CN112685978B - 一种适用于五次样条重构格式的自适应人工粘性控制方法 - Google Patents

一种适用于五次样条重构格式的自适应人工粘性控制方法 Download PDF

Info

Publication number
CN112685978B
CN112685978B CN202110264089.2A CN202110264089A CN112685978B CN 112685978 B CN112685978 B CN 112685978B CN 202110264089 A CN202110264089 A CN 202110264089A CN 112685978 B CN112685978 B CN 112685978B
Authority
CN
China
Prior art keywords
interface
physical quantity
unit
state value
artificial viscosity
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.)
Active
Application number
CN202110264089.2A
Other languages
English (en)
Other versions
CN112685978A (zh
Inventor
黄文锋
陈坚强
袁先旭
徐国亮
涂国华
刘旭亮
李辰
郭启龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202110264089.2A priority Critical patent/CN112685978B/zh
Publication of CN112685978A publication Critical patent/CN112685978A/zh
Application granted granted Critical
Publication of CN112685978B publication Critical patent/CN112685978B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Feedback Control In General (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供了一种适用于五次样条重构格式的自适应人工粘性控制方法，包括：步骤1、读取流场数据，求解五次样条重构方程组得到单元界面处物理量的状态值及其一阶导数；步骤2、计算出单元界面处的二阶至五阶导数；步骤3、根据波数识别方法，求出单元界面处流场的等效无量纲波数，从而确定该处流场的人工粘性系数；步骤4、根据单元界面处物理量的状态值计算无粘数值通量，并根据人工粘性系数添加六阶人工粘性项，最后采用相对应的时间离散方法在时间上进行推进。将局部瞬时流场特征与格式的解析能力结合起来，通过判断局部瞬时的流场特征是否在格式的可解析范围内，从而确定格式的人工粘性系数的取值，使格式的人工粘性控制更为合理。

Description

一种适用于五次样条重构格式的自适应人工粘性控制方法
技术领域
本发明涉及计算流体力学中数值计算方法领域，特别涉及一种适用于五次样条重构格式的自适应人工粘性控制方法。
背景技术
在计算流体力学中，数值格式的耗散起着非常重要的作用。数值耗散越大，计算越稳定，但对流场的刻画能力越弱；数值格式的耗散越小，流场的刻画能力就越强，计算越不稳定。因此如何调节格式的数值耗散一直是人们研究的重点。传统调节耗散的方法是调节格式的总体耗散，这种调节耗散的方式并不能完全满足人们的需求。理想的数值格式应当满足以下条件：在流动光滑区域数值耗散尽量小以保证格式精度和分辨率，在流动非光滑区或者梯度较大的区域数值耗散要足够大以保证计算稳定。因此发展根据局部流场信息确定格式耗散大小的自适应耗散方法具有广泛的科研以及应用前景。
发明内容
针对现有技术中存在的问题，提供了一种适用于五次样条重构格式的自适应人工粘性方法；该方法将流场的局部流动特征同格式的色散特性结合起来，由两者共同决定局部流场需要的人工粘性系数的大小，实现人工粘性系数的自适应改变。首先根据波数识别方法，由表征局部流场特征的密度的高阶导数求出该处流场的局部等效无量纲波数，再将等效无量纲波数与格式的修正无量纲波数进行对比。如果等效无量纲波数在格式的完全解析范围内，则人工粘性系数为零或者取一个很小的值；如果等效无量纲波数不在格式的完全解析范围内，则人工粘性系数根据一定的规则逐渐增大，直至最大值。这个判定过程即为自适应人工粘性判据。这种自适应人工粘性方法改变了以往的全流场使用统一的人工粘性系数且人工粘性系数的大小根据经验确定的方式，使得人工粘性系数的取值更为合理。
本发明采用的技术方案如下：一种适用于五次样条重构格式的自适应人工粘性控制方法，包括以下步骤：
步骤1、读取流场数据，求解五次样条重构方程组得到单元界面处物理量的状态值及其一阶导数；
步骤2、根据单元界面处物理量的状态值及一阶导数，计算出单元界面处的二阶至五阶导数；
步骤3、根据波数识别方法，由单元界面处流场密度的高阶导数求出单元界面处流场的等效无量纲波数，代入自适应人工粘性判据确定该处流场的人工粘性系数；
步骤4、根据单元界面处物理量的状态值计算无粘数值通量，并根据人工粘性系数添加六阶人工粘性项，最后采用相对应的时间离散方法在时间上进行推进。
进一步的，在结构网格有限体积方法的基本框架下，不计质量力和源项，在直坐标系下欧拉方程为：
Figure 159716DEST_PATH_IMAGE001
其中，
Figure 957907DEST_PATH_IMAGE002
为守恒变量，
Figure 563332DEST_PATH_IMAGE003
、
Figure 763369DEST_PATH_IMAGE004
、
Figure 163258DEST_PATH_IMAGE005
为直角坐标系
Figure 448746DEST_PATH_IMAGE006
下的无粘通量，具体表 达式为：
Figure 592282DEST_PATH_IMAGE007
其中，
Figure 318930DEST_PATH_IMAGE008
、
Figure 14353DEST_PATH_IMAGE009
、
Figure 662503DEST_PATH_IMAGE010
分别为流场的密度、速度矢量以及流场压力，
Figure 734365DEST_PATH_IMAGE011
为总能，其具体表达式为：
Figure 846677DEST_PATH_IMAGE012
其中
Figure 853947DEST_PATH_IMAGE013
为气体常数；
将欧拉方程在三维结构网格的控制体单元
Figure 848448DEST_PATH_IMAGE014
上进行积分，得到：
Figure 599366DEST_PATH_IMAGE015
其中
Figure 362923DEST_PATH_IMAGE016
为控制体单元的体积，
Figure 538165DEST_PATH_IMAGE017
是控制体面上的外法线 方向，
Figure 19962DEST_PATH_IMAGE018
是守恒变量在控制体单元上的平均值：
Figure 308992DEST_PATH_IMAGE019
Figure 458214DEST_PATH_IMAGE020
为通量张量，
Figure 807286DEST_PATH_IMAGE021
即为无粘通量，
Figure 776379DEST_PATH_IMAGE022
为控制 体的六个面，以
Figure 869100DEST_PATH_IMAGE023
表示为控制体的面标号，则面积分项写为控制体各面积分的加和：
Figure 607249DEST_PATH_IMAGE024
每个面上
Figure 127223DEST_PATH_IMAGE025
的具体表达式为：
Figure 583613DEST_PATH_IMAGE026
其中
Figure 948866DEST_PATH_IMAGE027
分别为单元界面处的密度、三个方向上的速度、 压力以及总能的状态值；
在求解过程中根据单元的平均值重构出单元界面处的物理量的状态值，五次样条重构方程组为：
Figure 479204DEST_PATH_IMAGE028
其中，
Figure 763555DEST_PATH_IMAGE029
Figure 379344DEST_PATH_IMAGE030
Figure 545359DEST_PATH_IMAGE031
Figure 930204DEST_PATH_IMAGE032
下标
Figure 651035DEST_PATH_IMAGE033
分别代表着
Figure 488541DEST_PATH_IMAGE034
三个方向，
Figure 461176DEST_PATH_IMAGE035
和
Figure 294003DEST_PATH_IMAGE036
分别代表
Figure 326681DEST_PATH_IMAGE037
单元和
Figure 979379DEST_PATH_IMAGE038
单元从单元左侧界面中心至单元右侧界面中心的距离，
Figure 755706DEST_PATH_IMAGE039
分别指的
Figure 443039DEST_PATH_IMAGE040
五个守恒物理量，
Figure 646618DEST_PATH_IMAGE041
和
Figure 786613DEST_PATH_IMAGE042
分别指在
Figure 428946DEST_PATH_IMAGE043
界面处第
Figure 846152DEST_PATH_IMAGE044
个物理量的界面状态值及其一阶导数，
Figure 345267DEST_PATH_IMAGE045
和
Figure 844994DEST_PATH_IMAGE046
分别指
Figure 87756DEST_PATH_IMAGE047
单 元和
Figure 359469DEST_PATH_IMAGE048
单元处的第
Figure 29485DEST_PATH_IMAGE049
个物理量的单元平均值，通过求解重构方程，可得到单元界面
Figure 19437DEST_PATH_IMAGE050
处物理量的状态值
Figure 331470DEST_PATH_IMAGE051
及其一阶导数
Figure 192110DEST_PATH_IMAGE052
；上述五次样条重构方 程组、界面状态值及一阶导数均为
Figure 33027DEST_PATH_IMAGE053
方向，
Figure 510276DEST_PATH_IMAGE054
方向和
Figure 94841DEST_PATH_IMAGE055
方向需要采用相同的方法求解。
进一步的，所述步骤2具体包括，根据单元界面
Figure 341145DEST_PATH_IMAGE056
处物理量的状态值
Figure 352964DEST_PATH_IMAGE057
及一阶导数
Figure 317509DEST_PATH_IMAGE042
计算单元中心
Figure 440186DEST_PATH_IMAGE058
处物理量的一阶至四阶导数：
Figure 540997DEST_PATH_IMAGE059
从而得到单元界面
Figure 989296DEST_PATH_IMAGE060
处物理量的二阶至五阶导数：
Figure 172628DEST_PATH_IMAGE061
Figure 771099DEST_PATH_IMAGE062
方向和
Figure 116630DEST_PATH_IMAGE055
方向采用相同的方法可求解得到高阶导数。
进一步的，波数识别方法为：
Figure 345617DEST_PATH_IMAGE063
其中
Figure 143809DEST_PATH_IMAGE064
即为单元界面
Figure 483654DEST_PATH_IMAGE065
处的等效无量纲波数，
Figure 683692DEST_PATH_IMAGE066
为一个 小量，防止分母为零，此处取
Figure 83580DEST_PATH_IMAGE067
，
Figure 634647DEST_PATH_IMAGE068
方向和
Figure 778184DEST_PATH_IMAGE055
方向采用相同的方法可求得两个 方向上的等效无量纲波数
Figure 832727DEST_PATH_IMAGE069
、
Figure 669096DEST_PATH_IMAGE070
。
进一步的，自适应人工粘性判据为：
Figure 441880DEST_PATH_IMAGE071
其中，
Figure 857949DEST_PATH_IMAGE072
为即为自适应人工粘性系数，
Figure 32579DEST_PATH_IMAGE054
方向和
Figure 771340DEST_PATH_IMAGE055
方向采用相 同的方法可求得两个方向上的
Figure 296999DEST_PATH_IMAGE073
、
Figure 782338DEST_PATH_IMAGE074
。
进一步的，所述步骤4具体包括：将单元界面处物理量的状态值带入控制方程计算无粘数值通量，无粘通量计算过程为：
Figure 811474DEST_PATH_IMAGE075
其中，
Figure 989646DEST_PATH_IMAGE076
界面上无粘通量
Figure 471443DEST_PATH_IMAGE077
为：
Figure 760473DEST_PATH_IMAGE078
其中，
Figure 378536DEST_PATH_IMAGE079
为
Figure 727609DEST_PATH_IMAGE080
界面处的外法向速度；
Figure 696702DEST_PATH_IMAGE081
为六阶人工粘性项，在计算中需要在
Figure 320581DEST_PATH_IMAGE082
后面添加此项来抑 制计算中的高频振荡其具体表达式为:
Figure 199676DEST_PATH_IMAGE083
其中，
Figure 844283DEST_PATH_IMAGE084
为自适应人工粘性系数；
Figure 176039DEST_PATH_IMAGE085
为尺度因子，对于各向 同性模型，其计算公式为：
Figure 931505DEST_PATH_IMAGE086
其中
Figure 131018DEST_PATH_IMAGE087
为界面处的声速，
Figure 946527DEST_PATH_IMAGE088
为界面
Figure 765579DEST_PATH_IMAGE089
处的面积矢 量；采用相同的方法可求得
Figure 590315DEST_PATH_IMAGE090
方向和
Figure 178422DEST_PATH_IMAGE055
方向的无粘通量。
与现有技术相比，采用上述技术方案的有益效果为：本发明的方案中人工粘性系数不再是全场统一的值，而是由自适应人工粘性判据根据局部流场的特征来确定,自适应人工粘性判据局部瞬时的流场特征和格式的色散特性共同来决定，此种方法将局部瞬时流场特征与格式的解析能力结合起来，通过判断局部瞬时的流场特征是否在格式的可解析范围内，从而确定格式的人工粘性系数的取值，使格式的人工粘性控制更为合理。
附图说明
图1是本发明的适用于五次样条重构格式的自适应人工粘性方法流程图。
图2是本发明中的等效无量纲波数与修正无量纲波数对比图。
图3是本发明中自适应人工粘性系数随波数变化曲线图。
图4是本发明中与自适应人工粘性系数对应的数值格式耗散曲线。
图5是采用一种固定人工粘性的二维黎曼问题密度等值线图。
图6是采用自适应人工粘性的二维黎曼问题密度等值线图。
图7是采用另一固定人工粘性的二维黎曼问题密度等值线图。
具体实施方式
下面结合附图对本发明做进一步描述。
本发明提出了一种适用于五次样条重构格式的自适应人工粘性方法，包括下述理论依据：
（1）波数识别方法
考虑傅里叶模态
Figure 899254DEST_PATH_IMAGE091
并假设网格均匀，其中
Figure 205601DEST_PATH_IMAGE092
为波数， 则傅里叶模态的一阶至五阶导数分别为：
Figure 302870DEST_PATH_IMAGE093
(1)
由式 (1)可得：
Figure 745484DEST_PATH_IMAGE094
(2)
将式（2）进行无量纲化，为防止导数为零造成分母为零，等效无量纲波数可写出如下形式：
Figure 637217DEST_PATH_IMAGE095
(3)
其中
Figure 227598DEST_PATH_IMAGE096
为一小量，防止分母为零。
（2）波数识别在五次样条重构格式中的验证
对于五次样条重构格式，采用傅里叶分析，单元界面
Figure 3924DEST_PATH_IMAGE097
处的一阶至五阶导数 分别为：
Figure 425678DEST_PATH_IMAGE098
(4)
上式中
Figure 363678DEST_PATH_IMAGE099
表示虚数。则由高阶导数推导出的等效无量纲波数为：
Figure 503673DEST_PATH_IMAGE100
(5)
而根据傅里叶分析得到的五次样条格式的修正无量纲波数为：
Figure 877498DEST_PATH_IMAGE101
(6)
其中实部为格式的色散特性，虚部为格式的耗散特性，
Figure 29124DEST_PATH_IMAGE102
为人工粘性系数，格式的 耗散大小可通过人工粘性系数调节。
图2显示了由高阶导数求出的等效无量纲波数与五次样条格式修正无量纲波数对比图，其中纵坐标表示色散特性，由图中可以看出在修正无量纲波数的可识别范围内，等效无量纲波数与修正无量纲波数符合的很好，说明由波数识别方法得到的等效无量纲波数在一定的波数范围内可以用来代表格式的修正无量纲波数。
（3）自适应人工粘性判据
由于等效无量纲波数在一定波数范围内可以代表格式的修正无量纲波数，因此可 以根据图2进行自适应人工粘性判据的设定：根据修正的无量纲波数的特性设一阈值
Figure 528239DEST_PATH_IMAGE103
， 当
Figure 296475DEST_PATH_IMAGE104
时，流场等效无量纲波数在五次样条重构格式的精确可识别范围内，人工粘 性系数取零或者很小；当
Figure 539237DEST_PATH_IMAGE105
时，此时五次样条重构格式无法精确识别流场波数， 需要增大人工粘性系数来抑制误差。前者根据Lele et al提出的色散误差容限
Figure 810950DEST_PATH_IMAGE106
来确定，当
Figure 480965DEST_PATH_IMAGE107
时，
Figure 470918DEST_PATH_IMAGE108
。后者满足Hu & Adams提出的色散耗散关系
Figure 986213DEST_PATH_IMAGE109
，由此得出的自适应判据如下：
Figure 174749DEST_PATH_IMAGE110
(7)
图3为人工粘性系数随波数的变化曲线，可以看出随着波数的增大，人工粘性随之变化。图4为自适应人工粘性系数对应的耗散曲线，其中纵坐标表示格式的耗散特性，可以看出自适应耗散在低波数区域耗散尽量小，不引入耗散误差；在高波数区域尽量大以抑制色散误差带来的波动。
基于上述理论依据，具体方案如下：
如图1所示，一种适用于五次样条重构格式的自适应人工粘性控制方法，包括以下步骤：
步骤1、读取流场数据，求解五次样条重构方程组得到单元界面处物理量的状态值及其一阶导数；
步骤2、根据单元界面处物理量的状态值及一阶导数，计算出单元界面处的二阶至五阶导数；
步骤3、根据波数识别方法，由单元界面处流场密度的高阶导数求出单元界面处流场的等效无量纲波数，代入自适应人工粘性判据确定此处流场的人工粘性系数；
步骤4、根据单元界面处物理量的状态值计算无粘数值通量，并根据人工粘性系数添加六阶人工粘性项，最后采用相对应的时间离散方法在时间上进行推进。
具体的，
步骤（1），
在结构网格有限体积方法的基本框架下，不计质量力和源项，在直坐标系下欧拉方程为：
Figure 625453DEST_PATH_IMAGE001
(8)
其中，
Figure 227336DEST_PATH_IMAGE002
为守恒变量，
Figure 952846DEST_PATH_IMAGE003
、
Figure 323785DEST_PATH_IMAGE111
、
Figure 43977DEST_PATH_IMAGE005
为直角坐标系
Figure 867577DEST_PATH_IMAGE112
下的无粘通量，具体 表达式为：
Figure 131199DEST_PATH_IMAGE007
(9)
其中，
Figure 559906DEST_PATH_IMAGE008
、
Figure 742626DEST_PATH_IMAGE009
、
Figure 460046DEST_PATH_IMAGE010
分别为流场的密度、速度矢量以及流场压力，
Figure 386414DEST_PATH_IMAGE011
为总能，其具体表达式为：
Figure 341731DEST_PATH_IMAGE012
(10)
其中
Figure 960932DEST_PATH_IMAGE013
为气体常数；
将欧拉方程在三维结构网格的控制体单元
Figure 962386DEST_PATH_IMAGE014
上进行积分，得到：
Figure 567811DEST_PATH_IMAGE015
(11)
其中
Figure 502268DEST_PATH_IMAGE016
为控制体单元的体积，
Figure 433315DEST_PATH_IMAGE017
是控制体面上的外法线 方向，
Figure 718803DEST_PATH_IMAGE018
是守恒变量在控制体单元上的平均值：
Figure 862340DEST_PATH_IMAGE019
(12)
Figure 120146DEST_PATH_IMAGE020
为通量张量，
Figure 81149DEST_PATH_IMAGE021
即为无粘通量，
Figure 729299DEST_PATH_IMAGE022
为控制 体的六个面，以
Figure 801160DEST_PATH_IMAGE023
表示为控制体的面标号，则面积分项写为控制体各面积分的加和：
Figure 848226DEST_PATH_IMAGE024
(13)
每个面上
Figure 714551DEST_PATH_IMAGE025
的具体表达式为：
Figure 177893DEST_PATH_IMAGE026
(14)
其中
Figure 663232DEST_PATH_IMAGE027
分别为单元界面处的密度、三个方向上的速度、 压力以及总能的状态值；
在求解过程中根据单元的平均值重构出单元界面处的物理量的状态值，五次样条重构方程组为：
Figure 957947DEST_PATH_IMAGE113
(15)
其中，
Figure 870540DEST_PATH_IMAGE029
Figure 617916DEST_PATH_IMAGE030
(16)
Figure 906946DEST_PATH_IMAGE031
Figure 993850DEST_PATH_IMAGE032
下标
Figure 467557DEST_PATH_IMAGE033
分别代表着
Figure 577596DEST_PATH_IMAGE034
三个方向，
Figure 529371DEST_PATH_IMAGE035
和
Figure 408465DEST_PATH_IMAGE036
分别代表
Figure 787494DEST_PATH_IMAGE037
单元和
Figure 447146DEST_PATH_IMAGE038
单元从单元左侧界面中心至单元右侧界面中心的距离，
Figure 77978DEST_PATH_IMAGE039
分别指的
Figure 201792DEST_PATH_IMAGE040
五个守恒物理量，
Figure 892667DEST_PATH_IMAGE041
和
Figure 836353DEST_PATH_IMAGE042
分别指在
Figure 333193DEST_PATH_IMAGE043
界面处第
Figure 918371DEST_PATH_IMAGE044
个物理量的界面状态值及其一阶导数，
Figure 904781DEST_PATH_IMAGE045
和
Figure 476708DEST_PATH_IMAGE046
分别指
Figure 573977DEST_PATH_IMAGE047
单 元和
Figure 344487DEST_PATH_IMAGE048
单元处的第
Figure 377165DEST_PATH_IMAGE049
个物理量的单元平均值，通过求解重构方程，可得到单元界面
Figure 295442DEST_PATH_IMAGE050
处物理量的状态值
Figure 806189DEST_PATH_IMAGE051
及其一阶导数
Figure 493523DEST_PATH_IMAGE052
；上述五次样条重构方 程组、界面状态值及一阶导数均为
Figure 697102DEST_PATH_IMAGE053
方向，
Figure 102675DEST_PATH_IMAGE068
方向和
Figure 682693DEST_PATH_IMAGE055
方向采用相同的方法可进行求解。
步骤（2）根据单元界面
Figure 224532DEST_PATH_IMAGE056
处物理量的状态值
Figure 926909DEST_PATH_IMAGE057
及一阶导数
Figure 429566DEST_PATH_IMAGE042
计算单元中心
Figure 406749DEST_PATH_IMAGE058
处物理量的一阶至四阶导数：
Figure 209620DEST_PATH_IMAGE114
Figure 879636DEST_PATH_IMAGE115
Figure 601080DEST_PATH_IMAGE116
(17)
Figure 54058DEST_PATH_IMAGE117
从而得到单元界面
Figure 570490DEST_PATH_IMAGE060
处物理量的二阶至五阶导数：
Figure 286773DEST_PATH_IMAGE118
(18)
步骤（3）：式(18)得到单元界面处物理量的二阶至五阶导数后，由式(5)可得到流场界面处的等效无量纲波数，此处求解过程中使用密度的高阶导数。根据自适应人工粘性判据式(7)，确定界面处的人工粘性系数。
步骤（4）：根据控制方程求解无粘数值通量
将单元界面处物理量的状态值带入控制方程计算无粘数值通量，无粘通量计算过程为：
Figure 623076DEST_PATH_IMAGE075
(19)
其中，
Figure 83008DEST_PATH_IMAGE076
界面上无粘通量
Figure 453946DEST_PATH_IMAGE077
为：
Figure 606710DEST_PATH_IMAGE119
(20)
其中，
Figure 164730DEST_PATH_IMAGE079
为
Figure 428352DEST_PATH_IMAGE080
界面处的外法向速度；
Figure 653797DEST_PATH_IMAGE081
为六阶人工粘性项，其具体表达式为:
Figure 977462DEST_PATH_IMAGE083
(21)
其中，
Figure 553937DEST_PATH_IMAGE084
为自适应人工粘性系数；
Figure 355671DEST_PATH_IMAGE085
为尺度因子，对于各向 同性模型，其计算公式为：
Figure 638885DEST_PATH_IMAGE120
(22)
其中
Figure 258085DEST_PATH_IMAGE087
为界面处的声速，
Figure 928713DEST_PATH_IMAGE088
为界面
Figure 658772DEST_PATH_IMAGE089
处的面积矢 量。
最后求出通量以后，采用相对应的时间离散方法在时间上进行推进。
在本实施例中，给出以下自适应人工粘性的数值验证，如图5、6、7为采用固定人工 粘性的二维黎曼问题的密度等值线图和使用自适应人工粘性的二维黎曼问题密度等值线 图，由图中可以看出：图5
Figure 734175DEST_PATH_IMAGE121
时的耗散太小，无法抑制高频振荡,且计算流场不 对称；图7
Figure 524277DEST_PATH_IMAGE122
时的结果耗散太大，流场结构不丰富；图6采用自适应人工粘性 后的流场具有对称且更为丰富的结构。
本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。如果本领域技术人员，在不脱离本发明的精神所做的非实质性改变或改进，都应该属于本发明权利要求保护的范围。
本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。
本说明书中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

Claims (1)

1.一种适用于五次样条重构格式的自适应人工粘性控制方法，其特征在于，包括以下步骤：
步骤1、读取流场数据，求解五次样条重构方程组得到单元界面处物理量的状态值及其一阶导数；
步骤2、根据单元界面处物理量的状态值及一阶导数，计算出单元界面处的二阶至五阶导数；
步骤3、根据波数识别方法，由单元界面处流场密度的高阶导数求出单元界面处流场的等效无量纲波数，代入自适应人工粘性判据确定该处流场的人工粘性系数；
步骤4、根据单元界面处物理量的状态值计算无粘数值通量，并根据人工粘性系数添加六阶人工粘性项，最后采用相对应的时间离散方法在时间上进行推进；
所述步骤1中求解过程具体包括：
在结构网格有限体积方法的基本框架下，不计质量力和源项，在直坐标系下欧拉方程为：
Figure 271659DEST_PATH_IMAGE001
其中，
Figure 51397DEST_PATH_IMAGE002
为守恒变量，
Figure 446606DEST_PATH_IMAGE003
、
Figure 198661DEST_PATH_IMAGE004
、
Figure 21124DEST_PATH_IMAGE005
为直角坐标系
Figure 288157DEST_PATH_IMAGE006
下的无粘通量，具体表达式为：
Figure 487057DEST_PATH_IMAGE007
其中，
Figure 357535DEST_PATH_IMAGE008
、
Figure 350899DEST_PATH_IMAGE009
、
Figure 105228DEST_PATH_IMAGE010
分别为流场的密度、速度矢量以及流场压力，
Figure 842240DEST_PATH_IMAGE011
为总能，其具体表达式为：
Figure 365625DEST_PATH_IMAGE012
其中
Figure 733153DEST_PATH_IMAGE013
为气体常数；
将欧拉方程在三维结构网格的控制体单元
Figure 974778DEST_PATH_IMAGE014
上进行积分，得到：
Figure 249902DEST_PATH_IMAGE015
其中
Figure 893373DEST_PATH_IMAGE016
为控制体单元的体积，
Figure 415490DEST_PATH_IMAGE017
是控制体面上的外法线方向，
Figure 144411DEST_PATH_IMAGE018
是守恒变量在控制体单元上的平均值：
Figure 223226DEST_PATH_IMAGE019
Figure 721203DEST_PATH_IMAGE020
为通量张量，
Figure 492850DEST_PATH_IMAGE021
即为无粘通量，
Figure 381172DEST_PATH_IMAGE022
为控制体的六个面，以
Figure 263677DEST_PATH_IMAGE023
表示为控制体的面标号，则面积分项写为控制体各面积分的加和：
Figure 616161DEST_PATH_IMAGE024
每个面上
Figure 558709DEST_PATH_IMAGE025
的具体表达式为：
Figure 183595DEST_PATH_IMAGE026
其中
Figure 604212DEST_PATH_IMAGE027
分别为单元界面处的密度、三个方向上的速度、压力以及总能的状态值；
在求解过程中根据单元的平均值重构出单元界面处的物理量的状态值，五次样条重构方程组为：
Figure 76781DEST_PATH_IMAGE028
其中，
Figure 190231DEST_PATH_IMAGE029
Figure 115461DEST_PATH_IMAGE030
Figure 277452DEST_PATH_IMAGE031
Figure 338949DEST_PATH_IMAGE032
下标
Figure 623300DEST_PATH_IMAGE033
分别代表着
Figure 35827DEST_PATH_IMAGE034
三个方向，
Figure 250777DEST_PATH_IMAGE035
和
Figure 432359DEST_PATH_IMAGE036
分别代表
Figure 622032DEST_PATH_IMAGE037
单元和
Figure 521855DEST_PATH_IMAGE038
单元从单元左侧界面中心至单元右侧界面中心的距离，
Figure 87965DEST_PATH_IMAGE039
分别指的
Figure 327317DEST_PATH_IMAGE040
五个守恒物理量，
Figure 687891DEST_PATH_IMAGE041
和
Figure 75010DEST_PATH_IMAGE042
分别指在
Figure 444812DEST_PATH_IMAGE043
界面处第
Figure 525288DEST_PATH_IMAGE044
个物理量的界面状态值及其一阶导数，
Figure 322342DEST_PATH_IMAGE045
和
Figure 931178DEST_PATH_IMAGE046
分别指
Figure 104671DEST_PATH_IMAGE047
单元和
Figure 53035DEST_PATH_IMAGE048
单元处的第
Figure 20991DEST_PATH_IMAGE049
个物理量的单元平均值，通过求解重构方程，可得到单元界面
Figure 117123DEST_PATH_IMAGE050
处物理量的状态值
Figure 828727DEST_PATH_IMAGE051
及其一阶导数
Figure 693915DEST_PATH_IMAGE052
；上述五次样条重构方程组、界面状态值及一阶导数均为
Figure 19723DEST_PATH_IMAGE053
方向；
所述步骤2具体计算过程为：根据单元界面
Figure 603151DEST_PATH_IMAGE054
处物理量的状态值
Figure 852867DEST_PATH_IMAGE055
及一阶导数
Figure 838140DEST_PATH_IMAGE042
计算单元中心
Figure 85582DEST_PATH_IMAGE056
处物理量的一阶至四阶导数：
Figure 156306DEST_PATH_IMAGE057
从而得到单元界面
Figure 209713DEST_PATH_IMAGE058
处物理量的二阶至五阶导数：
Figure 49493DEST_PATH_IMAGE059
；
波数识别方法为：
Figure 530153DEST_PATH_IMAGE060
其中
Figure 9544DEST_PATH_IMAGE061
即为单元界面
Figure 866642DEST_PATH_IMAGE062
处的等效无量纲波数，
Figure 560929DEST_PATH_IMAGE063
为一个小量，防止分母为零，此处取
Figure 212490DEST_PATH_IMAGE064
；
自适应人工粘性判据为：
Figure 195489DEST_PATH_IMAGE065
其中，
Figure 590698DEST_PATH_IMAGE066
为即为自适应人工粘性系数；
所述步骤4具体包括：将单元界面处物理量的状态值带入控制方程计算无粘数值通量，无粘通量计算过程为：
Figure 139491DEST_PATH_IMAGE067
其中，
Figure 227533DEST_PATH_IMAGE068
界面上无粘通量
Figure 494567DEST_PATH_IMAGE069
为：
Figure 880417DEST_PATH_IMAGE070
其中，
Figure 283717DEST_PATH_IMAGE071
为
Figure 542660DEST_PATH_IMAGE072
界面处的外法向速度；
Figure 296989DEST_PATH_IMAGE073
为六阶人工粘性项，其具体表达式为:
Figure 971684DEST_PATH_IMAGE074
其中，
Figure 495069DEST_PATH_IMAGE075
为自适应人工粘性系数；
Figure 924914DEST_PATH_IMAGE076
为尺度因子，对于各向同性模型，其计算公式为：
Figure 900960DEST_PATH_IMAGE077
其中
Figure 631543DEST_PATH_IMAGE078
为界面处的声速，
Figure 275014DEST_PATH_IMAGE079
为界面
Figure 610180DEST_PATH_IMAGE080
处的面积矢量。
CN202110264089.2A 2021-03-11 2021-03-11 一种适用于五次样条重构格式的自适应人工粘性控制方法 Active CN112685978B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110264089.2A CN112685978B (zh) 2021-03-11 2021-03-11 一种适用于五次样条重构格式的自适应人工粘性控制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110264089.2A CN112685978B (zh) 2021-03-11 2021-03-11 一种适用于五次样条重构格式的自适应人工粘性控制方法

Publications (2)

Publication Number Publication Date
CN112685978A CN112685978A (zh) 2021-04-20
CN112685978B true CN112685978B (zh) 2021-06-08

Family

ID=75458390

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110264089.2A Active CN112685978B (zh) 2021-03-11 2021-03-11 一种适用于五次样条重构格式的自适应人工粘性控制方法

Country Status (1)

Country Link
CN (1) CN112685978B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114117966B (zh) * 2021-12-09 2023-04-11 中国空气动力研究与发展中心计算空气动力研究所 一种物理属性与数据驱动耦合的流声模态分解及预测方法

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102460113A (zh) * 2009-05-11 2012-05-16 西安大略大学 监测介质的颗粒尺寸分布的超声法方法
CN102954910A (zh) * 2011-08-26 2013-03-06 航天科工防御技术研究试验中心 防热涂层拉伸剪切强度测试模型及其制作方法
CN104298797A (zh) * 2013-07-16 2015-01-21 中国石油化工股份有限公司 一种确定缝洞型油藏高导流通道圈闭剩余油的方法
CN108197072A (zh) * 2017-12-27 2018-06-22 中国空气动力研究与发展中心计算空气动力研究所 一种基于加权守恒变量阶跃的高精度间断Galerkin人工粘性激波捕捉方法
CN109032077A (zh) * 2018-09-05 2018-12-18 沈阳建筑大学 一种基于刀具姿态控制的五轴数控加工指令点插补方法
CN109614577A (zh) * 2018-12-10 2019-04-12 山东大学苏州研究院 一种Burgers方程求解方法及装置
US20200109740A1 (en) * 2018-10-06 2020-04-09 Toyota Jidosha Kabushiki Kaisha Spline telescopic shaft of vehicular propeller shaft
WO2020165059A1 (en) * 2019-02-11 2020-08-20 Miltenyi Biotec B.V. & Co. KG Generation of human pluripotent stem cell derived artificial tissue structures without three dimensional matrices
CN112052632A (zh) * 2020-07-27 2020-12-08 空气动力学国家重点实验室 一种高超声速流向转捩预测方法
CN112214869A (zh) * 2020-09-03 2021-01-12 空气动力学国家重点实验室 一种求解欧拉方程的改进型高阶非线性空间离散方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104317997A (zh) * 2014-10-17 2015-01-28 北京航空航天大学 一种高负荷风扇/压气机端壁造型优化设计方法
CN110209048A (zh) * 2019-05-20 2019-09-06 华南理工大学 基于动力学模型的机器人时间最优轨迹规划方法、设备

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102460113A (zh) * 2009-05-11 2012-05-16 西安大略大学 监测介质的颗粒尺寸分布的超声法方法
CN102954910A (zh) * 2011-08-26 2013-03-06 航天科工防御技术研究试验中心 防热涂层拉伸剪切强度测试模型及其制作方法
CN104298797A (zh) * 2013-07-16 2015-01-21 中国石油化工股份有限公司 一种确定缝洞型油藏高导流通道圈闭剩余油的方法
CN108197072A (zh) * 2017-12-27 2018-06-22 中国空气动力研究与发展中心计算空气动力研究所 一种基于加权守恒变量阶跃的高精度间断Galerkin人工粘性激波捕捉方法
CN109032077A (zh) * 2018-09-05 2018-12-18 沈阳建筑大学 一种基于刀具姿态控制的五轴数控加工指令点插补方法
US20200109740A1 (en) * 2018-10-06 2020-04-09 Toyota Jidosha Kabushiki Kaisha Spline telescopic shaft of vehicular propeller shaft
CN109614577A (zh) * 2018-12-10 2019-04-12 山东大学苏州研究院 一种Burgers方程求解方法及装置
WO2020165059A1 (en) * 2019-02-11 2020-08-20 Miltenyi Biotec B.V. & Co. KG Generation of human pluripotent stem cell derived artificial tissue structures without three dimensional matrices
CN112052632A (zh) * 2020-07-27 2020-12-08 空气动力学国家重点实验室 一种高超声速流向转捩预测方法
CN112214869A (zh) * 2020-09-03 2021-01-12 空气动力学国家重点实验室 一种求解欧拉方程的改进型高阶非线性空间离散方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Investigation of a Low-Toxicity Energetic Binder for a Solid Propellant: Curing, Microstructures, and Performance;Song Ma等;《ACS Omega》;20201201;第5卷(第47期);30538–30548 *
三维边界层流动失稳与Bypass转捩模式研究;徐国亮;《中国博士学位论文全文数据库础科学辑》;20121215(第12期);A004-6 *

Also Published As

Publication number Publication date
CN112685978A (zh) 2021-04-20

Similar Documents

Publication Publication Date Title
CN107272403A (zh) 一种基于改进粒子群算法的pid控制器参数整定算法
CN110851929B (zh) 基于自适应网格的二维叶型优化设计方法及装置
CN111079228B (zh) 一种基于流场预测的气动外形优化方法
CN104778327B (zh) 基于人工神经网络的飞机翼型优化设计方法
CN106126791B (zh) 一种考虑几何不确定性的高超声速机翼气动力/热分析方法
CN102540882B (zh) 一种基于最小参数学习法的飞行器航迹倾角控制方法
CN112685978B (zh) 一种适用于五次样条重构格式的自适应人工粘性控制方法
CN111006843B (zh) 一种暂冲式超声速风洞的连续变速压方法
CN113850008B (zh) 飞行器气动特性预测的自适应网格扰动域更新加速方法
CN112214869B (zh) 一种求解欧拉方程的改进型高阶非线性空间离散方法
CN107640183B (zh) 一种基于迭代学习控制的列车运行控制方法
CN113127976B (zh) 一种边界层分离诱导转捩预测方法、装置、设备及介质
CN113569360B (zh) 一种风力机叶片抗颤振翼型簇设计方法
CN104834772A (zh) 基于人工神经网络的飞机翼型/机翼反设计方法
CN112001109A (zh) 再生核粒子算法实现结构冲击动力学仿真方法
CN110765706A (zh) 基于ohngbm(1,1)的翼型非定常失速气动系数建模方法
CN106446419A (zh) 一种火力发电厂燃煤锅炉的建模方法及系统
CN107329131A (zh) 一种利用粒子滤波的雷达微弱目标检测跟踪方法
CN116373846A (zh) 一种基于bp神经网络优化的后轮转向车辆稳定性控制方法
CN110852888A (zh) 一种基于粒子滤波的证券投资组合优化方法
CN114185276B (zh) 一种基于多维泰勒网的非线性严格系统输出反馈控制方法
CN113977571B (zh) 一种柔性关节机器人输出力矩控制方法
CN113158339B (zh) 一种针对sst湍流模型的湍流长度尺度修正方法
Zhuo et al. Parameter Identification of Tire Model Based on Improved Particle Swarm Optimization Algorithm
CN111597741B (zh) 一种基于改进Hicks-henne算法的仿生蟹滑翔姿态下翼型优化设计方法

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