Skip to content

Commit

Permalink
Tests for TM
Browse files Browse the repository at this point in the history
  • Loading branch information
@chleh
chleh committed Nov 7, 2018
1 parent b4539be commit 4770eb8
Show file tree
Hide file tree
Showing 17 changed files with 2,085 additions and 0 deletions.
45 changes: 45 additions & 0 deletions ProcessLib/ThermoMechanics/Tests.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,3 +113,48 @@ AddTest(
ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu epsilon epsilon 1e-15 0
ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu displacement displacement 1e-16 1e-9
)

# Basic test that MFront models work for TM.
# Linear elastic, no internal state variables, but external temperature.
AddTest(
NAME ThermoMechanics_confined_thermal_expansion_mfront
PATH ThermoMechanics/LinearMFront
EXECUTABLE ogs
EXECUTABLE_ARGS cube_1e0_lin.prj
WRAPPER time
TESTER vtkdiff
REQUIREMENTS OGS_USE_MFRONT AND NOT (OGS_USE_LIS OR OGS_USE_MPI)
DIFF_DATA
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_1_t_1.000000.vtu sigma_1 sigma 1e-16 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_2_t_2.000000.vtu sigma_2 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_3_t_3.000000.vtu sigma_3 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_4_t_4.000000.vtu sigma_4 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_5_t_5.000000.vtu sigma_5 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_6_t_6.000000.vtu sigma_6 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_7_t_7.000000.vtu sigma_7 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_8_t_8.000000.vtu sigma_8 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_9_t_9.000000.vtu sigma_9 sigma 1e-8 0
cthex_ref.vtu cube_1e0_lin_pcs_0_ts_10_t_10.000000.vtu sigma_10 sigma 1e-8 0
)

# Test of a creep law.
AddTest(
NAME ThermoMechanics_BDT_mfront
PATH ThermoMechanics/BDT
EXECUTABLE ogs
EXECUTABLE_ARGS cube_1e0_bdt.prj
WRAPPER time
TESTER vtkdiff
REQUIREMENTS OGS_USE_MFRONT AND NOT (OGS_USE_LIS OR OGS_USE_MPI)
DIFF_DATA
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_51_t_300.000000.vtu epsilon_300 epsilon 1e-8 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_51_t_300.000000.vtu sigma_300 sigma 1e-3 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_151_t_900.000000.vtu epsilon_900 epsilon 1e-9 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_151_t_900.000000.vtu sigma_900 sigma 1e-3 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_251_t_1500.000000.vtu epsilon_1500 epsilon 1e-8 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_251_t_1500.000000.vtu sigma_1500 sigma 1e-3 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_501_t_3000.000000.vtu epsilon_3000 epsilon 1e-10 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_501_t_3000.000000.vtu sigma_3000 sigma 1e-3 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_1001_t_6000.000000.vtu epsilon_6000 epsilon 1e-7 0
bdt_ref.vtu cube_1e0_bdt_pcs_0_ts_1001_t_6000.000000.vtu sigma_6000 sigma 1e-3 0
)
148 changes: 148 additions & 0 deletions Tests/Data/ThermoMechanics/BDT/bdt.mfront
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
@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;
}
}
40 changes: 40 additions & 0 deletions Tests/Data/ThermoMechanics/BDT/bdt.mtest
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
@Behaviour<generic> 'src/libBehaviour.so' 'BDT';
@MaterialProperty<constant> 'YoungModulus' 4.4619E+04;
@MaterialProperty<constant> 'PoissonRatio' 3.0000E-01;
@MaterialProperty<constant> 'ThermalExpansion' 1.0000E-05;
@MaterialProperty<constant> 'f_c' 8.0750E+02;
@MaterialProperty<constant> 'm_0' 4.3100E+00;
@MaterialProperty<constant> 'alpha_p' 1.4200E-03;
@MaterialProperty<constant> 'n_exp_T' 4.7200E+00;
@MaterialProperty<constant> 'q_h0' 8.5900E-01;
@MaterialProperty<constant> 'chi_h' 0.0010;
@MaterialProperty<constant> 'alpha_d' 2.0000E-02;
@MaterialProperty<constant> 'h_d' 1.0000E+00;
@MaterialProperty<constant> 'Qact' 5.5400E+05;
@MaterialProperty<constant> 'A_creep' 5.3600E+10;
@MaterialProperty<constant> 'n_creep' 6.8000E+00;
@MaterialProperty<constant> 'El_1' -2.1693E-06;
@MaterialProperty<constant> 'El_2' 1.5710E-03;
@MaterialProperty<constant> 'El_3' 6.8680E-01;
@MaterialProperty<constant> 'at_1' 1.1522E-04;
@MaterialProperty<constant> 'at_2' -9.4003E-05;
@MaterialProperty<constant> 'at_3' -8.1075E-02;

// @Parameter 'Tref' 293.15;

@ExternalStateVariable 'Temperature' 1173.15;

@ImposedStrain 'EYY' {
0. : 0.,
750 : -0.0075,
6.1500E+03 : -6.1500E-02,
6.5260E+03 : -8.0300E-02,
6.7740E+03 : -1.0510E-01
};

@ImposedStress 'SXX' {0. : 0., 750 : -300, 6.7740E+03 : -300};

@ImposedStress 'SZZ' {0. : 0., 750 : -300, 6.7740E+03 : -300};

// @Times{0., 6.7740E+03 in 1000};
@Times{0., 6.0E+03 in 1000};
Loading

0 comments on commit 4770eb8

Please sign in to comment.