Adding plastic_damage_model for concrete (idaholab#184)

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

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

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>();
+
+/**
+ * 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<Real> & _D;
+  const MaterialProperty<Real> & _D_old;
+  const MaterialProperty<Real> & _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

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>();
+
+/**
+ * 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<Real> & _max_principal;
+   MaterialProperty<Real> & _min_principal;
+   MaterialProperty<Real> & _intnl0;
+   MaterialProperty<Real> & _intnl1;
+   MaterialProperty<Real> & _ele_len;
+   MaterialProperty<Real> & _gt;
+   MaterialProperty<Real> & _gc;
+
+   MaterialProperty<Real> & _tD;
+   MaterialProperty<Real> & _cD;
+   MaterialProperty<Real> & _D;
+   MaterialProperty<Real> & _min_ep;
+   MaterialProperty<Real> & _mid_ep;
+   MaterialProperty<Real> & _max_ep;
+   MaterialProperty<Real> & _sigma0;
+   MaterialProperty<Real> & _sigma1;
+   MaterialProperty<Real> & _sigma2;
+
+   Real ft(const std::vector<Real> & intnl) const;  /// tensile strength
+   Real dft(const std::vector<Real> & intnl) const; /// d(tensile strength)/d(intnl)
+   Real fc(const std::vector<Real> & intnl) const;  /// compressive strength
+   Real dfc(const std::vector<Real> & intnl) const; /// d(compressive strength)/d(intnl)
+   Real beta(const std::vector<Real> & intnl) const;
+   Real dbeta0(const std::vector<Real> & intnl) const;
+   Real dbeta1(const std::vector<Real> & intnl) const;
+   void weighfac(const std::vector<Real> & stress_params, Real & wf) const; /// weight factor
+   void dweighfac(const std::vector<Real> & stress_params, Real & wf, std::vector<Real> & dwf) const; /// d(weight factor)/d(stress)
+   Real damageVar(const std::vector<Real> & stress_params, const std::vector<Real> & intnl) const;
+
+  void computeStressParams(const RankTwoTensor & stress,
+                           std::vector<Real> & stress_params) const override;
+
+  std::vector<RankTwoTensor> dstress_param_dstress(const RankTwoTensor & stress) const override;
+
+  std::vector<RankFourTensor> d2stress_param_dstress(const RankTwoTensor & stress) const override;
+
+  void setEffectiveElasticity(const RankFourTensor & Eijkl) override;
+
+  virtual void preReturnMapV(const std::vector<Real> & trial_stress_params,
+                             const RankTwoTensor & stress_trial,
+                             const std::vector<Real> & intnl_old,
+                             const std::vector<Real> & yf,
+                             const RankFourTensor & Eijkl) override;
+
+  virtual void setStressAfterReturnV(const RankTwoTensor & stress_trial,
+                                     const std::vector<Real> & stress_params,
+                                     Real gaE,
+                                     const std::vector<Real> & intnl,
+                                     const yieldAndFlow & smoothed_q,
+                                     const RankFourTensor & Eijkl,
+                                     RankTwoTensor & stress) const override;
+
+  void yieldFunctionValuesV(const std::vector<Real> & stress_params,
+                            const std::vector<Real> & intnl,
+                            std::vector<Real> & yf) const override;
+
+  void computeAllQV(const std::vector<Real> & stress_params,
+                    const std::vector<Real> & intnl,
+                    std::vector<yieldAndFlow> & all_q) const override;
+
+   virtual void
+   flowPotential(const std::vector<Real> & stress_params,
+  							const std::vector<Real> & intnl,
+  							std::vector<Real> & r) const;
+
+   virtual void
+   dflowPotential_dstress(const std::vector<Real> & stress_params,
+  											 const std::vector<Real> & intnl,
+  											 std::vector< std::vector<Real> > & dr_dstress) const;
+
+
+   virtual void
+   dflowPotential_dintnl(const std::vector<Real> & stress_params,
+  											const std::vector<Real> & intnl,
+  											std::vector< std::vector<Real> > & dr_dintnl) const;
+
+	 virtual void
+   hardPotential(const std::vector<Real> & stress_params,
+								 const std::vector<Real> & intnl,
+								 std::vector<Real> & h) const;
+
+   virtual void
+   dhardPotential_dstress(const std::vector<Real> & stress_params,
+   												const std::vector<Real> & intnl,
+   												std::vector<std::vector<Real> > & dh_dsig) const;
+
+   virtual void
+   dhardPotential_dintnl(const std::vector<Real> & stress_params,
+												 const std::vector<Real> & intnl,
+												 std::vector<std::vector<Real> > & dh_dintnl) const;
+
+  void initialiseVarsV(const std::vector<Real> & trial_stress_params,
+                       const std::vector<Real> & intnl_old,
+                       std::vector<Real> & stress_params,
+                       Real & gaE,
+                       std::vector<Real> & intnl) const;
+
+  void setIntnlValuesV(const std::vector<Real> & trial_stress_params,
+                       const std::vector<Real> & current_stress_params,
+                       const std::vector<Real> & intnl_old,
+                       std::vector<Real> & intnl) const override;
+
+  void setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
+                            const std::vector<Real> & current_stress_params,
+                            const std::vector<Real> & intnl,
+                            std::vector<std::vector<Real>> & dintnl) const override;
+
+  virtual void consistentTangentOperatorV(const RankTwoTensor & stress_trial,
+                                          const std::vector<Real> & trial_stress_params,
+                                          const RankTwoTensor & stress,
+                                          const std::vector<Real> & stress_params,
+                                          Real gaE,
+                                          const yieldAndFlow & smoothed_q,
+                                          const RankFourTensor & Eijkl,
+                                          bool compute_full_tangent_operator,
+                                          const std::vector<std::vector<Real>> & dvar_dtrial,
+                                          RankFourTensor & cto) override;
+};
+
+#endif // PLASTICDAMAGESTRESSUPDATE_H

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<ComputeMultipleInelasticDamageStress>()
+{
+  InputParameters params = validParams<ComputeMultipleInelasticStress>();
+  return params;
+}
+
+ComputeMultipleInelasticDamageStress::ComputeMultipleInelasticDamageStress(const InputParameters & parameters)
+  : ComputeMultipleInelasticStress(parameters),
+    _D(getMaterialProperty<Real>("Damage_Variable")),
+    _D_old(getMaterialPropertyOld<Real>("Damage_Variable")),
+    _D_older(getMaterialPropertyOlder<Real>("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]);
+}