Fortran计算实例.pptVIP

13. 实例之一 实例之二 * * 1. 建模 如流体力学方程组 要点：作些简化后确定维数，注意诸项之取舍（如热传导、辐射、重力等） 激波管问题 初条 2. 无量纲化 对1D情况 2维流体力学方程组，*号已略去 3. 确定解域、划分网格及设置边界条件 可划分为上千个格点 模拟区域： 网 格： 均匀网格（最好用自适应网格） 边 条： 在激波传播至边界之前有效 4. 选取适当的差分格式 Lax格式或Friedrichs格式 精度： 稳定条件： 模型方程 对流体力学问题，其稳定条件为 5. 编程及调试 program shocktube parameter(im=1001) implicit double precision(a-h,o-z) dimension x(im),d(m),u(im),t(im) dimension d1(m),u1(im),t1(im) d0=1.67d-5 t0=3.d2 r=2.78d-2 v0=dsqrt(r*t0) L0=10 gamma=1.4d0 …… 可将3个变量一起放在一个2维数组 设定常数 d: 密度; u: 速度；t: 温度 dx=2.d0/dble(im-1) do 10 i=1,im x(i)=dx*dble(i-1)-1. continue do 20 i=1,(im+1)/2 d1(i)=1.d0 u1(i)=0.d0 t1(i)=1.d0 do 30 i=(im+1)/2+1,im d1(i)=1.d-1 u1(i)=0.d0 t1(I)=1.d0 划分网格 设定初值 10 20 30 dt=0.9d0*dx continue do 40 i=1,im d (i)= d1 (i) u(i)= u1 (i) t(i)= t1(i) do 50 i=2,im-1 d1(i)=0.5d0*(d(i-1)+d(i+1))-dt*(u(i)*(d(i+1)-d(i-1))/2.d0/dx+ d(i)*(u(i+1)-u(i-1))/2.d0/dx) u1(i)= 0.5d0*(u(i-1)+u(i+1))-dt*(u(i)*(u(i+1)-u(i-1))/2.d0/dx+ (t(i+1)-t(i-1))/2.d0/dx +t(i)/d(i)*(d(i+1)-d(i-1))/2.d0/dx) t1(i)=0.5d0*(t(i-1)+t(i+1))-dt*(u(i)*(t(i+1)-t(i-1))/2.d0/dx+ (gamma-1.d0)*t(i)*(u(i+1)-u(i-1))/2.d0/dx) 初始时间步长 赋新值 999 40 主体部分 + + + 50 可考虑加入人为耗散 tim=0.d0 初始时间 d1 (1)= d (1) u1(1)= u (1) t1(1)= t (1) d1 (im)= d (im) u1(im)= u (im) t1(im)= t (im) tim=tim+dt dt=1.d0 do 60 i=2,im-1 dt=dmin1(dt,dx/(dsqrt(gamma*t1(i))+dabs(u1(i))) Continue write(*,*)tim,dt 由柯朗条件定下一时间步长 固定边界条件 主体部分 60 局地无量纲声速 屏幕输出供调试 if(tim在某些时刻）call output(tim,d1,u1,t1) if(输出结果的次数等于tend/0.02)goto 888 If(dt.le.1.e-9) then输出”dt太小”并goto 888 goto 999 close(10) End Subroutine output(tim,d1,u1,t1) ….. If(tim.eq.0)open(10,…..form=‘unformatted’) write(10)d,u,t Return end 主体部分 如0.02, 0.04, 0.06, 0.08,……。即每0.02时间间隔则将结果输出到数据文件。 停止运算的无量纲时间 可采用另一计时变量，如tp, 其初值为0，tp=tp+dt，待tp快接近0.02，dt=0.02-tp，然后输出结果，并重新令tp=0. 888 磁场扩散 磁场扩散方程(仅1D) 1. 建模 初条 解析解 2. 无量纲化 3. 确定解域、划分网格及设置边界条件 可划分为上千个格点 模拟区域： 网 格： 均匀网格 边 条： 扩散至边界之前有效 4. 选取适当的差分格式 模型方程 全隐格式 精度：