Bifurcation analysis of quasi-periodic orbits of mechanical systems with 1:2 internal resonance via spectral submanifolds

We consider a periodically forced nonlinear mechanical system

where $\boldsymbol{x}\in\mathbb{R}^{n}$ is the generalized displacement vector; $\boldsymbol{M},\boldsymbol{C},\boldsymbol{K}\in\mathbb{R}^{n\times n}$ are the mass, damping and stiffness matrices; $\boldsymbol{f}(\boldsymbol{x},\dot{\boldsymbol{x}})$ is a $C^{r}$ smooth nonlinear function such that $\boldsymbol{f}(\boldsymbol{x},\dot{\boldsymbol{x}})\sim\mathcal{O}(|\boldsymbol{x}|^{2},|\boldsymbol{x}||\dot{\boldsymbol{x}}|,|\dot{\boldsymbol{x}}|^{2})$; and $\epsilon\boldsymbol{f}^{\mathrm{ext}}(\Omega t)$ denotes external harmonic excitation. We note that we allow for asymmetric damping and stiffness matrices to account for gyroscopic and follower forces, as well as nonlinear damping.

Solving the linear part of (1) leads to the generalized eigenvalue problem

where $\lambda_{j}$ is an eigenvalue, and $\boldsymbol{\phi}_{j}$ and $\boldsymbol{\theta}_{j}$ are associated right and left eigenvectors. We assume that the mechanical system (1) is in a 1:2 internal resonance; namely, there exist two pairs of complex conjugate modes that are 1:2 internally resonant. Without loss of generality, we let $\lambda_{1}$ and $\bar{\lambda}_{1}$ be the eigenvalues of the first pair of resonant modes, and $\lambda_{2}$ and $\bar{\lambda}_{2}$ be the eigenvalues of the second pair of resonant modes, we have

In addition, we assume that there are no resonances between the spectrum $\{\lambda_{1},\bar{\lambda}_{1},\lambda_{2},\bar{\lambda}_{2}\}$ and eigenvalues outside the spectrum.

We are interested in the forced dynamics when that the excitation frequency $\Omega$ is close to the natural frequency of the first pair of internally resonant modes, namely, $\mathrm{i}\Omega\approx\lambda_{1}$. Therefore, the first pair of modes is in external resonance. Due to the 1:2 internal resonance, energy transfer can occur between the two pairs of modes. We assume the origin of the system (1) is asymptotically stable, namely, $\mathrm{Re}(\lambda_{j})<0$ for $1\leq j\leq 2n$. Then the stable fixed point is perturbed as a stable small-amplitude periodic orbit when $\epsilon$ is sufficiently small guckenheimer2013nonlinear . However, when $\epsilon$ exceeds a critical value, the periodic orbit can become unstable via a secondary Hopf bifurcation (also referred to as Neimark-Sacker bifurcation), and quasi-periodic orbits appear. Our goal is to provide a thorough analysis of the bifurcations of the quasi-periodic orbits under the variations in $\epsilon$ and $\Omega$.

Solving quasi-periodic orbits of the mechanical system (1) is challenging for $n\gg 1$. We use reductions on spectral submanifolds (SSMs) to address this challenge. Next, we present a brief introduction to the theory of SSMs.

The second-order system (1) can be transformed into a first-order system as below

where

The associated generalized eigenvalue problem becomes (cf. (2))

where $\lambda_{j}$ is a generalized eigenvalue and $\boldsymbol{v}_{j}$ and $\boldsymbol{u}_{j}$ are the corresponding right and left eigenvectors, respectively. In particular, we have

In our setting, we take $\mathcal{E}=\mathrm{span}(\boldsymbol{v}_{1},\bar{\boldsymbol{v}}_{1},\boldsymbol{v}_{2},\bar{\boldsymbol{v}}_{2})$ as a master subspace for model reduction. It follows that we construct four-dimensional reduced-order models (ROMs) for the full system with a $2n$-dimensional phase space. A significant dimension reduction is obtained when $n$ is large. In the case of unforced system ($\epsilon=0$), the subspace $\mathcal{E}$ is invariant for the linear system $\boldsymbol{B}\dot{\boldsymbol{z}}=\boldsymbol{A}\boldsymbol{z}$, i.e., any trajectory of the linear system with an initial condition in $\mathcal{E}$ will stay in $\mathcal{E}$ for all times.

Under the addition of the nonlinear term $\boldsymbol{F}(\boldsymbol{z})$ in system (4), the master subspace $\mathcal{E}$ is perturbed into nonlinear normal modes (NNM), which are invariant manifolds tangent to $\mathcal{E}$ at the origin shaw1993normal . While there can be infinitely many NNMs associated with a given spectral subspace, there exists a unique, smoothest one among them under generically satisfied non-resonance conditions which is defined as the spectral submanifold (SSM) associated with $\mathcal{E}$. We denote this SSM as $\mathcal{W}(\mathcal{E})$. Slow SSMs are associated with the master subspace with the lowest linear decay rates and attract nearby full system trajectories exponentially fast. Hence, the reduced dynamics on slow SSMs provides a rigorous ROM. The SSM $\mathcal{W}(\mathcal{E})$ can be viewed as an embedding of an open set in reduced coordinates $\boldsymbol{p}=(p_{1},\bar{p}_{1},p_{2},\bar{p}_{2})\in\mathbb{C}^{4}$ via a map $\boldsymbol{z}=\boldsymbol{W}(\boldsymbol{p}):\mathbb{C}^{4}\to\mathbb{R}^{2n}$. In addition, the reduced dynamics on the SSM can be expressed as $\dot{\boldsymbol{p}}=\boldsymbol{R}(\boldsymbol{p})$.

Under further addition of the external periodic forcing $\epsilon\boldsymbol{F}^{\mathrm{ext}}(\Omega t)$, the autonomous SSM is perturbed as a unique, smoothest, periodic SSM, which we denote as $\mathcal{W}(\mathcal{E},\Omega t)$. Similarly to the autonomous SSM, the time-periodic SSM $\mathcal{W}(\mathcal{E},\Omega t)$ can be embedded via a map $\boldsymbol{W}_{\epsilon}(\boldsymbol{p},\Omega t):\mathbb{C}^{4}\times\mathbb{S}^{1}\to\mathbb{R}^{2n}$. The reduced dynamics on this time-periodic SSM is non-autonomous and can be expressed as jain2022HowComputeInvariant ; li2022NonlinearAnalysisForceda

We express the reduced coordinates in polar form as

and adopt a leading-order approximation to the non-autonomous part $\boldsymbol{S}(\boldsymbol{p},\Omega t)$  li2022NonlinearAnalysisForceda . Then the reduced dynamics is given as li2022NonlinearAnalysisForceda

where $\boldsymbol{\rho}=(\rho_{1},\rho_{2})$, $\boldsymbol{\theta}=(\theta_{1},\theta_{2})$, $f$ is a modal force, and $g_{i}(\boldsymbol{\rho},\boldsymbol{\theta})$ for $1\leq i\leq 4$ are nonlinear functions. Explicit expressions for $f$ and $g_{i}$ can be found in eqs. (37)-(38) of li2022NonlinearAnalysisForceda .

Due to the transformation (9), a fixed point of the vector field (10) corresponds to a periodic orbit $\boldsymbol{p}_{\mathrm{po}}(t)$ on the SSM. Therefore, we can extract the forced response of periodic orbits from the zero-level contours of the vector field (10). We can further infer bifurcations of periodic orbits of the full system from bifurcations of fixed points of the vector field (10) liNonlinearAnalysisForced2022c . In addition, a limit cycle of the SSM-based ROM (10) corresponds to a two-dimensional invariant torus of the full system liNonlinearAnalysisForced2022c . Hence, we can predict bifurcations of quasi-periodic orbits that stay on tori via the bifurcations of the limit cycles of the ROM.

As summarized in Fig. 1, SSM-based reductions enable effective bifurcation analysis of nonlinear dynamics of high-dimensional mechanical systems. Indeed, it achieves a significant dimension reduction because the $2n$-dimensional full system is reduced to a 4-dimensional ROM (10). Furthermore, thanks to the normal form parameterization style used in the polar ROM (10), a further simplification is obtained, i.e., quasi-periodic orbits of the full system are calculated as limit cycles of the ROM.

Inspired by the above discussion, we will present bifurcation analysis of quasi-periodic orbits in the Poincaré section of the full system (4). In particular, we take $\{(\boldsymbol{z},t):\mathrm{mod}(t,2\pi/\Omega)=0\}$ as the Poincaré section. Let $T=2\pi/\Omega$, the associated Poincaré map maps the state $\boldsymbol{z}(t)$ from $t=kT$ to $t=(k+1)T$ under the flow generated by (4). Therefore, we also call this map as period-$2\pi/\Omega$ map. We further define Poincaré intersection for an orbit/attractor as the intersection of the orbit/attrator with the Poincaré section.

We note that a chaotic attractor of the 4-dimensional SSM-based ROM (10) also corresponds to a chaotic attractor of the full system (4), as pointed out in Fig. 1. We will use the aforementioned Poincaré intersection as a tool to distinguish chaotic attractors from quasi-periodic attractors. In particular, the Poincaré intersection of a quasi-periodic attractor will be a closed curve, while the Poincaré intersection of a chaotic attractor will be a strange attractor (cf. Fig. 16). As an additional approach, we will also use the power spectral density (PSD) plot to distinguish them. The power spectral density (PSD) plot of a quasi-periodic attractor has isolated peaks (cf. Fig. 28), while the PSD plot of a chaotic attractor will have significant values for a broad frequency bandwidth (cf. the middle panel of Fig. 11).

We still need to determine the map $\boldsymbol{W}_{\epsilon}(\boldsymbol{p},\Omega t)$ and the unknown coefficients for the nonlinear functions $g_{i}$ in (10). For the consistency with (8), we approximate $\boldsymbol{W}_{\epsilon}(\boldsymbol{p},\Omega t)$ as

We further adopt a leading-order approximation to the non-autonomous part $\boldsymbol{X}(\boldsymbol{p},\Omega t)$, namely,

Here $\boldsymbol{X}_{0}(\Omega t)=\boldsymbol{x}_{0}e^{\mathrm{i}\Omega t}+\bar{\boldsymbol{x}}_{0}e^{-\mathrm{i}\Omega t}$ and the coefficient vector $\boldsymbol{x}_{0}$ is obtained by solving a system of linear equations jain2022HowComputeInvariant ; li2022NonlinearAnalysisForceda .

We use the computational procedure  jain2022HowComputeInvariant based on the parametrization method CabreIII ; Haro2016TheManifolds to solve for the autonomous map $\boldsymbol{W}(\boldsymbol{p})$ and as well the unknown coefficients in $g_{i}$ (see eq. (10)). The map $\boldsymbol{W}(\boldsymbol{p})$ is approximated by a truncated Taylor expansion in $\boldsymbol{p}$. Then the coefficients of the Taylor expansion, along with $g_{i}$, are determined from the invariance of SSM. In particular, they are solved from systems of linear equations in a recursive fashion, which allows us to automatically compute the SSM expansion up to arbitrary polynomial orders without change of coordinates using only the eigenvectors associated with the master subspace $\mathcal{E}$  jain2022HowComputeInvariant . In practice, we truncate the Taylor expansion based on the convergence of forced response with increasing expansion orders li2022NonlinearAnalysisForceda . We use notation $\mathcal{O}(q)$ with $q\in\mathbb{N}$ to denote the Taylor expansion is truncated at the $q$-th order in this paper.

The above parameterization method has been implemented in SSMTool ssmtool2 , an open-source package for the computation of invariant manifolds. We note that coco COCO ; dankowicz2013recipes ; ahsan2022methods , a continuation package, has been integrated into SSMTool li2022NonlinearAnalysisForceda ; liNonlinearAnalysisForced2022c . This integration enables the detection of bifurcations of fixed points of (10), parameter continuation of these bifurcated fixed points, as well as the continuation of (bifurcated) limit cycles of the ROM (10). In particular, we use the po-toolbox of coco to obtain the limit cycles. In the po-toolbox, boundary-value problems are formulated for the limit cycles, and collocation methods are used to solve for the limit cycles and their associated monodromy matrices. The stability and bifurcation of these limit cycles are further analyzed using Floquet theory. Therefore, we can easily perform bifurcation analysis of quasi-periodic orbits using SSMTool. Indeed, we will use this package to study the bifurcation of quasi-periodic orbits of a suite of 1:2 internally resonant mechanical systems in the coming sections, ranging from simple oscillators to high-dimensional FE models.

The governing equations of two coupled nonlinear oscillators are

The undamped natural frequencies of the linearized system are given as $\omega_{1}=1$ rad/s and $\omega_{2}=2$ rad/s. Thus, the system has a 1:2 internal resonance. In the following computations, system parameters are chosen as $c_{1}=0.005~{}\rm{N\cdot s/m}$, $c_{2}=0.01~{}\rm{N\cdot s/m}$, $b_{1}=b_{2}=1~{}\rm{N/m^{2}}$, $f_{1}=1~{}\rm{N}$, $f_{2}=0~{}\rm{N}$, and $\epsilon=0.01$, unless otherwise stated.

We first present the periodic and quasi-periodic orbits of the system under varying $\Omega$ obtained via reduction on SSM at $\mathcal{O}({3})$ truncation. We recall that these results for $x_{1}$ have been presented and also verified in liNonlinearAnalysisForced2022c . Here, we present the results for $x_{1}$ to make this paper self-contained. We will also present the forced response curve (FRC) of periodic and quasi-periodic orbits for $x_{2}$.

The FRCs of periodic orbits (FRC-PO) of the two oscillators for $\Omega\in$ [0.7, 1.1] are shown in Fig. 2. Here, $\|x_{1}\|_{\infty}$ and $\|x_{2}\|_{\infty}$ represent the infinite norm of $x_{1}(t)$ and $x_{2}(t)$, giving the amplitude of periodic or quasi-periodic response. There are four saddle-node (SN) bifurcation points and two secondary Hopf bifurcation (HB) points on the FRC-PO. In particular, the periodic orbits for $\Omega\in[\Omega_{\mathrm{HB1}},\Omega_{\mathrm{HB2}}]$ are unstable.

The secondary Hopf bifurcation marks the birth of quasi-periodic orbits. The FRCs of quasi-periodic orbits (FRC-QO) of the two oscillators under variations of $\Omega$ are shown in Fig. 3. We recall that these quasi-periodic orbits stay on some two-dimensional tori, which are computed as limit cycles of the ROM (10) (cf. Fig. 1). Therefore, the FRC-QO is obtained by parameter continuation of limit cycles that are born out of HB1. Along the FRC-QO, 20 period-doubling (PD) bifurcation points and 20 SN bifurcation points are observed, marking the complex bifurcations of quasi-periodic orbits. We note that the appearance of PD bifurcation points indicates that new families of quasi-periodic orbits with doubled periods can be obtained. We will analyze these new families of quasi-periodic orbits in Sect. 3.5.

Now we allow for the variations in the forcing amplitude $\epsilon$ to study how the PD and SN bifurcated quasi-periodic orbits in Fig. 3 evolve with varying $\epsilon$. We expect that they will disappear if $\epsilon$ is sufficiently small. We note that both the PD and SN bifurcations are codimension-one bifurcation. Thus we will have a curve of PD or SN bifurcated quasi-periodic orbits under the variations in $(\Omega,\epsilon)$. Since the PD and SN bifurcations of quasi-periodic orbits correspond to PD and SN bifurcations of the limit cycles of the ROM (10), we use parameter continuation of PD or SN bifurcated limit cycles of the ROM to extract the bifurcation curve. In particular, we initialize our continuation runs with the bifurcation points on FRC-QO shown in Fig. 3.

The continuation paths for the PD and SN bifurcation points with $(\Omega,\epsilon)\in$ [0.7, 1.1] $\times$ [0.0001, 0.03] are shown in the left and right panels of Fig. 4, respectively. As the upper bound of $\epsilon$ is increased to 0.03, we use $\mathcal{O}(7)$ truncation for SSM reduction to ensure converged predictions. Here, the lower panels are the projection of the corresponding upper panels. We observe from the figure that multiple branches exist for the PD or SN continuation path. The amplitudes of SN and PD bifurcated quasi-periodic orbits decrease as $\epsilon$ decreases, as seen in the upper two panels. Moreover, some PD or SN branches merge and then disappear as $\epsilon$ decreases, resulting in a dynamic change in the number of PD or SN bifurcation points on FRC-QO. For example, when $\epsilon$ is equal to 0.01, there are 20 SN and 20 PD bifurcation points, while at $\epsilon$ = 0.005, there are only 10 SN and 10 PD bifurcation points, and at $\epsilon$ = 0.0135, both numbers are increased to 22.

We observe from the right panel of Fig. 3 that the amplitude of the quasi-periodic orbits is close to that of the unstable periodic orbits for $\Omega\approx 1$. If the forcing amplitude $\epsilon$ is increased, it may occur that the quasi-periodic orbits intersect with the periodic orbits. In this subsection, we will discuss how the intersection of quasi-periodic and periodic orbits leads to the emergence of homoclinic bifurcations of quasi-periodic orbits.

We first compute the FRC-PO and FRC-QO with varying $\Omega$ but fixed $\epsilon=0.03$. The obtained results are plotted in Fig. 5. Here, the continuation run of quasi-periodic orbits starting at HB1 is terminated at point A of FRC-QO shown in Fig. 5, where the amplitude $||x_{2}||_{\infty}$ barely changes. We note that each quasi-periodic orbit on the FRC-QO consists of two incommensurable base frequencies: one is the excitation frequency $\Omega$, and the other one is much slower and given by $2\pi/T_{s}$, where $T_{s}$ is the period of the corresponding limit cycle of the SSM-based ROM (10) liNonlinearAnalysisForced2022c . As seen in the lower panel of Fig. 5, $T_{s}\to\infty$ as the continuation approaches point A. Therefore, the quasi-periodic orbit undergoes a infinite-period bifurcation Strogatz2018Nonlinear when $\Omega\to\Omega_{\mathrm{A}}$. As mentioned in Strogatz2018Nonlinear , there are two kinds of infinite-period bifurcation: saddle-node on invariant circle bifurcation (SNIC) and homoclinic bifurcation. Next, we show that our observed infinite-period bifurcation is the second kind.

To have a close look at the infinite-period bifurcation, we compute the intersection of the period-$2\pi/\Omega$ Poincaré map for some of these quasi-periodic orbits as $\Omega\to\Omega_{\mathrm{A}}$. We note that these quasi-periodic orbits stay on some two-dimensional invariant tori, and the intersection of these tori with the Poincaré section results in closed curves, which can be obtained by mapping the limit cycles of the SSM-based ROM (10) back to physical coordinates liNonlinearAnalysisForced2022c . The obtained results are presented in Fig. 6, from which we see that the limit cycle with an infinite period has a kink. This kink corresponds to a saddle fixed point. Therefore, we infer that the infinite-period bifurcation corresponds to a homoclinic bifurcation of quasi-periodic orbits.

A saddle fixed point on the Poincaré section corresponds to a periodic orbit of the full system with period-$2\pi/\Omega$. We plot the projection of this periodic orbit shown in the upper panel of Fig. 7. In fact, this periodic orbit is the one marked by B on the FRC-PO given by the upper panel of Fig. 5. The two markers A and B in the upper panel of Fig. 5 have not coincided because the vertical axis denotes the amplitude of periodic or quasi-periodic orbits. Indeed, as seen in the upper panel of Fig. 7, the intersection point for the periodic orbit B at the Poincaré section does not coincide with the peak position of the orbit. We plot the torus A along with the periodic orbit B in the lower panel of Fig. 7 for a direct visualization of the intersection, which illustrates clearly the global bifurcation, i.e., homoclinic bifurcation of quasi-periodic orbit.

We conclude this subsection by investigating how this global bifurcation evolves as $\epsilon$ changes. We expect that the homoclinic bifurcation will disappear when $\epsilon$ is less than a critical value. We perform parameter continuation of the homoclinic orbit under the variations in $(\Omega,\epsilon)$ liNonlinearAnalysisForced2022c ; li2023nonlinear . In particular, we fix the period of the associated limit cycle of the SSM-based ROM to be 5000, perform continuation of the limit cycle, and then map these limit cycles back to homoclinic bifurcated quasi-periodic orbits of the original system. The obtained continuation path of the homoclinic bifurcations is shown in Fig. 8 (marked as Tinf curve). We observe that the continuation path starting at point A terminates at point C, suggesting that there are two homoclinic bifurcations when $\epsilon=0.03$. In fact, this point C can also obtained from FRC-QO starting at HB2, as seen in Fig.  5. We also find that there is no limit cycle with an infinite period when $\epsilon$ is less than 0.0232. Thus, the critical value of $\epsilon$ for the homoclinic bifurcation is 0.0232, below which there is no homoclinic bifurcations.

Now we revisit Fig. 3 and discuss the periodic doubling bifurcation of quasi-periodic orbits. In particular, we extract the new families of quasi-periodic orbits born out of these PD bifurcations. For compactness, we focus on the region where $\Omega\in[\Omega_{\mathrm{PD9}},\Omega_{\mathrm{PD10}}]$. We switch to the continuation of limit cycles born out of PD9 of the ROM (10), and then map them back to quasi-periodic orbits of the original system. The FRC of these quasi-periodic orbits is plotted in red lines in Fig. 9 and is denoted by ’Period2’. As expected, it also intersects with the main FRC-QO (magenta lines) at PD10. Along this FRC, two new PD bifurcation points are detected at $\Omega\approx 1.0025$ and $\Omega\approx 1.003$, marking the existence of a new family of periodic orbits of doubled period. We switch the parameter continuation to this new family of quasi-periodic orbits and obtain the FRC plotted in blacked lines in Fig. 9. The internal period-$T_{s}$ gets further doubled for these quasi-periodic orbits, and hence they are denoted as ’Period4’ in the figure. Again, we detect two PD bifurcation points along the FRC of ’Period4’ quasi-periodic orbits. In short, the system displays a cascade of period-doubling bifurcations for the quasi-periodic orbits between the frequency internal $\Omega\in[\Omega_{\mathrm{PD9}},\Omega_{\mathrm{PD10}}]$.

A cascade of period-doubling bifurcations is a typical route to chaos. To have a close look at the transition from quasi-periodic to chaotic attractors, we present the bifurcation diagram for the attractors of the SSM-based ROM (10) for $\Omega\in[1.0023,1.0032]$ (cf. Fig. 9) in Fig. 10. This diagram is obtained via performing forward simulations of the ROM (10) at given sampled frequencies $\Omega$, extracting corresponding attractors in steady state, and locating the intersections of these attractors with appropriately constructed Poincaré map. This bifurcation diagram displays the transition to chaos via the cascade of period-doubling bifurcations. This bifurcation diagram is consistent with that of Fig. 9 but provides a more complete picture of the complex dynamics in this parameter region.

To further check these strange attractors are indeed chaotic, we randomly select one representative at $\Omega$ = 1.0027 to conduct a detailed analysis. We map the attractor back to physical coordinates and then plot the intersection of this attractor with the Poincaré section induced by the period-$2\pi/\Omega$ map in the upper panel of Fig. 11. As seen in the middle panel of the figure, the power spectral density of the intersected trajectory has significant values for a broad frequency bandwidth. We observe from the right panel of Fig. 11 that the maximum Lyapunov exponent of the attractor is 0.001. These three panels together indicate that the attractor is indeed chaotic. Here, the Lyapunov exponent is obtained via a discrete method detailed in geist1990comparison . In this method, variational equations are integrated, and QR decomposition is utilized to perform reorthonormalization.

To verify the effectiveness of the SSM-based prediction, we use numerical integration of the original system to generate reference results for the purpose of comparison. We take points D and E in Fig. 9 as representatives of Period2 and Period4 quasi-periodic orbits and plot their Poincaré intersections in Fig. 27 of Appendix section. The associated power spectral density (PSD) plots for these two quasi-periodic orbits are presented in Fig. 28. We observe that the results from the numerical integration match well with that of SSM-based prediction. As for the chaotic attractor, we follow a similar validation procedure using numerical integration, and the results are presented in Fig. 11. The original system indeed admits a chaotic attractor that coincides well with the SSM-based prediction.

Let $Z_{0}(x)$ be the initial, imperfect configuration of the curved beam, $w$ be the deflection relative to the initial configuration, the governing equation of the curved beam is given by oz2006TwotooneInternalResonances

along with boundary conditions

Here, $EI$ is flexural rigidity, $\rho$ is density, $A$ is the area of the cross-section, and $L$ is the length of the beam. We define the following dimensionless quantities wang2012DynamicsSimplySupported

and then (14) can be rewritten as

along with boundary conditions

We apply a Galerkin approach to discretize the partial-integral differential equation (17). Specifically, we substitute

into (17) and then apply a Galerkin projection, yielding a system of second-order ordinary differential equations below

where $\boldsymbol{q}\in\mathbb{R}^{n}$ is a generalized coordinate vector $\boldsymbol{M}$, $\boldsymbol{C}$, $\boldsymbol{K}$, $\boldsymbol{H}$, and $\boldsymbol{B}\in\mathbb{R}^{n\times n}$ are constant matrices, $\boldsymbol{d},\boldsymbol{f}\in\mathbb{R}^{n}$ and are constant vectors. Detailed expressions for these constant matrices and vectors are presented in (22) in the Appendix.

We use a parabolic function, namely $Z_{0}(x)=4a_{0}x\\ (L-x)$, to characterize the initial imperfect configuration. The associated $\varphi_{0}(\xi)=4a_{0}\xi(1-\xi)$. We tune $a_{0}$ such that the ratio of the first two undamped natural frequencies is 1:2. We take a five-mode truncation, i.e., $n=5$, because it is sufficient to generate converged solutions in our setting. Thus, the state space of the full system is 10, which is reduced to 4 via the SSM-based ROM (10). Numerical experiments show that the natural frequencies of the first two modes are given as $\omega_{1}=39.4784$ and $\omega_{2}=78.9580\approx 2\omega_{1}$ when $a_{0}$ = 14.2365. In the following computations, other system parameters are chosen as $\bar{\mu}=0.01$, $\kappa=1$, $\epsilon=0.01$, and $\bar{F}=5$, unless otherwise stated.

We first compute the forced response curve of periodic orbits for the system with $\Omega\in[39.2,39.7]$. Similar to the previous example, this FRC-PO is obtained by parameter continuation of the fixed point of the SSM-based ROM (10). In this example, an $\mathcal{O}(3)$ expansion is used because it is sufficient to yield converged solutions for the selected parameters. The obtained FRC-PO is shown in Fig. 13. Here $\|q_{1}\|_{\infty}$ and $\|q_{2}\|_{\infty}$ represent the infinite norm of the modal coordinates $q_{1}(\tau),q_{2}(\tau)$, giving the amplitude of the periodic response. Similar to the system of two coupled oscillators, there are two HB points (HB1 and HB2) and four SN bifurcation points on the FRC-PO. Moreover, the periodic motion for $\Omega\in[\Omega_{\mathrm{HB1}},\Omega_{\mathrm{HB2}}]$ is unstable, and we will focus on quasi-periodic motions born out of these two HB points of the system.

We switch to the parameter continuation of limit cycles of the SSM-based ROM (10) at HB1 and then extract the forced-response curve of quasi-periodic orbits (FRC-QO) under variations in $\Omega$. This continuation run proceeds until it reaches HB2. The obtained FRC-QO along with FRC-PO are shown in Fig. 14. 20 PD bifurcation points and 20 SN bifurcation points are detected on the FRC-QO, leading to the coexistence of stable and unstable quasi-periodic motions. In addition, we also observe the coexistence of periodic and quasi-periodic attractors in Fig. 14. For example, coexisting periodic attractor (Attractor1) and quasi-periodic attractor (Attractor2) can be found at $\Omega$ = 39.503. The validation of the SSM-based predictions for the FRC-PO and FRC-QO is provided in Appendix B, as seen in Fig. 29 and Figs. 30-32.

Since PD bifurcation points were detected, we expect that this curved beam system also has a cascade of PD bifurcations for the quasi-periodic orbits, leading to chaotic motion. This is indeed the case here. We focus on the region $\Omega\in[\Omega_{\rm{PD9}},\Omega_{\rm{PD10}}]$ and $[\Omega_{\rm{PD11}},\Omega_{\rm{PD12}}]$ to illustrate the cascade of PD bifurcations. We switch to the parameter continuation of limit cycles of the ROM (10) with doubled period (Period2) at PD9 and obtain the continuation path in red lines in the upper panel of Fig. 15. Along this Period2 FRC-QO, two new PD bifurcation points are detected, indicating the appearance of quasi-periodic motions with a quadrupled period (Period4). We switch to the continuation of these Period4 tori and obtain the Period4 FRC-QO plotted as black lines in the figure. We once again detect two new PD points on this black line, which illustrates the cascade of PD bifurcations. This cascade of PD bifurcations is also observed for $\Omega\in[\Omega_{\rm{PD11}},\Omega_{\rm{PD12}}]$, as can be seen in the lower panel of Fig. 15, where the red and black lines again denote the FRCs of Period2 and Period4 quasi-periodic orbits.

The above-observed cascades of PD bifurcations are expected to lead to chaotic motions. To confirm it, we check whether the system (20) exists chaotic attractors at two sampled excitation frequencies $\Omega=39.4745$ and $\Omega=39.4834$. Specifically, we perform forward simulations to the SSM-based ROM (10) and then map the generated trajectories back to the full system. Interestingly, we find that the steady state of the forward simulations can be either chaotic or quasi-periodic, depending on the choice of initial conditions. This differs from the previous example’s case (cf. Sect. 3.5). We present the Poincaré intersections of these coexisting quasi-periodic and chaotic attractors with period-$2\pi/\Omega$ map in Fig. 16, where the upper two panels show the two coexisting attractors at $\Omega=39.4745$, while the lower two panels display the two coexisting attractors at $\Omega=39.4834$. The two strange attractors in the right panels of Fig. 16 are indeed chaotic because their associated maximum Lyapunov exponents are 0.0038 and 0.0041, as detailed in Fig. 34.

We use numerical integration of the full system to validate the predicted cascade of PD bifurcations and the coexisting quasi-periodic and chaotic attractors. We take points A and C in Fig. 15 as representative of Period2 quasi-periodic orbits, and points B and D in Fig. 15 as representative of Period4 quasi-periodic orbits. We plot the Poincaré section of these quasi-periodic orbits in Fig. 33 in Appendix, from which we see that SSM-based predictions match perfectly with the reference solutions obtained via the numerical integration. In addition, as seen in Fig. 16, the predicted coexisting quasi-periodic and chaotic attractors are also successfully validated using the numerical integration of the full system. Further, we observe from Fig. 34 in Appendix that the predicated maximum Lyapunov exponents match well with the reference solutions as well. These demonstrate the effectiveness of the SSM-based ROM (10).

We conclude this subsection with bifurcation diagrams of attractors on the region $\Omega\in[\Omega_{\rm{PD9}},\Omega_{\rm{PD10}}]$ and $[\Omega_{\rm{PD11}},\Omega_{\rm{PD12}}]$. In particular, we follow the same method to produce Fig. 10 to extract the bifurcation diagram of attractors of the SSM-based ROM (10). The obtained bifurcation diagrams are shown in Fig. 17, from which we see the complete picture of the transition to chaos via the cascade of period-doubling bifurcations. These bifurcation diagrams are also consistent with that of Fig. 15.

The coexisting quasi-periodic and chaotic attractors at $\Omega=39.4745$ and $\Omega=39.4834$ indicate that the curved beam admits stable quasi-periodic orbits at these two sampled excitation frequencies. As seen in Fig. 15, the quasi-periodic orbits for $\Omega=39.4745$ and $\Omega=39.4834$ on the main branch of FRC-QO and its bifurcated branches with Period2 and Period4 quasi-periodic orbits (red and black lines) are all unstable. This implies that isolated branches of quasi-periodic orbits exist. To locate the isolas, we recall that each quasi-periodic orbit corresponds to a limit cycle of the SSM-based ROM (10) (cf. Fig. 1). So we perform two continuation runs of limit cycles of the ROM (10) under the variation in $\Omega$. These continuation runs are initialized with the limit cycles at $\Omega=39.4745$ and $\Omega=39.4834$ respectively. Indeed, we locate an isola in each of these continuation runs and then map the isolas of the two continuation runs back to the full system, yielding the isolas of quasi-periodic orbits that are marked as Isoa1 and Isola2 in Fig. 15. We observe that each of these two isolas has 2 PD bifurcation points and 4 SN bifurcation points. Moreover, the two quasi-periodic attractors in Fig. 16 correspond to point E on Isola1 and point F on Isola2, as seen in Fig. 15.

Next, we study how these isolas evolve under the variations in the forcing amplitude $\epsilon$. We perform the continuation of SN bifurcated quasi-periodic orbits with SN1 or SN7 as the initial solution. We restrict to the computational domain $\epsilon\in[0,0.015]$. The projection of the SN bifurcation curves of these two continuation runs are shown in Fig. 18. Here, the upper and lower panels show the SN bifurcation curves with SN1 and SN7 respectively. Interestingly, we observe from the upper panel that there exist two more SN points, i.e., SN5 and SN6, other than the four SN points on Isola1 (cf. the upper panel of Fig. 15), at $\epsilon=0.01$. This indicates the existence of other isolas at $\epsilon=0.01$. Indeed, we take SN5 as an initial point and perform continuation of quasi-periodic orbits with fixed $\epsilon$ and varying $\Omega$, resulting in the Isola3 in the upper panel of Fig. 15. Similarly, we observe from the lower panel that there are two more SN points, i.e., SN11 and SN12, other than the four SN points on Isola2 (cf. the lower panel of Fig. 15.), at $\epsilon=0.01$. We follow the same procedure and locate the Isola4 in the lower panel of Fig. 15.

We can further infer cusp, simple, and isola bifurcations from the continuation path shown in Fig. 18. Indeed, the kink points in Fig. 18 mark cusp bifurcations where two SN points merge and then disappear. Meanwhile, we see the emergence of two SN points at which $\partial\epsilon/\partial\Omega=0$ along the continuation path in Fig. 18. These critical points where $\partial\epsilon/\partial\Omega=0$ holds are either isola bifurcation points on which isolas are born or simple bifurcation points where two isolas merge. To illustrate them more clearly, we plot these continuation paths in $(\Omega,\epsilon,||q_{2}||_{\infty})$ along with some sampled FRC-QO in Fig. 19. We infer from these plot that Isola1 (Isola3) is born out of the isola bifurcation point IBA (IBB). In addition, Isola1 and Isola3 merge as a single isola via the simple bifurcation point SBA, and then the single isola persists as $\epsilon$ increases, as seen in the upper panel of Fig. 19. We observe similar bifurcation patterns for the lower panels of Fig. 18 and Fig. 19. In particular, Isola2 and Isola4 are born out of IBC and IBD and then merge via SBB.

Following li2022NonlinearAnalysisForceda , the geometric parameters of the shell in Fig. 20 are set to be: the length $L$ = 2 m, width $H$ = 1 m, thickness $t$ = 0.01 m, and the curvature parameter $w$ = 0.041 m. Material properties are specified with the density $\rho$ = 2700 kg/m³, Young’s modulus $E$ = 70 $\times$ 10⁹ Pa, and Poisson’s ratio $\nu$ = 0.33. The forcing amplitude is given as $F$ = 100 N.

The shallow shell model is discretized using flat, triangular shell elements. The FE model here contains 400 elements, and each node in the elements has six DOFs, resulting in $n$ = 1320 DOF li2022NonlinearAnalysisForceda . Thus, the dimension of state space is reduced from 2640 for the full system to 4 of SSM-based ROM (10). The ODEs of the system obtained by the finite element method can be expressed as

where $\boldsymbol{z}\in\mathbb{R}^{1320}$ is a displacement vector, $\boldsymbol{M},\boldsymbol{C},\boldsymbol{K}$ are the mass, damping, and stiffness matrices, $\boldsymbol{f}(\boldsymbol{z})$ is an analytic vector-valued function characterizing nonlinear internal force and $\boldsymbol{f}^{\mathrm{ext}}$ is the external load vector. Here we adopt a Rayleigh damping model such that the damping matrix can be expressed as $\boldsymbol{C}=\alpha\boldsymbol{M}+\beta\boldsymbol{K}$. The natural frequencies of the first two modes are given as $\omega_{1}=149.22$ and $\omega_{2}=298.78\approx 2\omega_{1}$. We choose $\alpha$ and $\beta$ such that the damping ratios of the first two modes are equal to 0.001. In the following computations, we set $\epsilon=0.03$ unless otherwise stated. We take the transverse displacements at two nodes with coordinates ($x$, $y$) = (0.25$L$, 0.5$H$) and (0.5$L$, 0.5$H$) to characterize the nonlinear vibration of the shell. They are denoted as $z_{\rm{out1}}$ and $z_{\rm{out2}}$ in this example.

We first compute the FRC-PO with $\Omega\in[144,153]$. This FRC-PO is obtained by parameter continuation of the fixed point of the SSM-based ROM (10). In this example, an $\mathcal{O}(5)$ expansion is used because it is sufficient to yield converged solutions for the selected parameters. Here, $\|z_{\rm{out1}}\|_{\infty}$ and $\|z_{\rm{out2}}\|_{\infty}$ represent the infinite norm of $z_{\rm{out1}}$, $z_{\rm{out2}}$, giving the amplitude of periodic or quasi-periodic response of the system. The FRC-PO is shown in Fig. 21. Similar to the previous two examples, there are two Hopf bifurcation (HB) points and four saddle-node (SN) bifurcation points on the FRC-PO. Since the periodic motion for $\Omega\in[\Omega_{\rm{HB1}},\Omega_{\rm{HB2}}]$ is unstable, we next analyze the quasi-periodic motion of the system between these two HB points.

The FRC-QO under varying $\Omega$ is shown in Fig. 22. Along this curve of quasi-periodic orbits, 12 periodic doubling (PD) bifurcation points and 6 saddle-node (SN) bifurcation points are detected. Similar to the previous two examples, the coexistence of stable and unstable quasi-periodic motions can be found. Moreover, we observe the coexistence of periodic and quasi-periodic attractors for $\Omega<\Omega_{\mathrm{HB1}}$ and $\Omega>\Omega_{\mathrm{HB2}}$ in Fig. 22.

We provide a validation of the SSM-based predictions for the FRC-PO is provided in Appendix C. In particular, we use a COCO-based shooting toolbox to extract the FRC-PO. As seen in Fig. 35, the FRC of periodic orbits obtained from SSM-based reduction matches well with that of the full system obtained via the shooting method. As detailed in Appendix C, a significant speed-up gain is obtained using the SSM-based reduction. Specifically, the computational times for extracting the FRC-PO of the SSM-based prediction and the shooting method are about 1 minute and 2 days, respectively.

We note that it is challenging to extract quasi-periodic orbits of this high-dimensional full system. As detailed in Appendix C, we take two representative stable tori, i.e., points A and B in Fig. 22, for the sake of verification. Specifically, we take a point on the invariant torus A or B obtained from SSM-based prediction as an initial condition, perform a forward simulation of the full system with the initial condition, and extract the associated quasi-periodic attractor in steady state. As seen in the upper and middle panels of Fig. 36, the SSM-based predictions match well with that of the reference solutions of the full system. Here, we again achieve a significant speed-up gain. Indeed, a forward simulation of the full system with 4500 cycles to ensure that the steady state arrived took more than 4 days of computational time. In contrast, we can obtain the whole FRC-QO shown in Fig. 22 in less than 4 minutes.

Similar to the case of two coupled oscillators, we perform parameter continuation of PD and SN bifurcated quasi-periodic orbits on the FRC-QO shown in Fig. 22 to study how these PD and SN bifurcated quasi-periodic orbits evolve with varying $\epsilon$. In particular, we take the PD and SN points on the FRC-QO in Fig. 22 as the initial solutions of the continuation runs. The projection of the continuation paths for the PD and SN bifurcation points with $(\Omega,\epsilon)\in[144,153]\times[0.0001,0.06]$ are shown in the left and right panels of Fig. 23, respectively. Here, the lower panels are the projection of the corresponding upper panels onto the parameter plane $(\Omega,\epsilon)$. The amplitudes of SN and PD bifurcated quasi-periodic orbits is reduced as $\epsilon$ decreases, as seen in the upper two panels. We observe that there are 6 branches for the PD continuation path and 3 branches for the SN continuation path. Further, the number of PD or SN bifurcation points changes dynamically as $\epsilon$ varies. In particular, these bifurcations disappear for sufficiently small forcing amplitude $\epsilon$. These features are similar to those in the system of two coupled oscillators.

Similarly to the previous two mechanical systems, we expect that there is a cascade of PD bifurcations for quasi-periodic orbits because PD bifurcation points were detected in Fig. 22. Indeed, the cascade of PD bifurcations is found on the region $\Omega\in[\Omega_{\rm{PD5}},\Omega_{\rm{PD6}}]$. By switching to parameter continuation of quasi-periodic orbits with doubled internal period (Period2) at PD5, we obtain the continuation path shown as the red line in Fig. 24. Along this Period2 FRC-QO, two new PD bifurcation points are detected. By switching to parameter continuation of quasi-periodic orbits with a quadrupled internal period (Period4) at one of the PD points on the Period2 FRC-QO, we obtain the continuation path of quasi-periodic orbits shown as a black line in Fig. 24. We once again detect two new PD points on this black curve of FRC-QO, which illustrates the cascade of PD bifurcations.

In the previous two mechanical systems, a cascade of PD bifurcations finally leads to chaotic motions. To test whether this is a common feature for mechanical systems with 1:2 internal resonance, we perform forward simulations of the SSM-based ROM (10) to extract attractors for $\Omega\in[149.388,149.438]$. We present the obtained bifurcation diagram for these attractors of the SSM-based ROM in Fig. 25 (cf. Fig. 24). This diagram is obtained following a similar procedure as that of Fig. 10. This bifurcation diagram displays the transition to chaos via the cascade of period-doubling bifurcations. This bifurcation diagram is consistent with that of Fig. 24 but provides a more complete picture of the complex dynamics in this parameter region. Interestingly, we observe a large window of period-3 limit cycles of the SSM-based ROM (10). These period-3 limit cycles correspond to period-3 tori of the full system.

We again have a close look at a representative chaotic attractor at $\Omega=149.41$. We map the simulated attractor back to physical coordinates. The intersections of the simulation results with the Poincaré section induced by the period-$2\pi/\Omega$ map are shown as a blue line on the upper panel of Fig. 26. The associated plot of PSD is shown in the lower panel of Fig. 26. We follow the same computational method as in the right panel of Fig. 11 to compute the maximum Lyapunov exponent of this strange attractor. The obtained maximum Lyapunov exponent is 0.0534. Therefore, a chaotic attractor is indeed observed in steady state.

To verify the effectiveness of the SSM-based predictions for the cascade of PD bifurcations, we take representative Period2 and Period4 quasi-periodic orbits, namely, points C and D on the FRC-QO shown in Fig. 24, and also the chaotic attractor in Fig. 26 to perform the verification. Specifically, we apply forward simulations of the full system to obtain reference solutions. More details on the forward simulations can be found in Appendix C. As seen in the lower two panels of Fig. 36, SSM-based predictions for the Period2 quasi-periodic orbit C and the Period4 quasi-periodic orbit D match well with that of the reference solutions. Likewise, as seen in Fig. 26, the predicted chaotic attractor is confirmed by the numerical integration of the full system. Here, computing the maximum Lyapunov exponent of the chaotic attractor via a long-time forward integration of the variational equations for the 2046-dimensional full system (21) is impractical, and thus, we do not report the reference result for this maximum Lyapunov exponent.