CN112364544B - 再入气动环境致结构热力响应有限元求解优化方法 - Google Patents

Abstract

本发明公开了一种再入气动环境致结构热力响应有限元求解优化方法，包括：通过有限元方法对基于热传导与材料热弹性动力学耦合控制方程进行离散并给出相应的算法流程；其中，在算法流程中，对于依赖于时间的偏微分方程，有限元方法先对空间区域进行离散，并得到求解区域的网格剖分，然后对时间项进行差分离散，按照迭代耦合松弛计算原理，逐步推进以捕捉服役期满大型航天器离轨再入强气动力热环境致结构材料在空间的振动以及热力响应变形非线性行为。本发明提供一种基于有限元方法的热力耦合响应求解的优化方法，有利于分析与研究材料结构在承受强气动力/热环境下的热力耦合响应，有利于开展对飞行器以及航天器结构性能预测与模拟。

Description

再入气动环境致结构热力响应有限元求解优化方法

技术领域

本发明属于飞行器空气动力技术领域，主要涉及力学、热学和材料科学问题，特别是涉及一种基于有限元方法对服役期满航天器离轨再入气动力热环境致结构变形热力耦合响应求解的优化方法。

背景技术

近地轨道运行的大型航天器服役期满或寿命末期，往往采用离轨控制，使其脱离当前运行轨道而自然陨落，属于无控飞行状态。该类大型航天器在轨任务完成后，将面临再入坠毁、安全评估预示、陨落飞行航迹确定问题。这类航天器再入大气层时所遇强气动力/热环境，使航天器本身承受严重的热流与过载，再入稠密大气环境超高速致高温热化学非平衡气流对航天器复合材料热解烧蚀累积效应、金属桁架结构热力耦合响应、变形软化熔融，会使航天器结构失效解体。因此需要发展适于航天器再入气动力热环境导致的结构热力耦合响应、变形软化熔融与热解烧蚀问题预测分析计算方法。

气动力是指大气压力和表面摩擦力，分别对飞行器产生升力和阻力等力/矩，而气动热直接为结构所感受成为热载荷。气动热使结构材料的力学性能降低，伴随应力作用以致发生蠕变，而结构部件之间的相互约束，在热载荷作用下，又将在结构中产生位移变形，从而使变形加剧并造成翘曲和蠕变特性的变化，同时温度的交替变化也会激起结构的热振动以至颤振。这些情况表明，热结构不仅关系到力学问题，也关系到热学和材料科学问题。金属材料的传热计算通过固体热传导方程进行求解，而结构的力学响应则需要利用计算固体力学等方法进行分析与模拟，外部的高超声速气动力/热通过边界条件传入结构内部，对于高温导致的材料热解烧蚀与金属桁架结构变形软化熔融失效问题，则需要分析材料的相变以及结构响应的非线性行为。

发明内容

本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。

为了实现根据本发明的这些目的和其它优点，提供了一种再入气动环境致结构热力响应有限元求解优化方法，包括：

通过有限元方法对热力耦合控制方程进行离散并给出相应的算法流程；

其中，在算法流程中，对于依赖于时间的偏微分方程，有限元方法先对空间区域进行离散，并得到求解区域的网格剖分，然后对时间项进行差分离散，按照迭代耦合松弛计算原理，逐步推进以得到材料内部的瞬态热传导参数，材料外表面位移膨胀量以及温度变化参数，进而捕捉材料在空间的振动以及热力响应非线性行为。

优选的是，所述算法流程被配置为包括：

S1、基于材料的相应参数、参考温度T₀、Newmark方法参数η、ω，以在确定求解区域后实现网格剖分，进而形成对应的有限元空间；

S2、设置时间积分区间

以及时间积分步N_t，以得到对应的时间步长

S3、设置初始迭代步n从0开始；

S4、分别对3N_v×3N_v整体弹性质量M_u与刚度矩阵K_u，N_v×N_v经修正的整体热质量矩阵M_θ与刚度矩阵K_θ，以及3N_v×N_v耦合矩阵L分别进行组装，以得到相应的整体矩阵；

S5、组装初始整体载荷向量F⁰，确定初始情况下的位移d⁰、速度

温度增量ξ⁰所对应的变量值，以得到第0步的初始加速度

S6、通过迭代计算得到第n+1项的位移、温度增量、热流、应变和热应力；

S7、对迭代步数的值是否收敛进行判断，若是则输出结果并退出，否则返回至S4。

优选的是，在S1中，通过在区域Ω上进行网格剖分，形成有限元空间V_h，其中h表示网格尺寸，则刚度项在有限元空间V_h中的离散形式表示为：

对于温度梯度，在代入离散温度增量表示式θ_h(x,t)后，表达式为

在有限元空间V_h中，将动力学方程中的耦合项

热传导方程的耦合项

与代入热力耦合方程的弱形式进行推导整理，由测试函数

与

的任意性，得到耦合弱形式的空间离散方程为：

优选的是，在S2中，时间导数项的离散推导，是将时间积分区别

划分为N个等距离区间，对于热弹性动力学方程，使用Newmark隐式方法来进行时间上的离散推进，在取ω＝1/2、η＝1/4的情况下，得到可保持格式无条件稳定的以下表达式：

对于耦合热传导方程，使用高精度的Crank-Nicolson格式得到以下的瞬态热传导方程表达式：

将热弹性动力学方程和瞬态热传导方程统一写成矩阵形式为：

优选的是，在S4中，所述整体矩阵被配置为：

优选的是，在S5中，所述初始加速度

被配置为采用以下表达式得到：

优选的是，在S6中，所述迭代计算包括：

S61、组装第n+1步向量Fⁿ⁺¹和Gⁿ⁺¹，以得到整体载荷向量：

S62、开始迭代计算，得到第n+1项的位移dⁿ⁺¹和温度增量ξⁿ⁺¹，并基于第n+1步的速度

加速度

计算相应的热流

应变

与热应力

优选的是，所述

的表达式分别为：

本发明至少包括以下有益效果：本发明通过发展热力耦合响应有限元计算方法，有利于分析与研究材料结构在承受强气动力/热环境下的热力耦合响应非线性行为，有利于开展对飞行器以及航天器结构性能预测与模拟，进而能够深入航天器结构内部，利用非线性有限元方法探索航天器再入阶段强气动力热环境致材料结构变形软化与熔融失效毁坏行为，最终建立起流/热/固一体化耦合计算方法，将具有重要的现实意义。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1为本发明的热力耦合响应有限元算法分类示意图；

图2为本发明的瞬态热力耦合问题计算模型与网格划分示意图(其中，结点个数231、单元个数400)；

图3为本发明在不同耦合状态下，各时刻的温度分布变化情况示意图；

图4为本发明在两种耦合状态r＝1、1.5，材料外表面的热膨胀随时间变化情况计算比较的示意图；

图5为本发明在不同耦合状态下，各时刻应力σ_y沿材料厚度y方向变化情况的示意图；

图6为本发明在热力耦合状态r＝1.5时，各时刻区域变形图(其中，位移放大尺度为M/β)；

图7为本发明在考虑振动的热力耦合状态r＝1.5时，材料上边界位移u₂变化示意图；

图8为本发明在热力耦合状态r＝1.5时，各时刻温度分布图；

图9为本发明在材料底部(y＝0)处的温度变化示意图；

图10为本发明在不同耦合状态下，包含振动的直接耦合法计算比较的示意图；

图11为本发明的四种热力耦合有限元算法计算比较示意图；

图12为本发明的四种热力耦合有限元方法，在不同时刻沿厚度方向应力σ_y计算比较中的其中一部分示意图；

图13为本发明的四种热力耦合有限元方法，在不同时刻沿厚度方向应力σ_y计算比较中的另外一部示意图；

图14为本发明在设置不同阻尼系数条件下，垂直位移以及温度振荡行为计算比较示意图；

图15为本发明的GKUA计算再入近连续过渡区(Kn_∞＝0.01、Ma_∞＝8.37)的竖直平板外部流场云图；

图16为本发明的平板结构热力耦合问题计算区域与网格划分局部图(其中，网格总结点数906，总单元数1500)；

图17为本发明5s时平板内温度增量分布的等值线云图；

图18为本发明平板内部温度与外部温度关联变化计算曲线图；

图19为本发明使用直接耦合法计算得到各时刻的平板变形图及稳态变形图(其中，变形放大尺度为50倍)；

图20为本发明在Kn_∞＝0.01,Ma_∞＝8.37情况下，稳态绕流平板结构场的应力分布图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

基于材料热弹性动力学与热传导方程的热力耦合模型的弱形式，构造设计了四种热力耦合状态下问题的有限元计算格式，建立了一套不考虑物体弹性振动与结构振动耦合瞬态热传导影响的有限元耦合算法。通过对一维瞬态热传导所产生的结构耦合问题计算检验，分析了在是否考虑振动的情况下，直接耦合法与顺序耦合法之间的相互区别与联系，因为Newmark算法在模拟位移动力学问题上表现出的卓越性能，使得本发明所设计的有限元算法能准确捕捉材料在空间的振动以及衰减行为。通过考虑振动的一维瞬态热传导耦合计算结果，证实了本发明提出的数值方法用于计算物体在温度与应力相互耦合情况下的温度分布与位移变化的可靠性。

在验证本发明热力耦合响应有限元计算方法正确性基础上，结合求解Boltzmann模型方程统一算法外部高超声速流场计算与本发明结构热力耦合问题有限元算法，开展了二维竖直平板在外部再入高超声速极端力热环境下的结构力-热耦合响应行为一体化计算，分析表明，在热力响应耦合计算过程中，平板左侧迎风区的物面温度、压力较来流值增加大于右侧背风区物面温度、压力，平板在此热力耦合双重作用下，产生变形，致使向平板背风面方向弯曲，揭示了飞行器在强气动力热绕流环境作用下，致物体发生弯曲变形的非定常变化过程与影响机制。

为了发展适于航天器再入气动力热环境致结构变形热力耦合响应有限元计算方法在一维、二维结构场的初步研究，有关三维结构场的算法应用还有待下一步研究展开，体会到，流场与结构场的耦合分析与计算是一个复杂过程，为简单起见，在利用求解Boltzmann模型方程的统一算法计算得到流场边界关于温度与压力等数据的基础上，开展结构的耦合响应问题模拟。计算中对于流场边界的处理是将边界的温度与压力值在各时间步上均取常值，这是我们进行流固耦合一体化计算的初始尝试，与真实情况有一定差距，如热力耦合方程边界条件所描述的，应该考虑的是材料外部复杂的热流边界条件，此时将流场边界的温度作为外部温度T_ext，同时考虑流场与结构的对流换热以及外部环境的辐射换热问题，由此产生的温度分布以及结构变形受力情况则更加符合实际，成为未来进一步研究发展的方向。

实施例：

在热弹性动力学方程中，设在空间区域

中，不考虑阻尼的影响，材料热弹性动力学方程可表示为：

其中，u＝[u₁ u₂ u₃]^T表示位移向量，f＝[f₁ f₂ f₃]^T为体积力向量，上标T表示矩阵转置，ρ为材料密度，τ表示材料阻尼系数，C_ijkl为材料四阶弹性张量，对于均匀各向同性材料，C_ijkl表示为：C_ijkl＝λδ_ijδ_kl+G(δ_ikδ_jl+δ_ilδ_jk).

其中，δ_ij为Kronecker符号，λ与G为材料Lamé常数，用材料杨氏模量E与泊松比ν表示为：

用θ＝T-T₀表示温度增量，T为绝对温度值，T₀为未变形状态下的温度值。β_ij为热弹性模量，可表示为：β_ij＝C_ijklα_kl.；

其中，α_kl为材料热膨胀系数，对于各向同性材料，α_kl＝α₀δ_kl，物体的热弹性模量可表示为：

其中，

为材料体弹性模量，说明温度的变化使得物体体积发生膨胀或收缩进而影响结构应力变化，物体的应变ε_ij表示为：

热应力σ_ij表示为：σ_ij＝C_ijklε_kl-β_ijθ.。

在瞬态热传导方程中，考虑位移场对温度场的影响，则耦合热传导方程为：

其中，c表示材料比热，k_ij为材料热传导张量，h为热源项，对于各向同性体，根据β_ij的表达式，方程中的第二项也可以表示为：

这里，ε_ii＝ε₁₁+ε₂₂+ε₃₃为物体体积应变，方程式说明物体的体积应变变化率会对温度场产生影响；

对于动力学方程式，初始条件为区域Ω上的初始位移与初始速度，即

对温度场，由于未知量选取的是温度增量，则初始温度增量为0，所以我们有θ(x,0)＝0.；

设物体的边界为光滑边界

在位移场中，边界条件分为位移边界Γ₁与力边界Γ₂，在位移边界Γ₁上给定位移值，一般取固定边界条件u＝0.。

在力边界Γ₂上，物体承受外部面力为p＝[p₁ p₂ p₃]^T，对于热弹性物体，边界条件为σ_ijn_j＝(C_ijklε_kl-β_ijθ)n_j＝p_i,

其中，n＝[n₁ n₂ n₃]^T为边界Γ₂单元外法向矢量。两种类型的边界不能重合，即有

温度场边界可分为温度边界Γ′₁与热流边界Γ′₂，Γ′₁上给定温度增量值

即

Γ′₂上的边界条件与热流值相关，如果同时考虑外边界的热流，对流与辐射边界条件则有

其中，q表示边界热流值，T_ext表示物体外部温度值，H为材料对流换热系数，ε是材料表面辐射发射率，σ＝5.67×10^-8W·m^-2·K^-4是斯特藩-玻尔兹曼常数。为简单起见，我们只考虑热流边界条件，即：

同样地，两种类型的边界不能重合，即：

利用变分原理推导热力耦合方程的弱形式，定义如下函数空间

与

其中，H¹(Ω)为函数及其一阶导数平方可积(即属于L²(Ω))的Sobolev(索伯列夫)函数空间。位移u(x,t)可看作如下映射：

这表示对任意时刻t，

类似地，温度增量θ(x,t)也看作如下映射：

并且有

对于热弹性动力学方程式分别乘以向量测试函数

的各分量，

并在空间Ω上积分，有

利用散度定理，得到：

结合边界条件与，并考虑到测试函数

在Γ₁满足

利用该条件可得到

同理，定义温度测试函数

将热传导方程式，乘以

在Ω上积分，有

同样利用散度定理，得到：；

代入边界条件，由测试函数

的要求，在边界Γ′₁满足

利用该条件，整理(2.21)式，得到：

综上，关于u与θ的热力耦合方程弱形式为：

由上述积分等式，可以看出温度场与位移场相互依赖，所以在有限元计算过程中两方程需要联立求解。

在实际操作中，计算流体力学是通过数值方法求解流体力学控制方程，得到流场的离散的定量描述，并以此预测流体运动规律的学科。在CFD中，把流体运动控制方程中的积分、微分项近似地表示为离散的代数形式(大概意思是将连续的曲线转化为一个个点)，使得积分或微分形式的控制方程转化为代数方程组；然后，通过电子计算机求解这些代数方程组，从而得到流场在离散的时间/空间点上的数值解(numerical solution)。CFD有时也称流场的数值模拟、数值计算、或数值仿真等。

在流体力学控制方程的微分和积分项中包括时间/空间变量(自变量)以及物理变量(因变量)，这些变量分别对应着时间/空间求解域和定义在求解域上的流动问题的解。要把这些积分或者微分项用离散的代数形式代替，必须首先把求解域表示为离散形式。在不同的离散方法中，求解域或者被近似为一系列网格点(grid points)的集合，或者被划分为一系列控制体或单元体(control volume,cell)。

因变量定义在网格点上或者控制体的中心、顶点或其他特征点上。在每一个网格点或者控制体上，流体运动方程中的积分或者微分项被近似地表示为离散分布的因变量和自变量的代数函数，并由此得到作为微分或积分型控制方程近似的一组代数方程，这个过程称为控制方程的离散化(discretization)，其中所采用的离散化方法称为数值方法或者数值格式。

这组代数方程的解(即数值解)给出了离散点上流场的定量描述。显然，为了得到流场结构的比较精细的描述，网格点或者单元体的数量必须足够多。对于科学和工程中的一些常见的流动问题，网格的数目常常需数万到数百万甚至更多(有人曾经使用十亿以上的网格进行湍流的直接数值模拟)。因此，不借助计算机，我们就无法有效地求解这些代数方程组。在实用中，人们采用各种程序设计语言求解过程编制成计算机程序，通过在计算机上运行这些程序得到数值解。

在具体操作过程中，通过利用有限元方法对热力耦合控制方程进行离散并给出相应的算法流程，对于依赖于时间的偏微分方程，有限元方法先对空间区域进行离散，并得到求解区域的网格剖分，然后在时间项利用改进的Euler方法进行差分离散，按照迭代耦合松弛计算原理，逐步推进求解，具体来说包括：

首先，对耦合弱形式的空间离散是使用有限元方法先对空间区域进行离散，得到求解区域的网格剖分，便于计算机进行运算。

在区域Ω上进行网格剖分，形成有限元空间V_h，这里的h表示网格尺寸，设网格结点个数为n_V，在V_h上位移与温度增量的插值函数u_h(x,t)、θ_h(x,t)分别表示为：

其中，d(t)为结点位移向量，其元素d_ij(t)，i＝1,...,n_V，j＝1,2,3表示第i个结点位移的第j个分量,ξ_i(t)表示在结点i处的温度增量值，

为结点形函数。

同样测试函数

与

也可以表示为

分别用上标“·”、“··”表示变量对时间的一阶与二阶导数，则未知量对时间的各阶导数值表示为

在处理动力学方程弱形式的刚度项

中，使用工程应变向量：

ε＝[ε₁₁ ε₂₂ ε₃₃ 2ε₁₂ 2ε₂₃ 2ε₁₃]^T.

用位移表示应变向量为

这里，偏导数矩阵算子

定义为：

代入离散位移表达式，得到ε_h(x,t)＝B(x)d(t)；

其中，

为应变矩阵。

由式，各向同性体的弹性系数C_ijkl可写成如下矩阵D的形式:

于是刚度项在有限元空间V_h中的离散形式表示为

对于温度梯度

其中，向量a定义为:a＝[1 1 1 0 0 0]^T；

代入离散温度增量表示式θ_h(x,t)，有

其中，

为梯度矩阵。

对于动力学方程中的耦合项

在有限元空间V_h中，ψ_i取所有的形函数

代入θ_h(x,t)式、

式以及β_ij＝βδ_ij，得到：

对于热传导方程的耦合项

代入位移u_h(x,t)与

的表达式，有

将等式、与代入热力耦合方程的弱形式中，推导整理，由测试函数

与

的任意性，我们最终得到耦合弱形式的空间离散方程为：

各矩阵表达式分别为：

这里，M_u与K_u分别为3N_v×3N_v整体弹性质量与刚度矩阵，M_θ与K_θ分别为N_v×N_v经修正的整体热质量与刚度，因为除以了参考温度T₀。C_u为3N_v×3N_v整体阻尼矩阵，并与能量的耗散相关。在工程上，C_u常取为所谓的Rayleigh阻尼，即

C_u＝τ₁M_u+τ₂K_u,；

其中，τ₁与τ₂为参数，一般通过实验得到，表示阻尼是粘性阻尼与固体阻尼之和。L称为3N_v×N_v耦合矩阵。

与(Γ′₂)^h分别是边界Γ₂与Γ′₂在有限元空间V_h的离散逼近。

其次，引入Newmark方法和Crank-Nicolson格式，Newmark隐式方法作用：进行时间上的离散推进，即从第n步推进至n+1步；Crank-Nicolson格式作用：采用的离散化方法，将连续值转化为适于计算机运算的离散值因为Newmark算法在模拟位移动力学问题上表现出的卓越性能，使得所设计的有限元算法能准确捕捉材料在空间的振动以及衰减行为。

上节已将材料热弹性与热传导耦合控制方程的弱形式在空间区域离散化处理，推导得到耦合常微分方程组。下面进行时间导数项的离散推导，将

划分为N个等距离区间：

其中，时间步长

t_n＝nΔt，n＝0,1,2,...,N_t，t_n时刻各未知量用上标表示为：

对于热弹性动力学方程，使用在求解动力学方程中应用最为广泛的一种Newmark隐式方法来进行时间上的离散推进，该方法采用以下公式计算速度与位移：

其中，参数0≤ω≤1，0≤η≤1/2，得到n+1时刻的加速度为：

对在n+1时刻的动力学方程：

代入dⁿ⁺¹、

整理得到：

通常取ω＝1/2、η＝1/4，可保持格式的无条件稳定。

在起始时间步上需要知道d⁰、

初始位移与速度d⁰和

以及初始温度增量ξ⁰已知，可由起始方程

计算得到

对于耦合热传导方程，使用高精度的Crank-Nicolson格式，有

最终整理得到

这样，将热弹性动力学方程和瞬态热传导方程统一写成矩阵形式为：

这是一般情况下热力耦合方程的时间推进计算格式，工程上称之为直接耦合法(DC，Direct Coupling)，根据所考虑问题的不同，有以下三种简化形式：

(1)不考虑物体的弹性振动，上述迭代格式(2.2.36)简化为

此时对于热弹性控制方程可以进行线性化处理，是一个静态问题，若将上述方程写成增量形式，即

dⁿ⁺¹＝dⁿ+Δdⁿ，ξⁿ⁺¹＝ξⁿ+Δξⁿ，Fⁿ⁺¹＝Fⁿ+ΔFⁿ；

并令Gⁿ＝0，在温度场上使用全隐式格式，我们得到

(2)考虑弹性振动但不考虑变形对温度的影响，格式(2.2.36)简化为：

此时两方程解耦，实际计算时可以单独求解温度增量场，再将温度增量所产生的热应力作为载荷加载到位移场中，这两步分别为

这是工程上遵循的热应力分析过程，通常称之为顺序耦合法(SC，SequentiallyCoupling)或者间接耦合法。

(3)不考虑弹性振动，也不考虑变形对温度的影响，格式简化为：

此时两方程也解耦，由上述格式，实际计算可分两步进行：

K_udⁿ⁺¹＝Lξⁿ⁺¹+Fⁿ⁺¹.

这是工程中最常使用的顺序耦合法，线弹性方程此时也变成了与第(1)种简化情况相似的静态问题。

综上，我们得到了四种热力耦合方程式的有限元计算格式，图1显示了上述推导的四种不同的方法，根据是否考虑变形对温度的影响将算法分成直接耦合法与顺序耦合法，再根据是否考虑弹性振动，进一步将直接耦合法与顺序耦合法细分成包含振动与不包含振动两种形式。

再次，得到热力耦合方程算法流程，基于一般的包含振动的直接耦合法(DC，VibIncl.)给出整个热力耦合问题的计算流程。

1.输入程序初始条件，具体来说是设置材料各参数材料密度ρ、材料比热c、材料杨氏模量E与泊松比ν(或者材料Lamé常数λ与μ)、热弹性模量β(或者材料热膨胀系数α)、参与求解阻尼的参数τ(或者τ₁、τ₂)、参考温度T₀以及Newmark方法参数η、ω。确定求解区域Ω，对区域作单元剖分，得到合适的有限元网格，选择适当的单元类型，形成有限元空间V_h，设置时间积分区间

以及时间积分步N_t，得到时间步长

取初始迭代步从n＝0开始。

2.程序的初始化，即计算第0步的相关数据，通过分别组装矩阵M_u、K_u、L、M_θ、K_θ，并形成计算格式中的整体矩阵：

组装初始整体载荷向量F⁰，确定初始变量值d⁰、

ξ⁰，计算初始加速度

3.迭代计算开始，得到第n+1项的位移dⁿ⁺¹和温度增量ξⁿ⁺¹，通过组装第n+1步向量Fⁿ⁺¹和Gⁿ⁺¹，形成计算格式中的整体载荷向量：

对位移与温度增量施加相应的边界条件，计算得到n+1步dⁿ⁺¹和ξⁿ⁺¹。

其中，d是位移，

是速度，

是加速度。

4.计算所需要输出的相关参数，即第n+1步的输出参数，通过依次计算n+1步速度

与加速度

根据需要计算相应的热流

应变

与热应力

计算表达式依次为：

5.当n+1步的结果收敛，即变化幅度小于一定值，程序计算结束，并输出结果，即令n＝n+1，当n＞N_t时，结束计算，否则返回第3步。

数值算例与结果分析，为验证所提出算法模型的有效性，通过下面的数值算例进行验证。

1、一维瞬态热传导计算检验

如图2(a)，考虑单元厚度的热弹性材料层，在整个区域上处于绝对温度T₀，在层底(y＝0)刚性固定(u＝0)并且绝热，在t＞0时，材料外表面(y＝1)维持在绝对温度T＝T₀+θ₀,即温度增量为θ₀。不考虑结构振动，外表面无面力，区域内无热源与体积力，对该问题的瞬态热力耦合问题进行计算分析，由现有技术看出，该问题具有解析解，温度增量与垂直位移的具体表达式为：

其中，

与

分别是温度增量θ与垂直位移u₂的Laplace变换，

s是Laplace变换的变量，M与A是与绝对温度和等温条件相关的模量，分别为：

可使用反Laplace变换，求得温度增量与垂直位移的值。这里，比值r＝A/M可用来反映两物理场的耦合程度，r越大，耦合程度越高，r＝1时表示两物理场解耦。

由问题区域的特点，在有限元计算中我们使用二维区域来模拟一维瞬态热力耦合问题，三维的动力学方程转化为二维平面应变动力学问题，取材料层区域为：Ω＝[0,0.5]×[0,1]，网格剖分如图2(b)所示，其中结点个数231、单元个数400，

网格单元采用线性三角形单元，由问题设定，材料在x方向无限长，材料底部固定，在边界x＝0上固定位移，在x＝0.5、y＝0.5上分别固定其法向位移，其余边界无面力。在上边界y＝1上设置温度增量θ₀＝1，其余边界绝热。这里我们进行无量纲计算，取材料物性参数为：

ρ＝1.0,c＝1.0,k＝1.0,E＝2.0,ν＝0.3,β＝0.05，

其中，参考温度T₀的值由比值r确定。解耦(r＝1)时T₀取为0，此时温度增量方程的弱形式两边不用除以T₀，直接消去耦合项即可。因为没有结构的振动，我们所使用的有限元算法为不考虑振动的直接耦合法(DC，Vib Excl.)和顺序耦合法(DC，Vib Excl.)，显然r＝1使用顺序耦合法，其余情况使用直接耦合法。

图3(a)-(d)绘出了之前介绍的方法计算不同热力耦合程度r＝1、1.1、1.5条件下，温度增量沿材料厚度y方向各时刻t＝0.1、0.2、0.5、1.0的变化分布情况。图中曲线表示算法结果、符号来自参考现有技术的结果，不同时刻不同热力耦合程度得到的两套结果变化规律完全一致，证实用于计算材料内部瞬态热传导问题，揭示温度变化规律的准确可靠性。分析表明，由于温度的增加，物体产生热膨胀，而由此产生的体积变化导致温度的增加变得滞后，耦合程度越高，这种滞后程度就越明显。随着计算时间推进，瞬态热传导逐渐趋于稳定，材料内外表面温度增量趋于给定值。

为了揭示瞬态热传导致材料外表面位移变化规律，图4给出了两种热力耦合响应状态r＝1、1.5下材料外表面y方向位移随时间变化分布与现有技术计算比较。图中线型表示计算结果、符号为来自现有技术的数据，可看出，计算不同热力耦合度得到的材料外表面位移随时间变化分布与现有技术及精确解均吻合很好，说明本发明的数值方法用于计算瞬态热传导引起材料外表面位移膨胀量是正确有效的。而且随着时间的延长，材料外表面位移膨胀量呈现非线性增大，相比于解耦情况，在全耦合状态下物体的位移也相对产生了延后，这是由于温度升高致延后的表现。

图5(a)-(d)展示了与图5(a)-(d)所对应的三种耦合状态r＝1、1.1、1.5在不同时刻沿厚度y方向应力分布变化情况。由于顺序耦合方程忽略了热力耦合响应引起变形所产生的温度场变化，反映出整体应力水平高于双向耦合情况。随时间推移，应力非线性效应减弱，材料内外表面应力大小趋近，热力响应耦合的程度越高，对应力所产生的滞后也越明显。

结合图3-图5，分析表明，在温度场与位移场相互耦合状态下，物体的温度变化与位移变化以及带来的应变与应力变化相互影响，都产生了一定的滞后效应，这种特性、影响和Yang、Chen[21]的描述相吻合，由此证实了所揭示变化规律的准确性。

2、考虑振动的一维瞬态热传导耦合计算分析

算例中的位移场忽略了结构的振动，下面给出该算例考虑物体弹性振动的耦合形式的计算结果。计算中均使用该算例中的材料物性参数式，只是在位移场中增加了惯性项。由图1，此时使用的有限元方法是包含振动的直接耦合法(DC，Vib Incl.)。

图6绘出了在耦合状态r＝1.5各时刻的结构变形图，从中我们看到，由于初始状态在边界上的温度突然升高，温度增量作为体载荷施加到结构中，使得材料结构产生一个瞬时加速度，从图6(a)与(b)可观察到一个弹性纵波向负y方向传播，在t＝0.5到1.0这个时间段，波前到达区域底部后反弹，波进行回传，与此同时由于温度持续升高，位移u₂不断增加，与上节算例中t＝1时温度场与位移场基本达到平衡的情况不同，此时的位移场仍在振动。图7显示了材料外表面的持续振动曲线，由于考虑了变形对温度的影响，一部分弹性振动能转化成热能，随时间的推移产生了耗散，最后位移场逐渐趋于稳定，最终达到与不考虑振动的间接耦合法相同的稳定状态。

图8给出了各时刻温度的变化情况，受位移场的影响，温度在由外向内增加时也存在波动行为，与图6(a)-(b)中弹性波传播到底部反弹相对应，温度场在t＝0.5到1.0这个时间段也有振荡，这从下文图9中材料底部(y＝0)的温度增量出现的先降后升的变化表现出来，最后温度场在t>1.0后也仍然持续振动并渐近趋于稳态。

前面我们只计算了在耦合状态r＝1.5时位移场与温度场的振动变化情况，下面在图10中用各种线型给出了在各种耦合状态r＝1、1.1、1.5位移场与温度场的振荡关系，从中可以看到，在考虑振动的情况下，耦合程度越低时，变形产生的弹性振动能转化成热能的比例就越小，此时能量耗散得越慢，位移场持续振荡的时间越久，但随着弹性振动能的耗散，位移场与温度场的振幅都逐渐变小，最终趋于稳定。另一方面，当考虑解耦的情况r＝1.0时，温度场的变化与位移场无关，那么此时温度场的计算与不考虑振动的顺序耦合法计算结果相同，但温度增量作为体积力施加在位移场上，使区域产生了持续的振动，在无阻尼的情况下，弹性振动能不发生耗散，那么这种振动则会一直保持下去。

图11用不同线型给出了四种有限元算法的计算结果，从图11(a)中材料表面位移变化曲线看出，除了包含振动的顺序耦合法，在其余三种方法中，位移场最终都会趋于同一稳定状态，在不考虑振动的方法中，位移在趋于稳定前，直接耦合法计算得到的各个时刻的位移值都比顺序耦合法小，达到稳定的时间也比顺序耦合法晚。另一方面，考虑振动与不考虑振动会使计算结果产生明显的差别，一旦产生弹性振动能耗散，位移场都会振荡趋近于稳定状态。对于温度场来说，由图11(b)中各种线型表示的计算结果看出，所有的方法都会趋于一个最终的稳态，本质上是因为热传导过程是一个不可逆的能量耗散过程。

图12-13给出了我们所关心的使用四种热力耦合有限元算法计算得到各个时刻沿厚度方向的应力变化曲线。其中在考虑振动的方法中，弹性波沿底部传播再反弹的现象清楚地反映出来。最终除了考虑振动的顺序耦合法，因为没有弹性振动能的损失会不断振荡下去之外，其他方法计算得到的应力分布都会趋近于同一个状态。

分析上述计算结果，注意到即使在热弹性动力学方程中不存在阻尼矩阵C_u，一旦耦合矩阵L存在，则整个系统也终将会趋于稳态。现在我们将阻尼矩阵加入到系统中以及计算格式中，使得该模型更加接近于实际。为简单起见，我们取C_u的形式为。图14使用不同的线型给出了无量纲在阻尼系数τ＝0,0.1,0.5以及1的条件下，结构外表面以及底部温度的振荡曲线。在图14(a)中的解耦状态下，热致振动将永不衰减，除非阻尼项存在，即τ＞0，而在耦合状态下，由图14(b)与(c)，系统的能量衰减同时来自于弹性动力学方程中的阻尼效应以及热传导方程中的耦合项，使得结果的收敛速度变得更加快。另一方向，虽然引入阻尼效应会增加额外的计算时间，但算法的稳定性得到了一定程度的改进。根据1-2步中的数值算例结果分析，验证了数值方法用于计算瞬态热传导致物体温度、位移变化引起应变、应力相互影响及材料热力耦合响应的正确可靠与适应性，这为下一步发展航天飞行器受到较高的热力耦合响应致物面材料内部温度分布与材料结构变形、热损毁数值模拟方法奠定了基础。

3、高超声速再入气动力热环境平板结构热力耦合响应计算分析，运用本发明建立的热力耦合响应计算方法，分析飞行器结构外表面在承受气动加热时的热传导与热变形响应行为。为简单起见，假设表面的气动力、气动热已由跨流域空气动力学数值方法]计算给出，因此直接把气动热与气动力代入飞行器热结构的热力耦合计算公式，作为边界条件，由此分析结构内部温度分布与受力变形行为。

为了将本发明有限元方法与求解Boltzmann模型方程统一算法(GKUA)相结合，这里我们研究了图15所示无限长竖直平板绕流模型问题，平板厚度0.015m，高度0.5m，平板外部受再入近连续稀薄过渡流区高超声速绕流冲击，来流马赫数Ma_∞＝8.3666，来流克努森数Kn_∞＝0.01，γ＝1.4。使用现有技术提出含转动非平衡效应Boltzmann模型方程统一算法求解外流场，因再入高超声速近连续过渡流，在竖直平板附近绕流场中形成了厚厚的脱体激波层，如图15(b)与(c)，平板表面由此受到强气动热与气动力作用，在迎风区承受高温与高压，若不考虑平板的振动，根据GKUA求解获得平板表面温度与压力作为本发明材料结构热力耦合算法界面条件，对平板内部的温度分布与变形情况进行计算与分析。

对平板截面内部结构进行网格剖分，得到如图16(b)的计算网格。固定平板下表面(y＝0)。由于用于绕流物体外部流场计算的网格与用于本发明热力耦合响应计算的平板内部结构网格的不一致性，这里将GKUA计算得到的外部流场中网格边界的温度值与压力值进行线性插值施加到平板内部热力耦合算法所用网格边界中，实现外部流场计算与平板内部热力响应的耦合计算。假设平板处于初始温度T₀，外部流场处于稳定状态，平板外表面温度达到T＝T₀+θ，其中T为GKUA计算竖直平板外流场稳定时获得的温度值，由此得到平板物面边界的温度增量值θ＝T-T₀。

在外部流场计算中，因飞行器再入不同高度跨流域绕流多尺度效应，通常采用无量纲模型化计算，而在物体内部结构场的热力耦合有限元计算中则需将其转换成更符合实际的有量纲温度与压力值。例如在流场中计算得到的温度与压力值分别为T^*与p^*，则实际温度与压力值T与P分别为：

T＝T^*T_∞，

其中，T_∞、ρ_∞、V_∞分别为参考温度、来流密度与来流速度。在热力耦合计算中选择平板为钢板，材料物性参数为：

E＝2.06×10¹¹Pa,ν＝0.3,α＝1.5×10^-6/℃

ρ＝7850kg/m³,c_v＝460J/(kg·℃),k＝50W/(m·℃)

参考温度T₀＝500℃，取时间步Δt＝0.05s，计算时间为5s。

图17(a)、(b)分别给出了使用不含振动效应的顺序耦合法与直接耦合法计算得到平板内部结构场在5s时的温度分布。可以看出，使用直接耦合法，温度还未达到完全稳态，而顺序耦合法计算的温度场在5s时已经到达了稳态。图18绘出了区域不同的部位外部流场与平板内部温度场相关联变化曲线，由图看出，平板外部使用求解Boltzmann模型方程统一算法计算得到迎风面y＝0.06、0.26、0.46处温度与背风面温度，均分别与本发明平板内部热力耦合计算得到的平板表面温度吻合一致，证实本发明关于物体内部热力响应耦合算法模型与程序设计实现的正确性。从图17绘出的平板内部温度分布看出，物体内部靠近平板迎风面驻点(坐标原点)的温度比附近其他点的温度高，平板表面最高的温度点则出现在平板的左上拐点，此处边界曲率变化最大，流场热流在此达到最大，由稳态热传导方程的极大值原理可知，材料内部温度在此也达到最大,故在平板的迎风边界由下到上温度出现一个先降后升的变化，如图18分别绘出y＝0.06、0.26、0.46三条曲线所示，这种变化在材料内部也清楚地反映出来。

图19给出了使用直接耦合法计算得到的平板各时刻的整体变形情况，可以看出，在受外部高超声速强气动力热绕流迎风面与背风面较大温度与压力差双重作用下，平板首先产生了较大幅度的弯曲(图19(a))，这是迎风面压力的作用结果，随后气动热使得平板内部温度逐渐升高，在内部形成了非均匀的热应力分布，高温使得平板在迎风面产生较背风面更大的热膨胀，所产生的热应力方向与迎风面压力方向相反，内外两种力相互作用使得平板弯曲幅度变得相比于计算初始阶段的弯曲幅度更小，如前所述，在5s时，使用直接耦合法计算的平板内部温度分布还未稳定，所以热膨胀仍在持续，即5s时平板的弯曲变形还将继续，最终达到图19(f)的稳态变形分布。两种方法得到的最终稳态分布是一致的，根据所取材料物性参数，对该算例定义的热力耦合度比值为r＝M/A≈1.034，可以看出，两种耦合情况所得到的计算结果较为一致、相差甚微，使用直接耦合法能更真实反应结构在各时刻的热力耦合状态，但若只关注结构最终的稳定状态，使用顺序耦合法进行热应力分析也是合理的，且以其快速与高效在工程中广泛使用。

图20给出Kn_∞＝0.01,Ma_∞＝8.37稳态绕流平板内结构应力分布，可以看到正应力σ_x与σ_y都与温度分布相关。最大剪应力及Von Mises应力(如图20(c)、(d)与(e))位于平板底部背风区处。由图19对于变形的计算分析，我们认为由于来流的压力作用，平板迎风区的热膨胀受到了抑制，而在平板背风区承受了最大的温差并且外部压力最小，使得该区域更容易产生非线性变形而失效。

以上方案只是一种较佳实例的说明，但并不局限于此。在实施本发明时，可以根据使用者需求进行适当的替换和/或修改。

这里说明的设备数量和处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言，可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (7)

1.一种再入气动环境致结构热力响应有限元求解优化方法，其特征在于，包括：

通过有限元方法对再入气动环境导致的结构热力耦合控制方程进行离散并给出相应的算法流程；

其中，在算法流程中，对于依赖于时间的偏微分方程，有限元方法先对空间区域进行离散，并得到求解区域的网格剖分，然后对时间项进行差分离散，按照迭代耦合松弛计算原理，逐步推进以得到材料内部的瞬态热传导参数，材料外表面位移膨胀量以及温度变化参数，进而捕捉材料在空间的振动以及热力响应非线性行为；

所述算法流程被配置为包括：

S1、基于材料的相应参数、参考温度T₀、Newmark方法参数η、ω，以在确定求解区域后实现网格剖分，进而形成对应的有限元空间；

S2、设置时间积分区间

以及时间积分步N_t，以得到对应的时间步长

S3、设置初始迭代步n从0开始；

S4、分别对3N_v×3N_v整体弹性质量M_u与刚度矩阵K_u，N_v×N_v经修正的整体热质量矩阵M_θ与刚度矩阵K_θ，以及3N_v×N_v耦合矩阵L分别进行组装，以得到相应的整体矩阵；

S5、组装初始整体载荷向量F⁰，确定初始情况下的位移d⁰、速度

温度增量ξ⁰所对应的变量值，以得到第0步的初始加速度

S6、通过迭代计算得到第n+1项的位移、温度增量、热流、应变和热应力；

S7、对迭代步数的值是否收敛进行判断，若是则输出结果并退出，否则返回至S4。

2.如权利要求1所述的再入气动环境致结构热力响应有限元求解优化方法，其特征在于，在S1中，通过在区域Ω上进行网格剖分，形成有限元空间V_h，其中h表示网格尺寸，则刚度项在有限元空间V_h中的离散形式表示为：

对于温度梯度，在代入离散温度增量表示式θ_h(x,t)后，表达式为

在有限元空间V_h中，将动力学方程中的耦合项

热传导方程的耦合项

与代入热力耦合方程的弱形式进行推导整理，由测试函数

与

的任意性，得到耦合弱形式的空间离散方程为：

3.如权利要求1所述的再入气动环境致结构热力响应有限元求解优化方法，其特征在于，在S2中，时间导数项的离散推导，是将时间积分区别

划分为N个等距离区间，对于热弹性动力学方程，使用Newmark隐式方法来进行时间上的离散推进，在取ω＝1/2、η＝1/4的情况下，得到可保持格式无条件稳定的以下表达式：

对于耦合热传导方程，使用高精度的Crank-Nicolson格式得到以下的瞬态热传导方程表达式：

将热弹性动力学方程和瞬态热传导方程统一写成矩阵形式为：

4.如权利要求1所述的再入气动环境致结构热力响应有限元求解优化方法，其特征在于，在S4中，所述整体矩阵被配置为：

5.如权利要求1所述的再入气动环境致结构热力响应有限元求解优化方法，其特征在于，在S5中，所述初始加速度

被配置为采用以下表达式得到：

6.如权利要求1所述的再入气动环境致结构热力响应有限元求解优化方法，其特征在于，在S6中，所述迭代计算包括：

S61、组装第n+1步向量Fⁿ⁺¹和Gⁿ⁺¹，以得到整体载荷向量：

S62、开始迭代计算，得到第n+1项的位移dⁿ⁺¹和温度增量ξⁿ⁺¹，并基于第n+1步的速度

加速度

计算相应的热流

应变

与热应力

7.如权利要求6所述的再入气动环境致结构热力响应有限元求解优化方法，其特征在于，所述

的表达式分别为：

Images