3. Fundamentals on Mortar Methods for Computational Contact Mechanics 23
3.5. Time integration and global solution scheme
3.5.1. Time integration for computational contact mechanics
No matter which spatial discretization scheme is employed, computational contact mechanics introduce an increased complexity with regard to time discretization than standard solid mechan-ics. This is due to the non-smooth characteristic of contact phenomena, which can be identified by discontinuities of the interface velocities in the event of an impact, see Laursen and Love
3.5. Time integration and global solution scheme [153]. Standard time integration methods, such as the One-Step-θ scheme or the generalized-α scheme, which was introduced in Section 2.2.3, implicitly assume that the time derivatives of the displacement unknownsdare continuous. Consequently, they lack in accuracy for some me-chanical conservation laws when applied to computational contact mechanics. However, in this thesis, the introduced generalized-α scheme is exclusively employed for time integration of the contact problems since non-smooth time integration schemes are beyond the scope of this work.
Thus, the final space and time discretized version of the frictional contact problem reads:
Kmassan+1−αm +Kdampvn+1−αf +fint(dn+1−αf) +fc(dn+1−αf,λn+1
−αf)−fext,n+1−αf =0. (3.57) Herein, the contact term arises in complete analogy to the internal forces, i.e. with linear inter-polation rules. Additionally, all constraints are enforced at the time tn+1. This time integration scheme provides a sufficient level in accuracy and robustness for all tested numerical examples later on. For a detailed discussion about conservation properties for linear and angular momen-tum as well as for energy conservation, the interested reader is referred to Popp [210]. Conserva-tion properties of the algorithm for contact of vertices, edges and surfaces developed in Chapter 4 in combination with the generalized-αscheme are discussed in Section 4.8.3.
3.5.2. Global solution scheme
To consider the additional contact nonlinearities stemming from the inequality constraints for normal direction (3.41) and tangential direction (3.44), the idea of aprimal-dual active set strat-egy (PDASS) is employed. It is well-known from literature on constrained optimization (cf. Qi and Sun [221]). When applying a PDASS, the set of all slave nodes S is split into an active setA for nodes with positive Lagrange multiplier values and an inactive setI for nodes which are currently not in contact. In addition, frictional contact requires the splitting of all active nodes A into a slip set SL and a stick set ST. These introduced sets allow for a reformula-tion of the inequality constraints into equality constraints for the corresponding nodal sets. To solve the fully discretized system of nonlinear equations within each time step, the well-known primal-dual active set strategy is re-interpreted as a semi-smooth Newton method, see Chris-tensen et al. [42] and Hintermüller et al. [109]. For this purpose, nodalnonlinear complementar-ity functions(NCP) are introduced, which are non-smooth but equivalently express the inequality constraints in (3.41) and (3.44).
First, the normal contact constraints in (3.41) are considered, see Figure 3.4. The comple-mentarity function for the normal direction has been discussed in Hüeber and Wohlmuth [115]
and Popp et al. [211] and is defined as
Cn,j(λj,d) =λn,j−max(0, λn,j−cn˜gn,j), cn >0. (3.58) The distinction between the active set A and the inactive set I is implicitly contained in this complementarity function through its two different solution branches for the non-smooth max-function. The normal contact constraints are exactly fulfilled when
Cn,j(λj,d) = 0. (3.59)
Figure 3.4: Nodal complementarity function Cn,j for the normal direction (here with cn = 1).
Normal contact constraints are fulfilled forCn,j = 0. The figure is taken from Farah et al. [74].
The parameter cndoes not affect the accuracy of results, but may influence the convergence be-havior of the semi-smooth Newton method depending on the specific problem setup. The choice of cn has been suggested to be at the order of Young’s modulusE of the involved contacting bodies to obtain optimal convergence, see Hüeber and Wohlmuth [115]. Numerical investiga-tions in Popp et al. [211, 212] have shown only a small influence on the convergence behavior of the semi-smooth Newton method. However, unless otherwise stated, all numerical examples in this thesis are realized withcn = 1.
Second, the frictional sliding constraints in (3.44) are reformulated to again apply a semi-smooth Newton scheme. There is not only one possibility to define a suitable complementarity function for this set of constraints, see for example Alart and Curnier [3] and Christensen et al.
[42]. In the following, the formulation introduced in Hüeber et al. [117] and successfully vali-dated in Gitterle et al. [88] is employed. It is visualized in Figure 3.5 and reads
Cτ,j(λj,d) =max(F(λn,j−cng˜n,j),||λτ,j+ctu˜τ,rel,j||)λτ,j
−Fmax(0, λn,j−cn˜gn,j)(λτ,j+ctu˜τ,rel,j), cn, ct >0. (3.60) Herein, the nodal weighted relative slip increment ˜uτ,rel directly arises by multiplication of the relative tangential velocity ˜vτ,relwith a time increment ∆t. In analogy to the NCP function corresponding to the non-penetration constraints (3.59), the solution of Coulomb’s friction law is equivalently expressed as
Cτ,j(λj,d) = 0. (3.61)
Again, the parameters cn as well as ct do not influence the accuracy of the results, but may control the convergence behavior. Considering (3.60), ct can be interpreted as a parameter that balances the different scales of the tangential part of the Lagrange multiplierλτ and the weighted
3.5. Time integration and global solution scheme
Figure 3.5: Nodal complementarity function Cτ,j for the tangential direction (here withcn = 1 andct = 1). Tangential contact constraints are fulfilled forCτ,j = 0. Note that this is a 1D visualization of the tangential contact constraints and thus valid for 2D contact problems only. The figure is taken from Farah et al. [74].
relative tangential slip incrementu˜τ,rel, see also Hüeber et al. [117]. Consequently, it depends on the current deformation and load. Thus, it is suggested to initially choose ct at the order of cn and adapt it by the ratio of the tangential part of the Lagrange multiplier λτ and the weighted relative tangential slip incrementu˜τ,relof the completed previous time step. Eventually, the global problem formulation to be solved consists of equations (3.39), (3.59) and (3.61). Herein, the nodal setsA,I,SLandST can be updated after each Newton iteration.
3.5.3. Algebraic form
Finally, an algebraic representation of the linearized system to be solved within each semi-smooth Newton step is provided. It is based on the work in Gitterle [87] and Popp [210]. The resulting system of equations is of saddle-point type and looks in an abstract form as follows:
The system of equations in (3.62) is of increased system size compared to classical structural problems, as both displacements and Lagrange multipliers show up as primary unknowns. Again, the solution vector contains increments of discrete displacements ∆d and Lagrange multipli-ers∆λ. Furthermore, the displacement unknowns are distinguished between interior nodes(·)N, slave nodes(·)S and master nodes(·)M. The matrix blocks denoted withKcontain terms from linearization of the internal force vector fint, damping terms Kdamp and mass matrix contribu-tions Kmass. The upper tilde symbol (˜·) indicates additional contributions from the linearized
Lagrange multiplier contact force vector fc. The Cmatrices represent the linearization of the complementarity functions in (3.58) and (3.60). One of the main advantages of dual mortar methods is the node-wise decoupling of the slave side contact virtual work resulting in a diagonal slave mortar matrixD, see the explanation in Section 3.4.1.2. Thus, it allows for a computation-ally efficient elimination of the additional Lagrange multiplier unknowns as done in Farah et al.
[74], Gitterle et al. [88], Popp et al. [211, 212] and Wohlmuth [289]. For this purpose, the third row of (3.62) is utilized to express the Lagrange multiplier increment in terms of displacement increments:
∆λi+1 =−Di,−T(KiSN∆di+1N + ˜KiSM∆di+1M + ˜KiSS∆di+1S +riS). (3.63) Here, the inversion of the diagonal slave side mortar matrix D is of negligible cost and all other operations are computationally cheap vector matrix multiplications. By inserting (3.63) into (3.62), the final condensed system of equations becomes:
KN N KN M KN S
KMN +PTKSN K˜MM+PTK˜SM K˜MS+PTK˜SS
−CλD−TKSN CM−CλD−TK˜SM CS−CλD−TK˜SS
i
∆dN
∆dM
∆dS
i+1
=
−hrN | rM+PTrS | rc−CλD−TrSii,T, (3.64) with the well-known mortar projection operator
P=D−1M. (3.65)
The final system of equations in (3.64) is of constant system size and the only remaining degrees of freedom are the displacement unknowns. The discrete Lagrange multipliers can be easily obtained by performing a recovery step based on (3.63).
The solution of the linear saddle-point type system in (3.62) or the condensed system in (3.64) with iterative solution techniques is beyond the scope of this thesis. Details on algebraic and geometric multigrid methods for this purpose can be found in Brunßen et al. [34], Krause [141], Wiesner [286] and Wohlmuth and Krause [292].