diff --git a/MaterialLib/SolidModels/MFront/DruckerPrager.mfront b/MaterialLib/SolidModels/MFront/DruckerPrager.mfront index bbb53778598a0f5deb768fd93b76d3f6d1e51136..fb10700c197b92d10c6f71d15ac042f5679c5f40 100644 --- a/MaterialLib/SolidModels/MFront/DruckerPrager.mfront +++ b/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; - } -} diff --git a/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj b/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj index 37a43241e6dd7abe8426e7e5b183e3fcfecbf0e7..8693f035896687a466af5d5437ab0cc8f2818609 100644 --- a/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj +++ b/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.