Line data Source code
1 : /****************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* BlackBear */
4 : /* */
5 : /* (c) 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by Battelle Energy Alliance, LLC */
9 : /* Under Contract No. DE-AC07-05ID14517 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* See COPYRIGHT for full restrictions */
13 : /****************************************************************/
14 : #include "DamagePlasticityStressUpdate.h"
15 : #include "libmesh/utility.h"
16 :
17 : registerMooseObject("BlackBearApp", DamagePlasticityStressUpdate);
18 :
19 : InputParameters
20 250 : DamagePlasticityStressUpdate::validParams()
21 : {
22 250 : InputParameters params = MultiParameterPlasticityStressUpdate::validParams();
23 500 : params.addParam<Real>(
24 : "yield_function_tolerance",
25 : "If the yield function is less than this amount, the (stress, internal parameters) are "
26 : "deemed admissible. A std::vector of tolerances must be entered for the multi-surface case");
27 750 : params.addRangeCheckedParam<Real>("biaxial_uniaxial_compressive_stress_factor",
28 500 : 0.1,
29 : "biaxial_uniaxial_compressive_stress_factor < 0.5 & "
30 : "biaxial_uniaxial_compressive_stress_factor >= 0",
31 : "Material parameter that relate biaxial and uniaxial "
32 : "compressive strength, i.e., \alfa = (fb0-fc0)/(2*fb0-fc0)");
33 500 : params.addRequiredParam<Real>("dilatancy_factor", "controls the dilation of concrete");
34 750 : params.addRangeCheckedParam<Real>("stiff_recovery_factor",
35 500 : 0.,
36 : "stiff_recovery_factor <= 1. & stiff_recovery_factor >= 0",
37 : "stiffness recovery factor");
38 500 : params.addRangeCheckedParam<Real>(
39 : "ft_ep_slope_factor_at_zero_ep",
40 : "ft_ep_slope_factor_at_zero_ep <= 1 & ft_ep_slope_factor_at_zero_ep >= 0",
41 : "slope of ft vs plastic strain curve at zero plastic strain");
42 :
43 500 : params.addRequiredParam<Real>(
44 : "tensile_damage_at_half_tensile_strength",
45 : "fraction of the elastic recovery slope in tension at 0.5*ft0 after yielding");
46 500 : params.addRangeCheckedParam<Real>("yield_strength_in_tension",
47 : "yield_strength_in_tension >= 0",
48 : "Tensile yield strength of concrete");
49 500 : params.addRangeCheckedParam<Real>("fracture_energy_in_tension",
50 : "fracture_energy_in_tension >= 0",
51 : "Fracture energy of concrete in uniaxial tension");
52 500 : params.addRangeCheckedParam<Real>("yield_strength_in_compression",
53 : "yield_strength_in_compression >= 0",
54 : "Absolute yield compressice strength");
55 500 : params.addRequiredParam<Real>("compressive_damage_at_max_compressive_strength",
56 : "damage at maximum compressive strength");
57 500 : params.addRequiredParam<Real>("maximum_strength_in_compression",
58 : "Absolute maximum compressive strength");
59 500 : params.addRangeCheckedParam<Real>("fracture_energy_in_compression",
60 : "fracture_energy_in_compression >= 0",
61 : "Fracture energy of concrete in uniaxial compression");
62 :
63 500 : params.addRequiredRangeCheckedParam<Real>(
64 : "tip_smoother",
65 : "tip_smoother>=0",
66 : "Smoothing parameter: the cone vertex at mean = cohesion*cot(friction_angle), will be "
67 : "smoothed by the given amount. Typical value is 0.1*cohesion");
68 500 : params.addParam<bool>("perfect_guess",
69 500 : true,
70 : "Provide a guess to the Newton-Raphson proceedure "
71 : "that is the result from perfect plasticity. With "
72 : "severe hardening/softening this may be "
73 : "suboptimal.");
74 250 : params.addClassDescription("Damage Plasticity Model for concrete");
75 250 : return params;
76 0 : }
77 :
78 192 : DamagePlasticityStressUpdate::DamagePlasticityStressUpdate(const InputParameters & parameters)
79 : : MultiParameterPlasticityStressUpdate(parameters, 3, 1, 2),
80 192 : _f_tol(getParam<Real>("yield_function_tol")),
81 :
82 384 : _alfa(getParam<Real>("biaxial_uniaxial_compressive_stress_factor")),
83 384 : _alfa_p(getParam<Real>("dilatancy_factor")),
84 384 : _s0(getParam<Real>("stiff_recovery_factor")),
85 :
86 384 : _Chi(getParam<Real>("ft_ep_slope_factor_at_zero_ep")),
87 384 : _Dt(getParam<Real>("tensile_damage_at_half_tensile_strength")),
88 384 : _ft(getParam<Real>("yield_strength_in_tension")),
89 384 : _FEt(getParam<Real>("fracture_energy_in_tension")),
90 384 : _fyc(getParam<Real>("yield_strength_in_compression")),
91 384 : _Dc(getParam<Real>("compressive_damage_at_max_compressive_strength")),
92 384 : _fc(getParam<Real>("maximum_strength_in_compression")),
93 384 : _FEc(getParam<Real>("fracture_energy_in_compression")),
94 :
95 192 : _ft0(_ft),
96 192 : _fc0(_fyc),
97 192 : _at(1.5 * std::sqrt(1 - _Chi) - 0.5),
98 192 : _ac((2. * (_fc / _fyc) - 1. + 2. * std::sqrt(Utility::pow<2>(_fc / _fyc) - _fc / _fyc))),
99 192 : _zt((1. + _at) / _at),
100 192 : _zc((1. + _ac) / _ac),
101 192 : _dt_bt(log(1. - _Dt) / log((1. + _at - std::sqrt(1. + _at * _at)) / (2. * _at))),
102 192 : _dc_bc(log(1. - _Dc) / log((1. + _ac) / (2. * _ac))),
103 :
104 384 : _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
105 :
106 192 : _sqrt3(std::sqrt(3.)),
107 384 : _perfect_guess(getParam<bool>("perfect_guess")),
108 192 : _eigvecs(RankTwoTensor()),
109 192 : _intnl0(declareProperty<Real>("damage_state_in_tension")),
110 192 : _intnl1(declareProperty<Real>("damage_state_in_compression")),
111 192 : _ele_len(declareProperty<Real>("element_length")),
112 192 : _gt(declareProperty<Real>("fracture_energy_in_tension")),
113 192 : _gc(declareProperty<Real>("fracture_energy_in_compression")),
114 192 : _tD(declareProperty<Real>("tensile_damage")),
115 192 : _cD(declareProperty<Real>("compression_damage")),
116 192 : _D(declareProperty<Real>("damage_variable")),
117 192 : _min_ep(declareProperty<Real>("min_ep")),
118 192 : _mid_ep(declareProperty<Real>("mid_ep")),
119 192 : _max_ep(declareProperty<Real>("max_ep")),
120 192 : _sigma0(declareProperty<Real>("damaged_min_principal_stress")),
121 192 : _sigma1(declareProperty<Real>("damaged_mid_principal_stress")),
122 384 : _sigma2(declareProperty<Real>("damaged_max_principal_stress"))
123 : {
124 192 : }
125 :
126 : void
127 448 : DamagePlasticityStressUpdate::initQpStatefulProperties()
128 : {
129 : /* There are multiple ways to determine the correct element length (or the characteristic length)
130 : used, the following commented lines show several different options. Some other options are
131 : still being considered. In this code, we define the element length as the cube root of the
132 : element volume */
133 :
134 : // if (_current_elem->n_vertices() < 3)
135 : // _ele_len[_qp] = _current_elem->length(0, 1);
136 : // else if (_current_elem->n_vertices() < 5)
137 : // _ele_len[_qp] = (_current_elem->length(0, 1) + _current_elem->length(1, 2)) / 2.;
138 : // else
139 : // _ele_len[_qp] =
140 : // (_current_elem->length(0, 1) + _current_elem->length(1, 2) + _current_elem->length(0, 4))
141 : // / 3.;
142 :
143 448 : _ele_len[_qp] = std::cbrt(_current_elem->volume());
144 :
145 448 : _gt[_qp] = _FEt / _ele_len[_qp];
146 448 : _gc[_qp] = _FEc / _ele_len[_qp];
147 :
148 448 : _min_ep[_qp] = 0.;
149 448 : _mid_ep[_qp] = 0.;
150 448 : _max_ep[_qp] = 0.;
151 448 : _sigma0[_qp] = 0.;
152 448 : _sigma1[_qp] = 0.;
153 448 : _sigma2[_qp] = 0.;
154 448 : _intnl0[_qp] = 0.;
155 448 : _intnl1[_qp] = 0.;
156 448 : _tD[_qp] = 0.;
157 448 : _cD[_qp] = 0.;
158 448 : _D[_qp] = 0.;
159 448 : MultiParameterPlasticityStressUpdate::initQpStatefulProperties();
160 448 : }
161 :
162 : void
163 589530 : DamagePlasticityStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
164 : {
165 : std::vector<Real> eigstrain;
166 589530 : _plastic_strain[_qp].symmetricEigenvalues(eigstrain);
167 589530 : _min_ep[_qp] = eigstrain[0];
168 589530 : _mid_ep[_qp] = eigstrain[1];
169 589530 : _max_ep[_qp] = eigstrain[2];
170 589530 : }
171 :
172 : void
173 1215970 : DamagePlasticityStressUpdate::computeStressParams(const RankTwoTensor & stress,
174 : std::vector<Real> & stress_params) const
175 : {
176 1215970 : stress.symmetricEigenvalues(stress_params);
177 1215970 : }
178 :
179 : std::vector<RankTwoTensor>
180 589530 : DamagePlasticityStressUpdate::dstress_param_dstress(const RankTwoTensor & stress) const
181 : {
182 : std::vector<Real> sp;
183 : std::vector<RankTwoTensor> dsp;
184 589530 : stress.dsymmetricEigenvalues(sp, dsp);
185 589530 : return dsp;
186 : }
187 :
188 : std::vector<RankFourTensor>
189 0 : DamagePlasticityStressUpdate::d2stress_param_dstress(const RankTwoTensor & stress) const
190 : {
191 : std::vector<RankFourTensor> d2;
192 0 : stress.d2symmetricEigenvalues(d2);
193 0 : return d2;
194 : }
195 :
196 : void
197 589530 : DamagePlasticityStressUpdate::setEffectiveElasticity(const RankFourTensor & Eijkl)
198 : {
199 : // Eijkl is required to be isotropic, so we can use the
200 : // frame where stress is diagonal
201 2358120 : for (unsigned a = 0; a < _num_sp; ++a)
202 7074360 : for (unsigned b = 0; b < _num_sp; ++b)
203 5305770 : _Eij[a][b] = Eijkl(a, a, b, b);
204 589530 : _En = _Eij[2][2];
205 589530 : const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
206 2358120 : for (unsigned a = 0; a < _num_sp; ++a)
207 : {
208 1768590 : _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
209 3537180 : for (unsigned b = 0; b < a; ++b)
210 1768590 : _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
211 : }
212 589530 : }
213 :
214 : void
215 589530 : DamagePlasticityStressUpdate::preReturnMapV(const std::vector<Real> & /*trial_stress_params*/,
216 : const RankTwoTensor & stress_trial,
217 : const std::vector<Real> & /*intnl_old*/,
218 : const std::vector<Real> & /*yf*/,
219 : const RankFourTensor & /*Eijkl*/)
220 : {
221 : std::vector<Real> eigvals;
222 589530 : stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
223 589530 : }
224 :
225 : void
226 589530 : DamagePlasticityStressUpdate::setStressAfterReturnV(const RankTwoTensor & /*stress_trial*/,
227 : const std::vector<Real> & stress_params,
228 : Real /*gaE*/,
229 : const std::vector<Real> & intnl,
230 : const yieldAndFlow & /*smoothed_q*/,
231 : const RankFourTensor & /*Eijkl*/,
232 : RankTwoTensor & stress) const
233 : {
234 : // form the diagonal stress
235 589530 : stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
236 : // rotate to the original frame
237 589530 : stress = _eigvecs * stress * (_eigvecs.transpose());
238 : // _dir[_qp] = _eigvecs;
239 589530 : Real D = damageVar(stress_params, intnl);
240 589530 : _sigma0[_qp] = (1. - D) * stress_params[0];
241 589530 : _sigma1[_qp] = (1. - D) * stress_params[1];
242 589530 : _sigma2[_qp] = (1. - D) * stress_params[2];
243 589530 : _intnl0[_qp] = intnl[0];
244 589530 : _intnl1[_qp] = intnl[1];
245 589530 : _D[_qp] = D;
246 589530 : }
247 :
248 : void
249 1215970 : DamagePlasticityStressUpdate::yieldFunctionValuesV(const std::vector<Real> & stress_params,
250 : const std::vector<Real> & intnl,
251 : std::vector<Real> & yf) const
252 : {
253 1215970 : Real I1 = stress_params[0] + stress_params[1] + stress_params[2];
254 1215970 : Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
255 1215970 : .secondInvariant();
256 3647910 : yf[0] = 1. / (1. - _alfa) *
257 2431940 : (_alfa * I1 + _sqrt3 * std::sqrt(J2) +
258 1935834 : beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) -
259 1215970 : fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]);
260 1215970 : }
261 :
262 : void
263 2339744 : DamagePlasticityStressUpdate::computeAllQV(const std::vector<Real> & stress_params,
264 : const std::vector<Real> & intnl,
265 : std::vector<yieldAndFlow> & all_q) const
266 : {
267 2339744 : Real I1 = stress_params[0] + stress_params[1] + stress_params[2];
268 2339744 : Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
269 2339744 : .secondInvariant();
270 2339744 : RankTwoTensor d_sqrt_3J2;
271 :
272 : Real fcbar, dfcbar;
273 2339744 : fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]);
274 2339744 : dfcbar = dfbar_dkappa(_fc0, _ac, 1. - _dc_bc, intnl[1]);
275 :
276 2339744 : if (J2 == 0)
277 0 : d_sqrt_3J2 = RankTwoTensor();
278 : else
279 2339744 : d_sqrt_3J2 = 0.5 * std::sqrt(3.0 / J2) *
280 2339744 : RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
281 2339744 : .dsecondInvariant();
282 :
283 4679488 : all_q[0].f = 1. / (1. - _alfa) *
284 4679488 : (_alfa * I1 + _sqrt3 * std::sqrt(J2) +
285 3739900 : beta(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) -
286 : fcbar;
287 :
288 9358976 : for (unsigned i = 0; i < _num_sp; ++i)
289 7019232 : if (J2 == 0)
290 0 : all_q[0].df[i] =
291 0 : 1. / (1. - _alfa) * (_alfa + beta(intnl) * (stress_params[2] < 0. ? 0. : (i == 2)));
292 : else
293 7019232 : all_q[0].df[i] =
294 7019232 : 1. / (1. - _alfa) *
295 12638308 : (_alfa + d_sqrt_3J2(i, i) + beta(intnl) * (stress_params[2] < 0. ? 0. : (i == 2)));
296 2339744 : all_q[0].df_di[0] =
297 3739900 : 1. / (1. - _alfa) * (dbeta0(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2]));
298 2339744 : all_q[0].df_di[1] =
299 3739900 : 1. / (1. - _alfa) * (dbeta1(intnl) * (stress_params[2] < 0. ? 0. : stress_params[2])) -
300 : dfcbar;
301 :
302 2339744 : flowPotential(stress_params, intnl, all_q[0].dg);
303 2339744 : dflowPotential_dstress(stress_params, all_q[0].d2g);
304 2339744 : dflowPotential_dintnl(stress_params, intnl, all_q[0].d2g_di);
305 2339744 : }
306 :
307 : void
308 10312480 : DamagePlasticityStressUpdate::flowPotential(const std::vector<Real> & stress_params,
309 : const std::vector<Real> & /* intnl */,
310 : std::vector<Real> & r) const
311 : {
312 10312480 : Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
313 10312480 : .secondInvariant();
314 10312480 : RankTwoTensor d_sqrt_2J2;
315 10312480 : if (J2 == 0)
316 0 : d_sqrt_2J2 = RankTwoTensor();
317 : else
318 10312480 : d_sqrt_2J2 = 0.5 * std::sqrt(2.0 / J2) *
319 10312480 : RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
320 10312480 : .dsecondInvariant();
321 :
322 41249920 : for (unsigned int i = 0; i < _num_sp; ++i)
323 30937440 : r[i] = (_alfa_p + d_sqrt_2J2(i, i));
324 10312480 : }
325 :
326 : void
327 4217408 : DamagePlasticityStressUpdate::dflowPotential_dstress(
328 : const std::vector<Real> & stress_params, std::vector<std::vector<Real>> & dr_dstress) const
329 : {
330 4217408 : Real J2 = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
331 4217408 : .secondInvariant();
332 4217408 : RankTwoTensor dII = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0, 0, 0)
333 4217408 : .dsecondInvariant();
334 4217408 : RankTwoTensor d_sqrt_2J2;
335 4217408 : RankFourTensor dfp;
336 : RankTwoTensor d2J2_dsigi_dsigj =
337 4217408 : RankTwoTensor(2. / 3., 2. / 3., 2. / 3., -1. / 3., -1. / 3., -1. / 3.);
338 4217408 : std::vector<Real> dJ2_dsigi(3);
339 16869632 : for (unsigned i = 0; i < 3; ++i)
340 12652224 : dJ2_dsigi[i] =
341 12652224 : (2 * stress_params[i] - stress_params[(i + 1) % 3] - stress_params[(i + 2) % 3]) / 3;
342 :
343 4217408 : if (J2 == 0)
344 : {
345 0 : for (unsigned i = 0; i < 3; ++i)
346 0 : for (unsigned j = 0; j < 3; ++j)
347 0 : dr_dstress[i][j] = 0.0;
348 : }
349 : else
350 : {
351 16869632 : for (unsigned i = 0; i < 3; ++i)
352 50608896 : for (unsigned j = 0; j < 3; ++j)
353 37956672 : dr_dstress[i][j] = 0.5 * (std::sqrt(2.0 / J2) * d2J2_dsigi_dsigj(i, j) -
354 37956672 : (1 / std::sqrt(2)) * 1.0 / std::sqrt(Utility::pow<3>(J2)) *
355 37956672 : dJ2_dsigi[i] * dJ2_dsigi[j]);
356 : }
357 4217408 : }
358 :
359 : void
360 2339744 : DamagePlasticityStressUpdate::dflowPotential_dintnl(
361 : const std::vector<Real> & /* stress_params */,
362 : const std::vector<Real> & /* intnl */,
363 : std::vector<std::vector<Real>> & dr_dintnl) const
364 : {
365 9358976 : for (unsigned i = 0; i < _num_sp; ++i)
366 21057696 : for (unsigned j = 0; j < _num_intnl; ++j)
367 14038464 : dr_dintnl[i][j] = 0.;
368 2339744 : }
369 :
370 : void
371 4217408 : DamagePlasticityStressUpdate::hardPotential(const std::vector<Real> & stress_params,
372 : const std::vector<Real> & intnl,
373 : std::vector<Real> & h) const
374 : {
375 : Real wf, ft, fc;
376 4217408 : weighfac(stress_params, wf);
377 4217408 : std::vector<Real> r(3);
378 4217408 : ft = f(_ft0, _at, intnl[0]);
379 4217408 : fc = f(_fc0, _ac, intnl[1]);
380 4217408 : flowPotential(stress_params, intnl, r);
381 4217408 : h[0] = wf * ft / _gt[_qp] * r[2];
382 4217408 : h[1] = -(1. - wf) * fc / _gc[_qp] * r[0];
383 4217408 : }
384 :
385 : void
386 1877664 : DamagePlasticityStressUpdate::dhardPotential_dstress(const std::vector<Real> & stress_params,
387 : const std::vector<Real> & intnl,
388 : std::vector<std::vector<Real>> & dh_dsig) const
389 : {
390 : Real wf, ft, fc;
391 1877664 : std::vector<Real> dwf(3);
392 1877664 : dweighfac(stress_params, wf, dwf);
393 1877664 : ft = f(_ft0, _at, intnl[0]);
394 1877664 : fc = f(_fc0, _ac, intnl[1]);
395 :
396 1877664 : std::vector<Real> r(3);
397 1877664 : flowPotential(stress_params, intnl, r);
398 1877664 : std::vector<std::vector<Real>> dr_dsig(3, std::vector<Real>(3));
399 1877664 : dflowPotential_dstress(stress_params, dr_dsig);
400 :
401 7510656 : for (unsigned i = 0; i < _num_sp; ++i)
402 : {
403 5632992 : dh_dsig[0][i] = (wf * dr_dsig[2][i] + dwf[i] * r[2]) * ft / _gt[_qp];
404 5632992 : dh_dsig[1][i] = -((1. - wf) * dr_dsig[0][i] - dwf[i] * r[0]) * fc / _gc[_qp];
405 : }
406 3755328 : }
407 :
408 : void
409 1877664 : DamagePlasticityStressUpdate::dhardPotential_dintnl(
410 : const std::vector<Real> & stress_params,
411 : const std::vector<Real> & intnl,
412 : std::vector<std::vector<Real>> & dh_dintnl) const
413 : {
414 : Real wf, dft, dfc;
415 1877664 : weighfac(stress_params, wf);
416 :
417 1877664 : dft = df_dkappa(_ft0, _at, intnl[0]);
418 1877664 : dfc = df_dkappa(_fc0, _ac, intnl[1]);
419 :
420 1877664 : std::vector<Real> r(3);
421 1877664 : flowPotential(stress_params, intnl, r);
422 :
423 1877664 : dh_dintnl[0][0] = wf * dft / _gt[_qp] * r[2];
424 1877664 : dh_dintnl[0][1] = 0.;
425 1877664 : dh_dintnl[1][0] = 0.;
426 1877664 : dh_dintnl[1][1] = -(1 - wf) * dfc / _gc[_qp] * r[0];
427 1877664 : }
428 :
429 : void
430 0 : DamagePlasticityStressUpdate::initialiseVarsV(const std::vector<Real> & trial_stress_params,
431 : const std::vector<Real> & intnl_old,
432 : std::vector<Real> & stress_params,
433 : Real & /* gaE */,
434 : std::vector<Real> & intnl) const
435 : {
436 0 : setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl);
437 0 : }
438 :
439 : void
440 2339744 : DamagePlasticityStressUpdate::setIntnlValuesV(const std::vector<Real> & trial_stress_params,
441 : const std::vector<Real> & current_stress_params,
442 : const std::vector<Real> & intnl_old,
443 : std::vector<Real> & intnl) const
444 : {
445 : RankTwoTensor sigma_trial = RankTwoTensor(
446 2339744 : trial_stress_params[0], trial_stress_params[1], trial_stress_params[2], 0, 0, 0);
447 2339744 : Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2];
448 2339744 : RankTwoTensor Identity_tensor = RankTwoTensor(1, 1, 1, 0, 0, 0);
449 2339744 : RankTwoTensor sigmadev_trial = sigma_trial - (I1_trial / 3.) * Identity_tensor;
450 2339744 : Real norm_sigmadev_trial = sigmadev_trial.L2norm();
451 :
452 2339744 : Real G = 0.5 * (_Eij[0][0] - _Eij[0][1]); // Lame's mu
453 2339744 : Real K = _Eij[0][1] + 2. * G / 3.; // Bulk modulus
454 2339744 : Real C3 = 3. * K * _alfa_p;
455 :
456 : RankTwoTensor denominator_tensor =
457 2339744 : (2 * G / norm_sigmadev_trial) * sigmadev_trial + C3 * Identity_tensor;
458 :
459 4679488 : RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0],
460 4679488 : trial_stress_params[1] - current_stress_params[1],
461 4679488 : trial_stress_params[2] - current_stress_params[2],
462 4679488 : 0.,
463 4679488 : 0.,
464 2339744 : 0.);
465 :
466 : // Implementing Eqn. (21) of Lee & Fenves, 2001, with backsubstituting eqn. (22)
467 2339744 : Real lam = dsig.L2norm() / denominator_tensor.L2norm();
468 :
469 2339744 : std::vector<Real> h(2);
470 2339744 : hardPotential(current_stress_params, intnl_old, h);
471 :
472 2339744 : intnl[0] = intnl_old[0] + lam * h[0];
473 2339744 : intnl[1] = intnl_old[1] + lam * h[1];
474 2339744 : }
475 :
476 : void
477 1877664 : DamagePlasticityStressUpdate::setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
478 : const std::vector<Real> & current_stress_params,
479 : const std::vector<Real> & intnl,
480 : std::vector<std::vector<Real>> & dintnl) const
481 : {
482 : RankTwoTensor sigma_trial = RankTwoTensor(
483 1877664 : trial_stress_params[0], trial_stress_params[1], trial_stress_params[2], 0, 0, 0);
484 1877664 : Real I1_trial = trial_stress_params[0] + trial_stress_params[1] + trial_stress_params[2];
485 1877664 : RankTwoTensor Identity_tensor = RankTwoTensor(1, 1, 1, 0, 0, 0);
486 1877664 : RankTwoTensor sigmadev_trial = sigma_trial - (I1_trial / 3.) * Identity_tensor;
487 1877664 : Real norm_sigmadev_trial = sigmadev_trial.L2norm();
488 :
489 1877664 : Real G = 0.5 * (_Eij[0][0] - _Eij[0][1]); // Lame's mu
490 1877664 : Real K = _Eij[0][1] + 2. * G / 3.; // Bulk modulus
491 1877664 : Real C3 = 3. * K * _alfa_p;
492 :
493 : RankTwoTensor denominator_tensor =
494 1877664 : 2 * G * (sigmadev_trial / norm_sigmadev_trial) + C3 * Identity_tensor;
495 3755328 : RankTwoTensor dsig = RankTwoTensor(trial_stress_params[0] - current_stress_params[0],
496 3755328 : trial_stress_params[1] - current_stress_params[1],
497 3755328 : trial_stress_params[2] - current_stress_params[2],
498 3755328 : 0.,
499 3755328 : 0.,
500 1877664 : 0.);
501 :
502 1877664 : Real lam = dsig.L2norm() / denominator_tensor.L2norm();
503 :
504 1877664 : std::vector<Real> dlam_dsig(3);
505 7510656 : for (unsigned i = 0; i < _num_sp; ++i)
506 9490674 : dlam_dsig[i] = dsig.L2norm() == 0. ? 0.
507 3857682 : : -(trial_stress_params[i] - current_stress_params[i]) /
508 3857682 : (dsig.L2norm() * denominator_tensor.L2norm());
509 :
510 1877664 : std::vector<Real> h(2);
511 1877664 : hardPotential(current_stress_params, intnl, h);
512 1877664 : std::vector<std::vector<Real>> dh_dsig(2, std::vector<Real>(3));
513 1877664 : dhardPotential_dstress(current_stress_params, intnl, dh_dsig);
514 1877664 : std::vector<std::vector<Real>> dh_dintnl(2, std::vector<Real>(2));
515 1877664 : dhardPotential_dintnl(current_stress_params, intnl, dh_dintnl);
516 :
517 5632992 : for (unsigned i = 0; i < _num_intnl; ++i)
518 15021312 : for (unsigned j = 0; j < _num_sp; ++j)
519 11265984 : dintnl[i][j] = dlam_dsig[j] * h[i] + lam * dh_dsig[i][j];
520 3755328 : }
521 :
522 : Real
523 34064582 : DamagePlasticityStressUpdate::fbar(const Real & f0,
524 : const Real & a,
525 : const Real & exponent,
526 : const Real & kappa) const
527 : {
528 34064582 : Real phi = 1. + a * (2. + a) * kappa;
529 34064582 : Real sqrt_phi = std::sqrt(phi);
530 : Real v = sqrt_phi;
531 34064582 : Real u = (1 + a) / a - sqrt_phi / a;
532 34064582 : Real cal_fbar = f0 * std::pow(u, exponent) * v;
533 34064582 : return (u > 0) ? cal_fbar : 1.E-6; // The minimum value for the strength parameter is 1.E-6, as a
534 : // zero value can cause numerical instabilities due to
535 : // derivatives and use of logarithmic functions.
536 : }
537 :
538 : Real
539 7019232 : DamagePlasticityStressUpdate::dfbar_dkappa(const Real & f0,
540 : const Real & a,
541 : const Real & exponent,
542 : const Real & kappa) const
543 : {
544 7019232 : Real phi = 1. + a * (2. + a) * kappa;
545 : Real dphi_dkappa = a * (2. + a);
546 7019232 : Real sqrt_phi = std::sqrt(phi);
547 : Real v = sqrt_phi;
548 7019232 : Real u = (1 + a) / a - sqrt_phi / a;
549 7019232 : Real dv_dphi = 1. / (2 * v);
550 7019232 : Real du_dphi = -(1 / (2 * a)) * (1 / sqrt_phi);
551 : Real cal_dfbar_dkappa =
552 7019232 : f0 * (std::pow(u, exponent) * dv_dphi + exponent * std::pow(u, exponent - 1) * v * du_dphi) *
553 7019232 : dphi_dkappa;
554 7019232 : return (u > 0) ? cal_dfbar_dkappa : 0.;
555 : }
556 :
557 : Real
558 12190144 : DamagePlasticityStressUpdate::f(const Real & f0, const Real & a, const Real & kappa) const
559 : {
560 12190144 : Real phi = 1. + a * (2. + a) * kappa;
561 12190144 : Real sqrt_phi = std::sqrt(phi);
562 : Real v = phi;
563 12190144 : Real u = (1 + a) * sqrt_phi;
564 12190144 : Real cal_f = (f0 / a) * (u - v);
565 12190144 : return (u > v) ? cal_f : 1.E-6; // The minimum value for the strength parameter is 1.E-6, as a
566 : // zero value can cause numerical instabilities due to derivatives
567 : // and use of logarithmic functions.
568 : }
569 :
570 : Real
571 3755328 : DamagePlasticityStressUpdate::df_dkappa(const Real & f0, const Real & a, const Real & kappa) const
572 : {
573 3755328 : Real phi = 1. + a * (2. + a) * kappa;
574 : Real dphi_dkappa = a * (2. + a);
575 3755328 : Real sqrt_phi = std::sqrt(phi);
576 : Real v = phi;
577 3755328 : Real u = (1 + a) * sqrt_phi;
578 : Real dv_dphi = 1.;
579 3755328 : Real du_dphi = (1 + a) / (2 * sqrt_phi);
580 3755328 : Real cal_df_dkappa = (f0 / a) * (du_dphi - dv_dphi) * dphi_dkappa;
581 3755328 : return (u > v) ? cal_df_dkappa : 0.;
582 : }
583 :
584 : Real
585 10574946 : DamagePlasticityStressUpdate::beta(const std::vector<Real> & intnl) const
586 : {
587 : Real fcbar, ftbar;
588 10574946 : fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]);
589 10574946 : ftbar = fbar(_ft0, _at, 1. - _dt_bt, intnl[0]);
590 :
591 10574946 : return (1. - _alfa) * (fcbar / ftbar) - (1. + _alfa);
592 : }
593 :
594 : Real
595 2339744 : DamagePlasticityStressUpdate::dbeta0(const std::vector<Real> & intnl) const
596 : {
597 : Real fcbar, ftbar, dftbar;
598 2339744 : fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]);
599 2339744 : ftbar = fbar(_ft0, _at, 1. - _dt_bt, intnl[0]);
600 2339744 : dftbar = dfbar_dkappa(_ft0, _at, 1. - _dt_bt, intnl[0]);
601 2339744 : return -(1. - _alfa) * fcbar * dftbar / Utility::pow<2>(ftbar);
602 : }
603 :
604 : Real
605 2339744 : DamagePlasticityStressUpdate::dbeta1(const std::vector<Real> & intnl) const
606 : {
607 : Real fcbar, ftbar, dfcbar;
608 2339744 : fcbar = fbar(_fc0, _ac, 1. - _dc_bc, intnl[1]);
609 2339744 : ftbar = fbar(_ft0, _at, 1. - _dt_bt, intnl[0]);
610 2339744 : dfcbar = dfbar_dkappa(_fc0, _ac, 1. - _dc_bc, intnl[1]);
611 2339744 : return dfcbar / ftbar * (1. - _alfa);
612 : }
613 :
614 : void
615 6684602 : DamagePlasticityStressUpdate::weighfac(const std::vector<Real> & stress_params, Real & wf) const
616 : {
617 : Real Dr = 0.;
618 : Real Nr = 0.;
619 26738408 : for (unsigned i = 0; i < _num_sp; ++i)
620 : {
621 20053806 : if (stress_params[i] > 0.)
622 : {
623 7515726 : Nr += stress_params[i];
624 7515726 : Dr += stress_params[i];
625 : }
626 : else
627 : Dr += -stress_params[i];
628 : }
629 6684602 : wf = Nr / Dr;
630 6684602 : }
631 :
632 : void
633 1877664 : DamagePlasticityStressUpdate::dweighfac(const std::vector<Real> & stress_params,
634 : Real & wf,
635 : std::vector<Real> & dwf) const
636 : {
637 1877664 : std::vector<Real> dNr(3, 0.), dDr(3, 0.);
638 : Real Dr = 0.;
639 : Real Nr = 0.;
640 7510656 : for (unsigned i = 0; i < _num_sp; ++i)
641 : {
642 5632992 : if (stress_params[i] > 0.)
643 : {
644 2643278 : Nr += stress_params[i];
645 2643278 : dNr[i] = 1.;
646 2643278 : Dr += stress_params[i];
647 2643278 : dDr[i] = 1.;
648 : }
649 : else
650 : {
651 : Dr += -stress_params[i];
652 2989714 : dDr[i] = -1.;
653 : }
654 : }
655 1877664 : wf = Nr / Dr;
656 :
657 7510656 : for (unsigned i = 0; i < _num_sp; ++i)
658 5632992 : dwf[i] = (dNr[i] - wf * dDr[i]) / Dr;
659 1877664 : }
660 :
661 : Real
662 589530 : DamagePlasticityStressUpdate::damageVar(const std::vector<Real> & stress_params,
663 : const std::vector<Real> & intnl) const
664 : {
665 589530 : Real sqrtPhi_t = std::sqrt(1. + _at * (2. + _at) * intnl[0]);
666 589530 : if (_zt > sqrtPhi_t / _at)
667 562986 : _tD[_qp] = 1. - std::pow(_zt - sqrtPhi_t / _at, _dt_bt);
668 : else
669 26544 : _tD[_qp] = 1. - 1.E-6;
670 :
671 : Real wf;
672 589530 : weighfac(stress_params, wf);
673 589530 : Real s = _s0 + (1. - _s0) * wf;
674 :
675 : Real sqrtPhi_c;
676 589530 : if (intnl[1] < 1.0)
677 562986 : sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * intnl[1]);
678 : else
679 26544 : sqrtPhi_c = std::sqrt(1. + _ac * (2. + _ac) * 0.99);
680 :
681 589530 : _cD[_qp] = 1. - std::pow((_zc - sqrtPhi_c / _ac), _dc_bc);
682 589530 : return 1. - (1. - s * _tD[_qp]) * (1. - _cD[_qp]);
683 : }
684 :
685 : void
686 100906 : DamagePlasticityStressUpdate::consistentTangentOperatorV(
687 : const RankTwoTensor & /* stress_trial */,
688 : const std::vector<Real> & /* trial_stress_params */,
689 : const RankTwoTensor & /*stress*/,
690 : const std::vector<Real> & /* stress_params */,
691 : Real /*gaE*/,
692 : const yieldAndFlow & /*smoothed_q*/,
693 : const RankFourTensor & elasticity_tensor,
694 : bool /* compute_full_tangent_operator */,
695 : const std::vector<std::vector<Real>> & /* dvar_dtrial */,
696 : RankFourTensor & cto)
697 : {
698 100906 : cto = elasticity_tensor;
699 100906 : return;
700 : }
|
