Rootfinding and Optimization Techniques for Solving Nonlinear Systems of Equations Arising from Cohesive Zone Models

TL;DR

The paper addresses solving highly nonlinear systems arising from cohesive zone models (CZMs) in fracture simulations by systematically comparing root-finding, optimization, and hybrid nonlinear solvers on a 1D finite element testbed. The authors formulate the problem as $\mathbf{r}(\mathbf{u})=0$ with a CZM-induced residual and examine solvers ranging from Picard and Newton to ADAM, BFGS, and dogleg/Steihaug trust-region methods, including exact and approximate Jacobians. Key findings show that fixed-point methods excel for partially damaged (ItP) problems on coarse meshes, while hybrid/trust-region approaches are more robust where fixed-point iterations struggle (ItC); on finer meshes, solver performance becomes highly sensitive to Jacobian treatment, with Broyden variants and Newton-type methods offering different trade-offs. The study highlights the potential of restarts and method switching to leverage region-specific strengths and points to future work on 3D CZMs and XFEM to broaden applicability. The results provide practical guidance for selecting and adapting nonlinear solvers in CZM-based fracture simulations, with implications for larger-scale and more complex materials problems.

Abstract

While approaches to model the progression of fracture have received significant attention, methods to find the solution to the associated nonlinear equations have not. In general, nonlinear solution methods and optimization methods have a rich body of work spanning back to at least the first century, providing the opportunity for advancement in the field of computational discrete damage modeling. In this paper, we explore the performance of established methods when applied to problems involving cohesive zone models to identify promising methods for further improvement in this specialized application. We first use a simple 1D example problem with low degrees of freedom (DoF) to compare nonlinear solution methods, thereby allowing for both straightforward and intuitive visualization of the residual space and reasoning about the cause for each method's performance. We then explore the impact of higher DoF discretizations of the same problem on the performance of the solution methods. Finally, we discuss techniques to improve performance or to overcome limitations of the various methods.

Figures (5)

• Figure 1: Residual norms as a function of iteration number for case ItP (left) and case ItC (right).
• Figure 2: Estimated and exact Jacobian condition number during iterations with Broyden's method for ItC case.
• Figure 3: Residual surfaces for case ItP (left) and case ItC (right).
• Figure 4: Example of instability across iterations of the L-BFGS-B method for the ItC case. The red line is the location of the cohesive zone.
• Figure 5: Residual norms as a function of iteration number for case ItP (left) and case ItC (right) with a finer mesh.