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