流体力学基礎講座（第２回午後）
名古屋工業大学大学院 創成シミュレーション工学専攻 後藤 俊幸流体数値計算のための準備
5
拡散方程式
時刻t = 0に原点x = 0にインクを１滴、静止流体中に落としたとしよう．するとインクは次第 に周りの流体中に広がって行く．この現象は拡散と呼ばれ、場所(x, y, z)、時刻tにおけるインク の濃度をu(x, y, z, t)とすると方程式 ∂u ∂t = κ∆u (68) に従う．κは拡散係数である．以下では１次元有限空間[0, 1]における拡散現象を扱う．初期条件は u(x, 0) = √ A 4πκt0 exp ( −(x− 1/2)2 4κt0 ) , A = 0.1, t0= 0.001. (69) とおこう．これはx = 1/2に鋭いピークをもつ釣り鐘型をしている．境界条件はu(x, t) = 0 for (x = 0, 1) (Dirichlet condition) (70)
∂u
∂x = 0 for (x = 0, 1) (Neumann condition) (71)
とする．２番目の境界条件は境界で質量流束が0（熱の場合には断熱）である条件である． もし、境界が無限の遠方にあるならば、(68)の解析解は u(x, t) = ∫ ∞ −∞ 1 √ 4πκtexp ( −(x− y)2 4κt ) u(y, 0)dy (72) とあたえられる．従って初期条件が(69)の様にx = x0に鋭いピークを持つ関数はy = x0以外で はほぼ0と見なして良いから(u(y, 0) = δ(y− x0)、ここでδ(x)はディラックのデルタ関数であ る）、u(x, t)は u(x, t) = √1 4πκtexp ( −(x− x0)2 4κt ) ≡ G(x − x0, t) (73) と書ける． いま区間をN + 1個の格子点で分割し、xi = i∆x, i = 0, 1,· · · , N、時間についても同様に tn= n∆t, n = 0, 1,· · ·とする．このときu(x, t)をu(xi, tn) = uni で表す．方程式の離散化は、時 間については前進差分、空間については中心差分を用いると un+1i = uni + κ ∆t ∆x2 ( uni+1− 2uni + uni−1) (74) となる．時間に関して前進差分を用いる方法を単純オイラー法とよぶ．
0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 u x
Evolution of conentration by Euler method (Dirichlet)
tn=0.04n 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 u x
Evolution of conentration by Euler method (Neumann)
tn=0.04n 図9: ディリ ク レ 条 件 の 下 で の 拡 散 の 様 子 ． tn = 0.04n. 図10: ノイマン条件の下での拡散の様子．tn= 0.04n.
5.1
解析例
上の問題に対して得られた結果は、中央部のuは次第に減少し周辺部はだんだん大きくなる．釣 り鐘形状のグラフが次第に平らになっていくのが見て取れる（図9と10参照）．ノイマン条件で は両端で質量流束が0であるのでuの曲線は境界上で水平になっている．そしてuの総量 Q(t) = ∫ 1 0 u(x, t)dx (75) は、方程式を区間[0,1]で積分して dQ dt = ∫ 1 0 ∂u(x, t) ∂t dx = κ ∫ 1 0 ∂2u ∂x2dx = κ ( ∂u ∂x x=0 − ∂u ∂x x=1 ) = 0 (76) となるから、時間について一定である．これは、プログラムのチェックに使える．ディリクレ条件 の場合には両端でu = 0であるので、uは時間とともに減少する．即ちt→ ∞のときu(x, t)→ 0．5.2
プログラム
この方程式の数値計算プログラムは以下のようになる．c Diffusion equation in 1D by Simple Euler
c
implicit real*8 (a-h,o-z) integer pout,gdraw,nx,tmax real*8 kappa
parameter(nx=50, kappa=0.1, dt=0.001, ibs=0)
（x軸の格子点数、拡散係数、 dt、 境界条件のスイッチ） parameter(tmax=400, gdraw=40) （時間ステップ数、結果を書き出す時間幅） dimension u(0:nx),uold(0:nx) c *** パラメータ設定 dx=1.d0/dfloat(nx) dx2i=1.d0/dx/dx pi=4.d0*atan(1.d0) amp=1.d-1 t_init=0.001 c *** 初期条件
x0=dx*(nx/2) do i=0,nx x=i*dx u(i)=amp/sqrt(2.*pi*t_init)*exp(-(x-x0)*(x-x0)/4/kappa/t_init) if(u(i).lt.1.d-20) u(i)=0.d0 end do c c *** 境界条件 if(ibs.eq.0) then c * データファイルのオープン
write(6,*) ’ Dirichlet Condition’ open(10, file=’diffusion_d.dat’)
c * ディリクレ条件
u(0)=0.d0 u(nx)=0.d0 else
write(6,*) ’ Neumann Condition’ open(10, file=’diffusion_n.dat’) c * ノイマン条件 u(0)=u(1) u(nx)=u(nx-1) endif c c *** Write u(x=1/2,t) c t=0.d0 write(6,1001)
1001 format(/,’ Evolution of u(1/2,t) ’, /)
write(6,1000) t,u(nx/2) do i=0,nx
write(10,1000) i*dx, u(i) end do write(10,1000) c c *** 単純オイラー法による時間発展 do n=1,tmax c c * uをu_oldにコピー do i=0,nx uold(i)=u(i) end do c * 時間を１つ進める do i=1,nx-1 diff=kappa*(uold(i+1)-2.d0*uold(i)+uold(i-1))*dx2i u(i)=uold(i)+diff*dt end do c *** 境界条件を課す if(ibs.eq.0) then c * ディリクレ条件 u(0)=0.d0 u(nx)=0.d0 else c * ノイマン条件 u(0)=u(1) u(nx)=u(nx-1) endif c c *** データを書き出す c if(mod(n,gdraw).eq.0) then do i=0,nx
write(10,1000) i*dx, u(i) end do
write(10,1000) c
write(6,1000) t,u(nx/2) endif
end do
1000 format(1x, 1pe10.3, 2x, 1pe12.6)
stop end
また、このgnuplotはgxdiffusion.pltであり、以下の内容である．
ips=0 ：出力先の選択 、ips=0はスクリーン, ips=1はepsfile
boundary=0 ：境界条件のスイッチ
if(ips==1) set terminal postscript eps plus color "Times-Italic" 22
set nologscale x ：線形の目盛り
set nologscale y
set xrange [0.:1.0] ：x軸の範囲
set yrange [0:1.5] set xlabel "x" set ylabel "u"
if(ips==1 && boundary==0) set output "diffusion_d.eps"
if(boundary==0) set title "Evolution of conentration by Euler method (Dirichlet)";\ plot ’diffusion_d.dat’ using 1:2 title ’$t_n=0.04n$’ with lines 1;\
1:2->xに対してyをプロット
pause -1 "Dirichelt condition"
if(ips==1 && boundary==1) set output "diffusion_n.eps"
if(boundary==1) set title "Evolution of conentration by Euler method (Neumann)";\ plot ’diffusion_n.dat’ using 1:2 title ’$t_n=0.04n$’ with lines 1;\
pause -1 "Neumann condition"
5.3
例題
以下の問題を解いてみよう． 1. t = 0でu(x, 0) = 1とする．このとき、拡散方程式をディリクレおよびノイマン境界条件の もとで解く． ノイマン境界条件の場合には、時間とともにuが両端から減少していくことを確認すること． 2. 無限に長い平板が水平においてあり（これをx軸にとる)、この上方（y ≥ 0）に粘性流体 が静止している．流体は無限に広がっているとする．平板が時刻t > 0から瞬間的に速度 U0 = 1で右側に移動しはじめた場合、流体の運動は以下の方程式に従う． ∂u ∂t = ν ∂2u ∂y2, y > 0, t > 0 (77) ここで、u(y, t)はx軸方向の速度であり、ν = µ/ρは流体の動粘性率であり、ここではν = 0.1 にとる．境界条件と初期条件は u(y = 0, t) = U0, u(y, t)−→ 0 y→ ∞のとき (78) u(y, 0) = 0. (79) これはレイリー（Rayleigh）の問題とよばれ、このときの流れをレイリーの流れという．y の無限遠方を計算機上で実現するには-0.2 0 0.2 0.4 0.6 0.8 1 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 y u Rayleigh problem tn=0.04n 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 u x
Evolution of conentration by Euler method (Dirichlet)
tn=0.04n 図11: レイリー問題における速度の時間変化． ν = 0.1．時刻は tn = 0.04n． 図12: 単純オイラー法による解の不安定性． ∆t = 0.002, tn= 0.04n. 少し工夫がいる．ここでは時間を余り長くまでとらないことにして、実質上y = 1のままで 計算することにする．計算結果は図11のようになる．時間とともに平板の速度の影響が流体 内部にまで次第に広がって（拡散）していくのが見て取れる．拡散の目安は、速度がU0e−1 になる距離u(y∗(t), t) = U0e−1 として定められる．理論的にはy∗(t) = √ 4νtとなる．すな わち、拡散は√tに比例して広がって行く．
5.4
安定性
上で用いた単純オイラー法は簡単であるが、数値的に不安定なスキームである．即ち、少しで も誤差があるとそれが時間とともに指数関数的に増大して、数値解が発散したり振動したりして 有用な解が得られなくなる．図12は∆t = 0.002で上の問題を解いたものである． 一般に放物型偏微分方程式における数値的安定性のためには条件 ∆x2 2κ ≤ ∆t (80) を満たす必要がある．上の条件を数学的に導くのは参考文献2, 6) に譲る事にして、直感的な議論 でおおよそのことを理解しよう． いまある時刻tから次の時刻t + ∆tへの解は u(x, t + ∆t) = ∫ ∞ −∞ 1 √ 4π∆texp ( −(x− y)2 4κ∆t ) u(y, t)dy (81) = ∫ ∞ −∞G(x− y, ∆t)u(y, t)dy (82) と書ける．ここで、G(x, ∆t)は(73)で与えられる．積分記号の中にある部分は、時刻tにおいて yとy + dyの区間にある物質量u(y, t)dyが、時間∆tの間にどれほどxに達するかを表している．Dx t=0 t=Dt t=2Dt -Dx<x<Dx 以内にQの大部分が含まれているためには y Dt Dx Dt/4 Dx/2 x = 2k t y 図13: 点 y にあった物質（総量 Q = 1）が時間 ∆t, 2∆tの間に拡散する．ある与えられ た ∆x に対して時間幅が大きすぎると、 区間 [−∆x, ∆x] にあるべき物質がたくさ ん外に漏れ出す。 図14: 時刻 ∆t における [−∆x, ∆x] の区間内に ある物質の総量と、[−∆x/2, ∆x/2] にあ る物質の総量が等しくなる時刻は ∆t/4． 物質の保存法則から Q(t) = ∫ ∞ −∞u(x, t + ∆t)dx =一定 (83) 1 = ∫ ∞ −∞G(x− y, ∆t)dx (84) であることがわかる．2番目の式の差分近似を考えてみると, ∆tと∆x = x− yが十分小さい時に はあまり物質は拡散しておらず、その多くはyの近傍にあって 1≈ 1 2[G(−∆x, ∆t) + 2G(0, ∆t) + G(∆x, ∆t)] ∆x (85) と近似してよいだろう．このところでこの近似が意味をなすためには、区間[−∆x, ∆x]にその大 部分の物質が含まれていなければならない（図13）．これは関数Gが|x − y| = ∆xのところで 十分小さくなっていることを意味するので、(73)の指数部が1より大きくなっていなければなら ない． ∆x2 4κ∆t > 1 (86) ここで係数が4ではなく2になっているが、これはより詳しい解析によるもので、ここでの議論 ではこだわらなくてよい．これが求める数値的安定性の条件である． いま、(86)を満たしてはいるけれども、しかし∆x2 / (2κ∆t)≈ 1となるようなぎりぎりの空間 と時間の格子間隔をとったとしよう．次に、(∆x/2, ∆t/2)としたとすると、(∆x/2)2 / 2κ(∆t/2) = ∆x2 / 2κ2∆t < 1となり、不安定になる．これは、(86)からみて空間格子幅が半分になったとき、時 間幅が∆t/2では大きすぎて、[−∆x/2, ∆x/2]の領域にそのほとんどが含まれているはずの物質が、 すでにこの領域から漏れ出てしまっていることを意味する．すなわち、この時間幅では物質の保存法 則が満たされない（図14）．保存法則を満たすようにするには、時間幅をもっと小さいものにして、物
質があまり遠くまでたくさん拡散してないようにする必要がある．これが、安定性条件についての一 つの物理的な見方である．上の問題では∆x = 1/50 = 0.02, ∆x2/(2κ) = 4×10−4/2/0.1 = 2×10−3 である．従って∆t = 0.001の時には計算ができ、∆t = 0.002の時には不安定になっている． このような不安定性を回避するために、いくつかの安定なスキームが考案されている．ひとつ はクランク・ニコルソン法（Crank-Nicolson Scheme）と呼ばれるものである．これは、拡散項を 時刻nとn + 1での平均でおきかえたものである． un+1i = uni + κ ∆t 2∆x2 ( un+1i+1 − 2un+1i + un+1i−1 + uni+1− 2uni + uni−1 ) (87) この式は右辺にun+1i が入るため、un+1について３重対角行列を係数とする連立方程式を解かな ければならない．このような解法を陰解法（Implicit Scheme）とよぶ．３重対角行列をひっくり 返すのは難しくはないが面倒である．しかし、数値的安定性は良いので、多少手間がかかっても 時間刻み幅を大きくとれるのでやる価値がある． 式(87)を整理して行列の形にすると、 2 + b −1 0 · · · 0 0 −1 2 + b −1 0 0 · · · 0 0 −1 2 + b −1 0 · · · ... .. . ... . .. . .. ... . .. 0 0 0 −1 2 + b −1 0 0 · · · 0 0 −1 2 + b un+10 un+11 .. . un+1N−1 un+1N = f0n f1n .. . fNn−1 fn N (88) fin= uni+1− (2 − b)uni + uni−1, b = 2∆x 2 κ∆t . (89) となる．３重対角なので直接解く方が反復解法よりも効率がよい．これを解くには漸化式を用い る．今、解を un+1i−1 = αi−1un+1i + βi−1 (90) とおいて式(88)に代入する．結果は un+1i = 1 2 + b− αi−1u n+1 i+1 + fin+ βi−1 2 + b− αi−1 (91) となるので、式(90)と見比べれば、係数αi, βiについての漸化式が得られる． αi = 1 2 + b− αi−1, βi = fin+ βi−1 2 + b− αi−1 (92) ここで境界条件を考える．Dirichlet条件u0 = uN = 0の時には un+10 = α0un+11 + β0 = 0−→ α0 = 0, β0= 0 (93) となる．したがってi = 0からi = N− 1まで式(92)によって順にαi, βi を計算する（前進消去）． その後un+1N = 0を用いればun+1N−1 = αN−1un+1N + βN−1より順次un+1i−1 が計算できることになる （後退代入）．このプロセスはGaussの消去法のプロセスそのものである．条件(80)を満たさない 大きな時間幅∆t = 0.004で計算して得られた解を図15に示す．安定な解が正しく得られている ことが分かる．面倒であるが、安定性が良くなったおかげで時間幅を4倍にまで大きくすること ができ、約４倍程度計算効率が向上している．
0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 u x
Evolution of conentration by Crank-Nicolson (Dirichlet) dt=0.004
tn=0.04n 図 15: クランク・ニコルソン法による解．∆t = 0.004, tn= 0.04n.
6
移流方程式
流れの方程式では、文字通りある速度でもってある物理量を流すという作用をもった項が現れ る．代表的なものは ∂u ∂t + c ∂u ∂x = 0, (94) という形をしている．波動方程式のところでみたように、この方程式の解もまた任意の関数f (x) を用いてそれぞれu(x, t) = f (x− ct)と表されることは代入してみれば確かめられる．それは、 ξ = x− ctとおくとu(x, t) = f (ξ)だから ∂f ∂t = df dξ ∂ξ ∂t =−c df dξ, ∂f ∂x = df dξ ∂ξ ∂x = df dξ, ∂u ∂t + c ∂u ∂x = df dξ(−c + c) = 0 (95) となるからである．このことからξの値が変化しない(x, t)の組は一次関数ξ = x− ctで与えら れ、この直線上ではf (ξ) = f (x− ct) = u(x, t)が一定であることがわかる．このことは時刻t = 0 でu(x, 0) = f (x)と与えられるなら、後の時刻t > 0でのu(x, t)はf (x)をctだけ正の側に平行 移動して得られることを意味している．即ち、方程式(94)は速度cで初期状態を変形せずに移動 （流）して行く現象を記述している（図6参照）．6.1
移流方程式の数値解
これからみると、移流方程式の解の理解はとても易しい．しかし、その数値計算となるとそれ ほど簡単ではない．まずは時間に関して前進差分、空間については中心差分をとってみると、 un+1j − unj ∆t =−c unj+1− unj−1 2∆x (96) となり、これを数値積分してみよう．初期条件は Type A u(x, 0) = { 1 f or |x − x0| < 0.05 0 otherwise (97)t x x-ct=0 0 t1 t2 A B C D x-ct=ξ-1 x-ct=ξ-2 図 16: 移流方程式の解．x− ct = ξ（一定）の直線上では u(x, t) は時刻 t = 0での初期値と同じ値をとる． である．境界条件は、変動が境界まで達する前に計算を終えるので、ここでは気にしなくともよ い．結果は図17である．初期に台形をしていた関数が、時間の経過とともに右に移動して行くが、 激しい振動をともなっており、これでは役にたたない． なぜこのようなことになったのだろうか？ それを理解するには移流方程式の物理的内容に立ち 返る必要がある．速度はc > 0であるから、全てはx軸の負から正の方向に流れる．即ち、全ての 変動は左から伝搬し、負の方向に向かって伝搬することは起こらない．従って、空間微分c∂u/∂x は点xj において変動は左から流れてくる（移流）と考えて、 c∂u ∂x = c uj− uj−1 ∆x , c > 0 (98) ととる方が理にかなっている（図18参照)．また、もしc < 0であるならば、全てはx軸の負の方 向に向かって変動は流れていくので c∂u ∂x = c uj+1− uj ∆x , c < 0 (99) ととる．このように、流す速度cの正負に応じて、常に考えている点の上流（風上）側での差分を とるようなスキームを上流（風上）差分スキームという．方程式(94)の上流差分式は un+1j − unj ∆t =−c unj − unj−1 ∆x , c > 0 (100) un+1j − unj ∆t =−c unj+1− unj ∆x , c < 0 (101) となる．しかし、cの正負をまとめて、 un+1j − un j ∆t =−c unj+1− unj−1 2∆x + |c|∆x 2 unj+1− 2unj + unj−1 ∆x2 (102) と表す方が便利であろう．右辺第２項は拡散係数が |c|∆x2 であたえられる拡散項に対応している． 即ち、上流差分は移流項に中心差分と小さい数値拡散係数を持った擬似的な拡散項を加えたもの と理解できる．図19は(100)による結果である．初期の台形が時間とともに正の側に流されてお り、激しい振動はみられない．しかし、台形の角は丸くなりしかも横に伸びて来ている．
-0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 u x
j-1
j
j+1
n+1
n+2
n
c
∆
t
Flow direction 図17: 移流方程式の解．Type A の初期条件で、 空間に中心差分 (96) を用いた結果．激 しい振動現象を伴った数値解得られる． 図18: 流れが正の方向に向いているとき、時刻 n + 1、点 xjには時刻 n、点 x− c∆t で の変動が速度 c の流れにのって運ばれて くる．6.2
例題
1. 上流差分法の例題プログラムはc > 0の場合について書いてある．これをc = 2,−1の場合 について解いて、どのように振る舞うかを調べよう． 2. 上の問題でc < 0の時にも正しく計算できるようにプログラムを修正し、実際に計算してみ よう． 方程式(94)の両辺を区間[0, L]にわたって積分すると ∂ ∂t ∫ L 0 u(x, t)dx =−c ∫ L 0 ∂u ∂xdx =−c(u(L, t) − u(0, t)) (103) となる．いま境界条件として周期境界条件u(0, t) = u(L, t)あるいは境界上で0となることを仮定 すれば(102)の右辺は0である．即ち ∂Q(t) ∂t = 0, Q(t) = ∫ L 0 u(x, t)dx =一定 (104) となる．これは簡単な計算のチェックとして使える．6.3
ラックス・ベンドルフの方法
上でみた中心差分は余り芳しい結果を生まなかった．これは、流れに方向性があるにも関わら ず空間微分を左右対称にとったため、小さな誤差を拡大してしまったからである．しかし、小さ な誤差が増幅しないように擬似的な減衰作用を付け加えることにより、安定化させて使えるよう にすることができる．以下に示すのはラックス・ベンドルフの方法(Lax-Wendroff)とよばれるも のである． un+1j − unj ∆t =−c unj+1− unj−1 2∆x + c2∆t 2∆x2 ( unj+1− 2unj + unj−1 ) (105)-0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 u x -0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 u x 図19: Type Aの初期条件で、空間に上流差分 (100)を用いた結果． 図20: Type Aの初期条件でラックス・ベンド ルフの方法 (101) を用いた結果．細かい 振動は抑制されている． -0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 u x -0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 u x 図21: Type Bの初期条件で中心差分による結 果． 図22: Type B の初期条件で上流差分の方法に よる結果． -0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 u x 図 23: 移流方程式の解．Type B の初期条件でラックス・ベンドルフの 方法結果．
最後の項は擬似的な拡散項c2∆t/2∂2u/∂x2に対応しており、数値拡散係数はc2∆t/2である．こ のことにより、細かい振動は抑制される．図20は同じ条件の下での計算結果である． これまでみると、上流差分に軍配が上がりそうであるが、これはType Aのような急激に変化す る衝撃波構造を持つかなり極端な初期条件についての比較である．通常の非圧縮流れでは、この ように極端な変化が作られることは少ない．そこで、Type Bのような比較的おとなしい場合につ いて３つのスキームでの計算結果を図21, 22, 6.2に示す． Type B u(x, 0) = 1 10√2πσ2 exp ( −(x− x0)2 σ2 ) , σ = 0.05 (106) いずれも、激しい振動はみられずおとなしく右に移動している．しかし、中心差分法では、時間 の経過とともに山の頂上が増大するとともに幅が狭まり、上流差分法では減衰と幅の広がりがみ られるのに対し、ラックス・ベンドルフ法では山の頂上の値と幅はほぼ一定のまま移動している． このようにみていくと、対象となる流れ場の力学的性質と予想される場の変化の様子を考えに 入れながら、数値計算スキームを選択する必要があることが理解される．もちろん、どんな流れ 場についても正しく解ける頑丈なスキームがあれば望ましいが、実際にはがまの油のような万能 薬の様なスキームはそう簡単には手に入らない．
6.4
プログラム
ラックス・ベンドルフの方法による数値計算プログラムは以下のようになる．c Convective equation in 1D by Lax-Wendroff or Central Difference
c
c *** Choice of the numerical scheme 差分法の選択の解説
c switch=0: Central difference
c switch=1: Lax-Wendroff
c
c *** Choice of the initial condition 初期条件の選択の解説
c init_type=0: Top hat
c init_type=1: Gaussian
c
implicit real*8 (a-h,o-z)
integer kint,kdraw,nx,nt,tmax,switch
c parameter(nx=100, tmax=40, kint=2, kdraw=20)
parameter(nx=100, tmax=100, kint=2, kdraw=20)
初期条件台形、中心差分を選択
parameter(c=1.d0, dt=0.002, width=0.05, init_type=0, switch=0) dimension u(0:nx),uold(0:nx) character*20 file(0:2) c *** Make parameters パラメータ設定 dx=1.d0/dfloat(nx) pi=4.d0*atan(1.d0) sigma2=width*width amp=1.d-1/sqrt(2.*pi*sigma2) r=c*dt/dx file(0)=’ Central_difference’ file(1)=’ Lax-Wendroff ’ c 出力ファイルを設定 if(switch==0) open(10,file=’Central_difference.dat’) if(switch==1) open(10,file=’Lax-Wendroff.dat’)
write(6,1000) file(switch),r
1000 format(1x,A20,1x,’Courant Number=’,1pe10.3,/)
write(6,*) ’ Conservation of Q’
c 初期条件を設定
c 保存量を計算（計算チェックのため）
c *** Initial condition and
c computation of total Q(t)=\int u(x,t)dx ***
x0=dx*(nx/2) do i=0,nx
x=i*dx u(i)=0.d0
if(init_type==0 .and. abs(x-x0)<=width) u(i)=1.d0 if(init_type==1) u(i)=amp*exp(-0.5d0*(x-x0)**2/sigma2) end do
c
c *** Output data for t=0
c
call output(u,dx,dt,nx,0,kint,kdraw) c
c *** Time advancing by simple Euler
do nt=1,tmax c c *** Copy u to uold do i=0,nx uold(i)=u(i) end do c
c *** Choice of the numerical scheme
c switch=1: Lax-Wendroff
c switch=0: Central difference
c 差分法の選択 do i=1,nx iplus=mod(i+1,nx) utmp=uold(i)-0.5d0*r*(uold(iplus)-uold(i-1)) diff=0.5d0*r*r*(uold(iplus)-2.d0*uold(i)+uold(i-1)) u(i)=utmp+switch*diff end do c
c *** Enforce periodic boundary condition
c 境界条件を課す
u(0)=u(nx) c
c *** Output data for t>0
c call output(u,dx,dt,nx,nt,kint,kdraw) end do stop end c ******************************************************************** subroutine output(u,dx,dt,nx,nt,kint,kdraw)
implicit real*8 (a-h,o-z) real*8 u(0:nx) save Q_init c **** write Q if(mod(nt,kint).eq.0) then Q=0.d0 do i=0,nx Q=Q+u(i) end do Q=Q*dx if(nt==0) Q_init=Q write(6,1000) nt*dt,Q/Q_init endif
c *** write u(x,t) if(mod(nt,kdraw).eq.0) then do i=0,nx write(10,2000) i*dx,u(i) end do write(10,2000) endif 1000 format(1x,2(1pe10.3,1x)) 2000 format(1x,2(1pe12.5,1x)) return end
