一种求解定常无黏流场的混合延拓方法

乔磊, 白俊强, 邱亚松, 华俊,2, 徐家宽

(1.西北工业大学 航空学院, 陕西 西安 710072; 2.中国航空工业集团公司 中国航空研究院, 北京 100012)

在飞行器气动特性的评估和优化设计中，通常需要对大量类似的几何构型反复进行定常流动模拟。因此，为使计算流体力学方法在飞行器设计中发挥更有力的作用，迫切需要研究具有更高计算效率的定常流场数值模拟计算方法。由于方程的非线性性质，其求解必须通过迭代法实现。根据迭代方法构造方式的不同，可以分为显式方法和隐式方法。由于定常问题的求解不受时间精度的限制，隐式解法可以充分发挥稳定性强、时间步长大的优点，从而在定常问题的求解中得到了广泛的应用。

在非线性问题的隐式解法中，需要通过牛顿法对问题进行线性化处理。牛顿迭代由于可以达到快速的平方收敛而成为应用广泛的非线性求解方法，但是牛顿迭代要求初始误差足够小才能避免发散。不像非定常计算可以利用已知时间步的信息构造质量较高的初值[1]，在定常空气动力学问题中，一般很难得到足够精确的初始解。为解决这一矛盾，人们一般通过对需要求解的问题进行修改，以改善问题的数学性质，使初始解得以满足牛顿迭代的收敛条件。这种扩大牛顿迭代可接受初值条件范围的操作，就称为牛顿迭代的全局化(globalization)。

根据具体操作对象的不同，牛顿迭代全局化方法可以被分为2种类型。第一类方法不涉及非线性问题本身，而是试图通过修改离散方式或解的修正量来保证稳定性。比较常见的有网格序列法[2]和线性搜索法[3]。这类方法不是本文研究的重点，故而不在这里作详细介绍。因为这一类方法独立于非线性问题本身，所以它们实际上可以与下面要介绍的第二类方法结合使用。第二类方法需要修改所要求解的问题本身来提高稳定性和计算效率。对于飞行器流场模拟问题，这类方法的典型代表有边界松弛(boundary condition relaxation)法和方程延拓(equation continuation)法[4]。由于定常流动问题的初场一般根据来流条件设置为均匀场，故而初始误差将集中在物面边界处。因此，通过边界松弛法逐步施加边界条件以改善问题的稳定性就成为一种顺理成章的选择。Lyra等人[5]对边界松弛法有较多的研究和应用。

方程延拓法是通过修改控制方程本身来改善非线性迭代收敛性的方法。方程延拓法中最为常见的就是经典的伪时间延拓法(PTC，pseudo time continuation)。为避免延拓项污染方程的解，方程延拓法需要建立一个延拓项逐渐减小并最终消失的方程序列，通过逐次求解近似方程最终得到原方程的解。伪时间推进法作为方程延拓法的一种，它从非定常方程的物理意义出发，在定常控制方程中引入了一个伪时间项，随着求解的收敛，时间项的作用逐渐消失，从而得到定常问题的解。作为一种经典方法，伪时间推进法得到了非常广泛的应用，也受到了研究者较多的关注。为提高迭代的稳定性，Kelley等人[6]发展了一种带约束的伪时间推进法。

除加入时间项外，引入人工耗散来提高解的稳定性也是一种常用的方法。早在速势方程的时代，Young等人提就出了用于加速收敛的黏性松弛法[7]。近年来，Hicken等人对比研究了黏性延拓和伪时间推进法在定常RANS方程求解中的效果，认为“黏性延拓法是伪时间推进法的一种具有相当鲁棒性和效率的替代方法”[8]，但是没有研究将二者进行结合的方法。最近，针对拟线性的标量对流输运方程，Pollock发表了一种通过拉普拉斯算子改善病态雅克比的方法[9]。不过，根据Hicken等人[8]的研究结论，由于NS方程的连续方程不存在黏性项，所以这种方法增稳作用比较有限。

受Pollock和Hicken等人工作的启发，本文对拉普拉斯算子和伪时间推进法相结合的混合延拓方法进行了研究，包括延拓方程的构造、延拓参数的初值选择、推进策略以及终止条件等。鉴于欧拉方程可以比较充分地体现控制方程的非线性特征，也可以反映激波这一引起流场求解不稳定的关键元素，本文在欧拉方程有限元求解器中实现了所构造的混合延拓方法，并通过GAMM鼓包和NACA0012翼型算例对所构造的混合延拓方法的效率进行了验证。

1 基本问题和求解方法

定常可压缩欧拉方程的隐式求解涉及一整套复杂的计算方法，其中主要包括方程延拓、欧拉方程的稳定化以及牛顿迭代中线性问题的求解，而各部分之间又互相关联，影响算法的整体效率。故虽然本文工作重点不在于欧拉方程的稳定化和线性问题的求解，为表述的清晰和完整起见，此处仍须对相关关键信息进行必要介绍。

1.1 基本控制方程

本文所研究问题的控制方程是无量纲的定常可压缩欧拉方程，表示为

R(w)=·Fc(w)=0

(1)

式中，w=(ρ,U,p)表示密度、速度和压力等流场自变量,R(w)表示残差,Fc(w)为无黏对流通量,具体为

(2)

式中，⊗表示直积,I为单位阵,E表示流体总能。

1.2 稳定化和激波捕捉

众所周知,在流动中出现不光滑结构时,使用常规有限元方法求解欧拉方程将会产生非物理的震荡,因而需要通过一些特殊处理解决这一问题。本文通过引入受人工黏性系数控制的人工黏性通量来完成方程的稳定化以及激波捕捉。方程的具体表述如下:

Rflux(w)=·[Fc(w)-Fav(w)]=0

(3)

虽然最初的问题是求解欧拉方程,但最终实际上被求解的,是加入黏性通量后的稳定化欧拉方程(3)。其中人工黏性通量的具体形式仿照Navier-Stokes方程定义如下:

(4)

式中，τij为应力张量,T为温度,μE和κE为待定义的人工黏性系数和人工导热系数。人工导热系数与人工黏性系数之间的关系也依物理关系确定

κE=PrμE

(5)

式中，Pr为普朗特数,本文取0.72。至此,唯一的待定参数就是人工黏性系数μE。

本文应用参考文献[10]中的方法来确定人工黏性系数。思路是,对于流动中的比较不光滑的区域,有限元解中的最高阶分量将占据较多的成分,反之亦然。根据这一指标构造相应的指示因子,就可以在流场中需要的部分加入合适的人工黏性。由于参考文献[10]中的方法是针对采用正交基函数的间断有限元方法,而本文采用的是基于拉格朗日基函数的连续有限元方法,故此处对与参考文献有所不同的一些细节进行交代。

在单元ΩK上,震荡指示因子FK定义如下

(6)

得到FK后,通过如下分段函数计算各单元上的人工黏性系数μE,K

(7)

式中，μE,l,Fl和Fu都是经验参数。相对于参考文献[10],此处引入了背景黏性系数μE,l。μE,l在大部分情况下取0,在某些特别光滑的情形中,震荡指示因子不足以给出适当的人工黏性,才需要引入非0的背景黏性。Fl和Fu是震荡指示因子的上下限,在无特别说明时,本文中二者的取值分别为-6.5和-3.5。参数μE,u,K是当地流场震荡严重时的人工黏性系数,这里与参考文献[10]中相同,取经典的一阶黏性:

(8)

此处cE是经验系数,本文取0.1。‖u‖∞,K和‖a‖∞,K为单元ΩK上速度和声速的无穷范数,在计算中取为单元高斯积分节点中的最大值。参数he,K为有效单元尺度,定义为

he,K=min(hK,rhhmin)

(9)

式中，hK为单元实际尺度,hmin为全局最小网格尺度,rh为约束系数,本文取值为3.0。对有效网格尺度进行约束的原因是为了避免在外流问题远场处相对巨大的网格中产生病态的人工黏性系数。

至此,通过人工黏性方法进行稳定处理的控制方程介绍完毕。为便于后文叙述,这里定义引入人工黏性后的通量为

Fflux≡Fc-Fav

(10)

1.3 加入延拓项

在方程延拓法求解非线性方程的过程中,除最后一步外,实际上被求解的都是延拓方程而非原始方程本身。这里,为了表述方便,引入形式上的延拓项,其具体定义在第2节中再进行详细讨论。由于延拓方程构成一个方程序列,我们记第n个延拓项为Cn(w),则延拓方程及其残差项可记为

(11)

1.4 有限元离散及数值求解

本文应用Q1线性拉格朗日单元对控制方程进行离散。对于给定计算域Ω和网格划分TK,由Q1单元可以建立解函数空间Vh,以任意测试函数zh∈Vh乘以延拓方程(11)并分部积分,即可得到弱形式离散方程

(12)

由于下文将只讨论离散计算问题,故在不引起混淆的情况下,不再以下标h区分连续方程和离散方程。边界条件通过界面通量Fflux(wh)∂Ω以弱形式施加给离散方程。界面通量的计算需要根据边界条件重构界面自变量。在本文中,物面处使用无穿透条件,远场使用黎曼边界条件。在压力入口处指定压力分量,外插其他自变量。边界条件的实现细节可以参考文献[11]。

建立完整的离散问题后,每个延拓方程都仍然是非线性方程,需要通过牛顿迭代进行求解。但是,由于延拓方程的解一般不是原始方程的解,没有必要精确求解各个延拓方程。在延拓方程求解时,延拓方程的残差范数或原始方程的残差范数降低一个数量级后立即停止当前延拓方程的求解。在最后一步延拓迭代中求解原始方程时,计算收敛的标准设为残差小于1×10-11。

牛顿迭代中的线性系统通过PETSc库[12]进行程序实现,相应的矩阵方程通过MUMPS[13]并行直接求解器进行求解。有限元离散到线性系统的数据接口通过deal.II有限元库[14]实现。

2 混合延拓方法的构造

2.1 拉普拉斯算子与伪时间项的对比分析

在定常流动的隐式解法中,线化处理的核心是牛顿迭代法。对于(12)式表示的问题略去延拓序号上标n,牛顿迭代可以表示为:

fork=0,1,2,…

(13)

牛顿迭代的收敛对方程的性质和初始解有较苛刻的要求:必须满足雅克比矩阵非奇异,初始解足够接近方程解。所谓“足够接近”,具体是指对于任意第k步牛顿迭代,须满足

(14)

式中，R″和R′分别为非线性算子的海森矩阵和雅克比矩阵,ε为当前解的误差。因而,改善牛顿解法收敛性的一般思路是降低雅克比矩阵的奇异性,以及增强问题的光滑性。作为一种常见的特例,伪时间推进法通过引入伪时间项,所求问题转化为

(15)

在实际计算中,Q(w)可以取为控制方程守恒变量,也可以简单取w本身。出于简化分析的考虑,本文取w。另外,τ是伪时间。这样,牛顿迭代转化为

(16)

式中，dτ为伪时间步长。由于对角阵的引入,所解问题的线性增强,使牛顿迭代更容易避免发散。

不过,由(14)式可知,要保证高牛顿迭代的收敛,在给定初始解的前提下,根本途径是降低雅克比矩阵的奇异性,增加问题的线性。在这些方面,任意非奇异的线性算子,都能起到降低海森矩阵范数,同时避免雅克比矩阵奇异的作用。正因为如此,拉普拉斯算子才可以成为伪时间项的一个有竞争力的替代方案。进一步引入拉普拉斯项的控制方程如(17)所示。其中cLP是拉普拉斯项的缩放系数,2表示拉普拉斯算子。

(17)

这样,牛顿迭代的公式变为

(18)

虽然拉普拉斯算子和伪时间项引入的对角阵都是对称正定矩阵,但是针对通常航空空气动力学问题中初始残差仅分布在物面附近的特点,时间项和拉普拉斯算子的作用之间的差异就比较大了。

由于伪时间项在雅克比矩阵中反映为一个对角阵,其逆算子仍然是对角阵,显然它对残差矢量的响应是当地的、局部的。而拉普拉斯算子则不同。通过格林函数法,知道对于残差向量R(w(x)),在三维拉普拉斯逆算子作用下的响应是

(19)

这说明任意一点的残差都会对解向量产生全局性的影响;同时任意一点的解向量变化,都包含了整个求解域中的残差信息。这正是拉普拉斯算子椭圆性的体现。在非线性迭代的初始阶段,由于残差仅在物面处非零,故拉普拉斯算子向流场内部传递边界信息的效果要远远强于伪时间项,解的残差也就可以更快地降低。这样,由(14)式可知,后续牛顿迭代的稳定性也就更容易保证,从而延拓项也可以更快速地衰减,最终得到更快的收敛速度。

但是,相对于伪时间推进法,拉普拉斯延拓法还有其不足之处。其中一点就是必须特别设计一种机制,使拉普拉斯项随迭代推进而逐步消失,否则将无法得到原始问题解。但是,如果使拉普拉斯项衰减过快,则有可能因为稳定性不足而导致发散。这就要求拉普拉斯延拓法必须在保守推进策略与激进推进策略之间寻求一个合适的平衡点,使得其应用难度高于时间推进法。另外,当流场中存在激波时,较大的拉普拉斯项将由于过度耗散而抹平激波,给收敛速度带来较大的不利影响;而过小的拉普拉斯项会使牛顿迭代稳定性不足而发散。在这种情况下,时间推进法可以通过足够小的时间步长可靠地使迭代收敛。

基于上述分析,可以看到,拉普拉斯项主要的优势在于迭代的最初始阶段。因此,可以构造一个混合延拓方法,拉普拉斯项在初始阶段占主导,而后逐渐过渡为伪时间推进法,以求综合二者的优势,达到更高的计算效率。

2.2 混合延拓项的构造

这里将混合延拓项定义为拉普拉斯算子的和伪时间项的线性组合。对于第n步延拓中的第k步牛顿迭代,延拓项的计算方法为

(20)

式中，λ为延拓项整体控制因子,ω为混合权重,两者的初始值都需要通过输入给定。上节分析中提到,当流场中出现激波结构时,应当提高伪时间项的权重。而在本文计算方法中,人工黏性系数是现成的流场光滑性指示因子。因此,对于第一步延拓之后的计算,定义混合权重ω的计算方法为

∀n>1

(21)

2.3 迭代策略的设计

本节给出用于计算各步λn的延拓项推进算法,包括递推算法、延拓参数的初始值和终止值。除关于混合延拓方法的的参数外,还有作为对照方法的时间推进法和单纯拉普拉斯延拓法的延拓策略。

首先是延拓参数初值的选择。对于拉普拉斯延拓而言,它的主要目的是耗散初始残差。根据参考文献[15]的结论,无量纲化的速度散度的值与马赫数的平方及压力梯度成正比,故此处依据来流马赫数的平方来设置拉普拉斯延拓的初始系数。

(22)

式中，M∞为来流马赫数,cc经验系数本文取0.1。在混合延拓方法的初期,拉普拉斯算子占主导地位,其初始延拓系数也通过上述方法设置。不过,时间推进法的初始CFL数难以找到先验的设定依据。由于时间推进法在本文中作为对照方法,所以其初始CFL是通过多次尝试得到的不至于发散的最大CFL数,这样其他延拓方法所取得较高计算效率才具有客观性。

延拓项系数的推进算法主要借鉴成熟的SER(switched evolution relaxation)[8]CFL数递增策略。本文对延拓项系数λn和CFL数η的倒数采用同样的递减模式,如(23)式所示

(23)

式中，上标n表示迭代步数,‖Rflux‖L2是残差向量的L2-范数。延拓参数κ0,β∈[0.5,1.5]和Rmax∈[0.5,1.0]是根据经验确定的常数。由于本文所使用的非线性迭代方法中没有保证残差单调降低的机制,因此需要引入强制缩减因子Rmax防止收敛停滞的发生。

最后是延拓的终止。由于指数递减的渐近性质,当残差无法下降到0时,延拓参数也不可能彻底为0。为避免残存的延拓项消耗不必要的计算资源,本文对“足够小”的延拓项系数采取截断处理。对于拉普拉斯延拓,“足够小”的判定是基于人工黏性实现的。由于本文利用人工黏性实现欧拉方程的稳定化,当延拓项的拉普拉斯算子贡献不足人工黏性贡献的1/10时,有理由认为延拓项已经没有存在的必要,可以直接将其置零,即

(24)

式中，cr=0.1。时间推进法中的延拓项性质与黏性不同,无法直接与人工黏性进行比较,对于伪时间延拓,通过时间项对残差的贡献来确定延拓项是否已经“足够小”。混合延拓法最末期由时间项主导,因此其终止判据沿用时间推进法的判据:

(25)

同样ct=0.1。

3 算例验证

由于本文使用稀疏矩阵直接求解器,其计算代价只与矩阵的规模和稀疏形式(sparsity patterm)有关,而这些参数都不受延拓方式影响。故本文通过各方法所需要的线性求解步数来衡量其计算代价。

3.1 算法参数汇总

由于本文用到的人工黏性和方程延拓方法中涉及较多的参数,为清晰起见,这里给出一个汇总说明。下文某些算例对模型参数有所调整,将会相应地给出明确说明。其中人工黏性涉及的经验参数如表1所示,延拓算法的参数如表2所示。

表1 无特别说明情况下的人工黏性参数

表2 无特别说明情况下的延拓算法参数

3.2 GAMM鼓包

本节给出GAMM鼓包[15]算例中3种延拓方法的效率对比。本文选取来流马赫数0.675的一个有激波存在的算例,用来验证混合延拓法的计算效率。图1给出了GAMM鼓包算例的计算网格。计算网格在鼓包的前后缘点以及激波可能出现的区域进行了预先加密。其中计算网格包含19 754个四边形单元。

图1 GAMM鼓包算例超临界状态的网格

本算例收敛解的马赫数云图在图2中给出，从图中可以看出,收敛解中没有明显的非物理震荡,证明本文所选用的人工黏性方法起到了预期的作用。

图2 GAMM鼓包算例的马赫数云图

图3给出了GAMM鼓包算例通量残差Rflux随线性求解步数Niter的收敛历史。图中PTC代表伪时间延拓法,LC代表拉普拉斯延拓法,BC代表混合延拓(blended continuation)法。其中,PTC方法所用的初始CFL数分别为η0=5。在小本节算例中,从整体上来看,PTC方法的收敛效率最低。在迭代初期,LC和延拓项受拉普拉斯算子主导的BC方法都能较快地使方程残差收敛,随着迭代的收敛和激波结构的清晰化,BC方法的收敛效率逐渐超越了LC方法。而由于PTC方法的耗散仅由人工黏性提供,其收敛速度随着马赫数增加不断降低。但是为了保证解的准确性,人工黏性只能保持在尽可能低的水平。拉普拉斯算子的引入正好解决了这一矛盾。

图4给出了延拓参数的变化过程。图中“CL-LC”表示拉普拉斯延拓法中拉普拉斯项的系数,“CT-BC”表示混合延拓方法中的时间项系数,“CL-BC”表示混合延拓方法中的拉普拉斯项系数,“CC-BC”表示混合延拓方法中这个延拓项的系数。从图中可以看出LC和BC两者效率有所差异的原因。随着迭代的进行,伪时间项逐步主导延拓项。而在迭代的后期,一方面由于激波这一尖锐结构的出现,另一方面由于解逐渐靠近问题的最终解,拉普拉斯算子的耗散和全局作用特性都对收敛不再有促进作用。因此,由于BC法通过把延拓项的主导成分变成伪时间项,快速脱离了拉普拉斯算子,达到了较快的收敛效率。

图5给出了本文计算得到的上下壁面马赫数分布与Žaloudek和Janda等[16-17]所给出结果的对比,可见符合程度较好,说明本文所发展的计算方法与程序实现是正确的。

图3 GAMM鼓包算例通量残差的收敛历史 图4 GAMM鼓包算例延拓项系数的变化历史 图5 GAMM鼓包算例马赫数分布的对比

3.3 NACA0012翼型

本小节NACA0012算例[18]对混合延拓方法在外流问题中的表现进行测试。算例中翼型的弦长L=1,远场的设置以翼型前缘点为原点,流向和纵向的范围都是[-25L,25L]。测试算例的来流马赫数为0.8,迎角为1.25°。图6给出了本节算例所使用的计算网格。与GAMM鼓包算例类似,本算例计算网格也在流动状态变化可能比较剧烈的地方进行了预先的局部加密处理。计算网格包含192 100个四边形网格单元。计算域的外围边界应用自由来流边界条件,翼型表面应用无黏壁面边界条件。

图6 NACA0012算例计算网格示意图

图7给出了NACA0012翼型算例的压力系数云图。从结果中可以看出,即使是在有激波的情况下,翼型附近的流场也没有出现非物理震荡,人工黏性较好地发挥了作用。不过,在远离翼型的区域,解存在一些不连续的现象。这是由于这些地方网格比较稀疏,以及单元界面悬点(hanging point)的间断处理方式导致的。

在本小节的算例中,LC和BC延拓的计算迭代参数均使用表1和表2中的值。为保证迭代稳定,PTC方法只能使用比较保守的迭代参数,即η0=1,β=1.1,Rmax=1/1.5。

图7 NACA0012算例的压力系数云图

图8给出了NACA0012算例的非线性残差收敛历史。从图中可见,BC方法显示出比LC更高的效率。由于本算例PTC方法收敛太慢,图中并未将其收敛历史显示到完全收敛。并且,在该算例条件下,也无法找到合适的参数使得PTC方法能达到与拉普拉斯延拓方法相当的收敛效率。

图9给出了LC方法和BC方法延拓参数的变化历史。从图中可以看出,混合延拓方法在迭代的后期延拓参数都是由伪时间项主导的。尽管如此,在光滑条件下,BC方法仍能达到与LP方法相当的收敛效率。这在一定程度上说明了,在非线性迭代的初始阶段给予足够的耗散是加速收敛的关键。

在图8显示的残差收敛历史中,值得注意的是BC方法的收敛单调性虽然较差,但是残差的上扬没有超出PTC方法的范围。由此,可以认为,按照BC方法的混合思想,综合LC和PTC 2种方法的优点,构造出在光滑问题中的收敛效率不弱于LC方法,且在激波问题中稳定性不弱于PTC方法的非线性迭代策略。

图8 NACA0012算例通量残差收敛历史 图9 NACA0012算例延拓系数的变化历史 图10 NACA0012算例翼型压力系数分布的对比

图10给出了本文计算得到的翼型压力分布与Vassberg等人[18]通过4 096×4 096规模的网格得到的压力分布的对比。从图中可见,本文计算结果与参考结果符合较好,并且激波压力无振荡。说明本文工作在提高效率的同时保证了计算结果的精度。

4 结 论

航空气动定常数值模拟中初始解通常为来流状态的均匀场,初始残差集中在壁面附近。针对这一特点,对伪时间延拓法和拉普拉斯延拓法进行了分析,得到了在迭代初始阶段伪时间延拓法存在耗散不足问题的结论。对于解中包含激波的情况,单纯的拉普拉斯延拓法存在过度耗散的问题,对非线性迭代收敛的效率有不利影响。综合上述两点,构造了一种结合伪时间延拓法和拉普拉斯延拓法优点的混合延拓方法。在混合延拓方法中,时间项和拉普拉斯项通过线性权重进行组合。在非线性迭代的初期,混合延拓项由拉普拉斯项主导。在组合权重的计算中,以人工黏性系数为作为解光滑程度的指标,当解中出现激波结构时,混合延拓项将偏重于伪时间延拓法。

通过GAMM鼓包和NACA0012 2个算例对混合延拓方法进行了测试。数值结果表明,在亚声速的光滑问题中,混合方法可以自动退化为拉普拉斯延拓法,并保持其计算效率;考虑到混合延拓法在后期基本由伪时间推进法主导,这说明非线性迭代初期的耗散是加速收敛的关键。在存在激波的跨声速问题中,混合延拓方法可以取得比拉普拉斯延拓法更高的效率。在本文算例中,由于伪时间延拓法耗散不足,且数值求解中应用的直接求解器不受矩阵对角占优的约束,故伪时间延拓法表现出的计算效率最低。

本文工作表明针对定常流动问题,通过在非线性迭代的初期引入额外耗散来提高计算效率是可行的。在更复杂的问题中,需要引入迭代法求解线性问题,这种情况下拉普拉斯项对整体计算效率的影响还需进一步研究。

致谢本文部分思路和工作是第一作者在美国Texas A&M 大学数学系联合培养期间形成的。在此作者特别感谢国家留学基金委的资助,以及联合培养期间Wolfgang Bangerth教授在计算数学理论和方法以及deal.II有限元库的使用等方面给予的指导和帮助。

[1] 贡伊明, 张伟伟, 刘溢浪. 非定常求解的内迭代初值对计算效率的影响研究[J]. 西北工业大学学报, 2016, 34(1):11-17

Gong Yiming, Zhang Weiwei, Liu Yilang. Reaserching How Initial Value if Internal Iteration Impacts on Computational Efficieny in Unsteady Flow Solving[J]. Journal of Northwestern Polytechnical University, 2016, 34(1):11-17 (in Chinese)

[2] Geuzaine P. Newton-Krylov Strategy for Compressible Turbulent Flows on Unstructured Meshes[J]. AIAA Journal, 2001, 39(3): 528-531

[3] Dennis J, Schnabel R. Numerical Methods for Unconstrained Optimization and Nonlinear Equations[M]. Society for Industrial and Applied Mathematics, 1996

[4] Brown D A, Zingg D W. Advances in Homotopy Continuation Methods in Computational Fluid Dynamics[R]. AIAA-2013-2944

[5] Lyra P R M, Morgan K. A Review and Comparative Study of Upwind Biased Schemes for Compressible Flow Computation. PartⅢ: Multidimensional Extension on Unstructured Grids[J]. Archives of Computational Methods in Engineering, 2002, 9(3): 207-256

[6] Kelley C T, Liao L Z, Qi L, et al. Projected Pseudotransient Continuation[J]. SIAM Journal on Numerical Analysis, 2008, 46(6): 3071-3083

[7] Young D P, Melvin R G, Bieterman M B, et al. Global Convergence of Inexact Newton Methods for Transonic Flow[J]. International Journal for Numerical Methods in Fluids, 1990, 11(8): 1075-1095

[8] Hicken J, Buckley H, Osusky M, et al. Dissipation-Based Continuation: a Globalization for Inexact-Newton Solvers[R]. AIAA-2011-3237

[9] Pollock S. A Regularized Newton-Like Method for Nonlinear PDE[J]. Numerical Functional Analysis and Optimization, 2015, 36(11): 1493-1511

[10] Persson P O, Peraire J. Sub-Cell Shock Capturing for Discontinuous Galerkin Methods[R]. AIAA-2006-0112

[11] Kuzmin D, Moller M, Gurris M. Flux-Corrected Transport: Principles, Algorithms,and Applications[M]. 2nd Ed. Springer, 2012: 193-238

[12] Balay S, Abhyankar S, Adams M F, et al. PETSc Users Manual[R]. ANL-95/11 - Revision 3.6, 2015

[13] Amestoy P R, Duff I S, L′Excellent J-Y, et al. A Fully Asynchronous Multi-Frontal Solver Using Distributed Dynamic Scheduling[J]. SIAM Journal on Matrix Analysis and Applications, 2001, 23(1):15-41

[14] Bangerth W, Hartmann R, Kanschat G. DealⅡ——a General Purpose Object Oriented Finite Element Library[J]. ACM Trans on Mathematical Software, 2007, 33(4): 1-27

[15] Delchini M O, Ragusa J C, Berry R A. Entropy-Based Viscous Regularization for the Multi-Dimensional Euler Equations in Low-Mach and Transonic Flows[J]. Computers & Fluids, 2015, 118 : 225-244

[16] Žaloudek M, FoǐJ, Fürst J. Numerical Solution of Compressible Flow in a Channel and Blade Cascade[J]. Flow, Turbulence and Combustion, 2006, 76(4) : 353-361

[17] Janda M, Kozel K, Liska R. Hyperbolic Problem: Theory, Numnerics, Applications[C]∥8th International Conference on Magdeburg, 2001: 563-572

[18] Vassberg J, Jameson A. In Pursuit of Grid Convergence, Part I: Two Dimensional Euler Solutions[R]. AIAA-2009-4114