Active-set algorithm for finite strain elastoplasticity and frictional contact

Francesco · December 10, 2024, 6:57am

Dear all,

I would like to raise some points that are not so clear to me.

I have been perusing the documentation of GetFEM with my colleagues and we found in the Python documentation, in particular, the following function:

add_finite_strain_elastoplasticity_brick(mim, lawname, unknowns_type, varnames=None, *args)

Could you please confirm that this Simo-Miehe finite-strain elastoplasticity model has actually been tested and implemented? Which algorithm is employed for the numerical discretization?

Also, concerning frictional contact: is the active-set algorithm already implemented in GetFEM, or would it represent a novel contribution?

If the active-set method is already implemented, it could be used also for finite-strain elastoplasticity. Also, in the general documentation, I read

In GetFEM we propose both the return mapping strategy and also an alternative strategy developed below which is mainly inspired from [PO-NI2016], [SE-PO-WO2015] and [HA-WO2009] (Active Set) and allow more simple tangent moduli.

Does it mean that the active-set method is available both for frictional contact and finite-strain elastoplasticity?

Would it be usable also via the Python interface?

Thank you very much for your feedback.

Konstantinos.Poulios · December 10, 2024, 11:09am

Active set method is equivalent to properly applied Newton-Raphson on semi-smooth problems.

https://epubs.siam.org/doi/10.1137/S1052623401383558

“Active set” is a good marketing term in the literature but it is not really an as distinct method as it gives the impression of. I prefer to refer to these methods as semi-smooth Newton or generalized Newton.

So regarding your question, the contact and plasticity methods implemented in GetFEM are equivalent to those that you will find under the “active set” term in the literature. There is nothing new to add there.

The finite strain elastoplasticity brick in GetFEM as far as I remember is based on

Simo, J., 1992. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory.

more detailed described in Chapter 14 of “Computational Methods for Plasticity” 2008.

If you only need linear hardening, rather than using the finite strain elastoplasticity brick, I would recommend you to use the implementation found in the examples:

https://git.savannah.nongnu.org/cgit/getfem.git/tree/contrib/continuum_mechanics

where you can see the actual equations been solved rather than using a black box.

Konstantinos.Poulios · December 10, 2024, 11:44am

if you want to use nonlinear hardneing laws, you can use the following implementation of the governing equations, to avoid black-box solutions. You can choose whether you use the pressure field as an extra variable or not.

md = gf.Model("real")

md.add_fem_variable("u", mfu)       # displacements
if mixed_up:
   md.add_fem_variable("p", mfp)    # pressure
md.add_internal_im_variable("dksi", mimd1)   # plastic multiplier increase in the current step

md.add_im_data("ksi0", mimd1)          # plastic multiplier value at previous time step
md.add_im_data("invCp0vec", mimd6)     # Components 11, 22, 33, 12, 13, 23 of the plastic part
md.set_variable("invCp0vec",           # of the inverse right Cauchy Green tensor at the
                np.tile([1,1,1,0,0,0], # previous step.
                mimd6.nbpts()))        #

md.add_initialized_data("gamma", [0.])          # Current load level (in percent)
md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
md.add_initialized_data("mu", E/(2.*(1.+nu)))   # Shear modulus

md.add_macro("Mat3x3ToVec6", "[[[1,0,0,0,0,0]  ,[0,0,0,0.5,0,0],[0,0,0,0,0.5,0]],"\
                             " [[0,0,0,0.5,0,0],[0,1,0,0,0,0]  ,[0,0,0,0,0,0.5]],"+\
                             " [[0,0,0,0,0.5,0],[0,0,0,0,0,0.5],[0,0,1,0,0,0]  ]]")
md.add_macro("Vec6ToMat3x3", "[[[1,0,0],[0,0,0],[0,0,0]],"\
                             " [[0,0,0],[0,1,0],[0,0,0]],"\
                             " [[0,0,0],[0,0,0],[0,0,1]],"\
                             " [[0,1,0],[1,0,0],[0,0,0]],"\
                             " [[0,0,1],[0,0,0],[1,0,0]],"\
                             " [[0,0,0],[0,0,1],[0,1,0]]]")

md.add_macro("F", "Id(3)+Grad(u)")
md.add_macro("J", "Det(F)")
md.add_macro("invCp0", "Vec6ToMat3x3.invCp0vec")
md.add_macro("devlogbetr", "Deviator(Logm(F*invCp0*F'))")
md.add_macro("devlogbe", "(1-2*dksi)*devlogbetr")
md.add_macro("tauD", "mu*pow(J,-2/3)*devlogbe")
if mixed_up:
   md.add_macro("tauH", "-p*J")
else:
   md.add_macro("tauH", "K*log(J)")
md.add_macro("tau", "tauH*Id(3)+tauD")
md.add_nonlinear_term(mim, "((tauH*Id(3)+tauD)*Inv(F)'):Grad(Test_u)")

if mixed_up:
   md.add_nonlinear_term(mimpress, "(p/K+log(J)/J)*Test_p")

md.add_macro("ksi", "ksi0+dksi*Norm(devlogbetr)")
if plast_law == "Ramberg Osgood":
   md.add_macro("Y", "{coef}*pow(ksi+1e-12,{exponent})"\
                     .format(coef=pow(E/pl_alpha * pow(2./3., (pl_n+1)/2.) * pow(pl_sigma_0, pl_n-1),
                                      1./pl_n),
                             exponent=1./pl_n))
elif plast_law == "Linear hardening":
   md.add_macro("Y", "{A}+{B}*ksi".format(A=np.sqrt(2./3.)*pl_sigma_0, B=2./3.*pl_H))
elif plast_law == "Linear and exponential hardening":
   md.add_macro("Y", "{A}+{B}*(1-exp(-{C}*ksi))+{D}*ksi"\
                     .format(A=np.sqrt(2./3.)*pl_sigma_0, B=np.sqrt(2./3.)*(pl_sigma_inf-pl_sigma_0),
                             C=np.sqrt(2./3.)*pl_delta, D=2./3.*pl_H))
else:
   raise NameError(plast_law)

md.add_macro("sigmaD", "mu*devlogbe")
md.add_macro("sigmaDtr", "mu*devlogbetr")
md.add_nonlinear_term(mim, "(Norm(sigmaD)-min(Y,Norm(sigmaDtr)+1e-12*dksi))*Test_dksi")

Konstantinos.Poulios · December 17, 2024, 4:33pm

@Francesco

There was a mistake that I have just fixed in the code from my previous message.

I also have an update regarding which version is implemented in GetFEM, although it is called Simo-Miehe, the actual implementation is with logarithmic strains according to Eterovic-Bathe.

The confusion came from the fact that Simo had also used logarithmic strains but in the other paper from 1992, that I mentioned above, not the one with Miehe. Simo’s model with logarithmic strains is equivalent to the Eterovic-Bathe, which is implemented in GetFEM.

if your want to compare the two models

logarithmic strain based
Simo-Miehe

you can have a look at this code which implements both

md = gf.Model("real")

md.add_fem_variable("u", mfu)       # displacements
if mixed_up:
   md.add_fem_variable("p", mfp)    # pressure
md.add_internal_im_variable("ksi", mimd1)   # plastic multiplier increase in the current step

md.add_im_data("gamma0", mimd1)        # accumulated plastic strain at previous time step
md.add_im_data("invCp0vec", mimd6)     # Components 11, 22, 33, 12, 13, 23 of the plastic part
md.set_variable("invCp0vec",           # of the inverse right Cauchy Green tensor at the
                np.tile([1,1,1,0,0,0], # previous step.
                mimd6.nbpts()))        #

md.add_initialized_data("loadfactor", [0.])     # Current load level (in percent)
md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
md.add_initialized_data("G", E/(2.*(1.+nu)))    # Shear modulus

md.add_macro("Mat3x3ToVec6", "[[[1,0,0,0,0,0]  ,[0,0,0,0.5,0,0],[0,0,0,0,0.5,0]],"\
                             " [[0,0,0,0.5,0,0],[0,1,0,0,0,0]  ,[0,0,0,0,0,0.5]],"+\
                             " [[0,0,0,0,0.5,0],[0,0,0,0,0,0.5],[0,0,1,0,0,0]  ]]")
md.add_macro("Vec6ToMat3x3", "[[[1,0,0],[0,0,0],[0,0,0]],"\
                             " [[0,0,0],[0,1,0],[0,0,0]],"\
                             " [[0,0,0],[0,0,0],[0,0,1]],"\
                             " [[0,1,0],[1,0,0],[0,0,0]],"\
                             " [[0,0,1],[0,0,0],[1,0,0]],"\
                             " [[0,0,0],[0,0,1],[0,1,0]]]")

md.add_macro("F", "Id(3)+Grad(u)")
md.add_macro("J", "Det(F)")
md.add_macro("invCp0", "Vec6ToMat3x3.invCp0vec")
if plast_version == "EterovicBathe":
   md.add_macro("devlogbetr", "Deviator(Logm(F*invCp0*F'))")
   md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
   md.add_macro("tauD", "G*devlogbe")
elif plast_version == "SimoMiehe":
   md.add_macro("betr", "F*invCp0*F'")
   md.add_macro("devbetr", "Deviator(betr)")
   md.add_macro("devbe", "(1-2/3*Trace(betr)*ksi)*devbetr")
   md.add_macro("tauD", "G*pow(J,-2/3)*devbe")
else:
   raise NameError(plast_version)

if mixed_up:
   md.add_macro("tauH", "-p*J")
else:
   md.add_macro("tauH", "K*log(J)")
md.add_macro("tau", "tauH*Id(3)+tauD")
if integration_degree_vol == integration_degree:
   md.add_nonlinear_term(mim, "((tauH*Id(3)+tauD)*Inv(F)'):Grad(Test_u)")
else:
   md.add_nonlinear_term(mim, "(tauD*Inv(F)'):Grad(Test_u)")
   md.add_nonlinear_term(mim_vol, "tauH*Inv(F)':Grad(Test_u)")

if mixed_up:
   md.add_nonlinear_term(mimpress, "(p/K+log(J)/J)*Test_p")

if plast_hardening == "Ramberg Osgood":
   md.add_macro("Y(gamma)",
                "{coef}*pow(gamma+1e-12,{exponent})"\
                .format(coef=pow(E/pl_alpha * pow(2./3., (pl_n+1)/2.) * pow(pl_sigma_0, pl_n-1),
                                 1./pl_n),
                        exponent=1./pl_n))
elif plast_hardening == "Linear hardening":
   md.add_macro("Y(gamma)", "{A}+{B}*gamma".format(A=np.sqrt(2./3.)*pl_sigma_0, B=2./3.*pl_H))
elif plast_hardening == "Linear and exponential hardening":
   md.add_macro("Y(gamma)",
                "{A}+{B}*(1-exp(-{C}*gamma))+{D}*gamma"\
                .format(A=np.sqrt(2./3.)*pl_sigma_0, B=np.sqrt(2./3.)*(pl_sigma_inf-pl_sigma_0),
                        C=np.sqrt(2./3.)*pl_delta, D=2./3.*pl_H))
else:
   raise NameError(plast_hardening)
if plast_version == "EterovicBathe":
   md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
   md.add_macro("tauDtr", "G*devlogbetr")
   #md.add_macro("sigmaD", "G/J*devlogbe")
   #md.add_macro("sigmaDtr", "G/J*devlogbetr")
elif plast_version == "SimoMiehe":
   md.add_macro("gamma", "gamma0+ksi*Norm(devbetr)")
   md.add_macro("tauDtr", "G*pow(J,-2/3)*devbetr")
   #md.add_macro("sigmaD", "G*pow(J,-5/3)*devbe")
   #md.add_macro("sigmaDtr", "G*pow(J,-5/3)*devbetr")
md.add_nonlinear_term(mim, "(Norm(tauD)-min(Y(gamma),Norm(tauDtr)+1e-12*ksi))*Test_ksi")

Francesco · December 18, 2024, 7:42am

Thank you very much, I will check this out!

Topic		Replies	Views
Implementation of the Gonzalez time scheme for dynamic hyperelasticity Questions and issues	16	93	February 13, 2025
Coupling matrix for (frictional) contact problems General	5	29	June 5, 2025
Self-contact only — Issue with augmented Lagrangian method General	2	28	May 12, 2025
Simple 2D isotropic linear elasticity problem on the unit square with a mesh of 2 triangles General	8	29	June 23, 2025
Any examples of XFEM in action? Questions and issues	2	31	November 26, 2024

Active-set algorithm for finite strain elastoplasticity and frictional contact

Related topics