Resonant three-dimensional nonlinear sloshing in a square base basin. Part 5. Three-dimensional non-parametric tank forcing

Assuming an inviscid incompressible liquid (with irrotational flows) partly filling a square base tank, which performs a small-amplitude sway/surge/pitch/roll periodic motion whose frequency is close to the lowest natural sloshing frequency, a nine-dimensional Narimanov–Moiseev-type (modal) system of ordinary differential equations with respect to the hydrodynamic generalised coordinates was derived in the Part 1 (Faltinsen et al., J. Fluid Mech., vol. 487, 2003, pp. 1–42). Constructing and analysing asymptotic periodic solutions of the system made it possible to classify steady-state resonant sloshing and its stability for the harmonic reciprocating (longitudinal, diagonal and oblique) forcing. The results were supported by experimental observations and measurements. The present paper finalises the case studies by considering the three-dimensional non-parametric (combined sway, pitch, surge, roll and yaw, but no heave) cyclic tank motions. It becomes possible after establishing an asymptotic equivalence of the associated periodic solutions of the modal system to those for a suitable horizontal translatory elliptic forcing so that, as a consequence, resonant steady-state waves and their stability can be considered versus angular position, semi-axis ratio $|\unicode[STIX]{x1D6FF}_{1}|$ and direction (counter- or clockwise) of the equivalent orbits. The circular orbit causes stable swirling waves (co-directed with the orbit) but may also excite stable nearly standing waves. The orbit direction does not affect the response curves for wall-symmetric (canonic) and diagonal orbit positions. This is not true for the oblique-type elliptic forcing. When the semi-axis ratio $|\unicode[STIX]{x1D6FF}_{1}|$ changes from 0 to 1, the response curves exhibit astonishing metamorphoses significantly influencing the frequency ranges of stable nearly standing/swirling waves and ‘irregular’ sloshing. For the experimental input data by Ikeda et al. (J. Fluid Mech., vol. 700, 2012, pp. 304–328), the counter-directed swirling disappears as $0.5\lesssim |\unicode[STIX]{x1D6FF}_{1}|$ but the frequency range of irregular waves vanishes for $0.75\lesssim |\unicode[STIX]{x1D6FF}_{1}|$.


Introduction
Suggesting an inviscid liquid with irrotational flows, a finite liquid depth, a forcing frequency σ close to the lowest natural sloshing frequency σ 1 , a forcing amplitude asymptotically smaller than the tank width/breadth (the five non-dimensional sway/surge/roll/pitch/yaw amplitudes are O( ) 1) and assuming that secondary resonance phenomena can be neglected, Faltinsen, Rognebakke & Timokha (2003, Part 1) derived a Narimanov-Moiseev-type nonlinear modal system (of ordinary differential equations), which effectively approximates resonant sloshing in a square base tank. The system couples the (hydrodynamic) generalised coordinates governing the free-surface elevation with nine degrees of freedom. Two generalised coordinates have the lowest asymptotic order O( 1/3 ), they are responsible for the primary excited natural sloshing modes. Three generalised coordinates have the second asymptotic order, O( 2/3 ), but the other four have the third asymptotic order O( ). This number of hydrodynamic generalised coordinates and their asymptotic order on the O( )-scale are a consequence of classical mathematical results by Moiseev (1958) and Narimanov (1957). Part 1 explains why the remaining (infinite number) degrees of freedom of the free surface are of the higher asymptotic order and can be neglected within the framework of the Narimanov-Moiseev theory. Under the assumptions stated above, the derived nonlinear modal system should theoretically be applicable to arbitrary (almost) periodic sway/surge/pitch/roll/yaw tank motions. However, being driven by available model tests and because of analytical problems (partly solved in the present paper), the forthcoming parts are exclusively focusing on the harmonic (longitudinal, oblique and diagonal) resonant tank reciprocation. The theoretical results on steady-state waves and their stability were supported by existing (including the authors') experimental observations and measurements. Raynovskyy & Timokha (2018a) employed an analogous Narimanov-Moiseev-type asymptotic approximation of the steady-state resonant sloshing for a circular base tank forced by sway, roll, surge, pitch and yaw. They showed that each resonant steady-state wave is then asymptotically equivalent to a steady-state wave excited by a suitable orbital elliptic horizontal translatory tank motion. The original and asymptotically equivalent solutions of the Narimanov-Moiseev-type modal equations differ only by the highest-order asymptotic components, which do not affect the waves stability. As a matter of the fact, studying the resonant steady-state sloshing due to periodic sway, surge, pitch, roll and yaw can be reduced to analysing the steady-state wave modes and their stability for tanks performing either reciprocating (a particular case; surveys are given by Royon-Lebeaud, Hopfinger & Cartellier (2007), Faltinsen, Lukovsky & Timokha (2016)) or orbital (relevant to bioreactors; see reviews by Weheliye, Yianneskis & Ducci (2013), Reclari et al. (2014), Bouvard, Herreman & Moisy (2017), Raynovskyy & Timokha (2018b) and ) tank motions.
Sloshing in a square base basin exposed to combined periodic sway, surge, roll, pitch and yaw is frequently associated with marine applications, storage and operational containers of, e.g. nuclear plants. Even though there exist numerous publications devoted to applied mathematical and engineering aspects of the problem (Ardakani & Bridges 2011;Ardakani 2019;Lyu et al. 2019;Zhang et al. 2019), studies on the resonant sloshing due to three-dimensional excitations are rare (Hiramitsu & Funakoshi 2015;Zhang, Ning & Teng 2018). There is a lack of fundamental knowledge on what types of resonant (steady-state) waves are realisable for this kind of three-dimensional forcing and in which frequency ranges. The Narimanov-Moiseev-type modal system from Part 1 allows us to fill this gap, A rigid square base container moves cyclically with a small amplitude in surge, sway, roll, pitch and yaw so that its translatory velocity is v O (t) = (v O1 (t), v O2 (t), 0) = (η 1 (t),η 2 (t), 0) and the instant angular velocity is ω(t) = (ω 1 (t), ω 2 (t), ω 3 (t)) = (η 4 (t), η 5 (t),η 6 (t)). The Oxyz coordinate system is rigidly fixed with the container. The mean (hydrostatic) free surface Σ 0 belongs to the Oxy plane and the origin is in the geometric centre of Σ 0 . In the non-dimensional statement, the tank breadth and width are equal to the unit but the non-dimensional mean liquid depth equals to h.
Each case in § § 4-7, especially, in § 4, which is associated with orbitally shaken containers, deserves an additional thorough theoretical study, perhaps, an independent publication. In order to restrict the paper volume, we report only a generic description of what kind of the steady-state waves are possible by using the reliable approximate mathematical model, which was supported by model tests for other classes of the resonant forcing. The present paper should also guide future experimental studies with a focus on co-existing stable resonant waves established for all forcing types.

The Narimanov-Moiseev-type modal equations
A cyclically moving rigid square base container is partially filled with a finite liquid depth by a perfect incompressible (irrotational flows) liquid. Figure 1 introduces necessary notation including the prescribed time-periodic translatory v O (t) = (η 1 (t),η 2 (t), 0) and instant angular ω(t) = (η 4 (t),η 5 (t),η 6 (t)) velocities of the container, which satisfy the periodicity condition v O (t + T) = v O (t) and ω(t + T) = ω(t), where T = 2π/σ is the forcing period. The developed nonlinear sloshing theory below requires that the heave container component is zero (η 3 (t) ≡ 0) to avoid parametric resonances in the hydrodynamic system (see a recent review on the parametric waves by Ibrahim (2015)).
Oscillatory liquid motions (sloshing) are considered in the non-inertial (tank-fixed) coordinate system Oxyz whose Oxy-plane coincides with the mean (hydrostatic) free surface Σ 0 and Oz runs through the centre of Σ 0 . Mathematically, the sloshing problem consists of finding the two unknowns: the free surface Σ(t) (described by the single-value representation z = f (x, y, t)) and the absolute velocity potential Φ(x, y, z, t). The corresponding free-surface problem and its variational analogy are presented in Part 1 (see more mathematical and physical details in chap. 2 of Faltinsen & Timokha (2009)).
To facilitate interested applied mathematicians and engineers, the supplementary materials A, available at https://doi.org/10.1017/jfm.2020.253, report computational formulas for the hydrodynamic coefficients and tables them for some particular values of h. An interesting mathematical fact is that η 6 (t) (yaw) is not presented in the Narimanov-Moiseev-type equations (2.6)-(2.8). This means that the yaw-type resonant forcing of the lowest natural sloshing frequency σ 1 = σ 0,1 = σ 1,0 exclusively affects the linear sloshing component in representation (2.2), i.e. only sway, pitch, roll and surge matter when σ is close to σ 1 . As was discussed in Part 1, the third-order generalised coordinates are 'driven', i.e. equations (2.6)-(2.7) do not depend on a 3 , b 3 , c 12 and c 21 but (2.8) contain only linear terms by these generalised coordinates. These generalised coordinates and (2.8) are a necessary mathematical component of the Narimanov-Moiseev theory, which requires to account for all the O( )-quantities in derivations. Part 4 also showed that contribution by a 3 , b 3 , c 12 and c 21 is necessary to get a good agreement with experimental measurements. Because damping was important in Part 4 to fit experimental data by Ikeda et al. (2012), the modal system is equipped with the underlined linear damping terms. The damping ratios ξ i,j express a summarised ('integral') effect of different physical factors, which include the laminar boundary layer on the wetted tank surface (Faltinsen & Timokha 2009, equations (6.139)). Part 4 also discusses when the linear damping terms can be utilised by the Narimanov-Moiseev-type modal equations, whose derivation suggests an ideal liquid with irrotational flows. A requirement is that the viscous vortical flows are localised near the wetted bottom/walls and/or the free surface. This includes the aforementioned boundary layer and, possibly, the effect of breaking waves. However, the damping ratios ξ i,j should remain asymptotically small values. Part 4 concludes that, mathematically, ξ 1,0 = ξ 0,1 = O( 2/3 ) to match the Narimanov-Moiseev asymptotic relations. The same order for the 'higher' damping rates means that the linear terms in (2.7) and (2.8) are O( 4/3 ) and O( 5/3 ), respectively, and, therefore, these terms can, as in Part 4, be neglected in the forthcoming derivations of the periodic a 2 , b 2 , c 1 and a 3 , b 3 , c 12 , c 21 .
Using the Narimanov-Moiseev-type modal system (2.6)-(2.8) facilitates the so-called classification of the steady-state resonant wave regimes (modes). The classification suggests: (i) derivation of the asymptotic (T = 2π/σ )-periodic solutions and examining their stability, (ii) specifying a link between these solutions and the corresponding steady-state wave modes and, finally, (iii) investigation of the wave-amplitude response curves identifying, in a parallel way, which stable steady-state wave mode(s) correspond to each point on the curves and how the wave modes change versus the forcing frequency σ /σ 1 and other typical input parameters. In the studied case, these parameters are the semi-axis ratio 0 |δ 1 | 1, angular position and direction of the horizontal translatory elliptic tank motion. The theoretical classification for the semi-axis ratio δ 1 = 0 was successfully done in Part 4. These results were supported by experiments by Ikeda et al. (2012).

Asymptotic periodic solution
The asymptotic derivation scheme from Parts 1 and 4 can be modified to get a general asymptotic T-periodic solution of (2.6)-(2.8) for three-dimensional non-parametric tank excitations. The scheme should, as usual, start with the lowest-order component which, according to the Narimanov-Moiseev theory, describes the asymptotically dominant contribution of the primary excited Stokes modes, f (1) 1 (x) and f (2) 1 (y), in the modal representation (2.2).
The first Fourier harmonics in a 1 (t), b 1 (t) are the only lowest-order asymptotic components of the steady-state solution within the framework of the Narimanov-Moiseev asymptotic theory. Inserting (3.1) into (2.7) derives the O( 2/3 )-components of a 2 (t), b 2 (t) and c 1 (t), which are quadratic functions of a,ā,b, b (these expressions are identical to those in Part 4). Furthermore, substituting (3.1) and a 2 (t), b 2 (t), c 1 (t) into (2.6) and (2.8) makes it possible to find the periodic solution of the entire modal system (2.6)-(2.8) up to the O( ) terms. Interested readers are referred to the supplementary materials B for more analytical details. The main result is the solvability (secularity) condition with regard to a,ā, b,b, Three-dimensional sloshing in a square. Part 5 894 A10-9 which contains for the three-dimensional forcing the four non-dimensional forcing amplitudes x ,¯ x , y ,¯ y = O( ) associated with the lowest Fourier harmonics on the right-hand side of (2.6a) and (2.6b), i.e.
where K x (t) = −P 1,0 (η 1 − S 1,0η5 − gη 5 ), K y (t) = −P 0,1 (η 2 + S 0,1η4 + gη 4 ) and ξ = 2ξ 0,1 = 2ξ 1,0 . Because we focus on resonant excitations of the lowest natural sloshing frequency, at least one of x ,¯ x ,¯ y and y should not be zero; specifically, all four nondimensional forcing amplitudes could be non-zero in the most general case. The nondimensional coefficients m 1 , m 2 and m 3 are functions of the non-dimensional liquid depth h. Computational formulas for m 1 , m 2 and m 3 are given in the supplementary materials B. Their values satisfy where h * = 0.3368 . . . is the critical depth (Hermann & Timokha 2008). Part 4 derived a similar secular system for the oblique forcing, in which¯ x = y = 0. As we will show, the more complicated right-hand side with y = 0 significantly changes the analytical properties of (3.2) and its solutions. This complicated right-hand side makes inapplicable analytical and numerical results from Part 4, which had to be modified to study the three-dimensional forcing.
The secular system (3.2) determines the lowest-order non-dimensional amplitude parameters a,ā, b,b as functions of the non-dimensional frequency parameter Λ and, therefore, it computes the first-order component of steady-state waves expressed in terms of the hydrodynamic generalised coordinates (3.1). Furthermore, the second-order wave component is associated with the generalised coordinates a 2 , b 2 and c 1 , which are governed by the homogeneous modal equations (2.7). As shown in Part 4, exact analytical expressions for a 2 , b 2 and c 1 are uniquely quadratic functions of a,ā, b,b. The third-order generalised coordinates a 3 (t), b 3 (t), c 12 (t) and c 21 (t) do not appear in the subsystem (2.6)-(2.7); they follow from (2.8). This explains why the O( )-order asymptotic component (in contrast to the O( 2/3 )-order contribution) of the constructed asymptotic periodic solution does not affect the lowest-order non-dimensional amplitude parameters a,ā, b,b. Finally, the supplementary materials B reproduce another analytical result from Part 4, which shows, in particular, that Lyapunov's stability of the constructed asymptotic periodic solution exclusively depends on a,ā, b,b and Λ. In summary, solving the secular system (3.2) for the fixed forcing amplitude parameters x ,¯ x ,¯ y , y gives the lowest-order wave-amplitude parameters a,ā, b,b versus the forcing frequency parameter Λ, determines the first-and second-order components of the corresponding steady-state waves and the waves stability. The higher-order Fourier harmonics in η i (t), i = 1, 2, 4, 5, 6 only contribute to the highest-order components of the steady-state wave flows.

Asymptotically equivalent orbital tank motion
Concentrating only on the first-and second-order approximations and Lyapunov's stability of the corresponding steady-state wave regimes reduces the external excitation effect to the four forcing amplitude parameters x ,¯ x ,¯ y and y appearing on the right-hand side of (3.2). Taking them instead of the original periodic sway, surge, roll, pitch and yaw introduces the prescribed out-of-phase sway-and-surge tank motion η * 1 (t) = e x cos σ t +ē x sin σ t, η * 2 (t) =ē y cos σ t + e y sin σ t, η * 4 (t) = η * 5 (t) ≡ 0, x = P 1 e x ,¯ x = P 1ēx ,¯ y = P 1ēy , y = P 1 e y ; P 1 = P 1,0 = P 0,1 < 0. (3.5) When using the Narimanov-Moiseev asymptotic theory, the resonant translatory forcing (3.5) yields the same secular system (3.2), derives the same first-and second-order components of the asymptotic periodic solution and keeps, specifically, the same stability properties. Under these circumstances, the combined out-of-phase sway-and-surge excitation (3.5) could be considered as an asymptotically equivalent tank forcing to the given periodic sway, surge, roll, pitch and yaw. The equivalence provides an efficient parameter study of the steady-state waves for a given combined periodic sway, surge, roll, pitch and yaw. The three-dimensional tank 'trajectories' are then characterised by an infinite set of input parameters, but the equivalent orbit (3.5) operates with the four real numbers e x ,ē x ,ē y , e y , only three of which, as we will show, are independent.
Keeping the lowest Fourier harmonics in (3.5) of the generic periodic sway, surge, roll, pitch and yaw is a serious advantage of the Narimanov-Moiseev asymptotic theory if the task consists of classifying the stable/unstable resonant waves due the resonant forcing of the lowest natural sloshing frequency. For example, the horizontal tank motion along the Bernoulli lemniscate of the non-dimensional width l a by η 1 (t) = l a cos σ t/(1 + sin 2 σ t) = 2l a ( √ 2 − 1) cos σ t + · · ·, η 2 (t) = l a sin σ t cos σ t/ (1 + sin 2 σ t) = 2l a (3 − 2

√
2) sin 2σ t + · · · is, according to this analytical result, equivalent to the longitudinal excitation by η 1 (t) = 2l a ( √ 2 − 1) cos t, which was considered in Part 1. For an oblique position of the Bernoulli lemniscate, one should apply results from Part 4. Difference between the steady-state sloshing due to original (lemniscate) and longitudinal trajectories is the O( )-order wave component, which, in addition, does not affect the hydrodynamic stability.
Another important consequence of the asymptotic equivalence is that it mathematically argues why the suspended container (Cooker 1994), which yields the purely harmonic forcing (3.5) for small-magnitude pendulum motions, could be employed (see, e.g., Turner, Bridges & Ardakani (2015a,b)) for modelling the non-parametric resonant vessel behaviour in a periodic sea, which can deal with rather complicated periodic sway, surge, roll, pitch and yaw.
We must stress that the proven equivalence requires the resonant forcing, for which the sloshing amplitude is of a lower asymptotic order than the forcing. Mathematically, it is not possible in the linear (non-resonant) approximation. However, the Narimanov-Moiseev theory is only a very particular case when the equivalent orbits exist. More complex resonances require the adaptive modal theory for which the orbital equivalence is, intuitively, allowed but requires a dedicated study.
When e x e y −ē xēy = 0, equation (3.5) implies a harmonic (longitudinal, oblique or diagonal) horizontal tank reciprocation. This particular case was studied in Part 4. If e x e y −ē xēy = 0, the tank moves translatory along a centred elliptic orbit, which is governed by The centred elliptic tank orbit is defined by its semi-axes |e x | and |e y | (without loss of generality, 0 |e y | |e x | = 0), angle γ between the major axis and Ox as well as by the angular direction (clockwise or counterclockwise) as schematically illustrated in (a). Panel (b) depicts the circular orbit (e x = e y ). Panels (c) and (d) show the two limiting cases, for which the ellipse axes are parallel to the symmetry planes associated with either vertical walls or diagonals, respectively. The two-sided arrows indicate the reciprocating harmonic tank excitations whose effect on the steady-state sloshing was studied in Part 4.
(readers can easily prove that by computing the fundamental invariants). Any centred elliptic trajectory is fully determined by (i) its semi-axis lengths |e x | and |e y | (without loss of generality, 0 |e y | |e x | = 0), (ii) the angle γ between the major axis and Ox and, finally, (iii) the trajectory direction (clockwise or counterclockwise). The trajectory is schematically depicted in figure 2(a). Assuming |e x |, |e y | and γ are known (these can be computed from (3.6)) leads to the generic parametric representation η * 1 (t) = [e x cos γ ] cos σ t − [e y sin γ ] sin σ t, η * 2 (t) = [e x sin γ ] cos σ t + [e y cos γ ] sin σ t, (3.7a,b) which should, to within the phase lag by t, be equivalent to (3.5). Here, e x = 0 and e y are not necessarily positive numbers to define the tank trajectory direction as follows e x e y > 0 -counterclockwise, e x e y < 0 -clockwise, e x e y = 0 -reciprocation.
3.3.1. Alternative forms of (3.13) An alternative secular system can be obtained after excluding the phase lags ψ and ϕ. The first (rather obvious) step of the derivation consists of taking 1 2 + 3 2 and 2 2 + 4 2 , which leads to (3.17) Further, we substitute ϕ = ψ + α into 2 and 4 and insert expressions for ( x cos ψ) and ( x sin ψ) from 1 and 3 . The result is the linear algebraic homogeneous whose solvability (zero-determinant condition) yields the third equation Even though this derivation step is based on analytical ideas from Part 4, the resulting (3.21) and (3.21c) differ from their analogy in Part 4 (equations (5.4) and (5.6)). This is especially clearly seen for (3.19), which depends now on three forcing amplitude parameters and has much more complicated analytical structure. This means that the three-dimensional forcing makes inapplicable numerical schemes from Part 4. The system (3.17), (3.19) (of three equations) couples the four unknowns Λ, A 2 , B 2 and α. Obviously, if the tetrad (Λ, A, B, −π α = ϕ − ψ π) satisfies (3.13), it automatically satisfies (3.17) and (3.19). However, because (in contrast to (3.13)) EQ1, EQ2, EQ3 are the π-periodic functions by α, the system (3.17), equation (3.19) cannot identify the original phase lag difference α but outputs α = α 0 + πk, k ∈ Z, where, without loss of generality, −π/2 α 0 π/2. The problem on restoring the original α from α 0 can be solved after noting that F(α) and D(α) on the left-hand side of (3.13) are also the π-periodic functions of α. This means that substituting α = α 0 + πk, k ∈ Z into F(α) and D(α) of (3.13) correctly computes their left-hand sides, and, because we assume x Υ = 0, one can restore the phase lags ϕ and ψ from (3.13) as well as the original phase lag difference −π α = ϕ − ψ π. Bearing in mind the latter analytical manipulation, one can conclude that the alternative secular system (3.17), (3.19) is mathematically equivalent to (3.13).
In order to find A, B and α, Part 4 (the reciprocating forcing, equations (5.4) and (5.6)) developed a specific numerical scheme, which is based on excluding the trigonometric terms with α and replacing them by C in (3.16). As we stated above, the three-dimensional forcing makes these secular equations and their solutions much more complicated. Numerical tests established situations when several different solutions exist with small but not zero C and, simultaneously, for the same input parameters, the system has solutions with finite and very large C. An accurate detection of these different C becomes a numerically consuming task. Thus, analytical ideas from Part 4 are not applicable. Computations have to adopt another parameter instead of the trigonometrical terms with α. Hereafter, we introduce 3.3.2. Numerical solutions of (3.21) The secular system (3.21) will be used for analytical and numerical studies. These suggest considering the following two cases. The first case EQ1 and, therefore, the unknowns A, β and B belong to the three-dimensional curvilinear trapezoidal (3.23a−c) where, as we see, the lower and upper bounds of A and β are fixed (do not depend on B) but the upper bound for B depends on A and β. Now, we can analytically solve EQ1 ± = 0 with respect to Λ for (A 2 , β; B 2 ) belonging to (3.23), Substituting (3.24) into (3.21b) and (3.21c) yields the system with respect to the three unknowns A 2 , β and B 2 , which belong to (3.23). The system (3.25) can be solved to find A 2 and B 2 as functions of β.
When considering the second case, Three-dimensional sloshing in a square. Part 5 894 A10-15 the unknowns belong to the curvilinear three-dimensional trapezoidal so that the problem reduces to the system ) with respect to the three unknowns A 2 , β and B 2 from (3.27). The system also determines A 2 and B 2 as functions of β.
One must note that the proposed numerical scheme provides, in fact, exhaustive search for all admissible solutions, which, as we showed, belong to the curvilinear trapezoidals (3.23) and (3.27). Thinking in terms of the wave-amplitude response curves, this means that it effectively detects and does not miss out both continuous branches existing in the entire resonant zone and locally distributed loop-type branches, which, as our computations will show, may emerge from a single point for certain sets of the input parameters.
Part 4 demonstrated applicability of the Narimanov-Moiseev-type modal system for the non-dimensional input parameters (3.30) which are associated with experimental model tests by Ikeda et al. (2012) who studied the steady-state resonant sloshing for the reciprocating (longitudinal, diagonal and oblique) harmonic tank motions. The input parameters are rather typical when the Narimanov-Moiseev theory is applicable. This includes the mean liquid depths 0.5 h (m 1 , m 2 and m 3 weakly change for those depths, see the supplementary materials B), relatively small forcing amplitudes (based on the authors experience) and the damping coefficient ξ , at least, for the laboratory tanks whose horizontal dimensions are less than one meter. For larger tanks, the viscous boundary layer effect on the wetted tank surface decreases but the 'integral' damping coefficient ξ (we assumed that the damping rates integrally account for all dissipative factors in the hydrodynamic system) become then significantly contributed by the local free-surface phenomena, which are discussed, e.g., in Part 1. Numerical examples in the forthcoming sections will, therefore, adopt (3.30).

Steady-state wave modes
Solving (3.25) and (3.29) and utilising (3.24) and (3.28), respectively, outputs (Λ, A 2 , B 2 , β). Later on, by using (3.22) and (3.26) to compute D and F, we restore the phase lags ψ and ϕ from (3.13) and compute α = ϕ − ψ. Finally, equation (3.14) returns the wave-amplitude parameters a,ā,b, b, which determine the steady-state wave Using (3.31) makes it possible to define the types of the corresponding steady-state waves. The first type implies the standing wave mode (particular cases are planar, diagonal and square-like; see details and discussions in Part 1 in the supplementary materials B). It occurs when surface patterns by S(x, y; a,b) and S(x, y;ā, b) in (3.31) are geometrically similar, i.e.
and, therefore, the lowest-order component in (3.31) admits the following form When (3.32) is not satisfied, the single-value representation z = S(x, y; a,b) cos σ t + S(x, y;ā, b) sin σ t determines the swirling wave mode, i.e. an azimuthally propagating wave running in either clockwise or counterclockwise direction around Oz. Swirling (rotary wave) for the square base tank slightly differs from what is swirling for the circular base tank, when the lowest-order approximation alike (3.31) may be re-written in a mathematical form, which clearly identifies an azimuthal wave propagation. Part 1 defines the swirling wave mode for the square base tank as where 'an almost flat crest travels around each of the four sides with an almost flat trough on the opposite side'; it also clarifies this definition by photos from the dedicated authors model tests, which are reproduced in the supplementary materials B. A mathematical treatment of the swirling wave mode can be done by considering the nodal line in the lowest-order approximation, i.e. the time-dependent curve on Σ 0 implicitly defined by the equation S(x, y; a,b) cos σ t + S(x, y;ā, b) sin σ t = 0. The motionless nodal line corresponds to a standing wave. When (3.32) is not fulfilled, the nodal line changes in time. Moreover, an elementary analysis shows that the nodal line steadily rotates (with deformations) around the Oz axis and the rotation direction depends on the sign of (3.34) so that, because A, B > 0, one can discriminate counterclockwise/clockwise swirling and standing waves as follows sin α > 0 -counterclockwise, sin α < 0 -clockwise, sin α = 0 -standing.
(3.35a−c) To explain experimental observations by Ikeda et al. (2012) and account for the asymptotic character of (3.31), Part 4 had to introduce the so-called nearly standing wave mode. Mathematically, it occurs when S(x, y; a,b) and S(x, y;ā, b) determine similar surface patters on the adopted asymptotic scale O( 1/3 ), i.e.
) and, as a consequence, the lowest-order component in (3.31) rewrites in the form Accounting for the latter wave type, the resonant steady-state waves fall apart into the purely and nearly standing wave modes, the swirling (rotary) wave mode, which . Symbolic classification of the steady-state wave modes by using (3.39). The criterion (3.40) makes it possible to distinguish co-and counter-directed (with respect to the orbital tank motions) swirling waves. The introduced symbols are further used as markers on the wave-amplitude response curves. Visualisation of the steady-state wave modes is given in Part 1 and the book by Faltinsen & Timokha (2009). We incorporate some of them into the supplementary materials.
occurs either clockwise or counterclockwise as well as irregular waves (all steadystate wave regimes are unstable and chaotic, modulated, etc. waves are realised). The irregular waves were detected in experiments from Part 1; they are also numerically confirmed by using the Narimanov-Moiseev and adaptive modal theories in Part 2 (Faltinsen, Rognebakke & Timokha 2005) for the longitudinal and diagonal forcing cases. The following criterion helps to identify them all steady-state waves are unstable ⇔ irregular (chaotic/modulated) sloshing.
The steady-state wave modes are symbolically classified in figure 3 where, due to (3.12) and (3.39), we can introduce swirling, which is co-or counter-directed with the elliptic orbit. The appropriate criterion reads as The forthcoming steady-state analysis operates with the wave-amplitude response curves in the (σ /σ 1 , A, B)-space. Different line types mark stable and unstable solutions but the criteria (3.39) and (3.40) are employed to specify/analyse the swirling wave direction. Figure 3 supplements the wave type classification.

Circular orbit
The circular tank forcing (figure 2b schematically depicts that) implies in (3.10) and (3.21a), (3.21b), where |e x | is the orbit radius, and the third secular equation (3.21c) takes the form One can see that (4.1) makes the secular system (3.21a), (3.21b), (4.2) independent of the sign of δ 1 . This means that its solutions (the tetrads (Λ, A 2 , B 2 , β)) and, therefore, the wave-amplitude response curves in the (σ /σ 1 , A, B) space do not depend on the circular direction (counter-or clockwise). However, using (3.13), (4.1) to restore the phase lags ϕ and ψ shows that the phase lag difference α is affected by the sign. This means that one can adopt, without loss of generality, δ 1 = 1 in calculations but we should use the criterion (3.40) to distinguish co-and counter-directed stable swirling waves.
Adopting (4.1) in the secular system (3.21a), (3.21b), (4.2) derives the following particular analytical solution inserting this solution into (3.13) (to get ϕ and ψ) deduces sin α = δ 1 . According to (3.39) and (3.40), (4.3) implies the swirling wave mode, whose angular propagation coincides with the circular orbit. Numerical analysis of the secular system (3.21a), (3.21b), (4.2) within (4.1) and β = 0 detects the steady-state sloshing with A = B and β = 0. This is opposite to the circular base tank considered by Raynovskyy & Timokha (2018a) when m 2 = m 3 in the corresponding secular equations. Mathematically, one can show that this solution is impossible in the considered case, too, if the identity m 2 = m 3 holds true. However, table 3 in the supplementary materials shows that this never happens for the square cross-section and finite liquid depths. The exhaustive search detected the aforementioned solution, using a very thick mesh with a refinement procedure, in a narrow interval of β, where 0.995 < β 1. In this interval, 1 − β 0.995 O( 2/3 ) on the adopted asymptotic scale (3.30) and, therefore, the found solutions imply either standing (when β = 1) or nearly standing wave modes. Specifically, Hiramitsu & Funakoshi (2015) did not find out those solutions in their numerical studies devoted to the circular tank forcing; they only numerically identified our analytical solution (4.3). Because of a very narrow range where these solutions exist, this could be a computational problem. Figure 4 demonstrates the wave-amplitude response curves for the circular forcing with the input data (3.30). Panel (a) shows a three-dimensional view but panels (b) and (c) are projections of the three-dimensional branching on the (σ /σ 1 , A) and (σ /σ 1 , B) planes, respectively. Each point of the response curves corresponds to FIGURE 4. The steady-state wave-amplitude response curves associated with the circular orbital tank motion schematically depicted in figure 2(b). The solid lines mark stable solutions but the dashed (deep-blue) lines are used to indicate the instability. Panel (a) shows the response curves in the (σ /σ 1 , A, B) space but (b) and (c) demonstrate their projections on the (σ /σ 1 , A) and (σ /σ 1 , B) planes, respectively. The branch P l S 0 WP r implies swirling whose angular wave propagation coincides with direction of the circular orbit. The loop-type response curves (symmetric with respect to the A = B plane) imply the nearly standing (square-like) wave mode but points E,Ē and P 0 ,P 0 correspond to the (purely standing) square-like waves. Symbols from figure 3 and criteria (3.39), (3.40) are used to specify stable standing, nearly standing and swirling waves. The input data (3.30) are adopted in the computations.
one (unique) steady-state wave whose lowest-order amplitude components in the Ox and Oy directions are equal to A and B, respectively. The solid lines denote stable steady-state solutions but the dashed lines are used to mark the instability. The stable steady-state wave modes are marked by symbols from figure 3. The swirling wave mode (analytical solution (4.3)) is represented by points on the continuous curve P l S 0 WP r . The curve belongs to the plane A = B in the (σ /σ 1 , A, B) space. It demonstrates the well-known damped hard-spring-type behaviour, that is, when the damping coefficient ξ decreases, point S 0 goes to infinity. As it typically happens for the hard-spring type branching, stable solutions (here, stable swirling waves) correspond to points on P l S 0 and WP r . This means that two different stable swirling waves are theoretically possible in the frequency range determined by W (the left bound) and S 0 (the right bound). According to (3.40), swirling on the branch P l S 0 WP r is co-directed with the circular tank orbit.
Standing (points P 0 ,P 0 and E,Ẽ where β = 1) and nearly standing waves are represented in figure 4 by the two loop-type branches EU P 0 D 0 E (B < A) and EŨ P 0D 0Ẽ (A < B). These branches are symmetric relative to the plane A = B in the (σ /σ 1 , A, B) space so that substitution A → B, B → A transforms EU P 0 D 0 E tõ EŨ P 0D 0Ẽ and visa-versa. This symmetry reflects invariance of the steady-state wave modes relative to the π/2-rotation of the Oxy coordinate system around Oz. There exist either four or eight (panels (b) and (c) make it more precise) solutions of the secular system, which imply the nearly standing (nearly square-like) wave mode in a certain frequency range with σ /σ 1 < 1. However, only two (on D 0 U andD 0Ũ ) of these four/eight solutions are stable. Positions of P 0 (P 0 ) and D 0 (D 0 ) depend on the damping coefficient ξ . They go to infinity as ξ → 0.
In summary, the circular orbital forcing always excites the stable resonant steady-state swirling whose angular propagation (around Oz) coincides with the forcing direction (clockwise or counterclockwise). For σ /σ 1 > 1, two different swirling waves may coexist in a certain frequency range whose upper bound strictly depends on damping. For σ /σ 1 < 1, the swirling wave co-exists (in a local frequency range) with two stable nearly standing square-like waves. The square-like waves are characterised by a relatively large amplitude relative to the wave amplitude of the co-existing stable swirling. Part 1 (figure 10) discovered a similar situation for the longitudinal forcing when the square-like waves of larger amplitude co-existed with stable planar wave of lower amplitude. Attempts to find out these square-like waves in experimental model tests (with almost zero initial scenarios) have failed except for h = 0.34. The latter became possible since smaller liquid depths lead to a narrow frequency range where only square-like waves are stable (planar wave is unstable). Unfortunately, such a frequency range is not possible for the circular forcing, even for smaller depths. Stable swirling always co-exists with the square-like waves for any resonant frequencies with σ /σ 1 < 1. As a consequence, the only way to experimentally detect the nearly standing square-like waves for circular excitations is to adopt exotic non-zero initial scenario. This is rather difficult to organise. From this point of view, the stable nearly standing waves on EU andẼŨ may be treated as almost unphysical; probability of their experimental occurrence looks very low.

Wall-symmetric elliptic orbits
Dealing with the so-called canonic position of the elliptic orbits in figure 2(c) mathematically implies in (3.10) and (3.21a), (3.21b), where |e x | is the major semi-axis length; (3.21c) ⇒ For (5.1), the secular system (3.21a), (3.21b), (5.2) becomes depending on the modulus |δ 1 | (the semi-axis ratio) but independent of the sign of δ 1 . This means that the tetrads (Λ, A 2 , B 2 , β) and the corresponding response curves in the (σ /σ 1 , A, B) space do not depend on the forcing direction (counterclockwise or clockwise) but, since the actual phase lags ϕ, ψ and their difference α are affected by the sign, results on the swirling wave mode should be interpreted by using the criterion (3.40), i.e. one must discriminate counter-directed and co-directed swirling waves. Symbols from figure 3 will be adopted to mark that feature on the response curves.
A transformation ('metamorphosis') of the response curves versus the semi-axis ratio 0 |δ 1 | < 1 is shown in figures 5 and 6. The calculations used the input data (3.30). The solid lines imply stable steady-state solutions but the dashed (deep-blue) lines mark instability. All notation is taken (with minor additions) from Part 4. Symbols from figure 3 specify the stable steady-state sloshing modes.
The starting point in figure 5 is the wave-amplitude response curves for the longitudinal harmonic tank excitations (δ 1 = 0). Part 4 derived analytical solution for that case, which is reported in the supplementary materials C. According to this result, the longitudinal forcing causes planar standing (β = 1, B = 0), swirling and nearly standing (nearly square-like) waves; a non-empty frequency range of irregular waves exists. The wave-amplitude branching (relevant to experiments by Ikeda et al. (2012)) is reproduced in figure 5(a). It is numerically (but not visually) identical to that by Faltinsen & Timokha (2017, figure 4a). Figure 5(a) contains a three-dimensional view in the (σ /σ 1 , A, B) space and two projections on (σ /σ 1 , A) and (σ /σ 1 , B), respectively. The response curves consist of the branch P l TEP 0 WP r (it belongs to the (σ /σ 1 , A) plane) and the arc-type branch ED 0 UVS 0 E. They are linked with each other at bifurcation points E and W. The branch P l TEP 0 WP r is responsible for planar waves; it exhibits the damped soft-spring type behaviour.
Stable solutions on the arc-type branch ED 0 UVS 0 W imply either swirling (VS 0 ) or nearly standing (D 0 U) steady-state waves. An important note is that each point on this branch corresponds to the two Ox-symmetric waves (in contrast to P l TEP 0 WP r ). For VS 0 , this symmetry means that there exist two physically identical rotary waves propagating in opposite (clock-or counterclockwise) directions. The two stable Ox-symmetric nearly standing (square-like) waves for each point on D 0 U (calculations give 1 = 1 − β 0.95 O( 2/3 ) with (3.30)) could be interpreted as two asymptotic solutions (3.37) with γ 1 , ±γ 2 (and, generally, different t 0 ). The planar (standing) waves are unstable between abscissas of T and W where stable swirling (associated with VS 0 ) and irregular waves are expected. The planar wave mode coexists with swirling in the frequency range defined by W and S 0 but the nearly square-like waves are realisable between abscissas of D 0 and T (theory and experiments in Part 1 showed that U moves into the zone between T and W for h 0.4). When the damping coefficient ξ tends to zero, points D 0 and S 0 go to infinity. This may theoretically enlarge the aforementioned zones of stability for the nearly standing and swirling wave modes. Figure 5(b) illustrates a metamorphosis of the response curves from column (a) when the semi-axis ratio |δ 1 | is relatively small but not zero. We were not able to find analytical (even particular) solutions for that case. The numerical scheme from § 3 was implemented to solve the secular system (3.21a), (3.21b), (5.2) equipped with (5.1). Computations suggested that β = 0 as 0 < |δ 1 | < 1 (substituting β = 0 into (5.2) ⇒ A = B ⇒ |δ 1 | = 1 from (3.21a), (3.21b) with (5.1)) but β = 1 (standing waves) was allowed. The response curves in figure 5(b) are drawn for the input data (3.30). They show that considering |δ 1 | = 0.05 splits the connected curve from the column (a) into the two non-connected branches, P l TD 0 UVS 0 WP r continuously going from P l to P r and the loop-type branch EP 0 S 0 V U D 0 E. A difference between (a) and (b) is that each The steady-state wave-amplitude response curves in the (σ /σ 1 , A, B) space for the input parameters (3.30). The solid lines specify stable solutions but the dashed (deepblue) lines mark their instability. The column (a) is devoted to the longitudinal harmonic forcing (δ 1 = 0, the results are taken from figure 4a by Faltinsen & Timokha (2017)). The branch P l TEP 0 WP r corresponds to the planar wave mode (stable on P l T and WP r ). Each point on the arc-type branch ED 0 UVS 0 W implies two steady-state waves including the stable nearly standing waves (D 0 U) and swirling (VS 0 ), which differ from each other only by the angular propagation direction. There are no stable solutions in the frequency range between T and V where an irregular sloshing is expected. The column (b) shows the response curves for |δ 1 | = 0.05 (relatively small semi-axis ratio). The response curves from (a) split for the case (b) into the two non-connected branches P l TD 0 UVS 0 P r and (loop type) EP 0 S 0 V D 0 E. The standing wave mode disappears except at points E and P 0 (unstable). Stable nearly planar waves (because of the relatively small |δ 1 |) are associated with points on P l T and WP r . The stable swirling from (a) splits into counter-(V S 0 ) and co-(VS 0 ) directed swirling modes in (b). Stable nearly square-like waves (different, in contrast to the case (a)) are associated with points on D 0 U and D 0 U . The frequency range of irregular waves is still determined by T and V.  point on P l TD 0 UVS 0 WP r and EP 0 S 0 V U D 0 E in (b) corresponds to a unique periodic solution and, therefore, implies only one steady-state wave. Indeed, each point on the arc-type branch ED 0 UVS 0 E in column (a) implies a pair of the Ox-symmetric steady-state waves but the non-zero |δ 1 | (passage to a wall-symmetric elliptic tank orbit) breaks this symmetry. As a consequence, the column (b) deals with two different stable nearly standing (1 − β 0.95 = O( 2/3 ) in calculations with (3.30)) waves on D 0 U and D 0 U and, of course, two stable swirling waves on VS 0 and V S 0 , respectively. These swirling waves have not only different angular propagation directions but also different wave amplitudes. Similar splitting into co-and counter-directed swirling was detected in computations by Hiramitsu & Funakoshi (2015) for a square base tank as well as described by Raynovskyy & Timokha (2018a) and  for the circular base container. However, sloshing in the circular base basin does not allow for the stable nearly standing wave mode. Points on P l T and WR r were associated with stable planar standing waves in column (a) but they become responsible for the nearly standing (nearly planar) waves in column (b). Positions of T and W are weakly affected by the non-zero |δ 1 | = 0.05 and, therefore, the frequency range of irregular waves remains practically the same. The limit ξ → 0 will shift points P 0 , D 0 (D 0 ) and S 0 (S 0 ) far away from the primary resonant zone. Increasing the semi-axis ratio to |δ 1 | = O(1) (on the chosen asymptotic scale) converts the nearly standing waves on P l T and WR r to the co-directed swirling. Figure 6 contains four panels, which demonstrate this and other metamorphoses of the wave-amplitude response curves. Only three-dimensional views are shown. The values |δ 1 | = 0.3, 0.5, 0.7 and 0.95 are used as representing most impressive peculiarities of the response curves. Notation, mathematical and physical interpretations of the 'stable' branches and the corresponding steady-state wave modes are adopted from figure 5(b). Two response curves in (a-c) look topologically identical to those in the column (b) of figure 5, i.e. there are the continuous branch P l TD 0 UVS 0 WP r going through the entire resonant zone and the loop-type branch D 0 EP 0 U D 0 . Differences between these panels become clear when monitoring types and stability ranges of the steady-state waves. As we stated above, nearly standing (planar) waves on P l T and WR r become impossible when |δ 1 | = O(1). Periodic solutions on P l T and WP r as well as D 0 U (nearly standing waves) covert then to the stable swirling, which is co-directed with the elliptic tank orbit.
The counter-directed swirling disappears when the semi-axis ratio |δ 1 | exceeds 0.5 for the input data (the same critical value was estimated in numerical experiments by Raynovskyy & Timokha (2018a) for the circular-base basin). A particular consequence of this fact is that the frequency range of irregular waves (TV) narrows with increasing |δ 1 |; the range becomes empty in panel (d) drawn with |δ 1 | = 0.95 (the tank orbit is close to the circle). The branching is then topologically identical to that in figure 4. Its principal 'novelty' with regard to (a-c) is the loop-type branchD 0ẼP 0Ũ responsible for the nearly standing wave mode (stable for points onD 0Ũ ). This branch emerges 'from a single point' at |δ 1 | = 0.8 (for the input data (3.30)). Hiramitsu & Funakoshi (2015) report a few calculations for the wall-symmetric elliptic orbits and demonstrate how an energetic parameter √ E changes versus the frequency parameter. Qualitatively, these calculations lead to similar branching (for √ E) but differ from our results with |δ 1 | approaching the unit. This difference was discussed for the circular orbit as, possibly, insufficient accuracy in calculations by Hiramitsu & Funakoshi (2015) since the nearly standing wave solution occurs in a very narrow range by β.
In summary, passage from the harmonic longitudinal reciprocation to elliptic tank (canonic) orbits in the Oxy plane (figure 2b), first, breaks away the Ox symmetry of nearly standing and swirling wave modes so that, in particular, swirling waves have different amplitudes depending on their angular propagation direction. Second, a further increase of the semi-axis ratio |δ 1 | eliminates, starting with a critical value of |δ 1 | (estimated at ≈0.5 for the tested input data), swirling whose angular propagation is opposite to the tank orbit (similar to the circular base basin studied by Raynovskyy & Timokha (2018a)). Third, the co-directed swirling waves become stable in the entire resonant zone and, as a consequence, irregular waves are not predicted when the tank trajectory approaches a circle. The latter limit topologically changes the wave-amplitude branching by yielding an extra loop-type branch emerging from a 'single point' in the (σ /σ 1 , A, B) space. Finally, the stable nearly standing (square-like) waves are possible for all values 0 |δ 1 | 1 and this wave mode co-exists with others. This means that experimental detections of the nearly standing waves need specific initial scenarios, which allow to discriminate it from other co-existing stable wave mode.

Diagonal-type elliptic orbits
The square base tank translatory moves along an elliptic orbit whose major axis belongs to the tank diagonal as shown in figure 2(d). Taking γ = π/4 reduces (3.21) to the form Three-dimensional sloshing in a square. Part 5 894 A10-25 In contrast to the cases in § § 4 and 5, the secular system depends on the sign of δ 1 , or, in other words, the orbit direction (clockwise or counterclockwise) affects the waveamplitude response curves. On the other hand, the secular system (6.1) possesses the following special kind of symmetry: if (A 2 , B 2 , β; δ 1 , ±) is a solution of the secular system, then (B 2 , A 2 , β; −δ 1 , ∓) is also its solution. This means that, as long as we know the response curves for, e.g. the diagonal-type counterclockwise elliptic forcing, substitution A → B, B → A outputs the response curves for the diagonal-type clockwise elliptic forcing. Due to this symmetry (relative to the diagonal plane), the forthcoming analysis may, without loss of generality, focus on the counterclockwise (0 < δ 1 1) tank orbits. There are two limiting cases of the diagonal-type elliptic forcing, which are associated with δ 1 = 0 (the diagonal harmonic tank excitation was studied in Part 4) and δ 1 = 1 (the circular forcing was analysed in § 4).
A purpose of Part 4 has been to resolve contradictions between experimental measurements by Ikeda et al. (2012) and theoretical results from Part 1. The undamped analysis in Part 1 established two physically identical (to within the angular propagation direction) swirling waves, which have, because of the four symmetry planes (x = 0, y = 0, y = x and y = −x), equal amplitudes on the perpendicular walls, i.e. A = B in terms of the lowest-order approximation. However, experiments by Ikeda et al. (2012) did not confirm the latter equality. According to the experimental measurements, B > A for the counterclockwise swirling but A < B when the rotary wave propagates clockwise. Part 4 explains and quantifies these inequalities by introducing a non-zero damping into the Narimanov-Moiseev modal theory.
Results from Part 4 (figure 8b) on the steady-state wave-amplitude response curves in the (σ /σ 1 , A, B) space with projections on the (σ /σ 1 , A) and (σ /σ 1 , B) planes are reproduced (almost identically) in figure 7(a). The same input data (3.30) and notation are adopted. The wave-amplitude branching consists of the connected curve with two bifurcation points U 1 and W where the branch P l D 0 U 1 WP r laying in the A = B plane is linked with the two branches U 1 D 1 V 1 S 0 R 1 W and U 1 P 0 V 2 R 2 W. Points on the branch P l D 0 U 1 WP r imply diagonal (standing) waves (stable and unstable). The branch was drawn by using the analytical solution The branches U 1 D 1 V 1 S 0 R 1 W and U 1 P 0 V 2 R 2 W are computed with the secular system (6.1). Stable periodic solutions on these branches imply counterclockwise and clockwise swirling, respectively. Irregular waves are expected in the frequency range between U 1 and V 1 .
The steady-state wave-amplitude response curves for the diagonal reciprocating (γ = π/4, δ 1 = 0, (a), reproduction from Part 4) and diagonal-type elliptic (γ = π/4, δ 1 = 0.05, (b)) forcing. The physical input parameters are defined in (3.30). In the column (a), P l D 0 U 1 WP r implies the standing diagonal wave mode (stable and unstable) with A = B but U 1 D 1 V 1 S 0 R 1 W and U 1 P 0 V 2 R 2 W correspond to the counterclockwise and clockwise swirling waves, which, due to the non-zero damping, are characterised by different wave amplitudes at the perpendicular tank walls (the latter phenomenon was observed by Ikeda et al. (2012) and theoretically described in Part 4). Irregular sloshing is expected for the forcing frequencies between U 1 and V 1,2 . In the column (b), the non-zero but relatively small semi-axis ratio δ 1 = 0.05 causes a splitting of the continuous branching and yields the branch P 1 D 0 E E D 1 U 1 V 1 S 0 R 1 P r (stable counterclockwise swirling on V 1 S 0 , nearly standing waves on P l T, D 0 U 1 and R 1 P r , and standing (square-like) waves at E and E ) as well as the loop-type branch WS 0 V 2 P 0 W (stable clockwise swirling on V 2 S 0 ). The irregular waves are then expected between U 1 and V 1 . Figure 7(b) demonstrates what happens with the wave-amplitude branching in the column (a) when δ 1 is not zero but remains relatively small (δ 1 = 0.05 in these computations). The column (b) shows a splitting of the originally continuous branching in the column (a). The first branch P l D 0 E E D 1 U 1 V 1 S 0 R 1 P r is a continuous curve going through the entire resonant zone but the second one, P 0 WS 0 V 2 P 0 , has a loop-type shape. Intervals P l T, D 0 U 1 and WP r in the column (a) are responsible for the diagonal (standing) waves. They convert to the nearly standing (nearly diagonal) waves in the column (b) (except at points E and E implying the (purely) square-like sloshing). Interval V 1 S 0 on P l D 0 E E D 1 U 1 V 1 S 0 R 1 P r is responsible for the stable steady-state solutions, which describe the counterclockwise swirling; V 2 S 0 on the loop-type curve WS 0 V 2 P 0 W implies the stable clockwise swirling. The frequency range of irregular waves (determined by U 1 and V 1 ) saves approximately fixed location and width in (a) and (b).
When the semi-axis ratio δ 1 increases to become of the order O(1) on the adopted asymptotic scale (3.30), the response curves are a subject of metamorphoses, which are illustrated in figures 8 and 9. In figure 8, computations were done with δ 1 = 0.3 (a) and 0.5 (b). Both three-dimensional views and projections on the planes (σ /σ 1 , A) and (σ /σ 1 , B) are depicted. The figure establishes the following global trends: disappearance of the loop-type branch P 0 V 2 S 0 so that stable clockwise swirling becomes impossible (0.38 δ 1 for the given input data), increase of the frequency range of irregular waves (the lower bound is associated with T), appearance of the stable nearly standing wave mode (on the small 'island' U 1 U 1 ) and, finally, conversion of the nearly standing wave mode to counterclockwise swirling on P 1 T and R 1 P r . The wave-amplitude branching in figure 8(b) is typical for the semi-axis ratios 0.38 δ 1 0.75 (at least, for the input parameters (3.30)). The branching consists of the single response curve P l TD 0 E U 1 U 1 U 1 E D 1 V 1 S 0 R 1 P r . Stable and unstable zones periodically alternate each other along this curve. Hiramitsu & Funakoshi (2015) publishes two computational examples in the (forcing frequency, non-dimensional energy)-plane for the diagonal-type forcing corresponding to δ 1 = 0 and 0.364. The results are qualitatively consistent with our analysis. Figure 9 demonstrates geometric and topologic metamorphoses of the response curves for the semi-axis ratios 0.75 δ 1 < 1. An emphasis is placed on how the curve P l TD 0 E U 1 U 1 U 1 E D 1 V 1 S 0 R 1 P r breaks away at δ 1 = 0.785 to yield one continuous curve in the entire resonant zone plus a loop-type branch. Details of this transformation are explained in panels (a-c). Panel (a) shows how the stability 'island' U 1 U 1 (responsible for nearly standing waves) disappears but the stable counterclockwise swirling mode becomes possible for points on the interval V 1 V 1 . The branching still consists of a single continuous curve. A further increase of δ 1 to 0.785 causes intersection of V 1 V 1 and TD 0 , which yields bifurcation point V 0 (panel b). Parameter β is far from the unit about V 1 but β ≈ 0.93 at V 1 so that the stable steady-state sloshing changes its type along the small interval V 1 V 0 V 1from the nearly standing (square-like) wave mode at V 1 to the counterclockwise swirling at V 1 . Panel (c) demonstrates a splitting of the connected response curves into the branch P l V 0 V 1 V 1 S 0 R 1 P r (counterclockwise swirling) and the loop-type branch D 1 E U 1 E D 0 V 0 V 1 D 1 (nearly standing waves). For those values of δ 1 , irregular waves are only possible in a narrow frequency range between V 1 and V 1 . This wave type fully disappears for δ 1 = 0.95 in panel (d). The latter panel also contains the extra loop-type branchD 0Ẽ Ũ 1Ẽ D 0 , which emerges from a single point at δ 1 = 0.88.
In summary, passage from reciprocating to elliptic diagonal-type excitations causes qualitatively similar metamorphoses of the wave-amplitude response curves and The steady-state wave-amplitude response curves for the diagonal-type (γ = π/4) elliptic counterclockwise orbital forcing with the semi-axis ratios δ 1 = 0.3 (a) and δ 1 = 0.5 (b). Physical input parameters, notation and symbols are taken from figure 7(b). Points E and E correspond to the square-like (standing) waves but points on D 0 U 1 and (relatively small) interval U 1 U 1 imply the stable nearly standing wave mode. Solutions on P l T and R 1 P r convert to the stable counterclockwise swirling, which is also represented by points on V 1 S 0 . Irregular waves are possible between U 1 and V 1 as well as T and U 1 . The clockwise swirling is possible for δ 1 = 0.3 but the semi-axis ratio = 0.5 makes this wave mode impossible.
associated steady-state wave regimes to those described in § 5 for the wall-symmetric position of the elliptic orbit. This includes conversion of stable (diagonal) standing waves to the swirling wave mode (co-directed with the orbit direction), a narrowing and vanishing frequency range for irregular waves, disappearance of the counterdirected swirling as well as a constantly existing zone of the stable nearly standing wave mode. On the other hand, original (δ 1 = 0) and intermediate (0 < δ 1 < 1) shapes of the response curves during their metamorphoses with increasing the semi-axis ratio differs from those in § 5.  indicates the vanishing stability 'island' U 1 U 1 (nearly standing waves) and appearance of the 'island' V 1 V 1 (stable counterclockwise swirling). Zone of irregular waves does not increase. Panels (b) and (c) show an emerge of one from two loop-type branches, which are expected in the limit δ 1 → 1. For δ 1 = 0.785 in (b), the 'island' V 1 V 1 touches TD 0 to constitute bifurcation point V 0 . A further increase of δ 1 (panel (c), δ 1 = 0.786 in computations) divides the branching at this bifurcation point and leads to the continuous response curve P l TV 0 V 0 V 1 S 0 R 1 P r (counterclockwise swirling) and the loop-type branch D 1 E U 1 E D 0 V 0 D 1 (the nearly standing sloshing except at E where a (purely) standing wave is expected). Panel (d) shows that the branching becomes close to that in figure 4 as δ 1 tends to the unit (the circular tank orbit).

Oblique-type elliptic tank orbits
After examining longitudinal (γ = 0, § 5) and diagonal (γ = π/4, § 6) positions of elliptic orbits and clarifying what are typical transformations of the wave-amplitude response curves versus the semi-axis ratio 0 |δ 1 | 1, one can finally focus on the most general angular position of the elliptic orbit in figure 2(a), which will be called the oblique-type elliptic forcing. This case turned out to be much more complicated, requiring us to consider, independently, co-and counter-directed elliptic orbits and, moreover, we have to present results for two trial 'oblique' angles γ = π/6 and π/12. The results are reported in figures 10-16.
Let us first focus on mathematical and physical aspects and admissible solutions of the considered secular equations (3.21). By choosing an appropriate position of the Oxy frame and, thereby, using symmetry relative to the diagonal planes, the angle γ in figure 2(a) may be assumed belonging to the interval −π/4 γ −π/4. Furthermore, Υ 2 y and, therefore, EQ1 ± and EQ2 ± in (3.21a) and (3.21b), respectively, are independent of the signs of γ and δ 1 . However, analysing EQ3 ± in (3.21c) shows that the wave-amplitude response curves and associated steady-state wave modes are different for (γ δ 1 ) > 0 (the case implies either counterclockwise trajectory for positive γ or clockwise trajectory with negative γ ) and (γ δ 1 ) < 0 (either clockwise trajectory for positive γ or counter-clockwise trajectory with negative γ ). This property is a consequence of the symmetry relative to the vertical coordinate planes Oxz and Oyz, which keeps invariant the sign of (γ δ 1 ) (just consider the Ox symmetric ellipses in figure 2a for different orbital directions!). There are no other symmetries of the square, thus, in order to describe all possible steady-state wave motions for the oblique-type elliptic forcing, one must independently analyse (γ δ 1 ) > 0 and (γ δ 1 ) < 0. In the present section, we adopt positive 0 < γ < π/4 but consider counterclockwise (0 < δ 1 < 1) and clockwise (−1 < δ 1 < 0) orbital tank motions.
The authors were not able to construct any particular analytical solutions of the secular system (3.21) for the oblique-type tank forcing and, therefore, the forthcoming analysis will be done by utilising the computational scheme in § 3. Numerical examples will adopt the trial angle γ = π/6 and π/12 and the input data (3.30). A requirement to consider the two different trial angles and build-up many illustrative graphs in figures 10-16 is caused, as our numerical experiments showed, by significantly different metamorphoses (with changing 0 < |δ 1 | < 1) of the response curves for 0 < γ π/8 and π/8 γ < π/4. Restricting to one trial 0 < γ < π/4 would make the steady-state analysis incomplete. Figure 10(a) reproduces results from Part 4 on the wave-amplitude response curves for the oblique harmonic (reciprocating) tank forcing with γ = π/6. The figure keeps unchanged basic notations for branches and bifurcation points by Faltinsen & Timokha (2017, figure 6b). The solid lines imply stable solutions (wave regimes) but the dashed (deep-blue) lines mark unstable ones. For that case, applicability of the Narimanov-Moiseev-type modal theory and the corresponding asymptotic steady-state analysis were supported by observations and measurements in Ikeda et al. (2012). The experiments also confirmed different maximum wave elevations at perpendicular walls (inequality A = B) for counterclockwise and clockwise swirling with the same forcing frequency σ /σ 1 .

Counterclockwise elliptic forcing
The wave-amplitude branching in the column (a) contains the continuous response curve P l TD 0 UVS 0 WW 0 P r running through the entire resonant zone and the loop-type curve P 0 S 0 V P 0 existing only in the primary resonant zone. The counterclockwise stable swirling (points on VS 0 and WW 0 ) is associated with the aforementioned continuous response curve but the clockwise stable swirling (V S 0 ) -with the loop-type branch P 0 S 0 V P 0 . All other stable steady-state waves are of the nearly standing type (points on P l T, W 0 P r and D 0 U). An irregular sloshing is expected in the frequency range between abscissas of U and V. The response curves do not contain points responsible for the (purely) standing wave mode (β = 1 or sin α = 0) even away from the primary resonant zone on P l T and W 0 P r where the linear sloshing theory could be applicable. Part 4 theoretically clarified this and other specific peculiarities by a non-zero damping in the hydrodynamic system. Indeed, when ξ = 0, points on P l T and W 0 P r imply the square-like standing wave mode.
Attention should be paid to the stability interval WW 0 P r . When approaching P r (a point away from the primary resonance zone), β becomes very close but not equal to unity (β ≈ 0.9999 for the input data (3.30)), namely, as we mentioned above, the nearly standing sloshing occurs. However, when going from P r to W, β rapidly decreases to achieve O( 1/3 ) 1 − β in the vicinity of W 0 and, later on, 1 − β = O(1), sin α > 0 in a local neighbourhood of W. This means that the nearly standing waves for points on P r W 0 change to a counterclockwise stable swirling on W 0 W.
FIGURE 10. The steady-state wave-amplitude response curves for the oblique-type elliptic forcing with γ = π/6 and counterclockwise tank orbits with the semi-axis ratio δ 1 = 0 (column a corresponds to the limiting case -the oblique harmonic reciprocation) and 0.45 (column b). The input data from (3.30). Notations are taken from Part 4, the solid lines imply the stable wave regimes but the dashed (deep-blue) lines mark the waves instability. The column (a) is a reproduction (with minor changes) of figure 6(b) from Part 4; it shows that the oblique reciprocation with positive 0 < γ < π/4 yields the continuous branch P l TD 0 UVS 0 WW 0 P r implying the counterclockwise stable swirling (on VS 0 and WW 0 ) in the primary resonant zone and the nearly standing waves away from the resonance zone (on P l T, W 0 P r and D 0 U). The stable clockwise stable swirling corresponds to points on V S 0 belonging to the loop-type branch P 0 S 0 V P 0 . Irregular waves are expected in the frequency range UV. Passage to positive 0.45 = δ 1 = O(1) in the column (b) saves P l TD 0 UVS 0 WW 0 P r but the clockwise stable swirling on P 0 S 0 V P 0 disappears. The column (b) shows that, when δ 1 = O(1) on the chosen asymptotic scale, the pieces P l T and W 0 P r become responsible for the stable counterclockwise swirling.
Passage from zero to positive but relatively small value of δ 1 keeps a topologically identical branching. As in figure 10(a), the branching consists of the continuous response curve P l TD 0 UVS 0 WW 0 P r and the loop-type branch P 0 S 0 V P 0 . However, the loop-type branch reduces with increasing δ 1 and, finally, disappears for the semi-axis ratio exceeding a critical value (≈ 0.5 for the input data (3.30)). When approaching this critical value from above, the loop-type branch P 0 S 0 V P 0 gets rid of points responsible for the stable clockwise swirling. In the present computations, this happens approximately at δ 1 = 0.4. Figure 10(b) demonstrates the wave-amplitude response curves for the oblique-type elliptic counterclockwise forcing with δ 1 = 0.45. They confirm the aforementioned feature. Comparison of the columns (a) and (b) in figure 10 shows that points on P l T and W(W 0 )P r become corresponding to the stable counterclockwise swirling (here, co-directed with the orbit) as it was observed before for other positions of the elliptic orbit when δ 1 = O(1). Irregular waves are still predicted in the frequency range defined by T and V. Furthermore, the column (b) demonstrates a decrease of D 0 U (nearly standing waves), an increase of VS 0 (counterclockwise swirling) and a drift of W to the right. Figure 11(a) shows that the loop-type branch indeed disappears for δ 1 = 0.5 (as in figure 8b for the diagonal-type forcing) and the branch P l TD 0 UVS 0 WP r keeps invariant its structure and consequence of the stable/unstable intervals. The clockwise (counterdirected to the tank orbit) stable swirling is impossible but the counterclockwise (codirected with the orbital forcing) swirling occurs in the frequency ranges, which are similar to those in figure 10(b). Figure 11(b-f ) demonstrate metamorphoses of the wave-amplitude response curves when 0.75 δ 1 0.95 (the cases are similar to those in figure 9 but for different γ ). The branching should reach the geometric shape in figure 4 as δ 1 approaches the unit. Comparing panels (a) and (b) (δ 1 = 0.5 and 0.75, respectively) shows generally negligible geometric changes of the continuous curve P l D 0Ẽ UẼ VS 0 WP r ; the only difference is two pointsẼ andẼ , which imply the square-like (purely standing) wave mode. An astonishing thing is appearance of the small loop-type branch D 0 , which is born from a 'single point' at, approximately, δ 1 = 0.746.
The limiting case with δ 1 = 1 also requires the second loop-type branch. Three panels (c-e) illustrate how this branch (D 0Ẽ UṼ V 0 D 0 ) 'buds' from the continuous response curve P l D 0Ẽ UẼ VS 0 WP r . Computations in the panels are done with δ 1 = 0.796, 0.7972 and 0.798, respectively. Evolution of the branching is analogous to what we observed in figure 9(a-c). It has the same three compulsory stages. These are, first, a new (and relatively small) 'island' V (V 0 )V , which emerges on P l D 0Ẽ UẼ VS 0 WP r (panel c) so that β rapidly changes along it, causing the nearly standing wave mode in a neighbourhood of the endpoint V and the counterclockwise swirling in a neighbourhood of the endpoint V . Second, the piece V (V 0 )V touches TD 0 to constitute the bifurcation point V 0 (panel (d), δ 1 = 0.7972). Finally, the loop-type branch D 0Ẽ UṼ V 0 D 0 breaks away with a further increase of δ 1 as it is shown in panel (e). For 0.8 δ 1 1, the stable counterclockwise (co-directed with the forcing trajectory) swirling exists in the entire resonant zone. Two stable nearly standing waves co-exist with the stable swirling in a certain frequency domain. Irregular sloshing (all steady-state regimes are unstable) is not predicted. Figure 12 demonstrates changes of the response curves when γ differs from the trial value π/6 (γ = π/12 in these calculations). Qualitatively, the wave-amplitude branching is similar to those in figures 10 and 11. The principal difference is the loop-type branch P 0 S 0 V U D 0 P 0 (P 0 D 0 U P 0 ). For γ = π/6, this branch vanished at δ 1 = 0.45 but, later on, it emerged from a single point at δ 1 = 0.746. Figure 12 shows that this branch exists (does not vanish for certain values of δ 1 ) for all 0 < δ 1 < 1 as γ = π/12. Physically, this fact implies two stable nearly standing waves (D 0 U and D 0 U ) in figure 12 for arbitrary semi-axis ratios.    figure 10(b) shows that the loop-type branch P 0 S 0 disappears for δ 1 = 0.5. In the range 0.47 δ 1 0.75, only stable counterclockwise swirling, nearly square-like waves and irregular sloshing are predicted. Panels (b-f ) show what happens for 0.75 δ 1 < 1 (the semi-axis ratio approaches the unit). A loop-type branch with B < A emerges from a 'single point' and becomes visible in (b). The branch grows with increasing δ 1 to yield a frequency range for the stable nearly standing waves. Another loop-type branch (A < B) is a consequence of touching V V (stable swirling at V and nearly standing wave at V ) and TD 0 (unstable sloshing). This kind of intersection and related bifurcation were described in figure 9 for the diagonal-type elliptic forcing. Irregular waves are not predicted in the cases (e) and ( f ), i.e. for the semi-axis ratios 0.8 δ 1 1. Hiramitsu & Funakoshi (2015) presented one numerical branching in the frequencyenergy plane for the counterclockwise orbit with γ = π/8 and δ 1 = 0.364. The results are qualitatively consistent with our calculations. They detect one continuous branch going through the entire resonant zone and a loop-type branch with an interval responsible for the counter-directed swirling mode.
curve P l D 0 UVS 0 WW 0 P r and the loop-type branch P 0 V S 0 V (U D 0 )P 0 . The first branch contains points responsible for the stable counterclockwise swirling (VS 0 ) and the second one has the interval V S 0 , which corresponds to the stable clockwise swirling mode. In other words, the oblique harmonic forcing cannot lead to physically identical but oppositely propagating swirling waves as happened for the longitudinal excitation.
To reach the limiting wave-amplitude branching in figure 4 as |δ 1 | → 1, the response curves should go through very specific metamorphoses, which are accompanied, in particular, by a removal of the stable counterclockwise from the P l D 0 UVS 0 WW 0 P r branch. How this may happen is demonstrated in the two illustrative examples for γ = π/6 and π/12 with −1 < δ 1 < 0. Let us consider the oblique harmonic tank reciprocation with the angle γ = π/6 and the corresponding wave-amplitude response curves in figure 10(a). After that, let us slightly perturb δ 1 < 0, which physically means clockwise elliptic tank motions with a small semi-axis ratio. Figure 13(a) presents the first interesting result for δ 1 = −0.015. The response curves for this semi-axis ratio look indeed rather interesting. One can see that the two originally (when δ 1 = 0) non-connected branches P l D 0 UVS 0 WW 0 P r and P 0 S 0 UD 0 touch each other and yield the intersection (bifurcation) point W 0 .
FIGURE 13. The steady-state wave-amplitude response curves caused by the oblique-type clockwise elliptic forcing with γ = π/6. Here, δ 1 = −0.015 (a) and δ 1 = −0.025 (b). The input data are taken from (3.30) but notation was introduced in figure 10(a) (oblique harmonic reciprocation, δ 1 = 0). The two columns demonstrate qualitative changes in the wave-amplitude branching when the tank orbit has a relatively small semi-axis ratio. The reciprocating tank forcing causes two non-connected branches in figure 10(a). The column (a) shows that these response curves join at bifurcation point W 0 , which implies the standing (square-like) wave mode. This happens about δ 1 = −0.015. Another standing wave is associated with point E, which coincides with the turning point T. When the semi-axis ratio continues increasing, this branching breaks away at W 0 (the bifurcation point 'jumps' to P 0 S 0 ). The column (b) (δ 1 = −0.025) demonstrates that the branching consists then of the single continuous curve P l TED 0 UVS 0 W 0 P 0 V S 0 W 0 P r . Disposition of frequency ranges for stable swirling, nearly standing and irregular waves is practically identical to that in figure 10(a) except for a narrow zone W 0 W , which defines the stable clockwise (co-directed with the forcing direction) swirling.
Specifically, the bifurcation point W 0 implies the standing (square-like) wave (β = 1, see figure 3). There also exists the second standing wave, which is associated with E, which coincides with the turning point T. Because the semi-axis ratio |δ 1 | = 0.015 is O( ) on the asymptotic scale (3.30), the frequency ranges of stable steady-state waves remain almost identical to those in figure 10(a) where δ 1 = 0. An exception is the piece (interval) W 0 W , which emerges from the bifurcation point W 0 . Points on W 0 W imply the stable nearly standing waves. Figure 13(b) illustrates modifications of the wave-amplitude response curves due to a further increase of the semi-axis ratio (computations were done with δ 1 = −0.025). After disconnecting at W 0 , the wave-amplitude branching becomes the single continuous curve P l TED 0 UVS 0 W 0 P 0 V S 0 W 0 P r where E and W 0 remain responsible for the unstable (purely) standing waves. Specifically, the point E moves then along this curve from 'periphery' to the primary resonance zone but W 0 'jumps' onto the interval P 0 S 0 . The semi-axis ratio |δ 1 | = 0.025 remains O( ) and, therefore, the steady-state wave amplitudes and the frequency ranges of the stable steady-state wave modes (see projections on (σ /σ 1 , A) and (σ /σ 1 , B)) negligibly change in the column (b). However, there are qualitative differences in the vicinity of the former bifurcation point W 0 . In panel (a), W 0 W corresponds to stable counterclockwise swirling and W 0 W implies stable nearly standing waves. Because W 0 jumps onto P 0 S 0 in panel (b), the interval W 0 W vanishes but, because β rapidly decreases along W 0 W , the stable clockwise (co-directed with the tank orbit) swirling wave is predicted in a local neighbourhood of W .
A further metamorphosis of the wave-amplitude response curves for the oblique-type clockwise elliptic forcing with increasing |δ 1 | = O(1) on the asymptotic scale (3.30) is demonstrated in figures 14 and 15. Figures 13(b, δ 1 = −0.025) and 14(a, δ 1 = −0.275) display the topologically equivalent branching, which consists of a single continuous curve running from P l to P r though the primary resonant zone. However, figure 14(a) exhibits the following differences: (i) points E and W 0 continuously drift into the primary resonance zone along the response curve, (ii) points on P l T and P r W become responsible for stable clockwise swirling (co-directed with the elliptic tank orbits), (iii) the interval VS 0 , on which stable counterclockwise swirling occurs, diminishes but V S 0 (clockwise, co-directed swirling) enlarges, (iv) a new interval U U appears, along which β rapidly changes to cause the nearly standing wave mode at U and counterclockwise swirling in the vicinity of U .
An increase of the semi-axis ratio |δ 1 | in figure 14(b-d) 'extrudes' the stable counterclockwise (counter-directed) swirling mode from the resonant zone. The associated topological metamorphoses look rather attractive. First, the branch piece U U touches P 0 S 0 and constitutes bifurcation point V as shown in panel (b) (δ 1 = −0.28). Second, a further increase of the semi-axis ratio causes a break away of the loop V S 0 VV (panel c) accompanied by a rapid reduction of VS 0 (stable counterclockwise swirling) with increasing |δ 1 |. The loop-type branch vanishes at δ 1 ≈ −0.3 (for the given input data). Panel (d, δ 1 = −0.5) demonstrates the branching, which is similar to that in figure 11(a). An exception is existence of the (purely) standing waves associated with E and W 0 . Similar disposition and sizes of the effective (for stable steady-state waves) frequency ranges and irregular waves are expected for both directions (clockwise and counterclockwise) of the oblique-type elliptic tank orbits when the semi-axis ratio is close to ≈0.5. Figure 15 examines the wave-amplitude response curves and associated steady-state waves when the clockwise oblique-type elliptic forcing trajectory tends to the circular shape. The calculations were done with δ 1 = −0.702 (a), δ 1 = −0.7029 (b), FIGURE 14. The steady-state wave-amplitude response curves caused by the oblique-type clockwise elliptic forcing for γ = π/6 with δ 1 = −0.275 (a), δ 1 = −0.280 (b), δ 1 = −0.290 (c) and δ 1 = −0.5 (d). The input data are defined in (3.30). Basic notation is taken from figure 13. The figure demonstrates a diminishing of the counterclockwise (counter-directed to the elliptic forcing) swirling through forming the bifurcation point V in (b), a further break away of the loop-type branch VS 0 (b,c) and, finally, disappearance of the latter branch (d).  δ 1 = −0.704 (c) and δ 1 = −0.95 (d) and γ = π/6. Formation of the loop-type branch D 0 EUW 0 D 0 (B > A) is similar to that in figure 11(c-e). In both cases, the co-directed stable swirling mode exists in the entire resonance zone (irregular waves become impossible) starting with the semi-axis ratio |δ 1 | = 0.75.    Figure 16 shows the wave-amplitude response curves for the oblique-type clockwise elliptic forcing when γ = π/12. The starting situation is depicted in figure 12(a), which corresponds to the harmonic oblique reciprocation (δ 1 = 0). Figure 16(a) demonstrates how the two non-connected branches in figure 12(a) join each other and constitute bifurcation point W 0 (as in figure 13b). Figure 16(b,c) demonstrate formation of a single response curve whose stable and unstable intervals alternate with each other. As in the previous examples, the counterclockwise swirling disappears for the semi-axis ratio ≈0.5. Indeed, the corresponding interval S 0 V is relatively small in cases (c-e) and it fully vanishes as 0.5 |δ 1 |. Figure 16(d,e) shows, first, intersection of the response curve at point B with the forthcoming break of at the bifurcation point, which leads to the loop-type branch W 0 UED 0 W 0 . Another loop-type branch W 0ŨẼD0W0 emerges from a single point when |δ 1 | tends to the unit.
In summary, the oblique-type elliptic forcing causes a large variety of the wave-amplitude response curves whose topology and geometry strongly depend on the semi-axis ratio and, In particular, on the tank orbit direction. Metamorphoses of the wave-amplitude response curves are very attractive when (γ δ 1 ) < 0 (either clockwise trajectory for positive γ or counter-clockwise trajectory with negative γ ). In particular, when 0 < |δ 1 | 0.5 (the upper bound depends on the physical input). When |δ 1 | ≈ 0.5, the counter-directed (to the tank orbit) swirling, normally, disappears. Furthermore, disposition and frequency ranges of the stable nearly standing and co-directed swirling waves as well as irregular waves remain similar for the semi-axis ratios 0.5 |δ 1 | 0.75. A further increase of the semi-axis ratio (0.75 |δ 1 |) transforms the wave-amplitude response curves to those for the circular orbital forcing.

Conclusions
Bearing in mind a systematic study of the resonant nonlinear sloshing in a square base container, steady-state and transient waves, the authors derived in Part 1 a (modal) system of nonlinear ordinary differential equations within the framework of the Narimanov-Moiseev asymptotic theory. The theory is valid for non-parametric tank forcing under special conditions to the forcing amplitude and frequency. The derived Narimanov-Moiseev-type modal system suggested a resonant excitation of the two lowest (degenerating) natural sloshing modes, a finite liquid depth and periodic tank motions by sway, surge, roll, pitch and, generally, yaw. The system was extensively exploited in the forthcoming parts to study the steady-state resonant wave regimes for the longitudinal (along parallel walls), diagonal and oblique reciprocations of the tank. Focusing exclusively on the reciprocating excitations was partly driven by existing (including the authors') experimental observations and measurements (all these have supported the theoretical results) but also it was caused by a series of mathematical problems in getting an analytical solution of the so-called secular system. Experiments by Ikeda et al. (2012) and the forthcoming theoretical investigations by Faltinsen & Timokha (2017, Part 4) showed that damping, even being relatively small, is important to explain some resonant wave paradoxes. Incorporating the linear damping terms in the Narimanov-Moiseev-type modal system made it impossible to get a general analytical solution of the secular system for the reciprocating forcing, more precisely, the authors were able to find analytical solutions only in very particular cases. They had to develop a numerical method in their studies. After introducing several modifications, the numerical method is used in the present paper to classify the steady-state sloshing due to the three-dimensional non-parametric cyclic tank motion (combined periodic sway, surge, roll, pitch and yaw). As matter of the fact, the present Part 5 realises the main goal of the case studies, fully revealing the potential of the Narimanov-Moiseev-type modal theory from Part 1.
The classification of the steady-state wave regimes became possible because, as we prove in the present paper, there exists an asymptotic equivalence between the steady-state resonant waves due to sway, surge, roll, pitch and yaw and those excited by a suitable horizontal (translatory) elliptic orbital tank motion. The asymptotically equivalent waves are characterised by the same first-and second-order asymptotic components (within the framework of the Narimanov-Moiseev asymptotic theory) and possess the same stability properties. This means that, instead of making an exhausting search of all possible periodic non-parametric excitation scenarios, one can investigate the steady-state wave regimes and their stability versus the three input parameters: 0 γ π/4 (characterises angular position of the elliptic orbit), 0 |δ 1 | 1 (the semi-axis ratio of the orbit) and the orbit direction (clockwise or counterclockwise). In view of this fact, the previous paper parts concentrated on the limiting case δ 1 = 0. The present paper independently examines another limiting case |δ 1 | = 1 (the circular tank forcing) and, furthermore, concentrates on the particular cases with canonic (wallsymmetric, γ = 0, 0 < |δ 1 | < 1) and diagonal (γ = π/4, 0 < |δ 1 | < 1) positions of the elliptic orbits. Finally, the oblique-type elliptic forcing is numerically analysed for γ = π/12 and π/6 with −1 < δ 1 < 1 where positive δ 1 denotes counterclockwise orbits but negative δ 1 -clockwise ones. Because analytical solutions were found only in very particular cases, the steady-state wave-amplitude analysis is done, numerically, with the input parameters (3.30), for which Part 4 showed a good agreement between the Narimanov-Moiseev-type theoretical predictions and experiments by Ikeda et al. (2012). Within the framework of the lowest-order approximation, the resonant steady-state waves fall apart into the purely and nearly standing wave modes and the swirling wave mode, which may occur either clockwise or counterclockwise, as well as an irregular (chaotic, modulated, etc.) sloshing. Introducing A and B, which measure the lowest-order wave-amplitude components in the Ox and Oy directions, respectively, and the corresponding phase lags, ψ and ϕ, makes it possible to systematically analyse the wave-amplitude response curves in the (σ /σ 1 , A, B) space specifying, in a parallel way, stable and unstable wave modes for each point on these curves. The latter uses a criterion based on the sin(ϕ − ψ)-values (proven in the present paper) and the Lyapunov-type stability analysis from Part 4.
The circular orbital tank forcing (|δ 1 | = 1) excites the stable resonant steady-state swirling whose angular propagation (around Oz) coincides with the forcing direction (clockwise or counterclockwise). This wave mode exists for all resonant forcing frequencies. Irregular waves are not predicted. The swirling wave amplitude is governed by the damped hard-spring-type response curves and, therefore, two different swirling waves may coexist in a certain frequency range where the forcing frequency is slightly larger than the lowest natural sloshing frequency. A novelty with regard to the same damped hard-spring type response curves for the orbitally excited swirling in an upright circular tank (Raynovskyy & Timokha 2018b) consists of two stable nearly standing waves occurring in a relatively narrow frequency range for the forcing frequency, which is slightly lower than the lowest natural sloshing frequency. The nearly standing waves co-exist with swirling. They are characterised by a relatively large amplitude, that is, their occurrence requires special initial scenarios, which, in our opinion, are difficult to organise in model tests.
Passage from the harmonic longitudinal reciprocation to elliptic tank orbits (γ = 0, 0 |δ 1 | 1) breaks away the Ox symmetry of nearly standing and swirling waves so that, e.g., co-and counter-directed swirling modes have different wave amplitudes. There is a critical value of |δ 1 | (≈ 0.5 for the tested input data) when the counter-directed swirling disappears (the same happened for upright circular base containers in Raynovskyy & Timokha (2018a)). A further increase of the semi-axis ratio makes the co-directed swirling stable in the entire resonant zone and, as a consequence, irregular sloshing may not occur when the tank orbit tends to a circle. The stable nearly standing wave mode is realisable for all |δ 1 |; it always co-exists with other stable steady-state waves. As we noted above, this co-existence makes difficult to detect this wave mode in experimental studies.
When the diagonal harmonic excitation acquires its non-small elliptic component (γ = π/4, |δ 1 | = O(1)), the wave-amplitude response curves and associated steady-state wave regimes become similar to those discovered for the wall-symmetric position of the elliptic orbit. This includes transformation of stable (diagonal) standing waves to the swirling wave mode (co-directed with orbits), narrowing and vanishing of the frequency range for irregular waves, disappearance of the counter-directed swirling as well as a constantly existing zone of the stable nearly standing wave mode.
The oblique-type elliptic forcing (0 < γ < π/4, −1 < δ 1 < 1, the δ 1 sign defines the orbit direction) leads to a variety of astonishing response curves and bifurcations, which strongly depend on δ 1 . However, when the semi-axis ratio is O(1) (on the adopted asymptotic scale), the counter-directed swirling disappears starting with a critical value (0.5 |δ 1 | in our calculations). The same happened for γ = 0 and γ = π/4. Disposition and frequency ranges of stable nearly standing and co-directed swirling wave modes as well as irregular waves are also similar.
The present paper exhausts potential of the Narimanov-Moiseev-type modal equations in the steady-state sloshing analysis for a (nearly) square base tank and, thereby, finalises the case studies originated in Part 1. Extension to other forcing types, e.g. with other resonant frequencies or by heave as well as for shallow liquid depth, requires their serious modification. However, the authors believe, the present wave classification results exhibit a useful 'forecast' of steady-state wave types and their effective frequency ranges for a broad variety of three-dimensional resonant tank motions and, thereby, it could both facilitate engineers who face these tank motions in everyday practice to model, for instance, a floating ship tank in oblique sea and, even more importantly, motivate mechanicians and mathematicians to pay more attention to the corresponding experimental and numerical studies. The present paper could also be a guidance for future experimental studies, which should focus on classification of the steady-state wave, i.e. establishing the frequency ranges where wave modes from figure 3 are stable. The classification may become rather difficult when two or more stable wave modes co-exist in a frequency range. This problem was already discussed in § 4 in the context of the nearly standing waves on the branches U E andŨ Ẽ and stable swirling. Their experimental detection needs a new strategy with non-zero initial scenarios.