www.mooseframework.org
JIntegral.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 calculates the J-Integral
8 //
9 #include "JIntegral.h"
10 #include "RankTwoTensor.h"
11 #include "libmesh/string_to_enum.h"
12 #include "libmesh/quadrature.h"
13 
14 template <>
15 InputParameters
17 {
18  InputParameters params = validParams<ElementIntegralPostprocessor>();
19  params.addRequiredParam<UserObjectName>("crack_front_definition",
20  "The CrackFrontDefinition user object name");
21  params.addParam<unsigned int>(
22  "crack_front_point_index",
23  "The index of the point on the crack front corresponding to this q function");
24  params.addParam<bool>(
25  "convert_J_to_K", false, "Convert J-integral to stress intensity factor K.");
26  params.addParam<unsigned int>("symmetry_plane",
27  "Account for a symmetry plane passing through "
28  "the plane of the crack, normal to the specified "
29  "axis (0=x, 1=y, 2=z)");
30  params.addParam<Real>("poissons_ratio", "Poisson's ratio");
31  params.addParam<Real>("youngs_modulus", "Young's modulus of the material.");
32  params.set<bool>("use_displaced_mesh") = false;
33  params.addParam<unsigned int>("ring_index", "Ring ID");
34  params.addParam<unsigned int>("ring_first", "First Ring ID");
35  MooseEnum q_function_type("Geometry Topology", "Geometry");
36  params.addParam<MooseEnum>("q_function_type",
37  q_function_type,
38  "The method used to define the integration domain. Options are: " +
39  q_function_type.getRawNames());
40  return params;
41 }
42 
43 JIntegral::JIntegral(const InputParameters & parameters)
44  : ElementIntegralPostprocessor(parameters),
45  _crack_front_definition(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
46  _has_crack_front_point_index(isParamValid("crack_front_point_index")),
47  _crack_front_point_index(
48  _has_crack_front_point_index ? getParam<unsigned int>("crack_front_point_index") : 0),
49  _treat_as_2d(false),
50  _Eshelby_tensor(getMaterialProperty<RankTwoTensor>("Eshelby_tensor")),
51  _J_thermal_term_vec(hasMaterialProperty<RealVectorValue>("J_thermal_term_vec")
52  ? &getMaterialProperty<RealVectorValue>("J_thermal_term_vec")
53  : NULL),
54  _convert_J_to_K(getParam<bool>("convert_J_to_K")),
55  _has_symmetry_plane(isParamValid("symmetry_plane")),
56  _poissons_ratio(isParamValid("poissons_ratio") ? getParam<Real>("poissons_ratio") : 0),
57  _youngs_modulus(isParamValid("youngs_modulus") ? getParam<Real>("youngs_modulus") : 0),
58  _ring_index(getParam<unsigned int>("ring_index")),
59  _q_function_type(getParam<MooseEnum>("q_function_type"))
60 {
61  if (_q_function_type == "TOPOLOGY")
62  _ring_first = getParam<unsigned int>("ring_first");
63 }
64 
65 void
67 {
69 
70  if (_convert_J_to_K && (!isParamValid("youngs_modulus") || !isParamValid("poissons_ratio")))
71  mooseError("youngs_modulus and poissons_ratio must be specified if convert_J_to_K = true");
72 }
73 
74 Real
76 {
77  Real scalar_q = 0.0;
78  RealVectorValue grad_of_scalar_q(0.0, 0.0, 0.0);
79 
80  const std::vector<std::vector<Real>> & phi_curr_elem = *_phi_curr_elem;
81  const std::vector<std::vector<RealGradient>> & dphi_curr_elem = *_dphi_curr_elem;
82 
83  for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i)
84  {
85  scalar_q += phi_curr_elem[i][_qp] * _q_curr_elem[i];
86 
87  for (unsigned int j = 0; j < _current_elem->dim(); ++j)
88  grad_of_scalar_q(j) += dphi_curr_elem[i][_qp](j) * _q_curr_elem[i];
89  }
90 
91  RankTwoTensor grad_of_vector_q;
92  const RealVectorValue & crack_direction =
94  grad_of_vector_q(0, 0) = crack_direction(0) * grad_of_scalar_q(0);
95  grad_of_vector_q(0, 1) = crack_direction(0) * grad_of_scalar_q(1);
96  grad_of_vector_q(0, 2) = crack_direction(0) * grad_of_scalar_q(2);
97  grad_of_vector_q(1, 0) = crack_direction(1) * grad_of_scalar_q(0);
98  grad_of_vector_q(1, 1) = crack_direction(1) * grad_of_scalar_q(1);
99  grad_of_vector_q(1, 2) = crack_direction(1) * grad_of_scalar_q(2);
100  grad_of_vector_q(2, 0) = crack_direction(2) * grad_of_scalar_q(0);
101  grad_of_vector_q(2, 1) = crack_direction(2) * grad_of_scalar_q(1);
102  grad_of_vector_q(2, 2) = crack_direction(2) * grad_of_scalar_q(2);
103 
104  Real eq = _Eshelby_tensor[_qp].doubleContraction(grad_of_vector_q);
105 
106  // Thermal component
107  Real eq_thermal = 0.0;
109  {
110  for (unsigned int i = 0; i < 3; i++)
111  eq_thermal += crack_direction(i) * scalar_q * (*_J_thermal_term_vec)[_qp](i);
112  }
113 
114  Real q_avg_seg = 1.0;
116  {
117  q_avg_seg =
120  2.0;
121  }
122 
123  Real etot = -eq + eq_thermal;
124 
125  return etot / q_avg_seg;
126 }
127 
128 Real
130 {
131  Real sum = 0;
132 
133  // calculate phi and dphi for this element
134  FEType fe_type(Utility::string_to_enum<Order>("first"),
135  Utility::string_to_enum<FEFamily>("lagrange"));
136  const unsigned int dim = _current_elem->dim();
137  UniquePtr<FEBase> fe(FEBase::build(dim, fe_type));
138  fe->attach_quadrature_rule(_qrule);
139  _phi_curr_elem = &fe->get_phi();
140  _dphi_curr_elem = &fe->get_dphi();
141  fe->reinit(_current_elem);
142 
143  // calculate q for all nodes in this element
144  _q_curr_elem.clear();
145  unsigned int ring_base = (_q_function_type == "TOPOLOGY") ? 0 : 1;
146 
147  for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i)
148  {
149  Node * this_node = _current_elem->get_node(i);
150  Real q_this_node;
151 
152  if (_q_function_type == "GEOMETRY")
154  _crack_front_point_index, _ring_index - ring_base, this_node);
155  else if (_q_function_type == "TOPOLOGY")
157  _crack_front_point_index, _ring_index - ring_base, this_node);
158 
159  _q_curr_elem.push_back(q_this_node);
160  }
161 
162  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
163  sum += _JxW[_qp] * _coord[_qp] * computeQpIntegral();
164  return sum;
165 }
166 
167 Real
169 {
170  gatherSum(_integral_value);
172  _integral_value *= 2.0;
173 
174  Real sign = (_integral_value > 0.0) ? 1.0 : ((_integral_value < 0.0) ? -1.0 : 0.0);
175  if (_convert_J_to_K)
176  _integral_value = sign * std::sqrt(std::abs(_integral_value) * _youngs_modulus /
177  (1 - std::pow(_poissons_ratio, 2)));
178 
179  return _integral_value;
180 }
virtual Real getValue()
Definition: JIntegral.C:168
virtual void initialSetup()
Definition: JIntegral.C:66
Real sign(Real x)
Definition: MathUtils.h:24
Real _poissons_ratio
Definition: JIntegral.h:42
virtual Real computeIntegral()
Definition: JIntegral.C:129
InputParameters validParams< JIntegral >()
Definition: JIntegral.C:16
JIntegral(const InputParameters &parameters)
Definition: JIntegral.C:43
Real _youngs_modulus
Definition: JIntegral.h:43
const RealVectorValue & getCrackDirection(const unsigned int point_index) const
virtual Real computeQpIntegral()
Definition: JIntegral.C:75
Works on top of NodalNormalsPreprocessor.
const std::vector< std::vector< Real > > * _phi_curr_elem
Definition: JIntegral.h:48
unsigned int _ring_index
Definition: JIntegral.h:44
Real DomainIntegralTopologicalQFunction(unsigned int crack_front_point_index, unsigned int ring_index, const Node *const current_node) const
unsigned int _ring_first
Definition: JIntegral.h:45
const std::vector< std::vector< RealGradient > > * _dphi_curr_elem
Definition: JIntegral.h:49
const CrackFrontDefinition *const _crack_front_definition
Definition: JIntegral.h:34
const unsigned int _crack_front_point_index
Definition: JIntegral.h:36
Real DomainIntegralQFunction(unsigned int crack_front_point_index, unsigned int ring_index, const Node *const current_node) const
Real getCrackFrontBackwardSegmentLength(const unsigned int point_index) const
bool _convert_J_to_K
Definition: JIntegral.h:40
std::vector< Real > _q_curr_elem
Definition: JIntegral.h:47
const MaterialProperty< RankTwoTensor > & _Eshelby_tensor
Definition: JIntegral.h:38
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
bool _has_symmetry_plane
Definition: JIntegral.h:41
MooseEnum _q_function_type
Definition: JIntegral.h:46
Real getCrackFrontForwardSegmentLength(const unsigned int point_index) const
const MaterialProperty< RealVectorValue > * _J_thermal_term_vec
Definition: JIntegral.h:39
bool _treat_as_2d
Definition: JIntegral.h:37