PSI - Issue 39

Daniele Amato et al. / Procedia Structural Integrity 39 (2022) 582–598 Author name / Structural Integrity Procedia 00 (2019) 000–000 • Stiffening regime : The contact pressure varies quadratically for penetrations in the range to , while the penalty stiffness increases linearly from to . The default final penalty stiffness, , is equal to 100 times the representative underlying element stiffness. The default value of is 3% of the same characteristic length used to compute (discussed above). • Constant final penalty stiffness regime : The contact pressure varies linearly, with a slope equal to for penetrations greater than .

588

7

Figure 5 Linear penalty pressure-overclosure relationship.

The low initial penalty stiffness typically results in better convergence of the Newton iterations and better robustness, whereas the higher final stiffness keeps the overclosure at an acceptable level as the contact pressure builds up. Advantages of the penalty method include numerical softening and no use of Lagrange multipliers. In general, the penalty enforcement methods sometimes provide more efficient solutions (generally due to reduced calculation costs per iteration and a lower number of overall iterations per analysis) at some (typically small) sacrifice in solution accuracy [32]. The result of the numerical problem, the stress and strain fields are computed for all specimen nodes, in particular for the mid-nodes of the collapsed brick elements. These stresses are used during the post-processor phase to obtain the K-factors. 3.3 Post-processor The post-processor phase aims to determine the Stress Intensity Factors for all crack front nodes and to extend them towards their propagation direction. The procedure for the calculation of the propagation direction and the crack growth rates is explained in detail in the next section. The stress analysis results for the cracked model, coming from the resolution of the numerical problem, are used to compute the SIFs in FRANC3D. The Stress Intensity Factors (SIFs) are computed at mid-side nodes along the crack front by using the M-integral formulation. In the 3-D domain version, the M -integral integration takes place over a volume, according to Eq. (2): ( 1 , 2 ) = ∫ � ( 1 ) ( 2 ) 1 + ( 2 ) ( 1 ) 1 − ( 1 , 2 ) 1 � = 1,2,3 = 1,2 (2) In this equation, is the Kronecker delta and is a function that is 1 at the stress intensity evaluation point and zero on the outer boundary of the domain of integration. To evaluate the M -integral numerically, two solutions are