Rubber roller compacting analysis: contact pressure plot about felupe HOT 19 CLOSED

import felupe as fem
import numpy as np
import pypardiso

r = 0.4
R = 1.0
thickness = 0.5

mesh = (
    fem.mesh.Line(a=r, b=R, n=6)
    .revolve(37, phi=180)
    .rotate(90, axis=2)
    .expand(n=6, z=thickness / 2)
)
mesh.update(points=np.vstack([mesh.points, [0, -1.1, 0]]))
x, y, z = mesh.points.T

region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])

inner = np.isclose(np.sqrt(x**2 + y**2), r)
outer = np.isclose(np.sqrt(x**2 + y**2), R)
boundaries = {
    **fem.dof.symmetry(field[0], axes=(True, False, True)),
    "move": fem.dof.Boundary(field[0], mask=inner, value=0),
    "bottom-x": fem.dof.Boundary(field[0], fy=-1.1, value=0, skip=(0, 1, 1)),
    "bottom-y": fem.dof.Boundary(field[0], fy=-1.1, value=0, skip=(1, 0, 1)),
    "bottom-z": fem.dof.Boundary(field[0], fy=-1.1, value=0, skip=(1, 1, 0)),
}

umat = fem.Hyperelastic(fem.mooney_rivlin, C10=0.4, C01=0.1)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
bottom = fem.MultiPointContact(
    field,
    points=np.arange(mesh.npoints)[outer],
    centerpoint=-1,
    skip=(1, 0, 1),
)
move = fem.math.linsteps([0, 0.1, 0.3], num=[1, 5])
step = fem.Step(
    items=[solid, bottom], ramp={boundaries["bottom-y"]: move}, boundaries=boundaries
)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["bottom-y"])
job.evaluate(tol=1e-1, solver=pypardiso.spsolve, parallel=True)
fig, ax = job.plot(
    x=move.reshape(-1, 1),
    yaxis=1,
    xlabel="Displacement (Bottom) $u_2$ in mm",
    ylabel="Reaction Force (Center) $F_2$ in N",
)

# Create a region on the boundary and project the Cauchy Stress, located at the
# quadrature points, to the mesh-points.
boundary = fem.RegionHexahedronBoundary(mesh, mask=outer)
stress = fem.project(solid.evaluate.cauchy_stress(field), region).reshape(-1, 9)

# Create a new displacement field only on the boundary and link it to the existing
# displacement field.
displacement = fem.FieldContainer([fem.Field(boundary, dim=3)])
displacement.link(field)

# Obtain the deformation gradient on the surface-quadrature points and evaluate the
# surface normals on the deformed configuration.
F = displacement.extract()[0]
N = boundary.normals
n = fem.math.dot(F, N, mode=(2, 1))
n /= fem.math.norm(n, axis=0)

# Interpolate the Cauchy-stress to the surface-quadrature points and evaluate the
# Cauchy stress-normal component (=contact pressure).
σ = fem.Field(boundary, values=stress).interpolate().reshape(F.shape)
p = -fem.math.ddot(σ, fem.math.dya(n, n, mode=1))

# Project the contact pressure from the surface-quadrature points to the mesh points
# for points which are in contact.
points_in_contact = bottom.results.force.reshape(-1, 3).nonzero()[0]
contact_pressure = np.zeros(mesh.npoints)
contact_pressure[points_in_contact] = fem.project(p, boundary)[points_in_contact]

# View the contact pressure
view = field.view(point_data={"Contact Pressure": contact_pressure})
view.plot("Contact Pressure", show_undeformed=False).show()

import felupe as fem
import numpy as np
import pypardiso

r = 0.4
R = 1.0
thickness = 0.5

mesh = (
    fem.mesh.Line(a=r, b=R, n=6)
    .revolve(37, phi=180)
    .rotate(90, axis=2)
    .expand(n=6, z=thickness / 2)
)

# add control points to the mesh and delete them from "points-without-cells"
# normally, displacements at these points are excluded from the active degrees-of-freedom
mesh.update(points=np.vstack([mesh.points, [0, -1.001, 0], [0, 0, 0]]))
mesh.points_without_cells = mesh.points_without_cells[:0]

x, y, z = mesh.points.T

region = fem.RegionHexahedron(mesh)
field = fem.FieldContainer([fem.Field(region, dim=3)])

inner = np.isclose(np.sqrt(x**2 + y**2), r)
outer = np.isclose(np.sqrt(x**2 + y**2), R)
mask = np.zeros(mesh.npoints, dtype=bool)
mask[-1] = True
boundaries = {
    **fem.dof.symmetry(field[0], axes=(True, False, True)),
    "center": fem.Boundary(field[0], mask=mask, skip=(0, 1, 0)),
    "bottom": fem.dof.Boundary(field[0], fy=-1.001, value=0),
}

umat = fem.Hyperelastic(fem.mooney_rivlin, C10=0.4, C01=0.1)
solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
bottom = fem.MultiPointContact(
    field,
    points=np.arange(mesh.npoints)[outer],
    centerpoint=-2,
    skip=(1, 0, 1),
)
center = fem.MultiPointConstraint(
    field,
    points=np.arange(mesh.npoints)[inner],
    centerpoint=-1,
    skip=(0, 0, 0),
)
load = fem.PointLoad(field, points=[-1])
force = -fem.math.linsteps([0, 1e-10, 0.01, 0.15], num=[1, 1, 5], axis=1, axes=3)
step = fem.Step(
    items=[solid, bottom, load, center], ramp={load: force}, boundaries=boundaries
)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["center"])
job.evaluate(tol=1e-2, solver=pypardiso.spsolve, parallel=True)
fig, ax = job.plot(
    y=force,
    xaxis=1,
    yaxis=1,
    xlabel="Displacement (Center) $u_2$ in mm",
    ylabel="External Force (Center) $F_2$ in N",
)

# Show the minimal-principal value of cauchy stress
ax = solid.imshow(
    "Principal Values of Cauchy Stress",
    component=2,
    show_undeformed=False,
    project=fem.project,
)

# Create a region on the boundary and project the Cauchy Stress, located at the
# quadrature points, to the mesh-points.
boundary = fem.RegionHexahedronBoundary(mesh, mask=outer)
stress = fem.project(solid.evaluate.cauchy_stress(field), region).reshape(-1, 9)

# Create a new displacement field only on the boundary and link it to the existing
# displacement field.
displacement = fem.FieldContainer([fem.Field(boundary, dim=3)])
displacement.link(field)

# Obtain the deformation gradient on the surface-quadrature points and evaluate the
# surface normals on the deformed configuration.
F = displacement.extract()[0]
N = boundary.normals
n = fem.math.dot(F, N, mode=(2, 1))
n /= fem.math.norm(n, axis=0)

# Interpolate the Cauchy-stress to the surface-quadrature points and evaluate the
# Cauchy stress-normal component (=contact pressure).
σ = fem.Field(boundary, values=stress).interpolate().reshape(F.shape)
p = -fem.math.ddot(σ, fem.math.dya(n, n, mode=1))

# Project the contact pressure from the surface-quadrature points to the mesh points
# for points which are in contact.
points_in_contact = bottom.results.force.reshape(-1, 3).nonzero()[0]
contact_pressure = np.zeros(mesh.npoints)
contact_pressure[points_in_contact] = fem.project(p, boundary)[points_in_contact]

# View the contact pressure
view = field.view(point_data={"Contact Pressure": contact_pressure})
view.plot("Contact Pressure", show_undeformed=False).show()

import numpy as np
import matplotlib.pyplot as plt
import pypardiso
import felupe as fem

# geometry
r = 30 / 2  # mm
R = 70 / 2  # mm
t = 2.5  # mm
length = 100  # mm

# load
F = 15000  # N
nsubsteps = 10  # number of substeps

# coating
E = 20  # MPa
poisson = 0.3

# rubber
C10 = 0.4  # MPa
C01 = 0.1  # MPa
K = 5000 * 2 * (C10 + C01)  # MPa

# meshing
nradial = (9, 4)  # number of radial points (rubber, coating)
ntheta = 181  # number of circ. points

# revolve two line-meshes, one for the rubber, one for the coating with equal `ntheta`
meshes = [
    fem.mesh.Line(a=r, b=R - t, n=nradial[0]).revolve(ntheta, phi=360),
    fem.mesh.Line(a=R - t, b=R, n=nradial[1]).revolve(ntheta, phi=360),
]

# add both meshes to a mesh-container, merge the points
container = fem.MeshContainer(meshes, merge=True)
ax = container.imshow(colors=["green", "red"])

# add control points to the meshes
for m in container.meshes:
    m.update(points=np.vstack([m.points, [0, -R * 1.001], [0, 0]]))

# take the container-points from the first mesh and create a stacked-mesh (all cells)
container.points = container.meshes[0].points
mesh = container.stack()

# delete control-points them from "points-without-cells"
# (otherwise the degrees-of-freedoms of these points would be ignored)
mesh.points_without_cells = mesh.points_without_cells[:-2]

# create a top-level region & field (for boundary conditions)
region = fem.RegionQuad(mesh)
field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)])

# create related sub-regions and -fields (for solid bodies)
regions = [fem.RegionQuad(m) for m in container.meshes]
fields = [fem.FieldContainer([fem.FieldPlaneStrain(r, dim=2)]) for r in regions]

# boundary conditions
x, y = mesh.points.T
inner = np.isclose(np.sqrt(x**2 + y**2), r)
outer = np.isclose(np.sqrt(x**2 + y**2), R)
mask = np.zeros(mesh.npoints, dtype=bool)
mask[-1] = True
boundaries = {
    "center": fem.Boundary(field[0], mask=mask, skip=(0, 1)),
    "bottom": fem.dof.Boundary(field[0], fy=-R * 1.001),
}

# hyperelastic solid bodies
rubber = fem.SolidBodyNearlyIncompressible(
    umat=fem.Hyperelastic(fem.mooney_rivlin, C10=C10, C01=C01),
    field=fields[0],
    bulk=K,
)
coating = fem.SolidBody(
    umat=fem.LinearElasticLargeStrain(E=E, nu=poisson),
    field=fields[1],
)

# constraints (rigid bottom-plate and multi-point-c. from inner-radius to center-point)
bottom = fem.MultiPointContact(
    field,
    points=np.arange(mesh.npoints)[outer],
    centerpoint=-2,
    skip=(1, 0),
)
mpc = fem.MultiPointConstraint(
    field,
    points=np.arange(mesh.npoints)[inner],
    centerpoint=-1,
    skip=(0, 0),
)

# point load
load = fem.PointLoad(field, points=[-1])

# force per thickness
force = -fem.math.linsteps(
    [0, 1e-9, F / length],
    num=[1, nsubsteps],
    axis=1,
    axes=2,
)

step = fem.Step(
    items=[rubber, coating, bottom, mpc, load],
    ramp={load: force},
    boundaries=boundaries,
)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["center"])
job.evaluate(x0=field, tol=1e-4, solver=pypardiso.spsolve)
fig, ax = job.plot(
    y=force,
    xaxis=1,
    yaxis=1,
    yscale=length,
    xlabel="Displacement (Center) $u_2$ in mm",
    ylabel="External Force (Center) $F_2$ in N",
)
ax = field.imshow("Principal Values of Logarithmic Strain", project=fem.topoints)

# Create a region on the boundary and shift the Cauchy Stress, located at the
# quadrature points, to the mesh-points
# (instead of `fem.topoints(values, region)`, one could also use `fem.project(values, region)`).
boundary = fem.RegionQuadBoundary(container.meshes[1], mask=outer, ensure_3d=True)
stress = fem.topoints(coating.evaluate.cauchy_stress(fields[1]), regions[1]).reshape(
    -1, 9
)

# Create a new displacement field only on the boundary and link it to the existing
# displacement field.
displacement = fem.FieldContainer([fem.FieldPlaneStrain(boundary, dim=2)])
displacement.link(field)

# Obtain the deformation gradient on the surface-quadrature points and evaluate the
# surface normals on the deformed configuration at the quadrature points.
F = displacement.extract()[0]
N = boundary.normals
n = fem.math.dot(F, N, mode=(2, 1))
n /= fem.math.norm(n, axis=0)

# Interpolate the Cauchy-stress to the surface-quadrature points and evaluate the
# Cauchy stress-normal component (=contact pressure).
σ = fem.Field(boundary, values=stress).interpolate().reshape(F.shape)
p = -fem.math.ddot(σ, fem.math.dya(n, n, mode=1))

# Project the contact pressure from the surface-quadrature points to the mesh points
# for points which are in contact.
points_in_contact = bottom.results.force.reshape(-1, 2).nonzero()[0]
contact_pressure = np.zeros(mesh.npoints)
contact_pressure[points_in_contact] = fem.project(p, boundary)[points_in_contact]

# View the contact pressure
view = field.view(point_data={"Contact Pressure": contact_pressure})
view.plot("Contact Pressure", show_undeformed=False).show()

# Evaluate and plot the contact pressure (length)
mesh_deformed = mesh.copy(points=mesh.points + field[0].values)
contact_coords = mesh_deformed.x[points_in_contact[:-1]]
contact_length = contact_coords.max() - contact_coords.min()

# Ensure that the mesh-points are already sorted
assert np.allclose(np.argsort(contact_coords), np.arange(len(contact_coords)))

fig, ax = plt.subplots()
ax.plot(
    mesh_deformed.x[points_in_contact[:-1]],
    contact_pressure[points_in_contact[:-1]],
    marker="o",
    label=f"Contact Length $\Delta x$ = {np.round(contact_length, 2)} mm",
)
ax.set_xlabel(r"Position $x$ in mm")
ax.set_ylabel(r"Contact Pressure in MPa")
ax.legend()

import numpy as np
import matplotlib.pyplot as plt
import felupe as fem

# geometry
r = 30 / 2  # mm
R = 70 / 2  # mm
t = 2.5  # mm
length = 100  # mm

# load
F = 15000  # N
nsubsteps = 10  # number of substeps

# coating
E = 20  # MPa
poisson = 0.3

# rubber
C10 = 0.4  # MPa
C01 = 0.1  # MPa
K = 5000 * 2 * (C10 + C01)  # MPa

# meshing
nradial = (9, 4)  # number of radial points (rubber, coating)
ntheta = 181  # number of circ. points

# revolve two line-meshes, one for the rubber, one for the coating with equal `ntheta`
meshes = [
    fem.mesh.Line(a=r, b=R - t, n=nradial[0]).revolve(ntheta, phi=360),
    fem.mesh.Line(a=R - t, b=R, n=nradial[1]).revolve(ntheta, phi=360),
]

# add both meshes to a mesh-container, merge the points
container = fem.MeshContainer(meshes, merge=True, decimals=8)
ax = container.imshow(colors=["green", "red"])

# add control points to the meshes
for m in container.meshes:
    m.update(points=np.vstack([m.points, [0, -R * 1.001], [0, 0]]))

# take the container-points from the first mesh and create a stacked-mesh (all cells)
container.points = container.meshes[0].points
mesh = container.stack()

# delete control-points them from "points-without-cells"
# (otherwise the degrees-of-freedoms of these points would be ignored)
mesh.points_without_cells = mesh.points_without_cells[:-2]

# create a top-level region & field (for boundary conditions)
region = fem.RegionQuad(mesh)
field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)])

# create related sub-regions and -fields (for solid bodies)
regions = [fem.RegionQuad(m) for m in container.meshes]
fields = [fem.FieldContainer([fem.FieldPlaneStrain(r, dim=2)]) for r in regions]

# boundary conditions
x, y = mesh.points.T
inner = np.isclose(np.sqrt(x**2 + y**2), r)
outer = np.isclose(np.sqrt(x**2 + y**2), R)
mask = np.zeros(mesh.npoints, dtype=bool)
mask[-1] = True
boundaries = {
    "center": fem.Boundary(field[0], mask=mask, skip=(0, 1)),
    "bottom": fem.dof.Boundary(field[0], fy=-R * 1.001),
}

# hyperelastic solid bodies
rubber = fem.SolidBodyNearlyIncompressible(
    umat=fem.Hyperelastic(fem.mooney_rivlin, C10=C10, C01=C01),
    field=fields[0],
    bulk=K,
)
coating = fem.SolidBody(
    umat=fem.LinearElasticLargeStrain(E=E, nu=poisson),
    field=fields[1],
)

# constraints (rigid bottom-plate and multi-point-c. from inner-radius to center-point)
bottom = fem.MultiPointContact(
    field,
    points=np.arange(mesh.npoints)[outer],
    centerpoint=-2,
    skip=(1, 0),
)
mpc = fem.MultiPointConstraint(
    field,
    points=np.arange(mesh.npoints)[inner],
    centerpoint=-1,
    skip=(0, 0),
)

# point load
load = fem.PointLoad(field, points=[-1])

# force per thickness
force = -fem.math.linsteps(
    [1e-10, 1e-9, 5e-8, F / length],
    num=[1, 1, nsubsteps],
    axis=1,
    axes=2,
)

step = fem.Step(
    items=[rubber, coating, bottom, mpc, load],
    ramp={load: force},
    boundaries=boundaries,
)

job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["center"])
job.evaluate(x0=field, tol=1e-4)
fig, ax = job.plot(
    y=force,
    xaxis=1,
    yaxis=1,
    yscale=length,
    xlabel="Displacement (Center) $u_2$ in mm",
    ylabel="External Force (Center) $F_2$ in N",
)
ax = field.imshow("Principal Values of Logarithmic Strain", project=fem.topoints)

# Create a region on the boundary and shift the Cauchy Stress, located at the
# quadrature points, to the mesh-points
# (instead of `fem.topoints(values, region)`, one could also use `fem.project(values, region)`).
boundary = fem.RegionQuadBoundary(container.meshes[1], mask=outer, ensure_3d=True)
stress = fem.topoints(coating.evaluate.cauchy_stress(fields[1]), regions[1]).reshape(
    -1, 9
)

# Create a new displacement field only on the boundary and link it to the existing
# displacement field.
displacement = fem.FieldContainer([fem.FieldPlaneStrain(boundary, dim=2)])
displacement.link(field)

# Obtain the deformation gradient on the surface-quadrature points and evaluate the
# surface normals on the deformed configuration at the quadrature points.
F = displacement.extract()[0]
N = boundary.normals
n = fem.math.dot(F, N, mode=(2, 1))
n /= fem.math.norm(n, axis=0)

# Interpolate the Cauchy-stress to the surface-quadrature points and evaluate the
# Cauchy stress-normal component (=contact pressure).
σ = fem.Field(boundary, values=stress).interpolate().reshape(F.shape)
p = -fem.math.ddot(σ, fem.math.dya(n, n, mode=1))

# Project the contact pressure from the surface-quadrature points to the mesh points
# for points which are in contact.
points_in_contact = bottom.results.force.reshape(-1, 2).nonzero()[0]
contact_pressure = np.zeros(mesh.npoints)
contact_pressure[points_in_contact] = fem.project(p, boundary)[points_in_contact]

# View the contact pressure
view = field.view(point_data={"Contact Pressure": contact_pressure})
view.plot("Contact Pressure", show_undeformed=False, style="points").show()

# Evaluate and plot the contact pressure (length)
mesh_deformed = mesh.copy(points=mesh.points + field[0].values)
contact_coords = mesh_deformed.x[points_in_contact[:-1]]
contact_length = contact_coords.max() - contact_coords.min()

# Ensure that the mesh-points are already sorted
assert np.allclose(np.argsort(contact_coords), np.arange(len(contact_coords)))

fig, ax = plt.subplots()
ax.plot(
    mesh_deformed.x[points_in_contact[:-1]],
    contact_pressure[points_in_contact[:-1]],
    marker="o",
    label=f"Contact Length $\Delta x$ = {np.round(contact_length, 2)} mm",
)
ax.set_xlabel(r"Position $x$ in mm")
ax.set_ylabel(r"Contact Pressure in MPa")
ax.legend()