matlab-流体Sod激波管问题.zip

%常数参数 gamma = 1.4; CFL = 0.5; Time = 0.14; %时间 a = 0; %横坐标范围 b = 1; r=0.1; %网格比 dx = 0.005; %距离间隔 dt = dx*r; x = a:dx:b; t = linspace(0,Time,401); maxn = (b-a)/dx+1; nc = 2:maxn-1; %初始时刻的值，密度(rho)、速度(u)、压强(p)、能量(E) rho0 = 1.*(x<0.5)+0.125.*(x>=0.5); u0 = 0.*(x<0.5)+0.*(x>=0.5); p0 = 1.*(x<0.5)+0.1.*(x>=0.5); m0 = u0.*rho0; E0 = p0./((gamma-1).*rho0)+u0.^2/2; w0 = [rho0; m0; E0]; f0 = [u0.*w0(1,:) ; u0.*w0(2,:) ;u0.*w0(3,:)] + [u0-u0 ; p0 ; p0.*u0]; U0=[rho0 ;u0 ;p0]; U = U0; for i = 1:Time/dt Up = U; Un(1,:,1)=Up(1,1);Un(2,:,1)=Up(2,1);Un(3,:,1)=Up(3,1); for j = 2:maxn-1 A1 = U2A_sgn(Up, j ,-1);%非正特征值对应的矩阵A- A2 = U2A_sgn(Up, j ,1); %非负特征值对应的矩阵A+ A=A1 + A2 U( : ,j) = Up( : ,j) - r * (A1*(Up( : ,j+1)-Up( : ,j))... +A2*(Up(:,j)-Up(:,j-1)));%一阶迎风格式 Un(1,j,i)=U(1,j);Un(2,j,i)=U(2,j);Un(3,j,i)=U(3,j); end Un(:,maxn,i)=U(:,maxn); end %0.14S之后的参数值 rho = U(1,:); u = U(2,:); p = U(3,:); m = rho.*u; E = p/(gamma-1)+0.5*rho.*u.^2; %画图 figure(1); subplot(5,1,1)%密度rho plot(x,rho0,'b-') hold on plot(x,rho,'ro--') set(gca,'YLim',[-0.2 1.2]) set(gca,'XLim',[a b]) %set(gca, 'YTick', [0 1]) hold off title(['dx=',num2str(dx),' CFL=',num2str(CFL),' t=',num2str(Time)]); xlabel('长度x'); ylabel('密度rho'); subplot(5,1,2)%速度u plot(x,u0,'b-') hold on plot(x,u,'ro--') set(gca,'YLim',[-0.2 1.2]) set(gca,'XLim',[a b]) %set(gca, 'YTick', [0 1]) hold off xlabel('长度x'); ylabel('速度u'); subplot(5,1,3)%压力p plot(x,p0,'b-') hold on plot(x,p,'ro--') set(gca,'YLim',[-0.2 1.2]) set(gca,'XLim',[a b]) %set(gca, 'YTick', [0 1]) hold off xlabel('长度x'); ylabel('压力p'); subplot(5,1,4)%质量流m plot(x,m0,'b-') hold on plot(x,m,'ro--') set(gca,'YLim',[-0.5 2.5]) set(gca,'XLim',[a b]) %set(gca, 'YTick', [0 1]) hold off xlabel('长度x'); ylabel('质量流m'); subplot(5,1,5)%能量E plot(x,E0,'b--') hold on plot(x,E,'ro--') set(gca,'XLim',[a b]) hold off xlabel('长度x'); ylabel('能量E'); leg0 = legend('初始状态','末状态');%图例 横排 set(leg0, 'Position', [0.1425 0.88 0.3500 0.0256]);%修改图例位置 %画动图 %之后的参数值 rhon = Un(1,:,:); %密度rhon rhon = reshape(rhon,maxn,i); un = Un(2,:,:); un = reshape(un,maxn,i); pn = Un(3,:,:); pn = reshape(pn,maxn,i); mn = rhon.*un; mn = reshape(mn,maxn,i); En = pn/(gamma-1)+0.5*rhon.*un.^2; En = reshape(En,maxn,i); figure(2) i1=i; M1=moviein(i1); for ii1=1:i1 plot(x,rhon(:,ii1),x,rhon(:,ii1),'ro--') set(gca,'YLim',[-0.2 1.2]) set(gca,'XLim',[a b]) xlabel('长度x'); ylabel('密度rho'); M1(i1)=getframe; end %movie(M1); figure(3)%速度u i2=i1; M2=moviein(i2); for ii2=1:i2 plot(x,un(:,ii2),x,un(:,ii1),'ro--') set(gca,'YLim',[-0.2 1.2]) set(gca,'XLim',[a b]) xlabel('长度x'); ylabel('速度u'); M2(:,i2)=getframe; end %movie(M2,1); figure(4)%压力p i3=i2; M3=moviein(i3); for ii3=1:i3 plot(x,pn(:,ii3),x,pn(:,ii1),'ro--') set(gca,'YLim',[-0.2 1.2]) set(gca,'XLim',[a b]) xlabel('长度x'); ylabel('压力p'); M3(:,i3)=getframe; end %movie(M3,1); figure(5)%质量流m i4=i3; M4=moviein(i4); for ii4=1:i4 plot(x,mn(:,ii4),x,mn(:,ii1),'ro--') set(gca,'YLim',[0 1.2]) set(gca,'XLim',[a b]) xlabel('长度x'); ylabel('质量流m'); M4(:,i4)=getframe; end %movie(M4,1); figure(6)%能量E i5=i4; M5=moviein(i5); for ii5=1:i5 plot(x,En(:,ii5),x,En(:,ii1),'ro--') set(gca,'YLim',[0 1.2]) set(gca,'XLim',[a b]) xlabel('长度x'); ylabel('能量E'); M5(:,i5)=getframe; end %movie(M5,1); figure(2); %函数U2A_sgn function [A_sgn] = U2A_sgn(U, j, sgn) g = 1.4 ;%gamma rho = U(1,j); u = U(2,j); p = U(3,j); a = (g*p/rho)^0.5; A = [ u rho 0 0 u 1/rho 0 g*p u ]; L = [ 0 -rho*a 1 a^2 0 -1 0 rho*a 1]; D = [ u-a 0 0 0 u 0 0 0 u+a];%A矩阵对角化之后 A_abs = inv(L)*(sign(D).*D)*L; A_sgn = 0.5*(A + sgn*A_abs); end