www.mooseframework.org
Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Types | Private Attributes | List of all members
InteractionIntegral Class Reference

This postprocessor computes the Interaction Integral. More...

#include <InteractionIntegral.h>

Inheritance diagram for InteractionIntegral:
[legend]

Public Member Functions

 InteractionIntegral (const InputParameters &parameters)
 
virtual Real getValue ()
 

Static Public Member Functions

static MooseEnum qFunctionType ()
 
static MooseEnum sifModeType ()
 

Protected Member Functions

virtual void initialSetup ()
 
virtual Real computeQpIntegral ()
 
virtual Real computeIntegral ()
 
void computeAuxFields (RankTwoTensor &aux_stress, RankTwoTensor &grad_disp)
 
void computeTFields (RankTwoTensor &aux_stress, RankTwoTensor &grad_disp)
 

Protected Attributes

unsigned int _ndisp
 
const CrackFrontDefinition *const _crack_front_definition
 
bool _has_crack_front_point_index
 
const unsigned int _crack_front_point_index
 
bool _treat_as_2d
 
const MaterialProperty< RankTwoTensor > * _stress
 
const MaterialProperty< RankTwoTensor > * _strain
 
std::vector< const VariableGradient * > _grad_disp
 
const bool _has_temp
 
const VariableGradient & _grad_temp
 
Real _K_factor
 
bool _has_symmetry_plane
 
Real _poissons_ratio
 
Real _youngs_modulus
 
unsigned int _ring_index
 
const MaterialProperty< RankTwoTensor > * _total_deigenstrain_dT
 
std::vector< Real > _q_curr_elem
 
const std::vector< std::vector< Real > > * _phi_curr_elem
 
const std::vector< std::vector< RealGradient > > * _dphi_curr_elem
 
Real _kappa
 
Real _shear_modulus
 
Real _r
 
Real _theta
 

Private Types

enum  QMethod { QMethod::Geometry, QMethod::Topology }
 
enum  SifMethod { SifMethod::KI, SifMethod::KII, SifMethod::KIII, SifMethod::T }
 

Private Attributes

const QMethod _q_function_type
 
const SifMethod _sif_mode
 

Detailed Description

This postprocessor computes the Interaction Integral.

Definition at line 24 of file InteractionIntegral.h.

Member Enumeration Documentation

enum InteractionIntegral::QMethod
strongprivate
Enumerator
Geometry 
Topology 

Definition at line 65 of file InteractionIntegral.h.

66  {
67  Geometry,
68  Topology
69  };
enum InteractionIntegral::SifMethod
strongprivate
Enumerator
KI 
KII 
KIII 

Definition at line 73 of file InteractionIntegral.h.

74  {
75  KI,
76  KII,
77  KIII,
78  T
79  };

Constructor & Destructor Documentation

InteractionIntegral::InteractionIntegral ( const InputParameters &  parameters)

Definition at line 67 of file InteractionIntegral.C.

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 }
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
const MaterialProperty< RankTwoTensor > * _strain
const SifMethod _sif_mode
const MaterialProperty< RankTwoTensor > * _total_deigenstrain_dT
const VariableGradient & _grad_temp

Member Function Documentation

void InteractionIntegral::computeAuxFields ( RankTwoTensor &  aux_stress,
RankTwoTensor &  grad_disp 
)
protected

Definition at line 283 of file InteractionIntegral.C.

Referenced by computeQpIntegral().

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 }
const SifMethod _sif_mode
Real InteractionIntegral::computeIntegral ( )
protectedvirtual

Definition at line 244 of file InteractionIntegral.C.

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 }
const std::vector< std::vector< Real > > * _phi_curr_elem
const unsigned int _crack_front_point_index
const QMethod _q_function_type
const CrackFrontDefinition *const _crack_front_definition
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()
std::vector< Real > _q_curr_elem
Real InteractionIntegral::computeQpIntegral ( )
protectedvirtual

Definition at line 152 of file InteractionIntegral.C.

Referenced by computeIntegral().

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 }
void computeAuxFields(RankTwoTensor &aux_stress, RankTwoTensor &grad_disp)
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
const std::vector< std::vector< Real > > * _phi_curr_elem
const unsigned int _crack_front_point_index
const MaterialProperty< RankTwoTensor > * _stress
const CrackFrontDefinition *const _crack_front_definition
std::vector< const VariableGradient * > _grad_disp
void computeTFields(RankTwoTensor &aux_stress, RankTwoTensor &grad_disp)
const std::vector< std::vector< RealGradient > > * _dphi_curr_elem
Real getCrackFrontBackwardSegmentLength(const unsigned int point_index) const
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
Real getCrackFrontForwardSegmentLength(const unsigned int point_index) const
const VariableGradient & _grad_temp
std::vector< Real > _q_curr_elem
void InteractionIntegral::computeTFields ( RankTwoTensor &  aux_stress,
RankTwoTensor &  grad_disp 
)
protected

Definition at line 345 of file InteractionIntegral.C.

Referenced by computeQpIntegral().

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 }
Real InteractionIntegral::getValue ( )
virtual

Definition at line 139 of file InteractionIntegral.C.

140 {
141  gatherSum(_integral_value);
142 
144  _integral_value +=
147 
148  return _K_factor * _integral_value;
149 }
Real getCrackFrontTangentialStrain(const unsigned int node_index) const
const unsigned int _crack_front_point_index
const CrackFrontDefinition *const _crack_front_definition
const SifMethod _sif_mode
void InteractionIntegral::initialSetup ( )
protectedvirtual

Definition at line 133 of file InteractionIntegral.C.

134 {
136 }
const CrackFrontDefinition *const _crack_front_definition
MooseEnum InteractionIntegral::qFunctionType ( )
static

Definition at line 18 of file InteractionIntegral.C.

Referenced by validParams< InteractionIntegral >().

19 {
20  return MooseEnum("Geometry Topology", "Geometry");
21 }
MooseEnum InteractionIntegral::sifModeType ( )
static

Definition at line 24 of file InteractionIntegral.C.

Referenced by validParams< InteractionIntegral >().

25 {
26  return MooseEnum("KI KII KIII T", "KI");
27 }

Member Data Documentation

const CrackFrontDefinition* const InteractionIntegral::_crack_front_definition
protected

Definition at line 41 of file InteractionIntegral.h.

Referenced by computeIntegral(), computeQpIntegral(), getValue(), and initialSetup().

const unsigned int InteractionIntegral::_crack_front_point_index
protected

Definition at line 43 of file InteractionIntegral.h.

Referenced by computeIntegral(), computeQpIntegral(), and getValue().

const std::vector<std::vector<RealGradient> >* InteractionIntegral::_dphi_curr_elem
protected

Definition at line 58 of file InteractionIntegral.h.

Referenced by computeIntegral(), and computeQpIntegral().

std::vector<const VariableGradient *> InteractionIntegral::_grad_disp
protected

Definition at line 47 of file InteractionIntegral.h.

Referenced by computeQpIntegral(), and InteractionIntegral().

const VariableGradient& InteractionIntegral::_grad_temp
protected

Definition at line 49 of file InteractionIntegral.h.

Referenced by computeQpIntegral().

bool InteractionIntegral::_has_crack_front_point_index
protected

Definition at line 42 of file InteractionIntegral.h.

bool InteractionIntegral::_has_symmetry_plane
protected

Definition at line 51 of file InteractionIntegral.h.

Referenced by computeQpIntegral().

const bool InteractionIntegral::_has_temp
protected

Definition at line 48 of file InteractionIntegral.h.

Referenced by computeQpIntegral(), and InteractionIntegral().

Real InteractionIntegral::_K_factor
protected

Definition at line 50 of file InteractionIntegral.h.

Referenced by getValue().

Real InteractionIntegral::_kappa
protected

Definition at line 59 of file InteractionIntegral.h.

Referenced by computeAuxFields(), and InteractionIntegral().

unsigned int InteractionIntegral::_ndisp
protected

Definition at line 40 of file InteractionIntegral.h.

Referenced by InteractionIntegral().

const std::vector<std::vector<Real> >* InteractionIntegral::_phi_curr_elem
protected

Definition at line 57 of file InteractionIntegral.h.

Referenced by computeIntegral(), and computeQpIntegral().

Real InteractionIntegral::_poissons_ratio
protected
std::vector<Real> InteractionIntegral::_q_curr_elem
protected

Definition at line 56 of file InteractionIntegral.h.

Referenced by computeIntegral(), and computeQpIntegral().

const QMethod InteractionIntegral::_q_function_type
private

Definition at line 71 of file InteractionIntegral.h.

Referenced by computeIntegral().

Real InteractionIntegral::_r
protected

Definition at line 61 of file InteractionIntegral.h.

Referenced by computeAuxFields(), computeQpIntegral(), and computeTFields().

unsigned int InteractionIntegral::_ring_index
protected

Definition at line 54 of file InteractionIntegral.h.

Referenced by computeIntegral().

Real InteractionIntegral::_shear_modulus
protected

Definition at line 60 of file InteractionIntegral.h.

Referenced by computeAuxFields(), and InteractionIntegral().

const SifMethod InteractionIntegral::_sif_mode
private

Definition at line 81 of file InteractionIntegral.h.

Referenced by computeAuxFields(), computeQpIntegral(), and getValue().

const MaterialProperty<RankTwoTensor>* InteractionIntegral::_strain
protected

Definition at line 46 of file InteractionIntegral.h.

Referenced by computeQpIntegral().

const MaterialProperty<RankTwoTensor>* InteractionIntegral::_stress
protected

Definition at line 45 of file InteractionIntegral.h.

Referenced by computeQpIntegral().

Real InteractionIntegral::_theta
protected

Definition at line 62 of file InteractionIntegral.h.

Referenced by computeAuxFields(), computeQpIntegral(), and computeTFields().

const MaterialProperty<RankTwoTensor>* InteractionIntegral::_total_deigenstrain_dT
protected

Definition at line 55 of file InteractionIntegral.h.

Referenced by computeQpIntegral(), and InteractionIntegral().

bool InteractionIntegral::_treat_as_2d
protected

Definition at line 44 of file InteractionIntegral.h.

Referenced by getValue(), and initialSetup().

Real InteractionIntegral::_youngs_modulus
protected

Definition at line 53 of file InteractionIntegral.h.

Referenced by computeTFields(), and InteractionIntegral().


The documentation for this class was generated from the following files: