摘要
随着制造业数字化时代到来，其软件的重要性日趋显现，并起着举足轻重的作用。国产CAE领域快速发展，作者及团队也围绕发动机多物理场耦合问题开展了一些工作，采用非结构化网格有限体积方法进行了多场耦合计算软件开发及相关应用研究。从前处理网格读取、格式转换，求解器开发各关键技术进行介绍，并重点介绍几个相关应用问题，展示该软件的应用能力和适用范围。
关键词
多物理场耦合；GTEA软件；有限体积方法；非结构化网格。
引言
随着制造业数字化时代到来，其软件的重要性日趋显现，并起着举足轻重的作用，计算科学发展成为影响国家利益与国家安全的战略性问题。一些国家将计算机建模与仿真列为优先发展的服务于国家利益的关键技术。有国内学者指出应将自主CAE软件的开发提高到战略发展高度。国内外开发了大量的通用和专用的CAE软件，哈尔滨工程大学自2003年起，启动CFD软件开发，最初研究目标是针对航空发动机燃烧室内的流动、传热与燃烧过程分析，柴油机缸内工作过程模拟，随着工作的深入和发展需要，发展了船舶水动力学分析模块。后续在结构声辐射、水动力噪声以及结构声耦合方法和软件模块开发方面进行初步工作。
网格前处理与转换技术
对于数值计算程序开发人员来说，划分求解域网格，得到求解域网格信息，例如：网格坐标，网格组成网格面系列等。这是一项具有重复性、复杂性的工作。网格生成技术及软件开发仍然存在挑战，国内相关单位开展了卓有成效的工作，但尚未形成通用、稳定和可靠的网格生成软件，也是当前CAE发展的短板之一。有经验的研究人员能够自己写出各种网格信息，但只限于几何结构简单，网格数量较少的场合；对于专业人员进行CFD求解器开发的研究人员，自己写出复杂问题的非结构化网格信息并保证不出错误是有一定难度的。因此采用一种折中的方法，团队最初开发了CGNS(CFD General Notation System)标准格式文件的输入输出，自编程序读取CGNS格式文件的网格信息，处理输出，应用于本课题组自行开发的求解器，即提供一个标准文件接口。可以采用通用的网格生成器，如ICEM CFD生成CGNS格式网格文件。为了进一步支持主流的CFD软件的前处理格式文件，提高软件的通用性，先后开发了前处理软件GAMBIT生成mesh格式文件和Star-CCM+软件生成网格文件。
CGNS 系统是一种良好的标准，核心求解器研发与商用前处理和后处理软件联系起来。核心求解器开发人员采用CGNS系统作为标准接口，则可以解除网格生成和计算结果后处理问题。作者也采用了CGNS标准作为输入输出格式，读取非结构网格和混合网格，包括二维三角形或四边形网格，三维四面体、六面体、金字塔形网格。
随着网格自适应和直角网格技术发展，如商用软件Star-CCM+和开源软件Openfoam等广泛应用，对多边形和多面体网格前后处理的需求越来越大。作者开发的软件采用了非结构化网格有限体积法，而且基于面循环的方法形成系数矩阵，该方法适用于多边形和多面体网格。问题在于目前的多边形网格和多面体网格生成器难以转化为CGNS标准格式，作者开发了直接读取Openfoam网格文件的接口程序，新版tecplot在原有基于单元的结果输出基础上，新增了基于面输出方法，包括多边形(FEPOLYGON)和多面体单元(FEPOLYHEDRON)。
GTEA求解器算法
3.1 基本方程
作者开发GTEA软件前提是连续介质假设，连续介质力学思想下统一流体力学、固体力学、声波动控制方程，详细过程可以参考相关文献。软件分为计算流体力学、结构应力波求解、流声耦合和声固耦合四大功能模块，采用非结构化网格离散控制方程，基本思想如图1所示。
3.2 并行计算流程
GTEA基本流程图如图2所示，在串行计算程序的基础上，进行了并行设计。与串行程序相比，并行计算的增加分为：区域分解与结果重建两个部分。程序开始运行后，读取网格，根据计算环境判断是否存在多个进程，存在多个进程时在主进程中根据进程数目划分网格并映射到各进程，读取设置参数后进行离散过程形成线性方程组，利用并行解法进行计算，方程组循环收敛后，对变量通过通信后进行更新，为下一步计算准备。在计算过程收敛后输出结果，对于存在多进程时，会增加输出子区计算结果，并在主进程中调用结果重建的环节。程序线性方程并行计算过程中，各相邻区域之间通过MPI进行通信，线性方程组求解过程需要全局通信。
图2 GTEA程序计算基本流程图
正如前面提及存在多进程计算时，需要将网格根据进程数目划分成相应的子区。本文调用了Metis软件系列的partdmesh工具在主进程区域分解部分，Metis软件包含了不同的工具，如partnmesh和partdmesh，对不同的进程数二者在分区时效率不同。为了支持任意混合网格，作者提出了先将网格文件转化成图形文件，然后调用kmetis划分图形。两者的区别是：前者以网格节点为节点转化成图形文件，而后者则以单元为节点。选择partdmesh是因为它可以使各相邻区域间通信两最小，而且各进程间负载较partnmesh平衡。采用了基于图论方法划分网格的partnmesh工具。
区域分解模块工作过程中，程序会通过METIS工具根据进程数将计算域分解，采用了主从模式，区域分解模块在主程序调用，对区域内的每个单元进行标识所在的进程标识，据此形成子区内部单元与ghost单元。因为子区域网格的构建依赖于通信模式的选择，在介绍子区域网格构造之前，先对通信模式做一下说明。通信模式主要有两种，以交界面和相邻单元为基础的交换模式。
a分解前区域
b通信模式1
c 通信模式2
图3 并行计算区域分解和通信模式示意图
仿真算例
4.1 自然对流与辐射传热耦合
在计算流体力学标准算例测试基础上，对软件在流动与辐射传热方面的应用进行了测试。主要算例为自然对流与辐射耦合分析。重点分析倾斜角a对流体流动和传热的影响。参数设置为Ra数为5×106，壁面发射率为e1,2,3,4 = 1，四种情况：t = 0 和t = 1, w = 0, 0.5, 1。腔体的长宽比AR是常数1。倾斜角以每步30度从-60度变化到90度。结果见图4。
图4a, b描述了倾斜角为负时，等温线随着倾斜角和散射反照率的改变的变化，正角度部分在图4d、f中描述。当t = 0时，无论倾斜角a如何变化，方腔中心部分的等温线都垂直于重力方向。对于负角度，等温壁面附近的温度梯度相比于正角度要大。从图4 a, b可看出，w=0和0.5时的等温线结构非常相似；而且，当w从0.5变化到1时，介质的散射明显改变了方腔内的热分布。然而，对于正倾斜角度，方腔内存在热分层并且中心区域的等温线是水平的。实际上，这个现象意味着少量的热传递。从图4 d, f可看出，散射反照率对方腔中等温线的结构略有影响。
(a) (b)
(c) (d)
(e) (f)
图 4 对t = 0和t = 1, w = 0, 0.5, 1，倾斜角对等温线的影响：
(1) t = 0; (2) t = 1, w = 0; (3) t = 1, w = 0.5; (4) t = 1, w = 1
4.2 发动机缸内燃油喷射过程仿真
无论是航空发动机还是柴油机，燃油雾化过程模拟是发动机流动、燃烧和传热模拟的重要环节之一，下面给出了燃油雾化过程与实验对比结果，验证软件在燃油雾化模型的准确性。
图5 不同液滴破碎模型下，燃油喷射模拟结果
4.3柴油机缸内工作过程模拟
为进一步说明，图6所示为层动区域及滑移边界的应用示意图。红线为静止的不匹配网格边界，主要有两个，一个是气道区域与气阀的区域的衔接处，一个是燃烧室区域与气缸区域的衔接处；蓝线为滑移边界，主要在气阀内部区域及外部区域的交接处以及气阀外部区域与气缸区域的交界处。彩色区域除区域2外均为层动区域，其中区域2为跟随气阀运动区域，区域1为适应阀杆运动的层动区域，区域3为适应气阀上表面运动的层动区域，区域5为适应活塞运动的气缸层动区域，区域4为适应气阀底部及活塞运动的层动区域。计算过程不同时刻的网格和流场如图7和图8所示。
图6层动区域及滑移边界的应用示意图
图7 动网格运动示意图
图8 缸内流动迹线图
4.4 燃气轮机燃烧过程仿真
利用GTEA软件对燃烧室内的流动和燃烧进行了数值模拟，数值模拟斜切径向旋流器环形燃烧室内冷、热态流场，并对所得计算结果进行分析。部分的计算结果如图9，计算所得燃烧室三维整体热态流场内通过旋流器中心截面的速度矢量图。由放大图A图可知，与双级轴向旋流器燃烧室一样，在扩压器的上下突扩段有明显的旋涡，这两个旋涡作为气动壁面根据进口状态的变化可自动进行调整，保证火焰筒进口气流稳定，旋流器出口高速旋转射流与主燃孔射流相互作用，在火焰筒头部产生了上下两个较强的旋涡，形成中心回流区如放大图B所示。C区放大图可以看到，主燃孔后，气流逐渐平稳，冷却孔形成局部的高速区。
图9 中心截面速度分布图
图9为计算所得中心截面气流温度场，从图中可以看到，从火焰筒气膜冷却孔流出的冷却空气对火焰筒壁面起到了很好的冷却作用，同时在主燃孔后，由于未燃尽燃气在火焰筒主燃孔后低压区继续燃烧，因此主燃孔后火焰筒壁面温度比周围的高。
4.5 水下航行体SUBOFF水动力仿真
确立的潜艇模型是以SUBOFF 潜艇模型为基础进行建立的。由美国国防预研规划署(DARPA)确立的SUBOFF 项目，SUBOFF 潜艇模型的艇体部分是用回转体进行模型确立的。关于SUBOFF 项目的提出，为了方便我们对潜艇模型进行模拟计算，分析其水动力性能和潜艇周围流场。自从1989 年美国国防预研规划署(DARPA)提出这个模型，已经掌握了大量的实验数据，为很多国内外的学者提供了研究验证的对象。
图10 SUBOFF模型流向速度速度和压力分布图
图11 SUBOFF模型中纵剖面上压力系数分布
4.6结构砰击水动力仿真
计算采用了文献中模型，如图12所示，计算域大小为3.22m×1m×1m，障碍物位于矩形水箱底部，水体前侧，大小为0.16m×0.4m×0.16m，距离布置水体侧壁面2.48m。初始水体大小为1.228m×1.0m×0.55m。在矩形水箱底部H4(0.582,0.5)设置波高仪监测自由面高度，在障碍物表面P3(2.4,0.475,0.10)设置压力传感器监测压力大小。
图12 计算模型示意图
计算模型网格取为结构化网格，网格数为161×50×50。图13为自由面高度监测点H4(0.582,0.5)处自由面高度随时间变化曲线，图14压力监测点P3(2.4,0.475,0.10)处的压力随时间变化曲线，可以看到计算结果和实验结果吻合良好。
图13 H4监测点处自由面高度随时间变化曲线
图14 P3处的压力随时间变化曲线
4.7复杂消声器声学性能仿真
如图15该消声器的出口管改到侧面，则该消声器结构的进出口边界由轴对称变为只有一个对称面，这就导致计算网格会成倍的增加，同时为测试FVM方法在大型计算中的能力，划分网格数为778579，利用8进程进行仿真。对于在本节提到的计算机上，FEM是难以完成的，而FVM方法由于采用显格式推进，并且不需要计算大型的刚度矩阵，所以可以比较快捷的处理较大型计算问题。为出口管位置不同时消声器的传递损失对比，可以看出，当出口改为侧面后在共振频率处出现多个峰值：
第一个共振峰是由于侧面出口管与消声器端部之间形成一个声学共振腔引起的；
后面的几个共振峰是由于侧面出口使消声器变为非对称结构而激发周向模态，消声器的高阶模态提前激起而形成的，出口改到侧面后消声器的整体声学性能得到改善，特别是1500Hz以下的声学性能得到明显提高。
图15 复杂消声器网格图与消声器传递损失比较
4.8 流场动力噪声
计算模型和网格划分见图16，采用相同的模型和网格，在圆柱壁面附近进行加密。圆柱直径为D=0.019m，流场马赫数Ma=0.2，雷诺数Re=200，介质为空气。时间步长为dt=2×10-5s，总计算时间为t=0.15s，当不可压流场计算趋于周期变化时，即当t=0.05 s时，开始计算声场。进行流场计算时，左边界为速度进口，右边界为压力出口，上下为对称边界，域内正方形为内部面，在声场计算中，域内正方形到外边界为PML区域。监测点设于坐标为A(0 ，10D)、B(0 ，35D)、C(0，44.9D)和D(0 ，45.1D)。
图16 计算模型和网格划分
图17 瞬时声压云图和y轴的声压分布
图17(a)为计算区域声压分布云图，从图中可以看出，圆柱绕流噪声是个偶极子声源，PML区域吸收良好，图17(b)是y轴的声压分布，可以看出圆柱上下方声压幅值相反。
4.9结构声耦合算例
考虑如图18所示的大坝-水库结构声耦合系统。大坝的介质属性为：弹性模量E = 3.437×109 Pa，泊松比μ = 0.25，密度ρs = 2000kg/m3，纵波波速为cp = 1436m/s；水的介质属性为：声速ca = 1436m/s，密度ρa = 1000kg/m3。水面高度H = 50m，均匀分布的正弦载荷f(t) = 200×sin(18t) N/m作用在大坝顶部。大坝底部采用固支边界条件，水库上侧水面为绝对软边界，下侧固体壁面为绝对硬边界，右侧为吸收边界。在结构子域、声学子域均采用显格式求解。
图18 大坝-水库系统和计算网格
(1) 四边形单元处理改进前后结果对比分析
网格模型如图14所示，改进前采用常应变型四边形单元模型，改进后采用双线性四边形单元求解结构声耦合问题。边界上网格尺寸为5m，与文献相同，结构子域中网格的最大边长Lmax ≈ 9.038m，最小边长为Lmin ≈ 2.744m，声学，子域中Lmax = Lmin = 5m。基于显-显耦合格式求解，依据稳定性条件，结构子域的时间步长需要满足∆ts ≤ 1.91×10-3s，声学子域的时间步长需要满足∆ta ≤3.48×10-3s，最终结构子域与声学子域取相同时间步长为1.5×10-3s。
(a)A点位移uy (b)B点声压
图19 四边形单元计算结果
图所示为A点位移uy的时间响应与B点压力的时间响应。将计算结果与文献中FEM程序的计算结果对比如图19所示，发现对于该非点源问题，改进处理前的四边形单元不能得到正确的结果，而改进后的计算结果与文献结果吻合更好，说明对四边形单元的双线性改进处理，提高了传统常应变单元方法在处理结构声耦合问题时的正确性。
4.10 结构模态分析
潜艇是一种既能在水上又能在水下航行的舰艇，根据其性能指标需求其外形有不同构形，但主体结构基本相似。对于单层壳潜艇，其耐压壳体可以看做是由圆柱壳、圆锥壳和半球壳组合组合结构，本节将采用此简化结构分析艇体主体结构振动特性，计算结果与实验值及FEM结果进行对比，FVM与FEM分析采用同一套网格。简化结构实物模型及FVM分析模型如下图所示：
图20简化结构实物模型 图21仿真模型
图22 前6阶模态结构
4.11 水下结构声固耦合仿真分析
考虑壳体弹性对声场的影响，壳体厚度为0.015m，结构为钢，介质属性为密度ρs = 7700.0kg/m3，杨氏模量E = 1.44×1011Pa，泊松比v = 0.35，假设为平面应变问题，纵波波速为cp = 5478.5m/s。壳体结构子域908个三角形单元划分，最小边长为Lmin ≈ 0.0106m，依据稳定性条件，选择的时间步长Δtso = 1.36×10-6s；声学子域采用62084个三角形单元划分，最小边长为Lmin ≈ 0.0144m，依据稳定性条件，选择的时间步长Δtao = 8.22×10-6s；最终采用不同时间步长技术∆ta = 5∆ts = 5.0×10-6s。
图23为水下弹性结构时，不同时刻声压云图分布，可以看到点声源的辐射场到达壳体后，在t = 0.0010s时已经激起结构振动并产生声辐射，t = 0.0020s后在水下结构头部可以明显看到结构振动产生的辐射声场，这主要是由于声源在尾部激起结构振动，而结构中纵波的波速远大于水中的声速，由此产生波传递的速度差，结构中纵波先到达头部，并且振动产生新的声辐射。
图23 弹性结构不同时刻声压(Pa)云图
5、 结论
文章介绍了明平剑教授团队围绕多物理场耦合问题进行计算研发的基本思想，计算软件开发的一些实际问题及解决方法，包括介绍了网格格式及并行计算问题等，基于非结构化网格有限体积方法研究可以参见明平剑教授及团队成员发表的文章。
往期回顾：
扫一扫
8852
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
