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 282 : SolidMechanicsPlasticJ2::validParams()
21 : {
22 282 : InputParameters params = SolidMechanicsPlasticModel::validParams();
23 564 : params.addRequiredParam<UserObjectName>(
24 : "yield_strength",
25 : "A SolidMechanicsHardening UserObject that defines hardening of the yield strength");
26 846 : params.addRangeCheckedParam<unsigned>(
27 564 : "max_iterations", 10, "max_iterations>0", "Maximum iterations for custom J2 return map");
28 564 : params.addParam<bool>("use_custom_returnMap",
29 564 : true,
30 : "Whether to use the custom returnMap "
31 : "algorithm. Set to true if you are using "
32 : "isotropic elasticity.");
33 564 : params.addParam<bool>("use_custom_cto",
34 564 : true,
35 : "Whether to use the custom consistent tangent "
36 : "operator computations. Set to true if you are "
37 : "using isotropic elasticity.");
38 282 : params.addClassDescription("J2 plasticity, associative, with hardening");
39 :
40 282 : return params;
41 0 : }
42 :
43 141 : SolidMechanicsPlasticJ2::SolidMechanicsPlasticJ2(const InputParameters & parameters)
44 : : SolidMechanicsPlasticModel(parameters),
45 141 : _strength(getUserObject<SolidMechanicsHardeningModel>("yield_strength")),
46 282 : _max_iters(getParam<unsigned>("max_iterations")),
47 282 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
48 423 : _use_custom_cto(getParam<bool>("use_custom_cto"))
49 : {
50 141 : }
51 :
52 : Real
53 6426488 : SolidMechanicsPlasticJ2::yieldFunction(const RankTwoTensor & stress, Real intnl) const
54 : {
55 6426488 : return std::sqrt(3.0 * stress.secondInvariant()) - yieldStrength(intnl);
56 : }
57 :
58 : RankTwoTensor
59 18432 : SolidMechanicsPlasticJ2::dyieldFunction_dstress(const RankTwoTensor & stress, Real /*intnl*/) const
60 : {
61 18432 : Real sII = stress.secondInvariant();
62 18432 : if (sII == 0.0)
63 0 : return RankTwoTensor();
64 : else
65 18432 : return 0.5 * std::sqrt(3.0 / sII) * stress.dsecondInvariant();
66 : }
67 :
68 : Real
69 23176 : SolidMechanicsPlasticJ2::dyieldFunction_dintnl(const RankTwoTensor & /*stress*/, Real intnl) const
70 : {
71 23176 : return -dyieldStrength(intnl);
72 : }
73 :
74 : RankTwoTensor
75 13824 : SolidMechanicsPlasticJ2::flowPotential(const RankTwoTensor & stress, Real intnl) const
76 : {
77 13824 : return dyieldFunction_dstress(stress, intnl);
78 : }
79 :
80 : RankFourTensor
81 6404600 : SolidMechanicsPlasticJ2::dflowPotential_dstress(const RankTwoTensor & stress, Real /*intnl*/) const
82 : {
83 6404600 : Real sII = stress.secondInvariant();
84 6404600 : if (sII == 0)
85 0 : return RankFourTensor();
86 :
87 6404600 : RankFourTensor dfp = 0.5 * std::sqrt(3.0 / sII) * stress.d2secondInvariant();
88 6404600 : Real pre = -0.25 * std::sqrt(3.0) * std::pow(sII, -1.5);
89 6404600 : RankTwoTensor dII = stress.dsecondInvariant();
90 25618400 : for (unsigned i = 0; i < 3; ++i)
91 76855200 : for (unsigned j = 0; j < 3; ++j)
92 230565600 : for (unsigned k = 0; k < 3; ++k)
93 691696800 : for (unsigned l = 0; l < 3; ++l)
94 518772600 : dfp(i, j, k, l) += pre * dII(i, j) * dII(k, l);
95 6404600 : return dfp;
96 : }
97 :
98 : RankTwoTensor
99 23176 : SolidMechanicsPlasticJ2::dflowPotential_dintnl(const RankTwoTensor & /*stress*/,
100 : Real /*intnl*/) const
101 : {
102 23176 : return RankTwoTensor();
103 : }
104 :
105 : Real
106 32159224 : SolidMechanicsPlasticJ2::yieldStrength(Real intnl) const
107 : {
108 32159224 : return _strength.value(intnl);
109 : }
110 :
111 : Real
112 19302112 : SolidMechanicsPlasticJ2::dyieldStrength(Real intnl) const
113 : {
114 19302112 : return _strength.derivative(intnl);
115 : }
116 :
117 : std::string
118 0 : SolidMechanicsPlasticJ2::modelName() const
119 : {
120 0 : return "J2";
121 : }
122 :
123 : bool
124 6430368 : 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 6430368 : if (!(_use_custom_returnMap))
136 21280 : 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 6409088 : yf.resize(1);
148 :
149 6409088 : Real yf_orig = yieldFunction(trial_stress, intnl_old);
150 :
151 6409088 : yf[0] = yf_orig;
152 :
153 6409088 : if (yf_orig < _f_tol)
154 : {
155 : // the trial_stress is admissible
156 9096 : trial_stress_inadmissible = false;
157 9096 : return true;
158 : }
159 :
160 6399992 : trial_stress_inadmissible = true;
161 6399992 : 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 6399992 : Real trial_equivalent_stress = yf_orig + yieldStrength(intnl_old);
166 : Real residual;
167 : Real jac;
168 6399992 : dpm[0] = 0;
169 : unsigned int iter = 0;
170 : do
171 : {
172 12878944 : residual = 3.0 * mu * dpm[0] - trial_equivalent_stress + yieldStrength(intnl_old + dpm[0]);
173 12878944 : jac = 3.0 * mu + dyieldStrength(intnl_old + dpm[0]);
174 12878944 : dpm[0] += -residual / jac;
175 12878944 : if (iter > _max_iters) // not converging
176 : return false;
177 12878944 : iter++;
178 12878944 : } while (residual * residual > _f_tol * _f_tol);
179 :
180 : // set the returned values
181 6399992 : yf[0] = 0;
182 6399992 : returned_intnl = intnl_old + dpm[0];
183 6399992 : 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 6399992 : returned_stress = 2.0 / 3.0 * nn * yieldStrength(returned_intnl);
188 6399992 : returned_stress.addIa(1.0 / 3.0 * trial_stress.trace());
189 6399992 : delta_dp = nn * dpm[0];
190 :
191 : return true;
192 : }
193 :
194 : RankFourTensor
195 6399992 : 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 6399992 : if (!_use_custom_cto)
203 0 : return SolidMechanicsPlasticModel::consistentTangentOperator(
204 : trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
205 :
206 6399992 : Real mu = E_ijkl(0, 1, 0, 1);
207 :
208 6399992 : Real h = 3 * mu + dyieldStrength(intnl);
209 6399992 : RankTwoTensor sij = stress.deviatoric();
210 6399992 : Real sII = stress.secondInvariant();
211 6399992 : Real equivalent_stress = std::sqrt(3.0 * sII);
212 6399992 : Real zeta = cumulative_pm[0] / (1.0 + 3.0 * mu * cumulative_pm[0] / equivalent_stress);
213 :
214 12799984 : return E_ijkl - 3.0 * mu * mu / sII / h * sij.outerProduct(sij) -
215 12799984 : 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 6399992 : SolidMechanicsPlasticJ2::useCustomCTO() const
226 : {
227 6399992 : return _use_custom_cto;
228 : }
|
