1、精選優質文檔-傾情為你奉上1 提出問題 問題描述 Sod激波管問題是典型的一類Riemann問題。如圖所示，一管道左側為高溫高壓氣體，右側為低溫低壓氣體，中間用薄膜隔開。t=0 時刻，突然撤去薄膜，試分析其他的運動。 Sod模型問題：在一維激波管的左側初始分布為：，右側分布為：，兩種狀態之間有一隔膜位于處。隔膜突然去掉，試給出在時刻Euler方程的準確解，并給出在區間這一時刻的分布圖。2 一維Euler方程組分析可知，一維激波管流體流動符合一維Euler方程，具體方程如下：矢量方程：分量方程：連續性方程、動量方程和能量方程分別是：其中 對于完全氣體，在量綱為一的形式下，狀態方程為：在量綱為一的

2、定義下，定容熱容為：聯立(1.2)，(1.3)，(1.4)消去溫度和定容比熱，得到氣體壓力公式為：上式中為氣體常數，對于理想氣體。3 Euler方程組的離散3.1 Jacibian矩陣特征值的分裂Jacibian矩陣A的三個特征值分別是，依據如下算法將其分裂成正負特征值：3.2 流通矢量的分裂這里對流通矢量的分裂選用Steger-Warming分裂法，分裂后的流通矢量為其中：為量綱為一的聲速：聯立(1.3)，(1.9)式，消去來流馬赫數得：3.3 一階迎風顯示格式離散Euler方程組得到 算法如下： 已知初始時刻t=0的速度、壓力及密度分布，則可得到特征值分裂值，從而求出流通矢量； 應用一階迎

3、風顯示格式可以計算出時刻的組合變量，從而得到時刻的速度、壓力及密度分布； 利用時刻的速度、壓力及密度分布可得特征值分裂值，從而求出流通矢量； 按照步驟2的方法即可得到時刻的速度、壓力及密度分布； 循環以上過程即可得到時刻的速度、壓力及密度分布。4 計算結果分析實際編程中，空間步長取0.001，空間網格數為1001，時間步長取0.00001，計算到終點時刻0.14s耗費機時137s，計算時間還是可以接受的。分析圖4-14-3，可以觀察到在隔膜附近流動參數變化劇烈，與初始條件相比，可以看出激波的影響范圍有限，始終在區間內變化。圖4-1是0.14時刻的密度分布圖，觀察可知，在密度波的傳播過程中，間斷

4、面上會出現了兩次“沉降”，說明密度在沉降位置發生了劇烈變化。圖4-2是0.14時刻的壓力分布圖，在壓力波的傳播過程中，在間斷面上出現了一個“壓力沉降”現象，說明壓力在沉降位置突降。圖4-3是0.14時刻的速度分布圖，在間斷面處產生一個向兩邊運動的速度，并且只有在隔膜附近才有氣體流動，其他地方靜止。圖4-1激波管內密度分布圖(0.14s)圖4-2激波管內壓力分布圖(0.14s)圖4-3激波管內速度分布圖(0.14s)源程序代碼： 分量和矩陣結合編寫的源程序：function sobtubing_SW()tic;close allee=1e-8;%劃分時空網格%delta_t=0.00001;Nt

5、=round(0.14/delta_t);delta_x=0.001;N_left=round(0.5/delta_x+1);N_right=round(0.5/delta_x+1)N=N_left+N_right-1;%1%初始條件%P=ones(1,N_left-1) 0.1*ones(1,N_right);Den=ones(1,N_left-1) 0.125*ones(1,N_right);u=zeros(1,N);gama=1.4;Den_u=Den.*u;E=P./(gama-1) +(0.5*u.2).*Den; %計算特征值分裂%for j=1:Ntepso=1e-8*ones(

6、1,N);gama=1.4;C=sqrt(gama*P./Den);lamta=ones(3,N);lamta_p=ones(3,N);lamta_n=ones(3,N);lamta(1,:)=u; lamta(2,:)=u-C; lamta(3,:)=u+C;for i=1:3 lamta_p(i,:)=0.5*(lamta(i,:)+sqrt(lamta(i,:).2+epso.2); lamta_n(i,:)=0.5*(lamta(i,:)-sqrt(lamta(i,:).2+epso.2);end %計算正通量%gama=1.4;C=sqrt(gama*P./Den);Ftran_p=

7、ones(3,N);f1_Pos=0.5/gama*Den.*(2*(gama-1)*lamta_p(1,:)+lamta_p(2,:)+lamta_p(3,:);f2_Pos=0.5/gama*Den.*(2*(gama-1)*lamta_p(1,:).*u+lamta_p(2,:).*(u-C)+lamta_p(3,:).*(u+C);f3_Pos=0.5/gama*Den.*(gama-1)*lamta_p(1,:).*u.2+0.5*lamta_p(2,:).*(u-C).2 .+0.5*lamta_p(3,:).*(u+C).2+(0.5*(3-gama)/(gama-1)*(lam

8、ta_p(2,:)+lamta_p(3,:).*C.2); %計算負通量%gama=1.4;C=sqrt(gama*P./Den);f1_Neg=0.5/gama*Den.*(2*(gama-1)*lamta_n(1,:)+lamta_n(2,:)+lamta_n(3,:);f2_Neg=0.5/gama*Den.*(2*(gama-1)*lamta_n(1,:).*u+lamta_n(2,:).*(u-C)+lamta_n(3,:).*(u+C);f3_Neg=0.5/gama*Den.*(gama-1)*lamta_n(1,:).*u.2+0.5*lamta_n(2,:).*(u-C).2

9、 .+0.5*lamta_n(3,:).*(u+C).2+(0.5*(3-gama)/(gama-1)*(lamta_n(2,:)+lamta_n(3,:).*C.2); %計算流動參數% for i=2:N-1 %密度計算 temp1(1,i) = (f1_Pos(1,i) - f1_Pos(1,i-1) / delta_x; temp2(1,i) = (f1_Neg(1,i+1) - f1_Neg(1,i) / delta_x; Den(1,i) = Den(1,i) - delta_t*(temp1(1,i) + temp2(1,i); % 密度、速度乘積計算 temp1(1,i) =

10、(f2_Pos(1,i) - f2_Pos(1,i-1) / delta_x; temp2(1,i) = (f2_Neg(1,i+1) - f2_Neg(1,i) /delta_x; Den_u(1,i) = Den_u(1,i) - delta_t*(temp2(1,i) + temp1(1,i); % 速度計算 u(1,i) = Den_u(1,i) / Den(1,i); % 能量計算 temp1(1,i) = (f3_Pos(1,i) - f3_Pos(1,i-1) /delta_x; temp2(1,i) = (f3_Neg(1,i+1) - f3_Neg(1,i) /delta_x

11、; E(1,i) = E(1,i) - delta_t*(temp2(1,i) + temp1(1,i); % 壓強計算 P(1,i) = (gama - 1)*(E(1,i) - 0.5*Den(1,i)*u(1,i)2); endend%結果顯示x=0:0.001:1;axis(0 1 0 1);plot(x,u,'LineWidth',2);xlabel('fontsize14x');ylabel('fontsize14速度V');figure(2);plot(x,P,'LineWidth',2);xlabel('fontsi