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