-
Notifications
You must be signed in to change notification settings - Fork 8
This issue was moved to a discussion.
You can continue the conversation there. Go to discussion →
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Rubber roller compacting analysis: contact pressure plot #800
Comments
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. |
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! |
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() |
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 ? |
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() |
No no, but it's much easier to start with. Please see the second example. |
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). |
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 |
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? |
Hello |
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. |
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: ***@***.***>
|
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. |
Sorry, I wanted to say ‘plain strain’! Here i also a view of what I want to achieve: |
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 = 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() 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. |
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): |
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? |
@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: |
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. |
This issue was moved to a discussion.
You can continue the conversation there. Go to discussion →
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.
Is it possible with Felupe lib, and if yes could you please indicate me the main steps ?
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 !!
very helpful
The text was updated successfully, but these errors were encountered: