Dense interpolations of LPT data in the presence of generic solid objects

Three-dimensional Lagrangian particle tracking measurements with helium filled soap bubbles provide quantitative flow visualizations in large measurement volumes up to the cubic metre scale. However, the instantaneously available fluid information density is severely restricted by the finite spatial resolution of the measurements. Therefore, the use of experimental data assimilation approaches are utilized to exploit the temporal information of the flow measurements, along with the governing equations of the fluid motion, to increase the measurement spatial resolution. Nevertheless, only in the last years, attempts to apply data assimilation methods to enhance the Lagrangian particle tracking (LPT) resolution in proximity of solid boundaries have been performed. Thus, in order to handle generic solid body intrusions within the densely interpolated LPT data, two different approaches based on the computational fluid-structure interaction frameworks are proposed. The introduced variants of the state of the art physics-driven data assimilation methods are assessed with a high fidelity numerical test case of flow over periodic hills. The accuracy superiority of the flow field reconstructions with the proposed approaches are denoted especially in close proximity of the interaction surface. An experimental application of the introduced methods is demonstrated to compute the pressure distribution over an unsteadily moving elastic membrane surface, revealing the time-resolved interaction between the flow structures and the membrane deformations.


Introduction
Tomographic particle image velocimetry (Tomo-PIV) (Elsinga et al 2006) allows acquisition of flow information * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. over a three dimensional measurement domain. However, Tomo-PIV measurements are accompanied with stringent constrains on the maximum achievable measurement volume due to the dispersed intensity of light illumination, light scattering performance of tracer particles and resolution characteristics of recorded images. Consequently, there exists an inverse relationship between the scattered and captured light intensity magnitudes and the depth of the measurement volume which puts strict limitations on the maximum allowable measurement volume for Tomo-PIV applications (Tokarev et al 2013). Accordingly, the first applications of Tomo-PIV were performed within volumes of of 13 cm 3 . Further applications by Schröder et al (2009), Violato et al (2011)and Ghaemi and Scarano (2011) were only able to reach 20 cm 3 while Fukuchi (2012) documented a measurement volume of 16 × 22 × 8 cm 3 for which the acceptable signal-to-noise ratio distribution was limited to a depth of 5 cm. Then, the maximum measurement volume was attained by Henningsson et al (2015) with a size of 20 × 24 × 5 cm 3 .
Hence, in order to enable large scale applications of Tomo-PIV, the use of neutrally buoyant helium filled soap bubbles (HFSB) as tracers for PIV experiments was proposed (Scarano et al 2015, Caridi et al 2016. Their controllable larger size and density, provided a great opportunity of achieving suitable flow tracing qualities and superior light scattering characteristics. Nevertheless, the limited production and seeding rate of the HFSB restrict the concentration of tracer particles within the measurement domain which in turn alleviates the instantaneously available spatial resolution. Thus, the reliability and accuracy of volumetric reconstruction might be compromised. Considering the spatial scarcity of instantaneously available HFSB particles, employing a Lagrangian particle tracking (LPT) approach enables a greater level of accuracy in terms of determining particle location, velocity and acceleration (Maas et al 1993). Accordingly, the 'Shake-The-Box' algorithm introduced by Schanz et al (2016), which performs an image-matching-based correction of particle distribution predictions in consecutive time-steps, yielded considerable improvements in accuracy and resolution of LPT information in comparison to the standard Tomo-PIV approaches. However, the data obtained from LPT is scattered over the measurement domain. Although straightforward interpolation of scattered data might attenuate the problem of reduced spatial resolution, it does not allow to resolve length scales smaller than the distance between two adjacent particles. Therefore, appropriate data assimilation methods are required for increasing the spatial resolution of the available experimental data which would provide greater accuracy and physical coherence of flow field interpolations.
The state-of-the-art data assimilation frameworks can be identified in two main groups depending on the basis of their algorithms. Based on statistical formulations, the Adaptive Gaussian Windowing (AGW) proposed by Agui and Jimenez (1987) captures the scattered information in a measurement volume on a predefined grid using a Gaussian weight averaging over the windows of Euclidean distances defined with respect to the grid location of interest (figure 1, top-left). Furthermore, the computed Gaussian average is scaled with a standard deviation corresponding to the desired spatial resolution. A more advanced algorithm based on Gaussian windowing is developed by Casa and Krueger (2013), where a Gaussian radial basis function (RBF) derived by Weiler et al (2005) is utilized with an iterative optimization procedure to obtain the best fit of data interpolation on a grid structure with respect to the original scattered dataset.
In order to incorporate the physical description of fluid dynamics within the data assimilation procedure, Gesemann et al (2016) introduced the use of 3D cubic B-splines for increasing signal-to-noise ratio of particle information reconstruction over particle velocity and accelerations. The method includes two main steps referred as TrackFit and FlowFit. TrackFit initiates the reconstruction algorithm with a noise reduction of the particle intensity signals similar to Wiener-/Kalman filter (Tenoudji 2016). This is performed by calculating B-splines over particle tracks which also provides the first and second order gradients of particle velocities since the least order of the spline is considered to be 2 using three particles of different time instants. Then, 3D B-splines are calculated for velocity and acceleration (or pressure) information by FlowFit where an iterative minimization of a cost function is performed enforcing various physical constraints such as the divergence freedom of velocity and velocity gradients accounting for conservation of mass and momentum (Navier-Stokes equation).
The ability of using the temporal information to enhance the measurement spatial resolution and accuracy was exploited by Schneiders et al (2015), who developed the VIC+ method based on vortex-in-cell (VIC) model presented by Christiansen (1973). The VIC+ method provides an alternative algorithm that allows the spatial resolution characteristics of raw LPT data to be increased while enabling the representation of the fluid behavior on a global sense over the prescribed computational domain (figure 1, top-right). Moreover, the vorticity transport equation enables the computation of the local time derivative component for the Lagrangian acceleration term only by means of velocity and vorticity fields via eliminating the contribution of pressure gradients. Hence, VIC+ also provides an opportunity to increase the resolution of material derivative information as its optimization procedure is based on maximizing the proximity of velocity and acceleration vectors to the measured values .
Furthermore, in order to take advantage of the recorded temporal history of particle tracks, the optimization procedure of VIC+ based on the velocity and vorticity distributions is expanded over Lagrangian particle trajectories reconstructed using multiple time instants . The method introduced is referred as VIC++ where the main objective is set to utilize the full potential of the temporal range of two dimensional measured data to increase the density of velocity information over the specified measurement domain in consistence with the temporal evolution of particle tracks. Further on, an extended version of VIC++ is introduced as time-segment assimilation (TSA) where the optimization problem of reconstructing a dense field of velocity information is performed on three dimensional measurement domains (Schneiders andScarano 2018, Gonzalez Saiz et al 2019). Based on a minimum of four time consecutive velocity distributions, TSA performs an iterative optimization procedure using the velocity information of a particle along with its propagation over a chosen time segment.
Although VIC+ and its variants are proven to successfully increase the spatial resolution of LPT data on uniformly structured computational grids, their applicability is restricted to purely fluidic domains. However, most practical engineering applications involve dynamic fluid-structure interactions which contain non-uniform solid boundaries that might deform unsteadily. Thus, computation of fluid information in presence of curved walls or solid objects (stationary or moving) introduces significant drawbacks in terms of numerical accuracy and appropriate surface definition for implementation of adequate boundary conditions. Without proper definition of solid boundaries, accurate determination of fluid properties along the structural surfaces becomes significantly problematic as boundary conditions and the information transfer between fluid and solid domains are strongly dependent on the interface description (Dowell 2004). More recently proposed methods for managing non-uniform solid intrusions, separate the computational domain in to multiple regions. In these approaches, the incapacity of available data assimilation algorithms to capture near surface information yields loss of valuable information (Gonzalez Saiz et al 2021). Although these methods provide accurate results in comparison with the extrapolation techniques, the improvements are associated with increased levels of complexity and computational cost. Therefore the current study introduces two approaches for providing VIC+ algorithm with the capability of performing data assimilation in presence of solid boundaries by implementing the well known computational fluid-structure interaction (FSI) frameworks of Arbitrary Lagrangian-Eulerian (Noh 1963) and immersed boundary treatment methods (Peskin 1982).

Non-uniform boundary treatments for VIC+
Two main approaches, named ALE-VIC+ and ImVIC+, are proposed to handle the solid boundary effects for VIC+ applications. Accordingly, section 2.1 gives brief background information on the boundary treatment approaches utilized in computational fluid dynamics simulations (CFD) which form the basis for the proposed methods. Then, sections 2.2 and 2.3 discuss the numerical implementation those approaches within the frameworks of ALE-VIC+ and ImVIC+ respectively.

Treatment of non-uniform boundary deformations in CFD
The ALE-VIC+ approach utilizes the Arbitrary Lagrangian-Eulerian developed for computational FSI simulations as outlined in section 2.1.1, accompanied with boundary fitted coordinate systems and mesh adaption procedures (figure 1, bottom-left). In ImVIC+ method, instead, the immersed boundary treatment is utilized to satisfy the appropriate boundary conditions while preserving uniform mesh formations which is discussed in section 2.1.2 (figure 1, bottom-right).

Arbitrary Lagrangian-Eulerian method.
In order to avoid the individual shortcomings of Eulerian and Lagrangian perspectives while leveraging from their respective advantages, a technique referred as the Arbitrary Lagrangian-Eulerian (ALE) method was introduced by Noh (1963), where both approaches are utilized in a coupled manner to handle appropriate fluidic and structural domains simultaneously. As the necessity of a coupled method emerges from the motion of boundaries enclosing the fluidic domain, the corresponding mesh formations are required to be modified accordingly. Hence, a natural approach can be considered as completely regenerating the grid structure at each time step. However, this process generally requires considerable user interaction and immense computational resources (Luke et al 2012).
Instead, grid deformation algorithms provide a valuable solution categorized under two main groups of physical analogy (Löhner and Yang 1996, Farhat et al 1998, Helenbrook 2003 and interpolation based schemes (Liu et al 2006, Jones and Samareh-Abolhassani 2012, Wang and Przekwas 2012. Utilizing these various methodologies, there exists multiple applications of the ALE method for vortex simulations based on the VIC model, for cost-effective high-fidelity numerical simulations of FSI problems. In this regard, Cottet and Poncet (2004) performed conformal mapping of the fluidic domain for three dimensional direct numerical simulations of wall bounded flows where the grid structure in close proximity of a sphere was fitted to the surface of the solid object for application of a hybrid particle-mesh vortex method. Furthermore, Kudela and Kozlowski (2009) employed a boundary fitted coordinate system for flow simulations around arbitrary shaped objects using the VIC framework during which fourth order interpolation kernels used by Cottet and Koumoutsakos (2000) are modified for particle-mesh switching of vorticity distributions in the near wall regions.

Immersed boundary treatment in CFD.
The VIC framework introduced by Christiansen (1973) allows the fast Fourier transform (FFT) based Poisson solvers to be employed over a predefined computational grid (Wu and JaJa 2013) allowing higher fidelity and resolution capabilities, while preserving the computational efficiency. As the fast Poisson solvers are utilized to characterize the rotational component of velocity vectors, physical intrusions within the flow field are required to be handled employing additional velocity or forcing terms. Therefore, the need for boundary fitted coordinate systems and introducing transformation operations between computational and physical coordinate systems is obliterated. Accordingly, Walther and Morgenthal (2002) and Cottet and Poncet (2004) implemented integral boundary equations on the VIC method in order to impose no-through and no-slip boundary conditions by means of surface singularities. However, defining no-through boundary condition over the penetrating velocity components results in an integral boundary equation which corresponds to an ill-posed Fredholm integral equation of first kind. Instead, the vortex sheet strengths are determined via the tangential velocity components yielding a Fredholm integral equation of second kind (Beale and Greengard 1994). Furthermore, in order to obtain a unique solution for multiply connected regions, Kelvins theorem of circulation conservation or Kutta condition is introduced as an additional constraint (Morgenthal and Walther 2007).

Arbitrary Eulerian-Lagrangian approach for VIC+
The ALE-VIC+ (figure 2, black) comprises of two main steps in which a computational grid is generated according to the surface information (section 2.2.1) and the flow governing equations within the VIC framework are solved over the generated boundary fitted coordinates (section 2.2.2). Then the problem closure is described by the definition of a costfunction for the optimization procedure (section 2.2.3).

Boundary fitted coordinate system adaptation.
Depending on the availability of the surface information, two different approaches can be followed to reconstruct the corresponding boundary shapes. If the immersed boundary shape is known a priori, it is directly provided as an input to the grid generation schemes along with the appropriate volume dimensions for the boundary fitted coordinate system to be created. In case where the unsteady motion of the solid boundary is captured by means of optical measurement methods, the captured structural motion is instantaneously provided to the grid generation algorithms. After an initial boundary fitted grid structure is generated, the mesh is deformed according to the boundary motion so that the exact interaction between the fluidic and solid domains can be expressed in a time-resolved manner. In order to achieve a continuous boundary conformation, the choice of deformation scheme is determined as RBF based mesh deformations (de Boer et al 2007) due to their proven mesh deformation accuracy and cost-efficiency (Smith et al 2000, Beckert andWendland 2001).
The formulations employed for RBF based mesh deformation method utilized are described by de Boer et al (2007) where a comparative study over many variants of RBFs was conducted using two different mesh quality metrics based on mesh size preservation and mesh skewness. Among 14 different RBFs with compact and global supports, thin plate spline and CP C 2 RBFs stood out by keeping the size alterations and skewness to a minimum whilst not compromising mesh adaptation efficiency. According to the presented results, CP C 2 RBF is utilized to update the computational grid due to a superior performance of mesh deformations especially in close proximity of the immersed boundary surface.

VIC+ on boundary fitted coordinates.
The flow governing equations on the updated conformal grid are solved, by performing a one-to-one mapping between the physical and a computational coordinate systems. The mapping relates to the transformation variables necessary for accurate description of mathematical operators constructing the link between the two grids. The physical coordinate system is selected as the Cartesian one since the flow field properties are obtained from the PTV/LPT measurements are already defined on a Cartesian coordinate system. Hence, the computational grid locations in any mesh form can be expressed as functions of the physical coordinates, where subscripts P and C, refer to the physical and computational coordinates respectively. Therefore, the flow governing equations within the VIC+ method are described on the computational coordinates by transforming the vector variables utilizing the chain rule and a transformation matrix that defines the link between the two coordinate systems. . The LPT measurements for the flow is provided (green box) and an initial velocity estimate made (black-2 & grey-1) to be introduced into VIC+ iterative procedure (blue boxes). In case of ImVIC+, the velocity fields are corrected with immersed boundary treatment (grey-3).
The two governing equations of the VIC framework described as the velocity(u)-vorticity(ω) formulation in equation (3) and the vorticity transport equation in equation (4), The viscous terms are neglected in equation 4. as the inviscid Navier-Stokes equations are considered to be valid due to fact that on large scale applications with HFSB tracers, particle concentrations near the walls are significantly reduced and the tracers are observed to be too large for entering the viscous sublayers (Faleiros et al 2018). Vector operators (∇ P , ∇ P × and ∇ 2 P referring to the gradient, curl and Laplacian respectively) can be expressed on the computational coordinates through equations (5)- (7), in which J T is the Jacobian of the transformation matrix and I is the identity matrix. Hence the velocity-vorticity relation in equation (3) where C is the cofactor of the transformation matrix. Thus, the corresponding the vorticity transport equation described by equation (4) is derived as follows.

Cost-function computation and adjoint gradient optimization.
As the orientation of velocity and acceleration vectors are preserved on the Cartesian descriptions, the resultant flow properties of velocity and material acceleration values at the original particle locations are used for calculating the error between the dense flow field interpolation and the measurement data. Following the exact procedure of adjoint gradient computation introduced for VIC+ with, additional modifications are implemented to take into account the link between the physical and computational grid structures. Hence, the adjoint gradient for the proceeding steps of the optimization procedure is computed and the optimization procedure is performed until a specified convergence criteria of J F /J I = 10 −3 of the initial value of the cost function is achieved.

Immersed boundary treatment for VIC+
Implementation of the immersed boundary treatment for VIC+ method is based on the theory of vector decomposition provided via the Helmholtz theorem (figure 2, gray). First, surface singularity elements are distributed over the solid boundary surface to satisfy the appropriate boundary conditions throughout the data assimilation procedure (section 2.3.1). The VIC method allows rotational component of the flow behavior to be characterized by means of distributed vortices and their respective velocity induction over a prescribed domain (section 2.3.2). Hence, the potential flow component of the velocity field due to the presence of a non-uniform surface intrusion is discretized by superimposing an additional velocity field due to the introduced surface singularity distribution (section 2.3.3). Then the problem closure is provided by the definition of a cost-function for the optimization procedure (section 2.3.4).

Solid boundary characterization and singularity attachment.
The numerical implementation of boundary integral equations for adequate description of boundary conditions over the interface surfaces is performed using the panel method introduced by Hess and Smith (1967). Considering that the structural information is obtained by means of either one of the aforementioned approaches described for ALE-VIC+, the solid boundary surface is characterized by means of quadrilateral panels each equipped with singularity elements of sources and doublets to introduce a scalar potential influence of the physical intrusion (figure 3). With singularity elements attached, the potential inductions of the quadrilateral panels are computed over the complete computational domain employing the formulations provided by Katz and Plotkin (2001).

Vortex-in-cell method application for rotational
component calculation. The initial estimation for the velocity field information is provided in terms of a freestream velocity value. Hence, the computational grid locations are equipped with pointwise vortices and the vorticity distribution over computational domain is utilized to compute the velocity distribution using the velocity-vorticity relation for the rotational velocity vector component in equation (10). (10)

Determination of potential flow component.
The calculated velocity field penetrates the physical boundaries. Thus, to satisfy the no-through boundary condition over the solid boundaries, the normal velocity components to each panel are equated to the relative source strengths (σ i ) of the corresponding panels. Hence, the strengths of the dipoles (µ i ) are reduced for ensuring numerical uniqueness (Tarafder et al 2010).
Furthermore, since the self induced scalar potential should vanish at collocation points over a surface of singularities (Lewis 1991), the dipole strengths are calculated by constructing a linear system of equations imposing that total sum of potential induction at the central locations of the quadrilateral panels equals to zero.
In order to establish a unique solution, the corresponding vortex rings (dipoles) are imposed to comply with the conversation of circulation while their strengths are being determined. The resultant overdetermined linear system is solved employing a least squares method (Soifer 2013).
Moreover, the scalar potential field induced by the surface singularities is differentiated in three dimensions to obtain the velocity vector components, which are then superimposed to calculate the resultant velocity field distributions.
Then, the corresponding velocity and vorticity fields are utilized to calculate the material derivative distributions over the computational domain in accordance to the inviscid Navier-Stokes formulation.

Cost function determination and optimization
procedure. The cost function for the optimization procedure is evaluated based on the errors of velocity and material acceleration values at the original particle track locations in comparison to the measured data. For the optimization procedure, it is assumed, due to the dependence of the potential flow component to the rotational elements, that the errors directly relate to the vortex-in-cell base which is dictated by the vorticity distributions. Hence, for each step of the optimization procedure, the gradient is calculated in terms of the vorticity strengths and the corresponding potential flow component is calculated to correct the velocity field distributions in order to satisfy the physical boundary condition of no penetration through the solid surfaces.

Test case and data processing
Validation studies of the proposed ALE-VIC+ method are performed with a direct numerical simulation (DNS) of flow over periodic hills. The simulations are conducted with a computational domain composed of two consecutive hill forms connected over a region of 144 h 3 , where h represents the nondimensional heights of the hills. Hence, the length dimensions of the domain in Cartesian coordinates are provided over normalized values with h. The numerical simulations are performed with periodic boundary conditions for inflow and outflow boundaries while the non-uniform surface contours of the hills are treated by means of an immersed boundary method to account for their influence over the fluidic domain imposing the no-slip boundary condition over the intruding hill surface (Chen et al 2014). The inflow conditions are initialized with a non-dimensional velocity distribution of u = 1 which corresponds to a hill height based Reynolds number of Re h = 10 595 as the non-dimensional kinematic viscosity is prescribed to be ν = 9.45 × 10 −5 .
In order to provide an accurate representation of scattered particle track information, the DNS results are randomly downsampled with particle concentrations of C = 25 par/h 3 , C = 125 par/h 3 and C = 250 par/h 3 . The corresponding particle image densities over the projected volumes obtained as N P = 0.01 ppp, N P = 0.05 ppp and N P = 0.1 ppp allowing a sensitivity investigation of the proposed data assimilation approaches to available particle concentrations. Simulation of LPT data structure is achieved via a pseudo-particle tracking approach applied over the instantaneously extracted flow field data. The procedure is initiated with random downsampling of high density flow information of the DNS data field in accordance to the prescribed particle concentrations for a selected time instant. Then a Runge-Kutta 4 time integration procedure is used to compute the particle propagation in time (Zheng and Zhang 2017). In order to minimize the truncation error propagation due to the numerical approximations, the integration procedure is applied in both forward and backward directions in time. For accurate extraction of velocity and acceleration values, the tracks are selected to be composed of seven time-steps. Then a Gaussian noise of 0.2 voxels in all three dimensions is artificially introduced to particle locations along the tracks to account for measurement and reconstruction errors (Sciacchitano et al 2022). Finally, the corresponding particle tracks are regularized with polynomials of order 2 to compute the velocity and material acceleration.

Results
To be able to demonstrate the improvements obtained with the proposed approaches over the state of the art of data assimilation, the standard VIC+ method is also employed to perform the dense interpolation of the velocity and material accelerations. However, since the base algorithm of VIC+ is not capable of handling the non-uniform boundary, the grid locations corresponding to the solid domain are artificially modified to have zero velocity and material acceleration values to prevent numerical divergence of the optimization procedure. Moreover, the spatial resolutions for the dense flow field reconstruction procedures are selected to have comparable grid spacing with the DNS simulations (∆h = 2 × 10 −3 h), to achieve better comparison of the reconstructed flow properties with the reference flow field. The corresponding computational grids for VIC+ and ImVIC+ (uniform structured grid with the immersed boundary surface) and for ALE-VIC+ (curvilinear grid boundary-fitted on the hill surface) are illustrated in figure 4 (left) and figure 4 (right) respectively.
To start with the pure LPT approach where the AGW is used to capture flow properties in close surroundings of each grid locations, the results failed to provide almost any relevant fluidic information due to the lack of particles. Hence, the corresponding need for an interpolation approach to reconstruct a coherent flow field description is addressed by linear interpolation, VIC+, ALE-VIC+ and ImVIC+ approaches. The initial comparisons performed qualitatively over the streamwise velocity magnitude distributions reveal similar results in terms of the overall coherence characteristics for both linear interpolation and VIC+ based approaches. Both methods identify the major flow structures of high and low velocity regions in relation to the reference simulations (figure 5). Accordingly, the reconstructed streamwise velocity distributions using both approaches reveal accelerated flow behavior over the hill form and separation regions with recirculatory flow characteristics downstream of the hill. Nevertheless, the detailed flow structures of local high and low velocity magnitude variations are depicted with increased accuracy using the VIC+ based methods as the separation effects are captured with greater agreement to the reference solution. Furthermore, the superiority of VIC+ variants becomes even more apparent when the comparisons are performed over the particle concentrations N P ⩾0.05 ppp (figure 5, middle & right). As the achievable resolution limits using data assimilation methods are strictly dominated by the interparticle distance of the available LPT data , low particle concentrations do not favor the excessive spatial resolution where the flow information beyond ∆h min = 0.25C −1/3 is optimized without a proper cost function definition. Thus, increasing the particle concentrations while keeping the dense interpolation grid constant  allows reconstruction of finer structures which demonstrates the ability of VIC+ based approaches to resolve the fluidic behavior with greater detail as well as preventing overestimation and underestimation errors by means of appropriate physical definitions.
Major differences between the original VIC+ and the proposed variants are observed especially in close proximity of the hill surface. The peak fluctuation magnitudes are resolved with a greater agreement to the reference simulations by employing both ALE-VIC+ and ImVIC+; conversely, the independence of VIC+ from the relevant surface structure caused overshoots and undershoots closer to the interface, which also propagated towards the regions further away from the surface ( figure 5, 4th row). Although the dominant effect of the boundary treatments for VIC+ is observed over the streamwise velocity, mitigation of velocity reconstruction errors also denotes the modifications influencing both normal and spanwise components (figure 6). The differences in the error distributions tend to decrease significantly in regions away from the hill surface, because the influence of singularity elements decay with the Euclidean distance. However, the corrections due to immersed boundary treatment implemented by ImVIC+ have no influence on the spanwise velocity components. Considering the orientation of the surface elements for the hill form and the corresponding normal vector directions, there are two velocity components (streamwise and normal) of the flow motion that induce penetration of fluid particles through the hill surface. Accordingly, the application ImVIC+ provides modifications to the flow properties within those directions to prevent unphysical flow penetration through the surface by ensuring the satisfaction of no-through boundary condition. Since the spanwise velocity components do not affect the satisfaction of no-through boundary condition, the resultant superposition of potential and rotational velocity fields does not impose any alterations on the spanwise velocity magnitudes.
Moreover, the velocity error distributions at various heights from the hill surface demonstrate the error reduction achieved via ImVIC+ and ALE-VIC+ compared to linear interpolation and VIC+ (figure 7). The differences in error reductions become more distinguishable at higher seeding concentrations, with the ALE-VIC+ approach yielding the minimum errors in the near-wall region. The differences in grid locations between the uniform and conformal grids (figure 4) causes the path of the optimization process to vary between the two grid formations. Due to the extreme scarcity of particles in the low particle concentration case (N P = 0.01 ppp), this difference causes variations of the reconstructed velocity field also away from the wall. The no-slip boundary condition defined over the conformed hill surface restricts the capability of ALE-VIC+ to properly determine the velocity variations within the separation region, thus yielding underestimations of streamwise velocity at the crest of the hill (figure 5, left-6th row). At this low ppp value, only ImVIC+ reveals a significant improvement of the reconstruction accuracy over VIC+. Nevertheless, as the particle concentration is increased (N P = 0.05 ppp), the perfectly conformed solid surface allows exact boundary conditions over the hill surface to be defined, thus enabling ALE-VIC+ to surpass the reconstruction accuracy of both VIC+ and ImVIC+ ( figure 8, left). For the highest N P value, the density of LPT information closer to the hill surface is considerably increased, which causes the accuracy gap between the ImVIC+ and VIC+ methods to collapse within a shorter range from the hill surface (figure 8, middle & right).
Finally, it should be noted that, whilst the modifications achieved by the proposed approaches are apparent both qualitatively and quantitatively, they are observed to be relatively mild. This is due two main reasons. Firstly, the periodic hill structure is located in close proximity of the computational domain boundaries. Secondly, the potential component of the flow is inferior in magnitude to the turbulent structures represented by means of the rotational component over the Helmholtz decomposition owing to the separation effects downstream of the hill. Accordingly, the corrections induced for the secondary flow structures of normal and spanwise flow components is relatively small which can be deduced from the correlation of velocity reconstruction errors correlating with the streamwise velocity components especially closer to the hill surface ( figure 7). Furthermore, separation effects correspond to a turbulent flow behavior which cannot be resolved via potential flow definitions solely. Nonetheless, the appropriate combination of vortex-in-cell approach with the immersed boundary treatment provides the superposition of necessary potential and rotational flow components to reconstruct the flow field properties with increased coherence levels to the physical flow behavior in close proximity of the hill surface.

Experimental assessment
As a part of the Holistic Optical Metrology for Aero-Elastic Research (HOMER) European Union H2020 project, the experimental setup employed in this work is designed to investigate fluid-structure kinematics of turbulent boundary layer-unsteady panel interactions where the experiments are conducted in a low-speed wind tunnel of TU Delft High Speed Laboratory.

Experimental setup
The experimental setup is designed with a square elastic membrane exposed to turbulent boundary layer conditions and equipped with a 3D-LPT image acquisition system (figure 9). The panel to be deformed is a square elastic membrane with sufficiently high moment of inertia to prevent any aeroelastic deformations and have full control over the membrane shape. It is clamped on an aluminum frame of 60 × 60 cm 2 from all edges to restrict its motion to only elastic deformations and avoid any transnational motion. A DC motor is connected to the center of the membrane by means of a gear and rod mechanism (figure 10, top-right). The membrane is steadily deformed at three different stages of upward, middle and downward positions, and unsteady deformed while being actuated at three different frequencies of 1 Hz, 3 Hz and 5 Hz with an amplitude of 30 mm from valley to crest (figure 11). The corresponding reduced frequencies (k = ω × b/V ∞ ) obtained with the membrane motion of 1 Hz and 3 Hz (b = 0.25 m & V ∞ = 12 m/s) are in the order of k 1 Hz,3 Hz ⩽ 0.06 where a quasi-steady state can be assumed. However, the reduced frequency at 5 Hz reaches to a level of k ≈ 0.2 and shall be considered within the unsteady aerodynamics regime (Leishman 2016). Moreover, a black foil with a regular grid of light-grey dots (0.8 mm diameter, 10 mm distance between adjacent dots, 36 × 36 dots) is applied to the upper face of the elastic membrane to enable the structural displacement measurements by means of LPT (figure 10, bottom-left and figure 11, top). A 120 cm long rigid plate with Lego blocks of 9.6 mm height located at 10 cm from its leading edge is installed upstream of the model to ensure well-developed turbulent boundary conditions at the test section (figure 9).  The processing of raw images and the pressure tab data acquired throughout the experimental campaign are performed with the commercial software package DaVis version 10.0.5 from LaVision GmbH and LabVIEW from NI respectively. For reconstructing fluid motion, the acquired images are provided to STB algorithm for LPT (Schanz et al 2016). Accordingly, the particle tracks are reconstructed with a minimum trace of seven particles which are regularized with 3th order polynomials for the computation of velocity and acceleration information. Furthermore, in order to extract the motion and deformation information of the membrane surface, the superior light scattering intensities of the surface markers are exploited. Since the individual analysis of the membrane motion requires the flow tracers to be removed from the images, a high-pass filter is applied to isolate the illuminated surface tracers. The corresponding pre-processed images are processed using the STB algorithm as well with the same minimum track length and regularization constraints.

Instantaneous flow characterization
Referring to the aforementioned discussion on the complete unsteady nature of FSI problems, the local closure of Collar's triangle of aeroelasticity demands a fully time-resolved characterization of the flow state. However, the employed largescale experimental setup resulted in a tracer particle concentration of C = 100 par/δ 3 99 where the thickness of the fully developed turbulent boundary layer just upstream of the membrane was recorded as δ 99 = 0.06 m (figure 13). Hence, with the corresponding image particle density of N P = 0.002 ppp, the instantaneous flow field data suffers from scarcity of tracer particles which are also scattered over the measurement domain. Although the measurement volume has a size of 20 × 20 × 15 cm 3 , in order to reconstruct the flow properties where accurate description of boundary conditions is required, the region of interest (ROI) is restricted to a domain of 10 × 10 × 8 cm 3 , where sufficient particle concentration is obtained.The relative position of the ROI with respect to the elastic panel is demonstrated in figure 11 for different membrane deformation stages.
The structural motion of the membrane captured via the LPT system and STB algorithm (as mentioned in section 4.1) is employed to reconstruct the membrane shape and deformation level both for steady (figure 11) and unsteady (figure 12) actuations. In order to quantify the accuracy of the structural motion reconstruction, the deformation of the central location of the membrane captured via the optical measurement setup is  compared against the known actuation at the pressure tap location. Accordingly, figure 12 shows the comparison of measured and induced motion at three different frequencies (1 Hz, 3 Hz and 5 Hz) over one period of each motion cycle. The RMS error levels of the membrane central position over the actuated motion are below ∼0.8, ∼0.9 and ∼0.8 mm, respectively. The instantaneous reconstruction characteristics are analyzed over the velocity and pressure information, by processing the time-resolved LPT data obtained for the steady and unsteady membrane deformations. For the steady deformation stages of the membrane, sampling rate for the instantaneous densely interpolated flow fields is determined based on the turbulent boundary layer profiles just upstream of the elastic membrane. The integral time scale is approximated by T I = δ 99 /U inf . In accordance to the boundary layer properties obtained prior to this experimental campaign, the integral time scale is T I = 0.06/10 = 6 ms. Hence, the corresponding the optimal sampling rate is SR = 1/2T I ≈ 83 Hz. As the measurements with steady membrane deformation are conducted Figure 13. Number of particles captured over the central location of the membrane at downward deformation position for the ensemble of STB output (left). Available particle tracking information (color-coded by streamwise velocity) for phase averaged ensemble of STB output (middle) and instantaneous STB output (right). The grey particles denote the track history of the captured LPT data. for 3 cycles with 3 kHz acquisitions of 1 seconds each governed by the temporary storage allowance on the cameras, a total of ∼240 statistically independent time instants are post-processed with the data assimilation algorithms for statistical analysis of dense flow field interpolations. Moreover, in order to provide a reference for assessment of the reconstruction accuracy achieved with various data assimilation approaches, the available STB data computed from a total of 6000 time instants is ensemble-averaged over bins of r STB = 2.5 mm (r STB = ∆x). The corresponding particle counts with a minimum of ∼250 tracks is regarded as a reliable source of information to represent the time averaged flow characteristics of the TBL interactions with the deformed membrane ( figure 13, left). However, this can only be achieved for the steady cases due to the lower number of total measurements for the unsteady deformations. Accordingly, the number of time instants corresponding to the upward, neutral and downward position of the membrane deformation stages for each frequency were determined by the membrane motion frequency for constant duration of measurement. Hence, the number of measurements are varied for different frequencies to match the total number of samples for each frequency. Therefore, a total 15 time instants are obtained for each membrane deformation stage for all three different motion frequencies of the membrane where the number of particles were too low for to perform a reliable ensemble average.

Velocity field reconstruction.
Initial investigations of the various data assimilation approaches are performed by extracting the time averaged streamwise velocity profiles within the turbulent boundary layers for steady and unsteady membrane deformations at upward, neutral and downward positions. Starting with the profiles obtained using the AGW method, the overestimation of the streamwise velocity values in close proximity to the membrane surface is apparent ( figure 14). This is mainly attributed to the low availability of the particle tracks in the vicinity of the membrane surface. Due to the large size of spherical bins (r AGW = 10 mm), the flow field information assigned to the grid locations close to the surface of the membrane are dominated by the particle tracks further away from the central location of the bins. Therefore, even though the influence of particle tracks are inversely related with Euclidean distance, frequent absence of any particles in the vicinity of the membrane causes the velocity values further away from the surface to dominate the resultant flow field data in close proximity to the surface.
On the other hand, linear interpolation provides the opposite characteristics with severe underestimations of the velocity profiles in close proximity to the membrane surface. This is again an influence of the scarcely present particle tracking data in close proximity of the membrane surface. Therefore, without any particles to provide the necessary velocity data in the lower layers of the turbulent boundary layers, the streamwise velocity magnitudes within the boundary layer profiles are significantly downgraded. Nevertheless, towards the locations further away from the membrane surface, the instantaneous availability of the particle tracking data increases which allows both linear interpolation and AGW to capture the boundary layer profiles with considerable accuracy. The improving accuracy towards the freestream is also a result of the physical characteristicsracteristics of the turbulent boundar of the turbulent boundary layers where the velocity gradients within the wall normal direction are concentrated in the vicinity of the membrane surface and possess considerably low amplitudes throughout the outer layers (Clauser 1956).
Furthermore, as the ALE-VIC+ approach is used to reconstruct the dense flow field information, a boundary fitted grid structure is created where the membrane surface is conformed by the mesh formation. Since the boundary for the dense interpolation domain is defined over the membrane surface, a no-slip condition is required to be implemented as a boundary condition. Hence, the neglected viscosity effects within the vorticity transport equation and the absence of particles in close proximity of the membrane surface, the implementation of a no-slip boundary condition also caused the ALE-VIC+ approach to underestimate the velocity magnitudes closer to the membrane. Nevertheless, the underestimations of the streamwise velocity components are considerably alleviated compared to the linear interpolation by means of the physics based interpolations of the VIC+ framework. This allowed the velocity gradients in close proximity of the membrane to be calculated with greater intensity compared to the linear interpolation which is also referred to as the more accurate resolution of turbulent boundary layers characteristics.
Moreover, the application of standard VIC+ and ImVIC+ method involves a similar approach of treating the lower boundary of the dense interpolation region which corresponds to the interior locations of the elastic membrane. Thus, the boundary conditions attributed to those locations do not necessarily refer to a no-slip boundary condition implementation over the membrane surface and allows the flow properties to vary according to the locally available particle tracking information. Therefore, especially throughout locations over the membrane where the deformation displacements are large, the underestimations of streamwise velocity magnitudes are further diminished compared to ALE-VIC+ using the standard VIC+ approach ( figure 14, top and figure 16, top). As the displacement amplitudes are reduced at the downward position of the membrane, the difference between the standard VIC+ and ALE-VIC+ is minimized since the grid layouts of the two approaches became almost identical. Hence, the resultant velocity boundary layer profiles revealed significant similarities between the standard VIC+ and ALE-VIC+ for steady and unsteady membrane deformation stage of downward position ( figure 14, bottom and figure 16, bottom).
Furthermore, a similar comparison is also valid for between the ALE-VIC+ and ImVIC+ approaches in terms of the influence of a boundary condition imposition directly over the membrane surface. Nonetheless, the ImVIC+ approach provided the necessary acceleration on the streamwise velocity profile within the close proximity of the membrane surface. This is achieved by enforcing the no-through boundary condition which allows the immersed boundary treatment to incorporate the relevant physical impact of a solid boundary presence on the VIC+ framework. The intensity of the induced acceleration is elevated as the surface deformation is increased to create a greater curvature of the membrane, cambered towards the incoming flow field (figure 14, top and figure 15, top). Thus, the corresponding dense flow field interpolations are well adapted with ImVIC+ to the modifications on turbulent boundary layers due to the changing membrane shape. The significance of the implementation of the immersed boundary treatment is further highlighted as the modifications on the dense flow field over the standard VIC+ approach concentrated in the close vicinity of the surface where the absence of the particle tracks is considered as a major source of error that prevents capturing accurate boundary layer characteristics as deduced by .
Finally, the physical characteristics of the turbulent boundary layers reconstructed at different membrane positions both with steady and unsteady deformations of the membrane are analyzed. The local cross sectional area of the flow streamlines are contracted as the intrusion of the membrane into to the flow field is increased from downward to upward deformation stages. Moreover, due the relatively low speed of the incoming flow (M ∞ = 0.035<< 0.3), the validity of incompressibility assumption is well preserved (Marchioro and Pulvirenti 1994). Thus, the conservation of mass imposes an acceleration of the flow over the membrane surface, governed by the positive gradients of membrane displacement within the streamwise direction. Accordingly, the freestream velocity values are observed to be increasing with increasing deformation levels as expected for both steady and unsteady deformation cases (figures 14 and 16). Moreover, with the increasing velocities of the freestream values over the boundary layers, the boundary layers are contracted in size which is in agreement with the turbulent boundary layer theory (Schlichting and Gersten 2017). Furthermore, the effect of elevated freestream velocity values and the reduced thickness of boundary layers are captured both between different deformation stages and locally over the membrane surface as the curvature of the membrane yields spatial variations of intrusion amplitudes towards the flow field.
A major difference is observed in terms of the velocity gradients especially in close proximity of the membrane surface between the steady and unsteady membrane deformation. As the membrane motion frequency is increased, the particle residence times within the measurement domain become comparable with the membrane speed which is also observed as the reduced frequency transitions from quasisteady to unsteady regimes as motion frequency is increased from 1 Hz to 5 Hz. Hence, for the 1 Hz case the effect of various membrane deformation stages becomes more profound which reduces as the motion frequency of membrane is increased where the difference between the boundary layer thicknesses at different stages of membrane motion also decreases. The gradual reduction of streamwise velocity gradients in the surface normal direction emphasizes this behavior especially for the upward deformation position where the variation of the velocity gradients and the corresponding boundary layer thickness is the greatest ( figure 15).

Pressure field reconstruction.
The pressure fields are computed by relating the material accelerations to the static pressure distributions over the Navier-Stokes equation (Marchioro and Pulvirenti 1994) while neglecting the viscous diffusion terms due to the low order of magnitude, O(10 −5 ), of their influence for turbulent flow conditions (Murai et al 2007). While the VIC+-based data assimilation approaches already provided material acceleration information as an output of the dense interpolation scheme, the STB and linear interpolation based pressure information is computed by utilizing the raw material acceleration. For the STB, the acceleration terms are binned similar to the velocity, thus yielding the pressure gradient field from equation (15). For the linear interpolation, the raw acceleration data computed for the instantaneously available tracks is linearly interpolated within the spatial domain, which also yielded the pressure gradient information over the reconstruction volume. Then, the computation of pressure distribution over the measurement domain is performed by spatially integrating the pressure gradients by means of solving a Poisson equation accompanied with the proper boundary conditions due to its superior accuracy (Charonko et al 2010).
The reconstruction of surface pressure over the nonuniformly deformed elastic membrane is performed utilizing two different approaches depending on the data-assimilation framework in use for elevating the spatial resolution characteristics. For the VIC+ and ImVIC+ employment, the rectangular uniform mesh formation is preserved throughout the data assimilation approach where the surface locations over the elastic membrane left hanging resulting in an intruding formation over the mesh structure. Hence, the pressure reconstruction is performed over the uniform Cartesian grid with masking the internal region of the solid intrusion. Then, the static pressure values over the exact surface locations are computed with an omni-directional integration procedure similar to the approach proposed by Jux et al (2020). In case of linear interpolation and ALE-VIC+ application, the resultant flow field information is represented over a boundary fitted grid structure conforming the membrane surface. Therefore, the pressure distribution over the membrane surface is computed using a reconstruction scheme with curvilinear transformations proposed by Cakir (2020). The letter approach enabled direct computation of the surface pressure alongside the global pressure information without the need of external extrapolations.
A significant improvement of accuracy and resolution is achieved by employing the VIC+ variants for the reconstruction of pressure information over the membrane surface. All VIC+ variants revealed similar properties in terms of capturing a lower pressure region towards the membrane central location. This corresponds to an accurate representation of the underlying physics of the problem where for all three deformation stages an upward curvature is present and the aforementioned acceleration in the direction of positive surface displacement gradients is accompanied with a reduction in static pressure (figure 11). Difference between the standard VIC+ approach and the proposed variants are observed to be amplified as the membrane deformation level increased. Nevertheless, even for the downward membrane position, a slight improvement in the reconstruction of pressure gradients is observed with both ALE-VIC+ and ImVIC+ over the standard VIC+ approach (figure 16, 3rd-5th columns-top). However, both proposed approaches are observed to be underestimating the pressure gradients over the membrane surface in comparison to the ensemble averaged STB data which yielded the static pressure at the center of the membrane to be overestimated with a value of ∼1.6 Pa (table 2). Moreover, as the membrane's intrusion into the flow field is increased, amplitudes of the spatial pressure gradients over the membrane due to the variations in surface deformation levels are captured with elevated accuracy with both ALE-VIC+ and ImVIC+ methods (figure 16, 3rd-5th columns-middle & bottom). Owing to the presence of boundary layer profiles, the streamwise velocity distributions reveal to lower amplitudes over the central region of the membrane. The underestimations of streamwise velocity gradients in the surface normal direction in close proximity to the membrane surface associated with the application of linear interpolation and standard VIC+ deviate the flow field characteristics from the condition of zero pressure gradients in the surface normal direction which is valid for turbulent boundary layers (Fernholz and Finleyt 1996). The violation of this assumption intensifies in the case of linear interpolations promoted by the increasing amplitude of velocity reconstruction errors. Hence, as the velocity gradients are computed locally from these velocity profiles the corresponding pressure gradients are completely inaccurate to capture the correct pressure distribution over the membrane surface. Nevertheless, employment of a surface treatment approach compensates the absence of particle tracks in the vicinity of the surface. In this regard, ImVIC+ provides lower deceleration over the central region of the membrane due to the modifications induced over the velocity field in order to satisfy the no through boundary condition.
Similarly when the VIC+ variants are employed, velocity and vorticity gradients are discretized by means of locally defined finite difference approximations which compose the temporal and convective components of the material acceleration terms. Hence, as the governing equations of the VIC+ framework do not possess any terms to handle turbulence or viscosity, the corresponding material derivative computations become vulnerable against the inaccuracy issues related to the Eulerian approach whilst reconstructing dense flow field data within the turbulent boundary layers. Therefore, the accuracy of acceleration field reconstructions with VIC+ variants is dictated by the balance between the achievable mesh resolution and the dominant scale of the turbulent structures. Deriving an analogy to the large eddy simulations where the local mesh resolution is determined by the dominant large scales within the turbulence cascade (Schmitt et al 1986), VIC+ framework responds to the variations in flow velocity in a similar manner. In absence of any turbulent treatment, the mesh resolution of data assimilation acts as a low-pass filter removing the influence of sub-grid scale flow field variations from the densely interpolated velocity and acceleration data. However, the achievable mesh fidelity for the VIC+ variants is dictated by the instantaneous density of particle tracks which is a function of the experimental setup parameters and independent of the flow behavior. Therefore, in case of large scale vortical structures, the variations in the velocity field can be captured by the direct description of the vorticity transport equation which also yields more accurate reconstruction of the material acceleration terms (van Gent et al 2017).
Nonetheless, the influence of acceleration terms on the cost function is controlled by a scaling factor determined according to the relative standard deviations of velocity and acceleration terms. This is referred to as the balance between the aforementioned balance of large and small scale flow field variation as well as an indication of their respective uncertainty levels. Therefore, the scaling factor strongly favors the velocity terms especially in case of turbulent boundary layers. Hence, even though the corresponding pressure gradients obtained from the VIC+ variants are substantially more accurate than the ones computed over an Eulerian approach, the influence of the Eulerian framework can be observed.
Comparing the surface pressure information over unsteady deformations of the membrane, the variations throughout the overall pressure profiles are observed to be relatively confined. Nevertheless, both ALE-VIC+ and ImVIC+ enabled a greater agreement with the pressure tab measurements by minor improvements (figure 17). Three main reasons can be deduced for these minor variations in correlation with the validation studies. First of all, the spatial gradients of both velocity and pressure over the measurement domain are dominated by the controlled motion of the elastic membrane. Hence, as the motion frequency of 1 Hz and 3 Hz corresponds to reduced frequencies of k ⩽ 0.05 (within the quasi-steady regime), these gradients can be reconstructed even with the linear interpolation without severe loss of accuracy. Secondly, considering the assumption of zero pressure gradients within the turbulent boundary layers and the membrane motion, the Dirichlet boundary condition dictates the majority of pressure variations in time. Although there exists non-zero pressure gradients throughout the realistic boundary layers, these relate to pressure modulations in the order of 6-8 Pa while the absolute pressure difference captured at highest deformation state of the membrane is ∼30 Pa. Hence, the pressure variations aimed to be reconstructed with a superior accuracy by the proposed algorithms correspond to ∼25%. Lastly, the boundary conditions for VIC+, ALE-VIC+ and ImVIC+ methods are determined via linearly interpolating the STB data over the corresponding boundaries. Considering the close proximity of the Dirichlet boundary condition of pressure reconstruction to the computational domain boundaries of the VIC+ variants, the differences between the various approaches in terms of the absolute pressure values are further alleviated. Finally, specific time instants are identified where the differences between VIC+ and the proposed variants are amplified with reconstruction errors up to ∼70% of the pressure fluctuation amplitudes obtained using VIC+. This is related to the extreme sensitivity of the standard VIC+ method to the spatial distribution of the available particle tracking information. In the absence of a proper boundary condition definition in close proximity of the elastic membrane, the VIC+ algorithm becomes strongly dependent on the presence of particles in the near wall regions to drive the optimization procedure towards an accurate solution which is a considerably rare situation in practical cases. Therefore, the ALE-VIC+ and ImVIC+ approaches provide the capability of accurately and consistently reconstructing the flow properties even in the cases of complete absence of particles in the near wall regions by providing a kinematic characterization of the solid boundary intrusions. The resultant accuracy improvements achieved by the implementation of ALE-VIC+ and ImVIC+ over linear interpolation and VIC+ are also observed over the cumulative error levels over the pressure reconstruction profiles (table 3).

Conclusions
Available spatial resolution characteristics for time-resolved flow field measurements with large scale LPT techniques are restricted due to the tracer particle specifications of HFSB tracers. However, the available governing equation based data assimilation techniques such as FlowFit  and VIC+  enable dense volumetric interpolations of flowfield information for regions of uniform rectangular computational domains with sole fluid presence. In this regard, the introduced variants of ALE-VIC+ and ImVIC+ provide the standard VIC+ algorithm with the capability of incorporating appropriate boundary condition definitions for dense flow field characterization in close proximity of solid objects with generic geometries. The ALE-VIC+ method implements the ALE method of fluid-structure interaction frameworks with a RBF based mesh deformation scheme to ensure the adaptability of the grid formations to the unsteady deformations of the FSI interface. On the other hand, ImVIC+ approach preserves FFT based Poisson solvers to increase computational efficiency using uniform predefined computational grids where immersed boundary treatments are utilized to satisfy the boundary conditions by means of additional flow components.
The validation studies of the proposed methods are performed with a numerical test case of flow over periodic hills, where the DNS datasets are manipulated to stimulate realistic experimental conditions. Even though, both linear interpolation and VIC+ variants resulted in coherent flow behaviors with the hill form, the local variations of velocity components favored the latter in terms of greater reconstruction accuracy. With the application of ALE-VIC+ and ImVIC+, reconstruction accuracy improvements over the standard VIC+ method are achieved especially in close proximity of the hill surfaces in terms of streamwise and normal velocity components composing the surface flow penetrations. Accuracy of the modifications are observed to be favoring the immersed boundary treatment approach with ALE-VIC+ containing additional numerical errors involved within the mesh transformation procedure (Lee and Tsuei 1992) while preservation of uniform grid structures in the case of ImVIC+ enabled minimization of the computational cost. Furthermore, ImVIC+ achieved better compatibility with the standard VIC+ approach with modifications being confided to the close proximity of the hill surface where the scarcity of particle tracking information prevented the construction of an accurate objective function for the optimization procedure.
Finally, both data assimilation approaches are applied to enable instantaneous flow field characterization for the measurements of turbulent boundary layer interactions with a dynamically deforming elastic membrane. The densely reconstructed flow field properties are then employed to compute the pressure distribution over membrane surface, revealing the time-resolved interaction between the flow structures and the membrane deformations. The ALE-VIC+ and ImVIC+ methods enabled the kinematic discretization of the unsteadily deforming membrane motion. Thus, the non-intrusive reconstruction of instantaneous velocity and pressure field computations are achieved compared to the trilinear interpolation and standard VIC+ approach. In comparison to surface pressure data, ImVIC+ approach also revealed a slight superiority of accuracy over ALE-VIC+ benefiting from the absence of numerical errors associated with curvilinear grids and the enforced no-slip boundary condition over the membrane surface. Nonetheless, this issue can be addressed by including the viscous terms in the NS equations provided that the particle concentrations and the corresponding spatial resolution is sufficient. In absence of these modifications, due to its versatile applicability to very complex geometries without tremendous increase of computational cost, makes ImVIC+ a potential method of choice over the ALE approach for providing physics-driven VIC+ method with the capability to accurately increasing the spatial resolution of the LPT data in the presence of generic solid boundaries.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.