An Efficient and Robust Topology Optimization Method for Thermoelastically Damped Microresonators

Abstract

The challenges of computational cost and robustness are critical obstacles in topology optimization methods, particularly for the iterative process of optimizing large-scale multiphysical structures. This study proposes an efficient and robust topology optimization method for minimizing the thermoelastic damping of large-scale microresonators. An evolutionary structural optimization method is adopted to passively determine the search direction of optimizing large-scale thermoelastic structures. To efficiently reduce the computational cost of the iterative process of an optimizing process, a model reduction method is developed based on the projection-based model reduction method whose reduced basis is generated within the Neumann series subspace. However, the projection-based model reduction method may be unstable when topology modifications are made during an iteration optimization process. To ensure robustness, a modal validation technique is first implemented in the iterative process to stabilize the iteration and narrow down the search domain, and a posterior evaluation of the Neumann series expansion is then developed to retain the convergence of the projection-based model reduction method. Furthermore, the efficiency and accuracy of the proposed topology optimization method are validated through numerical examples. Two large-scale numerical models are also used to demonstrate the advantage of the proposed method. It is found that large-scale thermoelastic structures with a phase-lag heat conduction law can be designed passively, precisely, and efficiently by using the proposed topology optimization method.

1. Introduction

Micro-electromechanical system (MEMS) resonators are currently dominating the billion-dollar market, particularly in the fields of biomedicine, aerospace, and others [1,2]. High-frequency resonators are widely used as the key components of MEMS. Two critical dynamic characteristics for microresonators are known as resonance frequency and energy dissipation (inverse quality factor). Frequency affects the sensitivity of the response of resonators, while the energy dissipation rate affects their accuracy and noise. Enlarging the quality factor of microresonators can, therefore, reduce their signal loss and improve stability [3]. The attainment of both high sensitivity and high accuracy is a vital requirement for most MEMS resonators; unfortunately, resonance frequency and energy dissipation are generally mutually exclusive.

Although the quest continues for higher frequency resonators, these have little-to-no use as MEMS resonators without appropriate accuracy (qualify factor). One of the major damping sources of microresonators is thermoelastic damping [4], which cannot be reduced by external measurement-and-control methods. Thermoelastic damping determines the upper limit of the performance of high-frequency MEMS resonators [5,6], although their resonance frequency is generally high [7]. For these reasons, the development of high-frequency yet low-dissipation microresonators is an optimization problem for an exercise in the compromise between frequency and quality.

As depicted in Figure 1, a thermoelastic structure undergoes vibration, resulting in the emergence of compressed and stretched regions within the structure. Since temperature increases in the compressed region but decreases in the stretched region, the temperature gradient-derived heat conduction will lead to irreversible thermoelastic damping. It has been examined by some works [8,9] that microstructural features can affect the thermoelastic damping of microresonators, and structural topology design methods can be used to enlarge the quality factor of microresonators. On the other hand, advancements in lithography manufacturing processes and 3D printing technology have gradually loosened the constraints of manufacturing processes on the design of thermoelastic structures. As a result, it has become possible to fabricate microresonators with intricate internal configurations [10].

The fundamental theory of the structural topology optimization method is defined by combining homogenization theory with numerical optimization methods. This method achieves optimal structural performance by modifying the micro-structural features within the design domain [11]. There have been extensive studies in this field, including the classic solid isotropic material with the penalization (SIMP)-based optimization criteria method [12], the evolutionary and bi-directional evolutionary structural optimization [13], the Smooth-edged material distribution for optimizing topology algorithm [14], and the floating projection topology optimization [15].

The heat transfer mechanism plays a significant role in the analysis and design of thermoelastic structures. It has been demonstrated through a large number of experiments that the assumption of infinite heat conduction velocity, which is required for linear Fourier heat conduction, is no longer valid for devices at the micron or even nanoscale. Instead, nonlinear non-Fourier heat conduction gradually dominates as the size of thermoelastic structures is down from the micrometer to the nanometer [16]. For instance, extreme temperatures in a micro-cantilever not only affect its Young’s modulus but also change the tension in the cantilever due to differential thermal expansion, resulting in the variation of both frequency and dissipation [17]. To better characterize these nonlinear heat conduction effects, scholars have extensively studied the phenomenon since the classical Maxwell–Cattaneo–Vernotte equation [18,19]. Among the non-Fourier heat conduction models, there is a specific branch dedicated to modeling non-Fourier phase-lag behavior. This branch attributes nonlinear non-Fourier heat conduction effects to the hysteresis effect between heat conduction and temperature gradients. One of the features of the non-Fourier phase-lag model is its ability to be degenerated by setting specific lag values [20]. According to the heat wave theory, the lag of heat conduction in solid structures usually affects the heat flux. In other words, the change of heat flux is driven by the changing temperature gradient [21]. Additionally, the single-phase lag and the dual-phase lag can be converted to each other [22]. For these reasons, we choose the dual-phase lag heat conduction model as the non-Fourier heat conduction law in this paper due to its universal nature.

The thermoelastic structures with the dual-phase lag heat conduction model will cause a nonlinear thermoelastic eigenproblem, which is very different to handle. Often, the linearization technique based on the state–space method [23,24] is used for solving the nonlinear eigenproblem. When the parametric partial differential equations of thermoelastic structures are discretized using commonly-used numerical techniques, the degrees of freedom of the governing equations will be very large. However, the dimension of the system matrix of the state–space method is much larger than the original dimension of the nonlinear thermoelastic eigenproblem, thereby resulting in a huge computational challenge.

There has been extensive research on accelerating the solution of large-scale linear eigenproblems, which can, of course, be used for handling the state–space method of thermoelastic eigenproblems with the dual-phase lag heat conduction model. The reduced basis method, which finds the relationship between the design variable and the system information (such as the displacement field) to solve the problem, has been employed to reduce the dimensionality of the problem and computational cost [25]. Furthermore, the matrix operator method, when combined with other optimization algorithms, has been utilized to improve stability and accelerate the solution of partial differential equations [26], and has been widely applied to solve finite element eigenvalue problems in multi-field coupled dynamical systems [27]. Currently, there are many types of matrix operators, including diagonal matrix operators, reduced-basis matrix operators, and differential matrix operators, and these operators are often necessary for subspace solutions of partial differential equations [28,29]. Among the types of matrix operators, diagonal and triangular matrix operators are known for their simplicity and fast convergence [30]. These efficient methods, although general, do not capture the essential features of thermoelastic structures. For the sake of further computational efficiency, Fu et al. [31] constructed diagonal systemic matrices of the lumped mass and heat generation matrices. With these diagonal systemic matrices, an efficient Krylov subspace method is proposed for the thermoelastic eigenproblem. Note that the analysis of thermoelastic eigenproblems with the Fourier heat conduction model has been provided recently for commercial software COMSOL. However, these works mentioned above can only be used for thermoelastic eigenproblems with the Fourier heat conduction model, rather than the non-Fourier heat conduction model of interest. Furthermore, the repeated updating of topology optimization designs is computationally challenging, especially for large-scale linear eigenproblems of thermoelastic eigenproblems with the dual-phase lag heat conduction model.

To the literature above, current system solutions of large-scale thermoelastic structures are difficult and time-consuming. Furthermore, the iterative optimization of thermoelastics requires an eigenvalue analysis in each iteration; as a result, the consumption of computational time and resources could be unbearable. Although an efficient design method of thermoelastic structures by topology optimization is proposed in our previous research [31], it could not resolve thermoelastic problems in the framework of the phase-lag heat conduction law and may also suffer from an iterative procedure oscillation. Thus, an effective and robust design method is urgently needed for fully-coupled thermoelastic microstructures under non-Fourier phase-lag heat conduction law. To address the need for an efficient design method for thermoelastic microstructures with phase-lag heat conduction law, and considering that existing design methods are not qualified to solve large-scale problems with complex geometry, this paper proposes an efficient and robust topology optimization method for minimizing the thermoelastic damping of large-scale microresonators with phase-lag heat conduction law.

This article is organized as follows. First, a state–space method is developed for large-scale microresonators with phase-lag heat conduction law in the framework of commonly-used numerical techniques, and the corresponding optimization method is defined to enlarge the quality factor of the large-scale microresonators. Then, to efficiently reduce the computational cost of the iterative process of optimizing process, a model reduction method is developed based on the projection-based model reduction method whose reduced basis is generated based on the Neumann series expansion technique. However, the model reduction method may be unstable when topology modifications are made during an iteration optimization process. Thus, to ensure robustness, a modal validation technique is first implemented in the iterative process to stabilize the iteration and narrow down the search domain, and a posterior evaluation of the Neumann series expansion is then developed to retain the convergence of the parameterized model reduction method. Finally, several numerical results are discussed to validate the efficiency and accuracy of the proposed method.

2. Statement of the Topology Optimization Problem

The objective of this topology optimization method is to enhance the quality factor of thermoelastic microresonators with the phase-lag heat conduction law. As per the definition of the phase-lag heat conduction model, the dissipation resulting from the non-Fourier form of the entropy increases because of the irreversible process related to both internal energy and heat flux and can be converted into a time-dependent lag form [32]. The introduction of higher-order time-dependent terms significantly increases the complexity of solving the system equations.

In an optimal design method, the search direction, search domain, efficiency, and convergence are the crucial properties [33]. The search direction determines the ability to obtain optimal solutions to a design problem, while the search domain influences the speed of the method. Moreover, the optimization method must be convergent to obtain an appropriate solution.

In this section, we propose a modified evolutionary structural optimization method to determine the search direction of the design method.

2.1. Thermoelastic Damping with Phase Lag Heat Conduction Law

In the case of a small energy dissipation (thermoelastic damping) of the microresonator, the quality factor of the i-th mode shape, $Q_i$, is similar to the dissipation angle [34], which can be obtained by the following formula [35]:

$$\begin{array}{c}\hfill Q_i^{-1}=2\left|\frac{\Im \left({\lambda }_i\right)}{\Re \left({\lambda }_i\right)}\right|,\end{array}$$

(1)

where ${\lambda }_i$ is the i-th eigenvalue of the system, and $\Im (•)$ and $\Re (•)$ denote the imaginary and real part of $(•)$, respectively. As the system eigenvalues have to be solved by the governing equations of the thermoelastic system, the fully-coupled equation will be derived through the law of momentum and heat conduction.

2.2. Governing Equations of Thermoelstic Structures

For a thermoelastic system constructed with linear isotropic thermoelastic material, the equation of motion can be obtained through the conservation law of momentum [36]:

$$\begin{array}{c}\hfill \frac{E}{1+\nu }\mathsf{ϵ}·\nabla +\frac{E\nu }{(1+\nu )(1-2\nu )}\nabla \left[\mathrm{Tr}\left(\mathsf{ϵ}\right)\right]-\frac{E\alpha }{1-2\nu }\nabla \theta +\tilde{\mathbf{b}}=\rho \ddot{\mathbf{u}}\end{array}$$

(2)

where $\theta$ is the temperature field, E is the elastic modulus, $\alpha$ is the coefficient of thermal expansion, $T_0$ is the reference temperature, and $\nu$ is Poisson’s ratio.

This paper mainly considers the dynamic characteristics of structures, i.e., free vibration analysis (modal analysis), so the body force $\tilde{\mathbf{b}}$, external force, and corresponding load terms are zero. Combining the finite element method and classical global matrix assembly method, and considering boundary conditions, the global equation of motion can be expressed as

$$\begin{array}{c}\hfill \mathbf{M}\ddot{\mathbf{U}}+\mathbf{KU}-\mathbf{G}\mathsf{\Theta }=\mathbf{0}\end{array}$$

(3)

where $\mathbf{M}$, $\mathbf{K}$, and $\mathbf{G}$ are the global mass, stiffness, and thermal stress matrices of the whole structure, $\mathbf{U}$ and $\mathsf{\Theta }$ are the global vectors of displacement and temperature, respectively. Note that the thermal stress term ($\mathbf{G}\mathsf{\Theta }$) is related to the temperature field of the structure. Therefore, the solution of the structural equation of motion requires the use of the structural heat conduction equation.

On the other hand, the dual-phase lag heat conduction model takes the following form [20,22]:

$$\begin{array}{c}\hfill \mathbf{q}(\mathbf{r},t)+{\tau }_q\dot{\mathbf{q}}(\mathbf{r},t)=-\kappa \left[\nabla T(\mathbf{r},t)+{\tau }_T\nabla \dot{T}(\mathbf{r},t)\right],\end{array}$$

(4)

where $\mathbf{q}$ and $\nabla T$ are the heat flux and temperature gradient, respectively. $\mathbf{r}$ and t are, respectively, the position vector and time. ${\tau }_q$ and ${\tau }_T$ are the phase lags on heat flux and temperature, respectively. Overdot denotes the derivative with respect to time.

The following equation describes the behavior of heat conduction and is known as the governing equation of heat conduction:

$$\begin{array}{c}\hfill \kappa \Delta \theta =C_v\dot{\theta }+\frac{E\alpha T_0}{1-2\nu }\mathrm{Tr}\left(\dot{\mathsf{ϵ}}\right)\end{array}$$

(5)

where $C_v=\rho C_p$ represents the specific heat capacity per unit volume, $\Delta$ represents the Laplacian operator, $\kappa$ is the thermal diffusivity, $\dot{\theta }$ is the rate of change of the temperature field, and $\dot{\mathsf{ϵ}}$ is the strain rate tensor.

Similar to the thermal stress term in the governing equation, the heat generation term ($\mathbf{F}\dot{\mathbf{U}}$) caused by strain in the thermal conduction equation is related to the structural displacement field, and therefore, the solution of the structural thermal conduction equation requires the use of the motion equation.

Combining Equations (4) and (5), the global equation of heat conduction considering the non-Fourier heat conduction dual-phase lag model can be expressed as [36]:

$$\begin{array}{c}\hfill {\mathbf{S}}_{}{\mathsf{\Theta }}_{}+{{\mathbf{S}}_{}}_T{\dot{\mathsf{\Theta }}}_{}+{\mathbf{H}}_{}{\dot{\mathsf{\Theta }}}_{}+{\mathbf{F}}_{}{\dot{\mathbf{U}}}_{}+{{\mathbf{H}}_{}}_q{\ddot{\mathsf{\Theta }}}_{}+{{\mathbf{F}}_{}}_q{\ddot{\mathbf{U}}}_{}=\mathbf{0}\end{array}$$

(6)

where

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\mathbf{S}}_T={\tau }_T{\mathbf{S}}_{}\hfill \\ {\mathbf{H}}_q={\tau }_q{\mathbf{H}}_{}\hfill \\ {\mathbf{F}}_q={\tau }_q{\mathbf{F}}_{}\hfill \end{array}\right.\end{array}$$

(7)

The governing heat conduction equation can take various forms depending on the chosen phase lag. In single-phase lag non-Fourier heat conduction models, the phase lag is usually applied to the heat flux [21]. Therefore, only the heat flux lag form is considered. When the temperature lag ${\tau }_T$ is equal to zero, Equation (6) degenerates to a single-phase lag form.

$$\begin{array}{c}\hfill {\mathbf{S}}_{}{\mathsf{\Theta }}_{}+{\mathbf{H}}_{}{\dot{\mathsf{\Theta }}}_{}+{\mathbf{F}}_{}{\dot{\mathbf{U}}}_{}+{\mathbf{H}}_{}{\ddot{\mathsf{\Theta }}}_{}+{\mathbf{F}}_q{\ddot{\mathbf{U}}}_{}=\mathbf{0}\end{array}$$

(8)

When both the temperature and heat flux lags are zero (${\tau }_q=0,{\tau }_T=0$), Equation (6) reduces to the classical Fourier heat conduction form [31]:

$$\begin{array}{c}\hfill \mathbf{S}\mathsf{\Theta }+\mathbf{H}\dot{\mathsf{\Theta }}+\mathbf{F}\dot{\mathbf{U}}=\mathbf{0}\end{array}$$

(9)

2.3. State–Space Method for Large-Scale Microresonators with Phase-Lag Heat Conduction Law

Combining Equations (3) and (6), we obtain

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\mathbf{M}\ddot{\mathbf{U}}+\mathbf{KU}-\mathbf{G}\mathsf{\Theta }=\mathbf{0}\hfill \\ {\mathbf{S}}_{}{\mathsf{\Theta }}_{}+{{\mathbf{S}}_{}}_T{\dot{\mathsf{\Theta }}}_{}+{\mathbf{H}}_{}{\dot{\mathsf{\Theta }}}_{}+{\mathbf{F}}_{}{\dot{\mathbf{U}}}_{}+{{\mathbf{H}}_{}}_q{\ddot{\mathsf{\Theta }}}_{}+{{\mathbf{F}}_{}}_q{\ddot{\mathbf{U}}}_{}=\mathbf{0}\hfill \end{array}\right.\end{array}$$

(10)

Equation (10) represents a unified thermoelastic governing equation that takes into account different heat conduction models. The full coupling between the temperature field and displacement field in the equation of motion and the heat conduction equation requires a simultaneous solution of the governing equations of heat conduction and motion. Moreover, the introduction of high-order terms related to the phase-lag time, such as ${\mathbf{S}}_T\dot{\mathsf{\Theta }}$, ${\mathbf{H}}_q\ddot{\mathsf{\Theta }}$, and ${\mathbf{F}}_q\ddot{\mathbf{U}}$, further increases the difficulty of solving this fully-coupled system equation.

There are several ways to solve the combined Equation (10). Direct decoupling of (10) is usually achieved by ignoring the coupling terms in the motion equation or the heat conduction equation, but this destroys the original equation form and neglects the impact of thermoelastic behavior on energy dissipation. Iterative methods such as the Newton iteration method are generally inefficient for solving fully-coupled high-order equations. Therefore, based on the idea of expanded complexity and descended order, this paper constructs a low-order state–space model that retains the thermoelastic influence while preserving the possibility of a subsequent accelerated solution. Unit matrices related to displacement rate and temperature rate are added:

$$\begin{array}{c}\hfill -\mathbf{I}\dot{\mathbf{U}}+\mathbf{I}\dot{\mathbf{U}}=\mathbf{0}\end{array}$$

(11)

$$\begin{array}{c}\hfill -\mathbf{I}{\dot{\mathsf{\Theta }}}_{}+\mathbf{I}{\dot{\mathsf{\Theta }}}_{}=\mathbf{0}\end{array}$$

(12)

Filling in the zero terms, we obtain the state–space method of large-scale microresonators with the phase-lag heat conduction law as

$$\left[\begin{array}{cccc}\mathbf{H}& {\mathbf{H}}_q& \mathbf{F}& {\mathbf{F}}_q\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{M}\\ \mathbf{I}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{I}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}\dot{\mathsf{\Theta }}\\ \ddot{\mathsf{\Theta }}\\ \dot{\mathbf{U}}\\ \ddot{\mathbf{U}}\end{array}\right]+\left[\begin{array}{cccc}\mathbf{S}& {\mathbf{S}}_T& \mathbf{0}& \mathbf{0}\\ -\mathbf{G}& \mathbf{0}& \mathbf{K}& \mathbf{0}\\ \mathbf{0}& -\mathbf{I}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& -\mathbf{I}\end{array}\right]\left[\begin{array}{c}\mathsf{\Theta }\\ \dot{\mathsf{\Theta }}\\ \mathbf{U}\\ \dot{\mathbf{U}}\end{array}\right]=\mathbf{0}$$

(13)

Equation (13) represents the governing equation of the system taking into account non-Fourier phase-lag heat conduction. In the 2D thermoelastic eigenvalue problem, the number of nodes $N_n$ directly affects the degrees of freedom of the system, where the basic degrees of freedom $N_d$ are given by $N_d=3N_n$ (including two displacement degrees of freedom and one temperature degree of freedom per node). However, the state–space method increases the degrees of freedom, and for the non-Fourier dual-phase lag heat conduction combined Equation (20), the degrees of freedom of the state–space solution increase to $N_s=6N_n$. Similarly, for the 3D thermoelastic eigenvalue problem, the basic degrees of freedom are $N_d=4N_n$, and for the non-Fourier dual-phase lag heat conduction combined Equation (20), the dimension of the state–space solution increases to $N_s=8N_n$.

With the help of the state–space method, for the free vibration considered herein, the solution is assumed to take the form:

$$\begin{array}{c}\hfill \mathsf{\Theta }=\tilde{\mathsf{\Theta }}{e}^{{\lambda }_{i}t},\mathbf{U}=\tilde{\mathbf{U}}{e}^{{\lambda }_{i}t}\end{array}$$

(14)

where ${\lambda }_i$ is the i-th eigenvalue of the thermoelastic system, ${\left[\tilde{\mathsf{\Theta }}\right]}_{N_n\times 1}$ and ${\left[\tilde{\mathbf{U}}\right]}_{2N_n\times 1}$ or ${\left[\tilde{\mathbf{U}}\right]}_{3N_n\times 1}$ are the amplitudes of temperature and 2D/3D displacement, respectively.

Then, Equation (13) can be expressed as

$$\begin{array}{c}\hfill \left({\lambda }_i\mathbf{B}-\mathbf{A}\right){\mathsf{\Phi }}_i=\mathbf{0}\end{array}$$

(15)

where, for the 2D non-Fourier dual-phase lag heat conduction model, ${\left[{\mathsf{\Phi }}_i\right]}_{6N_n\times 1}={\left[\begin{array}{c}\ddot{\mathsf{\Theta }},\dot{\mathsf{\Theta }},\ddot{\mathbf{U}},\dot{\mathbf{U}}\end{array}\right]}^{\mathrm{T}}$, and in the case of 3D, ${\left[{\mathsf{\Phi }}_i\right]}_{8N_n\times 1}={\left[\begin{array}{c}\ddot{\mathsf{\Theta }},\dot{\mathsf{\Theta }},\dot{\mathbf{U}},\ddot{\mathbf{U}}\end{array}\right]}^{\mathrm{T}}$.

Since the degrees of freedom of the system matrix reach $6N_n$ and $8N_n$ for the 2D and 3D thermoelastic problems, respectively, solving the system typically incurs heavy computational costs. Moreover, for large-scale design problems, the intermediate variables of the expanded system may exceed the storage capacity of computing devices. Once the thermoelastic system matrices are constructed, an optimal design method must be employed to solve the quality factor design problem.

2.4. Definition of Optimization Problem

The objective of the topology optimization process is to minimize the thermoelastic dissipation of the microbeam resonator. There are several methods for material structure topology optimization design, which can be roughly classified into gradient-based and non-gradient methods. The gradient-based method determines the search direction by computing the sensitivity gradient of the design variables. One specific branch of the gradient-based method is evolutionary structural optimization and its derivative algorithms, where the element pseudo-densities are the design variables. In this approach, a value of “1” for an element pseudo-density signifies that the element is filled with material, while a value of “0” denotes that the element is empty. By iteratively adjusting the element pseudo-densities based on the sensitivity of the design variable and the optimization evolution rate, the material in the structure is either added or removed, thereby altering its dynamical performance. The evolutionary iterative optimization process converges when the optimized structure reaches the target volume fraction within the constraints, and the objective function converges.

The advantage of evolutionary structural optimization lies in its clear and concise concept of element sensitivity. Elements with high design variable sensitivity make a greater contribution to the objective function and are retained, while those with low sensitivity are removed [37]. In addition, the structure obtained by its optimization has no intermediate density and is closer to the “0-1” topology configuration. Furthermore, since the iterative algorithm considers the relative values of the design variable sensitivity, there is no strict requirement for the positive or negative signs of the sensitivity, making it more convenient in handling multiphysics optimization problems than other optimization methods such as density-based methods. Thus, the evolutionary structural optimization method is chosen as the design method to determine the search direction of the design problem.

For the optimization of microstructures with the non-Fourier phase-lag heat conduction law, the thermoelastic optimization model can be expressed in a specific structural topology optimization form by combining topology optimization-related constraints, design variables $p_e$, and objective function $\mathbb{C}$. Taking the i-th quality factor $Q_i$ of the microresonator as the objective function, the matrix form of the thermoelastic optimization problem can be obtained as follows:

$$\begin{array}{cc}\hfill \mathrm{Max}& \mathbb{C}=Q_i\hfill \end{array}$$

$$\begin{array}{cc}\hfill s. t.& \mathbf{A}{\mathsf{\Phi }}_i={\lambda }_i\mathbf{B}{\mathsf{\Phi }}_i\hfill \end{array}$$

(16a)

$$\begin{array}{cc}& {\displaystyle \sum _{e=1}^N}{\mathsf{\Omega }}_e{p}_e-l_v{\mathsf{\Omega }}_0\le 0\hfill \end{array}$$

(16b)

$$\begin{array}{cc}& l_f{f}^{\left(\mathrm{origin}\right)}-f\le 0\hfill \end{array}$$

(16c)

$$\begin{array}{cc}& p_e=p_{min}\mathrm{or}1,(e=1,2,\cdots N)\hfill \end{array}$$

(16d)

The eigenvalue equation in (16a) is the main constraint of the optimization problem, where $\mathbf{A}$ and $\mathbf{B}$ are the system matrices mentioned earlier. Through finite element analysis, both free vibration boundary conditions and adiabatic boundary conditions are included in (16a). In addition, volume constraint (16b) and element volume constraint (16d) are imposed on the optimization process to obtain a lightweight structure and ensure that the optimization process is not affected by the singularity of the overall system matrix. Here, $p_e$ is the pseudo-density of the e-th element, which takes the value of 1 or a small positive number close to zero $p_{\min }$. N is the number of finite elements in the entire structure, $l_v$ is the volume fraction of the optimized structure, ${\mathsf{\Omega }}_e$ is the volume of each element, and ${\mathsf{\Omega }}_0$ is the initial volume of the structure. Furthermore, a frequency constraint (16c) is imposed to ensure the stability of the optimization process and prevent excessive reduction of structural frequency.

As mentioned before, this paper mainly focuses on the evolutionary structural optimization method to improve the quality factor of the structure. The sensitivity of the i-th order eigenvalue ${\lambda }_i$ of the system to the design variable $p_e$ can be expressed as:

$$\frac{\partial {\lambda }_i}{\partial p_e}={\mathbf{\Psi }}_i^{\mathrm{T}}\left(\frac{\partial \mathbf{A}}{\partial p_e}-{\lambda }_i\frac{\partial \mathbf{B}}{\partial p_e}\right){\mathsf{\Phi }}_i$$

(17)

As shown in Equation (17), the right-hand side can be obtained by solving the system equation. Therefore, this equation reveals an explicit relationship between system eigenvalue sensitivity and design variables. As mentioned before, we choose an evolutionary structural topology optimization method based on finite element analysis to solve the fully-coupled thermoelastic optimization problem in Equation (16). Since gradient-based evolutionary structural optimization methods require the sensitivity of the objective function with respect to the design variables, and considering that the displacement field, temperature field, and characteristic frequency of the system may be functions of the design variables themselves, according to the chain rule, the sensitivity of the objective function with respect to the design variables can be expressed as follows. Consequently, the sensitivity of the quality factor Q with respect to the element pseudo-density $p_e$ can be obtained using the chain rule and the relationship between the system’s eigenvalues and quality factor:

$$\begin{array}{c}\hfill \frac{\partial Q^{-1}}{\partial p_e}=2{\displaystyle \frac{\Im \left({\displaystyle \frac{\partial \lambda i}{\partial pe}}\right)\Re \left({\lambda }_i\right)-\Im \left({\lambda }_i\right)\Re \left({\displaystyle \frac{\partial \lambda i}{\partial pe}}\right)}{{\left[\Re \left({\lambda }_i\right)\right]}^2}}\end{array}$$

(18)

This equation represents the explicit expression of the sensitivity of the objective function Q with respect to the design variable $p_e$, and it contains terms related to the hysteresis time (such as ${\mathbf{S}}_T$, ${\mathbf{H}}_q$, and ${\mathbf{F}}_q$) as indicated by the definition of ${\lambda }_i$.

The optimization process of an ESO method can be discretized into several steps, as shown in Figure 2. For a large-scale design problem, the finite element analysis step of the fully-coupled system, as marked in red, requires heavy computational cost, while the repeated reconstruction of the structure further exacerbates the computational burden.

According to the definition and standard process of the ESO method, the convergence of the optimization process can be acquired when the structure reaches the design volume fraction, and the objective function changes little between two iterations [13].

As mentioned before, the search direction, search domain, efficiency, and convergence are vital to the design process. For instance, the gradient-based ESO method determines the search direction. The search domain needs to be narrowed down. The efficiency of the method needs to be improved, and the convergence of the optimization process needs to be validated. In the next section, we will propose an efficient and robust topology optimization method to solve this optimization problem with efficiency and accuracy.

3. Efficient and Robust Topology Optimization Method

In this section, a modal analysis method based on a Neumann series expansion is proposed to efficiently solve the combined system equation (Equation (13)). Then, a modal validation process based on the modal assurance criteria is implemented to narrow down the search domain, along with the posterior evaluation technique to retain the convergence of the design method.

3.1. Reconstruction of the System Matrices

According to the definition of thermoelastic damping, the system matrices become complex when considering the effects of thermoelastic damping. For the governing equation of motion, let $\mathbf{M}$ be a positive definite mass matrix, $\mathbf{K}$ be a non-negative definite symmetric stiffness matrix, and $\mathbf{G}$ be the thermal expansion matrix in the equation of motion. As mentioned above, the asymmetry of the combined system matrices $\mathbf{A}$ and $\mathbf{B}$ leads to the complex eigenvalue problem, which is difficult and time-consuming to solve. The key difficulty in solving the eigensolution is the decomposition and inversion step of the combined matrix. To facilitate the decomposition and inversion of the combined terms in the system equations, we construct the lumped mass matrix by making its diagonal elements proportional to the corresponding diagonal elements of the original consistent mass matrix. The lumped mass matrix is a diagonally positive definite matrix obtained from the traditional consistent mass matrix, which is widely used in finite element dynamic analysis. Its accuracy has been verified through numerical examples [38,39].

It is assumed that the total sum of all elements in both types of mass matrices is equal, that is

$$\begin{array}{c}\hfill {\left[\mathbf{M}\right]}_{ii}^{\mathrm{lumped}}=\frac{{\left[\mathbf{M}\right]}_{ii}^{\mathrm{consistent}}}{\sum {\left[\mathbf{M}\right]}_{ii}^{\mathrm{consistent}}}\sum \sum {\left[\mathbf{M}\right]}_{ij}^{\mathrm{consistent}}\end{array}$$

(19)

where ${\left[\mathbf{M}\right]}_{ij}$ represents the element in the i-th row and j-th column of matrix $\mathbf{M}$, and ∑ is the summation operator for the elements in the matrix. For convenience, we use $\mathbf{M}$ to denote the lumped mass matrix in the following text. Similarly, the heat generation matrix $\mathbf{H}$ in the heat conduction control equation is also constructed as a positive definite diagonal matrix.

For thermal-elastic coupling systems that consider non-Fourier dual-phase lag heat conduction models, the full coupling between the motion equation and the heat conduction equation requires a simultaneous solution of both equations. In the previous chapters, the order of the thermoelastic system was reduced by applying a state–space transformation to the system equations due to the influence of high-order time derivatives in the motion and heat conduction system equations. However, this increases the dimensionality of the system equations, and the classic direct combination method for control equations shown in Equation (13) results in two asymmetric matrices for the inherent terms of the system equations. This leads to a slow direct solution of the system equations due to the inversion and decomposition process of the global matrix.

To avoid the extremely time-consuming decomposition and inversion steps in solving the asymmetric matrix of the system equation, this paper combines the system matrices $\mathbf{M}$ and $\mathbf{H}$ into a single matrix to reduce computation time. Moreover, for the thermoelastic system with the non-Fourier dual-phase lag heat conduction model considered in this paper, the second-order term related to time cannot be directly transformed. Therefore, the first term on the left side of the equation is constructed as a triangular matrix to simplify the inversion and decomposition steps:

$$\left[\begin{array}{cccc}{\mathbf{H}}_q& \mathbf{0}& {\mathbf{F}}_q& \mathbf{0}\\ \mathbf{0}& \mathbf{I}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{M}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}\end{array}\right]\left[\begin{array}{c}\ddot{\mathsf{\Theta }}\\ \dot{\mathsf{\Theta }}\\ \ddot{\mathbf{U}}\\ \dot{\mathbf{U}}\end{array}\right]+\left[\begin{array}{cccc}{\mathbf{S}}_T+\mathbf{H}& \mathbf{S}+\mathbf{F}& \mathbf{0}& \mathbf{0}\\ -\mathbf{I}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& -\mathbf{G}& \mathbf{0}& \mathbf{K}\\ \mathbf{0}& \mathbf{0}& -\mathbf{I}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}\dot{\mathsf{\Theta }}\\ \mathsf{\Theta }\\ \dot{\mathbf{U}}\\ \mathbf{U}\end{array}\right]=\mathbf{0}$$

(20)

Thus, Equation (13) can be restructured into a form similar to the standard characteristic problem while preserving the original structure of the control equation to the greatest extent possible.

Equation (20) represents the eigenvalue problem of thermoelastic systems under free vibration considering the non-Fourier phase-lag heat conduction model. Since the first term is a triangular matrix, the inversion step is computationally efficient. However, the second term in the equation is asymmetric and may not be positive definite, which increases the computational cost of solving the eigenvalue problem of the combined equation. In the following text, we adopt a subspace method to approximate the solution in a small degree of freedom system to reduce the computational cost.

For the reconstructed system equations considering different heat conduction models, the combination matrix $\mathbf{B}$ is constructed as a triangular matrix, which greatly simplifies the inversion and decomposition of this matrix compared to the original system equations.

3.2. Reduced Basis in Neumann Series Subspace

During the iteration of topology optimization, structural configuration changes often occur by changing its material distribution to improve the mechanical system’s performance. Therefore, the system matrix change caused by changes in structural design parameters can be represented as follows:

$$\begin{array}{c}\hfill {\mathbf{A}}_m=\mathbf{A}+\Delta \mathbf{A},{\mathbf{B}}_m=\mathbf{B}+\Delta \mathbf{B}\end{array}$$

(21)

where ${\mathbf{B}}_m$ and ${\mathbf{A}}_m$ are the system matrices after structural configuration changes, and $\Delta \mathbf{B}$ and $\Delta \mathbf{A}$ represent the changes of the system matrices. The eigenvalue problem after structural configuration changes can be expressed as

$$\begin{array}{c}\hfill \left({\lambda }_{mi}{\mathbf{B}}_m-{\mathbf{A}}_m\right){\mathsf{\Phi }}_{mi}=\mathbf{0}\end{array}$$

(22)

where ${\lambda }_{mi}$ and ${\mathsf{\Phi }}_{mi}$ are, respectively, the ith eigenvalue and eigenvector of the modified viscoelastic structures.

Recalling Equation (21) and using some matrix manipulations, the modified eigenproblem (22) can be rewritten as

$$\begin{array}{c}\hfill \left({\mathbf{I}}_N-{\xi }_i{\mathbf{B}}_m^{-1}\Delta \mathbf{A}\right){\mathsf{\Phi }}_{mi}={\xi }_i{\mathbf{B}}_m^{-1}\mathbf{A}{\mathsf{\Phi }}_{mi}\end{array}$$

(23)

where

$$\begin{array}{c}\hfill {\xi }_i=\frac{1}{{\lambda }_{mi}}\end{array}$$

(24)

Applying Neumann series expansion to implement the inverse matrix of the coefficient matrix, ${\left({\mathbf{I}}_N-{\xi }_i{\mathbf{B}}_m^{-1}\Delta \mathbf{A}\right)}^{-1}$, one obtains

$$\begin{array}{c}\hfill {\mathsf{\Phi }}_{mi}={\xi }_i\sum _{k=0}^{\infty }{\left({\xi }_i{\mathbf{B}}_m^{-1}\Delta \mathbf{A}\right)}^k{\mathbf{B}}_m^{-1}\mathbf{A}{\mathsf{\Phi }}_{mi}\end{array}$$

(25)

Using the first r expansion terms and the initial eigenvector ${\mathsf{\Phi }}_i$ to approximate the modified eigenvector ${\mathsf{\Phi }}_{mi}$, we can obtain an approximate approach for calculating the modified eigenvector ${\mathsf{\Phi }}_{mi}$ as

$$\begin{array}{c}\hfill {\mathsf{\Phi }}_{mi}={\xi }_i\sum _{k=0}^r{\left({\xi }_i{\mathbf{B}}_m^{-1}\Delta \mathbf{A}\right)}^k{\mathbf{B}}_m^{-1}\mathbf{A}{\mathsf{\Phi }}_i\end{array}$$

(26)

Definition 1

(Neumann series subspace). The Neumann series subspace ${\mathcal{R}}_r$ is formed by an r-dimensional recursive sequence $\left\{r_i^{\left(0\right)},r_i^{\left(1\right)},r_i^{\left(2\right)},\cdots ,r_i^{\left(n-1\right)}\right\}$ generated by original system matrices $\mathbf{A}$ and $\mathbf{B}$, the changes in original system matrices $\Delta \mathbf{A}$ and $\Delta \mathbf{B}$, as well as the ith eigenvector ${\mathsf{\Phi }}_i$ of the original eigenproblem defined in Equation (18). That is,

$$\begin{array}{c}\hfill {\mathcal{R}}_r(\mathbf{A},\mathbf{B};{\mathbf{A}}_m,{\mathbf{B}}_m;{\mathsf{\Phi }}_i)\doteq \mathrm{span}\left\{{\mathbf{r}}_i^{\left(0\right)},{\mathbf{r}}_i^{\left(1\right)},\cdots ,{\mathbf{r}}_i^{(r-1)}\right\}\end{array}$$

(27)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{r}}_i^{\left(0\right)}={\left(\mathbf{B}+\Delta \mathbf{B}\right)}^{-1}\mathbf{A}{\mathsf{\Phi }}_i\hfill \\ {\mathbf{r}}_i^{\left(k\right)}={\left(\mathbf{B}+\Delta \mathbf{B}\right)}^{-1}\Delta {\mathbf{Ar}}_i^{(k-1)},\forall k=2,3\cdots ,r\hfill \end{array}\end{array}$$

(28)

Using a Neumann series subspace ${\mathcal{R}}_r$, it is convenient to use known modes to approximate the modified eigenvectors as the following form according to the approximate approach defined in Equation (26).

$$\begin{array}{c}\hfill {\mathsf{\Phi }}_{mi}\approx \sum _{k=0}^r{\xi }_i^{-(k+1)}{\mathbf{r}}_i^{\left(k\right)}\end{array}$$

(29)

where the r-dimensional recursive sequence has already been defined in (28).

Remark 1.

In the classic subspace method for solving the characteristic equation, the most computationally intensive part is often the inversion step, i.e., the process of solving ${{\mathbf{B}}_{\mathbf{m}}}^{-1}$. However, for the reconstructed system matrix ${\mathbf{B}}_m$ obtained from the lumped mass matrix $\mathbf{M}$ and heat generation matrix $\mathbf{H}$, which were previously discussed, the matrix decomposition and inversion steps have very low computational cost since ${\mathbf{B}}_m$ is a diagonal matrix.

3.3. Projection-Based Model Reduction Method

To maintain the robustness of the modified eigenproblem while reducing the degrees of freedom of the system, this section uses the projection-based model reduction method based on the Neumann series subspace ${\mathcal{R}}_r$.

Based on the r-th order Neumann subspace basis vectors, the modified eigenvector can be expressed as a linear combination of the original high-fidelity subspace in the Neumann series subspace ${\mathcal{R}}_r$:

$$\begin{array}{c}\hfill {\mathsf{\Phi }}_{mi}=\sum _{k=0}^{r-1}{\alpha }_i^{\left(k\right)}\mathbf{q}{}_i^{\left(k\right)}={\mathbf{R}}_i{\mathbf{\alpha }}_i\end{array}$$

(30)

where ${\mathbf{R}}_i$ is the space spanned by the orthonormalization basis ${\mathbf{q}}_i^{\left(k\right)}$ in the Neumann series subspace ${\mathcal{R}}_r$ and ${\mathbf{\alpha }}_i$ is the unknown parameters to be solved. Note that the previous equation used the following relation:

$${\mathbf{R}}_i=\left[\begin{array}{cccc}{\mathbf{q}}_i^{\left(0\right)},& {\mathbf{q}}_i^{\left(1\right)},& \dots ,& {\mathbf{q}}_i^{(r-1)}\end{array}\right],{\mathbf{\alpha }}_i={\left[\begin{array}{cccc}{\alpha }_i^{\left(0\right)},& {\alpha }_i^{\left(1\right)},& \dots ,& {\alpha }_i^{(r-1)}\end{array}\right]}^{\mathrm{T}}$$

(31)

According to the projection-based model reduction method, the degrees of freedom of the system matrices can be reduced as the following form [25]:

$$\begin{array}{c}\hfill {\mathbf{B}}_R={\mathbf{R}}_i^{\mathrm{T}}{\mathbf{B}}_m{\mathbf{R}}_i,{\mathbf{A}}_R={\mathbf{R}}_i^{\mathrm{T}}{\mathbf{A}}_m{\mathbf{R}}_i\end{array}$$

(32)

Because $N_s\gg r$, the reduced system matrices are very small.

Remark 2.

Because ${\mathbf{R}}_i$ is a reduced-basis matrix operator and ${\mathbf{A}}_m$ is not necessarily a self-adjoint matrix, the robustness of the rectangular basis matrix ${\mathbf{R}}_i$ must be improved by using the Gram–Schmidt orthogonalization algorithm. Thus, ${\mathbf{R}}_i$ is the space spanned by the orthonormalization basis ${\mathbf{q}}_i^{\left(k\right)}$ and can be generated by using the r-dimensional recursive sequence already defined in (28) based on the Gram–Schmidt orthogonalization algorithm. That is, the rectangular basis matrix ${\mathbf{R}}_i$ should be normalized such that each column becomes linearly independent.

Since the approximate eigenvectors should satisfy the eigenproblem (22), the unknown coefficients ${\mathbf{\alpha }}_i$ can be determined by substituting the approximate eigenvectors into Equation (22). Combining Equations (22) and (32), the standard eigenvalue problem of the system can be replaced by a reduced eigenvalue problem:

$$\begin{array}{c}\hfill \left({\lambda }_{mi}{\mathbf{B}}_R-{\mathbf{A}}_R\right){\mathbf{\alpha }}_i=\mathbf{0}\end{array}$$

(33)

The dimension of the reduced system matrix is $r\times r$, which is much smaller than the original system matrix $N_s\times N_s$. Therefore, it is only necessary to solve the reduced eigenvalue problem, which reduces the size of the eigenvalue problem from $N_s\times N_s$ to $r\times r$. Thus, the projection-based model reduction method can save the computational cost significantly. After solving for coefficients ${\mathbf{\alpha }}_i$, the modified eigenvector can be obtained using Equation (30).

As shown in Algorithm 1, the construction of the r order subspace of the Neumann series expansion combined with the reduced-basis solution improves the efficiency of the modal analysis. Numerous studies have shown that the computational complexity of the traditional method for solving the combined eigenvalue equation is $O\left(N_s^3\right)$ [40,41]. Since the calculation scale of solving the basis vectors is relatively large compared to other steps, we mainly focus on the computational complexity of solving the basis vectors. For the system matrix $\mathbf{A}$, the computational complexity of solving the left and right multiplication of the basis vectors is $4Lr{N}_{s}b$, where b is the half-bandwidth of the combination matrix, r is the number of basis vectors, and L is the order of eigenvalues to be solved; the half-bandwidth b is approximately proportional to $N_s^{0.5}$ [41]. Therefore, the computational complexity of the complex modal analysis process for thermoelastic structures is $O\left(4Lr{N}_s^{1.5}\right)$.

Thus, for large-scale problems with a large $N_s$, since $r\ll N_s$ and $L\ll N_s$, the projection-based model reduction method based on the Neumann series subspace ${\mathcal{R}}_r$ has lower computational cost compared to conventional solving methods. However, the projection-based model reduction method based on the Neumann series subspace ${\mathcal{R}}_r$ cannot be directly used in the optimization process due to the instability and convergent problems.

3.4. Modal Validation Process

According to the definition of the quality factor in the case of small dissipation (Equation (1)), the quality factor of the structure is dependent on the eigenvalues and mode shapes of the structure. Thus, during the iterative optimization process, a different mode shape may result in a different quality factor; the same order of mode shape has to be ensured in the optimization process.

For this reason, we use a modal assurance criteria-based (MAC-based) validation process to identify the first-order mode in the iterative optimization process. First, the first-order left and right eigenvectors ${\mathbf{\Psi }}_1^1$ and ${\mathsf{\Phi }}_1^1$ of a solid original structure without cavities are used as the reference group for the first iteration in the optimization process.

Then, the modal assurance criteria are introduced to the process as follows [42]:

$$\begin{array}{c}\hfill {\mathrm{MAC}}_1^{j,j+1}=\frac{{\left[{\left({\mathbf{\Psi }}_1^{j+1}\right)}^{\mathrm{T}}{\mathsf{\Phi }}_1^j\right]}^2}{\left[{\left({\mathbf{\Psi }}_1^j\right)}^{\mathrm{T}}{\mathsf{\Phi }}_1^j\right]\left[{\left({\mathbf{\Psi }}_1^{j+1}\right)}^{\mathrm{T}}{\mathsf{\Phi }}_1^{j+1}\right]}\end{array}$$

(34)

In the above equation, ${\mathbf{\Psi }}_1^j$ and ${\mathsf{\Phi }}_1^j$ represent the first-order left and right eigenvectors in the j-th iterative loop, respectively, while ${\mathbf{\Psi }}_1^{j+1}$ and ${\mathsf{\Phi }}_1^{j+1}$ represent the first-order left and right eigenvectors in the $(j+1)$-th iterative loop. ${\mathrm{MAC}}_1^{j,j+1}$ represents the similarity between the first-order mode in the $(j+1)$-th iterative loop and the original structure’s first-order mode, which ranges from 0 to 1. A larger MAC value indicates a higher similarity between the two modes, while a smaller value indicates a lower similarity. By searching through all modes in the $(j+1)$-th iteration, the mode most similar to the reference group’s left and right eigenvectors ${\mathbf{\Psi }}^j$ and ${\mathsf{\Phi }}^j$ is selected, and the corresponding eigenfrequency is used to calculate the quality factor of the current iteration step to maintain the stability of the structure modes during the optimization process.

Furthermore, as mentioned in the previous section, the solution of the fully-coupled system means heavy computational cost. While in the case of topology optimization problems, the solution of multiple eigenvalues and local modes will significantly affect the search domain and accuracy of the optimization method. To narrow down the search domain and improve the efficiency of the optimization process, the resonance frequency of the working mode shape in the j-th iterative loop will be used as the reference frequency in the next loop and only the system mode shapes with nearby resonance frequencies will be solved in the $(j+1)$-th iterative loop. Then, the working mode shape of the structure can be obtained precisely and efficiently during the optimization process.

From the literature above, the MAC-based validation process can reduce the search domain and stabilize the process of the iterative optimization method.

The computational efficiency of topology optimization methods based on finite element analysis and sensitivity gradients is heavily influenced by the number of degrees of freedom of the system. However, the fully-coupled nature of thermoelastic damping makes the system equations and sensitivity analysis more difficult to solve. Particularly for large and complex thermoelastic structures, the repetitive analysis process in iterative optimization further exacerbates the issue of low computational efficiency. On the other hand, current commercial software mainly focuses on decoupling and analyzing multi-field coupling systems, with optimization design modules only including statics and weak coupling analysis, lacking corresponding design tools for thermoelastic structures. Therefore, this paper proposes combining efficient modal analysis and modal validation techniques, along with the posterior evaluation technique to ensure the convergence of the optimization process, as a solution to the problem of difficulty and low efficiency in solving thermoelastic structures.

3.5. Posterior Evaluation Technique

The convergence of the Neumann series expansion depends on the eigenvalues of its non-identity terms [43]. Therefore, for the convergent Neumann series (26), the necessary and sufficient condition can be obtained as follows:

$$\begin{array}{c}\hfill \rho ({\xi }_i{\mathbf{B}}_m^{-1}\Delta \mathbf{A})<1\end{array}$$

(35)

where $\rho (•)$ denotes the spectral radius of the matrix $(•)$, i.e., it requires that the absolute value of all eigenvalues of the coefficient matrix are less than 1.

The maximum eigenvalue of the matrix ${\mathbf{B}}_m^{-1}\Delta \mathbf{A}$ can be obtained by solving the following eigenvalue problem for its minimum eigenvalue:

$$\begin{array}{c}\hfill {\mathbf{B}}_m{\mathbf{y}}_i=s_i\Delta \mathbf{A}{\mathbf{y}}_i\end{array}$$

(36)

As the thermal-elastic damping of the system is small, the known undamped (real part) frequencies can be used to approximate. Combining the minimum eigenvalue $s_{\min }$ and Equation (24) yields

$$\begin{array}{c}\hfill \left|s_{min}\right|<\left|{\lambda }_i\right|\end{array}$$

(37)

Note that the original eigenvalue ${\lambda }_i$ is used herein to approximate the modified eigenvalue ${\lambda }_{mi}$.

As shown in Equations (24) and (35), the convergence of the Neumann series expansion depends on the system eigenvalues ${\lambda }_i$, system matrix $\mathbf{B}$, and the system matrix perturbation $\Delta \mathbf{A}$. However, due to the thermoelastic damping in the system, the perturbation $\Delta \mathbf{A}$ is uncertain as it may change during the evolutionary optimization process when materials are added or removed from different locations. Therefore, the convergence of the Neumann series expansion cannot be verified by a priori methods. The convergence of the Neumann series expansion directly affects the efficiency and accuracy of the subspace dimensionality reduction method, and thus affects the optimization process. If the series diverges, it may cause oscillations in the iterative optimization process and cannot guarantee the convergence of the optimization algorithm.

As the technique of accelerating the solution of eigenproblems serves the optimization algorithm, and the convergence of the Neumann series cannot be verified in an a priori way, we adopt a posterior evaluation technique to ensure the convergence of the optimization algorithm. As shown in Algorithm 2, the parameter $evalnum$ is the number of iterations the posterior process has taken effect on. In the first iteration of evaluation, the evolution rate is adjusted, while in later iterations, the sensitivities are adjusted to ensure the convergence of the Neumann series expansion. The adjustment of sensitivities is repeated until the Neumann series expansion converges to ensure the robustness of the topology optimization process.

The evolution rate in gradient-based structural optimization, which is the overall volume reduction rate per iteration, directly affects the change of system matrix $\Delta \mathbf{A}$. Therefore, when the Neumann series expansion does not converge, the optimization process first returns to the structure obtained in the previous iteration and adjusts the evolution rate based on it to ensure convergence:

$$\begin{array}{c}\hfill e{r}^{(k+1)}={\mathcal{P}}_{er}\times e{r}^{\left(k\right)}\end{array}$$

(38)

where $e{r}^{\left(k\right)}$ denotes the evolution rate in the k-th iteration, and ${\mathcal{P}}_{er}$ is the penalty function for the evolution rate. After adjusting the evolution rate, a new structure and the corresponding system matrix $\Delta \mathbf{A}$ are obtained through optimization. If the new structure satisfies the convergence of the Neumann series expansion, the optimization process continues.

If the new structure obtained after adjusting the evolution rate still does not satisfy the convergence of the Neumann series expansion, the optimization process returns to the structure from the previous iteration. Then, the weights of the sensitivities of each element are adjusted based on the structure retained from the previous optimization. For the inefficient elements that were discarded in the previous optimization ($e\left(\mathrm{low}\right)$) and the efficient elements that were added in the current optimization ($e\left(\mathrm{high}\right)$), their penalty factors are increased to encourage the algorithm to choose other regions for optimization in the next iteration.

$$\begin{array}{c}\hfill {\left(\frac{\partial {\lambda }_{mi}}{\partial p_{e\left(\mathrm{low}\right)}}\right)}^{(k+1)}={\mathcal{P}}_{e\left(\mathrm{low}\right)}{\left(\frac{\partial {\lambda }_{mi}}{\partial p_{e\left(\mathrm{low}\right)}}\right)}^{\left(k\right)},{\left(\frac{\partial {\lambda }_{mi}}{\partial p_{e\left(\mathrm{high}\right)}}\right)}^{(k+1)}={\mathcal{P}}_{e\left(\mathrm{high}\right)}{\left(\frac{\partial {\lambda }_{mi}}{\partial p_{e\left(\mathrm{high}\right)}}\right)}^{\left(k\right)}\end{array}$$

(39)

where $\mathcal{P}e\left(\mathrm{low}\right)$ and $\mathcal{P}e\left(\mathrm{high}\right)$ are penalty functions for inefficient and efficient elements, respectively. After adjusting the element sensitivities, the optimization obtains a new structure and system matrix change $\Delta \mathbf{A}$, and optimizes for a new structure. If the new structure obtained after adjusting the sensitivities satisfies the convergence condition of the Neumann series expansion, the optimization process continues. Otherwise, based on the optimized elements in the current structure, the process returns to the previous structure and adjusts the element sensitivities to obtain a new structure. This process is repeated until the Neumann series expansion of the new structure converges.

By combining the traditional ESO convergent checks with the efficient modal analysis and posterior evaluation technique, the efficiency and robustness of the topology optimization method can be assured. The flowchart of the whole process is shown in Figure 3.

Step 1: Initialize mesh, boundary condition, evolutionary parameters, and penalty factors.

Step 2: Construct global system matrices using the reconstruction technique in Section 3.1.

Step 3: Efficient modal analysis can be implemented based on Algorithm 1.

Step 4: Use posterior evaluation for Neumann series subspace in Algorithm 2 to ensure convergence so that the proposed method is robust.

Step 5: Validate the mode shape of the system using Equation (34). With the help of the modal assurance criteria, the same order of mode shape is ensured, and the search domain of the optimization method is reduced.

Step 6: Analyze the global sensitivities through Equation (18) and construct a new design according to global sensitivity.

Step 7: Repeat Steps 2–6 until the objective function converges.

4. Results and Discussion

In this section, numerical examples are conducted to compare the iterative process with and without acceleration techniques in order to validate the accuracy and computational cost of the accelerated method. The accelerated method adopts the modified modal assurance criteria and the reduced basis method to solve the complex eigenfrequencies, the first example is conducted to verify its accuracy and computational cost. Then, two large-scale examples are discussed to show the accuracy and efficiency of the method. The accuracy of the finite element analysis of the proposed accelerated method is verified by comparing the results with the commercial software.

Both the non-accelerated and accelerated methods use the same material parameters, that is, the Young’s modulus $157\mathrm{GPa}$, density $2330{\mathrm{kg}/\mathrm{m}}^3$, Poisson’s ratio 0.22, thermal expansion coefficient $2.6\times {10}^{-6}{\mathrm{K}}^{-1}$, specific heat per unit mass $700\mathrm{J}/(\mathrm{kg}·\mathrm{K})$, thermal conductivity $90\mathrm{W}/(\mathrm{m}·\mathrm{K})$.

Moreover, the first-order thermoelastic quality factor ($Q_1$) of the objective function is normalized for better comparison.

4.1. Validation of Efficiency and Accuracy

The dual-phase lag heat conduction model is more representative of the actual operating conditions of microresonators at the microscale and below. Therefore, this section discusses the optimization results of different solutions considering the non-Fourier dual-phase lag heat conduction model.

The clamped–clamped microbeam resonator is often used as the key component of micro devices, and its quality factor will directly affect the performance of these devices [44,45]. Thus, the clamped–clamped microbeam resonator is chosen as the example in the validation process. Considering a clamped–clamped microbeam with the size of 8 × 1 µm, the entire design domain is discretized using a 2D finite element method into a grid mesh of $320\times 40$ 2D isoparametric thermoelastic finite elements, with a total system degree of freedom of 39,483. The evolutionary optimization parameters were set as $l_v=0.8$, $er=0.008$, $r_{\min }=4$, and the material penalty factors were set as $p_E=2$, $p_{\rho }=1$, $p_{\kappa }=3$, $p_{\alpha }=1$, $p_{C_v}=1$. The temperature phase-lag was set to ${\tau }_T=0.1$ µs, while the phase-lag of heat flux was set to ${\tau }_q=0.005$ µs.

As shown in Figure 4, for the case of non-accelerated method, the normalized quality factor of the microbeam after optimization was 130.07% of the initial structure. The optimization result for the accelerated method with the order of reduced basis $r=10$ was $134.39\%$ of the initial structure. The material distribution and contour map of temperature are similar for the optimization results of non-accelerated (Figure 4b) and accelerated (Figure 4c) methods, with some differences at the border of the vacations, this might be caused by the truncation error of mode shape and reduced basis process of the accelerated method.

On the other hand, the error of the normalized quality factor between these two methods is also acceptable, this indicates that the accelerated method can solve the optimization problem with accuracy.

As shown in

Table 1. The comparison of eigensolution and overall computational times of non-accelerated and accelerated design process of the clamped–clamped microbeam resonator.

	Clamped–Clamped Microbeam ($320\times 40$ Mesh)
	Average Eigenvalue Solution (s)	Iteration Loop	Overall Time (s)
non-accelerated method	26.3567	39	1453.64
accelerated method	8.1932	45	859.32
Time saved (%)	$68.91\%$		$40.88\%$

, the average computation time per iteration of the traditional subspace method for solving eigenvalue problems without acceleration is 26.3567 s, while the same step takes 8.1933 s with the accelerated method. As shown in Figure 4d, the non-accelerated method converges after 39 iterations, with a total computation time of 1453.64 s, including other operations such as sensitivity analysis. As shown in Figure 4e, the accelerated method converges after 45 iterations, with a total computation time of 859.32 s. As the posterior evaluation process takes effect in the iterative optimization process, which can be observed in iterations 28 to 32, the accelerated method takes more total iterations compared to the non-accelerated method. Therefore, the acceleration method proposed in this chapter saves 68.91% of the computational time of the eigenvalue solution step. Although in order to maintain the robustness of the optimization process, the accelerated method takes more iterations to achieve an optimized structure than the non-accelerated method, the accelerated method saves approximately 40.88% of the total computation time of the optimization process compared to the non-accelerated method.

Note that for traditional state–space solution of finite element analysis problems of the thermoelastic system, the complexity of the problem will increase exponentially with the degrees of freedom. Thus, for larger scale problems, the percentage of the computational time saved will increase with the scale of the problem, while for some problems, the time and storage consumed may be unbearable and exceed the maximum of many computation devices.

4.2. Optimal Design of Grid Beams

In microelectromechanical systems, periodic grid beam structures are often used as key components of actuators, sensors, and energy-absorbing structures [46,47]. Among them, grid beams considering in-plane vibrations are widely used, such as the capacitive load sensor and actuator structure. For the in-plane vibrating grid beam sensor, its quality factor affects the measurement accuracy and system stability, and its resonance frequency affects the transmission sensitivity. Therefore, the quality factor is vital to the grid beam system, and improving the quality factor will affect the overall performance of the system.

Thus, in this section, an in-plane vibrating grid beam with phase-lag heat conduction law is considered as an example of the design problem of thermoelastic damping. The sketch map of the structure is shown in Figure 5a, it can be observed that the grid structure is symmetric. Therefore, for the convenience of solving, the left side of the symmetrical structure is selected, and an elemental structure of five beams is chosen as an example in this research, as shown in Figure 5b. In the whole structure, the gray part is the non-design region, these parts of the beams have to be unchanged to maintain the connection with the support end, the finite element grid with a size of 5 µm × 132 µm and a discretization of 20 × 528 2D isoparametric thermoelastic finite elements. The yellow part is the design, where each individual beam has a size of 95 µm × 10 µm and a mesh discretized to 380 × 40 2D isoparametric thermoelastic finite elements. For the 2D microgrid beams, the system has a total 267,027 degrees of freedom. For the classic 2D space–state system solution, the total degrees of freedom is 534,054, due to the iterative characteristic of the optimization process and the medium parameters of the solution, the computation cost is inevitably high. In terms of boundary conditions, the displacement boundary condition is considered as the left fixed support constraint, and the thermal boundary condition is considered an insulated thermal boundary.

The material parameters of single crystal silicon are used as reference. Furthermore, as the existing commercial software lack the solution function considering non-Fourier phase lag heat conduction, the Fourier heat conduction is considered in the validation step of the structure (${\tau }_q=0,{\tau }_T=0$). As the energy dissipation rate is crucial for the system, the optimization objective of the thermoelastic is to maximize the system’s quality factor, subject to a volume constraint of 0.8.

First, the proposed 2D efficient thermoelastic finite element analysis solver is validated by comparing the structure’s frequencies and modes with the COMSOL software’s finite element thermoelastic analysis module. As shown in Figure 6, the original full structure’s frequency obtained using this software was $1.4207\times {10}^6$ Hz, while the frequency obtained through the commercial software was $1.4144\times {10}^6$ Hz, with an error of less than 1%. The original structure’s quality factor obtained using the accelerated method was 7456, and the quality factor obtained using the COMSOL software was 7481, with an error also less than 1%. Therefore, it can be concluded that the finite element system equation solver in this paper is stable and feasible.

Due to the lack of thermoelastic damping optimization module in existing commercial software packages, the optimization result is compared with the original structure solved by the proposed finite element method. The non-Fourier dual-phase lag heat conduction law is considered in the optimization step of the structure (${\tau }_q=0.01$ µs, ${\tau }_T=0.1$ µs), and a normalization process is applied to the objective function (quality factor) of the optimized structure to clearly show the improvement [36]. The temperature cloud-chart of the optimized structure is shown in Figure 6. The optimized structure shows an improved quality factor by 343.28%, compared to the original structure.

4.3. Optimal Design of MEMS Uniaxial Mirror

There are many thermoelastic MEMS applications related to microresonators, such as MEMS switches, energy harvesters, and actuators [48]. Among them, the MEMS switch is widely used in medical devices, mobile phones, radio frequency, and other devices [49]. Take the micro-electromechanical single-axis mirror structure in Figure 7 as an example [1], this mirror can be used as a switch-like component for high-speed switching, which requires high accuracy and stability, that is, the thermoelastic damping is small.

Its xoy-plane can be expressed as shown in Figure 7a, where the left side is two vibrating driving beams, the vibrating direction extends along the z-axis, and the right side is the working mirror. Therefore, it can be simplified as a three-dimensional structural optimization problem, as shown in Figure 7b, where the gray part is the non-design domain of the working mirror, the yellow part is the design domain, and the left side is the fixed constraint. The quality factor is considered the objective function and frequency constraints are imposed to reduce the excessive reduction of the structural frequency during the optimization process. The material parameters are the same as in former examples, and considering the degrees of freedom of the structure, only the optimization design in the range of Fourier heat conduction is considered. As shown in Figure 7b, the structural non-design domain mirror size is 400 µm × 800 µm × 40 µm and is discretized into a 80 × 160 × 8 mesh with 3D isoparametric thermoelastic finite elements. The design domain is two driving beams with fixed support on the left side, and the right side is rigidly connected with the mirror surface, the size of a single beam is 400 µm × 100 µm × 40 µm and the beam is discretized into a $80\times 20\times 8$ mesh with 3D isoparametric thermoelastic finite elements. The total degrees of freedom of the 3D system are 590,436. For the classic 3D state–space solution, the total degrees of freedom reaches the number of 1,180,872, which will significantly impair the efficiency of the solution. Moreover, considering the impact of intermediate variables of the solution, there is a high probability that it exceeds the storage limit of conventional computers.

We validate the solution of the 3D thermoelastic finite element analysis of the proposed method, and compare the structural frequency and mode shape with the finite element thermoelastic analysis module in the COMSOL software. As shown in Figure 8a,b, the original structure frequency obtained by using the accelerated method in this paper is $3.9225\times {10}^4$ Hz, and the frequency obtained by commercial software is $3.9105\times {10}^4$ Hz, the error is less than 1%. The quality factor of the original structure obtained by the accelerated solution is 14,164, and the quality factor of the COMSOL software solution is 14,196, and the error is acceptable. Therefore, it can be considered that the finite element system equation solving module in the software of this paper is feasible.

Similarly, as commercial software lacks a thermoelastic optimization module, the optimized structure is compared with only the original structure obtained by finite element analysis of this paper. The temperature contour map of the optimized result is shown in Figure 8c, with grid lines retained and only half of the structure displayed in the z-axis direction on the right side to highlight the internal configuration of the optimized structure. The optimized structural quality factor is 30,110, which is improved by $212.58\%$ compared to the original structure.

Thus, the proposed method saves significant computation time in comparison with the traditional linear method in the sensitivity analysis of the thermodynamically coupled system, and have the capability of the solution of complex systems with a large scale of degrees of freedom, with little loss of accuracy of the optimization result.

5. Conclusions

With the reduction in structural scales, the influence of thermoelastic damping becomes gradually dominant, determining the upper limit of the performance of high-frequency MEMS microresonators. Thus, enlarging thermoelastic damping based on the topology optimization design plays a critical role in improving the dynamical performance of microresonators. However, the challenges of computational cost and robustness are critical obstacles in the topology optimization design of microresonators. This is due to the iterative process of optimizing large-scale viscoelastic structures.

This study proposes an efficient and robust topology optimization method for minimizing the thermoelastic damping of large-scale microresonators. An evolutionary structural optimization method is adopted to passively determine the search direction of optimizing large-scale thermoelastic structures. To efficiently reduce the computational cost of the iterative process of optimizing process, a model reduction method is developed based on the projection-based model reduction method whose reduced basis is generated within the Neumann series subspace. However, the projection-based model reduction method may be unstable when topology modifications are made during an iteration optimization process. To ensure robustness, a modal validation technique is first implemented in the iterative process to stabilize the iteration and narrow down the search domain, and a posterior evaluation of the Neumann series expansion is then developed to retain the convergence of the projection-based model reduction method. Furthermore, the efficiency and accuracy of the proposed topology optimization method are validated through numerical examples.

Three numerical examples are discussed, results show that the proposed method saves significant computation time in comparison with the traditional linear method in the sensitivity analysis of the thermodynamic structures, with little loss of accuracy of the optimization result. Furthermore, the proposed method has the capability of designing complex large-scale structures.

The feasibility of the efficient and robust method is demonstrated in the gradient-based topology optimization problem for thermoelastic structures. Thus, the research work can be implemented in designing large-scale thermoelastic structures under phase-lag heat conduction law. Furthermore, the topology optimization method proposed in this work combines the matrix reconstruction and parameterized model reduction techniques for thermoelastic systems under the non-Fourier heat conduction law. Thus, it is hopeful for the eigensolution analysis of thermoelastic structures under the phase-lag heat conduction law and expedite the application of thermoelastic structures arising in different subjects, including conceptual design problems, damage analysis, and smart or adaptive structures, and further research in these areas is worth pursuing.