### Imposing torsion in a cylinder

239

views

0

Hi!

I'm working on a problem of non linear elasticity and I have an hollow cylinder to which I would like to impose a finite torsion on the top.

So, the code would like to be

from __future__ import print_function

from fenics import *

import mshr

import numpy as np

domain1 = mshr.Cylinder(Point(0,0,0),Point(0,0,100),10,10,100)

domain2 = mshr.Cylinder(Point(0,0,0),Point(0,0,100),5,5,50)

domain = domain1 - domain2

mesh = mshr.generate_mesh(domain,30)

Now, in order to impose that the displacement u_\Theta is equal to an angle, how can do it?

I know that I have to create a sort of map but I don't know how to do it

Thank you very much for the help

Giulia Bevilacqua

I'm working on a problem of non linear elasticity and I have an hollow cylinder to which I would like to impose a finite torsion on the top.

So, the code would like to be

from __future__ import print_function

from fenics import *

import mshr

import numpy as np

domain1 = mshr.Cylinder(Point(0,0,0),Point(0,0,100),10,10,100)

domain2 = mshr.Cylinder(Point(0,0,0),Point(0,0,100),5,5,50)

domain = domain1 - domain2

mesh = mshr.generate_mesh(domain,30)

Now, in order to impose that the displacement u_\Theta is equal to an angle, how can do it?

I know that I have to create a sort of map but I don't know how to do it

Thank you very much for the help

Giulia Bevilacqua

Community: FEniCS Project

### 3 Answers

0

The displacement $u_{\theta}$

$u_{\theta}\left(r\right)=\frac{r}{R}u_{\mathrm{max}}$

Now you transform this function into a Cartesian representation using

$r=\sqrt{x^2+y^2}$

$\theta=\atan\left(\frac{y}{x}\right)$

and

$u_x=-u_{\theta}\sin\theta$

$u_y=u_{\theta}\cos\theta$

assuming the third axis is the one you want to rotate around.

These components can then be used in the standard straightforward way of implementing DirichletBC in FEniCS.

`u`_{θ}cannot be equal to an angle. It is a displacement after all. You can prescribe a linear function in the radial coordinate $r$`r`with a maximum value of $u_{\mathrm{max}}$`u`_{max}on the mantle $r=R$`r`=`R`, i.e.:$u_{\theta}\left(r\right)=\frac{r}{R}u_{\mathrm{max}}$

`u`_{θ}(`r`)=`r``R``u`_{max}.Now you transform this function into a Cartesian representation using

$r=\sqrt{x^2+y^2}$

`r`=√`x`^{2}+`y`^{2},$\theta=\atan\left(\frac{y}{x}\right)$

`θ`=atan(`y``x`) ,and

$u_x=-u_{\theta}\sin\theta$

`u`_{x}=−`u`_{θ}sin`θ`,$u_y=u_{\theta}\cos\theta$

`u`_{y}=`u`_{θ}cos`θ`,assuming the third axis is the one you want to rotate around.

These components can then be used in the standard straightforward way of implementing DirichletBC in FEniCS.

0

i dont have answer to your exact question but to have an idea of how to apply boundary condition causing torsion you can see the following snippet. (it applies twist on a unit cube)

```
def mesh(self):
n = 8
return UnitCube(n, n, n)
def dirichlet_values(self):
clamp = Expression(("0.0", "0.0", "0.0"))
twist = Expression(("0.0",
"y0 + (x[1] - y0) * cos(theta) - (x[2] - z0) * sin(theta) - x[1]",
"z0 + (x[1] - y0) * sin(theta) + (x[2] - z0) * cos(theta) - x[2]"),
y0=0.5, z0=0.5, theta=pi/3)
return [clamp, twist]
def dirichlet_boundaries(self):
left = "x[0] == 0.0"
right = "x[0] == 1.0"
return [left, right]
```

0

Please login to add an answer/comment or follow this question.

However, if I use your advice, the code doesn't run because it did the first iteration and then I see -nan.

I think that the mistake is in the BC but I don't know what it is.

Thanks a lot

Giulia