**火火大火大火火火大火EQUIVALENCE(F(1,1,1),U(1,1)).(F(1,1,2),V(1,1)).(F(1,1,3),PC(1,1))&, (F(1,1,4),T(1,1))EQUIVALENCE(F(1,1,11),P(1,1).(F(1,1,12),RHO(1,1).(F(1,1,13)&,GAM(1,1),(F(1,1,14),CP(1,1)EQUIVALENCE(COF(1,1,1),CON(1,1)),(COF(1,1,2),AIP(1,1)&(COF(1,1,3),AIM(1,1)),(COF(1,1,4),AJP(1,1)),&(COF(1,1,5),AJM(1,1)),.(COF(1,1,6),AP(1,1))REAL*8,DIMENSION(NI)::TH,THU,THDIF,THCV,THCVSREAL*8THLEQUIVALENCE(X,TH),(XU,THU),(XDIF,THDIF),(XCV,THCV)&(XCVS,THCVS),(XL,THL)DATALSTOP,LSOLVE,LPRINT/.FALSE.,NFX4*.FALSE.,NFX4*.FALSE./DATALBLK/NFX4*.TRUE./DefaultDATAMODE.LAST.TIME.ITER/1,5,0..0/value!!DATARELAX.NTIMES/NFX4*1..NFX4*1/DATADT,IPREF,JPREF,RHOCON,CPCON/1.E+30,1,1,1.,1./ENDMODULEMODULEmodulename(1)Packagingdata（封装数据）；Module name is composed of two parts,(2）Initializingdata（数据初始化）；with a hyphen(-) at bottom in between.(3)Declaringtypeofdata（声明数据类型）11/79ENDMODULE
C****************************************** EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) &, (F(1,1,4),T(1,1)) EQUIVALENCE(F(1,1,11),P(1,1)),(F(1,1,12),RHO(1,1)),(F(1,1,13) &,GAM(1,1),(F(1,1,14),CP(1,1)) EQUIVALENCE(COF(1,1,1),CON(1,1)),(COF(1,1,2),AIP(1,1)), &(COF(1,1,3),AIM(1,1)),(COF(1,1,4),AJP(1,1)), &(COF(1,1,5),AJM(1,1)),(COF(1,1,6),AP(1,1)) REAL*8,DIMENSION(NI)::TH,THU,THDIF,THCV,THCVS REAL*8 THL EQUIVALENCE(X,TH),(XU,THU),(XDIF,THDIF),(XCV,THCV), &(XCVS,THCVS),(XL,THL) DATA LSTOP,LSOLVE,LPRINT/.FALSE.,NFX4*.FALSE., NFX4*.FALSE./ DATA LBLK/NFX4*.TRUE./ DATA MODE,LAST,TIME,ITER/1,5,0.,0/ DATA RELAX,NTIMES/NFX4*1.,NFX4*1/ DATA DT,IPREF,JPREF,RHOCON,CPCON/1.E+30, 1,1,1.,1./ END MODULE 11/79 Default value!!
Some explains to this mostimportant moduleREAL*8.DIMENSION(NINJ.NFX4)::FReal variable 3-D array, array title F, F(NI,NJ,NFX4);Variable number in three coordinates are NI,NJ and NFX4respectively,REAL*8 SMAX,SSUMReal variable of SMAX and SSUM, with length of eight digits:INTEGER*4NE,NP,NRHOInteger variable of NF,NP,NRHO, with length of four digits;EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1))Making the 1st variable of the 3D arrayF identical to the 2Darray U; the same for (F(1,1,2), V(1,1)12/79
REAL*8 SMAX,SSUM Real variable of SMAX and SSUM, with length of eight digits; INTEGER*4 NF,NP,NRHO Integer variable of NF,NP,NRHO, with length of four digits; REAL*8,DIMENSION(NI,NJ,NFX4)::F Real variable 3-D array, array title F，F(NI,NJ,NFX4); Variable number in three coordinates are NI,NJ and NFX4 respectively; EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)) Making the 1st variable of the 3D array F identical to the 2D array U; the same for (F(1,1,2), V(1,1))) Some explains to this most important module 12/79
1062.2PROGRAMMAIN(2)-Explained**求***水***水求水********水*******水水***水*水**水*水******求***水求**in detailMAIN**水*水**水**水**水水*水水*水水水*水***水水术*水水水*水**水水水*水PROGRAMMAIN!SetupentireflowchartCallingUSESTART LSharethe variables definedin theMODULEModuleIMPLICITNONE****************************************************OPEN(8,FILE='RESULT.txt')!ResultfileforoutputOnlyCALLGRID!Grid generation (setupinterface positions)Executed!Set up1-D array of geometric parametersCALLSETUP1CALLSTART!Set up initial fieldonceDO WHILE (NOT.LSTOP)! IfLSTOPis .F., then do loopCALLDENSE!Setupfluid densityCALLBOUND!Setup boundary conditionExecutedCALLOUTPUT!Print out representativeresultsManytimesCALLSETUP2!Keymodule:setcoefficientsand solveABEqs.ENDDOCALLOUTPUT!Print out some results!Simulationcompleted closefileRESULT.TXTinChannel8CLOSE(8)STOP!TerminatecomputationEND!Endofmainprogram13/79CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C*********************************************************** C-MAIN - C*********************************************************** PROGRAM MAIN USE START_L IMPLICIT NONE C********************************************************** OPEN(8,FILE='RESULT.txt') CALL GRID CALL SETUP1 CALL START DO WHILE (.NOT.LSTOP） CALL DENSE CALL BOUND CALL OUTPUT CALL SETUP2 ENDDO CALL OUTPUT CLOSE(8) STOP END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !Grid generation (setup interface positions) !Set up 1-D array of geometric parameters !Set up initial field ！If LSTOP is .F., then do loop !Set up fluid density !Set up boundary condition !Print out representative results !Key module: set coefficients and solve ABEqs. !Simulation completed close file RESULT.TXT in Channel 8 !Terminate computation ! Result file for output !Print out some results Share the variables defined in the MODULE 10.6.2.2 PROGRAM MAIN 13/79 !End of main program ! Set up entire flow chart Only Executed once Executed Many times Calling Module (2)-Explained in detail
For ABEqs, we determine coefficients using Subroutine DIFLOWFor all other scalar variables (including T, etc)(p'0) + div(p'oU) = div(T,grad)+ S;atapp,=ape +awdw+ard +asds+baw = DwA(|Pw+[Fw, 0ag = DA(PeD+-F,0ae/De:SchemeCentral differenceUpwinddifference1-0.5 Pae1+ -PAe, 0DefinitionHybridPower-lawExponentialFDPAeDePe,1--p0, (1 --0.1PA)5 + 0, -Pseexp(Pe)-1214/75
14/75 ( ) ,0 E e e e a D A P F = + − ( ) ,0 W w w w a D A P F = + P P E E W W N N S S a a a a a b = + + + + * * * ( ) div U div grad S ( ) ( ) t + = + For all other scalar variables (including T, etc) For ABEqs, we determine coefficients using Subroutine DIFLOW e e e F P D = aE / De :
10.6.2.3SUBROUTINEDIELOW(3)-ExplainedCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCin detailSUBROUTINE DIFLOW ! Determine D.A(P) of power law scheme i!TheinputvariablesareDIFFand FLOWUSESTARTL!Share thevariables defined in theMODULE START LIMPLICITNONEREAL*8TEMP!temporalvariable火欢ACOF=DIFF ! D.A(P)= D (ACOF represents D. A(|P)!1st return forIF(FLOW==o.)RETURN!Noflow,onlydiffusiondiffusion caseTEMP=DIFF-ABS(FLOW)*0.1 ! D-0.1|F= D(1-0.1|P DIP/>100ACOF=0. A(P D=max[0,(1-0.1|P D](1-0.1|Pbs5/P/<10!2nd return forIF(TEMP<=0.) RETURN!P|>10[P/>10TEMP-TEMP/DIFF 1-0.1|PACOF=DIFF*TEMP**5 ! D.(1-O'.1|PD' =D.A(PD!3rd return forRETURN[P/<10ae = D,A(|PD+[0,-F]!In SETUP2:$END15/79
USE START_L IMPLICIT NONE REAL*8 TEMP C********************************************************* ACOF=DIFF ！ IF(FLOW== 0.) RETURN ! TEMP=DIFF-ABS(FLOW)*0.1 ！ ACOF=0. IF(TEMP <= 0.) RETURN！ TEMP=TEMP/DIFF ！ ACOF=DIFF*TEMP**5 ！ RETURN END D F D P 0.1 (1 0.1 ) − = − 5 5 0 max[0,(1 0.1 ) ] ( 10 ( ) 1 0.1 ) 10 P A P P P P = − − ！ P 10 5 D P D A P (1 0.1 ) ( ) − = 1 0.1 − P D A P D ( ) = CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE DIFLOW ! Determine ( ) [0, ] E e e !In SETUP2: a = + − D A P F No flow, only diffusion of power law scheme ! Share the variables defined in the MODULE START_L 10.6.2.3 SUBROUTINE DIFLOW 15/79 ! temporal variable ! 1st return for diffusion case ! 2nd return for P 10 ! 3rd return for P 10 D A P( ) (3)-Explained in detail （ACOF represents ） !The input variables are DIFF and FLOW D A P( )