Result not conform to FEA (elastic roller deformation)

[yvanblanchard]

Hello,

I used a script for computing compaction of an elastic (rubber) roller on a flat surface, but results are too different compared to reference results (FEA and experimental measures).

Results of reference (with load force = 380 N) :

• Contact width = 30 mm
• Max pressure = 0.62 MPa
  NB: diff between FEA and experiments are mainly caused by inaccuracies of measure technique.
  

Felupe analysis results:

• Contact width = 49 mm
• Max pressure = 0.43 MPa
  (with 9.5mm of radial displacement).
  

NB:

• roller diameter = 70mm, inner diam = 35mm, width = 28mm.
• Elastic modulus = 0.275 MPa, and poisson = 0.475
• Use 'SolidBodyNearlyIncompressible' with 'NeoHooke' material model (rubber).

Please find attached the Python script.
roller_ShiminLu-script.zip

Thank you.

[adtzlr]

Hi @yvanblanchard!

This is generally out-of-scope because you're asking for some kind of engineering consulting instead of FElupe support. Anyway, what I would do:

In a first step, compare the force-displacement curves with your referenced article before looking at contact pressure and length. If these curves don't fit, create the model they used (3d, not 2d). If there are still differences, they could have added friction to the contact or they use a different material formulation (linear-elasticity with large strains is not unique) or the material parameters are simply wrong (by mistake, can happen). In Table 3 of your reference, the modulus is simply called "modulus". What if this was the shear and not young's modulus? I created a 3d-model, and obtained the 28 HA result with the 40 HA parameters. Something is strange here. I don't know.

code

import felupe as fem
import numpy as np
import pypardiso

rect = fem.Rectangle(a=(0, 35 / 2), b=(14, 70 / 2), n=(11, 16))
mesh = fem.MeshContainer(
    [
        rect.revolve(n=16, phi=90).rotate(90, axis=0),
        rect.revolve(n=31, phi=90).rotate(180, axis=0),
    ],
    merge=True,
    decimals=4,
).stack()
mesh.update(points=np.vstack([mesh.points, [0, 0, -70 / 2]]))

region = fem.RegionHexahedron(mesh)
field = fem.FieldsMixed(region, n=3)
umat = fem.ThreeFieldVariation(
    fem.LinearElasticLargeStrain(E=2 * 0.275 * (1 + 0.475), nu=0.475)
)
solid = fem.SolidBody(umat, field)

mask = np.isclose(mesh.y**2 + mesh.z**2, 35**2 / 4)
boundaries = dict(
    inner=fem.Boundary(field[0], mask=mask, skip=(0, 0, 1)),
    move=fem.Boundary(field[0], mask=mask, skip=(1, 1, 0)),
    **fem.dof.symmetry(field[0], axes=(1, 1, 0)),
)

points = np.where(np.isclose(mesh.y**2 + mesh.z**2, 70**2 / 4))[0][:-1]
contact = fem.MultiPointContact(field, points=points, centerpoint=-1, skip=(1, 1, 0))

table = fem.math.linsteps([0, -6], num=6)
step = fem.Step(
    [solid, contact], ramp={boundaries["move"]: table}, boundaries=boundaries
)
job = fem.CharacteristicCurve([step], boundary=boundaries["move"])
job.evaluate(solver=pypardiso.spsolve, tol=1e-2, parallel=True)

field.plot("Displacement", component=None, plotter=contact.plot()).show()
field.plot("Principal Values of Logarithmic Strain").show()
solid.plot("Principal Values of Cauchy Stress").show()
fig, ax = job.plot(
    xaxis=2,
    yaxis=2,
    yscale=4,
    xlabel="Displacement in mm",
    ylabel="Force in N",
)

[adtzlr]

And to politely remind you: Please use the discussions-feature for these types of questions, see CONTRIBUTING.md. Thanks 🙏🏻

[yvanblanchard]

Thank you very much @adtzlr
Sorry for this request but I wanted to know if the FElupe lib had a limitation of this case, or if maybe to help identifying an issue in the lib itself (trying to contribute to improve the lib).
I know that that the lib does not manage friction (that's why my reference from litterature don't use neither).
The elastic modulus is a Young equivalent modulus, not shear.
And I was able to reproduce the results of this reference at my side (using FEA and Abaqus), so no typos/errors in the ref paper.
Thank you again for your help, and I will still continue using FEA method for this case.

[adtzlr]

Thats great! I tried to reproduce this with another popular commercial software and the results were identical to FElupe. Interesting! What element type did you use in Abaqus?

[adtzlr]

Anyway, that is a great use case. I'll try to figure out what is the difference. I have some ideas 💡

[adtzlr]

@yvanblanchard the docs of Abaqus: Don't use the linear-elastic material for large-strain analyses; elastic strains should be smaller than 5%. Use a hyperelastic material instead.