PDEModels Best Practice

Conservative versus Non-Conservative Convection

Generally speaking a convection can be modeled in a conservative and non-conservative form. The conservative form is expressed as a ConservativeConvectionPDETerm and models . The alternative is to use a ConvectionPDETerm, that models . Note the different position of

Discontinuity of secondary unknown quantities

Sometimes, it is necessary to compute secondary quantities from the primary unknown variable. These secondary quantities are often proportional to the first derivative of the primary unknown variable. Unlike the primary unknown variable, secondary variables can become discontinuous across element boundaries.

The following example shows discontinuities that arise when taking the derivative of a dependent variable . Also, it shows how the secondary quantity becomes more accurate as the number of elements increases.

The equation to solve is a diffusion term, with coefficient , and a source term that varies quadratically in space. The equation is solved in a 1D domain from to . For this problem, an analytical solution exists.

The source term is represented by the following equation:

where is a constant and is the independent variable.

Load the finite element package:
Define variables:
Define the parameters to use:

Three different meshes will be used. The first will have 2 line elements, the second 4 and the last one 20.

Create the meshes:
Define the equation:

Two DirichletCondition will be specified at both ends of the line. At the left, a value of 1 will be specified, and at the right, a value of 0.

Define the Dirichlet conditions:
Solve the PDE for the different domains:

This problem has an analytical solution given below. Note that the solution if of 4th order. This means that the quadratic elements can used in the interpolating function can not perfectly interpolate the solution.

Define an analytical solution:
Visualize the different solutions, with element boundaries marked by red points:

We can see that in the primary quantity, the solution is more accurate when using more elements.

The interpolating function of the solution is quadratic. Now, when we compute the derivative of the solution the quadratic interpolating function effectively becomes linear. Additionally, it is in the nature of the finite elements that the linear elements will introduce discontinuity at the element boundaries. However, these discontinuities can be reduced by refining the mesh, as will be shown in the next.

Compute the derivative of the analytical expression:
Visualize the derivative of the analytical expression:
Compute the derivative of the numerical expressions:
Visualize the derivative of the 2-elements solution:
Compute the difference of the numerical derivative and the analytical derivative:
Evaluate the numerical derivative close by the discontinuity:
Visualize the derivative of the 4-elements solution:
Compute the difference of the numerical derivative and the analytical derivative:
Evaluate the derivative close at the discontinuity:

One can see a discontinuity at the intersection of each element, but the discontinuities are less pronounced than in the first case.

Visualize the derivative of the 20-elements solution:

With 20 elements, we no longer sees discontinuities.

Compute the difference of the numerical derivative and the analytical derivative:
Evaluate the derivative close at the discontinuity:
Plot how the error decays as more elements are used to approximate the derivative of the solution:

Dirac delta functions

Regularized delta functions are used to implement monopole sources in various fields of physics, such as acoustics.

To make use of a monopole source, the monopole source is located at . Then the monopole source term may be written as:

where is a Dirac delta function at the source location .

The Dirac delta function, however, poses a problem in numerical simulations as it can not be resolved in the discretized spatial domain. This is because the Dirac delta function is singular at the source location . Hence, an approximation to the Dirac delta function is needed. The process of approximating the Dirac delta function is called regularization.

There are various regularized delta functions available. In this tutorial we choose:

where is the regularization parameter that controls the support of the regularized delta functions . Typically should have a size comparable to the mesh spacing .

Set up the regularized delta functions centered at :
Inspect the regularized delta function centered at . Note that the width of support equals to :

Modifying will change the width of the regularized delta function, however, for the spatial integral is always 1.

Integrate over the regularized delta function:

Smoothed step functions

In some examples, we use a smoothed step function to prescribe a time profile for a transient parameter, for example the heat flux or the surface temperature . The smoothed step function is defined as follows:

Here the minimum value and the maximum value the function can reach are denoted by and . The location of the step is controlled by and the smoothed steepness is controlled by .

Define a smoothed step function ::
Visualize the smoothed step function :

Solving Memory-Intensive PDEs

The finite element method documentation has a section on Solving Memory-Intensive PDEs that may alo be useful to know about in the context of PDE solving.


Geometry Units

The units of a geometry can be rescaled. This is explained in the ElementMesh generation tutorial.

Material Units

Should the units of the geometry be different from the material units, then the material units can be scaled.

Set up a material and scale the units:

Internally, all material data units are converted to "SIBase" units. As a consequence the default unit of length is "Meters". If the units of the geometry are also in meters then nothing needs to be changed. If the units of the geometry are not in meters then either the PDE and material properties need to be scaled to the units the geometry or the geometry needs to be scaled to "Meters". To scale the units of the PDE and material parameters the parameter "ScaleUnits" can be given. If not explicitly stated otherwise, examples in this notebook use the default "SIBase" units.

The procedure works in the same way when the material parameters are given as a Quantity.

Set up a material and scale the units:
Create a scaled heat transfer PDE from the parameters: