Numerical solution method of thermal-elastic damping for MEMS beam based on interval piecewise infinite state representation

Abstract

This invention discloses a numerical solution method for thermoelastic damping of MEMS beams based on interval-segmented infinite state representation, belonging to the field of damping analysis technology for micro / nano resonators. The method first establishes a thermoelastic model of a micro / nano beam containing fractional-order heat conduction. After dimensionless transformation and Galerkin order reduction, it creatively divides the historical integration domain of the fractional-order derivative into short and long intervals: the short interval uses high-precision L1 interpolation, while the long interval uses infinite state representation technology to transform global convolution into local state recursion, and seamlessly connects the information of the two intervals using a sliding time window. Finally, the modal response is obtained and the loss factor is calculated within the PECE framework. This invention improves upon traditional algorithms... O(N2) The computational complexity is reduced to O(N), which greatly improves the computational efficiency of long-term simulation while maintaining high accuracy, providing a powerful tool for the rapid evaluation of the quality factor of MEMS resonators and the low-damping optimization design.

Description

Technical Field

[0001] This invention belongs to the field of damping analysis technology for micro / nano resonators, specifically involving a numerical solution method for thermoelastic damping of MEMS beams based on interval segmented infinite state representation. Background Technology

[0002] Micro-nano electromechanical systems (MEMS / NEMS) resonators are widely used in inertial sensing, frequency control, radio frequency communication, precision measurement, and biochemical detection. For these devices, thermoelastic damping is one of the fundamental sources of internal dissipation affecting quality factor and frequency stability. As device feature sizes shrink further from the micrometer scale to the nanometer scale and operating frequencies increase from kHz to MHz and even GHz, the impact of limited heat propagation velocity, long memory effect, and complex thermo-coupling behavior on device performance becomes increasingly prominent.

[0003] The diffusion-type heat transfer described by the classical Fourier heat conduction law can no longer accurately predict actual damping. To address this, non-Fourier heat conduction models such as the Cattaneo-Vernotte model, which contain fractional derivatives, have been introduced to capture the historical memory effect and finite propagation speed of heat flow.

[0004] However, when the governing equations include Caputo fractional derivatives, the system exhibits significant global history dependency during time progression. If a global convolutional L1 discretization method is directly used, solving at the nth time step requires backtracking through all historical states from the initial time to the current time, resulting in a total time complexity typically O(N). 2 The storage complexity increases linearly with the number of time steps. This computational cost presents a significant bottleneck for extracting long-term steady-state responses.

[0005] While fast convolution algorithms, convolution quadrature methods, and general fractional prediction-correction methods can improve solution efficiency in some scenarios, they often fail to reduce the time complexity to a linear level of O(N) for the thermoelastic damping problem of micro / nano beams, or they suffer from insufficient stability and significant accuracy loss when facing long-term simulations. Summary of the Invention

[0006] To address the aforementioned technical challenges, this invention proposes a numerical solution method for the thermoelastic damping of MEMS beams based on interval-segmented infinite state representation. This method achieves a fast solution that balances high accuracy, high stability, and linear time complexity, breaking through the computational bottleneck in the current analysis and optimization design of micro / nano resonators.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: A numerical solution method for thermoelastic damping of MEMS beams based on interval piecewise infinite state representation includes the following steps: S1. Obtain the geometric parameters, material parameters, excitation parameters and boundary conditions of the target micro / nano beam resonator. Based on the Euler-Bernoulli beam theory, uniaxial stress assumption and one-dimensional non-Fourier heat conduction assumption in the thickness direction, establish the thermo-elastic coupling control equation containing fractional-order Catteneo-Vernotte heat conduction relationship and dissipation nonlinear term to characterize the thermo-mechanical coupling behavior of the resonator. S2. The thermo-elastic coupling control equations are made dimensionless by defining dimensionless temperature, dimensionless time, dimensionless thickness coordinates, dimensionless driving frequency, and dimensionless thermal relaxation time, and a single-cycle dissipation energy dissipation model is constructed. With maximum energy storage The expression for the loss factor is: ; S3. Using the Galerkin modal projection method, the dimensionless partial differential equation system is reduced to a finite-dimensional modal equation system, which includes Caputo fractional derivative terms, and the modal cutoff order is determined according to the preset accuracy. S4. On a uniform time grid, the historical integration domain involved by the Caputo fractional derivative at the current time is divided into short-interval historical integration and long-interval historical integration; wherein, the short interval covers the recent period containing the current time and the kernel function exhibits weak singularity, and the long interval covers the distant historical period far beyond the recent period and the kernel function exhibits smooth decay. S5. Discretize the short-interval historical integral using the L1 piecewise linear interpolation method to obtain the recent memory term, which represents the influence of recent history on the modal variables at the current moment. S6. The long-interval historical integral is reconstructed using an infinite state representation method, transforming the far-historical contribution term in the form of global convolution into a local recursive form based on a set of internal state variables. S7. The internal state variables are updated using the Crank-Nicolson difference scheme, and the historical information that moves out of the short interval window over time is incorporated into the internal state variables in real time through a preset short interval history queue, so as to achieve seamless connection between recent and distant historical contributions. S8. Under the PECE prediction-evaluation-correction-evaluation time-progression framework, the modal temperature response at the current moment is iteratively solved by integrating recent memory terms and recursively updated internal state variables. S9. Once the system reaches a steady state, calculate the loss factor and output the loss factor for use in the quality factor evaluation and structural design of the micro / nano beam resonator.

[0008] The thermo-elastic coupling control equations established in S1 include fractional-order time derivative terms and thermal relaxation time constants for describing the non-instantaneous response of heat flow to temperature gradient at the micro-nano scale, as well as nonlinear heat source terms for describing the interaction between the temperature field and the strain field.

[0009] The specific steps in S4 are as follows: Set the current time... Historical integral domain Divided into long intervals and short interval ,in The preset short interval length is defined as the time step Δ. t Integer multiples of.

[0010] The specific steps in S6 are as follows: S61. The far-history convolution term is represented by a diffusion integral, and the integral is discretized into P nodes in the pseudo-frequency domain using the generalized Gauss-Laguerre orthogonal quadrature rule. S62. Define corresponding internal state variables for each quadrature node. , ,in The decay parameters for the corresponding quadrature nodes, For known driving terms related to the modal variables; The recursive solution of the first-order ordinary differential equations replaces the repeated convolution calculations for all distant histories.

[0011] The short interval history queue in S7 is a first-in-first-out queue. During time stepping, newly calculated recent historical states are pushed into the head of the queue. At the same time, historical state items that cross the short interval boundary and are removed from the tail of the queue are converted into known excitation sources to drive the update of the internal state variables.

[0012] The PECE time-progression framework in S8 specifically includes: Prediction: Based on the state at the previous time step, explicitly predict the initial state value at the current time step; Evaluation: Using the initial state value, calculate the nonlinear term of the current step, the recent memory term obtained in step S5, and the internal state variable updated in step S7; Correction: Substitute the evaluation results into the discrete equations and solve for the correction value of the state at the current time. Re-evaluation: The correction value is used to re-evaluate, in order to iteratively improve the solution accuracy of coupled nonlinear problems.

[0013] This invention also provides a design method for micro / nano beam resonators. Using the above method, a wideband scanning simulation is performed on a set of resonator design schemes with different geometric or material parameters to obtain the peak value of the loss factor corresponding to each scheme. Based on the correlation between the peak value of the loss factor and the resonator quality factor, the design scheme that meets the target performance is selected.

[0014] Compared with the prior art, the beneficial effects of the present invention are: 1. By using long-interval ISR recursion, the computational bottleneck of global history convolution was successfully overcome, reducing the simulation time complexity from O(N) to O(N). 2 The computation time is reduced to O(N). For long-term simulations requiring tens of thousands or even millions of time steps, the computation time is reduced from hours to minutes or even seconds, making rapid industrial design iteration possible.

[0015] 2. By retaining the short interval L1 high-precision discretization, the weak singularity effect of recent history is accurately captured, avoiding the accuracy loss common in far approximation algorithms.

[0016] 3. The Crank-Nicolson semi-implicit update of internal state variables and the PECE nonlinear processing framework ensure the stability of the algorithm in long-term simulations. This method is not only applicable to silicon beam resonators, but can also be extended to other materials, structures, and even a wider range of fractional-order multiphysics simulation problems.

[0017] 4. This method can directly output the peak value and position of the loss factor, which is directly related to the Q value of the resonator. It provides engineers with quantitative and efficient numerical basis for adjusting the size, material and boundary conditions of the resonator to avoid the high damping frequency band and achieve ultra-high Q value design. Attached Figure Description

[0018] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments are briefly introduced below.

[0019] Figure 1 This is a flowchart illustrating the overall steps of the rapid simulation method of the present invention.

[0020] Figure 2 This is a schematic diagram of the geometric model of the micro / nano beam resonator in an embodiment of the present invention.

[0021] Figure 3 This is a schematic diagram illustrating the principle of the short-interval window sliding and long-interval state update mechanism of the L1-ISR hybrid solver in this invention.

[0022] Figure 4 This is a detailed flowchart of the solver, which includes nonlinear iteration and steady-state determination, in a preferred embodiment of the present invention.

[0023] Figure 5 For traditional global memory convolution algorithms (O(N) 2 Comparison of the theoretical computation time of the linear complexity fast algorithm (O(N)) of this invention with that of the present invention.

[0024] Figure 6 The image shows a comparison of the steady-state temperature evolution waveforms of the global L1 algorithm, the L1-ISR algorithm of this invention, and the reference solution under typical operating conditions.

[0025] Figure 7 A comparison of the convergence of the dimensionless temperature response of microbeams under different modal truncation orders M.

[0026] Figure 8 A comparison of the convergence of the dimensionless temperature response under different numbers of ISR integration nodes P.

[0027] Figure 9 This is a comparison chart showing the actual computation time of the global L1 algorithm and the L1-ISR algorithm of this invention under the same simulation task. Detailed Implementation

[0028] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment is based on data from 12 key sensors collected from a ground test stand for a liquid rocket engine, but the invention is not limited to this specific application scenario.

[0029] Example 1 This embodiment uses the thermoelastic damping simulation of a silicon micro / nano beam resonator as an example. (See also...) Figure 1 The specific steps of the present invention will be described below.

[0030] S1. Construction of the thermo-elastic coupling model like Figure 1 As shown, this embodiment uses a micro / nano-scale beam as the object, with a beam length of L, width of b, and thickness of h. The analysis is performed under pure bending conditions, with simply supported boundary conditions on both sides.

[0031] The geometric parameters, material parameters, excitation parameters and boundary conditions of micro- and nano-scale beams are obtained. Under the Euler-Bernoulli beam theory, the uniaxial stress assumption and the one-dimensional non-Fourier heat conduction assumption in the thickness direction, a thermo-elastic coupling control equation containing fractional-order Catteneo-Vernotte heat conduction relations and dissipation nonlinear terms is established.

[0032] In this step, the heat flux, temperature gradient, and strain gradient are linked through a thermo-elastic coupling term; the fractional derivative is used to characterize the heat flux history memory effect, the thermal relaxation time parameter is used to characterize the finite heat flux propagation velocity, and the dissipative nonlinear term is used to characterize the nonlinear interaction between the temperature field, heat source term, and strain field at the micro-nano scale. This step provides a unified governing equation foundation for subsequent dimensionless transformation, order reduction, and time advancement.

[0033] S2. Dimensionlessness and Construction of Loss Factor Expression The governing equations established in S1 are dimensionless, constructing dimensionless temperature, dimensionless time, dimensionless thickness coordinates, dimensionless driving frequency, and dimensionless thermal relaxation time. Dimensionlessness reduces the impact of differences in the dimensions of different physical quantities on the accuracy and stability of numerical discretization, and facilitates unified comparison of results under different scale parameters.

[0034] Simultaneously, based on the energy definition of thermoelastic damping, a loss factor is established. The calculation expression, where This represents the energy dissipated by the system during a single steady-state cycle. This represents the maximum energy stored within that period. By integrating the energy expression with the output of the subsequent solver, the loss factor at the target frequency can be extracted after time progression is complete.

[0035] S3, Galerkin modal reduction Modal expansion of the continuous temperature field is performed using orthogonal basis functions that satisfy boundary conditions and thermal constraints, and a finite-dimensional modal equation system is obtained through the weighted residual method or orthogonal projection. After this step, the original partial differential equation system is transformed into a finite-dimensional modal dynamic equation system containing Caputo fractional derivatives, which facilitates efficient time propagation in the subsequent process.

[0036] The modal truncation order can be determined comprehensively based on accuracy requirements and computational cost. In a preferred embodiment, a suitable modal order can be selected by comparing the temperature response curves and loss factor convergence under different truncation orders, thus balancing accuracy and solution efficiency.

[0037] S4, Historical Integral Domain Partition For the current moment Caputo fractional derivative history integral field It is divided into long interval historical integrals. Historical integrals over short intervals Among them, the short-interval historical integrals retain the weak singular contributions near the current time, while the long-interval historical integrals correspond to the smoothly decaying distant historical memories.

[0038] To preset the short interval length, the preferred time step is... Multiples of integers. By setting a short interval history queue of fixed length, the accuracy of near history can be guaranteed while the distant history portion is handled by internal state variables. Recursive processing is used to avoid repeatedly backtracking all historical information at each time step.

[0039] S5: L1 Discretization of Short Interval History Integrals For short interval historical integrals The L1 piecewise linear interpolation discretization method is used to approximate the first derivative of the modal variables by the state difference quotient between adjacent time nodes, combined with fractional order. Time step The distances between each historical node and the current node are used to construct locally discrete weights. Since the short interval length is fixed, the summation scale of each step in this part remains within a controllable range, which can more accurately preserve the weak singular memory effect near the current time.

[0040] This step yields a discrete summation of recent historical memory items, which serves to provide a high-precision description of the most sensitive historical contributions at the current moment. Unlike simple truncation or completely equal-interval coarsening, this invention maintains a high-resolution representation of the recent historical portion, thus avoiding a significant sacrifice of numerical accuracy in recent memory in pursuit of reduced complexity.

[0041] S6: ISR Reconstruction of Long-Interval Historical Integrals For long interval historical integrals The infinite-state representation ISR is used for diffusion integral reconstruction. Specifically, the far-history convolution term is first represented by diffusion integral, then the weak singularity near the lower bound of integration is weakened by generalized nonlinear variable substitution, and finally the continuous pseudo-frequency domain is discretized using generalized Gauss-Laguerre quadrature or combined Gaussian quadrature. Each node is a quadrature node, and its internal state variables are constructed. .

[0042] In a preferred embodiment, each internal state variable Satisfying first-order ordinary differential equations ,in The decay parameters for the corresponding quadrature nodes, This is the driving term related to the modal variables. Through this internal state group, the long-interval history integral, which originally depended on the global history convolution, can be rewritten into a local recursive form that only depends on the state of the previous time step, thus creating conditions for reducing complexity.

[0043] S7: Short-interval historical queue sliding and internal state update like Figure 3As shown, the present invention synchronously performs short interval historical queue sliding and long interval internal state update at each time step. Figure 3 Left side of the middle The region represents the long-interval historical integral, processed recursively from the ISR state; the right side... The region represents the historical integral over a short interval, processed by L1 high-precision interpolation; as time progresses from... Advance to The entire short interval window slides to the right.

[0044] When a historical state moves out of the end of a short interval window, that historical state no longer participates in subsequent short interval L1 discretization, but is incorporated into the internal state variables as a known excitation source. The update process. Subsequently, the internal state variables... The Crank-Nicolson difference scheme is used for updates, enabling local recursion to be performed while maintaining numerical stability in long-range memories. This achieves a seamless connection between high-precision representation of near history and compressed recursion of far history.

[0045] S8: PECE Integration Promotion After obtaining the L1 discrete results of the short-interval history integral and the ISR recursive results of the long-interval history integral, the PECE prediction-evaluation-correction-evaluation time-progression framework is used to solve for the modal temperature response and vibration response at the current moment. Specifically, this includes: making explicit predictions based on the known state at the previous moment; evaluating the nonlinear term, short-interval history contribution, and long-interval history contribution using the predicted values; correcting based on the evaluation results and obtaining the state at the current moment; and repeating the evaluation and correction as necessary to improve convergence accuracy.

[0046] This step ensures that the short-interval L1 results and the long-interval ISR results are integrated within a unified time-progression framework. This not only helps to handle the nonlinear terms in the thermoelastic coupling equations, but also helps to maintain the continuity of the fractional memory effect throughout the entire time history.

[0047] S9: Steady-state determination, loss extraction, and parameter scan output After the system has progressed for a sufficiently long time, the difference between key state variables in adjacent cycles can be used to determine whether a steady state has been reached. Once a steady state is reached, the energy dissipation in a single cycle within a stable oscillation period can be used to determine this. With maximum energy storage Calculate the loss factor Furthermore, it outputs the steady-state temperature results at the target frequency.

[0048] In a preferred embodiment, such as Figure 4As shown, the core numerical method of this invention can be embedded into the specific solver operation flow: first, the solver operation input is performed; then, Picard iteration is used in the inner nonlinear convergence loop, and TDMA can be used to solve the discrete equation system; then, the solution is performed according to the step size. The timeline is adjusted, and steady-state criteria are applied at vibration cycle checkpoints. Once steady-state is reached, data output and post-processing are performed to complete the process. The integral and loss factor are calculated.

[0049] Algorithm Validation and Effect Description like Figure 5 As shown, the theoretical computation time of the traditional global memory convolution algorithm increases quadratically with the number of simulation time steps, while the theoretical computation time of the linear complexity fast algorithm increases approximately linearly with the number of simulation time steps. Figure 5 The red curve in the figure represents the traditional global convolution method, and the blue curve represents the reference curve for a fast algorithm with linear complexity. This figure is used to illustrate the necessity of the interval segmentation and ISR recursion proposed in this invention.

[0050] like Figure 6 As shown, under typical operating conditions, the efficient solution method adopted in this invention maintains consistency with the reference numerical algorithm in the main trend of steady-state temperature evolution. This figure illustrates that even after introducing interval segmentation and internal state variable recursion, the solver can still effectively preserve the main characteristics of the temperature response waveform, thus providing a reliable basis for loss factor extraction.

[0051] like Figure 7 As shown, as the modal truncation order gradually increases, the dimensionless temperature response curve of the microbeam and its steady-state local magnification gradually converge, indicating that the Galerkin modal reduction step has a clear truncation convergence law. Therefore, the modal truncation order in this invention can be reasonably selected according to accuracy requirements without disrupting the overall algorithm framework.

[0052] like Figure 8 As shown, with the increase of the number of ISR integration nodes, the dimensionless temperature local response curve gradually approaches the stable reference result, indicating that the ISR orthogonal quadrature discretization of long-interval history integrals has good node convergence characteristics. This result indirectly verifies the effectiveness of long-interval history integral reconstruction and internal state variable design in S6.

[0053] like Figure 9 As shown, under the same computational task, there is a significant time difference between the global L1 algorithm and the L1-ISR algorithm proposed in this invention. The L1-ISR algorithm has better computational efficiency in long-term simulation and wideband scanning scenarios.

[0054] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not describe all details exhaustively, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification.

Claims

1. A numerical solution method for thermoelastic damping of MEMS beams based on interval-segmented infinite state representation, characterized in that, Includes the following steps: S1. Obtain the geometric parameters, material parameters, excitation parameters and boundary conditions of the target micro / nano beam resonator. Based on the Euler-Bernoulli beam theory, uniaxial stress assumption and one-dimensional non-Fourier heat conduction assumption in the thickness direction, establish the thermo-elastic coupling control equation containing fractional-order Catteneo-Vernotte heat conduction relationship and dissipation nonlinear term to characterize the thermo-mechanical coupling behavior of the resonator. S2. The thermo-elastic coupling control equations are made dimensionless by defining dimensionless temperature, dimensionless time, dimensionless thickness coordinates, dimensionless driving frequency, and dimensionless thermal relaxation time, and a single-cycle dissipation energy dissipation model is constructed. With maximum energy storage The expression for the loss factor is: ； S3. Using the Galerkin modal projection method, the dimensionless partial differential equation system is reduced to a finite-dimensional modal equation system, which includes Caputo fractional derivative terms, and the modal cutoff order is determined according to the preset accuracy. S4. On a uniform time grid, the historical integration domain involved by the Caputo fractional derivative at the current time is divided into short-interval historical integration and long-interval historical integration; wherein, the short interval covers the recent period containing the current time and the kernel function exhibits weak singularity, and the long interval covers the distant historical period far beyond the recent period and the kernel function exhibits smooth decay. S5. Discretize the short-interval historical integral using the L1 piecewise linear interpolation method to obtain the recent memory term, which represents the influence of recent history on the modal variables at the current moment. S6. The long-interval historical integral is reconstructed using an infinite state representation method, transforming the far-historical contribution term in the form of global convolution into a local recursive form based on a set of internal state variables. S7. The internal state variables are updated using the Crank-Nicolson difference scheme, and the historical information that moves out of the short interval window over time is incorporated into the internal state variables in real time through a preset short interval history queue, so as to achieve seamless connection between recent and distant historical contributions. S8. Under the PECE prediction-evaluation-correction-evaluation time-progression framework, the modal temperature response at the current moment is iteratively solved by integrating recent memory terms and recursively updated internal state variables. S9. Once the system reaches a steady state, calculate the loss factor and output the loss factor for use in the quality factor evaluation and structural design of the micro / nano beam resonator.

2. The numerical solution method for the thermoelastic damping of MEMS beams based on interval-segmented infinite state representation according to claim 1, characterized in that, The thermo-elastic coupling control equations established in S1 include fractional-order time derivative terms and thermal relaxation time constants for describing the non-instantaneous response of heat flow to temperature gradient at the micro-nano scale, as well as nonlinear heat source terms for describing the interaction between the temperature field and the strain field.

3. The numerical solution method for the thermoelastic damping of MEMS beams based on the interval-segmented infinite state representation according to claim 1, characterized in that, The specific steps in S4 are as follows: Set the current time... Historical integral domain Divided into long intervals and short interval ,in The preset short interval length is defined as the time step Δ. t Integer multiples of.

4. The numerical solution method for the thermoelastic damping of MEMS beams based on the interval-segmented infinite state representation according to claim 1, characterized in that, The specific steps in S6 are as follows: S61. The far-history convolution term is represented by a diffusion integral, and the integral is discretized into P nodes in the pseudo-frequency domain by the generalized Gauss-Laguerre orthogonal quadrature rule. S62. Define corresponding internal state variables for each quadrature node. , ,in The decay parameters for the corresponding quadrature nodes, For known driving terms related to the modal variables; The recursive solution of the first-order ordinary differential equations replaces the repeated convolution calculations for all distant histories.

5. The numerical solution method for the thermoelastic damping of MEMS beams based on the interval-segmented infinite state representation according to claim 1, characterized in that, The short interval history queue in S7 is a first-in-first-out queue. During time stepping, newly calculated recent historical states are pushed into the head of the queue. At the same time, historical state items that cross the short interval boundary and are removed from the tail of the queue are converted into known excitation sources to drive the update of the internal state variables.

6. The numerical solution method for the thermoelastic damping of MEMS beams based on the interval-segmented infinite state representation according to claim 1, characterized in that, The PECE time-progression framework in S8 specifically includes: Prediction: Based on the state at the previous time step, explicitly predict the initial state value at the current time step; Evaluation: Using the initial state value, calculate the nonlinear term of the current step, the recent memory term obtained in step S5, and the internal state variable updated in step S7; Correction: Substitute the evaluation results into the discrete equations and solve for the correction value of the state at the current time. Re-evaluation: The correction value is used to re-evaluate, in order to iteratively improve the solution accuracy of coupled nonlinear problems.

7. A design method for a micro / nano beam resonator, characterized in that, Using the method described in any one of claims 1 to 6, a wideband scanning simulation is performed on a set of resonator design schemes with different geometric or material parameters to obtain the peak value of the loss factor corresponding to each scheme. Based on the correlation between the peak value of the loss factor and the resonator quality factor, a design scheme that meets the target performance is selected.

Images