季 奕 邢誉峰

* (北京航空航天大学固体力学研究所,北京 100083)

† (北京航空航天大学高等理工学院,北京 100083)

引言

时间积分方法被广泛应用在瞬态响应的求解中.根据递推公式的格式,其可被划分为两类,分别是显式方法和隐式方法.比较而言,显式方法效率高但条件稳定,隐式方法效率低但多为无条件稳定.对于非线性系统,确定显式方法的时间步长具有挑战性.因此,本文致力于发展对线性和非线性系统都是无条件稳定的隐式方法.

对瞬态热传导问题,广义梯形法则[1]是一种较为流行的方法,由它可衍生出Crank-Nicolson 方法、Galerkin 方法和后向差分方法(backward difference formula,BDF) 等.在这些方法中,只有Crank-Nicolson 方法具有二阶精度,但其在高频段无法提供数值耗散,从而会诱发出虚假的数值振荡.一阶精确的BDF 具有极强的数值耗散从而能迅速地过滤掉高频信息,但其无法准确地保留重要的低频信息.相比之下,Galerkin 方法能够有效地平衡低频精度和高频耗散.随着人们对精度、效率、耗散和稳定性的追求,在过去的几十年中,不断涌现出新的时间积分方法[2-31].

除减小步长外,提高算法阶次[2-6]也是一种提高精度的有效手段.针对一阶线性常微分方程,Fung借助加权残量法[2]、最小二乘法[3]和微分求积法[4]构造出一系列具有任意高阶精度的无条件稳定时间积分方法,但这些方法在求解过程中需要对系统矩阵进行扩维处理,从而大幅增加计算量.为了同时获得高精度和高效率,一种精细积分方法[7]被提出.该方法利用Taylor 级数来近似解析的传递矩阵,利用分段线性技术处理外载荷.此外,该方法采用了2m算法和增量矩阵存贮技术,大幅提高了效率、减少了舍入误差.虽然精细积分方法是条件稳定的,但由于在计算中采用了极小步长,因此条件稳定性并不影响其应用.精细时间积分方法已被广泛地应用在热传导[8-9]、双边值[10-11]和结构动力学[12-13]等问题中.

上面提到的优秀方法仅能求解线性问题,下面回顾一下能够求解一般非线性瞬态热传导问题的方法[14-31].Runge-Kutta 方法[14-16]是具有无条件稳定的高阶格式.但与单步方法相比,多级结构使其计算量更高.为了在不牺牲计算效率的前提下提高精度,线性多步方法[17-21]和多分步方法[22-28]被提出.一般来说,线性多步方法是二阶精确的,但其不能自启动.与单步方法相比,线性多步方法的使用略显繁琐.多分步方法是自启动的,但各个分步的有效刚度矩阵可能不同,致使其计算量略高于单步方法.瞬态热传导问题普遍存在于航空航天、土木和冶金等领域中,构造一种具有无条件稳定性且能够灵活处理该类问题的高精度、高效率时间积分方法是必要的.在这种背景下,本文提出了一种能处理一般瞬态热传导问题的无条件稳定单步时间积分方法,其具备二阶精度、高频耗散和自启动特性.利用拉格朗日插值函数分别近似温度场及其一阶导数场,并利用加权残量法建立二者之间的关系.理论分析和数值测试均显示出了本文方法在精度、耗散和稳定性方面的优势.

1 算法描述

瞬态热传导问题的非线性控制方程[32]可以写为

式中,C和K分别是N×N的热容矩阵和热传导矩阵,T和Q分别是N×1 的温度列向量和外部施加的热源列向量.当C和K是与温度T无关的定常矩阵时,式(1)退化为线性方程,即

为求解式(1)和式(2)中给出的瞬态热传导问题,本文首先把待求的时间域[0,tTotal]等间距地划分为n个时间单元,各离散点为0=t0

式中,ψj(t)是与第j个节点关联的拉格朗日插值函数,Tj和分别描述了与第j个节点关联的温度及其一阶导数.拉格朗日插值函数的定义为

其中,τ=(t−tk)/Δt,τj用来确定第j个节点的位置即tj=tk+τjΔt.在本文方法中,τ0=0,τ1=0.5,τ2=1.在式(3)中,利用拉格朗日函数作为基函数分别独立得到了T和,因此二者之间不存在关系T˙ =dT/dt.下面利用加权残量法来建立T和T˙ 之间的关系.定义如下加权残量r(t)

对于解析方法,上述残量在任意时刻均为0.为了最小化r(t),本文方法要求

式中,w(t)为权函数.将式(5)代入式(6)可得

选用不同的权函数w(t),可以获得不同数值特性的时间积分方法.但这种方式很难设计出具有无条件稳定和理想耗散的方法,为此本文没有直接使用权函数w(t),而是利用如下参数θk来控制算法的性能.

将式(4)和式(8)代入式(7)整理可得

本文方法的数值性能仅由参数θ1和θ2控制,在之后的章节中,它们将按所希望的性能要求进行设计.

2 改进Hughes 理论

经典Hughes[33]理论可用来评估已有时间积分方法在非线性瞬态热传导系统中的稳定性特性.为了能够设计出对该类非线性问题无条件稳定的新方法,本文对经典Hughes 理论进行了改进.

但大量的数值实验表明对线性系统无条件稳定的方法对非线性系统可能是失效的.在用时间积分方法求解非线性方程(1)时,由于不同时刻的C和K是不同的,因此利用线性化方法得到的不同时刻的热特征值和特征向量是不同的.为了能够考虑非线性系统的这个特征,Hughes[33]建立了如下模型

其中 λtk+1为tk+1=tk+Δt=(k+1)Δt时刻的热特征值.若系统是线性的,则 λtk+1不随时间变化,恒为常数,即 λtk+1=λ,此时模型(14)退化为经典线性模型(13).将时间积分方法的递推公式代入模型(14)可推出时变的传递因子A.若|A|≤ 1,则时间积分方法对单自由度非线性系统是无条件稳定的.对于非线性多自由度系统,对其控制方程(1)在不同时刻分别解耦时,同一阶次的φ在不同时刻是不同的.不同时刻的温度坐标和特征向量的顺序是根据对应时刻特征值λ的从小向大的顺序依次排列的.值得指出的是,由模型(14)推出的稳定性条件对单自由度非线性系统是充分且必要的,而对多自由度非线性系统是必要但不充分的.借助方程(14),Hughes[33]详细分析了广义梯形法则在非线性瞬态热传导系统中的稳定性.广义梯形法则[32]的递推公式为

将式(15)代入式(14)中可得传递因子为

将式(16)代入|A|≤ 1 中可得

从式(17)中可以发现,在广义梯形法则中只有BDF (α=1)对非线性系统是无条件稳定的,而对线性系统无条件稳定的Crank-Nicolson 方法(α=1/2)对非线性系统则是条件稳定的.另外,注意到对于显式方法(α=0),失稳原因是时间步长尺寸Δt不够小,而隐式方法(0 <α≤ 1)失稳的原因与热特征值λ演化规律相关.为了能够借助Hughes 理论设计出对非线性瞬态热传导系统是无条件稳定的隐式时间积分方法,本文在Hughes 理论中引入一个用来刻画λ演化规律的函数,其定义为

函数δt可以刻画λ在一个时间单元内[tk,tk+Δt]的变化情况.将函数δt代入模型(14)可得改进的模型为

利用改进的模型(19)得到的广义梯形法则的稳定条件为

比较式(20)和式(17)可以发现,由改进模型得到的稳定性条件与Hughes 理论得到的结果完全一致,从而说明改进理论在稳定性分析方面与Hughes理论是等价的.下面以广义梯形法则的稳定性条件(20)为例,简述如何利用改进Hughes 理论设计无条件稳定方法.从式(20)可以看出,当(− δtk+1α-α+1)≤0 成立时,不等式(20)恒成立.为了消除正的 δtk+1对稳定性的影响,α的取值区间应为α≥1.当α≥1 时,广义梯形法则是无条件稳定的,其中α=1 对应的是BDF.这种通过消除 δtk+1对稳定性影响的方法可以用来设计无条件稳定方法,并被用于本文方法的参数设计中.Hughes 理论原模型(14)虽然可以用于分析方法的稳定性,但难以用于设计无条件稳定方法.

3 算法设计与分析

将本文方法的递推公式(9)代入式(19)可得递推关系为

其中,时变的传递因子A和载荷算子L分别为

式中,Ω= λtkΔt,函数gi,和hi(i=1,2)为

其中 δtk+1/2= λtk+1/2/ λtk,δtk+1= λtk+1/ λtk,tk+1/2=tk+Δt/2=(k+1/2)Δt.

首先考虑本文方法的高频耗散特性.若一个无条件稳定时间积分方法的传递因子满足A|Ω→∞=0,则其被认为是L 型耗散方法.为了能够快速过滤掉高频振荡,L 型耗散是被期望的,为此要求

从而可以建立θ1和θ2的约束关系为

进而式(22)中的传递因子A和式(23)中的载荷算子L分别被更新为

下面讨论本文方法的稳定性.这里假设1/2≤θ1≤3/4,从而使得式(27)中的传递因子A的分母恒为正值,进而不等式|A|≤ 1 可等价表示为

式中,Ω∈[0,+∞),δtk+1/2∈(0,+∞),δtk+1∈(0,+∞).从式(29)中可以观察到,当0 < δtk+1<1 时,算法有可能失稳.为了消除这个不稳定因素,令θ1=3/4.另外可以注意到,对于线性系统(δtk+1/2≡ δtk+1≡ 1),不等式(29)是恒成立的.再次说明造成隐式方法对非线性问题失稳的原因与热特征值的时变性有关.至此,两个自由参数θ1和θ2确定完成,它们使得本文方法具有无条件稳定性和L 型耗散.

最后,利用局部截断误差σ[32]来分析本文方法的精度阶次.局部误差定义为

将式(27)和式(28)代入上式整理可得

可以看到对线性系统(δtk+1/2≡ δtk+1≡ 1),s0=s1=0,故本文方法具有严格的二阶精度.对于非线性系统,当 δtk+1/2和 δtk+1均趋近于1 时,s0和s1则趋于0,此时本文方法是近似二阶精确的.

为便于读者使用,下面给出本文方法对线性系统(2)的求解流程.求解过程中涉及tk+1/2和tk+1两个时刻的平衡方程,即

首先将递推公式(9)代入平衡方程(33)中可解出tk+1和tk+1/2两个时刻的温度,即

将求得的Ttk+1和Ttk+1/2回代到式(9) 即可得到,从而获得所有未知的变量信息.计算步骤如下:

第1 步准备工作:

(1) 构造热容矩阵C,热传导矩阵K,热源向量Q.

(2) 初始化T0和T˙0.

(3) 选择步长Δt和计算算法参数:

(4) 构造系数矩阵:

第2 步迭代运算:

(1) 计算tk+Δt时刻的有效载荷向量:

(2) 求解tk+1和tk+1/2时刻的温度:

(3) 求解tk+Δt时刻温度的一次导数:

对非线性系统(1),利用递推公式(9)、平衡方程(33)和Newton 迭代法,可求解出未知变量Ttk+1,Ttk+1/2和.

4 数值性能

本节对本文方法的无条件稳定性、L 型耗散和二阶精度进行数值验证和讨论.

4.1 传递因子

图1(a)给出了本文方法和一些现有方法在线性系统中的传递因子[32]曲线,其中解析传递因子的表达式为Aexact=exp(−λΔt).与Crank-Nicolson 方法和Galerkin 方法相比,本文方法具有高频耗散优势.与BDF 相比,本文方法有低频精度优势.图1(b)给出了本文方法在非线性系统下的传递因子曲线,可以看到对任意 δtk+1/2和 δtk+1的组合,在不同Ω值下|A|≤1 始终成立,并且几乎无数值振荡区.随着Ω值增加,A逐渐趋于0,说明本文方法对非线性问题具备较强的高频耗散能力.

图1 传递因子Fig.1 Amplification factor

4.2 收敛率

收敛率可检测时间积分方法的精度阶次,其定义为在指定时刻物理量的相对误差随步长减小而减小的速率.先测试线性系统[34],考虑如下方程

图2(a)给出了t=100 时的相对误差,可以看出本文方法与Crank-Nicolson 方法具有相同的斜率.从而说明对线性系统,本文方法具有严格的二阶精度.再测试非线性系统[34],考虑如下方程

该问题的解析解为T(t)=T0/(3σT03t+1)1/3.图2(b)给出了t=10 的相对误差,可以看出对非线性问题,本文方法也能够达到二阶精度,此时s0和s1趋于0.

图2 收敛率Fig.2 Convergence rate

5 数值测试

本节利用3 个算例来测试本文方法在不同类型瞬态热传导问题中的性能.本文方法的求解过程中涉及tk+1和tk+1/2两个时刻的平衡方程,而广义梯形法则仅求解tk+1时刻平衡方程.为保证精度比较在相同计算量下执行,在所有的算例中,本文方法的时间步长取为单步方法的两倍.在本节,除广义梯形法则外,具有二阶精确和L 型耗散的BDF2 方法[21]也被考虑,由于其不具备自启动特性,本文采用Crank-Nicolson 方法作为启动算法.此外,所有问题的参考解均由使用小步长的Crank-Nicolson 方法得到.

5.1 矩形功能梯度板中的二维热传导问题

考虑一个功能梯度板中的二维热传导问题[35],如图3 所知,假设导热系数kx1=kx2=1+(x1− 1)2+(x2− 1)2,密度ρ=1,比热容cv=1.整个域内施加单位初始温度,所有边界施加T=0 的温度场.该矩形域利用20×20 个4 结点矩形单元进行空间离散,广义梯形法则采用步长Δt=0.01.

图3 算例1 的几何尺寸和边界条件Fig.3 Geometry and boundary conditions for Ex.1

图4 给出了中点A(1,1)和角点B(1.9,1.9)的温度−时间曲线.可以看到在中点A处,所有方法均与参考解吻合得较好,但在边界点B处,Crank-Nicolson 方法存在明显的数值振荡.

图4 点A 和点B 的温度−时间曲线Fig.4 Temperatures of point A and point B versus time

为了进一步比较这些方法在域内的精度表现,图5 提供了所有方法在t=0.1 时刻温度误差的分布结果,可以看出本文方法不会像Crank-Nicolson 方法一样在边界处产生剧烈的数值振荡,且在所有方法中,本文方法的域内精度最高.BDF 和Galerkin 方法虽然没有数值振荡,但域内精度不令人满意.

图5 温度误差分布(t=0.1)Fig.5 Errors of temperature at t=0.1

5.2 一维非线性热传导问题

考虑一个材料参数与温度有关的一维杆[36],杆长l假设为1,材料参数为k=γT2+βT+η(γ,β,η∈R)和ρcv=1.在杆的两端施加T=0 的温度场,域内施加单位初始温度.利用20 个两结点单元对杆进行空间离散.

首先考虑γ=0,β=η=1 的情况,广义梯形法则的Δt=0.000 25.图6 给出了中点处(x=0.5)和边界处(x=0.95)的温度−时间曲线,可以看到本文提出的方法展示出明显的精度优势.

图6 温度−时间曲线(γ=0,β=η=1)Fig.6 Temperature-time history

图6 温度−时间曲线(γ=0,β=η=1) (续)Fig.6 Temperature-time history (continued)

为了展示本文方法的稳定性优势,另一种情况β=100,γ=η=1 被考虑.此时广义梯形法则的Δt=0.005.从图7 可以看到Crank-Nicolson 方法在计算一段时间后就失去了稳定性,而本文方法不仅没有失去稳定性,且在所有收敛的方法中具有最高的精度.

图7 温度−时间曲线(β=100,γ=η=1)Fig.7 Temperature-time history

5.3 具有外部热源的二维热传导问题

考虑一均质矩形板[9],如图8 所示,外部热源f(t)=cos(10t)从右边界持续不断地向域内输入,底部边界施加T=0 的温度场,其余边界为绝热状态.利用50×50 个4 结点矩形单元进行划分,广义梯形法则的Δt=0.005.

图8 算例3 的几何尺寸和边界条件Fig.8 Geometry and boundary conditions for Ex.3

图9 给出了点A(0.98,0.98)处的温度和温度一阶导数的结果,从绝对误差的数值上可以看出:本文方法要比二阶Crank-Nicolson 方法和BDF2 更加精确.总体来说,这3 种二阶精度算法的精度明显优于一阶精度的Galerkin 方法和BDF.

图9 点A 处温度和其一阶导数与时间的关系曲线Fig.9 Temperature and its derivative -time histories of point A

6 结论

本文提出了一种适用于求解一般瞬态热传导问题的无条件稳定单步时间积分方法.该方法具有无条件稳定性、二阶精度、L 型数值耗散和自启动特性.需要强调的是,与大多数现有时间积分方法不同,本文方法对线性和非线性瞬态热传导问题均是无条件稳定的.

线性和非线性数值比较结果表明,与著名的二阶Crank-Nicolson 方法相比,在计算量相同的前提下,本文方法具有更理想的精度和稳定性,并且不会在高频段诱发出虚假的数值振荡.与具有L 型耗散的二阶BDF2 相比,本文方法具有更高的精度且无需其他算法来启动计算.

鉴于本文方法具有优秀的数值性能和易执行性,故推荐其用于一般瞬态热传导问题的求解.