Hi,
I am following the following tutorial from @bleyerj Axisymmetric formulation for elastic structures of revolution — Numerical tours of continuum mechanics using FEniCS master documentation and I was wondering, if I wanted to find the principal stresses of the hemisphere when I am done, how can I do that?
I am trying to combine the examples of @bleyerj from indentation and axisymmetric calculations in order to get a matching condition of a sphere, like asked about in Why FEniCS values are so different than theoretical ones? - #11 by bleyerj.
from dolfinx import mesh, fem
import numpy as np
import ufl
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py.PETSc import Options as PETScOptions
from petsc4py.PETSc import ScalarType
N = 200
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N)
domain.geometry.x[:, 0] = domain.geometry.x[:, 0] ** 2
domain.geometry.x[:, 1] = 1 - domain.geometry.x[:, 1] ** 2
V = fem.functionspace(domain, ("CG", 1, ()))
V1 = fem.functionspace(domain, ("CG", 1, (domain.geometry.dim,)))
V2 = fem.functionspace(domain, ("CG", 1, (domain.geometry.dim + 1, domain.geometry.dim + 1)))
V0 = fem.functionspace(domain, ("DG", 0, ()))
def top(x):
return np.isclose(x[1], 1)
def bottom(x):
return np.isclose(x[1], 0)
def symmetry_x(x):
return np.isclose(x[0], 0)
bottom_edge = mesh.locate_entities_boundary(domain, 1, bottom)
top_edge = mesh.locate_entities_boundary(domain, 1, top)
sym_x = mesh.locate_entities_boundary(domain, 1, symmetry_x)
marked = np.hstack([bottom_edge, top_edge, sym_x])
mark_edges = np.hstack([np.full_like(bottom_edge, 1), np.full_like(top_edge, 2), np.full_like(sym_x, 3)])
sort_edges = np.argsort(marked)
edge_tags = mesh.meshtags( domain, 1, marked[sort_edges], mark_edges[sort_edges])
u_zero = np.array((0,) * domain.geometry.dim, dtype=ScalarType)
bc_bot = fem.dirichletbc( u_zero, fem.locate_dofs_topological(V1, 1, edge_tags.find(1)), V1)
sym_x = fem.locate_dofs_topological(V1.sub(0), 1, edge_tags.find(3))
bc_sym_x = fem.dirichletbc(ScalarType(0), sym_x, V1.sub(0))
bc = [bc_bot, bc_sym_x]
ds = ufl.Measure("ds", domain=domain, subdomain_data=edge_tags)
d = 0.02
R = 0.5
def obstacle(x):
return -d + (pow(x[0], 2)) / 2 / R
indent = fem.Function(V, name="indenter")
indent.interpolate(obstacle)
u = fem.Function(V1, name="u")
du = ufl.TrialFunction(V1)
u_ = ufl.TestFunction(V1)
p = fem.Function(V0, name="pressure")
E = fem.Constant(domain, 10.0)
nu = fem.Constant(domain, 0.3)
mu = E / 2 / (1 + nu)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
r = ufl.SpatialCoordinate(domain)[0]
def eps(v):
return ufl.sym(ufl.as_tensor([[v[0].dx(0), 0, v[0].dx(1)], [0, v[0] / r, 0], [v[1].dx(0), 0, v[1].dx(1)]]))
def sigma(v):
return lmbda * ufl.tr(eps(v)) * ufl.Identity(3) + 2.0 * mu * eps(v)
def ppos(x):
return (x + abs(x)) / 2.0
pen = fem.Constant(domain, 1e4)
F = ufl.inner(sigma(u), eps(u_)) * r * ufl.dx + pen * ufl.dot( u_[1], ppos(u[1] - indent)) * r * ds(2)
J = ufl.derivative(F, u, du)
problem = NonlinearProblem(F, u, bc, J=J)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
ksp = solver.krylov_solver
opts = PETScOptions()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}pc_type"] = "hypre"
solver.solve(u)
stress = fem.Function(V2, name="stress")
stress.interpolate( fem.Expression(sigma(u), V2.element.interpolation_points()))
My code does not run with an error, and I do get the correct value of force, but what about the stresses? is this the correct way to calculate the cauchy stress? And also, how would I take these values and calculate the principal stresses?
Principal stresses are the eigenvalues of the stress tensor.
So, please refer to the following discussion:
Hi, thanks for the reply! I have implemented the stress calculations, but when I do so, I am getting NaN values at the symmetry axis (as seen in the tutorial Axisymmetric formulation for elastic structures of revolution — Numerical tours of continuum mechanics using FEniCS master documentation, there are 1/r terms in the strain, leading to 1/r terms in the stress calcs). How can I fix that?
how about using L’Hopital’s rule?
\lim_{r \to 0} \frac{u_r}{r} = \lim_{r \to 0} \frac{du_r}{dr}
I don’t understand how I would implement this in FEniCSx. I want a whole stress field, which I interpolate over the functionspace’s interpolation points, so how would I add in L’Hopital’s rule?
How about making use of the ufl.conditional when you define the strain tensor (eps)?
For instance,
# strain
condition = ufl.ge(x[0], 1e-10)
def eps(v):
gradu = [[v[0].dx(0), 0, v[0].dx(1)],
[0, v[0]/x[0], 0],
[v[1].dx(0), 0, v[1].dx(1)]]
gradu_ = [row[:] for row in gradu] # copy gradu to gradu_
gradu_[1][1] = v[0].dx(0) # L'Hopital's rule
return ufl.conditional(condition,
ufl.sym(ufl.as_tensor(gradu)),
ufl.sym(ufl.as_tensor(gradu_)))
Thank you for this suggestion, it works great! One small difference I had to use was instead of the copy(), I recreated the tensor, as the copy() only copies the highest level list, so both gradu_ and gradu were modified.
Thanks, I revised the code above.