Confined jets in co-flow: effect of the flow rate ratio and lateral position of a return cannula on the flow dynamics

Co-axial tubes have been used to produce a co-flowing confined jet similar to that found in an Extracorporeal Membrane Oxygenation return cannula flow configuration. Particle Image Velocimetry was used to investigate the flow rate ratio between jet and co-flow as well as changes in flow characteristics due to cannula position. The flow was found to be dominated by three main structures: lateral flow entrainment, shear layer induced vortices and backflow along the wall. An increase in cannula flow rate amplified entrainment and recirculation, resulting in a decrease in length required to reach a fully developed flow. Changing cannula position relative the outer cylinder induced a significant reduction in recirculation zone as well as vortex formation on the side to which the cannula was tilted towards, whereas on the other side, the recirculating flow region was enhanced. Proper Orthogonal Decomposition demonstrated that the dominating structure found in the flow is the backflow, composing of several structures having different oscillation frequencies. The significance of the observed and measured flow structures is in enhancing mixing.


Introduction
In several medical applications, e.g. extracorporeal membrane oxygenation (ECMO), a life-saving therapy in the critically ill with cardio-respiratory failure, blood needs to be drained and returned to the patient's circulation requiring accesses via cannulae. It is important to understand the underlying flow characteristics and interaction of cannula and vessel blood flows to improve our knowledge regarding the fluid forces acting on the blood constituents (red and white blood cells, platelets, macromolecules) as well as overall features such as mixing of the two flows. In ECMO, the configuration of the return cannula in the vessel is similar to a confined round jet surrounded by a co-flow.
During the past decades, confined jets in different applications have been studied [13]. In most cases, confined jets for mixing applications have been investigated, focusing on mean velocity and turbulence characteristics. Kandakure et al. [9] used computational fluid dynamics (CFD) to study the effect of the enclosure size, i.e. confinement, over a wide range (6 ≤ D/d ≤ 50). They observed an increase in turbulent dissipation as well as a reduction in turbulent viscosity, entrainment rates and jet spreading angle with decreasing enclosure size. Thus, the jet decayed faster when the surrounding enclosure size became smaller.
Moreover, in comparison to a free jet, the presence of a co-flow affects the spreading and the mixing of a round jet. A co-flowing jet has two asymptotic behaviours, namely strongly advected region and weakly advected region [2,4]. The initial momentum excess per unit density M e0 = U 0 − U a U 0 d 2 p ∕4 determines the jet behaviour. Close to the nozzle tip, the jet is strong as compared to the co-flow and the flow is similar to that of a jet discharging into a still ambient fluid. The jet weakens with downstream distance as its velocity approaches that of the co-flow. Thus, the length scale defined by M 1∕2 e0 ∕U a can be used to estimate the location of the transition between strong and weak jets which occurs where ∼ 1 < x M 1∕2 e0 ∕U a <∼ 100 [4].
The flow rate ratio between the co-flow and the jet is another determinant factor regarding jet development and mixing. Rehab et al. [15] identified two flow regimes depending on the velocity ratio as compared to a critical value. The transition was explained by the turbulent entrainment rate and mass conservation. Such configuration characterised by high co-flow is common in chemically reacting flows, as in burners or in applications where substances are added to the co-flow via an internal jet. However, the case where the jet is faster than the co-flow is less documented in literature. Razinsky and Brighton [14] provided an extensive study on the velocity ratio effect for a faster jet, including turbulence measurements for configurations without recirculation zone. Henzler [8] established that a recirculation zone develops behind the nozzle near the outer cylinder (mixer) walls if the condition D∕d > G +V D ∕V d is valid. Here, D and d are the inner diameters of the outer and inner tubes respectively, V D and V d the corresponding flow rates and G ≈ 1 for equally dense mixed media. Barchilon and Curtet [3] established the dependence of the front boundary and centre positions of the recirculation zone to the flow rate ratio. For an increasing flow ratio by increasing the co-flow flow rate, the recirculation zone moves downstream and the recirculation flow rate q r , defined as the integral of the negative axial velocity over the cross-section, decreases. However, the boundary of the recirculation decay is invariant and remains at a distance of x/D ≈ 3 from the nozzle. A case without recirculation zone was also reported although the condition D∕d > G +V D ∕V d was fulfilled close to its limit of validity, showing that this condition may need to be adjusted. Moreover, Barchilon and Curtet [3] presented yet another parameter based on Euler's theorem to determine the presence of recirculation zones: the Craya-Curtet parameter m or Ct = 1∕ √ m . However, this parameter requires pressure measurements in the flow.
Confined co-flowing jets are often studied for mixers characterised by large diameter ratio (D/d) between the mixer and the nozzle. Razinsky and Brighton [14], Guiraud et al. [7] and Zhdanov et al. [22] presented three of the few studies found in literature about jet characteristics in a confined mixer having a low diameter ratio (D/d ≤ 6). Zhdanov et al. [22] highlighted the existence of a backflow composed by vortices with different length scales and improving mixing properties. However, coherent structures of the flow corresponding to the highlighted frequencies were not studied. Semeraro et al. [16] used Proper Orthogonal Decomposition (POD) and Dynamic Mode Decomposition (DMD) to determine the spatial and temporal structures composing a confined co-flowing jet. The oscillating motion of the jet appeared as two large structures located downstream (x/d > 10) of the nozzle in the first two POD modes. Higher modes showed the coupling between the recirculation zone near the inlet and the shear-layer oscillation.
Data of instantaneous flow structures and corresponding frequencies are lacking for highly confined mixers. Moreover, literature is also lacking regarding the effects due to nozzle lateral position for confined jets. For cannula flows, this is an important aspect that needs to be elucidated as in the case of ECMO, where the cannulae are flexible and the tip position is uncertain within the vessel section. The cannula tip may be centered or close to the vessel walls depending on internal factors (e.g. intravascular fluid status, breath volume, etc.) or external factors such as patient position. This work focuses on the effects of flow rate ratio and nozzle lateral position on the flow structures. The geometrical setup as well as flow rates considered corresponds to that of a (smaller) ECMO return cannula placed in inferior vena cava (IVC). A particular interest is brought to the spatial and temporal characteristics of the recirculation zone. The aim of this work is to investigate the details of cannula flows, but also to provide validation data for CFD, which can handle also the complex rheology of the blood.

Experimental rig
The experimental set-up consisted of two coaxial glass tubes: an outer tube (D = 18.3 mm) and an inner tube, the cannula (inner diameter d = 3.5 mm, outer diameter d o = 5 mm) (Fig. 1). This corresponds to a diameter ratio of D/d = 5.2. At the tip of the cannula there was a small asymmetry on the inner diameter due to fabrication process. The outer tube was 1 m long and connected to a buffer volume to suppress the flow fluctuations induced by the pump. The cannula was 500 mm long and localised at half height of the outer tube. A centrifugal (blood) pump was used to drive the flow in the cannula. The cannula inlet port was connected to a convergent section downstream of a larger tube that was 400 mm long to ensure undisturbed exit conditions. Tube lengths were chosen in order to guarantee a fully-developed flow in each tube prior the flow exiting at the tip of the cannula. Moreover, to investigate the effect of cannula position on flow characteristics, the cannula was fixed to a rotation system at its base, allowing the cannula to be tilted with respect to its position in the outer tube. Due to the dimensions of the tubes, the angle at which the cannula could be tilted was restricted to α < 0.5°. Thus, the tubes could still be considered coaxial although the cannula was tilted. For this reason, the cannula displacement is designed as a lateral shift and its tip position is defined by r c the position of the cannula centre as compared to the outer tube centre. The entire setup was placed into a rectangular glass tank filled with water to reduce optical deformations caused by the outer tube curvature.
The geometry is defined by the axial axis x and the radial axis r. U, u and u′ represent the mean velocity, the instantaneous velocity and the root-mean-square (rms) fluctuations respectively along x-axis. V, v and v′ represent the mean velocity, the instantaneous velocity and the rms fluctuations respectively along r-axis. The axis origin is located on the outer tube centreline at the cannula tip cross-section.
Despite the fact that blood is a non-Newtonian fluid characterized by shear-thinning, the fluid medium used in these experiments was water. Using blood analogies in the experiments would result in flow-dependent viscosity, which implies further uncertainties in the experimental data. By using water, at well-defined temperature, the current experimental data can be used for validation of flow simulations with higher degree of confidence.

Measurement technique
Particle Image Velocimetry (PIV) was used to measure the 2D velocity field in a vertical plane along the cannula axis from the cannula tip to 60 mm downstream, indicated in Fig. 1. A PIV system composed of a high-speed Nd:YAG laser and a Dantec high-speed double frame CCD camera was used. The laser was operating at 200 Hz and delivering pulses of 4 mJ at 532 nm with a duration of 10 ns. The camera (SpeedSense Miro 120 12-bit pixel depth with a sensor 1920 × 1200 pix) was mounted with a Nikkor 105 mm (f/ d max = 1.8) lens (Nikon) associated to a 20 mm extension ring to reach a magnification factor of 31 pix/mm. As it was not possible to place a calibration image with precision inside the tube, the magnification factor was determined through the image of the cannula whose outer diameter is precisely known. The laser sheet had a thickness of 1 mm, measured with a laser profiler camera, and was formed by a convergent and a divergent lens combination. Micrometric translation stages were used to position the laser sheet accurately in the centre of the tube. The velocity range of the flow varied depending on the case investigated. The time between frames was therefore determined so that the particle displacements was limited to ¼ of the interrogation window. The PIV system was synchronized through the Dantec software, acquiring 1000 couples of images at a sampling rate of 200 Hz.
Solid zirconium oxide particles having a diameter around 5 µm, mixed in an aqueous solution with liquid soap to avoid agglomeration, were used as PIV tracers. Based on the flow characteristics, the Stokes number was of the order 10 −2 , indicating the tracer particles correctly follow the flow. Prior to starting the measurement, the loop was run for 5 min to allow for the tracer particle to spread evenly throughout the circuit. The particle concentration was progressively adjusted to reach a particle density on the images high enough for the post-processing of the PIV data to be accurate.
The PIV calculation was performed using an in-house software developed by Bertrand Lecordier [11]. The method uses a multi-pass algorithm with an original iterative deformation of the particle images in order to gradually suppress the bias induced by the velocity gradients. The first iteration was performed with an interrogation window size of 32 × 32 pix progressively reduced to reach the final size of 16 × 16 pix (0.52 × 0.52 mm 2 ). Between each iteration the velocity field was validated by three criteria: Signal to Noise Ratio (SNR), maximum magnitude of the velocity and median filter [20]. To improve the quality of the images, a background correction was applied to rectify the laser sheet inhomogeneity and the background diffusion before the cross-correlation was performed. Finally, each velocity field was validated using the same criteria described above to detect and remove the spurious vectors, which never exceed 0.5% of the total number of vectors. The spatial resolution for the velocity measurement was estimated at 185 µm by using the cut-off wavenumber k c = 2.8/X, where X was the interrogation window size in pixel [5].

Case set-up
In this study, six different cases were studied, see Table 1. The parameters changed were the cannula flow rate and its position. The flow rate in the outer tube was kept constant at 1.3 L/min, corresponding to approximately U co = 0.1 m/s. Thus, the co-flow is laminar for all cases. Two different flow ratios (Q = Q o /Q c ) were investigated: Q = 1 and Q = 0.5. This resulted in a turbulent jet exiting from the cannula represented by a Reynolds number of Re jet = 7880 and Re jet = 15,760, respectively. Due to the diameter of the cannula, 2.6 L/min was the maximum flow rate that could be obtained. For both values of Q, three lateral positions of the cannula were compared. From a numerical simulation point of view, these experimental conditions can be modelled by applying a plug profile for the velocity and a zero pressure gradient condition at the inlet, and zero pressure and zero velocity gradient conditions at the outlet.
In the following, "jet" refers to the cannula jet, "interface" to the interface between the jet and the co-flow and "axis" to the axis defined by the cannula position. Moreover, when the cannula is shifted to the left, the left side of the jet will be referred to as the "narrow side" and the right side as the "wide side" (Cases 3-6).

Uncertainty assessment
One of the objectives of this study was to provide validation data for numerical simulations. Thus, careful attention has been given to the installation of the setup as well as to the quality of the post-processing of the PIV images. However, experimental measurements always include uncertainties due to systematic and random errors.
Geometrical characterisation revealed that the uncertainty associated with the cannula position relative to the outer tube in the plane perpendicular to the laser sheet was estimated to be less than 200 µm. Regarding the lateral position of the cannula in the plane of the laser sheet, measurements by imaging showed an uncertainty of 2 pix, i.e. 65 µm. The uncertainty of the measured inner cannula diameter d was 37 µm. The outer diameter of the cannula d o and the inner diameter of the outer tube D were directly measured using a calliper with an associated uncertainty of 20 µm. Regarding measurement technique related uncertainties, the laser sheet thickness was measured with a beam profiler and had an uncertainty of 8 µm. The uncertainty associated to the position of the laser sheet relative to the centreline of the tubes was estimated to 200 µm. Thus, the uncertainty of the laser sheet profile was considered negligible and the uncertainty associated to the laser sheet thickness was small enough to ensure all particles in the centre plane of the tube to be illuminated.
Regarding acquisition related uncertainty, the delay inter-image due to trigger and laser emission jitters was 76 ns. Applying the improved PIV algorithms (sub-pixel shift) [1], the particle displacement was determined with an accuracy of 0.05 pix. Using the Taylor series, the magnification factor uncertainty was found to be 0.3 pix/ mm corresponding to 1%. Thus, the velocity uncertainty varied from 1.2% (8 pix displacement) to 5.1% (1 pix displacement).
The flow meters had an accuracy of 1%. However, this uncertainty was considered to be included in the uncertainty arising from the reading on the display estimated at 0.05 L/min. Furthermore, a change in room temperature (± 2 °C) was noticed between experiments, leading to small variations of the water properties.
The associated statistical error in the mean value, due to the finite number of samples, was estimated by = ∕ √ N s where σ is the standard deviation and N s the number of independent samples. For the investigated cases, ε ≈ 0.5%.
Thus, the largest uncertainties concerning the geometry were associated to the positioning outside the observation plan, whereas the positioning of the different elements in the observation plan was very accurate and enabled a precise velocity measurement having an uncertainty of less than 5%. The most important source of error in this experiment was the flow rate associated uncertainty of 4%, in particular concerning the velocity in the jet.

Results and discussion
Due to laser reflections on the outer tube and the cannula tip, velocity measurements were carried out starting at a position 1 mm downstream of the cannula tip and at a distance not closer than 0.6 mm from the wall.
To assess the symmetry of the experimental set-up, i.e. study the effect of the side to which the cannula is shifted, measurements were carried out for Q = 1. Figure 2 shows the radial profiles of U and u′ when the cannula is shifted 1.32 mm to the right and to the left relative the centreline of the outer tube (curves for the cannula shifted to the right flipped to ease comparison). As shown in Fig. 2, the radial profiles of U and u′ remain the same for the crosssections where the jet is present (x/d = 0.6 and x/d = 5.7). A difference in the velocity magnitude and rms is observed further downstream (x/d = 17.1). The uncertainties associated with the cannula position and the flow rate measurements are possible reasons of this asymmetry. However, the evolution of the flow downstream of the cannula remains globally identical, independent of the side to which the cannula is shifted. Thus, in the following, results are presented for cannula shifted only to one side. represents the entrainment efficiency of the mass flow from the co-flow in the jet relative to the cannula flow rate. Figure 4 shows the behaviour of this coefficient for the different cases. Close to the cannula (x/d < 1.5) and for both sides of the jet, the entrainment coefficient shows similar trend independent of the case considered. Further downstream, K e is slightly sensitive to the flow rate ratio. However, the cannula position is found to be the most important parameter. Case 5 shows a significant decrease of K e on the narrow side of the jet while K e increases on the wide side for Cases 5 and 6, possibly due to the proximity of the wall limiting the available fluid to be entrained, preventing it to be set into motion by the jet, especially in case of lower jet flow rates. Due to Kelvin-Helmholtz instabilities, vortices are generated in the shear layer developing at the interface (2). In general, the vortices of a free jet are well defined appearing at a given frequency. However, the jet confinement complicates this behaviour, and in particular the recirculation zones found on the sides of the jet (3) are involved in changing jet characteristics. These recirculation zones, indicated by the streamlines, arise from the backflow propagating along the walls, characterised by a pulsatile motion. In this study, D/d = 5.2 and Q = 0.5-1. Considering these values, the condition D∕d > G + Q is largely fulfilled, confirming the backflow presence. The vortex location can be determined by calculating the Γ 2 field. Γ 2 criterion is based on Γ 1 criterion integrating the angular momentum over a domain around a point. Contrary to Γ 1 criterion, Γ 2 criterion takes the convection velocity of the vortices into account [6]. Instantaneous field of the Γ 2 criterion for Case 1, Fig. 5, shows that more vortices are formed on the right side of the jet at that instant even for a jet centred in the outer tube. This is explained by the backflow found on the right side of the jet at that instant. In the same way, when the backflow is present on the left side of the jet, more vortices are formed on the left side. Thus, the presence of a backflow enhances vortex formation. Moreover, when the cannula position is shifted, the vortices cannot develop on the narrow side of the jet due to the proximity of the wall.

Flow structures
For Case 1, the recirculation zone is located downstream of the cannula between x/d = 3 and x/d = 18.3 on the right side of the jet and between x/d = 5.8 and x/d = 18.3 on the left side. This result is similar to that of Zhdanov et al. [22], observing a recirculation zone at 2.5 < x/d < 21 for a similar flow configuration (D/d = 5 and Q = 1.3). The recirculation zones are non-symmetric due to the cannula asymmetric restriction at its tip mentioned in "Experimental rig". However, shifting the cannula, no recirculation zone is formed on the narrow side for Case 3. Increasing cannula flow rate (Case 4) forces the formation of a recirculation zone (R) on the narrow side of the jet although significantly smaller than that of the symmetric configuration. Furthermore, for Cases 5 and 6 (not shown in Fig. 3), the recirculation zone cannot form on the left side, possibly due to the insufficient space between the jet and the wall, and a similar flow as found for Case 3 is observed.
Considering Fig. 3a, the upstream (N) and downstream (P) positions of the recirculation zone boundary are characterised by zero axial velocity, and C denotes the recirculation zone centre (P is located outside the field of view in Fig. 3a). Figure 6 shows axial location for points N, C and P for all cases investigated, deduced from the mean velocity fields. Increasing the cannula flow rate moves the point N upstream and also slightly shifts C upstream. However, the position of the downstream limit P of the recirculation a b Fig. 4 Entrainment coefficient for the cases investigated a on the left side and b on the right side of the jet  Figure 7 shows the radial distributions of the axial velocity, both mean and rms fluctuations for the investigated cases ( Table 1). Note that the data sets for all measured velocity profiles in this study are provided in supplementary materials. Regarding the centred cannula, Fig. 7a, the jet velocity is one order of magnitude higher as compared to the co-flow at the cannula tip. The shear layer between the jet and the co-flow is responsible for the widening of the velocity profiles and its width can be estimated as the distance between the jet core boundary (r = d/2) and the radius where U = U co . Moreover, the widening of the jet velocity profile shows an early entrainment of the co-flow into the jet (x/d = 1.1). Downstream x/d = 5.7, the sign of the velocity changes in the near-wall region of the outer tube. This is due to the recirculation zones found on each side of the jet. From x/d = 20, the axial velocity is positive over the entire cross-section, confirming the recirculation zone ends around this location. Further downstream (x/d = 28.6), a flat profile typical to that expected for a turbulent pipe flow is observed (Re = 3000-4500). The establishment of the turbulent pipe flow can be noted further downstream in the centreline profile of the axial velocity, Fig. 8, where the axial velocity reaches an asymptotic value. The axial velocity decay is similar to that reported in previous studies and clearly faster as compared to a free jet.

Velocity field characteristics
When the cannula position is shifted, Fig. 7e and i, the kinetic energy of the jet maintains its shape until x/d = 2.8 and prevents it from deviating by influence of the coflow, corresponding to the weakly advected region as described by Davidson and Wang [4]. Beyond x/d = 2.8, the jet momentum dissipation through the shear layer that develops between the jet and the co-flow makes the dissymmetry of the flow evident. This is referred to the strongly advected region in which the co-flow controls the jet evolution. When altering the cannula position, the profiles no longer display negative axial velocity on the narrow side of the jet, except in Case 4, in accordance with what was observed in Fig. 3. However, the recirculation zone on the wide side of the jet is larger and has a higher backflow velocity. Figure 9 presents the maximal recirculation flow rate q r,max for the different cases. On the wide side of the jet, q r,max increases with the cannula shift, indicating a stronger backflow. Tilting the cannula, the resulting q r,max is always greater as compared to the centred cannula case even if the contributions from both recirculation zones (Cases 1 and 2) are considered. Furthermore, the maximum of the velocity profile is slightly shifted to the left of the cannula axis showing a deflection of the jet towards the wall. This deviation may be caused by the low-pressure zone created by the water entrained into the jet and not replaced fast enough by the co-flow as observed for a jet along a wall [19]. However, in the case investigated here, this may be further influenced by both sides of the jet being constrained by a wall along with the development of the recirculation zones. The recirculation zone development on the wide side of the jet prevents the jet to be deviated towards the right wall, but instead the recirculation zone acts in order to push the jet towards the left wall of the outer tube.
In case of the lower flow rate ratio (Q = 0.5), Fig. 7c, g and k show that the mean flow maintains similar characteristics as compared to Q = 1. However, the absolute velocity of the jet is increased by a factor of 2. This also leads to an increase in recirculation flow rate, Fig. 9.
In Fig. 7b, minor fluctuations are observed close to the cannula tip except in the shear layer between the jet and the co-flow which is highly unsteady. Fluctuations are found over the entire cross-section downstream of x/d > 2.8, displaying an increasing level of fluctuations further downstream, up to 50% at x/d = 22.9. Thus, maximum fluctuations on the jet axis are found just downstream of the recirculation zone. For a similar geometry (D/d = 5), Zhdanov et al. [22] found that the fluctuations were maximal at 25 nozzle diameters downstream of the nozzle tip for a confined jet with recirculation zone. Moreover, for a confined jet without recirculation zone and D/d = 6, Razinsky and Brighton [14] found the maximal fluctuations at 18 nozzle diameters. Thereafter, the  displayed in (a, b, e, f, i, j) and Q = 0.5 in (c, d, g, h, k, l)  fluctuation level decreases and becomes homogeneous over the cross-section at x/d > 28.6 (the distance where the flow is found to be fully developed, indicated by the mean velocity profile, Fig. 7a). Thus, approximately 30 cannula diameters are needed for the perturbations introduced by the jet to completely dissipate. However, turbulence level in the recirculation zone remains higher than 90%, Fig. 10. Turbulence even tends to be infinite at the front of the recirculation zone (point N). This behaviour is in agreement with Barchilon and Curtet [3]. Thus, in the coflow, the initially low turbulence level increases drastically entering the recirculation zone, being dominant, before being reduced in the fully developed flow. A shift in the cannula position does not have a significant effect on the fluctuations prior x/d = 5.7. However, after this cross-section, an increase in the fluctuations is noticed (Fig. 7f and j). Moreover, although the fluctuation level globally increases downstream of x/d = 5.7, the turbulence level in the recirculation zone decreases, having a minimum of 55% for Case 3 and 45% for Case 5.
A reduction in flow rate ratio Q enhances the fluctuations for x/d ≥ 11.4. Thus, a decrease in flow rate ratio Q leads to a decrease in the length for which a fully developed pipe flow is achieved, where the minimal distance required is approximately of 25 cannula diameters for Q = 0.5.

Spatio-temporal representation of the flow structures
Since the flow is strongly unsteady, it is complex to describe the flow only in terms of averaged structures. To capture the most energy containing structures, POD was applied using the snapshot method with 1000 snapshots [17]. The integral time scale of the flow is estimated in two points: one in the middle of the jet and the other in the recirculation zone. Estimating the integral time scale for each velocity component u and v, T u = 2.7 ms and T v = 3 ms in the jet and T u = 34 ms and T v = 12 ms in the recirculation zone. Thus, focusing on the jet, all the samples are independent. However, if considering only the recirculation zone, this yields approximately 150 independent samples which will restrain the range of frequency that can be studied. Figure 11 shows the energy distribution of the modes focusing on the first 15 modes. POD was calculated after substracting the mean value from the instantaneous measurements. Thus, the contribution of the mode 0 corresponding to the mean velocity field and representing around 90% of the flow energy is removed. So, among the fluctuating structures, it is evident that the first four POD modes contain most energy in all investigated cases. Mode 5 and higher represent each less than 2% of the total fluctuating flow energy. The main difference between the various cases is found in modes 1 and 2. In Case 1, the energy content is only slightly higher for modes 1 and 2 as compared to the higher modes. However, for Case 6, modes 1 and 2 contain 15.3% and 8.9% of the total energy, respectively, whereas mode 3 only contains 2.7% of the total energy. Thus, decreasing Q enhances energy content for modes 1 and 2. The effect of the position of the cannula is less clear. For Q = 0.5, modes 1 and 2 become increasingly dominant with increasing lateral shifting of the cannula, whereas for Q = 1, the maximum energy percentage is reached when the cannula is shifted by only 1.32 mm (Case 3). Thus, modes 1 and 2 are dominant when the cannula is shifted, however the position maximizing this effect is flow rate ratio dependent. The axial and radial components of the velocity of the first four POD modes for Case 1 are presented in Fig. 12. These modes present coherent structures for x/d > 6 corresponding to the location of the recirculation. At each time, the pulsating backflow is a combination of these modes. Thus, the backflow is a superimposition of an anti-symmetric structure (mode 1), a symmetric structure (mode 2) and a helicoidal structure (mode 3). Mode 1 could be interpreted as jet bending in the measurement plane according to the radial velocity. However, the high velocity area of this mode also includes the recirculation zone. Thus, the backflow is slowed down on the side towards the jet is bending and accelerated on the other side. Mode 1 is therefore a combination of the jet bending and the varying backflow velocity. Moreover, it is likely that mode 2 represents a similar type of variation but out of the measurement plane. However, only having access to the velocity measurements in a plane, this remains a hypothesis that cannot be confirmed. Other structures are also present in the flow although not significant for the flow energy. It has to be noted here that the flow being strongly 3D, the POD modes in 2D do not completely describe the flow. Thus, the anti-symmetric structure in 2D actually corresponds to a wobbling motion of the front of the backflow and the symmetric structure probably to a ring of fluid around the jet moving up and down.  In accordance with observations from the velocity profiles, decreasing Q does not modify the most energetic modes. However, the modes are impacted by the position of the cannula. Mode 1 for Cases 3 and 5 is not presented here due to the similarity with Case 1. However, modes 2, 3 and 4 for Cases 3 and 5 presented Fig. 13 are different from those of Case 1. Mode 2 for Cases 3 and 5 does not represent a symmetric backflow on each side of the jet as found for Case 1. In Case 3, mode 2 represents a recirculation zone along the entire tube width while in Case 5, it represents a longitudinal oscillation of the jet accompanied by a flow along the wall in the opposite direction of the jet oscillation. A helicoidal structure as observed for Case 1 is not found for modes 3 and 4 of Case 3, instead these modes represent a recirculating flow behaviour. Mode 3 of Case 5 is similar to Case 1. However, the structures are restricted to the side where the jet is located.
Moreover, this mode forms a pair with mode 4, which was not observed for modes 3 and 4 of Case 1.
Lawson and Davidson [10] and Uzol and Camsi [18] studied a flapping jet in a thin regular cavity and in a fluidic oscillator, respectively. The jet oscillations, accompanied by recirculation zones on each side of the jet, appeared to have a constant Strouhal number for a given geometry for Re jet ranging between 1000 and 80,000. Thus, having a Re jet ranging between 7880 and 15,760, a similar behaviour is expected for the recirculation zone in the present study. Based on the PIV data, a point was selected in the recirculation zone (wide side of the jet for Cases 3-6) to evaluate the power spectral density (PSD) of the backflow in this region. PSD spectrum as a function of Strouhal number  confirms the independence of the Strouhal number to Re jet noticed in previous studies for a given geometry. Thus, the oscillation frequency of the backflow increases proportionally to the cannula flow rate.
However, if the cannula is shifted closer to the wall (Cases 5 and 6), Strouhal number changes as compared to the previous configurations. For Q = 1 (Case 5), Strouhal number slightly increases to 0.0045, whereas, for Q = 0.5 (Case 6), Strouhal number decreases from 0.003 to 0.0018 as compared to the Cases 2 and 4. Thus, when the cannula is close enough to the wall, the backflow frequency changes depending on Q. This is probably due to the absence of the recirculation zone on the narrow side, cancelling its effect on the jet and on the recirculation zone on the wide side. Furthermore, the small value of St indicates that flow oscillations build up and dissipate rapidly due to viscosity preventing large propagation of those oscillations.
A spectral analysis of the POD time coefficients is carried out to compare the frequency pulsation of the most energetic POD modes with the backflow frequency based on the raw velocity fields. The pwelch function in Matlab was applied with overlapping windows on the signal. The length of the window was 90% of the acquisition time length and the overlap 80% of the acquisition time length. Figure 15 shows the frequency spectra of the modes matching the peaks of the spectrum from the streamwise velocity component for Cases 1-4. The POD mode 1 and 2 displays the best match with the main peaks of the u-component spectrum for all cases. For Q = 1, mode 4 also shows a good match with the low frequency peaks (St ≤ 0.001). For Q = 0.5, the match of mode 1 is even better and covers most of the peaks. In Case 2, mode 2 represents the flow structure that seems to pulsate at the secondary frequencies. About Case 5 (not presented here), a very good match of the mode 4 is noticed with the main peak at St = 0.0045. This confirms the different behaviour of the flow for Case 5 noticed in the previous results. Thus, the good agreement between the u-component spectrum and the POD mode frequencies validates the pulsation frequency found for the backflow and shows that the backflow is mostly composed by the symmetric and the anti-symmetric flow structures and pulses at their frequency. A contribution of the other modes at secondary frequencies is also present especially if the cannula flow rate is low.

Conclusion
Velocity measurements have been carried out using PIV on a confined co-flowing jet in a narrow mixer, corresponding to return cannula flows in ECMO. Different • Over the investigated distance range 0.3 ≤ x/d ≤ 17, two different types of interactions between jet and co-flow were observed. Up to x/d ≈ 3, the high kinetic energy of the jet as compared to that of co-flow results in a flow that shows similar development independent of cannula position. Further downstream, the jet is deviated by the asymmetric co-flow. • The flow rate ratios used in this study generates a mixing regime with recirculation zones on each side of the jet for a centred cannula in the mixer. • A lateral shift of the cannula influences the development of the recirculation zones. The recirculation zone does not develop on the narrow side of the jet when the cannula is positioned closer to the wall whereas it expands on the wide side. The recirculation flow rate increases with a shift in cannula position and with an increasing flow rate ratio. • Vortex formation in the shear layer is enhanced by the presence of a recirculation zone.
• For the flow rates investigated here, 25-30 diameters in downstream direction are required to obtain a fully developed flow for Q = 0.5 and Q = 1, respectively.
Spatial decomposition by POD demonstrated that the backflow near the mixer walls is the most energetic structure in the flow. It can be symmetric, anti-symmetric or helicoidal structure with different oscillation frequencies. In the measurement plan, the anti-symmetric mode is always the most dominant especially when the cannula is shifted. However, the 3D structures at the origin of these modes could possibly have a different energy distribution. Flow structures described by the first four modes are not impacted by the flow rate ratio. Moreover, except for the anti-symmetric mode, the main oscillating structures are modified by the cannula position. Frequency analysis showed that the backflow is a structure defined by a Strouhal number between 0.002 and 0.0045. In all cases, the main backflow frequency originates from a combination of the symmetric and the anti-symmetric modes while the secondary frequencies a b c d result from a more complex combination of structures leading to less organized oscillating structures.

Additional material
Files ESM_1-ESM_6 of the additional material provide raw data for the velocity and rms fluctuation profiles, respectively for Cases 1-6. The file ESM_7 contains the summary of the location of the recirculation zone characteristic points (N, C and P).
Acknowledgements Open access funding provided by Royal Institute of Technology. We thank Dr Nils Tillmark, Prof Laszlo Fuchs and the workshop staff at the Fluid Physics Laboratory (KTH Mechanics) for their expertise regarding the design and construction of the experimental rig. We also thank Dr Bertrand Lecordier for sharing his PIV algorithms and its corresponding advices. Funding from Region Stockholm through Grant HMT2018 is greatly acknowledged.
Funding This study was funded by Region Stockholm (Grant Number HMT2018).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.