Hi everyone,
My codes runs but I have get the wrong solution. I just follow the link as provided.
https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html
It runs but I think at my problem there should be a wave. It is about a 3d transient elasto-dynamics problem with BC at Y 0.05.
Here is my code.
Thanks.
#include FEniCS headers
from __future__ import print_function
from fenics import *
from ufl import nabla_div, nabla_grad, div
import numpy as np
# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
#Problem parameters variables
L = 1 #Beam length
W = 0.1 #Beam width
H = 0.1 #Beam height
E = 1000 #Young's modulus
nu= 0.3 #Poisson's ratio
lambda_ = (E*nu)/((1+nu)*(1-2*nu)) #Lame parameters
mu = E/(2*(1+nu)) #Lame parameters
rho=1.0 #density
a=1.0 # Rayleigh damping coefficients
b=0.001 # Rayleigh damping coefficients
gamma=1/2 # Newmark method parameters
beta=1/4 # Newmark method parameters
#Create mesh
#The points represent the diagonal ends of the box, followed by the number of elments along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(W, H, L ), 8, 8, 80) #3D mesh
# Define function space for displacement, velocity and acceleration
V = VectorFunctionSpace(mesh, "CG", 1)
# Define boundary condition
# Define function to identify the left boundary that is clamped
def clamped_boundary(x, on_boundary):
return on_boundary and x[2] < DOLFIN_EPS
# Create the boundary condition
bc_clamped = DirichletBC(V, Constant((0.00,0.00,0.00)), clamped_boundary) #3D
# Define function to identify the right boundary
def on_right(x, on_boundary):
return on_boundary and x[2] > 1 - DOLFIN_EPS
#Create the boundary condition
bc_right = DirichletBC(V.sub(1), Constant(0.05), on_right) #3D
# Combine the boundary condition
bc = [ bc_clamped, bc_right]
# Time-stepping parameters
T = 0.05 #Total T
Nsteps = 100 #Number of steps
dt = Constant(T/Nsteps)#delta t
# Creat output files
vtkfile = File('newmark/newmark.pvd')
# Test and trial functions
du = TrialFunction(V) #trial function
w = TestFunction(V) #test function
# Current (unknown) displacement
u = Function(V, name="Displacement")
# Fields from previous time step (displacement, velocity, acceleration)
u_old = Function(V) #the pervious displacement
v_old = Function(V) #the pervious velocity
a_old = Function(V) #the pervious acceleration
d = u.geometric_dimension() # space dimension
# Define expressions for strain and stress function
#strain
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#stress (St.Venant Kirchhoff model)
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# Mass form
def m(u, w):
return rho*inner(u, w)*dx
# Elastic stiffness form
def k(u, w):
return inner(sigma(u), epsilon(w))*dx
# Rayleigh damping form
def c(u, w):
return a*m(u, w) + b*k(u, w)
# Update formula for acceleration
# a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
beta_ = beta
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
# Update formula for velocity
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
def update_v(a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)
# Update formula for displacement
def update_u(u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return u_old + dt_*v_old+dt_*dt_/2*(1-2*beta)*a_old
# Define initial condition(I.C)(displacement)Previous
u_0=Expression(('0', '0', '0'), degree=1)
# Interpolate it
u_old=interpolate(u_0,V)
# Define initial condition(I.C)(velocity)Previous
v_0=Expression(('0', '0', '0'), degree=1)
# Interpolate it
v_old=interpolate(v_0,V)
# Define the weak form
a_new = update_a(du, u_old, v_old, a_old, ufl=True) #define the updated a
v_new = update_v(a_new, u_old, v_old, a_old, ufl=True) #define the updated v
u_new = update_u(u_old, v_old, a_old, ufl=True)
res = m(a_new, w) + c(v_new,w) + k(du, w) #here is the weak form(KMC equation)
a_form = lhs(res) #left hand side
L_form = rhs(res) #right hand side
# Compute solution by solving the weak formulation
u = Function(V)
t=0 #initial t
# Do time stepping
for n in range(Nsteps):
# Update current time
t += dt
# Compute solution
solve(a_form== L_form, u, bc) #solve the weak form
# Save to file and plot solution
vtkfile << (u, t)
# Update previous solution
u_old.assign(u) #update new u