## Documentation Center |

This section describes the solution of some parabolic PDE problems. The problems are solved using both the PDE app and the command-line functions. The topics include:

On this page… |
---|

This example shows how to solve a heat equation that describes the diffusion of heat in a body. The heat equation has the form:

Consider a metal block containing rectangular crack or cavity.
The left side of the block is heated to 100 degrees centigrade. At
the right side of the metal block, heat is flowing from the block
to the surrounding air at a constant rate. All the other block boundaries
are isolated. This leads to the following set of boundary conditions
(when proper scaling of *t* is chosen):

*u*= 100 on the left side (Dirichlet condition)∂

*u*/∂*n*= –10 on the right side (Neumann condition)∂

*u*/∂*n*= 0 on all other boundaries (Neumann condition)

Also, for the heat equation we need an initial value: the temperature
in the metal block at the starting time *t*_{0}.
In this case, the temperature of the block is 0 degrees at the time
we start applying heat.

Finally, to complete the problem formulation, we specify that the starting time is 0 and that we want to study the heat distribution during the first five seconds.

Once you have started the PDE app and selected the **Generic
Scalar** mode, drawing the CSG model can be done very quickly:
Draw a rectangle (R1) with the corners in `x = [-0.5 0.5 0.5
-0.5]` and `y = [-0.8 -0.8 0.8 0.8]`. Draw
another rectangle (R2) to represent the rectangular cavity. Its corners
should have the coordinates `x = [-0.05 0.05 0.05 -0.05] `and` y
= [-0.4 -0.4 0.4 0.4].` To assist in drawing the narrow rectangle
representing the cavity, open the Grid Spacing dialog box from the **Options** and
enter *x*-axis extra ticks at `-0.05` and `0.05`.
Then turn on the grid and the "snap-to-grid" feature.
A rectangular cavity with the correct dimensions is then easy to draw.

The CSG model of the metal block is now simply expressed as
the set formula `R1-R2`.

Leave the draw mode and enter the boundary mode by clicking
the ∂Ω button, and continue by selecting boundaries and
specifying the boundary conditions. Using the **Select
All** option from the **Edit** menu
and then defining the Neumann condition

for all boundaries first is a good idea since that leaves only the leftmost and rightmost boundaries to define individually.

The next step is to open the PDE Specification dialog box and enter the PDE coefficients.

The generic parabolic PDE that Partial Differential Equation Toolbox™ functions solve is

with initial values *u*_{0} = *u*(*t*_{0})
and the times at which to compute a solution specified in the array `tlist.`

For this case, you have *d* = 1, *c* =
1, *a* = 0, and *f* = 0.

Initialize the mesh by clicking the Δ button. If you want,
you can refine the mesh by clicking the **Refine** button.

The initial values *u*0 = 0, and the list of
times is entered as the MATLAB^{®} array [`0:0.5:5]`.
They are entered into the Solve Parameters dialog box, which is accessed
by selecting **Parameters** from the **Solve** menu.

The problem can now be solved. Pressing the **=** button
solves the heat equation at 11 different times from 0 to 5 seconds.
By default, an interpolated plot of the solution, i.e., the heat distribution,
at the end of the time span is displayed.

A more interesting way to visualize the dynamics of the heat
distribution process is to *animate* the solution.
To start an animation, check the **Animation** check
box in the Plot selection dialog box. Also, select the colormap` hot`.
Click the **Plot** button to start a recording
of the solution plots in a separate figure window. The recorded animation
is then "played" five times.

The temperature in the block rises very quickly. To improve
the animation and focus on the first second, try to change the list
of times to the MATLAB expression `logspace(-2,0.5,20)`.

Also, try to change the heat capacity coefficient `d` and
the heat flow at the rightmost boundary to see how they affect the
heat distribution.

This example shows how to solve for the heat distribution in the metal block with cavity using command-line functions. First, create geometry and boundary condition files. The files used here were created using the PDE app. `crackg.m` describes the geometry of the metal block, and `crackb.m` describes the boundary conditions.

To create an initial mesh, call `initmesh`:

```
[p,e,t]=initmesh('crackg');
```

The heat equation can now be solved using the `parabolic` function. The generic parabolic PDE that `parabolic` solves is

with initial value *u* 0 = *u* ( *t* 0), and the times at which to compute a solution specified in the array `tlist`. For this case, you have *d* = 1, *c* = 1, *a* = 0, and *f* = 0. The initial value *u* 0 = 0, and the list of times, `tlist`, is set to the array `0:0.5:5`.

d = 1; c = 1; a = 0; f = 0; u0 = 0; tlist = 0:0.5:5;

To compute the solution, call `parabolic`:

```
u = parabolic(u0,tlist,'crackb',p,e,t,c,a,f,d);
```

153 successful steps 0 failed attempts 308 function evaluations 1 partial derivatives 28 LU decompositions 307 solutions of linear systems

The solution `u` is a matrix with 11 columns, where each column corresponds to the solution at the 11 points in time 0, 0.5, . . . , 4.5, 5.0.

Plot the solution at *t* = 5.0 seconds using interpolated shading and a hidden mesh. Use the `hot` colormap:

pdeplot(p,e,t,'xydata',u(:,11),'mesh','off','colormap','hot')

This example shows how to solve a 3-D parabolic PDE problem by reducing the problem to 2-D using coordinate transformation. For a step-by-step command-line solution, see Heat Distribution in a Circular Cylindrical Rod.

Consider a cylindrical radioactive rod. At the left end, heat is continuously added. The right end is kept at a constant temperature. At the outer boundary, heat is exchanged with the surroundings by transfer. At the same time, heat is uniformly produced in the whole rod due to radioactive processes. Assume that the initial temperature is zero. This leads to the following problem:

where *ρ* is the density, *C* is
the rod's thermal capacity, *k* is the thermal conductivity,
and *f* is the radioactive heat source.

The density for this metal rod is 7800 kg/m^{3},
the thermal capacity is 500 Ws/kgºC, and the thermal conductivity
is 40 W/mºC. The heat source is 20000 W/m^{3}.
The temperature at the right end is 100 ºC. The surrounding temperature
at the outer boundary is 100 ºC, and the heat transfer coefficient
is 50 W/m^{2}ºC. The heat flux at the
left end is 5000 W/m^{2}.

But this is a cylindrical problem, so you need to transform the equation,
using the cylindrical coordinates *r*, *z*,
and θ. Due to symmetry, the solution is independent of θ,
so the transformed equation is

The boundary conditions are:

= 5000 at the left end of the rod (Neumann condition). Since the generalized Neumann condition in Partial Differential Equation Toolbox software is +

*qu*=*g*, and*c*depends on*r*in this problem (*c*=*kr*), this boundary condition is expressed as = 5000*r*.*u*= 100 at the right end of the rod (Dirichlet condition).= 50(100-

*u*) at the outer boundary (generalized Neumann condition). In Partial Differential Equation Toolbox software, this must be expressed as + 50*r*·*u*= 50*r*· 100.The cylinder axis

*r*= 0 is not a boundary in the original problem, but in our 2-D treatment it has become one. We must give the*artificial*boundary condition here.

The initial value is *u*(*t*_{0})
= 0.

Solve this problem using the PDE app. Model the rod as a rectangle
with its base along the *x*-axis, and let the *x*-axis
be the *z* direction and the *y*-axis
be the *r* direction. A rectangle with corners in
(-1.5,0), (1.5,0), (1.5,0.2), and (-1.5,0.2) would then model a rod
with length 3 and radius 0.2.

Enter the boundary conditions by double-clicking the boundaries
to open the Boundary Condition dialog box. For the left end, use Neumann
conditions with `0` for `q` and `5000*y` for `g`.
For the right end, use Dirichlet conditions with `1` for `h` and `100` for `r`.
For the outer boundary, use Neumann conditions with `50*y` for `q` and `50*y*100` for `g`.
For the axis, use Neumann conditions with `0` for `q` and `g`.

Enter the coefficients into the PDE Specification dialog box: `c` is `40*y`, `a` is
zero, `d` is `7800*500*y`, and `f` is `20000*y`.

Animate the solution over a span of 20000 seconds (computing
the solution every 1000 seconds). We can see how heat flows in over
the right and outer boundaries as long as *u * <
100, and out when *u * > 100. You can also open
the PDE Specification dialog box, and change the PDE type to **Elliptic**.
This shows the solution when *u* does not depend
on time, i.e., the steady state solution. The profound effect of cooling
on the outer boundary can be demonstrated by setting the heat transfer
coefficient to zero.

Was this topic helpful?