Experimental velocity data estimation for imperfect particle images using machine learning
EExperimental velocity data estimation for imperfect particle images using machinelearning
Masaki Morimoto, ∗ Kai Fukami, † and Koji Fukagata ‡ Department of Mechanical Engineering, Keio University, Yokohama 223-8522, Japan (Dated: May 5, 2020)We propose a method using supervised machine learning to estimate velocity fields from particleimages having missing regions due to experimental limitations. As a first example, a velocity fieldaround a square cylinder at Reynolds number of Re D = 300 is considered. To train machine learningmodels, we utilize artificial particle images (APIs) as the input data, which mimic the images of theparticle image velocimetry (PIV). The output data are the velocity fields, and the correct answersfor them are given by a direct numerical simulation (DNS). We examine two types of the inputdata: APIs without missing regions (i.e., full APIs) and APIs with missing regions (lacked APIs).The missing regions in the lacked APIs are assumed following the exact experimental situation inour wind tunnel setup. The velocity fields estimated from both full and lacked APIs are in greatagreement with the reference DNS data in terms of various statistical assessments. We furtherapply these machine learned models trained with the DNS data to experimental particle imagesso that their applicability to the exact experimental situation can be investigated. The velocityfields estimated by the machine learned models contain approximately 40 folds denser data thanthat with the conventional cross-correlation method. This finding suggests that we may be able toobtain finer and hidden structures of the flow field which cannot be resolved with the conventionalcross-correlation method. We also found that even the complex flow structures are hidden due tothe alignment of two square cylinders, the machine learned model is able to estimate the field inthe missing region reasonably well. The present results indicate a great potential of the proposedmachine learning based method as a new data reconstruction method for PIV. I. INTRODUCTION
High-resolution fluid data in both space and time are essential to advance our understanding on complex fluid flowphenomena. To meet this requirement, a vast range of measurement techniques have been proposed and improved todate, such as hot-wire anemometry, Schlieren photography, laser Doppler velocimetry, and particle image velocimetry(PIV). Among these, PIV has been widely used due to its capability of obtaining a velocity field in a two dimensionalplane (i.e., planar PIV) or even in a three dimensional volume (i.e., tomographic PIV). For PIV, the cross-correlationmethod [1] has been used as the major technique to estimate velocity fields from particle images. The mean velocityof the particles in a finite sized region called an interrogation window is determined by searching the region whichhas the maximum correlation of intensity distribution among two consecutive images in time. Another method isthe optical flow method proposed by Horn and Schunck [2]. In contrast to the concept of pattern matching in thecross-correlation method, the optical flow method is able to estimate velocity vectors by a spatio-temporal gradientof intensity distribution and to obtain higher resolution data than the cross-correlation method. However, one ofthe big issues here is that the error would increase significantly when the velocity magnitude increases because thetime differential term is included in the optical flow formulation. Therefore, proper experimental equipments, i.e., ahigh-speed camera and a laser emitter, are indispensable to apply the optical flow method in order to obtain particleimages in a short stride of time.Not only for the methods introduced above, estimation of velocity fields in PIV is directly related to the intensitydistribution of particle images. Therefore, the experimental images need to be as noiseless and appropriately bright-ened as possible. However, we often encounter cases where we cannot avoid regions with insufficient particle images,especially in the cases with immersed objects. For instance, the laser sheet is irregularly reflected on the surface of anobject. Also, an object having edges causes extremely white regions due to halation and shadow regions. An exampleof such a particle image around a square cylinder, taken by the authors, is shown in Fig. 1(a). Although the squarecylinder is made of transparent acrylic resin, we can observe white and dark regions below the cylinder due to itscorners. Also, it can be seen in Fig. 1(b) that the correlation coefficient is significantly low in the regions aroundand under the cylinder, which suggests that velocity estimation using the cross-correlation method is poor in thoseregions. ∗ masaki.morimoto@kflab.jp † [email protected] ‡ [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] M a y FIG. 1. (a) A typical particle image for a flow around a square cylinder made of transparent acrylic resin. Illumination isdone by a laser sheet from the top. The halation light and shadow regions significantly affect the velocity estimation. (b)Distribution of correlation coefficient. Correlation is quite low in the regions where the illumination is not appropriate.
To overcome such issues, several techniques have been proposed to date. One of the well-known techniques is topaint an object surface in flat black so that the laser light can be absorbed [3]. This can reduce the effect of irregularreflection and improve the velocity estimation. However, because the reflection cannot be prevented perfectly, it isstill difficult to acquire the particle images near the surface. Moreover, since the model is painted, information in theshaded region is totally lost. The other candidate is the use of laser induced fluorescence (LIF) [4]. In this case, tracerparticles are colored by fluorescent paint that emits light at a different frequency from that of the irradiated lasersheet, and by observing only the light emitted from particles the influence of reflection and halation can be mitigatedalmost perfectly. Although the LIF method performs very well, its use for gas flow experiments is burdensome becausethe fluorescent paint, e.g., rhodamine B, usually includes carcinogenetic material [5]. Therefore, a new method whichcan accurately estimate a velocity field without the aforementioned issues is eagerly desired.As a proof of concept of surrogates for the aforementioned conventional methods, some methods have been developedfor flow data estimation and reconstruction from limited measurements. These efforts can be seen not only inexperimental but also in computational fluid dynamics due to the high demand for high-resolution flow data. Gappyproper orthogonal decomposition (Gappy POD) proposed by Everson et al. [6] has been considered as one of thecandidates to reconstruct a flow field from incomplete data. Bui-Thanh et al. [7] applied the Gappy POD to a flowaround an airfoil and showed that the flow field can be successfully reconstructed from snapshots with 30% datamissing. However, the Gappy POD assumes that a missing portion differs in different snapshots so that it cannotbe applied to the situations where the missing region is fixed like Fig. 1. Use of the Kalman filter [8] and thelinear stochastic estimation [9] was also investigated for state estimation of turbulent channel flow. More recently,applications of machine learning have also emerged to take into account nonlinearities in the schemes [10–15]. Henryde Frahan et al. [16] utilized machine learning for recovering the computational fluid data of a cylinder wake andhomogeneous isotropic turbulence. Capability of machine learning based super resolution, which reconstructs high-resolution flow field from its low-resolution counterpart, was demonstrated by Fukami et al. [17, 18] for turbulent flows:the proposed model can substantially recover the energy spectra in the higher wavenumber range of two-dimensionaldecaying turbulence. The concept of super resolution has recently been extended to not only computational [19–21]but also experimental data sets [22].Of particular interest from the recent trends in reconstructing and estimating flow data using machine learning isits application to the PIV data processing procedure. To the best of our knowledge, the first application of machinelearning to PIV was demonstrated by Chen et al. [23]. They proposed a multi-layer perceptron based estimationmethod for PIV of a uniform flow and reported that the estimated mean velocity showed reasonable agreement withthe result of conventional PIV. Rabault et al. [24] utilized a convolutional neural network (CNN) for the first time toestimate velocity fields from synthetic particle images. Moreover, Cai et al. [25] have recently used a machine learningmodel called FlowNetS architecture [26] to estimate velocity fields of a cylinder wake, a backward facing step flow,and isotropic homogeneous turbulence from synthetic particle images and demonstrated that it can estimate higherresolution data than the conventional PIV. They have also applied the proposed method to experimental particleimages of turbulent boundary layer, and the results imply its great potential. Similarly, Grayver et al. [27] utilizeda supervised machine learning model to predict a velocity field from streak images of particles on both syntheticand experimental data. Gim et al. [28] have applied a shallow neural network to detect three-dimensional particlelocations for real-time particle detection in three-dimensional particle tracking velocimetry.Despite these recent trends, there are still only few studies focusing on applications of machine learning to ex-perimental data, and none of them handles cases with an object immersed in a flow which makes a fixed missingregion in PIV as exemplified in Fig. 1. This paper aims at estimating complete velocity fields from experimentalparticle images with a fixed missing region due to laser reflection etc. using a supervised machine learning method.An autoencoder-type CNN, which is robust for translation or rotation of images, is utilized for the present study.
FIG. 2. Overview of the machine learning based experimental flow data estimation. (a) Preparation of training data set. Theartificial particle image (API) at a certain time instant q r is generated by randomly distributing particles. The API with atime increment q (cid:48) can be obtained using the velocity data u DNS . (b) Schematic of the API used for the input to the machinelearning model, q = q r + q (cid:48) . (c) Training of machine learning model F for full data. We utilize artificial particle images (APIs) generated with a direct numerical simulation (DNS) to train the machinelearning model and apply the trained model to experimental images. The present paper is organized as follows. InSec. II, we introduce the overview, the training data, and the machine learning model of the present study. Theresults of estimation for APIs and experimental data are presented and discussed in Sec. III. At last, concludingremarks are offered in Sec. IV with some outlooks on experimental data estimation using data-oriented methods.
II. METHODSA. Overview of the present study
Let us introduce an overview of the present study in Fig. 2. For a training process of supervised machine learning,data sets comprising of pairs of an input q and the solution u must be prepared. Then, the model F is trainedsuch that u ≈ F ( q ). In this study, we aim to establish a machine learning model F to output the velocity fields u = { u, v } from input particle images q . Although our final objective is to apply this model to experimental data,pairs of particle images and the corresponding velocity fields obtained by the conventional particle image velocimetry(PIV) data processing method, e.g., the cross-correlation method, are not appropriate as the training data in termsof reliability due to the reflection of a laser light as mentioned above. To prepare clean training data, instead ofexperimental data with missing regions, we use synthetic particle images generated with the aid of direct numericalsimulation (DNS). Details for numerical setup of the present DNS can be found in Appendix.As the first step for preparing the training data set, we generate an artificial particle image (API) at a certain timeinstant q r as shown in Fig. 2(a). The particles are randomly distributed here. Details on the generation of APIswill be offered in Sec. II B. We utilize the DNS data to calculate the locations of these particles to obtain the APIwith a small time increment q (cid:48) . Then, the images at two time steps q r and q (cid:48) are summed up to obtain the API q as illustrated in Fig. 2(b). This procedure enables us to emulate the experimental particle images, which are notexact instantaneous snapshots because each particle would move a bit during the exposure time of high-speed camerawith real PIV experiment. Due to such movements of the particles, strictly speaking, the images contain very short FIG. 3. Examples of artificial particle images (APIs) around a square cylinder: (a) full API q = q r + q (cid:48) ; (b) lacked API (cid:98) q = (cid:98) q r + (cid:98) q (cid:48) . Lacked APIs (cid:98) q are prepared by applying a mask to the full APIs. trajectories of the particles [27]. The APIs at different time instants are generated likewise.As shown in Fig. 2(c), the next step is the construction of the machine learning model F for the full APIs q . Themodel F attempts to predict velocity fields u = { u, v } from the APIs q . Hence, the objective of training process forthe model F is to obtain optimized weights w so as to estimate the velocity fields from the APIs by minimizing thedifference between estimated velocity fields and reference DNS data used to generate APIs, such that w = argmin w || u DNS − F ( q ; w ) || . (1)In the present study, we choose L error norm as the loss function.Subsequently, we construct the other machine learning model (cid:98) F to estimate velocity fields u from the lacked APIs (cid:98) q , which is analogous to Fig. 2(c) but has a data missing region corresponding to the low correlation region in Fig.1(b). The lacked APIs (cid:98) q are generated by applying a mask to the full APIs q , and they include a lacked portionaround a square cylinder. Similar to the case of full APIs, the training process for the model (cid:98) F can be formulated as (cid:98) w = argmin (cid:98) w || u DNS − (cid:98) F ( (cid:98) q ; (cid:98) w ) || , (2)where (cid:98) w represents weights of the machine learning model (cid:98) F . At last, we evaluate the applicability of the machinelearned model (cid:98) F to experimental particle images. B. Artificial particle image
As introduced in Sec. II A, we use the artificial particle images (APIs) as the input for training the machine learningmodels. Figure 3 shows samples of full and lacked APIs. The APIs in the present study are generated based on theformulation proposed by Okamoto et al. [29]. The intensity I ∈ [0 ,
1] at location ( x, y ) is defined as I ( x, y ) = N (cid:88) i =1 I ,i exp (cid:18) − ( x − x p,i ) + ( y − y p,i ) ( d p / (cid:19) , (3) I ,i = 0 .
06 exp (cid:32) − z p,i σ l (cid:33) , (4)where ( x p , y p , z p ) is the location of a particle, d p is the particle diameter, σ l is the laser sheet thickness, and N is thenumber of particles, respectively. In this study, we set d p = 3 and σ l = 2 by referencing the experimental images ofour experimental setup. The intensity around each particle is defined with Gaussian distribution, whose maximumvalue should be tuned so as to roughly mimic the experimental images. In this study, we set it at 0 .
06 followingour preliminary tests. For the lacked APIs, we set a lacked region so as to mask the low correlation regions in ourexperimental situation and substitute zero there.Because the APIs must emulate experimental particle images, we further adjust the intensity distribution of theAPIs to that of the experimental images. For this adjustment, we utilize imhistmatch in MATLAB function, whichcan account for the first four statistical moments (i.e., average, variance, skewness factor, and flatness factor) of theintensity distribution. The values of those moments of the APIs with and without post-processing are summarized inTab. I. As the result of the histogram fitting, the intensity histogram of the post-processed API shows great agreementwith that of the experimental image. Although the variance is not improved, the post-processed API shows betterresults than without processing as a whole, since the average value is drastically improved. We also present in Fig. 4
TABLE I. Intensity fitting of APIs to the experimental image.Experimental image Artificial particle imagew/o processing w/ processingAverage 34.8 53.0 (52.6%) 34.5 (0.85%)Variance 8.24 8.60 (4.37%) 9.04 (9.71%)Skewness 2.34 2.53 (8.15%) 2.28 (2.66%)Flatness 6.44 7.74 (20.1%) 6.04 (6.23%)FIG. 4. An example of intensity fitting for APIs. (a) Experimental particle image, (b) API without and (c) with post-processing. the samples of experimental image and fitted APIs. It is found that the adjusted API can emulate the experimentalimage better than the original image. For readers’ information, we used a PIV laser G1000 (KATO KOKEN) anda high-speed camera K5-USB 8GB (KATO KOKEN) with AF Nikkor 85 mm F/1.6D (Nikon) lens to capture theparticle images of flow fields.
C. Machine learning model
In the present study, a model structure similar to a convolutional neural network (CNN) based autoencoder is usedto estimate velocity field u from the input particle image q . CNNs [30] have often been utilized in the field of imagerecognition, and recently, use of CNNs has also been propagated in the fluid dynamics community [31–39].The CNN is mainly composed of convolutional layers and pooling layers. The convolutional layer extracts keyfeatures of input data by filtering operation, c ( l ) ijm = ψ (cid:32) K − (cid:88) k =0 H − (cid:88) p =0 H − (cid:88) q =0 c ( l − i + p,j + q,k h pqkm + b m (cid:33) , (5)where c ( l − ijm and c ( l ) ijm are the input and output data at layer l , h pqkm denotes a filter of the size of ( H × H × K ) and b m is a bias. The output from the filtering operation is eventually multiplied by an activation function ψ . Using thenonlinear activation function here, a machine learning model can take into account nonlinearity in its structure. Inthe pooling layer, representative values are downsampled by pooling operations, e.g., maximum values (max pooling)or average values (average pooling). It is widely known that by incorporating the pooling operations CNN models canacquire robustness against the variance of input data due to decrease of spatial sensitivity [30]. Moreover, denoisingeffect can also be expected by several pooling operations, which is crucial for dealing with experimental images. FIG. 5. Autoencoder-based convolutional neural network used in the present study.
TABLE II. Details of the proposed autoencoder-based convolutional neural network.Layer Filter size β of Adam 0.9 β of Adam 0.999Number of epochs 3000 Early stopping patience 20 One of the recent hottest uses of CNNs is as an autoencoder [40], which is trained to output the same data asthe input via a low dimensional space called the latent space. Although the output of the present study is notthe same as the input as mentioned above, we here borrow the network structure of the CNN autoencoder. Letus present in Fig. 5 the proposed machine learning model for experimental data estimation. The details of theproposed model are summarized in Tab. II. As can be seen in both Fig. 5 and Tab. II, the input vector q is mappedinto a low dimensional latent space β with max pooling operations before the dimension is recovered at the outputlayer. As mentioned above, lower spatial sensitivity and robustness against noise can be acquired using the poolingoperations. As the pooling operation downsamples the representative values, e.g., max and average values, it lowersthe spatial sensitivity and generalizes the model. This feature is especially crucial for the present problem setting,since the location of each particle varies randomly in actual experimental images. Also, the denoising effect can beexpected through the pooling operations. Hence, the CNN autoencoder-like model which has the pooling procedureis suitable for our problem setting since we will apply a model (cid:98) F trained with lacked data to experimental data.We have checked that the present CNN autoencoder-like model outperforms a CNN model without the pooling andupsampling operations in our preliminary test. As the activation function, we use a rectified linear unit (ReLU),which has widely been utilized because of its stability in updating weights of a training process [41].For training of the present machine learning model, we use the Adam optimizer [42] and early stopping criteria[43] to avoid overfitting. The details of the hyperparameters are summarized in Tab. III. Note in passing that theaccuracy of the present model can be improved utilizing some theoretical optimization methods, e.g., hyperopt [44],Bayesian optimization [45–47], and random search [47], although we do not consider here. In the present study,baseline machine learning models of both single and staggered square cylinders are trained using 10000 x − y crosssectional snapshots, which correspond to 500 time steps per twenty positions in the z direction. You can see thedetails of preparation for training data set in Appendix. Note that dependence of the estimation accuracy on thenumber of the snapshots n snapshot will be investigated later. FIG. 6. Velocity fields estimated from artificial particle images (APIs) of flow around a square cylinder: (a) full APIs q ; (b)lacked APIs (cid:98) q .FIG. 7. (a) Mean centerline velocity profile and (b) probability density function (PDF) of streamwise velocity estimated withfull and lacked APIs of a flow around a square cylinder. III. RESULTS AND DISCUSSIONA. Example 1: Wake around a square cylinder
Let us first present in Fig. 6 the estimation using the full and lacked artificial particle images (APIs), q and (cid:98) q ,respectively. As shown here, the estimated fields with both attributes show great agreement with reference directnumerical simulation (DNS) data. The L error norm (cid:15) = || u DNS − u ML || / || u DNS || for both full and lacked APIsare approximately 10%, which indicates that the reconstruction ability is not influenced by the presence of missingregions in the input API. Note that the distribution of local L error norm of velocities are concentrated in the vortexshedding region rather than the missing region where the flow is almost steady in this particular demonstration.For the statistical evaluation, the mean streamwise velocity profile at y = 0 (i.e., the centerline velocity) and theprobability density function (PDF) of the streamwise velocity are shown in Fig. 7. The centerline velocities of the FIG. 8. Dependence on the number of snapshots and training time with a single cylinder.FIG. 9. Robustness against image intensities of the test data. fields estimated from the full and lacked APIs show almost perfect agreement with that of DNS. It also indicates thatpresence of the missing region has no influence on the estimation ability for this particular demonstration. The sametrends can be seen in the PDF, as shown in Fig. 7(b). Although there is a mismatch in the low probability region,for both the full and lacked APIs, the result of machine learned model shows great agreement with that of DNS.Since the machine learning models are trained to minimize the loss function, i.e., L error, it is more difficult for themachine learned model to estimate the low probability region.Let us examine the dependence of the L error norm on the number of snapshots used for training, as presentedin Fig. 8. Note again that the results presented above are obtained by the machine learned model trained with n snapshot = 10000. The errors of the estimated fields are relatively high in the case of lower n snapshot . This trend canbe seen with n snapshot = 50 in Fig. 8(a). The estimated fields of both attributes are blurry and not well matched withthe reference DNS. But, noteworthy here is that it can still estimate the large scale motion of the flow from the APIwith as little as n snapshot = 50. The error of the estimated field suddenly decreases around at n snapshot = 1000 andconverges approximately towards (cid:15) = 1 × − with n snapshot = 5000 to 10000. It suggests that although the trainingtime here is approximately 20 minutes on the NVIDIA TESLA V100 graphics processing unit (GPU) even in the casewith largest number of n snapshot , users should care the trade-off relationship between the computational burden andthe accuracy of the model.Next, we check the applicability and robustness to experimental situations. For the discussions above, the input APIsare generated with a single intensity I G, train . Since the intensity of experimental images varies on each experimentalsetup, machine learning models are required to be robust to the intensity variance unless the amount and kinds oftraining data are increased. Here, let us assess the robustness against the intensity of test images, as presented in Fig.9. The intensity of each image I G is defined as an average of intensity for each pixel of images. The estimation usingthe test images with the same intensity as the training data, i.e., I G, test /I G, train = 1, results in the lowest mean L error norm, as expected. As the test images become darker or brighter than the training images, the L error normincreases, but the maximum difference in the L error norm is less than 2%, as can be seen in Fig. 9, which suggeststhat the present machine learned model is robust against intensity variations in the test data.We then evaluate robustness against a noisy input q noise = q + κ n , where n denotes the Gaussian noise and κ is itsamplitude, as presented in Fig. 10. This assessment can be one of the benchmarks for applications to experiments.Note here that the scale of the vertical axis in Fig. 10 is normalized by the error in the case without noise. Up to κ = 0 .
10, the flow fields show no significant difference against the reference field presented on the upper-left portionof the error plot. However, with κ = 0 .
20, the finer scale motion in the wake region cannot be estimated, although
FIG. 10. Robustness against noise of the proposed machine learned model for a square cylinder flow.FIG. 11. Output of each hidden layer with an example of a single square cylinder flow. the large scale motion around the bluff body is still captured.Before applying our machine learned model to experimental data, we investigate in which regions the machinelearned model has strong interests for estimating the velocity fields. This is done by visualizing the output of eachhidden layer. Example outputs of each layer are shown in Fig. 11. Since each hidden layer has 16 or 32 outputsaccording to the number of filters, the output data with the largest average intensity is shown as the representatives.Note that the outputs of each hidden layer have only positive values because the rectified linear unit (ReLU) functionis utilized as the activation function in the present study, and therefore, larger values indicate stronger interests. Theoutputs of 2nd and 4th layer clearly show that the machine learned model has a strong interest on the location ofthe square cylinder. In contrast, for the layers near the output, e.g. 12th and 13th layers, the machine learned model
FIG. 12. Visualization at 14th linear layer. Orders above each figure indicate their contributions to the output. FIG. 13. Application of the machine learned model to experimental particle images. (a) Mean centerline velocity; (b) Reynoldsshear stress u (cid:48) v (cid:48) .FIG. 14. Comparison of velocity fields between the present model and the cross-correlation method. focuses on the wake region. This observation implies that the machine learned model first finds out the location of thesquare cylinder and then focuses on the regions of velocity fluctuations so as to output appropriate velocity fields. Weare also able to see clearly the second trend from the output of 13th layer, as summarized in Fig. 12. We show here thefour most dominant outputs. A noteworthy observation is that the positive values and negative values of velocity fieldare separated in different channels. This is perhaps because the machine learned model needs to represent negativevalue of velocity field at 14th layer, since we use ReLU — which has only positive range — as the activation functionin the upstream layers. Although we visualized the filters in convolutional neural network (CNN) such as Figs. 11and 12, different tools, e.g., Grad-CAM [48–50], can also be considered to interpret the results obtained by machinelearning models.We then apply the machine learned model to estimate the velocity field from our experimental images. To evaluatethe result, the centerline velocity and distribution of the Reynolds shear stress u (cid:48) v (cid:48) are assessed in Fig. 13, since thecorrect answers corresponding to the instantaneous experimental images do not exist. As shown in Fig. 13(a), themean velocity profile of the estimated fields is in excellent agreement with the reference DNS data. For the Reynoldsshear stress in Fig. 13(b), the structure can be captured well, which confirms that the fluctuation components canalso be estimated successfully, although the absolute value of the estimated field is slightly lower than that of DNS.An instantaneous flow field estimated using the machine learned model is also compared to the conventional cross-correlation method [1], as shown in Fig. 14. The flow field can be recovered successfully by using the machine learnedmodel. On the other hand, the result with the cross-correlation method is grossly affected by the missing regions. Itsuggests that the present machine learned model can retain the reconstruction ability by explicitly giving the lackedportion to training data. What is also noteworthy here is the density of velocity vectors obtained with the machinelearning based estimation compared to the conventional cross-correlation method. The machine learning based modelcan provide much finer resolution than the cross-correlation method, which can provide only one velocity vector froman interrogation window including some particles. Because the size of interrogation window is determined based onthe experimental setup, such as the particle diameter, the particle number density, and the spatial resolution of thecamera, it is usually tougher to obtain finer structures of flow field than a numerical simulation. In contrast, themachine learning model is trained to estimate vectors on every single computational grid points, it is able to estimatefiner structures of the field [25]. Particularly in our cases, the machine learned model can provide approximately 40folds denser field than the cross-correlation method. This advantage enables us to find small-scale structures which1 FIG. 15. Velocity fields estimated from APIs of flow around two staggered square cylinders. : (a) full APIs q ; (b) lacked APIs (cid:98) q .FIG. 16. Six dominant POD modes of a vorticity field around two staggered square cylinders. The value shown under eachfigure represents the contained energy of each mode. cannot be captured with the conventional methods. B. Example 2: Wakes around two staggered square cylinders
The problem setting of example 1 can be regarded as a relatively easy task for the machine learning model toaugment the velocity field of a lacked region because the flow there was nearly steady. For further investigation onthe capability of the machine learning model, let us consider a flow around two staggered square cylinders [51].The flow fields estimated from APIs are shown in Fig. 15. With both full and lacked inputs, the estimated flow2
FIG. 17. (a) Mean centerline velocity and (b) probability density function (PDF) of streamwise velocity estimated from fulland lacked APIs of flow around two staggered square cylinders.FIG. 18. Application of the machine learned model to experimental images of flow around two staggered square cylinders. (a)Mean centerline velocity; (b) Reynolds shear stress u (cid:48) v (cid:48) . fields are in reasonable agreement with the reference DNS data. The L error norms listed underneath the figuresshow approximately 10%, which are almost the same as that of example 1. However, we should note that the error isconcentrated in the region below the second square cylinder, where the complex structure is hidden as shown in Fig.15. It can also be seen from the proper orthogonal decomposition (POD) analysis with vorticity field as assessed inFig. 16. It is striking that the complex structures below the second square cylinder can be seen on POD modes 4 and5, a sum of which contains approximately 15% of the kinetic energy. Since there is no significant difference in the L error norm between full and lacked APIs, this POD analysis tells us the fact that the machine learned model for thelacked APIs (cid:98) F was able to estimate these complex structures as well as the model for full APIs F even the structuresare masked in advance.For further assessment, the mean centerline velocity profile and the PDF are presented in Fig. 17. The overallprofile of the mean centerline velocity agrees well with the DNS data, although the slight mismatch can be seen onthe vortex shedding region. The reason for this mismatch can be also seen in the PDF, which suggests that themachine learned model fails to estimate the low probability events. This is, again, likely because the training processis designed to minimize the L error, to which lower probability events have less contributions.For evaluation of the result with the experimental image inputs of a flow around two staggered square cylinders,the mean centerline velocity profile and the Reynolds shear stress u (cid:48) v (cid:48) distribution are assessed in Fig. 18. Similar tothe result on the APIs, the mean centerline velocity profile of the estimated field is in reasonable agreement with thatof the reference DNS, although there is a slight mismatch near the vortex shedding region. This trend can be seenalso in the Reynolds shear stress u (cid:48) v (cid:48) in Fig. 18(b), which shows underestimation. This is likely because the structureis substantially more complex than the first example and the low accuracy region corresponds to the high-energeticportion of POD modes 4 and 5, as discussed above.For this example, too, finer fields can be obtained by the present machine learning model, as presented in Fig.19. Similar to the single cylinder example, denser field can be seen with the use of the machine learned model. It isstriking that the differences in the velocity magnitude are clearer than that with the conventional cross-correlationdue to the finer structures. Although we have no correct answer for this experiment, the results suggest again thatthe machine learned model enables us to extract the hidden structures of flow fields.Finally, we perform three kinds of tests to investigate the influence of the alignment of bluff bodies on the machinelearning based estimation as follows:3 FIG. 19. Comparison of velocity fields around two staggered square cylinders between the present model and the cross-correlation method.FIG. 20. Estimation of flows around (a) a square cylinder and (b) two staggered square cylinders by a pre-trained model withthe other field and a model trained with both fields.
1. Apply the machine learned model pre-trained with flows around two staggered square cylinders to a flow arounda single square cylinder.2. Apply the machine learned model pre-trained with flows around a single square cylinder to a flow around twostaggered square cylinders.3. Apply the machine learned model pre-trained with both flows to both test data.Through these tests, we can see whether the present model can be utilized to a flow around unseen bluff bodyalignments or not. Figure 20 summarizes the results of these tests. For Fig. 20(a), APIs of flow around a squarecylinder are used as the test data. The model pre-trained with two square cylinders fails to estimate the correct field— the estimation is highly affected by the alignment of square cylinders of training data. The same trends can be seenin the case of estimating a flow around two square cylinders by using a model trained with a single square cylinder,as shown in Fig. 20(b). Since the machine learned model trained with a single flow configuration has relatively lowapplicability to flows around unseen bodies, we further assess the ability of the model trained using both types of flowfields as stated as procedure 3 above. As presented in Fig. 20, both estimated flow fields are in nice agreement withthe reference DNS data. The L error norm shown underneath each estimated field, is approximately 10%, whichis almost the same as those in the previous examples presented in Figs. 6 and 15. These results indicate that weneed to be careful in preparing the training data for machine learning based particle image velocimetry (PIV) dataaugmentation, especially in the case of estimation on flows around bluff bodies. IV. CONCLUSIONS
We presented a machine learning based data estimation method for particle image velocimetry (PIV) consideringa typical experimental limitations, i.e., missing regions. The training data set was prepared using a direct numerical4simulation (DNS) and artificial particle images (APIs). Two types of data, i.e., full APIs (no data missing region) andlacked APIs (with a fixed data missing region), were considered to examine the applicability of the machine learningmodel to experimental situations. The lacked regions were set by referencing to our exact wind tunnel setup. Weutilized a model similar to a convolutional neural network based autoencoder as the machine learning model.As the first example, we considered a flow around a single square cylinder at Re D = 300. The estimated flow fieldswith both full and lacked APIs showed great agreement with the reference DNS data. It was found that the machinelearned model can estimate the flow field from experimental images which has a lacked region, while the conventionalcross correlation method is heavily affected by the lacked region. These trends could be seen also in the statisticalassessment using the mean centerline velocity profile and the Reynolds shear stress distribution. In addition, by usingthe machine learning model, finer velocity fields could be obtained than that with a cross correlation method.For further assessment on more complex flow fields, we applied the proposed method to a flow around two staggeredsquare cylinders. Despite that the flow field contains finer and complex structures, the machine learned model couldestimate the velocity fields from APIs with L error norm less than 10%. We also found that even the complexstructures, with the energy rate of 15%, are hidden by the alignment of two square cylinders, the machine learnedmodel for lacked data is able to estimate the field as well as the model for full data.We also investigated the applicability of pre-trained model to a flow around unlearned alignment of bluff bodies.The results indicated that estimation with the machine learned model is affected by the alignment of bluff bodies. Wealso found that the machine learned model could be able to establish the function in not only estimating the field,but also recognizing the alignment of bluff bodies if proper training data sets are given. For more practical uses, caresshould be taken for the choice of training data set, which should be further examined in future.Recently, Hasegawa et al. [52, 53] have demonstrated that it is able to estimate the time series of flow field arounda cylinder by using a machine learned reduced order model at unlearned Reynolds numbers. This result encourages usto extended the present method to accommodate various Reynolds numbers. Also by properly choosing the trainingdata set, as mentioned above, we believe that the present concept, which has the ability to estimate denser fieldsthan a conventional cross correlation method, can be applied to more complex flows which contain various scales ofstructures, e.g., flows around an airfoil and turbulent flows. ACKNOWLEDGEMENTS
This work was supported through JSPS KAKENHI Grant Number 18H03758 by Japan Society for the Promotionof Science. The authors acknowledge Mr. Takaaki Murata, Mr. Taichi Nakamura (Keio University) and Mr. KazutoHasegawa (Keio University, Polimi) for fruitful discussions. The authors also thank Mr. Hikaru Murakami (KeioUniversity) for sharing his DNS data.
APPENDIX: DETAILS OF DNS AND DATA PREPARATION
We consider a flow around a square cylinder at the Reynolds number Re D = 300 based on the free-stream velocityand the side length of the cylinder D . The training data set is prepared by a direct numerical simulation (DNS) withnumerically solving of the incompressible Navier–Stokes equations with a penalization term [54], ∇ · v = 0 (6) ∂ v ∂t + ∇ · ( vv ) = − ∇ p + 1Re D ∇ v + λχ ( v b − v ) (7)where v = { u, v, w } and p and Re D are the non-dimensionalized velocity vector and pressure, respectively. Thepenalization term represents an object with a penalty parameter λ , a mask value χ , and a velocity vector of a flowinside the object v b which is 0. The mask value χ takes 0 outside of an object and 1 inside of an object.As shown in Fig. 21(a), the size of the computational domain here is ( L x × L y × L z ) = (20 D × D × D ). Thecomputational time step is ∆ t = 5 . × − , which results in the maximum Courant number around 0 .
3. Our DNScode has been verified with Franke et al. [55] and Robichaux et al. [56].For the training data of the first case with a single square cylinder, we focus on the volume around the square cylindersuch that (7 D × D × . D ) shown by red lines in Fig. 21(a). The number of grid points of the extracted region is( N (cid:93)x × N (cid:93)y × N (cid:93)z ) = (140 × × D = 300 [57], twenty x − y cross sections in z direction are used for the training data to account for the three-dimensionality. For the trainingdata we use 10000 x − y cross sectional snapshots, which include 500 time steps per twenty positions in z direction,as illustrated in Fig. 22. For the test data, we extract 160 x − y cross sections from five different instantaneous5 FIG. 21. Computational domain and vorticity field of flow around (a) a square cylinder and (b) two staggered square cylinders.FIG. 22. Preparation for training data in the present study. Twenty x − y sections of three-dimensional data with 500 timesteps are extracted. three-dimensional snapshots; namely, 800 snapshots are utilized for evaluation. Note that each instantaneous field ischosen to be distanced from each other, which enables us to evaluate the estimation ability by various states of theflow.The training data for the two staggered square cylinders are prepared similarly to the first example. The computa-tional domain with the size and grid points are summarized in Fig. 21(b). The Reynolds number is set to Re D = 300and the time step for the DNS is ∆ t = 0 . [1] R. J. Adrian, Twenty years of particle image velocimetry, Exp. Fluids (2005).[2] B. Horn and B. Schunck, Determining optical flow, Artif. Intell. , 185 (1981).[3] E. Paterna, P. Moonen, V. Dorer, and J. Carmeliet, Mitigation of surface reflection in PIV measurements, Meas. Sci.Technol. , 057003 (2013).[4] N. R. Zare, My life with LIF: A personal account of developing laser-induced fluorescence, Annu. Rev. Anal. Chem. , 1(2012).[5] R. Jain, M. Mathur, S. Sikarwar, and A. Mittal, Removal of the hazardous dye rhodamine B through photocatalytic andadsorption treatments, J. Environ. Manage. , 956 (2007).[6] R. Everson and L. Sirovich, Karhunen?lo`eve procedure for gappy data, J. Opt. Soc. Am. A , 8 (1995). [7] T. Bui-Thanh, M. Damodaran, and K. Willcox, Aerodynamic data reconstruction and inverse design using proper orthog-onal decomposition, AIAA J. , 8 (2004).[8] M. Chevalier, J. Hœpffner, T. R. Bewley, and D. S. Henningson, State estimation in wall-bounded flow systems. part 2.turbulent flows, J. Fluid Mech. , 167 (2006).[9] T. Suzuki and Y. Hasegawa, Estimation of turbulent channel flow at R e τ = 100 based on the wall measurement using asimple sequential approach, J. Fluid Mech. , 760 (2017).[10] M. P. Brenner, J. D. Eldredge, and J. B. Freund, Perspective on machine learning for advancing fluid mechanics, Phys.Rev. Fluids , 100501 (2019).[11] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, Machine learning for fluid mechanincs, Annu. Rev. Fluid Mech. ,477 (2020).[12] K. Fukami, K. Fukagata, and K. Taira, Assessment of supervised machine learning methods for fluid flows, Theor. Comput.Fluid Dyn. available online, https://doi.org/10.1007/s00162-020-00518-y (2020).[13] J. Rabault, M. Kuchta, A. Jensen, U. R´eglade, and N. Cerardi, Artificial neural networks trained through deep reinforce-ment learning discover control strategies for active flow control, J. Fluid Mech. , 281 (2019).[14] L. Guastoni, M. P. Encinar, P. Schlatter, H. Azizpour, and R. Vinuesa, Prediction of wall-bounded turbulence from wallquantities using convolutional neural networks, arXiv preprint, arXiv:1912.12969 (2019).[15] P. A. Srinivasan, L. Guastoni, H. Azizpour, P. Schlatter, and R. Vinuesa, Predictions of turbulent shear flows using deepneural networks, Phys. Rev. Fluids , 054603 (2019).[16] M. Henry de Frahan and R. Grout, Data recovery in computational fluid dynamics through deep image priors, arXivpreprint, arXiv:1901.11113 (2019).[17] K. Fukami, K. Fukagata, and K. Taira, Super-resolution reconstruction of turbulent flows with machine learning, J. FluidMech. , 106 (2019).[18] K. Fukami, K. Fukagata, and K. Taira, Super-resolution analysis with machine learning for low-resolution flow data, in , Paper 208 (2019).[19] R. Onishi, D. Sugiyama, and K. Matsuda, Super-resolution simulation for real-time prediction of urban micrometeorology,Scientific Online Letters on the Atmosphere (SOLA) , 178 (2019).[20] B. Liu, J. Tang, H. Huang, and X.-Y. Lu, Deep learning methods for super-resolution reconstruction of turbulent flows,Phys. Fluids , 025105 (2020).[21] K. Fukami, K. Fukagata, and K. Taira, Machine learning based spatio-temporal super resolution reconstruction of turbulentflows, arXiv preprint, arXiv:2004.11566 (2020).[22] Z. Deng, C. He, Y. Liu, and K. C. Kim, Super-resolution reconstruction of turbulent velocity fields using a generativeadversarial network-based artificial intelligence framework, Phys. Fluids , 125111 (2019).[23] P. H. Chen, J. Y. Yen, and J. L. Chen, An artificial neural network for double exposure PIV image analysis, Exp. Fluids , 373 (1998).[24] J. Rabault, J. Kolaas, and A. Jensen, Performing particle image velocimetry using artificial neural networks: a proof-of-concept, Meas. Sci. Technol. , 125301 (2017).[25] S. Cai, S. Zhou, C. Xu, and Q. Gao, Dense motion estimation of particle images via a convolutional neural network, Exp.Fluids , 60 (2019).[26] A. Dosovitskiy, P. Fischer, E. Ilg, P. Hausser, C. Haz?rbas, V. Golkov, P. Smagt, D. Cremers, and T. Brox, Flownet:Learning optical flow with convolutional networks, Proc. IEEE Int. Conf. on Computer Vision (2015).[27] A. V. Grayver and J. Noir, Particle streak velocimetry using ensemble convolutional neural networks, Exp. Fluids (2020).[28] Y. Gim, D. K. Jang, D. K. Sohn, H. Kim, and H. S. Ko, Three-dimensional particle tracking velocimetry using shallowneural network for real-time analysis, Exp. Fluids , 26 (2020).[29] K. Okamoto, S. Nishio, T. Saga, and T. Kobayashi, Standard images for particle-image velocimetry, Meas. Sci. Technol. , 685 (2000).[30] Y. A. LeCun, L. Bottou, Y. Bengio, and P. Haffner, Automatic early stopping using cross validation: quantifying thecriteria, Proc. IEEE , 2278 (1998).[31] K. Fukami, Y. Nabae, K. Kawai, and K. Fukagata, Synthetic turbulent inflow generator using machine learning, Phys.Rev. Fluids , 064603 (2019).[32] S. Lee and D. You, Data-driven prediction of unsteady flow fields over a circular cylinder using deep learning, J. FluidMech. , 217 (2019).[33] H. Salehipour and W. R. Peltier, Deep learning of mixing by two ‘atoms’ of stratified turbulence, J. Fluid Mech. , R4(2019).[34] T. Murata, K. Fukami, and K. Fukagata, Nonlinear mode decomposition with convolutional neural networks for fluiddynamics, J. Fluid Mech. , A13 (2020).[35] K. Hasegawa, K. Fukami, T. Murata, and K. Fukagata, Machine-learning-based reduced order modeling for unsteady flowsaround bluff bodies of various shapes, Theor. Comp. Fluid Dyn. (accepted) (2020).[36] J. Huang, H. Liu, and W. Cai, Online in situ prediction of 3-D flame evolution from its history 2-D projections via deeplearning, J. Fluid Mech. (2019).[37] J. Kim and C. Lee, Prediction of turbulent heat transfer using convolutional neural networks, J. Fluid Mech. (2020).[38] S. Lee and D. You, Mechanisms of a convolutional neural network for learning three-dimensional unsteady wake flow, arXivpreprint, arXiv:1909.06042 (2019).[39] R. Maulik, B. Lusch, and P. Balaprakash, Reduced-order modeling of advection-dominated systems with recurrent neuralnetworks and convolutional autoencoders, arXiv preprint, arXiv:2002.00470 (2020). [40] G. E. Hinton and R. R. Salakhutdinov, Reducing the dimensionality of data with neural networks, Science , 504 (2006).[41] V. Nair and G. E. Hinton, Rectified linear units improve restricted boltzmann machines, Proc. Int. Conf. Mach. Learn. ,807 (2010).[42] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint, arXiv:1412.6980 (2014).[43] L. Prechelt, Automatic early stopping using cross validation: quantifying the criteria, Neural Netw. , 761 (1998).[44] J. Bergstra, D. Yamins, and D. Cox, Making a science of model search: Hyperparameter optimization in hundreds ofdimensions for vision architectures, International conference on machine learning , 115 (2013).[45] E. Brochu, V. Cora, and de N. Freitas, A tutorial on bayesian optimization of expensive cost functions, with applicationto active user modeling and hierarchical reinforcement learning, Technical Report TR-2009-023, University of BritishColumbia (2009).[46] R. Maulik, A. Mohan, B. Lusch, S. Madireddy, and P. Balaprakash, Time-series learning of latent-space dynamics forreduced-order model closure, Physica D: Nonlinear Phenomena , 132368 (2020).[47] H. F. S. Lui and W. R. Wolf, Construction of reduced-order models for fluid flows using deep feedforward neural networks,J. Fluid Mech. , 963 (2019).[48] R. R. Selvaraju, A. Das, R. Vedantam, M. Cogswell, D. Parikh, and D. Batra, Grad-CAM: Why did you say that?, arXivpreprint, arXiv:1611.07450 (2016).[49] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, Grad-CAM: Visual explanations from deepnetworks via gradient-based localization, in Proceedings of the IEEE international conference on computer vision (2017)pp. 618–626.[50] E. Jagodinski, X. Zhu, and S. Verma, Uncovering dynamically critical regions in near-wall turbulence using 3D convolutionalneural networks, arXiv preprint, arXiv:2004.06187 (2020).[51] M. M. Alam, H. Bai, and Y. Zhou, The wake of two staggered square cylinders, J. Fluid Mech. , 475 (2016).[52] K. Hasegawa, K. Fukami, T. Murata, and K. Fukagata, Data-driven reduced order modeling of flows around two-dimensionalbluff bodies of various shapes, ASME-JSME-KSME Joint Fluids Engineering Conference, San Francisco, USA (2019).[53] K. Hasegawa, K. Fukami, T. Murata, and K. Fukagata, CNN-LSTM based reduced order modeling of two-dimensionalunsteady flows around a circular cylinder at different Reynolds numbers, (in Review) (2020).[54] J. P. Caltagirone, Sur linteraction fluide-milieu poreux: application au calcul des efforts excerses sur un obstacle par unfluide visqueux, C. R. Acad. Sci. Paris , 571 (1994).[55] R. Franke, W. Rodi, and B. Schonung, Numerical calculation of laminar vortex shedding flow past cylinders, J. Wind Eng.Ind. Aerodyn. , 237 (1990).[56] J. Robichaux, S. Balachandar, and S. P. Vanka, Three-dimensional floquet instability of the wake of square cylinder, Phys.Fluids , 560 (1999).[57] H. Bai and M. M. Alam, Dependence of square cylinder wake on Reynolds number, Phys. Fluids30