Free slip boundary condition in elasticity

General
1

I would like to have a mesh representing an elastic material to be able to “slide” or slip freely along a boundary. Eg., in this example:

from dolfin import *


class AllBoundaries(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < -49


def solve_problem(mesh, mfunc, force=[]):
    V = VectorFunctionSpace(mesh, "CG", 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    displacement = Function(V)

    bc = [DirichletBC(V, Constant((0, 0)), mfunc, 1)]

    F = Constant(force)
    E = Constant(5000)
    nu = Constant(0.3)
    mu = E / (2.0 * (1 + nu))
    lmbda = E * nu / ((1.0 + nu) * (1 - 2 * nu))
    sigma = 2.0 * mu * sym(grad(u)) + lmbda * tr(sym(grad(u))) * Identity(2)
    solve(inner(sigma, grad(v)) * dx == inner(F, v) * dx, displacement, bc)
    displacement.set_allow_extrapolation(True)
    return displacement


#################################################################################
mesh = RectangleMesh(Point(-50.0, -50.0), Point(50.0, 50), 20, 20)

mfunc = MeshFunction("size_t", mesh, 1, mesh.domains())
mfunc.set_all(0)
allb = AllBoundaries()
allb.mark(mfunc, 1)

F = [20, 10]
displacement = solve_problem(mesh, mfunc, F)
ALE.move(mesh, displacement)

### PLOT
import vedo
varrow = vedo.Arrow2D([0, 0], F).z(1).c("red4")
vmesh = vedo.Mesh([mesh.coordinates(), mesh.cells()]).c("k4").lc("k5")
vline = vedo.Line([(-50, -50), (-50, 50)]).c("red4").linewidth(5)
vedo.show(vmesh, vline, varrow, axes=1)

Screenshot from 2023-09-28 20-25-33
Screenshot from 2023-09-28 20-25-33541×596 93.7 KB

I would expect to see the mesh to “slide” a bit up the red boundary line.
How can I achieve that, if that’s possible at all?
Thanks!

2

If the slip boundary aligns with the coordinate system, you can apply a bc on the x component of the displacement (set it to 0), and let y be free.

3

If I try this I still get the same solution…:

from dolfin import *

class BoundarySource(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = 0
        # values[1] = 0 # free slip?

    def value_shape(self):
        return (2,)


def solve_problem(mesh, force=()):
    V = VectorFunctionSpace(mesh, "CG", 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    displacement = Function(V)

    subd = CompiledSubDomain("on_boundary && x[0] < -49")
    bnd_src = BoundarySource(degree=1)
    dbc = DirichletBC(V, bnd_src, subd)

    F = Constant(force)
    E = Constant(5000)
    nu = Constant(0.3)
    mu = E / (2 * (1 + nu))
    lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
    sigma = 2 * mu * sym(grad(u)) + lmbda * tr(sym(grad(u))) * Identity(2)
    solve(inner(sigma, grad(v)) * dx == inner(F, v) * dx, displacement, [dbc])
    return displacement


#############################################################################
mesh = RectangleMesh(Point(-50, -50), Point(50, 50), 20, 20)

F = [20, 10]
displacement = solve_problem(mesh, F)
ALE.move(mesh, displacement)

import vedo
varrow = vedo.Arrow2D([0, 0], F).z(1).c("red4")
vmesh = vedo.Mesh([mesh.coordinates(), mesh.cells()]).c("k4").lc("k5")
vline = vedo.Line([(-50, -50), (-50, 50)]).c("red4").linewidth(5)
vedo.show(vmesh, vline, varrow, axes=1)

what am I doing wrong?

Also, how can I generalize it to a boundary which is not aligned along a cartesian axis?

thanks for the help!

4

This is not how to set a sub condition, see: Help a starter! - #9 by dokken

5

Sorry - I tried the linked example but i’m still getting something wrong:

from dolfin import *

def clamped_boundary(x, on_boundary):
    return on_boundary and near(x[0], -50)

def solve_problem(mesh, force=()):
    V = VectorFunctionSpace(mesh, "CG", 1)
    u_ = TrialFunction(V)
    v = TestFunction(V)
    u = Function(V)

    # bc = DirichletBC(V, Constant([0,0]), clamped_boundary)
    bc = DirichletBC(V.sub(1), Constant(0), clamped_boundary)

    F = Constant(force)
    E = Constant(5000)
    nu = Constant(0.3)
    mu = E / (2 * (1 + nu))
    lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
    sigma = 2*mu * sym(grad(u_)) + lmbda * tr(sym(grad(u_))) * Identity(2)
    solve(inner(sigma, grad(v)) * dx == inner(F, v) * dx, u, [bc])
    return u

#############################################################################
mesh = RectangleMesh(Point(-50, -50), Point(50, 50), 20, 20)

F = [20, 10]
u = solve_problem(mesh, F)
ALE.move(mesh, u)

import vedo
varrow = vedo.Arrow2D([0, 0], F).z(1).c("red4")
vmesh = vedo.Mesh([mesh.coordinates(), mesh.cells()]).c("k4").lc("k5")
vline = vedo.Line([(-50, -50), (-50, 50)]).c("red4").linewidth(5)
vedo.show(vmesh, vline, varrow, axes=8)

Screenshot from 2023-09-29 17-43-03
Screenshot from 2023-09-29 17-43-03714×259 12.4 KB

6

Ouch! that was me being stupid… yes it works with:

    bc = DirichletBC(V.sub(0), Constant(0), clamped_boundary)
...
F = [20, 0]

Screenshot from 2023-09-29 18-13-20

Question is now if one could generalize it for boundaries which are not aligned in x or y.

7

Then you have to either:

  1. weakly enforce it with Nitsche/DG type methods
  2. use dolfinx MPC https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_stokes.py
1 Like