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