Experimental validation of computational fluid dynamics for solving isothermal and incompressible viscous fluid flow

In order to validate a computational method for solving viscous fluid flows, experiments are carried out in an eccentric cylindrical cavity showing various flow formations over a range of Reynolds numbers. Especially, in numerical solution approaches for isothermal and incompressible flows, we search for simple experimental data for evaluating accuracy as well as performance of the computational method. Verification of different computational methods is arduous, and analytic solutions are only obtained for simple geometries like a channel flow. Clearly, a method is expected to predict different flow patterns within a cavity. Thus, we propose a configuration generating different flow formations depending on the Reynolds number and make the experimental results freely available in order to be used as an assessment criterion to demonstrate the reliability of a new computational approach.


Introduction
The computational study of a physical problem is accepted as reliable, if the employed implementation is validated by using the experimental results. Especially, in fluid dynamics simulations, since the computational methods are still being developed, the experimental results with clearly identified initial and boundary conditions are of paramount importance. Experimental data are used for validating the numerical results, specifically for cases where analytic solutions are not known. Ample studies are available in the literature comparing the experimental and numerical results for Newtonian fluids like water [1,2], blood [3,4], even for very thick fluids [5], during casting [6], and also for geometries with an obstacle creating vortex shedding [7][8][9] or even for a mixture [10]. Especially, in engineering science, a commercial software or an in-house code is employed for a numerical analysis of a simplified system.
Assumptions have to be made for the simplification; thus, the reliability of the used software is challenging to justify. Furthermore, we stress that different types of partial differential equations make the coupled system difficult to solve [11]. Hence, various methods and implementations exist and are emerging; examining their reliability is of great interest. In benchmarking new numerical schemes [12], the reference results are generated by using direct numerical simulations in order to design a computationally challenging configuration [13][14][15][16] for testing the robustness of the numerical implementation. For some special applications, experimental data are utilized for a specific geometry and used to validate the numerical results [17,18]. For systems with lacking the experimental results, direct numerical simulations are used for generating test cases to verify the accuracy of a numerical implementation [19,20].
Generally speaking, a numerical implementation needs to be robust and the used method is expected to deliver reliably the accurate results. A computationally challenging configuration examines the robustness of the numerical implementation. An experimental setup provides a realistic setting such that we can quantify the accuracy of the method. We propose a simple yet effective geometry generating different flow formations that are challenging for a numerical implementation to simulate reliably. Experimental data are recorded, processed, and made public in [21] to be used for testing a new or existing implementation. The setup is simple, a cavity flow of a Newtonian fluid, and we compare the computed velocities with the experimental results directly. By changing the inlet velocity, the flow pattern changes completely such that developers may find these experimental results useful for examining performance of their numerical implementations.
In this work, we collect the velocity data of a Newtonian (linear viscous) fluid in the case of an incompressible flow. We refer to [22][23][24] for studies with a similar experimental setup for different systems. Herein, we present an experimental study on an eccentric cylindrical cavity, in which the flow is studied by means of particle image velocimetry (PIV) [25]. This standard technique is often employed for laminar flows [26][27][28], in order to study the transient response [29], or behind an obstacle [30,31]. In this work, fluorescent particles are used both for flow visualization (FV) and PIV at the central plane of a symmetric flow cavity. The 2-D measurements at this plane are compared to the results obtained from 3-D computations. As long as the computational results match the experimental results on the symmetry plane, we expect the implementation being reliable. We examine the method proposed by Abali [32] based on the finite element method by using opensource libraries developed under the FEniCS project [33]. The proposed approach is already verified by exact solutions for simple geometries, and the method attains a local monotonic convergence for a three-dimensional computational domain. Herein, we validate the method by using the experimental results obtained in this work.
We propose a simple experimental setup providing different flow formations depending on the Reynolds number. This experimental data are novel for examining a computation implementation since the same configuration generates various flow formations within the laminar regime by just changing the inflow velocity. A detailed comparative analysis is undertaken between the experimental data and computational results. Computational fluid dynamics (CFD) uses often the finite volume method (FVM), although we use herein the finite element method (FEM) that allows incorporating multiphysics more efficiently.
There are mainly three approaches like projection methods [34][35][36][37], stabilization with limiting strategies [38][39][40], and using enriched or mixed formulations for solving pressure and velocity fields [41][42][43][44]. FEM is still emerging for fluid dynamics, and we propose to use the experimental results for validating new computational approaches. We demonstrate qualitative and quantitative studies and make the experimental results and numerical implementation-CAD model, used mesh, and code in FEniCS to be used under GNU public license [45]-publicly available (see Abali [21]) in order to encourage scientific exchange.

Flow loop
The experimental setup is sketched in Fig. 1, and the camera from the top records the flow within the cavity through a high-pass filter as shown in pictures.
The flow geometry is a short aspect ratio cylindrical cavity, which is cast out of clear silicone. The cavity is manufactured in two steps. First, a mold assembly is machined of four aluminum parts to form the flow cavity. A low-temperature metal alloy (Cerrolow 117, melting point 47.5 • C ) is cast into this cavity. Then, the metal model is removed from the aluminum mold assembly and placed in a silicon gel bath (Sylgard 184, Dow Corning Inc., USA), hardened after curing, and the alloy is removed from the model by melting, leaving a transparent flow model.
The cavity is of diameter D = 19 mm and length L = 25.4 mm with a right-handed reference system (xyz) as illustrated in Fig. 2. The flow in the feed and discharge tubes is in the x-direction. The laser light sheet is in the xyplane that is the symmetry plane of the geometry and the camera looks into the flow model along the z-direction.
It is fed and discharged through eccentric in-line cylindrical pipes; both of diameter d = 6.35 mm and tangent to the cylindrical wall of the main flow cavity. The reason for the choice of this particular geometry is that it provides a combination of challenging boundary conditions, even for laminar flows. For example, the two-dimensional backward facing step geometry, both unbounded and bounded forms above, have been usually taken as a reference geometry in flow simulations (see, for example [46]). The current geometry presents bounded, three-dimensional forward and backward facing steps at once.
The working fluid is isopropyl alcohol of density = 0.785 g/cm 3 and kinematic viscosity of = 0.031 cm 2 ∕s . Its index of refraction is n = 1.378 , which is near, but not a perfect match, to that of the silicone casting ( n = 1.41 ) out of which the flow model is built. The fluid is seeded with 6 μm fluorescent particles (Duke Scientific Corporation, Red Fluorescent Polymer Microspheres) for flow visualization and PIV measurements. The flow pumping system described by Tsai and Savaş [47], which is comprised of a piston and a gear pump in series, is employed here to drive the fluid in the system. In the current application, the piston pump is turned off and the servomotor of the gear pump is replaced with a microstepper motor (25 000 pulses per revolution) to achieve stable operation at very low speeds. The gear pump capacity is 10.0 cm 3 /revolution. The gear pump speed ranged from 0.005 (200 s per revolution) to 1 revolution per second (rps). Flows are identified by their Reynolds number based on the flow rate in the inlet pipe, where Q is the volumetric flow rate. At very slow rotation rates of the pump, i.e., 0.005 and 0.010 rps, the flow rate was actually determined from the PIV data, as the static leakage of the gear pump introduces large uncertainties in flow rates determined from the displacement volume of the gear pump. For the steady inflow experiments presented here, the Reynolds numbers ranged from Re d = 2.7 -650. Later on for comparison, we will use Reynolds numbers all computed from the PIV data.

Optics
The optical components include a digital camera Kodak Megaplus ES1.0 for imaging (1008 × 1018 pixel resolution), a cemented spherical-cylindrical lens pair for laser light sheet generation, and a removable flat mirror for laser beam selection. The cemented pair consists of a + 400 mm FL plano-convex lens and − 20 mm FL cylindrical lens. Light beams form the PIV lasers, NewWave Gemini-30, dual-pulsed YAG laser, wavelength = 0.532 μm , and flow visualization laser, a continuous wave (CW) argon-ion laser, operated at green wavelength = 0.5145 μm , share the same optical path into the flow cavity. The removable mirror shown in the figure is put into position when the PIV lasers are needed. The thickness of the light sheet at the flow cavity is estimated to be about 50 μm . The camera, fitted with a Canon FD 50 mm, f1:1.4 lens, a 5-mm extension tube, and a C-mount adapter, is used for both FV and PIV image acquisitions. The Stokes's shift of the fluorescent seed particles is 50 nm. Therefore, a high-pass filter with a cutoff wavelength of 0.560 μm was used to record the light emitted from the fluorescent seed particles. Experimental  Table 1. Image acquisition and processing techniques are similar to those described and used in [22,48,49]. In particular, the PIV algorithm developed in house is used for correlation velocimetry. Postprocessing is done using routines written in IDL (Interactive Data Language).

Computational method
Isothermal flow of viscous fluid is given by balance equations of mass and momentum: respectively, where x i indicate the physical coordinates with index i from 1 to 3 and we employ the usual comma notation for spatial derivatives () ,i = ()∕ x i as well as Einstein's summation convention over repeated indices. All field variables are expressed in Cartesian coordinates. The unknowns to be solved are the mass density, , and the velocity, v i , of fluid particles occupying x i in space at the instant t. The specific supply term, g i , is because of gravitational forces, and its value is known, but this term is negligible in the simulations, as we know that the flow is predominantly caused by the pressure difference. We assume that isopropyl alcohol is a Newtonian fluid and model the flow by using a linear constitutive equation for Cauchy's stress: with the material constants , , and p, the hydrostatic pressure. This model is also called Navier-Stokes equation, as it has been proposed by both scientists in a phenomenological way, we refer to [50] for historical remarks. Under the usual assumption of an incompressible flow and the fluid has been homogeneous initially, the mass density remains constant in time and space. Then, from Eq. (2), we obtain In this case, namely for an incompressible flow, the mechanical pressure − 1 3 ii is identical to the hydrostatic pressure, p; hence, p corresponds to the (mechanical) pressure generated by a pump. Since the mass density is given, we use Eq. (4) 1 as the governing equation of p. We emphasize that usually Eq. (4) 1 is inserted into Eq. (4) 2 , which we circumvent for the sake of numerical stability. Various computational implementations exist in the literature for solving the governing equations (4) resulting in velocity and pressure fields. A straight-forward solution strategy in a two-dimensional continuum by using the finite element method is described in [51,Sect. 1.7]. The strategy relies on generating a so-called weak form out of governing equations. Trying to use this weak form in the case of a three-dimensional continuum starts causing numerical instabilities [52][53][54]. For circumventing these numerical problems, a great effort has been undertaken, for a brief review, we refer to [55]. As a remedy to these numerical problems, stabilization methods have been proposed in [56][57][58], and these methods are advanced [59][60][61][62] and implemented [63][64][65]. Also different numerical methods have been suggested for a robust implementation [66][67][68][69][70]. Another approach has been taken by changing the weak formulation [71][72][73][74]. Novel approaches have been proposed for advancing meshing algorithms as well as developing better solvers for fluid dynamics [75][76][77][78]. Also novel developments exist in hybrid methods as in [79][80][81] leading to research codes as in [82] as well as incorporating semi-analytical methods [83] or by modeling jump conditions on element boundaries in [84][85][86]. All these different theoretical approaches and various numerical implementations indicate that a computational method accepted by all communities is unlikely to occur in the near future. We employ the method developed and verified by Abali [32] in order to generate a weak form for computing velocity and pressure as follows: with The boundary term, n i ij , vanishes in the case of an incompressible fluid since the prescribed pressure p against the plane outward, t i = −pn i , leads to the formulation on the boundary, n i ij = n i ( ij + p ij ) = t j + n j p = 0.
We use continuous linear elements being the standard formulation for the finite element method [87]. For each node, four variables {p, v 1 , v 2 , v 3 } are computed in threedimensional space belonging to where [H n ] 4 is a four-dimensional Hilbertian Sobolev space of class C n . We use n = 1 , i.e., linear form functions for pressure and also for velocity. We circumvent any numerical problems affected by the Ladyzhenskaya-Babuska-Brezzi condition [88]. The additional governing equation restricting the pressure gradient enables using the same space for both fields (velocity and pressure). No stabilizing methods or limiting procedures are used in the implementation such that the approach attains monotonous convergence, in other words, using a finer mesh improves the accuracy without affecting the stability. We emphasize that g, s, mm are chosen as units leading to μN = g mm/s 2 and Pa = μN/mm 2 . The preprocessing has been realized in Salome 1 by using NetGen algorithms 2 for the triangulation. The mesh is of tetrahedron elements and is transformed as explained in [51][Appendix A.3] in order to implement in a Python code using packages developed under the FEniCS project, 3 which is wrapped in C++, and run on a Linux machine using Ubuntu. 4 Figure 3 shows the flow patterns for the cases listed in Table 1. The traces of the fluorescent particles recorded at long exposure under continuous illumination show the steady flow pattern at the central plane of the cavity. Each picture is constructed by adding subsequent images. High-speed regions show long streak, whereas the lowspeed recirculation zones show dotted trajectories. All flows are steady, with the possible exception of that at Re = 650 where the corresponding video sequence shows signs of instability at the unconstrained side of the inflow jet. Blurred signatures of these instabilities are marked with arrows in Fig. 3h. At the lowest two Reynolds numbers, the flow spreads into the cavity filling it, and, then, converges to the exit port. At the higher Reynolds number, the flow attains a ballistic regime where the jet from the inlet port heads straight to the exit port with little distortion in the velocity profile. A large circulation zone is set up in the cavity driven by the ballistic jet. At intermediate Reynolds numbers (16 and, to some extent, 32), the transition from cavity filling regime to the ballistic regime is apparent.

Results and discussion
Computational solution of the same geometry has been obtained on a mesh as in Fig. 4 with 53,208 degrees of freedom, which is a coarse mesh for such a study, selected on purpose, demonstrating the robustness of the implementation.
Then, we used a fine mesh as depicted in Fig. 5 resulting in 345,972 degrees of freedom and obtained very similar results; we present herein the results obtained with the coarse mesh.
We stress that the local monotonous convergence as demonstrated in [32] assures an improved accuracy for refined mesh choices. Because of this monotonous convergence, verification of the numerical results by an error quantification as in [89] is not necessary since we validate the computational results by comparing them to the experimental results. An important challenge is to maintain the incompressibility; hence, we use the ratio of mass flow at the inlet per outlet as a measure to quantify the accuracy of the incompressibility condition. Space and time discretizations affect this measure, by choosing the Fig. 3 Short time-averaged images for the cases compiled in Table 1. The arrows in (h) show blurred signatures of instabilities developing at the free edge of the jet aforementioned degrees of freedom with a time step of 1 ms, we have obtained an approximate error of 0.25%. The simulation lasts nearly 30 h on the coarse mesh and 600 h on the fine mesh by using 40 cores 5 in parallel by utilizing mpirun using the default backend PETSc [90] with the mesh partitioner SCOTCH [91] and by sharing vertices in separating the mesh to threads. Especially, the solver parameters are the key for a successful computation in incompressible flow simulations; therefore, we emphasize that the code is made public as well in [21]. A standard Newton-Raphson solver with consistent linearization is used. Because of using small time steps and the computed values for initial values for the algorithm, less than five iterations have been observed to reach a tolerance of 10 −3 in the aforementioned units, and tolerance means the Euclidian norm of changes in each iteration. For each iteration, Krylov subspace method is used for the linearized problem, and we have selected the generalized minimal residual method gmres for the underlying non-symmetric operator. We have used an algebraic multigrid amg preconditioner [92] as implemented in FEniCS by using PETSc algorithms. Specifically, we use an unsmoothed aggregation. The unrestarted gmres with a tolerance of 10 −4 has converged for coarse and fine meshes roughly the same manner such that the preconditioner amg bounds the condition number of the problem independent on the mesh size. The success of the preconditioner depends on the chosen boundary conditions [93]. We emphasize that on all boundaries, Dirichlet conditions of either unknowns are used. The experimental setup is modeled by setting the velocity profile at the inlet of radius r = d∕2 along the inlet pipe axis, providing a parabolic velocity distribution directed inward the cylindrical cavity. The maximum velocity, v m , at the center is given by, using the Reynolds number, Re, up to 600 linearly increased in 600 s under the assumption that the change is slow enough to attain the steady-state solution for each time step. We stress that the computation is transient and we compare the results by choosing the instants depending on the maximum velocity obtained at the same location from the PIV data and computation. Specifically, at the symmetry plane, z = 0 , we have chosen the center of inlet pipe, x∕L = −0.1 and y∕L = 0.625 such that the velocity is directed along the x-axis. The time values chosen for the comparison read 3 s, 6 s, 17 s, 34 s, 65 s, 145 s, 290 s, and 550 s, and we present in Table 2 maximum velocities, On all walls, the velocity is assumed to vanish-simply the no-slip and no-penetration condition is applied by setting zero velocities-throughout the simulation. At the outlet, the reference pressure of 1 bar is set as a Dirichlet boundary condition. Initially, zero velocity and reference pressure are defined in the whole domain. Figures 6, 7, 8, 9, 10, 11, 12, and 13 show side-by-side comparison of measurements and calculations at the symmetry plane of the flow cavity corresponding to the flows in Fig. 3.
In the figures, the experimental results are shown on left side and the numerical ones on the right side. Selected velocity profiles are shown on the top row and streamline patterns on the bottom row. For the computational postprocessing, the line integral convolution technique in ParaView has been employed [94]. Therefore, we compare streamlines on the symmetry plane where measurements were taken. Two important phenomena are captured in these simulations. First, the pressure drop because of the cavity results in a circulation beyond a threshold value between Re = 6 − 16 , and this phenomenon is visible by comparing streamlines of cases (b) and (c). Second, the circulation center moves from left to right with an increasing Reynolds number, see streamlines for cases (c) to (h). Both phenomena are predicted by simulations; hence, we conclude that, qualitatively, the agreement is adequate for all Reynolds numbers. We stress that the chosen time values are slightly changing the streamlines. Therefore, by choosing the instants we introduce an uncertainty to the comparison. We also denote that in higher Reynolds numbers even if the instant is chosen where the maximum velocity of the computation is lower than the experimental, plots show that the velocity is overestimated in computations. The exact location of the circulation center is not captured precisely, we repeat that the boundary conditions are set to increase the inlet velocity slowly so the  The pressure distribution from the computation is reasonable-despite the fact that the mesh is coarseproviding a numerical robustness difficult to obtain in an FEM computation. For a quantitative comparison of the results, we demonstrate in Fig. 15 velocity profiles as the Reynolds number increases.
Specifically, along selected lines along y at x∕L = {0.2, 0.5, 0.8} , velocities have been extracted from the PIV data and directly compared to the computational results. For postprocessing the computational results, velocity profiles are normalized by using the maximum velocity, v max , given in Table 2 in order to scale all the results in a single plot. The agreement provides a measure for reliability of the computation within the great range of Reynolds numbers presented herein. Especially, in low Reynolds numbers, the computed velocity profile is slightly higher than the experimental profile. This result relies on the chosen instant indicating a higher Reynolds number as well, see Table 2. In high Reynolds numbers, the computation overestimates the experimental results, we stress that experimental results are obtained by evaluating steady flows in the cavity. At Reynolds numbers beyond 650, the highest considered here, the flow in the cavity becomes unsteady, and turbulent at even higher Reynolds number. Even the case of Re = 650 discussed here shows occasional signs of unsteadiness. The quantitative agreement at x∕L = 0.5 demonstrates that the computational results approximate the motion of the recirculation center appropriately. However, the exact position of the recirculation center is not captured leading to a less accurate agreement at x∕L = 0.2 and x∕L = 0.8 . This issue is partly related to the chosen time values for the evaluation. Computational implementation from [32] has been used to obtain all the results. By selecting the numerical solver with the aforementioned settings, a transient computation of the nonlinear and coupled system has been performed without any numerical instabilities. The flow formation phenomena have been correctly predicted such as recirculation starting beyond a threshold velocity as well as its center moving from left to right. The exact location of the recirculation center is not detected by the computation. For judging the reliability of the numerical implementation, the overall agreement is adequate.

Closing remarks
An experimental validation of the viscous fluid flow computation has been presented in order to assess the reliability of a new approach or implementation. The setup configuration is designed for reducing the comparison to a (symmetry) plane, in which a flow pattern formation is observed in the steady state for Reynolds numbers of a sizeable range. By using this setup, validating a numerical implementation is simple; we have demonstrated such a procedure by using a numerical implementation based on the finite element method. For a qualitative analysis, streamline patterns are examined; for a quantitative verification, velocity profiles are compared along selected lines on the symmetry plane. Geometry and configuration are given in detail, and all experimental results and also its computational implementation used herein are publicly available in [21] to be used under GNU public license [45]. We hope that the data presented here may prove useful for testing new computational ideas.  Table 2