本文程序，使用有限元计算2D稳态热传导方程。
本程序只使用2D三角形单元求解，没有四边形单元求解。不过，也非常容易就可做到。
主程序为:
Heat_Conduction_steady.m
需要用到的子程序有：
coord3.m 用来生成节点坐标，节点坐标与单元的链接。边界条件等等。
生成文件
coordinates.dat 包含两列，第一列节点x坐标，第二列，节点y坐标
elements3.dat 包含三列：分别是各个单元的i,j,k节点的总体编号。
dirichlet.dat 包含两列：分别是第一类边界条件的两端点节点的总体编号。
另外，
coord4.m与coord3相同，但是，生成的是四边形单元。
即：
elements4.dat
f.m 生成由内生热产生的节点载荷qv/K,其中，qv为体积生热率，K为热导率
g.m 生成由表面热流产生的节点载荷。在有第二类边界条件的时候才使用。
stima3.m 计算三角形单元的单元刚度矩阵
stima4.m 计算四边形单元的单元刚度矩阵
u_d.m 狄利克雷边界条件。
先将各个函数列于下方：
----------------------------------------------------------------------------------------
% coord3.m
% the coordinate index is from 1~Nx(from left to right) for the bottom line
% and then from Nx+1~2*Nx (from left to right) for the next line above
% and then next
xmin=0;xmax=1;
ymin=0;ymax=1;
Nx=13;Ny=13;
x=linspace(xmin,xmax,Nx);
y=linspace(ymin,ymax,Ny);
k=0;
for i1=1:Ny
for i2=1:Nx
k=k+1;
Coord(k,:)=[x(i2),y(i1)];
end
end
save coordinates.dat Coord -ascii
% the element index
k=0;
vertices=zeros((Nx-1)*(Ny-1)*2,3);
for i1=1:Ny-1
for i2=1:Nx-1
k=k+1;
ijm1=i2+(i1-1)*Nx;
ijm2=i2+1+(i1-1)*Nx;
ijm3=i2+1+i1*Nx;
vertices(k,:)=[ijm1,ijm2,ijm3];
k=k+1;
ijm1=i2+1+i1*Nx;
ijm2=i2+i1*Nx;
ijm3=i2+(i1-1)*Nx;
vertices(k,:)=[ijm1,ijm2,ijm3];
end
end
save elements3.dat vertices -ascii
% The direchlet boundary condition (index of the two end nodes for each boundary line)
boundary=zeros(2*(Nx+Ny-2),2);
temp1=1:Nx-1;
temp2=2:Nx;
temp3=1:Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Nx*(1:Ny-1);
temp2=temp1+Nx;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Ny*Nx:-1:Ny*Nx-Nx+2;
temp2=temp1-1;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Nx*(Ny-1)+1:-Nx:Nx+1;
temp2=temp1-Nx;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
save dirichlet.dat boundary -ascii
----------------------------------------------------------------------------------------
% coord4.m
% the coordinate index is from 1~Nx(from left to right) for the bottom line
% and then from Nx+1~2*Nx (from left to right) for the next line above
% and then next
xmin=0;xmax=1;
ymin=0;ymax=1;
Nx=13;Ny=13;
x=linspace(xmin,xmax,Nx);
y=linspace(ymin,ymax,Ny);
k=0;
for i1=1:Ny
for i2=1:Nx
k=k+1;
Coord(k,:)=[x(i2),y(i1)];
end
end
save coordinates.dat Coord -ascii
% the element index
k=0;
vertices=zeros((Nx-1)*(Ny-1),4);
for i1=1:Ny-1
for i2=1:Nx-1
k=k+1;
il1=(i1-1)*Nx+i2;
il2=il1+1;
il3=il2+Nx;
il4=il1+Nx;
vertices(k,:)=[il1,il2,il3,il4];
end
end
save elements4.dat vertices -ascii
% The direchlet boundary condition (index of the two end nodes for each boundary line)
boundary=zeros(2*(Nx+Ny-2),2);
temp1=1:Nx-1;
temp2=2:Nx;
temp3=1:Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Nx*(1:Ny-1);
temp2=temp1+Nx;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Ny*Nx:-1:Ny*Nx-Nx+2;
temp2=temp1-1;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
temp1=Nx*(Ny-1)+1:-Nx:Nx+1;
temp2=temp1-Nx;
temp3=temp3+Nx-1;
boundary(temp3',:)=[temp1',temp2'];
save dirichlet.dat boundary -ascii
----------------------------------------------------------------------------------------
% 主程序Heat_Conduction_steady.m
%*****************************************************************************80
%
%% Aapplies the finite element method to steady heat conduction equation,
% that is Poission's equation
%
%
% The user supplies datafiles that specify the geometry of the region
% and its arrangement into triangular or quadrilateral elements, and
% the location and type of the boundary conditions, which can be any
% mixture of Neumann and Dirichlet.
%
% The unknown state variable U(x,y) is assumed to satisfy
% Poisson's equation (steady heat conduction equation):
% -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega
% with Dirichlet boundary conditions
% U(x,y) = U_D(x,y) on Gamma_D
% and Neumann boundary conditions on the outward normal derivative:
% Un(x,y) = G(x,y) on Gamma_N
% If Gamma designates the boundary of the region Omega,
% then we presume that
% Gamma = Gamma_D + Gamma_N
% F(x,y) is qv/K,here qv means volumetric heat generation rate
% but the user is free to determine which boundary conditions to
% apply.
%
% The code uses piecewise linear basis functions for triangular elements,
% and piecewise isoparametric bilinear basis functions for quadrilateral
% elements.
%
%
%
clear
%
% Read the nodal coordinate data file.
%
load coordinates.dat;
%
% Read the triangular element data file.
%
load elements3.dat;
%
% Read the quadrilateral element data file.
%
% load elements4.dat;
%
% Read the Neumann boundary condition data file.
% I THINK the purpose of the EVAL command is to create an empty NEUMANN array
% if no Neumann file is found.
%
eval ( 'load neumann.dat;', 'neumann=[];' );
%
% Read the Dirichlet boundary condition data file.
%
load dirichlet.dat;
A = sparse ( size(coordinates,1), size(coordinates,1) );
b = sparse ( size(coordinates,1), 1 );
%
% Assembly.
%
for j = 1 : size(elements3,1)
A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...
+ stima3(coordinates(elements3(j,:),:));
end
%{
for j = 1 : size(elements4,1)
A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) ...
+ stima4(coordinates(elements4(j,:),:));
end
%}
% Volume Forces.
%
% from the center of each element to Nodes
% Notice that the result of f here means qv/K instead of qv
for j = 1 : size(elements3,1)
b(elements3(j,:)) = b(elements3(j,:)) ...
+ det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...
f(sum(coordinates(elements3(j,:),:))/3)/6;
end
%{
for j = 1 : size(elements4,1)
b(elements4(j,:)) = b(elements4(j,:)) ...
+ det([1,1,1; coordinates(elements4(j,1:3),:)'] ) * ...
f(sum(coordinates(elements4(j,:),:))/4)/4;
end
%}
% Neumann conditions.
%
if ( ~isempty(neumann) )
for j = 1 : size(neumann,1)
b(neumann(j,:)) = b(neumann(j,:)) + ...
norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...
g(sum(coordinates(neumann(j,:),:))/2)/2;
end
end
%
% Determine which nodes are associated with Dirichlet conditions.
% Assign the corresponding entries of U, and adjust the right hand side.
%
u = sparse ( size(coordinates,1), 1 );
BoundNodes = unique ( dirichlet );
u(BoundNodes) = u_d ( coordinates(BoundNodes,:) );
b = b - A * u;
%
% Compute the solution by solving A * U = B for the remaining unknown values of U.
%
FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );
u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
%
% Graphic representation.
%
% show ( elements3, elements4, coordinates, full ( u ) );
figure
trisurf ( elements3, coordinates(:,1), coordinates(:,2), (full ( u )));
view ( -67.5, 30 );
shading interp
----------------------------------------------------------------------------------------
%stima3.m
function M = stima3 ( vertices )
%*****************************************************************************80
%
%% STIMA3 determines the local stiffness matrix for a triangular element.
%
% Discussion:
%
% Although this routine is intended for 2D usage, the same formulas
% work for tetrahedral elements in 3D. The spatial dimension intended
% is determined implicitly, from the spatial dimension of the vertices.
%
% Parameters:
%
% Input, real VERTICES(1:(D+1),1:D), contains the D-dimensional
% coordinates of the vertices.
%
% Output, real M(1:(D+1),1:(D+1)), the local stiffness matrix
% for the element.
%
d = size ( vertices, 2 );
D_eta = [ ones(1,d+1); vertices' ] \ [ zeros(1,d); eye(d) ];
M = det ( [ ones(1,d+1); vertices' ] ) * D_eta * D_eta' / prod ( 1:d );
----------------------------------------------------------------------------------------
%stima4.m
function M = stima4 ( vertices )
%*****************************************************************************80
%
%% STIMA4 determines the local stiffness matrix for a quadrilateral element.
%
%
% Parameters:
%
% Input, real VERTICES(1:4,1:2), contains the coordinates of the vertices.
%
% Output, real M(1:4,1:4), the local stiffness matrix for the element.
%
D_Phi = [ vertices(2,:) - vertices(1,:); vertices(4,:) - vertices(1,:) ]';
B = inv ( D_Phi' * D_Phi );
C1 = [ 2, -2; -2, 2 ] * B(1,1) ...
+ [ 3, 0; 0, -3 ] * B(1,2) ...
+ [ 2, 1; 1, 2 ] * B(2,2);
C2 = [ -1, 1; 1, -1 ] * B(1,1) ...
+ [ -3, 0; 0, 3 ] * B(1,2) ...
+ [ -1, -2; -2, -1 ] * B(2,2);
M = det ( D_Phi ) * [ C1 C2; C2 C1 ] / 6;
----------------------------------------------------------------------------------------
%f.g
function value = f ( u )
%*****************************************************************************80
%
%% F evaluates the right hand side of Laplace's equation. NOtice that F is qv/K instead of qv.
%
%
% This routine must be changed by the user to reflect a particular problem.
%
%
% Parameters:
%
% Input, real U(N,M), contains the M-dimensional coordinates of N points.
%
% Output, VALUE(N), contains the value of the right hand side of Laplace's
% equation at each of the points.
%
n = size ( u, 1 );
value(1:n) = ones(n,1);
%2.0 * pi * pi * sin ( pi * u(1:n,1) ) .* sin ( pi * u(1:n,2) );
----------------------------------------------------------------------------------------
%g.m
function value = g ( u )
%*****************************************************************************
%
%% G evaluates the outward normal values assigned at Neumann boundary conditions.
%
%
% This routine must be changed by the user to reflect a particular problem.
%
% Parameters:
%
% Input, real U(N,M), contains the M-dimensional coordinates of N points.
%
% Output, VALUE(N), contains the value of outward normal at each point
% where a Neumann boundary condition is applied.
%
value = zeros ( size ( u, 1 ), 1 );
----------------------------------------------------------------------------------------
%u_d.m
function value = u_d ( u )
%*****************************************************************************80
%
%% U_D evaluates the Dirichlet boundary conditions.
%
%
% The user must supply the appropriate routine for a given problem
%
%
% Parameters:
%
% Input, real U(N,M), contains the M-dimensional coordinates of N points.
%
% Output, VALUE(N), contains the value of the Dirichlet boundary
% condition at each point.
%
value = zeros ( size ( u, 1 ), 1 );
value(1:13)=1;
----------------------------------------------------------------------------------------
上面的程序计算了内生热为均匀qv/K=1，在底边界温度为1，在其它边界温度为0的计算的结果：
结果如下图
对于四边形单元，将主函数Heat_Conduction_steady.m修改如下
--------------------------------------------------------------------------------------
%Heat_Conduction_steady.m
%*****************************************************************************80
%
%% Aapplies the finite element method to steady heat conduction equation,
% that is Poission's equation
%
%
% The user supplies datafiles that specify the geometry of the region
% and its arrangement into triangular or quadrilateral elements, and
% the location and type of the boundary conditions, which can be any
% mixture of Neumann and Dirichlet.
%
% The unknown state variable U(x,y) is assumed to satisfy
% Poisson's equation (steady heat conduction equation):
% -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega
% with Dirichlet boundary conditions
% U(x,y) = U_D(x,y) on Gamma_D
% and Neumann boundary conditions on the outward normal derivative:
% Un(x,y) = G(x,y) on Gamma_N
% If Gamma designates the boundary of the region Omega,
% then we presume that
% Gamma = Gamma_D + Gamma_N
% F(x,y) is qv/K,here qv means volumetric heat generation rate
% but the user is free to determine which boundary conditions to
% apply.
%
% The code uses piecewise linear basis functions for triangular elements,
% and piecewise isoparametric bilinear basis functions for quadrilateral
% elements.
%
%
%
clear
%
% Read the nodal coordinate data file.
%
load coordinates.dat;
%
% Read the triangular element data file.
%
% load elements3.dat;
%
% Read the quadrilateral element data file.
%
load elements4.dat;
%
% Read the Neumann boundary condition data file.
% I THINK the purpose of the EVAL command is to create an empty NEUMANN array
% if no Neumann file is found.
%
eval ( 'load neumann.dat;', 'neumann=[];' );
%
% Read the Dirichlet boundary condition data file.
%
load dirichlet.dat;
A = sparse ( size(coordinates,1), size(coordinates,1) );
b = sparse ( size(coordinates,1), 1 );
%
% Assembly.
%{
for j = 1 : size(elements3,1)
A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...
+ stima3(coordinates(elements3(j,:),:));
end
%}
%%{
for j = 1 : size(elements4,1)
A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) ...
+ stima4(coordinates(elements4(j,:),:));
end
%%}
% Volume Forces.
%
% from the center of each element to Nodes
% Notice that the result of f here means qv/K instead of qv
%{
for j = 1 : size(elements3,1)
b(elements3(j,:)) = b(elements3(j,:)) ...
+ det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...
f(sum(coordinates(elements3(j,:),:))/3)/6;
end
%}
%%{
for j = 1 : size(elements4,1)
b(elements4(j,:)) = b(elements4(j,:)) ...
+ det([1,1,1; coordinates(elements4(j,1:3),:)'] ) * ...
f(sum(coordinates(elements4(j,:),:))/4)/4;
end
%%}
% Neumann conditions.
%
if ( ~isempty(neumann) )
for j = 1 : size(neumann,1)
b(neumann(j,:)) = b(neumann(j,:)) + ...
norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...
g(sum(coordinates(neumann(j,:),:))/2)/2;
end
end
%
% Determine which nodes are associated with Dirichlet conditions.
% Assign the corresponding entries of U, and adjust the right hand side.
%
u = sparse ( size(coordinates,1), 1 );
BoundNodes = unique ( dirichlet );
u(BoundNodes) = u_d ( coordinates(BoundNodes,:) );
b = b - A * u;
%
% Compute the solution by solving A * U = B for the remaining unknown values of U.
%
FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );
u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);
%
% Graphic representation.
%
% show ( elements3, elements4, coordinates, full ( u ) );
trisurf ( elements4, coordinates(:,1), coordinates(:,2), u')
shading interp
xlabel('x')
--------------------------------------------------------------------------------------
结果如下，和三节点完全相同：