[发明专利]基于FVM-TLBFS方法的飞行器热防护结构热传导计算方法

权利要求书

1.基于FVM-TLBFS方法的飞行器热防护结构热传导计算方法，其特征在于，针对热防护结构的传热特性问题，将LBE方法推广到结构热传导计算中；TLBFS方法可以应用于非均匀网格中；采用FVM方法对热防护结构进行空间离散，并利用热格子Boltzmann方程模型的局部解重构单元界面的数值通量，将D3Q6或D2Q4模型模型应用于TLBFS方法用于求解三维或二维结构热传导方程；采用基于格心格式的有限体积方法进行空间离散，采用显示Runge-Kutta方法进行时间推进。

2.根据权利要求1所述的基于FVM-TLBFS方法的飞行器热防护结构热传导计算方法，其特征在于，飞行器热防护结构二维结构热传导计算方法具体如下：

1)基于Boltzmann方程的结构热传导方程推导

对于各项同性材料，结构热传导控制方程可写成如下积分形式：

其中，ρ_s为固体结构密度，c_s为定压比热容，T为固体温度，k_s为固体传热系数，为温度梯度，n为单位法向矢量，S为控制体表面面积，Ω为控制体体积，Q为热流源项；

对于不可压缩流动传热问题，传热方程还可以写成如下形似

其中，u为结构传热方程，κ为热扩散系数，为梯度算子；

当热流源项为零时Q＝0，对于结构传热方程u＝0，方程(1)与方程(2)是等价的；方程(1)中的热传导系数k_s与方程(2)中的热扩散系数κ之间的关系如下，

为了还原宏观传热方程(2)，我们通过构造如下Boltzmann方程来实现：

式中，h为温度分布函数，h^eq是温度分布函数h在碰撞时间尺度τ_κ内通过粒子碰撞接近的平衡态，ξ＝(ξ₁,ξ₂)是粒子速度空间中的粒子速度，▽h为为温度分布函数梯度；

2)热格子Boltzmann模型D2Q4模型

下面将对D2Q4模型进行详细介绍：

由于采用离散的格子LBE模型，Boltzmann方程(4)可以重新写成如下形式：

在方程(4)中，粒子速度ξ可以写成e_α；在二维传热方程计算中，D2Q4-LBE模型可应用于离散的格子速度空间；在D2Q4模型中，平衡分布函数和离散粒子速度可以通过如下式子计算：

通过Chapman-Enskog展开分析，可以分别得到能量方程的分布函数与温度、法向通量的关系如下：

式中，为单元界面法向守恒变量，为单元界面法向通量，u₁为一维方向宏观速度，e_α,1为当地坐标系下粒子速度的第一个分量；另外，方程(9)中热扩散系数κ与碰撞时间尺度τ_κ可以写成如下关系式：

τ_κ＝2κ (10)

3)能量方程通量计算

对于热传导方程，单元界面的温度分布函数也可以表示为

式中h_α(0,t)为x＝0位置处的分布函数，为x＝0位置处的平衡分布函数，为x＝0位置处的非平衡分布函数，为无量纲碰撞时间；

将式(11)代入方程(9)中，可以的得到单元界面的法向能量通量为，

由此可见，法向能量通量也由两个部分组成：其一是F^s,I，代表界面平衡分布函数产生的贡献；其二是F^s,II，代表单元界面周围平衡分布函数的贡献；F^s,I与宏观流动速度直接相关；当u＝0时，F^s,I也为零，因此对于固体结构传热u＝0，只需要对第二部分通量F^s,II进行计算；

计算F^s,II时，应先求出单元界面附近点的平衡分布函数对于热流场中任意变量它在单元界面附近点的平衡分布函数科可以通过如下式子计算，

式中，和分别表示单元界面左右两边值和的梯度值；和则分别表示左右单元值和的梯度的平均值，即，和

联立方程(12)和(13)，对于二维传热方程计算，采用D2Q4模型可以求得

进一步地，将式(14)代入方程(12)中，同时令u＝0，可得到完整二维固体传热能量方程通量计算表达式为

至此，上式中最后未确定的参数为迁移时间步长δt，基于避免数值外插的原则，δt由如下式子进行计算，

2D：

式中，Δl和Δr分别为界面左右单元最短的边长；

4)计算流程

上面已经详细描述了FVM-TLBFS方法的整个推导过程和计算方法，下面将给出该求解器实施的具体步骤，如下：

(1)计算各网格单元守恒变量的导数，通过线性重构得到单元界面两侧的守恒变量；

(2)将单元界面两侧的守恒变量代入方程(6)、(7)、(13)，计算出平衡分布函数；

(3)通过方程(8)计算出单元界面上的守恒变量；

(4)通过方程(9)、(14)、(15)计算能量通量；

(6)通过方程(16)计算迁移时间步长δt；

(7)通过时间推进方法求解控制方程(1)，这一步能够计算出下一时间步网格单元中心新的守恒变量值；

(8)重复步骤(1)--(7)，直到得到满足条件的温度分布收敛解。

3.根据权利要求1所述的基于FVM-TLBFS方法的飞行器热防护结构热传导计算方法，其特征在于，飞行器热防护结构三维结构热传导计算方法具体如下：

1)基于Boltzmann方程的结构热传导方程推导

对于各项同性材料，结构热传导控制方程可写成如下积分形式：

其中，ρ_s为固体结构密度，c_s为定压比热容，T为固体温度，k_s为固体传热系数，为温度梯度，n为单位法向矢量，S为控制体表面面积，Ω为控制体体积，Q为热流源项；

对于不可压缩流动传热问题，传热方程还可以写成如下形式

其中，u为结构传热方程，κ为热扩散系数，为梯度算子；

当热流源项为零时Q＝0，对于结构传热方程u＝0，方程(1)与方程(2)是等价的；方程(1)中的热传导系数k_s与方程(2)中的热扩散系数κ之间的关系如下，

为了还原宏观传热方程(2)，通过构造如下Boltzmann方程来实现：

式中，h为温度分布函数，h^eq是温度分布函数h在碰撞时间尺度τ_κ内通过粒子碰撞接近的平衡态，ξ＝(ξ₁,ξ₂)是粒子速度空间中的粒子速度，▽h为温度分布函数梯度；

2)热格子Boltzmann模型D3Q6模型

下面将对D3Q6模型进行详细介绍：

由于采用离散的格子LBE模型，Boltzmann方程(4)可以重新写成如下形式：

在方程(4)中，粒子速度ξ可以成e_α；在二维传热方程计算中，D2Q4-LBE模型可应用于离散的格子速度空间；在D3Q6模型中，平衡分布函数和离散粒子速度通过如下式子计算：

e_α＝(±1,0,0),(0,±1,0),(0,0,±1),α＝1,2,3,4,5,6

通过两元扩散数Chapman-Enskog展开分析，将分布函数，时间导数和空间导数展开为如下形式

式中，为平衡分布函数，ε是一个与Knudsen克努森数成比例的小参数，h_α为分布函数，为ε级展开分布函数，为ε²级展开分布函数，为时间导数，为t1级时间导数其中t1＝εt，为t2级时间导数其中t2＝εt²，与均为空间导数；将式(7)代入方程(5)，得到如下三个表达式：

O(ε⁰):

O(ε¹):

O(ε²):

其中，O(ε⁰)为ε⁰级截断误差，O(ε¹)为ε¹级截断误差，O(ε²)为ε²级截断误差；

通过D3Q6模型将方程(9)和方程(10)关于指标α求和，可以得到

根据Boussinesq近似，∏⁽¹⁾可以表示为

这里，Ma为马赫数；联立方程(11)和方程(12)，可以得到

通过对比，发现方程热扩散系数κ与碰撞时间尺度τ_κ的关系如下：

τ_κ＝3κ (14)

此外，恢复能量方程的误差项为对于不可压缩流动的模拟，这个误差项忽略不计；

由此，分别得到能量方程的分布函数与温度、法向通量的关系如下：

式中，为单元界面法向守恒变量，为单元界面法向通量，u₁为一维方向宏观速度，x₁表示一维方向，e_α,1为当地坐标系下粒子速度的第一个分量；

3)能量方程通量计算

对于热传导方程，单元界面的温度分布函数也可以表示为

式中，h_α(0,t)为x＝0位置处的分布函数，为x＝0位置处的平衡分布函数，为x＝0位置处的非平衡分布函数，为无量纲碰撞时间；

将式(17)代入方程(16)中，可以的得到单元界面的法向能量通量为，

由此可见，法向能量通量也由两个部分组成：其一是F^s,I，代表界面平衡分布函数产生的贡献；其二是F^s,II，代表单元界面周围平衡分布函数的贡献；由于F^s,I与宏观流动速度直接相关；当u＝0时，F^s,I也为零，因此对于固体结构传热u＝0，只需要对第二部分通量F^s,II进行计算；

计算F^s,II时，应先求出单元界面附近点的平衡分布函数对于热流场中任意变量它在单元界面附近点的平衡分布函数科可以通过如下式子计算，

式中，和分别表示单元界面左右两边值和的梯度值；和则分别表示左右单元值和的梯度的平均值，即，和

联立方程(18)和(19)，对于三维传热方程计算，采用D3Q6模型可以求得

将式(20)代入方程(19)中，同时令u＝0，可得到完整三维固体传热能量方程通量计算表达式为

至此，上式中最后未确定的参数为迁移时间步长δt，基于避免数值外插的原则，δt由如下式子进行计算，

3D：

式中，Δl和Δr分别为界面左右单元最短的边长；

4)计算流程

上面已经详细描述了FVM-TLBFS方法的整个推导过程和计算方法，下面将给出该求解器实施的具体步骤，如下：

(1)计算各网格单元守恒变量的导数，通过线性重构得到单元界面两侧的守恒变量；

(2)将单元界面两侧的守恒变量代入方程(6)-(19)，计算出平衡分布函数；

(3)通过方程(15)-(17)计算出单元界面上的守恒变量；

(4)通过方程(18)、(20)、(21)计算能量通量；

(6)通过方程(22)计算迁移时间步长δt；

(7)通过时间推进方法求解控制方程(1)，这一步能够计算出下一时间步网格单元中心新的守恒变量值；

(8)重复步骤(1)--(7)，直到得到满足条件的温度分布收敛解。