激波管兩端封閉 內部充滿氣體 直管中有一薄膜將激波管隔開 在薄膜兩側充有均勻理想氣體

薄膜兩側氣體的壓力密度不同當t<0時 氣體處于靜止狀態

當t>=0時薄膜瞬時突然破裂氣體從高壓端沖向低壓端

控制方程有一維歐拉方程和理想氣體方程：

邊界條件：

初始條件：

MacCormack兩步差分迭代格式：

人工黏性項：

----------------------------------------------------------------------------

用到七個函數

主函數 shockwave_1D

子函數

timestep、MacCormack、U2F、boundarycondition、output

用來后處理的函數 post

----------------------------------------------------------------------------

function shockwave_1D()

global GAMA J CN

%全局變量

GAMA=1.4; %氣體常數

L=2.0; %計算區域

options={'網格數J';'總時間TT';'時間步長因子CN'};

inp={'200';'0.4';'0.8'};

dlg=inputdlg(options,'input',1,inp);

J=eval_r(dlg{1});?%網格數

TT=eval_r(dlg{2});?%總時間

CN=eval_r(dlg{3});?%時間步長因子

%輸入網格數、總時間、時間步長因子

rou1=1.0;

u1=0.0;

p1=1.0;

rou2=0.125;

u2=0.0;

p2=0.1;

%兩側的密度速度壓力

dx=L/J;

U=zeros(J+1,3);

for i=1:round(J/2)

U(i,1)=rou1;

U(i,2)=rou1*u1;

U(i,3)=p1/(GAMA-1)+rou1*u1*u1/2;

end

for i=round(J/2)+1:J+1

U(i,1)=rou2;

U(i,2)=rou2*u2;

U(i,3)=p2/(GAMA-1)+rou2*u2*u2/2;

end

%初始化矩陣U

t=0;?%時間變量

fp=zeros(3,1);

fp(1)=fopen('density.txt','a');

fp(2)=fopen('velocity.txt','a');

fp(3)=fopen('pressure.txt','a');

%建立三個文本文件 用來記錄每步的結果

for i=1:3

fprintf(fp(i),'?time?');

for

j=0:J

fprintf(fp(i),'?�',j*dx);

end

end

%初始化三個結果文件

%開始迭代計算

while t

dt=timestep(U,dx);

%求時間步長

t=t+dt;

U=MacCormack(U,dx,dt);

%調用maccormack迭代格式 并將新值賦給U

U=boundarycondition(U);

%邊界條件

output(U,t,fp);

%將計算結果輸出到文本文件

end

for i=1:3

fclose(fp(i));

end

%關閉結果文本文件

----------------------------------------------------------------------------

function

dt=timestep(U,dx)

%函數 用來計算時間步長

global GAMA J CN

maxvel=1e-100;

for i=2:J;

u=U(i,2)/U(i,1);

p=(GAMA-1)*(U(i,3)-0.5*U(i,1)*u*u);

vel=sqrt(GAMA*p/U(i,1))+abs(u);

if

vel>maxvel

maxvel=vel;

end

end

%獲得局部最大速度

dt=CN*dx/maxvel;

----------------------------------------------------------------------------

function U2=

MacCormack(U,dx,dt)

%兩步差分的MacCormack迭代格式 Ut和Ft 用來臨時變量存儲

global J

r=dt/dx;

eta=0.25;

%帶有密度控制開關的人工黏性方法 其中參數 eta theta

for i=2:J

theta=abs(abs(U(i+1,1)-U(i,1))-abs(U(i,1)-U(i-1,1)))...

/(abs(U(i+1,1)-U(i,1))+abs(U(i,1)-U(i-1,1))+1e-100);

%開關函數的參數theta

for

k=1:3

Ft(i,k)=U(i,k)+0.5*eta*theta*(U(i+1,k)-2*U(i,k)+U(i-1,k));%人工黏性項

end

end

for k=1:3

for

i=2:J

U(i,k)=Ft(i,k);

end

end

for i=1:J+1

Ft(i,:)=U2F(U(i,:));

end

for i=1:J

for

k=1:3

Ut(i,k)=U(i,k)-r*(Ft(i+1,k)-Ft(i,k));

end

end

%U(n+1/2)(i+1/2)

for i=1:J

Ft(i,:)=U2F(Ut(i,:));?end

%F(n+1/2)(i+1/2)

for i=2:J

for

k=1:3

U(i,k)=0.5*(U(i,k)+Ut(i,k))-0.5*r*(Ft(i,k)-Ft(i-1,k));

end

end

%U(n+1)(i)

U2=U;

----------------------------------------------------------------------------

function F=U2F(U)

%由矩陣U得到矩陣F的函數

global GAMA

u=U(2)/U(1);

p=(GAMA-1)*(U(3)-0.5*U(2)*U(2)/U(1));

F(1)=U(2);

F(2)=U(1)*u*u+p;

F(3)=(U(3)+p)*u;

----------------------------------------------------------------------------

function

U2=boundarycondition(U)

global J

U2=U;

for i=1:3

U2(1,i)=U(2,i);

U2(J+1,i)=U(J,i);

end

%自由流入流出 U2=U1 Un=Un-1

----------------------------------------------------------------------------

function output(U,t,fp)

%函數 當前步的計算結果 輸出到文本文件

global J GAMA

for i=1:3

fprintf(fp(i),'\nf?',t);

end

%輸出當前時間t

for i=1:J+1

rou=U(i,1);

u=U(i,2)/rou;

p=(GAMA-1)*(U(i,3)-0.5*U(i,1)*u*u);

fprintf(fp(1),'?.5e',rou);

fprintf(fp(2),'?.5e',u);

fprintf(fp(3),'?.5e',p);

end

%輸出對應內容 密度 速度 壓力

----------------------------------------------------------------------------

function post()

%后處理文件 調用輸出的文本文件

m=menu('選擇要輸出的內容','density','velocity','pressure')

if m==1

option='density';

end

if m==2

option='velocity';

end

if m==3

option='pressure';

end

%選擇要輸出的內容

filename=strcat(option,'.txt');

A=importdata(filename);

B=A.rowheaders;

A=A.data;

[a,b]=size(A);

for t=2:a

plot(A(1,:),A(t,:));

%畫圖

xlab='?t=';

ind=int2str(t-1);

while

length(ind)<3

ind=strcat('0',ind);

end

xlabel(strcat('x',xlab,B(t)));

ylabel(option);

frame =

getframe(1);

im =

frame2im(frame);

filename=strcat(option,ind);

filename=strcat(filename,'.bmp');

imwrite(im,filename,'bmp');

%將圖片逐幀輸出

end

----------------------------------------------------------------------------

計算結果如下圖

輸出t=0.4s時刻 密度 速度 壓強

總結

以上是生活随笔為你收集整理的激波管的matlab程序,matlab中用MacCormack法求解激波喷管问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。