I would like to use the Crank-Nicolson method to simulate the dynamics of a linear elastic system. As an example problem I set up a simple cantilever structure. In a first step I obtain the deflection for a static point load applied to the free end of the cantilver. The solution to the static part should then serve as the initial condition for the dynamic problem which is the mechanical oscillation from this deflection in the absence of any load.
To obtain the time-dependent solution I rewrite the initial PDE \nabla \cdot \pmb{\sigma} = \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} as a system of two coupled equations ( \mathbf{u}_1 = \mathbf{u}, \mathbf{u}_2 = \dot{\mathbf{u}} ) and derive the corresponding linear forms.
However, I am not sure if I implemented the forms as I am supossed to. To solve this problem I introduce an mixed function space with
V1 = fe.VectorElement('Lagrange', mesh.ufl_cell(), 3)
V2 = fe.VectorElement('Lagrange', mesh.ufl_cell(), 3)
Vt = fe.FunctionSpace(mesh, V1 * V2)
Just executing this command takes a couple of seconds. Even longer takes the assembly of the system using
A, b = fe.assemble_system(a, L, [bc1, bc2])
The final step solution step does not converge and I obtain after a long computation the error message
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
You find the complete code of my example below. This code only computes the first time-step of the time-dependent solution.
Is my problem related to an unproper definition of the function spaces or am my forms not correctly defined? I would be very happy if someone could point me in the right direction. Thanks in advance.
import fenics as fe
import numpy as np
import matplotlib.pyplot as plt
# Geometry parameters
l = 100e-6
h = 0.1 * l
w = 0.1 * l
# Mesh parameters
nl = 50
nw = 5
nh = 5
# Material parameters
E = 130e9
nu = 0.22
rho = 2330 # density kg/m**3
# Lame parameters
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
# Load parameters
f_ = fe.Constant((0, 0, 0))
delta = fe.PointSource(V.sub(2), fe.Point(l, w / 2, h / 2), -1)
mesh = fe.BoxMesh(fe.Point(0, 0, 0), fe.Point(l, w, h), nl, nw, nh)
V = fe.VectorFunctionSpace(mesh, 'Lagrange', 1)
def boundary(x, on_boundary):
return on_boundary and fe.near(x[0], 0)
bc = fe.DirichletBC(V, fe.Constant((0, 0, 0)), boundary)
def epsilon(u):
return 0.5 * (fe.grad(u) + fe.grad(u).T)
def sigma(u):
return 2 * mu * epsilon(u) + lmbda * (fe.div(u)) * fe.Identity(3)
u = fe.TrialFunction(V)
v = fe.TestFunction(V)
a = fe.inner(sigma(u), epsilon(v)) * fe.dx
L = fe.inner(f_, v) * fe.dx
A, b = fe.assemble_system(a, L, bc)
delta.apply(b)
uh = fe.Function(V)
fe.solve(A, uh.vector(), b)
uh.rename('displacement', 'displacement')
f = fe.File('static.pvd')
f << uh
# ## Time-Dependent Problem
dt = 0.0001
V1 = fe.VectorElement('Lagrange', mesh.ufl_cell(), 3)
V2 = fe.VectorElement('Lagrange', mesh.ufl_cell(), 3)
Vt = fe.FunctionSpace(mesh, V1 * V2)
def boundary(x, on_boundary):
return on_boundary and fe.near(x[0], 0)
bc1 = fe.DirichletBC(Vt.sub(0), fe.Constant((0, 0, 0)), boundary)
bc2 = fe.DirichletBC(Vt.sub(1), fe.Constant((0, 0, 0)), boundary)
v1, v2 = fe.TestFunctions(Vt)
u1, u2 = fe.TrialFunctions(Vt)
u10, u20 = fe.Function(Vt).split()
u10.assign(uh)
a1 = fe.inner(u1, v1) * fe.dx - fe.Constant(0.5 * dt) * fe.inner(u2, v1) * fe.dx
a2 = fe.inner(u2, v2) * fe.dx - fe.Constant(0.5 * dt / rho) * fe.inner(sigma(u1), epsilon(v2)) * fe.dx
a = a1 + a2
L1 = fe.inner(u10, v1) * fe.dx + fe.Constant(0.5 * dt) * fe.inner(u20, v1) * fe.dx
L2 = fe.inner(u20, v2) * fe.dx + fe.Constant(0.5 * dt / rho) * fe.inner(sigma(u10), epsilon(v2)) * fe.dx
L = L1 + L2
uht = fe.Function(Vt)
fe.solve(a==L, uht, [bc1, bc2])
A, b = fe.assemble_system(a, L, [bc1, bc2])