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
numberOfPDE = 1; pdem = createpde(numberOfPDE);
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. gdm = [3 4 0 width width 0 0 0 height height]'; g = decsg(gdm, 'S1', ('S1')'); % Convert the DECSG geometry into a geometry object % on doing so it is appended to the PDEModel geometryFromEdges(pdem,g); % Plot the geometry and display the edge labels for use in the boundary % condition definition. figure; pdegplot(pdem,'EdgeLabels','on'); axis([-.1 1.1 -.1 1.1]); title 'Geometry With Edge Labels Displayed';
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 = @(~,state) 2*hCoeff + 2*emiss*stefanBoltz*state.u.^3; f = 2*hCoeff*ta + 2*emiss*stefanBoltz*ta^4; d = thick*rho*specificHeat; specifyCoefficients(pdem,'m',0,'d',0,'c',c,'a',a,'f',f);
The bottom edge of the plate is set to 1000 degrees-Kelvin.
The boundary conditions are defined below. 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,
The initial guess is defined below.
Create the triangular mesh on the square with approximately ten elements in each direction.
hmax = .1; % element size msh = generateMesh(pdem,'Hmax',hmax); figure; pdeplot(pdem); axis equal title 'Plate With Triangular Element Mesh' xlabel 'X-coordinate, meters' ylabel 'Y-coordinate, meters'
Because the a and f coefficients are functions of temperature (due to the radiation boundary conditions),
solvepde automatically picks the nonlinear solver to obtain the solution.
R = solvepde(pdem); u = R.NodalSolution; figure; pdeplot(pdem,'XYData',u,'Contour','on','ColorMap','jet'); title 'Temperature In The Plate, Steady State Solution' xlabel 'X-coordinate, meters' ylabel 'Y-coordinate, meters' axis equal p = msh.Nodes; plotAlongY(p,u,0); title 'Temperature As a Function of the Y-Coordinate' xlabel 'Y-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 = 449.2 degrees-K
Include the "d" coefficient
specifyCoefficients(pdem,'m',0,'d',d,'c',c,'a',a,'f',f); 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; % Set the initial temperature on the bottom edge E1 to the value % of the constant BC, 1000 K setInitialConditions(pdem,1000,'edge',1); % Set solver options pdem.SolverOptions.RelativeTolerance = 1.0e-3; pdem.SolverOptions.AbsoluteTolerance = 1.0e-4; % |solvepde| automatically picks the parabolic solver to obtain the solution. R = solvepde(pdem,tlist); u = R.NodalSolution; figure; plot(tlist,u(3, :)); grid on title 'Temperature Along the Top Edge of the Plate as a Function of Time' xlabel 'Time, seconds' ylabel 'Temperature, degrees-Kelvin' % figure; pdeplot(pdem,'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' axis equal; % fprintf('\nTemperature at the top edge of the plate(t = %5.1f secs) = %5.1f degrees-K\n', ... tlist(1,end), u(4,end));
Temperature at the top edge of the plate(t = 5000.0 secs) = 441.3 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.