fixed - 结论 - 计数 - dp

本文探讨了一类矩阵计数问题,重点在于寻找满足特定条件的强联通有向图的数量。通过对图论的理解和数学技巧的应用,给出了详细的解题思路与步骤。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题目大意:给定n,mn,m,问有多少nnn列的矩阵AA,满足:
1)i,j[1,n],Ai,j[0,m)
2)i,j[1,n],p,q0,(Ap)i,j=(Aq)i,j∃i,j∈[1,n],∀p,q≥0,(Ap)i,j=(Aq)i,j。此时称(i,j)(i,j)为矩阵AA的不动点。
请注意,A0=I,即(A0)i,j[i=j](A0)i,j⇐⇒[i=j].
题解:
先把n=1判掉。
首先不要把它当成线性代数题,将其视为邻接矩阵。
那么就是要找到一对点,使得其两点间不同长度路径方案相同。
首先如果这个图不弱连通,那么一定存在一对点,永远到不了。
然后显然这个可以加强为不强联通。

然后过了一会,恩,也就是两个小时吧,意识到这个条件其实是充要条件(我会告诉你其间我还想到了什么exgcd之类的么)。

考虑反证,假设(x,y)(x,y)是其不动点。
1)若xyx≠y,则(Ap)x,y=(A0)x,y=0(Ap)x,y=(A0)x,y=0,但是因为其强联通,所以存在一个时刻xx能到达y
2)否则x=y,(Ap)x,y=(A0)x,x=1x=y,(Ap)x,y=(A0)x,x=1,但这也是不可能的,因为其强联通并且n2n≥2,所以存在一个时刻,xx能走一个与其相连的点zx然后从zz走到x;同时也可以一直走自环,因此方案数就至少是22了。

问题转化为n个点的有向强联通图计数,考虑主旋律的做法(顺便就当时主旋律的题解了吧)。显然自环不影响图的强联通性因此下文暂时不考虑。

先考虑一个有点问题的做法。
fnfn表示答案,gngn表示若干个强联通分量(下文简记SCC)的方案数,也就是gn=ni=1(n1i1)fignign=∑i=1n(n−1i−1)fign−i
考虑计算gngn,一个想法是,考虑所有不合g法的情况,缩点之后一定存在边。
钦定一些点是出度为0的点,然后剩下的点内部和从剩下的点到钦定的点间可以随便连,钦定的点就是g,钦定的点不能连向剩下的点,算出这个从总数中减去:
gn=hnn1i=1(ni)mi(ni)hnigign=hn−∑i=1n−1(ni)mi(n−i)hn−igi
其中hnhn表示nn个点的有向图

但这么算是肯定有问题的,但是我们想知道到底问题在哪里:

考虑真实情况是钦定的点有k个,缩点后有pp个,会被计算几次。

首先若k<n
考虑枚举的gigi只能恰好包含pp的一个非空子集,枚举这个的大小t:,其会被计算pt=1(pt)∑t=1p(pt)这么多次。
我们希望其被计算1次,其实也很简单:给它配上一个容斥系数就可以了:
pt=1(pt)(1)t+1=1∑t=1p(pt)(−1)t+1=1,也就是说这个系数应该配到连通块个数上,即gngn应当定义为,一个每个弱联通分量都是SCC的图,若其有tt个联通块,则对gn贡献是(1)t+1(−1)t+1
考虑在配上系数后k=nk=n的情况,此时由于i<ni<n所以tt不能枚举到p的全集:t=1p1(pt)(1)t+1=(1)p+1
因而最终k<nk<n会被计算11=01−1=0次,k=nk=n会被计算1((1)p+1)=(1)p+11−((−1)p+1)=(−1)p+1次,发现正好就是gngn的定义。
gngn推导出fnfn依然可以用最开始的式子,也就是与那个1−1的系数无关(因为本质上只不过是奇偶分开讨论而已)。
这样就做完了!开森。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<queue>
#include<utility>
#define lint long long
#define p 1000000007
#define gc getchar()
#define mp make_pair
#define fir first
#define sec second
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define clr(a,n) memset(a,0,sizeof(int)*((n)+1))
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
#define N 3010
using namespace std;
typedef pair<int,int> pii;
inline int inn()
{
    int x,ch;while((ch=gc)<'0'||ch>'9');
    x=ch^'0';while((ch=gc)>='0'&&ch<='9')
        x=(x<<1)+(x<<3)+(ch^'0');return x;
}
#define H(x) mi[(x)*((x)-1)]
int f[N],g[N],mi[N*N],c[N][N],h[N],fac[N],facinv[N];
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%p) (k&1)?ans=(lint)ans*x%p:0;return ans; }
inline int prelude(int n,int m)
{
    fac[0]=1;rep(i,1,n) fac[i]=(lint)fac[i-1]*i%p;facinv[n]=fast_pow(fac[n],p-2),mi[0]=1;
    for(int i=n-1;i>=0;i--) facinv[i]=facinv[i+1]*(i+1ll)%p;rep(i,1,n*n) mi[i]=(lint)mi[i-1]*m%p;return 0;
}
int main()
{
    int n=inn(),m=inn();if(n==1&&m==2) return !printf("1\n");
    prelude(n,m),f[1]=g[1]=1;rep(i,1,n) h[i]=H(i);rep(i,0,n) c[i][0]=1;
    rep(i,1,n) rep(j,1,i) c[i][j]=c[i-1][j-1]+c[i-1][j],(c[i][j]>=p?c[i][j]-=p:0);
    rep(i,2,n)
    {
        g[i]=h[i];
        rep(j,1,i-1) g[i]-=(lint)c[i][j]*h[i-j]%p*g[j]%p*mi[j*(i-j)]%p,(g[i]<0?g[i]+=p:0);
        f[i]=g[i];rep(j,1,i-1) f[i]+=(lint)c[i-1][j-1]*f[j]%p*g[i-j]%p,(f[i]>=p?f[i]-=p:0);
    }
    return !printf("%lld\n",(lint)mi[n]*((h[n]-f[n]+p)%p)%p);
}
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
06-03
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值