CN107729691B - 一种稀薄连续统一的气体流动特性数值模拟方法 - Google Patents

Abstract

本发明涉及一种稀薄连续统一的气体流动特性数值模拟方法，在解决非线性耦合本构模型的守恒方程时，引入了直接关联守恒变量的空间导数S，其次将守恒变量U和附加变量S表达成一组基函数的线性表示，以对其进行离散化。将原始非结构网格下的守恒量转换到标准网格下，在控制体上积分，其每一步迭代可以得到本构关系方程中的黏性应力和热传导项初值，将其带入转化为微分方程初值问题的非线性代数本构方程组迭代得到黏性应力和热传导项，再返回代入守恒方程，根据具体问题满足误差条件后判断计算结束，输出流场流动特性。

Description

一种稀薄连续统一的气体流动特性数值模拟方法

技术领域

本发明属于计算流体力学的数值模拟算法，涉及一种稀薄连续统一的气体流动特性数值模拟方法，是基于计算流体力学领域中稀薄连续统一气体动力学非线性耦合本构模型的气体流动特性数值模拟算法。

背景技术

临近空间和微尺度流动的流动特性计算中均涉及连续稀薄耦合流动，且无法严格区分连续区域和稀薄区域。此时，气体的流动性质发生变化，不再适用流体力学中的连续性假设。因此发展连续稀薄耦合气体动力学的流动特性数值算法非常必要，以解决临近空间中飞行器气体动力学及微机电系统中微尺度流动特性的预测技术难题。

连续、动量、能量三大守恒方程耦合牛顿黏性定律和傅里叶热传导定律，即纳维-斯托克斯-傅里叶方程(NSF方程)，在过去近两百多年里不断推动着流体力学的发展。但是在微尺度流动和稀薄气体流动环境下，NSF方程的连续性假设失效。为了研究稀薄领域的流动现象，G.A.Bird.在“The DSMC Method”一书中研究了直接蒙特卡罗模拟方法，在稀薄领域取得了成功，直接推动了稀薄气体动力学的发展，但同时其研究表明DSMC在连续流和过渡流区域计算代价极大，模拟近连续气体流动所消耗的时间资源太多。即使强制性耦合单独解决连续流动的NSF方程和单独解决稀薄流动的DSMC方法，也很难解决连续稀薄耦合流动问题。

Eu在“Kinetic Theory and Irreversible Thermodynamic”一书中对不可逆热力学过程进行了研究，对Boltzmann方程进行了模型化处理，从Boltzmann方程的定解条件H定理出发，依据从非平衡态到平衡态熵增的特点，构造的分布函数只是一种分布函数的形态定义而不是严格的具体表达式。根据熵增耗散的物理概念构建了非平衡态到平衡态的熵增耗散数学模型，其在趋近平衡态时熵增数学模型收敛为Rayleigh-Onsager耗散函数，以此为基础建立了非平衡态到平衡态统一的非线性熵增模型，进而处理了Boltzmann方程的碰撞项。由上述原则，Eu从Boltzmann方程导出了黏性应力、热传导等高阶量的本构方程，耦合同样由Boltzmann方程导出的守恒方程，完成了气体动力学方程的封闭。随后Myong在一篇名为“A generalized hydrodynamic computational model for rarefied andmicroscale diatomic gas flow”的文献中处理了本构方程高阶项并以此为基础构建了非线性耦合本构关系。由于非线性耦合本构关系方程的强非线性关系，传统的有限体积法求解流动守恒方程和非线性耦合本构关系模型受到限制。因此，发展稀薄连续统一气体动力学非线性耦合本构模型的气体流动特性数值模拟算法是非常有必要和有意义的。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种稀薄连续统一的气体流动特性数值模拟方法，避免现有的有限体积法求解流动守恒方程和非线性耦合本构关系模型受到限制的问题。主要目的是构造一种求解非线性耦合本构模型的稀薄连续统一的气体流动特性数值模拟方法。非线性耦合本构模型包括守恒方程和本构关系方程，前者是时间空间上的微分方程，而后者是非线性代数方程，基于其不同的性质分别用不同的数值方法来求解，进而建立稀薄连续统一的气体流动特性数值模拟方法。

技术方案

一种稀薄连续统一的气体流动特性数值模拟方法，其特征在于步骤如下：

步骤1：对于流动守恒方程引入附加变量S，S为黏性应力和热传导等流动守恒参数高阶量的空间导数，定义附加变量S的系列方程组；

流动守恒方程：

其中，t表示时间，U表示守恒变量，F_inv(U)表示非黏性项，F_vis(U)表示黏性项，

表示求偏导，

表示速度梯度。ρ表示密度，u表示流体速度，E表示能量，p表示压强，γ表示比热比，I表示单位张量，Π表示黏性应力，Q表示热传导，Ec＝(γ-1)Ma²，Ma表示马赫数。Pr＝Cpη/λ表示普朗特数，η表示剪切黏度，Cp表示定压比热容，λ表示导热率；

附加变量S的系列方程组：

其中，X表示位置向量；

步骤2、对附加变量S的系列方程组进行离散化：

首先将守恒变量U和附加变量S表达成一组基函数的线性表示：

其中，A_i,B_i表示基函数系数，n的取值与精度相关，

表示非结构网格下的基函数，标准网格下的

为6个：

步骤3：将三角形非结构化网格的原始网格转化为标准三角形网格：

(r,s)-(x,y)：

(x,y)-(r,s)：

其中：D是原始三角形非结构网格的面积；

步骤4：继续将标准三角形网格转化为标准正方形网格

(a,b)-(r,s):

(r,s)-(a,b):

步骤5：将U_h，S_h代入附加变量S的系列方程组，方程两边同时乘以

对网格积分得：

其中，υ表示控制体，Γ表示控制体边界，

表示通量项，

表示积分项；

步骤6：基于基函数

有性质：

步骤5的方程组变为：

步骤7：在每个网格内，方程组每一步迭代可得到附加参数S的近似表达式S_h中的系数B_i，根据位置坐标就能求解任一点S_h值；

步骤8：在每一步迭代中根据

和固体边界高斯积分点坐标求解出固体边界高斯积分点的S_h值，继而由

求出对应的Π₀,Q₀,Δ₀，Δ代表附加体积正应力；

步骤9：得到Π₀,Q₀,Δ₀代入黏性力和热传导非线性本构方程

其中：q()为函数q()＝sinh()/()，f_b为附加应力相对黏性系数，c为分子模型参数，R为无量纲Rayleigh-Onsager耗散函数，

()^T为转置，tra()为迹，双点积︰计算为

步骤10：令Π,Q,Δ对应x，Π₀,Q₀,Δ₀为初值，本构方程对应F(x)，

设H(x,s)＝F(x)+(s-1)F(x₀)＝0,s∈[0,1],x∈D，D为n维向量空间Rⁿ上的区域；当s＝0时，H(x,0)＝F(x)-F(x₀)＝0的解是初值x₀。当s＝1时，H(x,0)＝F(x)＝0的解就是非线性本构方程组的解；

步骤11：做变换s＝1-e^-t，得H(x,t)＝F(x)-e^-tF(x₀)＝0,t∈[0,+∞),x∈D，t为时间；

步骤12：对步骤11的方程两边取时间导数，得微分方程初值问题

当t＝0,H(x,0)＝0的解是初值x₀，当t→+∞,H(x,+∞)＝0的解就是F(x)＝0的解；

步骤13：若计算中出现奇异点，则采用H₁(x,t)＝H(x,t)+η(t)G(x)计算，G(x)连续可微，且η(+∞)＝0,G(x₀)＝0，则H₁(x₀,0)＝H(x₀,0)＝0,H₁(x^*,+∞)＝H(x^*,+∞)＝0,

且对应H₁(x,t)＝0的微分方程初值问题为

步骤14：令η(t)＝Ae^-t,G(x)＝x-x₀，其中A是调节参数；采用四级四阶经典Runge-Kutta迭代求积公式求解得到Π,Q,Δ；

步骤15：完成迭代得到黏性应力和热传导Π,Q,Δ代入步骤1的守恒方程，迭代得守恒变量U的近似表达式U_h的系数A_i后，根据位置坐标求解任一点的流动特性参数ρ,u,T。

有益效果

本发明提出的一种稀薄连续统一的气体流动特性数值模拟方法，在解决非线性耦合本构模型的守恒方程时，引入了直接关联守恒变量的空间导数S，其次将守恒变量U和附加变量S表达成一组基函数的线性表示，以对其进行离散化。将原始非结构网格下的守恒量转换到标准网格下，在控制体上积分，其每一步迭代可以得到本构关系方程中的黏性应力和热传导项初值，将其带入转化为微分方程初值问题的非线性代数本构方程组迭代得到黏性应力和热传导项，再返回代入守恒方程，根据具体问题满足误差条件后判断计算结束，输出流场流动特性。

本发明的有益效果是：本发明和现有的解耦方法、牛顿迭代法解决本构关系方程相比，不影响本构关系本身的特性，而是基于整个本构关系同时考虑，算法作为一种大范围的收敛法，简单高效，其中的初值不必尽可能地接近根，只要满足一定的条件，就能得到非线性代数方程组的解。与解决守恒方程的有限体积法相比，算法能够保证气体动力学守恒方程的黏性应力、热传导等高阶量和守恒变量在同一精度，同时解决了高阶量边界条件赋值的难题，在边界上只需要定义守恒变量的边界条件，进而通过在每个网格内迭代出相应的方程系数，根据位置坐标就可求解任一点的守恒变量。

附图说明

图1为非线性耦合本构关系封闭气体动力学方程基本思路

图2为稀薄连续统一气体动力学方程计算流程图

图3为非结构三角形网格转化标准三角形网格原理示意图

图4为标准三角形网格转化标准正方形网格原理示意图

图5为黏性应力和热传导非线性本构方程计算流程图

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本实施例是一种连续稀薄统一的流动特性数值方法，具体分为以下步骤：

步骤1：结合图1，可以看到由Boltzmann方程，结合H定理及熵增推导出质量守恒方程、动量守恒方程、能量守恒方程，完成了气体动力学方程的封闭。

步骤2：结合图2，对于流动守恒方程：

其中，t表示时间，U表示守恒变量，F_inv(U)表示非黏性项，F_vis(U)表示黏性项，

表示求偏导，

表示速度梯度。ρ表示密度，u表示流体速度，E表示能量，p表示压强，γ表示比热比，I表示单位张量，Π表示黏性应力，Q表示热传导，Ec＝(γ-1)Ma²，Ma表示马赫数。Pr＝Cpη/λ表示普朗特数，η表示剪切黏度，Cp表示定压比热容，λ表示导热率。

引入附加变量S，S为黏性应力和热传导等流动守恒参数高阶量的空间导数，定义附加变量S的系列方程组：

其中，X表示位置向量。

步骤3：对方程组进行离散化，首先将守恒变量U和附加变量S表达成一组基函数的线性表示：

其中，A_i,B_i表示基函数系数，n的取值与精度相关，

表示非结构网格下的基函数，标准网格下的

为6个：

步骤4：为减少后续的计算量，使基函数满足正交性，需要进行网格转换。如图3所示，由于原始网格是三角形非结构化网格，需要先将其转化为标准三角形网格。

(r,s)-(x,y)：

(x,y)-(r,s)：

其中D是原始三角形非结构网格的面积。

步骤5：如图4所示，继续将标准三角形网格转化为标准正方形网格。

(a,b)-(r,s):

(r,s)-(a,b):

步骤6：将U_h，S_h代入方程组，方程两边同时乘以

对网格积分得：

其中，υ表示控制体，Γ表示控制体边界，

表示通量项，

表示积分项。

步骤7：由于基函数

有性质：

方程组变为：

步骤8：在每个网格内，方程组每一步迭代可得到附加参数S的近似表达式S_h中的系数B_i，根据位置坐标就能求解任一点S_h值，如

所示，等号右边第一项为通量，只需给出边界上的守恒变量值，第二项为体积分(二维为面积分)，根据网格内的高斯积分点求出。

步骤9：在每一步迭代中根据

和固体边界高斯积分点坐标求解出固体边界高斯积分点的S_h值，继而由

求出对应的Π₀,Q₀,Δ₀，Δ代表附加体积正应力。

步骤10：得到Π₀,Q₀,Δ₀代入黏性力和热传导非线性本构方程：

其中q()为函数q()＝sinh()/()，f_b为附加应力相对黏性系数，c为分子模型参数，R为无量纲Rayleigh-Onsager耗散函数，

()^T为转置，tra()为迹，双点积:计算为

步骤11：如图5，令Π,Q,Δ对应x，Π₀,Q₀,Δ₀为初值，本构方程对应F(x)，设H(x,s)＝F(x)+(s-1)F(x₀)＝0,s∈[0,1],x∈D，D为n维向量空间Rⁿ上的区域。当s＝0时，H(x,0)＝F(x)-F(x₀)＝0的解是初值x₀。当s＝1时，H(x,0)＝F(x)＝0的解就是非线性本构方程组的解。

步骤12：做变换s＝1-e^-t，得H(x,t)＝F(x)-e^-tF(x₀)＝0,t∈[0,+∞),x∈D，t为时间。

步骤13：对方程两边取时间导数，得微分方程初值问题，

当t＝0,H(x,0)＝0的解是初值x₀，当t→+∞,H(x,+∞)＝0的解就是F(x)＝0的解。

步骤14：若计算中出现奇异点，则采用

计算，G(x)连续可微，且η(+∞)＝0,G(x₀)＝0，则

且对应

的微分方程初值问题为

步骤15：令η(t)＝Ae^-t,G(x)＝x-x₀，其中A是调节参数。采用四级四阶经典Runge-Kutta迭代求积公式求解，

K₁＝Y(x_n)

步骤16：如图5完成迭代得到黏性应力和热传导Π,Q,Δ代入守恒方程，迭代得守恒变量U的近似表达式U_h的系数A_i后，可根据位置坐标求解任一点的U_h值。满足误差条件后，即得各坐标位置的流动特性参数ρ,u,T。

Claims (1)

1.一种稀薄连续统一的气体流动特性数值模拟方法，其特征在于步骤如下：

步骤1：对于流动守恒方程引入附加变量S，S为黏性应力和热传导等流动守恒参数高阶量的空间导数，定义附加变量S的系列方程组；

流动守恒方程：

其中，t表示时间，U表示守恒变量，F_inv(U)表示非黏性项，F_vis(U)表示黏性项，

表示求偏导，

表示速度梯度，ρ表示密度，u表示流体速度，E表示能量，p表示压强，γ表示比热比，I表示单位张量，Π表示黏性应力，Q表示热传导，Ec＝(γ-1)Ma²，Ma表示马赫数，Pr＝Cpη/λ表示普朗特数，η表示剪切黏度，Cp表示定压比热容，λ表示导热率；

附加变量S的系列方程组：

其中，X表示位置向量；

步骤2、对附加变量S的系列方程组进行离散化：

首先将守恒变量U和附加变量S表达成一组基函数的线性表示：

其中，A_i,B_i表示基函数系数，n的取值与精度相关，

表示非结构网格下的基函数，标准网格下的

为6个：

步骤3：将三角形非结构化网格的原始网格转化为标准三角形网格：

(r,s)-(x,y)：

(x,y)-(r,s)：

其中：D是原始三角形非结构网格的面积；

步骤4：继续将标准三角形网格转化为标准正方形网格

(a,b)-(r,s):

(r,s)-(a,b):

步骤5：将U_h，S_h代入附加变量S的系列方程组，方程两边同时乘以

对网格积分得：

其中，υ表示控制体，Γ表示控制体边界，

表示通量项，

表示积分项；

步骤6：基于基函数

有性质：

步骤5的方程组变为：

步骤7：在每个网格内，方程组每一步迭代可得到附加参数S的近似表达式S_h中的系数B_i，根据位置坐标就能求解任一点S_h值；

步骤8：在每一步迭代中根据

和固体边界高斯积分点坐标求解出固体边界高斯积分点的S_h值，继而由

求出对应的Π₀,Q₀,Δ₀，Δ代表附加体积正应力；

步骤9：得到Π₀,Q₀,Δ₀代入黏性力和热传导非线性本构方程

其中：q()为函数q()＝sinh()/()，f_b为附加应力相对黏性系数，c为分子模型参数，R为无量纲Rayleigh-Onsager耗散函数，

()^T为转置，tra()为迹，双点积︰计算为

步骤10：令Π,Q,Δ对应x，Π₀,Q₀,Δ₀为初值，本构方程对应F(x)，

设H(x,s)＝F(x)+(s-1)F(x₀)＝0,s∈[0,1],x∈D，D为n维向量空间Rⁿ上的区域；当s＝0时，H(x,0)＝F(x)-F(x₀)＝0的解是初值x₀，当s＝1时，H(x,0)＝F(x)＝0的解就是非线性本构方程组的解；

步骤11：做变换s＝1-e^-t，得H(x,t)＝F(x)-e^-tF(x₀)＝0,t∈[0,+∞),x∈D，t为时间；

步骤12：对步骤11的方程两边取时间导数，得微分方程初值问题

当t＝0,H(x,0)＝0的解是初值x₀，当t→+∞,H(x,+∞)＝0的解就是F(x)＝0的解；

步骤13：若计算中出现奇异点，则采用H₁(x,t)＝H(x,t)+η(t)G(x)计算，G(x)连续可微，且η(+∞)＝0,G(x₀)＝0，则H₁(x₀,0)＝H(x₀,0)＝0,H₁(x^*,+∞)＝H(x^*,+∞)＝0,

且对应H₁(x,t)＝0的微分方程初值问题为

步骤14：令η(t)＝Ae^-t,G(x)＝x-x₀，其中A是调节参数；采用四级四阶经典Runge-Kutta迭代求积公式求解得到Π,Q,Δ；

步骤15：完成迭代得到黏性应力和热传导Π,Q,Δ代入步骤1的守恒方程，迭代得守恒变量U的近似表达式U_h的系数A_i后，根据位置坐标求解任一点的流动特性参数ρ,u,T。

Images