当前位置：首页 > news >正文

手把手教你修改MFiX源代码：扩展Sutherland公式支持多种气体粘度计算

news 2026/6/14 8:43:22

手把手教你修改MFiX源代码：扩展Sutherland公式支持多种气体粘度计算

在计算流体动力学（CFD）模拟中，气体粘度的准确计算对结果可靠性至关重要。MFiX作为开源CFD软件，默认采用Sutherland公式计算空气粘度，但实际工程常涉及多种气体混合场景。本文将深入解析如何通过修改源代码实现多气体粘度计算，并提供完整的参数配置方法与验证方案。

1. Sutherland公式原理与MFiX实现分析

Sutherland公式是描述气体粘度随温度变化的经典半经验模型，其通用形式为：

μ = μ_ref * (T/T_ref)^(3/2) * (T_ref + S)/(T + S)

其中：

μ_ref为参考温度下的粘度
T_ref为参考温度（通常取273K）
S为Sutherland常数（气体特性参数）

MFiX 19.3.1版本中，相关实现位于calc_mu_g.f文件的136行：

MU_G(IJK) = to_SI*1.7D-4 * (T_G(IJK)/273.0D0)**1.5D0 * (383.D0/(T_G(IJK)+110.D0))

关键参数说明：

1.7D-4：空气在273K时的动力粘度（Pa·s）
383.D0：空气的Sutherland参考温度
110.D0：空气的Sutherland常数

注意：公式中的to_SI是单位转换系数，确保最终输出为国际单位制

2. 源代码改造：实现多气体支持

2.1 修改前的准备工作

备份原始文件：

cp mfix-19.3.1/model/calc_mu_g.f calc_mu_g.f.bak

创建气体参数数据库：建议在文件开头添加气体参数声明块：

! Gas property database MODULE GAS_PROPERTIES IMPLICIT NONE ! Format: gas_name, mu_ref, T_ref, S_const CHARACTER(LEN=20), PARAMETER :: GAS_LIST(3) = & (/'AIR ', 'CO2 ', 'H2O '/) DOUBLE PRECISION, PARAMETER :: MU_REF(3) = & (/1.7D-4, 1.37D-4, 1.12D-4/) DOUBLE PRECISION, PARAMETER :: T_REF(3) = & (/273.0D0, 273.0D0, 350.0D0/) DOUBLE PRECISION, PARAMETER :: S_CONST(3) = & (/110.0D0, 222.0D0, 1064.0D0/) END MODULE GAS_PROPERTIES

2.2 核心算法改造

替换原136行代码为以下多气体支持版本：

USE GAS_PROPERTIES !... INTEGER :: gas_type ! Get gas type from user input (1=AIR, 2=CO2, 3=H2O) gas_type = 1 ! Default to air MU_G(IJK) = to_SI * MU_REF(gas_type) * & (T_G(IJK)/T_REF(gas_type))**1.5D0 * & ((T_REF(gas_type)+S_CONST(gas_type))/(T_G(IJK)+S_CONST(gas_type)))

改进说明：

通过gas_type参数控制气体选择
所有气体参数集中管理，便于扩展
保持原计算结构，确保数值稳定性

3. 参数配置与扩展方法

3.1 常见气体Sutherland参数

气体	μ_ref (Pa·s)	T_ref (K)	S常数
空气	1.7×10⁻⁴	273	110
CO₂	1.37×10⁻⁴	273	222
H₂O	1.12×10⁻⁴	350	1064
N₂	1.66×10⁻⁴	273	107
O₂	1.92×10⁻⁴	273	132

提示：新增气体只需扩展GAS_PROPERTIES模块中的数组即可

3.2 动态气体选择实现

对于需要运行时切换气体的场景，建议添加输入参数：

! 在输入模块中添加 INTEGER :: GAS_SELECTOR NAMELIST /GAS_INPUT/ GAS_SELECTOR ! 修改计算逻辑 gas_type = MIN(MAX(GAS_SELECTOR, 1), SIZE(GAS_LIST))

对应mfix.dat输入文件需新增：

&GAS_INPUT GAS_SELECTOR = 2 ! 1=Air, 2=CO2, 3=H2O /

4. 验证与调试指南

4.1 单元测试方法

创建测试函数验证不同温度下的计算结果：

SUBROUTINE TEST_SUTHERLAND() USE GAS_PROPERTIES IMPLICIT NONE DOUBLE PRECISION :: temp, mu INTEGER :: i PRINT *, "T(K) AIR CO2 H2O" DO i = 1, 10 temp = 273.0 + i*100.0 mu = MU_REF(1)*(temp/T_REF(1))**1.5*(T_REF(1)+S_CONST(1))/(temp+S_CONST(1)) PRINT "(F6.1,3E12.4)", temp, mu, & MU_REF(2)*(temp/T_REF(2))**1.5*(T_REF(2)+S_CONST(2))/(temp+S_CONST(2)), & MU_REF(3)*(temp/T_REF(3))**1.5*(T_REF(3)+S_CONST(3))/(temp+S_CONST(3)) END DO END SUBROUTINE

典型输出验证：

温度(K)	空气粘度(Pa·s)	CO₂粘度(Pa·s)	H₂O粘度(Pa·s)
300	1.85×10⁻⁵	1.48×10⁻⁵	0.99×10⁻⁵
500	2.71×10⁻⁵	2.45×10⁻⁵	1.72×10⁻⁵

4.2 常见问题排查

编译错误：
- 确保USE GAS_PROPERTIES语句位置正确
- 检查数组越界问题（gas_type应在有效范围内）
计算结果异常：
- 验证温度单位是否为开尔文
- 检查to_SI系数是否应用正确
扩展新气体无效：
- 确认GAS_PROPERTIES模块已更新
- 检查气体索引是否正确传递

5. 高级应用：混合气体粘度计算

对于混合气体场景，可采用Wilke混合规则：

FUNCTION MU_MIX(T, X, N) RESULT(mu_mix) DOUBLE PRECISION, INTENT(IN) :: T, X(N) INTEGER, INTENT(IN) :: N DOUBLE PRECISION :: mu_mix, phi(N,N), mu_pure(N) INTEGER :: i, j ! 计算各组分纯粘度 DO i = 1,N mu_pure(i) = MU_REF(i)*(T/T_REF(i))**1.5*(T_REF(i)+S_CONST(i))/(T+S_CONST(i)) END DO ! 计算相互作用参数 DO i = 1,N DO j = 1,N phi(i,j) = (1.0 + SQRT(mu_pure(i)/mu_pure(j))*(MW(j)/MW(i))**0.25)**2 / & SQRT(8.0*(1.0 + MW(i)/MW(j))) END DO END DO ! 计算混合粘度 mu_mix = 0.0 DO i = 1,N mu_mix = mu_mix + X(i)*mu_pure(i)/SUM(X(1:N)*phi(i,1:N)) END DO END FUNCTION

注意：此实现需要额外定义分子量数组MW和组分摩尔分数X

6. 性能优化建议

查表法加速：

! 预计算温度-粘度关系表 DOUBLE PRECISION, ALLOCATABLE :: MU_TABLE(:,:) INTEGER, PARAMETER :: N_TEMP = 1000 DOUBLE PRECISION :: TEMP_MIN = 200.0, TEMP_MAX = 2000.0 ! 初始化查表 ALLOCATE(MU_TABLE(N_TEMP, SIZE(GAS_LIST))) DO i = 1, N_TEMP temp = TEMP_MIN + (i-1)*(TEMP_MAX-TEMP_MIN)/(N_TEMP-1) DO j = 1, SIZE(GAS_LIST) MU_TABLE(i,j) = MU_REF(j)*(temp/T_REF(j))**1.5*(T_REF(j)+S_CONST(j))/(temp+S_CONST(j)) END DO END DO ! 查表计算 idx = INT((T_G(IJK)-TEMP_MIN)/(TEMP_MAX-TEMP_MIN)*(N_TEMP-1)) + 1 idx = MIN(MAX(idx, 1), N_TEMP) MU_G(IJK) = to_SI * MU_TABLE(idx, gas_type)