Fortran计算实例培训讲学

可采用另一计时变量，如tp, 其初值为0， tp=tp+dt，待tp快接近0.02，dt=0.02-tp， 然后输出结果，并重新令tp=0.
实例之二
1. 建模
磁场扩散
磁场扩散方程(仅1D)
初条
B(x,0) BB 0， 0， xx00
B 2B
t
解析解 B(x,t)2B 0 x/ 4teu2du 0
2. 无量纲化
B B0B*
x L0x*
t
L2 0
t*
B 2B t
3. 确定解域、划分网格及设置边界条件
B
可划分为上
1
千个格点
x
0
1
模拟区域： 网 格： 边 条：
0x1 x(i)(i1)/N (1)
B n(1 )0 ； B n(N )B 0(N )
扩散至边界之前有效
均匀网格
4. 选取适当的差分格式
b
s
5. 编程及调试
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 ……
t
x cs | u |
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 ……
+
(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+
50 +
(gamma-1.d0)*t(i)*(u(i+1)-u(i-1))/2.d0/dx)
Fortran计算实例
3. 确定解域、划分网格及设置边界条件
可划分为上 千个格点
1
模拟区域： 网 格：
边 条：
0
1x1
x
1
x(i)2 (i 1 )/N ( 1 ) 1 均匀网格（最好用自适应网格）
fn(1 )f0 (1 )； fn(N )f0 (N )，f其 ,u ,T 中
在激波传播至边界之前有效
a, b,c 系数
可将3个变量一起 放在一个2维数组
设定常数
dx=1./dble(im-1)
do 10 i=1,im
x(i)=dx*dble(im-1)
10
continue
do 20 i=2, im
20
b1(i)=1.d0
b1(1)=0.d0
划分网格 设定初值
tim=0.d0
dt=(x(2)-x(1))**2
主 体
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+ 可考虑
Continue write(*,*)tim,dt
屏幕输出供调试
如0.02, 0.04, 0.06, 0.08,……。即每0.02 时间间隔则将结果输出到数据文件。
主
体
部
if(tim在某些时刻）call output(tim,d1,u1,t1)
分
if(输出结果的次数等于tend/0.02)goto 888
4. 选取适当的差分格式
模型方程
u c u 0
t
x
Lax格式或Friedri来自百度文库hs格式
u i n 1 1 2 u i n 1 u i n 1 2 c x tu i n 1 u i n 1
精度：O(t, x2 ,x2) 稳定条件： t x
t
|c |
对流体力学问题，其稳定条件为
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
固定边界条件
局地无量 纲声速
由柯朗条件定 下一时间步长
60
dt=dmin1(dt,dx/(dsqrt(gamma*t1(i))+dabs(u1(i)))
If(dt.le.1.e-9) then输出”dt太小”并goto 888
goto 999
888 close(10) End
停止运算的 无量纲时间
Subroutine output(tim,d1,u1,t1)
…..
If(tim.eq.0)open(10,…..form=‘unformatted’)
write(10)d,u,t Return end
模型方程
u t
2u x 2
1
全隐格式
u i n 1 u i n x t 2u i n 1 1 2 u i n 1 u i n 1 1
精度：O(t,x2) 稳定条件：恒稳
u in 1 1 (2 x2 t)u in 1 u in 1 1 x2 tu in
c
a
d: 密度; u: 速度；t: 温度
可将3个变量一起 放在一个2维数组
设定常数
dx=2.d0/dble(im-1)
do 10 i=1,im
x(i)=dx*dble(i-1)-1.
10
continue
do 20 i=1,(im+1)/2
d1(i)=1.d0
u1(i)=0.d0
20
t1(i)=1.d0
do 30 i=(im+1)/2+1,im
d1(i)=1.d-1
u1(i)=0.d0
30
t1(I)=1.d0
划分网格 设定初值
tim=0.d0
初始时间
dt=0.9d0*dx
初始时间步长
999 continue
do 40 i=1,im
d (i)= d1 (i) u(i)= u1 (i)
赋新值
40
t(i)= t1(i)