domain: phys spi

text:
Enrico Stragiotti
email: enrico.stragiotti@onera.fr

Franc ¸ois-Xavier Irisarri

Cédric Julien

Joseph Morlier
email: joseph.morlier@isae-supaero.fr

Enhanced truss topology optimization (TTO) applied to a cellular wing box

Cellular architectured structures' most appealing characteristic is that they are able to create large structures assembling small repetitive components. Thanks to their modular nature, they bring interesting features among which reduced tooling, fast assembly, and short repair time. In the first part of the article, we formulate an optimizing method that minimizes the mass of a cellular architectured structure taking into account internal stresses, local buckling, and pattern repetition constraints. The proposed solving algorithm mitigates the appearance of local minima, a critical problem for discrete trusses. In the second part, the optimization method is applied to a real-size aeronautic application: the minimization of the mass of a cellular 3D wing box subject to lift and torsion loads. Compared to classic cell topologies, the proposed method found a cell 10-times lighter, at the cost of increased manufacturing difficulty.

I. Introduction

Assembled cellular ultralight structures present an opportunity to greatly improve the performance and cost efficiency of modern aerostructures [START_REF] Cramer | Elastic shape morphing of ultralight structures by programmable assembly[END_REF] . The repetitive nature brings various interesting features among which reduced tooling, fast assembly, and short repair time. Additionally, as the mechanical performance of the structure is greatly influenced by the topology and the materials of the repetitive pattern, cellular structures are naturally prone to optimization for various scenarios.

In the structure optimization field periodic materials are often modeled by means of asymptotic homogenization [START_REF] Zhou | Design of graded two-phase microstructures for tailored elasticity gradients[END_REF] . The heterogeneous cell topology (or Representative Volume Element (RVE)) is treated as homogeneous material with associated mechanical properties i.e., equivalent elastic tensor, shear modulus, etc. The homogenization approach is valid only if the RVE contains enough information about the heterogeneous material and if the structure presents significant periodicity. When designing cellular structures (and not materials), no scale separation is assumed between the repetitive pattern and the structure itself. Consequently, the hypotheses of asymptotic homogenization are not always verified. For that reason, full-scale approaches [START_REF] Wu | Topology optimization of multi-scale structures: a review[END_REF] have been developed. Among all the structure optimization methods, topology optimization allows achieving any shape within the design space instead of dealing with predefined configurations. Topology optimization considers the design domain as a continuum, in which each location may or may not have a material assigned to it. Even if the freedom of the design space offered by the continuum meshes is high, it is known that at very low volume fractions (ultralight structures) and, especially if buckling constraints and manufacturing consideration (minimum plate thickness) are taken into account, the optimal topology is a truss structure [START_REF] Sigmund | On the (non-)optimality of Michell structures[END_REF] . In addition to that, truss structure design naturally relies on constraints on maximum allowable stress, buckling, and maximum slenderness, which are all known for being difficult to implement on topology optimization using a Nested Analysis and Design (NAND) formulation [START_REF] Rozvany | On design-dependent constraints and singular topologies[END_REF] . On top of that, topology optimization proved to be extremely computation intensive, especially when optimizing large-scale 3D structures [START_REF] Aage | Giga-voxel computational morphogenesis for structural design[END_REF] with low volume fractions, even with adaptive meshes [START_REF] Salazar De Troya | Adaptive mesh refinement in stress-constrained topology optimization[END_REF] .

Building on a Simultaneous Analysis and Design (SAND) formulation, Truss Topology Optimization (TTO) is a structural optimization algorithm that finds a near-optimal truss structure with respect to a given load case. Instead of working with a continuous domain as topology optimization does, TTO optimizes the cross-sectional areas and the connectivity of the members of a discrete ground structure. In the early works, the TTO problem was formulated in terms of member forces 8 with plastic material formulation, ignoring the kinematic compatibility to obtain a Linear Programming (LP) problem. Formulated in that way, the problem becomes inexpensive to solve but loses in generality. One can solve only single load case problems for elastic design and is restricted to statically determinate structures or mechanism [START_REF] Rozvany | Layout Optimization of Structures[END_REF] . Aiming at mitigating these criticalities, TTO has been enriched with the kinematic compatibility constraint [START_REF] Reinschmidt | Applications of linear programming in structural layout and optimization[END_REF] , and also local buckling [START_REF] Achtziger | Local stability of trusses in the context of topology optimization Part I: Exact modelling[END_REF] and global buckling [START_REF] Kočvara | On the modelling and solving of the truss design problem with global stability constraints[END_REF] constraints.

The objective of the present research is to develop an optimizing algorithm to minimize the mass of an assembled cellular structure. The proposed formulation builds on the classic LP TTO optimization problem, adding the aforementioned kinematic compatibility and local buckling constraints. Furthermore, we enrich the formulation by imposing pattern repetition, implemented using the full-scale approach called variable linking [START_REF] Zhang | Scale-related topology optimization of cellular materials and structures[END_REF] . The proposed formulation performance is evaluated on a simplified 3D wing box load case and benchmarked against classic cell topologies.

II. Cellular truss topology optimization

In this section, we propose a TTO formulation developed specifically to tackle the difficulties in optimizing 3D ultralight cellular structures with no scale separation between the repeating pattern and the structure. The objective we choose is the minimization of the mass constrained by stress, local buckling, and pattern repetition constraints. The optimized structure is achieved by modifying the cross-sectional area and the connectivity of a discrete mesh. TTO acts on a ground structure that is defined as a set of bar members that connect a grid of M pin-jointed points.

The proposed formulation is based on the SAND approach, where, in addition to the member cross-sectional areas, the member forces and the nodal displacements are used as state variables. Consequently, the equations of structural mechanics are imposed as constraints of the optimization problem and not explicitly solved as for NAND approaches. One of the advantages of SAND approach is that the state variables are independent of each other and, thus, the sensitivity calculation of the constraints functions is simpler.

The problem formulation is stated in terms of linked cross-sectional member areas κ κ κ, member forces q q q and nodal displacements U U U as follows:

P 1 : min κ κ κ,q q q,U U U V = ℓ ℓ ℓ * T a a a (Volume minimization) s.t. B B Bq q q = f f f (Force equilibrium) q q q = a a aE ℓ ℓ ℓ b b b T U U U (Compatibility constraints)
q q q ≥ -sa a a 2 ℓ ℓ ℓ 2 (Euler's buckling constraints)

a a a = G G G ⊗ κ κ κ (Variable linking) -σ C a a a ≤ q q q ≤ σ T a a a (Stress constraints) κ κ κ ≥ 0, (Positive cross-sectional member areas)

where B B B is a 2M × N el matrix containing the direction cosines of j-th member with respect to the i-th degree of freedom to calculate the nodal force equilibrium. M is the number of nodes and N el the number of members of the ground structure. q q q = [q 1 , q 2 , . . . , q N el ] T is the vector containing the internal member forces caused by the external load

f f f = [ f 1 , f 2 , . . . , f 2M ] T . The state variable κ κ κ = [κ 1 , κ 2 , .
. . , κ n el ] T represents the cross-sectional area of the n el members of the repetitive cell and is linked to the structure cross-sectional member areas a a a = [a 1 , a 2 , . . . , a N el ] T through the cellular mapping vector G G G. σ C and σ T are the compressive and tensile maximum allowable stresses of the material, respectively. b b b represents the i-th column of B B B, while s is a constant used for Euler buckling calculation that depend on the cross-section topology. Our optimization takes into account only the linear behavior of the structure. The resolution of Problem P 1 produces complex structures made up of a multitude of small members that resemble to Michell structures. Substituting in the volume evaluation function V = ℓ ℓ ℓ T a a a the member physical length ℓ ℓ ℓ with ℓ ℓ ℓ * = [ℓ 1 + j c max(ℓ ℓ ℓ), ℓ 2 + j c max(ℓ ℓ ℓ), . . . , ℓ N el + j c max(ℓ ℓ ℓ)] T , one would penalize the appearance of small members [START_REF] Parkes | Joints in optimum frameworks[END_REF] . ℓ ℓ ℓ * is called augmented member length and j c is the joint cost. This approach permits avoiding the apparition of structures with tiny features when a fine ground structure is adopted.

The contributions of our studies is to revisit and enrich the classic LP TTO formulation [START_REF] Dorn | Automatic design of optimal structures[END_REF] . The cellular nature of lattice structures is taken into account by implementing the full-scale method called variable linking. These constraints, together with local buckling and kinematic compatibility constraints are regrouped in a unique formulation, with the idea of minimizing the need for post-processing between the optimized structure and the manufacturing phase. The resulting formulation proved to be extremely prone to local minima, and therefore we suggest a two-step optimization algorithm to limit that negative behavior. All these additions, together with the resolution algorithm, will be detailed in the next sections.

A. Cell structure optimization

To take into account the cellular topology of the structure, the variable linking scheme is implemented on a structured grid. This approach is usually found in the topology optimization literature and is adapted here to work with TTO. The variable linking approach is a full-scale method with a repeated pattern i.e. no scale separation is assumed between the repetitive cell and the structure.

Let us take a 2D truss structure for convenience. It is made by M joints and partitioned into N c = N c,x × N c,y repetitive cells (see Figure 1). Every cell is composed by m = m x × m y joints connected by n el = (m(m -1)/2) elements in the case of a fully-connected cell ground structure. We define κ κ κ = [κ 1 , κ 2 , . . . , κ n el ] T as the linked cross-sectional area of the members defined in the design cell Ω c . The optimization is conducted only on the repetitive design cell Ω c , while the structural response and the sensitivity analysis are evaluated over the whole structure Ω.

F N c,x N c,y m x m y Ω c Ω N c,x Ω Figure 1:
The variable linking scheme is applied to a cellular 2D structure. The design cell Ω c represents the design space of the optimizer, and it is repeated over the structure domain Ω.

Variable κ κ κ is linked to the structure cross-sectional area a a a over the entire structure Ω using the following equation:

g g g cell := a a a = G G G ⊗ κ κ κ =    G 1 κ κ κ . . . G N c κ κ κ    , (2)
where ⊗ stands for the Kroneker product. In the case of a single cell optimization,

G G G = [G 1 , G 2 , . . . , G N c ] ∈ R N c is defined as an all-ones vector of size N c , G G G = [1, 1, . . . , 1] ∈ R N c
. The Kroneker product operation described above tiles the variable κ κ κ N c -times. The Kronecker product is a linear operator and suitable for the sensitivity analysis of problem P 1 .

B. Kinematic compatibility constraints

To optimize cellular structures, that are usually statically indeterminate, or to solve multiple load cases, it is necessary to impose the kinematic compatibility constraints [START_REF] Cheng | Some aspects of truss topology optimization[END_REF] . These equations restrict the displacement field U U U in such a way that strains and internal stresses comply with Hooke's law σ j = E j ε j with j ∈ [1, 2, . . . , N el ]. Recalling that in a truss the relationship between nodal displacements and member deformations is b b b T j U U U = ℓ j ε j , we can formulate the kinematic compatibility constraints g g g comp as:

g j, comp := a j E j ℓ j b b b T j U U U -q j = 0 ∀ j ∈ [1, 2, . . . , N el ], (3)
An advantage of using the kinematic compatibility constraints to impose the linear mechanical behavior with respect to solving a Finite Element Analysis (FEA) inside a NAND formulation is that the problem stays wellposed even if the section variables go to 0. As the linear system K K KU U U = f f f is never explicitly solved during the optimization, it is not necessary to impose a lower bound on these variables to avoid a singular stiffness matrix. Kinematic compatibility constraints are non-linear as the design variables q q q are dependent on a a a and U U U.

C. Local buckling constraints

Adding local buckling constraints is fundamental as ultralight truss structures are often dominated by this mode of failure [START_REF] Sigmund | On the (non-)optimality of Michell structures[END_REF] . Using the notation of Achtziger 11 , the Euler's buckling constraints g g g buck are stated as:

g j, buck := q j + sa 2 j ℓ 2 j ≥ 0 ∀ j ∈ [1, 2, . . . , N el ], (4)
where s is a parameter dependent of the member material and section topology. In this work we choose to use circular sections with s = πE/4. Equation 4 is non-linear but still convex, and it would not create too many difficulties for a modern non-linear optimizer. The criticality is that by imposing local buckling constraints over a topology optimization problem (where the lower bound for the member section is 0), the optimization domain becomes disjointed [START_REF] Cheng | Some aspects of truss topology optimization[END_REF] . This problem is made even worse as we introduce the non-linearities of the kinematic compatibility constraints.

D. Optimization algorithm

Directly solving Problem P 1 using a Non-Linear Programming (NLP) solver systematically yields unsatisfying local minima. Getting rid of a member once the local buckling constraints are active proved to be challenging even for modern NLP solvers. The solutions are, thus, usually made by many intersecting thin bars. This behavior becomes even worse when a high number of cells is used.

Because of that, we propose a two-step optimization algorithm that aims at mitigating the proneness of the optimizer to get stuck on local minima. The steps are:

1. Problem P 1 is relaxed: kinematic compatibility constraints are omitted and local buckling constraints are linearized. We call this relaxed problem P 2 . The solution to this Problem P 2 proved to be less influenced by local minima as the response surface become smoother. Subsequently, the initial ground structure is simplified, removing all the members that do not appear in the relaxed problem solution.

2. The kinematic compatibility and the local buckling constraints are restored. Problem P 1 is solved in its original form on the simplified ground structure. The initial values for the cross-sectional areas are the solutions of Problem P 2 .

The algorithm proposed to solve the ultralight cellular structure optimization is presented in Fig. 2. The proposed two-step optimization algorithm. In the first part, the solution of a relaxed problem is used to simplify the fully connected cell ground structure. In the second part, the complete problem with the restored constraints is solved on the simplified ground structure.

While there is no mathematical guarantee that the optimized point found following this method is the global optimum of the NLP problem, we systematically find better results using this approach instead of directly solving P 1 .

Let us break down the algorithm presented in Fig. 2 and explicate all the phases in detail. Taking as input the boundary conditions (physical dimensions) and the cell parameters (number of repetitive cells and number of nodes per dimension), the first step is to generate the initial ground structure. The optimal structure is searched Association for Structural and Multidisciplinary Optimization in the UK (ASMO-UK) in a subset of this discretization. This step is critical as the solution is highly dependent on the initial ground structure [START_REF] Kočvara | On the modelling and solving of the truss design problem with global stability constraints[END_REF] .

In our work, we consider a fully-connected ground structure for every cell of the structure. This promotes the exploration of new cell topologies. These distinct ground structures are then assembled based on the cell parameters of the structure. Therefore, in the final structure, each cell contains its own members, and there are no members spanning over multiple cells. Overlapped members on the cell edges are kept during the optimization and their area is summed at the end.

Once the ground structure is defined, it is possible to iteratively solve the relaxed problem P 2 with linearized buckling constraints. The local buckling formula is linearized for every j member in the neighborhood of a 0 j using a first order Taylor expansion [START_REF] Schwarz | Efficient size and shape optimization of truss structures subject to stress and local buckling constraints using sequential linear programming[END_REF] : qcr, j = q cr, j (a 0 j ) + (a ja 0 j )

∂ q cr,i ∂ a a=a 0 j ( 5
)
where q cr, j (a 0 j ) represents the Euler's critical load of the j-th bar with cross-sectional area a 0 j . The linearized Euler buckling constraints the Euler's buckling constraints g g g buck are stated as:

g j buck := q j ≥ qcr,i = - sa 0 j 2a j -a 0 j ℓ 2 j ∀ j ∈ [1, 2, . . . , N el ], (6)
where superscript ∼ indicates linearized variables.

We can now state the relaxed Problem P 2 in terms of linked cross-sectional member areas κ κ κ and member forces q q q as: P 2 : min κ κ κ,q q q V = ℓ ℓ ℓ * T a a a (Volume minimization)

s.t. B B Bq q q = f f f (Force equilibrium)
q q q ≥ -sa a a 0 2a a aa a a 0 ℓ ℓ ℓ 2 (Euler's buckling constraints)

a a a = G G G ⊗ κ κ κ (Variable linking) -σ C a a a ≤ q q q ≤ σ T a a a (Stress constraints) κ κ κ ≥ 0. (Positive sections) (7)
As the objective function and the constraints are all linear, the problem can be iteratively solved by the means of Sequential Linear Programming (SLP) algorithm until convergence. The logic of the SLP algorithm is presented in Fig. 3.

Relaxed NLP problem a 0 , q 0 Linearize buckling constraint g g g buck around point a 0 The gradient of Eq. 6 goes to zero when a member is removed from the topology i.e., member which crosssectional area tend to zero are heavily penalized by the optimizer. In that way, even when the topology evolves during the iterations of the SLP, removed members are hardly reintroduced in the optimization. To allow suppressed member to reappear in the optimized structure we apply a reinjection heuristic strategy. The a 0 j value is updated every 30 SLP iterations inside Eq. 6 as follows:

a 0 j = max a 0 j , 1 10 max (a a a 0 0 0 ) . (8)
This step reintroduces periodically many suppressed members and helps avoid local minima as shown later in the document.

Once convergence of Problem P 2 is reached, the solution is used to simplify the initial ground structure, removing the elements in which a < max(a a a)×tol. The physical cross-sectional member areas ã * are retrieved using Eq. 2. The tol parameter is set to 10 -3 in our analysis.

The next step is to set up the FEA of the updated ground structure to evaluate Ũ * . The linear system to solve for a truss is then:

f f f = B B B T q q q = B B B T D D Dε ε ε = B B B T D D DB B BU U U = K K KU U U, (9)
where D D D = diag (E ã * i /ℓ i ) and the stiffness matrix K K K = B B B T D D DB B B. Later, ã * , q * and Ũ * are used as the starting point of the full NLP formulation where the kinematic compatibility and the exact local buckling formulation are restored. The solver finally outputs the optimized structure variables a * , q * , U * .

In this work, we used the python library CVXPY [START_REF] Diamond | CVXPY: A Python-Embedded Modeling Language for Convex Optimization[END_REF] with the MOSEK solver to solve the relaxed LP problem P 2 . The full NLP problem P 1 is instead solved using IPOPT [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for largescale nonlinear programming[END_REF] , a large-scale nonlinear optimization package. The Jacobian and the Hessian matrices of the problem P 1 are calculated at every optimization step to allow faster convergence. As every state variable of the optimization is independent of the others, these responses are easy to calculate and will not be detailed there.

III. Application to a cellular wing box

The proposed method is evaluated on a simplified 3D wing box. We simulate a drone wing box with an Aspect Ratio (AR) of 11 and dimensions of L L L = [5500 mm, 500 mm, 125 mm]. The applied aerodynamic loads are the simultaneous application of lift and distributed torque that simulate the pitching moment. The lift is assumed as an elliptical distribution on the span-wise direction x and uniform distribution on the chord-wise direction y as such :

q(x) = q max 1 - x L x /2 2 , ( 10
)
where q max = 4n z Mg πL x = 1130 N/m, if we assume n z = 1, and M = 1000 kg. The pitching torque is calculated using the panel method assuming inviscid and incompressible flow. The torque is written as a linear distribution on the span-wise direction x:

t(x) = t max 1 - x L x , (11)
where t max = 0.208 [N/m]. We notice that for a wing with a such AR, the pitching moment is almost negligible compared to the lift. We also assume that the aerodynamic loads are not influenced by the deflection of the structure i.e., we neglect the aeroelastic interaction between fluid and structure. All the loads are applied statically. The load case is graphically represented, in Fig. 4a. The material used for the optimization is an aluminum alloy whose properties are schematized in Table 1.

The wing box is discretized in three different ways to investigate how the solution varies depending on the input parameters. We first conduct an optimization subdividing the wing box into 11×4×1 equal cells with dimensions of 500×125×125 mm each. Subsequently, the structure is discretized by 44×4×1 cubic cells with an edge length of 125 mm to investigate how the shape and the size of the repeating cell influence the optimized solution. Ultimately, the wing box is discretized by 44×4×1 octet-truss cells with fixed topology, comparing the advantages and disadvantages of the two approaches.

Every cell being optimized by the TTO algorithm is discretized by m = 3x3x3 equally spaced nodes that support a fully connected ground structure of 351 members (see Fig. 4b for a cubic representation of the cell ground structure). For simplicity, no influence of skin is taken into account in this study. The loads are applied on the center node of the upper face of the cells.

A. 11x4x1 cells TTO wing box

The wing box presented in Fig. 4a is subdivided into 11×4×1 equal cells with dimensions of 500×125×125 mm each. Doing that way, we obtain a span-wise distribution of cells that is comparable with the span-wise ribs distance of a classic wing design.

The optimized structure can be seen in Fig. 5, while Fig. 6 shows the topology of one cell. The suggested design is made by long spar-like elements well visible in Fig. 5c orthogonally cut by ribs formed by diagonally connected squares (see Fig. 5b and Fig. 6b). The lower part of the structure works, as expected, in tension, while the upper part in compression. As the optimization is buckling constrained, one can note that the compression members are more massive than tension ones (see Fig. 6a) with a maximum cross-sectional area of 325 mm 2 . The final mass of the structure is 24.940 kg. Figure 7 shows the convergence history of the analysis, pointing out the difference between the LP and the NLP iterations. It is interesting to see that the reinjection strategy for local buckling used in the LP optimization permit to "jump" across a local minimum to land on a more interesting solution. We see that all the compressed bars of the web of the first topology are substituted by tension bars and a square with a tension diagonal in the second topology (see Fig. 6). Prioritizing the bars in tension, the second topology achieves a lower mass. In that example, the NLP only slightly modify the mass (+0.2%) of the LP solution because the cell topology is not statically indeterminate.

B. 44x4x1 cells TTO wing box

To investigate the effect of the cell size on the solution, we run another optimization subdividing the structure into 44×4×1 cubic cells with an edge length of 125 mm.

The optimized topology shows long 2D spars made by the repeating cell drawn in Fig. 8a. The topology is Association for Structural and Multidisciplinary Optimization in the UK (ASMO-UK) similar to the 11×4×1 case, with a big compression member in the upper part of the cell, connected by thin tension bars to the lower part of the cell. The square sub-structure on yz plane of the previous solution (shown in Fig. 6b) is here rearranged into a single 2D column. The 2D topology of the proposed solution should not be able to withstand the torque loads. However, the force-equilibrium constraint of the optimization is fully satisfied. The 2D structure is generated by a numeric artifact. To suppress the numeric noise of the solution at the end of the optimization, the members are thresholded with a value of r = max (a a a) * 10 -5 . It turns out that, due to the low value of the applied torque with respect to the lift, the bars responsible for equilibrium in the y direction are suppressed. Figure 9 presents the convergence history of the optimizing algorithm. We can see here again that the reinjection strategy used permits to avoid a local minimum.

C. 44x4x1 cells octet-truss wing box

The TTO algorithm is here benchmarked against one of the most popular cell topologies present in the literature: the octet-truss (see Fig. 10b). The octet-truss is a cell known for having very good effective mechanical properties that attain about half the theoretical values of the upper Hashin-Shtrikman bounds [START_REF] Deshpande | Effective properties of the octet-truss lattice material[END_REF] .

To perform the benchmark, the wing box of Figure 4a is meshed using 44×4×1 octet-truss cells of 125 mm of edge size (see Fig. 10). The cross-sectional areas of the cell members are all equal. The section value is evaluated by performing a parametric optimization on this section value, constrained by stress and local buckling for every member of the structure. The optimization is performed using Altair OptiStruct.

D. Discussion

The numerical results of the three tested cases are presented in Table 2. Some interesting conclusions can be drawn if we first compare the 11×4×1 TTO with the 44×4×1 TTO case, and later the 44×4×1 TTO with the octet-truss case. We start observing again the topologies of the two TTO-optimized structures in Fig. 6d and Fig. 8a. The main difference between the two structures is that the section of the upper and the lower members, charged in compression and tension, is now comparable. This is due to the smaller length of the members of the 44×4×1 case and the subsequent change of failure mode of the compression member from local buckling to stress. Because of this more efficient use of the material, the final weight of the 44×4×1 structure is less than half of the 11×4×1 structure.

Cell type

When using the variable linking scheme, the optimization problem becomes more and more constrained the more cells we put in the domain. One would expect to have better performances when using fewer cells. But in our case, this is not true. This is due to the fact that, even if the cell discretization remained the same for the two optimizations (3×3×3 nodes with a fully connected ground structure), the final ground structures are not the same. This is especially important as we saw that the 44×4×1 case ground structure is made up of smaller bars that help maximize the structure efficiency with respect to the buckling constraint.

One shortcoming of the two TTO-optimized cell topologies is that they are mechanisms on their own, acting as a structure only once connected with the other cells.

Comparing now the weight of the TTO and the octet-truss topologies is lighter by almost one order of magnitude. This big gap could be explained by two distinct observations: 1. We notice that in the octet-truss there are no members orientated exactly along the x axis, while in the TTO optimized cell they are the most massive. This clearly tells us that this is the most efficient direction to put the material to get a strong cell. On top of that, the upper and lower faces of the cell present a cross design (see Fig. 10b) that works well for torque but not for tension and compression loading. A new study exploring what happens if we rotate the cell could be interesting.

2. The octet-truss is a cell that exhibits good homogenized elastic properties in all directions thanks to its numerous plane of symmetry. It is, thus, less suitable for structural applications where all the cells present similar loading conditions, as the cell will be as stiff and strong in every direction and not aligned to the principal stress directions.

Nevertheless, the octet-truss is a stretching-dominated cell that is statically determined, and its manufacturability causes fewer issues than TTO-optimized cells.

It is important to note that the comparison that we are presenting here does not take into account the weight of fasteners and joints necessary to lint the cells together. In the case of the solution with 44 cells span-wise that weight would be multiplied by at least a factor of 4.

Finally, the deflection of the three wing boxes is represented in Fig. 11. Even being more massive, the wing tip deflection of the octet-truss case is almost double of the TTO. This is again due to the inefficient cross topology shown by the octet-truss cell in the upper and lower faces.

IV. Concluding remarks

This paper presents a comprehensive algorithm to optimize ultralight cellular structures. The classic LP formulation is extended to take into account an important mode of failure of trusses: local buckling. The repetitive aspect of the structure is implemented using the variable linking scheme as we suppose no scale separation between the structure and the repetitive pattern.

The optimization problem is solved using a two-step algorithm: we first solve a relaxed formulation by the means of SLP that is used to approximate the highly non-linear solution. Later, a full NLP formulation is solved where all the original constraints are restored. A reinjection numerical technique has been proposed to mitigate the appearance of local minima as well.

The proposed algorithm is applied to a real-size aerospace application, a cellular 3D wing box with a span of over 5 m. Compared to a classic cellular topology, the octet-truss, we saw a 10-times weight improvement, clearly showing the interest in cell optimization.

However, in the present form, TTO-optimized cells could present non-trivial shapes and questionable kinematic (mechanisms) that worsen the manufacturing phase, suggesting an interesting area of further research.

Figure 2 :

 2 Figure2: The proposed two-step optimization algorithm. In the first part, the solution of a relaxed problem is used to simplify the fully connected cell ground structure. In the second part, the complete problem with the restored constraints is solved on the simplified ground structure.

Figure 3 :

 3 Figure 3: Relaxed NLP resolution flowchart using an SLP formulation.

Figure 4 :

 4 Figure 4: (a) represents the 3D wing box discretized by 44×4×1 cells with lift (red forces) and torque (blue forces) applied. (b) shows the 351 members of a 3×3×3 fully connected ground structure. The cubic cell shown here can be stretched to accommodate non-cubic shapes.

Figure 5 :

 5 Figure 5: Orthogonal and perspective views of the 11×4×1 cells optimized topology. (a) shows the xz view (b) shows the yz view (c) shows the xy view and (d) shows the perspective view. Blue bars are loaded in compression, red ones are loaded in tension, while gray ones are unloaded. The cells create a complex structure that mimics the design of a classic wing with 4 spars-and 11 ribs-like substructures. (b) is enlarged 4 times to better display the topology to the reader.

Figure 6 :

 6 Figure 6: Orthogonal and perspective views of the optimized cell for the 11×4×1 wing box subdivision. (a) shows the xy view (b) shows the xy view (c) shows the xy view and (d) shows the perspective view.

Figure 7 :

 7 Figure 7: The volume of the structure is plotted against the optimizer iterations for the 11×4×1 TTO case. The plot shows the LP run in red, followed by the NLP in blue. The graph points out the efficiency of the reinjection strategy used in the LP at iteration 30 to avoid local minima.

Figure 8 :

 8 Figure 8: Orthogonal and perspective views of the 44×4×1 cells optimized topology. (a) shows the perspective view of the optimized cell. (b) shows the perspective view of the complete wing box. (c) shows the xy view of the complete wing box. Blue bars are loaded in compression, red ones are loaded in tension, while gray ones are unloaded.

Figure 9 :

 9 Figure9: The volume of the structure is plotted against the optimizer iterations for the 44×4×1 TTO case. The plot shows the LP run in red, followed by the NLP in blue. The graph points out the efficiency of the reinjection strategy used in the LP at iteration 30 to avoid local minima.

Figure 10 :

 10 Figure 10: (a) shows the OptiStruct mesh for the 44×4×1 octet-truss cells. (b) displays the unit octet-truss cell with blue struts representing the cuboct lattice bounded by 8 edge tetrahedrons (in red). The optimized section value found is 85.027 mm, that is almost the same values of the maximum of the 44×4×1 TTO solution. The final optimized weight of the octet-truss case is 103.893 kg.

Figure 11 :

 11 Figure 11: Comparison of the span-wise wing deflection of the three candidate structures. The TTO optimized solutions are more stiff and lighter than the octet-truss thanks to a more efficient topology of the cell.

Table 1 :

 1 Elastic material properties used for the optimization.

Table 2 :

 2 Performance of the optimized designs with a different tassellation strategy.

		Repetitions	Cell size	Mass	Compliance Tip deflection Maximum area
			[mm]	[kg]	[J]	[mm]	[mm2]
	TTO	11×4×1	500×125×125 24.940	585.248	356	325.327
	TTO	44×4×1	125×125×125 10.867	782.138	467	87.093
	Octet-truss	44×4×1	125×125×125 103.893	876.269	977	85.027

Association for Structural and Multidisciplinary Optimization in the UK (ASMO-UK)


A-S Rollier

B Legrand

D Collard

L Buchaillot

The stability and pull-in voltage of electrostatic parallel-plate actuators in liquid solutions

This paper deals with parallel-plate electrostatic actuators in liquids. We study the stability conditions of such actuators and show that the pull-in effect can be shifted beyond one-third of the gap, and can even be suppressed. We demonstrate that the insulating layers of the actuator plates, which are originally designed to avoid any current leakages or short-circuits, play a major role in this phenomenon. Experiments are performed on fabricated devices; silicon nitride layers are used to completely encapsulate the actuator plates. The voltages required to close the actuator gap are measured in various liquids and compared to the values obtained by analytical calculations. This study gives guidelines for the design of parallel-plate actuators featuring in liquids either a binary-state operation when the pull-in effect occurs, or a continuous displacement within the full gap.

Introduction

Air or vacuum operation of electrostatic actuators is well controlled as well as their technological fabrication. Some systems have been commercialized, such as the digital micromirror device (DMD) invented in 1987 by Texas Instruments, which is a fast reflective digital light electrostatic switch [START_REF] Gove | DMD Display systems: the impact of an all-digital display[END_REF].

However, one of the most ambitious applications would be to make them usable in a liquid environment and thus to consider original biomedical applications such as cellular handling and characterization, DNA manipulation or device motion in liquid environment. Only a few devices have been actuated in liquid solutions, such as electrostatic comb-drives insulated using native silicon dioxide (SiO 2 ) [START_REF] Sounart | Electrostatic actuation without electrolysis in microfluidic MEMS Transducers[END_REF].

A great advantage of electrostatic actuation in a liquid environment comes from the relative dielectric constant of liquid, e.g., 80.1 for deionized water and 21 for isopropyl alcohol [START_REF] Janz | Nonaqueous Electrolytes Handbook[END_REF], which is expected to magnify the produced forces. This is quite interesting for low voltage CMOS compatible applications.

In this paper, we focus on the parallel-plate electrostatic actuators that exhibit the well-known pull-in effect in air, which limits the stable positioning of the actuator to one-third of the gap prior to collapse. Some authors have extended the stable operating range of such actuators by using a current pulse drive method [START_REF] Castaner | Analysis of the extended operation range of electrostatic actuators by current-pulse drive Sensors[END_REF] or adding external series capacitors [START_REF] Chan | Electrostatic micromechanical actuator with extended range of travel[END_REF]. In this work, we demonstrate that the travel range of the actuator in liquids is also affected by the plate insulating layers, even leading to the suppression of the pull-in effect.

Test devices are fabricated with standard silicon microtechnologies. Only a few materials are suitable to insulate the actuator structural parts: silicon nitride, silicon dioxide and various types of resists. The resists cannot be selected because the fabrication of electrostatic actuators requires high temperatures annealing steps. Furthermore, the releasing of the devices is performed by a wet etching of a sacrificial oxide layer. Thus insulating the structural parts with oxide is impossible, except if the oxide growth occurs after the releasing, as in [START_REF] Sounart | Electrostatic actuation without electrolysis in microfluidic MEMS Transducers[END_REF]. However, silicon dioxide presents a strong compressive stress so that it is difficult to produce flat devices. Hence, we chose a fabrication process based on a low stress silicon nitride (Si x N y ) encapsulation of electrostatic actuators.

Since the encapsulation is performed with low stress layers, the structure exhibits no deflection and the topography is very flat. Finally, the breakdown field of Si x N y is higher than that of SiO 2 (respectively, 7 MV cm -1 against 1.7 MV cm -1 , measured in our laboratory) which enables the use of thinner encapsulation thicknesses.

In order to have actuation in liquids and avoid charged bilayer formation and potential screening, an ac drive signal is preferred to the dc polarization of the actuator [START_REF] Sounart | Electrostatic actuation without electrolysis in microfluidic MEMS Transducers[END_REF].

The voltages required to close the actuator gap are measured in various liquids and compared to the ones obtained through analytical calculations. This study gives guidelines to design parallel-plate actuators featuring either a binarystate operation when the pull-in effect occurs, or a continuous displacement within the full gap in liquids.

Analytical modelling of actuation

The parallel-plate actuator can be typically considered as a capacitor having a fixed electrode while the second electrode is connected to a mechanical spring and is free to move vertically once an electrostatic force is applied; see figure 1. The electrostatic force acts to reduce the gap between the two electrodes. The restoring mechanical force that is exerted by the spring on the mobile electrode is proportional to the spring stiffness and to its deformation [START_REF] Senturia | Microsystem Design[END_REF].

By using the principle of virtual work, the electrostatic force acting on the upper electrode is given by

F elect = 1 2 × ε 0 ε × S × V 2 t1ε ε1 + t2ε ε2 + (g -x) 2 x ∈ [0,g[( 1 )
where k is the spring constant, g the gap between the electrodes, V the supply voltage, x the displacement of the mobile part, S the surface of the mobile part, t 1 the thickness of insulating layer 1, ε 1 the relative permittivity of layer 1, t 2 the thickness of insulating layer 2, ε 2 the relative permittivity of layer 2 and ε the permittivity of the medium between the two electrodes.

The equilibrium position of the actuator is reached once the electrostatic and mechanical forces compensate. By using (1),

F elect + F mecha = 0 ⇔ 1 2 × ε 0 ε × S × V 2 t1ε ε1 + t2ε ε2 + (g -x) 2 = k × x. ( 2
)
Solving equation ( 2) makes its possible to determine the position x of the actuator according to the voltage V, under the condition that a stable equilibrium solution exists. As we do not introduce any dynamic and inertial effects in [START_REF] Sounart | Electrostatic actuation without electrolysis in microfluidic MEMS Transducers[END_REF], this modelling is applicable, from the mechanical point of view, to quasi-static cases, where mobile plate velocity and acceleration induce no resultant forces. This limits the range of working frequencies to those much lower than the natural resonant frequency of the actuator.

From the electrical point of view, the previous calculation of the electrostatic force in equation ( 1) is valid for a nonconducting liquid (i.e. without ions). For a conducting fluid, the voltage excitation has to be alternative in order to avoid electrode potential screening, as presented in [START_REF] Sounart | Electrostatic actuation without electrolysis in microfluidic MEMS Transducers[END_REF].

Assuming an ac driving voltage

V ac = V rms √ 2cos(ωt). (3)
The electrostatic force becomes

F elect-ac = 1 2 × ε 0 ε × S × V 2 rms t1ε ε1 + t2ε ε2 + (g -x) 2 × [1 + cos(2ωt)]. (4)
The induced electrostatic force is then a combination of a static force and of a force at the pulsation 2ω.

In this case, a quasi-static displacement of the actuator is achieved if the actuator mechanical response at pulsation 2ω is negligible compared to the actuator resonance frequency, f resonance ,thatistosay , 2ω ≫ 2π × f resonance .

(

) 5
Under this assumption [START_REF] Chan | Electrostatic micromechanical actuator with extended range of travel[END_REF], the modelling we propose remains valid and can be used in the case of ac driving voltage by simply replacing in equation (2) V by the rms value of the ac voltage V rms .

Stability study

To determine the range of the stable actuation, the force equilibrium condition in (2) can be rewritten as

h(x) = f(V). (6)
With for x ∈ [0,g[:

h(x) = t 1 × ε ε 1 + t 2 × ε ε 2 + (g -x) 2 × x 0( 7 ) f(V) = 1 2 × ε 0 ε × S × V 2 k 0. ( 8
)
Solving equation ( 6) yields the equilibrium position of the actuator mobile part provided that a solution exists. The stable equilibrium actuation zone is defined by the whole actuator position x in which a stable equilibrium occurs between the electrostatic force and the mechanical force.

Graphically, in figure 2, the function h(x) versus the displacement exhibits a maximum x 1 and a minimum 6). The stable equilibrium position extends from x = 0tox = min(g, x 1 ) according to [START_REF] Ho | Micro-electro-mechanical-systems (MEMS) and fluid flows[END_REF]. Arrows show the displacement of the actuator when a voltage is applied.

x 2 Displacement x h(x) x 1 ) (V f
defined by

x 1 = 1 3 t 1 ε ε 1 + t 2 ε ε 2 + g (9) x 2 = t 1 ε ε 1 + t 2 ε ε 2 + g. ( 10
)
For a fixed voltage V, two cases must be considered for the resolution of equation ( 6

). If x ∈ [0,g[, f(V) > h max (x),
there is no solution because F elect >F mecha and the mobile part collapses onto the buried electrode. There is a new equilibrium in x = g as the reactive force (equivalent to infinite spring stiffness in figure 1) has to be taken into account.

If x ∈ [0; g[, f(V) < h max (x), as the displacement is limited to x = g,andasg<x 2 (see equation ( 10)), there are two equilibrium solutions for equation ( 6), a stable one before x = x 1 and an unstable one after x = x 1 .

Thus, solving equation ( 6) yields the extension of the stable actuation zone corresponding to one-third of the initial air-gap. In dielectric liquids, the thickness and permittivity of insulating layers must be taken into account to determine the stable actuation zone.

Then this stable actuation zone ranges from x = 0t o x = min(g, x 1 ). Since the parallel-plate is voltage controlled, by using ( 7), ( 8), ( 9) and ( 10), the voltage stable actuation zone can be rewritten as

x ∈ [0, min(g, x 1 )] ⇔ V ∈ [0,V(min(g, x 1 ))] (11) V(g)= lim x→g V(x)= 2 × k × g × t1ε ε1 + t2ε ε2 2 ε 0 ε × S (12) V(x 1 ) = V pull-in = 8 × k × t1ε ε1 + t2ε ε2 + g 3 27 × ε 0 ε × S . (13)
Equation ( 12) corresponds to the voltage required to close the gap g if we are in a stable equilibrium configuration (g<x 1 ). In contrast, if g>x 1 the equilibrium is unstable and the pull-in effect occurs (see equation ( 13)).

The stable actuation zone in ( 8) can be modified to allow controlled motion in the entire gap and to overcome the pull-in effect after initial actuation gap so that pull-in effect does not occurs and the actuator motion is well controlled. [START_REF] Degani | Pull-in study of an electrostatic torsion microactuator[END_REF]. The relative permittivity is plotted versus the normalized thickness. t eq and ε eq are, respectively, the thickness and the permittivity of the insulating layers. With equation ( 9), we can obtain stable actuation within the whole gap, if

V pull-in V pull-out

x 1 >g ⇔ t 1 ε 1 + t 2 ε 2 > 2 ε × g. ( 14
)
If the condition ( 14) is not observed, then the equilibrium becomes unstable and the moving part collapses; it is the pull-in effect (figure 2). Indeed, if by construction, the previous condition ( 14)is respected, then the displacement of the moving part can be controlled within the entire gap.

Figure 3 represents the stability study with respect to the parameters of equation [START_REF] Degani | Pull-in study of an electrostatic torsion microactuator[END_REF]. To simplify the representation, the thickness and the permittivity of the two insulating layers are chosen to be equal; thus t 1 = t 2 = t eq and ε 1 = ε 2 = ε eq .

Motion-voltage characteristic

Solving equation (6) yields the x position of the actuator according to the applied voltage V.

Voltage-displacement characteristic of the actuator is schematically shown in figure 4(a) when the condition given in equation ( 14) is not respected. For a voltage lower than the pull-in voltage, V pull-in , the actuator motion can be controlled until the x 2 position. At V = V pull-in , the upper electrode collapses on the lower one. To release the contact, the voltage has to be reduced until the V pull-out value, which is the second threshold voltage of this hysteresis cycle. Then the actuator recovers a position in the stable actuation zone.

The pull-out voltage depends on the quality of the contact between the mobile part and the buried electrode. In fact, an intimate gap is achieved between the two layers depending on surface roughness, and interactive forces such as Van der Waals. All of these added parameters have to be taken into account to determine the expression of V pull-out [START_REF] Rebeiz | RF MEMS: Theory, Design, and Technology[END_REF].

Figure 4(b) represents the case where the condition ( 14)is fulfilled, and the pull-in effect does not occur. Consequently, the motion of the mobile part can be controlled in the entire gap and there is no hysteresis cycle.

For example, a parallel-plate actuator with g = 2 µmand 300 nm thick of Si x N y insulating layers (ε eq = 8) does not exhibit the pull-in effect in water (ε = 80.1). In fact, V pull-in occurs for a 2.6 µm displacement (see equation ( 13)) which exceeds the physical gap g (see point A in figure 3).

Whereas, according to point B in figure 3,2 4µmo f low stress silicon nitride would have been required in order to actuate the same parallel-plate actuator in air under stable equilibrium condition, which is unrealistic.

Consequently, we chose to encapsulate the structural part of the actuator with 300 nm Si x N y in order to control the motion in the entire gap during actuation in water.

Fabrication process

Three types of parallel-plate actuators have been fabricated using the same technological process: a simple parallel-plate device, a torsional actuator with a square electrostatic plate for actuation, two torsional beams and a lever to increase the deflection, and finally a membrane clamped in its four corners; see figure 5.

Starting from a silicon wafer, a 350 nm dry oxidation is performed and a 100 nm Si x N y (the exact stoichiometry has been optimized to reduce residual tensile stress) layer is deposited at 800 • C to insulate the substrate; then, a 350 nm n-doped polycrystalline silicon layer is deposited by low pressure chemical vapour deposition (LPCVD), transferred by photolithography and SF 6 reactive ion etching (RIE) to pattern buried electrodes. To insulate the electrodes during actuation, a 300 nm Si x N y layer is used; see figure 6(a). A2µm sacrificial layer of low temperature oxide (LTO) is deposited by LPCVD and anisotropically etched in a CF 4 and CHF 3 plasma; see figure 6(b).

Since actuators are operated in a liquid environment, to prevent any short-circuits and electrical leakages through the liquid medium, an encapsulation with Si x N y is achieved.

As the breakdown field of Si x N y deposited in our laboratory is 7 MV cm -1 , the 300 nm thick Si x N y layer is able to withstand a breakdown voltage of 200 V.

The first level of encapsulation is made of a 300 nm thick Si x N y layer anisotropically etched in order to open the electric contacts; figure 6(c). The structural material is a 2 µm n-doped polycrystalline silicon, which is anisotropically etched by SF 6 RIE; see figure 6(d). The second 300 nm silicon nitride for the front side insulation is deposited and patterned resulting in the complete encapsulation of the device; see figure 6(e).

In order to relax the mechanical stress in the stacked layers, a thermal annealing step is performed at 1100 • C during 3 h. The last step consists in the release of the structure in a hydrofluoric (HF) aqueous solution; see figure 6

(f ).

A parallel-plate actuator is shown in the scanning electron microscopy (SEM) picture in figure 7 of Si x N y . In order to validate the fabrication process, an interferometrical profile image has been performed along the cantilever (figure 7(c)). We can see that the profile topography is flat; thus no deflection appears which means that there is no notable stress gradient in the structure.

In order to precisely see the Si x N y encapsulation, a crosssection is made on a very small structure, a 2 µm width and 2 µm thick beam (figure 8

Experimental characterization

Experiments have been performed in different solutions with various permittivity and conductivity such as air, pure deionized water, isopropyl alcohol (IPA) and tap water. The three types of actuators have been tested, parallel-plate, torsional actuator and membrane. The devices have been tested using the experimental set-up shown in figure 9. A liquid reservoir is used to introduce the fluid into the sample, whilst keeping electrical contacts in a dry environment. The contacts are made by using two micromanipulators. A probe station enables the displacement visualization by way of an optical microscope. The dc and the ac polarization are generated by a waveform generator, connected to the device contacts and in parallel with an oscilloscope to monitor the voltage.

Electrode potential screening

As actuation occurs in dielectric liquid environment, during actuation ions accumulate in a electrical double layer (EDL) at the electrode interface; thus there is a phenomenon of potential screening [START_REF] Sounart | Frequency-dependent electrostatic actuation in microfluidic MEMS[END_REF]. Then, as deionized water presents an ionic strength around 10 -7 mol l -1 , potential screening is also possible and has been estimated from the Gouy-Chapman model [START_REF] Hunter | Zeta Potential in Colloid Science: Principles and Applications[END_REF]b y [START_REF] Sounart | Frequency-dependent electrostatic actuation in microfluidic MEMS[END_REF]. To avoid EDL, actuation must be driven by ac signals. The frequency of the signal to be employed depends on the dimensions of the structure (the dielectric layer thickness t eq and the gap g), and the liquid conductivity. Sounart et al [START_REF] Sounart | Frequency-dependent electrostatic actuation in microfluidic MEMS[END_REF] give us a reference chart for determining the actuation frequency. For deionized water, the minimum applied frequency must be 13 kHz. In the case of IPA, the a Value from [START_REF] Sounart | Frequency-dependent electrostatic actuation in microfluidic MEMS[END_REF] with a 300 nm thickness of dielectric layer and a 2 µm gap between the mobile part and the buried electrode. b Value from [START_REF] Janz | Nonaqueous Electrolytes Handbook[END_REF]. a Values from Si x N y [START_REF]Physics of Semiconductor Devices[END_REF].

driving frequency has to be higher than 17 kHz; in the case of tap water, the critical frequency is 1 MHz (table 1). Thus, the devices have been tested at 1 MHz, in the different fluids, in order to avoid charged bilayers and potential screening, which is likely to decrease the electrostatic force and thus lose the actuation [START_REF] Pomes | Fluorine migration in a soil bed submitted to an electric field: influence of electric potential on fluorine removal Colloids Surf[END_REF][START_REF] Ho | Micro-electro-mechanical-systems (MEMS) and fluid flows[END_REF].

Parallel-plate actuator

A parallel-plate actuator with the dimensions given in table 2 has been fabricated to test the out-of-plane displacement; the optical image in figure 10 was taken during actuation in tap water. Figure 10(a) corresponds to the upper position at V = 0 V. In figure 10(b) the 2 µm gap between the mobile part and the buried electrode is closed.

The structure has been firstly tested in air with and without the Si x N y encapsulation. Since condition of equation ( 14)i s not respected, the pull-in effect occurs, respectively, at 7 V and 8 V either without or with the nitride encapsulation (table 3). We can see that the two 300 nm thick dielectric layers induce an increasing in the actuation voltages by 15%.

If we raise the relative permittivity of the medium by actuating the parallel-plate in IPA, the condition ( 14) is not a Values from Si x N y [START_REF] Janz | Nonaqueous Electrolytes Handbook[END_REF].

respected any longer and the pull-in effect appears at 4 V (table 3). The actuation voltage is reduced by 50% compared to the one in air.

The last two tests occur in deionized and in tap water; thus the condition ( 14) is satisfied and the motion of the mobile part is controlled within the entire gap, the voltage to close the 2 µm gap is around 6 V (table 3).

These preliminary tests show that the electrostatic actuation principle can be applied to a large range of fluids with various permittivity and conductivity.

Other structures

To study electrostatic actuation in liquid medium, the other types of actuators have been tested with the same experimental set-up than in figure 9. An interesting type of actuation in liquid environment is the torsional electrostatic actuation [START_REF] Toshiyoshi | Electrostatic micro torsion mirrors for an optical switch matrix[END_REF] which could be used to make a micro-gripper. To highlight this principle, parallelplate structures are maintained by fixed beams; thus during actuation, these beams are subjected to mechanical moments of torsion [START_REF] Degani | Pull-in study of an electrostatic torsion microactuator[END_REF][START_REF] Satter | Macromodeling of an electrostatic torsional actuator Transducers[END_REF]. The picture in figure 11 shows a torsional actuator with a square electrostatic plate for actuation and two torsional beams; the deflection is observed on a lever beam, dimensioned in table 4.

This device has been tested in air, in tap water and in IPA (table 6). The V pull-in voltages are lower than the ones in the previous design, because of the higher surface of the actuation (see equation ( 13)).

As actuation occurs in a liquid environment, actuation involves fluid displacement. To test the robustness of the

V pull-in = 0.25 V V pull-in = 0.12 V V (2 µm) = 0.19 V Experimental voltage c V pull-in = 0.3 V V pull-in = 0.16 V V (2 µm) = 0.23 V Membrane d
Simulated voltage f V pull-in = 44.68 V V pull-in = 16.37 V V (2 µm) = 33.4 V Experimental voltage c V pull-in = 45 V V pull-in = 15.9 V V (2 µm) = 33.9 V structures, a membrane clamped in its four corners was fabricated with a 300 nm Si x N y encapsulation. The dimensions of the membrane are given in table 5. Figure 12 presents a microscope image of the membrane. The membrane has been actuated in air, in tap water and in IPA; the results are presented in table 6. We note a good agreement between experimental and modelled values. This type of design presents a strong effective stiffness, which explains the high values of the pull-in voltages (see equation ( 13)).

In figures 11 and 12, some holes were bored in the structural part to allow a better release of the membrane and the torsional actuator. Indeed, this release is carried out in liquid medium and these holes support the flow of HF.

Conclusion

In this paper, we have presented a novel optimized fabrication technique of parallel-plate actuators with a complete Si x N y encapsulation, allowing electrostatic actuation in a liquid environment. Analytical calculations show that the pull-in effect can be shifted far beyond one-third of the gap and even suppressed. The thicknesses and dielectric constants of the insulating layers play a major role in this phenomenon. Experiments have been performed on fabricated devices and actuation has been successfully observed on cantilevers, torsional actuators and membranes. The voltages required to close the actuators gap have been measured in IPA, DI water and tap water. They are in good agreement with the theoretical values. Further investigations will concern the dynamic behaviour of the actuators in liquids, which is strongly affected by a much higher damping coefficient in liquids than in air [START_REF] Lin | Miniature heart cell force transducer system implemented in MEMS technology[END_REF][START_REF] Maali | Hydrodynamics of oscillating atomic force microscopy cantilevers in viscous fluids[END_REF][START_REF] Bouhacina | Oscillation of the cantilever in atomic force microscopy: Probing the sample response at the microsecond scale[END_REF].

Figure 1 .

 1 Figure 1. Schematic view of a parallel-plate electrostatic actuator.

Figure 2 .

 2 Figure 2. Graphic representation of equation (6). The stable equilibrium position extends from x = 0tox = min(g, x 1 ) according to[START_REF] Ho | Micro-electro-mechanical-systems (MEMS) and fluid flows[END_REF]. Arrows show the displacement of the actuator when a voltage is applied.

Figure 3 .

 3 Figure 3. Graphical representation of condition[START_REF] Degani | Pull-in study of an electrostatic torsion microactuator[END_REF]. The relative permittivity is plotted versus the normalized thickness. t eq and ε eq are, respectively, the thickness and the permittivity of the insulating layers.

Figure 4 .

 4 Figure 4. Displacement versus voltage characteristic for electrostatic actuation with pull-in effect (a) and without pull-in effect (b).

Figure 5 .

 5 Figure 5. Schematic view of (a) a parallel-plate actuator; (b)the membrane clamped in its 4 corners and (c) the torsional actuator.

Figure 6 .

 6 Figure 6. Process flow: (a) buried electrodes patterning; (b) deposition and etch of the sacrificial layer; (c) first level of nitride encapsulation; (d) deposition of the structural layer; (e) second level of nitride encapsulation; (f ) releasing in a HF solution.

Figure 7 .Figure 8 .

 78 Figure 7. (a) SEM image of a parallel-plate actuator. (b)Close view of the nitride step covering. (c) Profile topography corresponding to dashed line in figure 4(a).

Figure 9 .

 9 Figure 9. Experimental set-up. A probe station is used to keep the contact and see the motion of the actuator. A test box enables the introduction and the removal of liquid, the electrical contact is taken out of the box. silicon nitride is presented in figure 8(b); we clearly see the 2 µm polycrystalline silicon surrounded by a 300 nm Si x N y layer.

Table 2 .

 2 Parameters of the cantilever type actuator. Parameter Value Actuator length LL = 250 µm Actuator width ww = 30 µm Actuator thickness tt = 2 µm Actuator mobile part surface SS = 6000 µm 2 Actuator effective stiffness kk = 0.8 N m -1 Gap between the actuator and the g = 2 µm electrode g Insulating layer 1 thickness t 1 t 1 = 300 nm Insulating layer 1 relative ε 1 = 8 permittivity a ε 1 Insulating layer 2 thickness t 2 t 2 = 300 nm Insulating layer 2 relative ε 2 = 8 permittivity a ε 2 Dielectric permittivity ε 0 ε 0 = 8.85 × 10 -12 Fm -1

Figure 10 .Table 4 .

 104 Figure 10. (a) Optical image during actuation in water when the gap is open. (b) Optical image during actuation in water when the 2 µm gap is closed. White arrows show actuation location.Table 4. Parameters of the torsional actuator tested. Parameter Value Actuator length LL = 300 µm Actuator width ww = 300 µm Actuator thickness tt = 2 µm Actuator mobile part surface SS = L × w = 0.09 mm 2 Actuator effective stiffness kk = 3.26 × 10 -9 Nmrad -1 Torsional beam length L t L t = 490 µm Torsional beam width w t w t = 5 µm Lever beam length L l L l = 1100 µm Lever beam width w l w l = 90 µm Gap between the actuator and the g = 2 µm electrode g Insulating layer 1 thickness t 1 t 1 = 300 nm Insulating layer 1 relative ε 1 = 8 permittivity a ε 1 Insulating layer 2 thickness t 2 t 2 = 300 nm Insulating layer 2 relative ε 2 = 8 permittivity a ε 2 Dielectric permittivity ε 0 ε 0 = 8.85 × 10 -12 Fm -1

Figure 11 .Table 5 . 1 a

 1151 Figure 11. (a) Optical image of a torsional actuator during actuation in water when the gap is open. (b) Optical image of a torsional actuator during actuation when the 2 µm gap is closed. White arrows show actuation location. Table 5. Parameters of the membrane tested. Parameter Value Actuator length LL = 378 µm Actuator width ww = 243 µm Actuator thickness tt = 2 µm Actuator mobile part surface SS = L × w = 91854 µm 2 Gap between the actuator and g = 2 µm the electrode g Insulating layer 1 thickness t 1 t 1 = 300 nm Insulating layer 1 relative ε 1 = 8 permittivity a ε 1 Insulating layer 2 thickness t 2 t 2 = 300 nm Insulating layer 2 relative ε 2 = 8 permittivity a ε 2 Dielectric permittivity ε 0 ε 0 = 8.85 × 10 -12 Fm -1 a Values from Si x N y [3].

Figure 12 .

 12 Figure 12. (a) Optical image of a membrane during actuation in water when the gap is open. (b) Optical image of a membrane during actuation when the 2 µm gap is closed. White arrows show actuation location.

Table 1 .

 1 Minimum driving frequency in dielectric liquid.

	Liquid	Conductivity b Frequency a
	Isopropyl alcohol 0.09 µScm -1 17 kHz
	Deionized water	0.07 µScm -1 13 kHz
	Tap water	5 µScm -1	1MHz

Table 3 .

 3 Theoretical and experimental voltage for the parallel-plate tested in different liquid medium with the dimension in table 1.

		Relative	Theoretical voltage to	Experimental voltage
	Actuation	permittivity a to close the 2 µm gap (V) to close the 2 µm gap (V) Difference (%)
	Air (without Si x N y )1 Air (with Si x N y )1 Isopropyl alcohol (1 MHz) 21.3 Deionised water (1 MHz) 80.1 Tap water (1 MHz) 80.1	V pull-in = 7.31 V pull-in = 7.73 V pull-in = 3.78 V (2 µm) = 6.37 V (2 µm) = 6.37	V pull-in = 74 V pull-in = 83 . 4 V pull-in = 45 . 8 V (2 µm) = 65 . 8 V (2 µm) = 65 . 8

a Values from Si x N y

[START_REF] Janz | Nonaqueous Electrolytes Handbook[END_REF]

.

Table 6 .

 6 Simulated and experimental voltage for torsional actuator and membrane in different liquid mediums with the dimension in table 3.

	Type of structure	Actuation	Air (with Si x N y )	Isopropyl alcohol (1 MHz) Water (1 MHz)
		Relative permittivity e	1	21.3	80.1
	Torsional actuator a Theoretical voltage b

Acknowledgment

The authors would like to thank the valuable work of the IEMN technical staff, and the 'Délégation Générale pour l'Armement' for the financial support of A-S Rollier.