《数值传热学》研究生课程教学资源（教案）流动传热问题数值仿真多层次案例

ovner =n+n, nnox(7ovXDIF()1aOu(n.OxOxXDIF()(GAM(i,)(G+1,)-u(,)xcv(i)(i-1,)()YCVO)u(i,j)-u(i-1,j)GAM(i-1 )xcv(i-1)XCVG-DXCV)Theaboveterm is takenas S,ofu-equation!图1-3动量方程源项的处理方法4.程序解读与说明对于教学程序的运用，主程序不需要改变，我们只需要对用户子程序进行改写。用户子程序包含六个模块，MODULEGRID，MODULESTART，MODULEDENSE，MODULEBOUND，MODULEOUTPUT，MODULEGAMSOR。针对这个端流流动换热问题，我们需要对这六个模块进行定义。(I)MODULEGRID在这个模块中，我们定义需要输出的二维场，需要求解的变量及其松弛因子，以及区域的大小和网格节点数，最后调用UGRID函数来生成网格。在该程序中可以求解10+4个变量，其中第1~4个变量默认为u,v,p,T,而最后四个变量分别为p,p,「,Cp。由于该问题还需要求解k，6方程，我们分别将其定义为第5个变量和第6个变量，求解了k，6方程后，用其来确定流粘性系数，并将其作为第7个变量。各变量的松弛因子如表1所示。该例子中，我们需要观察速度u,v，流函数，温度，k,6，，压力、速度这些二维场信息，因而将相应的TITLE,LSOVE，及LPRINT都设置为.TRUE.，如图1-4所示。RELAXO-0.SRELAX(2)-0.8RELAX(5)=0.6RELAX(6)-0.6!NF=7 for turbulent viscosity RELAX(7)-0.6LSOLVE-TRUELSOLVE(S)=TRUELSOLVE(O-.TRUELPRINT(I)-.TRUE.LPRINT(2)-.TRUEAll logical values forLPRINT(3)-.IRUE.solving and printingLPRINT(4)-.IRUE.LPRINT(5)-.TRUELPRINT()-.TRUELPRINT()-.TRUELPRINT(II)-.TRUELAST-100Regarding hasXL=1.the 7u element ofYL-4.F(LJ,NF)L1-7M=9!C,intheTexpressionfortemperatureCPCON=1000CALLUGRIDI,==nc,/rRETURN图1-4MODULEGRID定义

图 1-3 动量方程源项的处理方法 4. 程序解读与说明 对于教学程序的运用，主程序不需要改变，我们只需要对用户子程序进行改写。用户 子程序包含六个模块，MODULE GRID, MODULE START, MODULE DENSE, MODULE BOUND, MODULE OUTPUT, MODULE GAMSOR。针对这个湍流流动换热问题，我们需 要对这六个模块进行定义。 (1) MODULE GRID 在这个模块中，我们定义需要输出的二维场，需要求解的变量及其松弛因子，以及区 域的大小和网格节点数，最后调用 UGRID 函数来生成网格。在该程序中可以求解 10+4 个 变量，其中第 1~4 个变量默认为 u, v, p’, T, 而最后四个变量分别为 p, , , Cp。由于该问题 还需要求解 k, 方程，我们分别将其定义为第 5 个变量和第 6 个变量，求解了 k, 方程后， 用其来确定湍流粘性系数t，并将其作为第 7个变量。各变量的松弛因子如表 1 所示。该例 子中，我们需要观察速度 u, v, 流函数，温度，k, ，t，压力、速度这些二维场信息，因而 将相应的 TITLE, LSOVE, 及 LPRINT 都设置为.TRUE., 如图 1-4 所示。 图 1-4 MODULE GRID 定义

(2)MODULESTART在这个模块中，我们定义一些初始值和不变的边界条件，由于入口边界条件不变，因为壁面上的速度和温度已知，都可以在这个模块中定义。如图1-5所示，我们分别对速度初场、温度初场、k初场，初场及其相应的边界条件进行处理。ENTRYSTARTDO 100 J=1,M11%of inletkinetic energy,DO 101I=1,L1U(L,J-0.initial value,alsoB.C.for inletV(LJ)-10V(1,J)-0.V(I,2)-10n.:detemied fromIF(L.GI.4) V(I,2)-100.T(,J)=100.pV(2Lm)=100Re.T(1,J=0.nIF(L.GT.4) T(I,1)=400.AKE(LJ)-0.005*V(L.2)**21×100×1.0DIS(LJ)-0.1*AKE(LJ)**2100=.n=1.0101ENDDOn100ENDDOInitial value,also8=Cupkz/n,=0.09x1xk20.1kB.C.for inlet!图1-5MODULESTART定义(3）MODULEDENSE在这个模块中，需要对流体的密度进行定义。由于密度为常数，所以这个模块，我们无须进行定义，按主程序中的默认值给定。虽然该模块为空块，但必须保留，不能删除，否则会引起程序错误。这是因为主程序中调用了这个模块。(4）MODULEBOUND在这个模块里面，我们对该物理问题的边界条件进行设置。我们拟打算用出口速度提升的方法，来防止送代过程中出口法向流速出现负值对计算区域内的影响。具体实施过程如图1-6所示。ENTRYBOUNDIF(ITER-- 0) THENFLOWIN-DO 310 1-2,L2FLOWIN=FLOwIN+RHO(,)*V(,2)*xCV()!Flowrateat inlet310ENDDOFLOWINELSEFACTORFL=0.TAFL=0.Z[(V +V*RHOu*XCV(0]VMIN=0.ENDIFDO301 I2,L2IF(V(,M2)<0.) yMIN-DMAXI(VMIN,-V(,M2) 1 SearchforVminAFL-AFL+RHO(LMI)-XCV()FL-FL+RHO(LMI)*V(LM2)*XCV()FACTOR=FLOWIN/(FL+AFL*VMIN)30LENDDODO3021-212V(LMI)-(V(.M2)+VMIN)*FACTORy=FACTOR·(yM2+mmD302ENDDODO303J=2,M2AKE(L1,-AKE(L2,)! symmetry;decorationforprintoutDIS(LI1,J-DIS(L2,J)303ENDDORETURN图1-6MODULEBOUND定义

(2) MODULE START 在这个模块中，我们定义一些初始值和不变的边界条件，由于入口边界条件不变，因 为壁面上的速度和温度已知，都可以在这个模块中定义。如图 1-5 所示，我们分别对速度 初场、温度初场、k 初场，初场及其相应的边界条件进行处理。 图 1-5 MODULE START 定义 (3) MODULE DENSE 在这个模块中，需要对流体的密度进行定义。由于密度为常数，所以这个模块，我们 无须进行定义，按主程序中的默认值给定。虽然该模块为空块，但必须保留，不能删除， 否则会引起程序错误。这是因为主程序中调用了这个模块。 (4) MODULE BOUND 在这个模块里面，我们对该物理问题的边界条件进行设置。我们拟打算用出口速度提 升的方法，来防止迭代过程中出口法向流速出现负值对计算区域内的影响。具体实施过程 如图 1-6 所示。 图 1-6 MODULE BOUND 定义

(5)MODULEOUTPUT在这个模块中，我们可以输出一些重要的信息来观察程序是否达到收敛。针对这个问题，我们拟输出每个迭代步的SMAX，SSUM（质量守恒残差）以及速度v，温度T，k在某一位置的值来观察送代是否收敛。由于这个问题，速度场和温度场没有耦合，我们可以先求解速度场，等温度场收敛后，我们再进行温度场的求解，因而需要对求解变量进行切换，如图1-7所示。ENTRY OUTPUTIF(ITER-O)THENPRINT401WRITE(8,401)401FORMAT(1X,ITER',6X,SMAX,6X,SSUM,5X,V(6,6)14X,T(5,6),4X,KE(5,6))ELSEPRINT403,ITER, SMAX, SSUM, V(6,6),T(5,),AKE(5,6)WRITE(8,403)ITER,SMAX,SSUM,V(6,6),T(5,6),AKE(5,6)403FORMAT(1X,I6,1P5E11.3)ENDIFIF(ITER>=55)THEN! Switch off the solution variables:LSOLVE(4)-.TRUEFlow is not coupled withLSOLVE(1)-.FALSE.temperature! After obtainingLSOLVE(5)-.FALSE.converged flow field, temperatureLSOLVE(6)-.FALSE.issolvedENDIFIF (ITER-LAST)CALLPRINTRETURN图1-7MODULEOUTPUT定义(6)MODULEGAMSOR在这个模块中，我们分别对控制方程中的当量扩散系数「和源项进行设置。如图1-8所示，我们分别对每个控制方程的当量扩散系数进行了设置，其中流粘性系数可由求解的k和来获得，这里我们忽略了分子粘性系数，因为其值相对于瑞流粘性系数相差好几个数量级。ENTRYGAMSORIF(NF3)RETURNIF(NF-1)THEN! NF=7 for turbulent viscosityREL-1.-RELAX(7)DO500J-1,M1capk2DO 501 I=1,L1.AMT=CMU*RHO(LJ)*AKE(L,J)**2/(DIS(I,J)+1.E-30)6IFTER-0)AMUT(,J)-AMTInitial valuesAMUT(L,J)-RELAX(7)*AMT+REL*AMUT(LJ)501ENDDO!Underrelaxationforturbulentviscosity500ENDDOFACTOR=1ELSEPr=nc,12,A=nc,/PrIF(NF-4)FACTOR=CPCON/PRTIF(NF==5)FACTOR=1./PRK+)-fork.(n +)-foreIF(NF==6)FACTOR=1./PRDn0.a.DO 520 J-1,M1DO 521I=1,L1Laminarpartis omittedGAM(LJ-AMUT(LJ)*FACTORIF(NF/-1) GAM(L1,J)-0. ! Symmetric line, u-0GAM(,M1)-0. ! Local oneway for outlet521ENDDO520ENDDO-图1-8MODULEGAMSOR中当量粘性系数确定

(5) MODULE OUTPUT 在这个模块中，我们可以输出一些重要的信息来观察程序是否达到收敛。针对这个问 题，我们拟输出每个迭代步的 SMAX，SSUM（质量守恒残差）以及速度 v, 温度 T，k 在某 一位置的值来观察迭代是否收敛。由于这个问题，速度场和温度场没有耦合，我们可以先 求解速度场，等温度场收敛后，我们再进行温度场的求解，因而需要对求解变量进行切换， 如图 1-7 所示。 图 1-7 MODULE OUTPUT 定义 (6) MODULE GAMSOR 在这个模块中，我们分别对控制方程中的当量扩散系数和源项进行设置。如图 1-8 所 示，我们分别对每个控制方程的当量扩散系数进行了设置，其中湍流粘性系数可由求解的 k 和来获得，这里我们忽略了分子粘性系数，因为其值相对于湍流粘性系数相差好几个数 量级。 图 1-8 MODULE GAMSOR 中当量粘性系数确定

壁面处的当量扩散系数需要应用壁面函数法来确定，同时需要运用壁面函数法对壁面处各变量的边界条件进行设置，具体实施过程如图1-9所示。WFM implementation!WFDO530J-2,M2!Foru,p',k,e-SELECTCASE(NE)MFor u, p'and kadiabatic; For &, set up its valueCASE (1,3,5,)of lst innernode,then cut the connection toitsGAM(1,J)-0. boundary.All lead to F-0mCASE (2)!For velocity y, WFM should be used!pGAM(1,J)=AMU ! First, laminar viscosityis given for the left wallXPLUS(J-RHO(2J)*SQRT(AKE(2.)*CMU4*XDIF(2)/AMUemIF(XPLUS(J>11.5)GAM(1,J-AMU*XPLUS(J)e1 (ALOG(9.*XPLUS()*2.5) ! Turbulence viscositynCASE (4) ! Fortemperature, WFM for temperature1GAM(1,J)-AMU*CPCON/PR!First, laminarthermal conductivityatIF(XPLUS(J)>11.5) GAM(1,J=AMU*CPCON/PRT*XPLUS(J)1 /(2.5*ALOG(9.*XPLUS()+PFN) ! Turbulence thermal conductivityoENDSELECTpx(C1/4k1/2)nPr+530ENDDOn图1-9MODULEGAMSOR中的壁面函数法针对端流流动问题，由于源项较为复杂，因而程序的很大一部分是对各方程的源项进行处理。以速度方程u源项为例，图1-10介绍如何数值处理u方程源项。对于方程，由于靠近壁面的第一节内节点其值给定，我们采用大源项法进行处理，具体处理过程见图1-11。aaOuOvXDIF(O)IF(NF=-1) THENS.(ll,Cl(+)-(iOxOxOyOx(GAM(I.J)DO590J=2,M2xer(0)DO591I=3L2GAM(I-1,))-(-LD)CON(LJ)-(GAMLJ)*(U(I+1,J-U(L)/XCVM)xer(i-1)1-GAM(I-1,J)*(U(ILJ-U(I-1,J)/XCV(I-1)/XDIF()SourceGAMP=GAM(I,J+1)*GAM(I-1,J+1)/(GAM(/J+1)+GAM(I-1,J+1)+1.E-30)termforuGAMP=GAMP+GAM(I,J)*GAM(I-1,J)/(AM(I,J)+GAM(I-1,J)+1.E-30)eqGAMM-GAM(L,J-1)*GAM(I-1,J-1)/(GM(L,J-1)+GAM(I-1,J-1)+1.E-30)GAMM-AMM+GAM(I,J)*GAM(/,J/(GAM(I,)+GAM(I-1,J)+1.E-30)CON(L,J=4ON(I,J)+(GAMP*(V(I,J+1)-V(I-1,J+1)1-GAMM*(V-V(I-1,J)/(YCV(J)*XDIF())AP(L,J)=0.591ENDDOn(i-l)n(i)n(i-1j+1)n(i.j+1)tne590ENDDOn.(i-1, j)+n.(i,j) n.(i-1.j+1)+n.(i.j+1)RETURNRefertotextbookpage358图1-10MODULEGAMSOR中的u源项处理

壁面处的当量扩散系数需要应用壁面函数法来确定，同时需要运用壁面函数法对壁面 处各变量的边界条件进行设置，具体实施过程如图 1-9 所示。 图 1-9 MODULE GAMSOR 中的壁面函数法 针对湍流流动问题，由于源项较为复杂，因而程序的很大一部分是对各方程的源项进 行处理。以速度方程 u 源项为例，图 1-10 介绍如何数值处理 u 方程源项。对于方程，由于 靠近壁面的第一节内节点其值给定，我们采用大源项法进行处理，具体处理过程见图 1-11。 图 1-10 MODULE GAMSOR 中的 u 源项处理

S,=Apghen,S, =-A,Gen,GS=GenGC2peCp808A=1020~1030kk飞Large source term methodDO 600 J-2,M2DO601I-2L2CON(I,)=C1*GEN(I,J)*CMU*RHO(I,J)*AKE(I,J)?AP(L,J--C2*RHO(L,J)*DIS(L,J)/(AKE(L,J)+1.E-30)601ENDDOC314k,32600ENDDOKYPDO 602 J=2,M2DISS=CMU*AKE(2,J)**1.5/(0.4*CMU4*XDIF(2)福CON(2,J)-1.E30*DISSAP(2,J)--1.E30Adopt large source602ENDDOterm method for istSource termRETURNinnernodewherei-2;calculation forENDEpsilon eq.图1-11MODULEGAMSOR中的e源项处理能力提高实践案例二：歧管微通道中流动换热案例(二)1.问题描述随着电子设备集成化和小型化的快速发展，高集成电路和微电子器件得到广泛应用，导致电子设备的热流通量不断增加。调查显示，超过50%的电气设备故障是由过热引起的。微电子器件中巨大的散热对电子设备的热管理提出了严峻的挑战。然而，传统的空气冷却技术难以满足如此巨大的需求。为了解决这一问题，微通道散热器的概念被提出，因其优异的传热性能成为了最有前途的冷却技术之一吸引了大量的研究。本节通过一种多射流微通道(MJMC)散热器实例，介绍Fluent中基于有限体积法(FVM)采用SIMPLE算法对MIMC内的三维单相流动和换热过程进行数值模拟的工作步骤，MJMC单元的整体结构如图2-1(a)所示，选取MJMC单元单元的1/4区域图2-1(b)，进行流体流动和换热特性的数值研究。模型的边界条件设置如下：（1）由于计算域相对于单元整体的对称结构，1/4单元域左端面和前、后端面均采用对称边界条件：（2）基于单元体周期性重复的进出口结构，1/4MJMC单元中方存在5个交替分布的射流入口和出口，如图2-1a)所示：（3）通道单元上盖板的导热系数较低，模拟中将盖板与其他部件的接触面近似假设为绝热面，即计算域的上端面和右端面设置为绝热：（4）计算域底面通过恒定热流加热以模拟电子元件的发热，模拟中将热流密度设定为100Wcm-2。图2-1(c)为MJMC的截面视图，模型的具体几何尺寸参数设置如图所示

图 1-11 MODULE GAMSOR 中的 源项处理 (二) 能力提高实践案例二：歧管微通道中流动换热案例 1. 问题描述 随着电子设备集成化和小型化的快速发展，高集成电路和微电子器件得到广泛应用， 导致电子设备的热流通量不断增加。调查显示，超过 50%的电气设备故障是由过热引起的。 微电子器件中巨大的散热对电子设备的热管理提出了严峻的挑战。然而，传统的空气冷却 技术难以满足如此巨大的需求。为了解决这一问题，微通道散热器的概念被提出，因其优 异的传热性能成为了最有前途的冷却技术之一吸引了大量的研究。 本节通过一种多射流微通道(MJMC)散热器实例，介绍 Fluent 中基于有限体积法(FVM) 采用 SIMPLE 算法对 MJMC 内的三维单相流动和换热过程进行数值模拟的工作步骤。 MJMC 单元的整体结构如图 2-1(a)所示，选取 MJMC 单元单元的 1 / 4 区域图 2-1(b)，进行 流体流动和换热特性的数值研究。模型的边界条件设置如下：(1) 由于计算域相对于单元整 体的对称结构，1/4 单元域左端面和前、后端面均采用对称边界条件；(2) 基于单元体周期 性重复的进出口结构，1/4MJMC单元中方存在5个交替分布的射流入口和出口，如图2-1(a) 所示；(3) 通道单元上盖板的导热系数较低，模拟中将盖板与其他部件的接触面近似假设为 绝热面，即计算域的上端面和右端面设置为绝热；(4) 计算域底面通过恒定热流加热以模拟 电子元件的发热，模拟中将热流密度设定为 100 W cm−2。图 2-1(c)为 MJMC 的截面视图， 模型的具体几何尺寸参数设置如图所示