FEATool Multiphysics
v1.10 Finite Element Analysis Toolbox |

Mathematical Modeling

To model physical phenomena such as heat transfer, structural stresses, fluid dynamics, chemical reactions, wave propagation and more, a common approach is the continuum model which formulates the problems as partial differential equations. A typical model PDE can be derived by examining the rate of change of the involved quantities in an infinitesimal volume. This chapter introduces some mathematical theory for to derive FEM approach used by FEATool.

Consider a prototypical property or quantity , existing in an arbitrary domain . For this quantity to be conserved within , the rate of change of must be equal to the normal flux crossing the total surface (let denote the unit normal to ). With flux meaning the transport of , in and out of , with a vector field . Any generation from directional surface sources and scalar internal sources also contributes to an increased (or decreased) amount of . The following figure illustrates these phenomena for the region .

A balance equation explicitly stating the conservation of will have the form

which is easy to translate into mathematical vector notation

Application of Gauss' divergence theorem

to the surface terms results in the following reformulation

Let us consider the conservation of mass by substituting with the density . Under the assumption that the system is isolated it is possible to disregard any generation from internal and surface sources. What is left is the following balance equation

which must hold for a volume of any size. Thus, in the limit the well known mass conservation or continuity equation will appear, which reads

Moving on to conservation of linear momentum which, from Newton's second law, is the amount of mass transported with a certain velocity. Specifically, momentum transport per unit volume is therefore expressed as the density transported with the velocity field . If now the product of these quantities is substituted for in the general conservation equation then the flux of momentum over the boundaries will accordingly be represented by the quantity . Moreover, momentum is created and destroyed by forces acting internally and on the surface of . Inside of momentum is generated by the presence of forces, here labeled , while on the surface momentum is created by both pressure and shear stresses, which combined are represented by the stress tensor . Having made these observations enables the following formulation for the conservation of momentum

This should again be valid for a volume of any size or shape with the consequential removal of the integral operations, thus yielding the partial differential equation

The stress tensor can, as stated above, be subdivided into two parts, the static pressure and the shear stresses . Consider for the moment in the form of a unit cube where the coordinate axes are aligned with the surface normals. For a given coordinate axis, the pressure term will now only have one component which is acting directly against the surface normal, while the shear stress terms will have three different components, one in each coordinate direction. This enables a splitting of into the two terms and , that is

where denotes the diagonal unit matrix.

The Poisson equation is a very good model equation to study since it can represent many different physical processes such as diffusion, thermal conduction, and electric potential. The following example illustrates how to set up and solve the Poisson equation on a unit sphere ( ) with homogeneous boundary conditions, that is .

To derive the finite element discretization of the general Poisson equation first multiply the equation with an arbitrary function (called test function), and integrate over the whole domain , that is . By applying the Gauss theorem or partial integration to the left side we get

The boundary term can be neglected assuming that we only have prescribed or fixed value (Dirichlet) boundary conditions. In 2D the equation will look like

By discretizing the solution (trial function space) variable as and similarly test function , where are the finite element shape or basis functions (usually taken the same for and in the Galerkin approximation), and is the value of the solution at node . Inserting this we get

which after discretization and integration with numerical quadrature gives us a matrix system to solve

The modeled applications and processes are usually very specific in nature and the studies can thus be confined to a smaller spatial sub region or domain and a specific temporal extent . Both these restrictions will usually greatly simplify the modeling and also reduce the effort required to obtain a solution. It is particularly advantageous to utilize all existing symmetry axes to further shrink the computational domain. Sometimes it is also possible to transform a three dimensional model to two dimensions if the problem can be considered axisymmetric. The drawback of restricting the models spatially is that the equations now have to be supplied with suitable boundary conditions, which goal is to describe all interactions between and the rest of the non-modeled environment.

Boundary conditions usually come in two distinct variants, Dirichlet conditions which fixes the value of a quantity

and Neumann, also called natural, boundary conditions which specify the in- or outward-flux

Initial conditions, in the form of specified , must be prescribed in addition to the boundary conditions for applications where the temporal evolution is of interest. For non-linear stationary calculations it is also advantageous to specify a good initial guess so that the nonlinear solver will converge faster.

A tutorial example for the Poisson equation on a circular domain can be found in the Classic Equation - Poisson Equation on a Circle tutorial section. This example is split up into four parts showing how set up and solve the model in the GUI, as well as on the command line using physics modes, and the Poisson Equation using the CLI and core FEM library functions FEM core library functions such as matrix assembly.