How to set confining pressure with 3D cylinder in FEniCS?

1

Hi,
I am now performing deformation analysis on a three-dimensional cylindrical specimen. When performing a triaxial compression test, it is necessary to set a confining pressure on the surface of the cylinder, and the direction of the confining pressure is perpendicular to the cylinder surface inward. I don’t know how to set this force on the cylindrical surface（confine_boundary.mark(boundary_facet_function,1)）. Thank you so much for someone who can help me solve this problem.

from dolfin import *
from mshr import *
from ufl import replace
import numpy as np
import matplotlib.pyplot as plt

#Making a cylindrical geometry
#cylinder = Cylinder('coordinate of center of the top circle', 'coordinate of center of the bottom circle', 'radius of the circle at the top', 'radius of the circle at the bottom')
r1 = 0.025
r2 = 0.025
h1 = 0.0
h2 = 0.1
cylinder = Cylinder(Point(0, 0, h2), Point(0, 0, h1), r1, r2)
geometry = cylinder
# Making Mesh (40 corresponds to the mesh density)
mesh = generate_mesh(geometry, 30)
plot(mesh)

### Boundary Conditions ###
class ConfineBoundary(SubDomain):
    def inside(self, x, on_boundary):
        r = np.sqrt(x[0] * x[0] + x[1] * x[1])
        return near(r, r1) and between(x[0], (-r1, r1)) and between(x[1], (-r1, r1)) and on_boundary
confine_boundary = ConfineBoundary()
    
class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2],h1) and on_boundary
bottom_boundary = BottomBoundary()

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2],h2) and on_boundary
top_boundary = TopBoundary()

boundary_facet_function = MeshFunction('size_t', mesh, 2)
boundary_facet_function.set_all(0) # from 0 
confine_boundary.mark(boundary_facet_function,1)
bottom_boundary.mark(boundary_facet_function,2)
top_boundary.mark(boundary_facet_function,3)

ds = Measure("ds")(subdomain_data=boundary_facet_function)
2

In your variational linear form, something like

p = Constant(1.0)
n = FacetNormal(mesh)
L = dot(-p*n, u)*ds(1)
3

Thank you for your detailed answer, I will try it in the code

4

Hi，
I have changed my code with" dot(-p*n, v)*ds(1)", but the stress or strain results not change. Maybe there are some problem in my code.

p = Constant(2e7)
F = Constant((0, 0, -5e7))
n = FacetNormal(mesh)
res =  inner(sigma_new, epsilon(v))*dx - inner(F, v)*ds(3) - dot(-p*n, v)*ds(1)
a, L = lhs(res), rhs(res)