FEATool  v1.7
Finite Element Analysis Toolbox
 All Files Functions Pages
Physics Modes

In Equation mode partial differential equations (PDEs) together with equation coefficients must be chosen to accurately describe the physical phenomena to be simulated. Furthermore, in Boundary mode suitable boundary conditions must be prescribed in order to account for how the model interacts with its surroundings (outside of the modeled geometry and grid).

Equation and Boundary Modes

After creating a suitable geometry and grid one can switch between equation and boundary modes by using the and mode buttons. The corresponding Equation and Boundary Settings dialog boxes will also automatically open. Each boundary and equation setting tab corresponds to a different physics mode present in the model and allow for specifying the equation and boundary coefficients, initial conditions, and finite element shape functions. The coefficients and boundary coefficients will vary depending on the chosen physics mode and are explained in the following sections. The equation settings also allow for different equation and initial conditions to be set on a per subdomain basis by selecting the subdomains in the left hand side Subdomain: listbox.

Moreover, a convenient way to set up models is to use Model Constants and Expressions. The button opens the corresponding dialog box where one can enter constants and define expressions. One entered they can be used in stead of entering numerical values in equation coefficients and postprocessing expressions. More space for coefficients and expressions can also be added with the Add Row button.

phys_exprdlg_50.png


Multiphysics

The FEATool GUI also makes it easy to add and couple multiphysics equations and complex expressions to your models.

After starting to work with a physics mode it is simple to add one or more other modes by going to the last tab, with a + plus sign, in the Equation Settings dialog box. There you can simply choose an additional physics mode from the drop down combo box, select the dependent variable names you want to use (or keep the default ones), and press Add Physics >>> to add the mode. The new physics mode will now show up as a new tab with its short abbreviated name on the tab handle. Once the new mode has been added you can switch between the modes by clicking on the corresponding tabs in the Equation and Boundary Settings dialog boxes.

phys_addphys_50.png

The physics modes can be coupled by simply using the dependent variable names and derivatives in the coefficient expression dialog boxes. For example, the Navier-Stokes equations physics mode shown below uses the temperature variable T from the heat transfer mode in the source term for the y-direction.

featool-multiphysics-navier-stokes_50.jpg

Moreover, here below we can see how to make a three way multiphysics coupling. The convective velocities u and v are coupled from the Navier-Stokes equations physics mode and at the same time the temperature T and its two derivatives Tx and Ty are simultaneously coupled to the reaction rate source term in the convection and diffusion mode.

featool-multiphysics-convection-diffusion_50.jpg

As we have seen, it is very simple to set up multiphysics models in FEATool. This is made possible by the expression parsing functionality that allows you to enter and use complex expressions of dependent variables (for example u, v, T, c), their first derivatives (by just appending x or y to the variable names like Tx and Ty), the space dimensions x and y, as well as all common Octave and Matlab expressions and constants like pi, sin, cos, sqrt, ^2 etc.

Physics Modes

By using the predefined FEATool physics modes it is possible to easily and quickly implement models which simulate different physical effects such as fluid flow, structural stresses, chemical reactions, and heat transfer. This section describes how the various physics modes are defined, set up and used. The available physics modes are listed in the following table

Physics ModeDescriptionDefinition Function
Poisson EquationPoisson equationpoisson
Convection and DiffusionMass transport through convection and diffusionconvectiondiffusion
Conductive Media DCElectric potentialconductivemediadc
ElectrostaticsElectrostaticselectrostatics
MagnetostaticsMagnetostaticsmagnetostatics
Heat TransferHeat transport through convection and conductionheattransfer
Euler-Bernoulli BeamEuler-Bernoulli beam theory (1D)eulerbeam
Plane StressStructural mechanics (2D plane stress approximation)planestress
Plane StrainStructural mechanics (2D plane strain approximation)planestrain
Axisymmetric Stress-StrainStructural mechanics (2D cylindrical coordinates)axistressstrain
Linear ElasticityStructural mechanics (3D solid linear elasticity)linearelasticity
Darcy's LawDarcy's law porous media flowdarcyslaw
Brinkman EquationsBrinkman equations porous media flowbrinkmaneqns
Navier-Stokes EquationsIncompressible fluid flownavierstokes
Custom EquationUser defined equationcustomeqn

The equation parameters and coefficients for each selected physics mode are defined in the corresponding Equation Settings dialog boxes described below. In addition to also displaying the PDE equation, it is possible to change or edit the equation definitions with the edit equation button, and activate/deactivate physics modes in specific subdomains with the active button. Furthermore, the used finite element shape functions can also be selected either from the drop down box or directly entering a space separated list of shape functions for each dependent variable. Physics modes with convective effects allowing for numerical artificial stabilization also feature a button opening the stabilization settings.



Poisson Equation

The Poisson equation physics mode solves the classic elliptic Poisson equation for the scalar dependent variable $u$

\[ d_{ts}\frac{\partial u}{\partial t} + \nabla\cdot(-D\nabla u) = f \]

where $d_{ts}$ is a time scaling coefficient, $D$ is a diffusion coefficient, and $f$ is a scalar source term. In the Poisson Equation Settings dialog box shown below these equation coefficients, initial value $u_0 = u(t=0)$ can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

poi_eqdlg_50.png

The physics mode allows for both Dirichlet and Neumann (flux) boundary conditions. Dirichlet conditions prescribe a fixed value $r$ of the dependent variable $u = r$ on a boundary segment, while a Neumann condition will prescribe the normal flux to a boundary segment, that is $-\mathbf{n}\cdot D\nabla u = g$, where $\mathbf{n}$ is the outward directed normal, and $g$ therefore represents the value of the inward directed flux. The available boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
DirichletPrescribed value, $u = r$r
NeumannPrescribed flux, $-\mathbf{n}\cdot D\nabla u = g$g



Convection and Diffusion

The convection and diffusion physics mode models mass transport and reaction of a chemical species $c$. The governing equation for convection and diffusion reads

\[ d_{ts}\frac{\partial c}{\partial t} + \nabla\cdot(-D\nabla c) = R - \mathbf{u}\cdot\nabla c \]

where $d_{ts}$ is a time scaling coefficient, $D$ is a diffusion coefficient, $R$ is the reaction rate source term, and $\mathbf{u}$ a vector valued convective velocity field. In the Equation Settings dialog box shown below these equation coefficients, initial value $c_0 = c(t=0)$ can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

cd_eqdlg_50.png

For convection and diffusion problems with dominating convective effects it is advisable to use artificial stabilization. Pressing the lower Artificial Stabilization opens the corresponding dialog box which allows adding both isotropic artificial diffusion and anisotropic streamline diffusion. Turning coefficients are also provided to control the strength of the introduced artificial diffusion.

The convection and diffusion physics mode allows for four different boundary conditions; a prescribed concentration boundary condition, a convective flow (outflow condition), an insulation/symmetry condition which prescribes zero flux (or flow), and a prescribed flux condition. These boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
ConcentrationPrescribed concentration, $c = c_0$c0
Convective fluxOutflow, $-\mathbf{n}\cdot D\nabla c = 0$
Insulation/symmetryZero flux/flow, $-\mathbf{n}\cdot (D\nabla c + \mathbf{u}c) = 0$
Flux conditionPrescribed flux, $-\mathbf{n}\cdot (D\nabla c + \mathbf{u}c) = N_0$N0



Conductive Media DC

Electric potential and current can be modeled through the conductive media DC physics mode. To describe the flow of current and electric potential $V$ an electric field is firstly defined as $\mathbf{E} = \nabla V$. Furthermore, the current density $\mathbf{J}$ is related to the electric field by $\mathbf{J} = \sigma\mathbf{E}$ where $\sigma$ is the conductivity. By assuming conservation of the current one can pose a continuity equation for the current density, that is $\nabla\cdot\mathbf{J} = Q$ where $Q$ is a current source, which after expansion results in the following equation

\[ d_{ts}\frac{\partial V}{\partial t} + \nabla\cdot(-\sigma\nabla V) = Q \]

In the Equation Settings dialog box shown below the equation coefficients, initial value for the electric potential $V_0 = V(t=0)$ can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

dc_eqdlg_50.png

The conductive media DC physics mode features two boundary conditions, one prescribing the electric potential at a boundary $V = V_0$, and the other prescribing the current flow into the boundary. These boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
Electric potentialPrescribed electric potential, $V = V_0$V0
Current flowPrescribed current flow (flux), $-\mathbf{n}\cdot (\sigma\nabla V) = J_0$J0



Electrostatics

The electrostatics physics mode is an extension of the conductive media DC mode allowing for polarization vector, P, to be taken into account. The electrostatics equation for the electric potential V then reads

\[ -\nabla\cdot(\epsilon\nabla V - \mathbf{P}) = \rho \]

where $\epsilon$ is the permittivity and $\rho$ the charge density source term.

In the Equation Settings dialog box shown below the equation coefficients, initial value for the electric potential $V_0 = V(t=0)$ can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

es_eqdlg_50.png

The electrostatics physics mode features two boundary conditions, one prescribing the electric potential at a boundary $V = V_0$, and the other prescribing the current flow into the boundary. These boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
Electric potentialPrescribed electric potential, $V = V_0$V0
Ground/antisymmetryZero electric potential, $V = 0$
Surface chargePrescribed current flow (flux), $-\mathbf{n}\cdot (-\epsilon\nabla V + \mathbf{P}) = \rho_s$rhos
Insulation/symmetryZero current flow (flux), $-\mathbf{n}\cdot (-\epsilon\nabla V + \mathbf{P}) = 0$



Magnetostatics

The magnetostatics physics mode models simplified Maxwell's equations and in two dimensions solves

\[ -\nabla\cdot( \frac{1}{\mu}\nabla A_z - \mathbf{M} ) = J_z \]

for the magnetic vector potential $A_z$ in the z-direction where $\epsilon$ is the permeability, $\mathbf{M}$ the magnetization vector, and $J_z$ the current density. In three dimensions the equations are simplified to to exclude currents which reduces them to solving for a scalar magnetic potential $V_m$

\[ -\nabla\cdot( \mu\nabla V_m - \mu\mathbf{M} ) = 0 \]

In the Equation Settings dialog box shown below the equation coefficients, initial value for the magnetic potential can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

ms_eqdlg_50.png

The magnetostatics physics mode features four boundary conditions, one prescribing the magnetic potential at a boundary, one for magnetic insulation zero potential, surface current, and magnetic insulation/symmetry. These boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
Magnetic potentialPrescribed magnetic potential, $A_z/V_m = A_{z,0}/V_{m,0}$Az0/Vm0
Magnetic insulation/antisymmetryZero magnetic potential, $A_z/V_m = 0$
Surface currentPrescribed current flow (flux), $-\mathbf{n}\cdot (\frac{1}{\mu}\nabla A_z - \mathbf{M}/\mu\nabla V_m - \mu\mathbf{M}) = \mathbf{J}_s$Jx/y/z,s
Electric insulation/symmetryZero current flow (flux), $-\mathbf{n}\cdot (\frac{1}{\mu}\nabla A_z - \mathbf{M}/\mu\nabla V_m - \mu\mathbf{M}) = 0$



Heat Transfer (Convection and Conduction)

The heat transfer physics mode models heat transport through convection and conduction and heat generation through the following governing equation

\[ \rho C_p\frac{\partial T}{\partial t} + \nabla\cdot(-k\nabla c) = Q - \rho C_p\mathbf{u}\cdot\nabla T \]

where $\rho$ is the density, $C_p$ the heat capacity, $k$ is the thermal conductivity, $Q$ is the heat source term, and $\mathbf{u}$ a vector valued convective velocity field. In the Equation Settings dialog box show below the equation coefficients, initial value for the temperature $T_0 = T(t=0)$ can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

ht_eqdlg_50.png

For heat transfer problems with dominating convective effects it is advisable to use artificial stabilization. Pressing the lower Artificial Stabilization opens the corresponding dialog box which allows adding both isotropic artificial diffusion and anisotropic streamline diffusion. Turning coefficients are also provided to control the strength of the introduced artificial diffusion.

The heat transfer physics mode allows for four different boundary conditions; prescribed temperature, convective flow (outflow condition), an insulation/symmetry condition which prescribes zero flux (or flow), and a prescribed flux boundary condition.

The heat flux boundary condition involves several parameters. Firstly, an arbitrary expression for the heat flux may be prescribed with the coefficient q0. The second term, h*(Tinf-T), represents natural convection between the boundary and the surroundings. Here h is the convective heat transfer coefficient and Tinf a reference bulk temperature. The final term, Const*(Tamb^4-T^4), represents a radiation flux boundary condition where Const is the product between the emissivity of the boundary $\epsilon$ and the Stefan-Bolzmann constant $\sigma$, and Tamb is the surrounding ambient temperature.

The boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
TemperaturePrescribed temperature, $T = T_0$T0
Convective fluxOutflow, $-\mathbf{n}\cdot k\nabla T = 0$
Thermal insulation/symmetryZero heat flux, $-\mathbf{n}\cdot (k\nabla T + \rho C_p\mathbf{u}c) = 0$
Heat fluxPrescribed heat flux, $-\mathbf{n}\cdot (k\nabla T + \rho C_p\mathbf{u}T) = q_0 + h\cdot(T_{inf}-T) + Const\cdot(T_{amb}^4-T^4)$q0, h, Tinf, Const, Tamb



Euler-Bernoulli Beam

The Euler-Bernoulli beam physics mode models displacements, stresses and strains in a one-dimensional representations of beams and bars. The Euler-Bernoulli equations for the displacement $v$ reads

\[ \rho A \frac{\partial^2 v}{\partial t^2} + \frac{\partial^2}{\partial x^2}(EI\frac{\partial^2 v}{\partial x^2}) = q \]

where $\rho$ is the beam material density, A the cross-sectional area, E modulus of elasticity, I cross section moment of inertia, and q any distributed loads. In the Euler beam Equation Settings dialog box shown below these equation coefficients, initial value $v_0 = v(t=0)$ can be specified. The FEM shape function space can also be selected from the drop-down combobox (Hermite C1 shape functions by default).

eb_eqdlg_50.png

The boundary conditions for the Euler beam physics mode allows for any combination of prescribing displacements and edge loads. These boundary conditions are summarized in the table below.

Boundary ConditionDefinitionBoundary Coefficient
Fixed displacementsPrescribed displacements, $ u = u_0, v = v_0$v0
Edge loadsPrescribed edge loads, $ F_x = f_{x,0}, F_y = f_{y,0} $fy,0
eb_bdrdlg_50.png

Plane Stress

The plane stress physics mode models how structural stresses form in thin structures where the planar component of the stress can be neglected or considered zero. In this case the stress-strain relations can be written as

\[ \left( \begin{array}{c} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{array} \right) = \frac{E}{1-\nu^2} \left(\begin{array}{ccc} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{array} \right) \left( \begin{array}{c} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{array} \right) \]

where $E$ is the elastic or Young's modulus, and $\nu$ is the Poisson's ratio of the material. The strains are related to the material displacements ( $u$, $v$) as

\[ \epsilon_x = \frac{\partial u}{\partial x}, \qquad \epsilon_y = \frac{\partial v}{\partial y}, \qquad \gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \]

Balance equations for the stresses finally give the resulting governing equation system as

\[ \left\{\begin{aligned} - \frac{\partial\sigma_x}{\partial x} - \frac{\partial\tau_{xy}}{\partial y} = F_x \\ - \frac{\partial\tau_{xy}}{\partial x} - \frac{\partial\sigma_y}{\partial y} = F_y \end{aligned}\right. \]

where $F_x$ and $F_y$ are volume (body) forces in the x and y-directions, respectively. In the Equation Settings dialog box shown below the equation coefficients, initial value for the displacements can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

pss_eqdlg_50.png

Optional stress-strain temperature dependence can also be added, which takes the form

\[ \mathbf{\sigma} = \mathbf{D\epsilon} - \frac{E\alpha T}{1-\nu}[1\ \ 1\ \ 0]^T \]

where $\alpha$ is the coefficient of thermal expansion and $T$ either a prescribed temperature field, a dependent variable name from another physics mode that represents the temperature, or a combination such as $T-T_{ref}$.

The boundary conditions for the plane stress physics mode allows for any combination of prescribing displacements and edge loads. These boundary conditions are summarized in the table below.

Boundary ConditionDefinitionBoundary Coefficient
Fixed displacementsPrescribed displacements, $ u = u_0, v = v_0$u0, v0
Edge loadsPrescribed edge loads, $ F_x = f_{x,0}, F_y = f_{y,0} $fx,0, fy,0



Plane Strain

Like the plane stress physics mode, the plane strain mode models how structural stresses form in structures but where the z-component of the displacements can be neglected or considered zero. In this case the stress-strain relations can be written as

\[ \left( \begin{array}{c} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{array} \right) = \frac{E}{(1+\nu)(1-2\nu)} \left(\begin{array}{ccc} {1-\nu} & \nu & 0 \\ \nu & {1-\nu} & 0 \\ 0 & 0 & \frac{1-2\nu}{2} \end{array}\right) \left( \begin{array}{c} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{array} \right) \]

where $E$ is the elastic or Young's modulus, and $\nu$ is the Poisson's ration of the material. The strains are related to the material displacements ( $u$, $v$) as

\[ \epsilon_x = \frac{\partial u}{\partial x}, \qquad \epsilon_y = \frac{\partial v}{\partial y}, \qquad \gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \]

Balance equations for the stresses finally give the resulting governing equation system as

\[ \left\{\begin{aligned} - \frac{\partial\sigma_x}{\partial x} - \frac{\partial\tau_{xy}}{\partial y} = F_x \\ - \frac{\partial\tau_{xy}}{\partial x} - \frac{\partial\sigma_y}{\partial y} = F_y \end{aligned}\right. \]

where $F_x$ and $F_y$ are volume (body) forces in the x and y-directions, respectively. In the Equation Settings dialog box shown below the equation coefficients, initial value for the displacements can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

psn_eqdlg_50.png

Optional stress-strain temperature dependence can also be added, which takes the form

\[ \mathbf{\sigma} = \mathbf{D\epsilon} - \frac{E\alpha T}{(1+\nu)(1-2\nu)}[1\ \ 1\ \ 0]^T \]

where $\alpha$ is the coefficient of thermal expansion and $T$ either a prescribed temperature field, a dependent variable name from another physics mode that represents the temperature, or a combination such as $T-T_{ref}$.

The boundary conditions for the plane strain physics mode allows for any combination of prescribing displacements and edge loads. These boundary conditions are summarized in the table below.

Boundary ConditionDefinitionBoundary Coefficient
Fixed displacementsPrescribed displacements, $ u = u_0, v = v_0$u0, v0
Edge loadsPrescribed edge loads, $ F_x = f_{x,0}, F_y = f_{y,0} $fx,0, fy,0



Axisymmetric Stress-Strain

The axisymmetric stress-strain physics mode models how stresses and strains in cylindrical and rotationally symmetric geometries behave. The model equations are reduced from full 3D to a two-dimensional slice. In this case the balance equations for stress-strain (valid for the half plane r>0, with r=0 being the symmetry axis) can be written as

\[ \left\{\begin{aligned} - \frac{\partial\sigma_r}{\partial r} - \frac{\partial\tau_{rz}}{\partial z} - \frac{\sigma_r - \sigma_\theta}{r} = F_r \\ - \frac{\partial\tau_{rz}}{\partial r} - \frac{\partial\sigma_z}{\partial z} - \frac{\tau_{rz}}{r} = F_z \end{aligned}\right. \]

where $F_r$ and $F_z$ are volume (body) forces in the r and z-directions, respectively. The balance equations together with the constitutive relations

\[ \left( \begin{array}{c} \sigma_r \\ \sigma_\theta \\ \sigma_z \\ \tau_{rz} \end{array} \right) = \frac{E\nu}{(1+\nu)(1-2\nu)} \left(\begin{array}{cccc} \frac{1-\nu}{\nu} & 1 & 1 & 0 \\ 1 & \frac{1-\nu}{\nu} & 1 & 0 \\ 1 & 1 & \frac{1-\nu}{\nu} & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2\nu} \end{array} \right) \left( \begin{array}{c} \epsilon_r \\ \epsilon_\theta \\ \epsilon_z \\ \gamma_{rz} \end{array} \right) \]

uniquely define the equations to solve. Here $E$ is the elastic or Young's modulus, and $\nu$ is the Poisson's ratio of the material. The strains are related to the material displacements ( $u$, $w$) as

\[ \epsilon_r = \frac{\partial u}{\partial r}, \qquad \epsilon_\theta = \frac{u}{r}, \qquad \epsilon_z = \frac{\partial w}{\partial z}, \qquad \gamma_{rz} = \frac{\partial u}{\partial z} + \frac{\partial w}{\partial r} \]

In the PDE equation formulations the dependent variable u is replaced by u/r in order to avoid divisions by zero on the symmetry line (Thus custom postprocessing expressions including u should instead use r*u).

In the Equation Settings dialog box shown below the equation coefficients, initial value for the displacements can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

css_eqdlg_50.png

Optional stress-strain temperature dependence can also be added, which takes the form

\[ \mathbf{\sigma} = \mathbf{D\epsilon} - \frac{E\alpha T}{1-\nu}[1\ \ 1\ 1\ \ 0]^T \]

where $\alpha$ is the coefficient of thermal expansion and $T$ either a prescribed temperature field, a dependent variable name from another physics mode that represents the temperature, or a combination such as $T-T_{ref}$.

The boundary conditions for the plane stress physics mode allows for any combination of prescribing displacements and edge loads. These boundary conditions are summarized in the table below.

Boundary ConditionDefinitionBoundary Coefficient
Fixed displacementsPrescribed displacements, $ u = u_0, w = w_0$u0, w0
Edge loadsPrescribed edge loads, $ F_r = f_{r,0}, F_z = f_{z,0} $fr,0, fz,0



Linear Elasticity

The linear elasticity physics mode models how structural stresses form in solid structures where the stress-strain relations can be written as

\[ \left( \begin{array}{c} \sigma_x \\ \sigma_y \\ \sigma_z \\ \tau_{xy} \\ \tau_{yz} \\ \tau_{xz} \end{array} \right) = \frac{E}{(1+\nu)(1-2\nu)} \left(\begin{array}{cccccc} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \end{array} \right) \left( \begin{array}{c} \epsilon_x \\ \epsilon_y \\ \epsilon_z \\ \gamma_{xy} \\ \gamma_{yz} \\ \gamma_{xz} \end{array} \right) \]

where $E$ is the elastic or Young's modulus, and $\nu$ is the Poisson's ratio of the material. The strains are related to the material displacements ( $u$, $v$, $w$) as

\[ \epsilon_x = \frac{\partial u}{\partial x}, \qquad \epsilon_y = \frac{\partial v}{\partial y}, \qquad \epsilon_z = \frac{\partial w}{\partial z}, \qquad \gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}, \qquad \gamma_{yz} = \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}, \qquad \gamma_{xz} = \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \]

Balance equations for the stresses finally give the resulting governing equation system as

\[ \left\{\begin{aligned} - \frac{\partial\sigma_x}{\partial x} - \frac{\partial\tau_{xy}}{\partial y} - \frac{\partial\tau_{xz}}{\partial z} = F_x \\ - \frac{\partial\tau_{xy}}{\partial x} - \frac{\partial\sigma_y}{\partial y} - \frac{\partial\tau_{yz}}{\partial z} = F_y \\ - \frac{\partial\tau_{xz}}{\partial x} - \frac{\partial\tau_{yz}}{\partial y} - \frac{\partial\sigma_z}{\partial z} = F_z \end{aligned}\right. \]

where $F_x$, $F_y$, and $F_z$ are volume (body) forces in the x, y, and z-directions, respectively. In the Equation Settings dialog box shown below the equation coefficients, initial value for the displacements can be specified. The FEM shape function space can also be selected from the drop-down combobox (1st through 5th order conforming P1/Q1 shape functions), or further specified in the corresponding edit field.

le_eqdlg_50.png

Optional stress-strain temperature dependence can also be added, which takes the form

\[ \mathbf{\sigma} = \mathbf{D\epsilon} - \frac{E\alpha T}{1-2\nu}[1\ \ 1\ \ 1\ \ 0\ \ 0\ \ 0]^T \]

where $\alpha$ is the coefficient of thermal expansion and $T$ either a prescribed temperature field, a dependent variable name from another physics mode that represents the temperature, or a combination such as $T-T_{ref}$.

The boundary conditions for the linear elasticity physics mode allows for any combination of prescribing displacements and edge loads. These boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
Fixed displacementsPrescribed displacements, $ u = u_0, v = v_0, w = w_0$u0, v0 , w0
Edge loadsPrescribed edge loads, $ F_x = f_{x,0}, F_y = f_{y,0}, F_z = f_{z,0} $fx,0, fy,0, fz,0



Darcy's Law

Darcy's law models the pressure p in porous media flow through the relation

\[ d_{ts}\frac{\partial p}{\partial t} + \nabla\cdot(-\frac{\kappa}{\eta}\nabla p) = F \]

where $\kappa$ represents permeability, $\eta$ viscosity, F sources and sinks, and $d_{ts}$ is a time scaling coefficient. In the Darcy's Law Equation Settings dialog box shown below the equation coefficients, initial values for the pressure can be specified. The FEM shape function spaces can also be selected from the drop-down combobox, or further specified in the corresponding edit field.

dl_eqdlg_50.png

The Darcy's law physics mode allows prescription of three different boundary conditions. The prescribed pressure, flux, and insulation boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
PressurePrescribed pressure, $\mathbf{p} = p_0$p0
Insulation/symmetryZero flux/flow, $-\mathbf{n}\cdot(-\frac{\kappa}{\eta}\nabla p) = 0$
FluxPrescribed flux/flow, $-\mathbf{n}\cdot(-\frac{\kappa}{\eta}\nabla p) = N_0$N0



Brinkman Equations

The Brinkman equations physics mode models porous media flows and can be seen as a combination of Darcy's law with the Navier-Stokes equations. The Brinkman equations are defined as follows

\[ \left\{\begin{aligned} \rho\frac{\partial\mathbf{u}}{\partial t} - \eta\Delta\mathbf{u} + \frac{\eta}{\kappa}\mathbf{u} + \nabla p &= \mathbf{F} \\ \nabla\cdot\mathbf{u} &= 0 \end{aligned}\right. \]

which is to be solved for the unknown velocity field $\mathbf{u}$ and pressure $p$. In these equations $\rho$ represents the density of the fluid, $\eta$ the viscosity, and $\kappa$ permeability, moreover $\mathbf{F}$ represents body forces acting on the fluid. In the Equation Settings dialog box shown below the equation coefficients, initial values for the velocities and pressures can be specified. The FEM shape function spaces can also be selected from the drop-down combobox (P2P1/Q2Q1 and P1P-1/Q2P1 shape functions), or further specified in the corresponding edit field.

be_eqdlg_50.png

The Brinkman equations physics mode allows prescription of several different boundary conditions. Firstly, the no-slip (zero velocity) boundary condition which is appropriate for stationary walls. Moreover, a prescribed velocity condition can be prescribed to both in and outflows as well as moving walls. Prescribed pressure and neutral (zero viscous stress) conditions are both appropriate for outflows. Lastly, symmetry or slip conditions sets the flow to zero in one coordinate direction so as to prevent flow normal to the boundary. These boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
Wall/no-slipZero velocity, $\mathbf{u} = 0$
Inlet/velocityPrescribed velocity, $\mathbf{u} = (u_0, v_0, w_0)$u0, v0, w0
Neutral outflow/stress boundaryZero stress, $-\mathbf{n}\cdot (-p\mathbf{I} + \eta\nabla\mathbf{u}) = 0$
Outflow/pressurePrescribed pressure, $p = p_0$p0
Symmetry/slip, xi-directionZero velocity in xi-direction, $u_i = 0$



Navier-Stokes Equations (Incompressible Fluid Flow)

The Navier-Stokes equations physics mode models flows of incompressible fluid flows and which is described by

\[ \left\{\begin{aligned} \rho\left( \frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}\right) - \nabla\cdot(\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^T)) + p &= \mathbf{F} \\ \nabla\cdot\mathbf{u} &= 0 \end{aligned}\right. \]

which is to be solved for the unknown velocity field $\mathbf{u}$ and pressure $p$. In these equations $\rho$ represents the density of the fluid and $\mu$ the dynamic viscosity, moreover $\mathbf{F}$ represents body forces acting on the fluid. In the Equation Settings dialog box shown below the equation coefficients, initial values for the velocities and pressures can be specified. The FEM shape function spaces can also be selected from the drop-down combobox (P2P1/Q2Q1 and P1P-1/Q2P1 shape functions), or further specified in the corresponding edit field.

ns_eqdlg_50.png

For flow problems with dominating convective effects it is advisable to use artificial stabilization. Pressing the lower Artificial Stabilization opens the corresponding dialog box which allows adding both isotropic artificial diffusion and anisotropic streamline diffusion. Turning coefficients are also provided to control the strength of the introduced artificial diffusion.

The Navier-Stokes equations physics mode allows prescription of several different boundary conditions. Firstly, the no-slip (zero velocity) boundary condition which is appropriate for stationary walls. Moreover, a prescribed velocity condition can be prescribed to both in and outflows as well as moving walls. Prescribed pressure and neutral (zero viscous stress) conditions are both appropriate for outflows. Lastly, symmetry or slip conditions sets the flow to zero in one coordinate direction so as to prevent flow normal to the boundary. These boundary conditions are summarized in the table below

Boundary ConditionDefinitionBoundary Coefficient
Wall/no-slipZero velocity, $\mathbf{u} = 0$
Inlet/velocityPrescribed velocity, $\mathbf{u} = (u_0, v_0, w_0)$u0, v0, w0
Neutral outflow/stress boundaryZero stress, $-\mathbf{n}\cdot (-p\mathbf{I} + \mu(\nabla\mathbf{u}+\nabla\mathbf{u}^T)) = 0$
Outflow/pressurePrescribed pressure, $p = p_0$p0
Symmetry/slip, xi-directionZero velocity in xi-direction, $u_i = 0$



Custom Equation

User defined equations can be prescribed by using the custom equation physics mode.

ce_eqdlg_50.png

The equation specification can be accessed by either pressing the edit button next to the equation description, or the equation description itself. A dialog box will appear showing one edit field for equation and corresponding dependent variable.

ce_eqeditdlg_50.png

Note that the other physics mode equations are defined similarly and can also be edited by pressing the edit button in the corresponding tabs of the equation settings dialog box.

The syntax for equation specifications tries to as close as possible look like how one would write a partial differential equation with pen and paper. If for example the dependent variable is u like in the example above, then u' corresponds to the time derivative, ux the derivative in the x-direction, ux_x a second order derivative in the x-direction (to which partial integration will be applied according to the standard finite element derivation of the weak formulation), and u_t is the variable u multiplied with the fem test function, thus it will be assembled to the iteration matrix instead of the right hand side. The following table describes the syntax and legal operators

SyntaxDescriptionFormula
uDependent variable name$u$
xSpace dimension name$x$
uxDerivative in x-direction, rhs$\frac{\partial u}{\partial x}$
u_tDep. var multiplied with test function$u\cdot v$
u_xDerivative of test function$u\cdot \frac{\partial v}{\partial x}$
ux_tDerivative in x-direction$\frac{\partial u}{\partial x}\cdot v$
ux_x2nd derivative in x-direction$\frac{\partial u}{\partial x}\cdot\frac{\partial v}{\partial x}$
uxx_x3rd derivative in x-direction$\frac{\partial^2 u}{\partial x^2}\cdot\frac{\partial v}{\partial x}$
ux_xx3rd derivative in x-direction$\frac{\partial u}{\partial x}\cdot\frac{\partial^2 v}{\partial x^2}$
uxx_xx4th derivative in x-direction$\frac{\partial^2 u}{\partial x^2}\cdot\frac{\partial^2 v}{\partial x^2}$
+Addition
-Subtraction
*Multiplication
/Division
sqrt()Square root
^Power
()Delimit by enclosing in parentheses

where v is a test function. The equation syntax parser accepts numeric constants and coefficients defined in the fea.coef field. Higher dimensions work analogously and the space dimension x can be substituted with the others arbitrarily. For a more complicated example look at Custom Equation - Black-Scholes model equation in the tutorials section.

Artificial Stabilization

For problems with dominating convective effects such as can be found in convection and diffusion, heat transfer problems with convection, and fluid flow problems one can employ artificial stabilization if the grid size is too coarse to allow convergence.

The Artificial Stabilization dialog box allows control of both isotropic artificial diffusion and anisotropic streamline diffusion options. Turning coefficients are also provided to control the strength of the introduced artificial diffusion.

artstab_dlg_50.png

Isotropic artificial diffusion adds diffusion of magnitude delta*h_grid*|u| in all directions, where delta is the tuning coefficient, h_grid the local mean diameter of a grid cell, and |u| the magnitude of the convective velocity. Streamline diffusion modifies the finite element test function space and only adds a stabilization coefficient of the form delta*h_grid/|u| in the direction of the flow so as to minimize changes to the original problem.