Aurélien Doitrand
email: aurelien.doitrand@grenoble-inp.fr

Eric Martin
email: martin@lcts.u-bordeaux.fr

Dominique Leguillon
email: dominique.leguillon@upmc.fr

Numerical implementation of the coupled criterion: matched asymptotic and full finite element approaches

Keywords: Coupled criterion, Crack initiation, Brittle fracture, Finite element method

come

Numerical implementation of the coupled criterion: matched asymptotic and full finite element approaches

Introduction

Among various methods available for brittle fracture modeling such as, e.g., cohesive zone models [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Tvergaard | The relation between crack growth resistance and fracture process parameters in elastic-plastic solids[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF], models based on damage mechanics [START_REF] Giry | Stress-based nonlocal damage model[END_REF][START_REF] Pijaudier-Cabot | Nonlocal Damage Theory[END_REF] or phase-field models for fracture [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF][START_REF] Molnár | 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF], the coupled criterion (CC) [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF][START_REF] Carpinteri | A finite fracture mechanics approach to structures with sharp V-notches[END_REF] in the framework of Finite Fracture Mechanics (FFM) [START_REF] Hashin | Finite thermoelastic fracture criterion with application to laminate cracking analysis[END_REF][START_REF] Nairn | Exact and variational theorems for fracture mechanics of composites with residual stresses, traction-loaded cracks, and imperfect interfaces[END_REF][START_REF] Taylor | The fracture mechanics of finite crack extension[END_REF] proved its ability and efficiency to predict crack initiation in a broad range of materials and configurations. A detailed review of the CC applications was presented by Weissgraeber et al. [START_REF] Weißgraeber | A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers[END_REF]. The principle of the CC is the simultaneous fulfillment of a stress and an energy condition for crack initiation prediction. Coupling both conditions allows determining the initiation loading and crack length. A main advantage of the CC is that it basically requires only two fracture parameters, namely the material strength and fracture toughness which are physical quantities that can be determined experimentally.

The CC can be implemented through two different methods. The first one, initially proposed by Leguillon [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF], relies on the matched asymptotic (MA) expansion procedure proposed by Leguillon and Sanchez-Palencia [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF]. The MA approach consists in focusing only on a domain around the initiation crack location so as to compute the energy and the stress conditions without considering the whole specimen or structure. The boundary conditions imposed to the domain around the crack location are asymptotic displacement or stress fields corresponding to the stress concentrator under investigation (e.g. sharp or blunted crack, V-notch, hole).

Therefore, it relies on the assumption that the initiation crack location is far enough from the structure boundary or other features that may influence these displacement/stress fields. In the second approach, referred to as the full finite element (FFE) approach, the whole structure under consideration is modeled.

Compared to the FFE approach, the MA one is usually more computationally efficient since it requires less calculations. For instance, in the case of a sharp V-notch loaded in mode I, only one finite element (FE) calculation is needed [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF][START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavities size effects[END_REF][START_REF] Yosibash | A failure criterion for brittle elastic materials under mixed-mode loading[END_REF] whereas several calculations are necessary for the FFE approach. However, the asymptotic approach is only valid if the crack length remains small compared to the characteristic structure dimensions [START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF]. Moreover, the MA approach allows computing the initiation configuration but is unable to deal with unstable crack propagation that possibly occurs after initiation, which can be accounted for using the FFE approach [START_REF] Weißgraeber | A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers[END_REF][START_REF] Doitrand | Mesoscale analysis of damage growth in woven composites[END_REF][START_REF] Doitrand | Comparison between cohesive zone and coupled criterion modeling of crack initiation in rhombus hole specimens under quasi-static compression[END_REF][START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF][START_REF] Martin | Energetic conditions for interfacial failure in the vicinity of a matrix crack in brittle matrix composites[END_REF][START_REF] Sapora | Crack onset and propagation stability from a circular hole under biaxial loading[END_REF]. Some authors compared the MA and FFE implementation of the CC. For instance Yosibash et al. [START_REF] Yosibash | A failure criterion for brittle elastic materials under mixed-mode loading[END_REF] showed that the potential energy variations computed by both approaches for a sharp V-notch were close to each other (difference less than 1%) for small enough cracks, and that this difference was larger for larger cracks. Priel et al. [START_REF] Priel | Failure initiation at a blunt V-notch tip under mixed mode loading[END_REF] obtained similar results in the case of a blunted V-notch. Martin et al. [START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF] showed that the range of validity of the asymptotic approach may be determined by ensuring that the material characteristic length L mat = EGc (where E is the material Young's modulus, G c the fracture energy and σ c the tensile strength) is small enough compared to the characteristic dimension of the structure (L mat < 0.05h for the studied case of a bimaterial specimen submitted to a four point bending test, where h is the substrate thickness).

The theory of the CC can be found in several papers in the case of MA (e.g., [START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavities size effects[END_REF][START_REF] Yosibash | A failure criterion for brittle elastic materials under mixed-mode loading[END_REF][START_REF] Priel | Failure initiation at a blunt V-notch tip under mixed mode loading[END_REF][START_REF] García | Mixed-mode crack initiation at a v-notch in presence of an adhesive joint[END_REF][START_REF] Leguillon | An attempt to extend the 2D coupled criterion for crack nucleation in brittle materials to the 3D case[END_REF]) and FFE approaches (e.g., [START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF][START_REF] Doitrand | Mesoscale analysis of damage growth in woven composites[END_REF][START_REF] Doitrand | Damage onset modeling in woven composites based on a coupled stress and energy criterion[END_REF][START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Martin | A coupled strength and toughness criterion for the prediction of the open hole tensile strength of a composite plate[END_REF][START_REF] Martin | Understanding the tensile strength of ceramics in the presence of small critical flaws[END_REF]). However, contrary to the other approaches for brittle fracture modeling cited previously (such as cohesive zone or phase field models), dedicated tools are not available for the CC application in commercial FE codes. The objective of this paper is to fill this gap and to provide a FE procedure for the implementation in Abaqus/Standard of both MA and FFE applications of the CC, and to offer practical recommendations for their use. In Section 2, the CC theory is recalled. Details on FE implementation of both approaches are described in Section 3. Numerical results are given in section 4.

The coupled criterion

In this section, the CC theory is recalled in the case of MA (cf. Section 2.1) and FFE (cf. Section 2.2) approaches. For the sake of pedagogy, it will only be presented in the bi-dimensional case of an isotropic homogeneous linear elastic material loaded in pure mode I. The MA approach applied to a V-notch in the general mixed mode case is described in the Appendix. As stated in introduction, the CC requires the simultaneous fulfillment of two separate conditions to ensure crack initiation.

The first condition requires that the stress σ normal to the crack plane overcomes the material tensile strength σ c on the whole presupposed crack path of length l and described by its arc length parameter s:

σ(s) ≥ σ c , 0 ≤ s ≤ l. (1)
This stress condition is only a necessary condition for crack initiation. In the case of a monotonically decreasing stress, which is classicaly encountered at singular or stress raisers, the distance l s on which stress condition is fulfilled increases with increasing imposed loading. Moreover, if the stress condition is satisfied over a distance l s for a certain imposed loading, it is also verified for all distances l such that l ≤ l s (Fig. 1a). Therefore, the stress condition allows determining an upper bound for all the admissible crack lengths at initiation.

The second condition originates from a balance of the kinetic (δW k ), potential (δW p ) and crack surface (G c l, where G c is the material fracture toughness) energies between the state prior to and after the nucleation of a crack of length l:

δW k + δW p + G c l = 0. ( 2
)
A quasi-static initial state implies that δW k ≥ 0, which allows deriving the energy condition:

G inc (l) = -δW p l ≥ G c . (3)
G inc is called the incremental energy release rate and reverts to the Griffith definition of the energy release

Stress criterion fulfilled

Energy criterion fulfilled Initiation configuration

Crack length lower bound Crack length upper bound rate G when l → 0. Note that G can be expressed as a function of G inc [START_REF] Weißgraeber | A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers[END_REF][START_REF] Doitrand | Mesoscale analysis of damage growth in woven composites[END_REF][START_REF] Doitrand | Comparison between cohesive zone and coupled criterion modeling of crack initiation in rhombus hole specimens under quasi-static compression[END_REF][START_REF] Martin | Energetic conditions for interfacial failure in the vicinity of a matrix crack in brittle matrix composites[END_REF], which can be useful to study crack propagation after initiation. However, this point will not be treated in this work which is exclusively dedicated to the prediction of crack initiation. We consider here a positive geometry which implies that the crack propagates in an unstable manner after initiation. This corresponds to a monotonically increasing variation of G inc . Different situations (corresponding to negative geometries for which the crack propagation is stable after initiation) can also be encountered (cf. e.g., [START_REF] Weißgraeber | A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers[END_REF][START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF][START_REF] Martin | Energetic conditions for interfacial failure in the vicinity of a matrix crack in brittle matrix composites[END_REF][START_REF] Doitrand | Damage onset modeling in woven composites based on a coupled stress and energy criterion[END_REF][START_REF] García | Transverse crack onset and growth in cross-ply laminates under tension: application of a coupled stress and energy criterion[END_REF][START_REF] Martin | A twofold strength and toughness criterion for the onset of freeedge shear delamination in angle-ply laminates[END_REF]) . Therefore, for a given imposed loading, if the energy criterion is attained for a crack length l e , it is also attained for all crack lengths l such that l ≥ l e (Fig. 1a). The energy condition thus allows determining a lower bound for all the admissible initiation crack lengths. For a too small imposed loading (lower than the initiation loading), the lower bound l e for admissible initiation crack length provided by the energy criterion is higher than the upper bound l s provided by the stress condition hence crack initiation is not possible (cf. Fig. 1a). With increasing loading, l s increases whereas l e decreases and the initiation loading and crack length l c can be determined as the loading for which l c = l e = l s which ensures that both the stress and the energy conditions are simultaneously fulfilled (cf. Fig. 1b). In the following, as a matter of example, the case of a V-notch specimen (whose dimensions are depicted in Fig. 2) under three-point bending is studied. For the sake of simplicity, a sharp V-notch will be considered. The scripts for the MA and FFE approaches provided as supplementary materials (section 6.4 in the Appendix) also allow to analyze the case of a blunted V-notch [START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavities size effects[END_REF][START_REF] Priel | Failure initiation at a blunt V-notch tip under mixed mode loading[END_REF][START_REF] Leguillon | Crack onset at a v-notch. Influence of the notch tip radius[END_REF].

h t 2L x y β
Figure 2: Geometry of the V-notch specimen submitted to three-point bending.

MA approach

In the example depicted in Fig. 2, a crack initiates at the V-notch in mode I in y-axis direction. In the vicinity of the V-notch tip, under plane elasticity assumption, the outer asymptotic expansion of the far displacement field is given by:

U (x, y) = U (0, 0) + kr λ u(θ) + ... (4)
(x, y) and (r, θ) respectively hold for the cartesian and polar coordinates, U (0, 0) denotes the rigid translation of the origin located at the tip of the V-notch, λ is the singularity exponent and u(θ) is an angular shape function. They are solutions to an eigenvalue problem [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF] and depend on the V-notch opening angle β. For V-notch opening angles comprised between 0 (crack) and π (straight edge), the exponent λ varies between 0.5 and 1. k is the critical stress intensity factor (GSIF) and is proportional to the applied displacement or load.

The computation of the proportionality factor will be explained at the end of the section and illustrated in Section 3. The application of the CC through the MA approach consists in performing a zoom on a domain surrounding the crack initiation location. Dilating the space variables provides a (theoretically) unbounded domain in which the length of the crack is now equal to 1. This domain is classically denoted as the inner one and is used to define the inner asymptotic expansion of the near field. As shown by previous authors [START_REF] Leguillon | Fracture in heterogeneous materials -Weak and strong singularities[END_REF],

matching the inner expansion (near field) and the outer expansion (far field), the change δW p in potential energy due to crack initiation over a length l is:

-δW p (l) = k 2 l 2λ B + ..., (5)
where B is a geometrical coefficient depending only on the V-notch opening angle β. The energy condition (cf. Eq. ( 3)) then rewrites:

G inc = - δW p l = k 2 l 2λ-1 B ≥ G c . (6)
From Eq. ( 4) and using an appropriate normalization of the eigenvector u(θ) [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF][START_REF] Leguillon | Crack onset at a v-notch. Influence of the notch tip radius[END_REF], Eq. ( 1) rewrites as it is indicated in the Appendix (section 6.1):

σ(l) = kl λ-1 ≥ σ c . (7)
By combining Eqs. ( 6) and ( 7), the initiation crack length l c and the critical stress intensity factor k c can be derived:

l c = G c Bσ 2 c = 1 EB L mat , with L mat = EG c σ 2 c , (8)
k c = ( G c B ) 1-λ σ 2λ-1 c . (9)
It is worth noting that tabulated values of the coefficient B can be found in a previous paper [START_REF] Leguillon | Crack onset at a v-notch. Influence of the notch tip radius[END_REF]. The MA approach allows computing the generalized stress intensity factor k c and crack length l c at initiation without considering the whole structure (as depicted, e.g., in Fig 2). Eq. ( 8) indicates that l c is proportional to the material length L mat . The smallest value is obtained for small opening angles [START_REF] Leguillon | Crack onset at a v-notch. Influence of the notch tip radius[END_REF] with:

l c ≥ l min c = L mat 2π . (10)
Eq. ( 9) includes the usual Griffith criterion for a crack (β = 0, λ = 0.5) and the strength condition for a straight edge (β = π, λ = 1).

The link between the load applied on the structure and the generalized stress intensity factor k EF can be established by computing the following contour integral from a FE analysis on the uncracked structure:

k EF = Ψ(U EF , r -λ u -(θ)) Ψ(r λ u ( θ), r -λ u -(θ)) , (11)
with

Ψ(U a , U b ) = Γ σ(U a ) • N • U b -σ(U b ) • N • U a dl,
where Γ is any closed contour surrounding the V-notch and the crack extension and ending on the two faces of the V-notch and N the unit normal to Γ pointing toward the V-notch tip. The term r -λ u -(θ) is the so-called dual function to the singular mode r λ u(θ) [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF]. An example of calculation of the stress intensity factor using Eq. ( 11) is provided as a Supplementary material.

In the case of a blunted V-notch, for which the CC solution is given in Appendix (cf. Section 6.3), the inner domain is obtained by dilatation of the coordinates with respect to the blunted V-notch radius d.

Hence, contrary to the sharp V-notch for which the crack length in the inner domain is 1 (all quantities being dilated with respect to it), a non-dimensional crack length η = l/d is defined.

FFE approach

In the FFE approach, the whole structure under investigation is considered rather than only a domain surrounding the crack. The stress criterion can thus be computed using only one calculation on the undamaged configuration. Under linear elasticity and the assumption of small deformations, the local stress is proportional to the applied displacement U 0 so that Eq. ( 1) rewrites:

σ(s) = ξ(s)E U 0 w ≥ σ c , 0 ≤ s ≤ l. (12)
ξ(s) is a dimensionless function depending only on the structure geometrical features and w a characteristic length of the studied structure (selected here as the sample height as schematised in Fig. 2). Eq. ( 12) can be rewrited so as to determine the minimum imposed displacement U stress 0 necessary to fulfill the stress criterion on a distance l:

U stress 0 (l) = w σ c E 1 ξ(l) . (13)
The energy criterion requires several calculations with varying crack length in order to compute the potential energy variation, and hence the incremental energy release rate. Under linear elasticity and the assumption of small deformations, the potential energy is proportional to the square applied displacement so that Eq.

(3) provides:

G inc (l) = A(l)Ew( U 0 w ) 2 ≥ G c . (14)
A(l) is a dimensionless function depending only on the structural geometry. Eq. ( 14) supplies the minimum imposed displacement U energy 0 necessary to fulfill the energy criterion for a crack length l:

U energy 0 (l) = w G c A(l)Ew = w σ c E L mat w 1 A(l) , (15)
where Lmat w is similar to the brittleness number introduced by Mantic [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF].

The initiation displacement U c and crack length l c can thus be determined as the minimum applied displacement for which both criteria are fulfilled, i.e.:

l c = argmin l (max(U energy 0 (l), U stress 0 (l))), (16)
U c = max(U energy 0 (l c ), U stress 0 (l c )). ( 17
)
In the case of monotonically decreasing stress and monotonically increasing energy release rate, the initiation loading predicted by the stress and the energy conditions are equal so that the initiation length l c can be estimated by combining and solving Eqs. ( 13) and ( 15):

ξ(l c ) 2 A(l c ) = σ 2 c w G c E = w L mat . (18)
The applied displacement at initiation is then obtained with:

U c w σc E = L mat w 1 A(l c ) = 1 ξ(l c ) . (19)
Eq. ( 18) reveals that the initiation length l c depends on the material length L mat and on the structural geometry. In the case of a small notch defect submitted to traction, it was shown [START_REF] Martin | Understanding the tensile strength of ceramics in the presence of small critical flaws[END_REF] that the upper bound of the initiation length is:

l c ≤ l max c = L mat 2 . ( 20
)
It must be noted that the value of this upper bound was obtained while considering a monotonically increasing behaviour of the incremental energy release rate versus the crack length. Configurations for which B(l) exhibits a maximum are qualified as negative geometries [START_REF] Weißgraeber | Crack nucleation in negative geometries[END_REF]. For these configurations, such as for instance debonding initiation along an interface, l max c is a structural value which only depends on the geometry [START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF].

Finite element procedure -The case of a V-notch specimen

The numerical implementations of the MA and FFE approaches are described in this section, based on the three-point bending geometry depicted in Fig. 2.

Implementation of the MA approach

The asymptotic approach of the CC requires the coefficient B to be estimated so as to determine the initiation length and stress intensity factor at initiation with the help of Eqs. ( 8) and ( 9). The computation of B is based on the potential energy change due to the crack initiation (cf. Eq. ( 6)). The potential energies in presence of and without a crack of length unity are therefore assessed in the inner domain. An example of mesh of the inner domain is presented in Fig. 3. The displacement field related to the V-notch (cf. Eq. ( 4)

and Eq. ( 24) in the Appendix) is imposed as a Dirichlet condition on the domain boundary and the V-notch faces are stress-free. It is recommended to use the exact same mesh topology to compute the potential energy change and to double the node lying on the crack path in the cracked case since it reduces the numerical errors made on the potential energy change. Note that another possibility to calculate the potential energy change consists in using the contour integral Ψ defined by Eq. ( 11) [START_REF] Yosibash | A failure criterion for brittle elastic materials under mixed-mode loading[END_REF].

The inner domain is artificially bounded at a large distance and boundary conditions are necessary to prescribe a behavior at infinity (matching conditions), thus it is recommended to ensure that the inner domain is large enough. Following the convergence analysis presented in next section, we recommend to use an inner domain radius at least 200 times larger than the largest characteristic length (e.g., crack length, blunt notch radius) in the inner domain. Furthermore, the minimum mesh size m MA at the crack tip must be sufficiently small with m MA ≤ l 100 . The MA procedure also requires that the initiation length remains small compared to the defect length. Using Eq. ( 20), which gives an estimate of the upper bound of the initiation length, leads to:

L mat ≤ h 5 . (21)
As will be shown in next section, Eq. ( 21) must be satisfied to avoid any artefact introduced by the MA approach. In the present example, the stress variation can be derived analytically (Eq. ( 26) in the Appendix).

In other cases for which no analytical solutions exist, it is also possible to compute the stress condition using a FE calculation in the inner domain without crack. It does not require any supplementary calculations and stress intensity factor k c using Eqs. ( 8) and ( 9).

It can be noted that the FE calculations in the inner domain can be done once and for all with given values of (E,ν)=(E 1 ,ν 1 ) for which a value of the coefficient B = B 1 is computed. Indeed, it can be shown that the coefficient B 2 corresponding to another couple (E 2 ,ν 2 ) can be expressed as a function of B 1 and

(E 1 , ν 1 , E 2 , ν 2 ) with B 2 = B 1 1-ν 2 2 E2 E1 1-ν 2 1 .
We provide as supplementary data (section 6. are used. Then, a second script solves the CC (step 4) (CC Solution from FE calculations.py). Those two scripts can also be used to consider a mixed-mode loading. In this case, the crack may initiate with an orientation γ as schematized in Fig. 10. Minimizing k c I (γ) offers an estimate of the initiation angle γ c . The relevant formulation of the CC is given by Eqs. ( 31) and (32) in section 6.2 of the Appendix. The calculations for several V-notch openings β and crack angles γ have been performed (with E= 1 MPa and ν=0.36) and the results were stored in a database file (Data Vnotch r0.txt) provided as a Supplementary material. Similar results (only with γ = 0) were already given in a previous paper [START_REF] Leguillon | Crack onset at a v-notch. Influence of the notch tip radius[END_REF]. A python script (Solution Vnotch r0 modeIandII.py) allows solving the CC (step 4) in the case of a sharp V-notch (under mode I and mixed mode loading) using the MA without any FE calculations. It only requires these precalculated data and is therefore computationally very efficient.

Eqs. ( 43) and (44) in Section 6.3 of the Appendix give the formulation of the CC for a blunt V-notch under mixed-mode loading. The same scripts (ScriptV notch.py) and (CC Solution from FE calculations.py) can be used to solve the CC with the MA approach. As already mentioned, the inner domain is obtained by dilatation of the spatial coordinates with respect to the blunted V-notch radius d and a non-dimensional crack length η = l/d is introduced. Therefore, an estimate of an upper bound to the initiation crack length must be obtained before generating the FE mesh of the inner domain including a blunted V-notch radius of dimension 1 and a crack of non dimensional length η. In practice, the following procedure must be followed so as to compute the maximum non-dimensionnal crack length in the inner domain:

1) Choose the minimum blunted V-notch radius d min for which the CC solution has to be calculated, 2) Estimate the maximum initiation length with the condition l c ≤ Lmat 2 ,

3) Compute the non-dimensional upper bound for the initiation length Lmat 2dmin that can be employed to generate the inner domain mesh.

Then the MA approach can be used as already described in Section 2.1.

Implementation of the FFE approach

Crack initiation length l c and loading level U c are determined using Eqs. ( 18) and [START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF], which involves the dimensionless functions A(l) and ξ(s). The computation of ξ only requires one calculation on the structure without crack, whereas several calculations with increasing crack lengths are requested so as to obtain A.

Similarly to the MA approach, it is recommended to use the exact same mesh topology to compute the potential energy change by successively releasing the nodes lying on the crack path. A refined mesh of length M FFE with the mesh size m FFE (Fig. 4) must be introduced in the vicinity of the notch in order to capture accurately the initiation length l c . Using the lower bound l min c defined by Eq. ( 10) and the upper bound l max c defined by Eq. ( 20), a convergence analysis presented in next section leads to select: If the node distribution is regular, the number of calculations necessary to evaluate the incremental energy release rate A is N ≥ MFFE mFFE = 6π ∼ = 20. It will be shown in next section that N can be reduced by introducing a controlled variation of the mesh size within the refined zone. The steps to follow for a FE implemention of the FFE approach are the following:

m FFE ≤ l min c 6 = L mat 12π , (22)
M FFE ≥ l max c = L mat 2 . ( 23
)
1) Generation of the geometry and mesh of the studied structure, 2) Generation of several meshes with varying crack length, 3) Introduction of the material properties (E, ν) and computation of the scale factors ξ and A, 4) Introduction of the material properties (σ c , G c ) and calculation of the initiation crack length and loading level using Eqs. ( 18) and [START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF].

In practice, the dimensionless functions ξ and A do not depend on Young's modulus E and reveal weakly dependant on Poisson's ratio ν so that the material properties are only needed for the last step. We provide as supplementary data (section 6.4 in the Appendix) a python source code (F3P Vnotch.py) for steps 1 to 3 (using Abaqus/Standard) and a complementary one (post F3P.py) for step 4. Standard linear quadratic elements (plane strain conditions) are still used.

Numerical results

This section includes several convergence analysis and a comparison between MA and FFE approaches.

Convergence analysis to estimate the minimum mesh size and the radius of the inner domain (MA approach)

Convergence analyses are performed in order to determine the best meshing conditions of the inner domain. Fig. 5 shows the influence of the minimum mesh m MA size at the crack tip to crack length l ratio on the coefficient B. The condition mMA l ≤ 0.01 must be satisfied to reach convergence. The effect of the radius R of the inner domain on the initiation stress intensity factor and crack length is depicted on Fig. 6.

Both quantities are overestimated using a too small inner domain radius and converge to a constant value for larger domain radii with R ≥ 200 l. It is thus recommanded to keep the radius of the inner domain 200 times larger than the largest characteristic length. results were obtained for another angles β and V-notch depth h, which leads to recommand a mesh size

m FFE ≤ l min c 6 = Lmat 12π .
Employing a non-uniform mesh is useful to reduce the number of nodes within the refined zone. This is illustrated in Fig. 8b: starting from a mesh size m FFE =1.5 µm, a bias is introduced in order to verify the condition [START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF] close to the notch tip while allowing progressively a larger mesh size away from the notch tip. This procedure does not degrade or improve the accuracy but allows to reduce the number of calculations required to evaluate the scale factor A.

U c (m FFE )-U c (m REF ) U c (m REF )
in terms of applied displacement at initiation versus the characteristic length (FFE approach) m REF =0.5 µm : a) The mesh is uniform, b) The mesh is non-uniform and obtained after introducing a bias.

The three-point bending geometry with a sharp V-notch is analysed with 2L = 40 mm, w = 2 mm, h = 0.3 mm, β = 40deg. and E=400 GPa, ν = 0.2 and M FFE =200µm.

Comparison between MA and FFE approaches

MA and FFE approaches are now used to estimate the influence of the V-notch angle on crack initiation for the three-point bending geometry (Fig. 2) with 2L = 40 mm, w = 2 mm and h = 0.3 mm. The minimum mesh size and the radius of the inner domain employed with the MA method are respectively m MA =10µm

and R=200mm with l=1mm. The minimum mesh size and the length of the refined zone employed with the FFE method are respectively m FFE = m REF =0.5 µm and M FFE =200 µm. The MA approach requires the calculation of the GSIF for a given applied displacement with the contour integral defined by Eq. ( 11) (which can be obtained using the script Compute K2D.py provided in supplementary materials). Results displayed in Fig. 9a show that the applied critical displacement increases with the notch angle. Both procedures give similar results for lower values of the characteristic length but a significant difference appears with increasing L mat . The explanation is provided by Fig. 9b which plots the initiation length at crack initiation versus the notch angle: l c increases with L mat and β. As the scaling factor linking l c and L mat reveals smaller for the FF approach for higher values of L mat , smaller initiation lengths are estimated with the FF method for larger values of L mat . To remain within the framework of the asymptotic analysis, the initiation length must be at least one order of magnitude smaller than the notch length with l c ≤ h 10 =30 µm. This condition is no more satisfied for every notch angle for increasing values of L mat for which the MA results exhibits (in Fig. 9a) a non-physical minimum corresponding to a "critical angle" as already shown by previous authors [START_REF] Carpinteri | The problem of the critical angle for edge and center V-notched structures[END_REF].

Conclusion

To predict crack initiation in the vicinity of a stress concentration, the CC can be implemented with the help of MA or FFE methods. Both procedures are detailed in this paper in order to facilitate their use with a FE code. In addition, python source codes are provided to illustrate their implementation in Abaqus/Standard in the case of a sharp or a blunted V-notch. The supplied tools apply the MA approach for mixed mode loading and the FFE approach for a specimen submitted to three-point bending. Results of convergence analysis are analysed to link the characteristic length of the studied material with the needed mesh refinement at the notch tip for the FFE method. The MA approach allows the initiation loading and crack length to be determined without considering the whole structure under investigation by performing a zoom around the crack initiation location. It is computationally more efficient than the FFE approach that requires several calculations with different crack lengths to be computed on the whole model. The domain in which the calculations are performed must be large enough to prescribe a behavior "at infinity" (see Section 3.1). The asymptotic approach is restricted by the assumption that the initiation length is small compared to the structure characteristic dimensions, and that the crack location is far enough from any boundary or features that may affect the displacement/stress field around the crack. If these assumptions are not fulfilled, some differences on initiation crack length and loading level are obtained between both approaches. This is directly linked to the material characteristic length L mat which has an influence on the predicted initiation length (that also depends on the geometrical features of the studied structure). An additional study is needed to quantify more precisely these differences, nevertheless condition (21) provides a validity range estimate for the MA approach.

Appendix

The CC solution in mixed mode using the asymptotic approach is explained in this appendix for either sharp or blunted V-notch, which of course includes the particular case of a sharp V-notch specimen loaded under pure mode I loading presented in the paper for the sake of pedagoy. The scripts provided as supplementary data allow studying general cases of blunted V-notch under mixed mode loading.

Singularity exponents and asymptotic displacement field

The displacements in the vicinity of a traction free V-notch can be expressed by an asymptotic series, for which the first three terms expressed in polar coordinates are:

u(r, θ) =    u r u θ    = u(0, 0) + k I r λ1    u I r (θ) u I θ (θ)    + k II r λ2    u II r (θ) u II θ (θ)    (24)
with

u I r (θ) = [cos((1 + λ 1 )θ * ) + λL+3µL-λ1(λL+µL) (λL+µL)(1-λ1) sin(ω(1+λ1)/2) sin(ω)(1-λ1)/2) cos((1 -λ 1 )θ * )]/(2µ L λ 1 σ I-0 θθ ) u I θ (θ) = [-sin((1 + λ 1 )θ * ) -λL+3µL+λ1(λL+µL) (λL+µL)(1-λ1) sin(ω(1+λ1)/2) sin(ω)(1-λ1)/2) sin((1 -λ 1 )θ * )]/(2µ L λ 1 σ I-0 θθ ) u II r (θ) = [sin((1 + λ 2 )θ * ) + λL+3µL-λ2(λL+µL) (λL+µL)(1+λ2) sin(ω(1+λ2)/2) sin(ω)(1-λ2)/2) sin((1 -λ 2 )θ * )]/(2µ L λ 2 σ II-0 rθ ) u II θ (θ) = [cos((1 + λ 2 )θ * ) + λL+3µL+λ2(λL+µL) (λL+µL)(1+λ2) sin(ω(1+λ2)/2) sin(ω)(1-λ2)/2) cos((1 -λ 2 )θ * )]/(2µ L λ 2 σ II-0 rθ ) a
where λ L , µ L are the Lamé coefficient, θ * = θ -π/2, ω = 2π -β and

σ I-0 θθ = (1+λ1)sin(ω(1+λ1)/2) (1-λ1)sin(ω(1-λ1)/2) -1 σ II-0 rθ = 1 -(1-λ2)sin(ω(1+λ2)/2) (1+λ2)sin(ω(1-λ2)/2)
The singularity exponents λ 1 and λ 2 associated with a V-notch with opening angle β can be computed as the smallest roots of the following characteristic equations:

   sin(λ 1 (2π -β)) + λ 1 sin(2π -β) = 0 sin(λ 2 (2π -β)) -λ 2 sin(2π -β) = 0 (25)
The stress field in the vicinity of the V-notch can be expressed in polar coordinates as:

σ(r, θ) =          σ rr σ θθ σ rθ          = k I r λ1-1          σ I rr (θ) σ I θθ (θ) σ I rθ (θ)          + k II r λ2-1          σ II rr (θ) σ II θθ (θ) σ II rθ (θ)          (26)
where γ defines the crack angle with respect to the V-notch bisector (Fig. 10). Therefore the energy condition of the CC writes:

G inc = - δW p l = k 2 I l 2λ1-1 (B 11 (γ) + m(l)B 12 (γ) + m(l) 2 B 22 (γ)) G c (28)
The stress condition of the CC writes:

σ θθ = k I l λ1-1 (σ I θθ (γ) + m(l)σ II θθ (γ)) σ c (29)
Therefore, the generalized stress intensity factor k I can be computed by combining Eqs. ( 28) and ( 29

For a given V-notch opening angle β, the critical GSIF is estimated for several crack angles γ, the initiation GSIF and crack angle being given by the configuration minimizing k c I . In practice, the solution of the CC requires the knowledge of the functions B ij (i,j=1 or 2) and σ k θθ (k=I or II). The latter can be computed analytically in the case of a sharp V-notch (cf. Eq. ( 26)) while the former have to be computed by FE computations. First, FE calculations are performed for (k I =1, k II =0) in order to estimate the potential energies with (δW I p (l)) and without (δW I p (0)) a crack in pure mode I (and hence their difference -δW I p =δW I p (0) -δW I p (l)) from the stress σ I and strain I fields.

Figure 1 :

 1 Figure 1: Applied stress to strength ratio σ σc and incremental energy release rate to fracture toughness ratio G inc Gc as a function of the crack length l for a loading (a) lower than the initiation loading (the lower bound of admissible initiation length le is higher than the upper bound ls) and (b) equal to the initiation loading (the lower bound of admissible initiation length le is equal to the upper bound ls and to the initiation length lc).

1 )

 1 than those already performed to compute the potential energy change, but only a post-processing of the FE calculation in the inner domain with no crack. The steps to follow in a FE implementation of the MA approach are the following: Generation of the inner domain geometry and mesh, 2) Computation of the displacement boundary conditions imposed on the inner domain boundary, 3) Computation of the potential energy difference and thus B coefficient, B = WP (0) -WP (1) where WP denotes the potential energy computed in the inner domain and where 0 and 1 denote respectively the uncracked and cracked cases, 4) Introduction of the material properties (E, ν, σ c , G c ) and calculation of the initiation crack length l c

Figure 3 :

 3 Figure 3: Inner domain of radius R in the case of a V-notch with angle β containing a unity length crack and an example of mesh employed for the computation of the energy difference (MA approach).

Figure 4 :

 4 Figure 4: Refined mesh at the V-notch tip (FFE approach): the minimum mesh size and the length of the refined zone are respectively denoted m FFE and M FFE .

Figure 5 :

 5 Figure 5: Coefficient B (MA approach) versus the minimum mesh size to crack length ratio m MA l for a sharp V-Notch with

4. 2 .

 2 Fig. 7 depicts the crack length and applied displacement at initiation versus the ratio Lmat w for given values of h and β. For large values of the minimum mesh size, l c and U c are overestimated but decreasing m FFE towards the reference value m REF =0.5 µm clearly demonstrates a convergence. Fig. 8a shows that a value m FFE ≤ 1 µm is necessary to maintain the relative error with the reference solution below 1%. Similar

Figure 6 :Figure 7 :

 67 Figure 6: Initiation (a) critical stress intensity factor and (b) crack length determined with the MA approach as a function of the inner domain radius to crack length ratio for a sharp V-Notch with β = 90deg., E = 3300MPa , ν =0.36 , σc =72MPa and Gc =0.250MPa.mm. The minimum mesh size to crack length ratio is m MA l =0.01.

m=1. 5 (Figure 8 :

 58 Figure 8: Relative error

Figure 9 :

 9 Figure 9: Comparison between MA (square symbols and full lines) and FFE (disk symbols and dotted lines) approaches: a) Applied displacement and b) Initiation length at crack onset versus the notch angle for various values of the characteristic length Lmat. A three-point bending specimen with a sharp V-notch is considered with 2L = 40 mm, w = 2 mm, h = 0.3 mm and E=400 GPa, ν = 0.2, σ c = 100 MPa and Gc=1, 3 or 10 J/m 2 corresponding to Lmat=40, 120, or 400 µm.

Figure 10 :

 10 Figure 10: Inner domain of radius R containing a V-notch with opening angle β and tip radius d and a crack oriented with an angle γ with respect to Oy axis.

rθ 6 . 2 .-δW p = k 2 I

 622 with σ I rr (θ) = [cos((1 + λ 1 )θ * ) + 3-λ1 1-λ1 sin(ω(1+λ1)/2) sin(ω)(1-λ1)/2) cos((1 -λ 1 )θ * )]/σ I-0 θθ σ I θθ (θ) = [-cos((1 + λ 1 )θ * ) + 1+λ1 1-λ1 sin(ω(1+λ1)/2) sin(ω(1-λ1)/2) cos((1 -λ 1 )θ * )]/σ I-0 θθσ I rθ (θ) = [-sin((1 + λ 1 )θ * ) + sin(ω(1+λ1)/2) sin(ω(1-λ1)/2) sin((1 -λ 1 )θ * )]/σ I-0 θθ σ II rr (θ) = [sin((1 + λ 2 )θ * ) + 3-λ2 1+λ2 sin(ω(1+λ2)/2) sin(ω(1-λ2)/2) sin((1 -λ 2 )θ * )]/σ II-0 θθ σ II θθ (θ) = [-sin((1 + λ 2 )θ * ) + sin(ω(1+λ2)/2) sin(ω(1-λ2)/2) sin((1 -λ 1 )θ * )]/σ II-0 θθ σ II rθ (θ) = [cos((1 + λ 2 )θ * ) -1-λ2 1+λ2 sin(ω(1+λ2)/2) sin(ω(1-λ2)/2) cos((1 -λ 2 )θ * )]/σ II-0The sharp V-notch under mixed mode loadingIn the context of matched asymptotic expansions, the inner domain is dilated with respect to the crack length l. In this case the potential energy difference for a V-notch with opening angle β is expressed as: l 2λ1 (B 11 (γ) + m(l)B 12 (γ) + m(l) 2 B 22 (γ)) m(l) = kII kI l λ2-λ1

  ):k I = ( G c B 11 (γ) + m(l)B 12 (γ) + m(l) 2 B 22 (γ) ) 1-λ1 ( σ I θθ (γ) + m(l)σ II θθ (γ) σ c ) 1-2λ1(30)Using Eq. (28), the initiation crack length l c can be obtained by solving the following equation which is implicit because m depends on l c :l c = G c B 11 (γ) + m(l c )B 12 (γ) + m(l c ) 2 B 22 (γ) ( σ I θθ (γ) + m(l c )σ II θθ (γ) 11 (γ) + m(l c )B 12 (γ) + m(l c ) 2 B 22 (γ) ) 1-λ1 ( σ c σ I θθ (γ) + m(l c )σ II θθ (γ) ) 2λ1-1

  It provides B 11 for a given crack angle γ (the procedure being repeted for crack angles γ varying between -90 and 0 deg.), from Eq. (27):B 11 (γ) = -δW I p k 2 I l 2λ1 = δW I p (0) -δW I p (l) since k I = 1 and l = 1(34)Then, FE calculations are performed for (k I =0, k II =1) in order to evaluate the potential energies with (δW II p (l)) and without (δW II p (0)) a crack in pure mode II (and hence their difference -δW II p =δW II p (0) -δW II p (l)) from the stress σ II and strain II fields. A 22 for a given crack angle γ (the procedure being repeted for crack angles γ varying between -90 and 0 deg.), from Eq. (27):B 22 (γ) = -δW II p k 2 II l 2λ2 = δW II p (0) -δW II p (l) since k II = 1 and l = 1(36)Finally, the knowledge of B 11 and B 22 gives access to B 12 by computing the potential energy difference -δW I+II p for any GSIFs so that (k I ,k II ) =(0,0). For the sake of simplicity, we choose (k I =1, k II =1). B 12 (γ) = 1 m(l) ( -δW I+II p k 2 I l 2λ1 -A 11 (γ) -m(l) 2 A 22 (γ))

Actually, no extra FE calculations than those in pure mode I and pure mode II are required to compute B 12 since the potential energy in mode I+II can be computed by combining the stress and strain fields obtained in pure mode I and II:

2 Ω (σ I + σ II ) : ( I + II )dV [START_REF] Carpinteri | The problem of the critical angle for edge and center V-notched structures[END_REF] hence simplifying Eq. ( 37) into:

6.3. The blunted V-notch under mixed mode loading

In the case of a V-notch with a blunted tip of radius d (cf. Fig. 10), the inner domain is dilated with respect to either the crack length l or the notch radius d (leading to the same results). For practical reasons (ease of FE mesh generation), we choose a dilatation with respect to d hence defining the normalized crack length η = l/d. The potential energy difference reads:

Therefore the energy condition of the CC writes:

The stress condition of the CC writes:

where σ kb θθ differs from σ k θθ (k = I or II) defined in Section 6.2 because of the blunted notch tip. These functions can be computed by FE calculations on uncracked configuration. By combining Eqs. (41) and (42), we obtain the equation that must be solved so as to determine the initiation crack length l c = dη c for a given value of γ:

It can be noted that the solution of Eq. ( 43) depends on the actual value of the the blunted notch radius d (expressed in length unit) and not on a ratio between d and another characteristic length. From (41) and

(42) we deduce the critical GSIF:

The crack initiation angle γ c is solved as the one minimizing the critical GSIF k c I . In practice, the solution of the CC for the blunted V-notch requires the knowledge of the functions B b ij (i,j=1 or 2) and σ K θθ b (k=I or II). Contrary to the sharp V-notch case, the latter are not computed analytically but by means of FE calculations, as well as the former. For a given V-notch opening, B b ij must be computed for several crack 17 angles γ and normalized lengths η c in order to solve Eq. ( 43) for all these angles and find the one minimizing k c I . In fact, it was shown in [START_REF] Priel | Failure initiation at a blunt V-notch tip under mixed mode loading[END_REF] that the same crack initiation angle as for a sharp V-notch was found for a blunted V-notch. Therefore, the computation of B b ij can directly be performed for the crack angle γ c that minimizes the GSIF in the case of a sharp V-notch, which reduces the number of calculations required to solve the CC. The same approach as described in Section 6.2 is adopted to compute the B b ij functions.


Hafiz Ahmed
email: hafiz.h.ahmed@ieee.org

Rosane Ushirobira

Denis Efimov

Leonid Fridman

Yongqiang Wang

Output Global Oscillatory Synchronization of Heterogeneous Systems

Keywords: Output synchronization, Sliding mode observation, Van der Pol oscillator, Real-time control

The global output synchronization problem for heterogeneous nonlinear systems having relative degree 2 or higher is studied. The proposed approach consists in two steps. First, a partial projection of individual subsystems into the Brockett oscillators is performed using a sliding-mode control. Second, the network of these oscillators is synchronized using the global synchronization results of a particular second order nonlinear oscillator model from Ahmed et al. (2019). Our approach is based on output feedback and uses a higher order sliding mode observer to estimate the states and perturbations of the synchronized nonlinear systems. Along with numerical simulations, the performance of the proposed synchronization scheme is experimentally verified on a network of Van der Pol oscillators.

Introduction

Over the last decade, the synchronization of complex dynamical systems and/or network of systems has attracted a great deal of attention from multidisciplinary research communities due to their pervasive presence in nature, technology and human society [START_REF] Blekhman | Synchronization in science and technology[END_REF]; [START_REF] Pikovsky | Synchronization: a universal concept in nonlinear sciences[END_REF]; [START_REF] Strogatz | Sync: How order emerges from chaos in the universe, nature, and daily life[END_REF]; [START_REF] Osipov | Synchronization in oscillatory networks[END_REF]]. Among potential application domains of synchronization, it is worth to mention the smooth operations of microgrid [START_REF] Efimov | Robustness of delayed multistable systems with application to droop-controlled inverter-based microgrids[END_REF]; [START_REF] Schiffer | A consensus-based distributed voltage control for reactive power sharing in microgrids[END_REF]], secure communication [START_REF] Martínez-Guerra | Secure communications via synchronization of liouvillian chaotic systems[END_REF]; Fradkov & Markov (1997)], deployment of mobile sensor networks [START_REF] Wang | Energy-efficient pulse-coupled synchronization strategy design for wireless sensor networks through reduced idle listening[END_REF]], formation control [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF]], chaos synchronization [START_REF] Rodriguez | Synchronization in reduced-order of chaotic systems via control approaches based on high-order sliding-mode observer[END_REF]], genetic oscillators [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF]], etc.

Significant progress have been made in the past decade in the area of control design for synchronization, consensus or motion coordination, the existing literature is huge and covers a wide area of topics [START_REF] Gazi | Swarm stability and optimization[END_REF]; [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF]; [START_REF] Panteley | Synchronisation and emergent behaviour in networks of heterogeneous systems: A control theory perspective[END_REF]]. Until now, a large number of works are available on the problem of synchronization of networks with identical nodes, particularly when the nodes are linear time-invariant systems [START_REF] Scardovi | Synchronization in networks of identical linear systems[END_REF]; [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF]]. However, most physical systems are often not identical and frequently they are nonlinear in nature. The behavior of dynamical networks with heterogeneous nodes is much more complicated than the identical-node case. Usually, no common equilibrium for all nodes exists even if each isolated system has an equilibrium, the same for other invariant solutions, which can be destroyed or created by synchronization protocols.

The study of synchronization of dynamical networks with heterogeneous nodes is complicated and very few results have been reported by now. Some attempts have been made recently to propose output synchronization of heterogeneous systems and most of them employ the internal model principle [START_REF] Wieland | An internal model principle is necessary and sufficient for linear output synchronization[END_REF]; [START_REF] Isidori | Robust output synchronization of a network of heterogeneous nonlinear agents via nonlinear regulation theory[END_REF]; De Persis & Jayawardhana (2014); [START_REF] Liu | Robust decentralized output regulation with single or multiple reference signals for uncertain heterogeneous systems[END_REF]; [START_REF] Bidram | Synchronization of nonlinear heterogeneous cooperative systems using input-output feedback linearization[END_REF]]. The main idea is to assign each agent a (identical) local reference generator. The control algorithm then consists of two layers: a protocol, synchronizing the (identical) reference generators and local model-matching controllers, synchronizing the agents to their generators. Since agents share the same internal model, global network information is needed to implement the distributed controllers. Some other results are also available in the literature for particular classes of systems. For example, in [START_REF] Ahmed | Robust synchronization for multistable systems[END_REF], the authors have proposed robust synchronization for homogeneous/heterogeneous multistable systems. However, the systems are assumed to admit a decomposition without cycles (neither homoclinic nor heteroclinic orbits). Recently, the results of [START_REF] Ahmed | Robust synchronization for multistable systems[END_REF] have been applied to a multi-stable oscillator model [Ahmed et al. (2019)].

The goal of this work is to address the issue of synchronization of heterogeneous nonlinear systems using output feedback only, and an additional sub-goal is to have an oscillatory behavior in the synchronized state. Since many engineering systems have relative degree 2 or higher (e.g., pendulum systems, oscillators, robot manipulators, DC motor), the particular focus is put on this class of systems. Studying heterogeneous systems in general setting, we will assume that neither an equilibrium for each isolated node nor a synchronization manifold exists, so to synchronize them it is necessary to apply a feedback transformation [START_REF] Khalil | Nonlinear control[END_REF]] that projects all subsystems to a common (not necessarily identical) dynamics that can be synchronized next (the internal model principle). In this paper, the Brockett oscillator model is selected for this purpose. This is motivated by a global synchronization control recently proposed for such systems in Ahmed et al. (2019). Then higher order sliding mode (HOSM) observer is applied to estimate the unmeasurable states and perturbations using the idea presented in [START_REF] Fridman | Higher-order sliding-mode observer for state estimation and input reconstruction in nonlinear systems[END_REF]. In short, the main idea is to compensate the nonlinearities of individual systems followed by a nonlinear injection converting some parts of the systems into the Brockett oscillator form. The only restriction is that the individual systems should have relative degree 2 or higher (see Appendix for the definition), but most popular nonlinear benchmarks satisfy this criteria provided that the output signal is properly selected.

In this work, we have considered the output oscillatory synchronization of heterogeneous nonlinear systems. Potential applications for such a kind of synchronization include grid synchronization and load-sharing among DC/AC Inverters, to name a few. Grid current controllers working in synchronous reference frame require the phase of the grid voltage signal. This signal cannot be measured but can only be estimated.

Using the idea of master-slave output synchronization of oscillatory system, [Ahmed et al. (2019); [START_REF] Pay | Modeling and tuning of circular limit cycle oscillator fll with pre-loop filter[END_REF]; Ahmed et al. (2019)] estimated the phase of the grid voltage signal. In this approach, the grid voltage is considered as the output of the master oscillator while the slave oscillator tracks the master oscillator to estimate the parameters with oscillating dynamics. Similarly in the context of multiple oscillators, Virtual Oscillator Control (VOC) [START_REF] Sinha | Uncovering droop control laws embedded within the nonlinear dynamics of van der pol oscillators[END_REF]] for load sharing among multiple inverters can be considered as another potential application of the class of synchronization studied in this paper.

Main contribution: this paper studies the global output synchronization problem of nonlinear SISO (affine in control) heterogeneous multi-agent systems in a general setting. The proposed distributed synchronizing control law does not require any global network (or leader) information and uses the coupling with neighboring agents only. It can be applied to a network of systems having different orders. Moreover, precise information about the system parameters, uncertainties/disturbances etc., are also not needed as they are locally estimated using HOSM observer.

A preliminary conference version of this article has been presented in [START_REF] Ahmed | Oscillatory global output synchronization of nonidentical nonlinear systems[END_REF] where the authors consider only the case when all the agents have same relative degree. In this article, the results are extended for systems having different relative degrees and dimensions. All the omitted technical proofs of [START_REF] Ahmed | Oscillatory global output synchronization of nonidentical nonlinear systems[END_REF] are included in the current manuscript. Moreover, experimental results are provided to demonstrate the feasibility of the proposed controllers for real-time application. The results presented in this work use feedback transformation to project the subsystems dynamics to a common dynamics of Brockett oscillator [Ahmed et al. (2019)]. The results of Ahmed et al. (2019) are only applicable to the Brockett oscillator model which is of second order. However, the current work considers the global output synchronization problem of nonlinear SISO (affine in control) heterogeneous multi-agent systems in a general setting. This is a significant development with respect to Ahmed et al. (2019).

The rest of the article is organized as follows: Section 2 gives the problem statement followed by the synchronizing control design in Section 3. In Section 4, simulation and experimental studies are given, and finally Section 5 concludes this article. Preliminaries on relative degree and a summary of the result of Ahmed et al. (2019) can be found in the Appendix.

Problem Statement

The following family of nonlinear SISO systems (affine in control) is considered in this work for i = 1, N = 1, . . . , N with N > 1:

ẋi = f i (x i ) + g i (x i )u i , y i = h i (x i ), (1)
where x i ∈ R ni is the state, u i ∈ R (u i : R + → R is locally essentially bounded and measurable signal) is the input, y i ∈ R is the output; f i : R ni → R ni , h i : R ni → R and g i : R ni → R ni are sufficiently smooth functions. Denote the common state vector of (1) as

x = [x T 1 , . . . , x T N ] T ∈ R n with n = N i=1 n i , y = [y 1 , .
. . , y N ] T ∈ R N as the common output, and u = [u 1 , . . . , u N ] T ∈ R N as the common input. The relative degree condition (see Appendix) imposed on system (1) is summarized by the following assumptions: Assumption 1. For all i = 1, N , the systems in (1) have global uniform relative degree r i ∈ [2, n i ] and a globally defined normal form [START_REF] Marino | Nonlinear control design: geometric, adaptive and robust[END_REF]].

Under this assumption, for each subsystem in (1) there is a diffeomorphic transformation of coordinates T i : R n → R n such that [START_REF] Marino | Nonlinear control design: geometric, adaptive and robust[END_REF]; [START_REF] Khalil | Nonlinear control[END_REF]]:

η i ξ i = T i (x i ),
where ξ i ∈ R ri and η i ∈ R ni-ri are new components of the state, and for all i = 1, N the i th subsystem of (1) can be represented in the normal form:

ηi = ϕ i (η i , ξ i ), (2) ξi = A ri ξ i + b ri [α i (ξ i ) + β i (ξ i )u i ],
(3)

y i = c ri ξ i ,
where ϕ i : R ni → R ni-ri , α i : R ri → R and β i : R ri → R are smooth functions, β i is separated from zero, and

A ri =       0 1 0 . . . 0 0 0 0 1 . . . 0 0 . . . . . . . . . . . . 0 0 0 . . . 0 1 0 0 0 . . . 0 0       , b ri =       0 0 . . . 0 1       , c ri = 1 0 . . . 0
are in the canonical form. The subsystem (2) is called the zero dynamics of i th subsystem in (1), which we assume to be robustly stable:

Assumption 2. For all i = 1, N , the systems in (2) are input-to-state stable (ISS) with respect to the inputs ξ i [START_REF] Sontag | Smooth stabilization implies coprime factorization[END_REF]; [START_REF] Dashkovskiy | Input to state stability and allied system properties[END_REF]; [START_REF] Angeli | Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF]].

Concerning the definitions of ISS property used in this work, we will not distinguish ISS with respect to a set in the conventional sense [START_REF] Dashkovskiy | Input to state stability and allied system properties[END_REF]] or for a multistable system [START_REF] Angeli | Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF]], the only property we need here is the boundedness of the variables η i for bounded ξ i . More detailed analysis of the possible asymptotic behavior in (2) for the latter scenario is presented in [START_REF] Forni | Input-to-state stability for cascade systems with decomposable invariant sets[END_REF].

Then the synchronization problem consists in finding a control u such that the members of the family (1) perform synchronous (in phase) oscillations. Since the states of the subsystems in (1) may have different dimensions n i , a state synchronization error x i -x j cannot be defined in general (i.e. the states of the subsystems in (1) cannot follow their neighbors), but an output synchronization can be formulated:

Definition 1. The family (1) exhibits a global output synchronization if lim t→∞ (y i (t) -y j (t)) = 0, ∀i, j = 1, N for any initial conditions x i (0) ∈ R ni , i = 1, N .
Note that under Assumption 1 an additional requirement can be imposed on synchronization of derivatives:

lim t→∞ { ẏi (t) -ẏj (t)} = 0 ∀i, j = 1, N ,
and an auxiliary restriction for synchronization is

y i (t) = constant ∀i = 1, N ,
i.e. the systems perform a kind of oscillations in the synchronous mode.

An output feedback controller has to be designed to achieve the global output synchronization for (1).

Synchronization control design

The idea of this work is to design a feedback controller that will convert a part of subsystems (3) into the form of the Brockett oscillator [START_REF] Brockett | Synchronization without periodicity[END_REF]] through nonlinearity injection. Then global synchronization results can be easily obtained using the control proposed in Ahmed et al. ( 2019) (a summary is given in Appendix). However, this controller requires all components of the state vector to be available, which limits its implementation. Therefore, to overcome this difficulty, a high-order sliding-mode observer is used.

To simplify the presentation of the forthcoming synchronization protocol design, let us assume that

u i + d i = α i (ξ i ) + β i (ξ i )u i ,
where d i ∈ R is a new disturbance signal in (3) for each i = 1, N (since β i is not singular, such a representation always exists), and it may represent a parametric mismatch in the functions α i (ξ i ) and β i (ξ i ), but mainly it is introduced to model the estimation errors for the unmeasured state components ξ i . In such a case, if the further proposed observer for ξ i ensures its estimation, d i is always bounded and converging to zero or its vicinity. Finite-time convergence of d i using observer structure similar to the proposed one can also be found in the literature, cf. [START_REF] Ríos | An adaptive sliding-mode observer for a class of uncertain nonlinear systems[END_REF]. These useful properties are reflected in the following assumption: Assumption 3. For all i = 1, N , the unknown input d i : R + → R is continuously differentiable for almost all t ≥ 0, and there is a constant 0

< ν + < +∞ such that ess sup t≥0 | ḋi (t)| ≤ ν + .
The condition imposed by Assumption 3 is also satisfied by many engineering systems, e.g. stepper motor [START_REF] Defoort | A third-order sliding-mode controller for a stepper motor[END_REF], Assumption 1], DC/AC power inverter [START_REF] Gadelovits | Uncertainty and disturbance estimator based controller equipped with a multipletime-delayed filter to improve the voltage quality of inverters[END_REF], Eq. ( 2)], hydraulic actuators [START_REF] Ruderman | Virtual sensing of load forces in hydraulic actuators using second-and higher-order sliding modes[END_REF], Eq. ( 12)], to name a few. In many applications the disturbances are harmonic signals (coming from vibration of a part of the plant), Assumption 3 is also valid in those cases.

Observer design

Following the ideas presented in [START_REF] Fridman | Higher-order sliding-mode observer for state estimation and input reconstruction in nonlinear systems[END_REF] and [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF], let us first consider decoupling the input and disturbance. For that purpose, let us consider first, for all i = 1, N a Luenberger observer [START_REF] Luenberger | Observing the state of a linear system[END_REF]] for (3):

ζi = A ri ζ i + b ri u i + l i (y i -c ri ζ i ), (4)
where ζ i ∈ R ri is an auxiliary variable (an estimate of ξ i through a Luenberger observer), and l i ∈ R ri is the observer gain designed such that the matrix A ri -l i c ri is Hurwitz. The estimation error e i = ξ i -ζ i yields the following differential equation:

ėi = (A ri -l i c ri )e i + b ri d i ,
and to estimate also the unknown input d i , we consider an extended error vector:

e i = [e T i d i ] T .
Then from the available measurement output signal ψ i = c ri e i , we obtain:

˙ e i = A ri+1 e i -l i ψ i + b ri+1 ḋi ,
where [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF], the following high order sliding mode differentiator can be applied to estimate the error e i :

l i = [l T i 0] T . Based on
żi,1 = ν i,1 = -λ i,1 |z i,1 -ψ i | r i r i +1 sign(z i,1 -ψ i ) +z i,2 -l i,1 ψ i , żi,j = ν i,j = -λ i,j |z i,j -ν i,j-1 | r i -j+1 r i -j+2 sign(z i,j -ν i,j-1 ) +z i,j+1 -l i,j ψ i , j = 2, r i , (5) żi,ri+1 = -λ i,ri+1 sign(z i,ri+1 -ν i,ri ),
where λ i = [λ i,1 . . . λ i,ri+1 ] T ∈ R ri+1 is the vector of tuning parameters. The system (1) equipped with the observer (4), ( 5) is discontinuous due to the presence of sign function. The classical theory of differential equations is now not applicable since Lipschitz assumptions are employed to guarantee the existence of unique solutions. The solutions of the system (1) equipped with the observer (4), (5) are to be understood in the Filippov sense [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF]]. The solution concept proposed by Filippov for a differential equation with discontinuous right-hand-side is constructed as the "average" of the solutions obtained from approaching the point of discontinuity from different directions.

Denote by

ξ i = ζ i +    z i,1 . . . z i,ri    , d i = z i,ri+1
the estimates of ξ i and d i , respectively, provided by the observers (4) and ( 5). Then the following result can be proven:

Proposition 2. Let assumptions 1 and 3 be satisfied, the matrices A ri -l i c ri be Hurwitz. Then there exist λ i ∈ R ri+1 (λ i,ri+1 > ν + for all i = 1, N ) and T i > 0 such that for the system in (1) (or in (3)) and the observer in (4), ( 5) for all t ≥ T i :

y (j) i (t) = ξ i,j (t), j = 1, r i , d i (t) = d i (t).
Proof. Under the introduced assumptions, the system (1) can be transformed to the form (2), (3). In addition, the disturbance d i can be introduced with a bounded derivative for almost all instants of time. Next, for a properly selected λ i the result follows Lemma 8 in [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF], where the finite-time convergence and boundedness of the estimation error for (5) was proven, while the linear observer (4) serves to decouple the external disturbance d i and the control u i , which appear in the same equation.

The designed observers (4), ( 5) only use local input-output information u i and y i for each node i = 1, N , thus the proposed estimator is completely decentralized.

Control design

Once we have the estimates of ξ i and d i for all i = 1, N , i.e. the estimates for the states and the disturbances of the system (3), we are in position to design the synchronization control law.

The relative degree 2 case

First, assume that r i = 2. Then the following synchronizing control law can be proposed for all i = 1, N :

u i = -d i part 1 -ξ i,1 -b i ξ i,2 ξ 2 i,1 + ξ 2 i,2 -1 part 2 (6) + a i k ξ i-1,2 -2 ξ i,2 + ξ i+1,2 part 3
, where a i > 0, b i > 0 and k i > 0 are tuning parameters. The control law (6) has three parts: part 1 annihilates the nonlinearity of the original system (d i is dependent on α i and β i ); while part 2 injects additional nonlinearities to convert the system (3) into the form of the Brockett oscillator; finally, part 3 guarantees synchronization, since it contains the information of the left and right neighbors in the form of (17) (given in Appendix). In (6), part 1 and part 2 use only local information about the estimates calculated into the node ( ξ i and d i ), and only part 3 is based on signals sent over the network in (1). Thus, the control ( 6) is also decentralized, as in the observer (4), (5), and just the variables ξ i,2 (t) have to be communicated.

Theorem 3. Let assumptions 1, 2 and 3 be satisfied, the matrices A ri -l i c ri be Hurwitz, r i = 2 and λ i ∈ R ri+1 for all i = 1, N be selected as in Proposition 2. Consider the system (1) with the observers (4), ( 5) and the synchronizing feedback control (6). If there is an index 1 ≤ i ≤ N such that 2a i k < b i , then all trajectories in the closed-loop system are bounded, and for almost all initial conditions they converge to the largest invariant set, where the following restrictions are satisfied for all i = 1, N :

ẏi-1 + ẏi+1 = (2 + b i a i k (y 2 i + ẏ2 i -1)) ẏi , (7)
y 2 i + ẏ2 i = constant = 0, (y i -y i+1 ) 2 + ( ẏi -ẏi+1 ) 2 = constant .
Proof. Since all conditions of Proposition 2 are satisfied, the observer (4), ( 5) is finitetime convergent and for t ≥ T = max 1≤i≤N T i the estimates y (j) i (t) = ξ i,j (t) for j = 1, r i and d i (t) = d i (t) are valid for all i = 1, N . Therefore, due to the structure of the control (6), for t ≥ T the family (3) become a family of Brockett oscillators which implies the synchronization result stated in Theorem 9 thanks to part 3 of the controller (6). Let us prove that there is no finite-time escape of trajectories on the interval [0, T ]. The convergence of observers is independent of the form of control and the estimation errors ξ i -ξ i and d i -d i stay bounded for all t ≥ 0. Then the system (3) with the control ( 6) can be presented in the following form:

ξi,1 = ξ i,2 ξi,2 = 3 s=1 δ i,s -ξ i,1 + b i ξ i,2 , where δ i,1 = a i k(ξ i-1,2 -2ξ i,2 + ξ i+1,2 ), δ i,2 = ξ i,1 -ξ i,1 + d i -d i -b i (ξ i,2 -ξ i,2 ) +a i k{ ξ i-1,2 -ξ i-1,2 -2( ξ i,2 -ξ i,2 ) + ξ i+1,2 -ξ i+1,2 }, δ i,3 = b i (ξ i,2 -ξ i,2 )| ξ i | 2 -b i ξ i,2 | ξ i | 2
are auxiliary inputs, δ i,1 contains coupling signals, δ i,2 is bounded since it is composed by the observer estimation errors, and δ i,3 contains all nonlinear terms. Let us consider a Lyapunov function U i (ξ i ) = 0.5|ξ i | 2 for this subsystem, whose derivative takes the form (in the calculations below we will use convention for the indexes as N + 1 = 1):

Ui = ξ i,2 3 s=1 δ i,s + b i ξ 2 i,2 ≤ ξ i,2 3
s=1 δ i,s + b i U i ,
and the following upper estimates can be calculated using Young's inequality [START_REF] Khalil | Nonlinear control[END_REF]]:

ξ i,2 δ i,1 ≤ 0.5a i k(ξ 2 i-1,2 + ξ 2 i+1,2 ) ≤ a i k(U i-1 + U i+1 ), ξ i,2 δ i,2 ≤ 0.5(ξ 2 i,2 + δ 2 i,2 ) ≤ U i + 0.5δ 2 i,2 , ξ i,2 δ i,3 = b i (ξ i,2 -ξ i,2 )ξ i,2 | ξ i | 2 -b i ξ 2 i,2 | ξ i | 2 ≤ b i | ξ i,2 -ξ i,2 ||ξ i,2 || ξ i | 2 -b i ξ 2 i,2 | ξ i | 2 = b i | ξ i,2 -ξ i,2 | -|ξ i,2 | |ξ i,2 || ξ i | 2 ≤ b i 0 | ξ i,2 -ξ i,2 | ≤ |ξ i,2 | | ξ i,2 -ξ i,2 | 2 | ξ i | 2 | ξ i,2 -ξ i,2 | > |ξ i,2 | .
Note that

| ξ i | 2 = | ξ i -ξ i + ξ i | 2 ≤ 2| ξ i -ξ i | 2 + 2|ξ i | 2 ,
then finally we obtain:

ξ i,2 δ i,3 ≤ 2b i | ξ i,2 -ξ i,2 | 2 [| ξ i -ξ i | 2 + 2U i ] Therefore Ui ≤ a i k(U i-1 + U i+1 ) + (4b i | ξ i,2 -ξ i,2 | 2 + b i + 1)U i + 0.5δ 2 i,2 + 2b i | ξ i,2 -ξ i,2 | 2 | ξ i -ξ i | 2
There are constants i > 0 and i > 0 (dependent on initial conditions) such that

4b i | ξ i,2 -ξ i,2 | 2 + b i + 1 ≤ i , | ξ i,2 -ξ i,2 | 2 | ξ i -ξ i | 2 ≤ i , consequently, Ui < a i k(U i-1 + U i+1 ) + i U i + 0.5δ 2 i,2 + 2b i i
and considering a common Lyapunov function for (3)

U (ξ 1 , . . . , ξ N ) = N i=1 U i (ξ i )
we obtain that

U < N i=1 [a i k(U i-1 + U i+1 ) + i U i + 0.5δ 2 i,2 + 2b i i ] = σU + 0.5 N i=1 δ 2 i,2 + 4b i i ,
where σ = max 1≤i≤N i + max 1≤i≤N a i k. Since δ i,2 is bounded for all i = 1, N according to Proposition 2, then the above inequality implies that U admits an upper estimate exponential in time and the system is forward complete (Angeli & Sontag, 1999, Corollary 2.11).

Remark 4. In control (6), it is assumed that the agents are connected through a N -cycle graph1 [START_REF] Pemmaraju | Cycles, stars, and wheels[END_REF]] similar to Ahmed et al. (2019). Other network topology may also be considered.

The higher relative degree case

Now, consider the general case with r i ≥ 2, then the parts 2 and 3 of the control ( 6) form a reference signal ξ d i,3 for the variable ξ i,3 :

ξ d i,3 = -ξ i,1 -b i ξ i,2 ξ 2 i,1 + ξ 2 i,2 -1 +a i k ξ i-1,2 -2 ξ i,2 + ξ i+1,2 ,
where the parameters a i > 0, b i > 0 and k i > 0 save their meaning, and next this reference signal has to be propagated over chain of integrators, and the part 1 of the control ( 6) has to be applied on the last step to annihilate

d i . Let us denote ξi,s = ξ i,s -ξ d i,s
as the error or realization of a desired signal ξ d i,s by ξ i,s for s = 3, r i . For t ≥ T = max 1≤i≤N T i (the time when the estimates y

(j) i (t) = ξ i,j (t) for j = 1, r i and d i (t) = d i (t)
are valid for all i = 1, N ) consider the Lyapunov function from Ahmed et al. (2019):

V (ξ) = N i=1 α i 4 (ξ 2 i,1 + ξ 2 i,2 -1) 2 + 1 2 N i=1 (ξ i,1 -ξ i+1,1 ) 2 + (ξ i,2 -ξ i+1,2 ) 2 ,
where α i = bi kai , then

V = N i=1 -b i α -0.5 i z 2 (ξ, i) -a i z(ξ, i) ξi,3 and z(ξ, i) = ξ i-1,2 +ξ i+1,2 -ξ i,2 {2+α i (ξ 2 i,1 +ξ 2 i,2 -1)}.
Next, applying the backstepping approach [START_REF] Krstic | Nonlinear and adaptive control design[END_REF]], we can derive:

ξ d i,4 = ξd i,3 -κ ξi,3 + a i z(ξ, i), ξ d i,s+1 = ξd i,s -κ ξi,s -ξi,s-1 , s = 4, r i -1, u i = -d i + ξd i,ri -κ ξi,ri -ξi,ri-1 . ( 8
)
With such a control algorithm, the Lyapunov function

W (ξ) = V (ξ) + 0.5 N i=1 ri s=3 ξ2 i,s
admits the time derivative

Ẇ = - N i=1 b i α -0.5 i z 2 (ξ, i) + κ ri s=3 ξ2 i,s ≤ 0.
In order to implement the control u i given in ( 8) it is necessary to calculate the deriva-

tives d j ξi-1,2(t)
dt j and d j ξi+1,2(t) dt j

for j = 1, r i -3. Therefore, they have to be either transmitted by the network or reconstructed on-line by additional differentiators of a kind presented in (5). The latter solution is more preferable for a distributed synchronization protocol, i.e. for a = -1, +1:

˙ a i,0 = ν i,0 = -µ i,0 | a i,0 -ξ i+a,2 | r i -3 r i -2 sign( a i,0 -ξ i+a,2 ) + a i,1 , ˙ a i,j = ν i,j = -µ i,j | a i,j -ν i,j-1 | r i -j-3 r i -j-2 sign( a i,j -ν i,j-1 ) + a i,j+1 , j = 1, r i -4, (9)
˙ a i,ri-3 = -µ i,ri-3 sign( a i,ri-3 -ν i,ri ),
then by Lemma 8 in [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF] for a proper selection of

µ i = [µ i,0 , . . . , µ i,ri-3 ] T > 0 there is T i > 0 such that a i,j (t) = ξ (j)
i+a,2 (t) ∀t ≥ T i for all j = 0, r i -3 and a = -1, +1. Note that from (5), ξ i+a,2 has only r i+a -1 continuous derivatives, consequently

r i ≤ r i+a + 2 ∀i = 1, N , a = -1, +1. (10)
The following result has been proven:

Theorem 5. Let assumptions 1, 2 and 3 be satisfied, the matrices A ri -l i c ri be Hurwitz, r i ≥ 2 under (10), and λ i ∈ R ri+1 , µ i ∈ R ri-2 be properly selected for all i = 1, N .

Consider the system (1) with the observers (4), ( 5), ( 9) and the synchronizing feedback control (8). If there is an index 1 ≤ i ≤ N such that 2a i k < b i , then all trajectories in the closed-loop system are bounded, and for almost all initial conditions they converge to the largest invariant set, where the restrictions (7) are satisfied.

Proof. For t ≥ max 1≤i≤N T i + max 1≤i≤N T i all arguments presented above are valid, and the forward completeness can be proven similarly to the case of Theorem 3.

Remark 6. Theorem 5 doesn't assume that the relative degrees of the individual systems in a network are same. It shows that global oscillatory synchronization is possible even when the individual systems in a network have different relative degrees (see Sec. 4.1 for numerical example).

Instead of the conventional backstepping [START_REF] Krstic | Nonlinear and adaptive control design[END_REF]], the command filtered backstepping [START_REF] Farrell | Command filtered backstepping[END_REF]] can also be applied, and since it does not need derivatives of the virtual controls, the observers (9) can be avoided in this case. A drawback in this case is that exact synchronization becomes impossible.

Results and discussions

Numerical simulation results

In this Section, numerical simulations are considered to show the effectiveness of the proposed synchronization scheme. For this purpose, let us consider the synchronization of three heterogeneous nonlinear oscillatory systems composed of two Van der Pol oscillators and one oscillator with dynamics like Van der Pol oscillator. The dynamics of Van der Pol oscillator is given by

ẋ1i = x 2i , i = 1, 2 ẋ2i = -x 1i + ς i {1 -x 1i } x 2i + u i y i = x 1i (11)
where i represents the oscillator number, x 1i , x 2i are the states, y i is the output, ς i is the model parameter and u i is the control input of oscillator i. The dynamics of the third system is given by

ẋ1i = x 2i , ẋ2i = x 3i ẋ3i = -x 2i + ς i 1 -(x 2i ) 2 x 3i + u i y i = x 1i (12)
where i = 3 represents the oscillator number and the rest of the variables retain the same meaning as used for eq. ( 11). Although (12) represents one oscillator but we have used the notation i for uniformity. The parameters ς i are selected as ς i = 0.1 × i for simulation. So, here the synchronization of systems with different relative degree and order is considered. With the output y i , the system (11) has relative degree r i = n = 2, while with the output y i , the system (12) has relative degree r i = n = 3. Then the family ( 11), ( 12) can be written in the form (3) as for i = 1, 2

˙ x 1i x 2i = 0 1 0 0 x 1i x 2i + 0 1      u i + -x 1i + ς i 1 -(x 1i ) 2 x 2i d i      (13)
and

Figure 1. Evolution of the oscillator states with control ( 6) and ( 8).

  x 1i x 2i x 3i   =   0 1 0 0 0 1 0 0 0     x 1i x 2i x 3i   +   0 0 1       u i + -x 2i + ς i 1 -(x 2i ) 2 x 3i di     (14)
Since the states of Van der Pol oscillators are bounded, then the disturbances d i in ( 13) and ( 14) satisfy Assumption 3. Then, for systems ( 13) and ( 14), the design of observers ( 4) and ( 5) is straightforward. For the Van der Pol oscillators, the observer parameters are selected as L = 20.2 102 T , λ 1 = 3M , λ 2 = 1.5M and λ 3 = 1.1M

for M = 50 while for Van der Pol like oscillator, the observer parameters are selected as L = 30.6 312.1 1061.1 T , λ 1 = 5M , λ 2 = 3M , λ 3 = 1.5M and λ 4 = 1.1M

for M = 70. With the estimated states, the synchronizing controller (6) can be easily designed for Van der Pol oscillators represented by (11). System ( 12) is of relative degree 3. So, the controller (8) is to be designed for ( 14). The parameters of the controllers are selected as:

a 1 k 1 = 0.4, a 2 k 2 = 0.8, a 3 k 3 = 1.2, b 1 = 1, b 2 = 2, b 3 = 3 and κ = 10.
With these values of the parameters, the conditions of Theorem 5 are satisfied. The evolution of the oscillator states with controls ( 6) and ( 8) is given in Fig. 1. For this simulation, the oscillators were initialized as (2, -2), (0, 2) and (-2, 0, 2). From Fig. 1, it is clear that the control (6) successfully converted the Van der Pol oscillators in a finite-time to Brockett oscillators through nonlinearity injection while the control (8) converted part of system ( 14) into the form of Brockett oscillator. Then the result of global synchronization follows directly from Theorem 5. The phase plot of the individual systems and the the graph of (y i -y i+1 ) 2 + ( ẏi -ẏi+1 ) 2 as mentioned in Theorem 3 are given in figures 2 and 3 respectively. Numerical simulation results demonstrate the effectiveness the of the proposed synchronization scheme. The impact of perturbations are studied in detail in the experimental study section and skipped here to avoid repetition.

Experimental Study

To verify the practical feasibility and performance of the synchronization scheme developed in Section 3, let us consider the synchronization of a network of the Van der Pol oscillators with N = 3 as given by eq. ( 13). To implement the proposed synchronizing control, a dSPACE 1104 controller board was used. The control algorithms were implemented using Simulink. The solver was the fixed-step fourth-order Runge-Kutta and the time-step was 0.2 milliseconds. An overview of the experimental setup and the circuit diagram of an autonomous Van der Pol oscillator can be seen in Fig. 4. For experiments ς i = 0.1 was selected in (11), and the circuit parameters of Fig. 4 are chosen accordingly (the values of the circuit parameters are skipped here for brevity). We have selected identical values of circuit parameters for the experimental realization of the network of Van der Pol oscillators. However, nonidealities and nonlinearities of the practical devices, make the oscillators heterogeneous in a practical settings. For example the resistors that are used in this experiment have a tolerance limit of ±10%. So, even if we take the resistors of the same nominal values for individual oscillators, the final values are different due to the tolerance limit.

Van der Pol oscillator is a benchmark model for second order nonlinear system. Global synchronization of model ( 13) is a very interesting problem. However, it is difficult to prove analytically the global synchronization of model ( 13) with heterogeneous nodes. As an alternative way, it is possible to globally synchronize family (13) by transforming individual oscillators into Brockett form as mentioned in Section 3.

For system (13), the design of observers (4) and ( 5 To check the robustness of the proposed control strategy, two tests were done. In the first case, successive perturbations (in the form of additional DC voltage inputs) were added to x 23 . The evolution of the state variables and the control inputs in this case is given in Fig. 7. It is evident from Fig. 7 that the proposed control strategy is very robust with respect to perturbations. Although the oscillators lose synchrony at the beginning, they return to the synchronous state in a very short period of time due to global nature of the control (6). To further check the robustness, large amplitude perturbations were added to state x 12 (note that the control influences directly x 2i , but not x 1i , which is why it is an interesting case to consider). The results in this case are shown in Fig. 8. The experimental results again demonstrate the robustness of the proposed control. However, in this case, the oscillators need more time to synchronize because of large amplitude of perturbations and also due to the fact that x 1i do not depend directly on u i .

Conclusion

This paper studied the problem of global synchronization of nonlinear systems with relative degree 2 and higher using an output feedback. The nonlinear systems were first converted to a normal canonical form. Then higher order sliding mode observers were used to reconstruct the states and the perturbations in a finite time. Using these information, individual systems were projected to Brockett oscillator dynamics through dynamic output feedback control. The synchronization result was obtained by applying the results of Ahmed et al. (2019). Experimental study demonstrated the effectiveness of our method using a network of heterogeneous Van der Pol oscillators.

B. Synchronization of Brockett oscillators

The following family of Brockett oscillators [Brockett (2013)] is considered in this section for some N > 1:

ẋ1i = x 2i , ẋ2i = a i u i -x 1i -b i x 2i |x i | 2 -1 , i = 1, N , (16)
where a i , b i > 0 are the parameters of an individual oscillator, the state x i = [x 1i x 2i ] T ∈ R 2 and the control u i ∈ R (u i : R + → R is locally essentially bounded and measurable signal). Denote the common state vector of ( 16) by x = [x T 1 , . . . , x T N ] T ∈ R 2N and the common input by u = [u 1 , . . . , u N ] T ∈ R N .

The following synchronizing control is selected for family ( 16):

u = kM     
x 21 . . .

x 2(N -1) x 2N      , ( 17
)
where k > 0 is the coupling strength and

M =       -2 1 0 • • • 1 1 -2 1 • • • 0 0 1 -2 • • • 0 . . . . . . . . . 0 1 • • • 0 1 -2      
is the interconnection matrix. From a graph theory point of view, the oscillators are connected through a N -cycle graph [START_REF] Pemmaraju | Cycles, stars, and wheels[END_REF]] (each oscillator needs only the information of its left and right neighbors). Define the synchronization error among the various states of the oscillators as e 2i-1 = x 1i -x 1(i+1) , ė2i-1 = x 2i -x 2(i+1) = e 2i and e 2N -1 = x 1N -x 11 , ė2N-1 = x 2N -x 21 = e 2N . Then the main results of Ahmed et al. (2019) can be summarized as below:

Proposition 8. [Ahmed et al. (2019)] For any k > 0 in the system ( 16), (17) all trajectories are bounded and converge to the largest invariant set in

Ω ∞ = x ∈ R 2N : |x i | = const, e 2 2i-1 + e 2 2i = const, x 2(i-1) + x 2(i+1) = (2 + b i a i k (|x i | 2 -1))x 2i , i = 1, N .
Theorem 9. [Ahmed et al. (2019)] For any k > 0, if there is an index 1 ≤ i ≤ N such that 2a i k < b i , then in the system ( 16), (17) all trajectories are bounded and almost

Figure 2 .

 2 Figure 2. Phase plot i.e. y i vs ẏi of the three systems.

Figure 3 .

 3 Figure 3. Graph of (y i -y i+1 ) 2 + ( ẏi -ẏi+1 ) 2 for the three systems.

  ) are straightforward. The parameters of observers used in the experiments are L = 20.2 102 T , λ 1 = 3M ,

Figure 4 .

 4 Figure 4. Schematic overview of the practical implementation and circuit diagram of an autonomous Van der Pol oscillator.

Figure 5 .

 5 Figure 5. Convergence of the observer. A zoomed version of the the measured state variable (after analog to digital conversion) can be seen in the inset.

Figure 6 .

 6 Figure 6. Evolution of the oscillator states. Solid linesx 2i and dashed linesx 1i .

Figure 7 .

 7 Figure 7. Evolution of the oscillator states and control inputs in the case of perturbations added to x 23 . Solid linesx 2i and dashed linesu i . Arrow indicates the time when perturbations were added.

Figure 8 .

 8 Figure 8. Evolution of the oscillator states and control inputs in the case of large amplitude perturbations added to x 12 . Solid linesx 2i and dashed linesx 1i . Arrow indicates the time when the first perturbation was added.

A cycle graph C N is a graph on N nodes containing a single cycle through all nodes, or in other words, N number of vertices connected in a closed chain.

Funding This work was partially supported by Consejo Nacional de Ciencia y Tecnologia, Mexico [grant no. 282013]; Government Council on Grants, Russian Federation [Grant 08-08]; Ministry of Education and Science of the Russian Federation [Project 14.Z50.31.0031]; Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica (UNAM), Mexico 115419; HoTSMoCE INRIA associate team program.

Appendix

A. Relative Degree

Consider the following nonlinear system

where x ∈ R n is the state vector, u ∈ R is the input, y ∈ R is the output variable of the system, f and g are smooth vector fields. A vector field is said to be forward complete if all solutions to ẋ = f (x) are defined for all t ≥ 0 [ [START_REF] Khalil | Nonlinear control[END_REF]].

Definition 7. (Global Uniform Relative Degree [START_REF] Marino | Nonlinear control design: geometric, adaptive and robust[END_REF]]) The global uniform relative degree r of ( 15) is defined as the integer such that

We say that r = ∞ if

all of them converge to the largest invariant set in

In the set Ω ∞ we have for all i = 1, N :

for some r i ∈ R + and ρ i ∈ R + , and

for

kai and some c i ∈ R. Corollary 10. [Ahmed et al. (2019)] Let all conditions of Theorem 9 be satisfied, and all solutions of the following equations

with r i = 1 admit the restriction:

for some 1 ≤ i ≤ N . Then for almost all initial conditions the system ( 16), ( 17) is synchronized.