五阶WENO格式的基本原理及二维情况的代码

Detonation6

已于 2024-12-11 15:40:35 修改

阅读量3.5k

点赞数

分类专栏：计算流体力学里的若干问题介绍文章标签：算法

于 2023-07-31 20:01:16 首次发布

本文链接：https://blog.csdn.net/qq_35152062/article/details/131917946

版权

计算流体力学里的若干问题介绍专栏收录该内容

1 篇文章

订阅专栏

该文详细介绍了五阶WENO格式的基本原理，包括ENO重构模板的凸组合，权重因子的计算以及二维情况下求解界面两侧的源代码实现。通过权重优化，该方法旨在减少数值震荡，提高数值模拟的精度。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

1. 五阶WENO格式的基本原理

WENO的全称为加权本质无振荡方法，是ENO方法的衍生形态。通过将所有的ENO重构模板进行凸组合，获得较为逼近的参数空间分布。假设 $k$ 个模板的表达式如下所示，

$S_{r}(i) = \begin{Bmatrix} x_{i-r} & x_{i-r+1} & ... & x_{i-r+k-1} \end{Bmatrix}$

其中 $r$ 的取值如下所示，

$r = 0, 1, ... , k-1.$

那么可以获得 ${Q}_{i+1/2}^L$ 的 $k$ 种重构方式，具体数学表达式如下所示，

$Q_{i+1/2}^{(r)}=\sum_{j=0}^{k-1}Const_{rj}\cdot Q_{i-r+j}$

可能使用到的格心点参数组合有如下三种形式，

$S_0=\begin{Bmatrix} Q_{i-2} & Q_{i-1} & Q_{i} \end{Bmatrix}$

$S_1=\begin{Bmatrix} Q_{i-1} & Q_{i} & Q_{i+1} \end{Bmatrix}$

$S_2=\begin{Bmatrix} Q_{i} & Q_{i+1} &Q_{i+2} \end{Bmatrix}$

进一步地， $U_{i+1/2}^{(r)}$ 的三种拟合情况如下所示，

$\left\{\begin{matrix} Q_{i+1/2}^{(0)}= \frac{1}{3} Q_{i-2} -\frac{7}{6} Q_{i-1}+\frac{11}{6}Q_{i}\\ Q_{i+1/2}^{(1)}=-\frac{1}{6} Q_{i-1}+\frac{5}{6} Q_{i}+\frac{1}{3}Q_{i+1}\\ Q_{i+1/2}^{(2)}=\frac{1}{3} Q_{i}+\frac{5}{6} Q_{i+1}-\frac{1}{6}Q_{i+2}\\ \end{matrix}\right.$

加权重构公式如下所示，

$Q_{i+1/2}^L= w_{0}Q_{i+1/2}^{(0)}+w_{1}Q_{i+1/2}^{(1)}+w_{2}Q_{i+1/2}^{(2)}$

权重因子 $w_{k}$ 的表达式如下所示，

$\left\{\begin{matrix} w_{k} = \frac{\alpha _k}{\sum_{i=0}^{2}\alpha _k}\\ a_k=\frac{d_k}{(\beta _k+\varepsilon)^2} \end{matrix}\right.$

其中权系数 $d_k$ 有如下表达式，

$d_0=0.3, d_1=0.6, d_2=0.1$

参数 $\varepsilon$ 的作用是避免分母为零，

$\varepsilon = 1\times 10^{-6}$

此外光滑度度量系数 $\beta _k$ 的表达式为，

$\left\{\begin{matrix} \beta _2 = \frac{13}{12} (Q_i-2Q_{i+1}+Q_{i+2})^2+\frac{1}{4}(3Q_i-4Q_{i+1}+Q_{i+2})^2\\ \beta _1 =\frac{13}{12}(Q_{i-1}-2Q_i+Q_{i+1})^2+\frac{1}{4}(Q_{i-1}-Q_{i+1})^2\\ \beta_0 =\frac{13}{12}(Q_{i-2}-2Q_{i-1}+Q_i)^2+\frac{1}{4}(Q_{i-2}-4Q_{i+1}+3Q_{i})^2 \end{matrix}\right.$

二维情况下，用于求解界面+和-两侧的WENO格式m源代码如下所示，

function [Fs,Gs]=WENO_Scheme2D(Fp,Fn,Gp,Gn)
%统计总数
[~,N1,N2]=size(Fp);
%五点模板
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Fsp
%Fp
Fp0= 1/3*Fp(:,1:N1-5,4:N2-3)...,
    -7/6*Fp(:,2:N1-4,4:N2-3)...,
   +11/6*Fp(:,3:N1-3,4:N2-3);
Fp1=-1/6*Fp(:,2:N1-4,4:N2-3)...,
    +5/6*Fp(:,3:N1-3,4:N2-3)...,
    +1/3*Fp(:,4:N1-2,4:N2-3);
Fp2= 1/3*Fp(:,3:N1-3,4:N2-3)...,
    +5/6*Fp(:,4:N1-2,4:N2-3)...,
    -1/6*Fp(:,5:N1-1,4:N2-3);
%bFp
bFp0=13/12*(Fp(:,1:N1-5,4:N2-3)-2*Fp(:,2:N1-4,4:N2-3)+1*Fp(:,3:N1-3,4:N2-3)).^2+...,
       1/4*(Fp(:,1:N1-5,4:N2-3)-4*Fp(:,2:N1-4,4:N2-3)+3*Fp(:,3:N1-3,4:N2-3)).^2;
bFp1=13/12*(Fp(:,2:N1-4,4:N2-3)-2*Fp(:,3:N1-3,4:N2-3)+1*Fp(:,4:N1-2,4:N2-3)).^2+...,
       1/4*(Fp(:,2:N1-4,4:N2-3)-1*Fp(:,4:N1-2,4:N2-3)).^2;
bFp2=13/12*(Fp(:,3:N1-3,4:N2-3)-2*Fp(:,4:N1-2,4:N2-3)+1*Fp(:,5:N1-1,4:N2-3)).^2+...,
     1/4*(3*Fp(:,3:N1-3,4:N2-3)-4*Fp(:,4:N1-2,4:N2-3)+1*Fp(:,5:N1-1,4:N2-3)).^2;
%aFp
d0=1/10;
d1=3/5;
d2=3/10;
aFp0=d0./(eps+bFp0).^2;
aFp1=d1./(eps+bFp1).^2;
aFp2=d2./(eps+bFp2).^2;
%wFp
wFp0=aFp0./(aFp0+aFp1+aFp2);
wFp1=aFp1./(aFp0+aFp1+aFp2);
wFp2=aFp2./(aFp0+aFp1+aFp2);
%Fsp
Fsp=wFp0.*Fp0+wFp1.*Fp1+wFp2.*Fp2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Fsn            
%Fn
Fn0= 1/3*Fn(:,6:N1,4:N2-3)...,
    -7/6*Fn(:,5:N1-1,4:N2-3)...,
   +11/6*Fn(:,4:N1-2,4:N2-3);
Fn1=-1/6*Fn(:,5:N1-1,4:N2-3)...,
    +5/6*Fn(:,4:N1-2,4:N2-3)...,
    +1/3*Fn(:,3:N1-3,4:N2-3);
Fn2=1/3*Fn(:,4:N1-2,4:N2-3)...,
   +5/6*Fn(:,3:N1-3,4:N2-3)...,
   -1/6*Fn(:,2:N1-4,4:N2-3);
%bFn
bFn0=13/12*(Fn(:,6:N1,4:N2-3)-2*Fn(:,5:N1-1,4:N2-3)+1*Fn(:,4:N1-2,4:N2-3)).^2+...,
       1/4*(Fn(:,6:N1,4:N2-3)-4*Fn(:,5:N1-1,4:N2-3)+3*Fn(:,4:N1-2,4:N2-3)).^2;
bFn1=13/12*(Fn(:,5:N1-1,4:N2-3)-2*Fn(:,4:N1-2,4:N2-3)+1*Fn(:,3:N1-3,4:N2-3)).^2+...,
       1/4*(Fn(:,5:N1-1,4:N2-3)-1*Fn(:,3:N1-3,4:N2-3)).^2;
bFn2=13/12*(Fn(:,4:N1-2,4:N2-3)-2*Fn(:,3:N1-3,4:N2-3)+1*Fn(:,2:N1-4,4:N2-3)).^2+...,
     1/4*(3*Fn(:,4:N1-2,4:N2-3)-4*Fn(:,3:N1-3,4:N2-3)+1*Fn(:,2:N1-4,4:N2-3)).^2;
%aFn
aFn0=d0./(eps+bFn0).^2;
aFn1=d1./(eps+bFn1).^2;
aFn2=d2./(eps+bFn2).^2;
%wFn
wFn0=aFn0./(aFn0+aFn1+aFn2);
wFn1=aFn1./(aFn0+aFn1+aFn2);
wFn2=aFn2./(aFn0+aFn1+aFn2);
%Fsn
Fsn=wFn0.*Fn0+wFn1.*Fn1+wFn2.*Fn2;
Fs=Fsp+Fsn; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gsp
%Gp
Gp0= 1/3*Gp(:,4:N1-3,1:N2-5)...,
    -7/6*Gp(:,4:N1-3,2:N2-4)...,
   +11/6*Gp(:,4:N1-3,3:N2-3);
Gp1=-1/6*Gp(:,4:N1-3,2:N2-4)...,
    +5/6*Gp(:,4:N1-3,3:N2-3)...,
    +1/3*Gp(:,4:N1-3,4:N2-2);
Gp2= 1/3*Gp(:,4:N1-3,3:N2-3)...,
    +5/6*Gp(:,4:N1-3,4:N2-2)...,
    -1/6*Gp(:,4:N1-3,5:N2-1);
%bGp
bGp0=13/12*(Gp(:,4:N1-3,1:N2-5)-2*Gp(:,4:N1-3,2:N2-4)+1*Gp(:,4:N1-3,3:N2-3)).^2+...,
       1/4*(Gp(:,4:N1-3,1:N2-5)-4*Gp(:,4:N1-3,2:N2-4)+3*Gp(:,4:N1-3,3:N2-3)).^2;
bGp1=13/12*(Gp(:,4:N1-3,2:N2-4)-2*Gp(:,4:N1-3,3:N2-3)+1*Gp(:,4:N1-3,4:N2-2)).^2+...,
       1/4*(Gp(:,4:N1-3,2:N2-4)-1*Gp(:,4:N1-3,4:N2-2)).^2;
bGp2=13/12*(Gp(:,4:N1-3,3:N2-3)-2*Gp(:,4:N1-3,4:N2-2)+1*Gp(:,4:N1-3,5:N2-1)).^2+...,
     1/4*(3*Gp(:,4:N1-3,3:N2-3)-4*Gp(:,4:N1-3,4:N2-2)+1*Gp(:,4:N1-3,5:N2-1)).^2;
%aGp
aGp0=d0./(eps+bGp0).^2;
aGp1=d1./(eps+bGp1).^2;
aGp2=d2./(eps+bGp2).^2;
%wGp
wGp0=aGp0./(aGp0+aGp1+aGp2);
wGp1=aGp1./(aGp0+aGp1+aGp2);
wGp2=aGp2./(aGp0+aGp1+aGp2);
%Gsp
Gsp=wGp0.*Gp0+wGp1.*Gp1+wGp2.*Gp2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gsn
%Gn
Gn0=1/3*Gn(:,4:N1-3,6:N2)...,
    -7/6*Gn(:,4:N1-3,5:N2-1)...,
   +11/6*Gn(:,4:N1-3,4:N2-2);
Gn1=-1/6*Gn(:,4:N1-3,5:N2-1)...,
    +5/6*Gn(:,4:N1-3,4:N2-2)...,
    +1/3*Gn(:,4:N1-3,3:N2-3);
Gn2=1/3*Gn(:,4:N1-3,4:N2-2)...,
    +5/6*Gn(:,4:N1-3,3:N2-3)...,
    -1/6*Gn(:,4:N1-3,2:N2-4);
%bGn
bGn0=13/12*(Gn(:,4:N1-3,6:N2)-2*Gn(:,4:N1-3,5:N2-1)+1*Gn(:,4:N1-3,4:N2-2)).^2+...,
       1/4*(Gn(:,4:N1-3,6:N2)-4*Gn(:,4:N1-3,5:N2-1)+3*Gn(:,4:N1-3,4:N2-2)).^2;
bGn1=13/12*(Gn(:,4:N1-3,5:N2-1)-2*Gn(:,4:N1-3,4:N2-2)+1*Gn(:,4:N1-3,3:N2-3)).^2+...,
       1/4*(Gn(:,4:N1-3,5:N2-1)-1*Gn(:,4:N1-3,3:N2-3)).^2;
bGn2=13/12*(Gn(:,4:N1-3,4:N2-2)-2*Gn(:,4:N1-3,3:N2-3)+1*Gn(:,4:N1-3,2:N2-4)).^2+...,
     1/4*(3*Gn(:,4:N1-3,4:N2-2)-4*Gn(:,4:N1-3,3:N2-3)+1*Gn(:,4:N1-3,2:N2-4)).^2;
%aGn
aGn0=d0./(eps+bGn0).^2;
aGn1=d1./(eps+bGn1).^2;
aGn2=d2./(eps+bGn2).^2;
%wGn
wGn0=aGn0./(aGn0+aGn1+aGn2);
wGn1=aGn1./(aGn0+aGn1+aGn2);
wGn2=aGn2./(aGn0+aGn1+aGn2);
%Gsn
Gsn=wGn0.*Gn0+wGn1.*Gn1+wGn2.*Gn2;
Gs=Gsp+Gsn;
end