包含对流环热，热流边界，等温边界的稳态热传导方程的FEM求解。
以下面的问题为例：对于如图所示的平面传热问题，
若上端有给定的热流-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 % % 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 coefficient tf=20;% ambient temperature q=-2;% Surface heat flux at Neumann boundary condition K=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) 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) == 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 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:Ny x=linspace(xmin(i1),xmax(i1),Nx); 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(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
划分的单元如下：
计算的温度场分布如下：
