Added Lubby2 model based on MFront implementation.

MaterialLib/SolidModels/MFront/CMakeLists.txt

@@ -18,7 +18,8 @@ mfront_behaviours_check_library(OgsMFrontBehaviour
                                 MohrCoulombAbboSloan
                                 MohrCoulombAbboSloanOrtho
                                 StandardElasticityBrick
-                                StandardElasticityBrickOrtho)
+                                StandardElasticityBrickOrtho
+                                Lubby2)
 
 target_include_directories(MaterialLib_SolidModels_MFront
                            PUBLIC ThirdParty/MGIS/include)

MaterialLib/SolidModels/MFront/Lubby2.mfront

+
+@DSL Implicit;
+@Behaviour Lubby2;
+@Author Thomas Nagel;
+@Description{Lubby2 model for stationary
+and transient creep based on Burgers
+rheological model};
+
+@Algorithm NewtonRaphson;
+@MaximumNumberOfIterations 100;
+
+@Brick StandardElasticity;
+
+@Epsilon 1.e-14;
+@Theta 1.0;
+@Parameter local_zero_tolerance = 1.e-14;
+
+@ModellingHypotheses{".+"};
+@RequireStiffnessTensor<UnAltered>;
+
+// Intercept of yield function
+@MaterialProperty real GK0;
+GK0.setEntryName("KelvinShearModulus");
+@MaterialProperty real etaK0;
+etaK0.setEntryName("KelvinViscosity");
+@MaterialProperty real etaM0;
+etaM0.setEntryName("MaxwellViscosity");
+@MaterialProperty real mK;
+mK.setEntryName("KelvinElasticParameter");
+@MaterialProperty real mvK;
+mvK.setEntryName("KelvinViscoParameter");
+@MaterialProperty real mvM;
+mvM.setEntryName("MaxwellViscoParameter");
+
+@StateVariable Stensor epsK;
+epsK.setEntryName("KelvinStrain");
+
+@StateVariable Stensor epsM;
+epsM.setEntryName("MaxwellStrain");
+
+@Integrator
+{
+    // Implicit system
+    constexpr const auto id4 = Stensor4::Id();
+    const auto Pdev = Stensor4::K();
+
+    const auto sigma_eff =
+        std::max(sigmaeq(sig), local_zero_tolerance * D(0, 0));
+
+    const auto s = deviator(sig);
+
+    const auto etaK = etaK0 * std::exp(mvK * sigma_eff);
+    const auto etaM = etaM0 * std::exp(mvM * sigma_eff);
+    const auto GK = GK0 * std::exp(mK * sigma_eff);
+
+    // residuals
+    feel += depsM + depsK;
+
+    // calculate Kelvin strain residual
+    fepsK = depsK - dt / (2.0 * etaK) * (s - 2.0 * GK * (epsK + depsK));
+
+    // calculate Maxwell strain residual
+    fepsM = depsM - dt / (2.0 * etaM) * s;
+
+    // Jacobian
+    const auto dsigma_eff_ddeel = 3. / (2. * sigma_eff) * s | (Pdev * D);
+    const auto detaK_ddeel = etaK * mvK * dsigma_eff_ddeel;
+    const auto detaM_ddeel = etaM * mvM * dsigma_eff_ddeel;
+    const auto dGK_ddeel = GK * mK * dsigma_eff_ddeel;
+
+    dfeel_ddeel = id4;
+    dfeel_ddepsK = id4;
+    dfeel_ddepsM = id4;
+
+    dfepsK_ddeel = (dt / (2.0 * etaK * etaK) * (s - 2.0 * GK * (epsK + depsK)) ^
+                    detaK_ddeel) -
+                   (dt / (2.0 * etaK) * Pdev * D) +
+                   (dt / etaK * (epsK + depsK) ^ dGK_ddeel);
+    dfepsK_ddepsK = (1. + dt * GK / etaK) * id4;
+    // dfepsK_ddepsM is all zero --> nothing to do
+
+    dfepsM_ddeel = (dt / (2.0 * etaM * etaM) * s ^ detaM_ddeel) -
+                   (dt / (2.0 * etaM) * Pdev * D);
+    // dfepsM_ddepsK is all zero --> nothing to do
+    dfepsM_ddepsM = id4;
+}

MaterialLib/SolidModels/MFront/Lubby2.py

+import mtest
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+GM0  = 9.54e3
+KM0  = 2.78e4
+GK0  = 6.27e4
+mK   = -0.254
+mvK  = -0.327
+mvM  = -0.267
+etaK0= 1.66e5
+etaM0= 4.03e7
+tmax = 15.
+
+
+m    = mtest.MTest()
+mtest.setVerboseMode(mtest.VerboseLevel.VERBOSE_QUIET)
+m.setMaximumNumberOfSubSteps(10)
+m.setBehaviour('generic', 'src/libBehaviour.so', 'Lubby2')
+m.setExternalStateVariable("Temperature", 293.15)
+m.setImposedStress('SXX', 0.0)
+m.setImposedStress('SYY', 0.0)
+m.setImposedStress('SZZ', 0.0)
+m.setImposedStress('SXZ', 0.0)
+m.setImposedStress('SYZ', 0.0)
+m.setImposedStress('SXY', {0: 0, 1.e-5 : 5 * np.sqrt(2.), 15 : 5 * np.sqrt(2.)}) #Kelvin Mapping of shear components!!
+
+m.setMaterialProperty('YoungModulus', (9 * GM0 * KM0) / (3 * KM0 + GM0)) #this was the wrong equation
+m.setMaterialProperty('PoissonRatio', (3 * KM0 - 2 * GM0) / (2 * GM0 + 6 * KM0)) #this was the wrong equation
+m.setMaterialProperty('KelvinShearModulus', GK0)
+m.setMaterialProperty('KelvinViscosity', etaK0)
+m.setMaterialProperty('MaxwellViscosity', etaM0)
+m.setMaterialProperty('KelvinElasticParameter', mK)
+m.setMaterialProperty('MaxwellViscoParameter', mvM) #these were switched
+m.setMaterialProperty('KelvinViscoParameter', mvK) #these were switched
+
+sig_xy = 5.
+sig_eff = np.sqrt(3.)*sig_xy
+
+GK=GK0*np.exp(mK*sig_eff)
+etaK=etaK0*np.exp(mvK*sig_eff)
+etaM=etaM0*np.exp(mvM*sig_eff)
+
+
+eps_xy = lambda t: ((1./GM0 + t/etaM)*sig_xy + 1./GK*(1.-np.exp(-GK/etaK*t))*sig_xy)/2.
+
+s = mtest.MTestCurrentState()
+wk = mtest.MTestWorkSpace()
+m.completeInitialisation()
+m.initializeCurrentState(s)
+m.initializeWorkSpace(wk)
+
+ltime = np.append(np.linspace(0,1e-5,10),np.linspace(2.e-5,tmax,200))
+
+numerical = np.array([0.])
+
+#run sim
+for i in range(len(ltime)-1):
+    m.execute(s, wk, ltime[i], ltime[i + 1])
+    numerical = np.append(numerical,s.e1[3]/np.sqrt(2.)) #Kelvin mapping backwards!
+
+fig,ax = plt.subplots()
+ax.plot(ltime,numerical*100.,label='numerical')
+ax.plot(ltime,eps_xy(ltime)*100.,label='analytical')
+ax.set_xlabel('$t$ / d')
+ax.set_ylabel('$\\epsilon_{xy}$ / %')
+fig.legend()
+fig.savefig('lubby_mfront.pdf')