Fortran计算实例

null13. 实例之一13. 实例之一1. 建模如流体力学方程组要点：作些简化后确定维数，注意诸项之取舍（如热传导、辐射、重力等）激波管问题初条null2. 无量纲化对1D情况2维流体力学方程组，*号已略去null3. 确定解域、划分网格及设置边界条件可划分为上千个格点模拟区域：网 格：均匀网格（最好用自适应网格）边 条：在激波传播至边界之前有效null4. 选取适当的差分格式Lax格式或Friedrichs格式精度： 稳定条件：模型方程对流体力学问题，其稳定条件为 null5. 编程及调试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: 温度nulldx=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划分网格设定初值102030nulldt=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)初始时间步长赋新值99940主体部分+++50可考虑加入人为耗散tim=0.d0初始时间nulld1 (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局地无量纲声速屏幕输出供调试nullif(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. 建模初条解析解null2. 无量纲化null3. 确定解域、划分网格及设置边界条件可划分为上千个格点模拟区域：网 格：均匀网格边 条：扩散至边界之前有效null4. 选取适当的差分格式模型方程全隐格式精度： 稳定条件：恒稳null5. 编程及调试program diffusion parameter(im=1001) implicit double precision(a-h,o-z) dimension x(im),bb(im) dimension b1(im) Dimension a(im),b(im),c(im),s(im) b0=1.67d-5 t0=1.d5 L0=1.d5 ……可将3个变量一起放在一个2维数组设定常数a, b,c 系数nulldx=1./dble(im-1) do 10 i=1,im x(i)=dx*dble(im-1) continue do 20 i=2, im b1(i)=1.d0 b1(1)=0.d0划分网格设定初值1020nulldt=(x(2)-x(1))**2 continue do 30 i=1,im bb (i)= b1 (i) do 50 i=2,im-1 a(i)=2.d0+dx*dx/dt b(i)= -1.d0 c(i)=-1.d0 s(i)= dx*dx/dt*bb(i) a(1)=1.d0 b(1)=0.d0 c(1)=0.d0 s(1)=0.d0 a(im)=1.d0 b(im)=0.d0 c(im)=0.d0 s(im)=1.d0初始时间步长赋新值99930主体部分50tim=0.d0初始时间计算三对角矩阵系数null此过程由上依次向下的递推过程俗称“追”。此过程由下依次向上的回代过程俗称“赶”。主体部分dt=1.d0 do 60 i=2,im-1 dt=5.d0*dmin1(dt,dx*dx/2.d0) Continue tim=tim+dt write(*,*)tim,dt60暂取步长为显格式的5倍nullif(tim在某些时刻）call output(tim,d1,u1,t1) if(输出结果的次数等于tend/0.02)goto 888 goto 999 close(10) End Subroutine output(tim,d1,u1,t1) ….. If(tim.eq.0)open(10,…..) 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