ComputeFiniteStrain.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7
8 #include "ComputeFiniteStrain.h"
9 #include "Assembly.h"
10
12 #include "libmesh/utility.h"
13
14 MooseEnum
16 {
17  return MooseEnum("TaylorExpansion EigenSolution", "TaylorExpansion");
18 }
19
20 template <>
21 InputParameters
23 {
24  InputParameters params = validParams<ComputeIncrementalStrainBase>();
26  "Compute a strain increment and rotation increment for finite strains.");
29  "Methods to calculate the strain and rotation increments");
30  return params;
31 }
32
33 ComputeFiniteStrain::ComputeFiniteStrain(const InputParameters & parameters)
34  : ComputeIncrementalStrainBase(parameters),
35  _Fhat(_fe_problem.getMaxQps()),
36  _decomposition_method(getParam<MooseEnum>("decomposition_method").getEnum<DecompMethod>())
37 {
38 }
39
40 void
42 {
43  RankTwoTensor ave_Fhat;
44  Real ave_dfgrd_det = 0.0;
45  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
46  {
54
57
59  A -= Fbar;
60
61  // Fbar = ( I + gradUold)
63
64  // Incremental deformation gradient _Fhat = I + A Fbar^-1
65  _Fhat[_qp] = A * Fbar.inverse();
67
69  {
70  // Calculate average _Fhat for volumetric locking correction
71  ave_Fhat += _Fhat[_qp] * _JxW[_qp] * _coord[_qp];
72
74  ave_dfgrd_det += _deformation_gradient[_qp].det() * _JxW[_qp] * _coord[_qp];
75  }
76  }
77
79  {
80  // needed for volumetric locking correction
81  ave_Fhat /= _current_elem_volume;
83  ave_dfgrd_det /= _current_elem_volume;
84  }
85
86  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
87  {
89  {
90  // Finalize volumetric locking correction
91  _Fhat[_qp] *= std::cbrt(ave_Fhat.det() / _Fhat[_qp].det());
93  }
94
96  }
97 }
98
99 void
101 {
102  RankTwoTensor total_strain_increment;
103
104  // two ways to calculate these increments: TaylorExpansion(default) or EigenSolution
105  computeQpIncrements(total_strain_increment, _rotation_increment[_qp]);
106
107  _strain_increment[_qp] = total_strain_increment;
108
109  // Remove the eigenstrain increment
111
112  if (_dt > 0)
113  _strain_rate[_qp] = _strain_increment[_qp] / _dt;
114  else
115  _strain_rate[_qp].zero();
116
117  // Update strain in intermediate configuration
119  _total_strain[_qp] = _total_strain_old[_qp] + total_strain_increment;
120
121  // Rotate strain to current configuration
122  _mechanical_strain[_qp] =
123  _rotation_increment[_qp] * _mechanical_strain[_qp] * _rotation_increment[_qp].transpose();
124  _total_strain[_qp] =
125  _rotation_increment[_qp] * _total_strain[_qp] * _rotation_increment[_qp].transpose();
126 }
127
128 void
129 ComputeFiniteStrain::computeQpIncrements(RankTwoTensor & total_strain_increment,
130  RankTwoTensor & rotation_increment)
131 {
132  switch (_decomposition_method)
133  {
135  {
136  // inverse of _Fhat
137  RankTwoTensor invFhat(_Fhat[_qp].inverse());
138
139  // A = I - _Fhat^-1
140  RankTwoTensor A(RankTwoTensor::initIdentity);
141  A -= invFhat;
142
143  // Cinv - I = A A^T - A - A^T;
144  RankTwoTensor Cinv_I = A * A.transpose() - A - A.transpose();
145
146  // strain rate D from Taylor expansion, Chat = (-1/2(Chat^-1 - I) + 1/4*(Chat^-1 - I)^2 + ...
147  total_strain_increment = -Cinv_I * 0.5 + Cinv_I * Cinv_I * 0.25;
148
149  const Real a[3] = {invFhat(1, 2) - invFhat(2, 1),
150  invFhat(2, 0) - invFhat(0, 2),
151  invFhat(0, 1) - invFhat(1, 0)};
152
153  Real q = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) / 4.0;
154  Real trFhatinv_1 = invFhat.trace() - 1.0;
155  const Real p = trFhatinv_1 * trFhatinv_1 / 4.0;
156
157  // cos theta_a
158  const Real C1 =
159  std::sqrt(p + 3.0 * Utility::pow<2>(p) * (1.0 - (p + q)) / Utility::pow<2>(p + q) -
160  2.0 * Utility::pow<3>(p) * (1.0 - (p + q)) / Utility::pow<3>(p + q));
161
162  Real C2;
163  if (q > 0.01)
164  // (1-cos theta_a)/4q
165  C2 = (1.0 - C1) / (4.0 * q);
166  else
167  // alternate form for small q
168  C2 = 0.125 + q * 0.03125 * (Utility::pow<2>(p) - 12.0 * (p - 1.0)) / Utility::pow<2>(p) +
169  Utility::pow<2>(q) * (p - 2.0) * (Utility::pow<2>(p) - 10.0 * p + 32.0) /
170  Utility::pow<3>(p) +
171  Utility::pow<3>(q) * (1104.0 - 992.0 * p + 376.0 * Utility::pow<2>(p) -
172  72.0 * Utility::pow<3>(p) + 5.0 * Utility::pow<4>(p)) /
173  (512.0 * Utility::pow<4>(p));
174  const Real C3 =
175  0.5 * std::sqrt((p * q * (3.0 - q) + Utility::pow<3>(p) + Utility::pow<2>(q)) /
176  Utility::pow<3>(p + q)); // sin theta_a/(2 sqrt(q))
177
178  // Calculate incremental rotation. Note that this value is the transpose of that from Rashid,
179  // 93, so we transpose it before storing
180  RankTwoTensor R_incr;
182  for (unsigned int i = 0; i < 3; ++i)
183  for (unsigned int j = 0; j < 3; ++j)
184  R_incr(i, j) += C2 * a[i] * a[j];
185
186  R_incr(0, 1) += C3 * a[2];
187  R_incr(0, 2) -= C3 * a[1];
188  R_incr(1, 0) -= C3 * a[2];
189  R_incr(1, 2) += C3 * a[0];
190  R_incr(2, 0) += C3 * a[1];
191  R_incr(2, 1) -= C3 * a[0];
192
193  rotation_increment = R_incr.transpose();
194  break;
195  }
196
198  {
199  std::vector<Real> e_value(3);
200  RankTwoTensor e_vector, N1, N2, N3;
201
202  RankTwoTensor Chat = _Fhat[_qp].transpose() * _Fhat[_qp];
203  Chat.symmetricEigenvaluesEigenvectors(e_value, e_vector);
204
205  const Real lambda1 = std::sqrt(e_value[0]);
206  const Real lambda2 = std::sqrt(e_value[1]);
207  const Real lambda3 = std::sqrt(e_value[2]);
208
209  N1.vectorOuterProduct(e_vector.column(0), e_vector.column(0));
210  N2.vectorOuterProduct(e_vector.column(1), e_vector.column(1));
211  N3.vectorOuterProduct(e_vector.column(2), e_vector.column(2));
212
213  RankTwoTensor Uhat = N1 * lambda1 + N2 * lambda2 + N3 * lambda3;
214  RankTwoTensor invUhat(Uhat.inverse());
215
216  rotation_increment = _Fhat[_qp] * invUhat;
217
218  total_strain_increment =
219  N1 * std::log(lambda1) + N2 * std::log(lambda2) + N3 * std::log(lambda3);
220  break;
221  }
222
223  default:
224  mooseError("ComputeFiniteStrain Error: Pass valid decomposition type: TaylorExpansion or "
225  "EigenSolution.");
226  }
227 }
ComputeFiniteStrain(const InputParameters &parameters)
const Real & _current_elem_volume
const MaterialProperty< RankTwoTensor > & _total_strain_old
const DecompMethod _decomposition_method
InputParameters validParams< ComputeFiniteStrain >()
virtual void computeQpStrain()
MaterialProperty< RankTwoTensor > & _strain_increment
MaterialProperty< RankTwoTensor > & _mechanical_strain
virtual void computeQpIncrements(RankTwoTensor &e, RankTwoTensor &r)
static MooseEnum decompositionType()
MaterialProperty< RankTwoTensor > & _strain_rate
void subtractEigenstrainIncrementFromStrain(RankTwoTensor &strain)
const MaterialProperty< RankTwoTensor > & _mechanical_strain_old