中国科技核心期刊      中国指挥与控制学会会刊     军事装备类重点期刊
指挥控制与仿真

指挥控制与仿真 >

2022 , Vol. 44 >Issue 1: 58 - 62

DOI: https://doi.org/10.3969/j.issn.1673-3819.2022.01.008

理论研究

基于Richardson外推法的CFD离散误差分析

  • 左玲玉 1 ,
  • 司海青 1 ,
  • 李耀 1 ,
  • 仇静轩 1 ,
  • 李根 1 ,
  • 吴晓军 2 ,
  • 赵炜 2
展开
  • 1.南京航空航天大学通用航空与飞行学院, 江苏 南京 211100
  • 2.中国空气动力研究与发展中心, 四川 绵阳 621000

左玲玉(1997—),女,硕士研究生,研究方向为计算流体力学。

司海青(1976—),男,博士,教授。

Copy editor: 张培培

收稿日期: 2021-10-15

  要求修回日期: 2021-11-09

  网络出版日期: 2022-05-19

基金资助

南京航空航天大学研究生创新基地(实验室)开放基金(kfjj20200702)

江苏省研究生科研与实践创新计划(KYCX20_0214)

版权

版权所有,未经授权,不得转载、摘编本刊文章,不得使用本刊的版式设计。
收起

Analysis of CFD Discrete Error Based on Richardson Extrapolation Method

  • ZUO Ling-yu 1 ,
  • SI Hai-qing 1 ,
  • LI Yao 1 ,
  • QIU Jing-xuan 1 ,
  • LI Gen 1 ,
  • WU Xiao-jun 2 ,
  • ZHAO Wei 2
Expand
  • 1. College of General Aviation and Flight,Nanjing University of Aeronautics and Astronautics,Nanjing 211100
  • 2. China Aerodynamics Research and Development Center, Mianyang 621000, China

Received date: 2021-10-15

  Request revised date: 2021-11-09

  Online published: 2022-05-19

Copyright

Copyright reserved © 2022
Fold

摘要

为研究数值模拟中网格类型及网格尺度变化对离散误差产生的影响,对NACA2412三维翼型在0°、4°、8°、12°攻角下进行了计算流体力学(CFD)仿真,得到了对应条件下的升力系数及阻力系数,利用Richardson外推法研究多套连续变化的结构网格及非结构网格数值解的空间离散误差,并结合网格收敛指标(GCI)法对误差进行分析,最后,通过网格收敛性分析得到对应网格研究组修正后的精确解。结果表明:各个攻角条件下结构网格模型的离散误差均明显小于非结构网格,网格类型的改变对升力系数的影响更大,且细网格模型具有更好的网格收敛性。

关键词： 计算流体力学; NACA2412; 离散误差; Richardson 外推法; 网格收敛指标

本文引用格式

左玲玉 , 司海青 , 李耀 , 仇静轩 , 李根 , 吴晓军 , 赵炜 . 基于Richardson外推法的CFD离散误差分析[J]. 指挥控制与仿真, 2022 , 44(1) : 58 -62 . DOI: 10.3969/j.issn.1673-3819.2022.01.008

Abstract

In order to study the influence of grid type and grid scale changes on the discrete error in the numerical simulation, CFD simulation is performed on the NACA2412 three-dimensional airfoil at 0°, 4°, 8° and 12° angles of attack, and the lift coefficient and drag coefficient under corresponding conditions are obtained. The Richardson interpolation method is used to study the spatial discretization errors of multiple sets of continuously changing structured grids and unstructured grids. The grid convergence index (GCI) method is used to analyze the error. Finally, the modified exact solution of the corresponding grid research group is obtained by grid convergence analysis. The results show that the discrete errors of the structured grid model are significantly less than those of the unstructured grid under the conditions of various angles of attack. The change of the mesh type has a greater impact on the lift coefficient, and the fine mesh model has better mesh convergence.

Key words： computational fluid dynamics; NACA2412; discrete error; Richardson extrapolation method; grid convergence index

近年来,随着流体力学理论、数值方法以及高性能计算机技术的飞速发展,气动数值模拟技术(CFD)以其花费低、耗时短、损耗小等优点在航空航天领域发挥着越来越大的作用。然而,CFD技术要获得更广泛的应用及效益,对其可信度的要求也越来越高,所以,对CFD计算结果开展可信度研究十分必要,而可信度研究的基本内容和方法就是验证与确认(V&V)[1]。
验证主要是对数值模拟过程中的误差进行分析,以验证仿真结果的准确性;确认主要是将模拟结果与实验数据进行比较,以评估数值模型与实际模型之间的相似性[2]。国内外对CFD的验证与确认工作一直高度重视,迄今为止,大多数验证研究都是基于Richardson外推法[3]获得插值解从而得到数值计算的误差估计;Reed等[4]通过对高超声速飞行器二维和三维几何模型的机理识别以及CFD网格收敛性研究,最大程度地降低高超声速飞行器设计中各方面产生的不确定性;Lance等[5]对混合对流进行了实验研究,比较了实验数据和模型输出,为CFD模型提供了完整的验证数据;Islam等[6]对四种不同的船型进行验证与确认研究,使用基于安全系数和修正系数的方法估计算例的不确定性,虽然使用了相似的网格拓扑,但所有四个仿真计算都得到了不同的结果;于沛从算例的角度对验证确认的层次结构、算例的分类、选择等方面进行了分析说明[7];陈坚强等[8]利用Richardson插值法对典型高超声速流动问题进行了网格收敛性研究,然后,利用不确定度分析方法开展了该问题的确认;康顺[9]等总结了CFD模拟中的误差及估计准则,并以二维翼型为例介绍了如何进行误差估计和结果修正。通过总结分析各位专家学者们的研究发现,在验证与确认过程中针对数值模拟误差的研究占据着相当重要的位置,并且网格收敛性研究是必不可少的一部分。

1 CFD数值模拟误差

在采用CFD软件进行数值模拟时,模拟结果与真实结果之间存在一定误差,主要包括以下几个方面[10]。1)物理建模误差:由于模型简化、控制方程、边界条件等引入的误差。2)离散误差:主要源于对数学模型和定义域的离散,包括离散格式的截断误差、边界条件设置不合理、网格质量较差等情况导致的离散误差增大。3)迭代误差:CFD软件数值模拟过程中采用松弛迭代法求解各类代数方程时随时可能产生的误差,包括不完全迭代或迭代的步数太少等。4)几何外形:包括数模与真实模型外形之间的差异和生成网格时几何修改或近似引入的误差。5)来流条件变化:CFD仿真模拟中由于来流参数变化所导致的不确定性会使计算结果产生一定的偏差。
在各种误差中,空间离散误差是其中最重要且最难估计的部分,在误差源中占较大比重,所以,离散误差估计是我们的研究重点,而网格收敛性研究是目前空间离散误差估计最为普遍的方法[11]。

2 Richardson外推法

本文采用Richardson外推法对离散误差进行分析,该方法利用离散步长h将离散解与精确解的关系用泰勒展开的级数形式表示:
f k = f exact + g p h k p + o ( h k p + 1 )
其中,fk指第k网格上的值,gp指不随网格变化的第p阶误差项系数, f ^ exact是数值基准解,根据细化网格的数值解代入公式可以得到对应的基准解:
f ^ exact = f 1 - ε 21 r p - 1
其中,收敛精度p由下式确定:
p = ln ( ε 32 / ε 21 ) ln ( r )
同时可以得到网格上的误差估计:
g p = f 2 - f 1 r p - 1
而网格细化比r由下式定义:
r = N 1 N 2 1 / d
式中:Nk是网格k的网格单元数,d为空间维数。由于求解过程中会忽略高阶项,所以,误差往往偏大,对式(2)添加一个修正因子K[12]来提高对误差估计的精度,具体为
K = r p - 1 r p e - 1
其中,pe是对高阶项的估计。所以,修正后可以得到精确解Fexact:
F exact = f 1 - r p - 1 r p e - 1 ε 21 r p - 1
Roache将网格数值解的网格收敛指标(GCI)[13]定义为
GCI = F s r p - 1 | f 2 - f 1 |
式中,Fs为安全系数,对于三套网格,取Fs=1.25,最后定义收敛率R为
R = ε 21 ε 32
当0<R<1时,属于单调收敛,当R<0时为振荡收敛,而当R>1时为发散,Richardson插值法进行离散误差计算时只对单调收敛问题有效。

3 CFD离散误差分析

3.1 NACA2412翼型CFD仿真计算

本文对NACA2412三维翼型进行CFD数值仿真,其0°、4°、8°、12°攻角下升力系数和阻力系数的具体实验数据[14]如表1所示。
表1 升力系数和阻力系数不同攻角下的实验值
气动参数 0° 4° 8° 12°
升力系数 0.114 9 0.359 5 0.628 2 0.877 9
阻力系数 0.022 8 0.035 4 0.060 7 0.099 0
对模型进行了6种不同数量的网格划分,其中,前三种采用结构网格,后三种采用非结构网格,图1展示了采用O型网格划分情况,图2为非结构网格的划分情况,两种类型网格均在机翼前缘和后缘进行网格加密以便更好地捕捉流动特性。
图1 NACA2412结构网格划分细节图
图2 NACA2412非结构网格划分细节图
对于细网格U1、中网格U2、粗网格U3,相邻两套网格的细化比为 2 ,其中的U1网格单元数为248.20万,网格节点数为241.75万。对于非结构网格对应的细网格U4、中网格U5、粗网格U6,选择相同的网格细化比,其中,U4的网格单元数为1 598.95万,网格节点数为310.54万。
本文采用Fluent软件进行定常计算,湍流方程采用Spalart-Allmaras模型,马赫数Ma=0.23,动量方程及湍流耗散率方程均采用二阶迎风离散格式,计算残差设置为10-5,经过约1200步迭代计算后各个计算模型均达到收敛条件,表2和表3对应于NACA2412翼型不同网格数量在0°、4°、8°、12°攻角时CFD计算得到的升力系数和阻力系数。将不同攻角条件下的CFD计算值与实验数据进行对比,发现随着网格尺度的变化计算结果单调收敛且处于渐进收敛域内,可以进行网格收敛性研究。
表2 不同攻角下升力系数的CFD计算值
网格 0° 4° 8° 12°
U1 0.116 3 0.361 9 0.629 9 0.879 8
U2 0.119 2 0.364 3 0.636 9 0.882 5
U3 0.122 8 0.369 3 0.648 1 0.885 3
U4 0.116 0 0.362 6 0.632 1 0.879 4
U5 0.121 3 0.366 7 0.645 0 0.884 1
U6 0.127 8 0.371 6 0.664 9 0.885 6
表3 不同攻角下阻力系数的CFD计算值
网格 0° 4° 8° 12°
U1 0.022 0 0.035 0 0.060 1 0.099 0
U2 0.019 2 0.034 1 0.058 4 0.095 3
U3 0.015 7 0.032 9 0.055 1 0.096 0
U4 0.021 2 0.034 9 0.059 7 0.098 9
U5 0.017 4 0.033 3 0.057 0 0.096 5
U6 0.012 9 0.030 8 0.053 5 0.096 8

3.2 Richardson 外推法结果分析

为了展示网格类型以及网格粗密程度对误差估计精度的影响,将6套网格分为两组进行研究,其中,网格U1、U2和U3为网格研究组1,网格U4、U5和U6为网格研究组2,表4为利用公式(9)计算得到的两组网格不同攻角条件下的网格收敛判断。针对升力系数的计算,在0°、4°及8°时,网格研究组1(U1U3)和网格研究组2(U4U6)均为单调收敛,可以采用Richardson 插值法来估计其升力系数的计算误差,但网格研究组2(U4U6)在12°攻角时发散,不能进行误差估计。对于各个攻角的阻力系数的收敛性计算,两个网格研究组合在小于最大升力攻角12°时均为单调收敛,可进行误差估计,在12°攻角时都属于振荡收敛,无法进行阻力系数的误差估计。
表4 升阻力系数的收敛性判断
攻角 网格 升力系数 阻力系数
R 收敛性 R 收敛性
0° 1 <1 单调收敛 <1 单调收敛
2 <1 单调收敛 <1 单调收敛
4° 1 <1 单调收敛 <1 单调收敛
2 <1 单调收敛 <1 单调收敛
8° 1 <1 单调收敛 <1 单调收敛
2 <1 单调收敛 <1 单调收敛
12° 1 <1 单调收敛 <0 振荡收敛
2 >1 发散 <0 振荡收敛
针对0°、4°、8°攻角条件下升阻力系数,计算得到网格研究组1、2对应的离散误差估计gp以及网格收敛指标GCI值(下标l1代表研究组1升力系数对应的计算结果,下标d1代表研究组1的阻力系数对应的计算结果,依此类推),具体见表5。
表5 不同攻角下网格组合升阻力系数的分析结果
气动参数 内容 0° 4° 8°
升力系数 gpl1 0.012 5 0.002 5 0.012 0
gpl2 0.023 4 0.012 0 0.023 8
GCIl1 0.015 6 0.003 1 0.015 0
GCIl2 0.029 3 0.026 7 0.02 9
阻力系数 gpd1 0.012 7 0.001 6 0.000 9
gpd2 0.011 6 0.006 5 0.014 1
GCId1 0.015 9 0.002 0 0.002 2
GCId2 0.014 5 0.008 2 0.017 6
结果表明,两个网格研究组的网格误差和网格收敛指标计算结果接近,但结构网格的网格误差以及网格收敛指标整体上均小于非结构网格,同时网格类型的改变对升力系数的影响更大。
对0°、4°和8°攻角下升阻力系数计算值求解得到收敛精度p如表6所示。
表6 不同攻角下网格组合升阻力系数的收敛精度
网格 气动参数 0° 4° 8°
研究组1 升力系数 0.601 9 1.995 0 1.333 6
阻力系数 0.580 2 1.275 4 1.922 6
研究组2 升力系数 0.589 2 0.505 0 1.251 3
阻力系数 0.502 8 1.218 9 0.706 7
根据计算结果发现,同一攻角条件下升阻力系数收敛精度相差不大,说明CFD计算结果较为准确,同时网格研究组1的收敛精度大于网格研究组2的收敛精度,可见结构网格的离散解能以更快的速率接近精确解。所以,求解修正系数K时,对于网格研究组1根据经验pe 取2,而对于网格研究组2,假定结构网格组合达到对应的网格无关解,所以其各个攻角下的pe取其对应的网格研究组1计算得到的收敛精度p,由不同攻角下各网格组合升阻力系数收敛精度和对应的修正系数进一步求解得到修正后的精确解。
图3和图4给出了0°、4°、8°攻角下不同网格研究组合随网格数量变化的CFD计算值、修正后的精确解及其实验数据的对比图。其中,黑色实线为实验值,蓝色点划线为6套网格的CFD计算值,红色双点划线和绿色划线对应于网格研究组修正后的精确解F1、F2,这两条点划线所覆盖的网格长度对应于各自进行误差估计时所采用的网格研究组。
图3 不同攻角下升力系数对比图
图4 不同攻角下阻力系数对比
根据图3中的结果可见,在0°、4°、8°这三种攻角条件下网格研究组1计算的升力系数的精确解与实验值的误差均在2%以内,说明该种类型网格进行计算时随着网格数量的增加对机翼尾迹以及边界层处的流动情况模拟效果较为精确,结果趋于网格无关值,网格尺度对计算结果的影响逐渐消失。而网格研究组2由于CFD计算值随网格单元数的变化有较大差异,所以计算所得的精确解与实验值的误差远大于网格研究组1。
由图4可知,网格研究组1对于阻力系数的计算值也随着网格数量的增加逐渐达到网格无关解,但是本文主要考虑修正网格尺度变化所引起的误差,忽略了可能涉及的不完全迭代误差、时间步长误差等情况,所以所得的精确解与实验值仍有差异。对于网格研究组2,在4°攻角时CFD计算所得的阻力系数随着网格数量的增加更加趋近于网格无关解,而0°和8°攻角下得到的精确解与实验值的误差较大,说明对于网格研究组2,目前的网格单元数仍未达到最优效果,需要对网格尺寸进行调整,通过更加精细的网格来获得更加精确的计算结果。除此之外,可能还需要对模型本身或计算条件、边界条件等进行进一步的改进。

4 结束语

本文为评估网格类型及网格数量对NACA2412机翼CFD离散误差的影响,对其三维模型分别进行结构、非结构两种类型共6套网格的划分,并进行了三维定常流场数值模拟,结合Richardson外推法和GCI网格收敛指标法对升力系数、阻力系数进行了离散误差估计,对比分析了CFD计算结果、修正后的精确值、实验数据,完成了网格收敛性研究并得到结论:结构网格与非结构网格模型的离散误差有较大差异,网格类型的改变对升力系数的影响更大,两种方法计算所得的非结构网格升力系数的误差为结构网格误差的两倍左右;随着网格单元数的增加,各个攻角条件下结构网格的升阻力系数均比非结构网格更接近于网格无关值。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
邓小刚, 宗文刚, 张来平, 等. 计算流体力学中的验证与确认[J]. 力学进展, 2007, 37(2):279-288.

[2]
Oberkampf W L. Verification, Validation, and Predictive Capability in Computational Engineering and Physics[J]. Applied Mechanics Reviews, 2004, 57(5):345-384.

DOI

[3]
Roy C J. Grid Convergence Error Analysis for Mixed-Order Numerical Schemes[J]. AIAA Journal, 2003, 41(4):595-604.

DOI

[4]
Reed H L, Perez E, Kuehl J. Verification and Validation Issues in Hypersonic Stability and Transition Prediction[J]. Journal of Spacecraft and Rockets, 2015, 52(1):29-37.

DOI

[5]
Lance B W, Harris J R, Smith B L. Experimental Validation Benchmark Data for Computational Fluid Dynamics of Mixed Convection on a Vertical Flat Plate[J]. Journal of Verification Validation & Uncertainty Quantification, 2016, 1(2):021005.

[6]
Islam H, Soares C G. Uncertainty Analysis in Ship Resistance Prediction Using OpenFOAM[J]. Ocean Engineering, 2019, 191(2):105805.

DOI

[7]
于沛. CFD软件的算例验证确认技术[C]// 计算流体力学研究进展,第十二届全国计算流体力学会议论文集,北京, 2004.

[8]
陈坚强, 张益荣. 基于Richardson插值法的CFD验证和确认方法的研究[J]. 空气动力学学报, 2012, 30(2):176-183.

[9]
康顺, 石磊, 戴丽萍, 范忠瑶. CFD模拟的误差分析及网格收敛性研究[J]. 工程热物理学报, 2010, 31(12):2009-2013.

[10]
赵训友, 林景松, 童晓艳. 基于Richardson外推法的CFD中离散不确定度估计[J]. 系统仿真学报, 2014, 26(10):2315-2320.

[11]
陈坚强, 张益荣, 郭勇颜. 高超声速流动数值模拟方法及应用[M]. 北京: 北京科学出版社, 2019.

[12]
Wilson R V, Stern F, Coleman H W, et al. Comprehensive Approach to Verification and Validation of CFD Simulations-Part 2: Application for Rans Simulation of a Cargo/Container Ship[J]. J Fluid Eng, 2001, 123(4):803-810.

DOI

[13]
Roache P J. Perspective: A Method for Uniform Reporting of Grid Refinement Studies[J]. Journal of Fluids Engineering, 1994, 116(3):405-413.

DOI

[14]
Prabhakar A, Ohri A. Application of CFD Simulation in the Design of a Parabolic Winglet on NACA 2412[J]. Lecture Notes in Engineering & Computer Science, 2014, 2212(1):10-18.

Options
文章导航

/

〈 〉