Example 2: Coupled Duffings in 1:3 internal resonance

MMS example on coupled Duffing oscillators in 1:3 internal resonance subject to harmonic forcing. This configuration was studied by Nayfeh and Mook in Nonlinear Oscillations (1995), sections 6.6 and 6.6.2.

System description

Nonlinear system.

Illustration of two forced, nonlinearly coupled Duffing oscillators in 1:3 internal resonance, provided \(\omega_1 \approx 3 \omega_0\).

The system’s equations are

\[\begin{split}\begin{cases} \ddot{x}_{0} + 2 \mu_{0} \dot{x}_{0} + \omega_{0}^{2} x_{0} + \alpha_{1} x_{0}^{3} + \alpha_{2} x_{0}^{2} x_{1} + \alpha_{3} x_{0} x_{1}^{2} + \alpha_{4} x_{1}^{3} & = \Gamma_0 F \cos(\omega t), \\ \ddot{x}_{1} + 2 \mu_{1} \dot{x}_{1} + \omega_{1}^{2} x_{1} + \alpha_{5} x_{0}^{3} + \alpha_{6} x_{0}^{2} x_{1} + \alpha_{7} x_{0} x_{1}^{2} + \alpha_{8} x_{1}^{3} & = \Gamma_1 F \cos(\omega t), \end{cases}\end{split}\]

where

  • \(x_0,\; x_1\) are the oscillators’ coordinates,

  • \(t\) is the time,

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

  • \(\mu_0,\; \mu_1\) are the linear viscous damping coefficients,

  • \(\omega_0,\; \omega_1\) are the oscillator’s natural frequencies,

  • \(\alpha_i,\; i\in\{1,\cdots,8\}\) are nonlinear coefficients,

  • \(\Gamma_0,\; \Gamma_1\) are forcing coefficients,

  • \(F\) is the forcing amplitude,

  • \(\omega\) is the forcing frequency.

The oscillators are in 1:3 internal resonance. Taking \(\omega_0\) as the reference frequency, this implies

\[\omega_1 = 3\omega_0 + \sigma_1\]

where \(\sigma_1\) is the detuning of oscillator 1 wrt the 1:3 internal resonance condition.

A response around \(3\omega_0 \approx \omega_1\) is sought so the frequency is set to

\[\omega = 3\omega_0 + \epsilon \sigma\]

where

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

  • \(\sigma\) is the detuning wrt the frequency \(3\omega_0\).

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

  • \(\mu_i = \epsilon \tilde{\mu}_i\) indicates that damping is weak,

  • \(F = \epsilon \tilde{F}\) indicates that forcing is weak,

  • \(\sigma_1 = \epsilon \tilde{\sigma}_1\) indicates that detuning is small,

  • \(\alpha_i = \epsilon \tilde{\alpha}_i\) indicates that nonlinearities are weak.

Code description

The script below allows to

  • Construct the dynamical system.

  • Apply the MMS to the system,

  • Evaluate the MMS results at steady state,

  • Compute the backbone curve when only oscillator 1 responds.

 1# -*- coding: utf-8 -*-
 2
 3#%% Imports and initialisation
 4from sympy import symbols, Function
 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
10omega0, omega1, F, mu0, ma0 = symbols(r'\omega_0, \omega_1, F, \mu_0, \mu_1', real=True, positive=True)
11alphas = [symbols(r'\alpha_{{{}}}'.format(ii), real=True) for ii in range(1,9)] # Nonlinear coefficients
12Gam0, Gam1 = symbols(r"\Gamma_0, \Gamma_1", real=True) # Forcing coefficients
13
14t  = symbols('t')
15x0 = Function(r'x_0', real=True)(t)
16x1 = Function(r'x_1', real=True)(t)
17
18# Dynamical system
19Eq0 = x0.diff(t,2) + omega0**2*x0 + 2*mu0*x0.diff(t) + sum([alphas[ii]   * x0**(3-ii)*x1**ii for ii in range(4)])
20Eq1 = x1.diff(t,2) + omega1**2*x1 + 2*ma0*x1.diff(t) + sum([alphas[4+ii] * x0**(3-ii)*x1**ii for ii in range(4)])
21fF  = [Gam0, Gam1]
22dyn = MMS.Dynamical_system(t, [x0, x1], [Eq0, Eq1], [omega0,omega1], fF=fF, F=F)
23
24# Initialisation of the MMS sytem
25eps            = symbols(r"\epsilon", real=True, positive=True) # Small parameter epsilon
26Ne             = 1      # Order of the expansions
27omega_ref      = omega0 # Reference frequency
28ratio_omegaMMS = 3      # Look for a solution around omega1
29
30ratio_omega_osc = [1, 3] # Ratio between the oscillators' frequencies and the reference frequency
31sigma1          = symbols(r"\sigma_1", real=True) # Detuning of oscillator 1
32detunings       = [0, sigma1] # Detuning of the oscillators' frequency
33
34param_to_scale = (*alphas,          F, mu0, ma0, sigma1)
35scaling        = (*[1]*len(alphas), 1, 1,   1,   1)
36param_scaled, sub_scaling = MMS.scale_parameters(param_to_scale, scaling, eps)
37
38kwargs_mms = dict(ratio_omegaMMS=ratio_omegaMMS, ratio_omega_osc=ratio_omega_osc, detunings=detunings)
39mms = MMS.Multiple_scales_system(dyn, eps, Ne, omega_ref, sub_scaling, **kwargs_mms)
40
41# Application of the MMS
42mms.apply_MMS()
43
44# Evaluation at steady state
45ss = MMS.Steady_state(mms)
46ss.solve_bbc(solve_dof=1, c=param_scaled[9:11])
47
48# %%