A two-scale model for one-dimensional arrays of cantilevers and its verification

We present a simpliﬁed model of the mechanical behavior of large arrays of cantilevers in the dynamic operating regime. The supporting base is assumed to be elastic thus cross-talk effect between the canti-levers is taken into account. Beforehand, the model has been mathematically justiﬁed starting from a thin plate model, using the two-scale approximation theory issued from homogenization theory and taking into account the strong heterogeneities of the system. The resulting model is not standard, so in this paper dedicated to its veriﬁcation, we focus on some of its features in particular those related to the structure of its eigenmodes by both a qualitative and a quantitative analysis. parallelepiped cantilevers the un-clamped end of cantilevers, two one for free ends and one for ends with as in Atomic Force


Introduction
Cantilever arrays are used in a variety of application [1][2][3][4] including arrays of Atomic Force Microscope (AFM) as the Millipede from IBM dedicated to data storage [5]. Their direct numerical simulation, based on classical methods like Finite Element Methods (FEM), is prohibitive for today's computers, at least in a time compatible with the time constraints of a designer. The group of Bamieh has published a model of arrays of cantilevers [6] which takes into account electrostatic coupling between cantilevers and which derivation is based on purely phenomenological arguments. Besides, one of the authors has published a model for arrays AFMs in the elastostatic regime [7] and preliminary results for modeling their elastodynamic regime [8]. Both models are built by a rigorous mathematical approach based on an asymptotic method related to one of those used in the homogenization theory, see Refs. [9][10][11]. The derivation combines the concept of two-scale convergence (also called the unfolding method) together with this of strongly heterogeneous media to take into account the large difference between stiffness of the frame and this of the cantilevers having relatively large local motions. This type of two-scale modeling used for arrays of moving parts is not standard so that before use, one must check its validity because its mathematical justification alone cannot guarantee its relevance.
The main expected advantage of this two-scale model is that it should require little computing resources compared to classical simulations. In fact, we observe that using this model yields a very significant speedup of the computation, which is even true for large matrices. This point will not be much discussed along the paper which is focused on the approximation problem and not on the computation performance. Here, the presentation focuses first on special qualitative properties of the model and of its solutions from the point of view of the modal structure. In particular, investigations have been conducted to understand how the eigenmodes of the microstructures (the cantilevers) and those of the macrostructure (the base) are similarly organized in the two-scale model and in a finite element model of the three-dimensional elasticity problem. In addition to quantitative observations we provide an analysis based on quantitative comparisons that allows to point out some strengths and weaknesses of the two-scale model. In the same way that the mathematical method of derivation is quite general, we believe that the conclusions we draw for cantilever arrays are in fact more general and could apply to other arrays of micro and nanosystems.
The rest of the paper starts with the two-scale model statement in Section 2, followed by Section 3 that explains some features of the modal structure of the two-scale model. Finally, the comparison with the finite element model are detailed in Section 4.

2.
The two-scale model of cantilever arrays

Geometry of the problem
We consider a one-dimensional array of cantilevers, see Fig. 1a. It is comprised of a rectangle parallelepiped base crossing the array in the x 1 direction which is clamped to its ends and in which rectangle parallelepiped cantilevers are fixed. Concerning the unclamped end of cantilevers, two cases are considered, one for free ends and one for ends equipped with tips as in Atomic Force Microscopes. The array is a mere periodic repetition of a same cell in the direction x 1 . We assume that the number N of cells is sufficiently large in a sense to be specified in the next Section. Note that along the paper the superscripts B and C are referring to base and to cantilevers.

Two-scale approximation
We introduce the parameter e ⁄ = a/N inversely proportional to N, proportional to a parameter adefined hereafter, and assumed to be sufficiently small. Each point of the three-dimensional space, We consider the distributed field u(x) of elastic deflections in the array and we introduce its two-scale transform defined bŷ u e ðx 1 ; yÞ¼uðx c þ yÞ; for any x = x c + y which is constant, with respect to its first variable x 1 over each cell. Since it depends on the ratio e ⁄ , then it may be approximated by an asymptotic field, denoted by u 0 , obtained for a large number N of cells or equivalently when e ⁄ approaches (mathematically) 0 where O(e ⁄ ) is a function tending to zero when e ⁄ vanishes. The approximation u 0 is called the two-scale approximation of u. In general the partial function x 1´u 0 (x 1 ,.) is continuous unlike x 1 #û e which is piecewise constant. Now, we consider that the field of elastic deflections u is a solution of the Love-Kirchhoff thin elastic plate equation in the whole mechanical structure, including the base and the cantilevers. We assume that the ratio of cantilever thickness h C to base thickness h B is sufficiently small, namely This assumption implies that the ratio e ⁄4 of cantilever stiffness to base stiffness is sufficiently small. It is under this special choice that the asymptotic model constitutes an approximated model for thin cantilevers with relatively large displacements (but still in the framework of linear elasticity) acting on the base in a nonnegligible : The asymptotic anal-ysis shows that u 0 does not depend on the cell variable y in the base and so depends only on the spatial variable x 1 . Next, we remark that u 0 (x 1 ,y) is a two-scale field, and therefore cannot be directly used as an approximation of the field u(x) in the actual array of cantilevers. So, an inverse two-scale transform is to be applied to u 0 . However, we remark that x 1´u 0 (x 1 , y) is continuous, and so u 0 does not belong to the range of the two-scale transform operator and it has no pre-image. Hence we introduce an approximated inverse of the two-scale transform, vðx 1 ; yÞ# vðxÞ; in the sense that for any sufficiently regular one-scale function u(x) and two-scale function v(x 1 , y), Several choices are possible for vðxÞ, we define it only for functions v which restriction v B to the base is independent of the variable y.
in each cantilever where h:i x 1 represents the average in x 1 over the width of a cantilever. In total, we retain u 0 as an approximation of u in the actual physical system. Note that for the model in dynamics, the deflection u(t, x) is a time-space function. In our analysis we do not introduce a two-scale transform in time, so the time variable t acts as a simple parameter.

Model description
Now, we describe the model satisfied by the two-scale approximation u 0 (t, x 1 , y)o fu(t, x). Remark that as the deflection u in the Kirchhoff-Love model is independent of x 3 , thus u 0 is independent of y 3 . For further simplicity, we neglect cantilever torsion effect i.e. the variations of y 1´u 0 (t, x 1 , y).We denote by L B and L C⁄ the length of the base and the length of the cantilevers scaled by e ⁄ so the approximated model is posed in the rectangle X ={(x 1 , y 2 ) 2 (0, L B ) Â (0, L C⁄ )}, filled by an infinite number of cantilevers, see Fig. 2.  Thus, the cantilever motion is governed by a classical Euler-Bernoulli beam equation, in the microscopic variable y 2 , where ' C⁄ is the scaled width of the cantilevers, m C is a linear mass, E C the cantilever elastic modulus, I C the second moment of cantilever section, and f C a load per unit length in the cantilevers. This model represents the motion of an infinite number of cantilevers parameterized by x 1 . The base deflection are solution to the Euler-Bernoulli equation in the macroscopic variable x 1 posed on the lower boundary y 2 =0ofX, and f B are a linear mass, the base elastic modulus, the second moment of section of the base, a cantilever-base coupling coefficient and the load per unit length in the base. The term Àd B @ 3 y 2 ...y 2 u 0 jjunction is a distributed load originating from shear forces exerted by cantilevers on the base at base-cantilever junctions. In the model, the cantilevers appear as clamped in the base. So at the base-cantilever junctions, because @ y 2 u 0 ¼ 0 in the base. The equations of the free ends are and those of the ends equipped with a rigid part (usually a tip in Atomic Force Microscopes) are at junctions between elastic parts and rigid parts. Here, J R is a matrix of moments of the rigid part about the junction-plane, f R 3 is a load in the y 3 direction, F R 3 is a first moment of loads about the junction-plane, and F R 2 the first moment of loads in the y 2 direction about the beam neutral plane. Finally, the base clamping conditions are The loads f C , f B and f R in the model are some asymptotic loads which, in general, cannot be built from the physical problem. In a practical computation, they are simply replaced by the two-scale transformŝ f C ;f B andf R .

Structure of eigenmodes
There is a countable infinite number of eigenvalues k A and eigenvectors u A (x 1 ,y 2 ) associated to the model which equation has been stated in [8]. For convenience, we parameterize them by two independent indices i and j, both varying in the infinite countable set N.

Qualitative properties of the modal structure of cantilever arrays
We consider a one-dimensional silicon array of N = 10 cantilevers, with base dimensions 500 lm Â 16.7 lm Â 10 lm, and cantilever dimensions 41.7 lm Â 12.5 lm Â 1.25 lms oa = 2.1, see Fig. 3 for the two possible geometries, with or without tips. Then, for arrays of 15 and 20 cantilevers, only the cantilever width is changed so that to keep the same characteristic values of k A ij ; as it is seen in Fig. 4.
We have carried out our numerical study on both cases, but we limit the following comparisons to cantilevers without tips, because configuration including tips yields comparable results.
We restrict our attention to a finite number n B of eigenvalues k B i in the base. Computing the eigenvalues k A , we observe that they are grouped in bunches of size n B accumulated around a clamped-free cantilever eigenvalues. A number of other eigenvalues are isolated far from the bunches. It is remarkable that the eigenelements in a same bunch share a same cantilever mode shape even if they correspond to different indices j as discussed in the special case of a 10-cantilever array. This is why, these modes will be called ''dominated by cantilever modes''. Isolated eigenelements share also a common cantilever shape, which looks like a first clamped-free cantilever mode shape excepted that the clamped side is shifted far from zero. The induced global mode u A is then dominated by base deformations and therefore will be We discuss the comparison with the modal structure of the three-dimensional linear elasticity system for the cantilever array discretized by a standard Finite Element Method. The eigenvalues of the three-dimensional elasticity equations constitute also an increasing positive sequence that accumulates at infinity. As for the two-scale model, its density distribution exhibits a number of concentration points and also some isolated values. Here bunch sizes equal the number N of cantilevers, see sub-figs. 1, 3 and 5 in Fig. 4 representing eigenmode distributions for N = 10, 15 and 20. Extrapolating this observation shows that when the number of cantilevers increases to infinity bunch size increases proportionally. Since the two-scale model is an approximation in the sense of an infinitely large number of cantilevers, this explains why the two-scale model spectrum exhibit mode concentration with infinite number of elements. This remark provides guidelines for operating mode selection in the two-scale model. In order to determine an approximation of the spectrum for an N-cantilevers array, we suggest to operate a truncation in the mode list so that to retain a simple infinity of eigenvalues k A ij i¼1;...;Nand j2N : To conclude with the qualitative observation, we have remarked that the missing eigenmodes in the two-scale model correspond to physical effects not taken into account in the two Euler-Bernoulli models for the base and the cantilevers.

Quantitative verifications
The presentation is focused on the case of a 10-cantilever silicon array, i.e. for N = 10. The results relate to the first 40 modes in the FEM model and to the eigenelements k A and j 2 {1,2,3}, the latter being listed in Table 1. Note that the computation time is 0.76s for the modes of the two-scale model implemented in a nonoptimized MATLAB Ò code versus 88.14s for the finite element modes using COMSOL Ò with 20,859 quadratic elements with a regular laptop. We stress the fact that the N-eigenvalue bunches are not corresponding to a single row in Table 1 i.e. not corresponding to a single j. This is because the modes dominated by the deformation of base are interposed between the clusters of modes dominated by the deformations of cantilevers. The counterpart in terms of base modes is that they follow each other on consecutive columns but with possible line breaks.
To conduct a quantitative comparison of eigenvalues, it is required to match the modes of the two-scale model with those of the finite element model. Because of the proximity of many eigenvalues, a tool like the conventional Modal Assurance Criterion (MAC) is necessary to discriminate them, see [12].
it is equal to one if the shapes are identical and to zero when they are orthogonal in the sense of the inner product hÁ,Ái. Each subspace of eigenvectors u ref corresponding to a quasi-multiple eigenvalue is rotated so as to optimize the MAC matrix. The results are in Fig. 5 where the modes u A are arranged in the order such that the index i varies faster than the index j. The inner product is based on a sum over 300 points distributed along six parallel lines in the base and over 6 Â 10 Â N points along six lines in each cantilever. In both cases the six lines are along the four edges and along the central axes of the upper and lower faces and the points are regularly spaced. The FEM computation has been carried out with 20,859 elements. All modes u ref from the finite element model which are not sufficiently correlated with a mode u A i.e. with a MAC lower than 0.5 are not considered for comparison because they correspond to physical effects not modeled by the Euler-Bernoulli models. Some modes u A seem to correlate well with several modes u ref , like the eigenmodes 2, 11 and 12 so an additional criterion for selection should be applied. The most general method would be to add more points in the inner product, but here it was enough to eliminate the    Fig. 7b. Note that errors are far from being uniform among eigenvalues. In fact, the main error source resides in a poor precision of the Euler-Bernoulli model for representing base deformations in few particular cases. A careful observation of Finite Element modes shows that base torsion can be predominant for some modes, such as in the first mode of the first cantilever mode bunch.
In Fig. 7a, the distinction between the base modes and the cantilever modes is also marked. Their distinction could be done from the ratio of the amplitudes of deformation in the base and in the cantilevers. An equivalent way is to use the sensitivities with respect to characteristic parameters of the two modes of deformations. To find the influential parameters, the sensitivities of the model through parameter variations is established using a first-order finite difference method applied to the eigenvalues. The results are presented in Fig. 8 where all parameters have been tested, i.e. the Young's modulus, the volume mass, the thicknesses, the lengths and the widths.
Their values are denoted by E, q, h B , L B , l B , h C , L C , and l C where the superscripts B and C stand for base and cantilevers. The eigenvalues are mainly sensitive to the thickness h B of the base, to the length L C of the cantilevers, and for a lesser extent to the thickness h C of the cantilevers. Most of the eigenvalues are sensitive to only one of the two parameters h B or L C then they can be identified as a base mode or as a cantilever mode. The cantilever modes are clearly organized in clusters of N = 10 modes separated by base modes. At their interfaces some modes are almost equally sensitive to base and to cantilever parameters, they are referred as mixed mode in Fig. 9. However, for simplicity they are considered as base modes in Fig. 7. To illustrate the distinction between the three kind of modes, the Fig. 9a and c present a base mode and a cantilever mode when the Fig. 9b and d show two mixed modes.

Conclusion
A two-scale model of cantilever arrays in dynamic regime has been presented. Its derivation, previously carried out, uses a theory of strongly heterogeneous homogenization in which the cantilevers play the role of soft parts. In the resulting model, only the transverse displacement was retained. We analyzed it and compared it to ordinary finite element simulations from the viewpoint of modal structure. A special emphasis was placed on the distinction between modes dominated by base deformation and those dominated by cantilever deformation. We observe that this concept can be met in various kind of arrays of coupled systems, so the analysis methodology could be re-used in other applications. Globally, the two-scale model and the direct finite element model provide comparable results but some modes are not absolutely correct. A possible way to improve the current model would be to take into account the three mechanical displacements rather than the transverse displacement only.