www.mooseframework.org
InteractionIntegral.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 // This post processor returns the Interaction Integral
8 //
9 #include "InteractionIntegral.h"
10 #include "MooseMesh.h"
11 #include "RankTwoTensor.h"
12 #include "libmesh/string_to_enum.h"
13 #include "libmesh/quadrature.h"
14 #include "DerivativeMaterialInterface.h"
15 #include "libmesh/utility.h"
16 
17 MooseEnum
19 {
20  return MooseEnum("Geometry Topology", "Geometry");
21 }
22 
23 MooseEnum
25 {
26  return MooseEnum("KI KII KIII T", "KI");
27 }
28 
29 template <>
30 InputParameters
32 {
33  InputParameters params = validParams<ElementIntegralPostprocessor>();
34  params.addRequiredCoupledVar(
35  "displacements",
36  "The displacements appropriate for the simulation geometry and coordinate system");
37  params.addCoupledVar("temp",
38  "The temperature (optional). Must be provided to correctly compute "
39  "stress intensity factors in models with thermal strain gradients.");
40  params.addRequiredParam<UserObjectName>("crack_front_definition",
41  "The CrackFrontDefinition user object name");
42  params.addParam<unsigned int>(
43  "crack_front_point_index",
44  "The index of the point on the crack front corresponding to this q function");
45  params.addParam<Real>(
46  "K_factor", "Conversion factor between interaction integral and stress intensity factor K");
47  params.addParam<unsigned int>("symmetry_plane",
48  "Account for a symmetry plane passing through "
49  "the plane of the crack, normal to the specified "
50  "axis (0=x, 1=y, 2=z)");
51  params.addParam<Real>("poissons_ratio", "Poisson's ratio for the material.");
52  params.addParam<Real>("youngs_modulus", "Young's modulus of the material.");
53  params.set<bool>("use_displaced_mesh") = false;
54  params.addParam<unsigned int>("ring_index", "Ring ID");
55  params.addParam<MooseEnum>("q_function_type",
57  "The method used to define the integration domain. Options are: " +
58  InteractionIntegral::qFunctionType().getRawNames());
59  params.addRequiredParam<MooseEnum>("sif_mode",
61  "Stress intensity factor to calculate. Choices are: " +
62  InteractionIntegral::sifModeType().getRawNames());
63 
64  return params;
65 }
66 
67 InteractionIntegral::InteractionIntegral(const InputParameters & parameters)
68  : ElementIntegralPostprocessor(parameters),
69  _ndisp(coupledComponents("displacements")),
70  _crack_front_definition(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
71  _has_crack_front_point_index(isParamValid("crack_front_point_index")),
73  _has_crack_front_point_index ? getParam<unsigned int>("crack_front_point_index") : 0),
74  _treat_as_2d(false),
75  _stress(hasMaterialProperty<RankTwoTensor>("stress")
76  ? &getMaterialPropertyByName<RankTwoTensor>("stress")
77  : nullptr),
78  _strain(hasMaterialProperty<RankTwoTensor>("elastic_strain")
79  ? &getMaterialPropertyByName<RankTwoTensor>("elastic_strain")
80  : nullptr),
81  _grad_disp(3),
82  _has_temp(isCoupled("temp")),
83  _grad_temp(_has_temp ? coupledGradient("temp") : _grad_zero),
84  _K_factor(getParam<Real>("K_factor")),
85  _has_symmetry_plane(isParamValid("symmetry_plane")),
86  _poissons_ratio(getParam<Real>("poissons_ratio")),
87  _youngs_modulus(getParam<Real>("youngs_modulus")),
88  _ring_index(getParam<unsigned int>("ring_index")),
89  _total_deigenstrain_dT(hasMaterialProperty<RankTwoTensor>("total_deigenstrain_dT")
90  ? &getMaterialProperty<RankTwoTensor>("total_deigenstrain_dT")
91  : nullptr),
92  _q_function_type(getParam<MooseEnum>("q_function_type").getEnum<QMethod>()),
93  _sif_mode(getParam<MooseEnum>("sif_mode").getEnum<SifMethod>())
94 {
95  if (!hasMaterialProperty<RankTwoTensor>("stress"))
96  mooseError("InteractionIntegral Error: RankTwoTensor material property 'stress' not found. "
97  "This may be because solid mechanics system is being used to calculate a SymmTensor "
98  "'stress' material property. To use interaction integral calculation with solid "
99  "mechanics application, please set 'solid_mechanics = true' in the DomainIntegral "
100  "block.");
101 
102  if (!hasMaterialProperty<RankTwoTensor>("elastic_strain"))
103  mooseError("InteractionIntegral Error: RankTwoTensor material property 'elastic_strain' not "
104  "found. This may be because solid mechanics system is being used to calculate a "
105  "SymmTensor 'elastic_strain' material property. To use interaction integral "
106  "calculation with solid mechanics application, please set 'solid_mechanics = true' "
107  "in the DomainIntegral block.");
108 
110  mooseError("InteractionIntegral Error: To include thermal strain term in interaction integral, "
111  "must both couple temperature in DomainIntegral block and compute "
112  "total_deigenstrain_dT using ThermalFractureIntegral material model.");
113 
114  // plane strain
115  _kappa = 3.0 - 4.0 * _poissons_ratio;
116  _shear_modulus = _youngs_modulus / (2.0 * (1.0 + _poissons_ratio));
117 
118  // Checking for consistency between mesh size and length of the provided displacements vector
119  if (_ndisp != _mesh.dimension())
120  mooseError("InteractionIntegral Error: number of variables supplied in 'displacements' must "
121  "match the mesh dimension.");
122 
123  // fetch gradients of coupled variables
124  for (unsigned int i = 0; i < _ndisp; ++i)
125  _grad_disp[i] = &coupledGradient("displacements", i);
126 
127  // set unused dimensions to zero
128  for (unsigned i = _ndisp; i < 3; ++i)
129  _grad_disp[i] = &_grad_zero;
130 }
131 
132 void
134 {
136 }
137 
138 Real
140 {
141  gatherSum(_integral_value);
142 
144  _integral_value +=
147 
148  return _K_factor * _integral_value;
149 }
150 
151 Real
153 {
154  Real scalar_q = 0.0;
155  RealVectorValue grad_q(0.0, 0.0, 0.0);
156 
157  const std::vector<std::vector<Real>> & phi_curr_elem = *_phi_curr_elem;
158  const std::vector<std::vector<RealGradient>> & dphi_curr_elem = *_dphi_curr_elem;
159 
160  for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i)
161  {
162  scalar_q += phi_curr_elem[i][_qp] * _q_curr_elem[i];
163 
164  for (unsigned int j = 0; j < _current_elem->dim(); ++j)
165  grad_q(j) += dphi_curr_elem[i][_qp](j) * _q_curr_elem[i];
166  }
167 
168  // In the crack front coordinate system, the crack direction is (1,0,0)
169  RealVectorValue crack_direction(0.0);
170  crack_direction(0) = 1.0;
171 
172  // Calculate (r,theta) position of qp relative to crack front
173  Point p(_q_point[_qp]);
175 
176  RankTwoTensor aux_stress;
177  RankTwoTensor aux_du;
178 
180  computeAuxFields(aux_stress, aux_du);
181  else if (_sif_mode == SifMethod::T)
182  computeTFields(aux_stress, aux_du);
183 
184  RankTwoTensor grad_disp((*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
185 
186  // Rotate stress, strain, displacement and temperature to crack front coordinate system
187  RealVectorValue grad_q_cf =
189  RankTwoTensor grad_disp_cf =
191  RankTwoTensor stress_cf =
193  RankTwoTensor strain_cf =
195  RealVectorValue grad_temp_cf =
197 
198  RankTwoTensor dq;
199  dq(0, 0) = crack_direction(0) * grad_q_cf(0);
200  dq(0, 1) = crack_direction(0) * grad_q_cf(1);
201  dq(0, 2) = crack_direction(0) * grad_q_cf(2);
202 
203  // Calculate interaction integral terms
204 
205  // Term1 = stress * x1-derivative of aux disp * dq
206  RankTwoTensor tmp1 = dq * stress_cf;
207  Real term1 = aux_du.doubleContraction(tmp1);
208 
209  // Term2 = aux stress * x1-derivative of disp * dq
210  RankTwoTensor tmp2 = dq * aux_stress;
211  Real term2 = grad_disp_cf(0, 0) * tmp2(0, 0) + grad_disp_cf(1, 0) * tmp2(0, 1) +
212  grad_disp_cf(2, 0) * tmp2(0, 2);
213 
214  // Term3 = aux stress * strain * dq_x (= stress * aux strain * dq_x)
215  Real term3 = dq(0, 0) * aux_stress.doubleContraction(strain_cf);
216 
217  // Term4 (thermal strain term) = q * aux_stress * alpha * dtheta_x
218  // - the term including the derivative of alpha is not implemented
219  Real term4 = 0.0;
220  if (_has_temp)
221  {
222  Real sigma_alpha = aux_stress.doubleContraction((*_total_deigenstrain_dT)[_qp]);
223  term4 = scalar_q * sigma_alpha * grad_temp_cf(0);
224  }
225 
226  Real q_avg_seg = 1.0;
228  {
229  q_avg_seg =
232  2.0;
233  }
234 
235  Real eq = term1 + term2 - term3 + term4;
236 
238  eq *= 2.0;
239 
240  return eq / q_avg_seg;
241 }
242 
243 Real
245 {
246  Real sum = 0.0;
247 
248  // calculate phi and dphi for this element
249  FEType fe_type(Utility::string_to_enum<Order>("first"),
250  Utility::string_to_enum<FEFamily>("lagrange"));
251  const unsigned int dim = _current_elem->dim();
252  UniquePtr<FEBase> fe(FEBase::build(dim, fe_type));
253  fe->attach_quadrature_rule(_qrule);
254  _phi_curr_elem = &fe->get_phi();
255  _dphi_curr_elem = &fe->get_dphi();
256  fe->reinit(_current_elem);
257 
258  // calculate q for all nodes in this element
259  _q_curr_elem.clear();
260  unsigned int ring_base = (_q_function_type == QMethod::Topology) ? 0 : 1;
261 
262  for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i)
263  {
264  Node * this_node = _current_elem->get_node(i);
265  Real q_this_node;
266 
269  _crack_front_point_index, _ring_index - ring_base, this_node);
272  _crack_front_point_index, _ring_index - ring_base, this_node);
273 
274  _q_curr_elem.push_back(q_this_node);
275  }
276 
277  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
278  sum += _JxW[_qp] * _coord[_qp] * computeQpIntegral();
279  return sum;
280 }
281 
282 void
283 InteractionIntegral::computeAuxFields(RankTwoTensor & aux_stress, RankTwoTensor & grad_disp)
284 {
285  RealVectorValue k(0.0);
286  if (_sif_mode == SifMethod::KI)
287  k(0) = 1.0;
288  else if (_sif_mode == SifMethod::KII)
289  k(1) = 1.0;
290  else if (_sif_mode == SifMethod::KIII)
291  k(2) = 1.0;
292 
293  Real t = _theta;
294  Real t2 = _theta / 2.0;
295  Real tt2 = 3.0 * _theta / 2.0;
296  Real st = std::sin(t);
297  Real ct = std::cos(t);
298  Real st2 = std::sin(t2);
299  Real ct2 = std::cos(t2);
300  Real stt2 = std::sin(tt2);
301  Real ctt2 = std::cos(tt2);
302  Real ct2sq = Utility::pow<2>(ct2);
303  Real ct2cu = Utility::pow<3>(ct2);
304  Real sqrt2PiR = std::sqrt(2.0 * libMesh::pi * _r);
305 
306  // Calculate auxiliary stress tensor
307  aux_stress.zero();
308 
309  aux_stress(0, 0) =
310  1.0 / sqrt2PiR * (k(0) * ct2 * (1.0 - st2 * stt2) - k(1) * st2 * (2.0 + ct2 * ctt2));
311  aux_stress(1, 1) = 1.0 / sqrt2PiR * (k(0) * ct2 * (1.0 + st2 * stt2) + k(1) * st2 * ct2 * ctt2);
312  aux_stress(0, 1) = 1.0 / sqrt2PiR * (k(0) * ct2 * st2 * ctt2 + k(1) * ct2 * (1.0 - st2 * stt2));
313  aux_stress(0, 2) = -1.0 / sqrt2PiR * k(2) * st2;
314  aux_stress(1, 2) = 1.0 / sqrt2PiR * k(2) * ct2;
315  // plane stress
316  // Real s33 = 0;
317  // plane strain
318  aux_stress(2, 2) = _poissons_ratio * (aux_stress(0, 0) + aux_stress(1, 1));
319 
320  aux_stress(1, 0) = aux_stress(0, 1);
321  aux_stress(2, 0) = aux_stress(0, 2);
322  aux_stress(2, 1) = aux_stress(1, 2);
323 
324  // Calculate x1 derivative of auxiliary displacements
325  grad_disp.zero();
326 
327  grad_disp(0, 0) = k(0) / (4.0 * _shear_modulus * sqrt2PiR) *
328  (ct * ct2 * _kappa + ct * ct2 - 2.0 * ct * ct2cu + st * st2 * _kappa +
329  st * st2 - 6.0 * st * st2 * ct2sq) +
330  k(1) / (4.0 * _shear_modulus * sqrt2PiR) *
331  (ct * st2 * _kappa + ct * st2 + 2.0 * ct * st2 * ct2sq - st * ct2 * _kappa +
332  3.0 * st * ct2 - 6.0 * st * ct2cu);
333 
334  grad_disp(0, 1) = k(0) / (4.0 * _shear_modulus * sqrt2PiR) *
335  (ct * st2 * _kappa + ct * st2 - 2.0 * ct * st2 * ct2sq - st * ct2 * _kappa -
336  5.0 * st * ct2 + 6.0 * st * ct2cu) +
337  k(1) / (4.0 * _shear_modulus * sqrt2PiR) *
338  (-ct * ct2 * _kappa + 3.0 * ct * ct2 - 2.0 * ct * ct2cu -
339  st * st2 * _kappa + 3.0 * st * st2 - 6.0 * st * st2 * ct2sq);
340 
341  grad_disp(0, 2) = k(2) / (_shear_modulus * sqrt2PiR) * (st2 * ct - ct2 * st);
342 }
343 
344 void
345 InteractionIntegral::computeTFields(RankTwoTensor & aux_stress, RankTwoTensor & grad_disp)
346 {
347 
348  Real t = _theta;
349  Real st = std::sin(t);
350  Real ct = std::cos(t);
351  Real stsq = Utility::pow<2>(st);
352  Real ctsq = Utility::pow<2>(ct);
353  Real ctcu = Utility::pow<3>(ct);
354  Real oneOverPiR = 1.0 / (libMesh::pi * _r);
355 
356  aux_stress.zero();
357  aux_stress(0, 0) = -oneOverPiR * ctcu;
358  aux_stress(0, 1) = -oneOverPiR * st * ctsq;
359  aux_stress(1, 0) = -oneOverPiR * st * ctsq;
360  aux_stress(1, 1) = -oneOverPiR * ct * stsq;
361  aux_stress(2, 2) = -oneOverPiR * _poissons_ratio * (ctcu + ct * stsq);
362 
363  grad_disp.zero();
364  grad_disp(0, 0) = oneOverPiR / (4.0 * _youngs_modulus) *
365  (ct * (4.0 * Utility::pow<2>(_poissons_ratio) - 3.0 + _poissons_ratio) -
366  std::cos(3.0 * t) * (1.0 + _poissons_ratio));
367  grad_disp(0, 1) = -oneOverPiR / (4.0 * _youngs_modulus) *
368  (st * (4.0 * Utility::pow<2>(_poissons_ratio) - 3.0 + _poissons_ratio) +
369  std::sin(3.0 * t) * (1.0 + _poissons_ratio));
370 }
void computeAuxFields(RankTwoTensor &aux_stress, RankTwoTensor &grad_disp)
static MooseEnum qFunctionType()
virtual Real computeIntegral()
Real getCrackFrontTangentialStrain(const unsigned int node_index) const
void calculateRThetaToCrackFront(const Point qp, const unsigned int point_index, Real &r, Real &theta) const
calculate r and theta in the crack front polar cooridnate
InputParameters validParams< InteractionIntegral >()
const std::vector< std::vector< Real > > * _phi_curr_elem
const unsigned int _crack_front_point_index
const QMethod _q_function_type
const MaterialProperty< RankTwoTensor > * _stress
const CrackFrontDefinition *const _crack_front_definition
std::vector< const VariableGradient * > _grad_disp
Works on top of NodalNormalsPreprocessor.
void computeTFields(RankTwoTensor &aux_stress, RankTwoTensor &grad_disp)
Real DomainIntegralTopologicalQFunction(unsigned int crack_front_point_index, unsigned int ring_index, const Node *const current_node) const
const std::vector< std::vector< RealGradient > > * _dphi_curr_elem
Real DomainIntegralQFunction(unsigned int crack_front_point_index, unsigned int ring_index, const Node *const current_node) const
virtual Real computeQpIntegral()
Real getCrackFrontBackwardSegmentLength(const unsigned int point_index) const
virtual void initialSetup()
RealVectorValue rotateToCrackFrontCoords(const RealVectorValue vector, const unsigned int point_index) const
const MaterialProperty< RankTwoTensor > * _strain
const SifMethod _sif_mode
const MaterialProperty< RankTwoTensor > * _total_deigenstrain_dT
static MooseEnum sifModeType()
Real getCrackFrontForwardSegmentLength(const unsigned int point_index) const
const VariableGradient & _grad_temp
std::vector< Real > _q_curr_elem
InteractionIntegral(const InputParameters &parameters)