Compaction rubber roller analysis issue (with linear isotropic material assumption)

[yvanblanchard]

Hello,

I tried to model the compaction of a rubber roller but I face a convergence issue (max number of iterations reached).

Material: rubber , assumed as isotropic linear material (without any coating/cover material), considering that deformation is small enough (less than 20%).

Here is the python notebook file:
RollerCompaction-FEA_Felupe-Jiang.zip

Paper/article (input values and ref results):
jiang2021.pdf

Thank you for your help!

[adtzlr]

Hi @yvanblanchard,

I won't give support and won't have a look on notebooks - sorry. Please use text-based code-blocks for your questions.

```python
# your code goes here...
```

Nearly-incompressible materials

For nearly incompressible materials like the silicone rubber of the provided paper, you have to use SolidBodyNearlyIncompressible to avoid volumetric locking (Abaqus equivalent should be CPE4H, not CPE4RH) - even in case of relatively small strains. Of course there are other ways to overcome this issue but this is definitely the easiest one. SolidBodyNearlyIncompressible requires a strain energy density function only for the distortional part of the deformation. The volumetric part is added in the solid body by a piecewise (per element) constant pressure field. In the provided paper, a poisson ratio $\nu=0.4995$ results in a ratio of the bulk modulus and shear modulus $K/G = 1000$ (see Lame parameters).

# rubber
E = 1.4  # MPa
nu = 0.4995

# ...

rubber = fem.SolidBodyNearlyIncompressible(
    umat=fem.NeoHooke(mu=E / (2 * (1 + nu))),
    field=field,
    bulk=E / (3 * (1 - 2 * nu)),  # MPa
)

# ...

Here's the complete script.

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

# geometry
r = 20  # mm
R = 40  # mm
length = 1  # mm

# load
initial_compression = (R - r) / 4  # mm
F = 5  # N

# rubber
E = 1.4  # MPa
nu = 0.4995

# meshing
nradial = 11  # number of radial points
ntheta = 181  # number of circ. points

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

# add control points to the mesh
mesh.update(points=np.vstack([mesh.points, [0, -R * 1.1], [0, 0]]))

# remove 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 region & field (for boundary conditions)
region = fem.RegionQuad(mesh)
field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)])

# masks
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)

# hyperelastic solid body
rubber = fem.SolidBodyNearlyIncompressible(
    umat=fem.NeoHooke(mu=E / (2 * (1 + nu))),
    field=field,
    bulk=E / (3 * (1 - 2 * nu)),  # MPa
)

# 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),
)

# initial boundary conditions
boundaries_pre = {
    "center": fem.Boundary(field[0], fx=0, fy=0, mode="and", skip=(0, 1)),
    "move": fem.Boundary(
        field[0],
        fx=0,
        fy=0,
        mode="and",
        skip=(1, 0),
        value=-R * 0.1 - initial_compression,
    ),
    "bottom": fem.dof.Boundary(field[0], fy=mesh.y.min()),
}


# initial compression load case
step_pre = fem.Step(
    items=[rubber, bottom, mpc],
    boundaries=boundaries_pre,
)

job_pre = fem.CharacteristicCurve(steps=[step_pre], boundary=boundaries_pre["center"])
job_pre.evaluate()

# boundary conditions
boundaries = {
    "center": fem.Boundary(field[0], fx=0, fy=0, mode="and", skip=(0, 1)),
    "bottom": fem.dof.Boundary(field[0], fy=mesh.y.min()),
}

# point load
load = fem.PointLoad(field, points=[-1], values=[0, -F / length])

# load case
step = fem.Step(
    items=[rubber, bottom, mpc, load],
    boundaries=boundaries,
)

job = fem.Job(steps=[step], boundary=boundaries["center"])
job.evaluate()
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(mesh, mask=outer, ensure_3d=True)
stress = fem.topoints(rubber.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.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",
    plotter=field.plot(color="white", plotter=bottom.plot(line_width=3, opacity=0.1)),
).show()

# Evaluate and plot the contact pressure (length)
mesh_deformed = mesh.copy(points=mesh.points + field[0].values)
points_in_contact = [
    points_in_contact[0] - nradial,
    *points_in_contact[:-1],
    points_in_contact[-2] + nradial,
]
contact_coords = mesh_deformed.x[points_in_contact]
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)))

# Print displacements of the center of the roller
# a) the displacement field is the first field of the field container
# b) the center point is the last point of the mesh (manually added to the mesh points)
# c) the initial gap is subtracted from the displacement field
u = field[0].values - np.array([[0, -0.1 * R]])
print(f"Δd(center) = {u[-1, 1]:.2f} mm")

# Print deformed coordinates on the outer radius
# a) the displacement values are added to the mesh point coordinates
# b) the outer-mask is applied
X = mesh.points
x = X + u
# print(f"x(outer radius) = {x[outer]} mm")

fig, ax = plt.subplots()
ax.plot(
    mesh_deformed.x[points_in_contact],
    contact_pressure[points_in_contact],
    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()

100%|██████████| 1/1 [00:01<00:00,  1.13s/substep]
100%|██████████| 1/1 [00:02<00:00,  2.08s/substep]
Δd(center) = -2.05 mm

This time, the results from the paper

• center displacement $\Delta d = 2.01$ mm,
• contact length and
• max. contact pressure $p_\text{max} = 0.308$ MPa

are very similar to the ones from this script.

I slightly modified the contact pressure plot to include the zero values at the start and the end for a more precise contact length value.

P.S.: This simulation model is still stable if you refine the mesh.

nradial = 21  # number of radial points
ntheta = 721  # number of circ. points

[yvanblanchard]

Thank you very much for your help.
And sorry for the notebook, I missed that.

[adtzlr]

Thank you very much for your help.
And sorry for the notebook, I missed that.

You're welcome! Sorry for not having a look on the notebook, it's just to make future references to / re-reading of this discussion easier.