Nonlinear tension-field elastic membrane#
\(\newcommand{\bA}{\boldsymbol{A}} \newcommand{\bB}{\boldsymbol{B}} \newcommand{\bC}{\boldsymbol{C}} \newcommand{\bE}{\boldsymbol{E}} \newcommand{\bCel}{\boldsymbol{C}_\text{el}} \newcommand{\bF}{\boldsymbol{F}} \newcommand{\bG}{\boldsymbol{G}} \newcommand{\bI}{\boldsymbol{I}} \newcommand{\bS}{\boldsymbol{S}} \newcommand{\bT}{\boldsymbol{T}} \newcommand{\bU}{\boldsymbol{U}} \newcommand{\bV}{\boldsymbol{V}} \newcommand{\bW}{\boldsymbol{W}} \newcommand{\bX}{\boldsymbol{X}} \newcommand{\bY}{\boldsymbol{Y}} \newcommand{\bZ}{\boldsymbol{Z}} \newcommand{\be}{\boldsymbol{e}} \newcommand{\bu}{\boldsymbol{u}} \newcommand{\bx}{\boldsymbol{x}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bz}{\boldsymbol{z}} \newcommand{\T}{^\text{T}} \newcommand{\CC}{\mathbb{C}} \newcommand{\Pp}{\mathcal{P}} \newcommand{\Qq}{\mathcal{Q}} \newcommand{\dOm}{\text{d}\Omega} \newcommand{\dS}{\text{d}S}\)
This demo explores a nonlinear elastic membrane formulation using the so-called “tension-field” theory developed by [REI38], [PIP86]. This theory enables to account for wrinkles occuring in compressed regions of a thin elastic membrane in an effective manner. Indeed, a unilateral constitutive law is formulated which assumes to stiffness of the elastic membrane in compressed regions.
Contents#
Hyperelastic membrane#
We consider the computation of the displacement field \(\bu\) of an hyperelastic membrane subject to imposed tractions \(\bT\) on the Neumann part of the boundary \(\partial\Omega_\text{N}\) and to fixed displacement \(\bu = 0\) on the Dirichlet part of the boundary \(\partial \Omega_\text{D}\). In the finite-strain setting, the hyperelastic potential of a 3D material can be expressed as a function of some nonlinear strain measures. For instance, let us consider the Cauchy-Green strain tensor \(\bC=\bF\T\bF\) where \(\bF=\bI+\nabla \bu\) is the deformation gradient. The free energy hyperelastic potential is then \(\psi(\bC)\) and the displacement field can be obtained as the solution to the following minimum principle:
In the general case, such a problem is not convex which makes the analysis of hyperelastic materials rather complex. In particular, some equilibrium positions can be unstable and lead to buckling phenomena.
Tension-field elastic membrane formulation#
As regards thin hyperelastic membranes, local buckling (or wrinkling) will occur at very low load levels in compressed regions. In the limit of infinitely thin membranes, compressed stress states cannot be supported at all. Tension field theory [REI38] has been proposed in order to simplify the analysis of thin membranes.
In the finite-deformation case, the tension-field theory has been first formalized by [PIP86] by introducing a relaxed strain energy functional. More precisely, introducing the following quasi-convexification of \(\psi\):
the tension field variational principle is obtained when replacing \(\psi\) with \(\psi_\text{memb}\) in (1):
The above relaxed potential then provides a tension-field constitutive equation in terms of the second Piola-Kirchhoff stress \(\bS\) as follows:
where \(\bC_\text{w}\) can be seen as an inelastic wrinkling strain accounting for the occurrence of wrinkles in compressed regions. As a consequence, the resulting stress is always tensile.
Conic reformulation#
Finally, Pipkin showed that when \(\psi\) is a convex function of \(\bC\), \(\psi_\text{memb}\) turns out to be a convex function of \(\bF\) or, equivalently, of the displacement gradient \(\bG=\nabla \bu\). Indeed, if \(\psi(\bCel)\) is convex, the relaxed minimum principle (2) is a convex program due to the following conic reformulation of the SDP constraint (see also [KAN11]):
Let us first recall the *Schur complement lemma for a PSD block-matrix:
A symmetric block-matrix
\[\begin{equation*} \bZ = \begin{bmatrix} \bU & \bV \\ \bV\T & \bW \end{bmatrix} \end{equation*}\]is positive semi-definite if and only if \(\bW\succeq 0\) and \(\bU-\bV\bW^{-1}\bV\T \succeq 0.\)
Let us then consider the following symmetric matrix:
we have \(\bW=\bI \succ 0\) and \(\bU-\bV\bW^{-1}\bV\T = \bCel - (\bI+\bG\T)(\bI+\bG) = \bCel - \bI - \bG - \bG\T - \bG\T\bG=\bC\). As a result, using the Schur complement lemma, \(\bZ\succeq 0\) if and only if \(\bCel \succeq \bC\).
In conclusion, if \(\psi(\bC)\) admits a convex conic representation in terms of \(\bC\), problem (2) is a convex conic program.
Ogden-type hyperelastic materials#
In the following, we consider an incompressible Ogden-type material (with only one term for simplicity) for which the 3D potential reads:
where \(\mu\) is the shear modulus, \(\alpha\) is some power-law exponent (which we assume in the following to be larger than 2) and \(\lambda_1,\lambda_2,\lambda_3\) denote the principal stretches. The latter are the positive square roots of the Cauchy-Green tensor eigenvalues. Note that the case \(\alpha=2\) corresponds to a neo-Hookean material.
In the incompressible case, we have \(J=\lambda_1\lambda_2\lambda_3=1\) which yields a reduced energy:
From the following relation for the principal stretches in 2D:
one can see that the reduced energy can be equivalently expressed as a convex function of \(\bC\) as follows:
where we assume that \(\beta=\alpha/2 \geq 1\).
The above expression can further be re-expressed introducing additional auxiliary variables as:
and where the last constraints can be formulated using a suitable quadratic and power cones as follows:
Surfaces embedded in \(\mathbb{R}^3\)#
The above presentation essentially considered planar surfaces. In general, membranes are embedded in a 3D space. The only modification to account for is the expression of the deformation gradient \(\bF\). In this setting, since we consider out-of-plane displacements, the displacement and deformation gradients \(\bG\) and \(\bF\) are in fact \(3\times 2\) matrices such that:
where \(\be_1\), \(\be_2\) denotes orthonormal basis vectors of the membrane initial reference plane. In the following implementation, \(\overline{\bI}\) will be refered to as the metric in the reference configuration.
The strain tensors \(\bC,\bCel\) are still \(2\times 2\) matrices. The only modification to apply to section \ref{sec:membrane-conic-reformulation} concerns the SDP constraint. We now have:
where \(\bI_n\) is the identity matrix of dimension \(n\). The latter constraint can now be equivalently reformulated as follows:
Implementation#
#!/usr/bin/env python
# coding: utf-8
"""Implementation of nonlinear membrane behaviours."""
from ufl import as_vector, Identity, shape
from fenics_optim import (
ConvexFunction,
Pow,
Quad,
SDP,
to_mat,
to_vect,
block_matrix,
)
class OgdenMembrane(ConvexFunction):
"""
Ogden model of a thin membrane.
The 2D potential is given by:
:math:`\\psi(\\lambda_1,\\lambda_2) = \\frac{2\\mu}{\\alpha^2}\\left(
\\lambda_1^\\alpha + \\lambda_2^\\alpha +
\\lambda_1^{-\\alpha}\\lambda_2^{-\\alpha}
\right)`
Parameters
----------
expr : UFL expression
displacement gradient
mu : float
shear modulus
alpha : float
Ogden exponent, should be > 2
metric : matrix, optional
metric tensor for embedded membranes in 3D
"""
def __init__(self, expr, mu, alpha, metric=None, **kwargs):
ConvexFunction.__init__(self, expr, parameters=(mu, alpha, metric), **kwargs)
def conic_repr(self, X, mu, alpha, metric):
assert float(alpha) > 2, "Exponent should be larger than 2 for convexity."
beta = alpha / 2
t = self.add_var(name="t")
s = self.add_var(dim=3, cone=Quad(3), name="s")
w = self.add_var(dim=3, cone=Quad(3))
x = self.add_var(dim=3, cone=Pow(3, 1 / beta))
y = self.add_var(dim=3, cone=Pow(3, 1 / beta))
z = self.add_var(dim=3, cone=Pow(3, 1 / (1 + 2 * beta)))
self.add_eq_constraint(x[1], b=1)
self.add_eq_constraint(y[1], b=1)
self.add_eq_constraint(z[2], b=1)
self.add_eq_constraint(x[2] - t - s[0])
self.add_eq_constraint(y[2] - t + s[0])
self.add_eq_constraint(w[0] - t)
self.add_eq_constraint(w[2] - s[0])
self.add_eq_constraint(w[1] - z[1])
Cel = to_mat(as_vector([t + s[1], t - s[1], s[2]]))
self.set_linear_term(2 * mu / alpha ** 2 * (x[0] + y[0] + z[0]))
if metric is None:
d1, d2 = 2, 2
metric = Identity(2)
else:
d2, d1 = shape(metric)
dtot = d1 + d2
Z = to_mat(self.add_var(dim=dtot * (dtot + 1) // 2, cone=SDP(dtot)))
Id2 = Identity(d2)
Id1 = Identity(d1)
A = Z + block_matrix([[-Cel, X.T], [X, 0 * Id2]])
B = block_matrix([[0 * Id1, -metric.T], [-metric, Id2]])
self.add_eq_constraint(to_vect(A), b=to_vect(B), name="Stress")
Acknowledgements#
We acknowledge Alexander Niewiarowski (Princeton University) for stimulating discussions on this topic.
References#
[REI38] E. Reissner, in Proceedings of the Fifth International Congress on Applied Mechanics (1938) pp. 88–92.
[PIP86] A. C. Pipkin, IMA Journal of Applied Mathematics, 36, 85 (1986)
[KAN11] Kanno, Y. (2011). Nonsmooth mechanics and convex optimization. Crc Press.