Bilayer Wrinkling

The problem provided in this example is a bilayer wrinkling problem. The domain is a thin stiff film resting on top of a substrate. A small point-force perturbation is applied to the substrate to instigate the instability. The domain is compressed well after the onset of instability.

%matplotlib inline
from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import os
from arc_length.displacement_control_solver import displacement_control # import displacement control formulation of arc-length solver


parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

Define domain and function spaces

We also split the domain into substrate and film. We To make sure all the assumptions of a bilayer system are valid, we construct the doamin width and depth based on the analytical wavelength.

# Stiffness ratio
SR = 100

# Elasticity parameters (Plane strain assumptions)
E, nu = 10.0, 0.45
# Lame Constants
mu_f, lmbda_f = Constant(SR*E/(2*(1 + nu))), Constant(SR*E*nu/((1 + nu)*(1 - 2*nu)))
mu_s, lmbda_s = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# define mesh based on analytical wavelength
Hf = 1.0 # film height
ana_wavelength = float(2*np.pi*Hf*((mu_f)/(3*mu_s))**(1/3))
Hs = ana_wavelength*4.0
L = ana_wavelength*3 # length of the domain
SR = 100 # stiffness ratio
d = Hf/2 # mesh size

nx = int(L/d)
ny = int((Hs+Hf)/d)
mesh = RectangleMesh(Point(0,0), Point(L,Hs+Hf), nx, ny, "right")
V = VectorFunctionSpace(mesh, "Lagrange", 2)
plt.figure(figsize=(7,7))
plot(mesh)

[<matplotlib.lines.Line2D at 0x7ff3adc0da80>,
 <matplotlib.lines.Line2D at 0x7ff3adc0dd50>]

# Define Variational Form
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u = Function(V)                 # Solution

# Mark boundary subdomians
domain_markers = MeshFunction("size_t", mesh, mesh.topology().dim())
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

def Film(x, on_boundary):
    return x[1] >= Hs
    
def Substrate(x, on_boundary):
    return x[1] < Hs
    

AutoSubDomain(Film).mark(domain_markers,0)
AutoSubDomain(Substrate).mark(domain_markers,1)

# Check Domain
plt.figure(figsize=(7,7))
plt.title('Domain')
plot(domain_markers)

<matplotlib.collections.PolyCollection at 0x7ff3ad3c1db0>

Define Dirichlet Boundary Conditions

Note that for the case of displacement control, the FEniCS expression for the applied displacement mujst have be positive to prevent convergence issues.

For example:

apply_disp = Expression("t", t = 0.0, degree = 0)

is valid and will not have convergence issues while

apply_disp = Expression("-t", t = 0.0, degree = 0)

can cause convergence issues.

The direction of applied loading will be determined by the initial load step.

# Dirichlet boundary conditions
def Left(x, on_boundary):
    return on_boundary and near(x[0], 0, 1e-6)

def Right(x, on_boundary):
    return on_boundary and near(x[0], L, 1e-6 )

def Bottom(x, on_boundary):
        return on_boundary and near(x[1], 0, 1e-6)

#-----------------Applied Displacement-----------------------#
apply_disp = Expression("t",t = 0.0, degree = 0) 
#--------------Nonhomogenous Dirichlet Boundary Conditions------------------------------#
bc1 = DirichletBC(V.sub(0), Constant(0), Left)
bc2 = DirichletBC(V.sub(0), apply_disp, Right)
bc3 = DirichletBC(V.sub(1), Constant(0), Bottom)
bcs = [bc1, bc2, bc3]

Define Function for Point Load

Note that this is an approximation of a point load since the FEniCS UserExpression will be projected into a discontinuous space. To approach a point load the area near the point load will need to have fine mesh. In our case it is not neccessary since the point load is just a perturbation.

# Function for Point 
class PointLoad(UserExpression):
    def __init__(self, x0, f, tol,**kwargs):
        super().__init__(**kwargs)
        self.x0 = x0
        self.f = f
        self.tol = tol
    def eval(self, values, x):
        if near (x[0], self.x0[0],self.tol) and near(x[1], self.x0[1],self.tol):
            values[0] = self.f[0]
            values[1] = self.f[1]
        else:
            values[0] = 0
            values[1] = 0
    def value_shape(self):
        return (2,)

Kinematics and Weak form

# Kinematics
I = variable(Identity(mesh.topology().dim()))  # Identity tensor
F = variable(I + grad(u))                        # Deformation gradient
C = variable(F.T*F)                              # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Define Variational Form
dx = Measure('dx', domain=mesh, subdomain_data=domain_markers)

# Stored strain energy density (compressible neo-Hookean model)
psi_f = (mu_f/2)*(Ic - 2) - mu_f*ln(J) + (lmbda_f/2)*(ln(J))**2
psi_s = (mu_s/2)*(Ic - 2) - mu_s*ln(J) + (lmbda_s/2)*(ln(J))**2

# Applied Traction and Body Force
pert = 8e-2
T = Constant((0,0)) # Traction
B = Constant((0,0)) # Body Force
P = PointLoad(x0 = [L/2,Hf+Hs], f = [0,pert], tol = 1e-4, degree = 0, element=V.ufl_element()) # Point load

F_int = derivative(psi_f*dx(0) + psi_s*dx(1), u, v)
F_ext = derivative(dot(B, u)*dx + dot(T,u)*ds + dot(-P,u)*ds , u, v)
residual = F_int-F_ext
J = derivative(residual, u, du)

Solver

To use our solver we first have to define the type of solver (i.e. displacement control or force control) and solver parameters before using the solver. Note that the correct type of solver has to first be imported (see first cell).

Solver parameters

Here the parameters for both types of solvers:

• psi : the scalar arc-length parameter. When psi = 1, the method becomes the spherical arc-length method and when psi = 0 the method becomes the cylindrical arc-length method

• abs_tol (optional) : absolute residual tolerance for the linear solver (default value: 1e-10)

• rel_tol (optional) : relative residual tolerance for solver; the relative residual is defined as the ration between the current residual and initial residual (default value: DOLFIN_EPS)

• lmbda0 : the initial load parameter

• max_iter : maximum number of iterations for the linear solver

• solver (optional): type of linear solver for the FEniCS linear solve function – default FEniCS linear solver is used if no argument is used.

Aside from these solver parameters, the arguments need to solve the FEA problem must also be passed into the solver:

• u : the solution function

• F_int : First variation of strain energy (internal nodal forces)

• F_ext : Externally applied load (external applied force)

• J : The Jacobian of the residual with respect to the deformation (tangential stiffness matrix)

• displacement_factor : The incremental displacement factor

The solver can be called by:

solver = force_control(psi,abs_tol,rel_tol,lmbda0,max_iter,u,F_int,F_ext,bcs,J,displacement_factor,solver)

Using the solver

1. Initialize the solver by calling solver.initialize()

2. Iteratively call solver.solve() until desired stopping condition

# Solver Parameters
psi = 1.0
tol = 1.0e-6
lmbda0 = -0.07
max_iter = 30

# Set up arc-length solver
solver = displacement_control(psi=psi, abs_tol=tol, lmbda0=lmbda0, max_iter=max_iter, u=u,
                       F_int=F_int, F_ext=F_ext, bcs=bcs, J=J, displacement_factor=apply_disp)

disp = [u.vector().copy()]
lmbda = [0]
# Function space to compute reaction force at each iteration
v_reac = Function(V)
bcRx = DirichletBC(V.sub(0), Constant(1.0), Left) # take reaction force from the left boundary
f_reac = [0.0]

strain_crit = float((1/4)*((3*mu_s/mu_f))**(2/3))# critial strain
while np.abs(apply_disp.t) < strain_crit*L*1.2 and solver.converged:
    solver.solve()
    if solver.converged:
        # Store whole displacement field
        disp.append(u.vector().copy())
        # Store displacement factor
        lmbda.append(apply_disp.t)
        # Compute and store reaction force
        bcRx.apply(v_reac.vector())
        f_reac.append(assemble(action(residual,v_reac)))

Initializing solver parameters...
Starting initial Displacement Control Control with Newton Method:
Iteration 0: | 
Absolute Residual: 1.1603e+03| Relative Residual: 1.0000e+00
Iteration 1: | 
Absolute Residual: 8.5740e+02| Relative Residual: 7.3896e-01
Iteration 2: | 
Absolute Residual: 1.7790e+02| Relative Residual: 1.5333e-01
Iteration 3: | 
Absolute Residual: 1.3741e+02| Relative Residual: 1.1843e-01
Iteration 4: | 
Absolute Residual: 1.0448e+01| Relative Residual: 9.0050e-03
Iteration 5: | 
Absolute Residual: 1.3875e-01| Relative Residual: 1.1959e-04
Iteration 6: | 
Absolute Residual: 2.1193e-04| Relative Residual: 1.8265e-07
Iteration 7: | 
Absolute Residual: 2.8398e-11| Relative Residual: 2.4475e-14

Arc-Length Step 1 :
Iteration: 0 
|Total Norm: 1.1983e-10 |Residual: 2.8398e-11 |A: 1.1642e-10| Relative Norm : 1.0000e+00

Arc-Length Step 2 :
Iteration: 0 
|Total Norm: 2.1232e+03 |Residual: 4.2808e-02 |A: 2.1232e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.1557e+00 |Residual: 8.1457e-04 |A: -3.1557e+00| Relative Norm : 1.4863e-03
Iteration: 2 
|Total Norm: 1.4395e-02 |Residual: 1.2102e-10 |A: -1.4395e-02| Relative Norm : 6.7802e-06
Iteration: 3 
|Total Norm: 2.8894e-05 |Residual: 5.3262e-12 |A: 2.8894e-05| Relative Norm : 1.3609e-08
Iteration: 4 
|Total Norm: 5.7858e-08 |Residual: 5.1867e-12 |A: -5.7858e-08| Relative Norm : 2.7251e-11

Arc-Length Step 3 :
Iteration: 0 
|Total Norm: 2.1234e+03 |Residual: 3.6197e-02 |A: 2.1234e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5900e+00 |Residual: 4.2396e-08 |A: -3.5900e+00| Relative Norm : 1.6907e-03
Iteration: 2 
|Total Norm: 7.2107e-03 |Residual: 6.5145e-12 |A: 7.2107e-03| Relative Norm : 3.3958e-06
Iteration: 3 
|Total Norm: 1.4474e-05 |Residual: 6.4492e-12 |A: -1.4474e-05| Relative Norm : 6.8165e-09
Iteration: 4 
|Total Norm: 2.9453e-08 |Residual: 6.0589e-12 |A: 2.9453e-08| Relative Norm : 1.3871e-11

Arc-Length Step 4 :
Iteration: 0 
|Total Norm: 2.1233e+03 |Residual: 3.6085e-02 |A: 2.1233e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5862e+00 |Residual: 4.4008e-08 |A: -3.5862e+00| Relative Norm : 1.6889e-03
Iteration: 2 
|Total Norm: 7.2025e-03 |Residual: 9.8994e-12 |A: 7.2025e-03| Relative Norm : 3.3921e-06
Iteration: 3 
|Total Norm: 1.4458e-05 |Residual: 9.8700e-12 |A: -1.4458e-05| Relative Norm : 6.8088e-09
Iteration: 4 
|Total Norm: 2.9337e-08 |Residual: 1.0135e-11 |A: 2.9337e-08| Relative Norm : 1.3816e-11

Arc-Length Step 5 :
Iteration: 0 
|Total Norm: 2.1232e+03 |Residual: 3.5975e-02 |A: 2.1232e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5823e+00 |Residual: 4.6964e-08 |A: -3.5823e+00| Relative Norm : 1.6872e-03
Iteration: 2 
|Total Norm: 7.1942e-03 |Residual: 1.0703e-11 |A: 7.1942e-03| Relative Norm : 3.3883e-06
Iteration: 3 
|Total Norm: 1.4440e-05 |Residual: 1.1024e-11 |A: -1.4440e-05| Relative Norm : 6.8009e-09
Iteration: 4 
|Total Norm: 2.8405e-08 |Residual: 1.0788e-11 |A: 2.8405e-08| Relative Norm : 1.3378e-11

Arc-Length Step 6 :
Iteration: 0 
|Total Norm: 2.1231e+03 |Residual: 3.5864e-02 |A: 2.1231e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5784e+00 |Residual: 5.2303e-08 |A: -3.5784e+00| Relative Norm : 1.6854e-03
Iteration: 2 
|Total Norm: 7.1859e-03 |Residual: 1.2159e-11 |A: 7.1859e-03| Relative Norm : 3.3846e-06
Iteration: 3 
|Total Norm: 1.4422e-05 |Residual: 1.1911e-11 |A: -1.4422e-05| Relative Norm : 6.7928e-09
Iteration: 4 
|Total Norm: 2.8405e-08 |Residual: 1.2259e-11 |A: 2.8405e-08| Relative Norm : 1.3379e-11

Arc-Length Step 7 :
Iteration: 0 
|Total Norm: 2.1230e+03 |Residual: 3.5753e-02 |A: 2.1230e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5746e+00 |Residual: 6.1754e-08 |A: -3.5746e+00| Relative Norm : 1.6837e-03
Iteration: 2 
|Total Norm: 7.1776e-03 |Residual: 1.7571e-11 |A: 7.1776e-03| Relative Norm : 3.3808e-06
Iteration: 3 
|Total Norm: 1.4405e-05 |Residual: 1.8844e-11 |A: -1.4405e-05| Relative Norm : 6.7852e-09
Iteration: 4 
|Total Norm: 2.8638e-08 |Residual: 1.7147e-11 |A: 2.8638e-08| Relative Norm : 1.3489e-11

Arc-Length Step 8 :
Iteration: 0 
|Total Norm: 2.1229e+03 |Residual: 3.5643e-02 |A: 2.1229e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5707e+00 |Residual: 7.8105e-08 |A: -3.5707e+00| Relative Norm : 1.6820e-03
Iteration: 2 
|Total Norm: 7.1691e-03 |Residual: 1.9816e-11 |A: 7.1691e-03| Relative Norm : 3.3770e-06
Iteration: 3 
|Total Norm: 1.4387e-05 |Residual: 2.0020e-11 |A: -1.4387e-05| Relative Norm : 6.7768e-09
Iteration: 4 
|Total Norm: 2.8405e-08 |Residual: 1.9399e-11 |A: 2.8405e-08| Relative Norm : 1.3380e-11

Arc-Length Step 9 :
Iteration: 0 
|Total Norm: 2.1228e+03 |Residual: 3.5533e-02 |A: 2.1228e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5670e+00 |Residual: 1.0585e-07 |A: -3.5670e+00| Relative Norm : 1.6803e-03
Iteration: 2 
|Total Norm: 7.1607e-03 |Residual: 2.0475e-11 |A: 7.1607e-03| Relative Norm : 3.3733e-06
Iteration: 3 
|Total Norm: 1.4369e-05 |Residual: 1.9841e-11 |A: -1.4369e-05| Relative Norm : 6.7689e-09
Iteration: 4 
|Total Norm: 2.7823e-08 |Residual: 2.0318e-11 |A: 2.7823e-08| Relative Norm : 1.3107e-11

Arc-Length Step 10 :
Iteration: 0 
|Total Norm: 2.1226e+03 |Residual: 3.5422e-02 |A: 2.1226e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5633e+00 |Residual: 1.5262e-07 |A: -3.5633e+00| Relative Norm : 1.6787e-03
Iteration: 2 
|Total Norm: 7.1522e-03 |Residual: 2.1801e-11 |A: 7.1522e-03| Relative Norm : 3.3695e-06
Iteration: 3 
|Total Norm: 1.4351e-05 |Residual: 2.2095e-11 |A: -1.4351e-05| Relative Norm : 6.7609e-09
Iteration: 4 
|Total Norm: 2.8173e-08 |Residual: 2.1516e-11 |A: 2.8173e-08| Relative Norm : 1.3272e-11

Arc-Length Step 11 :
Iteration: 0 
|Total Norm: 2.1225e+03 |Residual: 3.5312e-02 |A: 2.1225e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5599e+00 |Residual: 2.3238e-07 |A: -3.5599e+00| Relative Norm : 1.6772e-03
Iteration: 2 
|Total Norm: 7.1440e-03 |Residual: 2.2486e-11 |A: 7.1440e-03| Relative Norm : 3.3658e-06
Iteration: 3 
|Total Norm: 1.4334e-05 |Residual: 2.3483e-11 |A: -1.4334e-05| Relative Norm : 6.7534e-09
Iteration: 4 
|Total Norm: 2.7241e-08 |Residual: 2.1829e-11 |A: 2.7241e-08| Relative Norm : 1.2835e-11

Arc-Length Step 12 :
Iteration: 0 
|Total Norm: 2.1223e+03 |Residual: 3.5200e-02 |A: 2.1223e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5569e+00 |Residual: 3.7247e-07 |A: -3.5569e+00| Relative Norm : 1.6759e-03
Iteration: 2 
|Total Norm: 7.1361e-03 |Residual: 2.3503e-11 |A: 7.1361e-03| Relative Norm : 3.3624e-06
Iteration: 3 
|Total Norm: 1.4318e-05 |Residual: 2.3736e-11 |A: -1.4318e-05| Relative Norm : 6.7462e-09
Iteration: 4 
|Total Norm: 2.9570e-08 |Residual: 2.3117e-11 |A: 2.9569e-08| Relative Norm : 1.3932e-11

Arc-Length Step 13 :
Iteration: 0 
|Total Norm: 2.1222e+03 |Residual: 3.5088e-02 |A: 2.1222e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5548e+00 |Residual: 6.2993e-07 |A: -3.5548e+00| Relative Norm : 1.6751e-03
Iteration: 2 
|Total Norm: 7.1290e-03 |Residual: 2.7705e-11 |A: 7.1290e-03| Relative Norm : 3.3593e-06
Iteration: 3 
|Total Norm: 1.4301e-05 |Residual: 2.8306e-11 |A: -1.4301e-05| Relative Norm : 6.7388e-09
Iteration: 4 
|Total Norm: 2.7707e-08 |Residual: 2.9772e-11 |A: 2.7707e-08| Relative Norm : 1.3056e-11

Arc-Length Step 14 :
Iteration: 0 
|Total Norm: 2.1220e+03 |Residual: 3.4973e-02 |A: 2.1220e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5545e+00 |Residual: 1.1330e-06 |A: -3.5545e+00| Relative Norm : 1.6751e-03
Iteration: 2 
|Total Norm: 7.1236e-03 |Residual: 3.4956e-11 |A: 7.1236e-03| Relative Norm : 3.3570e-06
Iteration: 3 
|Total Norm: 1.4290e-05 |Residual: 3.6204e-11 |A: -1.4290e-05| Relative Norm : 6.7340e-09
Iteration: 4 
|Total Norm: 2.8289e-08 |Residual: 3.6479e-11 |A: 2.8289e-08| Relative Norm : 1.3331e-11

Arc-Length Step 15 :
Iteration: 0 
|Total Norm: 2.1218e+03 |Residual: 3.4855e-02 |A: 2.1218e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5579e+00 |Residual: 2.1973e-06 |A: -3.5579e+00| Relative Norm : 1.6768e-03
Iteration: 2 
|Total Norm: 7.1214e-03 |Residual: 3.7540e-11 |A: 7.1214e-03| Relative Norm : 3.3562e-06
Iteration: 3 
|Total Norm: 1.4286e-05 |Residual: 3.6748e-11 |A: -1.4286e-05| Relative Norm : 6.7329e-09
Iteration: 4 
|Total Norm: 2.9453e-08 |Residual: 3.8109e-11 |A: 2.9453e-08| Relative Norm : 1.3881e-11

Arc-Length Step 16 :
Iteration: 0 
|Total Norm: 2.1216e+03 |Residual: 3.4728e-02 |A: 2.1216e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.5699e+00 |Residual: 4.6974e-06 |A: -3.5699e+00| Relative Norm : 1.6826e-03
Iteration: 2 
|Total Norm: 7.1233e-03 |Residual: 3.9741e-11 |A: 7.1233e-03| Relative Norm : 3.3575e-06
Iteration: 3 
|Total Norm: 1.4286e-05 |Residual: 3.8248e-11 |A: -1.4286e-05| Relative Norm : 6.7334e-09
Iteration: 4 
|Total Norm: 2.9337e-08 |Residual: 3.8661e-11 |A: 2.9337e-08| Relative Norm : 1.3827e-11

Arc-Length Step 17 :
Iteration: 0 
|Total Norm: 2.1214e+03 |Residual: 3.4584e-02 |A: 2.1214e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.6030e+00 |Residual: 1.1468e-05 |A: -3.6030e+00| Relative Norm : 1.6984e-03
Iteration: 2 
|Total Norm: 7.1208e-03 |Residual: 5.4688e-11 |A: 7.1208e-03| Relative Norm : 3.3566e-06
Iteration: 3 
|Total Norm: 1.4281e-05 |Residual: 4.0893e-11 |A: -1.4281e-05| Relative Norm : 6.7319e-09
Iteration: 4 
|Total Norm: 2.9570e-08 |Residual: 4.1504e-11 |A: 2.9569e-08| Relative Norm : 1.3939e-11

Arc-Length Step 18 :
Iteration: 0 
|Total Norm: 2.1211e+03 |Residual: 3.4398e-02 |A: 2.1211e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.6962e+00 |Residual: 3.4000e-05 |A: -3.6962e+00| Relative Norm : 1.7426e-03
Iteration: 2 
|Total Norm: 7.0135e-03 |Residual: 1.9105e-10 |A: 7.0135e-03| Relative Norm : 3.3065e-06
Iteration: 3 
|Total Norm: 1.4066e-05 |Residual: 3.8800e-11 |A: -1.4066e-05| Relative Norm : 6.6314e-09
Iteration: 4 
|Total Norm: 2.9220e-08 |Residual: 3.9831e-11 |A: 2.9220e-08| Relative Norm : 1.3776e-11

Arc-Length Step 19 :
Iteration: 0 
|Total Norm: 2.1208e+03 |Residual: 3.4092e-02 |A: 2.1208e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 3.9966e+00 |Residual: 1.3823e-04 |A: -3.9966e+00| Relative Norm : 1.8845e-03
Iteration: 2 
|Total Norm: 5.7595e-03 |Residual: 1.5407e-09 |A: 5.7595e-03| Relative Norm : 2.7157e-06
Iteration: 3 
|Total Norm: 1.1572e-05 |Residual: 4.2656e-11 |A: -1.1572e-05| Relative Norm : 5.4567e-09
Iteration: 4 
|Total Norm: 2.2817e-08 |Residual: 4.1945e-11 |A: 2.2817e-08| Relative Norm : 1.0759e-11

Arc-Length Step 20 :
Iteration: 0 
|Total Norm: 2.1201e+03 |Residual: 3.3338e-02 |A: 2.1201e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 5.3120e+00 |Residual: 1.0296e-03 |A: -5.3120e+00| Relative Norm : 2.5056e-03
Iteration: 2 
|Total Norm: 1.3499e-02 |Residual: 3.3193e-08 |A: -1.3499e-02| Relative Norm : 6.3673e-06
Iteration: 3 
|Total Norm: 2.7670e-05 |Residual: 4.3386e-11 |A: 2.7670e-05| Relative Norm : 1.3051e-08
Iteration: 4 
|Total Norm: 5.7160e-08 |Residual: 4.4566e-11 |A: -5.7160e-08| Relative Norm : 2.6961e-11

Arc-Length Step 21 :
Iteration: 0 
|Total Norm: 2.1180e+03 |Residual: 2.9969e-02 |A: 2.1180e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 1.8493e+01 |Residual: 3.9738e-02 |A: -1.8493e+01| Relative Norm : 8.7311e-03
Iteration: 2 
|Total Norm: 1.1711e+00 |Residual: 1.3147e-05 |A: -1.1711e+00| Relative Norm : 5.5291e-04
Iteration: 3 
|Total Norm: 8.3386e-03 |Residual: 4.9781e-09 |A: 8.3386e-03| Relative Norm : 3.9369e-06
Iteration: 4 
|Total Norm: 1.6346e-05 |Residual: 4.2348e-11 |A: -1.6346e-05| Relative Norm : 7.7173e-09
Iteration: 5 
|Total Norm: 3.2480e-08 |Residual: 4.2963e-11 |A: 3.2480e-08| Relative Norm : 1.5335e-11

Arc-Length Step 22 :
Iteration: 0 
|Total Norm: 2.1017e+03 |Residual: 1.0620e-01 |A: 2.1017e+03| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 1.3920e+02 |Residual: 1.4753e+00 |A: 1.3920e+02| Relative Norm : 6.6234e-02
Iteration: 2 
|Total Norm: 8.2069e+01 |Residual: 2.1121e+00 |A: -8.2042e+01| Relative Norm : 3.9049e-02
Iteration: 3 
|Total Norm: 5.0584e+02 |Residual: 1.2935e+01 |A: -5.0567e+02| Relative Norm : 2.4068e-01
Iteration: 4 
|Total Norm: 5.7802e+02 |Residual: 2.0530e-02 |A: -5.7802e+02| Relative Norm : 2.7503e-01
Iteration: 5 
|Total Norm: 3.2218e+02 |Residual: 6.1237e-01 |A: -3.2218e+02| Relative Norm : 1.5329e-01
Iteration: 6 
|Total Norm: 3.3445e+01 |Residual: 1.4510e-02 |A: 3.3445e+01| Relative Norm : 1.5913e-02
Iteration: 7 
|Total Norm: 1.0291e+01 |Residual: 5.5474e-04 |A: 1.0291e+01| Relative Norm : 4.8964e-03
Iteration: 8 
|Total Norm: 9.2142e-02 |Residual: 8.4985e-08 |A: 9.2142e-02| Relative Norm : 4.3842e-05
Iteration: 9 
|Total Norm: 4.8233e-05 |Residual: 4.5119e-11 |A: 4.8233e-05| Relative Norm : 2.2949e-08
Iteration: 10 
|Total Norm: 1.0012e-08 |Residual: 4.6327e-11 |A: 1.0012e-08| Relative Norm : 4.7637e-12

Arc-Length Step 23 :
Iteration: 0 
|Total Norm: 8.7006e+02 |Residual: 4.7428e+01 |A: 8.6876e+02| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 1.7038e+03 |Residual: 1.3156e-01 |A: -1.7038e+03| Relative Norm : 1.9582e+00
Iteration: 2 
|Total Norm: 5.8139e+02 |Residual: 5.3843e-01 |A: 5.8139e+02| Relative Norm : 6.6823e-01
Iteration: 3 
|Total Norm: 3.9780e+01 |Residual: 6.6247e-03 |A: 3.9780e+01| Relative Norm : 4.5721e-02
Iteration: 4 
|Total Norm: 3.5177e+00 |Residual: 2.7393e-05 |A: 3.5177e+00| Relative Norm : 4.0430e-03
Iteration: 5 
|Total Norm: 5.1063e-03 |Residual: 3.1942e-10 |A: 5.1063e-03| Relative Norm : 5.8689e-06
Iteration: 6 
|Total Norm: 9.3877e-07 |Residual: 5.0740e-11 |A: 9.3877e-07| Relative Norm : 1.0790e-09

Arc-Length Step 24 :
Iteration: 0 
|Total Norm: 1.5132e+02 |Residual: 2.9476e+01 |A: -1.4842e+02| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 6.9015e+02 |Residual: 1.9051e-01 |A: -6.9015e+02| Relative Norm : 4.5610e+00
Iteration: 2 
|Total Norm: 6.3485e+02 |Residual: 4.2235e-01 |A: 6.3485e+02| Relative Norm : 4.1955e+00
Iteration: 3 
|Total Norm: 3.9017e+01 |Residual: 3.4629e-03 |A: 3.9017e+01| Relative Norm : 2.5785e-01
Iteration: 4 
|Total Norm: 1.5919e+00 |Residual: 3.6610e-06 |A: 1.5919e+00| Relative Norm : 1.0520e-02
Iteration: 5 
|Total Norm: 8.6827e-04 |Residual: 5.9728e-11 |A: 8.6827e-04| Relative Norm : 5.7380e-06
Iteration: 6 
|Total Norm: 1.2259e-07 |Residual: 5.3523e-11 |A: 1.2259e-07| Relative Norm : 8.1012e-10

Arc-Length Step 25 :
Iteration: 0 
|Total Norm: 1.3318e+02 |Residual: 1.4316e+01 |A: -1.3241e+02| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 2.0608e+02 |Residual: 9.3473e-02 |A: -2.0608e+02| Relative Norm : 1.5473e+00
Iteration: 2 
|Total Norm: 1.9715e+02 |Residual: 3.2933e-02 |A: 1.9715e+02| Relative Norm : 1.4803e+00
Iteration: 3 
|Total Norm: 3.2879e+00 |Residual: 1.9443e-05 |A: 3.2879e+00| Relative Norm : 2.4687e-02
Iteration: 4 
|Total Norm: 8.3677e-03 |Residual: 2.0951e-10 |A: 8.3677e-03| Relative Norm : 6.2828e-05
Iteration: 5 
|Total Norm: 1.0462e-06 |Residual: 6.0425e-11 |A: 1.0462e-06| Relative Norm : 7.8555e-09
Iteration: 6 
|Total Norm: 9.3324e-10 |Residual: 5.9786e-11 |A: 9.3132e-10| Relative Norm : 7.0072e-12

Arc-Length Step 26 :
Iteration: 0 
|Total Norm: 1.1688e+02 |Residual: 9.4880e+00 |A: -1.1650e+02| Relative Norm : 1.0000e+00
Iteration: 1 
|Total Norm: 8.9529e+01 |Residual: 4.8316e-02 |A: -8.9529e+01| Relative Norm : 7.6598e-01
Iteration: 2 
|Total Norm: 8.4197e+01 |Residual: 4.8963e-03 |A: 8.4197e+01| Relative Norm : 7.2036e-01
Iteration: 3 
|Total Norm: 5.0529e-01 |Residual: 4.6408e-07 |A: 5.0529e-01| Relative Norm : 4.3230e-03
Iteration: 4 
|Total Norm: 2.1869e-04 |Residual: 6.2127e-11 |A: 2.1869e-04| Relative Norm : 1.8710e-06
Iteration: 5 
|Total Norm: 2.2468e-08 |Residual: 6.2099e-11 |A: 2.2468e-08| Relative Norm : 1.9223e-10

Optional: Stop at at specific criterion

To strongly impose the stop condition, we add Newton solver at the end with a specific prescribed displacement and replace the final solution with the solution from the Newton Solver. In this case, we enforced the stop criterion with the applied displacement is slightly over the analytical critical strain.

apply_disp.t = -strain_crit*L*1.2 # stop simulation after bifurcation
problem = NonlinearVariationalProblem(residual, u, bcs, J)
newton_solver = NonlinearVariationalSolver(problem)
newton_solver.solve()

disp.pop()
lmbda.pop()
f_reac.pop()
disp.append(u.vector())
lmbda.append(apply_disp.t)
f_reac.append(assemble(action(residual,v_reac)))

Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 2.151e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 2.765e-01 (tol = 1.000e-10) r (rel) = 1.285e-03 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 7.750e-05 (tol = 1.000e-10) r (rel) = 3.602e-07 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 8.076e-09 (tol = 1.000e-10) r (rel) = 3.754e-11 (tol = 1.000e-09)
  Newton solver finished in 3 iterations and 3 linear solver iterations.

Post Processing

Here we plot the final deformed shape, equilibrium path, and the wavelength. We obtain the analytical solutions from: https://royalsocietypublishing.org/doi/epdf/10.1098/rsta.2016.0163 and https://groups.seas.harvard.edu/hutchinson/papers/WrinklingPhenonmena-JAMF.pdf

plt.figure(figsize=(7,7))
plot(u, mode = 'displacement', cmap='Reds', edgecolor = 'k', linewidth = 0.2)
plt.title('Final Deformation')

Text(0.5, 1.0, 'Final Deformation')

test_y = np.diff(np.array(f_reac))

test_x = np.diff(-np.array(lmbda)/L)

fea_soln = np.nanargmax(np.abs((np.diff(np.diff(((np.array(f_reac[1:])/(Hs+Hf))/(-np.array(lmbda[1:]))/L)))))) # find the inflection point from FEA solution and use that as critical strain

plt.figure(figsize=(7,5))
plt.plot(-np.array(lmbda)/L, np.array(f_reac)/(Hs+Hf), c='k', marker = 'o', label = 'Equilibrium Path')
plt.plot(-np.array(lmbda[fea_soln+2])/L, np.array(f_reac[fea_soln+2])/(Hs+Hf), marker = 'd', c =(0.4,0,0) , ls = 'None', markersize = 10, label='FEA critical strain')
plt.axvline(strain_crit, ls = ':', c = (0.7,0,0), label = 'Analytical Critical Strain')
plt.xlabel('Applied Strain')
plt.ylabel('Stress per unit depth')
plt.title('Equilibrium path')
plt.legend()

percent_diff_crit_strain = ((strain_crit-(-np.array(lmbda[fea_soln+2])/L))/strain_crit) * 100
print('Percent Error between analytical critical strain and FEA critical strain:',percent_diff_crit_strain,'%')

Percent Error between analytical critical strain and FEA critical strain: 1.5933169703227463 %

Here we extract and compare the post-bifurcation wrinkling wavelength from FEA and analytical solution.

post_bif = np.argwhere(-np.array(lmbda)/L > strain_crit).reshape(-1) # get index after bifurcation
# get solution after onset of bifurcation:
disp_bif = disp[post_bif[0]]
u_bif = Function(V) 
u_bif.vector()[:] = disp_bif[:]


# Get dof coordinates:
x_dofs = V.sub(0).dofmap().dofs()
y_dofs = V.sub(1).dofmap().dofs()
theta_dofs = V.sub(1).dofmap().dofs()
dofs = V.tabulate_dof_coordinates()
dof_coords = dofs.reshape((-1, 2))

x_nodal_coord = dof_coords[x_dofs][:,0]
y_nodal_coord = dof_coords[y_dofs][:,1]
# Get nodal values 

top_layer = np.argwhere(np.abs(y_nodal_coord-(Hs+Hf)) < 1e-6)

# Plot displacement field
disp_x = x_nodal_coord[top_layer].reshape(-1)#(x_nodal_coord + u_bif.vector()[x_dofs])[top_layer].reshape(-1)
disp_y = (y_nodal_coord + u_bif.vector()[y_dofs])[top_layer].reshape(-1)

# Fit to cosine functions
from scipy import optimize

def fit_func(x, amp, wave, phi, offset):
    return amp * np.cos((2*np.pi)/wave * x + phi) + offset

params, params_covariance = optimize.curve_fit(fit_func, disp_x, disp_y,
                                               p0=[1.0,ana_wavelength,0.0,Hf+Hs])

plt.figure(figsize = (7,7))

plt.title('FEA wavelength vs Analytical Wavelength')
plt.scatter(disp_x,disp_y, marker = '.', c = 'r', s = 50, label='FEA data')
plt.scatter(disp_x,fit_func(disp_x, params[0], params[1], params[2], params[3]), c = 'None', marker = 's',edgecolor = 'k', s = 50, linewidth = 2, label = 'Fitted FEA data')
plt.scatter(disp_x,fit_func(disp_x, params[0], ana_wavelength, params[2], params[3]), c = 'None', marker = 'd', edgecolor = (0.5,0,0), s = 50, linewidth = 2, label = 'Analytical Wavelength')
plt.xlabel('x displacement')
plt.ylabel('y displacement')
plt.legend(loc = (1.01,0.5))

print('Percent Error between analytical wavelength and fitted FEM wavelength:',((params[1]-ana_wavelength)/ana_wavelength)*100, '%')

Percent Error between analytical wavelength and fitted FEM wavelength: 0.26494519096416136 %

Optional: Creating an animation from solution snapshots

For aesthethics reason, we only look at the thin fim part of the bilayer.

from matplotlib import animation, rc

plt.rcParams["animation.html"] = "jshtml"

u_plot = Function(V)

fig = plt.figure(figsize=(12,5))
ax = fig.gca()

# Make sure to only look at surface
ax.set_xlim([0,L])
ax.set_ylim([Hs,Hs+Hf+3])

def drawframe(n):
    fig.clf()
    ax = fig.gca()
    ax.set_xlim([0,L])
    ax.set_ylim([Hs,Hf+Hs+5])
    u_plot = Function(V)
    u_plot.vector()[:] = disp[n][:]
    return plot(u_plot, mode = 'displacement', cmap = 'Reds'),
# blit=True re-draws only the parts that have changed.
anim = animation.FuncAnimation(fig, drawframe, frames=len(lmbda), interval=40, blit=True)
anim

Once Loop Reflect