### Application of Traction BC

153

views

0

Hello everyone,

I'm trying to understand the application of traction boundary condition on specified face. For example, I want to fixed displacement on left side of the rectangular box and I want to apply vertical traction on right face of the box. For that purpose, I prepared a sample example on elasiticty code. But it seems like I am doing something wrong. Do you have any suggestion ? Thanks.

I'm trying to understand the application of traction boundary condition on specified face. For example, I want to fixed displacement on left side of the rectangular box and I want to apply vertical traction on right face of the box. For that purpose, I prepared a sample example on elasiticty code. But it seems like I am doing something wrong. Do you have any suggestion ? Thanks.

```
from __future__ import print_function
from fenics import *
%matplotlib inline
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 2)
# 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)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and abs(x[0] - L) < tol
right_boundary = RightBoundary()
boundary_part = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_part.set_all(0)
right_boundary.mark(boundary_part, 1)
ds = Measure("ds")[boundary_part]
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
u = Function(V)
d = 3 # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 0, 50))
Psi=(lambda_/2)*(tr(epsilon(u)))**2+mu*tr(epsilon(u)*epsilon(u))
Pi=Psi*dx-dot(f, u)*dx - dot(T, u)*ds(1)
F=derivative(Pi,u,v)
solve(F==0,u,bc)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
u_magnitude.vector().array().min(),
u_magnitude.vector().array().max())
```

Community: FEniCS Project

### 3 Answers

4

You mean this?

```
from __future__ import print_function
from fenics import *
L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'CG', 2)
# Define boundary condition
tol = 1E-14
d = 3 # space dimension
class ClampBound(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and on_boundary
class RightBound(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], L)
clamped_bc = ClampBound()
right_boundary = RightBound()
boundary_part = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_part.set_all(0)
clamped_bc.mark(boundary_part, 2)
right_boundary.mark(boundary_part, 3)
ds = Measure("ds")[boundary_part]
# plot bcs
ds_mark = File("ds_mark.pvd")
ds_mark << boundary_part
# boundary
bc = DirichletBC(V, Constant((0, 0, 0)), clamped_bc)
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
u = Function(V)
v = TestFunction(V)
T = Constant((0.05, 0., 0.))
Psi = (lambda_/2)*(tr(epsilon(u)))**2+mu*tr(epsilon(u)*epsilon(u))
Pi = Psi*dx - inner(T, u)*ds(3)
F = derivative(Pi, u, v)
solve(F == 0, u, bc)
# u_magnitude = sqrt(dot(u, u))
# u_magnitude = project(u_magnitude, V)
disp = File("disp.pvd")
disp << u
# plot(u_magnitude, 'Displacement magnitude')
# print('min/max u:',
# u_magnitude.vector().array().min(),
# u_magnitude.vector().array().max())
```

Yes exactly, thanks a lot.

written
4 months ago by
Christian

0

Hello,

For displacement traction you can select the dof of your bc model and apply the displacement, like this:

For displacement traction you can select the dof of your bc model and apply the displacement, like this:

```
class HorizBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and x[1] == W
hz_boundary = HorizBoundary()
bc_disp = DirichletBC(V.sub(1), Constant(2.5e-2), hz_boundary)
list_bcs = [bc, bc_disp]
....
solve(F==0, u, list_bcs)
```

0

Hello Ricardo,

Actually I am trying to apply force on the right face as a traction not displacement.

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