添加过氧化氢RS计算算法

1 year ago · b503b52d02
5 changed files with 242 additions and 0 deletions
--- a/h2o2_computer/selftest.sh
+++ b/h2o2_computer/selftest.sh
@ -0,0 +1,5 @@
+#!/bin/bash
+set -e
+gcc zh2o2_computer.c zh2o2_computer_selftest.c -lm -o selftest.out
+./selftest.out
+rm selftest.out
--- a/h2o2_computer/zh2o2_computer.c
+++ b/h2o2_computer/zh2o2_computer.c
@ -0,0 +1,198 @@
+#include "zh2o2_computer.h"
+
+#include <stdio.h>
+/**
+ * @brief 水蒸气饱和压力(Pws)
+ *
+ * @param pressure
+ * @return float
+ *
+ * Demo:
+ *  waterVaporSaturationPressure(20+273.15) = 2339.25
+ *
+ */
+static double fn05_Pws(double Tk) {
+  // pws(Tk) = a0* exp(a1/Tk + a2 + a3*Tk + a4*Tk^2 + a5*ln(a6*Tk))
+
+  const double a0 = 1.000000;
+  const double a1 = -6.0969385 * 1000;
+  const double a2 = 2.12409642 * 10;
+  const double a3 = -2.711193 * 0.01;
+  const double a4 = 1.673952 * 0.00001;
+  const double a5 = 2.433502;
+  const double a6 = 1.000000;
+
+  double pws = a0 * exp(a1 / Tk + a2 + a3 * Tk + a4 * Tk * Tk + a5 * log(a6 * Tk));
+
+  return pws;
+}
+/**
+ * @brief
+ */
+
+/**
+ * @brief H2O2蒸气饱和压力Phs
+ *
+ * @param Tk
+ * @return double
+ *
+ * Formula:
+ *          Phs(Tk) = 10^(a0+a1/Tk+a2*log10(Tk)+a3*Tk) / (7.5006*0.001)
+ *
+ * Demo   :
+ *          h2o2VaporSaturationPressure(20+273.15) = 181.295
+ *
+ */
+static double fn06_Phs(double Tk) {
+  const double a0 = 44.576;
+  const double a1 = -4025;
+  const double a2 = -12.996;
+  const double a3 = 0.004605;
+
+  double phs = pow(10, a0 + a1 / Tk + a2 * log10(Tk) + a3 * Tk) / (7.5006 * 0.001);
+  return phs;
+}
+
+/**
+ * @brief fn01_Ph
+ *
+ * @param p
+ * @param pws
+ * @return double
+ *
+ */
+static double fn01_Ph(double ppm, double atmospheric_pressure) { return ppm * atmospheric_pressure / 1000000.0; }
+
+/**
+ * @brief fn10_Yw Eq10(Page7)
+ *
+ * @param Xw
+ * @return double
+ */
+static double fn10_Yw(double Xw, double Tk) {
+  // Yw = exp{ (1-Xw)^2 * (1/RTk) * [B0+B1(1-4Xw)+B2(1-2Xw)*(1-6Xw)] }
+  // R  = Universal gas constant 1.986 cal∙K-1
+
+  double       b0 = -1017 + 0.97 * Tk;
+  const double b1 = 85;
+  const double b2 = 13;
+  const double R  = 1.986;
+
+  double Yw = exp((1 - Xw) * (1 - Xw) * (1 / Tk) * (b0 + b1 * (1 - 4 * Xw) + b2 * (1 - 2 * Xw) * (1 - 6 * Xw)) / R);
+
+  return Yw;
+}
+
+/**
+ * @brief fn11_Yh Eq11(Page7)
+ *
+ * @param Xw
+ * @param Tk
+ * @return double
+ */
+static double fn11_Yh(double Xw, double Tk) {
+  // Yh = exp{ Xw^2/RTk * [B0+B1(3-4Xw)+B2(1-2Xw)*(5-6Xw)] }
+
+  double       b0 = -1017 + 0.97 * Tk;
+  const double b1 = 85;
+  const double b2 = 13;
+  const double R  = 1.986;
+
+  double Yh = exp((Xw * Xw / (Tk * R)) * (b0 + b1 * (3 - 4 * Xw) + b2 * (1 - 2 * Xw) * (5 - 6 * Xw)));
+  return Yh;
+}
+
+static double fn8_Pw(double Pws, double Phs, double Ph, double Yw, double Yh) { return Yw * Pws * (1 - Ph / (Yh * Phs)); }
+static double fn17_Pms(double Pws, double Xw, double Yw) { return Yw * Xw * Pws; }
+
+static void _iteration_pms_and_pw(double xh_avg, double h2o2ppm, double Tk, double AirPressure, double* Pms, double* Pw) {  //
+  double Xw = 1 - xh_avg;
+  double Yw = fn10_Yw(Xw, Tk);
+  double Yh = fn11_Yh(Xw, Tk);
+
+  double Pws = fn05_Pws(Tk);
+  double Phs = fn06_Phs(Tk);
+
+  double Ph = fn01_Ph(h2o2ppm, AirPressure);
+
+  *Pw  = fn8_Pw(Pws, Phs, Ph, Yw, Yh);
+  *Pms = fn17_Pms(Pws, Xw, Yw);
+}
+
+static double t2k(double t) { return t + 273.15; }
+
+/*******************************************************************************
+ *                                   Extern                                    *
+ *******************************************************************************/
+
+double zh2o2_compute_pms(double ppm, double Tk, double AirPressure) {
+  double Xmin = 0;
+  double Xmax = 1;
+  double Pms, Pw;
+
+  int iteration_times = 0;
+  while (true) {
+    iteration_times++;
+
+    if (iteration_times > 100) {
+      printf("%s:%d warninging......,iteration_times > 100\n", __FILE__, __LINE__);
+      break;
+    }
+
+    double xh_avg = (Xmin + Xmax) / 2;
+    if (fabs(Xmax - Xmin) < 0.00001) {
+      break;
+    }
+
+    _iteration_pms_and_pw(xh_avg, ppm, Tk, AirPressure, &Pms, &Pw);
+    // printf("%d xh_avg=%lf,Pms=%lf,Pw=%lf\n", iteration_times, xh_avg, Pms, Pw);
+
+    if (Pw > Pms) {
+      Xmax = xh_avg;
+    } else {
+      Xmin = xh_avg;
+    }
+  }
+
+  return Pms;
+}
+
+double zh2o2_compute_rs(double ppm, double Tk, double rh100, double AirPressure) {
+  double Pms   = zh2o2_compute_pms(ppm, Tk, AirPressure);
+  double Pw    = fn05_Pws(Tk) * rh100 / 100;
+  float  rs100 = Pw / Pms * 100.0;
+  if (rs100 > 100) rs100 = 100;
+  if (rs100 < 0) rs100 = 0;
+  return rs100;
+}
+
+double zh2o2_t2k(double temperatureC) { return temperatureC + 273.15; }
+
+#define ZCHECK_EQ(val0, val1, permision)                                                                                    \
+  {                                                                                                                         \
+    double val0_ = val0;                                                                                                    \
+    double val1_ = val1;                                                                                                    \
+    if (fabs(val0_ - val1_) > permision) {                                                                                  \
+      printf("%s:%d ZCHECK_EQ(%s,%s,%s) failed: %lf != %lf\n", __FILE__, __LINE__, #val0, #val1, #permision, val0_, val1_); \
+      return -1;                                                                                                            \
+    }                                                                                                                       \
+  }
+
+int zh2o2_computer_self_test() {
+  ZCHECK_EQ(fn05_Pws(t2k(20)), 2339.25, 0.01);
+  ZCHECK_EQ(fn06_Phs(t2k(20)), 181.295, 0.01);
+  ZCHECK_EQ(fn01_Ph(500, 101325), 50.6625, 0.01);
+  ZCHECK_EQ(fn10_Yw(0.5, t2k(20)), 0.703911, 0.01);
+  ZCHECK_EQ(fn11_Yh(0.5, t2k(20)), 0.757218, 0.01);
+
+  double Pms, Pw;
+  _iteration_pms_and_pw(0.5, 500, t2k(20), 101325, &Pms, &Pw);
+  ZCHECK_EQ(Pw, 1038.94585, 0.01);
+  ZCHECK_EQ(Pms, 823.311334, 0.01);
+
+  ZCHECK_EQ(zh2o2_compute_pms(500, t2k(20), 101325), 1063.920083, 0.01);
+  ZCHECK_EQ(zh2o2_compute_rs(500, t2k(20), 42, 101325), 92.345718, 0.01);
+
+  printf("Test passed\n");
+  return 0;
+}
--- a/h2o2_computer/zh2o2_computer.h
+++ b/h2o2_computer/zh2o2_computer.h
@ -0,0 +1,35 @@
+#pragma once
+#include <math.h>
+#include <stdbool.h>
+#include <stdint.h>
+
+
+/**
+ * @brief 计算 water vapor saturation pressure  in an H2O2 vapor environment
+ *
+ * @param ppm 过氧化氢浓度
+ * @param Tk  温度(K)
+ * @param AirPressure 大气压
+ * @return double
+ */
+double zh2o2_compute_pms(double ppm, double Tk, double AirPressure);
+/**
+ * @brief 计算 过氧气化氢相对饱和度(RS)
+ *
+ * @param ppm 过氧化氢浓度
+ * @param Tk  温度(K)
+ * @param rh100 相对湿度(RH%)
+ * @param AirPressure 大气压
+ * @return double
+ */
+double zh2o2_compute_rs(double ppm, double Tk, double rh100, double AirPressure);
+/**
+ * @brief 计算 摄氏度转开尔文
+ *
+ * @param temperatureC
+ * @return double
+ */
+double zh2o2_t2k(double temperatureC);
+
+
+int zh2o2_computer_self_test();
--- a/h2o2_computer/zh2o2_computer_selftest.c
+++ b/h2o2_computer/zh2o2_computer_selftest.c
@ -0,0 +1,3 @@
+#include "zh2o2_computer.h"
+
+int main(int argc, char* const argv[]) { zh2o2_computer_self_test(); }
--- a/src/api/cmds/device_state_cmd_impl.cpp
+++ b/src/api/cmds/device_state_cmd_impl.cpp
@ -68,6 +68,7 @@ json DeviceStateCmdImpl::createSensorDataJson() {
 #if (defined PROJECT_TYPE_DRAW_BAR_BOX)
  report["draw_bar_box_disinfection"]["pressure"][0]                 = m_deviceIoControlService->DBDB__readPressureSensor(1);        // 液位
  report["draw_bar_box_disinfection"]["pressure"][1]                 = m_deviceIoControlService->DBDB__readPressureSensor(2);        // 空压机压力
+  report["draw_bar_box_disinfection"]["pressure"][2]                 = m_deviceIoControlService->DBDB__readPressureSensor(3);        // 空压机压力
  report["draw_bar_box_disinfection"]["heatingStrip"]["current"]     = m_deviceIoControlService->DBDB__heaterReadElectricCurrent();  // 电流
  report["draw_bar_box_disinfection"]["heatingStrip"]["temperature"] = m_deviceIoControlService->DBDB__heaterReadTemperatureData();  // 温度
 #endif