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