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 "BilinearMixedModeCohesiveZoneModel.h"
11 : #include "MooseVariableFE.h"
12 : #include "SystemBase.h"
13 : #include "MortarUtils.h"
14 : #include "MooseUtils.h"
15 : #include "MathUtils.h"
16 :
17 : #include "MortarContactUtils.h"
18 : #include "FactorizedRankTwoTensor.h"
19 :
20 : #include "ADReal.h"
21 : #include "CohesiveZoneModelTools.h"
22 : #include <Eigen/Core>
23 :
24 : registerMooseObject("ContactApp", BilinearMixedModeCohesiveZoneModel);
25 :
26 : InputParameters
27 90 : BilinearMixedModeCohesiveZoneModel::validParams()
28 : {
29 90 : InputParameters params = CohesiveZoneModelBase::validParams();
30 :
31 90 : params.addClassDescription("Computes the bilinear mixed mode cohesive zone model.");
32 :
33 : // Input parameters for bilinear mixed mode traction.
34 180 : params.addParam<MaterialPropertyName>("GI_c",
35 : "Critical energy release rate in normal direction.");
36 180 : params.addParam<MaterialPropertyName>("GII_c",
37 : "Critical energy release rate in shear direction.");
38 180 : params.addParam<MaterialPropertyName>("normal_strength", "Tensile strength in normal direction.");
39 180 : params.addParam<MaterialPropertyName>("shear_strength", "Tensile strength in shear direction.");
40 180 : params.addParam<Real>("power_law_parameter", "The power law parameter.");
41 180 : MooseEnum criterion("POWER_LAW BK", "BK");
42 180 : params.addParam<Real>("viscosity", 0.0, "Viscosity for damage model.");
43 180 : params.addParam<MooseEnum>(
44 : "mixed_mode_criterion", criterion, "Option for mixed mode propagation criterion.");
45 180 : params.addParam<bool>(
46 : "lag_displacement_jump",
47 180 : false,
48 : "Whether to use old displacement jumps to compute the effective displacement jump.");
49 180 : params.addParam<bool>("set_compressive_traction_to_zero",
50 180 : false,
51 : "Zero compressive traction (set to true, allowing the use of standard "
52 : "zero-penetration mortar contact constraints in "
53 : "the normal direction).");
54 180 : params.addParam<Real>(
55 180 : "regularization_alpha", 1e-10, "Regularization parameter for the Macaulay bracket.");
56 180 : params.addRangeCheckedParam<Real>(
57 : "penalty_stiffness", "penalty_stiffness > 0.0", "Penalty stiffness for CZM.");
58 180 : params.addParamNamesToGroup(
59 : "GI_c GII_c normal_strength shear_strength power_law_parameter viscosity "
60 : "mixed_mode_criterion lag_displacement_jump regularization_alpha "
61 : "penalty_stiffness",
62 : "Bilinear mixed mode traction");
63 : // End of input parameters for bilinear mixed mode traction.
64 :
65 90 : return params;
66 90 : }
67 :
68 45 : BilinearMixedModeCohesiveZoneModel::BilinearMixedModeCohesiveZoneModel(
69 45 : const InputParameters & parameters)
70 : : WeightedGapUserObject(parameters),
71 : PenaltyWeightedGapUserObject(parameters),
72 : WeightedVelocitiesUserObject(parameters),
73 : CohesiveZoneModelBase(parameters),
74 45 : _set_compressive_traction_to_zero(getParam<bool>("set_compressive_traction_to_zero")),
75 90 : _normal_strength(getMaterialProperty<Real>("normal_strength")),
76 90 : _shear_strength(getMaterialProperty<Real>("shear_strength")),
77 90 : _GI_c(getMaterialProperty<Real>("GI_c")),
78 90 : _GII_c(getMaterialProperty<Real>("GII_c")),
79 90 : _penalty_stiffness_czm(getParam<Real>("penalty_stiffness")),
80 90 : _mix_mode_criterion(getParam<MooseEnum>("mixed_mode_criterion").getEnum<MixedModeCriterion>()),
81 90 : _power_law_parameter(getParam<Real>("power_law_parameter")),
82 90 : _viscosity(getParam<Real>("viscosity")),
83 135 : _regularization_alpha(getParam<Real>("regularization_alpha"))
84 : {
85 :
86 : // Checks for bilinear traction input.
87 315 : if (!isParamValid("GI_c") || !isParamValid("GII_c") || !isParamValid("normal_strength") ||
88 315 : !isParamValid("shear_strength") || !isParamValid("power_law_parameter") ||
89 135 : !isParamValid("penalty_stiffness"))
90 0 : paramError("GI_c",
91 : "The CZM bilinear mixed mode traction parameters GI_c, GII_c, normal_strength, "
92 : "shear_strength, and power_law_parameter are required. Revise your input and add "
93 : "those parameters if you want to use the bilinear mixed mode traction model. ");
94 45 : }
95 :
96 : void
97 26886 : BilinearMixedModeCohesiveZoneModel::computeQpProperties()
98 : {
99 26886 : CohesiveZoneModelBase::computeQpProperties();
100 :
101 26886 : _normal_strength_interpolation = _normal_strength[_qp] * _JxW_msm[_qp] * _coord[_qp];
102 26886 : _shear_strength_interpolation = _shear_strength[_qp] * _JxW_msm[_qp] * _coord[_qp];
103 26886 : _GI_c_interpolation = _GI_c[_qp] * _JxW_msm[_qp] * _coord[_qp];
104 26886 : _GII_c_interpolation = _GII_c[_qp] * _JxW_msm[_qp] * _coord[_qp];
105 26886 : }
106 :
107 : void
108 53772 : BilinearMixedModeCohesiveZoneModel::computeQpIProperties()
109 : {
110 53772 : CohesiveZoneModelBase::computeQpIProperties();
111 :
112 : // Get the _dof_to_weighted_gap map
113 53772 : const auto * const dof = static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i));
114 :
115 : // TODO: Probably better to interpolate the deformation gradients.
116 107544 : _dof_to_normal_strength[dof] += (*_test)[_i][_qp] * _normal_strength_interpolation;
117 107544 : _dof_to_shear_strength[dof] += (*_test)[_i][_qp] * _shear_strength_interpolation;
118 :
119 107544 : _dof_to_GI_c[dof] += (*_test)[_i][_qp] * _GI_c_interpolation;
120 107544 : _dof_to_GII_c[dof] += (*_test)[_i][_qp] * _GII_c_interpolation;
121 53772 : }
122 :
123 : void
124 7087 : BilinearMixedModeCohesiveZoneModel::initialize()
125 : {
126 : // instead we call it explicitly here
127 7087 : CohesiveZoneModelBase::initialize();
128 :
129 : // Avoid accumulating interpolation over the time step
130 : _dof_to_mode_mixity_ratio.clear();
131 : _dof_to_normal_strength.clear();
132 : _dof_to_shear_strength.clear();
133 : _dof_to_GI_c.clear();
134 : _dof_to_GII_c.clear();
135 : _dof_to_delta_initial.clear();
136 : _dof_to_delta_final.clear();
137 : _dof_to_delta_max.clear();
138 7087 : }
139 :
140 : void
141 26886 : BilinearMixedModeCohesiveZoneModel::computeCZMTraction(const Node * const node)
142 : {
143 : using std::max, std::min;
144 :
145 : // First call does not have maps available
146 26886 : const bool return_boolean = _dof_to_weighted_gap.find(node) == _dof_to_weighted_gap.end();
147 26886 : if (return_boolean)
148 0 : return;
149 :
150 26886 : computeModeMixity(node);
151 26886 : computeCriticalDisplacementJump(node);
152 26886 : computeFinalDisplacementJump(node);
153 26886 : computeEffectiveDisplacementJump(node);
154 26886 : computeDamage(node);
155 :
156 : // Split displacement jump into active and inactive parts
157 : const auto interface_displacement_jump =
158 26886 : normalizeQuantity(_dof_to_interface_displacement_jump, node);
159 :
160 26886 : const ADRealVectorValue delta_active(max(interface_displacement_jump(0), 0.0),
161 : interface_displacement_jump(1),
162 : interface_displacement_jump(2));
163 26886 : const ADRealVectorValue delta_inactive(min(interface_displacement_jump(0), 0.0), 0.0, 0.0);
164 :
165 : // This traction vector is local at this point.
166 53772 : _dof_to_czm_traction[node] =
167 53772 : -(1.0 - libmesh_map_find(_dof_to_damage, node->id()).first) * _penalty_stiffness_czm *
168 26886 : delta_active -
169 80658 : _penalty_stiffness_czm * (_set_compressive_traction_to_zero ? 0.0 : delta_inactive);
170 : }
171 :
172 : void
173 26886 : BilinearMixedModeCohesiveZoneModel::computeModeMixity(const Node * const node)
174 : {
175 : using std::sqrt;
176 :
177 : const auto interface_displacement_jump =
178 26886 : normalizeQuantity(_dof_to_interface_displacement_jump, node);
179 :
180 26886 : if (interface_displacement_jump(0) > _epsilon_tolerance)
181 : {
182 : const auto delta_s =
183 20525 : sqrt(interface_displacement_jump(1) * interface_displacement_jump(1) +
184 20525 : interface_displacement_jump(2) * interface_displacement_jump(2) + _epsilon_tolerance);
185 :
186 20525 : _dof_to_mode_mixity_ratio[node] = delta_s / interface_displacement_jump(0);
187 : }
188 : else
189 6361 : _dof_to_mode_mixity_ratio[node] = 0;
190 26886 : }
191 :
192 : void
193 26886 : BilinearMixedModeCohesiveZoneModel::computeCriticalDisplacementJump(const Node * const node)
194 : {
195 : using std::sqrt;
196 :
197 : const auto interface_displacement_jump =
198 26886 : normalizeQuantity(_dof_to_interface_displacement_jump, node);
199 :
200 26886 : const auto mixity_ratio = libmesh_map_find(_dof_to_mode_mixity_ratio, node);
201 :
202 : const auto delta_normal_knot =
203 26886 : normalizeQuantity(_dof_to_normal_strength, node) / _penalty_stiffness_czm;
204 : const auto delta_shear_knot =
205 26886 : normalizeQuantity(_dof_to_shear_strength, node) / _penalty_stiffness_czm;
206 :
207 26886 : _dof_to_delta_initial[node] = delta_shear_knot;
208 :
209 26886 : if (interface_displacement_jump(0) > _epsilon_tolerance)
210 : {
211 20525 : const auto delta_mixed = sqrt(delta_shear_knot * delta_shear_knot +
212 20525 : Utility::pow<2>(mixity_ratio * delta_normal_knot));
213 :
214 41050 : _dof_to_delta_initial[node] = delta_normal_knot * delta_shear_knot *
215 61575 : sqrt(1.0 + mixity_ratio * mixity_ratio) / delta_mixed;
216 : }
217 26886 : }
218 :
219 : void
220 26886 : BilinearMixedModeCohesiveZoneModel::computeFinalDisplacementJump(const Node * const node)
221 : {
222 : using std::sqrt, std::pow;
223 :
224 : const auto interface_displacement_jump =
225 26886 : normalizeQuantity(_dof_to_interface_displacement_jump, node);
226 :
227 26886 : const auto mixity_ratio = libmesh_map_find(_dof_to_mode_mixity_ratio, node);
228 :
229 26886 : const auto normalized_GI_c = normalizeQuantity(_dof_to_GI_c, node);
230 26886 : const auto normalized_GII_c = normalizeQuantity(_dof_to_GII_c, node);
231 :
232 53772 : _dof_to_delta_final[node] =
233 53772 : sqrt(2.0) * 2.0 * normalized_GII_c / normalizeQuantity(_dof_to_shear_strength, node);
234 :
235 26886 : if (interface_displacement_jump(0) > _epsilon_tolerance)
236 : {
237 20525 : if (_mix_mode_criterion == MixedModeCriterion::BK)
238 : {
239 14842 : _dof_to_delta_final[node] =
240 7421 : 2.0 / _penalty_stiffness_czm / libmesh_map_find(_dof_to_delta_initial, node) *
241 7421 : (normalized_GI_c +
242 7421 : (normalized_GII_c - normalized_GI_c) *
243 22263 : pow(mixity_ratio * mixity_ratio / (1 + mixity_ratio * mixity_ratio),
244 14842 : _power_law_parameter));
245 : }
246 13104 : else if (_mix_mode_criterion == MixedModeCriterion::POWER_LAW)
247 : {
248 : const auto Gc_mixed =
249 26208 : pow(1.0 / normalized_GI_c, _power_law_parameter) +
250 26208 : pow(mixity_ratio * mixity_ratio / normalized_GII_c, _power_law_parameter);
251 39312 : _dof_to_delta_final[node] = (2.0 + 2.0 * mixity_ratio * mixity_ratio) /
252 13104 : _penalty_stiffness_czm /
253 13104 : libmesh_map_find(_dof_to_delta_initial, node) *
254 26208 : pow(Gc_mixed, -1.0 / _power_law_parameter);
255 : }
256 : }
257 26886 : }
258 :
259 : void
260 26886 : BilinearMixedModeCohesiveZoneModel::computeEffectiveDisplacementJump(const Node * const node)
261 : {
262 : using std::sqrt;
263 :
264 : const auto interface_displacement_jump =
265 26886 : normalizeQuantity(_dof_to_interface_displacement_jump, node);
266 :
267 : const auto delta_normal_pos =
268 26886 : MathUtils::regularizedHeavyside(interface_displacement_jump(0), _regularization_alpha) *
269 : interface_displacement_jump(0);
270 :
271 26886 : _dof_to_delta_max[node] = sqrt(Utility::pow<2>(interface_displacement_jump(1)) +
272 26886 : Utility::pow<2>(interface_displacement_jump(2)) +
273 53772 : Utility::pow<2>(delta_normal_pos) + _epsilon_tolerance);
274 26886 : }
275 :
276 : void
277 26886 : BilinearMixedModeCohesiveZoneModel::computeDamage(const Node * const node)
278 : {
279 26886 : const auto delta_max = libmesh_map_find(_dof_to_delta_max, node);
280 26886 : const auto delta_initial = libmesh_map_find(_dof_to_delta_initial, node);
281 26886 : const auto delta_final = libmesh_map_find(_dof_to_delta_final, node);
282 :
283 26886 : auto & [damage, damage_old] = _dof_to_damage[node->id()];
284 26886 : if (delta_max < delta_initial)
285 22976 : damage = 0;
286 3910 : else if (delta_max > delta_final)
287 0 : damage = 1.0;
288 : else
289 3910 : damage = delta_final * (delta_max - delta_initial) / delta_max / (delta_final - delta_initial);
290 :
291 26886 : if (damage < damage_old)
292 : // Irreversibility
293 491 : damage = damage_old;
294 :
295 : // Viscous regularization
296 53772 : damage = (damage + _viscosity * damage_old / _dt) / (_viscosity / _dt + 1.0);
297 26886 : }
298 :
299 : void
300 7087 : BilinearMixedModeCohesiveZoneModel::finalize()
301 : {
302 7087 : CohesiveZoneModelBase::finalize();
303 :
304 7087 : const bool send_data_back = !constrainedByOwner();
305 :
306 14174 : Moose::Mortar::Contact::communicateRealObject(
307 7087 : _dof_to_normal_strength, _subproblem.mesh(), _nodal, _communicator, send_data_back);
308 :
309 14174 : Moose::Mortar::Contact::communicateRealObject(
310 7087 : _dof_to_shear_strength, _subproblem.mesh(), _nodal, _communicator, send_data_back);
311 :
312 14174 : Moose::Mortar::Contact::communicateRealObject(
313 7087 : _dof_to_GI_c, _subproblem.mesh(), _nodal, _communicator, send_data_back);
314 :
315 14174 : Moose::Mortar::Contact::communicateRealObject(
316 7087 : _dof_to_GII_c, _subproblem.mesh(), _nodal, _communicator, send_data_back);
317 7087 : }
318 :
319 : Real
320 1890 : BilinearMixedModeCohesiveZoneModel::getModeMixityRatio(const Node * const node) const
321 : {
322 1890 : const auto it = _dof_to_mode_mixity_ratio.find(_subproblem.mesh().nodePtr(node->id()));
323 :
324 1890 : if (it != _dof_to_mode_mixity_ratio.end())
325 1114 : return MetaPhysicL::raw_value(it->second);
326 : else
327 : return 0.0;
328 : }
329 :
330 : Real
331 1890 : BilinearMixedModeCohesiveZoneModel::getCohesiveDamage(const Node * const node) const
332 : {
333 1890 : const auto it = _dof_to_damage.find(node->id());
334 :
335 1890 : if (it != _dof_to_damage.end())
336 1114 : return MetaPhysicL::raw_value(it->second.first);
337 : else
338 : return 0.0;
339 : }
340 :
341 : Real
342 1890 : BilinearMixedModeCohesiveZoneModel::getLocalDisplacementNormal(const Node * const node) const
343 : {
344 1890 : const auto it = _dof_to_interface_displacement_jump.find(_subproblem.mesh().nodePtr(node->id()));
345 1890 : const auto it2 = _dof_to_weighted_gap.find(_subproblem.mesh().nodePtr(node->id()));
346 :
347 1890 : if (it != _dof_to_interface_displacement_jump.end() && it2 != _dof_to_weighted_gap.end())
348 2228 : return MetaPhysicL::raw_value(it->second(0) / it2->second.second);
349 : else
350 : return 0.0;
351 : }
352 :
353 : Real
354 1890 : BilinearMixedModeCohesiveZoneModel::getLocalDisplacementTangential(const Node * const node) const
355 : {
356 1890 : const auto it = _dof_to_interface_displacement_jump.find(_subproblem.mesh().nodePtr(node->id()));
357 1890 : const auto it2 = _dof_to_weighted_gap.find(_subproblem.mesh().nodePtr(node->id()));
358 :
359 1890 : if (it != _dof_to_interface_displacement_jump.end() && it2 != _dof_to_weighted_gap.end())
360 2228 : return MetaPhysicL::raw_value(it->second(1) / it2->second.second);
361 : else
362 : return 0.0;
363 : }
|
