第 12 章 非圧縮性流れの解法|渦度輸送方程式と流れ関数方程式を解く方法
1960 年代のはじめに，Fromm は流れに平行に置かれた角柱まわりの非定常２次元平面流れをコンピュー
タで解析，角柱から発生する渦列の様子を可視化しアニメーションを作っている．この先駆的成果は，日本のテレビでも，コンピュータ・アプ リケーションの可能性を示す一例として紹介された．その解析には後に Fromm の方法と呼ばれる渦度輸送方程式と流れ関数方程式を解く方法 (vorticity transport equationand streamfunction equation method) のひとつが用いられていた．流れ関数は定常または非定常の２次元

流れ，または定常２次元圧縮性流れに対して定義できるもので，非定常２次元圧縮性流れ，３次元流れを扱うことはできない．この章には，定常または非定常の２次元非圧縮性の層流を解析するための，渦度輸送方

程式と流れ関数方程式を解く方法，すなわち渦度−流れ関数法 (!

;  法または  ; ! 法) について説明する．またこの種の解法の問題点とその対策についても述べる．

12.1

基礎方程式

非圧縮性流れの基礎方程式は次の Navier-Stokes(NS) 方程式と連続方程式である．

duu @uu
 +u ru = ;rp +  r2u
dt @t
r u = 0(12.1)(12.2)

 @=@t +u r は実質微分 (substantial derivative)，諸量は無次元化されており，u は流速，pu=dt = @uu=@t+ru 2=2;u  なるは静圧/密度， = 1=R，R は Reynolds 数である．NS 方程式 (12.1) は du
関係1を用い，渦度   ru を導入すれば
ただし d=dt@uu+ ru 2=2 ;u  = ;rp +  r2u@t(12.3)r(r) = 0 r(u  ) = (  r) u ;(u  r)  なる関係2を用いれば次となる．更にこの式の回転を取り，の渦度輸送方程式 (vorticity transport equation) が導かれる．

d= (  r) u +  r2
dt
なおこの式の左辺は渦度  の輸送項，右辺第 1 項はその発生項，第 2 項は拡散項である．
渦度−流れ関数法の基礎方程式を導く．渦度輸送方程式は２次元流れの場合には次のようになる．

d! @!
 @t +u  r! =  r2 !
dt



あるいは

; 
r   rr ; rrr

@2! @2! @! @! @!
+
u +v = +@t@x @y@x2 @y2rrr

r  ; r
r r
r
r r

(12.4)
r

1 ベクトル三重積 a (b c) = (c a) b (a b) c から u (u) = (u0 u) (u) u ，右辺第 1 項の下添字 0 はが u0 には演算されないことを示す．したがって (u 0 u ) u u + v v + w w = u 2=2 となる．
2 ベクトル三重積から(u  ) = ( ) u ( u )  ，右辺第 1 項の は  だけでなく u の演算子でもあることに注意すべき ) u = () u +u ( )，また=(u ) = 0．右辺第 2 項についても同様のことが言える．
である．したがって (

r

r

1

2

非圧縮性流れの解法 | !

;法ただし ! は渦度の紙面に垂直な成分で

!=

@v @u
;
@x @y

(12.5)

である．この式には発生項がなく，２次元流れの渦度は対流する間に生成消滅せず，拡散により最大値が減
少し最小値が増加し一様化していくことになる．ここで(x ) = (x0 ) +Zxi  u  dxx
zx0

ただし i z は紙面に垂直な単位ベクトル，あるいは

@
= u
@y

@
= ;v
@x

(12.6)

で定義される流れ関数 (streamfunction)  を導入する．２次元流れの渦度の式 (12.5) に流れ関数の定義式(12.6) を代入すれば，次の流れ関数の方程式 (streamfunction equation) が導かれる．

r2   @@x2 + @@y2 = ;!
22

(12.7)

渦度輸送方程式と流れ関数方程式を解く方法は，基本的には非定常２次元流れの解法で，解を時間ステップ

t だけ前進させるための計算手順は次のようになる．渦度輸送方程式 (12.4) の初期値問題から ! を計算，
流れ関数方程式 (12.7) の境界値問題から  を計算，式 (12.6) から流速成分 u v を計算，渦度の境界値 !
Bを計算．また定常流れは時間の十分経過した時点での収束解として求められる．

12.2

定常流れの解法 (Dennis-Chang 法)

Dennis-Chang 法は，直接 定常流れの渦度輸送方程式を解くものである．定常流れの場合には式 (12.4) はu r! =  r2 !

u

あるいは

 @2! @2! 
@! @!
+
v =
+
@x @y
@x2 @y2(12.8)になる．左辺の対流項は，1 次上流差分を用いれば等間隔長方形格子では

u ;ju j ! +1 ; !
u + ju j ! ; ! ;1
+
2
x
2
x
v ;jv j ! +1 ; !
v + jv j ! ; ! ;1
+ 2
+ 2
y
y

u r!  C (!) 
UD

ij

ij

ij

i

ij

j

ij

ij

ij

ij

ij

ij

i

j

ij

ij

ij

ij

(12.9)

また不等間隔長方形格子では

u ;ju j ! +1 ; !
+ u +2ju j !x; ! ;1
2
x +1 2
;1 2
v ;jv j ! +1 ; !
v + jv j ! ; ! ;1
+ 2
+ 2
y +1 2
y ;1 2u r!  C (!) 
UD

ij

i

ij

j
i

ij

ij

ij

ij

ij

ij

=

ij

i

ij

j

i

ij

ij

=

ij

j

=

ij

j

=

(12.10)となる．また 2 次中心差分を用いれば等間隔長方形格子では; (u!) ;1 + (v!) +1 ; (v!) ;1r (u !)  C (!)  (u!) +1 2x
2y
i

CD

j

i

j

ij

ij

(12.11)

Dennis - Chang 法

3

また不等間隔長方形格子ではu!)  CCD (!)

r (



 xi+1=2 xi;1=2 (u! )i+1j ; (u! )ij +(u! )ij ; (u! )i;1jxi;1=2 + xi+1=2 xi+1=2
xi;1=2yj ;1=2
yj +1=2 
1+(v! )ij +1 ; (v! )ij +(v! )ij ; (v! )ij ;1
yj ;1=2 + yj +1=2 yj +1=2
yj ;1=2
1

(12.12)となる．なお xi+1=2 = xi+1 ; xi である．2 次中心差分の式は保存形になっている．また式 (12.8) の右辺の拡散項は等間隔長方形格子では

 r 2 !  D (! ) 

1

 1 ;

R x



2 !i;1j ; 2!ij + i+1j +不等間隔長方形格子では

 r 2 !  D (! ) 1y

2 !ij ;1 ; 2!ij + ij +1



!ij ; !i;1j 
R xi;1=2 + xi+1=2 xi+1=2
xi;1=22!ij +1 ; !ij !ij ; !ij ;1 +;yj ;1=2 + yj +1=2 yj +1=2
yj ;1=2
1

 !i+1j !ij

;

2

;

(12.13);(12.14)となる．良く知られているように，1 次上流差分の式は十分安定であるが精度不足で，他方 2 次中心差分の式はレ イノルズ数が小さくないときには不安定性のため計算できなくなる．Dennis-Chang 法は次式により
32 次精度の解を安定に得るものである ．

CUD (!(n+1) ) ; D(!(n+1) ) = CUD (!(n) ) ; CCD (!(n) )(n)(12.15)(n+1)

この差分方程式は右辺に既知の !ij を入れ，!ijの連立 1 次方程式として解かれる．したがって計算は
(n+1)
(n)
1 次上流差分の式の安定さで行われ，解が収束し !ij= !ij になれば，両辺の 1 次上流差分の項は相殺されるので 2 次中心差分の式を解いたことになる．
次にこの解法の残りの部分すなわち流れ関数  ，流速成分 u v ，渦度の境界値 !B の計算について説明する．流れ関数の方程式 (12.7) の差分方程式は等間隔長方形格子の場合には次のようになる．
1

x2

;



i;1j ; 2ij + i+1j +

1

y2

;



ij ;1 ; 2ij + ij +1 = ;!ij(12.16)

また不等間隔長方形格子の場合には次のようになる．

ij ; i;1j xi;1=2 + xi+1=2 xi+1=2
xi;1=2ij +1 ; ij ij ; ij ;1 
2;= ;!ij+yj ;1=2 + yj +1=2 yj +1=2
yj ;1=22

 i+1j ij;;(12.17)

これらの方程式の境界値問題は大形疎行列の連立 1 次方程式になり，流れ関数 ij はこの連立 1 次方程式を解くことによって求められる．また流速 u v は次の式 (12.6) の差分式から求められる．等間隔長方形格子の場合には

uij = (ij +1 ; ij ;1 )=2y 
vij = ;(i+1j ; i;1j )=2x(12.18a)(12.18b)

3 Dennis, S.C.R. and Chang, G.Z., Numerical solutions for steady ow past a circular cylinder at Reynolds numbers up
to 100, J. Fluid Mech., Vol.42(1970), Part 3, 471{489.

4

非圧縮性流れの解法 | !

;法また不等間隔長方形格子の場合には

yj ;1=2y1(ij +1 ; ij ) + j +1=2 (ij ; ij ;1 ) yj ;1=2 + yj +1=2 yj +1=2
yj ;1=2xi;1=2
x
1
vij = ;(
i+1j ; ij ) + i+1=2 (ij ; i;1j )xi;1=2 + xi+1=2 xi+1=2
xi;1=2(12.19a)uij =(12.19b)

最後に，固体壁境界における渦度の値は次の Woods の条件から求められる．y = y0 境界がそれ自身に平行に速度 U = (y )i0 で運動している場合には

3 U + 3 ( ;  ) + O(y2 )1(12.20)2yy2 i0 i1
この式は Taylor 展開
11
i1 = i0 + y(y )i0 + y2 (yy )i0 + y3 (yyy )i0 +O(y4 )2!
3!
に (y )i0 = U ，式 (12.7) から得られる関係 (yy )i0 = ;!i0 ，これを微分した (yyy )i0 = ;(!i1 ;!i0 )=y +O(y) を代入し整理すれば導かれる．またほかの Woods の条件は次のようになる．
1
3 V + 3 ( ;  ) + O(x2 )!0j = ; !1j ;2x
x2 0j 1j1
3 V + 3 ( ; !Ij = ; !I ;1j +) + O(x2 )2
x
x2 Ij I ;1j1
3
3!iJ = ; !iJ ;1 ; U + 2 (iJ ; iJ ;1 ) + O(y2 )2yyDennis-Chang 法の計算の手順は，基本的には渦度の式 (12.15) とその境界条件 (12.20) から ! を求める
計算と，流れ関数の式 (12.7) とその境界条件から  を求める計算を反復するものである．その際に次式に!i0 = ; !i1 +より ! (n+1) に強いダンピングをかけることが必要である．

!ij(n+1) = !ij(n) + (!ij(n+1) ; !ij(n) )(12.21)

ここに  はダンピング係数でその適当な値は例えば 0.25 である．つまり前回の値 !ij をわずかに修正し(n)たところに先の解があるということで，ダンピングをかけなければ ! はたちまち発散してしまう．数値計算が行われる以前に また行われるようになってからも，各種の流れが渦度輸送方程式と流れ関数方程式を摂動法で解くことによって求められている．しかしながら摂動法では第 3 近似値までつまり n = 3 までが精々で，またダンピングを掛けるという考えは全くなかった．剥離を伴う流れの計算では，n = 50 程度で漸く剥離域が現れ，解の収束には更に多くの反復が必要である．摂動法の論文には奇妙で無意味な結果が数多く見られる．
ここでひとつの計算例として正方形キャビティ内の循環流れを解くことにする．この流れは解法の精度の検証のために広く用いられてきたものである．正方形キャビティの一辺の長さは 1 で，上壁面 (y = 1) は速

度 u = 1 で運動し，残る３つの壁面 (x = 0 1 y = 0) は固定されている．計算格子は 20

20 の不等間隔長方形格子である．次にこの流れを求める FORTRAN 95/90 Free source form プログラムを示す．
PROGRAM MAIN!**************************************************************! Problem: Square Cavity Flow
! Numerical Method: omega-psi Eqns, Dennis-Chang Method, SOL!**************************************************************

Dennis - Chang 法

5

PARAMETER(if=20,jf=20)
DIMENSION x(0:if),y(0:jf),h(0:if),g(0:jf),c(if,jf,8),omg(0:if,0:jf),psi(0:if,0:jf), &u(0:if,0:jf),v(0:if,0:jf),q(0:if,0:jf),iq(0:if)
CHARACTER*15 z1
COMMON Re,na,resomg,respsi
DATA psi,u,v/441*0.,441*0.,441*0./ Re,naf/500.,1000
CALL GRIGEN(x,y,h,g,c,if,jf)!grid generation
OPEN(20,FILE='OUTPUT.dat')
WRITE(20,'(A8, 2X 21F7.3)')'x or y =',x
FORALL(i=0:if)u(i,jf)=1.
na=0 10 na=na+1
CALL CALOMG(h,g,c,omg,psi,u,v,if,jf)!calculation of vorticity
CALL CALPSI(c,omg,psi,u,v,if,jf)!calculation of stream function
CALL CPU_TIME(sec)WRITE(20,60)na,resomg,respsi,sec
60 FORMAT(1H 'na =', I3, 2X 'resomg =' F7.3, 2X 'respsi =', F7.3, 2X 'cputime =', F5.2)
IF(resomg/100.+respsi>5.E-4.AND.na<naf)
GOTO 10
FORALL(i=0:if)iq(i)=i
DO k=1,4
SELECT CASE(k)CASE(1) z1='
vorticity' FORALL(i=0:if,j=0:jf)q(i,j)=omg(i,j)
CASE(2) z1='stream function' FORALL(i=0:if,j=0:jf)q(i,j)=psi(i,j)
CASE(3) z1=' velocity u ' FORALL(i=0:if,j=0:jf)q(i,j)=u(i,j)
CASE(4) z1=' velocity v ' FORALL(i=0:if,j=0:jf)q(i,j)=v(i,j)
ENDSELECTWRITE(20,'(//1H 5X A6,A15,A6/)')'***** ', z1, ' *****'
DO j=jf,0,-1
WRITE(20,'(I3, 21F9.3)')j,(q(i,j),i=0,if)
ENDDO
WRITE(20,'(21I9 )')iq
ENDDO
CLOSE(20)
CALL SCGRAPH(x,y,omg,psi,if,jf)
STOP
END PROGRAM MAIN
! ********** Grid Generation
SUBROUTINE GRIGEN(x,y,h,g,c,if,jf)
DIMENSION x(0:if),y(0:jf),h(0:if),g(0:jf),c(if,jf,8)
s=.08 pi=3.14159
DO i=0,if
x1=i/FLOAT(if) x(i)=x1-s*SIN(2.*pi*x1)
ENDDO
DO j=0,jf
y1=j/FLOAT(jf) y(j)=y1-s*SIN(2.*pi*y1)
ENDDO
FORALL(i=0:if-1)h(i)=x(i+1)-x(i)FORALL(j=0:jf-1)g(j)=y(j+1)-y(j)! Setup coefs of first- and second-derivatives
DO i=1,if-1
hm=h(i-1) h0=h(i)DO j=1,jf-1 gm=g(j-1) g0=g(j)c(i,j,1) = hm/h0/(hm+h0) c(i,j,2) = h0/hm/(hm+h0)c(i,j,3) = gm/g0/(gm+g0) c(i,j,4) = g0/gm/(gm+g0)c(i,j,5) = 2./h0/(hm+h0) c(i,j,6) = 2./hm/(hm+h0)c(i,j,7) = 2./g0/(gm+g0) c(i,j,8) = 2./gm/(gm+g0)
ENDDO ENDDO
END SUBROUTINE GRIGEN!1st-deriv coefs!2nd-deriv coefs

6

非圧縮性流れの解法 | !

;法! ********** Calculate vorticity transport eqn by SOR
SUBROUTINE CALOMG(h,g,c,o,p,u,v,if,jf)
DIMENSION h(0:if),g(0:jf),c(if,jf,8),o(0:if,0:jf),p(0:if,0:jf),u(0:if,0:jf),v(0:if,0:jf), &o0(0:20,0:20),c1(0:20,0:20,0:4),rhs(0:20,0:20)
COMMON Re,na,resomg,respsi
alpha=.2 IF(na>10)alpha=.4! Steup coefficients and rhs of vorticity transport eqn
FORALL(i=0:if,j=0:jf,k=1:4)c1(i,j,k)=0.
FORALL(i=0:if,j=0:jf)c1(i,j,0)=1.
FORALL(i=0:if,j=0:jf)o0(i,j)=o(i,j)
DO i=1,if-1 DO j=1,jf-1
up=(u(i,j)+ABS(u(i,j)))/2./h(i) um=(u(i,j)-ABS(u(i,j)))/2./h(i-1)
vp=(v(i,j)+ABS(v(i,j)))/2./g(j) vm=(v(i,j)-ABS(v(i,j)))/2./g(j-1)
c1(i,j,1) = um-c(i,j,5)/Re c1(i,j,2) =-up-c(i,j,6)/Re!lhs coefs of omegac1(i,j,3) = vm-c(i,j,7)/Re c1(i,j,4) =-vp-c(i,j,8)/Re
c1(i,j,0) = -c1(i,j,1)-c1(i,j,2)-c1(i,j,3)-c1(i,j,4)rhs(i,j) = um*(o(i+1,j)-o(i,j))-c(i,j,1)*(u(i+1,j)*o(i+1,j)-u(i,j)*o(i,j)) &+up*(o(i,j)-o(i-1,j))-c(i,j,2)*(u(i,j)*o(i,j)-u(i-1,j)*o(i-1,j)) &+vm*(o(i,j+1)-o(i,j))-c(i,j,3)*(v(i,j+1)*o(i,j+1)-v(i,j)*o(i,j)) &+vp*(o(i,j)-o(i,j-1))-c(i,j,4)*(v(i,j)*o(i,j)-v(i,j-1)*o(i,j-1))
ENDDO ENDDO! Consider Woods conds
dx=h(0) dy=h(0)FORALL(j=1:jf-1)
c1( 0,j,1)=.5 rhs( 0,j)=-3.*p(1,j)/dx/dx!on left side wallc1(if,j,2)=.5 rhs(if,j)=-3.*p(if-1,j)/dx/dx!on right side wall
ENDFORALLFORALL(i=1:if-1)
c1(i, 0,3)=.5 rhs(i, 0)=-3.*p(i,1)/dy/dy!on bottom wallc1(i,jf,4)=.5 rhs(i,jf)=-3.*p(i,jf-1)/dy/dy-3.*u(i,jf)/dy !on top wall
ENDFORALL! Solve vorticity transport eqn by SOR
n=0 20 n=n+1 reso=0.
ie=if IF(MOD(n,2)==0)ie=ie-2
DO ii=0,ie i=ii IF(MOD(n,2)==0)i=if-1-ii
im1=MAX0(i-1,0) ip1=MIN0(i+1,if)
je=jf IF(MOD(n,2)==0)je=je-2
DO jj=0,je j=jj IF(MOD(n,2)==0)j=jf-1-jj
jm1=MAX0(j-1,0) jp1=MIN0(j+1,jf)
res=c1(i,j,1)*o(ip1,j)+c1(i,j,2)*o(im1,j) &+c1(i,j,3)*o(i,jp1)+c1(i,j,4)*o(i,jm1) &+c1(i,j,0)*o(i,j)-rhs(i,j)
cor=-1.2*res/c1(i,j,0) o(i,j)=o(i,j)+cor
reso=AMAX1(reso,ABS(res))
ENDDO ENDDO
IF(reso>1.E-3.AND.n<4)
GOTO 20
FORALL(i=0:if,j=0:jf)o(i,j)=o0(i,j)+alpha*(o(i,j)-o0(i,j))
ENDSUBROUTINE CALOMG! ********** Solve boundary value problem of stream function eqn by SOR
SUBROUTINE CALPSI(c,o,p,u,v,if,jf)
DIMENSION c(if,jf,8),o(0:if,0:jf),p(0:if,0:jf),u(0:if,0:jf),v(0:if,0:jf)
COMMON Re,na,reso,respn=0 30 n=n+1 resp=0.
ie=if-1 IF(MOD(n,2)==0)ie=ie-2
DO ii=1,ie i=ii IF(MOD(n,2)==0)i=if-1-ii

7

対流差分法
je=jf-1 IF(MOD(n,2)==0)je=je-2
DO jj=1,je j=jj IF(MOD(n,2)==0)j=jf-1-jj
res=c(i,j,5)*(p(i+1,j)-p(i,j))-c(i,j,6)*(p(i,j)-p(i-1,j)) &+c(i,j,7)*(p(i,j+1)-p(i,j))-c(i,j,8)*(p(i,j)-p(i,j-1))+o(i,j)
cor=1.5*res/(c(i,j,5)+c(i,j,6)+c(i,j,7)+c(i,j,8))p(i,j)=p(i,j)+cor
resp=AMAX1(resp,ABS(res))
ENDDO ENDDO
IF(resp>1.E-5.AND.n<4)
GOTO 30! Compute velocities
FORALL(i=1:if-1,j=1:jf-1)
u(i,j)= c(i,j,3)*(p(i,j+1)-p(i,j))+c(i,j,4)*(p(i,j)-p(i,j-1))v(i,j)=-(c(i,j,1)*(p(i+1,j)-p(i,j))+c(i,j,2)*(p(i,j)-p(i-1,j)))ENDFORALLEND SUBROUTINE CALPSI

このプログラムについて少し 説明する．記号の意味は分かると思うが念のために示せば ，x = x, y = y ,omg = ! , psi =  , u = u, v = v また c は 1 階と 2 階差分の係数である．Re =j j

R は Reynolds 数，na = napprox

j jは反復数である．resomg = max r! は式 (12.15) と Woods の条件の最大残差値，また respsi = max r

は式 (12.17) の最大残差値である．このプログラムではサブルーチン GRIGEN で格子形成と係数 c の設定を行った後に，サブルーチン CALOMG の ! の計算とサブルーチン CALPSI の  の計算を反復し ，解の収束(maxjr! j=100: + maxjr j  0:0005) または所定の反復数 (na=naf) で計算を終了し ，!  u v の計算結果を打ち出している．サブルーチン CALOMG では式 (12.15) に Woods の条件を考慮したものを SOR 法で解いている．またサブルーチン CALPSI では式 (12.17) に Dirichlet 条件  = 0 を考慮したものを SOR 法で解き，また流速 u v を式 (12.19) から計算している．なお SOR 法の計算では掃引方向を毎回反転し，境界のところで点過緩和を続けて行わないようにしている．

計算に用いた格子と結果を図 12.1 に示す．計算格子の間隔は壁面のところで 0.025，キャビティ中心で0.075，また計算の結果は，解が十分収束した時点の渦度分布と流線である．最後に解の収束状況について述べる．このプログラムでは SOR 法の計算における続けての最大掃引数を n = 4 としているが，その根拠は下表による．なお掃引数 n = 2，ダンピング係数  = 0:35 では計算開始時に少し不安定になる．

n
2

0.30
0.35
0.40
4 0.300.35
0.40
10 0.15
0.20

12.3

napprox

139
126

発散

128
115

発散

183

CPU-time(sec)
0.66
0.44
0.66
0.39
1.10

発散

対流差分法

実質微分の演算子 d=dt

 @=@t + u r は流跡線に沿う内微分 (interior derivative) を表している．したがって式 (12.4) を図 12.2 に示す x t 空間内の１つの流跡線に沿って積分すれば次式が得られる．

!ijn+1 = ! +Z tn+1tn

r2 ! dt(12.22)非圧縮性流れの解法 | !

8;法図 12.1: 正方形キャビティ流れの計算格子，渦度分布と流線

r2 = @ 2=@x2 + @ 2=@y2 ，また上添字 * はこの流跡線と t = tn 面の交点すなわち線積分の起点を表す．
式 (12.22) は厳密な式で，流体粒子が時間間隔 t = tn+1 ; tn に１つの流跡線に沿って点 (tn x  y ) から

ただし

n+1 に変化するという物理現点 (tn+1 xi  yj ) まで移動する間に，流体粒子の渦度が拡散によって ! から !ij
象を直接的に記述している．式 (12.22) は対流差分 (convective-dierence) 法によって数値的に計算することができる．

!ijn+1 = ! +  r2 !t(12.23)

上添字 * の付いている起点の量はその量の時間 t = tn における格子点の値から補間によって求めることができ，またその際に必要な起点の位置 x は次式で与えられる．

Z tn+1x = x ij ; n u dt
また x(12.24)t;xij =  も対流差分法によって計算することができる．

 = ; ut(12.25)流跡線

H!ijn+1

t6
y

;
; -x

;;
;;
;;
;;
;;! ;; ;x;;;H t 面;; x ;;ijn

図 12.2: x t 空間内の積分路

9

対流差分法

対流差分の式 (12.23)，(12.25) はここでは予測子修正子法で計算することにする．これらの式中の−の付いている量は次のように置かれる．

n+1 )=2
修正子： f = (f + fijn
予測子： f = fij

この対流差分法の基本的計算手順を次に示す．まず予測子の計算では，式 (12.25) から起点の位置  (p) ，下n+1(p) を計算する．次に修正子の計算

記の 1 次補間式を用い起点の渦度 ! ，式 (12.23) から渦度の予測値 !ij

では，1 次補間で起点の流速 u ，式 (12.25) から起点の位置  (c) ，2 次補間で起点の渦度 ! と ( 2 ! ) ，式rn+1(c) を計算する．更に修正子の計算を繰り返すというものである．なお下記(12.23) から渦度の修正値 !ijのプログラムでは予測値はすでに求めた結果を基に時間に関し外挿し求めている．
補間公式については補間と関数近似の章で詳し く述べたので，ここではそのいくつかの式を挙げるにとどめる．まず等間隔長方形格子の場合には，補間公式は一般に次のようになる．

u =

m2
X

Ck ( )k=m1

m2
X

Cl ( ) ui+kj+l(12.26)l=m1

ただし Ck ( ) Cl ( ) は補間係数， = =x==y で，この式は一般に２次元の m2 ; m1 次の補間式である．1 次上流補間の式は

1
u = xf; ; ui;1 +(x ;j j)ui +  +ui+1 g= ; ; ui;1 +(1 ;j j)ui + + ui+1
= ui + ui1=2( =(12.27))

ただし   = (j j)=2 j j =  +; ; x = xi+1;xi  = =x  =  =x j j = +; ; ui+1=2 =ui+1 ;ui である．また補間係数は C;1 = ; ;  C0 = 1;j j C1 = + となる．1 次上流補間は常に安定で，1 次上流差分に比べ精度がよい．
次に 2 次補間の式はnu = x1 2 ; 12  (x ;  )ui;1 +(x +  )(x ;  )ui + 21 (x +  )ui+11
1(1 ; )ui;1 +(1 ; 2 )ui + (1+ )ui+1
2
2
1= ui + ui1=2 ;(1 ;j j)(1 ; r  )ui1=2
2

o

=;

ただし r  = ui1=2 =ui1 である．また補間係数は C;1 =(12.28)( =); (1 ; )=2 C0 = 1 ; 2  C1 = (1+ )=2である．この 2 次補間は，精度は十分であるが周知のように対流差分法を不安定にするおそれがあり，何らかの対策をこうじることが必要である．上流側の格子点を用いればある程度安定化することができる．この 2 次上流補間の式は，

I=(;1 ( < 0)
1 (  0)(12.29a)で定義される I を用い，上記の 2 次補間式でi ! i + I  !  ; Ix

! ;Iのような置換えを行えば得られる．

(12.29b)非圧縮性流れの解法 | !

10

;法

TVD の考えに基づく補間式を用いればより完全に安定化できる．2 次の TVD 補間式はn ou = ui + ui 1=2 + (jj; 1) 1 ; minmod r 2+1  2 ui 1=2 ( =  )ただし minmod は次式で定義される制限関数である．

(12.30)minmod(x y) = sign(x) max 0 minfjxj sign(x)yg]
この補間式は，2 次補間式 (12.28) を 1 次上流補間とそれを 2 次にする補正項の和で表し，補正項の大きさを u が区間;1 <  < 1 で極値を取らないように minmod 関数を用い制限したものである．

不等間隔長方形格子の場合には上記の補間公式は以下のようになる．その 1 次上流補間の式は形式的に式 (12.27) と全く同じである．

u = ;; ui;1 +(1 ;jj)ui + + ui+1 = ui + ui 1=2
ただし ( =  )(12.31)=  =xi 1=2  xi+1=2 = xi+1 ; xi である．2 次補間の式は次のように書くことができる．

 (xi+1=2 ;  )(x 1=2 +  )(xi+1=2 ;  )ui;1 + i;x
ui
(
x
+
x
)
i;1=2
i;1=2
i+1=2
i;1=2 xi+1=2(xi;1=2 +  )
+ (x
ui+1
(12.32)i;1=2 + xi+1=2 )xi+1=2
1 ; r u= ui + ui 1=2 ; (1 ;jj) 1+( =  )m i 1=2
ただし m = x 1=2 =x1=2 は格子間隔の比，r = (u 1=2 =x 1=2 )=(u1=2 =x1=2 ) は勾配の比である．また 2 次上流補間の式は式 (12.29a) の I を用い，上の 2 次補間式でu = ; x

i ! i + I  !  ; Ixi 1=2   !  ; I(12.33)のように置換えたものになる．また制限関数を用いる 2 次補間式は次のようになる．

n

 



o

u = ui + ui 1=2 + (1 ;jj) 1 ; minmod r1++mm  2 ui 1=2( =  )(12.34)! の値を 1 次補間で求め渦度輸送方程式を解く方法は Leith 法 (1965) として知られている4 ．この数値
スキームは安定で，その精度は 1 次上流差分よりも良く時空間で 2 次と言われている．もう少し詳し く言えば，1 次補間と 1 次上流差分の誤差は  = 0 の近傍では同程度で jj とともに増加するが，jj が更に増加すると 1 次上流差分では普通 増加の一途をたどるのに対し，1 次補間では途中で減少に転じ jj = 1 で再び誤差 0 になる．また Crowley(1967) は !  の値を 2 次補間で求め精度改善を図っている5 ．粘性の効果が弱く流速がほぼ一様の場合には，このように !  の値を精度良く求めることによって流れを精度良く解析することができる．しかしながら，流れの状態によらず精度を向上させるためには，補間公式のみでは不十分

で，上記の対流差分法のように起点の位置 x  を正確に求め直し ，修正子まで計算することが必要である．
修正子まで計算すれば 1 次補間でも精度をかなり改善できる．更に安定性の良い 2 次補間を用いて修正子を計算すれば精度を一層改善できる．

! ;  を用いる解法は，方程式の数が少なく，かつては良く用いられた．しかしながら 渦度 ! は流速 uを微分した量で u に比べ激しく変動し扱いにくい．乱流の計算に適さず，また 3 次元に拡張できない．
4
5

Leith, C.E., Numerical simulation of the earth's atmosphere, Methods in computational Physics, Vol.4, 1{28, 1965,
Crowley, W.P., Second-order numerical advection, J. of Comput. Physics, Vol.1(1967), 471{484.

11流れ関数の境界値12.4

流れ関数の境界値

本節には流れ関数の境界値の設定について述べる．始めに流路内に静止物体の置かれている２次元流れ;問題を考える．下壁面の流れ関数を 0 上壁面の流れ関数を 1 とすれば， 1
0 はこの流路の流量である．
ただし流路の高さは単位長さとする．静止物体表面の流れ関数 b の瞬間値は一定で 1b は物体の上側を通る流量， b;; 0 は下側の流量になる．この問題で，流路と物体の形状， 0  1 の値，レ イノルズ数がb の値は決まるはずである．しかしながら流れ関数方程式 (12.7) の境界

与えられれば，物理現象としては

b をどのように与えても成立する．平行平板間流路の中央に長方形柱が対称に置かれている問題

値問題は
を!

; 法で解いた論文はかなりあるが，その多くは長方形柱の上下を通る流量を同じとしている．物体後方に渦列の発生するときには当然 物体の上下を通る流量は時間とともに変動する．一様流れの中に円柱や角柱の置かれている問題を !

;法で解いた論文もかなりの数に上るが，そのほとんど すべては遠方の流

1 (x y) と共に b の値を固定しており，渦列を伴う流れを正しく計算しているとは言えない．本

れ関数

節ではまず

12.4.1b の値を決める条件を提示し，次いでこの条件を満足する内部流れと外部流れの計算法を示す．

圧力一価の条件

Navier-Stokes 方程式 (12.3) を２次元流れ場内に取られた一つの閉曲線 ; に沿って積分すれば次式が得られる．図 12.3 参照．

I

dP = Pf ; Ps = P = ;I=;(ut ; v! + !y )dx ;II(ut ; u  ! +  r ! )  dss(vt + u ! ; !x )dy(12.35)

ただし P = p + u 2=2 は全圧，下添字 s f は線積分の始点と終点で同一の点，! =



dx dy





0

0



! ，dss =0 6 ．流れ場内の圧力は一価であるからこの式の左辺は 0 となる．したがって，物理的に許される流れ場においては，次の条件が満足されなければならない．

II(ut ; v! + !y )dx + (vt + u ! ; !x )dy = 0(12.36)

ここではこの条件を圧力一価の条件 (one-valued
 pressure condition) と呼ぶことにする．
閉曲線上の外向き単位法線ベクトルを n =Za  dss =Z(ax dx + ay dy ) =

Z

nx ny

0 で定義すれば，流面内の任意ベクトル a に対し

Z(;ny ax + nx ay )ds =  n  a ds

図 12.3: 圧力一価の条件
6 ベクトル三重積より

r = ;r(ru) = ;r!
2u

非圧縮性流れの解法 | !

12

;

法

となる．今ベクトル a に NS 方程式 (12.3) を代入し ，ガウスの定理を適用すれば次式が得られる7 ．

Ps ; Pf =

ZZ 


d!
;
 r2 ! dxdydt

右辺の被積分関数は流れ場内では式 (12.4) から 0 になり，閉曲線の内部がすべて流体で満たされているときには積分値も 0 になる．したがってこのときには Pf = Ps となり式 (12.36) は満足される．しかしながら，閉曲線の内部に物体があるときには，物体表面上の流れ関数の値が適正に与えられなければ，式 (12.36) は一般に満足されないことになる．図 12.3 に示す流れの中を左回りに一巡する曲線と物体の表面を右回りに一巡する曲線を繋いだ ;1 に対しては式 (12.36) が成立するから，この式の流れの中を一巡する積分値と物体の表面を一巡する積分値は大きさ等しく符号反対になる．これより，圧力一価の条件は物体を一巡する一つの閉曲線に対して満足させれば十分であることが分かる．
物体に作用する流体力は次式で表される．

F=

I(;n p +  n ru )ds(12.37)

ただし積分路は物体表面を一巡するもので，n は物体表面の外向き単位法線ベクトルである．この式の右辺
第 1 項は物体表面の法線方向に作用する圧力，第 2 項は接線方向に作用する粘性力である．ここでは物体に作用する揚力と抗力は F の y または x 方向成分とする．なお諸量は無次元化されているので，これらの揚抗力は揚抗力係数と見ることもできよう．
揚抗力はこれらの式から直接求めることもできるが，流路内の一つの閉曲線に沿う積分から求めること; x 0  r で定義されを導入し ，対流項 u ru を保存形 (発散形)ruu

も可能でその方が計算し易いともいえる．NS 方程式 (12.1) は，u =





y

= 0 0
で表し，また粘性項を上記のように渦度 ! を用いて表せば，次のように書き換えることができる．
る紙面に垂直方向の流れ関数のベクトルr t + ru u + rp +  r! = 0
この式を流体で満たされた閉領域にわたって積分し ，ガウスの定理を適用すれば次式が得られる．

I

(n 

t+

n u u + n p +  n ! )ds = 0(12.38)

今閉領域を物体回りの環状領域に取れば，上式から物体に作用する流体力の式は次のようになる．

I

F =;

(n 

t+

n u u + n p +  n ! )ds(12.39)

この式の右辺は流れの中に任意に取られた物体を一巡する閉曲線に沿う積分で，n は閉曲線上の外向き単位
法線ベクトルである．式 (12.39) の意味すなわち揚抗力の具体的計算の仕方については後述のプログラムを参照されたい．
7ガウスの定理は任意のスカラー量 ，ベクトル量 a に対し次のように書くことができる．

ZZ
ZZ

ZZ

r dxdy =Ira dxdy =Ira dxdy =n  dsIn a ds
n a ds流れ関数の境界値

12.4.213

流路内にある物体まわりの流れ

図 12.4 に示す流路内に一つの物体の置かれている問題をとり上げる．物体上の流れ関数の値を b とする．その予測値は他の変数と同様に時間に関し外挿することによって，また修正値は圧力一価の条件が満足されるように反復計算によって求められる．以下には修正値の求め方について述べる．式 (12.35) を tn+1=2
に関する差分式に書き換えれば次式が得られる．

I
I;P n+1=2 = (ut ; v! + !y )n+1=2 dx + (vt + u ! ; !x )n+1=2 dy
n+1=2

ただし un+1=2 = (un + un+1 )=2 ut(12.40)= (un+1 ; un )=t である．式 (12.40) の右辺に既に求められている!n  u n の値と !n+1  u n+1 の予測値を入れ，物体を左回りに一巡する閉曲線 ; に沿って積分し P n+1=2を求めれば，その値は一般に 0 にならず圧力残差が生じる．物体表面の流れ関数 b の値はこの残差値が 0になるように修正される．
ところで式 (12.40) の被積分関数の各項は，NS 方程式の非定常項，対流項，粘性項に由来するが， b の値を変えたときに直接影響を受けるのは非定常項である．敷衍すれば， b の値を増減すれば，物体回りの循環が増減するので，循環流が加減速されることになり，これが圧力残差を減増することになる．一方，渦度 ! は物体表面近くで発生し流れに乗って対流し拡散によって減衰していくもので， b の渦度に及ぼす影
響は小さく，また第 2 項の;u ! は流れに垂直に作用する力で，その変化は閉曲線 ; 方向の力の線積分にあまり影響しない．今第 N +1 近似計算の式 (12.40) から第 N 近似計算の式 (12.40) を引き算した式を作り，
第 N +1 近似計算の圧力残差 P n+1=2(N +1) = 0 と置き，また b の値を変えたときに非定常項の変化に較
べ他の 2 項の変化は小さいものとして無視すれば次式が得られる．

P

n+1=2(N )

1

 t

nI

n+1(N +1)

(u dx + v dy )Io; (u dx + v dy)n+1(N )(12.41)

この式の積分は ; 上の循環で，この式は循環流の加減速によって圧力残差を 0 にできることを示している．
次に b と ; 上の循環の関係について考える．上記問題と同一流路内の流れ関数 ~ の境界値問題

r2 ~ = 0(~=(12.42)1
0(物体表面)(流路壁面，上流下流境界)の解は，物体まわりの非粘性ポテンシャル循環流れを表している．この流れの ; 上の循環は

Ib = ;

I @~
@n

ds

(12.43)図 12.4: 流路内にある物体まわりの流れ非圧縮性流れの解法 | !

14
ただし @=@n = n;

法

r，n は閉曲線 ; 上の外向き単位法線ベクトルである．この循環はもとの粘性流れの物体

表面の流れ関数の値を単位量増やしたときの ; 上の循環の増分を近似するものである．これより式 (12.41)(N +1) を乗じたものになり，この補正値は次式の右辺の循環の変動は Ib に物体表面の流れ関数の補正値  bから求めることができる．

Ib  b(N +1) = P n+1=2(N ) t

(12.44)

このポテンシャル流れの重畳は，渦度が b の値の変化の影響を受けないという上記の仮定と両立し，おおむね妥当なものといえる．しかしながら壁面近傍の流れの減速を無視していることは，循環流の増分 Ib(N +1) を過大に見積もることになる．それゆえ 1 よりも小さい係数を過小に評価し， bn+1(N +1)

0 を用い

= n+1(N ) + 0  (N +1)(12.45)のように置けば収束性を改善できる．なお実際の計算では，あらかじめ式 (12.42)，(12.43) から Ib を求め，各時間ステップの計算では予測値は時間に関する外挿によって，また修正値は式 (12.44)，(12.45) から求めることになる．

12.4.3一様流中にある物体まわりの流れ

次に図 12.5 に示す一様流中に複数の物体の置かれている問題を採り上げる．物体の数を M 個としその番めの物体表面上の流れ関数の値を  とする．その予測値は時間に関し外挿することによって，また修正
値は圧力一価の条件が満足されるように反復計算によって求められる．ここでは上記の修正値の求め方を一様流中にある複数物体の場合に拡張する．各物体を左回りに一巡する閉曲線 ; ( = 1 2 : : :  M ) に対し，
n+1=2(N )式 (12.40) の右辺に既に求められている ! n  u n と ! n+1(N )  u n+1(N ) の値を入れ，圧力残差値 Pを求める．これらの残差値はすべて 0 になるべきもので，それは物体表面の  ( = 1 2 : : :  M ) の値を適正に選ぶことによって実現できる．圧力一価の条件を満足させるための上記の循環を修正する式 (12.41)

は，この場合にも各 ; に対し有効である．

次に  と ; 上の循環の関係について考える．M 個の物体が同様に配置されている流れ場において，次の流れ関数 ~ の境界値問題を解く．

r2 ~ = 0(
~ =( = 1 2 : : :  M )1( 番めの物体表面)

0( 番め以外の物体表面，外側境界)

図 12.5: 一様流中にある複数物体まわりの流れ(12.46)

15

外部流れの外側境界
その解は番めの物体まわりの非粘性ポテンシャル循環流れを表している．この流れにおける各閉曲線 ;

上の循環は次のようになる．

I = ;

I @~



 @n

ds

( = 1 2 : : :  M  = 1 2 : : :  M )(12.47)

式 (12.41) の右辺の循環の変動は，I の 1 次結合で表され，その係数は物体表面の流れ関数の補正値  になる．したがって次式が成立する．
M
X
=1

I  (N +1) = Pn+1=2(N ) t( = 1 2 : : :  M )(N +1) はこの連立 1 次方程式を解くことによって求められる．また修正値は適当な係数

補正値  (12.48)0 を考慮した次式から求められる．
n+1(N +1)=

n+1(N )+
0  (N +1)(12.49)

なお実際の計算では，あらかじめ式 (12.46) の境界値問題を解いて式 (12.47) から I (   = 1 2 : : :  M )を求め，各時間ステップの計算では予測値は時間に関する外挿によって，また修正値は式 (12.48)，(12.49)から求めることになる．

12.5

外部流れの外側境界

一様流中に置かれた翼型など の物体周りの流れを解く場合には，外側境界は十分に遠方に取らなければならない．しかしながらかつてはメモリーや計算量との兼ね合いであまり遠方に取ることができなかった．
以下に述べる方法は外側境界の外部領域に非粘性ポテンシャル流れを補い外側境界を物体に近付けようというものである．粘性効果はもちろん流れ場全体に及ぶが，NS 方程式の粘性項は境界層と後流の外側では無視でき，この領域では NS 方程式をポテンシャル流れの方程式で置き換えても何ら差し支えない．ここに提案する方法は，外側境界の内部領域は粘性流れとして !

;法で，外部領域は非粘性ポテンシャル流れ

として解析的手法により解くものである．外側境界の位置は，実際にはこれら２つの領域の解の接続を完璧に実行することが難しく，物体からある程度離さざるをえない．プログラムは多少複雑になるが，少ない計算時間で精確な解が得られる．

図 12.6: 外部流れの外側境界非圧縮性流れの解法 | !

16

;法

図 12.6 に示すように，物体近傍と後流の外側に外側境界 ;0 を取り，外部流れ場を内部領域と外部領域に分ける．外部領域のポテンシャル流れの支配方程式は極座標 r  を用い次のように書くことができる．

r2  = rr + 1r r + r12  = 0

この方程式の十分遠方の一様流の条件 r(12.50)! 1 : u ! (U 0) を満足する解  =  (x y) は変数分離法で次のように求められる．

 = Uy + c0 0 +

L
X
l=1

cl l

(12.51a)ただし0 = log(r1 =r)(12.51b)2m;1 = r;m cos m 2m = r;m sin m
(m = 1 2 : : : )この式の右辺第 1 項は一様流れ，第 2 項は外部流れの循環を表している．なお式 (12.50) の解には定数項も含まれるが，ここでは遠方の流れ関数を固定し物体表面の流れ関数を圧力一価の条件を満足するように修正する立場を取っているので，定数項は除かれている．物体が一つの場合には，物体表面の流れ関数の値を固定し圧力一価の条件を満足するようにこの定数の値を決定し一様流の流れ関数の値を修正することもできる．また第 3 項の 1 = (1=r ) cos  の項は物体の前後で循環流を打ち消す作用，また 2 = (1=r ) sin の項は物体近傍の流線を膨らます作用をする．第 3 項は循環項の後始末と物体や境界層，後流による排除効果の項である．固有関数 l の値は次式で計算するのが便利である．

1 = x=r2  2 = y=r2 2m;1 = (x2m;3 ; y2m;2 )=r22m = (y2m;3 + x2m;2 )=r2

(m = 2 3 : : : )(12.52)式 (12.51a) の未定係数 cl の値は，境界 ;0 上の内部流れとの接続条件

=
@ = @@n @n(12.53a)(12.53b)を満足するように決定される．ラプラス方程式で支配される場では滑らかな境界上で関数自身とその法線微分値が連続であれば関数はこの境界を越えて解析接続8 される．これらの式に式 (12.51a) を代入すれば ;0上の k 番めの格子点に対して次式が得られる．

k = Uyk +

L
X
l=0

cl lk(12.54a)Xf^k = cl a^lk
L

(12.54b)

l=0

ただし

(

f^ = vU ;u(a^l = ;@l =@x;@l =@y(12.54c)

8 解析接続とは関数自身とすべての階数の微分値が Lipschitz 連続になることをいう．これを滑らか (smooth) ともいう．十分滑らか (suciently smooth) とは関数自身と必要なある階数までの微分値が連続なことをいう．

外部流れの外側境界

17

上側の式は前方境界 (x = const.)，下側の式は上下境界 (y = const.) に適用される．a^l の値は次式で計算するのが便利である．

(

(

(

2
2
2
4
4
a^0 = x=r2  a^1 = (x ; y4 )=r  a^2 = 2xy=r

y=r
2xy=r
(y2 ; x2 )=r4
(12.55)a^2m;1 = (m ;m1) r2 (xa2m;3 ; ya2m;2 )
(m = 2 3 : : : )a^2m = (m ;m1) r2 (ya2m;3 + xa2m;2 )内外領域の解の接続は，式 (12.54a) を境界 ;0 上のすべての格子点で満足させることは容易であるが，同時に式 (12.54b) も満足させることは困難である．ここでは反復計算の過程で，まず式 (12.54b) により内部
流れの u の値をもとに外部流れの未定係数 cl を求め，次に式 (12.54a) により内部流れの計算に必要な kの値を外部流れの  に等しいものとして与えることにする．次に cl の求め方について説明する．固有関数の数 L は ;0 上の格子点の数 K よりもかなり少なく，f^ の近似に必要最低なものが選ばれる．そのために
cl の値は，式 (12.54b) から直接ではなく，最小２乗法によってこの式が平均的に満足されるように求められることになる．式 (12.54b) の残差の２乗の和

R

K X
L
X
k=1 l=02
cl a^lk ; f^k = minが最小になるように，R を cm (m = 0 1 2 : : :  L) で微分し 0 と置く．すなわち@@cm R(c0  c1  : : :  cL ) = 0(m = 0 1 2 : : :  L)

これより次式が導かれる．

L
X
l=0

aml cl = fm

(m = 0 1 2 : : :  L)

(12.56a)

ただし

aml =
fm =

K
X

a^mk a^lk(12.56b)a^mk f^k(12.56c)

k=1
K
X

k=1

未定係数 cl (l = 0 1 2 : : :  L) はこの連立 1 次方程式を解くことによって求められる．
筆者らは，この外側境界の外側にポテンシャル流れを補う方法が，外側境界をかなり物体に近寄せても流れを正確に計算できるものであることを揚抗力の変化も含めて確かめている．なお cl の計算で，最小２乗法によらず簡単に選点の数 K を減らすことも考えられるが，この場合には隣り合う選点間で  の値が激しく変動する危険性が高い．

18

非圧縮性流れの解法 | !

12.6! ;  法のプログラム

12.6.1

流路内にある角柱まわりの流れ;法

幅 3.0 の平行平板間流路の中央に置かれた一辺の長さ 0.4 の正方形柱まわりの流れを上記の方法で解析する．上流の境界の流れは十分発達した層流で，放物線速度分布を持ちその平均流速 um = 1 とする．Re = 200．
角柱の辺の長さと平均流速を基にしたレ イノルズ数は 80 である．メインプログラムは３つの部分から構成

される．その始めの部分は準備段階で，格子形成，角柱まわりの循環 I の計算，初期定常流れの計算である．中間部分は各時間ステップにおける渦度 ! ，角柱表面の b ，流れ関数  ，渦度の境界値 !b の反復計算と揚抗力の計算である．また終わりの部分は渦度，流れ関数など の予測値の計算とグラフィックスである．

グラフィックスに関する部分はかなり長くなるのでここでは割愛した．次に FORTRAN 95/90 Free source

form で書かれたプログラムを示す．

program MAIN!*************************************************************************! Problem: Unsteady Flow past a Square Cylinder in Parallel Walled Duct
! Numerical Method: omega-psi Eqns, Convective-Difference Scheme!*************************************************************************PARAMETER(if=88,jf=37)
DIMENSION x(if),y(jf),h(if),g(jf),c(if,jf,8),omgm(if,jf),omg(if,jf),omg1(if,jf),e(if,jf), &d(if,jf),psim(if,jf),psi(if,jf),psi1(if,jf),u(if,jf),u1(if,jf),ub(if,jf),us(if,jf),&v(if,jf),v1(if,jf),vb(if,jf),vs(if,jf),p(if,jf)!xy:coordinates, h,g:grid spaces, c:coefficients of first- and second-derivatives!omgm:omega^n-1, omg:omega^n, omg1:omega^n+1, d:nabla^2omega, uv:velocities!u:u^n, u1:u^n+1, ub:baru, us:u^*, p:total pressure
COMMON i1,i2,j1,j2,dt,Re,n,na
TEXT(x,f1,f2,f3)=-.5*x*(1.-x)*f1+(1.-x*x)*f2+.5*x*(1.+x)*f3
DATA i1,i2,j1,j2/16,22,16,22/ dt,Re,nf/.04, 200., 600/
!i1,i2,j1,j2:square cylinder, dt:time increment, Re:Reynolds number, nf:final nOPEN(20,FILE='OUTPUT.dat')! Grid generation
CALL GRIGEN(x,y,h,g,c,if,jf)! Compute circulations
n=-1
FORALL(i= 1:if,j= 1:jf)psi(i,j)=0.
FORALL(i=i1:i2,j=j1:j2)psi(i,j)=1.
CALL CIRC(h,g,c,psi,if,jf,circ1)! Compute initial steady flow
n=0
CALL STEADY(x,y,h,g,c,omg,omg1,d,psi1,u1,v1,if,jf) GOTO 121!********************
101 n=n+1 t=t+dt
IF(MOD(n,20)==1)WRITE(20,61)
na=0 102 na=na+1! Compute vorticity
IF(n<=10.OR.na<=3)CALL CALOMG(h,g,omgm,omg,omg1,e,d,u,u1,ub,us,v,v1,vb,vs,if,jf)! Compute stream function
IF(na==1)THEN dpsib=psib1-psib dP=0.
GOTO 111 ENDIF
CALL CALDP(h,g,c,omg,omg1,u,u1,v,v1,p,if,jf,dP)!correct psi on square cylinder
dpsib=.7*dt*dP/circ1 psib1=psib1+dpsib
FORALL(i=i1:i2,j=j1:j2)psi1(i,j) = psib1
111 CONTINUECALL CALPSI(h,g,c,omg1,psi1,u1,v1,if,jf,respsi)! Compute bound values of omega and diffusion term of omegaプ ログ ラムIF(n<=10.OR.na<=3)CALL OMGB(h,g,omg1,psi1,if,jf)
IF(n<=10.OR.na<=3)CALL LAPOMG(c,omg1,d,if,jf)
IF(na==1)WRITE(20,62)n,tIF(na<=10.OR.MOD(na,10)==0)WRITE(20,63)na,respsi,psib1,dpsib,dP
IF(n==1)
GOTO 121
IF(na<10)
GOTO 102!********************CALL CALDP(h,g,c,omg,omg1,u,u1,v,v1,p,if,jf,dP)
CALL CALLD(h,g,omg,omg1,psi,psi1,u,u1,v,v1,p,if,jf,Cl,Cd)
CALL CPU_TIME(sec)WRITE(20,64)Cl,Cd,sec
IF(MOD(n,10)==0.AND.n>=500.AND.n<=560)CALL OUTPUT(omg1,psi1,u1,v1,if,jf)
121 CONTINUE! Advance time and get predicted values of omg, psi, etc by extrapolation
DO i=2,if DO j=1,jf
omgmm=omgm(i,j) omgm(i,j)=omg(i,j) omg(i,j)=omg1(i,j)
IF(n>3)omg1(i,j)=TEXT(2.,omgmm,omgm(i,j),omg(i,j))
psimm=psim(i,j) psim(i,j)=psi(i,j) psi(i,j)=psi1(i,j)
IF(n>3)psi1(i,j)=TEXT(2.,psimm,psim(i,j),psi(i,j))u(i,j)=u1(i,j) v(i,j)=v1(i,j)
ENDDO ENDDO
IF(n==1.AND.na<40)
GOTO 102
psibmm=psibm psibm=psib psib=psib1 IF(n>3)psib1=TEXT(2.,psibmm,psibm,psib)
IF(n<nf)
GOTO 101
CLOSE(20)STOP
61 FORMAT(/5X 'nt ', 4X 'na respsi
psib
4Pdpsib', 6X &
' dP
Cl
Cd
CPU-time')
62 FORMAT(5X I4,F8.2)
63 FORMAT(20X I2, F8.4, 2X F8.4, 4PF8.3, 4X 0PF8.4)
64 FORMAT(62X 2F8.4, 4X F8.2)
END PROGRAM MAIN
! ********** Generate nonuniform rectangular grid
SUBROUTINE GRIGEN(x,y,h,g,c,if,jf)
DIMENSION x(if),y(jf),h(if),g(jf),c(if,jf,8),iq(if)pi=3.14159
DO i= 1,15 xx=FLOAT(i-16)/15. x(i)=2.+xx*(xx*xx+1.) ENDDOFORALL(i=16:22)x(i)=2.+FLOAT(i-16)/15.
DO i=23,37 xx=FLOAT(i-22)/15. x(i)=2.4+.5*xx*xx+xx ENDDOFORALL(i=38:if)x(i)=3.9+FLOAT(i-37)/7.5
FORALL(j=1:jf)y(j)=FLOAT(j-19)/15.
DO j= 1,15 yy=FLOAT(j-16)/15. y(j)=y(j)-.6*yy*yy*yy-.9*yy*yy ENDDO
DO j=23,37 yy=FLOAT(j-22)/15. y(j)=y(j)-.6*yy*yy*yy+.9*yy*yy ENDDOFORALL(i=1:if-1)h(i)=x(i+1)-x(i) h(if)=h(if-1)
FORALL(j=1:jf-1)g(j)=y(j+1)-y(j)! Setup coefs of first- and second-order derivatives
FORALL(i=1:if,j=1:jf,k=1:8)c(i,j,k)=0.
DO i=2,if-1 DO j=1,jf hm=h(i-1) h0=h(i)c(i,j,1)=hm/h0/(hm+h0) c(i,j,2)=h0/hm/(hm+h0)!1st-deriv coefsc(i,j,5)=2./h0/(hm+h0) c(i,j,6)=2./hm/(hm+h0)!2nd-deriv coefs
ENDDO ENDDO
DO i=1,if DO j=2,jf-1 gm=g(j-1) g0=g(j)c(i,j,3)=gm/g0/(gm+g0) c(i,j,4)=g0/gm/(gm+g0)!1st-deriv coefsc(i,j,7)=2./g0/(gm+g0) c(i,j,8)=2./gm/(gm+g0)!2nd-deriv coefs
ENDDO ENDDO! Output x,y

19

20

非圧縮性流れの解法 | !

;法FORALL(i=1:if)iq(i)=i
WRITE(20,'(//5X A30/)')'***** grid point x *****'
DO m=1,3 ib=30*(m-1)+1 ie=MIN(ib+29,if)
WRITE(20,'(2X 30F8.4)')(x(i),i=ib,ie) WRITE(20,'(30I8)')(iq(i),i=ib,ie)
ENDDO
WRITE(20,'(//5X A30/)')'***** grid point y *****'
DO m=1,2 jb=30*(m-1)+1 je=MIN(jb+29,jf)
WRITE(20,'(2X 30F8.4)')(y(j),j=jb,je) WRITE(20,'(30I8)')(iq(j),j=jb,je)
ENDDO
ENDSUBROUTINE GRIGEN! ********** Compute circulation
SUBROUTINE CIRC(h,g,c,psi,if,jf,circ1)
DIMENSION h(if),g(jf),c(if,jf,8),psi(if,jf),o(if,jf),u(if,jf),v(if,jf)
DATA i5,i6,j5,j6/14,24,14,24/
FORALL(i=1:if,j=1:jf)o(i,j)=0.
CALL CALPSI(h,g,c,o,psi,u,v,if,jf,resp,dpsi)
circ1=0.
DO i=i5,i6-1 circ1 = circ1+( u(i,j5)+u(i+1,j5)-u(i,j6)-u(i+1,j6))/2.*h(i) ENDDO
DO j=j5,j6-1 circ1 = circ1+(-v(i5,j)-v(i5,j+1)+v(i6,j)+v(i6,j+1))/2.*g(j) ENDDO
WRITE(20,'(//5X A30//16X A8, F7.3//)')'***** circulations *****','circ1 =',circ1
ENDSUBROUTINE CIRC
! ********** Compute initial steady flow using convect-differnce method
SUBROUTINE STEADY(x,y,h,g,c,omg,omg1,d,psi,u,v,if,jf)
DIMENSION x(if),y(jf),h(if),g(jf),c(if,jf,8),omg(if,jf),omg1(if,jf),d(if,jf), &psi(if,jf),u(if,jf),v(if,jf)
COMMON i1,i2,j1,j2,dt,Re,n,na! Initial values
um=1.
DO i=1,if cycle_1: DO j=1,jf
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_1
omg(i,j) = um*2*y(j)/1.5 omg1(i,j) = um*2*y(j)/1.5
psi(i,j) = um*(1.5*y(j)-y(j)*y(j)*y(j)/3./1.5)
u(i,j) = um*(1.5*1.5-y(j)*y(j))/1.5 v(i,j) = 0.
ENDDO cycle_1 ENDDO! Compute steady flow
na=0 10 na=na+1 IF(na==1)
GOTO 11
CALL CALOMG0(h,g,c,omg,omg1,d,psi,u,v,if,jf,domg)!calculation of vorticity
omg(i1-3,j1+2)=10.
!trigger for asymmetric flow
11 CALL CALPSI(h,g,c,omg,psi,u,v,if,jf,respsi)!calculation of stream function
IF(na<=10.OR.MOD(na,10)==0)WRITE(20,60)na,domg,respsi
60 FORMAT(5X 'na =', I4, 4X '2Pdomg =' 2PF8.4, 4X 'respsi =', 0PF8.4)
IF(na<200)
GOTO 10ENDSUBROUTINE STEADY! ********** Calculate vorticity transport eqn using 1st-order convective-differences
SUBROUTINE CALOMG0(h,g,c,o,o1,d,p,u,v,if,jf,domg)
DIMENSION h(if),g(jf),c(if,jf,8),o(if,jf),o1(if,jf),d(if,jf),p(if,jf),u(if,jf),v(if,jf)
COMMON i1,i2,j1,j2,dt,Re,n,na! Compute vorticity using 1st-order convective-differences
CALL LAPOMG(c,o,d,if,jf)
CALL UPWIND(h,g,u,v,o,o1,if,jf)
DO i=2,if ip1=MIN(i+1,if) cycle_1: DO j=2,jf-1
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_1
o1(i,j) = o1(i,j)+d(i,j)*dt/Re
ENDDO cycle_1 ENDDO

プ ログ ラム

21

! Compute vorticity on boundaries by Woods conds
CALL OMGB(h,g,o1,p,if,jf)
domg=0.
DO i=2,if DO j=2,jf-1
dog = o1(i,j)-o(i,j) domg = AMAX1(domg,ABS(dog))
o1(i,j) = o(i,j)+.5*dog o(i,j)=o1(i,j)
ENDDO ENDDO
ENDSUBROUTINE CALOMG0! ********** Solve initial value prob of vorticity transport eqn by convective-differnece method
SUBROUTINE CALOMG(h,g,om,o,o1,e,d,u,u1,ub,us,v,v1,vb,vs,if,jf)
DIMENSION h(if),g(jf),om(if,jf),o(if,jf),o1(if,jf),e(if,jf),d(if,jf),u(if,jf), &u1(if,jf),ub(if,jf),us(if,jf),v(if,jf),v1(if,jf),vb(if,jf),vs(if,jf)
COMMON i1,i2,j1,j2,dt,Re,n,na
IF(na>1)
GOTO 10
FORALL(i=i1:i2)d(i,j1) = (o1(i,j1)-om(i,j1))/(2.*dt)*Re
FORALL(i=i1:i2)d(i,j2) = (o1(i,j2)-om(i,j2))/(2.*dt)*Re
FORALL(j=j1:j2)d(i1,j) = (o1(i1,j)-om(i1,j))/(2.*dt)*Re
FORALL(j=j1:j2)d(i2,j) = (o1(i2,j)-om(i2,j))/(2.*dt)*Re
FORALL(i= 1:if)d(i, 1) = (o1(i, 1)-om(i, 1))/(2.*dt)*Re
FORALL(i= 1:if)d(i,jf) = (o1(i,jf)-om(i,jf))/(2.*dt)*Re
FORALL(i=1:if,j=1:jf)e(i,j) = o(i,j)+dt*d(i,j)/2./Re!e=omg^n+dt*(nabla^2omg)^n/2/Re
RETURN!na=0
10 CONTINUECALL LAGINT(h,g,ub,vb,u,us,if,jf)
CALL LAGINT(h,g,ub,vb,v,vs,if,jf)
FORALL(i=2:if,j=2:jf-1)ub(i,j) = (us(i,j)+u1(i,j))/2.
FORALL(i=2:if,j=2:jf-1)vb(i,j) = (vs(i,j)+v1(i,j))/2.
CALL LAGINT(h,g,ub,vb,e,o1,if,jf)
FORALL(i=2:if,j=2:jf-1)o1(i,j) = o1(i,j)+dt*d(i,j)/2./Re!corrector
ENDSUBROUTINE CALOMG! Compute boundary values of vorticity by Woods cond
SUBROUTINE OMGB(h,g,o,p,if,jf)
DIMENSION h(if),g(jf),o(if,jf),p(if,jf)
COMMON i1,i2,j1,j2
OB(h,o1,p0,p1)=-o1/2.-3.*(p1-p0)/h/h
DO j=j1,j2
o(i1,j) = OB(h(i1-1),o(i1-1,j),p(i1,j),p(i1-1,j))!front wallo(i2,j) = OB(h(i2) ,o(i2+1,j),p(i2,j),p(i2+1,j))!rear wall
ENDDO o11=o(i1,j1) o12=o(i1,j2)
DO i=i1,i2
o(i,j1) = OB(g(j1-1),o(i,j1-1),p(i,j1),p(i,j1-1))!side wallo(i,j2) = OB(g(j2) ,o(i,j2+1),p(i,j2),p(i,j2+1))!side wall
ENDDO o(i1,j1)=(o11+o(i1,j1))/2. o(i1,j2)=(o12+o(i1,j2))/2.
DO i=1,ifo(i, 1) = OB(g(1),o(i,2),p(i, 1),p(i,2))!bottom wallo(i,jf) = OB(g(jf-1),o(i,jf-1),p(i,jf),p(i,jf-1))!top wall
ENDDO
ENDSUBROUTINE OMGB!********** Compute Laplacian omega
SUBROUTINE LAPOMG(c,o,d,if,jf)
DIMENSION c(if,jf,8),o(if,jf),d(if,jf)
COMMON i1,i2,j1,j2
DO i=2,if cycle_1: DO j=2,jf-1
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2) CYCLE cycle_1

22

非圧縮性流れの解法 | !

;法d(i,j) = c(i,j,7)*(o(i,j+1)-o(i,j))+c(i,j,8)*(o(i,j-1)-o(i,j))IF(i==if)CYCLE cycle_1
d(i,j) = c(i,j,5)*(o(i+1,j)-o(i,j))+c(i,j,6)*(o(i-1,j)-o(i,j))+d(i,j)
ENDDO cycle_1 ENDDO
ENDSUBROUTINE LAPOMG!********** Compute pressure difference along a closed curve
SUBROUTINE CALDP(h,g,c,o,o1,u,u1,v,v1,p,if,jf,dP)
DIMENSION h(if),g(jf),c(if,jf,8),o(if,jf),o1(if,jf), &u(if,jf),u1(if,jf),v(if,jf),v1(if,jf),p(if,jf)
COMMON i1,i2,j1,j2,dt,Re
i5=i1-2 i6=i2+2 j5=j1-2 j6=j2+2
p(i5,j5)=0. j=j5 10 CONTINUE
DO i=i5,i6 f1=f2
v0=(v(i,j)+v1(i,j))/2. om=(o(i,j-1)+o1(i,j-1))/2.
o0=(o(i,j)+o1(i,j))/2. op=(o(i,j+1)+o1(i,j+1))/2.
f2 = -(u1(i,j)-u(i,j))/dt+v0*o0-(c(i,j,3)*(op-o0)+c(i,j,4)*(o0-om))/Re
IF(i/=i5)p(i,j) = p(i-1,j)+(f1+f2)/2.*h(i-1)
ENDDO
i=i5 IF(j==j6)THEN i=i6 dP=p(i,j) ENDIF
DO j=j5,j6 f1=f2
u0=(u(i,j)+u1(i,j))/2. om=(o(i-1,j)+o1(i-1,j))/2.
o0=(o(i,j)+o1(i,j))/2. op=(o(i+1,j)+o1(i+1,j))/2.
f2 = -(v1(i,j)-v(i,j))/dt-u0*o0+(c(i,j,1)*(op-o0)+c(i,j,2)*(o0-om))/Re
IF(j/=j5)p(i,j) = p(i,j-1)+(f1+f2)/2.*g(j-1)
ENDDO
j=j6 IF(i==i5)
GOTO 10dP = p(i,j)-dP
ENDSUBROUTINE CALDP! ********** Solve boundary value problem of stream function eqn using SOR
SUBROUTINE CALPSI(h,g,c,o,p,u,v,if,jf,resp)!o=omg, p=psi
DIMENSION h(if),g(jf),c(if,jf,8),o(if,jf),p(if,jf),u(if,jf),v(if,jf)
COMMON i1,i2,j1,j2,dt,Re,n,na
ni=0 10 ni=ni+1 resp=0.
DO i=2,if cycle_1: DO j=2,jf-1
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_1
res = c(i,j,7)*(p(i,j+1)-p(i,j))+c(i,j,8)*(p(i,j-1)-p(i,j))+o(i,j)
IF(i/=if)res = c(i,j,5)*(p(i+1,j)-p(i,j))+c(i,j,6)*(p(i-1,j)-p(i,j))+resp(i,j) = p(i,j)+1.5*res/(c(i,j,5)+c(i,j,6)+c(i,j,7)+c(i,j,8))
IF(i<=58)resp = AMAX1(resp,ABS(res))
ENDDO cycle_1 ENDDO
IF(ni<7)
GOTO 10IF(n==-1.AND.ni<100)
GOTO 10IF((n==0.OR.n==1).AND.na==1.AND.ni<100)
GOTO 10DO i=1,if cycle_2: DO j=1,jf! Compute velocities
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_2
IF(j== 1)THEN u(i,j) = (p(i,j+1)-p(i,j))/g(j)ELSEIF(j==jf)THEN u(i,j) = (p(i,j)-p(i,j-1))/g(j-1)
ELSE u(i,j) = c(i,j,3)*(p(i,j+1)-p(i,j))+c(i,j,4)*(p(i,j)-p(i,j-1)) ENDIF
IF(i== 1)THEN v(i,j) =-(p(i+1,j)-p(i,j))/h(i)ELSEIF(i==if)THEN v(i,j) =-(p(i,j)-p(i-1,j))/h(i-1)
ELSE v(i,j) =-c(i,j,1)*(p(i+1,j)-p(i,j))-c(i,j,2)*(p(i,j)-p(i-1,j)) ENDIF
ENDDO cycle_2 ENDDO
ENDSUBROUTINE CALPSI

プ ログ ラム!********** Compute lift and drag coefficients
SUBROUTINE CALLD(h,g,omg,omg1,psi,psi1,u,u1,v,v1,p,if,jf,Cl,Cd)
DIMENSION h(if),g(jf),omg(if,jf),omg1(if,jf),psi(if,jf),psi1(if,jf), &u(if,jf),u1(if,jf),v(if,jf),v1(if,jf),p(if,jf)
COMMON i1,i2,j1,j2,dt,Re
i5=i1-2 i6=i2+2 j5=j1-2 j6=j2+2
j=j5 11 Cl1=0. Cd1=0.
DO i=i5,i6 f1=f2 d1=d2
u0=(u(i,j)+u1(i,j))/2. v0=(v(i,j)+v1(i,j))/2. omg0=(omg(i,j)+omg1(i,j))/2.
f2 = p(i,j)-u0*u0/2.+v0*v0/2.
d2 = (psi1(i,j)-psi(i,j))/dt+v0*u0+omg0/Re
IF(i>i5)Cl1 = Cl1+(f1+f2)*h(i-1)/2.
IF(i>i5)Cd1 = Cd1+(d1+d2)*h(i-1)/2.
ENDDO
i=i6 IF(j==j6)i=i5 Cl2=0. Cd2=0.
DO j=j5,j6 f1=f2 d1=d2
u0=(u(i,j)+u1(i,j))/2. v0=(v(i,j)+v1(i,j))/2. omg0=(omg(i,j)+omg1(i,j))/2.
f2 = (psi1(i,j)-psi(i,j))/dt-u0*v0+omg0/Re
d2 = -p(i,j)-u0*u0/2.+v0*v0/2.
IF(j>j5)Cl2 = Cl2+(f1+f2)*h(i-1)/2.
IF(j>j5)Cd2 = Cd2+(d1+d2)*h(i-1)/2.
ENDDO
IF(i==i6)THEN Cl = Cl1+Cl2 Cd = Cd1+Cd2 j=j6 GOTO 11
ELSE Cl = Cl-Cl1-Cl2 Cd = Cd-Cd1-Cd2 ENDIF
ENDSUBROUTINE CALLD!********** Output computational results
SUBROUTINE OUTPUT(o,p,u,v,if,jf)
DIMENSION o(if,jf),p(if,jf),u(if,jf),v(if,jf),q(if,jf),iq(if)CHARACTER*15 z1
FORALL(i=1:if)iq(i)=i
DO k=1,2
SELECT CASE(k)CASE(1) z1='
vorticity' FORALL(i=1:if,j=1:jf)q(i,j)=o(i,j)/100.
CASE(2) z1='stream function' FORALL(i=1:if,j=1:jf)q(i,j)=p(i,j)
ENDSELECTWRITE(20,'(//1H 5X A6,A15,A6/)')'***** ', z1, ' *****'
DO m=1,3 ib=30*(m-1)+1 ie=MIN(ib+29,if)
DO j=jf,1,-1 WRITE(20,'(I3, 30F9.3)')j,(q(i,j),i=ib,ie) ENDDO
WRITE(20,'(30I9 /)')(iq(i),i=ib,ie)
ENDDO ENDDO
ENDSUBROUTINE OUTPUT!***** 1st-order upwind interpolation
SUBROUTINE UPWIND(h,g,u,v,f,fs,if,jf)
DIMENSION h(if),g(jf),u(if,jf),v(if,jf),f(if,jf),fs(if,jf)
COMMON i1,i2,j1,j2,dt
DO i=2,if ip1=MIN(i+1,if) cycle_1: DO j=2,jf-1
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_1
ap=-(u(i,j)-ABS(u(i,j)))/2.*dt/h(i) am=-(u(i,j)+ABS(u(i,j)))/2.*dt/h(i-1) aa=ap-am
bp=-(v(i,j)-ABS(v(i,j)))/2.*dt/g(j) bm=-(v(i,j)+ABS(v(i,j)))/2.*dt/g(j-1) ba=bp-bm
fs(i,j) = -am* bp*f(i-1,j+1)+(1.-aa)* bp*f(i,j+1)+ap* bp*f(ip1,j+1) &-am*(1.-ba)*f(i-1,j )+(1.-aa)*(1.-ba)*f(i,j )+ap*(1.-ba)*f(ip1,j ) &+am* bm*f(i-1,j-1)-(1.-aa)* bm*f(i,j-1)-ap* bm*f(ip1,j-1)
ENDDO cycle_1 ENDDO
ENDSUBROUTINE UPWIND

23

24

非圧縮性流れの解法 | !

;法!***** 2nd-order Lagrangian interpolation
SUBROUTINE LAGINT(h,g,u,v,f,fs,if,jf)
DIMENSION h(if),g(jf),u(if,jf),v(if,jf),f(if,jf),fs(if,jf)
COMMON i1,i2,j1,j2,dt
DO i=2,if cycle_1: DO j=2,jf-1
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2) CYCLE cycle_1
ii=i xi=-u(i,j)*dt et=-v(i,j)*dt
IF(i>=3.AND.i<=i1+1.OR.i==if)THEN ii=i-1 xi=xi+h(ii) ENDIF
hm=h(ii-1) h0=h(ii) gm=g(j-1) g0=g(j)cxm=-xi*(h0-xi)/hm/(hm+h0) cx0=(hm+xi)*(h0-xi)/hm/h0 cxp=xi*(hm+xi)/h0/(hm+h0)
cym=-et*(g0-et)/gm/(gm+g0) cy0=(gm+et)*(g0-et)/gm/g0 cyp=et*(gm+et)/g0/(gm+g0)
fm = cym*f(ii-1,j-1)+cy0*f(ii-1,j)+cyp*f(ii-1,j+1)
f0 = cym*f(ii ,j-1)+cy0*f(ii ,j)+cyp*f(ii ,j+1)
fp = cym*f(ii+1,j-1)+cy0*f(ii+1,j)+cyp*f(ii+1,j+1)
fs(i,j) = cxm*fm+cx0*f0+cxp*fp
ENDDO cycle_1 ENDDO
ENDSUBROUTINE LAGINT

上記のプログラムを少し詳しく説明する．サブルーチン GRIGEN では格子形成と，１階と 2 階の中心差分の係数の設定が行われる．角柱は xi1xxi2  yj1yyj2 に置かれている．長方形格子の不等間隔の部分は 3 次多項式を基に格子間隔が滑らかに変化するように作られる．h(i) = xi+1=2 ，g(j) = yj +1=2 ．

また配列 c は中心差分の係数で，これらの係数は式 (12.14)，(12.17)，(12.19) など の計算に用いられる．

n = ;1 では，b の修正計算に必要な角柱回りの循環 circ1 = Ib が求められる．サブルーチン circ では
式 (12.42) のポテンシャル流れの境界値問題が解かれ，Ib の値が式 (12.43) から計算される．
n = 0 では，非定常流れの計算を確実に出発させるために初期定常流れが求められる．サブルーチン
STEADY では，初期定常流れの初期値として平行平板間流路の十分発達した流れが与えられ，定常流れが! ! !b !  の反復計算によって求められる．この計算はどこから始めても良いが， の初期値の不連続を除くべく  の計算から始めるのが適当である．! の計算はサブルーチン CALOMG0 を引用し ，対流差分法の式 (12.23) に 1 次上流補間の式 (12.31) を用いたものによって行われる．ここでは精度よりも安定性が重視される． と u v の計算はサブルーチン CALPSI の中で行われる．この反復計算は解が適当に収束した

時点で打ち切られる．なおここではトリガー渦を付加し周期的渦列の発生を促進しているために解は収束しないが，この渦を付加しないときの解の収束状況は次のようになる．

ただし

N

k

! k

k

 k

100
200

0.00130.00050.00280.0016k ! k= max j!ij(N ) ; !ij(N ;1) j k  k= max jij(N ) ; ij(N ;1) j である．この計算はこの程度の反復で十分である．
時間ステップ n

 1 の非定常流れの計算では，各時間ステップごとに反復計算が行われ，その第 N 近似

計算では

!n+1 ! bn+1 ! n+1 ! !bn+1 ! r2 !n+1の順に計算が行われる．! の計算はサブルーチン CALOMG の中で行われる．その予測値は前ステップの最後に時間に関し外挿されたものがそのまま用いられ，また修正値は対流差分法の式 (12.23) を書換えた次式から求められる．
;
n+1

!ij


= ! + 12  r2 !t + 12  (r2 !)nij+1 t

25

プ ログ ラムN = 1 では ! の計算は不要で，次に備え角柱表面の d = (r2 !)n と e = (! +  r2 !t=2)n の計算が行われる．N = 2 3 では上式から ! の修正値が求められる． の値は式 (12.32) の 2 次補間によって求められるが，角柱前方領域では安定化のため式 (12.33) の 2 次上流補間によって求められる．なおレ イノルズ数が大
きく不安定性が問題になるときには 2 次 TVD 補間式 (12.34) が用いられ，更に物体前面の隣接格子点には指数関数を用いた補間式が併用される．修正値の計算は第 3 近似までで十分である．

 の計算は SOR 法によって行われ，その掃引の反復の途中で角柱表面の b の修正も行われる．b の値
は，N = 1 では前ステップの最後に時間外挿されたものがそのまま用いられ，N  2 ではサブルーチン(N +1)CALDP を引用して式 (12.40) から圧力残差 dP= P n+1=2 を求め，式 (12.44) から補正値 dpsib= b
を求め，b の値は式 (12.45) によりダンピング係数 0.7 を考慮して決定される．領域内部の  の計算は定常流れの場合と同じで，サブルーチン CALPSI の中で 2 次中心差分の式 (12.17) が SOR 法で解かれる．このサブルーチンでは続いて流速 u v の値が式 (12.19) から求められる．

サブルーチン OMGB では角柱と流路壁境界の !b の値が Woods の条件 (12.20) から求められる．またサブ

ルーチン LAPOMG では式 (12.23) の拡散項のr2 ! の値が式 (12.14) から求められる．以上要するに各時間ス

テップにおける上記の反復計算は正確には次の手順で行われる．

! n+1(1) ! !bn+1(1) ! r2 !n+1(1)! !n+1(2) ! bn+1(2) ! n+1(2) ! !bn+1(2) ! r2 !n+1(2)! !n+1(3) ! bn+1(3) ! n+1(3) ! !bn+1(3) ! r2 !n+1(3)! bn+1(4) ! n+1(4)
n+1(10)! b! n+1(10)( ! LD )enj j

各時間ステップにおける解の収束状況は  の差分方程式 (12.17) の最大残差値 respsi = max r から判定される．ただしここではこの最大残差値の大きさによって反復を止めるのではなく 10 回ずつ反復している．
各時間ステップの解が収束した時点で，揚抗力 Cl = L，Cd = D が求められる．諸量は無次元化されているのでこれらの揚抗力は揚抗力係数でもある．サブルーチン CALLD では揚抗力が式 (12.37) から求められる．
各時間ステップの最後に次の時間ステップの ! n+1   n+1  bn+1 の予測値が時間に関する 2 次外挿式によって求められる．omgmm = ! n;2 ，omgm = ! n;1 ，omg = ! n ，omg1 = ! n+1 である．

図 12.7 はこの流れの計算結果で，渦度分布と流線のフローパターン，計算格子，角柱に作用する揚抗力と角柱表面の流れ関数の値の時間変化を示す．格子は十分滑らかに変化するように作られている．フローパ

ターンは流れが十分規則的に周期変化するようになった後のもので，1 周期をやや超える時間にわたるものである．渦度分布は強力な周期渦の発生を示している．ここでは;2:2  x  8:5 の範囲を計算し x  4:5の範囲を示した．なお下流側境界の近傍では壁面近傍の流れが剥離し一部逆流が見られる．抗力は時間的にほぼ一定になるが，揚力や角柱表面の流れ関数の値は渦の発生に同期してかなり変動する．
最後にこの計算の収束性について述べる．この計算では上記のように予測値を前時間ステップの終わりに 2 次式で外挿して与えているが，この予測値の精度はかなりのものである．! の収束状況は第 2 近似でk= maxj!ijn+1(N ) ; !ijn+1(N ;1) j < 0:3，第 3 近似で k ! k< 0:025 になり，! の大きさを考えればこれで十分と言える． の差分方程式 (12.17) の最大残差値 respsi = max jrij j < 0:001，この値はこの式の各項が渦度 ! の大きさであることを考えれば非常に小さい．また圧力残差値 jdPj = jP j < 0:0004，角(N )(N ;1)
j < 3  10;6 である．なお下流側境界近傍には多少の問題があ
柱表面の jdpsibj = jb j = jb ; b
り，これらの値はすべて x  4:5 の範囲で計算したものである．
kdomg = !

26

非圧縮性流れの解法 | !

;法図 12.7: 平行平板間流路内の正方形柱まわりの流れ27

プ ログ ラム12.6.2一様流中にある 2 個の角柱まわりの流れ

一様流中に一辺の長さ d = 0:4 の正方形柱が串形に置かれている場合を解析する．一様流の速度 U = 1とする．Re=200．d と U を基にしたレ イノルズ数は 80 である．内部領域と外部領域を区切る外側境界は前方が x =

;2 02，側方が = 1 2 に取られる．始めにこの計算のフローチャートを示す．計算の手順は，:y:

上記の内部流れの場合とだいたい同じであるが，外部領域のポテンシャル流れの計算が加わり少し複雑に^l  a;1 の計算と，各時間ステップごとの反復に cl  cl   f^なっている．すなわち準備段階に l  a

k

lm

k

k

k

の計算が加わることになる．次にこの計算に用いたプログラムを示す．なお CG に関する部分はかなり長くなるのでこのプログラムでは除かれている．

?

n = n +1

START N =0e

Grid Generationn = ;1

N = N +1

Compute I 
Compute l  a^l  a;lm1
n=0

Initial Steady Flow

Compute !

PPPYesPPn=1?
P

No
Compute P  b  k b ke

Compute cl  k cl k  k f^kCompute  k  k uCompute !b?
eCompute r2 !

Predict ! PPPPPYes


n=1
P
PPP & N<40?

PNo
Predict cl  b
Pf ?PP
Yes

P
Pn<n
PNo

STOP  WRITE
PPP
Pn=1?
e Yes
P
PNo
6PPPYes - e
PPN<10?
PNo
6Compute P L D

 WRITE

図 12.8: 外部流れのフローチャート

program MAIN!***********************************************************************! Problem: Unsteady Flow past Tandem Square Cylinders in Uniform Flow
! Numerical Method: omega-psi Eqns, Convective-Difference Scheme!***********************************************************************

28

非圧縮性流れの解法 | !

;法PARAMETER(if=90,jf=27,kf=205,lf=17)
DIMENSION x(if),y(jf),h(if),g(jf),c(if,jf,8),omgm(if,jf),omg(if,jf),omg1(if,jf),e(if,jf), &d(if,jf),psim(if,jf),psi(if,jf),psi1(if,jf),u(if,jf),u1(if,jf),ub(if,jf),us(if,jf),&v(if,jf),v1(if,jf),vb(if,jf),vs(if,jf),p(if,jf),cm(lf),c0(lf),c1(lf),phi(lf,kf),&a(lf,lf),ah(lf,kf)!xy:coordinates, h,g:grid spaces, c:coefficients of first- and second-derivatives!omgm:omega^n-1, omg:omega^n, omg1:omega^n+1, d:nabla^2omega, uv:velocities!u:u^n, u1:u^n+1, ub:baru, us:u^*, p:total pressure, cm:c^n-1, c0:c^n, c1:c^n+1!phi:eigenfunctions, a:a or a^-1, ah:hata
COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,Re,n,na
TEXT(x,f1,f2,f3) = -.5*x*(1.-x)*f1+(1.-x*x)*f2+.5*x*(1.+x)*f3
DATA i1,i2,j1,j2, i3,i4,j3,j4/8,14,11,17, 44,50,11,17/ &dt,Re,nf/.04, 200., 400/ ufs,pi/1., 3.14159/!i1,i2,j1,j2:front square cylinder, i3,i4,j3,j4:rear square cylinder!dt:time increment, Re:Reynolds number, nf:final n, ufs:u of free stream
OPEN(20,FILE='OUTPUT.dat')! Grid generation
CALL GRIGEN(x,y,h,g,c,if,jf)! Compute circulations
n=-1
FORALL(i= 1:if,j= 1:jf)psi(i,j)=0.
FORALL(i=i1:i2,j=j1:j2)psi(i,j)=1.
CALL CIRC(h,g,c,psi,if,jf,circ11,circ12)
FORALL(i= 1:if,j= 1:jf)psi(i,j)=0.
FORALL(i=i3:i4,j=j3:j4)psi(i,j)=1.
CALL CIRC(h,g,c,psi,if,jf,circ21,circ22)
WRITE(20,'(//5X A30//4(4X A8, F7.3))')'***** circulations *****', &'circ11 =',circ11,'circ12 =',circ12,'circ21 =',circ21,'circ22 =',circ22! Set up phi,ah,a
CALL OUTER0(x,y,if,jf,phi,ah,a,lf,kf)! Compute initial steady flow
n=0
CALL STEADY(x,y,h,g,c,omg,omg1,d,psi1,u1,v1,if,jf) GOTO 121!********************
101 n=n+1 t=t+dt
IF(MOD(n,20)==1)WRITE(20,61)
na=0 102 na=na+1! Compute vorticity
CALL CALOMG(h,g,omgm,omg,omg1,e,d,u,u1,ub,us,v,v1,vb,vs,if,jf)! Compute stream function
IF(na==1)THEN dpsib1=psib11-psib1 dP1=0.
dpsib2=psib21-psib2 dP2=0.
GOTO 111
ENDIF
CALL CALDP(h,g,c,omg,omg1,u,u1,v,v1,p,if,jf,dP1,1)!correct psi on obstacles
CALL CALDP(h,g,c,omg,omg1,u,u1,v,v1,p,if,jf,dP2,2)
dpsib1=.7*dt*(dP1*circ22-dP2*circ21)/(circ11*circ22-circ21*circ12) !psib^N-psib^N-1
dpsib2=.7*dt*(dP2*circ11-dP1*circ12)/(circ11*circ22-circ21*circ12)
psib11=psib11+dpsib1 psib21=psib21+dpsib2
FORALL(i=i1:i2,j=j1:j2)psi1(i,j) = psib11
FORALL(i=i3:i4,j=j3:j4)psi1(i,j) = psib21
111 CONTINUECALL OUTER(y,psi1,u1,v1,if,jf,ufs,phi,ah,a,c0,c1,lf,kf,dcl,dfh) !correct psi on outer boundary
CALL CALPSI(h,g,c,omg1,psi1,u1,v1,if,jf,respsi)! Compute bound values of omega and diffusion term of omega
CALL OMGB(h,g,omg1,psi1,if,jf)
CALL LAPOMG(c,omg1,d,if,jf)
IF(na==1)WRITE(20,62)n,t

プ ログ ラムIF(na<=10.OR.MOD(na,10)==0)WRITE(20,63)na,respsi,psib11,dpsib1,psib21,dpsib2,dcl,dfh,dP1,dP2
IF(n==1)
GOTO 121
IF(na<10)
GOTO 102
CALL CALDP(h,g,c,omg,omg1,u,u1,v,v1,p,if,jf,dP1,1)
CALL CALLD(h,g,omg,omg1,psi,psi1,u,u1,v,v1,p,if,jf,Cl1,Cd1,1)
CALL CALDP(h,g,c,omg,omg1,u,u1,v,v1,p,if,jf,dP2,2)
CALL CALLD(h,g,omg,omg1,psi,psi1,u,u1,v,v1,p,if,jf,Cl2,Cd2,2)
CALL CPU_TIME(sec)WRITE(20,64)Cl1,Cd1,Cl2,Cd2,sec!********************
121 CONTINUE! Advance time and get predicted values of omg, psi, etc by extrapolation
DO i=1,if DO j=1,jf
omgmm=omgm(i,j) omgm(i,j)=omg(i,j) omg(i,j)=omg1(i,j)
IF(n>3)omg1(i,j)=TEXT(2.,omgmm,omgm(i,j),omg(i,j))
psimm=psim(i,j) psim(i,j)=psi(i,j) psi(i,j)=psi1(i,j)
IF(n>3)psi1(i,j)=TEXT(2.,psimm,psim(i,j),psi(i,j))u(i,j)=u1(i,j) v(i,j)=v1(i,j)
ENDDO ENDDO
IF(n==1.AND.na<40)
GOTO 102
DO l=1,lf
cmm=cm(l) cm(l)=c0(l) c0(l)=c1(l) IF(n>3)c1(l)=TEXT(2.,cmm,cm(l),c0(l))
ENDDO
psib1mm=psib1m psib1m=psib1 psib1=psib11 IF(n>3)psib11=TEXT(2.,psib1mm,psib1m,psib1)
psib2mm=psib2m psib2m=psib2 psib2=psib21 IF(n>3)psib21=TEXT(2.,psib2mm,psib2m,psib2)
IF(n<nf)
GOTO 101
CLOSE(20)
FORALL(i=1:if,j=1:jf)p(i,j) = p(i,j)-(u1(i,j)*u1(i,j)+v1(i,j)*v1(i,j))/2.
STOP
61 FORMAT(/5X 'nt ', 4X 'na respsi
psib1 4Pdpsib1
psib2 4Pdpsib2
2Pdcl7', &'dfh ', 6X ' dP1
Cl1
Cd1
dP2
Cl2
Cd2
CPU-time')
62 FORMAT(5X I4,F8.2)
63 FORMAT(20X I2, F8.4, 2(2X 0PF8.4, 4PF8.3), 2X 2PF8.4, 0PF8.4, 4X F8.4, 18X F8.4)
64 FORMAT(96X 2F8.4,10X 2F8.4, 4X F8.2)
END PROGRAM MAIN
! ********** Generate nonuniform rectangular grid
SUBROUTINE GRIGEN(x,y,h,g,c,if,jf)
DIMENSION x(if),y(jf),h(if),g(jf),c(if,jf,8),iq(if)pi=3.14159
FORALL(i= 1: 7)x(i)=.1*(i-22)-.1*SIN(pi*(i-8)/10.)
FORALL(i= 8:50)x(i)=.2/3.*(i-29)
FORALL(i=51:55)x(i)=.1*(i-36)-.1*SIN(pi*(i-50)/10.)
FORALL(i=56:if)x(i)=1.8+.1*(i-55)
FORALL(j= 1:10)y(j)=.1*(j-13)-.1*SIN(pi*(j-11)/10.)
FORALL(j=11:17)y(j)=.2/3.*(j-14)
FORALL(j=18:jf)y(j)=.1*(j-15)-.1*SIN(pi*(j-17)/10.)
FORALL(i=1:if-1)h(i)=x(i+1)-x(i) h(if)=h(if-1)
FORALL(j=1:jf-1)g(j)=y(j+1)-y(j)! Setup coefs of first- and second-order derivatives

以下省略
ENDSUBROUTINE GRIGEN! ********** Compute circulation
SUBROUTINE CIRC(h,g,c,psi,if,jf,circ1,circ2)
DIMENSION h(if),g(jf),c(if,jf,8),psi(if,jf),o(if,jf),u(if,jf),v(if,jf)
DATA i5,i6,j5,j6/6,16,9,19/ i7,i8,j7,j8/42,52,9,19/

29

30

非圧縮性流れの解法 | !

;法FORALL(i=1:if,j=1:jf)o(i,j)=0.
CALL CALPSI(h,g,c,o,psi,u,v,if,jf,resp)
circ1=0. circ2=0.
DO i=i5,i6-1 circ1 = circ1+( u(i,j5)+u(i+1,j5)-u(i,j6)-u(i+1,j6))/2.*h(i) ENDDO
DO j=j5,j6-1 circ1 = circ1+(-v(i5,j)-v(i5,j+1)+v(i6,j)+v(i6,j+1))/2.*g(j) ENDDO
DO i=i7,i8-1 circ2 = circ2+( u(i,j7)+u(i+1,j7)-u(i,j8)-u(i+1,j8))/2.*h(i) ENDDO
DO j=j7,j8-1 circ2 = circ2+(-v(i7,j)-v(i7,j+1)+v(i8,j)+v(i8,j+1))/2.*g(j) ENDDO
ENDSUBROUTINE CIRC
! ********** Set up phi,ah,a
SUBROUTINE OUTER0(x,y,if,jf,phi,ah,a,lf,kf)
DIMENSION x(if),y(jf),phi(lf,kf),ah(lf,kf),a(lf,lf),wa(lf),iw(lf)k1=if k2=k1+jf-1 !kf=k2+if-1
DO k=1,kf
i=k1-k+1 IF(k>k1)i=1
 IF(k>k2)i=k-k2+1
j=jf IF(k>k1)j=k2-k+1 IF(k>k2)j=1
x1=x(i) y1=y(j) rr=x1*x1+y1*y1 x2=x1/rr y2=y1/rr r=SQRT(rr)phi(1,k) = x2 phi(2,k) = y2 phi(lf,k) = ALOG(2./r)!phi_lf:phi_0
IF(k>k1.AND.k<k2) &THEN ah(1,k) = x2*x2-y2*y2 ah(2,k) = 2.*x2*y2
 ah(lf,k) = x2
ELSE ah(1,k) = 2.*x2*y2
 ah(2,k) = y2*y2-x2*x2 ah(lf,k) = y2 ENDIF
DO l=4,lf-1,2 al=l/FLOAT(l-2)
phi(l-1,k) = x2*phi(l-3,k)-y2*phi(l-2,k)  phi(l,k) = y2*phi(l-3,k)+x2*phi(l-2,k)ah(l-1,k) = (x2*ah(l-3,k)-y2*ah(l-2,k))*al ah(l,k) = (y2*ah(l-3,k)+x2*ah(l-2,k))*alENDDO ENDDO
DO m=1,lf DO l=1,lf a(m,l) = 0. DO k=1,kf
a(m,l) = a(m,l)+ah(m,k)*ah(l,k)
ENDDO ENDDO ENDDO
det=1. epsl=0. CALL MATINV(a,wa,iw,lf,det,epsl)
WRITE(20,'(//5X A40//5X A8/(5X 17E12.4)/)') &'***** Coefs for outer bound cond *****','A^-1 =',a
END SUBROUTINE OUTER0! ********** Compute initial steady flow
SUBROUTINE STEADY(x,y,h,g,c,omg,omg1,d,psi,u,v,if,jf)
DIMENSION x(if),y(jf),h(if),g(jf),c(if,jf,8),omg(if,jf),omg1(if,jf),d(if,jf), &psi(if,jf),u(if,jf),v(if,jf)
COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,Re,n,na! Initial values
pi=3.14159 ufs=1.
rr0=(x(i2)-x(i1))*(y(j2)-y(j1))/pi
FORALL(i=1:if,j=1:jf)psi(i,j)=0.
DO i=1,if cycle_1: DO j=1,jf
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_1
IF(i>=i3.AND.i<=i4.AND.j>=j3.AND.j<=j4)CYCLE cycle_1
x1=x(i)-x(14) y1=y(j) rr=x1*x1+y1*y1 r=SQRT(rr)!in order to avoid rr=0
rra=rr0/rr xa=x1/r ya=y1/r
psi(i,j) = ufs*y1*(1.-rra)u(i,j) = ufs*(1.-rra+2.*ya*ya*rra)v(i,j) =-ufs*2.*xa*ya*rra
ENDDO cycle_1 ENDDO! Compute steady flow

以下省略
ENDSUBROUTINE STEADY! ********** Solve vorticity transport eqn by 1st-order convective-difference method
SUBROUTINE CALOMG0(h,g,c,o,o1,d,p,u,v,if,jf,domg)

プ ログ ラム
DIMENSION h(if),g(jf),c(if,jf,8),o(if,jf),o1(if,jf),d(if,jf),p(if,jf),u(if,jf),v(if,jf)
COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,Re,n,na! Compute vorticity using 1st-order convective-differences
CALL UPWIND(h,g,u,v,o,o1,if,jf)
DO i=2,if ip1=MIN(i+1,if) cycle_1: DO j=2,jf-1
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_1
IF(i>=i3.AND.i<=i4.AND.j>=j3.AND.j<=j4)CYCLE cycle_1
o1(i,j) = o1(i,j)+d(i,j)*dt/Re

以下省略
ENDSUBROUTINE CALOMG0! ********** Solve vorticity transport eqn by 2nd-order convective-differnece method
SUBROUTINE CALOMG(h,g,om,o,o1,e,d,u,u1,ub,us,v,v1,vb,vs,if,jf)
DIMENSION h(if),g(jf),om(if,jf),o(if,jf),o1(if,jf),e(if,jf),d(if,jf),u(if,jf), &u1(if,jf),ub(if,jf),us(if,jf),v(if,jf),v1(if,jf),vb(if,jf),vs(if,jf)
COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,Re,n,na
IF(na>3)
RETURNIF(na>1)
GOTO 10
DO k=1,2
IF(k==1)THEN ib=i1 ie=i2 jb=j1 je=j2 ENDIF
IF(k==2)THEN ib=i3 ie=i4 jb=j3 je=j4 ENDIF
FORALL(i=ib:ie)d(i,jb) = (o1(i,jb)-om(i,jb))/(2.*dt)*Re
FORALL(i=ib:ie)d(i,je) = (o1(i,je)-om(i,je))/(2.*dt)*Re
FORALL(j=jb:je)d(ib,j) = (o1(ib,j)-om(ib,j))/(2.*dt)*Re
FORALL(j=jb:je)d(ie,j) = (o1(ie,j)-om(ie,j))/(2.*dt)*Re
ENDDO
FORALL(i=2:if,j=2:jf-1)e(i,j) = o(i,j)+dt*d(i,j)/2./Re

以下省略
ENDSUBROUTINE CALOMG! ********** Compute boundary values of vorticity by Woods cond
SUBROUTINE OMGB(h,g,o,p,if,jf)
DIMENSION h(if),g(jf),o(if,jf),p(if,jf)
COMMON i1,i2,j1,j2,i3,i4,j3,j4
OB(h,o1,p0,p1) = -o1/2.-3.*(p1-p0)/h/h
DO k=1,2
IF(k==1)THEN ib=i1 ie=i2 jb=j1 je=j2 ENDIF
IF(k==2)THEN ib=i3 ie=i4 jb=j3 je=j4 ENDIF
DO j=jb,jeo(ib,j) = OB(h(ib-1),o(ib-1,j),p(ib,j),p(ib-1,j))!front wallo(ie,j) = OB(h(ie) ,o(ie+1,j),p(ie,j),p(ie+1,j))!rear wall
ENDDO obb=o(ib,jb) obe=o(ib,je)
DO i=ib,ieo(i,jb) = OB(g(jb-1),o(i,jb-1),p(i,jb),p(i,jb-1))!bottom wallo(i,je) = OB(g(je) ,o(i,je+1),p(i,je),p(i,je+1))!top wall
ENDDO o(ib,jb)=(obb+o(ib,jb))/2. o(ib,je)=(obe+o(ib,je))/2.
ENDDO
ENDSUBROUTINE OMGB
! ********** Compute Laplacian omega
SUBROUTINE LAPOMG(c,o,d,if,jf)

省 略
ENDSUBROUTINE LAPOMG! ********** Compute difference of total pressure along a closed curve
SUBROUTINE CALDP(h,g,c,o,o1,u,u1,v,v1,p,if,jf,dP,k)DIMENSION h(if),g(jf),c(if,jf,8),o(if,jf),o1(if,jf), &

31

32

非圧縮性流れの解法 | !

;法u(if,jf),u1(if,jf),v(if,jf),v1(if,jf),p(if,jf)
COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,ReIF(k==1)THEN i5=i1-2 i6=i2+2 j5=j1-2 j6=j2+2 ENDIF
IF(k==2)THEN i5=i3-2 i6=i4+2 j5=j3-2 j6=j4+2 ENDIF
p(i5,j5)=0. j=j5 10 CONTINUE

以下省略
ENDSUBROUTINE CALDP! ********** Correct psi on outer boundary
SUBROUTINE OUTER(y,p,u,v,if,jf,ufs,phi,ah,a,c,c1,lf,kf,dcl,dfh)
DIMENSION y(jf),p(if,jf),u(if,jf),v(if,jf),phi(lf,kf),ah(lf,kf),a(lf,lf), &c(lf),c1(lf),fh(kf),f(lf)COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,Re,n,na
k1=if k2=k1+jf-1 !kf=k2+if-1
DO k=1,kf!compute hatf_ki=k1-k+1 IF(k>k1)i=1
 IF(k>k2)i=k-k2+1
j=jf IF(k>k1)j=k2-k+1 IF(k>k2)j=1
IF(k>k1.AND.k<k2)THEN fh(k) = v(i,j)!upstream bound
ELSE fh(k) = ufs-u(i,j) ENDIF!both side bounds
ENDDO
DO m=1,lf f(m)=0. DO k=1,kf!compute f_m
f(m) = f(m)+ah(m,k)*fh(k)ENDDO ENDDO
10 dcl=0. dfh=0.
DO l=1,lf!compute
dcl=sum|c^n-c^n-1|
c10=c(l) IF(na==1)
GOTO 11c10=c1(l) c1(l)=0.
!compute c_l, dcl=sum|c^N-c^N-1|
DO m=1,lf c1(l)=c1(l)+a(l,m)*f(m) ENDDO
c1(l) = c10+.30*(c1(l)-c10)
11 dcl = dcl+ABS(c1(l)-c10)
ENDDO
DO k=1,kf!compute Psi_k & Psi_,nk(=fhk)
i=k1-k+1 IF(k>k1)i=1
IF(k>k2)i=k-k2+1
j=jf
IF(k>k1)j=k2-k+1 IF(k>k2)j=1
fhk=0. IF(na/=1)p(i,j)=ufs*y(j)cycle_1: DO l=1,lf
fhk
= fhk
+c1(l)* ah(l,k) IF(na==1)CYCLE cycle_1
p(i,j) = p(i,j)+c1(l)*phi(l,k)
ENDDO cycle_1 dfhk=fhk-fh(k)dfh = dfh+dfhk*dfhk/kf
ENDDO
dfh = SQRT(dfh)!dfh=sqrt(sum(fhk-hatf_k)^2/kf)
ENDSUBROUTINE OUTER! ********** Solve boundary value problem of stream function eqn using SOR
SUBROUTINE CALPSI(h,g,c,o,p,u,v,if,jf,resp)!o=omg, p=psi
DIMENSION h(if),g(jf),c(if,jf,8),o(if,jf),p(if,jf),u(if,jf),v(if,jf)
COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,Re,n,na
ni=0 10 ni=ni+1 resp=0.
DO i=2,if cycle_1: DO j=2,jf-1
IF(i>=i1.AND.i<=i2.AND.j>=j1.AND.j<=j2)CYCLE cycle_1
IF(i>=i3.AND.i<=i4.AND.j>=j3.AND.j<=j4)CYCLE cycle_1

以下省略
ENDSUBROUTINE CALPSI! ********** Compute lift and drag coefficients
SUBROUTINE CALLD(h,g,omg,omg1,psi,psi1,u,u1,v,v1,p,if,jf,Cl,Cd,k)プ ログ ラム
DIMENSION

h(if),g(jf),omg(if,jf),omg1(if,jf),psi(if,jf),psi1(if,jf), &u(if,jf),u1(if,jf),v(if,jf),v1(if,jf),p(if,jf)
COMMON i1,i2,j1,j2,i3,i4,j3,j4,dt,ReIF(k==1)THEN i5=i1-2 i6=i2+2 j5=j1-2 j6=j2+2 ENDIF
IF(k==2)THEN i5=i3-2 i6=i4+2 j5=j3-2 j6=j4+2 ENDIF
j=j5 11 Cl1=0. Cd1=0.

以下省略
ENDSUBROUTINE CALLD! ********** Output computational results
SUBROUTINE OUTPUT(o,p,u,v,if,jf)

省 略
ENDSUBROUTINE OUTPUT!***** Matrix inversion by Gauss-Jordan reduction
SUBROUTINE MATINV(a,wa,iw,n,det,e)
DIMENSION a(n,n),wa(n),iw(n)det=1.
e=MAX(e,1.E-5)
FORALL(k=1:n)iw(k)=k
DO k=1,nw=a(k,k) IF(ABS(w)>e)
GOTO 11l=k 12 l=l+1 IF(l>n)STOP 5555
w=a(l,k) IF(ABS(w)<=e)
GOTO 12
FORALL(j=1:n)wa(j)=a(k,j) a(k,j)=a(l,j) a(l,j)=wa(j)ENDFORALL
iw0=iw(k) iw(k)=iw(l) iw(l)=iw0
det=-det 11 det=w*det
a(k,k)=1.
FORALL(j=1:n)a(k,j)=a(k,j)/w
cycle_1: DO i=1,n IF(i==k)CYCLE cycle_1
w=a(i,k) a(i,k)=0.
FORALL(j=1:n)a(i,j)=a(i,j)-w*a(k,j)
ENDDO cycle_1
ENDDO
cycle_2: DO j=1,n IF(iw(j)==j)CYCLE cycle_2
l=j 13 l=l+1 IF(iw(l)/=j)
GOTO 13
FORALL(i=1:n)wa(i)=a(i,j) a(i,j)=a(i,l) a(i,l)=wa(i)ENDFORALL
iw(l)=iw(j)ENDDO cycle_2
END!***** 1st-order upwind interpolation
SUBROUTINE UPWIND(h,g,u,v,f,fs,if,jf)

省 略
ENDSUBROUTINE UPWIND!***** Lagrangian interpolation polynomial
SUBROUTINE LAGINT(h,g,u,v,f,fs,if,jf)

省 略
ENDSUBROUTINE LAGINT

33

非圧縮性流れの解法 | !

34

;法

以下にはこのプログラムについて前記の内部流れと重複する部分は省略し異なる部分のみ詳しく述べる．
サブルーチン GRIGEN では格子形成と中心差分の係数の設定が行われる．サブルーチン circ では式 (12.46)のポテンシャル流れの境界値問題が解かれ，角柱周りの循環が式 (12.47) から計算される．circ mn = Imn ．
1^l ，a = aml  = a;
サブルーチン OUTER0 では，外部ポテンシャル流れの計算のための，phi = l ，ah = alm の値の設定が行われる．外側境界は内部領域の前方と上下の境界からなり，k はこの境界上の格子点に右下か^l ，式 (12.56b) から行列 aml ，更にら順にうたれた格子点番号である．式 (12.52) から l ，式 (12.55) から a1
補助サブルーチン MATINV を用い逆行列 a;lm が求められる．準備段階の最後はサブルーチン STEADY による

初期定常流れの計算で，その初期値は一様流中にある図 12.6 に示すような等価な一つの円柱周りのポテン
シャル流れとして与えられる．

r2
(x y) = Uy 1 ; 02r
ufs = U は一様流の流速，r は前方角柱の中心に原点を持つ極座標，r0 は等価な円柱の半径である．ここでは内部領域の流れのみ定常粘性流れとして計算される．
非定常流れの各時間ステップにおける第 N 近似計算は次の手順で行われる．

!

! b !  !  ! !b ! r !
2! に関する計算は内部流れと同様である． に関する計算では，SOR 法による計算の過程で角柱表面と外側境界の境界値の修正が行われる．角柱表面の  の値は，N = 1 では前ステップの最後に時間外挿された予測値がそのまま用いられ，N 2 ではサブルーチン CALDP より各角柱まわりの圧力残差値 Pn = ( = 1 2)，+1 2

式 (12.48) から  ( = 1 2)，式 (12.49) から修正値  が決定される．外側境界の  の値はサブルーチ+1
ン OUTER の中で計算され，N = 1 では前ステップの最後に時間外挿された c1(l) = cn
l を基に求められる．
P
N 2 では式 (12.54c) から fh(k) = f^k ，式 (12.56c) から f(m) = fm ，cn+1 = L a;1 fm ，式 (12.51a)l

m=0 mln+1 が順に求められる．なお c1(lf) = cn+1 ．
から  の境界値 psi(i,j) = ij
0

図 12.6.2 は計算結果で，渦度分布と流線のフローパターン，２つの角柱に作用する揚抗力の時間変化を示す．フローパターンは内部領域のみを示したもので半周期をやや超える時間にわたるものである．渦度分
布は強力な周期的渦の発生を示している．角柱表面の流れ関数の値はかなり変化し (前方角柱では後方角柱では

0:041，

0:118 の範囲)，また流線は角柱の排除効果によって外側上下境界のところでかなり膨らんでいることも分かる．下流境界近傍を除けば不自然さは見られない．前方角柱の抗力は一定であるが，後方角柱の抗力は変動し ，また後方角柱の揚力も大きく変動しその振幅は抗力にほぼ匹敵する．

j j < 0:005 で，まも j j < 0:00004 j j < 0:00010 で非常に小さい．
た角柱表面の流れ関数の補正値  =  ;
PL
N N ; j は 0.001 の
外部領域の  の固有関数展開における係数 cl の補正値のノルム dcl =kclk= l jcl ;cl
大きさでこの係数は十分収束している．外側境界における流れ関数の接続については，関数自身は連続である
 PK
=
，
がその法線微分値は不連続になる．その不連続の程度を表すノルムを dfh =k f^k= f k (f^k ) g=K
PL^^^lk ; fk で定義すれば，その値は 0:16  0:29 となりかなり大きいものになる．なお
ただし fk = l cl a
k の総和の範囲を 51  155 に取り下流部分を除けばこのノルムの大きさは 0:040  0:043 となる．下流部
最後に計算の収束性について述べる．! と  に関しては前項と同様である．圧力残差値 P(N )(N ;1)

1

2

=0

( )

(

=1

=0

分を除けば外側境界での流れ関数の接続は良好といえる．

1)

2

1 2

プ ログ ラム

図 12.9: 一様流中の２つの正方形柱まわりの流れ

35

36

非圧縮性流れの解法 | !

;法+1 の計算のダンピング係数 0.3 に関連して少し駄弁を弄したい．ま以上でこの項の説明は終わるが，cnl

ず手始めに簡単な２つの代数方程式 f (x y ) = 0 g (x y ) = 0 を共に満足する解 x y を求める問題を考える．これらの方程式が xy 面上で図 12.10(a) のようになる場合には，x0 を f (x y ) = 0 に入れて y を求め，

これを g (x y ) = 0 に入れて x を求める手順を繰り返しても解は得られない．x0 を先に g (x y ) = 0 に入れy を求めるべきで，なお f (x y) = 0 から求めた y に f (x y) = 0 と g(x y) = 0 から求めた y の差 y に 1よりも大きい適切な係数を掛けたものを加えれば少ない反復で収束解を得ることができる．(b) 図の場合には，x0 を先に f (x y ) = 0 に入れて y を求めれば良いが，これに差 y に 1 よりも小さい適切な係数を掛けたものを加えれば，効率よく収束解を得ることができる．なおこの種の問題で，初期値 a(0) = 1 から収
束解 a = 2 に持って行くときにダンピング係数 0.5 をかすのは，a(1) = 1:5 a(2) = 1:75 a(3) = 1:875 : : :のように持って行こうというのではなく，a(1) = 3:2 になるときにダンピングをかけて a(1) = 2:1 にしようというのである．

y- -
J 



J
?
g (x y ) = 0

f (x y ) = 0

x0

(a)

y

@@
@
@
-@@ ?
@J@@@g (x y ) = 0

f (x y ) = 0x

図 12.10: 加速係数とダンピング係数

x0

(b)

x

本プログラムの場合には，f (x y ) = 0 と g (x y ) = 0 は内部領域の  のポアソン方程式と外部領域の のラプラス方程式に対応し ，また x y は外側境界における  とその勾配 @=@n に対応する．解曲面の様子を頭の中に描けばあるいは紙上に書けばダンピング係数の必要性は分かっても，その値は予め与えることは困難で試行錯誤によって決めらる．

37

曲線座標格子への拡張! ;  法の曲線座標格子への拡張

12.7

前節のプログラムは，1970 年代末に開発したものを本文の説明に合わせ分かりやすく書換えたもので，計算結果もパソコンで計算し直したものである．開発当時はまだコンピュータ環境が悪く，大型計算機センターのコンピュータを利用し，CPU-time 1 時間程度ずつに分割して計算しなければならず，turn-around-time
も数日からときに数週間かかり，この程度の計算が限界であった．本節にはこの解法の曲線座標格子への拡張について述べる．
物理空間のデカルト座標系を xy ，一般曲線座標系をとする．本節ではまずデカルト座標系から一般曲線座標系への変換について述べる．曲線座標格子は物体形状に合わせ，格子間隔が滑らかに変化するように，できるだけ直交に近いように，また流れの急激に変化するところが密になるように作られる．数値計算
は通常すべて，写像空間内の =  = 1 の正方形格子に対して行われる．

12.7.1 デカルト 座標系から一般曲線座標系への変換
デカルト座標と一般曲線座標の変換の測度 (metric) x  : : : と

"

x x
y y

# "

#

x
x

"

1 0
=
0 1
y

#

x

 : : : の間には次の関係が成立する9 ．

y(12.57)

また変換のヤコビアン J はx x 

J = = x y ; x y y y (12.58)となる．これらの関係から

x = J y  x = ;J y  y = ;J x  y = J x

x に関する微分を に関する微分に書き換えれば次のようになる10 ．
@= @ j @  r = (r j ) @@@xi @xi @ j
j
具体的には

@
1 y @ ; y @ =@x J  @@(12.59)(12.60)@1 ; x @ + x @ 
=


@y J
@
@

また一般に

@ @j
J
= 0
@ j @xi@
Jr j = 0@jであるから次の関係が導かれる．

@ @@  @j 
@
=
J j
=
J
@xi
@xi @ j @ j @xi@
Jr =J (r j )@j@  jk @ 
Jg
J r2 =@j@k
J(12.61a)(12.61b)(12.61c)= ( ) = ( ( ) ( )) であるから，この式を y で微分すれば 0 = y x +y x なる関係が得られ，式 (12.57)
Einstein(Einstein cnovention) が用いられる．すなわち同じ添字が繰り返されるときには ai bi = a1 b1+a2 b2

9 例えば xx  x  x y   x yはこれらの関係を纏めたものである．
10 ここでは
総和規約のように解釈する．

非圧縮性流れの解法 | !

38

;法

ただし (g ij ) は測度テンソルで，その成分は次のように定義される．

gij = ri rj =@i @j@xk @xk(12.62)流線に沿う微分は次のように書き換えられる．

u r = ui

@
@j @
@
= ui
= Uj
@xi
@xi @j@j(12.63)

ただし Uj は反変速度 (contravariant velocity) の j 方向成分である．反変速度 U は  空間の流速で 11 ，物理空間の流速 u との間には次の関係がある．

@j@xu  ui = i Uj
@xi i
@j

Uj =

(12.64)

流れの渦度 ! = vx ; uy は，式 (12.61a)，(12.64) を用いれば

!=

1 @ ;

J @



J (x y ; y x )U + J (x y ; y x )V +

1 @ ;

J @

J (x y ; y x )U + J (x y ; y x )V

更に式 (12.59) を用い J (x y ; y x ) = y y + x x  g21 のように置けば次式が得られる．

@@(g21 U + g22 V ) ;(g U + g12 V )@@ 11
ただし (gij ) は測度テンソルで，(gij );1 = (g ij )，その成分は次式で定義される．
  J! =gij =(12.65)@xx @xx @xk @xk
=@i @j @i @j(12.66)

ヤコビアン J は 物理空間の局所面積]/写像空間の相当面積] を表すので， は写像面上の流れの渦度である．物理空間の流れは２次元で高さ方向に一定であるが，ここで高さ方向の間隔 1 の２平面間の流れに着目することにする．写像面上の流れに流れ関数を定義するために，２つの流れの相当する面分を通る流量を合わせることにする．このとき写像面に付加される高さはヤコビアン J になり12 ，流れ関数は次のように定義できる．

@
= JU
@

12.7.2

@
= ;JV
@

(12.67)一般曲線座標系の基礎方程式

連続方程式 J r u = 0 は式 (12.61b) と (12.64) を用いれば次のようになる．

@@
JU + JV = 0@@(12.68)

また NS 方程式は，Ui = (ri ) u であるから，式 (12.1) を書換えた NS 方程式@uu+ r uu = ;rp ;  r@t

(12.69)+11 反変速度が写像空間の流速であることは次のように説明できる．ある時間に点 x にある流体粒子は単位時間後に点 x u に移動
する．点 x の曲線座標を  とすれば点 x u の曲線座標は  u r となる．したがって反変速度 U  u r は写像空間で流体粒子が単位時間に移動する距離すなわち流速であることが分かる．
12 このように写像面の流れに高さを与えることは，渦度の強さが流体の伸縮に比例して変化するという物理的性質，渦度  の定義
とも両立する．

++曲線座標格子への拡張

39の左から J ri  を演算すれば，対流項は式 (12.61b) を用い

J ri 

 @u
u

@t



@Ui
@
+ ri 
JU u
@t
@j j
@Ui
@
@
=J
+
JUj Ui ; u  JUj ri@t @j@j+ r uu = J

また圧力項と粘性項は

@p

@j

J ri r = J ri  !yJ ri rp = Jgij

;

!x



J r r = J (x !y ; y !x ) = y !y + x !x = ! 
J r r = J (x !y ; y !x ) = ;y !y ; x !x = ;!
結局 写像空間の反変速度の NS 方程式は次のようになる．
; @
; @
@
@
@
@
JUU + JV U = Ju U + V
 + Jv U + V

@
@
@
@ x
@
@ y
11
;Jg
p ; Jg12 p ;   
; @
; @
@
@
@
@
JVt + JUV + JV V = Ju U + V
x + Jv U + V

@
@
@
@
@
@ y
21
;Jg
p ; Jg22 p +  JUt +(12.70a)(12.70b)

渦度輸送方程式は式 (12.4) に式 (12.61c)(12.63) を用いればJ@!
@  jk @! 
@!
+ JUj
=
Jg
@t
@j
@j
@k

(12.71)となる．また対流差分法では，この式を  t 空間内の１つの流跡線に沿って積分した式

!

n+1
=
ij

! +

Z tn+1tn@  jk @! 
Jg
dt
J @j
@k
1

(12.72)が用いられる．ただし添字 は流跡線と t = tn 面の交点，線積分の起点を示す．流れ関数の方程式は渦度の式 (12.65) に流れ関数の定義式 (12.67) を代入することによって導かれる．

@  11 @@  @  21 @@ Jg+ Jg 12+Jg+ Jg 22
= ;
@
@
@
@
@
@

(12.73)

この式の導出には測度テンソル成分の関係 g11 = J 2 g 22  g12 = ;J 2 g 21  g21 = ;J 2 g 12  g22 = J 2 g 11 が用いられた．上式は式 (12.7) を式 (12.61c) を用い一般曲線座標系へ書き換えたものに他ならない．


次に圧力一価の条件を一般曲線座標系のものに書換える．x 空間の微小距離 dss = dx dy が  空間の
微小距離 d = d d に対応するものとする．式 (12.35) の各項は次のように書換えられる13 ．

u  dss =@xx@xx
U  d = gjk Uj dk = (g11 U + g21 V )d + (g12 U + g22 V )d@j j @k k=

13 念のため dssx =@k dk のようになるわけを説明する． 空間内のある点における x を x0 とする．この点から d だけ離れ@x

た点の x はテイラー展開により x0 dx
x0x d
x0 @x=@d @x=@d となる．y に関しても同様の関係が得られ，
これより上式の成立することが分かる．

+ = +grad= ++非圧縮性流れの解法 | !

40

;法

@y @x @x @y ;
U d = (x y ; y x ) ! (V d ; Ud)
@j @k @j @k j k@j @y  @!
j @x
r !  dss = !y dx ; !x dy = @;d
@y @k @x @k @j k= (y x ; x y )! +(y x ; x y )! d + (y x ; x y )! +(y x ; x y )! d= J (g 12 ! + g 22 ! )d ; J (g 11 ! + g 21 ! )d

u !  dss = !(vdx ; udy) = !

式 (12.35) にこれらの関係を用いれば一般曲線座標系における圧力一価の条件式が得られる．

I

P = ; (ut ;u ! +  r ! )  dssI=;g11 Ut + g21 Vt ; GV ! + J (g12 ! + g22 ! ) dI; g12 Ut + g22 Vt + GU! ; J (g11 ! + g21 ! ) d = 0
ただし G = x y(12.74);y x ，また積分は  空間内で，流れの中に取られる物体を一巡する閉曲線に沿って行われる．
次に流れの中にある物体に作用する流体力 F =

Fx Fy



を， 空間内で物体を一巡する閉曲線に沿

x 空間内の物体を左周りに一巡する閉曲線に沿う微小距離を上記のようにう積分から求める式を導く．
dss = dx dy とおけば，この曲線の外向き単位法線ベクトル n は，




@x
@y
;
n = dy ; dx  n ds =
dk
@k
@k
ds
ds
となる．また式 (12.39) の各項は次のように書換えられる．

n  t ds = ; dx dy t = ; x d + x d y d + y d tn  uu ds = (u dy ; v dx) u v = G u v (;V d + U d)
n p ds = dy ;dx p = y d + y d ;x d ; x d p
n ! ds = ; dx dy ! = ; x d + x d y d + y d !

式 (12.39) にこれらの関係を用いれば，次の物体に作用する流体力を  空間の閉曲線上で計算するための式が導かれる．

Fx =
Fy =

I

I

(x t + GuV ; y p + x ! )d +(y t + GvV + x p + y ! )d +II(x t ; GuU ; y p + x ! )d(12.75a)(y t ; GvU + x p + y ! )d(12.75b)


n

XXXX





X
y
X

s
ds

XX

41

曲線座標格子への拡張

12.7.3

差分法による解法

上記の一般曲線座標系の基礎方程式を  空間内の正方形格子に対し差分法で解く．ここでは単位正方形格子， =  = 1 を用いることにする．

式 (12.72) は対流差分法を適用すれば対流差分法による渦度輸送方程式の初期値問題の解法：!ijn+1 = ! +  r2 !t(12.76a)のようになり，予測子修正子法で解かれる．ただしの付いている量は流跡線に沿う積分領域の平均値で，n+1
n，修正子の計算では f = (f + fij)=2 のように置かれる．また予測子の計算では f = fijの付いている

起点の値は時間 t における格子点の値から補間によって求められ，起点の位置  はn

 =  ij ;Z tn+1
U dt =  ij ; U t

(12.76b)tn

から求められる．ラプラス差分演算子は  空間では次のようになる．

(J r2 f )ij = g~i11~i11+1=2j (fi+1j ; fij ) ; g;1=2j (fij ; fi;1j )+ g~i12~i12+1=2j (fi+1j +1 + fij +1 ; fi+1j ;1 ; fij ;1 )=4 ; g;1=2j (fij +1 + fi;1j +1 ; fij ;1 ; fi;1j ;1 )=4
21
21+ g~ij~ij+1=2 (fi+1j +1 + fi+1j ; fi;1j +1 ; fi;1j )=4 ; g;1=2 (fi+1j + fi+1j ;1 ; fi;1j ; fi;1j ;1 )=4

22
22+ g~ij~ij+1=2 (fij +1 ; fij ) ; g;1=2 (fij ; fij ;1 )(12.76c)

ただし g~ij = Jg ij ．予測子修正子法の計算手順と補間については 12.3 節と上記プログラムを参照されたい．

SOR 法による流れ関数方程式の境界値問題の解法： 式 (12.73) は 2 次中心差分を用いれば次のように書換えられる．

g~i11+1=2j (i+1j ; ij ) ; g~i11;1=2j (ij ; i;1j )22
22+ g~ij~ij+1=2 (ij +1 ; ij ) ; g;1=2 (ij ; ij ;1 ) =; g~i12+1=2j (i+1j+1 + ij+1 ; i+1j;1 ; ij;1 )=4 + g~i12;1=2j (ij+1 + i;1j+1 ; ij;1 ; i;1j;1 )=4
21
21; g~ij~ij+1=2 (i+1j +1 + i+1j ; i;1j +1 ; i;1j )=4 + g;1=2 (i+1j + i+1j ;1 ; i;1j ; i;1j ;1 )=4; ij

右辺の @ (~
g12  )=@ @ (~g21  )=@ に由来する項の値は SOR 法の計算の過程で随時新しいものに置換えら
れる．
流跡線

H




+1

! n
ij




t

6
;;


;;
;;
;;
;
;
;

;
;
;
;
!;
;
;
;;H t 面
;  ;










ij

図 12.11:  t 空間内の積分路n42

非圧縮性流れの解法 | !

;法Woods の条件による渦度の境界値の計算： x 空間内の平板境界上の点を添字 B，そこから法線方向に
n 離れた点を添字 1 を付して示すことにする．Woods の条件は境界上の渦度 !B の値を既知の !1  B  1の値から次式によって決定するものである．

3 u + 3 ( ;  ) + O(n2 )!B = ; 12 !1 + n
(12.77)B n2 B 1
ただし uB は平板境界がそれ自身に平行に移動する速度である．いま平板境界方向の格子線に沿って =const.
これらを横断する格子線に沿って
=const. とする．図 12.12 に示すように，点 B から隣接格子点までの距
離を dss = dx dy とすれば ，点 B における単位法線ベクトルは n = ;dy=ds dx=ds となる．また
法線方向の微小長さ n = j; (dy=ds)x+(dx=ds)y j，格子点 (i1) から法線までの格子線に沿う微小長さ
s = (dx=ds)x+(dy=ds)y となる．これより点 1 の !1  1 の値を格子点 (i;1 1) (i1) (i+1 1) から補間によって求める際に必要な s に対応する  の値は次のようになる14 ．
s  s = 1 q(x ; x )2 +(y ; y )2 = 1 ; i1 = si+11 i;11
2 i+11 i;11

図 12.12: 曲線座標格子の Woods の条件 (平板境界)
次にこの条件を滑らかな湾曲境界の場合に拡張する．点 B における接線方向の座標を s，法線方向の座標を n とする．Taylor 展開

1 = B +n(n )B + 21 n2 (nn )B + 16 n3 (nnn )B +   に静止壁面の条件 (n )B = 0，ここでの流れ関数の方程式 (nn )B = ;!B ; (ss )B ，その微分 (nnn )B =;(!1 ; !B )=n ;f(ss)1 ; (ss)B g=n +    を代入すれば，この場合の Woods の条件式が次のように得られる．

(12.78)!B = ; 12 !1 + n3 2 (B ; 1 ) ; (ss)B ; 21 (ss )1 +   
最後の 2 項は壁面曲率の効果を表わす．これらの項は下記のようにして求めることができる．なおこの場合には ds n  等は dx = (xi+10 ; xi;10 )=2 dy = (yi+10 ; yi;10 )=2 と置いて求められる．
(ss)1 の近似値は概略次のように求められる．まず点 (i+1 1) と点 (i ; 1 1) を結ぶ線と法線との交点 Aの座標 xA  yA を求める．点 A から点 1 までの距離 Dn1 = n1 ; nA を求める．nC = n1 +Dn1 ．次に法線上の (n) の値を 2 次式 (n) = B +cn2  c = (1 ; B )=n2 で近似し ，点 C における C を決定する．
D1 = C ; 1 ．点 1 における ss の近似値は(ss)1 = (i+11 ; 2i1 +i;11 +2D1 )=s2
14 dx

=

xi+10;xB  dy= +1 0 ;yiyB  ds=

p 2+ 2dx

dy  x=xi1;

xB  y

= 1;
yi

yB

．s  は正にも負にもなる．

曲線座標格子への拡張

43

図 12.13: 曲線座標格子の Woods の条件 (湾曲境界);1) と (i+1) 上の平均挙動に近いという前提で作られた．(ss )B の近似値も同様に決定できる．境界では (ss )B = 2DB =ds2 ，境界面が凸の場合には DB 正で曲率も正，凹面では共に負になる．
反変流速 U の計算： 式 (12.67) からJU = U~ = (ij+1 ; ij;1 )=2 JV = V~ = ;(i+1j ; i;1j )=2

この式は点 1 の曲率 ss が点 (i1) のものに近く，法線上の  の挙動が線 (i

閉曲線に沿う圧力差の計算：

物体を一巡する閉曲線は通常  空間内の格子線上に沿って取られる．し

たがって式 (12.74) の線積分は  か  のみの積分になる．物体に作用する流体力の計算に関しても同様のことが言える．
以上本項には デカルト座標格子の場合と相違する 主な曲線座標格子の差分式を示した．全体の計算の流れやプログラミングについてはデカルト座標格子の場合と同じになるのでここでは省略する．