学习中心

Modeling with PDEs: Poisson's and Laplace Equations

---

In Part 1 of this course on modeling with partial differential equations (PDEs), we begin with a quick introduction to using the general-purpose PDE interfaces in the COMSOL Multiphysics^®software. You will learn how to use the interfaces, which are based on the finite element method (FEM) and boundary element method (BEM), for modeling with Poisson's and Laplace equations, respectively. Built-in mathematics interfaces are available for modeling these specific PDEs (i.e., theLaplace Equationinterface andPoisson's Equationinterface), however, in order to demonstrate how to implement your own custom equations we will not use them. You will also see an example of how to solve for the Newtonian gravity of the Earth–Moon system, since this does not have any corresponding built-in options in the software. However, the techniques used here are applicable in a much wider context. In subsequent parts, we will learn more about the mathematics interfaces and how to model more general PDEs as well as systems of PDEs.

The Mathematics Interfaces in COMSOL Multiphysics^®

The mathematics interfaces can be found in the Model Wizard underMathematics, as seen in the figure below.

TheMathematicsinterfaces in the Model Wizard.

With these interfaces, you can solve various types of PDEs using different formulations and numerical methods. You can also use these interfaces to solve ordinary differential equations (ODEs) and other special types of equations. You can combine equations of different types to form a wide range of nonlinear equation systems, modeling all kinds of mathematical and physical processes. All of these interfaces, with the exception of some of those under theOptimization and SensitivityandMoving Interfacebranches, are available as part of the core functionality in the software.

The Finite Element and Boundary Element Methods

In COMSOL Multiphysics^®, the main method for solving PDEs is the finite element method and most of the mathematics and physics interfaces in the software are based on this method. However, in some cases it is beneficial to instead use BEM, frequently in combination with FEM. BEM uses a mesh on surfaces only, i.e., not in the computational volume. This is advantageous when the computational volume is large, which typically happens when the objects of interest are far apart. As a fun example of this, let us see how we can solve Laplace and Poisson's equation for the gravitational field of the Earth–Moon system (there is no corresponding built-in interface for this purpose). The Earth and the Moon are relatively far apart compared to their radii, which makes modeling the space that surrounds them suitable for BEM. To make things simple, we will assume classical stationary Newtonian gravity.

For electrostatics, corrosion, magnetostatics, and acoustics, there are dedicated physics interfaces that combine finite element and boundary element formulations. Since this is not the case for gravity, we will instead create custom equations and equation couplings by combining the finite-element-basedCoefficient Form PDEinterface with the boundary-element-basedPDE, Boundary Elementsinterface. An additional benefit with thePDE, Boundary Elementsinterface is that it automatically covers an infinite computational volume. Note that in our case we will ignore any other celestial objects, but that the model presented here can easily be extended by adding the Sun or other objects.

Newtonian Gravity

Let us first introduce some notation. We will model gravity using partial differential equations. For this purpose, we need the gravitational vector field,, having units of acceleration in SI units,. This field has three components in thex,y, andzdirection:

To make computations easier and be able to formulate a scalar-valued PDE, we will also need the gravitational potential, which is a scalar field,, with the somewhat unfamiliar unit of. We will see later how you can use theMathematicsinterfaces to specify custom units and integrate them with the automatic unit handling capabilities in the software.

We will use the following notation from vector calculus:

for the gradient of the scalar field.
for the divergence of the vector field.

Gauss's law for gravitationis:,
whereis the gravitational constant andis the mass density. We can include the varying density due to the layering of, for example, the Earth's interior if we want to. However, in our example we will model the Earth and Moon with their average mass density values: 5515and 3340, respectively. In the software, we will define these density values asVariables. This means you can change them to be functions of the coordinate variablesx,y, andz, or, if you prefer spherical coordinates,r,theta, andphi.

The relationship between the gravitational field and the gravitational potential is:

By combining this we get the classical scalar equation for gravity, which is the following version of Poisson's equation:

The gravitational constant has the unitand the density has the unit. Since these are multiplied together, the right-hand side of this equation has the unitor, equivalently,. Here, we will use the unit, which can be shown to also be equivalent to.

Since there is no mass density in free space,and we can use the Laplace equation:

Note that the gravity vector field,, and the gravity potential,, have the same type of relationship as the electric field,, and electric potential,, for electrostatics:

Also, there are dedicated physics interfaces for this that are available from the Model Wizard or theAdd Physicswindow in theAC/DCbranch, under theElectric Fields and Currentsbranch. In principle, we could use theElectrostaticsinterfaces for modeling gravity, but this would come with some drawbacks, such as unsuitable names for boundary conditions and the wrong units.

The Coefficient Form PDE Interface

One of the most important tools in COMSOL Multiphysics^®for equation-based modeling is theCoefficient Form PDEinterface. This equation template provides a powerful general interface for specifying linear and nonlinear equations, including the classical Poisson's and Laplace equations. Note that there areMathematicsinterfaces specifically for Poisson's and Laplace equations. These interfaces are available from the Model Wizard underMathematics>Classical PDEsand could have been used here. However, the goal of this article is to learn how to use the general-purpose PDE interfaces, such as theCoefficient Form PDE interface.

Using the notation of the user interface in the case where we have a scalar field, the equation for theCoefficient Form PDEinterface reads:

There are two important boundary conditions when using this interface. First, is the generalizedNeumannboundary condition:

,

withbeing the outward surface normal, which is used to model a boundary flux. The second is theDirichletboundary condition:

which is used to assign fixed values of the dependent variable on the boundary.

If we assume that we have just one dependent variable and that all materials are isotropic, then most of the coefficients are scalar quantities, but some are vector-valued as follows:

Scalars:
Vectors:

We will learn more about these coefficients in later articles. Note that for more general cases, the field quantities and coefficients can be vectors or matrices (tensors).

In the example of Newtonian gravity, we will now see how we can use theCoefficient Form PDEinterface to define a PDE by simply matching the coefficients.

Poisson's Equation

In this example, we will only need the stationary version of the Coefficient Form PDE. In other words, we assume:

for all times of interest. This also means that we ignore the influence of the coefficientsand.

We will not need many of the other coefficients. In our case, we have:

and end up with the equation (moving the negative sign in front):

To model our particular gravitational version of Poisson's equation, we will only use the coefficientsandand leave all other coefficients with their default values (zero). Specifically for modeling the equation

we will have:

the dependent variable:
the diffusion coefficient:
the source term:

We will shortly learn how to set this up from the beginning using the Model Wizard. To give you an idea on how this PDE will appear in the user interface, the correspondingCoefficient Form PDEinterface is shown in the figure below.

The gravity equation modeled using theCoefficient Form PDEinterface.

whereGandrhoare defined as a parameter and a variable for the gravitational constant and density elsewhere in the model, respectively. The value ofrhowill be that of the Earth, 5515; and the Moon, 3340, within the two spheres that we will use to represent these bodies.

The Model Wizard

Start by setting up a new model by, for example, selectingFile>New. In the Model Wizard, select3Das the space dimension, as shown in the figure below.

Selecting3Dfor the space dimension in the Model Wizard.

For the next step, in the Model Wizard, selectMathematics>PDE Interfaces>Coefficient Form PDE. Change both theField nameandDependent variables nametoV, as shown below. For more general models, theField namecan be used to refer to an array or vector of dependent variables.

The Select Physics page with the Added physics interfaces field on the left and the Review Physics Interface window on the right showing the dependent variables and units.
TheSelect Physics page in the Model Wizard.

As shown in the image above, we need to change to custom units for the dependent variable quantity, which is the gravity potential with unit, as well as the source term with unit.

Click theDefine Dependent Variable Unitbutton, available to the right of theDependent variable quantitytable in theUnitssection. TypeJ/kg, for theUnit, as shown here:

TheCustom unitsettings for theDependent variable quantity.

Click theDefine Source Term Unitbutton, available to the right of theSource term quantitytable in theUnitssection. TypeJ/(m^2*kg), for theUnit, as shown in the figure below.

TheCustom unitsettings for theSource term quantity.

For the next step, in the Model Wizard, selectStationaryfor the study type.

TheSelect Studypage of the Model Wizard.

ClickDoneto enter the Model Builder 3D workspace.

Using the Coefficient Form PDE for Poisson's Equation

We need to define a few parameters for the radii of the Earth,r_earth, and the Moon,r_moon, as well as the distance to the moon,d_moon. There is a built-in parameter for the gravitational constant,G_const. For convenience, we redefined this as a parameterG. This is all shown in the figure below.

The parameters used in the example model.

The geometry simply consists of two spheres at the distanced_moonapart, as shown in the figure below.

The COMSOL Multiphysics UI showing the Model Builder and Sphere settings window, with the geometry for the Earth–Moon system in the Graphics window.
The geometry model for the Earth–Moon system.

In the model, we have changed the length unit from m to km. Note that we do not need to enclose the two spheres in a volume mesh due to the fact that we will use a boundary element method for the gravity in space. This model could alternatively be defined as a 2D axisymmetric model.

The Earth–Moon system modeled as two spheres with a length unit set to km.

The average density of the Earth and the Moon are defined as two domain variables with the same name,rho, taking different values in domain 1 (Earth) and domain 2 (Moon), as shown in the figures below.

The variable for the average density of the Earth.

The variable for the average density of the Moon.

Note that you can make the expression for the varying density in space, as described earlier.

In the next step, define the coefficients for theCoefficient Form PDEsettings as described earlier (also shown in the figure below).

The COMSOL Multiphysics UI showing the Model Builder and the Coefficient Form PDE settings window, with the Earth–Moon system model in the Graphics window.
TheCoefficient Form PDE settings for the Earth–Moon system.

This concludes the setup of the finite-element-based Coefficient Form PDE. In the next step we add aPDE, Boundary Elementsinterface.

Using the PDE, Boundary Elements Interface for Laplace Equation

Since we will model the Earth–Moon system as an idealized and isolated system, we do not want to impose unnecessary and artificial boundary conditions. We would like to view the system as situated in infinite free space without external disturbances. For this purpose, we will use thePDE, Boundary Elementsinterface. This interface will more or less automatically handle the infinite computational domain outside (and between) the Earth–Moon system.

ThePDE, Boundary Elementsinterface covers a stationary PDE of the following form:

Notice that this form does not feature a source term, so we cannot use it to solve Poisson's equation. However, we can use it in free space since that is where the gravitational potential is governed by the Laplace equation:

In identifying the coefficients we see that:

the dependent variable:

the diffusion coefficient:

the absorption coefficient, or Helmholtz coefficient:

ThePDE, Boundary Elementsinterface is available in the list ofMathematicsinterfaces. To add this interface, you can right-clickComponent 1and select from theAdd Physicswindow. The interface is available at:Mathematics>PDE Interfaces>PDE, Boundary Elements, as shown below.

The COMSOL Multiphysics UI showing the Model Builder, Component settings, the Earth–Moon system in the Graphics window, and the Add Physics window with PDE, Boundary Elements highlighted.
The PDE, Boundary Elementsinterface is added to the model through the Add Physicswindow.

By default, domain 1 and 2 and theInfinite voiddomain are selected for thePDE, Boundary Elements. TheInfinite voidis COMSOL terminology for the space surrounding the spheres and all the way out to infinity. We need to remove domain 1 and 2 since we only want to apply thePDE, Boundary Elementsin theInfinite voiddomain, as shown below.

The domain selection for thePDE, Boundary Elementsinterface before deleting domains 1 and 2 from the geometry selection.

The domain selection for thePDE, Boundary Elementsinterface after deleting domains 1 and 2 from the geometry selection.

The next step is to change the units of thePDE, Boundary Elementsinterface to the same as that of theCoefficient Form PDEinterface, as shown in the figure below.

TheUnitssettings for thePDE, Boundary Elementsinterface.

After this step we can change the name of the dependent variable toV, as shown in the figure below.

TheDependent variablesettings for thePDE, Boundary Elementsinterface.

It is a feature in COMSOL Multiphysics^®that by using the same name of the dependent variable in both interfaces, we get a consistent coupling between the two PDEs on any common boundary. For this to be possible, the units need to be the same in the two interfaces.

Recall that our strategy is to define Poisson's equation using theCoefficient Form PDEinterface inside of the Earth and Moon, modeled as spheres, and then define aPDE, Boundary Elementsinterface to represent the gravity field in the surrounding space. In the figure below, you can see what this now looks like in the model tree.

The mathematics interfaces used in the Earth–Moon gravity model.

We now need to specify an asymptotic value for the gravity potential, or theCondition at infinity. We will set this value to zero in theSettingswindow of thePDE, Boundary Elementsinterface. You can see these settings in the figure below.

TheInfinite voidandCondition at infinitysettings for thePDE, Boundary Elementsinterface.

If we use the analogy between Newtonian gravitation and electrostatics, we can view this asymptotic condition as a "Ground" boundary condition at a sphere very far away from the Earth–Moon system. In this way, we do not need to specify any actual boundary conditions for theCoefficient Form PDEinterface. Instead, thisCondition at infinityis propagating — from theInfinite voidwhere thePDE, Boundary Elementsinterface is defined — inward via the coupling of the boundaries between the void and the two domains for the spheres. Again, the coupling is simply defined by assigning the same name of the dependent variable to the two equation interfaces.

With respect to the finite element mesh used in the spheres, a fixed element size of 500 km is used for the Earth, as shown in the figure below, and 200 km is used for the Moon. This is accomplished by setting the same value for the maximum and minimum element size. Inside of the spheres there will be tetrahedral finite elements. The boundary element mesh will consist of the triangular mesh at the surface of the spheres, adjacent to theInfinite voiddomain.

The COMSOL Multiphysics UI showing the Model Builder; the Mesh 1, Size 1 settings; and the Earth and Moon models in the Graphics window.
The settings for the finite element mesh of the Earth model.

Now, we can selectComputeto solve the PDE system. The computation time is a few minutes.

The Results

We can plot the magnitude of the surface gravity,

corresponding to the expressionsqrt(Vx^2+Vy^2+Vz^2), on the two planetary bodies, as shown in the figure below. You can type a mathematical expression like this one in the text fields for the expression to be visualized.

The expression used in a surface plot to visualize the magnitude of the surface gravity.

ASurfaceplot displaying the magnitude of surface gravity on the two planetary bodies.

Using a mouse, we can click on the surface of each sphere and read off the values 9.81for the Earth and 1.62for the Moon, which matches the expected values within less than a percent.

The table displaying the surface gravity, generated after clicking on the surface of each sphere in the results plot.

We can plot the magnitude of the gravity along a centerline in the space between the Earth and the Moon by plotting the variablepdebe.normgrad_u, a built-in variable used to efficiently compute the magnitude of the gravity vector field coming from the boundary element method solution. This gives us the plot shown below.

The settings for the plot where the built-in variable is used in the expression field.

The gravity on a line between the Earth and the Moon using logarithmic scaling.

We can see a minimum at about 345,000 km. Note, however, that this does not correspond to one of the so-calledLagrange points. To compute the location of the Lagrange points, we would need to include the counteracting acceleration effects of a rotating frame due to orbital motion. Note that theGriddata settings were increased to produce this plot, as shown in the figure below. TheGriddataset adds volumetric grid cells so that boundary element quantities can be visualized and evaluated in void domains.

The settings for theGrid 3Ddataset.

We can also visualize the gravitational field using colors, arrows, and contour lines, as shown in the figure below.

A slice plot of the Earth–Moon system compared to a blue color scale and Rainbow color scale with green arrows indicating the gravity direction.
In this visualization of the Earth–Moon system, the slice plot shows the gravity in space for values below 0.05 m/s². The surface plots shows the gravity on the surfaces, the contours show the gravity strength, and the arrows show the direction of the gravitational field.

From this plot we can, for example, see that the Earth's gravity dominates a large area also close to the Moon. Only when we come fairly close to the Moon does the Moon's gravity take over.

As a final plot, we can create a "gravity well" plot by using a deformation modification of a cut plane plot, as shown in the figure below.

A gravity well plot of the Earth–Moon system.

The Inverse-Square Law

Since the PDE for gravitation might not be familiar, let us investigate how this relates to the more familiar inverse-square law for the gravitational force. The inverse-square law is actually encoded in Poisson's equation. To see this, consider a point sourcehaving mass, located at the origin and the corresponding "point density":

For the mathematically inclined reader, this equation is meaningful only in the weak sense.

The analytical solution to this equation is:

whereis the distance from the point source (origin).

We get the inverse-square law by taking the gradient in spherical coordinates:

whereis a unit vector pointing outward in the radial direction.

For more information on spherical coordinates, seethis Wikipedia page.

We get the forceon a particle in this gravity field, having mass, by the relationship

which gives us the inverse square law:

Built-in Interfaces Based on the Boundary Element Method

For more information on using the built-in physics interfaces that are based on the boundary element method, please see the following blog posts:

• Acoustics:"How to Use the Boundary Element Method in Acoustics Modeling"
• Corrosion:"The Boundary Element Method Simplifies Corrosion Simulation"
• Magnetostatics:"Validating the Use of Boundary Elements for Magnetostatics Modeling"
• Electrostatics:"How to Create Electrostatics Models with Wires, Surfaces, and Solids"

Further Learning

The model used throughout this article is available for download. You can open and explore the model files to see exactly how the Laplace and Poisson's equations for the gravitational field of the Earth–Moon system were solved. In anotherLearning Center articlewe will examine the Coefficient Form PDE more.