A threescale domain decomposition method for the 3D analysis of debonding in laminates.
The prediction of the quasistatic response of industrial laminate structures requires to use fine descriptions of the material, especially when debonding is involved. Even when modeled at the mesoscale, the computation of these structures results in very large numerical problems. In this paper, the exact mesoscale solution is sought using parallel iterative solvers. The LaTInbased mixed domain decomposition method makes it very easy to handle the complex description of the structure; moreover the provided multiscale features enable us to deal with numerical difficulties at their natural scale; we present the various enhancements we developed to ensure the scalability of the method. An extension of the method designed to handle instabilities is also presented.
1 Introduction
Since the early 1980’s, a very large number of studies has been conducted on the prediction of debonding in composite laminates, resulting in better understanding of the failure processes of composites. As a result, micromechanicsbased models have been shown to enable accurate predictions of the debonding phenomenon.
The industrialists’ wish to replace expensive experiments by virtual tests for the design of their composite structures has brought about a new issue. Currently, the analysis of industrial problems using previously referenced micromodels is infeasible because the memory size and computing time requirements are far too large. Therefore, most applications to predict the initiation and propagation of debonding in laminates rely on mesoscale modeling [7, 18] which are also rooted in the analysis of micromechanics phenomena. In this paper we retain the model described in [2, 3]; it describes the behavior of the laminates distinctly in the plies, 3D entities, and in the interfaces, 2D entities, the debonding ability being localized in the interfaces and handled through a cohesive behavior.
Even for small test cases, the numerical problem resulting from the mesomodeling of laminate structures is huge. Nevertheless, the latest advances in domain decomposition and multiscale methods provide powerful tools enabling the calculation of laminate structures of reasonable size. We can distinguish two families of solvers able to handle such large numerical problems. The first one consists in using a nonlinear homogenized strategy [29, 11, 10] coupled with local enrichment methods [23, 16, 25, 24, 13, 26]. In this paper we focus on the prediction of debonding in process zones within the structure of the laminate, in which we consider that the debonding processes leading to final failure can be circumscribed. Though, within these potentially large process zones, no low gradient zone can be identified a priori. Thus, the enrichmentbased strategy might be difficult to use. Instead, we wish to use the second family of solvers for large numerical problems, which consists in computing the exact mesosolution everywhere in the process zones, using parallel iterative solvers [9, 22, 17, 14, 12, 8]. Coupling the solution with a reduced model such as a plate model [31, 4] will be the subject of further work, and is not dealt with in this paper.
The mixed domain decomposition strategy described in [21] uses an original concept which consists in splitting the structure in volume substructures separated by surface interfaces. Consequently, the reference problem resulting from the chosen mesomodeling is very naturally substructured, the cohesive interfaces of the model being handled within the interfaces of the domain decomposition method. This idea is developed in Section (2). Furthermore, the resolution of the substructured problem by a LaTIn iterative solver exhibits very interesting numerical properties: the nonlinearities are dealt with through local problems and very few reassembling steps are required. The incremental micromacro LaTIn algorithm is presented in Section (2.2) without any improvement. As shown in Section (2.3) several numerical difficulties are encountered when directly applying the method. The core of the paper (Sections 3 to 6) presents and assesses improvements and adaptations of the method to efficiently handle delamination computations.
The substructured LaTIn method is parametrized by two interface search direction operators. Optimal values have been determined and practical values have been assessed for many mechanical problems (perfect interfaces, contact with or without friction), though for cohesive interfaces an adaptation is required in order the method to be efficient (and in some case feasible), as explained in Section (3).
The method is granted a multiscale nature by the separation of a macroscopic part and a microscopic part of the interface fields. This separation, coupled with continuity conditions, leads to the construction of an automatic homogenization procedure. This concept has been successfully applied in the cases of crack propagation problems in homogeneous media in [15], of composite structures modeled at the microscale [19], and of multiscale problems in time and space [20]. Though, in the case of singularities resulting from the crack tips being localized on the interfaces of the domain decomposition strategy, the homogeneous solution is too poor to represent accurately the solution. It results in a loss of extensibility of the strategy. Therefore, we enhance the method with a specific technique for the calculation of the quantities in the process zone with increased accuracy, which results in a significant improvement of the convergence rate, which is the topic of Section (4).
A consequence to the substructuring naturally introduced to solve the reference problem is that the number of substructures required to solve large delamination problems becomes huge. Hence, the macroscopic problem can not be solved using direct solvers. The introduction of a thirdlevel problem is required to quickly propagate the veryhighwavelength part of the solution. This is achieved by solving the secondscale problem using the balancing domain decomposition method described in [22], as explained in Section (5).
When trying to predict the very final residual strength of the structure, which is of the industrialist’s interest, one has to deal with instabilities and limitpoints problems resulting from the local softening behavior. Hence, an adaptation of the threescale domain decomposition method to arclengthtype algorithms with local control of the loading amplitude has been developed, which is the subject of Section (6).
2 The reference problem and strategy
2.1 Reference problem and substructuring
The laminate structure occupying the domain is made out of adjacent plies occupying Domain , separated by cohesive interfaces (see Figure (1)). An external traction field (respectively a displacement field ) is applied to the structure on Part (respectively ) of the boundary of Domain . The normal to the boundary of Ply P, external to , is . The volume force is denoted . Let be the displacement field, the Cauchy stress tensor and the symmetric part of the displacement gradient in Ply .
The simulation is performed under the assumption of small perturbations and the evolution over time is considered to be quasistatic and isothermal. The problem is solved using an implicit time integration scheme. At each time step of the analysis, the reference equilibrium problem is:
Find , where , which verifies the following equations:

Kinematic admissibility on :
(1) 
Global equilibrium of Structure :
(2) 
Linear orthotropic behavior of the plies:
(3) 
Constitutive equation of the interfaces:
(4)
The gap of displacement of Interface such that has arbitrarily been introduced as . The operator establishes a relation between the primal interface unknown , and the dual interface unknown , which reads :
(5) 
The expression of the local stiffness operator of Interface can be made explicit in the basis (see Figure (2)):
is here the positive indicator function.
The local damage variables are introduced into the interface model in order to simulate its softening behavior when the structure is loaded. Their values range from (healthy interface point) to (completely damaged interface point). The parameters are related to the local energy release rates of the interface’s degradation modes. Denoting the surface strain energy of the cohesive interface,
(6) 
is here the surface strain energy of the cohesive interface.
We assume that the damage variables are functions of a single quantity: the maximum over time of a combination of the energy release rates , :
(7) 
Thus, the evolution laws are:
(8) 
This model has the advantage of using a single damage variable to handle several macroscopic delamination modes of the interface (traction along and shear along and ). However, when setting Parameters and to identified physical values such that , the energy dissipated due to the propagation of the crack is different for the three modes.
2.1.1 Substructured formulation of the reference problem
The laminates structure is decomposed into substructures and interfaces as represented in Figure (3). Each of these mechanical entities possesses its own kinematic and static unknown fields linked by its behavior. The substructuring is driven by the will to match domain decomposition interfaces with material cohesive interfaces, so that each substructure belongs to a unique ply and has a constant linear behavior. A substructure defined in Domain is connected to an adjacent substructure through an interface (Figure (4)). The surface entity applies force distributions , as well as displacement distributions , to and respectively. We write .
On a substructure such that , the boundary condition is applied through a boundary interface .
Let be the Cauchy stress tensor and the symmetric part of the displacement gradient in substructure .
Then, the substructured quasistatic problem consists in finding at a given step of the time integration scheme, where , which verifies the following equations:

Kinematic admissibility of Substructure :
(9) 
Static admissibility of Substructure :
(10) 
Linear orthotropic behavior of Substructure :
(11) 
Behavior of the interfaces :
(12) 
Behavior of the interface at the boundary:
(13)
The formal relation and named ”interface behavior” can be made explicit in the two cases that we have to handle:

Perfect interface:
(14) 
Cohesive interface:
(15) where Substructure (respectively ) belongs to Ply (respectively ), such that .
In the same manner, and for instance in the case of a prescribed displacement boudary interface , the formal relation reads:
(16) 
2.2 Twoscale iterative resolution of the substructured problem
2.2.1 Introduction of the macroscopic scale
In order the scalability of the method to be ensured, a global coarse grid problem is solved at each iteration of the solver. The definition of the macroscopic fields required to construct this problem is done on the interface only.
At each interface , the interface fields are split into a macro part and a micro part , the former belonging to a smalldimension subspace (e.g. 4 macro degrees of freedom per interface in 2D, 9 in 3D).
(17) 
Given the macrospaces and on Interface , the unicity of the decomposition of the interface fields in macro and micro data is ensured by uncoupling the interface virtual works:
(18) 
Usually, a common macro basis for both the kinematic and static interface macro fields is chosen. Numerical tests showed that in order to ensure the numerical scalability of the method the macro basis should extract at least the linear part of the interface forces (see Figures (5) and (7)). Indeed, this macro space contains the part of the interface fields with the highest wavelength. Consequently, according to the SaintVenant principle, the micro complement (which, as explained in next subsection, is found iteratively through the resolution of local problems) has only local influence.
2.2.2 The iterative algorithm
The iterative LaTIn algorithm, which enables the resolution of nonlinear problems, is here applied to the resolution of the substructured reference problem with nonlinearities localized in the interfaces.
The equations of the problem are split into two groups:

linear equations in substructure and macroscopic interface variables:

static admissibility of the substructures

kinematic admissibility of the substructures

linear behavior of the substructures

equilibrium of the macro interface forces


local equations in interface variables:

interface behavior

The interface solutions to the first set of equations belong to Space , while the interface solutions to the second set of equations belong to . The converged interface solution is such that:
(19) 
The resolution scheme consists in seeking the interface solution alternatively in these two spaces: first, one finds a solution in , then a solution in . In order for the two problems to be wellposed, search directions and linking the solutions and through the iterative process are introduced (see Figure (6)).
Hence, an iteration of the resolution scheme consists of two stages:

a local stage:
(20) 
a linear stage:
(21)
In the following sections, the subscript will be dropped.
Local stage
In the local stage, local problems are solved at each point of the interfaces :
(22) 
the two last equations of this system being the search direction . In the case of a cohesive interface, Problem (22) is nonlinear, and solved by a modified NewtonRaphson scheme.
Local linear problems are also solved at each point of the boundary interfaces :
(23) 
Linear stage
The linear stage consists in solving linear systems on each substructure under the constraint of macroscopic equilibrium of the interface forces.
(24) 
Macroscopic admissibility of displacements could also be enforced. In the case of perfect interfaces, it would be easy to derive. In the case of nonhomogeneous or nonlinear behavior at the interfaces, the macro condition would not be practically (and sometimes theoretically) feasible.
Condition (24) is incompatible with the monoscale search direction coupling the interface displacement and forces fields at the linear stage, which reads:
(25) 
Hence this search direction is weakened and verified at best under the macroscopic constraint. Technically this is realized using a Lagrangian whose stationarity leads to a modified local search direction:
(26) 
and to an equation expressing the continuity of the macroforces:
(27) 
where is unique for Interface and set to zero on .
A way to solve this set of equations consists in introducing a relation coupling and into (27). This relation is derived from the local equilibrium of each substructure (10) and from the local modified search direction (26). The problem to be solved for Substructure becomes:
(28) 
where .
Equation (28) is linear. Therefore, one can write a linear relation between the interface displacements and the loading:
(29) 
Operator is the dual Schur complement of Substructure modified by the search direction, while results from the condensation of the volumic loading on interface .
The corresponding interface forces are obtained through the modified search direction (26) and projected onto the macro space:
(30) 
where
is classically viewed as a homogenized behavior of Substructure and is calculated explicitly for each substructure by solving local subproblems (28) taking the vectors of the macro basis as boundary conditions on .
This relation is finally introduced into the equation expressing the admissibility of the macroforces (27), leading to the socalled macroproblem:
(31) 
The macroproblem is discrete by nature. Hence, its an algebraic form , where is the vector of the components of the Lagrange multiplier in the macro basis, is also used in the following.
The righthand side of Equation (31) can be interpreted as a macroscopic static residual obtained from the calculation of a singlescale linear stage. In order to derive this term, the problem (28) must be solved independently on each substructure, setting to zero. The resolution of the macroscopic problem (31) leads to the global knowledge of Lagrange multiplier , which is finally used as prescribed displacement to solve the substructure independent problems (28).
In order to perform the resolutions of (28) on the substructures, finite element method is used. Since the behavior of the substructures is linear, the stiffness operator of each substructure can be factorized once at the beginning of the calculation and reused without updating throughout the analysis, which gives the method high numerical performance.
Algorithm (1) sums up the iterative procedure which has been described in this section.
2.3 Delamination analysis example
A first example of quasistatic delamination analysis is shown in Figure (7). The problem consists in a double cantilever beam (DCB) case. The loading leading to mode I quasistatic crack’s propagation is increased linearly over 10 time steps. The first three of them correspond to the initiation of the delamination and the remainder to the crack’s propagation.
This assessment is realized with a C++ implementation of the mixed domain decomposition method capable of handling the quasistatic analysis of 3D nonlinear problems. The parallel computations use the MPI library to exchange data among several processors.
Each processor is assigned to a set of connected substructures (along with their interfaces); it calculates the associated operators and solves the local problems. This tends to reduce the number of interfaces duplicated among several processors (Figure (17)) and is achieved technically through a METIS routine.
This very simple test case already exhibits numerical difficulties:

The convergence rate of the LaTInbased strategy is highly dependant on the search direction parameters. In the case of cohesive interfaces the iterative solver can even stagnate when using too small values for the search direction parameters. Hence, we describe in next section the way we set them in order to ensure convergence.

The method loses its numerical scalability when the crack’s tip propagates. This phenomenon appears clearly in Figure (15) (”No sub resolution” label). When the delamination process propagates (Time steps 4 to 10), the number of LaTIn iterations to convergence becomes very large. A solution to this problem is developed in Section (4)
3 Analysis of the iterative algorithm parameters
The search direction parameters and are introduced as positive definite symmetric operators. It has been empirically shown in previous studies that an optimal set of these operators exists. Though, these optimal operators are known to be difficult to interpret theoretically, especially when using complex interface behaviors. In addition, they can be expensive to compute even in simplified cases (perfect interface behavior). Our goal is here to find an efficient local approximation of the search direction operators for debonding analysis.
3.1 Search direction
Using too small a value for Parameters can lead to the stagnation or divergence of the algorithm. Actually, in this case, the interface solution is seeked at the local stage in a truncated space , the solutions in the softening part of the local cohesive behaviour being unreachable. Figure (8) illustrates schematically this idea. The constitutive law of a cohesive interface and the search direction equation are projected on the normal direction to the interface, at a given integration point. The solution seeked at the local stage is the intersection of the two curves obtained. This solution is computed by a Modified NewtonRaphson scheme. It clearly appears here that the part of the constitutive law within the dotted frame cannot be reached when using this iterative procedure. Thus, the global minimum of the problem may not be found through the LaTIn iterations.
We thus choose to set to infinite values (see Figure 9) and focus on the choice of Search direction only. Although slightly improved convergence rate could be obtained using classical conjugate search directions and , the time costs of the local stage drop as the local problems can be solved directly.
3.2 Search direction : interpretation
The search direction must be separated into a macro part and a micro part in order to be interpreted. Equation (26) can be rewritten:
(32) 
(33) 
where is the space orthogonal to with respect to the inner product.
Perfect interfaces
Previous studies [17, 21] have focused on the choice of the micro search direction parameter . Classically, when dealing with perfect interfaces, the optimal parameter on interface can be linked to the Schur complement of the structure occupying the domain . As the introduced microscopic interface quantities have a local influence, this Schur complement can be calculated in the vicinity of Interface . In the case of isotropic and homogeneous materials, The scalar has been shown to be a good approximation of this local operator [21], being the Young’s modulus of the adjacent substructure and a characteristic length of the interface. As explained in paragraphs, computing this microscopic optimal search direction parameter is not of primal interest in our case.
The search direction parameter can be interpreted by considering an interface separating two adjacent substructures and . In order to simplify the explanations, we consider here a unique search direction parameters on for both substructures. Taking into account the macroscopic equilibrium of the interface forces along the interface , on can derive from (32) :
(34) 
This equation means that is a macroscopic stiffness of interface . In the case of perfect interface, setting to the infinity equals to enforcing the macroscopic displacement continuity through the interface. As the the macroscopic equilibrium of the interface forces is ensured by equation (27), this choice of the search direction parameter leads to the complete enforcement of the interface macroscopic behavior at the linear stage :
(35) 
Though, as the substructures are small compared to the size of the structure, the characteristic length introduced to approximate the microscopic optimal search direction parameter is very small. Hence, setting to is close to ensuring the macroscopic displacement continuity through perfect interfaces. Moreover, as it will be explained further on, the choice of a unique search direction paramater increases the simplicity and efficienty of the iterative strategy.
In the case of more complex interface behaviors, like softening, ignoring the influence of the macroscopic stiffness can lead to non physical solutions or even to the divergence of the iterative process. This idea is illlustrated by the results in Figure (10). The test case presented is a fourplies composite beam, the black parts of the cohesive interfaces corresponding to initial delamination. The third interface is artificially weakened. The solution labelled ”Reference solution” is obtained by a modified Newton algorithm, while the two following ones result from a LaTIn computation, two different set of search direction parameters being used.
Thus, an effort must be made to correctly set the interface conditions prescribed at the linear stage. In the case of simple interface behaviors (including the perfect interface behavior studied in the previous paragraph), a ”physical” macroscopic interface condition can be prescribed beetween substructures at the linear stage. We will derive the general case from the analysis of these specific studies.
Homogeneous isotropic elastic interface
Let us set the macroscopic search direction parameter of an elastic interface to , where is the local scalar stiffness of the interface. Equation (34) now reads :
(36) 
Hence, setting to the interface stiffness equals to prescribing the macroscopic interface behavior as an interface condition at the linear stage.
Delaminated interfaces under traction
In the case of a completely damaged cohesive interface loaded with traction (so that the gap of interface displacements is positive), the converged interface force fields are null. Thus, both micro and macro search direction parameters should optimaly be set to zero. Indeed, on interface as these quantities result from the local stage and verify the interface behavior. Consequently, setting leads to the enforcement of the relation (equations (32) and (33)).
General case
We can conclude from these three simplified cases that, as the macroscopic equilibrium of the interface forces is enforced at the linear stage through the choice of Admissibility space , one can ensure the verification of the complete linearized interface behavior by the macroscopic interface fields at this stage by setting the interface parameters to specific values. Thus, together with the choice of a rigid search direction, this choice leads to a hybrid iterative strategy :

a modified NewtonRaphson scheme on the macropart of the solution

a Latintype algorithm on the micro part, the search directions and being nonconjugated (which is not classical)
Though, when dealing with cohesive interfaces with nonconstant stiffness, with interfaces partially delaminated, or with delaminated interfaces loaded with both local traction and local compression, the interface behavior that should be verified by the macroscopic fields is not explicit. Yet a too coarse approximation of this behavior might drive the algorithm towards a local minimum of potential energy instead of a ”more physical” global minimum.
Moreover, the mechanical interpretation of the search direction parameters requires to introduce a micro search direction and a macro search direction . But when using this search direction separation, we found out that the CPU time increased, because of the large amount of projections of the interface fields in the macrospace required through the iterative process.
A practical way to choose a common micro and macro search direction parameter that ensures a good interface condition is the subject of next subsection.
3.3 Search direction : practical choice
The search direction parameters are set with respect to the interface behavior, as explained bellow:

perfect interfaces: are set to the optimal micro values described previsouly. As the characteristic length of the interface is very small, these values are high enough to ensure the continuity of the macrodisplacement through the perfect interfaces.

undamaged or partially damaged cohesive interfaces: are set to the optimal macro values for initially undamaged interfaces. In the case of an interface , we thus choose:
(37) even though it cannot be shown theoritically, this strategy has always led to the convergence of the iterative algorithm. Moreover, as the interface conditions between substructures are connected to mechanical properties, no branching to nonphysical solutions has been observed in our cases.

delaminated interfaces: are set to the same values used for undamaged interfaces, unless the contact status of the interface is known (i.e.: unless all integration points are in compression, or unless all integration points are in traction). In the last cases, the search direction parameters are updated as follows :

delaminated interface under compression
