以下面的问题为例：对于如图所示的平面传热问题，

若上端有给定的热流-2W/m²，即从下往上传输热量，结构下端有确定的温度100°，周围介质温度为20°，在两侧有换热，换热系数为α=100W/㎡/K，热导率200W/m/K,试分析其稳定温度场。

 

使用下面的程序包进行求解：

首先：

将区域划分为4个单元，各单元包含的节点在element3.dat中显示：

1 2 42 5 42 6 52 3 6

每个节点的横纵坐标在coordinates.dat文件中显示：

-10 00 010 0-5 100 105 10

边界包含三个：

恒温边界的节点编号，在dirichelet.dat文件中显示：

1 2

2 3

热流边界的节点编号在neumann.dat文件中显示：

4 5

5 6

对流环热边界的节点编号在convective.dat文件中显示：

3 6

1 4

 

下面介绍使用到的程序：

主程序Heat_conduction_steady.m

%*****************************************************************************%%    The unknown state variable U(x,y) is assumed to satisfy%    Poisson's equation (steady heat conduction equation):%      -K(Uxx(x,y) + Uyy(x,y)) = qv(x,y) in Omega%    here qv means volumetric heat generation rate,K is thermocal%    conductivity%    %      clear  close all%  Read the nodal coordinate data file.  load coordinates.dat;%  Read the triangular element data file.  load elements3.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;%  Read the Convective heat transfer boundary condition data file.  eval ( 'load Convective.dat;', 'Convective=[];' );alpha=100;% Convective heat transfer coefficienttf=20;% ambient temperatureq=-2;% Surface heat flux at Neumann boundary conditionK=200;% Thermal conductivity  A = sparse ( size(coordinates,1), size(coordinates,1) );  b = sparse ( size(coordinates,1), 1 );%%  Assembly.  if ( ~isempty(Convective) )    for j = 1 : size(elements3,1)        A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...      + stima3(coordinates(elements3(j,:),:),K) ...      + stima3_1(coordinates(elements3(j,:),:),elements3(j,:),Convective,alpha);    end  else    for j = 1 : size(elements3,1)        A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...      + stima3(coordinates(elements3(j,:),:),K);    end  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%  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,q)/2;    end  end%  Convective heat transfer boundary condition.  if ( ~isempty(Convective) )    for j = 1 : size(Convective,1)      b(Convective(j,:)) = b(Convective(j,:)) + ...        norm(coordinates(Convective(j,1),:) - coordinates(Convective(j,2),:)) * ...        h(sum(coordinates(Convective(j,:),:))/2,tf,alpha)/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.figure   trisurf ( elements3, coordinates(:,1), coordinates(:,2), (full ( u )));view(2)colorbarshading interpxlabel('x')

热传导确定的单元刚度矩阵的子程序stima3.m

function M = stima3 ( vertices ,K)%*****************************************************************************80%%% STIMA3 determines the local stiffness matrix for a triangular element.%%  Parameters:%%    Input, %   real VERTICES(3,2), contains the 2D coordinates of the vertices.%   K: thermal conductivity%    Output, %   real M(3,3), 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 );M=M*K;

由对流换热确定的单元刚度矩阵的子程序stima3_1.m

function M = stima3 ( vertices,elements ,Convective,alpha)%*****************************************************************************80%% STIMA3 determines the local stiffness matrix for a triangular element.%  Parameters:%    Input, %   real VERTICES(3,2), contains the 2D coordinates of the vertices.%   elements(3,1), the node index of local element%   Convective(N,2), node index of each boundaryline of convective heat%   transfer boundary%   alpha, convective heat transfer coefficient%    Output, %   real M(3,3), the local stiffness matrix  for the element.%M=zeros(3,3);for i1=1:length(Convective)    temp1=ismember(elements,Convective(i1,:));    if sum(temp1) == 2        temp2=find(temp1 == 0);        if temp2 == 1            ijm=23;        elseif temp2 == 2            ijm=31;        elseif temp2 == 3            ijm=12;        end        switch ijm        case 12            S=norm(vertices(2,:)-vertices(1,:));            M=M+alpha*S/6*[            2 1 0;            1 2 0;            0 0 0            ];        case 23        S=norm(vertices(3,:)-vertices(2,:));            M=M+alpha*S/6*[            0 0 0;            0 2 1;            0 1 2;            ];        case 31            S=norm(vertices(3,:)-vertices(1,:));            M=M+alpha*S/6*[            2 0 1;            0 0 0;            1 0 2;            ];        end         endend

由内生热导致的载荷项f.m：

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) = zeros(n,1);

由热流边界计算的载荷项g.m：

function value = g ( u,q )%*****************************************************************************80%% G evaluates the outward normal values assigned at Neumann boundary conditions.%  Parameters:%    Input, %    real U(N,2), contains the 2D coordinates of N points.%    q: surface heat flux at Neumann boundary%    Output, %   VALUE(N), contains the value of outward normal at each point%    where a Neumann boundary condition is applied.%  value = q*ones ( size ( u, 1 ), 1 );

由对流换热导致的载荷项h.m：

function value = h ( u,tf,alpha)%****************************************************************************%%  evaluates the Convective heat transfer force.%  Parameters:%    Input, %    real U(N,2), contains the 2D coordinates of N points.%    tf: ambient temperature at Convective boundary%    alpha: convective heat transfer coefficient%    Output, %   VALUE(N), contains the value of outward normal at each point%    where a Neumann boundary condition is applied.%  value = alpha*tf*ones ( 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 = ones ( size ( u, 1 ), 1 )*100;

 

上述程序，所得结果为：

 

 下面使用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 nextxminl=-10;xmaxl=10;xminu=-5;xmaxu=5;ymin=0;ymax=10;Nx=13;Ny=13;y=linspace(ymin,ymax,Ny);xmin=linspace(xminl,xminu,Ny);xmax=linspace(xmaxl,xmaxu,Ny);k=0;for i1=1:Ny    x=linspace(xmin(i1),xmax(i1),Nx);    for i2=1:Nx        k=k+1;        Coord(k,:)=[x(i2),y(i1)];        endendsave coordinates.dat Coord -ascii% the element indexk=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];            endendsave elements3.dat vertices -ascii% The direchlet boundary condition (index of the two end nodes for each boundary line)boundary=zeros(Nx-1,2);temp1=1:Nx-1;temp2=2:Nx;temp3=1:Nx-1;boundary(temp3',:)=[temp1',temp2'];save dirichlet.dat boundary -ascii% The Neumann boundary condition (index of the two end nodes for each boundary line)boundary=zeros(Nx-1,2);temp1=(1:Nx-1)+(Ny-1)*Nx;temp2=(2:Nx)+(Ny-1)*Nx;temp3=1:Nx-1;boundary(temp3',:)=[temp1',temp2'];save neumann.dat boundary -ascii% The Convective heat transfer boundary condition (index of the two end nodes for each boundary line)boundary=zeros((Ny-1)*2,2);temp1=Nx*(1:Ny-1);temp2=temp1+Nx;temp3=1:Ny-1;boundary(temp3',:)=[temp1',temp2'];temp1=Nx*(Ny-1)+1:-Nx:Nx+1;temp2=temp1-Nx;temp3=temp3+Ny-1;boundary(temp3',:)=[temp1',temp2'];save Convective.dat boundary -ascii

划分的单元如下：

计算的温度场分布如下：