CN109492240A - 一种基于二阶非线性本构模型的跨流域多尺度计算方法 - Google Patents

Abstract

本发明公开了一种基于二阶非线性本构模型的跨流域多尺度计算方法。该二阶非线性本构模型，从矩方法的模式出发，无需简化，直接采用Eu封闭方法对应力与热流输运方程左端进行封闭，然后在二阶本构模型求解中引进Rayleigh‑Onsager耗散函数将非线性代数方程组缩并成一个关于耗散函数R的方程，并采用最速下降方法求解单原子与双原子气体的R方程，最终得到气体对应的应力、热流与附加体积应力。弥补了NS方程在求解稀薄流及跨流域问题时失效的缺陷，为过渡流、稀薄流、及多尺度流动提供了有效的解决方案。在该方法的数值求解过程中采用了多种迭代方法，缩短收敛时间，提高了计算效率。

Description

一种基于二阶非线性本构模型的跨流域多尺度计算方法

技术领域

本发明涉及一种基于二阶非线性本构模型的跨流域多尺度计算方法，适用于不同努森数条件下的气体动力学数值计算方法，尤其适用于多尺度跨流域等计算状态。

背景技术

在20km～100km的高度范围被定义为临近空间，近年来引起了人们的广泛关注。在临近空间，气体的物理属性包括，密度、压力和温度，这些属性随高度不同而变化。由于临近空间介于连续稠密大气与稀薄近空间轨道之间，空间环境的特殊性决定了飞行器在此空域高速度、高机动、长航时飞行时必须考虑低密度大气与高温热环境对飞行器气动力、防隔热、通讯及控制的影响。除稀薄气体效应外，高马赫数飞行条件下气体头部滞止所产生的高温使得气体分子内部发生热力学能态激发、化学反应与辐射等复杂理化过程。因此，临近空间高超声速飞行器飞行环境存在时空多尺度与多物理效应耦合叠加，对传统空气动力学研究提出了新的课题与挑战。

为了准确模拟高超声速稀薄非平衡流动，几十年来国内外众多学者开展了大量研究，基于空气动力学与稀薄气体动力学的基本原理提出了多种类型的理论分析与数值仿真方法，包括考虑滑移边界条件修正的Navier-Stokes方程计算方法，蒙特卡洛直接数值模拟方法，线性化玻尔兹曼方程方法，玻尔兹曼方程BGK、ES-BKG或 Shakhov简化的模型方程方法，玻尔兹曼方程离散速度方法，以及Burnett方程、Grad 矩方法及其正则化理论和全局双曲化理论的扩展流体力学方法等。同时，近十年来兴起的UGKS和GKUA等统一气体动理论方法也为多尺度非平衡流动模拟提供了可发展的模拟手段。然而，基于原始或简化形式玻尔兹曼方程仿真所需的七维物理空间离散计算所带来的海量计算，以及稀薄气体流动粒子仿真类方法譬如DSMC方法的物理模拟本质所带来的分子自由程尺度网格需求、跨流域稀薄流仿真分子海量碰撞对计算、受碰撞频率约束的仿真分子自由运动时间等，使得这两类方法在高超声速飞行器跨流域稀薄流动仿真中的工程应用严重受限于计算机运算速度和内存的约束。由于物理化学模型的多样性与复杂性，在可见的时期内仍很难期望玻尔兹曼方程数值仿真方法与DSMC方法能够工程应用于复杂外形高超声速飞行器的热化学非平衡稀薄流动的仿真研究。另一方面，滑移物面边界条件虽然在物面附近努森层内一定程度上改善了宏观物理量分布准确度，然而其内部流场的应力与热流本构关系仍然是早已失效的Navier-Stokes-Fourier线性模型。采用一阶Maxwell-Smoluchowski (M/S)滑移边界条件的线性牛顿粘性/傅里叶热传导本构模型(NSF)在这些流域中的模拟能力十分有限，无法准确地刻画非平衡稀薄流中所表现出的强非线性特点。第三， Burnett方程与Grad十三矩方程由于与物理熵增定律相悖，数值计算的稳定性与结果正确性备受质疑。因此，上述稀薄气体动力学理论与计算方法近期难以克服的缺陷极大地降低了这些方法在临近空间复杂外形高超声速飞行器气动设计工程领域的可用性，建立可靠、高效的临近空间稀薄非平衡多尺度流动理论与预测解决方案就显得十分迫切。为了弥补线性NSF本构模型的精度不足，同时克服粒子类方法和 Boltzmann求解方法计算效率较低的缺陷，特别是为了充分利用近几十年来连续介质计算流体力学的宏大研究成果，近年来一套非线性耦合本构关系模型(Nonlinear coupled constitutive relations，NCCR)在由加拿大化学家B.C.Eu提出的严格满足热力学第二定律的广义流体力学方程(Generalized hydrodynamic equations,GHE)的基础上被发展起来。这套理论的核心在于巧妙地构建了一个形态定义非平衡态分布函数，作为一座桥梁把从非平衡态到平衡态演化的熵增特性和宏观非守恒量的耗散演化的过程紧密联系起来，使得这套理论从一开始就强制确保其满足热力学熵增定律及Boltzmann-H定理。韩国庆尚国立大学R.S.Myong教授在Eu方程基础上发展了一套有效的多维计算动力学模型，并为之提出一套行之有效的解耦求解算法，即非线性耦合本构关系模型(NCCR)。该模型由Eu的绝热假设和封闭假设条件下，通过对高阶非守恒量时间项和对流项的简化处理得到。通过解耦求解应力与热流二阶矩非线性耦合代数方程并有效地结合双曲守恒律控制方程，实现对流动的数值模拟。随后，该模型考虑了和分子内能非平衡与体积粘性有关的附加体积应力，通过引入附加体积应力这个高阶非守恒量演化方程，拓展到了双原子气体流动问题的模拟方面上来，并成功应用到了二维高超稀薄钝头绕流的模拟问题的研究当中。非线性耦合本构关系模型自提出之后得到了初步的研究和验证，表现出在单原子与双原子气体近平衡态连续流附近回归到NS方程解、非平衡流区域优于NS方程解的特点。目前，广义流体力学方程已经成功地被运用到一维激波结构、准一维库特流动、泊肃叶流努森层描述、声波吸收散布问题、二维平板边界层、稀薄圆柱绕流及方腔流动、三维高空高马赫数返回舱与类HTV-2飞行器绕流等问题研究当中，取得了显著的进展。其中，针对非线性本构关系耦合求解模型、间断有限元高精度算法及高精度物面滑移边界条件等问题，国内浙江大学陈伟芳、西北工业大学肖洪等均开展了初步的研究工作。已有研究成果表明，Eu方程与非线性本构关系理论具有巨大的研究潜力与应用价值。虽然前期广义流体力学方程的理论与计算研究已取得了一些突破，但临近空间高温稀薄气体的高超声速热化学非平衡流动，无论在理论计算还是试验模拟方面均是相当难于处理的多物理场耦合多尺度流动。为了克服多种计算方法的理论缺陷与技术不确定性，申请人首次提出了一种基于二维非线性本构模型的耦合计算方法，为直接求解三维多尺度跨流域问题提供了可靠的依据与广泛的应用前景。

发明内容

本发明的目的在于依据二阶非线性耦合模型，为多尺度跨流域问题提供一种的可靠有效的预测方法。

本发明提出的方法主要技术方案如下：

由Boltzmann-Curtiss方程分布函数在速度空间积分求矩，可以得单双原子三大守恒方程和非守恒量输运方程：

Boltzmann-Curtiss方程：

单原子气体守恒量控制方程：

非守恒量方程：

双原子气体守恒量控制方程：

非守恒量方程：

其中ρ,u,E,T,t,Π,Δ,Q,p,η,λ,I分别表示密度、速度矢量、总能、温度、时间、应力张量、附加正应力、热流、压力、粘性系数、热传导系数与单位矩阵，ψ⁽²⁾＝ψ⁽³⁾＝ψ⁽⁴⁾ψ⁽⁴⁾,ψ⁽⁵⁾,ψ⁽⁶⁾表示高阶矩，c_p,η_b,q(κ)分别表示定压比热、体粘性系数与非线性耗散因子。γ′＝(5-3γ)/2其中γ为比热比，由气体属性决定。

不进行简化，直接采用不同于Grad封闭的Eu封闭方法对应力与热流输运方程左端进行封闭：

ψ⁽²⁾＝ψ⁽³⁾＝ψ⁽⁴⁾＝0. (9)

最终可以得到一个更加精确的适用于单、双原子单一气体的一般化的流体动力学计算模型，无量纲化的演化本构方程被转换成如下形式，即构建获得二阶非线性本构模型：

其中c,f_b分别表示常系数与耗散函数、体积粘性和剪切粘性的比值(对于氮气f_b＝0.8)，Pr为普朗特数，ε＝1/(Pr(γ-1))，上标∧表示无量纲化，下标0表示NS方程应力与热流，以下公式如无特殊说明，符号均如前所述。

由于耦合本构方程的复杂性，可以通过引进Rayleigh-Onsager耗散函数将以上的方程缩并成一个方程：

其中

采用最速下降法，实现耗散函数的收敛，由于对于二阶非线性本构模型，可以表示成如下非线性方程组的一般形式：

F(x)＝0 (15)

其中考虑到剪切应力的无迹和对称性，是八维空间自变量，其中所有下标均表示对应张量的相应元素。构造目标函数G(x)

因此，求解非线性方程组(16)的解可转化求目标函数G(x)的极小值问题最速下降法的迭代思路是从目标函数的负梯度方向上选取合适的步长α_k，使得目标函数值最大程度的减少，从而逼近函数极小值的最优解，迭代过程如下：

如何选取合适的步长α_k，使得，

是最速下降法的关键。最直接的解析方法是构造有关步长α的单值函数，

然后通过对α求导解一元方程的根来寻找函数极小值。但是这种做法计算成本太高，本发明采用一维搜索近似得到每步较优步长来进行迭代，选取三个步长α₁＝0＜α₂＝α₃/2＜α₃，使得h(α₃)＜h(α₁)，利用，

h₁＝(h(α₂)-h(α₁))/α₂

h₂＝(h(α₃)-h(α₂))/(α₃-α₂) (20)

h₃＝(h₂-h₁)/α₃

构造二次多项式P(α)＝h(α₁)+h₁α+h₃α(α-α₂)在这三个点曲线拟合去近似单值函数 h(α)。在[α₁,α₃]定义近似步长使得作为该范围的最小值来代替h(α)的极小值，此时近似步长可作为k步迭代的合适步长α_k。近似步长由二次多项式P(α)临界点α₀＝(α₂-h₁/h₃)/2和α₃两者之间单值函数h(α)最小的一个确定。

本发明的有益效果是：

本发明的方法无需简化，直接采用Eu封闭方法对应力与热流输运方程左端进行封闭，然后在二阶本构模型求解中引进Rayleigh-Onsager耗散函数将非线性代数方程组缩并成一个关于耗散函数R的方程，并采用最速下降方法求解单原子与双原子气体的R方程，最终得到气体对应的应力、热流与附加体积应力。本发明提出的基于二阶非线性本构模型的跨流域多尺度计算方法方法能够有效的克服传统NS求解方法在过度流与稀薄流计算失效的缺陷，在各种计算状态下获得准确、稳定的连续流、过渡流与稀薄流数值计算结果。

附图说明

图1 返回舱几何模型图；

图2 计算域结构网格截面图；

图3 95km计算结果速度场对比图；

图4 95km计算结果压力场对比图；

图5 95km计算结果温度场对比图；

图6 物面压力系数分布曲线对比；

图7 物面摩阻系数分布曲线对比；

图8 物面热流系数分布曲线对比。

具体实施方式

下面结合附图详细描述本发明，本发明的目的和效果将变得更加明显。

本发明的基于二阶非线性本构模型的跨流域多尺度计算方法，是以Boltzman方程为基础，导出包含粘性应力、热传导、附加体积应力非守恒量的输运方程，并且根据Eu分布函数对由耗散引起的熵增，建立熵输运方程，从而构建获得二阶非线性本构模型，用于跨流域多尺度计算。

在具体的计算实现中包括流场初始化、无粘项离散与通量计算、非线性耦合本构模型迭代与粘性通量计算、新型边界条件赋值、时间步长计算等几个典型CFD过程。其中边界条件计算技术是按照结构化、模块化的思路进行的。该计算技术包含众多的空间离散格式、时间离散格式等先进的计算方法，以及限制器、熵修正、时间步计算、非物理解修正等多种选择，从而保证了计算功能的多种选择。

基于二阶非线性本构模型的跨流域多尺度计算技术，具体步骤与技术实现包括控制方程有限体积离散、无粘项采用了MUSCL插值与AUSMPW+通量差分方法、粘性项离散均采用二阶中心差分格式以及采用一阶Maxwell-Smoluchowski滑移跳跃边界条件求解出物面滑移速度和跳跃温度，为稀薄流动模拟提供准确的边界值。

采用本发明所述的方法针对Apollo返回舱着陆过程的绕流进行数值模拟与分析。来流参数具体如下：无穷远来流速度为U_∞＝8290m/s；采用等温壁模型，壁温为500K；当来流努森数为0.098时，来流压力为0.0145Pa；来流气体为空气，比热比1.4；气体常数296.9157；普朗特数0.72。在该流域气体已经非常稀薄，由于缺乏足够的气体粒子碰撞，非平衡效应非常强烈。图1为Apollo返回舱的计算模型。图 2为结构网格计算区域的截面图，由于来流方向物理量变化量比较大，进行了局部加密处理。从图3中至图5中可以看出，采用改进的二阶非线性本构模型耦合求解方法得到的计算结果，与传统NSF方法有显著的差异。其激波层厚度明显高于NSF的结果，更能体现稀薄气体效应。从图6压力系数、图7摩阻系数与图8热流系数的对比图中也可以观察发现，采用改进方法得到的计算结果与直接模拟蒙特卡洛方法(DSMC)为接近。

Claims (6)

1.一种基于二阶非线性本构模型的跨流域多尺度计算方法，其特征在于，以Boltzman方程为基础，导出包含粘性应力、热传导、附加体积应力非守恒量的输运方程，并且根据Eu分布函数对由耗散引起的熵增，建立熵输运方程，从而构建获得二阶非线性本构模型，用于跨流域多尺度计算。

2.一种基于二阶非线性本构模型的跨流域多尺度计算方法，其特征在于，包括如下步骤：

步骤(1)将针对单原子和双原子的Boltzmann统一动力学方程函数在速度空间积分求矩，最终得到质量、动量、能量三大守恒方程；

步骤(2)定义高阶速度矩，并对时间求导及分布函数积分，最终得到应力、热流、附加体积正应力非守恒量的输运方程；

步骤(3)对应力与热流输运方程左端进行封闭，并且引进Rayleigh-Onsager耗散函数将二阶非线性代数模型方程缩并成一个关于耗散函数R的方程；

步骤(4)在每一时间步计算中，采用最速下降方法，求解上述耗散函数R方程，根据所需的收敛准则得到最终的收敛解；并在此过程中计算非守恒量参数，应力、热流与附加体积应力，最终完成残差计算，得到流场结果。

3.如权利要求2所述的基于二阶非线性本构模型的跨流域多尺度计算方法，其特征在于，所述的步骤(3)具体为：

采用Eu封闭方法对应力与热流输运方程左端进行封闭：

ψ⁽²⁾＝ψ⁽³⁾＝ψ⁽⁴⁾＝0. (1)

得到适用于单、双原子单一气体的一般化的流体动力学计算模型，转换成如下形式，即构建获得二阶非线性本构模型：

其中c,f_b分别表示常系数与耗散函数、体积粘性和剪切粘性的比值，Pr为普朗特数，上标^表示无量纲化，下标0表示NS方程应力与热流,且ε＝1/(Pr(γ-1))；

引进Rayleigh-Onsager耗散函数将以上的方程缩并成一个方程：

其中

4.如权利要求2所述的基于二阶非线性本构模型的跨流域多尺度计算方法，其特征在于，所述的步骤(4)中采用最速下降法，实现耗散函数的收敛，具体如下：

将二阶非线性本构模型表示成如下非线性方程组的一般形式：

F(x)＝0 (7)

其中是八维空间自变量，构造目标函数G(x)

将求解非线性方程组(8)的解转化为求目标函数G(x)的极小值问题从目标函数的负梯度方向上选取合适的步长α_k，使得目标函数值最大程度的减少，从而逼近函数极小值的最优解，迭代过程如下：

选取合适的步长α_k，使得：

5.根据权利要求4所述的基于二阶非线性本构模型的跨流域多尺度计算方法，其特征在于，所述的步长α_k的选取方法为：

构造有关步长α的单值函数，

然后通过对α求导解一元方程的根来寻找函数极小值。

6.根据权利要求4所述的基于二阶非线性本构模型的跨流域多尺度计算方法，其特征在于，所述的步长α_k的选取方法为：

采用一维搜索近似得到每步较优步长来进行迭代：选取三个步长α₁＝0＜α₂＝α₃/2＜α₃，使得h(α₃)＜h(α₁)，利用，

构造二次多项式P(α)＝h(α₁)+h₁α+h₃α(α-α₂)在这三个点曲线拟合近似单值函数h(α)；在[α₁,α₃]定义近似步长使得作为该范围的最小值来代替h(α)的极小值，此时近似步长可作为k步迭代的合适步长α_k；近似步长由二次多项式P(α)临界点α₀＝(α₂-h₁/h₃)/2和α₃两者之间单值函数h(α)最小的一个确定。