Elastodynamic problem

dolfinx
1

Hi, I’m trying to solve a simple elastodynamic problem with a 2d beam, that’s my code:


import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
import time


# --------------------
# Functions and classes
# --------------------
def left(x, on_boundary):
    return (on_boundary and fe.near(x[0], 0.0))


# Strain function
def epsilon(u):
    return 0.5*(fe.nabla_grad(u) + fe.nabla_grad(u).T)


# Stress function
def sigma(u):
    return lmbda*fe.div(u)*fe.Identity(2) + 2*mu*epsilon(u)


# --------------------
# Parameters
# --------------------
# Young's modulus and Poisson's ratio
E = 0.02e9
nu = 0.1

# Lame's constants
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
rho = 200.0

# Time-stepping
t_start = 0.0  # start time
t_end = 1.0e-1  # end time
t_steps = 100  # number of time steps

t, dt = np.linspace(t_start, t_end, t_steps, retstep=True)
dt = float(dt) 

# --------------------
# Geometry
# --------------------
mesh = fe.RectangleMesh(fe.Point(0.,0.),fe.Point(5.,1.),100,10)

fe.plot(mesh)
plt.show()

# --------------------
# Function spaces
# --------------------
V = fe.VectorFunctionSpace(mesh, "CG", 1)
u_tr = fe.TrialFunction(V)
u_test = fe.TestFunction(V)

# --------------------
# Boundary conditions
# --------------------

bc1 = fe.DirichletBC(V, fe.Constant((0.0, 0.0)), left)
bc = [bc1]

# --------------------
# Initialization
# --------------------
u = fe.Function(V)
u_bar = fe.Function(V)
du = fe.Function(V)
ddu = fe.Function(V)
ddu_old = fe.Function(V)

file = fe.XDMFFile("FEniCS_explicit_output.xdmf")  # XDMF file

# --------------------
# Weak form
# --------------------
A_form = fe.inner(sigma(u_tr), epsilon(u_test))*fe.dx + 4*rho/(dt*dt)*fe.dot(u_tr - u_bar, u_test)*fe.dx

# --------------------
# Time loop
# --------------------
for ti in t:
    u_D.t = ti
    u_bar.assign(u + dt*du + 0.25*dt*dt*ddu)

    fe.solve(fe.lhs(A_form) == fe.rhs(A_form), u, bc)

    ddu_old.assign(ddu)
    ddu.assign(4/(dt*dt)*(u - u_bar))
    du.assign(du + 0.5*dt*(ddu + ddu_old))

    file.write(u, ti)

    #print(ti)
file.close()

The code works without errors, however when I plot the displacement I get a zero value. I have just write


fe.plot(u,mode='displacement')

I don’t understand if there is something wrong in the code or I have ton change the plot command.
Could you help me?
Thank you

2

Not sure if I am reading this right, but at the first time ti it seems to me that u_bar is assigned (a complicated expression which result is still) zero. In that case, accounting also for homogeneous Dirichlet BC, the system does not have any forcing, hence the solution u at the first time ti is still zero. Once the first u is zero, ddu and du will still be zero, and the next time steps will still have zero forcing term, hence zero solution.

3

Thank you very much. Initially indeed there was an imposed displacement that I have removed, I’ll try to fix it.