Imposing torsion in a cylinder

5 months ago by

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

5 months ago by
The displacement  $u_{\theta}$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}}$umax  on the mantle  $r=R$r=R , i.e.:
$u_{\theta}\left(r\right)=\frac{r}{R}u_{\mathrm{max}}$uθ(r)=rR umax .
Now you transform this function into a Cartesian representation using
$\theta=\atan\left(\frac{y}{x}\right)$θ=atan(yx ) ,
  $u_x=-u_{\theta}\sin\theta$ux=uθsinθ  ,
   $u_y=u_{\theta}\cos\theta$uy=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.
#Boundary condition
ux = Expression(("-((pow(x[0]*x[0]+x[1]*x[1],0.5)/r))*g*l*sin(pow(tan(x[1]/x[0]),-1))"), r = Ro, g=gamma, l=L, degree = 5)
uy = Expression(("((pow(x[0]*x[0]+x[1]*x[1],0.5)/r))*g*l*cos(pow(tan(x[1]/x[0]),-1))"), r = Ro, g=gamma, l=L, degree = 5)
u0 = Constant((0.0,0.0,0.0))
bctx = DirichletBC(V.sub(0).sub(0), ux, top)
bcty = DirichletBC(V.sub(0).sub(1), uy, top)
bcb = DirichletBC(V.sub(0), u0, bottom)
bcs = [bctx,bcty,bcb]

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

written 13 days ago by GiuliettaBevi  
5 months ago by
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]​
5 months ago by
Thank you very much

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

Similar posts:
Search »