计算程序_计算流体力学_对流方程_有限差分法_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时刻的计算结果') %%% 标题。