Solve Partial Differential Equation of Nonlinear Heat Transfer

Heat Transfer Equations for the Plate

The plate has planar dimensions 1 m by 1 m and is 1 cm thick. Because the plate is relatively thin compared with the planar dimensions, the temperature can be assumed to be constant in the thickness direction, and the resulting problem is 2-D.

Assume that convection and radiation heat transfers take place between the two faces of the plate and the environment with a specified ambient temperature.

The amount of heat transferred from each plate face per unit area due to convection is defined as

$$Q_c=h_c(T-T_a)$$,

where $$T_a$$ is the ambient temperature, $$T$$ is the temperature at a particular $$x$$ and $$y$$ location on the plate surface, and $$h_c$$ is a specified convection coefficient.

The amount of heat transferred from each plate face per unit area due to radiation is defined as

$$Q_r=ϵ\sigma (T^4-T_a^4)$$,

where $$ϵ$$ is the emissivity of the face and $$\sigma$$ 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 $$\rho$$ is the material density of the plate, $$C_p$$ is its specific heat, $$t_z$$ is its plate thickness, $$k$$ is its thermal conductivity, and the factors of two account for the heat transfer from both of its faces.

Define PDE Parameters

Set up the PDE problem by defining the values of the parameters. The plate is composed of copper, which has the following properties.

kThermal = 400;             % thermal conductivity of copper, W/(m-K)
rhoCopper = 8960;           % density of copper, kg/m^3
specificHeat = 386;         % specific heat of copper, J/(kg-K)
thick = 0.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)
tAmbient = 300;             % the ambient temperature
emiss = 0.5;                % emissivity of the plate surface

Extract PDE Coefficients

Define the PDE in symbolic form with the plate temperature as a dependent variable $$T(t,x,y)$$.

syms T(t,x,y)
syms eps sig tz hc Ta rho Cp k
Qc = hc*(T - Ta);
Qr = eps*sig*(T^4 - Ta^4);
pdeeq = (rho*Cp*tz*diff(T,t) - k*tz*laplacian(T,[x,y]) + 2*Qc + 2*Qr)

pdeeq(t, x, y) = 
$$2\hspace{0.17em}\mathrm{eps}\hspace{0.17em}\mathrm{sig}\hspace{0.17em}\left({T\left(t,x,y\right)}^4-{\mathrm{Ta}}^4\right)-k\hspace{0.17em}\mathrm{tz}\hspace{0.17em}\left(\frac{{\partial }^2}{\partial x^2}\mathrm{ }T\left(t,x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }T\left(t,x,y\right)\right)-2\hspace{0.17em}\mathrm{hc}\hspace{0.17em}\left(\mathrm{Ta}-T\left(t,x,y\right)\right)+\mathrm{Cp}\hspace{0.17em}\rho \hspace{0.17em}\mathrm{tz}\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }T\left(t,x,y\right)$$

Now, create coefficients to use as inputs in the PDE model as required by Partial Differential Equation Toolbox. To do this, first extract the coefficients of the symbolic PDE as a structure of symbolic expressions using the pdeCoefficients function.

symCoeffs = pdeCoefficients(pdeeq,T,'Symbolic',true)

symCoeffs = struct with fields:
    m: 0
    a: 2*hc + 2*eps*sig*T(t, x, y)^3
    c: [2x2 sym]
    f: 2*eps*sig*Ta^4 + 2*hc*Ta
    d: Cp*rho*tz

Next, substitute the symbolic variables that represent the PDE parameters with their numeric values.

symVars = [eps sig tz hc Ta rho Cp k];
symVals = [emiss stefanBoltz thick hCoeff tAmbient rhoCopper specificHeat kThermal];
symCoeffs = subs(symCoeffs,symVars,symVals);

Finally, since the fields symCoeffs are symbolic objects, use the pdeCoefficientsToDouble function to convert the coefficients to the double data type, which makes them valid inputs for Partial Differential Equation Toolbox.

coeffs = pdeCoefficientsToDouble(symCoeffs)

coeffs = struct with fields:
    a: @makeCoefficient/coefficientFunction
    c: [4x1 double]
    m: 0
    d: 3.4586e+04
    f: 1.0593e+03

Specify PDE Model, Geometry, and Coefficients

Now, using Partial Differential Equation Toolbox, solve the PDE problem using finite element analysis based on these coefficients.

First, create the PDE model with a single dependent variable.

numberOfPDE = 1;
model = createpde(numberOfPDE);

Specify the geometry for the PDE model—in this case, the dimension of the square.

width = 1; 
height = 1;

Define a geometry description matrix. Create the square geometry using the decsg (Partial Differential Equation Toolbox) function. For a rectangular geometry, the first row contains 3, and the second row contains 4. The next four rows contain the $$x$$-coordinates of the starting points of the edges, and the four rows after that contain the $$y$$-coordinates of the starting points of the edges.

gdm = [3 4 0 width width 0 0 0 height height]';
g = decsg(gdm,'S1',('S1')');

Convert the DECSG geometry into a geometry object and include it in the PDE model.

geometryFromEdges(model,g);

Plot the geometry and display the edge labels.

figure; 
pdegplot(model,'EdgeLabels','on'); 
axis([-0.1 1.1 -0.1 1.1]);
title('Geometry with Edge Labels Displayed');

Next, create the triangular mesh in the PDE model with a mesh size of approximately 0.1 in each direction.

hmax = 0.1; % element size
msh = generateMesh(model,'Hmax',hmax);
figure; 
pdeplot(model); 
axis equal
title('Plate with Triangular Element Mesh');
xlabel('X-coordinate, meters');
ylabel('Y-coordinate, meters');

Specify the coefficients in the PDE model.

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

Find Transient Solution

Apply the boundary conditions. Three of the plate edges are insulated. Because a Neumann boundary condition equal to zero is the default in the finite element formulation, you do not need to set the boundary conditions on these edges. The bottom edge of the plate is fixed at 1000 K. Specify this using a Dirichlet condition on all nodes on the bottom edge (edge E1).

applyBoundaryCondition(model,'dirichlet','edge',1,'u',1000);

Specify the initial temperature to be 0 K everywhere, except at the bottom edge. Set the initial temperature on the bottom edge E1 to the value of the constant boundary condition, 1000 K.

setInitialConditions(model,0);
setInitialConditions(model,1000,'edge',1);

Define the time domain to find the transient solution of the PDE problem.

endTime = 10000;
tlist = 0:50:endTime;

Set the tolerance of the solver options.

model.SolverOptions.RelativeTolerance = 1.0e-3; 
model.SolverOptions.AbsoluteTolerance = 1.0e-4;

Solve the problem using solvepde. Plot the temperature along the top edge of the plate as a function of time.

R = solvepde(model,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 (s)'
ylabel 'Temperature (K)'

Based on the plot, the transient solution starts to reach its steady state value after 6000 seconds. The equilibrium temperature of the top edge approaches 450 K after 6000 seconds.

Show the temperature profile of the plate after 10,000 seconds.

figure;
pdeplot(model,'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;

Show the temperature at the top edge at 10,000 seconds.

u(3,end)

ans = 
449.7914