Line data Source code
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 :
10 : #include "TensileStressUpdate.h"
11 : #include "libmesh/utility.h"
12 :
13 : registerMooseObject("TensorMechanicsApp", TensileStressUpdate);
14 :
15 : InputParameters
16 192 : TensileStressUpdate::validParams()
17 : {
18 192 : InputParameters params = MultiParameterPlasticityStressUpdate::validParams();
19 384 : params.addRequiredParam<UserObjectName>(
20 : "tensile_strength",
21 : "A TensorMechanicsHardening UserObject that defines hardening of the tensile strength");
22 384 : params.addParam<bool>("perfect_guess",
23 384 : true,
24 : "Provide a guess to the Newton-Raphson procedure "
25 : "that is the result from perfect plasticity. With "
26 : "severe hardening/softening this may be "
27 : "suboptimal.");
28 192 : params.addClassDescription(
29 : "Associative, smoothed, tensile (Rankine) plasticity with hardening/softening");
30 192 : return params;
31 0 : }
32 :
33 144 : TensileStressUpdate::TensileStressUpdate(const InputParameters & parameters)
34 : : MultiParameterPlasticityStressUpdate(parameters, 3, 3, 1),
35 144 : _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
36 288 : _perfect_guess(getParam<bool>("perfect_guess")),
37 288 : _eigvecs(RankTwoTensor())
38 : {
39 144 : }
40 :
41 : void
42 7936 : TensileStressUpdate::computeStressParams(const RankTwoTensor & stress,
43 : std::vector<Real> & stress_params) const
44 : {
45 : // stress_params[0] = smallest eigenvalue, stress_params[2] = largest eigenvalue
46 7936 : stress.symmetricEigenvalues(stress_params);
47 7936 : }
48 :
49 : std::vector<RankTwoTensor>
50 5568 : TensileStressUpdate::dstress_param_dstress(const RankTwoTensor & stress) const
51 : {
52 : std::vector<Real> sp;
53 : std::vector<RankTwoTensor> dsp;
54 5568 : stress.dsymmetricEigenvalues(sp, dsp);
55 5568 : return dsp;
56 : }
57 :
58 : std::vector<RankFourTensor>
59 0 : TensileStressUpdate::d2stress_param_dstress(const RankTwoTensor & stress) const
60 : {
61 : std::vector<RankFourTensor> d2;
62 0 : stress.d2symmetricEigenvalues(d2);
63 0 : return d2;
64 : }
65 :
66 : void
67 5440 : TensileStressUpdate::preReturnMapV(const std::vector<Real> & /*trial_stress_params*/,
68 : const RankTwoTensor & stress_trial,
69 : const std::vector<Real> & /*intnl_old*/,
70 : const std::vector<Real> & /*yf*/,
71 : const RankFourTensor & /*Eijkl*/)
72 : {
73 : std::vector<Real> eigvals;
74 5440 : stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
75 5440 : }
76 :
77 : void
78 5440 : TensileStressUpdate::setStressAfterReturnV(const RankTwoTensor & /*stress_trial*/,
79 : const std::vector<Real> & stress_params,
80 : Real /*gaE*/,
81 : const std::vector<Real> & /*intnl*/,
82 : const yieldAndFlow & /*smoothed_q*/,
83 : const RankFourTensor & /*Eijkl*/,
84 : RankTwoTensor & stress) const
85 : {
86 : // form the diagonal stress
87 5440 : stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
88 : // rotate to the original frame
89 5440 : stress = _eigvecs * stress * (_eigvecs.transpose());
90 5440 : }
91 :
92 : void
93 12080 : TensileStressUpdate::yieldFunctionValuesV(const std::vector<Real> & stress_params,
94 : const std::vector<Real> & intnl,
95 : std::vector<Real> & yf) const
96 : {
97 12080 : const Real ts = tensile_strength(intnl[0]);
98 : // The smoothing strategy means that the last yield function
99 : // gives the smallest value when returning to a line/point where,
100 : // without smoothing, the yield functions are equal.
101 : // Therefore, i use the smallest eigenvalue in this yield function
102 12080 : yf[0] = stress_params[2] - ts; // use largest eigenvalue
103 12080 : yf[1] = stress_params[1] - ts;
104 12080 : yf[2] = stress_params[0] - ts; // use smallest eigenvalue
105 12080 : }
106 :
107 : void
108 26720 : TensileStressUpdate::computeAllQV(const std::vector<Real> & stress_params,
109 : const std::vector<Real> & intnl,
110 : std::vector<yieldAndFlow> & all_q) const
111 : {
112 26720 : const Real ts = tensile_strength(intnl[0]);
113 26720 : const Real dts = dtensile_strength(intnl[0]);
114 :
115 : // yield functions. See comment in yieldFunctionValuesV
116 26720 : all_q[0].f = stress_params[2] - ts;
117 26720 : all_q[1].f = stress_params[1] - ts;
118 26720 : all_q[2].f = stress_params[0] - ts;
119 :
120 : // d(yield function)/d(stress_params)
121 106880 : for (unsigned yf = 0; yf < _num_yf; ++yf)
122 320640 : for (unsigned a = 0; a < _num_sp; ++a)
123 240480 : all_q[yf].df[a] = 0.0;
124 26720 : all_q[0].df[2] = 1.0;
125 26720 : all_q[1].df[1] = 1.0;
126 26720 : all_q[2].df[0] = 1.0;
127 :
128 : // d(yield function)/d(intnl)
129 106880 : for (unsigned yf = 0; yf < _num_yf; ++yf)
130 80160 : all_q[yf].df_di[0] = -dts;
131 :
132 : // the flow potential is just the yield function
133 : // d(flow potential)/d(stress_params)
134 106880 : for (unsigned yf = 0; yf < _num_yf; ++yf)
135 320640 : for (unsigned a = 0; a < _num_sp; ++a)
136 240480 : all_q[yf].dg[a] = all_q[yf].df[a];
137 :
138 : // d(flow potential)/d(stress_params)/d(intnl)
139 106880 : for (unsigned yf = 0; yf < _num_yf; ++yf)
140 320640 : for (unsigned a = 0; a < _num_sp; ++a)
141 240480 : all_q[yf].d2g_di[a][0] = 0.0;
142 :
143 : // d(flow potential)/d(stress_params)/d(stress_params)
144 106880 : for (unsigned yf = 0; yf < _num_yf; ++yf)
145 320640 : for (unsigned a = 0; a < _num_sp; ++a)
146 961920 : for (unsigned b = 0; b < _num_sp; ++b)
147 721440 : all_q[yf].d2g[a][b] = 0.0;
148 26720 : }
149 :
150 : void
151 5440 : TensileStressUpdate::setEffectiveElasticity(const RankFourTensor & Eijkl)
152 : {
153 : // Eijkl is required to be isotropic, so we can use the
154 : // frame where stress is diagonal
155 21760 : for (unsigned a = 0; a < _num_sp; ++a)
156 65280 : for (unsigned b = 0; b < _num_sp; ++b)
157 48960 : _Eij[a][b] = Eijkl(a, a, b, b);
158 5440 : _En = _Eij[2][2];
159 5440 : const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
160 21760 : for (unsigned a = 0; a < _num_sp; ++a)
161 : {
162 16320 : _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
163 32640 : for (unsigned b = 0; b < a; ++b)
164 16320 : _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
165 : }
166 5440 : }
167 :
168 : void
169 5440 : TensileStressUpdate::initializeVarsV(const std::vector<Real> & trial_stress_params,
170 : const std::vector<Real> & intnl_old,
171 : std::vector<Real> & stress_params,
172 : Real & gaE,
173 : std::vector<Real> & intnl) const
174 : {
175 5440 : if (!_perfect_guess)
176 : {
177 0 : for (unsigned i = 0; i < _num_sp; ++i)
178 0 : stress_params[i] = trial_stress_params[i];
179 0 : gaE = 0.0;
180 : }
181 : else
182 : {
183 5440 : const Real ts = tensile_strength(intnl_old[0]);
184 5440 : stress_params[2] = ts; // largest eigenvalue
185 5440 : stress_params[1] = std::min(stress_params[1], ts);
186 5440 : stress_params[0] = std::min(stress_params[0], ts);
187 5440 : gaE = trial_stress_params[2] - stress_params[2];
188 : }
189 5440 : setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl);
190 5440 : }
191 :
192 : void
193 32160 : TensileStressUpdate::setIntnlValuesV(const std::vector<Real> & trial_stress_params,
194 : const std::vector<Real> & current_stress_params,
195 : const std::vector<Real> & intnl_old,
196 : std::vector<Real> & intnl) const
197 : {
198 32160 : intnl[0] = intnl_old[0] + (trial_stress_params[2] - current_stress_params[2]) / _Eij[2][2];
199 32160 : }
200 :
201 : void
202 20944 : TensileStressUpdate::setIntnlDerivativesV(const std::vector<Real> & /*trial_stress_params*/,
203 : const std::vector<Real> & /*current_stress_params*/,
204 : const std::vector<Real> & /*intnl*/,
205 : std::vector<std::vector<Real>> & dintnl) const
206 : {
207 20944 : dintnl[0][0] = 0.0;
208 20944 : dintnl[0][1] = 0.0;
209 20944 : dintnl[0][2] = -1.0 / _Eij[2][2];
210 20944 : }
211 :
212 : Real
213 44240 : TensileStressUpdate::tensile_strength(const Real internal_param) const
214 : {
215 44240 : return _strength.value(internal_param);
216 : }
217 :
218 : Real
219 26720 : TensileStressUpdate::dtensile_strength(const Real internal_param) const
220 : {
221 26720 : return _strength.derivative(internal_param);
222 : }
223 :
224 : void
225 128 : TensileStressUpdate::consistentTangentOperatorV(const RankTwoTensor & stress_trial,
226 : const std::vector<Real> & trial_stress_params,
227 : const RankTwoTensor & /*stress*/,
228 : const std::vector<Real> & stress_params,
229 : Real /*gaE*/,
230 : const yieldAndFlow & /*smoothed_q*/,
231 : const RankFourTensor & elasticity_tensor,
232 : bool compute_full_tangent_operator,
233 : const std::vector<std::vector<Real>> & dvar_dtrial,
234 : RankFourTensor & cto)
235 : {
236 128 : cto = elasticity_tensor;
237 128 : if (!compute_full_tangent_operator)
238 0 : return;
239 :
240 : // dvar_dtrial has been computed already, so
241 : // d(stress)/d(trial_stress) = d(eigvecs * stress_params * eigvecs.transpose())/d(trial_stress)
242 : // eigvecs is a rotation matrix, rot(i, j) = e_j(i) = i^th component of j^th eigenvector
243 : // d(rot_ij)/d(stress_kl) = d(e_j(i))/d(stress_kl)
244 : // = sum_a 0.5 * e_a(i) * (e_a(k)e_j(l) + e_a(l)e_j(k)) / (la_j - la_a)
245 : // = sum_a 0.5 * rot(i,a) * (rot(k,a)rot(l,j) + rot(l,a)*rot(k,j)) / (la_j - la_a)
246 128 : RankFourTensor drot_dstress;
247 512 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
248 1536 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
249 4608 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
250 13824 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
251 41472 : for (unsigned a = 0; a < _num_sp; ++a)
252 : {
253 31104 : if (trial_stress_params[a] == trial_stress_params[j])
254 10368 : continue;
255 20736 : drot_dstress(i, j, k, l) +=
256 20736 : 0.5 * _eigvecs(i, a) *
257 20736 : (_eigvecs(k, a) * _eigvecs(l, j) + _eigvecs(l, a) * _eigvecs(k, j)) /
258 20736 : (trial_stress_params[j] - trial_stress_params[a]);
259 : }
260 :
261 128 : const RankTwoTensor eT = _eigvecs.transpose();
262 :
263 128 : RankFourTensor dstress_dtrial;
264 512 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
265 1536 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
266 4608 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
267 13824 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
268 41472 : for (unsigned a = 0; a < _num_sp; ++a)
269 31104 : dstress_dtrial(i, j, k, l) +=
270 31104 : drot_dstress(i, a, k, l) * stress_params[a] * eT(a, j) +
271 31104 : _eigvecs(i, a) * stress_params[a] * drot_dstress(j, a, k, l);
272 :
273 128 : const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
274 512 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
275 1536 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
276 4608 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
277 13824 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
278 41472 : for (unsigned a = 0; a < _num_sp; ++a)
279 124416 : for (unsigned b = 0; b < _num_sp; ++b)
280 93312 : dstress_dtrial(i, j, k, l) +=
281 93312 : _eigvecs(i, a) * dvar_dtrial[a][b] * dsp_trial[b](k, l) * eT(a, j);
282 :
283 128 : cto = dstress_dtrial * elasticity_tensor;
284 : }
|
