I am trying to solve simple linear elasticity problem like here. I am using 1D beam element in 3D space.
from dolfin import *
import numpy as np
import meshio
length = 10.0
N = 40 # number of elements
points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
cells = np.array([[i, i+1] for i in range(N)])
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
def left_end(x, on_boundary):
return near(x[0], 0,0.1)
def right_end(x, on_boundary):
return near(x[0], 10, 0.1)
facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(left_end).mark(facets, 1)
AutoSubDomain(right_end).mark(facets, 2)
ds = Measure("ds", domain=mesh, subdomain_data=facets)
def eps(v):
return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def sigma(v):
return lmbda*tr(eps(v))*Identity(3) + 2.0*mu*eps(v)
rho_g = 1
f = Constant((0,0,-rho_g))
T = Constant((0,0,0))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=3, dim=3)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx +inner(T, u_)*dx
def left(x, on_boundary):
return near(x[0], 0.)
bc = DirichletBC(V, Constant((0.,0.,0)), left)
u = Function(V, name="Displacement")
solve(a == l, u, bc)
The deflection shape seems wrong
I also tried with mixed element with 6 DOFs following this example.
from dolfin import *
from ufl import Jacobian, diag
#from plotting import plot
length = 10.0
N = 40 # number of elements
points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
cells = np.array([[i, i+1] for i in range(N)])
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
def tangent(mesh):
t = Jacobian(mesh)
return as_vector([t[0,0], t[1, 0], t[2, 0]])/sqrt(inner(t,t))
t = tangent(mesh)
ez = as_vector([0, 0, 1])
a1 = cross(t, ez)
a1 /= sqrt(dot(a1, a1))
a2 = cross(t, a1)
a2 /= sqrt(dot(a2, a2))
# We now define the material and geometrical constants which will be used in the constitutive relation. We consider the case of a rectangular cross-section of width $b$ and height $h$ in directions $a_1$ and $a_2$. The bending inertia will therefore be $I_1 = bh^3/12$ and $I_2=hb^3/12$. The torsional inertia is $J=\beta hb^3$ with $\beta\approx 0.26$ for $h=3b$. Finally, the shear areas are approximated by $S_1=S_2=\kappa S$ with $\kappa=5/6$.
# In[3]:
thick = Constant(1)
width = thick#/3
E = Constant(1e5)
nu = Constant(0.3)
G = E/2/(1+nu)
rho = Constant(1e-3)
g = Constant(1)
S = thick*width
ES = E*S
EI1 = E*width*thick**3/12
EI2 = E*width**3*thick/12
GJ = G*0.26*thick*width**3
kappa = Constant(5./6.)
GS1 = kappa*G*S
GS2 = kappa*G*S
# We now consider a mixed $\mathbb{P}_1/\mathbb{P}_1$-Lagrange interpolation for the displacement and rotation fields. The variational form is built using a function `generalized_strains` giving the vector of six generalized strains as well as a function `generalized_stresses` which computes the dot product of the strains with the above-mentioned constitutive matrix (diagonal here). Note that since the 1D beams are embedded in an ambient 3D space, the gradient operator has shape (3,), we therefore define a tangential gradient operator `tgrad` by taking the dot product with the local tangent vector $t$.
#
# Finally, similarly to Reissner-Mindlin plates, shear-locking issues might arise in the thin beam limit. To avoid this, reduced integration is performed on the shear part $Q_1\widehat{\gamma}_1+Q_2\widehat{\gamma}_2$ of the variational form using a one-point rule.
# In[6]:
Ue = VectorElement("CG", mesh.ufl_cell(), 1, dim=3)
W = FunctionSpace(mesh, Ue*Ue)
u_ = TestFunction(W)
du = TrialFunction(W)
(w_, theta_) = split(u_)
(dw, dtheta) = split(du)
def tgrad(u):
return dot(grad(u), t)
def generalized_strains(u):
(w, theta) = split(u)
return as_vector([dot(tgrad(w), t),
dot(tgrad(w), a1)-dot(theta, a2),
dot(tgrad(w), a2)+dot(theta, a1),
dot(tgrad(theta), t),
dot(tgrad(theta), a1),
dot(tgrad(theta), a2)])
def generalized_stresses(u):
return dot(diag(as_vector([ES, GS1, GS2, GJ, EI1, EI2])), generalized_strains(u))
Sig = generalized_stresses(du)
Eps = generalized_strains(u_)
dx_shear = dx(scheme="default",metadata={"quadrature_scheme":"default", "quadrature_degree": 1})
k_form = sum([Sig[i]*Eps[i]*dx for i in [0, 3, 4, 5]]) + (Sig[1]*Eps[1]+Sig[2]*Eps[2])*dx_shear
l_form = Constant(-rho*S*g)*w_[2]*dx
# Clamped boundary conditions are considered at the bottom $z=0$ level and the linear problem is finally solved.
# In[7]:
def bottom(x, on_boundary):
return near(x[0], 0.)
bc = DirichletBC(W, Constant((0, 0, 0, 0, 0, 0)), bottom)
u = Function(W)
solve(k_form == l_form, u, bc)
ffile = XDMFFile("shell-results.xdmf")
ffile.parameters["functions_share_mesh"] = True
v = u.sub(0, True)
v.rename("Displacement", "")
ffile.close()
The curvature still looks wrong looks too straight.
