An integrated measurement approach for the determination of the aerodynamic loads and structural motion for unsteady airfoils

The structural motion and unsteady aerodynamic loads of a pitching airfoil model that features an actuated trailing edge flap are determined experimentally using a single measurement and data processing system. This integrated approach provides an alternative to the coordinated use of multiple measurement systems for simultaneous position and flow field measurements in large-scale fluid–structure interaction experi- ments. The measurements in this study are performed with a robotic PIV system using Lagrangian particle tracking. Flow field measurements are obtained by seeding the flow with helium-filled soap bubbles, while the structural measurements are performed by tracking fiducial markers on the model surface. The unsteady position and flap deflection of the airfoil model are determined from the marker tracking data by fitting a rigid body model, that accounts for the motion degrees of freedom of the airfoil model, to the measurements. For the determination of the unsteady aerodynamic loads (lift and pitching moment) from the flow field measurements, two different approaches are evaluated, that are both based on unsteady potential flow and thin airfoil theory. These methods facilitate an efficient non-intrusive load determination on unsteady airfoils and produce results that are in good agreement with reference measurements from pressure transducers.


Introduction
The characterization of unsteady fluid-structure interactions (FSIs) is relevant for the design of deformable lifting structures in unsteady conditions, such as flexible flapping wings (Platzer et al., 2008;Heathcote and Gursul, 2007;Bleischwitz et al., 2017), rotor blades with unsteady periodic inflow (Wei et al., 2019;Strangfeld et al., 2016), or fixedwing aircraft subjected to gusts (Perrotta and Jones, 2017;Tang et al., 2010). The ongoing development of sophisticated computational methods for simulating these FSI problems requires validation data from wind tunnel experiments. This data is challenging to obtain because it requires the simultaneous measurement of a variety of quantities of different nature, which may explain why only a limited number of parameters has been typically used for validation purposes in the past (Tang and Dowell, 2001). The complete experimental characterization of a fluid-structure interaction requires the measurement of the structural response of the model as well as the unsteady flow around it. A particular focus is placed on the determination of the instantaneous geometry and aerodynamic load, which serve as coupling interfaces between fluid and structure in numerical FSI solvers (Kamakoti and Shyy, 2004).

Nomenclature c
Airfoil chord length, m C ℓ Section lift coefficient, C m Section pitching moment coefficient, Optical measurement techniques have received considerable attention in recent experimental FSI research because they do not require instrumentation of the model with installed sensors, which can be difficult or even impossible, for instance when investigating micro air vehicles with thin, lightweight, and/or highly flexible wings. However, until now, the optical measurement of the structural and the aerodynamic response in FSI experiments requires the simultaneous application of either multiple diagnostic techniques (Marimon Giovannetti et al., 2017;Zhang et al., 2019), or multiple data processing tools (Pereira Gomes and Lienhart, 2010;Liu et al., 2011). This complicates the experimental setup and increases the data processing effort and therefore limits the amount of information that can be extracted from the experiment. The present study aims to relieve this complication by using a single data acquisition and processing system to obtain non-intrusive position measurements of an unsteady airfoil with actuated flap as well as measurements of the unsteady flow and the aerodynamic response in terms of the lift, pitching moment and surface pressure distribution.
The state-of-the-art technique for unsteady flow field measurements is Particle Image Velocimetry (PIV, Raffel et al., 2018). The measurements in this study are obtained with a robotic volumetric PIV system , a recent development in PIV technology which facilitates large-scale volumetric flow measurements by effortlessly combining several flow field measurements from a compact Coaxial Volumetric Velocimeter (CVV, Schneiders et al., 2018) mounted on a robotic arm. The large-scale experiments are enabled by using helium-filled soap bubbles (HFSB, Scarano et al., 2015) as high-intensity flow tracers and the efficient Lagrangian particle tracking algorithm Shake-The-Box (STB, Schanz et al., 2016) for image processing. A limitation of this technique is the relatively low flow tracer density, which is typically on the order of one particle per cubic centimeter (Caridi et al., 2016). To augment the information density of the measurement, a phase-averaged approach can be adopted when conducting flow field measurements of (quasi-)periodic flow phenomena, as demonstrated by Martínez Gallar et al. (2019).
Of particular relevance in aeroelastic investigations of wings is the determination of the unsteady lift force. The classical approach to determine the force on an object from flow field (velocity) measurements, as introduced by Unal et al. (1997), is to apply the momentum equation in its integral form to a control volume enclosing the object. This approach has been applied successfully in several studies of airfoils in steady and unsteady flow conditions (van Oudheusden et al., 2007;David et al., 2009;Villegas and Diez, 2014a;Gharali and Johnson, 2014), but it requires the explicit determination of the pressure on the control volume boundary, which is not directly measurable with PIV. To avoid the source of error that this additional data processing step implicates, a vorticity-based formulation of the momentum equation was derived by Noca et al. (1997Noca et al. ( , 1999, where the pressure does not appear explicitly. This method has been applied in unsteady airfoil studies as well (Ferreira et al., 2011;Sterenborg et al., 2014b), however, its application requires the evaluation of several additional terms and gradient operations, hence, it is not evident if this method is superior (Rival and van Oudheusden, 2017).
An alternative method to determine the lift from the flow field is by application of the Kutta-Joukowski circulation theorem (see e.g. Anderson Jr. (2011), chapter 3). As this method requires no spatial derivation of the PIV velocity measurements, it is less prone to measurement errors than the control volume methods. The lift determination using PIV measurements based on the Kutta-Joukowski circulation theorem has been applied successfully in different studies of static airfoils (Lee and Su, 2012;Lind et al., 2014). Although the circulation theorem was originally derived under the conditions of steady, incompressible and irrotational flow, it was found to yield results that are in good agreement with validation data also for flows with separation and on unsteady airfoils (Sharma and Deshpande, 2012;Sterenborg et al., 2014a). Yuan and Olinger (2005) showed that the lift determination for an unsteady airfoil based on experimental circulation measurements can be improved compared to the quasi-steady formulation by adopting an unsteady approach, which accounts for the flow acceleration effects. The present study adopts a similar approach to determine the unsteady lift from circulation measurements.
Another load-related parameter of relevance in the aeroelastic analysis is the pitching moment. A circulation-based approach for lift determination does not permit a direct measurement of the pitching moment. In this study, it is proposed to determine the pitching moment by estimating the center of pressure of the lift force, which is derived from the circulation distribution along the airfoil chord.
A more direct approach to measure the pitching moment and the lift is to determine the load distribution by calculating the surface pressure from the PIV flow field data. The most popular method for the PIV-based pressure determination is by solving the Poisson equation for the pressure, after the pressure gradient in the flow field has been calculated from the velocity measurements using the momentum equation in its differential form (van Oudheusden, 2013). With this approach, it is possible to obtain the surface pressure distribution on objects immersed in the flow (Fujisawa et al., 2005) by prescribing Dirichlet boundary conditions for the pressure at some distance from the object, where Bernoulli's equation is typically invoked in regions not affected by rotation effects (Kurtulus et al., 2006). It was shown by Villegas and Diez (2014b), however, that the use of the Bernoulli equation in a quasi-steady sense is unsuited to determine the loads on an unsteady airfoil using this approach. Ragni et al. (2009) demonstrated (for a static airfoil) that a good estimation of the surface pressure can also be obtained using only velocity measurements near the airfoil surface. A similar approach is followed in this study, where relations from unsteady thin airfoil theory are applied to determine the surface pressure, including flow acceleration effects.
For the investigation of the structural dynamics of moving and deforming wings, Digital Image Correlation (DIC, Pan, 2018) and photogrammetric point tracking (Baqersad et al., 2017) are common non-intrusive approaches to replace the use of installed strain gauges and accelerometers (Black et al., 2010). The DIC technique provides a high measurement density of the structural displacement and deformation but requires the application of a dense speckle pattern on the investigated surface. Point tracking techniques, in contrast, provide discrete displacement information at locations where camera targets are placed. Because of the reduced model preparation effort, this approach can be viewed as more practical for investigating the large-scale structural motion of aircraft wings, where the observed deformations are relatively gradual so that a high spatial resolution of the measurement is not required (Liu et al., 2012).
Recently, Mitrotta et al. (2019) suggested that the STB algorithm is suitable for tracking circular structural markers in PIV images, and can therefore be used as a photogrammetric point tracking technique. In their study, marker tracking with STB was performed on the same PIV images that were used for the flow field analysis, showing a good agreement of STB marker position measurements and reference data obtained using a laser Doppler vibrometer. In the present study, the same approach is followed to obtain simultaneous measurements of the unsteady flow field and structural marker positions on a pitching airfoil with actuated flap, with a single optical data acquisition system. Subsequently, the airfoil position and flap deflection are determined from the measured STB marker positions using a regression approach based on the two degrees of freedom of the rigid-body airfoil/flap motion.

Wind tunnel setup
The experiments were conducted in the Open Jet Facility (OJF) at Delft University of Technology. The OJF is an open test section, closed return wind tunnel with an octagonal outlet (2.85 × 2.85 m 2 ), a maximum freestream velocity of 35 m/s  and a nominal freestream turbulence intensity of 0.5%. The wind tunnel setup is shown in Fig. 1, with an indication of all relevant components.
For the experiment, the investigated airfoil model is equipped with two circular end plates, where the top plate is attached to a support frame and the bottom plate is mounted on a steel frame that also holds the pitching mechanism. The latter controls the angle of attack oscillations α(t) of the model around the pitch axis, which is located at x r /c = 0.3125. The pitching motion is defined as: α(t) = α m + α a sin(2π t/T ), where α m is the mean angle, α a is the amplitude and T is the period of the pitching motion. The setup is mounted on a rotation table that allows to set the mean angle. The main parameters of the experimental setup are summarized in Table 1.

Airfoil model
The experimental model is a rigid rectangular wing with a NACA 0018 airfoil featuring a movable trailing edge flap that is hinged at x/c = 0.75. The airfoil model geometry is illustrated in Fig. 2.
The model is equipped with 24 Honeywell HSCS-RRN-1.6MD-SA5 differential pressure transducers (measurement range ±160 Pa, nominal accuracy ±4 Pa) installed inside the model, that are connected to 48 surface pressure taps. The taps are located along a diagonal line oriented at 15 • with respect to the inflow direction to reduce the effect of the upstream sensors on the downstream measurements. The sensors are used to directly measure the pressure difference between the model top and bottom surface at specific chordwise positions. To prevent laminar boundary layer separation near the airfoil leading edge, a zig-zag turbulator strip (thickness 2 mm, width 12 mm, pitch 60 • ) is installed on both sides of the model at x/c = 0.1. The airfoil model is identical to that used in a study by Raiola et al. (2018), to which the reader is referred for a more detailed discussion of the model design and instrumentation.
For this study, the airfoil model is painted black with a rectangular grid of 56 white markers (diameter 1.5 mm, spacing 35 mm) located between 0.35 < x/c < 1 and below the line of pressure sensors at mid-span, as shown in Fig. 2. The location of the pitch axis and the flap hinge are indicated in Fig. 2 as well, along with the positive direction of rotation for the angle of attack α and the flap deflection angle δ. The design of the present test model has been motivated in relation to the general framework of the ''Smart Rotor'' concept (Barlas and van Kuik, 2010), which aims at reducing construction, maintenance, and operation costs of wind energy systems by reducing the peak loads in wind turbine operation via active shape modification of the rotor blades. The active load reduction capabilities of the setup in this study are provided by a controller that acts based on the aerodynamic load measured using only 8 of the 24 differential pressure sensors to improve the response time, as discussed in detail in Fernández Barrio et al. (2020). The sensor signals measured are integrated over the chord using a trapezoidal integration scheme to provide an estimate of the instantaneous lift force. Based on the measured lift, the controller sends a signal to a JR DS8711HV servomotor, that is mounted on the bottom end plate and attached to the movable trailing edge flap, to counteract the unsteady loads induced by the pitching of the model, so as to maintain a constant value of the lift.
Results from two different experimental test cases are analyzed in this study, one where the flap is fixed with respect to the airfoil and one with actuated trailing edge flap. The lift coefficient over the period, calculated from the phase-averaged differential pressure transducer signals over 45 periods, is shown in Fig. 3 for the pitching airfoil model with fixed flap and with the flap actuated by the controller. Error bars indicate the standard deviation over the 45 periods, as obtained during the phase-averaging procedure. With the controller settings employed in this study, a reduction in the magnitude of the unsteady loads of about 80% is achieved, while the standard deviation at each phase is increased by about 40% due to the period-to-period variation in the flap motion, that is introduced by the controller as a result of the measurement noise in the pressure sensors signals.
Flow field measurements with PIV are conducted for both test cases. For the test case with actuated flap, a separate test was performed on the pitching airfoil model in absence of PIV measurements, where the controller output, which is the servomotor input signal, is stored. The recorded controller output is then phase-averaged and imposed on the servomotor during the PIV experiments, which means that when the PIV measurements are acquired for the test case with actuated flap, the controller itself is not active. This procedure improves the periodicity of the flap motion and prevents spurious flap deflection behavior when the pressure taps would become contaminated by the PIV seeding.
The pressure transducer measurements from both test cases are used as a reference for the PIV-based load determination. The procedure to generate this reference data from the pressure transducer measurements is described in Appendix A.

PIV measurements
The PIV system used for image recording is a LaVision Minishaker CVV probe mounted on a Universal Robots UR5 robotic arm, as described by Jux et al. (2018). The flow is seeded with HFSB that are generated by a 500 mm × 1000 mm seeding rake with 200 bubble-producing nozzles, which nominally produce 30 000 bubbles per second each, as described by Faleiros et al. (2019). The seeding rake was placed upstream in the wind tunnel settling chamber and the effective cross-section of the seeded flow region was ca. 200 mm × 200 mm. The illumination of the HFSB flow tracers was achieved  using two LaVision LED-Flashlight 300 illumination units mounted side by side to illuminate a chordwise region of 800 mm length, corresponding to −0.5 < x/c < 1.5 of the airfoil chord. This results in a total measurement volume size of about 32 liters. The spanwise region for the flow measurements is identical to the spanwise region with surface markers, as shown in Fig. 2. It is expected that at this position, the spanwise component of the flow velocity is negligible compared to the other two flow components which allows the flow field to be analyzed in a spanwise-averaged manner.
Three different stations of the CVV probe with respect to the airfoil model were used, as illustrated in a twodimensional view from the top in Fig. 4. The CVV position near the trailing edge is located about 50% further away from the airfoil than the other two CVV positions, to capture all surface markers over the full pitching motion and flap deflection range. The distances d 1 and d 2 indicated in Fig. 4 are approximately 0.4 m, whereas d 3 is approximately 0.6 m. These distances are sufficient to limit the effect of the presence of the CVV probe on the flow velocity around the model to less than 1% within a radius of 0.5 chord lengths around the airfoil . The larger distance d 3 partially impaired the quality of the flow field measurements from that CVV position, which exhibits larger values of random noise due to the low tomographic aperture angle of the CVV, as discussed in detail by Schneiders et al. (2018). The flow field measurements in this study cover only one side of the airfoil, in view of the symmetry of the pitching geometry.
PIV images are acquired at 700 Hz and processed with the LaVision DaVis 10 software. The procedure for obtaining the flow field measurements from the PIV images consists of applying a temporal high-pass filter to the images to remove light reflections (Sciacchitano and Scarano, 2014), performing a volume self-calibration (Wieneke, 2008(Wieneke, , 2018, generating a non-uniform optical transfer function (Schanz et al., 2013) and then applying the STB algorithm. The Lagrangian particle tracks that result from this procedure for the PIV image data from 45 recorded pitching motion cycles per CVV position are allocated to 100 temporal bins, each spanning 1% of the pitching motion period. Subsequently, the particle track data is ensemble-averaged spatially to a two-dimensional Cartesian grid with 5 mm grid spacing (i.e. 1.25% of the chord), based on a bin size of 20 mm × 20 mm with 75% overlap, hence producing phase-resolved flow field measurements in a spanwise-averaged sense. The ensemble averaging is performed using a top-hat filtering approach, as described by Agüera et al. (2016). Afterwards, a temporal sliding-average filter with a window size of 5% of the period is applied to each grid point to further reduce random errors in the unsteady flow field data. The measurement uncertainty of the flow fields is estimated from comparing the data in the overlapping regions of the adjacent measurement volumes, analogous to the procedure in Jux et al. (2018). The root-mean-square (RMS) of the discrepancies in the measured velocity magnitude are around 3% of the local velocity magnitude for the CVV positions 1 and 2, and around 7% for the CVV positions 2 and 3.
The PIV image data processing for the model position determination begins with removing the flow tracer information from the images by applying a temporal low pass filter to the acquired images. In the current study, a sliding minimum filter in time over a kernel of three images is used, as suggested by Mitrotta et al. (2019). In view of the relatively low pitching frequency and to reduce the processing time, only every seventh image from the PIV acquisition is used in the subsequent position determination procedure. The reduced sampling frequency of f s = 100 Hz results in a number of f s × T = 255 images per pitching cycle. After the image filtering, the STB algorithm is used to perform Lagrangian particle tracking of the fiducial markers on the airfoil surface.

Determination of the model position and flap deflection
The measurements of the marker positions are processed to determine the pitch and flap angles of the moving airfoil model. The processing method is illustrated in the following: first, measurements of the static airfoil are analyzed to determine the model position in laboratory coordinates and to illustrate the outlier filtering approach. In the second step, measurements of the pitching airfoil with fixed flap are used to illustrate the processing method for the moving airfoil. Results are discussed in detail for the pitching airfoil with actuated flap. This case is more challenging as the number of available marker position measurements per motion degree of freedom is smaller, compared to the other configurations (i.e., static airfoil and fixed flap).

Method
For the determination of the airfoil model position and flap deflection, the measured marker positions obtained in the coordinate system of the UR5 robot are fitted to the reference marker grid that was applied on the model. For the airfoil with fixed flap, this procedure involves finding the six degrees of freedom (DOFs) for the coordinate transformation from the laboratory reference frame to the airfoil reference frame. This is achieved by minimizing the overall distance of the marker measurements to the nearest reference marker position. The fitting procedure is described by the objective function f : where N is the number of tracked markers, ⃗ p i is the position vector of the ith marker track in laboratory coordinates, X , Y , and Z are the translations of the coordinate transformation, and ⃗ r is the position vector of the nearest reference marker in airfoil coordinates. ϕ, Θ and Ψ are the rotation angles of the coordinate transformation around the x, y, and z axes respectively, with rotation matrix ] . (2) The rotation Ψ around the vertical z-axis is directly related to the angle of incidence by α = −Ψ + α lab . Note that α is defined as nose-up (clockwise) positive. The angle of the laboratory coordinate system with respect to the inflow α lab is determined by analyzing position measurements of the static airfoil model at α = 0 • , as described in the following. The position corresponding to α = 0 • for the symmetric airfoil is determined by placing the airfoil model so that C ℓ (α) = 0, which, in this case, is verified by the measurements from the differential pressure transducers.
The result of the six DOFs fitting procedure is illustrated in Fig. 5 for marker tracking data from 100 images of the static airfoil model at α = 0 • . The fitting procedure is effective despite the presence of outliers which mainly occur due to reflections at the trailing edge and the wing-flap junction. The average number of marker tracks per time step including the outliers is 54.9, which is about 98% of the number of reference markers on the model. The low tomographic aperture angle of the CVV results in a relatively large measurement noise along the direction of the CVV viewing axis . The marker position measurements plotted as dots in Fig. 5 therefore appear as short line segments crossing the reference grid marker positions. The value for α lab that is found from this analysis is Ψ α=0 After an initial application of the fitting procedure, outliers are removed from the marker tracking data by setting a threshold for the distance to the nearest reference marker. When the optimization is performed for the second time for the static airfoil model after removing the outliers, the average number of marker tracks per time step is 48.2, which is about 86% of the reference markers. The remaining average distance from the track to the nearest marker is given by f min /N, which is 3.01 mm before and 1.21 mm after the outlier filtering.
For the analysis of the pitching airfoil, the fitting procedure is performed in two steps. In the first step, the fitting is performed for 500 images of the time series using the six DOFs optimization procedure using Eqs. (1) and (2) as described above for the static test case. The origin of the airfoil coordinate system is placed at the rotation point x r . The result of this first step of the analysis is shown in Fig. 6. It can be observed that all degrees of freedom in the optimization, except for the  rotation Ψ around vertical axis, are practically constant, with standard deviations of less than 1 mm for the translational DOFs and less than 0.1 • for the rotational DOFs. Therefore, the results obtained for the five remaining DOFs X , Y , Z , ϕ, and Θ are averaged and then prescribed in the second step of the fitting procedure, such that in this step only Ψ is determined. The second step is performed with two iterations per time step, with the outlier filtering applied after the first iteration.
For the analysis of the pitching airfoil model with actuated flap, the procedure is adjusted to allow for the additional degree of freedom corresponding to the flap deflection δ. This is achieved by deforming the reference marker grid in the optimization allowing for the variable flap deflection. Similar to the analysis of the pitching airfoil with fixed flap, the analysis is performed in two steps; after the determination of the five constant DOFs in the first step, the remaining two DOFs of the airfoil model motion Ψ and δ are determined in the second step.

Results
The results of the marker-based model position measurement for the pitching airfoil with actuated flap in terms of α and δ are shown in Fig. 7 in comparison with the corresponding reference data. The pitching motion that is imposed by the pitching mechanism is used as a reference for α, which is calculated from the measured rotation Ψ and the value that was found for α lab . The measured flap deflection δ is compared to the reference flap deflection angle which corresponds to the value of the control signal that is received by the servomotor which actuates the trailing edge flap. The reference flap deflection is calculated based on the kinematics of the flap hinge. The results from the marker-based position determination are sorted in phase for the comparison with the reference. Each data point in Fig. 7 corresponds to the α and δ obtained from the position measurement of the markers for a single time step.
The agreement between the determined angle of attack and the reference based on the motion kinematics imposed by the pitching mechanism is excellent, the RMS of the difference is 0.03 • . The RMS of the difference between the determined flap deflection angle and the reference based on the servomotor input signal is 0.44 • . Systematic differences to the reference and a higher level of random noise are expected in this case, because the servomotor that actuates the flap was operated closer to its power limit than the AC motor that drives the pitching mechanism. Furthermore, the actual rotation imparted by the servomotor was not measured but instead inferred from the actuator input signal. The largest systematic difference occurs at the beginning of the pitching cycle because the signal sent to the servomotor to generate the flap motion was updated for every pitching cycle at t/T = 0, causing a small delay in the servomotor motion here.
Apart from that, relatively large differences occur when the servomotor is moving slowly near the motion turning points,  where the flap motion was affected by friction between the airfoil model and the end plates, and when the servomotor is actuated at its maximum angular velocity, around t/T = 0.4 and t/T = 0.8. For the rest of the pitching cycle with 0.5 < t/T < 0.65 and 0.85 < t/T < 1, the agreement between the determined flap motion and the reference is similar to the agreement between the angle of attack α and its reference.
The obtained results for the airfoil model position and flap deflection over time are phase-averaged with a phase bin spanning 1% of the period, corresponding to the parameters used during the flow field processing. This correspondence facilitates the synchronous analysis of the phase-averaged model position and flow fields. The results thus obtained are shown for both test cases with fixed and actuated flap at four phase instants in Figs. 8 and 9, respectively.

Determination of the unsteady aerodynamic loads
Two methods are considered to determine the unsteady aerodynamic loads from the flow field measurements. The first method is based on a determination of the surface pressure distribution, while the second method is based on measurements of the circulation around the airfoil. The results are compared to the reference data from the differential pressure transducers installed in the airfoil model (see Appendix A).

Pressure-based load determination
For the determination of the lift and pitching moment on slender objects like airfoils, the primary interest is placed on the distribution of the aerodynamic load normal to the surface, as tangential forces are known to have little contribution The first term in Eq.
(3) directly relates the instantaneous value of the bound vortex strength distribution γ b (x, t) to the pressure difference and is therefore denoted as the quasi-steady pressure difference ∆p qs (x, t) = ρU ∞ γ b (x, t). The second term in Eq. (3) constitutes the unsteady contribution.
The aerodynamic load determination procedure begins with the calculation of the quasi-steady surface pressure using Bernoulli's equation: where V is the local velocity magnitude. The quasi-steady surface pressure cannot be determined directly from the velocity measurements due to the measurement errors and viscous effects that affect the near-wall region, and is instead determined using a linear extrapolation of the flow pressure obtained from applying Eq. (4) in the vicinity of the airfoil, similar to the approach of Ragni et al. (2009). The details of the implementation of the quasi-steady surface pressure determination are described in Appendix B.
Once the quasi-steady surface pressure C p,qs surf is determined for one side of the airfoil over the entire pitching motion period, the quasi-steady pressure difference on the airfoil ∆C p,qs is determined by subtracting measurements of C p,qs surf with a phase difference of ∆t/T = 0.5. This step is required as PIV measurements were conducted only on one side of the airfoil (see Section 2.3). The implications and underlying assumptions of this approach are discussed in Appendix A. When the quasi-steady pressure difference is known, the unsteady surface pressure difference can be calculated with Eq. (3). First, the vortex-sheet strength distribution along the chord is determined from the measurement of the quasi-steady pressure difference: With this result for γ b (x, t), the unsteady surface pressure difference follows from Eq. (3). From the pressure difference, the aerodynamic loads are determined by spatial integration. If the considered angles of attack α are small, so that cos(α) ≈ 1 and the influence of the tangential force on the lift is negligible, the unsteady lift per unit span L ′ (t) and the lift coefficient C ℓ (t) are obtained by chordwise integration of the pressure difference: where q ∞ = 1 2 ρU 2 ∞ is the dynamic pressure of the freestream. The nose-up pitching moment per unit span M ′ (t) and the moment coefficient C m (t) are similarly obtained by spatial integration of the pressure difference: In this study, the reference axis for the moment x ref is taken to correspond to the rotation axis of the pitching airfoil.

Circulation-based load determination
To circumvent the explicit determination of the surface pressure, the aerodynamic loads can be expressed in terms of the circulation Γ instead. The circulation is defined as the contour integral of velocity around the object: and is related to the integral of the enclosed vorticity according to Stokes' theorem (see e.g. Anderson Jr. (2011), chapter 2).
In the considered case, the circulation is due to the bound vortex strength on the airfoil γ b (x, t), as well as the vorticity shed in the wake, which is modeled by a vortex sheet of strength γ w (x, t). The bound vortex strength relates to the circulation bound to the airfoil as is the partial circulation on the airfoil between the leading edge and the chordwise position x, so that the total circulation bound to the airfoil is Γ c (t) = Γ b (c, t). The wake vorticity is produced as a result of the changes in bound circulation following Kelvin's theorem for the conservation of circulation. In the unsteady thin-airfoil theory it is assumed that vorticity is shed into the planar wake at the trailing edge and convected downstream at freestream velocity, so that U ∞ γ w (c, t) = − dΓc (t) dt (see Leishman (2006), chapter 8). It follows that the lift can be expressed in terms of the circulation when a chordwise integral is applied to the pressure difference according to Eq. (3): The first term in Eq. (9) is equivalent to the lift according to the Kutta-Joukowski theorem, which applies for steady flows, where as in the unsteady case, the bound circulation is time-dependent. The first term in Eq. (9) can therefore be considered as the quasi-steady lift contribution L ′ qs (t) = ρU ∞ Γ c (t). The second term in Eq. (9) relates to the flow acceleration effect, and may therefore be labeled as L ′ The unsteady lift L ′ (t) can hence be calculated from the bound circulation distribution, which is obtained from line integrals of the measured flow velocity. Similar as for the pressure-based approach, flow fields separated half a period are combined to form a complete flow field around the wing. The implementation of the circulation determination from the measured flow fields is described in Appendix C, where also the influence of the integration path on the obtained result for the circulation is discussed in detail.
The pitching moment can be calculated from the lift when the center of pressure x cp is known: with the center of pressure defined as Because Eq. (11) would require prior knowledge of the pressure distribution, it is convenient to express the pitching moment by introducing a separate center of pressure for each of the two lift terms in Eq. (9), where it was observed that where the center of pressure of the quasi-steady lift contribution is and the center of pressure of the lift due to flow acceleration effects is To calculate x cp,fa (t) with Eq. (14), knowledge of the temporal behavior of the circulation Γ b (x, t) is required. The computation of x cp,qs (t) with Eq. (13) would additionally require a further processing step to obtain the bound vortex strength distribution as the gradient of the bound circulation γ b (x, t) = ∂ ∂x Γ b (x, t). Therefore, Eq. (13) is further modified with an integration by parts: The pitching moment is hence determined from the circulation with Eqs. (12), (14) and (15).

Pressure-based load determination
Pitching airfoil with fixed flap. In Fig. 10, the surface pressure difference along the chord is shown for four phase instants.
The largest differences between the PIV-based pressure and the reference data appear in the airfoil nose region x/c < 0.1.
It is assumed that these larger differences are the result of insufficient spatial resolution of the flow field measurement to accurately capture the flow around the nose, which is characterized by relatively large spatial velocity gradients.
Downstream of x/c = 0.1, the PIV-based results are in good qualitative agreement with the reference data. The somewhat larger differences that appear around x/c = 0.2 and around x/c = 0.7 are considered as measurement artifacts due to the merging of the flow field data from the three different measurement volumes. The agreement of the PIV-based results with the reference data is improved compared to the quasi-steady formulation when the unsteady pressure formulation is adopted. The unsteady effects are the largest when the magnitude of the pitch rate is maximal, thus around t/T = 0 and t/T = 0.5. For example, at t/T = 0.05, the RMS of the difference to the reference data downstream of x/c = 0.1 is reduced by about 65% when the unsteady term is included. The aerodynamic loads as determined by a chordwise integration of the pressure difference are shown in Fig. 11. The qualitative behavior of the quasi-steady and unsteady loads over the period is similar, but the inclusion of the unsteady term introduces a phase shift of the signal. The phase shift is negative for the lift and positive for the moment, which can be explained by the fact that the unsteady term affects the load distribution mostly in the region downstream of the reference axis for the pitching moment computation that is located at x r /c = 0.3125, as it is observed in Fig. 10. This means that an increase in lift due to the inclusion of the unsteady term corresponds to a negative nose-up pitching moment and vice-versa, hence the opposing phase shift. For the lift and the moment, including the unsteady term produces a significant improvement in the agreement with the reference data in comparison to the quasi-steady loads. For the lift coefficient C ℓ , including the unsteady term results in a phase shift of ∆t/T = −2.2%, and the RMS error ϵ RMS , calculated from the difference to the reference, is reduced by 33%. For the moment around the pitch axis C m,r , the phase difference between the quasi-steady result and the unsteady result is ∆t/T = 5.2%, and ϵ RMS of the unsteady result is 54% lower than ϵ RMS of the quasi-steady result.
Pitching airfoil with actuated flap. The results for the chordwise distribution of the pressure difference at four phase instants for the pitching airfoil with actuated flap are shown in Fig. 12. The effect of the actuated flap is evident, as the pressure difference becomes negative downstream of the mid-chord, with a local extremum at approximately x/c = 0.75, which is the flap hinge point. Furthermore, the magnitude of the suction peak near the leading edge is reduced as a result of the change in the effective camber of the airfoil due to the flap deflection.
Only eight reference data points from the installed pressure transducers are available for comparison with the PIVbased results due to the particular setup of the flap actuation controller (see Section 2.2). Similar to the test case with fixed flap, the largest differences are observed for the largest magnitudes of ∆C p . The effect of the unsteady term, which is calculated as the temporal gradient of the chordwise integral of bound circulation (see Eq. (3)), is smaller than for the test case with fixed flap. That is because the controller is designed to keep the lift constant, therefore the flap motion tends to suppress the circulation variation and, hence, compensates an increase in the suction peak near the leading edge with a corresponding negative peak in the pressure difference further downstream, which means that the unsteady effects are only observed in the center region of the airfoil around x/c = 0.5. In this region, the agreement with the reference data is marginally improved when the unsteady pressure formulation is used, compared to the quasi-steady formulation.
The aerodynamic loads from the pressure-based approach for the pitching airfoil with actuated flap are shown in Fig. 13. No large differences are observed between the results when using the quasi-steady or the unsteady pressure formulation, which is in line with the observations for the pressure distributions themselves.

Circulation-based load determination
Pitching airfoil with fixed flap. The results for the circulation-based aerodynamic loads on the pitching airfoil with fixed flap are shown in Fig. 14. The quasi-steady result for the lift is in good agreement with the reference data, which is further improved when the unsteady formulation is used in which the result is shifted in phase by ∆t/T = −2.2%. For the pitching moment, the amplitude is underestimated by the circulation-based approach, so that only 75% of the amplitude of the reference data is reached, whereas for the lift 101% of the load amplitude is observed. This means that the center of pressure that was calculated with Eq. (13) is located too far downstream, closer to the pitch axis, that was used as the reference axis for the calculation of the moment. The phase shift in the moment signal of ∆t/T = 5.7%, that is produced by the unsteady formulation with respect to the quasi-steady result, is not large enough to reduce the remaining lag error with respect to the reference data below ϵ τ = −3.1% of the period.  Pitching airfoil with actuated flap. The aerodynamic loads for the pitching airfoil with actuated flap are shown in Fig. 15.
The quasi-steady and the unsteady formulation for the lift both provide results that are in very good agreement with the reference. While including the flow acceleration effects has an evident effect on the results, the agreement with the  reference is not significantly improved with this method. As discussed previously, the unsteady effects in this case are expected to be negligible, as the controller aims to maintain the circulation constant.
For the pitching moment, the amplitude of the circulation-based results is smaller than the amplitude of the reference, similar to the pitching airfoil with fixed flap, but less pronounced. In this case, the amplitude reaches 91% of the reference. In contrast to the observations for the pitching airfoil with fixed flap, no significant lag error is observed in the pitching moment results. The unsteady formulation is shifted in phase by ∆t/T = 0.4% with respect to the quasi-steady result, which reduces the RMS of the difference to the reference by 15%.

Comparison of the methods
The agreement of the results for the aerodynamic loads that were presented in the previous section with the respective reference obtained from the surface pressure transducer measurements is summarized in Table 2. The RMS error ϵ RMS is computed from the difference between measurement and reference. It is additionally given as a percentage of the maximum absolute value of the corresponding reference signal. The lag error ϵ τ is determined by cross-correlation with the reference signal and given both in seconds and as a percentage of the period T . The amplitude measurement error ϵ amp is determined by calculating the RMS error from the difference between the measurement and the reference after the measured signal has been shifted in phase according to ϵ τ .
The examination of the presented results yields the following observations: • The consideration of the flow acceleration effects in an unsteady formulation for the loads is highly effective for increasing the agreement with the reference data for the pitching airfoil with fixed flap, compared to the results from using a quasi-steady formulation. Similar improvements are observed for both load determination methods and originate primarily from a shift of the quasi-steady result in phase.
• For the pitching airfoil with actuated flap, the improvements due to the consideration of the flow acceleration effects are negligible. In this case, the flap actuation is designed to alleviate the loads, which also diminishes the unsteady effects on the load variation. Therefore, both quasi-steady and unsteady approaches yield accurate and comparable results, with RMS errors of less than 0.01 in lift and moment coefficients.
• For the pitching airfoil with fixed flap, the circulation-based approach yields more accurate results for the lift, when compared with the results from the pressure-based approach. In contrast to the pressure-based approach, the circulation-based approach does not use measurements near the airfoil surface, where measurement resolution problems are expected due to the presence of large velocity gradients and a turbulent boundary layer. The low measurement resolution was found to be most problematic for the pressure-based approach in the airfoil nose region, which is in accordance with the observations reported in the studies of Yuan and Olinger (2005) and Tagliabue et al. (2016), where a similar approach was used. Furthermore, the sensitivity of the circulation-based approach to individual flow velocity measurements is reduced by varying the integration path and averaging the result, so that the level of random error is decreased.
• For the pitching airfoil with fixed flap, the pressure-based approach yields more accurate results than the circulationbased approach for the evaluation of the pitching moment. The circulation-based approach underpredicts the pitching moment magnitude and exhibits a relatively large phase error in this case. It is noted that this lag is not a result of the selected approach to include the unsteady term, because the phase shift due to the consideration of flow acceleration effects is very similar in the pressure-based approach. Instead, a larger lag error is already observed in the quasi-steady circulation-based result for the pitching moment. The determination of the pitching moment with the circulation approach requires the integration of measurements near the surface so that the low measurement resolution can cause an accumulation of error here.
• For the pitching airfoil with actuated flap, the absolute measurement error is smaller for both methods, compared to the pitching airfoil with fixed flap. This can be explained with the reduced magnitude of the suction peak, which is the main cause of measurement resolution issues due to the large local velocity gradients near the leading edge of the airfoil. Compared to the test case with fixed flap, the presence of a negative peak in the pressure distribution and the associated measurement resolution problems cause a reduction of the error for the lift and an increase of the error for the moment, when considering the pressure-based approach. In comparison with the circulation-based approach, the pressure-based approach is slightly more accurate for the lift and the circulation-based approach is slightly more accurate for the pitching moment. However, the differences and overall errors are small enough so that both methods can be considered equivalent in this case.
Overall, it is established that both unsteady aerodynamic load determination approaches produce accurate results, in the two considered test cases with large and small lift variations, respectively. Neither of the two approaches appears as clearly superior to the other. Instead, it appears that the accuracy of the results for both approaches is compromised in different ways by deficiencies in the flow field measurement resolution. The pressure-based approach strongly depends on individual flow measurements near the airfoil surface, where the measurement resolution of this study is insufficient to accurately capture the flow behavior. Furthermore, the presence of the boundary layer on the airfoil increases the complexity of the implementation of the pressure-based method. In the circulation-based approach, flow measurements near the surface are not used in the determination of loads. Moreover, it is possible to reduce the sensitivity of the obtained loads to individual flow measurements by averaging the circulation result obtained with several integration paths. On the other hand, it is noted that the direct insights into the distribution of the loads from the pressure-based approach provide information on possible measurement issues that are not directly made available with the circulation-based approach.

Conclusions
An experimental approach was introduced to obtain the structural motion and aerodynamic loads on an unsteady airfoil using an integrated non-intrusive measurement approach. The integrated approach consists of simultaneously conducted optical measurements of flow and structure with a volumetric flow measurement device and the application of a Lagrangian particle tracking algorithm for the image data processing of both flow and structure.
The position of the unsteady airfoil was measured with the optical system by using structural markers on the experimental model. The image processing and marker tracking analysis was performed with the same particle tracking algorithm that is used for the analysis of the flow tracer particles. Next, the airfoil position in space was found by fitting the measurements of the marker positions to the structural marker grid that was painted on the model. Based on the experimental test case of a pitching airfoil with actuated flap, it can be concluded that the position determination method is suitable for reliably measuring the position and shape of deforming objects, as long as the deformation can be described in terms of a few degrees of freedom, as it is the case for rigid body motion with the two degrees of freedom in this study.
The unsteady aerodynamic loads in terms of the lift and the pitching moment are determined from the flow field measurement in this study using two novel implementations of an aerodynamic loads determination method that makes use of existing unsteady aerodynamics theory. The pressure-based and the circulation-based load determination approaches were both adapted with formulations from unsteady potential flow and thin airfoil theory to extract the unsteady aerodynamic loads. It was shown that adopting the unsteady formulations reduces the difference between the PIV-based load measurements and reference pressure data by more than 50% for the pitching airfoil without flap actuation, compared to results using quasi-steady formulations for the aerodynamic loads. It was also shown that both methods correctly predict the reduction of unsteady effects when the flap is actuated. While it was observed that the two different methods are affected differently by measurement resolution issues, no method was identified as providing superior results compared to the other. However, it is noted that the agreement of the results from both methods with the reference is very good, with typical errors of less than 0.01 in lift and moment coefficient.
The integrated approach for the determination of the unsteady fluid dynamic loads and the structural motion that was introduced in this study is useful for large-scale fluid-structure interaction experiments because the use of only a single measurement and data processing system in combination with the minimal model instrumentation requirements represent a significant reduction of the setup complexity, data processing effort and costs for these experiments. Furthermore, the two novel approaches for the aerodynamic loads determination that were introduced are connected to less data processing effort and pose lower measurement requirements on the flow field measurements than existing alternative methods that are based on the momentum equation in integral form, which require the evaluation of several more terms with multiple temporal derivatives and gradient operations (Rival and van Oudheusden, 2017). It should be noted that these methods can also be applied when the flow field data is obtained with a different measurement approach, such as a classical two-dimensional PIV setup, when the structural motion is known or measured with a separate system.
Although the analysis of the pitching airfoil with actuated flap was performed in a two-dimensional sense in this study, the measurement data was obtained in three-dimensional space. Therefore, the measurement approach is not inherently restricted to two-dimensional problems and a wide range of potential applications exists. However, it should be noted that with the sequential measurement procedure of the robotic PIV system used in the present implementation, applications are limited to periodic or repeatable phenomena. This limitation can be overcome by using an optical setup with multiple cameras instead, allowing to capture the entire flow field simultaneously. Furthermore, considering that the aerodynamic loads determination methods in this study were derived based on the assumption of unsteady potential flow, further research is necessary to assess the performance and possible necessary adaptations for the application of these methods in flow situations where the viscous effects are more dominant than in this study, such as in the case of separated flow over the airfoil.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments
The authors gratefully acknowledge the help of Jorge Fernández Barrio in preparing and conducting the experiment. This work has been carried out in the context of the HOMER (Holistic Optical Metrology for Aero-Elastic Research) project that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 769237.

Electronic supplementary material
The data that is presented in the figures of this article are available for download at: https://doi.org/10.4121/13810721.

Appendix A. Reference data from installed pressure transducer measurements
The pressure transducer measurements are acquired over 45 complete pitching cycles at a sampling rate of f s = 32 Hz.
However, because the PIV measurements were conducted only on one side of the airfoil model, a comparison between the loads derived from the flow field measurements and the differential pressure transducer data cannot be done directly.
Therefore, PIV measurements with a phase difference of ∆t/T = 0.5 have been combined, so as to reconstruct a flow field measurement that covers both sides of the airfoil in the optical load determination procedures. To produce consistent reference data for the PIV-based loads, the pressure transducer measurements from two phase instants separated by ∆t/T = 0.5, which correspond to one PIV-based load result, are combined to generate one reference data point. Based on the assumption that the aerodynamic load scales linearly with the angle of attack α for small values of α, it can be assumed that the PIV-based load, that is obtained from the two separate flow field measurements, will be equivalent to the mean value of the aerodynamic load at the two phase instants where the flow field was measured.
The procedure to generate the reference data for the lift coefficient C ℓ from the differential pressure transducer measurements is illustrated in The reference data is used for the comparison with the results of the PIV-based lift determination methods. Following a similar procedure as for the lift, sets of reference data are also determined for the pitching moment and surface pressure difference.

Appendix B. Determination of the quasi-steady surface pressure
The method to determine the quasi-steady surface pressure on the airfoil from flow velocity measurements in the vicinity of the airfoil is described in the following. In view of the presence of the boundary layer on the model, Bernoulli's principle (Eq. (4)) is unsuited to determine the pressure close to the surface. Therefore, the surface pressure is determined in this study from extrapolation of pressure data obtained at some distance from the surface, as illustrated in Fig. B.18.
For the chordwise positions with x/c ≤ 0.1, it is assumed that the boundary layer region is thin, compared to the spatial bin size that is used during the ensemble averaging, so that the flow measurements are not significantly affected by the presence of the boundary layer and Eq. (4) can be applied even close to the airfoil. Therefore, the surface pressure is determined in three steps that are applied for each vertical line of the Cartesian grid of the flow field measurement where x/c < 0.1: 1. Application of Eq. (4) to the flow measurement at the first grid point outside of the airfoil y 1 , to determine C p,qs 1 at a distance of ∆y 1 from the airfoil surface.
2. Application of Eq. (4) to the next following grid point away from the surface in y-direction, y 2 with C p,qs 2 , to approximate the pressure gradient ∂C p,qs 1 ∂y with a backward finite difference as ∂C p,qs 1 ∂y ≈ C p,qs 1 − C p,qs 2 y 1 − y 2 .
3. Extrapolation of the pressure on the surface of the airfoil C p,qs surf from the nearest grid point with C p,qs surf = C p,qs 1 − ∂C p,qs 1 ∂y · ∆y 1 .    For the determination of C p,qs surf downstream of x/c = 0.1, it is assumed that the application of Eq. (4) close to the airfoil surface yields erroneous results due to the presence of a thickened boundary layer. It is therefore necessary to exclude these data points from the analysis. The development of the boundary layer thickness on the airfoil can be roughly estimated using an empirical expression for the turbulent boundary layer on a flat plate with zero pressure gradient δ BL (x) (see White (2006) For the flow field measurements from the most downstream position of the CVV, where x/c > 0.65, more measurement noise is observed. To reduce the sensitivity of the surface pressure result to individual data points, the pressure gradient is not determined from a finite difference in this region, but with a linear regression of multiple data points instead, in this case over 10 data points. This approach also diminishes the modeling errors associated with the application of Eq. (B.1) to a region of non-zero pressure gradient, which would affect the measurements at y 1 and y 2 . The result of the linear regression is used to extrapolate the surface pressure in the chordwise region with x/c > 0.65, as illustrated in Fig. B.18.

Appendix C. Determination of the bound circulation
The distribution of the bound circulation on the airfoil is determined in this study by performing line integrals of the measured flow velocity, as described in the following. The procedure begins with a systematic analysis of the effect of the position of the integration path for determining the circulation.
Three rectangular integration paths are shown in Fig. C.19a, and the corresponding result from the line integration over the period is shown in Fig. C.19b. Note that the shown integration paths do not form a closed integration path, as flow field data is only available for one side, and therefore the obtained result from the corresponding partial line integral is labeled as pseudo-circulation Γ ′ . As a result of the open integration paths, it is necessary to subtract the mean value,Γ ′ , from the result of the line integration, in order to facilitate a direct comparison of the results obtained from the different integration paths. Path I is placed at a distance of 0.05 chord lengths from the airfoil surface, which means that the upstream segment of the integration path is located at x/c = −0.05, the downstream segment at x/c = 1.05 and the streamwise segment is located at y/c = −0.14, to additionally account for the thickness of the airfoil. Integration paths II and III are placed at a distance of 0.15 and 0.25 chord lengths from the airfoil, respectively.
Overall, the sensitivity of the pseudo-circulation to the position of the integration path observed in Fig. C.19 is small. A phase lag is observed with increasing distance to the airfoil, such that the result obtained with path III is shifted in phase by about 1% of the period with respect to the result obtained with path I. Furthermore, it appears that the amplitude of the pseudo-circulationΓ ′ is reduced with increasing distance of the integration path to the airfoil so thatΓ ′ is approximately 3.5% smaller for path III when compared to the result obtained with path I.
In the next step, the effect of varying the position of the three individual line segments of the integration path is investigated. In the following analysis, the distance to the airfoil for one of the three line segments is varied, while the other two distances remain constant at 0.15 chord lengths to the airfoil. The results of the line integrals with the varying integration paths are expressed in terms of the amplitudeΓ ′ and the phase lag τ , which is determined by cross-correlation and therefore supports the assumption of the wake vortex sheet being shed at the trailing edge and convected downstream at U ∞ . In unsteady potential flow, the circulation is independent of the positions of the upstream and streamwise segments