adding damage plasticity model for concrete ref #184

doc/content/bib/blackbear.bib

@@ -244,3 +244,50 @@ @article{poyet2009temperature
 	  year={2009},
 	  publisher={Elsevier}
 	}
+
+	@article{lee1998plastic,
+  title={Plastic-damage model for cyclic loading of concrete structures},
+  author={Lee, Jeeho and Fenves, Gregory L},
+  journal={Journal of engineering mechanics},
+  volume={124},
+  number={8},
+  pages={892--900},
+  year={1998},
+  publisher={American Society of Civil Engineers}
+}
+
+@book{lee1996theory,
+  title={Theory and implementation of plastic-damage model for concrete structures under cyclic and dynamic loading},
+  author={Lee, Jeeho},
+  year={1996},
+  publisher={University of California, Berkeley}
+}
+
+@article{wilkins2020method,
+  title = {A method for smoothing multiple yield functions},
+ author = {Andy Wilkins and Benjamin W. Spencer and Amit Jain and Bora Gencturk},
+   year = {2020},
+journal = {International Journal for Numerical Methods in Engineering},
+ volume = {121},
+ number = {3},
+  pages = {434--449},
+    doi = {10.1002/nme.6215}
+}
+
+@article{lubliner1989plastic,
+  title={A plastic-damage model for concrete},
+  author={Lubliner, Jacob and Oliver, Javier and Oller, Sand and Onate, EJIJos},
+  journal={International Journal of solids and structures},
+  volume={25},
+  number={3},
+  pages={299--326},
+  year={1989},
+  publisher={Elsevier}
+}
+
+@article{krabbenhoft2002basic,
+  title={Basic computational plasticity},
+  author={Krabbenh{\o}ft, KRISTIAN},
+  journal={University of Denmark},
+  year={2002}
+}

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/damage_plasticity_model/uniaxial_test.i block=Materials/stress
+
+!syntax parameters /Materials/ComputeMultipleInelasticDamageStress
+
+!syntax inputs /Materials/ComputeMultipleInelasticDamageStress
+
+!syntax children /Materials/ComputeMultipleInelasticDamageStress
+
+!bibtex bibliography

include/materials/ComputeMultipleInelasticDamageStress.h

@@ -0,0 +1,42 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/*                       BlackBear                              */
+/*                                                              */
+/*           (c) 2017 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#pragma once
+
+#include "ComputeMultipleInelasticStress.h"
+
+class ComputeMultipleInelasticDamageStress : public ComputeMultipleInelasticStress
+{
+public:
+  static InputParameters validParams();
+  ComputeMultipleInelasticDamageStress(const InputParameters & parameters);
+
+protected:
+  /// damage parameter for DamagePlasticityStressUpdate 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;
+};

include/materials/DamagePlasticityStressUpdate.h

@@ -0,0 +1,184 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/*                       BlackBear                              */
+/*                                                              */
+/*           (c) 2017 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#pragma once
+
+#include "MultiParameterPlasticityStressUpdate.h"
+
+class DamagePlasticityStressUpdate : public MultiParameterPlasticityStressUpdate
+{
+public:
+  static InputParameters validParams();
+  DamagePlasticityStressUpdate(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;
+
+private:
+  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 _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;
+};

src/materials/ComputeMultipleInelasticDamageStress.C

@@ -0,0 +1,73 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/*                       BlackBear                              */
+/*                                                              */
+/*           (c) 2017 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#include "ComputeMultipleInelasticDamageStress.h"
+#include "StressUpdateBase.h"
+
+registerMooseObject("BlackBearApp", ComputeMultipleInelasticDamageStress);
+
+InputParameters
+ComputeMultipleInelasticDamageStress::validParams()
+{
+  InputParameters params = ComputeMultipleInelasticStress::validParams();
+  return params;
+}
+
+ComputeMultipleInelasticDamageStress::ComputeMultipleInelasticDamageStress(
+    const InputParameters & parameters)
+  : ComputeMultipleInelasticStress(parameters),
+    _D(getMaterialProperty<Real>("elemental_damage_variable")),
+    _D_old(getMaterialPropertyOld<Real>("elemental_damage_variable")),
+    _D_older(getMaterialPropertyOlder<Real>("elemental_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];
+}
+
+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_older[_qp]);
+}