2019.04.16--紫金VS玉泉.-应用二阶TVD格式和五阶WENO格式计算一维无粘激波管流场

课程名称：计算空气动力学姓名：紫金VS玉泉指导老师：Stone专业：飞行器设计与工程日期：2019年4月16日应用二阶TVD格式和五阶WENO格式计算一维无粘激波管流场摘要：控制方程为一维无粘欧拉方程，二阶TVD格式采用空间VanLeer二阶迎风格式结合vanleer限制器，时间一阶前向差分，求解一维无粘激波管流场的二阶数值解。另外，应用五阶WENO格式求解一维无粘激波管流场的高阶数值解。①出发方程为：𝜕𝑈𝜕𝑡+𝜕𝐹𝜕𝑥=0其中：𝑈=[𝜌𝜌𝑢𝐸]，𝐹=[𝜌𝑢𝜌𝑢2+𝑝(𝐸+𝑝)𝑢]，E=𝑝𝛾−1+12𝜌𝑢2，γ=1.4②初始条件是：𝑈(𝑥,0)={𝑈𝐿𝑥0𝑈𝑅𝑥0其中：𝑈𝐿=(102.5)，𝑈𝑅=(0.12500.25)另外，时间步长∆𝑡=0.0001，空间步长∆𝑥=0.005③二阶TVD格式的具体形式如下：𝑀𝑎≤−1时，𝐹+=0，𝐹−=𝐹𝑀𝑎≥1时，𝐹+=𝐹，𝐹−=0|𝑀𝑎|1时，𝐹±=[𝑓1±𝑓1±𝛾[(𝛾−1)𝑢±2𝑐]𝑓1±2(𝛾2−1)[(𝛾−1)𝑢±2𝑐]2]其中：𝑓1±=±𝜌𝑐[𝑀𝑎±12]2二阶迎风格式的空间差分算式为：𝜕𝐹𝜕𝑥=3𝑓𝑖𝑛+−4𝑓𝑖−1𝑛++𝑓𝑖−2𝑛+−𝑓𝑖+2𝑛−+4𝑓𝑖+1𝑛−−3𝑓𝑖𝑛−2∆𝑥，一阶精度的时间差分算式为：𝜕𝑈𝜕𝑡=𝑈𝑖+1−𝑈𝑖∆𝑡，Vanleer限制器：𝛹(𝑟)=𝑟+|𝑟|1+𝑟，𝑟𝑖=𝑈𝑖−𝑈𝑖−1𝑈𝑖+1−𝑈𝑖最终的计算格式为：𝑈𝑖𝑛+1=𝑈𝑖𝑛−∆𝑡∆𝑥(𝑓𝑖𝑛+−𝑓𝑖−1𝑛++𝑓𝑖+1𝑛−−𝑓𝑖𝑛−)−𝛹(𝑟𝑖)[∆𝑡2∆𝑥(3𝑓𝑖𝑛+−4𝑓𝑖−1𝑛++𝑓𝑖−2𝑛+−𝑓𝑖+2𝑛−+4𝑓𝑖+1𝑛−−3𝑓𝑖𝑛−)−∆𝑡∆𝑥(𝑓𝑖𝑛+−𝑓𝑖−1𝑛++𝑓𝑖+1𝑛−−𝑓𝑖𝑛−)]④五阶WENO格式的具体形式如下：𝑈𝑖𝑛+1=𝑈𝑖𝑛−∆𝑡∆𝑥(𝐹𝑗+12𝑊𝐸𝑁𝑂−𝐹𝑗−12𝑊𝐸𝑁𝑂)正通量的情况：𝐹𝑗+12𝑊𝐸𝑁𝑂=𝜔1+𝐹𝑗+12(1)++𝜔2+𝐹𝑗+12(2)++𝜔3+𝐹𝑗+12(3)+𝜔𝑘+=𝛼𝑘+𝛼1++𝛼2++𝛼3+，𝑘=1,2,3𝛼𝑘+=𝐶𝑘(𝜀+𝐼𝑆𝑘+)𝑝，𝑘=1,2,3，p=2，𝜀=10−6，𝐶1=110，𝐶1=610，𝐶1=310𝐼𝑆1+=14(𝐹𝑗−2+−4𝐹𝑗−1++3𝐹𝑗+)2+1312(𝐹𝑗−2+−2𝐹𝑗−1++𝐹𝑗+)2𝐼𝑆2+=14(𝐹𝑗−1+−𝐹𝑗+1+)2+1312(𝐹𝑗−1+−2𝐹𝑗++𝐹𝑗+1+)2𝐼𝑆3+=14(3𝐹𝑗+−4𝐹𝑗+1++𝐹𝑗+2+)2+1312(𝐹𝑗+−2𝐹𝑗+1++𝐹𝑗+2+)2𝐹𝑗+12(1)+=13𝐹𝑗−2+−76𝐹𝑗−1++116𝐹𝑗+𝐹𝑗+12(2)+=−16𝐹𝑗−1++56𝐹𝑗++13𝐹𝑗+1+𝐹𝑗+12(3)+=13𝐹𝑗++56𝐹𝑗+1+−16𝐹𝑗+2+负通量的情况：𝐹𝑗−12𝑊𝐸𝑁𝑂=𝜔1−𝐹𝑗−12(1)−+𝜔2−𝐹𝑗−12(2)−+𝜔3−𝐹𝑗−12(3)−𝜔𝑘−=𝛼𝑘−𝛼1−+𝛼2−+𝛼3−，𝑘=1,2,3𝛼𝑘−=𝐶𝑘(𝜀+𝐼𝑆𝑘−)𝑝，𝑘=1,2,3，p=2，𝜀=10−6，𝐶1=110，𝐶1=610，𝐶1=310𝐼𝑆1−=14(𝐹𝑗+2−−4𝐹𝑗+1−+3𝐹𝑗−)2+1312(𝐹𝑗+2−−2𝐹𝑗+1−+𝐹𝑗−)2𝐼𝑆2−=14(𝐹𝑗+1−−𝐹𝑗−1−)2+1312(𝐹𝑗+1−−2𝐹𝑗−+𝐹𝑗−1−)2𝐼𝑆3−=14(3𝐹𝑗−−4𝐹𝑗−1−+𝐹𝑗−2−)2+1312(𝐹𝑗−−2𝐹𝑗−1−+𝐹𝑗−2−)2𝐹𝑗−12(1)−=13𝐹𝑗+2−−76𝐹𝑗+1−+116𝐹𝑗−𝐹𝑗−12(2)−=−16𝐹𝑗+1−+56𝐹𝑗−+13𝐹𝑗−1−𝐹𝑗−12(3)−=13𝐹𝑗−+56𝐹𝑗−1−−16𝐹𝑗−2−其中，正负通量𝐹±的表达式如下：𝑀𝑎≤−1时，𝐹+=0，𝐹−=𝐹𝑀𝑎≥1时，𝐹+=𝐹，𝐹−=0|𝑀𝑎|1时，𝐹±=[𝑓±𝑓±𝛾[(𝛾−1)𝑢±2𝑐]𝑓±2(𝛾2−1)[(𝛾−1)𝑢±2𝑐]2]，其中：𝑓±=±𝜌𝑐[𝑀𝑎±12]2⑤边界条件（N个空间网格点），应满足如下关系：𝑖=1时，𝑈𝑖−1𝑛=𝑈𝑖𝑛，𝑈𝑖−2𝑛=𝑈𝑖𝑛，𝐹𝑖−1𝑛=𝐹𝑖𝑛，𝐹𝑖−2𝑛=𝐹𝑖𝑛𝑖=2时，𝑈𝑖−2𝑛=𝑈𝑖𝑛，𝐹𝑖−2𝑛=𝐹𝑖𝑛𝑖=N时，𝑈𝑖+1𝑛=𝑈𝑖𝑛，𝑈𝑖+2𝑛=𝑈𝑖𝑛，𝐹𝑖+1𝑛=𝐹𝑖𝑛，𝐹𝑖+2𝑛=𝐹𝑖𝑛𝑖=N−1时，𝑈𝑖+2𝑛=𝑈𝑖𝑛，𝐹𝑖+2𝑛=𝐹𝑖𝑛本文中，N=200⑥图1(a)、图1(b)、图1(c)分别为t=0.2时刻，用二阶TVD格式计算出的密度、压力和速度分布及其同精确解的比较。图2(a)、图2(b)、图2(c)分别为t=0.2时刻，用五阶WENO格式计算出的密度、压力和速度分布及其同精确解的比较。图1可见，二阶TVD格式对于原二阶迎风格式在接触间断、激波以及膨胀波附近的非物理震荡有着良好的修正效果。但TVD格式受其构造方式所限，在解的局部极值点附近只能达到一阶精度。图2可见，五阶WENO格式不仅对于接触间断、激波以及膨胀波附近的非物理震荡有着良好的修正效果，而且对于接触间断、激波以及膨胀波的捕捉精度要比原二阶迎风格式精度要高。⑦源代码（以密度分布为例）：【基于VanLeer二阶TVD格式的一维激波管非定常流动数值解(t=0.2)】N=200;%网格节点数N_left=N./2;%Sod激波管膜片【左侧】的网格节点数N_right=N./2;%Sod激波管膜片【右侧】的网格节点数density=[ones(1,N_left)0.125*ones(1,N_right)];%各网格节点【t=0时刻密度】u=zeros(1,N);%各网格节点【t=0时刻速度】p=[ones(1,N_left)0.1*ones(1,N_right)];%各网格节点【t=0时刻压力】gamma=1.4;%比热比delta_t=0.0001;%时间步长delta_x=0.005;%空间步长lambda=delta_t./delta_x;U_1=zeros(1,N);%守恒量预赋值U_2=zeros(1,N);U_3=zeros(1,N);f_1_plus=zeros(1,N);%通量预赋值f_1_minus=zeros(1,N);f_2_plus=zeros(1,N);f_2_minus=zeros(1,N);f_3_plus=zeros(1,N);f_3_minus=zeros(1,N);Ma=zeros(1,N);%各网格节点处的马赫数预赋值e=zeros(1,N);%能量预赋值F_1_first=zeros(1,N);%【一阶精度】通量预赋值F_2_first=zeros(1,N);F_3_first=zeros(1,N);F_1_second=zeros(1,N);%【二阶精度】通量预赋值F_2_second=zeros(1,N);F_3_second=zeros(1,N);limiter_1=zeros(1,N);%【限制器】limiter_2=zeros(1,N);limiter_3=zeros(1,N);fort=0:0.0001:0.2forj=1:1:200%各网格节点处的马赫数赋值Ma(1,j)=u(1,j)./sqrt(gamma.*p(1,j)./density(1,j));endforj=1:1:200%能量赋初值e(1,j)=p(1,j)./(gamma-1)+0.5.*density(1,j).*u(1,j).*u(1,j);end%-----守恒量赋初值-----%forj=1:1:200U_1(1,j)=density(1,j);U_2(1,j)=density(1,j).*u(1,j);U_3(1,j)=p(1,j)./(gamma-1)+0.5.*density(1,j).*u(1,j).*u(1,j);endforj=1:1:200%-----第1个通量-----%ifMa(1,j)-1&&Ma(1,j)1%马赫数介于-1和1之间f_1_plus(1,j)=density(1,j).*sqrt(gamma.*p(1,j)./density(1,j)).*(0.5*(Ma(1,j)+1)).*(0.5*(Ma(1,j)+1));f_1_minus(1,j)=-density(1,j).*sqrt(gamma.*p(1,j)./density(1,j)).*(0.5*(Ma(1,j)-1)).*(0.5*(Ma(1,j)-1));elseifMa(1,j)=-1%马赫数=-1f_1_plus(1,j)=0;f_1_minus(1,j)=density(1,j).*u(1,j);elseifMa(1,j)=1%马赫数=1f_1_plus(1,j)=density(1,j).*u(1,j);f_1_minus(1,j)=0;end%-----第2个通量-----%ifMa(1,j)-1&&Ma(1,j)1%马赫数介于-1和1之间f_2_plus(1,j)=(density(1,j).*sqrt(gamma.*p(1,j)./density(1,j)).*(0.5*(Ma(1,j)+1)).*(0.5*(Ma(1,j)+1))).*((gamma-1).*u(1,j)+2.*sqrt(gamma.*p(1,j)./density(1,j)))./gamma;f_2_minus(1,j)=-density(1,j).*sqrt(gamma.*p(1,j)./density(1,j)).*(0.5*(Ma(1,j)-1)).*(0.5*(Ma(1,j)-1)).*((gamma-1).*u(1,j)-2.*sqrt(gamma.*p(1,j)./density(1,j)))./gamma;elseifMa(1,j)=-1%马赫数=-1f_2_plus(1,j)=0;f_2_minus(1,j)=density(1,j).*u(1,j).*u(1,j)+p(1,j);elseifMa(1,j)=1%马赫数=1f_2_plus(1,j)=density(1,j).*u(1,j).*u(1,j)+p(1,j);f_2_minus(1,j)=0;end%-----第3个通量-----%ifMa(1,j)-1&&Ma(1,j)1%马赫数介于-1和1之间f_3_plus(1,j)=(density(1,j).*sqrt(gamma.*p(1,j)./density(1,j)).*(0.5*(Ma(1,j)+1)).*(0.5*(Ma(1,j)+1))).*((gamma-1).*u(1,j)+2.*sqrt(gamma.*p(1,j)./density(1,j))).*((gamma-1).*u(1,j)+2.*sqrt(gamma.*p(1,j)./density(1,j)))./(2.