Modeling of the Effect of Process Variations on a Micromachined Doubly-Clamped Beam

In the fabrication of micro-electro-mechanical systems (MEMS) devices, manufacturing process variations are usually involved. For these devices sensitive to process variations such as doubly-clamped beams, mismatches between designs and final products will exist. As a result, it underlies yield problems and will be determined by design parameter ranges and distribution functions. Topographical changes constitute process variations, such as inclination, over-etching, and undulating sidewalls in the Bosch process. In this paper, analytical models are first developed for MEMS doubly-clamped beams, concerning the mentioned geometrical variations. Then, finite-element (FE) analysis is performed to provide a guidance for model verifications. It is found that results predicted by the models agree with those of FE analysis. Assigning process variations, predictions for performance as well as yield can be made directly from the analytical models, by means of probabilistic analysis. In this paper, the footing effect is found to have a more profound effect on the resonant frequency of doubly-clamped beams during the Bosch process. As the confining process has a variation of 10.0%, the yield will have a reduction of 77.3% consequently. Under these circumstances, the prediction approaches can be utilized to guide the further MEMS device designs.


Introduction
Precise processing control has turned into an issue, owing to the mass production of micro-electro-mechanical systems (MEMS) devices and their increasingly complicated manufacturing processes. Discrepancies between initial designs and products deteriorate quickly with the feature size reductions. Even with the state-of-art fabrication techniques, process variations occur inevitably [1][2][3]. The process variations mainly include misalignment, footing as well as critical dimension (CD) loss [4], manifested as inclination, over-etching and undulating sidewalls in the Bosch process [5]. Typically, the effects of relative tolerances in MEMS devices are more severe than macro-scale products [1,6,7]. Relative manufacturing tolerances are alternatives to performance uncertainties. In addition, microstructures are commonly performed with nonlinear parallel plate electrostatic forces, which contributes to the complexity of problems. For such a reason, adequate methods are required due to the limitations in design rules and linear theories of MEMS [8].
Studies on process variations have achieved abundant results [9] in the scope of integrated circuits (ICs). However, these IC achievements cannot meet the whole needs of MEMS technologies. The conventionally used trial-and-error approach severely relies on the design-test cycle, which not only postpones the development cycle, but is also costly and time-consuming. Therefore, efforts have been made in the domain of MEMS devices analysis, for a better understanding of the impacts of process variations and also for reductions of performance variabilities at the design stage [10][11][12][13][14][15].
Islam et al. [10] conducted simulations and stress analysis on a fixed-fixed beam in electrostatic situations. Their results have reflected that changes in length and thickness tend to be more strictly controlled. Microbeam resonators are commonly utilized to detect or filter signals in MEMS. Due to manufacturing uncertainties, microbeam resonators undergo significant variability from initial designs. For example, Liu et al. [11] achieved tradeoff designs regarding multiple and conflicting design criteria, while Rong et al. [12] focused on multilayer structures while considering the first and second-order sensitivities of frequency. Mawardi et al. [13] utilized enumeration search and input-output relationships to get the governing parameters as well as a wide range for operating resonant frequency. In addition, magnetometers adopted multiphysics-based optimization and nonlinear situations [14], and gyroscopes focused on the packaging with double yield [15]. Except for unique device analysis, methods that are generally applicable have advanced the processing improvement further [3,[16][17][18][19][20][21][22][23][24][25][26]. Mirzazaden et al. [16][17][18] investigated morphology uncertainties with reduced-order models through on-chip tests. Moreover, Shavezipur et al. [3,19], and Allen et al. [20] proposed the first-order second-moment (FOSM) and advanced FOSM reliability method, respectively, in a probabilistic way to obtain a linearized feasible region and maximize the yield. For those non-linear actuated MEMS devices, high fidelity optimization schemes have also been realized. To avoid the brute-force Monte Carlo (MC) scheme, Pfingsten et al. [21] considered a Bayesian Monte Carlo approach for yield estimation, with 90.0% computational savings but the same accuracy, compared with MC schemes. Vudathu et al. [22,26] applied a sensitivity analyzer for MEMS (SAM) by worst-case analysis, revealing the effects of parametric variations on performance and yield. Achievements have been reached to get a balance between precision and calculation-for example, the Sigma-Point approach applied to MEMS resonators with four orders of magnitude faster than MC [23], the generalized polynomial chaos (GPC) framework to handle stochastic coupled electromechanical analysis with the same precision and one order of magnitude faster compared with MC [24], and the Taguchi parameter design and statistical process-control method to minimize variability in performance response to fluctuations [25].
However, little work has been carried out for detailed analysis of specific processing steps. This paper intends to explore the effect of process variations on the resonance frequency of doubly-clamped beams, under the Bosch target processing environment. The commercial code ANSYS (11.0) [27] guides the verifications of the presented methods. The results suggest that with assigned process variations, structure performance and yield can be predicted. On the other hand, given design specifications, reasonable suggestions can be made for parameter error ranges under process variations.

Process Variations
Marked as a highly anisotropic etching process with high aspect ratios, deep reactive-ion etching (DRIE) is employed to create deep penetration, steep-sided holes and trenches in wafers or substrates. The Bosch process is one of the high-rate DRIE technologies, capable of fabricating 90 • vertical walls theoretically [28][29][30]. The Bosch process alternates between isotropic plasma etching and deposition of a passivation layer, also called pulsed or time-multiplexed etching. The etching-deposition procedures will be repeated until all of the demands are satisfied. However, it is hard to obtain a sidewall precisely vertical to the substrate. Morphology features like inclinations, undulating ripples as well as over-etching are inevitable and critical, which are called the trapezium effect, the ripple effect and footing effect in the following, respectively. The variations induced by these effects manifest roughly as planar sizes (dominated by photolithography and etching processes), planar position offset (dominated by alignment) and vertical sizes (dominated by thickness variations of thin films or the substrate). An ideal beam is a cuboid structure with a length dominated as l, a width as w, and the thickness as h. The cross section was supposed to be a standard rectangle, while, in fact, it appeared as the side-view given in Figure 1a  To reveal the significance of manufacturing process variations, a simplified doubly-clamped beam is illustrated in Figure 2. The thickness and length of the beam are assumed to undergo the same manufacturing process variation as 0.05 μm. Thus, for a 200 μm long and 2 μm thick beam, the relative error for the length equals 0.05%, while it is 5.0% for the thickness case. The parameter thickness is obviously more sensitive to process variations. Combined with the usually quadruple relationship of length in a beam's frequency, the relative error diminishes to 0.0006%. This finding reveals that the length variation can be ignored in certain cases to simplify the analysis models.

Problem Solution
Studies on MEMS devices have pointed out that manufacturing process variations have a close relationship with performance drift and device failure [22,25,31]. As the basic element in MEMS, the resonant frequency of the doubly-clamped beam underlies the majority of engineering designs. Assuming the section as a plane, the doubly-clamped beam can be treated as an Euler-Bernoulli beam. Without initial buckling, the differential equation for lateral oscillation can be expressed as where ��� is the bending stiffness, ���� is the linear density, ���� is the axial load, and ( , ) is the displacement along the z-axis. Ignoring the residual stress, the resonant frequency of the doubly-clamped beam approximates as [32,33]  To reveal the significance of manufacturing process variations, a simplified doubly-clamped beam is illustrated in Figure 2. The thickness and length of the beam are assumed to undergo the same manufacturing process variation as 0.05 µm. Thus, for a 200 µm long and 2 µm thick beam, the relative error for the length equals 0.05%, while it is 5.0% for the thickness case. The parameter thickness is obviously more sensitive to process variations. Combined with the usually quadruple relationship of length in a beam's frequency, the relative error diminishes to 0.0006%. This finding reveals that the length variation can be ignored in certain cases to simplify the analysis models. To reveal the significance of manufacturing process variations, a simplified doubly-clamped beam is illustrated in Figure 2. The thickness and length of the beam are assumed to undergo the same manufacturing process variation as 0.05 μm. Thus, for a 200 μm long and 2 μm thick beam, the relative error for the length equals 0.05%, while it is 5.0% for the thickness case. The parameter thickness is obviously more sensitive to process variations. Combined with the usually quadruple relationship of length in a beam's frequency, the relative error diminishes to 0.0006%. This finding reveals that the length variation can be ignored in certain cases to simplify the analysis models.

Problem Solution
Studies on MEMS devices have pointed out that manufacturing process variations have a close relationship with performance drift and device failure [22,25,31]. As the basic element in MEMS, the resonant frequency of the doubly-clamped beam underlies the majority of engineering designs. Assuming the section as a plane, the doubly-clamped beam can be treated as an Euler-Bernoulli beam. Without initial buckling, the differential equation for lateral oscillation can be expressed as where ��� is the bending stiffness, ���� is the linear density, ���� is the axial load, and ( , ) is the displacement along the z-axis. Ignoring the residual stress, the resonant frequency of the doubly-clamped beam approximates as [32,33]

Problem Solution
Studies on MEMS devices have pointed out that manufacturing process variations have a close relationship with performance drift and device failure [22,25,31]. As the basic element in MEMS, the resonant frequency of the doubly-clamped beam underlies the majority of engineering designs. Assuming the section as a plane, the doubly-clamped beam can be treated as an Euler-Bernoulli beam. Without initial buckling, the differential equation for lateral oscillation can be expressed as where EI is the bending stiffness, ρA is the linear density, σA is the axial load, and z(x, t) is the displacement along the z-axis. Ignoring the residual stress, the resonant frequency of the doubly-clamped beam approximates as [32,33] in which k i l stands for the coefficient of the ith mode of vibration, and the first three values as k 1 l = 4.730, k 2 l = 7.853, k 3 l = 10.996.

Effect of a Single Factor
Apart from geometrical size errors that can be presented in the behavior equations (refer to Appendixs A and B), morphology changes play an important role in the variability of devices' performance and yield. Appropriate models are needed to reflect the main causes that result from manufacturing uncertainties. As stated in Section 2, process variations are mainly rooted in the trapezium, footing and ripple effect. These effects primarily occur in the Bosch process, Reactive ion etching (RIE) for silicon-on-insulator (SOI) structures, as well as time-multiple-deep deposition (TMDE), respectively. More details on the model building are shown in Appendix A.
The cross section of the beam has been transformed from an ideal rectangle to a trapezoid profile due to process variations. This phenomenon is defined as the trapezium effect, which can be divided into the positive trapezium effect (the trapezoid angle θ > 0) and the negative trapezium effect (the trapezoid angle θ < 0). Figure 3 represents the latter, where Figure 3a denotes the section of a doubly-clamped beam. Figure 3b extracts its cross section models and parameters in the coordinate system. In such a case, the resonant frequency changes into: in which stands for the coefficient of the ith mode of vibration, and the first three values as 1 = 4.730, 2 = 7.853, 3 = 10.996.

Effect of a Single Factor
Apart from geometrical size errors that can be presented in the behavior equations (refer to Appendixes A and B), morphology changes play an important role in the variability of devices' performance and yield. Appropriate models are needed to reflect the main causes that result from manufacturing uncertainties. As stated in Section 2, process variations are mainly rooted in the trapezium, footing and ripple effect. These effects primarily occur in the Bosch process, Reactive ion etching (RIE) for silicon-on-insulator (SOI) structures, as well as time-multiple-deep deposition (TMDE), respectively. More details on the model building are shown in Appendix A.
The cross section of the beam has been transformed from an ideal rectangle to a trapezoid profile due to process variations. This phenomenon is defined as the trapezium effect, which can be divided into the positive trapezium effect (the trapezoid angle > 0) and the negative trapezium effect (the trapezoid angle < 0). Figure 3 represents the latter, where Figure 3a denotes the section of a doubly-clamped beam. Figure 3b extracts its cross section models and parameters in the coordinate system. In such a case, the resonant frequency changes into: (3) The footing effect usually shows up in the RIE/DRIE process of SOI structures. It causes inhomogeneous distributions of the mass and stiffness, even making the structure part collapse. The over-etching height ℎ and horizontal over-etching width are viewed as the key elements, as demonstrated in Figure 4. Without careful processing control or wide enough width, the structure is prone to crash down, inferred from Figure 4a. The footing effect changes the resonant frequency into: The footing effect usually shows up in the RIE/DRIE process of SOI structures. It causes inhomogeneous distributions of the mass and stiffness, even making the structure part collapse. The over-etching height h f and horizontal over-etching width w f are viewed as the key elements, as demonstrated in Figure 4. Without careful processing control or wide enough width, the structure is prone to crash down, inferred from Figure 4a. The footing effect changes the resonant frequency into:  Ripples are presented on the rough side walls of structures with high aspect ratios, owing to the TMDE technology. It is defined as the ripple effect, referring to the simplified model in Figure 5a,b. Given the ripple arc ranging from 30° to 180° and single ripple height from 0.1 µm to 1 µm, simulations in ANSYS have suggested that the deciding element in the ripple effect is single ripple height, rather than ripple arc with an error within 2.0%. Therefore, single ripple height, , can be treated as a key when dealing with the ripple effect. In such a case, the resonant frequency changes into: where the meaning of the symbols is marked in Figure 5. Ripples are presented on the rough side walls of structures with high aspect ratios, owing to the TMDE technology. It is defined as the ripple effect, referring to the simplified model in Figure 5a,b. Given the ripple arc ranging from 30 • to 180 • and single ripple height from 0.1 µm to 1 µm, simulations in ANSYS have suggested that the deciding element in the ripple effect is single ripple height, rather than ripple arc with an error within 2.0%. Therefore, single ripple height, t, can be treated as a key when dealing with the ripple effect.  Ripples are presented on the rough side walls of structures with high aspect ratios, owing to the TMDE technology. It is defined as the ripple effect, referring to the simplified model in Figure 5a,b. Given the ripple arc ranging from 30° to 180° and single ripple height from 0.1 µm to 1 µm, simulations in ANSYS have suggested that the deciding element in the ripple effect is single ripple height, rather than ripple arc with an error within 2.0%. Therefore, single ripple height, , can be treated as a key when dealing with the ripple effect. In such a case, the resonant frequency changes into: where the meaning of the symbols is marked in Figure 5. In such a case, the resonant frequency changes into: where the meaning of the symbols is marked in Figure 5.

Effect of Multiple Factors
Sensitivity analyses on the effects mentioned above are conducted, shown in Figure 6. With accurate processing control and mature techniques, single arc height can be restricted to be less than 0.01 µm so that the ripple effect can be diminished, as shown in Figure 6a. Under these circumstances, models can be simplified into two critical effects: the trapezium effect and footing effect. Thus, when side walls are thought to be smooth, the resonant frequency of doubly-clamped beams can be expressed as:

Effect of Multiple Factors
Sensitivity analyses on the effects mentioned above are conducted, shown in Figure 6. With accurate processing control and mature techniques, single arc height can be restricted to be less than 0.01 µm so that the ripple effect can be diminished, as shown in Figure 6a. Under these circumstances, models can be simplified into two critical effects: the trapezium effect and footing effect. Thus, when side walls are thought to be smooth, the resonant frequency of doubly-clamped beams can be expressed as: However, side walls cannot be treated as smooth all the time. Models containing the three effects simultaneously are essential. The corresponding coordinate is illustrated as Figure 7, where 1 = + 2ℎ , 2 = 1 − 2 . When the gap between beams is wide enough, like 6 µm or wider, the cross section can be assumed to be continuous, repeated and symmetrical. Assumptions can be raised that the central axis of the cross section equals the central axis of the trapezoid ABCD, which means = , and that the ripple is a semicircle. The moment of inertia, relative to central axis = 1 , equals the condition of axis 1 , approximately. Axis 1 is in the direction of the semicircle radius and perpendicular to the side waist of trapezoid ABCD, when the semicircle vibrates longitudinally. However, side walls cannot be treated as smooth all the time. Models containing the three effects simultaneously are essential. The corresponding coordinate is illustrated as Figure 7, where b 1 = b + 2h tan θ, b 2 = b 1 − 2w f . When the gap between beams is wide enough, like 6 µm or wider, the cross section can be assumed to be continuous, repeated and symmetrical. Assumptions can be raised that the central axis of the cross section equals the central axis of the trapezoid ABCD, which means y = y c , and that the ripple is a semicircle. The moment of inertia, relative to central axis y = y 1 , equals the condition of axis z 1 , approximately. Axis z 1 is in the direction of the semicircle radius and perpendicular to the side waist of trapezoid ABCD, when the semicircle vibrates longitudinally. The above results suggest that the resonant frequency changes into (refer to Appendix B for more details on model building):

Analysis and Results
Geometric features of MEMS devices usually do not comply with the design value, with the typical error around 5.0% [34,35] during the manufacturing processes. Design parameters for the doubly-clamped beam are listed in Table 1. Multi-field models are always complicated for the complex mechanism of MEMS devices. The situation deteriorates with stochastic manufacturing uncertainties. Leaving the cost alone, with repeated adjustment or tape-out of test structures, only a small percent of the data is acceptable. Due to the lack of manufacturing data, the trial-and-error method is not optimal. Adequate models underlying process variations should be developed.
Simulations are conducted in ANSYS WORKBENCH 14.5 (ANSYS, Pittsburgh, PA, USA), with the solver SOLID45 (ANSYS, Pittsburgh, PA, USA) [27]. The mesh method uses tetrahedrons with patch conforming. Lateral vibrations are considered to perform model verifications, as shown in Figure 8. The number of ripples is supposed to be 5 and 10, in order to simplify the validations. Assuming small variation ranges in key elements of the trapezium and footing effect, analyses of curves in Figure 8 are carried out. The resonant frequency is found to be in direct proportion to the number of ripples N, with improvement in stable behavior along with larger N. Morever, the frequency manifests a reverse proportion to the angle, undergoing serious shifts in the wake of deteriorative footing effect. The above results suggest that the resonant frequency changes into (refer to Appendix B for more details on model building):

Analysis and Results
Geometric features of MEMS devices usually do not comply with the design value, with the typical error around 5.0% [34,35] during the manufacturing processes. Design parameters for the doubly-clamped beam are listed in Table 1. Multi-field models are always complicated for the complex mechanism of MEMS devices. The situation deteriorates with stochastic manufacturing uncertainties. Leaving the cost alone, with repeated adjustment or tape-out of test structures, only a small percent of the data is acceptable. Due to the lack of manufacturing data, the trial-and-error method is not optimal. Adequate models underlying process variations should be developed.
Simulations are conducted in ANSYS WORKBENCH 14.5 (ANSYS, Pittsburgh, PA, USA), with the solver SOLID45 (ANSYS, Pittsburgh, PA, USA) [27]. The mesh method uses tetrahedrons with patch conforming. Lateral vibrations are considered to perform model verifications, as shown in Figure 8. The number of ripples is supposed to be 5 and 10, in order to simplify the validations. Assuming small variation ranges in key elements of the trapezium and footing effect, analyses of curves in Figure 8 are carried out. The resonant frequency is found to be in direct proportion to the number of ripples N, with improvement in stable behavior along with larger N. Morever, the frequency manifests a reverse proportion to the angle, undergoing serious shifts in the wake of deteriorative footing effect. In light of the above development, comparisons between modified models and ANSYS are conducted, with consideration of single effects, respectively. Figure 9a describes errors under the trapezium effect, with the bias less than 2.5%. It turns out that the model and simulations share a similar trend. Figure 9b states the situation for the ripple effect. As the situation for a high value of the arc height rarely occurs during processing, the result with an error of 10.1% is not accurate or applicable for further study. The first two results are receivable in general, limiting the errors within 2.0%. Confining errors within 3%, over-etching longitudinally introduces more variabilities in resonant frequency, referred to in Figure 9c. Complicated structures will result in larger differences. Analyses are conducted in Figure 9d to explain the footing and trapezium effects. In pursuit of less variabilities in frequency, the negative trapezium effect is proved to be effective. In addition, the footing effect occupies the dominant position, compared with the trapezium effect, according to Figure 6b,c.
Furthermore, comparisons have been established between modified models and FE analysis. The number of ripples N is assigned to 10, while = ℎ = 0.2 µm in the footing effect. The curves share an error within 2.6%, according to Figure 10. The outliers in the upper right corner of Figure 10 suggest biases of ANSYS simulations. This occurrence is attributed to the unpractical assumptions that angles in the positive trapezium effect can be 20°, in which case the principles of Timoshenke beams cannot be applied directly. However, the assumptions give credit to the negative case for the existence of the footing effect. The two curves trend similarly in general, which confirms the acceptability of the modified models.
Direct Monte Carlo (MC) simulations are performed based on the modified models. Yield is defined as a factor of the proportion falling into the same distribution range. Doubly-clamped beams with 400 μm length, 10 μm width and 4 μm thickness are raised as an example in these simulations. Hypotheses are proposed that all the parameters concerned comply with the Gaussian's distribution, along with the same process variation ±0.5 μm. The sampling numbers for MC simulations range from 100,000 to 1,000,000. The relative error for frequency turns out to be around 7.9% and an angle of around ±7° when considering the trapezium effect. The resonant frequency reduces from 22.7% to 33.0% while the yield decreases to nearly 67.0% under the footing effect. When doubling the numbers of ripples, relative errors for the resonant frequency can be improved about 1.0%. In light of the above development, comparisons between modified models and ANSYS are conducted, with consideration of single effects, respectively. Figure 9a describes errors under the trapezium effect, with the bias less than 2.5%. It turns out that the model and simulations share a similar trend. Figure 9b states the situation for the ripple effect. As the situation for a high value of the arc height rarely occurs during processing, the result with an error of 10.1% is not accurate or applicable for further study. The first two results are receivable in general, limiting the errors within 2.0%. Confining errors within 3%, over-etching longitudinally introduces more variabilities in resonant frequency, referred to in Figure 9c. Complicated structures will result in larger differences. Analyses are conducted in Figure 9d to explain the footing and trapezium effects. In pursuit of less variabilities in frequency, the negative trapezium effect is proved to be effective. In addition, the footing effect occupies the dominant position, compared with the trapezium effect, according to Figure 6b,c.
Furthermore, comparisons have been established between modified models and FE analysis. The number of ripples N is assigned to 10, while w f = h f = 0.2 µm in the footing effect. The curves share an error within 2.6%, according to Figure 10. The outliers in the upper right corner of Figure 10 suggest biases of ANSYS simulations. This occurrence is attributed to the unpractical assumptions that angles in the positive trapezium effect can be 20 • , in which case the principles of Timoshenke beams cannot be applied directly. However, the assumptions give credit to the negative case for the existence of the footing effect. The two curves trend similarly in general, which confirms the acceptability of the modified models.
Direct Monte Carlo (MC) simulations are performed based on the modified models. Yield is defined as a factor of the proportion falling into the same distribution range. Doubly-clamped beams with 400 µm length, 10 µm width and 4 µm thickness are raised as an example in these simulations. Hypotheses are proposed that all the parameters concerned comply with the Gaussian's distribution, along with the same process variation ±0.5 µm. The sampling numbers for MC simulations range from 100,000 to 1,000,000. The relative error for frequency turns out to be around 7.9% and an angle of around ±7 • when considering the trapezium effect. The resonant frequency reduces from 22.7% to 33.0% while the yield decreases to nearly 67.0% under the footing effect. When doubling the numbers of ripples, relative errors for the resonant frequency can be improved about 1.0%.

Conclusions
This paper has considered the problems existing in the Bosch process and their negative influences on MEMS doubly-clamped beam performance. Modified models of doubly-clamped beams were built, with consideration for the trapezium, footing, and ripple effects respectively and simultaneously. The relative performance error was restricted to 10.0%, with a yield of about 77.3% if process variations were assumed to be in the same range. FE verifications have been performed to validate the models built in this study, indicating that the heavy simulation work can be substituted in some cases by applying the models.
The model results can be viewed as the guidance for design cycle optimizations. Designers can directly figure out the key elements in the etching process by reconsidering the design sizes and shapes, and eventually compensate the errors brought by process variations to improve the yield.
Other critical elements such as residual stress, gaps between beams, or variations in Young's modulus were not considered in this paper and will be discussed in future work. Moreover, diversified distribution forms such as quasi-Gaussian can be applied in MC methods and will be the focus of future research.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Single Factor Effect
Appendix A. 1

. The Trapezium Effect
During the DRIE processes, the ideal rectangle profile of a beam can be transferred to a trapezoid profile due to process variations. Usually, the cases are divided into the positive trapezium effect (the trapezoid angle θ > 0) and the negative trapezium effect (the trapezoid angle θ < 0), shown in Figure 3. For the ideal doubly-clamped beam, the resonance frequency can be expressed as: Ew 2 3ρl 4 .
Inducing x-y coordinates (the same way for the positive), the upper width b 1 is: The neutral axis could be denoted by y na , along with the linear equations of x 1 and x 2 , that is: Here, S z is the static moment, and A t is the current cross-sectional area. Then, the moment of inertia is expressed as: Thus, the resonant frequency of the doubly-clamped beam for the modified model considering the trapezoid effect changes into: f t = (k i l) 2 6π Finally, the frequency of the doubly-clamped beam for the modified model considering the ripple effect changes into: ht . (A14)

Appendix B. Multiple Factors Effect
In addition, it is easy to get: Models should be developed to analyze arc coordinates for the solution of I arc when the trapezoid angle θ > 0 or θ < 0. The key to the arc coordinates is the position of the first arc (see Figure A2). Therefore, Then, the frequency of the doubly-clamped beam for the modified model considering three main effects changes into In Figure A2a, z 1 and z 2 are the neutral axes of the semicircle itself, by the point of their intersection o, and y 1 is parallel to y c . Consequently, the coordinate of y 1 is y 1 = h − r cos θ − d tan θ. For the condition in Figure A2b, y 1 = h − r cos θ + (d + r tan θ) sin θ. With the assumptions above, the moment of inertia for the semicircle relative to neutral axis z 1 equals I 0 = 1 2 · 1 64 πr 4 , and then the moment of inertia relative to the cross section for the ith arc becomes Therefore,