Two-Phase Approach for Fast Topology Optimization of Multi-Resonant MEMS Involving Model Order Reduction

TL;DR

This paper introduces a two-phase optimization method to speed up the design of multi-resonant MEMS devices using model order reduction.

Contribution

A novel two-phase optimization approach combining BESO and density-based methods with model order reduction for faster MEMS design.

Findings

The two-phase approach achieves optimal MEMS designs within 200 iterations.

Model order reduction reduces computation time by 50% in the examples tested.

Abstract

In this work, we propose a two-phase approach for a fast topology optimization of multi-resonant MEMSs. The approach minimizes the computation effort required to achieve an optimal design. In the first step, we perform a pre-optimization using bi-directional evolutionary structural optimization (BESO). We found in previous research that BESO can achieve optimal MEMS designs in a significantly lower number of iterations when compared to classical density-based methods. However, we encountered convergence issues with BESO towards the end of the optimization. Therefore, we introduced a second, density-based optimization phase to circumvent this issue. Finally, we introduced model order reduction to reduce the optimization time further. The novel approach is benchmarked with the design task of two common multi-resonant MEMS devices: a linear gyroscope and a micromirror. We show that the two-phase approach can achieve an optimal design within 200 iterations. With the addition of MOR, the computation of the goal function can be further reduced by 50% in our examples.

Figures

Comparison of the optimized geometries obtained by (a) two-phase optimization, (b) BESO, and (c) density-based TO with identical objective functions and constraints. The frequencies of the respective structures are compared in Table 3.

Full text

1. Introduction

Microelectromechanical systems (MEMSs) play an important part in our everyday life in different applications across a number of fields, e.g., telecommunications, consumer electronics, mobility, building automation, or healthcare. Their importance is expected to further increase at a growing rate due to the current development of the Internet of Things and portable/wearable electronics [1]. Many MEMS devices are resonant structures that are operated in resonance to realize their functions. These functions include sensing [2], actuation [3], timing [4], or filtering [5], with the very first application being reported in [6]. Since then, MEMS resonators have become prevalent in these aforementioned applications. This is mostly due to the possibility of fabricating these structures in batches with low-cost fabrication, similar to techniques used for integrated circuit manufacturing. Moreover, this similarity enables solutions for combining and integrating MEMSs with integrated circuits into a system on a chip. This leads to systems with superior form factors, robustness w.r.t. electromagnetic compatibility, and lower parasitic capacitances in the electrical connections between the involved components [7].

The design of MEMS resonators addresses their resonant modes and corresponding frequency characteristics. This includes tuning desired resonant modes to desired frequencies and suppressing undesired resonant modes. However, the devices’ resonance frequencies can often not be chosen independently, and geometric design alterations affect multiple resonance frequencies at once. In general, the field of MEMSs is an interdisciplinary field. MEMS devices intrinsically involve multiple physical domains, requiring different types of representation with varying levels of abstraction, making the design of MEMSs a highly complex matter. As reported in [8], the MEMS design process is still mostly performed in a “trial-and-error fashion”. It often relies on an expert’s experience with similar devices to find a suitable layout. This motivates a more systematic approach for identifying design geometries with desired properties, to increase the effectiveness and efficiency of the design process.

A systematic approach to geometry alteration is provided by mathematical optimization. It alters a design based on a user-defined objective function under given constraints. Typically, the layout of MEMSs is composed of simple geometrical substructures, such as rectangular mass blocks and flexible springs in the form of cantilevers or folded beams etc. Therefore, design optimization of MEMSs mainly consists of altering the geometric dimensions of these substructures. Due to the complex design space inherited from the complexity of the MEMS design process, heuristic approaches are preferred to traditional optimization algorithms [9]. Ref. [9] provides an overview of evolutionary design optimization approaches, which are predominant in the field of MEMSs. In [10,11], the authors propose an efficient, optimization-based framework for the design of MEMS gyroscopes and other inertial devices. The framework incorporates a combination of sub-structuring and static model reduction techniques. Recent studies on the optimization of scanning mirrors can be found in, e.g., Refs. [12,13], focusing on different aspects of the design, such as driving torque and crosstalk during actuation or dynamic deformation.

However, optimizations based on geometrical dimensions are strongly limited by the initial design, and the composition of simple geometric substructures will remain, although in optimized dimensions. More design freedom can be provided for topology optimization (TO). Instead of altering geometric dimensions, TO optimally distributes material onto a fixed design space. Thus, in theory, any shape or form can be achieved. This material distribution method was first introduced in [14]. In [15], the popular solid isotropic material with penalization (SIMP) method was introduced. The SIMP method assumes an isotropic material that is to be distributed over a finite-element (FE) domain. It introduces elemental pseudo-densities for each element relative to the density of the isotropic material as a design variable. This design variable determines the Young’s modulus and density via the respective power-law interpolation scheme. For SIMP-based TO, the method of moving asymptotes [16] or the optimality criteria (OC) method (e.g., [17]) is typically incorporated to find an optimum. An alternative approach to TO is evolutionary structural optimization (ESO) methods. The first examples of ESO can be found in [18,19] or [20]. ESO achieves optimal design by removing inefficient elements according to some predefined measure. The method was later extended to bi-directional ESO (BESO) in [21], which allows elements to be reintroduced. An overview of most state-of-the-art TO approaches is presented in [22,23], in which they are reviewed and compared. Some of the most recent studies have also introduced artificial intelligence (AI) to the field of TO. A comprehensive review of TO via machine and deep learning is provided in [24].

Even though TO was originally designed for compliance minimization (as shown in, e.g., [14]), it has proved effective in a plethora of applications across all fields of engineering, e.g., in aerospace engineering [25] or in the field of meta-materials [26]. The first relevant works, w.r.t. the field of MEMSs in the form of multi-frequency optimization, can be found in, e.g., [27,28,29] (SIMP) or [30] (BESO). A review discussing TO for vibration problems can be found in [31]. Early application of TO in actual MEMS devices can be found in [32,33]. More recent works have shown successful TO of many related devices: In [34,35], the authors consider three different optimization formulations for the optimization of MEMS gyroscopes to achieve targeted resonance frequencies; Ref. [36] proposes a TO to improve the dynamic stability of a resonant MEMS scanner used for LiDAR applications, and in [37], TO is performed on a multi-resonant piezoelectric energy harvester to maximize its power output for a predefined frequency range.

Because of the optimization problems posed by MEMS design, TO may require very high computational effort, especially when broadband frequency responses or transient behaviors (due to parasitic effects) and fine resolutions are required. This is why model order reduction (MOR) has been introduced to remedy this issue. MOR replaces the original high-fidelity model with a significantly lower-dimensional but still highly accurate surrogate. Thereby, it can reduce the computation time required for model evaluation by several magnitudes [38]. An up-to-date overview of classical MOR in TO is provided in [39]. Some recent studies also suggest using AI to reduce the computational burden of TO. Ref. [40] proposed a deep learning-based MOR that allows TO to be performed on coarser meshes while obtaining similar results compared to traditional methods on finer meshes.

In this work, we explore a two-phase TO approach for multi-resonant MEMS applications. The two phases combine an initial evolutionary TO with a subsequent density-based TO. Aside from its compatibility with commercial FE software [23], we found in our prior work that BESO can achieve an optimal design, requiring a significantly smaller number of iterations compared to density-based methods. However, in our previous work, we encountered convergence issues during BESO of resonant structures. This is why we introduce a second, density-based optimization phase to circumvent this issue. Such an approach has been discussed in the community; however, to our knowledge, it has yet to be implemented and evaluated. The only other work combining SIMP with BESO we are aware of is [41]. The aforementioned work is published in the field of multi-scale TO for lightweight cellular material, where BESO is applied to the optimization of the microstructure while SIMP is simultaneously applied to the microstructure.

The two-phase workflow is designed in such a way that the major optimization task is carried out during the BESO phase in a small number of iterations. Thereafter, a second optimization phase using density-based methods is performed to achieve algorithmic convergence. Since we aim to reduce the computational effort to a minimum, we also introduce MOR techniques as an additional measure to achieve this goal. Note that in our implementation, MOR is only considered in the density-based TO phase. This is because the commercial FE solver can provide all the required information during the BESO phase.

For the evaluation of the two-phase approach, we choose simple and well-established MEMS structures: an MEMS gyroscope and an MEMS scanning mirror. For the design of the MEMS gyroscope, we compare the results obtained by the combined approach with the results obtained by both BESO and the density-based approaches, individually, based on their performance given a limited number of iterations.

This paper is organized as follows. In Section 2, we briefly introduce the mathematical models of MEMS resonators and their reduced-order modeling. Section 3 defines the optimization problem we consider for the design of this class of devices to achieve the desired resonant frequency and amplitudes and avoid undesired ones. This includes the definition of an objective function, constraints, and sensitivities. In Section 4, the novel two-phase approach is proposed, before it is applied to the example models in Section 5. In Section 6, we conclude the paper and give an outlook on possible future research.

2. Numerical Modeling of Mechanical Resonators

In this section, we briefly introduce the numerical models involved in this paper. We first describe the model we obtain from the FE method. Subsequently, we provide a short outline of the projective model order reduction that is implemented in this work to speed up the optimization.

2.1. Finite Element Modeling

In this work, our geometries are modeled with industrial FE software (cf Section 5), using linear 3D hexahedron elements. The spatial discretization of the governing partial differential equation results in a second-order linear dynamical system, which is described by the following system of ordinary differential equations (see, e.g., [38,42]):

$$[eqn]$$

Here, in (1), $[eqn]$ is the state vector, i.e., a vector containing all nodal displacements, and $[eqn]$ are, respectively, the mass, damping, and stiffness matrices of the structure. As a simplification, we set $[eqn]$ , making the dynamics undamped. The remaining system matrices are assembled from respective element matrices $[eqn]$ and $[eqn]$ :

$$[eqn]$$

where the operator $[eqn]$ indicates the standard FE assembly process. u is the input of the system, which is the mechanical force in this case. $[eqn]$ is the input vector, which distributes the force upon the structure nodes. $[eqn]$ is the user-defined output of the system, gathered via the input matrix $[eqn]$ .

2.2. Model Order Reduction

There are many different ways to approach the model reduction of dynamical systems. In this work, we focus on the projective approaches. Their common idea is to achieve complexity reduction by approximating the state vector with a surrogate in a significantly lower dimensional basis and then projecting all governing equations onto a compatible low-dimensional subspace [43,44]. In detail, let the columns of $[eqn]$ be the basis of an r-dimensional subspace $[eqn]$ , with $[eqn]$ . An approximation of the state vector is given by $[eqn]$ , with $[eqn]$ . The residual $[eqn]$ of this approximation with reference to Equation (1) is given by

$$[eqn]$$

and the output equation is transformed into

$$[eqn]$$

By enforcing the Petrov–Galerkin condition,

$$[eqn]$$

we obtain the following reduced system:

$$[eqn]$$

where the reduced matrices are computed as

$$[eqn]$$

Different projective MOR approaches differ in their choice of projection matrices V and W. There are mainly three classes of approaches, approaches based on singular value decomposition (SVD), approaches based on Krylov subspaces, and approaches based on a combination of the first two. A detailed introduction of the approaches exceeds the scope of this work. Therefore, we will only give a brief overview of approaches widely used for the reduction of dynamical systems. A detailed introduction of MOR is provided in [43].

SVD-Based MOR

The class of SVD-based MOR methods includes widely used MOR methods such as proper orthogonal decomposition (POD), modal truncation, and truncated balanced realization (TBR). All these approaches utilize SVD to obtain their respective projection matrices—i.e., the singular vectors corresponding to the most dominant singular values are used as a projection basis. For more details, please refer to [43].

The most established SVD-based approach for the reduction of linear time-invariant dynamical systems is modal truncation. Instead of snapshot samples (POD) or Gramians (TBR), an eigenvalue decomposition is performed on the matrix pencil generated by the system matrices $[eqn]$ . The projection with the corresponding eigenbasis leads to a reduced-order system, which only retains the most dominant poles, i.e., the largest eigenvalues, of the system.

Note that when an r-dimensional modal truncation is applied, i.e., when the projection basis is chosen as $[eqn]$ , where $[eqn]$ is the i-th eigenvector of the system normalized by the mass matrix, i.e., $[eqn]$ , the reduced system matrices become

$$[eqn]$$

where $[eqn]$ gives the structure’s i-th resonance frequencies.

Note that there are two established normalization methods for eigenvectors, unity and mass normalization. Typically, eigenvectors are normalized by the mass matrix in industrial software. Therefore, mass normalization is considered throughout this work. However, we do not expect different normalization to affect the outcome of the optimization.

Due to the simple test cases, we implemented modal truncation in this contribution. This is only to establish a procedure and showcase the potential of MOR in combination with the proposed two-phase method. We expect state-of-the-art MOR methods to provide better performance when more complex optimization problems are considered, and especially when transient behaviors need to be considered due to parasitic effects.

3. Optimization Problem

The major design goal of multi-resonant MEMSs is to achieve predefined resonant frequencies for desired modes and suppress undesired or parasitic resonance effects. This is mostly achieved by adjusting the design in such a way that the resonance frequencies of parasitic modes appear at higher values. The secondary objective is mostly aimed at maximizing the amplitudes of the desired modes, e.g., in the case of a sensor, the sensor response, and for energy harvesting, the power output of the device. Therefore, we formulate the general optimization problem for this class of devices as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Here, we have defined the objective function (9a) as the maximization of static compliance of the structure, when a force is applied to the structure, e.g., actuation forces. The corresponding static force equilibrium is solved in (9d). We found that an energy-based objective function stabilizes the convergence of BESO. The optimization is then set up in such a way that the objective function becomes less influential compared to the frequency constraints, especially during the second optimization phase. Equations (9b)–(9f) define the constraints of the problem. (9b) enforces the frequencies of usable resonance modes $[eqn]$ to match the desired values $[eqn]$ , while (9c) imposes a lower bound for undesired modes, often referred to as parasitic modes. It is defined as a ratio $[eqn]$ between a desired mode i and an undesired mode j. (9e) is the eigenvalue problem that needs to be solved during optimization. Finally, (9f) restricts the volume of the structure to a predefined value $[eqn]$ and $[eqn]$ ranges from 0 to 1. For SIMP-based approaches, we enforce $[eqn]$ by introducing $[eqn]$ to avoid singularities.

In order to solve the optimization posed in (9), we use the well-established Lagrangian method. This extends the objective function (9a) to a Lagrangian function as follows:

$$[eqn]$$

with slack variables $[eqn]$ for the respective frequency constraints, turning them from inequality constraints into equality ones. $[eqn]$ denotes the Lagrange multiplier for the respective constraint.

The sensitivity of the Lagrangian function (10) with respect to the design variables $[eqn]$ is given by

$$[eqn]$$

where $[eqn]$ denotes the displacement of element e’s nodes corresponding to mode i or j’s respective mode shapes. The derivatives of the frequencies are computed using Rayleigh quotients:

$$[eqn]$$

where $[eqn]$ contains the component of the $[eqn]$ -th eigenvector corresponding to the nodes of element e. This leads to

$$[eqn]$$

To solve the optimization problem, the solutions of the additional variables of the Lagrangian function are also required. Therefore, they are treated like design variables. The derivatives w.r.t. the Lagrange multipliers are the respective constraints—i.e.,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

and the derivatives with respect to the slack variables result in

$$[eqn]$$

At an optimum, it is a necessary condition (Karush–Kuhn–Tucker) that all these derivatives (14)–(17) vanish. Therefore, with their help, an update scheme can be deployed for the search of these additional variables. In this work, we follow the update scheme suggested by [45,46]. As the domain of the additional variables introduced by the Lagrangian function is $[eqn]$ , ref. [46] introduces a function to reduce the search space. They suggest mapping each of the Lagrange multipliers from $[eqn]$ to $[eqn]$ :

$$[eqn]$$

making them more compatible with the range of $[eqn]$ . Note that with a searching scheme, the Slack variables $[eqn]$ only serve as indicators for the respective constraint, indicating whether it is active or inactive. Therefore, they do not need to be determined [46].

4. Topology Optimization

Topology optimization methods that are based on FE analysis mainly differ in their respective updating strategy. In [47], we applied BESO design multi-resonant MEMSs. Even though the optimizations resulted in promising designs, we experienced convergence issues with BESO. The resonant frequencies oscillate around the desired value during later iterations and sometimes even diverge (see Figure 1).

To remedy this issue, we suggest the introduction of a second optimization phase, in which a density-based approach is implemented to achieve algorithmic convergence. In this section, we illustrate the two TO methods applied in this work and propose a consecutive combination of both methods.

4.1. Bi-Directional Evolutionary Structural Optimization

The element updating rule of BESO deactivates elements with the lowest sensitivity values in each iteration. The deactivation threshold $[eqn]$ is defined in every iteration and is derived from the dynamic evolutionary rate ER^(i)^, which specifies the number of elements remaining after each design update:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the predefined upper and lower bounds for the evolutionary rate during optimization. $[eqn]$ and $[eqn]$ are parameters that control the location of the arctangent function’s inflection point and the function’s steepness around this point. Figure 2 shows the evolutionary rate ER versus volume fraction for some common $[eqn]$ and $[eqn]$ values in Figure 2. The update process is illustrated in Figure 3.

The sensitivity of an element w.r.t. compliance equals its strain energy [45], and the elemental sensitivity w.r.t. the resonance frequencies is derived in [48]. Coming from (11) and (13), the same sensitivities can be obtained by replacing the derivatives of the elemental matrices by

$$[eqn]$$

as elemental matrices vanish upon deactivation. The volume constraint is naturally considered in the updating scheme, so the corresponding part can be omitted, without affecting the order of the ranking:

$$[eqn]$$

To equal out the differences of the partial sensitivities induced by, e.g., the choice of units, we normalize the sensitivity values and project them onto the interval $[eqn]$ (similar to [45]):

$$[eqn]$$

where $[eqn]$ refers to the part of the sensitivity corresponding to the i- or j-th constraint. The superscript max or min marks the maximum or minimum value of the partial sensitivity over all elements e. For example, let the objective function be maximizing the compliance of the first mode shape. In addition, the second resonance frequency should be optimized towards $[eqn]$ , and the third resonance frequency should be higher than $[eqn]$ . Then, according to (21) and (11), the sensitivity is given by the following:

$$[eqn]$$

with partial sensitivities $[eqn]$ to $[eqn]$ . The normalization is applied to each of the partial sensitivities such that

$$[eqn]$$

It prevents partial sensitivities from being over-dominant and makes the overall sensitivity ill-conditioned.

Subsequently, we apply a filter to the sensitivities [49,50] to ensure the existence of solutions to the TO problem and to prevent commonly experienced issues, e.g., “checkerboard patterns” (e.g., [51,52]):

$$[eqn]$$

where $[eqn]$ is the set of adjacent elements in the neighborhood of element k and $[eqn]$ is the weight function based on a user-specified filter radius $[eqn]$ and the Euclidean distance between the elements k and l ( $[eqn]$ and $[eqn]$ denote the coordinates of the center point of the respective element). This sensitivity filtering is effectively a blurring of the values using a hat function. Therefore, this operation enables sensitivity values to be assigned to an inactive element, allowing it to become active again.

The reactivation of elements is further encouraged by the final modification. In [53], possible convergence issues were reported for evolutionary TO, such as jumps or the oscillation of objective values. One reason for these issues is the inaccurate sensitivities of the inactive elements, which may not be present in the FE-Model. Therefore, the authors suggested relaxing the difficulty of convergence by introducing history-averaging that smoothens the change sensitivity values from the $[eqn]$ -th to the $[eqn]$ -th TO iteration:

$$[eqn]$$

A major challenge during BESO is to conserve the structural integrity of the layout, i.e., to ensure that the layout does not become disconnected. Many studies have been published in recent years tackling this issue, e.g., [54]. The current implementation uses a similar approach to the aforementioned publication. Instead of a connectivity matrix, we use image recognition to check for the number of bodies present in each iteration. In case there is more than one body present during the check, a test identifies the elements causing the loss of connection. These elements (and their neighbors) are then excluded from optimization. This approach can also be used to control the minimum length scale present in the structure if a minimum number of adjacent elements is required for each element.

4.2. Density-Based Approach

For the second phase, the traditional simple isotropic material with penalization (SIMP) is implemented for the interpolation of material properties:

$$[eqn]$$

$[eqn]$ and $[eqn]$ denote the Young’s modulus and density of the material and $[eqn]$ and $[eqn]$ are the respective values assigned to void elements, as $[eqn]$ would lead to singular stiffness and mass matrices. p is the penalty factor, which is usually set to three by default (see [55] for a more in-depth discourse on the SIMP approach and the choice of the penalty factor).

For the density-based approach, we adopt the projection scheme from [56]. Here, the blurring filter is applied to the design variable instead of the sensitivity:

$$[eqn]$$

This allows us to introduce a subsequent threshold projection using a smoothed Heaviside function [56] to suppress gray elements, i.e., elements with intermediate density $[eqn]$ :

$$[eqn]$$

where $[eqn]$ is the threshold, e.g., an element with $[eqn]$ will be pushed towards being a void element and vice versa. The “aggressiveness” of the filter is controlled by the sharpness parameter $[eqn]$ , which determines the steepness of the slope in the Heaviside function. The Heaviside function for different values of $[eqn]$ and $[eqn]$ is plotted in Figure 4.

For the sake of consistency, the Lagrangian function (10) is also considered for the density-based TO. Therefore, the initial sensitivity is identical to (11), with

$$[eqn]$$

due to the SIMP interpolation scheme. This, then, has to be modified with the terms introduced by the applied filters and the chain rule:

$$[eqn]$$ $$[eqn]$$

where (31) corresponds to the blurring filter and (32) corresponds to the threshold projection. In total, the sensitivity of an element e is determined by

$$[eqn]$$

Furthermore, we also adopted the normalization introduced in (22). However, as opposed to BESO, not only does the relative ranking of the elemental sensitivities matter, but also their respective signs are required for density-based TO. Therefore, the normalization is adjusted to

$$[eqn]$$

with $[eqn]$ again being the part of the sensitivity corresponding to the i- or j-th constraint. $[eqn]$ is the infinity norm of the vector $[eqn]$ , which contains all $[eqn]$ of all elements.

4.3. Combining BESO with SIMP and MOR

The combined approach is defined as follows: at the beginning, BESO is initiated with a full design space. The discrete optimization is executed until the target volume fraction is achieved. Subsequently, we switch to the density-based approach. In this phase, we also introduce projection-based MOR to minimize the computational effort required during iterations. As the resonance frequencies (and therefore, the eigenvectors) need to be computed in every iteration, one could opt for a modal reduction of the model without any additional efforts required. Due to the pre-optimization with BESO, we do not expect large changes in structure; therefore, the projection basis V can be reused over a large number of iterations. This would make, e.g., Krylov subspace-based MOR or balanced truncation feasible for static objective functions.

As suggested in Section 4.1 and Section 4.2, the Lagrangian function is kept throughout both optimization phases for consistency. The Lagrange multipliers are determined with (18), where the supplement variables $[eqn]$ are initialized with small values and subsequently increased throughout both phases as long as the corresponding constraint is violated. A flowchart of the proposed TO approach is shown in Figure 5.

Note that we do not impose any geometric or gray penalization constraints as suggested in [34]. In our combined approach, the related issues are naturally taken care of by the BESO implementation. In our test cases, thin structures are exclusively generated during BESO optimization. Therefore, no “gray connections” (see [34]) need to be prevented. A threshold projection is still required to suppress gray areas. However, due to the BESO pre-optimization, the filter can be initialized at a higher $[eqn]$ value.

5. Optimization Results

In this section, we apply the proposed combined TO approach to two resonant MEMS structures. Both design spaces are implemented in Ansys^®^ Mechanical Enterprise Academic Research, Release 2023.1 [57]. TO is performed using pyAnsys, Release 2023.1.3 [58]. Both models are discretized using linear elements with an element size of two (μm). A detailed overview of optimization parameters used for each benchmark will be presented in the respective subsection. The optimal layouts are presented and discussed. Additionally, a comparison of the two-phase approach to BESO-only and purely density-based TO, as well as with the reference source, is provided for the first benchmark.

5.1. Benchmark: MEMS Gyroscope

In the first benchmark, we apply our combined approach to a simplified academic test case proposed in [34]. The considered structure is an in-plane linear MEMS gyroscope. Gyroscopes are inertial sensors for angular acceleration or velocities, making use of the Coriolis effect. For more details on MEMS gyroscopes, please refer to, e.g., [59]. Adopted from the reference paper, we target a drive frequency of $[eqn]$ , a sensing frequency of $[eqn]$ , and a frequency for the parasitic mode beyond $[eqn]$ . Instead of formally maximizing the sense response, we choose to maximize the compliance of the structure for the load case depicted in Figure 6. Low-frequency motions of MEMS structures with a concentrated proof mass and light suspending structures are dominated by the motion of the proof mass (a similar assumption is applied in [34]. Therefore, it is sufficient to impose a single force onto the mass’ center of gravity. In this way, we aim to maximize the deformation of the structure caused by the Coriolis force, which serves as a simplified estimation for the sense response. This is a deliberate choice, as compliance-based objective functions are applicable for a wide range of MEMS devices, whereas the sense response due to Coriolis force is very specific to gyroscope design.

The two-dimensional design space of this benchmark is depicted in Figure 7. The design space and the material parameters are both taken from the reference paper. The square proof mass has a dimension of $[eqn]$ . The design space on both sides of the proof mass is sized at $[eqn]$ . As described in the introduction, the element size is set to two. The element size of the model is increased compared to the reference paper due to instability in the communication between Python 3.9 and Ansys. To avoid undesired interaction between boundaries and filters, boundary blocks and a buffer zone are introduced [60]. A unitary isotropic material is assumed with Young’s modulus $[eqn]$ , Poisson’s ratio $[eqn]$ , and a density of $[eqn]$ (see Table 1).

For BESO, we set the upper and lower limits of the evolutionary rate to $[eqn]$ and $[eqn]$ . The parameters $[eqn]$ and $[eqn]$ are chosen as 40 and −20, respectively. The desired final volume fraction of the structure is set to $[eqn]$ , adopted from the reference paper. The Lagrange multipliers are initialized with $[eqn]$ . Upon violation of the constraint i or j, $[eqn]$ is increased to one—i.e., the corresponding auxiliary parameters $[eqn]$ are set to 0.5 and then increased by 0.004 after each iteration, as long as the constraint remains violated with a tolerance of $[eqn]$ . For filter parameters, the filter radius of the blurring filter used in this experiment is $[eqn]$ (reference: $[eqn]$ ). This is a consequence of the increased element size. In the threshold projection, we set $[eqn]$ and $[eqn]$ is initialized with $[eqn]$ , which is then increased linearly during the second phase. The increase is set to 2 every 20 iterations up to $[eqn]$ .

To reduce computation time, we exploit the symmetry of the design space and only perform TO on one of the quadrants (here, we use the bottom left quadrant). Additionally, a modal projection base, composed of the first three eigenvectors of the dynamical system is used for model order reduction during the SIMP phase. This means that only three linear equations need to be solved during each optimization loop (see Section 5.3). Finally, the modified optimality criteria method proposed in [61] is implemented as an optimizer. The change tolerance for convergence is set to 1% and the maximum number of iterations performed with the density-based method is restricted to 120.

The optimized layout is presented in Figure 8. A comparison of the achieved resonance frequencies versus the target values is listed in Table 2. Figure 9 shows the mode shapes of the optimized structure.

Comparison of Algorithms

In this section, we compare the results of the combined method with the results of BESO and density-based TO. For this comparison, we set up all optimizations with identical objective functions and constraints. The settings for both algorithms are identical to those deployed in the two-phase algorithm if possible. For density-based TO, we replaced the BESO pre-optimization by initially setting the filter parameter to $[eqn]$ . $[eqn]$ is then doubled every 20 iterations until $[eqn]$ . After that, $[eqn]$ is increased linearly, identical to the two-phase approach. With this configuration, the maximum number of iterations is identical across all optimizations. Additionally, we introducd the gray penalization constraints suggested by [34],

$$[eqn]$$ $$[eqn]$$

to ensure the integrity of the structure. Here, $[eqn]$ and $[eqn]$ are the drive and sense frequencies of the structure derived from the same set of intermediate design variables (before filtering) but with a more strict Heaviside filter, i.e., $[eqn]$ in our case. $[eqn]$ is adopted for this comparison.

The BESO-only approach resulted in a layout with some thin (folded) flexural beam structures, which closely resemble those utilized in classical MEMS layouts (see Figure 10b). Unfortunately, most of these structures disappear during the second phase, when combining both methods. The optimized geometry obtained by a purely density-based method resembles the structure from [34], which was also obtained by a purely density-based method. For comparison, the layouts in [34] are achieved with numbers of iterations ranging from 1232 to 1703. We capped the total number of iterations at 180 for all optimizations. This is due to the aforementioned issues concerning communication instability when interfacing with the commercial software. Consequently, the imposed tolerances for the target frequencies were increased from $[eqn]$ to 0.025 and the filter parameters also increased more rapidly. We do not claim and never expected the proposed, generic approach to outperform the reference paper, even with a similar number of iterations. However, this comparison with a well-tuned TO, specifically designed for this kind of device, clearly shows the potential and deficit of the approach.

Comparison of the optimized geometries obtained by (a) two-phase optimization, (b) BESO, and (c) density-based TO with identical objective functions and constraints. The frequencies of the respective structures are compared in Table 3.

More specifically, the comparison shows that BESO exhibits faster convergence towards optimal designs. We mainly attribute this to BESO not requiring complex filtering techniques. However, as stated in [62], BESO tends to converge towards inefficient local optimal designs, which is also apparent in this comparison. On a positive note, we can report that the inclusion of a second optimization phase successfully eliminated the convergence issues we encountered during BESO-only optimization. However, further investigation is required into whether a higher number of iterations with the density-based methods would lead to a recovery from the local optimum achieved by BESO for more efficient designs.

5.2. Benchmark: MEMS Scanning Mirror

In this benchmark, we apply our combined approach to a second, more realistic MEMS device, with actual dimensions and material properties. The device is an MEMS scanning mirror, often also referred to as a micromirror, with two rotational degrees of freedom. The tilting motion of the mirror, as shown in Figure 11, can be utilized to control the direction of reflected light rays. This feature is utilized in laser scanner projection applications. The typical performance parameters of an MEMS scanning mirror are its maximum scan angle, resonant frequency, resolution, and surface quality, i.e., smoothness and flatness [63]. In this benchmark, our approach is used to design a micromirror for image projection at 60 $[eqn]$ and a scanning resolution of 1080 horizontal lines. According to the desired performance parameters, we set the objective function to maximize the amplitude of both resonance modes to maximize their respective deflection angle. The objective function is set as the sum of compliance over the static load case depicted in Figure 12. In this way, we aim to maximize the horizontal deflection angle during operation. Note that a resonance frequency of 60 $[eqn]$ is not feasible for this considered design space due to its dimensions and the involved material. As a consequence, the vertical scanning motion will be actuated outside of resonance. Therefore, we minimize the frequency, which minimizes the energy required for actuation. In contrast, the horizontal scanning motion is driven in resonance; thus, we set the target frequency for the horizontal scanning mode to $[eqn]$ = 64,800 Hz.

The design space considered for this optimization is depicted in Figure 13. The mirror in the center of the structure has a dimension of $[eqn]$ similar to the structure proposed in [64]. The flexible structure surrounding the mirror is defined as the design space and has a dimension of $[eqn]$ . Also here, boundary blocks and buffer areas are introduced to avoid undesired boundary interactions. The material parameters are Young’s modulus $[eqn]$ , Poisson’s ratio $[eqn]$ , and density $[eqn]$ (see Table 4).

For the BESO, we reduced the limits of the evolutionary rate to $[eqn]$ and $[eqn]$ . $[eqn]$ and $[eqn]$ are kept from the gyroscope benchmark. The final volume fraction is set to 0.5. Again, the Lagrange multipliers are initialized with $[eqn]$ , and upon constraint violation, they and their auxiliary parameters $[eqn]$ are increased to 1 and 0.5, respectively. Here, the subsequent increase is changed to 0.01 after each iteration until the constraint becomes inactive within a tolerance of $[eqn]$ . All filter parameters remain identical to those used for the first benchmark, i.e., $[eqn]$ , $[eqn]$ , and $[eqn]$ . The same applies to the optimizer, maximum number of iterations, and convergence tolerance. Finally, quarter symmetry is again assumed, and the same “three-modes” modal reduction is performed to speed up the SIMP phase (see Section 5.3).

The optimized layout is presented in Figure 14. Table 5 shows the resonance frequencies of the optimized layout and the corresponding objective for the respective modes, and Figure 15 depicts the mode shapes of the optimized structure.

The optimized geometry achieved a layout satisfying the objectives and constraints defined in the optimization problem. To minimize the first resonance frequency, the material is concentrated at the top and bottom of the layout, and thin flexural beams are formed extending from the fixed support. These resemble features present in published designs, e.g., [65]. However, it is also evident that, if feasible, increasing the size of the design space would greatly benefit the optimization for the given objectives. Additional mass and an increasing length of the thin flexural structure would allow the optimizer to further reduce the resonance frequency of vertical scanning. Alternatively, an increase in thickness may also benefit the frequency characteristics. It would also allow for a connection between the mirror and design space with a reduced width, without compromising on the resonance frequency for horizontal scanning. Reducing the width of the connection would improve the planarity of the mirror during operation.

5.3. Model Order Reduction

As stated in Section 5.1 and Section 5.2, model order reduction is introduced to alleviate the computational burden. The model is reduced in such a way that the reduced-order model (ROM) only contains four degrees of freedom (DoFs). This is considerably lower in dimension compared to the full-order model (FOM) originating from FE modeling. The computation times required for the TO of both benchmark devices with and without the inclusion of MOR are shown in Table 6. For this comparison, we only consider the net computation time required for the second phase of the optimization process. This excludes, i.e., all file input and output and the overhead generated by the interaction between Python and Ansys. Table 6 shows a significant reduction in computation time when MOR is involved. The computation time required for the optimization of the gyroscope is reduced by half. The reduction in computation time achieved during the optimization of the micromirror is 66%. This demonstrates the potential of MOR-aided TO.

6. Conclusions and Outlook

In this work, we investigated the feasibility of a two-phase approach for a fast topology optimization of multi-resonant MEMS applications. The approach includes an initial optimization using BESO and a subsequent SIMP-based TO to achieve convergence. To further reduce the computation time, we introduced modal-based MOR during the second phase.

The proposed approach was applied to the design of two different benchmark models. In the first benchmark, we applied the approach to a single-mass, in-plane MEMS gyroscope. The optimized layout presented in Section 5.1 achieved the desired target resonance frequencies within the given tolerances. The comparison with BESO-only and purely density-based TO shows the following:

• 1.The convergence issues of BESO can be circumvented with the inclusion of a subsequent, density-based optimization;
• 2.BESO requires fewer iterations to converge to optimal designs …;
• 3.However, BESO tends towards local, inefficient optima that may not be recovered with a low number of density-based iterations.

In a second benchmark, we applied our two-phase approach to the design of a micromirror with two rotational degrees of freedom. The optimized layout was designed for imaging projection, with a resolution of $[eqn]$ at 60 $[eqn]$ . The proposed approach converged to a layout fulfilling the posed requirements. Depending on the spatial availability in the given application, a dimensional increase of the design space may further improve the design. Note that both models are generic test models. For specific applications, there are typically additional performance indicators or constraints to be considered, e.g., rigidity for specific disturbances or fabrication constraints. These, however, can be easily considered during optimization through additional terms in the objective function or as additional constraints.

The inclusion of MOR significantly reduced the computation time required for the second phase. We expect the advantage of MOR to further increase when the goal function becomes more complex. This is because MOR only affects the objective function. However, the optimization also includes many processes that do not benefit from MOR, such as the computation of the mode shapes and the resonance frequencies. Therefore, an optimized implementation of the method will also increase the advantage of MOR, as this would reduce the overall computation time, particularly the time required by the processes not affected by MOR.

During our experiments, we did not experience any change in convergence or final optimized layout due to MOR for both of the presented benchmarks. However, this needs further investigation. Future research may also exploit AI-based prediction for objective values such as resonance frequencies or compliance to achieve even faster optimizations: In [66], a neural network is trained to predict the resonant frequency, thermoelastic quality factor, and other quality measures of an MEMS resonator that can be relevant for design optimization. The prediction is significantly faster than traditional numerical simulation while delivering similarly accurate results.

Another relevant issue that shall be addressed in future research is the preservation of the structural integrity of the layout. The current implementation requires significant computational effort. The restore and exclude mechanism designed to prevent the layout from becoming disconnected makes the BESO process very sensitive to the change of optimization parameters, and the final layout is dependent on the course or history of the optimization. A more sophisticated approach could be, e.g., describing the FE mesh as a weighted graph, with elemental sensitivities as weights. Then, one could implement a path-finding algorithm such as A* [67] to find the shortest path along elements w.r.t. these sensitivities in each iteration. This would be a more deterministic approach and ensure structural integrity at minimal cost w.r.t. the optimization problem.