基于模态衰减的人工黏性方法及系统与流程

1.本发明涉及计算流体力学的技术领域，具体地，涉及基于模态衰减的人工黏性方法及系统，尤其涉及基于模态衰减的人工黏性模型。


背景技术：

2.在计算流体力学中，高阶格式对于复杂流动的数值模拟具有重要意义。以旋翼流动的模拟为例，对于该类流动，翼梢涡是其主要流场结构之一，采用通常的结构网格计算。而高阶格式对于这类问题则显示出了绝对的优势，采用高阶格式得到的旋翼气动特性和噪声的计算精度都有了显著的提高，并且其计算效率也比在低阶格式的基础上加密网格的效率要高。
3.高阶格式面临的一大挑战是：在可压缩流动中，以极小的数值耗散来分辨小尺度脉动，在光滑区保持格式的高精度特性不被破坏，同时光滑地捕捉流场中出现的各种间断。混合使用不同的高阶格式来计算，可以解决这一难题。混合格式是将在光滑流场中具有低耗散特性的高阶格式，和可以识别流场间断的激波捕捉格式相组合，在保证低耗散和高精度的同时，也能够准确地捕捉流场中的激波。激波捕捉方法一直是cfd研究的重点和难点，在有限差分和有限体积框架下，研究者们发展了大量经典的激波捕捉方法，包括限制器类方法、重构类方法、人工黏性方法和混合方法。
4.人工黏性方法是通过添加人工黏性项来抑制非物理振荡，其关键点是构造合适的光滑度量因子和选取适当的人工黏性系数。一种思路是采用某些变量(比如密度、速度等)的导数及其派生变量来作为光滑度量因子。该类方法特点是鲁棒性较强，并且由于针对不同的流动结构设计了不同的度量因子，因此能够灵活的适用于复杂流动。然而，对于更高阶数的导数可以更好的保证光滑区域的计算精度，但这会使得算法复杂度和计算代价大大增加。
5.而通过模态系数的衰减速度来构造人工黏性，方法构造形式简单，易于理解和实现。其基本思想是光滑流动区域高阶模态系数，应该迅速衰减至较小量级。激波捕捉方法可以在平滑区域保持相应计算格式的精度，对于非平滑的流动，该方法可以很好地抑制激波间断处的伪振荡。


技术实现要素：

6.针对现有技术中的缺陷，本发明的目的是提供一种基于模态衰减的人工黏性方法及系统。
7.根据本发明提供的一种基于模态衰减的人工黏性方法，所述方法包括如下步骤：
8.步骤s1：根据所模拟的具体流动，建立控制方程；
9.步骤s2：由模态衰减的间断指示器，构造模态分量函数然后建立关于衰减率τ的目标函数，并使用最小二乘的方法计算衰减率τ；
10.步骤s3：使周围网格守恒变量q组成新的守恒变量q
ex
代替原守恒变量q，同时由衰
减率τ得到新的人工黏性模型，并计算出黏性系数μ；
11.步骤s4：通过得到的守恒变量q
ex
和黏性系数μ，得到粘性通量g，基于间断有限元的理论求出随时间t变化的守恒变量q，进而得到整个流场的结果。
12.优选地，所述步骤s1包括如下步骤：
13.步骤s1.1：求解一般的双曲方程：
[0014][0015]
其中t为时间变量，q为守恒变量，f为对流通量，g为粘性通量；
[0016]
步骤s1.2：粘性通量g表示为：
[0017][0018]
其中μ为人工黏性系数；
[0019]
步骤s1.3：将守恒变量q分解为时间项和空间项：
[0020][0021]
其中n
p
和r分别表示单元网格中节点个数和坐标向量，m为单元网格中节点的序号，模态分量函数，为一般正交基函数；
[0022]
步骤s1.4：根据模态衰减的间断指示器，构造模态分量函数为：
[0023][0024]
其中c为已知参数，τ为衰减率。
[0025]
优选地，所述步骤s2包括如下步骤：
[0026]
步骤s2.1：构造一个最优的模态衰减系数bm：
[0027][0028]
步骤s2.2：将原模态分量函数与最优模态衰减系数bm组合，构造新的模态系数为
[0029][0030]
其中，表示为在网格单元k中取守恒变量q的标准二范数l2；
[0031]
步骤s2.3：通过如下公式，模态系数依次单调：
[0032][0033]
其中，n表示为网格节点的一组序号；
[0034]
步骤s2.4：使用最小二乘的方法计算衰减率τ，目标函数选取：
[0035][0036]
优选地，所述步骤s3包括如下步骤：
[0037]
步骤s3.1：将该网格单元守恒变量q和左右两个网格单元的守恒变量q
l
、qr构造一个新的守恒变量q
ex
：
[0038]qex
(r,t)＝γ
lql
(r,t) γq(r,t) γrqr(r,t)
ꢀꢀꢀ
(9)
[0039]
其中，γ
l
，γ，γr表示网格单元守恒变量对应的加权系数；
[0040]
步骤s3.2：构造的守恒变量值光滑变化，加权系数γ
l
，γ，γr选择如下所示：
[0041][0042]
其中参数r
*
为：
[0043][0044]
其中，p为计算的精度，p
*
和p
max
表示计算中所需的精度参数；
[0045]
步骤s3.3：在计算时如果耗散过度，会破坏高阶格式的精度，计算黏性系数μ的方法为：
[0046][0047]
其中μ
max
为最大的黏性系数，h为网格尺度，p为计算的精度。
[0048]
本发明还提供一种基于模态衰减的人工黏性系统，所述系统包括如下模块：
[0049]
模块m1：根据所模拟的具体流动，建立控制方程；
[0050]
模块m2：由模态衰减的间断指示器，构造模态分量函数然后建立关于衰减率τ的目标函数，并使用最小二乘的系统计算衰减率τ；
[0051]
模块m3：使周围网格守恒变量q组成新的守恒变量q
ex
代替原守恒变量q，同时由衰减率τ得到新的人工黏性模型，并计算出黏性系数μ；
[0052]
模块m4：通过得到的守恒变量q
ex
和黏性系数μ，得到粘性通量g，基于间断有限元的理论求出随时间t变化的守恒变量q，进而得到整个流场的结果。
[0053]
优选地，所述模块m1包括如下模块：
[0054]
模块m1.1：求解一般的双曲方程：
[0055][0056]
其中t为时间变量，q为守恒变量，f为对流通量，g为粘性通量；
[0057]
模块m1.2：粘性通量g表示为：
[0058]
[0059]
其中μ为人工黏性系数；
[0060]
模块m1.3：将守恒变量q分解为时间项和空间项：
[0061][0062]
其中n
p
和r分别表示单元网格中节点个数和坐标向量，m为单元网格中节点的序号，模态分量函数，为一般正交基函数；
[0063]
模块m1.4：根据模态衰减的间断指示器，构造模态分量函数为：
[0064][0065]
其中c为已知参数，τ为衰减率。
[0066]
优选地，所述模块m2包括如下模块：
[0067]
模块m2.1：构造一个最优的模态衰减系数bm：
[0068][0069]
模块m2.2：将原模态分量函数与最优模态衰减系数bm组合，构造新的模态系数为
[0070][0071]
其中，表示为在网格单元k中取守恒变量q的标准二范数l2；
[0072]
模块m2.3：通过如下公式，模态系数依次单调：
[0073][0074]
其中，n表示为网格节点的一组序号；
[0075]
模块m2.4：使用最小二乘的系统计算衰减率τ，目标函数选取：
[0076][0077]
优选地，所述模块m3包括如下模块：
[0078]
模块m3.1：将该网格单元守恒变量q和左右两个网格单元的守恒变量q
l
、qr构造一个新的守恒变量q
ex
：
[0079]qex
(r,t)＝γ
lql
(r,t) γq(r,t) γrqr(r,t)
ꢀꢀꢀ
(9)
[0080]
其中，γ
l
，γ，γr表示网格单元守恒变量对应的加权系数；
[0081]
模块m3.2：构造的守恒变量值光滑变化，加权系数γ
l
，γ，γr选择如下所示：
[0082][0083]
其中参数r
*
为：
[0084][0085]
其中，p为计算的精度，p
*
和p
max
表示计算中所需的精度参数；
[0086]
模块m3.3：在计算时如果耗散过度，会破坏高阶格式的精度，计算黏性系数μ的系统为：
[0087][0088]
其中μ
max
为最大的黏性系数，h为网格尺度，p为计算的精度。
[0089]
与现有技术相比，本发明具有如下的有益效果：
[0090]
1、该人工黏性模型可以很好的捕捉流场中的激波位置，从而有效抑制激波间断处的伪振荡，减小计算误差；
[0091]
2、通过结合当前网格及其相邻网格的变量，在不增加过多计算量下，增强用于估计衰减率的信息，适用于低阶的计算。
附图说明
[0092]
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
[0093]
图1为本发明方法完整算法流程图；
[0094]
图2为本发明方法算例结果图。
具体实施方式
[0095]
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
[0096]
本发明的目的是将应用模态衰减模型，包含周围网格单元的信息，扩展到高阶计算精度，构造新的人工黏性。本发明使用了形式简单的模态衰减模型，通过加权函数将周围流场信息组合。本发明不仅能保证准确捕捉流场中的激波位置，也适用与低阶格式的计算。
[0097]
本发明提出了一种基于模态衰减的人工黏性模型，其流程图如附图1和算例结果如附图2所示，它包含以下四个步骤：
[0098]
步骤s1：根据所模拟的具体流动，建立控制方程，具体包括如下子步骤：
[0099]
步骤s1.1：求解一般的双曲方程：
[0100][0101]
其中t为时间变量，q为守恒变量，f为对流通量，g为粘性通量；
[0102]
步骤s1.2：粘性通量g表示为：
[0103][0104]
其中μ为人工黏性系数；
[0105]
步骤s1.3：将守恒变量q分解为时间项和空间项：
[0106][0107]
其中n
p
和r分别表示单元网格中节点个数和坐标向量，m为单元网格中节点的序号，模态分量函数，为一般正交基函数；
[0108]
步骤s1.4：根据模态衰减的间断指示器，构造模态分量函数为：
[0109][0110]
其中c为已知参数，τ为衰减率。
[0111]
步骤s2：由模态衰减的间断指示器，构造模态分量函数然后建立关于衰减率τ的目标函数，并使用最小二乘的方法计算衰减率τ，具体包括如下子步骤：
[0112]
步骤s2.1：构造一个最优的模态衰减系数bm：
[0113][0114]
步骤s2.2：将原模态分量函数与最优模态衰减系数bm组合，构造新的模态系数为
[0115][0116]
其中，表示为在网格单元k中取守恒变量q的标准二范数l2；
[0117]
步骤s2.3：通过如下公式，模态系数依次单调：
[0118][0119]
其中，n表示为网格节点的一组序号；
[0120]
步骤s2.4：使用最小二乘的方法计算衰减率τ，目标函数选取：
[0121][0122]
步骤s3：使周围网格守恒变量q组成新的守恒变量q
ex
代替原守恒变量q，同时由衰减率τ得到新的人工黏性模型，并计算出黏性系数μ，具体计算步骤如下：
[0123]
步骤s3.1：将该网格单元守恒变量q和左右两个网格单元的守恒变量q
l
、qr构造一个新的守恒变量q
ex
：
[0124]qex
(r,t)＝γ
lql
(r,t) γq(r,t) γrqr(r,t)
ꢀꢀꢀ
(9)
[0125]
其中，γ
l
，γ，γr表示网格单元守恒变量对应的加权系数；
[0126]
步骤s3.2：构造的守恒变量值光滑变化，加权系数γ
l
，γ，γr选择如下所示：
[0127][0128]
其中参数r
*
为：
[0129][0130]
其中，p为计算的精度，p
*
和p
max
表示计算中所需的精度参数；
[0131]
步骤s3.3：在计算时如果耗散过度，会破坏高阶格式的精度，计算黏性系数μ的方法为：
[0132][0133]
其中μ
max
为最大的黏性系数，h为网格尺度，p为计算的精度。
[0134]
步骤s4：将得到的守恒变量q
ex
和黏性系数μ，代入方程(2)可以得到粘性通量g。之后基于间断有限元的理论可以求出随时间t变化的守恒变量q，进而可以得到整个流场的结果。
[0135]
本发明还提供一种基于模态衰减的人工黏性系统，所述系统包括如下模块：
[0136]
模块m1：根据所模拟的具体流动，建立控制方程。
[0137]
模块m1.1：求解一般的双曲方程：
[0138][0139]
其中t为时间变量，q为守恒变量，f为对流通量，g为粘性通量；
[0140]
模块m1.2：粘性通量g表示为：
[0141][0142]
其中μ为人工黏性系数；
[0143]
模块m1.3：将守恒变量q分解为时间项和空间项：
[0144][0145]
其中n
p
和r分别表示单元网格中节点个数和坐标向量，m为单元网格中节点的序号，模态分量函数，为一般正交基函数；
[0146]
模块m1.4：根据模态衰减的间断指示器，构造模态分量函数为：
[0147][0148]
其中c为已知参数，τ为衰减率。
[0149]
模块m2：由模态衰减的间断指示器，构造模态分量函数然后建立关于衰减率τ的目标函数，并使用最小二乘的系统计算衰减率τ。
[0150]
模块m2.1：构造一个最优的模态衰减系数bm：
[0151][0152]
模块m2.2：将原模态分量函数与最优模态衰减系数bm组合，构造新的模态系数为
[0153][0154]
其中，表示为在网格单元k中取守恒变量q的标准二范数l2；
[0155]
模块m2.3：通过如下公式，模态系数依次单调：
[0156][0157]
其中，n表示为网格节点的一组序号；
[0158]
模块m2.4：使用最小二乘的系统计算衰减率τ，目标函数选取：
[0159][0160]
模块m3：使周围网格守恒变量q组成新的守恒变量q
ex
代替原守恒变量q，同时由衰减率τ得到新的人工黏性模型，并计算出黏性系数μ。
[0161]
模块m3.1：将该网格单元守恒变量q和左右两个网格单元的守恒变量q
l
、qr构造一个新的守恒变量q
ex
：
[0162]qex
(r,t)＝γ
lql
(r,t) γq(r,t) γrqr(r,t)
ꢀꢀꢀ
(9)
[0163]
其中，γ
l
，γ，γr表示网格单元守恒变量对应的加权系数；
[0164]
模块m3.2：构造的守恒变量值光滑变化，加权系数γ
l
，γ，γr选择如下所示：
[0165][0166]
其中参数r
*
为：
[0167][0168]
其中，p为计算的精度，p
*
和p
max
表示计算中所需的精度参数；
[0169]
模块m3.3：在计算时如果耗散过度，会破坏高阶格式的精度，计算黏性系数μ的系统为：
[0170]
[0171]
其中μ
max
为最大的黏性系数，h为网格尺度，p为计算的精度。
[0172]
模块m4：通过得到的守恒变量q
ex
和黏性系数μ，得到粘性通量g，基于间断有限元的理论求出随时间t变化的守恒变量q，进而得到整个流场的结果。
[0173]
算例：本算例是检验基于模态衰减构造高阶的激波捕捉方法数值模拟性能的。计算区域为-5＜x＜5，计算网格量为200个，格式的精度阶数p＝3，边界条件为远场边条，计算终止时间为t＝1.8。初始的密度ρ、速度u、压强p为：
[0174][0175]
图2为本实施例模拟结果的密度沿x轴分布曲线。图2中可以看出，本发明清晰地捕捉了流场中的激波(x＜0.5)，说明本发明能够稳定、高精度地捕捉流场中的强间断；对于多尺度湍流结构嵌套区域0.5＜x＜2.5，本发明能够较好地捕捉几个极值点，说明本发明对多尺度湍流结构的分辨能力较好。
[0176]
本发明提出了一种基于模态衰减的人工黏性模型，而通过模态系数的衰减速度来构造人工黏性，方法构造形式简单，易于理解和实现。包含周围网格单元的信息，扩展到高阶计算精度，构造新的人工黏性。使用了形式简单的模态衰减模型，通过加权函数将周围流场信息组合。不仅能保证准确捕捉流场中的激波位置，也可以提高计算的鲁棒性。构造的人工黏性模型可以在平滑区域保持相应计算格式的精度，对于非平滑的流动，该方法可以很好地抑制激波间断处的伪振荡。
[0177]
本发明的核心算法是将应用模态衰减模型，包含周围网格单元的信息，扩展到高阶计算精度，构造新的人工黏性。使用了形式简单的模态衰减模型，通过加权函数将周围流场信息组合。在低阶格式的计算时，不会因为模板信息过少而对流场激波捕捉的不够准确，因而构造的人工黏性可以很好的适用于低阶格式的计算。如果对本发明的各种改动和变形不脱离本发明范围，仍然属于本发明的权利要求和同等技术范围之内。
[0178]
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统及其各个装置、模块、单元以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统及其各个装置、模块、单元以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同功能。所以，本发明提供的系统及其各项装置、模块、单元可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置、模块、单元也可以视为硬件部件内的结构；也可以将用于实现各种功能的装置、模块、单元视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
[0179]
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本技术的实施例和实施例中的特征可以任意相互组合。