1 Introduction↩︎

Nonlinear dynamic responses are commonly observed in micro-electro-mechanical systems (MEMS). On the one hand, such responses may cause the systems to fail to work properly. For example, in some high-precision sensors, such as gyroscopes [1], [2], pressure sensors [3], nonlinear responses reduce the accuracy of the sensor and even cause malfunctions [4]. On the other hand, nonlinearities can also be utilized to provide advantages over linear designs. For example, nonlinear dynamic responses can be used to improve the performance of energy harvesters [5]–[7], and their characteristic hardening or softening behavior can be utilized in the design of filters [8] and binary counters [9]. Thus, it is essential to tailor the dynamic responses of nonlinear systems if one aims to suppress the amplitude of structural vibration or to utilize them to optimally design devices.

To utilize nonlinear dynamic responses, several studies have focused on tailoring the hardening and softening behavior of nonlinear systems. One approach involves adjusting the hardening behavior of clamped-clamped beams by maximizing or minimizing the coefficient of cubic nonlinearity in the normal form [10]. Another method optimizes both the Q-factor and Duffing nonlinearity in finite element-based reduced-order models using derivative-free optimization techniques [11]. Additionally, backbone curves have been optimized on a reduced-order model (ROM) constructed via Lyapunov subcenter manifolds [12], [13] and spectral submanifolds (SSMs) [14].

Tailoring the hardening and softening behavior of a nonlinear system typically involves analyzing its autonomous dynamics and does not require consideration of external forcing. However, if the objective is to suppress or enhance the vibration amplitude under harmonic excitation, it is necessary to tailor the forced response curve (FRC). For low-dimensional nonlinear systems, the FRC of nonlinear systems can be computed by the harmonic balance method [15]–[17] combined with continuation techniques [18], [19]. Based on these techniques, [20] applied shape optimization to a clamped-clamped beam to minimize or maximize the amplitude of primary resonance. For high-dimensional nonlinear systems, such as finite element models (FEM) in topology optimization [21], computing and tailoring the FRC becomes challenging due to the increased computational cost.

In addition to amplitude control, tailoring the bifurcations of the FRC is also essential to nonlinear systems. Under harmonic excitation, nonlinear systems often exhibit saddle-node (SN) bifurcations on the FRC as the forcing amplitude increases. These bifurcations lead to abrupt jumps in the steady-state response, known as jump or hysteresis phenomena, which are particularly detrimental to devices intended to operate within a linear or predictable regime. Therefore, suppressing the occurrence of SN bifurcations becomes critical. In previous work, shape optimization along with numerical continuation is used to tailor the distance of two SN bifurcations on the nonlinear FRC [22].

Here, we aim to use topology optimization to tune the FRC of high-dimensional nonlinear systems, including minimizing or maximizing the peak of FRC, manipulating the hardening/softening behavior, and suppressing the occurrence of SN bifurcations. As a method to design continuum structures [23], topology optimization can produce more refined structures than shape optimization. However, the number of degrees of freedom (DOFs) of the FEM and the associated design variables in topology optimization is considerable.

To address the challenge of high-dimensionality, reduction on SSMs has been used to analytically extract the FRC of FEM with high DOFs [24]–[27]. In particular, SSM reduction transforms forced periodic orbits into fixed points of the associated reduced-order model (ROMs). In addition, the sensitivity of coefficients in SSM-based ROMs that govern the peak of FRC and hardening/softening behavior has been discussed in [28]. These sensitivities and their efficient computation [12] are essential in the iteration process of topology optimization. Based on these contributions, we will perform topology optimization of FRCs.

In this work, we adopt a density-based topology optimization algorithm to determine the optimal material layout. Specifically, the Solid Isotropic Material with Penalization (SIMP) scheme [29] is used to interpolate material properties and promote 0–1 designs. To enhance numerical stability and obtain clear structural boundaries, both density filtering and density projection techniques are employed [30]. The design variables are updated using the Method of Moving Asymptotes (MMA) [31]. All these techniques follow the implementation described in [12].

The rest of this paper is organized as follows. In Sec. 2, we introduce the SSM theory and the problem description about minimizing/maximizing the peak of FRC, tailoring the hardening/softening behavior, and suppressing the occurrence of SN bifurcations. In Sec. 3, we derive the sensitivity of coefficients used in the iteration process of topology optimization. Three numerical examples are presented in Sec. 5. Finally, the conclusions are presented in Sec. 6.

2 Problem formulation↩︎

2.1 Setup↩︎

We consider a periodically forced nonlinear mechanical system \[\label{eq:eom-second-full} \boldsymbol{M}(\boldsymbol{\mu})\ddot{\boldsymbol{x}}+\boldsymbol{C}(\boldsymbol{\mu})\dot{\boldsymbol{x}}+\boldsymbol{K}(\boldsymbol{\mu})\boldsymbol{x}+\boldsymbol{f}_2(\boldsymbol{\mu},\boldsymbol{x},\boldsymbol{x})+\boldsymbol{f}_3(\boldsymbol{\mu},\boldsymbol{x},\boldsymbol{x},\boldsymbol{x})=\epsilon \boldsymbol{f}^{\mathrm{ext}}(\Omega t),\quad 0<\epsilon\ll1\tag{1}\] where \(\boldsymbol{x}\in\mathbb{R}^n\) is the displacement vector; \(\boldsymbol{\mu}\in\mathbb{R}^q\) is the design variable representing the spatial distribution of material density in the structure; \(\boldsymbol{M},\boldsymbol{C},\boldsymbol{K}\in\mathbb{R}^{n\times n}\) are the mass, damping and stiffness matrices; \(\boldsymbol{f}_2(\boldsymbol{x},\boldsymbol{x})\) and \(\boldsymbol{f}_3(\boldsymbol{x},\boldsymbol{x},\boldsymbol{x})\) are internal nonlinear elastic forces arising from geometric nonlinearities; and \(\epsilon \boldsymbol{f}^{\mathrm{ext}}(\Omega t)\) denotes external periodic forcing with small amplitude \(\epsilon\ll1\).

Solving the linear part of 1 , we have the following generalized eigenvalue problem \[\label{eq:eom-second-eig} (\lambda_j^2\boldsymbol{M}+\lambda_j\boldsymbol{C}+\boldsymbol{K})\boldsymbol{\phi}_j=\boldsymbol{0},\quad (\lambda_j^2\boldsymbol{M}+\lambda_j\boldsymbol{C}+\boldsymbol{K})^{\mathrm{T}}\boldsymbol{\psi}_j=\boldsymbol{0}\tag{2}\] where \(\lambda_j\) is an eigenvalue, and \(\boldsymbol{\phi}_j\) and \(\boldsymbol{\psi}_j\) are associated right and left eigenvectors. In this work, we do not allow for internal resonance \((\lambda_i \approx k\lambda_j,\; k \in \mathbb{Z}^+)\) such that the high-dimensional system 1 can be reduced on a two-dimensional spectral submanifold (SSM). When internal resonance occurs during the optimization iteration process, additional constraints are introduced to eliminate the internal resonance. This will be discussed in Sec. 2.3.1.

In this work, we consider positive definite mass and stiffness matrices as well as Rayleigh damping, namely, \[\label{eq:Rayleigh32damping} \boldsymbol{K} \boldsymbol{\phi}_j = \omega_j^2 \boldsymbol{M} \boldsymbol{\phi}_j, \;\boldsymbol{\phi}_j^\mathrm{T}\boldsymbol{M}\boldsymbol{\phi}_j = 1, \; \boldsymbol{C} =\alpha\boldsymbol{M} + \beta \boldsymbol{K},\tag{3}\] where \(\omega_j\) is natural frequency, \(\alpha\) and \(\beta\) are the stiffness-proportional and mass-proportional damping constants. The values of \(\alpha\) and \(\beta\) are computed by \(\boldsymbol{\phi}_j^\mathrm{T}\boldsymbol{C}\boldsymbol{\phi}_j =\alpha + \beta \omega_j^2 = 2 \xi \omega_j\), where \(\xi\) is the damping ratio of the system. Further, taking \(j=1,2\), yields \[\label{eq:Rayleigh32damping32coefficient} \alpha = \xi \frac{2\omega_1 \omega_2}{\omega_1 + \omega_2}, \quad \beta = \xi \frac{2}{\omega_1 + \omega_2}.\tag{4}\]

2.2 Reduction on spectral submanifolds↩︎

Computing and tuning the FRC of the mechanical system 1 is challenging for \(n \gg 1\). To address this challenge, we briefly introduce SSMs reduction theory in this subsection, which transforms the full system 1 into an SSM-based ROM.

The second-order system 1 can be transformed into a first-order system as below \[\label{eq:full-first} \boldsymbol{B}\dot{\boldsymbol{z}} =\boldsymbol{A}\boldsymbol{z}+\boldsymbol{F}(\boldsymbol{z})+\epsilon\boldsymbol{F}^{\mathrm{ext}}({\Omega t}),\tag{5}\] where \[\begin{gather} \label{eq:zABF} \boldsymbol{z}=\begin{pmatrix}\boldsymbol{x}\\\dot{\boldsymbol{x}}\end{pmatrix},\,\, \boldsymbol{A}=\begin{pmatrix}-\boldsymbol{K} & \boldsymbol{0}\\\boldsymbol{0} & \boldsymbol{M}\end{pmatrix},\,\, \boldsymbol{B}=\begin{pmatrix}\boldsymbol{C} & \boldsymbol{M}\\\boldsymbol{M} & \boldsymbol{0}\end{pmatrix},\nonumber\\ \boldsymbol{F}(\boldsymbol{z})=\begin{pmatrix}{-\boldsymbol{f}_2(\boldsymbol{x},\boldsymbol{x})-\boldsymbol{f}_3(\boldsymbol{x},\boldsymbol{x},\boldsymbol{x})}\\\boldsymbol{0}\end{pmatrix},\,\, \boldsymbol{F}^{\mathrm{ext}}(\Omega t) = \begin{pmatrix}\boldsymbol{f}^{\mathrm{ext}}(\Omega t)\\\boldsymbol{0}\end{pmatrix}. \end{gather}\tag{6}\] The associated generalized eigenvalue problem becomes \[\boldsymbol{A}\boldsymbol{v}_j=\lambda_j\boldsymbol{B}\boldsymbol{v}_j,\quad \boldsymbol{u}_j^\ast \boldsymbol{A}=\lambda_j \boldsymbol{u}_j^\ast \boldsymbol{B},\] where \(\lambda_j\) is a generalized eigenvalue and \(\boldsymbol{v}_j\) and \(\boldsymbol{u}_j\) are the corresponding right and left eigenvectors, respectively. In particular, we have \[\boldsymbol{v}_j=\begin{pmatrix}\boldsymbol{\phi}_j \\ \lambda_j\boldsymbol{\phi}_j\end{pmatrix},\quad \boldsymbol{u}_j=\begin{pmatrix}\bar{\boldsymbol{\psi}}_j \\ \bar{\lambda}_j \bar{\boldsymbol{\psi}}_j\end{pmatrix}.\]

In this study, we construct the two-dimensional SSM associated with the master subspace \(\mathcal{E}=\mathrm{span}(\boldsymbol{v}_1,\bar{\boldsymbol{v}}_1)\) for model reduction. Specifically, let \(\boldsymbol{p} = (p, \bar{p})\) denote the reduced coordinates. We seek the autonomous SSM map \(\boldsymbol{z} = \boldsymbol{W}(\boldsymbol{p})\) and its associated reduced dynamics \(\dot{\boldsymbol{p}} = \boldsymbol{R}(\boldsymbol{p})\) at \(\epsilon = 0\). When external periodic forcing \(\epsilon \boldsymbol{f}^{\mathrm{ext}}(\Omega t)\) is considered, the SSM becomes time-dependent, yielding a periodic SSM map \(\boldsymbol{W}_\epsilon(\boldsymbol{p},\Omega t)\) and reduced dynamics \(\dot{\boldsymbol{p}}=\boldsymbol{R}_\epsilon(\boldsymbol{p},\Omega t)\). The map and reduced dynamics on the periodic SSM can be expressed as [25], [28] \[\begin{align} \label{eq:red-reps} \boldsymbol{W}_\epsilon(\boldsymbol{p},\Omega t)&=\boldsymbol{W}(\boldsymbol{p})+\epsilon\boldsymbol{X}(\boldsymbol{p},\Omega t)+\mathcal{O}(\epsilon^2), \\ \dot{\boldsymbol{p}}=\boldsymbol{R}_\epsilon(\boldsymbol{p},\Omega t)&=\boldsymbol{R}(\boldsymbol{p})+\epsilon\boldsymbol{S}(\boldsymbol{p},\Omega t)+\mathcal{O}(\epsilon^2). \end{align}\tag{7}\]

Plugging 7 into 5 , we obtain the invariant equation as \[\begin{align} \boldsymbol{B}\partial_{\boldsymbol{p}} \boldsymbol{W}_\epsilon \boldsymbol{R}_\epsilon + \partial_{\Omega t} \boldsymbol{W}_\epsilon \Omega = A \boldsymbol{W}_\epsilon + \boldsymbol{F}(\boldsymbol{W}_\epsilon) + \epsilon \, \boldsymbol{F}^{\mathrm{ext}} \cos \Omega t. \label{eq:invar-balance} \end{align}\tag{8}\] The autonomous SSM map \(\boldsymbol{z} = \boldsymbol{W}(\boldsymbol{p})\) and its associated reduced dynamics \(\dot{\boldsymbol{p}} = \boldsymbol{R}(\boldsymbol{p})\) can be obtained by solving the invariant equation 8 at \(\epsilon = 0\), using a Taylor expansion in \(\boldsymbol{p}\) truncated at cubic order, namely \[\begin{align} \boldsymbol{W}(\boldsymbol{p}) = \sum_{1 \leq m+n \leq 3} \boldsymbol{W}_{mn}p^m \bar{p}^n, \\ \boldsymbol{R}(\boldsymbol{p}) = \sum_{1 \leq m+n \leq 3} \boldsymbol{R}_{mn}p^m \bar{p}^n. \label{eq:Taylor32Expansion32W32and32R} \end{align}\tag{9}\] Cohomological equations can be derived, and the expansion coefficients can be obtained following a normal-form-style parameterization. More details can be found at [28].

The non-autonomous part of the periodic SSM mapping \(\boldsymbol{X}(\boldsymbol{p}, \Omega t)\) can be approximated by the leading-order terms as \[\begin{align} \boldsymbol{X}(\boldsymbol{p},\Omega t) = \boldsymbol{x}_0 e^{\mathrm{i} \Omega t} + \bar{\boldsymbol{x}}_0 e^{- \mathrm{i} \Omega t} \end{align}.\] Accordingly, the non-autonomous part of the reduced dynamics truncated at leading order can be written as \[\boldsymbol{S}(\boldsymbol{p},\Omega t)= \boldsymbol{s}_0^+ e^{\mathrm{i} \Omega t} + \boldsymbol{s}_0^- e^{- \mathrm{i} \Omega t}.\] One can again use the invariance equation to solve for these unknown coefficients. Interested readers can refer to [28] for more details.

We express the reduced coordinates in polar form as \[\label{eq:transform} p=\rho e^{\mathrm{i}(\theta+\Omega t)}\tag{10}\] and plug it into the reduced dynamics in 7 . Then the third-order reduced dynamics can be expressed as [28] \[\begin{align} \dot{\rho} &= \mathrm{Re}(\lambda) \rho + \mathrm{Re}(\gamma) \rho^3 + \epsilon ( \mathrm{Re}(\tilde{f}) \cos \theta + \mathrm{Im}(\tilde{f}) \sin \theta), \\ \dot{\theta} &= \mathrm{Im}(\lambda) - \Omega + \mathrm{Im}(\gamma) \rho^2 + \epsilon ( \mathrm{Im}(\tilde{f}) \cos \theta - \mathrm{Re}(\tilde{f}) \sin \theta)/\rho, \end{align} \label{eq:reduc-nonauto}\tag{11}\] where \(\lambda\) is the eigenvalue of master subspace, i.e., \(\lambda_1\); \(\gamma\) is the third-order backbone coefficient; \(\tilde{f}\) is a modal force. Explicit expressions for \(\gamma\) and \(\tilde{f}\) are given as \[\gamma = -\boldsymbol{\psi}^{\mathrm{T}}\boldsymbol{f}_{21}, \quad \tilde{f} = 0.5\boldsymbol{\psi}^\mathrm{T}\boldsymbol{f}^{\mathrm{ext}}, \label{eq:compute95gammma95force}\tag{12}\] where \(\boldsymbol{f}_{21}\) is associated with cubic nonlinearity, given in the Eq.(19) of [28]. It is noted that a fixed point of the ROM 11 corresponds to a forced periodic orbit of the full system 1 , and the forced response curve (FRC) of the full system can be analytically extracted via the SSM-based ROM [25], [28].

2.3 Optimization problem statement↩︎

In this paper, we consider a few optimization problems associated to the tuning of the FRC. As we will see, these optimization problems enable the tuning of the peak and the hardening/softening behavior of FRC, and also the SN bifurcations on the FRC.

2.3.1 Minimize the peak of FRC with target backbone↩︎

By using SSM reduction, we transform the optimization of the nonlinear finite element system 1 into the optimization of a reduced system 11 . In this subsection, we will present the formulation for minimizing the peak of FRC in reduced coordinates.

Since the fixed point of the reduced dynamics 11 corresponds to a periodic orbit of the full system [25], [28], setting \((\dot{\rho}, \dot{\theta}) = \boldsymbol{0}\) and eliminating \(\theta\), we obtain [28] \[\left(\mathrm{Re}(\lambda)\rho + \mathrm{Re}(\gamma)\rho^3 \right)^2 + \left(\mathrm{Im}(\lambda) - \Omega + \mathrm{Im}(\gamma)\rho^2 \right)^2 \rho^2 = \epsilon^2 |\tilde{f}|^2. \label{eq:reduce-FRC}\tag{13}\] This algebraic equation gives the FRC in reduced coordinates \(\rho(\Omega)\), which can be further mapped to the FRC of periodic orbits of the full system via the periodic SSM map \(\boldsymbol{W}_\epsilon(\boldsymbol{p},\Omega t)\) in 7 .

Furthermore, the peak point \((\rho_{\rm{max}}, \Omega_{\rm{max}})\) on the reduced FRC \(\rho(\Omega)\) corresponds to the peak of FRC in physical coordinates [28], [32]. Let \(d\rho/d\Omega = 0\), the peak point in reduced coordinates can be expressed as [28] \[\begin{align} &\Omega_{\max} = \mathrm{Im}(\lambda) + \mathrm{Im}(\gamma)\, \rho_{\max}^2, \tag{14} \\ &\mathrm{Re}(\lambda)\, \rho_{\max} + \mathrm{Re}(\gamma)\, \rho_{\max}^3 = \epsilon |\tilde{f}|\, \mathrm{sign}(\mathrm{Re}(\lambda)). \tag{15} \end{align}\] Equation 14 provides the analytical expression of the backbone curve in reduced coordinates. Specifically, when \(\mathrm{Im}(\gamma) > 0\), we obtain \(\Omega_{\rm{max}} > \mathrm{Im}(\lambda)\), corresponding to a hardening behavior; when \(\mathrm{Im}(\gamma) < 0\), we obtain \(\Omega_{\rm{max}} < \mathrm{Im}(\lambda)\), corresponding to a softening behavior. The corresponding amplitude \(\rho_{\rm{max}}\) can be computed by solving Equation 15 .

Now, we aim to minimize the peak of FRC by setting \(\rho_{\rm{max}}\) as the objective function. Meanwhile, we add a constraint to \(\gamma\) such that the backbone curve is also tuned. The corresponding optimization problem is formulated as \[\begin{align} &\underset{\boldsymbol{\mu}}{\min} \quad c_{\mathrm{nl}} = \rho_{\max} \\ &\text{s.t.:} \\ &\quad \left(\frac{\mathrm{Im}(\gamma)}{\gamma_{\text{target}}} -1\right)^2 \leq \varepsilon^2\\ &\quad (\omega_Y/\omega_{Y,\text{target}} - 1)^2 \leq \varepsilon^2 \\ &\quad \omega_X \geq \omega_{X,\text{target}} \\ &\quad A \leq A_{\max}. \\ \end{align} \label{eq:opt-rho95max}\tag{16}\] Here, the constraint on \(\mathrm{Im}(\gamma)\) is added to control the hardening/softening behavior. We also add constraints on \(\omega_{\text{Y}}\) (i.e., \(\omega_1\)) and \(\omega_{\text{X}}\) (i.e., \(\omega_2\)) to prevent the occurrence of internal resonance in the first and second modes. Specifically, the condition \(\omega_{X,\text{target}} > 3\omega_{Y,\text{target}}\) needs to be satisfied since the highest-order nonlinear term in the system is cubic. Further, the structural area fraction \(A\) is subject to an upper bound.

We note that the optimization problem 16 provides a concurrent design of the FRC, namely, the peak of the FRC and the softening or hardening of the associated backbone curve. To illustrate the effectiveness of the formulation 16 , we also consider two reference optimization problems: an optimization problem on minimizing the peak of FRC without tuning the backbone, and an optimization problem solely on tuning the backbone.

For the first reference problem, we consider a linear optimization formulation for periodic response [20], [33], namely \[\begin{align} &\underset{\boldsymbol{\mu}}{\min} \quad c_\mathrm{lin} = |\boldsymbol{U}|^{\mathrm{T}}\boldsymbol{L}|\boldsymbol{U}| \\ &\text{s.t.:} \\ &\quad (\omega_Y/\omega_{Y,\text{target}} - 1)^2 \leq \varepsilon^2 \\ &\quad \omega_X \geq \omega_{X,\text{target}} \\ &\quad A \leq A_{\max}, \\ \end{align} \label{eq:opt-linear}\tag{17}\] where \(\boldsymbol{U}\) denotes the linear displacement vector obtained from \((\boldsymbol{K} + \mathrm{i} \Omega \boldsymbol{C} - \Omega^2 \boldsymbol{M}) \boldsymbol{U} = \boldsymbol{f}^{\text{ext}}\), and the symbol \(|\cdot|\) represents taking the absolute value of each component of the vector. Here, \(\boldsymbol{L}\) is a diagonal selection matrix with ones at DOFs corresponding to the node, line, or domain of interest [33]. To maintain consistency, 17 has the same constraints as 16 , except for the constraint on \(\mathrm{Im}(\gamma)\), since linear optimization cannot capture or control the hardening/softening behavior. We will give an example to compare the optimization result of formulations 16 and 17 in Sec. 5.1.

As for the second reference problem, we consider the following optimization problem in which a target value of \(\mathrm{Im}(\gamma)\) is included in the objective function: \[\begin{align} &\underset{\boldsymbol{\mu}}{\min} \quad c_{\gamma} = \left(\frac{\mathrm{Im}(\gamma)}{\gamma_{\text{target}}} -1\right)^2 \\ &\text{s.t.:} \\ &\quad (\omega_Y/\omega_{Y,\text{target}} - 1)^2 \leq \varepsilon^2 \\ &\quad \omega_X \geq \omega_{X,\text{target}} \\ &\quad A \leq A_{\max}. \\ \end{align} \label{eq:opt-gamma}\tag{18}\] It is clear from the objective function that \(c_\gamma = 0\) is achieved when \(\mathrm{Im}(\gamma) = \gamma_{\text{target}}\). The constraints in the formulation 18 are consistent with those in 16 , ensuring a fair comparison between the two formulations. We will give an example to compare the result of formulations 16 and 18 in Sec. 5.2.

2.3.2 Control of SN bifurcations on FRC↩︎

As discussed in [25], [34], the SN bifurcations of the fixed points of the reduced dynamics 11 correspond to those of the periodic orbits of the full system 1 . This correspondence allows us to control the occurrence of SN bifurcations in reduced coordinates. In this subsection, we derive the expression of the SN bifurcation points in the reduced coordinates and eliminate their occurrence by enforcing the degeneracy condition of a cusp bifurcation. In the final step, we formulate an optimization problem to suppress the occurrence of SN bifurcations.

We first need to find the fixed point of SN bifurcation \((\rho_{\mathrm{SN}}, \Omega_{\mathrm{SN}})\) in reduced coordinates. The Jacobian of the reduced dynamics 11 can be expressed as [28] \[\boldsymbol{A}(\rho,\Omega) = \begin{pmatrix} \mathrm{Re}(\lambda) + 3\mathrm{Re}(\gamma) \rho^2 & -\left(\mathrm{Im}(\lambda) - \Omega + \mathrm{Im}(\gamma) \rho^2\right) \rho \\ 2\mathrm{Im}(\gamma) \rho + \left(\mathrm{Im}(\lambda) - \Omega + \mathrm{Im}(\gamma) \rho^2 \right)/\rho & \mathrm{Re}(\lambda) + \mathrm{Re}(\gamma) \rho^2 \end{pmatrix}.\] At the SN bifurcation of a fixed point, the Jacobian \(\boldsymbol{A}\) has a simple zero eigenvalue [35]. Thus, the determinant of \(\boldsymbol{A}\) is 0 at SN bifurcation, namely \[\begin{align} \left(\mathrm{Re}(\lambda) + 3\mathrm{Re}(\gamma) \rho^2\right)\left(\mathrm{Re}(\lambda) + \mathrm{Re}(\gamma) \rho^2\right) + \left(\mathrm{Im}(\lambda) - \Omega + \mathrm{Im}(\gamma) \rho^2\right) \\ \left(2\mathrm{Im}(\gamma) \rho^2 + \left(\mathrm{Im}(\lambda) - \Omega + \mathrm{Im}(\gamma) \rho^2 \right)\right) =0. \end{align} \label{eq:SN-deteq0}\tag{19}\] By combining 13 and 19 , and eliminating \(\Omega\), we derive the following algebraic equation for the response amplitude \(\rho_{\mathrm{SN}}\) in the reduced coordinates \[\begin{align} \left( 4\mathrm{Re}(\gamma)^4 + 4\mathrm{Im}(\gamma)^2 \mathrm{Re}(\gamma)^2 \right) \rho_{\mathrm{SN}}^{12} + 8\mathrm{Re}(\lambda)\mathrm{Re}(\gamma) \left( \mathrm{Re}(\gamma)^2 + \mathrm{Im}(\gamma)^2 \right) \rho_{\mathrm{SN}}^{10} \\ + 4\mathrm{Re}(\lambda)^2 \left( \mathrm{Re}(\gamma)^2 + \mathrm{Im}(\gamma)^2 \right) \rho_{\mathrm{SN}}^8 + 4\epsilon^2 |\tilde{f}|^2 \left( \mathrm{Re}(\gamma)^2 - \mathrm{Im}(\gamma)^2 \right) \rho_{\mathrm{SN}}^6 \\ + 4\epsilon^2 |\tilde{f}|^2 \mathrm{Re}(\lambda)\mathrm{Re}(\gamma) \rho_{\mathrm{SN}}^4 + \epsilon^4 |\tilde{f}|^4 = 0. \end{align} \label{eq:rho95SN}\tag{20}\] Substituting the value of \(\rho_{\mathrm{SN}}\) into 13 , we obtain the corresponding frequency \(\Omega_{\mathrm{SN}}\) as: \[\begin{align} \Omega_{\mathrm{SN}} &= \mathrm{Im}(\lambda) + \mathrm{Im}(\gamma) \rho_{\mathrm{SN}}^2 + \mathrm{sign}(\mathrm{Im}(\gamma)) k_{\mathrm{SN}}, \end{align}\] where \[\begin{align} k_{\mathrm{SN}} = \sqrt{{\epsilon^2 |\tilde{f}|^2}/{\rho_{\mathrm{SN}}^2} - \left( \mathrm{Re}(\lambda(\mu)) + \mathrm{Re}(\gamma(\mu)) \rho_{\mathrm{SN}}^2 \right)^2}. \end{align}\]

Note that there are usually two SN bifurcation points on the FRC. Next, we aim to suppress the occurrence of these two SN bifurcation points by enforcing the cusp bifurcation degeneration condition, namely [35] \[\begin{align} b = \langle \boldsymbol{\psi}_{\mathrm{SN}}, \;\boldsymbol{B}_{\mathrm{SN}}(\boldsymbol{\phi}_{\mathrm{SN}},\boldsymbol{\phi}_{\mathrm{SN}})\rangle = 0. \label{eq:compute95b} \end{align}\tag{21}\] where \(\boldsymbol{\psi}_{\mathrm{SN}}\) and \(\boldsymbol{\phi}_{\mathrm{SN}}\) are the left and right eigenvectors of \(\boldsymbol{A}_{\mathrm{SN}}=\boldsymbol{A}(\rho_{\mathrm{SN}},\Omega_{\mathrm{SN}})\), \(\boldsymbol{B}_{\mathrm{SN}} \in \mathbb{R}^{2 \times 2 \times 2}\) is the second derivative tensor of reduced dynamics 11 at \((\rho_{\mathrm{SN}},\Omega_{\mathrm{SN}})\). Specifically, we have \[\label{eq:comp95A95B} \boldsymbol{A}_{\mathrm{SN}}\boldsymbol{\phi}_{\mathrm{SN}}=0, \quad \boldsymbol{A}^\mathrm{T}_{\mathrm{SN}}\boldsymbol{\psi}_{\mathrm{SN}}=0,\quad \boldsymbol{B}_{\mathrm{SN}}(\boldsymbol{\phi}_{\mathrm{SN}}, \boldsymbol{\phi}_{\mathrm{SN}}) = \begin{pmatrix} \boldsymbol{\phi}_{\mathrm{SN}}^{\mathrm{T}} \boldsymbol{B}_1 \boldsymbol{\phi}_{\mathrm{SN}} \\ \boldsymbol{\phi}_{\mathrm{SN}}^{\mathrm{T}} \boldsymbol{B}_2 \boldsymbol{\phi}_{\mathrm{SN}} \end{pmatrix},\tag{22}\] where \[\boldsymbol{B}_1 = \begin{pmatrix} 6\mathrm{Re}(\gamma)\rho_{\mathrm{SN}} & 0 \\ 0 & \mathrm{Re}(\lambda)\rho_{\mathrm{SN}}+\mathrm{Re}(\gamma) \rho_{\mathrm{SN}}^3 \end{pmatrix},\, \boldsymbol{B}_2 = \begin{pmatrix} 2\mathrm{Im}(\gamma) - 2b_{2,1}/\rho_{\mathrm{SN}}^2 & b_{2,2} \\ b_{2,2} & b_{2,1} \end{pmatrix}\] with \(b_{2,1}=\mathrm{Im}(\lambda)-\Omega_{\mathrm{SN}}+\mathrm{Im}(\gamma)\rho_{\mathrm{SN}}^2\) and \(b_{2,2}=-(\mathrm{Re}(\lambda)+\mathrm{Re}(\gamma)\rho_{\mathrm{SN}}^2)/\rho_{\mathrm{SN}}\). Since \(\det(\boldsymbol{A}) = 0\), the eigenvectors \(\boldsymbol{\psi}_{\mathrm{SN}}\) and \(\boldsymbol{\phi}_{\mathrm{SN}}\) are not uniquely determined. To ensure uniqueness, we impose the following normalization conditions \[\label{eq:norm95vector} \boldsymbol{\psi}_{\mathrm{SN}}^{\mathrm{T}}\boldsymbol{\psi}_{\mathrm{SN}}=1, \quad\boldsymbol{\psi}_{\mathrm{SN}}^{\mathrm{T}}\boldsymbol{\phi}_{\mathrm{SN}}=1.\tag{23}\]

Many MEMS sensors (e.g., accelerometers, gyroscopes) are designed to operate within the linear regime to avoid hysteresis and excess noise caused by nonlinear effects [2]. To ensure such behavior, we aim to control the system’s bifurcation structure while enhancing its performance. Specifically, we formulate an optimization problem that maximizes the response amplitude \(\rho_{\max}\), subject to a constraint on the parameter \(b\) to limit the distance between the two SN bifurcation points. The resulting optimization problem is given as \[\begin{align} &\underset{\boldsymbol{\mu}}{\max} \quad c_{\mathrm{SN}} = \rho_{\max} \\ &\text{s.t.:} \\ &\quad 0 \leq |b| \leq b_\text{target} \\ &\quad \omega_X \geq \omega_{X,\text{target}} \\ &\quad (A/A_\text{target} - 1)^2 \leq \varepsilon^2, \\ \end{align} \label{eq:opt-SNbifur}\tag{24}\] where the constant \(b_{\text{target}}\) regulates the distance between the two SN bifurcation points. As \(b_{\text{target}}\) approaches zero, these two points coalesce, resulting in a cusp bifurcation. In addition, we enforce a frequency constraint on the X mode to ensure connectivity [36]. A two-sided constraint is also imposed on the structural area to improve the stability of the optimization process.

3 Sensitivity analysis↩︎

In the previous section, we formulated several optimization problems. Since the Method of Moving Asymptotes (MMA) is used to solve it, sensitivity analysis with respect to the design variables is essential for updating the design. In this section, we derive the sensitivity expressions for the key quantities \(\omega_j\), \(\lambda\), \(\gamma\), \(\tilde{f}\), \(\rho_{\max}\), and \(b\), which are critical for the optimization process.

The sensitivity of the eigenfrequencies \(\omega_j\) and eigenvalues \(\lambda\) has already been derived in Appendix.C of [28], namely \[\begin{align} \omega_j' = \frac{\boldsymbol{\phi}^{\text{T}}(\boldsymbol{K}' - \omega_j^2 \boldsymbol{M}')\boldsymbol{\phi}}{2\omega_j}, \\ \lambda' = -\frac{\xi' \omega \lambda + (\xi \lambda + \omega) \omega'}{\lambda + \xi \omega}. \end{align}\] where \(\omega\) is the eigenfrequency of the master subspace, i.e., \(\omega_1\). Here and throughout this paper, the apex \('\) denotes the derivative with respect to one of the system parameters contained in the vector \(\boldsymbol{\mu}\) introduced in 1 .

Following [28], the variation of \(\gamma\) can be obtained by the adjoint method, which can be expressed as \[\begin{align} \delta \gamma = & - \boldsymbol{\lambda}_{21}^{\mathrm{T}} \Big( \partial_{\boldsymbol{\mu}} \boldsymbol{f}_2(\boldsymbol{\mu}, \boldsymbol{\phi}, \boldsymbol{W}_{20}^{(1)}) + \partial_{\boldsymbol{\mu}} \boldsymbol{f}_2(\boldsymbol{\mu}, \boldsymbol{W}_{20}^{(1)}, \boldsymbol{\phi}) \\ & + \partial_{\boldsymbol{\mu}} \boldsymbol{f}_2(\boldsymbol{\mu}, \boldsymbol{\phi}, \boldsymbol{W}_{11}^{(1)}) + \partial_{\boldsymbol{\mu}} \boldsymbol{f}_2(\boldsymbol{\mu}, \boldsymbol{W}_{11}^{(1)}, \boldsymbol{\phi}) \delta \boldsymbol{\mu} \Big) \\ & - 3 \boldsymbol{\lambda}_{21}^{\mathrm{T}} \partial_{\boldsymbol{\mu}} \boldsymbol{f}_3(\boldsymbol{\mu}, \boldsymbol{\phi}, \boldsymbol{\phi}, \boldsymbol{\phi}) \delta \boldsymbol{\mu}\\ & + \boldsymbol{\lambda}_{20}^{\mathrm{T}} \Big( (4\lambda^2 + 2\lambda \alpha)\delta \boldsymbol{M} + (2\lambda \beta + 1)\delta \boldsymbol{K} \Big) \boldsymbol{W}_{20}^{(1)} \\ & + \boldsymbol{\lambda}_{11}^{\mathrm{T}} \Big( (4 [\mathrm{Re}(\lambda)]^2 + 2\mathrm{Re}(\lambda) \alpha) \delta \boldsymbol{M} + (2\mathrm{Re}(\lambda)\beta + 1)\delta \boldsymbol{K} \Big) \boldsymbol{W}_{11}^{(1)} \\ & + (\boldsymbol{\lambda}_{20} + \boldsymbol{\lambda}_{11})^{\mathrm{T}} \partial_{\boldsymbol{\mu}} \boldsymbol{f}_2(\boldsymbol{\mu}, \boldsymbol{\phi}, \boldsymbol{\phi}) \delta \boldsymbol{\mu} \\ & + \boldsymbol{\lambda}_{\boldsymbol{\phi}}^{\mathrm{T}} \Big( \delta \boldsymbol{K} - \omega^2 \delta \boldsymbol{M} \Big) \boldsymbol{\phi} + \lambda_{\mathrm{norm}} \boldsymbol{\phi}^{\mathrm{T}} \delta \boldsymbol{M} \boldsymbol{\phi}. \end{align}\] In the above equation, the damping coefficients \(\alpha\) and \(\beta\) are fixed, while the damping ratio \(\xi\) varies throughout the optimization process. This differs from the formulation in [28], where the damping ratio remains constant. We follow the tensorial approach in [12] to ensure efficient implementation.

The variation of \(\tilde{f}\) can also be obtained using the adjoint method. As shown in Appendix 7, the sensitivity of \(\tilde{f}\) is given by \[\begin{align} \delta \tilde{f} = \boldsymbol{\eta}_{\boldsymbol{\phi}}^{\mathrm{T}} \bigl(\delta\boldsymbol{K}\boldsymbol{\phi} - \omega^2 \delta \boldsymbol{M}\boldsymbol{\phi}\bigr) + \eta_{\mathrm{norm}} \boldsymbol{\phi}^{\mathrm{T}} \delta \boldsymbol{M} \boldsymbol{\phi} \end{align},\] where \(\boldsymbol{\eta}_{\boldsymbol{\phi}}\) and \(\eta_{\mathrm{norm}}\) are obtained from 25 . The sensitivity of \(\rho_{\max}\) can be obtained by taking derivative of 15 , resulting in \[\begin{align} \rho_{\max}' = \frac{\epsilon|\tilde{f}|' \mathrm{sign}(\mathrm{Re}(\lambda)) - \mathrm{Re}(\lambda')\rho_{\max} - \mathrm{Re}(\mathrm{\gamma'})\rho_{\max}^3}{\mathrm{Re}(\lambda) + 3\mathrm{Re}(\gamma)\rho_{\max}^2} \end{align}.\] Similarly, by taking derivative of 21 , the sensitivity of \(b\) can be expressed as \[\begin{align} b'=\boldsymbol{\psi}_{\mathrm{SN}}'\boldsymbol{B}_{\mathrm{SN}}(\boldsymbol{\phi}_{\mathrm{SN}},\boldsymbol{\phi}_{\mathrm{SN}}) + \boldsymbol{\psi}_{\mathrm{SN}}\boldsymbol{B}_{\mathrm{SN}}'(\boldsymbol{\phi}_{\mathrm{SN}},\boldsymbol{\phi}_{\mathrm{SN}})+ 2\boldsymbol{B}_{\mathrm{SN}}(\boldsymbol{\phi}_{\mathrm{SN}},\boldsymbol{\phi}_{\mathrm{SN}}') \end{align},\] where the sensitivities of \(\boldsymbol{\psi}_{\mathrm{SN}}\), \(\boldsymbol{\phi}_{\mathrm{SN}}\), and \(\boldsymbol{B}_{\mathrm{SN}}\) are derived in Appendix 8.

4 SIMP Model for Topology Optimization↩︎

In this work, the Solid Isotropic Material with Penalization (SIMP) method [29] is employed as the foundation for density-based topology optimization. In the SIMP-based framework, the material properties of the structure are interpolated element-wise based on the physical density variables \(\hat{\mu}_e\). These densities are used to assemble the global matrices and nonlinear force vectors. Specifically, the global mass and stiffness matrices, as well as the quadratic and cubic stiffness tensors, are computed as follows [12]: \[\begin{align} \boldsymbol{M} &= \bigcup_{e=1}^{N_{\mathrm{el}}} \hat{\mu}_e \boldsymbol{M}_e, \quad \boldsymbol{K} = \bigcup_{e=1}^{N_{\mathrm{el}}} \hat{\mu}_e \boldsymbol{K}_e, \\ F_2 &= \bigcup_{e=1}^{N_{\mathrm{el}}} \hat{\mu}_e F_{2e}, \quad F_3 = \bigcup_{e=1}^{N_{\mathrm{el}}} \hat{\mu}_e F_{3e}, \end{align}\] where the symbol \(\bigcup\) denotes the standard finite element assembly process, and \(N_{\mathrm{el}}\) is the total number of finite elements. The quadratic tensor \(F_2\) and cubic tensor \(F_3\) are used to compute the nonlinear internal force via Einstein notation, e.g., \[\begin{align} f_2^i(\boldsymbol{a},\boldsymbol{b})=F_2^{ijk}a^jb^k, \quad f_3^i(\boldsymbol{a},\boldsymbol{b},\boldsymbol{c})=F_3^{ijkl}a^jb^kc^l, \end{align}\] where the vectors \(\boldsymbol{a}\), \(\boldsymbol{b}\), and \(\boldsymbol{c}\) are the inputs of the nonlinear force. The physical densities \(\hat{\mu}_e\) are obtained by applying filtering and projection operations to the design variables \(\mu_e\), as described below.

To ensure mesh-independent and well-posed topology optimization results, the design variable \(\mu \in [0, 1]\) is first smoothed using a density filter [30], defined as \[\tilde{\mu}_e = \frac{\sum_{j \in \mathcal{N}_e} w_{j,e} \mu_j}{\sum_{j \in \mathcal{N}_e} w_{j,e}},\] where \(\tilde{\mu}_e\) is the filtered density of element \(e\), \(\mathcal{N}_e\) is the set of neighboring elements within a radius \(R\), and \(w_{j,e} = R - \lvert \mathbf{x}_j - \mathbf{x}_e \rvert\) is the weight based on centroid distance. Unless otherwise specified, the filter radius is set to \(R = 4\).

Next, a Heaviside-type projection [30] is applied to obtain the projected density \(\bar{\mu}_e\), which sharpens the transition between solid and void regions: \[\bar{\mu}_e = \frac{\tanh(\sigma \eta) + \tanh(\sigma (\tilde{\mu}_e - \eta))}{\tanh(\sigma \eta) + \tanh(\sigma (1 - \eta))},\] where \(\sigma\) is the projection steepness parameter and \(\eta \in (0, 1)\) is the threshold. In this work, we initially set \(\sigma = 10\) and \(\eta = 0.5\), and increase \(\sigma\) progressively during the optimization, unless otherwise specified.

Finally, the physical density \(\hat{\mu}_e\) used to interpolate material properties is computed via the SIMP scheme: \[\hat{\mu}_e = \hat{\mu}_0 + (1 - \hat{\mu}_0)\bar{\mu}_e^p,\] where \(\hat{\mu}_0 \ll 1\) is a small lower bound to avoid singular stiffness in void regions, and \(p\) is the penalization exponent to promote 0-1 designs. In this work, the value of \(\hat{\mu}_0\) is set to \(10^{-6}\), while \(p\) is initially set to 1 and gradually increased during the optimization process.

5 Numerical examples↩︎

In this section, three numerical examples are presented to show the result of minimizing the peak of FRC, controlling hardening/softening behavior, and suppressing the occurrence of SN bifurcations. In all the examples, the 2-dimensional structures are discretized using 4-node, square elements. Following [12], material properties are specified with the density \(\rho = 2330 \times 10^{-6} \;\rm{ng/um^3}\), Young’s modulus \(E=148 \times 10^9 \;\rm{Pa}\), and Poisson’s ratio \(\nu = 0.23\). The plane stress approximation is used, with a thickness of 24 \(\rm{\mu m}\).

The matrix assembly of the finite element model and the solution of the optimization problems are carried out using the YetAnotherFEcode package [37]. The FRCs of the optimized layouts are computed using SSMTool [38], an open-source software package for the computation of invariant manifolds and their reduced dynamics. In the Sec. 2.2 and 2.3.1, we employed a third-order reduced-order model 11 to obtain the expression of FRC in the reduced coordinates 13 . To assess the reliability of the optimized result, we further compute the FRC using SSMTool at higher-order truncations for the optimized structure. If the FRC at higher-order truncations agrees well with that at the third-order truncation, it indicates that the predicted FRC has converged and is hence reliable.

5.1 Minimizing the peak of FRC: linear vs. nonlinear optimizations↩︎

We consider a Messerschmitt-Bölkow-Blohm (MBB) beam with a fixed segment located at its center. Due to structural symmetry, the optimization is performed only on the left half of the beam [12], as illustrated in Fig. 1. The half MBB beam measures 800 \(\mathrm{\mu m}\) in length and 100 \(\mathrm{\mu m}\) in height, including a fixed region with a length of 100 \(\mathrm{\mu m}\). The left end of the beam is fully fixed, and the right end is constrained in the \(x\)-direction. A periodic excitation is applied at the midpoint of the right edge.

A finite element mesh of \(160 \times 20\) elements is employed, resulting in a system with \(n = 6699\) degrees of freedom in 1 . The damping ratio of the initial layout is set to \(\xi\) = 0.1%, which yields the damping constants \(\alpha = 2.99\) and \(\beta = 1.72 \times 10^{-7}\). During the optimization process, the forcing amplitude is given as \(f^{\mathrm{ext}} = 5 \times 10^{9} \;\mathrm{ng \cdot \mu m/ms^2}\). Unless otherwise specified, the parameter \(\epsilon\) is set to \(0.01\) throughout the computation. The target values for other constraints are set as follows: \(\omega_{Y,\text{target}} = 300 ~\mathrm{kHz}\), \(\omega_{X,\text{target}} = 1200 ~\mathrm{kHz}\), and \(A_{\max} = 50\%\). Based on these parameter settings, we compare the results obtained from the nonlinear optimization formulation 16 with those from the linear optimization formulation 17 .

With \(\gamma_{\mathrm{target}} = 5 \times 10^{-5}\), the iterative process for the nonlinear optimization formulation 16 is provided in Fig. 2. The objective \(\rho_{\max}\) exhibits an initial increase, followed by a rapid decrease until convergence. Meanwhile, the constraint on \(\mathrm{Im}(\gamma)\) is violated initially, and the value of \(\mathrm{Im}(\gamma)\) gradually converges and ultimately reaches the target \(\mathrm{Im}(\gamma) = \gamma_{\mathrm{target}}\). The abrupt changes observed in the objective \(\rho_{\max}\) and \(\mathrm{Im}(\gamma)\) during the iteration are caused by the scheduled updates of the penalization exponent \(p\) and the projection parameter \(\beta\) in the SIMP algorithm. The obtained optimal layout is given in the upper panel of Fig. 3. The computational performance of this iterative process is summarized in Table 1.

We solve the same nonlinear optimization formulation 16 but with varied \(\gamma_{\mathrm{target}}\). The optimal layouts obtained for different values of \(\gamma_{\mathrm{target}}\), as well as from the linear optimization formulation 17 , are shown in Fig. 3. As \(\gamma_{\mathrm{target}}\) varies, the optimal layouts from the nonlinear optimization exhibit significant changes. Moreover, all layouts obtained from the nonlinear optimization differ noticeably from the linear layout.

We now compare the FRCs of the optimal layouts shown in Fig. 3. The vertical axis \(\Vert z_{\rm{out1}} \Vert_\infty\) represents the vibration amplitude of the node located at (\(x\),\(y\))=(800, 50) in this example. As shown in the left panel of Fig. 4, the FRCs corresponding to the nonlinear optimization exhibit distinct hardening or softening behaviors depending on the value of \(\gamma_{\mathrm{target}}\). To make these nonlinear characteristics more evident, we increase the forcing amplitude to \(f^{\mathrm{ext}} = 2 \times 10^{10} \;\mathrm{ng \cdot \mu m/ms^2}\) and recompute the FRCs for the optimized structures. The results, presented in the right panel of Fig. 4, show that the FRC behavior gradually shifts from softening to hardening as \(\gamma_{\mathrm{target}}\) increases. Thus, while both linear and nonlinear optimization approaches succeed in controlling the FRC amplitude, only the nonlinear optimization is capable of adjusting the hardening/softening behavior. A posteriori computation shows that the structure obtained from the linear optimization yields \(\mathrm{Im}(\gamma) = 3.89 \times 10^{-5}\). This value is close to \(5\times10^{-5}\), which explains the FRC for the \(\gamma_{\mathrm{target}}=5\times10^{-5}\) case looks similar to the linear one, as seen in Fig. 4.

Since third-order SSM reduction is employed to compute the FRCs of the optimal layouts, it is necessary to assess the convergence of the SSM-based model reduction. The convergence analysis is presented in Appendix 9, which shows that the FRC obtained using third-order SSM reduction agrees well with that obtained using seventh-order SSM reduction under the forcing amplitude \(f^{\mathrm{ext}} = 2 \times 10^{10} \;\mathrm{ng \cdot \mu m/ms^2}\).

5.2 Tailoring hardening/softening behavior: concurrent design vs. backbone design only↩︎

Tailoring the hardening/softening behavior of the backbone curve in an MBB beam has been discussed in [12]. In this subsection, we aim to minimize the response amplitude while simultaneously tailoring the softening–hardening behavior, namely, the concurrent design 16 . We further compare this concurrent design with the design given by the reference problem 18 that tunes the backbone only. As shown in Fig. 5, the design domain of the beam has a length of 500 \(\mathrm{\mu m}\) and a height of 100 \(\mathrm{\mu m}\). A fixed region, occupying 20% of the total beam length, is located at the center and serves as a proof mass. The left end of the beam is fully fixed, while the right end is constrained in the \(x\)-direction. A periodic excitation is applied at the midpoint of the right edge.

The FE model uses a mesh of \(100 \times 20\) elements, leading to a system with \(n\)= 4179 DOF in 1 . The damping ratio of initial layout \(\xi\) = 0.1%, which leads to the values of damping constants \(\alpha = 6.90\) and \(\beta = 7.41 \times 10^{-8}\). During the optimization iteration process, the forcing amplitude is given as \(f^{\mathrm{ext}}= 5 \times 10^{9} \;\mathrm{ng \cdot \mu m/ms^2}\). In the following computation, we set \(\epsilon = 0.01\) unless otherwise stated. The projection parameter is fixed at \(\sigma = 10\) throughout the optimization, and the filter radius is set to \(R = 3\). In this example, no constraint is imposed on \(\omega_X\) since \(\omega_X > 3\omega_Y\) is naturally satisfied throughout the optimization iterations. Furthermore, the two-sided constraint on \(\gamma\) is replaced by a one-sided constraint. Specifically, \(\gamma > \gamma_{\mathrm{target}}\) is enforced when \(\gamma_{\mathrm{target}} > 0\), and \(\gamma < \gamma_{\mathrm{target}}\) is enforced when \(\gamma_{\mathrm{target}} < 0\). The target values for other constraints are set as follows: \(\omega_{Y,\text{target}} = 600 ~\mathrm{kHz}\), and \(A_{\max} = 50\%\). Based on these parameter settings, we compare the optimization results obtained from formulations 16 and 18 for \(\gamma_{\mathrm{target}} = \pm 1 \times 10^{-3}\).

5.2.1 Hardening optimization↩︎

At \(\gamma_{\mathrm{target}} = 1 \times 10^{-3}\), the optimal layouts for \(c_{\mathrm{nl}}\) in 16 and \(c_{\gamma}\) in 18 are shown in Fig. 6. The optimal layouts of \(c_{\mathrm{nl}}\) and \(c_{\gamma}\) are both symmetrical structures with respect to the axis \(y = 50\). Although both designs exhibit global symmetry, noticeable local differences can be observed between the two layouts. These differences reflect the influence of the chosen optimization objective on the final layout.

We then compare the FRCs of two optimal structures shown in Fig. 6 using SSMTool. The lower-left panel of Fig. 6 presents the FRCs computed under a forcing amplitude of \(f^{\mathrm{ext}} = 5 \times 10^{9}\), where the vertical axis \(\Vert z_{\rm{out1}} \Vert_\infty\) denotes the vibration amplitude at node (\(x\), \(y\)) = (500, 50). At this excitation level, the FRC peak of \(c_\mathrm{nl}\) is lower than that of \(c_{\gamma}\). Moreover, the hardening behavior of the FRCs is not prominent, as the forcing is too weak to activate significant nonlinear effects in the optimized layouts. Therefore, we increase the forcing amplitude to \(f^{\mathrm{ext}} = 2 \times 10^{10}\) and recompute the FRCs for both designs. The resulting curves are shown in the lower-middle panel of Fig. 6, where the FRC peak of \(c_\mathrm{nl}\) remains lower than that of \(c_{\gamma}\). Additionally, as the forcing amplitude increases, two SN bifurcation points are observed on the FRCs.

Since one of the SN bifurcation points lies close to the peak of the FRC, the variation of the FRC peak with respect to the forcing amplitude \(\epsilon\) can be studied via the SN bifurcation curves. As seen in the lower-right panel of Fig. 6, we allow for the variations in the forcing amplitude with \(\epsilon = [0.01, 1.2]\epsilon_o\) with \(\epsilon_o= 0.01\) to investigate how the SN bifurcations evolve with varying \(\epsilon\). We can see that the variation amplitude of SN bifurcations of optimal layout \(c_\mathrm{nl}\) is lower than that of \(c_{\gamma}\). Consequently, as the forcing amplitude increases, the peak of FRC of \(c_\mathrm{nl}\) remains lower than that of \(c_{\gamma}\). To verify the convergence of the reduced-order model, we recomputed the SN bifurcation curves using an SSM reduction truncated at \(\mathcal{O}(7)\). As shown in the same panel, the results confirm the convergence. Finally, the computational performance for \(\gamma_{\mathrm{target}} = 1 \times 10^{-3}\) is summarized in Table 2.

5.2.2 Softening Optimization↩︎

Now we consider the case of \(\gamma_{\mathrm{target}}<0\), which corresponds to softening behavior of FRCs. At \(\gamma_{\mathrm{target}} = -1 \times 10^{-3}\), the optimal layouts of the optimization problems 16 and 18 are shown in Fig. 7. Unlike the case of hardening optimization, the optimal layouts of \(c_{\gamma}\) and \(c_\mathrm{nl}\) are both asymmetrical structures. Further, the optimal layout of \(c_\mathrm{nl}\) is different from that of \(c_{\gamma}\).

We then compare the FRCs of these two optimal layouts. The FRCs under forcing amplitude \(f^{\mathrm{ext}}= 5 \times 10^{9}\) are shown in the lower-left panel of Fig 7, from which we can see that the peak of FRC of \(c_\mathrm{nl}\) is lower than that of \(c_{\gamma}\). To observe the softening behavior more clearly, we increase the forcing amplitude to \(f^{\mathrm{ext}} = 2 \times 10^{10}\). As seen in the lower-middle panel of Fig 7, the peak of FRC of \(c_\mathrm{nl}\) is still lower than that of \(c_{\gamma}\). Also, the softening behavior can be observed clearly for both optimal layouts. Further, two SN bifurcation points emerge as the forcing amplitude increases.

Similar to the case of hardening optimization, we perform parameter continuation of the SN bifurcation points with varying forcing amplitude \(\epsilon = [0.01,1.2]\epsilon_o\). The projections of the obtained SN bifurcation curves are shown in the lower-right panel of Fig 7, from which we can see that the vibration amplitude of SN bifurcation points of \(c_\mathrm{nl}\) is lower than that of \(c_{\gamma}\). Therefore, the optimized structure \(c_\mathrm{nl}\) demonstrates consistently enhanced resistance to external disturbances compared to \(c_{\gamma}\) under various forcing amplitudes. Finally, we examine the computational performance for the case of \(\gamma_{\mathrm{target}} = -1 \times 10^{-3}\), as summarized in Table 3.

5.3 Tuning saddle-node (SN) bifurcations on FRC↩︎

The final example is a microbeam in a mass sensor [39]. The goal is to maximize the vibration amplitude of the microbeam to improve its measurement accuracy. Meanwhile, to prevent excessive amplitude fluctuations, we control the SN bifurcation points on the FRC to maximize the region without bistability. As seen in Fig. 8, the design domain of the beam has a length of 600 \(\mathrm{\mu m}\) and a height of 100 \(\mathrm{\mu m}\). A fixed region in the middle accounts for 20% of the design domain. The left end of the beam is fixed, while the right end is constrained in the \(x\)-direction. A periodic excitation is applied at the center of the right edge.

Similar to the case of the MBB beam, a mesh of \(120 \times 20\) elements is used in the FE model, leading to a system with \(n\)= 5019 DOF in 1 . The damping ratio of initial layout \(\xi\) = 0.1%, which leads to the values of damping constants \(\alpha = 1.97\) and \(\beta = 2.77 \times 10^{-7}\). During the optimization iteration process, the forcing amplitude is given as \(f^{\mathrm{ext}}= 2 \times 10^{8} \;\mathrm{ng \cdot \mu m/ms^2}\). In the following computation, we set \(\epsilon = 0.01\) unless otherwise stated. The target values for constraints in formulation 24 are set as \(\omega_{X,\text{target}} = 1500 ~\mathrm{kHz}\), and \(A_{\mathrm{target}} = 40\%\). Based on these parameter settings, we present the optimization result of optimization formulation 24 with different values of \(b_{\mathrm{target}}\).

We consider three cases \(b_{\mathrm{target}}\in\{2,1,0.1\}\). The obtained results are shown in Fig. 9. From the first row, we observe that the optimal layouts evolve with decreasing \(b_{\mathrm{target}}\), showing both changes in topology and clearer material boundaries. The second row shows the variation of \(|b|\) during the iterative process. At the beginning of the iteration, the value of \(|b|\) is 0, which indicates that the SN bifurcation point has not yet occurred. As the number of iteration steps increases, the vibration amplitude of the system becomes larger, the SN bifurcation points emerge, the value of \(|b|\) is greater than 0, and it finally tends to the value of \(b_{\mathrm{target}}\). Further, as the value of \(b_{\mathrm{target}}\) decreases, more iterations are required for the convergence of \(|b|\). Finally, we consider the FRCs of the three optimal layouts, as shown in the third row of Fig. 9. In these plots, the vertical axis \(\Vert z_{\rm{out1}} \Vert_\infty\) denotes the vibration amplitude at the node located at \((x, y) = (600, 50)\). We observe that two SN bifurcation points get closer as the value of \(b_{\mathrm{target}}\) decreases.

Now we present a close look at the iterative process for the case of \(b_{\mathrm{target}}\) = 0.1. The evolution of the objective function \(\rho_{\max}\) in optimization formulation 24 is shown in the left panel of Fig. 10. It can be observed that the value of the objective function increases during the initial iterations and gradually stabilizes as the optimization progresses. The corresponding evolution of the FRCs is shown in the right panel of Fig.10, where the peak amplitude of the FRCs also increases with the number of iterations. The computational performance of the above iterative process is summarized in Table 4.

Similar to the previous two examples, it is necessary to examine the convergence of the SSM-based model reduction. In Appendix 10, we compute the FRC of the optimized structure shown in Fig. 9 using a seventh-order SSM-reduced model. The results indicate that the third-order SSM model has already achieved convergence.

6 Conclusion↩︎

In this paper, we have performed topology optimization to tune the forced response curves (FRCs) of high-dimensional nonlinear systems. By employing reduction on spectral submanifolds (SSMs), we have enabled efficient analysis and sensitivity computation of FRCs for large-scale finite element models. Several optimization objectives were considered, including minimizing the FRC peak, tailoring hardening/softening behavior, and tuning saddle-node (SN) bifurcations on FRC.

To tune SN points on the FRC, we have derived an explicit expression for the coefficient governing the distance between two SN bifurcation points, along with its sensitivity, which allows for the effective control of jump phenomena in nonlinear MEMS devices. Numerical examples demonstrate the effectiveness of the proposed approach in designing structures with desired dynamic characteristics. This framework offers a practical approach for incorporating nonlinear dynamic behavior into topology optimization, with potential applications in devices such as MEMS sensors and energy harvesters.

In future work, we aim to extend this framework to the control of internal resonances, which often arise in nonlinear systems and significantly influence complex dynamic responses. In particular, internal resonances may induce Hopf bifurcations, which lead to quasi-periodic motions on invariant tori [34], [40]. We plan to explore how such modal interactions can be systematically controlled by manipulating reduced-order dynamics, including the regulation of Hopf bifurcations and the resulting quasi-periodic responses, to enable advanced control of multi-frequency behaviors in nonlinear structures.

Acknowledgements↩︎

HL and ML acknowledge the financial support of the National Natural Science Foundation of China (No. 12302014) and State Key Laboratory of Structural Analysis, Optimization and CAE Software for Industrial Equipment (GZ24117). MP and JM acknowledge the financial support of STMicroelectronics (award number 4000614871).

Data availability↩︎

The data used to generate the numerical results included in this paper are available from the corresponding author on request.

Conflict of interest↩︎

The authors declare that they have no conflict of interest.

7 The sensitivity of \(\tilde{f}\)↩︎

We begin by listing all relevant constraints \[\begin{align} \tilde{f}-0.5\boldsymbol{\psi}^\mathrm{T}\boldsymbol{f}^{\mathrm{ext}}=0, \\ \boldsymbol{K}\boldsymbol{\phi}=\omega^2\boldsymbol{M}\boldsymbol{\phi}, \boldsymbol{\phi}^\mathrm{T}\boldsymbol{M}\boldsymbol{\phi}=1, \\ \boldsymbol{\psi}=\kappa\boldsymbol{\phi}, \\ \kappa \omega=\frac{-\mathrm{i}}{2\sqrt{1-\xi^2}}, \\ 2\xi\omega=\alpha+\beta\omega^2. \end{align}\] Now we define a Lagrangian as \[\begin{align} L=&\tilde{f} + \eta_{\tilde{f}}(\tilde{f}-0.5\boldsymbol{\psi}^\mathrm{T}\boldsymbol{f}^{\mathrm{ext}}) +\eta_{\boldsymbol{\phi}}(\boldsymbol{K}\boldsymbol{\phi}-\omega^2\boldsymbol{M}\boldsymbol{\phi})+\eta_{\mathrm{norm}}(\boldsymbol{\phi}^\mathrm{T}\boldsymbol{M}\boldsymbol{\phi}-1)\\ &+\eta_{\boldsymbol{\psi}}(\boldsymbol{\psi}-\kappa\boldsymbol{\phi}) +\eta_{\kappa}\bigl(\kappa \omega+\tfrac{\mathrm{i}}{2\sqrt{1-\xi^2}}\bigr) +\eta_{\xi}\bigl(2\xi\omega-(\alpha+\beta\omega^2)\bigr). \end{align}\] The variation of the Lagrangian is given by \[\begin{align} \delta L =& \delta \eta_{\tilde{f}} (\tilde{f} - 0.5\boldsymbol{\psi}^{\mathrm{T}}\boldsymbol{f}^{\mathrm{ext}}) + \delta \boldsymbol{\eta}_{\boldsymbol{\phi}}^{\mathrm{T}} (\boldsymbol{K}\boldsymbol{\phi} - \omega^2\boldsymbol{M}\boldsymbol{\phi}) \\ &+ \delta \eta_{\mathrm{norm}} (\boldsymbol{\phi}^{\mathrm{T}}\boldsymbol{M}\boldsymbol{\phi} - 1) + \delta \boldsymbol{\eta}_{\boldsymbol{\psi}}^{\mathrm{T}} (\boldsymbol{\psi} - \kappa\boldsymbol{\phi}) \\ &+ \delta \eta_{\kappa} \bigl(\kappa\omega + \tfrac{\mathrm{i}}{2\sqrt{1-\xi^2}}\bigr) + \delta \eta_{\xi} \bigl(2\xi\omega - (\alpha + \beta\omega^2)\bigr) \\ &+ \bigl(1 + \eta_{\tilde{f}}\bigr)\delta \tilde{f} + \bigl(\boldsymbol{\eta}_{\boldsymbol{\psi}}^{\mathrm{T}} -0.5\eta_{\tilde{f}}{\boldsymbol{f}^{\mathrm{ext}}}^\mathrm{T}\bigr) \delta \boldsymbol{\psi} \\ &+ \boldsymbol{\eta}_{\boldsymbol{\phi}}^{\mathrm{T}} \bigl(\delta\boldsymbol{K}\boldsymbol{\phi} - \omega^2 \delta \boldsymbol{M} \boldsymbol{\phi}\bigr) \\ &+ \bigl( \boldsymbol{\eta}_{\boldsymbol{\phi}}^{\mathrm{T}}(\boldsymbol{K} - \omega^2 \boldsymbol{M}) + 2 \eta_{\mathrm{norm}} \boldsymbol{\phi}^{\mathrm{T}} \boldsymbol{M} - \kappa \boldsymbol{\eta}_{\boldsymbol{\psi}}^{\mathrm{T}} \bigr) \delta \boldsymbol{\phi} \\ &+ \bigl(\eta_{\kappa} \kappa - 2\omega\boldsymbol{\eta}_{\boldsymbol{\phi}}^{\mathrm{T}}\boldsymbol{M}\boldsymbol{\phi} + 2\xi\eta_{\xi} - 2\eta_{\xi} \beta \omega \bigr) \delta \omega \\ &+ \eta_{\mathrm{norm}}\boldsymbol{\phi}^{\mathrm{T}}\delta \boldsymbol{M}\boldsymbol{\phi} + \bigl(\eta_{\kappa}\omega - \boldsymbol{\eta}_{\boldsymbol{\psi}}^{\mathrm{T}}\boldsymbol{\phi}\bigr)\delta \kappa \\ &+ \bigl(\eta_{\kappa}\tfrac{\xi\mathrm{i}}{2\sqrt{(1-\xi^2)^3}} + 2\eta_{\xi}\omega \bigr)\delta \xi. \end{align}\] By collecting the coefficients of the independent variations, we obtain the following adjoint equations \[\begin{align} \delta \tilde{f}:&\quad 1 + \eta_{\tilde{f}} = 0,\\ \delta \boldsymbol{\psi}:&\quad \boldsymbol{\eta}_{\boldsymbol{\psi}} - 0.5\eta_{\tilde{f}}{\boldsymbol{f}^\mathrm{ext}} = 0,\\ \delta \boldsymbol{\phi}:&\quad \bigl(\boldsymbol{K} - \omega^2 \boldsymbol{M}\bigr)\boldsymbol{\eta}_{\boldsymbol{\phi}} + 2\eta_{\mathrm{norm}}\boldsymbol{M}\boldsymbol{\phi} - \kappa \boldsymbol{\eta}_{\boldsymbol{\psi}} = 0,\\ \delta \omega:&\quad \eta_\kappa\,\kappa - 2\omega\boldsymbol{\eta}_{\boldsymbol{\phi}}^{\mathrm{T}}\boldsymbol{M}\boldsymbol{\phi} + 2\eta_\xi(\xi - \beta\omega) = 0,\\ \delta \kappa:&\quad \eta_\kappa\,\omega - \boldsymbol{\eta}_{\boldsymbol{\psi}}^{\mathrm{T}}\boldsymbol{\phi} = 0,\\ \delta \xi:&\quad \eta_{\kappa}\tfrac{\xi\mathrm{i}}{2\sqrt{(1-\xi^2)^3}} + 2\eta_{\xi}\omega = 0. \end{align}\] Thus, we have \(\eta_{\tilde{f}}=-1\), \(\boldsymbol{\eta}_{\boldsymbol{\psi}}=0.5\eta_{\tilde{f}}\boldsymbol{f}^{\mathrm{ext}}\), \(\eta_{\kappa}={\boldsymbol{\eta}_{\boldsymbol{\psi}}}^{\mathrm{T}}\boldsymbol{\phi}/\omega\), \(\eta_{\xi}=-\frac{\eta_{\kappa}\xi\mathrm{i}}{4 \omega \sqrt{(1-\xi^2)^3}}\). Let \(\boldsymbol{b}_{\mathrm{norm}}=\kappa \boldsymbol{\eta}_{\boldsymbol{\psi}}\), \(b_{\boldsymbol{\phi}}=\eta_\kappa \kappa + 2 \eta_{\xi}(\xi-\beta \omega)\), \(\boldsymbol{\eta}_{\boldsymbol{\phi}}\) and \(\eta_{\mathrm{norm}}\) can be obtained by \[\begin{align} \begin{pmatrix} 2 \omega \boldsymbol{\phi}^{\mathrm{T}}\boldsymbol{M} & 0 \\ \boldsymbol{K}-\omega^2 \boldsymbol{M} & 2 \boldsymbol{M} \boldsymbol{\phi} \end{pmatrix} \begin{pmatrix} \boldsymbol{\eta}_{\boldsymbol{\phi}} \\ \eta_{\mathrm{norm}} \end{pmatrix} = \begin{pmatrix} b_{\boldsymbol{\phi}} \\ \boldsymbol{b}_{\mathrm{norm}} \end{pmatrix}. \end{align} \label{eq:solve95lambda}\tag{25}\] Thus, we obtain the following expression for the sensitivity of \(\tilde{f}\) \[\begin{align} \delta \tilde{f} = \boldsymbol{\eta}_{\boldsymbol{\phi}}^{\mathrm{T}} \bigl(\delta\boldsymbol{K}\boldsymbol{\phi} - \omega^2 \delta \boldsymbol{M}\boldsymbol{\phi}\bigr) + \eta_{\mathrm{norm}} \boldsymbol{\phi}^{\mathrm{T}} \delta \boldsymbol{M} \boldsymbol{\phi}. \end{align}\]

8 The sensitivity of \(b\)↩︎

We first compute the sensitivity of \(\rho_{\mathrm{SN}}\) and \(\Omega_{\rm{SN}}\). Similar to compute \(\rho_{\max}\), we compute the sensitivity of \(\rho_{\rm{SN}}\) and \(\Omega_{\rm{SN}}\) by taking derivative of 20 , namely \[\begin{align} \rho_{\mathrm{SN}}' = \frac{c_1}{c_2}, \Omega_{\mathrm{SN}}' = c_3 -\frac{c_4}{c_5}, \end{align}\] where \[\begin{align} c_1 = & 8\mathrm{Im}(\gamma)\mathrm{Im}(\gamma)'\rho^6(\epsilon^2|\tilde{f}|^2 - (\mathrm{Re}(\lambda) \rho + \mathrm{Re}(\gamma) \rho^3)^2) \\ & + 4\mathrm{Im}(\gamma)^2\rho^6(2\epsilon^2|\tilde{f}||\tilde{f}|' - 2(\mathrm{Re}(\lambda)\rho + \mathrm{Re}(\gamma)\rho^3)(\mathrm{Re}(\lambda)'\rho + \mathrm{Re}(\gamma)'\rho^3)) \\ & - 2(\epsilon^2|\tilde{f}|^2 + 2 \mathrm{Re}(\lambda) \mathrm{Re}(\gamma) \rho^4 + 2\mathrm{Re}(\gamma)^2 \rho^6) \\ & (2\epsilon^2 |\tilde{f}||\tilde{f}|' + 2\mathrm{Re}(\lambda)'\mathrm{Re}(\gamma)\rho^4 + 2\mathrm{Re}(\lambda)\mathrm{Re}(\gamma)'\rho^4 + 4\mathrm{Re}(\gamma)\mathrm{Re}(\gamma)'\rho^6), \\ c_2 =& 2(\epsilon^2|\tilde{f}|^2 + 2\mathrm{Re}(\lambda)\mathrm{Re}(\gamma)\rho^4 + 2\mathrm{Re}(\gamma)^2\rho^4 + 2\mathrm{Re}(\gamma)^2\rho^6) \\ &(8\mathrm{Re}(\lambda)\mathrm{Re}(\gamma)\rho^3 + 12\mathrm{Re}(\gamma)^2\rho^5) \\ & + 24\mathrm{Im}(\gamma)^2\rho^5((\mathrm{Re}(\lambda)\rho + \mathrm{Re}(\gamma)\rho^3)^2 - \epsilon^2|\tilde{f}|^2) \\ & + 8\mathrm{Im}(\gamma)^2(\mathrm{Re}(\lambda)\rho + \mathrm{Re}(\gamma)\rho^3)(\mathrm{Re}(\lambda) +3\mathrm{Re}(\gamma)\rho^2)\rho^6, \\ c_3 = &\mathrm{Im}(\lambda)' + \mathrm{Im}(\gamma)'\rho^2 + 2\mathrm{Im}(\gamma)\rho\rho', \\ c_4 = & 2\epsilon^2|\tilde{f}||\tilde{f}|' - 2(\mathrm{Re}(\lambda)\rho + \mathrm{Re}(\gamma)\rho^3) \\ &(\mathrm{Re}(\lambda)'\rho + \mathrm{Re}(\lambda)\rho' + \mathrm{Re}(\gamma)'\rho^3+3\mathrm{Re}(\gamma)\rho^2\rho') \\ &-2(\mathrm{Im}(\lambda) - \Omega_{\mathrm{SN}} + \mathrm{Im}(\gamma)\rho^2)^2\rho\rho', \\ c_5 = &2(\mathrm{Im}(\lambda)-\Omega+\mathrm{Im}(\gamma)\rho^2)\rho^2. \end{align}\] Then we consider the sensitivity of \(\boldsymbol{A}_{\mathrm{SN}}\) and \(\boldsymbol{B}_{\mathrm{SN}}\). By using the chain rule to the expressions in 22 , we obtain \[\begin{align} \boldsymbol{A}_{\mathrm{SN}}' &= \partial_{\boldsymbol{\mu}} \boldsymbol{A}_{\mathrm{SN}} + \partial_{\Omega_{\mathrm{SN}}} \boldsymbol{A}_{\mathrm{SN}} \Omega_{\mathrm{SN}}' + \partial_{\rho_{\mathrm{SN}}} \boldsymbol{A}_{\mathrm{SN}} \rho_{\mathrm{SN}}', \\ \boldsymbol{B}_1' &= \partial_{\boldsymbol{\mu}} \boldsymbol{B}_1 + \partial_{\Omega_{\mathrm{SN}}} \boldsymbol{B}_1 \Omega_{\mathrm{SN}}' + \partial_{\rho_{\mathrm{SN}}} \boldsymbol{B}_1 \rho_{\mathrm{SN}}', \\ \boldsymbol{B}_2' &= \partial_{\boldsymbol{\mu}} \boldsymbol{B}_2 + \partial_{\Omega_{\mathrm{SN}}} \boldsymbol{B}_2 \Omega_{\mathrm{SN}}' + \partial_{\rho_{\mathrm{SN}}} \boldsymbol{B}_2 \rho_{\mathrm{SN}}', \\ \end{align}\] where \[\begin{align} \partial_{\boldsymbol{\mu}} \boldsymbol{A} &= \begin{pmatrix} \mathrm{Re}(\lambda)' + 3\mathrm{Re}(\gamma)' \rho_{\mathrm{SN}}^2 & -\left(\mathrm{Im}(\lambda)' + \mathrm{Im}(\gamma)' \rho_{\mathrm{SN}}^2\right) \rho_{\mathrm{SN}} \\ 2\mathrm{Im}(\gamma)' \rho_{\mathrm{SN}} + \left(\mathrm{Im}(\lambda)' + \mathrm{Im}(\gamma)' \rho_{\mathrm{SN}}^2 \right)/\rho_{\mathrm{SN}} & \mathrm{Re}(\lambda)' + \mathrm{Re}(\gamma)' \rho_{\mathrm{SN}}^2 \end{pmatrix}, \\ \partial_{\Omega_{\mathrm{SN}}} \boldsymbol{A} &= \begin{pmatrix} 0 & \rho_{\mathrm{SN}} \\ -1/\rho_{\mathrm{SN}} & 0 \end{pmatrix}, \\ \partial_{\rho_{\mathrm{SN}}} \boldsymbol{A} &= \begin{pmatrix} 6\mathrm{Re}(\gamma)\rho_{\mathrm{SN}} & -\mathrm{Im}(\lambda) + \Omega_{\mathrm{SN}} - 3\mathrm{Im}\gamma\rho_{\mathrm{SN}}^2 \\ 2 \mathrm{Im}(\gamma) - (\mathrm{Im}(\lambda) - \Omega_{\mathrm{SN}})/\rho_{\mathrm{SN}}^2 + \mathrm{Im}(\gamma) & 2\mathrm{Re}(\gamma)\rho_{\mathrm{SN}} \end{pmatrix}, \\ \partial_{\boldsymbol{\mu}} \boldsymbol{B}_1 &= \begin{pmatrix} 6\mathrm{Re}(\gamma)'\rho_{\mathrm{SN}} & 0 \\ 0 & \mathrm{Re}(\lambda)'\rho_{\mathrm{SN}} + \mathrm{Re}(\gamma)' \rho_{\mathrm{SN}}^3 \end{pmatrix}, \\ \partial_{\Omega_{\mathrm{SN}}} \boldsymbol{B}_1 &= \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix}, \\ \partial_{\rho_{\mathrm{SN}}} \boldsymbol{B}_1 &= \begin{pmatrix} 6\mathrm{Re}(\gamma) & 0 \\ 0 & \mathrm{Re}(\lambda)+3\mathrm{Re}(\gamma) \rho_{\mathrm{SN}}^2 \end{pmatrix}, \\ \partial_{\boldsymbol{\mu}} \boldsymbol{B}_2 &= \begin{pmatrix} 2\mathrm{Im}(\gamma)' - 2(\mathrm{Im}(\lambda)'+\mathrm{Im}(\gamma)'\rho_{\mathrm{SN}}^2)/\rho_{\mathrm{SN}}^2 & -(\mathrm{Re}(\lambda)'+\mathrm{Re}(\gamma)'\rho_{\mathrm{SN}}^2)/\rho_{\mathrm{SN}} \\ -(\mathrm{Re}(\lambda)'+\mathrm{Re}(\gamma)'\rho_{\mathrm{SN}}^2)/\rho_{\mathrm{SN}} & \mathrm{Im}(\lambda)'+\mathrm{Im}(\gamma)' \rho_{\mathrm{SN}}^2 \end{pmatrix}, \\ \partial_{\Omega_{\mathrm{SN}}} \boldsymbol{B}_2 &= \begin{pmatrix} 2/\rho_{\mathrm{SN}}^2 & 0 \\ 0 & -1 \end{pmatrix}, \\ \partial_{\rho_{\mathrm{SN}}} \boldsymbol{B}_2 &= \begin{pmatrix} 4(\mathrm{Im}(\lambda)-\Omega_{\mathrm{SN}})/\rho_{\mathrm{SN}}^3 & \mathrm{Re}(\lambda)/\rho_{\mathrm{SN}}^2-\mathrm{Re}(\gamma) \\ \mathrm{Re}(\lambda)/\rho_{\mathrm{SN}}^2-\mathrm{Re}(\gamma) & 2\mathrm{Im}(\gamma) \rho_{\mathrm{SN}} \end{pmatrix}. \end{align}\] Thus, the sensitivity of b with respect to the design variable \(\boldsymbol{\mu}\) is given by \[\begin{align} b'=\boldsymbol{\psi}_{\mathrm{SN}}'\boldsymbol{B}_{\mathrm{SN}}(\boldsymbol{\phi}_{\mathrm{SN}},\boldsymbol{\phi}_{\mathrm{SN}}) + \boldsymbol{\psi}_{\mathrm{SN}}\boldsymbol{B}_{\mathrm{SN}}'(\boldsymbol{\phi}_{\mathrm{SN}},\boldsymbol{\phi}_{\mathrm{SN}})+ 2\boldsymbol{B}_{\mathrm{SN}}(\boldsymbol{\phi}_{\mathrm{SN}},\boldsymbol{\phi}_{\mathrm{SN}}'), \end{align}\] where \(\boldsymbol{\psi}'_{\mathrm{SN}}\) and \(\boldsymbol{\phi}'_{\mathrm{SN}}\) are obtained by \[\begin{align} \begin{pmatrix} \boldsymbol{A}^{\mathrm{T}} & \boldsymbol{\psi}_{\mathrm{SN}} \\ \boldsymbol{\psi}_{\mathrm{SN}}^{\mathrm{T}} & 0 \end{pmatrix}\begin{pmatrix} \boldsymbol{\psi}'_{\mathrm{SN}} \\ \xi \end{pmatrix}&= \begin{pmatrix} -\boldsymbol{A}'^{\mathrm{T}}\boldsymbol{\psi}_{\mathrm{SN}} \\ 0 \end{pmatrix}, \\ \begin{pmatrix} \boldsymbol{A} & \boldsymbol{\psi}_{\mathrm{SN}} \\ \boldsymbol{\psi}_{\mathrm{SN}}^{\mathrm{T}} & 0 \end{pmatrix}\begin{pmatrix} \boldsymbol{\phi}'_{\mathrm{SN}} \\ \xi \end{pmatrix}&= \begin{pmatrix} -\boldsymbol{A}'\boldsymbol{\phi}_{\mathrm{SN}} \\ -\boldsymbol{\phi}_{\mathrm{SN}}^{\mathrm{T}} \boldsymbol{\psi}_{\mathrm{SN}}' \end{pmatrix}. \end{align}\]

9 The convergence validation of SSM-based reduction in Sec. 5.1↩︎

To assess the convergence of the SSM-based model reduction, we compute the FRCs of the optimized structures using seventh-order SSM reduction and compare them with those obtained from third-order SSM reduction in Fig. 4. As shown in Fig. 11, the FRCs for both cases exhibit excellent agreement between the two SSM orders under a forcing amplitude of \(f^{\mathrm{ext}} = 2 \times 10^{10} \;\mathrm{ng \cdot \mu m/ms^2}\). This confirms that the third-order SSM reduction is sufficiently accurate for capturing the nonlinear dynamic behavior of the optimized structures in this study.

10 The convergence validation of SSM-based reduction in Sec. 5.3↩︎

To validate the convergence of the SSM-based reduction in the SN bifurcation control example presented in Sec. 5.3, we compute the FRCs of the optimized structures using both third-order and seventh-order reduced models. As shown in Fig.12, consistent agreement is observed across all three target values of \(b_{\mathrm{target}}\) (2, 1, and 0.1), indicating that the third-order SSM approximation provides sufficient accuracy for capturing the key nonlinear response features in this study.

References↩︎