Sylvain Chateau
email: sylvain.chateau@usherbrooke.ca

Jean Mercat

Zhe Li

Sébastien Poncet

Julien Favier

Umberto D'ortona

EMERGENCE OF METACHRONAL WAVES IN CILIA ARRAYS: A HYDRODYNAMIC MECHANISM

The present work reports the emergence of metachronal waves in cilia arrays immersed in a two-fluid environment using a coupled lattice Boltzmann -Immersed Boundary method. The periciliary layer (PCL) is confined between the wall and the mucus layer. Its depth is chosen in such a way that the tips of the cilia can penetrate the mucus layer. The cilia are initially set in a random state but quickly synchronize with their immediate neighbors with a phase shift giving birth to sympleptic or antipleptic metachronal waves, depending on the strength of the fluid retroaction onto the cilia. Antiplectic waves are found to be the most efficient to transport mucus compared to other random or synchronised cilia motions.

INTRODUCTION

Fluid propulsion by moving cilia is an universal phenomenon that can be found everywhere in Nature from the locomotion of micro-organisms to multiple processes in biological organisms at the cellular scale. For examples, in the early human embryonic development, cilia are responsible for the heart placement on the left side of the chest. Cilia move also the ovules in the Fallopian tubes, or the nutriments in the brain. In the specific problem of mucociliary clearance, cilia are responsible for transporting the mucus outside the lungs. In severe respiratory diseases, such as cystic fibrosis or asthma, the number of cilia clusters is decreased compared to the case of a healthy patient, and many cilia may beat in an abnormal way resulting in a less efficient transport of mucus. It has been experimentally observed [START_REF] Sleigh | The biology of Cilia and Flagella[END_REF] that cilia beat in a synchronised way with a phase shift between two neighboring cilia, forming a metachronal wave which appears to greatly enhance the transport of mucus. Then, it is of great importance to understand how this wave emerges and which type of waves provides the most efficient transport. This work could result in a better understanding of the mucociliary clearance process.

NUMERICAL METHOD

The fluid part is solved on a Cartesian grid using the lattice Boltzmann method with a D3Q19 scheme. The model of Porter et al. [START_REF] Porter | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios[END_REF] is used to model the fluid-fluid repulsion forces. The cilia are modeled with a set of 20 Lagrangian points, and the immersed boundary method is used to ensure the no-slip condition along the cilia [START_REF] Li | An improved explicit immersed boundary method to couple with Lattice Boltzmann model for single-and multi-component fluid flows[END_REF]. The equations of motion for the cilia are taken from Chatelin [START_REF] Chatelin | Méthodes numériques pour l'écoulement de Stokes 3D : fluides à viscosité variable en géométrie complexe mobile; application aux fluides biologiques[END_REF] and the motion of each cilium is decomposed into a finite number of steps during a period. If needed, an interpolation between two steps can be done in order to have the right velocity values along the cilia. The code is suitable for High Performance Calculation using a MPI parallelization. Cilia are set in arrays in a two-fluid layer, with PCL at the bottom and mucus above it. The PCL depth is set in such a way that the tips of the cilia emerge into the mucus layer. Periodic boundary conditions are used in the X and Y directions. A Bounce-Back Rule condition is imposed at the bottom of the domain and a free-slip velocity condition at the top using the "mirror method". The fluid retroaction onto the cilia is evaluated by projecting the immersed boundary forces calculated at each time step for every Lagrangian points onto the corresponding velocity vectors, and then the torques are computed. By doing so, only the norm of the velocity vector, but not its direction, is modified. A new parameter α controls then the intensity of the retroaction:

|| V || = || V || ± α|| dV ||

RESULTS

Starting from an initially random state, cilia quickly synchronize with their immediate neighbors with a phase shift. With time, metachronal waves -antipleptic or sympleptic-emerge. The sign of the parameter α, and the ratio h/L cil of the PCL depth with regards to the cilia length are two key parameters that determine the kind of emerging waves. Positive (resp. negative) α-values correspond to sympleptic (resp. antipleptic) waves. A proper value of h P CL /L cil is needed for the metachronal waves to emerge. As it can be seen on figures 1(a) and 1(b), cilia move along a plane and so can be subjected to collisions. But, since a 3D configuration is considered, one can suppose that cilia would actually slip onto each other in reality. Figure 2 shows the mean mucus velocity evolving with time. It can be seen that, for the same set of parameters, antipleptic metachronal coordination appears to greatly enhance the transport of mucus compared to random or other synchronized Spacing between two cilia is 7 lattice units, and h P CL /L cil = 0.9. Size of the domain in X, Y, Z is 449, 8, 50. (b) Antipleptic metachronal wave emerging from an array of 16 cilia disposed along the X direction, with α = -3.5. Spacing between two cilia is 5 lattice units, and h P CL /L cil = 0.6. Size of the domain in X, Y, Z is 81, 6, 50. In both cases (a) and (b), L cil = 22 lattice units and the viscosity ratio is set to 15. collective motions. As for sympleptic metachronal coordination, the mean mucus transport is weaker that in the case of cilia beating randomly. We explain this by the fact that this kind of synchronization induces large recirculations in the mucus phase. Moreover, the cilia tips are always in the PCL phase and never emerge into the mucus. Confirming the simulations of Ding et al. [START_REF] Ding | Mixing and transport by ciliary carpets: a numerical study[END_REF] using the regularized Stokeslet method, antiplectic waves are the most efficient for mucus transport.

CONCLUSIONS

An efficient and validated 3D lattice Boltzmann solver coupled to the immersed boundary method has been applied for the first time to investigate the emergence of metachronal waves and their superiority to transport mucus in a two-phase flow environment. A simple and efficient method has been applied successfully to take into account the fluid retroaction onto the cilia. Metachronal waves emerge quickly before the flow reaches a steady state. The efficiency of the metachronal waves for the mucus transport has been calculated, and antipleptic waves are found to be the most efficient ones compared to other random or synchronised motions.

Figure 1 :

 1 Figure 1: (a) Sympleptic metachronal wave emerging from an array of 64 cilia disposed along the X direction, with α = 3.Spacing between two cilia is 7 lattice units, and h P CL /L cil = 0.9. Size of the domain in X, Y, Z is 449, 8, 50. (b) Antipleptic metachronal wave emerging from an array of 16 cilia disposed along the X direction, with α = -3.5. Spacing between two cilia is 5 lattice units, and h P CL /L cil = 0.6. Size of the domain in X, Y, Z is 81, 6, 50. In both cases (a) and (b), L cil = 22 lattice units and the viscosity ratio is set to 15.

Figure 2 :

 2 Figure 2: Mean mucus velocity for different kind of coordinations.

The authors acknowledge the support of the Physio-Assist company. The PhD of Sylvain Chateau is funded through the NSERC discovery grant program.


T Hirschler
email: thibaut.hirschler@insa-lyon.fr

R Bouclier

A Duval

T Elguedj

J Morlier

The Embedded Isogeometric Kirchhoff-Love Shell: from Design to Shape Optimization of non-conforming stiffened multipatch Structures

Keywords: Shape optimization, Shell, Mortar coupling, Isogeometric Analysis, Multipatch, Stiffened structures

Isogeometric shape optimization uses a unique model for the geometric description and for the analysis. The benefits are multiple: in particular, it avoids tedious procedures related to mesh updates. However, although the analysis of complex multipatch structures now becomes tractable with advanced numerical tools, isogeometric shape optimization has not yet been proven to be applicable for designing such structures. Based on the initial concept of integrating design and analysis, we develop a new approach that deals with the shape optimization of non-conforming multipatch structures. The model is built by employing the Free-Form Deformation principle. Introducing NURBS composition drastically simplifies the imposition of the shape updates in case of a non-conforming multipatch configuration. In the case of stiffened structures, the use of embedded surfaces enables to tackle the geometric constraint of connecting interfaces between the panel and the stiffeners during shape modifications. For the analysis, we introduce the embedded Kirchhoff-Love shell formulation. The NURBS composition defines the geometry of the shell while the displacement field is approximated using the same spline functions as for the embedded surface. We also formulate a new mortar method to couple non-conforming Kirchhoff-Love shells which intersect with any angle. We apply the developed method on different examples to demonstrate its efficiency and its potential to optimize complex industrial structures in a smooth manner.

Introduction

Structural shape optimization requires a suitable mix of an accurate geometric description and an efficient analysis model. Even more importantly, a close link between both the geometric and the analysis models is highly sought since they repeatedly communicate during the resolution. Thus, shape optimization is a great application of IsoGeometric Analysis since the latter integrates Computer-Aided Design and analysis [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]. As a result, isogeometric analysis optimization is an active research topic and it has been successfully applied in solid mechanics; see for example Wall et al. [START_REF] Wall | Isogeometric structural shape optimization[END_REF], Qian [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF], Nagy et al. [START_REF] Nagy | Isogeometric sizing and shape optimisation of beam structures[END_REF][START_REF] Nagy | Isogeometric design of elastic arches for maximum fundamental frequency[END_REF][START_REF] Nagy | Isogeometric design of anisotropic shells: Optimal form and material distribution[END_REF], Kiendl et al. [START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF], Taheri and Hassani [START_REF] Taheri | Simultaneous isogeometrical shape and material design of functionally graded structures for optimal eigenfrequencies[END_REF], Wang and Turteltaub [START_REF] Wang | Isogeometric shape optimization for quasi-static processes[END_REF], Fußeder et al. [START_REF] Fußeder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF], Kang and Youn [START_REF] Kang | Isogeometric shape optimization of trimmed shell structures[END_REF], Lian et al. [START_REF] Lian | Shape optimization directly from CAD: An isogeometric boundary element approach using T-splines[END_REF], Wang et al. [START_REF] Wang | Isogeometric shape optimization of smoothed petal auxetic structures via computational periodic homogenization[END_REF], Choi and Cho [START_REF] Choi | Constrained isogeometric design optimization of lattice structures on curved surfaces: computation of design velocity field[END_REF], Lei et al. [START_REF] Lei | Shape Optimization for Natural Frequency with Isogeometric Kirchhoff-Love Shell and Sensitivity Mapping[END_REF], Hirschler et al. [START_REF] Hirschler | Isogeometric sizing and shape optimization of thin structures with a solid-shell approach[END_REF], Ding et al. [START_REF] Ding | A gradientbased shape optimization scheme via isogeometric exact reanalysis[END_REF], Weeger et al. [START_REF] Weeger | Isogeometric shape optimization of nonlinear, curved 3D beams and beam structures[END_REF]. A general procedure, which has been improved over the years [START_REF] Daxini | Parametric shape optimization techniques based on Meshless methods: A review[END_REF], is commonly adopted. It is based on a multilevel design concept which consists in choosing different refinement levels of the same NURBS-based geometry to define both optimization and analysis spaces [START_REF] Nagy | Isogeometric design of anisotropic shells: Optimal form and material distribution[END_REF][START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF][START_REF] Wang | Isogeometric shape optimization for quasi-static processes[END_REF][START_REF] Hirschler | Isogeometric sizing and shape optimization of thin structures with a solid-shell approach[END_REF]. Shape updates are represented by altering the spatial location of the control points, and in some case the weights [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF][START_REF] Nagy | Isogeometric design of anisotropic shells: Optimal form and material distribution[END_REF], on the coarse level. The finer level defines the analysis model and is set to ensure good quality of the solution. The optimization and analysis refinement levels are independently chosen which provides a problem-adapted choice of the spaces.

From this overview, it appears that isogeometric shape optimization has not yet been applied to realworld structures. More precisely, mainly simple geometries modeled as single patch structures have been considered. However, it is well-known that a single NURBS patch cannot represent large and complex ge-ometries with for example holes and kinks. For instance, we are interested here in stiffened structures which are omnipresent in aeronautics [START_REF] Totaro | Recent advance on design and manufacturing of composite anisogrid structures for space launchers[END_REF][START_REF] Shroff | Design, analysis, fabrication, and testing of composite grid-stiffened panels for aircraft structures[END_REF][START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust nonintrusive coupling[END_REF]. These structures are obtained by the assembly of a main part (usually denoted by the skin or the panel) with subparts as stringers and ribs (the stiffeners). The fastening of the stiffeners with the skin gives a final structure with a high stiffness to weight ratio. Since stiffened structures are built with multiple parts, it seems obvious that multipatch modeling will be required to design them. Even more, the position and the shape of the stiffeners can be arbitrary which makes it difficult to model the whole structure with multiple conforming NURBS patches. Preprocessing procedures seem, in this case, inevitable to generate analysis-suitable conforming NURBS models [START_REF] Al Akhras | Towards an automatic isogeometric analysis suitable trivariate models gen-erationApplication to geometric parametric analysis[END_REF]. In this work we are interested, for example, in finding the optimal position of the stiffeners along the skin [START_REF] Zhang | A Moving Morphable Component Based Topology Optimization Approach for Rib-Stiffened Structures Considering Buckling Constraints[END_REF][START_REF] Wang | Streamline stiffener path optimization (SSPO) for embedded stiffener layout design of non-uniform curved grid-stiffened composite (NCGC) structures[END_REF]. In this case, the geometry of the skin remains unchanged during the whole optimization, and therefore it would be attractive not to modify its discretization. In addition, the generation of conforming meshes is not straightforward and may require time consuming procedures. This kind of framework has already been developed in the FEM community, see for example: Mulani et al. [START_REF] Mulani | EBF3PanelOpt: An optimization framework for curvilinear blade-stiffened panels[END_REF], Singh et al. [START_REF] Singh | An Optimization Framework for Curvilinearly Stiffened Composite Pressure Vessels and Pipes[END_REF], Zhao and Kapania [START_REF] Zhao | BLP Optimization of Composite Flying-wings with SpaRibs and Multiple Control Surfaces[END_REF]. It can be noted that these frameworks lack of compactness in the sense that several softwares are used during the different steps of optimization: geometry parametrization, meshing, analysis. Thus, we are looking for an optimization method that enables to use non-conforming multipatch models in order to suppress the tedious generation of conforming parametrization models. However, the imposition of the shape updates in case of non-conforming patches requires specific attention since there are no explicit links between the patches. We will show in this work how the use of Free-From Deformation techniques [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF] helps to deal with this issue. Especially, it transfers the shape modification between the patches in a master-slave way.

Isogeometric analysis of non-conforming multipatch structures is an active field of research. The main interest is to be able to perform the analysis directly on CAD models which almost systematically present trimming and non-conforming parametrization [START_REF] Marussig | A Review of Trimming in Isogeometric Analysis: Challenges, Data Exchange and Simulation Aspects[END_REF][START_REF] Teschemacher | Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis[END_REF][START_REF] Guo | Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format[END_REF]. The coupling can be addressed by different techniques. In case of arbitrary non-conforming interfaces, the coupling conditions are imposed in a weak sense. There are different approaches that can be employed; namely penalty methods [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF][START_REF] Breitenberger | Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures[END_REF][START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF], mortar methods [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF][START_REF] Brivadis | Isogeometric mortar methods[END_REF][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][START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: Geometric details, local refinement, inclusion, fracture[END_REF][START_REF] Dornisch | Dual and approximate dual basis functions for B-splines and NURBS Comparison and application for an efficient coupling of patches with the isogeometric mortar method[END_REF][START_REF] Sommerwerk | Reissner-Mindlin shell implementation and energy conserving isogeometric multi-patch coupling[END_REF][START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF][START_REF] Wunderlich | Biorthogonal splines for optimal weak patch-coupling in isogeometric analysis with applications to finite deformation elasticity[END_REF], and Nitsche methods [START_REF] Nguyen | Bonisoli, Nitsche's method for two and three dimensional nurbs patch coupling[END_REF][START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF][START_REF] Guo | A parameter-free variational coupling approach for trimmed isogeometric thin shells[END_REF][START_REF] Nguyen-Thanh | Isogeometric analysis of largedeformation thin shells using RHT-splines for multiple-patch coupling[END_REF][START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF]. Because stiffened structures are built with thin shells, we seek to consider isogeometric Kirchhoff-Love shell formulation. Moreover, this rotation-free shell formulation is very attractive for shape optimization because of its low computational cost [START_REF] Nagy | Isogeometric design of anisotropic shells: Optimal form and material distribution[END_REF][START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF][START_REF] Hirschler | Isogeometric sizing and shape optimization of thin structures with a solid-shell approach[END_REF]. In the specific context of Kirchhoff-Love shells, penalty [START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF] and Nitsche [START_REF] Guo | Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format[END_REF] methods have been successfully applied for large and complex industrial structures, even for non-linear analysis. Mortar coupling of Kirchhoff-Love shells has been done for geometric continuum surfaces [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF] and for non-linear analysis [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF][START_REF] Duong | A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries[END_REF]. However, the formulation in the context of linear elasticity analysis of multipatch Kirchhoff-Love shells which connect with arbitrary angles does not seem to have been performed yet. A solution can be found in [START_REF] Duong | A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries[END_REF] but the applicability in the context of linear analysis may not appear obvious. We present in this work a new mortar method applicable in this context. Differences with the mentioned approaches mainly concern the kinematic constraints imposed by the Lagrange multipliers. From an implementation point of view, the major benefit of the mortar method is its close link with domain decomposition techniques [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method?part I: A faster alternative to the two-level FETI method[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF][START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF][START_REF] Kleiss | IETI Isogeometric Tearing and Interconnecting[END_REF]. We discuss in this paper the advantages brought by such techniques to optimize multipatch structures.

Based on these observations, we develop in this work an approach which integrates the design and analysis into a single model accordingly to the isogeometric concept. Our approach deals with non-conforming multipatch models which are inevitable for designing stiffened structures. We offer a smooth way to impose the shape updates during the optimization even for complex geometries with varying curvatures. Neither during the design, the analysis nor the optimization steps we introduce approximations, reparametrizations, and other tedious procedures. To do so, we first tackle the geometric challenges in section 2. We present the idea of embedded surface related to Free-Form Deformation (FFD) methods [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF]. This will lead us to formulate, in section 3, the embedded isogeometric Kichhoff-Love shell. Then we present our mortar approach for the coupling and we remind basics on domain decomposition methods for the resolution. In section 4, we highlight the benefits and the potentials of using NURBS composition to perform shape optimization. Finally, section 5 presents advanced numerical examples that demonstrate the applicability and the efficiency of our approach to design innovative stiffened structures. This brings us to section 6 where concluding remarks and perspectives are discussed.

Geometric modeling

We first need to deal with the geometric difficulties of imposing shape updates to complex geometries. In par- ticular, this work focuses on stiffened structures which are obtained by the assembly of a main part (the panel) with one or more specific substructure (the stiffener). The main part and its substructure are coupled at a common interface. We introduce in this section the use of embedded surfaces to tackle this geometric constraint of connecting interfaces during shape modifications.

Geometric challenges

The simple case of a portion of a cylinder with one stiffener helps to highlight the geometric challenges. A first configuration can be defined as depicted in figure 1. Two classical NURBS surfaces describe the portion of cylinder and the stiffener. During the design, or more specifically, the optimization process, the designer may want to modify the shape of the different parts. In the IGA-based framework, it consists in changing the control point coordinates and, in few cases, their associated weights [3-8, 11, 14, 17]. Figure 2 shows two difficulties that arise with non-conforming multipatch configurations. The first difficulty occurs when trying to modify the shape of the stiffener. We want the bottom of the stiffener to perfectly lie on the cylinder. This geometric constraint is not easily fulfilled as shown in figure 2(a). In other words, sliding the control points along the edge of the cylinder disconnects the surfaces. The second difficulty is to modify the global shape (i.e., in this case the portion of cylinder) while keeping the substructure connected. More precisely, the question here is: how to transfer the shape update from the global part to its substructure? An illustration of the problem is shown in figure 2(b). If we move one control point of the main surface, the stiffener does not lie anymore on it. In this simple case, one could formulate an appropriate condition in order to glue the stiffener but, in general, it is not a trivial task. In case of curvilinear stiffeners lying on a curved surface, reparametrizations and approximations appear inevitable.

Embedded surface

Main idea

We introduce the concept of embedded surface in order to overcome the geometric challenges highlighted in the previous section. The idea comes from Bauer et al. [START_REF] Bauer | Embedded structural entities in NURBS-based isogeometric analysis[END_REF] who considered univariate embedded entities. More precisely, they embed a curve into a surface. This technique allows them to model cables in membrane structures or to apply line loads and supports to a shell structure without being limited to the edges of the surface.

A similar procedure is undertaken to define the stiffener for the stiffened cylinder depicted in figure 1. We do not use classical NURBS surfaces to describe the stiffeners. Instead of a simple surface, we introduce the composition of a NURBS volume and a NURBS surface. In other words, we immerse a surface into the parametric space of a volume. The volume can be seen as a mapping that smartly transforms the embedded surface in order to create a final surface with the desired properties. Figure 3 illustrates this process.

In general, embedded entities (curve, surface or volume) can be introduced when the shape modification of a given structure is controlled by another one or by sophisticated geometric rules. In this work, we restrict ourselves to embedded surfaces without loss of gener- ality. The embedded entities can be interpreted as slave parts since their shapes are driven by external geometric properties. Conversely, the master parts are those which dictate specific constraints during the shape update. For example, in the case depicted in figure 1, the portion of cylinder is the master part and the stiffener constitutes the slave part: modifying the master part imposes to transfer shape modification to the slave part and, in turn, the possible shape modifications of the slave part are restricted to follow the master geometry.

Mathematical description

Let us now introduce the following quantities. We define the mapping V as a NURBS volume:

V(ξ 1 , ξ 2 , ξ 3 ) = i jk R i jk (ξ 1 , ξ 2 , ξ 3 )P i jk , (1)
where R i jk are trivariate NURBS basis functions and P i jk are the related control points. The parametric coordinates ξ 1 , ξ 2 and ξ 3 are defined over the parametric domain Ω. The embedded surface S is a NURBS surface defined as:

S(θ 1 , θ 2 ) = ab Rab (θ 1 , θ 2 )Q ab , (2)
with Rab some bivariate NURBS basis functions and Q ab their associated control points. Under the consideration that the codomain of S is included in domain Ω, we finally compose the two previous mapping (1) and

(2) to get the phyical surface:

S(θ 1 , θ 2 ) = V S(θ 1 , θ 2 ) . (3)
The surface S is parametrized by the same set of variables (θ 1 , θ 2 ) as the embedded surface S and it returns value in the physical space R 3 .

Design benefits brought by embedded surfaces

The potential of using embedded entities is multiple from the design point of view. Our first motivation is to design stiffened structures found in aeronautics but, of course, it should not be limited to this purpose.

To demonstrate the flexibility of the proposed approach, we apply it to design the internal substructure of a curved wing. The wing is described by its external geometry as shown in figure 4(a). The outer skin of the wing can, for example, be defined by two classical NURBS surfaces. In order to build the internal substructure, we firstly introduce a NURBS volume that fills the inside of the wing. The two outer surfaces compose two opposite faces of the volume. Degree one is chosen between these faces. Thus, the control points associated to the volume are none other than the union of those associated to the both outer surfaces. Figure 4(b) shows the defined volume. Once the volume is generated, embedded surfaces are introduced into its parametric space. The composition of these embedded surfaces and the volume gives surfaces that are located inside the wing and that are perfectly coincident at some of their edges with both outer skins. One can simply modify the shape of the substructure by changing the location of some control points of the embedded surfaces. As shown in figure 5, a uni-directional translation of the embedded surface in the parametric space of the mapping leads to a smooth translation of the substructure which follows the curvature of the wing. It also automatically takes The skin of the wing can be defined by two NURBS surfaces (red and blue patches). These surfaces are used to define a NURBS volume that fills the inside of the wing. Then, one can embed surfaces into its parametric domain to describe the internal substructure (see for example figure 5).

1D translation

Complex motion

Initial Surfaces New Surfaces into account the cross-section variations of the wing. The geometric constraints given by the curvature and the cross-section variations could be really tedious, if not impossible without approximation, to be imposed with classical NURBS patches. The example of the wing is only one example of the capabilities of using embedded entities from the design point of view. We believe that such techniques can be successfully applied to a large range of design problems. The examples tackled in this article illustrate some of the applications (see more specifically section 5).

Shell analysis with embedded surfaces

The analysis of non-conforming multipatch structures is of high interest in the IGA community (see [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF][START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF][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][START_REF] Dornisch | Dual and approximate dual basis functions for B-splines and NURBS Comparison and application for an efficient coupling of patches with the isogeometric mortar method[END_REF][START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF][START_REF] Guo | A parameter-free variational coupling approach for trimmed isogeometric thin shells[END_REF][START_REF] Coox | A robust patch coupling method for NURBS-based isogeometric analysis of non-conforming multipatch surfaces[END_REF] to name a few). In this section, we firstly adapt the popular isogeometric Kirchhoff-Love shell formulation to the case of embedded surfaces. Secondly, we propose a mortar approach for the coupling of non-conforming parametrizations and we remind the basics of domain decomposition methods which enable to efficiently solve mixed problems.

Embedded Kirchhoff-Love formulation

The surface obtained by the composition of NURBS functions, as introduced in section 2, is not a classical NURBS surface. Thus, the standard isogeometric Kirchhoff-Love shell requires some adjustments to be used with embedded surfaces.

Looking at Kiendl et al. [START_REF] Kiendl | Isogeometric shell analysis with KirchhoffLove elements[END_REF], it can be seen that the starting point of the Kirchhoff-Love element formulation is the definition of the covariant base vectors. Taking equation ( 3), we observe that the surface is parametrized and, therefore, it has curvilinear coordinates which are θ 1 and θ 2 . As a consequence, the definition of the covariant base vectors is straightforward and is given by the chain rule:

a α = ∂S ∂θ α = k=1,3 S,θ α • e k V ,ξ k S α = 1, 2 , (4)
where the subscript (•) ,ξ = ∂(•)/∂ξ indicates the partial derivative with respect to variable ξ. The unit normal vector a 3 is obtained using the relation:

a 3 = 1 J a 1 × a 2 , where J = |a 1 × a 2 |. (5)
The Kirchhoff kinematic assumptions are wellknown: straight lines normal to the mid-surface are characterized by rigid-body motions and remain normal to the mid-surface after deformation. Such assumptions enable to describe the deformation of the shell body only by the displacement u of the mid-surface. Therefore, the overall displacement field U of the shell body can be written as: [START_REF] Nagy | Isogeometric design of elastic arches for maximum fundamental frequency[END_REF] where ζ ∈ [-t 2 ; t 2 ], t being the thickness of the shell. From Echter et al. [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF], the linearized rotation vector Φ is defined as a function of the covariant base vectors and rotation angles:

U(θ 1 , θ 2 , ζ) = u(θ 1 , θ 2 ) + ζ [Φ × a 3 ] (θ 1 , θ 2 ),
Φ = ϕ 1 a 1 + ϕ 2 a 2 , (7)
where the rotation angles ϕ 1 and ϕ 2 are:

ϕ 1 = 1 J u ,θ 2 • a 3 and ϕ 2 = - 1 J u ,θ 1 • a 3 . (8)
It follows that the linearized strain tensor of the shell body is found to be of the form

ε αβ = e αβ + ζκ αβ , (9)
where membrane strains e and bending strains κ are given by

e αβ = 1 2 u ,θ α • a β + u ,θ β • a α , (10)
κ αβ = -u ,θ α θ β • a 3 (11) + 1 J u ,θ 1 • a α,θ β × a 2 + u ,θ 2 • a 1 × a α,θ β + a 3 • a α,θ β J u ,θ 1 • a 2 × a 3 + u ,θ 2 • a 3 × a 1 .
Greek indices take the values 1 and 2. In fact, the kinematic assumptions make the transverse shear strains vanish (.i.e. ε α3 = 0). Looking at equation [START_REF] Fußeder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF] reveals that the derivatives of the covariant base vectors versus the curvilinear coordinates are required. In case of the composition of NURBS, the expressions of these derivatives are obtained starting from equation ( 4) and applying one more time the chain rule. It gives:

a α,θ β = k=1,3 S,θ α θ β • e k V ,ξ k S ( 12
)
+ k=1,3 l=1,3 S,θ α • e k S,θ β • e l V ,ξ k ξ l S .
One can notice that second derivatives are employed for the bending term. More specifically, the physical surface obtained by NURBS composition has to be C 1 (see equation ( 12)). This requirement can be satisfied with NURBS discretizations of higher continuity. If both the embedded surface S and the volume mapping V are C 1 then the NURBS composition is de facto C 1 . Applying Galerkin's method raises the question of the choice of the approximation space V h for the displacement field. In fact, the isoparametric concept may not be applicable in case of a NURBS composition. We need to define the solution space as all linear combinations of a given set of basis functions. Since second derivatives of the displacement field are involved in the bending term (see equation [START_REF] Fußeder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF]), the approximated displacement field has to be C 1 as the geometry. Two choices emerge:

1. based on the trivariate NURBS discretization of the mapping:

u h (θ 1 , θ 2 ) = i jk R i jk S(θ 1 , θ 2 ) U i jk , (13)
2. based on the bivariate NURBS discretization of the embedded surface:

u h (θ 1 , θ 2 ) = ab Rab (θ 1 , θ 2 ) U ab . (14)
The first choice provides an immersed-like approach where the deformation of the shell is prescribed by the surrounding volume [START_REF] Rank | Geometric modeling, isogeometric analysis and the finite cell method[END_REF][START_REF] Legrain | A NURBS enhanced extended finite element approach for unfitted CAD analysis[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][START_REF] Hansbo | Cut Finite Element Methods for Linear Elasticity Problems[END_REF]. The second choice is close to the initial Kirchhoff-Love isogeometric shell.

The mapping only plays a role from the geometric point of view since it modifies the shape of the embedded surface but it is not involved in the solution space. Depending on the application, one can adopt either one

θ 2 θ 1 ξ 3 ξ 1 ξ 2 x y z x y z FFD Step V(ξ 1 ,ξ 2 ,ξ 3 ) S(θ 1 ,θ 2 ) + U(θ 1 ,θ 2 )
Figure 6: Overview of the embedded isogeometric Kirchhoff-Love shell elements. The displacement field is approximated by the same parametrization as the embedded surface. However, the embedded surface does not describe the geometry of the shell. An additional mapping step further modifies the geometry of the shell. The deformed configuration is obtained by adding the displacement field to this NURBS composition.

or the other approach. The first approach can be interesting if the volume which surrounds the surface is a structural part. In case of void, a lot of degrees of freedom may be inactive leading to an ill-conditioned system, which is a well-known drawback of immersed method. Special treatments, such as penalization techniques [START_REF] Hansbo | Cut Finite Element Methods for Linear Elasticity Problems[END_REF][START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method[END_REF], have to be performed. It is particularly true for higher continuity discretizations since the supports of the NURBS basis functions become larger. The second approach leads to smaller system because it uses a bivariate parametrization (instead of a trivariate one for the first approach). No regularization is required and higher continuity discretizations (cubic and quartic for example) can be used. We notice that the second approach is also easier from the implementation point of view, especially during the stiffness matrix assembly.

In this work, we apply the second methodology ( 14). An overview of the approach is depicted in figure 6. Compared to the classical isogeometric Kirchhoff-Love shell, the difference lies in the additional trivariate mapping step that modifies the shape of the undeformed configuration. We name it the FFD step since it reminds the Free Form Deformation technique [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF]. Finally, introducing the approximation u h into (10) and [START_REF] Fußeder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF] gives the discrete membrane and bending strains. The equilibrium configurations of the shell followed from the principle of minimum potential energy leads to the typical linear system:

K U = F. ( 15
)
Further details on the expression of the stiffness matrix K and the vector of external force F can be found in Kiendl [START_REF] Kiendl | Isogeometric Analysis and Shape Optimal Design of Shell Structures[END_REF] and for example in the following papers; Hirschler et al. [START_REF] Hirschler | Isogeometric sizing and shape optimization of thin structures with a solid-shell approach[END_REF], Cirak et al. [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF], Kiendl et al. [START_REF] Kiendl | Isogeometric Kirch-hoffLove shell formulations for general hyperelastic materials[END_REF]. Here, one has to replace the expressions of the covariant vectors and their derivatives as given by equations ( 4) and (12).

Non-conforming coupling

We now introduce a new mortar method for the coupling of non-conforming isogeometric Kirchhoff-Love shells. Similarities with already published papers in the standard FEM community can be found, especially with those of Bernadou et al. [START_REF] Bernadou | Numerical analysis of junctions between plates[END_REF][START_REF] Bernadou | Numerical analysis of junctions between thin shells Part 1: Continuous problems[END_REF] where theoretical studies are carried out. In the isogeometric framework, mortar coupling of Kirchhoff-Love shells has been recently studied by Apostolatos et al. [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF] and Duong et al. [START_REF] Duong | A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries[END_REF]. Differences with our approach mainly lie in the kinematic constraints imposed by the Lagrange multiplier. In Apostolatos et al. [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF], the rotation vector is not clearly defined and the continuity of the rotation is imposed in the three directions of the global basis, which is only true in case of coplanar domains [START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF]. In Duong et al. [START_REF] Duong | A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries[END_REF], the chosen kinematic constraints deal with arbitrary shell junctions but they are only suited and usable for non-linear problems.

Mortar approach for Kirchhoff-Love shells

For the sake of simplicity, the method is presented in the case of two non-overlapping subdomains Ω 1 and Ω 2 without loss of generality. At the common interface denoted by Γ, the continuity of the displacement has to be ensured. Furthermore, in case of a shell junction with rigid hinge, an additional constraint is required which enforces the continuity of the rotation in the tangential direction associated to the interface curve. Thus, we formulate the kinematic constraints as:

u 1 = u 2 on Γ [3 displacements], (16) Φ 1 • t = Φ 2 • t on Γ [1 rotation], ( 17
)
where t is a unit tangent vector associated to the interface curve. The constraints are ensured in a weak sense by introducing two Lagrange multipliers λ ∈ L d and µ ∈ L r (L d and L r being ad-hoc spaces for the displacements and the rotation respectively). It results the formulation of the following Lagrangian for the coupled problem:

L (u 1 , u 2 , λ, µ) = 1 2 a 1 (u 1 , u 1 ) -l 1 (u 1 ) + 1 2 a 2 (u 2 , u 2 ) -l 2 (u 2 ) + λ, u 1 -u 2 + µ, Φ 1 • t -Φ 2 • t , (18)
where bilinear forms a s and linear forms l s constitute the standard variation forms of the elasticity problem on each subdomain. Operator •, • is a bilinear form defined such that:

v, w = Γ v • w dl . (19)
Then, the variational principle written in the discrete form gives the coupled linear system to be solved:

                         K 1 0 C T 1 Z T 1 0 K 2 C T 2 Z T 2 C 1 C 2 0 0 Z 1 Z 2 0 0                                                   U 1 U 2 Λ M                          =                          F 1 F 2 0 0                          , (20)
where vectors U 1 , U 2 , Λ, and M collect the Degrees Of Freedom corresponding to the discretizations of u 1 , u 2 , λ, and µ respectively. The displacement coupling matrices C s and the rotation coupling matrices Z s consist in sparse rectangular operators. As for the stiffness matrices K s , they are obtained by assembling local matrices C loc and Z loc that can be expressed at the control point level as:

C loc = ± Γ R λ R I d dl, (21)
Z loc = ± Γ R µ a 3 J R ,θ 2 (a 1 • t) -R ,θ 1 (a 2 • t) dl . ( 22
)
Here, we voluntarily skip the indices to ease the reading. R λ and R µ correspond to some basis functions of the displacement and of the rotation Lagrange multipliers respectively. R is a basis function of one displacement field. I d is the 3-by-3 identity matrix. Finally, C loc is a 3-by-3 matrix and Z loc is a 1-by-3 matrix. The symbol ± indicates whether the domain is the master (+) or the slave (-). The coupled linear system ( 20) is a saddle point problem. A special care may be required for the construction of the approximation subspaces of the Lagrange multipliers to avoid undesirable energy-free oscillations due to the non-satisfaction of the inf-sup condition. Let p denote the smaller degree of both subdomain displacement fields. We adopt the following strategy:

• for the displacement constraint ( 16), a B-Spline function λ h with degree p -1 is defined since it is mainly related to traction forces,

• for the rotation constraint ( 17), a B-Spline function µ h with degree p -2 is defined because it transfers a bending moment,

• same mesh refinement is chosen for both Lagrange multipliers λ h and µ h . We discretize these fields using as many elements as the coarsest of the domains Ω 1 and Ω 2 over the interface.

At this stage, it is important to say that this strategy is only based on numerical experiments. With such a choice, we never encountered instabilities in our computations. To illustrate this point, two classical test cases, the Scoderlis-Lo roof and the T-beam, are presented in the next section. Although the numerical studies validate numerically the formulation, we notice that a strong mathematical investigation of the choices of the discrete fields could be of interest from a theoretical point of view. For example, based on the work of Brivadis et al. [START_REF] Brivadis | Isogeometric mortar methods[END_REF] and Wunderlich et al. [START_REF] Wunderlich | Biorthogonal splines for optimal weak patch-coupling in isogeometric analysis with applications to finite deformation elasticity[END_REF] one could propose appropriate spaces to reach an optimal convergence rate for the Kirchhoff-Love shells. One may note that concerning stability issues, the penalty method [START_REF] Breitenberger | Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures[END_REF][START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF] and the Nitsche approach [START_REF] Guo | Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format[END_REF][START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF] are thus attractive since they involve simpler stability parameters.

Numerical investigation of the coupling

We shortly present two test cases to assess the presented coupling strategy. The goal is to underline the applicability of the approach for well-known problems.

The Scoderlis-Lo roof is discretized with three nonconforming embedded patches. Starting from the configuration depicted in figure 7 -3.005925 × 10 -1 (see Herrema et al. [START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF]). The results are given in figure 7(c). Good convergence of the solution is obtained. For the T-beam problem, we take the same numerical setting as in Herrema et al. [START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF]. Two planar patches are connected in a non-matching way.

The interface crosses the middle of the upper patch. We use here the classical Kirchhoff-Love shell to show that the proposed coupling strategy is fully applicable in this context. We obtain a similar deformation as in Herrema et al. [START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF] where penalty coupling was performed. In order to show that the rotation constraint is well prescribed, we plot in figure 8(c) the resulting angle between the patches along the interface. It can be seen that it remains equal to 90 • after deformation. Without the rotation constraint, the angle between the patches was found to be equal to 86.6 • at the end of the beam because the interface acts as a hinge [START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF]. The presence of small oscillations for the coarser parametrization highlights the non-optimality of the chosen spaces for the Lagrangian fields. It could be interesting to investigate the inf-sup condition to clarify the stability of the method [START_REF] Brivadis | Isogeometric mortar methods[END_REF]. However, these oscillations vanish with a finer parametrization.

In addition to these two examples, the following studies of section 5 will further confirm the potential of the proposed coupling approach.

Domain decomposition

A great benefit of mortar coupling is its link with domain decomposition methods [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method?part I: A faster alternative to the two-level FETI method[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF][START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF][START_REF] Kleiss | IETI Isogeometric Tearing and Interconnecting[END_REF]. The coupled problem (20) can be decomposed in order to solve only local and smaller systems. This is a particularly interesting technique to analyze large and complex structures characterized by the assembly of many substructures.

The starting point is to split the coupled system (20) into the following coupled equations:

U s = K -1 s F s -C T s Λ -Z T s M ∀s (23)
s

C s U s = 0 and s Z s U s = 0. (24)
The dual interface problem is obtained by the condensation of the problem at the interface. It consists in substituting equation ( 23) into equation [START_REF] Al Akhras | Towards an automatic isogeometric analysis suitable trivariate models gen-erationApplication to geometric parametric analysis[END_REF]. We get:

S d Λ M = T c T z (25)
where the condensed right hand side T and the so-called dual Shur complement operator S d read:

S d = s C s Z s K -1 s C T s Z T s , (26)
T c T z = s C s Z s • K -1 s F s . (27)
Now, the analysis of the non-conforming multipatch structure is done in two main steps: first we solve the interface problem [START_REF] Zhang | A Moving Morphable Component Based Topology Optimization Approach for Rib-Stiffened Structures Considering Buckling Constraints[END_REF], and secondly we deduce the local solutions from equation [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust nonintrusive coupling[END_REF]. One can see that these equations involve only local systems. Finally, a last ingredient comes into play to entirely decompose the problem (and make it parallelizable). The Shur complement is neither computed nor assembled. Indeed, the condensed symmetric linear system (25) is solved using a Krylov iterative solver, and more precisely a preconditioned conjugate gradient solver. With such a solver in hand, only matrix-vector products are required. As operator S d is expressed as a sum on subdomains, these products are computed as follows:

S d v = s C s Z s • y s , y s = K -1 s C T s Z T s • v. ( 28
)
The vectors y s are solutions of local systems. Thus, during the whole resolution, only local quantities are involved making the approach parallelizable. We do not go further into the details regarding domain decomposition in the present contribution. Other key aspects as efficient preconditionner, floating domain, and practical implementation of the method are not of interest here. The goal is to present the main concepts behind the domain decomposition methods in order to highlight the potential of using mortar coupling to analysis non-conforming multi-patch structures. For structural optimization, the domain decomposition method enables, in a preprocessing step, to factorize once for all the stiffness matrices of patches located in non-design regions. Following the invoked strategy, the solution at every iteration of the optimization process is therefore computed at a very competitive cost compare to the monolithic resolution [START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF].

Optimal design

Using IGA to perform structural optimization is attractive since it integrates design and analysis. We show in this section that introducing embedded surfaces into the IGA-based shape optimization framework offers new possibilities. It enables to tackle challenging optimization problems in a proper and accurate manner.

IGA-based shape optimization

A commonly adopted procedure is used in this work, and thus only the basics of IGA-based shape optimization are reminded. Interested readers can find more information in the following papers [START_REF] Nagy | Isogeometric design of anisotropic shells: Optimal form and material distribution[END_REF][START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF][START_REF] Fußeder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF][START_REF] Wang | Isogeometric shape optimization of smoothed petal auxetic structures via computational periodic homogenization[END_REF][START_REF] Hirschler | Isogeometric sizing and shape optimization of thin structures with a solid-shell approach[END_REF].

The main ingredient is the multilevel approach which allows to define the optimization model and the analysis model in a versatile manner. Both optimization and analysis models are initially obtained through krefinement of the initial geometry of the structure and represent the exact same geometry during the optimization process. The shape variations are imposed on the optimization model. Therefore, depending on the complexity of the admitted shape variations, one can choose an appropriate discretization level. A coarser optimization model may provide a simpler optimal shape than a finer optimization model. Then, the level of refinement for the analysis model is defined so as to ensure good quality of the solution. Thus, the analysis model is finer than the optimization model. Finally, thanks to the refinement procedure of NURBS functions, the link between both models is kept. The shape modifications imposed on the optimization model can be explicitly transferred to the analysis model by using the linear relation between the control points of these two discretizations [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF].

In this work, we focus on the common optimization problem that consists of minimizing the compliance under a given volume constraint. The mathematical formulation of this problem is:

min 1 2 F(s) • U(s) subject to s l ≤ s ≤ s u , V(s) ≤ V 0 , ( 29
)
where vector s collects the n s design variables. These variables take real values limited by lower bound s l and upper bound s u . In IGA-based shape optimization, the design variables are related to the control point coordinates in order to parametrize the geometry of the structure. Other quantities involved in the problem ( 29) are the load vector F, the state variables U, the volume of the structure V, and a prescribed volume V 0 . There are plenty of numerical methods to solve such an optimization problem. We use in this work the gradient based algorithm SLSQP [START_REF] Kraft | A Software Package for Sequential Quadratic Programming[END_REF] available in many scientific packages, as for instance SciPy [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF]. Computation of the sensitivities are required for specification of the search direction. Accurate sensitivities are crucial to obtain trustful results since they control the shape updates of the iterative optimization process. In this work, we manage to evaluate discrete sensitivities in a fully analytical fashion. This is made possible because of the close link between design and analysis within the isogeometric framework [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF][START_REF] Nagy | Isogeometric design of anisotropic shells: Optimal form and material distribution[END_REF][START_REF] Taheri | Simultaneous isogeometrical shape and material design of functionally graded structures for optimal eigenfrequencies[END_REF][START_REF] Lei | Shape Optimization for Natural Frequency with Isogeometric Kirchhoff-Love Shell and Sensitivity Mapping[END_REF].

L h q Figure 9:
The parabolic arch problem taken from Kiendl et al. [START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF].

When subjected to constant vertical load q, the optimal shape is a parabola with ratio h/L = 0.54779.

Design capabilities

Taking the example of the parabolic arch described in Kiendl et al. [START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF], we illustrate the design capabilities offered by embedded surfaces. We compare the resolution of this problem when using the classical Kirchhoff-Love shell and the introduced embedded version. The description of this benchmark example is given in figure 9. Under the consideration of a small thickness, the bending stiffness vanishes, and the optimal shape is analytically found to be a parabola with a specific length to height ratio (see again [START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF]).

One can solve the arch optimization problem by modelling the arch with isogeometric Kirchhoff-Love shell elements. In this case of a classical NURBS surface, the design variations are imposed in the vertical direction through design variables which move the inner control points of the beam, as depicted in figure 10. With the parameters given in [START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF], one has to find the final solution of a parabola with a height h opt = 5.4779 (see figure 10). Now, we solve the arch optimization problem with the embedded Kirchhoff-Love shell elements. The difference lies in the imposition of the shape variations. In fact, it can be done in two ways. Firstly, the design variables can be related to the control points of the volume mapping. We call this approach the FFD-like shape update because it recalls Free-Form Deformation techniques [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF]. Secondly, the design variables can be linked to the control points of the embedded surface in order to modify the shape of the shell. This second ap-

Random Initial Shape

Fixed Mapping

Movable CPs

Fixed CPs i=3 i=6 i=9 i=20

Optimal Shape i=0 proach is referred to the embedded shape update. Figures 11 and 12 illustrate the resolution of the arch optimization problems with these two techniques. Obviously, the parametrizations used here are not the unique way to solve the problem. For the FFD shape update, we define the embedded surface as a parabola with an arbitrary length to height ratio, because we know, in this particular case, that the final volume we are looking for is a linear mapping that corrects the ratio in order to obtain the expected parabola. Figure 11 shows the optimized solution. The final volume has a rectangular shape as predicted. It gives a parabola with the correct height. Note that another initialization than a parabola could have been considered and also lead to the correct optimal design. In case of embedded shape updates, the volume is fixed. We define it sufficiently large to ensure that the design space contains the optimal solution. Figure 12 shows that the optimized shape corresponds once again to the analytical solution.

The potential of the two design techniques is significant. We present in next section a few examples that highlight inherent benefits. The FFD design helps to optimize the global shape of structures composed of multiple parts. Conversely, the embedded technique offers a smooth way to optimize specific parts of a global structure. In any case, main and specific parts remain perfectly connected through the optimization process.

Advanced numerical examples

We apply the presented approach to three problems. The stiffened roof illustrates the use of the FFD design technique. The global shape modifications of the roof are automatically transferred to the stiffeners. The curved wall problem shows how the design modifications of the embedded surfaces allows to optimize the shape and the position of stiffeners along curved and complex geometries. Finally, the wing problem proves the applicability of our approach to design innovative aeronautical structures.

Stiffened roof

The problem of the stiffened roof is derived from the initial shape optimization problem of a roof under vertical uniform load [START_REF] Bletzinger | Computational methods for form finding and optimization of shells and membranes[END_REF][START_REF] Kegl | Shape optimization of truss-stiffened shell structures with variable thickness[END_REF]. Two stiffeners are added to the initial square plate. We introduce two embedded surfaces to define these stiffeners as shown in figure 13. The construction can be divided into two steps. First, the volume mapping is generated by the extrusion of the square plate. We set the degree to one in the direction of extrusion (.i.e. the knot vector in direction ξ 3 is defined as W = {0 0 1 1}). Thus, the generation of the volume is straightforward. The volume shares half of its control points with those associated to the square plate. The additional control points are a copy of the previous ones where the z-component is set to an offset value. The second step consists in embedding surfaces into the parameter space of the created volume. These embedded surfaces, once composed with the volume mapping, give the final surfaces within the physical space that define the stiffeners. Here, the embedded surfaces are B-Splines of one single element with degree 2-by-1. Thus, each surface has six control points as depicted in figure 13. The coordinates of their control points can directly be interpreted from this figure considering the ratio In this problem, the embedded surfaces remain unchanged during the optimization. We apply the FFD design technique. The shape modifications of the stiffeners are only imposed by the design updates of the surface representing the roof. To this purpose, the design variables are associated to the control points of the mapping. The volume mapping is discretized in 4-by-4-by-1 elements. Degree 3 is taken in directions ξ 1 and ξ 2 while degree 1 is kept in direction ξ 3 (red points in figure 13). Thus, the optimization model has 45 design variables that move the control points of the mapping in the vertical direction z. The analysis model is defined such that the roof and both stiffeners are discretized into 1024 and 256 bi-cubic elements respectively. The area of the roof without the stiffeners is constrained to be lower than 110% of the initial area of the square plate (.i.e. lower than V 0 = 110). The shape evolution of the optimization model and the analysis model is presented in figure 14.

Figure 15 presents the optimal results. The final shape is obtained after 125 iterations of the optimization algorithm when setting the stopping criteria regard- ing the relative objective function to 1e-8. Figure 15(c) depicts the evolution of the relative compliance during the resolution of the optimization problem. The compliance is drastically reduced: final relative compliance is equal to C opt /C 0 = 2.217e-3. The optimized structure deforms much less than the initial stiffened plate. The final displacement field is presented in figure 15(b). One can notice that the angle along the interface between the roof and the stiffeners varies. Thanks to our coupling approach, this type of configuration can be dealt with.

Looking at the continuity of the displacement field at the interface further confirms the suitability of the presented approach (see again figure 15(b)). It is interesting to note that the final shape has two planes of symmetry (see figure 15(a)). This is expected since the problem presents these symmetries. Therefore, some design variables have identical optimal values. The fact that we obtain the symmetries without setting groups of design variables is a meaningful indication to validate the result. From an engineering point of view, this problem of the stiffened roof highlights an interesting issue. We manage to solve the optimization problem with the substructure. One could also simplify the problem by ignoring in the optimization process the stiffeners assuming that they only have small influence on the behaviour of the global structure. We try to solve in a first step the optimization problem without the stiffeners, and we build the substructure afterwords. We compute the compliance on this new global structure and we get a compliance approximately 15% higher than the one obtained when optimizing the structure with the stiffeners. It highlights that tacking into account, during the optimization, sub-parts as stiffeners, holes, and other geometric details can be essential to design even better structures.

Curved wall

The second problem addresses the optimization of the shape and the position of a stiffener lying on a curved surface that represents a wall. Contrary to the stiffened roof problem, we do not modify the geometry of the wall which is here the master part. For this problem, we focus on the substructure only. In order to design this substructure, we apply the embedded design technique presented in section 4.2, figure 12.

The geometric construction of the stiffened curved wall is done similarly as the stiffened plate presented in figure 13. Starting from the surface defining the wall, we create a volume by extrusion. Once the volume is obtained, a surface is embedded in the parameter space of the volume. The composition of the embedded surface with the volume mapping gives the surface defining the stiffener in the physical space. More specifically, the wall is defined by a quadratic NURBS surface. This surface originally represents a quarter cylinder. We move some control points to generate the final curved wall depicted in figure 16. The embedded surface is defined by a quadratic NURBS with 2-by-4 elements. Its initial shape is described in figure 17. In order to perform the shape optimization of the stiffener, a set of design variables is associated to the control points of the embed-ded surface. These design variables modify the crosssection and the position of the stiffener. With the chosen parametrization, six groups of four control points are spread out along the main direction of the stiffener. For each group of control points, we define four design variables

v i = [v i 1 , v i 2 , v i 3 , v i 4 ]
as explained in figure 17(c) (i identifies the group number). Thus, a total of 6 × 4 = 24 design variables is used for this example. The initial shape is obtained by setting the design variables of each group equal to v 0 = [0.3, 0.90, 0.06, 0.05]. The area of the stiffener is constrained to be lower than the initial one which is V 0 = 21.579. Additional inequality constraints are included to the optimization problem in order to prevent undesired shapes during the resolution:

v i 3 -v i 1 ≤ 0, v i 1 + v i 3 -1 ≤ 0, v i 4 -v i 3 ≤ 0. ( 30
)
The first two constraints ensure that the embedded surface remains in the parameter space of the volume mapping. The third constraint limits the curvature of the cross-section of the stiffener. The optimization results are given in figures 18 and 19. During the resolution, the stiffener moves along the wall until it is located at the middle. Due to the symmetry of the problem, this result was predictable. The optimal shape has this symmetry which demonstrates the quality of the result. Especially, the weak coupling between the non-conforming patches is well defined and apparently does not yield any significant inaccuracy. The cross-section of the stiffener is larger at the bottom than at the top of the wall. A large crosssection at the bottom improves the fixation of the overall structure. The cross-section also becomes larger where the deformations are critical. Reducing the deformation at the bottom decreases the global deformation of the structure 18(a). The compliance of the final multipatch shape is equal to c opt = 2.416e-3, which is 39% lower that for the initial configuration. More details on the optimal shape are given in figure 19.

Wing

The last example of the wing highlights the applicability of our approach to design real world structures. We apply the concept presented in figures 4 and 5. The goal is to optimize the substructure of a simplified wing. We build the substructure by embedding surfaces into the parameter space of the volume that fills the inside of the wing. More precisely, we define a stiffener as presented in figure 20. The optimization problem consists of finding the best position of the stiffeners in order to increase the stiffness of the whole wing. We study successively the case of 2, 4 and 8 stiffeners.

One border of the wing is clamped and a uniform pressure is applied on its top. The initial configurations are obtained with evenly distributed stiffeners. Figure 21(a) shows these initial configurations with 2, 4, and 8 stiffeners. The placement of the stiffener is crucial for the global behaviour of the complete structure. One can notice that for the initial configuration with 2 stiff-eners, the upper and lower skins deform locally under the loading. This deformation appears due to the complex curvature of the wing. A better placement of the stiffeners could limit this phenomena. Therefore, one can perform the optimization with the present approach. Moving the stiffener along the wing is simplified by the use of the embedded surfaces as explained previously (see figure 5). One single design variable is associated to each embedded surface. These variables modify the appropriate coordinate of the control points in order to impose uni-directional translation. Thus, we have the same number of design variables as the number of stiffeners (i.e. 2, 4, or 8). Simple inequality constraints are added to the optimization problem so as to maintain a minimal distance between the stiffeners and an ascending order. Finally, there is no volume constraint in this example. The optimal results are depicted in figure 21(b). The decrease of the compliance is significant, and particularly in case of a low number of stiffeners.

More precisely, the final values of the compliance are 14.9%, 15.1%, 7.69% lower than the initial ones for 2,4 and 8 stiffeners respectively.

Concluding remarks

We successfully applied the concept of Free-Form Deformation to design multipatch shell structures. It consists in geometrically defining the shell mid-surface by a NURBS composition of an embedded surface with a volume mapping. We illustrated this construction principle on multiple examples and we applied it for stiffened structures. We managed to directly use the NURBS composition for the analysis by introducing the embedded Kirchhoff-Love shell formulation. This means that the original idea behind isogeometric analysis is maintained since both the design and the analysis are performed on an common model. It is particularly attractive when it comes to structural optimization. This way, analytical sensitivity of the compliance is achieved in this work. Moreover, the Free-Form Deformation techniques offer a smooth and an effective way to impose complex shape variations. It eliminates sophisticated geometric rules, reparametrizations and other approximations one can introduce when designing large and complex multipatch structures. Thus, we were able, for example, to create a versatile model of a wing with its substructure. Not only the position of the stiffeners are easily modified but also the shape of the whole design of the wing can be parametrized. One can change the shape of the outer skins by acting on the control points associated to the volume mapping. Due to the use of embedded entities, the internal substructure of the wing automatically follows the update of the outer geometry. This can be a very useful tool in a Multidisciplinary Design Optimization (MDO) context [START_REF] Balesdent | A survey of multidisciplinary design optimization methods in launch vehicle design[END_REF][START_REF] Martins | Multidisciplinary Design Optimization: A Survey of Architectures[END_REF]. One can simultaneously optimize the outer skin regarding an aerodynamic criteria and the shape of the internal substructure regarding the structural behaviour. In the aero-structural optimization framework, it seems that mainly sizing variables such as the thickness are introduced into the design process [START_REF] Kenway | Multipoint High-Fidelity Aerostructural Optimization of a Transport Aircraft Configuration[END_REF][START_REF] Keye | Aero-Structural Optimization of the NASA Common Research Model[END_REF]. The use of an isogeometric model with embedded entities can be very attractive to enlarge the design space by adding shape variables as the position of the stiffeners [START_REF] Gillebaart | Geometrically consistent static aeroelastic simulation using isogeometric analysis[END_REF][START_REF] Dubois | Parameterization Framework for the MDAO of Wing Structural Layouts[END_REF].

Finally, we pointed out the advantage of adopting a mortar coupling for optimizing multipatch structures. We believe that domain decomposition techniques are of great interest to efficiently solve multipatch optimization problem. One can factorize once and for all the stiffness matrix associated to the unmodified patches and perform the analysis with the domain decomposition approach. In addition, the decomposition allows to distribute the patches amongst several processors running in parallel, at a very competitive cost.

Figure 1 :

 1 Figure 1: Construction of a portion of cylinder with one stiffener by using two non-conforming NURBS surfaces.

Figure 2 :

 2 Figure 2: Geometric difficulties to impose design modifications to a stiffened cylinder: (a) while moving, the stiffener disconnects with the cylinder, and (b) the stiffener does not follow the shape variation of the cylinder.

Figure 3 :

 3 Figure 3: Construction of a stiffened structure using embedded surfaces. The final surface describing the stiffener is obtained by the composition of a NURBS surface and a NURBS volume. This construction ensures that the stiffener is perfectly lying on the portion of cylinder.

Figure 4 :

 4 Figure 4: Example of a wing: (a) the geometry is defined by its outer surface. (b)The skin of the wing can be defined by two NURBS surfaces (red and blue patches). These surfaces are used to define a NURBS volume that fills the inside of the wing. Then, one can embed surfaces into its parametric domain to describe the internal substructure (see for example figure5).

Figure 5 :

 5 Figure 5: Design of the internal substructure of a wing. Embedding surfaces in the parametric domain of the internal volume of the wing enables to describe ribs. The designer can simply modify the shape of the internal substructure by modifying the embedded surfaces. Here, a one-dimensional translation of an embedded surface provides the motion of the substructure along the curvature of the wing.

Figure 7 :

 7 Figure 7: Scoderlis-Lo non-conforming coupling: (a) the geometry is defined by three surfaces (two cubics and one quartic) which are embedded into a volume describing the portion of cylinder. The non-conformity is generated by the insertion of the indicated knots. (b) The distribution of the z-displacement field obtained with the three times refined parametrization. (c) Convergence of the z-displacement at points A and B with the refinement of the mesh.

Figure 8 :

 8 Figure 8: T-Beam non-conforming coupling: (a) Description of the problem: the interface cuts the upper patch at the middle of some elements. The upper patch has one more element in the beam direction than the bottom patch. (b) Deformed configuration with scale factor of 10 and the magnitude of the displacement field. (c) Angle between the patches along the interface after deformation. Coarse case corresponds to the mesh depicted in (b). Fine mesh is three times refined.

Figure 10 :Figure 11 :

 1011 Figure 10: Isogeometric resolution of the arch optimization problem with classical Kirchhoff-Love shell elements. A set of design variables (3 in the current case) move couples of control points in the z-direction.

Figure 12 :

 12 Figure 12: Resolution of the arch optimization problem with embedded Kirchhoff-Love elements in case of embedded shape updates. The volume mapping is fixed and it limits the design space. The shape modifications are imposed by moving the control points associated to the embedded surface (i denotes the SQP iteration).

Figure 13 :

 13 Figure 13: Construction of the geometry of the stiffened roof: starting from the square plate, the first step consists in generating a volume by extrusion. Degree one is set in the direction of the extrusion. Once the volume is obtained, NURBS surfaces are embedded into the parameter space of the volume. Finally, the compositions of these embedded surfaces with the volume give the final surfaces that describe the stiffeners.

  d 1 /d 2 = 0.35. The other geometric and mechanical parameters of the problem are the following: Young's modulus E = 210 • 10 9 , Poisson's ratio ν = 0.30, plate length L = 10, shell thickness t = 0.1, stiffener height h = 0.6, vertical uniform load P = 1000. The 4 corners of the plate are fixed.

Figure 14 :

 14 Figure 14: Shape updates of the stiffened roof during the optimization. At left the optimization model where the color map indicates the amplitude of the shape variation. These shape variations are imposed through the control points of the volume mapping (red points) which move in the vertical direction. At the right are depicted the associate analysis model. The stiffeners are perfectly lying on the roof during the whole process.

Figure 15 :

 15 Figure 15: Optimization results for the stiffened roof problem: (a) final shape with the associated shape modification field, (b) the displacement magnitude for the optimal shape, and (c) the optimization history which depicts the evolution of the relative compliance and the relative volume during the optimization.

Figure 16 :

 16 Figure 16: Geometric construction of the curved wall: control points coordinates of the surface and the volume mapping.

Figure 17 :

 17 Figure 17: Initial configuration of the stiffened curved wall: (a) the geometry where the stiffener is obtained by NURBS composition, (b) the magnitude of the displacement, and (c) the definition of the design variables which are associated to the control points of the embedded surface.

Figure 18 :

 18 Figure 18: Optimization histories for the resolution of the curved wall problem: (a) evolution of the geometry and the displacement field, and (b) evolution of the compliance and the volume.

2 Figure 19 :

 219 Figure 19: Optimal configuration of the stiffened curved wall: (a) final geometry and (b) the magnitude of the displacement. (c) The final embedded surface: side, top views and 3D views. The 2D views give information on the coordinates of the control points.

Figure 20 :

 20 Figure 20: Definition of a single stiffener for the wing problem. A surface is embedded into the parameter space of the volume that fills the inside of the wing. The mapping transforms the initial surface into the one defining the stiffener in the physical space.

Figure 21 :

 21 Figure 21: Initial (a) and optimal (b) configurations for the wing problem. From left to right, case of 2, 4, and 8 stiffeners. The displacement color map is scaled differently for each number of stiffeners. The deformation factor is the same for each configuration. As expected, the wing deforms less with the optimal position of the stiffeners.