Numerical calculation of the turbulent flow past a surface mounted cube with assimilation of PIV data

The application of Data Assimilation (DA) methods in Computational Fluid Dynamics (CFD) problems is a concept actively being explored to couple CFD with Experimental Fluid Dynamics data. Here, Particle Image Velocimetry (PIV) data are assimilated in an OpenFOAM based CFD solver to calculate the velocity and pressure fields of the turbulent flow past a surface mounted cube inside an atmospheric boundary layer for three planes belonging to the symmetry plane of the flow. At first, the SIMPLE algorithm is used to correct both pressure and velocity fields, with the PIV data used to formulate the initial and boundary conditions. The Reynolds stresses are calculated directly from the PIV data instead of using a turbulence model. Next, we use two implementations of the nudging method and two formulations of the Kalman Filter in order to assimilate the PIV data into the iterative SIMPLE procedure. A grid independence study is performed, and the performance of the different methods is assessed. The CFD predicted pressure field is in good agreement with pressure measurements on the cube surface. The results also show that the SIMPLE based correction step already leads to a significant reduction of both the mean and the variance of the continuity errors as well as the difference between the original PIV data and the resulting velocity fields. The application of the DA methods, particularly the KF, leads to minor further improvement of the results but does improve convergence of the CFD solver.


Introduction
Both Experimental Fluid Dynamics (EFD) and Computational Fluid Dynamics (CFD) are ways to get invaluable insights regarding fluid flows for a wide range of scientific fields, ranging from aeronautics and wind energy to flows around buildings and microfluidics.EFD and CFD methods have been independently developed for a long time, but their interaction has been limited to the validation and/or calibration of CFD models using EFD results or complementing EFD data using CFD results.However, in the last several years more and more attention has been paid to coupling EFD and CFD to overcome their individual weaknesses [1][2][3].
One approach to achieve this integration is by incorporating the measurements into the CFD solvers through the use of state observers, a widespread method in control theory.For instance, in [4], a state observer was used to reconstruct the boundary conditions for the turbulent flow through a square duct using the SIMPLER algorithm.Similarly, in [5], an observer method for the pressure difference between CFD and experimental data was applied to the two-dimensional flow around a square cylinder with von Kármán vortex shedding.Recent advancements include the application of an observer based method to the pressure equation in studying the flow around a square cylinder at Re = 100 using synthetic data [6], as well as the utilization of two state observer approaches to analyze the velocity difference between URANS and experimental data for the 3D flow behind a generic car mirror model [7].
In addition to state observers, nudging based methods, a particular type of state observer, have gained traction in recent years [8,9].These methods offer a straightforward integration into existing computational codes.For example, [10] employed a nudging method to estimate the turbulent flow around a canonical square cylinder at Re = 22 000 by modifying the momentum equation with DNS data.
Alternatively, a variety of Data Assimilation (DA) methods have been used in order to account for the uncertainties of the measurements that are used as initial/boundary conditions or are integrated into the solver, as well as the uncertainties of the CFD modeling.With these methods, the respective error covariance matrices usually quantify the uncertainties of the discrete sources of information.DA can be viewed as combining observations with model output to improve the latter while considering the uncertainties in both [11].DA methods have been used widely in the field of weather forecasting and earth sciences [11][12][13], but have been used in fluid dynamics as a way to couple CFD solutions and experimental data.
The Kalman Filter (KF), a classic DA approach, has been combined with segregated solvers for two and threedimensional laminar [14] and two-dimensional turbulent flows using synthetic [14,15] or experimental [14,16] data.Even reduced-order KF have been used [17] to combine Particle Tracking Velocimetry (PTV) data and DNS to analyze a planar jet.More sophisticated methods, such as the Ensemble KF (EnkF), have been applied to two and three-dimensional transonic flows using experimental data in [18], while in [19], the EnKF coupled with a high-order CFD code was applied to one and two-dimensional problems.Similarly, in [20], the EnKF method was employed to incorporate synthetic velocity data and wall or near wall quantities (eddy viscosity, friction coefficient, and wall pressure), in order to reconstruct twodimensional turbulent flows.
In [21], an additional forcing term calculated from a variational procedure using a direct-adjoint approach, was added to the Navier-Stokes (NS) equations to integrate timeaveraged Particle Image Velocimetry (PIV) data of the twodimensional mean flow field around an airfoil at Re = 13 500.Variational DA methods have also been used to generate inflow conditions for Direct Numerical Simulations (DNS), using PIV data collected from wind tunnel and synthetic measurements [22].However, in [23], the reverse was applied, as the mean flow field of the two-dimensional flow around an infinite cylinder at Re = 150 was reconstructed by assimilating DNS data.
In this work, both nudging and KF are applied for the calculation of the turbulent flow quantities around a surface mounted cube, using PIV data within the same computational and measurement domains.Determining the pressure field both on surfaces of solid bodies and within flow fields is of great importance for a wide variety of applications.For example, cavitation and aeroacoustic phenomena, the calculation of power curves as well as lift and drag coefficients of wind turbines and wind loads on structures [24,25].
Usually, the surface pressure is experimentally obtained using pressure orifices in the model, connected to pressure transducers or microphones.However, determining the instantaneous, or even the time-averaged, pressure fields is significantly more challenging.For instance, widespread pressure measurement techniques such as multiple-hole pitot pressure probes or fast-response pressure sensors provide singlepoint measurements with limited ability to reconstruct a pressure field with adequate spatial resolution.Additionally, their intrusive nature hinders their ability to capture dynamic phenomena [24] without interfering with the flow.Non-intrusive surface pressure measurement techniques exist, such as using pressure-sensitive paint [26] or employing air bubbles to measure static pressures [27].However, the applicability of these methods is narrow, as they usually fail to capture the whole pressure field in the domain of interest and cannot provide velocity and pressure measurements simultaneously.
The PIV measurement technique is a non-intrusive experimental technique that can extract instantaneous velocity fields over the whole domain of interest [28,29].As a result, various efforts have been made toward developing PIV based pressure calculation methods.The scope of these methods is to combine fundamental conservation equations with the PIV velocity data to provide their respective pressure fields.Assuming that the PIV data are mass-conservative, the solution of a Poisson-type pressure equation may lead to the extraction of both instantaneous [30,31] and time-averaged [32] pressure fields.Similarly, the use of Taylor's frozen turbulence hypothesis combined with a Poisson equation has been proposed to extract instantaneous pressure fields [33].
In [34], the solution of a Poisson pressure equation and line integration coupled with velocity smoothing techniques were applied in time-dependent incompressible flows using Digital-PIV data.Similarly, both the direct solution of the NS equations and the solution of a Poisson pressure equation were investigated in [35] for the unsteady wake flow around a cylinder of square cross-section, using Stereo-PIV data.Finally, a more recent approach [36] proposes a mesh-free method to solve a Poisson equation by interpolating the PIV velocity data using Radial Basis Functions.
Several of the studies presented above determine pressure from PIV.Some of these (e.g.[21,22]) assimilate PIV measurements to determine velocity and pressure fields in computational domains significantly larger than the measurement domains.Others (e.g.[31,33,36]) utilize only the PIV data and calculate the corresponding pressure fields without altering the velocity field, with the implied assumption that mass conservation was inherent in the PIV data.Another approach that incorporates PIV data into a CFD solver and simultaneously reconstructs the respective pressure field, consists of SIMPLE based methods where the initial and boundary conditions of the velocity are constructed from PIV data.An added benefit of combining PIV data with a solution for the pressure as well as the velocity fields is that both the continuity and momentum conservation laws are ensured.This is often a problem area for PIV data [29].Pressure fields from two-dimensional PIV snapshots have been extracted using the SIMPLER algorithm on a non-staggered grid for laminar, incompressible, and steady flow in [37].Similarly, in [38], two-dimensional PIV data were utilized to compute pressure fields for laminar, incompressible, steady, and unsteady flows.The use of two-dimensional PIV data to calculate pressure fields for turbulent, incompressible time-averaged flows was recently achieved using the Reynolds stresses from the PIV statistics and including them in the source terms of the RANS equations [39].
In this work, a SIMPLE based numerical solver is developed using the OpenFOAM CFD toolbox [40][41][42] to correct the velocity data and extract time-averaged pressure fields for the atmospheric boundary layer turbulent flow past a surface mounted cube using two-dimensional PIV velocity data [43,44].The method is similar to the previously mentioned applications of the SIMPLE method [37,38], but here it is applied to a turbulent flow with direct insertion of the Reynolds stresses in the equations, as they were measured by the PIV method rather than applying a turbulence model.The application is similar to that in [39], but here, a different code, SIMPLE algorithm variant and boundary conditions are applied.After the solver is verified based on experimental data [43,44], both nudging and Kalman filtering methods are applied.The goal is to reduce the variance of the continuity errors and the discrepancies between the initial experimental velocity fields and the final corrected ones.The calculation of pressure fields from PIV data while taking into account both the experimental and CFD errors is not new, but here we incorporate the data directly into an iterative procedure that solves the momentum and continuity conservation equations.This approach ensures that the resulting fields are both mass and momentum conservative, while simultaneously, the discrepancies introduced by the original experimental data and the CFD modeling are reduced.To the authors' knowledge, this is the first time that pressure calculation using an iterative solution of the NS equations from PIV data has been coupled with Data Assimilation methods, in order to assimilate the PIV data to the CFD results of a turbulent flow field throughout the iterative process.
The paper's structure is as follows: firstly, in section 2, the methods used in the solver implementation are thoroughly described.Next, in section 3, the developed solver is verified, a grid independence study is performed and the impact of DA methods is investigated.Finally, the conclusions of this work are presented in section 4.

Governing equations
The finite volume method is applied to solve, using iterative methods, the two-dimensional steady-state Reynoldsaveraged Navier-Stokes (RANS) equations.These are the mass conservation equation (continuity equation) (1a) and the momentum conservation equation (1b), expressed for incompressible, steady-state flow: where ū and p are the Reynolds-averaged velocity vector and pressure respectively, T is the Reynolds stress tensor, ∆ is the Laplace operator, ρ is the density and ν is the kinematic viscosity.For brevity reasons, the overbar symbols are omitted throughout this paper.Usually, in RANS approaches, eddy viscosity models are used to calculate the Reynolds stress tensor T. This work uses a different approach by calculating T from the PIV data and introducing it as an explicit momentum source in (1b).Considering there are N PIV snapshots, T = T PIV is given by where u ′ k,i = u k,i − ūPIV , i = 1, 2 corresponds to the x or y direction respectively, ūPIV is the mean velocity field and k is the PIV snapshot index.The possible effect of this approach on the model uncertainty is discussed in 2.3.3.

SIMPLE algorithm.
The solution of the RANS equations is performed using the SIMPLE algorithm [45], an iterative procedure that solves the NS equations in a segregated way.Here, the algorithm will be applied for calculating the pressure field and correcting the velocity field.The OpenFOAM toolbox [40][41][42] formulation is used, which is slightly different than the one used in [39] and consists of a momentum predictor, shown in (3), that determines a velocity field u * which satisfies the momentum equation for a given pressure field p k , at the k + 1−th iteration: where u * P is the velocity vector at point P, p P the pressure value at P and, a * P the coefficient of u * in the discretized form.Additionally, E(u  * ) is a matrix that includes the contributions of the neighboring cells of P into the discretized momentum equation, and depends on the discretization schemes being used.E(u * ) also contains the external momentum sources such as the gravity force as well as the control forcing used here but presented later on.Following the solution of (3), a Poisson equation is solved in order to determine a pressure field p * that satisfies the continuity equation for u * and is given by where the f index refers to the cell faces, meaning that u * is interpolated onto the cell faces before solving the Poisson equation.Then, the pressure field is relaxed explicitly with β p to find the pressure field at iteration k + 1 The volumetric fluxes across the cell faces are corrected: where S f is the surface vector at the f face, given by S f = S f n f , where S f is the face area, and n f is the normal vector.Finally, the velocity field is corrected by In OpenFOAM, the Rhie-Chow interpolation [46] is used implicitly to avoid checkerboarding.The upwind scheme [47] is used for the discretization of the convective terms, and a central differencing scheme is employed for the diffusion terms and the pressure gradient.Furthermore, iterative methods are employed to solve the discretized equations in matrix form.The preconditioned conjugate gradient (PCG) method [48,49] coupled with the Diagonal based Incomplete Cholesky preconditioner [48,50], is used to solve the (symmetric matrix) pressure equation.Similarly, the Preconditioned bi-conjugate gradient [48], method coupled with the Diagonal based Incomplete LU preconditioner [48], is used for the solution of the (asymmetric matrix) momentum equation.

Poisson equation.
Alternatively, to the SIMPLE based method [39], a Poisson equation may be solved to calculate the pressure from uncorrected PIV data [24].The pressure Poisson equation is devised by taking the divergence of the momentum equation for the PIV velocity field and assuming that the PIV data conform to the continuity equation: where T PIV is the Reynolds stress tensor available from the PIV data.

Boundary conditions.
For both the SIMPLE based method and the pressure Poisson equation, there is a need for the specification of a computational domain and boundary conditions.In the present case, the finite volume solver uses the (uniform) grids formed by the points where PIV data are available.Using a simple bi-linear interpolation, the nodal values are mapped onto the cell centers, as the solver requires.
This approach ensures that only interpolations are performed, avoiding extrapolating the experimental data.Also, contrary to [39] it ensures that the PIV planes and CFD grids have precisely the same dimensions and boundaries.In figure 1, the cell nodes where the PIV data are known are colored in red, while the cell centers where the values are stored in OpenFOAM are colored in blue.Additionally, in figure 1(b) four cells are colored in red, indicating the PIV nodes used in each interpolation to each CFD cell center.In figure 1(a), the same four cell centers are connected with blue lines to indicate the CFD cell centers that a single node influences when the interpolation is performed.
As the PIV measurement planes contain a fraction of the total flow field [29], the choice of boundary conditions has to be considered.The velocity boundary conditions are selected, in line with previous work [39], as Dirichlet boundary conditions using the PIV data values.
These boundary conditions are expanded into a perimetric zone of two cells in the interior part of the grid, where the velocity field is considered constant and equal to the respective PIV data.For example, in figure 2, the implementation of this zone for the upper boundary of a grid is demonstrated.It is worth noting that this selection generally introduces a zone of high continuity residuals near the boundaries, but their impact on the method's convergence and stability was insignificant.Furthermore, contrary to [39] where the pressure equation is solved only in the internal domain using Neumann boundary conditions, for the pressure predictor in (4), zero-valued Dirichlet boundary conditions are applied on any computational domain boundaries not adjacent to a solid physical boundary, and zero gradient (Neumann) boundary conditions are applied on domain boundaries adjacent to solid physical boundaries.In addition, inside the perimetric zone of fixed velocity, the pressure values are calculated using the fixed velocity values and, consequently are not updated during the iterative procedure.Therefore, the pressure boundary conditions impact the calculation of the (fixed) pressure inside the perimetric zone that acts as a Dirichlet boundary condition for the rest of the computational domain.The final procedure for the SIMPLE based method is summarized in algorithm 1.

Data assimilation methods
The flow field calculation procedure described above may be supplemented with Data Assimilation (DA) methods.Here, we apply the nudging method and the KF method.In the examined application, the flow is steady-state, but because it is solved with an iterative method, dynamic DA approaches can be used, considering the iterations are pseudo-time steps.

Best linear unbiased estimator -BLUE
Before we present the DA methods utilized in this work, we present the notion of the Best Linear Unbiased Estimator-BLUE, because it serves as the first step to determine an optimal (in the sense of minimizing the estimation error variance) gain matrix, a necessary step for the DA methods we have utilized.
A discrete model for the evolution of a physical system from time t k to time t k+1 with no input control may be described by a linear dynamic state equation [11]: where x is the model's state vector with x ∈ R n and M is the corresponding dynamics operator that can be time-dependent.
The f index refers to the forecast state, i.e. the state that is predicted using the model Equations.w k represents the modeling error with an associated error covariance matrix Q k .
Observations or measurements y ∈ R m are defined as where the t index of the x t k vector refers to the unknown true state of the system and H k is an observation operator that maps the n dimensional model space to the m dimensional observation space (usually m ≪ n) and can be time-dependent.v k represents the measurement/observation errors that can be considered a statistical process with a zero mean and an associated covariance matrix R k .
The modeling of w k and v k is an important aspect of DA and they are commonly considered to be independent [11] and with multivariate Gaussian/normal probability distributions Then, an a posteriori estimate x a , can be calculated as a linear combination of the forecast state x f and a weighted difference between the measurement prediction y o k and the forecast state x f k , using the Gauss-Markov theorem [11]: where the difference (y o k − Hx f k ) is called the innovation term and the optimal gain matrix K k is given by the Gauss-Markov theorem where P k and R k are the forecast error covariance matrix and observation error covariance matrix, respectively.In the general formalism of the DA approach, the forecast error covariance includes the model errors (Q k ) and any errors accumulated through the sequential analysis steps.While in the BLUE approach, the forecast error covariance matrix P k is equal to the model error covariance matrix Q k , as the notion of error accumulation cannot be defined when calculating the BLUE estimator without an error projection mechanism, e.g. the KF.
The task of determining the model error covariance matrix is more complex, as discussed in 2.3.3.
In this work, the framework described above treats the PIV measurements as observations, and the state vector is considered to be the velocity vectors on the CFD grid.Regarding the model of the dynamic system, the NS equations are solved using the SIMPLE method, i.e. each iteration of SIMPLE corresponds to a new state of the dynamic system.Thus, an estimate state is reached by combining the observations, i.e. the PIV measurements and the forecast state, i.e. the velocity fields predicted by the CFD.

2.3.2.
Nudging.An alternative method called nudging, or Newtonian relaxation, aims to dynamically adjust the model toward the observations by simply inserting a feedback term into the model equation, which is proportional to the observation and model misfit (similar to the innovation term) and nudges the model state toward the observations [51].The main advantage of nudging is its simplicity, as there is only the need to implement an observation operator and a slight modification of the model equation to add a feedback term.
Here, the steady-state RANS equations are modified by adding the nudging forcing term that is proportional to the misfit between the velocity field provided by CFD and the PIV measurements: where the gain matrix K is non-dimensional as expected from (13), and η is a constant that ensures dimensional compatibility.Due to using the SIMPLE algorithm to solve the RANS equations, the innovation term ηK (ū PIV − Hū) divided by ρ, is added to the momentum source matrix E(u * ), as shown in (3).Thus η units can be determined: However, both K and H are non dimensional, and consequently, The nudging method was implemented in the developed OpenFOAM based solver.In this kind of framework, the specification of K is crucial; it can be selected to be diagonal or even scalar, but the effect of it depends on the structure of the observation errors, their correlation in space, and/or their variance in time.But its tuning is not automated and relies on numerical experimentation.The BLUE approach can be used if information about observation errors is available, or an optimal nudging approach [52,53] can be used to minimize, with respect to K, a cost function that represents the error between the model and observations.In this work, two different options for selecting the gain matrix K are implemented.

Identity nudging.
In this application, the values have already been interpolated from all the PIV points onto all the cells of the CFD grids, as seen in figure 1.Even when grid refinement is performed, the PIV data are again mapped to all the points of the finer grids.Thus, H is considered to be the identity operator: The easiest option is to avoid assumptions about the error structures and use the identity matrix as the gain matrix.In this approach, the numerical solution and observation misfit are constantly forced into the system equations without any other considerations.This approach is the least effective but also provides a baseline to see if nudging has improved the final solution.Thus, the gain matrix is given by 2.3.2.2.BLUE nudging.Next, as described in (13) and taking into account (16), the BLUE estimator is used to calculate a gain matrix that minimizes the error variance: where P and R are the forecast and measurement error covariance matrices, respectively.As noted in the previous section, in the BLUE approach, P is equal to the model error covariance matrix Q.
The measurement/observation error covariance matrix is determined using a simplified form of the Reynolds stress tensor, i.e. assuming the velocity measurement errors are uncorrelated in space [31], the measurement error covariance matrix is considered diagonal and is given by where δ ij is the Kronecker delta and the standard deviation is given by σ 2 PIV,i = T i,i .However, PIV data contain random noise, that cannot easily be quantified.The Reynolds stress tensor used in this work also contains physics based information that stems from the turbulent nature of the flow.No effort has been made to a posteriori decompose the calculated Reynolds stress tensor to its physics based and random noise components.Consequently, it may be considered that this random uncertainty is accounted for, although overestimated.A more complex way to estimate the model error covariance matrix would most probably affect the performance of the DA techniques, but it is considered beyond the scope of the present work.
2.3.3.Model/forecast error covariance matrix.The modeling of the model (and consequently forecast) error covariance matrix Q is crucial due to its size and the fact that the true state is, by definition, unknown.Thus, compromises must be made to produce a computationally viable model [11,54].One approach is to separate the information about forecast error statistics from the innovation statistics [55].Alternatively, statistics can be derived for a surrogate quantity [56].
Other low-cost approaches for modeling the Q matrix use reduced bases, factorization, spectral methods [11,17] or simple assumptions about the matrix structure based on the physics of the problem [14,57].In this work, a reduced-cost option is implemented by making assumptions about the matrix structure based on the knowledge from CFD modeling [14].The basic assumption is that Q is diagonal and depends only on the spatial (and temporal for transient problems) discretization.This assumption is reasonable, as the spatial and temporal discretizations, on the one hand, influence the flow features that can be resolved, e.g.turbulence scales, recirculation length, etc.On the other hand, they affect the order of discretization/truncation error of the numerical discretization schemes used.
The main disadvantages of this assumption are that it disregards any errors in the boundary conditions, the modeling of the physics, and the correlation between the errors observed at different grid points/cells.With regard to the boundary conditions, these are fixed at the values from the PIV measurement data and so can be considered to not diverge from them.
Next, regarding the physics of the problem, the conservation equations are assumed incompressible and steady state (i.e.valid assumptions for this flow not involving any numerical modeling) and turbulence is not modeled but included through the PIV data itself i.e. identical to the measurements.In general, in the numerical simulations of turbulent flows, it is the interaction between both numerical errors and turbulence modeling assumptions that will influence the uncertainty of the numerical results.However, the assumptions inherent in the widely used models of turbulence, (e.g.isotropy, turbulence production and dissipation terms), introduce an uncertainty that we here assume to be larger than the one introduced by the direct inclusion of the measured Reynolds stress tensor.Obviously, although it is expected to be more accurate than the turbulence modeling alternative, this approach will introduce an uncertainty of its own.However, quantifying this uncertainty requires more complex and computationally demanding methods, e.g.Monte Carlo simulations, which are beyond the scope of this work.We must therefore acknowledge the possible measurement uncertainties of the Reynolds stress terms but proceed with the present analysis considering only the numerical uncertainties.We consider that this present work can be used as the foundation of future research work, that utilizes Monte Carlo methods, both for determining the sensitivity of the RANS equations to the direct inclusion of the Reynolds stress tensor and determining a more complex (non-diagonal) model error covariance matrix.
Finally, with regard to the error correlations, since the pressure equation is a Poisson equation, an elliptical PDE, meaning that information is propagated instantly throughout the whole domain, we consider this adequate to represent the propagation of the model error without needing a more complex model for the correlation between neighboring cells [14].Thus, a compromise between computational/modeling demands and physics validity is reached.
Based on these assumptions, the model error covariance matrix is considered a power function of the average cell length, ∆x, and the time-step ∆t, as proposed in [14] and given by where the SO (Spatial Order) and TO (Temporal Order) exponents correspond to the minimum order of accuracy of the used discretization schemes.Additionally, C is a user-defined coefficient that ensures dimensional similarity.In this application, the problem is steady-state, and the minimum order of the used spatial discretization schemes is one (O(∆x)), thus the model error covariance matrix becomes

Kalman filter -KF
The KF [58] is usually applied to dynamic problems and aims to minimize the a posteriori error covariance.The KF consists of a prediction/forecast step and a correction/analysis step.At the time t k (which here corresponds to one SIMPLE iteration), the result of a previous forecast is available, x f k , as is a set of measurements/observations, y o k .Based on these two vectors, we perform an analysis that produces x a k .Then, the evolution model is used to obtain a prediction of the state at time t k+1 .Finally, the result of the forecast is denoted x f k+1 and becomes the background, or initial guess, for the next time step.
The goal is to compute an optimal a posteriori estimate, x a k , that is a linear combination of an a priori estimate, x f k , and a weighted difference between the actual measurement, y o k , and the measurement prediction, H k x f k .The weight used is the same as the BLUE estimator defined in (13).Thus, the a posteriori estimate is defined as in (12), while the Kalman gain matrix is given by the BLUE estimator (13) for each time step.Consequently, when the measurement error covariance, R k , approaches zero, the gain K k , weighs the innovation more heavily.On the other hand, as the forecast error covariance, P k , approaches zero, the gain, K k , weighs the innovation less heavily.
The KF aims to propagate the model uncertainty introduced by the model error covariance matrix Q and accumulated in the forecast error covariance matrix P, using M to model the correlation between state variables for consecutive iterations, as seen in algorithm 2. In the spirit of the approximation of the model error covariance matrix in BLUE, the state matrix is also selected to be diagonal and equal to the inverse of the diagonal coefficient matrix, similarly to [14]; thus, it is given by where A is the matrix that contains all the diagonal coefficients a P as shown in (3).It is noted that this approximation is considered to directly propagate the error covariance matrix P where the model error Q is accumulated.Project the error covariance ahead: Compute the Kalman gain: Update estimate with measurement: Update the error covariance: Thus, all the matrices that influence the calculation of K are initialized as diagonal matrices, the M matrix that propagates the forecast error covariance matrix through the iterative process is also diagonal, and consequently, K is also diagonal.This structure of K, on the one hand, is beneficial from a computational cost point of view, and it also allows for different KF formulations, ensuring that they do not disturb the diagonal dominance of the equation matrices used in the SIMPLE method.On the other hand, the calculated gain matrix still minimizes the a posteriori error covariance equation, but the a priori and measurement error structures are simplified.
Two formulations that combine the KF and the SIMPLE algorithms are investigated (i) The classic formulation of the KF (ii) An integrated formulation of the KF [15] 2.3.4.1.Classic formulation.Firstly, the classic formulation is implemented.Each iteration of the SIMPLE algorithm is used as the model that produces the forecast state for the KF.Then, the forecast error covariance matrix and the optimal gain matrix are calculated, and the analysis state is obtained using the innovation term.The procedure is summarized in algorithm 3.One of the main properties of the KF is that it converges, rather quickly, to a steady-state solution [11,58,59] under certain conditions and consequently requires significantly fewer iterations than the SIMPLE algorithm.Thus, it is reasonable to assume that the KF does not significantly affect the SIMPLE convergence.An extensive discussion about convergence and computational cost of the developed methods is held in section 3.
Update state: Update forecast error covariance matrix: end while End

Integrated formulation.
Next, a formulation that combines the two prediction/correction schemes used in this work, the SIMPLE and the KF algorithms, is implemented into a single prediction correction scheme [15].Initially, the velocity field is predicted using the momentum equation and then updated implicitly using the KF gain matrix.The r.h.s of ( 4) is modified by so the new pressure equation is given by: After the solution of ( 24), with the Kalman Gain implicitly included, the volumetric fluxes are corrected, and the velocity field is calculated again, using (7), without explicitly applying the KF.Thus, the augmented velocity field has zero divergence (the mass continuity is enforced explicitly) instead of the intermediate velocity field.This approach is of interest as the non-zero divergence of the PIV velocity field is taken into account only in the pressure calculation, and the asymptotic convergence of the KF reduces the effect of these discontinuities.The integrated procedure is summarized in algorithm 4.

Experimental test case
The methods described in the previous sections are applied to the atmospheric boundary layer flow past a surface mounted cube with vertical openings.Wind tunnel Stereo-PIV measurements are available [43,44] for this case.The measurement planes are along the cube center plane and located upwind, above, and downwind of the cube as shown in figure 3. The flow is turbulent, as the Reynolds number at cube height is Re ≈ 2.4 • 10 4 , with air density ρ = 1.21 kg m −3 and air kinematic viscosity ν = 1.48 • 10 −5 m 2 s −1 .Two different upstream ABL profiles (high and low shear with a power-law exponent of 0.22 and 0.12 respectively), as well as two different cube setups (open and closed cube) are available.In this work, the data from the high shear and open cube setup are used.Furthermore, the Reynolds stress tensor field is calculated from the PIV measurements.
In line with the methodology described previously, for the velocity fields, Dirichlet boundary conditions using the PIV data are enforced and expanded into a two-cell deep perimetric zone inside the computational domain.For the calculation of the pressure fields, zero Dirichlet boundary conditions are applied on the boundaries not adjacent to a solid boundary.Zero gradient boundary conditions are enforced on the right boundary of plane A, the bottom boundary of plane B, and the left boundary of plane C.
Furthermore, as described above in section 2.2.3, the fixed pressure values inside the perimetric zone are calculated using the fixed velocity values and act as Dirichlet boundary conditions for the rest of the computational domain.The bottom boundaries of the A and C planes are close to the ground, which is a solid boundary.However, on the one hand, the distance between their boundaries and the ground is around 0.1 (plane (C)) to 0.15 (plane (A)) times the cube height; on the other hand, the distances between the boundaries that neighbor the cube are 2 to 3 times smaller.Thus, for the pressure fields, zero-valued Dirichlet boundary conditions are applied on the bottom boundaries of planes A and C.

Metrics
Surface pressure measurements along the cube centerline are available from experiments [43,44] and will contribute to the verification of the method.Thus the pressure coefficient C p is given by where p ref is a reference pressure and p dyn is the dynamic pressure 1 2 ρU 2 cube , where U cube is the freestream velocity at cube height measured using a hot-wire system [44].For the case investigated in this work, the velocity at cube height is U cube ≈ 3.2 m/s, and the corresponding dynamic pressure is p dyn ≈ 6.2 Pa.Since the reference pressure used in the experiments was well above the boundary layer [44], and consequently, the cube height of h cube = 0.11 m, it is outside the PIV planes where the method is applied and so the selection of the reference pressure in the numerical cases requires special care.The flow is solved in the PIV measurement planes, where information is available, and consistency with the freestream needs to be determined.Thus, for the numerical solution, the reference pressure was chosen at points inside the planes where the flow can be considered as close as possible to the freestream.This means that the pressure coefficient fields within each figure are with respect to an in-plane reference pressure, which was chosen at the same point for all DA implementations.This is common practice with PIV derived pressure fields where a pressure measurement point accompanies the PIV measurements [33,35] in order to alleviate a 'shift/bias' error that would otherwise arise with respect to the experimental measurements.Specifically, for planes A and B, the reference pressure is at a point near the top-left corner of each plane.The case of plane C poses an additional challenge, as the flow cannot be considered freestream-like at any point of the plane.Thus, an approach similar to [39] is considered by using the reference from plane B in combination with the average pressure difference between the overlapping points of the two planes.
Additionally, for the validation of the quality of the PIV measurements and the CFD results, some metrics are defined that quantify the velocity field divergence or the failure to obey the mass conservation Equation and the change of the fields using the SIMPLE algorithm.Thus, the normalized local continuity residual CRN is given by where ⟨∥u PIV ∥⟩ is the mean magnitude of the PIV velocity field for each measurement plane.Furthermore, the relative local velocity calculation error is defined as The normalization is performed using the mean magnitude of the velocity field within a given plane to avoid obtaining extreme values near the stagnation point where field values are small.The normalized RMSE of the velocity magnitude is used to quantify the effects of the applied methods on the results and is defined as where N cells is the number of the grid cells.Additionally, the mean value of the normalized continuity residual field µ CRN and the mean value of the relative velocity error field µ ϵ r,Velocity are used.All the methods above are implemented in the OpenFOAM solver described in the previous sections.

Results
At first, the SIMPLE algorithm (algorithm 1) and its variants using the nudging method (14a) and Kalman filtering methods (alorithms 3 and 4) are applied.In figure 4 the velocity magnitude contours calculated using the SIMPLE based method are presented.The CFD produces results, shown in figures 4(d)-(f), that have similar flow patterns compared to the PIV data, shown in figures 4(a)-(c).For example, the stagnation point at y ≈ 0.7h cube for plane A, as shown in figures 4(a) and (d), or the separation point for plane B, as shown in figures 4(b) and (e) appear in the same positions for both the CFD results and PIV data.However, for plane A, the effect of the boundary conditions is illustrated as different gradients are observed in the half-plane above the stagnation point.In addition, discrepancies are observed in the bottom half of the upstream cube wall in plane C.
As shown in figure 4, the KF, did not affect the resulting contour plots in a visually perceivable way.The same was found for the other DA methods, hence their impact is presented only in the form of the pressure coefficient C p distributions around the cube and shown in figure 5.The distributions produced by the SIMPLE based method follow the experimental data, even though the experimental data refer to surface pressure measurements and the CFD results are produced from the plane boundaries, which are at a distance from the cube walls.However, due to the use of the zero-valued Dirichlet boundary conditions in the bottom boundary of plane A in combination with the stagnation present in the right boundary, a steep gradient appears near the bottom boundary, inside the fixed velocity zone, as shown by the horizontal dotted black line in figure 5(a).
The use of DA methods did not change the shape of the C p distributions but produced a minor shift in the curves when the classic KF formulation is employed.Specifically, both nudging implementation results are almost identical, while the integrated KF formulation produces results similar to the nudging ones.However, the classic KF formulation produced a shifted C p distribution curve that also has a different curvature for plane A, as shown in figure 5(a).
Furthermore, solving only the pressure Poisson equation produces results that have visible discrepancies with the SIMPLE results and the experimental data, most notably for plane C in figure 5(c).In order to better indicate the discrepancies between the Poisson equation and the SIMPLE based methods, the pressure coefficient contours are shown in figure 6.The Poisson equation produces results that are not only offset from the ones produced by SIMPLE, as is the case for plane A in figures 6(a) and (d) but also describe different physical behavior as in plane C, shown in figures 6(c) and (f).On the contrary, for plane B, the results shown in figures 6(b) (Poisson equation) and (e) (SIMPLE) are similar.This behavior emanates from the validity of (a) the zerodivergence assumption for the PIV velocity data and (b) the 2-dimensionality assumption for the PIV velocity field.These assumptions are investigated in the next section.

Error distributions
As stated above, the assumption of zero-divergence velocity fields is a crucial one in the framework of pressure calculation from PIV data.In figure 7, the normalized continuity residual fields (CRN) are shown.The PIV data of plane B, shown in figure 7(b) have lower continuity residuals compared to planes A (figure 7(a)) and C (figure 7(c)).The regions where the zerodivergence assumption is inadequate are the ones close to the stagnation point in plane A and the recirculation zones on the left side of plane C.These discontinuities may arise due to poor image resolution and laser reflections from the solid Plexiglas cube used in the experimental study [44].
Another factor that influences the quality of the solution one can obtain is the existence of a non-zero mean velocity component that is normal to the PIV measurement plane (w).Since the measurement planes are located along the cube centerline, this component, (w), should be zero.Non-zero values are expected to stem from measurement errors, such as the non-zero laser thickness and laser sheet positioning errors.In figure 8, the influence of the 3-dimensionality of the PIV data is showcased.It is obvious that the w component is negligible compared to the magnitude of the 2D components, except for certain regions of high shear near the domain boundaries (upstream roof leading edge and leeward wall in the recirculation zone).These regions are expected to be more sensitive to measurement errors due to large spatial gradients but, as they are boundary values, they are not altered during the numerical solution of the field and thus the nearby velocity field is not altered in a non-physical way.Thus, it is clear that some form of velocity correction must be used for the pressure calculation from PIV data that present non-negligible mass discontinuities.
Furthermore, high CRN values are present near the field boundaries of the CFD solution, shown in figures 7(d)-(f).This behavior emanates from the usage of the 2-cell deep zones of fixed velocity.The velocity values inside the cells of these zones are neither corrected nor assimilated during the numerical solution.Accordingly, the CRN values on the outermost cells of these zones, are equal to the CRN values of the PIV data (shown in figures 7(a)-(c)), and the observed high CRN absolute values of the cells that lie on the interface of the perimetric fixed velocity zones and the actual computational domain are produced by the contribution to the divergence calculation of adjacent cells where velocity is either corrected or fixed.
In addition to the CRN contours, histograms of the CRN distributions are shown in figure 9. On the horizontal axis, the bins of the CRN field values, excluding the perimetric zones of fixed velocity, are presented and on the vertical axis, the relevant volume fraction values (VF), which are defined as the fraction of the total volume (area in 2D) of grid cells, in which the CRN values belong to the respective horizontal axis CRN bin These clearly indicate how the velocity correction achieved through the application of the SIMPLE algorithm has  significantly reduced the fraction of the solution space that exhibits a non-negligible continuity residual, compared to the raw PIV measurements.
As described in the previous section, the effect of DA methods cannot be easily evaluated visually using contour plots or histograms.Thus, the mean values and the standard deviations    of the methods' performance metrics are used.In tables 1-3, the effect of the SIMPLE based method and attempts towards its enhancement with DA methods is illustrated.The SIMPLE based method significantly reduces the mean (by 81% to 92%) and standard deviation (by 28% to 55%) of the CRN fields compared to the original PIV data.On the other hand, the effects of the DA methods on the performance metrics are minor, although they do reduce the mean values and the variance of both CRN and relative velocity errors, especially the classic KF formulation.As expected, the identity nudging method does not improve the results, while the BLUE nudging method manages to have an impact similar to the integrated KF formulation.This happens as using the integrated formulation, the convergence of the KF leads to adding a term to the momentum corrector, similar to the one added by the BLUE nudging method, but with a different gain matrix.
In summary, the best performing DA method (i.e. the classic formulation of the KF) reduces the extent to which the original PIV data are altered, specifically reducing the mean of the velocity error up to 3.5% and its standard deviation up to 6.5%, without adversely affecting the mass continuity of the resulting field.It is notable that it also reduces both the mean (ranging from 0% to 25%) and standard deviation (ranging from 3% to 5%) of the CRN fields.Although this effect on accuracy is minor, there is also an impact of DA methods on the SIMPLE based method's convergence, which will be illustrated in 3.6.

Grid refinement study
In CFD, mesh quality is one of the most important factors influencing the final solution's quality.Furthermore, as seen   in (21), mesh quality and size affect the errors in the calculations.The grids used in this application are chosen to comply with the spatial resolution of the PIV measurements and are structured and uniform (and orthogonal).However, they might be considered to be relatively coarse compared to a typical CFD mesh.Thus, the initial CFD grids are refined, and the data are mapped from the coarse to the finer grids.Another motivation behind this endeavour is to investigate whether the denser numerical mesh can take advantage of the coarser original data to better capture the flow behavior.Specifically, each cell of the uniform and orthogonal grid is split into four similar cells.The PIV data are mapped again onto the refined grids using bi-linear interpolation.Thus two finer grids (×4) and (×16) are produced for each studied case, as shown in figure 10.Then, after the refinement process, the flow is solved using the SIMPLE based method on each refined grid and compared to the results obtained from the nominal (initial) grids.
In table 4, the grids produced from the grid refinement procedure are summarized.As a result, all planes have similar cell numbers, and a comparison of the results is valid.
In tables 5-7, the effect of the grid refinement process on the performance metrics is shown.Both the RMSE and relative velocity error mean value ( µ ϵ r,Velocity ) and standard deviation ( σ ϵ r,Velocity ) do not change significantly for finer grids.This behavior is expected, as the data are interpolated into the refined grids and thus no additional information is available that would produce different flow behavior and consequently different RMSE and relative velocity errors.On the other hand, both the CRN mean value and standard deviation are significantly reduced as the grid is refined.Specifically, both the mean value and the standard deviation are reduced by 50% per refinement level.This improvement is observed, as the effect of the fixed perimeter zone boundary is reduced when the relative number of cells where velocity is fixed is reduced.
In figure 11, the pressure coefficient distribution is shown for different gird sizes.For all planes, the increase in grid size produces a minor offset in the pressure distributions with the most important impact observed in the regions with high gradients, especially for plane B, as shown in figure 11(b).This behavior stems from the fact that the grid refinement results in better resolution of the flow near the boundaries where zero Neumann boundary conditions are applied.
For brevity, only plane A contours are used in order to illustrate the effect of grid refinement on the resulting fields.It is evident from figures 12 and 13, that the grid refinement does not result in substantial differences in the final velocity and pressure coefficient fields, except perhaps for subtle differences near the boundaries.The contours of relative velocity errors, shown in figure 14, have similar behavior.The most noteworthy difference is observed in the CRN behavior in figure 15.Due to the grid refinement, the fixed velocity boundary zones are reduced as a percentage of the total cell number as the refinement quadruples the number of cells, while the number of fixed boundary cells doubles.Thus, the discontinuity that is observed in figure 15(a) at the fixed velocity zone is significantly reduced as the grid size is increased by a factor of 4, shown in figure 15(b) and a factor of 16 shown in figure 15(c).
Finally, it is also important to note that by refining the CFD grid, the PIV data can still be used, the advantages of the CFD solution are available, and the overall error is reduced.

Convergence and computational cost
An important aspect of the performance assessment of numerical methods is their computational cost.Apart from the  baseline SIMPLE based method, the DA methods are expected to increase the demands on computational resources per iteration.The nudging based methods are cheaper by design compared to the KF methods, as the former requires, at worst, three additional matrices (K, P, R) stored in memory and one multiplication, calculating the innovation term of (14b), to be performed during each iteration.In contrast, the KF method requires, at best, four additional matrices (K, P, R, M) stored in memory, four multiplications (the last four steps of algorithm 3), and one matrix inversion (the gain matrix update step in algorithms 3 and 4) to be performed during each iteration.For this reason, a performance scalability study is performed.The developed OpenFOAM solver is compiled using v2206 and executed in parallel using 2 cores in a shared memory system with an Intel i7-7500U processor and 8 GB of RAM.The mesh is decomposed using the Scotch library [60] and the parallel simulation is performed using Open MPI [61].Finally, the performance statistics are gathered using version 1.9 of the GNU Time utility.
The effect of the additional matrix storage and the multiplications does not significantly impact the memory and processor requirements as shown in table 8. Additionally, one would expect that the required matrix inversion in the DA methods (( 13) and ( 18)), would lead to a big difference in code performance, but its effect is not that significant as shown in table 8.This occurs because the matrix that has to be inverted in this implementation is diagonal, and its inverse can be found in linear time (O(n)).The main performance hit comes from the solution of the two linear systems defined by ( 3) and (4) using iterative methods for each SIMPLE iteration, which scale at best as O(n 2 ) [48], and have to be used for all methods.Furthermore, the classic KF method requires fewer iterations to converge compared to the SIMPLE method.
In table 9, the required computational resources are shown for both the SIMPLE based method and the KF formulations for increasing grid sizes.As expected, the KF methods require more RAM, and this requirement increases with the grid size.However, it is worth noting that the KF methods require fewer iterations to converge and consequently require less overall computational time and cost compared to the base SIMPLE method.

Conclusions
In this work, a method for pressure calculation from PIV data using the SIMPLE algorithm has been presented and supplemented with Data Assimilation methods.The DA methods were applied with the intention to improve the solution resulting from applying the SIMPLE algorithm, which corrects and then adds to the data of the velocity measurements, providing the pressure field.The developed solver directly included the PIV-based Reynolds stress tensor into the RANS equations.Regarding the DA methods, the nudging method was implemented, and two options for its gain matrix were investigated (identity and BLUE).Furthermore, the KF algorithm was coupled with the SIMPLE algorithm in two different formulations (classic and integrated).The methods were applied using two-dimensional PIV data and surface pressure measurements for three different regions along the center plane of a surface mounted cube immersed in a turbulent boundary layer flow.
The proposed SIMPLE based method was validated using experimental surface pressure measurements, and its effect in significantly reducing the continuity residuals of the original PIV data was illustrated.Specifically, the reduction of mean continuity residual values ranges from 81% to 92% and the reduction of its standard deviation ranges from 28% to 55%.
Applying Data Assimilation methods introduced a minor improvement in the resulting velocity and pressure fields.The best performing classic KF formulation reduced the mean of the velocity error up to 3.5% and its standard deviation up to 6.5%, while reducing both the mean (ranging from 0% to 25%) and standard deviation (ranging from 3% to 5%) of the continuity residuals.However, these methods did improve the convergence of the SIMPLE based method and thus came at a low computational cost.
Finally, a grid refinement study was performed to investigate the discrepancies that the PIV spatial discretization introduces to the availability of information and the subsequent performance of the SIMPLE based method.The PIV datavelocity and Reynolds stresses-were mapped onto finer grids by simple interpolations.There was no significant effect on the final results except for reducing the effect of the boundary conditions in the near boundary regions, and consequently reducing the CRN mean and variance.However, computational cost and the required iterations for convergence increased.It is noted that the turbulence effects were directly included, via the measured Reynolds stresses, in the solution of the flow field equations.An interesting concept would be to investigate a (sub-grid) model to incorporate the effect of turbulence at scales smaller than the PIV spatial and temporal discretization.
The positive outcome of the present endeavour suggests that assimilation of PIV data into numerical solutions may be a viable approach as a method to improve velocity field data and as a pressure calculation method in more complex geometries and/or configurations, where optical access problems may affect PIV data availability.Here, integration of the data, including the measured Reynolds stresses, in the iterative procedure provided the largest benefit but simple Data Assimilation methods indicated an incremental improvement in accuracy and also helped with convergence, at low computational cost.Supplementing data from other sources of measurement, e.g.surface pressure measurements and/or including sparse data in the computational simulation, might further alleviate PIV data quality problems while accounting for their uncertainties.This work can be used as a foundation for more complex ensemble based methods, in order to identify the sensitivity of the model to the direct inclusion of the Reynolds stress tensor, to improve the understanding of the model error covariance matrix, and possibly yield further improvements in accuracy.

Figure 1 .
Figure 1.PIV nodes (red) & CFD grid and cell centers (blue): (a) PIV nodes on a CFD grid (dashed lines) with representative CFD cell centers (blue) surrounding a single PIV node, and (b) representative CFD cell centers (open blue nodes) surrounded by PIV nodes (red).

Figure 2 .
Figure 2. Fixed perimeter cell zone for the top boundary of a grid, as noted in [39] for the: (a) PIV nodes & (b) CFD grid and cell centers.

Algorithm 3 .
Pressure calculation using the SIMPLE algorithm and the classic Kalman Filter formulation.Initial guess of p, u = u PIV and P = Q Compute Reynolds stresses from PIV T = T PIV while not converged do u * ← Momentum predictor: Solve discretized momentum equation, (3) Fix boundary zones as in figure 2 p * ← Pressure predictor: Solve pressure equation, (4) p k+1 ← Under-relax pressure field, (5) Correct volumetric fluxes, (6) u * k+1 ← Momentum corrector: Solve velocity correction equation, (7) Update forecast error covariance matrix:

Figure 3 .
Figure 3. Planes along the cube center plane where PIV measurements are available.

Figure 4 .
Figure 4. Velocity magnitude contours.The PIV data are displayed in the top row for planes: (a) A, (b) B, and (c) C. On the middle row, the CFD results using the SIMPLE based method (algorithm 1) are displayed for planes: (d) A, (e) B, and (f) C. On the bottom row, the CFD results using the SIMPLE based method coupled with the classic KF formulation are displayed for planes: (g) A, (h) B, and (i) C.

Figure 5 .
Figure 5. Pressure coefficient distribution along the cube centerline for different DA methods for planes: (a) A, (b) B, and (c) C. The horizontal dotted black lines indicate the boundary of the domain where the NS equations are solved for planes A & C.

Figure 6 .
Figure 6.Pressure coefficient contours.The Poisson equation results are displayed on the top row for planes: (a) A, (b) B, and (c) C. On the bottom row, the CFD results using the SIMPLE based method are displayed for planes: (d) A, (e) B, and (f) C.

Figure 7 .
Figure 7. Normalized continuity residual contours.The CRN contours of the initial PIV data are displayed in the top row for planes: (a) A, (b) B, and (c) C. On the bottom row, the CRN of the converged CFD results using the SIMPLE based method are displayed for planes: (d) A, (e) B, and (f) C.

Figure 8 .
Figure 8.The influence of the 2D velocity field assumption on the CRN of PIV data.The w PIV velocity components normalized by the respective mean PIV velocity magnitude are displayed for planes: (a) A, (b) B, and (c) C. On the bottom row, the CRN of the PIV data are displayed for planes: (d) A, (e) B, and (f) C.

Figure 9 .
Figure 9. Distribution of normalized continuity residual fields.The CRN distributions of the initial PIV data are displayed in the top row for planes: (a) A, (b) B, and (c) C. On the bottom row, the CRN distribution of the converged CFD results using the SIMPLE based method are displayed for planes: (d) A, (e) B, and (f) C.

Figure 10 .
Figure 10.Illustration of the grid refinement strategy for a reference cell.

Figure 11 .
Figure 11.Pressure coefficient distribution along the cube centerline for different grid sizes for planes: (a) A (a), (b) B, and (c) C. The horizontal dotted black lines indicate the boundary of the domain where the NS equations are solved for the nominal grids of planes A & C.

Figure 12 .
Figure 12.Velocity magnitude contours for plane A. Results are shown for: (a) the nominal grid, (b) the ×4 grid, and (c) the ×16 grid.

Figure 13 .
Figure 13.Pressure coefficient contours for plane A. Results are shown for: (a) the nominal grid, (b) the ×4 grid, and (c) the ×16 grid.

Figure 14 .
Figure 14.Relative velocity error contours for plane A. Results are shown for: (a) the nominal grid, (b) the ×4 grid, and (c) the ×16 grid.

Figure 15 .
Figure 15.Normalized continuity residual contours for plane A. Results are shown for: (a) the nominal grid, (b) the ×4 grid, and (c) the ×16 grid.

Table 1 .
Data assimilation methods performance for plane A. Bold values denote the minimum value of each column.

Table 2 .
Data assimilation methods performance for plane B. Bold values denote the minimum value of each column.

Table 3 .
Data assimilation methods performance for plane C. Bold values denote the minimum value of each column.

Table 4 .
Number of cells for the grids used in the grid refinement study.

Table 5 .
Grid refinement effect on the method performance metrics for plane A.

Table 6 .
Grid refinement effect on the method performance metrics for plane B.

Table 7 .
Grid refinement effect on the method performance metrics for plane C.

Table 8 .
Computational performance of different DA methods for plane A using the nominal grid.

Table 9 .
Performance scalability study for plane A.