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 "InteractionIntegral.h"
11 : #include "MooseMesh.h"
12 : #include "RankTwoTensor.h"
13 : #include "Conversion.h"
14 : #include "libmesh/string_to_enum.h"
15 : #include "libmesh/quadrature.h"
16 : #include "DerivativeMaterialInterface.h"
17 : #include "libmesh/utility.h"
18 : #include "CrackFrontDefinition.h"
19 :
20 : registerMooseObject("SolidMechanicsApp", InteractionIntegral);
21 : registerMooseObject("SolidMechanicsApp", ADInteractionIntegral);
22 :
23 : template <bool is_ad>
24 : MooseEnum
25 1928 : InteractionIntegralTempl<is_ad>::qFunctionType()
26 : {
27 3856 : return MooseEnum("Geometry Topology", "Geometry");
28 : }
29 :
30 : template <bool is_ad>
31 : MooseEnum
32 1928 : InteractionIntegralTempl<is_ad>::sifModeType()
33 : {
34 3856 : return MooseEnum("KI KII KIII T", "KI");
35 : }
36 :
37 : template <bool is_ad>
38 : InputParameters
39 964 : InteractionIntegralTempl<is_ad>::validParams()
40 : {
41 964 : InputParameters params = ElementVectorPostprocessor::validParams();
42 1928 : params.addRequiredCoupledVar(
43 : "displacements",
44 : "The displacements appropriate for the simulation geometry and coordinate system");
45 1928 : params.addCoupledVar("temperature",
46 : "The temperature (optional). Must be provided to correctly compute "
47 : "stress intensity factors in models with thermal strain gradients.");
48 1928 : params.addParam<MaterialPropertyName>(
49 : "functionally_graded_youngs_modulus",
50 : "Spatially varying elasticity modulus variable. This input is required when "
51 : "using the functionally graded material capability.");
52 1928 : params.addParam<MaterialPropertyName>(
53 : "functionally_graded_youngs_modulus_crack_dir_gradient",
54 : "Gradient of the spatially varying Young's modulus provided in "
55 : "'functionally_graded_youngs_modulus' in the direction of crack extension.");
56 1928 : params.addRequiredParam<UserObjectName>("crack_front_definition",
57 : "The CrackFrontDefinition user object name");
58 1928 : MooseEnum position_type("Angle Distance", "Distance");
59 1928 : params.addParam<MooseEnum>(
60 : "position_type",
61 : position_type,
62 : "The method used to calculate position along crack front. Options are: " +
63 : position_type.getRawNames());
64 1928 : params.addRequiredParam<Real>(
65 : "K_factor", "Conversion factor between interaction integral and stress intensity factor K");
66 1928 : params.addParam<unsigned int>("symmetry_plane",
67 : "Account for a symmetry plane passing through "
68 : "the plane of the crack, normal to the specified "
69 : "axis (0=x, 1=y, 2=z)");
70 1928 : params.addRequiredParam<Real>("poissons_ratio", "Poisson's ratio for the material.");
71 1928 : params.addRequiredParam<Real>("youngs_modulus", "Young's modulus of the material.");
72 964 : params.set<bool>("use_displaced_mesh") = false;
73 1928 : params.addRequiredParam<unsigned int>("ring_index", "Ring ID");
74 2892 : params.addParam<MooseEnum>("q_function_type",
75 : InteractionIntegralTempl<is_ad>::qFunctionType(),
76 : "The method used to define the integration domain. Options are: " +
77 : InteractionIntegralTempl<is_ad>::qFunctionType().getRawNames());
78 2892 : params.addRequiredParam<MooseEnum>(
79 : "sif_mode",
80 : InteractionIntegralTempl<is_ad>::sifModeType(),
81 : "Stress intensity factor to calculate. Choices are: " +
82 : InteractionIntegralTempl<is_ad>::sifModeType().getRawNames());
83 1928 : params.addParam<MaterialPropertyName>("eigenstrain_gradient",
84 : "Material defining gradient of eigenstrain tensor");
85 1928 : params.addParam<MaterialPropertyName>("body_force", "Material defining body force");
86 1928 : params.addCoupledVar("additional_eigenstrain_00",
87 : "Optional additional eigenstrain variable that will be accounted for in the "
88 : "interaction integral (component 00 or XX).");
89 1928 : params.addCoupledVar("additional_eigenstrain_01",
90 : "Optional additional eigenstrain variable that will be accounted for in the "
91 : "interaction integral (component 01 or XY).");
92 1928 : params.addCoupledVar("additional_eigenstrain_11",
93 : "Optional additional eigenstrain variable that will be accounted for in the "
94 : "interaction integral (component 11 or YY).");
95 1928 : params.addCoupledVar("additional_eigenstrain_22",
96 : "Optional additional eigenstrain variable that will be accounted for in the "
97 : "interaction integral (component 22 or ZZ).");
98 1928 : params.addCoupledVar("additional_eigenstrain_02",
99 : "Optional additional eigenstrain variable that will be accounted for in the "
100 : "interaction integral (component 02 or XZ).");
101 1928 : params.addCoupledVar("additional_eigenstrain_12",
102 : "Optional additional eigenstrain variable that will be accounted for in the "
103 : "interaction integral (component 12 or XZ).");
104 1928 : params.addParamNamesToGroup(
105 : "additional_eigenstrain_00 additional_eigenstrain_01 additional_eigenstrain_11 "
106 : "additional_eigenstrain_22 additional_eigenstrain_02 additional_eigenstrain_12",
107 : "Generic eigenstrains for the computation of the interaction integral.");
108 964 : params.addClassDescription(
109 : "Computes the interaction integral, which is used to compute various "
110 : "fracture mechanics parameters at a crack tip, including KI, KII, KIII, "
111 : "and the T stress.");
112 964 : return params;
113 964 : }
114 :
115 : template <bool is_ad>
116 840 : InteractionIntegralTempl<is_ad>::InteractionIntegralTempl(const InputParameters & parameters)
117 : : ElementVectorPostprocessor(parameters),
118 840 : _ndisp(coupledComponents("displacements")),
119 840 : _crack_front_definition(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
120 840 : _treat_as_2d(false),
121 1680 : _stress(getGenericMaterialPropertyByName<RankTwoTensor, is_ad>("stress")),
122 840 : _strain(getGenericMaterialPropertyByName<RankTwoTensor, is_ad>("elastic_strain")),
123 840 : _fe_vars(getCoupledMooseVars()),
124 840 : _fe_type(_fe_vars[0]->feType()),
125 840 : _disp(coupledValues("displacements")),
126 840 : _grad_disp(3),
127 840 : _has_temp(isCoupled("temperature")),
128 840 : _grad_temp(_has_temp ? coupledGradient("temperature") : _grad_zero),
129 840 : _functionally_graded_youngs_modulus_crack_dir_gradient(
130 840 : isParamSetByUser("functionally_graded_youngs_modulus_crack_dir_gradient")
131 924 : ? &getMaterialProperty<Real>("functionally_graded_youngs_modulus_crack_dir_gradient")
132 : : nullptr),
133 840 : _functionally_graded_youngs_modulus(
134 840 : isParamSetByUser("functionally_graded_youngs_modulus")
135 924 : ? &getMaterialProperty<Real>("functionally_graded_youngs_modulus")
136 : : nullptr),
137 1680 : _K_factor(getParam<Real>("K_factor")),
138 1680 : _has_symmetry_plane(isParamValid("symmetry_plane")),
139 1680 : _poissons_ratio(getParam<Real>("poissons_ratio")),
140 1680 : _youngs_modulus(getParam<Real>("youngs_modulus")),
141 1680 : _fgm_crack(isParamSetByUser("functionally_graded_youngs_modulus_crack_dir_gradient") &&
142 924 : isParamSetByUser("functionally_graded_youngs_modulus")),
143 1680 : _ring_index(getParam<unsigned int>("ring_index")),
144 840 : _total_deigenstrain_dT(
145 840 : hasMaterialProperty<RankTwoTensor>("total_deigenstrain_dT")
146 936 : ? &getGenericMaterialProperty<RankTwoTensor, is_ad>("total_deigenstrain_dT")
147 : : nullptr),
148 840 : _has_additional_eigenstrain(false),
149 840 : _additional_eigenstrain_gradient_00(isCoupled("additional_eigenstrain_00")
150 864 : ? &coupledGradient("additional_eigenstrain_00")
151 : : nullptr),
152 840 : _additional_eigenstrain_gradient_01(isCoupled("additional_eigenstrain_01")
153 864 : ? &coupledGradient("additional_eigenstrain_01")
154 : : nullptr),
155 840 : _additional_eigenstrain_gradient_11(isCoupled("additional_eigenstrain_11")
156 864 : ? &coupledGradient("additional_eigenstrain_11")
157 : : nullptr),
158 840 : _additional_eigenstrain_gradient_22(isCoupled("additional_eigenstrain_22")
159 864 : ? &coupledGradient("additional_eigenstrain_22")
160 : : nullptr),
161 840 : _additional_eigenstrain_gradient_02(isCoupled("additional_eigenstrain_02")
162 840 : ? &coupledGradient("additional_eigenstrain_02")
163 : : nullptr),
164 1680 : _additional_eigenstrain_gradient_12(isCoupled("additional_eigenstrain_12")
165 840 : ? &coupledGradient("additional_eigenstrain_12")
166 : : nullptr),
167 1680 : _q_function_type(getParam<MooseEnum>("q_function_type").template getEnum<QMethod>()),
168 1680 : _position_type(getParam<MooseEnum>("position_type").template getEnum<PositionType>()),
169 1680 : _sif_mode(getParam<MooseEnum>("sif_mode").template getEnum<SifMethod>()),
170 840 : _x(declareVector("x")),
171 840 : _y(declareVector("y")),
172 840 : _z(declareVector("z")),
173 840 : _position(declareVector("id")),
174 3360 : _interaction_integral(declareVector("II_" + Moose::stringify(getParam<MooseEnum>("sif_mode")) +
175 : "_" + Moose::stringify(_ring_index))),
176 840 : _eigenstrain_gradient(nullptr),
177 840 : _body_force(nullptr)
178 : {
179 840 : if (_has_temp && !_total_deigenstrain_dT)
180 0 : mooseError("InteractionIntegral Error: To include thermal strain term in interaction integral, "
181 : "must both couple temperature in DomainIntegral block and compute "
182 : "total_deigenstrain_dT using ThermalFractureIntegral material model.");
183 :
184 840 : if ((!_functionally_graded_youngs_modulus &&
185 840 : _functionally_graded_youngs_modulus_crack_dir_gradient) ||
186 42 : (_functionally_graded_youngs_modulus &&
187 42 : !_functionally_graded_youngs_modulus_crack_dir_gradient))
188 0 : paramError("functionally_graded_youngs_modulus_crack_dir_gradient",
189 : "You have selected to compute the interaction integral for a crack in FGM. That "
190 : "selection requires the user to provide a spatially varying elasticity modulus "
191 : "that "
192 : "defines the transition of material properties (i.e. "
193 : "'functionally_graded_youngs_modulus') and its "
194 : "spatial derivative in the crack direction (i.e. "
195 : "'functionally_graded_youngs_modulus_crack_dir_gradient').");
196 :
197 : // plane strain
198 840 : _kappa = 3.0 - 4.0 * _poissons_ratio;
199 840 : _shear_modulus = _youngs_modulus / (2.0 * (1.0 + _poissons_ratio));
200 :
201 : // Checking for consistency between mesh size and length of the provided displacements vector
202 840 : if (_ndisp != _mesh.dimension())
203 0 : paramError("displacements",
204 : "InteractionIntegral Error: number of variables supplied in 'displacements' must "
205 : "match the mesh dimension.");
206 :
207 : // fetch gradients of coupled variables
208 2898 : for (std::size_t i = 0; i < _ndisp; ++i)
209 2058 : _grad_disp[i] = &coupledGradient("displacements", i);
210 :
211 : // set unused dimensions to zero
212 1302 : for (std::size_t i = _ndisp; i < 3; ++i)
213 462 : _grad_disp[i] = &_grad_zero;
214 :
215 1680 : if (isParamValid("eigenstrain_gradient"))
216 : {
217 24 : _eigenstrain_gradient =
218 24 : &getGenericMaterialProperty<RankThreeTensor, is_ad>("eigenstrain_gradient");
219 24 : if (_total_deigenstrain_dT)
220 0 : paramError("eigenstrain_gradient",
221 : "eigenstrain_gradient cannot be specified for materials that provide the "
222 : "total_deigenstrain_dT material property");
223 : }
224 1680 : if (isParamValid("body_force"))
225 48 : _body_force = &getMaterialProperty<RealVectorValue>("body_force");
226 :
227 1704 : _has_additional_eigenstrain = _additional_eigenstrain_gradient_00 &&
228 840 : _additional_eigenstrain_gradient_01 &&
229 24 : _additional_eigenstrain_gradient_11;
230 :
231 840 : if (_mesh.dimension() == 3 && _has_additional_eigenstrain)
232 0 : if (!(_additional_eigenstrain_gradient_22 && _additional_eigenstrain_gradient_12 &&
233 0 : _additional_eigenstrain_gradient_02))
234 0 : paramError("additional_eigenstrain_gradient_12",
235 : "If additional eigenstrains are provided for the computation of the interaction "
236 : "integral in three dimensions, make sure the six components are provided.");
237 :
238 840 : if (_has_additional_eigenstrain)
239 24 : mooseInfo("A generic eigenstrain provided by the user will be considered in the interaction "
240 : "integral (via domain integral action).");
241 840 : }
242 :
243 : template <bool is_ad>
244 : void
245 840 : InteractionIntegralTempl<is_ad>::initialSetup()
246 : {
247 840 : _treat_as_2d = _crack_front_definition->treatAs2D();
248 840 : _using_mesh_cutter = _crack_front_definition->usingMeshCutter();
249 840 : }
250 :
251 : template <bool is_ad>
252 : void
253 840 : InteractionIntegralTempl<is_ad>::initialize()
254 : {
255 : std::size_t num_pts;
256 840 : if (_treat_as_2d && _using_mesh_cutter == false)
257 : num_pts = 1;
258 : else
259 258 : num_pts = _crack_front_definition->getNumCrackFrontPoints();
260 :
261 840 : _x.assign(num_pts, 0.0);
262 840 : _y.assign(num_pts, 0.0);
263 840 : _z.assign(num_pts, 0.0);
264 840 : _position.assign(num_pts, 0.0);
265 840 : _interaction_integral.assign(num_pts, 0.0);
266 :
267 3042 : for (const auto * fe_var : _fe_vars)
268 : {
269 2202 : if (fe_var->feType() != _fe_type)
270 0 : mooseError("All coupled displacements must have the same type");
271 : }
272 840 : }
273 :
274 : template <bool is_ad>
275 : Real
276 1717168 : InteractionIntegralTempl<is_ad>::computeQpIntegral(const std::size_t crack_front_point_index,
277 : const Real scalar_q,
278 : const RealVectorValue & grad_of_scalar_q)
279 : {
280 : // If q is zero, then dq is also zero, so all terms in the interaction integral would
281 : // return zero. As such, let us avoid unnecessary, frequent computations
282 1717168 : if (scalar_q < TOLERANCE * TOLERANCE * TOLERANCE)
283 : return 0.0;
284 :
285 : // In the crack front coordinate system, the crack direction is (1,0,0)
286 : RealVectorValue crack_direction(1.0, 0.0, 0.0);
287 :
288 : // Calculate (r,theta) position of qp relative to crack front
289 97904 : _crack_front_definition->calculateRThetaToCrackFront(
290 97904 : _q_point[_qp], crack_front_point_index, _r, _theta);
291 :
292 97904 : RankTwoTensor aux_stress;
293 97904 : RankTwoTensor aux_du;
294 97904 : RankTwoTensor aux_strain;
295 97904 : RankTwoTensor aux_disp;
296 :
297 97904 : if (_sif_mode == SifMethod::KI || _sif_mode == SifMethod::KII || _sif_mode == SifMethod::KIII)
298 88352 : computeAuxFields(aux_stress, aux_du, aux_strain, aux_disp);
299 9552 : else if (_sif_mode == SifMethod::T)
300 9552 : computeTFields(aux_stress, aux_du);
301 :
302 97904 : auto grad_disp = RankTwoTensor::initializeFromRows(
303 : (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
304 :
305 : // Rotate stress, strain, displacement and temperature to crack front coordinate system
306 : RealVectorValue grad_q_cf =
307 97904 : _crack_front_definition->rotateToCrackFrontCoords(grad_of_scalar_q, crack_front_point_index);
308 195808 : RankTwoTensor grad_disp_cf =
309 97904 : _crack_front_definition->rotateToCrackFrontCoords(grad_disp, crack_front_point_index);
310 97904 : RankTwoTensor stress_cf = _crack_front_definition->rotateToCrackFrontCoords(
311 97904 : MetaPhysicL::raw_value((_stress)[_qp]), crack_front_point_index);
312 97904 : RankTwoTensor strain_cf = _crack_front_definition->rotateToCrackFrontCoords(
313 97904 : MetaPhysicL::raw_value((_strain)[_qp]), crack_front_point_index);
314 : RealVectorValue grad_temp_cf =
315 97904 : _crack_front_definition->rotateToCrackFrontCoords(_grad_temp[_qp], crack_front_point_index);
316 :
317 97904 : if (getBlockCoordSystem() == Moose::COORD_RZ)
318 1008 : strain_cf(2, 2) = 0.0;
319 :
320 97904 : RankTwoTensor dq;
321 97904 : dq(0, 0) = crack_direction(0) * grad_q_cf(0);
322 97904 : dq(0, 1) = crack_direction(0) * grad_q_cf(1);
323 97904 : dq(0, 2) = crack_direction(0) * grad_q_cf(2);
324 :
325 : // Calculate interaction integral terms
326 :
327 : // Term1 = stress * x1-derivative of aux disp * dq
328 97904 : RankTwoTensor tmp1 = dq * stress_cf;
329 97904 : Real term1 = aux_du.doubleContraction(tmp1);
330 :
331 : // Term2 = aux stress * x1-derivative of disp * dq
332 97904 : RankTwoTensor tmp2 = dq * aux_stress;
333 97904 : Real term2 = grad_disp_cf(0, 0) * tmp2(0, 0) + grad_disp_cf(1, 0) * tmp2(0, 1) +
334 97904 : grad_disp_cf(2, 0) * tmp2(0, 2);
335 :
336 : // Term3 = aux stress * strain * dq_x (= stress * aux strain * dq_x)
337 97904 : Real term3 = -dq(0, 0) * aux_stress.doubleContraction(strain_cf);
338 :
339 : // Term4 (thermal strain term) = q * aux_stress * alpha * dtheta_x
340 : // - the term including the derivative of alpha is not implemented
341 : Real term4 = 0.0;
342 97904 : if (_has_temp)
343 : {
344 : Real sigma_alpha =
345 1536 : aux_stress.doubleContraction(MetaPhysicL::raw_value((*_total_deigenstrain_dT)[_qp]));
346 1536 : term4 = scalar_q * sigma_alpha * grad_temp_cf(0);
347 : }
348 :
349 : Real term4a = 0.0; // Gradient of arbitrary eigenstrain material contribution
350 97904 : if (_eigenstrain_gradient)
351 : {
352 : // Thermal strain gradient term in Nakamura and Parks, IJSS, 1992:
353 : // alpha * dT/dx_k*aux_stress*scalar_q
354 : // Generalization to the gradient of an arbitrary eigenstrain:
355 : // d_eigenstrain/dx_k*aux_stress*scalar_q
356 : const RealVectorValue & crack_dir =
357 768 : _crack_front_definition->getCrackDirection(crack_front_point_index);
358 768 : RankTwoTensor eigenstrain_grad_in_crack_dir =
359 768 : crack_dir * MetaPhysicL::raw_value((*_eigenstrain_gradient)[_qp]);
360 768 : term4a = scalar_q * aux_stress.doubleContraction(eigenstrain_grad_in_crack_dir);
361 : }
362 :
363 : // Gradient of arbitrary eigenstrain variable contribution.
364 : // Note that, unlike approaches 4 and 4a, approach 4b can integrate any type of eigenstrain, i.e.
365 : // it does not have to be thermal or volumetric.
366 : Real term4b = 0.0;
367 : // Term4b (eigen strain term) = q * aux_stress * d{eigen_strain}_x
368 97904 : if (_has_additional_eigenstrain)
369 : {
370 768 : RankTwoTensor I_x(RankTwoTensor::initNone);
371 768 : I_x(0, 0) = 1.0;
372 768 : RankTwoTensor I_xy(RankTwoTensor::initNone);
373 768 : I_xy(0, 1) = 1.0;
374 768 : I_xy(1, 0) = 1.0;
375 768 : RankTwoTensor I_y(RankTwoTensor::initNone);
376 768 : I_y(1, 1) = 1.0;
377 768 : RankTwoTensor I_z(RankTwoTensor::initNone);
378 768 : I_z(2, 2) = 1.0;
379 :
380 : const RealVectorValue & crack_dir =
381 768 : _crack_front_definition->getCrackDirection(crack_front_point_index);
382 :
383 : const auto eigenstrain_gradient_x =
384 768 : I_x.mixedProductJkI((*_additional_eigenstrain_gradient_00)[_qp]);
385 : const auto eigenstrain_gradient_xy =
386 768 : I_xy.mixedProductJkI((*_additional_eigenstrain_gradient_01)[_qp]);
387 : const auto eigenstrain_gradient_y =
388 768 : I_y.mixedProductJkI((*_additional_eigenstrain_gradient_11)[_qp]);
389 :
390 768 : GenericRankThreeTensor<is_ad> eigenstrain_gradient_z;
391 768 : if (_additional_eigenstrain_gradient_22)
392 768 : eigenstrain_gradient_z = I_z.mixedProductJkI((*_additional_eigenstrain_gradient_22)[_qp]);
393 :
394 768 : RankTwoTensor eigenstrain_grad_in_crack_dir =
395 768 : crack_dir * (MetaPhysicL::raw_value(eigenstrain_gradient_x) +
396 768 : MetaPhysicL::raw_value(eigenstrain_gradient_xy) +
397 768 : MetaPhysicL::raw_value(eigenstrain_gradient_y) +
398 768 : MetaPhysicL::raw_value(eigenstrain_gradient_z));
399 :
400 768 : GenericRankThreeTensor<is_ad> eigenstrain_gradient_xz;
401 768 : GenericRankThreeTensor<is_ad> eigenstrain_gradient_yz;
402 : // Add terms for three dimensions
403 768 : if (_mesh.dimension() == 3)
404 : {
405 : // Add influence of xz and yz.
406 0 : RankTwoTensor I_xz(RankTwoTensor::initNone);
407 0 : I_xz(0, 2) = I_xz(2, 0) = 1.0;
408 0 : eigenstrain_gradient_xz = I_xz.mixedProductJkI((*_additional_eigenstrain_gradient_02)[_qp]);
409 0 : RankTwoTensor I_yz(RankTwoTensor::initNone);
410 0 : I_yz(1, 2) = I_yz(2, 1) = 1.0;
411 0 : eigenstrain_gradient_yz = I_yz.mixedProductJkI((*_additional_eigenstrain_gradient_12)[_qp]);
412 :
413 0 : eigenstrain_grad_in_crack_dir = crack_dir * (MetaPhysicL::raw_value(eigenstrain_gradient_xz) +
414 0 : MetaPhysicL::raw_value(eigenstrain_gradient_yz) +
415 0 : MetaPhysicL::raw_value(eigenstrain_gradient_x) +
416 0 : MetaPhysicL::raw_value(eigenstrain_gradient_xy) +
417 0 : MetaPhysicL::raw_value(eigenstrain_gradient_y) +
418 0 : MetaPhysicL::raw_value(eigenstrain_gradient_z));
419 : }
420 :
421 768 : term4b = scalar_q * aux_stress.doubleContraction(eigenstrain_grad_in_crack_dir);
422 : }
423 :
424 : Real term5 = 0.0; // Body force
425 97904 : if (_body_force)
426 : {
427 : // Body force term in Nakamura and Parks, IJSS, 1992:
428 : // b_i*aux_du*crack_dir*scalar_q
429 : const RealVectorValue & crack_dir =
430 768 : _crack_front_definition->getCrackDirection(crack_front_point_index);
431 1536 : term5 = -scalar_q * MetaPhysicL::raw_value((*_body_force)[_qp]) * aux_du * crack_dir;
432 : }
433 :
434 : Real term6 = 0.0;
435 97904 : if (_fgm_crack && scalar_q != 0)
436 : {
437 : // See Equation 49 of J.-H. Kim and G. H. Paulino. Consistent formulations of the interaction
438 : // integral method for fracture of functionally graded materials. Journal of Applied Mechanics,
439 : // 72(3) 351-364 2004.
440 :
441 : // Term 6_1 = Cijkl,j(x) epsilon_{kl}^{aux} disp_{i,1}
442 4448 : RankTwoTensor cijklj_epsilonkl_aux; // Temporary second order tensor.
443 4448 : RankFourTensor cijklj;
444 4448 : cijklj.fillSymmetricIsotropicEandNu(1.0, _poissons_ratio);
445 4448 : cijklj *= (*_functionally_graded_youngs_modulus_crack_dir_gradient)[_qp];
446 4448 : cijklj_epsilonkl_aux = cijklj * aux_strain;
447 :
448 4448 : const Real term6_a = grad_disp_cf(0, 0) * cijklj_epsilonkl_aux(0, 0) +
449 4448 : grad_disp_cf(1, 0) * cijklj_epsilonkl_aux(0, 1) +
450 4448 : grad_disp_cf(2, 0) * cijklj_epsilonkl_aux(0, 2);
451 :
452 : // Term 6_2 = -Cijkl,1(x) epsilon_{kl} epsilon_{ij}^{aux}
453 4448 : RankTwoTensor cijkl1_epsilonkl;
454 4448 : cijkl1_epsilonkl = -cijklj * strain_cf;
455 4448 : const Real term6_b = cijkl1_epsilonkl.doubleContraction(aux_strain);
456 :
457 4448 : term6 = (term6_a + term6_b) * scalar_q;
458 : }
459 :
460 : // Add terms for possibly RZ (axisymmetric problem).
461 : // See Nahta and Moran, 1993, Domain integrals for axisymmetric interface crack problems, 30(15)
462 : // Note: Sign problem in Eq. (26). Also Kolosov constant is wrong in Eq. (13)
463 : Real term7 = 0.0;
464 97904 : if (getBlockCoordSystem() == Moose::COORD_RZ)
465 : {
466 1008 : const Real u_r = (*_disp[0])[_qp];
467 1008 : const Real radius_qp = _q_point[_qp](0);
468 :
469 1008 : const Real term7a =
470 2016 : (u_r / radius_qp * aux_stress(2, 2) + aux_disp(0, 0) / radius_qp * stress_cf(2, 2) -
471 1008 : aux_stress.doubleContraction(strain_cf)) /
472 : radius_qp;
473 :
474 : // term7b1: sigma_aux_ralpha * u_alpha,r
475 1008 : const Real term7b1 = aux_stress(0, 0) * grad_disp_cf(0, 0) +
476 1008 : aux_stress(0, 1) * grad_disp_cf(1, 0) +
477 1008 : aux_stress(0, 2) * grad_disp_cf(2, 0);
478 :
479 1008 : const Real term7b = 1.0 / radius_qp *
480 1008 : (term7b1 - aux_stress(2, 2) * grad_disp_cf(0, 0) +
481 1008 : stress_cf(2, 2) * (aux_du(0, 0) - aux_disp(0, 0) / radius_qp));
482 :
483 1008 : term7 = (term7a + term7b) * scalar_q;
484 : }
485 :
486 : Real q_avg_seg = 1.0;
487 97904 : if (!_crack_front_definition->treatAs2D())
488 : {
489 52416 : q_avg_seg =
490 104832 : (_crack_front_definition->getCrackFrontForwardSegmentLength(crack_front_point_index) +
491 52416 : _crack_front_definition->getCrackFrontBackwardSegmentLength(crack_front_point_index)) /
492 : 2.0;
493 : }
494 :
495 97904 : Real eq = term1 + term2 + term3 + term4 + term4a + term4b + term5 + term6 + term7;
496 :
497 97904 : if (getBlockCoordSystem() == Moose::COORD_RZ)
498 : {
499 1008 : const Real radius_qp = _q_point[_qp](0);
500 1008 : eq *= radius_qp;
501 :
502 1008 : std::size_t num_crack_front_points = _crack_front_definition->getNumCrackFrontPoints();
503 1008 : if (num_crack_front_points != 1)
504 0 : mooseError("Crack front has more than one point, but this is a 2D-RZ problem. Please revise "
505 : "your input file.");
506 :
507 1008 : const Point * crack_front_point = _crack_front_definition->getCrackFrontPoint(0);
508 1008 : eq = eq / (*crack_front_point)(0);
509 : }
510 :
511 97904 : return eq / q_avg_seg;
512 : }
513 :
514 : template <bool is_ad>
515 : void
516 149248 : InteractionIntegralTempl<is_ad>::execute()
517 : {
518 : // calculate phi and dphi for this element
519 149248 : const std::size_t dim = _current_elem->dim();
520 149248 : std::unique_ptr<FEBase> fe(FEBase::build(dim, _fe_type));
521 149248 : fe->attach_quadrature_rule(const_cast<QBase *>(_qrule));
522 149248 : _phi_curr_elem = &fe->get_phi();
523 149248 : _dphi_curr_elem = &fe->get_dphi();
524 149248 : fe->reinit(_current_elem);
525 :
526 : // calculate q for all nodes in this element
527 149248 : std::size_t ring_base = (_q_function_type == QMethod::Topology) ? 0 : 1;
528 :
529 423224 : for (auto icfp = beginIndex(_interaction_integral); icfp < _interaction_integral.size(); icfp++)
530 : {
531 273976 : _q_curr_elem.clear();
532 1991144 : for (std::size_t i = 0; i < _current_elem->n_nodes(); ++i)
533 : {
534 1717168 : const Node * this_node = _current_elem->node_ptr(i);
535 : Real q_this_node;
536 :
537 1717168 : if (_q_function_type == QMethod::Geometry)
538 1717168 : q_this_node = _crack_front_definition->DomainIntegralQFunction(
539 1717168 : icfp, _ring_index - ring_base, this_node);
540 : else // QMethod::Topology
541 0 : q_this_node = _crack_front_definition->DomainIntegralTopologicalQFunction(
542 0 : icfp, _ring_index - ring_base, this_node);
543 :
544 1717168 : _q_curr_elem.push_back(q_this_node);
545 : }
546 :
547 1991144 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
548 : {
549 : Real scalar_q = 0.0;
550 : RealVectorValue grad_of_scalar_q(0.0, 0.0, 0.0);
551 :
552 13555952 : for (std::size_t i = 0; i < _current_elem->n_nodes(); ++i)
553 : {
554 11838784 : scalar_q += (*_phi_curr_elem)[i][_qp] * _q_curr_elem[i];
555 :
556 45456576 : for (std::size_t j = 0; j < _current_elem->dim(); ++j)
557 33617792 : grad_of_scalar_q(j) += (*_dphi_curr_elem)[i][_qp](j) * _q_curr_elem[i];
558 : }
559 :
560 1717168 : if (getBlockCoordSystem() != Moose::COORD_RZ)
561 1607008 : _interaction_integral[icfp] +=
562 1607008 : _JxW[_qp] * _coord[_qp] * computeQpIntegral(icfp, scalar_q, grad_of_scalar_q);
563 : else
564 110160 : _interaction_integral[icfp] +=
565 110160 : _JxW[_qp] * computeQpIntegral(icfp, scalar_q, grad_of_scalar_q);
566 : }
567 : }
568 149248 : }
569 :
570 : template <bool is_ad>
571 : void
572 840 : InteractionIntegralTempl<is_ad>::finalize()
573 : {
574 840 : gatherSum(_interaction_integral);
575 :
576 2484 : for (auto i = beginIndex(_interaction_integral); i < _interaction_integral.size(); ++i)
577 : {
578 1644 : if (_has_symmetry_plane)
579 780 : _interaction_integral[i] *= 2.0;
580 :
581 1644 : const auto cfp = _crack_front_definition->getCrackFrontPoint(i);
582 1644 : _x[i] = (*cfp)(0);
583 1644 : _y[i] = (*cfp)(1);
584 1644 : _z[i] = (*cfp)(2);
585 :
586 1644 : if (_position_type == PositionType::Angle)
587 0 : _position[i] = _crack_front_definition->getAngleAlongFront(i);
588 : else
589 1644 : _position[i] = _crack_front_definition->getDistanceAlongFront(i);
590 :
591 1644 : if (_sif_mode == SifMethod::T && !_treat_as_2d)
592 252 : _interaction_integral[i] +=
593 252 : _poissons_ratio * _crack_front_definition->getCrackFrontTangentialStrain(i);
594 :
595 1644 : _interaction_integral[i] *= _K_factor;
596 : }
597 840 : }
598 :
599 : template <bool is_ad>
600 : void
601 0 : InteractionIntegralTempl<is_ad>::threadJoin(const UserObject & y)
602 : {
603 : const auto & uo = static_cast<const InteractionIntegralTempl<is_ad> &>(y);
604 :
605 0 : for (auto i = beginIndex(_interaction_integral); i < _interaction_integral.size(); ++i)
606 0 : _interaction_integral[i] += uo._interaction_integral[i];
607 0 : }
608 :
609 : template <bool is_ad>
610 : void
611 88352 : InteractionIntegralTempl<is_ad>::computeAuxFields(RankTwoTensor & aux_stress,
612 : RankTwoTensor & grad_disp,
613 : RankTwoTensor & aux_strain,
614 : RankTwoTensor & aux_disp)
615 : {
616 : RealVectorValue k(0.0);
617 88352 : if (_sif_mode == SifMethod::KI)
618 : k(0) = 1.0;
619 45824 : else if (_sif_mode == SifMethod::KII)
620 : k(1) = 1.0;
621 22912 : else if (_sif_mode == SifMethod::KIII)
622 : k(2) = 1.0;
623 :
624 88352 : Real t = _theta;
625 88352 : Real t2 = _theta / 2.0;
626 88352 : Real tt2 = 3.0 * _theta / 2.0;
627 88352 : Real st = std::sin(t);
628 88352 : Real ct = std::cos(t);
629 88352 : Real st2 = std::sin(t2);
630 88352 : Real ct2 = std::cos(t2);
631 88352 : Real stt2 = std::sin(tt2);
632 88352 : Real ctt2 = std::cos(tt2);
633 : Real ct2sq = Utility::pow<2>(ct2);
634 : Real st2sq = Utility::pow<2>(st2);
635 : Real ct2cu = Utility::pow<3>(ct2);
636 88352 : Real sqrt2PiR = std::sqrt(2.0 * libMesh::pi * _r);
637 :
638 : // Calculate auxiliary stress tensor
639 : aux_stress.zero();
640 :
641 88352 : aux_stress(0, 0) =
642 88352 : 1.0 / sqrt2PiR * (k(0) * ct2 * (1.0 - st2 * stt2) - k(1) * st2 * (2.0 + ct2 * ctt2));
643 88352 : aux_stress(1, 1) = 1.0 / sqrt2PiR * (k(0) * ct2 * (1.0 + st2 * stt2) + k(1) * st2 * ct2 * ctt2);
644 88352 : aux_stress(0, 1) = 1.0 / sqrt2PiR * (k(0) * ct2 * st2 * ctt2 + k(1) * ct2 * (1.0 - st2 * stt2));
645 88352 : aux_stress(0, 2) = -1.0 / sqrt2PiR * k(2) * st2;
646 88352 : aux_stress(1, 2) = 1.0 / sqrt2PiR * k(2) * ct2;
647 : // plane stress
648 : // Real s33 = 0;
649 : // plane strain
650 88352 : aux_stress(2, 2) = _poissons_ratio * (aux_stress(0, 0) + aux_stress(1, 1));
651 :
652 88352 : aux_stress(1, 0) = aux_stress(0, 1);
653 88352 : aux_stress(2, 0) = aux_stress(0, 2);
654 88352 : aux_stress(2, 1) = aux_stress(1, 2);
655 :
656 : // Calculate auxiliary strain tensor (only needed to FGM)
657 88352 : if (_fgm_crack)
658 : {
659 : aux_strain.zero();
660 4448 : RankFourTensor spatial_elasticity_tensor;
661 4448 : spatial_elasticity_tensor.fillSymmetricIsotropicEandNu(
662 4448 : (*_functionally_graded_youngs_modulus)[_qp], _poissons_ratio);
663 4448 : aux_strain = spatial_elasticity_tensor.invSymm() * aux_stress;
664 : }
665 :
666 : // Calculate x1 derivative of auxiliary displacements
667 : grad_disp.zero();
668 88352 : grad_disp(0, 0) = k(0) / (4.0 * _shear_modulus * sqrt2PiR) *
669 88352 : (ct * ct2 * _kappa + ct * ct2 - 2.0 * ct * ct2cu + st * st2 * _kappa +
670 88352 : st * st2 - 6.0 * st * st2 * ct2sq) +
671 88352 : k(1) / (4.0 * _shear_modulus * sqrt2PiR) *
672 88352 : (ct * st2 * _kappa + ct * st2 + 2.0 * ct * st2 * ct2sq - st * ct2 * _kappa +
673 88352 : 3.0 * st * ct2 - 6.0 * st * ct2cu);
674 :
675 88352 : grad_disp(0, 1) = k(0) / (4.0 * _shear_modulus * sqrt2PiR) *
676 88352 : (ct * st2 * _kappa + ct * st2 - 2.0 * ct * st2 * ct2sq - st * ct2 * _kappa -
677 88352 : 5.0 * st * ct2 + 6.0 * st * ct2cu) +
678 88352 : k(1) / (4.0 * _shear_modulus * sqrt2PiR) *
679 88352 : (-ct * ct2 * _kappa + 3.0 * ct * ct2 - 2.0 * ct * ct2cu -
680 88352 : st * st2 * _kappa + 3.0 * st * st2 - 6.0 * st * st2 * ct2sq);
681 :
682 88352 : grad_disp(0, 2) = k(2) / (_shear_modulus * sqrt2PiR) * (st2 * ct - ct2 * st);
683 :
684 : // Calculate aux displacement field (used in RZ calculations) - Warp3D manual
685 : aux_disp.zero();
686 88352 : if (getBlockCoordSystem() == Moose::COORD_RZ)
687 : {
688 1008 : const Real prefactor_rz = std::sqrt(_r / 2.0 / libMesh::pi) / 2.0 / _shear_modulus;
689 1008 : aux_disp(0, 0) = prefactor_rz * ((k(0) * ct2 * (_kappa - 1.0 + 2.0 * st2sq)) +
690 1008 : (k(1) * st2 * (_kappa + 1.0 + 2.0 * ct2sq)));
691 1008 : aux_disp(0, 1) = prefactor_rz * ((k(0) * st2 * (_kappa + 1.0 - 2.0 * ct2sq)) -
692 1008 : (k(1) * ct2 * (_kappa - 1.0 - 2.0 * st2sq)));
693 1008 : aux_disp(0, 2) = prefactor_rz * k(2) * 4.0 * st2;
694 : }
695 88352 : }
696 :
697 : template <bool is_ad>
698 : void
699 9552 : InteractionIntegralTempl<is_ad>::computeTFields(RankTwoTensor & aux_stress,
700 : RankTwoTensor & grad_disp)
701 : {
702 :
703 9552 : Real t = _theta;
704 9552 : Real st = std::sin(t);
705 9552 : Real ct = std::cos(t);
706 : Real stsq = Utility::pow<2>(st);
707 : Real ctsq = Utility::pow<2>(ct);
708 : Real ctcu = Utility::pow<3>(ct);
709 9552 : Real oneOverPiR = 1.0 / (libMesh::pi * _r);
710 :
711 : aux_stress.zero();
712 9552 : aux_stress(0, 0) = -oneOverPiR * ctcu;
713 9552 : aux_stress(0, 1) = -oneOverPiR * st * ctsq;
714 9552 : aux_stress(1, 0) = -oneOverPiR * st * ctsq;
715 9552 : aux_stress(1, 1) = -oneOverPiR * ct * stsq;
716 9552 : aux_stress(2, 2) = -oneOverPiR * _poissons_ratio * (ctcu + ct * stsq);
717 :
718 : grad_disp.zero();
719 9552 : grad_disp(0, 0) = oneOverPiR / (4.0 * _youngs_modulus) *
720 9552 : (ct * (4.0 * Utility::pow<2>(_poissons_ratio) - 3.0 + _poissons_ratio) -
721 9552 : std::cos(3.0 * t) * (1.0 + _poissons_ratio));
722 9552 : grad_disp(0, 1) = -oneOverPiR / (4.0 * _youngs_modulus) *
723 9552 : (st * (4.0 * Utility::pow<2>(_poissons_ratio) - 3.0 + _poissons_ratio) +
724 9552 : std::sin(3.0 * t) * (1.0 + _poissons_ratio));
725 9552 : }
726 :
727 : template class InteractionIntegralTempl<false>;
728 : template class InteractionIntegralTempl<true>;
|
