Decomposition method for dynamic contact problems of several deformable bodies

The iterative method of Schwarz domain decomposition (DDS) was used as a dynamic contact algorithm for two and more contacting 3D bodies allowing the consideration of the solution of the stress-strain state problem in the standard statement for each of the contacting body and the high-level parallelization of the solution process. According to the offered iteration algorithm, two steps are performed in each iteration, where the contact conditions for the displacement and stresses are fulfilled in turn for the nodes on the contact surface. The specific feature of the realization of the dynamic contact problem was the combination of the iteration procedures for the refinement of the contact problem solution according to the Newton-Raphson algorithm with the iterations of the DDS method. As the test calculation results have shown, the use of the Schwarz contact algorithm and the HHT-α scheme in combination with the mass redistribution method on the contact boundary is a reliable, stable and sufficiently accurate way for the numerical integration of the contact problems at impact interaction. Finally, accuracy of the proposed method is verified by a conservation of momentum through three contact examples.


Introduction
In the case of dynamic contact problems, the using of the classical scheme of step-by-step time integration becomes complicated because the accurate fulfillment of the boundary conditions of the contact, when Newmark implicit scheme for time integration is applied, leads to unphysical oscillations at the significant deviation of the total energy in a system. The appearance of such oscillations is mainly connected with that, when the standard approximation of the mass matrix is used, all the nodes, including the nodes on the contact surface, are thought to have their own mass which contributes into the total energy of the system. At the point of time corresponding to the contact interaction, the nodes velocity on the interaction surface is zero and their kinetic energy is zero as well; this results in the redistribution of the energy contributions into the total energy of the system leading to considerable unphysical numerical fluctuations.
The conservative properties of the schemes are preserved only for the linear elasticity problem with a constant contact surface, whereas, for the case of the varying contact boundary, the appearance of a new node in the contact area decreases the energy of the discrete system, and the node departure from the contact increases the energy. The effective methods for solving 3D dynamic contact problems are characterized by the chosen scheme of the numerical time-integration which in much determines the presence of parasitic oscillations due to the impossibility of the sufficiently accurate reproduction of the contribution of the high-frequency forms of oscillations into the solution of the dynamic contact problem; the conservation of the total energy of the system of contacting bodies; the stability of solution. One of the ways to exclude parasitic oscillations in the contact region is the use of dissipative schemes. As the computing practice shows, it is effective to use generalized implicit methods based on Newmark scheme such as Hilber-Hughes-Taylor (HHT-) scheme [1], which is unconditionally stable, has the second order accuracy and is dissipative for high frequencies. However, the use of the above scheme for the considered mathematical model of the impact interaction has shown that the scheme is unstable.
The alternative approach for preventing contact oscillations is also the mass redistribution method on the contact boundary [2]. The method can be summarized as follows: the mass of nodes in the contact area is nullified and, as a result, the inertia characteristics of the nodes are excluded. The method significantly stabilizes an unknown contact boundary and can be used in combination with any numerical schemes.
In the present paper, the iterative method of Schwarz domain decomposition (DDS) [3] was used as a dynamic contact algorithm [4,5] for two and more contacting 3D bodies allowing the consideration of the solution of the stress-strain state problem in the standard statement for each of the contacting body and the high-level parallelization of the solution process. According to the offered iteration algorithm, two steps are performed in each iteration, where the contact conditions for the displacement and stresses are fulfilled in turn for the nodes on the contact surface. The specific feature of the realization of the dynamic contact problem was the combination of the iteration procedures for the refinement of the contact problem solution according to the Newton-Raphson algorithm with the iterations of the DDS method. As the test calculation results have shown, the use of the Schwarz contact algorithm and the HHT-scheme in combination with the mass redistribution method on the contact boundary is a reliable, stable and sufficiently accurate way for the numerical integration of the contact problems at impact interaction.
Finally, accuracy of the proposed method is verified by a conservation of momentum through three contact examples. The results of the DDS method are validated by comparison with the experimental results of direct central collisions of three identical elastic spheres. In this case, two successive impacts were observed among three spheres. The DDS method results are agreed very well with the experimental results. In the following example, five steel spheres arranged with a normal gap consistently collide. The initial velocity of the first body was half meter per second. Eight conditions of the contract were used.
The results showed that there is no rotation of the spheres and the contact algorithm ensures high accuracy of the results at large time intervals. The difference in the maximum values of the contact force for each pair contact is two percent.
Also, in this paper we present the validation of the three-dimensional finite-element model of a human head. This model is built with the using of voxel data and the computing algorithm which includes the solution of the contact problem by Schwarz domain decomposition method and the generalized Newmark implicit dissipative scheme with the mass redistribution on the contact boundary.
Finally, three examples were demonstrated to verify the accuracy of results using the proposed method.

Dynamic contact of deformed bodies
We consider contact problem of homogeneous and isotropic linearly-elastic bodies that occupy the domains Ω (1) , Ω (2) , … , Ω ( ) . These domains are bounded by the surfaces δΩ ( ) , ∈ [1, ], which are correspond to the undeformed state of the bodies. We denote the surfaces that bound the bodies Ω at time ∈ [0, ] as Γ ( ) = δΩ ( ) . Contact of bodies occurs over a part of these surfaces which we call contact areas Γ ⊂ ⋃ =1 Γ ( ) . We assume that contact areas change during deformation.
We will consider frictionless contact problem without adhesion. We use the notation and for the normal and tangential parts of the displacement and stress components on the contact boundary Γ . For the case of small increments of displacements, the correspondence bodies at time between the points on Γ ( ) is given by the following overlap function: ( ) is negative when there is no contact and positive in the case of bodies overlapping. We assume that for each of interacting bodies the equilibrium equations, Hooke's law and Cauchy's relations are fulfilled: On the boundaries of bodies kinematic and static boundary conditions are applied: The initial displacements and velocities are given: We use the overlapping function (1) to specify the conditions on the contact surface: The first condition in (5) provides non-penetration, the second is compressive contact stresses, the latter is the compatibility condition when normal stresses are nonzero and there is no overlap but the contact exists. In many cases, the overlap between the surface and the contact boundary is considered as part of the known data determined by the geometry of the bodies.

Mass redistribution method and implicit timing scheme
In the model of dynamic contact interaction we consider numerical realization of contact conditions by using step-by-step analysis of the loading process and changing of contact conditions.
Step-by-step analysis allows to correct the discrete contact boundary depending on the behaviour of interacting bodies, current state of the contact zone, level of the stress-strain state, etc. The step-to-step changing represents the contact nodes inclusion/exclusion conditions which significantly influence the stability of the solution.
One of the solutions for removing unphysical oscillations in the contact region is the use of dissipative integration scheme. The practical results show that using generalized dissipative methods based on the Newmark scheme such as the HTT-scheme which is unconditionally stable, second order of accuracy and dissipative for high frequencies is more effective.
Another approach to prevent false oscillations is associated with the redistribution of the mass belongs to the nodes on the contact boundary area. Thus the resulting inertial characteristics of the mesh nodes change. This method significantly stabilizes the contact domain, in addition it can be used in combination with other schemes. Usually the matrix of masses is chosen to be diagonal and the mass of the cells adjusting to the node is concentrated in the node itself. One of the drawback of scheme (7) is that the nodes on the contact boundary also have their own inertia, which leads to instability even for conservative schemes. During contact phase node loses its kinetic energy and stops. In schemes with energy conservation, the oscillating movement of the contact node provides conservation of the total system energy.
The equilibrium equations of dynamic system of contact bodies without damping according to HTT-method can be written in the matrix form + Δ + (1 + ) +Δ − = (1 + ) +Δ ( +Δ ) − + (1 + ) +Δ − (6) where + Δ is the current time step; , are the mass and stiffness matrices, respectively; is the vector of external loads; is the vector of contact forces; is the vector of nodal displacements; ̈ is the acceleration vector; ∈ [−1/3,0] is a dissipative parameter.
At each integration step + Δ we construct an iterative process for solving a nonlinear system with a right-hand side depending on the contact forces . [ In this paper the system of equations was solved in block form (9), although it is possible to represent and solve it in the form of a system of two equations with respect to the vector of unknown displacements on the boundary which perhaps provides other more economical ways.

Conjugate analysis of contacting bodies by the Schwartz method
As can be seen from (7), after the linearization of the acceleration in the system, there is a nonlinearity due to the presence of contact forces. One of the variants of solving this kind of equation is the application of iterative refinement of the solution by Schwarz alternating method [5,6] for the subdomains of contacting bodies. Construct a version of the Schwartz alternating method, at each iteration of which the approximation is sought, as the result of an independent solution of the equations (2) for each body Ω ( ) with varying boundary conditions [7]. In this case, the displacement ( ) of the body Ω ( ) is computed in each iteration of the Schwartz method in two steps. In the first step, the initial displacement approximation is given starting from the range of expected values for the contact domain on the contact surfaces of the bodies. After solving this problem, the kinematic condition from (5) on the contact surface will be fulfilled, but the computed contact pressures on the surfaces belonging to the interacting bodies are unbalanced. In the second step, contact stresses are corrected and their equality is achieved, but the resulting displacements are not satisfied the condition (5). Further, the displacements are adjusted by the contacting surfaces aligning. The iterative process is performed until a given error is achieved in the kinematic and force conditions (5)  In order to ensure the coherence of inclusions and exclusions of nodes in different sets Λ Γ ( ) , nodes whose relative displacement with the contact surface of the opposite body is greater than the specified allowable value were excluded from the contact domain . If ( ( ) − ( ( ) ) − ( ) ) /( ( ) − ( ) ) > then the node ( ) is excluded from the set. This operation made it possible to reduce the number of contact iterations. In addition, the exception of inadmissible loads obtained at the previous iteration was considered. Similarly, for the body Ω ( ) .
For each subdomain Ω ( ) , systems of equations are constructed and degrees of freedom, connected with inner and contact nodes, are separated. In the iterative algorithm constructing for solving systems, the residual of the displacements ̅ ( ) and the forces ̅ ( ) on contact surfaces are used and alternately applied on iterations for each body and ensure the execution of kinematic and static conditions on the boundary of the bodies. For the kinematic system (even step), the displacement projection residual is defined as The vector of contact forces on the boundary Γ ( ) is computed in the form 2 +1 using the standard procedure of obtaining nodal forces equivalent to the distributed load, construct the component of the vector of the right-hand side 2 +1 ( ) for a system of equations on odd iterations. To solve systems, the conjugate gradient method was applied with diagonal preconditioning and parallelization of computations based on OpenMP and CUDA technologies After the systems solving check the performance of force and kinematic contact conditions. The convergence of such algorithms was considered in the [7].

Numerical results
The presented numerical studies were performed on unstructured meshes with hexagonal cells with different parameters of the cell quality [8] and for several interpolation options between incompatible meshes.   The difference was not more than 10%, which, among other things, is explained by different points of measurement of the reference values of velocities. The deviation of the total energy from the initial kinetic energy did not exceed 3%.

Five spheres
The sequence of pair contacts (figure 3) and the contact of several spheres (figure 4) were considered. In the first case, the distance between the spheres was chosen in such a way that the previous contact terminated at the beginning of the next one. The parameters of the spheres, the meshes and the speed of the first sphere were taken from the previous problem.

Head impact
Modeling of biomechanical processes involves the study of complex shape objects which requires to construct boundaries that describe the inner structure of biological tissues. One of the known problems is the Nahum's head impact experiment (figure 5 a) [10]. To increase the realism of the computational experiment, the finite element model was constructed from computed tomography data and consisted of three main materials (figure 5 b): soft tissue, skull and brain using the algorithm [11].   8 Note that along with the use of standard procedures for optimizing of the mesh cells quality, special functions for adjusting the shape of the material area and removing the "noisy" areas of the real physical model were used: filling internal cavities, searching for hanging cells not connected with the main model, etc. The unstructured mesh contains 1363369 nodes and 1308341 hexagonal cells. The properties of the facial tissues, scalp and skull layer, considered as a single whole, as well as brain tissues, are linearly elastic, the material constants are borrowed from [12][13][14]. As a result of numerical modeling of Nahum's experiment, a maximum of 6882 N was obtained, which agrees well with the experimental results of 6900 N (figure 6 a). The deviation of the total energy from the initial kinetic energy of the impactor did not exceed 2%.

Conclusion
To solve dynamic contact problems, Schwarz's alternating algorithm is provided, which ensures that the conditions on the varying contact boundary are fulfilled with an implicit time integration scheme and the mass redistribution of the contact nodes. This approach does not require special contact finite elements, Lagrange multipliers or regularization. Computational experiments are shown that the applied approaches to the computation of contact forces are provided the same error in computations for a smaller number of iterations than known algorithms. The independence of the solution approximation for each body makes it possible to use efficient numerical schemes and multi-level parallelization.