1.2　Mechanical模块与CFD模块在热分析中的区别

通过前例可知，Mechanical模块与CFD模块均可以进行热分析，两者的区别如下。

1．计算原理不同

Mechanical模块采用有限单元法（FEM），Fluent模块采用有限体积法（FVM）的单元中心（cell-centered）方式，CFX模块采用有限体积法的节点中心（node-centered）方式。

注意

在FEM中将划分网格形成的块称为Element，在FVM中称为Cell。

下面以通量平衡基础方程说明FEM与FVM的区别：

式中， u 为任意守恒物理量， 为该物理量的通量，即单位面积、单位时间内通过 域控制面的流动量。

采用有限元形式求解此类偏微分方程，并不拘泥强形式微分方程解，但至少满足积分点上的条件，即弱形式微分方程解。因此对模型 域进行积分建立方程：

式中， 为有限元计算加权试函数， V 为 域体积。

为了控制通量与模型内部向量的联系，对 采用散度定律，则：

式中， 为 域的边界， n 为边界上的法向量， S 为 域内封闭曲面的面积。

同时

则

代入模型 域的积分方程，得：

该式即为FEM对偏微分求解的基本方程。

当试函数与形函数一致时，即为伽辽金法。形函数只与几何、物理量相关，一般设为多项式，也可以定义为常数，当 时，则上式为：

该式即为FVM对偏微分求解的基本方程。

由此可知，FVM只是FEM求解中的一个特例，不同之处在于，FEM选取离散后有限数量的形函数，FVM选取离散后有限数量的体积域，分别使其满足控制域内的基本方程。

图1-2-1表示：对于内部节点，只有周围6个网格剖面线的单元（Element）与 域的域积分有关联；对于边界上（ ）的节点，只有周围3个单元（Element）与 域的域积分有关联，只有周围2个斜线剖面线的单元（Element）与 的域积分有关联。

图1-2-2表示：每个单元（Cell）都被视为一个域，不管是内部还是边界上的单元（Cell），都与 域的域积分有关联。

图1-2-3表示：在内部以每条边的中点加上每个单元（Cell）的中心构成一个域；在边界处除了边的中点和单元（Cell）的中心，再以边界进行封闭构成一个域（如果存在边界层，且为三角形网格，对于直角和钝角三角形单元，则以三角形最长边的中点加边的中点构成域）。不管是内部还是边界上的域，都与 域的域积分有关联。

对比FEM和FVM离散原理可知，FEM可以保证全局通量守恒，但不能保证局部通量守恒。换言之，在域边界处通量保证平衡，但是难以控制内部局部通量。以热分析为例，纯传导分析易于计算，FEM可以采用较粗的二次单元完成计算且保证内部热通量平衡；但是针对对流、辐射为主导的热分析，虽然可以在边界处得到精度较高的结果，但是内部热量流动过程较难处理。

FVM可以保证局部通量守恒，但是在边界上的解难以保证精度，即便在边界处划分更多的网格也不行，除了降低效率，还会影响边界上整体与全域内局部的匹配。因此在热分析中，FVM在边界处通常都存在一些热不平衡量误差，但是针对对流、辐射为主导的热分析，可以很容易得到内部热量流动过程。

2．分析的侧重点不同

Mechanical模块利用单位等效类比方式将结构有限元分析扩展到热学有限元分析，即将结构分析中的位移自由度变为热学分析中的温度自由度，同理将结构中的其他物理量变为热学中的物理量，如表1-2-1所示。

热传导计算公式为：

式中， q 为热通量， Q 为热量， A 为截面积， k 为热传导率， d 为温差面之间的距离， T 为温度。

对流计算公式为：

式中， q 为热通量， h 为对流换热系数（W·m ⁻² ·K ⁻¹ ）， T 为温度。

辐射计算公式为：

式中， Q 为热量， 为辐射率（0～1）， 为斯特藩-玻尔兹曼常数（约为5.67×10 ⁻⁸ W·m ⁻² ·K ⁻⁴ ）， A 为截面积， F 为辐射角系数， T 为温度。

对于传导过程，采用FEM可以快捷、准确地完成计算，完全不需要CFD计算；对于对流过程，如果对流换热系数为已知条件，FEM比CFD更方便计算，但是对流换热系数往往比较难获得一个较为准确的平均值，即便定义对流换热系数相对于时间或温度的函数，也存在一定的误差，所以对于强制对流形式，对流换热系数相对容易获得，可以采用CFD为主FEM为辅的计算形式，对于自然对流形式，对流换热系数较难获得，可以采用CFD计算；对于辐射过程，关键难度在于获得任意表面的辐射角系数，FEM将其过程大大简化，仅定义辐射率和Hemicube分辨率控制辐射所有参数，因此如果需要考虑更多的辐射因素，例如光谱依赖性、半透明固体的吸收、辐射热通量的平滑分布、太阳辐射等，采用CFD更加可靠。

下面以自然对流形式简要说明CFD分析流程。已知一个地暖区域，其二维截面积如图1-2-4所示，下边温度为315 K，周围环境温度为290 K，求稳态条件下该区域的对流换热系数。

对流换热由高温模型边界外流体运动产生，将较热的流体从边界层表面输送出去，该区域空间由较冷的流体补充，以至循环。这种微观运移简化为宏观描述，即为对流。本例中没有其他外界因素促使该区域内空气流动，为自然对流。自然对流中流体循环的动力仅来源于热流和冷流的密度差异，重力起到了非常重要的作用。

在自然对流计算过程中，采用瑞利数作为流动特征评估参数，不采用雷诺数。瑞利数计算公式为：

（ ，层流； ，湍流）

式中， （格拉晓夫数），用于定义浮力与黏性力之比，其中 为热膨胀系数，单位为℃ ⁻¹ ， g 为重力加速度，单位为m·s ⁻² ， L 为垂直于温度梯度的长度，单位为m，Δ T 为温差，单位为℃， 为流体运动黏度，单位为m ² ·s ⁻¹ ； （普朗特数），用于定义动量扩散与温度扩散之比，其中 为流体热扩散系数，单位为m ² ·s ⁻¹ ，在计算 的这个公式中 k 为流体热传导系数，单位为W·m ⁻¹ ·K ⁻¹ ， 为流体密度，单位为kg·m ⁻³ ， C 为比热容，单位为J·kg ⁻¹ ·K ⁻¹ 。

本例的物理参数如表1-2-2所示。

解：

将数值代入，则 ，流动形式为层流。

Nu （努塞尔特数）计算公式为： ，式中参数如表1-2-3所示。

本例中取 C =0.27， n =1/4，则 Nu =34.93；对流换热系数计算公式为 ，其中 L 对于竖板模型为模型上垂直于温度梯度的长度，对于平板模型为面积与周长之比，本例中 L =0.167 m，则 h =5.062 W·m ⁻² ·K ⁻¹ 。

下面用ANSYS Workbench的Fluent模块进行计算，以验证其计算结果。

（1）建立分析流程

如图1-2-5所示，建立热分析流程。其中包括A框架结构的划分网格、B框架结构的Fluent热分析。

注意

本例采用A3（Mesh）选项与Fluent的B2（Setup）选项建立关联的分析流程，也可以采用Geometry组合Fluent（with Fluent Meshing）模块的分析流程。Fluent（with Fluent Meshing）模块可以在开始时定义划分网格和求解计算的硬件资源，在单独使用Fluent进行大型工程分析时建议采用。

（2）前处理

建模过程省略。

如图1-2-6所示划分网格，其中Details窗口中的Physics Preference处选择CFD，Solver Preference处选择Fluent；插入Sizing项对整个面定义网格尺寸，在Geometry处选择整个面，Type处选择Element Size，Element Size设置为0.005 m；插入Inflation定义边界层，在Geometry处选择整个面，Boundary处选择模型的4条边，Maximum Layers设置为10。

划分网格后，在Details窗口的Quality→Mesh Metric处查看网格质量。在Mesh Metric处选择Aspect Ratio（网格纵横比）。

注意

可以在Smoothing处选择High项进一步提高网格质量。

注意

对于自然对流分析，必须定义边界层，一般定义5～10层，网格质量重点为纵横比，必须小于40。

如图1-2-7所示定义命名选择，以便Fluent定义边界条件。其中上边线命名为“top”，下边线命名为“bottom”，左边线命名为“left”，右边线命名为“right”，整个面命名为“air-face”。

（3）Fluent热分析流程

在B2 Setup处双击鼠标左键，保留默认设置（2D模型），单击Start进入Fluent分析模块，如图1-2-8所示，在1区选中Gravity复选框，在Y方向定义重力加速度为−9.81（自然对流必须设置）；在2区设置Energy为On，Viscous为Laminar（层流）；在3区双击Materials→Fluid→air，在弹出的对话框中修改Density为boussinesq，设置数值为1.1405 kg·m ⁻³ ，设置Thermal Expansion Coefficient的数值为0.00344827 K ⁻¹ （ 切记不要遗漏 ），其余参数均为默认设置，修改完成单击Change/Create按钮保存；在4区Cell Zone Conditions处选择Fluid（定义为流体域，软件可能默认为固体域）；在5区Boundary Conditions处双击，在弹出的对话框中单击Operating Conditions…按钮，在出现的对话框中设置Operating Pressure为0 Pa，Operating Temperature为290 K。

如图1-2-9所示，在Boundary Conditions→wall项中依次双击bottom、left、right、top（依据命名选择自动创建边界条件项），对“bottom”边在Thermal选项卡中设置Temperature为315 K；对“left”和“right”边，在Thermal选项卡中设置Heat Flux为0；对“top”边，在Thermal选项卡中设置Temperature为290 K。其余全部采用默认设置。

注意

对于自然对流分析，必须定义重力。在定义流体密度参数时，一般采用Incompressible ideal gas（设定密度，推荐）或Boussinesq（单击Operating Conditions…按钮，在弹出的对话框的Operating Temperature处设置温度，温差在20%以内）形式，就计算精度而言，Boussinesq计算结果会导致羽流宽度变窄，存在一定误差。本例仅展示计算流程，所以采用Boussinesq。因为本计算模型为开放区域，所以Operating Pressure设置为0 Pa。

（4）Fluent热分析求解设置

如图1-2-10所示，在1区双击Solution→Methods，在弹出的对话框中，Scheme下拉列表框选择Coupled，Gradient下拉列表框选择Least Squares Cell Based，Pressure下拉列表框选择PRESTO！，Momentum下拉列表框选择First Order Upwind，Energy下拉列表框选择First Order Upwind，选中Pseudo Transient复选框；在2区双击Solution→Controls，在弹出的对话框中所有选项均保留默认设置；在3区双击Monitors→Residual，在弹出的对话框中，将Continuity的Absolute Criteria项设置为1E−5，其余项均为默认设置。

注意

对于自然对流分析，因为流场和能量为强耦合，所以Scheme下拉列表框选择Coupled，如果采用SIMPLE也可以计算收敛，但计算效率和精度均低于Coupled。同样对于自然对流，Pressure下拉列表框只能选择PRESTO!或者Body-Force Weighted，其中前者应用范围更广一些，主要用于压力梯度较大的计算中，且计算效率高于后者。First Order Upwind与Second Order Upwind的主要区别在于收敛速度，前者容易收敛，但精度略逊于后者。选中Pseudo Transient（伪瞬态）复选框主要为了提高较大纵横比模型的收敛性能。

Solution Controls用于定义求解稳定性的亚松弛因子，该值与最终收敛结果无关，Fluent默认值可以用于大多数场合，如有必要可以调节（调小，提高收敛稳定性；调大，提高收敛速度）。

Fluent稳态分析收敛判断基于3点：①计算残差小于设定的标准残差；②计算域内的质量、动量、能量和标量达到平衡；③某个需要监控的物理量达到稳定状态。对于本例而言，如果仅以小于标准残差或系统平衡判定收敛，则容易出现很大的计算误差。在自然对流中，“left”边和“right”边的温度保持稳定才是该稳态分析的关键。如图1-2-11所示，先定义监控项，右击Report Definitions项，在出现的快捷菜单中选择New→Surface Report→Area Weighted Average…，在弹出的对话框中进行设置，Name文本框采用默认命名，在Field Variable下方的下拉列表中分别选择Temperature…和Static Temperature，其中“report-def-0”对应的Surfaces选择为“left”，“report-def-1”对应Surfaces选择为“right”。

如图1-2-12所示定义收敛监控项，右击Monitors→Report plots项，在出现的快捷菜单中选择New，在弹出的对话框中进行设置，Name文本框采用默认命名，分别单击Add按钮添加“report-def-0”和“report-def-1”，其余项均为默认设置。

单击Run Calculation，将Number of Iteration设置为500，其余项采用默认设置，计算完成后的监控图如图1-2-13所示。可以看到，在400步左右，监控左右两边线温度达到稳定。

（5）Fluent热分析后处理

自然对流计算完成，先查看速度与温度分布云图，如图1-2-14所示。右击Contours项，在出现的快捷菜单中选择New，在弹出的对话框中进行设置，Contour Name文本框采用默认命名，在Contours of…下方的下拉列表框中选择“Velocity…”“Velocity Magnitude”和“Temperature…”“Static Temperature”，分别创建“contour-1”和“contour-2”，在Surfaces下方的列表框内选择“surface_body”，单击Compute按钮，最后单击Save/Display按钮即可。采用类似操作，生成速度向量图，如图1-2-15所示。

由图可见，从“bottom”边界到“top”边界的自然对流速度和热边界层的运移趋势，在“right”边界旁边形成浮力羽流，在“left”边界旁边形成浮射流。浮力羽流形成是因为流体初始动量很小（图中显示速度较慢，0.1 m·s ⁻¹ 左右），但是主要受浮力作用（密度差，图中显示温差跨度约为模型的整个高度）继续向上流动和扩散，其扩散面形如羽毛。浮射流形成是因为流体初始动量较大（图中显示速度较快，0.25 m·s ⁻¹ 左右），同时也受浮力作用（图中显示温差跨度约为模型高度的一半），在动量占优势的情况下，向下流动和扩散，最终动能耗散，又形成羽流，构成一个循环。

如图1-2-16所示，双击Results→Reports→Surface Integrals，在弹出的Surface Integrals对话框中，在Report Type下拉列表框中选择Area-Weighted Average，Field Variable下方的下拉列表框分别选择“Wall Fluxes…”和“Surface Heat Transfer Coef.”，在Surfaces下方的列表框内选择bottom、left、right、top，单击Compute按钮，可得对流换热系数为5.03856 W·m ⁻² ·K ⁻¹ ，与之前的计算结果5.062 W·m ⁻² ·K ⁻¹ 相差约0.46%。

综上所述，对流换热系数的获得虽然可以采用公式计算（或查表），但是对于复杂形貌、不能简单套用公式计算的流场区域，CFD模块可以很方便地计算自然对流状态的热传递过程并求出对流换热系数。但是该系数表现为区域内均值，不能反映区域内细节，即便FEM热计算时可将其作为输入已知边界条件，仍然对区域内的热分布细节估计不足，存在一定的系统误差。 t/YHhNqhChkUSyJopl7q8Mz8XglDDwKw2HANSyQZeDRDHPMsk9xhm3BmBHJ9g32N