Inspired by recent studies of a squid-like swimmer, we propose a three-dimensional jet propulsion system composed of an empty chamber enclosed within a deformable body with an opening. By prescribing the body deformation and jet velocity profile, we numerically investigate the jet flow field and propulsion performance under the influence of background flow during a single deflation procedure. Three jet velocity profiles, i.e., constant, cosine and half cosine, are considered. We find that the maximum circulation of the vortex ring is reduced at a higher background flow velocity. This is because stronger interaction between the jet flow and background flow makes it harder to feed the leading vortex ring. Regarding thrust production, our analysis based on conservation of momentum indicates that with the constant profile the peak thrust is dominated by the time derivative of the fluid momentum inside the body, while momentum flux related thrust accounts for the quasi-steady thrust. For the cosine profile, its peak is mainly sourced from momentum flux associated with the unsteady vortex ring formation. No prominent thrust peak exists with the half cosine profile whose thrust continuously increases during the jetting. For all the three jet velocity profiles, added-mass related thrust attributed to body deformation enhances the overall thrust generation non-negligibly. Under the present tethered mode, the background flow has negligible influence on the thrust attributed to momentum flux and momentum change of the fluid inside the body. However, it indeed affects the over pressure-related thrust but its effect is relatively small. The overall thrust declines due to the significantly increased drag force at large incoming flow speed despite the rise of added-mass related thrust. Unsteady thrust involving vortex ring formation becomes more important in the overall thrust generation with an increased background flow velocity, reflected by larger ratios of the unsteady impulse to jet thrust impulse.

## I. INTRODUCTION

Jet propulsion utilized by squids and other cephalopods is a unique aquatic locomotion method different from the commonly used fin oscillation and body undulation modes. Its simplicity and effectiveness (as demonstrated by the extremely high-speed bursting capacity of squids) have attracted many researchers in quest for bio-inspired propulsion systems of underwater robots or vehicles. A typical jetting cycle of a squid starts with the inflation of the mantle, a muscular organ that encloses the mantle cavity and the viscera (Bartol , 2008). During this inflation, the internal volume of the cavity increases so that water fills in. Afterward, the squid shrinks the mantle and water is ejected through the funnel tube. Through body deflation, effective escape locomotion and high maneuverability can be achieved (O'Dor, 2013). Some squids are even capable of bursting as fast as the fastest fish in terms of peak speed measured by body length per second (Gosline and DeMont, 1985; Wardle, 1975).

The studies of squid locomotion started with direct observation and measurement of live squids. For example, anatomical studies indicated that the squid mantle contains two predominant muscle groups, i.e., circumferential muscle fibers (or so-called circular muscles) that constitute the bulk of the mantle wall and radial muscles that extend from the inner to the outer surface of the mantle wall (Ward and Wainwright, 2009). The activity patterns of these two muscles during jet-propelled locomotion were recorded by Gosline (1983). The mantle muscle structure and related kinematics were further studied by Thompson and Kier (2006). These studies set the basis of our understanding of the morphology and physiology of squid jetting.

Beyond that, efforts have also been made to investigate the hydrodynamics and flow patterns around live squids for the purpose of designing an efficient jet swimmer. Based on recorded kinematic data of swimming squids by high-speed camera, Anderson and DeMont (2000) calculated jet velocity, thrust generation, and propulsive efficiency. However, their approximation based on the Bernoulli equation led to significant inaccuracy since important physical mechanisms related to vortex dynamics were omitted. In subsequent studies, more detailed flow structures of the jets behind squids were reported via digital particle imaging velocimetry (DPIV) techniques (Anderson and Grosenbaugh, 2005; Bartol , 2016). Through the visualization of the flow field around the squid in pulsed jetting, two principal jet structures were identified, mode I and mode II (Bartol , 2009). In mode I, the ejected fluid rolls up into an isolated vortex ring; whereas in mode II, the discharged fluid develops into a leading vortex ring followed by a trailing flow with distributed vorticity.

In addition to the above studies which focus on the hydrodynamics of jet flow around live squids, some more fundamental research of vortex ring formation may also be valuable in elucidating the underlying mechanism of jet propulsion. For instance, similar flow structures of the vortex ring to that observed behind live squid were generated through a piston/cylinder arrangement in the laboratory environment (Gharib , 1998). They found that the ratio of the piston stroke to diameter *L*/*D* determines the transition between the two distinct states. Namely, with a constant jet velocity profile, only an isolated vortex ring is formed at a small stroke ratio (*L*/*D* < 4), while the flow field at a large stroke ratio with *L*/*D* beyond 4 consists of a leading vortex ring followed by trailing flow, reminiscent of the mode I and mode II, respectively. This phenomenon was confirmed by subsequent numerical studies (Mohseni , 2001; Rosenfeld , 1998). Furthermore, Linden and Turner (2001) performed a theoretical analysis to elucidate the underlying physics. They compared the circulation, impulse, volume and kinetic energy of the jet plug with the counterparts of a finite-core vortex ring. Their results showed that a single vortex ring cannot retain these properties after the critical equivalent stroke ratio is reached. Besides, they suggested the limiting vortex was “optimal” in achieving maximum impulse for given energy input. Indeed, recent studies of squid-inspired jet propulsion showed optimal propulsive performance was reached near the critical stroke ratio (Bi and Zhu, 2018; Christianson , 2020).

Based on the above findings, the vortex ring formation was then taken into consideration in a few studies of cephalopod-inspired jet propulsion. For example, a novel cephalopod-inspired jet propulsion system was proposed by Bi and Zhu (2018). They focused on a single bursting cycle and found the optimal speed coincides with the critical stroke ratio. However, their numerical model was based on the potential-flow theory, where the viscosity effect was not included so that the accuracy of the results was compromised. A subsequent study of an axisymmetric squid jetting model took the fluid viscosity into account by solving the Navier-Stokes equations (Bi and Zhu, 2020). Their results highlighted the role of jet acceleration and viscous dissipation in the vortex ring evolution. This study was characterized by an axisymmetric model. Moreover, like many other studies (e.g., Advaith , 2017; Krueger , 2008; Nguyen , 2019), this axisymmetric study was conducted in a still fluid environment, a scenario different from real squid locomotion in steady swimming. According to the measurement by Anderson and Grosenbaugh (2005), squid swims at a forward velocity not much lower than the jet velocity. The ratio between swimming speed and the jet velocity of adult squid ranges from 0.33 to 0.69 based on their observation. The background flow would affect the vortex ring formation and, thus, the jet propulsion performance (Anderson and Grosenbaugh, 2005). This effect was considered in a numerical study by Jiang and Grosenbaugh (2006). Nevertheless, their model was based on a rigid body, distinct from cephalopod locomotion involving body deformation. Actually, recent studies indicated that dynamic shape change plays a non-negligible role in the jet process (Bi and Zhu, 2019; Steele , 2017). Another limitation of the work by Jiang and Grosenbaugh (2006) is that their numerical model is two dimensional.

In this work, we numerically examined a three-dimensional cephalopod-inspired jet swimmer focusing on propulsion performance. This system has a pressure chamber with a fixed nozzle serving as the outlet for discharged fluid. The deformation of the body is prescribed so that we will be able to investigate the effect of equivalent stroke ratio on jet flow and pulsed thrust production under given jet velocity profiles. With background flow being considered, we aim at examining the vortex ring formation and jet propulsion capabilities at different jet velocity profiles in a low jet-based Reynolds number scenario.

Compared with previous numerical studies of the squid-like jet propulsion system (Bi and Zhu, 2018, 2020; Steele , 2017), the novelty of the present work lies in three aspects. First, we take incoming flow into consideration, which has been neglected in most previous studies. The presence of incoming flow plays a non-negligible role in the hydrodynamics of squids (Anderson and Grosenbaugh, 2005). Second, our squid-inspired system with an empty pressure chamber involves body deformation, which is much closer to actual squids and some existing robotic design. Finally, a three-dimensional jet propulsor modeling that takes viscous effect into consideration is introduced in our study.

The rest of this paper is organized as follows. The problem formulation, including model geometry, deformation kinematics, and performance metric, is introduced in Sec. II. These are followed by the governing equations and numerical techniques in Sec. III. Afterward, simulation results are presented in Sec. IV. Finally, conclusions are summarized in Sec. V.

## II. PROBLEM FORMULATION

A squid-like swimmer consisting of a pressure chamber with an opening is considered in this work. This chamber mimics the mantle cavity of the squid, and its opening serves as both the inlet and outlet of fluid. During a typical deflation-inflation process, the chamber body shrinks so that the fluid inside the cavity is expelled through the nozzle, and then fluid is sucked back into the chamber for the next jetting. In this particular study, we focus on a single deflation deformation that produces most of the thrust during the deflation-inflation propulsion cycle.

*x*-axis, as shown in Fig. 1(a). The body length of the swimmer is

*L*, and the nozzle size is denoted by

*D*. The upper surface [above the upper dashed line in Fig. 1(a)] and lower surface [below the lower dashed line in Fig. 1(a)] are two halves of an ellipse, respectively. In a cylindrical coordinate system, the radial distance from the

*x*-axis to a certain point on the body surface is given by

As shown in Fig. 1(b), during deflation, the length of the body *L* and nozzle diameter *D* = 0.2 *L* remain unchanged, while the eccentricity of the ellipse increases so that the volume of the internal body is reduced. This design is similar to an octopus-inspired robot in Weymouth (2015) whose internal skeleton is rigid, and therefore the body length of the robot is constant. This is also in line with the jet locomotion of a squid, during which the body length remains almost unchanged, but the perimeter decreases as the body shrinks (Ward, 2009).

*h*which equals 0.02

*L*at its depleted state is varied so that the overall volume of the solid body (excluding the fluid inside it) remains a constant to avoid the mass change of the body itself. The two extreme states, i.e., the fully inflated and deflated, are denoted by $ e = e 0$ and $ e = e 1$, whose internal volumes (the volume of the fluid which fills in the chamber) are $ V ( e 0 )$ and $ V ( e 1 )$, respectively. The equivalent stroke ratio associated with the body deflation is defined as $ \Gamma = 4 \Lambda / ( \pi D 3 )$, where $ \Lambda = V ( e 0 ) \u2212 V ( e )$ (

*e*is the instantaneous eccentricity). To characterize each case, we define the maximum equivalent stroke ratio $ \Gamma m$, which is closely related to body deformation and reached at the end of the body shrinking, as $ \Gamma m = 4 ( V ( e 0 ) \u2212 V ( e 1 ) ) / ( \pi D 3 )$. During the deflation, the spatially averaged speed of the ejected jet flow $ V j$ at the nozzle can be determined by

$ e 0$ . | 0.904 . | 0.898 . | 0.883 . | 0.868 . | 0.844 . |
---|---|---|---|---|---|

$ \Gamma m$ | 3.31 | 4.66 | 7.60 | 10.59 | 15.07 |

$ e 0$ . | 0.904 . | 0.898 . | 0.883 . | 0.868 . | 0.844 . |
---|---|---|---|---|---|

$ \Gamma m$ | 3.31 | 4.66 | 7.60 | 10.59 | 15.07 |

^{2}/s.

*x*-direction, $\rho $ represents the fluid density, and $P$ is the instantaneous power expenditure of the body, which can be calculated as

## III. NUMERICAL MODEL AND APPROACH

*u*,

*v*, and

*w*being the three velocity components in the Cartesian coordinate system and

*E*being the total energy of the flow. $n$ represents the normal unit vector in the outward direction and $\Omega $ denotes the fluid control domain with the boundary as

*S*. $F$ is the flux vector which consists of two parts, i.e., the convective and pressure fluxes $ F c$ and the fluxes $ F d$ arising from the viscous shear stress and thermal diffusion, with

*p*is the pressure. During the deformation of the swimmer body, the arbitrary Lagrangian–Eulerian (ALE) formulation is implemented to handle the flow equations on a deformable mesh. It is achieved by defining fluxes relative to the motion of the control domain surfaces, which are expressed in terms of relative velocity

*u*,

_{r}*v*, and

_{r}*w*given by

_{r}*i*,

*j*,

*k*), a semi-discretized form of Eq. (8) can be obtained as

*i*,

*j*,

*k*), and $ R i , j , k$ is the residual measuring the convective and diffusive fluxes entering the hexahedral cell through its surface. $ D i , j , k$ represents an artificial viscosity term which is introduced to stabilize the numerical scheme and eliminate the spurious numerical oscillations (Jameson , 1981).

*n*denote the two previous time levels and this yields a second-order accuracy. A hybrid multistage Runge–Kutta scheme is then implemented to integrate Eq. (16).

Within the numerical flow solver, parallelization is enabled by domain decomposition via message passing interface (MPI) to achieve large-scale computation. Furthermore, the local time-stepping and multigrid methods are implemented to accelerate the convergence. Meanwhile, implicit residual smoothing is applied to increase the stability of the solution. More details about the fluid solver can be found in Sadeghi (2004).

It is worth noting that when simulating an incompressible flow problem considered here using the present compressible flow solver, it is necessary to ensure that the compressibility is negligibly small. As we did in the previous applications to other relevant incompressible flow simulations, the freestream Mach number, defined as $ M a \u221e = U 0 / a \u221e$ with $ U 0$ denoting the freestream flow velocity and $ a \u221e$ being the local speed of sound, is chosen as 0.06 − 0.1. This value is far below the critical value of 0.3, where the compressibility effect becomes pronounced, but still sufficiently large to ensure numerical stability. Given that the boundary is moving during the deformation of the body and flow is ejected with a certain speed, the actual Mach number may become larger than $ M a \u221e$ around the swimmer. Therefore, the Mach number distribution in the whole computational domain is monitored during computation to ensure that they are always below the critical value for computational accuracy.

This compressible flow solver has been successfully validated and applied to various incompressible flow simulations involving boundary motion in our previous work (Liu , 2016; Luo , 2020a; Shi , 2020; Xiao and Liao, 2010). The validations include simulations of flow over a flexible cantilever behind a square cylinder, the bending of a flexible plate in the uniform flow, the response of a flexible plate in forced heave motion (Luo , 2020b), and flow over a NACA 0012 airfoil (Luo , 2020c).

The computational domain and generated fluid mesh around the body surface are shown in Fig. 3. On the body surface, the no-slip condition is applied, while for the other boundaries, the non-reflective far-field boundary condition is imposed. A self-consistency study is performed to assess the appropriate mesh and time step resolution by solving the laminar flow for $ U 0 = 0.42 \u2009 V j m$ and $ \Gamma m = 10.59$. Specifically, three grids are generated: a coarse grid with 5 109 112 cells, a medium-size grid with 6 437 992 cells, and a fine grid with 8127368 cells. The minimum grid spacing for the three meshes is $ 3.5 \xd7 10 \u2212 3$ *L*. Like the other simulations in this study, the body deformation is started after the surrounding quasi-steady flow field is fully developed at the presence of incoming background velocity $ U 0$. First, the mesh convergence test is performed in which three meshes are used along with the dimensionless time step size, defined as $ \Delta t \xaf = \Delta t V j m / L$, $ \Delta t \xaf = 0.02$, as shown in Fig. 4(a). It can be seen that the medium mesh yields a result that is very close to that of the fine mesh. Thereafter, the effect of the time step size is also studied, as shown in Fig. 4(b). These tests show that the results are not sensitive to the mesh size and time step size if they are sufficiently small. Therefore, the following simulations hereafter are carried out with a medium-size mesh at $ \Delta t \xaf = 0.025$ to reduce computational cost and ensure accuracy.

## IV. RESULTS

Similar to Jiang and Grosenbaugh (2006), we initiate the body deformation after the external quasi-steady flow field is fully developed. To elaborate, the fluid field where the inflated body is at rest in the presence of incoming uniform flow is first simulated and then taken as the initial flow field for the following computations when the body deflation is included.

### A. Effect of maximum equivalent stroke ratio (body deformation)

To visualize the wake structure around the swimmer body at different maximum equivalent stroke ratios, we plot the *Z*-vorticity contour, streamline, and *Q* criterion distribution at the plane *z* = 0, as well as the wake structure by iso-surface of *Q* criterion, in Fig. 5. It can be seen from Figs. 5(a) and 5(g) that at a small maximum equivalent stroke ratio almost all of the ejected fluid is entrained into a single vortex ring. When $ \Gamma m$ is increased to 10.59, on the other hand, the wake consists of a leading vortex ring followed by a trailing flow, as shown in Fig. 5(b). This is reminiscent of the two jet modes observed behind the jetting of squid, i.e., in mode I the ejected fluid rolls up into an isolated vortex ring, while in mode II a leading vortex ring pinches off from the trailing jet (Bartol , 2008). The cases for $ \Gamma m$ = 3.31 and 10.59 show characteristics of these two modes, respectively.

Nevertheless, a close inspection of Figs. 5(f) and 5(h) reveals that there are no remarkable secondary vortices formed behind the leading vortex ring at $ \Gamma m$= 10.59, unlike what was observed in Gharib (1998). This difference may be mostly attributed to the fact that larger jet-based Reynolds number $ Re j$ (above 2500) was considered in the experiment, while at a small Reynolds number (below 1000), our results are consistent with numerical and experimental data in Bi and Zhu (2020) and Palacios-Morales and Zenit (2013). The interaction between the jet flow and background flow near the nozzle can be observed from the streamline distributions in Figs. 5(c) and 5(d), implying that the incoming flow may affect the wake flow field and vortex ring formation, which will be examined later.

For further insight into the evolution of the near-body flow field and vortex ring, the *Z*-vorticity contour and *Q* criterion distribution at the plane *z* = 0 and the wake structure at several instantaneous equivalent stroke ratios during the jetting at $ \Gamma m = 10.59$ are presented in Fig. 6. However, it is important to note that the detachment of the leading vortex ring from its trailing jet cannot be determined by the *Q*-criterion contours presented here, as they exclude the region of the free shear layer which is the source of vorticity for the leading vortex ring growth. Although the critical value of stroke ratio was found to be around 4 by Gharib (1998), some subsequent studies suggested that the critical stroke ratio may vary depending on specific conditions, e.g., velocity distribution at the nozzle plane and jet-based Reynolds number (Mohseni , 2001; Palacios-Morales and Zenit, 2013; Rosenfeld , 1998). For example, Bi and Zhu (2020) found that at $ Re j = 150$, the trailing flow was not strong enough to induce a secondary vortex ring, and it remained connected with the leading vortex ring after $ \Gamma = 4$ so that the trailing flow was still able to feed the leading vortex ring.

*Q*criterion to obtain a closed area where

*Q*> 0 to integrate the vorticity and calculate the vortex ring circulation. The circulation

*C*is obtained by

*Z*-vorticity at plane

*z*= 0 given by

It is interesting to examine the effect of body deformation on the actual jet velocity profile. In most previous studies involving a piston arrangement (Gao , 2020; Gharib , 1998; Jiang and Grosenbaugh, 2006), the nozzle features a tubular shape so that the flow is mainly ejected in the axial direction. In comparison, our system is more like a conical nozzle as the lateral dimension decreases during the body shrinking. For this reason, there may be a significant radial component in the jet velocity, as shown from the streamline distribution in Figs. 5(c) and 5(d). For further insight into the actual jet velocity profile at the exit plane, the distribution of axial and radial velocity is presented in Fig. 8. As shown in Fig. 8(b), a maximum radial velocity of around 0.3 $ V j m$ is observed near the nozzle surface. On the other hand, we also observe a small peak of the axial velocity, which is larger than the jet speed near the position *r*/*D* = 0.35 − 0.45.

*A*represents the nozzle plane. A comparison between Figs. 9(a) and 10(a) reveals that the jet momentum flux related thrust is not the main contributor to the total thrust during the initial stage when $ C T$ is much larger than $ C T j$, although at the quasi-steady phase $ C T j$ approaches $ C T$.

We further plot the excessive pressure-related thrust coefficient $ C T p = F p / ( 0.5 \rho V j m 2 D 2 )$ in Fig. 10(b), where $ F p = \u222b A \Delta p d S$ ( $ \Delta p = p \u2212 p \u221e$, with $ p \u221e$ being the pressure in the far-field) is the excessive pressure force at the nozzle plane. Nevertheless, by itself $ C T p$ is not able to bridge the difference between $ C T j$ and the total thrust $ C T$, especially at its peak. For instance, at $ \Gamma m = 15.09$, the peak value of $ C T$ is 3.37 while the maxima of $ C T p + C T j$ is only 1.81. The remaining thrust comes from $ F m$, whose normalized form $ C T m = F m / ( 0.5 \rho V j m 2 D 2 )$ is shown in Fig. 10(c). Taking $ \Gamma m = 15.07$ as an example, we find that the overall peak thrust at $ t = 0.13 \u2009 T d$ is mainly dominated by the horizontal momentum change of the fluid inside the chamber $ C T m$, as shown in Figs. 9(a) and 10(c). By comparing Fig. 9(a) with Fig. 10(a), the thrust component attributed to jet momentum flux $ C T j$ dominates thrust generation roughly after $ t = 0.4 \u2009 T d$. Although previous studies suggested that over pressure effect would enhance the generated thrust of the starting jet (Krueger and Gharib, 2003), we find its contribution $ C T p$ is relatively small during the whole procedure. In a comparison of Figs. 10(b) and 10(d), added-mass related thrust (which will be discussed later) attributed to body deformation even contributes more thrust than $ C T p$. The correlation of overall thrust production between the three thrust sources applies to other $ \Gamma m$ as well.

Besides, we further quantitatively compare the three thrust components with the total thrust at $ \Gamma m = 15.07$ in Fig. 12(a). It can be seen that the peak value of the sum of the three parts matches well with that of the total thrust. However, during the quasi-steady phase, the combination of the three components overestimates the thrust production. The difference is partially attributed to the drag force on the body caused by the incoming flow. At $ \Gamma m = 15.07$, our simulations of a rigid body at its inflated state with $ e 0 = 0.844$ in the incoming flow reveals that it sustains a drag coefficient $ C d = F x / ( 0.5 \rho V j m 2 D 2 ) = 0.77$, as shown in Fig. 12(c). This value may decrease when the body shrinks as both the surface area and the projected area in the flow direction decrease, and the shape becomes more streamlined (for example, $ C d$ drops to 0.54 when *e* = *e*_{1}= 0.92). Meanwhile, the added-mass-related thrust $ F a$ also contributes to the total thrust. We plot the dimensionless added-mass related thrust coefficient $ C T a = F a / ( 0.5 \rho V j m 2 D 2 )$ (its calculation is described in the next section) at different $ \Gamma m$ in Fig. 10(d). As can be seen, larger body deformation ( $ \Gamma m$) leads to a larger added-mass related thrust. Now the added-mass-related thrust and the viscous drag force generally bridge the gap between the other three thrust components and the total thrust at the constant thrust phase, as shown in Fig. 12(b).

Figure 13(a) presents the time history of the power expenditure coefficient $ C P$. It is found that there is an initial increase in power input, and it later approaches a plateau. In Fig. 13(b), we display the dependences of the time-averaged thrust coefficient and the propulsion factor of impulse generation upon $ \Gamma m$. In line with the previous study (Bi and Zhu, 2020), larger mean thrust production and propulsion factor are obtained at a small maximum equivalent stroke ratio $ \Gamma m$.

### B. Effect of background flow velocity

We first examine the effect of incoming flow on the starting process before the body deformation starts. The time histories of the drag force coefficient at several different incoming flow velocities are shown in Fig. 14(a). It can be found that after the initial phase, the drag force coefficient $ C d$ eventually approaches a constant. The constant drag force coefficient is then plotted against the incoming velocity in Fig. 14(b). It roughly presents a linear correlation with respect to $ U 0 / V j m$. An inspection on the *Z*-vorticity contour, presented in Fig. 15, reveals that stronger vorticity is formed near the nozzle plane of the swimmer at a larger $ U 0$. It is interesting to see how these enhanced vortices from incoming flow would influence the jet flow and thrust production when the body starts to deflate.

To demonstrate the effect of incoming velocity during deformation, we first plot the normalized vortex ring circulation as a function of equivalent stroke ratio $\Gamma $ at various $ U 0$ with $ \Gamma m = 15.07$ in Fig. 16. The maximum circulation decreases with a larger incoming flow velocity. This may be due to the stronger interaction of the incoming flow and the jet flow near the nozzle at larger $ U 0$ (as shown in Fig. 17), which consumes energy from the jet and makes it harder to feed the leading vortex ring. This is reminiscent of the phenomenon that the vortex ring observed in a squid jet is weaker than the formation of a vortex ring in still water (Anderson and Grosenbaugh, 2005). At a larger incoming velocity, more fluid ejected from the chamber is entrained into the vortices near the nozzle along with the background flow, instead of contributing to the growth of the vortex ring. Thus, stronger vortices are formed near the nozzle plane, while trailing vortices behind the leading vortex ring are more suppressed with an increased $ U 0$, as shown in Fig. 18.

Afterward, the effect of the background flow velocity on the pulsed jet propulsion is examined. Figure 19 presents the variation of thrust coefficient, mean thrust and propulsion factor under different background flow velocities. It is shown that they all decrease with increased $ U 0$. To explore the main cause of this total thrust reduction, we first compare $ C T j$, $ C T p$, and $ C T m$ at various $ U 0$ in Figs. 20(a)–20(c). As can be seen, the incoming flow velocity generally has negligible influence on the momentum flux related thrust $ C T j$ and the thrust $ C T m$ sourced from the change rate of the horizontal momentum of the fluid inside the chamber, especially at the quasi-steady phase. It has a more noticeable effect on over pressure-related thrust $ C T p$, mainly during the initial phase of the jetting near $ t = 0.2 \u2009 T d$.

Indeed, these quantities are calculated from the fluid field at the exit plane and inside the chamber, where they are less likely to be affected by the background flow. Even so, background flow with larger incoming velocity enhances both axial and radial velocity at the exit plane in the initial phase of the jetting, as shown in Fig. 21. As a result, we observe an increase in $ C T j$ for a larger $ U 0$ at the corresponding instant. Regarding the radial velocity that comes from body deformation, especially the shrinking in the radial direction, its initial increase may be caused by stronger induction of the recirculation formed by the interaction of background flow and pre-ejected flow near the nozzle. Afterward, the axial velocity at the nozzle gradually plateaus and reaches similar steady values despite the difference in $ U 0$.

Excluding $ C T j$, $ C T p$, and $ C T m$ as the main cause for the overall thrust reduction at high incoming flow velocity, we further examine the added mass related effect. As mentioned earlier, the lateral dimension of the swimmer decreases during the deflation, leading to added-mass related force, $ \u2212 d d t ( m a V ) = \u2212 m \u0307 a V \u2212 m a V \u0307$ with *V* being the speed of the swimmer (for the present tethered mode, we use $ U 0$ instead), in which the term $ F a = \u2212 m \u0307 a V$ contributes positively to thrust generation. By approximating the shape of the swimmer as a prolate ellipsoid, we are able to calculate the instantaneous added mass as $ m a = C a 4 3 \pi \rho a 3 ( 1 \u2212 e \u2032 2 )$ where $ e \u2032 = 1 \u2212 ( 1 \u2212 e 2 + D / L ) 2$. The added mass coefficient $ C a$ is given by $ C a = \alpha \alpha \u2212 2$ with $ \alpha = 2 ( 1 \u2212 e \u2032 2 ) e \u2032 3 ( 1 2 log \u2009 1 + e \u2032 1 \u2212 e \u2032 \u2212 e \u2032 )$ (Imlay, 1961). We plot the dimensionless added-mass related thrust coefficient $ C T a$ in Fig. 20(d). Under the current tethered mode, $ C T a$ becomes larger with increased incoming flow velocity, and it decreases with the body approaching the deflated state. Nevertheless, the increased positive contribution from the added-mass related force at larger $ U 0$ does not overcome the rise of viscous drag. For example, the drag force coefficient $ C d$ at the inflated static body shape is 0.77 when $ U 0 = 0.42 \u2009 V j m$, while it increases significantly to 1.79 at $ U 0 = 0.69 \u2009 V j m$, as shown in Fig. 14(b). Therefore, given the other three thrust components remain almost unchanged and $ F a$ increases, we conclude that the reduction of overall thrust at high incoming flow velocity is caused mostly by viscous drag.

Besides, we are curious about the effect of the incoming velocity upon the unsteady impulse $ I j$ and $ I p$ related to vortex ring formation (see Sec. IV A for the definition) as compared with the overall impulse production. For this purpose, we calculate the ratios between these unsteady impulses and the overall impulse $I$ at different values of $ U 0 / V j m$ shown in Fig. 22. With the increase in incoming flow velocity, the flux momentum and over-pressure associated with unsteady vortex ring dynamics become more crucial to the overall jet thrust production. Although the over-pressure effect at the nozzle plane would enhance the overall impulse, it is relatively small as $ U 0 / V j m$ is varied. In comparison, the flux momentum is of more importance in contributing to the overall thrust generation at a larger incoming flow velocity.

### C. Effect of the jet velocity profile

The time histories of the circulation of the vortex ring for the above two jet velocity patterns are presented in Fig. 24. Throughout the jetting procedure, for the half cosine profile, the circulation *C* with a smaller $ \Gamma m$ is larger than that with a larger $ \Gamma m$ at any instantaneous $\Gamma $. This is because the same amount of fluid mass (thus, the same $\Gamma $) is ejected with higher speed when $ \Gamma m$ is smaller. For the same reason, the maximum circulation for the two jet velocity profiles is smaller than that for a constant jet velocity profile at the same $ \Gamma m$. As suggested by Bi and Zhu (2020), circulation *C* of the vortex ring at a certain instantaneous $\Gamma $ is affected by the evolution of the jet speed before that instant. Therefore, even at the same equivalent stroke ratio $\Gamma $, the circulations for the two profiles may be different. In comparison, when a constant jet velocity is considered, the time histories of circulation *C* for different $ \Gamma m$ are close to each other, as shown in Fig. 7.

Figure 25 shows the time histories of $ C T$ for the two jet velocity profiles. For the cosine jet velocity, a more significant acceleration of jet velocity results in a faster increase in $ C T$ profile before $ t = 0.5 \u2009 T d$ compared with the half cosine profile. After the peak jet velocity is reached, the thrust production falls monotonically before a gradual increase at the late jetting. Negative thrust is yielded at the late jetting phase associated with the deceleration of the jet. In contrast, for the half cosine profile, $ C T$ generally increases throughout the jetting process. In general, the overall thrust production for the constant fashion almost always surpasses those of the other two profiles. The underlying mechanism will be discussed in the following.

To gain further insight into the sources of the overall thrust depicted in Fig. 25 for the two velocity profiles, we compare the time histories of $ C T j$, $ C T p$, $ C T m$, and $ C T a$ in Fig. 26. By comparing Fig. 25(a) with Fig. 26(a), we find that the peak thrust for the cosine profile is dominated by the jet momentum flux out of the nozzle plane. After the peak roughly at $ t = 0.6 \u2009 T d$, $ C T j$, $ C T p$, and $ C T a$ start to decline as the jet velocity decreases. The thrust source $ C T m$ even becomes negative as the deceleration of the velocity of the fluid inside the chamber at the late phase of the jetting, which accounts for the negative total thrust $ C T$ after $ t = 0.8 \u2009 T d$ along with the viscous drag force $ C d$.

Regarding the half cosine profile, all the values of the four thrust components keep positive due to the continuous acceleration of the jet velocity. Generally, the variation patterns of the $ C T j$, $ C T p$, and $ C T a$ are similar for large $ \Gamma m$. The main difference among cases with different $ \Gamma m$ is in $ C T m$, especially the remarkable peak at $ t = 0.5 \u2009 T d$ for $ \Gamma m = 4.70$. This large peak of $ C T m$ is because the maximum jet velocity $ V j m$ is reached faster for a smaller $ \Gamma m$ with the same half cosine profile. It leads to a larger acceleration of the fluid inside the chamber. For the same reason, we see that the value of $ C T m$ at $ \Gamma m = 10.60$ is larger than that at $ \Gamma m = 15.10$. Thus, the small peak of the overall thrust $ C T$ at $ t = 0.5 \u2009 T d$ for $ \Gamma m = 4.70$, as shown in Fig. 25(b), is dominated by $ C T m$. Regarding the maximum total thrust $ C T$ reached at the end of the deflation, they are all attributed to the thrust component $ C T j$ from jet momentum flux.

To interpret the different temporal variations of overall thrust for the three jet velocity patterns shown in Fig. 25, we replot their time histories of $ C T j$, $ C T p$, $ C T m$, and $ C T a$ at $ \Gamma m = 15.10$ in Fig. 27. Not surprisingly, the pattern of $ C T j$ is closely correlated with the specified jet velocity profile. In the constant jet profile, the large overall thrust $ C T$ in the late stage of deflation is attributed to the relatively large jet velocity in that phase. The more distinct peak of $ C T$ for the constant jet profile during the initial phase is closely related to $ C T m$ when the internal flow is quickly accelerated. Indeed, to reach the constant jet profile, the fluid inside the chamber has to experience considerable acceleration from rest to establish the specified jet speed. In contrast, the variations of $ C T m$ for the other two jet velocity profiles are more gentle due to the slow change of $ V j$. The maximum $ C T p$ are close for the three jet velocity profiles, but located at different time instants. Regarding the added-mass related thrust $ F a$, although different variation patterns of $ C T a$ are presented for the three jet velocity profiles, it plays a non-significant role in the overall thrust generation. Their contribution may increase if large body deformation is considered.

## V. CONCLUSIONS

Using a three-dimensional, unsteady, viscous and compressible Navier-Stokes fluid solver based on cell-centered finite volume method, we numerically study squid-like jet propulsion through body deformation. The cephalopod-inspired swimmer is idealized as an ellipsoid-like pressure chamber with an opening serving as the outlet of jet flow. With the increase in the eccentricity of the body, the chamber deflates so that internal fluid is ejected. We focus on a single deflation deformation and systematically investigate the near-body fluid field, vortex ring formation and jet propulsion performance with various background flow conditions and different jet velocity profiles.

For a constant jet velocity profile, the flow field visualization shows that no pronounced secondary vortices are observed behind the leading vortex ring at a relatively small jet-based Reynold number $ Re j$ considered in this study, which is different from the high Reynolds number cases. The incoming flow speed $ U 0$ of the background flow plays a role in the formation of the vortex ring. Larger $ U 0$ leads to stronger interaction between the background flow and jet flow near the nozzle. A pronounced recirculation zone is formed near the exit plane as the jet flow is entrained into these vortices, making it harder to feed the vortex ring and tailing vortices so that circulation of the vortex ring is reduced at high incoming velocity. The jet velocity profile is another significant factor determining the evolution of the vortex ring. In comparison with the cosine and half-cosine jet velocity profiles, a constant jet velocity profile produces larger peak values of *C* at a given equivalent stroke ratio. This is explained by the fact that the fluid is initially ejected at a higher speed under a constant jet velocity profile without an acceleration phase.

Based on conservation of momentum, three distinctive thrust sources are identified, i.e., momentum flux related thrust $ C T j$, excessive pressure (at the nozzle exit) related thrust $ C T p$, and thrust $ C T m$ due to time variation of the momentum of the fluid inside the pressure chamber. For a constant jet velocity pattern, total thrust peak is observed in the initial phase of the jetting and mostly dominated by $ C T m$, while $ C T j$ determines the thrust at the quasi-steady phase. For the cosine jet velocity profile, the jet momentum flux related thrust $ C T j$ accounts for the overall peak thrust near the middle of the jetting phase. Regarding the thrust production for the half cosine profile, $ C T$ generally experience a continuous increase throughout the jetting mainly attributed to the growing $ C T j$. By comparing thrust production for the three jet velocity profiles, we find that the thrust $ C T j$ from jet momentum flux, which is closely related to instantaneous jet speed, is essential to maintain a positive total thrust. For example, with the deceleration of the jet during the second half of the deflation period for a cosine jet velocity profile, the total thrust becomes negative. In terms of $ C T$, the constant profile outperforms the other two profiles during the later stage of deflation quasi-steady thrust phase due to larger jet velocity (thus, greater $ C T j$). Our evaluation of added-mass related thrust $ F a$ due to the body deformation suggests that it contributes a non-negligible portion in the overall thrust for all the three jet velocity profiles for the body kinematics considered here.

Although larger incoming flow velocity $ U 0$ would result in thrust reduction, the background flow is proven to have a negligible effect on momentum flux related thrust $ C T j$ and thrust $ C T m$ from the horizontal momentum change of the fluid inside the chamber. Increased incoming flow velocity results in greater added mass related thrust $ F a$ attributed to body deformation. Nevertheless, the rise of $ F a$ is overweighted by the significantly increased viscous drag force at larger $ U 0$. Therefore, the overall thrust declines at large incoming flow speed. We also find that larger background flow velocity would increase the portion of $ C T j$ and $ C T p$ in the overall thrust generation.

Finally, we note that a three-dimensional squid-like swimmer model is introduced in this study. It is then applied to a relatively simple scenario that only focuses on single-cycle deflation deformation. More sophisticated cases involving three-dimensional effect, e.g., symmetry-breaking instability demonstrated in 2D studies, will be explored in future work.

## ACKNOWLEDGMENTS

This research was supported by an EPSRC Supergen ORE Hub Flexible Fund Program Grant “Autonomous Biomimetic Robot-fish for Offshore Wind Farm Inspection” (No. EP/S000747/1). This work used the Cirrus UK National Tier-2 HPC Service at EPCC (http://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC (No. EP/P020267/1) and ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk) based at the University of Strathclyde. Yang Luo thanks China Scholarship Council (CSC) for financial support during his study in the UK.

## DATA AVAILABILITY

The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.15129/f47dd7f5-efaa-4c1f-9b89-22e167322cc9.

## REFERENCES

*Lolliguncula brevis*: Evidence of multiple jetmodes' and their implications for propulsive efficiency

*Loligo opalescens*

*The Complete Expressions for Added Mass of a Rigid Body Moving in an Ideal Fluid*

*Parallel Computation of Three-Dimensional Aeroelastic Fluid-Structure Interaction*