C ---------------------------------------------------------------- C ２次精度中心差分法による平行平板間乱流の直接シミュレーション C 公開版 Ver.1.2 （次バージョンでは不等間隔格子の補間を修正する予定） C ---------------------------------------------------------------- C このファイルは、プログラム名 DNS37AB2.F を使用するに当たって C 必要な情報をコメント行として書き込んだものです。 C *************************************************************** C * DNS37AB2 Ver.1 1994.09.13 (4th order accuracy) * C * Ver.3B 1997.07.12 (2nd order; corrected) * C * by Takeo KAJISHIMA, Dept. Mech. Eng., Osaka University * C * 2-1 Yamada-oka, Suita, Osaka 565-0871, Japan * C * TEL: +81 6 879 7249, FAX: + 81 6 879 7250 * C * E-mail: kajisima@mech.eng.osaka-u.ac.jp * C *************************************************************** C C 「平行平板間乱流の直接数値シミュレーション」 C 空間差分近似：２次精度中心差分法 C 時間発展解法：２次精度Adams-Bashforth（非線形項と粘性項） C １次精度後退Euler（圧力勾配項と連続の式） C C * エネルギー収支、レイノルズ応力収支、各種乱流統計量を算出する C サブルーチンは除去してあります。 C * このプログラムはベクトル機や並列機のために調整されていません。 C C * 多くのサンプルが必要な場合： C 初期値は一組（ある瞬時の u,v,w,p）だけを提供しますが、 C 計算領域が有限なので、瞬時の平均流れ場はゆったりと変動します。 C 多くのサンプルがほしい場合には、長い時間走らせて、別の状態を C 多数作って下さい。無次元時間で 1〜2（時間刻み DT=1.D-4 ならば C 10,000〜20,000ステップ）経過すれば、初期値とは相関のない状態に C 変わると思います。 C C *（あまり勧めませんが）格子数、格子配置を変更する場合： C その配置における統計的定常状態に達するまで、十分な時間進行が必要です。 C 格子配置をどの程度変更するかにもよりますが、無次元時間で5〜10 C （時間刻み DT=1.D-4 ならば 50,000〜100,000ステップ）は必要でしょう。 C C * 時間平均流れ場を求めようとする場合： C 最低次の統計量である 断面平均速度 UME が非常にゆったり変化するので C これが数回変動する時間で平均する必要があります。 C また、時間平均する長さが十分であるかどうかのチェックは、 C (1) せん断応力分布（粘性応力とレイノルズ応力の和）が -1 から +1 を C 結ぶ直線になっているか（与えた圧力勾配と釣り合っているか） C (2) 本来 0 になるべき横断方向速度変動(w)の３次モーメント(skewness) C が例えば他の方向のそれに比べて十分に小さいか C を調べればわかります。許容範囲は、主観的ですが誤差１％程度を目安と C してはどうでしょうか。経験的には 無次元時間で 5〜10 位の幅での時間 C 平均が必要でしょう。 C （計算機資源に余裕のある人は、格子点と計算領域を増やし、 C 瞬時のサンプル数を大きくすれば、より短い時間平均でも大丈夫です） PROGRAM DNS37AB2 PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*3 CT(0:100) CHARACTER*4 FILE1 COMMON /FILE1/FILE1 /RMS/URMSX,URMSC,VRMSC,WRMSC COMMON /RET/RET /DT/DT /STEP/ISTEP /TSTEP/TSTEP /START/ISTART COMMON /DIV/DIVX /COU/COMX /ENE/ENE /UME/UME,UMX COMMON /ERRP/POIERR /ITRP/ITRP DATA CT/'000', &'001','002','003','004','005','006','007','008','009','010', &'011','012','013','014','015','016','017','018','019','020', &'021','022','023','024','025','026','027','028','029','030', &'031','032','033','034','035','036','037','038','039','040', &'041','042','043','044','045','046','047','048','049','050', &'051','052','053','054','055','056','057','058','059','060', &'061','062','063','064','065','066','067','068','069','070', &'071','072','073','074','075','076','077','078','079','080', &'081','082','083','084','085','086','087','088','089','090', &'091','092','093','094','095','096','097','098','099','100'/ C FILE1 はファイル名の先頭の4文字を与えています。 C ファイル名は 'dns2'+CT+'u/v/w/p.d' です。 FILE1='dns2' C RET は流路幅 H と 壁面平均摩擦速度 U_tau に基づくレイノルズ数です。 C これを固定することは、圧力勾配を与えたことに対応します。 RET=300.D0 C DT は時間刻み（流路幅 H と 壁面平均摩擦速度 U_tau で無次元化）です。 C このプログラムでは格子配置を変更しない限り DT=2.0D-4 でも平気です。 C 非定常状態の解析には時間解像度を確保するために DT=1.0D-4 がよいでしょう。 DT= 1.0D-4 C 時間刻みは 壁指標では RET*DT で換算して下さい。 C DT=1.0D-4 のときには DT+=0.03 です。 C このプログラムは ISKIP=40 ステップ（壁指標での時間間隔は C DT*300*40 = 1.2）ごとに ISAVE=100 回データを保存するようになっています。 C 保存する時間間隔、保存回数は適当に調整してください。 ISKIP=40 ISAVE=100 C 時間進行させるステップ数 ISTOP は保存するステップ間隔と保存回数の積です。 ISTOP=ISKIP*ISAVE C プログラム実行後の時間ステップ数 ICOUNT と 保存回数 ICT を 0 にリセット。 ICOUNT=0 ICT=0 C 格子配置を決め、差分係数などを求めます。 CALL SBRMSH C 初期値(u,v,w,p)を読み込みます。 C CT(0)='000' だから 初期値のファイル名は 'dns4000u/v/w/p.d' となります。 C このファイルが同一ディレクトリーにないとプログラムは正常に動きません。 CALL SBRDTR(CT(ICT)) C 初期値の状態を出力しておきます。 WRITE( 6,1000) CALL SBRUMR CALL SBRST1 CALL SBRCHK WRITE( 6,1100) ISTART,TSTEP,ITRP,POIERR,DIVX,COMX & ,UME,UMX,INT(RET*UME),INT(RET*UMX),ENE,URMSX,URMSC,VRMSC,WRMSC 1000 FORMAT(3X,4HSTEP,7X,1HT,1X,4HITRP,4X,6HPOIerr,4X,6HDIVmax &,2X,6HCOUmax,3X,5HUmean,4X,4HUmax,4X,3HRem,4X,3HRex &,4X,6HEnergy,3X,5HU'max,5X,3HU'c,5X,3HV'c,5X,3HW'c) 1100 FORMAT(I7,F8.4,I5,2E10.2,3F8.3,2I7,F10.5,4F8.4) C================================================================= C ここが時間進行ループの先頭です。 10 CONTINUE C 非線形項（対流項）を計算します。 CALL SBRCON C 粘性項を計算します。 CALL SBRVIS C 部分段階速度を計算します。プログラムの実行開始直後には 現時点での C 流れ場しかわからないので 陽的Euler法、2回目以降は Adamas-Bashforth法 C を適用します。 IF(ICOUNT.GE.1) THEN CALL SBRPRE(1.5D0,-0.5D0) ELSE CALL SBRPRE(1.0D0, 0.0D0) ENDIF C 圧力方程式の右辺を計算します。 CALL SBRRHP C 圧力のポアソン方程式を逐次過緩和(SOR)法で解きます。 C 反復回数を IPSTOP=20 回としておきますが、div u の出力結果を見て C 収束精度に不安があるなら反復回数を増やしてください。 IPSTOP=20 CALL SBRSOR(IPSTOP) C 圧力勾配項を加えて時間進行を完了します。 CALL SBRCOR ICOUNT=ICOUNT+1 ISTEP=ISTART+ICOUNT TSTEP=TSTEP+DT C 最初の 10 ステップおよびそれ以降の 10 ステップおきに流れ場を C チェックするため、瞬時の統計量を出力します。 IF(MOD(ISTEP,10).eq.0.OR.ICOUNT.LE.10) THEN CALL SBRUMR CALL SBRST1 CALL SBRCHK WRITE( 6,1100) ISTEP,TSTEP,ITRP,POIERR,DIVX,COMX & ,UME,UMX,INT(RET*UME),INT(RET*UMX),ENE,URMSX,URMSC,VRMSC,WRMSC ENDIF C ISTEP: 時間ステップ数 C TSTEP: 無次元の経過時間 C ITRP: SOR法の反復回数 C POIERR:ポアソン式の誤差（誤差のノルム／右辺のノルム） C DIVX: 全格子における連続条件の誤差 div u の最大値 C COMX: 全格子におけるクーラン数の最大値 C UME: 断面平均流速 C UMC: 中央断面での平均速度 C RET*UME: バルク速度によるレイノルズ数。 C RET*UMC: 中心断面の速度によるレイノルズ数。 C ENE: 全流れ場での運動エネルギーの合計 C URMSX: 主流方向速度変動の最大値 C URMSC,VRMSC,WRMSC: 中心断面での各方向速度の変動強度 C ISKIP ステップ毎にデータ(u,v,w,p)をファイル 'dns4???u/v/w/p.d' に C 保存します。ファイルはサブルーチン内で新規に作成されます。 C 同一名のファイルが存在していると、プログラムの実行を停止します。 IF(MOD(ISTEP,ISKIP).EQ.0) THEN ICT=ICT+1 CALL SBRDTS(CT(ICT)) ENDIF C 指定のステップ数に達したら、時間進行ループを抜け出します。 IF(ICOUNT.GE.ISTOP) GO TO 20 C 開始後 10 ステップが過ぎた後、div u の最大値が大きすぎたら C 警告を発して実効を停止します。 IF(ISTEP.GE.10.AND.DIVX.GE.1.D+2) GO TO 30 GO TO 10 C================================================================= C ここが時間進行ループの最後尾です。 20 CONTINUE STOP 30 CONTINUE WRITE(6,*) 'Stopped because of the numerical instability' STOP END C ========================================================================= C 原則的には、これ以降のサブルーチンの内容を変更しないで下さい。 C 変更の余地のある箇所だけ（たとえば格子配置）にコメントを入れておきます。 C ========================================================================= C ********************************************************************* C * SBR. MSH : MESH PARAMETERS * C ********************************************************************* C 格子配置を決め、差分係数などを求めます。 SUBROUTINE SBRMSH PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /H/HX,HZ /Y/YV(0:NY),YP(0:NY1) COMMON /D/DX,DZ /DY/DY(NY) /RET/RET /DT/DT CALL MESHY CALL MESHXZ WRITE(6,1000) NX,HX,DX,RET*HX,RET*DX,NZ,HZ,DZ,RET*HZ,RET*DZ &,NY,1.,DY(1),DY(NY/2),RET*DY(1),RET*DY(NY/2),YP(1),RET*YP(1) &,INT(RET),DT 1000 FORMAT(/10(1H*)," IN SBR.MSH ",10(1H*)/ &/5X,"NX=",I4,5X,"HX=",F7.3,3X,"DX=",F7.5 &,5X,"HX+=",F7.1,3X,"DX+=",F7.3 &/5X,"NZ=",I4,5X,"HZ=",F7.3,3X,"DZ=",F7.5 &,5X,"HZ+=",F7.1,3X,"DZ+=",F7.3 &/5X,"NY=",I4,5X,"HY=",F7.3,3X,"DYmin=",F7.5,3X,"DYmax=",F7.5 &/30X,"DYmin+=",F6.3,3X,"DYmax+=",F6.2 &/5X,"YP(1)=",F7.5,3X,"YP(1)+=",F6.2 &/5X,"RET=",I6/5X,"DT=",F9.5) CALL MDQ1VN CALL MDQ1VD CALL MDQ1P CALL MDQ2V CALL MDQ2PV CALL MDQ2PP RETURN END C--------------------------------------------------------------------- C----- Mesh generation ..... for X, Z directions ------------------ SUBROUTINE MESHXZ PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /H/HX,HZ /D/DX,DZ /RET/RET COMMON /IP/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(NZ) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA COMMON /CMPX/CM1X,C00X,CP1X /CMPZ/CM1Z,C00Z,CP1Z C 主流(x)方向および横断(z)方向の格子幅は、それぞれ壁指標で C DX+=18, DZ+=9 となるように設定されています。 C したがって計算領域は 壁指標でそれぞれ 格子数を乗じて C HX+=1152, HZ+=576 となります。 C これで壁近傍のストリーク間隔 Lz+ = 約100 が再現されます。 C もっと粗くすると、ストリーク間隔がなまってしまいます。 C もっと細かくすると、格子数 64 では計算領域が不足し、周期条件の C 影響が大きくなってしまいます。 DX=18.D0/RET DZ= 9.D0/RET HX=NX*DX HZ=NZ*DZ C 計算機資源にもう少し余裕のある場合（例えば3倍くらい容量を確保 C できる場合には C NX= 96, (DX+=12.), DX=12.D0/RET C NZ=128, (DX+=4.5), DX=4.5D0/RET C とすることを勧めます。 C その際、時間刻みは DT=0.75D-4 くらいにするのが無難でしょう。 C 現在の格子設定では、２次モーメント（速度変動強度）までは大丈夫ですが、 C ３次（skewness）や４次（flatness）に対しては解像度が不足しています。 C 限られた計算機容量では、横断(z)方向の格子を細かくして解像度を高める C のが最も効果的で、壁に垂直(y)方向の格子を細かく、あるいは壁に寄せても C あまり効果がありません。 BXA=1.D0/DX BZA=1.D0/DZ TXA=BXA/2.D0 TZA=BZA/2.D0 C --- List vectors to identify the neighbouring point DO 50 I=1,NX IP(I)=I+1 IM(I)=I-1 IF(IP(I).GT.NX) IP(I)=IP(I)-NX IF(IM(I).LT. 1) IM(I)=IM(I)+NX 50 CONTINUE DO 60 K=1,NZ KP(K)=K+1 KM(K)=K-1 IF(KP(K).GT.NZ) KP(K)=KP(K)-NZ IF(KM(K).LT. 1) KM(K)=KM(K)+NZ 60 CONTINUE C --- Parameters for finite-difference of Poisson equation CM1X= 1.D0/DX**2 C00X=-2.D0/DX**2 CP1X= 1.D0/DX**2 CM1Z= 1.D0/DZ**2 C00Z=-2.D0/DZ**2 CP1Z= 1.D0/DZ**2 RETURN END C--------------------------------------------------------------------- C----- Mesh generation ..... for Y direction ---------------------- SUBROUTINE MESHY PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) /DY/DY(NY) /DYC/DYC(NY-1) C 壁に垂直な方向には、tanh 関数を使って、壁近傍で格子が密になるように C 配置をしています。(Moin & Kim, J.F.M.118, pp.341-377, 1982と同様) C ALG を大きくすると壁に集中、小さくすると均一な格子になります。 C 現在の設定(NY=64, ALG=0.95)では、壁に接する格子幅は DY+=0.9, C 流路中央では DY+=9（DZ+と同じ）となっています。 C スタガード格子なので、最初の速度点(u,w)は壁から y+=0.45 にあります。 C 0.90 < ALG < 0.98 が妥当と思いますが、陽解法では ALG=0.95 くらいに C 止めた方がよいでしょう。格子を壁に寄せすぎると 流路中央部でまばらに C なってしまいます。また、あまり壁に寄せても結果に著しい改善は期待でき C ません。前述のとおり、DZ+を細かくする方が効果的です。 ALG=0.95D0 AT=DLOG((1.D0+ALG)/(1.D0-ALG))/2.D0 YV(0)=0. DO 10 J=1,NY-1 ETA=AT*(-1.D0+2.D0*DBLE(J)/DBLE(NY)) YV(J)=(DTANH(ETA)/ALG+1.D0)/2.D0 10 CONTINUE YV(NY)=1.D0 DO 15 J=1,NY ETA=AT*(-1.D0+2.D0*(DBLE(J)-0.5D0)/DBLE(NY)) YP(J)=(DTANH(ETA)/ALG+1.D0)/2.D0 15 CONTINUE C ... Outer points (half mesh) YP(0)=2.D0*YV(0)-YP(1) YP(NY1)=2.D0*YV(NY)-YP(NY) DO 20 J=1,NY DY(J)=-YV(J-1)+YV(J) 20 CONTINUE DO 30 J=1,NY-1 DYC(J)=-YP(J)+YP(J+1) 30 CONTINUE RETURN END C----------------------------------------------------------------------- C---- FDM operators at V points using data on P points and INSIDE walls SUBROUTINE MDQ1VN PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) /D1VN/D1VNM(0:NY),D1VNP(0:NY) DO 10 J=1,NY-1 DYV=YP(J+1)-YP(J) D1VNM(J)=-1.D0/DYV D1VNP(J)= 1.D0/DYV 10 CONTINUE C For the Neumann B.C. at the wall (DP/DY=0) for Pressure D1VNM(0)= 0.D0 D1VNP(0)= 0.D0 D1VNM(NY)= 0.D0 D1VNP(NY)= 0.D0 RETURN END C----------------------------------------------------------------------- C---- FDM operators at V points using data on P points and AT walls ---- SUBROUTINE MDQ1VD PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) COMMON /D0V/D0VM(NY-1),D0VP(NY-1) /D1V/D1VM(0:NY),D1VP(0:NY) DO 10 J=1,NY-1 DYV=YP(J+1)-YP(J) D0VM(J)=(YP(J+1)-YV(J))/DYV D0VP(J)=(YV(J)-YP(J))/DYV D1VM(J)=-1.D0/DYV D1VP(J)= 1.D0/DYV 10 CONTINUE Caution! Stencils are shifted at the wall point DYP0=YP(2)-YP(1) D1VM(0)= YP(2)/YP(1)/DYP0 D1VP(0)=-YP(1)/YP(2)/DYP0 DYPN=YP(NY)-YP(NY-1) D1VM(NY)= (1.-YP(NY))/(1.-YP(NY-1))/DYPN D1VP(NY)=-(1.-YP(NY-1))/(1.-YP(NY))/DYPN RETURN END C----------------------------------------------------------------------- C--- FDM operators for velocity at P points using data on V points ---- SUBROUTINE MDQ1P PARAMETER (NY=64,NY1=NY+1,NYH1=NY-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) /DY/DY(NY) COMMON /D0P/D0PM(NY),D0PP(NY) /D1P/D1PM(NY),D1PP(NY) DO 10 J=1,NY D0PM(J)=(YV(J)-YP(J))/DY(J) D0PP(J)=(YP(J)-YV(J-1))/DY(J) D1PM(J)=-1.D0/DY(J) D1PP(J)= 1.D0/DY(J) 10 CONTINUE RETURN END C----------------------------------------------------------------------- C--- FDM operators for viscous terms at V points --------------------- SUBROUTINE MDQ2V PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /D2V/D2VM(NY-1),D2V0(NY-1),D2VP(NY-1) COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D1P/D1PM(NY),D1PP(NY) DO 10 J=1,NY-1 D2VM(J)=D1VM(J)*D1PM(J) D2V0(J)=D1VM(J)*D1PP(J)+D1VP(J)*D1PM(J+1) D2VP(J)=D1VP(J)*D1PP(J+1) 10 CONTINUE RETURN END C----------------------------------------------------------------------- C--- FDM operators for viscous terms at P points --------------------- SUBROUTINE MDQ2PV PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /D2P/D2PM(2:NY),D2P0(NY),D2PP(NY-1) COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D1P/D1PM(NY),D1PP(NY) J=1 D2P0(J)=D1PM(J)*D1VM(J-1)+D1PP(J)*D1VM(J) D2PP(J)=D1PM(J)*D1VP(J-1)+D1PP(J)*D1VP(J) DO 10 J=2,NY-1 D2PM(J)=D1PM(J)*D1VM(J-1) D2P0(J)=D1PM(J)*D1VP(J-1)+D1PP(J)*D1VM(J) D2PP(J)= D1PP(J)*D1VP(J) 10 CONTINUE J=NY D2PM(J)=D1PM(J)*D1VM(J-1)+D1PP(J)*D1VM(J) D2P0(J)=D1PM(J)*D1VP(J-1)+D1PP(J)*D1VP(J) RETURN END C----------------------------------------------------------------------- C--- FDM operators for pressure's Poisson equation ----------------- SUBROUTINE MDQ2PP PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /DPP/DPPM(2:NY),DPP0(NY),DPPP(NY-1) COMMON /D1VN/D1VNM(0:NY),D1VNP(0:NY) /D1P/D1PM(NY),D1PP(NY) DO 10 J=1,NY IF(J.GT. 1) DPPM(J)=D1PM(J)*D1VNM(J-1) DPP0(J)=D1PM(J)*D1VNP(J-1)+D1PP(J)*D1VNM(J) IF(J.LT.NY) DPPP(J)=D1PP(J)*D1VNP(J) 10 CONTINUE RETURN END C ********************************************************************* C * SBR. DTR : DATA READ * C ********************************************************************* SUBROUTINE SBRDTR(CT) PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*3 CT CHARACTER*4 FILE1 COMMON /FILE1/FILE1 /START/ISTART /TSTEP/TSTEP COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /P/P(NZ,NX,NY) OPEN(10,file=FILE1//CT//'u.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,2000) U CLOSE(10) OPEN(10,file=FILE1//CT//'v.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,2000) V CLOSE(10) OPEN(10,file=FILE1//CT//'w.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,2000) W CLOSE(10) OPEN(10,file=FILE1//CT//'p.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,3000) P CLOSE(10) 1000 FORMAT(I10,F12.6) 2000 FORMAT(8F10.6) 3000 FORMAT(8F10.5) WRITE(6,1200) ISTART,TSTEP 1200 FORMAT(/'Data READ',10X,'STEP=',I7,10X,'T=',F10.3) RETURN END C ********************************************************************* C * SBR. CON : NONLINEAR TERM * C ********************************************************************* SUBROUTINE SBRCON CALL SBRNLU CALL SBRNLV CALL SBRNLW RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLU PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /UF/UF(NZ,NX,NY) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W1/CX(NZ,NX,NY) /W2/CY(NZ,NX,0:NY) /W3/CZ(NZ,NX,NY) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D0P/D0PM(NY),D0PP(NY) COMMON /IP/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(NZ) C ... CX:-U^XU_X at P, CZ:-W^XU_Z at P DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ CX(K,I,J)=-TXA*(U(K,IM(I),J)+U(K,I,J))*(-U(K,IM(I),J)+U(K,I,J)) CZ(K,I,J)=-TZA*(W(K,I,J)+W(K,IP(I),J))*(-U(K,I,J)+U(KP(K),I,J)) 10 CONTINUE C ... CY:-V^XU_Y at UV DO 30 J=0,NY IF(J.GT.0.AND.J.LT.NY) THEN DO 31 I=1,NX DO 31 K=1,NZ CY(K,I,J)=-(D1VM(J)*U(K,I,J)+D1VP(J)*U(K,I,J+1)) & *(V(K,I,J)+V(K,IP(I),J))/2.D0 31 CONTINUE ELSE DO 32 I=1,NX DO 32 K=1,NZ CY(K,I,J)=0.D0 32 CONTINUE ENDIF 30 CONTINUE C ... UF:-(U^XU_X)^X-(V^XU_Y)^Y-(W^XU_Z)^Z at U DO 50 J=1,NY DO 50 I=1,NX DO 50 K=1,NZ UF(K,I,J)=(CX(K,I,J)+CX(K,IP(I),J)+CZ(KM(K),I,J)+CZ(K,I,J))/2.D0 & +D0PM(J)*CY(K,I,J-1)+D0PP(J)*CY(K,I,J) 50 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLV PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /VF/VF(NZ,NX,0:NY) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W1/CX(NZ,NX,NY) /W2/CY(NZ,NX,0:NY) /W3/CZ(NZ,NX,NY) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA COMMON /D0V/D0VM(NY-1),D0VP(NY-1) COMMON /D0P/D0PM(NY),D0PP(NY) /D1P/D1PM(NY),D1PP(NY) COMMON /IP/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(NZ) C ... CY:-(V^Y*V_Y)^Y at P DO 22 J=1,NY DO 22 I=1,NX DO 22 K=1,NZ CY(K,I,J)=-(D0PM(J)*V(K,I,J-1)+D0PP(J)*V(K,I,J)) & *(D1PM(J)*V(K,I,J-1)+D1PP(J)*V(K,I,J)) 22 CONTINUE C ... CX:-(U^Y*V_X) at UV, CZ:-(W^Y*V_Z) at VW DO 30 J=1,NY-1 DO 30 I=1,NX DO 30 K=1,NZ CX(K,I,J)=-(D0VM(J)*U(K,I,J)+D0VP(J)*U(K,I,J+1)) & *BXA*(-V(K,I,J)+V(K,IP(I),J)) CZ(K,I,J)=-(D0VM(J)*W(K,I,J)+D0VP(J)*W(K,I,J+1)) & *BZA*(-V(K,I,J)+V(KP(K),I,J)) 30 CONTINUE C ... VF:-(U^Y*V_X)^X-(V^Y*V_Y)^Y-(W^Y*V_Z)^Z at V DO 50 J=1,NY-1 DO 50 I=1,NX DO 50 K=1,NZ VF(K,I,J)=(CX(K,IM(I),J)+CX(K,I,J)+CZ(KM(K),I,J)+CZ(K,I,J))/2.D0 & +D0VM(J)*CY(K,I,J)+D0VP(J)*CY(K,I,J+1) 50 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLW PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /WF/WF(NZ,NX,NY) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W1/CX(NZ,NX,NY) /W2/CY(NZ,NX,0:NY) /W3/CZ(NZ,NX,NY) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D0P/D0PM(NY),D0PP(NY) COMMON /IP/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(NZ) C ... CX:-U^ZW_X at UW, CZ:-W^ZW_Z at P DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ CX(K,I,J)=-TXA*(U(K,I,J)+U(KP(K),I,J))*(-W(K,I,J)+W(K,IP(I),J)) CZ(K,I,J)=-TZA*(W(KM(K),I,J)+W(K,I,J))*(-W(KM(K),I,J)+W(K,I,J)) 10 CONTINUE C ... CY:-V^ZW_Y at VW DO 30 J=0,NY IF(J.GT.0.AND.J.LT.NY) THEN DO 31 I=1,NX DO 31 K=1,NZ CY(K,I,J)=-(D1VM(J)*W(K,I,J)+D1VP(J)*W(K,I,J+1)) & *(V(K,I,J)+V(KP(K),I,J))/2.D0 31 CONTINUE ELSE DO 32 I=1,NX DO 32 K=1,NZ CY(K,I,J)=0.D0 32 CONTINUE ENDIF 30 CONTINUE C ... WF:-(U^ZW_X)^X-(V^ZW_Y)^Y-(W^ZW_Z)^Z at W DO 50 J=1,NY DO 50 I=1,NX DO 50 K=1,NZ WF(K,I,J)=(CX(K,IM(I),J)+CX(K,I,J)+CZ(K,I,J)+CZ(KP(K),I,J))/2.D0 & +D0PM(J)*CY(K,I,J-1)+D0PP(J)*CY(K,I,J) 50 CONTINUE RETURN END C ********************************************************************* C * SBR. VIS : VISCOUS TERM * C ********************************************************************* SUBROUTINE SBRVIS PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) /RET/RET COMMON /UF/UF(NZ,NX,NY) /VF/VF(NZ,NX,0:NY) /WF/WF(NZ,NX,NY) COMMON /D2V/D2VM(NY-1),D2V0(NY-1),D2VP(NY-1) COMMON /D2P/D2PM(2:NY),D2P0(NY),D2PP(NY-1) COMMON /IP/IP(NX) /IM/IM(NX) /CMPX/CM1X,C00X,CP1X COMMON /KP/KP(NZ) /KM/KM(NZ) /CMPZ/CM1Z,C00Z,CP1Z DO 10 J=1,NY C00=C00X+C00Z+D2P0(J) IF(J.EQ.1) THEN DO 21 I=1,NX DO 21 K=1,NZ UF(K,I,J)=UF(K,I,J) &+(CM1X*U(K,IM(I),J)+CM1Z*U(KM(K),I,J)+C00*U(K,I,J) & +CP1Z*U(KP(K),I,J)+CP1X*U(K,IP(I),J)+D2PP(J)*U(K,I,J+1))/RET WF(K,I,J)=WF(K,I,J) &+(CM1X*W(K,IM(I),J)+CM1Z*W(KM(K),I,J)+C00*W(K,I,J) & +CP1Z*W(KP(K),I,J)+CP1X*W(K,IP(I),J)+D2PP(J)*W(K,I,J+1))/RET 21 CONTINUE ELSE IF(J.LT.NY) THEN DO 20 I=1,NX DO 20 K=1,NZ UF(K,I,J)=UF(K,I,J) &+(D2PM(J)*U(K,I,J-1)+CM1X*U(K,IM(I),J)+CM1Z*U(KM(K),I,J) & +C00*U(K,I,J) & +CP1Z*U(KP(K),I,J)+CP1X*U(K,IP(I),J)+D2PP(J)*U(K,I,J+1))/RET WF(K,I,J)=WF(K,I,J) &+(D2PM(J)*W(K,I,J-1)+CM1X*W(K,IM(I),J)+CM1Z*W(KM(K),I,J) & +C00*W(K,I,J) & +CP1Z*W(KP(K),I,J)+CP1X*W(K,IP(I),J)+D2PP(J)*W(K,I,J+1))/RET 20 CONTINUE ELSE DO 22 I=1,NX DO 22 K=1,NZ UF(K,I,J)=UF(K,I,J) &+(D2PM(J)*U(K,I,J-1)+CM1X*U(K,IM(I),J)+CM1Z*U(KM(K),I,J) & +C00*U(K,I,J)+CP1Z*U(KP(K),I,J)+CP1X*U(K,IP(I),J))/RET WF(K,I,J)=WF(K,I,J) &+(D2PM(J)*W(K,I,J-1)+CM1X*W(K,IM(I),J)+CM1Z*W(KM(K),I,J) & +C00*W(K,I,J)+CP1Z*W(KP(K),I,J)+CP1X*W(K,IP(I),J))/RET 22 CONTINUE ENDIF ENDIF 10 CONTINUE DO 30 J=1,NY-1 C00=C00X+C00Z+D2V0(J) DO 30 I=1,NX DO 30 K=1,NZ VF(K,I,J)=VF(K,I,J) &+(D2VM(J)*V(K,I,J-1)+CM1X*V(K,IM(I),J)+CM1Z*V(KM(K),I,J) & +C00*V(K,I,J) & +CP1Z*V(KP(K),I,J)+CP1X*V(K,IP(I),J)+D2VP(J)*V(K,I,J+1) & )/RET 30 CONTINUE RETURN END C ********************************************************************* C * SBR. PRE : PREDICTION STEP * C ********************************************************************* SUBROUTINE SBRPRE(AB,BB) PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /DT/DT /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /UF/UF(NZ,NX,NY) /VF/VF(NZ,NX,0:NY) /WF/WF(NZ,NX,NY) COMMON /UB/UB(NZ,NX,NY) /VB/VB(NZ,NX,NY-1) /WB/WB(NZ,NX,NY) DAB=AB*DT DBB=BB*DT DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ UF(K,I,J)=UF(K,I,J)+2.D0 U(K,I,J)=U(K,I,J)+DAB*UF(K,I,J)+DBB*UB(K,I,J) W(K,I,J)=W(K,I,J)+DAB*WF(K,I,J)+DBB*WB(K,I,J) UB(K,I,J)=UF(K,I,J) WB(K,I,J)=WF(K,I,J) 10 CONTINUE DO 20 J=1,NY-1 DO 20 I=1,NX DO 20 K=1,NZ V(K,I,J)=V(K,I,J)+DAB*VF(K,I,J)+DBB*VB(K,I,J) VB(K,I,J)=VF(K,I,J) 20 CONTINUE RETURN END C ********************************************************************* C * SBR. RHP : R.H.S. OF POISSON EQ. * C ********************************************************************* C --- R.H.S. FOR SCALER POTANTIAL ------------------------------------- SUBROUTINE SBRRHP PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W3/Q(NZ,NX,NY) /DT/DT /D/DX,DZ COMMON /D1P/D1PM(NY),D1PP(NY) /IM/IM(NX) /KM/KM(NZ) DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ Q(K,I,J)=((-U(K,IM(I),J)+U(K,I,J))/DX & +D1PM(J)*V(K,I,J-1)+D1PP(J)*V(K,I,J) & +(-W(KM(K),I,J)+W(K,I,J))/DZ )/DT 10 CONTINUE RETURN END C ********************************************************************* C * SBR. SOR : S.O.R. SCHEME FOR POISSON EQ. * C ********************************************************************* SUBROUTINE SBRSOR(ISOR) PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /DPP/DPPM(2:NY),DPP0(NY),DPPP(NY-1) COMMON /P/P(NZ,NX,NY) /W3/Q(NZ,NX,NY) /ERRP/POIERR /ITRP/ITRP COMMON /IP/IP(NX) /IM/IM(NX) /CMPX/CM1X,C00X,CP1X COMMON /KP/KP(NZ) /KM/KM(NZ) /CMPZ/CM1Z,C00Z,CP1Z ITRP=0 SUMS=0.D0 DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ SUMS=SUMS+Q(K,I,J)**2 10 CONTINUE 100 ITRP=ITRP+1 SUMR=0. DO 20 J=1,NY C00=C00X+C00Z+DPP0(J) DPC=-1.5D0/C00 IF(J.EQ.1) THEN DO 31 I=1,NX DO 31 K=1,NZ RESI=CM1X*P(K,IM(I),J)+CM1Z*P(KM(K),I,J)+C00*P(K,I,J) & +CP1Z*P(KP(K),I,J)+CP1X*P(K,IP(I),J)+DPPP(J)*P(K,I,J+1) & -Q(K,I,J) P(K,I,J)=P(K,I,J)+RESI*DPC SUMR=SUMR+RESI**2 31 CONTINUE ELSE IF(J.LT.NY) THEN DO 30 I=1,NX DO 30 K=1,NZ RESI=DPPM(J)*P(K,I,J-1)+CM1X*P(K,IM(I),J)+CM1Z*P(KM(K),I,J) & +C00*P(K,I,J) & +CP1Z*P(KP(K),I,J)+CP1X*P(K,IP(I),J)+DPPP(J)*P(K,I,J+1) & -Q(K,I,J) P(K,I,J)=P(K,I,J)+RESI*DPC SUMR=SUMR+RESI**2 30 CONTINUE ELSE DO 32 I=1,NX DO 32 K=1,NZ RESI=DPPM(J)*P(K,I,J-1)+CM1X*P(K,IM(I),J)+CM1Z*P(KM(K),I,J) & +C00*P(K,I,J)+CP1Z*P(KP(K),I,J)+CP1X*P(K,IP(I),J) & -Q(K,I,J) P(K,I,J)=P(K,I,J)+RESI*DPC SUMR=SUMR+RESI**2 32 CONTINUE ENDIF ENDIF 20 CONTINUE POIERR=DSQRT(SUMR/SUMS) C WRITE(6,*) 'ITRP=',ITRP,' ERR=',POIERR IF(ITRP.GE.ISOR.OR.POIERR.LT.1.D-5) GO TO 200 GO TO 100 200 CONTINUE RETURN END C ********************************************************************* C * SBR. COR : CORRECTION STEP * C ********************************************************************* SUBROUTINE SBRCOR PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /P/P(NZ,NX,NY) /DT/DT /D/DX,DZ /DY/DY(NY) COMMON /IP/IP(NX) /KP/KP(NZ) /D1VN/D1VNM(0:NY),D1VNP(0:NY) TX=DT/DX TZ=DT/DZ DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ U(K,I,J)=U(K,I,J)-TX*(-P(K,I,J)+P(K,IP(I),J)) W(K,I,J)=W(K,I,J)-TZ*(-P(K,I,J)+P(KP(K),I,J)) 10 CONTINUE DO 20 J=1,NY-1 DO 20 I=1,NX DO 20 K=1,NZ V(K,I,J)=V(K,I,J)-DT*(D1VNM(J)*P(K,I,J)+D1VNP(J)*P(K,I,J+1)) 20 CONTINUE RETURN END C ********************************************************************* C * SBR. UMR : MEAN & RMS VALUES * C ********************************************************************* SUBROUTINE SBRUMR PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /UM/UM(NY),VM(0:NY),WM(NY) /S12/S12L(0:NY),S12V(0:NY) COMMON /IP/IP(NX) /IM/IM(NX) /D0V/D0VM(NY-1),D0VP(NY-1) COMMON /BVU/BVU2(NY),BVU3(NY),BVU4(NY) COMMON /BVV/BVV2(0:NY),BVV3(0:NY),BVV4(0:NY) COMMON /BVW/BVW2(NY),BVW3(NY),BVW4(NY) ARXZ=1./DBLE(NX*NZ) DO 20 J=1,NY SUMU1=0.D0 SUMU2=0.D0 SUMW1=0.D0 SUMW2=0.D0 DO 10 I=1,NX DO 10 K=1,NZ SUMU1=SUMU1+U(K,I,J) SUMU2=SUMU2+U(K,I,J)**2 SUMW1=SUMW1+W(K,I,J) SUMW2=SUMW2+W(K,I,J)**2 10 CONTINUE UM(J)=ARXZ*SUMU1 BVU2(J)=ARXZ*SUMU2 WM(J)=ARXZ*SUMW1 BVW2(J)=ARXZ*SUMW2 20 CONTINUE DO 50 J=1,NY-1 SUMV1=0.D0 SUMV2=0.D0 SUMUV=0.D0 DO 40 I=1,NX DO 40 K=1,NZ SUMV1=SUMV1+V(K,I,J) SUMV2=SUMV2+V(K,I,J)**2 SUMUV=SUMUV+(D0VM(J)*U(K,I,J)+D0VP(J)*U(K,I,J+1)) & *(V(K,I,J)+V(K,IP(I),J)) 40 CONTINUE VM(J)=ARXZ*SUMV1 BVV2(J)=ARXZ*SUMV2 S12L(J)=ARXZ*SUMUV/2.D0 50 CONTINUE RETURN END C ********************************************************************* C * SBR. STT : TURBULENCE STATISTICS * C ********************************************************************* C --- Mean Velocity and Turbulence Intensity---------------------------- SUBROUTINE SBRST1 PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /UM/UM(NY),VM(0:NY),WM(NY) /PM/PM(NY) COMMON /UR/UR(NY),VR(0:NY),WR(NY) /PR/PR(NY) COMMON /ENE/ENE /UME/UME,UMX /RMS/URMSX,URMSC,VRMSC,WRMSC COMMON /DY/DY(NY) /DYC/DYC(NY-1) COMMON /BVU/BVU2(NY),BVU3(NY),BVU4(NY) COMMON /BVV/BVV2(0:NY),BVV3(0:NY),BVV4(0:NY) COMMON /BVW/BVW2(NY),BVW3(NY),BVW4(NY) COMMON /BVP/BVP2(NY) ENE=0.D0 UME=0.D0 UMX=0.D0 URMSX=0.D0 DO 10 J=1,NY UR(J)=DSQRT(BVU2(J)-UM(J)**2) WR(J)=DSQRT(BVW2(J)-WM(J)**2) PR(J)=DSQRT(BVP2(J)-PM(J)**2) ENE=ENE+DY(J)*(UR(J)**2+WR(J)**2) UME=UME+DY(J)*UM(J) UMX=DMAX1(UMX,UM(J)) URMSX=DMAX1(URMSX,UR(J)) 10 CONTINUE DO 20 J=1,NY-1 VR(J)=DSQRT(BVV2(J)-VM(J)**2) ENE=ENE+DYC(J)*VR(J)**2 20 CONTINUE URMSC=(UR(NY/2)+UR(NY/2+1))/2.D0 VRMSC=VR(NY/2) WRMSC=(WR(NY/2)+WR(NY/2+1))/2.D0 ENE=ENE/2.D0 RETURN END C ********************************************************************* C * SBR. CHK : CHECK of DIVERGENCE and COURANT-NUMBER * C ********************************************************************* SUBROUTINE SBRCHK PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /D/DX,DZ /DY/DY(NY) /DT/DT /DIV/DIVX /COU/COMX COMMON /D0P/D0PM(NY),D0PP(NY) /D1P/D1PM(NY),D1PP(NY) COMMON /UR/UR(NY),VR(0:NY),WR(NY) COMMON /IM/IM(NX) /KM/KM(NZ) DIVX=0.D0 COMX=0.D0 DO 40 J=1,NY DNOMAL=(DX*DY(J)*DZ)**(1.D0/3.D0)/UR(J) DO 40 I=1,NX DO 40 K=1,NZ DIV=(-U(K,IM(I),J)+U(K,I,J))/DX & +D1PM(J)*V(K,I,J-1)+D1PP(J)*V(K,I,J) & +(-W(KM(K),I,J)+W(K,I,J))/DZ UCP=(U(K,IM(I),J)+U(K,I,J))/2.D0 WCP=(W(KM(K),I,J)+W(K,I,J))/2.D0 VCP=D0PM(J)*V(K,I,J-1)+D0PP(J)*V(K,I,J) DIVX=DMAX1(DIVX,DNOMAL*DIV) COU=DT*(DABS(UCP)/DX+DABS(VCP)/DY(J)+DABS(WCP)/DZ) COMX=DMAX1(COMX,COU) 40 CONTINUE RETURN END C ********************************************************************* C * SBR. DTS : DATA SAVE * C ********************************************************************* SUBROUTINE SBRDTS(CT) PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*3 CT CHARACTER*4 FILE1 COMMON /FILE1/FILE1 /STEP/ISTEP /TSTEP/TSTEP COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /P/P(NZ,NX,NY) OPEN(10,file=FILE1//CT//'u.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,2000) U CLOSE(10) OPEN(10,file=FILE1//CT//'v.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,2000) V CLOSE(10) OPEN(10,file=FILE1//CT//'w.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,2000) W CLOSE(10) OPEN(10,file=FILE1//CT//'p.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,3000) P CLOSE(10) 1000 FORMAT(I10,F12.6) 2000 FORMAT(8F10.6) 3000 FORMAT(8F10.5) WRITE(6,1200) ISTEP,TSTEP 1200 FORMAT(/'Data SAVE',10X,'STEP=',I7,10X,'T=',F10.3/) RETURN END