Hello,
Fenics Beginner here. I am trying to simulate the deflection of a cantilever beam with orthotropic elasticity coefficients. This code (based on Jeremy Bleyers tutorial orthotropic linear elasticity.
The code gives accurate results when I use a second order Lagrange element, but the results are off by a factor of 1000 when I use first order Lagrange elements. Even if I change the size of the grid it still does not work, so mesh size does not seem to be the issue.
Please help…this is part of a bigger simulation problem and I would really like to do make it work with first order elements.
# Cantilever beam with Orthotropic Elasticity
#Non Scaled Version
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
C11, C22, C12, C13, C33, C44 =166.5,166.5,64.1, 64.1,166.5,79.6
C = ufl.as_matrix([[C11, C12, C13, 0.0, 0.0, 0.0],
[C12, C22, C13, 0.0, 0.0, 0.0],
[C13, C13, C33, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, C44, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, C44, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, C44]])
scale =1E-06
rho=2230
escale=1E9
L = 6000.0*scale
W = 1000.0*scale
H=45*scale
delta = W / L
g=9.8
deflection = (1.5*rho*g*(L**4))/(130E9*(H**2))#Theoretical value
print (deflection)
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],[10, 6, 4], cell_type=mesh.CellType.tetrahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
def clamped_boundary(x):
return np.isclose(x[0], 0.0)
def bottom_boundary(x):
return np.isclose(x[2],0.0)
def other_boundary(x):
return np.logical_or.reduce((np.isclose(x[0],L),
np.isclose(x[1],0.0),
np.isclose(x[1],W),
np.isclose(x[2],H)))
def strain(u,repr="vectorial"):
eps_t= ufl.sym(ufl.grad(u))
if repr == "vectorial":
return ufl.as_vector([eps_t[0,0],eps_t[1,1],eps_t[2,2],2*eps_t[0,1], 2*eps_t[1,2],2*eps_t[0,2]])
elif repr == "tensorial":
return eps_t
def stress(u,repr="vectorial"):
sigv=ufl.dot(C,strain(u))
if repr == "vectorial":
return sigv
elif repr == "tensorial":
return ufl.as_matrix([[sigv[0],sigv[4],sigv[6]],
[sigv[4],sigv[1],sigv[5]],
[sigv[6],sigv[5],sigv[2]]])
fdim = domain.topology.dim - 1
fixedboundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom_boundary)
other_facets = mesh.locate_entities_boundary(domain, fdim, other_boundary)
bvalues =np.full_like(bottom_facets,1)
ovalues = np.full_like(other_facets,2)
bfacets = np.hstack([bottom_facets,other_facets])
bvalues = np.hstack ([bvalues,ovalues])
sort_order = np.argsort(bfacets)
mtb = mesh.meshtags(domain, fdim, bfacets[sort_order],bvalues[sort_order])
ds=ufl.Measure("ds", domain=domain, subdomain_data=mtb)
T = fem.Constant(domain,0.0)
n=ufl.FacetNormal(domain)
#define problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g/escale)))
a = ufl.inner(stress(u), strain(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T*n, v) * ds(1)+ufl.dot(fem.Constant(domain,0.0)*n, v) * ds((2,3,4,5))
#Define BC
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, fixedboundary_facets), V)
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
'''