有限查分计算2D稳态热传导方程

下面的有限差分程序计算2D稳态热传导方程。

其模型和文章有限元计算2D稳态热传导方程的模型一样。从后面也可以看到，其结果也和该文章的结果一样。

即有限差分和有限元法，解决同一问题，得到了相同的解。

程序见下面：

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

% Solution of 2D steady heat conduction equation(Poisson equation): d2FI/dx2+d2FI/dy2=f
% with finite differences on a regular grid
%   here f=qv/K, where qv is volumetric heat production rate, K is thermal conductivity

% Clean all variables
clear;
% Clear all figures
clf;

% Numerical model parameters
% Model size, m
xsize   =   1; % Horizontal
ysize   =   1; % Vertical

% Numbers of nodes
xnum    =   13; % Horizontal
ynum    =   13; % Vertical
% Grid step
xstp    =   xsize/(xnum-1); % Horizontal
ystp    =   ysize/(ynum-1); % Vertical

 f=1;
%f=qv/K

% Matrix of coefficients initialization
L       =   sparse(xnum*ynum,xnum*ynum);
% Vector of right part initialization
R       =   zeros(xnum*ynum,1);

% Composing matrix of coefficients L()
% and vector (column) of right parts R()
% Boundary conditions: FI=1 for y==0, FI=0 for other boundaries
% Process all Grid points
for i=1:1:ynum
  for j=1:1:xnum
    % Global index for current node
    k=(j-1)*ynum+i;
    %Boundary nodes
    if(j==1 || j==xnum || i==xnum)
        L(k,k)=1;
        R(k,1)=0;
     elseif i == 1
        L(k,k)=1;
        R(k,1)=1;
    else
    %Internal nodes
        %d2FI/dx2+d2FI/dy2=f
        L(k,k-ynum) = 1/xstp^2;
        L(k,k-1   ) = 1/ystp^2;
        L(k,k     ) =-2/xstp^2-2/ystp^2;
        L(k,k+1   ) = 1/ystp^2;
        L(k,k+ynum) = 1/xstp^2;
        R(k,     1) = f;
    end
  end
end
    
%Obtaining vector of solutions S()
S       =   L\R;

% Reload solutions S() to 2D array FI()
FI      =   zeros(ynum,xnum);
% Process all Grid points
for i=1:1:ynum
  for j=1:1:xnum
    % Global index for current node
    k       =   (j-1)*ynum+i;
    FI(i,j) =   S(k);
  end
end

% Making vectors for nodal points positions
x=0:xstp:xsize; % Horizontal
y=0:ystp:ysize; % Vertical

% Making new figure
figure(1);
% Plotting solution
surf(x,y,FI);
shading interp;
xlabel('x, km')
ylabel('y, km')
 view ( -67.5, 30 );

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

其结果如下：

可见，和使用有限元法得到的结果一致。