Mousai: An Open-Source General Purpose Harmonic Balance Solver

ASME Dayton Engineering Sciences Symposium 2017
Joseph C. Slater, October 23, 2017

Overview

A wide array of contemporary problems can be represented by nonlinear ordinary differential equations with solutions that can be represented by Fourier Series:

• Limit cycle oscillation of wings/blades
• Flapping motion of birds/insects/ornithopters
• Flagellum (threadlike cellular structures that enable bacteria etc. to swim)
• Shaft rotation, especially including rubbing or nonlinear bearing contacts
• Engines
• Radio/sonar/radar electronics
• Wireless power transmission
• Power converters
• Boat/ship motions and interactions

• Cardio systems (heart/arteries/veins)
• Ultrasonic systems transversing nonlinear media
• Responses of composite materials or materials with cracks
• Near buckling behavior of vibrating columns
• Nonlinearities in power systems
• Energy harvesting systems
• Wind turbines
• Radio Frequency Integrated Circuits
• Any system with nonlinear coatings/friction damping, air damping, etc.

These can all be observed in a quick literature search on 'Harmonic Balance'.

Why (did I) write Mousai?

• The ability to code harmonic balance seems to be publishable by itself
  • It's not research- it's just application of a known family of technique
• A limited number of people have this knowledge and skill
  • Most cannot access this technique
  • "Research effort" is spent coding the technique, not doing research

Why write Mousai? (continued)

• Matlab command eig unleashed power to the masses
  • Very few papers are published on eigensolutions- they have to be better than eig
  • eig only provides simple access to high-end eigensolvers written in C and Fortran
  • Undergraduates with no practical understanding of the algorithms easily solve problems that were intractable a few decades ago.
  • Access and ease of use of such techniques enable greater science and greater research.
  • The real world is nonlinear, but linear analysis dominates because the tools are easier to use.
  • With Mousai, an undergraduate can solve a nonlinear harmonic response problem easier then a PhD can today.

Theory:

Linear Solution

• Most dynamics systems can be modeled as a first order differential equation \begin{equation} \ddot{\mathbf{z}}(t)=\mathbf{f}(\mathbf{z}(t),\mathbf{u}(t)) \end{equation}
  • Use finite differences
  • Use Galerkin methods (Finite Elements)
  • Of course discrete objects

• This is the common State-Space form: -solutions exceedingly well known if it is linear

• Finding the oscillatory response, after dissipation of the transient response, requires long time marching.

  • Without damping, this may not even been feasible.
  • With damping, tens, hundreds, or thousands of cycles, therefore thousands of times steps at minimum.

For a linear system in the frequency domain this is \begin{equation} j\omega\mathbf{Z}(\omega)=\mathbf{f}(\mathbf{Z}(\omega),\mathbf{U}(\omega)) \end{equation} \begin{equation} j\omega\mathbf{Z}(\omega)=A\mathbf{Z}(\omega)+B\mathbf{U}(\omega) \end{equation} where \begin{equation} A = \frac{\partial \mathbf{f}(\mathbf{Z}(\omega),\mathbf{U}(\omega))}{\partial\mathbf{Z}(\omega)},\qquad B = \frac{\partial \mathbf{f}(\mathbf{Z}(\omega),\mathbf{U}(\omega))}{\partial\mathbf{U}(\omega)} \end{equation} are constant matrices.

The solution is: \begin{equation} \mathbf{Z}(\omega) = \left(Ij\omega-A\right)^{-1}B\mathbf{U}(\omega) \end{equation} Where the magnitudes and phases of the elements of $\mathbf{Z}$ provide the amplitudes and phases of the harmonic response of each state at the frequency $\omega$.

Nonlinear solution

• For a nonlinear system in the frequency domain we assume a Fourier series solution \begin{equation} \mathbf{z}(t)=\lim_{N\to\infty}\sum_{n=-N}^{N}\mathbf{Z}_n e^{j n \omega t} \end{equation}

• $N=1$ for a single harmonic. $n=0$ is the constant term.

• This can be substituted into the governing equation to find $\dot{\mathbf{z}}(t)$: \begin{equation} \dot{\mathbf{z}}(t)=\mathbf{f}(\mathbf{z}(t),\mathbf{u}(t)) \end{equation}

• This is actually a function call to a Finite Element Package, CFD, Matlab function, - whatever your solver uses to get derivatives

• We can also find $\dot{\mathbf{z}}(t)$ from the derivative of the Fourier Series: \begin{equation} \dot{\mathbf{z}}(t)=\lim_{N\to\infty}\sum_{n=-N}^{N}j n \omega\mathbf{Z}_n e^{j n \omega t} \end{equation}

• The difference between these methods is zero when $\mathbf{Z}_n$ are correct.

\begin{equation} \mathbf{0} \approx\sum_{n=-N}^{N}j n\omega \mathbf{Z}_n e^{j n \omega t}-\mathbf{f}\left(\sum_{n=-N}^{N}\mathbf{Z}_n e^{j n \omega t},\mathbf{u}(t)\right) \end{equation}

• These operations are wrapped inside a function that returns this error
• This function is used by a Newton-Krylov nonlinear algebraic solver.
• Calls any solver in the SciPy family of solvers with the ability to easily pass through parameters to the solver and to the external derivative evaluator.

Examples:

Duffing Oscillator

$$\ddot{x}+0.1\dot{x}+x+0.1 x^3=\sin(\omega t)$$

# Define our function (Python)
def duff_osc_ss(x, params):
    omega = params['omega']
    t = params['cur_time']
    xd = np.array([[x[1]],
                   [-x[0] - 0.1 * x[0]**3 - 0.1 * x[1] + 1 * sin(omega * t)]])
    return xd

# Arguments are name of derivative function, number of states, driving frequency,
# form of the equation, and number of harmonics

t, x, e, amps, phases = ms.hb_so(duff_osc_ss, num_variables=2, omega=.1,
                                 eqform='first_order', num_harmonics=5)
print('Displacement amplitude is ', amps[0])
print('Velocity amplitude is ', amps[1])

Displacement amplitude is  0.946956354668
Velocity amplitude is  0.0946956354558

Mousai can easily recreate the near-continuous response

time, xc = ms.time_history(t, x)

pltcont()# abbreviated plotting function

Let's sweep through driving frequencies to find a frequency response function

omegal = np.arange(3, .03, -1 / 200) + 1 / 200
ampl = sp.zeros_like(omegal)
ampl[:] = np.nan
t, x, e, amps, phases = ms.hb_so(duff_osc_ss, num_variables=2,
                                 omega=3, eqform='first_order', num_harmonics=1)
for i, freq in enumerate(omegal):
    # Here we try to obtain solutions, but if they don't work,
    # we ignore them by inserting `np.nan` values.
    x = x - np.average(x)
    try:
        t, x, e, amps, phases =
        ms.hb_so(duff_osc_ss, x0=x,
                 omega=freq, eqform='first_order', num_harmonics=1)
        ampl[i] = amps[0]
    except:
        ampl[i] = np.nan
    if np.isnan(ampl[i]):
        break

plt.plot(omega,amp, label='Up sweep')
plt.plot(omegal,ampl, label='Down sweep')
plt.legend()
plt.title('Amplitude versus frequency for Duffing Oscillator')
plt.xlabel('Driving frequency $\\omega$')
plt.ylabel('Amplitude')
plt.grid()

Two degree of freedom system

$ \begin{bmatrix} 1&0\\ 0&1 \end{bmatrix} \begin{bmatrix} \ddot{x}_1\\ \ddot{x}_2\end{bmatrix} + \begin{bmatrix} 2&-1 \\ -1&2 \end{bmatrix} \begin{bmatrix} {x}_1\\{x}_2\end{bmatrix} +$ $ \begin{bmatrix} \alpha x_{1}^{3}\\ 0 \end{bmatrix} $ $= \begin{bmatrix} 0 \\ A \sin(\omega t) \end{bmatrix} $

def two_dof_demo(x, params):
    omega = params['omega']
    t = params['cur_time']
    force_amplitude = params['force_amplitude']
    alpha = params['alpha']
    # The following could call an external code to obtain the state derivatives
    xd = np.array([[x[1]],
                   [-2 * x[0] - alpha * x[0]**3 + x[2]],
                   [x[3]],
                   [-2 * x[2] + x[0]]] + force_amplitude * np.sin(omega * t))
    return xd

Let's find a response.

parameters = {'force_amplitude': 0.2}
parameters['alpha'] = 0.4
t, x, e, amps, phases = ms.hb_so(two_dof_demo, num_variables=4,
                                 omega=1.2, eqform='first_order', params=parameters)
amps

array([ 0.86696762,  0.89484597,  0.99030411,  1.04097851])

Or a parametric study of response amplitude versus nonlinearity.

alpha = np.linspace(-1, .45, 2000)
amp = np.zeros_like(alpha)
for i, alphai in enumerate(alpha):
    parameters['alpha'] = alphai
    t, x, e, amps, phases = ms.hb_so(two_dof_demo, num_variables=4, omega=1.2,
                                     eqform='first_order', params=parameters)
    amp[i] = amps[0]

plt.plot(alpha,amp)
plt.title('Amplitude of $x_1$ versus $\\alpha$')
plt.ylabel('Amplitude of $x_1$')
plt.xlabel('$\\alpha$')
plt.grid()

Two degree of freedom system with Coulomb Damping

$ \begin{bmatrix} 1&0\\ 0&1 \end{bmatrix} \begin{bmatrix} \ddot{x}_1\\ \ddot{x}_2\end{bmatrix} + \begin{bmatrix} 2&-1 \\ -1&2 \end{bmatrix} \begin{bmatrix} {x}_1\\{x}_2\end{bmatrix} +$ $ \begin{bmatrix} \mu |\dot{x}|_{1}\\ 0 \end{bmatrix} $ $= \begin{bmatrix} 0 \\ A \sin(\omega t) \end{bmatrix} $

def two_dof_coulomb(x, params):
    omega = params['omega']
    t = params['cur_time']
    force_amplitude = params['force_amplitude']
    mu = params['mu']
    # The following could call an external code to obtain the state derivatives
    xd = np.array([[x[1]],
                   [-2 * x[0] - mu * np.abs(x[1]) + x[2]],
                   [x[3]],
                   [-2 * x[2] + x[0]]] + force_amplitude * np.sin(omega * t))
    return xd

mu = np.linspace(0, 1.0, 200)
amp = np.zeros_like(mu)
for i, mui in enumerate(mu):
    parameters['mu'] = mui
    t, x, e, amps, phases = ms.hb_so(two_dof_coulomb, num_variables=4, omega=1.2,
                                     eqform='first_order', num_harmonics=3, params=parameters)
    amp[i] = amps[0]

Too much Coulomb friction can increase the response.

• Did you know that?
• This damping shifted resonance.

plt.plot(mu,amp)
plt.title('Amplitude of $x_1$ versus $\\mu$')
plt.ylabel('Amplitude of $x_1$')
plt.xlabel('$\\mu$')
plt.grid()

But can I solve an equation in one line? Yes!!!

Damped Duffing oscillator in one command.

out = ms.hb_so(lambda x, v,
               params: np.array([[-x - .1 * x**3 - .1 * v + 1 *
                                  sin(params['omega'] * params['cur_time'])]]),
               num_variables=1, omega=.7, num_harmonics=1)
out[3][0]

1.4779630014433971

OK - that's a bit obtuse. I wouldn't do that normally, but Mousai can.

How to get this?

• Install Scientific Python from SciPy.org
  • AFRL: See your tech support to get the Enthought distribution installed
• See the Mousai documents for installation instructions
  • pip install mousai
  • AFRL: Talk to me- install is easy if I send you the files.
• See Mousai on GitHub (https://github.com/josephcslater/mousai)

Conclusions

• Nonlinear frequency solutions are within reach of undergraduates
• Installation is trivial
• Already in use (GitHub logs indicate dozens of users)
• Custom special case and proprietary solvers such as BDamper can be replaced for free
• Research potential is about to be unleashed

Future

• Add time-averaging method
  • currently requires high number of harmonics for non-smooth systems
• Add masking of known harmonics (average is often fixed and known)
• Automated sweep control
• Branch following
• Condense the one-line method
• Evaluate on large scale problems
  • Currently attempting to hook to ANSYS
• Parallelize
• Leverage CUDA