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 416 : TensileStressUpdate::validParams()
17 : {
18 416 : InputParameters params = MultiParameterPlasticityStressUpdate::validParams();
19 832 : params.addRequiredParam<UserObjectName>(
20 : "tensile_strength",
21 : "A SolidMechanicsHardening UserObject that defines hardening of the tensile strength");
22 832 : params.addParam<bool>("perfect_guess",
23 832 : 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 416 : params.addClassDescription(
29 : "Associative, smoothed, tensile (Rankine) plasticity with hardening/softening");
30 416 : return params;
31 0 : }
32 :
33 312 : TensileStressUpdate::TensileStressUpdate(const InputParameters & parameters)
34 : : MultiParameterPlasticityStressUpdate(parameters, 3, 3, 1),
35 312 : _strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
36 624 : _perfect_guess(getParam<bool>("perfect_guess")),
37 624 : _eigvecs(RankTwoTensor())
38 : {
39 312 : }
40 :
41 : void
42 15128 : 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 15128 : stress.symmetricEigenvalues(stress_params);
47 15128 : }
48 :
49 : std::vector<RankTwoTensor>
50 11912 : TensileStressUpdate::dstress_param_dstress(const RankTwoTensor & stress) const
51 : {
52 : std::vector<Real> sp;
53 : std::vector<RankTwoTensor> dsp;
54 11912 : stress.dsymmetricEigenvalues(sp, dsp);
55 11912 : return dsp;
56 11912 : }
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 0 : }
65 :
66 : void
67 11656 : 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 11656 : stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
75 11656 : }
76 :
77 : void
78 11656 : 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 11656 : stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
88 : // rotate to the original frame
89 11656 : stress = _eigvecs * stress * (_eigvecs.transpose());
90 11656 : }
91 :
92 : void
93 23472 : TensileStressUpdate::yieldFunctionValuesV(const std::vector<Real> & stress_params,
94 : const std::vector<Real> & intnl,
95 : std::vector<Real> & yf) const
96 : {
97 23472 : 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 23472 : yf[0] = stress_params[2] - ts; // use largest eigenvalue
103 23472 : yf[1] = stress_params[1] - ts;
104 23472 : yf[2] = stress_params[0] - ts; // use smallest eigenvalue
105 23472 : }
106 :
107 : void
108 56072 : TensileStressUpdate::computeAllQV(const std::vector<Real> & stress_params,
109 : const std::vector<Real> & intnl,
110 : std::vector<yieldAndFlow> & all_q) const
111 : {
112 56072 : const Real ts = tensile_strength(intnl[0]);
113 56072 : const Real dts = dtensile_strength(intnl[0]);
114 :
115 : // yield functions. See comment in yieldFunctionValuesV
116 56072 : all_q[0].f = stress_params[2] - ts;
117 56072 : all_q[1].f = stress_params[1] - ts;
118 56072 : all_q[2].f = stress_params[0] - ts;
119 :
120 : // d(yield function)/d(stress_params)
121 224288 : for (unsigned yf = 0; yf < _num_yf; ++yf)
122 672864 : for (unsigned a = 0; a < _num_sp; ++a)
123 504648 : all_q[yf].df[a] = 0.0;
124 56072 : all_q[0].df[2] = 1.0;
125 56072 : all_q[1].df[1] = 1.0;
126 56072 : all_q[2].df[0] = 1.0;
127 :
128 : // d(yield function)/d(intnl)
129 224288 : for (unsigned yf = 0; yf < _num_yf; ++yf)
130 168216 : 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 224288 : for (unsigned yf = 0; yf < _num_yf; ++yf)
135 672864 : for (unsigned a = 0; a < _num_sp; ++a)
136 504648 : all_q[yf].dg[a] = all_q[yf].df[a];
137 :
138 : // d(flow potential)/d(stress_params)/d(intnl)
139 224288 : for (unsigned yf = 0; yf < _num_yf; ++yf)
140 672864 : for (unsigned a = 0; a < _num_sp; ++a)
141 504648 : all_q[yf].d2g_di[a][0] = 0.0;
142 :
143 : // d(flow potential)/d(stress_params)/d(stress_params)
144 224288 : for (unsigned yf = 0; yf < _num_yf; ++yf)
145 672864 : for (unsigned a = 0; a < _num_sp; ++a)
146 2018592 : for (unsigned b = 0; b < _num_sp; ++b)
147 1513944 : all_q[yf].d2g[a][b] = 0.0;
148 56072 : }
149 :
150 : void
151 11656 : 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 46624 : for (unsigned a = 0; a < _num_sp; ++a)
156 139872 : for (unsigned b = 0; b < _num_sp; ++b)
157 104904 : _Eij[a][b] = Eijkl(a, a, b, b);
158 11656 : _En = _Eij[2][2];
159 11656 : const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
160 46624 : for (unsigned a = 0; a < _num_sp; ++a)
161 : {
162 34968 : _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
163 69936 : for (unsigned b = 0; b < a; ++b)
164 34968 : _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
165 : }
166 11656 : }
167 :
168 : void
169 11656 : 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 11656 : 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 11656 : const Real ts = tensile_strength(intnl_old[0]);
184 11656 : stress_params[2] = ts; // largest eigenvalue
185 11656 : stress_params[1] = std::min(stress_params[1], ts);
186 11656 : stress_params[0] = std::min(stress_params[0], ts);
187 11656 : gaE = trial_stress_params[2] - stress_params[2];
188 : }
189 11656 : setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl);
190 11656 : }
191 :
192 : void
193 67728 : 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 67728 : intnl[0] = intnl_old[0] + (trial_stress_params[2] - current_stress_params[2]) / _Eij[2][2];
199 67728 : }
200 :
201 : void
202 43776 : 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 43776 : dintnl[0][0] = 0.0;
208 43776 : dintnl[0][1] = 0.0;
209 43776 : dintnl[0][2] = -1.0 / _Eij[2][2];
210 43776 : }
211 :
212 : Real
213 91200 : TensileStressUpdate::tensile_strength(const Real internal_param) const
214 : {
215 91200 : return _strength.value(internal_param);
216 : }
217 :
218 : Real
219 56072 : TensileStressUpdate::dtensile_strength(const Real internal_param) const
220 : {
221 56072 : return _strength.derivative(internal_param);
222 : }
223 :
224 : void
225 256 : 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 256 : cto = elasticity_tensor;
237 256 : 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 256 : RankFourTensor drot_dstress;
247 1024 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
248 3072 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
249 9216 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
250 27648 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
251 82944 : for (unsigned a = 0; a < _num_sp; ++a)
252 : {
253 62208 : if (trial_stress_params[a] == trial_stress_params[j])
254 20736 : continue;
255 41472 : drot_dstress(i, j, k, l) +=
256 41472 : 0.5 * _eigvecs(i, a) *
257 41472 : (_eigvecs(k, a) * _eigvecs(l, j) + _eigvecs(l, a) * _eigvecs(k, j)) /
258 41472 : (trial_stress_params[j] - trial_stress_params[a]);
259 : }
260 :
261 256 : const RankTwoTensor eT = _eigvecs.transpose();
262 :
263 256 : RankFourTensor dstress_dtrial;
264 1024 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
265 3072 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
266 9216 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
267 27648 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
268 82944 : for (unsigned a = 0; a < _num_sp; ++a)
269 62208 : dstress_dtrial(i, j, k, l) +=
270 62208 : drot_dstress(i, a, k, l) * stress_params[a] * eT(a, j) +
271 62208 : _eigvecs(i, a) * stress_params[a] * drot_dstress(j, a, k, l);
272 :
273 256 : const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
274 1024 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
275 3072 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
276 9216 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
277 27648 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
278 82944 : for (unsigned a = 0; a < _num_sp; ++a)
279 248832 : for (unsigned b = 0; b < _num_sp; ++b)
280 186624 : dstress_dtrial(i, j, k, l) +=
281 186624 : _eigvecs(i, a) * dvar_dtrial[a][b] * dsp_trial[b](k, l) * eT(a, j);
282 :
283 256 : cto = dstress_dtrial * elasticity_tensor;
284 256 : }
|
