Implementing elastoplasticty coupled to pore pressure with mgis.fenics

variational formulation
1

Hello,

I would like to simulate elastoplasticity coupled to pore pressure in 2D. I am implementing the Mohr-Coulomb constitutive model and so I intend to use mgis.fenics to make use of MFront for the Mohr-Coulomb model. The governing equation is

div(σ) - I p + ρ b = 0

where σ is the effective stress, p is the pore pressure, ρ is the density, b is the volumetric force and I is the identity. In a full FEniCS, my working formulation is:

#u_trial: = trial function for displacement u
#u_test: test function for displacement u
#p: pore pressure
A_form = inner(sigma(u_trial),epsilon(u_test))*dx
A_residual = -inner(Identity(2)*p,epsilon(u_test))*dx-density*dot(body_force,u_test)*dx

In mgis.fenics, this part is handled by MFrontNonlinearProblem class. Without pore pressure p, one has

problem = MFrontNonlinearProblem(u, material, quadrature_degree=2, bcs=boundary_conditions)
problem.set_loading(density*dot(body_force,u_test)*problem.dx)

Hence, my understanding is that I only need to add a variant of “-inner(Identity(2)*p,epsilon(u_test))*dx” to the argument of set_loading. Going over the tutorials for mgis.fenics, I tried this formulation

problem = MFrontNonlinearProblem(u, material, quadrature_degree=2, bcs=boundary_conditions)
problem.initialize()
strain = problem.get_gradient("Strain")
problem.set_loading(-inner(Identity(2)*p,strain)*problem.dx-density*dot(body_force,u_test)*problem.dx)

but this returns an error: ufl.log.UFLException: Shapes do not match

What is the right approach? I am only starting to explore MFront and mgis.fenics so I appreciate any comments and help.

2

Hi, I think the reason for this is that strain is returned as a vector of components not a 2D tensor. Besides, in you loading definition you don’t want to use strain but the corresponding variation with respect to a test function u_test. So you can just replace this with sym(grad(u_test)).

1 Like
3

This works. You are right, I thought I can also get the variation with respect to the test function from the problem.get_gradient() but your suggestion is indeed the correct way. Thank you!