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 : // MOOSE includes
11 : #include "ExplicitDynamicsContactConstraint.h"
12 : #include "FEProblem.h"
13 : #include "DisplacedProblem.h"
14 : #include "AuxiliarySystem.h"
15 : #include "PenetrationLocator.h"
16 : #include "NearestNodeLocator.h"
17 : #include "SystemBase.h"
18 : #include "Assembly.h"
19 : #include "MooseMesh.h"
20 : #include "MathUtils.h"
21 : #include "Executioner.h"
22 : #include "AddVariableAction.h"
23 : #include "ContactLineSearchBase.h"
24 : #include "ExplicitDynamicsContactAction.h"
25 :
26 : #include "libmesh/id_types.h"
27 : #include "libmesh/string_to_enum.h"
28 : #include "libmesh/sparse_matrix.h"
29 :
30 : registerMooseObject("ContactApp", ExplicitDynamicsContactConstraint);
31 :
32 : const unsigned int ExplicitDynamicsContactConstraint::_no_iterations = 0;
33 :
34 : InputParameters
35 360 : ExplicitDynamicsContactConstraint::validParams()
36 : {
37 360 : InputParameters params = NodeFaceConstraint::validParams();
38 360 : params += ExplicitDynamicsContactAction::commonParameters();
39 360 : params += TwoMaterialPropertyInterface::validParams();
40 :
41 720 : params.addRequiredParam<BoundaryName>("boundary", "The primary boundary");
42 720 : params.addParam<BoundaryName>("secondary", "The secondary boundary");
43 720 : params.addRequiredParam<unsigned int>("component",
44 : "An integer corresponding to the direction "
45 : "the variable this constraint acts on. (0 for x, "
46 : "1 for y, 2 for z)");
47 720 : params.addCoupledVar(
48 : "displacements",
49 : "The displacements appropriate for the simulation geometry and coordinate system");
50 720 : params.addCoupledVar("secondary_gap_offset", "offset to the gap distance from secondary side");
51 720 : params.addCoupledVar("mapped_primary_gap_offset",
52 : "offset to the gap distance mapped from primary side");
53 360 : params.set<bool>("use_displaced_mesh") = true;
54 720 : params.addParam<Real>("penalty",
55 720 : 1e8,
56 : "The penalty to apply. Its optimal value can vary depending on the "
57 : "stiffness of the materials");
58 720 : params.addParam<Real>("friction_coefficient", 0, "The friction coefficient");
59 720 : params.addParam<Real>("tangential_tolerance",
60 : "Tangential distance to extend edges of contact surfaces");
61 720 : params.addParam<bool>(
62 720 : "print_contact_nodes", false, "Whether to print the number of nodes in contact.");
63 :
64 360 : params.addClassDescription(
65 : "Apply non-penetration constraints on the mechanical deformation in explicit dynamics "
66 : "using a node on face formulation by solving uncoupled momentum-balance equations.");
67 360 : return params;
68 0 : }
69 :
70 270 : ExplicitDynamicsContactConstraint::ExplicitDynamicsContactConstraint(
71 270 : const InputParameters & parameters)
72 : : NodeFaceConstraint(parameters),
73 : TwoMaterialPropertyInterface(this, Moose::EMPTY_BLOCK_IDS, buildBoundaryIDs()),
74 270 : _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
75 540 : _component(getParam<unsigned int>("component")),
76 540 : _model(getParam<MooseEnum>("model").getEnum<ExplicitDynamicsContactModel>()),
77 270 : _update_stateful_data(true),
78 270 : _mesh_dimension(_mesh.dimension()),
79 270 : _vars(3, libMesh::invalid_uint),
80 270 : _var_objects(3, nullptr),
81 270 : _has_secondary_gap_offset(isCoupled("secondary_gap_offset")),
82 270 : _secondary_gap_offset_var(_has_secondary_gap_offset ? getVar("secondary_gap_offset", 0)
83 : : nullptr),
84 270 : _has_mapped_primary_gap_offset(isCoupled("mapped_primary_gap_offset")),
85 270 : _mapped_primary_gap_offset_var(
86 270 : _has_mapped_primary_gap_offset ? getVar("mapped_primary_gap_offset", 0) : nullptr),
87 540 : _penalty(getParam<Real>("penalty")),
88 540 : _print_contact_nodes(getParam<bool>("print_contact_nodes")),
89 270 : _residual_copy(_sys.residualGhosted()),
90 270 : _gap_rate(&writableVariable("gap_rate")),
91 486 : _neighbor_vel_x(isCoupled("vel_x") ? coupledNeighborValue("vel_x") : _zero),
92 486 : _neighbor_vel_y(isCoupled("vel_y") ? coupledNeighborValue("vel_y") : _zero),
93 756 : _neighbor_vel_z((_mesh.dimension() == 3 && isCoupled("vel_z")) ? coupledNeighborValue("vel_z")
94 270 : : _zero)
95 : {
96 270 : _overwrite_secondary_residual = false;
97 :
98 540 : if (isParamValid("displacements"))
99 : {
100 : // modern parameter scheme for displacements
101 2160 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
102 : {
103 810 : _vars[i] = coupled("displacements", i);
104 810 : _var_objects[i] = getVar("displacements", i);
105 : }
106 : }
107 :
108 540 : if (parameters.isParamValid("tangential_tolerance"))
109 0 : _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
110 :
111 540 : if (parameters.isParamValid("normal_smoothing_distance"))
112 0 : _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
113 :
114 540 : if (parameters.isParamValid("normal_smoothing_method"))
115 0 : _penetration_locator.setNormalSmoothingMethod(
116 : parameters.get<std::string>("normal_smoothing_method"));
117 :
118 270 : if (_model == ExplicitDynamicsContactModel::FRICTIONLESS_BALANCE)
119 : {
120 : bool is_correct =
121 648 : (isCoupled("vel_x") && isCoupled("vel_y") && _mesh.dimension() == 2) ||
122 1080 : (isCoupled("vel_x") && isCoupled("vel_y") && isCoupled("vel_z") && _mesh.dimension() == 3);
123 :
124 216 : if (!is_correct)
125 0 : paramError("vel_x",
126 : "Velocities vel_x and vel_y (also vel_z in three dimensions) need to be provided "
127 : "for the 'balance' option of solving normal contact in explicit dynamics.");
128 : }
129 270 : }
130 :
131 : void
132 35178 : ExplicitDynamicsContactConstraint::timestepSetup()
133 : {
134 35178 : if (_component == 0)
135 : {
136 11726 : updateContactStatefulData(/* beginning_of_step = */ true);
137 11726 : _update_stateful_data = false;
138 : _dof_to_position.clear();
139 : }
140 35178 : }
141 :
142 : void
143 11726 : ExplicitDynamicsContactConstraint::updateContactStatefulData(bool beginning_of_step)
144 : {
145 187810 : for (auto & [secondary_node_num, pinfo] : _penetration_locator._penetration_info)
146 : {
147 176084 : if (!pinfo)
148 0 : continue;
149 :
150 176084 : const Node & node = _mesh.nodeRef(secondary_node_num);
151 176084 : if (node.n_comp(_sys.number(), _vars[_component]) < 1)
152 0 : continue;
153 :
154 176084 : if (beginning_of_step)
155 : {
156 176084 : if (_app.getExecutioner()->lastSolveConverged())
157 : {
158 176084 : pinfo->_contact_force_old = pinfo->_contact_force;
159 176084 : pinfo->_accumulated_slip_old = pinfo->_accumulated_slip;
160 176084 : pinfo->_frictional_energy_old = pinfo->_frictional_energy;
161 176084 : pinfo->_mech_status_old = pinfo->_mech_status;
162 : }
163 0 : else if (pinfo->_mech_status_old == PenetrationInfo::MS_NO_CONTACT &&
164 0 : pinfo->_mech_status != PenetrationInfo::MS_NO_CONTACT)
165 : {
166 0 : mooseWarning("Previous step did not converge. Check results");
167 : // The penetration info object could be based on a bad state so delete it
168 0 : delete pinfo;
169 0 : pinfo = nullptr;
170 0 : continue;
171 : }
172 :
173 176084 : pinfo->_starting_elem = pinfo->_elem;
174 176084 : pinfo->_starting_side_num = pinfo->_side_num;
175 176084 : pinfo->_starting_closest_point_ref = pinfo->_closest_point_ref;
176 : }
177 176084 : pinfo->_incremental_slip_prev_iter = pinfo->_incremental_slip;
178 : }
179 11726 : }
180 :
181 : bool
182 337632 : ExplicitDynamicsContactConstraint::shouldApply()
183 : {
184 337632 : if (_current_node->processor_id() != _fe_problem.processor_id())
185 : return false;
186 :
187 : bool in_contact = false;
188 :
189 : std::map<dof_id_type, PenetrationInfo *>::iterator found =
190 337632 : _penetration_locator._penetration_info.find(_current_node->id());
191 337632 : if (found != _penetration_locator._penetration_info.end())
192 : {
193 337632 : PenetrationInfo * pinfo = found->second;
194 337632 : if (pinfo != nullptr)
195 : {
196 : // This computes the contact force once per constraint, rather than once per quad point
197 : // and for both primary and secondary cases.
198 337632 : if (_component == 0)
199 112544 : computeContactForce(*_current_node, pinfo, true);
200 :
201 337632 : if (pinfo->isCaptured())
202 : in_contact = true;
203 : }
204 : }
205 :
206 : return in_contact;
207 : }
208 :
209 : void
210 112544 : ExplicitDynamicsContactConstraint::computeContactForce(const Node & node,
211 : PenetrationInfo * pinfo,
212 : bool update_contact_set)
213 : {
214 : RealVectorValue distance_vec(node - pinfo->_closest_point);
215 :
216 : // Adding current velocity to the distance vector to ensure proper contact check
217 112544 : dof_id_type dof_x = node.dof_number(_sys.number(), _var_objects[0]->number(), 0);
218 112544 : dof_id_type dof_y = node.dof_number(_sys.number(), _var_objects[1]->number(), 0);
219 : dof_id_type dof_z =
220 112544 : _mesh.dimension() == 3 ? node.dof_number(_sys.number(), _var_objects[2]->number(), 0) : 0;
221 112544 : RealVectorValue udotvec = {(*_sys.solutionUDot())(dof_x)*_dt,
222 112544 : (*_sys.solutionUDot())(dof_y)*_dt,
223 112544 : _mesh.dimension() == 3 ? (*_sys.solutionUDot())(dof_z)*_dt : 0};
224 : distance_vec += udotvec;
225 :
226 112544 : if (distance_vec.norm() != 0)
227 112544 : distance_vec += gapOffset(node) * pinfo->_normal;
228 :
229 : const Real gap_size = -1.0 * pinfo->_normal * distance_vec;
230 :
231 : // This is for preventing an increment of pinfo->_locked_this_step for nodes that are
232 : // captured and released in this function
233 : bool newly_captured = false;
234 :
235 : // Capture nodes that are newly in contact
236 112544 : if (update_contact_set && !pinfo->isCaptured() &&
237 : MooseUtils::absoluteFuzzyGreaterEqual(gap_size, 0.0, 0.0))
238 : {
239 : newly_captured = true;
240 : pinfo->capture();
241 : }
242 :
243 112544 : if (!pinfo->isCaptured())
244 23680 : return;
245 :
246 88864 : switch (_model)
247 : {
248 2592 : case ExplicitDynamicsContactModel::FRICTIONLESS:
249 : {
250 2592 : const Real penalty = getPenalty(node);
251 : RealVectorValue pen_force(penalty * distance_vec);
252 2592 : pinfo->_contact_force = pinfo->_normal * (pinfo->_normal * pen_force);
253 : break;
254 : }
255 86272 : case ExplicitDynamicsContactModel::FRICTIONLESS_BALANCE:
256 86272 : solveImpactEquations(node, pinfo, distance_vec);
257 : break;
258 0 : default:
259 0 : mooseError("Invalid or unavailable contact model");
260 : break;
261 : }
262 :
263 88864 : if (update_contact_set && pinfo->isCaptured() && !newly_captured)
264 : {
265 : const Real contact_pressure = -(pinfo->_normal * pinfo->_contact_force);
266 88280 : if (-contact_pressure >= 0.0)
267 : {
268 : pinfo->release();
269 : pinfo->_contact_force.zero();
270 : }
271 : }
272 : }
273 :
274 : void
275 86272 : ExplicitDynamicsContactConstraint::solveImpactEquations(const Node & node,
276 : PenetrationInfo * pinfo,
277 : const RealVectorValue & /*distance_gap*/)
278 : {
279 : // Momentum balance, uncoupled normal pressure
280 : // See Heinstein et al, 2000, Contact-impact modeling in explicit transient dynamics.
281 :
282 86272 : dof_id_type dof_x = node.dof_number(_sys.number(), _var_objects[0]->number(), 0);
283 86272 : dof_id_type dof_y = node.dof_number(_sys.number(), _var_objects[1]->number(), 0);
284 : dof_id_type dof_z =
285 86272 : _mesh_dimension == 3 ? node.dof_number(_sys.number(), _var_objects[2]->number(), 0) : 0;
286 :
287 86272 : auto & u_dot = *_sys.solutionUDot();
288 :
289 : // Get lumped mass value
290 86272 : const auto & diag = _sys.getVector("mass_matrix_diag_inverted");
291 :
292 86272 : Real mass_node = 1.0 / diag(dof_x);
293 86272 : Real mass_face = computeFaceMass(diag);
294 :
295 86272 : Real mass_eff = (mass_face * mass_node) / (mass_face + mass_node);
296 :
297 : // Include effects of other forces:
298 : // Initial guess: v_{n-1/2} + dt * M^{-1} * (F^{ext} - F^{int})
299 86272 : Real velocity_x = u_dot(dof_x) + _dt / mass_node * -1 * _residual_copy(dof_x);
300 86272 : Real velocity_y = u_dot(dof_y) + _dt / mass_node * -1 * _residual_copy(dof_y);
301 86272 : Real velocity_z = u_dot(dof_z) + _dt / mass_node * -1 * _residual_copy(dof_z);
302 :
303 86272 : Real n_velocity_x = _neighbor_vel_x[0];
304 86272 : Real n_velocity_y = _neighbor_vel_y[0];
305 86272 : Real n_velocity_z = _neighbor_vel_z[0];
306 :
307 : RealVectorValue secondary_velocity(
308 86272 : velocity_x, velocity_y, _mesh.dimension() == 3 ? velocity_z : 0.0);
309 : RealVectorValue closest_point_velocity(
310 86272 : n_velocity_x, n_velocity_y, _mesh.dimension() == 3 ? n_velocity_z : 0.0);
311 : // Compute initial gap rate
312 86272 : Real gap_rate = pinfo->_normal * (secondary_velocity - closest_point_velocity);
313 :
314 : // Compute the force increment needed to set gap rate to 0
315 86272 : RealVectorValue impulse_force = pinfo->_normal * (gap_rate * mass_eff) / _dt;
316 86272 : pinfo->_contact_force = impulse_force;
317 :
318 : // recalculate velocity to determine gap rate
319 86272 : velocity_x -= _dt / mass_eff * impulse_force(0);
320 86272 : velocity_y -= _dt / mass_eff * impulse_force(1);
321 86272 : velocity_z -= _dt / mass_eff * impulse_force(2);
322 :
323 : // Recalculate gap rate for backwards compatibility
324 86272 : secondary_velocity = {velocity_x, velocity_y, _mesh.dimension() == 3 ? velocity_z : 0.0};
325 86272 : gap_rate = pinfo->_normal * (secondary_velocity - closest_point_velocity);
326 : // gap rate is now always near "0"
327 86272 : _gap_rate->setNodalValue(gap_rate);
328 86272 : }
329 :
330 : Real
331 0 : ExplicitDynamicsContactConstraint::computeQpSecondaryValue()
332 : {
333 : // Not used in current implementation.
334 0 : return _u_secondary[_qp];
335 : }
336 :
337 : Real
338 2399328 : ExplicitDynamicsContactConstraint::computeQpResidual(Moose::ConstraintType type)
339 : {
340 : // We use kinematic contact. But adding the residual helps avoid element inversion.
341 2399328 : PenetrationInfo * pinfo = _penetration_locator._penetration_info[_current_node->id()];
342 2399328 : Real resid = pinfo->_contact_force(_component);
343 :
344 2399328 : switch (type)
345 : {
346 266592 : case Moose::Secondary:
347 266592 : return _test_secondary[_i][_qp] * resid;
348 :
349 2132736 : case Moose::Primary:
350 2132736 : return _test_primary[_i][_qp] * -resid;
351 : }
352 :
353 : return 0.0;
354 : }
355 :
356 : Real
357 112544 : ExplicitDynamicsContactConstraint::gapOffset(const Node & node)
358 : {
359 : Real val = 0;
360 :
361 112544 : if (_has_secondary_gap_offset)
362 0 : val += _secondary_gap_offset_var->getNodalValue(node);
363 :
364 112544 : if (_has_mapped_primary_gap_offset)
365 0 : val += _mapped_primary_gap_offset_var->getNodalValue(node);
366 :
367 112544 : return val;
368 : }
369 :
370 : Real
371 2592 : ExplicitDynamicsContactConstraint::getPenalty(const Node & /*node*/)
372 : {
373 : // TODO: Include normalized penalty values.
374 2592 : return _penalty;
375 : }
376 :
377 : Real
378 86272 : ExplicitDynamicsContactConstraint::computeFaceMass(const NumericVector<Real> & lumped_mass)
379 : {
380 : // Initialize face mass to zero
381 : Real mass_face(0.0);
382 :
383 : // Get the primary side of the current contact
384 86272 : const auto primary_side = _current_primary->side_ptr(_assembly.neighborSide());
385 :
386 : // Get the dofs on the primary (face) side of the contact
387 : std::vector<dof_id_type> face_dofs;
388 86272 : _primary_var.getDofIndices(primary_side.get(), face_dofs);
389 :
390 : // Get average mass of face
391 431360 : for (const auto dof : face_dofs)
392 345088 : mass_face += 1.0 / lumped_mass(dof);
393 :
394 86272 : mass_face /= face_dofs.size();
395 :
396 86272 : return mass_face;
397 86272 : }
|
