From 3d21014a470359df5a49d25368d168beceb98e9c Mon Sep 17 00:00:00 2001
From: cbsilver <christian.silbermann@ifgt.tu-freiberg.de>
Date: Wed, 18 Jan 2023 11:48:01 +0100
Subject: [PATCH] Hypoelastic MCC models: absolute formulation with analytical
solution for pc
---
...ModCamClay_semiExplParaInitNLnu_abs.mfront | 63 +++++++++----------
1 file changed, 30 insertions(+), 33 deletions(-)
diff --git a/MaterialLib/SolidModels/MFront/ModCamClay_semiExplParaInitNLnu_abs.mfront b/MaterialLib/SolidModels/MFront/ModCamClay_semiExplParaInitNLnu_abs.mfront
index 5ff42dd92c8..9d37ea6635a 100644
--- a/MaterialLib/SolidModels/MFront/ModCamClay_semiExplParaInitNLnu_abs.mfront
+++ b/MaterialLib/SolidModels/MFront/ModCamClay_semiExplParaInitNLnu_abs.mfront
@@ -18,7 +18,6 @@
eel, deel (elastic strain (increment))
sig (stress)
dlp (plastic increment)
- dpc (pre-consolidation pressure increment)
Output: feel (strain residual depending on deel, dlp, dpc)
flp (yield function residual depending on deel, dpc)
@@ -67,10 +66,6 @@ pc0.setEntryName("InitialPreConsolidationPressure");
@StateVariable real lp;
lp.setGlossaryName("EquivalentPlasticStrain");
-// @StateVariable stress pc;
-// pc.setEntryName("ReducedPreConsolidationPressure");
-@IntegrationVariable strain rpc;
-
// An auxiliary state variable is a persistent variable but not an integration variable.
@AuxiliaryStateVariable stress pc;
pc.setEntryName("PreConsolidationPressure");
@@ -93,7 +88,6 @@ p0.setEntryName("InitialPressure");
@LocalVariable bool withinElasticRange;
@LocalVariable real M2;
@LocalVariable real young;
-@LocalVariable real rpcMin;
@Includes{
#ifndef MFRONT_PRESSUREDEPENDANTBULKMODULUS_IMPLEMENTATION
@@ -114,11 +108,10 @@ p0.setEntryName("InitialPressure");
const auto s = deviator(eel);
const auto r = 3 * (1 - 2 * nu) / (2 * (1 + nu));
- double K, G, p;
// explicit computation of the hydrostatic pressure
- p = p0 * exp(-v0_ka * e);
- K = v0_ka * p;
- G = r * K; // dG_de = r * dK_de = r * dK_dp * dp_de = - G * dK_dp;
+ const auto p = p0 * exp(-v0_ka * e);
+ const auto K = v0_ka * p;
+ const auto G = r * K; // dG_de = r * dK_de = r * dK_dp * dp_de = - G * dK_dp;
const auto dp_de = -K;
const auto dK_dp = v0_ka;
// Hooke's law
@@ -152,10 +145,7 @@ p0.setEntryName("InitialPressure");
p0 = p;
v = v0;
}
-
young = 3.0 * K * (1.0 - 2*nu);
- rpc = pc / young;
- rpcMin = pcMin / young;
// computation of the elastic prediction (does not work for plane stress!)
const auto eps_el = StrainStensor{eel + deto};
@@ -190,56 +180,63 @@ p0.setEntryName("InitialPressure");
const auto s = deviator(sig);
const auto q = std::sqrt(1.5 * s | s);
const auto p = -trace(sig) / 3;
- // update the internal (state) variables (rpc holds the old value!)
- const auto rpcNew = rpc + theta * drpc;
- const auto pcNew = rpcNew * young;
+ // update the internal (state) variables (pc holds the old value!)
+ const auto deelV = trace(deel);
+ const auto detoV = trace(deto);
+ auto deplV = detoV - deelV;
+ auto eplV = epl_V + theta * deplV;
+ const auto pcNew = (pc0 - pcMin) * exp(-the * eplV) + pcMin; // TODO only for eplV0=0!
// calculate the direction of plastic flow
- const auto f = (q * q + M2 * p * (p - pcNew));
+ const auto f = q * q + M2 * p * (p - pcNew);
const auto df_dp = M2 * (2 * p - pcNew);
+ const auto df_dpc = -M2 * p;
const auto df_dsig = eval(3 * s - df_dp * id2 / 3);
auto norm = std::sqrt(6 * q * q + df_dp * df_dp / 3); // = std::sqrt(df_dsig|df_dsig);
norm = std::max(norm, epsr * young);
const auto n = df_dsig / norm;
const auto ntr = -df_dp / norm;
// plastic strain and volumetric part
- auto depl = eval(dlp * n);
- const auto deplV = trace(depl);
+ const auto depl = eval(dlp * n);
+ deplV = trace(depl);
+ eplV = epl_V + theta * deplV;
const auto fchar = pc0 * young; // OR: young*young
// residual
feel = deel + depl - deto;
flp = f / fchar;
- frpc = drpc + deplV * the * (rpcNew - rpcMin);
// auxiliary derivatives
const auto dnorm_dsig = (9 * s - 2 * M2 / 9 * df_dp * id2) / norm;
const auto dn_ddeel = (3 * Pr4 + 2 * M2 / 9 * (id2 ^ id2) - (n ^ dnorm_dsig)) / norm * dsig_deel * theta;
- const auto dn_ddrpc = (id2 + df_dp * n / norm) * M2 / (3 * norm) * theta * young;
- const auto dfrpc_ddeplV = the * (rpcNew - rpcMin);
+ const auto dn_dpc = (id2 + df_dp * n / norm) * M2 / (3 * norm);
+ const auto dpc_deplV = -the * (pcNew - pcMin);
+ const auto ddeplV_ddlp = ntr;
+ const auto ddeplV_dn = eval(dlp * id2);
+ //const auto dpc_dn = dpc_deplV * theta * ddeplV_dn;
- // jacobian (all other parts are zero)
- dfeel_ddeel += dlp * dn_ddeel;
- dfeel_ddlp = n;
- dfeel_ddrpc = dlp * dn_ddrpc;
+ const auto dpc_ddlp = dpc_deplV * theta * ddeplV_ddlp;
+ const auto dpc_ddeel = dpc_deplV * theta * (ddeplV_dn | dn_ddeel);
+ //const auto dpc_ddeel = dpc_deplV * theta * dlp * (id2 | dn_ddeel);
- dflp_ddeel = (df_dsig | dsig_deel) * theta / fchar; // in case of problems with zero use:
- dflp_ddlp = strain(0); // (q<epsr) ? strain(1) : strain(0);
- dflp_ddrpc = -M2 * p * theta / fchar * young;
+ // jacobian (all other parts are zero)
+ const auto dfeel_dpc = dlp * dn_dpc;
+ dfeel_ddeel += dlp * dn_ddeel + (dfeel_dpc ^ dpc_ddeel);
+ dfeel_ddlp = n + dfeel_dpc * dpc_ddlp;
- dfrpc_ddlp = dfrpc_ddeplV * ntr;
- dfrpc_ddeel = dfrpc_ddeplV * dlp * (id2 | dn_ddeel);
- dfrpc_ddrpc = 1 + (deplV * the + dfrpc_ddeplV * dlp * trace(dn_ddrpc)) * theta;
+ const auto dflp_dpc = df_dpc / fchar;
+ dflp_ddeel = (df_dsig | dsig_deel) * theta / fchar + dflp_dpc * dpc_ddeel;
+ dflp_ddlp = strain(0) + dflp_dpc * dpc_ddlp;
}
// explicit treatment as long as change of v (or e) during time increment is small
@UpdateAuxiliaryStateVariables
{
Time += dt;
- pc += drpc * young;
const auto deelV = trace(deel);
const auto detoV = trace(deto);
epl_V += detoV - deelV;
+ pc = (pc0 - pcMin) * exp(-v0 / (la - ka) * epl_V) + pcMin;
v += v0 * detoV;
}
--
GitLab