A temperature-robust level-set approach for eigenfrequency optimization

The optimization of target eigenfrequencies is crucial for several engineering applications, including dynamical systems. Micro-electro-mechanical systems (MEMS) used in time-keeping applications, for example, require exceptional frequency stability. Most eigenfrequency structural optimization methods focus on a deterministic approach, often neglecting potential fluctuations in operational conditions. Among these, temperature variations have long been known to have a detrimental effect on the natural frequencies of a structure. In this work, we show how eigenfrequency optimization can be applied to the field of structural dynamics while minimizing the variance of natural frequencies caused by external temperature uncertainties. To accomplish this, we employ a level-set optimization algorithm, known for its computational efficiency and ability to define crisp interfaces.


Introduction
Several engineering applications require or greatly benefit from the optimization of natural frequencies. This is usually done for performance reasons, safety improvement, disturbance rejection, human comfort, or noise reduction.
A few examples of such applications in the field of structural dynamics are, e.g., resonators, bridges, energy harvesters, or car chassis (Nobari et al. 2015;Dunning et al. 2016;Giannini et al. 2022). Several techniques exist to optimize such structures, depending on the problem at hand and on the desired goal. In cases where the design is constrained to adopt a prescribed geometry, a parametric-or size-optimization (Choi and Kim 2005;Giannini et al. 2020) is usually recommended.
On the other hand, topology optimization techniques (Bendsoe and Sigmund 2004;Dunning and Kim 2015) can attain any suitable shape within the design space, thus increasing the overall freedom of the design process. There exist several approaches to topology optimization, such as Solid Isotropic Material with Penalization (SIMP) method (Bendsøe and Sigmund 1999), level-set method (Wang et al. 2003;Allaire et al. 2004), Evolutionary Structural Optimization (ESO) method (Xie and Steven 1993), and the more recent Moving Morphable Components (MMC) method (Guo et al. 2014). Among these, level-set optimization has been successfully adopted to solve a number of engineering problems, mainly due to its ability to define sharp interfaces without filtering stages, and to its computational efficiency derived from the fact that only the boundary is propagated, instead of the whole physical domain.
Despite the successful application of topology optimization in all these fields, the great majority of existing works focus on the deterministic performance of the structures at hand. While this leads to optimal structures under nominal conditions, their performance can be greatly affected by a change in loading conditions, working conditions, or environmental factors, among others.
Robust formulations have been adopted for structures subject to varying loads or uncertainties on the loading point. For instance, Dunning et al. (2011) proposed an approach to minimize the expected compliance, whereas Dunning and Kim (2013); Xu et al. (2019) presented algorithms that consider both the expected value and the variance of the compliance. On the other hand, Guest and Igusa (2008) applied a multi-load approach to truss structures subject to Gaussian uncertainties in the external loads and extended the method to consider uncertain nodal locations, whereas Chen et al. (2010) perform the optimization under random field uncertainties.
Similarly, there are works presenting formulations able to keep into account the variance of a manufacturing process, which may lead to over-or under-etching of a nominal structure geometry. For instance, Wang et al. (2011) introduced filtering stages to ensure mesh convergence at both global and local scales. These projection stages are used to optimize compliant mechanisms and heat conductors under a minimum length scale. In a similar way, Jang et al. (2012); Andreasen et al. (2020) exploit these filtering stages to implement a worst-case approach for compliance minimization, whereas Li et al. (2019) dealt with frequency maximization in the presence of possible over-etch. All these methods follow a worst-case approach, which sometimes leads to solutions that are too conservative. On the other hand, probabilistic methods, such as the Monte Carlo method (Rubinstein and Kroese 2016) or polynomial chaos expansion (Wiener 1938), are able to better describe the manufacturing uncertainty field, thus allowing the computation of the response statistics and their sensitivities. Among these are the works of Schevenels et al. (2011), Tootkaboni et al. (2012, Lazarov et al. (2012), Chen and Chen (2011) and Zhang and Kang (2017).
Concerning structural optimization in thermal environments, Chung et al. (2020) presented a level-set topology optimization technique for nonlinear thermoelasticity, whereas Kambampati et al. (2020) deals with stress-based topology optimization in presence of a spatially varying temperature field. Topology optimization approaches for thermal buckling criteria were presented by Wu et al. (2019), Gan and Wang (2022). Regarding dynamic problems, Li (2013, 2014) focused on dynamic compliance minimization.
The problem of eigenfrequencies stability in temperature is extremely relevant in several applications, with literature providing numerous theoretical and experimental studies of the phenomenon (Talebian et al. 2010;Wang et al. 2012). Perfect examples of this are microelectromechanical systems (MEMS). For instance, MEMS for time-keeping applications (Wu et al. 2020) are crucial in the electronics industry and represent a multi-billion dollar market (van Beek and Puers 2011). While resonators have historically been produced using quartz crystals, the recent advances in MEMS technology allowed us to obtain reliable, miniaturized, and cheap MEMS resonators and oscillators.
However, frequency drifts are among the biggest challenges that MEMS designers have to face, whether they are dealing with gyroscopes, oscillators, or micro-mirrors (Jiang et al. 2021;Ferguson et al. 2005;Wolter et al. 2005).
Even though specific MEMS architectures and/or feedback loops exist to keep the oscillating frequency at the desired value (Prache et al. 2016;Salvia et al. 2010), it would be beneficial to have layouts intrinsically robust to temperature variations. This, in fact, would keep the oscillating frequency as close as possible to the system's resonance, thus reducing the amount of forcing needed to keep the device in motion. However, to the authors' knowledge, there are no works in literature that present techniques able to reduce or minimize the variance of a system's natural frequencies in presence of a varying temperature.
The rest of the paper is organized as follows. Section 2 briefly presents the most common existing approaches for reducing the temperature effects on MEMS resonators. After that, Sect. 3 shows how the thermal behavior of the material properties and the thermal expansion phenomenon affect the natural frequencies and the mode shapes of a structure. Then, the temperature-robust eigenfrequency optimization method is described, along with the statistical analysis and the sensitivity derivation. Different numerical examples of the proposed approach are presented in Sect. 5, including the application to a 1 MHz MEMS resonator. Finally, Sect. 6 presents the conclusions.

Existing approaches for reducing temperature effects
In the field of MEMS resonators for timing applications, the natural frequency sensitivity with respect to temperature variations is expressed through the Temperature Coefficient of Frequency (TCF, Wu et al. 2020): where f 0 is the nominal value of the natural frequency. This term is usually measured in ppm/°C. Typically, silicon resonators exhibit a negative TCF around −30 ppm/°C, which is too high for industrial applications. Therefore, compensation techniques are required to achieve a low TCF and a highly stable timing reference. The most common temperature compensation techniques are based either on passive or active approaches (Wu et al. 2020). An example of active compensation is electrostatic tuning, which exploits the electrostatic stiffness effect obtained through a bias voltage in capacitive MEMS. The bias voltage is carefully tuned to compensate for temperature-induced frequency variations. A TCF of 39 ppm over a temperature range between 25 • C and 125 • C has been achieved by Sundaresan et al. (2007). However, this approach is not always feasible since it requires large bias voltages which may not be suitable for standard CMOS. Another approach is electronic compensation, which is typically implemented through a Phase-Locked Loop (PLL) circuit coupled with a temperature sensor. Even though the PLL circuit inevitably increases the phase noise and the jitter of the output clock, this approach has been used to design a MEMS-based programmable oscillator with a TCF less than 0.1 ppm over a temperature range between −45 and 105 °C (Roshan et al. 2016). A different approach is a micro-oven control, where a MEMS resonator is placed inside a thermally isolated micro-oven whose temperature is controlled with a feedback circuit. For instance, the Oven-Controlled MEMS Oscillator (OCMO,  is characterized by a TCF of around 0.3 ppm in the range of −45 and 85 • C. Compared to active techniques, passive compensation approaches do not require additional power consumption or dedicated circuits. For instance, composite structures exploit the positive TCF of some materials, such as SiO 2 , to compensate for the negative TCF of silicon. Similarly, doping (Csavinszky and Einspruch 1963) allows for a decrease in the TCF of Silicon through heavy doping. Besides these techniques, which rely on material engineering, another passive approach is geometry engineering. The idea is to design a layout where the temperature dependence of the material properties is compensated by the thermal stresses, yielding a device with a low TCF. For instance, Hsu and Nguyen (1998) applies this approach to a low-frequency foldedbeam MEMS device, whereas Hsu et al. (2000) designed a 10 MHz MEMS resonator with a TCF around −2.5 °C.
Typically, in industrial applications, a combination of active and passive compensation approaches is used. In particular, composite structures or doping is used to greatly reduce the initial TCF. Then, electronic compensation or micro-oven control is implemented to further improve the frequency robustness of such devices.

Temperature effects on the natural frequencies
In this Section, we show how considering only elasticity is not sufficient when optimizing a structure for frequency robustness in temperature, and it is necessary to take into account thermal stresses as well.

Influence of the material properties
The characteristic equation of a linear second-order vibrating system is where M and K are respectively the mass and stiffness matrices, q is the mode shape vector, and is the natural frequency. Under the assumption of linear elastic material, the matrices M and K linearly depend on the material properties. Given that, for the applications considered (i.e., MEMS resonators), the frequency drift resulting from temperature changes is primarily due to the thermal behavior of Young's modulus (Jiang et al. 2021) and that thermal stresses affect the K but not M (Wu et al. 2020), this study focuses on the changes related to the stiffness matrix. For most isotropic materials, there is a wide range of temperature over which the change of Young's modulus can be approximated to be linear (Ledbetter 1982), where E 0 is the initial value for E, ΔT is the temperature variation, and T is the linear thermal coefficient (Ledbetter 1982).
Since the matrix K linearly depends on Young's modulus, Eq.
(2) is rewritten as: where the coefficient accounts for the thermal behavior of Young's modulus. Therefore, a reference angular frequency 0 is affected by a change in temperature according to In particular, in the case of silicon, Young's modulus E tends to decrease (soften) with increasing temperatures (Shirai 2013), leading to a subsequent frequency decrease.
That is, if only Eq. (4) was considered, structures made of the same material but characterized by different natural frequencies at the same temperature will experience the same relative variation ∕ 0 . Moreover, the mode shape associated with the natural frequency remains unaffected by the temperature change. Therefore, in order to introduce a dependency on the shape of the structure, thermal stresses need to be considered.

Influence of the thermal stress
A common approach to include the thermal stresses is to add to Eq. (2) the geometric stiffness matrix (Crisfield 1996;Liu et al. 2013): The matrix depends on the nodal displacement vector u T , obtained from the solution of the static thermal expansion problem, where f T is the equivalent force vector, proportional to the thermal expansion coefficient k , that accounts for the thermal stresses.
It is possible to show that f T linearly depends on the temperature variation ΔT (Crisfield 1996;Liu et al. 2013): Consequently, it is possible to write the normalized thermal expansion problem: It can be demonstrated (Crisfield 1996;Liu et al. 2013) that linearly depends on ΔT as well: Moreover, linearly depends also on Young's modulus. Therefore, considering both the thermal behavior of Young's modulus (Eq. 3) and the thermal stresses (Eq. 9), the characteristic equation (Eq. 2) can be rewritten as In this work, the mass normalization conditions of the mode shapes are used: The formulation in Eq. (11) is suitable for a topology optimization approach and will be used in the rest of this work.

Level-set topology optimization
In this Section, we show the statistical analysis used to compute the eigenfrequency expected value and variance. After that, we illustrate the temperature-robust level-set topology optimization algorithm.

Statistical analysis
The temperature difference ΔT is assumed to follow a Gaussian distribution characterized by a mean value and variance 2 . The value of ΔT affects the eigenfrequency according to Eq. (11). By applying the Law of the Unconscious Statistician (LOTUS) (Allen 2006, Section 18.3.4), the expected value of the natural frequency [ ] can be computed as follows: where P(ΔT) is the Normal Probability Density Function associated with ΔT . Then, the variance is computed as To solve Eqs. (13) and (14), the natural frequency is linearized around ΔT = : A simple way to approximate T is through finite differences. However, this approach would be associated with a high computational cost since it requires solving additional eigenvalue problems. For instance, the first-order forward and backward schemes require an additional solution. On the other hand, if the central difference scheme is used, two additional eigenvalue problems are required (at ΔT = ± ). Therefore, in this work, an analytical approach is used.
The expression of T is computed by differentiating Eq. (11): Using the mass normalization conditions Eqs. (12) and (16) can be rewritten as: where the coefficient m T is defined as Finally, Eq. (17) allows to compute Using the previous results, Eqs. (13) and (14) lead:

Temperature-robust eigenfrequency optimization
Equations (20) and (21) allow to define the temperaturerobust eigenfrequency problem as follows: where the objective function is a linear combination of [ ] and 2 [ ] , A indicates the area fraction (i.e., the fraction of the available design domain that is occupied by the material), and A max refers to the maximum acceptable area fraction. This formulation aims at maximizing the expected natural frequency while minimizing its variance. The coefficients w e and w s are computed according to Asadpoure et al. (2011): where 0 ≤ w ≤ 1 , and 0 is the eigenfrequency of the initial design. If w = 1 , only the expected value is maximized and the robust problem falls back into the deterministic case. On the other hand, if w = 0 , only the variance is minimized. To track the desired mode shape during the optimization process, the Modal Assurance Criterion (MAC) is implemented (Kim and Kim 2000;Pastor et al. 2012).
Problem (22) is solved iteratively using a level-set optimization method (Algorithm 1). Just like in the classical Level-Set Topology Optimization (LSTO, Allaire et al. 2002;Wang et al. 2003;Allaire et al. 2004), the Level-Set Method (LSM, Osher and Sethian 1988;Sethian 1999;Osher and Fedkiw 2006) is exploited to track and evolve the boundary of the structure throughout the robust optimization process (Fig. 1). The level-set function Φ is propagated according to the Hamilton-Jacobi equation, where V n (x) indicates the velocity field, obtained by extending the optimal boundary points velocities to the entire design domain (Osher and Fedkiw 2006). At each iteration, the structure identified by the zero level-set of Φ is discretized using the Marching Squares algorithm (Lorensen and Cline 1987;Maple 2003). This allows to compute the boundary points z j , which are the design variables of the optimization problem, as well as the area fractions k of all the grid elements (Hyun and Kim 2021b). Figure 2 provides a graphical representation of the Marching Squares discretization. The area fractions are then used to define the structural matrices according to the Ersatz material approximation with a linear interpolation scheme for the material properties (Allaire et al. 2004).
The boundary point sensitivities are computed using the approach presented by Kambampati et al. (2021): the derivatives of [ ] and 2 [ ] with respect to the area fractions k of individual elements (Sect. 4.3) are combined with a local perturbation of the boundary points to obtain the shape sensitivities at the boundary points (Chen and Chen 2011;Wu et al. 2021). The sensitivities at the boundary points are used in the sub-optimization problem (Allaire et al. 2004) that is solved with IPOPT (Interior Point OPTimizer, Wächter and Biegler 2006), an open-source library for large-scale nonlinear constrained optimization.
The optimal velocities of the boundary points, which come from the solution to the sub-optimization problem, are then used to update the level-set function. This operation is performed by solving the level-set equation in its discrete form, where Φ is the level-set function, k refers to the iteration number, i indicates a grid node, Δt is the time step of the propagation, and V n,i indicates the velocity of node i (normal to the geometry boundary).
In this work, a fixed grid made of unitary square elements is used, along with the Hamilton-Jacobi Weighted Essentially Non-Oscillatory method (HJ-WENO, Osher and Fedkiw 2006) for spatial discretization, and a Forward Euler approach for time discretization. In addition, Δt = 1 is used, leading to the following Courant-Friedrichs-Lewy (CFL) condition for numerical stability (Osher and Fedkiw 2006): Throughout this work, CFL = 0.2 has proven to be a good trade-off between convergence and numerical stability.
Convergence is reached when all constraints are satisfied and the objective (or cost) variation between two successive iterations is sufficiently small in terms of relative difference.

Sensitivity analysis
In this Section, we present the analytical procedure used to to obtain the derivatives of the expected natural frequency [ ] and its variance 2 [ ] with respect to the area fractions k of the grid elements. These derivatives are then combined with local perturbation of the boundary points (Kambampati et al. 2021) to compute the boundary sensitivities of Problem (22).
While the derivative of the expected natural frequency [ ] k is derived by following an approach similar to that presented in Li (2013, 2014), the derivative of the frequency variance requires a novel formulation that is presented later in this Section. The derivative of [ ] with respect to k is written as where ũ T,i is the i-th element of ũ T that is obtained by solving Eq. (9). The term within square brackets in Eq.(28), can be rewritten in a compact form as where the operator is defined as To obtain the derivative of T with respect to k , we start by differentiating Eq. (17): where the term dd k is the total derivative with respect to k .
Computing q k for all the grid elements would be computationally expensive. In order to avoid this operation, the adjoint method is used by exploiting the Static Modal Derivatives (SMD, Idelsohn and Cardona 1985a, b;Slaats et al. 1995;Varona et al. 2019) of Eq. (11) with respect to k : Combining Eqs. (35) and (36): To avoid the computation of q k , we solve the adjoint equation which allows to simplify Eq. (37) as follows: where y = q + ΔT . Using the operator introduced in Eq. (30), it is possible to show that: As before, the adjoint method is used to avoid the computation of the term ũ T k . Equation (9)

Hence, Eq. (42) simplifies as
Finally, combining Eq. (44) and (39) and rearranging the terms, the sensitivity of T is

Numerical examples
In this Section, we present the results of the temperaturerobust eigenfrequency optimization. To show the validity of the approach, Problem (22) is solved at different temperature conditions ( and ) with different weights w . After that, we show how to slightly modify the formulation of Problem (22) to design a thermally robust resonator with a specific natural frequency. The parameters used in the following examples are reported in Table 1. The material considered is polysilicon, whose mechanical properties as a function of temperature have been intensely studied by Sharpe et al. (2001). The (43) K = − (y, q).
(44) problem settings, the boundary conditions, and the initial layout are shown in Fig. 3. All the results are obtained using OpenLSTO (Kambampati et al. 2018).

Thermally robust frequency maximization
The aim of this example is to maximize the natural frequency associated with the oscillation of the proof mass along the y axis (Fig. 3a). The approach is applied to the same problem with different combinations of mean temperature , temperature standard deviation , and weight w . The numerical results are reported in Table 2, whereas Fig. 4 shows two optimal layouts obtained at the same temperature condition but with two different weights.
By using a value of w close to 1, the algorithm converges to a solution with stiffer elements, thus increasing the expected natural frequency as well as its variance. On the other hand, if a low value of w is used, the optimal layout presents more flexible elements that allow the residual thermal stresses to compensate for the frequency variation. This results in a reduction of both its variance and the expected value. When w = 1 , there is no control over the frequency variance.

Design of a thermally robust resonator
In this Section, we show how the approach presented in Sect. 4 can be used to optimize a MEMS resonator with a prescribed frequency t while minimizing its variance. Problem (22) is rewritten as Fig. 3 Problem setting, boundary conditions, and initial layout. In (a) the black region is the proof mass that remains fixed, whereas the gray areas represent the design space which is initially seeded with holes (b)  In this formulation, the expected natural frequency is imposed using the third constraint, whereas its variance is minimized through the objective function. We solve this problem for t = 1 MHz, A max = 70% , = 0 • C , and = 100 • C . The desired frequency corresponds to the oscillations of the proof mass along the y axis (Fig. 3).
The results are reported in Fig. 5. This robust solution is compared to a layout obtained through a deterministic approach, i.e. imposing the natural frequency without minimizing its variance. Both solutions exhibit the desired natural frequency of 1 MHz. However, the robust layout presents a much lower variance (4.09 kHz) with respect to the deterministic solution (68.59 kHz). Concerning the TCF value (Sect. 2), the deterministic and the robust designs are characterized respectively by −82.83 ppm∕ • C and −20.22 ppm∕ • C . In particular, the TCF of the robust design is one-third lower than the typical TCF of silicon resonators ( −30 ppm∕ • C). (46) All the layouts in Fig. 5 have been validated using COM-SOL Multiphysics © . The geometry is discretized using a body-fitted mesh of triangular elements. The plane stress approximation is used along with a geometrically linear formulation. After enforcing the boundary conditions, the natural frequencies are computed by solving an eigenfrequency step. By following this procedure, we have computed the target natural frequencies (Table 3) at reference temperature ( ΔT = 0 • C ). In addition, Fig. 6 shows the target mode shapes of the resonators (eigenmode corresponding to the oscillation of the proof mass along the y axis). The eigenfrequencies in presence of thermal effects are simulated by coupling the mechanical and thermal domains. Firstly, given a temperature variation ΔT , Young's modulus is modified according to Eq.(3). After that, a thermal expansion is applied in the entire domain. Then, a static step is used to compute the residual thermal stresses. Finally, the natural frequencies are computed by solving an eigenfrequency step starting from the results of the previous step. With this procedure, we have computed the TCF values (Table 4) related to the target natural frequencies of the resonators.
The discrepancy between all the numerical results (Table 3) is below 1%, thus confirming the quality of the solutions obtained with OpenLSTO.

Effects of the boundary conditions
The method proposed in Sect. 4 involves residual thermal stresses, which depend on the boundary conditions. In this Section, we use the same problem formulation presented in Sect. 5.2 (Problem (46)) to show how the optimal solution changes in presence of different boundary conditions. In particular, we consider the two boundary conditions illustrated in Fig. 7. Figure 8 shows the results obtained with the partial fixed boundary conditions (Fig. 7a), whereas Fig. 9 shows the results obtained with the four corners boundary conditions (Fig. 7b). In all the examples, the target frequency of 1 MHz has been reached. In addition, the robust solutions exhibit lower values of frequency variance and TCF compared to their deterministic counterparts, thus confirming the validity of the proposed approach.

Conclusion
In this work, a novel approach that deals with temperature uncertainties has been described. Statistical analysis was used to obtain analytical expressions for the expected natural frequencies and their variance over temperature. These values were used to formulate the temperaturerobust eigenfrequency optimization, which aims at maximizing the expected natural frequency while minimizing its variance due to temperature variations. The optimization problem was solved using a level-set optimization algorithm, and the sensitivities were computed through the adjoint method by exploiting the definition of Static Modal Derivatives (SMD). Even though there exist more accurate approaches for the calculation of modal derivatives, the static approximation was chosen mainly due to its simplicity and computational efficiency. A possible limitation of the proposed approach lies in the first-order approximation of the natural frequencies. While higher-order terms would improve the accuracy of the proposed method, their inclusion would increase the complexity of the formulation, especially that of the statistical and sensitivity analysis. Numerical examples considering different temperature conditions, weights, and boundary conditions were presented to support the effectiveness of the proposed approach, along with an example of a 1 MHz temperaturerobust MEMS resonator. The proposed layout, which has been validated using COMSOL, exhibits a TCF coefficient comparable to that of real-world devices. It is worth mentioning that the desired TCF value may not always be achievable. Nonetheless, any attempt aimed at minimizing the TCF will reduce the effort required by other active compensation techniques as well as the amount of energy required by such systems.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement. STMicroelectronics (Award Number 4000614871).

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.

Replication of results
OpenLSTO (Kambampati et al. 2018), freely available at https:// github. com/ M2DOL ab/ OpenL STO, was used to obtain the results above. The sensitivities equations, as well as the implementation details, have been provided in the previous sections. The linear systems and the eigenvalue problems were solved using the libraries Eigen (Guennebaud et al. 2010) and Spectra (Qiu 2015). The solution to the sub-optimization problems was computed using the open-source library IPOPT (Interior Point OPTimizer, Wächter and Biegler (2006)).

Ethical approval
The manuscript does not contain clinical studies or patient data.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 9
Comparison between the deterministic layout (a), obtained with a deterministic level-set optimization, and the robust layout (b), obtained by solving Problem (46). Solutions obtained with the four corners boundary conditions (Fig. 7b)