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 "TensorMechanicsPlasticJ2.h"
11 : #include "RankFourTensor.h"
12 :
13 : registerMooseObject("TensorMechanicsApp", TensorMechanicsPlasticJ2);
14 :
15 : InputParameters
16 126 : TensorMechanicsPlasticJ2::validParams()
17 : {
18 126 : InputParameters params = TensorMechanicsPlasticModel::validParams();
19 252 : params.addRequiredParam<UserObjectName>(
20 : "yield_strength",
21 : "A TensorMechanicsHardening UserObject that defines hardening of the yield strength");
22 378 : params.addRangeCheckedParam<unsigned>(
23 252 : "max_iterations", 10, "max_iterations>0", "Maximum iterations for custom J2 return map");
24 252 : params.addParam<bool>("use_custom_returnMap",
25 252 : true,
26 : "Whether to use the custom returnMap "
27 : "algorithm. Set to true if you are using "
28 : "isotropic elasticity.");
29 252 : params.addParam<bool>("use_custom_cto",
30 252 : true,
31 : "Whether to use the custom consistent tangent "
32 : "operator computations. Set to true if you are "
33 : "using isotropic elasticity.");
34 126 : params.addClassDescription("J2 plasticity, associative, with hardening");
35 :
36 126 : return params;
37 0 : }
38 :
39 63 : TensorMechanicsPlasticJ2::TensorMechanicsPlasticJ2(const InputParameters & parameters)
40 : : TensorMechanicsPlasticModel(parameters),
41 63 : _strength(getUserObject<TensorMechanicsHardeningModel>("yield_strength")),
42 126 : _max_iters(getParam<unsigned>("max_iterations")),
43 126 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
44 189 : _use_custom_cto(getParam<bool>("use_custom_cto"))
45 : {
46 63 : }
47 :
48 : Real
49 3213376 : TensorMechanicsPlasticJ2::yieldFunction(const RankTwoTensor & stress, Real intnl) const
50 : {
51 3213376 : return std::sqrt(3.0 * stress.secondInvariant()) - yieldStrength(intnl);
52 : }
53 :
54 : RankTwoTensor
55 6144 : TensorMechanicsPlasticJ2::dyieldFunction_dstress(const RankTwoTensor & stress, Real /*intnl*/) const
56 : {
57 6144 : Real sII = stress.secondInvariant();
58 6144 : if (sII == 0.0)
59 0 : return RankTwoTensor();
60 : else
61 6144 : return 0.5 * std::sqrt(3.0 / sII) * stress.dsecondInvariant();
62 : }
63 :
64 : Real
65 8352 : TensorMechanicsPlasticJ2::dyieldFunction_dintnl(const RankTwoTensor & /*stress*/, Real intnl) const
66 : {
67 8352 : return -dyieldStrength(intnl);
68 : }
69 :
70 : RankTwoTensor
71 4608 : TensorMechanicsPlasticJ2::flowPotential(const RankTwoTensor & stress, Real intnl) const
72 : {
73 4608 : return dyieldFunction_dstress(stress, intnl);
74 : }
75 :
76 : RankFourTensor
77 3201726 : TensorMechanicsPlasticJ2::dflowPotential_dstress(const RankTwoTensor & stress, Real /*intnl*/) const
78 : {
79 3201726 : Real sII = stress.secondInvariant();
80 3201726 : if (sII == 0)
81 0 : return RankFourTensor();
82 :
83 3201726 : RankFourTensor dfp = 0.5 * std::sqrt(3.0 / sII) * stress.d2secondInvariant();
84 3201726 : Real pre = -0.25 * std::sqrt(3.0) * std::pow(sII, -1.5);
85 3201726 : RankTwoTensor dII = stress.dsecondInvariant();
86 12806904 : for (unsigned i = 0; i < 3; ++i)
87 38420712 : for (unsigned j = 0; j < 3; ++j)
88 115262136 : for (unsigned k = 0; k < 3; ++k)
89 345786408 : for (unsigned l = 0; l < 3; ++l)
90 259339806 : dfp(i, j, k, l) += pre * dII(i, j) * dII(k, l);
91 3201726 : return dfp;
92 : }
93 :
94 : RankTwoTensor
95 8352 : TensorMechanicsPlasticJ2::dflowPotential_dintnl(const RankTwoTensor & /*stress*/,
96 : Real /*intnl*/) const
97 : {
98 8352 : return RankTwoTensor();
99 : }
100 :
101 : Real
102 16065744 : TensorMechanicsPlasticJ2::yieldStrength(Real intnl) const
103 : {
104 16065744 : return _strength.value(intnl);
105 : }
106 :
107 : Real
108 9640338 : TensorMechanicsPlasticJ2::dyieldStrength(Real intnl) const
109 : {
110 9640338 : return _strength.derivative(intnl);
111 : }
112 :
113 : std::string
114 0 : TensorMechanicsPlasticJ2::modelName() const
115 : {
116 0 : return "J2";
117 : }
118 :
119 : bool
120 3215672 : TensorMechanicsPlasticJ2::returnMap(const RankTwoTensor & trial_stress,
121 : Real intnl_old,
122 : const RankFourTensor & E_ijkl,
123 : Real ep_plastic_tolerance,
124 : RankTwoTensor & returned_stress,
125 : Real & returned_intnl,
126 : std::vector<Real> & dpm,
127 : RankTwoTensor & delta_dp,
128 : std::vector<Real> & yf,
129 : bool & trial_stress_inadmissible) const
130 : {
131 3215672 : if (!(_use_custom_returnMap))
132 8480 : return TensorMechanicsPlasticModel::returnMap(trial_stress,
133 : intnl_old,
134 : E_ijkl,
135 : ep_plastic_tolerance,
136 : returned_stress,
137 : returned_intnl,
138 : dpm,
139 : delta_dp,
140 : yf,
141 : trial_stress_inadmissible);
142 :
143 3207192 : yf.resize(1);
144 :
145 3207192 : Real yf_orig = yieldFunction(trial_stress, intnl_old);
146 :
147 3207192 : yf[0] = yf_orig;
148 :
149 3207192 : if (yf_orig < _f_tol)
150 : {
151 : // the trial_stress is admissible
152 7002 : trial_stress_inadmissible = false;
153 7002 : return true;
154 : }
155 :
156 3200190 : trial_stress_inadmissible = true;
157 3200190 : Real mu = E_ijkl(0, 1, 0, 1);
158 :
159 : // Perform a Newton-Raphson to find dpm when
160 : // residual = 3*mu*dpm - trial_equivalent_stress + yieldStrength(intnl_old + dpm) = 0
161 3200190 : Real trial_equivalent_stress = yf_orig + yieldStrength(intnl_old);
162 : Real residual;
163 : Real jac;
164 3200190 : dpm[0] = 0;
165 : unsigned int iter = 0;
166 : do
167 : {
168 6431796 : residual = 3.0 * mu * dpm[0] - trial_equivalent_stress + yieldStrength(intnl_old + dpm[0]);
169 6431796 : jac = 3.0 * mu + dyieldStrength(intnl_old + dpm[0]);
170 6431796 : dpm[0] += -residual / jac;
171 6431796 : if (iter > _max_iters) // not converging
172 : return false;
173 6431796 : iter++;
174 6431796 : } while (residual * residual > _f_tol * _f_tol);
175 :
176 : // set the returned values
177 3200190 : yf[0] = 0;
178 3200190 : returned_intnl = intnl_old + dpm[0];
179 3200190 : RankTwoTensor nn = 1.5 * trial_stress.deviatoric() /
180 : trial_equivalent_stress; // = dyieldFunction_dstress(trial_stress, intnl_old) =
181 : // the normal to the yield surface, at the trial
182 : // stress
183 3200190 : returned_stress = 2.0 / 3.0 * nn * yieldStrength(returned_intnl);
184 3200190 : returned_stress.addIa(1.0 / 3.0 * trial_stress.trace());
185 3200190 : delta_dp = nn * dpm[0];
186 :
187 : return true;
188 : }
189 :
190 : RankFourTensor
191 3200190 : TensorMechanicsPlasticJ2::consistentTangentOperator(const RankTwoTensor & trial_stress,
192 : Real intnl_old,
193 : const RankTwoTensor & stress,
194 : Real intnl,
195 : const RankFourTensor & E_ijkl,
196 : const std::vector<Real> & cumulative_pm) const
197 : {
198 3200190 : if (!_use_custom_cto)
199 0 : return TensorMechanicsPlasticModel::consistentTangentOperator(
200 : trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
201 :
202 3200190 : Real mu = E_ijkl(0, 1, 0, 1);
203 :
204 3200190 : Real h = 3 * mu + dyieldStrength(intnl);
205 3200190 : RankTwoTensor sij = stress.deviatoric();
206 3200190 : Real sII = stress.secondInvariant();
207 3200190 : Real equivalent_stress = std::sqrt(3.0 * sII);
208 3200190 : Real zeta = cumulative_pm[0] / (1.0 + 3.0 * mu * cumulative_pm[0] / equivalent_stress);
209 :
210 6400380 : return E_ijkl - 3.0 * mu * mu / sII / h * sij.outerProduct(sij) -
211 6400380 : 4.0 * mu * mu * zeta * dflowPotential_dstress(stress, intnl);
212 : }
213 :
214 : bool
215 0 : TensorMechanicsPlasticJ2::useCustomReturnMap() const
216 : {
217 0 : return _use_custom_returnMap;
218 : }
219 :
220 : bool
221 3200190 : TensorMechanicsPlasticJ2::useCustomCTO() const
222 : {
223 3200190 : return _use_custom_cto;
224 : }
|
