Merge pull request #2643 from endJunction/MFrontDruckerPrager

Update of DruckerPrager implementation.

MaterialLib/SolidModels/MFront/DruckerPrager.mfront

@@ -7,86 +7,83 @@
  *
  */
 
-@Parser Implicit;
+@DSL Implicit;
 @Behaviour DruckerPrager;
-@Algorithm NewtonRaphson_NumericalJacobian;
-@Theta 1.;
+@Author Thomas Nagel;
+@Description{
+Non-associated DP (perfect plasticity)
+F = alpha * I_1 + sqrt(J_2) - K
+};
 
-@MaterialProperty real young;
-young.setGlossaryName("YoungModulus");
-@MaterialProperty real nu;
-nu.setGlossaryName("PoissonRatio");
-@MaterialProperty real kappa;
-@MaterialProperty real beta;
+@Algorithm NewtonRaphson;
+@MaximumNumberOfIterations 100;
 
-@StateVariable real pla_mult;
-@StateVariable real a;
+@Brick StandardElasticity;
 
-@LocalVariable stress lambda;
-@LocalVariable stress mu;
-@LocalVariable stress Fel;
-@LocalVariable StressStensor s0;
+@Epsilon 1.e-14;
+@Theta 1.0;
+@Parameter local_zero_tolerance = 1.e-14;
 
-@InitLocalVariables
-{
-    lambda = computeLambda(young, nu);
-    mu = computeMu(young, nu);
-    StressStensor sigel(lambda * trace(eel + deto) * Stensor::Id() +
-                        2 * mu * (eel + deto));
-    const auto s_dev = deviator(sigel);
-    const stress J2 = (s_dev | s_dev) / 2.;
-    const stress I1 = trace(sigel);
+@ModellingHypotheses{".+"};
+@RequireStiffnessTensor<UnAltered>;
 
-    Fel = sqrt(J2) + beta * I1 - kappa;
-}
+// Intercept of yield function
+@MaterialProperty real K;
+K.setEntryName("Cohesion");
+// coefficient of I1 in yield function
+@MaterialProperty real alpha_y;
+alpha_y.setEntryName("FrictionParameter");
+// coefficient of I1 in plastic potential
+@MaterialProperty real alpha_g;
+alpha_g.setEntryName("DilatancyParameter");
+
+@StateVariable strain lam;
+lam.setGlossaryName("EquivalentPlasticStrain");
 
-@ComputeStress
+@LocalVariable bool F;
+
+@InitLocalVariables
 {
-    sig = (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel);
+    const auto sig_el = computeElasticPrediction();
+    const auto I1_el = trace(sig_el);
+    const auto s_el = deviator(sig_el);
+    const auto J2_el = max((s_el | s_el) / 2., local_zero_tolerance);
+    const auto sqrt_J2_el = sqrt(J2_el);
+    F = alpha_y * I1_el + sqrt_J2_el - K > 0.;
 }
 
 @Integrator
 {
-    Stensor s_dev = deviator(sig);
+    constexpr const auto id = Stensor::Id();
+    constexpr const auto id4 = Stensor4::Id();
+    const auto Pdev = id4 - (id ^ id) / 3;
 
-    if (Fel > 0)
+    if (F)
     {
-        Stensor nq = sqrt(3.) / 2. / sigmaeq(sig) * s_dev;
-        Stensor np = 1.0 * Stensor::Id();
+        const auto I1 = trace(sig);
+        const auto s = deviator(sig);
+        const auto J2 = max((s | s) / 2., local_zero_tolerance);
+        const auto sqrt_J2 = sqrt(J2);
+        const auto i_sqrt_J2 = 1. / sqrt_J2;
+        // yield function
+        const auto Fy = alpha_y * I1 + sqrt_J2 - K;
 
-        const auto yield =
-            sqrt((s_dev | s_dev) / 2.) + beta * trace(sig) - kappa;
+        // flow direction
+        const auto n = eval(alpha_g * id + s * i_sqrt_J2 / 2.);
 
-        const auto dev_flow = 1.0;
-        const auto iso_flow = beta;
+        // yield function gradient
+        const auto nF = eval(alpha_y * id + s * i_sqrt_J2 / 2.);
+        const auto dn_dsig =
+            eval(i_sqrt_J2 / 2. * (id4 - (s ^ s) / (2. * J2)) * Pdev);
 
-        feel = deel - deto + dpla_mult * (nq * dev_flow + np * iso_flow);
-        fpla_mult = yield / young;
-        fa = da - sqrt(2. / 3. * (dpla_mult * (nq * dev_flow + np * iso_flow)) |
-                       (dpla_mult * (nq * dev_flow + np * iso_flow)));
-    }
-    else
-    {
-        feel = deel - deto;
+        // residuals
+        feel += dlam * n;
+        flam = Fy / D(0, 0);
+        // Jacobian
+        dfeel_ddeel += theta * dlam * (dn_dsig * D);
+        dfeel_ddlam = n;
+        dflam_ddlam = strain(0.);
+        dflam_ddeel = theta * (nF | D) / D(0, 0);
     }
 }
 
-@TangentOperator
-{
-    if ((smt == ELASTIC) || (smt == SECANTOPERATOR))
-    {
-        computeElasticStiffness<N, Type>::exe(Dt, lambda, mu);
-    }
-    else if (smt == CONSISTENTTANGENTOPERATOR)
-    {
-        StiffnessTensor De;
-        Stensor4 Je;
-        computeElasticStiffness<N, Type>::exe(De, lambda, mu);
-        getPartialJacobianInvert(Je);
-        Dt = De * Je;
-    }
-    else
-    {
-        return false;
-    }
-}

Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj

@@ -13,8 +13,9 @@
                 <material_properties>
                     <material_property name="YoungModulus" parameter="E"/>
                     <material_property name="PoissonRatio" parameter="nu"/>
-                    <material_property name="beta" parameter="beta"/>
-                    <material_property name="kappa" parameter="kappa"/>
+                    <material_property name="Cohesion" parameter="kappa"/>
+                    <material_property name="FrictionParameter" parameter="beta"/>
+                    <material_property name="DilatancyParameter" parameter="beta"/>
                 </material_properties>
             </constitutive_relation>
             <!-- Configuration of the Ehlers model used to compute the reference solution.