CFD入门: 有限差分和基本概念

upup，我试着写了一下BTCS的格式，但是不知道怎么加周期性边界条件[委屈][委屈]，可以帮我看看嘛[呲牙]

隐格式：

时间上采用后向差分，空间上采用中心差分：令CFL=Δt/Δx

格式：u(j,n+1)+CFL/2*[ u(j+1,n+1)- u(j-1,n+1)]= u(j,n)

方程组：A*U(n+1)=U(n)

A=diag(ones(1,Nx+1))+diag(-1/2*CFL*ones(1,Nx),-1)+diag(1/2*CFL*ones(1,Nx),1)

U(n+1)=[u(1,n+1), u(2,n+1)…u(N,n+1)]

U(n)=[u(1,n)+CFL/2*u(0,n+1),u(2,n)…u(N-1,n), u(N,n)- CFL/2*u(N+1,n+1)]

代码如下

clear

xa = 0; xb = 2*pi; % 求解域

Nx = 200; % 单元数

Nx1 = Nx + 1; % 格点数

hx = (xb - xa)/Nx; % Δx

Nt = 500; % 单元数

tend = 2*pi; % 终止时间

Nt1 = Nt + 1; % 格点数

ht = (tend - 0)/Nt; % Δx

% 真解

ureal = @(x,t) sin(x-t); 

bcL = 1; bcR = 1; 

CFL = ht/hx; % CFL数，Δt/Δx

X = xa:hx:xb; % 网格

% 系数矩阵

U=zeros(Nx1,Nt1);

U(:,1)=ureal(X,0);

A= zeros(Nx1,Nx1);

A = A + diag(ones(1,Nx + 1));

A = A + diag(-1/2*CFL*ones(1,Nx),-1);

A = A + diag(1/2*CFL*ones(1,Nx),1);

F = zeros(Nx1,1); % 右端项

for i=2:Nt1 

  F(1)=U(2,i-1)+1/2*CFL*U(1,i);

  F(end)=U(end-1,i-1)-1/2*CFL*U(end,i);

  for j=2:Nx1-1

    F(j)=U(j,i-1);

  end

  U(:,i)=A\F;

end

figure

plot(X,U(:,end),'ro',X,ureal(X,tend),'b-')

legend('num','exact');