Hello,
I am attempting to solve the following problem:
A 2D beam is clamped on its left end subjected a tensile load in positive x-direction on its right end. An inner square in the middle of the beam consists of a stiffer material.
The displacement solution is in line with expectation and other finite element software. The stress distribution however shows some unexpected behaviour. Below is a plot of the normal stress in x-dimension. The lower material boundary behaves as expected, with the inner material showing significantly higher stress values due to its higher stiffness. However, the upper boundary does not reflect this trend.
I am not sure what the reason for this is, but I would suspect that it has something to do with the interpolation of the stresses? Below is my code to reproduce the plot above. I appreciate any help.
import numpy as np
import ufl
import matplotlib.pyplot as plt
from dolfinx import fem, mesh, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
# structural parameters
Lx = 2
Ly = 1
# number of elements per unit length in each dimension
N = 100
# lower left and upper right corner of inner square
ll = [0.75, 0.25]
ur = [1.25, 0.75]
# force on right end
f_right = (0.5, 0.0)
tol = 0
element_degree = 1
quadrature_degree = 4
# material parameters
matmod = 'Hooke'
E1 = default_scalar_type(1.)
E2 = default_scalar_type(10.)
nu = default_scalar_type(0.3)
# mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [Lx, Ly]], [Lx*N, Ly*N], mesh.CellType.quadrilateral)
Q = fem.functionspace(domain, ("DG", 0))
V = fem.functionspace(domain, ("CG", element_degree, (domain.geometry.dim, )))
# left end of beam
def left(x):
return np.isclose(x[0], 0)
# right end of beam
def right(x):
return np.isclose(x[0], Lx)
# outer material
def domain1(x):
return np.logical_or(np.logical_or(x[0] <= ll[0] + tol, x[0] >= ur[0] - tol), np.logical_or(x[1] <= ll[1] + tol, x[1] >= ur[1] - tol))
# inner material
def domain2(x):
return np.logical_and(np.logical_and(ll[0] - tol <= x[0], x[0] <= ur[0] + tol), np.logical_and(ll[1] - tol <= x[1], x[1] <= ur[1] + tol))
# mark cells for bcs
fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, left)
right_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag_bc = mesh.meshtags(domain, domain.topology.dim - 1, marked_facets[sorted_facets], marked_values[sorted_facets])
# define mu and lambda differently for each material
cells1 = mesh.locate_entities(domain, domain.topology.dim, domain1)
cells2 = mesh.locate_entities(domain, domain.topology.dim, domain2)
mu = fem.Function(Q)
lam = fem.Function(Q)
mu.x.array[cells1] = np.full_like(cells1, E1 / (2 * (1 + nu)), dtype=default_scalar_type)
mu.x.array[cells2] = np.full_like(cells2, E2 / (2 * (1 + nu)), dtype=default_scalar_type)
lam.x.array[cells1] = np.full_like(cells1, E1 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)
lam.x.array[cells2] = np.full_like(cells2, E2 * nu / ((1 + nu) * (1 - 2 * nu)), dtype=default_scalar_type)
# dirchlet bc
u_zero = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
dirichlet_bc = [fem.dirichletbc(u_zero, left, V)]
# test and solution functions
v = ufl.TestFunction(V)
u = fem.Function(V)
# traction
t = fem.Constant(domain, default_scalar_type(f_right))
# spatial dimension
d = len(u)
# identity tensor
I = ufl.variable(ufl.Identity(d))
# lagrange strain
eps = ufl.sym(ufl.grad(u))
# strain energy
sig = 2 * mu * eps + lam * ufl.tr(eps) * I
psi = 0.5 * ufl.inner(eps, sig)
# set gauß quadrature degree
metadata = {'quadrature_degree': quadrature_degree}
# measure inner domain
dx = ufl.Measure('dx', domain=domain, metadata=metadata)
# measure right end
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag_bc, metadata=metadata)
# potential energy (integrate strain energy on inner and outer material, traction energy on right boundary)
Pi = psi*dx - ufl.dot(t, u)*ds(2)
# variation of Pi
dPidu = ufl.derivative(Pi, u, v)
# define problem F(u) = 0
problem = NonlinearProblem(dPidu, u, dirichlet_bc)
# set solver and options
solver = NewtonSolver(domain.comm, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
# solve
solver.solve(u)
# get stress
stress = 2 * mu * eps + lam * ufl.tr(eps) * I
# interpolate normal stress in x
V = fem.functionspace(domain, ("CG", element_degree))
sxx_expr = fem.Expression(stress[0, 0], V.element.interpolation_points())
sxx = fem.Function(V)
sxx.interpolate(sxx_expr)
# plot normal stress in x
x = np.array(V.tabulate_dof_coordinates())
plt.figure(1)
plt.scatter(x[:, 0], x[:, 1], c=sxx.vector.array, marker='s')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar(label='sigma_xx')
plt.show()
