基于特征正交分解的一类瞬态非线性热传导问题的新型快速分析方法1)

朱强华 杨 恺 梁 钰 高效伟

∗(南昌航空大学飞行器工程学院,南昌 330063)

†(大连理工大学航空航天学院,工业装备结构分析国家重点实验室学,大连 116024)

引言

瞬态非线性热传导问题在实际工程中普遍存在,这是因为材料的热物性参数通常是随温度变化的.由于这类问题的复杂性,很难得到其解析解,一般采用有限差分法[1-2]、有限元法[3-5]、有限体积法[6]等传统的数值方法得到其近似解.鉴于这类问题的重要性,研究者仍不断致力于发展新的求解手段.比如,顾元宪等[7]采用精细积分法求解瞬态非线性热传导问题,推导了热传导方程精细积分的具体列式;Tang 等[8]采用半边界法(half-boundary method,HBM)求解核燃料棒的一维瞬态非线性热传导问题;Yang等[9]采用径向积分边界元法(radial integration boundary element method,RIBEM)求解导热系数随温度变化形成的瞬态非线性热传导问题;Cui 等[10]采用Gao 等[11]近年提出的单元微分法(element differential method,EDM)求解多维瞬态非线性热传导问题.然而,这些数值方法在对工程中的大型、复杂系统进行求解时,得到的离散方程规模巨大,其自由度数量可达数百万甚至千万,这一方面对计算机的性能和存储容量提出了极高的要求,另一方面需要耗费大量的计算时间,无法满足一些需要对温度进行实时控制或快速预测的领域的要求.为此,特征正交分解(proper orthogonal decomposition,POD)作为一种有效的模型降阶方法,能够将高维系统转变为低维系统,从而大大缩短计算所需要的时间,因此是解决这个问题的重要途径.

POD 方法又被称为主成分分析(principal component analysis,PCA)或奇异值分解(singular value decomposition,SVD),在流体力学[12-13]、结构动力学[14-16]、最优控制[17-18]等诸多领域都得到广泛应用.在传热学领域,Biaecki 等[19]采用POD 方法对瞬态线性热传导问题的有限元模型进行降阶,使计算效率得到显著提高;Zhang 和Xiang[20]将POD 方法和无网格法结合建立了瞬态线性热传导问题的降阶模型,在保持计算精度的同时使计算时间大为减少;Gao 等[21]提出了一种多域POD 方法,能够对由多种常物性材料组成的复杂求解域的瞬态热传导过程进行快速分析;胡金秀和高效伟[22]建立了变系数瞬态热传导问题边界元格式的POD 降阶模型,解决了边界元法计算速度慢、算法稳定性差的问题;冯俞楷等[23]提出从温度场的梯度提取POD 模态建立瞬态线性热传导问题的降阶模型,不仅提高了计算效率,还具有更强的外推能力.这些研究展示出基于POD的模型降阶方法在处理线性热传导问题时的高效性及其与各种数值方法结合的广泛适用性.

虽然POD 是依赖信号来重构系统的一种线性处理方法,但研究表明它对于非线性系统也有一定的应用价值[24-25].目前国内外对瞬态非线性热传导问题的模型降阶研究尚非常缺乏.Fic 等[26]初步研究了对瞬态非线性热传导问题采用POD 方法建立其降阶模型的可行性,发现无需重新计算POD 模态;Han 等[27]针对瞬态变物性热传导问题发展了基于适体坐标(body-fitted coordinate,BFC)和有限体积法的POD 降阶模型.尽管降阶分析能够准确地预测非线性热传导过程,但是对计算效率的提升效果却远不及其对线性热传导问题的那样显著,通常要小1∼2个数量级.其原因是在瞬态非线性热传导问题的迭代计算过程中,每个迭代步都必须更新降阶模型的系数矩阵,其关联计算会耗费大量时间.于是POD降阶模型的计算效率成为制约其在瞬态非线性热传导问题中应用的关键,也成为了近年来的研究热点.Gaonkar 和Kulkarni[28]采用两级离散和多级格式相结合的方法(two level discretization-multilevel scheme,TLD-MS)有效提高了瞬态非线性热传导问题的POD降阶模型计算效率,但该方法局限于结构化网格且存在网格匹配性要求以及插值引起额外的精度损失;Feng 等[29]采用组合近似方法(combined approximations,CA),而Ding 等[30]则采用等几何分析和独立系数相结合的方法(isogeometric analysis-independent coefficients,IGA-IC),分别建立了瞬态非线性热传导问题的非POD 降阶模型,虽然提高了计算效率,但在全阶模型的自由度数量较少时效果亦不明显.

在前期的研究中[31],我们发现可以采用瞬态线性热传导分析中得到的POD 模态对相同求解域的瞬态非线性热传导问题建立降阶模型进行近似计算,这虽然简化并加快了降阶模型的构建,但并未涉及到如何解决非线性降阶模型计算耗时的关键问题.本文对此展开了深入系统的高效算法研究,结构安排如下:第二部分给出瞬态非线性热传导问题的有限元全阶模型及其计算方法;第三部分先采用POD 方法建立瞬态非线性热传导问题的降阶模型,并简要介绍其常规计算方法,然后重点阐述新型计算方法的基本理论;第四部分通过二维和三维算例验证本文提出的瞬态非线性热传导降阶模型新型计算方法的准确性和高效性;最后进行总结.

1 瞬态非线性热传导方程的有限元法

1.1 瞬态非线性热传导问题的控制方程

考虑材料的导热系数随温度变化形成的非线性热传导问题,在求解域Ω 内其非稳态、无内热源时的控制方程可表示为

式中,ρ 为材料的密度(kg/m3),c为材料的比质量热容(J/(kg·K)),T(x,t)表示t时刻空间坐标x处的温度(K),k(T)表示材料随温度T变化的导热系数(W/(m·K)),t0为初始时刻(s).

方程(1)是一类抛物型偏微分方程,其定解条件包括初值条件

和边值条件

式中,T0为初始时刻的温度场,Γ=Γ1+Γ2+Γ3为求解域Ω 的边界,其中Γ1是第一类边界(Dirichlet 条件),是在边界Γ1上给定的温度,Γ2是第二类边界(Neumann 条件),是在边界Γ2上给定的热流密度(W/m2),Γ3是第三类边界(Robin 条件),是在边界Γ3上给定的对流换热系数(W/(m2·K)),Ta为环境温度,ni(i=1,2,3)是边界Γ 外法线的方向余弦.当和不随时间变化时称为定常边界条件,否则称为时变边界条件.

这里需要说明的是,由于本文提出的模型降阶快速分析方法目前仅适用于解决第二类和第三类边界条件的瞬态非线性热传导问题,因此本文暂时不涉及第一类边界条件的相关问题.

1.2 瞬态非线性热传导方程的有限元离散格式

应用有限元法求解上述瞬态非线性热传导问题,方程(1)及其定解条件(2)和(3)的弱解形式为

式中,N为权函数.由于所构造的权函数需要满足正交的特性,这在整个求解域内是难以实现的.为此有限元法通过将求解域离散成若干个微元子域(即网格单元),在各子域内构造相应的权函数以求得其温度场的近似解

式中,Te(t)为单元e的近似温度,n为单元e所含的节点数量,和分别为单元e内节点i的权函数和温度,Ne和Te(t)分别为单元e内各节点的权函数向量和温度向量,nelem为求解域Ω 的离散单元数量.通过Galerkin 投影可以得到单元e的非线性代数方程组

式中,Ce,Ke和Fe分别是单元e的热容矩阵、热传导矩阵和温度载荷向量,其计算式如下

显然,Ce和Ke都是n×n方阵,而Fe是n×1 列阵.

于是,上述瞬态非线性热传导问题的有限元离散格式可由形如式(6)的求解域Ω 内各单元的非线性代数方程组组集而成,即

式中,T(t)为t时刻的节点温度向量

其中,nnode为求解域Ω 内的节点总量,而C,K和F分别是热传导系统的总体热容矩阵、总体热传导矩阵和总体温度载荷向量

其中,C和K都是nnode×nnode方阵,而F是nnode×1列阵.

由式(8)可以解得求解域Ω 的温度场.由于瞬态非线性热传导控制方程(1)的有限元模型的方程数量等于节点总量(即自由度数量),当求解域Ω 的规模较大、形状复杂时其内部须划分数量庞大的网格单元,导致求解需要很大的存储资源并耗费大量的计算时间.

1.3 瞬态非线性热传导有限元模型的计算方法

式(8)中总体热传导矩阵K是与温度T(t)相关的,因此该式是一个非线性微分方程组.本文采用无条件稳定的Euler 隐式格式

将其转化为非线性代数方程组,并应用N-R 方法进行迭代求解.

已知时刻tn的温度场Tn,第i次迭代后的残值可表示成

将残值R应用泰勒展开,令有

式中,KT为切向刚度矩阵,KT=∂R/(∂T),于是

当修正温度∆T或残值R达到收敛标准时,时刻tn+1的温度场为

2 瞬态非线性热传导方程的模型降阶法

2.1 瞬态非线性热传导问题的降阶模型

对于需要实时控制或快速计算的工程领域,有限元等常规的数值方法不能满足要求.因为实际问题的自由度数量高达上千万甚至更多,求解这种超大规模的非线性方程组需要很长的时间.为此需要建立降阶模型以减小计算量,从而达到实时控制或快速计算的目标.根据本课题组前期的研究[31],基于特征正交分解方法建立的降阶模型能够对瞬态非线性热传导过程进行准确预测,但是由于非线性问题的特殊性,降阶模型在计算效率的提升方面不够显著,本文致力于解决这一问题.为保持文章的完整性,下面简要介绍针对瞬态非线性热传导问题建立的降阶模型,而关于特征正交分解方法的详细内容可参阅文献[32].

首先,通过实验测量、数值模拟等方法获得瞬态非线性热传导问题在若干个时刻点上的温度场,通常称其为瞬像,将所有瞬像按照时间顺序排列成瞬像矩阵A=[T(t1)T(t2)···T(tl)]nnode×l.

其次,对瞬像矩阵A进行特征正交分解,获得一组正交基矢量,通常称其为POD 模态,截取前r阶能够捕捉绝大部分能量的模态构成模态矩阵Φ=

最后,利用模态矩阵Φ可以将任意时刻的温度场T(t)表示成

2.2 瞬态非线性热传导问题降阶模型的常规计算方法

采用隐式时间推进方法建立的瞬态非线性热传导降阶模型的离散格式可表示成

然而需要强调的是,由式(16)可见此降阶模型是不完备的,即仍然取决于原物理空间中的高维温度向量Tn+1,而不能直接通过当前时刻得到的低维向量来计算,使得在式(18)的计算过程当中每个迭代步都要出现下列额外的计算环节:

(2)根据新的温度场更新导热系数,利用式(7b)和式(10b)重新计算各单元的热传导矩阵并组集成非线性热传导系统的总体热传导矩阵Kn+1.

(3)通过式(17b)的矩阵运算将Kn+1转换成非线性热传导降阶模型的系数矩阵

由于有限元全阶模型的自由度数量较大,这些额外的数据转换、矩阵组集和矩阵运算将耗费不少的计算时间,使得降低模型的小规模非线性代数方程组的求解带来的时间收益被抵消.在后面的算例验证部分我们可以看到,在有限元全阶模型的规模不太大时,采用常规方法求解降阶模型对计算效率的提升效果并不理想.

2.3 瞬态非线性热传导问题降阶模型的新型计算方法

本小节将介绍瞬态非线性热传导降阶模型的一种新型高效计算方法,它是直接通过对上述常规计算方法中两个耗时环节的技术处理以达到提升计算效率的目的:一方面,采用单元预转换方法(element preconversation method,EPM)压缩降阶模型系数矩阵的计算时间;另一方面,采用多级线性化方法(multilevel linearization method,MLM)消除非线性方程组的迭代计算过程.这种新方法解决了导热系数随温度变化情况下EPM 和MLM 的结合问题,以及如何对换热边界单元的热传导矩阵进行相应计算.下面对此作详细阐述.

2.3.1 EPM 理论与计算

这要求在有限元全阶模型的单元层面就预先将单元热传导矩阵Ke进行转换,其计算式为

式中,Φe是模态矩阵Φ中单元e的对应部分,Φe是n×r矩阵.显然,与一样,也是r×r方阵.

将式(7b)改写成

于是,式(22)可以简化成

将上式依次代入式(21)、式(20),可得

当瞬态非线性热传导问题的有限元全阶模型和对流换热边界确定后,都可以预先计算并存储备用.在迭代过程中,采用式(27)计算降阶模型的系数矩阵时,只需要根据新的温度场计算出各个单元的平均导热系数即可.因此,与2.2 节的常规计算方法相比,EPM 能够有效压缩的计算时间,进而提高降阶模型的计算效率.

2.3.2 MLM 理论与计算

MLM 的主要思想是通过插值技术取代非线性热传导系统的非线性项,得到其近似的线性方程组,从而利用前面若干级时刻点上已知的温度场直接计算出当前时刻的温度场.方便起见,将式(8)表示的瞬态非线性热传导问题的有限元全阶模型改写成

式中,系数阵W=CK(T),载荷向量f=CF(t).

本文采用无条件稳定的Dupont II 多级线性化格式来求解式(28),此格式为[33]

式中,W∗表示插值温度T∗下的系数矩阵W,f∗表示插值时间t∗时的载荷向量f,T∗和t∗按下式计算

将式(15)代入式(29),然后方程两边同乘ΦT就可以得到瞬态非线性热传导降阶模型的多级线性化格式

式中,Ir为r阶单位矩阵,而分别为

由此降阶模型的多级线性化格式可见,在前两个时刻tn和tn+1的已知时,可先利用式(15)计算原物理空间的温度场Tn和Tn+1,然后根据式(30)计算插值温度T∗和插值时间t∗,继而求出W∗和f∗,并利用式(32)将它们转换为和就可以通过式(31)直接计算时刻tn+2的得到其温度场Tn+2.

实际计算时,初始时刻t0的温度场是已知的,而第二个时刻t1的温度场未知,需要通过2.2 节的常规方法或者基于前面EPM 方法改进的常规方法求出,此后各个时刻的温度场均可由多级线性化格式直接计算得到,消除了耗时的迭代过程,因而有利于提高降阶模型的计算效率.

2.3.3 EPM 和MLM 的结合技术

鉴于EPM 和MLM 分别是针对常规方法中两个不同环节做出的改进,这就为两者结合以获得降阶模型更高的计算效率提供了可能.然而,其结合的困难在于对MLM 方法中式(32)的处理,将W∗和f∗代入其中可见

该式要求进行总体热容矩阵的求逆和总体矩阵及其与模态矩阵之间的多重乘法运算,这与EPM 方法从单元层面的预处理方式是不协调的.

本文针对材料的比热容为常数、导热系数随温度变化这一类型的瞬态非线性热传导问题,提出有限元全阶模型中的各个单元对降价模型系数矩阵亦有其相应的贡献,贡献大小用表示,即

由式(7a)得到的单元热容矩阵Ce为协调质量热容矩阵,为了计算式(35)中的单元矩阵须将其转化为集中质量热容矩阵,这样由式(10a)得到的总体热容矩阵C为对角矩阵,可写为

由于比热容为常数,于是C是不随温度变化的恒常矩阵,其逆矩阵亦为对角矩阵,此时可由下式进行计算

将式(24)代入式(37),并补充单元处于对流换热边界上时存在的单元热传导矩阵修正项可得

式中,重复指标p表示执行爱因斯坦求和约定,对二维问题和三维问题p分别取2 和3.

将上式依次代入式(35)、式(34),可得

这样,就将EPM 和MLM 结合起来了.在用式(31)降阶模型多级线性化格式进行计算时,只需要根据前两个时刻已经得到的温度场计算出插值温度和插值时间即可,其余步骤和过程是与EPM 完全相同的.由此可见,这种针对瞬态非线性热传导降阶模型发展的新型EPM-MLM 复合方法,与EPM 具有良好的一致性,其区别仅在于将和替换成了和简单易行,便于编制规范统一的程序和软件.

3 算例分析

本部分通过二维和三维瞬态非线性热传导算例的计算和分析,验证本文所提出的降阶模型新型计算方法的可行性、准确性和高效性.为方便对比说明,下面将瞬态非线性热传导降阶模型的常规计算方法和新型计算方法分别命名为POD1 和POD2.

3.1 算例1:方板的瞬态非线性热传导问题

考虑如图1 所示边长为0.5 m 的二维方板,板材的密度ρ=7850 kg/m3、比热容c=460 J/(kg·◦C),其导热系数随温度线性变化,k=65-0.05TW/(m·◦C).方板的初始温度为20◦C,其右侧为对流换热边界,换热系数125 W/(m2·◦C)、环境温度1120◦C;左侧为放热边界,热流密度10 000 W/m2;上、下侧均为绝热边界.

图1 方板的几何模型及热边界条件Fig.1 Geometric model and thermal boundary conditions of the square plate

方板的有限元模型如图2 所示.计算域采用线性三角形单元进行离散,各侧边均布21 个节点,总计888 个单元、485 个节点.为便于后处理分析,取方板的左侧中点A、中心点B 和右侧中点C 等三个点作为参考点.

为建立方板瞬态非线性热传导问题的降阶模型,本文采取数值模拟的方法构造瞬像矩阵:时间步长取为50 s,通过有限元法计算0∼20 000 s 时间范围内方板在上述定常热边界条件下温度场的变化,选取t=1000,2000,···,20 000 s 共计20 个时刻点上方板的温度场构成瞬像矩阵A.对A进行特征正交分解分析,截取前4 个特征值(即r=4)对应的POD 模态构造模态矩阵Φ.利用式(15)可以得到仅包含4 个非线性方程的降阶模型,相比于原来由485 个非线性方程组成的有限元全阶模型,方程组的规模显著减小.

图2 方板的有限元模型Fig.2 Finite element model of the square plate

图3 给出了有限元全阶分析、常规降阶分析和新型降阶分析三种方法得到的方板在10 000 s 时刻的温度场对比,可以看到三者是高度一致的.图4 给出了方板上3 个参考点的温度时变曲线,它将常规降阶分析和新型降阶分析两种方法的计算结果与商业有限元软件ANSYS 对此全阶模型的计算结果(图中以实线表示)进行比较,可见新型降阶算法POD2 与常规的降阶算法POD1 一样,在所研究的整个时间范围内能够得到与有限元方法吻合的结果.因此,这些研究结果表明基于POD 的模型降阶方法能够用于瞬态非线性热传导问题的分析并得出与有限元方法相同的结果,而所提出的新型算法POD2 也可以得到正确的结果.

图3 三种算法10 000 s 时刻方板的温度场对比Fig.3 Comparison of temperature distributions in the square plate at t=10 000 s using different algorithms

图3 三种算法10 000 s 时刻方板的温度场对比(续)Fig.3 Comparison of temperature distributions in the square plate at t=10 000 s using different algorithms(continued)

图4 参考点处温度随时间的变化曲线Fig.4 Temperature evolutions at reference points

下面进一步对降阶模型新型算法的计算精度进行定量分析.首先,在方板的上边线上取等间距的6个点,将其10 000 s 时刻的计算温度和相对误差列于表1.相对误差是以有限元全阶模型的计算结果为基准的误差,其计算式为

可见,该时刻下POD2 算法的相对误差相比于POD1算法略微增加,这是因为POD2 算法不仅包含POD1算法的模态截断误差,还附加了线性化误差.但是POD2 算法的相对误差不超过0.2%,仍然保持了比较高的精度,因而这种新型算法的计算结果是准确的.

表1 方板上边线的计算温度及相对误差Table 1 Temperatures and relative errors on the top side of the square plate

其次,考虑所研究时间范围内各个时刻的方板整体温度的降阶模型计算误差,采用均方根误差(root mean square error,RMSE)这一指标来衡量其大小

图5 给出了RMSE 的时变曲线,可见降阶模型的计算误差在经过起始很短时间的脉动后迅速趋近于零,其峰值仅为1.1%.图中也显示出新型算法POD2 的误差略大于常规算法POD1,其原因如前所述.

图6 对有限元全阶分析、常规降阶分析和新型降阶分析3 种方法花费的计算时间进行了比较,可以直观看到POD1 用的计算时间较FEM 略有减少,而POD2 用的计算时间则大幅减少,这表明采用常规算法分析瞬态非线性热传导降阶模型时的计算效率并不比有限元全阶分析有明显改进,而采用新型算法对降阶分析计算效率的提高是非常显著的.

图5 方板温度的均方根误差随时间的变化曲线Fig.5 Temporal variation of the RMSE for Example 1

图6 各种算法所用计算时间的比较Fig.6 Comparison of computational times cost by various algorithms for Example 1

下面改变方板计算域的节点密度,使有限元全阶模型的自由度数量从129 逐渐增加到2941,分析这个因素对降阶模型POD1 和POD2 两种算法计算效率的影响.表2 中列出了不同自由度下各种算法所花费的计算时间,同时给出了POD1 和POD2 的倍率数值.倍率是有限元全阶分析与POD 降阶分析所用计算时间的比值,其大小可以衡量降阶分析对计算效率提升效果的强弱.由表可见,常规算法POD1 的计算效率提高得不多,尤其是在自由度数量较小时,降阶分析所用的计算时间几乎与全阶分析的相近.与之相比,新型算法POD2 的计算效率则提高了2 个数量级,而且在自由度数量较小时,降阶分析所用的计算时间也是大幅减少.此外,随着自由度数量的增大,新型算法POD2 对计算效率的增幅愈加显著.

表2 不同自由度时各种算法所用的计算时间Table 2 Computational times cost by various algorithms under different DOFs

最后,考虑到非线性热传导问题的有限元全阶分析往往要花费较长的计算时间,而且实际工程上面对着各种复杂的热边界条件,测试也比较的困难,导致获取相应问题的模态矩阵建立其降阶模型既耗时又不易.如果能够用简单的定常热边界条件下建立的降阶模型对相同计算域在复杂时变热边界条件下的瞬态非线性热传导过程进行分析和预测,将是非常有意义和应用价值的.为此,下面采用图1 中定常热边界条件下建立的方板瞬态非线性热传导降阶模型来分析其在新的时变热边界条件下的热传导过程,以检验这种方法的可行性.

保持方板右侧对流换热条件和上、下侧绝热条件不变,将其左侧改为时变的放热热流边界,如图7所示,有二种情况:

(a)时变边界1:热流密度按正弦波

的形式光滑连续变化,在0∼20 000 s 的时间范围内周期性变化4 次.

(b)时变边界2:热流密度按三角波的形式非光滑连续变化,波动幅度增大为时变边界1 的两倍,而在0∼20 000 s 的时间范围内周期性变化次数则减少到两次.

采用新型算法POD2 对前面定常热边界条件下建立的方板瞬态非线性热传导降阶模型进行计算,得到两种新的时变边界条件下方板的温度场,图8 所示为3 个参考点上的温度变化.由图可见,点A 处的温度明显呈现出波动状升温的过程,因为该点就处于左侧时变热边界上,甚至在时变边界2 的条件下,处于方板右侧边界上的点C 也由于左侧大幅度变化的热流边界而使其温度呈现出轻微的波动.为验证该降阶模型计算结果的准确性,图中同时以曲线形式给出了ANSYS 对相应边界条件下方板温度场的模拟结果,可见两者是完全吻合的,因此用定常热边界条件下建立的瞬态非线性热传导降阶模型能够准确有效的对相同计算域在各种复杂的时变热边界条件下的热传导过程进行分析与预测.

图7 方板左侧施加的时变热流Fig.7 Time-varying heat flux on the left side of the square plate

图8 时变边界下参考点处温度的时变曲线Fig.8 Temperature evolutions at reference points under time-varying boundary conditions

3.2 算例2:散热器的瞬态非线性热传导问题

考虑一个三维结构的散热器,如图9 所示,其底面为边长0.2 m 的方形,高为0.15 m.3 个翅片的厚度及其间距均为δ,δ=0.04 m,翅片高为0.1 m.材料的密度ρ=7800 kg/m3、比热容c=400 J/(kg·◦C)、导热系数随温度线性变化,k=60+0.05TW/(m·◦C).散热器的初始温度为20◦C,其底面是与发热体接触的加热面,热流密度50 000 W/m2;翅片与环境气体间进行对流换热,换热系数100 W/(m2·◦C)、环境温度20◦C;其余表面均为绝热面.

图9 散热器的几何模型及热边界条件Fig.9 Geometric model and thermal boundary conditions of the radiator

散热器的有限元模型如图10 所示.计算域采用六面体单元进行离散,总计4400 个单元、5796 个节点.为便于后处理,取中部翅片的角点A和B作为参考点,以及边线CD进行温度分析.

图10 散热器的有限元模型Fig.10 Finite element model of the radiator

与算例1 一样,为建立散热器的瞬态非线性热传导降阶模型,采取数值模拟的方法构造瞬像矩阵:时间步长取5 s,通过有限元法计算0∼500 s 时间范围内散热器温度场的变化,选取t=25,50,···,500 s 共计20 个时刻点上的温度场构成瞬像矩阵A.对A进行特征正交分解分析,截取前5 个特征值(即r=5)对应的POD 模态构造模态矩阵Φ.利用式(15)可以得到仅包含5 个非线性方程的降阶模型,相比于由5796 个非线性方程组成的有限元全阶模型,方程组的规模急剧减小.

图11 给出了有限元全阶分析、常规降阶分析和新型降阶分析3 种方法得到的散热器在末时刻500 s 时的温度场对比,可以看到三者基本上是相同的.图12 给出了散热器上A,B二个参考点的温度时变曲线,图中将POD1 和POD2 两种降阶模型算法的计算结果与商业有限元软件ANSYS 对全阶模型的计算结果进行了比较,同样可以看到新型算法POD2与常规算法POD1 一样能够得到与有限元方法吻合的结果.因此本文提出的新型算法POD2 可以对瞬态非线性热传导问题进行准确的分析.

图11 三种算法500 s 时刻散热器的温度场对比Fig.11 Comparison of temperature distributions in the radiator at t=500 s using different algorithms

图12 参考点处温度随时间的变化曲线Fig.12 Temperature evolutions at reference points

在散热器的边线CD 上取等间距的6 个点,将其500 s 时刻的计算温度和相对误差列于表3 中.由表可见,与算例1 类似,该时刻下POD2 算法的相对误差相比于POD1 算法略微增大.常规算法POD1 的最大相对误差约为1.12%,而新型算法POD2 的最大相对误差约为1.81%,因而这种新型算法仍然保持了比较高的精度,其计算结果是准确的.

表3 散热器边线CD 上的计算温度及相对误差Table 3 Temperatures and relative errors on the line CD of the radiator

图13 给出了用POD1 和POD2 两种降阶模型算法计算得到的散热器整体温度的均方根误差RMSE的时变曲线,可以看到降阶模型的计算误差经过起始时间的脉动,在约200 s 后缓慢稳定下降.POD1和POD2 的脉动峰值分别为0.01%和0.03%.降阶模型的计算误差在稳定段基本小于0.01%,新型算法POD2 的误差略大于常规算法POD1.

图14 对散热器算例采用有限元全阶分析、常规降阶分析和新型降阶分析3 种方法花费的计算时间进行了比较.由图可见,由于该算例的自由度数量比较大,有限元全阶分析需花费长达约862 s 的计算时间,此时常规降阶分析和新型降阶分析花费的计算时间均有不同程度减少,尤其是后者,用时仅为7 s,计算效率提高了将近123 倍.

图13 散热器温度的均方根误差随时间的变化曲线Fig.13 Temporal variation of the RMSE for Example 2

图14 各种算法所用计算时间的比较Fig.14 Comparison of computational times cost by various algorithms for Example 2

3.3 新型算法的特性分析

通过前面两个算例的计算分析,验证了本文提出的降阶模型新型算法是准确有效的,特别是在计算效率方面,较常规算法有大幅度的提升.

与近年来发展的针对瞬态非线性热传导问题的其它一些降阶算法相比,本文的新型算法也具有更高的计算效率.表4 中列有IGA-IC[30]、TLD-MS[28]和CA[29]等一些方法对自由度数量与本文算例尽量接近的算例所用的计算时间,二维问题和三维问题分栏给出.须说明的是,CA 方法的计算时间文献中未直接给出,是从其图上提取的,故标为约数.考虑到各个研究者所用计算机性能的不同,对计算时间是有影响的,故表中同时给出了倍率值(全阶与降阶计算时间的比值),以更客观的评价各种方法对降阶模型计算效率的提升效果.由表可见,本文的POD2 方法在计算效率上较其它方法是占优的,特别是三维问题.而二维问题中TLD-MS 方法的计算效率虽然接近POD2 方法,但该方法需要建立疏、密两套网格,算法繁杂,既存在网格匹配的要求,又会因不同网格间的插值造成精度损失,不利于通用软件的编制和实际工程的应用.

表4 新型算法和文献中其他方法的计算效率对比Table 4 Comparison of computational efficiency among POD2 and other methods in literatures

4 结论

本文针对导热系数随温度变化的一类瞬态非线性热传导问题,发展并提出了一种基于POD 模型降阶技术的新型高效快速分析方法.通过单元预转换方法和多级线性化方法的有效结合解决了非线性降阶模型的常规计算方法在系数矩阵更新和迭代等计算环节的费时问题,从而能够获得很高的计算效率.

通过二维和三维算例分析,演示了瞬态非线性热传导降阶模型新型计算方法的应用细节,通过与全阶模型有限元法计算结果的比较,验证了该新型算法的准确性,并通过与降阶模型常规计算方法以及近年来发展出的其他新方法在计算时间上的比较,验证了该新型算法的高效性.

研究结果表明,该新型算法具有良好的特性:其一,计算稳定性好,算法基于的隐式时间推进格式和多级线性化格式均为无条件稳定;其二,计算速度快,常规算法在全阶模型自由度数量较小时计算时间没有明显减少,而新型算法无论自由度数量是大是小都能够加快计算速度,并且自由度数量越多计算效率的提升效果越好,在本文自由度数量不超过6000 的小规模计算中其计算效率能提高2∼3 个数量级;其三,计算通用性好,便于编制规范、统一的程序和软件,有利于实际应用.

此外,采用常数热边界条件下得到的POD 模态建立的瞬态非线性热传导降阶模型可以对相同计算域在各种复杂的光滑/非光滑时变热边界条件下的非线性热传导过程进行快速、准确的分析和预测,这进一步拓展了本文所提降阶模型新型算法在工程应用中的价值.