中心差分法 matlab代码_笔记|有限差分法|CFD Ferziger

干脆直接上完这一章的笔记，内容参考Ferziger的经典之作<Computational Methods for Fluid Dynamics^[1]>的第三章，主要简述了有限差分法的基本思路及各类常见差分的构造，方程求解部分没有太多涉及（求解部分被归纳到了第五章中）。

因为懒得码字所以前一部分全是图片，清晰度应该没问题；后一部分附算例（一维对流扩散方程）个人写的的MATLAB代码及结果图片。

感觉FDM本身就很自然，无论是差分近似微分，还是迎风格式（流场上游决定/影响下游）

*另外书里Spectral method几乎没怎么讲，所以也没有做太多笔记，自己写的一点主要来自于< Spectral methods in MATLAB^[2]>一书

笔记

算例 1-D convection diffusion equation （with Dirichlet boundary）

main

clear all
close all

%% parameters
L=1;
rho=1;
u=1;
Gamma=0.02;
Pe=rho*u*L/Gamma;
phi0=0;
phiL=1;

%% setting number of nodes
num=11;
phi=zeros(num,1);

%% initializing discretization
% uniform grid
x=zeros(1,num);
for i=1:num
    x(i)=(i-1)*L/(num-1);
end

% non-uniform grid
% for i=1:num
%     x(i)=1-1/2^(i-1);
% end

%% diffusion term
% CDS
Ad=zeros(num,3);
for i=2:(num-1)
    [Ad(i,1),Ad(i,2),Ad(i,3)]=diffusionCDS(Gamma,x(i-1),x(i),x(i+1));
end

%% convection term
% UDS (upwind difference scheme)
Ac=zeros(num,3);
for i=2:(num-1)
	[Ac(i,1),Ac(i,2),Ac(i,3)]=convectionUDS(rho,u,x(i-1),x(i),x(i+1));
end

% CDS
% for i=2:(num-1)
% 	[Ac(i,1),Ac(i,2),Ac(i,3)]=convectionCDS(rho,u,x(i-1),x(i),x(i+1));
% end


%% construct coefficient matrix
A=zeros(num);

A(1,1)=1;
A(num,num)=1;
for i=2:(num-1)
    A(i,i-1)=Ac(i,1)+Ad(i,1);
    A(i,i)=Ac(i,2)+Ad(i,2);
    A(i,i+1)=Ac(i,3)+Ad(i,3);
end

%% boundary condition
% Dirichlet
Q=zeros(num,1);
Q(1)=phi0;
Q(num)=phiL;

%% solution
% discretization solution
phi=inv(A)*Q;
plot(x,phi,'LineWidth',1);
hold on

% exact solution
num_exact=101;
x_exact=zeros(1,num_exact);
for i=1:num_exact
    x_exact(i)=(i-1)*L/(num_exact-1);
end
phi_exact=phi0+(exp(x_exact*Pe/L)-1)/(exp(Pe)-1)*(phiL-phi0);
plot(x_exact,phi_exact,'--','LineWidth',2);
legend(['UDS  node ',num2str(num)],'Exact solution')

functions

function [AW,AP,AE]=diffusionCDS(Gamma,x1,x2,x3)
% CDS of diffusion term

AW=-2*Gamma/((x3-x1)*(x2-x1));
AE=-2*Gamma/((x3-x1)*(x3-x2));
AP=-AW-AE;

end

function [AW,AP,AE]=convectionUDS(rho,u,x1,x2,x3)
% UDS for convection term

AE=min(rho*u,0)/(x3-x2);
AW=-max(rho*u,0)/(x2-x1);
AP=-AE-AW;

end

function [AW,AP,AE]=convectionCDS(rho,u,x1,x2,x3)
% CDS for convection term

AE=rho*u/(x3-x1);
AW=-rho*u/(x3-x1);
AP=-AE-AW;

end

Solutions

中心差分 11节点

中心差分 51节点

迎风格式 11节点

迎风格式 51节点

参考

1. ^Ferziger J H, Perić M, Street R L. Computational methods for fluid dynamics[M]. Berlin: springer, 2002.
2. ^Trefethen L N. Spectral methods in MATLAB[M]. Society for industrial and applied mathematics, 2000.