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

流动传热问题数值仿真多层次案例*一、摘要本教学案例设计坚持“紧跟学科发展、突出能力提高”，面向能源、航天航空、化工等领域国家重大需求，科教相长，及时更新研究生课程案例教学内容，将流动传热问题数值仿真教学案例设计与国家重大需求高度匹配，设计了包含流动传热数值模拟自编程基础案例以及面向芯片热管理、电池热管理、太阳能高效利用等为背景的能力提高实践案例。通过上述多层次的案例，既夯实研究生基础理论知识，又培养研究生运用数值方法解决实际工程问题的能力，达到培养研究生的创新意识和创新能力的目的。关键词：流动传热、数值仿真、电池热管理、太阳能高效利用*负责人简介：陶文铨，男，西安交通大学教授，博士生导师，中科院院士。编制说明：本教学案例将流动传热问题数值仿真教学案例设计与国家重大需求高度匹配，设计了包含流动传热数值模拟自编程基础案例以及面向芯片热管理、电池热管理、太阳能高效利用等为背景的能力提高实践案例。夯实研究生基础理论知识，培养研究生运用数值方法解决实际工程问题的能力

流动传热问题数值仿真多层次案例* 一、 摘要 本教学案例设计坚持“紧跟学科发展、突出能力提高”，面向能源、航天航空、化工等 领域国家重大需求，科教相长，及时更新研究生课程案例教学内容，将流动传热问题数值 仿真教学案例设计与国家重大需求高度匹配，设计了包含流动传热数值模拟自编程基础案 例以及面向芯片热管理、电池热管理、太阳能高效利用等为背景的能力提高实践案例。通 过上述多层次的案例，既夯实研究生基础理论知识，又培养研究生运用数值方法解决实际 工程问题的能力，达到培养研究生的创新意识和创新能力的目的。 关键词：流动传热、数值仿真、电池热管理、太阳能高效利用 * 负责人简介：陶文铨，男，西安交通大学教授，博士生导师，中科院院士。 编制说明：本教学案例将流动传热问题数值仿真教学案例设计与国家重大需求高度匹配，设计了包含流动 传热数值模拟自编程基础案例以及面向芯片热管理、电池热管理、太阳能高效利用等为背景的能力提高实 践案例。夯实研究生基础理论知识，培养研究生运用数值方法解决实际工程问题的能力

二、背景信息流动传热过程是能源与环境领域的关键过程。流动传热过程的数值仿真及其强化控制是新时期我国经济社会可持续发展、国防安全的重要保证。1）随着学科的发展，针对极端条件下（高温、高压、高超声速、微纳尺度等）的流动换热问题的研究日益增多，极端条件给实验教学带来了难度。编写流动传热问题的数值仿真案例，具有降低实验危险、减少资源浪费、降低成本的优势。2）近20年来，计算机硬件工业的发展为数值传热学提供了坚实的物质基础。随着先进数值方法的不断发展，数值仿真技术和理论分析和实验一起，成为热流问题研究的三种重要手段。在某些极端条件下（如高温高压、微纳尺度等），数值仿真甚至是主要的研究手段。数值仿真具有成本较低和能模拟复杂或较理想过程等优点。数值仿真对流动与传热过程的研究发挥了重要的作用。因此，流动传热数值模拟已发展成为本领域毕业生必须熟练掌握的一项技术。3）最后特别需要指出的是，我国目前广泛使用的商业计算流体力学软件如Fluent、Comsol、Star-CD等产权均为国外所有，目前尚无具有自主知识产权的成熟CFD软件。随着我国科研中“卡脖子”技术负面影响的日益凹陷，非常有必要在本科教学中进行数值仿真的数值理论及实践教学，为我国自主研发CFD软件储备人才。综上所述，通过实体实验方法进行教学已经无法满足能源工程领域高速发展的潮流，通过数值仿真和传统实验相结合的方法对研究生教学实践环节进行改革势在必行。本案例依托我校优势学科，服务《数值传热学》、《计算传热学近代进展》、《强化传热原理与技术》等课程，综合具有自主软件著作权的开源计算流体力学和传热学程序并整合商用CFD软件，构建数值仿真方法理论基础及实践教学的多层级综合案例，通过案例教学旨在夯实学生专业知识及数值仿真理论基础并培养学生创新实践能力，培养基础知识扎实、实践能力强创新型人才，为我国重点战略领域如航天航空发动机设计、核能反应堆设计、新能源汽车热管理等储备高水平青年人才，提高我国核心领域的原始创新力和核心竞争力。三、案例正文（一）基础案例：二维椭圆型流动与传热问题通用程序案例为便于学生掌握数值求解流动与传热问题的方法，全面提高自主编程能力，本教学案例将二维圆型流动与传热问题的通用程序在西安交通大学数值传热学研究小组网页上公布，并附带八个教学案例，包括：（1）直角坐标中二维稳态无内热源的导热；（2）空心圆柱内的稳态热传导：（3）正方形管道内的充分发展对流换热；（4）内壁上有直肋的环形通道内充分发展对流换热：（5）二维突扩通道中的流动与换热：（6）方形通道内的复杂充分流动：（7）旋转圆盘上的冲击流动：（8）带中心射流通道内的瑞流换热。该教学程序，最初来源于美国明尼苏达大学的Patankar教授。原教学程序所采用语言为Fortran77，本教学团队将其改为Fortran95，并做了相应改进：采用了新的通用控制方程，便于处理流固耦合传热界面，同时可以加快收敛速度

二、 背景信息 流动传热过程是能源与环境领域的关键过程。流动传热过程的数值仿真及其强化控制 是新时期我国经济社会可持续发展、国防安全的重要保证。 1）随着学科的发展，针对极端条件下（高温、高压、高超声速、微纳尺度等）的流动 换热问题的研究日益增多，极端条件给实验教学带来了难度。编写流动传热问题的数值仿 真案例，具有降低实验危险、减少资源浪费、降低成本的优势。 2）近 20 年来，计算机硬件工业的发展为数值传热学提供了坚实的物质基础。随着先 进数值方法的不断发展，数值仿真技术和理论分析和实验一起，成为热流问题研究的三种 重要手段。在某些极端条件下（如高温高压、微纳尺度等），数值仿真甚至是主要的研究手 段。数值仿真具有成本较低和能模拟复杂或较理想过程等优点。数值仿真对流动与传热过 程的研究发挥了重要的作用。因此，流动传热数值模拟已发展成为本领域毕业生必须熟练 掌握的一项技术。 3）最后特别需要指出的是，我国目前广泛使用的商业计算流体力学软件如 Fluent、 Comsol、Star-CD 等产权均为国外所有，目前尚无具有自主知识产权的成熟 CFD 软件。随 着我国科研中“卡脖子”技术负面影响的日益凹陷，非常有必要在本科教学中进行数值仿真 的数值理论及实践教学，为我国自主研发 CFD 软件储备人才。 综上所述，通过实体实验方法进行教学已经无法满足能源工程领域高速发展的潮流， 通过数值仿真和传统实验相结合的方法对研究生教学实践环节进行改革势在必行。本案例 依托我校优势学科，服务《数值传热学》、《计算传热学近代进展》、《强化传热原理与技术》 等课程，综合具有自主软件著作权的开源计算流体力学和传热学程序并整合商用 CFD 软件， 构建数值仿真方法理论基础及实践教学的多层级综合案例，通过案例教学旨在夯实学生专 业知识及数值仿真理论基础并培养学生创新实践能力，培养基础知识扎实、实践能力强、 创新型人才，为我国重点战略领域如航天航空发动机设计、核能反应堆设计、新能源汽车 热管理等储备高水平青年人才，提高我国核心领域的原始创新力和核心竞争力。 三、 案例正文 (一) 基础案例：二维椭圆型流动与传热问题通用程序案例 为便于学生掌握数值求解流动与传热问题的方法，全面提高自主编程能力，本教学案 例将二维椭圆型流动与传热问题的通用程序在西安交通大学数值传热学研究小组网页上公 布，并附带八个教学案例，包括：(1) 直角坐标中二维稳态无内热源的导热; (2) 空心圆柱内 的稳态热传导；（3）正方形管道内的充分发展对流换热；（4）内壁上有直肋的环形通道 内充分发展对流换热；（5）二维突扩通道中的流动与换热；（6）方形通道内的复杂充分 流动；（7）旋转圆盘上的冲击流动；（8）带中心射流通道内的湍流换热。该教学程序， 最初来源于美国明尼苏达大学的 Patankar 教授。原教学程序所采用语言为 Fortran77, 本教学 团队将其改为 Fortran95，并做了相应改进：采用了新的通用控制方程，便于处理流固耦合 传热界面，同时可以加快收敛速度

该程序分为两大部分：主程序和用户子程序。其中主程序在课程学习期间，用户不需要对其做改动，主要功能是针对提供的网格节点信息生成网格及相应的数组，建立节点间的代数方程，以及对代数方程进行求解。用户子程序需要根据具体物理问题进行相应设置，为主程序提供网格信息，初始条件，流体密度，边界条件，当量扩散系数和源性等信息。该教学程序的主要特点如下：采用有限体积法来导出离散方程，基于SIMPLER算法处理压力和速度的耦合关系，以原始变量UV,P作为流场求解变量并采用交叉网格：对流-扩散项离散采用乘方格式，但对流项可方便地改为任一种三点格式：源项采用局部线性化方式处理，界面上当量扩散系数采用调和平均方式计算：非稳态问题采用全隐格式，第二、三类边界条件采用附加源项法进行处理：代数方程求解采用送代法，用交替方向隐式线送代方法辅以块修正技术求解代数方程。在课程结束而且学生对教学程序已经充分理解后，学生自已可对主程序进行修改，以适应求解更为复杂的问题，例如从二维改为三维。这是培养学生自主开发软件能力的一种很好的训练。自本课程开设以来，已经有不少研究生基于本教学程序完成了自己的学位论文。下面以带中心射流通道内的端流换热的案例为例，介绍如何运用教学程序，标准k-E模型和壁面函数法来求解流流动传热问题。标准k-c模型和壁面函数法是目前求解端流中最常用的模型之一，也是商业软件中经常采用的模型。通过教学程序，可以更加深入全面的认识这一物理过程。1.物理问题描述如图1-1所示，带一般中心射流的流体进入一平行板通道。流动为紊流，动力粘度n-10-6，Pr=0.7。请采用标准k-c模型和壁面函数法，确定突扩区域中的流场和温度场。1.04.0Tw=0T=100T=400V=100V=10图1-1带中心射流的通道内的端流换热2.数学描写

该程序分为两大部分：主程序和用户子程序。其中主程序在课程学习期间，用户不需 要对其做改动，主要功能是针对提供的网格节点信息生成网格及相应的数组，建立节点间 的代数方程，以及对代数方程进行求解。用户子程序需要根据具体物理问题进行相应设置， 为主程序提供网格信息，初始条件，流体密度，边界条件，当量扩散系数和源性等信息。 该教学程序的主要特点如下：采用有限体积法来导出离散方程，基于 SIMPLER算法处理压 力和速度的耦合关系，以原始变量 u,v,p 作为流场求解变量并采用交叉网格；对流-扩散项 离散采用乘方格式，但对流项可方便地改为任一种三点格式；源项采用局部线性化方式处 理，界面上当量扩散系数采用调和平均方式计算；非稳态问题采用全隐格式，第二、三类 边界条件采用附加源项法进行处理；代数方程求解采用迭代法，用交替方向隐式线迭代方 法辅以块修正技术求解代数方程。在课程结束而且学生对教学程序已经充分理解后，学生 自己可对主程序进行修改，以适应求解更为复杂的问题，例如从二维改为三维。这是培养 学生自主开发软件能力的一种很好的训练。自本课程开设以来，已经有不少研究生基于本 教学程序完成了自己的学位论文。 下面以‘带中心射流通道内的湍流换热’的案例为例，介绍如何运用教学程序，标准 k- 模型和壁面函数法来求解湍流流动传热问题。标准 k-模型和壁面函数法是目前求解湍流中 最常用的模型之一，也是商业软件中经常采用的模型。通过教学程序，可以更加深入全面 的认识这一物理过程。 1. 物理问题描述 如图 1-1 所示，带一般中心射流的流体进入一平行板通道。流动为紊流，动力粘度 ＝ 10-6 , Pr = 0.7。请采用标准 k-模型和壁面函数法，确定突扩区域中的流场和温度场。 图 1-1 带中心射流的通道内的湍流换热 2. 数学描写

完整的数学描写包括控制方程和边界条件，现在分别对控制方程和边界条件进行描述，假设流动为稳态，则可以忽略非稳态项，各控制方程可用如下统一的形式来表达：div(p up)=div(Tgradp)+S*(1-1)其中，分别为u,v,T,,c,p,p。各个变量的当量扩散系数如表1-1所示：表1-1各变量当量扩散系数和松弛因子选取1234NF=567811VPcTkU8PVariableGntn.cpn.n.<Q-ntntPr,Ok0eα0.80.81.00.60.60.6对于u控制方程，其源项为：aou)OvS.=-(Men ax+(Menaxaxay(1-2)对于控制方程，其源项为aauaav(1-3)S,=(enr1-(uenayayaxay对于k控制方程，其源项为S=n,G-ps(1-4)对于控制方程，其源项为S, =GenG_ Spe?(1-5)kk其中G为： (%+))(1-6)边界条件为：(1)入口处：k取来流平均动能的1%；ε按照ε=cupk2n,取值

完整的数学描写包括控制方程和边界条件，现在分别对控制方程和边界条件进行描述。 假设流动为稳态，则可以忽略非稳态项，各控制方程可用如下统一的形式来表达： * * div u div ( ) ( grad ) S    =  +  (1-1) 其中， 分别为 u, v, T, k, , p, p’。各个变量的当量扩散系数如表 1-1 所示： 表 1-1 各变量当量扩散系数和松弛因子选取 对于 u 控制方程，其源项为： eff eff ( ) ( ) u u v S x x y x       = +     (1-2) 对于 v 控制方程，其源项为 eff eff ( ) ( ) v u v S x y y y       = +     (1-3) 对于 k 控制方程，其源项为 k t S G = −   (1-4) 对于 控制方程，其源项为 2 1 t 2 c G c S k k    = − (1-5) 其中 G 为： 2 2 2 2 t u v u v G x y y x                 = + + +                       (1-6) 边界条件为： (1) 入口处：k 取来流平均动能的 1%； 按照 2 t c k     = 取值

(2）壁面处：采用壁面函数法确定。对于垂直于壁面的速度，其法向一阶导数为0；对于平行于壁面的速度，设置为0，壁面处的等效端流粘度系数和流导热系数则由壁面函数法确定：对于k，其法向一阶导数为0；对于，其第一个与壁面相邻近的内节点给定为C=C3l*kp3'2 /kyp.(3）出口处按照局部单向化处理；（4）对称线上，t=0,其余法向一阶导数为03.数值方法(1）k-a方程源项的处理将k-s方程源项局部线性化S=S。+Sp处理，如图1-2所示。引入一个负S，来增加数值计算的稳定性。S, =n,G- pe=n,G-(PE)kKScSpCpeS-SenGc,en,GCpe18kkkkScSp图1-2k-s方程源项局部线性化处理方法(2）出口速度的提升采用局部单向化条件时，为防止送代过程中出口法向流速出现负值对计算区域内的影响，采用提升出口界面流速的方法避免。FLOWINFACTOR=(1-7)L2Z[(V,M2 +Vmm)* RHOi,M * XCV(i)]=2V,m = FACTOR-(v,M2 +/vmin)(1-8)(3）动量方程源项的处理本程序采用交叉网格，速度存储在界面上，而其他物理信息存储在主节点上。在对动量方程源项处理时，需要考虑速度位置和节点间的关系，以u控制方程为例，其动量方程源项的处理方法如图1-3所示

(2) 壁面处：采用壁面函数法确定。对于垂直于壁面的速度，其法向一阶导数为 0；对于平 行于壁面的速度，设置为 0，壁面处的等效湍流粘度系数和湍流导热系数则由壁面函数 法确定；对于 k, 其法向一阶导数为 0；对于，其第一个与壁面相邻近的内节点给定为 3/4 3/2 C k y  P P   = 。 (3) 出口处按照局部单向化处理； (4) 对称线上，u=0, 其余法向一阶导数为 0. 3. 数值方法 (1) k-方程源项的处理 将 k-方程源项局部线性化 S = Sc + Sp处理，如图 1-2 所示。引入一个负 Sp 来增加数 值计算的稳定性。 图 1-2 k-方程源项局部线性化处理方法 (2) 出口速度的提升 采用局部单向化条件时，为防止迭代过程中出口法向流速出现负值对计算区域内的影 响，采用提升出口界面流速的方法避免。 min 2 , 2 , 1 2 [( ) ( )] L i M i M i FLOWIN FACTOR V RH V O XCV i = =  +   (1-7) , 1 , 2 min ( ) i M i M v FACTOR v v =  + (1-8) (3) 动量方程源项的处理 本程序采用交叉网格，速度存储在界面上，而其他物理信息存储在主节点上。在对动 量方程源项处理时，需要考虑速度位置和节点间的关系，以 u 控制方程为例，其动量方程 源项的处理方法如图 1-3 所示

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 的截面视图， 模型的具体几何尺寸参数设置如图所示