一种高超声速流动-传热与结构响应的多场耦合瞬态数值的方法

技术领域

本发明涉及现代高速飞行器设计技术领域，特别涉及一种高超声速流动-传 热与结构响应的多场耦合瞬态数值的方法。

背景技术

高超声速飞行器的快速发展给热防护设计带来了更为严峻的挑战。准确的预 测气动热/力环境、结构温度和应力状态，能够在提高飞行器安全性能的同时减 小热防护系统设计冗余，对提高飞行器的性能有着极为重要的意义。

在现有技术中，传统的高超声速飞行器的热/力载荷环境的预测与热防护结 构性能的分析基本还处于分离状态。现有技术中的方法一般都是首先在给定的 等温壁面条件下进行流体计算，得到壁面热流或传热系数；然后将得到的热载 荷作为边界条件加载到结构上进行热分析得到固体热分布；最后再根据该固体 热分布计算得到结构的应力和应变。

现有技术中的上述方法实际上是把多物理场耦合的事实人为的分割成多个 独立的物理场，并且也没有考虑各个物理场之间的相互作用。因此，在这种情 况下，既无法得到精确的气动热/力载荷环境，也无法正确地评价热防护材料及 其结构的服役特征。

高超声速飞行热防护设计是一个涉及到真实气体效应、耦合传热和结构热力 响应的复杂的多物理场耦合问题，必须采用多场耦合的方法求解。但是，由于 多物理场耦合问题的复杂性，还需要进一步开展分析方法研究，深刻把握防热 系统多场耦合规律及其效应。

发明内容

有鉴于此，本发明提供一种高超声速流动-传热与结构响应的多场耦合瞬态 数值的方法，从而可以实现高超声速非平衡流动求解器与结构热/力全耦合求解 器相耦合的多场耦合计算，使得对高超声速飞行器的气动热力环境和结构热力 响应的预测更符合物理实际，并可保证计算精度。

本发明的技术方案具体是这样实现的：

一种高超声速流动-传热与结构响应的多场耦合瞬态数值的方法，该方法包 括：

A、预先建立多物理场耦合模型并设置当前边界条件；

B、根据结构确定壁面温度和位移边界条件，在流体-固体耦合界面进行 数据交换，得到流体区域的当前温度和位移边界条件；

C、根据所述流体区域的当前温度和位移边界条件，在流体区域同时求 解预设的各个守恒方程的耦合解格式，在计算一个时间步Δt之后，得到当 前的热流和压力；

D、根据当前的热流和压力，在流体-固体耦合界面进行数据交换，得到 固体区域的边界条件；

E、根据所述固体区域的边界条件，在固体区域通过热力全耦合的方法 进行求解，计算一个时间步Δt之后，得到壁面温度T_w和结构位移u_s；

F、判读是否满足预设的停止条件，如果是，则停止整个流程；否则，返 回执行步骤B。

较佳的，所述设置当前边界条件包括：

由结构初始条件确定壁面温度和位移边界条件，进行高超声速稳态流动 计算，同时求解预设的各个守恒方程的耦合解格式，得到初始热流和初始压 力；

将所述初始热流和初始压力作为当前边界条件。

较佳的，所述数据交换包括：关联和插值。

较佳的，在进行数据交换时：

对于位移和温度采用最近邻居搜索方法计算；

对于压力和热流通量等载荷则采用守恒插值方法。

较佳的，所述停止条件为：当前的计算周期为最后一个计算周期。

较佳的，所述各个守恒方程包括：

连续守恒方程、动量守恒方程、能量等守恒方程和组分守恒方程。

如上可见，通过使用本发明的高超声速流动-传热与结构响应的多场耦合瞬 态数值的方法，可以解决传统方法中热/力载荷环境的预测与热防护结构性能分 析互相分离的问题，从而实现高超声速非平衡流动求解器与结构热/力全耦合求 解器相耦合的多场耦合计算，使得对高超声速飞行器的气动热力环境和结构热 力响应的预测更符合物理实际；而且，由于本发明中考虑了真实气体效应、耦 合传热和结构热力全耦合等复杂的效应，从而还可以大大提高高超声速飞行器 热/力载荷环境和结构响应预测的计算精度。

附图说明

图1为本发明实施例中的高超声速流动-传热与结构响应的多场耦合瞬 态数值的方法的流程示意图。

图2为本发明实施例中的多物理场耦合模型的示意图。

图3为本发明实施例中的耦合策略示意图。

图4为本发明实施例中的圆柱绕流计算模型。

图5为本发明实施例中计算得到的驻点温度随时间变化图。

图6为本发明实施例中计算得到的驻点热流随时间变化图。

图7为本发明实施例中计算得到的表面温度分布。

图8为本发明实施例中计算得到的结构Mises应力和位移图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下参照附图并举 实施例，对本发明进一步详细说明。

本实施例提供了一种高超声速流动-传热与结构响应的多场耦合瞬态数 值的方法。

图1为本发明实施例中的高超声速流动-传热与结构响应的多场耦合瞬 态数值的方法的流程示意图。如图1所示，本发明实施例中的高超声速流动- 传热与结构响应的多场耦合瞬态数值的方法可以包括如下所述的步骤：

步骤11，预先建立多物理场耦合模型并设置当前边界条件。

在本发明的技术方案中，首先需要预先建立相应的多物理场耦合模型。

在本发明的技术方案中，可以根据实际应用的需要，建立所需的多物理 场耦合模型。一般来说，可以使用多种方法来建立所述多物理场耦合。以下 将以其中的一种具体实施方式为例，对本发明的技术方案进行详细的介绍。

例如，图2为本发明实施例中的多物理场耦合模型的示意图。如图2所 示，本发明中的多物理场耦合模型具体可以包括：用于流体区域的流体分析 模型、用于固体区域的热-结构分析模型和用于数据交换的数据交换模型。

在所述用于流体区域的流体分析模型中，主要是利用同时求解连续、动 量、能量等守恒方程的耦合解格式进行流体区域中的气动热与气动力分析， 得到壁面热流q_w和壁面压力P_w。

在用于固体区域的热-结构分析模型中，主要是利用热力全耦合方法来进 行热-结构分析，得到壁面温度T_w和结构位移u_s。

在用于数据交换的数据交换模型中，主要是利用关联和插值等方法，将 流体分析模型的结果和热-结构分析模型的结果在网格间进行参数交换。

图3为本发明实施例中的耦合策略示意图。如图3所示，在本发明的技 术方案中，可以采用分区求解方法完成对高超声速流动-传热和结构响应的耦 合分析。其中，在流体区域和固体区域的求解器均为瞬态求解，每个求解器 所需要的数据在耦合界面上将进行反复交换。以下，将对上述耦合策略进行 详细的介绍。

另外，在本发明的技术方案中，在进行上述耦合策略之前，还需要预先 设置当前边界条件。

在本发明的技术方案中，可以根据实际应用的需要，使用多种方法来建 预先设置当前边界条件。以下将以其中的一种具体实施方式为例，对本发明 的技术方案进行详细的介绍。

较佳的，在本发明的具体实施例中，所述设置当前边界条件可以包括：

步骤11a，由结构初始条件确定壁面温度和位移边界条件，进行高超声 速稳态流动计算，同时求解预设的各个守恒方程的耦合解格式，得到初始热 流和初始压力。

在本步骤中，可以求解瞬态耦合分析的初始条件。即由结构初始条件确 定壁面温度和位移边界条件，进行高超声速稳态流动计算，并将计算果作为 瞬态耦合分析的初始条件，即初始热流和初始压力。

在本发明的技术方案中，可以使用多种具体实施方式来实现上述的步骤 11a，以下将以其中的一种具体实施方式为例，对本发明的技术方案进行详细 的介绍。

例如，在本发明的较佳实施例中，所述各个守恒方程包括：连续守恒方 程、动量守恒方程和能量等守恒方程。

考虑到高超声速飞行器周围的流场为化学非平衡、热力学平衡的粘性可 压缩连续流动，所述连续守恒方程、动量守恒方程和能量等守恒方程分别可 以表示为：

$\frac{\partial \rho }{\partial t}+▿·\left(\rho \overrightarrow{v}\right)=0---\left(1\right)$

$\frac{\partial \left(\rho \overrightarrow{v}\right)}{\partial t}+▿·\left(\rho \overrightarrow{v}\overrightarrow{v}\right)+▿p=▿·\overline{\stackrel{‾}{\tau }}---\left(2\right)$

$\frac{\partial \left(\rho E\right)}{\partial t}+▿·\left[\overrightarrow{v}\right(\rho E+p\left)\right]=▿·[\lambda ▿T-\underset{i}{\Sigma }{h}_i{\overrightarrow{J}}_i+(\overline{\stackrel{‾}{\tau }}·\overrightarrow{v}\left)\right]---\left(3\right)$

其中，ρ为密度，v为速度，p为压力，E为比总能量，J_i为扩散通量，为应力张量。

较佳的，在本发明的具体实施例中，对于包含组分混合或反应的流动， 所述各个守恒方程中还进一步还需要包括：组分守恒方程。该组分守恒方程 可以表示为：

$\frac{\partial \left({\mathrm{\rho Y}}_i\right)}{\partial t}+▿·\left(\rho \overrightarrow{v}{Y}_i\right)=-▿·{\overrightarrow{J}}_i+R_i---\left(4\right)$

其中，Y_i为组分质量分数，R_i为组分i的产物。

另外，在本发明的较佳实施例中，组分的分压可以由组分密度和混合气 体的温度求得，而混合气体的压力则可由Dalton定律给出：

p_i＝ρ_iR_iT(5)

$p=\underset{i}{\Sigma }{p}_i---\left(6\right)$

其中，p_i为组分i的压力。p为混合气体的总压力。

另外，在本发明的较佳实施例中，在高超声速化学非平衡流动计算中， 对于上述的连续守恒方程、动量守恒方程、能量守恒方程和组分守恒方程采 用有限体积法同时求解；而空间离散格式可采用AUSM+格式，时间上采用 交替方向隐式求解法。

此外，在本发明的技术方案中，由于存在动能耗散和激波，高超声速飞 行器周围的空气会达到极高的温度，而高温将使得气体离解甚至电离。因此， 在本发明的技术方案中，化学非平衡假设即特征化学反应时间与流动的特征 时间相当。可以采用Park5组分(O,N,NO,O2,N2)17反应化学动力学模型， 并考虑第三体效应。反应机制包括三个分解反应和两个交换反应。

其中，第r个反应的一般形式可以表示为：

其中，N表示化学组分数量，表示反应物化学计量系数，表示生成 物化学计量系数，M_i表示组分i。

另外，正反应速率常数k_f,r可以用Arrhenius形式表示，而逆反应速率常 数k_b,r则可通过正反应速率常数求出：

$k_{f,r}=A_r{T}^{{\beta }_r}\exp [-E_r/\mathrm{RT}]---\left(8\right)$

$k_{b,r}=\frac{k_{f,r}}{K_r}---\left(9\right)$

其中，A_r表示指前因子，β_r表示温度指数，E_r表示活化能，T表示温度， R表示通用气体常数，K_r表示平衡常数。

单一组分的粘性系数、热导率和扩散系数可由气体分子动力论给出：

${\mu }_i=\frac{2.67\times {10}^{-6}\sqrt{M_{i}T}}{{\sigma }_i^2{\Omega }_{\mu i}}---\left(10\right)$

${\lambda }_i=\frac{15}{4}\frac{R}{M_i}\mu [\frac{4}{15}\frac{C_{pi}{M}_i}{R}+\frac{1}{3}]---\left(11\right)$

$D_{ij}=0.00188\frac{{\left[T^3\left(\frac{1}{M_i}+\frac{1}{M_j}\right)\right]}^{1/2}}{{\mathrm{p\sigma }}_{ij}^2{\Omega }_{D_{ij}}}---\left(12\right)$

其中，μ_i为粘性系数，λ_i为热导率，D_ij为扩散系数。

而混合气体的系数可以由半经验Wilke公式给出，即：

$\mu =\underset{i}{\Sigma }\frac{X_i{\mu }_i}{{\Sigma }_j{X}_j{\phi }_{ij}}---\left(13\right)$

$\lambda =\underset{i}{\Sigma }\frac{X_i{\lambda }_i}{{\Sigma }_j{X}_j{\phi }_{ij}}---\left(14\right)$

式中，X_i表示组分i的摩尔分数，M_i表示组分i 的分子质量。

混合气体的多组分扩散系数则可由以下近似表达式计算：

$D_{im}=\frac{1-X_i}{\underset{j}{\Sigma }(X_i/D_{ij})}---\left(15\right)$

有限速率化学反应流动的计算需要每种单一组分的热力学属性(比热C_p,i和焓h_i)。本文中只考虑流动为化学非平衡、热力学平衡状态的反应气体， 每种组分的热力学属性由当地静温计算。高温气体的比热C_p,i和焓h_i可以通过 温度的多项式曲线拟合函数计算：

C_p,i＝R(A₁+A₂T+A₃T²+A₄T³+A₅T⁴)(16)

$h_i=RT(A_1+\frac{A_2}{2}T+\frac{A_3}{3}{T}^2+\frac{A_4}{4}{T}^3+\frac{A_5}{5}{T}^4+\frac{A_6}{T})---\left(17\right)$

则混合气体的属性可以由单一组分的属性确定：

$C_p=\underset{i=1}{\overset{NS}{\Sigma }}{C}_i{C}_{p,i},h=\underset{i=1}{\overset{NS}{\Sigma }}{C}_i{h}_i---\left(18\right)$

另外，根据傅里叶定律，热流可以表示为

$q=-\lambda \frac{\partial T}{\partial x}---\left(19\right)$

其中，λ是热导率。

因此，根据上述的公式(19)，通过对上述各个公式的求解，即可得到 初始热流；根据上述的公式(6)，并通过对上述各个公式的求解，即可得到 初始压力。

步骤11b，将所述初始热流和初始压力作为当前边界条件。

通过上述的步骤11a～11b，即可通过计算得到当前边界条件。

步骤12，根据结构确定壁面温度和位移边界条件，在流体-固体耦合界 面进行数据交换，得到流体区域的当前温度和位移边界条件。

在本步骤中，在获得当前边界条件之后，即可根据结构确定壁面温度和 位移边界条件，并通过在流体-固体耦合界面进行数据交换，从而可得到流体 区域的当前温度和位移边界条件，即经过数据交换之后的当前温度T_w和位移 u_s。其中，T_w表示在边界处的温度，而u_s则表示在边界处的位移。

在本发明的技术方案中，可以根据实际应用的需要，使用多种方法来进 行数据交换。以下将以其中的一种具体实施方式为例，对本发明的技术方案 进行详细的介绍。

在本发明的技术方案中，设置了两个区域：流体区域和固体区域。在上 述两个区域中分别进行相关计算时，需要将这两个区域分别进行网格划分， 各个数据均被定义在某种类型的网格上。然而，由于这两个区域的网格划分 是分别进行的，因此在两个区域的交界处，虽然都是使用网格描述同一个几 何体，但两个区域中的单元尺寸和网格节点的位置一般都不相同，这被称之 为“不匹配网格”。

因此，在将流体区域的数值作为固体区域的边界条件进行计算时，需要 先将流体区域中的边界条件(即网格中的数值)进行数据交换，转换成固体 区域中的网格中的数值，然后才能在固体区域中进行相应的计算。同理，在 流体区域中进行相应计算之前，也需要先将固体区域中的边界条件通过数据 交换转换成流体区域中的网格中的数值。

较佳的，在本发明的具体实施例中，在进行数据交换时，对于位移和温 度等场变量可采用最近邻居搜索方法计算；而对于压力和热流通量等载荷则 需要满足守恒性，即需要采用守恒插值方法。

较佳的，在本发明的具体实施例中，上述的所述数据交换可以包括：关 联和插值。对于每一个目标网格中的点，可以在源网格内搜索离其最近的单 元，产生一个节点-单元关系；而邻居搜索则可基于Kd-tree进行。

另外，在本发明的较佳实施例中，可以采用形函数算法完成上述的关联 和插值。如果使用的是线性单元，那么可以精确的映射线性函数，而精确的 映射二次函数则需要使用二次单元。

步骤13，根据流体区域的当前温度和位移边界条件，在流体区域同时求 解预设的各个守恒方程的耦合解格式，在计算一个时间步Δt之后，得到当 前的热流和压力。

由于在步骤12中可以得到流体区域的当前温度和位移边界条件，因此在 本步骤中即可根据流体区域的当前温度和位移边界条件进行高超声速稳定流 动计算，通过同时求解预设的各个守恒方程的耦合解格式的方式，在计算一 个时间步Δt之后，得到当前的热流和压力。

较佳的，在本发明的较佳具体实施例中，本步骤13中的“同时求解预设 的各个守恒方程的耦合解格式”的具体实现方式与上述步骤11a中的具体实 现方式可以是相同或相类似的，因此在此不再赘述。

步骤14，根据当前的热流和压力，在流体-固体耦合界面进行数据交换， 得到固体区域的边界条件。

在本步骤中，在获得当前的热流和压力之后，即可通过在流体-固体耦合 界面对当前边界条件(即当前的热流和压力)进行数据交换，从而可得到固 体区域(即结构求解器)的边界条件，即经过数据交换之后的热流q_w和压力 P_w。其中，q_w表示在边界处的热流，而P_w则表示在边界处的压力。

较佳的，在本发明的较佳具体实施例中，本步骤14中的数据交换与步骤 12中的数据交换可以使用相同的具体实现方式，因此在此不再赘述。

步骤15，根据所述固体区域的边界条件，在固体区域通过热力全耦合的 方法进行求解，计算一个时间步Δt之后，得到壁面温度T_w和结构位移u_s。

由于在上述步骤14中获得了固体区域的边界条件(即经过数据交换后的 热流q_w和压力P_w)，因此在本步骤中，即可根据所述固体区域的边界条件， 在固体区域(即热流、压力联合载荷下的结构场)中采用热力全耦合的方法 进行求解，在计算一个时间步Δt后，即可得到壁面温度T_w和结构位移u_s。

此外，在本发明的技术方案中，可以根据实际应用的需要，使用多种方 法来求解得到壁面温度T_w和结构位移u_s。以下将以其中的一种具体实施方式 为例，对本发明的技术方案进行详细的介绍。

例如，在本发明的技术方案中，基于能量守恒定律和Fourier定律，可以 得到结构瞬态热传导方程：

${\mathrm{\rho C}}_p\frac{\partial T}{\partial t}={\lambda }_s\frac{\partial }{\partial x_i}\left(\frac{\partial T}{\partial x_i}\right)---\left(20\right)$

当固体结构受到加热时，其体内温度将发生变化，此时在固体结构内部 将会有温差存在。由于热膨胀以及结构的约束作用，固体结构将产生热应力， 而热应力又导致了热变形的发生。对于二维结构的响应方程为：

$\frac{\partial }{\partial x}\left(\begin{array}{c}-{\sigma }_x\\ -{\tau }_{xy}\end{array}\right)+\frac{\partial }{\partial y}\left(\begin{array}{c}-{\tau }_{xy}\\ -{\sigma }_y\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)---\left(21\right)$

[σ]＝[D][B]δ^e＝[S]δ^e(22)

σ_x、σ_y和τ_xy为固体结构的应力分量；[D]是弹性矩阵，[B]和[S]分别为应 变、应力矩阵；δ^e为位移矩阵，其中，位移矩阵中的元素为结构位移u_s在x、 y、z方向上的分量。因此，根据该位移矩阵δ^e即可获知相应的结构位移u_s； 然后，根据上述的公式(20)～(22)，即可求解得到壁面温度T_w。

在固体区域内，可以采用有限元法进行热-结构分析。考虑应力与温度分 布之间的双向耦合关系，进行热力全耦合分析。温度用向后差分格式积分， 非线性耦合系统采用牛顿法求解。因此，在本发明的技术方案中，热力与结 构问题是同时求解的。

步骤16，判读是否满足预设的停止条件，如果是，则停止整个流程；否 则，返回执行步骤12。

在本发明的技术方案中，所述停止条件可以根据实际应用的需要而预先 设置。例如，较佳的，在本发明的具体实施例中，所述停止条件为：

当前的计算周期为最后一个计算周期。

在本发明的技术方案中，可以根据实际应用需要预先设置N个计算周期， 并在每一个计算周期中都执行上述的步骤12～15。因此，在步骤16中，即可 判断当前的计算周期是否为最后一个计算周期，如果是，则停止整个流程； 如果不是最后一个计算周期(即不是第N个计算周期)，则可以返回执行步 骤12，执行下一个计算周期。

另外，在本发明的技术方案中，所述时间步Δt的大小可以根据实际应 用需要而预先设置。例如，在本发明的较佳实施例中，所述时间步Δt的值 可以是：0.01秒(s)、0.001s或0.0001s等。

为验证本发明中所提出的高超声速流动-传热与结构响应的多场耦合瞬 态数值的方法的正确性，以下将以经典的圆柱绕流实验为例，对本发明中的 高超声速流动-传热与结构响应的多场耦合瞬态数值的方法进行实验。

在本发明中的上述实验中，所使用的来流马赫数、压力和温度分别为 6.47、648.1Pa和241.5K(如图4所示)；实验中所使用的不锈钢管的长度、 直径和厚度分别为0.6096m、0.0762m和0.0127m。流体区域和固体区域的网 格均为结构网格，但在耦合界面处不匹配。对边界层网格进行加密以保证计 算结果的网格无关性，并使其具有足够的分辨率。

图4为本发明实施例中的圆柱绕流计算模型。本发明中可以使用如图4 所示的流场和圆柱的计算模型。在本发明的技术方案中，计算过程中可以采 用SSTk-ω湍流模型进行计算。

图5为本发明实施例中计算得到的驻点温度随时间变化图。图5中提供 了不同耦合时间步长时驻点处的温度随时间变化曲线。

图6为本发明实施例中计算得到的驻点热流随时间变化图。图6中所示 为驻点热流随时间变化的趋势，其中图6的纵坐标q_stag即为边界处的热流q_w。 如图5和图6所示，随着加热时间的增长，驻点温度也将逐渐升高，而热流 则趋势相反。另外，温度和热流都在初始的一小段时间里变化剧烈，而后变 化趋于缓和。

在本发明的技术方案中，基于不同耦合时间步长得到了不同的温度分布 结果。其中，Δt＝0.0001s时得到t＝2s时的结果(384.9K)与实验结果(即图 5中所示的实验点处的数值，约为388.72K)最为接近，而对于Δt＝0.001s和 Δt＝0.01s，得到了略偏高的温度结果。因此可以得出结论，随着时间步长偏 离某个合适的特定值逐渐增大，预测的结构温度相应会逐渐偏高。

图7为本发明实施例中计算得到的表面温度分布。图7给出了Δt＝0.0001s 时不同时间点的温度分布。其中，图7中的横坐标为曲线长度(CurveLength)。 如图7所示，在气动加热开始后的前0.5s内，温度变化比较明显。

图8为本发明实施例中计算得到的结构Mises应力和位移图。如图8所 示，在热力联合载荷下，驻点处的位移(displacement)和Mises应力都随时 间逐渐增大，但是在当前工况下位移量较小，因此不足以明显影响流场。

通过将上述的实验结果与风洞试验的结果进行对比，验证了本发明中的 高超声速流动-传热与结构响应的多场耦合瞬态数值的方法的正确性和有效 性。

综上可知，通过使用本发明的高超声速流动-传热与结构响应的多场耦合 瞬态数值的方法，可以解决传统方法中热/力载荷环境的预测与热防护结构性 能分析互相分离的问题，从而实现高超声速非平衡流动求解器与结构热/力全 耦合求解器相耦合的多场耦合计算，使得对高超声速飞行器的气动热力环境 和结构热力响应的预测更符合物理实际；而且，由于本发明中考虑了真实气 体效应、耦合传热和结构热力全耦合等复杂的效应，从而还可以大大提高高 超声速飞行器热/力载荷环境和结构响应预测的计算精度。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本 发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在 本发明保护的范围之内。