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 "ComputeSimoHughesJ2PlasticityStress.h"
11 :
12 : registerMooseObject("SolidMechanicsApp", ComputeSimoHughesJ2PlasticityStress);
13 :
14 : InputParameters
15 56 : ComputeSimoHughesJ2PlasticityStress::validParams()
16 : {
17 56 : InputParameters params = DerivativeMaterialInterface<ComputeLagrangianStressPK1>::validParams();
18 56 : params += SingleVariableReturnMappingSolution::validParams();
19 56 : params.addClassDescription("The Simo-Hughes style J2 plasticity.");
20 112 : params.addParam<MaterialPropertyName>(
21 : "elasticity_tensor", "elasticity_tensor", "The name of the elasticity tensor.");
22 112 : params.addRequiredParam<MaterialName>("flow_stress_material",
23 : "The material defining the flow stress");
24 56 : return params;
25 0 : }
26 :
27 42 : ComputeSimoHughesJ2PlasticityStress::ComputeSimoHughesJ2PlasticityStress(
28 42 : const InputParameters & parameters)
29 : : DerivativeMaterialInterface<ComputeLagrangianStressPK1>(parameters),
30 : SingleVariableReturnMappingSolution(parameters),
31 42 : _elasticity_tensor_name(_base_name + getParam<MaterialPropertyName>("elasticity_tensor")),
32 42 : _elasticity_tensor(getMaterialProperty<RankFourTensor>(_elasticity_tensor_name)),
33 84 : _F_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
34 42 : _ep_name(_base_name + "effective_plastic_strain"),
35 42 : _ep(declareProperty<Real>(_ep_name)),
36 42 : _ep_old(getMaterialPropertyOldByName<Real>(_ep_name)),
37 42 : _be(declareProperty<RankTwoTensor>(_base_name +
38 : "volume_preserving_elastic_left_cauchy_green_strain")),
39 42 : _be_old(getMaterialPropertyOldByName<RankTwoTensor>(
40 42 : _base_name + "volume_preserving_elastic_left_cauchy_green_strain")),
41 42 : _Np(declareProperty<RankTwoTensor>(_base_name + "flow_direction")),
42 42 : _flow_stress_material(nullptr),
43 42 : _flow_stress_name(_base_name + "flow_stress"),
44 42 : _H(getMaterialPropertyByName<Real>(_flow_stress_name)),
45 42 : _dH(getDefaultMaterialPropertyByName<Real, false>(
46 168 : derivativePropertyName(_flow_stress_name, {_ep_name}))),
47 42 : _d2H(getDefaultMaterialPropertyByName<Real, false>(
48 252 : derivativePropertyName(_flow_stress_name, {_ep_name, _ep_name})))
49 : {
50 126 : }
51 :
52 : void
53 42 : ComputeSimoHughesJ2PlasticityStress::initialSetup()
54 : {
55 42 : ComputeLagrangianStressBase::initialSetup();
56 :
57 42 : _flow_stress_material = &getMaterial("flow_stress_material");
58 :
59 : // Enforce isotropic elastic tensor
60 42 : if (!hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC))
61 0 : mooseError("ComputeSimoHughesJ2PlasticityStress requires an isotropic elasticity tensor");
62 42 : }
63 :
64 : void
65 160 : ComputeSimoHughesJ2PlasticityStress::initQpStatefulProperties()
66 : {
67 160 : ComputeLagrangianStressPK1::initQpStatefulProperties();
68 160 : _be[_qp].setToIdentity();
69 160 : _ep[_qp] = 0;
70 160 : }
71 :
72 : void
73 160600 : ComputeSimoHughesJ2PlasticityStress::computeQpPK1Stress()
74 : {
75 : usingTensorIndices(i, j, k, l, m);
76 160600 : const Real G = ElasticityTensorTools::getIsotropicShearModulus(_elasticity_tensor[_qp]);
77 : const Real K = ElasticityTensorTools::getIsotropicBulkModulus(_elasticity_tensor[_qp]);
78 : const auto I = RankTwoTensor::Identity();
79 160600 : const auto Fit = _F[_qp].inverse().transpose();
80 160600 : const auto detJ = _F[_qp].det();
81 :
82 : // Update configuration
83 160600 : RankTwoTensor f = _inv_df[_qp].inverse();
84 160600 : RankTwoTensor f_bar = f / std::cbrt(f.det());
85 :
86 : // Elastic predictor
87 160600 : _be[_qp] = f_bar * _be_old[_qp] * f_bar.transpose();
88 160600 : RankTwoTensor s = G * _be[_qp].deviatoric();
89 321008 : _Np[_qp] = MooseUtils::absoluteFuzzyEqual(s.norm(), 0) ? std::sqrt(1. / 2.) * I
90 160600 : : std::sqrt(3. / 2.) * s / s.norm();
91 160600 : Real s_eff = s.doubleContraction(_Np[_qp]);
92 :
93 : // Compute the derivative of the strain before return mapping
94 160600 : if (_fe_problem.currentlyComputingJacobian())
95 59632 : _d_be_d_F = _F_old[_qp].inverse().times<l, m, i, j, k, m>(
96 29816 : (I.times<i, k, j, l>(f_bar * _be_old[_qp].transpose()) +
97 29816 : I.times<j, k, i, l>(f_bar * _be_old[_qp])) /
98 29816 : std::cbrt(f.det()) -
99 89448 : 2. / 3. * _be[_qp].times<i, j, l, k>(_inv_df[_qp]));
100 :
101 : // Check for plastic loading and do return mapping
102 160600 : Real delta_ep = 0;
103 160600 : if (computeResidual(s_eff, 0) > 0)
104 : {
105 : // Initialize the derivative of the internal variable
106 132163 : if (_fe_problem.currentlyComputingJacobian())
107 : {
108 : _d_deltaep_d_betr.zero();
109 25825 : if (MooseUtils::absoluteFuzzyEqual(s.norm(), 0))
110 0 : _d_n_d_be.zero();
111 : else
112 25825 : _d_n_d_be = G / std::sqrt(6) / s.norm() *
113 25825 : (3 * I.times<i, k, j, l>(I) - 2 * _Np[_qp].times<i, j, k, l>(_Np[_qp]) -
114 77475 : I.times<i, j, k, l>(I));
115 : }
116 :
117 132163 : returnMappingSolve(s_eff, delta_ep, _console);
118 :
119 : // Correct the derivative of the strain after return mapping
120 132163 : if (_fe_problem.currentlyComputingJacobian())
121 : _d_be_d_F -=
122 25825 : 2. / 3. *
123 25825 : (_be[_qp].trace() * _Np[_qp].times<i, j, k, l>(_d_deltaep_d_betr) +
124 103300 : delta_ep * _Np[_qp].times<i, j, k, l>(I) + delta_ep * _be[_qp].trace() * _d_n_d_be) *
125 51650 : _d_be_d_F;
126 : }
127 :
128 : // Update intermediate and current configurations
129 160600 : _ep[_qp] = _ep_old[_qp] + delta_ep;
130 160600 : _be[_qp] -= 2. / 3. * delta_ep * _be[_qp].trace() * _Np[_qp];
131 160600 : s = G * _be[_qp].deviatoric();
132 160600 : RankTwoTensor tau = (K * (detJ * detJ - 1) / 2) * I + s;
133 160600 : _pk1_stress[_qp] = tau * Fit;
134 :
135 : // Compute the consistent tangent, i.e. the derivative of the PK1 stress w.r.t. the deformation
136 : // gradient.
137 160600 : if (_fe_problem.currentlyComputingJacobian())
138 : {
139 29816 : RankFourTensor d_tau_d_F = K * detJ * detJ * I.times<i, j, k, l>(Fit) +
140 59632 : G * (_d_be_d_F - I.times<i, j, k, l>(I) * _d_be_d_F / 3);
141 29816 : _pk1_jacobian[_qp] = Fit.times<m, j, i, m, k, l>(d_tau_d_F) - Fit.times<k, j, i, l>(tau * Fit);
142 : }
143 160600 : }
144 :
145 : Real
146 313211 : ComputeSimoHughesJ2PlasticityStress::computeReferenceResidual(const Real & effective_trial_stress,
147 : const Real & scalar)
148 : {
149 313211 : const Real G = ElasticityTensorTools::getIsotropicShearModulus(_elasticity_tensor[_qp]);
150 313211 : return effective_trial_stress - G * scalar * _be[_qp].trace();
151 : }
152 :
153 : Real
154 473811 : ComputeSimoHughesJ2PlasticityStress::computeResidual(const Real & effective_trial_stress,
155 : const Real & scalar)
156 : {
157 473811 : const Real G = ElasticityTensorTools::getIsotropicShearModulus(_elasticity_tensor[_qp]);
158 :
159 : // Update the flow stress
160 473811 : _ep[_qp] = _ep_old[_qp] + scalar;
161 473811 : _flow_stress_material->computePropertiesAtQp(_qp);
162 :
163 473811 : return effective_trial_stress - G * scalar * _be[_qp].trace() - _H[_qp];
164 : }
165 :
166 : Real
167 313211 : ComputeSimoHughesJ2PlasticityStress::computeDerivative(const Real & /*effective_trial_stress*/,
168 : const Real & scalar)
169 : {
170 313211 : const Real G = ElasticityTensorTools::getIsotropicShearModulus(_elasticity_tensor[_qp]);
171 :
172 : // Update the flow stress
173 313211 : _ep[_qp] = _ep_old[_qp] + scalar;
174 313211 : _flow_stress_material->computePropertiesAtQp(_qp);
175 :
176 313211 : return -G * _be[_qp].trace() - _dH[_qp];
177 : }
178 :
179 : void
180 181048 : ComputeSimoHughesJ2PlasticityStress::preStep(const Real & scalar, const Real & R, const Real & J)
181 : {
182 181048 : if (!_fe_problem.currentlyComputingJacobian())
183 151184 : return;
184 :
185 : const auto I = RankTwoTensor::Identity();
186 29864 : const Real G = ElasticityTensorTools::getIsotropicShearModulus(_elasticity_tensor[_qp]);
187 :
188 : // Update the flow stress
189 29864 : _ep[_qp] = _ep_old[_qp] + scalar;
190 29864 : _flow_stress_material->computePropertiesAtQp(_qp);
191 :
192 29864 : _d_R_d_betr =
193 29864 : G * _Np[_qp] - G * scalar * I - (G * _be[_qp].trace() + _dH[_qp]) * _d_deltaep_d_betr;
194 29864 : _d_J_d_betr = -G * I - _d2H[_qp] * _d_deltaep_d_betr;
195 29864 : _d_deltaep_d_betr += -1 / J * _d_R_d_betr + R / J / J * _d_J_d_betr;
196 : }
|
