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