11-6-3ProgramreadingCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCMODULEUSERL***************************************************************INTEGER*4IJREAL*8 GBR,DPDZ,PR,AMU,FRE,WBAR,TMENDMODULECCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCSUBROUTINEUSER*************水****************************水*********************USE START LUSE USER LIMPLICITNONEC*****************************************************************CPROBLEMSEVEN-cComplexfullydevelopedlaminarfluidflowandheattransferinachorizontal squareductC*************************************************************11/98
11/98 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MODULE USER_L C*************************************************************** INTEGER*4 I,J REAL*8 GBR, DPDZ, PR, AMU, FRE, WBAR, TM END MODULE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE USER C**************************************************************** USE START_L USE USER_L IMPLICIT NONE C***************************************************************** C-PROBLEM SEVEN- - C Complex fully developed laminar fluid flow and heat transfer in a C horizontal square duct C**************************************************************** 11-6-3 Program reading
ENTRYGRIDTITLE(1)-VEL U:TITLE(2)=VEL VTITLE(3)=°STR FN:TITLE(4)=.TEMP:wistreated asfifthvariable!TITLE(5)=W/WBAR.TITLE(11)-PRESSURERELAX(1)=0.8! Notfor w,Withknown u, y,the w eq is linearRELAX(2)=0.8LSOLVE(1)=TRUE! u, y, p, Tare solved firstLSOLVE(4)-.TRUELPRINT(1)=.TRUE!In SIMPLERcode.whenthe1st variableis setLPRINT(2)=.TRUEto be solved, the 2nd and 3rd ones (v and pe) areLPRINT(3)=.TRUELPRINT(4)=.TRUEautomatically regarded as variablestobe solvedLPRINT(5)=.TRUE.LPRINT(11)-TRUE.LAST=-25XL=1.!ComputationfortheentireregionYL=1.L1=7M1=7CALLUGRIDRETURN12/98
12/98 ENTRY GRID TITLE(1)=‘ .VEL U.' TITLE(2)=‘ .VEL V.' TITLE(3)=‘ .STR FN.' TITLE(4)=‘ . TEMP .' TITLE(5)=‘ .W/WBAR.' TITLE(11)='PRESSURE' RELAX(1)=0.8 RELAX(2)=0.8 LSOLVE(1)=.TRUE. LSOLVE(4)=.TRUE. LPRINT(1)=.TRUE. LPRINT(2)=.TRUE. LPRINT(3)=.TRUE. LPRINT(4)=.TRUE. LPRINT(5)=.TRUE. LPRINT(11)=.TRUE. LAST=25 XL=1. YL=1. L1=7 M1=7 CALL UGRID RETURN ！Not for w; With known u, v, the w eq is linear. ！Computation for the entire region ! u, v, p, T are solved first ！In SIMPLER code, when the 1st variable is set to be solved, the 2nd and 3rd ones (v and pc ) are automatically regarded as variables to be solved. ！w is treated as fifth variable!
ENTRYSTARTpgβGBR=1.E4拉DPDZ--3000DO 100 J=1,M1DO 101 I=1,L1免热U(I,J)=0!Initial temperature and someV(IL,J)=0.绝热T(I,J)=0.boundary conditionsT(L1,J)=1F(I,J,5)=100.! Initial field for axial velocity wIF (I==1.OR.I==L1) F(I,J,5)=0.IBoundary cond. of w at four walls w=0IF (J==1.OR.J==M1) F(I,J,5)=0101ENDDO100ENDDO! CPCON-1, default valuePR=0.7!Pr=uc,/2;=uc,/PrAMU=1AMUP=AMU*CPCON/PR！afor temperatureRETURNBossinesq assumption: It ignores density differencesENTRYDENSEexcept where they appear in terms pgRETURN13/98
13/98 ENTRY START GBR=1.E4 DPDZ=-3000. DO 100 J=1,M1 DO 101 I=1,L1 U(I,J)=0. V(I,J)=0. T(I,J)=0. T(L1,J)=1. F(I,J,5)=100. IF (I==1.OR.I==L1) F(I,J,5)=0. IF (J==1.OR.J==M1) F(I,J,5)=0. 101 ENDDO 100 ENDDO PR=0.7 AMU=1. AMUP=AMU*CPCON/PR RETURN ENTRY DENSE RETURN ！ g !Boundary cond. of w at four walls w=0 Bossinesq assumption: It ignores density differences except where they appear in terms g !Initial temperature and some boundary conditions ! Initial field for axial velocity w ! for temperature Pr / ; / Pr p p ！ = = c c ! CPCON=1, default value
ENTRYBOUNDFRE-0.!wisnotsolvedwhenITER<20IF(ITER<20)RETURN!Switchofthesolvedvariables,onlyexecutedonceIF(NOTLSOLVE(5))THENThedefaultvalueofLSOLVEisFALSE.WhenLSOLVE(1)=.FALSEITER=-20,.NOTLSOLVE(5)isTRUE:WhenLSOLVE(5)-.TRUEITER>=21,itis.FALSEENDIF!OnceSthvariableissolved,thesetwoWBAR=0statements arenotneededtoexecutedanymoreDO302J=2.M2DO 303 I=2.L2WBAR=WBAR+F(IJ.5)*XCV()*YCV(J)!Forcomputingaveragevelocity303ENDDO!Computing(fRe)accordingtodefinitionShown in the next page302ENDDOFRE=-DPDZ*2.*4.*(XL*YL)**3/(XL+YX)**2/(WBAR*AMU)RETURN!WBAR=WBAR/(XL*YL)BecausebothXLandYL=l,thiscalculationisignored!14/98
14/98 ENTRY BOUND FRE=0. IF(ITER<20) RETURN IF(.NOT.LSOLVE(5)) THEN LSOLVE(1)=.FALSE. LSOLVE(5)=.TRUE. ENDIF WBAR=0. DO 302 J=2,M2 DO 303 I=2,L2 WBAR=WBAR+F(I,J,5)*XCV(I)*YCV(J) 303 ENDDO 302 ENDDO FRE=-DPDZ*2.*4.*(XL*YL)**3/(XL+YL)**2/(WBAR*AMU) RETURN ！w is not solved when ITER<20 ！Once 5th variable is solved, these two statements are not needed to executed any more. ！Switch of the solved variables, only executed once. The default value of LSOLVE is FALSE. When ITER=20, .NOT.LSOLVE(5) is TRUE; When ITER>=21, it is .FALSE. ！Computing (fRe) according to definition; Shown in the next page. ! For computing average velocity ! WBAR=WBAR/(XL*YL) Because both XL and YL =1, this calculation is ignored!
FRE=-DPDZ*2.*4.*(XL*YL)**3/(XL+YL)**2/(WBAR*AMUf Re=-[(dp / da)D, / /pwaj Pw.D,n-2dp / dzf Re = -2[(dp / dz)D, / wmn] =n(Zwi,jA,A-2dp / dznZwi,jAA,-2dp / dz4 XL*YL·XL*YL2(XL + YL)nZwijAij,j-2dp / dz 4(XL* YL)3nZwi,jAAi,j(XL+YL)15/98
15/98 2 , , 2 / 4( ) i j i j dp dz A A w A P − = • 2 , , 2 / 4 * ( ) * 2( ) i j i j dp dz XL YL XL YL w A XL YL − = • • + 3 2 , , 2 / 4( * ) i j i j ( ) dp dz XL YL w A XL YL − = • + FRE=-DPDZ*2.*4.*(XL*YL)**3/(XL+YL)**2/(WBAR*AMU) 2 , , 2 / 4( ) ( ) i j i j dp dz A w A A P = • − 1 2 Re [( / ) / ] 2 m h h m w D f dp dz D w = − 2 Re 2[( / ) / ] h m f dp dz D w = −