Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "JIntegral.h"
11 : #include "RankTwoTensor.h"
12 : #include "Conversion.h"
13 : #include "libmesh/string_to_enum.h"
14 : #include "libmesh/quadrature.h"
15 : #include "libmesh/utility.h"
16 : #include "CrackFrontDefinition.h"
17 : registerMooseObject("TensorMechanicsApp", JIntegral);
18 :
19 : InputParameters
20 547 : JIntegral::validParams()
21 : {
22 547 : InputParameters params = ElementVectorPostprocessor::validParams();
23 1094 : params.addRequiredCoupledVar(
24 : "displacements",
25 : "The displacements appropriate for the simulation geometry and coordinate system");
26 1094 : params.addRequiredParam<UserObjectName>("crack_front_definition",
27 : "The CrackFrontDefinition user object name");
28 1094 : MooseEnum position_type("Angle Distance", "Distance");
29 1094 : params.addParam<MooseEnum>(
30 : "position_type",
31 : position_type,
32 547 : "The method used to calculate position along crack front. Options are: " +
33 547 : position_type.getRawNames());
34 1094 : MooseEnum integral_vec("JIntegral CIntegral KFromJIntegral", "JIntegral");
35 1094 : params.addRequiredParam<MooseEnum>("integral",
36 : integral_vec,
37 547 : "Domain integrals to calculate. Choices are: " +
38 547 : integral_vec.getRawNames());
39 1094 : params.addParam<unsigned int>("symmetry_plane",
40 : "Account for a symmetry plane passing through "
41 : "the plane of the crack, normal to the specified "
42 : "axis (0=x, 1=y, 2=z)");
43 1094 : params.addParam<Real>("poissons_ratio", "Poisson's ratio");
44 1094 : params.addParam<Real>("youngs_modulus", "Young's modulus of the material.");
45 547 : params.set<bool>("use_displaced_mesh") = false;
46 1094 : params.addRequiredParam<unsigned int>("ring_index", "Ring ID");
47 1094 : MooseEnum q_function_type("Geometry Topology", "Geometry");
48 1094 : params.addParam<MooseEnum>("q_function_type",
49 : q_function_type,
50 547 : "The method used to define the integration domain. Options are: " +
51 547 : q_function_type.getRawNames());
52 547 : params.addClassDescription("Computes the J-Integral, a measure of the strain energy release rate "
53 : "at a crack tip, which can be used as a criterion for fracture "
54 : "growth. It can, alternatively, compute the C(t) integral");
55 547 : return params;
56 547 : }
57 :
58 428 : JIntegral::JIntegral(const InputParameters & parameters)
59 : : ElementVectorPostprocessor(parameters),
60 428 : _crack_front_definition(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
61 856 : _integral(getParam<MooseEnum>("integral").getEnum<IntegralType>()),
62 428 : _J_thermal_term_vec(hasMaterialProperty<RealVectorValue>("J_thermal_term_vec")
63 1284 : ? &getMaterialProperty<RealVectorValue>("J_thermal_term_vec")
64 : : nullptr),
65 856 : _Eshelby_tensor(_integral != IntegralType::CIntegral
66 832 : ? &getMaterialProperty<RankTwoTensor>("Eshelby_tensor")
67 : : nullptr),
68 428 : _Eshelby_tensor_dissipation(
69 428 : _integral == IntegralType::CIntegral
70 452 : ? &getMaterialProperty<RankTwoTensor>("Eshelby_tensor_dissipation")
71 : : nullptr),
72 428 : _fe_vars(getCoupledMooseVars()),
73 428 : _fe_type(_fe_vars[0]->feType()),
74 856 : _has_symmetry_plane(isParamValid("symmetry_plane")),
75 880 : _poissons_ratio(isParamValid("poissons_ratio") ? getParam<Real>("poissons_ratio") : 0),
76 880 : _youngs_modulus(isParamValid("youngs_modulus") ? getParam<Real>("youngs_modulus") : 0),
77 856 : _ring_index(getParam<unsigned int>("ring_index")),
78 856 : _q_function_type(getParam<MooseEnum>("q_function_type").getEnum<QMethod>()),
79 856 : _position_type(getParam<MooseEnum>("position_type").getEnum<PositionType>()),
80 428 : _x(declareVector("x")),
81 428 : _y(declareVector("y")),
82 428 : _z(declareVector("z")),
83 428 : _position(declareVector("id")),
84 428 : _j_integral(declareVector((_integral == IntegralType::KFromJIntegral
85 : ? "K_"
86 452 : : (_integral == IntegralType::JIntegral ? "J_" : "C_")) +
87 856 : Moose::stringify(_ring_index)))
88 : {
89 428 : }
90 :
91 : void
92 423 : JIntegral::initialSetup()
93 : {
94 423 : _treat_as_2d = _crack_front_definition->treatAs2D();
95 423 : _using_mesh_cutter = _crack_front_definition->usingMeshCutter();
96 435 : if (_integral == IntegralType::KFromJIntegral &&
97 471 : (!isParamValid("youngs_modulus") || !isParamValid("poissons_ratio")))
98 0 : mooseError("youngs_modulus and poissons_ratio must be specified if integrals = KFromJIntegral");
99 423 : }
100 :
101 : void
102 423 : JIntegral::initialize()
103 : {
104 : std::size_t num_pts;
105 423 : if (_treat_as_2d && _using_mesh_cutter == false)
106 : num_pts = 1;
107 : else
108 117 : num_pts = _crack_front_definition->getNumCrackFrontPoints();
109 :
110 423 : _x.assign(num_pts, 0.0);
111 423 : _y.assign(num_pts, 0.0);
112 423 : _z.assign(num_pts, 0.0);
113 423 : _position.assign(num_pts, 0.0);
114 423 : _j_integral.assign(num_pts, 0.0);
115 :
116 1419 : for (const auto * fe_var : _fe_vars)
117 : {
118 996 : if (fe_var->feType() != _fe_type)
119 0 : mooseError("displacements", "All coupled displacements must have the same type");
120 : }
121 423 : }
122 :
123 : Real
124 1174248 : JIntegral::computeQpIntegral(const std::size_t crack_front_point_index,
125 : const Real scalar_q,
126 : const RealVectorValue & grad_of_scalar_q)
127 : {
128 1174248 : RankTwoTensor grad_of_vector_q;
129 : const RealVectorValue & crack_direction =
130 1174248 : _crack_front_definition->getCrackDirection(crack_front_point_index);
131 1174248 : grad_of_vector_q(0, 0) = crack_direction(0) * grad_of_scalar_q(0);
132 1174248 : grad_of_vector_q(0, 1) = crack_direction(0) * grad_of_scalar_q(1);
133 1174248 : grad_of_vector_q(0, 2) = crack_direction(0) * grad_of_scalar_q(2);
134 1174248 : grad_of_vector_q(1, 0) = crack_direction(1) * grad_of_scalar_q(0);
135 1174248 : grad_of_vector_q(1, 1) = crack_direction(1) * grad_of_scalar_q(1);
136 1174248 : grad_of_vector_q(1, 2) = crack_direction(1) * grad_of_scalar_q(2);
137 1174248 : grad_of_vector_q(2, 0) = crack_direction(2) * grad_of_scalar_q(0);
138 1174248 : grad_of_vector_q(2, 1) = crack_direction(2) * grad_of_scalar_q(1);
139 1174248 : grad_of_vector_q(2, 2) = crack_direction(2) * grad_of_scalar_q(2);
140 :
141 : Real eq;
142 :
143 1174248 : if (_integral != IntegralType::CIntegral)
144 1161448 : eq = (*_Eshelby_tensor)[_qp].doubleContraction(grad_of_vector_q);
145 : else
146 12800 : eq = (*_Eshelby_tensor_dissipation)[_qp].doubleContraction(grad_of_vector_q);
147 :
148 : // Thermal component
149 : Real eq_thermal = 0.0;
150 1174248 : if (_J_thermal_term_vec && _integral != IntegralType::CIntegral)
151 : {
152 4645792 : for (std::size_t i = 0; i < 3; i++)
153 3484344 : eq_thermal += crack_direction(i) * scalar_q * (*_J_thermal_term_vec)[_qp](i);
154 : }
155 :
156 : Real q_avg_seg = 1.0;
157 1174248 : if (!_crack_front_definition->treatAs2D())
158 : {
159 949200 : q_avg_seg =
160 1898400 : (_crack_front_definition->getCrackFrontForwardSegmentLength(crack_front_point_index) +
161 949200 : _crack_front_definition->getCrackFrontBackwardSegmentLength(crack_front_point_index)) /
162 : 2.0;
163 : }
164 :
165 1174248 : Real etot = -eq + eq_thermal;
166 :
167 1174248 : return etot / q_avg_seg;
168 : }
169 :
170 : void
171 73964 : JIntegral::execute()
172 : {
173 : // calculate phi and dphi for this element
174 73964 : const std::size_t dim = _current_elem->dim();
175 73964 : std::unique_ptr<FEBase> fe(FEBase::build(dim, _fe_type));
176 73964 : fe->attach_quadrature_rule(const_cast<QBase *>(_qrule));
177 73964 : _phi_curr_elem = &fe->get_phi();
178 73964 : _dphi_curr_elem = &fe->get_dphi();
179 73964 : fe->reinit(_current_elem);
180 :
181 : // calculate q for all nodes in this element
182 73964 : std::size_t ring_base = (_q_function_type == QMethod::Topology) ? 0 : 1;
183 :
184 248084 : for (std::size_t icfp = 0; icfp < _j_integral.size(); icfp++)
185 : {
186 174120 : _q_curr_elem.clear();
187 1348368 : for (std::size_t i = 0; i < _current_elem->n_nodes(); ++i)
188 : {
189 1174248 : const Node * this_node = _current_elem->node_ptr(i);
190 : Real q_this_node;
191 :
192 1174248 : if (_q_function_type == QMethod::Geometry)
193 1149288 : q_this_node = _crack_front_definition->DomainIntegralQFunction(
194 1149288 : icfp, _ring_index - ring_base, this_node);
195 : else // QMethod::Topology
196 24960 : q_this_node = _crack_front_definition->DomainIntegralTopologicalQFunction(
197 24960 : icfp, _ring_index - ring_base, this_node);
198 :
199 1174248 : _q_curr_elem.push_back(q_this_node);
200 : }
201 :
202 1348368 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
203 : {
204 : Real scalar_q = 0.0;
205 : RealVectorValue grad_of_scalar_q(0.0, 0.0, 0.0);
206 :
207 9693384 : for (std::size_t i = 0; i < _current_elem->n_nodes(); ++i)
208 : {
209 8519136 : scalar_q += (*_phi_curr_elem)[i][_qp] * _q_curr_elem[i];
210 :
211 33201696 : for (std::size_t j = 0; j < _current_elem->dim(); ++j)
212 24682560 : grad_of_scalar_q(j) += (*_dphi_curr_elem)[i][_qp](j) * _q_curr_elem[i];
213 : }
214 :
215 1174248 : _j_integral[icfp] +=
216 1174248 : _JxW[_qp] * _coord[_qp] * computeQpIntegral(icfp, scalar_q, grad_of_scalar_q);
217 : }
218 : }
219 73964 : }
220 :
221 : void
222 423 : JIntegral::finalize()
223 : {
224 423 : gatherSum(_j_integral);
225 :
226 1272 : for (std::size_t i = 0; i < _j_integral.size(); ++i)
227 : {
228 849 : if (_has_symmetry_plane)
229 738 : _j_integral[i] *= 2.0;
230 :
231 849 : Real sign = (_j_integral[i] > 0.0) ? 1.0 : ((_j_integral[i] < 0.0) ? -1.0 : 0.0);
232 :
233 849 : if (_integral == IntegralType::KFromJIntegral)
234 12 : _j_integral[i] = sign * std::sqrt(std::abs(_j_integral[i]) * _youngs_modulus /
235 12 : (1.0 - Utility::pow<2>(_poissons_ratio)));
236 :
237 849 : const auto cfp = _crack_front_definition->getCrackFrontPoint(i);
238 849 : _x[i] = (*cfp)(0);
239 849 : _y[i] = (*cfp)(1);
240 849 : _z[i] = (*cfp)(2);
241 :
242 849 : if (_position_type == PositionType::Angle)
243 126 : _position[i] = _crack_front_definition->getAngleAlongFront(i);
244 : else
245 723 : _position[i] = _crack_front_definition->getDistanceAlongFront(i);
246 : }
247 423 : }
248 :
249 : void
250 0 : JIntegral::threadJoin(const UserObject & y)
251 : {
252 : const auto & uo = static_cast<const JIntegral &>(y);
253 :
254 0 : for (auto i = beginIndex(_j_integral); i < _j_integral.size(); ++i)
255 0 : _j_integral[i] += uo._j_integral[i];
256 0 : }
|
