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