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

**URL:** <https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306>\
**Category:** General\
**Created:** [November 15, 2024, 12:51pm UTC](https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306 "2024-11-15T12:51:18Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![TomM](https://avatars.discourse-cdn.com/v4/letter/t/c0e974/32.png) [@TomM](https://fenicsproject.discourse.group/u/TomM)\
**Post date:** [November 15, 2024, 12:51pm UTC](https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306/1 "2024-11-15T12:51:18Z")

</div>

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?

 ![cylinder](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/ec56f99056335ebcad7603641b7cd806d81ecfdc.png)  
If I remove the second cylinder it looks like this (as expected)  
 ![single_cylinder](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/b7f09deaf8fbe7f037dd13fe587caa1b6354bb1c.png)  
If I move them further away from each other this effect increases as well as the elongation only in x direction  
 ![cylinder_extrem](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d2e4c2941026d31a61bcddda4c4e8da505b04ddb.png)  
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.  
 ![cylinder_z](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/752bf3353ac9d16349b835ccaa1b3be2cd1cbd15.png)

This is my code to reproduce this:

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 17, 2024, 1:21pm UTC](https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306/2 "2024-11-17T13:21:54Z")

</div>

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:

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/8e975f6a8071711a261a4fb24cdc6b6185beecee.png)

---

<div class="post-metadata">

**Author:** ![TomM](https://avatars.discourse-cdn.com/v4/letter/t/c0e974/32.png) [@TomM](https://fenicsproject.discourse.group/u/TomM)\
**Post date:** [November 18, 2024, 9:28am UTC](https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306/3 "2024-11-18T09:28:37Z")

</div>

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

```auto
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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 18, 2024, 11:57am UTC](https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306/4 "2024-11-18T11:57:59Z")

</div>

> [@TomM](#):
>
> ```auto
> 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.

---

<div class="post-metadata">

**Author:** ![TomM](https://avatars.discourse-cdn.com/v4/letter/t/c0e974/32.png) [@TomM](https://fenicsproject.discourse.group/u/TomM)\
**Post date:** [November 18, 2024, 12:51pm UTC](https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306/5 "2024-11-18T12:51:45Z")

</div>

My full sigma is this:

```auto
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?

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

```

 ![better_mesh](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/7ac219ffaa4fa6d867a099a65fa80a6210807db8.png)  
Thanks for helping me with this!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 18, 2024, 12:57pm UTC](https://fenicsproject.discourse.group/t/why-do-non-interacting-cylinders-bend-in-unintuitive-direction-under-identical-forces/16306/6 "2024-11-18T12:57:57Z")

</div>

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