Mixing Neuman and Dirichlet BC on cylinder (linear elasticity)

5 months ago by
Dear FEniCS community,

I am quite new to fenics but I already enjoy it. I've been trying to adapt the linear elasticity problem that can be found in the fenics tutorial. Here is what I want to do:

1) I use a cylinder geometry : domain = Cylinder(Point(cylinder_length, 0, 0), Point(0, 0, 0), W, W, 16)
2) As in the tutorial I applied the Dirichlet boundary conditions on the x = 0 face of the cylinder
3) Now, I would like to apply a Neuman boundary condition on the rest of the boundary such that :
        i) For the boundary corresponding to the face x = cylinder_length, there is a constant Neuman boundary condition (pressure) T = (0,F,0) (Note that T should be a vector)
        ii) For the remaining part of the boundary (boundary of the cylinder minus the face at x=0 minus the face at x=cylinder_length), the Neuman boundary condition (pressure) is T = (0,0,0)

My question is how do I write an expression for point 3) to impose this Neuman condition ? I had a look on how to write Expressions in fenics but its not clear to me how to achieve this.

Note that everything else should be correct. It's mainly a copy paste from the tutorial. For example my functional is :
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

Thank you in advance for your time !!!

In python I would write an expression like :
tol = 1e-14
def MyExpression(x[0],x[1],x[3]):
      if near(x[0],cylinder_length, tol):
             return [0,F,0]
             return [0,0,0]

Can I do something similar ?

My code so far

#!/usr/bin/env python

from __future__ import print_function
from pylab import *
from fenics import *
from mshr import *

# Scaled variables
cylinder_length = 1.; W = 0.2
mu = 1
rho = 1
delta = W/cylinder_length
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
#domain = Circle(Point(0,0),1)
domain = Cylinder(Point(cylinder_length, 0, 0), Point(0, 0, 0), W, W, 16)
mesh = generate_mesh(domain, 16)
#mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 20, 6, 6)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))

#My Question : how do I write an Expression for T? Can I write it as a python function ?

a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)
Community: FEniCS Project

3 Answers

5 months ago by
The best way to do it is marking the top boundary of your geometry and then impose the Neumann BC. Try something like this:

from dolfin import *

class Top(SubDomain):
  def inside(self, x, on_boundary):
    return x[0] > L - DOLFIN_EPS  # You also could use near(...)


# Marker
m = MeshFunction('size_t', mesh, mesh.geometry().dim()-1)
Top().mark(m, 1)

# Measure
ds = Measure('ds', domain=mesh, subdomain_data=m)


# Variatonal formulation
F = 1.0
T = Constant((F, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds(1)
Thank you very much for your answer.
I got no problem to get my code working with your suggestion. However, I am getting strange results when I let
m = MeshFunction('size_t', mesh, mesh.geometry().dim()-1)​

Is it possible that I should change the dimension to

m = MeshFunction('size_t', mesh, mesh.geometry().dim()-2)

or dim()-3. Like this I have a solution that looks much more physical.

However as I said earlier, I am just starting with fenics so I should learn how the functions appearing in the code you sent me works.

Also, so far I wrote the functional (as in the tutorial) as

L = dot(f, v)*dx + dot(T, v)*ds(1)

but SICHENG SUN suggests

L = dot(f, v)*dx + inner(T, v)*ds(1)

Can my problem comes from that?

Anyway, that was very helpful ! Thanks a lot

PS : a small detail: as SICHENG SUN it should be

T = Constant((0.0, F, 0.0))

to be consistent with my question. But it was easy to solve !

written 5 months ago by Mike  
With the argument mesh.geometry().dim()-1 you are telling to dolfin that the MeshFunction is defined on entities of dimension 2 (facets), so it should be fine. Maybe the problem is that m is not initialized with zeros. Try:

m = MeshFunction('size_t', mesh, mesh.geometry().dim()-1)​
Top().mark(m, 1)​

To your second question, both definitions of L are equivalent (your problem is not there).

Finally, I think that you have a little issue with the definition of your axes (see your uploaded figure and check their orientation)
written 5 months ago by Hernán Mella  
Thank you again for your help and for the explanations. That's very useful.
Actually, everything works even without the line
I think I have no issue with my axes. There are two situations:
1) setting, as you suggest, T = Constant((F, 0, 0)) is the correct expression for studying a force that acts along the axis to compress (F<0) or dilate (F>0) the beam (cylinder)
2) setting, as I did, T = Constant((0, F, 0)) is the correct expression for studying a force that acts on one side of the beam and is perpendicular to the axis of the cylinder. The other end of the cylinder is fixed to a wall by the Dirichlet BCs u=0. In this case, we will have a deflection of the beam. This is this problem that I am interested in.

I checked with my code that both situations are working.

The problem that I mentioned in my previous comment was coming from the fact that setting T =  (0,F=1,0) leads to a force that is so strong that it destroys the material (the mesh was completely deformed). Then I tried with T=(0,F=1e-2,0), and I obtained a beautiful result.

Note that for dilatation (situation 1) setting T = (F=1,0,0) works perfectly as the material is much more resistant to dilatation that to deflexion.

In a nutshell, that's you so much for your very helpful suggestions !

PS : result of the beam deflexion along y for T = (0,F=1e-2,0)

written 5 months ago by Mike  
5 months ago by
Could you please upload a figure to show your geometry and boundary condition? You can do it simply with Powerpoint, Paint or any other CAD software like Rhino, auto cad.
It will be much clear to know what's the problem.

Thank you very much for the answer. I updated my post with a scheme of the geometry and I also added the code that I have written so far. My question is how to write the expression for T that describes my BCs. Thank you in advance for your help.

written 5 months ago by Mike  
You do not need to define your Neuman condition where T = 0. So the code given by Mella should work. If it doesn't, try:
T = Constant((0.0, F, 0.0))

L = inner(T, v)*ds(1)​

BTW, you can also use other software like GMSH to define your geometry. It is useful particular with complicated geometry shape, and it is very friendly to the beginer.
written 5 months ago by SICHENG SUN  
I have edited the definition of T in my answer. Thank you!
written 5 months ago by Hernán Mella  

Thank you so much for your answer. Your answer combined with the one of Hernan were very helpful. I am still getting a small problem that I describe in my answer to Hernan. Other than that my code runs without errors.
Thanks for suggesting to have a look at GMSH. I was indeed wondering how to generate complicated meshes with mshr. I will have a look into it and also learn how to import the generated mesh into fenics.

written 5 months ago by Mike  
4 months ago by
maybe, have a look at my generalised boundary condition:
Please login to add an answer/comment or follow this question.

Similar posts:
Search »