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