Numerical study on the aerodynamic performance of the three-dimensional wing of a jellyfish-like flyer

The jellyfish-like flying machine is a new development direction of the future bionic flapping-wing aircraft besides the insect-mimic and bird-mimic micro air vehicles (MAVs). To better understand the underlying fluid mechanisms of the jellyfish-like flyer, we numerically simulated the aerodynamic forces of the three-dimensional flapping wings under different control parameters. The effects of flapping amplitude, vortex wake, up-flight speed, and wing–wing interaction on aerodynamic performance were investigated. The results show that, at hovering, the mean lift rises rapidly at first and then tends to be stable with the increase in flapping amplitude. The vortex wake can improve the lift at large flapping amplitudes, while it reduces the lift at very small flapping amplitudes. With the increase in up-flight speed, the lift decreases. However, the sources of lift reduction are different for different flapping amplitudes. When the two wings flap together and the distance between the wings is small enough, the wing–wing interaction can improve the lift by about 15% compared with that of a single wing, but much higher power is required, resulting in lower efficiency. The results of this study provide new insights into the flight mechanism of the jellyfish-like aircraft and have important guiding significance for the design and optimization of the jellyfish-like flying machine.


Introduction
After hundreds of millions of years of evolution, flying creatures can obtain astonishing aerodynamic performance by flapping wings, which has become the inspiration for designing bionic aircraft. The design of bionic aircraft is generally based on two typical types of motion: the lift-based stroke in a horizontal plane (also termed the normal mode) [1], like that of the hummingbird, and the drag-based paddling stroke in a vertical plane [2], like that of the dragonfly. Since the conventional steady-state aerodynamic theory is insufficient to explain the high aerodynamic force coefficients required by these small creatures [3], much endeavor has been put into obtaining the unsteady aerodynamic interpretations. The revealed unsteady aerodynamic mechanisms include delayed stall, rapid pitching rotation, rapid acceleration (or added mass), clap and fling, and wake capture [4][5][6][7]. Although the flapping-powered robotic flyer shows great potential in producing high aerodynamic forces and energy efficiency relying on these unsteady mechanisms, it also needs to be equipped with feedback control systems to maintain stable flight due to its inherent instability [8][9][10][11].
In an experimental study, Childress et al. [12] observed that a pyramid-shaped flexible body with an array of four wings can achieve a stable hovering within a vertically oscillating airflow. Motivated by this experiment, Ristroph and Childress [13] developed a prototype of the jellyfish-like flying machine. This flapping-wing air vehicle can achieve self-righting hovering by opening and closing four downward-pointing wings, without relying on additional aerodynamic surfaces and without feedback control. Each wing of the aircraft was hinged on the upper horizontal loop and rotates around a pivot point, which is very similar to the simple pitching motion. This provides a new way of movement for designing bionic aircraft.
This intrinsically stable jellyfish-like flying machine and the relatively simple actuation motion attracted the attention of researchers, prompting them to carry out more detailed research on the aerodynamics and flows of such flyers. Using twodimensional (2D) models, Fang et al. [14] and Zhang et al. [15] explored the aerodynamics and flight dynamics of the flying machine with theoretical and numerical methods, respectively. Both studies found that pitching amplitude and wing separation played a determining role in the locomotion state and aerodynamic performance. In addition, Zhang et al. [16] simulated the passive flight of a pair of 2D flexible pitching foils with a concave-down configuration in the vertically oscillating airflow. It was shown that flexibility can strongly affect the weight-supporting capability and stability of the flyer.
It is worth noting that all the aforementioned studies that were carried out on the jellyfish-like aircraft are 2D studies. More recently, Liu et al. [17] experimentally investigated the effect of the time asymmetric pitching motion on the lift generation and the flow-field characteristics of this kind of ornithopters, which is the only threedimensional (3D) research on hovering. They found that the fast downstroke and slow upstroke pitching pattern was superior to symmetric pitching. Although the flow-field was 3D in their study, the vortex structures were still displayed by the 2D PIV measurement. The evolution of the side-edge vortex was not considered, and its effect on the forces of the wing remains unclear so far. The experimental study by Buchholz and Smits [18] demonstrated that the streamwise vorticity formed at the 3D panel edges tended to dominate the wake formation. Other studies [19][20][21] have also demonstrated that the 3D effect plays a key role in the formation of vortices and the aerodynamic performance of the flapping wings. Therefore, with respect to this jellyfish-like flyer, it is still far from being possible to fully understand the details of 3D vortex evolution and aerodynamics, especially the vortex structures around all edges of the wing.
In contrast to the previous 2D studies, this paper numerically investigated the aerodynamics of 3D flapping wings employing the locomotion of the jellyfish-like flyer. The hovering was emphasized because it is an indispensable flight mode for MAVs to perform special tasks. This is also different from most other works, which studied the pitching wing in an incoming flow [18][19][20][21]. The aerodynamic performance under some key motion parameters was examined and related to the flow fields. In addition, the wing-wake interaction was also considered, which was neglected in previous studies. The remainder of the paper is organized as follows: the wing motion and the numerical method with its validation are described in Section 2. The simulated results are presented in Section 3. In Section 3.1, the influence of flapping amplitude on the aerodynamics of a single wing under hovering conditions is studied first, and then the vortex structures are examined in detail. The role of vortex wake is also emphasized in Section 3.2. Next, the effect of up-flight speed on aerodynamic performance is investigated in Section 3.3. Finally, in Section 3.4, the interference effect between wings is studied and comparisons with relevant biological data are presented. In Section 4, we summarize the conclusions and discuss future research avenues.

Wing model and motion
In this work, the wing was modeled as a flat plate with rounded edges. The wing shape is similar to the teardrop-shaped wing of the jellyfish-like aircraft manufactured by Ristroph and Childress [13] ( Figure 1A). The aspect ratio R/c of the wing is 2.25, where R and c are the wing length and mean chord length, respectively. The radius of gyration r 2 of the wing is 0.62R.
To facilitate the flow-field analysis, one opposing pair of wings was used instead of the bell-like configuration of four pieces of wings in the original prototype of the ornithopter. It should also be noted that the simple drive mechanism shown in Figure 1A does not close all four wings simultaneously but rather causes one opposing pair to lead the other by a quarter period. Figure 1B illustrates the pitching motion. The two wings flap about the wing root with a mirror-symmetric sinusoidal pitching motion in time. The origin of the Cartesian coordinates OXYZ coincides with the pivot of the right wing, with the Y-axis passing through the two wing roots and the Z-axis pointing vertically upward ( Figure 1B). The wings can move upward at a speed of V, and V 0 represents the hovering state. The separation distance between the two wings' hinge points is H. θ and Φ represent the opening angle and peak-to-peak flapping amplitude, respectively (see Figure 1B). The pitching angles of the right (ϕ 1 ) and left (ϕ 2 ) wings are expressed as follows: where n is the wingbeat frequency and t is the time. A dimensionless time (τ) is defined as that varies between 0 and 1 for each flapping cycle. One flapping cycle includes two strokes, with τ 0 and τ 1 representing the start of the first right-to-left stroke and the end of the second left-to-right stroke, respectively. This paper first studied the aerodynamics of a single wing (the right wing) at different flapping amplitudes and up-flight speeds, assuming that the wing-wing spacing is large enough and there is no interference between the left and right wings. The opening angle for the single wing was set at θ 90°to ensure that the resultant force of the wing points vertically upward [22]. Then, the interference effect of two wings at different distances was studied. It should be noted Frontiers in Physics frontiersin.org that to avoid collision between wings, the opening angle was no longer 90°when the wing-wing interaction was considered.

Numerical methods
The flows and aerodynamic forces were computed using the method of computational fluid dynamics (CFD), and the laminar flow model was used due to the low Reynolds number. The governing equations of the flow are the three-dimensional incompressible unsteady Navier-Stokes (N-S) equations: In the aforementioned equations, u, t, and p represent dimensional fluid velocity, time, and pressure, respectively. ∇ is the gradient operator, and ∇ 2 is the Laplacian operator. v and ρ are the kinematics viscosity and density of the fluid, respectively. In the non-dimensionalization, U, c, and c/U are taken as reference velocity, length, and time, respectively. U is the mean flapping velocity at the radius of gyration of the wing, U 2Φnr 2 . Thus, the non-dimensional parameters involve: u* u U , t* tU c , p* p ρU 2 , z zt* c U z zt , ∇* c∇. The non-dimensional wingbeat cycle is written as T* 2Φr 2 /c, which is proportional to the peak-to-peak amplitude. Then, the dimensionless form of the N-S equations can be deduced as follows: zu* zt* + u* · ∇*u* −∇p* + 1 Re ∇* 2 u*.
Here, Re is the Reynolds number, and Re cU/v. Its value is set at 200 in the following simulations. The method for solving the N-S equations is the same as that in the previous studies [23,24]. The algorithm is based on the method of artificial compressibility developed by Rogers et al. [25] and Rogers and Pulliam [26]. To discretize the time derivative of the momentum equations, a three-point backward differencing formula with second-order accuracy is used. Moreover, a pseudo-time and a pseudo-time derivative of pressure are introduced into the continuity equation. The resulting system of equations is iterated in pseudo-time. As long as the pressure derivative in pseudo-time approaches zero, the velocity divergence at the new time level approaches zero. The derivatives of the viscous fluxes in the momentum equation are approximated using second-order central differences. The upwind differencing method based on the flux-difference splitting technique is used for the derivatives of convective fluxes. For points located close to boundaries, the second-order upwind differencing is employed; for points located toward the interior, the third-order upwind differencing is used. This algorithm was described in detail in [25,26]. Under the far-field boundary conditions, at the inflow boundary, the relative velocity at the boundary is used to determine the velocity components, while  [13]. The prototype pulls in and pushes out four wings by using a motor (M) to rotate a crankshaft (CS) that connects via a link (L*) to each wing. The body of the flyer consists of two crossed vertical loops that support the motor below and a horizontal upper loop. More detailed information about this flyer can be found in [13]. (B) Simplified wing-pair model. Here, V is the up-flight speed; θ represents the wings' base opening angle; R and H denote the wing length and the distance between the two pivot points, respectively. Two rigid plates are driven to flap with prescribed mirror-symmetric sinusoidal pitching motions of peak-to-peak amplitude Φ. Frontiers in Physics frontiersin.org pressure is extrapolated from the interior; at the outflow boundary, the velocity is extrapolated from the interior, and pressure is set to be equal to the static pressure of the still air. On the wing surfaces, impermeable wall and no-slip boundary conditions are applied, and the pressure on the boundary is obtained through the normal component of the momentum equation written in the moving coordinate system. To facilitate the study of the wing-wing interaction, the technique of moving overset grids was employed. Figures 1C, D show the grid system that consists of two body-fitted curvilinear wing grids that extend a relatively short distance from the wing surfaces and a Cartesian background grid that extends to the domain's far-field boundaries. The wing grid captures features such as boundary layers, separated vortices, and vortex-wing interactions, and the background grid captures the vortex wake and carries the solution to the far field. Data are interpolated from one grid to another at the inter-grid boundary points using tri-linear interpolation. The wing grid (O-H type) was created using a Poisson solver, which was based on the work of Hilgenstock [27]. It allows fast generation of the threedimensional wing grid with full boundary control. The background grid was generated algebraically. The nondimensional time step (Δt*) is equal to T*/440. The wing grid has dimensions 43 × 71 × 70 in the normal direction, around the wing, and in the spanwise direction, respectively (the first grid layer thickness is 0.001c). The outer boundary of the wing grid is 2.5c from the wing spanwise axis in the normal direction and 1.5c from the wing root and tip in the spanwise direction. The background grid has dimensions 121 × 121 × 121 in the X, Y, and Z directions, respectively, and its boundaries are 20c away from the wing root. To accurately capture the vortex structure, the background grid is refined in the regions covered by the wing grid motion and the main development region of the vortex wake. In these regions, the background grid is evenly distributed, and the grid spacing is approximately the same as that of the outer part of the wing grid. A detailed study of the variables such as time step, domain size, and grid dimensions was conducted, and the values mentioned previously are appropriate for the present study.
After the pressure and velocity of the fluid are obtained, the aerodynamic forces are calculated by integrating the pressure and the viscous stresses on the wing surface. The lift (L) and drag force (D) of the wing are defined as the force components along the Z and Y axes, respectively. The pressure, lift, and drag coefficients are defined as C p (p − p ∞ )/(0.5ρ(U) 2 ), C L L/(0.5ρ(U) 2 S), and C D D/(0.5ρ(U) 2 S), respectively. Here, p is the pressure, p ∞ is the pressure at infinity, and S is the area of one wing. The aerodynamic power (P a ) required to overcome the aerodynamic moments is determined by the aerodynamic torque and angular velocity of the wing. The aerodynamic power coefficient is calculated as C P,a P a /(0.5ρ(U) 3 S). The propulsion efficiency η is given as η C L /C P,a . C L and C P,a are cycle-averaged C L and C P,a , respectively.
The code has been verified by comparing the calculated forces with those of experimental measurement of flapping insect wings [28]. It showed that the calculated results matched the experimental measurements closely. To further ensure that our code is correct, Figure 2 compares our calculations with the work of Li and Lu [29], which studied the forces of pitching and heaving flat plates in an incoming flow. Again, the aerodynamic data are in good agreement. In addition, according to the previous research of our group [24], this flow solver can accurately predict the wing-wing or wing-wake interactions.
3 Results and discussion 3.1 Effect of flapping amplitude As mentioned in Section 2.2, the non-dimensional wingbeat period (T* 2Φr 2 /c) is proportional to the peak-to-peak amplitude Φ. Thus, the mean angular velocity is the same for different amplitudes. The main purpose of changing Φ is to compare the difference in aerodynamic performance between smaller, faster motions and larger, slower motions.
In this section, the aerodynamics and flow structures of a single wing under the condition of hovering are simulated first because hovering is an important flight mode for the MAVs to perform specific tasks. The amplitude of pitching is varied from Φ 10°to 80°at intervals of 10°. It should be noted that the range of pitching amplitude investigated in this study is larger than those in previous works [18 -21]. The transient lift and drag FIGURE 2 Comparison of the present instantaneous lift coefficient (A) and drag coefficient (B) with previous data obtained by Li and Lu [29] for a rectangular plate with an aspect ratio of 0.5 at Re 500.
Frontiers in Physics frontiersin.org coefficients of the wing after the flow reached a stable periodic state are shown in Figure 3. It is seen that there are two positive crests and two negative troughs of lift in one cycle for each curve ( Figure 3A). The magnitude of positive crests is larger than that of the negative troughs, which ensures that the wing is subjected to an upward average lift. With the increase in Φ, the peak of lift increases and shifts slightly rightward in phase. The curve of drag for each flapping amplitude has only one crest and one trough with equal magnitude, which results in zero cycle-averaged drag ( Figure 3B). As shown in Figure 3B, the crest and trough of the drag increase as Φ decreases. C L , C P,a , and η of the wing as a function of Φ are shown in Figure 4. It shows that C L increases with the enlargement of Φ at first, and then it remains approximately stable when Φ is larger than 30°. Since C P,a decreases gradually, η shows a general trend of continuous growth as Φ increases. It indicates that a larger, slower flapping motion can improve the efficiency of the wing compared to a smaller, faster flapping motion. This is also consistent with the findings of Van Buren et al. [30] in the study of a heaving and pitching wing in an incoming flow.
Next, we further discuss the vortex structures around the wing under different flapping amplitudes. Figure 5 illustrates the evolution of the vortex structures around the wing at eight typical instants during one flapping cycle at Φ 40°. The vortex is visualized by the Q criterion, and the oblique view is used to facilitate the identification of different vortices. The vortex structures at τ 0 (the beginning of the right-to-left stroke) and τ 0.5 (the beginning of the left-to-right stroke) in Figure 5 can be viewed as mirror images of each other due to the symmetry in the wing kinematics. The "mirrored image" is of help to clearly demonstrate the vortex structures around the wing.
As shown in Figure 5, when the wing starts to flap from right to left (τ 0), the main vortex structures are two curved vortex rings. One is composed of the side-edge vortices (SV 2 ) and the wingtip vortex (TV 2 ) formed in the preceding left-to-right stroke, and the other is composed of SV 1 and TV 1 formed in the earlier right-to-left  gradually generated at the wing edge, forming a new vortex ring. At this instant, the wing as well as the newly formed SV 3 bypasses SV 2 , providing an opportunity for the wing-wake interaction. At τ 0.25

FIGURE 5
Formation and evolution of vortices at eight typical instants (A-H) during one flapping cycle for a single wing at Φ 40°. In all snapshots, the isosurface of the Q criterion with Q 2 is shown.

FIGURE 6
Instantaneous vortex structure of the wing for different amplitudes at Φ (A-H) τ 0.5. In all snapshots, the iso-surface of the Q criterion with Q 2 is shown.
Frontiers in Physics frontiersin.org 06 ( Figure 5C), SV 3 and TV 3 become stronger, and SV 2 is distorted due to the induction of the newly generated vortex ring. From τ 0.375 to 0.5, TV 3 detaches away from the wing, resulting in a gradually bent vortex ring. In addition, the vortex ring formed by SV 2 and TV 2 has been seriously distorted, mainly because SV 2 is further entrained and induced by TV 3 . It should be noted that SV 1 and TV 1 gradually dissipate from τ 0 to 0.5 and finally disappear in the following process. In the subsequent stroke (from τ 0.5 to 0.875), the formation and evolution of vortex structures are similar to those in the preceding stroke (from τ 0 to 0.375), that is, the formation of SV 4 and TV 4 and their entraining and inducing effect on SV 3 . By 2D PIV measurement, Liu et al. [17] showed that the flow field for the wing of the jellyfish-like flyer is characterized as the consecutive shedding of two vortices at the wingtip with an opposite rotation sense, which are generated during the downstroke and upstroke, respectively. This is similar to the generation and shedding of TV 3 and TV 4 in the present study, although the Reynolds numbers are different in the two studies.
The vortex structures under different flapping amplitudes at τ 0.5 are shown in Figure 6. It should be noted that the wing begins to flap from left to right at this moment. When Φ is smaller than 40°, the vortex structures become more complex, although they show similar topology. It is seen that SV 3 and TV 3 become larger gradually with the increase in Φ (Figures 6A-D). Due to the short dissipation time, the vortex structures generated by the previous strokes for small flapping amplitude can still be displayed. When Φ is larger than 40°, the strength of the side vortices changes little, which is consistent with the cycle-averaged lift force performance in Figure 4. Because of the long dissipation time at the large flapping amplitude, SV 2 and TV 2 , which can be observed for Φ 40°, have almost disappeared for Φ 70°and 80°.
The iso-surfaces of time-averaged velocity u Z in the flow field for different Φ cases are shown in Figure 7. The iso-surfaces in gray with a threshold of u Z 0.2 indicate the unfavorable upwash for propulsion, while those in blue with u Z −0.3 represent the favorable downwash for propulsion. The upwash is supposed to be generated by the part of the distorted vortex ring (SV 2 ), which is labeled in Figure 6. As shown in Figure 7, the downwash concentrates near the wing at hovering. This is different from the case with incoming flow, in which the wakes developed into bifurcated jets [21,31]. The downwash is distributed under and on both sides of the wing, while the upwash is mainly concentrated on both sides of the wing. As Φ increases from 10°to 40°, the downwash region becomes significantly wider in the Y-axis direction and longer in the Z-axis direction (Figures 7A-D). The longer downwash in the vertical direction indicates a larger lift. As Φ increases further, the width of the downwash keeps increasing, but the length changes little ( Figures 7E-H), so the overall downwash beneficial for lift generation is basically unchanged. This is consistent with the results in Figure 4, which shows that the mean lift increases rapidly first and then stabilizes approximately. Meanwhile, the upwash on both sides of the wing decreases continuously. As Φ reaches 80°, the upwash has nearly disappeared. This is because the distorted part of the vortex wake (SV 2 ) dissipates faster for larger Φ, which can also be observed in the circled part in Figure 6.
According to the vortex structures and downwash in Figure 6 and Figure 7, it can be summarized that with the increase in the flapping amplitude from 10°to 40°, SV and TV on the wing surface are significantly enhanced, forming a stronger vortex ring to induce stronger downwash, thus causing a larger lift. As the flapping amplitude increases further, the vortex wake, as well as the accompanying downwash, mainly widens in the horizontal direction, but the strength of the vortex and the downwash in the vertical direction changes little, resulting in basically no change in the lift.

Wake effect
For flapping wings at hovering, the vortices are convected away from the wing at a very small speed due to the lack of incoming flow. Hence, as a wing decelerates and reverses to start the next stroke, the wing inevitably captures and interacts with the remaining wake from the previous stroke. It is necessary to explore the wake effect under different amplitudes at hovering, which has not been examined for the jellyfish-like wing motion in previous studies. The contribution of wake capture to the aerodynamic force is isolated by comparing the averaged lift force of the first stroke when the wing starts to flap in a quiescent fluid and the wake effect is absent, with that of the ninth stroke after the periodic flow field has been established. The results are shown in Figure 8. It is seen that with the increase in Φ, the effect of wake capture changes from reducing lift to increasing lift, and the improvement of lift keeps increasing (see the solid and dashed black lines with delta symbols). Also displayed in Figure 8 are the two components of lift contributed by friction and pressure, respectively. For all Φ considered, whether there is wake or not, the friction-based lift (denoted as C L(F) ) is always negative, while the pressure-based lift (denoted as C L(P) ) is always positive. The wake effect is to increase C L(P) and decrease C L(F) . As Φ increases, the improvement in C L(P) caused by the wake effect becomes larger (compare the solid and dashed red lines with diamond symbols), while the reduction in C L(F) due to the wake effect becomes smaller (compare the solid and dashed red lines with gradient symbols). As a result, the combined effect of the wake becomes larger with the increase in Φ. It should be noted that at Φ 10°, the wake effect on C L(P) is negligibly small, but its reduction effect on C L(F) is much larger, so the wake effect is to decrease C L at this small flapping amplitude.
Next, we analyzed the case of Φ 40°in detail to explore the reason why the pressure-based lift was enlarged by the vortex wake. Figure 9 shows the spanwise component of vorticity (ω S ) and pressure contours of the wing at the 0.7R section at τ 0.125 for the first stroke and the ninth stroke. The vortex structures and the location of the section are also displayed in the figure to assist the analysis. Compared with the result in the first stroke, not only the positive pressure on the lower wing surface in the ninth stroke is larger, but the negative pressure on the upper wing surface is also lower, which makes the pressurebased lift larger (see Figures 9A, B).
For the first stroke, only SV 3 is generated at this time instant (Figures 9A, C). In contrast, for the ninth stroke in Figures 9B, D, besides the newly generated SV 3 on the upper wing surface, two edge vortices (SV 2 ) can be seen on the lower surface of the wing, and they are part of the vortex ring generated in the preceding stroke ( Figure 9D). The two edge vortices (SV 2 ) or the previous vortex ring can induce favorable fluid velocity toward the lower surface of the wing, resulting in higher pressure as shown in Figure 9B. Furthermore, it is noted that on both edges of the wing, the vorticity signs of the newly generated vortex and wake vortex are opposite ( Figure 9B). Therefore, as the wing approaches the vortex wake, narrow channels are formed between SV 3 and SV 2 , which will accelerate the fluid in the middle and promote the strength of SV 3 at the same time. As a result, larger side-edge vortices (SV 3 ) are generated on the upper wing surface, resulting in lower pressure. This wake capture effect is similar to the findings of Birch and Dickinson [32], and Lua et al. [33]. These findings indicated that the wake capture can improve the lift of insect wings. However, it is somewhat different from the study of Wu and Sun [34], which shows that the vortex wake can reduce the lift of the wing by generating a downwash.
It has been shown in Figure 6 that with the decrease in flapping amplitude, the strength of the vortex gets weaker. At Φ 10°, SV 2 and TV 2 generated in the preceding stroke are very small ( Figure 6A), so the surface pressure profit is small. Moreover, because the wing is more vertically placed, the contribution of the change in surface pressure to the lift can be almost ignored. At larger Φ, not only does the effective horizontal projection area of the wing increase, but the stronger SV 2 and TV 2 also prompt more FIGURE 8 Variation in the mean lift and the pressure-and friction-based lift of the wing as a function of the flapping amplitude Φ with/without wake effect. Subscripts "P" and "F" in the symbols represent the pressure-and friction-based components of the lift, respectively. Subscripts "W" and "N" in the symbols represent the cases with and without wake effect, respectively.

FIGURE 9
Pressure and spanwise vorticity (ω S ) contours at the 0.7 R wing section, as well as the instantaneous vortex structure (Q 2) of the wing at t*/T* 0.125 (A,C) and t*/T* 4.125 (B,D) (solid and broken lines indicate positive and negative vorticity, respectively).

Frontiers in Physics
frontiersin.org favorable changes in surface pressure, so the effect of the wake capture on enlarging the C L(P) becomes larger. The reason for the negative C L(F) is easy to understand. The flapping wing obtains upward reaction force by driving the fluid to convect downward. The downward movement of fluid causes a velocity gradient near the wing surface, which can result in the negative lift. When there is a vortex wake, the downward velocity will be enhanced by the downwash, leading to an increase in the friction force. Thus, the negative friction-based lift becomes smaller when there is a vortex wake.

The effect of up-flight speed
By changing the flapping amplitude, the jellyfish-like aircraft can switch between hovering, ascending, and forward flight states [13]. In this section, we further explored the aerodynamic performance of a single wing at different up-flight speeds (V, see Figure 1B). Figure 10 plots the mean lift coefficient C L of the wing flying upward at V 0, 0.25V 0 , 0.50V 0 , 0.75V 0, and 1.00V 0 , respectively. It should be noted that V 0 is defined as the mean velocity of the wingtip at Φ 40°. It is seen that for a specific ascending speed, C L rises rapidly at first and then tends to be stable with the increase in flapping amplitude. When the flapping amplitude is fixed, C L decreases gradually with the increase in up-flight speed. The smaller the flapping amplitude is, the faster the C L drops, and even negative average lift values appear in some cases.
The vortex structures of the wing with different up-flight speeds at τ 0 for three cases (Φ 10°, 40°, and 70°) are shown in Figure 11. To show more vortex structures, the criterion of Q 0.5 is used in these images. For Φ 10°(left column in Figure 11), the wake under the wing consists of two sets of vortex rings (shed a vortex ring in each half-cycle). This is similar to the work of Buchholz and Smits [18] in which a rectangular plate conducting pitching motion was studied and a series of "horseshoe vortices" was observed in the wake. With the increase in the up-flight speed, the vortex street becomes longer. As the vortices convect downward and dissipate gradually, the vortex ring breaks and forms two separate branches. When Φ increases to 40°, the vortex rings become larger, and they are no longer concentrated in the vertical direction but convect downstream at an oblique angle relative to the wake centerline. For Φ 70°, the vortex rings are further enhanced and the inclination angle becomes larger. Figure 12 plots the instantaneous lift coefficients and its two components (C L , C L(P) , and C L(F) ) of the wing for Φ 10°, 40°, and 70°cases at different up-flight speeds. At Φ 10°, the change in C L(P) with the increase in V is relatively small ( Figure 12B), but C L(F) obviously decreases ( Figure 12C). This leads to the rapid decrease in C L with the up-flight speed for the Φ 10°case, as shown in Figure 12A. The study by Floryan et al.
[35] on the heaving or pitching foils also showed that viscous drag added an approximately constant negative offset to the thrust coefficient. Our results show that this constant negative offset increases rapidly with the increase in the up-flight speed at a very small flapping amplitude. Considering that the flapping wing is approximately vertical for the Φ 10°case, the larger up-flight speed will cause a bigger velocity gradient in the normal direction at the wing surface, resulting in a greater frictional force.
When Φ increases to 40°and 70°, C L(F) is essentially unchanged at different up-flight speeds ( Figures 12F, I), while C L(P) decreases ( Figures 12E, H), resulting in the reduction of C L with the increase in V. The change in C L(P) is caused by the velocity change in the incoming flow relative to the wing surface. When the wing moves up at a certain speed, the lower wing surface will face a lower relative inflow velocity as the wing flaps downward, while the upper wing surface will face a higher relative inflow velocity as the wing flaps upward. This opposite change in the relative inflow velocity can reduce the magnitude of positive C L(P) at about τ 0.125 and increase the magnitude of negative C L(P) at around τ 0.375 ( Figures 12E, H). Since the effective horizontal projection area of the wing also becomes larger following the increase in the flapping amplitude, the reduction in C L(P) would be further amplified at larger Φ.

The effect of wing-wing interaction
Using the wing-pair model in Section 2.1, we explored the effect of the aerodynamic interaction between two wings at different separations (H, see Figure 1B) at the hovering state. The opening angle θ 80°and the flapping amplitude Φ 20°were selected to ensure that the two wings would parallel with each other at the end of the downstroke. Here, the process of two wings approaching is called downstroke, and the process of moving apart is called upstroke. Figure 13 displays the instantaneous lift and drag coefficients of the right wing at different wing-wing spacings, along with those of a single wing for comparison. It is seen that the magnitudes of the peak and trough of lift and drag of the wing pair are larger than those of the single wing, and they increase with the decrease in H. Since the opening angle is now 80°, the setting of the relative position of the pair of wings in this section is similar to that of Liu et al. [17] Therefore, the change trends of the forces in both studies for the symmetric pitching motion are also similar, that is, the downstroke mainly produces positive lift and the upstroke produces negative lift (see Figure 3A in Ref. [17]). The flow-field information of two wing-pair cases (H 0.1 c and H 1.2 c) and the single wing case was selected for analysis and comparison. The pressure distributions and the velocity vectors in the YOZ plane at τ 0.125 and τ 0.75 are plotted in Figure 14. At H 0.1 c, it is seen that as the wings move close to each other, the fluid between them is "squeezed" and the pressure significantly increases compared with that of the single wing ( Figures 14A, C), indicating that the larger force is caused by the squeezed air between the two wings. During the upstroke, the "moving apart" of the wings in the opposite direction increases the gap between the wings and creates a larger suction pressure there ( Figures 14B, D). This suction pressure produces a larger negative lift and drag than that of the single wing to keep the wings from moving apart. The large aerodynamic forces due to the wing-wing interaction are also similar to those of wings conducting clap and fling motion [36,37]. When the separation increases to H 1.2 c (Figures 14E, F), the pressure and the velocity vectors of the wing pair are not very different from those of the single wing, leading to almost the same aerodynamic force of the two cases. It should be noted that a favorable strong downward jet is generated under the two wings at H 0.1 c (Figure 14C), while no jet but downwash is generated in the other cases. Therefore, we can deduce that when multiple flapping wings are arranged at the distance of the prototype of the jellyfish-like aircraft, which is about 1.2c, the wing-wing interaction is very small. The variation in C L , C P,a , η, and the mean absolute value of the drag | C D | of the right wing as a function of H is displayed in Figure 15. It shows that C L , C P,a , and | C D | increases with the decrease in H, but η has the opposite changing trend. At H 1.2c, the average lift is almost the same as that of a single wing, which indicates that the interference effect between the two wings becomes negligibly small when H reaches 1.2 c or larger. Compared with the single wing, C L and | C D | of the wing at H 0.1 c are increased by 15.5% and 24.7%, respectively, but at the same time, C P,a is enlarged by 66.3%, resulting in a 30.1% decrease in efficiency. Although reducing the distance between wings can improve the lift to a certain extent, it requires more energy consumption and becomes less efficient, which is not beneficial to the economic flight of this jellyfish-like flyer. Dewey et al.
[38] found a thrust enhancement  of about 70% for the out-of-phase 2D foils in a side-by-side configuration compared to an isolated foil with no observable change in the propulsive efficiency. Fang et al. [14] found that efficiency increased as the wings get closer due to a mirror-image "ground effect" between the wings. The differences in the results between our current research and the previous two research studies may lie in the difference in wing models (2D vs 3D), and/or the incoming flow settings (with vs without incoming flow).
It is instructive to test our calculation results against biological observations. For real jellyfish, there are two main mechanisms for propulsion [39]. One is the jet propulsion mechanism used by the elongated jellyfish. During the contraction stroke of the elongated jellyfish, it drives the entire bell and creates a powerful jet in the subumbrellar cavity. Periodic bell contractions cause the ejected jets to form a series of uniform vortex rings with the same rotational direction in the wake. This motion can produce rapid acceleration and swimming speed but require high energy expenditure. The other is the "rowing mechanism" used by the oblate jellyfish [40]. For oblate jellyfish, contractions occur primarily at the bell margin, which is similar to using paddles (the bell margin) to row a boat (the

FIGURE 15
Variation in C L , C P,a , η, and the mean absolute value of the drag | C D | of the wing as a function of wing-wing separations (H). The symbol "∞" indicates the single wing case.
Frontiers in Physics frontiersin.org body). An oblate jellyfish produced periodic vortex wake instead of strong jets. A train of closely spaced vortex ring pairs was formed in the wake and convected downstream. This motion was more efficient in propulsion. For the numerical results in this study, the velocity fields at H 0.1 c and H 1.2 c correspond well to those of elongated and oblate jellyfish, respectively. At H 0.1 c, although wings can gain a larger aerodynamic force and generate a favorable jet, the efficiency is lower. At H 1.2 c, although the force on the wings is small and no jet is generated, the efficiency is higher. Therefore, the real jellyfish-like flyer is closer to a simplified rowing model without the interference effect between multiple wings.

Conclusion
By solving the incompressible Navier-Stokes equation, we numerically investigated the aerodynamic performance of the flapping wing of the jellyfish-like flyer under some control parameters, including the flapping amplitude, up-flight speed, and the spacing of the wing pair. The wake capture effect and wing-wing interaction were emphasized and discussed. Different from previous 2D models, the 3D wing model with a shape similar to that of the jellyfish-like aircraft wing was used in this study.
With the increase in the flapping amplitude, the lift of the wing first rises rapidly and then becomes approximately stable. As the flapping amplitude increases from 10°to 40°, the vortices generated at the wing edge are significantly enhanced, forming an intenser vortex ring to induce stronger downwash. When the flapping amplitude increases further, the vortex wake, as well as the accompanying downwash, mainly widens in the horizontal direction, but the strength of the vortex and the downwash in the vertical direction changes little, resulting in basically no change in lift.
The wake capture effect changes from reducing the lift at Φ 10°to increasing the lift when Φ is larger than 20°, and with the increase in Φ, the improvement of lift caused by the vortex wake keeps increasing. The vortex wake improves the lift mainly by inducing a larger fluid velocity relative to the lower wing surface and inducing larger edge vortices on the upper wing surface, making the positive pressure and suction force on both wing sides larger.
When the wing has an up-flight speed, the variation in lift with amplitude is similar to the results at hovering. At a specific Φ, the mean lift gradually decreases with the increase in up-flight speed. It is also found that the viscous drag accounts for the reduction in lift at small flapping amplitude, while the change in surface pressure dominates the force variation for large flapping amplitude.
The aerodynamic effect of the wing-wing interaction was also studied. It is found that when the distance between wings is small enough, the wing-wing interaction can improve the lift of the wing, but it requires more power consumption and results in lower efficiency, which is unfavorable to the aircraft. When the wings' spacing is the same as that of the prototype of the jellyfish-like flyer, the aerodynamic interference between multiple wings can be neglected.
The findings in this study may shed some light on further understanding of the aerodynamics of the jellyfish-like flapping wing and provide useful insights into the future design of jellyfishlike aircraft. While the flyer is described as being "jellyfish-like", there are some apparent differences between the kinematics of its wings and that of a jellyfish's bell. During swimming, jellyfish contracts and expands the soft deformable bell very complicatedly, with asymmetries in the contraction and expansion phases. Therefore, it is necessary to consider the effects of flexibility and time asymmetry motion while adopting the 3D wing model, as well as to evaluate the postural stability of the air vehicle.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.