Rubber roller compacting analysis: contact pressure plot #805
-
|
Hello, I would like to model a rubber roller, crushed on a flat plane (see image). The rroller has a rigid sleeve inside, with a load force applied in the center, and I would like to extract the contact length and pressure distribution in contact zone. There is also a coating (thin plastic cover) but I can ignore it for the moment.
Notice that I would like to use Mooney-Rivlin material model, which seems to be supported by Felupe (nice!). I noticed there is the 'morph-rubber-wheel' tutorial, but it is not exactly what I would like to do (and it misses contatc length and contact pressure): Thanks for sharing your amazing library !! |
Beta Was this translation helpful? Give feedback.
Replies: 26 comments 9 replies
-
|
Hi @yvanblanchard, node-to-segment or segment-to-segment contact formulations aren't implemented yet. All contact formulations in the examples are either based on frictionless node-to-rigid contact (which is in fact a compression-only multi-point-constraint, like your referenced example) or with additional solid-contact elements with very low stiffness. I could provide examples with both methods, would that be helpful for you? Both methods would be frictionless. A view on the contact pressure is definitely possible. |
Beta Was this translation helpful? Give feedback.
-
|
Hi @adtzlr Yes, this simplified friction-less contact is fine for me (at least for first trial), so interested by an example how to model it (with both approaches), with contact pressure plotting as well. Thanks again for your help! |
Beta Was this translation helpful? Give feedback.
-
|
Hi @yvanblanchard, here's a first example which uses Codeimport 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()
|
Beta Was this translation helpful? Give feedback.
-
|
Thank you very much for this first fantastic result @adtzlr !! very nice If I want to control the load force: do you advise me to check the curve plot and recation force value, and then adjust (with a second analysis run) the max applied displacement (here 0.3) in order to get the expected force ? |
Beta Was this translation helpful? Give feedback.
-
|
This one is load-controlled. A very small initial force has to be provided. Codeimport 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()
|
Beta Was this translation helpful? Give feedback.
-
No no, but it's much easier to start with. Please see the second example. |
Beta Was this translation helpful? Give feedback.
-
|
Thank you again @adtzlr could you also please tell me how to integrate a kind of coating to the roller ? There is a stiff thin plastic cover in reality, and I would like to model it using two or three rows of elements thick (0.2mm thick). So it’s just about how to build such a mesh around the cylinder (with coincident nodes at interface), and how apply a different material (here kind of isotropic plastic). |
Beta Was this translation helpful? Give feedback.
-
|
FElupe has Isotropic Plasticity implemented - but not for large rotations. Other than that, meshing is straight forward. You'll have to create a https://felupe.readthedocs.io/en/latest/howto/composite.html |
Beta Was this translation helpful? Give feedback.
-
|
Hi @yvanblanchard, are you interested in a 3D simulation or is 2D / Plane Strain fine, i.e. how is the thickness / diameter ratio of your problem? |
Beta Was this translation helpful? Give feedback.
-
|
Hello |
Beta Was this translation helpful? Give feedback.
-
|
Thanks for clarification. Yes, 2d plane strain is possible (and easier + faster). With the given node-to-rigid contact a non-flat plate is not possible, but with very soft elements in between it is. It could be that a displacement controlled approach would be much more stable. I'll check that out. However, I'm still unsure how to handle plasticity of the coating. I'll come up with an example after the weekend. |
Beta Was this translation helpful? Give feedback.
-
|
Great, thanks !
I am definitely interested by the same things you prepared for me, but for plane stress (2D).
About the plastic coating: there is a misunderstanding, i am sorry: I wanted to say that for my application, the silicone rubber has a thin ring of plastic material (ptfe I think), that we can consider as isotropie. so no need of plasticity behavior.
About the possibility to handle contact with a non-planar surface: that is more a best-have and not a must-have for I want to do. So I don’t want you to spend to much time on it. I can use prepomax and calculix solver for this particular case (even if your library, allowing to have everything in à python notebook, without files exchange, is just fantastic!).
Yvan Blanchard
Coriolis Composites
…________________________________
De : Andreas Dutzler ***@***.***>
Envoyé : Saturday, June 29, 2024 10:17:24 AM
À : adtzlr/felupe ***@***.***>
Cc : Yvan BLANCHARD ***@***.***>; Mention ***@***.***>
Objet : Re: [adtzlr/felupe] Rubber roller compacting analysis: contact pressure plot (Issue #800)
Do not click links or open attachments unless you recognize the sender and know the content is safe.
Thanks for clarification.
Yes, 2d plane strain is possible (and easier + faster). With the given node-to-rigid contact a non-flat plate is not possible, but with very soft elements in between it is. It could be that a displacement controlled approach would be much more stable. I'll check that out.
However, I'm still unsure how to handle plasticity of the coating. I'll come up with an example after the weekend.
—
Reply to this email directly, view it on GitHub<#800 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/APHYJYQRHU3TXNA27OJAJMLZJZURJAVCNFSM6AAAAABJ5IHK66VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCOJYGA2DGOJRGA>.
You are receiving this because you were mentioned.Message ID: ***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
|
Are you sure on plane stress? You said the length of the roller is rather long - I'd use plane strain. Anyway, I'll prepare an example. |
Beta Was this translation helpful? Give feedback.
-
|
Sorry, I wanted to say ‘plain strain’! Here i also a view of what I want to achieve:
|
Beta Was this translation helpful? Give feedback.
-
|
Hi @yvanblanchard! Here is the initial version of the script (typo included!).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()Here's the improved script.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 = 500 * 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()Here's a version with a displacement-controlled pre-compression (improves convergence).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
initial_compression = (R - r) / 2 # mm
F = 10000 # N
# coating
E = 20 # MPa
poisson = 0.3
# rubber
C10 = 0.18 # MPa
C01 = 0.05 # 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.1], [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]
# 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 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),
)
# 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, coating, bottom, mpc],
boundaries=boundaries_pre,
)
job_pre = fem.CharacteristicCurve(steps=[step_pre], boundary=boundaries_pre["center"])
job_pre.evaluate(x0=field)
# 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, coating, bottom, mpc, load],
boundaries=boundaries,
)
job = fem.Job(steps=[step], boundary=boundaries["center"])
job.evaluate(x0=field)
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)))
# 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"u(center) = {u[-1]} 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[:-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()Unfortunately, you can't see too much in the plot of the contact pressure. Hence, I added a plot for the contact pressure as a function of the deformed x-coordinate on the outer radius. The value of the initial (small) external force is crucial for the convergence of the simulation. To improve convergence, we could use two consecutive steps, one with e.g. 1-2 mm displacement-controlled compression and then continue with another force-controlled step on top of the initial compression. But the script will get longer and longer, so it depends on what you want to achieve. I hope I could help and the provided example is useful for your problem. Would it be ok for you if I convert this issue to a discussion? The nice thing would be that we can select one of our comments as the answer.
More Plots
|
Beta Was this translation helpful? Give feedback.
-
|
Thank you very much @adtzlr I agree that even if convergence is not optimized here, it is not a big deal for my need. Unfortunately, when running the script (through Jupyter), I have an error (newtonRaphson):
|
Beta Was this translation helpful? Give feedback.
-
|
No, it is not related to Jupyter. I've seen that too if I changed any of the parameters. Unfortunately, the simulation is super unstable for the initial node-to-rigid contact. We could try to improve this by an initial displacement controlled compression. The nearly-incompressible solid (mixed-field formulation) also makes the simulation more unstable. How is the bulk vs. shear modulus of your material? Nearly-incompressible or more like a compressible foam? |
Beta Was this translation helpful? Give feedback.
-
|
@adtzlr It's strange that keeping your python code as-it (same input parameters), I am not able to achieve convergence and your good results. Anyway, to answer your questions:
|
Beta Was this translation helpful? Give feedback.
-
|
Thanks for the details! Unfortunately I haven't implemented the anisotropic linear-elastic material yet. More important: Sorry, there was a typo in the script - it should be improved now. I had deactivated the x-direction of the MPC and added a contact at y=0 instead 🙅🏻 ...the model's center was not fixed in direction x. I edited the above answer. If there are rigid-body motions possible in the model, sometimes Pypardiso finds a solution on one platform (Win) and not on the other one (e.g. Linux). I switched to the SciPy solver because the performance is similar and you don't need an external dependency here. |
Beta Was this translation helpful? Give feedback.
-
|
@yvanblanchard I have moved this conversation to the discussion section. |
Beta Was this translation helpful? Give feedback.
-
|
Thank you @adtzlr |
Beta Was this translation helpful? Give feedback.
-
|
No, it was not in the definition of the boundaries but in the skip-arguments of mpc/contact instead. |
Beta Was this translation helpful? Give feedback.
-
|
Thanks. There seems to be a difference between Win/Linux or other dependencies. I tried it in Colab and reducing the bulk modulus to 500x the shear modulus helps. |
Beta Was this translation helpful? Give feedback.
-
|
Yes, I changed the bulk factor (K variable), using 'K = 500 * (C10+C01)' formula, and it converges. Thanks for the tip. But if I use other (smaller) values for C10 and C01 (i.e. C10 = 0.18, C01 = 0.05 MPa) with this same modified bulk formula, it does not converge anymore. NB: I used a formula linking the rubber 'shore' indicator to the equivalent elastic modulus, and then C01/C10 coefficients, see below.
|
Beta Was this translation helpful? Give feedback.
-
|
I just tested those values of C01/C10, with number of steps of 25 (instead of 10, initially), and it converges well. |
Beta Was this translation helpful? Give feedback.
-
|
@adtzlr |
Beta Was this translation helpful? Give feedback.
-
|
Hi @yvanblanchard, I added a third version of the script to the comment with the scripts #805 (comment). This one includes a displacement-controlled pre-compression of the roller which eliminates the convergence problems. Regarding the displacements of the roller: the initial gap between the rigid plate and the roller must be subtracted from the displacement field values. Also, the previously given # 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"u(center) = {u[-1]} 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") |
Beta Was this translation helpful? Give feedback.
-
|
@adtzlr |
Beta Was this translation helpful? Give feedback.
-
|
Sorry @adtzlr I notice that calculations do not converge for 'softer' material (foam). But maybe it's because I don't have the proper strain formulation/law for my rubber material. Or maybe it's because this 'M-R' formulation is not appropriate for 'foam' roller material. First problem: I tried such elastic material property for rubber, but the contact pressure is null: Secondly: I am supposed to get a radial deformation (from paper/experiemental testing) of bit less of 5mm (for 500 Newtons load force , 60mm roller length and E roller = 0.8 MPa), but I only get 3.2mm with calculations.
|
Beta Was this translation helpful? Give feedback.
-
|
In my opinion, the results of this simulation are quite similar to the ones from the paper. I quickly adjusted both The contact length of this simulation is shorter and hence, the contact pressure is higher compared to the paper. This seems to be a result of the anisotropic vs. isotropic modeling approach. Beside that, specific hyperelastic foam models exist, like the Blatz-Ko model but I have no practical experience with it. Scriptimport numpy as np
import matplotlib.pyplot as plt
import felupe as fem
# geometry
r = 30 / 2 # mm
R = 68.3 / 2 # mm
t = 2.5 # mm
length = 60 # mm
# load
initial_compression = (R - r) / 4 # mm
F = 750 # N
# coating
E = 20 # MPa
poisson = 0.3
# rubber foam
G = 0.65 # MPa
nu = 0.25
# meshing
nradial = (11, 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.1], [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]
# 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 bodies
# rubber = fem.SolidBodyNearlyIncompressible(
# umat=fem.Hyperelastic(fem.mooney_rivlin, C10=C10, C01=C01),
# field=fields[0],
# bulk=K,
# )
rubber = fem.SolidBody(
umat=fem.LinearElasticLargeStrain(E=2 * G * (1 + nu), nu=nu),
field=fields[0],
)
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),
)
# 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, coating, bottom, mpc],
boundaries=boundaries_pre,
)
job_pre = fem.CharacteristicCurve(steps=[step_pre], boundary=boundaries_pre["center"])
job_pre.evaluate(x0=field)
# 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, coating, bottom, mpc, load],
boundaries=boundaries,
)
job = fem.Job(steps=[step], boundary=boundaries["center"])
job.evaluate(x0=field)
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)))
# 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"u(center) = {u[-1]} 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[:-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()F = 500 N: F = 750 N: |
Beta Was this translation helpful? Give feedback.
-
|
Thank you @adtzlr , it's much clearer now Problem is that we don't have the same deformation value as in this paper. About Blatz-Ko model (for foams): I will try to implement it in order to get better contact length (and pressure values). Thank you again. |
Beta Was this translation helpful? Give feedback.
-
|
Blatz-Ko is definitely possible with Hyperelastic. We use a totally different material model (Mooney-Rivlin vs. orthotropic-linear-elastic) for the foam in this script compared to the paper, with nearly-incompressibility instead of a very compressible foam. Differences are to be expected. The compressible Neo-Hookean material is an isotropic material - so you can't copy/paste the material parameters from the paper. One thing what could still be a problem with this script is the bending behavior of the default quad element formulations in this example. |
Beta Was this translation helpful? Give feedback.
-
|
In the very last test/attempt, we used isotropic linear material ( and no longer Mooney-Rivlin). About your quad element formulation potential issue: |
Beta Was this translation helpful? Give feedback.
-
|
This is right, in the r-t plane this is possible. I've overlooked that - sorry. |
Beta Was this translation helpful? Give feedback.
Hi @yvanblanchard,
I added a third version of the script to the comment with the scripts #805 (comment). This one includes a displacement-controlled pre-compression of the roller which eliminates the convergence problems. Regarding the displacements of the roller: the initial gap between the rigid plate and the roller must be subtracted from the displacement field values. Also, the previously given
tol-argument ofJob.evaluate()isn't necessary anymore. (Of course it can still be specified to quit the Newton iterations earlier. This saves some computation time but produces more unbalanced, i.e. out-of-equilibrium, residual forces.)