《非结构网格通量重构算法下三种紧致WENO限制器对比研究》由会员分享，可在线阅读，更多相关《非结构网格通量重构算法下三种紧致WENO限制器对比研究（18页珍藏版）》请在金锄头文库上搜索。
1、 非结构网格通量重构算法下三种紧致WENO限制器对比研究 石京昶，严 红，*(1. 西北工业大学 长三角研究院，太仓 215400；2. 陕西省航空发动机内流动力学重点实验室，西安 710072)0 引 言高精度数值格式因日益强烈的工程需求被关注和研究。其中代表性的方法有加权本质无振荡（weighted essentially non-oscillatory, WENO）、间断伽辽金（discontinuous Galerkin methods, DG）、通量重构（flux reconstruction, FR）等。WENO 方法通过非线性重构保证了激波附近“本质无振荡”。DG 和FR 方法天
2、然适用于非结构网格处理复杂几何外形，然而DG、FR 方法本身不能处理激波，所以需要引入限制器。DG、FR 方法框架下激波捕捉方法的一般思路是：首先采用间断探测器找到需要进行重构的网格单元，然后对目标网格单元采用限制器重构出新的解多项式。如何检测间断有不同的思路，参见Qiu 等的对比研究1。其中KXRCF 间断检测器2被广泛用于高阶DG 格式的激波捕捉方法研究中。有限体积方法下发展出的TVD、TVB 限制器可以应用于高阶DG、FR 方法，但会丢失高阶方法的高阶特性。WENO 的思想可以用来构造适用于DG 和FR 方法的限制器。Qiu 等3和Luo 等4提出了Hermite WENO 限制器，但是
3、并不紧致，目标单元解的重构需要用到相邻网格的相邻网格中的解。这种非紧致特性造成DG、FR 方法并行扩展性损失。Zhong 等5首先提出了紧致WENO 限制器用于DG 方法，核心思想是将相邻网格单元中的解插值到目标单元，并将单元平均替换为目标单元的单元平均，形成的新多项式与目标单元原本的多项式以WENO 的方式加权平均得到最终的重构多项式解。此WENO 限制器思路清晰，实现简单，故称为简单WENO 限制器（简称SWENO 限制器）。Zhu 等6-7在此基础上加以改进，将简单WENO 限制器中目标单元的相邻单元中的解直接插值到目标单元的方式改进为相邻单元中的解通过最小二乘法投影到目标单元。最近，Z
4、hu 等8-9又进一步提出了新的多精度WENO 限制器（简称MWENO限制器），其相对于改进的简单WENO 限制器做出的核心改进是，目标单元中的原始高阶多项式分解成从一阶至高阶的多个不同阶次多项式，相邻单元中的解多项式经最小二乘法投影到目标单元，形成一阶多项式后，用于计算最低阶多项式的光滑因子，最终多个不同阶次多项式仍采用WENO 方式加权重构得到最终的解多项式。Li 等10也于最近提出了p阶加权WENO 限制器（简称PWENO 限制器）。核心思想与Zhu 等最新的多精度WENO 限制器类似，不过相邻网格单元中的解经最小二乘法投影到目标单元，形成一阶多项式，参与WENO 重构得到最终的解多项式
5、，而不是只用于计算最低阶多项式的光滑因子。以上WENO 限制器在极端低密度低压的局部区域会重构出包含负密度负压的非物理解，因此需要使用保正限制器避免此情形出现。Zhang 等11提出的保正限制器能够保持原始高阶精度无损失，被广泛应用于高阶DG 格式激波捕捉方法研究中。但值得一提的是，此保正保精度限制器证明过程中假设单元解点必须是Gauss-Lobatto 点，相比于Gauss-Legendre 点精度较低。高阶格式计算含有激波的稳态问题时较难收敛到机器零，WENO 格式存在此问题。Zhang 等12提出迎风型重构减弱激波后伪振荡。Zhu 等13构造了新的多精度WENO 格式。DG 方法结合WE
6、NO 限制器激波捕捉同样难以收敛到机器零。最近，Zhu 等14提出了一种新的间断探测器，结合其多精度WENO限制器能够在标准算例中获得良好的收敛特性。该间断探测器核心思想是将KXRCF 中基于单元面上解计算差值改为基于单元内解计算差值，本文称之为CKXRCF 间断探测器。本文在FR 框架下，结合间断探测器和保正保精度限制器，对比研究简单WENO 限制器、p阶加权WENO 限制器和多精度WENO 限制器在多个经典算例中的性能，并讨论其高精度的优势和收敛性问题。常见文献中WENO 限制器验证测试算例集中在无黏Euler 方程，本文将其扩展到激波与层流边界层相互干扰的算例中，针对N-S 方程测试以上
7、WENO 限制器的性能，为相关研究提供初步结果。1 FR 方法简介非定常可压缩Navier-Stokes 方程的微分形式可写成如下双曲型守恒律的形式：式中，通量F(Q)=Fc(Q)Fv(Q)，守恒变量Q、无黏通量和黏性通量分别如下所示：本文简介FR 方法二维情形下的基本思想，详细阐述见Wang 等15和Huynh 等16的综述。计算域离散成N个网格单元ViiN=1。方程（1）的加权弱形式可改写为以下形式：式中，Qi是网格单元Vi上的近似解，且属于k阶多项式空间，即QiPk(Vi)。F(Qi)是F(Qi)在多项式空间Pk上的投影。iPk(Vi)是Vi上的修正项，由式（3）定义：式中，W是权函数，
8、Fn=FcnomFn(Qi)是网格界面上黎曼通量与单元内部局部通量的差值。将自由度定义在网格单元内部的解点。为计算修正项，在单元界面处定义通量点，则解点上的修正项由式（4）计算：式中，j,f,l是独立于解的常数，Sf是单元界面面积，j是解点在当前网格单元内的编号，f是当前网格单元面的编号，l是通量点在当前面的编号。综上所述，FR 方法由以下离散方程近似求解原双曲律方程：本文所有算例均使用自主开发的N-S 方程非结构网格并行求解器NFR，时间推进方法均采用三阶强稳定性显式SSP-RK3 方法，无黏通量采用Roe 通量或LaxFriedrichs 通量，黏性通量均采用BR2 格式17。本文所有算例
9、的结果均将原始网格单元多项式阶次均等剖分为子单元，如一个四边形网格单元的P2 剖分为4 个四边形子单元。2 WENO 限制器简介本文采用的WENO 限制器的基本思路是使用间断探测器找出需要限制解的问题网格单元，然后对其应用WENO 限制器重构得到新解。对于欧拉方程系统，为了更好地避免激波附近的振荡，使用通量雅可比矩阵将守恒变量转换为特征变量，经WENO 限制器重构得到新的特征变量，再经通量雅可比矩阵转换回守恒变量。下述3 种WENO 限制器均只简介其针对一维标量方程的限制过程。2.1 简单WENO 限制器简单WENO 限制器的核心是将问题网格单元Ij及其相邻网格单元的解多项式加权组合得到重构解
10、，保证其与原始解有相同的单元平均和k阶精度。权重由解多项式的光滑因子确定。记网格单元Ij1、Ij、Ij+1的解为pL(x)、p0(x)、pR(x)，简单WENO 限制器重构的新解为：式中，L、R是归一化的非线性权重，(x)、(x)是相邻网格单元内解多项式pL(x)、pR(x)投影到Ij中的新多项式。为保证守恒性，相邻网格单元内的解向Ij中的投影形式为：式中，单元平均定义为：注意以上积分均为目标单元Ij内的积分。上述归一化非线性权重l计 算式如下：式中，其中，常数=1106，r=2。线性权重L=0.001，0=0.998,R=0.001。光滑因子l采用经典的WENO方式，定义如下：2.2 p 阶
11、加权WENO 限制器记网格单元Ij1、Ij、Ij+1的解为pL(x)、p0(x)、pR(x)。重构问题网格单元为Ij，其内各阶解多项式p0,s可通过将原始高阶解多项式p0转换为模态多项式获得。为保证守恒性且将相邻网格单元内的线性解多项式投影到问题网格单元，在问题网格单元上可定义来自相邻网格单元Ij1的线性多项式如下：来自相邻网格单元Ij+1的线性多项式可通过相似方式得到。最终的WENO 重构解多项式定义为目标单元内的解与相邻网格单元解的加权组合，保证了和原始解p0有相同的单元平均和k阶精度。上述归一化非线性权重l计算同公式（9），式中l计算式如下：线性权重l和参数l与简单WENO 限制器中的不
12、同，不再是常数，而是因为l对应的阶次不同而有区别，设定如下：式中，参数K取0.01。显然，线性权重l和参数l对 阶次的依赖关系弱化了相邻网格单元Ij1内的解投影到目标网格单元Ij中的线性多项式，强化了Ij中的高阶模态。另外，当低阶模态的振荡程度比高阶模态更严重时，令低阶模态的线性权重为0，去掉振荡程度更严重的部分模态。按照如下平均方式比较两者的光滑因子：式中，参数Ktrunc常取1，s,r是多项式p0,s的光滑因子s的分量，即因此，计算光滑因子s时存储s,r留待此处使用。方程（19）中右端高阶模态的光滑因子平均值对每个s阶模态是变化的。令s从最高阶k降序至2，对每个s检查；如果低阶模态的光滑因
13、子平均值更大，则令高阶模态之外的模态多项式对应的线性权重为0。2.3 多精度WENO 限制器记网格单元Ij1、Ij、Ij+1的解为pL(x)、p0(x)、pR(x)。重构问题网格单元为Ij，其内各阶解多项式p0,s可通过将原始高阶解多项式p0转换为模态多项式获得。基于各阶解多项式重构得到新的线性解多项式p和非线性解多项式p。上述非线性权重l1,l2计算如下：式中，光滑因子0,1的计算并不依赖于零阶多项式p0,0，而是依赖相邻网格单元中重构得到的多项式：式中，L、R为相邻网格单元重构解的光滑因子。为保证守恒性，且将相邻网格单元内的线性解多项式投影到目标单元，在目标单元上可定义来自相邻网格单元Ij
14、1的线性多项式如下：来自相邻网格单元Ij+1的线性多项式用相似方式得到。最终的WENO 重构解多项式定义为最高k阶非线性解多项式，保证了和原始解p0有相同的单元平均和k阶精度。2.4 保正保精度限制器尽管前述限制器抑制了FR 方法在激波附近的振荡，但无法完全避免出现负密度和压力，限制器对解多项式的重构并不考虑是否会重构出非物理的解。因此还需要保正保精度限制器配合激波捕捉方法。本文实现的保正保精度限制器11要求使用Legendre-Gauss-Lobatto 点，以保证采用显式RK 时间推进格式和一定的CFL（Courant-Friedrichs-Lewy）数下，可以从数学上证明该限制器保正保精
15、度且守恒。记目标网格单元为Ij，k阶多项式解分布在N(k)个Legendre-Gauss-Lobatto 点上，记守恒标量qj,i、密度 j,i、压力pj,i，i1,N(k)，=11015表征极小的正值。保正保精度限制器的算法流程如下：1）限制密度。找出当前网格单元所有解点上密度 j,i的最小值，若该值小于允许的最小正值，则将该网格单元中所有解点上的密度 j,i等比例缩小。本文所有算例均使用了保正保精度限制器。3 间断探测器简介3.1 KXRCF 间断探测器记网格单元Ij1、Ij、Ij+1的解为pL(x)、p0(x)、pR(x)。将网格单元Ij的所有面分为两组：流入面组Ij和流出面组I+j，则可由如下条件判断间断：其中，m=1，Ck=1，hj是网格单元Ij的特征半径。公式（33）的物理含义是，如果该网格单元与相邻网格单元在“流入面组”上的间断差值占比超过了Ck值，就判定其为问题网格单元。3.2 CKXRCF 间断探测器记网格单元Ij1、Ij、Ij+1的解为pL(x)、p0(x)、pR(x)。先计算目标网格单元及其相邻网格单元内解的积分，再由其差值判断间断。具体如下：其中，Ck=1，hj是网格单元Ij的特征半径。CKXRCF间断探测器与KXRCF 间断探测器的不同在于前者使用了单元内解的积分，而非网格单元交界面通量点处解的积分。一般来说，流场中存
