Dear joha2,
I also found this example after my last post here and did adopt it with some parts from other examples from over here (
https://github.com/FEniCS/dolfin/blob/m ... sticity.py for pvd export, ...).
In the provided example is 'gravity' induced force on the complete volume used and I tried to change this to the calculix-example used in FC.
But I struggled in the end with introducing a point load at the end of the beam.
For the 'gravity' or 'self weight' problem do match the results from calculix-FC and the Fenics script well (fine mesh with 'scalefactor = 1'), even if the solving time is much much higher for fenics.
Displacement from Fenics:
Code: Select all
Solving linear variational problem.
min/max u: -3.97926569654e-09 0.00207006247747
And FC gives:
2.24684 mm
Since I do not have time yet to further look into it, I post the last version of my file hoping it will help.
BR,
HoWil
Code: Select all
"""
FEniCS tutorial demo program: Linear elastic problem.
-div(sigma(u)) = f
The model is used to simulate an elastic beam clamped at
its left end and deformed under its own weight.
"""
from __future__ import print_function
from fenics import *
# Scaled variables
Len = 8; W = 1; H = 1 # dimensions in m
#mu = 1
rho = 7900 # kg/m^3 has to be adopted
# Elasticity parameters
E = 210.0e9 # Pa
nu = 0.3 # dimensionless
mu = E/(2.0*(1.0 + nu))
lambda_ = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
#delta = W/Len
#gamma = 0.4*delta**2
#beta = 1.25
#lambda_ = beta
g = 9.81 # m/s^2 gamma
# Create mesh and define function space
scalefactor = 1
mesh = BoxMesh(Point(0, 0, 0), Point(Len, W, H), 30/scalefactor, 15/scalefactor, 15/scalefactor)
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)
#return sym(nabla_grad(u))
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(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, -10000000)) # force vector # does not work yet!!!!!
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds(0) # self load + load on face
Fz = 100000.
#ps = PointSource(V.sub(0),Point(Len, 0., 0.), Fz) # !!!!!!! does not work
#delta = PointSource(V, Point(0., 0., 0.), 0.)
#delta.apply(b)
# Compute solution
u = Function(V)
#ps.apply(u.vector()) # !!!!!!! does not work
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./2*inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')
# Compute magnitude 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())
# Save solution to file in VTK format
File('elasticity/displacement.pvd') << u
File('elasticity/von_mises.pvd') << von_Mises
File('elasticity/magnitude.pvd') << u_magnitude
# Hold plot
interactive()
