Main Content

Inhomogeneous Heat Equation on Square Domain

This example shows how to solve the heat equation with a source term.

The basic heat equation with a unit source term is

ut-Δu=1

This equation is solved on a square domain with a discontinuous initial condition and zero temperatures on the boundaries.

Create a transient thermal model.

thermalmodel = createpde("thermal","transient");

Create a square geometry centered at x = 0 and y = 0 with sides of length 2. Include a circle of radius 0.4 concentric with the square.

R1 = [3;4;-1;1;1;-1;-1;-1;1;1];
C1 = [1;0;0;0.4];
C1 = [C1;zeros(length(R1) - length(C1),1)];
gd = [R1,C1];
sf = 'R1+C1';
ns = char('R1','C1')';
g = decsg(gd,sf,ns);

Append the geometry to the model.

geometryFromEdges(thermalmodel,g);

Specify thermal properties of the material.

thermalProperties(thermalmodel,"ThermalConductivity",1,...
                               "MassDensity",1,...
                               "SpecificHeat",1);

Specify internal heat source.

internalHeatSource(thermalmodel,1);

Plot the geometry and display the edge labels for use in the boundary condition definition.

figure
pdegplot(thermalmodel,"EdgeLabels","on","FaceLabels","on")
axis([-1.1 1.1 -1.1 1.1]);
axis equal
title("Geometry With Edge and Subdomain Labels")

Figure contains an axes object. The axes object with title Geometry With Edge and Subdomain Labels contains 11 objects of type line, text.

Set zero temperatures on all four outer edges of the square.

thermalBC(thermalmodel,"Edge",1:4,"Temperature",0);

The discontinuous initial value is 1 inside the circle and zero outside. Specify zero initial temperature everywhere.

thermalIC(thermalmodel,0);

Specify non-zero initial temperature inside the circle (Face 2).

thermalIC(thermalmodel,1,"Face",2);

Generate and plot a mesh.

msh = generateMesh(thermalmodel);
figure;
pdemesh(thermalmodel); 
axis equal

Figure contains an axes object. The axes object contains 2 objects of type line.

Find the solution at 20 points in time between 0 and 0.1.

nframes = 20;
tlist = linspace(0,0.1,nframes);

thermalmodel.SolverOptions.ReportStatistics ='on';
result = solve(thermalmodel,tlist);
97 successful steps
0 failed attempts
196 function evaluations
1 partial derivatives
21 LU decompositions
195 solutions of linear systems
T = result.Temperature;

Plot the solution.

figure
Tmax = max(max(T));
Tmin = min(min(T));
for j = 1:nframes
    pdeplot(thermalmodel,"XYData",T(:,j),"ZData",T(:,j));
    axis([-1 1 -1 1 0 1]);
    Mv(j) = getframe;
end

Figure contains an axes object. The axes object contains an object of type patch.

To play the animation, use the movie(Mv,1) command.