Stresses in Axisymmetric Calculations

General
2

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?