Partial Differential Equation Toolbox
This example shows how to perform a heat transfer analysis of a thin plate using the Partial Differential Equation Toolbox™.
The plate is square and the temperature is fixed along the bottom edge. No heat is transferred from the other three edges (i.e. they are insulated). Heat is transferred from both the top and bottom faces of the plate by convection and radiation. Because radiation is included, the problem is nonlinear. One of the purposes of this example is to show how to handle nonlinearities in PDE problems.
Both a steady state and a transient analysis are performed. In a steady state analysis we are interested in the final temperature at different points in the plate after it has reached an equilibrium state. In a transient analysis we are interested in the temperature in the plate as a function of time. One question that can be answered by this transient analysis is how long does it take for the plate to reach an equilibrium temperature.
The plate has planar dimensions one meter by one meter and is 1 cm thick. Because the plate is relatively thin compared with the planar dimensions, the temperature can be assumed constant in the thickness direction; the resulting problem is 2D.
Convection and radiation heat transfer are assumed to take place between the two faces of the plate and a specified ambient temperature.
The amount of heat transferred from each plate face per unit area due to convection is defined as
where is the ambient temperature, is the temperature at a particular x and y location on the plate surface, and is a specified convection coefficient.
The amount of heat transferred from each plate face per unit area due to radiation is defined as
where is the emissivity of the face and is the Stefan-Boltzmann constant. Because the heat transferred due to radiation is proportional to the fourth power of the surface temperature, the problem is nonlinear.
The PDE describing the temperature in this thin plate is
where is the material density, is the specific heat, is the plate thickness, and the factors of two account for the heat transfer from both plate faces.
It is convenient to rewrite this equation in the form expected by PDE Toolbox
The plate is composed of copper which has the following properties
k = 400; % thermal conductivity of copper, W/(m-K) rho = 8960; % density of copper, kg/m^3 specificHeat = 386; % specific heat of copper, J/(kg-K) thick = .01; % plate thickness in meters stefanBoltz = 5.670373e-8; % Stefan-Boltzmann constant, W/(m^2-K^4) hCoeff = 1; % Convection coefficient, W/(m^2-K) % The ambient temperature is assumed to be 300 degrees-Kelvin. ta = 300; emiss = .5; % emissivity of the plate surface
The expressions for the coefficients required by PDE Toolbox can easily be identified by comparing the equation above with the scalar parabolic equation in the PDE Toolbox documentation.
c = thick*k; % Because of the radiation boundary condition, the "a" coefficient % is a function of the temperature, u. It is defined as a MATLAB % expression so it can be evaluated for different values of u % during the analysis. a = sprintf('2*%g + 2*%g*%g*u.^3', hCoeff, emiss, stefanBoltz) f = 2*hCoeff*ta + 2*emiss*stefanBoltz*ta^4; d = thick*rho*specificHeat;
a = 2*1 + 2*0.5*5.67037e-08*u.^3
For a square, the geometry and mesh are easily defined as shown below.
width = 1; height = 1; % define the square by giving the 4 x-locations followed by the 4 % y-locations of the corners. gdmTrans = [3 4 0 width width 0 0 0 height height]; g = decsg(gdmTrans', 'S1', ('S1')'); % Create the triangular mesh on the square with approximately % ten elements in each direction. hmax = .1; % element size [p, e, t] = initmesh(g, 'Hmax', hmax); pdeplot(p,e,t); title 'Plate With Triangular Element Mesh' xlabel 'X-coordinate, meters' ylabel 'Y-coordinate, meters'
The bottom edge of the plate is set to 1000 degrees-Kelvin.
The boundary conditions are defined in the function below which is referred to in the PDE Toolbox documentation as a boundary file. Three of the plate edges are insulated. Because a Neumann boundary condition equal zero is the default in the finite element formulation, the boundary conditions on these edges do not need to be set explicitly. A Dirichlet condition is set on all nodes on the bottom edge, edge 1,
b=@boundaryFileThinPlate; type boundaryFileThinPlate
function [ q, g, h, r ] = boundaryFileThinPlate( p, e, u, time ) %BOUNDARYFILETHINPLATE Boundary conditions for heatTransferThinPlateExample % [ q, g, h, r ] = BOUNDARYFILETHINPLATE( p, e, u, time ) returns the % Neumann BC (q, g) and Dirichlet BC (h, r) matrices for the % heatTransferThinPlateExample example. % p is the point matrix returned from INITMESH % e is the edge matrix returned from INITMESH % u is the solution vector (used only for nonlinear cases) % time (used only for parabolic and hyperbolic cases) % % See also PDEBOUND, INITMESH % Copyright 2012 The MathWorks, Inc. % $Revision: 188.8.131.52 $ $Date: 2012/08/21 01:12:25 $ N = 1; ne = size(e,2); q = zeros(N^2, ne); g = zeros(N, ne); h = zeros(N^2, 2*ne); r = zeros(N, 2*ne); for i=1:ne ei = e(5,i); if(ei == 1) % Set the temperature at both vertices on the edge h(1,i) = 1; h(1,i+ne) = 1; r(1,i) = 1000; r(1,i+ne) = 1000; end end end
Because the a and f coefficients are functions of temperature (due to the radiation boundary conditions), the nonlinear solver pdenonlin must be used to obtain the solution.
u = pdenonlin(b,p,e,t,c,a,f, 'jacobian', 'lumped'); figure; pdeplot(p, e, t, 'xydata', u, 'contour', 'on', 'colormap', 'jet'); title 'Temperature In The Plate, Steady State Solution' xlabel 'X-coordinate, meters' ylabel 'Y-coordinate, meters' plotAlongY(p, u, 0); title 'Temperature As a Function of the Y-Coordinate' xlabel 'X-coordinate, meters' ylabel 'Temperature, degrees-Kelvin' fprintf('Temperature at the top edge of the plate = %5.1f degrees-K\n', ... u(4));
Temperature at the top edge of the plate = 448.9 degrees-K
endTime = 5000; tlist = 0:50:endTime; numNodes = size(p,2); % Set the initial temperature of all nodes to ambient, 300 K u0(1:numNodes) = 300; % Find all nodes along the bottom edge and set their initial temperature % to the value of the constant BC, 1000 K nodesY0 = abs(p(2,:)) < 1.0e-5; u0(nodesY0) = 1000; rtol = 1.0e-3; atol = 1.0e-4; % The transient solver parabolic automatically handles both linear % and nonlinear problems, such as this one. u = parabolic(u0, tlist, b,p,e,t,c,a,f,d,rtol,atol); figure; plot(tlist, u(3, :)); grid; title 'Temperature Along the Top Edge of the Plate as a Function of Time' xlabel 'Time, seconds' ylabel 'Temperature, degrees-Kelvin' % figure; pdeplot(p, e, t, 'xydata', u(:,end), 'contour', 'on', 'colormap', 'jet'); title(sprintf('Temperature In The Plate, Transient Solution( %d seconds)\n', ... tlist(1,end))); xlabel 'X-coordinate, meters' ylabel 'Y-coordinate, meters' % fprintf('\nTemperature at the top edge of the plate(t=%5.1f secs) = %5.1f degrees-K\n', ... tlist(1,end), u(4,end));
65 successful steps 0 failed attempts 95 function evaluations 1 partial derivatives 16 LU decompositions 94 solutions of linear systems Temperature at the top edge of the plate(t=5000.0 secs) = 447.2 degrees-K
As can be seen, the plots of temperature in the plate from the steady state and transient solution at the ending time are very close. That is, after around 5000 seconds, the transient solution has reached the steady state values. The temperatures from the two solutions at the top edge of the plate agree to within one percent.