diff --git a/MaterialLib/SolidModels/MFront/CMakeLists.txt b/MaterialLib/SolidModels/MFront/CMakeLists.txt index 6facf3fc1416de3699a491ec8643f4b289a04838..a765275d31db16fa27bcd54fc1d8bc9677cd99ce 100644 --- a/MaterialLib/SolidModels/MFront/CMakeLists.txt +++ b/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) diff --git a/MaterialLib/SolidModels/MFront/Lubby2.mfront b/MaterialLib/SolidModels/MFront/Lubby2.mfront new file mode 100644 index 0000000000000000000000000000000000000000..698afb3b357732fe49f97400e93a6422b44b8d92 --- /dev/null +++ b/MaterialLib/SolidModels/MFront/Lubby2.mfront @@ -0,0 +1,86 @@ + +@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; +} diff --git a/MaterialLib/SolidModels/MFront/Lubby2.py b/MaterialLib/SolidModels/MFront/Lubby2.py new file mode 100644 index 0000000000000000000000000000000000000000..92905601e5f54ca233b32983912ba916a0696810 --- /dev/null +++ b/MaterialLib/SolidModels/MFront/Lubby2.py @@ -0,0 +1,69 @@ +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')