How to define non-uniform traction boundary conditions at nodes

General
1

Hello everyone，as shown, for the elasticity problem in a 2D area (left side fixed, horizontal traction applied on the right side), if the traction force is constant, the code is as follows:

from fenics import *

mesh = RectangleMesh(Point(0., 0.), Point(5, 5), 10, 10)
V = VectorFunctionSpace(mesh, "P", 1)

def left(x, on_boundary):
    return near(x[0], 0) and on_boundary
boundary_mark = MeshFunction("size_t", mesh, 1)
class Rightboundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and abs(x[0]-5) < tol

right = Rightboundary()
right.mark(boundary_mark, 1)
ds = Measure('ds', domain=mesh, subdomain_data = boundary_mark)

E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

bc = DirichletBC(V, Constant((0, 0)), left)
T = 100
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_)) * dx
l = dot(Constant((T, 0)), u_) * ds(1)

u = Function(V, name="Displacement")
solve(a == l, u, bc)
File("result/U.pvd") << u

If the traction force at each node is different (assuming from bottom to top, they are (T0, 0), (T1, 0), (T2, 0)…), how should it be handled? I hope to receive some help, thank you.
faf490e84025b0b650b58b27e70e261

2

I assume that the traction should be continuous between the nodes, maybe linear? If so, create a first order Lagrange space, and assign the traction to the relevant node, using V.tabulate_dof_coordinates() and a distance computation from your input points as shown in

3

In fact, I’ve already made a similar attempt before:

from dolfin import *
import numpy as np

# Create mesh and function space
mesh = RectangleMesh(Point(0., 0.), Point(5, 5), 5, 5)
V = VectorFunctionSpace(mesh, "CG", 1)

T_data = [0, 50, 100, 120, 150, 0]
T_coordinates = [[5, 0], [5, 1], [5, 2], [5, 3], [5, 4], [5, 5]]

class Rightboundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 5) and on_boundary

right = Rightboundary()
vertex_mark = MeshFunction("size_t", mesh, 0, 0)
right.mark(vertex_mark, 1)
vertices = vertex_mark.where_equal(1)
vertex_to_dof_map = vertex_to_dof_map(V)

T_load = Function(V)

# Get all boundary coordinates
dof_coords = V.tabulate_dof_coordinates()
boundary_dofs = [vertex_to_dof_map[vertex] for vertex in vertices]
boundary_coords = np.array([dof_coords[boundary_dof] for boundary_dof in boundary_dofs])

points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(T_coordinates, 0)
distances = np.sum(np.square(points_A - points_B), axis=2)
is_close = distances < 1e2 * DOLFIN_EPS
positions = np.nonzero(is_close)
for row, col in zip(*positions):
    T_load.vector()[boundary_dofs[row]] = T_data[col]
    print(row, col, boundary_coords[row], T_coordinates[col], T_data[col])

E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

boundary_mark = MeshFunction("size_t", mesh,1)
right.mark(boundary_mark, 1)
ds = Measure('ds', domain=mesh, subdomain_data = boundary_mark)

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_)) * dx
l = dot(T_load, u_) * ds(1)

u = Function(V, name="Displacement")
solve(a == l, u)

However, the output results were missing some parts：

1 0 [5. 0.] [5, 0] 0
3 1 [5. 1.] [5, 1] 50
5 2 [5. 2.] [5, 2] 100
Solving linear variational problem.

On the other hand, as far as I know, traction should be defined in the form of a 2D vector because uniform constant traction is defined as T =Constant((0，0)). After I modified the form of the traction:
T_data = [[0,0], [50,0], [100,0], [120,0], [150,0], [0,0]]
the following error occurred:

Traceback (most recent call last):
  File "/home/dyfluid/Desktop/contact2D/node.py", line 34, in <module>
    T_load.vector()[boundary_dofs[row]] = T_data[col]
    ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
IndexError: Indices to set must be a 1D array

I hope you can provide some help, thank you.

4

Make T_load in a scalar space, i.e.
T_space = FunctionSpace(mesh, "Lagrange", 1),
T_load = Function(T_space), then in turn assign the T_data.

Then, instead of T_load in the variational form, use:

Traction = ufl.as_vector((T_load,0))
l = dot(Traction, u_)*ds(1)
5

Thank you for your guidance, my issue has been resolved. I have a question regarding the definition of traction: in the FEniCS tutorial, T is defined as the force along the boundary normal. So, for a vertical boundary (x=0), does the traction condition have to be set horizontally (T = Constant(a, 0))? If it is set as T = Constant(a, b), what would be the result? I hope you can answer my question, thank you.

1715326235229
1715326235229876×241 37.8 KB

6

\sigma \cdot n should be zero in the tangential direction if you are only applying normal stresses.
See for instance page 3 of:

In particular:

Stress at a point x can be mathematically defined as limA→{x} tA/A, where tA repre-
sents the internal force applied across the cutting surface A, which encloses x and shrinks
to the point x. The force tA is a vector and does not necessarily have to be orthogonal to
the surface. The component of this force that is normal to the surface is known as the nor-
mal stress (σn), and the component tangential to the surface is termed the shearing stress
(τn), with the subscript n indicating a normal vector to the cutting surface.

and Page 6:

  1. Neumann boundary conditions (also known as traction boundary conditions).
    These conditions specify the traction vector t (force per unit area) on the Neumann bound-
    ary ΓN . Mathematically: σ · n = t on ΓN , where σ is the stress tensor, n is the unit
    outward normal vector to the boundary, and t is a given function. Physically, this means
    that the boundaryΓN is subjected to a specified force or traction
7

Sure, I understand now. Thank you for your patient explanation!