Shape optimization of geometrically nonlinear modal coupling coefficients: an application to MEMS gyroscopes | Scientific Reports

Apr 04, 2025

Scientific Reports volume 15, Article number: 10957 (2025) Cite this article

394 Accesses

Metrics details

Micro- and nanoelectromechanical system (MEMS and NEMS) resonators can exhibit rich nonlinear dynamics as they are often operated at large amplitudes with high quality factors and possess a high mode density with a variety of nonlinear modal couplings. Their impact is strongly influenced by internal resonance conditions and by the strength of the modal coupling coefficients. On one hand, strong nonlinear couplings are of academic interest and promise novel device concepts. On the other hand, however, they have the potential to disturb the linear system behavior on which industrial devices such as gyroscopes and micro mirrors are based. In either case, being able to optimize the coupling coefficients by design is certainly beneficial. A main source of nonlinear modal couplings are geometric nonlinearities. In this work, we apply node-based shape optimization to tune the geometrically nonlinear 3-wave coupling coefficients of a MEMS gyroscope. We demonstrate that individual coupling coefficients can be tuned over several orders of magnitude by shape optimization, while satisfying typical constraints on manufacturability and operability of the devices. The optimized designs contain unintuitive geometrical features far away from any solution an experienced human MEMS or NEMS designer could have thought of. Thus, this work demonstrates the power of shape optimization for tailoring the complex nonlinear dynamic properties of MEMS and NEMS resonators.

Nonlinear dynamic phenomena are often observed in micro- and nanoelectromechanical system (MEMS and NEMS) resonators. They include simple Duffing-like behavior, nonlinear damping, parametric excitations, or multidmodal effects, just to name a few1,2,3,4. Nonlinear dynamics related to internal resonances between different modes5,6,7 are particularly relevant in MEMS and NEMS resonators, as shown, e.g., in8. Such effects depend on the strength of nonlinear coupling coefficients and resonance conditions between the modes.

It has been proposed to exploit nonlinearities for novel device concepts. For example, nonlinear modal couplings can be employed for angular rate detection in gyroscopes, enhancement of a gyroscope’s sensitivity, frequency stabilization in resonators, energy harvesting, or generation of phononic frequency combs9,10,11,12,13. In these cases, being able to increase the modal coupling coefficients would enhance favorable effects that arise due to the related nonlinearities.

On the other hand, in typical commercial applications, e.g., gyroscopes or micro mirrors, nonlinear effects related to 3-wave and 4-wave modal couplings can interfere with device functionality and even lead to failure14,15. The governing internal resonance conditions are highly sensitive to fabrication tolerances, whereas the involved nonlinear modal coupling coefficients are not16. Being able to reduce the modal coupling coefficients sufficiently by design could ensure device functionality despite internal resonances.

In MEMS and NEMS resonators, oscillation amplitudes are often large compared to their structural dimensions. In such a large mechanical deflection scenario, the quadratic displacement components of the mechanical strain become relevant and induce geometric nonlinearities which depend on the mechanical structure and its mode shapes17. Thus, geometric nonlinearities are one of the main sources of the aforementioned nonlinear phenomena, and structural optimization can be employed to tailor them. So far, the majority of publications on the optimization of MEMS and NEMS resonators was concerned with eigenfrequencies and quality factors18,19,20,21,22. Regarding nonlinear modal couplings, it has been shown that topology optimization can be used to avoid internal resonances in simple structures23. Recently, we developed a node-based shape optimization methodology and applied it to a MEMS gyroscope in order to avoid internal resonances24. Alternatively, shape optimization has been applied to tune the geometrically nonlinear Duffing coefficient and a 3-wave coupling coefficient in MEMS resonators25. Based on this approach, designs with optimized Duffing coefficients were fabricated and experimentally verified26. The approach in25 is limited to structures that can be approximated by beam elements and optimizes their width. For example, this implies that an initially straight beam can not become curved. Furthermore, it has been shown that one can also tune nonlinear coefficients of MEMS resonators via topological changes27.

In this work, we extend our previously developed node-based shape optimization methodology24 to the optimization of geometrically nonlinear 3-wave modal coupling coefficients and exemplify our approach on an industrial MEMS gyroscope. In contrast to our previous work, where we optimized eigenfrequencies, this requires optimization of functions that depend on mode shapes. The advantage of our node-based shape optimization is that it can generate significant shape changes, which go beyond variations in widths. Thus, the presented approach is applicable to arbitrarily complex geometries.

In this section, we introduce the underlying equations, which govern the dynamics of a mechanical resonator. In particular, we will focus on geometric nonlinearities and the resulting coupling of the structure’s eigenmodes.

It is useful to formulate the resonator’s dynamics in terms of its linear eigenmodes, since they form a basis in which the displacement can be expressed. For that purpose, we employ the finite element method (FEM) to simulate arbitrarily complex geometries. The structure is discretized and the mechanical displacement and geometry are interpolated by shape functions. This leads to the generalized eigenvalue problem

with linear stiffness matrix \(\varvec{K}\), mass matrix \(\varvec{M}\), eigenvector \(\varvec{\phi }_i\), angular eigenfrequency \(\omega _i=2\pi f_i\) and eigenfrequency \(f_i\) of mode i. A modal analysis is performed to obtain the relevant eigenmodes from Eq. (1). The eigenvectors are then mass-normalized, i.e., \(\varvec{\phi }_i^T\varvec{M}\varvec{\phi }_i=1\).

MEMS and NEMS resonators are typically driven to oscillate at relatively large amplitudes such that geometrically nonlinear effects are not negligible. These effects are solely determined by the geometry of the resonator and its material parameters. We assume linear material behavior, which is customary for standard MEMS and NEMS materials, such as, e.g., silicon. Starting from continuum mechanics, the strain energy U of an elastic body is given by

with the second Piola–Kirchhoff stress tensor \(\varvec{S}\), the Green-Lagrange strain tensor \(\varvec{E}\) and the undeformed volume of the body V28. Note that Eq. (2) contains terms that are quadratic, cubic and quartic in the displacement \(\varvec{u}\). The displacement can be written in the basis of the linear eigenmodes as

where N is the number of degrees of freedom in the discretized model and \(q_i\) is the modal amplitude of mode i. Inserting the modal superposition, Eq. (3), into the strain energy, Eq. (2), one obtains

where \(\alpha _{n,m,l}\) and \(\beta _{n,m,l,k}\) are modal 3-wave and 4-wave coupling coefficients that arise due to geometric nonlinearities. This work focuses on the 3-wave coupling coefficients \(\alpha _{n,m,l}\). Within the FEM approximation they can be calculated as

where e is the element index and the summation is performed over all elements, \(V^e\) is the undeformed volume of element e, \(\varvec{\varepsilon }^e_n\) is the linear strain of eigenvector n within element e, \(\varvec{D}\) is the elasticity matrix and \(\varvec{\eta }_{m,l}^e\) is the nonlinear strain of modes m and l within element e. Note that Voigt notation is used here by collecting the components of the symmetric strain tensors in vectors \(\varvec{\varepsilon }^e_n\) and \(\varvec{\eta }_{m,l}^e\). The linear and nonlinear strains can be obtained as

where \(\varvec{B}^e_{\varepsilon }\) contains shape function derivatives and is identical to the strain-displacement matrix of the linear theory and \(\varvec{B}^e_{\eta }\left( \varvec{\phi }^e_m\right)\) contains products of shape function derivatives and derivatives of \(\varvec{\phi }^e_m\), which is the eigenvector of mode m inside element e. The structures of \(\varvec{B}^e_{\varepsilon }\) and \(\varvec{B}^e_{\eta }\) can be found in literature28.

An equation of motion can be derived from Eq. (4) for each mode, which determines its time-dependent modal amplitude and reads

where \(Q_n\) is the quality factor of mode n, \(F_n\) is the modal force acting on mode n and we defined

Equation (8) contains quadratic and cubic nonlinearities, which couple mode n to all other modes. The significance of these couplings is determined by resonance conditions between the eigenfrequencies of the modes and the magnitude of the corresponding coupling coefficients. A detailed description of many of the nonlinear effects that can arise from Eq. (8) can be found, e.g., in5. To facilitate the numerical modelling of the underlying systems with potentially large FEM models, efficient approaches have been developed to simulate the dynamics arising due to Eq. (8) in a reduced order model (ROM)29.

In this section, we give an overview of the employed shape optimization methodology. Particular attention is being paid to the sensitivity analysis of eigenvector-dependent functions, such as the 3-wave coupling coefficients. For that purpose, the adjoint method is introduced.

We build upon the shape optimization methodology that we introduced previously24. A detailed explanation of the approach can be found therein and we will only give a brief summary here.

We apply node-based shape optimization to geometries that are extruded along the z-direction. For that purpose, parameters \(p_j\) are introduced, which shift the exterior nodes of the geometry along the surface normals in the xy-plane. The design parameters are collected in a vector \(\varvec{p}\). In order to avoid sharp kinks in the optimized design, the parameters are defined such that they also drag neighboring nodes along in a linearly decaying manner. Essentially, each parameter \(p_j\) corresponds to one exterior node of the geometry. Increasing the value of \(p_j\) will shift the corresponding node, as well as some neighboring nodes, along the geometry’s surface normal in the xy-plane. By that, the shape of the extruded geometry can be varied in a parametrized way and gradients with respect to such shape changes can be calculated analytically.

We are concerned with optimization problems where the objective function and the constraints depend on parameters explicitly as well as implicitly through eigenfrequencies and eigenvectors. These optimization problems are generally nonlinear and require an iterative solution procedure. The aim is to find the parameter set \(\varvec{p}\), which modify the shape of the geometry such that all specified constraints are fulfilled. Therefore, our optimization loop starts by performing a modal analysis of the current design. Subsequently, the objective function and constraints are evaluated and their sensitivities with respect to the design parameters are calculated. The method of moving asymptotes (MMA)30 is then employed to solve the nonlinear optimization problem and provide an updated set of design parameters, based on the function values and sensitivities. The design parameters are then applied to shift the positions of the node coordinates. By that, an updated design is obtained. The steps are repeated until all constraints are satisfied.

We wish to optimize functions

which have explicit and implicit dependencies on a given parameter \(p_j\). In general, a function c can depend on any number of parameters \(p_{j}\), eigenfrequencies \(f_{i}\) and eigenvectors \(\varvec{\phi }_i\). Evaluating the sensitivity \(\frac{\textrm{d}c}{\textrm{d}p_{j}}\) by calculating the eigenvector sensitivity \(\frac{\textrm{d}\varvec{\phi }_i}{\textrm{d}p_{j}}\), for potentially thousands of parameters, is inefficient. The adjoint method can be applied to evaluate \(\frac{\textrm{d}c}{\textrm{d}p_{j}}\), without having to obtain the eigenvector sensitivity, leading to a reduced computational complexity31. For that purpose, the Lagrangian \({\mathcal {L}}\) is defined by adding zeros to Eq. (11) as

where \(\varvec{\lambda }_{i}\) and \(\eta _i\) are the adjoint variables and the index i runs over all modes on which c depends. The adjoint variables can be chosen arbitrarily, since they are being multiplied by terms which equal zero. Therefore, Eqs. (11) and (12) have identical values. However, calculating the sensitivity of Eq. (12) is more efficient and can be obtained as

where the adjoint variable \(\eta _i\) has to be calculated according to

and the adjoint variable \(\varvec{\lambda }_{i}\) has to be obtained by solving the equation system

Equations (13)–(16) were obtained by taking the derivative of Eq. (12) and then choosing the adjoint variables such that all terms proportional to \(\frac{\textrm{d}\varvec{\phi }_i}{\textrm{d}p_{j}}\) and \(\frac{\textrm{d}f_i}{\textrm{d}p_{j}}\) cancel out. We solve Eq. (15) and Eq. (16) with Nelson’s method32,33. One can also solve an equation system for \(\eta _i\) and \(\varvec{\lambda }_{i}\) simultaneously25,31, but this would lead to a coeffcient matrix with an increased bandwidth and is less efficient. For all functions that depend on eigenvectors we will use the adjoint method according to Eqs. (13)–(16) for sensitivity analysis. The required partial derivatives \(\frac{\partial c}{\partial p_j}\), \(\frac{\partial c}{\partial \varvec{\phi }_i}\) and \(\frac{\partial c}{\partial f_i}\) are derived in Appendix A.1 for the 3-wave coupling coefficients. Details regarding the calculation of the stiffness and mass matrix sensitivities, which appear in Eq. (13), can be found in24.

In this section, we briefly describe the MEMS gyroscope, on which we will demonstrate our methodology. Furthermore, we formulate two optimization problems, to demonstrate the capability of tuning geometrically nonlinear 3-wave coupling coefficients by shape optimization. We choose an exemplary MEMS gyroscope, to highlight that, in contrast to previous work25, our method requires no simplifications of the model, such as using beam-theory, and is applicable to MEMS structures with a realistic complexity. The application to MEMS gyroscopes holds also practical relevance, as the occurence of nonlinear modal couplings can be detrimental to their operation14,15.

Top view of the single-axis MEMS gyroscope’s initial design. The dashed red lines show the symmetry axes, blue elements indicate the substrate anchors, green elements the springs and orange elements the comb electrodes.

Mode shapes of the initial design. Only the modes relevant to the optimizations are displayed. For simplicity, only the displacement of the surface is shown. The coloring indicates the sum of the three displacement components at each point, normalized to each mode’s maximum value. (a) Mode 1 with an eigenfrequency of 22.0 kHz. (b) Mode 3 (drive) with an eigenfrequency of 24.4 kHz. (c) Mode 4 (sense) with an eigenfrequency of 27.3 kHz. (d) Mode 7 with an eigenfrequency of 51.0 kHz. (e) Mode 9 with an eigenfrequency of 55.7 kHz. (f) Mode 13 with an eigenfrequency of 70.3 kHz.

The single-axis MEMS gyroscope, that we will shape-optimize here, has already been subject of previous publications, which dealt with the modeling of nonlinearities34,35,36,37,38 and the shape optimization of eigenfrequencies24, where a more detailed description of the model can be found. Figure 1 depicts a top view of the initial design, as the three-dimensional (3D) model is simply extruded along the out-of-plane direction. We assume that the entire mechanical structure is composed of linear isotropic polycrystalline silicon. The design exhibits a quarter symmetry in the xy-plane. In Fig. 1, the substrate anchors are colored in blue, the springs in green and the comb drive electrodes in orange. Several eigenmodes of the initial design, which will be relevant in this work, are shown in Fig. 2. The operation of the MEMS gyroscope is based on the drive mode, shown in Fig. 2b, and the sense mode, shown in Fig. 2c. The relevance of the remaining modes in Fig. 2 will be explained later on. An alternating voltage applied to the comb electrodes forces the drive mode to oscillate. An angular rate around the x-axis then excites the sense mode via the Coriolis force. The oscillation of the sense mode can be measured capacitively via electrodes underneath the masses.

We will demonstrate the shape optimization of geometrically nonlinear 3-wave coupling coefficients with two different optimization problems. Both optimizations will be performed on the MEMS gyroscope in Fig. 1. We apply the node-based shape optimization to the springs of the gyroscope. The design parameters are defined such that they shift the exterior nodes of the springs along the surface normals in the xy-plane. We formulate the optimization problems such that the squared norm of the design parameters \(\varvec{p}\cdot \varvec{p}\) is minimized. By that, the optimization converges towards the design which is most similar to the initial design. The actual goals of the optimization will be introduced as constraints. In both optimization problems, we will enforce that the eigenfrequencies of drive mode \(f_d\) and sense mode \(f_s\) remain fixed within a tolerance. Additionally, we will constrain the mode shape of the drive mode \(\varvec{\phi }_{d}\) to remain parallel to the y-axis. This is motivated by the fact that the initial layout of a gyroscope is usually already designed such that drive and sense mode have desired eigenfrequencies and mode shapes, i.e., basic sensor functionality is ensured. Furthermore, manufacturability constraints on the minimum width of structures and minimum distance between structures will be employed. These are customary for MEMS resonators, since typical etching processes can not release structures which are arbitrarily close to each other or too wide.

The aim of the first optimization problem is the reduction of two different geometrically nonlinear 3-wave couplings between drive mode and two spurious modes by factors of 1000. We found that 3-wave couplings between in-plane, out-of-plane and torsional modes are typically large in our design. Therefore, we will decrease the 3-wave coupling coefficients \({\bar{\alpha }}_{d,1,9}={\tilde{\alpha }}_{d,1,9}+{\tilde{\alpha }}_{d,9,1}\) and \({\bar{\alpha }}_{d,7,9}={\tilde{\alpha }}_{d,7,9}+{\tilde{\alpha }}_{d,9,7}\) between the drive mode \(d=3\), a torsional mode 9, and two distinct out-of-plane modes 1 and 7, see Fig. 2a,b,d and e. Note that \({\bar{\alpha }}_{d,1,9}\) and \({\bar{\alpha }}_{d,7,9}\) contain all nonlinear 3-wave couplings proportional to \(q_{d}q_1q_9\) and \(q_{d}q_7q_9\) in the potential energy, Eq. (4). We optimize their absolute values, as these determine the magnitude of the nonlinear coupling effects. The first optimization problem is formulated as

where the subscript 0 indicates quantities in the initial design, \(\varvec{\phi }_{d,x}\) and \(\varvec{\phi }_{d,y}\) refer to the x- and y-components of the drive mode’s eigenvector, \(V_{el}\) is the domain of the comb drive electrodes which is highlighted in orange in Fig. 1, \(d_j\) is the distance between structures, \(w_j\) is the width of structures and \(n_p=3413\) is the number of design parameters. The initial values of the 3-wave coupling coefficients were \(|{\bar{\alpha }}_{d,1,9}|_0=4.1 \times 10^{18}{\frac{1}{\sqrt{{\textrm{kg}}}\cdot \hbox {m}\cdot \hbox {s}^2}}\) and \(|{\bar{\alpha }}_{d,7,9}|_0=5.0 \times 10^{18}{\frac{1}{\sqrt{{\textrm{kg}}}\cdot \hbox {m}\cdot \hbox {s}^2}}\). The minimization in Eq. (17) ensures that the design is close to the initial design; the constraints on \(|{\bar{\alpha }}_{d,1,9}|\) and \(|{\bar{\alpha }}_{d,7,9}|\) reduce the magnitude of the 3-wave couplings by a factor of 1000 with respect to their initial values; the constraints on \(f_{d}\) and \(f_{s}\) keep the drive and sense mode eigenfrequencies fixed within \(\pm 1\%\) of their initial values; the constraint on the ratio between \(\varvec{\phi }_{d,x}\) and \(\varvec{\phi }_{d,y}\) ensures that the drive mode shape remains parallel to the y-axis within the comb drive electrodes, so that the sensor can be properly actuated; the constraints on \(d_j\) and \(w_j\) enforce that structures are at least \(2\,\upmu \hbox {m}\) apart and \(1.5\,\upmu \hbox {m}\) wide, to ensure the manufacturability with customary MEMS fabrication processes. The limits on \(p_j\) were chosen such that a converged solution could be obtained without being unnecessarily large. Details on the calculation of \(d_j\) and \(w_j\) can be found in24.

In contrast to the first optimization, the second optimization problem will demonstrate the enhancement of geometrically nonlinear 3-wave effects. We will enforce a 1:2 internal resonance and increase the corresponding coupling coefficient by a factor of 250. The factor of 250 is due to the fact that this was roughly the largest coupling increase that we could obtain with our methodology. We consider the 1:2 coupling between drive mode and mode 7. Their mode shapes are shown in Fig. 2b and d. The second optimization problem reads

where the constraint on \(|{\tilde{\alpha }}_{d,d,7}|\) leads to an increase of the 1:2 internal resonance coupling coefficient between drive mode and mode 7 by a factor of 250 with respect to the initial value \(|{\tilde{\alpha }}_{d,d,7}|_0=2.2 \times 10^{15}{\frac{1}{\sqrt{{\textrm{kg}}}\cdot \hbox {m}\cdot \hbox {s}^2}}\); the constraint on \(|2f_{d}-f_7|\) enforces that the frequencies of drive mode and mode 7 deviate no more than 100 Hz from the ideal 1:2 internal resonance condition; the constraint on \(|f_7-f_{13}|\) ensures that the eigenfrequencies of mode 7 and 13 remain at least 2 kHz apart to avoid mode mixing and the remaining constraints and minimization are as in Eq. (17). Without the constraint on \(|f_7-f_{13}|\), we found that modes 7 and 13 would mix. Looking at Fig. 2d and f, this meant that mode 7, which is initially an out-of-plane mode, also exhibited significant in-plane motion within the springs, similar to mode 13, leading to a large coupling \(|{\tilde{\alpha }}_{d,d,7}|\). However, this would be an undesirable solution, as such an increase, based on mode mixing, would likely be very sensitive to the specific values of \(f_7\) and \(f_{13}\) and therefore to fabrication tolerances. Note that \({\tilde{\alpha }}_{d,d,7}\) contains all coupling coefficients that appear in the potential energy, Eq. (4), proportional to \(q_{d}^2q_7\).

The sensitivities of the objective function and constraints in Eqs. (17) and (18), with respect to the design parameters, are required for the gradient-based shape optimization. The sensitivity of the objective function is trivially obtained. For all constraints that depend on eigenfrequencies but not on eigenvectors we calculate the sensitivities of the eigenfrequencies as in24 and obtain the sensitivities of the constraints via the chain rule. The sensitivities of the manufacturability constraints are also calculated as in24. For all constraints that depend on eigenvectors we apply the adjoint method, according to Eqs. (13)–(16). The partial derivatives, required for the adjoint method, are detailed in Appendix A. During the shape optimization, the initial order of the modes might change. In order to consistently refer each constraint to the correct mode, we track the modes as described in Appendix B. We exploit the quarter symmetry of the initial design during the optimization. By that, the optimized designs also exhibit a quarter symmetry. All simulations and calculations were performed in a self-written FEM code in MATLAB.

History of the optimized 3-wave coupling coefficients during the optimization process. The dashed lines indicate the target values. (a) First optimization problem, Eq. (17). (b) Second optimization problem, Eq. (18).

Final quarter model of the converged shape optimization according to Eq. (17). The full design has a quarter symmetry, which is exploited in the optimization procedure. The springs are colored in green.

Final quarter model of the converged shape optimization according to Eq. (18). The full design has a quarter symmetry, which is exploited in the optimization procedure. The springs are colored in green.

Dependence of the 3-wave coupling coefficients on process variations in the initial and optimized designs. Their values are shown for the nominal design, on which the optimization was based, as well as for the same design which was either dilated or eroded. The dilated design emulates an under-etching of the entire model such that all structures are thicker by 100 nm in the xy-plane. The eroded design emulates an over-etching of the entire model such that all structures are thinner by 100 nm in the xy-plane.

The convergence history of the optimized 3-wave couplings for both optimization problems is shown in Fig. 3. The optimized quarter designs, corresponding to one quadrant of the complete structure depicted in Fig. 1, are shown in Figs. 4 and 5 for the first and second optimization problem, respectively. The mode shapes of the optimized designs, corresponding to the initial modes in Fig. 2, are shown in Appendix C, for completeness. However, they do not deviate significantly from the initial mode shapes.

The first optimization problem, as stated in Eq. (17), converged after 34 iterations. It can be seen in Fig. 3a that the two depicted 3-wave couplings decreased continuously, except for a few iterations. Decreasing them by factors of 1000 took 34 iterations and the remaining constraints in Eq. (17) were already fulfilled in fewer iterations. We assume that the spikes in Fig. 3a are due to the strong nonlinear dependence of the 3-wave couplings on the design parameters. A line search algorithm could be used to determine a step size which ensures reduction of the 3-wave coupling coefficients in every iteration. Looking at the optimized quarter layout in Fig. 4, we see rather subtle design changes. Nevertheless, the optimized 3-wave coupling coefficients were significantly reduced.

The second optimization problem, as stated in Eq. (18), converged after 102 iterations. Again, most constraints were fulfilled in fewer iterations and the 3-wave coupling coefficient was the limiting factor in the optimization. Figure 3b shows a rather smooth convergence of the 3-wave coupling coefficient towards the target value. The optimized quarter layout in Fig. 5 exhibits significant morphing. Note that the shape changes are not just limited to variations in the width but also show significant bending of the initially straight springs.

It appears that, in our examples, decreasing an initially large 3-wave coupling coefficient requires fewer iterations and smaller design changes than increasing an initially small 3-wave coupling coefficient. In both optimizations we find unintuitive design changes which have a significant impact on the optimized 3-wave coupling coefficients. This demonstrates the usefulness of gradient-based shape optimization for tuning geometrically nonlinear modal couplings without manual topology adjustments.

To investigate the sensitivity of the optimized designs to fabrication tolerances, we also calculate the coupling coefficients of the initial and optimized designs for a dilated and eroded geometry. This is a common approach when modelling MEMS or NEMS resonators, where the tolerances of the employed etching processes can result in variations of the in-plane width of structures22,39. For this purpose, we apply a uniform shift of +50 nm (dilated design) or -50 nm (eroded design) to all exterior nodes along their outward surface normal direction in the xy-plane. For example, a +50 nm shift will widen a spring by +100 nm, due to it being applied on both sides of the spring, which is a typical order of magnitude for fabrication tolerances in MEMS. The results for the minimized coupling coefficients, as shown in Fig. 3a, and the maximized coupling coefficient, as shown in Fig. 3b, are summarized in Fig. 6. The results are shown for the initial and the optimized designs for a dilated (thicker), eroded (thinner) and nominal geometry. The initial nominal geometry corresponds to Fig. 1 and the optimized nominal geometries correspond to Figs. 4 and 5. One can see a clear impact of the fabrication tolerances on the optimized 3-wave coupling coefficients, whereas the effect in the initial design is rather small. The coupling coefficients that were reduced, as shown in Fig. 3a, increase by roughly one order of magnitude in the optimized design for the eroded and dilated geometry compared to the nominal design, as seen in Fig. 6. Nevertheless, they still remain around two orders of magnitude below the initial values. The coupling coefficient that was increased in the second optimization, as seen in Fig. 3b, changes by less than one order of magnitude for the eroded and dilated geometry compared to the nominal design, as seen in Fig. 6. Its value remains roughly two orders of magnitude larger than in the initial design.

Our results highlight that even with fabrication errors the nonlinear coupling coefficients can be tuned over various orders of magnitude. Note that we did not consider the fabrication errors during the optimization. A common approach is to evaluate the dilated, nominal and eroded design and to optimize the worst case out of the three designs in each iteration22,39. Fabricating and experimentally verifying such structures will be subject of future research. For that purpose, it would be advantageous to consider the fabrication errors during the optimization, to obtain a robust design.

In this work, we developed a framework for the node-based shape optimization of geometrically nonlinear 3-wave coupling coefficients, and applied it to the design of a MEMS gyroscope. It was demonstrated that 3-wave coupling coefficients can be increased or decreased by several orders of magnitude with our approach, while also fulfilling manufacturability constraints and constraints on eigenfrequencies. The 3-wave coupling coefficients between in-plane, out-of-plane and torsional modes could be reduced by 3 orders of magnitude with very subtle design changes. The 3-wave coupling coefficient leading to an 1:2 internal resonance between an in-plane and an out-of-plane mode could be increased by a factor of 250, while also enforcing the corresponding resonance condition. The optimizations lead to design changes which are non-intuitive even to experienced human designers. Furthermore, it was estimated that, even in the presence of fabrication errors, the optimized coupling coefficients were still around two orders of magnitude smaller or larger than the initial values. Our results highlight that shape optimization is a powerful design tool for MEMS and NEMS resonator applications which exploit nonlinearities as well as applications where nonlinearities lead to parasitic effects. In the future, it would be intriguing to fabricate and characterize test structures with optimized geometrically nonlinear 3-wave couplings, to experimentally prove the applicability of our approach. For that purpose it would also be benficial to consider fabrication errors during optimization to obtain robust designs.

All the data are available in the manuscript, or they will be provided upon reasonable request.

Bachtold, A., Moser, J. & Dykman, M. Mesoscopic physics of nanomechanical systems. Rev. Mod. Phys. 94, 045005. https://doi.org/10.1103/RevModPhys.94.045005 (2022).

Article ADS MathSciNet CAS MATH Google Scholar

Eichler, A. & Zilberberg, O. Classical and Quantum Parametric Phenomena (Oxford University Press, 2023).

Book MATH Google Scholar

Bukhari, S. A. R., Saleem, M. M., Hamza, A. & Bazaz, S. A. A novel design of high resolution mems gyroscope using mode-localization in weakly coupled resonators. IEEE Access 9, 157597–157608 (2021).

Article Google Scholar

Miao, T. et al. Nonlinearity-mediated digitization and amplification in electromechanical phonon-cavity systems. Nat. Commun. 13, 2352 (2022).

Article ADS CAS PubMed PubMed Central MATH Google Scholar

Nayfeh, A. H. & Mook, D. T. Nonlinear Oscillations (Wiley, 2008).

MATH Google Scholar

Kerschen, G. (ed.) Modal Analysis of Nonlinear Mechanical Systems Vol. 555 (Springer, 2014).

MATH Google Scholar

Touzé, C. & Amabili, M. Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures. J. Sound Vib. 298, 958–981. https://doi.org/10.1016/j.jsv.2006.06.032 (2006).

Article ADS MATH Google Scholar

Opreni, A., Vizzaccaro, A., Frangi, A. & Touzé, C. Model order reduction based on direct normal form: Application to large finite element MEMS structures featuring internal resonance. Nonlinear Dyn. 105, 1237–1272. https://doi.org/10.1007/s11071-021-06641-7 (2021).

Article MATH Google Scholar

Sarrafan, A., Azimi, S., Golnaraghi, F. & Bahreyni, B. A nonlinear rate microsensor utilising internal resonance. Sci. Rep. 9, 8648. https://doi.org/10.1038/s41598-019-44669-3 (2019).

Article ADS CAS PubMed PubMed Central Google Scholar

Nitzan, S. H. et al. Self-induced parametric amplification arising from nonlinear elastic coupling in a micromechanical resonating disk gyroscope. Sci. Rep. 5, 9036. https://doi.org/10.1038/srep09036 (2015).

Article CAS PubMed PubMed Central MATH Google Scholar

Antonio, D., Zanette, D. H. & Lopez, D. Frequency stabilization in nonlinear micromechanical oscillators. Nat. Commun. 3, 806. https://doi.org/10.1038/ncomms1813 (2012).

Article ADS CAS PubMed MATH Google Scholar

Nabholz, U., Lamprecht, L., Mehner, J. E., Zimmermann, A. & Degenfeld-Schonburg, P. Parametric amplification of broadband vibrational energy harvesters for energy-autonomous sensors enabled by field-induced striction. Mech. Syst. Signal Process. 139, 106642. https://doi.org/10.1016/j.ymssp.2020.106642 (2020).

Article MATH Google Scholar

Ganesan, A., Do, C. & Seshia, A. Phononic frequency comb via intrinsic three-wave mixing. Phys. Rev. Lett. 118, 033903. https://doi.org/10.1103/PhysRevLett.118.033903 (2017).

Article ADS PubMed MATH Google Scholar

Nabholz, U., Curcic, M., Mehner, J. E. & Degenfeld-Schonburg, P. Nonlinear Dynamical System Model for Drive Mode Amplitude Instabilities in MEMS Gyroscopes https://doi.org/10.1109/ISISS.2019.8739703 (2019).

Nabholz, U., Schatz, F., Mehner, J. E. & Degenfeld-Schonburg, P. Spontaneous parametric down-conversion induced by non-degenerate three-wave mixing in a scanning MEMS micro mirror. Sci. Rep. 9, 3997. https://doi.org/10.1038/s41598-019-40377-0 (2019).

Article ADS CAS PubMed PubMed Central MATH Google Scholar

Nabholz, U., Stockmar, F., Mehner, J. E. & Degenfeld-Schonburg, P. Validating the critical point of spontaneous parametric downconversion for over 600 scanning MEMS micromirrors on wafer level. IEEE Sensors Lett. 4, 1–4. https://doi.org/10.1109/LSENS.2020.2964384 (2020).

Article MATH Google Scholar

Lifshitz, R. & Cross, M. C. Nonlinear dynamics of nanomechanical and micromechanical resonators. Rev. Nonlinear Dyn. Complex. https://doi.org/10.1002/9783527626359.ch1 (2008).

Article MATH Google Scholar

Giannini, D., Bonaccorsi, G. & Braghin, F. Size optimization of MEMS gyroscopes using substructuring. Eur. J. Mech. A. Solids 84, 104045. https://doi.org/10.1016/j.euromechsol.2020.104045 (2020).

Article MathSciNet MATH Google Scholar

Giannini, D., Braghin, F. & Aage, N. Topology optimization of 2D in-plane single mass MEMS gyroscopes. Struct. Multidiscip. Optim. 62, 2069–2089. https://doi.org/10.1007/s00158-020-02595-3 (2020).

Article MathSciNet MATH Google Scholar

Giannini, D., Aage, N. & Braghin, F. Topology optimization of MEMS resonators with target eigenfrequencies and modes. Eur. J. Mech. A. Solids 91, 104352. https://doi.org/10.1016/j.euromechsol.2021.104352 (2022).

Article MathSciNet MATH Google Scholar

Shin, D. et al. Spiderweb nanomechanical resonators via Bayesian optimization: Inspired by nature and guided by machine learning. Adv. Mater. 34, 2106248. https://doi.org/10.1002/adma.202106248 (2022).

Article CAS Google Scholar

Hoj, D. et al. Ultra-coherent nanomechanical resonators based on inverse design. Nat. Commun. 12, 5766. https://doi.org/10.1038/s41467-021-26102-4 (2021).

Article ADS CAS PubMed PubMed Central MATH Google Scholar

Pedersen, N. Designing plates for minimum internal resonances. Struct. Multidiscip. Optim. 30, 297–307. https://doi.org/10.1007/s00158-005-0529-x (2005).

Article MathSciNet MATH Google Scholar

Schiwietz, D., Hörsting, M., Weig, E. M., Degenfeld-Schonburg, P. & Wenzel, M. Shape Optimization of Eigenfrequencies in MEMS Gyroscopes. https://doi.org/10.48550/ARXIV.2402.05837 (2024).

Dou, S., Strachan, B. S., Shaw, S. W. & Jensen, J. S. Structural optimization for nonlinear dynamic response. Philos. Trans. R. Soc.y A: Math. Phys. Eng. Sci. 373, 20140408. https://doi.org/10.1098/rsta.2014.0408 (2015).

Article ADS MathSciNet MATH Google Scholar

Li, L. L. et al. Tailoring the nonlinear response of MEMS resonators using shape optimization. Appl. Phys. Lett. 110, 081902. https://doi.org/10.1063/1.4976749 (2017).

Article ADS CAS MATH Google Scholar

Zega, V., Langfelder, G., Falorni, L. G. & Comi, C. Hardening, softening, and linear behavior of elastic beams in mems: An analytical approach. J. Microelectromech. Syst. 28, 189–198 (2019).

Article MATH Google Scholar

Wriggers, P. Nonlinear Finite Element Methods (Springer, 2008).

MATH Google Scholar

Touze, C., Vizzaccaro, A. & Thomas, O. Model order reduction methods for geometrically nonlinear structures: A review of nonlinear techniques. Nonlinear Dyn. 105, 1141–1190. https://doi.org/10.1007/s11071-021-06693-9 (2021).

Article MATH Google Scholar

Svanberg, K. The method of moving asymptotes-a new method for structural optimization. Int. J. Numer. Meth. Eng. 24, 359–373. https://doi.org/10.1002/nme.1620240207 (1987).

Article MathSciNet MATH Google Scholar

Lee, T. H. An adjoint variable method for structural design sensitivity analysis of a distinct eigenvalue problem. KSME International Journal 13, 470–476. https://doi.org/10.1007/BF02947716 (1999).

Article ADS MATH Google Scholar

Tcherniak, D. Topology optimization of resonating structures using SIMP method. Int. J. Numer. Meth. Eng. 54, 1605–1622. https://doi.org/10.1002/nme.484 (2002).

Article MATH Google Scholar

Nelson, R. B. Simplified calculation of eigenvector derivatives. AIAA J. 14, 1201–1205. https://doi.org/10.2514/3.7211 (1976).

Article ADS MathSciNet MATH Google Scholar

Putnik, M., Cardanobile, S., Nagel, C., Degenfeld-Schonburg, P. & Mehner, J. Simulation and modelling of the drive mode nonlinearity in MEMS-gyroscopes. Procedia Eng. 168, 950–953. https://doi.org/10.1016/j.proeng.2016.11.313 (2016).

Article MATH Google Scholar

Putnik, M. et al. Incorporating geometrical nonlinearities in reduced order models for MEMS gyroscopes https://doi.org/10.1109/ISISS.2017.7935656 (2017).

Article Google Scholar

Putnik, M. et al. A static approach for the frequency shift of parasitic excitations in MEMS gyroscopes with geometric nonlinear drive mode https://doi.org/10.1109/ISISS.2017.7935651 (2017).

Article Google Scholar

Putnik, M. et al. Predicting the resonance frequencies in geometric nonlinear actuated MEMS. J. Microelectromech. Syst. 27, 954–962. https://doi.org/10.1109/JMEMS.2018.2871080 (2018).

Article CAS MATH Google Scholar

Putnik, M. et al. Simulation methods for generating reduced order models of MEMS sensors with geometric nonlinear drive motion https://doi.org/10.1109/ISISS.2018.8358112 (2018).

Article Google Scholar

Sigmund, O. Manufacturing tolerant topology optimization. Acta. Mech. Sin. 25, 227–239 (2009).

Article ADS MATH Google Scholar

Download references

The IPCEI ME/CT project is supported by the Federal Ministry for Economic Affairs and Climate Action on the basis of a decision by the German Parliament, by the Ministry for Economic Affairs, Labor and Tourism of Baden-Württemberg based on a decision of the State Parliament of Baden-Württemberg, the Free State of Saxony on the basis of the budget adopted by the Saxon State Parliament, the Bavarian State Ministry for Economic Affairs, Regional Development and Energy and financed by the European Union - NextGenerationEU. The authors are thankful to Martin Putnik at Robert Bosch GmbH for providing the model of the initial design.

Open Access funding enabled and organized by Projekt DEAL.

Robert Bosch GmbH, Corporate Sector Research and Advance Engineering, 71272, Renningen, Germany

Daniel Schiwietz, Marian Hörsting, Matthias Wenzel & Peter Degenfeld-Schonburg

Technical University of Munich, School of Computation, Information and Technology, 80333, Munich, Germany

Daniel Schiwietz & Eva Maria Weig

Technical University of Munich, Center for Quantum Engineering (ZQE), 85748, Garching, Germany

Eva Maria Weig

You can also search for this author inPubMed Google Scholar

You can also search for this author inPubMed Google Scholar

You can also search for this author inPubMed Google Scholar

You can also search for this author inPubMed Google Scholar

You can also search for this author inPubMed Google Scholar

D.S. and P.D. crafted the original idea. D.S. implemented the method into code. D.S., M.H. and M.W. developed the underlying shape optimization tool. D.S. performed the optimization. D.S. wrote the manuscript. P.D. and E.W. supervised the work. All authors reviewed the manuscript.

Correspondence to Daniel Schiwietz.

The authors declare no competing interests.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

For the geometrically nonlinear 3-wave coupling coefficients, defined in Eq. (5), we obtain the partial derivatives for the adjoint method as

where \(\varvec{J}^e\) is the Jacobian matrix of element e, p is the integration point index, \(W_p\) is the integration point weight, \(\varvec{L}^e\) is a matrix which maps the element eigenvector to the global eigenvector as \(\varvec{\phi }_i^e=\varvec{L}^e\varvec{\phi }_i\) and \(\delta _{i,j}\) is the Kronecker delta. To derive Eq. (21), we used the property \(\varvec{B}^e_{\eta }\left( \varvec{\phi }^e_i\right) \varvec{\phi }^e_j=\varvec{B}^e_{\eta }\left( \varvec{\phi }^e_j\right) \varvec{\phi }^e_i\). The required sensitivities of \(\varvec{B}_\varepsilon ^e\) and \(\det (\varvec{J}^e)\) are defined as in24. Note that \(\varvec{B}^e_{\eta }\left( \varvec{\phi }^e_i\right) =\varvec{B}^e_{\eta }\left( \varvec{H}_i^e,\varvec{B}^e\right)\) is bilinear in \(\varvec{H}_i^e\) and \(\varvec{B}^e\), with \(\varvec{H}_i^e\) being the displacement gradient of the eigenvector of mode i inside element e and \(\varvec{B}^e\) containing shape function derivatives. The displacement gradient and its partial derivative can be obtained as

where \(\varvec{\varPhi }_i^e\) is a matrix in which each column corresponds to one node of element e and contains the 3 eigenvector components belonging to that node for mode i. The matrices \(\varvec{B}^e\) and \(\frac{\partial \varvec{B}^e}{\partial p_j}\) are defined as in24. Exploiting the bilinearity of \(\varvec{B}^e_{\eta }\left( \varvec{\phi }^e_i\right)\), its partial derivative is obtained via the product rule as

After evaluating Eqs. (1920)–(21) for all relevant \(\alpha _{n,m,l}\), the corresponding partial derivatives of the in Section 4.2 introduced \(|{\bar{\alpha }}_{d,1,9}|\), \(|{\bar{\alpha }}_{d,7,9}|\) and \(|{\tilde{\alpha }}_{d,d,7}|\) are obtained via the chain rule and their total sensitivities can be evaluated using the adjoint method.

The sensitivity of the constraint on the drive mode’s eigenvector is required. Note that \(\varvec{\phi }_{d,x}\) and \(\varvec{\phi }_{d,y}\) have the same dimensions as \(\varvec{\phi }_{d}\) but simply contain zeros everywhere, except for the x or y degrees of freedom inside \(V_{el}\), where they contain the corresponding values of \(\varvec{\phi }_{d}\). The sensitivity is calculated via the adjoint method and the required partial derivatives are

We employ a similar approach to mode tracking as in24. We define a modal assurance criterion (MAC) as

where \(\varvec{R}\) is the upper triangular matrix obtained from the Cholesky decomposition \(\varvec{M}=\varvec{R}^T\varvec{R}\), A and B are weights, the superscript \(\{k\}\) refers to quantities calculated in the current iteration k, \(\{k-1\}\) refers to the previous iteration and \(\{0\}\) to the initial design. The mode i for which \(\text {MAC}^{\{k\}}_{ij}\) has the largest magnitude determines the mode number in iteration k corresponding to mode j from the previous iteration. We use \(A=0.1\) and \(B=0.9\), which we found to ensure a robust identification of the initially defined modes in this work, even if modes mix and demix throughout the optimization.

Figures 7 and 8 show the modes of the two optimized designs, corresponding to the modes in Fig. 2. Note that the order of the modes has changed in the optimized designs, when sorted by eigenfrequencies. However, for comparability, we label them in the same order as in the initial design.

Mode shapes of the optimized design according to Eq. (17). For simplicity, only the displacement of the surface is shown. The coloring indicates the sum of the three displacement components at each point, normalized to each mode’s maximum value. We refer to the modes by their numbers in the initial design. (a) Mode 1 with an eigenfrequency of 21.4 kHz. (b) Mode 3 (drive) with an eigenfrequency of 24.7 kHz. (c) Mode 4 (sense) with an eigenfrequency of 27.1 kHz. (d) Mode 7 with an eigenfrequency of 49.8 kHz. (e) Mode 9 with an eigenfrequency of 52.5 kHz. (f) Mode 13 with an eigenfrequency of 105.5 kHz.

Mode shapes of the optimized design according to Eq. (18). For simplicity, only the displacement of the surface is shown. The coloring indicates the sum of the three displacement components at each point, normalized to each mode’s maximum value. We refer to the modes by their numbers in the initial design. (a) Mode 1 with an eigenfrequency of 15.8 kHz. (b) Mode 3 (drive) with an eigenfrequency of 24.7 kHz. (c) Mode 4 (sense) with an eigenfrequency of 27.1 kHz. (d) Mode 7 with an eigenfrequency of 49.4 kHz. (e) Mode 9 with an eigenfrequency of 45.3 kHz. (f) Mode 13 with an eigenfrequency of 51.4 kHz.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

Schiwietz, D., Hörsting, M., Weig, E.M. et al. Shape optimization of geometrically nonlinear modal coupling coefficients: an application to MEMS gyroscopes. Sci Rep 15, 10957 (2025). https://doi.org/10.1038/s41598-025-95412-0

Download citation

Received: 01 September 2024

Accepted: 19 March 2025

Published: 31 March 2025

DOI: https://doi.org/10.1038/s41598-025-95412-0

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative