Why do non-interacting cylinders bend in unintuitive direction under identical forces?

[TomM]

I have two cylinders which are not supposed to interact with each other (yet). I apply a pressure at the right face of both cylinders. Why do they curve like this instead of going straight?

If I remove the second cylinder it looks like this (as expected)

If I move them further away from each other this effect increases as well as the elongation only in x direction

It does not matter where I put the 2nd cylinder, the direction in which it gets deformed is always the same. Here I put the 2nd cylinder on the z axis and still the first one is deformed in negative y and negative z direction.

This is my code to reproduce this:

import numpy as np # type: ignore
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction # type: ignore
from mpi4py import MPI # type: ignore
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.fem import locate_dofs_geometrical, DirichletBC
from petsc4py import PETSc # type: ignore
import matplotlib.pyplot as plt # type: ignore
import gmsh

pressure_value = 8.0
E = 10 # Young's modulus
radius = 1
radius_2 = 1
length = 10

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

gdim = 3

# add first cylinder
gmsh.model.add("cylinder")
cylinder_tag = gmsh.model.occ.addCylinder(0, 0, 0, length, 0, 0, radius)
gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(gdim, [cylinder_tag], 0)

surfaces = gmsh.model.getBoundary([(3, cylinder_tag)], oriented=False, recursive=False)

right_face_tag = surfaces[1][1]
gms.model.addPhysicalGroup(gdim - 1, [right_face_tag], 1)
gmsh.model.setPhysicalName(gdim - 1, 1, "Right")

left_face_tag = surfaces[0][1]
gmsh.model.addPhysicalGroup(gdim - 1, [left_face_tag], 2)
gmsh.model.setPhysicalName(gdim - 1, 2, "Left")

lateral_surface_tags = [surf[1] for surf in surfaces[2:]]
gmsh.model.addPhysicalGroup(gdim - 1, lateral_surface_tags, 3)
gmsh.model.setPhysicalName(gdim - 1, 3, "Lateral")

# add second cylinder
second_cylinder_tag = gmsh.model.occ.addCylinder(0, 0,20*radius, length, 0, 0, radius_2)
gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(gdim, [second_cylinder_tag], 1)
surfaces = gmsh.model.getBoundary([(3, second_cylinder_tag)], oriented=False, recursive=False)

second_right_face_tag = surfaces[1][1]
gmsh.model.addPhysicalGroup(gdim - 1, [second_right_face_tag], 4)
gmsh.model.setPhysicalName(gdim - 1, 4, "secondRight")

second_left_face_tag = surfaces[0][1]
gmsh.model.addPhysicalGroup(gdim - 1, [second_left_face_tag], 5)
gmsh.model.setPhysicalName(gdim - 1, 5, "secondLeft")

second_lateral_surface_tags = [surf[1] for surf in surfaces[2:]]
gmsh.model.addPhysicalGroup(gdim - 1, second_lateral_surface_tags, 6)
gmsh.model.setPhysicalName(gdim - 1, 6, "secondLateral")

gmsh.model.mesh.generate(gdim)

domain, _, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim)

gmsh.finalize()

shape = (gdim,)

# Define the function space
V = fem.functionspace(domain, ("P", 2, shape))

u = TrialFunction(V)
v = TestFunction(V)
u_h = fem.Function(V, name="Displacement")

# Define strain and stress tensors
def epsilon(u):
return sym(grad(u))

def sigma(u):
return E * epsilon(u)

pressure_force = fem.Constant(domain, np.array([pressure_value, 0, 0]))

dx = Measure("dx", domain=domain) # Volume measure
ds = Measure("ds", domain=domain, subdomain_data=facets) # Surface measure for boundaries

right_face = 1
right_face_second = 4

# Define the boundary conditions
def left(x):
return np.isclose(x[0], 0)

left_dofs = locate_dofs_geometrical(V, left)
bcs = [
fem.dirichletbc(np.zeros(3,), left_dofs, V) # keep the left face fixed
]

# Define the weak form
a = (inner(sigma(u), epsilon(v)) * dx)
L = inner(pressure_force, v) * (ds(right_face) + ds(right_face_second))
# in case of a single cylinder:
# L = inner(pressure_force, v) * ds(right_face)

bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

problem = fem.petsc.LinearProblem(
a, L, u=u_h, bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)

vtk = io.VTKFile(domain.comm, "results/dynamic/dynamic_cylinder.pvd", "w")

# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)

assemble_vector(b, linear_form)

# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

# Solve linear problem
solver.solve(b, u_h.x.petsc_vec)
u_h.x.scatter_forward()
vtk.write_function(u_h)

vtk.close()

[dokken]

Please note that your code is not runnable (as there are typos and missing indentation).

Secondly, fixing those and inspecting your variational form, it doesn’t look like a standard elasticity law, as you only have a(u,v) = \int_\Omega E \frac{1}{2}(\nabla u + \nabla u^T)\frac{1}{2}(\nabla v + \nabla v^T)~\mathrm{d}x.

What does this stem from, and what are the physically correct units for E and the pressure on the rhs?

Please also note that as your problem is posed on a very coarse mesh, it is probably not radially symmetric, which could cause forces in x and y direction.
The cylinders are also not identical.
You can see this by looking at your mesh:

[TomM]

In the end I want to do this:

\int_{\Omega} \left( \sigma(u) \cdot \epsilon(v) + \rho \frac{1}{\Delta t} \, u \cdot v \right) \, \text{d}x = \int_{\Omega} \rho \frac{1}{\Delta t} \, u_n \cdot v \, \text{d}x + \int_{\Gamma_{\text{right}}} f \cdot v \, \text{d}s

or static:

\int_\Omega \sigma(u) \cdot \epsilon(v) \, \text{d}x = \int_{\Gamma_{\text{right}}} f \cdot v \, \text{d}s

lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)

a = (
    inner(sigma(u), epsilon(v)) * dx +
    (1 / dt) * rho * inner(u, v) * dx +
    friction_coefficient * inner(u, v) * ds(skin)
)

L = inner(pressure_force, v) * (ds(right_face)) + (1 / dt) * rho * inner(u_n, v) * dx

but for testing I set friction_coefficient = 0 and poisson ratio nu = 0 as well as dt to infinite (and therefore simplifying the equations). Did I make a mistake when defining these equations?

The pressure on the rhs and the Youngs modulus E are in pascal.

If I understand your reply correctly I should define my mesh in a way so that is it radially symmetric and this could fix my problem?

Thank you for your time.

[dokken]

@TomM

def sigma(u):
return E * epsilon(u)

@TomM

In the end I want to do this:

\int_{\Omega} \left( \sigma(u) \cdot \epsilon(v) + \rho \frac{1}{\Delta t} , u \cdot v \right) , \text{d}x = \int_{\Omega} \rho \frac{1}{\Delta t} , u_n \cdot v , \text{d}x + \int_{\Gamma_{\text{right}}} f \cdot v , \text{d}s ∫Ω(σ(u)⋅ϵ(v)+ρ1Δtu⋅v)dx=∫Ωρ1Δtun⋅vdx+∫Γrightf⋅vds

I.e. your sigma is not correct.

Secondly, my point is that your mesh is

1. Way too coarse to represent a cylinder
2. Is not symmetric in the way you are thinking.

[TomM]

My full sigma is this:

def sigma(u):
    return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)

Is this correct and if not what is the mistake?

When changing the mesh by setting these numbers I get better results. Is this the way to go or should I try something else?

gmsh.option.setNumber("Mesh.MeshSizeMax", 0.1) 
gmsh.option.setNumber("Mesh.MeshSizeMin", 0.0001)  

Thanks for helping me with this!

[dokken]

THe above sigma looks better, and also the mesh resolution look more sensible.