Solving Partial Differential Equations with Finite Elements

Introduction

The aim of this tutorial is to give an introductory overview of the finite element method (FEM) as it is implemented in NDSolve. The notebook introduces finite element method concepts for solving partial differential equations (PDEs). First, typical workflows are discussed. The setup of regions, boundary conditions and equations is followed by the solution of the PDE with NDSolve. The visualization and animation of the solution is then introduced, and some theoretical aspects of the finite element method are presented. Classical PDEs such as the Poisson and Heat equations are discussed. Coupled PDEs are also introduced with examples from structural mechanics and fluid dynamics.

Why Finite Elements?

Explicit closed-form solutions for partial differential equations (PDEs) are rarely available. The finite element method (FEM) is a technique to solve partial differential equations numerically. It is important for at least two reasons. First, the FEM is able to solve PDEs on almost any arbitrarily shaped region. Second, the method is well suited for use on a large class of PDEs. While it is almost always possible to conceive better methods for a specific PDE on a specific region, the finite element method performs quite well for a large class of PDEs. In summary, the finite element method is important since it can deal with:

What Is Needed for a Finite Element Analysis

To solve partial differential equations with the finite element method, three components are needed:

The Scope of the Finite Element Method as Implemented in NDSolve

The current version of the implementation of the finite element method supports the following features:

Regions

To perform a finite element analysis, a region over which the partial differential equation is to be solved needs to be defined. Rectangular regions may be specified using {v,vmin,vmax} for each of the spatial independent variables. The exact argument specification is detailed in NDSolveValue.

Solve a Poisson equation and visualize the solution over a rectangular region.
In[2]:=
Click for copyable input
Out[3]=

Regions of arbitrary shape may be specified using the notation varsΩ, where Ω is a region so that RegionQ[Ω] gives True.

One way to specify a nonrectangular region is by using Boolean predicates. A predicate is a function that returns either True or False. Since the term predicate has many meanings in mathematical literature, the terminology Boolean predicate is used to emphasize that the predicate returns a Boolean. To describe regions, Boolean predicates can be used as in functions like RegionPlot and RegionPlot3D, for example.

Region setup and visualization.
In[4]:=
Click for copyable input
Out[5]=

This tutorial concentrates on solving partial differential equations with the finite element method, without emphasis on the creation of regions and meshes. More detailed information on this topic can be found in "Element Mesh Generation".

Classical Partial Differential Equations

The Coefficient Form of Partial Differential Equations

What types of equations can be solved with the finite element method as implemented in NDSolve? Consider a single partial differential equation in :

The PDE is defined in . Here is the dependent variable for which a solution is sought. The coefficients , , and are scalars; , and are vectors; and is an × matrix.

What follows are some well-known PDEs and their corresponding coefficients. To illustrate the generality of (1), the components that are relevant to a specific equation are black, while the non-relevant components are gray.

The Laplace equation simply contains a diffusive term:

To model Poisson's equation, only a small modification is needed; add a load term :

Helmholtz's equation adds a reaction term :

Convection-diffusion-reaction type equations are another common class of PDEs. Compared to the previous examples, these have an additional convection term :

The PDEs considered so far are stationary, i.e. they have no time dependence. The heat equation adds time dependence to the Poisson equation. It has the following form:

Similarly, the wave equation is given as:

Equation (1) provides the components for modeling a range of different phenomena, since it provides spatial derivatives up to order 2.

Poissons Equation with Dirichlet Conditions

The Poisson equation is to be solved over a region with boundary conditions. First, a region needs to be defined where the equation will be solved. Then the equation and boundary conditions are defined. Finally, the equation is solved over the region.

One way to specify a region is by using Boolean predicates.

Region setup and visualization.
In[6]:=
Click for copyable input
Out[7]=

Next, define a Poisson operator.

PDE operator setup.
In[8]:=
Click for copyable input

The third step is to set up the boundary conditions. The first boundary condition enforces the dependent variable to a value of 0 wherever the evaluation of x08y10 yields True on the boundary of the region . In this case, this is the upper-left part of the boundary. This type of boundary condition is called a Dirichlet boundary condition. A second Dirichlet boundary condition specifies the value of the dependent variable whenever yields True on the boundary of region . Boundary conditions will be explained in further detail in a later section.

Dirichlet boundary conditions setup.
In[9]:=
Click for copyable input

In the solution step, NDSolve needs the PDE, the boundary conditions and the region.

Solve the PDE.
In[10]:=
Click for copyable input
Out[10]=

Internally, the region is converted to a finite element mesh and a finite element method is used to solve the equation automatically.

To visualize the solution, a contour plot can be used.

Visualize the solution.
In[11]:=
Click for copyable input
Out[11]=

Partial Differential Equations and Boundary Conditions

NDSolve and related functions allow for specifying three types of spatial boundary conditions: Dirichlet conditions, Neumann values and periodic boundary conditions.

Dirichlet boundary conditions prescribe a constraint on the dependent variable of value on some part of the boundary:

Generalized Neumann boundary values, also known as Robin values or Neumann boundary conditions, specify a value . The value prescribes a flux over the outward normal on some part of the boundary:

The finite element method is based on the weak form of the differential equation. This form is obtained by taking equation (1), multiplying it by a so-called test function , and integrating over the region :

Integration by parts gives:

This process is done internally. The integrand in the boundary integral in (11) can be replaced with (9) and yields:

specifies the value of the boundary integral integrand (11) in the weak form and thus the name NeumannValue.

Instead of making use of integration by parts to obtain equation (11), the divergence theorem and Green's identities can also be used.

Periodic boundary conditions make the dependent variables behave according to a given relation between two distinct parts of the boundary. Other boundary conditions are conceivable, but currently not implemented.

In most cases, Dirichlet boundary conditions need not be associated with a particular equation. They can be specified independently of the equation. Dirichlet boundary conditions are specified as equations. Generalized Neumann values, on the other hand, are specified by giving a value, since the equation satisfied is implicit in the value. Neumann values are mathematically tied to the PDE.

For practical reasons, in NDSolve and related functions, NeumannValue needs to be given as a part of the equation. Consider a system of coupled PDEs. One would like to be able to unambiguously specify any given Neumann value to any given single PDE of that system of PDEs. It is not possible to derive (unambiguously) from the Neumann value with which PDE equation the value should be associated. Making a NeumannValue part of a PDE equation solves this problem without ambiguity. DirichletCondition may also be given in a PDE equation as well.

{op[u[x1,]] ΓN,ΓD,ΓP}partial differential equation with differential operator op, Neumann boundary values ΓN, Dirichlet boundary condition equations ΓD and periodic boundary condition ΓP
{op1[u[x1,],v[x1,],] ΓN1,op2[u[x1,],v[x1,],] ΓN2,,ΓD,ΓP}system of coupled partial differential equation operators opj, Neumann boundary values ΓNj and Dirichlet and periodic boundary condition equations ΓD and ΓP

Specifying partial differential equations with boundary conditions.

DirichletCondition, NeumannValue and PeriodicBoundaryCondition all require a second argument that is a predicate describing the location on the boundary where the conditions/values are to be applied. Additionally, the PeriodicBoundaryCondition has a third argument specifying the relation between the two parts of the boundary.

A partial differential equation typically needs at least one Dirichlet boundary condition on some part of the region to be uniquely solvable. For this reason, Dirichlet boundary conditions are also called essential boundary conditions.

The use of generalized Neumann values arises from the boundary integral from the weak form, and so they are also commonly referred to as natural conditions.

A pure Neumann boundary equation is a generalized Neumann boundary equation (9) where :

If is also not specified, then a Neumann zero boundary value is automatically implied ():

Note that when the Neumann value is zero, the term drops out in (12). This is typically a case where there is no flux or stress at the boundary. When no boundary conditions are explicitly specified for a part of the region boundary , a Neumann zero value is assumed.

Poissons Equation with Neumann Values

Consider the solution of Poisson's equation with generalized Neumann boundary values.

Use that same region and the PDE operator.
In[12]:=
Click for copyable input

Poisson's equation is given as:

The corresponding Neumann boundary equation is:

Note how the normal derivative of the Neumann-type equation relates to the PDE. In the Neumann boundary equation, one specifies the values of and on the right-hand side. The left-hand side is implicitly specified through the PDE. If the coefficient in the PDE is given as an identity matrix, then the coefficient on the left-hand side of the Neumann equation is that same coefficient .

The Neumann boundary value is specified though the values and . Typically when , one speaks of a Neumann boundary value, and in the case , one speaks of a generalized Neumann or Robin boundary value. It is common to use the term Neumann boundary condition for both the generalized Neumann and Neumann boundary value.

Specify a generalized Neumann value and Dirichlet boundary condition.
In[14]:=
Click for copyable input
In[15]:=
Click for copyable input

It is important to observe that the constants in the term of the Neumann boundary equation exactly match the constants in from the weak form of the PDE (11). Because these terms are required to match, NeumannValue only requires the user to enter the right-hand side of the boundary equation. For example, to enter the boundary equation , it suffices to enter NeumannValue[1,]. The correct expression for the left-hand side of the boundary equation (i.e. ) is implied by the term in the weak form of the PDE (12). This close relation of the Neumann boundary value and the PDE is further accommodated by making the NeumannValue part of the equation.

Solve the PDE.
In[16]:=
Click for copyable input
Out[16]=
Plot the solution.
In[17]:=
Click for copyable input
Out[17]=

In the previous examples, not all boundaries had boundary conditions explicitly attributed; for those parts of the boundaries, the Neumann zero condition is used.

Inspect the values of the solution at .
In[18]:=
Click for copyable input
Out[18]=

Note that the values beyond 8 are zero, as required by the Dirichlet condition.

Poissons Equation with a Periodic Boundary Condition

The purpose of PeriodicBoundaryCondition is to relate the values of a solution of a PDE at two distinct parts of the boundary. PeriodicBoundaryCondition[a+b u,pred,f] specifies a periodic boundary condition where the predicate pred specifies a target at which the solution should relate to values from a distinct source, and the function f maps coordinates from the target locations (blue) to the source locations (green) where the dependent variable relation is evaluated. The resulting values will be set as a condition at the target.

68.gif

An example is the best way to understand how PeriodicBoundaryCondition works.

Specify and visualize a region marked with a blue target and green source boundary.
In[19]:=
Click for copyable input
Out[20]=

The PDE in this example is a Laplacian with a load term of 1 unit active in a rectangle specified by . Everywhere else there the load is 0.

Specify a PDE.
In[21]:=
Click for copyable input

PeriodicBoundaryCondition allows either periodic or antiperiodic or relational boundary conditions to be modeled. In this example, a periodic boundary condition is presented. A periodic condition is specified when the first argument is the dependent variable, u, and requires that the dependent variable u at the predicate (blue, target) is to have the same values as at position (green, source). During the solution process, the dependent variable u is to be evaluated at the source. For this, a mapping from the predicate (blue, target) to (green, source) is needed. The function FindGeometricTransform can construct such a mapping from a list of coordinates.

Compute a mapping for the PeriodicBoundaryCondition using FindGeometricTransform.
In[22]:=
Click for copyable input

Investigate to see that the mapping maps a point from the green source to the blue target.

In[23]:=
Click for copyable input
Out[23]=

The PeriodicBoundaryCondition is then set up by giving the relation of the dependent variable, the target predicate and a mapping. In the periodic case, the relation is .

In[24]:=
Click for copyable input

Additionally, on both curved sides DirichletCondition is placed. Note that the DirichletCondition is excluded at such that there is no DirichletCondition on the target boundary. Having Dirichlet boundary conditions there would not be consistent with the solution mapped to these places.

Specify a DirichletCondition at .
In[25]:=
Click for copyable input
Solve the PDE.
In[26]:=
Click for copyable input
Visualize the solution.
In[27]:=
Click for copyable input
Out[27]=
Inspect that the solution at is the same as at .
In[28]:=
Click for copyable input
Out[29]=

An antiperiodic boundary condition can be realized by using the negative of the dependent variable in the PeriodicBoundaryCondition. In the next example, everything except the antiperiodic boundary condition remains the same.

Specify a antiperiodic PeriodicBoundaryCondition with .
In[30]:=
Click for copyable input
Solve the PDE.
In[31]:=
Click for copyable input
Visualize the solution.
In[32]:=
Click for copyable input
Out[32]=
Inspect the sum of the solution at the two parts of the boundary.
In[33]:=
Click for copyable input
Out[34]=

Partial Differential Equations with Variable Coefficients

In the following example, the PDE coefficients depend on space and time. Here indicates the temporal variable and the spatial variables:

The example of a shielded microstrip illustrates the effect of variable coefficients. It is an electrostatics application that can be modeled with a Poisson equation where the coefficient is the relative permittivity and depends on space:

A geometry is constructed where the boundary separates two materials. In each material, has a different value.

First, a boundary mesh is created. For more information about creating element meshes, please consult ToBoundaryMesh.

Create a boundary mesh.
In[35]:=
Click for copyable input
Out[36]=
Visualize the boundary mesh.
In[37]:=
Click for copyable input
Out[37]=

The conversion of the boundary mesh to an element mesh is done with ToElementMesh.

Convert the boundary mesh to a full mesh and visualize the resulting mesh.
In[38]:=
Click for copyable input
Out[39]=

Next, the PDE is set up. The load term of the PDE is the charge density divided by the vacuum permittivity :

The lower part of the geometry is silicon, while the upper part is vacuum. Silicon has a relative permittivity of about 11.7, while vacuum has a value of 1.

In[40]:=
Click for copyable input
Out[40]=

The charge density and vacuum permittivity are given as follows.

In[41]:=
Click for copyable input
Set up the partial differential operator.
In[43]:=
Click for copyable input
Out[43]=

For the use of Inactive in PDEs, see the section Formal Partial Differential Equations.

The potential is set to 0 wherever . A second potential is set to 1000 in the inner strip.

Set up the boundary conditions.
In[44]:=
Click for copyable input
Solve the PDE.
In[45]:=
Click for copyable input
Out[45]=
Visualize the result.
In[46]:=
Click for copyable input
Out[46]=

In many cases when there are explicit discontinuities in the PDE coefficients, NDSolve automatically includes these in the mesh it generates to solve the problem. The discontinuity at is automatically detected for this example.

Create a symbolic specification for the region.
In[47]:=
Click for copyable input
Out[49]=
Solve the PDE.
In[50]:=
Click for copyable input
Out[50]=
Inspect the boundary mesh that has been auto-created.
In[51]:=
Click for copyable input
Out[51]=

As can be seen, the original input region did not have an internal boundary at the discontinuity, but during the mesh generation process, the internal line along the discontinuity was added.

Partial Differential Equations with Nonlinear Coefficients

Some PDE coefficients may, in addition to space and time , also depend on the dependent variable and the first derivatives . Here indicates the spatial variables , , . If a coefficient depends on the dependent variable , the equation is nonlinear. Consider the nonlinear equation:

As an example, to compute a minimal surface the following equation can be used:

The equation is nonlinear, since the derivative of the dependent variable appears in the diffusion coefficient.

Set up a disk as a region.
In[270]:=
Click for copyable input
Set up the partial differential operator.
In[271]:=
Click for copyable input
Out[272]=

To solve nonlinear stationary equations, an initial seed for the solution can be specified. This is done by making use of InitialSeedings. InitialSeedings specifies an expression to be evaluated over the region of interest. The expression can depend on the spatial variables. If no initial seed is specified, 0 is assumed as an initial seed.

Solve the nonlinear equation.
In[273]:=
Click for copyable input
Out[273]=
Visualize the result.
In[275]:=
Click for copyable input
Out[275]=

Partial Differential Equations with Nonlinear Variable Coefficients

The following section presents a magneto statics application example. Magneto static applications typically make use of several, possibly nonlinear, materials in different areas of the simulation domain. This example looks at a simplified model of a motor, for which the distribution of the magnetic potential and a vector density field of the magnetic flux will be computed.

In the preceding representation of the motor geometry, the stator sections are gray, the rotor is red, and the blue area represents air. Also, several coils are grouped and displayed in orange, light orange and yellow.

The PDE model is derived from Ampère's circuital law of Maxwell equations:

Here is the current density and the permeability. The magnetic flux density is connected to the magnetic potential via

Inserting 1 in 2, you obtain

To simplify the model to two dimensions, assume that currents only flow orthogonal to the the and plane: . Further assuming that is constant in the direction and and are constant, you can write:

Since a model is assumed that is homogenous in the direction, the components drop out and the following is used:

The equation is nonlinear, since the permeability is a function of the spatial coordinates and the derivative of the unknown . In the rotor and stator, the permeability of a ferromagnetic material is used, while in the rest of the domain, the permeability of air will be used. The right-hand side of the equation, , is the -component of a variable current density that flows through the coils.

The first step is to create a mesh that has markers for the different regions. These markers can then be used to activate, for example, in different regions of the simulation domain. Using markers for this is much simpler than specifying the formula for regions where different parts of the PDE are active. More information on markers and their generation in meshes can be found in the "Element Mesh Generation" tutorial.

Use a manually generated mesh region.
In[64]:=
Click for copyable input
Generate a boundary element mesh from the MeshRegion.
In[65]:=
Click for copyable input
Out[65]=

For markers to be attributed to different subregions, coordinates lie in those subregions need to be specified.

Set up the marker coordinates.
In[66]:=
Click for copyable input

Visualizing the boundary mesh and the marker coordinates helps one to understand the process of specifying marker coordinates. For the visualization, it is necessary to group the marker coordinates and set up colors that will go with the markers.

Group the marker coordinates and set up marker colors.
In[73]:=
Click for copyable input
Visualize the boundary mesh and the marker coordinates in different colors.
In[75]:=
Click for copyable input
Out[75]=

The coils, currently visualized with an orange point, are to carry different currents. To simplify this process, each coil will get a marker assigned such that coils carrying the same current are grouped together.

Specify the coil marker.
In[76]:=
Click for copyable input
Check that there are as many coil markers as there are coil coordinates.
In[77]:=
Click for copyable input
Out[77]=

Now you have a list of coordinates lying in the subregions and have attributed markers with those coordinates.

Specify the markers.
In[78]:=
Click for copyable input
Out[78]//Short=

Next, create the full mesh.

Generate the mesh with markers.
In[79]:=
Click for copyable input
Out[79]=
Visualize the mesh with the stator in gray, the rotor in red, air in blue and the coils in orange, light orange and yellow.
In[80]:=
Click for copyable input
Out[80]=

Now the material data is set up.

Set up the permeability of air.
In[81]:=
Click for copyable input

The change of permeability for ferromagnetic materials is typically given in BH-curves. These contain material data. For this data, an interpolating function can be created that can be used as a PDE coefficient. As a more efficient alternative, a fit can be generated that can be used as the PDE coefficient. Both approaches will be shown.

Specify values of a ferromagnetic material from the BH-curve data sheet.
In[82]:=
Click for copyable input
Set up the BH-curve data.
In[84]:=
Click for copyable input

While generating the interpolation, make sure that the smallest and largest values of the BH-curve are continued if the interpolating function is queried outside of its domain. Also instruct the interpolating function to not give warning messages when queried outside of the domain.

Generate an interpolating function from the BH-curve data.
In[85]:=
Click for copyable input
Compute the norm of the magnetic flux density and avoid it becoming zero.
In[86]:=
Click for copyable input

If the permeability coefficient is evaluated either in the rotor or the stator, then the ferromagnetic material will be used in all other regions the permeability for air is used.

Specify an inverse permeability based on the interpolated data.
In[87]:=
Click for copyable input

With the interpolating function approach comes great flexibility, and the coefficient can be used as it is. However, a more efficiently evaluatable approach is to find a fit for the BH-curve data.

Find a fit for the BH-curve data with a Gaussian model and generate a function.
In[88]:=
Click for copyable input
Out[91]=
Visualize the difference between the interpolated and the fitted BH-curve.
In[92]:=
Click for copyable input
Out[92]=
Specify an inverse permeability that is efficiently evaluatable.
In[93]:=
Click for copyable input

The current density in the coils is specified such that there is a positive current in elements with ElementMarker 4, a negative current in the elements with ElementMarker 6 and a zero current elsewhere.

Specify the current density.
In[94]:=
Click for copyable input
Specify the partial differential equation.
In[96]:=
Click for copyable input

The boundary condition is set up in such a way that the magnetic potential is zero at the out rim of the stator.

Set up the boundary condition.
In[97]:=
Click for copyable input
Solve the nonlinear PDE.
In[98]:=
Click for copyable input
Out[98]=
Find the minimum and maximum of the magnetic potential.
In[99]:=
Click for copyable input
Out[99]=
Visualize the magnetic potential and magnetic flux density B.
In[100]:=
Click for copyable input
Out[100]=

Partial Differential Equations and Nonlinear Boundary Conditions

NDSolve and related functions allow for specifying nonlinear generalized Neumann values. Dirichlet and periodic boundary conditions cannot be nonlinear. A nonlinear generalized Neumann value is given by:

This is very similar to the linear generalized Neumann value, except that the dependent variable is now a function of .

A common use case for generalized nonlinear Neumann value is to model radiating boundaries. As a one-dimensional example, a nonlinear radiating heat transfer is used:

The boundary condition on the right is a nonlinear radiation condition:

On the left boundary, a Dirichlet condition is made use of. Here is the thermal conductivity coefficient, is a cross-sectional area exposed to radiation, is the StefanBoltzmann constant, is the emissivity, is the ambient temperature, and is a heat source.

Set up the nonlinear equation.
In[2]:=
Click for copyable input
Set up the radiation boundary condition.
In[3]:=
Click for copyable input
Set up a Dirichlet condition.
In[4]:=
Click for copyable input
Set up material parameters.
In[5]:=
Click for copyable input
Set up the PDE with material parameters substituted.
In[6]:=
Click for copyable input
Out[6]=
Solve the equation.
In[8]:=
Click for copyable input
Out[8]=
Visualize the result.
In[9]:=
Click for copyable input
Out[9]=

For a verification of a similar example, please refer to the nonlinear Finite Element verification tests.

Heat Equation

The heat equation is a time-dependent Poisson equation

where the dependent variable depends on the spatial coordinates and time .

Both Dirichlet boundary conditions and Neumann boundary values may also depend on time.

The overall procedure to solve PDEs remains the same: a region needs to be specified and a PDE with boundary conditions needs to be set up.

Set up the spatial region.
In[65]:=
Click for copyable input

The dependent variable now depends on time and space.

Set up the PDE operator, being sure to include the time as an argument of .
In[66]:=
Click for copyable input

Also, the boundary conditions may depend on space and time.

Set up the boundary conditions.
In[67]:=
Click for copyable input
In[68]:=
Click for copyable input

To solve first-order time-dependent problems, NDSolve needs an initial condition and a time region.

Solve the PDE.
In[69]:=
Click for copyable input
Make an animation of the solution.
In[72]:=
Click for copyable input
Out[73]=

Note that it is more efficient to convert the mesh once and use the discretized mesh during the plotting instead of rediscretizing the region for every animation frame.

Wave Equation

The wave equation is a second-order time-dependent equation:

Define circular region .
In[59]:=
Click for copyable input
Define the wave equation PDE operator.
In[60]:=
Click for copyable input
Dirichlet boundary conditions are set to a fixed value of on all boundaries.
In[61]:=
Click for copyable input

The wave equation is a second-order time-dependent PDE, and initial conditions and the derivative of the initial condition must be specified.

Initial conditions for the wave equation.
In[62]:=
Click for copyable input
In[63]:=
Click for copyable input
Solve the PDE and have NDSolve find initial values.
In[64]:=
Click for copyable input
Out[64]=

To make an animation of the resulting InterpolatingFunction, use Plot3D and Animate.

Make an animation of the solution.
In[65]:=
Click for copyable input
Out[66]=

If no boundary conditions are specified, which implicitly means a Neumann zero boundary condition, then the wave is reflected at the boundary. This can be seen well in a 1D example.

Define line region .
In[67]:=
Click for copyable input
Define a wave operator.
In[68]:=
Click for copyable input
Initial conditions for the wave equation.
In[69]:=
Click for copyable input
Solve the wave equation with reflecting boundaries.
In[71]:=
Click for copyable input
Out[71]=

To make an animation of the resulting InterpolatingFunction, use Plot3D and Animate.

In[72]:=
Click for copyable input
Out[73]=

For time-dependent systems, boundary conditions may be specified for any time derivative of a dependent variable of order strictly less than the time derivative order of that dependent variable in the equation. Thus for equations that are first order in time, boundary conditions may only be given for the dependent variables. For equations that are second order in time, boundary conditions may be given for the dependent variables and their first derivative with respect to time.

To model a wave equation with absorbing boundary conditions, one can proceed by using a temporal derivative of a Neumann boundary condition. For this, the exact same wave equation as in the preceding is used, but a Neumann value is added to the equation.

Solve the wave equation with absorbing boundaries.
In[74]:=
Click for copyable input
Out[74]=

To make an animation of the resulting InterpolatingFunction, use Plot3D and Animate.

In[75]:=
Click for copyable input
Out[76]=

More information about the wave equation and absorbing boundary conditions can be found in the Acoustics in the time domain tutorial.

Modeling periodic boundary conditions of waves can be implemented with PeriodicBoundaryCondition. In this example, the derivative of the initial condition is not zero but chosen such that the initial wave travels to the right.

Initial conditions for the wave equation.
In[77]:=
Click for copyable input
Specify a periodic boundary condition that absorbs the wave on the right and reinserts it on the left.
In[78]:=
Click for copyable input
Solve the wave equation with absorbing boundaries.
In[79]:=
Click for copyable input
Out[79]=

To make an animation of the resulting InterpolatingFunction, use Plot3D and Animate.

In[80]:=
Click for copyable input
Out[81]=

Formal Partial Differential Equations

A PDE that is making use of Laplacian can also be formulated with Div and Grad. To do so, use the coefficient in . Here is an × matrix where is the number of dimensions in .

In[82]:=
Click for copyable input
Out[82]=

To prevent evaluation, Inactive can be used.

Set up the PDE operator.
In[83]:=
Click for copyable input
Out[84]=

The formal version of the PDE allows you to easily modify the coefficient matrix .

One scenario where the use of Inactive can be important is when modeling the relation between a PDE and a Neumann value. Consider the following PDE:

The corresponding Neumann boundary equation is

Consider a PDE where is the identity matrix and is given by . Specifying the PDE without the use of Inactive will cause the expression to evaluate.

In[85]:=
Click for copyable input
Out[87]=

Note, however, that the evaluation changes the PDE to the following:

is and is .

In[88]:=
Click for copyable input

These PDEs are equivalent.

In[90]:=
Click for copyable input
Out[90]=

Even though the PDEs are equivalent, the corresponding Neumann values are not:

For the second PDE, the Neumann value modeled is the following:

To prevent evaluation and model a Neumann value of the following form, the use of Inactivate is needed.

In[91]:=
Click for copyable input

A more in-depth discussion of formal PDEs can be found in "Finite Element Method Usage Tips".

Also, note that while it is possible to specify a coefficient as a function, this deprives the internal code of the ability to optimize the coefficient integration. In this specific case, the coefficient is constant, and with the optimization, the coefficient is only integrated once. A black box function would need continuous evaluations over the region. For example, a function that hides its internals is not ideal, since the FEM code cannot make use of the optimized diffusion operator for constant coefficients.

A black box function hiding its internals that should be avoided as a PDE coefficient.
Click for copyable input

Systems of Partial Differential Equations

The Coefficient Form of Systems of Partial Differential Equations

Consider a system of two equations with dependent variables and :

For a region in , each is an × matrix, and each of , and are vectors.

In this case, the generalized Neumann equations are:

One-Dimensional Coupled Partial Differential Equation

The following example is set up such that the parameters are easy to modify. First, the region is specified.

Specify a region.
In[92]:=
Click for copyable input
Out[93]=

Here is the PDE to be solved:

The PDE has two dependent variables and . The coupling is between two diffusion operators. Note that since this is a one-dimensional system, each coefficient is a 1×1 matrix, which is equivalent to a number.

Set up the PDE operators.
In[94]:=
Click for copyable input

Each equation has its own set of boundary conditions. In the most general form, Dirichlet boundary conditions are now systems of conditions:

Neumann values specify a system of the following form:

Next, the boundary conditions are set up.

Specify boundary conditions for the first equation.
In[97]:=
Click for copyable input
In[99]:=
Click for copyable input
Specify boundary conditions for the second equation.
In[101]:=
Click for copyable input
In[103]:=
Click for copyable input
Solve the system of equations.
In[105]:=
Click for copyable input

Here is the analytical solution of the system of PDEs computed with DSolve.

The analytical solution.
In[106]:=
Click for copyable input
Out[106]=

Note how in the DSolve formulation of the system of equations the generalized Neumann boundary equations have to be given explicitly, and how their coefficients match those in the PDE.

To compare the analytical result to the numeric result, make a plot of the difference between the two.

Plot the difference between the analytical and numerical solutions.
In[107]:=
Click for copyable input
Out[107]=

This example also makes apparent why the NeumannValue is part of the equation. The equations for the Neumann value are given by:

Consider the boundary value .

In[108]:=
Click for copyable input
Out[108]=

There is no way to know from the Neumann value with which equation it should be associated. Making the Neumann value part of the equation disambiguates with which equation the Neumann value is associated.

Structural Mechanics

In structural mechanics, the deformation of objects under load is computed. Structural mechanics provides good examples of coupled PDEs. As an example region, consider a beam five units long and one unit high.

The PDE that models the physics is a plane stress relation and is valid for thin objects. Here, "thin" means thin relative to the other dimensions of the object. The plane stress formulation depends on two dependent variables and the independent variables. The two dependent variables give the component-wise contribution of the deformed object. As material data, the Young's modulus and the Poisson ratio need to be specified.

Define a plane stress operator.
In[109]:=
Click for copyable input

Furthermore, a boundary load may be specified. In this case, the boundary load applies at the right edge of the beam and a force of unit is applied in the negative direction.

On the left-hand side of the region, both dependent variables are held fixed with 0 displacement.

Set up the boundary conditions so that the beam is held fixed at the left.
In[111]:=
Click for copyable input
Solve the equation.
In[112]:=
Click for copyable input

In order to show the deformed beam, an element mesh must be created for the beam at rest, which is deformed by the interpolating functions given.

Visualize the boundary of the beam at rest versus the deformed mesh of the beam.
In[113]:=
Click for copyable input
Out[114]=

More information about the usage of ElementMeshDeformation can be found in "Element Mesh Visualization".

The following contour plots show the deformation of the beam. The dependent represents the deformation in the direction. It can be seen that the top-right part of the beam experiences a shift in the positive direction. The lower-right part is moved in the negative direction.

In[115]:=
Click for copyable input
Out[115]=

The contour plot of the dependent represents the deformation in the direction. The right-hand side of the beam undergoes a deformation in the negative direction.

In[116]:=
Click for copyable input
Out[116]=

Note that the color scale of the two contour plots is not the same.

For completeness, a PDE that models the physics of a plane strain relation is presented. The plane strain formulation is, in contrast to the plane stress formulation, valid for thick objects. Here, "thick" means thick relative to the other dimensions of the object. The plane strain formulation depends on two dependent variables and the independent variables. The two dependent variables give the component-wise contribution of the deformed object. As material data, the Young's modulus and the Poisson ratio need to be specified.

Define a plane strain operator.
In[117]:=
Click for copyable input

The boundary conditions this time are prescribed displacements on both sides. On the left-hand side of the region, both dependent variables are held fixed with 0 displacement. On the right-hand side, a displacement of one unit in the direction of negative axes is given and the displacement on the direction is held fixed.

Set up the boundary conditions so that the beam is held fixed at the left and the right.
In[119]:=
Click for copyable input
Solve the equation.
In[120]:=
Click for copyable input

In order to show the deformed beam, an element mesh must be created for the beam at rest, which is deformed by the interpolating functions given.

Visualize the boundary of the beam at rest versus the deformed mesh of the beam.
In[121]:=
Click for copyable input
Out[122]=

The following contour plots show the additive deformation of the beam. The dependent represents the deformation in the direction, and the dependent represents the deformation in the direction.

In[123]:=
Click for copyable input
Out[123]=

Fluid flow

A linear model of a fluid flow is Stokes's equation. In this example, a flow through a narrowing tube is analyzed.

Specify and visualize a region that narrows.
In[124]:=
Click for copyable input
Out[125]=
The Stokes operator is set for a viscosity .
In[126]:=
Click for copyable input

Here represents the fluid's velocity in the direction and represents the velocity in the direction. stands for the pressure in the fluid.

Specify the PDE.
In[127]:=
Click for copyable input

The inflow velocity in the direction is prescribed with the boundary conditions. This is done by specifying an inflow profile for the dependent variable on the left-hand side of the region. On the bottom and upper part of the region, a no-flow condition is specified. This is done by setting both the and components on those parts to 0. An outflow condition is set on the right-hand side of the domain. Here the pressure is arbitrarily set to 0.

Setup of the boundary conditions.
In[128]:=
Click for copyable input

A stable solution can be found if the velocities are interpolated with a higher order than the pressure. NDSolve allows an interpolation order for each dependent variable to be specified.

Solve the equations while specifying a finer mesh than the default.
In[129]:=
Click for copyable input

An explanation of the usage of the finite element method option interpolation order is given in "Finite Element Method Usage Tips".

Visualize the flow in the direction.
In[130]:=
Click for copyable input
Out[130]=
Visualize the flow in the direction.
In[131]:=
Click for copyable input
Out[131]=
Visualize the pressure.
In[132]:=
Click for copyable input
Out[132]=
Visualize the flow field.
In[133]:=
Click for copyable input
Out[133]=

A nonlinear model of a fluid flow is NavierStokes's equation. In this example, a flow-through around a cylinder is analyzed by looking at a 2D cross section.

Specify and visualize a region around a cross section of a cylinder.
In[2]:=
Click for copyable input
Out[4]=
The NavierStokes operator is set for a viscosity and density .
In[248]:=
Click for copyable input

Here represents the fluid's velocity in the direction, and represents the velocity in the direction. stands for the pressure in the fluid.

Specify the PDE.
In[249]:=
Click for copyable input

The inflow velocity in the direction is prescribed with the boundary conditions. This is done by specifying an inflow profile for the -dependent variable on the left-hand side of the region. On the bottom and upper part of the region and around the cylinder, a no-flow condition is specified. This is done by setting both the and components on those parts to 0. An outflow condition is set on the right-hand side of the domain. Here the pressure is arbitrarily set to 0.

Setup of the boundary conditions.
In[250]:=
Click for copyable input

Suppose you would like to compute the pressure difference from a point before the cylinder and a point behind the cylinder . Additionally, the length of recirculation is to be computed. Also the drag and lift force on the cylinder are of interest. is the coordinate of the end of the cylinder and is the end of the recirculation area. In order to get a good result, the area around the cylinder is to be refined in a piriform shape, and a higher than default accuracy and precision are to be used for the mesh generation.

Set up a piriform refinement region around the cylinder.
In[251]:=
Click for copyable input
Visualize the refinement region and the simulation domain .
In[252]:=
Click for copyable input
Out[252]=
Create a refinement function based on the region of refinement.
In[253]:=
Click for copyable input

A stable solution can be found if the velocities are interpolated with a higher order than the pressure. NDSolve allows an interpolation order for each dependent variable to be specified.

Solve the equations, while specifying points to be included, increased accuracy and precision and a refinement function.
In[254]:=
Click for copyable input

An explanation of the usage of the finite element method option interpolation order is given in "Finite Element Method Usage Tips".

Compute the pressure difference before and after the cylinder.
In[255]:=
Click for copyable input
Out[255]=

The expected value for the pressure difference is between 0.1172 and 0.1176.

Find the recirculation length after the cylinder.
In[256]:=
Click for copyable input
Out[256]=

The expected value for the recirculation length is between 0.0842 and 0.0852.

Set up a function to compute the drag and lift force around the cylinder.
In[275]:=
Click for copyable input
Compute the drag and lift force on the cylinder.
In[276]:=
Click for copyable input
Out[276]=

The expected value for the drag force is between 5.57 and 5.59, and for the lift force between 0.0104 and 0.0110.

Visualize the and components of the velocity as , pressure contours
as black lines, and the velocity vector field.
In[246]:=
Click for copyable input
Out[246]=

The transient NavierStokes equation can also be solved. For this, the same geometry as in the stationary case is used. First the time-dependent NavierStokes equation is set up.

The transient NavierStokes operator is set for a viscosity and density .
In[247]:=
Click for copyable input

Suppose you would like to model an inflow on the left, set a pressure condition on the outflow on the right, and all remaining walls have a non-slip boundary condition. For the inflow boundary condition to be consistent, a helper function that smoothly ramps up the inflow over time is created.

Function to ramp up the inflow boundary condition over time.
In[56]:=
Click for copyable input
Out[58]=
Set up the boundary conditions.
In[59]:=
Click for copyable input
The initial conditions are set to 0.
In[23]:=
Click for copyable input

Because the Navier-Stokes equation is a high-index differential algebraic equation, options are specified to time-integrate the equations efficiently. The maximum difference order of the time integrator is reduced to minimize the number of rejected steps. The time-dependent boundary conditions are differentiated to create a consistent system of equations that are more efficient to time integrate. Lastly, a mesh somewhat finer than the default mesh is used. Solving this equation can, depending on the hardware used, take a few minutes.

Time-integrate the NavierStokes equation.
In[60]:=
Click for copyable input
Out[60]=
Out[61]=
Find the maximum and minimum of the velocity range.
In[62]:=
Click for copyable input
Out[62]=
For efficient visualization, extract the mesh from the solution.
In[63]:=
Click for copyable input
Out[63]=
Create contour plots of the velocity.
In[37]:=
Click for copyable input
Out[37]=
Visualize the velocity component.
In[64]:=
Click for copyable input
Out[64]=

As a next step, the NavierStokes equation is coupled to a heat equation to model energy transport.

First a region and some parameters are specified.

Specify parameters and a rectangular region with cutouts.
In[2]:=
Click for copyable input
Visualize the region.
In[4]:=
Click for copyable input
Out[4]=
Specify a Prandtl number and a Rayleigh number .
In[5]:=
Click for copyable input

Set up a viscous NavierStokes equation that is coupled to a heat equation making use of a Boussinesq approximation. Use the material parameters specified.

Set up a NavierStokes equation coupled to a heat equation.
In[6]:=
Click for copyable input

Boundary conditions and initial conditions need to be set up.

Set up no-slip boundary conditions for the velocities on all boundary walls.
In[7]:=
Click for copyable input
Set up a reference pressure point.
In[8]:=
Click for copyable input
Specify a temperature difference between the left and right walls.
In[9]:=
Click for copyable input
Replace parameters in the boundary conditions.
In[10]:=
Click for copyable input
Set up initial conditions such that the system is at rest.
In[11]:=
Click for copyable input

Monitor the time integration progress and the total time it takes to solve the PDE while using a refined mesh and interpolating the velocities and and the temperature with second order and the pressure with first order.

Monitor the time integration of the equations.
In[12]:=
Click for copyable input
Out[12]=

The last step is post-processing.

Construct a visualization of the region's boundary.
In[13]:=
Click for copyable input
Out[14]=
Visualize the pressure distribution.
In[15]:=
Click for copyable input
Out[15]=
Visualize the temperature distribution.
In[16]:=
Click for copyable input
Out[16]=
Visualize the velocity field.
In[17]:=
Click for copyable input
Out[17]=
Animate the change in temperature and the velocity streamlines.
In[18]:=
Click for copyable input
Out[18]=

For the curious, investigate the different velocity field produced by simply swapping the cut-out positions.

A different region to try.
In[24]:=
Click for copyable input
Out[24]=