A generalised path-following solver for robust analysis of material failure

When analysing complex structures with advanced damage or material models, it is important to use a robust solution method in order to trace the full equilibrium path. In light of this, we propose a new path-following solver based on the integral of the rate of dissipation in each material point, for solving problems exhibiting large energy dissipating mechanisms. The method is a generalisation and unification of previously proposed dissipation based path-following solvers, and makes it possible to describe a wider range of dissipation mechanisms, such as large strain plasticity. Furthermore, the proposed method makes it possible to, in a straightforward way, combine the effects from multiple dissipation mechanisms in a simulation. The capabilities of the solver are demonstrated on four numerical examples, from which it can be concluded that the proposed method is both versatile and robust, and can be used in different research domains within computational structural mechanics and material science.

mass-scaling to overcome long simulation times 1 . Implicit transient solvers, on the other hand, may encounter convergence issues near sudden force drops. Although a reduction of the time step often leads to improved convergence due to the increased effect of the mass matrix, such small time steps lead to long simulation times.
As a result, quasi-static simulation techniques are often favourable in instances where the failure event is slow and gradual (e.g. ductile failure) or the inertia and highfrequency oscillations are of secondary importance. In addition, whereas a quasi-static solver can be used to analyse the overall stability of a mechanical system, this is also not true for transient solvers for which potential structural instabilities are missed as the solution passes dynamically to another stable branch [2]. Quasi-static simulations will reveal these unstable points.
Conventional quasi-static solution procedures, where the applied external loads or prescribed displacements are controlled incrementally, includes their own sets of issues. For example, if the problem in question experience snap-back or snap-through behaviour, caused e.g. by unstable crack growth and/or buckling, they will fail to converge. In such cases, important information regarding limit points, equilibrium paths, damage tolerance, or other mechanical responses may be lost. An alternative quasi-static solution scheme that has been shown to provide robustness in cases of material or structural instability, are so-called path-following methods, where the solution is parametrised with a path-following constraint.
The most prominent group of path-following methods within the field of computational mechanics, are the arclength methods. Some of the pioneering works on these methods date back to 1970s by Wempner [3] and Riks [4]. Further improvements have been suggested by Ramm [5] and Crisfield [6] (for an extended historical review of the methods, the reader is referred to Geers [7]). The core idea in these arc-length methods is that the equilibrium path is parametrised with an additional constraint equation, which is expressed in terms of the norm of the incremental vector of displacement degrees of freedom. Such arc-length solvers have proven to be effective in solving problems with geometric non-linearities, such as buckling.
However, the classical arc-length methods encounter difficulties in problems where large local material instabilities are present, such as rapid crack growth or localised plastic deformation [8]. These difficulties can be explained by the observation that only a few degrees of freedom in the model may be controlling the physical instability, resulting in that a global norm fails to accurately accommodate for this locality. As a remedy, so-called local arc-length methods have been investigated, where only a subset of degrees of freedom near the material instability are included in the constraint equation, see e.g. [8] or [9] for application to strain-softening behaviour.
A drawback of these local arc-length methods is that the location of the failure process zone must be known in advance. This in turn has resulted in the investigation of algorithms for monitoring of internal variables in order to smartly select controlling degrees of freedom in the vicinity of the instability, see e.g Alfano et al. for cohesive zones [10], or Geers [7] or Pohl et al. [11] for work on the sub-plane method. To our knowledge there exists no general selection scheme that can accommodate all types of damage mechanisms, and selection algorithms have only been developed for some specific applications.
An alternative approach for tracing the equilibrium path in problems with material instabilities was suggested by Gutiérrez [12]. He derived a path-following constraint equation based on the energy release rate in geometrically linear solids with continuum damage, which was expressed solely in terms of the global displacements and force vectors. Such a dissipation based constraint leads to the advantage that no a priori information about the location of material instabilities is required. Furthermore, the energy release rate also fulfils the requirements needed for the path-following constraint; it is monotonically increasing and it is a function of the unknown degrees of freedom of the system.
The proposed dissipation-based path-following constraint was later extended by Verhoosel et al. [13], who derived equivalent path-following constraint for two additional types of problems, (i) a geometrically linear model with plasticity and (ii) a geometrically non-linear model with damage. As opposed to Gutiérrez, these constraints were not expressed only in terms of the global state and force vector, but included additional terms with new global vector quantities (that were assembled similarly to the internal force vector). The two new constraints were shown to robustly solve problems with significant material instabilities. A drawback of this method is that each type of problem has to be treated as a special case. Thereby, a general approach which accommodates different dissipating mechanisms can not be formulated in a unified way.
In the current contribution, we reformulate and generalise the dissipation based path-following equation proposed in [12,13], by determining the dissipation at material point level instead of being calculated with global quantities. Specifically, the dissipation is calculated in each continuum and/or interface quadrature point and then integrated numerically over the entire domain. This integral formulation has the benefit of unifying the path-following constraint for all types of dissipative mechanisms, which means that no adjustment of the solver is needed to accommodate a specific problem. Furthermore, there are also instances where it is difficult or impossible to formulate the dissipation constraint equation using global quantities, thereby making the computation of the dissipation on material point level the only feasible option.
The remainder of this paper is organised as follows. In the next section we present the general theory for solvers with path-following constraints, followed by details about the proposed approach. In Sect. 3 we show four numerical examples, where the proposed method is demonstrated and evaluated against alternative solutions methods. The numerical examples include (i) buckling driven damage, (ii) phase field damage modelling, (iii) multiple and unstable delamination growth in a laminated fibre-reinforced composite and (iv) plasticity with large deformation kinematics coupled with material interface (cohesive) failure. In Sect. 4 we summarise the paper with and present the main conclusions together with perspectives on future implication and future work.

A path-following constraint based on dissipation in material points
In this section we first summarise the theory related to pathfollowing solvers, and review some common path-following constraints found in the literature. Next, we describe the novel path-following constraint introduced in the work, and outline the corresponding solution procedure.

Preliminaries of path-following solvers
Consider the weak form of the governing equations for a body in quasi-static equilibrium (body forces ignored), see Fig. 1: where is the reference configuration of the body, S is the Second Piola-Kirchhoff stress, P is the first Piola-Kirchhoff stress, E(δu) is the (virtual) Green-Lagrange strain. Furthermore,ū are prescribed displacements on u , andt is the prescribed traction on the boundary t with outward normal n.
In a finite element setting, the displacement field u(X) and the virtual displacements δu(X) can be interpolated by: where N I (X) are vector-valued shape functions, → R dim , a I are displacement degrees of freedom, and N dof is the number of degrees of freedom in the model. Substituting these relations into the Eq. (1) and using that δa I is arbitrary, we obtain: Fig. 1 A body subjected to external loadst and prescribed displacementsū The left and right-hand side of Eq. (3) represent the elements of the internal force vector, f int (a), and the external force vector, f ext , respectively, which can be expressed as the following residual r: where a is the vector of the unknown degrees of freedom a I . Equation (4) is typically solved in either a displacement or force controlled setting together with Newton iterations. In a force controlled setting, the external force vector is typically expressed as f ext = λf +f , where λ is a load parameter, andf is a unit direction vector, andf is a force vector independent of λ.
If the solution of the problem in question involves either snap-back or snap-through behaviour, the displacement or force controlled solution approach will encounter convergence issues, and the full equilibrium path can not be traced. The reason behind this is explained and illustrated for a system with a single degree of freedom in Fig. 2. A possible way of tracing the entire equilibrium path is to make use of path-following methods, where the load parameter λ is reinterpreted as an independent unknown load variable. This leads to the new residual equation: Since the system of equations now has an additional unknown variable (i.e λ), it needs to be complemented with an additional equation, typically called the path-following constraint, ϕ(a, λ) = 0.
Note that ϕ(a, λ) can, in its most general form, be a function of both a and λ. Combining the path-following equation in Eq. (6) with the residual in Eq. (5), the new system of equations takes the form: The augmented system of Eq. (7) can be linearised and solved using a standard predictor-corrector scheme with Newton iterations. The solution at iteration k +1 is therefore obtained as where a is the incremental update of the degrees of freedom between two load steps, and K is the standard tangent stiffness matrix, which is defined as: and h (herein called the dissipation gradient) and w are defined as If one wants to use a displacement controlled solution procedure, the above equations need to be reformulated. First, the displacements are expressed as a combination of free and prescribed displacements [14]: whereâ is a unit direction vector (length N dof with N a prescribed degrees of freedom), a p accounts for boundary conditions independent of the load variable λ (length N dof with N p degrees of freedom with Dirichlet conditions) and a f contains all free degrees of freedoms (length The matrix T is a constraint matrix of size N dof × N f , defining dependence between the free and constrained degrees of freedom. With the absence of external forces, the residual in Eq. (5) can be reformulated in terms of the as: Combining and linearising the path-following constraint ϕ(a, λ) and Eq. (12), the linearised system of equations becomes: In the next section we will describe the path-following constraint proposed herein.

A path-following constraint based on the rate of dissipation in material points
Inspired by the integral formulation of the energy release rate constraint proposed in Verhoosel et al. [13], and later also seen in May et al. [15], we propose an alternative for how the dissipation can be used to control the solution. More specifically, the (specific) rate of dissipation,Ḋ (J/s/m 3 ), is computed directly in each material point, and integrated over the domain: where G is the total dissipation rate (J/s), and is the same domain as in the weak form in Eq. (1). In general,Ḋ will be dependent on the same quantities used to determine the stress-strain response. As such,Ḋ is most easily evaluated in the material routine directly. Furthermore, expressions forḊ will be different for each dissipative mechanisms and material models, but are often readily available if a thermo-dynamically consistent framework has been used. In scenarios where this is not the case, it is likely possible to derive an approximate measure for the rate of dissipation. Explicit expressions forḊ are given for some common models in Sect. 3. By time-discretising Eq. (14) (in this work we have used a backward Euler scheme), the path-following constraint equation ϕ(a) takes the form: where D is the incremental energy dissipation between two load steps, and τ is a parameter constraining the amount of dissipation in the load step. In addition, the path-following method also require the constraint derivatives h, and w. Returning to the introductory example in Eq. (1), the dissipation gradient, h, is obtained via the chain rule: where h I are the entries of h, and where ∂ E/∂a I was introduced in Eq. (3). As noted earlier, D is evaluated using the same quantities as the stress-strain response, so ∂ D/∂ E is most easily computed in the material routine as well. Finally, with the definition of ϕ(a) in Eq. (15), w will always be equal to zero as ϕ(a) is always independent of λ.
There are multiple advantages of using the (volume integrated) material point dissipation measure introduced in Eq. (15). Firstly, Eq. (15) is a unification of the multiple dissipation constraints previously proposed in the literature by Gutiérrez [12] and Verhoosel et al. [13]. Therein, the proposed constraint equations were derived as special cases of different dissipation mechanisms. Thus, with the proposed approach it is no longer necessary to adjust the solver to specific types of problems, since the complexity of describing the dissipation mechanism has been moved to the material point. This fact is also beneficial for implementation in FE-codes, since the same procedure is used for all possible dissipative mechanisms.
A second advantage of the proposed generalised dissipation measure is that it enables the inclusion of damage mechanisms for which it is impossible (or at least very difficult) to derive a global dissipation expression. One such example is the case of geometrically non-linear deformations with plasticity, where it is not possible to express the elastic energy rate for the solid in terms of global vectors (a step which is crucial in the methods proposed by Gutiérrez [12] and Verhoosel et al. [13]). On the other hand, dissipation measures in the material point is readily available for many large-strain materials, and will be further discussed in Sect. 3.4.
A third advantage is that the proposed formulation also provides a straightforward approach to combine multiple dissipation mechanisms. This is automatically achieved in the solver, by simply adding the dissipation measures from the different damage mechanisms (exemplified in Sect. 3.4).

Numerical examples
In this section we will demonstrate the capabilities of the proposed path-following solver in four numerical examples.
The numerical examples include (i) non-linear delamina-tion buckling (with dissapative interfaces) (ii) phase field damage (iii) multiple and unstable delamination growth in a laminated fibre-reinforced composite and (iv) plasticity with large deformation kinematics coupled with material interface (cohesive) failure.
Before examining the numerical examples, a few clarifications regarding the simulations are made.

Initiation of simulations
The simulations are in all four numerical examples initiated with the Riks arc-length method [4]. This is to accommodate for the fact that no dissipation exists at the beginning of each simulations, whereby the dissipation based path-following solver will not work. Note that the choice of initiating the simulation with a Riks method, as opposed to for example a standard force/displacement controlled increments, is made in order to traverse potential geometrical non-linearites occurring before dissipation has developed in the structure. An automatic switch to the dissipation based path-following solver is performed when the total dissipation rate in the simulation exceeds a user-defined limit, τ switch . 2

Automatic adjustment of 1
In order to avoid too large or too small load steps, the path parameter, τ , is automatically adjusted after each load step. To do so, the value of τ is adjusted according to the number of Newton iterations, I , required to fulfil the convergence criteria in the previous load step: where (n − 1) and (n) refers to the previous and next load step, respectively, I opt denotes the desired number of Newton iterations per load step, and α and β are user-defined parameters. In this work, α has been set to 1/2 and β to 1/4 (following the works of [13]). In addition, if the path parameter has been increased to a large value at the vicinity of a sharp or sudden force drop, the Newton corrector iterations will likely fail to converge. We accommodate for this by halving τ until convergence within the load step is achieved, where I F is the number of failed (i.e. non-converged) attempts at a load step.
In "Appendix B", an algorithm outlining the implementation of the solver is presented.

Delamination buckling of a double cantilever beam
In this first numerical example we will compare the performance between the arc-length solver proposed by Crisfield [6] (Crisfield), the solver proposed by Verhoosel et al. in Section 2.1.3 in [13] (Verhoosel), and the generalised dissipation solver proposed in this paper (Current). The problem to be analysed is replicated from [13], and consists of a double cantilever beam subjected to a compression force, see Fig. 3. This load case includes two non-linear failure mechanisms; one mechanism comes from the buckling behaviour, and the other from the progressive delamination of the interface.
The cantilever beam is L = 20 mm long and h = 0.4 mm thick, and has an initial crack length of a 0 = 10 mm. Furthermore, two axial compression forces are applied to the lower and upper half of the beam, together with two small vertical perturbation forces (10 −5 N) to trigger the buckling mode. The beam is discretised with quadratic quadrilateral elements (plane strain assumed), with 400 elements along the beam, and 4 elements through the thickness.
In Sect. 2, the path-following constraint was introduced in the context of large strain continuum. In this numerical example, the dissipation will be calculated in the cohesive zone elements, whereby the path-following constraint ϕ [Eq. (15)], and its gradient h [Eq. (16)], are instead calculated on the cohesive domain cz : where, J is the displacement jump vector along the cohesive interface.
The interface between the upper and lower part of the beam is modelled using a standard bi-linear cohesive zone law [16]. The constitutive law relates the interface cohesive traction, t, to the displacement jump between the bottom and top surface, J, via: where D defines the initial stiffness of the interface. Furthermore, d is the damage variable that depends on the ultimate traction and fracture toughness of the interface (see Turon et al. [16] for more details). The expression for theḊ along the cohesive interface is well-defined 3 (and derived in detail in "Appendix A"): For this example, the fracture toughness in both normal and transverse direction is set to G I = G I I = 0.2 N/mm. Two values for the ultimate traction in the normal direction, σ ult , and tangential direction, τ ult , will be used to assess the robustness of the solvers; σ ult = τ ult = 40 MPa and 75 MPa. Finally, the bulk material of the beam is modelled using a St. Venant elastic material (to accommodate large displacements and rotations), with a Young's modulus E =135 GPa, and Poisson ratio ν = 0.18.
The resulting relation between the applied external tip force and the resulting horizontal tip displacement is shown in Fig. 4, where the marks on each equilibrium curve represents a converged solution step. In practice, the curves from the three solvers fall on top of each other making it impossible to distinguish any significant differences. Therefore, the solution from the Verhoosel solver has been excluded since it performs very similar to the generalised dissipation solver. Furthermore, the total number of load steps from each solver is presented in Table 1.
From the results, it can be seen that in the case when σ ult = τ ult = 40 MPa (Fig. 4a), the Crisfield solver performs well, showing the ability to take large steps throughout the simulation. The two dissipation based solvers also trace the full path, but require more steps in order to complete the simulation, due to buckling being the dominating failure mode.
In the case when σ ult = τ ult = 75 MPa (Fig. 4b), the energy dissipation from the interface is to a larger extent controlling the deformation. This results in that the Crisfield Fig. 4 Resulting force-displacement curve for different values of ultimate traction in the double cantilever beam example. The markers indicate a converged load step. The zoom in (a) illustrates that the arc-length solver is able to take larger step when buckling is the dominating mode of failure, but also that the same solver experience difficulties when the dissipative interface failure is controlling the deformation (zoom in (b)) Generalised Crisfield Tip force [N] (b) t ult = 75MPa solver encounters difficulties when the interface damage starts to propagate, and eventually leads to that the simulation terminates prematurely. On the other hand, both dissipation based solvers are able to robustly trace the full equilibrium path. It can also be noted that the dissipation based solvers take very similar amount of steps in the simulations. The results confirm that the arc-length solver is favourable in problems involving geometric non-linearities, also that it has a problem to converge in problems involving material instabilities. On the contrary, both the dissipation solver proposed by Verhoosel et al. and the generalised dissipation solver are able to trace the path for both cases of σ ult /τ ult , thereby showing their robustness and ability to solve problems for a range of material parameters.

Phase field damage for brittle fracture
Path-following techniques specifically aimed for phase field damage formulations were first introduced by Singh et al. [14]. It was also later shown to perform well by May et al. [15]. We will here demonstrate that our proposed solution framework works well for simulation of phase field damage modelling.
As a preliminary, a brief overview of the phase-field formulation used is given here, including the definition of the energy dissipation needed for the path-following constraint. The governing equations related to phase field modelling can be stated as: u =ū on u ∇d · n = 0 on (22) where d and u are the unknown (scalar) damage phase field and (vector) displacement field, respectively, G c is the critical energy release rate, and l c is a length scale parameter for regularising the solution. Furthermore,t is an applied traction force on the boundary t whilstū are prescribed displacements on u ( t u = ). The (Cauchy) stress field, σ , and strain field, ε, are defined in this case as where 0 is the energy density for the undamaged solid. Furthermore, H is a history field for enforcing damage irreversibility. H can be defined in several ways, but in this work we follow [14] and use: The specific dissipation rate,Ḋ, is computed as: which emanates from the concept of fracture surface area, used in phase field damage modelling [14]. The loadcase considered is depicted in Fig. 5, where the geometry is inspired from [18] (however we use modified boundary conditions). The geometry consists of a rectangular plate with width W = 65 mm, height H = 120 mm, and thickness t = 1 mm, which has a circular cutout with radius 20 mm near the centre of the plate (x = 36.5 mm, y = 51 mm), and two smaller cutouts with radii 5 mm located near the top and bottom (x = 20 mm, y = 20 mm) and (x = 20 mm, y = 100 mm), respectively. Furthermore, the plate has an initial notch with length a 0 = 10 mm on its left side at y = 65 mm (modelled as a small separation in mesh). The inner boundary of the top cutout is controlled by a prescribed displacement in the vertical direction, whilst the inner boundary on the bottom cutout is kept fixed.
The bulk material has a Young's modulus of E = 6 GPa and a Poisson's ratio of ν = 0.22. Furthermore, the energy The resulting force displacement response is shown in Fig. 6. The structure behaves linearly up until the first force drop at around 8 kN, after which a sharp snap-back behaviour can be identified when the damage field starts to form. This process is repeated when a new damage field is initiated on the right side of the hole. Note that due to the very brittle fracture process, small incremental load steps are required in order to trace the full path. However, the size of the load steps are automatically found by the method described in the beginning of this section. Also note that instabilities in the equilibrium curve can be identified when the crack propagates into the hole. The instability is likely caused by a lowering of the dissipation rate and the use of an unstructured mesh closer to the rim of the hole. However, it can be seen that the solver manages to trace this part as well. Finally, the results from this example confirms that the solver can handle very brittle fractures events with significant snap-back behaviour (Fig. 7).

Multiple delamination growth a in curved composite beam
In the following example, we show how the generalised dissipation solver can be applied to problems related to failure prediction and delamination in fibre reinforced polymer composites. Fibre composites often possess brittle interfacial material properties, potentially leading to large snapback behaviour in quasi-static simulations. The post-failure behaviour in these types of structures is highly influenced by the damage propagation, and for a good understanding of the system the full equilibrium path should be traced. Furthermore, standard arc-length solvers are typically not applicable to these types of problem due to local material instabilities in the interfaces. Therefore, a path-following method based The composite structure to be analysed is a beam with a 90 • bend, shown in Fig. 8. This type of structure has previously been tested experimentally in [19], and a similar model setup is replicated here. The radius at the 90 • bend is R = 6.65 mm, while the length of the arms of the beam is L = 25.4 mm. The total thickness of the beam is h = 3.31 mm, and the width is 20 mm. The left vertex (see figure) of the beam is kept fixed, while the right vertex is constrained in the y-direction, with an applied force in the x-direction. Furthermore, the beam is modelled using 151 linear quadrilateral elements (plane strain) running along the length of the beam, and with 2 elements per layer (total 8 through the thickness).
The fibre composite consists of four 0 • layers (main fibre direction running along the length of the beam). Each layer is modelled as a transversely isotropic material with parameters E 11 = 141 GPa, E 22 = 11 GPa, G 12 = 5.8 GPa, G 23 = 5.8 GPa, ν 12 = 0.3, ν 23 = 0.4. The ultimate traction in shear and normal direction are set to τ ult = 24 MPa and σ ult = 39 MPa, which are within the range of ultimate trac-tion reported in [19]. Furthermore, the critical energy release rate in mode I and mode II are set to G I = 0.0906 N/mm and G I I = 0.943 N/mm.
The interfaces between the layers are modelled with an exponential cohesive zone law proposed by Kolluri et al. [20]. The interface traction in the normal, t n , and tangential, t t , directions, are related to the the interface separations J n and J t , via the following constitutive equations: where δ n and δ t are the critical separations in the normal and tangential directions, respectively, and H denotes the Heaviside-function (that takes into account the noninterpenetration condition). The cohesive material model includes four damage variables; two variables describing the damage progression in tension (d n ) and shear (d t ), along with two damage variables coupling the two modes (d cn and d ct ) (for more detail these parameters, see [20]). A measure forḊ for the exponential cohesive material was derived in Auth et al. [21]: The resulting interfacial damages are shown in Fig. 9, and the equilibrium path obtained with the generalised dissipation solver is shown in Fig. 10. It can be seen that the structure behaves elastically up till the first force drop at around 95 N, where the first crack starts to initiate at the centre of the structure (between the second and third layer). After the unloading and propagation of the first crack, a second delamination starts to grow at the bottom layer. A similar process is followed by a final crack forming at the top interface. It can be concluded that the solver is able to robustly trace the progression of the propagation of all the three delaminations. On the contrary, solutions with both the Crisfield and a Newton solver failed to converge after the formation of the first crack.

Micro structure including two competing damage mechanisms
In this example a simplified composite micro structure in the form of a reduced fibre-matrix unit cell will be analysed, see Fig. 11. The problem includes two dissipative mechanisms; the debonding of the fibre-matrix interface, and the Symmetric bc. Fig. 11 The loadcase and mesh for the fibre-matrix micro structure elasto-plastic deformation of the matrix material. It will be shown that both these mechanisms can be accounted for in the solver framework presented herein. The fibre-matrix volume (square) has a side length of L = 0.045 mm, with the radius of the strand being R f = 0.03 mm. Symmetry boundary conditions are used for the bottom and left face, while a prescribed displacement is applied to the right boundary. The top surface is assumed to be traction free.
Large plastic strain is expected in the matrix, and it is therefore modelled using a standard large strain plasticity material model with isotropic hardening (see e.g. works by Simo [22]). This material model assumes a multiplicative split of the deformation gradient, F = F e F p , where F e is the elastic part of deformation gradient and F p is the plastic part. The Second Piola-Kirchhoff stress, S, and the right Cauchy-Green deformation tensor, C, are related via the following constitutive equation: Here,S is the Second Piola-Kirchhoff stress in an intermediate configuration, evaluated from an Neo-Hookian hyper elastic law: where C e = F T e F e is the elastic right Cauchy-Green deformation tensor.
Since the material model is formulated in a thermodynamically consistent framework, an expression forḊ is readily available: For the fibre material, only small deformations are expected, whereby a simple (isotropic) elastic St Venant model is used with the following material parameters: Young's modulus E = 10 GPa and Poisson's ratio ν = 0.35. For the matrix material, the following parameters are used; Young's modulus E = 1 GPa, Poisson's ratio ν = 0.25, yield stress σ y = 65 MPa and hardening stress H = 20 MPa. Finally, the same exponential cohesive zone law as presented in Sect. 3.3, is used for the fibre-matrix interface. The critical energy release of the interface is set to G I I = G I = 0.005 N/mm. Different values for the ultimate traction, t ult = (σ ult = τ ult ), will be used to assess the robustness of the solver (ultimate traction in shear and normal directions is assumed to be equal in all cases). Figure 12 shows the load-displacement curve and deformed mesh when t ult = 70 MPa. For this case, the interface is the limiting factor against failure, represented by interface debonding between the matrix and fibre. It is shown that the generalised dissipation solver is able to find the solution for the complete propagation of the fibre-matrix interface. Furthermore, since the dissipation from the interface is dominating the solution, the Verhoosel dissipation solver is equally applicable for these specific set of parameters (even though it does not incorporate information about the dissipation from the matrix material). Figure 12 shows that the Verhoosel solver traces the same path as the generalised solver in a similar fashion.
Next we investigate the behaviour of the micro structure when the matrix material accounts for the majority of the energy dissipation in the solution. We do this by increasing the ultimate traction to t ult = 95 MPa. Figure 13 shows the load-displacement curve and deformed mesh. As before, the material-point based dissipation solver is able to trace the full path since it incorporates information from both dissipa-tion mechanisms. However, the Verhoosel dissipation solver cannot be used for this set of parameters, since no global expression for the matrix dissipation can be obtained. Also note that for this case the equilibrium path is relatively simple and does not show any snap-back behaviour, whereby a standard Newton solver could also be used to trace the path (not shown in the figure).
For the final investigation, we will analyse what happens if simultaneous energy dissipation from both the fibre-matrix interface and matrix plasticity is present (and competing) in the solution. We do this by choosing the ultimate traction to t ult = 84 MPa. To provide a reference solution, an explicit dynamic (quasi-static) solution has been performed for this case. The densities of the fibre and matrix was assumed to be 2×10 −3 g/mm 3 and 1×10 −3 g/mm 3 , respectively. However, due to the small size of the micro-structure, these densities cased very small time steps, whereby they were scaled by a factor of 10. A quasi-static event was simulated by setting the loading rate to 0.1 mm/ms.  The resulting load-displacement curve is presented in Fig. 14. It is shown that the generalised dissipation solver is able to trace the full path in accordance with the explicit reference solution, while the Verhoosel dissipation solver enters a different (incorrect) branch. This branch matches the solution from when t ult = 95 (Fig. 13), which means that it enters a branch where the interface remains intact. The reason it enters the incorrect branch (in relation to the reference solution), is due to it lacking information about the complete dissipation state of the load step. Finally, it is also worth to note that a standard Newton solver is also not able to trace the correct path as it is not well-equipped to accommodate for snap-back behaviour.

Concluding remarks
In this contribution, we have proposed a new path-following solver for tracing the equilibrium path in quasi-static problems exhibiting energy dissipating mechanisms and potential snap-back or snap-through instabilities. The proposed constraint equation is based on the integrated rate of dissipation in each material point, and allows to include both continuum and interface dissipative mechanisms in a unified way.
An important aspect of path-following solvers based on the total dissipation rate (or energy release rate), is that they are only applicable if the structure being analysed exhibit some energy dissipation. However, this is not always the case, for example at the start of a simulation. In this work, this was resolved by initiating the simulations with a standard arclength solver, and then switch to the dissipation based solver when the dissipation exceeded a pre-defined limit. However, for more complicated structures/problems with alternating elastic and dissipating parts of the equilibrium path, a more general switching algorithm could prove to be more efficient.
In such cases, hybrid methods of dissipation and geometrical path-following constraints can be considered, see e.g Bellora et al. [23].
A consequence of the generalised formulation of the pathfollowing constraint is that the gradient h I in Eq. (16) needs to be updated and assembled each Newton iteration. This is commonly not required in other path-following constraints, where h I only needs to be assembled once every load step. As such, the proposed method requires more computational operations than other dissipation based constraints, which will slightly increase the computational cost compared to other solvers.
Furthermore, the proposed method requires that an expression forḊ of the material is available. This requirement does not exist for the path-following constraints derived in Verhoosel et al, where only the material tangent is required to approximate the dissipation. Since the material tangent is always available in a finite element setting, their method is able to simulate materials where an expression for theḊ is missing.
While the solver presented herein has been shown to be robust, it may still fail to converge at specific load steps where the equilibrium state is difficult to trace (such as sudden force drops). As a remedy, we introduced a scaling (halving) of the path-following parameter (Eq. (18)) which appears to be an effective method to handle such events. Furthermore, it was also found that the solver was sensitive to the parameters α, β (Eq. (17)) and τ switch , which affect the solver's ability obtain convergence. A trial and error approach for choosing these parameters is therefore sometimes needed in order to find the full equilibrium path.
The capabilities of the proposed dissipation based pathfollowing solver is tested in four numerical examples, where it is demonstrated that our novel formulation enables robust solution of a wide range of problems. The result from the numerical examples are summarised below.
The first example involves a double cantilever beam exhibiting both material and geometrical instabilities (interface failure and buckling), where it is demonstrated that a dissipation based path-following solver outperforms traditional arc-length solvers when the material instabilities (in this case interface crack growth) are dominating the problem.
The second example shows that the proposed solver can be directly used for solving phase field damage problems with a monolithic solver scheme. It is further demonstrated that the solver is stable enough to follow the entire load history even with severe snap-back instabilities.
The third example demonstrates that the method can be used to simulate multiple and unstable crack growth in layered fibre reinforced structures. These types of structures often possess brittle interface properties which can lead to local and global snap-back behaviour. There is thus an obvious need for robust solvers to analyse these types of structures, such as the dissipation based path following solver proposed herein.
The last example involves a simplified micro-structure of a fibre-matrix interface where dissipative mechanisms from both the fibre-matrix interface and the matrix plasticity is included. It is demonstrated that our solver can incorporate the the total dissipation rate from both the interface and matrix plasticity, which is required for solving the full parameter space.
As an outlook, we believe that the numerical results in this paper indicate that the path-following method presented herein could be beneficial also when analysing other types of problems where dissipation is present e.g. problems with rate dependent materials or gradient damage models. Thus, further investigation in these areas would be interesting to see in the future.
Funding Open access funding provided by Chalmers University of Technology.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Dissipation measure for bi-linear cohesive zone law
The second law of thermodynamics for a cohesive zone describing interface failure can be expressed as: where t is the traction vector at the interface,J is the time derivative of the interface separation, and is the Helmholtz's free energy associated with the interface deformation and failure. For a standard bi-linear cohesive zone material, can be expressed as [16]: where d is a damage variable that can vary between 0 and 1, and D is a diagonal tensor representing the stiffness of the interface. Inserting Eq. (32) into Eq. (31) we obtain: By definition, t = (1 − d) D · J, whereby leading to that the reduced dissipation inequality becomes: An evolution law for d is presented in Turon [16].

B Proposed algorithm for solver
Algorithm 1 presents a possible algorithmic scheme for how the solver presented herein can be implemented in a finite element code.