Nonlinear FEM
Nonlinear Systems
Nonlinear problems in computational modeling are typically formulated by discretizing a field \(\phi\) into a vector \(\boldsymbol{\Phi}_n\) and defining a nonlinear residual \(\boldsymbol{r}\left(\boldsymbol{\Phi}_n\right)\) which for a viable solution returns
Let’s consider an introductory example
where \(\boldsymbol{f}\) is a constant right hand side. There is no general method to solve the above problem in closed form, meaning iterative methods are necessary.
Newton method and Picard iteration for General Nonlinear Systems
Start from the previous iteration \(i\) with values \(\boldsymbol{\Phi}_{i,n}\) and find updated solution \(\boldsymbol{\Phi}_{i+1,n}\) by approximating the nonlinear system with a Taylor expansion \(\boldsymbol{r}\) around \(\boldsymbol{\Phi}_{i,n}\) of first order and solve the local approximation
where \(\boldsymbol{J}\) is the jacobian matrix of the nonlinear residual \(\boldsymbol{r}\) with respect to \(\boldsymbol{\Phi}_{n}\):
\(\boldsymbol{J}\) can be found in terms by using the chain rule
where \(\nabla_{\boldsymbol{\Phi}_n} \boldsymbol{K}\left(\boldsymbol{\Phi}_{i,n}\right)\) is a third order tensor denoting entry-wise gradient of the matrix \(\boldsymbol{K}\) with respect to \(\boldsymbol{\Phi}_n\) which also acts on said vector via the single contraction \(\cdot\). We will see later that in many numerical schemes this simplifies greatly. Write residual at current iteration
and rearrange to
By repeatedly solving this linear problem, we will converge to the solution where \(\boldsymbol{r} \approx \boldsymbol{0}\). In practice this form is rarely chosen and further simplifications are done. We define the incremental update \(\boldsymbol{\Phi}_{incr,n}\)
after which the update of our solution changes to
and the linear problem simplifies to
which gets rid of the matrix-vector-product \(\boldsymbol{J} \boldsymbol{\Phi}_{i,n}\) and saves time. The procedure so far is the classical Newton method. An easier method often used in multiphysics is the Picard iteration where \(\boldsymbol{J}\) is further simplified to
It is often done as the element-wise gradient does not have to be calculated \(\nabla_{\boldsymbol{\Phi}_n} \boldsymbol{K}\left(\boldsymbol{\Phi}_{i,n}\right)\) and it often preserves the symmetry of \(\boldsymbol{K}\). It can basically be considered freezing \(\boldsymbol{K}\) at the current intermediate \(\boldsymbol{\Phi}_{i,n}\). The second part of the Jacobian is also often called the Newton correction
Newton method and Picard iteration in Numerical Schemes like FEM
In many common numerical schemes, \(\boldsymbol{K}\) as a function of the field \(\phi\) which in turn is a function of the discretized values \(\boldsymbol{K}\left(\phi \left(\boldsymbol{\Phi}_n\right)\right)\). Since in many common methods, \(\phi\left(\boldsymbol{\Phi}_n\right)\) is linear with respect to \(\boldsymbol{\Phi}_n\). Thus by applying the chain rule to the Newton correction (where the right hand side represents the already contracted form.)
we may rewrite the Newton correction
with \(\left(\nabla_{\boldsymbol{\Phi}_n} \phi \right)^T\) given by the discretization technique at hand and \(\frac{\partial \boldsymbol{K}\left(\phi\left(\boldsymbol{\Phi}_{i,n}\right)\right)}{\partial \phi}\) which e.g. in nonlinear heat conduction this would be the derivative of the heat conductivity tensor with respect to temperature. This is just a material property that has to be known/measured. For FE we remember the interpolation with shape functions \(\boldsymbol{N}\) from the previous chapter
therefore
Substituting this into Newton correction, we arrive at the final expression for FE:
Nonlin. Elasticity, stress and strain measures
In nonlinear mechanics, we distinguish between the (undeformed) reference configuration and the (deformed) current configuration. The motion \(\varphi\) maps material points from the reference position \(\boldsymbol{x} \in \Omega_0\) to the current position \(\varphi(\boldsymbol{x}) \in \Omega\). Depending on which configuration the balance laws are expressed in, different stress measures and their thermodynamic conjugate strain measures are used.
Stress measure |
Work-conjugate strain measure |
|---|---|
First Piola–Kirchhoff \( \boldsymbol{P} \) |
deformation gradient \(\boldsymbol{F}\) |
Second Piola–Kirchhoff \( \boldsymbol{S} \) |
Green–Lagrange strain \(\boldsymbol{E}\) |
Cauchy stress \( \boldsymbol{\sigma} \) |
The stress measures are related by
and
(a) Reference configuration with First Piola–Kirchhoff stress \(\boldsymbol{P}\)
Here \(\nabla_0\) is the gradient w.r.t. the reference coordinate \(\boldsymbol{x}\), \(\boldsymbol{b}_0\) is the body force per reference volume \(V_0\), and \(\bar{\boldsymbol{t}}_0\) is the nominal traction (per reference area).
(b) Reference configuration with Second Piola–Kirchhoff stress \(\boldsymbol{S}\)
We can rewrite the previous weak form also in terms of the 2nd Piola-Kirchhoff stress tensor \(\boldsymbol{S}\) using \(\boldsymbol{P}=\boldsymbol{F}\boldsymbol{S}\)
Equivalently, using the Green–Lagrange strain \(\boldsymbol{E}=\tfrac{1}{2}(\boldsymbol{F}^\mathrm{T}\boldsymbol{F}-\boldsymbol{I})\) and its variation \(\delta\boldsymbol{E}=\operatorname{sym}(\boldsymbol{F}^\mathrm{T}\nabla_0 \boldsymbol{w})\):
(c) Current configuration with Cauchy stress \(\boldsymbol{\sigma}\)
We can write weak form also in terms of the Cauchy stress tensor \(\boldsymbol{\sigma}\) using \(\boldsymbol{P}=\boldsymbol{F}\boldsymbol{S}\)
where \(\nabla\) is the spatial gradient w.r.t. \(\phi(x)\), \(\boldsymbol{b}\) is body force per current volume \(V\), and \(\bar{\boldsymbol{t}}=\boldsymbol{\sigma}\boldsymbol{n}\) is the Cauchy traction on the Neumann boundary \(\Gamma_N\) (with outward normal \(\boldsymbol{n}\)). On the Dirichlet boundary \(\Gamma_D\), \(\boldsymbol{w}=\boldsymbol{0}\).