Investigations of the Effects of Geometric Imperfections on the Nonlinear Static and Dynamic Behavior of Capacitive Micomachined Ultrasonic Transducers

In order to investigate the effects of geometric imperfections on the static and dynamic behavior of capacitive micomachined ultrasonic transducers (CMUTs), the governing equations of motion of a circular microplate with initial defection have been derived using the von Kármán plate theory while taking into account the mechanical and electrostatic nonlinearities. The partial differential equations are discretized using the differential quadrature method (DQM) and the resulting coupled nonlinear ordinary differential equations (ODEs) are solved using the harmonic balance method (HBM) coupled with the asymptotic numerical method (ANM). It is shown that the initial deflection has an impact on the static behavior of the CMUT by increasing its pull-in voltage up to 45%. Moreover, the dynamic behavior is affected by the initial deflection, enabling an increase in the resonance frequencies and the bistability domain and leading to a change of the frequency response from softening to hardening. This model allows MEMS designers to predict the nonlinear behavior of imperfect CMUT and tune its bifurcation topology in order to enhance its performances in terms of bandwidth and generated acoustic power while driving the microplate up to 80% beyond its critical amplitude.


Introduction
Microresonators have been used in several applications such as mass sensing [1], micropumps [2], gas sensors [3], gyroscopes [4] and accelerometers [5]. These microsystems are commonly fabricated with an initial deflection due to the mechanical stress caused by the fabrication process.
One of the most famous micro-systems that have been developed is the capacitive micromachined ultrasonic transducer (CMUT), which is a flexible microplate that can vibrate under an applied voltage. CMUTs have been used as transducers by emitting ultrasonic waves, also as a sensors to capture ultrasounds. Despite the excessive use of piezoelectric transducers in industrial applications, the capacitive micro-machined ultrasonic transducers (CMUTs) are more advantageous than the piezoelectric ones. First they can be fabricated using micromachining techniques [6][7][8] with tight parameter specifications and with complex shapes (circular, square, hexagonal...) [9][10][11], which makes them operate in different applications. Batch fabrication enables the micromachining of large numbers of CMUTs at the same time with good uniformity. Another important feature of CMUTs is their compatibility with integrated circuit (IC) technology. Compared to piezoelectric transducers, CMUTs have a larger bandwidth, which can result in good sensitivity (when the device is operating in the receive mode), better transmission (in the case of transmitting mode) and low noise [12].
However, the microfabrication techniques reveal several imperfections such as stress gradients [13] and initially curved shapes [14]. Even a small initial deflection, with respect to the dimensions of the resonator, has a critical effect on its static and dynamic behavior. These effects have been studied on several microresonators; for instance, the static response of a microbeam has been investigated by Ouakad et al. [15,16]. In the case of the downward initial deflection, the increase of the initial deflection leads to a decrease in the pull-in voltage due to the fact that the microbeam becomes closer to the bottom electrode. The upward initial deflection results in a snap through behavior [17,18] where the initial deflection increases the snap through voltage and decreases the pull-in voltage.
The initial imperfection also has a major effect on the resonance frequencies of a resonator. An increase of the initial deflection leads to a divergence of the odd resonance frequencies also known as frequency veering [19][20][21][22]. This phenomenon was reported in several studies such as the effect of varying the ratio of the side lengths of a rectangular plate [23] and the effect of varying the sagging levels on cables [24].
The initial deflection dose not only have an effect on the static behavior of microsystems but also affects nonlinear dynamic behavior. Ruzziconi et al. [25,26] developed a nonlinear model for a clamped microbeam obtained using Galerkin method. With only two degrees of freedom, the numerical model is able to predict with satisfactory accuracy the nonlinear dynamic behavior of the microbeam, since it correctly simulates both the primary and the 1/3 superharmonic resonances. The dynamic pull-in and dynamic snap-through of an arch microbeam was also investigated by Younis et al. [15] by exciting the microresonator with high DC and AC voltages around the fundamental resonance frequency. The increase in the shallow arch has two effects on the nonlinear behavior of the microbeam: first it reduces its vibration amplitude due to the increase of the distance between the two electrodes; and second it reduces the hardening behavior and increases its softening behavior. Experimental investigation was also conducted on the microstructure around the first and the third resonance frequency [27,28]. For low DC and AC voltages, the microsystem has a softening behavior at the first resonance frequency and a hardening behavior at the third resonance frequency. The experimental tests showed that with high voltages, we can obtain dynamic snap-through cross-well motions. Also, the microbeam can exhibit large oscillations of a continuous band of snap-through motion between the first primary resonance regime and the super-harmonic resonance regime. These characteristics of arches resonators were used in developing bandpass filters [29], low-powered switches [30] and logic devices [31].
The effect of initial imperfection was also studied for the case of microplate. Saghir et al. [32,33] derived the equation of motion of a rectangular microplate excited with an electrostatic force. The governing equations, derived using the von Kármán plate theory, are transformed into a Reduced Order Model (ROM) using the Galerkin discretization procedure. The static and the dynamic behavior of the microsystem were studied experimentally and the obtained results were compared to the numerical results obtained with the developed model. Near primary resonance, the microplate has a hardening type behavior when the excitation force is low and the increase in the AC voltage increases the gap between the two stable solutions. The dynamic behavior of the microresonator was also investigated at ω/3 super-harmonic excitation. In this case, the increase in the AC voltage does not only increase the vibration amplitude but also changes the nonlinear dynamic behavior from hardening to softening.
Several studies have been conducted on the nonlinear behavior of rectangular plate with imperfection at the macroscale. Celep et al. [34,35] investigated the effect of an initial imperfection on the static and dynamic behavior of a rectangular plate subjected to four different boundary conditions. Chena et al. [36] derived the nonlinear differential equation of the rectangular plate with initial imperfection and under arbitrary initial stress. The governing equations are discretized using the Galerkin procedure and the partial differential equations are transformed into a set of ordinary differential equations. The obtained equations were solved using the Runge-Kutta method to obtain the linear and nonlinear vibration frequencies.
Medina et al. [37,38] developed a reduced order model which has been validated using finite element analysis in order to study the axisymmetric snap-through of an initially curved circular plate subjected to an electrostatic actuation. The analysis revealed that due to the nonlinearity of the electrostatic force, the snap-through occurs at a lower displacement compared to the one induced by mechanical loads. For the case of CMUTs, few studies have been conducted on nonlinear behavior. Vogl et al. [39,40] used the von Kármán plate theory and the Galerkin method to determine the nonlinear differential equations of a circular CMUT. These equations have been solved using the method of multiple scales [41]. This analytical method is used to determine the nonlinear frequency and force responses of the CMUT around the first resonance frequency. The microsystem can display a softening or a hardening behavior depending on the applied DC voltage.
The nonlinear behavior of circular plates was also investigated at the macroscale by Touzé and Thomas [42,43]. The circular plate, with free edge boundary condition, was forced with a frequency around the asymmetric mode by taking into account the effects of damping and small geometric imperfection. The nonlinear frequency response of the system can exhibit different nonlinear behavior, softening or hardening, depending on the design parameters. From these solutions they showed that traveling waves are generic when coupled-mode solutions arise. The analytical model has been validated with respect to experimental results obtained by measuring the deflections of an antinode of each preferential configuration. By estimating the input parameter of the system, the developed model was able to predict the nonlinear phenomena of the plate.
The differential quadrature method is a numerical technique that has been used for discretization of a mathematical problem. It was derived by Richard Bellman in the early 1970s [44]. This method has been employed as an alternative technique for finite difference method and FEM. We can mention the work of Najar et al. [45] in modeling microbeams. The static and dynamic behavior of the microbeam were studied by discretizing the partial differential equation using DQM. This method is accurate and efficient. With only nine grid points, the error in the pull-in voltage is about 0.2%. The Generalized Differential Quadrature Finite Element Method (GDQFEM) was also used to solve the static and dynamic behavior of plates with complex shapes [46].
In this paper, the nonlinear behavior of initially curved circular CMUTs is modeled. The partial differential equations of the upper microplate of the CMUT are derived based on the von Kármán plate theory, while including the initial curvature, the geometric and the electrostatic nonlinearities. First, the differential quadrature method (DQM) is used to transform the governing equations into a system of ordinary differential equations with respect to time. Then, the time derivatives are dropped and the resulting coupled algebraic equations are solved in order to investigate the static behavior of the CMUT with respect to the DC voltage and the initial deflection. The obtained results are validated with respect to the experimental data. Moreover, the nonlinear frequency and force responses are investigated by solving the nonlinear differential equations using the harmonic balance method (HBM) coupled with the asymptotic numerical method (ANM). Finally, the effects of the actuation voltages and the initial deflection on the dynamic behavior of the CMUT are numerically investigated.

Mathematical Model
In this section, we consider a circular microplate with a uniform cross-section h and a homogeneous material with a density ρ and a Young's modulus E, as shown in Figure 1. The vibrating electrode is modeled as a perfectly clamped circular microplate with radius R excited with an electric voltage V(t) = V dc + V ac cos(ωt).
The bottom electrode, with radius R e , is placed below the microplate at a distance d. We suppose that the microplate is subjected to an axisymetric initial deflection w 0 (r) defined as follows: where w max represent the initial deflection of the membrane at the center. We have chosen this expression of the initial deflection because it satisfies the boundary conditions of a clamped circular microplate and also it is a good approximation to the initial deflection determined experimentally [47]. The microplate equations of motion are derived from the von Kármán plate theory. Luan and Robbins [48] proved that, at the microscale, the continuum mechanics theory is still valid. We consider the following hypothesis of the von Kármán plate theory [49]: • The material of the plate is elastic, homogeneous, and isotropic. • A straight line (filament), initially normal to the midplane, remains straight and normal to the surface during the deformation.

•
The stress normal to the midplane, σ z , is small compared to the other stress components and may be neglected in the stress-strain relations.

Problem Formulation
We follow the work of Nayfeh [50] to derive the equations of motion for an axisymmetric circular plate. We denote u(r, t) and w(r, t) the displacement component of a point in the plate with respect to the r, θ, z coordinate system. The displacement vector D is given by: where w ,r is the derivative of w with respect to r. The von Kármán nonlinear strains for a plate in large displacement are given by: where w ,rr is the second derivative of w with respect to r. ε 0 represents the initial strain due to the residual stress N 0 . e 1 and e 2 are the midplane strain components and are given as: where ν is the Poisson's ratio and w 0,r is the first derivative of the initial deflection w 0 with respect to r. In order to obtain the equations of motion and the boundary conditions, we use the extended Hamilton principle which can be written as: where δW nc , δΠ, δT and δW e denote the variations of the nonconservative, kinetic, elastic and electrostatic potential energies (W nc , Π, T and W e ), expressed, respectively, as: The variation of the electrostatic energy for a parallel plate capacitor with an initial deflection w 0 is defined as: where A is the undeformed area of the reference plane. σ and ε are the Jaumann stresses and strains, respectively. 0 is the dielectric constant of the plate material and V(t) is the applied voltage. The second time derivative and the spatial variation of the displacement vector can be written as: Substituting Equation (9) into Equation (7) and integrating over the thickness z ∈ [−h/2, h/2], yields: Substituting Equation (3) into Equation (6), yields Integrating over the thickness where M θθ and M rr are the plate moments. N θθ and N rr are the plate stress resultants: We perform partial integration for the following terms: where w 0,rr is the second derivative of the initial deflection w 0 with respect to r. Substituting (15) into (12), we obtain: Substituting (16), (10) and (8) into (5) and setting each of the coefficients of δu and δw equal to zero, we obtain where N rr,r M θθ,r and M rr,r are the first derivatives of N rr , M θθ and M rr with respect to r. M rr,rr is the second derivative of M rr with respect to r. The corresponding boundary conditions at r = 0 and r = R are defined as: rM rr = 0 or δw ,r = 0

Dimensional Equations of Motion
For thin plates, we adopt the plane stress assumption. The stress-strain relations for an linear elastic material can be stated as [49]: Therefore is the shear modulus. Substituting Equation (21) into Equations (13) and (14) and performing the integration yields By substituting Equations (22) into the Equations of motion (17) and (18), neglecting the inplane and rotatory inertia terms (since they have negligible effect on the transverse motion because the transverse natural frequencies are negligible to the inplane natural frequency [49]) and adding the damping force, the governing equations can be written as: This system of differential equations takes into consideration the effects of both the initial deflection and the residual stress. The last term between two brackets in Equation (23) represents the mid-plane stretching effects which are caused by the large displacement of the microplate. The axisymmetric boundary conditions are defined with a horizontal tangent at r = 0 and a lateral displacement equal to zero u = 0. For the clamped boundary condition (r = R), the in-plane and the out-of-plane displacements u and w are equal to zero. Therefore the boundary conditions of the clamped circular microplate are: where ∇ 4 is the bi-harmonic operator expressed as is the plate flexural rigidity and f d = cẇ is the viscous damping force. The electric forcing term in Equation (23) is a parallel-plate approximation to the capacitance where the fringing field effect is ignored for a small aspect ratio d R of the capacitor [51].

Nondimensional Equations of Motion
In order to reduce the complexity of the physical parameters, we rewrite Equations (23) and (24) in terms of nondimensional variables, which are defined as follow: Substituting these definitions into (23) and (24) and dropping the tilde (∼) in the result, we get: where Γ is defined as: wherew 0 = d w 0 and the tilde (∼) represents the dimensional quantities. For convenience, we rewrite the boundary conditions in its nondimensional form:

Differential Quadrature Method
In order to solve the equation of motion of the system, we reduce the PDEs into a set of ODEs. One of the methods that has been used in the literature is the differential quadratic method. The aim of this method is to discretize the space domain into a grid of sampling points and approximating the derivative with respect to r as a weighted sum of the function at the grid points. To this end, the derivatives can be written as: where r i are the grid points defined by Chebyshev-Gauss-Lobatto [52]: u i , w i are the in-plane and out-of-plane displacements of plate at the grid point r = r i . A

(k)
ij are the DQM weighting coefficients of the kth order derivative which are given by [53][54][55]: where A (k) ij are the respective weighting coefficients that depend only on the sampling points. The advantage of using this configuration is that it can converge only with a few number of grid points. Also, it reduces the computational time as a result of the properties of the weighting matrix A (k) [56,57].

Discretized Equations of Motion
If we apply the DQM to the equation of motion (27) and (28), we obtain: whereΓ i is the discretization ofΓ.
and the boundary conditions (30) can be rewritten as: At r = 0 : where w i 0 , w i 0,rr and w i 0,rr are respectively the initial displacement and its associated first and second derivatives with respect to r at the i th grid points.

Static Response
In order to investigate the static behavior of the microplate actuated by a DC voltage, we drop the derivative with respect to time and let w i (t) and u i (t) be time independent. Therefore the equations of motion can be written as: where w st i = w st (r i ) and u st i = u st (r i ) are the static displacements of the microplate at the i th grid point.

Implementation of the Arclength Method
In this section, we solve the system of Equations (38) and (39) using the arclength continuation method as defined in the Appendix A, where: X(s) = w st 1 (s), w st 2 (s), · · · , w st n (s), u st 1 (s), u st 2 (s), · · · , u st n (s) and α(s) = V dc (s) (40) and At this stage, it is crucial to check the convergence of the DQM solution. In Figure 2, we plot the total displacement of the microplate center (r = 0) with respect to V dc voltage for different grid numbers n. While assuming that the residual stress N 0 is negligible, the physical parameters used in the simulation are presented in Table 1. The stars marked in this graph ( * ) denote the bifurcation points where the solution changes from stable to unstable and it is called pull-in point. At these points, the slope of the curve becomes infinite. To determine the bifurcation points, we solve the following equation: At pull-in (s = s pi ), the pull-in voltage and amplitude are defined as (V dc (s pi ) , w i (s pi )). We can clearly see that the solution converges as we increase the number of grid points n. The advantage of using the DQM is that the solution converges with a few number of points, compared to the finite element method FEM [58].  In our case, the radius of the electrode is smaller than the radius of the membrane. Therefore, the electrostatic force is defined as a Heaviside function: The problem with using this type of function is that it requires an important number of points to converge because it is not a continuous function. To increase the convergence of the DQM we use two different approximations of the Heaviside function: the first one is using a smooth logistic function defined as: k is the steepness of the curve (k = 30). The second one is by approximating the Heaviside function as the derivative of a ramp function defined as: In Figure 3 we compare the Heaviside function with its approximation. The logistic function is a good approximation of the Heaviside function except near the discontinuity point (r = R e R ) where the transition is much smoother. However, for the second approximation, using the derivative of the ramp function, the solution is oscillatory especially near the singular point (r = R e R ). This effect has been reported earlier by Eftekhari [59] and Jung [60] in the case of the Dirac function. In Figure 4, we compare the pull-in voltages between the Heaviside function and its approximation for different number of grid points n. For small number of points n ≤ 12, the Heaviside function and the logistic function have the same amplitude and by increasing the number of points the logistic function converges much faster than the Heaviside function. For the second approximation, where we used the derivative of the ramp function, the solution amplitude converges to the same solution as the logistic function. Comparing the three curves, we can conclude that the approximation using the derivative of the ramp function leads to a fast convergence with high accuracy compared to the other functions. For n = 10, the static solution of the 3rd approximation converges and the error between the pull-in amplitude for n = 10 and n = 50 is less than 0.2%.  In order to validate the accuracy of the proposed model, we plot the total displacement of the plate at its center with respect to V dc and we compare it to the results presented in [47]. As it is shown in Figure 5, we obtain an excellent agreement between the experimental and the numerical results. For the numerical simulation we used DQM with only n = 10 points since 10 points are sufficient to ensure the convergence of the numerical solution. Actually, the static displacement increases monotonously and it is characterized by a vertical slope at the pull-in voltage equal to 81 V and for which the microplate touches the bottom electrode and collapses. Around this excitation level, the plate displacement becomes very sensitive to voltage variations. Clearly, the microplate has no snap-through motion since its static profile shows only one stable mode of operation [18,61,62]. This is because the initial deflection is small compared to the radial dimension of the microplate (w max << R) and the residual stress is negligible. For the lower branch solution, the two analytical models based on DQM and Galerkin method, converge with a few degrees of freedom. However, the use of an even number of modes in the ROM leads to the divergence of the upper branch solution [45]. Then, we investigate theoretically the effect of the initial deflection on the pull-in voltage of circular microplates. Figure 6a displays the variation of the static displacement with respect to the applied DC load for microplates with different initial deflections. The stars ( * ) marked in this graph denote the bifurcation points where the solution becomes unstable and goes to pull-in, which increases with respect to the initial deflection. The pull-in voltage increases up to 45% and this is due to the fact that the effective gap between the two electrodes is larger than the expected one d = 750 nm, so that a larger electrostatic force is required to reach the pull-in. In order to determine the effective gap distance, we identify the equivalent electrostatic force in the case of initially curved microplates, by integrating over z ∈ [0, R e ]. We obtain the following expression [47]: The square of the effective gap distance is approximated as the mean of the square of the total gap distance. In Figure 6b, we compare the pull-in voltage of the circular CMUT calculated using the presented model and approximated using the effective gap distance d e f f . The proposed approximation has a good agreement with the numerical model especially at small deflections. In our case, the maximum error between the two curves in the interval [−300 nm, 300 nm] is about 5%. Therefore, the effective gap distance is a good approximation to determine the pull-in voltage of a CMUT with initial deflection. However, this approximation is only valid when the initial deflection is small (less than 40% of the gap distance).

Eigenfrequency Analysis
The initial deflection does not only affect the static response of the microplate, but also it modifies its resonance frequencies. The linear undamped eigenvalue problem of a clamped circular plate excited with a DC voltage is obtained by linearizing (34) and (35) where {φ i , ψ i } are the mode shapes of w and u at r = r i . Substituting (47) into Equations (34) and (35), we obtain:ẅ The electrostatic force for a circular microplate with initial deflection w max can be expanded into Taylor series around the static solution: For n grid points, the total number of equations is equal to 2n. At the boundary nodes, we use the boundary Equation (37) to determine the boundary displacement {w 1 , w n−1 , w n , u 1 , u n−1 } as function of w i and u j where i = 2, · · · , (n − 2) and j = 2, · · · , (n − 1). This simplification reduces the total number of equations to 2n − 5.
The resonance frequencies of the structure and the mode shapes are determined by solving the eigenvalue problem of the system of Equation (48), which can be written as: where M, C and K are the mass, damping and stiffness matrices; K NL (Y(t),Ẏ(t)) is the vector with nonlinear parameters and Y(t) is the vector of unknowns defined as: By dropping the nonlinear and damping terms and by replacing w dy i (t) = φ dy i e jwt and u dy i (t) = ψ dy i e jwt , the linear natural frequencies are determined using the following characteristic equation.
As we can notice in Equation (48), the nondimensional eigenfrequencies depend mainly on the axial force N 0 , the DC voltage V dc , the static displacement w st i and the initial deflection w 0 . As shown in Figure 7, the resonance frequencies decrease as we increase the DC voltage. At the pull in voltage the first resonance frequency drops to zero due to the increase of the electrostatic force. Also, we show the influence of the initial deflection on the first two resonance frequencies.
As displayed in Figure 6b, the pull in voltage increases as we increase the initial deflection. At a constant DC voltage, the initial deflection shifts the resonance frequency of the microplate. This is due to the increase of the gap distance which reduces the electrostatic force, i.e., an important actuation voltage is needed to bend the microplate.

Nonlinear Dynamic Analysis
For n grid points, the total number of differential equations to solve is equal to 2n − 5, which is increases the computational time. To reduce the number of equations, we write the radial displacement u i as a function of w i using Equation (35). Therefore, Equation (34) can be transformed into: (54) where Equation (54) is solved using DQM coupled with the HBM and the ANM, as it is defined in [63,64]. Doing so, one can plot the nonlinear frequency and force responses of the CMUT under primary resonance.

Frequency Response
The frequency response of the microplate is obtained by solving Equation (54) at different excitation frequencies. In order to examine the convergence of the DQM with respect to the number of grid points we simulate the frequency response of a CMUT actuated with a DC voltage V dc = 50 V and an AC voltage V ac = 1 V for a quality factor Q = 100. In Figure 8, we present the effect of the number of the grid points on the frequency response of the CMUT around its first resonance frequency. The use of a low number of points (n = 6) leads to a shift in the resonance frequency and an error in the response of the microplate. By increasing n, the nonlinear response of the CMUT converges, and the use of 10 grid points is sufficient to describe the nonlinear dynamic response of the device. The maximum error between the two curves for n = 10 and n = 17 is less than 1%.  Figure 8. Convergence of the dynamic solution of the microplate actuated with V dc = 50 V and V ac = 1 V in the case of a quality factor Q = 100.
In order to investigate the effects of the actuation voltages, we plot in Figure 9 the frequency response of the CMUT actuated with V ac = 1 V and for several DC voltages. For a low DC voltage (V dc = 20 V), the vibration amplitude of the microplate is low and the system exhibits a linear behavior. By increasing the DC voltage, the out-of plane displacement w(r, t) increases and the curve bends to the left and displays a softening behavior due to the increase in the electrostatic nonlinearities. Moreover, the resonance frequency decreases due to the negative stiffness term which is proportional to V 2 dc and inversely proportional to the cube of the gap distance as it is shown in Equation (50).  Figure 9. Effects of changing the DC voltage on the frequency responses of the microplate for a V ac = 1 V and with a quality factor Q = 100.
Unlike the effect of V dc on the CMUT frequency response, Figure 10 shows that the variation of the AC voltage has a significant impact on the excitation amplitude and a negligible one on the negative electrostatic stiffness. Actually, increasing V ac leads to an increase in the vibration amplitude and enlarges the bistability domain without any remarkable frequency shift. This phenomenon is displayed in Figure 10a,b for the case of low and high DC voltage.
In Figure 10a, the CMUT is actuated with a DC voltage V dc = 20 V. In this case, the increase of the AC voltage from V ac = 3 V up to V ac = 4 V leads to an increase in the hardening effect on the nonlinear behavior of the CMUT, while vibrating slightly below the critical amplitude [65][66][67][68]. This is due to the fact that the geometric nonlinearities dominate the microplate dynamics. By increasing more the AC voltage up to (V ac = 5 V), the CMUT exhibits a mixed behavior [69] due to the increase of the electrostatic force [70], which increases the CMUT vibration up to 80% beyond the critical amplitude.
Remarkably, for V dc = 50 V (Figure 10b), the increase of the AC voltage leads to a softening behavior beyond the critical amplitude, due to the dominance of the electrostatic nonlinearities over the geometric nonlinearities.
The proposed model enables the capture of the main dynamical phenomena of the CMUT in the frequency domain, which describes the competition between mechanical and electrostatic nonlinearities.
In Figure 11, we present the effect of the initial deflection on the frequency response of a microplate actuated with the same AC and DC voltage (V dc = 20 V and V ac = 4 V. While a flat microplate has a softening behavior, an increase in the initial deflection leads to a change in the CMUT behavior from softening to hardening. The initial deflection increases the gap distance between the two electrodes, which decreases the electrostatic force. In this case, the geometric nonlinearity becomes larger compared to the electrostatic nonlinearity. Consequently, the frequency response curve bends to the right and exhibits a hardening behavior.     Figure 11. The effect of the initial deflection on the frequency response for the case of microplate excited with V dc = 15 V and V ac = 0.25 V.

Force Response Analysis
For a better understanding of the effect of the excitation amplitude, we display in Figure 12 the force response curves of a circular CMUT with different initial deflections. The force response curves are constructed by varying the excitation voltage V ac while the excitation frequency ω and the bias voltage V dc are set to constant values. For a null AC voltage, the vibration amplitude is equal to zero, because in this case the CMUT is excited only with a DC voltage. For a low V ac , the solution amplitude has one possible solution. By raising the AC voltage, the vibration amplitude increases and the dynamic solution becomes nonlinear and displays bistability.
In Figure 12, we display the effect of an initial deflection on the force response curve in the case of a softening behavior. It is shown that as the initial deflection increases, the bifurcation point is shifted to the right, which enlarges the bistability domain. This is due to the increase in the effective gap distance, which means that a higher AC voltage is needed to reach the bifurcation point.

Conclusions
In this paper, the mathematical model of a circular CMUT has been derived while taking into account the effects of an axisymmetric initial deflection and the main sources of nonlinearities. The partial differential equations have been spatially discretized using the DQM. The resulting equations are then transformed into coupled ODEs. For the CMUT static behavior, we presented the good agreement between the experimental results and the numerical solutions, which proves that DQM, with only 10 grid points, is an efficient tool to study the static response of such a device. We showed that a positive initial imperfection increases the pull-in voltage of the CMUT. Also, we presented an approximation of the effective gap distance that can be used with a flat plate model to predict the pull-in voltage.
The numerical model has been also used to investigate the effect of the initial deflection on the resonance frequencies of a CMUT at several DC voltages. For the nonlinear frequency responses, the ODEs have been solved using the HBM coupled with the ANM. It has been shown that the CMUT can exhibit a hardening, a softening or a mixed behavior depending on the forcing parameter V dc . Also, we have highlighted that the initial deflection can change the frequency response of the CMUT from softening to hardening, by increasing the effective gap distance, which decreases the electrostatic nonlinearities. Finally, we have shown that the initial deflection shifts the bifurcation points to the right while enlarging the bistability domain.
Future work will include the investigation of the effects of the squeezed air film and the fluidic medium on the nonlinear response of the microplate.