Example 10: Subharmonic response of parametrically excited coupled centrifugal pendulums in 1:1 internal resonance

MMS example on two coupled centrifugal pendulums in 1:1 internal resonance subject to direct and parametric forcings, triggering a subharmonic response. The pendulums are indirectly and weakly coupled through the rotor that supports them. This configuration was studied in V. Mahé et al., Subharmonic centrifugal pendulum vibration absorbers allowing a rotational mobility, MSSP (2022).

System description

Nonlinear system.

Illustration of the modes of a CPVA made of a rotor and 2 pendulums. The rotor is in grey, the pendulums in black, they path in blue, and their position on that path in red. Their rotation with respect to the rotor is not represented here. The direction of the rotor’s motino with respect to the pendulums’ is shown through the black arrow.

Nonlinear system.

Illustration of CPVA modes 1 (phase opposition) and 2 (unison) as oscillators. This representation can be adopted when solving for the rigid body mode, and then injecting this solution in the two other modes. This reduces the system’s degrees of freedom by one, and produces the 2 coupled (reduced) modes system. Note that the excitation order \(n\) and the rotor’s angular position \(\vartheta\) in this sketch stand for the excitation frequency \(\omega\) and time \(t\), as is usual in the context of rotating machines.

The system’s equations in modal space are

\[\begin{split}\begin{cases} \ddot{\zeta}_1 + n_p^2 \zeta_1 & = f_1(\boldsymbol{\zeta}, \dot{\boldsymbol{\zeta}}, \ddot{\boldsymbol{\zeta}}, t) \\ \ddot{\zeta}_2 + n_p^2 \zeta_2 & = f_2(\boldsymbol{\zeta}, \dot{\boldsymbol{\zeta}}, \ddot{\boldsymbol{\zeta}}, t), \end{cases}\end{split}\]

where

  • \(\zeta_1,\; \zeta_2\) are the (modal) oscillators’ coordinates,

  • \(\boldsymbol{\zeta}\) is the vector of oscillator coordinates,

  • \(t\) is the angular position of the rotor supporting the pendulums, analogous to time

  • \(\dot{(\bullet)} = \mathrm{d}(\bullet)/\mathrm{d}t\) is a angle derivative,

  • \(n_p\) is the pendulums’ tuning order and their natural frequency when they are uncoupled from the rotor,

  • \(f_1(\boldsymbol{\zeta}, \dot{\boldsymbol{\zeta}}, \ddot{\boldsymbol{\zeta}}, t)\) and \(f_2(\boldsymbol{\zeta}, \dot{\boldsymbol{\zeta}}, \ddot{\boldsymbol{\zeta}}, t)\) are (small) functions involving linear damping, an additional linear stiffness on oscillator 2, nonlinear terms, and both direct and parametric forcings. They take the form:

    \[\begin{split}\begin{aligned} \begin{split} f_1(\boldsymbol{\zeta}, \dot{\boldsymbol{\zeta}}, \ddot{\boldsymbol{\zeta}}, t) = & -\Lambda_m^{-1} \Big[ (\Lambda_m \dot{\zeta}_1 - n_t^2(1+n_t^2)\zeta_1\zeta_2)(\mu n_p^2 \Lambda_c \zeta_2 + T_1\cos(\omega t)) \\ & + 2 \mu n_t^2 \Lambda_m \dot{\zeta}_1 (\zeta_1\dot{\zeta}_1 + \zeta_2 \dot{\zeta}_2) \\ & + 6 \eta \alpha_1 \alpha_3 (\zeta_1 \dot{\zeta}_1^2 + 2 \zeta_2 \dot{\zeta}_1 \dot{\zeta}_2 + \zeta_1 \dot{\zeta}_2^2 + \zeta_1^2 \ddot{\zeta}_1 + 2 \zeta_1 \zeta_2 \ddot{\zeta}_2 + \zeta_2^2 \ddot{\zeta}_1) \\ & - 2 x_4 (\zeta_1^3 + 3\zeta_1\zeta_2^2) + b \dot{\zeta}_1 \Big], \end{split} \\[ 0.5em] \begin{split} f_2(\boldsymbol{\zeta}, \dot{\boldsymbol{\zeta}}, \ddot{\boldsymbol{\zeta}}, t) = & -\Lambda_m^{-1} \Bigg[ \left(\Lambda_c + \Lambda_m \dot{\zeta}_2 - \frac{n_t^2(1+n_t^2)}{2}(\zeta_1^2 + \zeta_2^2)\right)(\mu n_p^2 \Lambda_c \zeta_2 + T_1\cos(\omega t)) \\ & + 2 \mu n_t^2 (\Lambda_c + \Lambda_m \dot{\zeta}_2) (\zeta_1 \dot{\zeta}_1 + \zeta_2 \dot{\zeta}_2) \\ & - \Lambda_c \mu n_p^2\frac{n_t^2(1+n_t^2)}{2}(3\zeta_1^2\zeta_2 + \zeta_2^3) \\ & + \Lambda_c \mu n_t^2(1+n_t^2)(2\zeta_1 \dot{\zeta}_1 \dot{\zeta}_2 + \zeta_2\dot{\zeta}_1^2 + \zeta_2 \dot{\zeta}_2^2) \\ & + 6 \eta \alpha_1 \alpha_3 (2\zeta_1 \dot{\zeta}_1 \dot{\zeta}_2 + \zeta_2 \dot{\zeta}_1^2 + \zeta_2 \dot{\zeta}_2^2 + 2\zeta_1 \ddot{\zeta}_1 \zeta_2 + \ddot{\zeta}_2 \zeta_1^2 + \ddot{\zeta}_2 \zeta_2^2) \\ & - 2\tilde{x}_{4}(3\zeta_1^2\zeta_2 + \zeta_2^3) + \tilde{b}\dot{\zeta}_2 \Bigg]. \end{split} \end{aligned}\end{split}\]

    where

    • \(\Lambda_m\) is an equivalent pendulum mass,

    • \(\Lambda_c\) is a linear coupling coefficient among the pendulums,

    • \(\mu\) is a linear coupling coefficient between the pendulums and the rotor that supports them,

    • \(n_t\) and \(x_4\) are linear and nonlinear pendulums’ path coefficients, respectively,

    • \(\alpha_1\) and \(\alpha_3\) are linear and nonlinear pendulums’ rotation coefficients, respectively,

    • \(b\) is the linear viscous damping coefficient,

    • \(T_1\) is the forcing amplitude,

    • \(\omega\) is the forcing frequency (often written \(n\) as it is actually a forcing order).

A response around \(n_p\) is sought so the frequency is set to

\[\omega = n_p + \epsilon \sigma\]

where

  • \(\epsilon\) is a small parameter involved in the MMS,

  • \(\sigma\) is the detuning wrt the pendulums’ tuning order.

The parameters are then scaled to indicate how weak they are:

  • \(\mu = \epsilon \tilde{\mu}\) indicates that pendulums are weakly coupled,

  • \(b = \epsilon \tilde{b}\) indicates that damping is weak,

  • \(T_1 = \epsilon \tilde{T}_1\) indicates that forcing is weak,

  • \(x_4 = \epsilon \tilde{x}_4, \; \alpha_3 = \tilde{\alpha}_3\) indicates that nonlinearities are weak.

Code description

The script below allows to

  • Construct the dynamical (modal) system,

  • Apply the MMS to the system,

  • Evaluate the MMS results at steady state,

  • Introduce compact notations,

  • Compute the backbone curve and forced response when only oscillator (mode) 1 responds,

  • Evaluate the stability of the forced solution using cartesian coordinates,

  1# -*- coding: utf-8 -*-
  2
  3#%% Imports and initialisation
  4from sympy import symbols, Function, Rational, sin, cos, sqrt, solve
  5from sympy.physics.vector.printing import init_vprinting
  6init_vprinting(use_latex=True, forecolor='White') # Initialise latex printing 
  7from oscilate import MMS
  8
  9# Parameters and variables
 10nP, nt, mu, eta, Lm, b, T1  = symbols(r"n_p, n_t, \mu, \eta, \Lambda_m, b, T_1", real=True, positive=True)
 11Lc, alpha1, alpha3, x4 = symbols(r"\Lambda_c, \alpha_1, \alpha_3, x_4", real=True)
 12t = symbols('t') # Replaces the rotor's angular position \vartheta.
 13z1 = Function(r'\zeta_1', real=True)(t) # Mode 1 coordinate (phase-opposition mode)
 14z2 = Function(r'\zeta_2', real=True)(t) # Mode 2 coordinate (unison mode)
 15
 16# Dynamical system
 17dz1 = z1.diff(t)
 18dz2 = z2.diff(t)
 19ddz1 = z1.diff(t,2)
 20ddz2 = z2.diff(t,2)
 21
 22f1 = -Lm**-1 * \
 23     (
 24         (Lm*dz1 - nt**2*(1+nt**2)*z1*z2) * (mu*nP**2*Lc*z2)
 25         + 2*mu*nt**2*Lm*dz1*(z1*dz1 + z2*dz2)
 26         + 6*eta*alpha1*alpha3*(z1*dz1**2 + 2*z2*dz1*dz2 + z1*dz2**2 + z1**2*ddz1
 27                                + 2*z1*z2*ddz2 + z2**2*ddz1)
 28         - 2*x4*(z1**3 + 3*z1*z2**2)
 29         + b*dz1
 30       )
 31
 32f2 = -Lm**-1 * \
 33     (
 34         (Lc + Lm*dz2 - Rational(1,2)*nt**2*(1+nt**2)*(z1**2 + z2**2)) * (mu*nP**2*Lc*z2)
 35         + 2*mu*nt**2*(Lc + Lm*dz2)*(z1*dz1 + z2*dz2)
 36         - Lc*mu*nP**2*Rational(1,2)*nt**2*(1+nt**2)*(3*z1**2*z2 + z2**3)
 37         + Lc*mu*nt**2*(1+nt**2)*(2*z1*dz1*dz2 + z2*dz1**2 + z2*dz2**2)
 38         + 6*eta*alpha1*alpha3*(2*z1*dz1*dz2 + z2*dz1**2 + z2*dz2**2 
 39                                + 2*z1*ddz1*z2 + ddz2*z1**2 + ddz2*z2**2)
 40         - 2*x4*(3*z1**2 *z2 + z2**3)
 41         + b*dz2
 42       )
 43
 44fF1 = -Lm**-1 * (Lm*dz1 - nt**2*(1+nt**2)*z1*z2)
 45fF2 = -Lm**-1 * (Lc + Lm*dz2 - Rational(1,2)*nt**2*(1+nt**2)*(z1**2 + z2**2))
 46
 47Eq1 = ddz1 + nP**2*z1 - f1
 48Eq2 = ddz2 + nP**2*z2 - f2
 49
 50kwargs_dyn = dict(F=T1, fF=[fF1, fF2])
 51dyn = MMS.Dynamical_system(t, [z1, z2], [Eq1, Eq2], [nP, nP], **kwargs_dyn)
 52
 53# Initialisation of the MMS sytem
 54eps            = symbols(r"\epsilon", real=True, positive=True) # Small parameter epsilon
 55Ne             = 1  # Order of the expansions
 56omega_ref      = nP # Reference frequency
 57ratio_omegaMMS = 2
 58
 59param_to_scale = (mu, alpha3, x4, b, T1)
 60scaling        = (1 , 1,      1,  1, 1)
 61param_scaled, sub_scaling = MMS.scale_parameters(param_to_scale, scaling, eps)
 62
 63mms = MMS.Multiple_scales_system(dyn, eps, Ne, omega_ref, sub_scaling,
 64                                 ratio_omegaMMS=ratio_omegaMMS)
 65
 66# Application of the MMS
 67mms.apply_MMS()
 68
 69# Evaluation at steady state
 70ss = MMS.Steady_state(mms)
 71
 72# Substitute parameters
 73cc = symbols(r"c_c", real=True, positive=True)
 74ct, cp = symbols(r"c_t, c_p", real=True)
 75mut, alpha3t, x4t = param_scaled[0:3]
 76def_param = [(cp, 3*(x4t + 2*nP**2*eta*alpha1*alpha3t)),
 77             (ct, Rational(1,2) * Lc*mut*nP**2*nt**2*(1+nt**2)),
 78             (cc, mut*nt**4)
 79             ]
 80sub_nt = [(nt, nP*sqrt(Lm))]
 81sub_param = [(x4t, solve(def_param[0][1] - def_param[0][0], x4t)[0]),
 82             (Lc*mut*nt, solve(def_param[1][1] - def_param[1][0], Lc*mut*nt)[0]),
 83             (mut*nt, solve(def_param[2][1] - def_param[2][0], mut*nt)[0])]
 84
 85a, beta = ss.coord.a, ss.coord.beta
 86collect_f = [a[0]*a[1]**2*sin(beta[0]-beta[1]),
 87             a[0]**2*a[1]*sin(beta[0]-beta[1]),
 88             a[0]*a[1]**2*cos(beta[0]-beta[1]),
 89             a[0]**2*a[1]*cos(beta[0]-beta[1]),
 90             a[0]**3, a[1]**3, a[0]*a[1]**2, a[0]**2*a[1]**2]
 91ss.sol.fa    = [fa_ix.subs(sub_param+sub_nt).simplify().expand().collect(collect_f) for fa_ix in ss.sol.fa]
 92ss.sol.fbeta = [fbeta_ix.subs(sub_param+sub_nt).simplify().expand().collect(collect_f) for fbeta_ix in ss.sol.fbeta]
 93
 94# Solve the evolution equations for a given dof
 95solve_dof = 0 # dof to solve for
 96ss.solve_forced(solve_dof=solve_dof)
 97ss.solve_bbc(solve_dof=solve_dof, c=param_scaled[-2])
 98
 99# Stability analysis 
100kwargs_stab = dict(coord="cartesian", eigenvalues=True, bifurcation_curves=True, analyse_blocks=True, kwargs_bif=dict(var_a=True, var_sig=True, solver=solve))
101ss.stability_analysis(**kwargs_stab)
102
103# %%