Mixing Neuman and Dirichlet BC on cylinder (linear elasticity)
5 months ago by
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
if near(x,cylinder_length, tol):
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 < 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
5 months ago by
from dolfin import * class Top(SubDomain): def inside(self, x, on_boundary): return x > L - DOLFIN_EPS # You also could use near(...) . . . # Marker m = MeshFunction('size_t', mesh, mesh.geometry().dim()-1) m.set_all(0) 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)
5 months ago by
It will be much clear to know what's the problem.
Please login to add an answer/comment or follow this question.