Stress continuity using mixed-dimensional branch

2

Hi Hernan,

To impose the stress continuity condition, I guess you could use a second Lagrange multiplier, as you did for the displacement continuity condition.
You would then have 4 spaces in your MixedFunctionSpace :

V_I_disp  = FunctionSpace(mesh_I,FE)   # Lagrange multiplier - displacement continuity
V_I_stress  = FunctionSpace(mesh_I,FE)   # Lagrange multiplier -stress continuity
W = MixedFunctionSpace(V_L, V_R, V_I_disp, V_I_stress)

# Function and test functions
(ul, ur, l_disp, l_stress) = TrialFunctions(W)
(wl, wr, e_disp, e_stress) = TestFunctions(W)

Regarding your variational formulation, I suspect a few mistakes :
(1) in a_left : the condition on the neumann boundary \sigma(u_l)n_l = g ~\text{in}~ \Gamma_N (note : I guess you meant u_l here since u_r is not defined on that boundary) uses ds_L(3) where 3 is the marker for the interface. I think it should be ds_L(2) instead (and it has to match with the term involving g).
(2) I am not sure we need the term inner(sigma(ul)*FacetNormal(mesh_R), wr)*ds_R(3) in a_right.
(3) To be coherent with the term inner(ur-ul, e_disp)*ds_I (1st Lagrange multiplier), shouldn’t we have inner(l_disp, wr-wl)*ds_I ?

Imposing the stress continuity condition with the second Lagrange multiplier, the variational formulation should look like :

# Variational problem
a_left = inner(sigma(ul),sym(grad(wl)))*dx_L + inner(sigma(ul)*FacetNormal(mesh_L), wl)*ds_L(2)
L_left = inner(g, wl)*ds_L(2)

a_right = inner(sigma(ur),sym(grad(wr)))*dx_R # - inner(sigma(ul)*FacetNormal(mesh_R), wr)*ds_R(3)
L_right = inner(Constant((0,0)), wr)*dx_R

nL = Constant((1,0))
nR = Constant((-1,0))

a = a_left + a_right \
  + inner(l_disp, wr-wl)*ds_I + inner(ur-ul, e_disp)*ds_I \
  + inner(l_stress, sigma(wl)*nL + sigma(wr)*nR)*ds_I + inner(sigma(ul)*nL + sigma(ur)*nR, e_stress)*ds_I
L = L_left + L_right

Note : FacetNormal cannot be used directly in cell type integrals (as ds_I which is of type dx). Here, the domain is simple and we know the expression of the normal at the interface. When it is not the case, you can define nL as the projection of FacetNormal(mesh_L) on V_L (same for nR).

Hope that helps!

4 Likes