计算程序_计算流体力学_对流方程_有限差分法_Lax格式_迎风格式_FTCS格式

% 一维对流方程迎风格式、Lax格式、FTCS格式差分法计算

% 潭花林清华大学航天航空学院

% FTCS格式对于一维对流方程不稳定，最好不用

clc

clear all

% 1.参数定义

dx=1;

x1=-18;

x2=18;

x=x1:dx:x2;

L1=length(x);

% dt=0.5*dx; % 收敛

dt=2*dx; % 不收敛

t1=0;

t2=t1+80*dt;

t=t1:dt:t2;

L2=length(t);

alpha=1;

lambda=alpha*dt/dx;

geshi=1; % 迎风格式

% geshi=2; % Lax格式

% geshi=3; % FTCS格式

% 2.显式求解

zeta=zeros(L1,L2);

for kk=1:3

geshi=kk;

for ii=1:L1

if x(ii)>0

zeta(ii,1)=1;

else

if x(ii)==0

zeta(ii,1)=1/2;

else

if x(ii)<0

zeta(ii,1)=0;

end

end

end

end

if geshi==1

for ii=2:L1

for jj=1:(L2-1)

zeta(ii,jj+1)=zeta(ii,jj)-lambda*(zeta(ii,jj)-z eta(ii-1,jj));

end

zeta(1,jj+1)=zeta(2,jj+1);

end

zeta1=zeta;

else if geshi==2

for ii=2:(L1-1)

for jj=1:(L2-1)

zeta(ii,jj+1)=(zeta(ii+1,jj)+zeta(ii-1,jj))/2-. ..

lambda/2*(zeta(ii+1,jj)-zeta(ii-1,jj));

end

zeta(1,jj+1)=zeta(2,jj+1);

zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));

end

zeta2=zeta;

else if geshi==3

for ii=2:(L1-1)

for jj=1:(L2-1)

zeta(ii,jj+1)=zeta(ii,jj)-lambda/2*(zeta(ii+1,j j)-zeta(ii-1,jj));

end

zeta(1,jj+1)=zeta(2,jj+1);

zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));

end

zeta3=zeta;

end

end

end

end

% 3.绘图

% 3.1 t=0

figure(1)

n=1;

plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图

% axis equal %%% 是否要求x、y坐标间距相等

% grid on %%% 是否要求画网格

xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义

%text(1,2,'f(x)') %%% 图中文字标识

legend('迎风格式','Lax格式','FTCS格式

') %%% 不同曲线的线型区分

title('t=0时刻的计算结果') %%% 标题axis([-18,18,-0.2,1.2])

% 3.2 t=10

figure(2)

n=(10-t(1))/dt;

plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图

% axis equal %%% 是否要求x、y坐标间距相等

% grid on %%% 是否要求画网格

xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义

%text(1,2,'f(x)') %%% 图中文字标识

legend('迎风格式','Lax格式','FTCS格式

') %%% 不同曲线的线型区分

title('t=10s时刻的计算结果') %%% 标题

% 3.3 t=20

figure(3)

n=(20-t(1))/dt;

plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图

% axis equal %%% 是否要求x、y坐标间距相等

% grid on %%% 是否要求画网格

xlabel('x/m'),ylabel('t/s') %

%% x,y轴表示的变量含义

%text(1,2,'f(x)') %%% 图中文字标识

legend('迎风格式','Lax格式','FTCS格式

') %%% 不同曲线的线型区分

title('t=20s时刻的计算结果') %%% 标题

% 3.4 t=40

figure(4)

n=(40-t(1))/dt;

plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图

% axis equal %%% 是否要求x、y坐标间距相等

% grid on %%% 是否要求画网格

xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义

%text(1,2,'f(x)') %%% 图中文字标识

legend('迎风格式','Lax格式','FTCS格式

') %%% 不同曲线的线型区分

title('t=40s时刻的计算结果') %%% 标题