autotoc_md63

jupyter: jupytext: text_representation: extension: .md format_name: markdown format_version: '1.3' jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python

name: python3

Part 1: Theory

Introduction and motivation for the mixed finite element formulation

The theory discussed in this notebook has been developed as part of the finite element library, MoFEM [Kaczmarczyk, 2020], For further information visit our website here.

We will start this notebook by introducing the idea of the mixed formulation and its difference from the coupled problem formulation:

In the case of the coupled problem, we consider the interplay between two (or more) sub-problems governed by different physical equations and therefore described by two (or more) fields. Examples of coupled problems include thermoelasticity, electroelasticity, fluid-structure interaction and many others.
Mixed problem formulation is obtained upon introduction of one (or more) auxiliary variable(s) in a problem governed by one physical process.

Examples of classical problems permitting (and benefiting from) mixed formulation include, but are not limited to:

Incompressible elasticity (auxiliary variable: hydrostatic pressure)
Stokes/Navier-Stokes problem for viscous incompressible flow (auxiliary variable: fluid pressure)
Convection-Reaction-Diffusion problem (auxiliary variable: flux)
Mechanical contact problem (auxiliary variable: contact pressure).

There are multiple reasons for using mixed formulation, see [Boffi, 2013] for more details:

Presence of constraints in a problem under study (incompressible elasticity/fluid flow, contact problem)
Importance of new variables appearing in the formulation (accurate computation of stresses in elastic problem and fluxes in diffusion problem)
Possibility to obtain weaker formulation with less requirements on the regularity of the solution
Embedded reliable and efficient a posteriori error estimates.

Derivation of the mixed weak form for large strain elasticity

For a more detailed discussion of this formulation see this video.

Conservation law and boundary conditions

The strong form for quasi-static large strain elasticity is the classical conservation of linear momentum defined as:

\begin{align} -\textrm{DIV}\:(\mathbf P) = 0 & \quad \textrm{in}\; \Omega_0 \ \mathbf u - \bar {\mathbf u} = 0 & \quad \textrm{on}\; \Gamma_u, \ \bar{\mathbf t} - \mathbf P \cdot \mathbf N = 0 & \quad \textrm{on}\; \Gamma_t, \end{align}

where we assume no body forces acting on the domain, $\mathbf P$ represents the first Piola stress, $\mathbf u$ is the spatial displacement, $\bar{\mathbf u}$ is the prescribed spatial displacement boundary condition, $\bar{\mathbf t}$ is the applied traction and $\mathbf N$ is the normal vector of the respective surface.

Kinematics

To fully understand the mixed formulation in this notebook we assume a prior knowledge of non-linear continuum mechanics, but we will first introduce some lesser-used additional kinematics. Based on the polar decomposition of the deformation gradient $\mathbf{F} = \mathbf{R}\mathbf{U}$ where $\mathbf{U}$ is the symmetric stretch tensor and $\mathbf{R}$ is the orthonormal rotation tensor. The rotation can be expressed using the exponential map: $$ \mathbf{R} = \exp(\mathbf{A}) = \mathbf{1} + \frac{\sin(\omega)}{\omega} \mathbf{A} + \frac{2\sin^2(\omega/2)}{\omega^2} \mathbf{A}^2 $$

In index notation: $$ R_{iI} = \delta_{iI} + \frac{\sin(\omega)}{\omega} \varepsilon_{iIJ} \omega_J + \frac{2\sin^2(\omega/2)}{\omega^2} \varepsilon_{iKJ} \varepsilon_{KIL} \omega_J \omega_L $$

with: $$ \omega = |\boldsymbol{\omega}|, \quad \overline{\boldsymbol{\omega}} = \frac{\boldsymbol{\omega}}{\omega}, \quad \mathbf{A} = \boldsymbol{\varepsilon} \cdot \overline{\boldsymbol{\omega}} $$

Finally, the variation of spatial rotation is: $$ \delta R_{iI} = \varepsilon_{KIJ} R_{jK} \delta \omega_J $$

Energy Functional

To derive the mixed formulation we start with an energy functional below expressed in terms of the four independent fields as displacement vector $\boldsymbol{u}$, rotation tensor $\boldsymbol{R}$, stretch tensor $\boldsymbol{U}$ and first Piola-Kirchhoff stress tensor $\boldsymbol{P}$: \begin{equation} \Pi\left(\boldsymbol{u},\boldsymbol{R},\boldsymbol{U},\boldsymbol{P}\right) = W\left(\boldsymbol{R}\boldsymbol{U}\right) + \int_{\Omega_0} \left( \frac{\partial x_i}{\partial X_J} - R_{iK} U_{KJ} \right) P_{iJ} \, dV - \int_{\Gamma_0} \bar{u}_i t_i \, dS, \end{equation}

Consistency equation

We first obtain our consistency equation by taking variation with respect to first Piola-Kirchhoff stress as: \begin{align} \text{D } \Pi \left[\delta P_{iJ} \right] &= \int_{\Omega_0}\delta P_{iJ}\frac{\partial x_{i}}{\partial X_J}\text{ dV} - \int_{\Omega_0}\delta P_{iJ}\left(R_{iI} U_{IJ}\right) \text{ dV};\ &= \int_{\Omega_0}\delta P_{iJ}\frac{\partial u_i}{\partial X_J} \text{ dV} - \int_{\Omega_0} \delta P_{iJ}\left(R_{iI} U_{IJ}-\delta_{iJ}\right) \text{ dV}, \end{align}

Applying integration by parts to $\int_{\Omega_0}\delta P_{iJ}\frac{\partial u_i}{\partial X_J}\text{ dV}$: \begin{equation} \begin{split} \text{D } \Pi \left[\delta P_{iJ} \right] = - \int_{\Omega_0} \delta P_{iJ} \left( R_{iI} U_{IJ} - \delta_{iJ}\right) \text{ dV} - \int_{\Omega_0} \frac{\partial\delta P_{iJ}}{\partial X_J} u_i \text{ dV} +\int_{\Gamma_0} N_J \delta P_{iJ} \overline{u}_i\text{ dS}. \end{split} \end{equation}

Notice that the displacement $\bar u$ is now applied in the weak sense, which is typical for the mixed formulation. Traction boundary condition is applied a priori. This equation is used to enforce the consitency between our calculated displacements and solved displacements.

Physical equation

Next, the physical/constitutive equations are enforced in integral strong sense as follows: \begin{equation} \begin{split} \text{D } \Pi \left[(R\delta U)_{iJ} \right] &=\frac{\partial W \left(U\right)}{\partial \left( RU\right)_{iJ}} (R\delta U)_{iJ} - \int_{\Omega_0} (R\delta U)_{iJ}P_{iJ} \,\text{ dV} \ \text{D } \Pi \left[R_{iK}\delta U_{KJ} \right] &= R_{iK} \frac{\partial W (U)}{\partial U_{KJ}} R_{iK}\delta U_{KJ} - \int_{\Omega_0} R_{iK}\delta U_{KJ}P_{iJ} \,\text{ dV} \ \text{D } \Pi \left[\delta U_{KJ} \right] &=\frac{\partial W}{\partial U_{KJ}}\delta U_{KJ} - \int_{\Omega_0} R_{iK}\delta U_{KJ}P_{iJ} \,\text{ dV} \ \end{split} \end{equation} Note it is weighted only by bases of stretch tensor $\mathbf U$ since the rotations should be stress free. This equation enforces that the solved for first Piola-Kirchhoff stress is consistent with with calculated stress from the constitutive equation we choose, defined later.

Angular momentum equation

Next, angular momentum is satisfied with taking into consideration the rotation kinematics defined previously, as follows:

\begin{equation} \text{D } \Pi \left[\delta R_{iK} \right] = -\int_{\Omega_0} \left( \delta R_{iK} U_{KJ} P_{iJ}\right) \text{ dV} \end{equation} This equation ensures the conservation of angular momentum enforcing the symmetry of the cauchy stress tensor in a weak sense.

Linear momentum equation

Finally, the conservation of linear momentum equation is enforced as follows: \begin{equation} \text{D } \Pi \left[\delta u_i \right] = \int_{\Omega_0} \frac{\partial P_{iJ}}{\partial X_J}\delta u_i - \delta u_i f_i \, \text{ dV} \end{equation} This equation solves for the displacements of our problem.

Element equations

Combining these four equations gives a system of equations as follows:

\begin{equation*} \begin{split}

\int_{\Omega_0} \delta P_{iJ} \left( R_{iI} U_{IJ} - \delta_{iJ}\right) \text{ dV} - \int_{\Omega_0} \frac{\partial\delta P_{iJ}}{\partial X_J} u_i \text{ dV} +\int_{\Gamma_0} N_J \delta P_{iJ} \overline{u}_i\text{ dS} &= \boldsymbol{0};\ \frac{\partial W}{\partial U_{KJ}}\delta U_{KJ} - \int_{\Omega_0} R_{iK}\delta U_{KJ}P_{iJ} \,\text{ dV} &= \boldsymbol{0};\ -\int_{\Omega_0} \left( \delta R_{iK} U_{KJ} P_{iJ}\right) \text{ dV} &= \boldsymbol{0}; \ \int_{\Omega_0} \frac{\partial P_{iJ}}{\partial X_J}\delta u_i - \delta u_i f_i \, \text{ dV} &= \boldsymbol{0}; \end{split} \end{equation*}

Spaces

With above equations at hand, we choose approximation spaces as follows: $$ \left\{ \begin{array}{l} \delta\mathbf{P}\in U_0^h
\subset H^\textrm{div}_0(\Omega_0^h) :=
\delta\mathbf{P} \in H^\textrm{div}(\Omega_0^h) :
N_{j}\delta P_{ij} = 0\;\textrm{on}\;\Gamma_0^{h,u}