Non-linear elasticity

Algorithm to solve the non-linear elasticity problem using a minimization technique.

Problem

Let $\Omega \subset \mathbb{R}^{3}$ denotes a hyperelastic material and $\partial \Omega= \Gamma_{0} \cup \Gamma_{1}$ its boundary, $\Gamma_{0}$ and $\Gamma_{1}$ denote respectively the disjoint parts of the boundary ($\Gamma_{0} \cap \Gamma_{1}= \emptyset$) where a null displacement and a surface traction $\mathbf{t}$ are applied.

The problem is to find the displacement field $\mathbf{u}$ of the body $\Omega$, which minimizes the total potential energy $\mathcal{E}$ given by:

$ \displaystyle{ \mathcal{E}(\mathbf{v})= \int_{\Omega} \Psi \, dV - \int_{\Gamma_{1}} \mathbf{t}.\mathbf{v} \, dA } $

Therefore:

$ \displaystyle{ \mathbf{u}= \underset{\mathbf{v} \in \mathcal{A}}{\text{argmin}}(\mathcal{E}(\mathbf{v})) } $

Where $\Psi$ is the strain energy function and $\mathcal{A}$ is the admissible displacements set defined by:

$ \displaystyle{ \mathcal{A}=\left\lbrace \mathbf{v} \in \left( H^{1}(\Omega)\right)^{3} \; ; \; \mathbf{v}=0 \; \text{on} \; \Gamma_{0} \right\rbrace } $

Algorithms

Material parameters

File NeoHookean.idp

real c10 = 1.2345;
real c01 = 0.;
real nu = 0.499;

real K = 6.*(c01+c10)/(3.*(1.-2.*nu));

// The energy density function is equal to : W(I1,I2,I3) = c1.I1.I3^(-1/3) + c2.I3 + c3.I3^(1/2) + c4

real c1 = c10;
real c2 = K/2.;
real c3 = -K;
real c4 = -3.*c10 + K/2.;

// The energy density
macro W(I)
(
	c1*I[0]*(I[2]^(-1./3.)) + c2*I[2] + c3*(I[2]^(1./2.)) + c4
) //

// The first differential of the energy density
macro dW(I, dI)
(
	c1*(I[2]^(-1./3.))*dI[0] + c1*I[0]*(-1./3.)*(I[2]^(-4./3.))*dI[2]	+ c2*dI[2]
		+ c3*(1./2.)*(I[2]^(-1./2.))*dI[2]
) //

// The second differential of the energy density
macro ddW(I, dI, dII)
(
	  c1*(-1./3.)*(I[2]^(-4./3.))*dI[0]*dII[2]
	+ c1*(-1./3.)*(I[2]^(-4./3.))*dI[2]*dII[0]
	+ c1*I[0]*(-1./3.)*(-4./3.)*(I[2]^(-7./3.))*dI[2]*dII[2]
	+ c3*(1./2.)*(-1./2.)*(I[2]^(-3./2.))*dI[2]*dII[2]
) //

Elasticity law

Definition of the elasticity law in 2D.

File ElasticLaw2d.idp

macro C2(d)
[
	1. + 2.*dx(d[0]) + dx(d[0])*dx(d[0]) + dx(d[1])*dx(d[1]),
	1. + 2.*dy(d[1]) + dy(d[0])*dy(d[0]) + dy(d[1])*dy(d[1]),
	 dy(d[0]) + dx(d[1]) + dx(d[0])*dy(d[0]) + dx(d[1])*dy(d[1])
] //
macro dC2(d, dd)
[
	2.*dx((dd)[0]) + 2.*dx((dd)[0])*dx(d[0]) + 2.*dx((dd)[1])*dx(d[1]),
	2.*dy((dd)[1]) + 2.*dy((dd)[0])*dy(d[0]) + 2.*dy((dd)[1])*dy(d[1]),
	dy((dd)[0]) + dx((dd)[1]) + dx((dd)[0])*dy(d[0]) + dx((dd)[1])*dy(d[1])
		+ dx(d[0])*dy((dd)[0]) + dx(d[1])*dy((dd)[1])
] //
macro ddC2(dd, ddd)
[
	2.*dx((dd)[0])*dx((ddd)[0]) + 2.*dx((dd)[1])*dx((ddd)[1]),
	2.*dy((dd)[0])*dy((ddd)[0]) + 2.*dy((dd)[1])*dy((ddd)[1]),
	dx((dd)[0])*dy((ddd)[0]) + dx((dd)[1])*dy((ddd)[1])
		+ dx((ddd)[0])*dy((dd)[0]) + dx((ddd)[1])*dy((dd)[1])
] //

macro I2C(C)
[
	C[0] + C[1] + 1.,
	C[0]*C[1] + C[1] + C[0] - C[2]*C[2],
	C[0]*C[1] - C[2]*C[2]
] //EOM

macro dI2C(C, dC)
[
	dC[0] + dC[1] ,
	dC[0]*C[1] + dC[1] + dC[0] - 2.*dC[2]*C[2] + C[0]*dC[1],
	dC[0]*C[1] + C[0]*dC[1] - 2.*C[2]*dC[2]
] //

macro  ddI2C(dC, ddC)
[
	0.*dC[0]*ddC[0],
	dC[0]*ddC[1] - 2.*dC[2]*ddC[2] + ddC[0]*dC[1],
	dC[0]*ddC[1] + ddC[0]*dC[1] - 2.*ddC[2]*dC[2]
] //

macro I2d(d) I2C(C2(d))  //
macro dI2d(d, dd) dI2C(C2(d), dC2(d, dd)) //
macro ddI2d(d, dd, ddd) (ddI2C(dC2(d, dd), dC2(d, ddd)) + dI2C(C2(d), ddC2((dd), (ddd)))) //
macro W2d(d) W(I2d(d)) //
macro dW2d(d, dd) dW(I2d(d), dI2d(d, dd)) //
macro ddW2d(d, dd, ddd) (ddW(I2d(d), dI2d(d, dd), dI2d(d, ddd)) + dW(I2d(d), ddI2d(d, dd, ddd))) //

Minimization algorithm

load "ff-Ipopt"

include "ElasticLaw2d.idp"
include "NeoHookean.idp"

// Dimension constants
real L = 0.5;
real l = 1.;
real cx = 0;
real cy = -l/2;

// Discretization constants
int Nx = 5;
int Ny = 2*Nx;
real f1 = 0, f2 = -0.876;
int[int] labs = [1, 2, 3, 4];
mesh Th = square(Nx, Ny, [L*x+cx, l*y+cy], label=labs, flags=1);

fespace Wh(Th, [P1, P1]);

//The total potential energy
func real iW (real[int] &D) {
	Wh [DDx, DDy];
	DDx[] = D;

	real res = int2d(Th)(W2d([DDx, DDy])) - int1d(Th,3)([DDx, DDy]'*[f1, f2]);
	return res;
}

//The gradient of the energy
func real[int] idW (real[int] &D) {
	Wh [DDx, DDy];
	DDx[] = D;

	varf dWW ([Dx, Dy], [Vx, Vy])
		= int2d(Th)(
			  dW2d([DDx, DDy], [Vx, Vy])) - int1d(Th, 3)([Vx, Vy]'*[f1, f2]);
	real[int] idWW = dWW(0, Wh);
	return idWW;
}

//The Hessian of the energy
matrix iddWW;
func matrix iddW (real[int] &D) {
	Wh [DDx, DDy];
	DDx[] = D;

	varf ddWW ([Dx, Dy], [Vx, Vy])
		= int2d(Th)(
			  ddW2d([DDx, DDy], [Dx, Dy], [Vx, Vy])
		);
	iddWW = ddWW(Wh, Wh);
	return iddWW;
}

// Boundary conditions
Wh [ub1, ub2] = [1e19, 1e19]; // Unbounded in interior
Wh [lb1, lb2] = [-1e19, -1e19]; // Unbounded in interior
varf vGamma1([Dx, Dy], [v1, v2]) = on(1, Dx=1, Dy=1); // Boundary conditions
varf vGamma4([Dx, Dy], [v1, v2]) = on(4, Dx=1, Dy=1); // Boundary conditions
real[int] onGamma1 = vGamma1(0, Wh);
real[int] onGamma4 = vGamma4(0, Wh);
Wh [ubb1, ubb2] = [1e19, 0.];
Wh [lbb1, lbb2] = [-1e19, 0.];
Wh [ubb4, ubb5] = [0., 1e19];
Wh [lbb4, lbb5] = [0., -1e19];
ub1[] = onGamma1 ? ubb1[] : ub1[]; // enforcing the boundary condition
lb1[] = onGamma1 ? lbb1[] : lb1[];
ub1[] = onGamma4 ? ubb4[] : ub1[]; // enforcing the boundary condition
lb1[] = onGamma4 ? lbb4[] : lb1[];
for (int i = 0; i < ub1[].n; i++)
	if (onGamma1[i] && onGamma4[i]){ ub1[][i] = 0.; lb1[][i] = 0.; }

Wh [Dx, Dy] = [0., 0.];

IPOPT(iW, idW, iddW, Dx[], lb=lb1[], ub=ub1[]); // Minimize

mesh Thm = movemesh(Th, [x+Dx, y+Dy]);
[Dx, Dy] = [Dx, Dy];

plot([Dx, Dy]);
Result - Deformation and displacement field
Displacement

Validation

See FRICTIONLESS CONTACT PROBLEM FOR HYPERELASTIC MATERIALS WITH INTERIOR POINT OPTIMIZER, H. Houssein, S. Garnotel, F. Hecht (upcoming article)

Authors

Author: Houssam Houssein