Hybrid of monolithic and staggered solution techniques for the computational analysis of fracture, assessed on fibrous network mechanics

The computational analysis of fiber network fracture is an emerging field with application to paper, rubber-like materials, hydrogels, soft biological tissue, and composites. Fiber networks are often described as probabilistic structures of interacting one-dimensional elements, such as truss-bars and beams. Failure may then be modeled as strong discontinuities in the displacement field that are directly embedded within the structural finite elements. As for other strain-softening materials, the tangent stiffness matrix can be non-positive definite, which diminishes the robustness of the solution of the coupled (monolithic) two-field problem. Its uncoupling, and thus the use of a staggered solution method where the field variables are solved alternatingly, avoids such difficulties and results in a stable, but sub-optimally converging solution method. In the present work, we evaluate the staggered against the monolithic solution approach and assess their computational performance in the analysis of fiber network failure. We then propose a hybrid solution technique that optimizes the performance and robustness of the computational analysis. It represents a matrix regularization technique that retains a positive definite element stiffness matrix while approaching the tangent stiffness matrix of the monolithic problem. Given the problems investigated in this work, the hybrid solution approach is up to 30 times faster than the staggered approach, where its superiority is most pronounced at large loading increments. The approach is general and may also accelerate the computational analysis of other failure problems.


Introduction
The theory of strong discontinuities (e.g., see [1][2][3][4][5][6][7][8][9][10][11]) has made a significant impact on the computational analysis of fracture and extended the functionality of many commercial finite element method (FEM) packages.Strong discontinuities enrich the displacement solution space with jumps across the fracture surfaces.The embedded discontinuity finite element method (ED-FEM or E-FEM) and the extended finite element method (X-FEM) have gained the most popularity (e.g., see [12,13] for comparative studies).The E-FEM and the X-FEM represent strong discontinuities via elemental and nodal enrichments of the displacement solution space, respectively.It allows the approximation of fracture problems using coarse FEM meshes, no re-meshing, and results in mesh-independent simulation results.However, said approaches often require a crack-tracking algorithm (except for recent work [14] on E-FEM without crack tracking) and exhibit difficulties capturing crack branching and coalescence of multiple cracks, difficulties that are avoided with the newly emerged phase-field fracture method (e.g., see [15,16]).Unlike discrete crack approaches, it is a smeared crack approach like damage mechanics [17][18][19][20] and therefore requires a localization limiter or a characteristic length scale parameter.Such measure is difficult to determine in practice [21] and can be difficult to implement upon unstructured meshes [22].In addition, a sufficiently refined mesh is necessary to adequately describe the mechanics of the localization zone, leading to the high computational cost of smeared crack approaches.Thus, it prevents from an efficient analysis of failure in fiber networks where the fracture process zone is to be resolved in the individual fibers.
The tangent stiffness matrix of the two-field problem can be indefinite or negative definite [23][24][25][26].The energy potential to be minimized is therefore non-convex.Uncoupling the problem and using a staggered solution approach that solves for the displacement and the state of damage alternatingly, overcomes this issue and results in a convex alternating minimization problem.It is the preferred approach in phase-field modeling of fracture [27] and available in commercial software [28][29][30].It results in a positive definite tangent stiffness matrix, but does not represent the linearized residuum, and the solution then converges much slower as compared to Newton-Raphson iterations.The efficiency of the consistently linearized coupled (monolithic) problem is therefore sacrificed towards the robustness of the solution technique [31].
The E-FEM has received significant attention in the description of failure in truss-bars, beams, and structures thereof [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46].No crack-tracking algorithm is then needed [47].The E-FEM framework allows static condensation of the additional DOF associated with the fracture kinematics [33] directly at the element level, resulting in an operator-splitting method [34] for the evaluation of these DOF.As with the phase-field method, a staggered solution method may be used, and the continuous and discontinuous fields are minimized alternatively [39].In the present work, we expand these ideas and propose a hybrid method to optimize the performance and robustness of the E-FEM models in the analysis of fiber network failure [39,48].The hybrid solution approach is general and may also be applied to other fracture mechanics problems.

Enhanced finite element formulation
We consider a local Cartesian coordinate system {  ,   ,   } in the description of a 3D Timoshenko beam of the length  and the cross-section .The beam's neutral axis  ∈ [0, ] is aligned with the coordinate along which distributed  as well as concentrated  loads are applied.
where  is the axial strain,   ,   denote the shear strains along the  and  directions,   is the change in the angle of twist   around the beam's neutral axis and   ,   represent bending curvatures.The stress resultants  = [           ]  are conjugate to the generalized strain measures Here, the axial force , the shear forces   ,   along the  and  directions, the torsional moment   , and the bending moments   ,   along the  and  directions have been introduced.
The strong form of the local equilibrium equations [49] reads where   ,   ,   and   ,   ,   denote the components of the distributed forces and moments per unit length, respectively.
For simplicity, we consider a two-node beam element of the length   with the linear shape functions  1 = 1 − /  and  2 = /  , and a single Gauss point in the center is used to integrate the FEM equations.The interpolated generalized displacement field then takes the form where  ̅ 1 and  ̅ 2 are the six-dimensional generalized nodal displacement vectors.In addition, a softening hinge allows for the formation of localized failure in the center of the Timoshenko beam.We therefore enrich where  = [  1  2 ]  represents the element's nodal displacement vector, and denotes the corresponding interpolation matrix.The strain field then takes the form and the matrix interpolates the jump, where  = − 2 / = −1/  and  denotes the 6 × 6 identity matrix.In addition,    =    / is the Dirac delta function centered at   , and represents the standard 6 × 12 strain-displacement interpolation matrix, where  1 = where  = [            ]  represents a distributed load along the beam and  a concentrated load applied at the position   along the beam.In the derivation of (11) 2 , we used the identity ∫ (   −  2 ) = 0, a consequence of having the discontinuity in the middle of the beam,   =   /2.In addition, the second term of (11) 2 vanishes if   = 0 or   =  e , which we apply here.
From the principle of virtual work,  int −  ext = 0, and the independence of the admissible variations , , we obtain the two variational statements [33,34,50] where the condition  = ∫        has been used.The condition (12) 2 represents the local residual at the discontinuity,  =   , and ensures the equilibrium between the traction  acting on the discontinuity and the stress resultant vector  in the bulk material, − +  = .
The implementation of (12) results in the set Incremental formulations are used to express the constitutive relations of the bulk material and the fracture process zone, expressions that close the set of equations (13).We use Δ = Δ bulk ; Δ = Δ (14) to describe the development of the stress resultant  and the traction , where  and  denote the respective tangent constitutive tensors.Given a beam of the cross-section  that is made of a linearelastic material with Young's modulus  and shear modulus , the tangent determines the development of , where  is the shear correction factor,  is the polar moment of inertia, whilst  11 and  22 denote the area moments of inertia.Note that (15) is the simplest set of uncoupled linear elastic constitutive equations and based on the assumption that the beam cross-section possesses appropriate symmetries [51].
With the strain  bulk in the bulk material (8), the increment of the stress resultant and incremental local equilibrium across the discontinuity −Δ + Δ =  are to be enforced.
The linearization of the residual force equations (13) with respect to the unknown displacement Δ and the discontinuous displacement Δ yields the system with the sub-matrices The term  *  is explained in Section 3.
In the implementation of these equations, we distinguish between elastic and failure loading.We follow the framework of inelasticity [23,52] and the decision is based on the introduction of a failure surface Φ, see Section 3. Given elastic loading, or unloading, the increment Δ =  and the system reduces to the standard Timoshenko beam FEM model    Δ =   int −   ext .Different implementations concerning damage loading are discussed in the forthcoming sections.

Monolithic, staggered and hybrid FEM implementation
The monolithic implementation considers the consistently derived finite element stiffness as shown in (17), where off-diagonal terms couple the increment of the displacement Δ and the increment of the jump Δ.As with other embedded approaches, our model allows for the static condensation of Δ directly at the element level, which then results in a displacement-based FEM implementation.The second equation in (17), the internal equilibrium, is then used to substitute Δ through In the limit   → 0,    remains positive definite and invertible [6,33].In general, however, the condition (19) poses a limit on the size   of the finite element in the analysis of strain-softening materials.The substitution of Δ by (19) and the embedded formulation again allows to solve it directly at the element level.The system of FEM equations may then be assembled, and the solution of the global system yields the nodal displacements.Here, determines the corresponding finite element stiffness, and given it represents an inconsistently linearized residuum, which results in poor convergence of the staggered approach.
Towards optimizing performance and robustness of the finite element model, we propose the hybrid definition of the finite element stiffness, where  is a numerical parameter, chosen to ensure a positive definite finite element stiffness matrix  hyb  .It is set according to where because we only allow failure in tension.In (16),  *  is the tangent when Δ   = Δ   = Δ   = Δ   = Δ   = 0.
Towards closing the failure description, the evolution of the internal softening variable is to be linked to the evolution of the jump displacement, and Δ = 0 at  < 0 ; Δ = Δ at  = 0 (29) specifies said correspondence at the cases of elastic and failure loading, respectively.
With (16), the internal equilibrium − +  =  for the only non-trivial loading mode reads − +  = 0, resulting in Here,  bulk =  +  denotes the strain in the bulk material, and  is the jump displacement with the initial condition  = 0 at  =  ̅ .Given   and   from the previous time point, we can iteratively derive  and  at the current time point.We therefore expand (30) towards where  trial = ( +   ) and Δ =  −   denotes a consistency parameter that ensures (28); it eventually corrects the trial state by  corr .The implementation follows the classical concept of computational inelasticity [53], and once Δ is given,  =   + Δ sign() updates the solution.
At failure loading, Φ = 0, the algorithmic consistency parameter therefore reads Δ = Φ trial  −  sign() > 0 (33) and allows us to update the solution.The algorithmic consistency parameter Δ must be positive, a standard stability condition to be enforced to avoid snap-back behavior.Table 1 summarizes the predictor-corrector implementation.

Mesh independence
To illustrate the behavior of the finite element, we consider a 0.1  long cantilever beam that is loaded in tension.Young's modulus  = 1.0 , the cross-section area  = 1.0  2 , and the ultimate tensile force  ̅ = 1.0  further describe the problem.The cantilever beam is discretized with one and ten evenly spaced finite elements, where the ultimate tensile force is lowered by 1% in the most left element.Otherwise, the homogenous stress field of the problem would have resulted in an ambiguous failure pattern, see Section 3 elsewhere [39].Figure 1 presents the force-displacement response of our structural problem for three different fracture energies   , a property also expressed as the area under the curve.As expected from the fracture model, the force-displacement response is independent of the finite element discretization; it is identical between the one-element and ten-element discretization.
Fig. 1: Failure of a cantilever beam under tension using a one-element and ten-element discretization.
Axial force  as a function of the axial displacement  at the free end is shown for three different fracture energies   .The inset illustrates the investigated discretization with the strain localization in the red element.

Tensile test of a fibrous tissue specimen
A planar and random network of interconnected 3D beams describes the mechanics of a fiber network of the densities   = {300; 500; 1000} / 3 and with the fiber properties listed in Table 3. In-plane it covers the area  ×  = 18 × 6  2 and the displacement  0 = 9.0  is prescribed along one edge, while all six DOFs at the opposite edge are fixed.Out-of-plane displacements and rotations are prevented, and all information on how the mesh has been generated is reported elsewhere [55] together with the ANSYS input file [39].A quasi-static failure analysis was computed through the incremental application of the prescribed displacement  0 .
Figure 2 shows the stress-strain response of the fiber networks.Computational results refer to the staggered solutions.As expected, we observe a more dissipative response of the fiber networks with the higher fracture energy of the fibers.As an example, Figure 3 shows the evolution of the jump  for the densest network of fibers with the fracture energy   = 0.1 .The fiber network configurations refer to the points A, B and C in Figure 2. Given said fracture energy, the linear softening law (28) provides the jump   = 0.85  at the state of complete rupture.

Numerical stability
The numerical stability of the 1000 / 3 dense fibrous network made of fibers with the fracture energy of 0.2  was investigated with respect to the size of the displacement increment.The prescribed displacement was applied through 500, 2k, 4k, and 8k steps.In each sub-step a maximum of 500 equilibrium iterations were allowed (NEQIT=500 in ANSYS), before the next step was processed.
Whilst the staggered solution converged for all step sizes, the monolithic failed for some step sizes before the staggered solution.Figure 4 reports the results from the stability analysis and indicates the points of failure of the monolithic solution.Failure is most likely linked to the inability to minimize the related non-convex problem, and a further investigation of ANSYS' internal algorithm was not feasible.
Given the monolithic method did not fail, the monolithic and staggered approaches solve the same equilibrium equations (13) and therefore result in the same stress-stain response of the fibrous network.and staggered (dotted line) solutions are explored, where a prescribed displacement was applied through 500, 2k, 4k, and 8k steps, respectively.Crosses denote the point of termination of the (monolithic) computation.

Tensile test of a notched fibrous tissue specimen
The afore explored specimen geometry (see Section 5.2) is modified, and a sharp notch with an opening angle of 20° is introduced in the center of the tensile specimen.A 1000 / 3 dense fibrous network with the fiber properties listed in Table 3 is considered.Figure 5 shows the stress-strain response of the notched tensile specimen together with the evolution of the jump .Aligned with the previous problem, the monolithic approach is numerically unstable, given the prescribed displacement was applied through 2k and 10k steps, respectively.

Performance and stability of the hybrid solution technique
Among all our simulations, the notched fibrous tissue specimen made of beams with the fracture energy   = 0.1  was computationally most demanding.This case will therefore be considered to explore the performance of the proposed hybrid solution technique (25) and benchmarking it against the monolithic and staggered schemes.Towards the optimization of the computation time of large fiber networks with many DOFs, we limit the maximum number of steps to 500 in our benchmarking exercise.We emphasize that at small displacement increments even the staggered solution method converges within a low number of iterations (2 to 4) and then no gain in performance is to be achieved with the hybrid solution.Table 4 lists the cumulative iterations, the number of all iterations needed to compute the solution until the prescribed displacement  0 = 1.08  (corresponding to the endpoint in Fig. 5) is reached, a measure sensitive to the computational effort to solve the problem.In addition, Figure 8 shows the resulting stress-strain curves and the predicted end-configurations for the hybrid solution technique, similar to the results shown in Fig. 5.  Numerical effort is listed in Table 4 and compared to the monolithic and staggered solution techniques.
Right: Development of failure in the fibrous network as predicted through the prescription of 20, 100, 200 and 500 displacement increments, where the jump at complete rupture is   = 0.85 .Any  >   indicates the fiber has failed.

Summary and Conclusion
We studied the computational analysis of failure in fibrous materials, where the individual fibers are modeled as Timoshenko beams with embedded strong discontinuities.Representative benchmark examples have been used, where the recently proposed staggered solution method [39] has been tested against the monolithic solution strategy.Whilst the staggered approach is numerically robust, it does not use the consistent linearization of the nodal forces and therefore suffers from a poor convergence rate.This is especially the case for large displacement increments.The monolithic approach, in contrary, follows from the consistent linearization but can result in a non-positive definite element stiffness matrix, that then requires the solution of unstable equilibria.It is therefore practically not applicable to solve the benchmark problems studied in this work; it erratically fails, and step-size refinement is not always successful.We therefore proposed a novel hybrid solution method that forms the element stiffness through an adaptive 'mixing' of the stiffness of the monolithic and staggered approaches.It may also be seen as a matrix regularization technique to retain a positive definite element stiffness matrix while approaching the tangent stiffness matrix of the monolithic problem.The hybrid method results in a robust and computational efficient solution technique to explore failure in fibrous materials.
The approach is general and may also accelerate the computational analysis of other failure problems.

)
of non-linear equations at the element level, where   int = ∫       and   ext = ∫       +   (  ) are the internal and external element nodal force vectors, respectively.
the condition for  critical ⇒ det  hyb  = 0 is derived as follows.Excluding the six rigid body motion-related DOFs from the system, the effective stiffness matrices are of the dimension 6 × 6, and the only physical root of det  hyb  = 0 results in an expression for  critical and (26) then guarantees a positive definite finite element stiffness matrix.Although our stiffness matrices are sparsely populated, the direct solution of det  hyb  = 0 requires the eigenvalue analysis of one 12 × 12 matrix at every Gauss point for each solution step; a faster and tailored implementation for beam rupture is discussed in Section 4. 3 Predictor-corrector implementation of beam rupture Aiming at modeling the failure of soft fibers, we limit ourselves to the description of failure under tension.Therefore, only the component    =  of the jump  is allowed to evolve, whilst    =    =    =    =    = 0.The development of the traction according to (14) 2 is then determined by the tangent constitutive tensor component  11 =  < 0, whilst all the other components of  are 0. As a result,  * in (18) takes the form

1 . 4 Implementation
Compute elastic trial force and test for failure loading  trial = ( +   ) Φ trial =  trial − ( ̅ +   ) IF Φ trial ≤ 0 THEN  =  trial ;  =   ;  =   sign()  =  trial + Δ sign( trial )  =   + Δ sign( trial )  =   + Δ END DO END IF The description of beam rupture through failure under tension and the uncoupled constitutive model (15) result in an uncoupling of the condition (25).Consequently, the identification of the stability parameter (26) yields the single scalar equation and avoids then the eigenvalue analysis of the stiffness matrices  mono  and  stagg .Here,  min = ℎ tol /  > 0 determines a minimum stiffness, where the numerical tolerance level ℎ tol is determined by the precision level of the hardware/software realization, whilst   is the element length.To avoid illconditioning towards the development of complete fiber rupture, the corresponding elements in the global stiffness matrix are limited to be larger than the minimum stiffness  min , where again said stiffness threshold has been used.

Fig. 2 :
Fig. 2: Stress-strain response of fiber networks with the density of 300 / 3 (yellow), 500 / 3 (blue) and 1000 / 3 (green).Dotted lines represent results achieved with the staggered solution method.The fibers within the network are modeled as Timoshenko beams with the properties listed in Table 3 and the facture energy   = 0.1  (left column) and   = 0.2  (right column).

Fig. 3 : 2 .
Fig. 3: Illustration of the development of failure in the fibrous network of the density 1000 / 3 (left column) and 300 / 3 (right column).The fibers within the network are modeled as Timoshenko beams with the properties listed in Table 3 and the facture energy   = 0.1  (left column) and   = 0.2  (right column).The jump  [] in the displacement field is shown at points A, B, and C (left column) and D, E, and F (right column), highlighted in Fig. 2. At complete rupture, the jump   = 0.85  (left column) and   = 1.70  (right column).Any  >   indicates the fiber has failed.

Fig. 4 :
Fig. 4: Numerical stability with respect to the size of the displacement increment in the failure analysis of the 1000 / 3 dense fibrous network.The fibers within the network are modeled as Timoshenko beams with the properties listed in Table 3 and the facture energy   = 0.2 .Monolithic (solid line)

Fig. 5 :
Fig. 5: Left: Stress-strain response of fiber network with the density of 1000 / 3 and fibers modeled as Timoshenko beams with the properties listed in Table 3. Solid and dotted lines represent results achieved with the monolithic and the staggered solution method, with the monolithic approach failing for large displacement steps.For illustration purposes, the green-red and the blue-orange transition of the solid lines denote a decrease in the time step size for the monolithic solutions to to the staggered solutions.Right: Development of failure in the fibrous network with the jump   = 0.85  at complete rupture.Any  >   indicates the fiber has failed.

Fig. 6 :
Fig. 6: Left: Stress-strain response of fiber networks with the density of 1000 / 3 and fibers modeled as Timoshenko beams with the properties listed in Table 3 and the facture energy   = 0.1 .The curves present results from the hybrid solution technique through the application of 20, 100, 200 and 500 (denoted by a light-to-dark color transition) displacement increments, respectively.
]  with the subscripts ,  and  denoting the respective displacement and rotation components.The generalized strain measures of the 3D Timoshenko beam then read The displacements () and the rotations () are collectively represented by the generalized displacement vector  ̅ = [ denotes the Heaviside step function centered at the middle of the element   .With the boundary conditions  1 =  ̅ 1 and  2 =  ̅ 2 + , where  1 and  2 are the total nodal displacement vectors, the Towards reinforcing the robustness of the model, we may uncouple ,  and solve the problem in a staggered way.The internal equilibrium equation (13) 2 is then solved at the nodal displacement (17)7)results in the system of equations from the previous solution.It explicitly reads =  ′ ;   int ′ =     int ;   ′ =     . and the local {  ,   ,   } systems, respectively.

Table 2
summarizes the implementation in the commercial software package ANSYS through the user-Transform to global coordinates   ′ =      and   int ′ =     int 10.Return , ,   ′ ,   int ′ Benchmark exercises ,   ,   ,   ,   as  =  5. Element degrees of freedom (DOF) deletion IF (  >   ) THEN

Table 4 :
Cumulative iterations (ANSYS variable CUM ITER) to analyze the failure of the notched fibrous tissue specimen formed by beams with the fracture energy   = 0.1 .Failure to derive a computational result is indicated by 'f' in the table.