The accompanying Jupyter notebook can be obtained here 2_load_and_boundary_conditions

Loads and boundary conditions

The preceding tutorial focused on body forces. However, in situations where you need to apply traction force or a Dirichlet Boundary Condition (DBC) to specific parts of the mesh, this tutorial will address those scenarios.

from dolfin import *

Mesh entities

Conceptually, a mesh (modeled by the class Mesh), consists of a collection of mesh entities. A mesh entity is a pair (d, i), where d is the topological dimension of the mesh entity and i is a unique index of the mesh entity. Mesh entities are numbered within each topological dimension from 0 to nd − 1, where nd is the number of mesh entities of topological dimension d.

Entity Dimension
Vertex 0
Edge 1
Face 2
Facet D-1
Cell D 
length, depth = .6, .200
num_ele_along_depth = 10
ele_size = depth/num_ele_along_depth
mesh = RectangleMesh(Point(0, 0), Point(length, depth),
                     int(length/ele_size), int(depth/ele_size))

V = FunctionSpace(mesh, 'CG', 1)
U = VectorFunctionSpace(mesh, 'CG', 1)
T0 = TensorFunctionSpace(mesh, 'DG', 0)

dim = mesh.topology().dim()

Define the boundaries (Subdomains)

clamped_boundary = CompiledSubDomain("near(x[0],0)")
load_boundary = CompiledSubDomain("near(x[1],0.2) && x[0]>0.5")

Mark the boundaries on mesh

support_tag, load_tag = 1, 2
mf = MeshFunction("size_t", mesh, 1)
mf.set_all(0)
clamped_boundary.mark(mf,support_tag)
load_boundary.mark(mf,load_tag)

Visualize the boundaries

with XDMFFile("input/subdomains.xdmf") as outfile:
    outfile.write(mf)

bc = DirichletBC(U, Constant((0,)*dim), clamped_boundary)

E, nu = 2e11, 0.3
rho, g = 7800, 9.81
\[ \begin{align} \sigma &= \lambda\,\hbox{tr}\,(\varepsilon) I + 2\mu\varepsilon,\\ \varepsilon &= \frac{1}{2}\left(\nabla u + (\nabla u)^{\top}\right), \end{align} \]
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))


def epsilon(u):
    return 0.5*(grad(u) + grad(u).T)


def sigma(u):
    return lmbda*tr(epsilon(u))*Identity(dim) + 2*mu*epsilon(u)

# Define variational problem
u, v = TrialFunction(U), TestFunction(U)

f = Constant((0, -rho*g))

ds = Measure("ds",subdomain_data=mf)

a = inner(sigma(u), epsilon(v))*dx
t = Constant((0,100))
L = dot(f, v)*dx + dot(t,v)*ds(load_tag)

# Compute solution
u = Function(U)
solve(a == L, u, bc)

plt.figure(figsize=(18, 16))
# Plot solution
scale_factor = 1e5
plot(u*scale_factor, title='Displacement', mode='displacement')

<matplotlib.collections.PolyCollection at 0x7f38d4d042b0>

png

u.vector().min()

-1.8010419648112784e-06

plt.figure(figsize=(18, 16))
# Plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(dim)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Von Mises stress')

<matplotlib.tri.tricontour.TriContourSet at 0x7f38d4c3d1d0>

png

plt.figure(figsize=(18, 16))
# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
      u_magnitude.vector().vec().array.min(),
      u_magnitude.vector().vec().array.max())

min/max u: -5.71041349380966e-13 1.8360022678476485e-06

png

stress = project(sigma(u),T0)
plt.figure(figsize=(18, 16))
plot(stress[0,0], title='$\sigma_{xx}$')

<matplotlib.tri.tricontour.TriContourSet at 0x7f38d4af0b38>

png

plt.figure(figsize=(18, 16))
plot(stress[0,1], title='$\sigma_{xy}$')

<matplotlib.tri.tricontour.TriContourSet at 0x7f38d4a8ad68>

png

plt.figure(figsize=(18, 16))
plot(stress[1,1], title='$\sigma_{yy}$')

<matplotlib.tri.tricontour.TriContourSet at 0x7f38d4a100f0>

png