企业官网建设流程全解析

• 科学计算

本文围绕 CoolProp 的一项核心数值问题展开：如何在delta=0（即零密度）处安全地计算 Helmholtz 残差自由能的偏导数，并从中提取精确的维里系数 B(T)、C(T)。文档docs/superpowers/derivations/virial-axy.md给出了完整的符号推导（由 sympy 完成、机器校验），配套脚本为dev/derivations/virial_axy_derivations.py（输出见dev/derivations/virial_axy_output.txt）。读完本文，你将掌握：Axy无量纲导数群的定义与"按delta^y缩放形式累加、避免运行时计算delta*(something/delta)"的编码原则、理想气体主导项/广义指数项（GenExp）/非解析项（IAPWS-95 型）/Gao-B/SAFT 缔合/广义立方型六类 Helmholtz 项在delta=0处的有限性证明与闭式表达式、以及从c1/c2泰勒系数重建 B(T)、C(T)、dB/dT、dC/dT 的完整公式链。

1. 问题背景：delta=0处的0*Inf = NaN

CoolProp 采用无量纲化 Helmholtz 自由能模型：a = a0(τ,δ) + ar(τ,δ)，其中τ = Tc/T为对比温度的倒数、δ = ρ/ρ_r为对比密度。在主流实现中，残差部分的d(ar)/d(delta)是通过对数导数形式再乘1/delta编码的——参见 src/Helmholtz.cpp 中 GenExp 累加段（:124-133）与末尾的one_over_delta乘法（:274-289），以及理想气体主导项IdealHelmholtzLead::all（:1280-1290，其中dalphar_ddelta += 1.0/delta）。这类1/delta因子在delta=0处发散，而压力公式p = rho*R*T*(1 + delta*dalphar_dDelta)又在delta=0处乘上一个显式的delta，于是产生0*Inf = NaN。

文档中已给出可复现的失败现象：对任意 HEOS 流体执行update(DmolarT, 0, T)都会抛出"p is not a valid number"——该异常信息正是 src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp 中post_update()后置校验所抛出的ValueError（见:1730）。压力计算的底层逻辑在calc_pressure_nocache（:2723起），其核心行p_last = rhomolar * R_u * T * (1 + delta * d.dalphar_ddelta)（:2765）正是 0*Inf 的来源。

解决办法：定义无量纲导数群

Axy = tau^x * delta^y * d^{x+y}(a) / (dtau^x ddelta^y)

其中delta^y缩放因子在解析层面与1/delta相消，使Axy在delta=0处有限。文档明确了两条 API 路径：get_Ar(x,y)返回残差部分的Axy，get_Aig(x,y)返回理想气体部分的Axy（teqp 风格命名）；而维里系数走另一条"零密度专用路径"，仅从delta低幂次项提取c1/c2。

2. 核心原则与符号约定

• 编码原则：每一个项都直接以delta^y缩放后的形式累加Axy，使源码中每个1/delta^k都能与显式的delta^k相消——绝不在运行时计算delta * (something/delta)。
• 维里系数身份：ar在delta=0处可展开为ar = c1*delta + c2*delta^2 + ...，因此c1、c2分别是ar对delta的一阶、二阶导数在 0 处的取值（带系数）：

c1 = d(ar)/d(delta)|_0 B(T) = c1 / rho_r c2 = (1/2) d^2(ar)/d(delta^2)|_0 C(T) = 2*c2 / rho_r^2 dB/dT = (dc1/dtau) * dtau/dT / rho_r dtau/dT = -T_r / T^2 dC/dT = 2*(dc2/dtau) * dtau/dT / rho_r^2

• 关键观察：c1只接收delta幂次d=1项的贡献，c2只接收d∈{1,2}项的贡献。这正是 GH #2991 所要求的"零密度处维里系数的独立路径"与通用（非维里）导数路径get_Ar/get_Aig的分工。
• sympy 陷阱：对符号正指数d求limit(delta**d * ..., delta, 0)会错误地返回oo；因此有限性校验必须代入具体d值（d=1是关键最低幂次），而乘积形式本身对符号d是正确的。

3. 理想气体主导项：a0_lead = log(delta) + a1 + a2*tau

理想气体的典型零密度情形是d(a0)/d(delta) = 1/delta → Inf，但Aig(0,1)=1有限。sympy 推导出如下闭式（源码见 src/Helmholtz.cpp:1284-1289，当前实现正是log(delta) + a1 + a2*tau及1/delta、-1/delta/delta系列导数）：

其余理想气体项只含tau（如log(tau)类），对delta的导数恒为 0，仅贡献到A[itau][0] = tau^itau * d^itau(a0)/dtau^itau（对既有tau导数做机械缩放，本来就有限）。A00含有log(delta)（在 0 处非有限），但压力与维里系数从不读取A00。

4. 广义指数项（GenExp）：大多数流体的主力项

单个 GenExp 元素的形式为：

ar_elem = n*tau^t*delta^d*exp(u), u = -c*delta^l - omega*tau^m - eta1*(delta-eps1) - eta2*(delta-eps2)^2 - beta1*(tau-gam1) - beta2*(tau-gam2)^2

定义两个"不含1/delta"的括号：

Bdelta = d - c*l*delta^l - delta*(eta1 + 2*eta2*(delta-eps2)) (= delta * dln(ar_elem)/ddelta) Btau = t - omega*m*tau^m - tau*(beta1 + 2*beta2*(tau-gam2))

sympy 生成的Axy乘积形式（对d=1,2均验证在delta=0处有限），可直接作为 C 代码：

// A01 = delta*d(ar)/ddelta A01 = pow(delta,d)*n*pow(tau,t)*( -c*pow(delta,l)*l + d - delta*(eta1 + 2*eta2*(delta-epsilon2)) )*exp(u); // A02 = delta^2 * d2(ar)/ddelta2 A02 = -pow(delta,d)*n*pow(tau,t)*( c*pow(delta,l)*pow(l,2) - c*pow(delta,l)*l - d*(d-1) + 2*d*(c*pow(delta,l)*l + delta*(eta1+2*eta2*(delta-epsilon2))) + 2*pow(delta,2)*eta2 - pow(c*pow(delta,l)*l + delta*(eta1+2*eta2*(delta-epsilon2)), 2) )*exp(u); // A10 = tau*d(ar)/dtau A10 = pow(delta,d)*n*pow(tau,t)*( -m*omega*pow(tau,m) + t - tau*(beta1 - 2*beta2*(gamma2-tau)) )*exp(u);

A11、A12形式类似（完整 ccode 见脚本输出），同样在delta=0处有限。

维里系数。记u0 := u(tau, delta=0) = beta1*gam1 - beta1*tau - beta2*gam2^2 + 2*beta2*gam2*tau - beta2*tau^2 + eps1*eta1 - eps2^2*eta2 - omega*tau^m，则：

值得注意的细节：-c交叉项仅在l=1时进入c2（此时delta^{d+l}=delta^2）；l>=2时它不出现。dc1/dtau、dc2/dtau的闭式在脚本输出中给出，用于计算 dB/dT、dC/dT。

XiangDeiters项形如phi0 + acentric*phi1 + theta*phi2，其中每个phi_k都是 GenExp 子项——按上式逐子项应用再线性组合即可，无需新的推导。其源码结构见 src/Helmholtz.cppResidualHelmholtzXiangDeiters::all（:959-973），确实就是derivs = derivs + derivs0 + derivs1 * acentric + derivs2 * theta的线性叠加。

5. 非解析项（IAPWS-95 型；Water、CO2 等）

ar_elem = delta*n*DELTA^b*PSI theta = (1-tau) + A*((delta-1)^2)^(1/(2*beta)); PSI = exp(-C*(delta-1)^2 - D*(tau-1)^2) DELTA = theta^2 + B*((delta-1)^2)^a

该形式源码中没有1/delta（参见 src/Helmholtz.cpp:739-752的累加代码，alphar += delta * ni * DELTA_bi * PSI），因此原始导数在delta=0处本就有限；实现只需将既有贡献按delta^y * tau^x缩放。脚本以具体参数（n=0.7, a=3.5, b=0.85, beta=0.3, A=0.32, B=0.2, C=28, D=700, tau=1.1）做了数值验证：A01,A02,A10,A11,A12在delta→0时全部趋于 0 且有限。

维里系数：显式的delta因子给出

c1 = (n*DELTA^b*PSI)|_{delta=0} = n * ((1-tau+A)^2 + B)^b * exp(-C - D*(tau-1)^2) c2 = (1/2) d^2/ddelta^2|_0 (复用 delta=0 处的既有 d2alphar_ddelta2 贡献)

对测试参数c1 ≈ 1.35e-16——因为PSI(0)=exp(-C)，数值上指数级微小。这符合物理预期：非解析项局域在临界点附近，对第二维里系数 B 的贡献可忽略。有限且正确。

6. Gao-B 项

ar_elem = n*tau^t*exp(1/(b+beta*(tau-gam)^2)) * delta^d*exp(eta*(delta-eps)^2)

文档指出源码中/delta只是编码选择（旧版 src/Helmholtz.cpp 中该处的写法），功能上ar_elem在delta=0处是解析的——当前源码ResidualHelmholtzGaoB::all（:826起）正是按Ftau = pow(tau,t)*exp(1/(b+beta*pow(-gamma+tau,2)))与Fdelta = pow(delta,d)*exp(eta*pow(delta-epsilon,2))分离乘积实现（:839-840），不含1/delta。Axy 乘积形式（d=1,2有限）：

// A01 A01 = n*pow(tau,t)*( d*pow(delta,d) + 2*pow(delta,d+1)*eta*(delta-epsilon) ) * exp( (eta*(b+beta*pow(gamma-tau,2))*pow(delta-epsilon,2) + 1)/(b+beta*pow(gamma-tau,2)) ); // A02 A02 = pow(delta,d)*n*pow(tau,t)*( 4*d*delta*eta*(delta-epsilon) + d*(d-1) + 2*pow(delta,2)*eta*(2*eta*pow(delta-epsilon,2)+1) ) * exp( eta*pow(delta-epsilon,2) + 1.0/(b+beta*pow(gamma-tau,2)) );

（A10,A11,A12在脚本输出中；全部在 0 处有限。）

维里系数，记w := eta*eps^2 + 1/(b+beta*(tau-gam)^2)：

dc1/dtau、dc2/dtau同样在脚本输出中。

7. SAFT 缔合项

ar = m*a*(log(X) - X/2 + 1/2), X = 2/(sqrt(1+4*Deltabar*delta)+1) Deltabar = g(vbarn*delta)*(exp(epsbar*tau)-1)*kappabar, g(eta)=0.5*(2-eta)/(1-eta)^3

关键事实：X(delta=0) = 1（脚本中已验证），故ar(delta=0)=0。所有分母都是(2*Deltabar*X*delta+1)的幂，在delta=0处趋于 1，因此源码中既有的X导数闭式（src/Helmholtz.cpp:861-1110区间的Xdd/dXdd_d...系列，以及ResidualHelmholtzSAFTAssociating::all（:1228））在delta=0处全部有限；实现只需按delta^y*tau^x缩放。

脚本以具体参数（m=a=1, vbarn=0.45, epsbar=5, kappabar=0.002, tau=1.2）数值验证：A01,A02,A10,A11,A12全部有限且随delta→0趋于 0；c1=-0.402、c2=-0.129（有限非零——SAFT 确实对 B、C 有贡献）。这与非解析项形成鲜明对比：缔合作用在低密度极限并非可忽略。

8. 广义立方型项

立方型状态方程（PengRobinson 等）委托给m_abstractcubic->alphar(tau,delta,z,itau,idelta)返回原始混合偏导。立方 EOS 在rho=0处解析（alphar是delta线性式的对数形式），因此Axy = tau^x*delta^y * alphar(...,x,y)在delta=0处有限——没有1/delta需要相消，只需施加delta^y*tau^x缩放。立方型的调用路径与委托结构可结合 src/Helmholtz.cppResidualHelmholtzGeneralizedCubic::all（:776）理解。

9. 数值交叉验证：维里路径的基石

文档给出了关键的数值对拍：B*rho_r = c1 = lim_{delta→0} d(ar)/d(delta)，与 CoolProp 自带的dalphar_dDelta（内部存在1e-12hack）独立比较，T=300 K：

约 6 位有效数字的收敛确认了维里路径所依赖的泰勒系数身份。（Water 在δ=1e-3处有偏差——那是非解析项——但到δ=1e-6即收敛。）

10. 复现与验证方法

推导脚本可直接复现：python3 -u dev/derivations/virial_axy_derivations.py。它依赖 sympy，按项输出Axy闭式、ccode 及有限性检查，完整结果已固化在dev/derivations/virial_axy_output.txt。脚本内部逻辑要点：

• AXY = [(0,1),(0,2),(1,0),(1,1),(1,2)]枚举全部所需导数阶；
• GenExp 与 GaoB 走axy_symbolic（符号 + ccode）与virials_symbolic（对d∈{1,2,3}、l∈{1,2}的组合求c1/c2及dc1/dtau、dc2/dtau）；
• 非解析与 SAFT 因分数幂/根号会令符号化简爆炸，走numeric_check：在delta=1e-4/1e-7/1e-10三个尺度检查Axy有限性并数值求c1/c2；
• 关键注释明确提醒 sympy 对符号d求极限的错误行为。

文档还给出了实现阶段的三重验证门（oracle）：手推Axy与delta^y*tau^x*(原始导数)在远离 0 处吻合至 ~1e-9（快速路径 oracle），以及one_mcx多重复数自动微分（绝对 oracle）——后者与源码中#if ENABLE_CATCH保护的one_mcx方法（如 GaoB:917、XiangDeiters:975）对应。

11. 实现映射：从推导到源码改写

文档为实施分支（CoolProp-izv）给出了逐项的改写清单（下表行号对应目标实现修订版，本文已验证的当前源码位置在括号内标注）：

配套数据结构与 API 面（推导文档中的设计目标）：HelmholtzDerivatives存储A[itau][idelta]（itau∈{0,1}, idelta∈{0,1,2}）并暴露getA(itau,idelta)；后端提供calc_Ar/calc_Aig，公开get_Ar/get_Aig；压力改为calc_pressure = rho*R*T*(1+Ar(0,1))并对rho==0短路；维里系数直接从c1/c2重建。请注意：calc_Ar/calc_Aig/get_Ar/get_Aig属于推导文档描述的目标实现 API，在当前仓库主分支源码中尚未出现，属于实施分支（CoolProp-izv / GH #2991）的映射清单。

12. 总结

Axy方案的实质是用"代数先消去、再数值求值"替代"数值上 0*Inf"的病态组合：所有闭式均由 sympy 符号生成并被机器校验，有限性/维里主张全部可复现（dev/derivations/virial_axy_derivations.py+ 输出文件）；核心身份c1 = B*rho_r已与 CoolProp 自身导数引擎交叉验证。对想要在 CoolProp 中修复update(DmolarT,0,T)崩溃、或为自定义 Helmholtz 项推导零密度安全导数与维里系数的开发者而言，本文给出了从符号推导、闭式到源码改写的完整闭环。

• 科学计算