module constants
! 定义物理常数和模拟参数
implicit none
! 数学常数
real, parameter :: PI = 3.1415926535 ! 圆周率π
! 物理常数
real, parameter :: R_univ = 8.314 ! 通用气体常数 [J/(mol·K)]
real, parameter :: kB = 1.380649e-23 ! 玻尔兹曼常数 [J/K]
real, parameter :: NA = 6.02214076e23 ! 阿伏伽德罗常数
! 水特性参数
real, parameter :: M_water = 0.018 ! 水的摩尔质量 [kg/mol]
real, parameter :: R_vapor = R_univ/M_water ! 水蒸气气体常数 [J/(kg·K)]
real, parameter :: rho_l = 1000.0 ! 液态水密度 [kg/m?]
real, parameter :: sigma = 0.075 ! 水的表面张力 [N/m]
real, parameter :: L_evap = 2.5e6 ! 水的蒸发潜热 [J/kg]
real, parameter :: Pr_vapor = 0.71 ! 水蒸气的普朗特数
! 发动机固定参数
real, parameter :: De = 0.2 ! 喷口直径 [m]
real, parameter :: T0 = 1000.0 ! 喷口温度 [K]
real, parameter :: P0 = 50000.0 ! 喷口压力 [Pa]
real, parameter :: Ue = 2000.0 ! 喷口速度 [m/s]
! 喷口组分摩尔分数
real, parameter :: Xco2 = 0.16 ! CO2摩尔分数
real, parameter :: Xco = 0.15 ! CO摩尔分数
real, parameter :: Xh2o = 0.22 ! H2O摩尔分数
real, parameter :: Xn2_o2 = 1.0 - Xco2 - Xco - Xh2o ! N2和O2的摩尔分数
! 环境参数
real, parameter :: P_inf = 26500 ! 环境压力 [Pa]
real, parameter :: T_inf = 223.15 ! 环境温度 [K]
! 计算参数
real, parameter :: dt = 1e-3 ! 增大时间步长 [s]
real, parameter :: t_total = 100.0 ! 减少总模拟时间 [s]
real, parameter :: gamma = 1.4 ! 比热比
real, parameter :: cp_mix = 1000.0 ! 混合气体定压比热 [J/(kg·K)]
real, parameter :: cooling_time = 1.0 ! 温度下降特征时间 [s]
real, parameter :: pressure_time = 2.0 ! 压力下降特征时间 [s]
integer, parameter :: nsteps = nint(t_total/dt) ! 时间步数
! 混合气体摩尔质量计算
real, parameter :: M_mix = Xco2*0.044 + Xco*0.028 + Xh2o*0.018 + &
Xn2_o2*0.78*0.028 + Xn2_o2*0.22*0.032
real, parameter :: R_mix = R_univ / M_mix ! 混合气体气体常数
end module constants
program contrail_simulation
use constants
implicit none
! 定义状态变量
real :: time, dist, T, P, Yv, Yl, N_drop, r_drop, rho
integer :: i, output_counter
real :: p_sat_val, S
logical :: condensation_occurred
! 打开输出文件
open(unit=10, file='contrail_results.txt', status='replace')
write(10, '(a)') 'Time(s), Distance(m), Temperature(K), Pressure(Pa), Yv, Yl, N_drop(m^{-3}), r_drop(m), Saturation'
! 设置初始条件 (喷口处)
P = P0
T = T0
rho = P / (R_mix * T)
Yv = (Xh2o * M_water) / M_mix
Yl = 0.0
N_drop = 0.0
r_drop = 1e-9 ! 小初始值避免除零
time = 0.0
dist = 0.0
output_counter = 0
! 计算初始饱和比
p_sat_val = p_sat(T)
S = (Yv * P * M_mix / M_water) / p_sat_val
! 输出初始状态
write(10, '(9g16.8)') time, dist, T, P, Yv, Yl, N_drop, r_drop, S
print *, '初始状态: T=', T, 'K, P=', P, 'Pa, Yv=', Yv, 'S=', S
! 主时间循环
do i = 1, nsteps
time = time + dt
dist = dist + Ue * dt
! 1. 压力随时间逐渐降低到环境压力
P = P_inf + (P0 - P_inf) * exp(-time/pressure_time)
! 2. 温度随时间逐渐降低到环境温度
T = T_inf + (T0 - T_inf) * exp(-time/cooling_time)
! 3. 更新密度
rho = P / (R_mix * T)
! 4. 计算当前饱和比
p_sat_val = p_sat(T)
S = (Yv * P * M_mix / M_water) / p_sat_val
! 5. 冷凝计算标志
condensation_occurred = .false.
! 6. 只有在过饱和(S>1)时才进行冷凝计算
if (S > 1.0) then
call calculate_condensation(T, P, rho, S, Yv, Yl, N_drop, r_drop, condensation_occurred)
end if
! 7. 每100步输出一次详细信息
if (mod(i, 100) == 0) then
print *, 'Step ', i, ': T=', T, 'K, P=', P, 'Pa, Yv=', Yv, 'Yl=', Yl
print *, ' N_drop=', N_drop, 'm^-3, r_drop=', r_drop*1e6, 'μm, S=', S
print *, ' 冷凝发生:', condensation_occurred
end if
! 8. 每10步输出一次结果到文件
output_counter = output_counter + 1
if (output_counter >= 10) then
write(10, '(9g16.8)') time, dist, T, P, Yv, Yl, N_drop, r_drop, S
output_counter = 0
end if
end do
close(10)
print *, '模拟完成! 结果已保存到 contrail_results.txt'
print *, '最终状态: T=', T, 'K, P=', P, 'Pa, Yl=', Yl, 'r_drop=', r_drop*1e6, 'μm'
contains
! 冷凝计算子程序
subroutine calculate_condensation(T, P, rho, S, Yv, Yl, N_drop, r_drop, occurred)
real, intent(in) :: T, P, rho, S
real, intent(inout) :: Yv, Yl, N_drop, r_drop
logical, intent(out) :: occurred
real :: J_nuc, r_c, drdt, m_v
real :: rho_g, term1, term2, lambda_v, Kn, denominator, numerator
real :: initial_r_drop, initial_N_drop
real :: m_single, delta_mu, G_c, q_c = 1.0 ! 蒸发系数取1
! 记录初始状态
initial_r_drop = r_drop
initial_N_drop = N_drop
occurred = .false.
! 检查过饱和程度是否足够
if (S < 1.0) then
return
end if
! 1. 计算临界半径 - 使用公式(2.8)
if (log(S) > 1e-10) then
r_c = 2.0 * sigma / (rho_l * R_vapor * T * log(S))
else
r_c = 1e-6 ! 设置一个合理的最小值
end if
! 2. 计算成核率 - 使用公式(2.9)
m_single = M_water / NA ! 单个水分子的质量
rho_g = Yv * rho ! 水蒸气密度
! 计算临界吉布斯自由能 G_c = 4πr_c?σ/3
G_c = 4.0 * PI * r_c**2 * sigma / 3.0
! 成核率计算公式
term1 = q_c * rho_g**2 / rho_l
term2 = sqrt(2.0 * sigma / (PI * m_single**3))
J_nuc = term1 * term2 * exp(-G_c / (kB * T))
! 限制成核率范围
J_nuc = min(max(J_nuc, 0.0), 1e20)
! 3. 计算液滴生长率 - 使用公式(2.13)
lambda_v = 0.025 ! 蒸汽导热系数 [W/(m·K)]
! 计算Knudsen数
Kn = lambda_v / (P * r_drop) * sqrt(PI * R_mix * T / 2.0)
! 计算分子和分母
numerator = lambda_v * (T - T_inf) * (1.0 - r_c / max(r_drop, 1e-10))
denominator = 1.0 + (2.0 / (1.5 * Pr_vapor)) * (gamma / (1.0 + gamma)) * Kn
! 计算生长率
drdt = numerator / (rho_l * L_evap * denominator)
! 限制生长率范围
drdt = min(max(drdt, -1e-4), 1e-4)
! 4. 计算相变质量
m_v = J_nuc * (4.0/3.0 * PI * r_c**3 * rho_l) + &
N_drop * rho_l * 4.0 * PI * r_drop**2 * drdt
! 限制相变质量不超过可用水蒸气的20%
m_v = min(m_v, Yv * rho * 0.2)
! 5. 更新状态变量
if (m_v > 0.0) then
N_drop = N_drop + J_nuc * dt
Yl = Yl + (m_v / rho_l) * dt
Yv = Yv - (m_v / rho) * dt
! 6. 更新液滴半径
if (N_drop > 0.0) then
r_drop = (3.0 * Yl / (4.0 * PI * rho_l * N_drop))**(1.0/3.0)
end if
occurred = .true.
end if
! 输出显著变化
if (abs(r_drop - initial_r_drop) > 1e-12 .or. &
abs(N_drop - initial_N_drop) > 1e-12) then
print *, '冷凝更新: r_drop=', r_drop*1e6, 'μm, N_drop=', N_drop, 'm^-3, m_v=', m_v
end if
end subroutine calculate_condensation
! 饱和蒸汽压计算函数
function p_sat(T) result(p)
real, intent(in) :: T
real :: p, t_c
t_c = T - 273.15
if (t_c >= 0) then
p = 611.2 * exp(17.62 * t_c / (243.12 + t_c))
else
p = 611.2 * exp(22.46 * t_c / (272.62 + t_c))
end if
end function p_sat
end program contrail_simulation为什么液滴半径在不断减小,一直减小
最新发布