Novel Hybrid Aquatic-Aerial Vehicle to Survey in High Sea States: Initial Flow Dynamics on Dive and Breach

1. Introduction and Literature

From Leonardo Da Vinci’s concept of the ornithopter Primrose [1], to the design of Japan’s beak-nosed bullet train, biomimicry has been used to solve a host of problems and curiosities that would potentially continue to remain unresolved without natural inspiration to draw upon. One problem humanity is beginning to contend with is that single-use items are less of a commodity and more of a burden, and that future generations will require environmentally friendly technology to thrive. Hybrid Vehicle technology that can operate across multiple domains or environments [2-4] may offer some measure of a solution, while biomimicry may show how to make it a reality. Take an oceanographic survey ship attempting to study the properties of an ocean column; multiples of torpedo probes are potentially carried aboard to capture the thermoclines and haloclines, with few expected to return to the surface. Alternatively, consider a patrol aircraft attempting to locate with sonar a missing submarine, launching dozens of disposable sonobuoys without guaranteeing success. In the age of multi-billion transistor microchips, nature may still inspire a solution to designing and creating items to be reused or applied across multiple domains.

Numerous species of marine birds have evolved to develop the ability to dive from a cruise condition into water to catch fish below the wave action of an ocean surface. In the case of the Northern Gannet, the velocities at which they pierce the water’s surface can often exceed 50 m/s [5] and research of this Gannet action has inspired Hybrid Aquatic-Aerial Vehicle (HAAV) design [6-9]. Meanwhile, Emperor Penguins can dive from ice sheets of great heights before leaping up out of the water, breaching at short spurt speeds around 2-4 m/s as well as impressive ascent rates of up to 120 m per minute from depths as great as 60 m [10, 11]. Both Avians are trans-median, meaning their anatomies have evolved to allow them to operate in both air and water, admittedly to vastly different degrees - while the Gannet cannot swim, nor can the Penguin fly. This trans-median behaviour, the ability to comfortably operate above and below water while comfortably transitioning between the two, presents a novel method for solving several issues. Take the previous sonar search examples; if either were able to recover their sensory equipment autonomously, neither would have to carry as many units, ultimately saving space, consumables, and money. Moreover, if the scientific equipment could fly, dive, swim, breach, and fly again, then the operational coverage of either system would be significantly expanded, reducing the need for an entire ship for weeks or an aircraft for days. By observing the behaviours and abilities of two different bird species, we can conceive of a method to expand the capabilities of relevant materials and equipment to use them repeatedly, ultimately saving space and expenditure.

The concept of autonomous aerial vehicles able to land or dive into water, transit submerged below the waves, before breaching the water to take flight again – is not new. However, most studies involving prototypes and field trials, have occurred only in the last few years [3]. Such developments are primarily due to commercial electric-powered drone developments, including uncrewed underwater vehicles, supported by high-fidelity two-phase computational fluid dynamics (CFD) modelling [4], that is now achieving a consensus in the best approach (see methodology in section 2). Many of the novel designs arguably derive from major power tension [12], driving demand for autonomous robots across domains [13], including through commercialization into non-state actors [14]. Most of the published experimental and field work with novel designs has been the Chinese-based research into the Nezha HAAV demonstrator designs [15-19]. Narayanan, et al. [20] similarly highlighted numerous designs of trans-media vehicles, ranging from quad-copters that land and take off from the water’s surface to rocket-propelled vehicles that splash down into the water on entry.

The comprehensive literature review on hybrid aerial-aquatic vehicle (HAAV) designs and technologies by Yao, et al. [3] finds that the variable sweep designs are ‘the only HAAV that performs well in all four phases of motion, namely aerial flight, underwater cruising, water-entry and water-exit.’ Further, they find these have high load capacity, the greatest mission potential, and that the ‘sweptback-wing appears to be the more appropriate structure for future study.’ Hence, due to its morphing wings and propulsion, one of the most promising developments in HAAV design is the Nezha-SeaDart tail-sitting delta design shown in Error! Reference source not found.. This Nezha HAAV design underwent record-setting field trials at sea; however, due to its sub-metre size and its open propeller design and no sea state details, it has likely only achieved a sea state of one or 0.1 m wave height. These sea state trials were likely given confidence by the numerical and control research of Bi, et al. [21], into the effect of Sea State One waves on quad-copter propellers in a HAAV, attempting to transition. Notwithstanding such progress, Jin, et al. [18] realistically note, ‘Current HAUV prototypes cannot fully cope with the unstructured, harsh, and hazardous ocean environments considering the high-pressure and the strong hydrodynamic perturbations.’ Research into viable HAAV designs for sea states greater than one, especially with typical sea state periods of around eight to 11 seconds, is a gap likely requiring greater size, transition speed, protection for propulsion systems, and research funding than seen so far.

Figure 1. Nezha SeaDart design and trials (reproduced from [19]).

Figure 1. Nezha SeaDart design and trials (reproduced from [19]).

Our research on conceptual HAAV design began, like the earlier Nezha designs, with a surface landing and sinking — in our case, a submersible seaplane [4, 22, 23]. The conceptual designs made us quickly realize the transition was critically dependent on sea state. A second design effort in 2023 focused on plunge diving and accelerated breaching using sweeping wings and retractable propulsion, similar to the Nezha-SeaDart. A new addition is shown in Error! Reference source not found. whereby the HAAV would operate as a ‘child aircraft’ from the parent maritime patrol aircraft for electrical charging, range extension, and superior signal intelligence [24].

Figure 2. Recoverable parent-child HAAV for sonar search [24].

Figure 2. Recoverable parent-child HAAV for sonar search [24].

Three design features of the submersible seaplane were retained in the new HAAV design that is shown in Error! Reference source not found., with dimensions in Table 1. Two of the retained features come from Australian industry involvement in deep-sea research vessels. The first retained feature is using a proprietary syntactic foam, at least for the fuselage, for lightweight compressible strength at depth, but with good tensile strength compared to other commercial syntactic foams. The second retained feature is the use of a fast buoyancy engine for rapid ascent from depth to achieve a viable breach above high-sea-state waves. Finally, the third retained feature is electric ducted fans (EDF) as the proposed single source of propulsion for both air and subsea, retracted into the fuselage quickly before trans-media travel. EDF propulsion will likely be as limited as open propellers underwater compared to above water; however, according to [25], ‘EDFs have a series of advantages compared to open propellers due to their geometry and shroud including a physical safety barrier [22,23], thrust vectoring [22], reduced noise emissions [24,25], smaller size improving ground clearance [18] and superior static performance [26].’ The parent-child relationship was introduced to alleviate the autonomy demands of such a mission in the open ocean, and it is related to military research on swarming and ethical control of robotic and autonomous systems. Such parent-child aircraft use was envisaged last decade by researchers like Fresconi and Fermen-Coker [26] but only recently enabled by the safe reel-in technology of Sparrowhawk [27] and, more recently, X-61A Gremlin. Other parent-child docking possibilities for HAAVs can be surface-ship, subsea station, or land-based in archipelagic waters.

Figure 3. 3D printed vehicle models showing the penetration/breach (left) and glide/cruise (right) configurations, with key part locations identified.

Figure 3. 3D printed vehicle models showing the penetration/breach (left) and glide/cruise (right) configurations, with key part locations identified.

The Nezha HAAV designs have illustrated that two-phase CFD methods are accurate predictors for trans-media transition using both experimental methods and field trials [15-19]. Hence, our new HAAV design will initially be modelled and simulated based on these increasingly consistent methods. Qi, et al. [28], examined the forces experienced by an air-launched Autonomous Underwater Vehicle on initial entry into the water. Leveraging from similar research, Ma, et al. [29] looked into both the water entry and exit dynamics of a slender, axisymmetric vehicle, finding very good agreement between CFD modelling and simulation at low speeds of around 20-30 m/s. Such accuracy leverages years of research into the most appropriate CFD methods for aquatic-aerial trans-media modelling. A quantum leap in the accuracy of such CFD analysis for propeller-driven aircraft to transit water, such as quad-copters of the Nezha design, was the research by Wei, et al. [30] to appropriately account for each rotor in the two-phase medium. Building on that design-enabling numerical research base is perhaps exemplified best by Sun, et al. [31] who found numerical-experimental agreement in a slender HAAV during water exit. This work was extended by research still under review [32] where the HAAV is driven by an aft water propeller and a nose-mounted air propeller with careful scheduling through the transition stages to avoid propulsion overspeed and to achieve control. The numerical accuracy in the presence of multiple propellers, coupled with the aforementioned wave effect modelling [21], shows significant promise to digitally engineer HAAV designs much more before their prototyping, especially for sea state vulnerabilities and control.

Wu, et al. [33] researched numerically a swing-wing vehicles’ water entry characteristics under various entry conditions where the HAAV was 2 m in length, diving at 10 m/s into waves about Sea State Three (1 m high, 10 m apart). They mainly examined the effect of waves on the peak impact forces rather than the trajectory disturbance, and they did not vary the wave period or sea state. Similarly, Dong, et al. [34] researched numerically the plunge dive of a HAAV that was 1.62 m long with a swing-wing, nominally operating from a parent ship. They varied the dive speed from 2 m/s to 12 m/s and the angle of inclination in the dive, again examining the impact forces. Neither of these water entry researchers examined the penetration depth or time, which are critical for a HAAV buoyancy engine to establish underwater operation.

Liu, et al. [35] made an immersive numerical analysis of a HAAV’s exit behaviours, using various initial angles and relatively high breaching velocities around 13 m/s to 18 m/s. Their studies showed that wave action on the body of an existing axisymmetric vehicle significantly influenced its breaching angle, varying it by up to 25 deg. They were able to establish ‘a mathematical model for the motion stability of submersible aerial vehicles at the water–air interface’ … ‘so that the submersible aerial vehicle’s pitch angle and velocity at the end of vehicle’s water-exit process, corresponding to any initial inclination angle and velocity, can be solved.’ The model had sufficient accuracy to aid the control of such transitions. Their research determined that the ‘minimum water-exit velocity of submarine aerial vehicles should be greater than 10.8 m/s’ for a ‘level 5 wave’, which is presumed to be Sea State Five (wave heights 2.5 - 4 m). Another very recent approach by Yun, et al. [36] was to experiment with slender HAAV shapes exiting the water to build a relationship between ‘the pop-up height of the cylinder … to the release depth.’ Such relationships help characterize the design depth and buoyancy exchange necessary to attain sufficient height to establish air-based propulsion for continued launch above sea state waves.

While most studies of water-entry and water-exit have been numerical and focused on HAAV transition, a recent study by Lu, et al. [37] sought to minimize the time of exposure during both transitions by equations of motion. This research focused on breach velocities around 40 m/s and dives constrained to be between 0 m/s and 80 m/s with an optimum trajectory impact velocity of around 32 m/s. These values are well beyond the viability of reusable sonar HAAVs and more appropriate to submarine and air-launched missiles; however, the approach could be used to reduce the time of exposure of HAAVs during transition to irregular ocean waves of the type modelled by Liang, et al. [38].

Taken together, all these studies have shown that the forces experienced by a small air vehicle transitioning at speed into or out of water may be tolerable for our current level of technology, which leads to our study. In this paper, we begin to analyse the subject HAAV’s transitions numerically, focusing on one critical concern of whether the impact and breach loads could be withstood by the fuselage’s proposed structure of syntactic foam – a material primarily chosen to survive hydrostatic pressure in compression. Such analysis also helps bound the deceleration function of the dual-EDFs considered in balance with sufficient time to exchange buoyancy underwater and minimise the effect time of incident waves for higher sea state tolerance. This paper numerically analyses a three-dimensional model of their HAAV design transitioning to and from the air-water interface. Using a CFD software with a Volume of Fluid method and Overset Meshing of the region of interest, similar to [30], we can freely move the vehicle model along a singular axis and develop the splash motion of the multiphase domain. The simulation was initially done with the model descending, to examine the slam loads and moments on various parts of the vehicle at varying speeds as it entered the water. Then, in the ascending (breach) direction, to observe the effects of water shedding and surface deformation on vehicle dynamics. In both instances, a still column of seawater with a flat free surface was replicated to establish a confirmable baseline from which to later simulate the dynamic wave action of the vehicle's intended operational environment. The intent of this initial study is to provide an instrument for analysing the structural stresses encountered by a vehicle and the validity of construction/design materials, whilst generating impact and tensile load cases to optimise the HAAV’s reinforced syntactic foam construction.

Our research significance is in the novelty of the HAAV design to operate in diverse sea states from substantial ‘breach-enabling’ depths and by morphing both the wings and retracting the dual-EDFs during transition. We prepare a baseline analysis to confirm the requisite loads and timings for the dive and breach, which can then be validated with pool experiments before further numerically analysing the vehicle’s vulnerabilities in realistic sea state waves. The remainder of the research article covers the method, results, and discussion of this initial numerical loads analysis during the transitions.

2. Methodology

2.1. Numerical Setup

The numerical work was generally informed by the texts by Ferziger, et al. [39], Tu, et al. [40], Wilcox [41], and Versteeg and Malalasekera [42].

2.1.1. Governing Equations

Given that the interaction velocities in all cases are far below the acoustic speed of water, only the incompressible components of the Reynolds Averaged Navier-Stokes (RANS) equations are considered. In their Cartesian form, they respectively equate to:

$${\displaystyle \frac{\partial \mathit{\rho }}{\partial \mathit{t}}}+\nabla \bullet \left(\mathit{\rho }\mathbf{v}\right)=0$$

(1)

$${\displaystyle \frac{\partial \left(\mathit{\rho }{\mathit{u}}_{\mathit{i}}\right)}{\partial \mathit{t}}}+\nabla \bullet \left(\mathit{\rho }{\mathit{u}}_{\mathit{i}}\mathbf{v}\right)=\nabla \bullet {\mathit{\tau }}_{\mathit{i}}+\mathit{\rho }{\mathit{b}}_{\mathit{i}}$$

(2)

where ρ is the fluid density, t is the time scale, $\mathbf{v}$ is fluid velocity, u_i is the Cartesian components of the averaged velocity, τ_i is the Newtonian Stress Tensor Cartesian component, and b_i is the coordinate body forces per unit mass [19]. In their incompressible form, both can be simplified to:

$$\nabla \bullet \mathbf{v}=0$$

(3)

$${\displaystyle \frac{\partial {\mathit{u}}_{\mathit{i}}}{\partial \mathit{t}}}+\nabla \bullet \left({\mathit{u}}_{\mathit{i}}\mathbf{v}\right)=\nabla \bullet \left({\mathsf{\nu }\nabla \mathit{u}}_{\mathit{i}}\right)-{\displaystyle \frac{\nabla \bullet \left({\mathit{p}\mathit{i}}_{\mathit{i}}\right)}{\mathit{\rho }}}+{\mathit{b}}_{\mathit{i}}$$

(4)

where kinematic velocity, ν, equates to fluid viscosity over density (μ/ρ), p is the mean pressure, and i_i is the Cartesian unit vector in the direction of x.

To account for the turbulent eddies that would be experienced, Wilcox’s two-equation k-omega [41] (k-ω) with Menter’s Shear Stress Transport (SST) [43] models were used and are defined by:

$${\displaystyle \frac{\partial \left(\mathit{\rho }\mathit{k}\right)}{\partial \mathit{t}}}+\nabla \bullet \left(\mathit{\rho }\mathit{k}\overline{\mathbf{v}}\right)=\nabla \bullet \left[\left(\mathit{\mu }+{\mathit{\sigma }}_{\mathit{k}}{\mathit{\mu }}_{\mathit{t}}\right)\nabla \mathit{k}\right]+{\mathit{P}}_{\mathit{k}}-\mathit{\rho }{\mathit{\beta }}^{\mathrm{*}}\left(\mathit{\omega }\mathit{k}-{\mathit{\omega }}_0{\mathit{k}}_0\right)+{\mathit{S}}_{\mathit{k}}$$

(5)

$${\displaystyle \frac{\partial \left(\mathit{\rho }\mathit{\omega }\right)}{\partial \mathit{t}}}+\nabla \bullet \left(\mathit{\rho }\mathit{\omega }\overline{\mathbf{v}}\right)=\nabla \bullet \left[\left(\mathit{\mu }+{\mathit{\sigma }}_{\mathit{\omega }}{\mathit{\mu }}_{\mathit{t}}\right)\nabla \mathit{\omega }\right]+{\mathit{P}}_{\mathit{\omega }}-\mathit{\rho }\mathit{\beta }\left({\mathit{\omega }}^2-{\mathit{\omega }}_0^2\right)+{\mathit{S}}_{\mathit{\omega }}$$

(6)

where $\overline{\mathbf{v}}$ is the mean velocity, μ is the dynamic viscosity, ${\mathit{\beta }}^{\mathit{*}}$ are model coefficients, and ${\mathit{S}}_{\mathit{k}/\mathit{\omega }}$ are applicable source terms [44].

The model coefficients are defined as:

$${\mathit{\sigma }}_{\mathit{k}}={\mathit{F}}_1{\mathit{\sigma }}_{{\mathit{k}}_1}+(1-{\mathit{F}}_1){\mathit{\sigma }}_{{\mathit{k}}_2}$$

(7)

$${\mathit{\sigma }}_{\mathit{\omega }}={\mathit{F}}_1{\mathit{\sigma }}_{{\mathit{\omega }}_1}+(1-{\mathit{F}}_1){\mathit{\sigma }}_{{\mathit{\omega }}_2}$$

(8)

$$\mathit{\beta }={\mathit{F}}_1{\mathit{\alpha }}_1^{\mathrm{*}}+(1-{\mathit{F}}_1){\mathit{\beta }}_2$$

(9)

Turbulence production terms, ${\mathit{P}}_{\mathit{k}}$ and ${\mathit{P}}_{\mathit{\omega }}$, are defined as:

$${\mathit{P}}_{\mathit{k}}={\mathit{G}}_{\mathit{k}}+{\mathit{G}}_{\mathit{n}\mathit{l}}+{\mathit{G}}_{\mathit{b}}$$

(10)

$${\mathit{P}}_{\mathit{\omega }}={\mathit{G}}_{\mathit{\omega }}+{\mathit{D}}_{\mathit{\omega }}$$

(11)

where:

$$\mathrm{Turbulent}\mathrm{Production},{\mathit{G}}_{\mathit{k}}={\mathit{\mu }}_{\mathit{t}}{\mathit{f}}_{\mathit{c}}{\mathit{S}}^2-{\displaystyle \frac{2}{3}}\mathit{\rho }\mathit{k}\nabla \bullet \overline{\mathbf{v}}-{\displaystyle \frac{2}{3}}{\mathit{\mu }}_{\mathit{t}}(\nabla \bullet \overline{\mathbf{v}}{)}^2$$

(12)

$$\mathrm{Non}-\mathrm{linear}\mathrm{Turbulent}\mathrm{Production},{\mathit{G}}_{\mathit{n}\mathit{l}}=(\mathrm{T}\mathrm{RANS},\mathrm{NL})\text{:}\nabla \overline{\mathbf{v}}$$

(13)

$$\mathrm{Buoyancy}\mathrm{Production},{\mathit{G}}_{\mathit{b}}=\mathit{\beta }{\displaystyle \frac{{\mathit{\mu }}_{\mathit{t}}}{\mathit{P}{\mathit{r}}_{\mathit{t}}}}(\nabla \overline{\mathit{T}}\bullet \mathit{g})$$

(14)

$$\mathrm{Specific}\mathrm{dissipation}\mathrm{Production},{\mathit{G}}_{\mathit{\omega }}=\mathit{\rho }\mathit{\gamma }[{\mathit{S}}^2-{\displaystyle \frac{2}{3}}(\nabla \bullet \overline{\mathbf{v}}{)}^2)-{\displaystyle \frac{2}{3}}\mathit{\omega }\nabla \bullet \overline{\mathbf{v}}]$$

(15)

$$\mathrm{Cross}-\mathrm{diffusion}\mathrm{term},{\mathit{D}}_{\mathit{\omega }}=2\mathit{\rho }(1-{\mathit{F}}_1){\mathit{\sigma }}_{{\mathit{\omega }}_2}{\displaystyle \frac{1}{\mathit{\omega }}}\nabla \mathit{k}\bullet \nabla \mathit{\omega }$$

(16)

And, Turbulent eddy viscosity is calculated as:

$${\mathit{\mu }}_{\mathit{t}}=\mathit{\rho }\mathit{k}\mathit{T}$$

(17)

$$\mathrm{where}T=\mathrm{m}\mathrm{i}\mathrm{n}({\displaystyle \frac{{\alpha }^{\mathrm{*}}}{\omega }},{\displaystyle \frac{a_1}{S{F}_2}})$$

(18)

S is the modulus of the mean strain rate tensor and is related to the velocity gradient of the fluid, while F, f, and $\gamma$ are additional model coefficients.

The ability of the combined k-ω turbulence modelling with SST equations to predict flow separations in adverse pressure gradients makes their use in predicting cavitation development during water entry and exit attractive, against the standard k-ε model [40].

2.2. Simulation Setup

In all cases, the simulations were examined using the CFD software StarCCM+ with a second-order implicit unsteady solver – the same package and solver combination used by Dong, et al. [34] and Sun, et al. [45] to examine their own HAAV’s free fall into water. In contrast, other researchers have in the last two years effectively used the same CFD approach in ANSYS Fluent to perform two-phase HAAV transition analysis [32, 46]. StarCCM’s unsteady solver employs the SIMPLE scheme to resolve the pressure-velocity coupling within the multiphase flow, a technique further described by Versteeg and Malalasekera [42]. The motion within the simulation was defined by the Dynamic Fluid-Body Interaction (DFBI) module within the software, which incorporates a six-degree-of-freedom solver that integrates pressure and shear forces across the surfaces of embedded bodies. The resultant forces and moments across these faces are then used to solve the governing equations of motion to find the body’s new position. At the same time, the adjusted boundary conditions of the body surface feed into the computation of the surrounding flow field [44]. With a minimum Weber number far greater than 1 (We ≈ 547,000) when using D as the characteristic length, the inertial forces considerably dominate the surface tension forces of the water [47], to the degree that they were considered negligible.

2.2.1. Model Geometry and Parameters

UNSW’s HAAV Mk.1 is designed to operate as a child vehicle to a parent aircraft, such as the P-8, as shown in Error! Reference source not found.. During operation, the vehicle would detach from the mothership and glide to its area of operation. To limit the slam forces on dive entry, it would conduct numerous manoeuvres, including a reverse thrust descent, to reduce its velocity to a maximum water penetration speed of 20 m/s. Before water entry, both wings and twin EDF propulsors would retract inboard to present the slimmest profile possible. Once submerged, both propulsors would unfurl and work in concert with the buoyancy control system to descend the vehicle to its operational depth. On completion of operations, it would ascend with its EDFs supplementing its positive buoyancy, to breach at a velocity where it can deploy its wings and accelerate away to marry with the parent aircraft once again. It would then be winched back onto its mounting pylon, similar to the General Atomics Sparrowhawk vehicle [27], where it would upload its data and be recharged, ready for re-use.

Rather than scaling based on Reynolds numbers, the full-sized HAAV model was utilised and reconfigured using Fusion360 for the air-water transition by folding the wings rearward without scissoring over one another. The hinge point was subsequently located on the wing root’s trailing edge, though sufficiently outboard to not result in scalloping of the fuselage, where the hinged wing roots would intersect. Table 1 summarises the key vehicle dimensions and parameters in both its cruise and transition configurations.

The initial vehicle weight was estimated to be 218.6kg while airborne. However, given its volume and average seawater density, its weight during entry was increased linearly with decreasing velocity until it was neutrally buoyant, remaining consistent with the conceptual design intent of flooding compartments with seawater. Conversely, to account for its required buoyancy control and free flooding compartments, the vehicle’s weight on exit was rapidly decreased beyond its buoyancy means in the first meter of flight, to simulate expunging open compartments of water. The centre of gravity (CoG) was moved from its optimal point at the wing’s quarter chord to 10 cm behind the trailing edge to account for the sweep of the battery-laden wings. Noting that most underwater gliders have a mechanism to alter their CoG longitudinally between at least two positions [48], but without yet knowing the particular configuration of the HAAV’s buoyancy system, the centre of buoyancy was kept static and coexistent with the CoG throughout the mass range.

Because the model was moving through the larger domain, the local coordinate system was aligned with the universal system, with the vehicle's upper surface in the -X direction and long axis in the +Z direction as shown in Error! Reference source not found.. Due to the moment developed by the underslung wings and wing roots interacting with the water, some vehicle rotational motion was expected. However, the trailing flat surfaces and controls of the vehicle would be expected to stabilise its motion during transition. So, to simplify the model and gather peak data, vehicle motion was restricted to a single degree of motion along the Z axis, as opposed to allowing the model to freely move and rotate.

2.2.2. Grid Generation and Boundary Conditions

Two alternative full-scale domains were created for the differing scenarios as depicted in Error! Reference source not found.. Both domains were broken into four regions, consisting of background, overlap, and volume of interest regions, with a horizontal refined volume to model the water surface.

Rather than re-meshing the entire domain every timestep, the volume of interest region was cast as an overset region to model the vehicle’s movement. An alternative to the overset approach has been proposed and proven by Bhalla, et al. [49] but the benefits are not yet clear. The overset works in two parts: first, by meshing the background without any overset bodies, and second, by meshing only the overset region with the body present. During the simulation, the properties of the background cells covered by the body and the overlapping grid are deactivated, while a thin shell of cells on the overset remains active in both regions. These cells trade information on the properties inside the overset region and the outer domain [39].

The Overset volume measured 0.97 x 1.35 x 1.81L_S to achieve even and adequate mesh growth across the vehicle. An Overlap region of the same width and depth was also created; however, it was extended the entire length of the Background volume to facilitate the meshing interface between the Overset and Background volumes as the vehicle moved through the domain. Overset volume walls were set as interfaces while the vehicle surface was a non-slip wall with 6 body-fitted prism layers expanding 0.008L_S from the model surface.

For the entry cases, a Background volume measuring 2.7 x 2.7 x 9.59L_S was used, similar to [45], with the high-density (water) surface located at Z = 1.64L_S as measured from the top of the domain. Such a gap ensured an adequate low-density volume (Air), to precipitate dropping the HAAV without boundary effects influencing the fluid interface during vehicle entry. The vehicle was positioned nose down approximately 0.08L_S above the simulated fluid surface to limit its free air acceleration prior to water contact. The domain lower boundary was located 8.95L_S below the water surface to allow the vehicle to develop appropriate entry splash effects and trailing cavitation while observing trends encountered during its varying mass descent.

In the Breach case, the background volume was reduced in total height to 7.70L_S, with the high-density (water) surface located at Z =2.98L_S from the domain top to account for the vehicle leaping from the surface during breach. The vehicle was positioned 2.0L_S below the surface prior to release, to allow water flow-surface interactions to adequately develop.

Figure 4. Domain and meshing environment of the Penetration cases (left) and Breach cases (right).

Figure 4. Domain and meshing environment of the Penetration cases (left) and Breach cases (right).

3. Grid Convergence

Two refinement levels to the base grid size were used in both entry and breach cases to establish grid convergence and determine the influence of cell length, dl, on the numerical results. The timestep was determined for each grid convergence scenario by ensuring a Courant-Friedrichs-Lewy (CFL) value of 1 for the varying refinement levels. The penetration cell values across the fine, medium, and coarse refinements were on the order of 7.7x10⁶, 1.7x10⁶, and 5.2x10⁵, with corresponding minimum dl’s of 4.0, 7.5, and 15 mm, respectively. Error! Reference source not found. shows the grid independence results regarding the Velocity behaviour of a 15 m/s penetration. In the figure, we see that there is a good agreement between all three refinements and likely considerable scope to optimise solution efficiency if necessary. Notwithstanding, the medium case appeared optimal of all three for a two-second simulation.

Figure 5. Mesh refinement results comparison.

Figure 5. Mesh refinement results comparison.

4. Results and Discussion

All of the simulations and results were conducted on a computer using an eight-core Ryzen 7 3700X CPU at a standard 3.6 GHz processing speed, requiring up to 16 hours of continuous computing to determine a solution.

4.1. Penetration

Three initial velocities, V_i, were examined for the penetration cases using the medium base mesh and under a common time step of Δt = 0.0075 s. After a short settling period to allow the initialised values to stabilise, the vehicle was allowed to freely drop under gravity with a V_i of 10, 15, and 20 m/s. The reason for examining this spread of entry velocities was to account for the vehicle’s ability to enter at a specific speed and consider the vertical fluid motion of passing waves. In all cases, the vehicle’s Ogival nose cone would break the water surface, initially forcing fluid to jet outward to the sides. As it continued to descend, more fluid was pushed to the sides of the nose, plateauing the impact force at its peak. Rather than inducing a pressure wave in the fluid from the point of entry, incompressibility-induced splash dynamics fountained water up the sides of the vehicle in the beginnings of a splash curtain as a cavity is formed where the model occupied the space. The formation of these splash effects has been previously studied by Korobkin and Pukhnachov [50] and well outlined in Truscott, et al. [51]. Simultaneously, viscosity influenced the incoming water mass relative to the vehicle to decrease the overall dynamic pressure on the nose cone. The decrease in dynamic pressure allowed the vehicle to continue accelerating from gravity, initially at a significantly reduced rate, before reaccelerating with the decreasing dynamic pressure. The additional vehicle speed and flow of the splash curtain up the sides of the fuselage interacted with the fuselage wing roots, introducing the first additional drag source. However, the vehicle would continue to accelerate until the velocity peaked just before the root caps of the swept wings contacted the surface. Deceleration was sudden on entry of the wings and typically represented the point of maximum drag. The drag forces experienced by the wing surfaces are delayed due to the blunt body separation of the flat root caps creating their own cavities, until hydrostatic pressure and gravity reattach the flow to the wings. On reattachment, the trailing turbulent eddy interaction between the wings and horizontal stabiliser inhibited their drag force development until the entire body was finally submerged. The dynamics and associated forces for the water entry are visualised in Error! Reference source not found. and Error! Reference source not found. respectively.

Figure 6. Vehicle surface penetration dynamics showing surface streak lines and areas of dynamic pressure intensity.

Figure 6. Vehicle surface penetration dynamics showing surface streak lines and areas of dynamic pressure intensity.

Figure 7. Vehicle Dynamics Variation during water-entry.

Figure 7. Vehicle Dynamics Variation during water-entry.

Comparing the acceleration and velocity profiles of the different initial velocity cases, three key deceleration peaks throughout the transition through the water surface were observed. The peaks were subsequently associated with the nose cone, root cap surfaces, and water entry of the full horizontal stabiliser span. Vehicle deceleration ranged between 5G’s and 11G’s, peaking at the root cap transition, while the follow-on decelerative impulse appears more acute with increasing velocity values; a feature also apparent in the studies of Dong, et al. [34]. Both the peak values and their following impulses have important implications for the internal components, as it was determined that these parameters were central to defining both hardening and survivability requirements for internal components in the Z-axis. We also see two surprising features between the acceleration and velocity values. The first being the high deceleration of the 20 m/s entry case, slowing the vehicle to a similar velocity as the 10 m/s entry around t = 0.4 s. It should be noted, the 20 m/s case does not achieve depth faster than the 10 m/s case, but rather that the impact energies are transferred from the vehicle into the fluid over similar periods. The second observation is the delayed reduction in deceleration by the 10 m/s entry between t = 0.4 s and t = 0.6 s, possibly due to trailing cavity interactions.

Figure 8. Vehicle Force Variation during water-entry.

Figure 8. Vehicle Force Variation during water-entry.

Error! Reference source not found. shows that the vehicle’s nose cone force ranged from 1,500 to 4,500 Newtons (N). Not to be confused with the peak dynamic pressures (which approximately range between 7.75x10⁴ and 2.6x10⁵ Pascals), these values represent the local force experienced at the interface of the nose cone’s base, to be translated up the length of the fuselage barrel. Unsurprisingly, they mirror the acceleration values highlighted above and are similar in behaviour to the Zhang, et al. [52] study of a torpedo’s entry into water.

Once the swept wing bases contact the water surface, the impact forces peak between approximately 8,700 N and 21,000 N and are transmitted primarily through the wing hinges. Interestingly, the peak horizontal stabiliser force for the 10 m/s case has a significantly longer impulse than the other cases, and is more graduated. This is likely reflective of the turbulent attachment of flow to the tail cone as compared to the trailing cavity production of the higher speed cases. In all cases, we see the effect of turbulent flow interaction between the upper wing surfaces and the lower horizontal stabiliser surfaces, as the collapsed wing root cavity remnants are loaded and accelerated, consequently producing an oscillating source of shear drag on the affected surfaces.

4.2. Breaching

Prior to conducting any of the Breach cases, the drag coefficient of the vehicle (C_D) needed to be determined to ascertain the buoyancy component of the Breaching Velocity, V_B. This was done by setting the vehicle into a standalone volume, the same dimensions and mesh parameters as the Breach case, and running a steady state simulation with seawater flowing at 10 m/s as the medium. With the C_D determined, it was assumed that the drag induced by the terminal velocity of a buoyancy-driven ascent would balance the buoyancy force. Therefore, the remaining velocity to achieve V_B, combined with the C_D would give the residual propulsive force necessary to maintain velocity. For simplicity, the propulsor’s thrust was maintained throughout the simulation, including when the propulsors had notionally cleared the water’s surface, and reverted to their air-breathing performance parameters. In reality, both EDFs would shut down just prior to transitioning the surface and then restart once purged of water, to avoid overspeed.

Three breach cases were examined based on the buoyancy of the vehicle, those being the 10%, 20%, and 30% positive buoyancy. A value of 30% was chosen as the highest buoyancy due to the likelihood that the EDF stowage compartment would be free flooding and could not contribute to buoyancy control, leaving the remaining volume to contribute. As highlighted earlier, the weight in air of the vehicle is estimated at 218.6 kg, while the neutral buoyancy weight is 527 kg. To achieve a 10% buoyancy force, the vehicle must weigh 10% less than the neutral weight, or 52.7 kg lighter. A subsequent vehicle mass of 474.3 kg for the same volume would result in approximately 517 N of positive buoyancy, ${\mathit{F}}_{\mathit{B}}$. As outlined above, the total drag during an ascent at V_B is balanced by both the buoyancy and propulsive forces. Incorporating them into the drag equation, we get:

$$\mathit{D}={\sum \mathit{F}=\mathit{F}}_{\mathit{B}}+{\mathit{F}}_{\mathit{P}}={\displaystyle \frac{1}{2}}\mathit{\rho }{\mathit{V}}_{\mathit{B}}^2{\mathit{S}}_{\mathit{z}}{\mathit{C}}_{\mathit{D}}$$

(19)

where ρ is the average density of sea water, ${\mathit{S}}_{\mathit{z}}$ is the cross-sectional area of the vehicle facing its direction of travel, and C_D is the coefficient of drag equal to 0.3784. Inputting table 1 values into Equation 19 with a $F_B$ of 517 N, 1034 N, and 1551 N for the increasing buoyant cases, the thrust output of the propulsors, $F_P$, necessary to achieve a V_B of 5 m/s is 913 N, 396 N, and -121 N, respectively.

With the thrust values known, we were able to incorporate them into the simulations to achieve a steady velocity while submerged.

For the 30% buoyancy case, the initial acceleration and final deceleration of the vehicle are higher than in the other two cases due to the excessive amount of buoyancy necessary to achieve V_B. As the effect of the fluid force during breach could not be accounted for, zero propulsive force was applied, allowing the vehicle to transition “naturally” rather than attempting to regulate its V_B. As can be seen from Error! Reference source not found., in all three cases the vehicle is sufficiently able to clear the water surface, with the 20% and 30% cases only marginally able to do so. With the propulsors located so close to the nose of the vehicle, it is reasonable to conclude that in all cases, they are able to transition to an airbreathing operation in waves as high as 2.5 m. Further scrutiny in Error! Reference source not found. and Error! Reference source not found. also reveals that as the vehicle approaches the surface, there is an increased acceleration. The sudden increase in velocity is potentially a result of fluid viscosity at the free surface, allowing for momentum to be transferred from the vehicle to the water column, in the Z direction, rather than forcing fluid to the side. The result of the fluid’s upward momentum is the formation of a dome over the nose cone as it departs the water (Error! Reference source not found.), which is consonant with the simulations of Chen, et al. [53] and Moshari, et al. [54]. A similar effect is seen at the point where the root caps also breach the surface. There is also an initial variation in the deceleration values as the water from the propulsor stowage compartment is evacuated of water, within its first meter of exit, between 0.6 and 1.6m of height, easing its deceleration gradient (Error! Reference source not found.). Lastly, we see the effect of the vehicle’s remaining buoyancy and propulsive force, delaying the onset and eventually offsetting gravitational acceleration.

Figure 9. Vehicle Breach Behaviours.

Figure 9. Vehicle Breach Behaviours.

Unfortunately, with the lack of a retarding force to limit the vehicle’s velocity, the 30% buoyancy case can accelerate to a post-release V_B 3.3% greater than the remaining cases. Although a minor difference, this slight change in velocity has resulted in an upward left shift of the 30% case’s total force profile. The difference in the shift demonstrates the sensitivity of the vehicle’s breach loads with varying velocities.

For the remaining buoyancy cases, we see that they mirror each other closely up to the point where the nose cone fully departs the water at t =1.6 s, and drains the remaining surface dome fluid. The difference in buoyancy, as a fraction of the remaining vehicle weight, makes itself evident in the left shift of the 20% and 30% profiles post t = 1.75 s.

Following the departure of the root caps centred around t = 1.90s, we find that the key force inhibiting the upward motion of the vehicle disappears, resulting in a dramatic reversal in the total force, from an overall drag effect, to an upward propulsion, due to the buoyancy of the remaining vehicle portion, which is still submerged. There is a moment of pseudo-stability between t = 2.00s and t = 2.15s as the weight of the draining fluid from the emerging root caps and vehicle stifles the propulsive buoyancy (Error! Reference source not found.). The weight of water clinging to the vehicle to be shed is also noted in simulations by Chen, et al. [46], Sun, et al. [45] and accounted for by Yun, et al. [36], when predicting breach heights. Finally, there is an equivalent reduction in buoyant propulsive force as the horizontal stabilisers begin to emerge, culminating at the point where the wing tips exit the surface.

Figure 10. Vehicle surface breach dynamics showing surface distension and fluid carriage.

Figure 10. Vehicle surface breach dynamics showing surface distension and fluid carriage.

Figure 11. Vehicle forces during breach.

Figure 11. Vehicle forces during breach.

5. Conclusions

This research article has provided the first journal overview of a new hybrid aquatic-aerial vehicle design concept that aims to operate as a child aircraft to a maritime patrol parent aircraft, diving underwater with swept wings and retracted electric ducted fans, to conduct sonar search, before returning via a significant breach and return flight, with wings and fans out. The HAAV receives recharge and comprehensive data analysis from the parent aircraft. The vehicle’s significant design challenge is to pervasively operate in sea states that no HAAVs have yet achieved. To exploit new lightweight materials with compressive strength and to time the trans-media morphing for variable sea state, this article presents the first findings of the HAAV’s free dive entry and controlled breach, from a static column of seawater. We used numerical simulation software to determine structural load cases and to inform the next design iteration of the vehicle's morphing and propulsion. Combining RANS solvers with SST k-omega turbulence models, we were able to conduct simulations for three different entry velocity cases, three different breach buoyancy cases, and grid independence studies. The vehicle dynamics during the transition through the water surface and the forces acting on key parts of the vehicle were able to be determined and compared across all cases, for both scenarios. Additionally, we were able to demonstrate the validity of our propulsive force calculations to maintain breach velocity during a powered ascent, prior to determining the vehicle’s ability to sufficiently clear passing wave tops to climb away from the surface.

We found that the vehicle’s behaviours on entry were very similar, such that only the amplitude of the behaviours would vary with each case. Across all entry velocities, once the nose cone had fully submerged, viscosity and splash dynamics would accelerate the vehicle until interaction with the wing roots and the swept wing bases, resulting in significant decelerations on the order of 11G’s for the 20 m/s entry case. For the same case, maximum loads experienced at the nose-cone to fuselage interface and wing hinges were of the order of 4,500 and 21,000 N, respectively. Accelerated flow and turbulent eddy interaction between the folded wings and the overhead horizontal tail would result in the strong oscillatory moments toward each other, varying the total drag, which tapered with decreasing velocity. For the breach cases, the energy expended in creating an exit splash only marginally impacted the vehicle’s ability to leap from the water's surface. Importantly, the degree of buoyancy and buoyancy control appear to have the greatest impact on a vehicle’s exit, while form drag and water drainage largely affect the exit velocity behaviour. The velocity behaviour in turn determined the vehicle’s maximum leap height post breach, which ranged between 1.08 and 1.19L_S (4.0 and 4.4m) as measured from nose to water surface.

Further work using plunge pool experiments will be necessary to validate results and refine CFD parameters to reflect the HAAV’s nuanced specifications. Already, wind tunnel testing with the HAAV in its glide configuration has been conducted, with future work planned to examine its dive configuration characteristics, while other work will examine and refine the conformity of the wing fold system. The validated CFD from the pool experiments will be used with representative numerical sea waves to confirm the viability of the HAAV to pervasively operate in high sea states.