Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 :
10 : #include "CappedWeakInclinedPlaneStressUpdate.h"
11 : #include "RotationMatrix.h" // for rotVecToZ
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", CappedWeakInclinedPlaneStressUpdate);
15 :
16 : InputParameters
17 74 : CappedWeakInclinedPlaneStressUpdate::validParams()
18 : {
19 74 : InputParameters params = CappedWeakPlaneStressUpdate::validParams();
20 74 : params.addClassDescription("Capped weak inclined plane plasticity stress calculator");
21 148 : params.addRequiredParam<RealVectorValue>("normal_vector", "The normal vector to the weak plane");
22 74 : return params;
23 0 : }
24 :
25 56 : CappedWeakInclinedPlaneStressUpdate::CappedWeakInclinedPlaneStressUpdate(
26 56 : const InputParameters & parameters)
27 : : CappedWeakPlaneStressUpdate(parameters),
28 56 : _n_input(getParam<RealVectorValue>("normal_vector")),
29 56 : _n(declareProperty<RealVectorValue>("weak_plane_normal")),
30 112 : _n_old(getMaterialProperty<RealVectorValue>("weak_plane_normal")),
31 : _rot_n_to_z(RealTensorValue()),
32 : _rot_z_to_n(RealTensorValue()),
33 56 : _rotated_trial(RankTwoTensor()),
34 112 : _rotated_Eijkl(RankFourTensor())
35 : {
36 56 : if (_n_input.norm() == 0)
37 2 : mooseError("CappedWeakInclinedPlaneStressUpdate: normal_vector must not have zero length");
38 : else
39 54 : _n_input /= _n_input.norm();
40 :
41 54 : _rot_n_to_z = RotationMatrix::rotVecToZ(_n_input);
42 54 : _rot_z_to_n = _rot_n_to_z.transpose();
43 54 : }
44 :
45 : void
46 192 : CappedWeakInclinedPlaneStressUpdate::initQpStatefulProperties()
47 : {
48 192 : CappedWeakPlaneStressUpdate::initQpStatefulProperties();
49 192 : _n[_qp] = _n_input;
50 192 : }
51 :
52 : void
53 480 : CappedWeakInclinedPlaneStressUpdate::finalizeReturnProcess(const RankTwoTensor & rotation_increment)
54 : {
55 480 : CappedWeakPlaneStressUpdate::finalizeReturnProcess(rotation_increment);
56 480 : if (_perform_finite_strain_rotations)
57 0 : for (const auto i : make_range(Moose::dim))
58 : {
59 0 : _n[_qp](i) = 0.0;
60 0 : for (const auto j : make_range(Moose::dim))
61 0 : _n[_qp](i) += rotation_increment(i, j) * _n_old[_qp](j);
62 : }
63 480 : }
64 :
65 : void
66 480 : CappedWeakInclinedPlaneStressUpdate::initializeReturnProcess()
67 : {
68 480 : CappedWeakPlaneStressUpdate::initializeReturnProcess();
69 480 : if (_perform_finite_strain_rotations)
70 : {
71 0 : _rot_n_to_z = RotationMatrix::rotVecToZ(_n[_qp]);
72 0 : _rot_z_to_n = _rot_n_to_z.transpose();
73 : }
74 480 : }
75 :
76 : void
77 480 : CappedWeakInclinedPlaneStressUpdate::preReturnMap(Real /*p_trial*/,
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 480 : if (yf[1] >= 0)
86 384 : _stress_return_type = StressReturnType::no_compression;
87 96 : else if (yf[2] >= 0)
88 0 : _stress_return_type = StressReturnType::no_tension;
89 :
90 480 : _rotated_trial = stress_trial;
91 480 : _rotated_trial.rotate(_rot_n_to_z);
92 480 : _in_trial02 = _rotated_trial(0, 2);
93 480 : _in_trial12 = _rotated_trial(1, 2);
94 480 : _in_q_trial = q_trial;
95 :
96 480 : _rotated_Eijkl = Eijkl;
97 480 : _rotated_Eijkl.rotate(_rot_n_to_z);
98 480 : }
99 :
100 : void
101 672 : CappedWeakInclinedPlaneStressUpdate::computePQ(const RankTwoTensor & stress,
102 : Real & p,
103 : Real & q) const
104 : {
105 672 : RankTwoTensor rotated_stress = stress;
106 672 : rotated_stress.rotate(_rot_n_to_z);
107 672 : p = rotated_stress(2, 2);
108 672 : q = std::sqrt(Utility::pow<2>(rotated_stress(0, 2)) + Utility::pow<2>(rotated_stress(1, 2)));
109 672 : }
110 :
111 : void
112 480 : CappedWeakInclinedPlaneStressUpdate::setEppEqq(const RankFourTensor & /*Eijkl*/,
113 : Real & Epp,
114 : Real & Eqq) const
115 : {
116 480 : Epp = _rotated_Eijkl(2, 2, 2, 2);
117 480 : Eqq = _rotated_Eijkl(0, 2, 0, 2);
118 480 : }
119 :
120 : void
121 480 : CappedWeakInclinedPlaneStressUpdate::setStressAfterReturn(const RankTwoTensor & /*stress_trial*/,
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 480 : stress = _rotated_trial;
132 480 : 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 480 : stress(0, 0) -= _rotated_Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
136 480 : stress(1, 1) -= _rotated_Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
137 480 : if (_in_q_trial == 0.0)
138 0 : stress(2, 0) = stress(2, 1) = stress(0, 2) = stress(1, 2) = 0.0;
139 : else
140 : {
141 480 : stress(2, 0) = stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
142 480 : stress(2, 1) = stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
143 : }
144 :
145 : // rotate back to the original frame
146 480 : stress.rotate(_rot_z_to_n);
147 480 : }
148 :
149 : void
150 0 : CappedWeakInclinedPlaneStressUpdate::consistentTangentOperator(const RankTwoTensor & stress_trial,
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 : {
162 0 : TwoParameterPlasticityStressUpdate::consistentTangentOperator(stress_trial,
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 0 : }
174 :
175 : RankTwoTensor
176 480 : CappedWeakInclinedPlaneStressUpdate::dpdstress(const RankTwoTensor & /*stress*/) const
177 : {
178 480 : RankTwoTensor dpdsig = RankTwoTensor();
179 480 : dpdsig(2, 2) = 1.0;
180 480 : dpdsig.rotate(_rot_z_to_n);
181 480 : return dpdsig;
182 : }
183 :
184 : RankTwoTensor
185 480 : CappedWeakInclinedPlaneStressUpdate::dqdstress(const RankTwoTensor & stress) const
186 : {
187 480 : RankTwoTensor rotated_stress = stress;
188 480 : rotated_stress.rotate(_rot_n_to_z);
189 480 : RankTwoTensor dqdsig = CappedWeakPlaneStressUpdate::dqdstress(rotated_stress);
190 480 : dqdsig.rotate(_rot_z_to_n);
191 480 : return dqdsig;
192 : }
193 :
194 : RankFourTensor
195 0 : CappedWeakInclinedPlaneStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
196 : {
197 0 : RankTwoTensor rotated_stress = stress;
198 0 : rotated_stress.rotate(_rot_n_to_z);
199 0 : RankFourTensor d2qdsig2 = CappedWeakPlaneStressUpdate::d2qdstress2(rotated_stress);
200 0 : d2qdsig2.rotate(_rot_z_to_n);
201 0 : return d2qdsig2;
202 : }
|
