文章快速检索     高级检索
  空气动力学学报  2019, Vol. 37 Issue (3): 373-381  DOI: 10.7638/kqdlxxb-2017.0212

引用本文  

刘君, 陈泽栋, 陈洁, 等. 一种用于间断装配的流场结构辨识算法[J]. 空气动力学学报, 2019, 37(3): 373-381.
LIU J, CHEN Z D, CHEN J, et al. A surface extraction method of flow feature for shock-fitting scheme[J]. Acta Aerodynamica Sinica, 2019, 37(3): 373-381.

基金项目

国家自然科学基金(11872144, 11532016)

作者简介

刘君(1965-), 男, 博士, 教授, 研究方向:计算流体力学, 空气动力学.E-mail:liujun65@dlut.edu.cn

文章历史

收稿日期:2017-12-11
修订日期:2018-04-11
一种用于间断装配的流场结构辨识算法
刘君1 , 陈泽栋2 , 陈洁1 , 邹东阳1     
1. 大连理工大学 航空航天学院, 大连 116024;
2. 大连理工大学 工程力学系, 大连 116024
摘要:使用激波装配法时,初始激波是否准确将会对计算过程产生影响。为了确定初始激波的位置,提出了一种新的流场结构辨识算法。该算法以捕捉法计算得到的流场作为系统观测数据,根据密度、压力等参数从该数据中获取激波和接触间断等流动特征周围的网格节点作为离散点集。通过将该离散点集分割成若干子区域,在各子区域内进行分片拟合,最终将离散点集拟合成连续光滑的实体模型,并将此作为初始激波面。在二维方法的基础上,通过引入单位球模型成功将该辨识算法拓展到三维应用。结果表明,采用该方法获得的间断曲面(激波和接触间断)与捕捉法流场中的间断分布吻合较好,作为初始间断面用于装配法可快速得到收敛解。该方法解决了应用激波装配法时确定初始间断面的难题。此外,该方法还可用于网格自适应方法。选择不同流动参数,可以获得相应流场特征结构的空间曲面,在此曲面的基础上可进行网格局部加密或重剖分。该流场结构辨识算法用于网格自适应具有网格尺度自由设置的优势。
关键词激波探测    激波装配法    非结构网格自适应    曲面拟合    结构辨识    数值模拟    
A surface extraction method of flow feature for shock-fitting scheme
LIU Jun1 , CHEN Zedong2 , CHEN Jie1 , ZOU Dongyang1     
1. School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China;
2. Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, China
Abstract: A new surface extraction method named source-rays method is proposed to locate and shape initial shock wave surface for shock-fitting algorithms. In this method, the flow field data acquired by the shock-capturing schemes are used as the observed data, and the grid nodes with the features such as shock waves and contact discontinuities are obtained as a discrete point set according to the parameters of density, pressure and others. The discrete point set is finally assembled into a continuous smooth surface by being segmented into several subsets for fragment fitting. Based on the two-dimensional method, this identification algorithm is extended to three-dimensional application by introducing a unit sphere model. The discrete points are fitted to a series of triangular patches with determined connection relationship by a unit sphere, and a smooth continuous surface can be assembled. The results show that the feature surface (shock wave and contact discontinuity) obtained by the source-rays method coincides with the discontinuous distribution in the flow field calculated by the shock-capturing schemes. By using the extracted feature surface as the initial shock wave surface, the convergence solution can be quickly obtained in shock-fitting algorithm. In addition, this method can also be used in the adaptive mesh refinement procedure. According to different flow parameters, the spatial surfaces of the corresponding flow field structures can be obtained. Based on these surfaces, the mesh grid can be locally refined. This flow field feature surface extraction method has the advantage of grid scale set arbitrarily for mesh adaptation procedure.
Keywords: shock detection    shock-fitting    unstructured mesh adaptation    surface fitting    feature identification    numerical simulation    
0 引言

流体微团加速到超声速以后无法感知到达之前遇到的空间物体,只能在物体前通过急剧压缩的激波才能把速度降下来。因此,激波是可压缩流动的基本特征。早期计算流体力学(CFD)的算法无法模拟这一现象,直到1954年Lax提出处理包含有间断流动问题的弱解理论后[1-2],国内外众多科学家开始研究在计算过程中能自动分辨出激波的算法,提出TVD、NND、ENO、WENO等具有里程碑意义的激波捕捉格式[3],解决了航天航空等领域很多工程模拟问题,极大提高了CFD地位,使其发展为一个重要的学科分支。尽管建立这些格式时考虑了激波间断,但是本质上还是1950年von Neumann提出的把间断描述成流动参数连续变化过渡区的思路[4],构造出来的格式存在以下难以克服的理论困境。

数值计算的离散格式是对偏微分方程的近似,按照CFD收敛性理论,数值解随着网格尺度减小趋向理论解,CFD追求的极限是方程的理论解,原则上说对于方程不能描述的流动现象,数值解和实验值不具有可比性。海平面环境下马赫数小于10的激波厚度大约1×10-8 m,流体在这么小的空间内受到剧烈压缩,暂且不说准确模拟所需要的网格尺度及其带来的计算量,从连续介质假设出发建立的N-S方程是否适用于激波模拟也存在争议,因为Stoke假设(3λ+2μ=0)不再成立[5],计算时需要对体积膨胀黏性系数λ进行修正才能得到与实验相符的激波内部温度[6]。从事稀薄气体动力学研究的学者认为激波内部结构涉及到热力学非平衡,连续介质假设本身就不成立,“从理论上说,Navier-Stokes方程不能用来描述激波结构的[7]”。无黏流体力学把激波表示成间断,两侧参数和运动速度满足兰金-许贡纽(Rankine-Hugoniot, R-H)关系式,这是迄今为止几乎很少有争议的理论。所谓间断在数学上就是在同一空间位置对应至少两组(激波前后)流体变量,在多维情况下可能更多,例如二维马赫反射三岔点有四组流体变量。很多基于弱解理论构造的激波捕捉格式考虑了间断,但是用参数连续变化的思路描述激波注定不可能计算出真正的数学间断,捕捉到的激波主要是现象模拟。

按照上述观点,常规N-S方程及其理论解不能用来描述激波内部物理现象,那么基于离散格式的数值解至少在过渡区内没有物理意义,采用Euler方程产生的后果是无法评价过渡区内数值解的精度(激波前后参数相差很大),激波捕捉格式对于激波相交等现象更是无能为力。

20世纪60年代Moretti等人将激波装配方法和Lax-Wendroff格式结合,在CDC 600型计算机上进行,总耗时为360 s,给出了超声速钝头体问题的数值解,2002年Moretti回顾从业36年来CFD领域模拟激波的进展,针对超声速钝头体绕流比较了捕捉法和装配法计算结果,认为准确模拟激波还需要发展装配法[8-9]。对于50年前提出的装配法,文献[10]的评价具有代表性——“激波装配法具有计算精度高、物理概念清晰等优点,但是计算过程复杂、计算量大,并且必须事先知道激波的大概位置。一般来说,激波装配法仅适用于那些激波比较简单而事先可以估计激波形状和位置的定常流动问题。”

造成计算复杂的主要原因是求解过程中激波运动和变形。装配法减少了网格量,造成计算量大的原因是激波及其相交点等空间网格点需要判别、标记或特殊处理,增加代码编写复杂性和并行计算难度。随着CFD网格技术发展,这些难点也得到缓解。国外文献把最早Moretti提出的称为边界激波装配法,仅能处理形状较为规则和位移不大的脱体激波。后来Richtmyer等人[11]提出了浮动激波装配方法,在静止网格基础上嵌入R-H关系式,允许激波在背景网格上滑动,解决了内部激波装配难题。近年来Zhong等人将这种方法和高精度有限差分格式相结合有效提高内部光滑流场的计算精度,成功应用于许多高超声速流动的求解中[12-13]。基于结构网格的边界激波装配法和浮动激波装配法在模拟激波相交和马赫反射等复杂波系结构时面临一定困难。Paciorri和Bonfigolioli等人[14-16]采用非结构网格技术和局部网格重构提出了激波阵面追踪法,记录激波在背景网格上的位置,对激波附近的背景网格进行挖洞,然后再利用Delanay方法在洞内重构以激波阵面分界的新网格,时间推进过程中根据R-H关系式得到新的激波位置,重新判别、挖洞和填补背景网格。这种算法不需要考虑网格节点之间的连接关系,也就减少了对于特殊点进行处理的麻烦,增加了代码的普适性,但是在激波边界节点对应的网格单元只能采用一阶插值,像常用的多块重叠网格技术一样存在插值引起的精度降低问题。其次,不允许激波在一个时间步内跨越多个网格节点,通常选择足够小的时间步长,因而对计算效率有所影响。上述这3类装配法时间推进过程中需要判别激波位置,增加了计算量,而且激波前后网格数目是变化的,不利于并行计算。2015年刘君等人[17]提出了基于非结构动网格的激波(间断)装配方法,从能够描述网格变形的任意拉格朗日-欧拉坐标系(Arbitrary Lagrangian-Eulerian, ALE)下的控制方程出发,得到激波位置及其两侧的参数后,装配的激波就变成了一个移动边界,对于存在内部激波的计算域则变成类似于拼接网格的多个分区,激波运动由上下游参数驱动,采用成熟的动网格技术[18]进行移动和变形。与Paciorri及Bonfigolioli的装配法相比,在保留非结构网格的优点外,新的装配法不需要在每一步对计算网格进行挖洞重构,避免了插值误差,易于并行计算,并且时间步长不受限制。采用这种算法模拟了内部存在正规反射和马赫反射等复杂激波结构的算例[19-21],从效果看,较好地解决了传统装配法遇到的“计算过程复杂、计算量大”问题。

为了突破“激波装配法仅适用于那些激波比较简单而事先可以估计激波形状和位置的定常流动问题”的限制,刘君等人提出了在采用捕捉法计算的流场基础上辨识并装配主要激波结构、计算过程中自动装配内部细致结构和新生激波的解决思路。在技术途径上提出了新的数据结构,在激波边界通量计算时,如果采用Roe、van Leer、AUSM等格式就是捕捉法,如果采用R-H关系式就是装配法,实现两种方法之间任意转换。按照上述思路,为了解决装配法的应用难题,亟需发展激波、接触间断等流场特征结构的辨识算法。

目前,常用的激波辨识算法主要有三种[22]:第一种根据最大变量梯度确定激波面,该方法实现简单,但需要设置合适的过滤器才能消除错误的结果;第二种为法向马赫数法[24],该方法认为在激波面上流动的法向马赫数等于1,该方法与第一种方法相比适用性更强,但也需要设置适当的过滤器才能获得良好结果,且难以用于探测接触间断;第三种为基于特征线理论的方法[25-26],该方法与前两种方法相比更精确、更有效,但实施过程复杂且计算量大。

综上所述,这些激波辨识算法难以在复杂三维问题中广泛应用,且上述算法辨识出的激波通常为具有一定宽度的激波带,难以直接用于装配法。因此,需要发展新的辨识算法来推动激波装配法实用化进程。

本文提出了一种光源射线法,该方法可将捕捉法中探测到的具有一定宽度的激波带拟合成一系列离散线段(二维)或三角形面片(三维),利用曲线/曲面拟合算法进一步将其拟合成连续的曲线/曲面。该曲线/曲面可作为激波装配法中的激波面,并且可根据计算需要对其进行任意形式的网格划分。为便于表述,先用二维模型的实际操作过程进行原理说明,然后简单介绍三维推广过程,并展示实际工程问题的辨识结果。

1 流场结构辨识算法 1.1 二维辨识算法流程

算法的基本实现过程为:首先,通过流动特征探测技术辨识出符合流动特征要求的单元。如图 1(a)为超声速钝体扰流的流场结构,图 1(b)为利用激波探测技术在上述流场中辨识出的激波单元。上述激波单元的格心或格点构成集合T。其次,在钝体内部空间设置光源点S,设定角度间隔φ1φiφn-1,从光源处发射出来的n条光线l1liln将空间分割为n-1个子区域,射线lili+1构成区域Ai的边界,φi为射线lili+1夹角,如图 1(d)所示。单个区域Ai内的所有离散点构成子集Ti,满足T=$\bigcup\limits_{i = 1}^n {{T_i}} , {T_i}\bigcap\limits_{i \ne j} {{T_j}} = \mathit{\Phi } $。然后,统计每个区域Ai内包含的离散点个数Ni,如果某个区域内离散点数Ni少于限定数m(通常取m=3),则该区域内不进行线性拟合;如果出现相邻区域内点数都少于m,可以重新调整光源点位置或角度间隔,直到满足区域内离散点个数要求。对区域Ai内的Ni个离散点进行线性拟合(如最小二乘法)可获得线段Pi, iPi+1, i,其中Pi, j表示区域Ai内拟合出的直线与光源射线li的交点。对于内部区域,射线li分别与区域Ai-1Ai内的拟合直线产生交点Pi, i-1Pi, i。当这两个交点重合时,点Qi称为一个特征点,满足Qi=Pi, i-1=Pi, i;当这两个交点不重合时在射线li上寻找一个特征点Qi,使其满足Qi=f(Pi, i-1, Pi, i)。映射函数f可以为权函数,即Qi=k1Pi, i-1+k2Pi, i。本文取k1=k2=0.5,即特征点Qi为线段Pi, i-1Pi, i的中点。对于边界区域A1An,射线l1ln上的交点即为特征点,即Q1=P1, 1Qn=Pn, n-1。所有特征点Q1QiQn构成一个集合Q,各点之间的连接关系确定且唯一,如图 1(e)所示。最后,通过曲线拟合方法将集合Q中的元素拟合成一条连续曲线(如图 1c),可在该曲线上进行任意形式的网格划分。


图 1 光源射线法示意图 Fig.1 Sketch of source-rays method
1.2 脱体激波辨识

针对二维超声速钝头体绕流问题,采用1.1节所述方法进行激波面辨识。算例中,空气来流马赫数Ma=5.0、无量纲密度ρ=1、无量纲温度T=1,来流攻角为0°。钝头体中心位置坐标为(0, 0),钝头体半径为0.0254。网格为160×100的四边形网格,其中沿钝头体法向的节点数为100。计算网格和捕捉法流场压力云图分别如图 2(a)图 2(b)所示。


图 2 基于压力梯度的激波辨识结果 Fig.2 Shock wave surface extraction based on pressure gradient

考虑到很多高精度格式的限制器选择压力梯度来进行调节,在此也将其作为激波特征参数来探测脱体激波。对于捕捉法得到的压力场,用格林公式计算出每个网格单元的压力梯度矢量。作为间断面,激波厚度理论上处于分子自由程大小这一数量级,其两侧压力的梯度值为无穷大。但是,在捕捉法数值计算中,激波表现为带状分布的过渡区且具有较大压力梯度。对于不同的问题以及不同的计算网格,激波带上的压力梯度值都有所不同。因此,难以通过一个固定的压力梯度值将激波带提取出来。本文为了准确获得激波点集的分布,首先计算出流场中的最大压力梯度值|∇p|max,然后针对不同问题设置相应的阈值系数k,压力梯度阈值|∇p|thre=k|∇p|max,所有满足条件|∇p|≥|∇p|thre的单元格心点构成激波点集T。Case1取阈值系数k=0.1,获得离散点集,如图 2(c)所示。针对上述激波离散点集,采用光源射线法进行曲线拟合。光源S位置固定在(x, y)=(0.1, 0)处,相邻两条射线间的夹角取常值φi=4.5°。辨识出的结果如图 2(d)所示,从图中可以看出辨识出的激波曲线与捕捉法流场中的激波点集T吻合度很高,能准确描述激波形状,对激波曲线重新进行网格离散即可作为激波边界直接用于激波装配法。

离散点集的分布取决于所选择参数和阈值,为了比较不同压力梯度阈值下所辨识出的激波曲线,在捕捉法流场的基础上另外补充两个辨识状态:Case2取k=0.2,Case3取k=0.3。得到的点集及其最终辨识出来的脱体激波曲线如图 3所示。从3个状态的辨识结果可以看出,不同的压力梯度条件下所获取的离散点集有所不同,导致最终辨识出的激波曲线长度有所区别,但三条激波曲线都能与捕捉法流场中激波的位置和形状保持一致,将其作为激波装配法中的激波边界可快速得到收敛解。得到初始激波位置以后的装配法计算方法和过程参见文献[17, 19-21],这里不再赘述。


图 3 不同压力梯度阈值下的脱体激波辨识结果比较 Fig.3 Comparison of shock surface extraction results under different pressure gradient thresholds
1.3 内部激波及接触间断辨识

对波系更为复杂的马赫反射问题进行数值模拟,计算区域及边界条件如图 4所示。边长为LAF=0.4、LAB=0.1、LDE=0.312、LEF=0.85,斜楔BC倾角为14°;边界AF为超声速入口、DE为超声速出口,其余为滑移壁。计算网格为300×218的四边形网格,其中入口和出口边界上的网格节点数为300。来流马赫数Ma=1.9,攻角为0°。


图 4 计算域 Fig.4 Calculation domain

由于斜楔BC的偏转作用,自由来流在点B处产生一道入射激波(Shock1);激波达到壁面FE处时发生马赫反射,产生一道反射激波(Shock2)和一条马赫杆(Shock3);反射激波在壁面CD处再次发生反射,产生第二道反射激波(Shock4)。流场稳定后的密度云图以及波系结构如图 5所示。采用激波及接触间断探测方法可从流场中获得主要的四条激波点集和一条接触间断点集分布。为了检验不同的探测方法对最终结果的影响,激波点集采用密度梯度作为参考量进行探测,接触间断点集由密度梯度和压力梯度共同决定。同样,先计算出全场最大密度梯度值|∇ρ|max和最大压力梯度|∇p|max。激波探测的密度梯度阈值|∇ρ|threS=kρS|∇ρ|max,所有满足|∇ρ|≥|∇ρ|threS的单元节点构成激波点集,这里取kpS=0.03。接触间断单元由密度梯度和压力梯度共同决定,满足密度梯度较大且压力梯度较小,即|∇ρ|≥|∇ρ|threc、|∇p|≤|∇p|threc,其中|∇ρ|threc=kρc|∇ρ|max、|∇p|threc|=kpc|∇p|max。所有满足该条件的单元节点构成曲线拟合所需的接触间断点集。这里取kρc=0.05、kpc=0.035,其中在不同密度梯度阈值|∇ρ|threc下辨识出的接触间断长度有一定区别,但位置和形状基本一致。激波和接触间断离散点集分布如图 5所示。


图 5 捕捉法密度云图及间断点集分布 Fig.5 Density contours and discontinuity point set distribution under shock-capturing schemes

获得离散点集后,通过在区域内合理布置光源可拟合出每条激波和接触间断曲线。首先在Ⅱ区设置一个光源,光源发出等间距射线,采用上文光源射线法可辨识出表征激波1和激波2的曲线,如图 6(a)所示。同理,在Ⅲ区布置光源可辨识出激波4和接触间断曲线,如图 6(b)所示。在I区设置光源可辨识出马赫杆曲线,如图 6(c)所示。所有曲线的最终辨识结果与流场密度云图的比较如图 7所示。可以看出,辨识出的曲线与采用捕捉法获得的波系结构高度吻合,直接用于装配法有助于快速得到收敛解。得到初始激波位置以后的装配法计算方法和过程参见文献[17, 19-21],后文也会有简单介绍,这里不再赘述。


图 6 激波和接触间断辨识过程 Fig.6 Shock wave and contact discontinuity surface extraction process


图 7 激波及接触间断辨识结果与流场比较 Fig.7 Comparison between feature surface extraction result and flow field structure

值得注意的是,这里光源的位置分布并不唯一,可采用不同的分布方式达到相似的效果,比如可先在Ⅰ区设置光源同时拟合出激波1和激波3曲线,再在Ⅲ区设置光源得到激波2、激波4和接触间断曲线。

2 三维实现 2.1 三维辨识算法流程

三维问题的辨识过程与二维类似,首先通过流场探测技术获取捕捉法流场中相应流动特征的离散点集,然后在离散点集的基础上利用光源射线法进行特征曲面的拟合。三维问题中射线l的方向在直角坐标系中由两个角度决定,分别为与z轴的夹角θ以及lxy平面内的投影与x轴的夹角φ,如图 8所示。通过给定一系列角度{θi, φi},离散点集T所在的空间区域A被射线集{li}分割成n个子区域。每个子区域Ap由三条射线liljlk构成,满足$ A = \bigcup\limits_{p = 1}^n {{A_p}} $${A_p}\bigcap\limits_{p \ne q} {{A_q}} = \mathit{\Phi } $。子区域Ap中包含的所有离散点构成集合T的子集Tp,满足$ T = \bigcup\limits_{p = 1}^n {{T_p}} , {T_p}\bigcap\limits_{p \ne q} {{T_q}} = \mathit{\Phi }$。上述过程中,确定合适的射线连接关系,使子区域之间满足连续且不相交条件比较困难。


图 8 射线方向示意图 Fig.8 Ray direction

本文通过引入单位球可确定各射线之间的连接关系,如图 9所示。对半径为1的球面进行三角形离散,获得一系列三角形面片,各顶点连接关系确定,且三角形面片满足连续且不相交条件。


图 9 单位球 Fig.9 Unit sphere

将单位球的球心作为光源S,球心与单位球上三角形面片的三个顶点aiajak的连线liljlk构成一个子区域Ap,如图 10(a)所示。若Ap内离散点的个数不小于3,则对Ap内所有离散点进行平面拟合,拟合出的平面与三条射线的交点分别为bi, pbj, pbk, p;若该子区域内离散点的个数小于3,则跳过该区域。所有区域拟合完成后,每条射线li与不同子区域内的拟合平面形成若干交点。将各射线li上的交点取平均值获得特征点bi,则各子区域Ap的三条射线边界上的特征点bibjbk构成目标三角形面片,如图 10(b)所示。最终,整个离散点集T被拟合成一系列连续的三角形面片,三角形面片顶点的连接关系由单位球确定。可对上述三角形面片构成的空间曲面重新进行网格划分,获得可用于装配法的激波面网格。


图 10 子区域内三角形面片拟合过程 Fig.10 Triangular fitting process in sub-region
2.2 三维测试算例

通过以下两个算例验证本文辨识算法在三维问题中的有效性。

算例一:对三维半椭球面进行辨识。椭球面方程为$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1, (z \geqslant 0) $,其中a=b=2, c=3。在如图 8所示坐标系下,椭球面上的点A(x, y, z)满足以下关系式:

$ \left\{\begin{array}{l}{x=a \sin \theta \cos \varphi} \\ {y=b \sin \theta \cos \varphi} \\ {z=c \cos \theta}\end{array}\right. $

在椭球面附近随机生成100 000个点作为离散点集,随机点坐标B(x′, y′, z′)满足:

$ \left\{\begin{array}{l}{x^{\prime}=a^{\prime} \cos \theta^{\prime} \cos \varphi^{\prime}} \\ {y^{\prime}=b^{\prime} \cos \theta^{\prime} \sin \varphi^{\prime}} \\ {z^{\prime}=c^{\prime} \sin \theta^{\prime}}\end{array}\right. $

其中,

$ \left\{\begin{array}{l}{a^{\prime}=(0.97+0.06 k) a} \\ {b^{\prime}=(0.97+0.06 l) b} \\ {c^{\prime}=(0.97+0.06 m) c} \\ {\varphi^{\prime}=\alpha_{1}+i\left(\alpha_{2}-\alpha_{1}\right), \quad 0 \leqslant \alpha_{1}, \quad \beta_{2} \leqslant 2 \pi} \\ {\theta^{\prime}=\beta_{1}+j\left(\beta_{2}-\beta_{1}\right), \quad 0 \leqslant \beta_{1}, \quad \beta_{2} \leqslant 0.5 \pi}\end{array}\right. $

式中,ijklm为区间[0, 1]内的随机数。

获得的离散点集如图 11所示。单位球共有736个节点、1468个三角形单元(如图 9),球心位置坐标为(0, 0, 0)。单位球位置及椭球面辨识结果如图 12所示。从结果可以看出,辨识出的曲面光滑完整且形状与离散点分布形状一致。


图 11 椭球随机点集 Fig.11 Random point set of ellipsoid


图 12 单位球及拟合结果 Fig.12 Unit sphere and the extracted surface of ellipsoid

算例二:对三维双椭球面进行拟合。椭球面一满足方程$\frac{x^{2}}{2^{2}}+\frac{y^{2}}{5^{2}}+\frac{z^{2}}{2^{2}}=1, (y \leqslant 0) $,椭球面二满足方程$ \frac{x^{2}}{2.5^{2}}+\frac{y^{2}}{3^{2}}+\frac{z^{2}}{2.5^{2}}=1, (y \leqslant 0)$两椭球交集的外边界作为目标曲面。采用与算例一类似的随机函数分别获取两个椭球面附近的随机点集坐标,其中π≤α1, α2≤2π、0≤β1, β2≤π,得到目标离散点集如图 13所示。采用与算例一相同的单位球,球心位置坐标为(0, -1.5, 0)。双椭球面辨识结果如图 14所示,从图中可以看出区域边缘附近的拟合曲面不太完整,但其他部位的拟合曲面连续且光滑。


图 13 双椭球随机点集 Fig.13 Random point set of double-ellipsoid


图 14 双椭球面拟合结果 Fig.14 Extracted surface of double-ellipsoid
2.3 三维绕流激波辨识

对三维超声速球头绕流进行模拟,来流马赫数Ma=5.0, 无量纲密度ρ=1.0, 无量纲温度T=1.0,来流攻角和侧滑角均为0°。球头的球心位置坐标为(0, 0, 0),球头半径为0.0254。计算网格共有508 365个的六面体单元,球面网格如图 15(a)所示。图 15(b)为流场压力云图,图 15(c)为根据流场压力梯度获得的激波带离散点集。算例中,压力梯度阈值系数k=0.05离散点为满足|∇p|≥|∇p|thre的所有单元节点,共有23 640个。单位球共包含1734个节点,3464个三角形面片,球心坐标为(2, 0, 0)。最终辨识出的激波面如图 15(d)所示,除了边缘附近外的其他区域辨识结果与离散点集吻合较好,将边缘位置稍作处理即可作为初始激波面用于激波装配法。对于不同的初始位置的激波,最终收敛解相同,故本例不再展示装配法的计算过程,计算方法参见文献[17, 19-21]。


图 15 三维超声速球头绕流激波辨识过程 Fig.15 Shock wave surface extraction process of 3D blunt body flow

在球头绕流算例的基础上,对波系较为复杂的球柱锥组合体外形进行模拟,检验特征面提取方法处理三维复杂问题的能力。计算模型如图 16(a)所示,来流马赫数Ma=4.04,来流攻角为20°。波系结构如图 16(b)所示,在最外层形成一道弓形激波,在锥裙与柱体的结合处形成一个内嵌激波。探测出的激波点集和提取出的激波面分别如图 16(c)图 16(d)所示。


图 16 球柱锥组合体绕流激波辨识过程 Fig.16 Shock wave surface extraction process of cylinder with a hemispherical nose and a conical flow
3 装配法实例

利用上述方法获得初始激波面后,将其作为激波边界用于激波装配法,进行流动数值模拟,考查所述激波辨识方法在实际使用中的效果。

对于定常激波装配法,初始激波位置及形状会影响计算的收敛过程,但最终的收敛结果都相同。因此,本节直接采用文献[19]算例的部分结果说明激波装配法与捕捉法的区别。该算例为具有复杂波系结构的马赫反射问题,来流马赫数2.0,偏转角δ=14°。计算区域为L×L的正方形,边界上网格节点均匀分布,节点间距为0.01L,内部区域为三角形网格,单元数为22 429。流场波系结构及边界条件如图 17所示,入射激波IS和反射激波RS、马赫杆MS交于三波点TP,SS为接触间断。


图 17 边界条件及流场结构示意图[19] Fig.17 Boundary conditions and sketch of flow field structure[19]

首先采用激波捕捉法获得初始流场,在此基础上采用本文辨识算法辨识出激波及接触间断的位置和形状,将其作为装配法的初始间断面,结合R-H关系式进行装配法流场计算。得到初始激波位置以后的计算过程参见文献[19],计算结果完全一致。下面给出部分结果,如图 18所示,其中图 18(a)为采用捕捉法的计算结果,图 18(b)为采用装配法的计算结果。从图中可以看出,捕捉法流场中激波附近出现等值线的扭曲和震荡,且激波呈带状分布,具有很宽的过渡区。与此相比,采用装配法所获得的流场中各区域之间界线分明,激波是一个没有宽度的界面,激波两侧变量值严格满足R-H关系式。


图 18 装配法与捕捉法计算结果对比(马赫数云图)[19] Fig.18 Comparison between shock-capturing solution and shock-fitting solution[19]
4 结论

通常对于具有一定“厚度”的空间点集,采用传统的拟合方法难以将其拟合成形状一致的光滑曲面,本文提出的光源射线法很好地解决了此类问题。

对于从捕捉法流场中获取的具有一定“厚度”的激波带和接触间断带,采用该方法进行空间曲面拟合后可作为激波装配法中的初始间断面。尤其对于具有复杂波系结构的流场,该方法提供了一种装配初始流场的途径。本文详细介绍了该算法在二维问题中的实现过程,并将其成功推广应用于三维问题。算例结果表明,辨识出的激波曲面整体上与捕捉法流场中激波带分布相吻合,但是三维问题中激波曲面边缘附近的辨识效果有待进一步提高。对于复杂工程问题,往往存在较为复杂的三维波系结构,间断面的提取过程更为困难和繁琐,该方法还需进一步完善。

对于传统激波装配法,在初始激波面给定后需要经过若干步迭代才能收敛到准确的形状和位置。尤其对于三维问题,收敛过程更加漫长,甚至出现难以收敛的情况。采用本文方法可提供较为准确的初始激波,有助于缩短收敛过程。

本文所提方法除了用于激波装配法外,未来还可应用于自适应网格算法。根据不同的流动特征可拟合出相应的空间曲面,进而可以任意调整该曲面附近区域的网格分布,达到网格根据流动特征自动调整的目的。

参考文献
[1]
LAX P D. Weak solutions of nonlinear hyperbolic equations and their numerical computation[J]. Communications on Pure and Applied Mathematics, 1954, 7(1): 159-193. DOI:10.1002/(ISSN)1097-0312
[2]
LAX P, WENDROFF B. Systems of conservation laws[J]. Communications on Pure and Applied mathematics, 1960, 13(2): 217-237. DOI:10.1002/(ISSN)1097-0312
[3]
张涵信, 沈孟育. 计算流体力学-差分格式原理和应用[M]. 北京: 国防工业出版社, 2003.
ZHANG H X, SHEN M Y. Computational fluid dynamics-theory and application of difference scheme[M]. Beijing: National Defense Industry Press, 2003. (in Chinese)
[4]
VON NEUMANN J, RICHTMYER R D. A method for the numerical calculation of hydrodynamic shocks[J]. Journal of Applied Physics, 1950, 21(3): 232-237. DOI:10.1063/1.1699639
[5]
KUNDU P K. Fluid mechanics[M]. San Diego: Academic Press, 1990.
[6]
赵学瑞, 廖其奠. 黏性流体力学[M]. 北京: 机械工业出版社, 1983.
ZHAO X R, LIAO Q D. Viscous fluid mechanics[M]. Beijing: China Machine Press, 1983. (in Chinese)
[7]
吴其芬, 陈伟芳, 等. 稀薄气体动力学[M]. 湖南长沙: 国防科技大学出版社, 2004.
WU Q F, CHEN W F. Rarefied gas dynamics[M]. Changsha: National University of Defense Technology Press, 2004. (in Chinese)
[8]
ABBETT M, MORETTI G. A time-dependent computational method for blunt body flows[J]. AIAA Journal, 1966, 4(12): 2136-2141. DOI:10.2514/3.3867
[9]
MORETTI G. Thirty-six years of shock fitting[J]. Computers & Fluids, 2002, 31(4): 719-723.
[10]
张德良. 计算流体力学教程[M]. 北京: 高等教育出版社, 2010.
ZHANG D L. A Course in Computational Fluid Dynamics[M]. Beijing: Higher Education Press, 2010. (in Chinese)
[11]
RICHTMYER R D, MORTON K W. Difference methods for initial-value problems[M]. New York: Interscience Publisher, 1967.
[12]
RAWAT P S, ZHONG X L. On high-order shock-fitting and front-tracking schemes for numerical simulation of shock-disturbance interactions[J]. Journal of Computational Physics, 2010, 229(19): 6744-6780. DOI:10.1016/j.jcp.2010.05.021
[13]
PRAKASH A, PARSONS N, WANG X, et al. High-order shock-fitting methods for direct numerical simulation of hypersonic flow withchemical and thermal nonequilibrium[J]. Journal of Computational Physics, 2011, 230(23): 8474-8507. DOI:10.1016/j.jcp.2011.08.001
[14]
PACIORRI R, BONFIGLIOLI A. A shock-fitting technique for 2D unstructured grids[J]. Computers & Fluids, 2009, 38(3): 715-726.
[15]
PACIORRI R, BONFIGLIOLI A. Shock interaction computations on unstructured, two-dimensional grids using a shock-fitting technique[J]. Journal of Computational Physics, 2011, 230(8): 3155-3177. DOI:10.1016/j.jcp.2011.01.018
[16]
BONFIGLIOLI A, GROTTADAUREA M, PACIORRI R, et al. An unstructured, three-dimensional, shock-fitting solver for hypersonic flows[J]. Computers & Fluids, 2013, 73: 162-174.
[17]
刘君, 邹东阳, 徐春光. 基于非结构动网格的非定常激波装配法[J]. 空气动力学学报, 2015, 33(1): 10-16.
LIU J, ZOU D Y, XU C G. An unsteady shock-fitting technique based on unstructured moving grids[J]. Acta Aerodynamica Sinica, 2015, 33(1): 10-16. DOI:10.7638/kqdlxxb-2014.0109 (in Chinese)
[18]
刘君, 徐春光, 白晓征. 有限体积方法和非结构动网格[M]. 北京: 科学出版社, 2016.
LIU J, XU C G, BAI X Z. Finite volume methods and unstructured dynamic meshes[M]. Beijing: Science Press, 2016. (in Chinese)
[19]
ZOU D, XU C, DONG H, et al. A shock-fitting technique for cell-centered finite volume methods on unstructured dynamic meshes[J]. Journal of Computational Physics, 2017, 345: 866-882. DOI:10.1016/j.jcp.2017.05.047
[20]
刘君, 邹东阳, 董海波. 动态间断装配法模拟斜激波壁面反射[J]. 航空学报, 2016, 37(3): 836-846.
LIU J, ZOU D Y, DONG H B. A moving discontinuity fitting technique to simulate shock waves impinged on a straight wall[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(3): 836-846. (in Chinese)
[21]
刘君, 邹东阳, 董海波. 基于非结构变形网格的间断装配法原理[J]. 气体物理, 2017(1): 13-20.
LIU J, ZOU D Y, DONG H B. Principle of new discontinuity fitting technique based on unstructured moving grid[J]. Physics of Gases, 2017(1): 13-20. (in Chinese)
[22]
WU Z N, XU Y Z, WANG W B, et al. Review of shock wave detection method in CFD post-processing[J]. Chinese Journal of Aeronautics, 2013, 26(3): 501-513. DOI:10.1016/j.cja.2013.05.001
[23]
PAGENDARM H G, SEITZ B. An algorithm for detection and visualization of discontinuities in scientific data fields applied to flow data with shock waves[J]. Scientific Visualization:Advanced Software Techniques, 1993, 161-177.
[24]
LOVELY D, HAIMES R. Shock detection from computational fluid dynamics results[C]//14th Computational Fluid Dynamics Conference. 1999: 3285. https://www.researchgate.net/publication/2581120_Shock_Detection_from_Computational_Fluid_Dynamics_Results
[25]
KANAMORI M, SUZUKI K. Shock wave detection in two-dimensional flow based on the theory of characteristics from CFD data[J]. Journal of Computational Physics, 2011, 230(8): 3085-3092. DOI:10.1016/j.jcp.2011.01.007
[26]
KANAMORI M, SUZUKI K. Shock wave detection based on the theory of characteristics for CFD results[C]//20th AIAA Computational Fluid Dynamics Conference. Honolulu, Hawaii, 2011. https://www.researchgate.net/publication/268474447_Shock_Wave_Detection_based_on_the_Theory_of_Characteristics_for_CFD_Results