Accelerating the pace of engineering and science

# Documentation

## Parabolic PDEs

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:

### Heat Equation for Metal Block with Cavity

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

$d\frac{\partial u}{\partial t}-\Delta u=0.$

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 t0. 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.

#### Using the PDE App

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

$\frac{\partial u}{\partial n}=0$

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

$d\frac{\partial u}{\partial t}-\nabla \cdot \left(c\nabla u\right)+au=f,$

with initial values u0 = u(t0) 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 u0 = 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.

#### Metal Block Using Command-Line Functions

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')


### Heat Distribution in a Radioactive Rod

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:

$\rho C\frac{\partial u}{\partial t}-\nabla \text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(k\nabla u\right)=f,$

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/m3, the thermal capacity is 500 Ws/kgºC, and the thermal conductivity is 40 W/mºC. The heat source is 20000 W/m3. 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/m2ºC. The heat flux at the left end is 5000 W/m2.

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

$r\rho C\frac{\partial u}{\partial t}-\frac{\partial }{\partial r}\left(kr\frac{\partial u}{\partial r}\right)-\frac{\partial }{\partial z}\left(kr\frac{\partial u}{\partial z}\right)=fr.$

The boundary conditions are:

• $\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(k\nabla u\right)$ = 5000 at the left end of the rod (Neumann condition). Since the generalized Neumann condition in Partial Differential Equation Toolbox software is $\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(k\nabla u\right)$ + qu = g, and c depends on r in this problem (c = kr), this boundary condition is expressed as $\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\nabla u\right)$ = 5000r.

• u = 100 at the right end of the rod (Dirichlet condition).

• $\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(k\nabla u\right)$ = 50(100-u) at the outer boundary (generalized Neumann condition). In Partial Differential Equation Toolbox software, this must be expressed as $\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\nabla u\right)$ + 50r  ·  u = 50r · 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 $\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\nabla u\right)$ here.

The initial value is u(t0) = 0.

#### Using the PDE App

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.