基于RH输入自动匹配EAIM浓度与活度值的代码实现求助

架构老林

2026-8-23

需求说明

我是编程新手，以下是我的部分代码。需要实现：当设置特定RH值（范围0.2-0.9，步长0.1）时，自动从对应的EAIM矩阵中选取浓度和活度值。例如，当RH设为0.2时，需自动选取各矩阵第一行的数据。目前不知如何实现，恳请帮助！

原始代码

import math  # 原代码用到math.pi但未导入，此处补上

RH = 0.8 #input any value from 0.2-0.9 with 0.1 steps 
OrgMF = 0.1 #no units - range from 0.1 to 0.9 - organic to sulfate mass fraction
Rc = 200*10**-9 #m (nm) - inorganic particle radius
Vinorg = 4*math.pi*(Rc**3)/3 #m3 - volume inorganic 
Minorg = Vinorg*1.5 #mass inorganic
Morg = Minorg/(1/OrgMF-1) #mass organic 
Vorg = Morg/0.9 #m3 - volume organic 
Vtot = Vinorg + Vorg #m3 - volume total
Rtot = (3*Vtot/4/math.pi)**(1/3) #m - radius total
Lorg = Rtot - Rc #m - thickness of organic coating
Ma = 0.1191390 #kg/mol - mass of IEPOX
Sa = 4*3.1415*((Rc+Lorg)**2) #m2 - surface area of particles 
# 原代码中Ptot未定义，需补充实际值才能运行
Dg = 1e-4*0.1/Ptot #m2/s
a = 1 #no units - mass accomodation 
V = (4/3)*3.1415*((Rc+Lorg)**3) #m3 - particle vol 
T = 280.0 #K
R = 0.082057 #L*atm/k*mol 
R2 = 8.314 #j/k
Haq = 1*10**8 #M/atm - henrys law constant for IEPOX
Horg = 6*10**5 #M/atm - henry's law for org layer 
Dorg = 100*((1.38*10**-23)*T)/(6*3.1315*(5*10**5)*(315*10**-12)) #cm2/s - k*t/6*pi*viscosity*molecular radius - diffusion of IEPOX in org 
w = ((8*R2*T)/(3.1415*Ma))**0.5  #m/s - book equ for mean speed (281 in paper?)
kh = 3.6*(10**-2) #M-1s-1 - rate constant of H3O+ 
kso4 = 1*(10**-4) #M-1s-1 - ballpark - rate constant of nucleophile SO4-2 
khso4 = 7.3*(10**-4) #M-1s-1 rate consant of general acid HSO4-

#RH vs H concentration, activity
H_eaim = [[1.6*10**-5, 5.28*10**2], #RH 0.2 
          [1.67*10**-5, 1.6*10**2], #RH 0.3 
          [1.85*10**-5, 5.19*10**1], #RH 0.4
          [2.13*10**-5, 1.89*10**1], #RH 0.5 
          [2.53*10**-5, 7.41*10**0], #RH 0.6
          [3.07*10**-5, 3.1*10**0],  #RH 0.7
          [3.75*10**-5, 1.39*10**0], #RH 0.8
          [4.37*10**-5, 7.01*10**-1]] #RH 0.9
      
#RH vs HSO4 concentrations,  activity 
HSO4_eaim = [[8.4*10**-5, 1.97*10**0], 
             [8.33*10**-5, 1.89*10**0], 
             [8.15*10**-5, 1.81*10**0], 
             [7.87*10**-5, 1.7*10**0], 
             [7.47*10**-5, 1.55*10**0], 
             [6.93*10**-5, 1.37*10**0], 
             [6.25*10**-5, 1.14*10**0], 
             [5.63*10**-5, 8.49*10**-1]]

#RH vs SO4 concentrations, activity vs RH value 
SO4_eaim = [[1.6*10**-5, 6.09*10**-5], 
            [1.67*10**-5, 2.22*10**-4], 
            [1.85*10**-5, 6.25*10**-4], 
            [2.13*10**-5, 1.53*10**-3], 
            [2.53*10**-5, 3.02*10**-3], 
            [3.07*10**-5, 5.25*10**-3], 
            [3.75*10**-5, 8.3*10**-3], 
            [4.37*10**-5, 1.42*10**-2]]


H =  #M of H+ - EAIM CONCENTRATION - need to choose from the eaim values above for all M concentrations  
ah =  #EAIM ACTIVITY
         
SO4 = 3*10**-2 #M of nucleophile SO4-2 EAIM CONCENTRATION 
HSO4 = 1*10**0 #M of general acid HSO4- EAIM CONCENTRATION
         
Mt = 0.13615 #kg/mol mass of tetrol 
Ms = 0.216123 #kg/mol mass of sulfate product 
        
Htet = 1*10**8 #M/atm tetrol estimate  
Hs = 1*10**16 #M/atm sulfate product estimate
    
Kaq = (kh*H) + (kso4*SO4*ah) + (khso4*HSO4) #s-1 - reaction rate constant from gaston - using auto-selected eaim concentration/activity values to calculate
faqH = (kh*H)/Kaq
faqSO4 = ((kso4*SO4*ah) + (khso4*HSO4))/Kaq

实现方案

核心思路

RH值从0.2到0.9，步长0.1，正好对应矩阵的8个行索引（0到7）。通过计算RH对应的索引，直接提取对应行的浓度和活度值。为避免浮点数精度误差，先将RH转换为保留1位小数的精确值。

修改后的完整代码

import math

RH = 0.8 #input any value from 0.2-0.9 with 0.1 steps 
OrgMF = 0.1 #no units - range from 0.1 to 0.9 - organic to sulfate mass fraction
Rc = 200*10**-9 #m (nm) - inorganic particle radius
Vinorg = 4*math.pi*(Rc**3)/3 #m3 - volume inorganic 
Minorg = Vinorg*1.5 #mass inorganic
Morg = Minorg/(1/OrgMF-1) #mass organic 
Vorg = Morg/0.9 #m3 - volume organic 
Vtot = Vinorg + Vorg #m3 - volume total
Rtot = (3*Vtot/4/math.pi)**(1/3) #m - radius total
Lorg = Rtot - Rc #m - thickness of organic coating
Ma = 0.1191390 #kg/mol - mass of IEPOX
Sa = 4*3.1415*((Rc+Lorg)**2) #m2 - surface area of particles 
# 补充Ptot示例值，需根据实际情况修改
Ptot = 1.0
Dg = 1e-4*0.1/Ptot #m2/s
a = 1 #no units - mass accomodation 
V = (4/3)*3.1415*((Rc+Lorg)**3) #m3 - particle vol 
T = 280.0 #K
R = 0.082057 #L*atm/k*mol 
R2 = 8.314 #j/k
Haq = 1*10**8 #M/atm - henrys law constant for IEPOX
Horg = 6*10**5 #M/atm - henry's law for org layer 
Dorg = 100*((1.38*10**-23)*T)/(6*3.1315*(5*10**5)*(315*10**-12)) #cm2/s - k*t/6*pi*viscosity*molecular radius - diffusion of IEPOX in org 
w = ((8*R2*T)/(3.1415*Ma))**0.5  #m/s - book equ for mean speed (281 in paper?)
kh = 3.6*(10**-2) #M-1s-1 - rate constant of H3O+ 
kso4 = 1*(10**-4) #M-1s-1 - ballpark - rate constant of nucleophile SO4-2 
khso4 = 7.3*(10**-4) #M-1s-1 rate consant of general acid HSO4-

#RH vs H concentration, activity
H_eaim = [[1.6*10**-5, 5.28*10**2], #RH 0.2 
          [1.67*10**-5, 1.6*10**2], #RH 0.3 
          [1.85*10**-5, 5.19*10**1], #RH 0.4
          [2.13*10**-5, 1.89*10**1], #RH 0.5 
          [2.53*10**-5, 7.41*10**0], #RH 0.6
          [3.07*10**-5, 3.1*10**0],  #RH 0.7
          [3.75*10**-5, 1.39*10**0], #RH 0.8
          [4.37*10**-5, 7.01*10**-1]] #RH 0.9
      
#RH vs HSO4 concentrations,  activity 
HSO4_eaim = [[8.4*10**-5, 1.97*10**0], 
             [8.33*10**-5, 1.89*10**0], 
             [8.15*10**-5, 1.81*10**0], 
             [7.87*10**-5, 1.7*10**0], 
             [7.47*10**-5, 1.55*10**0], 
             [6.93*10**-5, 1.37*10**0], 
             [6.25*10**-5, 1.14*10**0], 
             [5.63*10**-5, 8.49*10**-1]]

#RH vs SO4 concentrations, activity vs RH value 
SO4_eaim = [[1.6*10**-5, 6.09*10**-5], 
            [1.67*10**-5, 2.22*10**-4], 
            [1.85*10**-5, 6.25*10**-4], 
            [2.13*10**-5, 1.53*10**-3], 
            [2.53*10**-5, 3.02*10**-3], 
            [3.07*10**-5, 5.25*10**-3], 
            [3.75*10**-5, 8.3*10**-3], 
            [4.37*10**-5, 1.42*10**-2]]

# 计算RH对应的矩阵索引，处理浮点数精度
rh_rounded = round(RH, 1)
index = int((rh_rounded - 0.2) / 0.1)

# 从EAIM矩阵提取对应数据
H, ah = H_eaim[index]
SO4, _ = SO4_eaim[index]  # 提取浓度，若需活度可替换为变量名
HSO4, _ = HSO4_eaim[index]  # 提取浓度，若需活度可替换为变量名
         
Mt = 0.13615 #kg/mol mass of tetrol 
Ms = 0.216123 #kg/mol mass of sulfate product 
        
Htet = 1*10**8 #M/atm tetrol estimate  
Hs = 1*10**16 #M/atm sulfate product estimate
    
Kaq = (kh*H) + (kso4*SO4*ah) + (khso4*HSO4) #s-1 - reaction rate constant from gaston - using auto-selected eaim concentration/activity values to calculate
faqH = (kh*H)/Kaq
faqSO4 = ((kso4*SO4*ah) + (khso4*HSO4))/Kaq

# 测试输出，验证提取结果
print(f"当前RH={RH}，提取的H浓度：{H}，活度：{ah}")
print(f"SO4浓度：{SO4}，HSO4浓度：{HSO4}")
print(f"计算得到的Kaq：{Kaq}")

注意事项

原代码中Ptot未定义，需补充实际业务值才能正常运行，示例中暂设为1.0。
浮点数处理：使用round(RH,1)确保RH值精确到一位小数，避免因浮点数存储误差导致索引计算错误。
若需要使用HSO4或SO4的活度，只需将代码中的_替换为对应的变量名即可，例如SO4_conc, SO4_act = SO4_eaim[index]，后续按需调用SO4_act。

内容的提问来源于stack exchange，提问作者livlh

基于RH输入自动匹配EAIM浓度与活度值的代码实现求助

核心思路 #

修改后的完整代码 #

注意事项 #

核心思路

修改后的完整代码

注意事项