1. Introduction
With the increasing number of passengers and the need for more freight transport around the world, logistics and transport companies are pressured to increase their number of flights. However, this has led to a major increase in global air traffic, and the aviation industry is under pressure to reduce noise and emissions. One of the potential solutions to this issue is to build airports away from populated areas in offshore locations. This means moving the take-off and landing paths over water. However, the cost of land reclamation and the need for new terminal buildings and pathways to be constructed are very high. An alternative to expensive infrastructure would be the use of a hybrid configuration that can reach high speeds while carrying a heavy payload. The introduction of high-speed planing hulls has led to the development of a hybrid configuration known as the seaplane, which is fast, agile, and capable of covering a wide range of applications because it combines the characteristics of aeroplanes and hydroplanes. The configuration of the seaplane is started as a ship, but with the advantage of using the aerodynamic phenomenon known as the ground effect to fly just above the surface of the water.
When a craft flies next to the surface of the water or ground, it is influenced by the surface effect of aero-hydrodynamics, which leads to a high lift-to-drag ratio. The deceleration of air trapped between the ground and the wing surface causes a significant increase in the pressure on the under-surface of the wing [
1]. Thus, the ground effect (GE) can be defined as the increase in the lift-to-drag ratio of a lifting surface flying at a small relative distance from the underlying surface [
2]. This technology was most notably used in the development of the “Ekranoplan” designed in the Soviet Union in the beginning of the 1960s [
1,
2,
3].
Seaplanes have a unique instability phenomenon called porpoising, which, according to Faltinsen [
4], can be defined as a periodic, bounded vertical motion that a craft might exhibit at certain speeds, as illustrated in
Figure 1. This phenomenon can be seen as an oscillatory motion in the heave and pitch axes and can cause severe damage to the structure of the craft [
5]. In some cases, if the hull leaves the water and returns with a negative trim angle, the craft will submarine [
6]. Porpoising can be fatal, even in calm water, because the loss of longitudinal stability can cause self-induced heave and pitch oscillations (porpoising) and the submergence of the bow area [
7].
It is of great importance to study the response of the coupled motions of an object when it is subjected to an external periodic force. Such a set-up finds applications in resonance vibrations, similar to what is seen when a seaplane is accelerating to take-off. This could lead to minimising unwanted resonance and avoiding chaos [
8]. In addition, the exact analytical solutions of the nonlinear, coupled, and forced equations of motion are very helpful for performing comprehensive analyses. These solutions can be used to capture the main nonlinear phenomena, such as nonlinear resonance and chaos, that characterise these systems with respect to the linear ones. In contrast to the numerical solutions, which can be useful for investigating a specific case, the analytical solution can be helpful for conducting a detailed parametric analysis aimed at understanding the whole system’s behaviour. This is needed in the pre-design of any mechanical structure that interacts with a moving fluid [
9]. Hence, understanding the role of coupling and nonlinearity on system behaviour could enhance the dynamic stability of seaplanes and provide insight into key design parameters. In fact, predicting the take-off performance by combining aerodynamic and hydrodynamic calculations is a major challenge in seaplane design [
10]. The structural design of a sea-craft requires the evaluation of its motions and wave-induced pressure distributions over the hull in wavy sea conditions. Thus, a nonlinear, time-dependent seakeeping solution must simulate, with sufficient accuracy and efficiency, the seakeeping behaviour of a sea-craft in wave trains that may be several hours long in duration [
11]. With this goal in mind, the development of a nonlinear periodic solution was initiated.
Most studies conducted to date on the dynamic stability of seaplanes have been experimental or numerical in nature. For instance, Guo et al. [
10] used the volume-of-fluid (VOF) approach to numerically simulate the air–water flow field to study the resistance characteristics and attitude of a seaplane at various take-off speeds. They concluded that a seaplane's trim angle should be operated in an ideal trim range, typically between 4° and 6°, to minimise hydrodynamic resistance. Song et al. [
12] conducted a numerical study to understand the characteristics of a seaplane’s take-off performance. The study only focused on the suction force, pressure, free surface, and motion of the seaplane model to understand the characteristics of the planing motion of the seaplane when moving through waves of different velocities, wavelengths, and heights. A two-dimensional numerical investigation was performed using the VOF approach proposed by Xie et al. [
13]. The latter used a gas–liquid two-phase model to simulate complex fluid fluxes, with the air–water interface captured using the VOF technique. Their findings elucidated the mechanism of complicated flows under the influence of the bottom cabin and inclined postures. After conducting a comprehensive review of the academic literature on seaplane hydrodynamics, Tavakoli et al. [
14] concluded that more work needs to be conducted on mathematical models, particularly those that simulate manoeuvring motions, to fully account for the challenges of operating in a real-world maritime environment. Zhou et al. [
15] stated that the majority of the present research on seaplane planing in waves is still concentrated on experiments and numerical simulations in calm water or under low-wave conditions. According to Wang et al. [
16], it is no longer difficult to forecast the porpoising conditions for a specific sea-craft, experimentally or numerically. The challenge is finding out how to solve the porpoising problem, which is achievable with an understanding of the nonlinear behaviour of sea-crafts accelerating on a wavy sea surface. Zheng et al. [
17] used a modified two-phase flow solver in OpenFOAM to simulate the porpoising phenomenon of a large seaplane with four turbines during take-off. They discovered that heaving and pitching oscillations are primarily determined by the hydrodynamic force, with the aerodynamic force playing an auxiliary role. A numerical study by Zha et al. [
18] investigated the water-landing behavior of seaplanes and concluded that careful control of seaplane motion is essential to prevent potential damage from overwashing events. After thoroughly reviewing the methods used in seaplane hydrodynamic design, Morabito [
19] stated that, because the majority of the research in this field dates back to the 1940s, studying the hydrostatic and hydrodynamic challenges associated with seaplane design is currently an important endeavour.
This study aimed to present an analytical tool that can be used to predict the nonlinear motion of seaplanes advancing through head-sea waves. The analytical tool is a system of coupled, nonlinear heaving and pitching equations of motion, which are assumed to be externally excited by a sinusoidal head-sea wave. These equations are presented in the form of a system of coupled Duffing equations with cubic nonlinearity. The focus is on heave and pitch motions because, at a speed-to-length ratio of less than 4 (take-off and landing speeds), the hydrodynamic coefficients associated with motions in other directions are neglected [
20]. Hence, the dynamic stability of seaplanes during take-off and landing was examined from heaving and pitching motions. The analytical tool was validated using computational fluid dynamics (CFD) simulations performed on Ansys Fluent, using the VOF method, as the flow was two-phase. Analytical data for aero-hydrodynamics and flight mechanics have a direct impact on safety and reliability and can assist in enhancing both. This provides an opportunity to investigate the nature of the nonlinearity derived directly from motion physics. This should provide a deeper understanding of the nonlinear phenomenon associated with motion across waves, by explaining the effect of the nonlinear and coupling coefficients on the magnitude of the external force/moment. This study is the first of its kind to use the coupled nonlinear Duffing equation to describe the heaving and pitching motions of seaplanes advancing through head-sea waves. A closed-form analytical solution to the equations of motion, describing the structural response of the seaplane hull, is also presented in this study.
2. Mathematical Formulation
Coupled Duffing oscillators were previously investigated analytically in [
8,
9,
21,
22,
23]. In all previous works, the attention has been focused on the route to chaos and the frequency response. No attempts have been made to obtain an exact, closed-form analytical solution representing the heaving and pitching responses of a seaplane advancing through head-sea waves, which was the aim of this work. The chaos theory is one of the most significant achievements in nonlinear science. The Duffing equation is associated with the mathematical chaotic behaviour. This chaotic behaviour exists in many natural systems, such as weather and climate systems [
24]. It also occurs spontaneously in systems with artificial components, such as road traffic. Chaos theory has been applied in several other disciplines, including meteorology, sociology, and environmental sciences. Chaos is defined as a periodic long-term behaviour in a deterministic system that exhibits a sensitive dependence on the initial conditions [
25]. Three properties must exist in a dynamical system for it to be classified as chaotic:
It must exhibit periodic long-term behaviour, meaning that the solution of the system settles into an irregular pattern as . The solution either does not repeat or oscillates in periodically.
It is sensitive to the initial conditions. This means that any small change in the initial condition can change the trajectory, which may result in significantly different long-term behaviour.
It must be “deterministic”, which means that the irregular behaviour of the system is due to the nonlinearity of the system, rather than external forces.
Thus, Duffing oscillators find applications in chaos theory, which is a field of study in mathematics that investigates the behaviour of dynamical systems that are highly sensitive to the initial conditions. Small differences in the initial conditions (such as rounding errors in the numerical computation) yield widely diverging outcomes for such dynamical systems, rendering long-term predictions generally impossible [
25]. There are several analytical methods that can be used to solve the Duffing equation in its various forms, such as the method of multiple scales [
26], the Bogolubov–Mitropolski method [
27], He’s energy balance method [
28], the global error minimisation method [
29], the variational iteration method [
30], the Jacobi elliptic functions method [
31], the homotopy perturbation method [
32], and the Poincaré–Lindstedt method [
33].
Several physical problems faced by engineers today demonstrate certain fundamental features, such as nonlinearities, that prevent exact analytical solutions. Thus, perturbation methods are used to obtain information regarding the solutions of nonlinear equations. The solution is presented in the form of convergent power series expansions with respect to a small dimensionless parameter that is of the same order of nonlinearity and amplitude of motion [
34]. The Poincaré–Lindstedt method is used to uniformly approximate periodic solutions to ordinary differential equations. This method eliminates the secular terms arising from the straightforward application of the regular perturbation theory to weakly nonlinear equations [
35]. Secular terms are those that grow without bounds. These terms have a singularity point at which a given mathematical object is not defined [
35].
The Poincaré–Lindstedt method was adopted in this research to find a closed-form solution to the coupled equations of motion because it is used to obtain periodic solutions to nonlinear differential equations. It describes the period of the unstable motion of seaplanes such that porpoising can be predicted analytically. The use of this technique leads to an expression for the desired solution in terms of a formal power series with a small parameter (
), known as a perturbation series, that quantifies the deviation from the exact solvable problem [
36,
37,
38].
2.1. Heave and Pitch Equations of Motion
The equations of motion for a seaplane advancing at a constant forward velocity with an arbitrary heading in regular sinusoidal sea waves are presented in this section. The performance of the seaplane is presented using a six-DOF (degrees of freedom) system. Considering that the responses are linear and harmonic, the six linear equations of motion can be written using the subscript notation as follows [
4]:
where
is the component of the generalised mass matrix of the craft in the direction owing to the motion.
is the added-mass coefficient in the direction due to the motion.
is the damping coefficient in the direction due to the motion.
is the hydrostatic restoring force coefficient in the direction due to the motion.
is the complex amplitude of the exciting forces and moments in the direction.
Heaving and pitching due to moving through head-sea waves are oscillatory motions along the vertical and sagittal axes, as illustrated in
Figure 2. Heaving can be defined as the translational (linear) motion in the vertical direction, while pitching is the rotational (angular) motion of the centre of gravity of the aircraft [
39].
The heaving and pitching system of seaplane motion behaves like a two-degrees-of-freedom spring-mass system. According to Faltinsen [
4], this assumption is clear when a craft model is given heave or pitch displacements from its equilibrium position; it rapidly oscillates several times before it comes to rest. Therefore, the resulting equations for the heave and pitch motions of the seaplane can be expressed as follows:
Subscripts 3 and 5 represent the heave and pitch motions, respectively. In addition,
is the heave translational displacement, and
is the pitch angular displacement.
is the external excitation force and
is the external excitation moment. Both quantities were assumed to be caused by sinusoidal waves. Hence, the excitation forces and moments can be expressed as follows:
Both the time-dependent sinusoidal force and moment functions oscillate at an external excitation frequency . Here, denotes time, represents the amplitude of the force, is the amplitude of the moment, and is the phase angle, indicating a phase shift in the moment relative to the force.
The hydrostatic restoring force coefficient (buoyancy), hydrodynamic coefficients (added mass and damping), and external forces and moments (amplitudes and frequency) acting on a planing hull were calculated using strip theory. This theory utilises the integration of the hydrodynamic pressure acting on the wetted surface of the planing hull to obtain the forces and moments. The calculation of the coefficients is discussed in detail in
Section 4.1.
The general equations of the heave and pitch motions of seaplanes were reduced to Duffing oscillator equations by assuming that the restoring forces behave according to Hooke’s law but have a nonlinear nature (cubic nonlinearity). The system is excited by a sinusoidal head-sea wave, and the damping effect is negligible because the hydrodynamic coefficients associated with damping are much smaller than the added mass and restoring force coefficients. This is because the damping coefficients are functions of frequency and speed, whereas the added mass and restoring force coefficients are functions of the seaplane mass [
40]. In addition, the pitch’s angular acceleration was neglected in the heave equation, and the heave translational acceleration was neglected in the pitch equation. Nevertheless, the elastic deformation of the rigid body was neglected because the system of equations of motion was formed using Lagrangian mechanics, in which the system is assumed to be conserved. Finally, a mathematical model was developed to predict porpoising when the craft is under maximum buoyancy support; hence, the hydrostatic restoring force is dominant. Taking this into consideration, this two-degrees-of-freedom system can be modelled using two coupled second-order nonlinear differential equations with the following forms:
where
is the pitch coupling term coefficient.
is the heave nonlinear term coefficient.
is the heave coupling term coefficient.
is the pitch nonlinear term coefficient.
is the heave amplitude.
and are the pitch amplitudes.
is the natural frequency.
is the external excitation frequency.
Equations (5) and (6) are harmonically excited Duffing oscillator equations. They can be solved using the Poincaré–Lindstedt method by introducing a new, small, non-dimensional parameter (
) into the equations. This was achieved by assuming that the nonlinearity, coupling, and excitation force/moment are small relative to the harmonic motion and that they appear in the same order in each equation. In addition, (
) can be used as a measure of the strength of the nonlinearity. Hence, a small parameter (
) was introduced, and Equations (5) and (6) can be rewritten as follows:
2.2. Amplitude Equations
We next turned to the question of determining an expression for the amplitude and frequency of porpoising oscillations in terms of the nonlinearity, coupling, and external forcing parameters. The frequency of the system
can be explicitly presented in the differential equation by introducing the following transformation:
The frequency of the forcing term was initially assumed to be unknown, but near a multiple of the natural frequency,
. To account for this,
was introduced to define the order of resonance, allowing the solution to cover various cases (i.e., different orders of resonance). By applying the chain rule, the equations were transformed as follows:
where the prime symbol indicates the derivative with respect to the re-scaled time variable
. An unknown frequency of the system appeared in the differential equations. Hence, approximate solutions were obtained for both equations in the form of power series expansions in terms of
. In other words, the following were introduced:
The first term in the expansion represents the linear solution of the equation, also referred to as the leading-order term. In this study, the leading-order frequency,
, was set to unity, representing the natural frequency of the system. Substituting Equations (12), (13), and (14) into Equations (10) and (11) yielded the following:
Separating the terms of the same order and eliminating the terms of order yielded the following equations:
The leading-order solutions of the two equations were assumed to have the following form:
The corrections to the linear solution were determined in the course of the analysis by requiring the expansion of
and
to be uniform for all
. However, the particular solution of the first-order terms of
and
contains secular terms that make the expansion non-uniform. To achieve uniform expansion, the secular terms in
,
,
, and
must be eliminated. To this end, the coefficients of the secular terms were set to zero to find expressions for the first-order frequency
. Thus, only the inhomogeneous terms in Equations (18), (19), (21), and (22) governing
,
,
, and
were inspected. This resulted in the expansions being uniform for all first- and second-order solutions because secular terms did not appear in them. Removing all secular terms up to the first order required satisfying the following system of equations:
while
satisfied the following cubic polynomial:
represents the amplitude of the external force, which is the amplitude of the sea wave. The nonlinear coefficients
and
were obtained from the polynomial, assuming in the first approximation that both coefficients are equal. The other terms were obtained from the following equations:
Hence, the amplitude equations obtained are as follows: