《中国科学院大学 计算流体力学 作业 6.1》由会员分享，可在线阅读，更多相关《中国科学院大学 计算流体力学 作业 6.1（11页珍藏版）》请在金锄头文库上搜索。
1、计算流体力学第五次作业6.1作业6.1针对如下Sod激波管问题红+也2=0 dt dx d(pu) g(pr?4-p)_0x0,1dt dx6(pE) , dEu + pu) _ q dt dxA. /、Jo,U X0 5川5阶WENO格式计算其数值解，加出t=0.14时刻密度、速度及压力的分布;并与粘 确解进行比拟(要求数值解与精确解画在同一张图上,便于比拟)。要求：空间网格数100,时间推进格式选用3阶Runge-Kutta,时间步长自选。本次作业包含两个局部：第一局部：第二次作业求sod激波管理论解：(详见第二次作业).主程序lyr_cfd_Sodshocktube (画图局部作了轻微改
2、动)1 .子程序(1) fcenter中心区压力函数(2) dfcenter中心区压力函数导函数(3) newton牛顿迭代法求解中心区压力px(4) drawpic画出t=0.14s时速度、压力、密度的分布曲线第二局部：求sod激波管数值解：本次作业基于流通矢量分裂方法，Steger-Warming分裂。方法如下PPT所示：wl=al/(al+a2+a3);w2=a2/ (al+a2+a3);w3=a3/(al+a2+a3);fweno(i,j)=wl*Ujl+w2*Uj 2+w3*Uj 3;end国处理边界条件句=1时取w1=w2=0fweno(i,1)=1/3*U(i,1)+5/6*U(
3、i,2)-1/6*U(i, 3);% j =2 时取 wl = 0IS2_2 = l/4*(U(i,1)-U(iA 3)A2 + 13/12*(U(i,1)-2*U(i,2)+U( i,37)A2;工S3_2 = 1/4*(3*U(i,2)-4*U(i,3)+U(i,4) )A2 + 13/12* (U(iA2)- 2夫U?i,3)+U(i,4)A2;a2_2=(6/10)/(10A(-6)+IS2_2)人2;a32=(3/10)/(10A(-6)+工S312)2w2_2=a2_2/(a2_2+a3_2);w3_2=a3_2/(a2_2+a3_2);fweno(i,2)=w2_2*(-1/6夫
4、U(i,1)+5/6*U(i,2)+1/3*U(i,3)+w 3_2夫(1/3夫U(i,5)+5/6*U(i,3)-l/6*U(i,4);%j=100 时取w3=0IS1_10 0 = 1/4*(U(i,98)-4*U(iz 99)+3*U(i,100)A2 + 13/12*(U (i,8)-2夫U(i,99)+U(i,100)八2;IS2_10 0 = 1/4*(U(iz99)-U(if101)A2 + 13/12*(U(iz 99)-2*U(i ,100)+U(i,101)A2;al_100=(1/10)/(10A(-6)+ISl_100)A2;a21100= (6/10) / (10人(
5、-6) + 工S2100) A2;wl2100=al_100/(al_100+a2_100);w2100=a2l00/(al100+a2100);fweno(i,100)=wl_100*(l/3*U(i,98)-7/6*U(i,99)+ll/6*U(i z100)+w2_100*(=1/6*U(i,99)+5/6*U(i,100)+1/3*U(i,101) )；%j=101 时取 w2=w3=0fweno(i,101)=l/3*U(i,99)-7/6*U(i,100)+11/6*U (i,101);for j =2:101fweno_pos(i,j)=(fweno(i,j)-fwno(i,j-
6、1)/dx; endfweno_pos(i,1)=0;endWENO_negfunction fweno_neg=WENO_neg(U) fwno_ng=zros(3,101);fweno=zeros(3,101);dx=0.01;2空间步长for i=l:3for j=3:99Ujl = l/3*U (iz j+2)-7/6*U(i,j+1)+11/6*U(i, j);Uj2=-l/6*U(i,j+l)+5/6*U(i,j)+l/3*U(i,j-1);Uj3=l/3*U(i,j)+5/6*U(i,j-1)-l/6*U(i,j-2);ISl = l/4*(U(i,j+2)-4*U(i,j+1)
7、+3*U (iA j) )A2 + 13/12* (U (i,j +2)-2*U(i,j+l)+U(i,j)人2;IS2 = l/4*(U(i,j+1)-U(i,j-1)A2 + 13/12* (U(i,j+1)-2*U (i,j )+U(i,j-l)A2;工S3=l/4*(3*U(i,j)-4*U(i,j-l)+U(i,j-2)A2 + 13/12*(U (i,j )-2*U(i,j-1)+U(i,j-2)人2;Cl=l/10;C2=6/10;C3=3/10;p=2;al=Cl/ (10A (-6)+IS1)Ap;a2=C2/(10A(-6)+IS2)Ap;a3=C3/ (10A (-6)+
8、IS3)Ap;wl=al/(al+a2+a3);w2=a2/(al+a2+a3);w3=a3/ (al+a2+a3);fweno(i,j)=wl*Ujl+w2*Uj 2+w3*Uj 3; end国处理边界条件%j =1 时取 w2=w3=0fweno(i,1)=1/3*U(i,3)一7/6夫U(i,2)+11/6*U(i, 1);% j =2 时取 w3=0ISl_2 = l/4*(U(i,4)-4*U(i,3)+3*U(i,2)2 + 13/12* (U(i,4)-2夫U?i, 3) +U (i, 2) )2工S2_2=l/4*(U(i,3)-U(i,1)A2+13/12*(U(i,3)-2
9、(i,2)+U( i,lT)A2;al_2=(1/10)/(10A(-6)+IS1_2)人2;322=(6/10)/(10A(-6)+工S212)八2;wl2=al_2/(al_2+a2_2); 一w2_2=a2_2/(al_2+a2_2);fweno(i,2)=wl_2*(1/3*U(i,4)-7/6夫U(i,3)+11/6*U(i,2)+w 2_2* (-1/6*U (匚3) +5/6*U (i, 2) +1/3夫U (i, 1);%j=100 时取wl=0IS2_10 0 = 1/4* (U(iz101)-U(iz99)A2 + 13/12*(U(i, 101) -2*U( i,100)
10、+U(i,99)A2;IS3_10 0 = 1/4*(3*U(iz100)-4*U(iA 99)+U(i,98)A2 + 13/12*(U (iA100)-2*U(i,99)+U(i,98)A2;a2_100=(6/10)/(10A(-6)+IS2_100)A2;a31100=(3/10)/(10人(-6)+工S3100)A2;w2100=a2_100/(a2_100+a3_100);w3100=a3100/(a2100+a3100);fweno(i,100)=w2_100*(-1/6*U(i,101)+5/6*U(i,100)+1/3*U (i,99)+w3_100*71/3*U(i,100
11、)+5/6*U(i,99)-l/6*U(i,98) 罚=101 时 取 wl=w2=0fweno(i,101)=l/3*U(i,101)+5/6*U(i,100)-1/6*U (iA 99);for j =1:100fweno_neg(i,j)=(fweno(i,j+1)-fweno(i,j)/dx; endfweno_neg(i,101)=0;end Steger-Warming具体步骤（以一维为例）己知 U = 0叫EV1)2)3)4)5)计兑 计修 计修 带入pmp4 = ，友=- C，4j = +C %, 4 (k = 123)利用不同的迎风格式，分别计算ar ardx dx（后差，前
12、差）6)、df ar sr计算十至十区diidu.F Cl= 0StSxsf+-dxdCdx(2 +2)122(7 -1)4 + & + 42(7- 1)Z+Z(一，)+ Z( +，)z(7-l)请+(4+母( + c+w(3-7)区+石1 H -2(7-D(1)f+ = f(x+),r = f(k-)7）时间推进箓黑。激波捕捉格式采取5阶WENO格式。方法如下：Jiang & Shu的5阶WENO格弋型+幽=o,仆)=皈0(1)正通量情况(a0)dx dx严二(刀缪。一 臂)/Ax鬻 =外力+电ZS，2=1/3/_2 _7/641 +11/6刀 以=一1/64T +5/64+1/347 学Z
13、3于于2-116于万3=- % = 攵= 1,23p = 2,e=10 电+(4,-2/ +3 科=% -4% +.%勖 2加 +加产 412412412j-2j-2j+1j+25阶WENO格式+ 要皿=o, f(ii) = au.a 0负通收情况(a0)dx dx=(惜-堂)z f=钠 +3注意是他.而不是j+1/2#3 = I%.2 -7/6加 +11/6力行/2 = T/6力.1+ 5/6 + 1/3加承2 =心力+ 5/6力.1/6加akCk%一 ,+。2+。3. =( + /S )P ”=1，2.3p=Z = l()fG =1/10,。2 =6/10,C3 .3/10电二；3“ -切
14、J3 +的)，电=2 -川+凯-2 +3珞。3力T& +川+浜 +人)2利用对称性，正通M差分格式中下标“j+k”改成“jk” 即得到负通地的差分格式、-,，WENO / ,WENO z*TTLEAY? / a - r *汪：fj =5+1/2 - /短）/ 不变空间网格数为100,即空间步长0.01。时间步长取O.OOlo时间推进格式采取3阶Runge-Kutta方法，如下列图：3阶（TVD型）。=Un + 加 0（U）。=3/4。 +1/4。 +1/440（。）Un+l= 1/3U +2/3U。）+ 2/340（。）程序包含：1 .主程序 WENO5_numerical_solution.子程序4个（1） FVS_pos（2） FVS_neg（3） WENO_pos（4） WENO_neg其中主程序主要进行时间推进和绘图。子程序采用流通矢量分裂(FVS)方法Steger-Warming分裂对通量进 行分裂，得到正通量和负通量。然后用5阶WENO格式对得到的通 量进行差分。得到的结果如图:1.10.110.90.80.70.50.40.30.2工0.6t=0.14s时刻密度RHO分布00.51x-0.51.51t=0.14s时亥ij速度U分布x1.210.8d 0.6
