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 "ADComputeFiniteStrain.h"
11 : #include "RankTwoTensor.h"
12 : #include "RankFourTensor.h"
13 : #include "SymmetricRankTwoTensor.h"
14 : #include "SymmetricRankFourTensor.h"
15 :
16 : #include "libmesh/quadrature.h"
17 : #include "libmesh/utility.h"
18 :
19 : registerMooseObject("SolidMechanicsApp", ADComputeFiniteStrain);
20 : registerMooseObject("SolidMechanicsApp", ADSymmetricFiniteStrain);
21 :
22 : template <typename R2, typename R4>
23 : MooseEnum
24 3328 : ADComputeFiniteStrainTempl<R2, R4>::decompositionType()
25 : {
26 6656 : return MooseEnum("TaylorExpansion EigenSolution", "TaylorExpansion");
27 : }
28 :
29 : template <typename R2, typename R4>
30 : InputParameters
31 3328 : ADComputeFiniteStrainTempl<R2, R4>::validParams()
32 : {
33 : InputParameters params = ADComputeIncrementalStrainBase::validParams();
34 3328 : params.addClassDescription(
35 : "Compute a strain increment and rotation increment for finite strains.");
36 6656 : params.addParam<MooseEnum>("decomposition_method",
37 : ADComputeFiniteStrainTempl<R2, R4>::decompositionType(),
38 : "Methods to calculate the strain and rotation increments");
39 3328 : return params;
40 0 : }
41 :
42 : template <typename R2, typename R4>
43 2496 : ADComputeFiniteStrainTempl<R2, R4>::ADComputeFiniteStrainTempl(const InputParameters & parameters)
44 : : ADComputeIncrementalStrainBaseTempl<R2>(parameters),
45 2496 : _Fhat(this->_fe_problem.getMaxQps()),
46 2496 : _decomposition_method(
47 4992 : this->template getParam<MooseEnum>("decomposition_method").template getEnum<DecompMethod>())
48 : {
49 2496 : }
50 :
51 : template <typename R2, typename R4>
52 : void
53 1492042 : ADComputeFiniteStrainTempl<R2, R4>::computeProperties()
54 : {
55 1492042 : ADRankTwoTensor ave_Fhat;
56 12785918 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
57 : {
58 : // Deformation gradient
59 11293876 : auto A = ADRankTwoTensor::initializeFromRows(
60 11293876 : (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
61 :
62 : // Old Deformation gradient
63 11293876 : auto Fbar = ADRankTwoTensor::initializeFromRows(
64 11293876 : (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);
65 :
66 : // A = gradU - gradUold
67 11293876 : A -= Fbar;
68 :
69 : // Fbar = ( I + gradUold)
70 11293876 : Fbar.addIa(1.0);
71 :
72 : // Incremental deformation gradient _Fhat = I + A Fbar^-1
73 11293876 : _Fhat[_qp] = A * Fbar.inverse();
74 11293876 : _Fhat[_qp].addIa(1.0);
75 :
76 : // Calculate average _Fhat for volumetric locking correction
77 11293876 : if (_volumetric_locking_correction)
78 6617616 : ave_Fhat += _Fhat[_qp] * _JxW[_qp] * _coord[_qp];
79 : }
80 :
81 1492042 : if (_volumetric_locking_correction)
82 827202 : ave_Fhat /= _current_elem_volume;
83 :
84 1492042 : const auto ave_Fhat_det = ave_Fhat.det();
85 12785912 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
86 : {
87 : // Finalize volumetric locking correction
88 11293872 : if (_volumetric_locking_correction)
89 : {
90 : using std::cbrt;
91 13235232 : _Fhat[_qp] *= cbrt(ave_Fhat_det / _Fhat[_qp].det());
92 : }
93 :
94 11293872 : computeQpStrain();
95 : }
96 1492040 : }
97 :
98 : template <typename R2, typename R4>
99 : void
100 13716816 : ADComputeFiniteStrainTempl<R2, R4>::computeQpStrain()
101 : {
102 13716816 : ADR2 total_strain_increment;
103 :
104 : // two ways to calculate these increments: TaylorExpansion(default) or EigenSolution
105 13716816 : computeQpIncrements(total_strain_increment, _rotation_increment[_qp]);
106 :
107 13716814 : _strain_increment[_qp] = total_strain_increment;
108 :
109 : // Remove the eigenstrain increment
110 13716814 : this->subtractEigenstrainIncrementFromStrain(_strain_increment[_qp]);
111 :
112 13716814 : if (_dt > 0)
113 13624226 : _strain_rate[_qp] = _strain_increment[_qp] / _dt;
114 : else
115 92588 : _strain_rate[_qp].zero();
116 :
117 : // Update strain in intermediate configuration
118 13716814 : _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
119 13716814 : _total_strain[_qp] = _total_strain_old[_qp] + total_strain_increment;
120 :
121 : // Rotate strain to current configuration
122 13716814 : _mechanical_strain[_qp].rotate(_rotation_increment[_qp]);
123 13716814 : _total_strain[_qp].rotate(_rotation_increment[_qp]);
124 :
125 13716814 : if (_global_strain)
126 0 : _total_strain[_qp] += (*_global_strain)[_qp];
127 13716814 : }
128 :
129 : template <typename R2, typename R4>
130 : void
131 13716816 : ADComputeFiniteStrainTempl<R2, R4>::computeQpIncrements(ADR2 & total_strain_increment,
132 : ADRankTwoTensor & rotation_increment)
133 : {
134 : using std::sqrt;
135 :
136 13716816 : switch (_decomposition_method)
137 : {
138 13711912 : case DecompMethod::TaylorExpansion:
139 : {
140 : // inverse of _Fhat
141 13711912 : const ADRankTwoTensor invFhat = _Fhat[_qp].inverse();
142 :
143 : // A = I - _Fhat^-1
144 13711912 : ADRankTwoTensor A(ADRankTwoTensor::initIdentity);
145 13711912 : A -= invFhat;
146 :
147 : // Cinv - I = A A^T - (A + A^T);
148 13711912 : ADR2 Cinv_I = ADR2::timesTranspose(A) - ADR2::plusTranspose(A);
149 :
150 : // strain rate D from Taylor expansion, Chat = (-1/2(Chat^-1 - I) + 1/4*(Chat^-1 - I)^2 + ...
151 27423824 : total_strain_increment = -Cinv_I * 0.5 + Cinv_I.square() * 0.25;
152 :
153 : const ADReal a[3] = {invFhat(1, 2) - invFhat(2, 1),
154 : invFhat(2, 0) - invFhat(0, 2),
155 : invFhat(0, 1) - invFhat(1, 0)};
156 :
157 13711912 : const auto q = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) / 4.0;
158 13711912 : const auto trFhatinv_1 = invFhat.trace() - 1.0;
159 27423824 : const auto p = trFhatinv_1 * trFhatinv_1 / 4.0;
160 :
161 : // cos theta_a
162 13711912 : const ADReal C1_squared =
163 27423824 : p + 3.0 * Utility::pow<2>(p) * (1.0 - (p + q)) / Utility::pow<2>(p + q) -
164 41135736 : 2.0 * Utility::pow<3>(p) * (1.0 - (p + q)) / Utility::pow<3>(p + q);
165 13711912 : if (C1_squared <= 0.0)
166 0 : mooseException(
167 : "Cannot take square root of a number less than or equal to zero in the calculation of "
168 : "C1 for the Rashid approximation for the rotation tensor. This zero or negative number "
169 : "may occur when elements become heavily distorted.");
170 :
171 13711912 : const ADReal C1 = sqrt(C1_squared);
172 :
173 : ADReal C2;
174 13711912 : if (q > 0.01)
175 : // (1-cos theta_a)/4q
176 76392 : C2 = (1.0 - C1) / (4.0 * q);
177 : else
178 : // alternate form for small q
179 54745792 : C2 = 0.125 + q * 0.03125 * (Utility::pow<2>(p) - 12.0 * (p - 1.0)) / Utility::pow<2>(p) +
180 13686448 : Utility::pow<2>(q) * (p - 2.0) * (Utility::pow<2>(p) - 10.0 * p + 32.0) /
181 : Utility::pow<3>(p) +
182 : Utility::pow<3>(q) *
183 54745792 : (1104.0 - 992.0 * p + 376.0 * Utility::pow<2>(p) - 72.0 * Utility::pow<3>(p) +
184 13686448 : 5.0 * Utility::pow<4>(p)) /
185 27372896 : (512.0 * Utility::pow<4>(p));
186 :
187 13711914 : const ADReal C3_test =
188 27423824 : (p * q * (3.0 - q) + Utility::pow<3>(p) + Utility::pow<2>(q)) / Utility::pow<3>(p + q);
189 13711912 : if (C3_test <= 0.0)
190 2 : mooseException(
191 : "Cannot take square root of a number less than or equal to zero in the calculation of "
192 : "C3_test for the Rashid approximation for the rotation tensor. This zero or negative "
193 : "number may occur when elements become heavily distorted.");
194 27423820 : const ADReal C3 = 0.5 * sqrt(C3_test); // sin theta_a/(2 sqrt(q))
195 :
196 : // Calculate incremental rotation. Note that this value is the transpose of that from Rashid,
197 : // 93, so we transpose it before storing
198 13711910 : ADRankTwoTensor R_incr;
199 13711910 : R_incr.addIa(C1);
200 54847640 : for (unsigned int i = 0; i < 3; ++i)
201 164542920 : for (unsigned int j = 0; j < 3; ++j)
202 246814380 : R_incr(i, j) += C2 * a[i] * a[j];
203 :
204 13711910 : R_incr(0, 1) += C3 * a[2];
205 13711910 : R_incr(0, 2) -= C3 * a[1];
206 13711910 : R_incr(1, 0) -= C3 * a[2];
207 13711910 : R_incr(1, 2) += C3 * a[0];
208 13711910 : R_incr(2, 0) += C3 * a[1];
209 13711910 : R_incr(2, 1) -= C3 * a[0];
210 :
211 27423822 : rotation_increment = R_incr.transpose();
212 : break;
213 : }
214 :
215 4904 : case DecompMethod::EigenSolution:
216 : {
217 : // Add a small perturbation to F for the case when F=I, which occurs with no deformation,
218 : // which commonly occurs in initialization. The perturbation to F does not affect the computed
219 : // stress, but prevents a singularity in the AD-computed material Jacobian.
220 4904 : if (this->_fe_problem.currentlyComputingJacobian() &&
221 5832 : _Fhat[_qp] == ADRankTwoTensor::Identity())
222 264 : _Fhat[_qp] +=
223 528 : ADRankTwoTensor(0.0, 5.0e-13, 5.0e-13, 5.0e-13, 0.0, 5.0e-13, 5.0e-13, 5.0e-13, 0.0);
224 :
225 4904 : FADR2 Chat = ADR2::transposeTimes(_Fhat[_qp]);
226 4904 : FADR2 Uhat = MathUtils::sqrt(Chat);
227 9808 : rotation_increment = _Fhat[_qp] * Uhat.inverse().template get<ADRankTwoTensor>();
228 9808 : total_strain_increment = MathUtils::log(Uhat).template get<ADR2>();
229 : break;
230 : }
231 :
232 0 : default:
233 0 : mooseError("ADComputeFiniteStrain Error: Pass valid decomposition type: TaylorExpansion or "
234 : "EigenSolution.");
235 : }
236 13716814 : }
237 :
238 : template class ADComputeFiniteStrainTempl<RankTwoTensor, RankFourTensor>;
239 : template class ADComputeFiniteStrainTempl<SymmetricRankTwoTensor, SymmetricRankFourTensor>;
|
