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 [2], 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_oscillator(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# %%