[发明专利]一种稀薄连续统一的气体流动特性数值模拟方法

权利要求书

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。