Elasticity problem

1

Hello everyone,
I’m trying to run the tutorials: elasticity problem
the code is following:
from dolfin import *

#Scaled variables
L=1; W= 0.2
mu=1
rho=1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_= beta
g = gamma

#Create mesh and define function space
mesh = BoxMesh(Point(0,0,0),Point(L,W,W),10,3,3)
V = VectorFunctionSpace(mesh,‘P’,1)

#Define boundary condition
tol = 1E-14

def clamped_boundary(x,on_boundary):
return on_boundary and x[0]<tol

bc = DirichletBC(V, Constant((0,0,0)),clamped_boundary)

#Define strain and stress

def epsilon(u):
return 0.5*(nabla_grad(u)+nabla_grad(u).T)

def sigma(u):
return lambda_nabla_grad(u)Identity(d)+2muepsilon(u)

#Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() #space dimension
v = TestFunction(V)
f = Constant((0,0,-rho*g))
T = Constant((0,0,0))
a = inner(sigma(u),epsilon(v))*dx
L = dot(f,v)*dx+dot(T,v)*ds

#Compute solution
u = Function(V)
solve(a==L,u,bc)

#plot solution
plot(u,title = ‘Displacement’, mode = ‘displacement’)

#Plot stress
s = sigma(u)-(1./3)*tr(sigma(u))Identity(d) #deviatoric stress
von_Mises = sqrt(3./2inner(s,s))
V = FunctionSpace(mesh,‘P’,1)
von_Mises = project(von_Mises,V)
plot(von_Mises,title=‘Stress intensity’)

#Compute magutude of displacement
u_magnitude = sqrt(dot(u,u))
u_magnitude = project(u_magnitude,V)
plot(u_magnitude, ‘Displacement magnitude’)

print(‘min/max u:’,
u_magnitude.vector().array().min(),
u_magnitude.vector().array().max())

The above doesn’t work. Can somebody help me ?
Thanks in advance!

2

Hi
It seems like you have missed multiplication symbol (*) in defining sigma(u) and von_Mises stress. This code should work fine:

from fenics import *

# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u, title='Displacement', mode='displacement')

s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  
von_Mises = sqrt(3./2*inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')

u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
      u_magnitude.vector().array().min(),
      u_magnitude.vector().array().max())
3

Thank you very much that really helped me a lot.