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