FEM计算2D瞬态热传导方程
本文程序，使用有限元计算2D瞬态热传导方程。
本程序只使用2D三角形单元求解，没有四边形单元求解。不过，也非常容易就可做到。
主程序为:
Heat_Conduction_transient.m
大部分子程序和这篇文章相同，但是多了下面的初始温度场的子程序。
u_init.m
现将各个程序列于下：
%function Heat_Conduction_transient( ) %*****************************************************************************80 % %% Applies the finite element method to the heat equation. % % The user supplies datafiles that specify the geometry of the region % and its arrangement into triangular elements, and the location and % type of the boundary conditions, which can be any mixture of Neumann and % Dirichlet, and the initial condition for the solution. % % Note that only uses triangular elements in this version. % % The unknown state variable U(x,y,t) is assumed to satisfy % the time dependent heat equation: % % dU(x,y,t)/dt = Uxx(x,y,t) + Uyy(x,y,t) + F(x,y,t) in Omega x [0,T] % % with initial conditions % % U(x,y,0) = U_INIT(x,y,0) % % with Dirichlet boundary conditions % % U(x,y,t) = U_D(x,y,t) on Gamma_D % % and Neumann boundary conditions on the outward normal derivative: % % Un(x,y) = G(x,y,t) on Gamma_N % % If Gamma designates the boundary of the region Omega, % then we presume that % % Gamma = Gamma_D + Gamma_N % % but the user is free to determine which boundary conditions to % apply. Note, however, that the problem will generally be singular % unless at least one Dirichlet boundary condition is specified. % % The code uses piecewise linear basis functions for triangular elements. % % The user is required to supply a number of data files and MATLAB % functions that specify the location of nodes, the grouping of nodes % into elements, the location and value of boundary conditions, and % the right hand side function in the heat equation. Note that the % fact that the geometry is completely up to the user means that % just about any two dimensional region can be handled, with arbitrary % shape, including holes and islands. % % % Local Parameters: % % Local, real DT, the size of a single time step. % % Local, integer NT, the number of time steps to take. % % Local, real T, the current time. % % Local, real T_FINAL, the final time. % % Local, real T_START, the initial time. % % 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. % There must be at least one Dirichlet boundary condition. % load dirichlet.dat; % % Determine the bound and free nodes. % BoundNodes = unique ( dirichlet ); FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes ); A = sparse ( size(coordinates,1), size(coordinates,1) ); B = sparse ( size(coordinates,1), size(coordinates,1) ); t_start = 0.0; t_final = 1; nt = 10;% nt means number of time intervals instead of timesteps t = t_start; dt = ( t_final - t_start ) / nt; U = zeros ( size(coordinates,1), nt+1 ); % % Assembly. % for j = 1 : size(elements3,1) A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ... + stima3(coordinates(elements3(j,:),:)); end for j = 1 : size(elements3,1) B(elements3(j,:),elements3(j,:)) = B(elements3(j,:),elements3(j,:)) ... + det ( [1,1,1; coordinates(elements3(j,:),:)' ] )... * [ 2, 1, 1; 1, 2, 1; 1, 1, 2 ] / 24; end % % Set the initial condition. % U(:,1) = u_init ( coordinates, t ); % % Given the solution at step I-1, compute the solution at step I. % for i = 1 : nt t = ( ( nt - i ) * t_start ... + ( i ) * t_final ) ... / nt; b = sparse ( size(coordinates,1), 1 ); u = sparse ( size(coordinates,1), 1 ); % % Account for the volume forces, evaluated at the new time T. % for j = 1 : size(elements3,1) b(elements3(j,:)) = b(elements3(j,:)) ... + det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ... dt * f ( sum(coordinates(elements3(j,:),:))/3, t ) / 6; end % % Account for the Neumann conditions, evaluated at the new time T. % for j = 1 : size(neumann,1) b(neumann(j,:)) = b(neumann(j,:)) + ... norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ... dt * g ( sum(coordinates(neumann(j,:),:))/2, t ) / 2; end % % Account for terms that involve the solution at the previous timestep. % b = b + B * U(:,i); % % Use the Dirichlet conditions, evaluated at the new time T, to eliminate the % known state variables. % u(BoundNodes) = u_d ( coordinates(BoundNodes,:), t ); b = b - ( dt * A + B ) * u; % % Compute the remaining unknowns by solving ( dt * A + B ) * U = b. % u(FreeNodes) = ( dt * A(FreeNodes,FreeNodes) + B(FreeNodes,FreeNodes) ) ... \ b(FreeNodes); U(:,i+1) = u; end show ( elements3, coordinates, U, t_start, t_final );
u_init.m
function value = u_init ( xy, t ) %*****************************************************************************80 % %% U_INIT sets the initial condition for the state variable. % % Discussion: % % The user must supply the appropriate routine for a given problem % % In many problems, the initial time is 0. However, the value of % T is passed, in case the user wishes to use this same routine to % evaluate, for instance, the exact solution. % % % Parameters: % % Input, real XY(N,M), contains the M-dimensional coordinates of N points. % N is probably the total number of points, and M is probably 2. % % Input, real T, the initial time. % % Output, VALUE(N), contains the value of the solution U(X,Y,T). % n = size ( xy, 1 ); value = zeros ( n, 1 ); middle = floor ( n / 2 ); value ( middle ) = 1.0;
u_d.m
View Code
function value = u_d ( u, t ) %*****************************************************************************80 % %% U_D evaluates the Dirichlet boundary conditions. % % Discussion: % % 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. % n = size ( u, 1 ); value =zeros(n,1);% 0.1 * ( u(:,1) + u(:,2) ); value(1:13)=1;
stima3.m
View Code
function M = stima3 ( vertices ) %*****************************************************************************80 % %% STIMA3 determines the local stiffness matrix for a triangular element. % % Discussion: % % Although this routine is intended for 2D usage, the same formulas % work for tetrahedral elements in 3D. The spatial dimension intended % is determined implicitly, from the spatial dimension of the vertices. % % Modified: % % 23 February 2004 % % Author: % % Jochen Alberty, Carsten Carstensen, Stefan Funken. % % Reference: % % Jochen Alberty, Carsten Carstensen, Stefan Funken, % Remarks Around 50 Lines of MATLAB: % Short Finite Element Implementation, % Numerical Algorithms, % Volume 20, pages 117-137, 1999. % % Parameters: % % Input, real VERTICES(1:(D+1),1:D), contains the D-dimensional % coordinates of the vertices. % % Output, real M(1:(D+1),1:(D+1)), 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 ); return end
show.m
View Code
function show ( elements3, coordinates, U, t_start, t_final ) %*****************************************************************************80 % %% SHOW displays the solution of the finite element computation. % % % Parameters: % % Input, integer ELEMENTS3(N3,3), the nodes that make up each triangle. % % Input, real COORDINATES(N,1:2), the coordinates of each node. % % Input, real U(N,NT+1), the solution, for each time step. % % Input, integer NT, the number of time steps. % % Input, real T_START, T_FINAL, the start and end times. % umin = min ( min ( U ) ); umax = max ( max ( U ) ); nt=size(U,2)-1; for i = 0 : nt % % Print the current time T to the command window, and in the plot title. % t = ( ( nt - i ) * t_start ... + ( i ) * t_final ) ... / nt; fprintf ( 1, 'T = %f\n', t ); % picture = trisurf ( elements3, coordinates(:,1), coordinates(:,2), ... U(:,i+1)', 'EdgeColor', 'interp', 'FaceColor', 'interp' ); % % We want all the plots to use the same Z scale. % axis ( [0 1 0 1 umin umax] ); % % Write some labels on the plot. % xlabel ( 'X axis' ); ylabel ( 'Y axis' ); s = sprintf ( 'T = %8f', t ); title ( s ); % drawnow end
f.m
View Code
function value = f ( xy, t ) %*****************************************************************************80 % %% F evaluates the right hand side of the heat equation. % % % This routine must be changed by the user to reflect a particular problem. % % % Parameters: % % Input, real XY(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 ( xy, 1 ); value(1:n) = 0;
g.m
View Code
function value = g ( u, t ) %*****************************************************************************80 % %% G evaluates the outward normal values assigned at Neumann boundary conditions. % % This routine must be changed by the user to reflect a particular problem. % % For this particular problem, we want to set the value of G(X,Y,T) to 1 % if X is 1, and to 0 otherwise. % % Parameters: % % Input, real U(N,M), contains the M-dimensional coordinates of N points. % % Input, real T, the current time. % % Output, VALUE(N), contains the value of the outward normal at each point % where a Neumann boundary condition is applied. % n = size ( u, 1 ); value = zeros ( n, 1 );
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 xmin=0;xmax=1; ymin=0;ymax=1; Nx=13;Ny=13; x=linspace(xmin,xmax,Nx); y=linspace(ymin,ymax,Ny); k=0; for i1=1:Ny 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(2*(Nx+Ny-2),2); temp1=1:Nx-1; temp2=2:Nx; temp3=1:Nx-1; boundary(temp3',:)=[temp1',temp2']; temp1=Nx*(1:Ny-1); temp2=temp1+Nx; temp3=temp3+Nx-1; boundary(temp3',:)=[temp1',temp2']; temp1=Ny*Nx:-1:Ny*Nx-Nx+2; temp2=temp1-1; temp3=temp3+Nx-1; boundary(temp3',:)=[temp1',temp2']; temp1=Nx*(Ny-1)+1:-Nx:Nx+1; temp2=temp1-Nx; temp3=temp3+Nx-1; boundary(temp3',:)=[temp1',temp2']; save dirichlet.dat boundary -ascii
