AN UNCONDITIONALLY STABLE METHOD FOR TRANSIENT HEAT CONDUCTION
-
摘要: 瞬态热传导问题普遍存在于航空航天、土木和冶金等领域中, 对这类问题精确、高效的数值求解方法一直备受关注. 为此, 本文提出了一种无条件稳定的单步时间积分方法. 在所提出的方法中, 拉格朗日插值函数被用来近似真实的温度场及其一次导数场, 之后, 加权残量法被用来建立二者之间的关系. 通过对算法参数的巧妙设计, 本文提出的方法具有二阶精度和L型数值耗散, 其中, L型耗散使得本文方法能够快速过滤掉虚假的高频振荡. 多数现有时间积分方法仅对线性系统具有无条件稳定性, 对非线性系统则是条件稳定的. 为此, 本文改进了Hughes针对一阶非线性热传导问题提出的时间积分方法稳定性评估理论, 并将改进的理论用于方法的参数设计中. 理论分析的结果表明本文方法对线性和非线性热传导系统均是无条件稳定的. 即使对于著名的Crank-Nicolson方法失稳的非线性热传导问题, 本文方法仍能给出稳定且精确的预测. 数值测试结果显示, 所提出的方法相较于当前流行的方法具有明显的精度、耗散和稳定性优势.Abstract: Time-dependent transient heat conduction problems are widely encountered in aerospace, civil engineering, metallurgical engineering, etc., and for such problems, accurate and fast numerical approaches have always attracted attention in the past decades. To achieve this goal, this paper proposes an unconditionally stable single-step time integration method for general transient heat conduction systems. In the proposed method, the temperature vector and its time derivative are formulated independently by the Langrage interpolation function, and then the relation between the temperature vector and its time derivative is defined with the weighted residual method. Theoretical analysis, including convergence rate and amplification factor, illustrates that the proposed method is strictly second-order accurate for the temperature vector and its time derivative, and it has the strong algorithmic dissipation (L-dissipation), meaning that it can quickly filter out the unwanted numerical oscillations in the high-frequency range. At present, most existing time integration methods, such as the Crank-Nicolson method and the Galerkin method, are unconditionally stable for linear transient heat conduction systems, but they are conditionally stable for nonlinear ones. To this end, this work improved the stability analysis theory for nonlinear transient heat conduction systems proposed by Hughes and used the improved stability analysis theory to design the free parameters of the proposed method. Because of this reason, the proposed method is unconditionally stable for both linear and nonlinear transient heat conduction problems. Due to the desirable algorithmic stability, the proposed method can still provide accurate and stable predictions for nonlinear transient heat conduction problems where the excellent Crank-Nicolson method fails. Some linear and nonlinear transient heat conduction problems are solved in this paper, and the results of these problems show that compared to the currently popular time integration methods, such as the Crank-Nicolson method and the backward difference formula, the proposed method enjoys noticeable advantages in accuracy, dissipation and stability.
-
引言
时间积分方法被广泛应用在瞬态响应的求解中. 根据递推公式的格式, 其可被划分为两类, 分别是显式方法和隐式方法. 比较而言, 显式方法效率高但条件稳定, 隐式方法效率低但多为无条件稳定. 对于非线性系统, 确定显式方法的时间步长具有挑战性. 因此, 本文致力于发展对线性和非线性系统都是无条件稳定的隐式方法.
对瞬态热传导问题, 广义梯形法则[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]可以写为
$$\left. \begin{array}{l} {\boldsymbol{C}}\left( {\boldsymbol{T}} \right){\dot{\boldsymbol T}} + {\boldsymbol{K}}\left( {\boldsymbol{T}} \right){\boldsymbol{T}} = {\boldsymbol{Q}} \\ {\boldsymbol{T}}\left( 0 \right) = {{\boldsymbol{T}}_0} \end{array} \right\}$$ (1) 式中, C和K分别是N × N的热容矩阵和热传导矩阵, T和Q分别是N × 1的温度列向量和外部施加的热源列向量. 当C和K是与温度T无关的定常矩阵时, 式(1)退化为线性方程, 即
$$\left. \begin{array}{l} {\boldsymbol{C\dot T}} + {\boldsymbol{KT}} = {\boldsymbol{Q}} \\ {\boldsymbol{T}}\left( 0 \right) = {{\boldsymbol{T}}_0} \end{array} \right\}$$ (2) 为求解式(1)和式(2)中给出的瞬态热传导问题, 本文首先把待求的时间域[0, tTotal]等间距地划分为n个时间单元, 各离散点为0 = t0 < t1 < ··· < tk < ··· < tn-1 < tn = tTotal, 其中tk = kΔt (k = 1, 2, ···, n), Δt = tTotal /n为时间单元长度. 在任一时间单元t∈[tk, tk + Δt]中, 温度场及其一阶导数场被近似为
$$\left. \begin{array}{l} {\dot{\boldsymbol T}}\left( t \right) = {\psi _0}\left( t \right){{{\dot{\boldsymbol T}}}_0} + {\psi _1}\left( t \right){{{\dot{\boldsymbol T}}}_1} + {\psi _2}\left( t \right){{{\dot{\boldsymbol T}}}_2} \\ {\boldsymbol{T}}\left( t \right) = {\psi _0}\left( t \right){{\boldsymbol{T}}_0} + {\psi _1}\left( t \right){{\boldsymbol{T}}_1} + {\psi _2}\left( t \right){{\boldsymbol{T}}_2} \\ \end{array} \right\}$$ (3) 式中, ψj(t)是与第j个节点关联的拉格朗日插值函数, Tj和
${{\dot{\boldsymbol T}}_j}$ 分别描述了与第j个节点关联的温度及其一阶导数. 拉格朗日插值函数的定义为$${\psi _j}\left( t \right) = \prod\limits_{k = 0 \atop \left( {k \ne j} \right) }^2 {\frac{{\tau - {\tau _k}}}{{{\tau _j} - {\tau _k}}}} $$ (4) 其中, τ = (t − tk)/Δt,
${\tau _j} $ 用来确定第j个节点的位置即tj = tk + τjΔt. 在本文方法中, τ0 = 0, τ1 = 0.5, τ2 = 1. 在式(3)中, 利用拉格朗日函数作为基函数分别独立得到了T和${\dot{\boldsymbol T}}$ , 因此二者之间不存在关系${\dot{\boldsymbol T}}$ = dT/dt. 下面利用加权残量法来建立T和${\dot{\boldsymbol T}}$ 之间的关系. 定义如下加权残量r(t)$$\begin{split} {\boldsymbol{r}}\left( t \right) =& {{\psi _0}\left( t \right){{{\dot{\boldsymbol T}}}_0} + {\psi _1}\left( t \right){{{\dot{\boldsymbol T}}}_1} + {\psi _2}\left( t \right){{{\dot{\boldsymbol T}}}_2}} -\\ &\left( {{{\dot \psi }_0}\left( t \right){{\boldsymbol{T}}_0} + {{\dot \psi }_1}\left( t \right){{\boldsymbol{T}}_1} + {{\dot \psi }_2}\left( t \right){{\boldsymbol{T}}_2}} \right) \end{split} $$ (5) 对于解析方法, 上述残量在任意时刻均为0. 为了最小化r(t), 本文方法要求
$$\int_0^{\Delta t} {w\left( t \right)\frac{{{{\rm{d}}^{i - 1}}{\boldsymbol{r}}\left( t \right)}}{{{\rm{d}}{t^{i - 1}}}}} {\rm{d}}t = {\boldsymbol{0}},\;\; {i = 1,2} $$ (6) 式中, w(t)为权函数. 将式(5)代入式(6)可得
$$\int_0^{\Delta t} {w\left( t \right)\left( {\sum\limits_{j = 0}^2 {\frac{{{{\rm{d}}^{i - 1}}{\psi _j}\left( t \right)}}{{{\rm{d}}{t^{i - 1}}}}} {{{\dot{\boldsymbol T}}}_j} - \sum\limits_{j = 0}^2 {\frac{{{{\rm{d}}^i}{\psi _j}\left( t \right)}}{{{\rm{d}}{t^i}}}{{\boldsymbol{T}}_j}} } \right)} {\rm{d}}t = {\boldsymbol{0}}$$ (7) 选用不同的权函数w(t), 可以获得不同数值特性的时间积分方法. 但这种方式很难设计出具有无条件稳定和理想耗散的方法, 为此本文没有直接使用权函数w(t), 而是利用如下参数θk来控制算法的性能.
$${\theta _k} = \frac{{\displaystyle\int_0^{\Delta t} {w\left( t \right){t^k}} {\rm{d}}t}}{{\Delta {t^k}\displaystyle\int_0^{\Delta t} {w\left( t \right)} {\rm{d}}t}},\;\; {k = 1,2} $$ (8) 将式(4)和式(8)代入式(7)整理可得
$$\left. \begin{array}{l} {{{\dot{\boldsymbol T}}}_{{t_{k + {1 / 2}}}}} = {\alpha _{11}}{{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}} + {\alpha _{12}}{{\boldsymbol{T}}_{{t_{k + 1}}}} + {\beta _1}{{\boldsymbol{T}}_{{t_k}}} + {\eta _1}{{{\dot{\boldsymbol T}}}_{{t_k}}} \\ {{{\dot{\boldsymbol T}}}_{{t_{k + 1}}}} = {\alpha _{21}}{{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}} + {\alpha _{22}}{{\boldsymbol{T}}_{{t_{k + 1}}}} + {\beta _2}{{\boldsymbol{T}}_{{t_k}}} + {\eta _2}{{{\dot{\boldsymbol T}}}_{{t_k}}} \\ \end{array} \right\}$$ (9) 式中
$$\begin{split} & {\boldsymbol{\alpha }}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{\alpha _{11}}}&{{\alpha _{12}}} \\ {{\alpha _{21}}}&{{\alpha _{22}}} \end{array}} \right]{\rm{ = }}\frac{1}{{4\Delta t\left( {2\theta _1^2 - \theta _2^{}} \right)}}\cdot\\ & \left[ {\begin{array}{*{20}{c}} { - 4\left( {8\theta _1^2 - 4\theta _1^{} - 4\theta _2^{} + 1} \right)}&{ {16\theta _1^2 - 4\theta _1^{} - 8\theta _2^{} + 1} } \\ { - 16\left( {4\theta _1^2 - 2\theta _1^{} - 2\theta _2^{} + 1} \right)}&{4\left( {8\theta _1^2 - 2\theta _1^{} - 4\theta _2^{} + 1} \right)} \end{array}} \right] \end{split}$$ (10) $${\boldsymbol{\beta }}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{\beta _1}} \\ {{\beta _2}} \end{array}} \right]{\rm{ = }}\frac{1}{{4\Delta t\left( {2\theta _1^2 - \theta _2^{}} \right)}}\left[ {\begin{array}{*{20}{c}} { {16\theta _1^2 - 12\theta _1^{} - 8\theta _2^{} + 3} } \\ {4\left( {8\theta _1^2 - 6\theta _1^{} - 4\theta _2^{} + 3} \right)} \end{array}} \right]$$ (11) $${\boldsymbol{\eta }}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{\eta _1}} \\ {{\eta _2}} \end{array}} \right]{\rm{ = }}\frac{1}{{4\left( {2\theta _1^2 - \theta _2^{}} \right)}}\left[ {\begin{array}{*{20}{c}} { {8\theta _1^2 - 4\theta _1^{} - 4\theta _2^{} + 1} } \\ {4\left( {2\theta _1^2 - 2\theta _1^{} - \theta _2^{} + 1} \right)} \end{array}} \right]$$ (12) 本文方法的数值性能仅由参数θ1和θ2控制, 在之后的章节中, 它们将按所希望的性能要求进行设计.
2. 改进Hughes理论
经典Hughes[33]理论可用来评估已有时间积分方法在非线性瞬态热传导系统中的稳定性特性. 为了能够设计出对该类非线性问题无条件稳定的新方法, 本文对经典Hughes理论进行了改进.
多自由度线性系统(2)可以借助
${\boldsymbol{T}}=\displaystyle\sum\limits_i^N {{T^i}{\varphi ^i}}$ 解耦为N个关于温度坐标Ti的方程${\dot T^i}$ + λiTi = ri, 其中λ为由(K−λC)φ = 0得到的热特征值, φ为热特征向量. 因此, 大多数已有时间积分方法的理论分析是基于如下线性模型[32]$$\dot T + \lambda T = r$$ (13) 但大量的数值实验表明对线性系统无条件稳定的方法对非线性系统可能是失效的. 在用时间积分方法求解非线性方程(1)时, 由于不同时刻的C和K是不同的, 因此利用线性化方法得到的不同时刻的热特征值和特征向量是不同的. 为了能够考虑非线性系统的这个特征, Hughes[33]建立了如下模型
$${\dot T_{{t_{k + 1}}}} + {\lambda _{{t_{k + 1}}}}{T_{{t_{k + 1}}}} = {r_{{t_{k + 1}}}}$$ (14) 其中
${\lambda _{{t_{k + 1}}}}$ 为tk+1 = tk + Δt = (k + 1)Δt时刻的热特征值. 若系统是线性的, 则${\lambda _{{t_{k + 1}}}}$ 不随时间变化, 恒为常数, 即${\lambda _{{t_{k + 1}}}}$ = λ, 此时模型(14)退化为经典线性模型(13). 将时间积分方法的递推公式代入模型(14)可推出时变的传递因子A. 若|A| ≤ 1, 则时间积分方法对单自由度非线性系统是无条件稳定的. 对于非线性多自由度系统, 对其控制方程(1)在不同时刻分别解耦时, 同一阶次的φ在不同时刻是不同的. 不同时刻的温度坐标和特征向量的顺序是根据对应时刻特征值λ的从小向大的顺序依次排列的. 值得指出的是, 由模型(14)推出的稳定性条件对单自由度非线性系统是充分且必要的, 而对多自由度非线性系统是必要但不充分的. 借助方程(14), Hughes[33]详细分析了广义梯形法则在非线性瞬态热传导系统中的稳定性. 广义梯形法则[32]的递推公式为$${T_{{t_{k + 1}}}}{\rm{ = }}{T_{{t_k}}} + \Delta t\left[ {\left( {1 - \alpha } \right){{\dot T}_{{t_k}}} + \alpha {{\dot T}_{{t_{k + 1}}}}} \right]$$ (15) 将式(15)代入式(14)中可得传递因子为
$$A = \frac{{1 - {\lambda _{{t_k}}}\Delta t\left( {1 - \alpha } \right)}}{{1 + {\lambda _{{t_{k + 1}}}}\Delta t\alpha }}$$ (16) 将式(16)代入|A| ≤ 1中可得
$$\left. \begin{array}{l} \alpha = 0:\Delta t \leqslant 2/{\lambda _{{t_k}}}\\ \alpha \in \left( {0,1} \right):\\ \;\;\;\;\;\;{\rm{stable}},\;\;{\lambda _{{t_{k + 1}}}} \geqslant \dfrac{{\left( {1 - \alpha } \right)}}{\alpha }{\lambda _{{t_k}}}\\ \;\;\;\;\;\;\Delta t \leqslant \dfrac{2}{{\left( {1 - \alpha } \right){\lambda _{{t_k}}} - \alpha {\lambda _{{t_{k + 1}}}}}},\;\;{\lambda _{{t_{k + 1}}}} < \dfrac{{\left( {1 - \alpha } \right)}}{\alpha }{\lambda _{{t_k}}}\\ \alpha = 1:{\rm{stable}} \end{array} \right\} $$ (17) 从式(17)中可以发现, 在广义梯形法则中只有BDF (α = 1)对非线性系统是无条件稳定的, 而对线性系统无条件稳定的Crank-Nicolson方法(α = 1/2)对非线性系统则是条件稳定的. 另外, 注意到对于显式方法(α = 0), 失稳原因是时间步长尺寸Δt不够小, 而隐式方法(0 < α
$\leqslant$ 1)失稳的原因与热特征值λ演化规律相关. 为了能够借助Hughes理论设计出对非线性瞬态热传导系统是无条件稳定的隐式时间积分方法, 本文在Hughes理论中引入一个用来刻画λ演化规律的函数, 其定义为$${\delta _t} = \frac{{{\lambda _t}}}{{{\lambda _{{t_k}}}}},t \in \left[ {{t_k},{t_k} + \Delta t} \right]$$ (18) 函数δt可以刻画λ在一个时间单元内[tk, tk + Δt]的变化情况. 将函数δt代入模型(14)可得改进的模型为
$${\dot T_{{t_{k + 1}}}} + {\delta _{{t_{k + 1}}}}{\lambda _{{t_k}}}{T_{{t_{k + 1}}}} = {r_{{t_{k + 1}}}}$$ (19) 利用改进的模型(19)得到的广义梯形法则的稳定条件为
$$\left( { - {\delta _{{t_{k + 1}}}}\alpha - \alpha + 1} \right){\lambda _{{t_k}}}\Delta t \leqslant 2,\;\;{\rm{ }}{\delta _{{t_{k + 1}}}} = \frac{{{\lambda _{{t_{k + 1}}}}}}{{{\lambda _{{t_k}}}}} \in \left( {0, + \infty } \right)$$ (20) 比较式(20)和式(17)可以发现, 由改进模型得到的稳定性条件与Hughes理论得到的结果完全一致, 从而说明改进理论在稳定性分析方面与Hughes理论是等价的. 下面以广义梯形法则的稳定性条件(20)为例, 简述如何利用改进Hughes理论设计无条件稳定方法. 从式(20)可以看出, 当(−
${\delta _{{t_{k + 1}}}}$ α-α + 1)$\leqslant $ 0成立时, 不等式(20)恒成立. 为了消除正的${\delta _{{t_{k + 1}}}}$ 对稳定性的影响, α的取值区间应为α$\geqslant$ 1. 当α$\geqslant $ 1时, 广义梯形法则是无条件稳定的, 其中α = 1对应的是BDF. 这种通过消除${\delta _{{t_{k + 1}}}}$ 对稳定性影响的方法可以用来设计无条件稳定方法, 并被用于本文方法的参数设计中. Hughes理论原模型(14)虽然可以用于分析方法的稳定性, 但难以用于设计无条件稳定方法.3. 算法设计与分析
将本文方法的递推公式(9)代入式(19)可得递推关系为
$${T_{{t_{k + 1}}}} = A{T_{{t_k}}} + L$$ (21) 其中, 时变的传递因子A和载荷算子L分别为
$$A = \frac{{{\varOmega }_{}^2{g_2} + {\varOmega }{g_1} + 2}}{{{\varOmega }_{}^2{h_2} + {\varOmega }{h_1} + 2}}\quad\;$$ (22) $$L = - \Delta t\frac{{{\varOmega }{{\bar g}_2} + {{\bar g}_1}}}{{{\varOmega }_{}^2{h_2} + {\varOmega }{h_1} + 2}}$$ (23) 式中, Ω =
${\lambda _{{t_k}}}$ Δt, 函数gi,${\bar g_i}$ 和hi (i = 1,2)为$$\begin{split} {g_1} = &1 - 4\theta _1^{} - 4\theta _2^{} + 8\theta _1^2 + \\ &\quad{\delta _{{t_{k + {1 / 2}}}}}\left( { - 3 + 6\theta _1^{} + 4\theta _2^{} - 8\theta _1^2} \right) \\ \end{split} \tag{24a}\;\;\;\;\;\;\;\;\;\;\;\, $$ $${g_2} = {\delta _{{t_{k + {1 / 2}}}}}\left( {1 + 2\theta _1^2 - 2\theta _1^{} - \theta _2^{}} \right) \tag{24b} \;\;\;\;\;\;\;\;\;\;\;\;$$ $$\begin{split} {{\bar g}_1} = &\left( {{r_{{t_k}}} + {r_{{t_{k + 1}}}}} \right)\left( {1 - 4\theta _1^{} - 4\theta _2^{} + 8\theta _1^2} \right) - \\ &\quad 4{r_{{t_{k + {1 / 2}}}}}\left( {1 - 2\theta _1^{} - 2\theta _2^{} + 4\theta _1^2} \right) \\ \end{split} \tag{24c}\;\;\;\;\;\;\,$$ $$\begin{split} {\bar g_2} =& {\delta _{{t_{k + {1 / 2}}}}}{r_{{t_k}}}\left( {1 - 2\theta _1^{} - \theta _2^{} + 2\theta _1^2} \right) +\\ &\quad{\delta _{{t_{k + {1 / 2}}}}}{r_{{t_{k + 1}}}}\left( {\theta _2^{} - 2\theta _1^2} \right) \end{split} \tag{24d}\;\;\;\;\;\;\;\;\,\,\,\,\,\,\,$$ $$\begin{split} {h_1} =& \left( {{\delta _{{t_{k + {1 / 2}}}}} - {\delta _{{t_{k + 1}}}}} \right)\left( {1 + 8\theta _1^2} \right) - \\ &\quad2{\delta _{{t_{k + {1 / 2}}}}}\left( {\theta _1^{} + 2\theta _2^{}} \right) + 4{\delta _{{t_{k + 1}}}}\left( {\theta _1^{} + \theta _2^{}} \right) \end{split} \tag{24e}\,\, $$ $${h_2} = {\delta _{{t_{k + {1 / 2}}}}}{\delta _{{t_{k + 1}}}}\left( {2\theta _1^2 - \theta _2^{}} \right) \tag{24f}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\,\,\,\,$$ 其中
${\delta _{{t_{k + {1 / 2}}}}}$ =${\lambda _{{t_{k + {1 / 2}}}}}$ /${\lambda _{{t_k}}}$ ,${\delta _{{t_{k + 1}}}}$ =${\lambda _{{t_{k + 1}}}}$ /${\lambda _{{t_k}}}$ , tk+1/2 = tk + Δt/2 = (k + 1/2)Δt.首先考虑本文方法的高频耗散特性. 若一个无条件稳定时间积分方法的传递因子满足A|Ω→∞ = 0, 则其被认为是L型耗散方法. 为了能够快速过滤掉高频振荡, L型耗散是被期望的, 为此要求
$${\left. A \right|_{{\varOmega } \to \infty }}{\rm{ = }}0 \Rightarrow {g_2}{\rm{ = 0}}$$ (25) 从而可以建立θ1和θ2的约束关系为
$$\theta _2^{} = 2\theta _1^2 - 2\theta _1^{} + 1$$ (26) 进而式(22)中的传递因子A和式(23)中的载荷算子L分别被更新为
$$\begin{split} A = &\left\{\varOmega \left[ {{\delta _{{t_{k + 1/2}}}}\left( {1 - 2\theta _1^{}} \right) - \left( {3 - 4\theta _1^{}} \right)} \right] + 2\right\}\Big/\left\{\varOmega _{}^2{\delta _{{t_{k + 1/2}}}}{\delta _{{t_{k + 1}}}}\right.\cdot\\ &\;\;\;\;\;\;\;\;\left.\left( {2\theta _1^{} \!-\! 1} \right) + \varOmega \left[ {{\delta _{{t_{k + 1}}}}\left( {3 - 4\theta _1^{}} \right) + 3{\delta _{{t_{k + 1/2}}}}\left( {2\theta _1^{} - 1} \right)} \right] \!+\! 2\right\} \end{split}$$ (27) $$\begin{split} L = &\Delta t\;\left\{ {\varOmega {\delta _{{t_{k + 1/2}}}}{r_{{t_{k + 1}}}}\left( {2{\theta _1} - 1} \right) + \left[ {3{r_{{t_k}}} - 4{r_{{t_{k + 1/2}}}} + 3{r_{{t_{k + 1}}}} - } \right.} \right.\\ &\quad \left. {\left. {4{\theta _1}\left( {{r_{{t_k}}} - 2{r_{{t_{k + 1/2}}}} + {r_{{t_{k + 1}}}}} \right)} \right]} \right\}/\left\{{\varOmega _{}^2{\delta _{{t_{k + 1/2}}}}{\delta _{{t_{k + 1}}}}\left( {2{\theta _1} - 1} \right) + } \right.\\ &\quad \left. {\varOmega \left[ {{\delta _{{t_{k + 1}}}}\left( {3 - 4\theta _1^{}} \right) + 3{\delta _{{t_{k + 1/2}}}}\left( {2\theta _1^{} - 1} \right)} \right] + 2} \right\}\\[-12pt] \end{split}$$ (28) 下面讨论本文方法的稳定性. 这里假设1/2
$\leqslant$ θ1$\leqslant $ 3/4, 从而使得式(27)中的传递因子A的分母恒为正值, 进而不等式|A|$\leqslant $ 1可等价表示为$$ \begin{split} & 0\leqslant {\varOmega }_{}^{2}{\delta }_{{t}_{k+1/2}}{\delta }_{{t}_{k+1}}\left(2{\theta }_{1}^{}-1\right)+2\varOmega {\delta }_{{t}_{k+1/2}}\left(2{\theta }_{1}^{}-1\right)+4+\\ &\qquad\underset{{\text{不稳定因素}}}{\underbrace{\varOmega \left({\delta }_{{t}_{k+1}}-1\right)\left(3-4{\theta }_{1}^{}\right)}}\end{split}$$ (29) 式中, Ω∈[0, +∞),
${\delta _{{t_{k + {1 / 2}}}}}$ ∈(0, +∞),${\delta _{{t_{k + 1}}}}$ ∈(0, +∞). 从式(29)中可以观察到, 当0 <${\delta _{{t_{k + 1}}}}$ < 1时, 算法有可能失稳. 为了消除这个不稳定因素, 令θ1 = 3/4. 另外可以注意到, 对于线性系统(${\delta _{{t_{k + {1 / 2}}}}}$ ≡${\delta _{{t_{k + 1}}}}$ ≡ 1), 不等式(29)是恒成立的. 再次说明造成隐式方法对非线性问题失稳的原因与热特征值的时变性有关. 至此, 两个自由参数θ1和θ2确定完成, 它们使得本文方法具有无条件稳定性和L型耗散.最后, 利用局部截断误差σ[32]来分析本文方法的精度阶次. 局部误差定义为
$$\sigma = {{\left( {{T_{{t_{k + 1}}}} - A{T_{{t_k}}} - L} \right)} / {\Delta t}}$$ (30) 将式(27)和式(28)代入上式整理可得
$$\sigma = {s_0}O\left( {\Delta {t^0}} \right) + {s_1}O\left( {\Delta t} \right) + {s_2}O\left( {\Delta {t^2}} \right)$$ (31) 式中
$$\begin{split} {s_0} = &4\left( {{{\dot T}_{{t_k}}} + {\delta _{{t_{k + {1 / 2}}}}}{\lambda _{{t_k}}}{T_{{t_k}}} - {r_{{t_k}}}} \right) + {\varOmega }\left[ \left( {3{\delta _{{t_{k + {1 / 2}}}}} - 2} \right){{\dot T}_{{t_k}}} +\right.\\ &\left. {\delta _{{t_{k + {1 / 2}}}}}{\delta _{{t_{k + 1}}}}{\lambda _{{t_k}}}{T_{{t_k}}} - {\delta _{{t_{k + {1 / 2}}}}}{r_{{t_k}}} \right] \end{split} \tag{32a}$$ $${s_1} = \frac{1}{2}{\varOmega }\left[ {\left( {3{\delta _{{t_{k + {1 / 2}}}}} - 1} \right){{\ddot T}_{{t_k}}} + 2{\delta _{{t_{k + {1 / 2}}}}}{\delta _{{t_{k + 1}}}}{\lambda _{{t_k}}}{{\dot T}_{{t_k}}} - 2{\delta _{{t_{k + {1 / 2}}}}}{{\dot r}_{{t_k}}}} \right] \tag{32b}$$ $${s_2} = \frac{1}{8}{\delta _{{t_{k + {1 / 2}}}}}{\varOmega }\left( {3{{\dddot T}_{{t_k}}} + 4{\delta _{{t_{k + 1}}}}{\lambda _{{t_k}}}{{\ddot T}_{{t_k}}} - 4{{\ddot r}_{{t_k}}}} \right) \tag{32c}$$ 可以看到对线性系统(
${\delta _{{t_{k + {1 / 2}}}}}$ ≡${\delta _{{t_{k + 1}}}}$ ≡ 1), s0 = s1 = 0, 故本文方法具有严格的二阶精度. 对于非线性系统, 当${\delta _{{t_{k + {1 / 2}}}}}$ 和${\delta _{{t_{k + 1}}}}$ 均趋近于1时, s0和s1则趋于0, 此时本文方法是近似二阶精确的.为便于读者使用, 下面给出本文方法对线性系统(2)的求解流程. 求解过程中涉及tk+1/2和tk+1两个时刻的平衡方程, 即
$$\left. \begin{array}{l} {\boldsymbol{C}}{{{\dot{\boldsymbol T}}}_{{t_{k + {1 / 2}}}}} + {\boldsymbol{K}}{{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}} = {\boldsymbol{Q}}\left( {{t_{k + {1 / 2}}}} \right) \\ {\boldsymbol{C}}{{{\dot{\boldsymbol T}}}_{{t_{k + 1}}}} + {\boldsymbol{K}}{{\boldsymbol{T}}_{{t_{k + 1}}}} = {\boldsymbol{Q}}\left( {{t_{k + 1}}} \right) \\ \end{array} \right\}$$ (33) 首先将递推公式(9)代入平衡方程(33)中可解出tk+1和tk+1/2两个时刻的温度, 即
$$\begin{split} {{\boldsymbol{T}}_{{t_{k + 1}}}}{\rm{ = }}&{\left( {{\alpha _{12}}{\boldsymbol{C}} - {{\boldsymbol{C}}_4}} \right)^{ - 1}}\left[ {\left( {{{\boldsymbol{Q}}_{{t_{k + 1/2}}}} - {{\boldsymbol{C}}_5}{{\boldsymbol{Q}}_{{t_{k + 1}}}}} \right) + } \right.\\ &\quad \left. {\left( {{{\boldsymbol{\beta}} _2}{{\boldsymbol{C}}_5} - {{\boldsymbol{\beta}} _1}} \right){\boldsymbol{C}}{{\boldsymbol{T}}_{{t_k}}} + \left( {{{\boldsymbol{\eta}} _2}{{\boldsymbol{C}}_5} - {\eta _1}} \right){\boldsymbol{C}}{{\mathop T\limits^. }_{{t_k}}}} \right] \end{split}\;\;\;$$ (34) $${{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}}{\rm{ = }}{{\boldsymbol{C}}_1}\left( {{{\boldsymbol{Q}}_{{t_{k + 1}}}} - {{\boldsymbol{C}}_3}{{\boldsymbol{T}}_{{t_{k + 1}}}} - {\beta _2}{\boldsymbol{C}}{{\boldsymbol{T}}_{{t_k}}} - {\eta _2}{\boldsymbol{C}}{{{\dot{\boldsymbol T}}}_{{t_k}}}} \right)$$ (35) 式中
$$\left. \begin{array}{l} {{\boldsymbol{C}}_1} = {\left( {{\alpha _{21}}{\boldsymbol{C}}} \right)^{ - 1}},{{\boldsymbol{C}}_2} = \left( {{\alpha _{11}}{\boldsymbol{C}} + {\boldsymbol{K}}} \right),{{\boldsymbol{C}}_3} = \left( {{\alpha _{22}}{\boldsymbol{C}} + {\boldsymbol{K}}} \right)\\ {{\boldsymbol{C}}_4} = {{\boldsymbol{C}}_2}{{\boldsymbol{C}}_1}{{\boldsymbol{C}}_3},{{\boldsymbol{C}}_5} = {{\boldsymbol{C}}_2}{{\boldsymbol{C}}_1} \end{array} \right\} $$ (36) 将求得的
${{\boldsymbol{T}}_{{t_{k + 1}}}}$ 和${{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}}$ 回代到式(9)即可得到${{\dot{\boldsymbol T}}_{{t_{k + 1}}}}$ 和${{\dot{\boldsymbol T}}_{{t_{k + {1 / 2}}}}}$ , 从而获得所有未知的变量信息. 计算步骤如下:第1步 准备工作:
(1) 构造热容矩阵C, 热传导矩阵K, 热源向量Q.
(2) 初始化T0和
${{\dot{\boldsymbol T}}_0}$ .(3) 选择步长Δt和计算算法参数:
θ2, α11, α12, α21, α22, β1, β2, η1, η2.
(4) 构造系数矩阵:
$$\begin{array}{l} {{\boldsymbol{C}}_1} = {\left( {{\alpha _{21}}{\boldsymbol{C}}} \right)^{ - 1}},{{\boldsymbol{C}}_2} = \left( {{\alpha _{11}}{\boldsymbol{C}} + {\boldsymbol{K}}} \right),{{\boldsymbol{C}}_3} = \left( {{\alpha _{22}}{\boldsymbol{C}} + {\boldsymbol{K}}} \right), \\ {{\boldsymbol{C}}_4} = {{\boldsymbol{C}}_2}{{\boldsymbol{C}}_1}{{\boldsymbol{C}}_3},{{\boldsymbol{C}}_5} = {{\boldsymbol{C}}_2}{{\boldsymbol{C}}_1} \\ \end{array} $$ (5) 构造有效刚度矩阵
${{\overset{\frown}{\boldsymbol{ K}}}}$ :${{\overset{\frown}{\boldsymbol{ K}}}} = \left( {{\alpha _{12}}{\boldsymbol{C}} - {{\boldsymbol{C}}_4}} \right)$ .第2步 迭代运算:
(1) 计算tk +Δt时刻的有效载荷向量:
$${{\overset{\frown}{\boldsymbol{ Q}}}} = \left( {{{\boldsymbol{Q}}_{{t_{k + {1 / 2}}}}} - {{\boldsymbol{C}}_5}{{\boldsymbol{Q}}_{{t_{k + 1}}}}} \right) + \left( {{\beta _2}{{\boldsymbol{C}}_5} - {\beta _1}} \right){\boldsymbol{C}}{{\boldsymbol{T}}_{{t_k}}} + \left( {{\eta _2}{{\boldsymbol{C}}_5} - {\eta _1}} \right){\boldsymbol{C}}{{\dot{\boldsymbol T}}_{{t_k}}}$$ (2) 求解tk+1和tk+1/2时刻的温度:
$$ {{\boldsymbol{T}}_{{t_{k + 1}}}} \!\!=\!\! {{{\overset{\frown}{\boldsymbol{ K}}}}^{ - 1}}{{\overset{\frown}{\boldsymbol{ Q}}}},\; {{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}}\!\! =\!\! {{\boldsymbol{C}}_1}\left( {{{\boldsymbol{Q}}_{{t_{k + 1}}}} - {{\boldsymbol{C}}_3}{{\boldsymbol{T}}_{{t_{k + 1}}}} - {\beta _2}{\boldsymbol{C}}{{\boldsymbol{T}}_{{t_k}}} - {\eta _2}{\boldsymbol{C}}{{{\dot{\boldsymbol T}}}_{{t_k}}}} \right)$$ (3) 求解tk +Δt时刻温度的一次导数:
$${{\dot{\boldsymbol T}}_{{t_{k + 1}}}}{\rm{ = }}{\alpha _{21}}{{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}} + {\alpha _{22}}{{\boldsymbol{T}}_{{t_{k + 1}}}} + {\beta _2}{{\boldsymbol{T}}_{{t_k}}} + {\eta _2}{{\dot{\boldsymbol T}}_{{t_k}}}$$ 对非线性系统(1), 利用递推公式(9)、平衡方程(33)和Newton迭代法, 可求解出未知变量
${{\boldsymbol{T}}_{{t_{k + 1}}}}$ ,${{\boldsymbol{T}}_{{t_{k + {1 / 2}}}}}$ 和${{\dot{\boldsymbol T}}_{{t_{k + 1}}}}$ .4. 数值性能
本节对本文方法的无条件稳定性、L型耗散和二阶精度进行数值验证和讨论.
4.1 传递因子
图1(a)给出了本文方法和一些现有方法在线性系统中的传递因子[32]曲线, 其中解析传递因子的表达式为Aexact = exp(−λΔt). 与Crank-Nicolson方法和Galerkin方法相比, 本文方法具有高频耗散优势. 与BDF相比, 本文方法有低频精度优势. 图1(b)给出了本文方法在非线性系统下的传递因子曲线, 可以看到对任意
${\delta _{{t_{k + {1 / 2}}}}}$ 和${\delta _{{t_{k + 1}}}}$ 的组合, 在不同Ω值下|A|$\leqslant $ 1始终成立, 并且几乎无数值振荡区. 随着Ω值增加, A逐渐趋于0, 说明本文方法对非线性问题具备较强的高频耗散能力.4.2 收敛率
收敛率可检测时间积分方法的精度阶次, 其定义为在指定时刻物理量的相对误差随步长减小而减小的速率. 先测试线性系统[34], 考虑如下方程
$$ \dot{T}+T=10\mathrm{cos}\left(0.1t\right),\;{T}_{0}=100$$ (37) 图2(a)给出了t = 100时的相对误差, 可以看出本文方法与Crank-Nicolson方法具有相同的斜率. 从而说明对线性系统, 本文方法具有严格的二阶精度. 再测试非线性系统[34], 考虑如下方程
$$ \dot{T}+\sigma {T}^{4}=0,\;{T}_{0}=200,\sigma =3.629{\rm{e}}^{-8}$$ (38) 该问题的解析解为T(t) = T0/(3σT03t + 1)1/3. 图2(b)给出了t = 10的相对误差, 可以看出对非线性问题, 本文方法也能够达到二阶精度, 此时s0和s1趋于0.
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.
图4给出了中点A(1,1)和角点B(1.9,1.9)的温度−时间曲线. 可以看到在中点A处, 所有方法均与参考解吻合得较好, 但在边界点B处, Crank-Nicolson方法存在明显的数值振荡.
为了进一步比较这些方法在域内的精度表现, 图5提供了所有方法在t = 0.1时刻温度误差的分布结果, 可以看出本文方法不会像Crank-Nicolson方法一样在边界处产生剧烈的数值振荡, 且在所有方法中, 本文方法的域内精度最高. BDF和Galerkin方法虽然没有数值振荡, 但域内精度不令人满意.
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)的温度−时间曲线, 可以看到本文提出的方法展示出明显的精度优势.
为了展示本文方法的稳定性优势, 另一种情况β = 100, γ = η = 1被考虑. 此时广义梯形法则的Δt = 0.005. 从图7可以看到Crank-Nicolson方法在计算一段时间后就失去了稳定性, 而本文方法不仅没有失去稳定性, 且在所有收敛的方法中具有最高的精度.
5.3 具有外部热源的二维热传导问题
考虑一均质矩形板[9], 如图8所示, 外部热源f(t) = cos(10t)从右边界持续不断地向域内输入, 底部边界施加T = 0的温度场, 其余边界为绝热状态. 利用50 × 50个4结点矩形单元进行划分, 广义梯形法则的Δt =0.005.
图9给出了点A(0.98, 0.98)处的温度和温度一阶导数的结果, 从绝对误差的数值上可以看出: 本文方法要比二阶Crank-Nicolson方法和BDF2更加精确. 总体来说, 这3种二阶精度算法的精度明显优于一阶精度的Galerkin方法和BDF.
![]()
图 9 点A处温度和其一阶导数与时间的关系曲线Figure 9. Temperature and its derivative -time histories of point A6. 结论
本文提出了一种适用于求解一般瞬态热传导问题的无条件稳定单步时间积分方法. 该方法具有无条件稳定性、二阶精度、L型数值耗散和自启动特性. 需要强调的是, 与大多数现有时间积分方法不同, 本文方法对线性和非线性瞬态热传导问题均是无条件稳定的.
线性和非线性数值比较结果表明, 与著名的二阶Crank-Nicolson方法相比, 在计算量相同的前提下, 本文方法具有更理想的精度和稳定性, 并且不会在高频段诱发出虚假的数值振荡. 与具有L型耗散的二阶BDF2相比, 本文方法具有更高的精度且无需其他算法来启动计算.
鉴于本文方法具有优秀的数值性能和易执行性, 故推荐其用于一般瞬态热传导问题的求解.
-
[1] Bathe KJ. Finite Element Procedures. Upper Saddle River: Prentice Hall, 1996.
[2] Fung TC. Weighting parameters for unconditionally stable higher-order accurate time step integration algorithms. Part 1: First-order equations. International Journal for Numerical Methods in Engineering, 1999, 45: 941-970 doi: 10.1002/(SICI)1097-0207(19990720)45:8<941::AID-NME612>3.0.CO;2-S
[3] Fung TC. Higher-order accurate least-squares methods for first-order initial value problems. International Journal for Numerical Methods in Engineering, 1999, 45: 77-99 doi: 10.1002/(SICI)1097-0207(19990510)45:1<77::AID-NME579>3.0.CO;2-5
[4] Fung TC. Solving initial value problems by differential quadrature method–Part 1: First-order equations. International Journal for Numerical Methods in Engineering, 2001, 50: 1411-1427 doi: 10.1002/1097-0207(20010228)50:6<1411::AID-NME78>3.0.CO;2-O
[5] Kim W, Reddy JN. A new family of higher-order time integration algorithms for the analysis of structural dynamics. Journal of Applied Mechanics-Transactions of the ASME, 2007, 84(7): 071008
[6] Kim W, Reddy JN. Effective higher-order time integration algorithms for the analysis of linear structural dynamics. Journal of Applied Mechanics-Transactions of the ASME, 2007, 84(7): 071009
[7] Zhong WX, Williams FW. A precise time step integration method. Part C: Journal of Mechanical Engineering Science, 1994, 208(6): 427-430 doi: 10.1243/PIME_PROC_1994_208_148_02
[8] Wu F, Gao Q, Zhong WX. Fast precise integration method for hyperbolic heat conduction problems. Applied Mathematics and Mechanics (English Edition)
, 2013, 34(7): 791-800 doi: 10.1007/s10483-013-1707-6 [9] Li QH, Chen SS, Kou GX. Transient heat conduction analysis using the MLPG method and modified precise time step integration method. Journal of Computational Physics, 2011, 230: 2736-2750 doi: 10.1016/j.jcp.2011.01.019
[10] 张文志, 富明慧, 蓝林华. 两端边值问题的通用精细积分法. 中山大学学报(自然科学版), 2010, 49(6): 15-19 (Zhang Wenzhi, Fu Minghui, Lan Linhua. General precise integration method for two-point value problems. Acta Scientiarum Baturalium Universitatis Sunyatseni, 2010, 49(6): 15-19 (in Chinese) [11] 彭海军, 高强. 线性非齐次常微分方程两端边值问题精细积分法. 大连理工大学学报, 2010, 50(4): 475-480 (Peng Haijun, Gao Qiang. Precise integration method for two-point boundary value problems of linear nonhomogeneous ordinary differential equations. Journal of Dalian University of Technology, 2010, 50(4): 475-480 (in Chinese) doi: 10.7511/dllgxb201004002 [12] 谭述君, 钟万勰. 非齐次动力方程Duhamel项的精细积分. 力学学报, 2007, 39(3): 374-381 (Tan Shujun, Zhong Wanxie. Precise integration method for Duhamel terms arising from non-homogenous dynamic systems. Chinese Journal of Theoretical and Applied Mechanics, 2007, 39(3): 374-381 (in Chinese) doi: 10.3321/j.issn:0459-1879.2007.03.011 [13] 汪梦甫, 周锡元. 结构动力学方程的更本文精细积分方法. 力学学报, 2004, 36(2): 191-195 (Wang Mengfu, Zhou Xiyuan. Renewal precise time step integration method of structural dynamic analysis. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(2): 191-195 (in Chinese) doi: 10.3321/j.issn:0459-1879.2004.02.009 [14] Kennedy CA, Carpenter MH. Diagonally implicit Runge-Kutta methods for stiff ODEs. Applied Numerical Mathematics, 2019, 146: 221-244 doi: 10.1016/j.apnum.2019.07.008
[15] Boom PD, Zingg DW. Optimization of high-order diagonally-implicit Runge–Kutta methods. Journal of Computational Physics, 2018, 1(15): 168-191
[16] Issa S, Wallin M, Ristinmaa M, et al. Diagonally implicit Runge–Kutta (DIRK) integration applied to finite strain crystal plasticity modeling. Computational Mechanics, 2018, 62: 1429-1441 doi: 10.1007/s00466-018-1572-y
[17] Dong S. BDF-like methods for nonlinear dynamic analysis. Journal of Computational Physics, 2010, 229(8): 3019-3045 doi: 10.1016/j.jcp.2009.12.028
[18] Gear CW. Numerical Initial Value Problems in Ordinary Differential Equations. New Jersey: Prentice Hall; 1971.
[19] Zhang HM, Zhang RS, Masarati P. Improved second-order unconditionally stable schemes of linear multi-step and equivalent single-step integration methods. Computational Mechanics, 2021, 67: 289-313 doi: 10.1007/s00466-020-01933-y
[20] Zhang J. A-stable two-step time integration methods with controllable numerical dissipation for structural dynamics. International Journal for Numerical Methods in Engineering, 2020, 121: 54-92 doi: 10.1002/nme.6188
[21] Zhang J. A-stable linear two-step time integration methods with consistent starting and their equivalent single-step methods in structural dynamics analysis. International Journal for Numerical Methods in Engineering, 2021, 122: 2312-2359 doi: 10.1002/nme.6623
[22] Malakiyeh MM, Shojaee S, Bathe KJ. The Bathe time integration method revisited for prescribing desired numerical dissipation. Computers & Structures, 2019, 212: 289-298
[23] Noh G, Bathe KJ. The Bathe time integration method with controllable spectral radius: the ρ∞-Bathe method. Computers & Structures, 2019, 212: 299-310
[24] Li JZ, Li XY, Yu KP. Enhanced studies on the composite sub-step algorithm for structural dynamics: the Bathe-like algorithm. Applied Mathematics Modelling, 2020, 80: 33-64 doi: 10.1016/j.apm.2019.11.033
[25] Wen WB, Deng SY, Wang NB, et al. An improved sub-step time-marching procedure for linear and nonlinear dynamics with high-order accuracy and high-efficient energy conservation. Applied Mathematics Modelling, 2021, 90: 78-100 doi: 10.1016/j.apm.2020.08.068
[26] Ji Y, Xing YF. An optimized three-sub-step composite time integration method with controllable numerical dissipation. Computers & Structures, 2020, 231: 106210
[27] Li JZ, Yu KP, He HN. A second-order accurate three sub-step composite algorithm for structural dynamics. Applied Mathematics Modelling, 2020, 77: 1391-1412 doi: 10.1016/j.apm.2019.08.022
[28] Xing YF, Ji Y, Zhang HM. On the construction of a type of composite time integration methods. Computers & Structures, 2019, 221: 157
[29] 朱强华, 杨恺, 梁钰等. 基于特征正交分解的一类瞬态非线性热传导问题的新型快速分析方法. 力学学报, 2020, 52(1): 124-138 (Zhu Qianghua, Yang Kai, Liang Yu, et al. A novel fast algorithm based on model order reduction for one class of transient nonlinear heat conduction problem. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(1): 124-138 (in Chinese) doi: 10.6052/0459-1879-19-323 [30] 李艾伦, 傅卓佳, 李柏伟等. 含肿瘤皮肤组织传热分析的广义有限差分方法. 力学学报, 2018, 50(2): 1198-1205 (Li Ailun, Fu Zhuojia, Li Bowei, et al. Generalized finite difference method for bioheat transfer analysis on skin tissue with tumors. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 1198-1205 (in Chinese) [31] 刘硕, 方国东, 王兵等. 近场动力学与有限元方法耦合求解热传导问题. 力学学报, 2018, 50(2): 339-348 (Liu Shuo, Fang Guodong, Wang Bing, et al. Study of thermal conduction problem using coupled peridynamics and finite element method. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 339-348 (in Chinese) doi: 10.6052/0459-1879-17-332 [32] Hughes TJR. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. New Jersey: Prentice Hall, 1987.
[33] Hughes TJR. Unconditionally stable algorithms for nonlinear heat conduction. Computer Methods in Applied Mechanics and Engineering, 1977, 10: 135-139 doi: 10.1016/0045-7825(77)90001-9
[34] Wang YZ, Maxam D, Tamma KK, et al. Design/analysis of GEGS4-1 time integration framework with improved stability and solution accuracy for first-order transient systems. Journal of Computational Physics, 2020, 422: 109763 doi: 10.1016/j.jcp.2020.109763
[35] Wang BL, Tian ZH. Application of finite element-finite difference method to the determination of transient temperature field in functionally graded material. Finite Elements in Analysis and Design, 2005, 41: 335-349 doi: 10.1016/j.finel.2004.07.001
[36] Chen BS, Gu YX, Guan ZQ, et al. Nonlinear transient heat conduction analysis with precise time integration method. Numerical Heat Transfer Part B-Fundamentals, 2001, 40: 325-341 doi: 10.1080/104077901317091712
-
期刊类型引用(2)
1. 冀佳. 建筑电气设备铝合金壳体的焊接技术. 焊接技术. 2023(10): 92-95 . 
百度学术
2. 顾崴，刘铖，安志朋，史东华. 一种基于Hamel形式的无条件稳定动力学积分算法. 力学学报. 2022(09): 2577-2587 . 本站查看
其他类型引用(1)
下载:
计量
- 文章访问数: 898
- HTML全文浏览量: 300
- PDF下载量: 124
- 被引次数: 3
