Pressure reconstruction of a planar turbulent flow field within a multiply connected domain with arbitrary boundary shapes

Over the past two decades, it has been demonstrated that the instantaneous spatial pressure distribution in a turbulent flow field can be reconstructed from the pressure gradient field non-intrusively measured by particle image velocimetry (PIV). Representative pressure reconstruction methods include the omni-directional integration,^1–9 the Poisson equation approach,^10,11 the least-square methods,^12,13 and, most recently, data assimilation methods,^14,15 which are essentially flow field computation with least-square optimization to reconcile the differences between the computational results and the experimental data, and using the Poisson equation to reconstruct the pressure. Please note, as demonstrated by Wang et al.,¹⁶ the least-square approach^12,13 is mathematically equivalent to the Poisson equation approach with Neumann boundary conditions.

Most of these previous pressure reconstruction examples, however, were applied to simply connected domains.¹⁷ None of these previous studies have discussed how to apply the pressure reconstruction procedures to multiply connected domains.¹⁷ However, the pressure reconstruction capability in multiply connected domains is critical for cavitation, multiphase flow, and porous medium flow research and applications. In view of this need, and to fill the gap, we want to acquire and establish the pressure reconstruction capability for cavitation and multiphase flow measurements. To achieve this goal, we first introduce procedures to implement the parallel-ray omni-directional integration algorithm for pressure reconstruction in multiply connected domains with arbitrary inner and outer boundary shapes, and then validate the pressure reconstruction code using two test cases, respectively, including the Johns Hopkins Turbulence Database (JHTBD) on isotropic turbulence DNS (Direct Numerical Simulation) data¹⁸ and the Caltech DNS data on shock-induced non-spherical bubble collapse.¹⁹ In addition, the performance of the parallel-ray omni-directional integration code is also compared with the conventional Poisson equation approach by using the JHTBD isotropic turbulence pressure gradient field with embedded random error distributions.

The pressure reconstruction algorithm used in the current study is the rotating parallel-ray omni-directional integration algorithm,⁶ which is an upgrade from its predecessor, the virtual boundary omni-directional integration algorithm.^1–3 As demonstrated in Liu and Moreto⁷ based on simply connected flow domains, omni-directional integration methods outperform the conventional Poisson equation approach in pressure reconstruction accuracy, with the parallel ray being the algorithm with the best performance in accuracy among all omni-directional integration methods.

The introduction of the omni-directional integration was motivated by the need of minimizing the influence of the pressure gradient measurement error on the final reconstructed pressure result. The design of the omni-directional integration algorithm was based on the observation of the scalar potential properties according to the field theory, which states that pressure is a scalar potential and its gradient forms a conservative field. Thus, spatial integration of the scalar potential gradient is independent of the integration path. With this understanding, Liu and Katz^1,2 introduced the original omni-directional integration algorithm, which first integrates the measured pressure gradient from all directions toward the point of interest, and then further averaging of the pressure values obtained along different integration paths allows the minimization of the effect of the errors embedded in the measured pressure gradient data so as to obtain the final reconstructed pressure.

Figure 1 shows the schematic of the process of pressure reconstruction from the measured pressure gradient in a simply connected domain by using the omni-directional integration algorithm. This method utilizes parallel rays as guidance paths for line integration of pressure gradient. By rotating the parallel rays, effectively omni-directional paths with equal weights coming from all directions toward the point of interest at any location within the computation domain can be generated. To reduce the effect of erroneous pressure along the boundary, the parallel ray omni-directional integration algorithm utilizes iterations of boundary pressure values in the integration process. Initially, integration along the perimeter of the measurement domain gives the initial boundary pressure values, but since each boundary node is crossed by numerous paths, the initial boundary values can be corrected/replaced with the results of the omni-directional integration. The integration process is then repeated using the new boundary values, and the iterations proceed until the results converge within an acceptable threshold level, thus achieving accurate Dirichlet conditions for boundary pressure values. The convergence of the boundary pressure calculation follows an exponential error decay fashion, as demonstrated theoretically in Liu and Moreto.⁷ Once the accurate boundary pressure values are achieved, the pressure values over the entire domain can then be determined by invoking the final round of the omni-directional integration over the entire domain. Please note, according to Liu and Moreto,⁷ the performance of the rotating parallel ray depends on the rotation angle increment Δα and the distance between rays Δd^* (normalized by the grid spacing). Consistent with Liu and Moreto,⁷ throughout the work presented in this paper, all pressure reconstruction results based on the parallel ray code are obtained with the parallel ray parameters of Δd^* = 0.4 and Δα = 0.20°.

The aforementioned key features of the rotating parallel ray omni-directional integration are preserved during the implementation of the same algorithm in a multiply connected domain with arbitrary domain boundary shapes. The procedures of the implementation are illustrated in Fig. 2. A bold parallel ray serving as guidance for the line integration penetrates the planar domain at the boundary point A. Along the way, it encounters the inner void boundaries at points B, C, D, and E and finally leaves the entire domain at the boundary point F. To be consistent with the zig-zag line integration pattern (i.e., no line integration along the diagonal direction between neighboring grids is involved in the pressure gradient integration)^1–3,5–7 in the omni-directional integration implementation, here, the boundary points are determined by following the four-connected criterion;²⁰ according to which, if any of the four neighboring nodal points (left, right, upper, and lower neighboring nodal points) is in a void region, the current domain point is then regarded as a boundary point. Pressure values at all the boundary points identified in this way, including both the inner and the outer domain boundary points, will then go through the boundary value iteration process, similar to that outlined in the procedures described in Liu and Moreto.⁷ After achieving convergence, the pressure at inner nodal points within the multiply connected domain will then be calculated. Please note, during the entire pressure reconstruction process, nodal points within the void areas are not involved in any form of calculations. While preserving the nature and advantage of the parallel ray omni-directional pressure reconstruction at places with flow data, the new implementation of the algorithm is capable of processing an arbitrary number of inner void areas with arbitrary boundary shapes, as shown in the following validation examples.

Validation of the multiply connected domain pressure reconstruction code is first conducted using a fully resolved time-varying isotropic turbulence database available at the JHTDB.¹⁸ The dataset is generated by direct numerical simulation of forced isotropic turbulence,^21,22 which was performed on a three-dimensional periodic grid with 1024 ×1024 ×1024 nodal points using a pseudo-spectral method for solving the Navier–Stokes equations. Fluctuation energy of the isotropic flow field at low wave number is injected at each step of simulation to maintain steady state conditions. The database provides 1024 instantaneous realizations of the isotropic turbulence flow in a 2π × 2π × 2π three components of the velocity and the pressure fully resolved, both temporally and spatially. Figure 3 shows the pressure distribution over a sample planar area of 256 ×256 grids selected from the dataset. Also shown in the figure is the pressure gradient components obtained by central finite difference from sample area. We further allocated arbitrary void areas over the pressure gradient fields and then use the parallel ray omni-directional integration code to reconstruct the pressure. Comparison of the reconstructed pressure with the DNS pressure is then made to characterize the error in the reconstructed pressure. Results of this validation are shown in Fig. 4, where the first column shows the reconstructed pressure with different inner void patterns, the second column shows the distribution of the difference between the reconstructed and the original DNS pressure nondimensionalized by the standard deviation of the DNS pressure, and the third column shows the corresponding probability density function (pdf) of the error. From Fig. 4, it can be seen that the error distributions of the reconstructed pressure are fairly consistent with various void patterns and arbitrary boundary shapes, involving sharp concave/convex corners on both inner and outer boundaries, as shown in the top three rows of plots in Fig. 4, with the characteristic standard deviation of the errors varying from 0.55% to 0.61% of the characteristic turbulence pressure fluctuation magnitude, i.e., the standard deviation of the turbulence pressure fluctuation. The significant increase in the standard deviation of the errors (increased to 14.8%) as well as the secondary peak on the pdf plot for the letter “SDSU” void pattern shown in the last row of plots in Fig. 4 is due to the isolated region within the letter “D”, according to an analysis (not shown here due to the length consideration). In that analysis, by excluding the isolated region within the letter “D” in pressure reconstruction, it turns out that the error statistics remains the same as other void pattern cases discussed before for Fig. 4.

To further evaluate the performance of the code, random noise was homogeneously superimposed onto the entire domain of the DNS isotropic turbulence pressure gradient field. The noise amplitude is chosen as 40% of the maximum DNS pressure gradient $∇pDNSmax$⁠. A total of 1000 statistically independent realizations of noise distribution were generated in this manner. We then further superimpose different void shapes onto those 1000 realizations of the noise-imbedded pressure gradient fields, with one sample realization shown in the first two columns in Fig. 5. Furthermore, we reconstruct the pressure (the third column of plots from left in Fig. 5) using the parallel ray omni-directional integration code and then compare the differences between the reconstructed and the original DNS pressure to get the error distributions (as shown in the fourth column of plots from left in Fig. 5). In addition, to compare the performance in the accuracy of the reconstructed pressure, we also reconstructed the pressure from the 1000 realizations of the error-embedded isotropic turbulence pressure gradient fields using the conventional Poisson equation approach with Neumann boundary condition and the Dirichlet boundary condition. Please note the Dirichlet boundary condition is provided by the converged boundary pressure values from the parallel ray omni-directional integration. Comparisons of the cumulative average errors of the reconstructed pressure based on those 1000 realizations of the error-embedded isotropic turbulence pressure gradient fields using the parallel ray and the Poisson equation methods are presented in Fig. 6. The statistics of those pressure reconstruction results are shown in Tables I–III. The cumulative average error of the reconstructed pressure is based on the error $ϵij$ in the reconstructed pressure $p̃ij$ computed as the difference between the reconstructed pressure $p̃ij$ and the exact (i.e., DNS) pressure $pij$⁠,

with its standard deviation defined as

where (and throughout the paper) the overbar represents the spatial averaging over the entire domain for a specific reconstructed pressure realization.

As can be seen from Fig. 6, convergence of the error statistics is achieved for most test cases after 500 realizations. It can also be seen that, as shown in Fig. 6(a), when void areas are introduced, the error of the reconstructed pressure in the multiply connected domain also increases in comparison with the simply connected domain for the reconstructed pressure obtained by parallel ray omni-directional integration code.

Comparing Figs. 6(a) and 6(b), it can be seen that, even the smallest error of the Poisson equation Neumann boundary condition approach, which takes a value of 0.57 as seen on Fig. 6(b), is significantly greater than the highest error of the parallel ray approach, which only takes a value of 0.27. This is a clear indication that the rotating parallel ray omni-directional integration is capable of providing a better reconstructed pressure accuracy than the Poisson equation Neumann boundary condition approach when dealing with noise-embedded pressure gradient.

Comparison of the cumulative error statistics between the rotating parallel ray [Fig. 6(a)] and the Poisson equation approach with Dirichlet condition [Fig. 6(c), where boundary pressure values are taken from the parallel ray calculation results] shows that the reconstructed pressure obtained by the Poisson equation Dirichlet is still not as good as that obtained by the parallel ray approach alone. This also can be seen by comparing the first columns of Tables I and III. Clearly, we can see that the rotating parallel ray omni-directional integration outperforms the Poisson equation approach.

For further validation, we also applied the parallel-ray omni-directional integration code to reconstruct the instantaneous pressure in a shock-induced bubble collapse in water based on the Caltech DNS data for that transient process.¹⁹ Basically, the Caltech DNS data document the evolution process in which an isolated spherical bubble with an initial diameter of 100 μm immersed in water with a pre-shock pressure of 1.0 atm is suddenly subjected to a shock pressure $ps$ with a magnitude of 353.0 atm swept from left to right in the field. Impacted by the shock, the bubble undergoes a quick collapse process. When the bubble reaches its smallest volume due to the effect of shock compression, it starts a deformed expansion, along which, a quick water jet is formed and is penetrating through the bubble. Figure 7 (Multimedia view) shows two sample instants of the bubble deformation process, one is at t = 244.36 ns, when the bubble is during the shrinking process [Fig. 7(a)]. The other is at t = 351.53 ns, when the bubble volume is bounced back from its minimum value while a jet is penetrating the deformed bubble. In Fig. 7, the DNS pressure, the reconstructed pressure, the spatial distribution, and the Probability Density Function of the reconstructed pressure error (represented by the difference between the reconstructed and the original DNS pressure) are presented for those two sample instants. The pressure gradient needed for the pressure reconstruction is calculated using a central difference scheme from the DNS pressure. As can be seen from Fig. 7, the standard deviation of the reconstructed error for both cases is all very small, with $prec−pdns/psstd=5.4×10−5$ for the t = 244.36 ns case and $1.53×10−4$ for the t =351.53 ns case.

In conclusion, the Rotating Parallel Ray Omni-directional Integration algorithm has been successfully applied for the first time to the pressure reconstruction in multiply connected domains with arbitrary boundary shapes. Validation using both DNS data on isotropic turbulence and shock-induced non-spherical bubble collapse indicates that the Parallel Ray Omni-directional Integration can accurately reconstruct the pressure in multiply connected domains, even when the pressure gradient field is embedded with error. Moreover, the parallel ray method once again outperforms the Poisson equation approach, consistent with the observation made in Liu and Moreto⁷ for simply connected domains.

The successful implementation of the parallel ray pressure reconstruction method to multiply connected domains paves the way for a variety of important applications, including, for example, experimental characterization of pressure field changes during the process of cavitation bubble inception, growth and collapse, non-intrusive unsteady aerodynamic force assessment for an arbitrary body shape immersed in flows, multi-phase flow investigations, etc. In particular, as an immediate follow-up effort, the parallel ray pressure code will be used for the instantaneous pressure distribution reconstruction in a turbulent flow surrounding cavitation inception bubbles occurring on top of a cavity trailing corner based on high-speed PIV measurements.²³ 

Please note that the procedures described in this paper can be readily extended to three-dimensional (3D) applications.²³ The major difficulty in 3D applications is the significant increase in computation cost due to the addition of the third dimension to the computation domain. However, this technical difficulty can be resolved by invoking GPU (Graphics Processing Unit) computation, which has been successfully demonstrated by Wang et al.⁸ and Agarwal et al.⁹ for 3D pressure reconstruction in simply connected rectangular cuboid computation domains. In 3D pressure reconstruction applications, the rotating parallel ray omni-directional integration should be adapted to icosahedron parallel ray omni-directional integration as adopted in Wang et al.,⁸ where the normal direction of the surfaces of a densely configured virtual icosahedron surrounding the 3D computation domain serves as the orientation guidance for uniformly distributed 3D parallel rays passing through the computation domain. The method of the identification of the boundary points of multiply connected 3D domain is similar to that for the two-dimensional (2D) multiply connected domains as discussed in this paper, except that the 6-connected criterion²⁰ instead of the 4-connected should be adopted for the 3D boundary point identification so as to ensure that the 3D boundary point determination is consistent with the 3D zig-zag fashion of the spatial line integration. The convergence rate of the 3D application should be similar to that for the 2D application by following an exponential error decay fashion with respect to the times of the iteration, which has been predicted theoretically by Liu and Moreto.⁷ As demonstrated in Liu and Moreto,⁷ a fine increment in parallel ray rotating angle is beneficial for a fast convergence rate and a high accuracy in reconstructed pressure. Similar behavior should be expected for applications of pressure reconstruction in multiply connected 3D domains, where the rotating angle for 2D is replaced by the dense spatial angel arrangement for 3D as represented by the dense icosahedron surfaces. All these will be verified in a follow-up paper in 3D applications of the technique in multiply connected domains.

This work was sponsored by the Office of Naval Research Grant No. N00014-21-1-2392 (Dr. K.-H. Kim is the Program Officer). The authors would like to thank Professor Tim Colonius for his permission of the use of his shock-induced non-spherical bubble collapse DNS data in this research. The help from Mr. Jose Rodolfo Chreim in retrieving the bubble collapse DNS data and the help from Dr. Minping Wan and Professor Charles Meneveau in retrieving the isotropic turbulence DNS data from JHTDB are all gratefully acknowledged.