Example elements

StationaryFourier(k)

For solving stationary linear heat conduction (which uses Fourier's law) with conductivity k, such that the heat flux is $\boldsymbol{q}=-k \nabla T$, where $T$ is the temperature field.

The strong form is,

\[ \nabla \cdot \boldsymbol{q} = h \quad \textbf{x} \in \Omega,\]

and the corresponding weak form is

\[ -\int_\Omega [\nabla \delta T] \cdot \boldsymbol{q}\ \mathrm{d}\Omega = - \int_\Gamma \delta T\ q_\mathrm{n}\ \mathrm{d}\Gamma + \int_\Omega \delta T\ h\ \mathrm{d}\Omega\]

where, on the right hand side, $q_\mathrm{n}$ is a heat flux normal to the boundary $\Gamma$, and $h$ is a volumetric heat supply. These contributions are not included in the element, and should be added with FerriteNeumann.jl

TransientFourier(k, c)

For solving the transient linear heat equation (which uses Fourier's law) with conductivity k, such that the heat flux is $\boldsymbol{q}=-k \nabla T$, where $T$ is the temperature field.

The strong form is,

\[ c\ \dot{T} + \nabla \cdot \boldsymbol{q} = h \quad \textbf{x} \in \Omega,\]

and the corresponding time-discretized weak form is

\[ \int_\Omega \delta T\ c\ \frac{T - {}^\mathrm{n}T}{\Delta t}\ \mathrm{d}\Omega - \int_\Omega [\nabla \delta T] \cdot \boldsymbol{q}\ \mathrm{d}\Omega = - \int_\Gamma \delta T\ q_\mathrm{n}\ \mathrm{d}\Gamma + \int_\Omega \delta T\ h\ \mathrm{d}\Omega\]

where ${}^\mathrm{n}T$ is the old temperature (in the previous timestep) and $\Delta t$ is the timestep. On the right hand side, $q_\mathrm{n}$ is a heat flux normal to the boundary $\Gamma$, and $h$ is a volumetric heat supply. These external contributions on the right hand side are not included in the element, and should be added with FerriteNeumann.jl

ElasticPlaneStrain(;E=2.e3, ν=0.3)

For solving linear elasticity for plane strain, where Young's modulus, E, and Poisson's ratio, ν, is used to construct the correct stiffness tensor, $\boldsymbol{\mathsf{C}}$, such that the stress, $\boldsymbol{\sigma}=\boldsymbol{\mathsf{C}}:\boldsymbol{\epsilon}$, where the strain tensor, $\boldsymbol{\epsilon}=[\boldsymbol{u}\otimes\nabla]^\mathrm{sym}$, is calculated from the displacement field, $\boldsymbol{u}(\boldsymbol{x},t)$.

The strong form of the mechanical quasi-static equilibrium is

\[ \boldsymbol{\sigma} \cdot \nabla + \boldsymbol{b} = 0 \quad \textbf{x} \in \Omega,\]

with the corresponding weak form,

\[ \int_\Omega [\boldsymbol{\delta u}\otimes\nabla]^\mathrm{sym} : \boldsymbol{\sigma} \mathrm{d}\Omega = \int_\Gamma \boldsymbol{\delta u} \cdot \boldsymbol{t} \mathrm{d}\Gamma + \int_\Omega \boldsymbol{\delta u} \cdot \boldsymbol{b} \mathrm{d}\Omega\]

The external loading on the right hand side is not included in the element, but can be implemented using FerriteNeumann.jl. (Note that this constructor returns LinearElastic, which just stores the correct stiffness tensor for the case of isotropic plane strain)

PoroElasticPlaneStrain(;E=2.e3, ν=0.3, k=0.05, α=1.0, β=1/2e3)

The strong forms are given as

\[\begin{aligned} \boldsymbol{\sigma}(\boldsymbol{\epsilon}, p) \cdot \boldsymbol{\nabla} &= \boldsymbol{0} \\ \dot{\Phi}(\boldsymbol{\epsilon}, p) + \boldsymbol{w}(p) \cdot \boldsymbol{\nabla} &= 0 \end{aligned}\]

where $\boldsymbol{\epsilon} = \left[\boldsymbol{u}\otimes\boldsymbol{\nabla}\right]^\mathrm{sym}$ The constitutive relationships are

\[\begin{aligned} \boldsymbol{\sigma} &= \boldsymbol{\mathsf{E}}:\boldsymbol{\epsilon} - \alpha p \boldsymbol{I} \\ \boldsymbol{w} &= - k \boldsymbol{\nabla} p \\ \Phi &= \phi + \alpha \mathrm{tr}(\boldsymbol{\epsilon}) + \beta p \end{aligned}\]

with $\boldsymbol{\mathsf{E}}=2G \boldsymbol{\mathsf{I}}^\mathrm{dev} + 3K \boldsymbol{I}\otimes\boldsymbol{I}$. The material parameters are then the shear modulus, $G$, bulk modulus, $K$, permeability, $k$, Biot's coefficient, $\alpha$, and liquid compressibility, $\beta$. The porosity, $\phi$, doesn't enter into the equations (A different porosity leads to different skeleton stiffness and permeability).

The weak forms are

\[\begin{aligned} \int_\Omega \left[\left[\boldsymbol{\delta u}\otimes\boldsymbol{\nabla}\right]^\mathrm{sym}: \boldsymbol{\mathsf{E}}:\boldsymbol{\epsilon} - \boldsymbol{\delta u} \cdot \boldsymbol{\nabla} \alpha p\right] \mathrm{d}\Omega &= \int_\Gamma \boldsymbol{\delta u} \cdot \boldsymbol{t} \mathrm{d} \Gamma \\ \int_\Omega \left[\delta p \left[\alpha \dot{\boldsymbol{u}} \cdot \boldsymbol{\nabla} + \beta \dot{p}\right] + \boldsymbol{\nabla}(\delta p) \cdot [k \boldsymbol{\nabla}]\right] \mathrm{d}\Omega &= -\int_\Gamma \delta p w_\mathrm{n} \mathrm{d} \Gamma \end{aligned}\]

where $\boldsymbol{t}=\boldsymbol{n}\cdot\boldsymbol{\sigma}$ is the traction and $w_\mathrm{n} = \boldsymbol{n}\cdot\boldsymbol{w}$ is the normal flux.

WeakForm(f::Function)

Solve the problem with primary variable, $u$, variation, $\delta u$, and a weak form

\[ \int_\Omega f(\delta u, \delta u \otimes\nabla, u, u \otimes\nabla, \dot{u}, \dot{u} \otimes\nabla) \mathrm{d}\Omega - \int_\Gamma \delta u\ h(x,t,n)\ \mathrm{d}\Gamma - \int_\Omega \delta u\ b(x,t)\ \mathrm{d}\Omega = 0\]

where the function f is given to the weak form, and h and b are given with FerriteNeumann.

Examples

Transient heat flow

Weak form

\[ \int_\Omega \delta u\ c\ \dot{u} + k\ [\nabla \delta u] \cdot [\nabla u]\ \mathrm{d}\Omega + \int_\Gamma \delta u\ q_\mathrm{n}\ \mathrm{d}\Gamma - \int_\Omega \delta u\ b\ \mathrm{d}\Omega = 0\]

Implementation This implementation is equivalent to TransientFourier, but it is also possible to add the body load directly in the weak form (if desired).

c = 1.0; k = 1.0; # heat capacity and heat conductivity (material parameters)
qn = 1.0; b=1.0;  # Normal boundary flux and internal heat source (external loading)
material = WeakForm((δu, ∇δu, u, ∇u, u_dot, ∇u_dot) -> δu*c*u_dot + k*(∇δu ⋅ ∇u))
nh = NeumannHandler(dh)
add!(nh, Neumann(:u, 2, getfacetset(dh.grid, "right"), (x,t,n)->qn))
add!(nh, BodyLoad(:c, 1, (x,t)->b))

Linear elasticity

This implementation is equivalent to ElasticPlaneStrain, but it is also possible to add the body load directly in the weak form (if desired). Weak form

\[ \int_\Omega [\boldsymbol{\delta u}\otimes\nabla]^\mathrm{sym} : \boldsymbol{\sigma}\ \mathrm{d}\Omega - \int_\Gamma \boldsymbol{\delta u} \cdot \boldsymbol{t}\ \mathrm{d}\Gamma - \int_\Omega \boldsymbol{\delta u} \cdot \boldsymbol{b}\ \mathrm{d}\Omega = 0\]

Implementation

G = 80e3; K = 160e3; # Shear and bulk modulus (material parameters)
tn = 1.0, b=Vec((0.0, 0.0, -1.0)); # Normal traction and body force (external loading)
material = WeakForm((δu, ∇δu, u, ∇u, u_dot, ∇u_dot) -> (∇δu ⊡ (2*G*dev(symmetric(∇u)) + 3*K*vol(∇u))))
nh = NeumannHandler(dh)
add!(nh, Neumann(:u, 2, getfacetset(dh.grid, "right"), (x,t,n)->tn*n))
add!(nh, BodyLoad(:c, 2, (x,t)->b))

J2Plasticity(;E, ν, σ0, H) <: MaterialModelsBase.AbstractMaterial

This plasticity model is taken from Ferrite.jl's plasticity example, and considers linear isotropic hardening, with Young's modulus, E, Poisson's ratio, ν, initial yield limit, σ0, and hardening modulus, H. It is defined as an AbstractMaterial following the MaterialModelsBase interface.