レイノルズ方程式をプログラミング
数値計算に関することを書いてきました。
本当は流体力学の研究をしてるんです。
（数値計算と二刀流は大変）
流体力学の中にレイノルズ方程式という有名な式があります。（大学の講義でもでてくると思う）
ひょっとしたらテストにも出てくるかも
そこで、レイノルズ方程式とそのサンプルプログラムを紹介したいと思います。
ナビエストークス方程式
レイノルズ方程式を説明する前に、ナビエストークス方程式という流体の式があります。
∂υ/∂t + (υ * ∇)*υ = -(1/ρ) * ∇p + ν∇²*υ + F
こんな式ですね。
この式を解くことが出来れば大体の流体現象が表せると云われてます。
レイノルズ方程式いらないじゃんと思いますがまあ色々あるんですよ（笑）
コンピュータの限界
ナビエストークス方程式は見ての通り複雑なんですよね。
今では、パソコンのメモリーは4G、８G等、Gオーダーが当たり前ですよね。
でも、20年くらい前だとメモリーはMオーダーだったんですよ。
それではナビエストークス方程式は解けない、そこで、ナビエストークスを核として式を簡単にする必要があったんです。
その一つがレイノルズ方程式なんです。
追記
今ではナビエストークス方程式でシミュレーションするのが一般的です。
レイノルズ方程式
レイノルズ方程式を適応するには条件があります。
① 非圧縮性、層流、ニュートン流体、粘度一定。
② 粘性力が強く、慣性力は無視できる。
③ 薄膜である。
④ 壁面と流体はすべらない
です。
そうすると
d/dx((h³/η) * (dp/dx)) = 6U * dh/dx
になります。
高さhが、分かれば圧力を出せます。
この条件に当てはまるものとしては軸受の潤滑等があります。
（流体潤滑、混合潤滑、境界潤滑など等）
この分野ではまだレイノルズ方程式は使われてますね。
次にレイノルズ方程式を中心差分します。
P(i)を求めたいので
h_x = (-h(i-1)+h(i+1))
p_x = (-p(i-1)+p(i+1))
p(i) = -1.5 * U * η * Δx / h(i)**3 * h_x
& +3.0/8.0 / h(i) * h_x * p_x
& + 0.5 * (p(i-1)+p(i+1))
こんな感じになります。
レイノルズ方程式を使った計算結果の例はこんな感じになりますね。
関連記事
最後にレイノルズ方程式のサンプルプログラムを紹介します。
プログラムコード
program raynol
IMPLICIT NONE
real p(0:100000)
real h(0:100000)
real x(0:1000)
real y(0:1000)
real p1(0:10000)
real avg(0:100000)
real RES (-1000:1000000)
real U,eta,delx, h_x, p_x,total,initial pressure
character dum
integer i,m,n,i_min,i_max,n_max
c———————————————————————–
n_max = 1000
i_min = 0
i_max = 20
total=0
c———————————————————————–
OPEN (11,FILE=’./kusabi1.dat’,STATUS=’unknown’)
OPEN (12,FILE=’./kusabi3.dat’,STATUS=’unknown’)
read(11,*) dum,dum
read(12,*) dum,dum
do i = i_min , i_max
read(11,*) x(i), y(i)
read(12,*) h(i), p(i)
end do
close(11)
close(12)
OPEN (13,FILE=’./parameter1.txt’,STATUS=’unknown’)
OPEN (14,FILE=’./initial pressure.txt’,STATUS=’unknown’)
read(13,*) U
read(13,*) eta
read(14,*) initial pressure
close(13)
close(14)
c———————————————————————-
delx = abs( x(1) – x(0) )
do n =1,n_max,+1
do i = i_min+1, i_max-1
h_x = (-h(i-1)+h(i+1))
p_x = (-p(i-1)+p(i+1))
p1(i) = -1.5 * U * eta * delx / h(i)**3 * h_x
& +3.0/8.0 / h(i) * h_x * p_x
& + 0.5 * (p(i-1)+p(i+1))
end do
p1(i_min) = 0
p1(i_max) = 0
RES(n)=-1.0e10
do i = i_min , i_max
RES(n)=AMAX1( RES(n),abs(p(i)-p1(i)))
end do
IF(MOD(N,10).EQ.0) THEN
write(*,*) n, RES(n)
end if
do i = i_min , i_max
p(i)=p1(i)
end do
end do
OPEN (12,FILE=’./pressure.dat’,STATUS=’unknown’)
write(12,'(A)’) ‘ x p ‘
do i = i_min , i_max , +1
write(12,'(2(e16.8e3,x))’)x(i),p(i)
end do
CLOSE(12)
stop
end program