Hi,

I’m trying to find the solution of displacements (u) and Lagrange multipliers (p) in 2D rectangle body with fixed left edge (u = 0) and traction specified on the right edge. All other edges are traction-free. I have never dealt with Lagrange multipliers and I’m struggling to set up the problem in Fenics.

Here’s the Euler-Lagrnage equations I’ve found:

∇∙σ=0 (1)

∇∙u=0 (2)

σ=2μϵ+pI (3)

and by definition of strain:

ϵ=1/2 (∇u+(∇u)^T ) (4)

To solve for u and p, I did the following:

# Define function space

V1 = VectorFunctionSpace(mesh, ‘P’, 1)

V2 = FunctionSpace(mesh, ‘R’, 0)

# Assign left/right boundaries

def left(x,on_boundary):

return near(x[0],0) and on_boundary

def right(x,on_boundary):

return near(x[0],3) and on_boundary

# Mark the boundary on the right

boundary_subdomains = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)

boundary_subdomains.set_all(0)

AutoSubDomain(left).mark(boundary_subdomains, 1)

AutoSubDomain(right).mark(boundary_subdomains, 2) # This is the boundary where the traction is applied

dss = ds(subdomain_data=boundary_subdomains)

# BC: u_x = u_y = 0 at the left

bc = DirichletBC(V,Constant((0,0)),left)

# Define strain

def epsilon(u):

return 0.5*(grad(u) + grad(u).T)

# Define stress

def sigma(u,p):

return 2*mu*epsilon(u) + p*Identity(d)

# Define variational problem

u = TrialFunction(V1)

p = TrialFunction(V2)

d = u.geometric_dimension() # space dimension of u

v = TestFunction(V1)

q = TestFunction(V2)

f = Constant((0,0)) # External force = 0

T = Constant((0,0)) # Traction free BC

I’m planning to add Petrov-Galerkin term in the weak form of the equation (Why I introduced Lagrange multiplier here). Am I going the right way so far? How should I derive the weak form to solve for u and p?