計(jì)算流體力學(xué)大作業(yè)

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

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

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

4、面上會(huì)出現(xiàn)了兩次“沉降”，說(shuō)明密度在沉降位置發(fā)生了劇烈變化。圖4-2是0.14時(shí)刻的壓力分布圖，在壓力波的傳播過(guò)程中，在間斷面上出現(xiàn)了一個(gè)“壓力沉降”現(xiàn)象，說(shuō)明壓力在沉降位置突降。圖4-3是0.14時(shí)刻的速度分布圖，在間斷面處產(chǎn)生一個(gè)向兩邊運(yùn)動(dòng)的速度，并且只有在隔膜附近才有氣體流動(dòng)，其他地方靜止。圖4-1激波管內(nèi)密度分布圖(0.14s)圖4-2激波管內(nèi)壓力分布圖(0.14s)圖4-3激波管內(nèi)速度分布圖(0.14s)源程序代碼： 分量和矩陣結(jié)合編寫(xiě)的源程序：function sobtubing_SW()tic;close allee=1e-8;%劃分時(shí)空網(wǎng)格%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; %計(jì)算特征值分裂%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 %計(jì)算正通量%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); %計(jì)算負(fù)通量%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); %計(jì)算流動(dòng)參數(shù)% for i=2:N-1 %密度計(jì)算 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); % 密度、速度乘積計(jì)算 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); % 速度計(jì)算 u(1,i) = Den_u(1,i) / Den(1,i); % 能量計(jì)算 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); % 壓強(qiáng)計(jì)算 P(1,i) = (gama - 1)*(E(1,i) - 0.5*Den(1,i)*u(1,i)2); endend%結(jié)果顯示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