Apply a normal or tangential force on a surface

dolfinx
1

Hello,

As things stand, I am applying a constant body and tension force on an imported mesh of a cube (not dependant of the shape of the mesh). Is it possible to apply a normal force or a tangential force depending on each node of the surface of application? I know it is possible to use Expression or Function but I don’t know how.

Thanks a lot for your help !

B = Constant(mesh, PETSc.ScalarType((0, 0, 0)))
T = Constant(mesh, PETSc.ScalarType((100, 0, 0)))


# Stored strain energy density (compressible neo-Hookean model)
psi = ufl.variable((mu / 2.0) * (Ic - 2) - mu * ufl.ln(J) + (lmbda / 2.0) * (ufl.ln(J))**2)

# dx = ufl.Measure("dx")

# Pi = psi*dx - ufl.inner(b, u)*dx
# First Piola-Kirchhoff Stress
P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=mesh, subdomain_data=ft, metadata=metadata)
dx = ufl.Measure("dx", domain=mesh, subdomain_data = ct, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P)*dx(10) - ufl.inner(v, T)*ds(14) - ufl.inner(v, B)*dx(10)
2

You can use the ufl.FacetNormal of mesh to create normal stresses, and similarly use a tangent projection of the normal to get a tangent vector, see Tangential gradient operator error: invalid ranks 1 and 1 in product - #2 by dokken

3

Thanks @dokken for your quick answer.
I succesfully applied a normal stress with a scalar value of 100 on each facet :

# F = ufl.inner(ufl.grad(v), P)*dx(5) - ufl.inner(100*ufl.FacetNormal(mesh), v)*ds(3)

I followed the link that you provided and tried to use the function:

def tangential_proj(u, n):
    return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u

I replaced the second inner product in the expression of F :

F = ufl.inner(ufl.grad(v), P)*dx(5) - tangential_proj(ufl.inner(ufl.FacetNormal(mesh), v), ufl.FacetNormal(mesh))*ds(3)

But it doesn’t seem to work :

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In [139], line 9
      6     return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u
      7 # Define form F (we want to find u such that F(u) = 0)
      8 # F = ufl.inner(ufl.grad(v), P)*dx(5) - ufl.inner(100*ufl.FacetNormal(mesh), v)*ds(3)
----> 9 F = ufl.inner(ufl.grad(v), P)*dx(5) - tangential_proj(ufl.inner(ufl.FacetNormal(mesh), v), ufl.FacetNormal(mesh))*ds(3)
     11 # Pi = psi*dx(7) - ufl.inner(u, T)*ds(5)
     12 
     13 # - ufl.inner(u, B)*dx(10)
     14 # F_res = ufl.derivative(Pi, u, v)
     15 # Jac = ufl.derivative(F_res, u)
     17 problem1 = dolfinx.fem.petsc.NonlinearProblem(F, u, bcs)

Cell In [139], line 6, in tangential_proj(u, n)
      5 def tangential_proj(u, n):
----> 6     return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u

IndexError: tuple index out of range

Can you help or provide an example on the use of this function ?

4

You need to provide a vector-valued expression as the first argument to tangential_proj, not a scalar.

I think perhaps you are misunderstanding the purpose of tangential_proj. To create a normal force, you can simply multiply a scalar value by the FacetNormal, since the normal direction is uniquely defined. However, you cannot create a tangent force in the same way, i.e. by multiplying the scalar value of the force by the “facet tangent”, since the tangent direction is not unique. Instead, you must define a force vector f and project it into the tangent space of the facet, i.e. remove the component of f that lies in the normal direction. tangential_proj(u, n) accomplishes exactly this, by projecting the vector u into the tangent space, i.e. the plane normal to n. For example:

f = dolfinx.fem.Constant(mesh, PETSc.ScalarType((1.0, 2.0, 3.0)))
F = ufl.inner(ufl.grad(v), P)*dx(5) - ufl.inner(tangential_proj(f, ufl.FacetNormal(mesh)), v)*ds(3)
3 Likes
5

Thank you a lot @conpierce8 for your answer ! It was very instructive and I indeed misunderstood how to use the tangential projection.

I can then use it to extract the normal projection of a vector:

def normal_proj(u,n):
    return u - tangential_proj(u,n)

Is it however possible to have more control on the value of the tangential vector, i.e. define directly its components in the same way than with the Constant function?

6

A normal projection can be achieved by simply calling

def projection(u, n):
    return ufl.inner(u, n) * n

n = FacetNormal(mesh)
vec = ufl.as_vector([1,2,3])
vec_n = projection(vec, n)

If you know the tangent of your surface, you can easily compute the projection of any vector into the tangent plane, as shown above with the normal vector. In general a tangent on a facet (and a facet normal) is not a constant value.

3 Likes