error #6649: This symbol has multiple INTENT statement/attribute declarations which is not allowed. [Y]
Compilation Aborted (code 1)
Compilation Aborted (code 1)
error #6406: Conflicting attributes or multiple declaration of name. [Y]
error #6406: Conflicting attributes or multiple declaration of name. [Y]
error #6649: This symbol has multiple INTENT statement/attribute declarations which is not allowed. [Y]修改以下程序! 凝结尾迹形成模拟程序
! 包含Schmidt-Appleman判据计算和Gyarmathy液滴生长模型
! 添加了完整的控制方程组求解
!======================================================================
module SchmidtAppleman
implicit none
!!定义常量
real, parameter :: a_water = 7.63 !!a,b为经验常数,水面取7.63和241.9 冰面取9.5和265.5
real, parameter :: b_water = 241.9
real, parameter :: Q = 35e3 !!燃料的燃烧热
real, parameter :: b = 265.5
real, parameter :: EIH2O = 1.25 !!水蒸气排放指数(假设为1.25)
real, parameter :: Cp = 1.004 !!空气的定压比热容
real, parameter :: epsilon = 0.6222 !!水与干空气的分子量之比
real, parameter :: a = 9.5
real, parameter :: E0 = 6.11 !!温度为0时,E0的数值,hPa(马格努斯经验公式)
real, parameter :: E0_water = 6.112
real, parameter :: EIT = 0.6 !!喷气发动机的平均推进效率
contains
subroutine check_contrail_formation()
real :: RHc !!温度T下的凝结尾形成的阈值湿度,计算可得
real :: RHi !!冰面饱和水汽,计算可得
real :: RHw !!相对于水面的饱和水汽压,kPa
real :: G !!计算可得
real :: Tc !!计算可得
real :: RH
real :: RH_0 = 80.0 !! 添加初始值80%相对湿度
real :: T=-50 !!大气环境温度 ℃
real :: P=26500 !!环境的大气压强
real :: e_T !!T温度下的饱和水汽压
real :: eTc !!温度T下的饱和水汽压,由马格努斯公式计算可得
integer :: have
RH = RH_0 / 100.0
RHw = E0_water*10**((a_water*T)/(b_water + T))
RHi = RHw*(6.061*exp(18.102*T/(249.52+T))/6.1162/exp(22.577*T/(237.78+T)))
G = (EIH2O * Cp * P) / (epsilon * Q * (1 - EIT))
Tc = -46.46 + 9.43*log(G-0.053) + 0.72*log(G-0.053)*log(G-0.053)
eTc = E0*10**(a*T/(b+T))
RHc = (G*(T-Tc) + eTc)/eTc
print *, "G, Tc, RHw =", G, Tc, RHw
If (RHc <= RH) then
If (RHc < 1) then
If (RHi < 1) then
have = 1
write (*,*) "emergence,能产生凝结尾迹"
else
have = 0
write (*,*) "will not emergence,不能产生凝结尾迹"
end if
end if
end if
end subroutine check_contrail_formation
end module SchmidtAppleman
!======================================================================
! 凝结尾迹综合模拟模块
! 包含气相方程(有限差分法)和液相方程(迎风格式)
!======================================================================
module contrail_simulation
use, intrinsic :: iso_fortran_env, only: dp => real64
implicit none
!------------------ 物理常数定义 ------------------!
real(dp), parameter :: pi = 3.1415926535897932_dp ! 圆周率
real(dp), parameter :: R_gas = 287.0_dp ! 空气气体常数 [J/(kg·K)]
real(dp), parameter :: R_v = 461.5_dp ! 水蒸气气体常数 [J/(kg·K)]
real(dp), parameter :: k_B = 1.380649e-23_dp ! 玻尔兹曼常数 [J/K]
real(dp), parameter :: NA = 6.02214076e23_dp ! 阿伏伽德罗常数 [mol⁻¹]
real(dp), parameter :: Cp_air = 1004.0_dp ! 空气定压比热容 [J/(kg·K)]
real(dp), parameter :: lambda_v = 0.026_dp ! 蒸汽导热系数 [W/(m·K)]
real(dp), parameter :: h_tg = 2.5e6_dp ! 汽化潜热 [J/kg]
real(dp), parameter :: Pr = 0.71_dp ! 普朗特数(空气)
real(dp), parameter :: lambda = 1.0_dp ! 分子自由程比
real(dp), parameter :: sigma = 0.072_dp ! 表面张力 [N/m]
real(dp), parameter :: rho_l = 1000.0_dp ! 液态水密度 [kg/m³]
real(dp), parameter :: mol_diameter = 3.0e-10_dp ! 水分子直径 [m]
real(dp), parameter :: mu = 1.8e-5_dp ! 空气动力粘度 [Pa·s]
real(dp), parameter :: k_thermal = 0.024_dp ! 空气热导率 [W/(m·K)]
!------------------ 环境参数 ------------------!
real(dp) :: Te = 223.15_dp ! 环境温度 [K]
real(dp) :: P_env = 26500.0_dp ! 环境总压 [Pa]
real(dp) :: Ue = 2000.0_dp ! 对流速度 [m/s] (修正为更合理的值)
!------------------ 射流参数 ------------------!
real(dp) :: T_jet = 1000.0_dp ! 射流温度 [K]
real(dp) :: U_jet = 2000.0_dp ! 射流速度 [m/s]
real(dp) :: P_jet = 50000.0_dp ! 射流压力 [Pa]
!------------------ 计算网格参数 ------------------!
integer, parameter :: nx = 100 ! 轴向网格数
integer, parameter :: nr = 50 ! 径向网格数
real(dp), parameter :: L = 50.0_dp ! 轴向计算域长度 [m]
real(dp), parameter :: R_max = 10.0_dp ! 最大径向距离 [m]
real(dp), parameter :: dx = L / nx ! 轴向步长 [m]
real(dp), parameter :: dr = R_max / nr ! 径向步长 [m]
real(dp), parameter :: dt = 1.0e-7_dp ! 时间步长 [s] (减小以适应高速射流)
!------------------ 场变量定义 ------------------!
! 气相场变量
real(dp), dimension(nx, nr) :: rho_g = 0.3_dp ! 气相密度 [kg/m³]
real(dp), dimension(nx, nr) :: u = 2000.0_dp ! 轴向速度 [m/s] (初始化为射流速度)
real(dp), dimension(nx, nr) :: v = 0.0_dp ! 径向速度 [m/s] (轴对称假设)
real(dp), dimension(nx, nr) :: p = 50000.0_dp ! 压力 [Pa] (初始化为射流压力)
real(dp), dimension(nx, nr) :: T_g = 1000.0_dp ! 温度 [K] (初始化为射流温度)
real(dp), dimension(nx, nr) :: h_g = 0.0_dp ! 焓 [J/kg]
! 液相场变量
real(dp), dimension(nx, nr) :: N_d = 1.0e-20_dp ! 液滴数目密度 [1/m³]
real(dp), dimension(nx, nr) :: Y_d = 0.0_dp ! 液滴质量分数 [kg/kg]
real(dp), dimension(nx, nr) :: r_d = 1.0e-9_dp ! 液滴半径 [m]
! 源项
real(dp), dimension(nx, nr) :: S_m = 0.0_dp ! 质量源项
real(dp), dimension(nx, nr) :: S_u = 0.0_dp ! 动量源项
real(dp), dimension(nx, nr) :: S_v = 0.0_dp ! 径向动量源项
real(dp), dimension(nx, nr) :: S_h = 0.0_dp ! 能量源项
real(dp) :: time = 0.0_dp ! 模拟时间 [s]
integer :: step_count = 0 ! 时间步计数
contains
!==================== 主计算程序 ====================!
subroutine solve()
real(dp) :: drdt, J ! 半径增长率 [m/s], 成核率
real(dp) :: mv_nucleation, mv_growth ! 成核与生长质量源项
real(dp) :: vapor_consumption ! 蒸气消耗量 [kg/m³]
real(dp) :: T_C, S, Kn, R_c, Delta_T ! 中间计算变量
real(dp) :: mean_free_path ! 分子平均自由程 [m]
real(dp) :: x_pos, r_pos ! 位置坐标
real(dp) :: conv_term, diff_term ! 对流项和扩散项
real(dp) :: rho_air ! 空气密度 [kg/m³]
real(dp) :: dudx, dudy, dvdx, dvdy ! 速度梯度
real(dp) :: dTdx, dTdy, dpdx, dpdy ! 梯度
real(dp) :: div_u, laplace_u ! 散度和拉普拉斯算子
integer :: i, j, step
real(dp) :: area_factor ! 轴对称面积因子
!=============== 初始化阶段 ================!
rho_air = P_env / (R_gas * Te) ! 计算空气密度
h_g = Cp_air * T_g ! 初始化焓值
! 设置射流初始条件 (整个射流区域)
do j = 1, nr
do i = 1, nx
rho_g(i,j) = P_jet / (R_gas * T_jet) ! 使用射流条件计算密度
u(i,j) = U_jet
v(i,j) = 0.0_dp
p(i,j) = P_jet
T_g(i,j) = T_jet
! 仅在中心区域初始化液滴
if (i <= nx/5 .and. j <= nr/2) then
N_d(i,j) = 1.0e15_dp
Y_d(i,j) = 0.01_dp
r_d(i,j) = 1.0e-6_dp
else
N_d(i,j) = 1.0e-20_dp
Y_d(i,j) = 0.0_dp
r_d(i,j) = 1.0e-9_dp
end if
end do
end do
! 打开输出文件
open(unit=10, file="contrail_evolution.csv", status="replace")
write(10, "(A)") "Step,Time[s],x[m],r[m],N[1/m3],Y[kg/kg],R[m],T[K],P[Pa],S"
! 打开云图输出文件
open(unit=20, file="temperature_field.csv", status="replace")
open(unit=30, file="pressure_field.csv", status="replace")
write(20, "(A)") "x[m],r[m],T[K]"
write(30, "(A)") "x[m],r[m],P[Pa]"
!=============== 时间推进循环 ================!
do step = 1, 200000 ! 增加时间步数以适应更小的时间步长
step_count = step_count + 1
! 0. 重置源项
S_m = 0.0_dp
S_u = 0.0_dp
S_v = 0.0_dp
S_h = 0.0_dp
! 1. 液相求解 (迎风格式)
do j = 2, nr-1
do i = 2, nx-1
! 计算位置坐标
x_pos = (i - 0.5_dp) * dx
r_pos = (j - 0.5_dp) * dr
area_factor = 1.0_dp / r_pos ! 轴对称因子
! 液滴数密度方程 (7)
conv_term = 0.0_dp
! 轴向对流 (迎风)
if (u(i,j) > 0.0_dp) then
conv_term = conv_term + u(i,j) * (N_d(i,j) - N_d(i-1,j)) / dx
else
conv_term = conv_term + u(i,j) * (N_d(i+1,j) - N_d(i,j)) / dx
end if
! 径向对流 (迎风)
if (v(i,j) > 0.0_dp) then
conv_term = conv_term + v(i,j) * (N_d(i,j) - N_d(i,j-1)) / dr
else
conv_term = conv_term + v(i,j) * (N_d(i,j+1) - N_d(i,j)) / dr
end if
! 添加轴对称项
conv_term = conv_term + v(i,j) * N_d(i,j) * area_factor
! 更新液滴数密度
N_d(i,j) = N_d(i,j) - dt * (conv_term - J)
! 液滴质量分数方程 (8)
conv_term = 0.0_dp
! 轴向对流 (迎风)
if (u(i,j) > 0.0_dp) then
conv_term = conv_term + u(i,j) * (Y_d(i,j) - Y_d(i-1,j)) / dx
else
conv_term = conv_term + u(i,j) * (Y_d(i+1,j) - Y_d(i,j)) / dx
end if
! 径向对流 (迎风)
if (v(i,j) > 0.0_dp) then
conv_term = conv_term + v(i,j) * (Y_d(i,j) - Y_d(i,j-1)) / dr
else
conv_term = conv_term + v(i,j) * (Y_d(i,j+1) - Y_d(i,j)) / dr
end if
! 添加轴对称项
conv_term = conv_term + v(i,j) * Y_d(i,j) * area_factor
! 更新质量分数
Y_d(i,j) = Y_d(i,j) - dt * (conv_term - mv_growth)
! 更新液滴半径 (9)
if (N_d(i,j) > 1e-10_dp) then
r_d(i,j) = (3.0_dp * Y_d(i,j) * rho_g(i,j) / (4.0_dp * pi * rho_l * N_d(i,j)))**(1.0_dp/3.0_dp)
end if
! 计算成核和生长源项
mean_free_path = k_B * T_g(i,j) / (sqrt(2.0_dp)*pi*(mol_diameter**2)*p(i,j))
T_C = T_g(i,j) - 273.15_dp
S = p(i,j) / (611.2_dp * exp(22.46_dp*T_C/(T_C + 272.62_dp)))
R_c = (2.0_dp*sigma) / (rho_l * R_gas * T_g(i,j) * log(S))
J = 1e-3_dp * (p(i,j)/(R_gas*T_g(i,j)))**2 / rho_l * &
sqrt(2.0_dp*sigma/(pi*(18e-3_dp/NA)**3)) * &
exp(-4.0_dp*pi*R_c**2*sigma/(3.0_dp*k_B*T_g(i,j)))
Delta_T = (T_g(i,j)**2 * R_v / h_tg) * log(S)
Kn = mean_free_path / r_d(i,j)
drdt = (lambda_v*Delta_T*(1.0_dp - R_c/r_d(i,j))) / &
(p(i,j)*h_tg*r_d(i,j)*(1.0_dp + &
(2.0_dp/(1.5_dp*Pr))*(r_d(i,j)/(1.0_dp+lambda))*Kn))
mv_nucleation = J * (4.0_dp/3.0_dp)*pi*R_c**3*rho_l
mv_growth = N_d(i,j) * 4.0_dp*pi*r_d(i,j)**2 * drdt * rho_l
! 设置源项 (4-6)
S_m(i,j) = -mv_growth
S_u(i,j) = u(i,j) * S_m(i,j)
S_v(i,j) = v(i,j) * S_m(i,j)
S_h(i,j) = h_tg * S_m(i,j)
end do
end do
! 2. 气相求解 (有限差分法)
do j = 2, nr-1
do i = 2, nx-1
! 计算位置坐标
x_pos = (i - 0.5_dp) * dx
r_pos = (j - 0.5_dp) * dr
area_factor = 1.0_dp / r_pos ! 轴对称因子
! 连续性方程 (1)
dudx = (u(i+1,j) - u(i-1,j)) / (2.0_dp*dx)
dvdy = (v(i,j+1) - v(i,j-1)) / (2.0_dp*dr)
div_u = dudx + dvdy + v(i,j)*area_factor
rho_g(i,j) = rho_g(i,j) - dt * (rho_g(i,j)*div_u + u(i,j)*(rho_g(i+1,j)-rho_g(i-1,j))/(2.0_dp*dx) &
+ v(i,j)*(rho_g(i,j+1)-rho_g(i,j-1))/(2.0_dp*dr)) + dt*S_m(i,j)
! x方向动量方程 (2)
dudx = (u(i+1,j) - u(i-1,j)) / (2.0_dp*dx)
dudy = (u(i,j+1) - u(i,j-1)) / (2.0_dp*dr)
dpdx = (p(i+1,j) - p(i-1,j)) / (2.0_dp*dx)
laplace_u = (u(i+1,j) - 2.0_dp*u(i,j) + u(i-1,j))/(dx*dx) + &
(u(i,j+1) - 2.0_dp*u(i,j) + u(i,j-1))/(dr*dr)
u(i,j) = u(i,j) - dt * (u(i,j)*dudx + v(i,j)*dudy + (1.0_dp/rho_g(i,j))*dpdx) &
+ dt * (mu/rho_g(i,j)) * (laplace_u + (1.0_dp/3.0_dp)*dudx) &
+ dt * S_u(i,j)/rho_g(i,j)
! r方向动量方程 (2) - 轴对称形式
dvdx = (v(i+1,j) - v(i-1,j)) / (2.0_dp*dx)
dvdy = (v(i,j+1) - v(i,j-1)) / (2.0_dp*dr)
dpdy = (p(i,j+1) - p(i,j-1)) / (2.0_dp*dr)
laplace_v = (v(i+1,j) - 2.0_dp*v(i,j) + v(i-1,j))/(dx*dx) + &
(v(i,j+1) - 2.0_dp*v(i,j) + v(i,j-1))/(dr*dr)
v(i,j) = v(i,j) - dt * (u(i,j)*dvdx + v(i,j)*dvdy + (1.0_dp/rho_g(i,j))*dpdy) &
+ dt * (mu/rho_g(i,j)) * (laplace_v - v(i,j)/(r_pos*r_pos) + (1.0_dp/3.0_dp)*dvdy) &
+ dt * S_v(i,j)/rho_g(i,j)
! 能量方程 (3)
dTdx = (T_g(i+1,j) - T_g(i-1,j)) / (2.0_dp*dx)
dTdy = (T_g(i,j+1) - T_g(i,j-1)) / (2.0_dp*dr)
h_g(i,j) = h_g(i,j) - dt * (u(i,j)*dTdx + v(i,j)*dTdy) * Cp_air &
+ dt * (k_thermal/(rho_g(i,j)*Cp_air)) * &
((T_g(i+1,j)-2.0_dp*T_g(i,j)+T_g(i-1,j))/(dx*dx) + &
((T_g(i,j+1)-2.0_dp*T_g(i,j)+T_g(i,j-1))/(dr*dr)) &
+ dt * S_h(i,j)/rho_g(i,j)
T_g(i,j) = h_g(i,j) / Cp_air
! 状态方程
p(i,j) = rho_g(i,j) * R_gas * T_g(i,j)
end do
end do
! 3. 边界条件
! 入口边界 (i=1) - 射流条件
i = 1
do j = 1, nr
rho_g(i,j) = P_jet / (R_gas * T_jet)
u(i,j) = U_jet
v(i,j) = 0.0_dp
p(i,j) = P_jet
T_g(i,j) = T_jet
if (j <= nr/2) then
N_d(i,j) = 1.0e15_dp
Y_d(i,j) = 0.01_dp
r_d(i,j) = 1.0e-6_dp
else
N_d(i,j) = 1.0e-20_dp
Y_d(i,j) = 0.0_dp
r_d(i,j) = 1.0e-9_dp
end if
end do
! 出口边界 (i=nx) - 零梯度
i = nx
do j = 1, nr
rho_g(i,j) = rho_g(i-1,j)
u(i,j) = u(i-1,j)
v(i,j) = v(i-1,j)
p(i,j) = P_env ! 出口压力设为环境压力
T_g(i,j) = T_g(i-1,j)
N_d(i,j) = N_d(i-1,j)
Y_d(i,j) = Y_d(i-1,j)
r_d(i,j) = r_d(i-1,j)
end do
! 对称轴 (j=1)
j = 1
do i = 1, nx
v(i,j) = 0.0_dp
! 对称条件
rho_g(i,j) = rho_g(i,j+1)
u(i,j) = u(i,j+1)
p(i,j) = p(i,j+1)
T_g(i,j) = T_g(i,j+1)
N_d(i,j) = N_d(i,j+1)
Y_d(i,j) = Y_d(i,j+1)
r_d(i,j) = r_d(i,j+1)
end do
! 外边界 (j=nr) - 环境条件
j = nr
do i = 1, nx
rho_g(i,j) = rho_air
u(i,j) = Ue
v(i,j) = 0.0_dp
p(i,j) = P_env
T_g(i,j) = Te
N_d(i,j) = 1.0e-20_dp
Y_d(i,j) = 0.0_dp
r_d(i,j) = 1.0e-9_dp
end do
! 4. 输出结果
if (mod(step,100) == 0) then
! 输出中心线数据
i = nx/2
j = 1
x_pos = (i - 0.5_dp) * dx
r_pos = (j - 0.5_dp) * dr
T_C = T_g(i,j) - 273.15_dp
S = p(i,j) / (611.2_dp * exp(22.46_dp*T_C/(T_C + 272.62_dp)))
write(10,'(I0,9(",",ES12.5))') step, time, x_pos, r_pos, &
N_d(i,j), Y_d(i,j), r_d(i,j), T_g(i,j), p(i,j), S
print *, "步数:",step," 时间:",time," 温度:",T_g(i,j)," 压力:",p(i,j)
! 输出云图数据 (每500步)
if (mod(step,500) == 0) then
rewind(20)
rewind(30)
write(20, "(A)") "x[m],r[m],T[K]"
write(30, "(A)") "x[m],r[m],P[Pa]"
do j = 1, nr
do i = 1, nx
x_pos = (i - 0.5_dp) * dx
r_pos = (j - 0.5_dp) * dr
write(20, '(3(ES12.5,","))') x_pos, r_pos, T_g(i,j)
write(30, '(3(ES12.5,","))') x_pos, r_pos, p(i,j)
end do
end do
end if
end if
time = time + dt
! 终止条件:超过0.05秒 (高速射流需要更少的时间)
if (time > 0.5_dp) exit
end do
!=============== 结束阶段 ================!
close(10)
close(20)
close(30)
! 输出最终云图
open(unit=40, file="final_temperature.csv", status="replace")
open(unit=50, file="final_pressure.csv", status="replace")
write(40, "(A)") "x[m],r[m],T[K]"
write(50, "(A)") "x[m],r[m],P[Pa]"
do j = 1, nr
do i = 1, nx
x_pos = (i - 0.5_dp) * dx
r_pos = (j - 0.5_dp) * dr
write(40, '(3(ES12.5,","))') x_pos, r_pos, T_g(i,j)
write(50, '(3(ES12.5,","))') x_pos, r_pos, p(i,j)
end do
end do
close(40)
close(50)
print *, "=== 模拟完成 ==="
print "(A,ES10.3)", "最终时间 [s]: ", time
print "(A,ES10.3)", "最大液滴半径 [m]: ", maxval(r_d)
print "(A,ES10.3)", "最高温度 [K]: ", maxval(T_g)
print "(A,ES10.3)", "最低压力 [Pa]: ", minval(p)
print "(A,ES10.3)", "最大速度 [m/s]: ", maxval(u)
end subroutine solve
end module contrail_simulation
!======================== 主程序 ========================
program main
use SchmidtAppleman
use contrail_simulation
implicit none
! 执行凝结尾迹形成判据
call check_contrail_formation()
! 执行综合模拟
call solve()
end program main
最新发布