Applying stress constraints to a point in two-dimensional geometry

Hello everyone, as shown in the diagram, I want to perform mechanical calculations on a two-dimensional circular region. How can I add stress constraints at the center of the circle in FEniCS: radial stress equals circumferential stress?(The red point represents the center of the circle)
3f0cc804d5306117a92a945fc5e4c8e

It’s hard to give some guidance without a set of equations and a MWE, but a general suggestion could be: if you are willing to impose a constraint, use a Lagrange multiplier. Since your constraint is pointwise, the correct multiplier will be in a finite element space that is known in legacy FEniCS as "Real"

I want to calculate the thermal expansion displacement of a circular region at 100°C. Here are my code and geo file. Could you provide specific guidance based on this? I want to apply a constraint at the center of the circle to make the radial stress equal to the circumferential stress

from dolfin import *

# mesh and function space
mesh = Mesh("single.xml")
boundary = MeshFunction("size_t", mesh, "single_facet_region.xml")
subdomains = MeshFunction("size_t", mesh, "single_physical_region.xml")
Vu = VectorFunctionSpace(mesh, 'CG', 2)

# parameters
E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = Constant(1e-5)

def eps(v):
    return sym(grad(v))
def sigma(v, dT):
    return (lmbda*tr(eps(v))- alpha*(3*lmbda+2*mu)*dT)*Identity(2) + 2.0*mu*eps(v)

# variational problem
du = TrialFunction(Vu)
u_ = TestFunction(Vu)
Wint = inner(sigma(du, 100), eps(u_))*dx
aM = lhs(Wint)
LM = rhs(Wint)

# Compute solution
u = Function(Vu, name="Displacement")
solve(aM == LM, u)
File("result/circle_u.pvd") << u

mesh:single.geo

// Gmsh project created on Tue Apr 09 10:25:53 2024
SetFactory("OpenCASCADE");

N = 50;
//+
Circle(1) = {0, 0, 0, 1, 0, 2*Pi};
//+
Transfinite Curve {1} = N Using Progression 1;
//+
Curve Loop(1) = {1};
//+
Plane Surface(1) = {1};
//+
Physical Curve("f", 2) = {1};
//+
Physical Surface("f", 3) = {1};

Unfortunately the idea I had in mind doesn’t work as shown eg in Assemble only on the vertices because a required feature is not implemented. I’ll let someone chime in in case there are other ideas or approaches.

Thank you for your reply!

Hi, what is the reason for which you would like to apply such a constraint ? You could achieve this with a Lagrange multiplier and a MixedElement formulation with the constraint acting only in the central region. But this seems weird to me in terms of physics…