diff --git a/MaterialLib/SolidModels/MFront/BDT.mfront b/MaterialLib/SolidModels/MFront/BDT.mfront new file mode 100644 index 0000000000000000000000000000000000000000..85ae6ce7a3c70d686797d5e4fd7284ec6236bba2 --- /dev/null +++ b/MaterialLib/SolidModels/MFront/BDT.mfront @@ -0,0 +1,164 @@ +@Parser Implicit; +@Behaviour BDT; +@Algorithm NewtonRaphson_NumericalJacobian; +@Theta 1.; + +@MaterialProperty real young; +young.setGlossaryName("YoungModulus"); +@MaterialProperty real nu; +nu.setGlossaryName("PoissonRatio"); +@MaterialProperty thermalexpansion alpha; +alpha.setGlossaryName("ThermalExpansion"); +@MaterialProperty real f_c; +@MaterialProperty real m_0; +@MaterialProperty real alpha_p; +@MaterialProperty real n_exp_T; +@MaterialProperty real q_h0; +@MaterialProperty real chi_h; +@MaterialProperty real alpha_d; +@MaterialProperty real h_d; +@MaterialProperty real Qact; +@MaterialProperty real A_creep; +@MaterialProperty real n_creep; +@MaterialProperty real El_1; +@MaterialProperty real El_2; +@MaterialProperty real El_3; +@MaterialProperty real at_1; +@MaterialProperty real at_2; +@MaterialProperty real at_3; + +@StateVariable real pla_mult; +@StateVariable real a; +@StateVariable real dam; +@StateVariable real Rv; +@StateVariable real triax_p; +@StateVariable real vp; + +@LocalVariable stress lambda; +@LocalVariable stress mu; +@LocalVariable stress Fel; +@LocalVariable StressStensor s0; + +// @Parameter Tref; +// Tref.setDefaultValue(293.15); + +@InitLocalVariables +{ + auto const T_Celsius = T - 273.15; + + const auto y2 = El_1 * power<2>(T_Celsius) + El_2 * T_Celsius + El_3; + lambda = computeLambda(young * y2, nu); + mu = computeMu(young * y2, nu); + StressStensor sigel(lambda * trace(eel + deto) * Stensor::Id() + + 2 * mu * (eel + deto)); // StresssStensor in Tutorial + const auto s_dev = deviator(sigel); + // J2 = (s_dev|s_dev)/2.; //defines double contraction + const stress seq = sigmaeq(s_dev); + const stress I1 = trace(sigel); + const auto aux = (seq + I1) / (3. * f_c); + + const auto aux_ex = (1. - 1. / n_exp_T); + const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10.), n_exp_T); + + Rv = max(at_1 * T_Celsius + at_2 * I1 / 3. + at_3, 1.0e-4); + auto Rpel = -(1. - q_h0) * pow((a + vp) / Rv, 2) + + 2. * (1. - q_h0) * (a + vp) / Rv + q_h0; + if ((a + vp) > Rv) + Rpel = 1.0; + + const auto q_h = Rpel / (pow(aux_arg, aux_ex)); + + Fel = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) + + (m_0 * aux - 1) * power<2>(q_h); +} + +@ComputeStress +{ + sig = (1 - dam) * (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel); +} + +@Integrator +{ + auto const T_Celsius = T - 273.15; + + Stensor nvp = Stensor(0.); + Stensor s_dev = deviator(sig / (1 - dam)); + const stress seq = sigmaeq(s_dev); + + if (seq > 1.e-15) + { + nvp = 1.5 * s_dev / seq; + } + fvp -= + dt * A_creep * exp(-Qact / 8.3144598 / T_Celsius) * pow(seq, n_creep); + + if (Fel > 0) + { + const real pla_mult_ = pla_mult + theta * dpla_mult; + Stensor nq = s_dev; + Stensor np = 1.0 * Stensor::Id(); + const auto aux = (seq + trace(sig / (1 - dam))) / (3. * f_c); + + const auto aux_ex = (1. - 1. / n_exp_T); + const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10), n_exp_T); + real a_; + real dam_; + real vp_; + + a_ = a + theta * da; + dam_ = dam + theta * ddam; + vp_ = vp + theta * dvp; + Rv = max(at_1 * T_Celsius + at_2 * trace(sig / (1 - dam)) / 3. + at_3, + 1.0e-4); + auto Rp_ = -(1. - q_h0) * pow((a_ + vp_) / Rv, 2) + + 2. * (1. - q_h0) * (a_ + vp_) / Rv + q_h0; + if ((a_ + vp_) > Rv) + Rp_ = 1.0; + const auto q_h = Rp_ / (pow(aux_arg, aux_ex)); + + const auto yield = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) + + (m_0 * aux - 1) * power<2>(q_h); + + const auto big_aux = power<2>(aux) * (1 - q_h) + seq / f_c; + const auto dev_flow = + power<2>(q_h) * m_0 / (2. * f_c * seq) + + 2. * ((1 - q_h) * aux / (f_c * seq) + 3. / (2. * f_c * seq)) * + big_aux; + const auto iso_flow = power<2>(q_h) * m_0 / (3. * f_c) + + 4. * (1. - q_h) * aux / (3. * f_c) * big_aux; + + feel = deel - deto + dpla_mult * (nq * dev_flow + np * iso_flow) + + dvp * nvp; + 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))); + triax_p = iso_flow / dev_flow; + + dam = max(min(1 - exp(-(((a + vp) - Rv) / alpha_d)), 1.0), 0.0); + sig = (1 - dam) * sig; + } + else + { + feel = deel - deto + dvp * nvp; + } +} + +@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/ThermoMechanics/BDT/bdt.mtest b/MaterialLib/SolidModels/MFront/BDT.mtest similarity index 96% rename from Tests/Data/ThermoMechanics/BDT/bdt.mtest rename to MaterialLib/SolidModels/MFront/BDT.mtest index d2c5b2b1a296aaa40abcf1f041b2a27391ad0b2e..4b058e896512b0e6729dd67e0c945ea61c21eafe 100644 --- a/Tests/Data/ThermoMechanics/BDT/bdt.mtest +++ b/MaterialLib/SolidModels/MFront/BDT.mtest @@ -1,4 +1,4 @@ -@Behaviour<generic> 'src/libBehaviour.so' 'BDT'; +@Behaviour<generic> 'libOgsMFrontBehaviour.so' 'BDT'; @MaterialProperty<constant> 'YoungModulus' 4.4619E+04; @MaterialProperty<constant> 'PoissonRatio' 3.0000E-01; @MaterialProperty<constant> 'ThermalExpansion' 1.0000E-05; diff --git a/MaterialLib/SolidModels/MFront/CMakeLists.txt b/MaterialLib/SolidModels/MFront/CMakeLists.txt index 4e1a425beb9a311ae11b891c4309499a4d4a6cb5..54df69a8d44b3b6ac41bde131ebbb00dc5e7c38a 100644 --- a/MaterialLib/SolidModels/MFront/CMakeLists.txt +++ b/MaterialLib/SolidModels/MFront/CMakeLists.txt @@ -16,7 +16,10 @@ target_link_libraries(MaterialLib_SolidModels_MFront set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/ThirdParty/MGIS") include(cmake/modules/tfel) -mfront_behaviours_check_library(OgsMFrontBehaviour Elasticity +mfront_behaviours_check_library(OgsMFrontBehaviour + BDT + DruckerPrager + Elasticity StandardElasticityBrick) target_include_directories(MaterialLib_SolidModels_MFront diff --git a/MaterialLib/SolidModels/MFront/DruckerPrager.mfront b/MaterialLib/SolidModels/MFront/DruckerPrager.mfront new file mode 100644 index 0000000000000000000000000000000000000000..5f138f4ea68b1ed99fe66288d2b182b4c05c630b --- /dev/null +++ b/MaterialLib/SolidModels/MFront/DruckerPrager.mfront @@ -0,0 +1,83 @@ +@Parser Implicit; +@Behaviour DruckerPrager; +@Algorithm NewtonRaphson_NumericalJacobian; +@Theta 1.; + +@MaterialProperty real young; +young.setGlossaryName("YoungModulus"); +@MaterialProperty real nu; +nu.setGlossaryName("PoissonRatio"); +@MaterialProperty real kappa; +@MaterialProperty real beta; + +@StateVariable real pla_mult; +@StateVariable real a; + +@LocalVariable stress lambda; +@LocalVariable stress mu; +@LocalVariable stress Fel; +@LocalVariable StressStensor s0; + +@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); + + Fel = sqrt(J2) + beta * I1 - kappa; +} + +@ComputeStress +{ + sig = (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel); +} + +@Integrator +{ + Stensor s_dev = deviator(sig); + + if (Fel > 0) + { + Stensor nq = sqrt(3.) / 2. / sigmaeq(sig) * s_dev; + Stensor np = 1.0 * Stensor::Id(); + + const auto yield = + sqrt((s_dev | s_dev) / 2.) + beta * trace(sig) - kappa; + + const auto dev_flow = 1.0; + const auto iso_flow = beta; + + 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; + } +} + +@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/DruckerPrager.mfront b/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.mfront deleted file mode 100644 index b6a5571748a8b3ed43b1f3b402da321050af8aab..0000000000000000000000000000000000000000 --- a/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.mfront +++ /dev/null @@ -1,72 +0,0 @@ -@Parser Implicit; -@Behaviour DruckerPrager; -@Algorithm NewtonRaphson_NumericalJacobian; -@Theta 1.; - -@MaterialProperty real young; -young.setGlossaryName("YoungModulus"); -@MaterialProperty real nu; -nu.setGlossaryName("PoissonRatio"); -@MaterialProperty real kappa; -@MaterialProperty real beta; - -@StateVariable real pla_mult; -@StateVariable real a; - -@LocalVariable stress lambda; -@LocalVariable stress mu; -@LocalVariable stress Fel; -@LocalVariable StressStensor s0; - -@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); - - Fel = sqrt(J2) + beta * I1 - kappa; -} - -@ComputeStress { - sig = (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel); -} - -@Integrator { - Stensor s_dev = deviator(sig); - - if (Fel > 0) { - Stensor nq = sqrt(3.)/2./sigmaeq(sig)*s_dev; - Stensor np = 1.0 * Stensor::Id(); - - const auto yield = sqrt((s_dev|s_dev)/2.) + beta * trace(sig) - kappa; - - const auto dev_flow = 1.0; - const auto iso_flow = beta; - - 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 ; - } -} - -@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/DruckerPrager.so b/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.so deleted file mode 100755 index a311106f008e899cd4c8fe9eb3146de01d08ee50..0000000000000000000000000000000000000000 Binary files a/Tests/Data/Mechanics/Ehlers/MFront/DruckerPrager.so and /dev/null differ diff --git a/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj b/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj index 936c3b4fbad43f6ce2affa1bf0feccb39afe8c05..37a43241e6dd7abe8426e7e5b183e3fcfecbf0e7 100644 --- a/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj +++ b/Tests/Data/Mechanics/Ehlers/MFront/cube_1e0_dp.prj @@ -9,7 +9,6 @@ <integration_order>2</integration_order> <constitutive_relation> <type>MFront</type> - <library>DruckerPrager.so</library> <behaviour>DruckerPrager</behaviour> <material_properties> <material_property name="YoungModulus" parameter="E"/> diff --git a/Tests/Data/ThermoMechanics/BDT/bdt.mfront b/Tests/Data/ThermoMechanics/BDT/bdt.mfront deleted file mode 100644 index 8ba480df97ec7010735d29d06f577e49dfd3b58a..0000000000000000000000000000000000000000 --- a/Tests/Data/ThermoMechanics/BDT/bdt.mfront +++ /dev/null @@ -1,148 +0,0 @@ -@Parser Implicit; -@Behaviour BDT; -@Algorithm NewtonRaphson_NumericalJacobian; -@Theta 1.; - -@MaterialProperty real young; -young.setGlossaryName("YoungModulus"); -@MaterialProperty real nu; -nu.setGlossaryName("PoissonRatio"); -@MaterialProperty thermalexpansion alpha; -alpha.setGlossaryName("ThermalExpansion"); -@MaterialProperty real f_c; -@MaterialProperty real m_0; -@MaterialProperty real alpha_p; -@MaterialProperty real n_exp_T; -@MaterialProperty real q_h0; -@MaterialProperty real chi_h; -@MaterialProperty real alpha_d; -@MaterialProperty real h_d; -@MaterialProperty real Qact; -@MaterialProperty real A_creep; -@MaterialProperty real n_creep; -@MaterialProperty real El_1; -@MaterialProperty real El_2; -@MaterialProperty real El_3; -@MaterialProperty real at_1; -@MaterialProperty real at_2; -@MaterialProperty real at_3; - -@StateVariable real pla_mult; -@StateVariable real a; -@StateVariable real dam; -@StateVariable real Rv; -@StateVariable real triax_p; -@StateVariable real vp; - -@LocalVariable stress lambda; -@LocalVariable stress mu; -@LocalVariable stress Fel; -@LocalVariable StressStensor s0; - -// @Parameter Tref; -// Tref.setDefaultValue(293.15); - -@InitLocalVariables { - auto const T_Celsius = T - 273.15; - - const auto y2 = El_1 * power<2>(T_Celsius) + El_2 * T_Celsius + El_3; - lambda = computeLambda(young * y2, nu); - mu = computeMu(young * y2, nu); - StressStensor sigel(lambda * trace(eel + deto) * Stensor::Id() + - 2 * mu * (eel + deto)); // StresssStensor in Tutorial - const auto s_dev = deviator(sigel); - // J2 = (s_dev|s_dev)/2.; //defines double contraction - const stress seq = sigmaeq(s_dev); - const stress I1 = trace(sigel); - const auto aux = (seq + I1) / (3. * f_c); - - const auto aux_ex = (1. - 1. / n_exp_T); - const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10.), n_exp_T); - - Rv = max(at_1 * T_Celsius + at_2 * I1 / 3. + at_3, 1.0e-4); - auto Rpel = -(1. - q_h0) * pow((a + vp) / Rv, 2) + - 2. * (1. - q_h0) * (a + vp) / Rv + q_h0; - if ((a + vp) > Rv) - Rpel = 1.0; - - const auto q_h = Rpel / (pow(aux_arg, aux_ex)); - - Fel = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) + - (m_0 * aux - 1) * power<2>(q_h); -} - -@ComputeStress { - sig = (1 - dam) * (lambda * trace(eel) * Stensor::Id() + 2 * mu * eel); -} - -@Integrator { - auto const T_Celsius = T - 273.15; - - Stensor nvp = Stensor(0.); - Stensor s_dev = deviator(sig / (1 - dam)); - const stress seq = sigmaeq(s_dev); - - if (seq > 1.e-15) { - nvp = 1.5 * s_dev / seq; - } - fvp -= dt * A_creep * exp(-Qact / 8.3144598 / T_Celsius) * pow(seq, n_creep); - - if (Fel > 0) { - const real pla_mult_ = pla_mult + theta * dpla_mult; - Stensor nq = s_dev; - Stensor np = 1.0 * Stensor::Id(); - const auto aux = (seq + trace(sig / (1 - dam))) / (3. * f_c); - - const auto aux_ex = (1. - 1. / n_exp_T); - const auto aux_arg = 1 + pow(alpha_p * (T_Celsius - 10), n_exp_T); - real a_; - real dam_; - real vp_; - - a_ = a + theta * da; - dam_ = dam + theta * ddam; - vp_ = vp + theta * dvp; - Rv = max(at_1 * T_Celsius + at_2 * trace(sig / (1 - dam)) / 3. + at_3, 1.0e-4); - auto Rp_ = -(1. - q_h0) * pow((a_ + vp_) / Rv, 2) + - 2. * (1. - q_h0) * (a_ + vp_) / Rv + q_h0; - if ((a_ + vp_) > Rv) - Rp_ = 1.0; - const auto q_h = Rp_ / (pow(aux_arg, aux_ex)); - - const auto yield = power<2>((1. - q_h) * power<2>(aux) + seq / f_c) + - (m_0 * aux - 1) * power<2>(q_h); - - const auto big_aux = power<2>(aux) * (1 - q_h) + seq / f_c; - const auto dev_flow = - power<2>(q_h) * m_0 / (2. * f_c * seq) + - 2. * ((1 - q_h) * aux / (f_c * seq) + 3. / (2. * f_c * seq)) * big_aux; - const auto iso_flow = power<2>(q_h) * m_0 / (3. * f_c) + - 4. * (1. - q_h) * aux / (3. * f_c) * big_aux; - - feel = - deel - deto + dpla_mult * (nq * dev_flow + np * iso_flow) + dvp * nvp; - 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))); - triax_p = iso_flow / dev_flow; - - dam = max(min(1 - exp(-(((a + vp) - Rv) / alpha_d)), 1.0), 0.0); - sig = (1 - dam) * sig; - } else { - feel = deel - deto + dvp * nvp; - } -} - -@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/ThermoMechanics/BDT/cube_1e0_bdt.prj b/Tests/Data/ThermoMechanics/BDT/cube_1e0_bdt.prj index 0719b16c569438c1321f7112e1117f9a3560a042..6efa350a3df47546c3bad3e518134bcfda47758b 100644 --- a/Tests/Data/ThermoMechanics/BDT/cube_1e0_bdt.prj +++ b/Tests/Data/ThermoMechanics/BDT/cube_1e0_bdt.prj @@ -9,7 +9,6 @@ <integration_order>2</integration_order> <constitutive_relation> <type>MFront</type> - <library>libBehaviour.so</library> <behaviour>BDT</behaviour> <material_properties> <material_property name="YoungModulus" parameter="E"/> diff --git a/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj b/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj index 1641b3dd5698f5f99177c49bf23d5da4962c4d15..cf3dcc9c588b8ce3fc497e154377a559e8d53513 100644 --- a/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj +++ b/Tests/Data/ThermoMechanics/LinearMFront/cube_1e0_lin.prj @@ -9,7 +9,6 @@ <integration_order>2</integration_order> <constitutive_relation> <type>MFront</type> - <library>libBehaviour.so</library> <behaviour>Elasticity</behaviour> <material_properties> <material_property name="YoungModulus" parameter="E"/> diff --git a/Tests/Data/ThermoMechanics/LinearMFront/elasticity.mfront b/Tests/Data/ThermoMechanics/LinearMFront/elasticity.mfront deleted file mode 120000 index 5b5479bafded1f2fccc30bd08d5805ef3dad2aaa..0000000000000000000000000000000000000000 --- a/Tests/Data/ThermoMechanics/LinearMFront/elasticity.mfront +++ /dev/null @@ -1 +0,0 @@ -../../Mechanics/Linear/MFront/disc_with_hole/elasticity.mfront \ No newline at end of file diff --git a/Tests/Data/ThermoMechanics/LinearMFront/libBehaviour.so b/Tests/Data/ThermoMechanics/LinearMFront/libBehaviour.so deleted file mode 100755 index 93632cd232f2d286a7400d749ef622ae8b81f758..0000000000000000000000000000000000000000 Binary files a/Tests/Data/ThermoMechanics/LinearMFront/libBehaviour.so and /dev/null differ