以下面的问题为例：对于如图所示的平面传热问题，
若上端有给定的热流-2W/m2，即从下往上传输热量，结构下端有确定的温度100°，周围介质温度为20°，在两侧有换热，换热系数为α=100W/㎡/K，热导率200W/m/K,试分析其稳定温度场。
使用下面的程序包进行求解：
首先：
将区域划分为4个单元，各单元包含的节点在element3.dat中显示：
1 2 4 2 5 4 2 6 5 2 3 6
每个节点的横纵坐标在coordinates.dat文件中显示：
-10 0 0 0 10 0 -5 10 0 10 5 10
边界包含三个：
恒温边界的节点编号，在dirichelet.dat文件中显示：
1 2
2 3
热流边界的节点编号在neumann.dat文件中显示：
4 5
5 6
对流环热边界的节点编号在convective.dat文件中显示：
3 6
1 4
下面介绍使用到的程序：
主程序Heat_conduction_steady.m
View Code %***************************************************************************** % % 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 % % clearclose 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 coefficient tf=20;% ambient temperature q=-2;% Surface heat flux at Neumann boundary condition K=200;% Thermal conductivityA = 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);endelsefor j = 1 : size(elements3,1)A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...+ stima3(coordinates(elements3(j,:),:),K);endend% Volume Forces. % % from the center of each element to Nodes % Notice that the result of f here means qv/K instead of qvfor 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;endend% 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;endend % % 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.figuretrisurf ( elements3, coordinates(:,1), coordinates(:,2), (full ( u ))); view(2) colorbar shading interp xlabel('x')
热传导确定的单元刚度矩阵的子程序stima3.m
View Code 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
View Code 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) == 2temp2=find(temp1 == 0);if temp2 == 1ijm=23;elseif temp2 == 2ijm=31;elseif temp2 == 3ijm=12;endswitch ijmcase 12S=norm(vertices(2,:)-vertices(1,:));M=M+alpha*S/6*[2 1 0;1 2 0;0 0 0];case 23S=norm(vertices(3,:)-vertices(2,:));M=M+alpha*S/6*[0 0 0;0 2 1;0 1 2;];case 31S=norm(vertices(3,:)-vertices(1,:));M=M+alpha*S/6*[2 0 1;0 0 0;1 0 2;];end end end
由内生热导致的载荷项f.m：
View Code 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：
View Code 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：
View Code 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：
View Code 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程序自动将计算区域划分单元，且节点编号，边界条件的分配等。
View Code % 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 xminl=-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:Nyx=linspace(xmin(i1),xmax(i1),Nx);for i2=1:Nxk=k+1;Coord(k,:)=[x(i2),y(i1)]; end endsave coordinates.dat Coord -ascii% the element index k=0; vertices=zeros((Nx-1)*(Ny-1)*2,3); for i1=1:Ny-1for i2=1:Nx-1k=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(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
划分的单元如下：
计算的温度场分布如下：
