A Lebée

K Sab

Homogenization of a space frame as a thick plate: application of the Bending-Gradient theory to a beam lattice

Keywords: Plate theory, Higher-order Models, Homogenization, Periodic plates, Space frame

The Bending-Gradient theory for thick plates is the extension to heterogeneous plates of Reissner-Mindlin theory originally designed for homogeneous plates. In this paper the Bending-Gradient theory is extended to in-plane periodic structures made of connected beams (space frames) which can be considered macroscopically as a plate. Its application to a square beam lattice reveals that classical Reissner-Mindlin theory cannot properly model such microstructures. Comparisons with exact solutions show that only the Bending-Gradient theory captures second order effects in both deflection and local stress fields.

Introduction

The classical theory of plates, known also as Kirchhoff-Love plate theory is based on the assumption that the normal to the mid-plane of the plate remains normal after transformation. This theory is also the first order of the asymptotic expansion with respect to the thickness [START_REF] Ciarlet | Justification Of The 2-Dimensional Linear Plate Model[END_REF]. Thus, it presents a good theoretical justification and was soundly extended to the case of periodic plates [START_REF] Caillerie | Thin elastic and periodic plates[END_REF]. It enables to have a first order estimate of the macroscopic deflection as well as local stress fields. In most applications the first order deflection is accurate enough. However, this theory does not capture the local effect of shear forces on the microstructure because shear forces are one higher-order derivative of the bending moment in equilibrium equations (Q α = M αβ,β ).

Because shear forces are part of the macroscopic equilibrium of the plate, their effect is also of great interest for engineers when designing structures. However, modeling properly the action of shear forces is still a controversial issue. Reissner [START_REF]The effect of transverse shear deformation on the bending of elastic plates[END_REF] suggested a model for homogeneous plates based on a parabolic distribution of transverse shear stress through the thickness (Reissner-Mindlin theory). This model performs well for homogeneous plates and gives more natural boundary conditions than those of Kirchhoff-Love theory. Thus, it is appreciated by engineers and broadly used in applied mechanics. However, the direct extension of this model to more complex microstructures raised many difficulties. Many suggestions were made for laminated plates [START_REF] Whitney | Stress Analysis of Thick Laminated Composite and Sandwich Plates[END_REF][START_REF] Reddy | On refined computational models of composite laminates[END_REF][START_REF] Noor | An assessment of five modeling approaches for thermo-mechanical stress analysis of laminated composite panels[END_REF][START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF][START_REF] Nguyen | Shear Correction Factors for Functionally Graded Plates[END_REF][START_REF] Nguyen | First-order shear deformation plate models for functionally graded materials[END_REF] as well as in-plane periodic plates [START_REF] Lewiski | Effective models of composite periodic platesI. Asymptotic solution[END_REF][START_REF] Cecchi | A homogenized ReissnerMindlin model for orthotropic periodic plates: Application to brickwork panels[END_REF], leading to more complex models.

Revisiting the approach from Reissner [START_REF]The effect of transverse shear deformation on the bending of elastic plates[END_REF] directly with laminated plates, Lebée and Sab [START_REF] Lebée | A Bending-Gradient model for thick plates. Part I: Theory[END_REF][START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF] showed that the transverse shear static variables which come out when the plate is heterogeneous are not shear forces Q α but the full gradient of the bending moment R αβγ = M αβ,γ . Using conventional variational tools, they derived a new plate theory -called Bending-Gradient theory -which is actually turned into Reissner-Mindlin theory when the plate is homogeneous. This new plate theory is seen by the authors as an extension of Reissner's theory to heterogeneous plates which preserves most of its simplicity. Originally designed for laminated plates, it was also extended to in-plane periodic plates using averaging considerations such as Hill-Mandel principle and successfully applied to sandwich panels [START_REF] Lebée | Homogenization of thick periodic plates: Application of the Bending-Gradient plate theory to a folded core sandwich panel[END_REF][START_REF] Lebée | Homogenization of cellular sandwich panels[END_REF].

In order to give a more comprehensive illustration of the features of this new theory, we extend its homogenization scheme to space frames (Section 2). Space frames are large roofings made of many identical unit-cells. Numerous illustration are given in Buckminster Fuller's achievements. In this work, a space frame is a unit-cell made of connected beams periodically reproduced in a plane and which "from far" can be considered as a plate. Many devices fall into this category: space trusses [START_REF] Deshpande | Collapse of truss core sandwich beams in 3-point bending[END_REF][START_REF] Deshpande | Effective properties of the octet-truss lattice material[END_REF][START_REF] Wallach | Mechanical behavior of a three-dimensional truss material[END_REF][START_REF] Wicks | Optimal truss plates[END_REF], tensegrity, nexorade [START_REF] Douthe | Design of nexorades or reciprocal frame systems with the dynamic relaxation method[END_REF], gridshells [START_REF] Douthe | Gridshell in composite materials: Towards wide span shelters[END_REF], lattices, expanded metal, gratings, etc.

Let's point out that very few methods exist in the literature for deriving a thick plate macroscopic model when the microstructure is made of structural elements. Clearly, theories for laminates are not suitable for the present microstructure and we need an approach which takes into account the periodicity of the plate. There is also a large literature dedicated to the homogenization of beam lattices (also referred as discrete medias [START_REF] Pradel | Homogenization of discrete media[END_REF][START_REF] Florence | A rigorous homogenization method for the determination of the overall ultimate strength of periodic discrete media and an application to general hexagonal lattices of beams[END_REF]). However, these approaches lead only to 3D or in-plane macroscopic models. Finally, the work from Lewiski [START_REF] Lewiski | Effective models of composite periodic platesI. Asymptotic solution[END_REF] and Cecchi and Sab [START_REF] Cecchi | A homogenized ReissnerMindlin model for orthotropic periodic plates: Application to brickwork panels[END_REF] suggest homogenization techniques for a thick plate which apply only to microstructures with the classical Cauchy's continuum. Hence, in this respect the present approach is innovative.

One can argue that from an engineering point of view, the full simulation of the lattice remains affordable in terms of computation and will directly lead to more accurate results. However, because these structures are periodic, their design is often based on the assessment of a single unit-cell. Engineers check if this unit-cell sustains highest loads in critical areas and then, set the design for the entire roof. Having a single unit-cell reproduced many times enables returns to scale and reduces costs. Such an approach requires a good knowledge of local stress fields in the unit-cell generated by macroscopic loadings, independently of the configuration of the plate, which is exactly the purpose of homogenization techniques. Hence these approaches are complementary of full simulations.

In Section 3 we consider a square beam lattice (Figure 1) and apply both Bending-Gradient and Cecchi and Sab [START_REF] Cecchi | A homogenized ReissnerMindlin model for orthotropic periodic plates: Application to brickwork panels[END_REF] Reissner-Mindlin homogenization scheme. This very simple pattern will enable the derivation of closed-form solutions of the auxiliary problems which are easy to interpret even if the approach can handle 3D geometries. Let us already point out that, because of the patterns symmetries, the Reissner-Mindlin shear forces stiffness F ∼ (Q α = F αβ γ β , where γ α is Reissner-Mindlin shear strain) is necessarily isotropic1 : F αβ = F δ αβ . However, the pattern must be somehow sensitive to the bending orientation. In order to check this prediction, a comparison with exact solutions of the cylindrical bending of the lattice in two orientations is performed in Section 4. It reveals that only the Bending-Gradient theory is able to capture second order effects both in terms of deflection and local stress fields.

Homogenization of a periodic space frame as a thick plate

Summary of the Bending-Gradient plate model

We consider a linear elastic plate which mid-plane is the 2D domain ω ⊂ R2 . Cartesian coordinates (x 1 , x 2 , x 3 ) in the reference frame (e -1 , e -2 , e -3 ) are used to describe macroscopic fields. At this stage, the microstructure of the plate is not specified.

The membrane stress N αβ , the bending moment M αβ , and shear forces Q α (α, β, γ... = 1, 2) are the usual generalized stresses for plates. Moreover, the main feature of the Bending-Gradient theory is the introduction of an additional static unknown: the gradient of the bending moment R αβγ = M αβ,γ . The 2D third-order 2 tensor R ⌢ complies with the following symmetry: R αβγ = R βαγ . It is possible to derive shear forces Q

-from R ⌢ with: Q α = R αββ . The full bending gradient R ⌢ has six components (taking into account symmetries of indices) whereas Q has two components. Thus, using the full bending gradient as static unknown introduces four additional static unknowns. More precisely: R 111 and R 222 are respectively the cylindrical bending part of shear forces Q 1 and Q 2 , R 121 and R 122 are respectively the torsion part of these shear forces and R 112 and R 221 are linked to strictly self-equilibrated stresses. Equilibrium equations and stress boundary conditions are detailed in Appendix A. They are very similar to those of Reissner-Mindlin theory where

Q α = M αβ,β is replaced by R αβγ = M αβ,γ .
The main difference between Reissner-Mindlin and Bending-Gradient plate theories is that the latter enables the distinction between each component of the gradient of the bending moment whereas they are mixed into shear forces with Reissner-Mindlin theory.

Generalized stresses N αβ , M αβ and R αβγ work respectively with the dual strain variables: e αβ , the conventional membrane strain, χ αβ the curvature and Γ αβγ the third order tensor related to generalized shear strains. These strain fields must comply with the compatibility conditions and boundary conditions detailed in Appendix A.

Finally, assuming uncoupling between (N ∼ , M ∼ ) and R ⌢ (see [START_REF] Lebée | A Bending-Gradient model for thick plates. Part I: Theory[END_REF][START_REF] Lebée | Homogenization of thick periodic plates: Application of the Bending-Gradient plate theory to a folded core sandwich panel[END_REF]), the Bending-Gradient plate constitutive equations are written as:

       N ∼ = A ∼ ∼ : e ∼ + B ∼ ∼ : χ ∼ (1a) M ∼ = T B ∼ ∼ : e ∼ + D ∼ ∼ : χ ∼ (1b) Γ ⌢ = f ⌢ ⌢ • • • R ⌢ (1c)
where A ∼ ∼ , B ∼ ∼ , D ∼ ∼ are conventional Kirchhoff-Love fourth order stiffness tensors and ":" is the double contraction product. The generalized shear compliance tensor f ⌢ ⌢ is a sixth order tensor3 and " • • • " denotes a triple contraction product:

f ⌢ ⌢ • • • R ⌢ = (f αβγδǫζ R ζǫδ ).
In some cases, the Bending-Gradient is turned into a Reissner-Mindlin plate model. This is the case for homogeneous plates. In order to estimate the difference between both plate models we defined the isotropic projection of the Bending-Gradient stress energy density on a Reissner-Mindlin one in [START_REF] Lebée | A Bending-Gradient model for thick plates. Part I: Theory[END_REF]. According to this projection, the Reissner-Mindlin part of f ⌢ ⌢ is:

f ⌢ ⌢ RM = 2 3 i ∼ ∼ • i ∼ ∼ • • • f ⌢ ⌢ • • • 2 3 i ∼ ∼ • i ∼ ∼ (2)
where i ∼ ∼ is the identity for in-plane elasticity tensors (i αβγδ = 1 2 (δ αγ δ βδ + δ αδ δ βγ )) and i ∼ ∼ •i ∼ ∼ = (i αβγη i ηδǫζ ) is a sixth order tensor. Consequently, we suggested the following relative distance between the Bending-Gradient and the Reissner-Mindlin stress energy densities:

∆ RM/BG = f ⌢ ⌢ -f ⌢ ⌢ RM f ⌢ ⌢ , where f ⌢ ⌢ = f αβγδǫζ f αβγδǫζ (3)
The scalar ∆ RM/BG gives an estimate of the pure warping fraction of the shear stress energy and is a criterion for assessing the need for the Bending-Gradient model: 0 ≤ ∆ RM/BG < 1. When the plate constitutive equation reduces to a Reissner-Mindlin one, we have exactly ∆ RM/BG = 0 and both theories are strictly identical. Finally, if f ⌢ ⌢ writes as:

f ⌢ ⌢ = i ∼ ∼ • f ∼ • i ∼ ∼ where f
∼ is a second order tensor, then ∆ RM/BG = 0 and the Reissner-Mindlin compliance is exactly f

∼ . The constitutive tensors A ∼ ∼ , B ∼ ∼ , D ∼ ∼ and f ⌢ ⌢
are function of the plate microstructure. Their derivation is based on the homogenization scheme discussed in the next section.

Homogenization scheme

When the plate is laminated, closed-form expression for

A ∼ ∼ , B ∼ ∼ , D ∼ ∼ and f ⌢ ⌢
where derived in [START_REF] Lebée | A Bending-Gradient model for thick plates. Part I: Theory[END_REF]. When the plate is periodic those tensors must be derived using an homogenization procedure. This procedure relies on the fundamental assumption that loadings as well as macroscopic variables (M ∼ , M ∼ , R ⌢ ) varies slowly compared to the size of the microstructure (the unit-cell). This means that at the scale of the unit-cell the macroscopic fields can be considered as uniform. Hence, it is possible to define auxiliary problems which enable the derivation of the actual microscopic fields as a linear superposition of macroscopic fields. Once this is done, it is straightforward to take the average of the elastic energy stored inside the unit-cell and equate it to the macroscopic elastic energy density (namely the macroscopic constitutive equation). This analogy between micro and macroscopic energy is usually called Hill-Mandel principle. When the plate microstructure is constituted of a 3D continuum, the auxiliary problems were introduced in [START_REF] Lebée | Homogenization of thick periodic plates: Application of the Bending-Gradient plate theory to a folded core sandwich panel[END_REF]. In the present case, we consider periodic plates which unit-cell is constituted of connected beams. This necessitates the adaptation of the homogenization scheme for 3D continuum. In this section, we first describe the unitcell and the related beam model. Then the Kirchhoff-Love auxiliary problem which enables the derivation of A ∼ ∼ , B ∼ ∼ , D ∼ ∼ is adjusted to this type of plate. Finally, the Bending-Gradient auxiliary problem is derived through a direct approach and defines the compliance f ⌢ ⌢ .

The unit-cell of the space frame

In this section, we give a general description of the space frame and set the main assumptions. We also introduce some useful definitions.

We consider an in-plane periodic plate which unit-cell is made of an assembly of connected beams (Figure 2). The unit-cell reference frame is (e -1 , e -2 , e -3 ) associated with the local coordinate y -= (y 1 , y 2 , y 3 ). There are n b beams constituting the unit-cell intersecting at nodes (see Figure 2). Each beam connects two nodes through a path γ k . The collection of γ k is Γ = ∪ n b k=1 γ k . Without loss of generality, we assume here a perfect connection at nodes (fully clamped). Furthermore, a subset of these nodes are pairs of periodicity nodes: they connect the unit-cell to its neighbors. Thus interior nodes are denoted I i (i = 1, n i ) and periodicity pairs P j and P j , (j = 1, n p ).

Local beam equations. The local reference frame along each beam is denoted (t -, n -, b -) associated with the coordinate s -= (s 1 , s 2 , s 3 ). The unit vector t -is oriented in the direction of the neutral axis of the beam. For the sake of simplicity, beams are assumed rectilinear and follows St Venant's approximation with uniform torsion. Consequently, the resultant force r -and the moments m -are defined as follows:

r -= S σ ∼ B • t -dS and m -= S s -× σ ∼ B • t -dS (4)
where S is the beam section, σ ∼ B is the beam's local 3D stress and × denotes the cross product. Each beam is subjected to distributed forces k -and torques l -. The corresponding equilibrium equations are:

r - ′ + k -= 0 m - ′ + t -× r -+ l -= 0 (5)
Where " ′ " denotes the derivative with respect to s 1 . The constitutive equation writes as follows:

r -= E ∼ • u - ′ + t -× θ - m -= G ∼ • θ - ′ (6)
where u -is the beam displacement and θ -its rotation. E ∼ and G ∼ are the beam stiffness tensors. For instance, when there are two axis of symmetry for the section, these tensors write in the local reference frame of the beam (t -, n -, b -), as:

E ∼ =   ES 0 0 0 GS s,2 0 0 0 GS s,3   and G ∼ =   GJ 0 0 0 EI 2 0 0 0 EI 3   ( 7
)
where E is the Young modulus, G the shear modulus of the constitutive material of the beam. S is the section area, S s,2 and S s,3 are the shear areas in each direction. J is the torsion constant, I 2 = S s 2 3 dS and I 3 = S s 2 2 dS are the second moments of inertia.

Kinematic and static constraints on the assembly. Since we assumed perfect connection, displacements and rotations must be continuous at interior nodes I i in all the following. Moreover, local equilibrium at each interior node must be fulfilled:

∀j = 1...n i , n bj k=1 ǫ k k r -(I j ) = 0 and n bj k=1 ǫ k k m -(I j ) = 0 (8)
where n bj is the number of beams connecting at node I j , k r -denotes the resultant of beam k connected to the node. Moreover, ǫ k = +1 if k t -, the tangent vector of beam k, is directed toward the node and

ǫ k = -1 if k t
-is directed away the node.

Additionally we define periodicity conditions as:

∀j = 1...n p , θ -(P j ) = θ -P j and u -(P j ) = u -P j (9)
Elastic energy. The complementary elastic energy per unit surface stored inside the unit-cell will be used in the homogenization scheme. It is defined as the sum along Γ of the beams stress energy density:

w * int r -, m -= 1 2A Γ T r -• E ∼ -1 • r -+ T m -• G ∼ -1 • m -dΓ ( 10
)
where A is the in-plane area of the unit-cell

Homogenization as Kirchhoff-Love plate

In this section, we recast the Kirchhoff-Love auxiliary problem -originally defined for 3D continuumin a suitable form for the unit-cell we just introduced.

In the case of a 3D continuum, the Kirchhoff-Love auxiliary problem consists in applying the global displacement field

U - KL = ê ∼ • y -+ y 3 χ ∼ • y -- 1 2 T y -• χ ∼ • y - e -3 (11)
on average on the unit-cell, taking into account upper and lower free faces of the plate [START_REF] Lebée | Homogenization of thick periodic plates: Application of the Bending-Gradient plate theory to a folded core sandwich panel[END_REF]. Here, we introduced the convenient notation:

Then, the elastic energy stored in the unit-cell is computed for each pair of component of e ∼ and χ ∼ and yields the Kirchhoff-Love stiffness tensors

A ∼ ∼ , B ∼ ∼ , D ∼ ∼ .
In the present case, there is an additional displacement field which carries energy in the unit-cell: the 3D rotation field θ -. Rotations are related to the skew-symmetric part of the gradient of the displacement. The rotation pseudo-vector related to the skew symmetric part of U - KL ⊗ ∇ -is :

Θ - KL = ω ∼ • χ ∼ • y - (13)
where

ω ∼ = e -2 ⊗ e -1 -e -1 ⊗ e -2
is the permutation operator. Thus we suggest to apply both U - KL and Θ - KL on average on the unit-cell in the case of a space frame. This leads to the following Kirchhoff-Love auxiliary problem P KL :

P KL :                                  r - KL ′ = 0 on Γ m - KL ′ + t -× r - KL = 0 on Γ r - KL = E ∼ • u - KL ′ + t -× θ - KL on Γ m - KL = G ∼ • θ - KL ′ on Γ u - KL = u - per + U - KL on Γ θ - KL = θ - per + Θ - KL on Γ ∀i = 1...n i , u -
per and θ - per continuous at node I i ∀i = 1...n i , r - KL and m - KL equilibrated at node I i ∀j = 1...n p , u - per , θ - per , r - KL and m - KL periodic at node P j [START_REF] Lebée | Homogenization of thick periodic plates: Application of the Bending-Gradient plate theory to a folded core sandwich panel[END_REF] where the subscript "per" refers to periodicity (u - per and θ - per are also named corrector fields). NB: The free faces boundary condition is implicitly included in the present auxiliary problem. Actually, from one node to another, beams are free. In addition, no vertical displacement is prescribed at nodes. However, the displacement of one node belonging to the reference plane of the plate must be blocked in order to prevent rigid motion.

Solving the problem for each individual component of e ∼ and χ ∼ leads to the stress localization tensors r e iαβ , m e iαβ and r χ iαβ , m χ iαβ . The complete local stress field can be reconstructed by linear combination:

r - KL = r - (e) + r - (χ) = r e iαβ e αβ + r χ iαβ χ αβ = r ∼ - e : e ∼ + r ∼ - χ : χ ∼ (15a) m - KL = m - (e) + m - (χ) = m e iαβ e αβ + m χ iαβ χ αβ = m ∼ - e : e ∼ + m ∼ - χ : χ ∼ (15b)
Kirchhoff-Love stiffness tensors are derived using Hill-Mandel principle. The Kirchhoff-Love plate energy density is identified to the energy stored in the unit-cell (Equation 10):

w * int r - KL , m - KL = 1 2 e ∼ : A ∼ ∼ : e ∼ + 2e ∼ : B ∼ ∼ : χ ∼ + χ ∼ : D ∼ ∼ : χ ∼ (16)
Kirchhoff-Love plate moduli are then evaluated as follows:

A ∼ ∼ = 1 A Γ T r ∼ - e • E ∼ -1 • r ∼ - e + T m ∼ - e • G ∼ -1 • m ∼ - e dΓ (17a)
B ∼ ∼ = 1 A Γ T r ∼ - e • E ∼ -1 • r ∼ - χ + T m ∼ - e • G ∼ -1 • m ∼ - χ dΓ (17b)
D ∼ ∼ = 1 A Γ T r ∼ - χ • E ∼ -1 • r ∼ - χ + T m ∼ - χ • G ∼ -1 • m ∼ - χ dΓ (17c)
Finally, using the inverted constitutive equation,

e ∼ = a ∼ ∼ : N ∼ + b ∼ ∼ : M ∼ (18a) χ ∼ = T b ∼ ∼ : N ∼ + d ∼ ∼ : M ∼ (18b)
we get the following local stress field generated by uniform membrane stress N ∼ and bending moment M ∼ on the unit-cell:

r - (N ) = r ∼ - N : N ∼ = r ∼ - e : a ∼ ∼ + r ∼ - χ : T b ∼ ∼ : N ∼ (19a) m - (N ) = m ∼ - N : N ∼ = m ∼ - e : a ∼ ∼ + m ∼ - χ : T b ∼ ∼ : N ∼ (19b) r - (M) = r ∼ - M : M ∼ = r ∼ - e : b ∼ ∼ + r ∼ - χ : d ∼ ∼ : M ∼ (19c) m - (M) = m ∼ - M : M ∼ = m ∼ - e : b ∼ ∼ + m ∼ - χ : d ∼ ∼ : M ∼ (19d)

Homogenization as Bending-Gradient plate

Once the Kirchhoff-Love homogenization is performed, it becomes possible to derive the Bending-Gradient localization field and constitutive equation. In [START_REF] Lebée | Homogenization of thick periodic plates: Application of the Bending-Gradient plate theory to a folded core sandwich panel[END_REF], the homogenization scheme was suggested for 3D Cauchy continuum. However, in the present case, the unit-cell is constituted of structural elements. We suggest here a direct approach similar to the one used for deriving laminated plates localization fields [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF].

Provided M ∼ is uniform in the unit-cell, the localization fields r - (M) , m - (M) in ( 19) follow the equilibrium equations of Problem 14. Let us assume a macroscopic linear variation of M ∼ through the unit-cell: M * αβ = M αβ + R αβγ y γ and define the Bending-Gradient localization field as:

r - BG = r - (N ) + r ∼ - M : M ∼ * + r - (R) (20)
and

m - BG = m - (N ) + m ∼ - M : M ∼ * + m - (R) (21)
Injecting this field in Problem 14 leads to the definition of the Bending-Gradient auxiliary problem:

                         r - (R) ′ + k - (R) = 0 on Γ m - (R) ′ + t -× r - (R) + l - (R) = 0 on Γ r - (R) = E ∼ • u - (R) ′ + t -× θ - (R) on Γ m - (R) = G ∼ • θ - (R) ′ on Γ ∀i = 1...n i , u - (R) and θ - (R) continuous at node I i ∀i = 1...n i , r - (R) and m - (R) equilibrated at node I i ∀j = 1...n p , u - (R) , θ - (R) , r - (R) and m - (R) periodic at node P j (22)
where

k - (R) = r ∼ - M : R ⌢ • t -= r M iαβ t γ R βαγ (23a) l - (R) = m ∼ - M : R ⌢ • t -= m M iαβ t γ R βαγ (23b)
are the distributed loadings generated by a uniform gradient of the bending moment R ⌢ .

Solving the problem for each individual component of R αβγ leads to the localization stress fields r R iγβα , m R iγβα . The complete local stress field can be reconstructed by linear combination:

r - (R) = r R iγβα R αβγ = r ⌢ - R • • • R ⌢ (24) m - (R) = m R iγβα R αβγ = m ⌢ - R • • • R ⌢ (25)
The generalized shear compliance is then evaluated as follows:

f ⌢ ⌢ = 1 A Γ T r ⌢ - R • E ∼ -1 • r ⌢ - R + T m ⌢ - R • G ∼ -1 • m ⌢ - R dΓ (26)
3. Homogenization of a square lattice as a Bending-Gradient plate

In order to illustrate this homogenization scheme, we consider a very simple structure: a square beam lattice (Figure 1). In this section, we derive the homogenized behavior of the lattice as a Bending-Gradient plate and also as a Reissner-Mindlin plate.

The unit-cell

Figure 3 shows the unit-cell of the lattice idealization. It is constituted of two beams ("h": horizontal, "v": vertical) which lie in the (e -1 , e -2 ) plane. Hence their local reference frames write as (e -1 , e -2 , e -3 ) for beam "h" and (e -2 , -e -1 , e -3 ) for beam "v". Both beams are identical with length b. We assume the beam section has two orthogonal axis of symmetry. One axis is e -3 so that the beam stiffness tensors write as Equation 7in the local reference frame.

Kirchhoff-Love auxiliary problem

Computation details of Kirchhoff-Love auxiliary problem are given in Appendix B. Here, we give directly the generalized stress field localization as function of N ∼ and M ∼ in the local reference frame of each beam:

h r - KL =   bN 11 bN 12 0   h and h m - KL =   -bM 12 bM 11 bN 12 s 1 -b 2   h (27a)
  v and v m - KL =   bM 12 bM 22 -bN 12 s 1 -b 2   v (27b)
This result is rather intuitive: membrane traction in Direction 1 generates only uniform traction of beam "h", as well as N 22 generates traction in beam "v". In plane shear N 12 puts both beams in non-uniform in-plane bending. Similarly, pure bending in Direction 1 generates pure out-of-plane bending of beam "h" as well as M 22 generates pure out-of-plane bending in beam "v". Finally, plate torsion M 12 generates pure torsion in both beams. NB : All localization fields are continuous except the one generated by N 12 . This is not incorrect, the generalized stress discontinuity being introduced by the node. Actually, the node is equilibrated.

The corresponding Kirchhoff-Love stiffness tensors are:

A ∼ ∼ =    ES b 0 0 0 ES b 0 0 0 b GSs,2 + b 3 12EI3 -1    , B ∼ ∼ = 0 and D ∼ ∼ =   EI2 b 0 0 0 EI2 b 0 0 0 GJ b   ( 28
)
where we introduced Voigt notation A ∼ ∼ for fourth order tensors A ∼ ∼ as:

A ∼ ∼ =   A 1111 A 2211 √ 2A 1211 A 2211 A 2222 √ 2A 1222 √ 2A 1211 √ 2A 1222 2A 1212   (29)
NB : the corresponding notation for N ∼ and M ∼ is:

N ∼ =   N 11 N 22 √ 2N 12   (30)
There is no coupling between N ∼ and M ∼ (B ∼ ∼ = 0) and there is no Poisson's effect (A 1122 = D 1122 = 0). 23. The beam stress localization field related to M ∼ is given in the local reference frame by:

h r - (M) = v r - (M) = 0 -, h m - (M) =   -bM 12 bM 11 0   h and v m - (M) =   bM 12 bM 22 0   v (31) First, we have directly h k - (R) = v k - (R) = 0 -. Then, for beam "h" we have R ⌢ • t -= (R αβ1 ). Replacing M αβ by this in h m - (M) leads to: h l - (R) =   -bR 121 bR 111 0   h (32)
Similarly for beam "v" we have:

v l - (R) =   bR 122 bR 222 0   v ( 33
)
The solution of Problem 22 is derived a very similar way as for Kirchhoff-Love auxiliary problem. It is not detailed here. The generalized stress localization in the local reference frame of each beam is:

h r - (R) =   0 0 b (R 111 + R 122 )   h and h m - (R) =   bR 121 s 1 -b 2 bR 122 s 1 -b 2 0   h (34a) v r - (R) =   0 0 b (R 121 + R 222 )   v and v m - (R) =   -bR 122 s 1 -b 2 bR 121 s 1 -b 2 0   v (34b)
Let recall that

Q 1 = R 111 + R 122 and Q 2 = R 121 + R 222 .
Thus the resultants h r - (R) and v r - (R) are directly function of the shear forces. However, the torsion part of shear forces (R 122 and R 121 ) generates additional torsion and bending in h m - (R) and v m - (R) . Thus, shear forces related to cylindrical bending and torsion do not generate the same local stress field.

Reissner-Mindlin localization

In order to compare in details the Bending-Gradient theory with Reissner-Mindlin theory, we also considered a Reissner-Mindlin homogenization scheme. We chose Cecchi and Sab [START_REF] Cecchi | A homogenized ReissnerMindlin model for orthotropic periodic plates: Application to brickwork panels[END_REF] method for deriving a Reissner-Mindlin constitutive equation. Following the work from Whitney [START_REF] Whitney | Stress Analysis of Thick Laminated Composite and Sandwich Plates[END_REF], Cecchi and Sab [START_REF] Cecchi | A homogenized ReissnerMindlin model for orthotropic periodic plates: Application to brickwork panels[END_REF] assume the plate is under cylindrical bending for deriving local stress fields and a Reissner-Mindlin constitutive equation. The derivation is not detailed here. However, the related Reissner-Mindlin field localization is easily obtained, setting

R 111 = Q 1 , R 222 = Q 2 and R 121 = R 122 = R 112 = R 221 = 0 in Equation 34: h r - (Q) =   0 0 bQ 1   h and h m - (Q) = 0 - (35a) v r - (Q) =   0 0 bQ 2   v and v m - (Q) = 0 - (35b)
Thus, there is no more correction of m -.

Constitutive equations

Bending-Gradient. The Bending-Gradient compliance is derived using Equation 26 and is:

f ⌢ ⌢ = b GS s,3 i ∼ ∼ • i ∼ ∼ + b 3 12 1 GJ + 1 EI 2 t ⌢ ⌢ (36)
where

i ∼ ∼ • i ∼ ∼ =         1 0 0 0 0 1/ √ 2 0 0 0 0 0 0 0 0 1/2 0 1/ √ 2 0 0 0 0 0 0 0 0 0 1/ √ 2 0 1 0 1/ √ 2 0 0 0 0 1/2         and t ⌢ ⌢ =        
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/2

        (37)
in the main reference frame. We used the following Voigt notation (see [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF] for full details):

f ⌢ ⌢ =         f 111111 f 111122 √ 2f 111121 f 111211 f 111222 √ 2f 111221 f 221111 f 221122 √ 2f 221121 f 221211 f 221222 √ 2f 221221 √ 2f 121111 √ 2f 121122 2f 121121 √ 2f 121211 √ 2f 121222 2f 121221 f 112111 f 112122 √ 2f 112121 f 112211 f 112222 √ 2f 112221 f 222111 f 222122 √ 2f 222121 f 222211 f 222222 √ 2f 222221 √ 2f 122111 √ 2f 122122 2f 122121 √ 2f 122211 √ 2f 122222 2f 122221         (38)
In order to go further in the interpretation, we rewrite f ⌢ ⌢ as:

f ⌢ ⌢ = b GS s,3 i ∼ ∼ • i ∼ ∼ + η t ⌢ ⌢ (39)
where

η = b 2 GS s,3 12
1 GJ + 1 EI 2 (40)
Reissner-Mindlin. The Reissner-Mindlin shear force compliance is derived using Equation 26 again where the Reissner-Mindlin localization field (35) is plugged:

f αβ = b GS s,3 δ αβ (41)
which is the surface average of individual beam shear force stiffness.

Difference between Reissner-Mindlin and Bending-Gradient constitutive equation

Setting η = 0 in f ⌢ ⌢
leads exactly to the Reissner-Mindlin constitutive equation (41). Thus the additional

term η t ⌢ ⌢ in f ⌢ ⌢
is a contribution purely related to the Bending-Gradient model. It reflects the additional torsion and bending in the beams when shear forces are not purely cylindrical in the main reference frame (eg: Q 1 = M 11,1 ). In contrast to i ∼ ∼ • i ∼ ∼ which is an isotropic tensor for the Bending-Gradient constitutive equation (eg. invariant to in-plane rotations [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF]), t ⌢ ⌢ is not isotropic which makes f

⌢ ⌢
sensitive to the orientation of the beam lattice.

Let us give some order of magnitude for η. For a circular cross section, we have: S s,3 = 9 40 πd 2 , I 2 = πd 4 64 and J = 2I 2 where d is the beam diameter. Thus:

η = 3b 2 5d 2 2 + ν 1 + ν ≈ b 2 d 2 (42)
The ratio b/d is typically larger than 4 (otherwise a 3D modelization of the unit-cell would be more suitable) and leads to η typically larger than 10. Clearly, the contribution of η t ⌢ ⌢ to f ⌢ ⌢ is dominant. This point is confirmed when deriving the relative distance between Bending-Gradient and Reissner-Mindlin model introduced in Subsection 2.1:

∆ RM/BG = 2 √ 2 3 η 8 + (1 + η) 2 (43) For large η, ∆ RM/BG = 2 √ 2 3 ≈ 94%. Furthermore, for η ≫ 1 the eigenvalues of f ⌢ ⌢ are approximately: b GS s,3 η 2 , η 2
, 1, 1, 0, 0 Thus, depending on which eigenvalue of f ⌢ ⌢ is activated by the plate's configuration, the results might be quite different both in deflection and in local stress fields as will be illustrated in the next section.

Cylindrical bending of a square beam lattice

In order to give a clear demonstration of the additional warping effects included in the Bending-Gradient theory we compare Kirchhoff-Love, Reissner-Mindlin and Bending-Gradient models to the exact solution in cylindrical bending configuration where the lattice direction is rotated by 0 and 45 with respect to the bending direction.

Lattice at 0

Figure 4 shows the lattice configuration when the bending direction is the same as the one of the lattice. It is bent in Direction x 1 along 2n unit-cells and is infinitely periodic in Direction x 2 . Longitudinal beams ("h") are transversely loaded with the sine shaped load k -= -k 0 sin κx 1 e -3 , where κ = π L is the wave number and L is the span: L = 2nb.

. Exact solution

The exact solution for this configuration is straightforward. Transverse beams ("v") in Direction 2 are not loaded and longitudinal beams are simply supported with the transverse load k -. The deflection along longitudinal beams is:

h u ex 3 = - k 0 κ 4 EI 2 1 + π 2 EI 2 L 2 GS s,3 sin κx 1 (44)
the bending moment is:

h m ex 2 = - k 0 κ 2 sin κx 1 (45)
and the shear force is:

h r ex 3 = - k 0 κ cos κx 1 (46)

Homogenization approximations

The plate surrogate for the configuration of the beam lattice given in Figure 4 is a simply supported plate under cylindrical bending. The plate distributed load is p -= k 0 /b sin κx 1 e -3 in order to have the same loading per unit surface.

The boundary value problem is stated as follows. We assume invariance in x 2 direction. The plate is simply supported in x 1 = (0, L): U 3 (0, L) = 0, M ∼ (0, L) = 0 and loaded with p -. We search for fields which are compatible (Equations (A.1) and (A.2)) and fulfill the constitutive equation [START_REF] Ciarlet | Justification Of The 2-Dimensional Linear Plate Model[END_REF]. Computations are fully detailed in Appendix D.1. The Bending-Gradient and the Reissner-Mindlin solutions are both strictly identical to the exact solution. The Kirchhoff-Love approximation yields a less detailed result. The macroscopic deflection is:

U KL 3 = - k 0 κ 4 EI 2 sin κx 1 (47)
There is no localization related to shear forces and then:

h r KL 3 = 0 (48) h m KL 2 = - k 0 κ 2 sin κx 1 (49)
As expected, the Kirchhoff-Love homogenization scheme underestimates the deflection because it does not take into account second order effects related to R ⌢ . Additionally it does not capture local stress generated by shear forces. However, it is common practice for laminated plates to use the 3D equilibrium equation integrated through the thickness in order to derive this field. Here, we can use the microstructure's beam equilibrium equation [START_REF] Reddy | On refined computational models of composite laminates[END_REF] for this: plugging m KL 2 into equation [START_REF] Reddy | On refined computational models of composite laminates[END_REF] gives the correct shear force in the beam.

Lattice at 45

Now we simply rotate the lattice of 45 with respect to the bending direction. Figure 5 shows the lattice in this new configuration. It is bent in Direction x 1 along 2n diamond shaped unit-cells and is infinitely periodic in direction x 2 . The beams are transversely loaded with the sine shaped load k -= -k 0 sin κx 1 e -3 . Now the span is: L = 2 √ 2nb (Figure 6). There are two sets of beams. The first set "+" is numbered 1..2n and is rotated by +π/4 in (e -1 , e -2 ) plane and the second set "-" is rotated by -π/4. It is symmetric to set + with respect to any (A j , e -1 , e -3 ) plane. Hence the focus is on set + in the following.

Beams are perfectly connected at nodes A j located at x 1 = x j = jb/ √ 2. Simple support condition occurs at node A 0 and A 4n : both beams are perfectly connected together, however the node is free to rotate with respect to the support.

Exact solution

The derivation of the exact solution of this problem, though cumbersome, does not raise any difficulty. It is detailed in Appendix C. The results are summarized here for the set of beams "+". The generalized stresses in the local reference frame of the beam for x j < x 1 < x j+1 are:

+ r - ex = - √ 2 k 0 κ cos κx 1 e -3 (50)
+ m - ex = -m *    α 1+α cosj -cosj+1 2 √ 3γ sin κx 1 + α 1+α cosj+1 -cosj 2 √ 3γ 0    + ( 51
)
where cos j = cos κx j , α = GJ EI2 is the ratio between torsion and bending of the beam, γ 2 = κ 2 b 2 24 is the unit-cell size versus span ratio and m * = 2k0 κ 2 is the maximum bending moment in the beams. The resultant + r - ex is rather similar to the 0 case. However, with + m - ex , there is an additional uniform torsion in each beam.

The deflection along beams "+" for x j < x 1 < x j+1 is:

+ u ex 3 = -U KL 1 + γ 2 (1 + α) ηα (1 + α) sin κx 1 - √ 3γα cos j+1 (x 1 -j) 2 -cos j (x 1 -j -1) 2 + cos j -cos j+1 1 -cos 1 ( 52
)
where

x 1 = √ 2x 1 /b is the normalized coordinate and U KL = 4k0 κ 4 EI2(1+α) is Kirchhoff-Love maximum deflection.

Homogenization approximations

The plate surrogate for the configuration in Figure 5 is a simply supported plate under cylindrical bending where the distributed load is now p -= 2k0 b sin κx 1 e -3 in order to have the same loading per unit surface. The boundary value problem is exactly the same as in the 0 case. Hence the approximation is derived following a similar procedure. However one have to take into account the rotation of the lattice. Full details are given in Appendix D.2 Deflection. The macroscopic deflection for the Bending-Gradient is:

U BG 3 = -U KL       1 + γ 2 α (1 + α) 2 η RM + α 2 1 + γ 2 BG       sin κx 1 (53)
Removing underbraced terms with "BG" leads to Reissner-Mindlin approximation and removing both "BG" and "RM" terms leads to Kirchhoff-Love deflection.

For most sections α ≈ 1 and γ is small. Since η is typically larger than 10, the "BG" term in Equation 53 is large compared to the "RM" one. Thus, the Reissner-Mindlin approach underestimates the actual second order deflection. This is confirmed when plotting the deflection error, U3-u ex 3 U KL along the span (Figure 7). Beams have the same circular cross section as in Section 3.3. The Poisson's coefficient is 0.3 and the ratio between the cell-size and the section diameter is b/d = 5, which leads to η ≃ 26. The span is made of four unit-cells (n = 2, 0 < x 1 < 4n). Only half the span is plotted versus the normalized coordinate x 1 = √ 2x 1 /b. Globally, the errors remain small compared to the exact solution. As expected Kirchhoff-Love approximation underestimates the deflection. However, the correction given by Reissner-Mindlin approximation is not large enough. Actually, only the Bending-Gradient is able to properly capture all corrections related to second order effects.

Local fields. Local stress fields in the beams with the Bending-Gradient approximation write as follow:

+ r - BG = - √ 2k 0 κ cos κx 1 e -3 RM (54) 0 0.5 1 1.5 2 2.5 3 3.5 4 -12 -10 -8 -6 -4 -2 0 2 x 10 -3
Span coordinate x 1 Deflection error,

u ex 3 -U3 U KL BG RM KL
+ m - BG = - m * 1 + α           α sin κx 1 - α 1 + γ 2 γ 2 sin κx 1 + √ 3γ (x 1 -j -1/2) cos κx 1 BG sin κx 1 - 1 1 + γ 2 -αγ 2 sin κx 1 + 2 √ 3αγ (x 1 -j -1/2) cos κx 1 BG 0           + (55)
where the same convention as in Equation 53 is used. Exactly as in the 0 case (Section 4.1), Kirchhoff-Love approximation does not capture the shear forces r 3 in the beams. In addition, using the beam equilibrium for recovering the effect shear forces as we did in Subsubsection (4.1.2) leads to + r -= -1 1+α √ 2k0

κ cos κx 1 e -3 which underestimates the actual shear force. This fact was already pointed out in [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF] when dealing with laminated plates. The FE simulation using integration through the thickness for deriving transverse shear stress distribution was leading to un-balanced fields when the direction of orthotropy was not the same as the direction of bending. Considering now the Reissner-Mindlin approximation, there is an additional correction for m - BG which is only carried by the Bending-Gradient solution.

In order to have a closer look at this contribution, the bending moment m 2 and the torsion m 1 along beams "+" are plotted respectively in Figures 8 and9 for the exact solution and homogenization approximations. Kirchhoff-Love and Reissner-Mindlin solutions are identical and plotted as the same continuous curve which follows the exact solution on average. By contrast, the Bending-Gradient approximation follows closely the exact solution. Especially, it captures the change of sign of the bending moment m 2 close to the support. In case the lattice is made of concrete, this kind detail is useful for an efficient design. This is not included in the other approximations.

The errors

m1-m ex 1 m * and m2-m ex 2 m *
are also plotted in Figures 10 and11. It shows that the maximum error is 1% for the Bending-Gradient and 9% for the other approaches.

Discussion

In the previous section, three homogenization approximations of the cylindrical bending of a beam lattice were derived. The first one, leading to a Kirchhoff-Love plate model is also known to be the first order of the asymptotic expansion [START_REF] Ciarlet | Justification Of The 2-Dimensional Linear Plate Model[END_REF][START_REF] Caillerie | Thin elastic and periodic plates[END_REF]. This approach enables to have a first order estimate of the macroscopic deflection as well as local stress fields. Except with sandwich panels and laminated plates, in most application the first order deflection is accurate enough. This is also the case in this paper (the error is already smaller than 10 -3 for a span of few unit-cells). However, this approach does not capture the local effect of shear forces on the microstructure and the common practice for deriving the actual effect of transverse shear appears clearly as not reliable in more general configurations. The Reissner-Mindlin homogenization approach used in this paper clearly misses some part of their effect. This is because of the assumption of cylindrical bending which is made in most approach for deriving a Reissner-Mindlin-like constitutive equation. This assumption enforces Q 1 = M 11,1 and precludes any possible effect of the other derivatives of the bending moment on the local stress on the microstructure.

In the present case, the Bending-Gradient approximation captures all second order effects and gives a very satisfying description of both the deflection and the local stress fields. Especially, when rotating the lattice, this accuracy is preserved. Moreover, it brought out two scaling parameters between the first order deflection and the second order deflection. The first one is related to shear forces in the beams and is the same as the one in the Reissner-Mindlin approach. It is proportional to γ 2 /η, namely (d/L) 2 . The second one is related to the second order torsion/bending of the beams and is proportional to γ 2 , namely (b/L) 2 . Thus the latter effect is dominant at second order provided it is activated by the plate configuration. This distinction between scaling parameters is completely discarded when using a Reissner-Mindlin approach.

Our approach still presents some limitations. First, in the present case, we chose supports conditions for the lattice so that edge effects remain limited. It is well known in homogenization that these effects might be not negligible in some configurations. This point is far out of the scope of this paper and raises the question -common to all models -regarding the mechanical meaning of plate boundary conditions. Second, the local load distribution on the microstructure was arbitrarily chosen. At first order, this choice does not affect the fields localization (Kirchhoff-Love case). However, different load distributions could generate different localization of stress fields when looking for a higher order approximation (see [START_REF] Boutin | Microstructural effects in elastic composites[END_REF] for instance).

Finally, in the present paper, the plate model was solved in a simple, closed-form solution. This raises the question regarding the finite element implementation of the Bending-Gradient theory. Actually, it is very close to the classical Reissner-Mindlin model (with its well-known limitations). The most straightforward implementation is to add the four C 0 generalized rotation fields which are coupled to Reissner-Mindlin FEM formulation. This way the computational burden remains perfectly acceptable.

Conclusion

In this paper, the Bending-Gradient homogenization scheme was extended to space frames. Its application to a square beam lattice led to closed form solutions of the auxiliary problems. It also revealed that this microstructure cannot be properly modeled with Reissner-Mindlin theory. This observation was confirmed when comparing homogenization approaches with the exact solution of the cylindrical bending of the lattice. Only the Bending-Gradient theory is able to capture second order effects both in macroscopic deflection and in local fields.

Many other periodic structures might not be correctly modeled with Reissner-Mindlin theory such as concrete waffle slabs, orthotropic decks or composite floor systems and requires more accurate approaches such as the Bending-Gradient model. Hence, let us point out that it is perfectly possible to model a microstructure which is a mix between 3D, shell or beam elements following the direct derivation of the Bending-Gradient auxiliary problem suggested in the present paper and that other patterns are currently under investigation.

One can ask also if the Bending-Gradient theory can be easily extended to more random microstructures. The answer is actually the same as for case of a 3D random medium: it might be possible to establish bounds and argue on the statistical distribution of the microstructure. However, there is no clear argument for expecting that the Bending-Gradient would turn back to a Reissner-Mindlin model.

Finally, provided the limitations pointed out in Section 5 are kept in mind, the Bending-Gradient plate theory is considered by the authors as a sound extension of the Reissner-Mindlin theory to the case of heterogeneous plates. More precisely, it is an attempt to concentrate most of the results given by the second order of the asymptotic expansion [START_REF] Lewiski | Effective models of composite periodic platesI. Asymptotic solution[END_REF] into a model which keeps most aspects of Reissner-Mindlin model, well understood by engineers.

Beam equilibrium. Equilibrium equations (5) applied to beam "h" lead to uniform resultant along the beam:

h r -= h r - 0 (B.2)
and linear bending moment:

h m -= -e -1 × h r -(s 1 -b/2) + h m - b 2 (B.3)
where h m - b 2 is the bending moment taken at s 1 = b/2. Similar results occur for beam "v". We focus now on the derivation of h r - 0 and h m - b 2 which fully determine the stress state of beam "h". Inserting this in node equilibrium B.1 leads to:

e -1 × h r - 0 + e -2 × v r - 0 = 0 (B.4)
Taking the cross product with e -1 enforces: h r 0 2 =v r 0 2 and h r 0 3 = v r 0 3 = 0, where the components are taken in the local reference frame of each beam.

Beam displacement. The constitutive equation 6 related to rotations leads to:

h G ∼ • h θ --h θ - 0 = e -1 × h r - 0 (s 1 -b) s 1 2 + h m - b 2 s 1 (B.5)
where h θ - 0 = h θ -(s 1 = 0), and periodicity for θ - per enforces: The constitutive equation for the translations leads to:

h G ∼ • ∆ h θ -= h m - b 2 b (B.6) where ∆ h θ -= h θ -(0) -h θ -(b). Since θ - per is periodic, ∆ h θ -is
h u --h u - 0 = h E ∼ -1 • h r - 0 s 1 -e -1 × h G ∼ -1 • -e -1 × h r - 0 s 1 3 - b 2
s 2 1 2 + h m - b 2 s 2 1 2 + h θ - 0 s 1 (B.8)
where h u - 0 = h u -(0). The periodicity of h u - per enforces:

∆ h u -+ b 2 e -1 × ∆ h θ -= h k ∼ • h r - 0 -b e -1 × h θ - 0 (B.9)
where

h k ∼ =    b ES 0 0 0 b GSs,2 + b 3 12EI3 0 0 0 b GSs,3 + b 3 12EI2    h (B.10) and ∆ h u -=   be 11 be 12 -b 2 /2χ 11   h (B.11)
The first line of Equation B.9 leads directly to the traction in beam "h":

h r 0 1 = ESe 11 (B.12)
Similar analysis for beam "v" with:

∆ v u -=   be 22 -be 12 -b 2 /2χ 22   v and ∆ v θ -=   bχ 12 bχ 22 0   v (B.13) leads to: v m - b 2 =   GJ χ 12 EI 2 χ 22 0   v (B.14) and v r 0 1 = ESe 22 (B.15)
Displacement continuity at node. There remain h r 0 2 =v r 0 2 components to derive. Rotation field continuity enforces: h θ - 0 = v θ - 0 . Using this with Equation B.9 and its counterpart for beam "v" leads to:

U BG 3 = - k 0 κ 4 1 EI 2 + π 2 L 2 GS s,3
sin κx 1 (D.4)

Macroscopic stress. Furthermore, the solution for the bending moment and the bending gradient from [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF] are:

M ∼ =   M 11 M 22 √ 2M 12   = -1 g ∼ -1 • g - p 0 κ 2 sin κx 1 (D.5) R ⌢ =         R 111 R 221 √ 2R 121 R 112 R 222 √ 2R 122         =       -1 g ∼ -1 • g - 0 0 0       p 0 κ cos κx 1 (D.6)
It becomes:

M ∼ =   -1 0 0   k 0 bκ 2 sin κx 1 and R ⌢ =         -1 0 0 0 0 0         k 0 bκ cos κx 1 (D.7)
Localization. We introduce these values in the localization Equations 27 and 34 and get the following local generalized stress fields:

h m BG 2 = - k 0 κ 2 sin κx 1 (D.8) h r BG 3 = - k 0 κ cos κx 1 (D.9)
Reissner-Mindlin approximation. The Reissner-Mindlin approximation is easily derived while repeating the previous steps and setting f

⌢ ⌢ = b GSs,3 i ∼ ∼ • i ∼ ∼
. This leads to exactly the same result as the Bending-Gradient approximation.

Kirchhoff-Love approximation. The Kirchhoff-Love approximation is derived from the Bending-Gradient approximation while setting f ⌢ ⌢ = 0. The approximate deflection becomes:

U KL 3 = - k 0 κ 4 EI 2 sin κx 1 (D.10)
Furthermore, there is no localization related to shear forces and then:

Bending-Gradient approximation

The approximation for the lattice at 45 is derived following a similar procedure as in the previous section. However one have to take into account the rotation of the lattice. This affects both the constitutive equation of the Bending-Gradient model but also the localization in the unit-cell. Full details about rotation of R ⌢ and its constitutive equation are provided in [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF].

Macroscopic fields. The flexural stiffness rotated by 45 is now: Inserting this in Equation D.1 leads to the following macroscopic deflection:

d ∼ ∼ = b 2EI 2 1 + α α   1 -1-α 1+α 0 -1-α
U 3 = -U KL 1 + γ 2 α (1 + α) 2 η + α 2 1 + γ 2
sin κx 1 (D. [START_REF] Lebée | Homogenization of cellular sandwich panels[END_REF] where U KL = 4k0 κ 4 EI2(1+α) is the Kirchhoff-Love maximum deflection and γ 2 = κ 2 b 2 24 is the ratio between the size of the unit-cell and the span.

The bending moment writes as:

M ∼ = -    1 1 1+γ 2 1-α 1+α + γ 2 0    2k 0 bκ 2 sin κx 1 (D.16)
and the bending gradient:

R ⌢ = -          1 1 1+γ 2 1-α 1+α + γ 2 0 0 0 0          h 2k 0 bκ cos κx 1 (D.

Figure 1 :

 1 Figure 1: A square beam lattice

Figure 2 :

 2 Figure 2: An example of unit-cell of a space frame with four pairs of periodicity nodes and one interior node.

Figure 3 :

 3 Figure 3: The beam lattice unit-cell

3. 3 .

 3 Bending-Gradient and Reissner-Mindlin auxiliary problems 3.3.1. Bending-Gradient localization Once the Kirchhoff-Love auxiliary problem is solved, it is possible to derive the distributed loads of the Bending-Gradient auxiliary problem, k - (R) and l -(R) from Equation

Figure 4 :

 4 Figure 4: Beam lattice configuration at 0 in cylindrical bending configuration

Figure 5 :Figure 6 :

 56 Figure 5: Beam lattice configuration at 45 in cylindrical bending configuration

Figure 7 :

 7 Figure 7: Macroscopic deflection error along the span for a 4 unit-cell lattice (only half the span is plotted)

Figure 8 :

 8 Figure 8: Bending moment m 2 along Beams "+" for a 4 unit-cell lattice (only half the span is plotted)

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Torsion m 1 along Beams "+" for a 4 unit-cell lattice (only half the span is plotted)

  fully determined by Θ - KL defined in Equation13and we have:∆ h θ -

+ b 3 12EI3(B. 16 ) 2 ) 3 )

 31623 Result. The preceding results are summarized in the following localization field: Love stiffness tensors are derived using Equation 17 which leads to the result in Equation (28). Finally, using the constitutive equations: N ∼ = A ∼ ∼ : e ∼ and M ∼ = D ∼ ∼ : χ ∼ leads to the field localization given in Equation (27).Macroscopic deflection. The macroscopic deflection from[START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF] is:The matrix d ∼ ∼ is Kirchhoff-Love flexural stiffness and the 3 × 3 matrix f∼ is the upper-left 3 × 3-block of f ⌢ ⌢.For the configuration at 0 , d ∼ ∼ is given by Equation 28 and fWith these values for d ∼ ∼ and f ∼ , we have:

κ 2 sin κx 1 (D. 12 )

 112 Appendix D.2. Lattice at 45 Appendix D.2.

17

 17

  )

Two orthogonal direction of orthotropy and invariance under 90 in-plane rotation

Vectors and higher-order tensors are boldfaced and different underlinings are used for each order: vectors are straight underlined, u -. Second order tensors are underlined with a tilde: M ∼ . Third order tensors are underlined with a parenthesis: R ⌢ . Fourth order tensors are are doubly underlined with a tilde: D ∼ ∼. Sixth order tensors are doubly underlined with a parenthesis: f ⌢ ⌢ .

f αβγδǫζ follows major symmetry: f αβγδǫζ = f ζǫδγβα and minor symmetry f αβγδǫζ = f βαγδǫζ . Thus there are only 21 independent components.

Appendix A. Compatible fields for the Bending-Gradient theory

Equilibrium equations and boundary conditions involving stress fields are gathered in the set of statically compatible fields:

where the derivation operator ∇ -is also formally represented as a vector: a ∼ • ∇ -= a αβ ∇ β = a αβ,β is the divergence and a ∼ ⊗∇ -= a αβ ∇ γ = a αβ,γ is the gradient. ∂ω s is the portion of edge on which static boundary conditions apply and n -the related outer normal:

is the force per unit length and M ∼ d the full bending moment enforced on the edge. The plate is loaded exclusively with the out-of-plane distributed force p -= p 3 e -3 . This set of equations is almost identical to Reissner-Mindlin equations where shear forces have been replaced by the bending gradient R ⌢ .

The compatibility conditions and boundary conditions for strain fields are gathered in the set of kinematically compatible fields:

where U -= (U i ) is the average through the thickness of the plate 3D displacement and Φ ⌢ = (Φ αβγ ) is the generalized rotation. Γ ⌢ and Φ ⌢ are 2D-third-order tensors with the following symmetry: Φ αβγ = Φ βαγ . Moreover, ∂ω k is the portion of edge on which kinematic boundary conditions apply:

is a symmetric second-order tensor related to a forced rotation on the edge. These fields are almost identical to Reissner-Mindlin kinematically compatible fields where the rotation vector is replaced by the generalized rotation

, where ϕ -is a vector representing rotations leads to Reissner-Mindlin-like kinematics:

Thus in the general case, Reissner-Mindlin kinematics can be interpreted as a restriction of

Appendix B. Computation details for Kirchhoff-Love homogenization

In this section we detail the application of the Kirchhoff-Love auxiliary problem to the unit-cell of the lattice given in Figure 3. For convenience, the "KL" superscript for the solution is dropped in this section.

Node equilibrium. The equilibrium of node P 1 leads to:

where ∆m -= m -(s 1 = 0)m -(s 1 = b) and the left superscript denotes the beam (see Figure 3).

Appendix C. Exact solution for the cylindrical bending of a 45 lattice

Effect of symmetries. The lattice configuration is symmetric with respect to (A 2n , e -2 , e -3 ) plane and any (A j , e -1 , e -3 ) plane. Thus, the results are derived for half the lattice (see Figure 6) and we have the following relation between fields in each set:

where + denotes the stress in the first beam set and-denotes the stress in the second beam set.

is the symmetry matrix expressed in the main reference frame (m -and θ -are pseudo vectors, which explains the sign). Moreover the symmetry vanishes some displacement at nodes A j . We have in the main reference frame u j 2 = θ j 1 = θ j 3 = 0 where u j 2 means u 2 (x j ) 22 Equilibrium at nodes. Cutting the lattice at node A j with (A j , e -2 , e -3 ) plane leads to the following equilibrium equations at nodes (see Figure C.12):

where + r - j = + r -(x j ) is the resultant on beam "+" at x j . r - is the resultant through cutting plane. Taking into account symmetry, this leads to:

where + r j i and -r j i are the components of + r - j and -r - j taken in the local reference frame of the beam and r j i are the components of r - j taken in the main reference frame.

Macroscopic equilibrium.

Using global equilibrium at node A j , is possible to derive resultants and moments r - j and m - j through the cutting plane. We have:

in the main reference frame, where cos j = cos κx j .

Individual beam problem to solve. Combining preceding equations, it is possible to derive boundary condition to apply to each beam of set "+" in the local reference frame:

To this must be added + u 0 3 = 0 for the simple support and + θ 2n 2 = 0 for symmetry condition at mid span. With these boundary conditions and the beam model introduced in Section 2.2.1 one gets the result given in Section 4.2.1 Appendix D.1. Lattice at 0

The cylindrical bending of a Bending-Gradient plate loaded with p -= p 0 sin κx 1 e -3 is solved in details in [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF].

Localization of the solution. In order to relocalize macroscopic fields, we rotate them in the reference frame of beams "+" (rotation of +45 , cf. [START_REF] Lebée | A Bending-Gradient model for thick plates, Part II: Closed-form solutions for cylindrical bending of laminates[END_REF]). We get:

We introduce these values in the localization Equations 27 and 34 and get the local generalized stress fields for beams "+" given in Equations 54 and 55.

Appendix D.2.2. Reissner-Mindlin approximation Macroscopic fields. Again, the Reissner-Mindlin approximation is derived taking f

in Equation D.14 and following previous steps. We have the macroscopic deflection:

The macroscopic bending moment is:

Shear forces are derived taking the divergence of this:

Localization. In order to relocalize macroscopic fields, we rotate them in the reference frame of beams "+" (rotation of +45 ).


Mohamad Hairol Jabbar

Dominique Houzet
email: dominique.houzet@grenoble-inp.fr

Omar Hammami
email: omar.hammami@ensta-paristech.fr

Impact of 3D IC on NoC Topologies: A Wire Delay Consideration

Keywords: 3D NoC Architecture, Exploration, Network on Chip, Partitioning, Physical design

In this paper, we perform an exploration of 3D NoC architectures through physical design implementation based on two tiers Tezzaron 3D technology. The 3D NoC partitioning is done by dividing the NoC's datapath component into two blocks placed in the two tiers. Two Stacked NoC architectures namely Stacked 3D-Mesh NoC and Stacked 2D-Hexagonal NoC developed based on this partitioning strategy are analyzed by comparing their performances with Stacked 2D-Mesh NoC and classical 2D-Mesh and 3D-Mesh NoC. In order to measure the impact of wire delay on performance, two technology libraries (130 nm and 45 nm) representing old and advanced technologies have been used for the performance analysis. Results from physical implementations show that in advanced technologies such as 45 nm and below, the performance of Stacked 2D NoC topologies with datapath partitioning method have better performances compared with traditional 2D/3D Mesh topologies and Stacked 3D Mesh topology. We advocate here that with stacking there is no need for 3D NoC topologies for advanced 2-tier 3D IC and this is also confirmed for multistage networks like butterfly.

I. INTRODUCTION

As moving to sub-20 nm CMOS technology poses great design and manufacturing challenges, 3D integration [START_REF] Knickerbocker | Three-dimensional silicon integration[END_REF] is increasingly seen as a solution to those challenges for designing complex system on chip [START_REF]ITRS Report Interconnect[END_REF]. Global interconnect wire is one of the primary concern for advanced process technology (65 nm and below) (Figure 1) that has substantial contribution to the wire delay as well as power consumption even with the repeaters. By stacking dies or wafers, the performance can be increased due to reduction in interconnect length and so is power consumption due to the reduction in number of repeaters along the wires.

Performance improvement of 3D integration can be more prominent when compared with shrinking transistor technology [START_REF] Nomura | Performance Analysis of 3D-IC for Multi-Core Processors in sub-65nm CMOS technologies[END_REF]. Implementing the 2D design in 3D architecture on the same process technology could provides higher performance benefit than CMOS migration to the next process technology. Stacking multiple dies will also reduce the total footprint of a chip making it very suitable for mobile devices. However, several challenges such as thermal, yield, cost, design tools and testing of 3D architecture must be overcome before 3D IC technology can be widely adopted as a mainstream technology [START_REF] Li | Is 3D integration an opportunity or just a hype?[END_REF].

This paper presents an exploration of 3D NoC architectures through physical design case studies. Our motivation is that we want to explore different partitioning strategies from the previous reported works for 3D NoC architecture and then evaluate their performance accurately based on layout-level routed netlist. The contributions of this work are as follows: II. RELATED WORKS Many issues in 2D NoC architecture and design have been studied over the past years covering various aspects such as design flow, implementation evaluation and design space exploration. However, research in 3D NoC is still new and many issues remain unexplored especially in real design and implementation. Design space exploration of 3D NoC topologies through cycle accurate simulation have been performed showing the benefits of 3D design in terms of throughput, latency and energy dissipation for mesh-based and tree-based NoC architecture [START_REF] Feero | Networks-on-Chip in a Three-Dimensional Environment: A Performance Evaluation[END_REF]. In [START_REF] Pavlidis | 3-D Topologies for Networks-on-Chip[END_REF], zero load latency and power consumption analytical models of various 3D NoC topologies have been evaluated proving the advantages of combining 3D IC with 3D NoC architecture. We base upon this literature to investigate further the results by doing analysis from physical design implementation results. Another work [START_REF] Kim | A novel dimensionally-decomposed router for on-chip communication in 3D architectures[END_REF] proposed a novel 3D router architecture by decomposing the router into different dimensions to provide better performance over other 3D NoC architectures. We differ from the previous reported works as we focus on partitioning 3D NoC architectures and evaluate their performance through layout-level netlist for more accurate analysis of wirelength, timing, area and power consumption.

Several experiments have been conducted investigating the performance of 3D architectures based on the results from physical design implementations. Work in [START_REF] Lee | Timing analysis and optimization for 3D stacked multi-core microprocessors[END_REF] has studied different partitioning styles for implementing 3D multicore architectures namely core level, block level and gate level showing that TSV capacitance, EDA tools and timing optimization methods have strong impact on the performance of the final 3D architecture. In [START_REF] Pathak | Through-siliconvia management during 3D physical design: When to add and how many?[END_REF], they showed that 3D architecture could lose or reduce its benefit due to the tools inability to perform 3D-aware optimization. On the other hand, larger circuits tend to gain more improvement from 3D architecture over its 2D counterpart for advanced technology such as 45 nm node. In [START_REF] Thorolfsson | Logic-onlogic 3D integration and placement[END_REF], the study of different 3D placement methods on the performance of three 3D architectures showed that true-3D placement method produces the highest performance improvement over other methods at old technology (130 nm) indicating the importance of 3D-aware tools to obtain maximum benefits of 3D integration. However, no previous work has been presented with detailed performance evaluation on various physical design metrics (wirelength, timing, impact of wire length) of 3D NoC architecture in particular with 3D Meshbased NoC architecture.

III. 3D IC TECHNOLOGY

The 3D integration technology we used is based on Tezzaron [START_REF] Patti | Homogeneous 3D Integration[END_REF] that uses TSV for peripheral IOs and microbumps for inter die connections. The two-tier 3D stacking method is based on wafer-to-wafer bonding, face-to-face method with via-first approach as illustrated in Figure 2. The inter-die microbumps provide high interconnection density up to 40,000 per mm 2 without interfering to FEOL (front-end-ofline) device or routing layers. It is also possible to implement four tiers by stacking through back-to-back using TSV of the two face-to-face stacking in order to have higher design complexity but it will not be covered in this paper. In order to analyze performance of 3D NoC architectures in advanced technology, we have chosen 45 nm standard library from ST Microelectronic [START_REF] Microelectronic | 45 nm ST Microelectronic standard library[END_REF]. We use similar 3D structure for inter-tier connections using microbumps as in Tezzaron technology but we replace the 130 nm technology of Global Foundries with 45 nm ST Microelectronic standard library. The 45 nm technology used in this study has seven metal layers where metal seven is used for bonding and the routing is limited until metal six.

A. Design Flow Based on 2D EDA Tools

The 3D design flow is developed based on the 2D EDA tools. This flow is made possible with the Tezzaron 3D technology using mirobumps for inter-tier connections. These microbumps have negligible delay for the inter-tier connection and thus we can perform 3D timing analysis at post-synthesis stage without any inaccurate delay estimation of inter-tier connection. Post-synthesis static timing analysis (STA) for each tier is done separately before 3D timing analysis is performed. In order to perform 3D timing analysis at post-synthesis stage, we create a top level netlist that instantiate both tiers and connect them using inter-tier wires that represent microbumps. Using the generated timing constraints, timing optimization is carried out using 2D place and route tool for each tier separately. For post-route 3D timing analysis, we create the top level netlist as in the postsynthesis step and feed the SPEF file of each tier into Synopsys PrimeTime for timing and power analysis. The parasitics for the microbumps are ignored due to their negligible delay.

IV. BASELINE NOC ARCHITECTURE A. Router and Network Interface Architecture

The router and network interface architectures are standard architectures used for mesh NoC. In this experiment, we did not include the processor core in order to make the experiment easier. Also, if we include the processor cores, the results will be more significant (the benefit of the partition will be higher) because of the increased inter-router links.

B. Baseline 3D Mesh NoC

In this architecture, the 3D NoC is implemented on two tiers where each tier has identical blocks as shown in Figure 3. This is the straightforward extension of 2D Mesh NoC architecture to 3D Mesh NoC where we just take a copy of a tile (a router and a network interface) and put it on top of each tile. Compared with the area of 2D Stacked Mesh NoC, this architecture has slightly more area due to the additional ports for vertical connections. This 4x2x2 Mesh NoC architecture is based on a 3D router architecture that has vertical links for inter-tier connections between routers. These physical vertical links shown in red color are based on the logical vertical links in each 3D router.

V. EXPLORATION OF 3D NOC ARCHITECTURES A. 3D NoC Partitioning

In this section, we describe the partitioning method to be used for the next 3D stacked NoC architectures. The FIFO buffer is dominating silicon area in the NoC architecture. Thus, it is a good approach to partition it into two tiers. Other datapath components are also partitioned into two tiers at bitlevel. For example, for the 32 bit FIFO size, the resulting implementation will be 16 bits per tier. For the non-datapath components such as routing logic, arbitration logic and FIFO control, we place them on each tier by trying to balance the area of both tiers. Figure 4 illustrates this partitioning method with respects to 2D and baseline 3D Mesh NoC architectures. Rather than using automatic tools such as HMetis to partition the design, we focus on dividing the datapath manually into two parts and place them into two tiers in order to preserve the homogeneous properties of tile block architecture. Another reason for not using this automatic tool is because the tool also tries to optimize the nets between gates in the netlist with no capability of 3D placement meaning that logic cells can be interchangeably partitioned into the two tiers which will eventually affect the 3D timing path. The first architecture designed with the partitioning method is depicted in Figure 4. Here, rather than stacking the tiles on top of each other, we map the 3D NoC on the 2D layout and then partition it into two tiers. As shown in Figure 6(a), the green links represent logical vertical connections between 3D routers while the physical vertical links in orange color are basically the 2D logical links. By doing this, the area is slightly increased compared with the 3D Mesh NoC but reduced compared with the 2D one. However, this partitioning method requires higher number of inter-tier connections than pure 3D Mesh NoC. One disadvantage of this structure is that the inter-router wire links are not equal between all routers because vertical wire links are longer than other links. 8 (a). We used the here the same datapath partitioning method (cf. section 5.1). Previous work [START_REF] Malony | Regular processor arrays[END_REF] had proved that hexagonal topology is the most efficient topology and theoretical exploration of addressing, routing and broadcasting in hexagonal mesh architecture has also been explored.

1) Packet Routing

Routing is illustrated in Figure 5. Basically the packets will be first routed through X direction and then to Y direction to reach the destination. However, in the case of a router with a diagonal link in the same direction as Y direction, the packets will be routed through this diagonal link instead of X axis link. Therefore, from the diagram, the packets will be routed from router 00 to 33 through router 11, 12 and 23. The diameter of the Hexagonal NoC can be formulated as d = (x-1) + (y-1) -(x/2), where x is the number of hops in X axis and y is the number of hops in Y axis. This hexagonal routing is a dimension ordered routing which is thus deadlock free.

2) Physical Implementation

The physical size of the tiles is determined by measuring the distance between tiles such that the distance between the six neighboring tiles is equal. This is to make sure that this square floorplan area is identical to the original hexagonal shape. Although it is possible to create hexagonal floorplan in SoC Encounter, we choose to adopt rectangular floorplan. As shown in Figure 7, the rectangular floorplan for the hexagonal architecture can be carefully arranged such that the interrouter links are equal and use only vertical and horizontal links, and thus we avoid the use of diagonal wires in the case of an hexagonal shape having four diagonal edges. We adopted the equation of (a/2) 2 + b 2 = c 2 , where a is the tile's height, b is tile's width and c is the physical direct distance between the two tiles to determine the size of each tile. We first fix the value of a and then find the value of b such that c is equal to a at the same time meeting initial target utilization. We also have derived mathematical formulation proving that the surface area of the square floorplan is identical to the original hexagonal structure. Let's say a equal to 579 µ m, following the equation above will obtain the value of b equal to 500 µm and c equal to 577 µm with the initial target utilization of 60%. To compare the diameter for both topologies, consider an example of a 4 x 5 network. The diameter for 3D Mesh NoC is 6 while for 3D Hexagonal NoC is 6 but Hexagonal NoC has shorter inter-router wire links benefit.

VI. EXPERIMENTAL RESULTS

For older technology such as 130 nm and above, wire length effect is not significant and the delay in the critical paths is mostly determined by the delay of the gates. As shown in Table I and Figure 9, the 3D NoC architectures do not benefit in terms of speed and power consumption. The power consumption is even higher in 3D architectures due to the additional gates as well as the increased wirelength. In this study, we used simple partitioning method to partition the 2D design into 2 tiers. However, some studies have shown that automatic partitioning tools could provide performance improvement over 2D architecture even using old technology such as 130 nm and 180 nm [START_REF] Thorolfsson | Reconfigurable five-layer three-dimensional integrated memoryon-logic synthetic aperture radar processor[END_REF]. Partitioning is very important in 3D design primarily for old technology. Using automatic partitioning tool such as hMetis [START_REF] Karypis | Multilevel hypergraph partitioning: applications in VLSI domain[END_REF] helps to improve the performance of 3D architecture although it is still not significant because the tool try to optimize the connections between gates in the synthesized netlist but is not able to perform in-place 3D optimization during place and route as in usual 2D optimization. At 45 nm, automatic partitioning tools can provide higher performance improvement for the 3D architecture than for old technology.

A. Analysis on the Impact of Wire Delay

As for wire delay, older technology nodes (such as 130 nm) do not have significant wire delay on the performance. The critical path for all designs in this study is located within the tile block (from bottom tier to top tier in 3D stacked architecture) except for 3D Mesh NoC architecture where its critical path is between two routers. Looking at the 3D critical paths for all 3D NoC architectures in this study, the ratio of wire delay is about 3% of the total critical path delay. For comparison, the wire delay in 2D architecture using the same process technology is about 5.7% of the total critical path delay and thus we can generally conclude that 3D architecture in this technology will not offer any benefit in terms of speed.

However, there is still an opportunity to gain benefit from 3D architecture by optimizing partitioning method as demonstrated by several works previously using older technology nodes [START_REF] Thorolfsson | Logic-onlogic 3D integration and placement[END_REF] although the results is not very significant compared with ideal improvement we should get. Additionally, analyzing the critical path delay for the 2D architecture using 45 nm technology indicates that wire delay is about 1% due to very small area. We expect to see larger portion of wire delay in the critical path for 2D architecture with larger design.

B. Extrapolation for Advanced Technology

For designs using 45 nm used in this study, the 3D architectures still do not provide any improvement over its 2D design as shown in Table II and Figure 10. However, it shows a reduction trend of the gap between 3D and 2D architectures compared with the results in Figure 9 using 130 nm technology. If we look at the area, we can see that this design consumes very small area (less than 1 mm 2 ) and this is the primary reason why there is no improvement obtained using 45 nm technology. Previous work have demonstrated for large designs (about 36 mm 2 in 2D architecture), substantial performance improvement (75% reduction in longest path delay) that could be achieved over 2D architecture using the same 45 nm technology because wirelength becomes significant [START_REF] Pathak | Through-siliconvia management during 3D physical design: When to add and how many?[END_REF]. Table III shows the extrapolation of wire delay for 22 nm technology based on the critical path wire delay in 45 nm and on the data from ITRS 2007 interconnect report for global wire without drivers (1.02ns (45nm), 3.3ns (22nm), 5.9ns (16nm)). This extrapolation is intended to show that when the design used is realistically large, we will see improvement for the proposed hexagonal NoC topology in stacked 3D architecture compared to the other solutions. The gate delay value for 22 nm is assumed to improve two times over 45 nm technology because it is two technology generations from 45 nm and the tile area (and thus the interrouter wire length) is assumes to be 3 mm x 3 mm for the 3D Mesh NoC considering the area of commercial grade LEON3 processor [START_REF] Pathak | Through-siliconvia management during 3D physical design: When to add and how many?[END_REF]. From the 3 mm inter-router wirelength of 3D Mesh NoC, we calculate the wire length for 2D Stacked Mesh and 3D Mesh as follows:

From the rectangular area equation for the hexagonal floorplan, where a is the inter-router length for 3D Stacked Hexagonal NoC, x and y is the new inter-router length for 2D Stacked Mesh and 3D Mesh respectively. The wire length for 3D Mesh and 3D Stacked Mesh is equal because 3D Stacked Mesh has half the area of 3D Mesh but has double the interrouter length for 3D logical vertical links due to the 2D mapping of the 3D Stacked Mesh. This extrapolation is simplified by ignoring the router area impact due to the different number of ports for different topologies which is 4, 5 and 6 ports for 2D Mesh, 3D Mesh (also 3D Stacked Mesh) and 2D Hexagonal respectively. As can be seen from the table, the wire delay is becoming more significant for 16 nm technology and thus it will have strong impact on the critical path delay especially for 3D Mesh NoC and 3D Stacked Mesh NoC (because of logical vertical links between routers) since it has longer inter-router wire links. The 2D Stacked Mesh NoC outperforms the 3D Mesh NoC and 3D Stacked Mesh NoC in the wire delay and eventually the total delay. However, the 3D Stacked Hexagonal NoC is shown to have better improvement than 2D Stacked Mesh in terms of network latency because it has a lower diameter compared with 2D Stacked Mesh and thus it will benefit for applications running on large networks.

C. Extending for Other NoC Topology

The results from this experiment can also be applied to other topologies such as multi-stage interconnection and butterfly network. Previous work has proposed to partition the butterfly network by folding it into several tiers [START_REF] Hammami | Design of 3D-IC for butterfly NOC based 64 PE-multicore: Analysis and design space exploration[END_REF]. As has been explained in previous sections for the hexagonal topology, the same method can also be used for butterfly topology, meaning that partitioning it into two tiers (which can be referred to 3D Stacked Butterfly) could improve its performance due to the shorten wire links between stages. VII. CONCLUSIONS This paper has presented an exploration of 3D NoC architectures through physical design case studies. The performance of these 3D NoC architectures, namely 3D Stacked Mesh NoC and 3D Stacked Hexagonal NoC has been analyzed by comparing with the 2D Mesh NoC and also traditional 3D Mesh NoC architecture. For advanced technologies such as 45 nm and beyond, the 3D NoC architectures based on this partitioning method show better performance than traditional 3D Mesh NoC architecture due to significance of wire delay effect. Future work will extend this study to manycore MPSoC and 3D EDA tools implementation.

Figure 1 :

 1 Figure 1: Global interconnect and gate delay trend for different process technologies from ITRS report • Propose a partitioning strategy to optimize NoC architecture performance based on stacking router's datapath components into two blocks on two tiers. • Study the wire delay impact on the 3D NoC architectures by performing physical implementation using two standard cell libraries, 130 nm and 45 nm.

Figure 2 :

 2 Figure 2: Cross section of Tezzaron 3D IC technology with corresponding parameters

Figure 3 :

 3 Figure 3: 3D Mesh NoC design a) block diagram b) floorplan c) routed layout

Figure 4 :

 4 Figure 4: Partitioning method for the 3D stacked NoC architecture C. 3DNoC2: 3D Stacked Hexagonal NoC Due to unequal inter-router wire link lengths in the 3D Stacked Mesh NoC architecture because of the logical vertical links (green lines in Figure 6(a)), we proposed a new topology having same length of inter-router physical links called hexagonal topology shown in Figure8(a). We used the here the same datapath partitioning method (cf. section 5.1). Previous work[START_REF] Malony | Regular processor arrays[END_REF] had proved that hexagonal topology is the most efficient topology and theoretical exploration of addressing, routing and broadcasting in hexagonal mesh architecture has also been explored.

Figure 5 :

 5 Figure 5: Hexagonal routing block diagram

Figure 7 :Figure 8 :

 78 Figure 7: Rectangular floorplan for hexagonal architecture showing the inter-router links



  And c = a (because equal inter-router links), thus,

Figure 10 :

 10 Figure 10: Performance comparison of 3D NoC architectures over 2D NoC in 45 nm technology

Table I :

 I Physical implementation results for 3D NoC architectures in 130 nm technology

	Parameters	3D Mesh NoC	3D Stacked Mesh NoC	3D Stacked Hexagonal NoC
	Vertical connections	1763	6261	7255
	Core area (mm 2 )	3.24	4.37	5.43
	Total wirelength (m)	12.48	14.01	17.03
	Longest path (ns)	4.20	4.60	4.73
	Power @ 333 MHz (W)	1.44	1.25	1.40
	Figure 9: Performance comparison of 3D NoC architectures
	over 2D NoC in 130 nm technology

Table II :

 II Physical implementation results of 3D NoC architectures in 45 nm technology

	Parameters	3D Mesh NoC	3D Stacked Mesh NoC	3D Stacked Hexagonal NoC
	Vertical connections	1763	6261	7255
	Core area (mm 2 )	0.79	0.91	1.01
	Total wirelength (m)	5.5	5.9	6.5
	Longest path (ns)	3.23	3.33	3.59
	Power @ 333 MHz (W)	0.23	0.24	0.26

Table III :

 III Extrapolation of delay from physical implementation result for 3D NoC architectures performance

	Technology / NoC architecture	Gate delay (ns)	Wire delay (ns)	Total delay (ns)	Diameter (8x9 network)	Max latency (ns)
		2D Stacked Mesh NoC	2.6	1.5	4.1	15	61.5
		3D Mesh NoC	2.6	3.0	5.6	12	67.2
	45 nm	3D Stacked Mesh NoC	2.6	3.0	5.6	12	67.2
		3D Stacked
		Hexagonal	2.6	1.59	4.19	11	46.1
		NoC
		2D Stacked Mesh NoC	1.3	4.95	7.55	15	113.25
		3D Mesh NoC	1.3	9.9	11.2	12	134.4
	22 nm	3D Stacked Mesh NoC	1.3	9.9	11.2	12	134.4
		3D Stacked
		Hexagonal	1.3	5.25	6.55	11	72.05
		NoC
		2D Stacked Mesh NoC	0.6	8.85	9.45	15	141.75
		3D Mesh NoC	0.6	17.7	18.3	12	219.6
	16 nm	3D Stacked Mesh NoC	0.6	17.7	18.3	12	219.6
		3D Stacked
		Hexagonal	0.6	9.38	9.98	11	109.78
		NoC