Marie Tirvaudey

Robin Bouclier

Jean-Charles Passieux

Ludovic

Ludovic Chamoin

Non-invasive implementation of nonlinear IsoGeometric Analysis in an industrial FE software

come

Introduction

Isogeometric Analysis (IGA) was introduced by Hughes et al. [Hughes et al. 2005] and later detailed by Cottrell et al. [Cottrel et al. 2009] in order to solve problems directly on the geometry extracted from a Computer-Aided-Design (CAD) model. The main idea is to use the same basis, e.g. Non-Uniform-Rational-B-Splines (NURBS) [Cohen et al. 1980, Piegl andTiller 1997] or T-Splines [START_REF] Bazilevs | [END_REF], for analysis as the one used to describe the geometry in CAD. Beyond the reinforced link between CAD and analysis, IGA turned out to be a superior computational mechanics technology [Evans et al. 2009, Lipton et al. 2010, Schillinger et al. 2013, Yin et al. 2015, Rauen et al. 2017]. However, despite its real enthusiasm in the computational mechanics community, the implementation of IGA in existing industrial codes still appears quite invasive, which limits its massive deployment in industry. Indeed, IGA comes with an intermediate space for the denition of the shape functions which questions the element-wise structure of most existing FE codes. Unlike the standard Finite Element Method (FEM) where each element has its own parametrization, the parametric space for NURBS is localized onto a patch that is composed of several elements. Some IGA implementations in commercial FE packages exist such as in LS-Dyna [Hartmann et al. 2011, Benson et al. 2013, Hartmann et al. 2016, Chen et al. 2016], Abaqus [Elguedj et al. 2012, Duval et al. 2015, Lai et al. 2017] or Radioss [Occelli et al. 2019], but it is still quite a few.

The concept of Bézier extraction has proved to be a milestone to tie IGA and FEM closer together. Originally introduced by Borden et al. [Borden et al. 2011] for B-Splines and NURBS, the concept has now been generalized to a large variety of advanced splines such as T-Splines [Scott et al. 2011, Scott et al. 2012], hierarchical B-Splines and NURBS [Schillinger et al. 2012, Hennig et al. 2016, Angella et al. 2018], hierarchical T-Splines [Evans et al. 2015, Chen andde Borst 2018], and LR B-Splines [Dokken et al. 2013]. Focusing on B-Splines and NURBS, the technology enables to formulate a smooth polynomial B-Spline (respectively NURBS) function in terms of C 0 polynomial Bernstein (resp. rational Bézier) functions. Nevertheless, it must be stressed, at this stage, that implementing this strategy still requires modications both at the elementary and at the assembly levels: (1) modify the standard FE shape functions subroutine to introduce a new element based on polynomial Bernstein (or rational Bézier) functions, (2) apply the extraction for each element and, (3) modify connectivity since it is based on the NURBS numbering.

An extension of the idea of Bézier extraction to standard nodal FE functions has been more recently proposed in Schillinger et al. [Schillinger et al. 2016]. It gave birth to the Lagrange extraction operator that directly links Lagrange nodal basis with smooth B-Spline basis. This new operator encompasses the Bézier operator and oers an alternative and simple implementation: it merely requires to evaluate the B-Spline basis functions at the nodal points associated to the Lagrange polynomials. The Lagrange extraction especially appears of great interest for geometries based on polynomial B-Splines. Indeed, it eliminates the need for introducing Bernstein polynomials as new shape functions, which makes it possible not to touch the FE subroutines of a standard FE code at the element level. The Lagrange extraction thus goes a step further to allow for the integration of IGA in an industrial FE code with minimal invasiveness. However, it has to be underlined from this overview that the procedure cannot be applied to the general situation of NURBS (or any rational splines). In this case, a family of rational Lagrange functions can be introduced and used to directly express the NURBS basis; the problem is that such functions (that integrate weights and a weight function) are not part of a typical FE code which is restricted to polynomial basis functions. Moreover, it may be noticed that only the potential of the methodology has been underlined in [Schillinger et al. 2016]. Up to now, it seems that only a very recent attempt that extends this framework to the real open source FE-code FEniCS can be listed [Kamensky and Bazilevs 2019].

If pure numerical eciency is the only criterion for qualifying an implementation of IGA, the best approach would certainly be to rewrite everything into a dedicated piece of code written using a low-level language (see, e.g., PetIGA [Dalcin et al. 2016] for high-performance IGA based on the PETSc library, [Chang et al. 2016] for an implementation using Fortran and [Nguyen et al. 2015] for global implementation aspects). In addition, let us mention that recent works to implement IGA eciently start with a modication of the standard looping over elements in the assembly process [Calabro et al. 2017], which is incompatible with the use of standard FE codes. All of these approaches may require signicant development eorts, especially whenever a new nonlinear constitutive law needs to be added. The approaches based on subroutines, that were mentioned before, represent also interesting alternatives, as constitutive laws of existing industrial codes can be reused. However, these routines depend on the considered code and still require a certain eort of implementation. This is why, in this study, it is chosen to reduce the human time required for the computer implementation (even if, obviously, this one may come with a slight increase of the computational cost).

In other words, this work aims at performing IGA in an available FE software in the least possible invasive manner. Starting with the Lagrange extraction technology, a global view is adopted and a full algebraic bridge that directly goes from Lagrange nodal polynomials to NURBS functions is formulated. It leads to a novel, alternative implementation procedure in which no modication of the whole FE routines is required. As a result, and compared to the current practice, no modication of the shape functions subroutine and of the assembly is performed. More precisely, the strategy simply consists of the projection of the FE linear system on an ad-hoc regular reduced basis whose algebraic construction is automatic and represents an objectively moderate implementation eort. In this way, the method solely uses features oered by most modern industrial FE codes, such as, for instance, the possibility to extract the (tangent) stiness operators. This appears of crucial interest regarding industrial FE environments since the access to the whole code is not needed and the FE routines that may be highly optimized are not touched. In that sense, the proposed method lies in the family of so-called non-invasive (also denoted non-intrusive) methods that has been recently applied to local-global coupling [Gendre et al. 2009, Bouclier et al. 2016, Bouclier et al. 2017], domain decomposition [Duval et al. 2016], contact problems [Oumaziz et al. 2017, Oumaziz et al. 2019], transient dynamics analysis [Bettinotti et al. 2014] and stochastic partial dierential equations [Chevreuil et al. 2013]. Their ultimate goal is then to tackle real industrial applications (see, e.g., [Guinard et al. 2018] for real aeronautical structures). The strategy is applied here to simply incorporate isogeometric capabilities in the industrial FE code Code_Aster, which is a familiar open source software package for numerical simulation in structural mechanics developed by the EDF R&D company [Code_Aster 2014].

The paper is organized as follows: Section 2 species the existing link between NURBS-based IGA and FEM by reviewing the Bézier and Lagrange extraction technologies. Section 3 is devoted to the development of the proposed implementation procedure. Then, for validation and demonstration purposes, some numerical experiments are rst carried out involving 2D and 3D elasticity in Section 4, before performing the more complex simulation of a nonlinear elastoplastic structure in Section 5. Finally, Section 6 concludes on this work.

2

The link between IGA and FEM

In this section, the fundamentals to tie IGA and FEM closer together are given. Although the Lagrange extractor encompasses the Bézier one, it is chosen here to start with the original Bézier version [Borden et al. 2011] before introducing the Lagrange transformation [Schillinger et al. 2016] given the importance of the concept of Bézier extraction in the state-of-the-art of IGA. For completeness, note that further details regarding NURBS and related algorithms can be found in [Cohen et al. 1980, Piegl andTiller 1997]. From a set of non-decreasing coordinates in the parametric space collected in the knot vector Ξ = ξ 1 , ξ 2 , ..., ξ n+p+1 , a basis of univariate B-Spline functions of order p, denoted by (N i,p ) i∈{1..n} , can be dened recursively using the Cox-de Boor formula [Cohen et al. 1980]. A B-Spline function of order p has p-1 continuous derivatives. Besides, if a knot ξ l has a multiplicity m l , the number of continuous derivatives at this knot will decrease by m l . A comparison between B-Spline and Lagrange shape functions is given in Figure 1. Then, a B-Spline curve C BS is constructed as follows:

C BS (ξ) = n i=1 N i,p (ξ)P i = P BS T N(ξ), (1)
where

P BS =   x 1 1 ..x d 1 ... x 1 n ..x d n   is the (n × d) matrix that collects the d coordinates of the n control points and N =   N 1 ... N n  
is the vector of the B-Spline basis functions. For a multi-dimensional space domain, the B-Spline basis functions are determined by performing the tensor product of 1D B-Spline functions, so that in 3D it yields:

N A = N i,p (ξ) × N j,q (η) × N k,r (ζ), (2)
for control point P A that corresponds to the i th , j th and k th control points in each univariate spatial direction.

NURBS.

The multivariate NURBS functions

R IGA = R IGA A n A A=1 are dened from a set of n A B-Spline functions {N A } n A
A=1 . In order to do so, it is necessary to introduce w IGA A n A A=1 as the weights of each of the control points P IGA = P IGA A n A A=1 associated to the NURBS entity. The rational functions read as follows:

R IGA A = N A w IGA A W IGA , where W IGA = n A A=1 N A w IGA A . (3)
Following the denition of the B-Spline entities in Eq. (1), NURBS objects are then dened such that:

V IGA = n A A=1 R IGA A P IGA A = P IGA T R IGA . ( 4
)
The NURBS functions exactly describe the geometry of conic sections. For illustration purpose, an example of a 2D circular beam is considered in Figure 1(a). Starting with a quadratic one-element mesh, one can perform k-renement [Cottrell et al. 2007] in the arc direction to obtain a quadratic mesh composed of two elements.

Remark 1 The present article is restricted to quadratic functions since almost all industrial FE codes do not go beyond second-order Lagrange nite elements. Nevertheless, we emphasize that the developed implementation could be directly applied to higher-order spline-based functions, provided that the corresponding higher-order nite elements are available in a standard FE code.

The Bézier extraction

The rst step to link IGA and FEM is to formulate the smooth polynomial B-Spline functions in terms of polynomials that are C 0 on the element edges. This is possible as the space generated by the B-Spline functions is included in the one generated by the C 0 functions. Such a link has been developed in [Borden et al. 2011] under the name of Bézier extraction which maps the Bernstein basis onto the B-Spline basis.

Bézier decomposition by knot insertion.

In order to form a structured C 0 Bernstein mesh from a smooth B-spline mesh, one simply needs to repeat all the inside knots of the knot vectors until they reach a multiplicity p (see Figure 1(b) for the example of a quadratic B-Spline curve with two elements). The geometry is preserved with the resulting C 0 Bernstein representation. It must be stressed at this stage that the advantage of Bernstein functions is that they exhibit an elementary structure which is similar to FEM (see Figure 1(c)).

Bernstein polynomial functions.

The Bézier transformation relies on the denition of Bézier curves [Bézier 1986] and Bernstein functions. Originally, a Bézier curve of order p is a linear combination of p + 1 univariate Bernstein functions

B( ξ) = B i,p ( ξ) p+1 i=1 associated with control points P BER = {P i,p } p+1 i=1
. Regarding the expression of the Bernstein functions, they are dened directly in the parent space since each Bernstein element has its own parametrization. As with B-Spline functions, multi-variate Bernstein functions are built by applying a tensor product between the univariate functions. From here on, the notations B and P BER are extended to denote the functions and control points of a multi-variate Bernstein mesh composed of several elements.

Bézier extraction operator.

By making use of the above knot-insertion process, a Bézier extractor matrix C can be constructed in order to ensure the link between the B-Spline functions and the Bernstein functions as follows:

N = CB.
(5)

In order to determine the positions of the Bernstein control points the equality between the expression of the B-Spline entity (see Eq. ( 1)) and the Bernstein one is used. By expressing the B-Spline functions using Eq. ( 5), it can be shown that P BER = C T P BS .

The Lagrange extraction

The Lagrange extraction operator directly maps a Lagrange nodal basis onto a B-Spline basis. In the following, the remaining link between Bernstein and Lagrange shape functions is rst established and then the direct construction of the Lagrange extractor is presented.

From Lagrange to Bernstein polynomials.

Denoting by L the classical FE Lagrange functions, the goal is here to build operator D LB that satises:

B = D LB L. (6)
Obviously, operator D LB is constructed from its elementary representation D e LB which is the same for each element. For illustration purpose, one-dimensional Bernstein and Lagrange shape functions of order 2 are plotted in Figures1(b) and 1(c).

To construct operator D e LB , one simply needs to express the Bernstein functions as a linear combination of the Lagrange functions at some interpolation points. Making use of the interpolatory property of the Lagrange functions, such an operator can be eciently constructed by evaluating the Bernstein functions at the nodal points associated to the Lagrange basis. Once operator D LB is built, the same treatment as with the Bézier operator can be performed to construct the FE nodes from the Bernstein mesh: P F E = D T LB P BER . Such nodes could be used to construct the input mesh for classical nite element codes. Obviously, these nodes interpolate the geometry. The B-Spline, Bernstein and Lagrange control meshes that generate the same two-element quadratic polynomial beam are plotted in Figures 2(b), 2(c) and 2(d) respectively. As expected, the Bernstein control points are far from the geometry in comparison with the Lagrange nodes that interpolate the geometry.

A direct link between Lagrange and B-Spline functions.

With previous operators C and D LB in hand, the computation of the Lagrange extraction operator D becomes straightforward:

N = DL with D = C D LB . (7)
Nevertheless, for better numerical eciency, the Lagrange extraction operator is never computed this way. Indeed, the same procedure as for the Lagrange-Bernstein operator D LB can be directly applied to the full Lagrange extraction operator D: it merely requires to evaluate the B-spline basis functions at nodal points. An ecient algorithm for this has been proposed in [Schillinger et al. 2016].

The extraction in case of NURBS

Until now, the link between IGA and FEM is restricted to geometries based on polynomial B-Splines. However, in general, spline discretizations consist of rational basis functions that integrate weights associated to control points. The Lagrange extraction can be extended to the case of NURBS. It consists of establishing a link between NURBS and rational Lagrange basis functions, as detailed in [Schillinger et al. 2016]. Briey, it can be extracted from Eqs. ( 3) and ( 7) that:

R IGA = W IGA DL W IGA , (8)
where W IGA is the diagonal matrix of NURBS weights and

W IGA = n A A=1 w IGA A
N A is the NURBS weight function. The NURBS weight function can be rewritten using the Lagrange basis as:

W IGA = n A A=1 w IGA A N A = (w IGA ) T N = (w IGA ) T DL = (D T w IGA ) T L = (w LAG ) T D = W LAG , (9)
where the weights associated to the rational Lagrange control points are:

w LAG = D T w IGA . ( 10
)
The rational Lagrange functions are then dened as follows:

R LAG = W LAG L W LAG , (11)
where W LAG is the diagonal matrix of the Lagrange weights. The link between NURBS functions and rational Lagrange functions is nally made using Eqs. ( 11) and (9) in Eq. ( 8). Consequently, a new extraction operator D W is created as described below:

R IGA = W IGA D(W LAG ) -1 R LAG = D W R LAG . ( 12
)
Therefore, the rational Lagrange control points depend on the NURBS control points:

P LAG = (D W ) T P IGA .
Let us emphasize at this stage that this link from rational Lagrange functions to NURBS is exact since both bases are rational and NURBS are of higher-order smoothness. However, one must keep in mind that the denition of these rational Lagrange functions requires, from standard Lagrange polynomials, the incorporation of Lagrange weights W LAG and of the Lagrange weight function W LAG . Those operations are not part of standard FE codes that are restricted to polynomial basis functions. As shown in Section 3, an additional work thus needs to be achieved to directly go from Lagrange polynomials to NURBS.

3

The non-invasive implementation

Let us recall that our interest is to implement IGA in an available FE software with a reduced level of invasiveness. The potential of the Lagrange extraction technology to allow for the implementation of IGA in standard FE codes has been underlined in [Schillinger et al. 2016]. Especially in case of geometries based on polynomial B-Splines, a strategy with minimal invasiveness has been derived. After exposing the current practice for setting up IGA using the previous extraction concepts, the proposed novel implementation procedure is presented.

3.1 A rst option 3.1.1 General case of NURBS.

The rational Lagrange geometry that has been created using D W (see Eq. ( 12)) has an elementary structure which makes it more likely to be implemented in a FE software. The resulting implementation scheme and the interaction between the initial NURBS mesh and the computed rational Lagrange mesh are described in Figure 3. On each branch, the quantities of interest needed to perform the next step are pointed out, for example the global and local rational Lagrange extraction operators D W and D e W , or the NURBS and rational Lagrange shape functions R IGA and R LAG . Moreover, the connectivity table used for the matrix assembly is indicated. From this procedure, the following modications to be done in the FE software can be listed:

1. Modify the standard FE shape functions subroutine to incorporate the Lagrange weights and the weight function and thus, construct the rational Lagrange shape functions from the existing Lagrange polynomials, as expressed in Eq. ( 11);

2. Apply the extraction for each element using operator D e W ; 3. Change the connectivity table in order to perform an IGA assembly.

Figure 3: General standard implementation procedure using the Lagrange extraction for NURBS: from a NURBS denition of the geometry and with the creation of operator D W , a rational Lagrange representation of the geometry is created. By modifying the shape functions subroutine, to express rational Lagrange functions, and applying the local Lagrange extraction, the IGA elementary stiness matrix is computed. Finally, the assembly subroutine needs to be modied to allow the use of the NURBS degrees of freedom to get the global NURBS stiness matrix.

Specic case of B-Splines.

In the particular case where the geometry is generated only with B-Spline polynomials (no rational basis functions), the standard FE subroutines do not need to be touched at the element level. In other words, point 1 above does not stand anymore. Indeed, after computing the standard FE elementary stiness matrices K e F E , only additional matrix-matrix products are required to obtain the corresponding IGA elementary matrices K e IGA :

K e IGA = D e K e F E D eT . ( 13
)
This is made possible since a direct link is established between the Lagrange polynomials that are present in standard FE software and the smooth B-Spline functions. As a result, part of our work now consists in elaborating a direct link between NURBS and Lagrange polynomials to meet such a non-invasiveness in the case of rational functions as well.

3.2 The proposed strategy 3.2.1 Principle.

As stated above, no modication of the whole FE routines is envisaged in this work (i.e., no modication into the shape functions subroutine (point 1 above) and concerning the assembly (point 2 above)). In order to do so, a global view of the extraction is adopted, as depicted in Figure 4. The path starting with a B-Spline mesh is exact and mainly consists of applying procedure exposed in [Schillinger et al. 2016], but in a global way. The path related to NURBS however requires the construction of an additional operator to go from polynomials to rational functions. Such a transformation cannot be exact since this is the space of the rational functions that includes the associated polynomials and not the other way around. A projection thus needs to be performed. For simplicity, the projection is performed on the Lagrange side, i.e., the idea is to project the rational Lagrange discretization onto the polynomial Lagrange space. Operator D LL is introduced as:

R LAG = D LL L, (14)
and its construction is explained in the next paragraphs.

Figure 4: Approach to link a NURBS mesh to a FE mesh using dierent global operators. Those operators are then used to recover the NURBS stiness matrix from the FE one computed using a classical FE software, taken as a black-box. Black arrows represent an exact link, gray arrows an approximation.

Once the nite element mesh is created from the IGA mesh in a pre-processing step, it is used as an input for the classical FE code in order to compute the FE stiness matrix. With the dierent elaborated operators, successive transformations are then performed to obtain the nal isogeometric stiness matrix, as follows:

K IGA = D W D LL K F E D T LL D T W . ( 15
)
This transformation is also applied on the right-hand side:

F IGA = D W D LL F F E . ( 16
)
Consequently, the system K IGA U IGA = F IGA can be solved to obtain the displacement U IGA . Finally, note that the resulting IGA displacement can be back-converted in terms of nodal displacements:

U F E = D T LL D T W U IGA , (17)
so that existing subroutines of the FE code can be used for post-processing.

Remark 2 Notice that Eqs. ( 15) and ( 16) can be interpreted as projecting the FE system onto the isogeometric basis. As a result, the developed algebraic bridge (D W D LL ) provides a new point of view on the relation between IGA and FEA: IGA can be interpreted as the projection of FEA on a specic regular reduced basis.

Simple approximation.

Acting at the Lagrange level for performing the projection between the rational and associated polynomial spaces oers the opportunity to follow a pragmatic yet accurate strategy. Indeed, it has to be emphasized that the control points of the rational Lagrange discretization interpolate the geometry. As a result, it is possible to simply consider that the position of the FE nodes P F E is exactly the same as the position of the rational Lagrange control points P LAG . In this case, D LL formally reads as the identity operator (denoted by I D in the following), so that the derived procedure does not add any extra-computational eort from the current practice. As an illustration, Figure 5 investigates the dierence in terms of geometry in case of a quarter circular beam. It can be observed that the approximation is already very accurate for a single element (Figure 5(a)) and, obviously, it is improved through the renement of the mesh since more interpolated control points are added (Figure 5(b)). The accuracy of the approximation technique can also be appreciated through Figure 5(d) that displays the geometry error in terms of the relative dierence in radius from the exact value evaluated along the circular annulus. Given this high accuracy, it is expected that the error related to the mapping D LL = I D appears largely insignicant compared to that associated to the nite element resolution of an underlying mechanical problem.

(a) Simple approximation for one element. (b) Simple approximation for two elements (in the arc direction).

(c) Geometry error in terms of the relative dierence in radius between the simple approximation and the exact circle.

(d) Weight evolution with renement. The weight of the rational Lagrange control point P LAG is equal to cos 2 (α/2), which quickly converges to 1 as the mesh is rened.

Figure 5: The simple approximation in case of a quarter circular beam. The dierence in terms of geometry quickly vanishes with the renement of the mesh.

An alternative interpretation can be given: the strategy simply consists of considering all the weights of the rational Lagrange discretization equal to one. Taking back the example of a quarter circular beam, the approximation can be illustrated as follows (see Figure 5(d)). As explained in the NURBS technology [Cottrel et al. 2009, Piegl andTiller 1997], the weight of the middle point of a quadratic circular one-element arc is equal to the cosine of half of the angle subtended by the arc (the weights of the two boundary control points being one). From here on, α denotes the angle of interest (see Figure 5(d) again). Making use of previous equality (10), then expressing operator D using (7), and nally noting that C = I D in case of a mesh composed of a single NURBS element, the weights associated to the rational Lagrange discretization of this circular arc can be computed as follows:

w LAG = D e 1D LB T w IGA , (18)
with

w IGA =   1 cos(α) 1   and D e 1D
LB the 1D version of D e LB (6), which ends up with the weight:

w LAG = 1 2 (1 + cos(α)) = cos 2 α 2 (19)
for the middle point of the rational Lagrange parametrization. Therefore, the good accuracy of the simple approximation is accounted for: when α decreases cos 2 (α/2) quickly tends towards one. Note nally that going from Bernstein polynomials to rational Bézier functions with the same treatment would lead to a convergence speed of cos(α) instead, and going from B-Spline polynomials to NURBS functions similarly would converge even less quickly.

Local least-square approximation.

As a more rigorous mathematical approach, a strategy based on a projection using a L 2 -norm is also proposed to approximate the rational Lagrange functions by means of linear combinations of Lagrange polynomials. The presentation of the method is restricted to the two-dimensional setting. Given the local structure of rational Lagrange functions as well as Lagrange polynomials, the methodology proceeds element-by-element for the construction of operator D LL , thus providing simplicity and limiting the computational cost. Note that similar works have been made in quite related contexts, such as to project C 0 discretizations to smoother basis (see, e.g., [START_REF] Thomas | [END_REF] and [Schillinger et al. 2016] that build Bézier and Lagrange projection, respectively). Here, the following simple treatment can be performed. To start with, the vertex Lagrange nodes of the considered element are taken equal to the corresponding vertex rational Lagrange control points. Then, in order to be directly compatible with the neighboring elements, the projection is formulated by writing a least-square problem on each element edge. As a result, the developed strategy does not require the assembly step encountered in common local least-square procedures (see, e.g., [Govindjee et al. 2012]) in order to combine in some way the dierent values obtained from the dierent elements for the same point. More precisely, let us consider an elementary rational Lagrange geometry expressed similarly as in Eqs ( 1) and ( 4):

S e LAG = P e LAG T R e LAG , ( 20
)
where P e LAG is a matrix that collects the positions of the control points related to the considered rational Lagrange element and R e LAG is the vector of the elementary rational Lagrange functions. The objective is to determine the positions of the equivalent Lagrange nodes P e F E that give an accurate approximation of the initial rational Lagrange geometry S e LAG . Considering an element edge dened by the same second coordinate η in the parent space, the quadratic functional J is minimized such that:

J = 1 2 1 ξ=-1 C e LAG ( ξ) -C e F E ( ξ) 2 , ( 21
)
to compute the equivalent Lagrange nodes related to the considered edge. Let us emphasize that the curves C e LAG ( ξ) and C e F E ( ξ) only involve the one-dimensional rational Lagrange functions and Lagrange polynomials, respectively. After dierentiation, we end up with the following linear system that relates each of the spacial coordinates of the edge Lagrange nodes (X e F E ) with that of the edge rational Lagrange control points (X e LAG ):

M e 1D X e F E = M e R X e LAG , ( 22
)
where M e 1D is simply the elementary mass matrix (with unit density) of the 1D Lagrange polynomials and M e R reads as follows:

           1 ξ=-1 L 1 ( ξ)R LAG 1 ( ξ) 1 ξ=-1 L 1 ( ξ)R LAG 2 ( ξ) 1 ξ=-1 L 1 ( ξ)R LAG 3 ( ξ) 1 ξ=-1 L 2 ( ξ)R LAG 1 ( ξ) 1 ξ=-1 L 2 ( ξ)R LAG 2 ( ξ) 1 ξ=-1 L 2 ( ξ)R LAG 3 ( ξ) 1 ξ=-1 L 3 ( ξ)R LAG 1 ( ξ) 1 ξ=-1 L 3 ( ξ)R LAG 2 ( ξ) 1 ξ=-1 L 3 ( ξ)R LAG 3 ( ξ)            . ( 23
)
From system ( 22), the position of the middle Lagrange node alone is actually picked since, as stated above, the vertex nodes are taken equal to the vertex Lagrange control points. The procedure is then repeated to all element edges. Unlike operator D e LB , note that D e LL is not the same for each element, as it depends on the weights of the control points of the considered rational Lagrange element. However, let us emphasize that those weights only appear on the right-hand side of the local least-square problem (through M e R , see Eq. ( 23)). As a consequence, M e 1D is analytically computed and inverted once and for all elements:

(M e 1D ) -1 =     9/2 -3/4 3/2 -3/4 9/8 -3/4 3/2 -3/4 9/2     . ( 24
)
The only requirement of the developed procedure is actually to compute M e R on each element, that is, to perform a numerical integration in the parent space of an element edge. This is straightforward using standard numerical analysis packages.

Still on a quarter circular arc, Figure 6 shows the dierence in terms of geometry encountered when applying the proposed local least-square procedure for a mesh composed of a single element. From a global point of view, this more rigorous alternative seems to produce equivalent results (at least in terms of geometry) as the previous pragmatic strategy that considers all the weights equal to one. This accounts for the validity of the previous simple approximation in practical (engineering) applications.

4

Examples with elasticity models

As a rst assessment of the method, numerical experiments with two and three-dimensional linear elasticity models are carried out in this section. For validation purpose, the rst test-case (for which a reference analytical solution is available) is computed using a Matlab homemade code. Then, for proof-of-concept, the open source package Code_Aster, which is developed by the EDF R&D company [Code_Aster 2014], is used as an industrial FE software for the second benchmark. It may be noticed that no unit is given for the following test-cases (as it is usually the case in the literature). The given values are those that are directly taken as input for the numerical process.

2D circular beam under end shear

The rst example consists of a 2D circular beam under plane stress subjected to end shear, as illustrated in Figure 7. The material properties are the following: Young Modulus, E = 10000 and Poisson ratio, ν = 0.25. A radial constant displacement u 0 = 0.01 is prescribed on the bottom edge of the beam. The closed form solution in terms of strain energy can be found in [Zienkiewicz et al. 2003].

In Figure 8, the evolution of the relative error regarding the strain energy is given as a function of the number of degrees of freedom. The relative error is computed as:

|E ex d -E f e d | |E ex d | , (25)
where E ex d denotes the reference exact strain energy and E f e d the strain energy of the discrete model. To start with, it can be noticed that the results for the 9-node element created with the simple or the leastsquare projections are similar. This behavior was expected because, as it has been previously illustrated from a geometrical viewpoint, the simple approximation was already good for one element in the arc direction. Here, with the rst renement used, there are already 6 elements in this direction; therefore, the error related to the mapping D LL = I D is largely insignicant compared to the associated FE error. The analysis performed with the FE mesh created by both of these approximation methods gives an error evolution similar to the FEM reference [Zienkiewicz et al. 2003]. Then, in order to obtain the IGA stiness matrix, and thus to recover the isogeometric Figure 8: Evolution of the relative energy error of the circular beam. Full FEM [Zienkiewicz et al. 2003] and IGA resolutions are used as references. On the one hand, the error is calculated by using the simple approximation between rational Lagrange functions and Lagrange polynomials (see section 3.2.2) for the resulting FE problem (curve "Simple projected FEM") and back-projected IGA problem (curve "Simple projected IGA"). On the other hand, this error is computed for the FE and IGA problems with a least-square approximation between rational Lagrange functions and Lagrange polynomials (see section 3.2.3, associated curves "Least-square projected IGA and FEM"). results, the FE stiness matrices are back-projected using the developed global bridge. The errors found with the simple and least-square projected IGA are similar to the error of a full IGA reference computation (the curves "IGA reference", "Simple projected IGA" and "Least-square projected IGA" are superimposed). This accounts for the accuracy of the developed global bridge in case of NURBS. Finally, one can notice that for a given mesh renement, the errors are about the same for FEM and IGA. The main dierence is the number of degrees of freedom which signicantly decreases when IGA is performed (the two curves are simply horizontally translated). This gap between the two curves illustrates the increased per-degree-of-freedom accuracy of IGA (reduction by a factor of more than 3 here) which has been deeply demonstrated in the literature (see, e.g., [Evans et al. 2009, Cottrell et al. 2006, Schillinger et al. 2013]).

Solid horseshoe

In this section, a solid horseshoe subjected to equal and opposite displacements on the top surfaces is used as a full 3D example (see, e.g., [Hughes et al. 2005] for similar computations). The geometry and the NURBS mesh associated with the solid horseshoe are given in Figure 9(a). Let us emphasize that it involves a multi-patch model which can be easily handled using the proposed methodology. The NURBS mesh is composed of 1152 elements of degree 2 with 7020 degrees of freedom. The material properties are the following: E = 3 × 10 7 and ν = 0.3. An initial displacement u 0 = y is applied to the top-left surface and a similar displacement u 0 = -y is applied to the top-right surface. In addition, the top surface displacement across the x-direction is locked whereas for the z-displacement only the corner points on each top surface are locked. This choice of boundary conditions makes the loading asymmetric. Using the developed global bridge, a corresponding nite element mesh, composed of hexahedral elements with 27 nodes, is created (see Figure 9(b)). The FE mesh is composed of 33507 degrees of freedom. It is used as an input mesh in the industrial FE software. After recovering and solving the IGA linear system from the FE one constructed in Code_Aster, the obtained IGA displacement eld can be back-converted in terms of nodal displacements (see Eq. ( 17)) and the pre-processing routines of the FE code can be used to plot the NURBS stress, see Figure 10(a). A zoom on the area of stress concentration of σ yy , due to the asymmetric loading, is made in Figure 10(b). Results are smooth and in good agreement with reference IGA computations [Hughes et al. 2005], which means that the developed global bridge does not introduce signicant errors for the projection.

Extension to nonlinear analysis

Making IGA more accessible for industrial applications requires not to touch optimized routines which perform the integration of the nonlinear mechanical behavior. In this part the case of a nonlinear problem is considered with minimal invasiveness (no modication into the nonlinear solver). The only requirement of the method is to be able to extract the tangent stiness and right-hand-side operators at each iteration of the nonlinear solver.

Implementation of the method

In the nonlinear case, the resolution of the tangent problem is performed within a Newton solver which is still realized by the commercial FE code. The method is thus non-invasive also with respect to the options of the nonlinear solver. Technically, it requires additional features of the FE code, which consist in being able to pause the nonlinear resolution in order to externalize only the resolution of the global tangent system; then, to re-inject the displacement eld solution; and nally, to restart the nonlinear resolution without any other external treatments. The method is implemented here using Code_Aster software which oers this possibility thanks to Python subroutines and the so-called STAT_NON_LINE in splitted commands solver. More precisely, Figure 11 shows how the method works in a nonlinear context. After the pre-processing that does not change with respect to the elastic case, the solver STAT_NON_LINE in splitted commands enters the Newton loop. The nonlinear behavior is integrated and the tangent problem assembled with the optimized Fortran Figure 11: Flowchart of the non-invasive and nonlinear implementation in an existing FE software, taken as a blackbox. Note that operator Π is introduced to dene the full bridge between IGA and FEA (i.e., D for B-Splines or D W D LL for NURBS). routines of Code_Aster. On the other hand, the resolution of the tangent linear system is replaced by a Python subroutine, which consists in solving its projection onto the subspace spanned by Π. Once this resolution is made, the IGA displacement solution vector is transferred on the FE space so that it can be re-introduced in Code_Aster in order to update the behavior. It is worth noticing that a second very short Python subroutine is needed for the calculation of the IGA residual from the FE one. Some elements of these subroutines are provided hereafter for the implementation in Code_Aster. For instance, the subroutine iga_solver.py is roughly organized as follows:

• get the Aster concept which corresponds to the matrix and the right-hand side and convert to numpy arrays: Kef = sparse.csr_matrix(STIFF.EXTR_MATR()) Fef = FORCE.EXTR_COMP().valeurs

• projection of the global system and resolution: Uiga = scipy.sparse.linalg.spsolve( Pi.T.dot(K.dot(Pi)) , Pi.T.dot(F) )

• reconstruction of the FE solution vector Uef = Pi.dot(Uiga)

This Python subroutine is called in the COMM le of Code_Aster using the command below: from iga_solver import * Uef=iga_solver(STIFF,FORCE)

Finally, in the COMM le, an Aster concept has to be built from the displacement Uef using CREA_TABLE(...) and CREA_CHAMP(...) functions.

Remark 3 The application programming interface of Code_Aster being Python, it is possible in this particular case to avoid reading/writing the linear operators on the drive, as it may be required with other industrial software. The operators are kept in memory in sparse format. The algorithm simply requires to project the FE tangent system and FE residual on an ad-hoc reduced basis. This operation is eciently performed using a dedicated sparse Library Scipy.Sparse. Note that this kind of projection is classically used in many model order reduction techniques (Proper Orthogonal Decomposition, Reduced Basis, see, e.g. [Kerfriden et al. 2012, Quarteroni et al. 2016, Chinesta et al. 2017]).

Remark 4 It may also be noticed that our method appears compatible with previous works performed regarding the parallel computation of nonlinear problems using Code_Aster (see, e.g, [Duval et al. 2016, Oumaziz et al. 2019]). As a result, the use of our strategy in combination with such works should allow the parallel computation of nonlinear IGA using Code_Aster.

Application to elastoplasticity

For illustration purpose, a 3D dog-bone sample in tension is considered in this section. Such an experimental test is often used to characterize the hardening of metallic materials [Mathieu et al. 2015]. The specimen is 100mm long, the ligament is 5 mm wide and the sample is 2.5 mm thick, as shown in Figure 12 (top). An elastoplastic constitutive behavior is considered with a Von Mises mixed hardening (see Figure 12 (bottom)). The parameters are the Young modulus E=22.13GPa and Poisson ratio ν=0.3. The Prager constant is set to 2200MPa and the nonlinear hardening is based on the tensile stress-strain curve given in Figure 12. The sample is subjected to a remote tension: the load goes from p = 0MPa to p = 80MPa within 10 (non-uniform) increments. With such parameters, the nonlinearity is easily computed using very few iterations of the Newton-based solver at each loading step.

Remark 5 To prescribe the Dirichlet conditions, Code_Aster uses the double Lagrange multiplier method (see [Charras et al. 1993]) to avoid losing the positive-deniteness of the operator. The associated Lagrange multipliers are used in the nonlinear solver of Code_Aster to update the behavior. Although this does not represent any theoretical limitation, for the sake of simplicity, a Neumann problem is considered here, so that the multipliers vanish.

The FE mesh is built from the CAD as described in Section 3. Even if, in this case, only one element was used along the sample thickness, the number of degrees of freedom is signicantly reduced: 6561 for the IGA versus 23409 for the associated nite element mesh, see Fig. 13.

The longitudinal component of the displacement eld solution is depicted in Figure 14 (left). The Von Mises stress eld obtained with the non-invasive IGA method (STAT_NON_LINE in splitted commands + subroutine Python) is also presented in Figure 14 (center). It can be seen that, in the region of the ligament, this eld has values above the yield stress (≈ 200MPa) which illustrates that in its central part, the specimen has undergone plastic deformations. This IGA stress eld is compared to the Von Mises stress eld obtained using the input FE mesh and the nonlinear FE solver STAT_NON_LINE of Code_Aster with the same mesh. Let us notice that no signicant dierence regarding the convergence speed of the nonlinear solver has been observed between the IGA and FE versions for this test-case. Since the converged solutions are very close, the relative discrepancy in Von Mises stress between the IGA and FE solutions is plotted in Figure 14 (right). Despite a reduction of more than 70% in the number of degrees of freedom, the solutions are very close (less than 2% of local mismatch while the two solutions come from dierent approximation subspaces), which conrms, in the nonlinear framework, the observations that had been made in elasticity. In order to quantitatively assess the proposed non-invasive IGA method with respect to the classical FE method, the tensile force-displacement curve is eventually plotted in Figure 15. Again, despite its much smaller dimension, the IGA subspace is able to provide a solution almost identical to that generated by FEM. This result is in line with an interpretation of the IGA as the projection of FEA on a specic regular reduced basis (see Remark 2).

Conclusions

In this work, a novel implementation procedure was developed that further simplies the integration of IGA in well-established FE environments, the ultimate goal being to make IGA more accessible for industrial applications. The attractive property of our approach is that it exhibits a reduced level of invasiveness. More precisely, in the proposed implementation, the whole FE routines remain completely untouched. As a result, the total FE code that may integrate complex, nonlinear numerical models and optimized routines can be used as a black box. The only requirements is that the commercial code (a) is able to output the tangent stiness and right-hand side operators and (b) has full quadratic elements (the 9-node quadrilateral element in 2D, or the 27-node cubic element in 3D) so that it becomes possible to recover and solve, in an external script, the isogeometric system from the FE one. In that sense, the proposed procedure lies in the family of non-invasive methods whose interest for the transfer of advanced technologies to the industrial community has been proved these last years [Duval et al. 2016].

The key aspect of our approach was to adopt a global vision regarding the link between IGA and FEM. Starting with a global formulation of the Lagrange extraction operator [Schillinger et al. 2016], it was shown that it is possible to perform non-invasive B-Spline based IGA in a real industrial FE code, which already constitutes a progress from a practical point of view. Then, to meet a similar reduced level of invasiveness in the case of rational functions as well, a global approximate link to go from standard nodal Lagrange polynomials to NURBS functions was developed. An additional operator that enables to take into account the weights of the rational Lagrange control points was constructed, which enabled to perform non-invasive NURBS-based IGA. The performance of the implementation was rst illustrated through a series of numerical experiments involving 2D and 3D elasticity problems and the use of the industrial FE software Code_Aster developed by the EDF R&D company [Code_Aster 2014]. According to the authors knowledge, this is the rst time that IGA has been implemented in Code_Aster. Then, the potential of the method was assessed in the nonlinear regime. None of the routines concerning the nonlinearity has been touched: at each iteration of the nonlinear solver, the FE tangent system and the FE nonlinear residual were simply projected on an ad-hoc reduced basis to obtain their isogeometric counterparts.

From a general point of view, this work further contributes to bridging the gap between IGA and standard FEM. Indeed, the developed full algebraic bridge going from Lagrange polynomials to (possibly rational) isogeometric functions oers a new lighting on the relation between IGA and FEA. IGA can be interpreted as the projection of FEA on a specic regular reduced basis. Pushing this interpretation a step further, the implementation procedure may also appear of interest in the context of reduced-order modeling [Kerfriden et al. 2012, Quarteroni et al. 2016, Chinesta et al. 2017]: it could be used to implement any reduced order model in a non-invasive way from a standard FE code.

Figure 1 :

 1 Figure 1: From B-Spline to Lagrange: (a) the initial B-Spline based discretization (four control points associated with four global quadratic B-Spline functions) ; (b) the Bernstein discretization (one control point is added throughout the Bézier decomposition and each element has three local Bernstein shape functions) and (c) Lagrange (FE) discretization (same number of functions as Bernstein, the control points are now on the curve). (top) the curve with the control points and (bottom) the univariate shape functions plots realized in the parameter space (ξ ∈ [0; 1]).

Figure 2 :

 2 Figure 2: Dierent representations of a two-element quadratic curved beam analysed in the example section. (a) the NURBS representation using rational functions and weights generates the exact circular geometry; (b) B-Spline, (c) Bernstein and (d) Lagrange control meshes and corresponding geometry which approximate the circles arcs as described in Section 3.2.2.

  (a) Local least-square approximation for one element. (b) Zoomed window. (c) Geometry error in terms of the relative dierence in radius.

Figure 6 :

 6 Figure 6: The local least-square approximation in case of a quarter circular beam. It consists of an element-level projection technology which nds the best position of the Lagrange nodes that generates a geometry close to the Rational Lagrange one.

Figure 7 :

 7 Figure 7: Problem description of the circular beam under end shear.

  (a) Initial NURBS denition of the horseshoe geometry (b) Associated constructed FE mesh.

Figure 9 :

 9 Figure 9: Example of coarse meshes of the horseshoe.

  (a) Stress σzz over the global surface (b) Detail on the localized stress σyy due to the asymmetric loading.

Figure 10 :

 10 Figure 10: Contour plot of the stress for the solid horseshoe (simple projected IGA).

Figure 12 :

 12 Figure 12: Problem description of the elastoplastic dog-bone sample in tension.

Figure 13 :

 13 Figure 13: Initial IGA mesh (left) and rened automatically generated FE mesh (right).

Figure 14 :

 14 Figure 14: Accuracy of the proposed non-invasive IGA implementation: Longitudinal displacement eld (left, amplication 100) ; Von Mises stress eld (center); and relative discrepancy in Von Mises stress from the FE reference solution computed using Code_Aster (right). The seventh increment, which corresponds to p = 70MPa, is considered for the plot.

Figure 15 :

 15 Figure 15: Force versus displacement curve. Comparison between Code_Aster FE reference (solid black line) and non-invasive IGA (blue crosses).

Acknowledgements

The authors would like to thank Paul Oumaziz, from Núcleo Cientíco Multidisciplinario-DI, Universidad de Talca, Chile, for his valuable help regarding the implementation of the method in Code_Aster.


Armelle Chabot
email: armelle.chabot@ifsttar.fr

M4-5n

M4-5n development for the structural study of Charging Box Inclusion in Electrified Roads

Keywords: Interface, multilayered structures, discontinuities, eRoad 1

may come

M4-5n development for the structural study of Charging Box Inclusion in Electrified Roads

Introduction

Under the objectives of the European Energy Agency to reduce the CO2 emissions for 2050, electrified Roads (eRoad) are proposed. Among all the proposed idea, one of them is to inductively charge the moving Electric Vehicle (EV). It leads to introduce Charging Units (CU) in one of the top layers of the existing multi-layered structures [START_REF] Rim | The development and deployement of online electric vehicules (OLEV)[END_REF] [START_REF] Chen | Electrification of roads: Opportunities and challenges[END_REF] [START_REF] Pérez | Implementing recharging inductive technology on heavy-duty pavement bringing unlimited autonomy to electrical vehicles[END_REF]). An additional thin layer is then placed to covert the all and to assure a good wheel/pavement contact (Figure 1). But due to severe different moving loads (depending on number of truck axles, speed, lateral wandering effect, etc.) and environmental conditions (temperature, water, freezing, etc), failures between the different materials involved may occur quickly (Petit and Chabot, 2017) [START_REF] Buttlar | State-of-the-Art Report of the RILEM Technical Committee 241-MCD Series[END_REF]. Depending on the ratio of modulus between the different materials, cracking and debonding phenomenon between the CU and the structure is one of the major distress parameters to consider for the durability of such new road concepts (Chabot and Deep, 2018). To first analyse the discontinuities in such composite structure, this paper presents a new development done in the dedicated 2D pavement tool, called M4-5nW [START_REF] Nasser | A half-analytical elastic solution for 2D analysis of cracked pavements[END_REF]. It aims to introduce the change of material (due to the CU location) in one (or several) layer(s) of such eRoad multi-layers. Then, a 2D electrified road application case with or not a vertical crack is illustrated.

M4-5n development for an eRoad with a built-in CU

The M4-5nW tool is based on the M4-5n for the Multi-Particle Model of Multilayer Materials (M4) with 5n of equilibrium equations (n: total number of layers) and Winkler's Springs (W) for the half space soil. The M4-5n is specially designed to analyse the delamination of composite materials for multilayer bending problems [START_REF] Chabot | Modèles Multiparticulaires des Matériaux Multicouches M4_5n et M4_(2n+1)M pour l'étude des effets de bord[END_REF]. The M4-5n can be viewed as a superimposition, following the vertical coordinate z, of n Reissner plates linked by shear stresses,

, and the normal stress, , at the interface between the consecutive layers i and i+1

(similarly i-1 and i, i ϵ{1,n-1}, α ϵ {1,2}) where ℎ " # and ℎ " $ are the coordinate of the higher and lower face of the layer i respectively (1).

𝜏 &

","#( (𝑥, 𝑦) = 𝜎 &0 (𝑥, 𝑦, ℎ " # ) = 𝜎 &0 (𝑥, 𝑦, ℎ "#( $ ) 𝜈 2,2#( (𝑥, 𝑦) = 𝜎 00 (𝑥, 𝑦, ℎ " # ) = 𝜎 &0 (𝑥, 𝑦, ℎ "#( $ )

This modelling reduces the real 3D problem to the determination of plane fields (x, y) per each layer i and at the interface i, i+1, (and i-1, i). Thus, the real 3D (2D) object is transforming into a 2D (1D) geometry. It allows convenient parametric calculations without singularity problems (Chabot and Ehrlacher, 1998) [START_REF] Chabot | Mechanical analysis of a mixed mode debonding test for composite pavements[END_REF]. Each i layer has its own behavioural laws, equilibrium equations and lateral boundary conditions. A vertical crack is easy to introduce into each layer separately [START_REF] Chabot | A simplified modeling for cracked pavements / Modèle simplifié de chaussées fissurées[END_REF].

In the so-called M4-5nW tool developed for 2D plane strain problems, the pavement is chosen equivalent to 3 material layers resting on a soil. The soil is assumed equivalent to a combination of a fictitious layer (shear soil layer) ensuring the transfer of shear stresses between the pavement multilayer and Winkler's Springs (W). In that case, a series of M4-5nW equation manipulations lead to the writing of a second-order differential system of twelve analytical equations depending only on one variable, "x" (Nasser and Chabot, 2018) (2). . 𝐴𝑋 55 (𝑥) + 𝐵𝑋 5 (𝑥) + 𝐶𝑋(𝑥) = 𝐷𝑌 ;,( < (𝑥) + 𝐸𝑌 >,? < (𝑥) + 𝐹𝑌 ;,( (𝑥) + 𝐺𝑌 >,? (𝑥)

τ α i,i+1 (x, y) ν i,i+1 (x, y) The matrices [A]12x12, [B]12x12, [C]12x12, [D]12x2, [E]12x2, [F] (2)
12x2 and [G]12x2 depend solely on the geometric and mechanical parameters of the equivalent elastic problem. The first order M4-5nW unknowns 𝑈 ( 2 (𝑥), 𝑈 0 2 (𝑥) and 𝜙 " (𝑥) are, respectively, the average in-plane displacements, average outof-plane displacements and average rotations of the layer i (𝑖 ∈ [START_REF] Buttlar | State-of-the-Art Report of the RILEM Technical Committee 241-MCD Series[END_REF][START_REF] Chabot | Mechanical analysis of a mixed mode debonding test for composite pavements[END_REF]) written as follows (3):

[𝑿] = K 𝑋 ( ⋮ 𝑋 M N 0M×( ; [𝑋] = Q
Following several previous works [START_REF] Chabot | A simplified modeling for cracked pavements / Modèle simplifié de chaussées fissurées[END_REF], the 2D tool solves semi-analytically the equations using the Newmark discretization (4).

W XYZ < #W X < U = W XYZ #W X [ XYZ ; 𝑘 ]#( = 𝑥(𝑗 + 1) -𝑥(𝑗) (4)
To test if the M4-5nW could be a convenient tool for the parametrical study of material inclusions in eRoads with possible macro cracks between materials, it is here proposed to analyse the elastic stress distribution around the charging box for the 2D plane strain example of Chen's PhD works [START_REF] Chen | Electrification of roads: Opportunities and challenges[END_REF][START_REF] Chen | Towards an understanding of the structural performance of future electrified roads: a finite element simulation study[END_REF]. Due to the symmetry of the problem, only the half cross geometry is studied (Figure 2). A "b" index notation, for the points of layers 2 and 3 containing the box inclusion, is used. In Zone 3 of Figure 2, the moving load q(x) describes the vertical and uniform load distribution from the tyre contact. The boundary condition at the upper interface (0,1 index) are then given in (5).

`𝜏( ;,( (𝑥) = 0 𝜈 ;,( (𝑥) = q(𝑥)

for x ∈ [𝑙 h + 𝑙 i , 𝑙 h + 𝑙 i + 𝑙 j ] (5)
Considering 𝐸 k , as the Young's modulus of the soil, the boundary condition at the bottom of the multilayer, interface 4,5 index, provides the connection between the shear layer and the Winkler foundation (with f, a correction factor, equal to 0.83 for a multilayer) (Nasser and Chabot, 2018) [START_REF] Chen | Electrification of roads: Opportunities and challenges[END_REF].

`𝜏( >,? (𝑥) = 0

𝜈 >,? (𝑥) = -𝑘𝑈 0 > (𝑥) 𝑤𝑖𝑡ℎ 𝑘 = 𝐸 k ℎ * 𝑎𝑛𝑑 ℎ * = r 𝑓𝑒 " u 𝐸 " 𝐸 k v "w> "w( (6)
According to [START_REF] Chabot | A simplified modeling for cracked pavements / Modèle simplifié de chaussées fissurées[END_REF], the M4-5n boundary conditions at x=0 and x=l describing the fixed boundary on the both ends (Figure 2), are provided through equation ( 7):

⎩ ⎪ ⎨ ⎪ ⎧ 𝑈 ( " (
0) = 𝑈 ( " (𝑙) = 0 ; 𝜙 ( " (0) = 𝜙 ( " (l) = 0 ⎩ ⎨ ⎧ 𝜙 ( " 0 + 𝑈 0,( " (0) + 1 + 𝜗 " 5𝐸 " •𝜏 ( "$(," (0) + 𝜏 ( ","#( (0)€ = 0 𝜙 ( " (𝑙) + 𝑈 0,( " (𝑙) + 1 + 𝜗 " 5𝐸 " •𝜏 ( "$(," (𝑙) + 𝜏 ( ","#( (𝑙)€ = 0

; 𝑖 ∈ [1,4] (7)
The Zone I, introduced in the M4-5nW tool, pertains to the CU location in layers nos. 2 and 3 of the structure (Figure 2). Since the CU box is made of a different material than the layer materials, the matrices of both Zone 1 (𝑥 ∈ [0, 𝑙 h ]) and outside (𝑥 ∈ [𝑙 h , 𝑙])) differ (with the "b" index notation) in the second-order differential system function of the variable "x" (2). The continuity equations for displacement and rotations at (𝑥 = 𝑙 h ) are then used according to equations ( 8), ( 9) and [START_REF] Nasser | A half-analytical elastic solution for 2D analysis of cracked pavements[END_REF].

𝑈 ( U• (𝑙 h ) = 𝑈 ( U (𝑙 h ) and 𝑈 ( 0• (𝑙 h ) = 𝑈 ( 0 (𝑙 h ) (8) 𝑈 0 U• (𝑙 h ) = 𝑈 0 U (𝑙 h ) and 𝑈 0 0• (𝑙 h ) = 𝑈 0 0 (𝑙 h ) (9) 𝜙 ( U• (𝑙 h ) = 𝜙 ( U (𝑙 h ) and 𝜙 ( 0• (𝑙 h ) = 𝜙 ( 0 (𝑙 h ) (10)
The two interfacial behaviour relationships at interface i, i+1, as given in Equations ( 11) and ( 12) [START_REF] Chabot | Modèles Multiparticulaires des Matériaux Multicouches M4_5n et M4_(2n+1)M pour l'étude des effets de bord[END_REF], are then used to derive the continuity relationship between the two interfacial unknowns present due to the material discontinuity of the CU in layers 2 and 3.

𝑒 "#( 12 𝑈 0,( "#( (𝑥) + 𝑒 " 12 𝑈 0,( " (𝑥) + 𝑈 ( "#( (𝑥) -𝑈 ( " (𝑥) -5𝑒 "#( 12 𝜙 ( "#( (𝑥) -5𝑒 " 12 𝜙 ( " (𝑥)

= -(1 + 𝑣 " )𝑒 " 12𝐸 " 𝜏 ( "$(," (𝑥) -(1 + 𝑣 "#( )𝑒 "#( 12𝐸 "#( 𝜏 ( "#(,"#U (𝑥) + ƒ (1 + 𝑣 " )𝑒 " 4𝐸 " + (1 + 𝑣 "#( )𝑒 "#( 4𝐸 "#( " 𝜏 ( ","#( (𝑥) [START_REF] Pérez | Implementing recharging inductive technology on heavy-duty pavement bringing unlimited autonomy to electrical vehicles[END_REF] 𝑈 0 "#( (𝑥) -𝑈 0 " (𝑥) = 9 70 𝑒 " 𝐸 " 𝜈 "$(," (𝑥) + 9 70

𝑒 "#( 𝐸 "#( 𝜈 "#(,"#U (𝑥) + 13 35 ƒ 𝑒 " 𝐸 " + 𝑒 "#( 𝐸 "#( " 𝜈 ","#( (𝑥) [START_REF] Rim | The development and deployement of online electric vehicules (OLEV)[END_REF] In subtracting Equation ( 12) at point x = lb and using Equation ( 9), Equations ( 13) to (15) are obtained for the three interfaces affected by the presence of the CU (between layers 1 and 2, then between layers 2 and 3 and lastly between layers 3 and 4) (Figure 2).

13 35 ƒ 𝑒 ( 𝐸 ( + 𝑒 U 𝐸 hU " 𝜈 h((,U) (𝑙 h ) + 9 70 𝑒 U 𝐸 •U 𝜈 h(U,0) (𝑙 h ) = 13 35 ƒ 𝑒 ( 𝐸 ( + 𝑒 U 𝐸 U " 𝜈 (,U (𝑙 h ) + 9 70 𝑒 U 𝐸 U 𝜈 U,0 (𝑙 h ) (13)
9 70

𝑒 U 𝐸 hU 𝜈 h((,U) (𝑙 h ) + 13 35 ƒ 𝑒 U 𝐸 hU + 𝑒 0 𝐸 h0 " 𝜈 h(U,0) (𝑙 h ) + 9 70 𝑒 0 𝐸 •0 𝜈 h(0,>) (𝑙 h ) = 9 70 𝑒 U 𝐸 U 𝜈 (,U (𝑙 h ) + 13 35 ƒ 𝑒 U 𝐸 U + 𝑒 0 𝐸 0 " 𝜈 U,0 (𝑙 h ) + 9 70 𝑒 0 𝐸 0 𝜈 0,> (𝑥𝑙 h ) (14)
13 35 ƒ 𝑒 0 𝐸 h0 + 𝑒 > 𝐸 > " 𝜈 h(0,>) (𝑙 h ) + 9 70 𝑒 0 𝐸 •0 𝜈 h(U,0) (𝑙 h ) = 13 35 ƒ 𝑒 0 𝐸 0 + 𝑒 > 𝐸 > " 𝜈 0,> (𝑙 h ) + 9 70 𝑒 0 𝐸 0 𝜈 U,0 (𝑙 h ) (15)
Similarly, by subtracting Equation [START_REF] Pérez | Implementing recharging inductive technology on heavy-duty pavement bringing unlimited autonomy to electrical vehicles[END_REF], at x = lb and using Equations ( 8) to (10), a system of the three additional equations is obtained to assure the continuity of the interface shear stresses of the three interfaces involved. The system can thus be solved numerically given that all the boundary, continuity and interface conditions have been formulated and defined for this eRoad inclusion problem.

2D case results

To illustrate the results, the interfacial normal stresses are drawn for several elastic pavement structure cases (Figure 3). The load characteristics as well as the box dimension and layer thicknesses are given on Figure 1. The solutions are obtained using a uniform discretization (mesh width of 20mm) of the variable "x" and approximations around the changes of Zone 1 and 2 (Chabot and Deep, 2018). The soil model has been adapted with E s =160MPa ( 6) and (E 4 =270MPa, 𝜐 > = 0.35) in order to fit Chen's result [START_REF] Chen | Towards an understanding of the structural performance of future electrified roads: a finite element simulation study[END_REF]. These equivalent soil values are kept for all the results presented here after. The prefabricated CU is supposed been made with a cement concrete material (E b2 =E b3 =30000MPa; 𝜐 •U = 𝜐 •0 = 0.25) (Figure 2). For each case, a flexible pavement (purple curves with 𝜐 ( = 𝜐 U = 𝜐 0 = 0.35 𝑎𝑛𝑑 𝐸 ( = 𝐸 U = 𝐸 0 ) with the assumption of the presence of a vertical crack between the CU and the layer is compared to a perfect bonded composite pavement made of a cement concrete for layers 1 and 2 with the same characteristic that the CU inclusion (𝐸 ( = 𝐸 U = 30000MPa, 𝜐 ( = 𝜐 U = 0.25 ) casted on an asphalt layer with 𝐸 0 ≠ 𝑎𝑛𝑑 𝜐 0 = 0.35. to -0.2, see Figure1), the higher the intensity of the interfacial normal stresses are between layers 1 and 2. In addition, the highest the asphalt modulus of the flexible pavement is (E BBSG (15°C)= 5400MPa, E GB3 (15°C)= 9300MPa), the highest the interfacial normal stress, between layers 1 and 2, is. Although it is reducing the ratio of modulus between layers 2 and 3 in case of the composite pavement structure.

Conclusion

In order to adapt pavement structures into eRoads with CUs that inductively charge the moving Electric Vehicle, pavement durable solutions in terms of material and structure need to be found. Due to change of materials between the inductive box inclusion and the layers, debonding phenomenon between those may occur quickly and under such condition, it is one of the major distress parameters to study for pavement structures. In that aim a new development that introduces material discontinuity in one or several layers have been done in an alternative modelling tool, called the M4-5nW. To contribute to this energy transition period, the tool offers quick parametric mechanical study possibilities.

FIG. 1 -

 1 FIG. 1 -A half cross-sectional geometry of an eRoad with an inductive CU (After Chen et al., 2019)

  𝑌 ;,( ] = S 𝜏 ;,( 𝜈 ;,( T U×( ; [𝑌 M,M#( ] = S 𝜏 M,M#( 𝑣 M,M#( T U×(𝑖 ∈[START_REF] Buttlar | State-of-the-Art Report of the RILEM Technical Committee 241-MCD Series[END_REF][START_REF] Chabot | Mechanical analysis of a mixed mode debonding test for composite pavements[END_REF]

FIG. 2 -

 2 FIG. 2 -M4-5nW equivalent eRoad model of half the cross-sectional geometry (After Chabot and Deep, 2018)

FIG. 3 -

 3 FIG. 3 -Interfacial normal stress comparison between the various cases: (a) BBSG (40°,𝐸 0 = 1000𝑀𝑃𝑎, w=0); (b) BBSG (40°, w=-0.2); (c) BBSG (15°C, w=-0.2) ; GB3 (15°C, w=-0.2)


Lucas Deniel
email: lucas.deniel@u-psud.fr

Mathilde Gay

Diego Perez-Galacho

Charles Baudot

Laurent Bramerie

Oskars Ozolins

Frédéric Boeuf

Laurent Viven

Christophe Peucheret

Delphine Marris-Morini

Diego Pérez-Galacho

Fréderic Boeuf

O-band DAC-less PAM-4 generation with a silicon dual-drive Mach-Zehnder Modulator

Keywords: Silicon, Mach-Zehnder modulator, PAM-4, O-band, DAC-less

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

INTRODUCTION

For the next generation of intra-datacenter and short-reach data communications, 4-level pulse amplitude modulation (PAM-4) has been chosen as a convenient format to fulfill the increasing spectral efficiency requirements. Silicon photonics is a promising technology for these applications as it provides economical, flexible, and low-power solutions, in addition to its mature complementary metal-oxide semiconductor (CMOS) fabrication process and the possibility of co-integration with electronics [START_REF]GbE Task Force[END_REF][START_REF] Shi | Silicon photonic modulators for PAM transmissions[END_REF].

Usually, PAM-4 signal generation is done in the electrical domain, and this signal is used to drive the electrooptical modulator. In this case, power hungry digital-to-analog converters (DAC) are needed to convert 2-level non-return-to-zero (NRZ) data-streams into 4-level ones. For this reason, the generation of PAM-4 signals in the optical domain was proposed and successfully demonstrated using various structures of silicon modulators on which two different NRZ binary sequences were directly applied. These structures are sometimes referred to as optical DACs, and most of these demonstrations were done in the C-band [START_REF] Roshan-Zamir | A 40Gb/s PAM4 optical DAC silicon microring resonator modulator transmitter[END_REF][START_REF] Patel | Silicon photonic segmented modulator-based electro-optic DAC for 100Gb/s PAM-4 Generation[END_REF][START_REF] Li | Silicon photonic ring-assisted MZI for 50Gb/s DAC-less and DSP-free PAM-4 transmission[END_REF]. Interestingly, the well-known dual-drive Mach-Zehnder modulator (MZM) structure is used for classical OOK modulation, but it was demonstrated that it could also be used for quadrature phase shift keying (QPSK) and PAM-4 generation, making it an interesting choice for the evolving transmission standards. [START_REF] Zheng | Silicon PAM-4 optical modulator driven by two binary electrical signals with different peak-to-peak voltages[END_REF][START_REF] Pérez-Galacho | QPSK modulation in the O-band using a single dual-drive Mach-Zehnder silicon modulator[END_REF] We demonstrate optical PAM-4 generation at 20 Gbit/s using an O-band silicon dual-drive MZM. In the first part, we introduce the MZM design and the phase shifter characteristics. A numerical model is then used to find out the peak-to-peak voltages to apply on the MZM arms to obtain 4 equally spaced intensity levels while benefiting from the transmission full dynamic range. Experimental results are then reported, including bit error rate measurements.

MODULATOR DESIGN AND PHASE SHIFTER CHARACTERISTICS

The device fabrication was carried out by ST-Microelectronics on a 300-mm silicon-on-insulator (SOI) technological platform [START_REF] Baudot | Daphne silicon photonics technological platform for research and development on WDM applications[END_REF]. The optical access of the chip are ensured by grating couplers. The combiner and splitter of the MZM structure are 3-dB 12 multimode interference (MMI) based couplers. Each arm contains a phase shifter which consists in a 1-mm-long lateral PN junction embedded in a rib waveguide, operating in depletion regime. The waveguide height is 300 nm and its width is 400 nm, allowing propagation of the fundamental TE mode, while the slab part is 50-nm thick. The targeted doping concentrations are 510 17 cm -3 for the P region and 4.510 17 cm -3 for the N region. Higher doping concentrations of about 10 19 cm -3 are used near the electrical contacts to reduce the access resistance, thus enhancing the modulator bandwidth. Traveling wave coplanar electrodes on the surface of the chip are used for radio-frequency (RF) waves propagation and metallic strips are used as resistors (heaters) close to each arm to thermally tune the static operating point of the MZM, by applying a DC voltage.

Static characterizations of the phase shifter efficiency were made experimentally by tracking the resonance wavelengths of a ring resonator, as it allows a precise measurement of the phase variation within the waveguide as a function of the applied voltage. The measured phase shifter efficiency varies from 1.1 V cm at -0.5 V to 1.4 V cm at -7 V. A 17-GHz electro-optical bandwidth was measured on the MZM using an Agilent 86030A 50 GHz Lightwave Component Analyzer.

The phase shifter insertion loss was estimated as a function of the voltage through a numerical model. First, carrier concentrations within the PN junction as a function of the applied voltage are obtained using an electrical simulation. Soref model is then used to calculate the corresponding local absorption coefficients, and a mode solver is finally used to estimate the total effective attenuation factor of the waveguide. Loss estimation varies from 1.8 dB mm -1 at 0V to 1.1 dB mm -1 at -7V.

PAM-4 MZM MODEL AND COMPARISON WITH EXPERIMENTAL RESULTS

As the effective index change is not linear with the applied voltage, as the insertion loss depends on the applied voltage, finding the dynamic voltage amplitudes to apply on each arm to obtain 4 equally-spaced intensity levels while maximizing the extinction ratio in the optical domain is not straightforward. For this reason, the following equations are numerically implemented to model the MZM:

( ) ( ) ( ) IN 2 exp ( ) exp j ( ) 2 n eff eff n eff eff n E E V L n n V L      = - +  +      (1) 2 2 1 2 OUT 2 H H j j E e E e E   - + = (2)
where EIN is the input electric field, En is the electric field at the output of arm labelled n, EOUT is the output electric field, L is the phase shifter length,  is the wavelength, H is the static phase shift introduced by thermal heaters, eff and neff are the effective attenuation factor and refractive index of the optical mode, respectively, eff(Vn) and

BER MEASUREMENTS

A BER measurement is conducted to evaluate the MZM PAM-4 performance under this dual-driving method. For this experiment, light is produced by an external cavity laser at 1310 nm and input to the device after being amplified with a semiconductor optical amplifier (SOA). At the output of the chip, the modulated beam is amplified with a praseodymium doped fiber O-band amplifier and split in two paths with a fibered directional coupler. One half of the power is sent to a 25-GHz oscilloscope displaying the eye-diagram in real time. The other half is sent to a 40-GHz photodiode and the converted electrical signal goes through a 35-GHz amplifier before being received by the error counter of a BER tester. Two complementary pseudo random binary sequences (PRBSs) are generated electrically at a 10-Gbit/s rate. One of them is delayed by 200 ps with an RF tunable delay line, and both signals are amplified with tunable gain RF amplifiers. A -4-V DC voltage is added on both lines with DC bias tees to bias the PN junction in reverse. The two driving signals are then injected into the RF electrodes of the MZM via RF probes before being absorbed by 50  terminations loaded after DC blocks on the electrical output. The peak-topeak voltages are set to 4 V on one arm and 7 V on the other one by tuning the RF amplifier gains. The optical power received on the photodiode is swept from -8.5 to 0.5 dBm for the BER measurement. Fig. 2 shows the obtained BER values as a function of the received power. It can be seen that the BER is under the 3.810 -3 hard-decision forward-error-correction (HD-FEC) threshold for a received power of -6 dBm and under the 2.210 -4 KP4 threshold for a received power of -3.4dBm.

CONCLUSIONS

A 20-Gbit/s PAM-4 signal was generated at 1310 nm using a silicon dual drive MZM driven by two different NRZ binary signals. After an introduction of this modulator design and characteristics, a numerical model was presented and used to find out the optimum voltages to apply to each arm of the MZM for PAM-4 generation. Finally, a BER under the HD-FEC threshold of 3.810 -3 was measured for a received power of -6 dBm using the silicon MZM as an optical PAM-4 generator.

  neff(Vn) are the effective attenuation factor change and refractive index change of the optical mode, respectively, under the application of the voltage Vn on the phase shifter labelled n. eff + eff(Vn) is estimated and neff(Vn) is calculated as described in the previous section. Fig. 1(a) shows the PAM-4 levels obtained with the MZM numerical model by applying peak-to-peak voltages of 4 V on one arm of the MZM and 7 V on the other. Fig. 1(b) shows an eye diagram obtained experimentally by applying the same voltages to the real silicon MZM, biased at the quadrature point. The matching of level spacings between the two figures demonstrates the numerical model accuracy.

Figure 1 .

 1 Figure 1. (a) Obtained PAM-4 levels by applying a 4-V peak-to-peak voltage on one arm and 7-V peak-to-peak voltage on the second arm on the silicon MZM numerical model. (b) Experimental eye diagram obtained after applying the same peak-to-peak voltages on the silicon MZM.

Figure 2 .

 2 Figure 2. Measured O-band PAM-4 BER with the silicon MZM driven by a 7-VPP and a 4-VPP PRBS as a function of the received optical power.

ACKNOWLEDGEMENTS

This work was supported by the European project Cosmicc (H2020-ICT-27-2015-688516).


Doping profile extraction in thin SOI films: application to A2RAM

Keywords: A2RAM, C-V characteristic, SOI, electrical characterization, doping profile. I

We propose for the first time a method based on C-V measurement to extract the bridge doping profile which governs the A2RAM performances. Assessed with TCAD simulation and simple extraction model adapted from bulk devices, this technique is validated with experimental data.

INTRODUCTION

The A2RAM (Fig. 1) is a 1T-DRAM [START_REF] Okhonin | Zero Capacitor Embedded Memory Technology for System on Chip[END_REF] with an N-doped layer (nMOS case) called 'the bridge' located at the bottom of the low-doped p-type body. The bridge allows short-circuiting the source and drain. The concept of A2RAM has been presented in [START_REF] Rodriguez | New concepts for 1T-DRAMS: covercoming the scaling limits[END_REF] and its experimental performance in [START_REF] Rodriguez | Experimental demonstration of capacitorless A2RAM cells on silicon-on-insulator[END_REF]. Recently, the scalability issues of the A2RAM have been studied by TCAD simulations [START_REF] Wakam | Optimization guidelines of A2RAM cell performance through TCAD simulations[END_REF]. We noticed some inconsistencies with measurements, presumably related to the mismatch of the bridge doping profile between TCAD simulations and experiments. As the A2RAM behavior is strongly related to these parameters, we need an accurate solution to evaluate them. The use of physical characterization such as Secondary-Ion Mass Spectroscopy (SIMS) is not reliable here because (i) the silicon film is too thin and (ii) the target doping is too low (around 10 18 cm -3 ). Thus we need a method based on electrical measurement. For this purpose, the C-V electrical characterization is emulated through TCAD simulations and validated with experiments. II. A2RAM: TCAD SIMULATION METHODOLOGY AND C-V CHARACTERISTIC We simulate the A2RAM structure (Fig. 1) as fabricated in [START_REF] Rodriguez | Experimental demonstration of capacitorless A2RAM cells on silicon-on-insulator[END_REF] with Synopsys TCAD tools [START_REF]TCAD Sentaurus software -Synopsys vM[END_REF] using the same models as in [START_REF] Wakam | Optimization guidelines of A2RAM cell performance through TCAD simulations[END_REF].

The main technological parameters are: abrupt Gaussian doping profiles, body thickness Tbody = 14 nm, bridge thickness Tbridge = 22 nm, bridge doping Nbridge = 10 18 cm -3 , Equivalent Oxide Thickness EOT = 3.1 nm and Buried Oxide thickness (BOX) TBOX = 10 nm. We performed small signal simulation and report on Fig. 2

III. BRIDGE DOPING PROFILE EXTRACTION METHOLOGY

WITH TCAD To extract the bridge doping profile, we revisit the method described in [START_REF] Schroeder | Semiconductor material and device characterization[END_REF][START_REF] Ghibaudol | Improved method for the oxide thickness extraction in MOS structures with ultra-thin gate dielectrics[END_REF][START_REF] Sareen | Effect of Si cap layer on parasitic channel operation in Si/SiGe metal-oxidesemiconductor structures[END_REF][START_REF] Chattopadhyay | C-V characterization of strained Si/SiGe multipleheterojunction capacitors as a tool for heterojunction MOSFET channel design[END_REF][START_REF] Lee | Non-homogeneous electrical transport through silicon-on-saphire thin films: evidence of the intermal stress influence[END_REF][11]. The depletion depth depends on the gate voltage: a small signal variation on the metal modifies the depletion depth and the charge density at bottom of the film.

A. Extraction of the 'doping profile of the bridge'

Having Cgg, and considering the depletion regime in the body and all dopants ionized, we need to compute the so-called 'doping function' Nnet (1) defined as an apparent density of charge in the silicon film [START_REF] Schroeder | Semiconductor material and device characterization[END_REF]: with q the elementary charge, Si the silicon permittivity and A the gate area. Fig. 3-(a) shows │Nnet│ versus Vg derived from the C-V characteristic of the A2RAM (Fig. 2-(a)). In the gate bias range between Vthb and Vthf, the 'doping function' presents a broad peak with the maximum value 9.10 17 cm -3 close to the doping value set in TCAD (Fig. 1).

N 2 qA ε (1)
In Fig. 4-(a), similar simulations are performed for various bridge doping value Nbridge (3.10 18 , 10 18 , 3.10 17 , and 10 17 cm -3 ) keeping constant the bridge thickness Tbridge (22 nm). With higher values of Nbridge (≥ 3.10 18 cm -3 ), the extracted bridge doping can be accurately evaluated (Fig. 4-a) because less and less portion of the bridge is depleted by the field effect. If the Nbridge is low (< 3.10 17 cm -3 ), the extracted value has no sense because the bridge is fully depleted. This observation is confirmed in Fig. 6-a where we compare the C-V of A2RAM as described in Fig. 1 but with Nbridge = 10 17 cm -3 , and 3.10 17 cm -3 , with the C-V of a SOI device transistor with the same silicon thickness. For a bridge doping higher than 3.10 17 cm -3 , the change in slope starts to become significant, but as the device is still fully depleted (including the bridge), that is why the extraction remains "inaccurate". Indeed, for Nbridge < 3.10 17 cm - 3 the shape of the C-V curve shows the same slope even if we have an Nbridge.

On the other hand, Fig. 5-(a) shows the impact of the bridge thickness Tbridge (28, 22, 14, 10, and 8 nm) while keeping constant the bridge doping (Nbridge =10 18 cm -3 ). As for the bridge doping, good resolution of the Nbridge extracted is guaranteed for thicker bridge (Tbridge ≥ 22 nm). But, the extracted Nbridge is less accurate for films thinner than 14 nm because the bridge is more and more depleted. For ultrathin Tbidge (< 8 nm), the bridge is totally depleted. In Fig. 6-b, we compare the C-V curves of an ARAM as described in Fig. 1 but with Tbridge = 8 nm and a SOI device transistor with same silicon thickness. They are similar with a smooth slope in the case of the A2RAM.

B. Extraction of the 'Body' thickness

We aim to extract the bridge doping Nbridge and the body thickness Tbody. At a gate voltage between Vthb and Vthf, the bridge is still formed, and the body is depleted: this means that the gate capacitance Cgg corresponds to Cox and Csi in series, with Cox and Csi the gate oxide and the depletion layer capacitance respectively. In this particular case, the depletion layer corresponds to the p-doped body layer, so we have:

depletion depth T ε C C (2)
As we precisely know the value of the total silicon thickness TSi, we can easily deduce Tbridge (Tbridge = TSi -Tbody). ). The value extracted, Tbody = 21.6 nm, is 8 nm higher than the 14 nm of Tbody defined in TCAD (Fig. 1). This difference is related to the relatively low value of the bridge doping. To confirm this hypothesis, we have extracted Tbody for different Nbridge, and the results are shown in Fig. 4-(b). As in Nbridge extraction, Tbody extraction is more accurate for higher doping: for Nbridge = 3.10 18 cm -3 , we obtain Tbody= 15.7 nm, just 2 nm thicker than the expected 14 nm. If Nbridge is too low (<3.10 17 cm -3 ), the silicon film is fully depleted (Fig. 6-a), including the bridge. This implies that Cgg capacitance depends on BOX capacitance and equation ( 2) cannot be used (this explains why Tbody>Tsi for low Nbridge). In Fig. 5-(b), the results achieved for different bridge thickness Tbridge values (28, 22, 14, 10, and 8 nm) at the same bridge doping Nbridge (10 18 cm -3 ) are shown. For Tbridge > 8 nm, the Tbody extracted are 16, 21, 27, and 32 nm instead of 8, 14, 22, and 26 nm thus inaccurate because the initial bridge doping is relative low as we have noticed in Fig. 3 and Fig. 4. For low Tbrige (< 8 nm) like for low values of Nbridge the value of Tbody extracted is meaningless because the bridge is depleted (Fig. 6-(b)).

C. Comments about the resolution of the technique

The observations derived from Fig. 4 could infer that the minimum value of the bridge doping that can be extracted is 3.10 17 cm -3 . However, Fig. 7 shows if the back gate is properly biased (Vback), it is possible to extract a doping value even for low doping (< 3.10 17 cm -3 ). Reducing the bridge doping value while keeping the same mid gap metal at the back gate changes the flat band condition, thus we need to apply a back bias Vback. Even in flat band condition at the back interface, if the value of Nbridge is low (< 3.10 17 cm -3 ), the value of Tbody extracted is less accurate. In fact, the depletion of the bridge due to the pn junction (body-bridge) is no longer negligible. Considering the standard expression of the depletion depth at thermal equilibrium [START_REF] Mathieu | Physique des semiconducteurs et des composants électroniques[END_REF], we can estimate the depletion depth in the ntype (3.10 17 cm -3 ) bridge equals 3.3 nm. Since the 'bridge thickness' is 22 nm, neglecting the depletion of the bridge due to the body junction is not realistic.

In the same way, from Fig. 5 one can deduce that below Tbody = 8 nm it is not possible to extract a bridge doping Nbridge and a body thickness Tbody. But in Fig. 8-(a) and 8-(b) biasing the back interface enhances the resolution of the extraction. As observed in Fig. 7 the shift of the flat band condition can be managed by biasing the back gate. However, there is a limit of the bias which can be applied on the back interface Vback. In Fig. 9, we compare the electron density extracted at gate voltage Vg = 0 V with the original 'doping concentration' for different structures at variable Vback. For Vback > 500 mV, the electron density at the back interface is higher than the concentration in the volume of the bridge. Its variation is just few nanometer above the BOX, due to "weak inversion". That is why for this work, we have limited Vback at 350 mV.

Previous simulations consider a metal electrode directly below the BOX with a mid-gap work-function. However, in real A2RAM device [START_REF] Rodriguez | Experimental demonstration of capacitorless A2RAM cells on silicon-on-insulator[END_REF] a silicon ground plane GP acts as the back gate. To determine if the presence of GP could have an impact on our extraction methodology, we performed simulation of A2RAM including a p-type GP (doping 10 18 cm -3 ). Because of the back gate work-function shift, C-V curves are not aligned (Fig. 10-a): we need to apply 500 mV on the GP to superpose the C-V plots for metal and silicon back gate (Fig. 10-(a,b)). In this condition, the bridge doping extracted is the same as the one without the GP (Fig. 10-c). Note also that Tbody extraction is not impacted by the GP (Fig. 10-d). 1 with Nbridge = 3.10 18 cm -3 , (b) A2RAM as defined in Fig. 1 with Nbridge = 10 17 cm -3 , (c) A2RAM as defined in Fig. 1 with Tbridge = 28 nm, and (d) A2RAM as defined in Fig. 1 with Tbridge = 8 nm; obtained TCAD simulations.

IV. EXPERIMENTAL RESULTS AND DISCUSSION

We extract the doping profile of the bridge on the samples fabricated at LETI with process flow described in [START_REF] Rodriguez | Experimental demonstration of capacitorless A2RAM cells on silicon-on-insulator[END_REF]. The experimental C-V characteristics with two Vback shown in Fig. 11-(a) are comparable to C-V from TCAD simulations.

The shift of the Vthb that can be evidenced on Fig. 11-(b) (at Vback = 0 V) like on Fig. 10-(b) is mainly due to the presence of the GP. In consequence, the bridge doping value of 4.10 17 cm -3 read on Fig. 11-(c) at Vback = 0 V is lower than the real value but the accuracy can be improved by biasing the back gate at Vback = 500 mV as evidenced by simulation. In this condition, the extracted doping value equals nearly 10 18 cm -3 . The Tbody extracted (Fig. 11-d) at Vback = 0 V (27 nm) or at Vback = 500 mV (24 nm) is overestimated because the value of Nbridge is low (as we have noticed and concluded in Fig. 3-b, Fig. 4-b).

V. CONCLUSION

In this paper, we have developed an electrical characterization method to extract the doping profile in thin SOI films with an application to the A2RAM. We have noticed that for low values of Nbridge and Tbridge, the accuracy can be enhanced by biasing the back gate. The extraction technique has been then implemented on the samples fabricated in LETI [START_REF] Rodriguez | Experimental demonstration of capacitorless A2RAM cells on silicon-on-insulator[END_REF] and the results are in agreement with the ones shown in simulation.

  -(a) the gate capacitance Cgg versus the gate voltage Vg for both the A2RAM and a standard SOI device transistor. Similarly to a back-biased FD SOI devices with an inverted back channel, the A2RAM C-V characteristic shows a double hump related to the presence of the inversion charge at back and front interfaces. It is possible to define the threshold voltage for the formation of these channels. Fig.2-(b) shows the first derivative of the C-V curves. When the front-gate bias increases, the first hump is related to the 'conduction of the bridge' at Vthb and the second one to the onset of the front channel at Vthf (Fig.2-(b)).

Figure 1 .

 1 Figure 1. Simulated A2RAM cell with parameters defined by the process flow in [3].

Figure 2 .

 2 Figure 2. (a) C(Vg) characteristic of the A2RAM and standard SOI transistor with same silicon thickness but no doped bridge. (b) First derivative of the C-V curves.

Fig. 3 -

 3 Fig.3-(b)shows the depletion depth versus gate bias Vg. The depletion depth corresponding to the electrical thickness of the 'body' must be determined at same gate voltage as the one used for the bridge doping extraction (in this case Vg=0). The value extracted, Tbody = 21.6 nm, is 8 nm higher than the 14 nm of Tbody defined in TCAD (Fig.1). This difference is related to the relatively low value of the bridge doping. To confirm this hypothesis, we have extracted Tbody for different Nbridge, and the results are shown in Fig.4-(b). As in Nbridge extraction, Tbody extraction is more accurate for higher doping: for Nbridge = 3.10 18 cm -3 , we obtain Tbody= 15.7 nm, just 2 nm thicker than the expected 14 nm. If Nbridge is too low (<3.10 17 cm -3 ), the silicon film is fully depleted (Fig.6-a), including the bridge. This implies that Cgg capacitance depends on BOX capacitance and equation (2) cannot be used (this explains why Tbody>Tsi for low Nbridge).

Figure 3 .

 3 Figure 3. (a) Doping function and (b)-depletion depth in A2RAM by TCAD simulations.

Figure 4 .

 4 Figure 4. (a) Doping function and (b) depletion depth for different Nbridge (3.10 18 , 10 18 , 3.10 17 , and 10 17 cm -3 ) obtained by A2RAM TCAD simulations.

Figure 5 .

 5 Figure 5. (a) Doping function and (b) depletion depth for different Tbridge (28, 22, 14, 10 and 8 nm) obtained by A2RAM TCAD simulations.

Figure 6 .

 6 Figure 6. Comparison of the C-V of a SOI device transistor with (a) A2RAM as defined in Fig.1 but with Nbridge = 3.10 17 cm -3 and Nbridge = 10 17 cm -3 , (b) A2RAM as defined in Fig.1 but with Tbridge = 8 nm, obtained by TCAD simulations.

Figure 7 .

 7 Figure 7. (a) Doping function and (b) depletion depth with back-gate bias (Vback = 350 mV) for different Nbridge (3.10 18 , 10 18 , 3.10 17 , and 10 17 cm -3 ) obtained by simulations.

Figure 10 .

 10 Figure 10. (a) C-V characteristic, (b) first derivative of the C-V, (c) doping function, and (d) depletion depth of the A2RAM with mid-gap back metal and with a p-type GP (TCAD simulations).

Figure 9 .

 9 Figure 9. comparison of the doping concentration with the electron density extracted at Vg = 0 V for different Vback for an: (a) A2RAM as defined in Fig.1with Nbridge = 3.10 18 cm -3 , (b) A2RAM as defined in Fig.1with Nbridge = 10 17 cm -3 , (c) A2RAM as defined in Fig.1with Tbridge = 28 nm, and (d) A2RAM as defined in Fig.1with Tbridge = 8 nm; obtained TCAD simulations.

Figure 8 .

 8 Figure 8. (a) Doping function and (b) depletion depth with the back gate (Vback) at 350 mV for different Tbridge (28, 22, 14, 10 and 8 nm) obtained by A2RAM TCAD simulations.

ACKNOWLEDGMENT

The research leading to these results has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 687931 REMINDER. This work has been partially supported by the LabEx Minos ANR-10-LABX-55-01.