以下面的問題為例：對于如圖所示的平面傳熱問題，

若上端有給定的熱流-2W/m²，即從下往上傳輸熱量，結構下端有確定的溫度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

劃分的單元如下：

計算的溫度場分布如下：

?

轉載于:https://www.cnblogs.com/heaventian/archive/2012/12/04/2801606.html

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的包含对流环热，热流边界，等温边界的稳态热传导方程的FEM求解。的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。