Application of Traction BC


153
views
0
4 months ago by
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.

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
4 months ago by
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
4 months ago by
Hello,

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
4 months ago by

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.

Similar posts:
Search »