From ae592b0aaae0fabc0140115519a7262d1ee16011 Mon Sep 17 00:00:00 2001 From: jain651 Date: Mon, 8 Mar 2021 11:01:00 -0700 Subject: [PATCH] Adding plastic_damage_model for concrete (#184) --- .../ComputeMultipleInelasticDamageStress.md | 21 + .../PlasticDamageStressUpdateOne_v3.md | 15 + .../ComputeMultipleInelasticDamageStress.h | 49 ++ .../PlasticDamageStressUpdateOne_v3.h | 197 ++++++ .../ComputeMultipleInelasticDamageStress.C | 71 ++ .../PlasticDamageStressUpdateOne_v3.C | 629 ++++++++++++++++++ .../analysis/plotting_SEM_v3.py | 159 +++++ .../input/multiaxial1_cmp_v3.i | 309 +++++++++ .../input/multiaxial1_ten_v3.i | 309 +++++++++ .../input/multiaxial2_cmp_v3.i | 314 +++++++++ .../input/multiaxial2_ten_v3.i | 315 +++++++++ .../input/shr_cyclic_v3.i | 309 +++++++++ .../input/uni_cmp_dila_ap25_v3.i | 499 ++++++++++++++ .../input/uni_cmp_dila_ap2_v3.i | 499 ++++++++++++++ .../input/uni_cmp_dila_ap3_v3.i | 499 ++++++++++++++ .../plastic_damage_model/input/uni_cmp_v3.i | 499 ++++++++++++++ .../input/uni_ten_msh_sen_1ele_v3.i | 317 +++++++++ .../input/uni_ten_msh_sen_2ele_v3.i | 367 ++++++++++ .../input/uni_ten_msh_sen_4ele_v3.i | 367 ++++++++++ .../plastic_damage_model/input/uni_ten_v3.i | 329 +++++++++ 20 files changed, 6073 insertions(+) create mode 100644 doc/content/source/materials/ComputeMultipleInelasticDamageStress.md create mode 100644 doc/content/source/materials/PlasticDamageStressUpdateOne_v3.md create mode 100644 include/materials/ComputeMultipleInelasticDamageStress.h create mode 100644 include/materials/PlasticDamageStressUpdateOne_v3.h create mode 100644 src/materials/ComputeMultipleInelasticDamageStress.C create mode 100644 src/materials/PlasticDamageStressUpdateOne_v3.C create mode 100644 test/tests/plastic_damage_model/analysis/plotting_SEM_v3.py create mode 100644 test/tests/plastic_damage_model/input/multiaxial1_cmp_v3.i create mode 100644 test/tests/plastic_damage_model/input/multiaxial1_ten_v3.i create mode 100644 test/tests/plastic_damage_model/input/multiaxial2_cmp_v3.i create mode 100644 test/tests/plastic_damage_model/input/multiaxial2_ten_v3.i create mode 100644 test/tests/plastic_damage_model/input/shr_cyclic_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_cmp_dila_ap25_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_cmp_dila_ap2_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_cmp_dila_ap3_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_cmp_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_ten_msh_sen_1ele_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_ten_msh_sen_2ele_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_ten_msh_sen_4ele_v3.i create mode 100644 test/tests/plastic_damage_model/input/uni_ten_v3.i diff --git a/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md b/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md new file mode 100644 index 00000000..6b739c75 --- /dev/null +++ b/doc/content/source/materials/ComputeMultipleInelasticDamageStress.md @@ -0,0 +1,21 @@ +# Compute Multiple Inelastic Damage Stress + +!syntax description /Materials/ComputeMultipleInelasticDamageStress + +## Description + +`ComputeMultipleInelasticDamageStress` computes the damage stress. + +## Example Input Files + +The input settings for the inelastic material model is as follows + +!listing test/tests/plastic_damage_model/input/uni_ten_v3.i block=Materials/stress + +!syntax parameters /Materials/ComputeMultipleInelasticDamageStress + +!syntax inputs /Materials/ComputeMultipleInelasticDamageStress + +!syntax children /Materials/ComputeMultipleInelasticDamageStress + +!bibtex bibliography diff --git a/doc/content/source/materials/PlasticDamageStressUpdateOne_v3.md b/doc/content/source/materials/PlasticDamageStressUpdateOne_v3.md new file mode 100644 index 00000000..34404a12 --- /dev/null +++ b/doc/content/source/materials/PlasticDamageStressUpdateOne_v3.md @@ -0,0 +1,15 @@ +# Plastic Damage Model + +!syntax description /Materials/PlasticDamageStressUpdateOne_v3 + +# Plastic damage model for concrete + +Documentation will be updated later. + +!syntax parameters /Materials/PlasticDamageStressUpdateOne_v3 + +!syntax inputs /Materials/PlasticDamageStressUpdateOne_v3 + +!syntax children /Materials/PlasticDamageStressUpdateOne_v3 + +!bibtex bibliography diff --git a/include/materials/ComputeMultipleInelasticDamageStress.h b/include/materials/ComputeMultipleInelasticDamageStress.h new file mode 100644 index 00000000..b8f3fb3f --- /dev/null +++ b/include/materials/ComputeMultipleInelasticDamageStress.h @@ -0,0 +1,49 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#ifndef COMPUTEMULTIPLEINELASTICDAMAGESTRESS_H +#define COMPUTEMULTIPLEINELASTICDAMAGESTRESS_H + +#include "ComputeMultipleInelasticStress.h" + +class ComputeMultipleInelasticDamageStress; + +template <> +InputParameters validParams(); + +/** + * ComputeMultipleInelasticDamageStress is a specialized version of + * ComputeMultipleInelasticStress for use only with the + * PlasticDamageStressUpdate model. + */ + +class ComputeMultipleInelasticDamageStress : public ComputeMultipleInelasticStress +{ +public: + ComputeMultipleInelasticDamageStress(const InputParameters & parameters); + +protected: + /// damage parameter for PlasticDamageStressUpdate model + const MaterialProperty & _D; + const MaterialProperty & _D_old; + const MaterialProperty & _D_older; + + virtual void computeQpJacobianMult() override; + + virtual void computeAdmissibleState(unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & inelastic_strain_increment, + RankFourTensor & consistent_tangent_operator) override; + + virtual void updateQpStateSingleModel(unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & combined_inelastic_strain_increment) override; +}; + +#endif // ComputeMultipleInelasticDamageStress_H diff --git a/include/materials/PlasticDamageStressUpdateOne_v3.h b/include/materials/PlasticDamageStressUpdateOne_v3.h new file mode 100644 index 00000000..6b74bd29 --- /dev/null +++ b/include/materials/PlasticDamageStressUpdateOne_v3.h @@ -0,0 +1,197 @@ +/****************************************************************/ +/* MOOSE - Multiphysics Object Oriented Simulation Environment */ +/* */ +/* All contents are licensed under LGPL V2.1 */ +/* See LICENSE for full restrictions */ +/****************************************************************/ +#ifndef PLASTICDAMAGESTRESSUPDATEONE_v3_H +#define PLASTICDAMAGESTRESSUPDATEONE_v3_H + +#include "MultiParameterPlasticityStressUpdate.h" +#include "TensorMechanicsHardeningModel.h" + +class PlasticDamageStressUpdateOne_v3; + +template <> +InputParameters validParams(); + +/** + * PlasticDamageStressUpdateOne_v3 implements rate-independent associative tensile failure + * ("Rankine" plasticity) with hardening/softening. + */ +class PlasticDamageStressUpdateOne_v3 : public MultiParameterPlasticityStressUpdate +{ +public: + PlasticDamageStressUpdateOne_v3(const InputParameters & parameters); + + /** + * Does the model require the elasticity tensor to be isotropic? + */ + bool requiresIsotropicTensor() override { return true; } + +protected: + virtual void initQpStatefulProperties() override; + virtual void finalizeReturnProcess(const RankTwoTensor & rotation_increment) override; + const Real _f_tol; + + const Real _alfa; + const Real _alfa_p; + const Real _s0; + + const Real _Chi; + const Real _Dt; + const Real _ft; + const Real _FEt; + + const Real _fyc; + const Real _Dc; + const Real _fc; + const Real _FEc; + + const Real _at; + const Real _ac; + const Real _zt; + const Real _zc; + const Real _dPhit; + const Real _dPhic; + const Real _sqrtPhit_max; + const Real _sqrtPhic_max; + const Real _dt_bt; + const Real _dc_bc; + const Real _ft0; + const Real _fc0; + const Real _small_smoother2; + + const Real _c; + const Real _eps; + const int _nGpt; + const Real _tol; + + const Real _sqrt3; + + /// Whether to provide an estimate of the returned stress, based on perfect plasticity + const bool _perfect_guess; + + /// Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to stress space + RankTwoTensor _eigvecs; + + MaterialProperty & _max_principal; + MaterialProperty & _min_principal; + MaterialProperty & _intnl0; + MaterialProperty & _intnl1; + MaterialProperty & _ele_len; + MaterialProperty & _gt; + MaterialProperty & _gc; + + MaterialProperty & _tD; + MaterialProperty & _cD; + MaterialProperty & _D; + MaterialProperty & _min_ep; + MaterialProperty & _mid_ep; + MaterialProperty & _max_ep; + MaterialProperty & _sigma0; + MaterialProperty & _sigma1; + MaterialProperty & _sigma2; + + Real ft(const std::vector & intnl) const; /// tensile strength + Real dft(const std::vector & intnl) const; /// d(tensile strength)/d(intnl) + Real fc(const std::vector & intnl) const; /// compressive strength + Real dfc(const std::vector & intnl) const; /// d(compressive strength)/d(intnl) + Real beta(const std::vector & intnl) const; + Real dbeta0(const std::vector & intnl) const; + Real dbeta1(const std::vector & intnl) const; + void weighfac(const std::vector & stress_params, Real & wf) const; /// weight factor + void dweighfac(const std::vector & stress_params, Real & wf, std::vector & dwf) const; /// d(weight factor)/d(stress) + Real damageVar(const std::vector & stress_params, const std::vector & intnl) const; + + void computeStressParams(const RankTwoTensor & stress, + std::vector & stress_params) const override; + + std::vector dstress_param_dstress(const RankTwoTensor & stress) const override; + + std::vector d2stress_param_dstress(const RankTwoTensor & stress) const override; + + void setEffectiveElasticity(const RankFourTensor & Eijkl) override; + + virtual void preReturnMapV(const std::vector & trial_stress_params, + const RankTwoTensor & stress_trial, + const std::vector & intnl_old, + const std::vector & yf, + const RankFourTensor & Eijkl) override; + + virtual void setStressAfterReturnV(const RankTwoTensor & stress_trial, + const std::vector & stress_params, + Real gaE, + const std::vector & intnl, + const yieldAndFlow & smoothed_q, + const RankFourTensor & Eijkl, + RankTwoTensor & stress) const override; + + void yieldFunctionValuesV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & yf) const override; + + void computeAllQV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & all_q) const override; + + virtual void + flowPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & r) const; + + virtual void + dflowPotential_dstress(const std::vector & stress_params, + const std::vector & intnl, + std::vector< std::vector > & dr_dstress) const; + + + virtual void + dflowPotential_dintnl(const std::vector & stress_params, + const std::vector & intnl, + std::vector< std::vector > & dr_dintnl) const; + + virtual void + hardPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & h) const; + + virtual void + dhardPotential_dstress(const std::vector & stress_params, + const std::vector & intnl, + std::vector > & dh_dsig) const; + + virtual void + dhardPotential_dintnl(const std::vector & stress_params, + const std::vector & intnl, + std::vector > & dh_dintnl) const; + + void initialiseVarsV(const std::vector & trial_stress_params, + const std::vector & intnl_old, + std::vector & stress_params, + Real & gaE, + std::vector & intnl) const; + + void setIntnlValuesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl_old, + std::vector & intnl) const override; + + void setIntnlDerivativesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl, + std::vector> & dintnl) const override; + + virtual void consistentTangentOperatorV(const RankTwoTensor & stress_trial, + const std::vector & trial_stress_params, + const RankTwoTensor & stress, + const std::vector & stress_params, + Real gaE, + const yieldAndFlow & smoothed_q, + const RankFourTensor & Eijkl, + bool compute_full_tangent_operator, + const std::vector> & dvar_dtrial, + RankFourTensor & cto) override; +}; + +#endif // PLASTICDAMAGESTRESSUPDATE_H diff --git a/src/materials/ComputeMultipleInelasticDamageStress.C b/src/materials/ComputeMultipleInelasticDamageStress.C new file mode 100644 index 00000000..efd9f392 --- /dev/null +++ b/src/materials/ComputeMultipleInelasticDamageStress.C @@ -0,0 +1,71 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "ComputeMultipleInelasticDamageStress.h" + +#include "StressUpdateBase.h" + +registerMooseObject("BlackBearApp", ComputeMultipleInelasticDamageStress); + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + return params; +} + +ComputeMultipleInelasticDamageStress::ComputeMultipleInelasticDamageStress(const InputParameters & parameters) + : ComputeMultipleInelasticStress(parameters), + _D(getMaterialProperty("Damage_Variable")), + _D_old(getMaterialPropertyOld("Damage_Variable")), + _D_older(getMaterialPropertyOlder("Damage_Variable")) +{ +} + +void +ComputeMultipleInelasticDamageStress::computeQpJacobianMult() +{ + ComputeMultipleInelasticStress::computeQpJacobianMult(); + _Jacobian_mult[_qp] = (1.0 - _D_older[_qp]) * _Jacobian_mult[_qp]; + // _Jacobian_mult[_qp] = (1.0 - _D[_qp]) * _Jacobian_mult[_qp]; +} + +void +ComputeMultipleInelasticDamageStress::updateQpStateSingleModel( + unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & combined_inelastic_strain_increment) +{ + ComputeMultipleInelasticStress::updateQpStateSingleModel(model_number, + elastic_strain_increment, + combined_inelastic_strain_increment); + _Jacobian_mult[_qp] = (1.0 - _D_older[_qp]) * _Jacobian_mult[_qp]; + // _Jacobian_mult[_qp] = (1.0 - _D[_qp]) * _Jacobian_mult[_qp]; +} + +void +ComputeMultipleInelasticDamageStress::computeAdmissibleState(unsigned model_number, + RankTwoTensor & elastic_strain_increment, + RankTwoTensor & inelastic_strain_increment, + RankFourTensor & consistent_tangent_operator) +{ + _models[model_number]->updateState(elastic_strain_increment, + inelastic_strain_increment, + _rotation_increment[_qp], + _stress[_qp], + _stress_old[_qp] / (1.0 - _D_older[_qp]), + // _stress_old[_qp] / (1.0 - _D[_qp]), + _elasticity_tensor[_qp], + _elastic_strain_old[_qp], + _tangent_operator_type == TangentOperatorEnum::nonlinear, + consistent_tangent_operator); + // _stress[_qp] *= (1.0 - _D[_qp]); + _stress[_qp] *= (1.0 - _D_older[_qp]); +} diff --git a/src/materials/PlasticDamageStressUpdateOne_v3.C b/src/materials/PlasticDamageStressUpdateOne_v3.C new file mode 100644 index 00000000..45009a35 --- /dev/null +++ b/src/materials/PlasticDamageStressUpdateOne_v3.C @@ -0,0 +1,629 @@ +/****************************************************************/ +/* MOOSE - Multiphysics Object Oriented Simulation Environment */ +/* */ +/* All contents are licensed under LGPL V2.1 */ +/* See LICENSE for full restrictions */ +/****************************************************************/ +#include "PlasticDamageStressUpdateOne_v3.h" +#include "libmesh/utility.h" + +registerMooseObject("BlackBearApp", PlasticDamageStressUpdateOne_v3); + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params.addParam("yield_function_tolerance", "If the yield function is less than this amount, the (stress, internal parameters) are deemed admissible. A std::vector of tolerances must be entered for the multi-surface case"); + + params.addRangeCheckedParam("factor_relating_biaxial_unixial_cmp_str", + 0.1, + "factor_relating_biaxial_unixial_cmp_str < 0.5 & factor_relating_biaxial_unixial_cmp_str >= 0", + "Material parameter that relate biaxial and uniaxial compressive strength, i.e., \alfa = (fb0-fc0)/(2*fb0-fc0)"); + params.addRequiredParam("factor_controlling_dilatancy", + "parameter for the dilation"); + params.addRangeCheckedParam("stiff_recovery_factor", + 0., + "stiff_recovery_factor <= 1. & stiff_recovery_factor >= 0", + "stiffness recovery parameter"); + + params.addRangeCheckedParam("ft_ep_slope_factor_at_zero_ep", + "ft_ep_slope_factor_at_zero_ep <= 1 & ft_ep_slope_factor_at_zero_ep >= 0", + "slope of ft vs plastic strain curve at zero plastic strain"); + params.addRequiredParam("damage_half_ten_str", + "Fraction of the elastic recovery slope in tension at 0.5*ft0 after yielding"); + params.addRangeCheckedParam("yield_ten_str", + "yield_ten_str >= 0", + "Tensile yield strength of concrete"); + params.addRangeCheckedParam("frac_energy_ten", + "frac_energy_ten >= 0", + "Fracture energy of concrete in uniaxial tension"); + + params.addRangeCheckedParam("yield_cmp_str", + "yield_cmp_str >= 0", + "Absolute yield compressice strength"); + params.addRequiredParam("damage_max_cmp_str", + "damage at maximum compressive strength"); + params.addRequiredParam("max_cmp_str", + "Absolute maximum compressive strength"); + params.addRangeCheckedParam("frac_energy_cmp", + "frac_energy_cmp >= 0", + "Fracture energy of concrete in uniaxial compression"); + + + params.addRequiredRangeCheckedParam("tip_smoother", + "tip_smoother>=0", + "Smoothing parameter: the cone vertex at mean = cohesion*cot(friction_angle), will be smoothed by the given amount. Typical value is 0.1*cohesion"); + params.addParam("perfect_guess", + true, + "Provide a guess to the Newton-Raphson proceedure " + "that is the result from perfect plasticity. With " + "severe hardening/softening this may be " + "suboptimal."); + params.addClassDescription( + "Plastic damage model for concrete"); + return params; +} + +PlasticDamageStressUpdateOne_v3::PlasticDamageStressUpdateOne_v3(const InputParameters & parameters) + : MultiParameterPlasticityStressUpdate(parameters, 3, 1, 2), + _f_tol(getParam("yield_function_tol")), + + _alfa(getParam("factor_relating_biaxial_unixial_cmp_str")), + _alfa_p(getParam("factor_controlling_dilatancy")), + _s0(getParam("stiff_recovery_factor")), + + _Chi(getParam("ft_ep_slope_factor_at_zero_ep")), + _Dt(getParam("damage_half_ten_str")), + _ft(getParam("yield_ten_str")), + _FEt(getParam("frac_energy_ten")), + + _fyc(getParam("yield_cmp_str")), + _Dc(getParam("damage_max_cmp_str")), + _fc(getParam("max_cmp_str")), + _FEc(getParam("frac_energy_cmp")), + + _at(1.5*std::sqrt(1-_Chi) - 0.5), + _ac((2.*(_fc/_fyc)-1.+2.*std::sqrt(std::pow((_fc/_fyc),2.)-_fc/_fyc))), + + _zt((1.+_at)/_at), + _zc((1.+_ac)/_ac), + _dPhit(_at*(2.+_at)), + _dPhic(_ac*(2.+_ac)), + _sqrtPhit_max((1.+_at+sqrt(1.+_at*_at))/2.), + _sqrtPhic_max((1.+_ac )/2.), + _dt_bt(log(1.-_Dt)/log((1.+_at-sqrt(1.+_at*_at))/(2.*_at))), + _dc_bc(log(1.-_Dc)/log((1.+_ac )/(2.*_ac))), + _ft0(0.5*_ft/((1.-_Dt)*pow((_zt-_sqrtPhit_max/_at),(1.-_dt_bt))*_sqrtPhit_max)), + _fc0(_fc/((1.-_Dc)*pow((_zc-_sqrtPhic_max/_ac),(1.-_dc_bc))*_sqrtPhic_max)), + _small_smoother2(std::pow(getParam("tip_smoother"), 2)), + + _c(2.2523), + _eps(1.E-6), + _nGpt(12), + _tol(1.E-3), + + _sqrt3(sqrt(3.)), + _perfect_guess(getParam("perfect_guess")), + _eigvecs(RankTwoTensor()), + _max_principal(declareProperty("max_principal_stress")), + _min_principal(declareProperty("min_principal_stress")), + _intnl0(declareProperty("damage_parameter_for_tension")), + _intnl1(declareProperty("damage_parameter_for_compression")), + _ele_len (declareProperty("element_length")), + _gt(declareProperty("fracture energy in tension for the element")), + _gc(declareProperty("fracture energy in compression for the element")), + _tD(declareProperty("tensile_damage")), + _cD(declareProperty("compression_damage")), + _D(declareProperty("Damage_Variable")), + _min_ep(declareProperty("min_ep")), + _mid_ep(declareProperty("mid_ep")), + _max_ep(declareProperty("max_ep")), + _sigma0(declareProperty("damaged_min_principal_stress")), + _sigma1(declareProperty("damaged_mid_principal_stress")), + _sigma2(declareProperty("damaged_max_principal_stress")) +{ +} + +void +PlasticDamageStressUpdateOne_v3::initQpStatefulProperties() +{ + if (_current_elem->n_vertices() < 3) + _ele_len[_qp] = _current_elem -> length (0,1); + else if (_current_elem->n_vertices() < 5) + _ele_len[_qp] = (_current_elem -> length (0,1) + + _current_elem -> length (1,2))/2.; + else + _ele_len[_qp] = (_current_elem -> length (0,1) + + _current_elem -> length (1,2) + + _current_elem -> length (0,4))/3.; + + _gt[_qp] = _FEt/_ele_len[_qp]; + _gc[_qp] = _FEc/_ele_len[_qp]; + + _min_ep[_qp] = 0.; + _mid_ep[_qp] = 0.; + _max_ep[_qp] = 0.; + _sigma0[_qp] = 0.; + _sigma1[_qp] = 0.; + _sigma2[_qp] = 0.; + _intnl0[_qp] = 0.; + _intnl1[_qp] = 0.; + _tD[_qp] = 0.; + _cD[_qp] = 0.; + _D[_qp] = 0.; + MultiParameterPlasticityStressUpdate::initQpStatefulProperties(); +} + +void +PlasticDamageStressUpdateOne_v3::finalizeReturnProcess( + const RankTwoTensor & /*rotation_increment*/) +{ + std::vector eigstrain; + _plastic_strain[_qp].symmetricEigenvalues(eigstrain); + _min_ep[_qp] = eigstrain[0]; + _mid_ep[_qp] = eigstrain[1]; + _max_ep[_qp] = eigstrain[2]; +} + +void +PlasticDamageStressUpdateOne_v3::computeStressParams(const RankTwoTensor & stress, + std::vector & stress_params) const +{ + stress.symmetricEigenvalues(stress_params); +} + +std::vector +PlasticDamageStressUpdateOne_v3::dstress_param_dstress(const RankTwoTensor & stress) const +{ + std::vector sp; + std::vector dsp; + stress.dsymmetricEigenvalues(sp, dsp); + return dsp; +} + +std::vector +PlasticDamageStressUpdateOne_v3::d2stress_param_dstress(const RankTwoTensor & stress) const +{ + std::vector d2; + stress.d2symmetricEigenvalues(d2); + return d2; +} + +void +PlasticDamageStressUpdateOne_v3::setEffectiveElasticity(const RankFourTensor & Eijkl) +{ + // Eijkl is required to be isotropic, so we can use the + // frame where stress is diagonal + for (unsigned a = 0; a < _num_sp; ++a) + for (unsigned b = 0; b < _num_sp; ++b) + _Eij[a][b] = Eijkl(a, a, b, b); + _En = _Eij[2][2]; + const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]); + for (unsigned a = 0; a < _num_sp; ++a) + { + _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom; + for (unsigned b = 0; b < a; ++b) + _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom; + } +} + +void +PlasticDamageStressUpdateOne_v3::preReturnMapV(const std::vector & /*trial_stress_params*/, + const RankTwoTensor & stress_trial, + const std::vector & /*intnl_old*/, + const std::vector & /*yf*/, + const RankFourTensor & /*Eijkl*/) +{ + std::vector eigvals; + stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs); +} + +void +PlasticDamageStressUpdateOne_v3::setStressAfterReturnV(const RankTwoTensor & /*stress_trial*/, + const std::vector & stress_params, + Real /*gaE*/, + const std::vector & intnl, + const yieldAndFlow & /*smoothed_q*/, + const RankFourTensor & /*Eijkl*/, + RankTwoTensor & stress) const +{ + // form the diagonal stress + stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0); + // rotate to the original frame + stress = _eigvecs * stress * (_eigvecs.transpose()); +// _dir[_qp] = _eigvecs; + Real D = damageVar(stress_params, intnl); + _sigma0[_qp] = (1.-D)*stress_params[0]; + _sigma1[_qp] = (1.-D)*stress_params[1]; + _sigma2[_qp] = (1.-D)*stress_params[2]; + _intnl0[_qp] = intnl[0]; + _intnl1[_qp] = intnl[1]; + _D[_qp] = D; +} + +void +PlasticDamageStressUpdateOne_v3::yieldFunctionValuesV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & yf) const +{ + Real I1 = stress_params[0] + stress_params[1] + stress_params[2]; + Real J2 = (pow(stress_params[0]-stress_params[1],2.) + + pow(stress_params[1]-stress_params[2],2.) + + pow(stress_params[2]-stress_params[0],2.))/6.+_small_smoother2; + Real sqrtJ2 = sqrt(J2); + yf[0] = 1./(1.-_alfa)*(_alfa*I1 + _sqrt3*sqrtJ2 + + beta(intnl)*(stress_params[2]<0.?0.:stress_params[2])) - fc(intnl); +} + +void +PlasticDamageStressUpdateOne_v3::computeAllQV(const std::vector & stress_params, + const std::vector & intnl, + std::vector & all_q) const +{ + Real I1 = stress_params[0] + stress_params[1] + stress_params[2]; + Real J2 = (pow(stress_params[0]-stress_params[1],2.) + + pow(stress_params[1]-stress_params[2],2.) + + pow(stress_params[2]-stress_params[0],2.))/6.+_small_smoother2; + Real sqrtJ2 = sqrt(J2); + std::vector DevSt(3); // vector of principal deviatoric stress + for (unsigned i = 0; i < 3; ++i) + DevSt[i] = stress_params[i]-I1/3.; + +// yieldFunctionValuesV(stress_params, intnl, all_q[0].f); + all_q[0].f = 1./(1.-_alfa)*(_alfa*I1 + _sqrt3*sqrtJ2 + + beta(intnl)*(stress_params[2]<0.?0.:stress_params[2])) - fc(intnl); + + for (unsigned i = 0; i < _num_sp; ++i) + all_q[0].df[i] = 1./(1.-_alfa)*(_alfa + _sqrt3*DevSt[i]/(2.*sqrtJ2) + beta(intnl)*(stress_params[2]<0.?0.:(i==2))); + all_q[0].df_di[0] = 1./(1.-_alfa)*(dbeta0(intnl)*(stress_params[2]<0.?0.:stress_params[2])); + all_q[0].df_di[1] = 1./(1.-_alfa)*(dbeta1(intnl)*(stress_params[2]<0.?0.:stress_params[2])) - dfc(intnl); + + flowPotential(stress_params, intnl, all_q[0].dg); + dflowPotential_dstress(stress_params, intnl, all_q[0].d2g); + dflowPotential_dintnl(stress_params, intnl, all_q[0].d2g_di); +} + +void +PlasticDamageStressUpdateOne_v3::flowPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & r) const +{ + Real J2 = (pow(stress_params[0]-stress_params[1],2.) + + pow(stress_params[1]-stress_params[2],2.) + + pow(stress_params[2]-stress_params[0],2.))/6.+_small_smoother2; + Real invsqrt2J2 = 1./sqrt(2.*J2); + std::vector DevSt(3); + DevSt[0] = (2.*stress_params[0]-stress_params[1]-stress_params[2])/3.; // dJ2/dsig0 + DevSt[1] = (2.*stress_params[1]-stress_params[2]-stress_params[0])/3.; // dJ2/dsig1 + DevSt[2] = (2.*stress_params[2]-stress_params[0]-stress_params[1])/3.; // dJ2/dsig2 + + Real D = damageVar(stress_params, intnl); + + for (unsigned int i = 0; i < _num_sp; ++i) + r[i] = (_alfa_p + (J2<_f_tol ? 0. : DevSt[i]*invsqrt2J2))*pow((1.-D),1); +} + +void +PlasticDamageStressUpdateOne_v3::dflowPotential_dstress(const std::vector & stress_params, + const std::vector & intnl, + std::vector< std::vector > & dr_dstress) const +{ + Real J2 = (pow(stress_params[0]-stress_params[1],2.) + + pow(stress_params[1]-stress_params[2],2.) + + pow(stress_params[2]-stress_params[0],2.))/6.+_small_smoother2; + Real invsqrt2J2 = 1./sqrt(2.*J2); + std::vector DevSt(3); + DevSt[0] = (2.*stress_params[0]-stress_params[1]-stress_params[2])/3.; // dJ2/dsig0 + DevSt[1] = (2.*stress_params[1]-stress_params[2]-stress_params[0])/3.; // dJ2/dsig1 + DevSt[2] = (2.*stress_params[2]-stress_params[0]-stress_params[1])/3.; // dJ2/dsig2 + + Real D = damageVar(stress_params, intnl); + + for (unsigned i = 0; i < _num_sp; ++i) + for (unsigned j = 0; j < (i+1); ++j) + { + if (i!=j) + { + dr_dstress[i][j] = J2<_f_tol ? 0. : invsqrt2J2*(-1./3.-DevSt[i]*DevSt[j]/(2.*J2))*pow((1.-D),2); + dr_dstress[j][i] = dr_dstress[i][j]; + } + else + dr_dstress[i][i] = J2<_f_tol ? 0. : invsqrt2J2*(2./3.-DevSt[i]*DevSt[j]/(2.*J2))*pow((1.-D),2); + } +} + +void +PlasticDamageStressUpdateOne_v3::dflowPotential_dintnl(const std::vector & /* stress_params */, + const std::vector & /* intnl */, + std::vector > & dr_dintnl) const +{ + for (unsigned i = 0; i < _num_sp; ++i) + for (unsigned j = 0; j < _num_intnl; ++j) + dr_dintnl[i][j] = 0.; +} + +void +PlasticDamageStressUpdateOne_v3::hardPotential(const std::vector & stress_params, + const std::vector & intnl, + std::vector & h) const +{ + Real wf; + weighfac(stress_params, wf); + std::vector r(3); + flowPotential(stress_params, intnl, r); + h[0] = wf*ft(intnl)/_gt[_qp]*r[2]; + h[1] = -(1.-wf)*fc(intnl)/_gc[_qp]*r[0]; +} + +void +PlasticDamageStressUpdateOne_v3::dhardPotential_dstress(const std::vector & stress_params, + const std::vector & intnl, + std::vector > & dh_dsig) const +{ + Real wf; + std::vector dwf(3); + dweighfac(stress_params, wf, dwf); + + std::vector r(3); + flowPotential(stress_params, intnl, r); + std::vector > dr_dsig(3, std::vector(3)); + dflowPotential_dstress(stress_params, intnl, dr_dsig); + + for (unsigned i = 0; i < _num_sp; ++i) + { + dh_dsig[0][i] = (wf*dr_dsig[2][i] + dwf[i]*r[2])*ft(intnl)/_gt[_qp]; + dh_dsig[1][i] = -((1.-wf)*dr_dsig[0][i] - dwf[i]*r[0])*fc(intnl)/_gc[_qp]; + } +} + +void +PlasticDamageStressUpdateOne_v3::dhardPotential_dintnl(const std::vector & stress_params, + const std::vector & intnl, + std::vector > & dh_dintnl) const +{ + Real wf; + weighfac(stress_params, wf); + std::vector r(3); + flowPotential(stress_params, intnl, r); + + dh_dintnl[0][0] = wf*dft(intnl)/_gt[_qp]*r[2]; + dh_dintnl[0][1] = 0.; + dh_dintnl[1][0] = 0.; + dh_dintnl[1][1] = -(1-wf)*dfc(intnl)/_gc[_qp]*r[0]; +} + +void +PlasticDamageStressUpdateOne_v3::initialiseVarsV(const std::vector & trial_stress_params, + const std::vector & intnl_old, + std::vector & stress_params, + Real & /* gaE */, + std::vector & intnl) const +{ + setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl); +} + +void +PlasticDamageStressUpdateOne_v3::setIntnlValuesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl_old, + std::vector & intnl) const +{ + Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2]; + Real J2_trial = (pow(trial_stress_params[0]-trial_stress_params[1],2.) + + pow(trial_stress_params[1]-trial_stress_params[2],2.) + + pow(trial_stress_params[2]-trial_stress_params[0],2.))/6.+_small_smoother2; + Real invsqrt2J2_trial = 1./sqrt(2.*J2_trial); + Real G = 0.5*(_Eij[0][0]-_Eij[0][1]); // Lame's mu + Real K = _Eij[0][1]+2.*G/3.; // Bulk modulus + Real C1= (2.*G*invsqrt2J2_trial); + Real C2 = -(I1_trial/3.*G*invsqrt2J2_trial - 3.*K*_alfa_p); + Real C3 = 3.*K*_alfa_p; + + RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0]-current_stress_params[0], + trial_stress_params[1]-current_stress_params[1], + trial_stress_params[2]-current_stress_params[2],0.,0.,0.); + RankTwoTensor fac = J2_trial< _f_tol ? + C3*RankTwoTensor(1.,1.,1.,0.,0.,0.) : + RankTwoTensor(C1*trial_stress_params[0]-C2, + C1*trial_stress_params[1]-C2, + C1*trial_stress_params[2]-C2,0.,0.,0.); + + Real lam = dsig.L2norm() / fac.L2norm(); + std::vector h(2); + hardPotential(current_stress_params, intnl_old, h); + + intnl[0] = intnl_old[0] + lam*h[0]; + intnl[1] = intnl_old[1] + lam*h[1]; +} + +void +PlasticDamageStressUpdateOne_v3::setIntnlDerivativesV(const std::vector & trial_stress_params, + const std::vector & current_stress_params, + const std::vector & intnl, + std::vector> & dintnl) const +{ + Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2]; + Real J2_trial = (pow(trial_stress_params[0]-trial_stress_params[1],2.) + + pow(trial_stress_params[1]-trial_stress_params[2],2.) + + pow(trial_stress_params[2]-trial_stress_params[0],2.))/6.; + Real invsqrt2J2_trial = 1./sqrt(2.*J2_trial); + Real G = 0.5*(_Eij[0][0]-_Eij[0][1]); // Lame's mu + Real K = _Eij[0][1]+2.*G/3.; // Bulk modulus + Real C1= (2.*G*invsqrt2J2_trial); + Real C2 = -(I1_trial/3.*G*invsqrt2J2_trial - 3.*K*_alfa_p); + Real C3 = 3.*K*_alfa_p; + + RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0]-current_stress_params[0], + trial_stress_params[1]-current_stress_params[1], + trial_stress_params[2]-current_stress_params[2],0.,0.,0.); + RankTwoTensor fac = J2_trial< _f_tol ? + C3*RankTwoTensor(1.,1.,1.,0.,0.,0.) : + RankTwoTensor(C1*trial_stress_params[0]-C2, + C1*trial_stress_params[1]-C2, + C1*trial_stress_params[2]-C2,0.,0.,0.); + + Real lam = dsig.L2norm() / fac.L2norm(); + + std::vector dlam_dsig(3); + for (unsigned i = 0; i < _num_sp; ++i) + dlam_dsig[i] = dsig.L2norm() == 0.? 0. : + -(trial_stress_params[i]-current_stress_params[i])/(dsig.L2norm() * fac.L2norm()); + + std::vector h(2); + hardPotential(current_stress_params, intnl, h); + std::vector > dh_dsig(2, std::vector(3)); + dhardPotential_dstress(current_stress_params, intnl, dh_dsig); + std::vector > dh_dintnl(2, std::vector(2)); + dhardPotential_dintnl(current_stress_params, intnl, dh_dintnl); + + for (unsigned i = 0; i < _num_intnl; ++i) + for (unsigned j = 0; j < _num_sp; ++j) + dintnl[i][j] = dlam_dsig[j]*h[i] + lam*dh_dsig[i][j]; +} + +Real +PlasticDamageStressUpdateOne_v3::ft(const std::vector & intnl) const +{ + Real sqrtPhi_t = sqrt(1.+_at*(2.+_at)*intnl[0]); + if (_zt > sqrtPhi_t/_at) + return _ft0 * pow(_zt-sqrtPhi_t/_at,(1.-_dt_bt))*sqrtPhi_t; + else + return _ft0 * 1.E-6; +} + +Real +PlasticDamageStressUpdateOne_v3::dft(const std::vector & intnl) const +{ + Real sqrtPhi_t = sqrt(1.+_at*(2.+_at)*intnl[0]); + if (_zt > sqrtPhi_t/_at) + return _ft0*_dPhit/(2*sqrtPhi_t)*pow(_zt-sqrtPhi_t/_at,-_dt_bt)*(_zt-(2.-_dt_bt)*sqrtPhi_t/_at); + else + return 0.; + +} + +Real +PlasticDamageStressUpdateOne_v3::fc(const std::vector & intnl) const +{ + Real sqrtPhi_c; + if(intnl[1]<1.0) + sqrtPhi_c = sqrt(1.+_ac*(2.+_ac)*intnl[1]); + else + sqrtPhi_c = sqrt(1.+_ac*(2.+_ac)*0.99); + return _fc0*pow((_zc-sqrtPhi_c/_ac),(1.-_dc_bc))*sqrtPhi_c; +} + +Real +PlasticDamageStressUpdateOne_v3::dfc(const std::vector & intnl) const +{ + if(intnl[1]<1.0) + { + Real sqrtPhi_c = sqrt(1.+_ac*(2.+_ac)*intnl[1]); + return _fc0*_dPhic/(2.*sqrtPhi_c)*pow(_zc-sqrtPhi_c/_ac,-_dc_bc)*(_zc-(2.-_dc_bc)*sqrtPhi_c/_ac); + } + else + return 0.; +} + +Real +PlasticDamageStressUpdateOne_v3::beta(const std::vector & intnl) const +{ + return (1.-_alfa)*fc(intnl)/ft(intnl) - (1.+_alfa); +} + +Real +PlasticDamageStressUpdateOne_v3::dbeta0(const std::vector & intnl) const +{ + return -(1.-_alfa)*fc(intnl)*dft(intnl)/pow(ft(intnl),2.); +} + +Real +PlasticDamageStressUpdateOne_v3::dbeta1(const std::vector & intnl) const +{ + return dfc(intnl)/ft(intnl)*(1.-_alfa); +} + +void +PlasticDamageStressUpdateOne_v3::weighfac(const std::vector & stress_params, Real & wf) const +{ + Real Dr = 0.; + Real Nr = 0.; + for (unsigned i = 0; i < _num_sp; ++i) + { + if (stress_params[i]>0.) + { + Nr += stress_params[i]; + Dr += stress_params[i]; + } + else + Dr += -stress_params[i]; + } + wf = Nr/Dr; +} + +void +PlasticDamageStressUpdateOne_v3::dweighfac(const std::vector & stress_params, Real & wf, std::vector & dwf) const +{ + std::vector dNr(3, 0.), dDr(3, 0.); + Real Dr = 0.; + Real Nr = 0.; + for (unsigned i = 0; i < _num_sp; ++i) + { + if (stress_params[i]>0.) + { + Nr += stress_params[i]; + dNr[i] = 1.; + Dr += stress_params[i]; + dDr[i] = 1.; + } + else + { + Dr += -stress_params[i]; + dDr[i] = -1.; + } + } + wf = Nr/Dr; + + for (unsigned i = 0; i < _num_sp; ++i) + dwf[i] = (dNr[i]-wf*dDr[i])/Dr; +} + +Real +PlasticDamageStressUpdateOne_v3::damageVar(const std::vector &stress_params, const std::vector & intnl) const +{ + Real sqrtPhi_t = sqrt(1.+_at*(2.+_at)*intnl[0]); + if (_zt > sqrtPhi_t/_at) + _tD[_qp] = 1. - pow(_zt-sqrtPhi_t/_at,_dt_bt); + else + _tD[_qp] = 1. - 1.E-6; + + Real wf; + weighfac(stress_params,wf); + Real s = _s0 + (1.-_s0)*wf; + + Real sqrtPhi_c; + if(intnl[1]<1.0) + sqrtPhi_c = sqrt(1.+_ac*(2.+_ac)*intnl[1]); + else + sqrtPhi_c = sqrt(1.+_ac*(2.+_ac)*0.99); + + _cD[_qp] = 1. - pow((_zc-sqrtPhi_c/_ac),_dc_bc); + return 1. - (1.-s*_tD[_qp])*(1.-_cD[_qp]); +} + +void +PlasticDamageStressUpdateOne_v3::consistentTangentOperatorV(const RankTwoTensor & /* stress_trial */, + const std::vector & /* trial_stress_params */, + const RankTwoTensor & /*stress*/, + const std::vector & /* stress_params */, + Real /*gaE*/, + const yieldAndFlow & /*smoothed_q*/, + const RankFourTensor & elasticity_tensor, + bool /* compute_full_tangent_operator */, + const std::vector> & /* dvar_dtrial */, + RankFourTensor & cto) +{ + cto = elasticity_tensor; + return; +} diff --git a/test/tests/plastic_damage_model/analysis/plotting_SEM_v3.py b/test/tests/plastic_damage_model/analysis/plotting_SEM_v3.py new file mode 100644 index 00000000..775d068b --- /dev/null +++ b/test/tests/plastic_damage_model/analysis/plotting_SEM_v3.py @@ -0,0 +1,159 @@ +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import matplotlib.image as mpimg +from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator) +from matplotlib.offsetbox import TextArea, DrawingArea, OffsetImage, AnnotationBbox, AnchoredOffsetbox +from numpy import * +import numpy as np +import pandas as pd +import math + +# Tension test +uni_ten_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_ten_v3.csv') + +multiaxial1_ten_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/multiaxial1_ten_v3.csv') +multiaxial2_ten_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/multiaxial2_ten_v3.csv') + +# Mesh sensitivity test +uni_ten_msh_sen_1ele_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_ten_msh_sen_1ele_v3.csv') +uni_ten_msh_sen_2ele_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_ten_msh_sen_2ele_v3.csv') +uni_ten_msh_sen_4ele_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_ten_msh_sen_4ele_v3.csv') + +# Compression test +uni_cmp_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_cmp_v3.csv') +# bi_cmp_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/bi_cmp_v3.csv') + +multiaxial1_cmp_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/multiaxial1_cmp_v3.csv') +multiaxial2_cmp_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/multiaxial2_cmp_v3.csv') + +# Dilatancy test +uni_cmp_dila_ap_02_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_cmp_dila_ap_02_v3.csv') +uni_cmp_dila_ap_025_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_cmp_dila_ap_025_v3.csv') +uni_cmp_dila_ap_03_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/uni_cmp_dila_ap_03_v3.csv') + +# Calculation of volumetric strain for Dilatancy test +vol_strain_ap_02_v3 = uni_cmp_dila_ap_02_v3['e_xx'] + uni_cmp_dila_ap_02_v3['e_yy'] + uni_cmp_dila_ap_02_v3['e_zz'] + uni_cmp_dila_ap_02_v3['max_ep'] + uni_cmp_dila_ap_02_v3['mid_ep'] + uni_cmp_dila_ap_02_v3['min_ep'] +vol_strain_ap_025_v3 = uni_cmp_dila_ap_025_v3['e_xx'] + uni_cmp_dila_ap_025_v3['e_yy'] + uni_cmp_dila_ap_025_v3['e_zz'] + uni_cmp_dila_ap_025_v3['max_ep'] + uni_cmp_dila_ap_025_v3['mid_ep'] + uni_cmp_dila_ap_025_v3['min_ep'] +vol_strain_ap_03_v3 = uni_cmp_dila_ap_03_v3['e_xx'] + uni_cmp_dila_ap_03_v3['e_yy'] + uni_cmp_dila_ap_03_v3['e_zz'] + uni_cmp_dila_ap_03_v3['max_ep'] + uni_cmp_dila_ap_03_v3['mid_ep'] + uni_cmp_dila_ap_03_v3['min_ep'] + +# Shear Cyclic Loading +shr_cyclic_v3 = pd.read_csv('./test/tests/plastic_damage_model/output/shr_cyclic_v3.csv') + +# Lee results +Lee_uni_ten = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_uni_ten.csv', header=None) +Lee_uni_cmp = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_uni_cmp.csv', header=None) +Lee_bi_ten = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_bi_ten.csv', header=None) +Lee_bi_cmp = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_bi_cmp.csv', header=None) +Lee_cyc_uni_ten = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_cyc_uni_ten.csv', header=None) +Lee_cyc_uni_cmp = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_uni_cmp.csv', header=None) +Lee_cyc_uni_ten_cmp = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_cyc_uni_ten_cmp.csv', header=None) +Lee_cyc_shr = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_cyc_shr.csv', header=None) +Lee_dila_ap_02 = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_dila_ap_02.csv', header=None) +Lee_dila_ap_025 = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_dila_ap_025.csv', header=None) +Lee_dila_ap_03 = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_dila_ap_03.csv', header=None) +Lee_uni_ten_msh_sen_2_elem = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_uni_ten_msh_sen_2_elem.csv', header=None) +Lee_uni_ten_msh_sen_4_elem = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_uni_ten_msh_sen_4_elem.csv', header=None) +Lee_uni_ten_msh_sen_8_elem = pd.read_csv('./test/tests/plastic_damage_model/analysis/Lee_uni_ten_msh_sen_8_elem.csv', header=None) + +# function to put image +def place_image(im, loc=3, ax=None, zoom=1, **kw): + if ax==None: ax=plt.gca() + imagebox = OffsetImage(im, zoom=zoom*0.72) + ab = AnchoredOffsetbox(loc=loc, child=imagebox, frameon=False, **kw) + ax.add_artist(ab) + +## plotting Single Element Model (SEM) response +############################################################################# +SEM, SEM_ax = plt.subplots(2,2, sharex=False, sharey=False, figsize = (10,8)) +plt.subplots_adjust(left=0.08, bottom=0.1, right=0.95, top=0.95, wspace = 0.25, hspace = 0.3) +# SEM.text(0.03, 0.5, 'Force [N]', va='center', rotation='vertical') +# SEM.text(0.5, 0.05, 'Deformation [mm]', ha='center') + +[i, j] = [0, 0] +SEM_ax[i][j].set_title('Multi-axial Tensile Loading', fontsize=11) +SEM_ax[i][j].plot(Lee_uni_ten[0], Lee_uni_ten[1], c = '0.6', linestyle ='solid', linewidth=1.0, label= 'Uniaxial (Lee)') +SEM_ax[i][j].plot(Lee_bi_ten[0] , Lee_bi_ten[1], c = '0.6', linestyle ='dashdot', linewidth=1.0, label= 'Biaxial (Lee)') +SEM_ax[i][j].plot(multiaxial1_ten_v3['displacement_x']/25.4, multiaxial1_ten_v3['s_xx'], c = '0', linestyle ='solid', linewidth=1.0, label= 'Uniaxial_v3') +SEM_ax[i][j].plot(multiaxial2_ten_v3['displacement_x']/25.4, multiaxial2_ten_v3['s_xx'], c = '0', linestyle ='dashdot', linewidth=1.0, label= 'Biaxial_v3') + +SEM_ax[i][j].grid(color="grey", ls = '--', lw = 0.5) +# SEM_ax[i][j].xaxis.set_minor_locator(MultipleLocator(0.2)) +# SEM_ax[i][j].yaxis.set_minor_locator(MultipleLocator(0.2)) +SEM_ax[i][j].set_xlabel("Strain [m/m]") +SEM_ax[i][j].set_ylabel("Stress [MPa]") +# SEM_ax[i][j].set_xlim(0, 0.001) +# SEM_ax[i][j].set_ylim(0,0.15) +# STT_image = mpimg.imread('./test/tests/plastic_damage_model_2017/analysis/STT.png') +# place_image(STT_image, loc='lower right', ax=SEM_ax[i][j], pad=0, zoom=0.18) +SEM_ax[i][j].legend(loc='upper right', facecolor = 'white', framealpha = 0.4, edgecolor ='none') + + +[i, j] = [0, 1] +SEM_ax[i][j].set_title('Multi-axial Compressive Loading', fontsize=11) +SEM_ax[i][j].plot(-1*Lee_uni_cmp[0], -1*Lee_uni_cmp[1], c = '0.6', linestyle ='solid', linewidth=1.0, label= 'Uniaxial (Lee)') +SEM_ax[i][j].plot(-1*Lee_bi_cmp[0] , -1*Lee_bi_cmp[1], c = '0.6', linestyle ='dashdot', linewidth=1.0, label= 'Biaxial (Lee)') +SEM_ax[i][j].plot(-1*multiaxial1_cmp_v3['displacement_x']/25.4, -1*multiaxial1_cmp_v3['s_xx'], c = '0', linestyle ='solid', linewidth=1.0, label= 'Uniaxial (Model)') +SEM_ax[i][j].plot(-1*multiaxial2_cmp_v3['displacement_x']/25.4, -1*multiaxial2_cmp_v3['s_xx'], c = '0', linestyle ='dashdot', linewidth=1.0, label= 'Biaxial (Model)') +SEM_ax[i][j].grid(color="grey", ls = '--', lw = 0.5) +# SEM_ax[i][j].xaxis.set_minor_locator(MultipleLocator(0.2)) +# SEM_ax[i][j].yaxis.set_minor_locator(MultipleLocator(0.2)) +SEM_ax[i][j].set_ylabel("Force [N]") +SEM_ax[i][j].set_xlabel("Deformation [mm]") +# SEM_ax[i][j].set_xlim(0, 0.01) +# SEM_ax[i][j].set_ylim(0,35) +SEM_ax[i][j].legend(loc='lower center', facecolor = 'white', framealpha = 0.4, edgecolor ='none') + + +# [i, j] = [0, 2] +# SEM_ax[i][j].set_title('Shear Cyclic Loading', fontsize=11) +# SEM_ax[i][j].plot(Lee_cyc_shr[0], Lee_cyc_shr[1], c = '0.6', linestyle ='solid', linewidth=1.0, label= 'Uniaxial (Lee)') +# SEM_ax[i][j].plot(shr_cyclic_v3['displacement_x']/25.4, shr_cyclic_v3['s_xy'], c = '0', linestyle ='solid', linewidth=1.0, label= 'Uniaxial_v3') +# SEM_ax[i][j].grid(color="grey", ls = '--', lw = 0.5) +# # SEM_ax[i][j].xaxis.set_minor_locator(MultipleLocator(0.2)) +# # SEM_ax[i][j].yaxis.set_minor_locator(MultipleLocator(0.2)) +# SEM_ax[i][j].set_xlabel("Strain [m/m]") +# SEM_ax[i][j].set_ylabel("Stress [MPa]") +# # SEM_ax[i][j].set_xlim(0, 0.001) +# # SEM_ax[i][j].set_ylim(0,0.15) +# # STT_image = mpimg.imread('./test/tests/plastic_damage_model_2017/analysis/STT.png') +# # place_image(STT_image, loc='lower right', ax=SEM_ax[i][j], pad=0, zoom=0.18) +# SEM_ax[i][j].legend(loc='lower right', facecolor = 'white', framealpha = 0.4, edgecolor ='none') + + + +[i, j] = [1, 0] +SEM_ax[i][j].set_title('Mesh sensitivity', fontsize=11) +SEM_ax[i][j].plot(Lee_uni_ten_msh_sen_2_elem[0], Lee_uni_ten_msh_sen_2_elem[1], c = '0.6', linestyle ='solid', linewidth=1.0, label= '1 element (Lee)') +SEM_ax[i][j].plot(Lee_uni_ten_msh_sen_4_elem[0], Lee_uni_ten_msh_sen_4_elem[1], c = '0.6', linestyle ='dashdot', linewidth=1.0, label= '2 elements (Lee)') +SEM_ax[i][j].plot(Lee_uni_ten_msh_sen_8_elem[0], Lee_uni_ten_msh_sen_8_elem[1], c = '0.6', linestyle ='dashed', linewidth=1.0, label= '4 elements (Lee)') +SEM_ax[i][j].plot(uni_ten_msh_sen_1ele_v3['displacement_x'], -1*uni_ten_msh_sen_1ele_v3['react_x'], c = '0', linestyle ='solid', linewidth=1.0, label= '1 element (Model)') +SEM_ax[i][j].plot(uni_ten_msh_sen_2ele_v3['displacement_x'], -1*uni_ten_msh_sen_2ele_v3['react_x'], c = '0', linestyle ='dashdot', linewidth=1.0, label= '2 elements (Model)') +SEM_ax[i][j].plot(uni_ten_msh_sen_4ele_v3['displacement_x'], -1*uni_ten_msh_sen_4ele_v3['react_x'], c = '0', linestyle ='dashed', linewidth=1.0, label= '4 elements (Model)') +SEM_ax[i][j].grid(color="grey", ls = '--', lw = 0.5) +# SEM_ax[i][j].xaxis.set_minor_locator(MultipleLocator(0.2)) +# SEM_ax[i][j].yaxis.set_minor_locator(MultipleLocator(0.01)) +SEM_ax[i][j].set_ylabel("Force [N]") +SEM_ax[i][j].set_xlabel("Deformation [mm]") +SEM_ax[i][j].set_xlim(0, 0.02) +SEM_ax[i][j].set_ylim(0,90) +SEM_ax[i][j].legend(loc='upper right', facecolor = 'white', framealpha = 0.4, edgecolor ='none') + + +[i, j] = [1, 1] +SEM_ax[i][j].set_title('Dilatancy check', fontsize=11) +SEM_ax[i][j].plot(Lee_dila_ap_02[0], -1*Lee_dila_ap_02[1], 'k', marker='o', ms=3, linewidth=1.0, label= '$\\alpha_p$ = 0.20 (Lee)') +SEM_ax[i][j].plot(Lee_dila_ap_025[0], -1*Lee_dila_ap_025[1], 'r', marker='o', ms=3, linewidth=1.0, label= '$\\alpha_p$ = 0.25 (Lee)') +SEM_ax[i][j].plot(Lee_dila_ap_03[0], -1*Lee_dila_ap_03[1], 'b', marker='o', ms=3, linewidth=1.0, label= '$\\alpha_p$ = 0.30 (Lee)') +SEM_ax[i][j].plot(vol_strain_ap_02_v3, -1*uni_cmp_dila_ap_02_v3['s_xx'], 'k', linewidth=1.0, label= '$\\alpha_p$ = 0.12 (Model)') +SEM_ax[i][j].plot(vol_strain_ap_025_v3, -1*uni_cmp_dila_ap_025_v3['s_xx'], 'r', linewidth=1.0, label= '$\\alpha_p$ = 0.15 (Model)') +SEM_ax[i][j].plot(vol_strain_ap_03_v3, -1*uni_cmp_dila_ap_03_v3['s_xx'], 'b', linewidth=1.0, label= '$\\alpha_p$ = 0.20 (Model)') +SEM_ax[i][j].grid(color="grey", ls = '--', lw = 0.5) +# SEM_ax[i][j].xaxis.set_minor_locator(MultipleLocator(0.2)) +# SEM_ax[i][j].yaxis.set_minor_locator(MultipleLocator(0.2)) +SEM_ax[i][j].set_ylabel("Force [N]") +SEM_ax[i][j].set_xlabel("Deformation [mm]") +SEM_ax[i][j].set_xlim(-0.001, 0.004) +# SEM_ax[i][j].set_ylim(0,0.15) +SEM_ax[i][j].legend(loc='upper right', facecolor = 'white', framealpha = 0.4, edgecolor ='none') + +plt.show() diff --git a/test/tests/plastic_damage_model/input/multiaxial1_cmp_v3.i b/test/tests/plastic_damage_model/input/multiaxial1_cmp_v3.i new file mode 100644 index 00000000..58379601 --- /dev/null +++ b/test/tests/plastic_damage_model/input/multiaxial1_cmp_v3.i @@ -0,0 +1,309 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 2 + ny = 2 + + xmin = -12.7 + xmax = 12.7 + + ymin = -12.7 + ymax = 12.7 + [] +[] + +# [Mesh] +# type = MeshGeneratorMesh +# [] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '-1E-3*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = ElementAverageValue + variable = stress_xx + [../] + # [./react_x] + # type = NodalSum + # variable = resid_x + # boundary = 'left' + # [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.1E4 + poissons_ratio = 0.18 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain] + type = ComputePlaneIncrementalStrain + [../] + [./pdm] + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.70 + damage_half_ten_str = 0.51 + frac_energy_ten = 12.3E-3 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.25 + frac_energy_cmp = 825E-3 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/multiaxial1_cmp_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/multiaxial1_ten_v3.i b/test/tests/plastic_damage_model/input/multiaxial1_ten_v3.i new file mode 100644 index 00000000..a97feb3e --- /dev/null +++ b/test/tests/plastic_damage_model/input/multiaxial1_ten_v3.i @@ -0,0 +1,309 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 2 + ny = 2 + + xmin = -12.7 + xmax = 12.7 + + ymin = -12.7 + ymax = 12.7 + [] +[] + +# [Mesh] +# type = MeshGeneratorMesh +# [] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '1E-4*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = ElementAverageValue + variable = stress_xx + [../] + # [./react_x] + # type = NodalSum + # variable = resid_x + # boundary = 'left' + # [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain] + type = ComputePlaneIncrementalStrain + [../] + [./pdm] + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.70 + damage_half_ten_str = 0.51 + frac_energy_ten = 12.3E-3 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 1750E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/multiaxial1_ten_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/multiaxial2_cmp_v3.i b/test/tests/plastic_damage_model/input/multiaxial2_cmp_v3.i new file mode 100644 index 00000000..7ea7d55e --- /dev/null +++ b/test/tests/plastic_damage_model/input/multiaxial2_cmp_v3.i @@ -0,0 +1,314 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 2 + ny = 2 + + xmin = -12.7 + xmax = 12.7 + + ymin = -12.7 + ymax = 12.7 + [] +[] + +# [Mesh] +# type = MeshGeneratorMesh +# [] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '-1E-3*t' + [../] + [./top_surface] + type = FunctionPresetBC + variable = disp_y + boundary = 'top' + function = '-1E-3*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = ElementAverageValue + variable = stress_xx + [../] + # [./react_x] + # type = NodalSum + # variable = resid_x + # boundary = 'left' + # [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.1E4 + poissons_ratio = 0.18 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain] + type = ComputePlaneIncrementalStrain + [../] + [./pdm] + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.70 + damage_half_ten_str = 0.51 + frac_energy_ten = 12.3E-3 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.25 + frac_energy_cmp = 825E-3 + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/multiaxial2_cmp_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/multiaxial2_ten_v3.i b/test/tests/plastic_damage_model/input/multiaxial2_ten_v3.i new file mode 100644 index 00000000..ab1ae4fd --- /dev/null +++ b/test/tests/plastic_damage_model/input/multiaxial2_ten_v3.i @@ -0,0 +1,315 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 2 + ny = 2 + + xmin = -12.7 + xmax = 12.7 + + ymin = -12.7 + ymax = 12.7 + [] +[] + +# [Mesh] +# type = MeshGeneratorMesh +# [] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '1E-4*t' + [../] + [./top_surface] + type = FunctionPresetBC + variable = disp_y + boundary = 'top' + function = '1E-4*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = ElementAverageValue + variable = stress_xx + [../] + # [./react_x] + # type = NodalSum + # variable = resid_x + # boundary = 'left' + # [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain] + type = ComputePlaneIncrementalStrain + [../] + [./pdm] + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.70 + damage_half_ten_str = 0.51 + frac_energy_ten = 12.3E-3 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 1750E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/multiaxial2_ten_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/shr_cyclic_v3.i b/test/tests/plastic_damage_model/input/shr_cyclic_v3.i new file mode 100644 index 00000000..fe937da3 --- /dev/null +++ b/test/tests/plastic_damage_model/input/shr_cyclic_v3.i @@ -0,0 +1,309 @@ +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 1 + ny = 1 + + xmin = -12.7 + xmax = 12.7 + + ymin = -12.7 + ymax = 12.7 + [] + construct_side_list_from_node_list = true +[] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./disp_x_zero] + type = FunctionPresetBC + variable = disp_x + # boundary = '1 4' + boundary = 'bottom' + function = '0.' + [../] + [./disp_y_zero] + type = FunctionPresetBC + variable = disp_y + boundary = 'top bottom' + # boundary = '1 2 3 4' + function = '0.' + [../] + [./disp_x_loading] + type = FunctionPresetBC + variable = disp_x + boundary = 'top' + # boundary = '2 3' + function = '1E-3*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'top' + # boundary = '2 3' + [../] + [./eps_xy] + type = ElementAverageValue + variable = strain_xy + [../] + [./s_xy] + type = ElementAverageValue + variable = stress_xy + [../] + # [./react_x] + # type = NodalSum + # variable = resid_x + # boundary = 'left' + # [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.1E4 + poissons_ratio = 0.18 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain] + type = ComputePlaneIncrementalStrain + [../] + [./pdm] + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.70 + damage_half_ten_str = 0.51 + frac_energy_ten = 12.3E-3 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.25 + frac_energy_cmp = 825E-3 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 0.2 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/shr_cyclic_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_cmp_dila_ap25_v3.i b/test/tests/plastic_damage_model/input/uni_cmp_dila_ap25_v3.i new file mode 100644 index 00000000..6f529c97 --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_cmp_dila_ap25_v3.i @@ -0,0 +1,499 @@ +[Mesh] + type = GeneratedMesh + dim = 3 + nx = 1 + ny = 1 + nz = 1 + xmin = -0.5 + xmax = 0.5 + ymin = -0.5 + ymax = 0.5 + zmin = -0.5 + zmax = 0.5 +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./disp_z] + [../] +[] + +[Kernels] + [./TensorMechanics] + displacements = 'disp_x disp_y disp_z' + save_in = 'resid_x resid_y resid_z' + [../] +[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./resid_z] + [../] + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./min_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./mid_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./max_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma0] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma1] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma2] + order = CONSTANT + family = MONOMIAL + [../] + [./f0] + order = CONSTANT + family = MONOMIAL + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_xz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xz + index_i = 0 + index_j = 2 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_yz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yz + index_i = 1 + index_j = 2 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + [./strain_xx] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_xz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xz + index_i = 0 + index_j = 2 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_yz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yz + index_i = 1 + index_j = 2 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_zz + index_i = 2 + index_j = 2 + [../] + [./f0_auxk] + type = MaterialStdVectorAux + property = plastic_yield_function + index = 0 + variable = f0 + [../] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./min_ep] + type = MaterialRealAux + property = min_ep + variable = min_ep + [../] + [./mid_ep] + type = MaterialRealAux + property = mid_ep + variable = mid_ep + [../] + [./max_ep] + type = MaterialRealAux + property = max_ep + variable = max_ep + [../] + [./sigma0] + type = MaterialRealAux + property = damaged_min_principal_stress + variable = sigma0 + [../] + [./sigma1] + type = MaterialRealAux + property = damaged_mid_principal_stress + variable = sigma1 + [../] + [./sigma2] + type = MaterialRealAux + property = damaged_max_principal_stress + variable = sigma2 + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] +[] + +[Postprocessors] + # active ='react_x react_y react_z displacement_x' + [./react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [../] + [./react_y] + type = NodalSum + variable = resid_y + boundary = 'left' + [../] + [./react_z] + type = NodalSum + variable = resid_z + boundary = 'left' + [../] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = PointValue + point = '0 0 0' + variable = stress_xx + [../] + [./s_xy] + type = PointValue + point = '0 0 0' + variable = stress_xy + [../] + [./s_xz] + type = PointValue + point = '0 0 0' + variable = stress_xz + [../] + [./s_yy] + type = PointValue + point = '0 0 0' + variable = stress_yy + [../] + [./s_yz] + type = PointValue + point = '0 0 0' + variable = stress_yz + [../] + [./s_zz] + type = PointValue + point = '0 0 0' + variable = stress_zz + [../] + [./e_xx] + type = PointValue + point = '0 0 0' + variable = strain_xx + [../] + [./e_xy] + type = PointValue + point = '0 0 0' + variable = strain_xy + [../] + [./e_xz] + type = PointValue + point = '0 0 0' + variable = strain_xz + [../] + [./e_yy] + type = PointValue + point = '0 0 0' + variable = strain_yy + [../] + [./e_yz] + type = PointValue + point = '0 0 0' + variable = strain_yz + [../] + [./e_zz] + type = PointValue + point = '0 0 0' + variable = strain_zz + [../] + [./max_ep] + type = PointValue + point = '0 0 0' + variable = max_ep + [../] + [./mid_ep] + type = PointValue + point = '0 0 0' + variable = mid_ep + [../] + [./min_ep] + type = PointValue + point = '0 0 0' + variable = min_ep + [../] + [./sigma0] + type = PointValue + point = '0 0 0' + variable = sigma0 + [../] + [./sigma1] + type = PointValue + point = '0 0 0' + variable = sigma1 + [../] + [./sigma2] + type = PointValue + point = '0 0 0' + variable = sigma2 + [../] + [./f0] + type = PointValue + point = '0 0 0' + variable = f0 + [../] + [./D] + type = PointValue + point = '0 0 0' + variable = D + [../] + [./intnl0] + type = PointValue + point = '0 0 0' + variable = intnl0 + [../] + [./intnl1] + type = PointValue + point = '0 0 0' + variable = intnl1 + [../] +[] + +[BCs] + [./x_r] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '-2E-5*x*t' + [../] + [./x_l] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0' + [../] + [./y_l] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0' + [../] + [./z_l] + type = FunctionPresetBC + variable = disp_z + boundary = 'back' + function = '0' + [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./strain] + type = ComputeIncrementalSmallStrain + displacements = 'disp_x disp_y disp_z' + [../] + [./pdm] + # type = PlasticDamageStressUpdateOne + # yield_function_tol = 1.E-5 + # alpha = 0.109 # fb0 = -19.47 + # alpha_p = 0.23 + # s0 = 0.001 + # + # at = 0.0001 # parameter to adjust yield strength + # Dt = 0.51 # degradation at half of maximum tensile strength + # ft = 3.48 # maximum yield strength + # gt = 4.8E-4 + # + # ac = 3.77 # parameter to adjust yield strength + # Dc = 0.40 # degradation at maximum compressive strength + # fc = 27.6 # maximum compressive strength of concrete + # gc = 6.889E-2 # + # + # tip_smoother = 1.E-6 + # smoothing_tol = 1.E-3 + + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.155 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.88 + damage_half_ten_str = 0.51 + frac_energy_ten = 4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + [../] +[] + +[Executioner] +# nl_abs_tol=1E-6 +# petsc_options_iname = '-pc_type' +# petsc_options_value = 'lu' + + end_time = 4000 + dt = 5 + type = Transient +[] + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_cmp_dila_ap_025_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_cmp_dila_ap2_v3.i b/test/tests/plastic_damage_model/input/uni_cmp_dila_ap2_v3.i new file mode 100644 index 00000000..82a06133 --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_cmp_dila_ap2_v3.i @@ -0,0 +1,499 @@ +[Mesh] + type = GeneratedMesh + dim = 3 + nx = 1 + ny = 1 + nz = 1 + xmin = -0.5 + xmax = 0.5 + ymin = -0.5 + ymax = 0.5 + zmin = -0.5 + zmax = 0.5 +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./disp_z] + [../] +[] + +[Kernels] + [./TensorMechanics] + displacements = 'disp_x disp_y disp_z' + save_in = 'resid_x resid_y resid_z' + [../] +[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./resid_z] + [../] + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./min_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./mid_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./max_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma0] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma1] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma2] + order = CONSTANT + family = MONOMIAL + [../] + [./f0] + order = CONSTANT + family = MONOMIAL + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_xz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xz + index_i = 0 + index_j = 2 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_yz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yz + index_i = 1 + index_j = 2 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + [./strain_xx] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_xz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xz + index_i = 0 + index_j = 2 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_yz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yz + index_i = 1 + index_j = 2 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_zz + index_i = 2 + index_j = 2 + [../] + [./f0_auxk] + type = MaterialStdVectorAux + property = plastic_yield_function + index = 0 + variable = f0 + [../] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./min_ep] + type = MaterialRealAux + property = min_ep + variable = min_ep + [../] + [./mid_ep] + type = MaterialRealAux + property = mid_ep + variable = mid_ep + [../] + [./max_ep] + type = MaterialRealAux + property = max_ep + variable = max_ep + [../] + [./sigma0] + type = MaterialRealAux + property = damaged_min_principal_stress + variable = sigma0 + [../] + [./sigma1] + type = MaterialRealAux + property = damaged_mid_principal_stress + variable = sigma1 + [../] + [./sigma2] + type = MaterialRealAux + property = damaged_max_principal_stress + variable = sigma2 + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] +[] + +[Postprocessors] + # active ='react_x react_y react_z displacement_x' + [./react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [../] + [./react_y] + type = NodalSum + variable = resid_y + boundary = 'left' + [../] + [./react_z] + type = NodalSum + variable = resid_z + boundary = 'left' + [../] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = PointValue + point = '0 0 0' + variable = stress_xx + [../] + [./s_xy] + type = PointValue + point = '0 0 0' + variable = stress_xy + [../] + [./s_xz] + type = PointValue + point = '0 0 0' + variable = stress_xz + [../] + [./s_yy] + type = PointValue + point = '0 0 0' + variable = stress_yy + [../] + [./s_yz] + type = PointValue + point = '0 0 0' + variable = stress_yz + [../] + [./s_zz] + type = PointValue + point = '0 0 0' + variable = stress_zz + [../] + [./e_xx] + type = PointValue + point = '0 0 0' + variable = strain_xx + [../] + [./e_xy] + type = PointValue + point = '0 0 0' + variable = strain_xy + [../] + [./e_xz] + type = PointValue + point = '0 0 0' + variable = strain_xz + [../] + [./e_yy] + type = PointValue + point = '0 0 0' + variable = strain_yy + [../] + [./e_yz] + type = PointValue + point = '0 0 0' + variable = strain_yz + [../] + [./e_zz] + type = PointValue + point = '0 0 0' + variable = strain_zz + [../] + [./max_ep] + type = PointValue + point = '0 0 0' + variable = max_ep + [../] + [./mid_ep] + type = PointValue + point = '0 0 0' + variable = mid_ep + [../] + [./min_ep] + type = PointValue + point = '0 0 0' + variable = min_ep + [../] + [./sigma0] + type = PointValue + point = '0 0 0' + variable = sigma0 + [../] + [./sigma1] + type = PointValue + point = '0 0 0' + variable = sigma1 + [../] + [./sigma2] + type = PointValue + point = '0 0 0' + variable = sigma2 + [../] + [./f0] + type = PointValue + point = '0 0 0' + variable = f0 + [../] + [./D] + type = PointValue + point = '0 0 0' + variable = D + [../] + [./intnl0] + type = PointValue + point = '0 0 0' + variable = intnl0 + [../] + [./intnl1] + type = PointValue + point = '0 0 0' + variable = intnl1 + [../] +[] + +[BCs] + [./x_r] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '-2E-5*x*t' + [../] + [./x_l] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0' + [../] + [./y_l] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0' + [../] + [./z_l] + type = FunctionPresetBC + variable = disp_z + boundary = 'back' + function = '0' + [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./strain] + type = ComputeIncrementalSmallStrain + displacements = 'disp_x disp_y disp_z' + [../] + [./pdm] + # type = PlasticDamageStressUpdateOne + # yield_function_tol = 1.E-5 + # alpha = 0.109 # fb0 = -19.47 + # alpha_p = 0.23 + # s0 = 0.001 + # + # at = 0.0001 # parameter to adjust yield strength + # Dt = 0.51 # degradation at half of maximum tensile strength + # ft = 3.48 # maximum yield strength + # gt = 4.8E-4 + # + # ac = 3.77 # parameter to adjust yield strength + # Dc = 0.40 # degradation at maximum compressive strength + # fc = 27.6 # maximum compressive strength of concrete + # gc = 6.889E-2 # + # + # tip_smoother = 1.E-6 + # smoothing_tol = 1.E-3 + + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.12 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.88 + damage_half_ten_str = 0.51 + frac_energy_ten = 4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + [../] +[] + +[Executioner] +# nl_abs_tol=1E-6 +# petsc_options_iname = '-pc_type' +# petsc_options_value = 'lu' + + end_time = 4000 + dt = 5 + type = Transient +[] + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_cmp_dila_ap_02_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_cmp_dila_ap3_v3.i b/test/tests/plastic_damage_model/input/uni_cmp_dila_ap3_v3.i new file mode 100644 index 00000000..e2ce7806 --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_cmp_dila_ap3_v3.i @@ -0,0 +1,499 @@ +[Mesh] + type = GeneratedMesh + dim = 3 + nx = 1 + ny = 1 + nz = 1 + xmin = -0.5 + xmax = 0.5 + ymin = -0.5 + ymax = 0.5 + zmin = -0.5 + zmax = 0.5 +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./disp_z] + [../] +[] + +[Kernels] + [./TensorMechanics] + displacements = 'disp_x disp_y disp_z' + save_in = 'resid_x resid_y resid_z' + [../] +[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./resid_z] + [../] + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./min_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./mid_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./max_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma0] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma1] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma2] + order = CONSTANT + family = MONOMIAL + [../] + [./f0] + order = CONSTANT + family = MONOMIAL + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_xz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xz + index_i = 0 + index_j = 2 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_yz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yz + index_i = 1 + index_j = 2 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + [./strain_xx] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_xz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xz + index_i = 0 + index_j = 2 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_yz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yz + index_i = 1 + index_j = 2 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_zz + index_i = 2 + index_j = 2 + [../] + [./f0_auxk] + type = MaterialStdVectorAux + property = plastic_yield_function + index = 0 + variable = f0 + [../] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./min_ep] + type = MaterialRealAux + property = min_ep + variable = min_ep + [../] + [./mid_ep] + type = MaterialRealAux + property = mid_ep + variable = mid_ep + [../] + [./max_ep] + type = MaterialRealAux + property = max_ep + variable = max_ep + [../] + [./sigma0] + type = MaterialRealAux + property = damaged_min_principal_stress + variable = sigma0 + [../] + [./sigma1] + type = MaterialRealAux + property = damaged_mid_principal_stress + variable = sigma1 + [../] + [./sigma2] + type = MaterialRealAux + property = damaged_max_principal_stress + variable = sigma2 + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] +[] + +[Postprocessors] + # active ='react_x react_y react_z displacement_x' + [./react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [../] + [./react_y] + type = NodalSum + variable = resid_y + boundary = 'left' + [../] + [./react_z] + type = NodalSum + variable = resid_z + boundary = 'left' + [../] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = PointValue + point = '0 0 0' + variable = stress_xx + [../] + [./s_xy] + type = PointValue + point = '0 0 0' + variable = stress_xy + [../] + [./s_xz] + type = PointValue + point = '0 0 0' + variable = stress_xz + [../] + [./s_yy] + type = PointValue + point = '0 0 0' + variable = stress_yy + [../] + [./s_yz] + type = PointValue + point = '0 0 0' + variable = stress_yz + [../] + [./s_zz] + type = PointValue + point = '0 0 0' + variable = stress_zz + [../] + [./e_xx] + type = PointValue + point = '0 0 0' + variable = strain_xx + [../] + [./e_xy] + type = PointValue + point = '0 0 0' + variable = strain_xy + [../] + [./e_xz] + type = PointValue + point = '0 0 0' + variable = strain_xz + [../] + [./e_yy] + type = PointValue + point = '0 0 0' + variable = strain_yy + [../] + [./e_yz] + type = PointValue + point = '0 0 0' + variable = strain_yz + [../] + [./e_zz] + type = PointValue + point = '0 0 0' + variable = strain_zz + [../] + [./max_ep] + type = PointValue + point = '0 0 0' + variable = max_ep + [../] + [./mid_ep] + type = PointValue + point = '0 0 0' + variable = mid_ep + [../] + [./min_ep] + type = PointValue + point = '0 0 0' + variable = min_ep + [../] + [./sigma0] + type = PointValue + point = '0 0 0' + variable = sigma0 + [../] + [./sigma1] + type = PointValue + point = '0 0 0' + variable = sigma1 + [../] + [./sigma2] + type = PointValue + point = '0 0 0' + variable = sigma2 + [../] + [./f0] + type = PointValue + point = '0 0 0' + variable = f0 + [../] + [./D] + type = PointValue + point = '0 0 0' + variable = D + [../] + [./intnl0] + type = PointValue + point = '0 0 0' + variable = intnl0 + [../] + [./intnl1] + type = PointValue + point = '0 0 0' + variable = intnl1 + [../] +[] + +[BCs] + [./x_r] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '-2E-5*x*t' + [../] + [./x_l] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0' + [../] + [./y_l] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0' + [../] + [./z_l] + type = FunctionPresetBC + variable = disp_z + boundary = 'back' + function = '0' + [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./strain] + type = ComputeIncrementalSmallStrain + displacements = 'disp_x disp_y disp_z' + [../] + [./pdm] + # type = PlasticDamageStressUpdateOne + # yield_function_tol = 1.E-5 + # alpha = 0.109 # fb0 = -19.47 + # alpha_p = 0.23 + # s0 = 0.001 + # + # at = 0.0001 # parameter to adjust yield strength + # Dt = 0.51 # degradation at half of maximum tensile strength + # ft = 3.48 # maximum yield strength + # gt = 4.8E-4 + # + # ac = 3.77 # parameter to adjust yield strength + # Dc = 0.40 # degradation at maximum compressive strength + # fc = 27.6 # maximum compressive strength of concrete + # gc = 6.889E-2 # + # + # tip_smoother = 1.E-6 + # smoothing_tol = 1.E-3 + + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.195 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.88 + damage_half_ten_str = 0.51 + frac_energy_ten = 4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + [../] +[] + +[Executioner] +# nl_abs_tol=1E-6 +# petsc_options_iname = '-pc_type' +# petsc_options_value = 'lu' + + end_time = 4000 + dt = 5 + type = Transient +[] + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_cmp_dila_ap_03_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_cmp_v3.i b/test/tests/plastic_damage_model/input/uni_cmp_v3.i new file mode 100644 index 00000000..1d71fba4 --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_cmp_v3.i @@ -0,0 +1,499 @@ +[Mesh] + type = GeneratedMesh + dim = 3 + nx = 1 + ny = 1 + nz = 1 + xmin = -0.5 + xmax = 0.5 + ymin = -0.5 + ymax = 0.5 + zmin = -0.5 + zmax = 0.5 +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./disp_z] + [../] +[] + +[Kernels] + [./TensorMechanics] + displacements = 'disp_x disp_y disp_z' + save_in = 'resid_x resid_y resid_z' + [../] +[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./resid_z] + [../] + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yz] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] + [./min_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./mid_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./max_ep] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma0] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma1] + order = CONSTANT + family = MONOMIAL + [../] + [./sigma2] + order = CONSTANT + family = MONOMIAL + [../] + [./f0] + order = CONSTANT + family = MONOMIAL + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_xz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xz + index_i = 0 + index_j = 2 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_yz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yz + index_i = 1 + index_j = 2 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + [./strain_xx] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_xz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_xz + index_i = 0 + index_j = 2 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_yz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_yz + index_i = 1 + index_j = 2 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = elastic_strain + variable = strain_zz + index_i = 2 + index_j = 2 + [../] + [./f0_auxk] + type = MaterialStdVectorAux + property = plastic_yield_function + index = 0 + variable = f0 + [../] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./min_ep] + type = MaterialRealAux + property = min_ep + variable = min_ep + [../] + [./mid_ep] + type = MaterialRealAux + property = mid_ep + variable = mid_ep + [../] + [./max_ep] + type = MaterialRealAux + property = max_ep + variable = max_ep + [../] + [./sigma0] + type = MaterialRealAux + property = damaged_min_principal_stress + variable = sigma0 + [../] + [./sigma1] + type = MaterialRealAux + property = damaged_mid_principal_stress + variable = sigma1 + [../] + [./sigma2] + type = MaterialRealAux + property = damaged_max_principal_stress + variable = sigma2 + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] +[] + +[Postprocessors] + active ='react_x react_y react_z displacement_x' + [./react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [../] + [./react_y] + type = NodalSum + variable = resid_y + boundary = 'left' + [../] + [./react_z] + type = NodalSum + variable = resid_z + boundary = 'left' + [../] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = PointValue + point = '0 0 0' + variable = stress_xx + [../] + [./s_xy] + type = PointValue + point = '0 0 0' + variable = stress_xy + [../] + [./s_xz] + type = PointValue + point = '0 0 0' + variable = stress_xz + [../] + [./s_yy] + type = PointValue + point = '0 0 0' + variable = stress_yy + [../] + [./s_yz] + type = PointValue + point = '0 0 0' + variable = stress_yz + [../] + [./s_zz] + type = PointValue + point = '0 0 0' + variable = stress_zz + [../] + [./e_xx] + type = PointValue + point = '0 0 0' + variable = strain_xx + [../] + [./e_xy] + type = PointValue + point = '0 0 0' + variable = strain_xy + [../] + [./e_xz] + type = PointValue + point = '0 0 0' + variable = strain_xz + [../] + [./e_yy] + type = PointValue + point = '0 0 0' + variable = strain_yy + [../] + [./e_yz] + type = PointValue + point = '0 0 0' + variable = strain_yz + [../] + [./e_zz] + type = PointValue + point = '0 0 0' + variable = strain_zz + [../] + [./max_ep] + type = PointValue + point = '0 0 0' + variable = max_ep + [../] + [./mid_ep] + type = PointValue + point = '0 0 0' + variable = mid_ep + [../] + [./min_ep] + type = PointValue + point = '0 0 0' + variable = min_ep + [../] + [./sigma0] + type = PointValue + point = '0 0 0' + variable = sigma0 + [../] + [./sigma1] + type = PointValue + point = '0 0 0' + variable = sigma1 + [../] + [./sigma2] + type = PointValue + point = '0 0 0' + variable = sigma2 + [../] + [./f0] + type = PointValue + point = '0 0 0' + variable = f0 + [../] + [./D] + type = PointValue + point = '0 0 0' + variable = D + [../] + [./intnl0] + type = PointValue + point = '0 0 0' + variable = intnl0 + [../] + [./intnl1] + type = PointValue + point = '0 0 0' + variable = intnl1 + [../] +[] + +[BCs] + [./x_r] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '-2E-5*x*t' + [../] + [./x_l] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0' + [../] + [./y_l] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom' + function = '0' + [../] + [./z_l] + type = FunctionPresetBC + variable = disp_z + boundary = 'back' + function = '0' + [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./strain] + type = ComputeIncrementalSmallStrain + displacements = 'disp_x disp_y disp_z' + [../] + [./pdm] + # type = PlasticDamageStressUpdateOne + # yield_function_tol = 1.E-5 + # alpha = 0.109 # fb0 = -19.47 + # alpha_p = 0.23 + # s0 = 0.001 + # + # at = 0.0001 # parameter to adjust yield strength + # Dt = 0.51 # degradation at half of maximum tensile strength + # ft = 3.48 # maximum yield strength + # gt = 4.8E-4 + # + # ac = 3.77 # parameter to adjust yield strength + # Dc = 0.40 # degradation at maximum compressive strength + # fc = 27.6 # maximum compressive strength of concrete + # gc = 6.889E-2 # + # + # tip_smoother = 1.E-6 + # smoothing_tol = 1.E-3 + + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.12 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.88 + damage_half_ten_str = 0.51 + frac_energy_ten = 4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.25 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + [../] +[] + +[Executioner] +# nl_abs_tol=1E-6 +# petsc_options_iname = '-pc_type' +# petsc_options_value = 'lu' + + end_time = 4000 + dt = 5 + type = Transient +[] + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_cmp_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_ten_msh_sen_1ele_v3.i b/test/tests/plastic_damage_model/input/uni_ten_msh_sen_1ele_v3.i new file mode 100644 index 00000000..c80d73f9 --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_ten_msh_sen_1ele_v3.i @@ -0,0 +1,317 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 1 + ny = 1 + + xmin = -25.4 + xmax = 25.4 + + ymin = -12.7 + ymax = 12.7 + [] + [./extra_nodeset1] + type = ExtraNodesetGenerator + input = gmg + new_boundary = 'bottom_left' + coord = '-25.4 -12.7' + [] + [./extra_nodeset2] + type = ExtraNodesetGenerator + input = extra_nodeset1 + new_boundary = 'top_left' + coord = '-25.4 12.7' + [] +[] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom_left' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '1E-4*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [../] +[] + +[Materials] + [./elasticity_tensor_right] + type = ComputeIsotropicElasticityTensor + block = '0' + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress_right] + type = ComputeMultipleInelasticDamageStress + block = '0' + inelastic_models = pdm_right + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain_right] + type = ComputePlaneIncrementalStrain + block = '0' + [../] + [./pdm_right] + type = PlasticDamageStressUpdateOne_v3 + block = '0' + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.40 + ft_ep_slope_factor_at_zero_ep = 0.90 + damage_half_ten_str = 0.27 + frac_energy_ten = 25E-3 #4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_ten_msh_sen_1ele_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_ten_msh_sen_2ele_v3.i b/test/tests/plastic_damage_model/input/uni_ten_msh_sen_2ele_v3.i new file mode 100644 index 00000000..9bec2b2b --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_ten_msh_sen_2ele_v3.i @@ -0,0 +1,367 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 2 + ny = 1 + + xmin = -25.4 + xmax = 25.4 + + ymin = -12.7 + ymax = 12.7 + [] + [./subdomains] + type = SubdomainBoundingBoxGenerator + input = gmg + bottom_left = '0 12.7 0' + top_right = '25.4 -12.7 0' + block_id = '1' + location = INSIDE + [] + [./extra_nodeset1] + type = ExtraNodesetGenerator + input = subdomains + new_boundary = 'bottom_left' + coord = '-25.4 -12.7' + [] + [./extra_nodeset2] + type = ExtraNodesetGenerator + input = extra_nodeset1 + new_boundary = 'top_left' + coord = '-25.4 12.7' + [] +[] + +# [Mesh] +# type = MeshGeneratorMesh +# [] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom_left' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '1E-4*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [../] +[] + +[Materials] + [./elasticity_tensor_left] + type = ComputeIsotropicElasticityTensor + block = '0' + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress_left] + type = ComputeMultipleInelasticDamageStress + block = '0' + inelastic_models = pdm_left + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain_left] + type = ComputePlaneIncrementalStrain + block = '0' + [../] + [./pdm_left] + type = PlasticDamageStressUpdateOne_v3 + block = '0' + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.90 + damage_half_ten_str = 0.27 + frac_energy_ten = 25E-3 #4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] + [./elasticity_tensor_right] + type = ComputeIsotropicElasticityTensor + block = '1' + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress_right] + type = ComputeMultipleInelasticDamageStress + block = '1' + inelastic_models = pdm_right + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain_right] + type = ComputePlaneIncrementalStrain + block = '1' + [../] + [./pdm_right] + type = PlasticDamageStressUpdateOne_v3 + block = '1' + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.40 + ft_ep_slope_factor_at_zero_ep = 0.90 + damage_half_ten_str = 0.27 + frac_energy_ten = 25E-3 #4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_ten_msh_sen_2ele_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_ten_msh_sen_4ele_v3.i b/test/tests/plastic_damage_model/input/uni_ten_msh_sen_4ele_v3.i new file mode 100644 index 00000000..4481ebd5 --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_ten_msh_sen_4ele_v3.i @@ -0,0 +1,367 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 4 + ny = 1 + + xmin = -25.4 + xmax = 25.4 + + ymin = -12.7 + ymax = 12.7 + [] + [./subdomains] + type = SubdomainBoundingBoxGenerator + input = gmg + bottom_left = '12.7 12.7 0' + top_right = '25.4 -12.7 0' + block_id = '1' + location = INSIDE + [] + [./extra_nodeset1] + type = ExtraNodesetGenerator + input = subdomains + new_boundary = 'bottom_left' + coord = '-25.4 -12.7' + [] + [./extra_nodeset2] + type = ExtraNodesetGenerator + input = extra_nodeset1 + new_boundary = 'top_left' + coord = '-25.4 12.7' + [] +[] + +# [Mesh] +# type = MeshGeneratorMesh +# [] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom_left' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '1E-4*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./react_x] + type = NodalSum + variable = resid_x + boundary = 'left' + [../] +[] + +[Materials] + [./elasticity_tensor_left] + type = ComputeIsotropicElasticityTensor + block = '0' + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress_left] + type = ComputeMultipleInelasticDamageStress + block = '0' + inelastic_models = pdm_left + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain_left] + type = ComputePlaneIncrementalStrain + block = '0' + [../] + [./pdm_left] + type = PlasticDamageStressUpdateOne_v3 + block = '0' + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.90 + damage_half_ten_str = 0.27 + frac_energy_ten = 25E-3 #4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] + [./elasticity_tensor_right] + type = ComputeIsotropicElasticityTensor + block = '1' + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress_right] + type = ComputeMultipleInelasticDamageStress + block = '1' + inelastic_models = pdm_right + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain_right] + type = ComputePlaneIncrementalStrain + block = '1' + [../] + [./pdm_right] + type = PlasticDamageStressUpdateOne_v3 + block = '1' + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.40 + ft_ep_slope_factor_at_zero_ep = 0.90 + damage_half_ten_str = 0.27 + frac_energy_ten = 25E-3 #4.8E-4 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 6.889E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_ten_msh_sen_4ele_v3 + exodus = true + [./csv] + type = CSV + [../] +[] diff --git a/test/tests/plastic_damage_model/input/uni_ten_v3.i b/test/tests/plastic_damage_model/input/uni_ten_v3.i new file mode 100644 index 00000000..76c18799 --- /dev/null +++ b/test/tests/plastic_damage_model/input/uni_ten_v3.i @@ -0,0 +1,329 @@ +# [MeshGenerators] +[Mesh] + [gmg] + type = GeneratedMeshGenerator + dim = 2 + nx = 1 + ny = 1 + + xmin = -12.7 + xmax = 12.7 + + ymin = -12.7 + ymax = 12.7 + [] + [./subdomains] + type = SubdomainBoundingBoxGenerator + input = gmg + bottom_left = '0 12.7 0' + top_right = '12.7 -12.7 0' + block_id = '1' + location = INSIDE + [] + [./extra_nodeset1] + type = ExtraNodesetGenerator + input = subdomains + new_boundary = 'bottom_left' + coord = '-12.7 -12.7' + [] + [./extra_nodeset2] + type = ExtraNodesetGenerator + input = extra_nodeset1 + new_boundary = 'top_left' + coord = '-12.7 12.7' + [] +[] + +# [Mesh] +# type = MeshGeneratorMesh +# [] + +[GlobalParams] + displacements = 'disp_x disp_y' + volumetric_locking_correction = true + out_of_plane_strain = strain_zz +[] + +[Variables] + [./disp_x] + [../] + [./disp_y] + [../] + [./strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[Kernels] + [./TensorMechanics] + use_displaced_mesh = true + save_in = 'resid_x resid_y' + [../] + + [./solid_z] + type = WeakPlaneStress + variable = strain_zz + use_displaced_mesh = true + [../] +[] + +#[Modules/TensorMechanics/Master] +# [./all] +# add_variables = true +# incremental = false #true +# strain = small +# generate_output = 'stress_xx stress_xy stress_xz stress_yy stress_yz stress_zz +# strain_xx strain_xy strain_xz strain_yy strain_yz strain_zz +# max_principal_stress mid_principal_stress min_principal_stress +# secondinv_stress thirdinv_stress vonmises_stress +# secondinv_strain thirdinv_strain +# elastic_strain_xx elastic_strain_xy elastic_strain_xz elastic_strain_yy elastic_strain_yz elastic_strain_zz' +## plastic_strain_xx plastic_strain_xy plastic_strain_xz plastic_strain_yy plastic_strain_yz plastic_strain_zz' +# save_in = 'resid_x resid_y' +# planar_formulation = plane_strain +# [../] +#[] + +[AuxVariables] + [./resid_x] + [../] + [./resid_y] + [../] + [./D] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl0] + order = CONSTANT + family = MONOMIAL + [../] + [./intnl1] + order = CONSTANT + family = MONOMIAL + [../] + + [./stress_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./stress_zz] + order = CONSTANT + family = MONOMIAL + [../] + + [./strain_xx] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_xy] + order = CONSTANT + family = MONOMIAL + [../] + [./strain_yy] + order = CONSTANT + family = MONOMIAL + [../] + [./material_strain_zz] + order = CONSTANT + family = MONOMIAL + [../] +[] + +[AuxKernels] + [./D_auxk] + type = MaterialRealAux + property = Damage_Variable + variable = D + [../] + [./intnl0_auxk] + type = MaterialRealAux + property = damage_parameter_for_tension + variable = intnl0 + [../] + [./intnl1_auxk] + type = MaterialRealAux + property = damage_parameter_for_compression + variable = intnl1 + [../] + + [./stress_xx] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xx + index_i = 0 + index_j = 0 + [../] + [./stress_xy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_xy + index_i = 0 + index_j = 1 + [../] + [./stress_yy] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_yy + index_i = 1 + index_j = 1 + [../] + [./stress_zz] + type = RankTwoAux + rank_two_tensor = stress + variable = stress_zz + index_i = 2 + index_j = 2 + [../] + + [./strain_xx] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xx + index_i = 0 + index_j = 0 + [../] + [./strain_xy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_xy + index_i = 0 + index_j = 1 + [../] + [./strain_yy] + type = RankTwoAux + rank_two_tensor = total_strain + variable = strain_yy + index_i = 1 + index_j = 1 + [../] + [./strain_zz] + type = RankTwoAux + rank_two_tensor = total_strain + variable = material_strain_zz + index_i = 2 + index_j = 2 + [../] +[] + +[BCs] + [./left_x] + type = FunctionPresetBC + variable = disp_x + boundary = 'left' + function = '0.' + [../] + [./y1] + type = FunctionPresetBC + variable = disp_y + boundary = 'bottom_left' + function = '0.' + [../] + [./right_surface] + type = FunctionPresetBC + variable = disp_x + boundary = 'right' + function = '1E-4*t' + [../] +[] + +[Postprocessors] + [./displacement_x] + type = AverageNodalVariableValue + variable = disp_x + boundary = 'right' + [../] + [./s_xx] + type = ElementAverageValue + variable = stress_xx + [../] + # [./react_x] + # type = NodalSum + # variable = resid_x + # boundary = 'left' + # [../] +[] + +[Materials] + [./elasticity_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 3.17E4 + poissons_ratio = 0.18 + [../] + [./stress] + type = ComputeMultipleInelasticDamageStress + inelastic_models = pdm + perform_finite_strain_rotations = false + tangent_operator = nonlinear + [../] + [./strain] + type = ComputePlaneIncrementalStrain + [../] + [./pdm] + type = PlasticDamageStressUpdateOne_v3 + factor_relating_biaxial_unixial_cmp_str = 0.109 # fb0 = -20.862 + factor_controlling_dilatancy = 0.23 + stiff_recovery_factor = 0.001 + + yield_ten_str = 3.48 + ft_ep_slope_factor_at_zero_ep = 0.70 + damage_half_ten_str = 0.51 + frac_energy_ten = 12.3E-3 + + yield_cmp_str = 18.30 + max_cmp_str = 27.60 + damage_max_cmp_str = 0.40 + frac_energy_cmp = 1750E-2 + + yield_function_tol = 1.E-5 + tip_smoother = 1.E-6 + smoothing_tol = 1.E-3 + [../] +[] + +[Preconditioning] + active = SMP + [./SMP] + type = SMP + full = true + [../] + [./FDP] + type = FDP + full = true + [../] +[] + +[Executioner] + solve_type = 'NEWTON' + nl_max_its = 100 + nl_abs_tol = 1.E-5 + nl_rel_tol = 1E-3 + + line_search = none + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + + petsc_options = '-snes_converged_reason' + + type = Transient + end_time = 200 + dt = 1 +[] + + +[Outputs] + file_base = ./test/tests/plastic_damage_model/output/uni_ten_v3 + exodus = true + [./csv] + type = CSV + [../] +[]