+ Statement of need
+ Quasi-static evolution problems arising in fracture are strongly
+ nonlinear
+ (Marigo,
+ 2023),
+ (Bourdin
+ et al., 2008). They can admit multiple solutions, or none
+ (León
+ Baldelli & Maurini, 2021). This demands both a functional
+ theoretical framework and practical computational tools for real case
+ scenarios. Due to the lack of uniqueness of solutions, it is
+ fundamental to leverage the full variational structure of the problem
+ and investigate solutions up to second order, to detect nucleation of
+ stable modes and transitions of unstable states. The stability of a
+ multiscale system along its nontrivial evolutionary paths in phase
+ space is a key property that is difficult to check: numerically, for
+ real case scenarios with several length scales involved, and
+ analytically, in the infinite-dimensional setting. Despite the concept
+ of unilateral stability is classical in the variational theory of
+ irreversible systems
+ (Mielke
+ & Roubíček, 2015) and the mechanics of fracture
+ (Francfort
+ & Marigo, 1998) (see also Nguyen
+ (2000)),
+ few studies have explored second-order criteria for crack nucleation
+ and evolution. Although sporadic, these studies are significant,
+ including
+ (Pham
+ et al., 2011),
+ (Pham
+ & Marigo, 2013),
+ (Sicsic
+ et al., 2014),
+ (León
+ Baldelli & Maurini, 2021), and
+ (Zolesi
+ & Maurini, 2024). The current literature in computational
+ fracture mechanics predominantly focuses on unilateral first-order
+ criteria, systematically neglecting the exploration of higher-order
+ information for critical points. To the best of our knowledge, no
+ general numerical tools are available to address second-order criteria
+ in evolutionary nonlinear irreversible systems and fracture
+ mechanics.
+ To fill this gap, our nonlinear solvers offer a flexible toolkit
+ for advanced stability analysis of systems which evolve with
+ constraints.
+
+
+ <bold>Functionality</bold>
+ We attack the following abstract problem which encodes a selection
+ principle:
+
+ 0, \text{ find an } \text{ irreversible-constrained evolution } y_t
+ ]]>
+ P(0): Given T>0, find an irreversible-constrained evolution yt
+
+
+ yt:t∈[0,T]↦Xt such that
+
+
+ [Unilateral Stability]E(yt)≤E(yt+z),∀z∈V0×K0+[1]
+ Above,
+
+ T
+ defines a horizon of events. The system is represented by its total
+ energy
+
+ E
+ and
+
+ Xt
+ is the time-dependent space of admissible states. A generic element of
+
+
+ Xt
+ contains a macroscopic field that can be externally driven (or
+ controlled, e.g. via boundary conditions) and an internal field (akin
+ to an internal degree of order). In the applications of fracture, the
+ kinematic variable is a vector-valued displacement
+
+
+ u(x)
+ and the degree of order
+
+ α(x)
+ controls the softening of the material. Irreversibility applies to the
+ internal variable, hence an irreversible-constrained evolution is a
+ mapping parametrised by
+
+ t
+ such that
+
+ αt(x)
+ is non-decreasing with respect to
+
+ t.
+ The kinematic variable is subject to bilateral variations belonging to
+ a linear subset of a Sobolev vector space
+
+
+ V0,
+ whereas the test space for the internal order parameter
+
+
+ K0+
+ only contains positive fields owing to the irreversibility constraint.
+ The main difficulties are to correctly enforce unilateral constraints
+ and to account for the changing nature of the space of variations.
+ HybridSolver (1)
+ BifurcationSolver, (2) and
+ StabilitySolver (3) address the solution of [1]
+ in three stages:
+
+
+ A constrained variational inequality; that is first order
+ necessary conditions for unilateral equilibrium.
+
+
+ A singular variational eigen-problem in a vector space; that is
+ a bifurcation problem indicating uniqueness (or lack thereof) of
+ the evolution path.
+
+
+ A constrained eigen-inequality in a convex cone; originating
+ from a second order eigenvalue problem indicating stabilty of the
+ system (or lack thereof).
+
+
+ These numerical tools can be used to study general evolutionary
+ problems formulated in terms of fully nonlinear functional operators
+ in spaces of high or infinite dimension. In this context, systems can
+ have surprising and complicated behaviours such as symmetry breaking
+ bifurcations, endogenous pattern formation, localisations, and
+ separation of scales. Our solvers can be extended or adapted to a
+ variety of systems described by an energetic principle formulated as
+ in [1].
+
+ Software
+ Our solvers are written in Python and are
+ built on DOLFINx, an expressive and
+ performant parallel distributed computing environment for solving
+ partial differential equations using the finite element method
+ (Baratta
+ et al., 2023). It enables us wrapping high-level functional
+ mathematical constructs with full flexibility and control of the
+ underlying linear algebra backend. We use PETSc
+ (Balay
+ et al., 2023), petsc4py
+ (Dalcin
+ et al., 2011), SLEPc.EPS
+ (Hernandez
+ et al., 2005), and dolfiny
+ (Habera
+ & Zilian, 2024) for parallel scalability.
+ Our solver’s API receives an abstract energy functional, a
+ user-friendly description of the state of the system as a dictionary
+ {"u": u, "alpha": alpha},
+ where the first element is associated to the reversible field and
+ the second to the irreversible component, the associated constraints
+ on the latter, and the solver’s parameters (see an example in the
+ Addendum).
+ Solvers can be instantiated calling
+ solver = {Hybrid,Bifurcation,Stability}Solver(
+ E, # An energy (dolfinx.fem.form)
+ state, # A dictionary of fields describing the system
+ bcs, # A list of boundary conditions
+ [bounds], # A list of bounds (lower and upper) for the order parameter
+ parameters # A dictionary of numerical parameters
+ )
+ where [bounds]=[lower, upper] are required
+ for the HybridSolver. Calling
+ solver.solve(<args>) triggers the
+ solution of the corresponding variational problem. Here,
+ <args> depend on the solver (see the
+ documentation for details).
+ HybridSolver solves a (first order)
+ constrained nonlinear variational inequality, implementing a
+ two-phase hybrid strategy which is ad hoc for
+ energy models typical of applications in damage and fracture
+ mechanics. The first phase (iterative alternate minimisation) is
+ based on a de-facto industry standard, conceived to exploit the
+ (partial, directional) convexity of the underlying mechanical models
+ (Bourdin
+ et al., 2000). Once an approximate-solution enters the
+ attraction set around a critical point, the solver switches to
+ perform a fully nonlinear step solving a block-matrix problem via
+ Newton’s method. This guarantees a precise estimation of the
+ convergence of the first-order nonlinear problem based on the norm
+ of the (constrained) residual.
+ BifurcationSolver is a variational
+ eigenvalue solver which uses SLEPc.EPS to explore the lower part of
+ the spectrum of the Hessian of the energy, automatically computed
+ performing two directional derivatives. Constraints are accounted
+ for by projecting the full Hessian onto the subspace of inactive
+ constraints
+ (Jorge
+ Nocedal, 1999). The relevance of this approach is typical of
+ systems with threshold laws. Thus, the solve
+ method returns a boolean value indicating whether the restricted
+ Hessian is positive definite. Internally, the solver stores the
+ lower part of the operators’ spectrum as an array.
+ StabilitySolver solves a constrained
+ variational eigenvalue inequality in a convex cone, to check whether
+ the (restricted) nonlinear Hessian operator is positive therein.
+ Starting from an initial guess
+
+ z0*,
+ it iteratively computes (eigenvalue, eigenvector) pairs
+
+
+ (λk,zk)
+ converging to a limit
+
+ (λ*,z*)
+ (as
+
+ k→∞),
+ by implementing a simple projection and scaling algorithm
+ (Moreau,
+ 1962),
+ (Pinto
+ da Costa & Seeger, 2010). The positivity of
+
+
+ λ*
+ (the smallest eigenvalue) allows to conclude on the stability of the
+ current state (or lack thereof), hence effectively solving P(0).
+ Notice that, if the current state is unstable
+ (
+
+ λ*<0),
+ the minimal eigenmode indicates the direction of energy
+ decrease.
+ We dedicate a separate contribution to illustrate how the three
+ solvers are algorithmically combined to solve problem P(0) in the
+ case of fracture.
+ [fig:convergence]
+ illustrates the numerical convergence properties of the
+ StabilitySolver in a 1d verification
+ test.
+ In a
+ supplementary
+ document, we perform a thorough verification of the code
+ through parametric benchmark for investigating the stability of a 1D
+ mechanical system, providing analytical expressions used for
+ comparison with numerical solutions, as well as all parameters
+ (numerical and physical) employed in the calculations. Accuracy and
+ reliability of the solvers is shown by the close agreement between
+ numerical and analytic solutions in a benchmark minimisation of (a
+ constrained) Rayleigh ratio, a key problem for applications in
+ structural mechanics and stability analysis.
+
+ Rate of convergence for
+ StabilitySolver in 1d (cf. benchmark
+ problem in the
+ Addendum).
+ Targets are the eigenvalue
+
+ limkλk=:λ*
+ (pink) and the associated eigen-vector
+
+
+ x*
+ (error curve in blue). Note that the residual vector (green) for
+ the cone problem need not be zero at a
+ solution.
+ ![]()
+
+
+
+ Acknowledgements
+ ALB acknowledges the students of MEC647 (Complex Crack
+ Propagation in Brittle Materials) of the
+ Modélisation Multiphysique Multiéchelle des Matériaux et des Structures
+ master program at ENSTA Paris Tech/École Polytechnique for their
+ contributions, motivation, and feedback; Yves Capdeboscq,
+ Jean-Jacques Marigo, Sebastien Neukirch, and Luc Nguyen, for
+ constructive discussions and one key insight that was crucial for
+ this project. The work of PC was supported by the JSPS Innovative
+ Area grant JP21H00102 and JSPS Grant-in-Aid for Scientific Research
+ (C) JP24K06797. PC holds an honorary appointment at La Trobe
+ University and is a member of GNAMPA.
+
+
+
+
+