Fixing dofs normal to direction vector used in traction boundary condition for linear elasticity

General
1

Hi,

I would like to partially fix the dofs where a traction force is applied on a circle on the surface of a body. Right now the way I am applying the traction is the following:

# Traction direction: push along +x (change sign if needed)
T = fem.Constant(domain, -p_traction * default_scalar_type((node_normal[0], node_normal[1], node_normal[2])))

a = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx
L = ufl.dot(f_vec, v) * ufl.dx + ufl.dot(T, v) * ds(2)

Is there a way of enforcing zero displacement normal to the node_normal vector for the nodes where the traction is applied?

Running code and a vtk mesh can be found in another post:

Thank you,
Alex

2

As you are using Lagrange elements, you can’t enforce this through a normal dirichlet condition if your normals do not align with the cartesian grid. You can however use DOLFINx MPC for it, as shown in: Stokes flow with slip conditions — DOLFINx-MPC: An extension to DOLFINx for multi point constraints in the setting of Stokes equation.
You could alternatively try to enforce it weakly with Nitsche’s method.

3

Thank you for the prompt reply. Could you give me an example of how to do it for a coordinate direction for the current code?

Do you have a post reference discussing how Nitsche’s method is applied to a similar problem in Fenicsx?

Thank you

4

Then you can take the DirichletBC in the appropriate sub space. For instance, if you want the component in the ith direction to be 0, you set a bc on V.sub(i).
In Mixed finite element problems — FEniCS Workshop there is a fairly comprehensive explaination.

Nitsche methods is for instance covered in: Weak imposition of Dirichlet conditions for the Poisson problem — FEniCSx tutorial
and asimov-contact/python/dolfinx_contact/one_sided/nitsche_ufl.py at main · Wells-Group/asimov-contact · GitHub goes through linear elasticity with nitsche/contact conditions, I.e. dot(u,n)-g=0

5

Thank you for the reply. I added this at the end of the build_primal_problem() function:

# Dirichlet BC to stop deflection normal to x (lock_num == 1) or y (lock_num == 2) axis   
if lock_num == 1:
    right_facets = facet_tags.find(2)
    right_dofs   = fem.locate_dofs_topological(V, fdim, right_facets)
    bc_1 = fem.dirichletbc(0.0, right_dofs, V.sub(1))
    bc_2 = fem.dirichletbc(0.0, right_dofs, V.sub(2))
    return V, a, L, [bc_left, bc_1, bc_2]
if lock_num == 2:
    right_facets = facet_tags.find(2)
    right_dofs   = fem.locate_dofs_topological(V, fdim, right_facets)
    bc_1 = fem.dirichletbc(0.0, right_dofs, V.sub(0))
    bc_2 = fem.dirichletbc(0.0, right_dofs, V.sub(2))
    return V, a, L, [bc_left, bc_1, bc_2]

return V, a, L, [bc_left]

Is this correct? From the output displacements it seems to be working

6

No, to be completely safe here, you should use:

# Dirichlet BC to stop deflection normal to x (lock_num == 1) or y (lock_num == 2) axis   
if lock_num == 1:
    right_facets = facet_tags.find(2)
    right_dofs1   = fem.locate_dofs_topological(V.sub(1), fdim, right_facets)
    bc_1 = fem.dirichletbc(0.0, right_dofs1, V.sub(1))
    right_dofs2   = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
    bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
    return V, a, L, [bc_left, bc_1, bc_2]
if lock_num == 2:
    right_facets = facet_tags.find(2)
    right_dofs0   = fem.locate_dofs_topological(V.sub(0), fdim, right_facets)
    bc_1 = fem.dirichletbc(0.0, right_dofs0, V.sub(0))
    right_dofs2   = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
    bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
    return V, a, L, [bc_left, bc_1, bc_2]

return V, a, L, [bc_left]
7

That works thank you! The previous case was deflecting less but it was still moving upwards.