李新磊 吴 坤 赵林英 范学军,

* (中国科学院力学研究所高温气体动力学国家重点实验室,北京 100190)

† (中国科学院大学工程科学学院,北京 100049)

** (西北工业大学机电学院,西安 710072)

†† (合肥中科重明科技有限公司,合肥 230601)

引言

超燃冲压发动机作为一种高超声速巡航飞行器,受气动加热和燃料的燃烧释热影响因而承受着巨大的热载荷[1-3].在发动机的燃烧室壁中,通常嵌入矩形直通道构成再生冷却系统,以保证发动机能够长时间稳定运行.然而,随着飞行马赫数的提升,传统的再生冷却结构面临着冷却能力不足、结构超温等问题[4-6].因此,亟需开展新型热防护方案的设计及优化.

国内外学者对此进行了大量研究,其中最为常见的一类做法是添置粗糙元件.如Sunden 等[7]通过在圆管内设置非对称的波浪形翅片,使得管道的换热系数显著提高,壁面温度显著降低.Li 等[8]对带有截断肋的矩形通道进行了研究,分析了截断肋周围的强化传热机制,并探究了肋高和肋间距等参数对通道流动换热性能的影响规律.杨泽南等[9]对由Kagome 型网格阵列填充的通道进行了研究,发现在相同工况下,错排网格填充的通道,不仅减重效果明显,网格后复杂的涡系结构也使得流体努塞尔数明显提高,壁面温度显著下降.此外,部分研究还通过改变再生冷却通道的流动模式来达到增强换热的目的.如Li 等[10]基于变分法生成了能够自适应几何和热边界条件的变管径冷却通道,在牺牲部分压降的条件下使得壁面平均温度和温度不均匀度均有所下降.Li 等[11]对嵌套式冷却通道进行了研究,发现同向平行通道有着更好的传热性能,且内通道的嵌入有助于克服热分层现象.Zhang 等[12]提出一种双向流动模式,该模式使得壁面温度分布更加均匀,正癸烷的转化率提升以及较大的化学热沉,但壁面最高温度和压降均有不同程度的增加.

除了以上研究,近年来以拓扑优化为代表的新兴优化方法,逐渐成为散热系统设计的新途径[13-16].该方法是以微结构的材料性质作为设计变量,在固定的设计空间内自动寻找最优的材料分布方案[17],因此在理论上拥有最高的设计自由度,被认为是最具前途的优化方法之一[18].Tang 等[14]曾对考虑材料变物性的导热热沉结构进行拓扑优化,得到了如同树枝枝干般的热沉结构,同时发现固体材料的变物性显著改变了拓扑优化构型,进而影响了其在真实热环境中的性能表现.Zhang 等[19]基于伪密度法对冲压发动机的再生冷却结构进行了拓扑优化,发现其设计域内生成了众多垂直于主流方向的微通道,显著改善了整体换热,但通道内局部的逆流和滞留区导致出现了局部高温.本文认为,针对发动机再生冷却通道的拓扑优化是结合了湍流流动、燃料物性剧烈变化等特征的复杂流热耦合问题[19],其计算量大且极易出现网格依赖性布局、棋盘格单元等构型缺陷,进而影响结构性能.针对上述问题,本文将开发一套考虑变热物性的流热耦合拓扑优化求解器,采用建表插值法求解热物性,用以缓解热物性求解计算量大的问题.此外,将完整再生冷却设计域简化为一个换热设计单元,通过适当的滤波技术和松弛求解策略,以尽量避免出现构型缺陷.

本文主要内容如下,首先建立拓扑优化模型,基于连续伴随法对考虑变物性的伴随方程和灵敏度进行推导,并基于开源CFD 平台OpenFOAM[20]构建拓扑优化求解器.针对流热耦合问题开展拓扑优化设计,重点考察能量耗散约束对优化构型的影响.为了评估拓扑优化构型的流动换热性能,本文还将对拓扑优化构型轮廓进行提取,开展三维共轭换热数值模拟分析.

1 拓扑优化模型

拓扑优化方法种类众多,包括伪密度法、水平集法、相场法以及进化结构优化法等[18].在各方法中,以伪密度法的应用最为广泛和成熟.因此,本节将基于伪密度法建立拓扑优化模型,并利用连续伴随法对考虑变物性的伴随方程及灵敏度进行推导.

1.1 控制方程

针对流热耦合的拓扑优化问题,其物理模型示意图如图1 所示,设计域承受来自壁面的热载荷q或由内部产生的体积热Q,流体在流经该设计域时与固体实现热交换,进而达到热平衡.该优化问题的核心是在满足压降或材料体积比等约束下,找到最优的流固材料单元的空间分布,从而使区域的平均温度或最高温度等目标降到最低.

图1 拓扑优化设计示意图Fig.1 Schematic of the topology optimization design

首先建立流热耦合拓扑优化模型,该模型的控制方程包含如下的质量方程、动量方程和能量方程

其中,R(w,γ) 表示等式约束向量,γ 为设计变量,取值范围为0～1,w=(p,u,T) 为隐式依赖于设计变量的状态向量,在该问题中分别对应压力、速度和温度;ρ,Cp,µ,λ 和h分别为密度、定压比热容、动力黏度、热导率和比焓,均考虑其温度的相关性.τ为黏性应力张量,其定义式为τ=2µeffS(u)=(µ+µt)·为有效动力黏度,表示分子动力黏度µ与湍流黏度µt之和,Prt为湍流普朗特数,默认值为0.85,Q为体积热源.

与常规形式N-S 方程不同的是,拓扑优化模型的动量方程中引入了体积力项 -α(γ)u,α (γ) 是与多孔介质渗透率的倒数相关的插值函数.引入设计变量 γ 对 α (γ) 进行插值,用以区分流体与固体材料单元,从而实现相似多孔介质流动的描述,因此该变量又称为“伪密度”.插值函数采用如下形式[21]

式中,θ 为形状参数,取值范围为0.005～0.1,该函数形状如图2 所示.αmax的取值则与达西数Da相关

图2 不同形状参数下的插值函数Fig.2 Interpolation function with different shape parameterθ

式中，l为几何特征尺寸,Da表征了多孔介质的渗透能力,其取值小于 10-5较为合适[15].

相应地,在能量方程中对热导率进行插值处理,λmin和 λmax分别为流体和固体材料的热导率.当设计变量 γ=0 时,α 取值为一个较大数 αmax,表明该处单元为渗透率较低的固体单元;当设计变量 γ=1 时,α 取值为0,表明该处为流体单元.

本文采用了比焓形式的能量方程,与定压比热容存在如下热力学关系

因此,在求解温度时,可利用如下形式的牛顿迭代法计算得到

式中,pold,Told,Cp(pold,Told) 和h(pold,Told) 分别为上个迭代步中的压力、温度、比热和比焓,hnew为利用式(1)计算得到的更新后的比焓,该迭代在相对误差小于10-4时终止.

除了上述控制方程,本文定义了如下形式的Pnorm 函数[22]作为该优化问题的目标函数

该目标有利于避免设计域内局部区域温度过高,式中n取值为30.此外,引入如下的不等式约束函数

其中,G(w,γ) 为不等式约束向量,分别用来约束设计域的材料体积比例以及进出口能量耗散.Vmax为材料体积比例的期望值,φmax为进出口能量耗散的期望值.

综上所述,一个完整的基于密度法的拓扑优化问题可以由以下数学表达式来概述.

1.2 伴随方程与灵敏度推导

本节将基于连续伴随法对流热耦合问题的伴随关系式及灵敏度公式进行推导,首先构造如下形式的拉格朗日函数

式中,Ψ 为目标函数,λ=(pa,ua,Ta) 为拉格朗日乘子,又称伴随向量.由于R(w,γ)=0,因此原目标函数关于设计变量的灵敏度等价于如下拉格朗日函数的全导数

依据偏微分方程约束的KKT (Karush–Kuhn–Tucker)条件[23],构造如下的伴随变量方程组

利用积分变换对上述方程组进行化简,便可得到求解伴随变量所需的控制方程,其具体形式及推导过程参见附录A.能够发现,由于考虑了温度相关的变热物性和输运性质,伴随变量的控制方程形式相比于常物性模型[15]复杂很多.而伴随变量控制方程与原始状态变量的控制方程存在形式上的相似性,因此可采用相似的数值算法对其进行求解.

此时拉格朗日函数的全导数,即目标函数的灵敏度,经推导后转化为如下形式

在求解得到原始状态变量与伴随变量解后,便可计算得到该灵敏度数值.

2 数值方法

本节将对拓扑优化过程中遇到的关键数值问题以及所采取的数值方法进行介绍,包括滤波与投影函数的使用、热物性的求解、拓扑优化求解器的构建以及共轭换热数值求解器的验证等.

2.1 滤波与投影

在拓扑优化模型中,由于对材料物性采用了特殊的惩罚策略,导致优化过程中容易出现一些不切实际的布局,如网格依赖性布局或棋盘格单元等[24-26],这严重影响了材料布局的重构和收敛.研究表明,一些正则化技术,如密度滤波器[27]可以有效地解决上述问题.本文采用亥姆霍兹型偏微分方程作为滤波器,其方程形式如下

式中 γ0为滤波前的伪密度,γf为滤波后的伪密度,Rmin为滤波半径.当网格单元尺寸为 ∆x时,滤波半径Rmin越大,则 γf的形状特征越容易被抹平,因此,本文控制各算例的最大滤波半径不超过 2 ∆x.

此外,采用滤波器通常会产生灰色单元,即处于流固材料单元之间的模糊单元.因此,可以配合使用投影函数[28]来获得清晰的流固边界.本文采用如下的双曲正切函数来对 γf进行投影

式中 γp为投影后的伪密度,β 和 η 为用来控制投影函数形状的形状因子.投影函数的 η 取值为Vmax,用来保证材料的体积比例,而控制斜率的 β 因子随优化过程进行而逐渐变大,本文中 β 初始值为1,最大值为32,投影后的优化变量 γp近乎收敛到接近于{0,1}的集合,可以得到较为清晰的流固边界.

在对设计变量 γ 进行滤波和投影操作后,此时的灵敏度形式为

为了获得原始灵敏度,需要进行如下的链式求导[15]

2.2 热物性求解

在发动机冷却通道内,碳氢燃料的热物性随温度与压力在很大的范围内变化[29].以正癸烷为例,作为一种典型的碳氢燃料,在再生冷却的流动吸热过程中,其自身温升高达数百摄氏度,热力学性质和输运性质会在临界点附近发生陡峭变化,容易引发方程求解过程中的数值积分刚性问题.此外,在利用连续伴随方法推导得到的伴随方程中,存在热力学性质及输运性质关于温度的偏导项,包括 ∂h/∂T,∂ρ/∂T,∂Cp/∂T,∂ µ/∂T和 ∂ λ/∂T,当采用解析形式的物性求解方程时,偏导项的求解将变得十分困难.

因此,本文对物性参数及各偏导项的计算采取了建表-插值[30]的方法,首先基于CoolProp 软件[31]建立数据表,在不考虑燃料高温裂解的前提下,该物性数据表以压力和温度作为基础变量,变量间隔分别为1 kPa 和2 K,模拟过程中采用线性插值计算得到相关的参数.图3 所示为正癸烷在3 MPa 压力下(Tpc=648.4 K)的密度、比热以及二者关于温度的热力学偏导性质,其中利用建表-插值得到的数据与利用CoolProp 直接计算得到的数据十分接近,最大相对误差均不超过0.2%,表明本文所采用的物性计算方法是准确的.

图3 热物理性质及热力学偏导Fig.3 Thermo-physical properties and derivatives

2.3 求解器构建

基于上述的模型与数值方法,本文利用开源平台OpenFOAM 搭建了拓扑优化求解器.图4 所示为求解器的主要求解框架,其主程序包含以下4 部分:(1)对原始状态变量方程求解;(2)对伴随变量方程求解;(3)预估目标函数与约束值;(4)对灵敏度和伪密度进行滤波和投影.其中对原始状态变量和伴随变量控制方程的求解模块均基于内置的标准可压缩流动求解器rhoSimpleFoam 进行改造.对热物性及热力学偏导的求解模块则按上节所述的建表-插值法编译为动态链接库.耦合了MMA (method of moving asymptotes)[32]算法作为优化工具,同样将其编译成库供主程序调用,对设计变量进行更新.求解器迭代优化的终止条件以目标函数的相对变化小于0.1%或固定迭代步数为准则.

图4 考虑变物性的拓扑优化求解流程图Fig.4 Flowchart of the topology optimization approach considering variable physical properties

2.4 共轭换热数值验证

本文所开展的拓扑优化设计均基于二维模型,为了对其实际三维构型的流动换热性能进行评估,将对拓扑优化构型进行提取,并利用OpenFOAM 内置的标准共轭换热求解器chtMultiRegionFoam 进行数值模拟.为了验证该共轭换热求解器的准确性,首先与Zhu 等[33]进行的流动换热实验进行对比分析.该实验中固体材料为不锈钢,燃料为正癸烷.按照对应的实验条件,入口给定质量流量0.6083 g/s,入口温度为625.93 K,出口给定背压4.19 MPa,测试段体积热源大小为1.758×108W/m3,连接段热流边界的热流密度大小为11550 W/m2.采用k-ϵ 湍流模型进行验证,其结果如图5 所示,显然,该数值计算结果与实验结果较为吻合,外壁温及冷却剂温度的最大相对误差均小于5%,证明该共轭换热求解器可以较好地预测流动换热性能.

图5 冷却剂及外壁面温度的数值与实验结果对比Fig.5 Comparison of predicted coolant and out wall temperatures against the experimental data

3 结果与分析

3.1 流热耦合拓扑优化设计

本节将针对发动机再生冷却通道开展流热耦合拓扑优化设计,为了节省计算成本,以图6 所示的二维换热单元作为设计域,即不考虑构型在z方向的优化.其几何及边界条件设置见图6,计算域进出口尺寸d=5 mm,厚度为2 mm,进、出口段为长度10d的非设计域,中间设计域长度为36d,流体材料为正癸烷,固体材料为GH3128,其热导率随温度的变化如下

图6 流-热耦合结构计算域示意图Fig.6 Schematic of the computational domain of the coupled thermal–fluid structure

设定入口速度为1 m/s,初始温度为300 K,对应的入口雷诺数为2500,计算域出口压力设为3 MPa,其余边界为对称边界,对设计域施加体积热源Q=5×108W/m3.利用均匀结构化网格对区域进行离散化,考虑到要尽可能保留小尺度特征,同时还要保证加工可制造性,本文将网格尺寸设定为0.1 mm×0.1 mm,网格总数为140000,采用k-ε湍流模型进行模拟.

该优化模型以P-norm 函数 Ψ 作为目标,施加材料体积约束gvol,Vmax取值为0.68,αmax取初始值为103,每步增长1.05 倍,最大值设为108.为了考察不同能量耗散约束值对拓扑优化构型的影响,本文首先对 γ0=0.5,αmax=103的初始分布构型进行计算,得到其能量耗散值约为0.015 k g·m/s2,随后将该数值作为参考值,对 φmax的取值范围在0.02～0.25 内的算例进行了优化,并从中提取出 φmax为0.03,0.05,0.14,0.19 和0.21 的5 个算例,分别命名为Case 1～Case 5.

图7 为优化得到的伪密度场分布及对应的速度分布与温度分布.为了方便观察,这里将计算结果按关于x轴对称形式显示.能够发现,在材料比例基本一致的条件下,Case 1～Case 5 的固体域均发生了不同程度的分裂,且各分裂的固体胞元块呈现出不同的排列方式.在Case 1 中,大部分固体域分布在通道两侧,流体流通面积较大,流速较低.而随着能量耗散约束值增大,固体域分裂成的块数逐渐增多,流体域内的微细通道数量增多,流体流通面积减小,流速增大,与固体域的接触面积明显增加.从Case 1～Case 5,流体的流动分离与再混合行为逐渐增多,多处区域存在局部的流动加速,同样促进了流体与固体接触面的热量交换.

图7 拓扑优化构型及状态变量云图Fig.7 Contours of the optimized layouts and state variables

图7 拓扑优化构型及状态变量云图 (续)Fig.7 Contours of the optimized layouts and state variables (continued)

通过观察各伪密度法场对应的温度云图能够发现,尽管采用了投影等数值策略,在优化构型中仍存在部分灰色单元区域,如在Case 5 的后段,由于灰色单元的存在,导致该区域中出现流动死区,其温度呈现出非连续变化,这在一定程度上干扰了对流动换热性能的评估.因此,为了评估真实的拓扑优化通道的流动换热性能,本文将进一步对拓扑优化构型进行提取,将其对应的三维结构进行共轭换热数值模拟.

3.2 拓扑冷却通道的流动换热性能

本节将对三维拓扑优化通道开展共轭换热数值模拟,并对其流动换热性能进行分析.首先将5 种拓扑优化构型中 γ <0.5 的区域进行提取,利用样条曲线对拓扑胞元进行建模,图8 所示为5 种拓扑通道(Case 1～Case 5)的三维模型,此外,加入了传统直通道构型Case 0 作为对比算例,其固体壁肋宽约为1.70 mm,高度为2 mm,冷热壁厚度均为0.75 mm.

图8 拓扑优化通道的三维几何视图Fig.8 Three-dimensional layouts of topology optimization channels

随后对上述模型进行计算,其入口速度仍为1 m/s,初始温度为300 K,出口压力为3 MPa.按能量守恒原则,将体积热源Q换算为在固体域底部面施加的面热流q,其大小约为1 MW/m2,上下层固体域左右面为对称边界,其余壁面为绝热边界,采用kε湍流模型并结合壁面函数,对通道内各固体壁面添加边界层网格,第一层网格高度为0.004 mm,以保证无量纲距离y+<10.

图9 所示为计算得到的各通道中心截面和底部被加热面的温度分布云图.能够发现,由于固体胞元的不规则排列,5 个通道的热壁面温度呈现非均匀分布,在某些较大的胞元间隙位置,如Case 3 的190 mm 和210 mm 位置处,由于流体流通面积较大,流速较低,导致底面加热位置的壁温较高,但从Case 1到Case 5,各通道内的高温区域显著减少.

图9 xy 截面的温度分布云图Fig.9 Contours of temperature on xy plane

图10 为中心xy截面上的速度云图、压力云图以及湿周底部面的努塞尔数云图.在Case 1～Case 5中,冷却剂平均流速逐渐增大,压力损失也逐渐增大.在Case 3～Case 5 中,通道前端位置的固体胞元首先起到了节流的作用,使冷却剂流速加快,随后在流动过程中发生了大量的流动分离和再混合,带来了较大的能量耗散和压力损失.但由于冷却剂流动边界层和热边界层的分离和再附着,使得传热热阻降低,增强了流体与固体间的换热,局部努塞尔数增大.此外,在多个固体胞元的钝前缘位置,局部努塞尔数也有明显提高,这表明流体对固体的冲击效应,同样使得胞元前缘位置的传热增强.

表1 所示为各拓扑冷却通道的流动和换热性能数据,图11 为热壁面平均温度和最高温度与通道压降的性能关系曲线.能够发现,从Case 1～Case 5,各通道的压降逐渐增加,且热壁面平均温度和最高温度均随压降升高而显著降低,其中Case 5 热壁面平均温度和最高温度分别降低了91.1 K 和55.1 K.各通道的平均努塞尔数逐渐增加,表明各通道整体的换热能力逐渐增强.但相比于Case 0,仅Case 3～Case 5 起到了强化换热的效果,其平均努塞尔数增益百分比分别为12.6%,16.0%和23.4%.

表1 6 种通道的流动换热性能Table 1 Flow and heat transfer performances of six channels

图11 各通道热壁面平均温度和最高温度随压降的性能关系曲线Fig.11 Relationship between the average and maximum temperature in terms of the pressure drop of the cooling channels

3.3 拓扑冷却通道的强化换热机理

图12 所示为各通道在流向不同位置处yz截面上的涡量云图.能够观察到,相比于Case 0～Case 2,Case 3～Case 5 通道内冷却剂在多个位置呈现二次涡形态,如在Case 4 和Case 5 的x=70～85 mm、Case 5 的105～115 mm 等位置.该涡系结构的形成是由于流体在流经固体胞元间隙时,受其他分支的流体掺混与扰动而形成,其涡系形态一直保持到胞元后缘,但沿流向涡强度逐渐减弱,直到下一个固体域分裂的区域,流体再次形成明显的二次涡结构.该涡系结构促使近壁面流体与主流发生掺混,促进热交换.

图12 不同流向位置yz 截面上的涡量云图Fig.12 Contours of vorticity on yz planes at different streamwise locations

图13 所示为各拓扑冷却通道在xy中心截面上冷却剂的湍动能分布云图.在Case 3～Case 5 通道内,多个固体胞元的侧壁以及尾部回流区等位置,由于局部的流动加速或流动混合,使流体的湍动能得以激发,换热能力增强.此外,流体湍动能增强区域与二次涡形成区域基本吻合,表明二次涡结构有助于激发湍动能,进而增强换热.

图13 xy 截面湍动能云图Fig.13 Contours of turbulence energy on xy plane

4 结论

本文针对发动机的再生冷却结构,开展了考虑变物性的流热耦合拓扑优化设计,并对三维拓扑结构的流动换热性能进行了数值模拟与分析,得到主要结论如下.

(1)在基于连续伴随法的求解框架中,考虑工质的变物性会使伴随方程的形式更加复杂,求解更加困难.利用建表-插值法对方程中热物性、输运性质及各偏导项进行求解,是一种有效且准确的方法.

(2)通过对流热耦合结构进行拓扑优化设计,发现其优化构型随着能量耗散约束值的增大而愈加复杂,且伴有多重的流动分离和再混合现象.

(3)通过对三维拓扑优化构型进行共轭换热数值模拟,结果表明,相比于直通道构型Case 0,拓扑优化冷却通道内的冷却剂流态复杂,存在二次涡等涡系结构,在带来更大压力损失的同时,各通道局部的换热效果得到增强.在5 个拓扑通道中(Case 1～Case 5),Case 3～Case 5 整体实现了强化换热的效果,其平均努塞尔数增益百分比分别为12.6%,16.0%和23.4%.

综上所述,面向发动机再生冷却的拓扑优化设计具有良好的应用前景,后续将继续开展更多三维流热耦合结构的设计工作,并逐步完善求解器的其他功能.

附录A

如1.2 节所述,伴随变量方程组满足正文关系式(12),在对其进行化简前,首先将目标函数的定义域分为边界和内部域两部分

则目标函数关于状态变量的偏导数为

式中wd表示目标关于状态变量的导数的方向,变量p,u和T对应的导数方向分别为pd,ud和Td. 下面对方程组(12)中的各积分项进行积分变换,结果如下

此时,拉格朗日函数关于状态变量的偏导数可以推导为如下形式

若要满足正文式(12),可得到伴随变量控制方程

以及对应的边界条件方程

对于方程(A12)～式(A14),本文将其分别作为伴随变量ua,pa和Ta的显式边界条件,其数值随优化迭代而更新.方程中所涉及的目标函数关于状态变量的各项偏导数表达式如表A1所示.此时化简后的灵敏度表达式如下所示

对于能量耗散约束函数的灵敏度,其推导过程与换热目标函数灵敏度相似,本文不予以赘述.

附表1 目标及约束函数关于状态变量的偏导数Table 1 Derivatives of the objectives and constraints w.r.t state variables