基于特征分裂有限元准隐格式的共轭传热整体耦合数值模拟方法1)

刘 瑜 邓家钰 王成恩,2) 苏红星

∗(上海交通大学机械与动力工程学院,上海 200240)

†(大连理工大学航空航天学院,辽宁大连 116024)

引言

共轭传热现象在科学和工程领域中大量存在[1],比如核反应系统[2],燃气轮机涡轮叶片冷却[3-4],航天器热防护[5-6],化工处理[7],生物医学[8-9],电子设备散热[10],MEMS 设备中的微槽道对流换热[11]等.

共轭传热问题除了实验研究以外,可以通过理论分析或数值模拟的方法求解.理论分析方法对简单的模型(如平板,圆管等)进行建模,通过简化得到精确或近似的理论解[12].理论方法对共轭传热问题的理解很重要,但只适用于简单的问题.数值模拟将流动的数值解与固体传热的数值解耦合起来,主要有两种耦合方法:分区耦合(partitioned)[13]和整体耦合(monolithic)[14]方法.分区耦合方法对流场和固体热分析采用各自相应的成熟的求解器,然后通过迭代方法满足流热耦合的界面条件,即在界面上温度和热流应是连续的.分区耦合方法的优点是可以采用现有的流动和固体热分析求解程序,但是为了满足界面条件需要进行迭代求解,在不同的求解器之间传递边界条件,算法实现复杂,并且还需要考虑分区耦合算法在流固耦合界面上的稳定性条件[15-20].整体耦合方法将流动方程和固体传热方程统一离散,然后求解单一的方程.整体耦合方法的优点是不需要处理流固耦合界面条件,缺点是如果流动特征时间和固体传热特征时间差别很大,非定常共轭传热的计算量将是巨大的[21-22].

目前存在大量的数值方法模拟共轭传热问题，如有限差分法[23],有限体积法[6,24-25],浸入式边界法[26],无网格法[27],格子玻尔兹曼方法[28-30],边界元法[31],有限元法等.有限元法在计算固体传热问题上具有优势，如果流场求解也采用有限元求解,可以降低程序实现的复杂性.采用有限元法求解共轭传热的文献相对较少.Misa 和Sarkar[32]采用整体耦合方法,利用满足inf-sup 条件的经典伽辽金(Galerkin)有限元方法对方腔自然对流的共轭传热进行了模拟.Malatip 等[33]采用Rice 和Schnipke[34]提出的交错有限元方法和整体耦合方法求解定常共轭传热问题.Malatip 等将Choi 和Yoo[35]提出的二阶时间精度分裂有限元方法推广到求解共轭传热问题[36]和流−热−结构耦合问题[37].Verstraete 和Scholl[19]以及Scholl 等[20]采用压力稳定有限元法计算流动,定常有限元法计算固体传热,对分区耦合方法计算定常共轭传热的稳定性问题进行了研究.Long 等[38]将光滑有限元方法(smoothed finite element method) 和改进的光滑粒子流体动力学方法耦合起来求解热−流−结构相互作用问题.

利用伽辽金有限元法求解不可压缩流动,需要有限单元满足所谓的inf-sup 条件,即速度和压力的插值函数不能具有相同的阶数[39].这增加了算法实现的难度,特别是对于三维问题更是如此.为了采用等阶元计算,Hughes 等[40]提出了用PSPG (Pressure-Stabilizing/Petrov-Galerkin) 稳定项来规避inf-sup 条件,为了实现对流主导流动的稳定计算还需要添加SUPG(Streamline-Upwind/Petrov-Galerkin) 稳定项[41].Zienkiewicz 等[42]在分数步法的基础上,提出了基于特征分裂的有限元法(CBS).CBS 方法有半隐格式(semi-implicit),准隐格式(quasi-implicit) 和基于人工压缩性的格式(CBS-AC)[42].准隐格式的时间步长只取决于对流稳定性,计算效率高于半隐格式和基于人工压缩性的格式[43].在作者阅读的有限文献范围内,还没有看到将CBS 方法的准隐格式用于共轭传热数值模拟的例子.CBS 方法既适用于不可压缩流动的求解,也适用于亚声速、超声速流动的求解;并且同样可以采用等阶插值函数;CBS 方法是分数步方法,每一个分裂步得到的方程是标量方程,对存储量要求小,方程的求解也相对容易.本文的主要目的是采用CBS 方法的准隐格式,发展一种整体耦合层流共轭传热数值模拟方法.

本文接下来给出共轭传热的控制方程,并介绍流场计算采用的基于特征分裂的有限元方法,随后给出整体耦合共轭传热数值模拟算法.为了验证本文所采用算法的准确性,对不可压缩流动问题和共轭传热问题进行了模拟.

1 计算方法

1.1 控制方程

非守恒形式的不可压缩Navier-Stokes 方程为

上式是无量纲化的方程,其中ui,(i=1,2,3)是流体速度分量,p是流体静压,Re为雷诺数.方程中变量下标满足爱因斯坦求和规则.

假设流体的密度是常量,且其他流体性质也与温度无关.将能量方程表示为温度的形式,即

其中,T为温度,Qf为流体中的热源项,Pr为普朗特数.在这种情形下,能量方程与连续性方程和动量方程解耦.温度只是简单的输运变量,在得到了速度场后,计算温度方程就可以得到温度场的分布.

固体传热方程为

其中κsf=κs/κf,cr=ρscps/ρfcpf.ρf,ρs分别是流体和固体的密度,κf,κs分别是流体和固体的导热系数,cps和cpf分别表示固体和流体的比热.推导以上无量纲方程采用的参考量见附录.

1.2 求解不可压缩流动的分数步法

方程(1)的半离散形式为

其中∆t是时间步长,θ1,θ2和θ3分别是控制对流项、压力项和黏性项时间离散的参数,其取值范围为[0,1].为了避免求解非线性方程,只考虑θ1=0 的情形.在推导CBS 格式之前,首先引入辅助变量∆u∗和∆u∗∗,且令

采用基于Helmholtz-Hodge 分解的2 级分数步方法,将方程分裂为

1.3 对流稳定项

对于对流占优问题,CBS 格式需要添加对流稳定项[43].本文在半离散方程(5) 中添加显式的稳定项,即

其中∆ts为对流稳定项需要的时间步长,可以认为是稳定性参数,与SUPG 中的稳定参数τ[41]的作用相似.在分数步方法中,稳定项也要进行分裂,在预测步式(7)中添加的稳定项为

即不考虑压力项.在修正步式(8) 中,稳定项仅包含压力项,即

类似可以得到带稳定项的能量方程半离散形式,即

其中最后一项是对流稳定项.

1.4 基于特征分裂有限元方法的矩阵形式

对半离散方程采用伽辽金有限元方法进行空间离散,根据参数的不同,可以得到不同的格式[43].

1.4.1 半隐格式

步骤1

其中Mu,i,C,K分别表示速度质量矩阵,离散对流算子和黏性算子.Fi表示黏性边界积分,fs,i是稳定项.

步骤2

步骤3

其中L,D,G分别表示离散拉普拉斯算子,离散散度算子和离散梯度算子,fp,i是稳定项.半隐格式是条件稳定的,稳定性由对流项和黏性项确定.

在计算得到了速度场后,即可以将速度场代入到能量方程计算温度场.

步骤4

其中,MT,KT分别是温度方程质量矩阵和温度扩散算子,FT表示温度扩散项边界积分,QT为热源项,fT为稳定项.

半隐格式的单元稳定时间步长为

其中上标e表示单元,h是单元的网格尺度,∆tcon,∆tvis分别是对流时间步长和黏性时间步长.全局时间步长取所有单元稳定时间步长的最小值,即

黏性稳定性在一定条件下会对时间步长产生严重的限制.为了克服这一缺点,提出了准隐格式.

1.4.2 准隐格式

与半隐格式不同,准隐格式对黏性项进行隐式处理[43],即1/2θ31.

步骤1

黏性项的边界积分仍然作显式处理.步骤2 和步骤3与半隐格式相同.

步骤3 还可以采用非投影形式,即对于速度修正步

步骤4

准隐格式的时间步长只受对流项限制,即

以上矩阵和源项的具体形式见附录.

说明1:半隐格式的时间步长受黏性稳定性条件影响;而基于人工可压缩性的格式(CBS-AC),根据文献[43]报道,准隐格式的计算效率比其高3 ∼100 倍.因此本文计算采用准隐格式.

说明2:CBS 方法允许等阶插值.本文选择等阶有限元,即速度、压力、温度都采用相同的有限元基函数,这大大方便了程序的实现.

说明3:准隐格式求解的四步方程的矩阵是对称正定的,故可以采用共轭梯度法求解.预处理矩阵本文采用雅克比方法计算.

说明4:本文计算所采用的参数为θ1=0,θ2=θ3=θ4=1.

1.5 共轭传热耦合算法

在进行流固耦合模拟时,需要划分流场网格和固体网格.本文设计了一种流固耦合网格和边界条件表达方式.出于简单考虑,流固耦合界面是适配的,即在流固耦合界面上流场网格点和固体网格点是一致的.为了方便,流体区域和固体区域作为一个统一的计算域划分单一的网格.为了区分流体和固体网格,需要用一个指标集来表示网格单元的材料属性,考虑用一个材料属性文件来存储网格单元的属性,即mateiral.txt.

对于温度场而言,流固耦合边界属于内面,因此在流固耦合边界上只需要指定流场变量边界条件.而在固体的其他边界上只需要指定温度边界条件,不需要指定流场变量的边界条件.

如果采用前面提及的流固耦合网格及边界条件,可以得到共轭传热耦合算法如Algorithm 1 所示.

固体传热方程离散后的形式为

温度方程在流体区域按照式(27) 离散,在固体区域按照式(29) 离散,最后将温度方程装配为统一的方程组.

该算法是一种整体耦合算法,与分区算法相比,不需要通过迭代就可以隐式满足流热耦合界面条件.

1.6 时间步长

1.6.1 流体时间步长

Bevan 等[43]研究了CBS 方法中的半隐格式、准隐格式和基于人工压缩性的格式的计算效率,发现准隐格式由于只需要考虑对流时间步长,相比于其他两种格式具有较高的计算效率.文献[43] 中没有对时间推进步长和稳定项中的时间步长进行区分,在稳定项中采用的也是时间推进步长,即∆ts=∆t,并且建议采用半隐格式和准隐格式时,时间步长应该选择全局时间步长.

作者认为稳定项中的时间步长∆ts可以视为稳定性参数,而稳定性参数应具有当地性质.当计算采用全局时间步长时,如果整个计算区域的稳定性时间步长都取同样的值,有可能不能较好地抑制数值不稳定性.因此本文采用准隐方法进行计算,对物理推进时间步长和稳定项中的时间步长做了区分.物理推进时间步长∆t按式(28)计算，取全场的最小时间步长作为全局时间步长.稳定项中的时间步长为局部时间步长∆ts，按照式(22)计算.

1.6.2 固体时间步长

在固体传热方程(4)中,固体密度与比热的乘积与流体密度与比热的乘积之比,cr,在一般情况下有cr≫1.这导致传热方程时间尺度远大于流动方程的时间尺度.He 等[44]估计了典型共轭传热问题的固体和流体的稳定时间步长之比为103∼105.当对共轭传热问题的流体温度方程和固体温度方程统一求解时,如果采用流体部分的稳定时间步长进行计算,则固体部分的温度场需要迭代很多步才能收敛.对于定常问题,为了提高收敛速度,在固体部分需要采用比流体时间步长更大的时间步长

其中∆tflu和∆tsol分别表示流动计算和固体传热计算采用的时间步长,n为放大因子.

式(30)中如果选择n=cr,则相当于在方程(4)中令cr=1,∆t=∆tflu=∆tsol.这是文献[23,45]采用的做法.本文定义名义上的cr值为cr′,即

计算时用cr′代替cr.名义上固体时间步长与流体时间步长相同,但是通过调整cr′实际上改变了固体时间步长.cr′可以足够小,在满足稳定性的前提下提高收敛速度.这将在下面的算例中得到验证.

2 算法验证

2.1 稳定项时间步长选取方法对计算结果的影响

首先对稳定项中的时间步长的选择方式对计算结果的影响进行比较.采用3 种方式计算稳定项中的时间步长,分别是:与全局时间步长相同,采用局部时间步长和取零值(即不添加稳定项).计算采用的例子是方腔顶盖驱动流.方腔的顶盖按指定的速度运动,方腔的其他边界固定.在顶盖的作用下,方腔中流体将发生运动.计算域Ω 为(0,1)×(0,1),边界条件都为Dirichlet 条件.顶盖的运动速度,按照文献[46],给定如下

当网格足够稠密时,稳定项的作用不明显,为此在这里特意选择比较稀疏的网格进行计算.网格包含800个线性三角形单元和441 个节点.网格如图1(a) 所示.图1(b)∼图1(d) 分别给出了由3 种方法计算的Re=5000 时x方向上的速度分量u的云图.图2 给出了在y=0.8 上u的分布.可见如果不加入稳定项,流场变量将出现很大的伪振荡;采用文献[43] 中的方法,由于稳定项中的时间步长全场都是相同的,不能完全抑制伪振荡的出现;而采用本文的方法流场变量分布光滑,抑制了伪振荡的出现.

图1 顶盖驱动流的计算网格和速度云图Fig.1 Computational mesh and contours of x-component of velocity

图2 由不同稳定项计算方法得到的速度的x 分量在y=0.8 上的分布Fig.2 Velocity profile at y=0.8 predicted by different stabilization schemes

2.2 层流计算验证

2.2.1 方腔顶盖驱动流

方腔顶盖驱动流是用来检验层流计算准确性的典型算例.计算和边界条件与上节相同.采用一阶三角形单元计算,单元数为3200,节点数为1681,在边界附近进行了加密.分别对Re=400,1000,5000,10 000 时的情形进行了计算.图3 比较了x=0.5 上速度分量v和y=0.5 上速度分量u的分布与Ghia等[47]的计算结果,可见符合很好.

2.2.2 层流后台阶流动

对于这个问题有实验数据可供比较,因而为很多作者所引用.本文选择雷诺数为229 的情形进行计算.问题的计算域和边界条件如图4 所示.进口速度u的分布的表达式为

在壁面处按无滑移边界条件处理,在出口速度指定Neumann 条件,压力指定为0,在其余边界上压力指定为Neumann 条件.计算网格包含15 200 个三角形单元,7841 个节点.

图5 分别给出了压力云图和速度云图.图(6)给出了不同截面上速度分布与实验数据[48]的比较,可见计算结果与实验符合很好.

图3 不同雷诺数条件下的速度型面Fig.3 Velocity profiles in the x and y directions

图4 层流后台阶流动计算域和边界条件Fig.4 Domain and boundary conditions of flow over a backward facing step

图5 层流后台阶流动的压力和速度云图Fig.5 Pressure and velocity contours of flow over a backward facing step

2.3 共轭传热数值模拟

本文对若干共轭传热问题进行模拟,以检验算法和程序的准确性.其中,反向流动热交换器在模拟中考虑了cr′的取值对收敛速度的影响,其余算例都设置cr′=1.

2.3.1 共轭传热库埃特流

共轭传热库埃特流[33]的示意图为图7.上壁以固定速度运动,温度为T2=0;下壁固定具有厚度,需要考虑热传导,下壁底部温度为T1=1;流固耦合界面上温度的精确解为

其中h1,h2分别为固体和流体区域的厚度.流体域中的温度和固体域中的温度在y方向上为线性分布.

图8 给出了取不同的κsf值时,由一阶三角形单元计算的温度沿y方向上的分布曲线,网格单元数为1200.可见计算值与理论值十分吻合.

图8 不同条件下共轭传热库埃特流温度分布曲线与理论解的比较Fig.8 Comparison of benchmark solutions for Couette flow condition

2.3.2 二维反向流动热交换器

二维反向流动热交换器的示意图如图9 所示.在热交换器中分布两个平行流道,流道被一金属平板分隔开.两个流道的厚度分别为δ1,δ3,隔板的厚度为δ2.流道的外壁是绝热的,流道入口处的速度和温度被认为是均匀的.计算采用的参数如下:几何参数δ1=δ2=δ3=0.1 m,L=1 m;流动参数U1=0.2 m/s,T1=800 K,U3=0.1 m/s,T3=300 K,ρf=1000 kg/m3,κf=10 W/(m·K),cpf=25 J/(kg·K),µ=0.15 kg(m·s);金属板的物性ρs=8000 kg/m3,κs=50 W/(m·K),cps=500 J/(kg·K).

图9 二维反向流动热交换器Fig.9 A conjugate counter flow heat exchanger

计算网格采用线性三角形单元,单元数为3360,节点数为1763.

在本例中cr的准确值为160.为了考虑cr′的取值对收敛速率的影响,分别对cr′采用准确值和cr′=1,0.001 的情形进行了模拟,图10 给出了温度残值的收敛曲线.可见cr′如果采用准确的值,收敛将非常缓慢.如果令cr′=1,收敛速度将大幅增大,继续减小cr′的值，收敛速度会继续增大.cr′的取值不影响最终的收敛解.

图10 采用不同的cr′ 的收敛速率比较Fig.10 Comparison of convergence rate for different cr′ values

计算的速度场和温度场分别如图11(a)和图11(b)所示.本文还计算了不同的固体和流体热传导系数之比条件下的共轭传热.在x=L/2 上温度的分布曲线如图12 所示,并与Malatip 等[37]计算的结果进行比较,可见与文献给出的结果符合良好.

图11 二维反向流动热交换器x 方向上的速度云图和温度云图Fig.11 Predicted x-component of velocity and temperature contours for a counter flow heat exchanger

图12 分别取κsf=0.1,1,5,10,二维反向流动热交换器的温度在x= L/2 上的分布与文献[37]的比较Fig.12 Temperature distributions at x= L/2 for κsf=0.1,1,5,10 of a counter flow heat exchanger compared with results from Ref.[37]

2.3.3 加热固体的强迫对流冷却

第3 个测试例子是流动流经两块平行平板间含有3 个矩形加热固体的强迫对流冷却问题.该问题首次由Davalath 和Bayazitoglu[23]提出用来模拟集成电路元件的冷却.该问题的计算域和边界条件如图13 所示,矩形块尺寸相同,且等间距分布.流体以完全发展的速度型面从左侧进入,从右侧流出.3 个固体块中的生成的热可以假设为具有均匀分布的热源项,本文中假设热源项Qs=8.计算时取Pr=0.71,κsf=10,分别计算了雷诺数为Re=100,500,1000 的情形.计算网格如图14 所示,网格中含有个10 900 个线性三角形单元,5708 个网格节点.图15 给出了雷诺数为100 时的压力、速度和温度云图.图16 给出了固体−流体界面上温度的分布,并与文献[37]的结果进行了比较,可见本文计算结果与文献符合很好.

图13 加热固体的强迫对流冷却问题的计算域和边界条件Fig.13 Domain and boundary conditions of forced convection cooling across rectangular blocks

图14 加热固体的强迫对流冷却问题的计算网格Fig.14 Computational mesh of forced convection cooling across rectangular blocks

图15 雷诺数为100 时,加热固体的强迫对流冷却的压力,速度幅度和温度云图Fig.15 Pressure,velocity magnitude and temperature contours for Re=100

图16 雷诺数为100,500,1000 时固体−流体界面上温度的分布与文献[37]的比较Fig.16 Comparison of wall temperature distributions along solid-fluid interface of rectangular blocks for Re=100,500,1000 with Ref.[37]

3 结论

本文通过将物理推进时间步长与稳定项中的时间步长进行区别,改进了基于特征分裂的有限元法的准隐格式,使其具有了更好的稳定性.在此基础上,利用改进了的准隐格式和整体耦合方法实现了不可压缩流体和固体间共轭传热的数值模拟并通过与理论、实验和文献中的数值解进行比较,验证了本文提出方法的正确性.还研究了固体区域时间步长对定常共轭传热问题数值模拟收敛性的影响,发现选择合适的固体时间步长可以显著提高收敛速率.

本文的一个不足之处是没有对稳定时间步长对基于特征分裂的有限元法准隐格式的稳定性的影响进行理论分析,在未来的工作中将对此进行研究;同时还可以研究物理推进时间步长采用局部时间步长对定常共轭传热问题计算收敛性的影响.未来还将采用本文提出的方法对自然对流和混合对流问题进行研究,进一步探索该方法的适用范围.

附录

控制方程的无量纲化

有量纲形式的控制方程为: