# Reduced order modelling of fully coupled electro-mechanical systems through invariant manifolds with applications to microstructures

Attilio Frangi, Alessio Colombo, Alessandra Vizzaccaro, Cyril Touzé

► To cite this version:

Attilio Frangi, Alessio Colombo, Alessandra Vizzaccaro, Cyril Touzé. Reduced order modelling of fully coupled electro-mechanical systems through invariant manifolds with applications to microstructures. International Journal for Numerical Methods in Engineering, 2024, 126 (3), pp.e7641. 10.1002/nme.7641. hal-04935025

HAL Id: hal-04935025

https://hal.science/hal-04935025v1

Submitted on 7 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Reduced order modelling of fully coupled electro-mechanical systems
through invariant manifolds with applications to microstructures

Attilio Frangi¹, Alessio Colombo¹, Alessandra Vizzaccaro², and Cyril Touzé³

¹Department of Civil and Environmental Engineering, Politecnico di Milano, P.za Leonardo da Vinci 32, 20133 Milano, Italy

²College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, United Kingdom

³Institute of Mechanical Sciences and Industrial Applications (IMSIA), ENSTA Paris - CNRS - EDF - CEA, Institut Polytechnique de Paris, 828 Boulevard des Maréchaux, 91761 Palaiseau, France

Abstract

This paper presents the first application of the direct parametrisation method for invariant manifolds to a fully coupled multiphysics problem involving the nonlinear vibrations of deformable structures subjected to an electrostatic field. The formulation proposed is intended for model order reduction of electrostatically actuated resonating Micro-Electro-Mechanical Systems (MEMS). The continuous problem is first rewritten in a manner that can be directly handled by the parametrisation method, which relies upon automated asymptotic expansions. A new mixed fully Lagrangian formulation is thus proposed which contains only explicit polynomial nonlinearities, which is then discretised in the framework of finite element procedures. Validation is performed on the classical parallel plate configuration, where different formulations using either the general framework, or an approximation of the electrostatic field due to the geometric configuration selected, are compared. Reduced-order models along these formulations are also compared to full-order simulations operated with a time integration approach. Numerical results show a remarkable performance both in terms of accuracy and wealth of nonlinear effects that can be accounted for. In particular, the transition from hardening to softening behaviour of the primary resonance while increasing the constant voltage component of the electric actuation, is recovered. Secondary resonances leading to superharmonic and parametric resonances are also investigated with the reduced-order model.

# 1 Introduction

Micro-electro-mechanical systems (MEMS) have increasingly attracted considerable interest due to their small size, high reliability, and low power consumption. MEMS accelerometers, gyroscopes, pressure sensors, micromirrors, magnetometers, microphones and many others are now essential components in many devices of everyday life [1]. The analysis of the dynamical behaviour of MEMS devices requires time-dependent, nonlinear, multiphysics models including electromagnetics, piezoelectricity, and fluid-structure interaction. Microdevice engineering faces intricate geometries and is burdened by uncertainties regarding material parameters and fabrication imperfections. This makes traditional Full-Order simulation strategies extremely expensive, if not infeasible.

In MEMS resonators, for instance, frequency is highly sensitive to the level of the oscillation. Indeed, while large actuation voltages might be desirable in a resonator to decrease the phase noise and facilitate the design of the electronic control circuit, on the other hand, they can excite non-linear phenomena, both in the mechanical structures and in the electrostatic readout. In particular, a non-linear frequency response causes a shift in the resonance frequency of the oscillator [2, 3] and a consequent loss in terms of frequency stability of the device. For this reason, it is fundamental to be able to model the dynamic behaviour of such devices by taking into account all the possible sources of non-linearities during the design process.

Apart from real-time clocks, the literature on MEMS is rich in applications associated with nonlinear effects, as recently reviewed in [4]. For instance, some MEMS like electrostatic micromirrors are intrinsically nonlinear [5] due to the large rotations of the reflecting surfaces. In [6, 7, 8], the behaviour of resonator arrays for mass sensing applications and of rate-integrating gyroscopes is addressed. Moreover, different sources of non-linearities can be combined to achieve an overall linear response of the sensor, to compensate for temperature effects, instabilities or for sensing applications [9, 10, 11]. Non-linear phenomena like bistability [12], internal resonances [13, 14, 15], self-induced parametric amplification [16, 17, 18] and frequency combs [19, 20] have also been investigated in recent years, both theoretically and experimentally.

One can identify multiple sources of nonlinearities when transformations are no longer infinitesimal: geometric effects refer to the nonlinear evolution and coupling of stress components within a structure; electrostatic forces depend in an intrinsically nonlinear manner on geometrical gaps, just like gas dissipation. The combination of geometrical and electrostatic nonlinearities has attracted many investigations like the semianalytical approach discussed by [4,
21, 22, 23, 24]. The application of Von Kármán beam theory and analytical expressions of the electrostatic forces to parallel plates capacitors allows one to simulate the classical transition from hardening to softening behaviour in the frequency response function, together with superharmonic and sub-harmonic resonances. However, the extension of these approaches to more general configurations is problematic.

In other simulations of coupled electromechanical problems [25, 26, 27, 28] the electrostatic field was computed resorting to the Boundary Element approach, with a heavy computational burden.

In recent years, several coupled electromechanical problems have been addressed in the literature using a reduction technique based on implicit condensation. This method, sometimes also denoted as the *applied force method*, relies on applying prescribed static forces to the structure to infer the nonlinear stiffness of a Reduced Order Model (ROM) [29]. It has been first applied to MEMS in [30] and then extended in [31] to simulate a MEMS gyroscope test-structure exhibiting 1:2 internal resonance, in [3] to present refined modelling of non-linearities in MEMS resonators, and in [32] to deal with coupled piezoelectric structures. However, this approach is known to have limitations that are intrinsically linked to the static underlying assumptions used in its derivation, see e.g. [33, 34, 35]. On the other hand, as explained for example in [36], nonlinear reduction techniques based on the invariant manifold theory do not have such limitations and thus should be preferred for dynamical problems.

Model order reduction techniques using nonlinear normal modes (NNM) defined as invariant manifolds of the phase space have now a long history since the first works on the subject [37, 38, 39, 40]. In recent years, the computation of invariant manifolds has been generalised to be directly applicable to large-scale systems in a series of different implementations, which will be referred to as *direct* in the remainder. Using the normal form approach [39, 40], direct computations of invariant manifolds have been proposed in [41, 42, 43] with special emphasis on using the method in a non-intrusive manner and tackling internal resonances. Based on the parametrisation method for invariant manifolds [44, 45], the uniqueness of the smoothest manifold tangent to the linear eigenspace for dissipative systems has been proven in [46] and named as Spectral submanifold (SSM) [46]. Then, a direct implementation of the parametrisation method has been proposed in [47, 48], publicly released in the code SSMtools [49]. An alternative direct implementation of the parametrisation method for invariant manifolds has also been proposed by the authors in [50, 51]. An original treatment of the non-autonomous term is also proposed in [52], that generalises previous developments used in [51, 47], and allows dealing with e.g. superharmonic resonances. Both the autonomous and non-autonomous versions of this algorithm are embedded in the open-source code MORFE (Model order reduction for finite element problems) [53]. In the following, we will refer to any *direct* implementation of the method as DPIM (direct parametrisation of invariant manifolds). The DPIM has shown great versatility and has been used to compute efficient ROMS for different applications, e.g. wing-like structure [47], MEMS structures (microbeams, micromirrors) [50, 51], rotating structures with centrifugal effect [54], viscoelastic beams with gyroscopic force [48], to name a few. For MEMS problems with weakly coupled physics, it has also been used to take into account piezoelectric actuation [55] for devices where the polarisation history is imposed.

The present article aims to extend the state of the art of direct implementations of the parametrisation method to the first truly multiphysics application, namely fully coupled electromechanical problems, and to electrostatically actuated MEMS structures in particular, which still represent the vast majority of such devices. Since the parametrisation method is best suited for polynomial nonlinearities, and thus works well with geometrical nonlinearities, the coupled electromechanical problem poses severe challenges. Indeed, the electrostatic problem is generally formulated in the current configuration where it is fully linear, and nonlinearities are *hidden* in the morphing of the mesh which is required to follow the structure deformation. To make them explicit, the problem can be rewritten in Lagrangian form on the initial configuration, see e.g. [25, 26, 28] where a Boundary Element approach has been used for the solution of the electrostatic field. However, this step introduces non-polynomial nonlinearities. Since our goal here is to propose an efficient reduction technique compatible with the parametrisation method, the equations of motion are reworked to express the semi-discretised finite element (FE) problem as differential-algebraic equations (DAE) containing only polynomial nonlinearities, by introducing a new field associated with the gradient of the electric scalar potential. The resulting mixed formulation is new and original to the best of the authors’ knowledge. Although this new formulation might not be more efficient than existing ones for simulating the full order model (FOM), its only purpose is to produce rigorous ROMs, which will be orders of magnitude faster than the FOM, regardless of the formulation. Since the parametrisation method can handle systems formulated as a DAE without important modification to the algorithm, as proven in [56] and subsequently used in [52, 57], thanks to the mixed formulation proposed in this paper, an efficient and accurate ROM strategy will be derived for arbitrary order solutions of forced electromechanical systems, and seamlessly embedded in the DPIM implementation.

The paper is organised as follows. Section 2 presents the general formalism for a generic electromechanical problem involving a deformable body and rigid electrodes. Particular attention is devoted to rewriting such a problem in the original configuration and introducing an auxiliary field in order to express the nonlinearities under polynomial form. Section 3 discusses the FE discretisation and summarises the semi-discrete problem to solve. The DPIM procedure is then presented starting from the semi-discrete DAE in Section 4. The configuration where the structure and the
electrode are two parallel plates is then addressed in Section 5, and the resulting simplifications that can be used to express the electrostatic forces in such case are recalled. Finally, numerical results for a two-dimensional parallel plate problem are collected in Section 6. Pull-in verifications, primary resonance solution with hardening/softening behaviour, and secondary resonances are investigated and comparisons between full-order simulations and ROMs are reported, showing the very good accuracy provided by the reduction technique.

## 2 Problem formulation

In this section, the equations of motion for the coupled nonlinear electromechanical problem are derived, and a particular emphasis is put on rewriting them in such a way that they can be readily solved by the direct parametrisation method for invariant manifold. This leads to adding new unknown fields in order to express the semi-discrete problem resulting from the finite element discretisation as a differential-algebraic equation (DAE) with only polynomial nonlinearities.

Let us consider a deformable solid, of initial configuration $\Omega_0$ as shown in Fig. 1. The body is subjected to an electric field due to the interaction with a set of rigid electrodes $\Omega_K$ and is otherwise immersed in the infinite empty space $A_0$. The solid is considered a perfect conductor, a condition which is very well respected in silicon microsystems, with known given potential $\varphi^D = V_{DC} + V_{AC} \sin \omega t$, where the constant part $V_{DC}$ and the alternating amplitude $V_{AC}$ are such that $V_{AC} \ll V_{DC}$. A uniform voltage $\varphi_K^D$ is imposed on each $\Omega_K$ (equipotential body) and it is here assumed, for simplicity, that $\varphi_K^D = 0$, although more general configurations can be easily considered.

To simplify the notation, $\varphi^D$ denotes the given voltage on $\partial\Omega_0 \cup S_0$, where $S_0$ is the collection of all the (fixed) surfaces of the electrodes.

<figure><img src="image_2.png" /><figcaption>Figure 1: Example of the layout considered for the electromechanical problem, shown in the initial configuration.</figcaption></figure>

The problem is classically formulated in terms of a potential function $\varphi$ such that the electric field is $-\text{grad}\varphi$, and $\varphi$ is governed by the following set of equations written in the current configuration:

$$ \text{div} \text{grad} \varphi = 0, \quad \text{in } A, \qquad (1) $$

$$ \varphi = \varphi^D, \quad \text{on } S_0, \qquad (2) $$

$$ \text{grad}\varphi \cdot \hat{n} = 0, \quad \text{on } S_\infty, \qquad (3) $$

where div and grad denote operators with respect to the current coordinates $\mathbf{y}$, A is the current air domain, $S_\infty$ is a sufficiently distant truncation boundary and the vector $\hat{n}$ is the unit normal vector oriented from the solid into the air. On the surface of a conductor, grad$\varphi$ is orthogonal to the surface. Alternatively, the problem can be formulated
in a weak manner as follows.

Find $\varphi \in C_{\varphi}(\varphi^{\mathrm{D}})$ such that:

$$ \int_A \mathrm{grad}\tilde{\varphi} \cdot \mathrm{grad}\varphi \, d\Omega = 0, \quad \forall \tilde{\varphi} \in C_{\varphi}(0), \qquad (4) $$

where $C_{\varphi}(\varphi^{\mathrm{D}})$ denotes the space of sufficiently continuous functions respecting the Dirichlet boundary conditions given by Eq. (2). The solution to Eq. (4) requires to mesh the air domain and to truncate it at $S_{\infty}$. Moreover, it is worth stressing that the domain A is evolving according to the deformation of the solid.

The electrostatic pressure (force per unit surface), exerted on the surface of the deformable body at an arbitrary point y reads:

$$ \hat{\mathbf{b}} = \frac{\hat{\sigma}^2}{2\varepsilon_0} \hat{\mathbf{n}}, \quad \hat{\sigma} = -\varepsilon_0 \mathrm{grad}\varphi \cdot \hat{\mathbf{n}}, \qquad (5) $$

where $\hat{\sigma}$ is the surface charge on a conductor (Coulomb theorem).

The electrostatic forces induce vibrations in the structure and the dynamical equilibrium is governed by the Principle of Virtual Power (PVP) enforced in the current unknown configuration:

$$ \int_{\Omega} \rho \ddot{\mathbf{u}} \cdot \mathbf{u} \, d\Omega + \int_{\Omega} \boldsymbol{\sigma} : \mathrm{grad} \tilde{\mathbf{u}} \, d\Omega = \int_{\partial \Omega} \hat{\mathbf{b}} \cdot \tilde{\mathbf{u}} \, dS, \quad \forall \tilde{\mathbf{u}} \in C_u(0), \qquad (6) $$

where $\boldsymbol{\sigma}$ is the Cauchy stress tensor and $\hat{\mathbf{b}}$ are surface forces in the current configuration. The test field $\tilde{\mathbf{u}}$ is defined over the space $C_u(0)$ of admissible functions, i.e. functions that are vanishing where Dirichlet boundary conditions are prescribed. The first integral in Eq.(6) expresses the power $P_a$ of inertia forces, while the second term is the opposite of the power of internal forces $P_i$, and the right-hand side is the power of external forces $P_e$, such that the PVP can also be simply rewritten as:

$$ P_a - P_i = P_e, \quad \forall \tilde{\mathbf{u}} \in C_u(0). \qquad (7) $$

It is important to remark that with this continuous weak formulation given by Eqs. (4) and (6), the nonlinearities of the problem, stemming from both the mechanics (geometric nonlinearity) and the electric coupling, do not explicitly appear, being hidden in the fact that the current configuration is unknown. A classical method to make them explicit is to rewrite the problem in the original configuration.

## 2.1 Formulation in the original configuration

The problem can be rewritten in the original configuration, using classical developments in continuum mechanics. For the sake of brevity, only the main results are reported here and the interested reader can find details of this derivation in appendix A. The PVP (6) in the original configuration reads:

$$ \int_{\Omega_0} \rho_0 \ddot{\mathbf{u}} \cdot \tilde{\mathbf{u}} \, d\Omega_0 + \int_{\Omega_0} \boldsymbol{\sigma}^{\mathrm{PK}}[\mathbf{u}] : (\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}]) \, d\Omega_0 = \int_{\partial \Omega_0} \mathbf{b} \cdot \tilde{\mathbf{u}} \, dS_0 \quad \forall \tilde{\mathbf{u}} \in C_u(0), \qquad (8) $$

where $\varepsilon$ and $\bar{\nabla}$, defined in Appendix A, stand for:

$$ \varepsilon[\tilde{\mathbf{u}}] = \frac{1}{2}(\nabla\tilde{\mathbf{u}} + \nabla^T\tilde{\mathbf{u}}) \quad \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}] = \frac{1}{2}(\nabla^T\mathbf{u} \cdot \nabla\tilde{\mathbf{u}} + \nabla^T\tilde{\mathbf{u}} \cdot \nabla\mathbf{u}), $$

$\nabla$ denotes the gradient with respect to the initial coordinates $\mathbf{x}$, and $\mathbf{y} = \mathbf{x} + \mathbf{u}$. In Eq. (8), $\boldsymbol{\sigma}^{\mathrm{PK}}$ denotes the second Piola-Kirchhoff stress tensor. Since for the class of problems at hand strains are small even if transformations are non-infinitesimal, the Kirchhoff constitutive law holds, yielding:

$$ \boldsymbol{\sigma}^{\mathrm{PK}} = \mathcal{A} : \mathbf{e} = \mathbf{e} : \mathcal{A}, \qquad (9) $$

where $\mathcal{A}$ is the tensor of elastic constants and $\mathbf{e}$ is the Green-Lagrange strain tensor:

$$ \mathbf{e}[\mathbf{u}] = \frac{1}{2}(\nabla\mathbf{u} + \nabla^T\mathbf{u} + \nabla^T\mathbf{u} \cdot \nabla\mathbf{u}). \qquad (10) $$

Let us now detail the derivation of the electrostatic force in the initial configuration. Since it is a follower force, the calculations are not straightforward and have been thus reported in Appendix C. In particular, using Eq. (78), one has:

$$ \begin{aligned} P_e &= \int_{\partial\Omega} \frac{\hat{\sigma}^2}{2\varepsilon_0} (\tilde{\mathbf{u}} \cdot \hat{\mathbf{n}}) dS \\ &= \frac{\varepsilon_0}{2} \int_{\partial\Omega_0} \sigma^2 \left( \mathbf{n} + \frac{1}{J_S} (\mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2}) \right) \cdot \tilde{\mathbf{u}} dS_0, \end{aligned} \qquad (11) $$
where $J_S$ is the surface jacobian, see Eq.(71), $\mathbf{u}_{,\alpha}$ denotes the derivative with respect to the surface coordinate $a_\alpha$ used to represent the FE surface elements, $\mathbf{g}_\alpha$ are surface covariant vectors and $\sigma(\mathbf{x}) = \hat{\sigma}(\mathbf{y}(\mathbf{x}))$. Eq. (11) indirectly provides an explicit expression for $\mathbf{b}$ in Eq. (8):

$$ \mathbf{b} = \sigma^2 \left( \mathbf{n} + \frac{1}{J_S} (\mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2}) \right). $$

In order to propose a unified formulation expressed in the initial configuration, the electrostatics equilibrium, Eq. (4), also needs to be rewritten, a step that is rarely done in the literature and is detailed next. One of the problems to be solved is linked to the fact that the deformation of the air domain A follows the displacement of the structure $\mathbf{u}$. Since this displacement is defined only for the structural nodes, the air domain deformation is strictly speaking only defined on $\partial\Omega$. One method to cope with this issue is to extend smoothly the structural displacement vector $\mathbf{u}$ to A to represent the deformation of the air domain. This extension, denoted $\mathbf{s}$, is indeed fictitious and not unique. In a FEM sense, one could a priori set the nodal values to zero in the air domain immediately outside the conductors, but this leads to strong gradients, thus generating a poor solution quality. The option adopted in the present investigation is to solve a Dirichlet problem for each displacement component, i.e. find $s_i \in C_s(u_i)$, $\forall i \in [1, d]$:

$$ \int_{A_0} \nabla \tilde{s} \cdot \nabla s_i \, d\Omega_0 = 0 \quad \forall \tilde{s} \in C_s(0), \qquad (12) $$

where the space $C_s(u_i)$ collects the functions such that $s_i = u_i$ on $\partial\Omega_0$, $s_i = 0$ on $S_0 \cup S_\infty$ and $d$ is the dimension of the problem, $d = 2$ if the problem is 2-dimensional, $d = 3$ if it is 3-dimensional. An example of this extension procedure is presented in Section 6, where a clamped-clamped beam faces a rigid electrode. Other choices that would be computationally more efficient will be commented on in Section 2.3.1. Once $\mathbf{s}$ has been defined, the electrical equilibrium in the initial configuration reads

$$ \int_{A_0} \nabla \tilde{\varphi} \cdot \mathbf{c}^{-1}[\mathbf{s}] \cdot \nabla \varphi J d\Omega_0 = 0, \quad \forall \tilde{\varphi} \in C_\varphi(0), \qquad (13) $$

where $J = \det \mathbf{f}[\mathbf{s}]$ is the Jacobian of the transformation:

$$ \mathbf{f}[\mathbf{s}] = \mathbf{1} + \nabla \mathbf{s}, \qquad (14) $$

$$ \mathbf{c}[\mathbf{s}] = \mathbf{f}^{\mathrm{T}}[\mathbf{s}] \cdot \mathbf{f}[\mathbf{s}] = \mathbf{1} + \nabla \mathbf{s} + \nabla^{\mathrm{T}} \mathbf{s} + \nabla^{\mathrm{T}} \mathbf{s} \cdot \nabla \mathbf{s}. \qquad (15) $$

The main concern with this formulation is that $\mathbf{c}^{-1}$ is not polynomial in $\mathbf{s}$. The next section aims to circumvent this difficulty by adopting a new formulation of the continuous problem.

## 2.2 Mixed formulation with polynomial nonlinearities

The goal of this Section is to show how the introduction of an additional field allows the expression of the problem to be solved with polynomial nonlinearities only. To remove the non-polynomial terms from Eq. (13), an additional vector field $\psi$ can be introduced as:

$$ \psi = \nabla \varphi \cdot \mathbf{f}^{-1}[\mathbf{s}] = \mathbf{f}^{-\mathrm{T}}[\mathbf{s}] \cdot \nabla \varphi, \quad \nabla \varphi = \mathbf{f}^{\mathrm{T}}[\mathbf{s}] \cdot \psi = \psi \cdot \mathbf{f}[\mathbf{s}]. \qquad (16) $$

Remark that $\psi$ corresponds to the gradient in actual coordinates and is strictly related to the electric field. The space of $\psi$ might be for instance that of piecewise (discontinuous) linear vectors (element by element). We enforce Eq.(16) in a weak form, using as test functions all possible $\tilde{\psi} \in C_\psi$:

$$ \int_{A_0} \tilde{\psi} \cdot (\mathbf{f}^{\mathrm{T}}[\mathbf{s}] \cdot \psi - \nabla \varphi) d\Omega_0 = 0, \quad \forall \tilde{\psi} \in C_\psi, \qquad (17) $$

while Eq. (13) now becomes:

$$ \int_{A_0} \nabla \tilde{\varphi} \cdot (J \mathbf{f}^{-1}[\mathbf{s}]) \cdot \psi \, d\Omega_0 = 0, \quad \forall \tilde{\varphi} \in C_\varphi(0). \qquad (18) $$

As $J\mathbf{f}^{-1}$ is the adjoint of $\mathbf{f}$, i.e. the transpose of the cofactor matrix of $\mathbf{f}$, we will use the notation $\mathrm{adj}(\mathbf{f})$ in what follows. Recall that the coefficient $(i, j)$ of the cofactor matrix of $\mathbf{f}$ is $(-1)^{i+j} M_{ij}$, where $M_{ij}$ is the minor of the coefficient. Hence, each coefficient of $\mathrm{adj}(\mathbf{f})[\mathbf{s}]$ is a quadratic polynomial in $\mathbf{s}$ in a 3D space. On the other hand, for 2D space problems, each coefficient of $\mathrm{adj}(\mathbf{f})[\mathbf{s}]$ is linear in $\mathbf{s}$. Such a simplification is helpful in solving planar problems and will be used in Section 6.
Since on the surface of the conductors grad$\varphi$ is parallel to $\hat{n}$, the following relationships hold:

$$ \sigma \hat{n} = -\epsilon_0 \nabla \varphi \cdot \mathbf{f}^{-1}, \quad \sigma^2 = \epsilon_0^2 \|\psi\|^2, $$

and the virtual power of the electrostatic pressure in the reference configuration can be written as a function of the additional vector field $\psi$ in a formulation without non-polynomial nonlinearities as:

$$ P_e = \frac{\epsilon_0}{2} \int_{\partial\Omega_0} \|\psi\|^2 \left( \mathbf{n} + \frac{1}{J_S} (\mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2}) \right) \cdot \tilde{\mathbf{u}} \, dS_0. \quad (19) $$

## 2.3 Summary of the General Formulation

To summarise, we have proposed a formulation of the coupled electromechanical problem consisting of a system of four equations which can be expressed as follows.

Find $\mathbf{u} \in C_u(0)$, $\mathbf{s} \in C_s(\mathbf{u})$, $\varphi \in C_\varphi(\varphi^D)$, $\psi \in C_\psi$ such that:

* the mechanical equilibrium:

$$ \int_{\Omega_0} \rho_0 \ddot{\mathbf{u}} \cdot \tilde{\mathbf{u}} \, d\Omega_0 + \int_{\Omega_0} \sigma^{\text{PK}}[\mathbf{u}] : (\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}]) \, d\Omega_0 \\ = \frac{\varepsilon_0}{2} \int_{\partial\Omega_0} \|\psi\|^2 \left( \mathbf{n} + \frac{1}{J_S} (\mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2}) \right) \cdot \tilde{\mathbf{u}} \, dS_0, \quad \forall \tilde{\mathbf{u}} \in C_u(0), \quad (20) $$

* the extension problem, $\forall i \in [1, d]$:

$$ \int_{A_0} \nabla \tilde{s} \cdot \nabla s_i \, d\Omega_0 = 0, \quad \forall \tilde{s} \in C_s(u_i), \quad (21) $$

* the compatibility equation:

$$ \int_{A_0} \tilde{\psi} \cdot (\psi + \nabla^T \mathbf{s} \cdot \psi - \nabla \varphi) d\Omega_0 = 0, \quad \forall \tilde{\psi} \in C_\psi, \quad (22) $$

* the electrostatics equilibrium:

$$ \int_{A_0} \nabla \tilde{\varphi} \cdot (\psi + \operatorname{adj}(\nabla \mathbf{s}) \cdot \psi) \, d\Omega_0 = 0, \quad \forall \tilde{\varphi} \in C_\varphi(0), \quad (23) $$

are all satisfied.

The system of equations (20)-(23) will be referred to as the Mixed General Formulation (MGF) in the rest of the paper, as it mixes standard kinematic unknowns like displacements and new force variables like $\psi$. The MGF falls in the category of mixed approaches that are vastly debated in the literature, see e.g. [58].

It provides a very general framework to solve for coupled electromechanical problems, which also fits the requirements of the parametrisation method for invariant manifolds [45]. Indeed, the problem with non-polynomial nonlinearities has been replaced by a DAE containing only quadratic and cubic terms, at the price of introducing an auxiliary field, which results in an added compatibility equation. Section 4 will show how the general reduction strategy provided by the parametrisation method can be applied to the semi-discrete problem resulting from a FE discretisation of these equations.

### 2.3.1 Remark on computational efficiency

From the inspection of the MGF it is clear that the introduction of the two fields $\mathbf{s}$ and $\psi$ is a considerable burden from the computational point of view. However, this can be drastically reduced by defining, in the original configuration, a boundary layer $L_0$ of suitable thickness around the deformable structure such that $\mathbf{s}$ decays smoothly to zero from $\partial\Omega_0$ to the outer boundary $\partial L_0$ of the layer. In $A_0 \setminus L_0$, $\mathbf{s}$ vanishes, $\mathbf{c}[\mathbf{s}] = \mathbf{1}$, there is no need to introduce $\psi$ and the original trivial electrostatic formulation is recovered.
### 2.3.2 Simplifications in 2D

Since the results analysed in Section 6 are obtained using a 2D implementation, the simplifications pertaining to this specific case are here detailed. Solids are extruded along the z direction, the first surface coordinate $a_1$ coincides with z and hence $\mathbf{g}_1 = \mathbf{e}_z$ and $\mathbf{u}_z = \mathbf{0}$. The outer surfaces of the solids are represented in 2D by the lines which are their intersection with the $z = 0$ plane. The virtual power $\mathcal{P}_e$ of electrostatic forces hence becomes:

$$ \mathcal{P}_e = \frac{\varepsilon_0}{2} \int_{\partial\Omega_0} \|\psi\|^2 \left( \mathbf{n} + \frac{\mathbf{e}_z \wedge \mathbf{u}_{,a}}{J_S} \right) \cdot \tilde{\mathbf{u}} \, ds, \quad (24) $$

where $a$ is an abscissa running along the line-boundary $\partial\Omega_0$, $\mathbf{u}_{,a}$ is the directional derivative of $\mathbf{u}$ with respect to $a$ and $ds = J_S da$. Moreover, $\mathrm{adj}(\nabla s)$ is only linear in $s$.

## 3 Solution procedure

In coupled electromechanical problems, vibrations develop around an initial static solution associated with the presence of the constant voltage $V_{DC}$, according to the definition of the applied voltage given in Section 2. This is known to have important effects on the nonlinear response of the structure, as will be evidenced in Section 6.

In the first step, the static solution is computed using a Newton-Raphson solver linearising the governing equations around the current guess of the solution and generating the so-called *tangent stiffness matrix* at every iteration. It is worth stressing that this matrix, at convergence, will also represent the linear part of the expansion in the DPIM procedure, see Section 4, that, together with the mass matrix, defines the eigenvalue problem.

We will henceforth denote with $\mathbf{u}_0, \mathbf{s}_0, \psi_0, \varphi_0$ the static solution and with $\mathbf{u}_0 + \mathbf{u}, \mathbf{s}_0 + \mathbf{s}, \psi_0 + \psi, \varphi_0 + \varphi$ the total solution without changing notation, for the sake of simplicity.

In the next section, we will inject this decomposition into the governing equations and partition all the terms into linear, quadratic and cubic contributions, in order to prepare the path for the discretisation and subsequent application of the DPIM procedure to obtain a reduced-order model. Constant terms only define the static solution and will be disregarded.

### 3.1 FEM discretisation

The discretised version of the equations is here provided according to standard FE procedures, by focusing on two-dimensional planar problems, for the sake of simplicity. Note that there is no difference between the 2D and 3D problems from the theoretical point of view, nevertheless, to obtain a rapid proof of concept of the reduction method, the presentation and the results shown in Section 6 will only consider 2D problems. All constant terms associated with the initial static solution are removed. In the present implementation, $\mathbf{u}$ and $\mathbf{s}$ are interpolated with quadratic triangles, as well as $\varphi$. On the contrary, $\psi$ is modelled as piecewise-linear, discontinuous from triangle to triangle. This ensures convergence at least in the linear case [58].

**Mechanical equilibrium** We analyse here the terms associated with the virtual power of electrostatic forces in Eq. (20), while we simply provide the final result for the more classical mechanical terms, since they have already been considered elsewhere, see e.g. [55], and are briefly discussed in Appendix A. If $\mathbf{n}_0$ denotes the normal to the surface deformed by the static solution $\mathbf{u}_0$, and $J_{S0}$ is the corresponding surface Jacobian,

$$ \mathcal{P}_e = \frac{\varepsilon_0}{2} \underbrace{\int_{\partial\Omega_0} \tilde{\mathbf{u}} \cdot \left( 2\mathbf{n}_0(\psi_0 \cdot \psi) + \|\psi_0\|^2 \frac{\mathbf{e}_z \wedge \mathbf{u}_{,a}}{J_{S0}} \right) ds}_{\text{linear}} + \underbrace{\frac{\varepsilon_0}{2} \int_{\partial\Omega_0} \tilde{\mathbf{u}} \cdot \left( \mathbf{n}_0 \|\psi\|^2 + 2 \frac{\mathbf{e}_z \wedge \mathbf{u}_{,a}}{J_{S0}} (\psi_0 \cdot \psi) \right) ds}_{\text{quadratic}} + \underbrace{\frac{\varepsilon_0}{2} \int_{\partial\Omega_0} \tilde{\mathbf{u}} \cdot \left( \frac{\mathbf{e}_z \wedge \mathbf{u}_{,a}}{J_{S0}} \|\psi\|^2 ds, \right)}_{\text{cubic}} \quad (25) $$

which generates the discretised contributions:

$$ \mathcal{P}_e = \tilde{\mathbf{U}}^{\mathrm{T}} (\mathbf{R}_u \mathbf{U} + \mathbf{R}_{\psi} \mathbf{\Psi} + \mathbf{R}_{\psi\psi} (\mathbf{\Psi}, \mathbf{\Psi}) + \mathbf{R}_{\psi u} (\mathbf{\Psi}, \mathbf{U}) + \mathbf{R}_{\psi\psi u} (\mathbf{\Psi}, \mathbf{\Psi}, \mathbf{U})). \quad (26) $$

Here $\mathbf{U}$ and $\tilde{\mathbf{U}}$ are vectors collecting the unknown and test nodal values, respectively, for the $\mathbf{u}$ and $\tilde{\mathbf{u}}$ fields. Similarly, $\mathbf{\Psi}$ collects the nodal values of the $\psi$ discretisation, element by element. Other upright bold symbols, like $\mathbf{K}$, denote
matrices, while italic bold symbols, like **R**, denote quadratic or cubic nonlinear operators. When added to the $\mathcal{P}_a$ and $\mathcal{P}_i$ terms, the final semi-discrete equation reads:

$$ M\ddot{U} + (\alpha M + \beta K)\dot{U} + (K - R_u)U - R_\psi\Psi + G(U, U) + H(U, U, U) - R_{\psi\psi}(\Psi, \Psi) - R_{\psi u}(\Psi, U) - R_{\psi\psi u}(\Psi, \Psi, U) = 0, \quad (27) $$

where a mechanical dissipation term has been added in the form of Rayleigh damping.

**Extension equation** The extension equation is linear and is the easiest one to discretise. It is written for all the components $s_i$ at once by collecting in **S** the unknown nodal values of the interpolation for **s**:

$$ E_s S + E_u U = 0, \quad (28) $$

where $E_s$ is a square positive definite matrix and $E_u U$ accounts for the Dirichlet boundary conditions imposing that $\mathbf{s} = \mathbf{u}$ on $\partial\Omega_0$.

**Compatibility equation** The compatibility equation, setting $f_0 = 1 + \nabla s_0$, becomes:

$$ \underbrace{\int_{A_0} \tilde{\psi} \cdot (f_0^T \cdot \psi + \nabla^T s \cdot \psi_0 - \nabla \varphi) d\Omega_0}_{\text{linear}} + \underbrace{\int_{A_0} \tilde{\psi} \cdot \nabla^T s \cdot \psi d\Omega_0}_{\text{quadratic}} = 0, \quad (29) $$

which generates the discretised form:

$$ C_u U + C_s S + C_\psi \Psi + C_\varphi \Phi + C_{\psi s}(\Psi, S) = C_V V_{AC} \sin(\omega t), \quad (30) $$

where $\Phi$ is the collection of the unknown nodal values for the $\varphi$ field. It is worth remarking that $C_u U$ stems from the Dirichlet boundary conditions on $\mathbf{s}$, while $C_V$ is due to the imposed voltage which fixes the nodal values of $\varphi$ on $\partial\Omega_0$. It is also important to stress that $V_{DC}$ appears in the above equations only through the static solution which has been computed as a first step.

**Electrostatics equilibrium** The last equation of the problem concerns the electrostatics equilibrium and reads:

$$ \underbrace{\int_{A_0} \nabla \tilde{\varphi} \cdot (\operatorname{adj}(f_0) \cdot \psi + \operatorname{adj}(\nabla s) \cdot \psi_0) d\Omega_0}_{\text{linear}} + \underbrace{\int_{A_0} \nabla \tilde{\varphi} \cdot \operatorname{adj}(\nabla s) \cdot \psi d\Omega_0}_{\text{quadratic}} = 0, \quad (31) $$

which then generates the discretised form:

$$ D_u U + D_s S + D_\psi \Psi + D_{\psi s}(\Psi, S) = 0. \quad (32) $$

# 4 Model order reduction

Now that the semi-discretised problem has been derived, the DPIM procedure, i.e. a direct implementation of the parametrisation method for invariant manifold, can be applied. This nonlinear reduction technique is a simulation-free method that operates directly from the physical space, and computes arbitrary order nonlinear mappings and reduced dynamics along the selected invariant manifold. It relies on a general procedure explained for instance in [45], which has been specified for mechanical problems for example in [46, 59, 47, 50, 51, 52]. To apply the DPIM algorithm by strictly following the procedure proposed in [52], the first step consists of writing the system as first-order DAE. To cope with the second-order time derivative in Eq. (27), an additional velocity field **V** is introduced via $M\ddot{U} = MV$. Combining this additional equation with Eq. (27), Eq. (28), Eq. (30), and Eq. (32), one obtains:

$$ M\dot{V} = -(\alpha M + \beta K)V + (-K + R_u)U + R_\psi\Psi - G(U, U) - H(U, U, U) + R_{\psi\psi}(\Psi, \Psi) + R_{\psi u}(\Psi, U) + R_{\psi\psi u}(\Psi, \Psi, U) \quad (33a) $$

$$ M\ddot{U} = MV \quad (33b) $$

$$ 0 = -E_s S - E_u U \quad (33c) $$

$$ 0 = -C_u U - C_s S - C_\psi \Psi - C_\varphi \Phi - C_{\psi s}(\Psi, S) + C_V V_{AC} \sin(\omega t) \quad (33d) $$

$$ 0 = -D_u U - D_s S - D_\psi \Psi - D_{\psi s}(\Psi, S), \quad (33e) $$
where the derivative terms have been brought to the left-hand side and the algebraic terms on the right-hand side of the equations.

We can collect all the variables of the problem into a single variable y defined as:

$$ y = \begin{bmatrix} V \\ U \\ S \\ \Psi \\ \Phi \end{bmatrix} \qquad (34) $$

and one gets a compact form for the problem in the sole y variable with nonlinearities up to cubic order, which reads:

$$ B\dot{y} = Ay + Q(y, y) + T(y, y, y) + \varepsilon(C_{+}e^{+i\omega t} + C_{-}e^{-i\omega t}). \quad (35) $$

In particular, the linear terms of the system are collected by means of the matrices A and B, defined as:

$$ B \doteq \begin{bmatrix} M & 0 & 0 & 0 & 0 \\ 0 & M & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix}, \quad A \doteq \begin{bmatrix} -(\alpha M + \beta K) & -K + R_u & 0 & R_\psi & 0 \\ M & 0 & 0 & 0 & 0 \\ 0 & -E_u & -E_s & 0 & 0 \\ 0 & -C_u & -C_s & -C_\psi & -C_\phi \\ 0 & -D_u & -D_s & -D_\psi & 0 \end{bmatrix}. \quad (36) $$

The polynomial nonlinearities are collected into two vectors separating the quadratic forces Q from the cubic ones T, which are defined as:

$$ Q(y, y) \doteq \begin{bmatrix} -G(U, U) + R_{\psi\psi}(\Psi, \Psi) + R_{\psi u}(\Psi, U) \\ 0 \\ 0 \\ -C_{\psi s}(\Psi, S) \\ -D_{\psi s}(\Psi, S) \end{bmatrix}, \quad T(y, y, y) \doteq \begin{bmatrix} -H(U, U, U) + R_{\psi\psi u}(\Psi, \Psi, U) \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \quad (37) $$

Lastly, the non-autonomous term has been defined as:

$$ \varepsilon \mathbf{C}_{\pm} \doteq \pm \frac{1}{2i} \begin{bmatrix} 0 \\ 0 \\ 0 \\ \mathbf{C}_V \\ 0 \end{bmatrix} V_{AC}, \quad \text{such that:} \quad \varepsilon (\mathbf{C}_+ e^{+i\omega t} + \mathbf{C}_- e^{-i\omega t}) = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \mathbf{C}_V \\ 0 \end{bmatrix} V_{AC} \sin(\omega t). \quad (38) $$

Note that the meaning of ε coincides with the amplitude of the non-autonomous term. Consequently, C_± has to be unitary in some sense. To meet this constraint, C_± will be normalised with respect to one of the eigenvectors of the system.

The compact form of the problem around the fixed point given by Eq. (35) is the starting point for the non-autonomous DPIM algorithm introduced in [52]. The autonomous version of the DPIM procedure [50], following a general strategy exposed for instance in [45], seeks a parametrisation of the smoothest invariant manifold (i.e. the spectral submanifold as shown in [46]) tangent to a given linear subspace, whereas the non-autonomous one [52] takes into account the time dependence of said manifolds under the action of the forcing term by seeking a parametrisation of the perturbed whisker [60, 61, 51]. The selected linear subspace is called the *master* subspace, and it is obtained from the eigenvectors of the system linearised around the fixed point. The first step of the DPIM procedure is then to compute the solution of the eigenproblem for a complex conjugate pair (or more) that will define the master subspace.

Let us define the triplet {$\lambda_k, \mathbf{Y}_k, \mathbf{X}_k$} to be one eigenvalue, right eigenvector, and conjugate left eigenvector of the system respectively. They fulfil the following properties:

$$ (\lambda_k \mathbf{B} - \mathbf{A}) \mathbf{Y}_k = \mathbf{0} \qquad (39a) $$

$$ (\lambda_k \mathbf{B}^T - \mathbf{A}^T) \mathbf{X}_k = \mathbf{0}. \qquad (39b) $$

The normalisation condition of the eigenvectors reads:

$$ \mathbf{X}_k^T \mathbf{B} \mathbf{Y}_k = b_k \delta_{kl}, \quad \forall k, l \qquad (40) $$
Therefore, the normalisation condition for the forcing reads:

$$ ||\mathbf{X}_k^T \mathbf{C}_\pm||_2 = b_k. \qquad (41) $$

The master subspace is typically selected to be the subspace spanned by the slowest modes of the system, i.e. the eigenvectors related to the eigenvalues with the smallest real part. In the non-autonomous context, we are interested in reducing the system around resonance. The forcing frequency $\omega$ will be close to a complex conjugate pair of eigenvalues of the system, thus exciting one of the modes: $\lambda_k \approx +i\omega$ and $\bar{\lambda}_k \approx -i\omega$, as the problem presented is lightly damped, therefore the real part of the eigenvalues is much smaller than the imaginary one.

We seek a parametrisation of the invariant manifold with respect to new normal coordinates $\mathbf{z}$ in the form:

$$ \mathbf{W}(\mathbf{z}) = \sum_{p=1}^{o} \sum_{k=1}^{m_p} \mathbf{W}^{(p,k)} \mathbf{z}^{\alpha(p,k)}, \qquad (42) $$

on the manifold, the dynamics is also described by a polynomial expansion as:

$$ \dot{\mathbf{z}} = \mathbf{f}(\mathbf{z}) = \sum_{p=1}^{o} \sum_{k=1}^{m_p} \mathbf{f}^{(p,k)} \mathbf{z}^{\alpha(p,k)}, \qquad (43) $$

The monomial associated to $\alpha(p, k)$ simply reads $\mathbf{z}^{\alpha(p,k)} = z_1^{\alpha_1} z_2^{\alpha_2} \dots z_d^{\alpha_d}$.

In the context of a non-autonomous problem, the normal coordinates are of two different types: the coordinates parametrising the autonomous manifold $\bar{\mathbf{z}}$, and those accounting for the time dependence of the invariant manifold $\tilde{\mathbf{z}}$:

$$ \mathbf{z} = \begin{bmatrix} \bar{\mathbf{z}} \\ \tilde{\mathbf{z}} \end{bmatrix}. \qquad (44) $$

In particular, the first set of variables $\bar{\mathbf{z}}$ represent an identity-tangent transformation of the modal coordinates of the selected master subspace. Typically, the forcing is resonant with only one mode, so a single complex conjugate pair of $\bar{\mathbf{z}}$ coordinates is enough to reproduce the system's behaviour. Regarding the second set of variables $\tilde{\mathbf{z}}$, they are book-keeping coordinates used in the implementation to automatically expand in the $\varepsilon$ parameter, thus avoiding first-order truncation on the amplitude of the forcing [47, 51], see [52] for a more detailed presentation. They are defined as:

$$ \tilde{\mathbf{z}} = \begin{bmatrix} \varepsilon e^{+i\omega t} \\ \varepsilon e^{-i\omega t} \end{bmatrix}. \qquad (45) $$

When the change of coordinates $\mathbf{W}(\mathbf{z})$ is inserted into the original equation in place of $\mathbf{y}$, the invariance equation is obtained:

$$ \mathbf{B} \nabla \mathbf{W}(\mathbf{z}) \mathbf{f}(\mathbf{z}) = \mathbf{A} \mathbf{W}(\mathbf{z}) + \mathbf{Q}(\mathbf{W}(\mathbf{z}), \mathbf{W}(\mathbf{z})) + \mathbf{T}(\mathbf{W}(\mathbf{z}), \mathbf{W}(\mathbf{z}), \mathbf{W}(\mathbf{z})) + [\mathbf{C}_+ \quad \mathbf{C}_-] \tilde{\mathbf{z}}, \qquad (46) $$

where the first term comes from the chain rule applied to $\dot{\mathbf{y}}$.

The invariance equation is solved recursively, order by order. At order one, the linear problem is retrieved, while higher orders from $p = 2$ give rise to the so-called homological equation of order $p$. These equations can be written at the level of an arbitrary monomial, resulting in an underdetermined, ill-conditioned problem. The underdeterminacy comes from the fact that two unknowns are input, the nonlinear mapping and the reduced dynamics. This problem is cured by introducing different styles of solution, the most important ones being the graph style and the normal form style [45]. The ill-conditioning comes from the nonlinear resonance relationships. In our implementation of the parametrisation method, this is cured by augmenting the size of the problems to solve with a bordering technique[50, 51, 52].

At the end of the process, the reduction method provides the nonlinear mapping (42) and the reduced dynamics (43) at arbitrary order. To study the convergence of the ROM, an expansion order needs to be selected. Since the two different parts of the master coordinate $\mathbf{z}$ in Eq. (44) have different meanings, the truncation of the ROM will be hereafter referred to as $\mathcal{O}(p, q)$. The first integer $p$ refers to the maximal order of the $\mathbf{z}$ coordinate, while $q$ refers to the maximal order of $\tilde{\mathbf{z}}$. As a matter of fact, $q$ directly represents the order $\varepsilon^q$ of the truncation of the non-autonomous part. In general, the order $q$ needed to reach convergence is much smaller than that needed for $p$, see the numerical results in Section 6 and those reported in previous works [52].

Lastly, note that, at the end of the procedure, a minimal ROM described by Eq. (43) is obtained. The size of this ROM coincides with the size of $\mathbf{z}$, which, for a forcing resonant with a single mode of the system, will be equal to 4. However, since the coordinates $\tilde{\mathbf{z}}$ are known from Eq. (45), they do not constitute free variables hence the effective size of the ROM is equal to 2. Moreover, different values of the forcing can be simulated from a single run of the DPIM procedure by simply varying the $\varepsilon$ parameter.
# 5 Simplified approach for large parallel plates

A model problem which provides an excellent benchmark for the proposed formulation is represented by a long clamped-clamped beam facing a parallel rigid electrode of the same length with initial gap $g$, small with respect to the length. In this case a very accurate analytical approximation of the electrostatic forces $\mathbf{b}$ is given by the *parallel plates* formula which assumes an electric vector field that is orthogonal to the mid-plane of the plates, and thus colinear with the $\mathbf{e}_2$ vector in our notations, yielding:

$$ \mathbf{b}(\mathbf{x}) = \frac{1}{2} \varepsilon_0 \frac{(V_{\text{DC}} + V_{\text{AC}} \sin \omega t)^2}{(g - u_2(\mathbf{x}))^2} \mathbf{e}_2, \qquad (47) $$

where $\mathbf{x}$ stands for the point position in the initial configuration and the force is exerted only on the upper surface $U_0$ facing the electrode. This approximation ignores fringing fields and components of the force tangent to the electrodes, but is known to be extremely accurate for small gaps between large electrodes [23, 4]. Fringing fields may in general affect the capacitance value, but have negligible effects on the force, which is associated with the capacitance variation. In real electromechanical MEMS, gaps between electrodes are generally much smaller than the other dimensions of the electrodes, so this simplification has many practical applications.

In this case, a simplified formulation of the MGF can be obtained by using the analytical expression for $\mathbf{b}$, yielding:

$$ \int_{\Omega_0} \rho_0 \ddot{\mathbf{u}} \cdot \tilde{\mathbf{u}} \, d\Omega_0 + \int_{\Omega_0} \sigma^{\text{PK}}[\mathbf{u}] : (\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}]) \, d\Omega_0 = \frac{\varepsilon_0}{2} \int_{U_0} \frac{(V_{\text{DC}} + V_{\text{AC}} \sin \omega t)^2}{(g - u_2(\mathbf{x}))^2} \tilde{u}_2 \, dS_0, \quad \forall \tilde{\mathbf{u}} \in \mathcal{C}_s(0). \qquad (48) $$

This model will be referred to as the parallel-plate formulation (PPF) in the rest of the paper. The analytical expression of $\mathbf{b}$ substantially simplifies the problem such that direct time integration routines can be directly applied to Eq. (48). The PPF will thus be used in the next developments to provide time-domain full-order reference solutions. However, Eq. (48) still contains non-polynomial nonlinearities, which represents a major obstacle for several solution schemes and for the parametrisation method for invariant manifold in particular. In analytical approaches, like the multiple scales technique used for instance in [21, 23], two solutions are adopted. The most direct approach consists in expanding $1/(g - u_2(\mathbf{x}))^2$ in Taylor series. However, the series converges very slowly and a very large number of terms is required when $u_2$ reaches significant fractions of the gap. As an alternative *exact* approach which is also strictly connected with the general formulation developed in the previous section, one can introduce an auxiliary field:

$$ \psi(\mathbf{x}) = \frac{V_{\text{DC}} + V_{\text{AC}} \sin \omega t}{g - u_2(\mathbf{x})}. \qquad (49) $$

Remark that $\psi$ has the meaning of a gradient of the potential, so that a notation similar to $\psi$ in the general formulation of the previous sections is adopted. The power of electric forces are quadratic in $\psi$ and Eq. (48) becomes:

$$ \int_{\Omega_0} \rho_0 \ddot{\mathbf{u}} \cdot \tilde{\mathbf{u}} \, d\Omega_0 + \int_{\Omega_0} \sigma^{\text{PK}}[\mathbf{u}] : (\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}]) \, d\Omega_0 = \frac{\varepsilon_0}{2} \int_{U_0} \psi^2(\mathbf{x}) \tilde{u}_2 \, dS_0, \quad \forall \tilde{\mathbf{u}} \in \mathcal{C}_s(0). \qquad (50) $$

Moreover, the strong form (49) is replaced with a weak form of the compatibility equation:

$$ \int_{U_0} \psi(\mathbf{x})(g - u_2(\mathbf{x}))\tilde{\psi}(\mathbf{x})dS_0 = (V_{\text{DC}} + V_{\text{AC}} \sin \omega t) \int_{U_0} \tilde{\psi}(\mathbf{x})dS_0, \quad \forall \tilde{\psi}(\mathbf{x}) \in \mathcal{C}_{\psi}, \qquad (51) $$

where $\tilde{\psi}$ and $\psi$ are chosen in the same space $\mathcal{C}_{\psi}$. In Section 6, the results presented are generated by an implementation where $\mathcal{C}_{\psi}$ is the space of piecewise constant shape functions.

Eqs. (50)-(51) represent the mixed semi-analytical formulation for parallel plates problems only, and is referred to as MPPF. It can be discretised following the same steps as for the MGF and then reduced using the DPIM procedure. Details are collected in Appendix B. The MPPF model will be used to validate numerical results in Section 6.

# 6 Numerical results

In this section, we collect some numerical tests intending to validate the proposed formulations.

First, a case of static pull-in instability is addressed in Section 6.1 using the mixed electrostatic formulation introduced in Section 2.1 and a graphical representation of the new associated fields is also provided for the sake of illustration. Next, in Section 6.2.1 the different ROM strategies are applied to the case of a Parallel Plate setting, where a clamped-clamped beam faces a parallel rigid electrode of the same length. Despite its simple geometrical layout, this configuration activates several nonlinear phenomena when computing either primary or secondary resonances, and hence represents a challenging testbed.
## 6.1 Pull-in verification

To validate the proposed mixed general formulation MGF for the electrostatic field, we examine the static pull-in phenomenon that occurs when an electrostatically actuated system reaches a critical voltage, causing its movable component to abruptly snap into contact with a fixed electrode due to electrostatic forces.

<figure><img src="image_3.png" /><figcaption>Figure 2: Static pull-in verification. Green mesh: deformable structure $\Omega_0$; orange mesh: air domain $A_0$. The lower deformable mass is anchored by two thin beams on either side to model linear spring constraints.</figcaption></figure>

A classical result concerns parallel plates in which a rigid plate, with given voltage $V_{DC}$ and constrained by linear springs, faces a parallel fixed electrode with zero voltage. When $V_{DC}$ increases quasi-statically, the former plate stably approaches the electrode until the displacement is equal to one-third of the gap. Above this value, the system reaches the pull-in instability, and a further increase of the voltage induces a dynamic evolution where the movable plate will snap to the electrode. To reproduce such a set-up numerically, the case shown in Fig. 2 is considered. The fixed upper electrode is shown in white, while the movable part, in green, consists of a central stiff mass, constrained by two thin beams on either side. The nonlinear mechanical effects in these beams are deliberately disabled to simulate linear springs. The problem is two-dimensional, and the mass has a length of 6 µm and a height of 3 µm, while the beam thickness is 0.5 µm for a length of 5 µm. The gap $g$ between the electrode and the mass is 1 µm. The fictitious elastic constants are $E = 154$ MPa and $\rho = 2330$ kg/m$^3$.

Figures 3 and 4 present the solution fields for the particular case of $V_{DC} = 10$ V. In particular, Figure 3 is the contour plot of the vertical displacement of the structure $u_{02}$, and its extension $s_{02}$ to the air domain $A_0$. The overall field is continuous and vanishes on $S_\infty$ and on the surface of the electrode.

<figure><img src="image_4.png" /><figcaption>Figure 3: Static pull-in verification. Static analysis with $V_{DC} = 10$ V. Vertical displacement $u_{02}$ in µm in the solid and its extension $s_{02}$ to the air domain. The two fields are continuous on the solid surface. The field is plotted in the deformed configuration for graphical purposes, but the formulation is fully Lagrangian.</figcaption></figure>

The plots of Figures 4a-4c correspond to the potential $\varphi_0$ and the two components of the Eulerian gradient $\psi_0$, Eq. (16), always for $V_{DC} = 10$ V. It can be appreciated that the electric field is almost constant and vertical between the plates, which justifies the comparison with parallel-plate formulas for the pull-in value.

Next, the voltage bias between the electrode and the top side of the mass is raised, and the static deformation at different voltage levels is computed. Figure 5 reports the evolution of the maximum displacement and the first eigenvalue $\omega_1$ as the potential varies. The plot clearly shows that the mass displacement increases with the actuating voltage in a stable manner; nevertheless, as the midspan displacement $u_m$ approaches the critical value of one-third of
<figure><img src="image_5.png" /><figcaption>Figure 4: Static pull-in verification. Static analysis with $V_{DC} = 10$ V. (a) potential $\varphi_0$. The values of $V_{DC}$ assigned on the electrode and the body, respectively 0 V and 10 V, correspond to $\varphi_0$ on the surface of the bodies. Remark the linear distribution between the plates; (b) Horizontal component $\psi_{01}$ of the Eulerian gradient showing singularities near the edges of the electrodes; (c) Horizontal component $\psi_{02}$ of the Eulerian gradient which is nearly constant within the parallel plates, corresponding to the presence of a uniform electric field.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

Figure 5: Static pull-in verification. Evolution of the relative displacement $u_m/g$ (left axis) of the lower rigid mass and of the first eigenfrequency (right axis) at different actuating voltages. As a consequence of the pull-in phenomenon, when the displacement approaches one-third of the initial gap, even a small voltage increase causes a large change in the displacement and the first eigenvalue vanishes. The proposed mixed MGF formulation is compared with the classical formulation utilising only the potential $\varphi$ showing negligible differences.

the initial gap, even small increments in the actuating voltage result in significant changes in deformation, following the typical behaviour observed in pull-in phenomena. Similarly, the first eigenfrequency gradually decreases due to the softening influence of the potential difference, then rapidly converges to 0 as the critical voltage threshold is approached, confirming that the instability phenomenon is correctly captured.

## 6.2 Reduced order modelling for parallel plates

We now turn to the analysis of the dynamical response of a clamped-clamped microbeam of length $L = 510$ µm, thickness $h = 1.5$ µm and air-gap $g = 1.18$ µm with respect to a parallel electrode that has the same length. The dimensions correspond to those of realistic microstructures, see for instance [21, 22, 23], and are chosen to trigger the wealth of nonlinear effects pointed out in what follows. The microbeam is excited by a voltage bias of the form $V_{DC} + V_{AC} \sin(\omega t)$. We also assume that the beam is subjected to an axial force of 0.0009 N. The elastic constants are $E = 154$ GPa and $\rho = 2330$ kg/m$^3$ while the Poisson coefficient is set to 0 to simulate a large thin plate.

In the following section, the various formulations proposed in this paper are used to generate the Frequency Response Curves (FRCs) for the mid-span displacement $u_m$ of the beam, i.e. the amplitude of the oscillation of the centre point when the frequency $\omega$ of the $V_{AC}$ input signal spans a specific range of frequencies.

In all the applications, a single master mode is selected, excluding a priori possible phenomena of internal resonance. The master mode corresponds to the first bending mode of the clamped-clamped beam interacting with the surrounding electric field. This implies that the size of the normal coordinates vector $z$ in the reduced order model $\dot{z} = f(z)$ in Eq.(43) is 4. The first two coordinates, $z_1$ and $z_2 = \bar{z}_1$, correspond to the normal coordinates representing the master mode, while the last two, ($z_3, z_4$), are the non-autonomous variables representing the harmonic forcing, such that $z_3 = e^{i\omega t}$ and $z_4 = e^{-i\omega t}$.

Both primary and secondary resonances will be investigated in the following sections. For convenience, we gather in Table 1 the values utilised in the different analyses for $V_{DC}$, $V_{AC}$ and for the actuating frequency $\omega$.

As a general remark for all the examples addressed, it should be recalled that in real electromechanical MEMS gaps between electrodes are generally very small to achieve strong actuating forces. Moreover, maximum displacements cannot exceed fractions of the gap, as they would induce other instability phenomena known as pull-in, which are carefully avoided in practice. This physical constraint would be met at lower amplitudes than the convergence limit of the reduction technique.
<table><caption>Table 1: Parameter values for the different test cases performed.</caption><thead><tr><th></th><th><i>V</i><sub>DC</sub> [V]</th><th><i>V</i><sub>AC</sub> [V]</th><th><i>ω</i> [Hz]</th></tr></thead><tbody><tr><td>Primary resonance - hardening</td><td>1.2</td><td>0.2</td><td>≈ <i>ω</i><sub>1</sub></td></tr><tr><td>Primary resonance - softening</td><td>4.0</td><td>0.045</td><td>≈ <i>ω</i><sub>1</sub></td></tr><tr><td>Superharmonic resonance</td><td>1.7</td><td>0.7</td><td>≈ 0.5 <i>ω</i><sub>1</sub></td></tr><tr><td>Parametric (subharmonic) resonance</td><td>1.7</td><td>0.7</td><td>≈ 2 <i>ω</i><sub>1</sub></td></tr></tbody></table>

### 6.2.1 Comparison between the reduced order models of the general MGF and parallel-plate MPPF formulations

This section aims to validate the parallel plate assumptions by comparing the results provided by the MGF with those of the MPPF, with the ultimate goal of using in the next section the more agile and cheaper MPPF to provide full-order reference solutions. As a consequence, we start by comparing the outcomes of the ROMs obtained from the two formulations. Since the DPIM procedure has already been shown to deliver very accurate results for purely mechanical vibrating systems, see e.g. [50, 52, 51], we will use this strategy to validate the fact that MPPF provides trustworthy solutions, in line with numerous previous studies showing the accuracy of the parallel plate analytical formula.

Two cases of actuation are selected to draw out this comparison. First, a small value of $V_{DC}$ is imposed to get a hardening behaviour for the vibration of the fundamental mode. Next, $V_{DC}$ is increased and the negative spring effect due to the electrostatic field becomes dominant, turning the overall beam response from hardening to softening, a well-known effect in electrostatic MEMS. For each formulation, MGF and MPPF, a ROM has been derived with the DPIM up to order $O(5,1)$ using the complex normal form style. The ROM is then solved by means of a continuation technique implemented in the software Matcont [62] to obtain the frequency-response curve (FRC) in the vicinity of the fundamental eigenfrequency.

Figures 6a and 6b collect the results for two different sets of parameters. The case of Figure 6a corresponds to $V_{DC} = 1.2$ V, $V_{AC} = 0.2$ V and $\alpha_{30} = 3 \times 10^{-3} s^{-1}$, $\alpha_{45} = 4.5 \times 10^{-3} s^{-1}$. Here the $V_{DC}$ bias is small and the beam retains the classical hardening behaviour of a purely mechanical clamped-clamped beam. In Figure 6b, on the contrary, we set $V_{DC} = 4$ V and $V_{AC} = 0.045$ V, a situation much closer to pull-in.

Note that, for the sake of simplicity, the stability analysis is not reported in the figures. Since the cases under study refer to classical cases in nonlinear vibrations, the stability of the different branches of solution can be easily inferred from the results, and are correctly predicted by the ROMs. The level of actuation has been set such that the maximum vibration amplitude reported in Fig. 6 corresponds to 1/3 of the thickness and more than 1/3 of the gap $g$, hence giving rise to important nonlinearities with respect to the setup.

In the numerical results, two effects should be noticed. First, there is a shift of the overall FRC towards smaller frequencies. Indeed, increasing the $V_{DC}$ value the negative spring effect due to the electrostatic field overcomes the expected mechanical hardening.

Secondly, the nonlinear response of the beam is now dominated by electrostatic forces and becomes softening. In general, the agreement between the two different approaches is very satisfactory, undoubtedly showing that for this application and the range of parameters used, the simplification of the electric field is perfectly well justified. Consequently, the MPPF formulation will be used in the rest of the paper to produce the FOM solutions.

### 6.2.2 Comparison between the PPF full-order model and the ROM

Since the reduced models obtained with the MGF and MPPF are fully equivalent, as shown in the previous section, we now compare the results obtained from the MPPF-ROM with those stemming from a time integration approach, directly applied to the full-order PPF model, Eq.(48), which will be used as FOM reference solution.

In particular, we analyse the same cases of hardening and softening behaviour as in the previous section. In Figure 6c, the voltages are set to $V_{DC} = 1.2$ V and $V_{AC} = 0.2$ V, while in Figure 6d, $V_{DC} = 4$ V and $V_{AC} = 0.045$ V. The structural losses are modelled with a mass-proportional Rayleigh damping $C = \alpha M$, and the Frequency Response Curves are computed for two different dissipation values, respectively $\alpha_{45} = 4.5 \times 10^{-3} s^{-1}$ and $\alpha_{30} = 3 \times 10^{-3} s^{-1}$.

The FRCs of the MPPF-ROM is obtained with the continuation package Matcont, while the FOM is integrated in time with a nonlinear implicit Newmark scheme, with $\beta = 1/4$, $\gamma = 1/2$. In the latter case, each cycle is partitioned in 120 steps and 3Q cycles are simulated for each value of $\omega$ always using a homogeneous state as an initial condition for the analysis, where Q is the quality factor. This time integration allows to follow only the stable branches, while the continuation technique can reproduce the whole Frequency Response Curve.

In all four cases, the ROM and FOM results are well aligned, providing additional confirmation of the very good accuracy provided by the reduction method.
<figure><img src="image_7.png" /><figcaption>Figure 6: Parallel plates. (a-b) Comparison between the responses of the Reduced Order Models MPPF-ROM and MGF-ROM. Frequency response curves for the mid-span displacement $u_m$. All the plots report the amplitude of the oscillation around the static solution. (a) Hardening case: $V_{DC} = 1.2$ V, $V_{AC} = 0.2$ V, with two different levels of dissipation $\alpha_{30} = 3 \times 10^{-3} s^{-1}$ and $\alpha_{45} = 4.5 \times 10^{-3} s^{-1}$. (b) Softening case $V_{DC} = 4.0$ V, $V_{AC} = 0.045$ V with the same dissipation levels. Increasing the $V_{DC}$ value, the negative spring effect due to the electrostatic field dominates and the overall beam response turns from hardening to softening. (c-d) Comparison of the full Full Order Model PPF, time-integrated, with the Reduced Model MPPF-ROM, solved with a numerical continuation procedure. The same actuation levels and dissipation values as in (a-b) are utilised.</figcaption></figure>

### 6.2.3 Convergence study

Having showcased the MPPF-ROM capability to yield accurate results, aligning closely with the solution provided by the full model, this section aims to assess the convergence rate of the ROM obtained thanks to the DPIM procedure towards the actual solution, while increasing the expansion order. As the analyses described in the previous section already made apparent, a rapid convergence has been observed when increasing the orders in the parametrisation method, in the case of primary resonance. Specifically, in the hardening test case, convergence is achieved as early as the $O(5,1)$ expansion order, as shown in Fig. 7a. However, when softening effects are more prominent and the pull-in instability is approached, the rate of convergence tends to slow down: Fig. 7b underscores that there remains a noticeable difference between order $O(5,1)$ and $O(7,1)$ expansions, indicating the necessity of high order expansions
<figure><img src="image_8.png" /><figcaption>Figure 7: Parallel plates. Convergence analysis of the ROM for the primary resonance. (a): hardening behaviour with low damping. $V_{\text{DC}} = 1.2$ V, $V_{\text{AC}} = 0.2$ V, $\alpha = 3 \times 10^{-3}$ s$^{-1}$; (b) softening behaviour with low damping. $V_{\text{DC}} = 4.0$ V, $V_{\text{AC}} = 0.045$ V, $\alpha = 3 \times 10^{-3}$ s$^{-1}$. The amplitudes computed with the PPF-FOM and obtained by direct time integration are also reported as black stars for comparison.</figcaption></figure>

near pull-in. Pushing the order further only yields minimal improvement, as the ROM with $O(9, 1)$ expansion almost overlaps with the $O(7, 1)$ one. Concerning the convergence of the non-autonomous development, which is set to 1 in the previous study, it has been numerically verified (and not shown here for the sake of brevity) that increasing the order does not change the FRC which is indeed converged with an order 1 on the forcing term. This result is in line with previous studies on primary resonances, which were limited to this order of accuracy on the forcing term and showed a good convergence rate, see for example the discussions and examples reported elsewhere [52, 51, 48]. The amplitudes predicted by the PPF-FOM, and computed with a direct numerical time integration, are also reported in Fig. 7 for the sake of completeness.

### 6.2.4 Secondary resonances

In this final Section, the nonlinear vibrating behaviour of the coupled electro-mechanical system when the forcing frequency $\omega$ is in the vicinity of the secondary resonances, corresponding to $\omega \simeq \omega_1/2$, and $\omega \simeq 2\omega_1$ (with $\omega_1$ the linear frequency of the fundamental bending mode), is investigated. Comparisons between MPPF-ROM and PPF-FOM will be given whenever possible, in order to show the accuracy of the reduction technique in a variety of adverse nonlinear cases.

The first case under study is that of the superharmonic resonance, when the frequency of the actuation voltage approaches $\omega_1/2$. As shown in [52], a first-order expansion on the non-autonomous term is not able to cope with such superharmonic resonance, and the ROM needs to include an $\epsilon$-order that is at least equal to the order of the superharmonic considered, 2 in this case. Fig. 8 reports the results found when the actuation is in the vicinity of half the eigenfrequency of the fundamental bending mode. Three ROMs with increasing orders are computed to achieve convergence: $O(3, 2)$, $O(5, 4)$ and $O(7, 6)$. The results are compared to a FOM simulation obtained by running a time integration directly on the PPF formulation using the same settings as in the previous cases (Newmark scheme, 120 points per period). Whereas the truncation order $O(3, 2)$ is not converged, the overlap of the two other curves underlines that an order $O(5, 4)$ is sufficient to reproduce this nonlinear secondary resonance. It is important to mention that orders of the truncation of the non-autonomous term lower than 2 are unable to capture this superharmonic behaviour. One can note that such behaviour on secondary resonances has already been reported [63, 22], with a simplified model for the beam and the electrostatic coupling. Here the model computes the structural vibrations with an FE model without kinematical assumptions, and fully takes into account the electromechanical coupling, offering quantitative predictions on the amplitude levels one can reach in such superharmonic regime that could be directly compared e.g. to experiments. The direct comparisons provided with the FOM solution undoubtedly illustrates the convergence of the ROM to the actual solution.
<figure><img src="image_9.png" /><figcaption>Figure 8: Parallel plates. Superharmonic resonance with $\omega \approx \omega_1/2$. $V_{DC} = 1.7$ V, $V_{AC} = 0.7$ V, and mass-proportional damping with $\alpha = 3 \times 10^{-3} s^{-1}$. Convergence of the MPPF-ROM with increasing orders in the DPIM, and comparison to direct numerical time integration of the PPF-FOM.</figcaption></figure>

The last case under study is the secondary resonance occurring when the forcing frequency is in the vicinity of two times the eigenfrequency of the fundamental mode. This case has already been analysed with von Kármán assumptions for the beam and first-order multiple scales solutions [63, 22], and was termed subharmonic resonance. However, the main driver of the instability is a parametric resonance, as shown in the solvability conditions derived in [63], that contains the same term as the classical parametric resonance [2, 64]. The resonance curves shown in [63, 22] are also typical of a parametric excitation, and not a subharmonic resonance of a single degree of freedom oscillator, which is characterised by a detached solution branch [2, 65]. As a matter of fact, this resonance scenario can be interpreted in light of the dynamical system equations (35), and appears to be a subharmonic resonance for the complete system, but a parametric resonance for the mechanical vibratory systems, once the other equations have been condensed into the dynamical problem. For this reason, we will refer to this case hereafter as a parametric excitation.

Before showing the numerical results, let us first highlight how the reduced dynamics provided by the ROM can also be used to analyse, understand, and predict the awaited behaviours. In this specific case of parametric resonance in a subharmonic excitation scenario, it provides a fruitful physical interpretation of the nonlinear dynamical phenomena at hand, in the line of the analysis reported for example in [66]. Truncating the ROM to order 3 for both the autonomous ($z_1, \bar{z}_1$) coordinates and the non-autonomous ones ($z_3, z_4$), the DPIM procedure using complex normal form [50, 52] provides the following reduced dynamics:

$$ \dot{z}_1 = \lambda_1 z_1 + f_1 z_1^2 \bar{z}_1 + f_2 \bar{z}_1 z_3 + f_3 z_1 z_3 z_4 + \mathcal{O}(|z|^3), \quad (52) $$

where $\lambda_1 = -\xi_1\omega_1 + i\omega_1\sqrt{1 - \xi_1^2}$ is the complex eigenvalue of the fundamental mode, and $\xi_1$ refers to the modal damping ratio. The second monomial corresponds to the trivially resonant cubic term, while the third one represents the parametric excitation term, since $z_3 = e^{i\omega t}$, such that $f_2$ is directly proportional to the amplitude of the forcing. Finally, the last monomial $z_1z_3z_4$ is a trivially resonant term of higher order in the non-autonomous term and needs an order 2 on the non-autonomous coordinates to be considered.

In Eq. (52), the driver of the parametric instability is the second-order term $f_2\bar{z}_1z_3$ which reads $f_2\bar{z}_1e^{i\omega t}$. This monomial is indeed resonant when $\omega$ is in the vicinity of $2\omega_1$. In particular, through a period-doubling bifurcation from the main branch of zero solution, a periodic response occurs in the form:

$$ z_1(t) = \rho e^{i\phi} e^{i\omega t/2}, \quad (53) $$

where a polar representation of the normal coordinate is introduced with an amplitude $\rho$. In order to show how the ROM predicts the bifurcation points, let us truncate Eq. (52) to its lowest-order significant terms by neglecting the
two monomials of order three with coefficients $f_1$ and $f_3$. Inserting the ansatz (53) into the normal form with only the $f_2$ term of the parametric instability, and assuming a steady-state response since we are interested in the solution branches of the permanent regime, one obtains:

$$ \rho(\dot{i}\omega/2 - \lambda_1 - f_2e^{-2i\phi}) = 0 + \mathcal{O}(\rho^2) \quad (54) $$

Hence, at the lowest order, a non-zero branch of solution with $\rho \neq 0$ exists if and only if the term into the bracket vanishes, which gives rise to the instability region of this problem. Separating real and imaginary parts of this term, and eliminating the angle $\phi$ thanks to the relationship $\cos^2 2\phi + \sin^2 2\phi = 1$, yields:

$$ \Omega_{a,b} = 2\omega_1\sqrt{1 - \xi_1^2} \pm 2\sqrt{f_2^2 - \xi_1^2\omega_1^2}. \quad (55) $$

This equation gives the two limit points $\Omega_{a,b}$ where the two branches of periodic solutions are connected to the solution at rest $\rho = 0$, i.e. the two period-doubling bifurcation points.

The analysis can be even pushed further by considering the term $f_3z_1z_3z_4$ in addition. Since $z_3 = e^{i\omega t}$ and $z_4 = e^{-i\omega t}$, this trivially resonant monomial simply reads $f_3z_1$ and will thus have a direct consequence on the linear term and on the oscillation frequency, which is perturbed. Considering an order 2 on the non-autonomous coordinate will thus change the prediction of the location of the bifurcation points given in Eq. (55), which needs to be rewritten by changing $\omega_1$ to $\omega_1 + f_3$. All those predictions will be assessed thanks to the numerical results, underlining the predictive capability and the wealth of the ROM.

<figure><img src="image_10.png" /><figcaption>Figure 9: Parallel plates. Parametric (subharmonic) resonance with $\omega \approx 2\omega_1$. $V_{DC} = 1.7$ V, $V_{AC} = 0.7$ V, and mass-proportional damping with $\alpha = 3 \cdot 10^{-3} s^{-1}$. (a) Convergence of the ROM with increasing orders. (b) Close-up view and comparison with direct numerical integration points from the PPF-FOM.</figcaption></figure>

The first numerical experiment considers mass-proportional damping as in the previous simulations. This choice of Rayleigh damping has been selected in the full-order PPF formulation Eq.(48) because the mass matrices in the PPF and MPPF formulations are identical. On the other hand, the stiffness terms are different, as in the current implementation of the FOM electrostatic terms are directly included in the stiffness matrix. Consequently, a stiffness proportional damping term would generate different FOM and ROM models, making the comparisons much more difficult to handle. Fig. 9(a) reports the convergence of the ROM obtained from the DPIM with increasing orders both on the autonomous and non-autonomous parts. As awaited from previous theoretical derivations using simpler models, the subharmonic resonance gives rise here to a parametric instability, and the typical shape of the bifurcated branches are retrieved, with a satisfactory convergence of the ROM obtained from order $\mathcal{O}(5, 2)$. In this case however, it has been found that the bifurcated branches do not close at high amplitudes. This behaviour is attributed to using mass-proportional damping in the simulations. Indeed, the modal damping coefficients of the slave modes are decreasing with increasing frequency. Consequently, the nonlinear damping coefficients, that gather all the losses contributions of the slave modes in the reduced dynamics (see e.g. [40, 41] for more details), are not important enough, leading to this behaviour in the bifurcated branches that do not close and stay open. This specific point can be addressed by using stiffness-proportional damping, as it will be shown next. Fig. 9(b) shows a close-up view in the
vicinity of the bifurcation points, and compares the outcomes of the ROM to direct time integration performed on the PPF-FOM. Again, a very good agreement is found in terms of the amplitudes of the bifurcated solutions, as well as the location of the bifurcation points. In particular, the analytical development shown before, highlighting that the period doubling point will shift as soon as an order 2 in the non-autonomous part is considered, is retrieved from the numerical simulation, underlining the accuracy of the ROM in predicting such details.

<figure><img src="image_11.png" /><figcaption>Figure 10: Parallel plates. Parametric resonance with Rayleigh damping. $V_{\text{DC}} = 1.7$ V, $V_{\text{AC}} = 0.7$ V, $\alpha = -3.81 \times 10^{-2} s^{-1}$, $\beta = 2 \times 10^{-1}$ s. (a) Convergence of the autonomous part of the expansion. (b) Convergence of the non-autonomous part of the expansion</figcaption></figure>

The numerical results reported in Figs. 10a-10b conclude this section by considering the same case of the parametric resonance, where a complete Rayleigh damping including a stiffness-proportional term is now included in the ROM simulation. The damping matrix used to build the ROM now reads $\alpha M + \beta K$, and numerical values have been set as $\alpha = -3.81 \times 10^{-2} s^{-1}$ and $\beta = 2 \times 10^{-1} s$. This case is here analysed to show that taking into account a more realistic loss representation which increases with frequencies, allows to retrieve closed bifurcating branches in this parametric resonance. This is the direct consequence of the fact that in the ROM derivation, the damping of all the slave modes is aggregated into nonlinear damping coefficients in the reduced dynamics, offering a refined representation of the losses when orbits are decaying along the invariant manifold [40, 41]. Note that for this particular case of stiffness-proportional damping, a FOM simulation using PPF formulation has not been set up for the reasons explained before. Since the focus of this work is on the derivation of ROMs using the DPIM, and not on the set-up of accurate FOM, only the predictions of the ROM are here analysed.

Fig. 10a shows the bifurcated solution branches of the parametric resonance for increasing odd orders (3, 5, 7 and 9) of the autonomous expansion, while keeping the non-autonomous truncation to 1. Interestingly, the y-axis represents the amplitude of the physical displacement of the centre point of the beam $u_m$. Note that this value is small but not vanishing. Indeed, since $\omega \approx 2\omega_1$, there is a residual non-resonant linear excitation, far from resonance, which makes the displacement solution non-zero here. This remark also holds for the mass-proportional damping case. Whereas the ROM coordinates $z$ exactly bifurcates from $\rho = 0$, it is the nonlinear mapping that relates physical coordinates to normal coordinates, which allows retrieving this non-zero (though small) value of the displacement. The bifurcated branches of solutions then show the classical portrayal of a parametrically excited system, with a hardening behaviour obtained with the parameters selected here. Again, stiffness-proportional damping is needed to reach a maximal amplitude in the bifurcated branches, which contrasts with the previous results shown in Fig. 9, and definitely highlights the importance of the damping of the slave modes and the fact that they play a crucial role for the nonlinear reduced dynamics.

Finally, Fig. 10b shows the convergence of the ROMs when increasing the order of the non-autonomous terms. As awaited from the analysis, an important shift of the solution is obtained when increasing the order of the non-autonomous coordinates to 2, since including the term $f_3 z_1 z_3 z_4$ in the analysis. As announced, taking this term into account leads to a slight frequency shift such that the location of the bifurcation points is slightly modified. Higher order terms then show minimal effects and the convergence is finally reached with order $O(9, 3)$.
## 6.3 Computational costs

All the computations presented in the preceding sections were performed using an Intel i9-12900K processor with 32 GB of RAM. The mesh employed for the parallel plate formulations, both for the reduced model and the full one, counts 5834 degrees of freedom. With this mesh, the PPF-FOM model with the Newmark scheme requires, on average, 67 minutes per computed point. The computational time needed to construct the ROM models depends on the selected order of expansion. We report in Table 2 the computational times for the converged ROM and for the FOM for each of the cases shown in the previous section.

<table><caption>Table 2: Computational times for the converged ROM and for the FOM in each of the presented cases. The total time of the FOM relates to the number of points reported in the figures. The FOM for the parametric with mass-proportional damping (open branches) shows a much slower convergence, taking upwards of 2000 periods to reach a steady state solution.</caption><thead><tr><th>Resonance type</th><th>Order</th><th>ROM</th><th>FOM (1 point, avg)</th><th>FOM (total)</th></tr></thead><tbody><tr><td>Primary - hardening</td><td>O(7, 1)</td><td>12s</td><td>67m</td><td>1008m</td></tr><tr><td>Primary - softening</td><td>O(9, 1)</td><td>35s</td><td>67m</td><td>472m</td></tr><tr><td>Superharmonic</td><td>O(7, 6)</td><td>45s</td><td>67m</td><td>801m</td></tr><tr><td>Parametric - open</td><td>O(7, 3)</td><td>24s</td><td>22h</td><td>133h</td></tr><tr><td>Parametric - closed</td><td>O(9, 4)</td><td>2m30s</td><td>-</td><td>-</td></tr></tbody></table>

These times include the computation of the static solution, employing a Newton-Raphson scheme, the solution to the generalised eigenvalue problem, and the solution to the homological equations for all of the monomials involved in the specific parametrisation. In particular, for this last step, a set of linear systems needs to be solved. Let $o$ be the maximum order and $n$ be the number of normal coordinates considered. The number of monomials at order $p$, $m_p$, and the total number of monomials, $M_o$, for which a linear system has to be solved are

$$m_p = \frac{1}{2} \frac{(p + n - 1)!}{p!(n - 1)!}, \qquad (56)$$

$$M_o = \frac{1}{2} \frac{(o + n)!}{o!n!}, \qquad (57)$$

where the factor one-half is a consequence of the fact that only one linear system needs to be solved for each pair of conjugate monomials. Note that the expression of $M_o$ given above only holds when the order is equal for the autonomous and non-autonomous parts of the expansion: if the latter is kept smaller, as it is done in some of the numerical results in the paper, the formula will change and the total number will decrease. For each of the $M_o$ monomials, a linear system of size $N$ (DOFs) is solved, so that the computational time of the overall computation scales with the number of DOFs as $O(N^3)$, exactly like a linear solve. To solve these systems, we exploited the Pardiso solver. It is important to note that, while a static solution and eigenvalue problem are also included, the bulk of the computational time is taken by the assemblage and solution of the $M_o$ linear systems. Moreover, since the entire assemblage procedure has been implemented in Julia, the native just-in-time compilation of the software is exploited to achieve a substantial speed-up over fully-interpreted languages.

Once the reduced model has been retrieved, the continuation procedure exploiting MatCont [62] for each curve takes from 10 seconds to 1 minute: the specific time within this range depends on the number of points computed on the curve.

## 7 Conclusions

In this paper, a direct implementation of the parametrisation method for invariant manifolds (DPIM) has been applied to a fully coupled multiphysics problem involving the nonlinear vibrations of deformable structures under the influence of an electrostatic field, hence significantly extending the realm of applications that can be dealt with this reduction technique.

The main difficulty was related to the expression of the nonlinear forces of the electromechanical problem, which are nonpolynomial, thus asking for a special treatment in order to be then automatically treated in the arbitrary order expansions issued from the DPIM procedure. To overcome this challenge, a new mixed fully Lagrangian formulation has been proposed, that exclusively incorporates explicit polynomial nonlinearities. This was achieved by rewriting
the coupled problem in its original configuration, and by introducing an auxiliary field strictly related to the electric field.

The governing equations have been semi-discretised using the standard finite element procedures, resulting in a system of differential-algebraic equations compatible with the general derivation of the DPIM for non-autonomous problems discussed in [52].

The proposed formulation is fully general and applicable to arbitrary geometries. For the purpose of validation, we have restrained ourselves to a 2D implementation of the standard, but challenging benchmark of a clamped-clamped beam facing a parallel electrode, where well-known and accurate simplified expressions for the electrostatic forces can be used for comparison.

This simplification enabled validation of the reduced-order models against full-order simulations obtained through direct time integration. Numerical results have shown a remarkable performance both in terms of accuracy and wealth of nonlinear effects that can be accounted for, including primary resonances, with the transition of the frequency response curves from hardening to softening, as well as secondary resonances. The convergence of the ROM has also been analyzed for a classical benchmark, and its very good predictive capacity is used to give insights into the nonlinear dynamical phenomena observed. Application of this nonlinear reduction technique paves the way for the generation of efficient and reliable ROMS with remarkable accuracy and efficient predictive capacities for electrostatically actuated resonating Micro Electro Mechanical Systems and provides a much-needed tool for the design and optimisation of a whole new family of devices with unprecedented performances. This research contributes significantly to advancing the understanding and application of the parametrisation method in tackling complex multiphysics problems, particularly those arising in the domain of MEMS.

## Acknowledgments

Alessio Colombo acknowledges the support of ST Microelectronics under the framework of the Joint Research Platform with the Politecnico di Milano, STEAM. Attilio Frangi acknowledges the PRIN 2022 Project "DIMIN- DIgital twins of nonlinear Microstructures with iNnovative model-order-reduction strategies" (No. 2022XATLT2) funded by the European Union - NextGenerationEU.

## References

[1] A. Corigliano, R. Ardito, C. Comi, A. Frangi, A. Ghisi, and S. Mariani. Mechanics of microsystems. John Wiley & Sons, 2018.

[2] A. H. Nayfeh and D. T. Mook. Nonlinear oscillations. John Wiley & sons, New-York, 1979.

[3] V. Zega, G. Gattere, S. Koppaka, A. Alter, G. D. Vukasin, A. Frangi, and T. W. Kenny. Numerical modelling of non-linearities in mems resonators. Journal of Microelectromechanical Systems, 29(6):1443-1454, 2020.

[4] A.Z. Hajjaj, N. Jaber, S. Ilyas, F.K. Alfosail, and M.I. Younis. Linear and nonlinear dynamics of micro and nano-resonators: Review of recent advances. International Journal of Non-Linear Mechanics, 119(103328):1-22, 2020.

[5] A. Frangi, A. Guerrieri, N. Boni, R. Carminati, M. Soldo, and G. Mendicino. Mode coupling and parametric resonance in electrostatically actuated micromirrors. IEEE Transactions on Industrial Electronics, 65(7):5962-5969, 2018.

[6] S. Baguet, V.N. Nguyen, C. Grenat, C.H. Lamarque, and R. Dufour. Nonlinear dynamics of micromechanical resonator arrays for mass sensing. Nonlinear Dynamics, 95(2):1203-1220, 2019.

[7] S.H. Nitzan, P. Taheri-Tehrani, M. Defoort, S. Sonmezoglu, and D.A. Horsley. Countering the effects of nonlinearity in rate-integrating gyroscopes. IEEE Sensors Journal, 16(10):3556-3563, 2016.

[8] C. Grenat, S. Baguet, C.-H. Lamarque, and R. Dufour. Mass sensing by symmetry breaking and localization of motion in an array of electrostatically coupled nonlinear mems resonators. International Journal of Non-Linear Mechanics, 140:103903, 2022.

[9] D.D. Shin, Y. Chen, I.B. Flader, and T.W. Kenny. Temperature compensation of resonant accelerometer via nonlinear operation. IEEE Micro Electro Mechanical Systems (MEMS), pages 1012-1015, 2018.

[10] Ben Rivlin and David Elata. Design of nonlinear springs for attaining a linear response in gap-closing electrostatic actuators. International Journal of Solids and Structures, 49(26):3816-3822, 2012.
[11] Jérôme Juillard, Pierre Prache, and Núria Barniol. Analysis of mutually injection-locked oscillators for differential resonant sensing. *IEEE Transactions on Circuits and Systems I: Regular Papers*, 63(7):1055–1066, 2016.

[12] M.H. Ghayesh and H. Farokhi. Bistable nonlinear response of mems resonators. *Nonlinear Dynamics*, 90(3):1627–1645, 2017.

[13] S. Houri, D. Hatanaka, M. Asano, R. Ohta, and H. Yamaguchi. Limit cycles and bifurcations in a nonlinear mems resonator with a 1:3 internal resonance. *Applied Physics Letters*, 114(10):103103, 2019.

[14] David A. Czaplewski, Changyao Chen, Daniel Lopez, Oriel Shoshani, Axel M. Eriksson, Scott Strachan, and Steven W. Shaw. Bifurcation generated mechanical frequency comb. *Phys. Rev. Lett.*, 121:244302, Dec 2018.

[15] David A. Czaplewski, Scott Strachan, Oriel Shoshani, Steven W. Shaw, and Daniel López. Bifurcation diagram and dynamic response of a MEMS resonator with a 1:3 internal resonance. *Applied Physics Letters*, 114(25):254104, 06 2019.

[16] O. Thomas, F. Mathieu, W. Mansfield, C. Huang, S. Trolier-McKinstry, and L. Nicu. Efficient parametric amplification in MEMS with integrated piezoelectric actuation and sensing capabilities. *Applied Physics Letters*, 102(16):163504, 2013.

[17] P.M. Polunin and S.W. Shaw. Mode coupling and parametric resonance in electrostatically actuated micromirrors. *International Journal of Non-Linear Mechanics*, 94:300–308, 2017.

[18] Shai Shmulevich, Inbar Hotzen Grinberg, and David Elata. A mems implementation of a classic parametric resonator. *Journal of Microelectromechanical Systems*, 24(5):1285–1292, 2015.

[19] A. Ganesan, C. Do, and A. Seshia. Phononic frequency comb via intrinsic three-wave mixing. *Physical Review Letters*, 118(033903), 2017.

[20] A. Ganesan, C. Do, and A. Seshia. Phononic frequency comb via three-mode parametric resonance. *Applied Physics Letters*, 112(021906), 2018.

[21] M. I. Younis and A. H. Nayfeh. A study of the nonlinear response of a resonant microbeam to an electric actuation. *Nonlinear Dyn.*, 31:91–117, 2003.

[22] A. H. Nayfeh and M. I. Younis. Dynamics of mems resonators under superharmonic and subharmonic excitations. *J. Micromech. Microeng.*, 15:1840–7, 2005.

[23] M. I. Younis. *MEMS linear and nonlinear statics and dynamics*. Springer, New-York, 2011.

[24] J. Juillard. Analysis of resonant pull-in of micro-electromechanical oscillators. *Journal of Sound and Vibration*, 350:123–139, 2015.

[25] G. Li and N.R. Aluru. A lagrangian approach for electrostatic analysis of deformable conductors. *J. Microelectromech. Syst.*, 11:245–254, 2002.

[26] S.K. De and N.R. Aluru. Full-lagrangian schemes for dynamic analysis of electrostatic MEMS. *J. Microelectromech. Syst.*, 13:737–758, 2004.

[27] R.C. Batra, M. Porfiri, and D. Spinello. Electromechanical model of electrically actuated narrow microbeams. *J. Microelectromech. Syst.*, 15:1175–1189, 2006.

[28] R.C. Batra, M. Porfiri, and D. Spinello. Review of modeling electrostatically actuated microelectromechanical systems. *Smart Mater. Struct.*, 16:R23–R31, 2007.

[29] J. J. Hollkamp and R. W. Gordon. Reduced-order models for non-linear response prediction: Implicit condensation and expansion. *Journal of Sound and Vibration*, 318:1139–1153, 2008.

[30] A. Frangi and G. Gobat. Reduced order modelling of the non-linear stiffness in MEMS resonators. *International Journal of Non-Linear Mechanics*, 116:211 – 218, 2019.

[31] G. Gobat, V. Zega, P. Fedeli, L. Guerinoni, C. Touzé, and A. Frangi. Reduced order modelling and experimental validation of a mems gyroscope test-structure exhibiting 1: 2 internal resonance. *Scientific Reports*, 11(1):16390, 2021.
[32] A. Givois, J.-F. Deü, and O. Thomas. Dynamics of piezoelectric structures with geometric nonlinearities: a non-intrusive reduced order modelling strategy. *Computers & Structures*, 253:106575, 2021.

[33] Y. Shen, N. Béreux, A. Frangi, and C. Touzé. Reduced order models for geometrically nonlinear structures: Assessment of implicit condensation in comparison with invariant manifold approach. *European Journal of Mechanics - A/Solids*, 86:104165, 2021.

[34] G. Haller and S. Ponsioen. Exact model reduction by a slow-fast decomposition of nonlinear mechanical systems. *Nonlinear Dynamics*, 90:617-647, 2017.

[35] E. Nicolaidou, T. L. Hill, and S. A. Neild. Indirect reduced-order modelling: using nonlinear manifolds to conserve kinetic energy. *Proc. R. Soc. A.*, 476:20200589, 2021.

[36] C. Touzé, A. Vizzaccaro, and O. Thomas. Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques. *Nonlinear Dynamics*, 105:1141-1190, 2021.

[37] S. W. Shaw and C. Pierre. Non-linear normal modes and invariant manifolds. *Journal of Sound and Vibration*, 150(1):170-173, 1991.

[38] S. W. Shaw and C. Pierre. Normal modes for non-linear vibratory systems. *Journal of Sound and Vibration*, 164(1):85-124, 1993.

[39] C. Touzé, O. Thomas, and A. Chaigne. Hardening/softening behaviour in non-linear oscillations of structural systems using non-linear normal modes. *Journal of Sound and Vibration*, 273(1-2):77-101, 2004.

[40] C. Touzé and M. Amabili. Non-linear normal modes for damped geometrically non-linear systems: application to reduced-order modeling of harmonically forced structures. *Journal of Sound and Vibration*, 298(4-5):958-981, 2006.

[41] A. Vizzaccaro, Y. Shen, L. Salles, J. Blahos, and C. Touzé. Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures. *Computer Methods in Applied Mechanics and Engineering*, 284:113957, 2021.

[42] A. Opreni, A. Vizzaccaro, A. Frangi, and C. Touzé. Model order reduction based on direct normal form: application to large finite element MEMS structures featuring internal resonance. *Nonlinear Dynamics*, 105:1237-1272, 2021.

[43] A. Opreni, M. Furlan, A. Bursuc, N. Boni, G. Mendicino, R. Carminati, and A. Frangi. One-to-one internal resonance in a symmetric mems micromirror. *Appl. Phys. Lett.*, 121(173501):1-6, 2022.

[44] X. Cabré, E. Fontich, and R. de la Llave. The parameterization method for invariant manifolds. III. Overview and applications. *J. Differential Equations*, 218(2):444-515, 2005.

[45] A. Haro, M. Canadell, J.-L. Figueras, A. Luque, and J.-M. Mondelo. *The parameterization method for invariant manifolds. From rigorous results to effective computations*. Springer, Switzerland, 2016.

[46] G. Haller and S. Ponsioen. Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction. *Nonlinear Dynamics*, 86(3):1493-1534, 2016.

[47] S. Jain and G. Haller. How to compute invariant manifolds and their reduced dynamics in high-dimensional finite-element models. *Nonlinear Dynamics*, 107:1417-1450, 2022.

[48] M. Li, S. Jain, and G. Haller. Nonlinear analysis of forced mechanical systems with internal resonance using spectral submanifolds – part I: Periodic response and forced response curve. *Nonlinear Dynamics*, 110:1005-1043, 2022.

[49] S. Jain, T. Thurnher, M. Li, and G. Haller. Computation of invariant manifolds and their reduced dynamics in high-dimensional mechanics problems. https://github.com/haller-group/SSMTool-2.1, 2021.

[50] A. Vizzaccaro, A. Opreni, L. Salles, A. Frangi, and C. Touzé. High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to large amplitude vibrations and uncovering of a folding point. *Nonlinear Dynamics*, 110:525-571, 2022.

[51] A. Opreni, A. Vizzaccaro, C. Touzé, and A. Frangi. High order direct parametrisation of invariant manifolds for model order reduction of finite element structures: application to generic forcing terms and parametrically excited systems. *Nonlinear Dynamics*, 111:5401-5447, 2023.
[52] A. Vizzaccaro, G. Gobat, A. Frangi, and C. Touzé. Direct parametrisation of invariant manifolds for forced non-autonomous systems including superharmonic resonances. submitted to Nonlinear Dynamics, 2023. arXiv:2306.09860.

[53] A. Opreni, A. Vizzaccaro, A. Martin, A. Frangi, and C. Touzé. MORFEInvariantManifold.jl: Nonlinear model order reduction of large dimensional finite element models using the direct parametrisation method for invariant manifolds. https://github.com/MORFEproject, 2022.

[54] A. Martin, A. Opreni, A. Vizzaccaro, M. Debeurre, L. Salles, A. Frangi, O. Thomas, and C. Touzé. Reduced order modeling of geometrically nonlinear rotating structures using the direct parametrisation of invariant manifolds. Journal of Theoretical, Computational and Applied Mechanics, 10430, 2023.

[55] A. Opreni, G. Gobat, C. Touzé, and A. Frangi. Nonlinear model order reduction of resonant piezoelectric micro-actuators: an invariant manifold approach. Computers and Structures, 289-107154:1–17, 2023.

[56] Mingwu Li, Shobhit Jain, and George Haller. Model reduction for constrained mechanical systems via spectral submanifolds. Nonlinear Dynamics, 111(10):8881–8911, May 2023.

[57] Aurélien Grolet, Alessandra Vizzaccaro, Marielle Debeurre, and Olivier Thomas. High order invariant manifold model reduction for systems with non-polynomial non-linearities: geometrically exact finite-element structures and validity limit. PREPRINT (Version 1), April 2024.

[58] A. Quarteroni and A. Valli. Numerical Approximation of Partial Differential Equations. Springer, New-York, 1994.

[59] S. Ponsioen, T. Pedergnana, and G. Haller. Automated computation of autonomous spectral submanifolds for nonlinear modal analysis. Journal of Sound and Vibration, 420:269 – 295, 2018.

[60] A. Haro and R. de la Llave. A parameterization method for the computation of invariant tori and their whiskers in quasi-periodic maps: numerical algorithms. Discrete & Continuous Dynamical Systems-B, 6(6):1261, 2006.

[61] T. Breunung and G. Haller. Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanical systems. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2213):20180083, 2018.

[62] A. Dhooge, W. Govaerts, and Y. A. Kuznetsov. Matcont: a matlab package for numerical bifurcation analysis of ODEs. ACM SIGSAM Bulletin, 38(1):21–22, 2004.

[63] E. M. Abdel-Rahman and A. H. Nayfeh. Dsecondary resonances of electrically actuated resonant microsensors. J. Micromech. Microeng., 13:491–501, 2003.

[64] J. J. Thomsen. Vibrations and stability. Advanced theory, analysis and tools. Springer, Berlin, Heidelberg, 2003. 2nd edition.

[65] M. Volvert and G. Kerschen. Phase resonance nonlinear modes of mechanical systems. Journal of Sound and Vibration, 511:116355, 2021.

[66] A. de Figueiredo Stabile, C. Touzé, and A. Vizzaccaro. Normal form analysis of nonlinear oscillator equations with automated arbitrary order expansions. Journal of Theoretical, Computational and Applied Mechanics, submitted, 2024.

# Appendix

## A Virtual power of internal stresses

In this section, we briefly summarise the steps left implicit in Section 8 that allows reformulating the terms associated to the power of internal stresses in Eq. (6). The Cauchy stress tensor $\sigma$ is replaced by the second Piola-Kirchhoff stress $\sigma^{PK}$ defined as

$$ \sigma^{PK} = Jf^{-1} \cdot \sigma \cdot f^{-T} $$

where $f[\mathbf{u}] = \mathbf{1} + \nabla \mathbf{u}$ is the transformation gradient and $J$ is the Jacobian of the transformation, so that

$$ \int_{\Omega} \sigma : \mathrm{grad} \tilde{\mathbf{u}} \, d\Omega = \int_{\Omega_0} (f[\mathbf{u}] \cdot \sigma^{\mathrm{PK}}[\mathbf{u}]) : \nabla^T \tilde{\mathbf{u}} \, d\Omega_0 $$
Applying standard rules of tensor calculus:

$$
\begin{align}
\int_{\Omega_0} (\mathbf{f}[\mathbf{u}] \cdot \boldsymbol{\sigma}^{\mathrm{PK}}[\mathbf{u}]) : \nabla^\mathrm{T} \tilde{\mathbf{u}} \, \mathrm{d}\Omega_0 &= \int_{\Omega_0} \boldsymbol{\sigma}^{\mathrm{PK}}[\mathbf{u}] : \mathrm{sym}(\nabla^\mathrm{T} \tilde{\mathbf{u}} \cdot \mathbf{f}) \, \mathrm{d}\Omega_0 \nonumber \\
&= \int_{\Omega_0} \boldsymbol{\sigma}^{\mathrm{PK}}[\mathbf{u}] : (\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}]) \, \mathrm{d}\Omega_0 \tag{58}
\end{align}
$$

where:

$$
\mathrm{sym}(\nabla^T \tilde{\mathbf{u}} \cdot \mathbf{f}) = \frac{1}{2} (\nabla \tilde{\mathbf{u}} + \nabla^T \tilde{\mathbf{u}} + \nabla^T \tilde{\mathbf{u}} \cdot \nabla \mathbf{u} + \nabla^T \mathbf{u} \cdot \nabla \tilde{\mathbf{u}}) \quad (59)
$$

and we have defined the operators

$$
\varepsilon[\mathbf{a}] = \frac{1}{2}(\nabla\mathbf{a} + \nabla^T\mathbf{a}) \qquad (60)
$$

$$
\bar{\nabla}[\mathbf{a}, \mathbf{b}] = \frac{1}{2}(\nabla^T \mathbf{a} \cdot \nabla \mathbf{b} + \nabla^T \mathbf{b} \cdot \nabla \mathbf{a}). \quad (61)
$$

Note also that:

$$
\mathbf{e}[\mathbf{u}] = \varepsilon[\mathbf{u}] + \frac{1}{2} \bar{\nabla}[\mathbf{u}, \mathbf{u}]
$$

so that:

$$
\begin{equation}
\begin{split}
\boldsymbol{\sigma}^{\mathrm{PK}}[\mathbf{u}] : (\boldsymbol{\varepsilon}[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}]) ={}& \boldsymbol{\varepsilon}[\tilde{\mathbf{u}}] : \mathcal{A} : \boldsymbol{\varepsilon}[\mathbf{u}] + \\
& \frac{1}{2} \boldsymbol{\varepsilon}[\tilde{\mathbf{u}}] : \mathcal{A} : \bar{\nabla}[\mathbf{u}, \mathbf{u}] + \bar{\nabla}[\tilde{\mathbf{u}}, \mathbf{u}] : \mathcal{A} : \boldsymbol{\varepsilon}[\mathbf{u}] + \\
& \frac{1}{2} \bar{\nabla}[\tilde{\mathbf{u}}, \mathbf{u}] : \mathcal{A} : \bar{\nabla}[\mathbf{u}, \mathbf{u}]
\end{split}
\tag{62}
\end{equation}
$$

The first term is linear, the second and the third are quadratic and the last term is cubic.

In many applications, the solution is sought as **u**₀ + **u** where **u**₀ is known and **u** is the unknown increment. Indeed, this is a slight abuse of notation, since now **u** denotes the increment and not the total displacement, but simplicity will benefit from this choice. It is then important to identify the new contributions in the integrand of Eq.(58):

$$
\sigma^{\mathrm{PK}}[\mathbf{u}] : (\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}])
$$

The factors can be expanded as:

$$
\mathbf{e}[\mathbf{u}_0 + \mathbf{u}] = \mathbf{e}[\mathbf{u}_0] + \varepsilon[\mathbf{u}] + \bar{\nabla}[\mathbf{u}_0, \mathbf{u}] + \frac{1}{2}\bar{\nabla}[\mathbf{u}, \mathbf{u}] \quad (63)
$$

$$
\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}_0 + \mathbf{u}, \tilde{\mathbf{u}}] = \varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}_0, \tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}] \quad (64)
$$

so that the product yields a constant term:

$$
\boldsymbol{\sigma}^{\mathrm{PK}}[\mathbf{u}_0] : (\boldsymbol{\varepsilon}[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}_0, \tilde{\mathbf{u}}])
$$

a linear term:

$$
\boldsymbol{\sigma}^{\mathrm{PK}}[\mathbf{u}_0] : \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}] + (\boldsymbol{\varepsilon}[\mathbf{u}] + \bar{\nabla}[\mathbf{u}_0, \mathbf{u}]) : \mathcal{A} : (\boldsymbol{\varepsilon}[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}_0, \tilde{\mathbf{u}}])
$$

a quadratic term:

$$
(\varepsilon[\mathbf{u}] + \bar{\nabla}[\mathbf{u}_0, \mathbf{u}]) : \mathcal{A} : \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}] + \frac{1}{2} \bar{\nabla}[\mathbf{u}, \mathbf{u}] : \mathcal{A} : (\varepsilon[\tilde{\mathbf{u}}] + \bar{\nabla}[\mathbf{u}_0, \tilde{\mathbf{u}}])
$$

and finally a cubic one:

$$
\frac{1}{2} \bar{\nabla}[\mathbf{u}, \mathbf{u}] : \mathcal{A} : \bar{\nabla}[\mathbf{u}, \tilde{\mathbf{u}}]
$$

The FEM discretisation of the internal power leads to:

$$
\tilde{\mathbf{U}}^{\mathrm{T}} (\mathbf{F} + \mathbf{KU} + \mathbf{G}(\mathbf{U}, \mathbf{U}) + \mathbf{H}(\mathbf{U}, \mathbf{U}, \mathbf{U})) \quad (65)
$$

where the matrices can be computed directly from the decomposition above.
## B Discretisation of the MPPF formulation

Let us consider the simplified approach for parallel plates of Section 5 formulated by the coupled continuous equations (50) and (51). Following the same procedure as for the MGF formulation, let us denote with $\mathbf{u}_0, \psi_0$ the static solution of the MPPF eqs.(50)-(51) when only $V_{DC}$ is acting and with $\mathbf{u} + \mathbf{u}, \psi_0 + \psi$ the total solution under the combined action of $V_{DC}$ and $V_{AC}$. Ignoring the constant terms which enter only in the computation of the static solution the virtual power of electrostatic forces in Eq.(50) becomes:

$$ P_e = \varepsilon_0 \underbrace{\int_{U_0} \psi_0 \psi \tilde{u}_2 \, dS_0}_{\text{linear}} + \frac{\varepsilon_0}{2} \underbrace{\int_{U_0} \psi^2 \tilde{u}_2 \, dS_0}_{\text{quadratic}} \qquad (66) $$

so that the semi-discretised equation for the mechanical equilibrium reads, also using Eq.(65):

$$ M\ddot{U} + KU + G(U, U) + H(U, U, U) = R_{\psi}\Psi + R_{\psi\psi}(\Psi, \Psi) \qquad (67) $$

The compatibility equation (51) has a right-hand side which generates the forcing while the left hand generates the terms:

$$ \underbrace{\int_{U_0} (\psi(g - u_{02}) - \psi_0 u_2) \tilde{\psi} \, dS_0}_{\text{linear}} + \underbrace{\int_{U_0} -\psi u_2 \tilde{\psi} \, dS_0}_{\text{quadratic}} \qquad (68) $$

The final form of the discretised compatibility equation is:

$$ C_u \mathbf{U} + C_{\psi} \mathbf{\Psi} + C_{\psi u}(\mathbf{\Psi}, \mathbf{U}) = \mathbf{F}_{\psi} V_{AC} \sin(\omega t) \qquad (69) $$

where the right-hand side vector comes from:

$$ \tilde{\mathbf{\Psi}}^{\mathrm{T}} \mathbf{F}_{\psi} = \int_{U_0} \tilde{\psi} \, dS $$

The system of Eqs. (67)-(69) is a Differential Algebraic System which can be reduced via the DPIM procedure as described in Section 4 for the equivalent DAE stemming from the MGF of Section 2.3.

## C Elements of differential geometry of a surface

Let us consider a surface $S_0$ embedded in a three-dimensional Euclidean space $x_1, x_2, x_3$ and mapped onto the parameter space $a_1, a_2$ (e.g. the parameters of a FEM discretisation).

We now introduce $\mathbf{g}_\alpha$, the two covariant base vectors:

$$ \mathbf{g}_{\alpha} = \mathbf{x}_{,\alpha} \left(= \frac{\partial \mathbf{x}}{\partial a_{\alpha}}\right) \qquad (70) $$

which are in general non-orthogonal, non-unitary vectors tangent to the surface. The vector product of the two covariant vectors yields the normal:

$$ \mathbf{g}_1 \wedge \mathbf{g}_2 = \mathbf{n} \cdot J_S \quad \text{where} \quad J_S = \|\mathbf{g}_1 \wedge \mathbf{g}_2\| \qquad (71) $$

The differential surface element can be expressed as:

$$ dS_0 = J_S da_1 da_2 \qquad (72) $$

and in particular, one has:

$$ \mathbf{n} \cdot dS_0 = \mathbf{g}_1 \wedge \mathbf{g}_2 \, da_1 \, da_2 \qquad (73) $$

Let us now suppose that the surface deforms due to displacements $\mathbf{u}$ with $\mathbf{y} = \mathbf{x} + \mathbf{u}$ being the actual coordinates. Let $\hat{\mathbf{g}}_\alpha$ denote the “actual” covariant vectors

$$ \hat{\mathbf{g}}_{\alpha} = \mathbf{y}_{,\alpha} = \mathbf{g}_{\alpha} + \mathbf{u}_{,\alpha} (= \mathbf{g}_{\alpha} + \nabla \mathbf{u} \cdot \mathbf{g}_{\alpha}) \qquad (74) $$
The vector product of the new covariant vectors gives the actual normal:

$$ \hat{\mathbf{g}}_1 \wedge \hat{\mathbf{g}}_2 = \mathcal{J}_S \hat{\mathbf{n}} \qquad (75) $$

We first compute the expansion of the vector product $\hat{\mathbf{g}}_1 \wedge \hat{\mathbf{g}}_2$:

$$ \begin{aligned} \hat{\mathbf{g}}_1 \wedge \hat{\mathbf{g}}_2 &= (\mathbf{g}_1 + \mathbf{u}_{,1}) \wedge (\mathbf{g}_2 + \mathbf{u}_{,2}) \\ &= \mathbf{g}_1 \wedge \mathbf{g}_2 + \mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2} \end{aligned} \qquad (76) $$

showing that the surface element $\mathcal{J}_S \hat{\mathbf{n}}$ admits an exact quadratic expansion in $\mathbf{u}$:

$$ \mathcal{J}_S \hat{\mathbf{n}} = J_S \mathbf{n} + \mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2} \qquad (77) $$

In particular, one has:

$$ \begin{aligned} \hat{\mathbf{n}} \mathrm{dS} &= \hat{\mathbf{g}}_1 \wedge \hat{\mathbf{g}}_2 \mathrm{d}a_1 \mathrm{d}a_2 = (J_S \mathbf{n} + \mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2}) \mathrm{d}a_1 \mathrm{d}a_2 \\ &= \left(\mathbf{n} + \frac{1}{J_S} \mathbf{g}_1 \wedge \mathbf{u}_{,2} + \mathbf{u}_{,1} \wedge \mathbf{g}_2 + \mathbf{u}_{,1} \wedge \mathbf{u}_{,2}\right) \mathrm{dS}_0 \end{aligned} \qquad (78) $$