An Analytical Model for CMUTs with Square Multilayer Membranes Using the Ritz Method

Capacitive micromachined ultrasonic transducer (CMUT) multilayer membrane plays an important role in the performance metrics including the transmitting efficiency and the receiving sensitivity. However, there are few studies of the multilayer membranes. Some analytical models simplify the multilayer membrane as monolayer, which results in inaccuracies. This paper presents a new analytical model for CMUTs with multilayer membranes, which can rapidly and accurately predict static deflection and response frequency of the multilayer membrane under external pressures. The derivation is based on the Ritz method and Hamilton’s principle. The mathematical relationships between the external pressure, static deflection, and response frequency are obtained. Relevant residual stress compensation method is derived. The model has been verified for three-layer and double-layer CMUT membranes by comparing its results with finite element method (FEM) simulations, experimental data, and other monolayer models that treat CMUTs as monolayer plates/membranes. For three-layer CMUT membranes, the relative errors are ranging from 0.71%–3.51% for the static deflection profiles, and 0.35%–4.96% for the response frequencies, respectively. For the double-layer CMUT membrane, the relative error with residual stress compensation is 4.14% for the central deflection, and −1.17% for the response frequencies, respectively. This proposed analytical model can serve as a reliable reference and an accurate tool for CMUT design and optimization.


Introduction
Capacitive Micromachined Ultrasonic Transducers (CMUTs) were introduced to the ultrasound community about three decades ago [1,2]. Compared to traditional piezoelectric transducers, CMUTs exhibit some attractive features such as lower mechanical impedance [3], broader bandwidth [4], lower sensitivity to temperature [5], stable device properties, ease of fabrication, and compatibility with Micro-Electro-Mechanical Systems (MEMS) devices [6]. Thus, they have become an excellent choice suitable for a wide range of applications including medical diagnostic imaging [7,8], nondestructive testing [9], and automotive collision avoidance applications such as parking assistance or blind-spot monitoring [10]. In order to facilitate the design and optimization of CMUTs, numerous investigations have been directed towards understanding, predicting and controlling their static and dynamic characterizations [11,12].
A typical CMUT is built with a square, rectangular, circular, or hexagonal membrane separated from a fixed substrate by a small airgap [13,14]. The working principles are illustrated in Figure 1. The vibrating membrane can generate or detect ultrasonic waves as a transmitter or a receiver. During transmit-mode operation, a DC voltage is applied between the two electrodes. The membrane is membrane is attracted to the fixed substrate by an electrostatic force. When an AC voltage is superimposed over the DC voltage, the membrane will move in response to the applied signal, generating an ultrasonic wave that is launched into the ambient environment. During the receive-mode operation, the CMUT membrane is subjected to an ultrasonic wave. The membrane deformation will occur and produce a capacitance variation. A microelectronic circuit is used to detect the capacitance variation and process this into a useful signal. Because the CMUT's performances depend on the change of the deflection before and after that the membrane is biased, accurate calculation of the deflection is crucial [15]. Besides, the response frequency, as an important dynamic characteristic, determines the specific application and operational bandwidth of the device [16]. Therefore, we have chosen the static deflection and the response frequency as the research focuses of this paper.
Since the exact shape of a deformed clamped membrane is not known, generalized plate theories have been applied by many authors to obtain a functional form of the membrane deformation curve. The deflection calculation must satisfy the boundary conditions, membrane geometry, and the specific loading condition. As is known that the boundary condition depends on specific membrane geometry (e.g., square, rectangular, circular, or hexagonal), the load deflection model and deflection functions are unique for each of the cases and each demands separate treatment. For clamped square and rectangular membranes, no simple exact deflection solution exists and approximate methods must be used. In this paper, only the square membrane has been considered.
The first CMUTs were built using a sacrificial release process that has become the standard CMUT fabrication method. Numerous variations of the sacrificial release process have been published, all based on the same basic principle [2,5]. The wafer-bonding method is another widely-used CMUT fabrication technique nowadays, based on a different approach to cavity formation that uses a combination of bulk and surface micromachining techniques [2,13]. Commonly, the membrane is fabricated directly using a thin conduction film (e.g., aluminum or polysilicon) or a composite of a non-conduction film (e.g., silicon nitride) with a conduction coating on top (e.g., aluminum or gold). To avoid electrical breakdown after the pull-in phenomenon, an insulation layer is added, either under the membrane or on top of the fixed substrate. Additionally, a passivation layer is fabricated on top of the membrane to protect the CMUT from environmental contamination [13]. The fabrication process described above determines the multilayer structure of the CMUT membrane.
Researchers have described a variety of approaches to determine the deflection curve of single-layer clamped square membranes used in MEMS devices. Timoshenko and Woinowsky-Krieger presented an approximate mathematical expression to determine the deflection curve for thin plate in the small deflection regime using a trigonometric series [17]. However, the model is computationally expensive because it requires extensive numerical calculations to determine a set of coefficients. Bao provided the mechanical analysis of the square membrane, from which the deflection profile and the vibration frequency could be obtained [6]. Ben Moussa used the Galerkin Because the CMUT's performances depend on the change of the deflection before and after that the membrane is biased, accurate calculation of the deflection is crucial [15]. Besides, the response frequency, as an important dynamic characteristic, determines the specific application and operational bandwidth of the device [16]. Therefore, we have chosen the static deflection and the response frequency as the research focuses of this paper.
Since the exact shape of a deformed clamped membrane is not known, generalized plate theories have been applied by many authors to obtain a functional form of the membrane deformation curve. The deflection calculation must satisfy the boundary conditions, membrane geometry, and the specific loading condition. As is known that the boundary condition depends on specific membrane geometry (e.g., square, rectangular, circular, or hexagonal), the load deflection model and deflection functions are unique for each of the cases and each demands separate treatment. For clamped square and rectangular membranes, no simple exact deflection solution exists and approximate methods must be used. In this paper, only the square membrane has been considered.
The first CMUTs were built using a sacrificial release process that has become the standard CMUT fabrication method. Numerous variations of the sacrificial release process have been published, all based on the same basic principle [2,5]. The wafer-bonding method is another widely-used CMUT fabrication technique nowadays, based on a different approach to cavity formation that uses a combination of bulk and surface micromachining techniques [2,13]. Commonly, the membrane is fabricated directly using a thin conduction film (e.g., aluminum or polysilicon) or a composite of a non-conduction film (e.g., silicon nitride) with a conduction coating on top (e.g., aluminum or gold). To avoid electrical breakdown after the pull-in phenomenon, an insulation layer is added, either under the membrane or on top of the fixed substrate. Additionally, a passivation layer is fabricated on top of the membrane to protect the CMUT from environmental contamination [13]. The fabrication process described above determines the multilayer structure of the CMUT membrane.
Researchers have described a variety of approaches to determine the deflection curve of single-layer clamped square membranes used in MEMS devices. Timoshenko and Woinowsky-Krieger presented an approximate mathematical expression to determine the deflection curve for thin plate in the small deflection regime using a trigonometric series [17]. However, the model is computationally expensive because it requires extensive numerical calculations to determine a set of coefficients. Bao provided the mechanical analysis of the square membrane, from which the deflection profile and the vibration frequency could be obtained [6]. Ben Moussa used the Galerkin method with a polynomial basis function to determine the deformation curve [18]. This model is limited in accuracy and exhibits slow convergence. Kerrour and Hobar improved the accuracy and convergence of the approach in [19] by replacing the polynomial basis function with a trigonometric basis function [19]. However, the model has not been verified against finite element method (FEM) or experimental results. Elgamel proposed a deflection shape function using a single-term Fourier approximation of the exact bending shape that incorporates a double cosine square term [20]. Because of its simpler form, this function is widely used for load-deflection analysis of clamped square membranes subject to large deflection. However, its accuracy is compromised by the truncation of higher order terms in the Fourier series. Rahman et al. developed a deflection shape function by modifying a polynomial-based function presents in [18] to include two empirically determined parameters to improve the accuracy, and the analytical results were compared with FEM and experiments for the single-layer membrane [15].
However, there is a lack of research focusing on multilayer characterization. Analytically, many current derivations treat the multilayer membrane as a monolayer plate/membrane, and are dependent mainly on the properties of the device layer [2,16,21,22]. In fact, when the multilayer membrane is under the action of external forces, each layer suffers deformation and internal forces developed between each layer. In such cases, material coupling will occur due to transverse bending and in-plane stretching. Obviously, such simplification will lead to inaccuracies and errors in the strain-displacement equations of the membrane. FEM has also been widely used to understand transducer characteristics and to optimize transducer response, since FEM results are highly accurate and reliable. However, in most applications, the thickness of the CMUT's membrane is small in comparison with its lateral dimensions. Since the FEM mesh size is limited to a minimum size, the finite element simulation may have a huge mesh generation with a poor element quality. Also, the simulation could take a long time to make a single change in one parameter. Since most CMUTs have regular shapes with limited boundary conditions, finite element simulations are not efficient for design optimization. In contrast, analytical models can provide better insight into the multilayer characteristics. Therefore, it is necessary to establish an accurate analytical model for the multilayer characterizations of CMUTs.
In this paper, we presented a new analytical model which focuses on the CMUT's multilayer characterization. The proposed model offers (1) a new, readily usable, simple, and accurate deflection shape function for uniformly loaded clamped square multilayer membrane used in typical CMUT design; (2) a response frequency calculation model derived under the external pressure variation; and (3) a residual stress compensation method based on energy functional using the Ritz method. The rest of the paper has been organized in the following order: Section 2 derives the proposed energy functional model and the compensation method based on the Ritz method and Hamilton's principle. The relationship between the external pressure, the static deflection, and the response frequency is obtained Section 3 provides the experimental and FEM validation of the model for both the three-layer and double-layer CMUT membranes. Analytical derivations with and without the residual stress compensation have also been compared for the double-layer CMUT membrane. Section 4 concludes the paper and assesses the validity of the proposed model. Figure 2 shows a schematic drawing of a multilayer membrane used in CMUTs. The vibrating section is a clamped square N-layer membrane that is excited by an electrostatic force. The origin O of the Cartesian coordinate system is in the center, and the xoy surface is in the middle layer (the "midsurface") of the multilayer membrane. 2a and h are the width and the total thickness. The top and bottom surfaces of the membrane are located at h/2 and´h/2 in the z direction. The ith layer is located between z i and z i+1 . Here, E i , ρ i , and µ i are Young's Modulus, the density, and Poisson's ratio of the ith layer, respectively. of the Cartesian coordinate system is in the center, and the xoy surface is in the middle layer (the "midsurface") of the multilayer membrane. 2a and h are the width and the total thickness. The top and bottom surfaces of the membrane are located at h/2 and −h/2 in the z direction. The ith layer is located between zi and zi+1. Here, Ei, ρi, and μi are Young's Modulus, the density, and Poisson's ratio of the ith layer, respectively.

Static Deflection Analysis
Under Kirchhoff's hypothesis, the inplane displacements u, v and the transverse deflection w of the membrane in the x, y, and z directions are: where u 0 px, yq, v 0 px, yq, and w 0 px, yq are the displacement components in the midsurface in the x, y, and z directions. λ x and λ y are the rotations of the midsurface about the x/y axis given by: The strain-displacement relations of the membrane are: where ε 0 x , ε 0 y , and ε 0 xy are the reference surface strains at z = 0 defined by: and k x , k y , and k xy are the membrane curvatures given by: Micromachines 2016, 7, 55 5 of 20 The nonlinear strain-displacement relations given by Equation (4) are those of von Kármán [23]. Due to the constitutive equations, the total stress components for the ith layer in Cartesian coordinates are: The strain energy for the ith layer is: where V i is the volume of the ith layer. Substituting Equations (3) and (6) into Equation (7) gives: where S is the area of the ith layer, as shown in Figure 2. The total potential energy Π 1 of the multilayer membrane as an elastic system is: where U is the total strain energy, and W is the potential energy of the uniform external load p to the N-layer membrane.
Since it is extremely difficult to directly solve the nonlinear Equation (9), the Ritz method [24] can be used to approach the analytical approximation under the CMUT's first order frequency. Both the Galerkin method [25,26] and the Ritz method are effective methods to solve elastic mechanical problems, yet they have different approaches. The Ritz method is a kind of the energy method and based on the principle of least potential energy, while the Galerkin method is a special form of the weighted residual method. For conservative elastic systems whose functional exist, the Ritz method is more efficient and practical. The Ritz method only requires the trail solution to meet the constraint boundary conditions, while the Galerkin methods also requires natural boundary conditions. In our case, the membrane shape is regular, the elastic system is conservative, the functional exist, and the trail solution is available. So we chose the Ritz method over the Galerkin method.
For a clamped square membrane, the boundary conditions are: Thus, the approximate solutions for the static displacement components can be written as: 3 V m , and W m are the ratios of the maximum static displacement components in the x, y, and z directions to the total thickness h. Equation (11) is then substituted into Equation (9). According to the Ritz method, the energy expression Π 1 is a function of three coefficients whose numerical values can be determined from the conditions that: Equation (12) yields the following set of algebraic equations: For the cubic equation in W m , the coefficients B 1 -B 3 are: where: where the coefficients α ij´i " 1, . . . , 6 j " 1, . . . , 6¯are elements in the compliance matrix of the multilayer membrane. For our case of an isotropic multilayer membrane clamped at its periphery, the compliance matrix can be reduced to a symmetric matrix with six independent elastic constants listed in Equation (15), and other elements are zero [24]. Using the material properties of the CMUT's layers, Equation (15) can be calculated to determine the static displacements in all three directions.

Response Frequency Analysis
Similarly, the time-dependent vibrating three-dimensional displacements are given by: where ω is the angular frequency of the vibrating multilayer membrane. u t 0 px, yq, v t 0 px, yq, and w t 0 px, yq are the vibrating displacement components in the midsurface in the x, y, and z directions. λ t x and λ t y are the vibrating rotations of the midsurface about the x/y axis given by Thus the total displacements in the x, y, and z directions can be written as the summation of the static and the vibrating displacement components: The strain energy U T i and the kinetic energy T i of the ith layer now are: Micromachines 2016, 7, 55 8 of 20 where: Then the total potential energy Π 2 of the multilayer membrane is: According to Hamilton's principle [27], Equation (22) should satisfy that for arbitrary t 1 and t 2 , Π 2 attains its minimum value. The static deflections u 0 , v 0 , and w 0 obtained from Equation (9) can be used for the calculations of the vibrating displacements u t 0 , v t 0 , w t 0 and their angular frequency ω. Again, the Ritz method is utilized to evaluate the desired extreme value of Π 2 .
Given the boundary conditions of a clamped square membrane as in Equation (10), the approximate solutions for the dynamic displacement components can be written as: and W t are the ratios of the maximum vibrating displacement components in the x, y, and z directions to the total thickness h.
Together with the Ritz method and Hamilton's principle, the energy functional Π 2 satisfies the following conditions: where: where β ij´i " 1, 2, 3 j " 1, 2, 3¯and γ ij´i " 1, 2, 3 j " 1, 2, 3¯are the reduced stiffnesses and flexural rigidities in the plane-stress constitutive matrix of the multilayer membrane, respectively. According to matrix Equation (25), if the operation frequency of the multilayer CMUT membrane is known, the relevant material and geometric parameters can be obtained. Similarly, the frequency of the multilayer membrane can be acquired by using the known material and geometric parameters. In order to acquire the nonzero solution of the matrix Equation (25), the matrix determinant should be zero. During this calculation process, the angular frequency ω of the CMUT membrane can be obtained.
Note that for an isotropic material having elastic symmetry with respect to the midsurface, the values of α 13 and γ 13 for the static and dynamic displacement calculations are zero, meaning that the coupling rigidities between transverse bending and in-plane stretching will disappear. Thus, it can be predicted that our proposed model will be sensitive to symmetric geometries.

Residual Stress Compensation
According to the plate theory [17], geometries can be divided into a "thin plate" category and a "thin membrane" category, according to their dimension ratios of thickness to minimum lateral dimensions. Between them two, thin membranes suffer more from the residual stress. For the isotropic thin plate without residual stress, the method mentioned above is fully capable to acquire the accurate solutions. For anisotropic thin plates, the stiffness matrix or the compliance matrix needs to be modified [25] (which will not be discussed in this paper). For geometries with residual stress, compensation needs to be added.
As one of the most important mechanical characterizations, residual stress in material is influential and inevitable. During the fabrication process, residual stress has been produced. It can be empirically predicted within a certain range, but the exact value cannot be predicted. On the contrary, it has to be measured by specialized instruments, and the measuring methods include the X-ray, the Raman spectrum, the infrared analysis, and the electron diffraction.
When the CMUT membrane vibrates, the residual stress works as a pre-tension. In order to obtain the accurate predictions of the static deflection and the frequency response, residual stress compensation has to be considered.
CMUT membranes can be regarded as isotropic thin plates/membranes, so the strain in the thickness direction can be ignored (ε 0 z i " 0), and the residual stress can be treated as plane stress. Due to the presence of the residual stress, an additional elastic potential energy can be created during the membrane's vibration. Let the residual stress for the ith layer be´σ 0 x i , σ 0 y i , σ 0 xy i¯, and the vibration displacement be: V px, y, z, tq " u px, y, z, tq i`v px, y, z, tq j`w px, y, z, tq k where u px, y, z, tq , v px, y, z, tq , w px, y, z, tq are the displacement components in x, y, z directions, and i, j, k are the unit vectors in x, y, z directions. The strain-displacement relationship of the membrane can be defined by: So the additional elastic potential energy due to the residual stress can be expressed as: For actual devices, the residual stress can be measured, and the additional elastic potential energy can be calculated using Equation (30). The energy functional for CMUT's static and vibrating displacements can be modified as: Through similar calculation processes mentioned above, the static deflection and frequency response of the CMUT membrane with the presence of residual stress will be obtained.
This paper focuses on the static deflection and the frequency response of the CMUT's multilayer membrane. In the next section, Equation (11) will be used for the central deflection and the deflection profiles under a variety of external pressures. Equation (25) will be used to calculated the frequency of the CMUT's multilayer membrane by f " ω 2π . For actual devices, the residual stress compensation will be done using Equations (28)-(31).

Validation and Discussion
In order to validate the theoretical analyses obtained above, their results were compared with those from simulations, previous literatures, and experiments. Due to the fact that the square membrane cannot be simplified as 2D axisymmetric models as the circular membranes, 1/4 substructure 3D FEM models were constructed. The cubic mesh was chosen for our case, since it gives higher mesh quality for a reduced computation time.
This section is organized as follows: Firstly, three-layer CMUT membranes are analyzed and modeled by the proposed model, simulations, and previous literatures. The analytical predictions of the static deflection profiles and the frequency responses are compared with simulations for different dimensions under different external pressure values.
Secondly, a double-layer CMUT membrane is designed, fabricated, and tested. The static deflection profiles and the frequency response under the standard ambient pressure are measured, using the Sensofar 3D Optical Profiler (Sensofar, Barcelona, Spain) and the Agilent Precision Impedance Analyzer (Agilent Technologies Inc., Palo Alto, CA, USA). The analytical results predicted by the proposed model with and without the residual stress compensation, the other literatures, and the simulations are compared with the experimental data.

Three-Layer Membrane Case
The three-layer membranes were chosen for their prevalence in the MEMS devices [1,14]. The geometries and dimensions of the three-layer CMUT membranes are listed in Table 1. Due to different fabrication processes, CMUT's device layer can be either conductive or non-conductive. For CMUTs A and B, the 1st layer (SiN x ) is the insulation layer, the 2nd layer (Al) serves as both the device layer and the electrode, and the 3rd layer (SiN x ) is the passivation layer. For CMUTs C and D, the 1st layer (SiN x ) is the device layer, the 2nd layer (Au) is the top electrode, and the 3rd layer (SiN x ) is the passivation layer. Note that CMUTs A and B share the same geometry, and their three layers are symmetrical about the midsurface. CMUTs C and D share the other geometry, which is not symmetrical about the midsurface. CMUTs A and C have the same side-length, while B and D have the same side-length. All layers within CMUTs A-D membranes are supposed to be ideally homogenous, isotropic, and without residual stress.
The material properties needed for calculations and simulations are listed in Table 2.

. Static Deflection Analysis
The CMUT's multilayer membrane can be deflected by the external pressure. The central deflections and the deflection profiles under the external pressure are calculated by our model, Bao [6], Rahman [15], and simulated by FEM. The vertical displacement distribution contours of CMUTs A and C for 0.5 MPa obtained from FEM are shown in Figure 3. The deflection profiles of the three-layer CMUT membranes are provided in Figure 4. The deflection profiles from the membrane center along the x-axis are chosen for comparison. Static deflection error analyses of three-layer CMUT membranes for different external pressures are listed in Table 3.    According to plate theory [18], the stiffness of a 70 µm-side-length membrane is bigger than that of a 120 µm-side-length membrane. Thus, for both analytical calculations and FEM simulations, the external pressures are 0.5/1 MPa for CMUTs A and C, and 0.1/0.2 MPa for CMUTs B and D.    The deflection distribution contours in Figure 3 indicate that for the same side-length, CMUTs A and C share a similar deflection tendency, and the largest vertical displacement is exhibited at the membrane center. However, due to the symmetrical geometry, A has a bigger central deflection and smoother deflection variation compared to C. When the external pressure is 0.5 MPa, A's central deflection is´0.207 µm, and C's is´0.186 µm. Also, for the same external pressure, a bigger average vertical displacement can be expected for A. For example, the vertical displacements between 15-25 µm for A is´0.15-´0.05 µm, while for C the range is about´0.12-´0.05 µm. Further discussion on this can be obtained by referring to Figure 4 and Table 3.
From Figure 4, it is worth noting that along the x-axis away from membrane center, the vertical displacements become smaller until they reach zero at the membrane edges in all cases. The proposed model's predictions are in strong agreement with the FEM simulations, while the other two show larger deviations. Data from Table 3 indicates that the relative errors of the center deflections between our model and FEM range from 0.71%-3.51%, which is fairly satisfactory. Therefore, it is quite necessary to treat CMUT's membranes as multilayer models, and our model for CMUT's multilayer membrane is much more accurate than other analytical models that treat the CMUT membranes as monolayer plates or membranes. Besides, monolayer models predict smaller deflections. Take CMUT A for 0.5 MPa, for example. Bao's and Rahman's models predict similar central deflections, which is around´0.15 µm. However, the FEM simulation is´0.207 µm (Figure 3), and the proposed model's prediction is´0.212 µm (Table 3). This deviation caused a relative error around 25%, which is not suitable for accurate determination of deflection profiles. Moreover, their predictions about the deflection profiles are similar, too, and Rahman's profile is smoother than Bao's due to the electrical coupling coefficients. Similar conclusions could be obtained from other comparisons for the same external pressure in Figure 4 and Table 3.
To make a further discussion on our model's validation, we focus on the comparisons between the proposed model's predictions and the FEM simulations.
From Figure 4a and Table 3, it can be observed that for CMUT A, when the external pressure is bigger, the relative error is smaller. At the membrane center (x = 0), our proposed model predicts slightly bigger vertical displacement than FEM. For the 0.5 MPa case, the relative error between our prediction and FEM is 2.41%. For the 1 MPa case, the relative error is reduced to 1.82%. Similar conclusions could be obtained from other sub-figures in Figure 4 and Table 3. Thus, the proposed model can maintain accurate predictions for a wide range of external pressures.
Take the deflection comparison between Figure 4a,c for 0.5 MPa, for example. It can be noticed that A has bigger vertical displacement than C, which was proved by Figure 3, too. Also, A's relative errors are bigger than C. For the 0.5 MPa case, A's relative error is 2.41%, while C's is 1.52%. Moreover, for the 1 MPa case, A's relative error is 1.82%, while C's is 1.05%. A and C have the same side-length but different geometry, yet their relative errors are similar and small. Thus, the proposed model can obtain stable predictions regardless of the membrane geometry. Comparisons between CMUTs B and D for the same external pressure can lead us to similar conclusions, too.
Comparisons between A and B in Figure 4a,b show that when the side-length is bigger, the relative error is slightly bigger within an acceptable range. For CMUT A, the relative errors are 1.82%-2.41%, while B's relative errors are 1.96%-3.51%. Due to the dimension error accumulation effect, this relative error deviation between A and B is acceptable. Given the ratios of the membrane thickness to the lateral dimension, both A and B can be categorized as "thin plates" [18]. Thus, the proposed model is suitable for "thin plates" category to obtain accurate deflection profiles. Similar results could be achieved by comparing CMUTs C and D.
For a clearer comparison, the absolute errors for the three-layer membranes in relation to the simulation data are shown in Figure 5. Together with Table 3, it can be seen that for symmetrical geometries (A and B), the absolute errors are ranging from 7.261-24.925 nm. For non-symmetrical geometries (C and D), the absolute errors are 6.649-25.265 nm. Besides, for the 70 µm-side-length (A and C) and the 120 µm-side-length three-layer membranes (B and D), the absolute errors are 6.649-15.712 nm, and 9.935-25.265 nm, respectively. These absolute errors are accurate and acceptable for CMUT design and modeling. In Figure 5, there is a peak in every curve. The peak is located around two-thirds of the half side-length away from the membrane center, indicating that the deflection shape function of the membrane can be further modified and improved.
Micromachines 2016, 7, 5 14 of 20 geometries (C and D), the absolute errors are 6.649-25.265 nm. Besides, for the 70 μm-side-length (A and C) and the 120 μm-side-length three-layer membranes (B and D), the absolute errors are 6.649-15.712 nm, and 9.935-25.265 nm, respectively. These absolute errors are accurate and acceptable for CMUT design and modeling. In Figure 5, there is a peak in every curve. The peak is located around two-thirds of the half side-length away from the membrane center, indicating that the deflection shape function of the membrane can be further modified and improved.

Frequency Response Analysis
Similarly, a frequency shift can be expected in response to external pressure variation. For the response frequency analysis, Bao's model [6] is taken as reference. The response frequencies for different pressures are calculated by the proposed model, Bao's model, and simulated by FEM. The response frequencies of the three-layer CMUT membranes are provided in Figure 6, and the error analyses are listed in Table 4. For the 70 μm-side-length membranes (CMUTs A and C), the external pressure range is 0-1 MPa. For the 120 μm-side-length membranes (CMUTs B and D), the external pressure range is 0.1-0.2 MPa.

Frequency Response Analysis
Similarly, a frequency shift can be expected in response to external pressure variation. For the response frequency analysis, Bao's model [6] is taken as reference. The response frequencies for different pressures are calculated by the proposed model, Bao's model, and simulated by FEM. The response frequencies of the three-layer CMUT membranes are provided in Figure 6, and the error analyses are listed in Table 4. For the 70 µm-side-length membranes (CMUTs A and C), the external pressure range is 0-1 MPa. For the 120 µm-side-length membranes (CMUTs B and D), the external pressure range is 0.1-0.2 MPa.
From Figure 6, it can be clearly observed that the proposed model's predictions are in strong agreement with the FEM simulations, while the predictions from Bao's model [6] show larger deviations. When the external pressures increase, the response frequencies increase, too. However, since the external pressures are not taken into consideration in Bao's model, the predictions from [6] are invariable, unlike the other two. Therefore, the proposed multilayer model for CMUT membrane is more accurate and suitable for response frequency determination.
Besides, compared to FEM simulations, the response frequencies of the monolayer models could be either smaller or bigger, due to different geometries. Take CMUT A, for instance. The response frequencies using FEM are 4.6703-4.8207 MHz, while Bao's prediction is 4.4799 MHz, which is smaller than the FEM simulations. On the other hand, for CMUT C featuring A's side-length but a different geometry, the Bao's prediction (3.9547 MHz) is higher than FEM simulations (3.3597-3.4868 MHz). Thus, the monolayer models like Bao's are not accurate enough for response frequency predictions for different external pressures. Similar conclusions could be obtained from comparisons between CMUTs B and D in Figure 6 and Table 4.

Frequency Response Analysis
Similarly, a frequency shift can be expected in response to external pressure variation. For the response frequency analysis, Bao's model [6] is taken as reference. The response frequencies for different pressures are calculated by the proposed model, Bao's model, and simulated by FEM. The response frequencies of the three-layer CMUT membranes are provided in Figure 6, and the error analyses are listed in Table 4. For the 70 μm-side-length membranes (CMUTs A and C), the external pressure range is 0-1 MPa. For the 120 μm-side-length membranes (CMUTs B and D), the external pressure range is 0.1-0.2 MPa.   According to Table 4, the relative errors between the proposed model and FEM for all cases are ranging from 0.35%-4.96%, which is quite satisfactory. To make a further discussion on our model's validation, we focus on the comparisons between the proposed model's predictions and the FEM simulations.
From Figure 6a and Table 4, it can be observed that when the external pressure increases, the response frequency and the relative error increase, too. When the external pressure is 0, the relative error between the proposed model and FEM is the smallest (1.31%). For the 1 MPa case, the relative error is increased to 2.75%, which is still within an acceptable range.
Take the response frequency comparison between Figure 6a,c, for example. A and C have the same side-length but different geometry, and A's relative errors (1.31%-2.75%) are slightly smaller than C's (1.77%-2.97%). Similar phenomenon can be observed from comparisons between CMUTs B and D. Thus, the proposed model can accurately predict the response frequencies of different geometries, and more accurate predictions can be expected for the symmetrical geometries.
Comparisons between CMUTs A and B show that when the external pressure is small, B's relative errors are smaller than A's; when the external pressure is big, B's relative errors are bigger than A's. CMUTs A and B have the same geometry but different side-length, yet their relative errors are within the same acceptable range. Similar conclusions can be obtained by comparing CMUTs C and D. Thus, the proposed model can produce accurate and stable predictions regardless of the dimensions.

Double-Layer Membrane Case
To verify the accuracy of the developed model for actual devices, a CMUT featuring a double-layer membrane has been designed, fabricated, and tested. The two layers are a 1 µm-thick, 70 µm-side-length square polysilicon membrane, and a 0.3 µm-thick, 70 µm-side-length square gold electrode located on top of the polysilicon membrane. The double-layer membrane is separated from the fixed substrate by a 1 µm-deep cavity. Using wafer-bonding technology, the CMUT was fabricated and sealed with a low-vacuum gap (10´2-10´3 mbar). The Young's modulus of polysilicon is 160 GPa, the Poisson's ratio is 0.22. The residual stress is around´200 MPa tested by the Renishaw inVia confocal Raman microscope, so in such plane stress condition, it can be obtained that σ 0 x i " σ 0 y i "´200 MPa and σ 0 xy i " 0 in the principal strain direction. These values will be used for the residual stress compensation in Equations (28)-(31). Inside the laboratory, the ambient pressure is around 1.013ˆ10 5 Pa, the temperature is 22˝C, and the granule concentration under 1 µm is less than 1000/m 3 , which are very suitable for delicate experiments. In this subsection, the analytical results from our proposed model with and without the residual stress compensation, the other literatures, and simulations are compared with the experimental data to verify the model's accuracy.

Static Deflection Analysis
The static deflection testing system is illustrated in Figure 7. The Sensofar 3D optical profiler (S neox Non-contact 3D Optical Profiler, SensoSCAN 5.3, Sensofar, Barcelona, Spain) is located on the Accurion active vibration isolation desktop unit (Accurion GmbH i4-OD-3173, Accurion GmbH, Goettingen, Germany). An image of the CMUT is obtained using a blue LED light source (the wavelength is 460 nm) for a high lateral resolution.

Static Deflection Analysis
The static deflection testing system is illustrated in Figure 7. The Sensofar 3D optical profiler (S neox Non-contact 3D Optical Profiler, SensoSCAN 5.3, Sensofar, Barcelona, Spain) is located on the Accurion active vibration isolation desktop unit (Accurion GmbH i4-OD-3173, Accurion GmbH, Goettingen, Germany). An image of the CMUT is obtained using a blue LED light source (the wavelength is 460 nm) for a high lateral resolution.  The deflection profiles of the double-layer CMUT membrane for 0.1 MPa predicted by the proposed model with and without the residual stress compensation (hereinafter referred as RSC), the other literatures, and the FEM simulations are compared with the experimental data, as shown in Figure 8. The deflection profiles from the membrane center along the x-axis are chosen for comparison. The central deflections and error analyses are listed in Table 5.
From Figure 8 and Table 5, it can be observed that the proposed model's predictions are in good agreement with the FEM simulations and the experimental data. The deviations between the proposed model and the experimental data are reducing along from center to membrane edge. After the residual stress compensation, the central deflection becomes nearer to the experimental data, and the relative errors between the proposed model and the experimental data have reduced from 5.42% to 4.14%. Bao's predictions are smaller than the experimental data, while the other four results are bigger. Rahman's deviation at membrane center is the smallest, and the deviations increase along the x-axis. The deflection profiles of the double-layer CMUT membrane for 0.1 MPa predicted by the proposed model with and without the residual stress compensation (hereinafter referred as RSC), the other literatures, and the FEM simulations are compared with the experimental data, as shown in Figure 8. The deflection profiles from the membrane center along the x-axis are chosen for comparison. The central deflections and error analyses are listed in Table 5.    The deviations between the proposed model and the experimental data may be due to the fabrication uncertainties such as manufacturing errors, residual stresses, and heavy doping. They can lower the accuracy of the analytical calculations, since the analytical calculations depend on the geometrical dimensions and material properties. For instance, the heavy doping is necessary for good conductivity, but it also changes the density of the material.

Frequency Response Analysis
The Agilent 4294A Precision Impedance Analyzer (Agilent Technologies Inc.) is used for the frequency response measurements. A 20 V DC voltage was applied and a 0.01 V AC variation was superimposed. The DC bias was chosen based on the following considerations. On one hand, a smaller DC bias voltage can avoid huge electrical coupling effects that can interfere with the validation of the proposed model. On the other hand, a larger DC bias voltage can guarantee accurate measurements and strong electrical signals. According to our calculations, the pull-in voltage of the 70 µm-side-length double-layer membrane is above 200 V, so the DC bias voltage was chosen to be less than 10% of the pull-in voltage.
The testing results and FEM simulations are illustrated in Figure 9. Since the proposed model directly gives an exact value for the response frequency, the analytical data cannot be illustrated in Figure 9. Instead, the analytical predictions are compared with other values in Table 6.
Note that the experimental data and FEM simulations are arranged using different horizontal and vertical coordinates for a clearer view. This is more suitable, since the impedance measurements represent the real part of the admittance values changing with respect to frequency, while the FEM simulations represent the vertical displacements of membrane center varying with frequency.
Observations from Figure 9 and Table 6 show that for the double-layer CMUT membrane, the proposed model's predictions strongly agree with the experimental data. With the residual stress compensation, the relative errors between the analytical prediction and the experimental data have reduced from´1.49% to´1.17%, which is quite satisfactory. The FEM simulations provide a similar response frequency and the relative error is´3.02%. Bao's prediction is bigger than the experimental data, producing a relative error of 6.76%, which also validates the multilayer membrane model. Comparisons and analyses above clearly indicate that the proposed model can offer accurate predictions for the response frequency of the double-layer CMUT membrane, and the residual stress compensation can improve the analytical predictions.
chosen to be less than 10% of the pull-in voltage.
The testing results and FEM simulations are illustrated in Figure 9. Since the proposed model directly gives an exact value for the response frequency, the analytical data cannot be illustrated in Figure 9. Instead, the analytical predictions are compared with other values in Table 6.
Note that the experimental data and FEM simulations are arranged using different horizontal and vertical coordinates for a clearer view. This is more suitable, since the impedance measurements represent the real part of the admittance values changing with respect to frequency, while the FEM simulations represent the vertical displacements of membrane center varying with frequency. Figure 9. The testing results and FEM simulations of the double-layer CMUT membrane. The black solid line represents the experimental data, arranged against the bottom x-axis and left y-axis in black. The blue solid line denotes the FEM simulations, arranged against the top x-axis and right y-axis in blue. The first peak of the experimental data is the first-order response frequency, and the second peak is the second-order response frequency.  The blue solid line denotes the FEM simulations, arranged against the top x-axis and right y-axis in blue. The first peak of the experimental data is the first-order response frequency, and the second peak is the second-order response frequency.

Conclusions and Further Study
As the crucial vibrating component, the CMUT's membrane plays an important role in CMUT performance, including the transmitting efficiency and the receiving sensitivity. Accurate determination of the CMUT's multilayer membrane is vital for design and optimization. Focused as it is on the CMUT's multilayer characterization, this paper presents a new analytical model for CMUT design and optimization. The theoretical analysis and FEM/experimental verification lead to the following conclusions: (1) A new analytical model for static deflection and the frequency response of the CMUT's multilayer membrane under pressure variation has been presented. The derivation of the proposed model and the relevant residual stress compensation method is based on the Ritz method and Hamilton's principle. The relationships between the external pressure, the static deflection, and the frequency response are obtained.
(2) For three-layer CMUT membranes, the static deflection profile and the frequency response under external pressures have been calculated by the proposed model, other literatures, and FEM simulations. The relative errors are ranging from 0.71%-3.51% for the static deflection profiles, and 0.35%-4.96% for the response frequencies, respectively.
(3) For the double-layer CMUT membrane, the static deflection profile and the frequency response under external pressure variations have been verified by the proposed model, other literatures, FEM simulations, and the experiments. The proposed model can provide accurate predictions. With the residual stress compensation, the relative errors for static deflection have reduced from 5.42% to 4.14%, and the relative errors for frequency response have reduced from´1.49% to´1.17%. This proposed analytical model can accurately and rapidly predict the static deflection and the frequency response under a wide range of external pressures of the CMUT's multilayer membrane, and can serve as a reliable reference and an accurate tool for CMUT design and optimization. At the present stage, the model only provides the mechanical characterization of CMUTs, which is basic yet not enough in practice. In our future work, equivalent circuit models will be further developed for the multilayer CMUTs. Besides, CMUTs with rectangular shape membranes will be fabricated and analyzed using our model.