Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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("SolidMechanicsApp", JIntegral);
18 :
19 : InputParameters
20 680 : JIntegral::validParams()
21 : {
22 680 : InputParameters params = ElementVectorPostprocessor::validParams();
23 1360 : params.addRequiredCoupledVar(
24 : "displacements",
25 : "The displacements appropriate for the simulation geometry and coordinate system");
26 1360 : params.addRequiredParam<UserObjectName>("crack_front_definition",
27 : "The CrackFrontDefinition user object name");
28 1360 : MooseEnum position_type("Angle Distance", "Distance");
29 1360 : params.addParam<MooseEnum>(
30 : "position_type",
31 : position_type,
32 680 : "The method used to calculate position along crack front. Options are: " +
33 680 : position_type.getRawNames());
34 1360 : MooseEnum integral_vec("JIntegral CIntegral KFromJIntegral", "JIntegral");
35 1360 : params.addRequiredParam<MooseEnum>("integral",
36 : integral_vec,
37 680 : "Domain integrals to calculate. Choices are: " +
38 680 : integral_vec.getRawNames());
39 1360 : 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 1360 : params.addParam<Real>("poissons_ratio", "Poisson's ratio");
44 1360 : params.addParam<Real>("youngs_modulus", "Young's modulus of the material.");
45 680 : params.set<bool>("use_displaced_mesh") = false;
46 1360 : params.addRequiredParam<unsigned int>("ring_index", "Ring ID");
47 1360 : MooseEnum q_function_type("Geometry Topology", "Geometry");
48 1360 : params.addParam<MooseEnum>("q_function_type",
49 : q_function_type,
50 680 : "The method used to define the integration domain. Options are: " +
51 680 : q_function_type.getRawNames());
52 680 : 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 680 : return params;
56 680 : }
57 :
58 534 : JIntegral::JIntegral(const InputParameters & parameters)
59 : : ElementVectorPostprocessor(parameters),
60 534 : _crack_front_definition(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
61 1068 : _integral(getParam<MooseEnum>("integral").getEnum<IntegralType>()),
62 534 : _J_thermal_term_vec(hasMaterialProperty<RealVectorValue>("J_thermal_term_vec")
63 1602 : ? &getMaterialProperty<RealVectorValue>("J_thermal_term_vec")
64 : : nullptr),
65 1068 : _Eshelby_tensor(_integral != IntegralType::CIntegral
66 1036 : ? &getMaterialProperty<RankTwoTensor>("Eshelby_tensor")
67 : : nullptr),
68 534 : _Eshelby_tensor_dissipation(
69 534 : _integral == IntegralType::CIntegral
70 566 : ? &getMaterialProperty<RankTwoTensor>("Eshelby_tensor_dissipation")
71 : : nullptr),
72 534 : _fe_vars(getCoupledMooseVars()),
73 534 : _fe_type(_fe_vars[0]->feType()),
74 1068 : _has_symmetry_plane(isParamValid("symmetry_plane")),
75 1100 : _poissons_ratio(isParamValid("poissons_ratio") ? getParam<Real>("poissons_ratio") : 0),
76 1100 : _youngs_modulus(isParamValid("youngs_modulus") ? getParam<Real>("youngs_modulus") : 0),
77 1068 : _ring_index(getParam<unsigned int>("ring_index")),
78 1068 : _q_function_type(getParam<MooseEnum>("q_function_type").getEnum<QMethod>()),
79 1068 : _position_type(getParam<MooseEnum>("position_type").getEnum<PositionType>()),
80 534 : _x(declareVector("x")),
81 534 : _y(declareVector("y")),
82 534 : _z(declareVector("z")),
83 534 : _position(declareVector("id")),
84 534 : _j_integral(declareVector((_integral == IntegralType::KFromJIntegral
85 : ? "K_"
86 566 : : (_integral == IntegralType::JIntegral ? "J_" : "C_")) +
87 1068 : Moose::stringify(_ring_index)))
88 : {
89 534 : }
90 :
91 : void
92 534 : JIntegral::initialSetup()
93 : {
94 534 : _treat_as_2d = _crack_front_definition->treatAs2D();
95 534 : _using_mesh_cutter = _crack_front_definition->usingMeshCutter();
96 550 : if (_integral == IntegralType::KFromJIntegral &&
97 598 : (!isParamValid("youngs_modulus") || !isParamValid("poissons_ratio")))
98 0 : mooseError("youngs_modulus and poissons_ratio must be specified if integrals = KFromJIntegral");
99 534 : }
100 :
101 : void
102 675 : JIntegral::initialize()
103 : {
104 : std::size_t num_pts;
105 675 : if (_treat_as_2d && _using_mesh_cutter == false)
106 : num_pts = 1;
107 : else
108 180 : num_pts = _crack_front_definition->getNumCrackFrontPoints();
109 :
110 675 : _x.assign(num_pts, 0.0);
111 675 : _y.assign(num_pts, 0.0);
112 675 : _z.assign(num_pts, 0.0);
113 675 : _position.assign(num_pts, 0.0);
114 675 : _j_integral.assign(num_pts, 0.0);
115 :
116 2245 : for (const auto * fe_var : _fe_vars)
117 : {
118 1570 : if (fe_var->feType() != _fe_type)
119 0 : mooseError("displacements", "All coupled displacements must have the same type");
120 : }
121 675 : }
122 :
123 : Real
124 2023056 : JIntegral::computeQpIntegral(const std::size_t crack_front_point_index,
125 : const Real scalar_q,
126 : const RealVectorValue & grad_of_scalar_q)
127 : {
128 2023056 : RankTwoTensor grad_of_vector_q;
129 : const RealVectorValue & crack_direction =
130 2023056 : _crack_front_definition->getCrackDirection(crack_front_point_index);
131 2023056 : grad_of_vector_q(0, 0) = crack_direction(0) * grad_of_scalar_q(0);
132 2023056 : grad_of_vector_q(0, 1) = crack_direction(0) * grad_of_scalar_q(1);
133 2023056 : grad_of_vector_q(0, 2) = crack_direction(0) * grad_of_scalar_q(2);
134 2023056 : grad_of_vector_q(1, 0) = crack_direction(1) * grad_of_scalar_q(0);
135 2023056 : grad_of_vector_q(1, 1) = crack_direction(1) * grad_of_scalar_q(1);
136 2023056 : grad_of_vector_q(1, 2) = crack_direction(1) * grad_of_scalar_q(2);
137 2023056 : grad_of_vector_q(2, 0) = crack_direction(2) * grad_of_scalar_q(0);
138 2023056 : grad_of_vector_q(2, 1) = crack_direction(2) * grad_of_scalar_q(1);
139 2023056 : grad_of_vector_q(2, 2) = crack_direction(2) * grad_of_scalar_q(2);
140 :
141 : Real eq;
142 :
143 2023056 : if (_integral != IntegralType::CIntegral)
144 1997456 : eq = (*_Eshelby_tensor)[_qp].doubleContraction(grad_of_vector_q);
145 : else
146 25600 : eq = (*_Eshelby_tensor_dissipation)[_qp].doubleContraction(grad_of_vector_q);
147 :
148 : // Thermal component
149 : Real eq_thermal = 0.0;
150 2023056 : if (_J_thermal_term_vec && _integral != IntegralType::CIntegral)
151 : {
152 7989824 : for (std::size_t i = 0; i < 3; i++)
153 5992368 : eq_thermal += crack_direction(i) * scalar_q * (*_J_thermal_term_vec)[_qp](i);
154 : }
155 :
156 : Real q_avg_seg = 1.0;
157 2023056 : if (!_crack_front_definition->treatAs2D())
158 : {
159 1576416 : q_avg_seg =
160 3152832 : (_crack_front_definition->getCrackFrontForwardSegmentLength(crack_front_point_index) +
161 1576416 : _crack_front_definition->getCrackFrontBackwardSegmentLength(crack_front_point_index)) /
162 : 2.0;
163 : }
164 :
165 2023056 : Real etot = -eq + eq_thermal;
166 :
167 2023056 : return etot / q_avg_seg;
168 : }
169 :
170 : void
171 141600 : JIntegral::execute()
172 : {
173 : // calculate phi and dphi for this element
174 141600 : const std::size_t dim = _current_elem->dim();
175 141600 : std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(dim, _fe_type));
176 141600 : fe->attach_quadrature_rule(const_cast<QBase *>(_qrule));
177 141600 : _phi_curr_elem = &fe->get_phi();
178 141600 : _dphi_curr_elem = &fe->get_dphi();
179 141600 : fe->reinit(_current_elem);
180 :
181 : // calculate q for all nodes in this element
182 141600 : std::size_t ring_base = (_q_function_type == QMethod::Topology) ? 0 : 1;
183 :
184 449160 : for (std::size_t icfp = 0; icfp < _j_integral.size(); icfp++)
185 : {
186 307560 : _q_curr_elem.clear();
187 2330616 : for (std::size_t i = 0; i < _current_elem->n_nodes(); ++i)
188 : {
189 2023056 : const Node * this_node = _current_elem->node_ptr(i);
190 : Real q_this_node;
191 :
192 2023056 : if (_q_function_type == QMethod::Geometry)
193 1976592 : q_this_node = _crack_front_definition->DomainIntegralQFunction(
194 1976592 : icfp, _ring_index - ring_base, this_node);
195 : else // QMethod::Topology
196 46464 : q_this_node = _crack_front_definition->DomainIntegralTopologicalQFunction(
197 46464 : icfp, _ring_index - ring_base, this_node);
198 :
199 2023056 : _q_curr_elem.push_back(q_this_node);
200 : }
201 :
202 2330616 : 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 16457808 : for (std::size_t i = 0; i < _current_elem->n_nodes(); ++i)
208 : {
209 14434752 : scalar_q += (*_phi_curr_elem)[i][_qp] * _q_curr_elem[i];
210 :
211 55989312 : for (std::size_t j = 0; j < _current_elem->dim(); ++j)
212 41554560 : grad_of_scalar_q(j) += (*_dphi_curr_elem)[i][_qp](j) * _q_curr_elem[i];
213 : }
214 :
215 2023056 : _j_integral[icfp] +=
216 2023056 : _JxW[_qp] * _coord[_qp] * computeQpIntegral(icfp, scalar_q, grad_of_scalar_q);
217 : }
218 : }
219 141600 : }
220 :
221 : void
222 675 : JIntegral::finalize()
223 : {
224 675 : gatherSum(_j_integral);
225 :
226 2010 : for (std::size_t i = 0; i < _j_integral.size(); ++i)
227 : {
228 1335 : if (_has_symmetry_plane)
229 1180 : _j_integral[i] *= 2.0;
230 :
231 1335 : Real sign = (_j_integral[i] > 0.0) ? 1.0 : ((_j_integral[i] < 0.0) ? -1.0 : 0.0);
232 :
233 1335 : if (_integral == IntegralType::KFromJIntegral)
234 20 : _j_integral[i] = sign * std::sqrt(std::abs(_j_integral[i]) * _youngs_modulus /
235 20 : (1.0 - Utility::pow<2>(_poissons_ratio)));
236 :
237 1335 : const auto cfp = _crack_front_definition->getCrackFrontPoint(i);
238 1335 : _x[i] = (*cfp)(0);
239 1335 : _y[i] = (*cfp)(1);
240 1335 : _z[i] = (*cfp)(2);
241 :
242 1335 : if (_position_type == PositionType::Angle)
243 189 : _position[i] = _crack_front_definition->getAngleAlongFront(i);
244 : else
245 1146 : _position[i] = _crack_front_definition->getDistanceAlongFront(i);
246 : }
247 675 : }
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 : }
|
