Classification of vortex patterns of oscillating foils in side-by-side configurations

The unsteady hydrodynamics of two in-phase pitching foils arranged in side-by-side (parallel) configurations is examined for a range of Strouhal number and separation distance. Three distinct vortex patterns are identified in the Strohual number-separation distance phase maps, which include separated wake, merged wake, and transitional-merged wake. Furthermore, a novel model is introduced based on fundamental flow variables including velocity, location, and circulation of dipole structures to quantitatively distinguish vortex patterns in the wake. The physical mechanism of wake merging process is also elucidated. When an oscillating foil experiences the jet deflection phenomenon, secondary structures shed from the primary street traverse in the other direction by making an angle with its parent vortex street. For parallel foils, secondary structures from the vortex street of the lower foil interact with the primary vortex street of the upper foil under certain kinematic conditions. This interaction triggers the wake merging process by influencing circulation of coherent structures in the upper part of the wake. It is unveiled that merging of the wakes leads to enhancements in propulsive efficiency by increasing thrust generation without a significant alteration in power requirements. These are attributed to the formation of a high-momentum jet by the merged vortex street, which possesses significantly larger circulation due to the amalgamation of the vortices, and major alterations in the evolution of leading edge vortices. Thus, flow physics that are thoroughly explored here are crucial in providing novel insights for future development of flow control techniques for efficient designs of bio-inspired underwater propulsors.


Introduction
Bio-mimicking is an innovative way of designing highly efficient robotic platforms for numerous engineering applications. Hence, understanding physical mechanisms employed by natural aquatic species is important to design next-generation autonomous swimming robots. An oscillating foil presents a generic model for the motion used by fish to propel (Anderson et al. 1998). Various studies on such dynamic systems were conducted by researchers in the past few decades, addressing different aspects of the associated wake mechanics and hydrodynamic performance of oscillating foils. Triantafyllou et al. (1993) argued that Strouhal number, defined as St = f A/U ∞ , was one of the governing parameter in fish-like swimming problems. Here, f denotes the oscillation frequency, A is the amplitude of oscillations, and U ∞ shows the free-stream flow velocity. They showed that Strouhal number of highest propulsive efficiency of an oscillating foil overlapped with that of natural swimming for many fish, cetaceans, and marine mammals.
At low Strouhal numbers, oscillating foils produce well-known Bénard-von Kármán (BvK) vortex street. The wake transitions to reverse Bénard-von Kármán vortex street with increasing St (Koochesfahani 1989). A further increase in the value of St triggers the symmetry breaking process in the wake, resulting in formation of deflected (asymmetric) Bénard-von Kármán streets (Jones et al. 1998). This phenomena has been extensively investigated in the literature for single foils (Liang et al. 2011;von Ellenrieder & Pothos 2008;Cleaver et al. 2012). Godoy-Diana et al. (2008) demonstrated that the flow parameters of flapping locomotion in nature coincides with the parameters of oscillating foils that produce deflected wakes. They further conjectured that natural swimmers and fliers could either utilize deflected wakes as their maneuvering strategy or avoid them during forward locomotion. Godoy-Diana et al. (2009) argued that even though three-dimensional (3D) effects influence vortex dynamics of oscillating foils, the wake deflection was a quasi twodimensional (Q2D) phenomenon. The wake deflection occurred when self-advection of the first shed dipole was strong enough to divert the main flow and subsequent dipoles away from the wake centerline. Godoy-Diana et al. (2008) further proposed a model that quantitatively predicted wake deflection, considering offset between dipolar velocity and advection velocity of the dipoles. Although St has a significance on asymmetric characteristics of oscillating foils, amplitude of the oscillation is observed to considerably influence the attributions of deflected wakes. Symmetry breaking is triggered in the wake of oscillating foils at noticeably high oscillation amplitude for a fixed St (Godoy-Diana et al. 2008). Further increase in the amplitude results in the transition from twodimensional wake to three-dimensional wake, which suggests that transition from reverse BvK to deflected BvK is required for the formation of three-dimensional instabilities in the wake (Deng et al. 2015). On the other hand, there are some conditions that inhibits the formation of deflected wakes. For example, interactions between the foil and shed vortices in the wake of an oscillating flexible foil may cancel the constitution of asymmetric vortex street. Furthermore, Calderon et al. (2014) showed that threedimensional effects in the wake of finite span foils hinder the wake deflection, which is observed for the effectively infinite span foils under the same flow conditions. Threedimensionality introduced by the tip vortex, which prevents the vortex coupling, and the symmetric circulation of interconnected vortex loops, which are due to the vortex topology of finite span foil, are two underlying reasons that were provided for cancellation of the deflection.
Efficient propulsion through effective vorticity control implementation has been an outstanding challenge in engineering community for decades. Under inviscid and incompressible flow assumption, trapping a free vortex on the upper surface of a twodimensional wing is theoretically capable of increasing lift generation through the introduction of a low pressure region (Huang & Chow 1982). An adequately stabilized spanwise vortex can enhance the coefficient of lift up to 10 fold, which can be beneficial for design of short takeoff and landing (STOL) aircrafts (Rossow 1978). Saffman & Sheffield (1977) calculated the exact solution for potential flow over a flat plate with a free line vortex positioned on the upper boundary and estimated highly improved lift generation. Leading edge vortices (LEVs) substantially impact, often dominate the wake of simultaneously heaving and pitching foils, where their development could extensively amplify the propulsive performance of foils depending on their formation and shedding. These all are influenced by the foil kinematics. For instance, amalgamation of an LEV and a trailing edge vortex (TEV), which coincides with high efficiency and improved thrust generation, occurs when vortical structures are controlled using various parameters of foil kinematics, such as phase angle between heave and pitch, St, amplitude of heave motion or maximum angle of attack (Anderson 1996;Anderson et al. 1998). Likewise, three distinct vortex patterns are formed behind the foil simulatenosuly heaving and pitching in the wake of a D-section cylinder (Gopalkrishnan 1993;Gopalkrishnan et al. 1994;Shao & Pan 2011). Implementation of active vorticity control by dictating the flow kinematics yields constructive interaction mode, destructive interaction mode and expanding wake modes, which correspond to trough, peak and mixed responses in efficiency.
Natural swimmers and flyers are known to exploit physical mechanisms for their best interest to achieve most efficient way to displace themselves in a fluid medium. Reattachment of LEVs, formed by the flow separation due to dynamic stall, on the upper surface of the wing of hawkmoths or fruit flies during the downstroke of flapping greatly contributes to lift production (Ellington et al. 1996;Birch & Dickinson 2001;Bomphrey et al. 2005). Although balancing the body weight with enhanced lift production plays a crucial role in insect flight, it constitutes insignificant adversity for aquatic animals, owing to the presence of strong bouyant forces. The main concern for aquatic swimmers is to overcome the drag exerted by water, which is three orders of magnitude denser than air. Borazjani & Daghooghi (2013) carried out numerical simulations on selfpropelled virtual swimmers with three different tail geometries inspired from mackerel body. They demonstrated that an attached LEV is formed on the body during locomotion settings that resemble natural swimming conditions for most fish. Evolution of the LEV is remarked to consequentially influence pressure distribution around the tail and the generated force for different tail shapes. In an experimental study, propulsive force of actively swimming bottlenose dolphins was calculated using digital particle image velocimetry (DPIV) measurements of vortex generated by the large amplitude fluke stroke of the dolphin (Fish et al. 2014). Effect of body shape (mackerel body or lamprey body) and swimming kinematics (anguilliform or carangiform) on the hydrodynamics of self-propelled virtual body/caudal fin swimmers was numerically examined by Borazjani & Sotiropoulos (2010) for a range of Reynolds number. It is noted that the form and kinematics of swimmers differently impact the swimming efficiency in viscous, transitional and inertial regimes. Liu et al. (2017) simulated a more complex model, which includes both fin-fin and body-fin interactions, by reconstructing body shape and kinematics of steady swimming crevalle jack using high-speed cameras. They demonstrated that posterior body vortices captured by the caudal fin strengthens LEVs around the fin, which produces most of the swimming thrust.
Fish schooling is defined as a behaviour seen in many fish species that appears as aggregation of a number of individuals and their collective navigation in the flow. Various reasons was propounded by evolutionary biologists and zoologists to explain this habit. These include but are not limited to improved success of finding a mate, effective defence strategy by confusing predators, and better chance of finding prospection. A fundamental question is raised in the minds of engineers interested in bio-inspiration: "Could coordinated swimming enhance the propulsive performance of an individual swimmer?" One of the pioneer works, which addresses this question, was presented by Weihs (1973). They argued, based on a highly idealized, two-dimensional, and inviscid model, that individuals in schooling formation may enjoy hydrodynamic benefits if the spacing and synchronization between swimmers is adequately adjusted. Daghooghi & Borazjani (2015) conducted large-eddy simulations of self-propelled synchronized mackerels in a variety of infinite rectangular schooling patterns. They observed that schooling fish enjoy significant enhancements in swimming speed without more power requirements. They achieve this through exploitation of the channeling effect. Ashraf et al. (2016Ashraf et al. ( , 2017 experimentally approached this problem by examining the swimming of two red nose tetra fish in a shallow-water tunnel with controlled velocity. By tracking the kinematics of fish using stereoscopic video recordings, it is demonstrated that the possibility of fish to locomote in side-by-side configuration and synchronize their tail beat frequencies is strongly correlated with the increasing speed of water. This suggests that in-phase or out-of-phase synchronization of the collectively swimming fish in parallel with a proper spacing provides intensified propulsive performance during demanding flow conditions. Fish schools are often modeled using multiple oscillating rigid hydrofoils arranged in different configurations due to simplicity that it offers. Boschitsch et al. (2014) carried out experiments on the propulsive performance and wake structures of two pitching foils in an in-line configuration for a range of separation distances and phase differences. They observed that both the performance and wake structures of the front foil were affected by the presence of the downstream foil only for considerably small separation distances. They distinguished two different wake modes: branched and coherent. In the coherent mode, a single vortex street is formed behind the follower foil whose time-averaged wake corresponds to a single high-momentum jet. The branched mode, on the other hand, have two angled high-momentum jets in its time-averaged wake. The peaks in the thrust, power and efficiency coincide with the coherent mode wakes while the branched mode wakes are associated with the troughs. Recently, Lagopoulos et al. (2020) focused more on the wake deflection and production of side force by simultaneously heaving and pitching foils in an in-line configuration. They identified three distinct vortex patterns in the wake and showed that wake deflection introduced by the upstream foil could be eliminated due to the presence of the downstream body. Dewey et al. (2014) qualitatively examined vortex patterns of in-phase, mid-phase, and out-of-phase pitching foils in side-by-side configurations for a fixed separation distance and Strouhal number. They then proposed models of wake development for each case. To this end, it was shown that in-phase, out-of-phase, and mid-phase pitching foils produced merging symmetric, diverging symmetric, and asymmetric wakes, respectively. More recently, Gungor & Hemmati (2020) reported on numerical studies of in-phase and out-of-phase pitching foils in parallel arrangements at different Strouhal numbers that maintain a constant gap between them. They showed that the two foils produced quasi-steady symmetric wakes for both phase difference at low St, whereas asymmetricto-symmetric and symmetric-to-asymmetric transitions were observed in the wake at high St for in-phase and out-of-phase oscillations, respectively. A similar symmetry breaking phenomenon in the wake of foils, performing out-of-phase oscillations in a parallel configuration, was observed by Bao et al. (2017) and Zhang et al. (2018). However, they demonstrated the asymmetric wake at a single time instant without examining the transient formation process of their asymmetry.
Although there are a few studies that demonstrate the development of vortex structures behind parallel foils, there are none that provide quantitative explanations for the vortex interactions in the wake. Furthermore, studies focused on explaining vortex patterns in the wake mostly overlooked unsteady interactions and their impact on propulsive performance, which are expected at high St. In this study, we examine merging of the vortex streets in the wake of in-phase pitching foils in side-by-side arrangement at a range of St and separation distance inspired from fish schools. The Reynolds number of the flow (Re = U ∞ c/ν, where, c is the chord length of the foil and ν is dynamic viscosity of the fluid) is fixed at Re = 4000 considering wake patterns of oscillating foils reach a plateau after Re 1000 (Das et al. 2016) and their propulsive performance exhibits negligible alteration after Re 4000 (Senturk & Smits 2019). Therefore, this paper aims to illuminate three novel points that are currently missing in literature: (i) quantification and classification of vortex patterns behind two parallel foils undergoing in-phase pitching oscillations, (ii) physical mechanisms, governing the wake merging phenomenon, and (iii) influence of the merger on propulsive performance of the system. For this purpose, this paper is structured as follows. A description of dynamic system composed of two pitching foils and our numerical setup is provided in section 2. Section 3 includes the results on the wakes of parallel foils and discussions concerning the vortex patterns and wake merging phenomena, which is followed by main conclusions in section 4.

Methodology
The flow around two oscillating, rigid teardrop foils in a side-by-side configuration is numerically simulated using OpenFOAM. For this purpose, Navier-Stokes equations are directly solved using pimpleDyMFOAM solver, which is an incompressible transient flow solver for systems requiring dynamic grids. The solver utilizes PIMPLE algorithm, which is a hybrid of PISO (Pressure-Implicit with Splitting Operators) and SIMPLE (Semi Implicit Method for Pressure Linked Equations). The time-step size is adequately selected to limit Courant number of the flow below 0.8 throughout the domain. It is achieved by using over 3500 time-steps per oscillation cycles. The divergence terms of Navier-Stokes equations are discretized using upwind-biased, second-order accurate "Linear Upwind" technique. Second-order, implicit backward time method is employed for temporal terms. The convergence criterion, which is the residual of velocity components and pressure in the momentum equations, is set to 10 −5 .
Both foils, Foil 1 (lower foil) and Foil 2 (upper foil), have chord-lengths of c and semicircular leading edges with radii of 0.05c. They perform pure pitching motion, which is mathematically defined as: Here, θ 0 is the pitching amplitude, t is time, and φ is the phase difference between the two foils. In-phase pitching condition are imposed by setting φ = 0, and θ 0 is fixed at 8 • . A schematic representation of the in-phase pitching is provided in Figure 1a. The grid is morphed by the solver in each timestep in order to ensure the pitching motion, while maintaining its quality. The separation distance between the foils (d) is varied from 0.5c to 2c with increments of 0.5c and is non-dimensionalized by c, i.e., 0.5 < d * = d/c < 2. Moreover, St ranges from 0.15 to 0.5 for each value of d * . This parameter space (d * extend) includes the range of separation distance between two red nose tetra fish, synchronously swimming in side-by-side arrangement, during the fish tank experiments by Ashraf et al. (2016Ashraf et al. ( , 2017. The St space covers the formation of BvK, reverse BvK, and reverse BvK regimes in the wake of single oscillating foils (Godoy-Diana et al. 2008) and natural swimming St of various fish species (Triantafyllou et al. 1993).
A computational domain, similar to the one reported in our previous work ) is employed in this study, which follows the experiments of Dewey et al. (2014). It extends 30c in the streamwise (x−) direction and 16c in the cross-flow (y−) direction. Also, the leading edge of the foils are placed 8c away from the inlet. Neumann condition for both pressure and velocity are applied at the outlet boundary, while a uniform velocity (u = U ∞ , v = 0, w = 0) is prescribed to the inlet boundary. Boundary conditions for the upper and lower walls, and foil surfaces are selected to be slip and no-slip, respectively.
A non-homogeneous spatial grid, consisting of 7.87 × 10 5 hexahedral elements, is generated to simulate the flow. The grid is most refined around the foils with 600 nodes on the surface of each foil, which is consistent with the numerical setup of Senturk & Smits (2019). The grid size expands towards the boundaries without exceeding the expansion ratio of 1.03 in the entire computational domain. Sensitivity analyses for grid, time-step, and domain sizes as well as a validation study of our computational methodology are provided in . More details of the presently utilized grid around the foils are presented in Figure 2.
Two-dimensional versus three-dimensional simulations are an important numerical complexity that can have implications on wake dynamics at high Re flow conditions. To this effect, we carried out three-dimensional sensitivity studies to confirm that underlying physics of coherent structures in the flow, including wake deflection, wake merging, and vortex interactions, follow a two-dimensional or Q2D mechanism (Godoy-Diana et al. 2008, 2009Dewey et al. 2014;Shoele & Zhu 2015;Lagopoulos et al. 2020). Deng et al. (2016) notes that 2D-to-3D transition in the wake of pure in-phase pitching foils occur at considerably high St, which excludes the parameter space employed here. Contour plots in Figure 3 compare coherent structures, and their interactions, along the center xy−plane of the wake at Re = 4000, St = 0.3, and d * = 1 with those from two-dimensional simulations. These results confirm that two-and three-dimensional simulations render very similar results in terms of coherent structures, wake dynamics and interactions related to merging. Moreover, it has been previously established that there are no significant variations observed between two-and three-dimensional cases in studying propulsive performance of pure pitching foils at moderate St, e.g., thrust, efficiency and power (Zurman-Nasution et al. 2020). This range comprises St of the flow examined in the current study. Note that the impact of three-dimensionality on wake structures (Deng et al. 2016) and performance (Zurman-Nasution et al. 2020) becomes remarkable at relatively lower St for pure heaving foils.
The cycle-averaged coefficients of thrust ( C T ) and power ( C P ) together with Froude efficiency (η) are calculated to discuss propulsive performance of the foils. These parameters are defined as: (2.5) Here, F x is the streamwise force applied by the foil to the fluid, M z is the moment in z direction applied to the foil, ρ is fluid density, s is span of the foil. F x and M z are averaged within each oscillation cycles over at least 3500 timesteps.

Results & Discussion
We begin our analysis with examining vortex dynamics and wake interactions of parallel pitching foils. In a previous study (Gungor & Hemmati 2020), we examined transient wake developments of foils, performing in-phase and out-of-phase pitching in side-by-side configurations for St = 0.25 − 0.5 and d * = 1 at Re = 4000. It was demonstrated that wake structures showed quasi-steady characteristics and were in perfect agreement with the findings of Dewey et al. (2014), i.e., merging symmetric wake for in-phase pitching and diverging symmetric wake for out-of-phase pitching foils. However, wake structures and propulsive performance of both in-phase and out-of-phase pitching foils were observed to be highly transient. The wake of in-phase pitching foils initially consisted of two deflected vortex streets parallel to each other. These streets merged after some time and formed a symmetric wake. The merging process coincides with the enhancement in time-averaged thrust and efficiency of the foils. The opposite phenomena was observed in the wake of out-of-phase pitching foils. The foils initially produced diverging symmetric wakes whose symmetry was broken after several oscillation cycles. Here, we expand the parameter space to introduce and build classifications of the vortex patterns, which elucidate active flow control techniques possibly employed by natural swimmers, to gain a desired hydrodynamic performance. Furthermore, we also present a quantitative explanation for underlying physical mechanisms of the wake merging phenomenon.

Classification of Vortex Patterns
We identify three distinct vortex patterns in the wake of in-phase parallel pitching foils (side-by-side configuration) for the given parameter space, i.e., 0.5 d * 2 and 0.15 St 0.5. Merged−separated characteristics of the wakes were taken into consideration when classifying the wake in Figure 4. Here, a merged wake corresponds to the vortex topology that involves the vortex streets shed by upper and lower foils merging in mid-wake and forming a single street, which constitutes a new flow configuration. In separated wakes, on the other hand, upper and lower vortex streets remain isolated without interacting with each other. Both merged and separated wakes demonstrate quasi-steady characteristics as the vortex patterns are formed within several pitching cycles and remain stable during the next oscillation cycles without altering significantly (compare t 1 = 14P and t 2 = 20P of Figure 4a and 4b, where P is period of the pitching cycle). Conversely, transitional-merged wakes undergo distinct separation and merging stages, transitioning from the former to the latter configuration. As explained earlier, oscillating foils produce deflected reverse BvK vortex streets at considerably high St. In the wake of parallel foils, interaction between vortex streets shed by each foil results in the constitution of the symmetric wake, in which upper and lower wakes amalgamate around the centerline. (see Figure 4c). The pitching cycle, in which the merging takes place, greatly depends on d * and St (see Table 1). For the sake of comparison, the merging process of the wakes occurs around 22 nd cycle for d * = 1 and St = 0.5, whereas more than 75 cycles are needed for this phenomenon to occur for d * = 2 and St = 0.5. These vortex patterns were gathered in a St − d * phase diagram in order to provide a thorough classification of the wakes of parallel foils in Figure 5. In this diagram, separated and merged wakes are observed in upper (d * 1.5) and lower (d * 1) regions of the diagram, respectively. On the other hand, transitional-merged wakes fall into the high St region. It is important to note that this type of wakes are formed only at sufficiently high St that facilitate the formation of deflected wakes. The relation between the deflection phenomena and wake merging will be further explained in the next part of this section.
Thereafter, we present quantification of important characteristics of the wakes and propose a model that disfperftinguishes emerging vortex patterns. At this point, it is important to recall the dipole model by Godoy-Diana et al. (2009), which presented a quantitative threshold for the wake deflection behind a single oscillating foil. This model also holds true for the wake of undulating foils (Khalid et al. 2020). The wake of oscillating foils that consist of a reverse BvK vortex street is dominated by shedding of a counterclockwise (positive sign) and a clockwise (negative sign) vortex per oscillation cycle that are located slightly above and below the centerline, respectively (Koochesfahani 1989). The structure that consists of these two vortices is called dipole. Circulations of the vortices in the dipole induce velocity normal to the line that connects the vortex centers as described by the two-dimensional Bio-Savart rule (Naguib et al. 2011). When the selfadvection velocity of the dipole is strong enough, it diverts the dipole from the centerline, which is followed by the consecutive dipoles. Therefore, the model was based on the offset between advection velocity of the propulsive wake, i.e., U phase , and self-induced translation velocity of the dipole, i.e., U dipole . They can be mathematically defined as follows: where X i is the x-coordinate of the spatial position of a vortex core, Γ denotes the circulation of counter-rotating vortices, and ξ represents the distance between the centers Separated wake circle Merged wake square Transitional-Merged wake diamond Figure 5: Classification of the wake patterns of in-phase pitching foils in side-by-side configuration at a range of separation distance and Strouhal number for Re = 4000. Dashed lines correspond to the boundary that distinguished merged and separated wakes.
of the vortices (see Figure 6a). Circulation (Γ ) is computed either from a line integral of the velocity field or from a surface integral of voricity over the area bounded by a closed curve. Godoy-Diana et al. (2009) used a rectangular frame, whose size was determined by Gaussian fit, to extract the boundary of each vortex towards calculating Γ . However, this method has a downside of potential numerical errors due to the possibility that rectangular frames may include counter-rotating vortices, particularly in the case of structures traversing in close proximity of one another. Therefore, we use an arbitrary closed curve to accurately capture each vortex proposed earlier by Khalid et al. (2020), which eliminates this error in computing Γ . The boundary of the curve is defined such that it only encompasses the region with magnitude of vorticty (ω) greater than 5% of its value in the flow field. Then, we determine the circulation of each vortex by calculating the line integral of the velocity field around the curve using the following definition: where α is the angle between U phase and U dipole as presented in Figure 6a.
Here, we present a model that distinguishes different classes of vortex patterns using U * p . Although the model of Godoy-Diana et al. (2009) successfully predicts whether the wake is deflected behind an isolated oscillating foil, but it cannot properly identify the nature of vortex patterns, i.e., merged or separated, for multiple parallel foils that exhibit wakes in close proximity of one another. In order to construct an effective mathematical model, our current work focuses on differentiating merged and separated wakes and supplying information about the direction of their deflections. To illustrate it further, transitional-merged wake at St = 0.5 and d * = 2 exhibits deflection during each stage of wake development. During the separated stage at t 1 = 30P , both top and bottom wakes are deflected downwards, whereas the wake fully transitions to that of a merged configuration at t 2 = 90P . In the latter stage, upwards deflected bottom vortex street  Figure 7a), which produce deflected vortex streets, and at lower St and 1.5 d * 2, which exhibit horizontal vortex streets (e.g. Figure 4b), are treated disparately, while they are all classified as separated wake given the wake deflection is present only in a single case. Thus we introduce the term sin α to the formulation to take the direction of deflection into account, because sin α yields positive values for upwards wakes and negative values for downwards and non-deflected wakes. Thus, the effective phase velocity of the coupled vortex system or U * p,sys can now be defined as follows: Here, U * p,upper and U * p,lower are the effective phase velocities computed for upper and lower vortex streets. This mathematical relation produces negative U * p,sys values for merging wakes due to the opposite deflection directions. In order to examine this model using our results, vortex dipoles are tracked as they move downstream the wake and U * p,sys is calculated with respect to the distance dipoles travel after shedding and full separation from the foils. The displacement of these coherent structures is given by r = (X 1 − X 0 ) 2 + (Y 1 − Y 0 ) 2 , where X 1 and Y 1 describe the instantaneous location of a vortex core. X 0 and Y 0 provide location of the vortex core just after its full detachment from the foils. Note that geometric mean of the respective quantities for the counterrotating vortices is used as the location of the dipole. It is evident from Figure 6b that the proposed model works perfectly for the given parameter space.
In this model, separated and merged stage of transitional-merged wakes are treated individually and marked with different colors, since these stages are contradictory to one another in terms of their vortex configuration. It is important to note that merged wake cases are tracked for relatively short radial displacement, i.e., r/c 3. This is because their upper and lower vortex streets merge at mid-wake, which inhibits further tracking. However, dipoles of the separated wakes are traceable until circulation of the vortices shrink to negligible values due to the viscous diffusion around r/c = 5. To clarify the working mechanism of the model, U * p of merged wakes (see Figure 4a or Figure 4c at t 2 = 90) yield positive values for both bottom and top vortex streets as they are deflected upwards and downwards, respectively. Furthermore, sin α for top and bottom wakes switch signs (sin α < 0 for top and sin α > 0 for bottom), which results in U * p,sys < 0. On the other hand, horizontal vortex streets in the separated wakes (see Figure 4b) have U * p < 0 and sin α < 0, thus leading to U * p,sys > 0. Finally, the separated wake, whose vortex streets are deflected (see Figure 7a), or transitional-merged wake at separated stage (see Figure 4c at t 1 = 30P ) yield U * p > 0 and sin α < 0, which translates to U * p,sys > 0.

Mechanism of Wake Merging
Now that we have established a mathematical model to quantitatively identify and characterize different wake patterns behind parallel foils, we focus our attention to identify and explain the mechanism of wake merging. To this effect, we analyze the wake merging phenomenon by associating it with the production and dynamics of secondary vortex structures. When an oscillating single foil produces a deflected wake, secondary structures appear from the primary vortex street to move away from the direction of deflection . Such secondary structures were also observed in the experiments of Godoy-Diana et al. (2008) and Jones et al. (1998) and numerical simulations of Liang et al. (2011). But no further analyses were conducted for this important feature of the wake dynamics. These structures are considerably weaker in their relative strength compared to those in the main street, which is why they have not received adequate attention in literature. We hypothesize that secondary structures play a key role in the merger of upper and lower vortex streets behind parallel oscillating foils. Figure 7a shows separated wakes of the foils for the cases with St = 0.4 and d * = 2 at t = 60P . Even though secondary structures are present in the wake, structures from the lower wake are convected in the downstream direction before reaching the upper wake. At St = 0.5 and d * = 2, in which wakes are merged, however, these secondary structures from the lower foil deflect upward to interact with primary street of the upper foil, traversing downward (see Figure 7b). These observations hint that this interaction triggers the wake merging process, because constructive or destructive interference of secondary vortices with the bigger coherent structures change their strengths in terms of circulation (Zhu et al. 2002;Akhtar et al. 2007;Khalid et al. 2021b,a). Furthermore, onset of the complete wake merging appears located very close to the point of interaction between secondary structures and primary street. For instance, secondary structures shed by the lower foil at St = 0.5 and d * = 1.5 catches the upper main vortex street at x/c ≈ 3, which coincides with the location for the merging of wakes (please see Figure 4c). Likewise, both the interactions between secondary structures and the upper wake as well as the merging of upper and lower wakes occur at x/c ≈ 4.65 for St = 0.5 and d * = 2 (please see the supplementary video and Table 1). The alignment of spatial merging location and that of vortex interactions strengthens our argument on the role of these smaller structures in initiating and facilitating the wake merging process. We proceed with providing another quantitative explanation to the impact of secondary structures on the overall wake dynamics. In this manner, locations and circulation of vortex dipoles are traced in the wake to provide evidence for the impact of secondary structures. Figure 8 exhibits the change of Γ for negative vortices associated with separated (St = 0.4 − d * = 2) and transitional-merged (St = 0.5 − d * = 2) wakes, as dipoles move downstream the wake. For the transitional-merged wake (Figure 8b), circulation of the negative vortices of upper and lower wakes overlap in the near wake region. However, proximity in this sense is broken in favor of the upper wake, after which there is constructive interference of secondary structures with negative vorticity at r/c ≈ 3.5. This observation remains valid at different time instants. Non-dimensional circulation at t 1 = 44P and t 2 = 55P is presented in Figure 8b. On the contrary, there exists no significant difference in circulation of the upper and lower wakes for the separated wake (Figure 8a), because the secondary structures have no influence on the upper wake. Similarly, wakes at t 1 = 40P and t 2 = 50P show that this trend is independent of the wake evolution and time. Furthermore, circulation of negative vortices of the upper and lower wakes is computed for other separated wakes and transitional merged wakes to support the set explanation for this mechanism. It is evident from Figure  Transitional-Merged Wake Figure 9: Magnitude of non-dimensional circulation of negative vorticity of upper and lower vortex streets for transitional-merged wakes before the merger and separated wakes. Table 1: Streamwise location (x/c) and time instant (t/P ) in which the wake merging occurs as well as the percent improvement in the cycle-averaged coefficient of thrust (∆ C T ) for separated and transitional-merged wake cases at St = 0.4 and St = 0.5. 9 that our illustration stays valid for all transitional-merged wake cases investigated in this study.

Effect of Wake Merging on the Propulsive Performance of the System
We now focus on the relationship between propulsive performance of parallel foils (sideby-side configuration) and the wake merging phenomena by assessing the cycle-averaged performance metrics. Figure 10 shows temporal variations of the system-averaged (Foil 1 and Foil 2) coefficients of thrust and power, as well as efficiency at a range of St and d * . The performance parameters of this dynamical system is examined for 65 oscillation cycles. The exception is the case of St = 0.5 and d * = 2, in which the wake merging process occurs later at the 78 th cycle (see Table 1). Thus, we examined this particular case using 90 cycles instead. Figure 10 covers all transitional-merged wakes observed in the current study as well as a separated wake, i.e. St = 0.4 and d * = 2. In this section, we aim at establishing the impact of unsteady alterations in wake dynamics on the propulsive performance of the system. Thus, parameters for cases with lower St, i.e., St 0.3, are not presented here as they display no significant wake transitions. Previously, it was demonstrated by Gungor & Hemmati (2020) that parallel foils generated highly quasi-steady performance and wake characteristics for lower St at d * = 1. With this background, It is further noticed in the present study that the same attributions persist for other separation distances examined here. Power requirements of the system marginally vary in time, which translates into resembling trends for the cycle-averaged coefficient of thrust and efficiency. The percent improvement in thrust calculated between 5 th and 65 th cycles (5 th and 90 th for St = 0.5 and d * = 2) is given in Table 1 together with the location and time instant of the wake merging. It is evident from Figure 10a that generated thrust by the transitional-merged wakes improves with time and reaches a steady-steady after the wake is fully merged. This indicates that the wake merging process could be a contributing factor for thrust enhancement. Furthermore, it is important to note that the separated wake presented here (St = 0.4 and d * = 2) experiences trivial alterations in propulsive performance parameters, which further supports our argument. Thus, we affirm that merging of these wakes improves propulsive thrust of the system by increasing thrust generation through amplification of the circulation associated with amalgamated vortices around the mid-wake. This results in the formation of high momentum jet on the centerline. There are two consequential inferences from Figure 10 and Table 1, which strengthens our argument. First, the time instance of peaks in thrust variation lags the instant of wake merging process. For example, thrust generation for the case of St = 0.5 and d * = 1 has its maximum at t/P = 29, whereas its wake merging occurs at t/P = 22. Second, thrust enhancement in the system decreases as the streamwise location of the onset of wake merging moves downtream. This is due to the reduction in circulation of the dipoles as they travel downstream the wake (see Figures 8 and 9). When upper and lower vortex streets merge in closer proximity to the foils, amalgamated vortices yield more circulation. This translates to an increased momentum jet that is induced by the vortex street. This is consistent with insignificant improvements in thrust for the case of St = 0.4 and d * = 1.5, whose wake merging occurs farther in the wake, i.e., increased distance from the foils. We now proceed with implementing the finite-core vortex array model to our wake data, following Naguib et al. (2011). This model suggests that velocity profiles in the wake can be determined by superimposing finite amount of vortex cores onto a uniform flow. Streamwise and transverse velocity profiles that induced by superposition of N vortices can be determined using: where r i is the radial distance from i th vortex center to the point of calculation, x ci and y ci are the streamwise and transverse location of the center of the i t h vortex, respectively. The model requires the number of vortices greater than or equal to 10 in order to converge (Naguib et al. 2011). To this end, locations and circulations of vortices within the range of 3 x/c 7 are measured for St = 0.5 and d * = 1 at t 1 = 13P (separated stage) and at t 2 = 50P (merged stage). The St and d * of the flow are selected, considering it yields the highest percent improvement in thrust production (see Figure 10 and Table 1). Moreover, the range of execution is determined considering that the wake merging occurs following x/c = 2.5 and circulation of the dipoles diminish after x/c 7. Mean streamwise velocity profiles calculated using the vortex array model is plotted for x/c = 4 and x/c = 6 in Figure 11. High velocity excess is observed around the centerline (y/c = 0) for the merged wake (t 1 = 50P ), whereas two distinct peaks associated with the jets created by the upper and lower vortex streets are visible in the separated wake (t 1 = 13P ). Excess momentum at the outlet profiles (x/c = 4 and x/c = 6) created by the single velocity peak is greater than combination of the two peaks. It clearly shows impact of the wake merging on the formation of high-momentum jets, which translates to improvements in thrust generation. Furthermore, we compare velocity profiles obtained from the vortex array model with contours of mean horizontal velocity from Figure 12. The model accurately captures the general trend and locations of the velocity peaks. However, it underestimates the magnitude of peaks in velocity profiles. This short-coming may be due to the sampling space given we only considered the region after wake merging in our calculations for the model. This focus on a particular region was driven by our emphasis on the influence of wake merging.
The formation, growth and interaction of LEVs can be an important mechanism that impact propulsive performance of oscillating foils (Anderson et al. 1998;Pan et al. 2012), while also impacting the wake development (Hemmati et al. 2019). To this effect, we now look at how alterations in the formation and growth of LEVs around the foils influence their propulsive performance during separated-to-merged wake transition. LEVs are formed when the angle of attack is high enough that a separation bubble is formed on the foils. Their presence and evolvement on the surface of a fin or wing is responsible for a large part of thrust and lift generation in aquatic locomotion (Borazjani & Daghooghi 2013;Bottom Ii et al. 2016;Liu et al. 2017;Xiong & Liu 2019) and insect flight (Ellington et al. 1996;Birch & Dickinson 2001;Bomphrey et al. 2005), respectively. Unsteady thrust variations of the foils in the transitional-merged wake (St = 0.5 and d * = 1) throughout the separated and merged time ranges is presented in Figure 13, which clearly demonstrates that higher peaks and lower troughs are achieved after the wake merging. Contours of vorticity focusing on surfaces of the foils at time instants that correspond to the highest (θ = 8 • ) and the lowest (θ = 0 • ) angles of attack is shown in Figure  14. This enables comparing this process to the evolution of LEVs around the foils. Note that these instants roughly overlap with the times at which the foils yield their highest and lowest thrust generation, as marked in Figure 13. Negative (clockwise rotating) and positive (counter clockwise rotating) LEV are formed on the upper and lower surface of the foils at the separated stage of the wake evolution, respectively. On the other hand, constituting positive LEVs on the lower surface of Foil 1 and negative LEVs on the upper surface of Foil 2 is significantly suppressed when the wake is fully merged. Furthermore, it is evident from contour plots in Figure 14 that stronger LEVs are formed after the wake merging, which hints at a potential factor for thrust enhancement. These observations are valid for the other transitional-merged wake cases with recognizable thrust enhancement as well, although they are not shown here for brevity. Note that profiles of unsteady thrust have two peaks and two troughs per oscillation cycle. The impact of LEVs on thrust generation can be explained through the low pressure region (suction) formed by vortices. LEVs attached close to the anterior part of the foil favorably affect thrust by dropping the pressure in this region, whereas their influence is adverse if located around the posterior part of the foil. For example, Foil 1 at t 7 = 50.25P have enhanced LEV around the front edge comparing to Foil 1 at t 3 = 13.25P . However, the distribution of LEVs on rear surfaces of the foils are matching, which translates to thrust enhancements for Foil 1. Likewise, thrust of Foil 2 at t 5 = 49.75P is considerably larger than that of Foil 1. This is due to stronger LEVs formed close to the forehead of Foil 2, while an LEV with large negative vorticity is present on the rear part of the Foil 1.
It has been previously demonstrated by Gungor & Hemmati (2020) that merging of vortex streets is accompanied by the restoration of wake symmetry for parallel foils. Both performance and wake characteristics exhibit symmetric behavior with a delay of half pitching period. This lag between the foils is due to the formation process of vortex dipoles. Although TEV 1 and TEV 2 are shed at the same time (e.g., t 5 = 49.75P ), TEV 2 establishes a dipole with TEV 0 that has been shed half a period prior to these TEVs, whereas the coupling of TEV 1 and TEV 3 are delayed by half a cycle. Therefore, distribution of LEVs around Foil 1 and Foil 2 at t 5 = 49.75 in Figures 14e is mirror image symmetric with switched signs of vorticity with that around Foil 2 and Foil 1 at t 7 = 50.25P in Figures 14g, respectively. Besides, shifting C T values for Foil 1 by a half period in either direction results in a perfect overlap with those for Foil 2 during the merged stage of the wake evolution, or vice versa (see Figure 13). Contrarily, there is neither a coherent similarity of the arrangement of LEVs around the foils between / / 1 = 12.75 5 = 49.75 2 = 13 6 = 50 3 = 13.25 7 = 50.25 4 = 13.5 8 = 50.5 Figure 13: Variations in unsteady thrust coefficient of Foil 1 and Foil 2 for St = 0.5 and d * = 1 (transitional-merged wake) during separated stage (12 t/P 14) illustrated in black and merged stage (49 t/P 51) of the wake evolution illustrated in red.

Conclusions
Numerical simulations on the flow around two in-phase pitching foils in side-by-side configuration are examined at a range of Strouhal number, 0.15 < St < 0.5 and separation distance, 0.5 < d * < 2, at Reynolds number of 4000. First, we classified the vortex patterns in the wake. Separated and merged wakes, which exhibit quasisteady performance and wake characteristics, were observed at lower Strouhal numbers. Small spacing between the foils yielded the constitution of merged wakes, while separated wakes are seen at higher separation distances. On the other hand, transitional-merged wakes, which are often observed at high Strouhal numbers, exhibited wake evolution in time. Two distinct and deflected vortex streets shed by each foil were observed at early stages of the oscillations. Upper and lower vortex streets approached each other in time, which eventually resulted in merging of the wakes. A novel mathematical model was proposed, which quantitatively established the threshold for the two set vortex patterns. This model utilized the locations and circulations of individual vortices in a dipole. It was further tested using the current parameter space and performed perfectly in determining if the wake is separated or merged. Then, we proceeded with evaluating and explaining the physical mechanism associated with the wake merging phenomenon. This analysis revealed a novel process in which secondary structures in the wake were responsible in part for the wake merging. The wake merging occurred when secondary structures from the lower vortex street were strong enough to form a constructive interaction with main vortex street of the upper wake. This interaction triggered the merging of wakes by increasing the circulation of negative vortex in the upper vortex street. In turn, this impacted the resultant induced velocity (flow) by the two vortex streets, which now do not match, leading to further deflection of wakes and their subsequent merger. Finally, it was observed that merging of the wakes enhances propulsive performance of the foils by combining circulations of amalgamated vortices. This process induces high-momentum jet around the centerline. Evolution of leading-edge vortices played a major role in the performance enhancement. Alterations in the distribution of leading-edge structures and the amplification in their strength, which occurred after the wake merging, was a contributing factor for the improvements in thrust generation.