CN107729691A - 一种稀薄连续统一的气体流动特性数值模拟方法 - Google Patents
一种稀薄连续统一的气体流动特性数值模拟方法 Download PDFInfo
- Publication number
- CN107729691A CN107729691A CN201711110884.6A CN201711110884A CN107729691A CN 107729691 A CN107729691 A CN 107729691A CN 201711110884 A CN201711110884 A CN 201711110884A CN 107729691 A CN107729691 A CN 107729691A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mtd
- mfrac
- mtr
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
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表示守恒变量,Finv(U)表示非黏性项,Fvis(U)表示黏性项,表示求偏导,表示速度梯度。ρ表示密度,U表示流体速度,E表示能量,p表示压强,γ表示比热比,I表示单位张量,Π表示黏性应力,Q表示热传导,Ec=(γ-1)Ma2,Ma表示马赫数。Pr=Cpη/λ表示普朗特数,η表示剪切黏度,Cp表示定压比热容,λ表示导热率;
附加变量S的系列方程组:
其中,X表示位置向量;
步骤2、对附加变量S的系列方程组进行离散化:
首先将守恒变量U和附加变量S表达成一组基函数的线性表示:
其中,Ai,Bi表示基函数系数,n的取值与精度相关,表示非结构网格下的基函数,标准网格下的为6个:
步骤3:将三角形非结构化网格的原始网格转化为标准三角形网格:
(r,s)-(x,y):
(x,y)-(r,s):
其中:D是原始三角形非结构网格的面积;
步骤4:继续将标准三角形网格转化为标准正方形网格
(a,b)-(r,s):
(r,s)-(a,b):
步骤5:将Uh,Sh代入附加变量S的系列方程组,方程两边同时乘以对网格积分得:
其中,υ表示控制体,Γ表示控制体边界,表示通量项,表示积分项;
步骤6:基于基函数有性质:
步骤5的方程组变为:
步骤7:在每个网格内,方程组每一步迭代可得到附加参数S的近似表达式Sh中的系数Bi,根据位置坐标就能求解任一点Sh值;
步骤8:在每一步迭代中根据和固体边界高斯积分点坐标求解出固体边界高斯积分点的Sh值,继而由求出对应的Π0,Q0,Δ0,Δ代表附加体积正应力;
步骤9:得到Π0,Q0,Δ0代入黏性力和热传导非线性本构方程
其中:q()为函数q()=sinh()/(),fb为附加应力相对黏性系数,c为分子模型参数,R为无量纲Rayleigh-Onsager耗散函数,()T为转置,tra()为迹,双点积︰计算为
步骤10:令Π,Q,Δ对应x,Π0,Q0,Δ0为初值,本构方程对应F(x),
设H(x,s)=F(x)+(s-1)F(x0)=0,s∈[0,1],x∈D,D为n维向量空间Rn上的区域;当s=0时,H(x,0)=F(x)-F(x0)=0的解是初值x0。当s=1时,H(x,0)=F(x)=0的解就是非线性本构方程组的解;
步骤11:做变换s=1-e-t,得H(x,t)=F(x)-e-tF(x0)=0,t∈[0,+∞),x∈D,t为时间;
步骤12:对步骤11的方程两边取时间导数,得微分方程初值问题
当t=0,H(x,0)=0的解是初值x0,当t→+∞,H(x,+∞)=0的解就是F(x)=0的解;
步骤13:若计算中出现奇异点,则采用H1(x,t)=H(x,t)+η(t)G(x)计算,G(x)连续可微,且η(+∞)=0,G(x0)=0,则H1(x0,0)=H(x0,0)=0,H1(x*,+∞)=H(x*,+∞)=0,
且对应H1(x,t)=0的微分方程初值问题为
步骤14:令η(t)=Ae-t,G(x)=x-x0,其中A是调节参数;采用四级四阶经典Runge-Kutta迭代求积公式求解得到Π,Q,Δ;
步骤15:完成迭代得到黏性应力和热传导Π,Q,Δ代入步骤1的守恒方程,迭代得守恒变量U的近似表达式Uh的系数Ai后,根据位置坐标求解任一点的流动特性参数ρ,u,T。
有益效果
本发明提出的一种稀薄连续统一的气体流动特性数值模拟方法,在解决非线性耦合本构模型的守恒方程时,引入了直接关联守恒变量的空间导数S,其次将守恒变量U和附加变量S表达成一组基函数的线性表示,以对其进行离散化。将原始非结构网格下的守恒量转换到标准网格下,在控制体上积分,其每一步迭代可以得到本构关系方程中的黏性应力和热传导项初值,将其带入转化为微分方程初值问题的非线性代数本构方程组迭代得到黏性应力和热传导项,再返回代入守恒方程,根据具体问题满足误差条件后判断计算结束,输出流场流动特性。
本发明的有益效果是:本发明和现有的解耦方法、牛顿迭代法解决本构关系方程相比,不影响本构关系本身的特性,而是基于整个本构关系同时考虑,算法作为一种大范围的收敛法,简单高效,其中的初值不必尽可能地接近根,只要满足一定的条件,就能得到非线性代数方程组的解。与解决守恒方程的有限体积法相比,算法能够保证气体动力学守恒方程的黏性应力、热传导等高阶量和守恒变量在同一精度,同时解决了高阶量边界条件赋值的难题,在边界上只需要定义守恒变量的边界条件,进而通过在每个网格内迭代出相应的方程系数,根据位置坐标就可求解任一点的守恒变量。
附图说明
图1为非线性耦合本构关系封闭气体动力学方程基本思路
图2为稀薄连续统一气体动力学方程计算流程图
图3为非结构三角形网格转化标准三角形网格原理示意图
图4为标准三角形网格转化标准正方形网格原理示意图
图5为黏性应力和热传导非线性本构方程计算流程图
具体实施方式
现结合实施例、附图对本发明作进一步描述:
本实施例是一种连续稀薄统一的流动特性数值方法,具体分为以下步骤:
结合图1,可以看到由Boltzmann方程,结合H定理及熵增推导出质量守恒方程、动量守恒方程、能量守恒方程,完成了气体动力学方程的封闭。
步骤2:结合图2,对于流动守恒方程:
其中,t表示时间,U表示守恒变量,Finv(U)表示非黏性项,Fvis(U)表示黏性项,表示求偏导,表示速度梯度。ρ表示密度,U表示流体速度,E表示能量,p表示压强,γ表示比热比,I表示单位张量,Π表示黏性应力,Q表示热传导,Ec=(γ-1)Ma2,Ma表示马赫数。Pr=Cpη/λ表示普朗特数,η表示剪切黏度,Cp表示定压比热容,λ表示导热率。
引入附加变量S,S为黏性应力和热传导等流动守恒参数高阶量的空间导数,定义附加变量S的系列方程组:
其中,X表示位置向量。
步骤3:对方程组进行离散化,首先将守恒变量U和附加变量S表达成一组基函数的线性表示:
其中,Ai,Bi表示基函数系数,n的取值与精度相关,表示非结构网格下的基函数,标准网格下的为6个:
步骤4:为减少后续的计算量,使基函数满足正交性,需要进行网格转换。如图3所示,由于原始网格是三角形非结构化网格,需要先将其转化为标准三角形网格。
(r,s)-(x,y):
(x,y)-(r,s):
其中D是原始三角形非结构网格的面积。
步骤5:如图4所示,继续将标准三角形网格转化为标准正方形网格。
(a,b)-(r,s):
(r,s)-(a,b):
步骤6:将Uh,Sh代入方程组,方程两边同时乘以对网格积分得:
其中,υ表示控制体,Γ表示控制体边界,表示通量项,表示积分项。
步骤7:由于基函数有性质:
方程组变为:
步骤8:在每个网格内,方程组每一步迭代可得到附加参数S的近似表达式Sh中的系数Bi,根据位置坐标就能求解任一点Sh值,如所示,等号右边第一项为通量,只需给出边界上的守恒变量值,第二项为体积分(二维为面积分),根据网格内的高斯积分点求出。
步骤9:在每一步迭代中根据和固体边界高斯积分点坐标求解出固体边界高斯积分点的Sh值,继而由求出对应的Π0,Q0,Δ0,Δ代表附加体积正应力。
步骤10:得到Π0,Q0,Δ0代入黏性力和热传导非线性本构方程:
其中q()为函数q()=sinh()/(),fb为附加应力相对黏性系数,c为分子模型参数,R为无量纲Rayleigh-Onsager耗散函数,()T为转置,tra()为迹,双点积:计算为A:
步骤11:如图5,令Π,Q,Δ对应x,Π0,Q0,Δ0为初值,本构方程对应F(x),设H(x,s)=F(x)+(s-1)F(x0)=0,s∈[0,1],x∈D,D为n维向量空间Rn上的区域。当s=0时,H(x,0)=F(x)-F(x0)=0的解是初值x0。当s=1时,H(x,0)=F(x)=0的解就是非线性本构方程组的解。
步骤12:做变换s=1-e-t,得H(x,t)=F(x)-e-tF(x0)=0,t∈[0,+∞),x∈D,t为时间。
步骤13:对方程两边取时间导数,得微分方程初值问题,
当t=0,H(x,0)=0的解是初值x0,当t→+∞,H(x,+∞)=0的解就是F(x)=0的解。
步骤14:若计算中出现奇异点,则采用计算,G(x)连续可微,且η(+∞)=0,G(x0)=0,则且对应的微分方程初值问题为
步骤15:令η(t)=Ae-t,G(x)=x-x0,其中A是调节参数。采用四级四阶经典Runge-Kutta迭代求积公式求解,
K1=Y(xn)
步骤16:如图5完成迭代得到黏性应力和热传导Π,Q,Δ代入守恒方程,迭代得守恒变量U的近似表达式Uh的系数Ai后,可根据位置坐标求解任一点的Uh值。满足误差条件后,即得各坐标位置的流动特性参数ρ,u,T。
Claims (1)
1.一种稀薄连续统一的气体流动特性数值模拟方法,其特征在于步骤如下:
步骤1:对于流动守恒方程引入附加变量S,S为黏性应力和热传导等流动守恒参数高阶量的空间导数,定义附加变量S的系列方程组;
流动守恒方程:
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>U</mi>
<mo>=</mo>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<mi>&rho;</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&rho;</mi>
<mi>u</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&rho;</mi>
<mi>E</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>F</mi>
<mrow>
<mi>i</mi>
<mi>n</mi>
<mi>v</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>U</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>&rho;</mi>
<mi>u</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&rho;</mi>
<mi>u</mi>
<mi>u</mi>
<mo>+</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msup>
<mi>&gamma;Ma</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mi>p</mi>
<mi>I</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<mi>&rho;</mi>
<mi>E</mi>
<mo>+</mo>
<mfrac>
<mn>1</mn>
<mrow>
<msup>
<mi>&gamma;Ma</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mi>P</mi>
<mo>)</mo>
<mi>u</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>F</mi>
<mrow>
<mi>v</mi>
<mi>i</mi>
<mi>s</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>U</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>Re</mi>
</mfrac>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>&Pi;</mo>
<mo>+</mo>
<msub>
<mi>f</mi>
<mi>b</mi>
</msub>
<mi>&Delta;</mi>
<mi>I</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<mo>&Pi;</mo>
<mo>+</mo>
<msub>
<mi>f</mi>
<mi>b</mi>
</msub>
<mi>&Delta;</mi>
<mi>I</mi>
<mo>)</mo>
<mi>u</mi>
<mo>+</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mi>E</mi>
<mi>c</mi>
<mi>Pr</mi>
</mrow>
</mfrac>
<mi>Q</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中,t表示时间,U表示守恒变量,Finv(U)表示非黏性项,Fvis(U)表示黏性项,表示求偏导,▽表示速度梯度。ρ表示密度,U表示流体速度,E表示能量,p表示压强,γ表示比热比,I表示单位张量,Π表示黏性应力,Q表示热传导,Ec=(γ-1)Ma2,Ma表示马赫数。Pr=Cpη/λ表示普朗特数,η表示剪切黏度,Cp表示定压比热容,λ表示导热率;
附加变量S的系列方程组:
其中,X表示位置向量;
步骤2、对附加变量S的系列方程组进行离散化:
首先将守恒变量U和附加变量S表达成一组基函数的线性表示:
其中,Ai,Bi表示基函数系数,n的取值与精度相关,表示非结构网格下的基函数,标准网格下的为6个:
步骤3:将三角形非结构化网格的原始网格转化为标准三角形网格:
(r,s)-(x,y):
<mrow>
<mi>x</mi>
<mo>=</mo>
<mo>-</mo>
<mfrac>
<mrow>
<mi>r</mi>
<mo>+</mo>
<mi>s</mi>
</mrow>
<mn>2</mn>
</mfrac>
<msub>
<mi>x</mi>
<mn>1</mn>
</msub>
<mo>+</mo>
<mfrac>
<mrow>
<mi>r</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mn>2</mn>
</mfrac>
<msub>
<mi>x</mi>
<mn>2</mn>
</msub>
<mo>+</mo>
<mfrac>
<mrow>
<mi>s</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mn>2</mn>
</mfrac>
<msub>
<mi>x</mi>
<mn>3</mn>
</msub>
<mo>,</mo>
</mrow>
<mrow>
<mi>y</mi>
<mo>=</mo>
<mo>-</mo>
<mfrac>
<mrow>
<mi>r</mi>
<mo>+</mo>
<mi>s</mi>
</mrow>
<mn>2</mn>
</mfrac>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mo>+</mo>
<mfrac>
<mrow>
<mi>r</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mn>2</mn>
</mfrac>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mo>+</mo>
<mfrac>
<mrow>
<mi>s</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mn>2</mn>
</mfrac>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mo>.</mo>
</mrow>
(x,y)-(r,s):
<mrow>
<mi>r</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mn>2</mn>
<mi>x</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mo>-</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mn>2</mn>
<mi>y</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mn>3</mn>
</msub>
<mo>-</mo>
<msub>
<mi>x</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mn>3</mn>
</msub>
<mo>-</mo>
<msub>
<mi>x</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mo>+</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mo>-</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mn>3</mn>
</msub>
<mo>+</mo>
<msub>
<mi>x</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mn>2</mn>
<mi>D</mi>
</mrow>
</mfrac>
<mo>,</mo>
</mrow>
<mrow>
<mi>s</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mn>2</mn>
<mi>x</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mn>2</mn>
<mi>y</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>x</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>x</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mn>3</mn>
</msub>
<mo>+</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mo>-</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mn>3</mn>
</msub>
<mo>+</mo>
<msub>
<mi>x</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mo>-</mo>
<mn>2</mn>
<mi>D</mi>
</mrow>
</mfrac>
<mo>.</mo>
</mrow>
其中:D是原始三角形非结构网格的面积;
步骤4:继续将标准三角形网格转化为标准正方形网格
(a,b)-(r,s):
<mrow>
<mi>r</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<mi>a</mi>
<mo>)</mo>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mi>b</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</mfrac>
<mo>-</mo>
<mn>1</mn>
<mo>,</mo>
<mi>s</mi>
<mo>=</mo>
<mi>b</mi>
</mrow>
(r,s)-(a,b):
<mrow>
<mi>a</mi>
<mo>=</mo>
<mn>2</mn>
<mfrac>
<mrow>
<mn>1</mn>
<mo>+</mo>
<mi>r</mi>
</mrow>
<mrow>
<mn>1</mn>
<mo>-</mo>
<mi>s</mi>
</mrow>
</mfrac>
<mo>-</mo>
<mn>1</mn>
<mo>,</mo>
<mi>b</mi>
<mo>=</mo>
<mi>s</mi>
<mo>;</mo>
</mrow>
步骤5:将Uh,Sh代入附加变量S的系列方程组,方程两边同时乘以对网格积分得:
其中,υ表示控制体,Γ表示控制体边界,表示通量项,表示积分项;
步骤6:基于基函数有性质:
步骤5的方程组变为:
步骤7:在每个网格内,方程组每一步迭代可得到附加参数S的近似表达式Sh中的系数Bi,根据位置坐标就能求解任一点Sh值;
步骤8:在每一步迭代中根据和固体边界高斯积分点坐标求解出固体边界高斯积分点的Sh值,继而由求出对应的Π0,Q0,Δ0,Δ代表附加体积正应力;
步骤9:得到Π0,Q0,Δ0代入黏性力和热传导非线性本构方程
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>&Pi;</mo>
<mi>q</mi>
<mrow>
<mo>(</mo>
<mi>c</mi>
<mi>R</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<msub>
<mi>f</mi>
<mi>b</mi>
</msub>
<mi>&Delta;</mi>
<mo>)</mo>
</mrow>
<msub>
<mo>&Pi;</mo>
<mn>0</mn>
</msub>
<mo>+</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<mo>&Pi;</mo>
<mo>&CenterDot;</mo>
<mo>&dtri;</mo>
<mi>u</mi>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&Delta;</mi>
<mi>q</mi>
<mrow>
<mo>(</mo>
<mi>c</mi>
<mi>R</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>&Delta;</mi>
<mn>0</mn>
</msub>
<mo>+</mo>
<mfrac>
<mn>3</mn>
<mn>2</mn>
</mfrac>
<msub>
<mi>f</mi>
<mi>b</mi>
</msub>
<mo>(</mo>
<mo>&Pi;</mo>
<mo>+</mo>
<msub>
<mi>f</mi>
<mi>b</mi>
</msub>
<mi>&Delta;</mi>
<mi>I</mi>
<mo>)</mo>
</mrow>
<mo>:</mo>
<mo>&dtri;</mo>
<mi>u</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>Q</mi>
<mi>q</mi>
<mrow>
<mo>(</mo>
<mi>c</mi>
<mi>R</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>+</mo>
<msub>
<mi>f</mi>
<mi>b</mi>
</msub>
<mi>&Delta;</mi>
<mo>)</mo>
</mrow>
<msub>
<mi>Q</mi>
<mn>0</mn>
</msub>
<mo>+</mo>
<mo>&Pi;</mo>
<mo>&CenterDot;</mo>
<msub>
<mi>Q</mi>
<mn>0</mn>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中:q()为函数q()=sinh()/(),fb为附加应力相对黏性系数,c为分子模型参数,R为无量纲Rayleigh-Onsager耗散函数,()T为转置,tra()为迹,双点积︰计算为
步骤10:令Π,Q,Δ对应x,Π0,Q0,Δ0为初值,本构方程对应F(x),
设H(x,s)=F(x)+(s-1)F(x0)=0,s∈[0,1],x∈D,D为n维向量空间Rn上的区域;当s=0时,H(x,0)=F(x)-F(x0)=0的解是初值x0。当s=1时,H(x,0)=F(x)=0的解就是非线性本构方程组的解;
步骤11:做变换s=1-e-t,得H(x,t)=F(x)-e-tF(x0)=0,t∈[0,+∞),x∈D,t为时间;
步骤12:对步骤11的方程两边取时间导数,得微分方程初值问题
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mfrac>
<mrow>
<mi>d</mi>
<mi>x</mi>
</mrow>
<mrow>
<mi>d</mi>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mo>-</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msup>
<mi>F</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mi>F</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>x</mi>
<mo>(</mo>
<mn>0</mn>
<mo>)</mo>
<mo>=</mo>
<msub>
<mi>x</mi>
<mn>0</mn>
</msub>
<mo>,</mo>
<mi>t</mi>
<mo>&Element;</mo>
<mo>&lsqb;</mo>
<mn>0</mn>
<mo>,</mo>
<mo>+</mo>
<mi>&infin;</mi>
<mo>)</mo>
</mtd>
</mtr>
</mtable>
</mfenced>
当t=0,H(x,0)=0的解是初值x0,当t→+∞,H(x,+∞)=0的解就是F(x)=0的解;
步骤13:若计算中出现奇异点,则采用H1(x,t)=H(x,t)+η(t)G(x)计算,G(x)连续可微,且η(+∞)=0,G(x0)=0,则H1(x0,0)=H(x0,0)=0,H1(x*,+∞)=H(x*,+∞)=0,
且对应H1(x,t)=0的微分方程初值问题为
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mfrac>
<mrow>
<mi>d</mi>
<mi>x</mi>
</mrow>
<mrow>
<mi>d</mi>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mo>-</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msup>
<mi>F</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>&eta;</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>G</mi>
<mo>&prime;</mo>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>&lsqb;</mo>
<mi>F</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>+</mo>
<mi>&eta;</mi>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
<mi>G</mi>
<mo>(</mo>
<mi>x</mi>
<mo>)</mo>
<mo>&rsqb;</mo>
<mo>,</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>x</mi>
<mo>(</mo>
<mn>0</mn>
<mo>)</mo>
<mo>=</mo>
<msub>
<mi>x</mi>
<mn>0</mn>
</msub>
<mo>,</mo>
<mi>t</mi>
<mo>&Element;</mo>
<mo>&lsqb;</mo>
<mn>0</mn>
<mo>,</mo>
<mo>+</mo>
<mi>&infin;</mi>
<mo>&rsqb;</mo>
</mtd>
</mtr>
</mtable>
</mfenced>
步骤14:令η(t)=Ae-t,G(x)=x-x0,其中A是调节参数;采用四级四阶经典Runge-Kutta迭代求积公式求解得到Π,Q,Δ;
步骤15:完成迭代得到黏性应力和热传导Π,Q,Δ代入步骤1的守恒方程,迭代得守恒变量U的近似表达式Uh的系数Ai后,根据位置坐标求解任一点的流动特性参数ρ,u,T。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711110884.6A CN107729691B (zh) | 2017-11-13 | 2017-11-13 | 一种稀薄连续统一的气体流动特性数值模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711110884.6A CN107729691B (zh) | 2017-11-13 | 2017-11-13 | 一种稀薄连续统一的气体流动特性数值模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107729691A true CN107729691A (zh) | 2018-02-23 |
CN107729691B CN107729691B (zh) | 2020-05-22 |
Family
ID=61215049
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201711110884.6A Active CN107729691B (zh) | 2017-11-13 | 2017-11-13 | 一种稀薄连续统一的气体流动特性数值模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107729691B (zh) |
Cited By (14)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108920872A (zh) * | 2018-07-26 | 2018-11-30 | 上海交通大学 | 针对dsmc方法的bcp粒子定位实现方法及系统 |
CN109002629A (zh) * | 2018-08-01 | 2018-12-14 | 苏州慧德仿真技术有限公司 | 一种多相流模拟仿真的卷积神经网络及快速可视化方法 |
CN109489745A (zh) * | 2018-11-23 | 2019-03-19 | 宁波水表股份有限公司 | 一种基于数据迭代的流量计量方法 |
CN110705008A (zh) * | 2019-08-16 | 2020-01-17 | 北京航空航天大学 | 一种等离子体涡旋驱动装置的附加磁场位型的优化方法 |
CN110705007A (zh) * | 2019-08-16 | 2020-01-17 | 北京航空航天大学 | 一种等离子体涡旋发生器的效率评估方法 |
CN111259319A (zh) * | 2020-01-21 | 2020-06-09 | 扬州大学 | 一种确定稻田温室气体通量基准、周期与年际趋势的方法 |
CN111354086A (zh) * | 2018-12-24 | 2020-06-30 | 中国空气动力研究与发展中心超高速空气动力研究所 | 一种适用于dsmc方法网格位置属性判断的双向三维扫描方法 |
CN112685973A (zh) * | 2020-12-28 | 2021-04-20 | 中国航天空气动力技术研究院 | 非平衡高温气体混合物黏性系数计算方法及电子设备 |
CN113792432A (zh) * | 2021-09-15 | 2021-12-14 | 沈阳飞机设计研究所扬州协同创新研究院有限公司 | 基于改进型fvm-lbfs方法的流场计算方法 |
CN114065659A (zh) * | 2021-10-22 | 2022-02-18 | 中国人民解放军国防科技大学 | 一种两相流动界面演化模拟方法及系统 |
WO2022067498A1 (zh) * | 2020-09-29 | 2022-04-07 | 中南大学 | 一种气液相变的介观模拟方法 |
CN114398844A (zh) * | 2022-01-25 | 2022-04-26 | 南京航空航天大学 | 基于连续水膜流动的稳态防冰仿真方法 |
CN115618171A (zh) * | 2022-06-06 | 2023-01-17 | 北京理工大学 | 一种基于同伦算法的推进剂燃烧平衡产物求解方法 |
CN117393062A (zh) * | 2023-12-13 | 2024-01-12 | 上海交通大学四川研究院 | 刚性化学反应流回退自适应半隐半显耦合时间的模拟方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102890751A (zh) * | 2012-09-18 | 2013-01-23 | 天津空中代码工程应用软件开发有限公司 | 求解二维黎曼问题模拟亚音速无粘流的数值方法 |
CN104091065A (zh) * | 2014-07-03 | 2014-10-08 | 南京信息工程大学 | 一种求解浅水问题模拟间断水流数值的方法 |
-
2017
- 2017-11-13 CN CN201711110884.6A patent/CN107729691B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102890751A (zh) * | 2012-09-18 | 2013-01-23 | 天津空中代码工程应用软件开发有限公司 | 求解二维黎曼问题模拟亚音速无粘流的数值方法 |
CN104091065A (zh) * | 2014-07-03 | 2014-10-08 | 南京信息工程大学 | 一种求解浅水问题模拟间断水流数值的方法 |
Non-Patent Citations (3)
Title |
---|
HONG XIAO等: "A Unifed Framework for Modeling Continuum and Rarefed Gas Flows", 《SCIENTIFIC REPORTS》 * |
JUN-LIN WU等: "Numerical study on rarefied unsteady jet flow expanding into vacuum using the Gas-Kinetic Unified Algorithm", 《COMPUTERS AND FLUIDS》 * |
肖洪等: "稀薄气体动力学的非线性耦合本构方程理论及验证", 《航空学报》 * |
Cited By (19)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108920872A (zh) * | 2018-07-26 | 2018-11-30 | 上海交通大学 | 针对dsmc方法的bcp粒子定位实现方法及系统 |
CN109002629A (zh) * | 2018-08-01 | 2018-12-14 | 苏州慧德仿真技术有限公司 | 一种多相流模拟仿真的卷积神经网络及快速可视化方法 |
CN109489745A (zh) * | 2018-11-23 | 2019-03-19 | 宁波水表股份有限公司 | 一种基于数据迭代的流量计量方法 |
CN111354086A (zh) * | 2018-12-24 | 2020-06-30 | 中国空气动力研究与发展中心超高速空气动力研究所 | 一种适用于dsmc方法网格位置属性判断的双向三维扫描方法 |
CN111354086B (zh) * | 2018-12-24 | 2023-04-14 | 中国空气动力研究与发展中心超高速空气动力研究所 | 一种适用于dsmc方法网格位置属性判断的双向三维扫描方法 |
CN110705008A (zh) * | 2019-08-16 | 2020-01-17 | 北京航空航天大学 | 一种等离子体涡旋驱动装置的附加磁场位型的优化方法 |
CN110705007A (zh) * | 2019-08-16 | 2020-01-17 | 北京航空航天大学 | 一种等离子体涡旋发生器的效率评估方法 |
CN111259319B (zh) * | 2020-01-21 | 2023-03-14 | 扬州大学 | 一种确定稻田温室气体通量基准、周期与年际趋势的方法 |
CN111259319A (zh) * | 2020-01-21 | 2020-06-09 | 扬州大学 | 一种确定稻田温室气体通量基准、周期与年际趋势的方法 |
WO2022067498A1 (zh) * | 2020-09-29 | 2022-04-07 | 中南大学 | 一种气液相变的介观模拟方法 |
CN112685973A (zh) * | 2020-12-28 | 2021-04-20 | 中国航天空气动力技术研究院 | 非平衡高温气体混合物黏性系数计算方法及电子设备 |
CN113792432A (zh) * | 2021-09-15 | 2021-12-14 | 沈阳飞机设计研究所扬州协同创新研究院有限公司 | 基于改进型fvm-lbfs方法的流场计算方法 |
CN114065659A (zh) * | 2021-10-22 | 2022-02-18 | 中国人民解放军国防科技大学 | 一种两相流动界面演化模拟方法及系统 |
CN114398844A (zh) * | 2022-01-25 | 2022-04-26 | 南京航空航天大学 | 基于连续水膜流动的稳态防冰仿真方法 |
CN114398844B (zh) * | 2022-01-25 | 2023-04-07 | 南京航空航天大学 | 基于连续水膜流动的稳态防冰仿真方法 |
CN115618171A (zh) * | 2022-06-06 | 2023-01-17 | 北京理工大学 | 一种基于同伦算法的推进剂燃烧平衡产物求解方法 |
CN115618171B (zh) * | 2022-06-06 | 2023-10-24 | 北京理工大学 | 一种基于同伦算法的推进剂燃烧平衡产物求解方法 |
CN117393062A (zh) * | 2023-12-13 | 2024-01-12 | 上海交通大学四川研究院 | 刚性化学反应流回退自适应半隐半显耦合时间的模拟方法 |
CN117393062B (zh) * | 2023-12-13 | 2024-02-23 | 上海交通大学四川研究院 | 刚性化学反应流回退自适应半隐半显耦合时间的模拟方法 |
Also Published As
Publication number | Publication date |
---|---|
CN107729691B (zh) | 2020-05-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107729691A (zh) | 一种稀薄连续统一的气体流动特性数值模拟方法 | |
Guo et al. | Discrete unified gas kinetic scheme for all Knudsen number flows. II. Thermal compressible case | |
Hidalgo et al. | ADER schemes for nonlinear systems of stiff advection–diffusion–reaction equations | |
Sun et al. | Development of a vapor–liquid phase change model for volume-of-fluid method in FLUENT | |
Yang et al. | A hybrid lattice Boltzmann flux solver for simulation of viscous compressible flows | |
Cécora et al. | Differential Reynolds stress modeling for aeronautics | |
Yang et al. | Circular function-based gas-kinetic scheme for simulation of inviscid compressible flows | |
Feng et al. | Preserving the film coefficient as a parameter in the compact thermal model for fast electrothermal simulation | |
Furst | Numerical simulation of transitional flows with laminar kinetic energy | |
Im et al. | Large eddy simulation of coflow jet airfoil at high angle of attack | |
Hosseini et al. | Lattice Boltzmann advection-diffusion model for conjugate heat transfer in heterogeneous media | |
Ranjan et al. | Partially Averaged Navier Stokes simulation of turbulent heat transfer from a square cylinder | |
Liang et al. | Large eddy simulation of compressible turbulent channel flow with spectral difference method | |
Fürst et al. | Comparison of several models of the laminar/turbulent transition | |
Gawas et al. | Rayleigh-Bénard type natural convection heat transfer in two-dimensional geometries | |
Mittal et al. | Hysteresis in flow past a NACA 0012 airfoil | |
Giles et al. | Progress in adjoint error correction for integral functionals | |
Davidson et al. | Hybrid LES-RANS: Computation of the flow around a three-dimensional hill | |
Chen et al. | Cartesian grid method for gas kinetic scheme on irregular geometries | |
Fang et al. | Design optimization of unsteady airfoils with continuous adjoint method | |
Kim et al. | Turbulent Transition Prediction Using Large-Eddy Simulation with the Stability Theory | |
Orang et al. | Computational study of incompressible turbulent flows with method of characteristics | |
Waters et al. | Application of a dynamic LES model with an H-adaptive FEM for fluid and thermal processes | |
Zhang et al. | An improved MRT-LBM and investigation to the transition and periodicity of 2D lid-driven cavity flow with high Reynolds numbers | |
Bassi et al. | DISCONTINUOUS GALERKIN SOLUTION OF THE REYNOLDS-AVERAGED NAVIER–STOKES AND KL-KT-ω TRANSITION MODEL EQUATIONS |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |