www.mooseframework.org
CappedWeakInclinedPlaneStressUpdate.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "RotationMatrix.h" // for rotVecToZ
12 #include "libmesh/utility.h"
13 
15 
18 {
20  params.addClassDescription("Capped weak inclined plane plasticity stress calculator");
21  params.addRequiredParam<RealVectorValue>("normal_vector", "The normal vector to the weak plane");
22  return params;
23 }
24 
26  const InputParameters & parameters)
27  : CappedWeakPlaneStressUpdate(parameters),
28  _n_input(getParam<RealVectorValue>("normal_vector")),
29  _n(declareProperty<RealVectorValue>("weak_plane_normal")),
30  _n_old(getMaterialProperty<RealVectorValue>("weak_plane_normal")),
31  _rot_n_to_z(RealTensorValue()),
32  _rot_z_to_n(RealTensorValue()),
33  _rotated_trial(RankTwoTensor()),
34  _rotated_Eijkl(RankFourTensor())
35 {
36  if (_n_input.norm() == 0)
37  mooseError("CappedWeakInclinedPlaneStressUpdate: normal_vector must not have zero length");
38  else
39  _n_input /= _n_input.norm();
40 
42  _rot_z_to_n = _rot_n_to_z.transpose();
43 }
44 
45 void
47 {
49  _n[_qp] = _n_input;
50 }
51 
52 void
54 {
57  for (const auto i : make_range(Moose::dim))
58  {
59  _n[_qp](i) = 0.0;
60  for (const auto j : make_range(Moose::dim))
61  _n[_qp](i) += rotation_increment(i, j) * _n_old[_qp](j);
62  }
63 }
64 
65 void
67 {
70  {
72  _rot_z_to_n = _rot_n_to_z.transpose();
73  }
74 }
75 
76 void
78  Real q_trial,
79  const RankTwoTensor & stress_trial,
80  const std::vector<Real> & /*intnl_old*/,
81  const std::vector<Real> & yf,
82  const RankFourTensor & Eijkl)
83 {
84  // If it's obvious, then simplify the return-type
85  if (yf[1] >= 0)
87  else if (yf[2] >= 0)
89 
90  _rotated_trial = stress_trial;
94  _in_q_trial = q_trial;
95 
96  _rotated_Eijkl = Eijkl;
98 }
99 
100 void
102  Real & p,
103  Real & q) const
104 {
105  RankTwoTensor rotated_stress = stress;
106  rotated_stress.rotate(_rot_n_to_z);
107  p = rotated_stress(2, 2);
108  q = std::sqrt(Utility::pow<2>(rotated_stress(0, 2)) + Utility::pow<2>(rotated_stress(1, 2)));
109 }
110 
111 void
113  Real & Epp,
114  Real & Eqq) const
115 {
116  Epp = _rotated_Eijkl(2, 2, 2, 2);
117  Eqq = _rotated_Eijkl(0, 2, 0, 2);
118 }
119 
120 void
122  Real p_ok,
123  Real q_ok,
124  Real gaE,
125  const std::vector<Real> & /*intnl*/,
126  const yieldAndFlow & smoothed_q,
127  const RankFourTensor & /*Eijkl*/,
128  RankTwoTensor & stress) const
129 {
130  // first get stress in the frame where _n points along "z"
131  stress = _rotated_trial;
132  stress(2, 2) = p_ok;
133  // stress_xx and stress_yy are sitting at their trial-stress values
134  // so need to bring them back via Poisson's ratio
135  stress(0, 0) -= _rotated_Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
136  stress(1, 1) -= _rotated_Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
137  if (_in_q_trial == 0.0)
138  stress(2, 0) = stress(2, 1) = stress(0, 2) = stress(1, 2) = 0.0;
139  else
140  {
141  stress(2, 0) = stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
142  stress(2, 1) = stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
143  }
144 
145  // rotate back to the original frame
146  stress.rotate(_rot_z_to_n);
147 }
148 
149 void
151  Real p_trial,
152  Real q_trial,
153  const RankTwoTensor & stress,
154  Real p,
155  Real q,
156  Real gaE,
157  const yieldAndFlow & smoothed_q,
158  const RankFourTensor & Eijkl,
159  bool compute_full_tangent_operator,
160  RankFourTensor & cto) const
161 {
163  p_trial,
164  q_trial,
165  stress,
166  p,
167  q,
168  gaE,
169  smoothed_q,
170  Eijkl,
171  compute_full_tangent_operator,
172  cto);
173 }
174 
177 {
178  RankTwoTensor dpdsig = RankTwoTensor();
179  dpdsig(2, 2) = 1.0;
180  dpdsig.rotate(_rot_z_to_n);
181  return dpdsig;
182 }
183 
186 {
187  RankTwoTensor rotated_stress = stress;
188  rotated_stress.rotate(_rot_n_to_z);
189  RankTwoTensor dqdsig = CappedWeakPlaneStressUpdate::dqdstress(rotated_stress);
190  dqdsig.rotate(_rot_z_to_n);
191  return dqdsig;
192 }
193 
196 {
197  RankTwoTensor rotated_stress = stress;
198  rotated_stress.rotate(_rot_n_to_z);
199  RankFourTensor d2qdsig2 = CappedWeakPlaneStressUpdate::d2qdstress2(rotated_stress);
200  d2qdsig2.rotate(_rot_z_to_n);
201  return d2qdsig2;
202 }
RealTensorValue _rot_n_to_z
Rotation matrix that rotates _n to "z".
RealVectorValue _n_input
User-input value of the normal vector to the weak plane.
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
auto norm() const -> decltype(std::norm(Real()))
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
RankFourTensor _rotated_Eijkl
Elasticity tensor rotated to the frame where _n points along "z".
RankTwoTensor _rotated_trial
Trial stress rotated to the frame where _n points along "z".
virtual void consistentTangentOperator(const RankTwoTensor &stress_trial, Real p_trial, Real q_trial, const RankTwoTensor &stress, Real p, Real q, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, RankFourTensor &cto) const
Calculates the consistent tangent operator.
CappedWeakInclinedPlaneStressUpdate(const InputParameters &parameters)
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
static constexpr std::size_t dim
virtual void consistentTangentOperator(const RankTwoTensor &stress_trial, Real p_trial, Real q_trial, const RankTwoTensor &stress, Real p, Real q, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, RankFourTensor &cto) const override
Calculates the consistent tangent operator.
virtual void setStressAfterReturn(const RankTwoTensor &stress_trial, Real p_ok, Real q_ok, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
Sets stress from the admissible parameters.
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
void addRequiredParam(const std::string &name, const std::string &doc_string)
RealTensorValue _rot_z_to_n
Rotation matrix that rotates "z" to _n.
void rotate(const RankTwoTensorTempl< Real > &R)
unsigned int _qp
TensorValue< Real > RealTensorValue
Real _in_trial02
trial value of stress(0, 2)
virtual void initializeReturnProcess() override
Derived classes may use this to perform calculations before any return-map process is performed...
virtual void preReturnMap(Real p_trial, Real q_trial, const RankTwoTensor &stress_trial, const std::vector< Real > &intnl_old, const std::vector< Real > &yf, const RankFourTensor &Eijkl) override
Derived classes may employ this function to record stuff or do other computations prior to the return...
CappedWeakInclinedPlaneStressUpdate performs the return-map algorithm and associated stress updates f...
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment) override
Derived classes may use this to perform calculations after the return-map process has completed succe...
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
virtual void computePQ(const RankTwoTensor &stress, Real &p, Real &q) const override
Computes p and q, given stress.
Real _Epp
elasticity tensor in p direction
registerMooseObject("SolidMechanicsApp", CappedWeakInclinedPlaneStressUpdate)
void rotate(const TypeTensor< T > &R)
const bool _perform_finite_strain_rotations
Whether to perform finite-strain rotations.
GenericRealTensorValue< is_ad > rotVecToZ(GenericRealVectorValue< is_ad > vec)
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment) override
Derived classes may use this to perform calculations after the return-map process has completed succe...
enum CappedWeakPlaneStressUpdate::StressReturnType _stress_return_type
IntRange< T > make_range(T beg, T end)
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
static InputParameters validParams()
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
virtual void setEppEqq(const RankFourTensor &Eijkl, Real &Epp, Real &Eqq) const override
Set Epp and Eqq based on the elasticity tensor Derived classes must override this.
const MaterialProperty< RealVectorValue > & _n_old
Old value of the normal.
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
CappedWeakPlaneStressUpdate performs the return-map algorithm and associated stress updates for plast...
virtual void initializeReturnProcess() override
Derived classes may use this to perform calculations before any return-map process is performed...
MaterialProperty< RealVectorValue > & _n
Current value of the normal.
virtual RankTwoTensor dpdstress(const RankTwoTensor &stress) const override
d(p)/d(stress) Derived classes must override this
Real _in_trial12
trial value of stress(1, 2)