A simple method for extracting material parameters of multilayered MEMS structures using resonance frequency measurements

Recently, multilayered structures have been utilized in MEMS. Researched applications include uncooled microcantilever infrared focal plane arrays [1], radio frequency (RF) components [2], micromachined mirrors [3], etc. It is well known that MEMS devices are highly dependent on factors such as the Young's modulus and residual stress of the multilayered films. These properties determine both the final shape and the functionality of released micromachined structures, and should therefore be accurately assessed. The Young's modulus and residual stress for the single layer film have been widely studied by the cantilever deflections, wafer curvatures, displacements of variously designed microstructures, buckling lengths, membrane deflections, resonance frequency, pull-in voltages and fixed-fixed beam deflections. However, these methods are not easily extended to multilayered films. The Young's modulus of multilayered films has been measured without the residual stress [4–7], while the residual stress for multilayered films has been studied but still requires the Young's modulus of each layer to be known [8–12]. It is more desirable to directly measure both the Young's modulus and residual stress for the multilayered films simultaneously. An attempt has been made to characterize these two mechanical properties by the pull-in voltages [13, 14] or resonant frequency [15]. The resonant frequency of the multilayered beam fixed at two ends varies with the dimensions of the beam. By designing beam test structures for each layer with different widths, the Young's modulus and residual stress for an individual layer may be extracted by the measured resonant frequency of the multilayered beam. Taking into account the buckling or no buckling problem of the multilayered beam, this paper first introduces a model for the resonant frequency of the beam for each layer with different widths. An approach to extract the Young's modulus and residual stress for the individual layer is then developed. The validity of this approach has been studied using finite element modeling (FEM). As a typical example, test structures for a metal (gold)/polysilicon bilayer beam were fabricated using MEMSCAP PolyMUMPs process. A Polytech scanning laser Doppler vibrometer (LDV) system with hardware modules of MSV-400-M2 was used to measure the resonant frequency of the beam. The results show the proposed method produces results that agree with the values most often given in the literature.

This method uses multilayered micromachined doubly clamped beams that are excited to resonance and are measured using a scanning LDV. The resonance frequency predicted by a composite-beam model is used to extract the Young's moduli and stresses of all the individual thin-film layers in the beam. To easily extract these material parameters, test structures for the doubly–clamed beams with different widths are designed. Figure 1 shows the cross-section of the beam with different widths. The length of the beam is l. The width, thickness, density, Young's modulus, residual stress and Poisson's ratio of the ith layer film are b_i, h_i, ρ_i, E_i, σ_i, ѵ_i, respectively. Here a positive stress represents the tensile and a minus stress represents the compressive.

Zoom In Zoom Out Reset image size

In the designed test structures, since the width of the beam is bigger than five times thickness of the beam, the beam is considered to be wide, and Young's modulus E_i is replaced by the effective Young's modulus $\overline{{{E}_{i}}}$, i.e.

$$\begin{eqnarray}&&\overline{{{E}_{i}}}=\frac{{{E}_{i}}}{1-\nu _{i}^{2}}\end{eqnarray} \tag{ 1 }$$

In this study, the effective Young's modulus $\overline{{{E}_{i}}}$ is used instead of the Young's modulus E_i. If the height direction of the beam is assumed to be along z axis, and the position of the beam bottom, z₀, is 0 while the position of the top of the ith layer is z_i and the distance between the neutral axis and the bottom of the beam is z_c. It yields

$$\begin{eqnarray}&&{{z}_{i}}=\underset{j=1}{{i \sum}} {{h}_{i}},\quad \left(i=1,2,\cdots ,n\right).\end{eqnarray} \tag{ 2 }$$

According to the static equilibrium condition of the forces, z_c can be written as [16]

$$\begin{eqnarray}&&{{z}_{c}}=\frac{\underset{i-1}{\overset{n}{\sum}} \overline{{{E}_{i}}}{{b}_{i}}\left(z_{i}^{2}-z_{i-1}^{2}\right)}{2\underset{i-1}{\overset{n}{\sum}} \overline{{{E}_{i}}}{{b}_{i}}\left({{z}_{i}}-{{z}_{i-1}}\right)}\end{eqnarray} \tag{ 3 }$$

Then the moment of inertia of ith layer, I_i, with respect to the neutral axis of the beam is given by

$$\begin{eqnarray}&&{{I}_{i}}=\underset{{{z}_{i-1}}}{\overset{{{z}_{i}}}{\int}} {{z}^{2}}\text{d}{{A}_{i}}=\frac{1}{3}{{b}_{i}}\left(z_{i}^{3}-z_{i-1}^{3}\right),\quad \left(i=1,2,\cdots ,n\right),\end{eqnarray} \tag{ 4 }$$

where A_i is the cross-section area of the ith layer. For a multilayered beam, the bending stiffness, $\overline{EI}$, the axial stiffness, $\overline{EA}$, the linear density, $\overline{\rho A}$, and the axial load, $\overline{\sigma A}$, are defined by the sum of the individual terms for each layer in the beam as follows [6]

$$\begin{eqnarray}&&\overline{EI}=\underset{i=1}{\overset{n}{\sum}} \overline{{{E}_{i}}{{I}_{i}}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\overline{EA}=\underset{i=1}{\overset{n}{\sum}} \overline{{{E}_{i}}{{A}_{i}}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\overline{\rho A}=\underset{i=1}{\overset{n}{\sum}} {{\rho}_{i}}{{A}_{i}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\overline{\sigma A}=\underset{i=1}{\overset{n}{\sum}} {{\sigma}_{i}}{{A}_{i}}\end{eqnarray} \tag{ 8 }$$

As is well known [17], in the beam theory there is an nth critical buckling load ${{P}^{\left(n\right)}}$defined as the axial force sufficient to maintain the beam in a slightly bent form. When the axial load P is greater than the critical value${{P}^{\left(n\right)}}$, the beam remains straight. And if $P<{{P}^{\left(n\right)}}$ (i.e., $\bigg|P\bigg|>\bigg|{{P}^{\left(n\right)}}\bigg|$), the beam is initially statically deflected and a postbuckling displacement $b{{\varphi}_{n}}\left(\hat{x}\right)$ is induced, where b is the dimensional scaling constant and ${{\varphi}_{n}}\left(\hat{x}\right)$ is the nth buckling mode shape. These two cases will be discussed below.

2.1. Initially flat beam

For an unbuckled beam, the axial load P = $\overline{\sigma A}$ is greater than the critical value ${{P}^{\left(n\right)}}$, so the beam remains flat. For the undamped beam with a uniform cross section along its length, the natural frequencies ω for the initially flat beam is given by the eigenvalues [18]

$$\begin{eqnarray}\underbrace {\left| {\begin{array}{*{20}{c}} 0&1&0&1\\ {{\lambda _1}}&0&{{\lambda _2}}&0\\ {\sin ({\lambda _1}l)}&{\cos ({\lambda _1}l)}&{\sinh ({\lambda _2}l)}&{\cosh ({\lambda _2}l)}\\ {{\lambda _1}\cos ({\lambda _1}l)}&{ - {\lambda _1}\sin ({\lambda _1}l)}&{{\lambda _2}\cosh ({\lambda _2}l)}&{{\lambda _2}\sinh ({\lambda _2}l)} \end{array}} \right|}_{\bf{\rm{S}}} = 0 \end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{\beta}_{1}}=\frac{\overline{\sigma A}}{\overline{EI}},\\ &&\nonumber {{\beta}_{2}}=\frac{\overline{\rho A}}{\overline{EI}}{{\omega}^{2}},\\ &&\nonumber {{\lambda}_{1}}={{\left(\frac{{{\left(\beta _{1}^{2}+4{{\beta}_{2}}\right)}^{1/2}}-{{\beta}_{1}}}{2}\right)}^{1/2}},\\ &&\nonumber {{\lambda}_{2}}={{\left(\frac{{{\left(\beta _{1}^{2}+4{{\beta}_{2}}\right)}^{1/2}}+{{\beta}_{1}}}{2}\right)}^{1/2}} \end{eqnarray} \tag{ 10 }$$

Here the matrix in equation (9) is considered to be S.

2.2. Buckled beam

For the buckled beam, $P<{{P}^{\left(n\right)}}$, the beam will buckle and a postbuckling displacement $b{{\varphi}_{n}}\left(\hat{x}\right)$ is induced. The natural frequencies ω for the buckled beam is given by the eigenvalues [17, 19]

$$\begin{eqnarray}\underbrace {\left| {\begin{array}{*{20}{c}} 0&1&0&1&1\\ {{\lambda _1}}&0&{{\lambda _2}}&0&0\\ {\sin ({\lambda _1}l)}&{\cos ({\lambda _1}l)}&{\sinh ({\lambda _2}l)}&{\cosh ({\lambda _2}l)}&1\\ {{\lambda _1}\cos ({\lambda _1}l)}&{ - {\lambda _1}\sin ({\lambda _1}l)}&{{\lambda _2}\cosh ({\lambda _2}l)}&{{\lambda _2}\sinh ({\lambda _2}l)}&0\\ {\frac{{\frac{{\alpha {\lambda _1}}}{{{l^2}}}(\cos ({\lambda _1}l) - 1)}}{{\lambda _1^2 - {{(\frac{{2\pi }}{l})}^2}}}}&{\frac{{\frac{{ - \alpha {\lambda _1}}}{{{l^2}}}\sin ({\lambda _1}l)}}{{\lambda _1^2 - {{(\frac{{2\pi }}{l})}^2}}}}&{\frac{{\frac{{\alpha {\lambda _2}}}{{{l^2}}}(1 - \cosh ({\lambda _2}l))}}{{\lambda _2^2 + {{(\frac{{2\pi }}{l})}^2}}}}&{\frac{{\frac{{ - \alpha {\lambda _2}}}{{{l^2}}}\sinh ({\lambda _2}l)}}{{\lambda _2^2 + {{(\frac{{2\pi }}{l})}^2}}}}&{\frac{{{\beta _2} - \frac{{\pi \alpha }}{l}}}{{2\pi }}} \end{array}} \right|}_{\bf{B}} = 0\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&\alpha =\frac{8{{\pi}^{3}}{{b}^{2}}}{{{l}^{3}}}\frac{\overline{EA}}{\overline{EI}}\\ &&\nonumber {{\beta}_{1}}=\frac{{{P}^{\left(1\right)}}}{\overline{EI}},\\ &&\nonumber {{\beta}_{2}}=\frac{\overline{\rho A}}{\overline{EI}}{{\omega}^{2}},\\ &&\nonumber {{\lambda}_{1}}={{\left(\frac{{{\left(\beta _{1}^{2}+4{{\beta}_{2}}\right)}^{1/2}}-{{\beta}_{1}}}{2}\right)}^{1/2}},\\ &&\nonumber {{\lambda}_{2}}={{\left(\frac{{{\left(\beta _{1}^{2}+4{{\beta}_{2}}\right)}^{1/2}}+{{\beta}_{1}}}{2}\right)}^{1/2}}. \end{eqnarray} \tag{ 12 }$$

Here the matrix in equation (11) is considered to be B.

2.3. Extraction of material properties

The matrices in equations (9) and (11) are defined as S and B respectively, equations (9) and (11) can be rewritten as

$$\begin{eqnarray}&&{{F}_{\text{S}}}=\left(\omega ,\overline{{{E}_{1}}},\overline{{{E}_{2}}},\bullet \bullet \bullet ,\overline{{{E}_{n}}},{{\sigma}_{1}},{{\sigma}_{2}},\bullet \bullet \bullet ,{{\sigma}_{n}}\right)=\det \left(\text{S}\right)=0\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{F}_{\text{B}}}=\left(\omega ,\overline{{{E}_{1}}},\overline{{{E}_{2}}},\bullet \bullet \bullet ,\overline{{{E}_{n}}},{{\sigma}_{1}},{{\sigma}_{2}},\bullet \bullet \bullet ,{{\sigma}_{n}}\right)=\det \left(\text{B}\right)=0\end{eqnarray} \tag{ 14 }$$

where $\overline{{{E}_{1}}},\overline{{{E}_{2}}},\bullet \bullet \bullet ,\overline{{{E}_{n}}},{{\sigma}_{1}},{{\sigma}_{2}},\bullet \bullet \bullet ,{{\sigma}_{n}}$ represent the effective moduli and residual stresses of each layer, respectively, and $\omega$ is the natural frequency. Since equations (13) and (14) indicate that the natural frequencies are dependent of the material properties and geometries of the beam, the Young's modulus and residual stress of each layer of the composite beam can be extracted by using natural frequencies of a set of beams with different geometries. For an n-layer composite beam, at least two kinds of beams with different lengths are needed, and each kind of beam with the same length needs at least n different types in widths. If the widths of each layer in the jth beam are assumed to be ${{b}_{j1}},{{b}_{j2}},\cdots ,{{b}_{jn}}\quad \left(1\leqslant j\leqslant n\right),$ and those in the kth beam with the same length as the jth beam are ${{b}_{k1}},{{b}_{k2}},\cdots ,{{b}_{kn}}\quad \left(1\leqslant k\leqslant n, k\ne j\right)$, it is required that the vector $\left({{b}_{j1}},{{b}_{j2}},\cdots ,{{b}_{jn}}\right)$ and vector $\left({{b}_{k1}},{{b}_{k2}},\cdots ,{{b}_{kn}}\right)$ must be linearly independent. It is well known that for a single-layer beam when n equals to 1, it is not restricted by this principle. So at least 2n kinds of beams are needed. The natural frequencies of these beams are required, and deflection profiles of these beams need to be observed so as to determine equation (13) or equation (14) to be used. Once the natural frequencies of the designed test structures are measured, the Young's modulus and residual stress of each layer for the composite beam can be extracted by numerical methods such as the Newton iteration method or least-square method. Generally, there are infinite solutions (Young's modulus E and residual stress σ) to the above equation. However, for an engineering problem, for example, the Young's modulus of the polysilicon has been reported in the range of 120 Gpa~201 GPa. Therefore, only one solution will be in this range. This solution is what we want. For the bilayered beam there are four unknowns (two Young's moduli and two residual stresses), four beams with different (length and width) dimensions are needed. If it is assumed that the Young's modulus and residual stress are the same for the each layer of the four beams because they are located close to each other, we have

$$\begin{eqnarray}&&\bigg\{\begin{array}{c} {{F}_{11}}=\left({{\omega}_{11}},\overline{E1},\overline{E2},\sigma 1,\sigma 2\right)=0 \\ {{F}_{21}}=\left({{\omega}_{21}},\overline{E1},\overline{E2},\sigma 1,\sigma 2\right)=0 \\ {{F}_{12}}=\left({{\omega}_{12}},\overline{E1},\overline{E2},\sigma 1,\sigma 2\right)=0 \\ {{F}_{22}}=\left({{\omega}_{22}},\overline{E1},\overline{E2},\sigma 1,\sigma 2\right)=0 \\ \end{array}.\end{eqnarray} \tag{ 15 }$$

In the experiments, we measured and obtained four different fundamental frequencies (ω₁₁, ω₂₁, ω₁₂, ω₂₂) for the four beams by a scanning LDV. Then, we can extract the two Young's moduli and two residual stresses for the bilayered beam according to the solving flow. In the experiments, by using a digital holographic microscope (DHM, R2200), initially flat or buckled beams were checked so that equation (13) for the initially flat or equation (14) for the buckled beams could be selected.

To verify the proposed approach, analysis of doubly clamped bilayer beams as an example were simulated by FEM. Figure 2 shows a model for the doubly clamped bilayer beam, where the bottom layer is polysilicon and the top layer is gold. Four different beams in length and width are taken so that the Young's moduli and residual stresses of both layers can be extracted. Input parameters of the beams for the simulation are listed in table 1. To distinguish between the initially flat and buckled conditions, the simulations are carried out three times with two different pairs of residual stresses acting on the beams. Using the parameters in table 1, the fundamental frequencies and deflection profiles of the beams are obtained in FEM simulations by Ansys. The results are then put back into equations (13) or (14) to figure out the Young's moduli and residual stresses. Table 2 lists the comparison between materials properties set to the FEM analysis and those calculated from FEM analysis. It is clear from table 2 that the relative errors are less than 23%.

Zoom In Zoom Out Reset image size

Test structures were fabricated using MEMSCAP PolyMUMPs process. Layer Poly2 and layer Metal were used as structural materials. Both bilayer and single-layer doubly clamped beams were fabricated. For the bilayer beams, two kinds of lengths were applied while for the beams with the same length but different widths of layer Metal. And additional single-layer beams were also provided, which can be treated as the 0 width of layer Metal. Parameters of the test structures are listed in table 3. Scanning electron microscopy (SEM) photographs of the test structures are shown in figure 3. All the measurements were performed at room temperature.

Zoom In Zoom Out Reset image size

To determine the deformation of the beams, a DHM (R2200) from Lyncée Tec was used. The setup of the DHM system is shown in figure 4. The DHM can produce an off-axis hologram through the interference between the wave reflected by the sample and a reference wave. Figure 5 shows a 3 D image of one buckled beam. The profiles of the beams extracted from the DHM measurement are depicted as shown in figure 6. It can be observed clearly that some profiles of the beams are similar to a cosine curve, while the other profiles gently bend like bows. Due to the step-up support of the beams introduced by the surface micro-machined fabrication process, the doubly clamped beams applied with a slightly compressive load before buckling will have a small deformation rather than keep straight [21]. Hence the beams with profiles like cosine curves are considered as buckled and others slightly bent are treated as initially flat beams. The observation results are summarized in table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After the deformation profiles were tested, resonant frequencies of the beams were then measured. A Polytech scanning LDV system with hardware modules of MSV-400-M2 was used. The setup of the LDV system is shown in figure 7. The test structures were excited by a piezoelectric plumbum zirconate titanate (PZT) ceramic slice with a periodical chip signal. By applying driven signal with frequencies ranging from 1 kHz to 1 MHz, the resonant frequencies were extracted from the frequency response functions [22]. Figure 8 shows the frequency response curve of a bilayer beam with length of 300 μm. The fundamental frequencies of six beams were measured. The measured fundamental frequencies are listed in table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With the measured results, the effective moduli and residual stresses were extracted simultaneously by equations (13) or (14). There are 12 types of combinations in total, made up of four beams out of the whole six ones, which forms 12 groups of solvable equations from equations (13) or (14). The Young's moduli and residual stresses of both layers are then extracted. The beam combinations and calculated results are listed in table 4. The results show the proposed method produces results that agree with the values most often given in the literature. For example, the Young's modulus of gold in this experiment coincides with 78 GPa as reported in [23].

There are some factors that may affect the accuracy of the extraction model. First of all, squeezed damping in air may lead to the measured resonant frequency to shift from the expected value [22]. To reduce this deviation, the amplification factor for the frequency response functions at points of the resonant frequencies are required to be more than 20, and in this situation the error caused by the frequency shift will be less than 0.01%. Secondly, compliance of step-up support fabricated by the surface micro-machined fabrication process may cause a relaxation of residual stresses. It is reported that under a compressive biaxial residual stress of 11 MPa, the relaxed stress is about 9 MPa as beam length is 1000 μm, while the relaxed stress decreases to 3 MPa as beam length decreases to 500 μm [25]. In this study, lengths of doubly clamped beams are less than 400 μm. Hence, the relaxation of residual stress in the test structures is small and acceptable.

In this paper, an approach is presented to measure the Young's modulus and residual stress of an individual layer for MEMS composite films. A theoretical model for extracting material properties by resonant frequencies of doubly clamped multilayered beams is analyzed. Material properties for the composite films can be obtained by a set of doubly clamped beams with different lengths and widths. Both initially flat and buckled beams have also been considered. The validity of the proposed method has been confirmed by the finite element analysis and experiments. This method is suitable for monitoring the MEMS fabrication process.

This project is partly supported by the National Science and Technology Major Project under contract no. 2011ZX02507.