How to extract (vector valued) displacement magnitude + questionable highest stress not in contact area

dolfinx
1

Hello everyone,

I am working on a hyperelastic contact model using an analytically defined circular indenter within the penalty approach.

Question 1: How can i automatically evaluate maximum displacements and stresses in general?
This Example in the FEniCS legacy used the dot product and a projection for accomplishing this task. How can it be done using FEniCSx?

Displacements

image
image781×226 16.6 KB

Question 2: How can the highest stress occour inside the bar? Paraview automatically visulizes the stress tensor using a maginitude (unusual?) - How can I avoid it? Should I only use the normal projection for contact stresses?

Cauchy Stress

image
image821×221 16.2 KB

Here’s the MWE as short as possible

import numpy as np
import ufl
from dolfinx import fem, mesh, log
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Options
from dolfinx import nls

# domain geometries
length = 0.25
height = 0.02

# dimensions of stiff indenter 
D = 0.15
center_x = length / 2
center_y = (height + (D/2)) -  0.01

# build rectangular domain 
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0],[length, height]], [40, 9])

x = ufl.SpatialCoordinate(domain)
n = ufl.FacetNormal(domain)

def bottom(x):
    return np.isclose(x[1], 0)

def top(x):
    return np.isclose(x[1], height)

tdim = domain.topology.dim
fdim = tdim - 1

bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
top_facets = mesh.locate_entities_boundary(domain,fdim,top)

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

metadata = {"quadrature_degree":4}
ds = ufl.Measure('ds',domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
V = fem.VectorFunctionSpace(domain, ("CG",1))

u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)

# fixed bottom as dirichlet bc
u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_tag.find(1)), V)

# definition of gap function - circular indenter
def gapp(u, x):
    R = ufl.sqrt(pow(((x[0] - center_x) + u[0]), 2) + pow(((x[1] - center_y) + u[1]), 2)) 
    return  R - (D/2)

# maculay bracket for contact enforcement
def maculay(x): 
    return (x+abs(x))/2

# test and trial functions
u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)

######### kinematic #########
# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))


# hyperelastic material properties, according to demo
E = ScalarType(1.0e4)
nu = ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
P = ufl.diff(psi, F)
S = 2 * ufl.diff(psi, C)

# create XDMF output
xdmf = XDMFFile(domain.comm, "ironing/output.xdmf", "w")
xdmf.write_mesh(domain)

# penalty stiffness for the contact
pen = fem.Constant(domain, ScalarType(1e6))

# weak formulation to solve for the nonlinear problem
form = ufl.inner(ufl.grad(u_), P)*dx + pen * ufl.inner(u_, n*maculay(-gapp(x, u)))*ds(2)

# set up nonlinear problem to solve
Jac = ufl.derivative(form, u, du)
problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=Jac)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-6
solver.rtol = 1e-6
solver.convergence_criterion = "incremental"
solver.report = True
ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# solver logging on
log.set_log_level(log.LogLevel.INFO)

# solver solve
w, converged = solver.solve(u)
assert(converged)
xdmf.write_function(u)
print(f"Number of interations: {w:d}")

# postprocessing: True stress / Cauchy stress
V_stress = fem.TensorFunctionSpace(domain, ("DG", 0))
cauchy_stress = (1./J) * F*S*F.T
stress_expr = fem.Expression(cauchy_stress, V_stress.element.interpolation_points())
cauchystress = fem.Function(V_stress, name="CauchyStress")
cauchystress.interpolate(stress_expr)
xdmf.write_function(cauchystress)

xdmf.close()

Thanks in advance!

2

See Electromagnetics example FEniCSx and

1 Like