Large-scale volumetric flow visualization of the unsteady wake of a flapping-wing micro air vehicle

The objective of this experimental investigation is the volumetric visualization of the near wake topology of the vortex structures generated by a flapping-wing micro air vehicle. To achieve the required visualization domain (which in the present experiments amounts to a size of 60,000 cm3), use is made of robotic particle image velocimetry, which implements coaxial illumination and imaging in combination with the use of helium-filled soap bubbles as tracer particles. Particle trajectories are determined via Lagrangian particle tracking and information of different phases throughout the flapping cycle is obtained by means of a phase-averaging procedure applied to the particle tracks. Experiments have been performed at different settings (flow speed, flapping frequency, and body angle) that are representative of actual flight conditions, and the effect of reduced frequency on the wake topology is investigated. Furthermore, experiments have been carried out in both tethered and free-flight conditions, allowing an unprecedented comparison between the aerodynamics of the two conditions.


Introduction
An increasing interest in the development of micro air vehicles (MAVs) has been seen in recent years due to the wide range of missions, where they can be employed, from search and rescue activities to security and surveillance applications. These MAVs typically fly at very low speeds or even in hovering conditions, but high manoeuvrability and flight efficiency are desirable for wind gusts rejection and increased flight endurance, respectively. Different classes of MAVs can be distinguished: fixed wing, rotatory wing, and flapping wing (Mueller 2001;Mueller and DeLaurier 2003). Well-established aircraft design concepts apply to fixed and rotatory wing configurations, but these suffer from a degradation in performance when scaled down in size (Jones and Babinsky 2010). On the other hand, flapping-wing propulsion of natural fliers has shown extraordinary flight capabilities at low Reynolds numbers (Cheng et al. 2016;Muijres et al. 2014;Karasek et al. 2018), thus serving as an inspiration to the development of flapping-wing MAV (FWMAV) designs (Shyy et al. 1999(Shyy et al. , 2010Mueller 2001). Multiple bio-inspired FWMAVs have been developed worldwide (e.g., Nakata et al. 2011;Keennon et al. 2012;Ren et al. 2013;Mazaheri and Ebrahimi 2011), among them also the DelFly of TU Delft (de Croon et al. 2012. The DelFly will be used in this study as subject to investigate the aerodynamics of flapping wings and the vehicle itself will be described in further detail in Sects. 1.2 and 2.1. Flapping-wing propulsion constitutes an area that holds interest in the fields of both biological and aerospace sciences and this shared interest has stimulated progress in flapping-wing aerodynamic research. Quantitative visualizations techniques, particle image velocimetry (PIV) in particular, have led to a better understanding of flapping flight, supporting the identification of several unsteady flow mechanisms (Ellington 1984;Birch and Dickinson 2001;Sane 2003;Spedding et al. 2003;Lehmann et al. 2005). From the analysis of both natural fliers (e.g., Spedding and Hedenström 2009) and mechanical systems, an understanding of the general flow characteristics of flapping wings and their relation to force generation has been obtained. Nevertheless, distinct features for specific animals/MAVs may arise from the particular wing arrangement and flapping kinematic patterns involved.
Different strategies may be distinguished regarding the experimental procedure followed in the aerodynamic investigation of flapping-wing propulsion. The most direct approach captures the animal (or MAV), as it passes in its natural flying mode through the observation area while performing measurement of the footprint of the vortex structures generated during the flapping flight (Langley et al. 2014;Gutierrez et al. 2016;del Estal Herrero et al. 2018). Although this puts relatively low demands on the test setup (no wind tunnel is needed, for example), it may require extensive training (animal) or accurate control (MAV) of the test object, to properly align the flight path with the observation area. Even then, the repeatability of the measurements is likely to be negatively affected by the inherent unpredictability of animal behaviour or the limited control authority on the MAV (del Estal Herrero et al. 2018). For these reasons, like in general aerospace research, the common practice to investigate flapping-wing flight has been to keep the test object at a fixed position relative to the laboratory frame (and hence the measurement setup) while representing the forward flight effect using a wind tunnel facility. Most often, the position and orientation of the test object (animal/MAV) are fixated with a mounting structure or a tether, preferably under conditions that are close to those of natural flight (Ellington et al. 1996;Srygley and Thomas 2002;Bomphrey et al. 2006;Nakata et al. 2011;Mazaheri and Ebrahimi 2011;Ren et al. 2013). This mounting results in a relatively simple setup, while, furthermore, the forces generated by the flapping wings can be measured with a balance simultaneously to the flow visualization. This allows valuable information to be obtained on the relation between flow behaviour and resulting force generation. When possible, one may alternatively approach the natural free-flight condition more closely by training/controlling the animal/ MAV to remain suspended without physical restraint in the test section (Spedding et al. 2003;Hedenström et al. 2006Hedenström et al. , 2009Johansson and Hedenström 2009;Muijres, et al. 2011;Karasek et al. 2019).

Flow characterization of flapping flight
Qualitative smoke visualisation has proved instrumental in revealing particular flow structures, such as the occurrence of leading edge vortices (LEVs) on insect wings, in, for example, hawkmoths (Ellington et al. 1996) and butterflies (Srygley and Thomas 2002). Planar particle image velocimetry (PIV) provided further quantitative information on the flow structures near the flapping wings as well as in the (near) wake, from observations on natural flyers and flapping mechanisms. Important aspects of flapping flight were thus quantified, such as the characterization of the LEV structure on insect-scale wings (Birch and Dickinson 2001) and the clap-and-fling interaction between two wings (Lehmann et al. 2005). These flow features were also found to occur in actual flyers, such as the vortical structures generated by birds (Hedenström et al. 2006) and flapping-wing MAVs (De Clercq et al. 2009). By combining the information obtained in several planes along the span, a reconstruction of the three-dimensional wake flow can be proposed (Spedding et al. 2003;Bomphrey et al. 2006;Groen et al. 2010;Ren et al. 2013). A more accurate representation of the 3D flow structure can be obtained from dense multi-plane 2C/3C data, by phase-synchronizing sequentially acquired planes, such as carried out for the flow generated by a flappingwing mechanism (Lu and Shen 2008;Carr et al. 2013) and also for reconstruction of the flow in the wake of the DelFly MAV (Percin et al. 2014(Percin et al. , 2017Deng and van Oudheusden 2016). Evidently, this procedure is laborious and, furthermore, depends critically on cycle repeatability, which does not make it likely a suitable procedure for animal studies, although it was successfully applied for the investigation of butterflies by Fuchiwaki et al. (2013).
With the availability of high-speed PIV, a simpler yet more approximate approach has become common for characterizing flapping-wing wakes. Here, a convection model based on Taylor's hypothesis is used to reconstruct the three-dimensional wake structure from the time-resolved data obtained in a single spanwise-oriented wake plane behind the flapping wings. This approach assumes that the flow structures that are generated by the flapping wings and measured in the wake plane are convected downstream without distortion and at a constant velocity, equal to that of the free stream, thus neglecting the diffusion of the wake structures as well as their deformation as a result of the induced velocities and wing motion. This procedure has been widely applied for the wake investigations of animals Johansson and Hedenström 2009;Hubel et al. 2010;Muijres et al. 2011) and also for MAVs (Percin et al. 2014;Deng and van Oudheusden 2016). This allowed to identify the specific wake structures of different animal species and their dependence on flight regime while further providing a means to quantify how the circulation and spanwise distribution of lift of the wings varies over the flapping cycle.
The accuracy of this wake reconstruction method evidently depends critically on the reliability of the convection model underlying the Taylor's hypothesis, and may be violated when the distorting influence of the induced velocities becomes significant. Bomphrey et al. (2012) used highspeed volumetric tomographic PIV to reconstruct the wake of a desert locust in tethered conditions. The measurement volume in the experiment (60 × 80 × 4 mm 3 ) had a limited thickness of 4 mm oriented in the longitudinal direction, to be able to cover the wake in width and height. This made the use of Taylor's hypothesis again necessary to propagate the flow information downstream to describe the full-wake length. Analysis of the results revealed that this procedure introduces errors on wake morphology, which may lead to erroneous force calculations based on such information. After correcting for this deformation by realigning the subsequent wake elements, the shape of the wake structures was found to be changed significantly and it was concluded that the volumes captured at subsequent time intervals do not reliably compose the actual wake shape. Instead, a large volumetric imaging technique would be desired to properly capture the longitudinal development of the wake.
A variety of volumetric flow measurement techniques have been reported in the literature, as reviewed recently by Discetti and Coletti (2018), but most techniques have important restrictions on the achievable volume, especially for experiments in air. Therefore, apart from generic flappingwing configurations (e.g., David et al. 2012;Cheng et al. 2014;van de Meerendonk et al. 2018), volumetric visualization studies of full-scale true flyers are quite scarce. Langley et al. (2014) applied synthetic aperture PIV (SA-PIV) for the investigation of the flow around a butterfly, in a single volume of size 330 cm 3 . Henningsson et al. (2015) performed a tomographic PIV (tomo-PIV) measurement of the wake of a desert locust in a domain of 200 × 240 × 58 cm 3 (length × span × height), corresponding to a volume of 2800 cm 3 , which allowed to capture the full extent of the wake cycle. To the authors' knowledge, this constitutes presently the largest single-volume flapping-wing visualization experiment in air.

Flow visualization research on the DelFly MAV
The DelFly MAV (Fig. 1) is a flapping-wing MAV concept with a tip-to-tip wingspan of 28 cm and a nominal mass of around 20 g. It features a biplane configuration, with wings moving alternately towards and away from each other when flapping (see Fig. 2) to generate both thrust and lift, taking advantage of the clap-and-fling wing-wing interaction effect, which in the case of flexible deforming wings is rather to be interpreted as a "clap-and-peel" (De Clercq et al. 2009;Armanini et al. 2016;Percin et al. 2016Percin et al. , 2017. of a DelFly model, flapping at 13 Hz in a quiescent environment representing a hover flight configuration, in a plane at 75% semi-span. The measurements provided a characterization of the wing deformations, including the occurrence of the clap-and-fling interaction. The presence of the LEV during instroke and outstroke was confirmed, although its visualization during most of the flapping motion was optically obstructed by the wings. To resolve this issue, Groen et al. (2010) performed phase-locked stereo-PIV measurements with data acquisition triggered when the wing surface is aligned normal to the illumination plane, for a given phase of the flapping motion. With this approach, the existence of an LEV during the flapping motion was explicitly revealed. Measurements carried out at different spanwise positions showed that the LEV had an approximately conical shape up to about 86% of the wingspan, where it was speculated to be connected to the tip vortex. To obtain more detailed three-dimensional information on the flow structure around the wings and in the near wake, phase-synchronized multiplane measurements were performed by (Percin et al. 2014(Percin et al. , 2017Deng and van Oudheusden 2016). In addition to the true spatial reconstruction of the wake, also the spatio-temporal wake reconstruction using Taylor's hypothesis was investigated, which may be assumed to be reasonable for a  (Percin et al. 2014;Deng and van Oudheusden 2016). This allowed to visualize the main structures in the wake, notably the tip vortices generated by the four wings and the starting trailing edge vortex connected to the tip vortex. With this approach, it was also possible to assess the effect of the reduced frequency on the interaction between the vortex structures. It was concluded that an increase in reduced frequency leads to a stronger interaction between the vortices, which is linked as well with the enhancement of the clap-and-fling mechanism. Comparing the spatio-temporal reconstruction with the true structure of the wake obtained from the spatial interpolation, very similar characteristics were observed. This was interpreted by the authors as a justification of the convection model, while some detail differences were attributed to the non-uniform convection velocity in the near wake. However, the use of Taylor's hypothesis fails to represent the wake at more realistic flight conditions, where the angle of attack is quite high and the interaction between the velocities induced by the wings and the free stream is too complex to be accounted for by a simple uniform convection model. Free-flight visualizations of the DelFly (del Estal Herrero et al. 2018;Karasek et al. 2019) give, furthermore, clear evidence that in slow forward flight, the wake is predominantly aligned along the body axis direction, making the assumption of convection in the free stream direction unrealistic.

Objectives of the present investigation
The above considerations motivate the (large-scale) volumetric imaging to properly capture the longitudinal development of the wake, under realistic conditions that are representative of actual flight. With the wing span of the Delfly being ~ 30 cm and typical wake structures ~ 20-25 cm in length (del Estal Herrero et al. 2018, see also Sect. 3.3), the minimum volume requirements would amount to around 20,000 cm 3 . This cannot be achieved by regular tomo-PIV, in view of the volume restrictions, which has been reported as being typically no more than 100 cm 3 for low-repetition-rate systems in air (Scarano 2013;Scarano et al. 2015). Exceptionally large tomo-PIV measurement volumes, using a powerful laser illumination of about 1 J/pulse, were reported in the studies of Fukuchi (2012) on a cylinder wake and Henningsson et al. (2015) on the wake of a desert locust, achieving volume sizes of around 2800 cm 3 . An important step towards substantially increasing the dimensions of the measurement volume has been achieved with the introduction of helium-filled soap bubbles (HFSB) as tracer particles, which are of sub-millimetre size and have a light-scattering capability of typically 10 4 times that of micron-sized droplet seeding (Scarano et al. 2015). Large-scale experiments using HFSB have been reported, amongst others, by Kühn et al. (2011) using tomo-PIV in a convection cell (56,000 cm 3 ) and Huhn et al. (2018) using Lagrangian particle tracking of a jet flow at 4 m/s (54,000 cm 3 ). Earlier experiments have shown that the use of HFSB is not detrimental to the operation of the DelFly (del Estal Herrero et al. 2018), although occasional cleaning of the wings is required, which, therefore, makes it a feasible option for the current purpose as well. In the present investigation, the particular strategy to achieve the required large-scale measurement volume to visualize the entire wake extent of the DelFly is combining the HFSB seeding with a robotic PIV system based on coaxial illumination and imaging Jux et al. 2018); more details are given in Sect. 2.3.
In addition to the realization of a full-wake visualization, a second objective of the current work is to compare the experiments carried out in a tethered configuration in a wind tunnel with free-flight visualization experiments. Earlier DelFly studies have compared the aerodynamic forces derived from flight tests to balance data taken under similar conditions in the wind tunnel (Caetano et al. 2015;Karasek et al. 2016). In general, good agreement was found, especially for the thrust force in the direction along the body longitudinal axis. However, larger differences were found in the wing stroke direction, which were attributed to structural vibrations of the MAV body, which affect the force balance measurements depending on the clamping position. For the free-flight measurements in the current study, the same procedure is followed as for the preliminary free-flight visualization experiments described in (Karasek et al. 2019), by combining the PIV visualization with a position-control system for autonomous flight in the wind tunnel environment.

Experimental procedures
The DelFly MAV was used as the experimental object in the present study. Wind tunnel experiments were performed in two different campaigns, one with the DelFly model in tethered conditions and another with the complete DelFly vehicle in actual free flight. For the first set of experiments, the selected configurations aim to represent conditions that the DelFly would have in free flight. This is achieved by setting appropriate values of the relevant parameters, which are the flight speed ( U ∞ ), the flapping frequency ( f flap ) and the angle of attack of the body ( b ) with respect to the free-stream direction, based on data from previous free-flight experiments (see Caetano et al. 2015;Karasek et al. 2016). The main parameters for each of the considered configurations are provided in Table 1, including the corresponding value of the reduced frequency k . The adopted nomenclature of the configurations is such that the first digit relates to the freestream flow speed, while TC and FC indicate the tethered or free-flight configuration, respectively. The experimental procedures of the investigation are presented in the following sections, addressing the MAV vehicle that has been used (Sect. 2.1), the different experimental facilities (Sect. 2.2) and the robotic PIV system (Sect. 2.3).

The DelFly MAV
The DelFly MAV used in the experiments is shown in Fig. 1, together with the layout and dimensions of one of the half-wings. The MAV has a configuration with two fullspan wings, arranged in cross configuration with a dihedral angle of 12°, that flap in counter motion. The wings have a span of 280 mm and mean chord of 80 mm. The wing kinematic motion is illustrated in Fig. 2, which shows that during one flapping cycle, the wings periodically move away and towards each other, which is referred to as the outstroke and instroke phases of the cycle, respectively. To assist the discussion of the results, the flapping cycle has been divided in eight phases (indicated as I, II, etc.) equally spaced in time over the cycle. The maximum angle between the wings (reached in the phase indicated as V) is 87°. The wings are made of 15 µm Mylar foil and are reinforced with carbonfibre leading edge spars and wing stiffeners. As a result, the wings are highly deformable, which is a characteristic that was found to enhance the aerodynamic performance significantly ). The DelFly is, furthermore, equipped with conventional horizontal and vertical tail planes and corresponding control surfaces, elevator, and rudder (see Fig. 1a), which make it a stable and controllable MAV platform relying on tail actuation. The weight of the DelFly is approximately 25 g and the distance from the leading edge of the wings to the trailing edge of the tail is about 215 mm. The tail has a span of 170 mm and is made from 2 mm-thick polystyrene sheet; it was painted black for the experiments to avoid unwanted reflections from its surfaces. This treatment was not applied to the wings of the DelFly, as the paint could change considerably the stiffness and the weight of the wings, and thus affect the flight behaviour. However, as the current interest lies primarily in the characterization of the wake structures, light reflections from the wings as well as their possible obstruction of the measurement domain are no direct concern.

Tethered experiments
The tethered experiments were performed in a low-speed wind tunnel at TU Delft, which has an open-exit test section with cross sectional dimensions of 60 × 60 cm 2 and a contraction ratio of 3.6. The measurements were performed at free-stream velocities ( U ∞ ) in the range of 0.5-2 m/s, with corresponding mean-chord-based Reynolds number ( Re = U ∞ c∕ ) ranging from 2800 to 11,300. The DelFly was positioned in the middle of the exit stream of the wind tunnel, supported by a structure that allows to set the desired body angle of attack (Fig. 3). A Hall-effect switch, installed next to the magnet-equipped output gear of the DelFly gearbox, provided the information about the wing position by delivering a pulse signal at every flapping cycle. By means of high-speed camera recordings, the system was calibrated to determine the position of the wings, and hence the phase in the flapping cycle, when the pulse signal was sent. The power input for the DelFly is regulated to achieve a constant flapping frequency, which is varied between 10 and 13 Hz, depending on the particular experiment (see Table 1).

Free-flight experiments
The second experimental campaign was performed in the open jet facility (OJF) of TU Delft, which is a relatively large-size wind tunnel with an open test section of 2.85 × 2.85 m 2 cross section, and a contraction ratio of 3.0. The main difference for these experiments is encountered with respect to the positioning of the DelFly, which in this case is flying in the open test section of the wind tunnel without any support (Fig. 4). The DelFly position is controlled by an on-board autopilot, which steers the vehicle based on feedback from the on-board IMU (used for attitude estimation) and from an external motion tracking system  that provides position and heading information with respect to ground (Karasek et al. 2019). The external motion tracking system employs 12 OptiTrack Flex 13 cameras (Natu-ralPoint, Inc., resolution 1280 px × 1024 px, 120 fps). The tracking system data are collected to determine the location of the DelFly as well as the instantaneous wing position throughout the measurement. This provides the information on the flapping phase of the wings required in the phaseaveraging procedure while also allowing to apply corrections on the measured data to account for the dynamic motion of the DelFly. More detailed information on the control system and the data synchronization with the PIV data set can be found in (Karasek et al. 2019).

Robotic PIV system
The investigation uses robotic volumetric velocimetry ) as visualization technique, which allows the measurement of all three components of the velocity in a large volume (of the order of several litres). The optical component of the measurement system is formed by the coaxial volumetric velocimeter (CVV, see Schneiders et al. (2018)) which makes use of a coaxial arrangement between imaging and illumination. It features a compact module containing four cameras and the laser light source (see Fig. 5).
In the robotic volumetric velocimetry system, the CVV module is mounted on a robotic arm to move it with millimetric precision. The advantages of this compact and robotically actuated device are notably the versatility in positioning and the time saving, as only one calibration is needed in view of the constant illumination and imaging configuration that is ensured by the CVV module. This facilitates combining subsequent measurements in adjacent volumes, such as to achieve an enlarged measurement domain. The robotic volumetric velocimetry system is used in combination with HFSB flow tracers (Scarano et al. 2015), which are sub-millimetre neutrally buoyant particles that scatter light with much higher intensity than conventional micron-sized droplet particles, thus allowing the volumetric visualization capability of the CVV approach to be exploited. The complete Robotic PIV system is, hence, composed of three major parts: the HFSB seeding generator, the CVV optical system, and the robotic arm (type Universal Robots-UR5) to which the CVV is mounted. The seeding generator consists of 10 wings with 20 seeding nozzles each, spanning a total cross section of 50 × 100 cm 2 and with a nozzle pitch between generators of 5 cm in both directions. Each nozzle produces neutrally buoyant helium-filled soap bubbles of 300 µm diameter at a nominal rate of ~ 30,000 bubbles per second (Faleiros et al. 2019). For the tethered experiments, the seeding generator is placed in the settling chamber of the wind tunnel to minimize flow intrusiveness. After the wind tunnel contraction, the seeded stream-tube covers a region of about 30 × 55 cm 2 around the DelFly. In the case of the free-flight experiments, the seeding generator is placed after the contraction, so that the seeded stream tube has a cross section of 50 × 100 cm 2 . The tracer concentration is between 1.5 and 2 tracers/cm 3 , resulting in a number of particles per pixel (ppp) between 0.07 and 0.09.
The CVV system is composed of four CMOS cameras (10-bit, 800 × 600 pixels, 4.4 µm pixel size, 471 Hz acquisition frequency at full resolution) assembled close together at low tomographic aperture ( 0 = 4 • ). The cameras mount lenses with a focal length f = 4 mm, whose optical aperture is set at f# = 8. The illumination is provided by a Quantronix Darwin Duo Nd:YLF laser (21 mJ/ pulse at 1 kHz, wavelength of 527 nm) via an optical fibre located in the middle of the cameras arrangement. The cameras recorded 8000 images per cone acquired at a rate of 471 Hz in continuous mode for all the measurements performed. The images need to be pre-processed before the particle-tracking algorithm can be applied due to the reflections generated by the body and wings of the DelFly. The approach followed to eliminate the unsteady reflections is to apply a Gaussian smoothing in space (kernel size of 3 × 3 pixels) to spread the reflection area along the different images, upon which it can be removed by a high-pass frequency Butterworth filter (Sciacchitano and Scarano 2014). The total measuring time per sequence is 17 s which corresponds to approximately 200 wingbeats, with 40 images recorded for each wingbeat, at an average flapping frequency of 12 Hz. The CVV system records images in a conical volume of about 25 cm in height and width, and 30 cm in depth. Five repeated measurements were taken for each conical volume, resulting in 40,000 images in total.
To capture the entire wake of the DelFly wings, for each configuration, multiple conical volumes are taken at different locations and the information obtained for the different cones is then combined together with an in-house Matlab code. The spatial merging of the different conical regions is achieved from the coordinate calibration that forms an inherent feature of the robotic system (see Jux et al. 2018). Four different conical volumes are taken in total, as depicted in the graphical representation in Fig. 6 (not to scale). The orientation of the volumes is selected to be along the longitudinal body axis, because under the present flow conditions, the wake alignment is dominated by the induced flow of the flapping wings del Estal Herrero et al. 2018). The depth direction of the cameras was placed along the spanwise direction of the DelFly, as shown in Figs. 3 and 4. The CVV was moved in the longitudinal plane, i.e., normal to the spanwise direction. One pair of cones was needed to visualize the upper and lower parts of the wake. In addition, two consecutive pairs of cones were taken along the body direction to capture the length of the wake. The dimensions of the combined measuring volume amounted to 50 × 30 × 40 cm 3 (in x, y, and z-directions, respectively), providing a total volume of 60,000 cm 3 .

Particle-tracking algorithm
The CVV system records the particles intensity in the images, and by means of particle-tracking velocimetry (PTV), the velocity information of the flow is extracted from the images. The information about the location of each particle in time is determined in a Lagrangian way, using the Shake-the-Box (STB) algorithm (Schanz et al. 2016) as implemented in LaVision DaVis version 8.4.0, to reconstruct the particle trajectories. With the information of the location of the particles in time, a second-order polynomial fit is used to compute the particle velocities. This technique is several orders of magnitude less demanding in terms of computational time than the tomographic PIV technique, due to the fact that the volume reconstruction is particle-based instead of voxel-based and that no cross-correlation analysis is conducted. The output of the Shake-the-Box processing is the information of the tracks that the particles have followed in time.
The number of particles captured at each time instant, and therefore the track density in the measurement domain, is usually not large enough to achieve instantaneous velocity fields with sufficient spatial resolution to resolve the main flow structures. Therefore, an averaging approach needs to be applied to increase the density of particle tracks in the measurement domain. In the present investigation, the periodic evolution of the flow is obtained by computing phaseaveraged flow fields that correspond to different phases in the flapping cycle. The temporal (i.e., phase) synchronization of the different measurement areas is achieved using the wing position as phase indicator, as explained in Sect. 2.2. This allows the velocity information obtained for every individual particle track to be uniquely associated with a location in the measurement domain and phase in the flapping cycle. Information about the wing position through time, correlated with the recording of the images, becomes crucial for the success of constructing the phase-averaged flow visualization. To achieve a better accuracy at each phase, the information from one-time step to the following one is interpolated to obtain the information at the exact time when the wings are at a specified location (i.e., phase). Finally, by means of a binning procedure (Agüera et al. 2016), the flow results that are originally obtained in a Lagrangian (unstructured) representation are converted into a Eulerian one to obtain a three-dimensional grid that relates to the measured volume. Before the spatial averaging (binning) is conducted, the particle track information is mapped onto a global coordinate system, such that no artefacts are generated at the boundaries of the overlapping conical regions. A cubic bin of linear size 30 mm was used with an overlap of 75% which results in a vector pitch of 7.5 mm.

Results
The current section reports the evolution of the wake vortex structures through the flapping cycle. To start, the velocity results for one of the configurations (3TC) are discussed in detail, to explain the main features of the wake topology of the DelFly. Wake planes that are oriented perpendicular to the body axis are addressed first (Sect. 3.1) and next the three-dimensional vortical topology (Sect. 3.2). Subsequently, the different configurations in tethered conditions are compared in Sect. 3.3 to analyse the effect of changing flight conditions and reduced frequency on the wake structures. Finally, in Sect. 3.4, a comparison is made between the results for tethered and free-flight configurations at similar conditions.

Wake topology during a flapping cycle
From the data measured for configuration 3TC, information on the longitudinal velocity is extracted in a single wake plane (Fig. 7). The wake plane is located at 10 mm downstream of the trailing edge of the tail and its particular position is indicated by the orange line in the top and side views depicted in Fig. 14. The information is given in the body reference frame, with the x-component set along the body of the DelFly and y the spanwise direction (see also Fig. 6). In the plots, the position of the DelFly is schematically indicated by a black outline to relate the flow field results with the position of the wings. The flapping cycle has been divided into eight phases to show the behaviour of the flow structures and their evolution over the flapping cycle (see Fig. 2). The flapping cycle is defined to start when the wings are closed and the outstroke is going to start (phase I). Note, however, that the flow measured in the observation plane is not corresponding to the instantaneous position of the physical wings at that precise moment, in view of the convective delay between the wing position and the visualization plane. Therefore, to assist the interpretation of these flow fields, a "virtual wing position" (hence, flapping phase) is introduced, which has been deduced from the information of the change in crossflow velocity along the wingbeat cycle. In the setup of these experiments, the CVV system was positioned at the starboard side of the DelFly (i.e., positive y, see Fig. 3) which leads to a better tracking of the particles on that side, because the particles that are closer to the laser source will scatter more light as the intensity decreases with the fourth power of the distance. As a consequence, having more particles tracked contributes to a higher accuracy in that measuring side.
Regarding the accelerated flow region, represented by the red areas of relatively high longitudinal velocity, its evolution along the flapping cycle can be observed in Fig. 7. The first phase (Fig. 7-I) of the cycle corresponds to the start of the outstroke. There is no clear acceleration region in the wake, for the reason that the wings are closed (see also Fig. 2). It is seen in Fig. 7-II that a high-velocity region has drastically increased in size with respect to the previous phase, due to the opening of the wings. The deformation of the wings plays a major role in this effect due to the "peeling" action. To illustrate the wing deformation during the clap-and-peel process, video images of the corresponding phases are depicted in Fig. 8, with the approximate chordwise shape of the two wings indicated by the orange lines.
The acceleration region keeps increasing during phases III and IV as the wings are still opening further. On phase V, when the instroke has started, the flow acceleration that is generated by the wings reaches its maximum, again due to the deformation of the wings and the "clap" effect that is created as the wings start moving together (see Fig. 8-left). During the following phases, from VI to VIII, the acceleration region reduces its size and it is concentrated into two different regions, one near the body of the DelFly and the other near the tip of the wings. The flow at the first location follows the direction of the body, while the second acceleration region has a lateral component of velocity directing the flow outward along the wings (Fig. 9). The outward velocity occurs as a consequence of the deformation of the wings. The deformation along the wing is restricted because of the presence of the stiffeners, as shown in Fig. 1; this may explain why the acceleration region at the end separates in two. In the most outward part of the wing, apart from the expected chordwise deformation, there is a spanwise Fig. 7 Velocity component in the body direction in a spanwise-oriented plane at 10 mm distance from the end of the tail at 8 phases along the flapping cycle (configuration 3TC) deformation of the wing which is the probable cause of the observed flow acceleration in the lateral direction. Figure 10 displays the longitudinal vorticity ω x (i.e., in the direction of the MAV body axis) for the subsequent phases of the flapping cycle. The first observation is that there are two main counter-rotating vortices originating from the tips of the DelFly during the entire flapping cycle. According to lifting line theory, this implies that positive lift is generated during the entire flapping cycle, which is crucial for the MAV to provide the necessary lift to maintain its altitude. Lift is defined as the force generated perpendicular to the free stream, so the vertical force in an absolute reference frame. The vortices are shed by the lower pair of wings during the outstroke and by the upper wings during the instroke. This finding leads to the conclusion that the main vortices are always shed by the wing that is moving in the direction opposite to that of the free stream when the relative velocity of the wing with respect to the flow is higher and the majority of the force is generated. Considering the wings moving along with the flow, the velocity experienced by the wing surface is reduced significantly. Assuming that the effective flapping-induced velocity produced by the wing is located at two-thirds of the wing span ( b = 14 cm ), it is estimated that the total velocity seen by the wing moving against the flow is more than three times larger than the one seen by the wing going in the free-stream direction.
Looking more in detail to the different phases, two strong counter-rotating tip vortices are retrieved in phase I of the cycle (starting of the outstroke) that have been generated by the upper wings during the previous instroke. These tip vortices (wake structures A in Fig. 10) interact with the wings while they are closed, and then with the tip vortices generated by the lower wings during the outstroke, resulting in a major tip vortex structure with the same sign for the entire cycle, as it can be clearly appreciated in the 3D views that will be shown in Sect. 3.2. The tip vortex structure in phase II is quite elongated, as is also the case for phase III. The shape of the tip vortices along the outstroke will be further discussed in Sect. 3.2. Several smaller vortices, indicated by the label B, start to appear in phase II, but they are not yet defined until phase III. The crossflow going over the wings from tip to root needs to be redirected once it reaches the root, generating these vortices.
In phase IV, the wings are completely open, indicating the transition from the outstroke to the instroke. According to wing theory, a lifting surface of finite span will generate tip and root vortices at the extremities of the surface, that are aligned with the free-stream direction. This same concept can be applied to flapping wings, similar to propeller theory. In this way when the upper wings start to generate lift, a spanwise-aligned starting vortex is produced that is followed by the generation of counter-rotating tip and root vortices. From Fig. 9, phases IV and V, the evolution of root and tip vortices can be confirmed (wake structures C and A, respectively). Visualization of the starting vortex is performed by analysis of the spanwise vorticity. In Fig. 11, such a plot is presented for phase V, where the starting vortex can be seen at the upper side of the DelFly. The vorticity in that region has a negative sign which means that positive circulation, and therefore, positive lift is being generated. The starting vortex can only be clearly seen at the beginning of the instroke, because at the beginning of the outstroke, a strong interaction between the wings occurs, so the change in circulation of the wings is not significant. On the lower side of Fig. 11, however, the opposite behaviour is found in a large area with high positive y-vorticity. It corresponds to the trailing edge vortices shed from the lower wings when they stop producing lift. The trailing edge vortex generated by a lifting surface is indeed of opposite sign to that of the starting vortex.
In phase V, the wings are completely open, which causes the crossflow along the lower wings to be maximum; hence, the secondary vortices formed on the upper part of the Del-Fly are enhanced. Then the peak vorticity of these vortices diminishes drastically when the wings come closer together, as the flow is now directed away from the root, and towards the tips. In the same way, the root vortices also experience a decrease in strength and they nearly disappear at phase VIII, as the wings are approximately closed and the DelFly resembles a fixed wing vehicle. The tip vortices, where the highest flapping-induced velocity can be found, still maintain their strength throughout the final phases of the cycle.

Three-dimensional wake topology
Considering again configuration 3TC, the vortex structures identified in the previous section are now discussed by means of their three-dimensional representation. First, a complete three-dimensional visualization of the entire wake is shown in Fig. 12, for one specific phase (an animation for the full cycle is provided as online supplementary material). It can be observed that the wake displays a good level of symmetry with respect to the centreline plane of the DelFly. Therefore, although data are available for the full span of the wake, in the following discussions only the structures at one side (corresponding to positive y) are shown, to facilitate visual interpretation.
The development of the topology obtained for different phases of the cycle is sketched to give an overview of the variation of the structures, while they are convected. In Fig. 13, the different vortices are identified, as well as the distinction between instroke and outstroke, indicated by the grey/white colour bands in the figure. The latter division has been obtained from visual inspection of the vortical shapes in accordance with the reasoning presented in Sect. 3.1. It may be reminded here again that the x-axis is aligned with the body axis, so the upward deflection of the wake that can be observed in the side view is the result of the body inclination with respect to the free-stream direction (body angle b = 31 • in this case).
The large tip vortex structures that were already mentioned in Sect. 3.1 can be clearly recognized in Fig. 13. The shape of the convected tip vortex seems to be straight with respect to the DelFly body during the instroke, while it is deflected downwards during the outstroke, which was also seen in Fig. 10 on phases II and III. During the outstroke, the lower wing is the one that is generating the tip vortex, which agrees with the downward movement of the vortex. Another reason for the shape of the vortex on the outstroke comes from the fact that the convection velocity is lower in this part of the cycle (Fig. 7), especially in the tip area, where the air is being redirected towards the root. Therefore, the vortex does not travel away as rapidly along the x-direction in the first stages of the outstroke. In addition, the vortex is moved slightly towards the centre of the wake due to the crossflow that appears. This also explains that the instroke structures are longer than the ones generated during the outstroke (the grey area is shorter than the white area), similar as in Percin et al. (2014).
The wake topology is highly affected by the different convection velocities that are experienced throughout the flapping cycle in different locations in the wake (Fig. 7). This makes the analysis of the wake evolution more complex. The vortices generated by the flapping wings are not convected in the direction of the body of the DelFly, neither in the direction of the freestream. There is a combination between the freestream and the effect of the varying induced flow generated throughout the cycle. The result is that the angle of the wake varies as the structures travel away from the wings. With downstream distance the importance of the induced flow acceleration decreases, and the wake convection velocity will become more uniform and closer to the free-stream velocity. However, this situation is not reached within the measuring volume studied and, instead, it will be assumed that over the measured region the angle of the wake is approximately constant. The value of this wake angle is computed by taking the information of the wake at two different spanwise-oriented planes that are spaced in the x-direction, such that they correspond to the same phase of subsequent flapping cycles. A correlation between the x-vorticity values of both planes provides a measure of the average displacement of the vortex structures during the length of a flapping cycle. This procedure is carried out for each of the eight phases, as shown in Sect. 3.1. The result for each phase varies, so only the correlations with a signal to noise ratio SNR > 2.9 were considered. For the current configuration (3TC), the mean angle of the wake with respect to the longitudinal body direction is established to be approximately 15°, which agrees with a visual inspection of the side view in Fig. 13, and clearly different from the body angle (31°). The angle of the wake in the lateral direction is not so easily verified visually, but the crosscorrelation procedure computes an outwards movement of the wake structures at an angle of 2°. The outward velocity that was presented in the previous section has an influence on the movement of the wake, which slightly moves away from the centre plane.
The evolution of the wake over three different phases during the flapping cycle can be seen in Fig. 14, where the wake plane location corresponding to the information depicted in Sect. 3.1 is also indicated (orange line). The assumption of a constant wake angle seems to be reasonable when looking at the side views for the different phases. Starting from phase II that has been explained in detail before, phase V shows clearly the starting vortex linked to the root and tip vortices during the instroke. The secondary vortex created during the outstroke moves upwards and inwards, as visible in, respectively, the side and top views. As can be expected, the vortex structures become more diffuse, as they convect Fig. 12 Perspective view of the full wake of configuration 3TC, vortical structures are visualized by isosurfaces of Q = 600 s −2 and coloured by ω x (animation is provided as online supplementary material) downstream, which can be seen when looking at the root vortex for phase II, which is not so prominent as in phase V. The peak vorticity of all vortex structures also decreases. On phase VII, the vortex structures from the instroke have moved downstream and the trailing edge vortex from the lower wings is seen below the root vortex. The secondary vortex seems to interact with small vortices generated before. This interaction probably originates before this point, but the reason behind this interaction was not studied in detail.

Effect of reduced frequency
Previous studies (Percin et al. 2014) have indicated that the interaction between the vortices in the wake of a flappingwing system is affected by the reduced frequency. However, in their investigation, the body angle was set to zero, resulting in symmetric conditions for top and bottom wings, while the investigated range of reduced frequencies spanned from 0.47 to 1.17. As such, the conditions investigated by Percin et al. (2014) are not representative of realistic slow forward flight, in contrast to the conditions considered in the present study (see Table 1). To assess the impact of the flight condition on the wake structure, the three investigated configurations (1TC, 2TC, and 3TC) are compared, where the reduced frequency of configuration 1TC is five times higher than that of configuration 3TC, being 6.5 and 1.3, respectively. It should be noted, however, that in this case, the configurations not only differ in the value of the reduced frequency, but also in the body inclination with respect to The 3D wake topology for each of the configurations is presented in Fig. 15, for a specific phase in the cycle. An animation of the entire cycle can be found in the supplementary material available online. There are large differences between the topologies shown, but some similarities can still be observed. The major vortex structure appreciated in all three configurations is the tip vortex. The tip vortex in case 1TC is no longer a smooth structure throughout the flapping cycle, but it has a distinct shape for instroke and outstroke. During the instroke, the tip vortex and the starting vortex cannot be distinguished, forming a single oblique vortex that starts close to the root of the DelFly and moves outwards as the instroke develops. The tip vortex generated during the outstroke is shorter than for the instroke, as it was the case for configuration 3TC, and it follows a straighter path with respect to the body of the MAV. A second tip vortex structure generated during the instroke on the lower part of the DelFly can be seen for configuration 1TC. This indicates that the lower wing is producing a vertical force. However, the x-vorticity of this structure is negative, which indicates that negative lift is being generated. This can be confirmed from the analysis of the y-vorticity which has positive values for the lower wing during the instroke (figure not shown as it is similar to Fig. 11). However, in this case, it is not the trailing vortex that is shed but the starting vortex related to a sudden change in lift generation. The strength of this vortex is smaller than the one measured for the tip vortex of the upper wing, which can be understood from their different relative velocities with respect to the flow, so the total resulting lift force remains positive. Another vortex structure that appears for configuration 1TC is located on the upper part of the instroke tip vortex. This vortex is not present close to the wings, but as the wake travels downstream, the vortices interact strongly with each other. The outstroke tip vortex seems to interact with the tip vortex of the following instroke. When looking at configuration 2TC, the wake topology is more similar to the one of configuration 3TC. There is a clear tip vortex, while a starting vortex can also be seen. The shape, however, is more curved rather that straight as in the case of configuration 3TC, which means that the wake of configuration 2TC starts to look more similar to that of configuration 1TC.
No root vortex is observed for configurations 1TC and 2TC for any phase of the flapping cycle. Due to the relatively high flapping frequency, the flow at the root of the wings is constantly changing and the dominant component is the crossflow. Therefore, it is conjectured that the root vortex is not present, as it is not generated as fast as the change in wing conditions happens. As the crossflow is augmented, the secondary vortex is strengthened, while also its length changes from one configuration to the other. It is Fig. 15 Three-dimensional views of the wake structures visualized by isosurfaces of Q = 500 s −2 and coloured by ω x for phases III of configurations 1TC (top), 2TC (middle) and 3TC (bottom) hypothesized that the deformation of the wings is an important contributing factor in this. As the flapping frequency increases, the deformation of the wings is enhanced and the phase lag between the leading edge and the trailing edge of the wings increases . This allows for the secondary vortex to last for a longer time, as the wings do not close completely until later in the flapping cycle. The secondary vortex seems to be generated closer to the root on configurations 1TC and 2TC than in configuration 3TC in both y-and z-directions. The fact that it is closer in the z-direction is likely due to the reduction of the free-stream component, so the convection of the vortices is much more influenced by the flow acceleration induced by the wings. This causes the angle of the wake with respect to the body to be smaller.
Although the observed wake topologies qualitative confirm these conclusions, an accurate quantitative characterisation can only partly be achieved, as the quality of the wake reconstructions reduces when moving from configuration 3TC towards 1TC. A number of characteristic parameters of the wake have been collected in Table 2: the wake length corresponding to a single flapping cycle (L), the maximum longitudinal velocity (U b MAX ) , and the maximum longitudinal vorticity ( x MAX ), both defined in the body frame of reference. The wake length of a flapping cycle has been computed by determining the difference between the location of the start of the instroke, which is defined as the time when the starting vortex is shed. This length is retrieved by visual means from the flow visualization, and should be treated as an approximate value. The values of maximum velocity and vorticity are computed, selecting only regions where more than 1200 particles have been tracked on the same slice as the contour plots, as shown in Sect. 3.1.
Based on the wake cycle length L observed in the visualizations, a mean convective velocity has been calculated, according to: U conv = Lf flap . As can be seen in Table 2, the effective convection velocity is much different from the free-stream speed, being substantially higher due to the flapping-induced flow. For configuration 3TC, the maximum induced velocity, U b MAX , is nearly twice the free-stream velocity, while the estimated mean convection speed is about 40% higher than the free-stream velocity, and close to the average value of U b MAX and U ∞ . The differences between the convection velocity and free-stream velocity become progressively pronounced for decreasing free speed (increasing reduced frequency), their ratio becoming about 2 (for 2TC) to 4 (for 1TC). Together with the wake alignment observed in Fig. 15, these results clearly illustrate the invalidity of Taylor's hypothesis for low flight speeds, which combine a high value of the reduced frequency with a large body angle. For comparison, also the reduced frequency k conv based on the convective velocity has been computed, which is seen to display a much smaller variation than the reduced frequency based on the free-stream speed. The maximum values of the longitudinal vorticity ( x ) occur always in the tip vortex structure across all phases analysed. The highest maximum value occurs for configuration 3TC, as more force is generated on the wings than for the other configurations. The overall maximum value along the flapping cycle is found at the beginning of the outstroke, when the wings are close together and coming apart from each other, which can relate to the interaction of the tip vortex from the previous instroke that combines its strength with the following vortex generated by the lower wing.

Free-flight vs. tethered conditions
To finalize the analysis, the comparison between tethered and free-flight experimental conditions is investigated using configurations 3TC and 3FC. The reduced frequency for both cases is very similar, albeit that the flapping frequency for case 3FC is about 10% higher than for case 3TC. This affects the length of the wake for a flapping cycle, as seen in Fig. 16, yielding approximately 27 cm and 24 cm for configurations 3TC and 3FC, respectively. Such result is consistent with the difference in flapping frequency, as very similar values of the convection velocity are obtained (see Table 2). The vortex structures that can be appreciated for free-flight are not captured with the same level of detail as for the tethered condition, and only the main structure, namely, the tip vortex, can be clearly seen. The root vortex is one of the structures that are not visible for configuration 3FC, whereas they are well defined for configuration 3TC.
The dynamic motion of the DelFly is inherent to the free-flight tests, as it is flying around a set point with small deviations from it. The largest deviations are found in the  Using the information from the flight tracking system, the position deviations are corrected for, to be able to process the data with the phase-averaging procedure. However, the motion of the DelFly, especially the lateral component which also involves the heading of the Delfly with respect to the free stream, changes the precise position and the orientation of the wake and to some extent the way the vortex structures are shed. This translates in fluctuations of the position of the individual wake structures that cause the smaller structures to be spread out over the volume, such that they are no longer detectable in the phase-averaged visualisation. The larger structures still remain visible but, as seen in Fig. 16, the core is larger due to this "blurring" of the vortex and the x-vorticity reached is smaller than for tethered conditions. The overall shape of the tip vortex, however, is still very similar to the vortex structure found in tethered conditions, indicating that the structures generated in both cases follow the same trend.

Conclusions
A large-scale volumetric flow visualization of the wake of the flapping-wing micro air vehicle "DelFly" was performed using robotic coaxial volumetric velocimetry in combination with helium-filled soap bubbles as tracers. A total measurement volume of about 50 cm × 30 cm × 40 cm (60,000 cm 3 ) was achieved, allowing the wake to be captured in both longitudinal and lateral extent. As such, these experiments constitute the first direct volumetric characterization of the wake of an FWMAV of this size. A phase-averaging procedure was applied to extract the wake structure at different phases along the flapping cycle. Experiments were carried out for different configurations, in tethered as well as free-flying conditions. The tethered experiments simulate free-flight scenarios along the envelope of the MAV in forward flight, by setting appropriate values of the flight parameters (flight speed, flapping frequency and body angle). This allows to compare the wake topologies at different conditions and to investigate the influence of the reduced frequency in particular.
The dominant vortex structures generated by the flapping wings have been identified to be the tip vortices, which form a continuous shape along the complete flapping cycle. The tip vortex starts to develop at the beginning of the instroke at the upper wings, and as the outstroke commences, it is maintained and deflected downwards by the lower wings. The fact that a continuous tip vortex structure is generated indicates that lift is produced throughout the complete flapping cycle. The creation, interaction, and development of the tip vortex and the other wake structures are strongly affected by the flow acceleration produced by the flapping wings. Due to the interaction between the freestream and the induced flow, the angle at which the wake is oriented is a combination of the body direction and the direction of the freestream, which will depend on the relative importance of each component. As such, the wake convection does not follow the velocity and direction of the freestream and Taylor's hypothesis was found to be (highly) invalid, notably for the conditions corresponding to the lowest flight velocities.
The length of the wake generated during the instroke is slightly longer than the ones generated during the outstroke because of the larger flow acceleration during the instroke phase. The deformation and the shape of the tip vortex depends to a great extent on the reduced frequency. When the reduced frequency is low (fast forward flight), the tip vortex is quite straight and aligned in the streamwise direction. With increase of the reduced frequency (decreasing flight speed), an oblique vortex structure can be seen as the tip vortex cannot be clearly distinguished from the starting vortex. The starting vortex is generated at the beginning of the instroke and is predominantly oriented perpendicular to the plane of symmetry. The close interaction between the tip vortex (streamwise direction) and the starting vortex (lateral direction) results in the formation of this oblique vortex. Due to the deformation of the wings, additional vortical structures appear during the final part of the outstroke and the beginning of the instroke. The crossflow going from the tip to the root of the wings along the span needs to be redirected causing secondary vortices to appear. The wings deformed during the outstroke creating a space that does not close until the instroke has already started, which allows the vortices to last until this gap disappears. When the reduced frequency increases, the deformation of the wings becomes stronger, and thus, the secondary vortices also become more prominent.
The second challenge addressed in this investigation was the flow visualization under free-flight conditions. Although the complexity involved in such measurements increased considerably, due to the unsteadiness of the DelFly position in the test area, meaningful results could, nevertheless, be achieved. Due to this inherent dynamic behaviour, however, only the larger and stronger vortex structures were captured. In addition, the vortex structures located closer to the body of the DelFly could not be visualized correctly. The shape of the tip vortex, being again the main vortex structure, and the traces of other structures suggest a very similar topology as the one found in tethered conditions.
In view of the scarcity of comparable volumetric flapping-wing studies in the literature, a tentative comparison is made to the 3D tomo-PIV visualization carried out by Henningsson et al. (2015) of the desert locust at a realistic flight condition. However, as the wing configuration of the desert locust (tandem wings) is quite different from that of the Del-Fly (counter-flapping wings), only a partial agreement can be expected. Common to both studies, strong tip vortices are observed as the main feature of the wake. In addition, for the locust, a more or less continuous and relatively linear tip vortex structure is observed to occur during an extensive part of the flapping cycle, in this case due to the downstroke of forewing and hindwing, respectively. In the DelFly, a similar prolonged tip vortex structure occurs, but then in relation to the top and bottom wings (see Sect. 3.2). For the locust, the hindwing tip vortex is connected to an arcshaped structure formed by starting vortex and root vortex, and for the DelFly, a comparable structure is observed at the start of the instroke (see Fig. 13). In the locust study, more additional small-scale structures are revealed, in view of the higher spatial resolution compared to the current experiment (approximately number of vectors per wing span length for the two experiments is 100 and 35, respectively).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.