工程建设科学技术奖申报网站,网站外链建设与文章发布规范,想要网站推广页面,免费加盟游戏代理以下面的问题为例#xff1a;对于如图所示的平面传热问题#xff0c; 若上端有给定的热流-2W/m2#xff0c;即从下往上传输热量#xff0c;结构下端有确定的温度100#xff0c;周围介质温度为20#xff0c;在两侧有换热#xff0c;换热系数为α100W/㎡/K#xff0c;热导…以下面的问题为例对于如图所示的平面传热问题 若上端有给定的热流-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 22 3 热流边界的节点编号在neumann.dat文件中显示 4 55 6 对流环热边界的节点编号在convective.dat文件中显示 3 61 4 下面介绍使用到的程序 主程序Heat_conduction_steady.m View Code %*****************************************************************************
%
% The unknown state variable U(x,y) is assumed to satisfy
% Poissons 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[]; );
alpha100;% Convective heat transfer coefficient
tf20;% ambient temperature
q-2;% Surface heat flux at Neumann boundary condition
K200;% 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,d1); vertices ] \ [ zeros(1,d); eye(d) ];
M det ( [ ones(1,d1); vertices ] ) * D_eta * D_eta / prod ( 1:d );
MM*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.
%
Mzeros(3,3);
for i11:length(Convective)temp1ismember(elements,Convective(i1,:));if sum(temp1) 2temp2find(temp1 0);if temp2 1ijm23;elseif temp2 2ijm31;elseif temp2 3ijm12;endswitch ijmcase 12Snorm(vertices(2,:)-vertices(1,:));MMalpha*S/6*[2 1 0;1 2 0;0 0 0];case 23Snorm(vertices(3,:)-vertices(2,:));MMalpha*S/6*[0 0 0;0 2 1;0 1 2;];case 31Snorm(vertices(3,:)-vertices(1,:));MMalpha*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 Laplaces 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 Laplaces% 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 Nx1~2*Nx (from left to right) for the next line above
% and then next
xminl-10;xmaxl10;
xminu-5;xmaxu5;
ymin0;ymax10;
Nx13;Ny13;
ylinspace(ymin,ymax,Ny);
xminlinspace(xminl,xminu,Ny);
xmaxlinspace(xmaxl,xmaxu,Ny);
k0;
for i11:Nyxlinspace(xmin(i1),xmax(i1),Nx);for i21:Nxkk1;Coord(k,:)[x(i2),y(i1)]; end
endsave coordinates.dat Coord -ascii% the element index
k0;
verticeszeros((Nx-1)*(Ny-1)*2,3);
for i11:Ny-1for i21:Nx-1kk1;ijm1i2(i1-1)*Nx;ijm2i21(i1-1)*Nx;ijm3i21i1*Nx;vertices(k,:)[ijm1,ijm2,ijm3];kk1;ijm1i21i1*Nx;ijm2i2i1*Nx;ijm3i2(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)
boundaryzeros(Nx-1,2);
temp11:Nx-1;
temp22:Nx;
temp31: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)
boundaryzeros(Nx-1,2);
temp1(1:Nx-1)(Ny-1)*Nx;
temp2(2:Nx)(Ny-1)*Nx;
temp31: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)
boundaryzeros((Ny-1)*2,2);
temp1Nx*(1:Ny-1);
temp2temp1Nx;
temp31:Ny-1;
boundary(temp3,:)[temp1,temp2];
temp1Nx*(Ny-1)1:-Nx:Nx1;
temp2temp1-Nx;
temp3temp3Ny-1;
boundary(temp3,:)[temp1,temp2];
save Convective.dat boundary -ascii 划分的单元如下 计算的温度场分布如下 转载于:https://www.cnblogs.com/heaventian/archive/2012/12/04/2801606.html
