E Lapina

P Oumaziz

R Bouclier

Immersed boundary-conformal isogeometric LaTIn method for multiple non-linear interfaces

Keywords: Splines, Non-conformal coupling, Robin conditions, Multiscale, Contact, Delamination

Computational micromechanics appears of the utmost importance, especially in the current context of digital twins in mechanics of materials. The objective here is to develop an efficient solver for the simulation of geometrically complex composite microstructures involving numerous inclusions connected with the matrix through various non-linear interface behaviors. To do so, we resort to IsoGeometric Analysis, which provides increased per-degree-of-freedom accuracy, and leverage the recently introduced immersed boundary-conformal method to retrieve conformal matrix/inclusion interfaces through the construction of conformal layers from it. Then, the approach is enhanced with the Large Time INcremental method that allows to separate the non-linear interface equations from those related to the subdomains, the latter being all linear and subdomain-wise independent. It results in an immersed hybrid mixed higher-order numerical scheme that is naturally parallelizable between the different subdomains and that is flexible to treat any non-linear interface behavior. The stabilization of the formulation occurs within the bulk equations where Nitsche couplings are performed. The accuracy and efficiency of the developed algorithm are demonstrated by solving a range of non-linear examples in 2D, including different numbers of inclusions in unilateral contact, frictional contact, and delamination with the matrix of the composite microstructure.

Introduction

With the advent of imaging tools on one side, and the desire for data assimilation and digital twins on the other side, computational micromechanics, that seeks to perform numerical simulation at the micro scale, definitely constitutes an emerging field nowadays in the mechanics of materials. Starting from images built, computational micromechanics enables to perform numerical simulation using the exact, imperfect geometry of the sample [START_REF] Legrain | An x-fem and level set computational approach for image-based modelling: application to homogenization[END_REF][START_REF] Verhoosel | Image-based goal-oriented adaptive isogeometric analysis with application to the micro-mechanical modeling of trabecular bone[END_REF][START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF][START_REF] Claus | Contact modeling from images using cut finite element solvers[END_REF]. Combining the methodology with in-situ experimental tests also opens the door for data assimilation at the micro scale [START_REF] Heinze | Experimental and numerical investigation of single pores for identification of effective metal foams properties[END_REF][START_REF] Wu | A combined experimental and numerical method to estimate the elastic modulus of single trabeculae[END_REF], by comparing measured fields coming, e.g., from digital image correlation [START_REF] Buljac | Digital volume correlation: review of progress and challenges[END_REF][START_REF] Rouwane | Architecture-driven digital volume correlation: application to the analysis of in-situ crushing of a polyurethane foam[END_REF], and simulated fields obtained from computational micromechanics. In other contexts, computational micromechanics offers the opportunity to feed artificial intelligence to create homogenized constitutive laws [START_REF] Fernández | Anisotropic hyperelastic constitutive models for finite deformations combining material theory and data-driven approaches with application to cubic lattice metamaterials[END_REF][START_REF] Masi | Multiscale modeling of inelastic materials with thermodynamics-based artificial neural networks (tann)[END_REF], or to perform virtual material design [START_REF] Herráez | Computational micromechanics evaluation of the effect of fibre shape on the transverse strength of unidirectional composites: an approach to virtual materials design[END_REF], i.e. to optimize the material microstructure to achieve a desired behavior. In line with this background, the purpose of this work is to develop an efficient and accurate numerical method able to compute multiscale and geometrically complex composite microstructures involving numerous inclusions connected with the matrix through various non-linear interface behaviors.

As a first step towards efficiency, the present work utilizes IsoGeometric Analysis (IGA). Introduced in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF], IGA proposes to use the higher-order and smooth spline bases encountered in Computer-Aided Design (CAD), e.g. the Non-Uniform Rational B-Splines (NURBS) or simpler B-splines [START_REF] Cohen | Discrete B-spline and subdivision techniques in computer aided geometric design and computer graphics[END_REF][START_REF] Piegl | The NURBS book. 2Nd[END_REF], not only for the representation of the geometry but also for the approximation of solution fields in numerical simulations (see, for instance, [START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF] for a recent review). While the original aim was to streamline the time-consuming process of creating mechanical models from CAD programs, the use of splines offers its own advantages from an analysis point of view. Indeed, spline functions can be C (p-1) regular between elements for a polynomial degree p, while Lagrange polynomials, which are used in the standard Finite Element Method (FEM), attain only a C 0 regularity at these locations. As a result, IGA can exhibit higher per-degree-of-freedom accuracy and robustness compared to the standard FEM for mechanical simulations [START_REF] Evans | n-Widths, sup-infs, and optimality ratios for the k-version of the isogeometric finite element method[END_REF], which makes this approach often seen as a High-Performance Computational (HPC) tool.

Nevertheless, modeling multiscale and geometrically complex objects, such as composite microstructures involving a great number of local inclusions, is challenging in IGA. Indeed, multivariate B-splines and NURBS come with a rigid tensor-product structure, therefore precluding the simple modeling of local behaviors. More advanced splines, such as locally refined splines [START_REF] Dokken | Polynomial splines over locally refined boxpartitions[END_REF], hierarchical splines [START_REF] Hennig | Bézier extraction and adaptive refinement of truncated hierarchical NURBS[END_REF][START_REF] Angella | Efficient extraction of hierarchical b-splines for local refinement and coarsening of isogeometric analysis[END_REF], T-splines and hierarchical T-splines [START_REF] Veiga | Analysis-suitable T-splines are dualcompatible[END_REF][START_REF] Evans | Hierarchical T-splines: Analysis-suitability, Bézier extraction, and application as an adaptive basis for isogeometric analysis[END_REF], have emerged to offer local mesh refinement, but the integration of arbitrary local models (e.g., inclusions) within spline patches (e.g., matrix) still remains an issue [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: Geometric details, local refinement, inclusion, fracture[END_REF][START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF][START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF][START_REF] Lapina | A fully non-invasive hybrid iga/fem scheme for the analysis of localized non-linear phenomena[END_REF]. From a general point of view, mainly two approaches can be followed. Firstly, one can try to invoke spline re-parametrization procedures [START_REF] Massarwi | Volumetric untrimming: Precise decomposition of trimmed trivariates into tensor products[END_REF], thereby leading to the splitting of the new geometry into several patches with C 0 regularity at the boundaries. This may entail a considerable modeling and computational effort which is often as complex and time-consuming as standard mesh generation and then, is opposed to the core idea of IGA. The second option consists in developing advanced numerical methods that allow to relax the need for geometrically conformal meshes. In this respect, it may appear appealing to resort to the family of immersed boundary methods that have been subjected to intensive studies in IGA these last years (see [START_REF] Rank | Geometric modeling, isogeometric analysis and the finite cell method[END_REF][START_REF] Schillinger | The finite cell method: A review in the context of higher-order structural analysis of CAD and image-based geometric models[END_REF] for the origin with the Finite-Cell Method (FCM), [START_REF] Verhoosel | Image-based goal-oriented adaptive isogeometric analysis with application to the micro-mechanical modeling of trabecular bone[END_REF][START_REF] Sai C Divi | Topology-preserving scan-based immersed isogeometric analysis[END_REF] for trabecular bones, [START_REF] Rouwane | Adjusting fictitious domain parameters for fairly priced image-based modeling: Application to the regularization of digital image correlation[END_REF] for cellular materials, [START_REF] Hoang | Skeleton-stabilized immersogeometric analysis for incompressible viscous flow problems[END_REF] for flow problems in porous media, [START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF] for composite materials, [START_REF] Casquero | The divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics[END_REF] for fluid-structure interactions in vesicles, and [START_REF] Patton | Cost-effective and accurate interlaminar stress modeling of composite kirchhoff plates via immersed isogeometric analysis and equilibrium[END_REF] for composite kirchhoff plates, to name a few). The idea is to simply use unfitted structured meshes for the interpolation of the mechanical fields, while the cut or trimmed domains are accurately captured by means of suitable quadrature rules [START_REF] Nagy | On the numerical integration of trimmed isogeometric elements[END_REF][START_REF] Kudela | Efficient and accurate numerical quadrature for immersed boundary method[END_REF][START_REF] Antolin | Isogeometric analysis on V-reps: first results[END_REF][START_REF] Legrain | Non-negative moment fitting quadrature rules for fictitious domain methods[END_REF][START_REF] Garhuom | Non-negative moment fitting quadrature for cut finite elements and cells undergoing large deformations[END_REF]. In this work, we propose to follow this second path to be completely free from the geometric complexity of the model.

The main task then lies in the formulation and implementation of a coupling scheme between the background spline mesh (describing, e.g., the composite matrix) and the multiple foreground inclusions. This appears far from trivial since the coupling interfaces are expected to be both non-conformal (the latter may cut in any way the elements of the background grid) and non-linear (such as incorporating contact or delamination). As a remedy, the first ingredient of our approach is to consider the recently introduced Immersed Boundary-Conformal Method (IBCM) [START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF] that has proved to be accurate and robust for perfect (linear) interfaces. It consists in a pragmatic strategy that transforms the initial interface, through the construction of conformal layers from it, into three different interfaces: the initial one between the matrix and the inclusion that becomes conformal, and two non-conformal ones that now lie within the matrix and inclusion. Therefore, the strategy leverages the geometric flexibility of the immersed methods with the advantages of conformal discretizations, which seems to us all the more relevant in our context. Indeed, (i) the solution and especially the stress fields can be properly described around the interfaces with controlled conformal discretizations, (ii) complex non-linear laws can be enforced in standard manners through the conformal interfaces, and (iii) the non-conformal aspect can be treated within the framework of linear elasticity, which is now well established in the field, in particular using Nitsche-type methods [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF][START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF][START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF]. Moreover, let us note that the difficulty of building conformal layers is generally mitigated when the method is applied to composite materials, as the local inclusions may have rather simple geometries (for instance, they merely consist of circular cylinders in case of fiber-reinforced composites).

As a result, the objective now is to extend the IBCM to account for multiple non-linear interfaces and to perform HPC. The treatment of this type of problem does not seem to be covered by the HPC methods developed in IGA, which are currently restricted mainly to the framework of linear partial differential equations (see, e.g., [START_REF] Calabro | Fast formation of isogeometric galerkin matrices by weighted quadrature[END_REF][START_REF] Mantzaflaris | Low rank tensor methods in galerkin-based isogeometric analysis[END_REF][START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF] for fast operator assembly and [START_REF] Hofreither | Robust multigrid for isogeometric analysis based on stable splittings of spline spaces[END_REF][START_REF] Hirschler | A dual domain decomposition algorithm for the analysis of non-conforming isogeometric kirchhoff-love shells[END_REF][START_REF] Micha L Bosy | A domain decomposition method for isogeometric multi-patch problems with inexact local solvers[END_REF][START_REF] Olof B Widlund | Bddc deluxe algorithms for twodimensional h (curl) isogeometric analysis[END_REF] for multigrid and domain decomposition solvers). Conversely, numerous works may be reported in the field of standard FEM to answer the issue, including the so-called Large Time INcremental (LaTIn) method introduced in [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation[END_REF]. Originally designed to consider material non-linearities such as plasticity and damage in a non-incremental approach, the method has quickly been applied to domain decomposition with non-linear interfaces, starting with unilateral and frictional contact [START_REF] Blanzé | Contact problems in the design of a superconducting quadrupole prototype[END_REF][START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF][START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF] and then going to delamination [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Violeau | Micromodel-based simulations for laminated composites[END_REF][START_REF] Fernandez | The effect of the mesh refinement on a multiscale domain decomposition method for the non-linear simulation of composite structures[END_REF]. More recently, the method has also been extended to the case of immersed interfaces for multiple contacts in FEM [START_REF] Claus | A stable and optimally convergent latin-cutfem algorithm for multiple unilateral contact problems[END_REF][START_REF] Claus | Contact modeling from images using cut finite element solvers[END_REF] by enhancing it with specific stabilization terms inspired by the CutFEM methodology [START_REF] Burman | Cutfem: discretizing geometry and partial differential equations[END_REF]. The principle of the LaTIn approach for domain decomposition is to separate the non-linear interface equations from those related to the subdomains, the latter being all linear and subdomain-wise independent. Therefore, the interest is that the subdomain equations can be solved efficiently in parallel, and the non-linear interface ones can be written locally in an easy fashion regardless of the treated non-linearity.

In this context, our idea is to develop an immersed boundary-conformal isogeometric LaTIn method in order to accurately and efficiently compute composite microstructures involving multiple inclusions. In our framework, the LaTIn method is therefore applied across conformal interfaces without any difficulty, which results in the local formulation of the non-linear interface equations at each interface Gauss point, and in the parallel solution of subdomain-wise Robin-type linear elastic problems. The stabilization of the formulation is moved within the bulk equations related to the matrix and the inclusions, where Nitsche couplings with minimal stabilization are performed [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF][START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF]. The implemented algorithm consists in an immersed, hybrid (bulk/interface), mixed (displacement/force,) higher-order, and smoother numerical scheme that appears accurate and scalable to treat multiple inclusions in unilateral contact, frictional contact, and delamination with the matrix of the composite microstructure.

The paper is organized as follows: after this introduction, the strong form of the reference multiscale coupling problem, as well as its discretization in IGA, are discussed in Section 2. Then, in Section 3, the proposed approach for dealing with both non-linear and non-conformal interfaces is presented, introducing the IBCM, stabilized Nitsche coupling, as well as the LaTIn approach. Next, a range of numerical experiments are carried out in 2D to assess the performance of our method with respect to accuracy and efficiency in Section 4, followed by conclusions and discussion in Section 5.

Reference mechanical coupling problem and discretization

This section establishes the context of the study and introduces the corresponding notations. The reference mechanical coupling problem is presented, before some remarks are added to highlight the difficulty of modeling multiple local behaviors within the IGA framework.

Governing equations

We undertake to solve a multi-domain problem characterized by a physical domain Ω ⊂ R d , d = 2 or 3 being the topological dimension of the domain. More precisely, since our target applications are composite materials such as fiber-reinforced ones, we consider the coupling of a global model (say the matrix) with multiple local models (say the fibers) through non-linear interfaces. For the sake of simplicity, let us take here only two subdomains. Ω is therefore divided into two disjoint, open and bounded subsets Ω m and Ω f joining at interface Γ, i.e. such that Ω = Ω m ∪ Ω f ∪ Γ and Ω m ∩ Ω f = ∅ (see Fig. 1 for illustration). Subscripts m and f may be viewed as referring to the matrix and fiber, respectively. Linear elasticity is assumed for the two bodies while the interface is expected to exhibit non-linear behavior, in particular unilateral or frictional contact, or even delamination.

Ω m and Ω f are subjected to body forces f g m and f g f , respectively. Furthermore, surface forces F g m and F g f are applied to boundaries Γ Fm and Γ F f , and displacements u g m and u g f are prescribed on boundaries Γ um and Γ u f , respectively (see Fig. 1 again). In each subdomain, the equilibrium equations, constitutive relations, and boundary conditions have to be verified. Using subscript i to denote a quantity that is valid over region Ω i , the corresponding governing equations read: ∀i ∈ {m, f }, where σ i denotes the Cauchy stress tensor, ε(u i ) the infinitesimal strain tensor, C i the Hooke tensor and n i the outward unit normal, all associated to Ω i . Finally, the interface conditions have to be added to complete the formulation of the boundary value problem. From a general point of view, the latter can be written as:

div(σ i ) + f g i = 0 in Ω i ; σ i = C i ε(u i ) in Ω i ; u i = u g i on Γ ui ; σ i n i = F g i on Γ Fi ; (1)
σ m n m + σ f n f = 0 on Γ ; (2a) g(u m , u f , σ m n m , σ f n f ) = 0 on Γ ; (2b)
where (2a) ensures the equilibrium of the two subdomains along Γ, and (2b) expresses the possible non-linear law between the interface displacements and tractions through a given function g. In the following (see Section 3.3.2), g will be detailed considering the case of contact (unilateral or frictional) and delamination.

Spline discretization and challenges

We recall that the aim of this study is to make use of IGA at the discretization level. IGA being now mature and relatively well-known in the scientific computing community, we do not enter into the details of the technology here, but rather clarify the difficulty of solving multiscale problems similar to (1)-(2) within this framework. For more information, the interested reader is referred to the pioneering contributions [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF] and, e.g., [START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF] for a recent review. Basically, in its common version, the principle of IGA is to use B-spline and NURBS functions to build the approximation spaces when applying Galerkin's method. As stated in the introduction, although this can provide higher accuracy and robustness compared to the standard FEM, this also unfortunately makes the local modeling of geometrically complex inclusions within a spline patch highly challenging. This point actually seems to be closely related to what is called the analysissuitable model issue in the field. Indeed, standard IGA requires a boundary fitted discretization for the analysis while in CAD programs, where the only matter is the rendering of the geometry, entities are described as collections of their boundary surfaces.

As an illustration, let us take the simple example of Fig. 1. In CAD, this object may consist of: (i) a one-patch B-spline surface for the whole plate (linear, 1 element) and (ii) a NURBS trimming curve (quadratic, 2 elements) that forms the boundary between the global and local domains Ω m and Ω f , respectively (see Fig. 2(a)). As a consequence, the underlying spline surface is unaffected by the trimming object and preserves its topology. Conversely, using standard IGA for the analysis of such a problem would require a delicate spline re-parametrization of the whole model to explicitly define the two domains. In practice, this would inevitably lead to the splitting of the geometry into several (tensorproduct) patches with C 0 continuity at the boundaries (see Fig. 2(b) for an example of a boundary fitted NURBS parametrization of the considered problem). In the case of multiple inclusions, the situation is getting even trickier since the resulting geometry of the matrix strongly differs topologically from a square. The new spline model of Fig. 2(b) is commonly referred to as an analysis-suitable model in the sense that it can be easily enhanced using classic spline refinement [START_REF] Lee | Degree elevation of nurbs curves by weighted blossom[END_REF][START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF][START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF] to compute the solution of a corresponding mechanical problem.

This simple example underlines the difficulties of generating analysis-suitable models for multiscale IGA. As mentioned in the introduction, one can perform as in Fig. 2(b); that is, one can strive to remove all trimmed regions by invoking spline re-parametrization strategies, which generally appears cumbersome in practice. Alternatively, one can implement specific numerical schemes adapted to the models with the mentioned defects coming from the geometric modeling [START_REF] Marussig | A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects[END_REF]. This work is heading in this second direction with the development of a novel immersed-like method (see Fig. 2(c)). The challenge then consists in formulating a proper coupling formulation that is adapted to a non-conformal interface (see Fig. 2(c) again). In our context, it is thus requested to handle interfaces that are not only non-conformal but also non-linear, which is the object of the next section.

Isogeometric immersed LaTIn method

As stated in the introduction, we propose to draw inspiration from the ICBM recently introduced in [START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF]. The idea is to separate the difficulty by making appear two different types of interface (similar mindset as in [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF] for instance): non-conformal but perfect interfaces (i.e., that ensure the equilibrium and the standard kinematic compatibility between the subdomains) and non-linear but conformal ones (i.e., where the boundaries of the subdomains are aligned between each other). The method is thoroughly described in this section, starting with the principle, and then focusing on the numerical schemes to address the non-conformal (but perfect) interfaces and finally the non-linear (but conformal) ones.

Principle

More precisely, from interface Γ depicted in Figs. 1 and2, the idea is to extrude it both in the matrix and fiber directions to obtain two conformal layers. The situation is illustrated in Fig. 3. The process results in a possibly non-linear interface Γ that becomes conformal, while the non-conformal aspect is moved inside the subdomains. As for notations, we introduce the superscript j ∈ {l, b} to precise whether we are in the layers (l) or in the bulks (b). Therefore, we denote Ω l m (resp. Ω l f ) the layer within the matrix (resp. fiber) and Ω b m (resp. Ω b f ) the complementary matrix subdomain (resp. fiber subdomain). Finally, the process makes appear two new interfaces, denoted Γ m and Γ f , that allows to couple the layer and the bulk within the matrix and the fiber, respectively, in a non-conformal but perfect manner. Eventually, given the conformal nature of Γ, the non-linear interface behavior can be enforced in a standard manner. In this work, it will be formulated locally, especially pointwise at each interface Gauss point, by means of the LaTIn method [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF]. Regarding the perfect but non-conformal couplings, a well-established Nitsche method will be used following [START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF].

Returning to the formulation, the subdomain-wise Eqs. ( 1) become (see again Fig. 3 for the notations

): ∀i ∈ {m, f } and ∀j ∈ {l, b}, div(σ j i ) + f g i = 0 in Ω j i ; σ j i = C i ε(u j i ) in Ω j i ; u j i = u g i on Γ u j i ; σ j i n j i = F g i on Γ F j i ; (3)
with the following perfect interface conditions: ∀i ∈ {m, f },

σ l i n l i + σ b i n b i = 0 on Γ i ; (4a) u l i -u b i = 0 on Γ i . (4b)
b i (i ∈ {m, f })
for the layer and the bulk, respectively. The possible non-linear coupling between the two layers is performed through interface Γ (still in red) that becomes conformal. The non-conformal aspect is moved inside the subdomains: for i ∈ {m, f }, non-conformal but perfect interfaces Γ i appear (depicted in green) to connect Ω l i and Ω b i . The LaTIn method will be applied for the coupling at interface Γ, while a Nitsche coupling scheme will be used at interfaces Γ i , thereby arriving at a robust, simple, and efficient strategy.

And lastly, Eqs. (2) now concern only the layers, which yields:

σ l m n l m + σ l f n l f = 0 on Γ ; (5a) g(u l m , u l f , σ l m n l m , σ l f n l f ) = 0 on Γ. (5b)
The object in what follows is to derive robust and efficient numerical schemes to solve problem (3)-( 4)-(5a)-(5b).

Non-conformal coupling with a Nitsche approach

Let us focus here on the enforcement of Eqs. ( 3)-( 4), i.e. the formulation and discretization related to the perfect and non-conformal couplings across Γ i , ∀i ∈ {m, f }. As stated above, we consider a Nitsche method for this purpose. Given its primal nature, variational consistency, positive definiteness, symmetry, and good conditioning properties, the Nitsche approach seems to have established itself over the years in the immersed boundary community for the perfect and non-conformal coupling of linear subdomains (see [START_REF] Hansbo | Nitsche's method for interface problems in computa-tional mechanics[END_REF][START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF][START_REF] Claus | A stable and optimally convergent latin-cutfem algorithm for multiple unilateral contact problems[END_REF][START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF] to name a few). More precisely, we follow here reference [START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF] that relies on the minimal stabilization scheme [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF].

Nitsche coupling formulation

Let us start by defining the functional spaces U j i and V j i that contain the displacement solution and test functions, respectively: ∀i ∈ {m, f } and ∀j ∈ {l, b},

U j i = {u j i ∈ [H 1 (Ω j i )] d , u j i | Γ u j i = u g i } ; V j i = {v j i ∈ [H 1 (Ω j i )] d , v j i | Γ u j i = 0}.
To account for the coupling in the Nitsche approach, both the interface conditions (4a) and (4b) are imposed weakly and a penalty-like stabilization term is added to ensure the ellipticity of the boundary value problem. The corresponding weak formulation can be written as: ∀i ∈ {m, f }, find

(u l i , u b i ) ∈ U l i × U b i such that, j∈{l,b} a j i (u j i , v j i ) - Γi {C i ε(v i )} n l i • [[u i ]]dΓ - Γi [[v i ]] • {σ i } n l i dΓ +ζ i Γi [[v i ]] • [[u i ]]dΓ = j∈{l,b} l j i (v j i ), ∀(v l i , v b i ) ∈ V l i × V b i ; (6)
where standard bilinear form a j i and linear form l j i associated to subdomain Ω j i read:

a j i u j i , v j i = Ω j i ε(v j i ) : C i ε(u j i )dΩ ; l j i (v j i ) = Ω j i v j i • f g i dΩ + Γ F j i v j i • F g i dΓ.
In Eq. ( 6), [[u i ]] = u l i -u b i and {σ i } n l i denote the displacement jump and the stress flux across Γ i , respectively, and ζ i is a stabilization parameter that shall depend on the mesh discretizations and material properties of the two subdomains to be coupled. Following [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF][START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF], we adopt the one-sided flux, i.e. that we consider only the term for the conformal subdomain:

{σ i }n l i = σ l i n l i = C i ε(u l i )n l i ;
and we take ζ i as:

ζ i = β i ((h l i ) -1 + (h b i ) -1 ) with β i = 6(p max i ) 2 × 8(E max i /(1 -2ν max i )) ;
where h j i is the maximum element size in Ω j i , and p max i , E max i and ν max i are the maximum polynomial degree, Young modulus, and Poisson ratio, respectively, of the subdomains (Ω j i , ∀j ∈ {l, b}) to be coupled. With such choices, the formulation is rather simple and has proved to be stable and accurate (see again [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF][START_REF] Wei | Immersed boundaryconformal isogeometric method for linear elliptic problems[END_REF]). Here, since the layer is of the same material as that of the corresponding bulk, we obviously have

E max i = E l i = E b i = E i and ν max i = ν l i = ν b i = ν i grouped in the Hooke tensor C i , ∀i ∈ {m, f }.
Finally, let us condense formulation (6) by writing it as follows: ∀i ∈ {m, f },

find u i ∈ U i such that, a i (u i , v i ) = l i (v i ), ∀v i ∈ V i , (7)
where

U i = U l i × U b i , V i = V l i × V b i
, and a i (u i , v i ) and l i (v i ) are the left-hand side and right-hand side of (6), respectively. We recall that Eq. ( 7) enforces ( 3)-( 4). It thus remains to ensure [START_REF] Blanzé | Contact problems in the design of a superconducting quadrupole prototype[END_REF], which is the object of the next Section 3.3.

Implementation care with immersed methods

The implementation of formulation [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF] requires special care to (i) evaluate integrals over pieces of d-variate elements (for subdomains Ω b i ), (ii) evaluate integrals over (d -1)-variate interfaces that cut the meshes (for computing the Nitsche coupling operators over Γ i ), and (iii) manage the possible ill-conditioning of the stiffness matrices over Ω b i (i.e., the operators associated to bilinear forms a b i ). For point (i), for simplicity and robustness, we use the recursive, quad-tree based quadrature approach employed in FCM [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Schillinger | The finite cell method: A review in the context of higher-order structural analysis of CAD and image-based geometric models[END_REF]. For sure, more geometrically faithful quadrature rules [START_REF] Nagy | On the numerical integration of trimmed isogeometric elements[END_REF][START_REF] Kudela | Efficient and accurate numerical quadrature for immersed boundary method[END_REF][START_REF] Antolin | Isogeometric analysis on V-reps: first results[END_REF][START_REF] Legrain | Non-negative moment fitting quadrature rules for fictitious domain methods[END_REF][START_REF] Garhuom | Non-negative moment fitting quadrature for cut finite elements and cells undergoing large deformations[END_REF] could also have been considered since our approach is generic in terms of integration schemes. Regarding point (ii), we pay attention to discretize interfaces Γ i as the intersection of the two meshes on both sides of it (see Fig. 4). This ensures that on each interface element, all the involved basis functions from the two subdomains to couple are polynomials rather than piecewise polynomials. Finally, for point (iii), we start by removing from the corresponding stiffness matrices all the degrees of freedom associated with basis functions whose support does not intersect with the integration domain. Then, Figure 4: Interface discretization for coupling non-conformal meshes using the Nitsche approach: to compute integrals over Γ m , the interface discretization is chosen as the intersection of two nonconformal meshes on both sides of it. Interface Γ m is shown, Γ f is of course discretized in a similar manner.

we make use of a diagonal scaling preconditioner as in [START_REF] Antolin | Isogeometric analysis on V-reps: first results[END_REF]. For a given (n × n) matrix K collecting terms K i,j for the ith line and jth column, the preconditioner is simply defined as:

D = diag 1/ K 1,1 , 1/ K 2,2 , . . . , 1/ K n,n ,
and is applied on the left and on the right of the linear system of interest. With the latter, we never encountered instabilities in our numerical experiments.

Enforcement of non-linear interface behaviors with the LaTIn method

Let us now move to the non-linear interface Γ. The object here is to add Eq. ( 5) to formulation [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF]. In order to do so, we recall that we consider the LaTIn method [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation[END_REF] which has the interest of leading to a naturally parallel algorithm in case of multiple local models. More precisely, the idea of the LaTIn approach is to separate the equations in two groups: in our case, the first one concerns the subdomains equations, which are all linear and subdomain-wise independent; and the second one gathers the interface equations, the latter being non-linear and local (i.e. pointwise independent). We then iterate between these two sets of equations by means of two search directions to obtain the solution of the problem; that is, the one that satisfies all the equations. In the following, we first outline the LaTIn numerical scheme from a general domain decomposition point of view and then briefly specify the method depending on the interface non-linear behavior (unilateral or frictional contact, and delamination). For further details, the interested reader is advised to consult, e.g., [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF][START_REF] Claus | A stable and optimally convergent latin-cutfem algorithm for multiple unilateral contact problems[END_REF] and the other works cited hereafter.

General LaTIn numerical scheme

Separation of the equations To formulate the method, let us start by introducing the following interface fields for subdomains Ω i : ∀i ∈ {m, f }, λ i represents a surface force applied to Ω i through boundary Γ, and w i is the trace of u l i (or directly of u i ) over Γ. With these new variables, we can define the first partial solution spaces A i that group the solutions of the linear equations set on the subdomains: ∀i ∈ {m, f },

A i : (λ i , w i ) such that        (7) σ l i n l i = λ i on Γ w i = u l i on Γ . ( 8
)
The spaces A i are affine spaces often called spaces of admissible fields. We denote their union:

A = i∈{m,f } A i .
Then, let us define additional force and displacement interface fields ( λi , ŵi ) i∈{m,f } that belong to the second partial solution space L, i.e. that verify the behavior of the interface:

L : ( λi , ŵi ) i∈{m,f } such that λm + λf = 0 on Γ g( ŵm , ŵf , λm , λf ) = 0 on Γ . (9)
L is a manifold referred to as the local space since it is defined by pointwise independent equations. Finally, it is required to add the search directions k + and k -to communicate between spaces A and L and thus close the problem. This is performed in a mixed way as follows: ∀i ∈ {m, f },

k + : λi -λ i -k + i ( ŵi -w i ) = 0 on Γ ; (10a) k -: λi -λ i + k - i ( ŵi -w i ) = 0 on Γ ; (10b)
where k + i and k - i represent interface stiffnesses. In practice, we choose k + i = k - i = k i > 0 which is the classical setting of the LaTIn algorithm [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation[END_REF]. Note at this stage that problem ( 8)-( 9)-( 10) does correspond to our problem of interest ( 7)-( 5) since Eqs (10a) and (10b) obviously lead to λi = λ i and ŵi = w i , ∀i ∈ {m, f }.

Iterative algorithm With the above separation of equations in hand, we can then perform a fixed point numerical scheme to solve the problem. More precisely, for the nth iteration, starting with initial guesses

S (0) i = (λ (0) i , w (0) i ) ∈ A i , ∀i ∈ {m, f }, we subsequently perform two steps: 1. Local (non-linear) stage: ∀i ∈ {m, f } , given S (n-1) i = (λ (n-1) i , w (n-1) i ) ∈ A i , find Ŝ(n-1) i = ( λ(n-1) i , ŵ(n-1) i ) ∈ L ∩ k + . (11)
2. Global (linear) stage:

∀i ∈ {m, f } , given Ŝ(n-1) i = ( λ(n-1) i , ŵ(n-1) i ) ∈ L, find S (n) i = (λ (n) i , w (n) i ) ∈ A i ∩ k -. (12)
To ensure the convergence of the algorithm, a relaxation step is also performed at the end of the linear stage:

∀i ∈ {m, f } , S (n) i ← θS (n) i + (1 -θ)S (n-1) i ,
where we choose θ = 0.5 following again [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation[END_REF]. The overall procedure is illustrated in Fig. 5. In addition, the search direction k i should represent, as close as possible, the condensed stiffness on Γ of the complementary part of the whole domain to minimize the number of iterations of the algorithm (see, e.g., [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF])). In 2D, we therefore take

k m = E f L f and k f = Em Lm
, where E f (resp. E m ) and L f (resp. L m ) are the Young modulus and characteristic length of subdomain Ω f (resp. Ω m ). Finally, an indicator of error is used to quantify the distance between A and L for two successive partial solutions. It is written in an energy norm such that:

η = i∈{m,f } Γ k i w (n) i - ŵ(n-1) i 2 + 1 k i λ (n) i - λ(n-1) i 2 dΓ i∈{m,f } Γ k i w (n) i 2 + ŵ(n-1) i 2 + 1 k i λ (n) i 2 + λ(n-1) i 2 dΓ .
Figure 5: Latin method: a schematic depiction of the iterative process ( S(i) is an approximation before relaxation).

Implementation Benefiting from the above LaTIn algorithm, it becomes clear that (i) the nonlinear stage [START_REF] Bouclier | Development of a new, more regular, mortar method for the coupling of nurbs subdomains within a nurbs patch: Application to a non-intrusive local enrichment of nurbs patches[END_REF] can be performed locally and (ii) the linear stage ( 12) can be performed in parallel between the different subdomains. This is the interest of the method when considering multiple local models connected through non-linear interfaces. We somehow naturally end up with a non-linear parallel domain decomposition algorithm.

From a practical point of view, the local equations are solved at each integration point of the conformal interface Γ. Section 3.3.2 will provide some details depending on the interface non-linear behavior considered. Note that the computations can be performed in parallel between the different integration points. The hat quantities ( λi , ŵi ) i∈{m,f } are therefore "heart" quantities in the sense that they are only defined at the interface Gauss points and they are not directly discretized with some basis functions. Obviously, to make it possible, the other interface fields (λ i , w i ) i∈{m,f } are evaluated at the corresponding integration points prior to performing the local stage.

Regarding the linear stage, the problems to be solved are actually subdomain-wise elastic problems subjected to generalized Robin boundary conditions. Replacing λ i in (8) by its expression from (10b), and returning to u l i or even u i by applying the trace over Γ, the linear problems reads at iteration n: ∀i ∈ {m, f }, find u

(n) i ∈ U i such that, a i (u (n) i , v i ) + Γ k i u (n) i • v i dΓ = l i (v i ) + Γ ( λ(n-1) i + k i ŵ(n-1) i ) • v i dΓ, ∀v i ∈ V i , (13)
where we reuse the notations of Section 3.2 for a i , l i , U i and V i . In the discrete setting, since we consider here a conformal interface Γ, problems (13) should be rather simple to solve and do not require additional stabilization terms as in the case of non-conformal LaTIn interfaces (see, e.g., [START_REF] Claus | A stable and optimally convergent latin-cutfem algorithm for multiple unilateral contact problems[END_REF]).

The stabilization is actually moved to the Nitsche coupling in Section 3.2 in our approach. Here, we simply take the existing (d -1)-variate (boundary) basis functions associated to subdomain Ω i that generate Γ to discretize w i . We, therefore, end up with an interface mass matrix for operator

Γ u i • v i dΓ.
Then, for the discretization of right-and sides Γ λi • v i dΓ and Γ ŵi • v i dΓ, we actually only need the evaluation of λi and ŵi at the interface Gauss points, which is the case at the end of the local stage. Note that this corresponds to performing a L 2 -projection of the heart quantities ( λi , ŵi ) i∈{m,f } onto the continuous space used for w i . Finally, once we have u i in hand, we can obtain w i by using the control variables of u i associated to the boundary control points generating Γ, and subsequently evaluate w i at the interface Gauss points, and λ i at the same locations using the descent search direction (10b).

Remark 1. It may be noticed at this stage that interface Γ in the proposed LaTIn strategy actually needs to be geometrically conformal but it can be non-matching. In other words, interface Γ needs to be aligned with the boundary edges of the subdomains to be coupled but the discretizations of these subdomains over Γ (mesh refinements and/or polynomial degrees) can differ (see, e.g., [START_REF] Bouclier | Development of a new, more regular, mortar method for the coupling of nurbs subdomains within a nurbs patch: Application to a non-intrusive local enrichment of nurbs patches[END_REF] for the nomenclature). With a non-matching interface, only additional care is required to define the integration points over Γ for the local equations: in this case, interface Γ is discretized as the intersection of the two meshes on both sides of it and:

n gp = max i∈{m,f } (p i ) + 1
Gauss points are taken per (d -1)-variate elements of the latter. That being said, let us underline that taking different discretizations for the two layers may appear strange since there is no reason to prioritize the accuracy of the interface fields of one of the subdomains to correctly represent the non-linear interface behaviors of interest.

Treatment of the non-linearity in specific cases

Now that the general LaTIn algorithm has been presented, let us give a few more details regarding the local stage [START_REF] Bouclier | Development of a new, more regular, mortar method for the coupling of nurbs subdomains within a nurbs patch: Application to a non-intrusive local enrichment of nurbs patches[END_REF]. For simplicity, we omit in what follows the superscript (n). ∀i ∈ {m, f }, starting with quantities (λ i , w i ), the aim is to find the hat quantities ( λi , ŵi ), at each interface Gauss point, such that Eqs. ( 9) and (10a) are satisfied. Once again, we consider contact (unilateral or frictional) and delamination (using a cohesive interface) for the non-linear behaviors. The object of this part is, therefore, to specify Eq. ( 9) (i.e., function g) in these cases, and to provide some insights regarding the resolution. In order to do so, we actually need to introduce a pseudo-time to account for the load history and then write the LaTIn method in a quasi-static framework. Consequently, the above LaTIn algorithm is slightly modified as the search directions now link, at each time step t, the surface forces λ

(t) i , λ (t)
i and the pseudo-velocities ẇ(t) i , ŵ

i such that:

k + : λ(t) i -λ (t) i -k v i ( ŵ(t) i - ẇ(t) i ) = 0 on Γ ; (14a) k -: λ(t) i -λ (t) i + k v i ( ŵ(t) i - ẇ(t) i ) = 0 on Γ. (14b)
Next, a standard implicit scheme is used to express the velocities as a function of the displacements:

ẇ(t) i = w (t) i -w (t-1) i ∆t and ŵ(t) i = ŵ(t) i - ŵ(t-1) i ∆t . (15)
Finally, contact and cohesive behaviors require the use of a local basis, at each interface Gauss point, defined by the tangent and normal vectors associated with the interface. In the following, the subscript τ (resp. n) refers to the tangential (resp. normal) component of a quantity expressed in this basis. The normal vector is arbitrarily defined from the matrix to the fiber.

Remark 2. Note that k v i in Eqs. (14a) and (14b) does not strictly have the same dimension as k i in Eqs. (10a) and (10b), due to the introduction of the time increment ∆t. However, let us emphasize that we can take in practice ∆t = 1 since we are in the case of a quasi-static transformation ( i.e., without real dynamic effects). Therefore, we take for our numerical experiments the same values for k v i as for k i , ∀i ∈ {m, f }.

Unilateral and frictional contact First, let us present the equations for a contact interface that follows the Coulomb friction law. The case of unilateral contact (i.e., contact without friction) will then be easily obtained by simplifying these equations, see remark 3. Denoting by j n and µ the possible initial normal gap between the two subdomains and the friction coefficient, respectively, the interface quantities need to satisfy in the case of friction contact the following equations at each interface Gauss point and time step:

λ(t) m + λ(t) f = 0, (Mechanical equilibrium) ŵ(t) mn - ŵ(t) fn + j n ≥ 0, (Non interpenetration) λ(t) m n - λ(t) fn ≥ 0, (Positive reaction forces) ( λ(t) mn - λ(t) f n )( ŵ(t) mn - ŵ(t) f n + j n ) = 0, (Signorini conditions)    If ∥ λ(t) iτ ∥ < µ| λ(t) in | then ŵ(t) mτ - ŵ(t) f τ = 0, If ∥ λ(t) iτ ∥ = µ| λ(t) in | then ∃α > 0, ŵ(t) mτ - ŵ(t) fτ = -α λ(t) mτ - λ(t) fτ . (Sticking)
(Sliding) [START_REF] Claus | A stable and optimally convergent latin-cutfem algorithm for multiple unilateral contact problems[END_REF] Remark 3. In case of unilateral contact, the last two equations are replaced by the following one:

λ(t) mτ - λ(t) fτ = 0,
and, as the equilibrium of the surface forces is always verified, this leads to null tangent surface forces:

λ(t) mτ = λ(t) fτ = 0.
As a consequence, the problem of the local stage consists here in solving Eqs. ( 16), (14a) and ( 15). This can actually be easily performed by defining two indicators C (t) n and G (t) τ , which are explicitly determined from known fields at time step t:

C (t) n = ŵ(t-1) mn - ŵ(t-1) fn + j n + ∆t ẇ(t) mn -ẇ(t) mn -∆t k v -1 m λ (t) mn -k v -1 f λ (t)
fn , and:

G (t) τ = λ (t) mτ -k v m ẇ(t) m - k v f k v m + k v f λ (t) mτ + λ (t) fτ -k v m ẇ(t) mτ -k v f ẇ(t) fτ . C (t)
n refers to the contact state of the corresponding interface Gauss point at time step t: it can be either in contact (C

(t) n ⩽ 0) or not (C (t) n > 0). In case of contact, the second indicator G (t)
τ indicates whether the solution is sticking or sliding. Then, the solution is explicit and can be determined without any difficulty. For more details regarding the solution algorithm, we urge the interested reader to consult, e.g., [START_REF] Blanzé | Contact problems in the design of a superconducting quadrupole prototype[END_REF][START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF].

Remark 4. Two frameworks can be used to take into account the quasi-static assumption. The first one is the incremental approach, where the LaTIn method can be viewed as a non-linear solver used to find the solution for one time step before determining the next time step. The second framework is the non-incremental approach. The LaTIn method is seen as a more sophisticated solver since all time steps are embedded in the LaTIn iterations. Indeed, in this variant, each local stage and linear stage consists in finding an approximation of the solution for all the time steps. This non-incremental approach is particularly adapted to multi-parametric studies [START_REF] Roulet | A parallel strategy for the multiparametric analysis of structures with large contact and friction surfaces[END_REF] or time/space decomposition with model order reduction [START_REF] Giacoma | Toward an optimal a priori reduced basis strategy for frictional contact problems with latin solver[END_REF], which is not our case here so we carry out the incremental approach.

Cohesive model Let us now specify the equations regarding the cohesive model for the interface and how it is managed within the LaTIn method. Cohesive models are based on the idea that the interface is defined by a stiffness that links the surface force to the normal displacement jump. Initially, the cohesive interface is undamaged, but during loading, the interface can be progressively damaged until it breaks completely. A fully damaged interface can be considered as a contact interface. In this part, we assume that the interface is opening and the normal displacement jump is positive. The case of a closed interface and a null displacement jump is considered as contact behavior in the normal direction and is enforced as we described previously. Therefore, to avoid cumbersome notations, we do not specify the positive part of the displacement jump in the following.

In general terms, the constitutive law can be written as follows:

λ(t) fn = k [[ ŵ(t) n ]] [[ ŵ(t) n ]],
where k is the stiffness of the interface which depends on the normal displacement jump [[ ŵ(t) n ]] = ŵ(t) mn -ŵ(t) fn . There are several laws that describe cohesive models. In our case, we consider a simple bilinear law [START_REF] Alfano | Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues[END_REF], as illustrated in Fig. 6. It begins with linear elastic behavior which initial stiffness is denoted as k 0 . The displacement jumps δ 0 and δ c , corresponding to the initiation and completion of damage are derived from the energy restitution rate G c and the critical stress σ c as δ c = 2 Gc σc and δ 0 = 0.1 • 2 Gc σc , respectively. As the displacement jump increases, the critical stress σ c is reached and the interface starts to be damaged, the surface force decreases and in the case of unloading, the elastic return follows a damaged stiffness k. In this case, the complete problem to be solved at the local stage consists in finding a solution that verifies the above constitutive law, the equilibrium of the interface, and the search direction equations, which read:

λ(t) fn = k [[ ŵ(t) n ]] [[ ŵ(t) n ]], (Constitutive law) λ(t) fn + λ(t) mn = 0, (Equilibrium of surface forces) λ(t) m -λ (t) m -k v m ( ŵ(t) m -ẇ(t) m ) = 0, (Search direction for the matrix) λ(t) f -λ (t) f -k v f ( ŵ(t) f - ẇ(t) f ) = 0. (Search direction for the fiber) (17)
At this step, it is often convenient to change the search direction and consider an infinite one (i.e., choosing k v m → ∞ and k v f → ∞), which leads to simply replace the last two equations of ( 17) by ŵ(t)

m = w (t) m and ŵ(t) f = w (t)
f . As a result, the computation of the surface forces are fully explicit and there is no need for a non-linear solver. We advise the interested reader to consult [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Violeau | Micromodel-based simulations for laminated composites[END_REF][START_REF] Fernandez | The effect of the mesh refinement on a multiscale domain decomposition method for the non-linear simulation of composite structures[END_REF] for more details on this topic.

Remark 5. So far we have written the problem for the normal component of the interface quantities. We can define the same kind of equation for the tangential ones. The main difference is that compressive behavior does not make much sense for the tangential component, shear damage will occur in any direction in the tangential plane. Therefore, the tangential displacement jump is considered as an absolute value.

Linear stage under a quasi-static assumption Finally, let us underline that the linear stage ( 13) is slightly modified to take into account the quasi-static assumption. To do so, we simply need to repeat the procedure to obtain [START_REF] Burman | Cutfem: discretizing geometry and partial differential equations[END_REF], but with ( 14) as search directions. Combining this with the time integration scheme [START_REF] Casquero | The divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics[END_REF], the problem at the linear stage now reads: ∀i ∈ {m, f }, at time step t find u

(t) i ∈ U i such that, a i (u (t) i , v i ) + Γ 1 ∆t k i u (t) i • v i dΓ = l (t) i (v i ) + Γ λ(t) i + k i ŵ(t) i + 1 ∆t k i w (t-1) i • v i dΓ, ∀v i ∈ V i .

Numerical results

To assess the performance of the developed isogeometric IBCM LaTIn scheme, we now present a series of numerical experiments in 2D that covers different geometries, discretizations and non-linear interface scenarios. More precisely, we investigate 4 test cases: (i) one inclusion in unilateral contact with the matrix, (ii) several inclusions in unilateral contact with the matrix, (ii) a two-body assembly with frictional contact inside, and (iv) one inclusion in delamination with the matrix.

Figure 7: Matrix with a central inclusion in frictionless contact: description and data of the problem.

Figure 8: Model decomposition: each part (matrix and inclusion) is composed of two subparts, a conformal layer, and a background mesh. The meshes used for the solution approximation in the matrix are colored in black while the meshes used for the solution in the inclusion are colored in red. grey color corresponds to the integration area for stiffness operators.

Matrix with a central inclusion in frictionless contact

To start with, we considered a single circular inclusion in a rectangular matrix, as described in Figure 7. Both matrix and inclusion materials were linear elastic, and we assumed plane stress state and small deformations. However, a large elastic contrast between the matrix and inclusion materials was prescribed: in particular, the Young modulus was chosen 80 times larger for the inclusion than for the matrix. Contact behavior without friction is considered on the matrix/inclusion interface. Due to the symmetry of the problem, only one quarter of the area was actually modeled with corresponding symmetrical boundary conditions (see again Figure 7). We decomposed the model into two independent subdomains corresponding to the matrix and to the inclusion. Then, following our approach, in each model we added two conformal layers on the matrix/inclusion interface (see Figure 8). The thickness of each of the conformal layers was 0.2mm. Thus, each subdomain consists of two subparts: a conformal layer and a regular background mesh. Each layer was coupled by Nitsche method with its corresponding model (through a non-conformal but perfect interface). The stiffness operators of these corresponding regular background meshes, following FCM, were built by integrating solely in the physical domain which is colored in grey on Figure 8. Together, the subdomains were connected by the developed LaTIn method through interface Γ (a nonlinear but conformal interface). We recall that this technique allows us to circumvent the difficulty of dealing with a coupling interface both non-conformal and incorporating contact. The following results were obtained for a rather fine mesh containing 32 × 32 NURBS quadratic elements in each layer and 80 × 80 B-spline quadratic elements in the matrix background mesh. We chose the same element size for the inclusion's background mesh. The search directions for the LaTIn method, due to the significant contrast between the properties of the materials, were as the Young modulus of the complement divided by its characteristic size (see Section 3.3). The LaTIn iterations were stopped once indicator η reaches 10 -5 . All the shown results correspond to converged LaTIn solutions. Figure 9 shows the obtained displacement on a deformed geometry with a scale factor of 10. Figure 10 depicts stresses on the deformed geometry. We observe that since σ yy is continuous on the top of the inclusion, materials are in contact in this area, while there is no contact on the right of the inclusion. In addition, Figure 11 (left) shows the status of the integration points which is defined by the sign of the contact indicator C (t) n . Normal tractions on the interface Γ are depicted on Figure 11 (right) where angle α is the angle between the position vector of a point and the positive direction of x axis. The absence of oscillations in these results numerically accounts for the stability of the numerical scheme.

Then, the convergence of our method with the mesh refinement was studied. For this purpose, four meshes were considered, see Figure 12. In addition, NURBS of degrees 2 and 3 were used. Finer meshes were obtained from the initial one in a hierarchical manner by consequently doubling the number of elements in each direction. As a reference solution, we took a very refined one following our approach (320 × 320 cubic elements in background meshes and 128 × 128 cubic elements in the annuli). We investigated the convergence rate in the energy norm, which is defined as:

Err = ∥u h -u ref ∥ ∥u ref ∥ = i∈{m,f } Ωi ε(u h i -u ref i ) : C i ε(u h i -u ref i ) dΩ i∈{m,f } Ωi ε(u ref i ) : C i ε(u ref i ) dΩ
, where u h and u ref denote the solution of interest and the reference solution, respectively. For the evaluation of the error, the coarser solutions were interpolated at the integration points of the reference solution. We show the convergence plots for 2-and 3-degree splines on Figure 13a with respect to size indicator h. It can be observed that when evaluating the error over the whole domain the rate of convergence for the quadratic solution is almost the one encountered in linear elasticity (h 2 ). This demonstrates the superior accuracy of our method with respect to the state-of-the-art of the LaTIn technique that usually restricts to linear FEM. Nevertheless, it can also be noticed that going to cubic splines does not improve much the convergence. This is expected due to the lack of regularity of the contact problem [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]. For completeness on this point, we finally computed the error only in the matrix background mesh. The associated convergence curves are given in Figure 13b. This time, we did obtain the convergence rates usually reached in linear elasticity for both the quadratic and cubic solutions. The solution in the matrix background mesh being quite distant from the contact interface is actually less affected by the lack of regularity of the problem.

Remark 6. This mesh refinement study suggests that it might be worthwhile to discretize the layers with splines of relatively low polynomial degrees (say 2) and to use splines of higher polynomial degrees (say 3 or more) for the bulk. This would allow to correctly catch the non-linear behavior at the interface while increasing the per-degree-of-freedom accuracy in the bulk. Note that this treatment is possible with our method without any effort, since it relies on a stabilized Nitsche scheme that allows to couple spline immersed subdomains of different degrees.

Matrix with multiple inclusions

We are now interested in the convergence of the LaTIn method for different number of inclusions. To this end, we considered models consisting of an 8mm × 8mm square-shaped matrix with 1, 2, 4, or 16 (uniformly distributed) inclusions, see Figure 14. Each inclusion was of a 1mm diameter, each annulus thickness was 0.1mm. We apply symmetry boundary conditions on the left and bottom edges of the matrix and a horizontal displacement on the right edge. As in the previous case, both matrix and inclusions materials were linear elastic and we assumed again plane stress state and small deformations. We considered the same material for all inclusions. Material parameters were the same as in 4.1, so that there was a significant contrast in the matrix and inclusions Young moduli. Furthermore, we considered once again frictionless contact behavior on all matrix/inclusion interfaces.

For the numerical simulations, we took rather fine meshes. Matrix background mesh contained 32 × 32 elements, each of its conformal layers being composed of 16 × 4 elements. All inclusions had identical meshes but, contrary to the previous test case, the inclusions' background meshes were eight times finer than the matrix one (see again Figure 14 for meshes representation).

Figure 15 shows stress distribution for the 4-inclusion problem on a deformed geometry with a scale factor of 10. As can be seen from the continuity of stresses σ yy the highest and the lowest parts of inclusions are in contact with the matrix while some detachments appear on the left and right parts of the inclusions. This seems consistent with the Poisson effect of the plate. For completeness on this point, Figure 16 depicts the contact status of interface integration points for this case. Eventually, the convergence of the LaTIn method for 1, 2, 4, and 16 inclusions problems was studied. Figure 17 shows the corresponding convergence curves. It can be observed that the convergence of the developed algorithm is almost independent of the number of inclusions. More precisely, the algorithm needs only about 30 iterations for different numbers of inclusions to attain an indicator η of 10 -5 . Recalling that the different interface and subdomain computations during the local and global stages, respectively, can be carried out in parallel, our approach can thus be used as an efficient non-linear domain decomposition solver. The matrix actually plays the role of the coarse problem of parallel domain decomposition approaches [START_REF] Hirschler | A dual domain decomposition algorithm for the analysis of non-conforming isogeometric kirchhoff-love shells[END_REF][START_REF] Micha L Bosy | A domain decomposition method for isogeometric multi-patch problems with inexact local solvers[END_REF][START_REF] Olof B Widlund | Bddc deluxe algorithms for twodimensional h (curl) isogeometric analysis[END_REF]. It allows to transmit directly the information all over the inclusions, which makes the algorithm scalable.

Frictional contact problem

In the next step, we applied our method to a model with frictional contact at interfaces. To this end, a model with two curvilinear interfaces was created (see Figure 18). This model was inspired by [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF] where interfaces were straight.

The model comprised two domains, Ω 1 and Ω 2 , where Ω 1 , analogous to the matrix model in the previous examples, was itself composed of 2 parts, physical and fictitious ones. Two conformal layers were built on curvilinear non-linear interfaces 1 and 2 to treat the non-conformal aspect of domain Ω 1 . Domain Ω 2 , in its turn, was a simple solely physical domain that had a conformal interface discretization with interfaces 1 and 2. Frictional contact between Ω 1 and Ω 2 was modeled by the Coulomb friction law with friction coefficients µ 1 and µ 2 on interfaces 1 and 2, respectively. In addition, interface 3 imposed boundary conditions on the right boundary of Ω 2 . It was a frictionless contact interface with initial gap j between the interface and domain Ω 2 . During the linear stage, the displacement w corresponding to interface 3 was forced to be zero as if there were a rigid body. Then, interface forces λ were calculated from w with the use of the search direction. Here, both materials are linear elastic with equal properties. The discretization used to solve the problem involved 16 × 48 quadratic B-spline elements for Ω 1 , 32 × 2 quadratic B-spline elements for the conformal layers, and 32 × 32 quadratic B-spline elements for Ω 2 (see again Figure 18).

The problem studied in this section is quasi-static, and it involved two load steps. In the first step, at time t 0 , a vertical load F 1 was applied on the top boundary of domain Ω 1 , causing interfaces 1 and 2 to come into contact and thus resulting in non-zero contact forces on the whole of these interfaces. During this step, there was no external load on Ω 2 , so F 2 = 0. In step two, at t = t 1 , load F 2 = F max 2 was applied on domain Ω 2 while F 1 is being kept constant. The aim was to investigate the effect of 20 1.

In this situation of frictional contact, we considered that the algorithm had converged when the criterion η < 10 -7 was achieved. The required number of iterations varies depending on the time step and on the values of µ 1 and µ 2 . However, it is of the order of a hundred iterations.

The contact zone at interface 3 can vary based on the chosen friction coefficients. It could be fully in contact, partially in contact, or not in contact. We computed the resulting reaction force R 3 on this interface to quantify the contact status. Results are depicted in Figure 19, where we can see 2 main zones. A zone for which the friction coefficients are too high and prevent contact. It results in a null reaction force (see zone close to point (d) in Figure 19). In the second zone, a partial or full contact is reached on interface 3: the size of the contact area depends on the friction coefficient, the smaller the friction coefficients, the higher the reaction force (see zone close to point (a) in Figure 19). Remark 7. It can be seen that contrary to the model with straight interfaces from [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF], the reaction force for µ 1 = µ 2 = 0 is bigger than the applied horizontal load F 2 . Indeed, when an interface is in contact, it results in nonzero normal forces. For the curvilinear interfaces, these normal forces have a nonzero horizontal component, which gives additional horizontal load on domain Ω 2 .

To give a closer look at the results, Figure 20 shows horizontal displacement on the deformed configuration for four particular pairs of (µ 1 , µ 2 ). It can be seen that interface 3 becomes completely in contact if there is no friction on interfaces 1 and 2 (Figure 20a). Then, Figure 20b and Figure 20c depict two cases where interface 3 is only partially in contact. Finally, for µ 1 = µ 2 = 0.6 (Figure 20d) there is no contact on interface 3.

Matrix with one inclusion and delamination

Finally, we investigated the case of a cohesive inclusion/matrix interface with the last example. The problem description is provided in Fig. 21(a), and we took the same material properties as in the first example (see Section 4.1) for the matrix and inclusion. Regarding discretization, the matrix contained 32 × 16 NURBS elements in the conformal layer and 32 × 32 B-spline elements in the background mesh. For the inclusion, we opted for a mesh two times finer than for the matrix, both in the bulk and in the layer. Additionally, we employed different splines degrees: p = 2 for the inclusion and p = 4 for the matrix (see Figure 21(b)). Therefore, in view of assessing the robustness of the developed approach (and especially regarding the LaTIn strategy), we placed ourselves in the case of a geometrically conformal but non-matching interface here (see remark 1).

Dirichlet boundary conditions were imposed on the top and bottom of the region. The simulation involved 40 time steps with linear increments, resulting in a final displacement of u g = 0.2mm. To simulate the delamination of the matrix/fiber interface, we employed the standard bi-linear law described in sec. 3.3.2. In our computations, we set G c = 0.5 mm•MPa, and σ c = 9 MPa. The reaction force at the bottom of the domain is shown in Figure 22 at each time step. More precisely, Figure 22 shows three configurations : (a) the case of an undamaged structure, for which the behavior of the interface remains in the first linear part (see again Fig. 6); (b) the structure with the cohesive interface starting from an undamaged one and (c) the structure with a fully damaged interface. As expected, the behavior of the standard structure falls between the two extreme configurations. It starts with an undamaged interface to reach the behavior of a fully damaged interface as the load increases. In fact, damage to the interface progressively reduces the overall stiffness of the structure.

Finally, Figure 23 presents the damage parameter at the integration points of the interface at the last time step. A parameter value of zero indicates an undamaged interface, while a value of one corresponds to complete interface damage. As it can be seen, the interface is fully damaged on the top and at the bottom, where jumps in normal displacement are maximum, while it stays intact on the left and on the right, where these jumps are insignificant. This can also be seen in Figure 24, as the stress at the top and bottom of the inclusion is zero. Since the interface is mostly broken, the inclusion is hardly loaded in the y-direction.

Conclusions

In this work, we developed an efficient and accurate solver for the simulation of composite microstructures composed of multiple inclusions connected to the matrix with unilateral contact, frictional contact, and delamination. To do so, three interconnected ingredients were employed: (i) IGA, which has established itself as an HPC tool nowadays, was used for the geometric representation of the different interfaces and for the problem discretization, (ii) the pragmatic immersed boundary-conformal method was considered to offer great geometric flexibility while retrieving conformal interfaces, and (iii) the LaTIn technique was implemented in this context to obtain a parallel domain decomposition algorithm capable of treating any non-linear interface behavior. In other words, the key point was to separate the difficulty of treating interfaces both non-conformal and non-linear, so as to obtain a simple numerical scheme which stability is ensured within the bulk equations through Nitsche coupling. The strategy was applied to a range of numerical examples in 2D and proved to be accurate regardless of the geometry, discretization, and non-linear interface scenarios. Although the focus in terms of applications was on composite microstructures, we underline that the proposed algorithm may also find success in other contexts, such as for large assemblies of interacting solids [START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF], or for bio-medical applications [START_REF] Claus | Contact modeling from images using cut finite element solvers[END_REF].

The extension of the implementation to the 3D case would obviously constitute an interesting prospect for this work. In this respect, the main issue may concern the surface integration over the non-conformal interfaces for the Nitsche coupling, which would require computing the intersection of a spline surface mesh with the interior of a spline volume mesh [START_REF] Antolin | Region extraction in mesh intersection[END_REF]. The extension of the formulation to mixed-dimensional problems, for instance, to incorporate beam models within volume meshes [START_REF] Kerfriden | A mixed-dimensional cutfem methodology for the simulation of fibre-reinforced composites[END_REF], may also appear desirable to cover a larger variety of structures, such as encountered in civil engineering. Finally, starting from images, e.g. obtained with X-ray computed tomography, and making use of specific image-based spline procedures [START_REF] Passieux | Image-based isogeometric twins of lattices with virtual image correlation for varying cross-section beams[END_REF], would be attractive to take into account the real inner geometries of materials.

Figure 1 :

 1 Figure 1: Schematic representation of the reference problem: two domains, Ω m and Ω f , are coupled through interface Γ.

Figure 2 :

 2 Figure 2: Illustration of the analysis-suitable model issue when applying IGA to the reference multiscale problem of Fig. 1. (a) Initial CAD parametrization; (b) Spline re-parametrization to obtain a boundaryfitted analysis-suitable model; (c) Discretizations considered in the case of immersed methods. In Fig. (c) the meshes are associated with the basis functions of the mechanical fields, while the gray areas correspond to the integration domains.

Figure 3 :

 3 Figure 3: Introducing conformal layers into the initial problem. Now, each subdomain consists of two, Ω li and Ω b i (i ∈ {m, f }) for the layer and the bulk, respectively. The possible non-linear coupling between the two layers is performed through interface Γ (still in red) that becomes conformal. The non-conformal aspect is moved inside the subdomains: for i ∈ {m, f }, non-conformal but perfect interfaces Γ i appear (depicted in green) to connect Ω l i and Ω b i . The LaTIn method will be applied for the coupling at interface Γ, while a Nitsche coupling scheme will be used at interfaces Γ i , thereby arriving at a robust, simple, and efficient strategy.

Figure 6 :

 6 Figure 6: Bilinear law considered for the cohesive interface to model delamination.

Figure 9 :

 9 Figure 9: Horizontal (left) and vertical (right) displacements plotted over the deformed geometry with a scale factor of 10.

Figure 10 :

 10 Figure 10: σ xx stress (on the left) and σ yy stress (on the right). Deformed geometry with a scale factor of 10.

Figure 11 :

 11 Figure11: Integration points on the interface: in contact (orange) and not in contact (grey) (left); normal traction on the interface, α is the angle between the position vector of a point and the positive direction of x axis (right).

Figure 12 :

 12 Figure 12: Meshes considered for convergence study: matrix model (black) and inclusion model (red).

( a )

 a Energy norm error for splines of second and third degrees computed in the whole domain.(b) Energy norm error for splines of second and third degrees computed only in the matrix background mesh (matrix without its conformal layer).

Figure 13 :

 13 Figure 13: Convergence plots for 2-and 3-degree splines.

Figure 14 :

 14 Figure 14: Description of the matrix with multiple inclusions problem: 1, 2, 4, and 16 inclusions. Black color corresponds to the matrix mesh, orange -to the inclusions meshes. Model parameters are given at the bottom.

Figure 15 :

 15 Figure 15: σ xx stress (on the left) and σ yy stress (on the right) for the model with 4 inclusions. Deformed geometry with a scale factor of 10.

Figure 16 :

 16 Figure 16: Integration points on the interfaces: in contact (orange) and not in contact (grey). Model with 4 inclusions.

Figure 17 :

 17 Figure 17: Convergence of the LaTIn algorithm for the multiple inclusions problem: models with 1, 2, 4 and 16 inclusions.

Figure

  Figure18: Model description for the frictional contact problem. Black colour corresponds to domain Ω 1 , orange -to domain Ω 2 . Note that due to the geometric simplicity of domain Ω 2 , the latter is directly made conformal to the frictional interfaces, so only one layer is built for Ω 1 at each of these interfaces.

Figure 19 :

 19 Figure 19: Reaction force (in N) on interface 3 for variation in friction coefficients µ 1 and µ 2 .

Figure 20 :

 20 Figure 20: Horizontal displacement at the end of the last time step. Deformed geometry with a scale factor of 50. Four pairs of coefficients considered. Note that for the two central cases (pairs (µ 1 = 0, µ 2 = 0.6) and (µ 1 = 0.6, µ 2 = 0)) the displacement varies throughout interface 3 (with the maximum equal to gap j) which means that only a part of interface 3 is in contact.

Figure 21 :

 21 Figure 21: Delamination: (a) problem description and geometrical parameters; (b) spline discretizations with a zoom on the non-matching interface.

Figure 22 :

 22 Figure 22: Reaction force versus applied displacement at each loading step.

Figure 23 :

 23 Figure 23: Damage parameter on the last time step on the interface integration points.

Figure 24 :

 24 Figure 24: σ yy stress in MPa for the last time step. Deformed geometry with a scale factor of 2.

Table 1 :

 1 18: Model description for the frictional contact problem. Black colour corresponds to domain Ω 1 , orange -to domain Ω 2 . Note that due to the geometric simplicity of domain Ω 2 , the latter is directly made conformal to the frictional interfaces, so only one layer is built for Ω 1 at each of these interfaces. Model parameters for the frictional contact problem. friction coefficients on the reaction forces at interface 3. We varied µ 1 and µ 2 from 0 to 0.6 with 20 values for each coefficient. The model parameters are listed in Table

	Parameters	Values
	E	210 GPa
	ν	0.3
	h	50 mm
	e	10 mm
	j	0.04 mm
	d	5 mm
	F max 1 F max 2	30 MPa 30 MPa
	µ 1 , µ 2	∈ [0, 0.6]
	∆t	1 s
	η	10 -7

Funding and/or Conflicts of interests/Competing interests

The authors have no relevant financial or non-financial interests to disclose.


W Daniau

R Salut

J.-M Friedt
email: jmfriedt@femto-st.fr

A Peña

J Debray

B Menaert

P Armand

P Papet

A Penarier

P Nouvel

Experimental demonstration of Surface Acoustic Wave propagation on α-GeO 2 for wireless, passive sensor design

Keywords: GeO2, surface acoustic wave, resonator, passive wireless sensor

On the path towards developing passive transducers for wireless sensing in high temperature (>560 • C) environments, we experimentally demonstrate the design, manufacturing and interrogation of resonators patterned on GeO2 substrates.

I. INTRODUCTION

Passive sensors using the linearity of piezoelectric effects for far-field wireless cooperative target short range RADAR measurement have been investigated for probing quantities in environments incompatible with a wired connection, whether moving targets or harsh environmental conditions. High temperature monitoring matches the latter context, with most investigations focusing on langasite family materials concluding to the long-term instability of this subtrate when exposed to high temperature. In this investigation, a material of the same crystal class (point group 32) as α-quartz and without phase transition up to its melting temperature (1116 • C) -α-GeO 2 -is considered as the substrate for designing and realizing a passive transducer compatible with short range RADAR interrogation.

II. SURFACE ACOUSTIC WAVE PROPAGATION: SIMULATION

The propagation of Surface Acoustic Waves has been investigated numerically previously [START_REF] Taziev | SAW properties in quartz-like α-GeO 2 single crystal[END_REF]. In this work, the mechanical and electrical constants extracted from [START_REF] Lignie | Growth of piezoelectric water-free GeO 2 and SiO 2 -substituted GeO 2 single-crystals[END_REF] are used to define the mixed matrix acoustic wave propagation and reflection parameters [START_REF] Ballandras | Simulations of surface acoustic wave devices built on stratified media using a mixed finite element/boundary integral formulation[END_REF]. Based on this analysis, and considering aluminum electrodes selected for their adequation to cleanroom manufacturing technologies for repeating multiple prototyping steps on a given substrate, the electromechanical coupling coefficient map, acoustic velocity map, and beam steering map are computed.

The electromechanical coupling coefficient K 2 is observed (Fig. 1) to remain below the 1% range, with a maximum at l = t = 0 in the IEEE Stds 176-1949 orientation standards. For this propagation direction, the velocity is slowest around The beamsteering is negligible in this propagation direction (Fig 3), easing the design of a resonator, the only geometry compatible with delaying the transducer response beyond clutter for a wireless interrogation using a short range RADAR system.

From this acoustic propagation modelling, the (Y X) cut of α-GeO 2 is selected for manufacturing tranducers.

III. TRANSDUCER DESIGN AND MANUFACTURING

From the previous analysis, the low electromechanical coupling coefficient prevents the design of a reflective delay line, and a resonator has been designed with 78 interdigitated electrode pairs in the transducer and 80 electrodes in the Bragg mirrors on each side, in either synchronous and asynchronous designs aimed at maximizing the quality factor Q at resonance frequency f . Indeed the objective of separating the sensor response decaying as Q/(πf ) from clutter is to reach a Q-factor high enough for clutter to have faded out, with an objective in the Q ≃ 5000 optimum between sensor identification and measurement speed to reach τ ≥ 3.5 µs assuming a short range RADAR distance detection limit of 150 m.

Due to the minute dimensions of the available substrates and the challenges of patterning 1.15 µm-wide electrodes separated by 1.15 µm, electron-beam (e-beam) lithography is employed to pattern interdigitated transducers with 180 µm acoustic aperture, mirrors and contact pads. A total of 11 resonators are patterned on the 6 × 10 mm 2 samples to assess the validity of the material constants used in the simulation. Etching the aluminum electrodes leads to increasing surface roughness and only a few patterning and etching steps aimed at scaling the electrodes to reach the 434-MHz ISM band can be performed before the sample needs to be polished again.

IV. WIRELESS MEASUREMENT

Multiple resonators with quality factors in the 2000-3000 range were manufactured and probed under a probe-station. The microwave probes were fitted with an helical monopole antenna and wireless interrogation in a sub-meter range was achieved using a dedicated short range RADAR system based on a frequency swept pulsed measurement for scanning the returned power as a function of frequency. The interrogation range is at the moment limited by the quality factor of the resonator due to excessive degradation of the surface roughness after multiple processing steps (Fig. 4).

Nevertheless, multiple temperature sweeps were performed in the ambiant to 200 • C temperature range, allowing for a first estimate of the first and second order Coefficient of Temperature with Frequency CTF (Fig. 5). V. CONCLUSION An experimental demonstration of surface acoustic wave propagation on α-GeO 2 was achieved, allowing for wireless measurement of the passive transducer resonance frequency as a function of temperature at sub-meter range.

  This work is supported by the French National Research Agency ANR under the OVERHEAT grant 21-CE08-0017. The RENATECH network is acknowledged for supporting the MIMENTO cleanroom facility.

Fig. 1 .

 1 Fig. 1. Electromechanical coupling coefficient (%) of Rayleigh waves for α-GeO 2 (Y Xlt)/l/t as a function of cut angle l and propagation angle t.

  2100 m/s (Fig.2), leading to lithography resolution challenges when designing Ultra-High Frequency (UHF) range transducers, e.g. meeting the radiofrequency emission regulations in the narrowband European Industrial, Scientific and Medical (ISM) band centered on 434 MHz.

Fig. 2 .

 2 Fig. 2. Velocity map (m/s) of Rayleigh waves for α-GeO 2 (Y Xlt)/l/t as a function of cut angle l and propagation angle t.

Fig. 3 .

 3 Fig. 3. Beam steering map ( • ) of Rayleigh waves for α-GeO 2 (Y Xlt)/l/t as a function of cut angle l and propagation angle t.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Wireless sensing of the GeO 2 resonator shown in inset.


Vincent Maurice

Clément Carlé

Shervin Keshavarzi

Kumar Ravinder

Samuel Chutani

Ludovic Queste

Jean-Marc Gauthier-Manuel

Rémy Cote

Moustafa Abdel Vicarini

Rodolphe Hafiz

Boudot

Wafer-level vapor cells filled with laser-actuated hermetic seals for integrated atomic devices

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.