import theme
theme.load_style()
This lecture by Tim Fuller is licensed under the Creative Commons Attribution 4.0 International License. All code examples are also licensed under the MIT license.
import numpy as np
import sympy as sp
%pylab inline
sp.init_printing("latex")
Populating the interactive namespace from numpy and matplotlib
Nonlinearity arises naturally in physical systems and can be modeled by finite elements in one (or more) of several ways.
They types of nonlinearity encountered in static analysis can be classified as either material nonlinearity or geometric nonlinearity.
To illustrate, consider the case of nonlinear elasticity, $$ \sigma = f\left(\epsilon\right) $$
where $\epsilon$ is a strain determined from the stretch $\lambda=L/L_0$
$$ \epsilon = g\left(\lambda\right) $$Here, $f$ and $g$ are nonlinear functions. Combined, we get
$$ \sigma = f\left(g\left(\lambda\right)\right) $$Suppose laboratory experiments show
lam = np.linspace(.1, 5)
sig = 38 * np.log(lam)
plt.xlabel(r"$\lambda$")
plt.ylabel(r"$\sigma$")
plot(lam, sig);
There are many ways we can express the $f$ and $g$ functions:
$$\begin{align} \sigma &= f\left(\epsilon\right) = 38 \epsilon \\ \epsilon &= g\left(\lambda\right) = \ln{\lambda} \end{align}$$or
$$\begin{align} \sigma &= f\left(\epsilon\right) = 38 \ln{\left(\epsilon + 1\right)} \\ \epsilon &= g\left(\lambda\right) = \lambda - 1 \end{align}$$or
$$\begin{align} \sigma &= f\left(\epsilon\right) = 19 \ln{\left(2\epsilon + 1\right)} \\ \epsilon &= g\left(\lambda\right) = \frac{1}{2}\left(\lambda^2 - 1\right) \end{align}$$etc.
When combined, all of these choices give
$$ \sigma = f\left(g\left(\lambda\right)\right) = 38 \ln{\lambda} $$Choice 1 has a linear stress-strain relationship but a nonlinear strain-stretch relationship. This is called geometric nonlinearity.
Choice 2 has a nonlinear stress-strain relationship but a linear strain-stretch relationship. This is called material nonlinearity.
Choice 3 has both material and geometric nonlinearity.
No matter the choice of where to put nonlinearity, the end result of stress as a function of stretch should be the same. Thus, the terms material nonlinearity and geometric nonlinearity are somewhat arbitrary. They do indicate differences in the FEM implementation, but not differences in the final answer.
In solid mechanics, static equilibrium requires that
$$ \frac{d\sigma\left[u\right]}{dx} + b = 0 $$where we have now that $\sigma\left[u\right]$ is a (not necessarily linear) function of the displacement $u$ and and $b$ is a force per unit volume.
The weak form is derived in the usual way
$$\begin{align} &\int_{\Omega}\left(\frac{d\sigma\left[u\right]}{dx} + b\right)w\,dx = 0, \qquad \forall w \\ &\int_{\Omega}\left(\sigma\left[u\right]\frac{dw}{dx} - bw\right)\,dx - \sigma\left[u\right] w\Bigg|_{\partial \Omega} = 0 \end{align} $$Because stress is not necessarily linear in $u$, $u$ must be found iteratively.
The finite element solution procedure for the nonlinear problem follows nearly identical that of the linear problem, namely, the isplacement field is approximated by interpolating polynomials possessing the property of compact support (shape functions):
$$ u = \sum_{j=1}^n u_j N_j $$and the weighting functions as
$$ w = \sum_{i}^n w_i N_i $$The stress corresponding to a given displacement is found by computing the strain
$$ \epsilon = \frac{du}{dx} = \sum_{j=1}^{n}u_j \frac{dN_j}{dx} $$and then passing the strain to the material model $\sigma = \sigma\left[\epsilon\right]$
Upon substitution of the preceding equations in to the weak form we arrive at
$$ \int_{\Omega}\left(\sigma\left[\epsilon(u_j)\right]\frac{dN_i}{dx} - bN_i\right)\,dx - \sigma\left[\epsilon(u_j)\right] N_i\Bigg|_{\partial \Omega} = 0, \quad \forall x_i \text{ on } \partial_\sigma \Omega $$$$ u_i = u_i^*\left(x_i\right), \quad \forall x_i \text{ on }\partial_u\Omega $$This is a set of $n$ equations for the $n$ unknown nodal displacements $u_i$. This set of equations is similar to those in the linear case, but because the the consititutive law relating stress and strain is not necessarily linear, it must be solved numerically.
One method of solving the finite element equations is by Newton-Raphson iteration whereby we linearize the finite element equations about some guess for the $u_j$ and obtain successfully better approximations by iteration. The linearization of $\sigma$ about $u_j$ is
$$ \sigma\left[u_k\right] \approx \sigma\left[u_k^0\right] + \frac{d\sigma}{du_k}du_k = \sigma\left[u_k^0\right] + \frac{d\sigma}{d\epsilon}\frac{d\epsilon}{du_k}du_k $$recalling that
$$ \epsilon = \sum_{k=1}^n u_k \frac{dN_k}{dx} \rightarrow \frac{d\epsilon}{du_j} = \sum_{k=1}^n \frac{du_k}{du_j} \frac{dN_k}{dx} = \sum_{k=1}^n \delta_{kj} \frac{dN_k}{dx} = \frac{dN_j}{dx} $$we can write
$$ \int_{\Omega}\left(\frac{d\sigma}{d\epsilon}\frac{dN_i}{dx}\frac{dN_j}{dx}\,dx\right) du_j + \int_{\Omega}\sigma\frac{dN_i}{dx}\, dx - \int_{\Omega}bN_i\,dx - \sigma N_i\Bigg|_{\partial \Omega} = 0 $$or
$$ K_{ij}du_{j} = F_{i} - R_{i} $$where $$ K_{ij} = \int_{\Omega}\frac{d\sigma}{d\epsilon}\frac{dN_i}{dx}\frac{dN_j}{dx}\,dx, \quad R_{i} = \int_{\Omega}\sigma\frac{dN_i}{dx}\, dx, \quad F_{i} = \int_{\Omega}bN_i\,dx + \sigma N_i\Bigg|_{\partial \Omega} $$
This expression is nearly identical to the linear case with the addition of two terms:
STEP 0 Initial guess for $u_j$ ($u_j=0$ is customary)
STEP 1 Compute global $K_{ij}$, $R_i$, and $F_i$
STEP 2 Enforce boundary conditions
STEP 3 Solve $K_{ij}du_j = F_i - R_i$
STEP 4 Set $u_j \leftarrow u_j + du_j$
STEP 5 If $\lVert du_j\rVert>$ TOL: GO TO STEP 1
STEP 6 Write output, postprocess, etc.
Let us investigate what additional information we will need to compute the element residual. The element residual is given by
$$ r_k^e = \int_{x_A}^{x_B}\sigma \frac{dN_k}{dx} dx $$Using Gauss-Legendre quadrature
$$ r_k^e = \int_{-1}^{1}\sigma \frac{1}{J} \frac{dN_k}{d\xi} Jd\xi = \sum_{l=1}^{n}\sigma\frac{dN_k\left(\xi_l\right)}{d\xi} w_l $$Computing the residual will require computing the stress at a Gauss point and, implicitly, the strain since $\sigma = f(\epsilon)$. Recall tha the strain is given by
$$ \epsilon = \sum_{i=1}^{N} u_i \frac{dN_i}{dx} = \sum_{i=1}^{N} \frac{1}{J} u_i \frac{dN_i}{d\xi} $$data
external functions
procedure
Step 1: Update nodal coordinates $$x_k=X_k+u_k$$
Step 2: Initialize the element residual $$r_k^e=0$$
Step 3: Perform Gauss integration
for $n$ from 1 to number of Gauss points do
Step a: Compute shape function derivative wrt to $\xi$ and the Jacobian $$\frac{dN_k}{d\xi} = \frac{dN_k\left(\xi_n\right)}{d\xi}$$ $$J = J\left(x_k, \xi_n\right)$$
Step b: Compute the strain at $\xi_n$ $$\epsilon = \sum_k \frac{1}{J} u_k \frac{dN_k\left(\xi_n\right)}{d\xi}$$
Step c: Compute the stress at $\xi_n$ $$\sigma = f\left(\epsilon\right)$$
Step d: Update element residual
for $j$ from 1 to number of element nodes do
$$r^e_j = r^e_j + \sigma\frac{dN_j}{d\xi} w_n$$
end for
end for
The element stiffness calculation follows as before, except that instead of a constant element stiffness, the stiffness may change with deformation
data
external functions
procedure
Step 1: Update nodal coordinates $$x_k=X_k+u_k$$
Step 2: Initialize the element stiffness $$k_{ij}^e=0$$
Step 3: Perform Gauss integration for $n$ from 1 to number of Gauss points do
Step a: Compute shape function derivative wrt to $\xi$ and the Jacobian $$\frac{dN_k}{d\xi} = \frac{dN_k\left(\xi_n\right)}{d\xi}$$ $$J = J\left(x_k, \xi_n\right)$$
Step b: Compute the strain at $\xi_n$ $$\epsilon = \sum_k \frac{1}{J} u_k \frac{dN_k\left(\xi_n\right)}{d\xi}$$
Step c: Compute the element stiffness at $\xi_n$ $$\frac{d\sigma}{d\epsilon} = f\left(\epsilon\right)$$
Step d: Update element stiffness
for $i$ from 1 to number of element nodes do
for $j$ from 1 to number of element nodes do
$$k_{ij}^e = k_{ij}^e +
\frac{1}{J}\frac{d\sigma}{d\epsilon}\frac{dN_i}{d\xi}\frac{dN_i}{d\xi} w_n$$
end for
end for
end for
An example of a nonlinear material model, widely used in the study of polymer-like materials subjected to moderate strains, is the single parameter Neo-Hookean material model whereby the stress is related to the stretch by
L, L0, lam, e, E = sp.symbols("L L0 lambda epsilon E")
s = E / 2 * (lam - 1 / lam)
s
where $E$ is Young's modulus. In terms of engineering strain
$$ \epsilon = \frac{L - L_0}{L_0} = \lambda - 1 \rightarrow \lambda = \epsilon + 1 $$the stress is
s = s.subs({lam: e + 1})
s.simplify()
The stiffness, $\frac{d\sigma}{d\epsilon}$ is given by
dsde = sp.diff(s, e)
dsde.simplify()
Let's make sure that the stiffness at $\epsilon=0$ is $E$
dsde.subs({e:0})
Implement a Newton solver using a linear elastic material model, verify it gives the same answer as the direct solver.
Doing so will require calling an external material model whenever the stress and/or material stiffness are needed. To reduce complications and repeated code, write a function like:
def update_material_state(props, strain, mat_id=1):
if mat_id == 1:
# linear elastic
E = props[0]
stress = E * strain
stiff = E
elif mat_id == ...
return stress, stiff
Different material models can then be implemented by giving them a mat_id, i.e., mat_id will now be a model input. The FEModel constructor should also make sure that the props array is set up correctly for the material.
When element stiffness and/or stress are needed, just call update_material_state as
el_stress, el_stiff = update_material_id(props, strain, mat_id)
Of course, you'll need to make sure that props, strain, and mat_id are correct.