On the prediction of the turbulent flow behind cylinder arrays via echo state networks

This study aims at the prediction of the turbulent flow behind cylinder arrays by the application of Echo State Networks (ESN). Three different arrangements of arrays of seven cylinders are chosen for the current study. These represent different flow regimes: single bluff body flow, transient flow, and co-shedding flow. This allows the investigation of turbulent flows that fundamentally originate from wake flows yet exhibit highly diverse dynamics. The data is reduced by Proper Orthogonal Decomposition (POD) which is optimal in terms of kinetic energy. The Time Coefficients of the POD Modes (TCPM) are predicted by the ESN. The network architecture is optimized with respect to its three main hyperparameters, Input Scaling (INS), Spectral Radius (SR), and Leaking Rate (LR), in order to produce the best predictions in terms of Weighted Prediction Score (WPS), a metric leveling statistic and deterministic prediction. In general, the ESN is capable of imitating the complex dynamics of turbulent flows even for longer periods of several vortex shedding cycles. Furthermore, the mutual interdependencies of the TCPM are well preserved. However, optimal hyperparameters depend strongly on the flow characteristics. Generally, as flow dynamics become faster and more intermittent, larger LR and INS values result in better predictions, whereas less clear trends for SR are observable.


Introduction
Machine learning (ML) provides the possibility of modeling turbulent flows without solving the complex underlying equations [1][2][3][4][5][6][7].The possible outcomes can range from the reconstruction of the collective response of the turbulent flow in a confined geometry (such as cumulative heat transfer or aerodynamic forces) to the complete reconstruction of the flow field itself.These reconstructions can be a modeling of the turbulent flow or a prediction.Predicting the exact values of the velocity is far more challenging (or even impossible) than modeling an output with similar dynamics.Therefore, a temporal prediction of the velocity field in a turbulent flow is the ultimate objective, one that might never be entirely achieved but is undeniably worthwhile to strive for.
Due to the inherent chaotic nature of turbulence, temporal predictions diverge rapidly in the course of time.One possible solution for that is to feed the ML algorithms with some limited information about the flow continuously.This can be achieved through limited sensor data in the field [8][9][10] or complete input of some secondary variables, such as temperature or density, to the algorithm [11,12].Bright et al [8].were able to reconstruct the flow from sparse pressure sensor measurements over the surface of a cylinder for the case of a von Kármán vortex street (KVS).Later, Raissi et al [11].conducted pioneering work by reconstructing the pressure and velocity fields behind a bluff body from available flow visualization data.
Despite the challenges of the independent temporal prediction of turbulence without any input information during the prediction phase, studies have already begun to address this issue to some extent.Sekar et al [13].conducted one of the pioneering works to predict the velocity field over an airfoil via a combination of convolutional neural networks and deep Multilayer Perceptrons (MLP).Meanwhile, Srinivasan et al [14].combined MLPs with long short term memory networks to predict turbulent shear flows.Later, physical principles were incorporated with ML in the form of physics-informed deep learning to assist the algorithms in the prediction of turbulence [15][16][17].
Complex deep learning algorithms combined with physical principles can certainly improve predictions if they are carefully chosen and adjusted.However, the complexity and computational needs of such algorithms increase dramatically.Furthermore, it has not been proven yet that more complex algorithms would necessarily result in better modeling and predictions of turbulence in general.In this regard, reservoir computing [18,19] or, more specifically, Echo State Networks (ESN) [20] are promising algorithms that feature simple structure, fast computation, and proven learning capabilities for time series.An ESN consists of a reservoir of sparsely connected neurons, in which only the output layer is trained, and the neurons retain a memory of their previous state [21][22][23].Several studies have implemented ESNs to model 2-dimensional Rayleigh-Bénard convection [24][25][26][27] and showed the capability of ESNs in imitating the complex dynamics of turbulent flows.
Recently, Ghazijahani et al [28] predicted the unsteady KVS behind a cylinder.There, data reduction was applied via Proper Orthogonal Decomposition (POD) to the flow, and then the ESN predicted the Time Coefficients of the POD Modes (TCPM).The results demonstrated the potential of the ESN for the prediction of unsteady dynamics of turbulent KVS for a relatively long period of time.However, there is still a lot to be learned about the application of ESNs to more complex turbulent flows.In this regard, the flow past an array of cylinders is an ideal test case [29].Varying the distance between the cylinders, the dynamics of the turbulent wake significantly changes and shifts from single bluff body flow to transient flow and finally to co-shedding flow.Furthermore, once POD is applied to the flow, the increasing complexity of the flow shows also in the higher relative values of kinetic energy for higher modes.
Hence, the current study aims to investigate the predictive capability of an ESN in such complex dynamics.Furthermore, it is interesting to investigate the possible impact of the array variation and, therefore, the dynamics on the performance of ESNs in general and using optimized hyperparameters in particular.Thus, three different arrangements of cylinder arrays (see Ghazijahani and Cierpka [29] for more details) with single bluff body, transient and co-shedding flow regimes are chosen for the current study.They are named in the following mVnH, where the first number stands for the vertical distance and the second for the horizontal distance in diameter of the cylinders, respectively.Then, the flow is reduced by POD, and the TCPM are predicted by an ESN with a reservoir of 3000 neurons.The three main hyperparameters of the ESN, namely input scaling (INS), spectral radius (SR), and leaking rate (LR), are varied to investigate their respective influence.Finally, the results are visualized and quantified to assess the predictions deterministically and statistically.

Experimental setup
Figure 1 shows a photo of the experimental setup for the Particle Image Velocimetry (PIV) measurements behind the cylinders.Seven cylinders were arranged in two columns, four in the rear and three in the front.The cylinders were rigid and had a diameter of D = 1 mm.The vertical distance of the cylinders in each column (V/D) and the horizontal distance between the two columns (H/D) were varied for V/D = [2,4,6] and H/D = [2,4,8], respectively.This resulted in nine different arrays of cylinders, each with its distinct flow features.The arrays were placed in a water channel with a cross-section of 50 × 50 mm 2 .The free stream velocity was set to V ∞ = 133 mm s −1 and thus the Reynolds number based on the cylinder diameter was Re = V ∞ D/ν ∼ 100.The flow was seeded with polyamide particles of 5 µm in diameter and hallow glass spheres of 10 µm in diameter.The latter was used only for arrays of V/D = 6 due to the necessity of larger particles because of the larger field of view in these cases.
A continuous-wave laser (Laserworld Green-200532) illuminated the vertical mid-plane of the channel with a light sheet with a thickness of 1 mm.Then, PIV images of the illuminated particles were captured using a high-speed camera (HS 4 M by LaVision GmbH) perpendicular to the laser sheet from the outside of the water channel.The images were calibrated against a calibration target.The PIV double-frame images were recorded at 200 Hz for 15 s.This corresponds to 330 vortex-shedding events for the wake of a single isolated cylinder, which is sufficient for data convergence.
An advanced cross-correlation evaluation was carried out for PIV processing via DaVis 10 (LaVision GmbH) for an initial rectangular interrogation window size of 64 × 64 pixels with 50% overlap, and a final circular Gaussian window weighting of 16 × 16 pixels again with 50% overlap.This yielded a final spatial  resolution of 0.069 × 0.069 mm 2 for the arrays with V/D = 2 and 4, and 0.12 × 0.12 mm 2 for the arrays with V/D = 6.For all arrays, the field of view consisted of 231 × 138 vectors.The number of outlier vectors remained under two percent for all measurements, applying a normalized median test [30].The dynamic spatial and velocity ranges were [0-14] pixels and [0-150] mm s −1 , respectively.Dynamic spatial range is the maximum length of the field of view divided by the vector spacing, whereas dynamic velocity range is the maximal velocity divided by the uncertainty (often taken as 0.1 or 0.05 pixel displacement).For properly adjusted experiments, it has usually been assumed that the absolute error is in the range of 0.1 pixels [31].This results in minimum relative uncertainties of 1% and 1.3% for V = 2, 4 and V = 6 in the free stream.

ESN
In this study, an ESN is used to predict the TCPM of the velocity fields behind three different cylinder arrays with different arrangements, namely 2V2H, 4V2H, and 6V2H.An ESN is a type of recurrent neural network where only the output layer is trained.Figure 2 shows a schematic sketch of the ESN.The ESN consists of a reservoir of N neurons that are randomly connected to each other with a weight matrix W (green arrows) and to the input signals with a weight matrix W in (blue arrows).The weight matrices of these connections (W and W in ) are randomly generated at the beginning of the training process, and their random state is fixed with a variable named as Random Seed (RS).The network thus takes the past values from previous time steps and uses a nonlinear and randomly generated reservoir state to learn the output connection to best predict the future time steps.Then each neuron state s(n) is updated by following equations: Here, LR stands for the LR, and it is the speed at which the neurons are updating their state.Next, the reservoir starts to learn some connections (W out , red arrows in figure 2) from the neurons to the output signals q(n) so that it can predict them accurately: where Here, β > 0 stands for the ridge regression parameter, which is used to prevent the amplification of small differences in state dimensions by large rows of W out .Furthermore, it prevents the algorithm from overfitting, which occurs when it learns the training data by heart and performs poorly on new unseen data.
For the current study, the ESN is written in Python using the easyesn library [32].It consists of N = 3000 neurons randomly connected to 20 percent of the other neurons in the reservoir.The INS limits the input weights to The network is trained for 700 time steps and tested for another 700 time steps to assess its capability in prediction of the TCPM.For more details of the flow physics and vortex shedding for the different cases the interested reader is referred to Ghazijahani and Cierpka [29].There are three main hyperparameters that are believed to have the most important role in the performance of an ESN, their values are varied to investigate how they affect the predictions.The LR, is responsible for the speed at which the reservoir is updating; thus, its optimum value can change with the dynamics of the data.The second main hyperparameter is the SR which is the maximum eigenvalue of W. In other words, it is the chaotic degree of the nonlinear interactions between the neurons.Thus, the higher the SR value, the more chaotic the interactions between the neurons.However, in general, SR is recommended to be kept below one in order to ensure the Echo State Property [21], which is the independence of the network from its initial state and dependence only on the past input.INS is the third main hyperparameter, which limits the input weights and, thus, the relative weight of the input against the reservoir's internal history (in the prediction, input will be the output of the reservoir itself).Finally, the random realization of the reservoir itself can impact the predictions.However, it has been shown previously that around 24 RS are enough to make a conclusion about the general performance of the reservoir in a specific hyperparameter set [28].Thus, for the current study, LR and SR are varied for 10 different values between [0,1], and INS is varied for ten values between [0, 20].Moreover, for each hyperparameter set, 24 RS are tested to draw conclusions independent of the effect of the random realizations of the ESN.

Results
The KVS behind a single cylinder has been used extensively in the literature for ML applications.The turbulent flow behind arrays of cylinders showed very rich and complex dynamics that were highly dependent on the distances between [29].Three different general flow regimes were observed: a single bluff body, transient flow, and co-shedding flow.
In order to provide examples of rich, complex, yet different turbulent flows for the current machine learning algorithm (ESN), one example out of each of the three aforementioned flow regimes is chosen.2V2H for the bluff body flow, 4V2H for the transient flow, and 6V2H for the co-shedding flow.Then, snapshot POD [33] was applied to the velocity fields to reduce the amount of data fed to the reservoir.Over the past decades, POD-based reduced order modeling has been widely used in the literature [34,35].Additionally, POD is an excellent tool for gaining a solid understanding of the nature of the flows as the modes are optimal in terms of kinetic energy.In this respect, figure 3 shows the cumulative percentage of the kinetic energy up to each POD mode for the three arrays that are used in the current study.Evidently, the complexity of the flow increases as one moves from 2V2H to 4V2H and 6V2H arrays.For 2V2H and 4V2H arrays, around 83% and 71% of the kinetic energy is included in the first ten modes, while this number reduces to only 64% for the 6V2H array.For the current study, the first hundred POD modes are used for predictions via ESN, and this corresponds to 97%, 89%, and 92% of the total kinetic energy for 2V2H, 4V2H, and 6V2H arrays, respectively.Accordingly, it can be assumed that the main characteristics of the dynamics of the flows are well represented for all three arrays.In addition, the amount is relatively similar, so that an influence on the prediction capabilities might not be strong.It is important to note before moving on to the optimization of the ESN and predictions that the purpose of the study is to investigate the possibility of predicting and/or modeling such complex dynamics via ESN.Furthermore, it is essential to observe how the performance of the predictions is affected by the nature of the flow when one changes from a single bluff body flow to transient or co-shedding flow.Lastly, as already stated, ESNs are a primary example of a ML algorithm, and the effect of variations in its hyperparameters on its predictions in different chaotic dynamics can give valuable insight for further application and to design more complex ML algorithms to overcome its shortcomings as well.

Optimization
Network optimization or, in other words, identification of the best prediction and thus the best set of hyperparameters, is an essential part of ML applications in turbulence.In this regard, the application's objective and the task's complexity will be heavily influential.For the current study, it is clear that the inherent chaotic temporal evolution of the turbulent flow will prevent the network from having perfect deterministic predictions.Therefore, the aim is to maintain general statistical similarity in the predictions and search for the maximum possible deterministic accuracy.While deterministic accuracy translates to perfect alignment with the ground truth, the statistical similarity is more than just the proximity of the probability density function (PDF) estimates of the time coefficients.In other words, two signals can have identical PDFs and very different oscillations.Then, once the optimization parameter for a single time coefficient has been identified, one must figure out how to combine the parameter values of all hundred modes to come up with a final value for the entire prediction set.The solution for this is actually quite straightforward since the ratio of the total kinetic energy attributed to a certain POD mode to the total kinetic energy in the flow can be used as a weighting coefficient.
In this regard, one might suggest the use of conventional measures such as mean squared error in a weighted form, as already discussed.However, as shown in figure 13 in the appendix, such measures end up identifying minimum fluctuations around the average values as the best predictions since any offset in the phase or domain of the fluctuations in the predictions will quickly diverge to higher error measures.Purely statistical measures such as standard deviation (σ) or Kullback-Leibler divergence (KLD) of the entire values in the prediction phase have the drawback that they will not necessarily result in predictions with similar dynamics.Although the results of optimization based on the weighted minimum KLD in the appendix (figure 13) shows relatively promising outcomes, one needs to consider that it is always possible to have two signals with very different dynamics but similar PDF distributions, and the particular promising outcome of the KLD for the current example is mainly due to the very effective grasp of the general dynamic of the flow by the ESN and the absence of predictions with similar PDFs but very different dynamics such as dominant frequency.Furthermore, as already mentioned, the current study aims for predictions with maximum statistical and deterministic similarity simultaneously, and therefore, we choose to not rely on purely statistical measures such as KLD for optimization.It has been shown in the appendix in figure 14 that the optimization parameter of the current study WPS (which will be later explained) has similar performance in terms of statistical convergence to the ground truth compared to KLD, although it is not aiming for purely statistical convergence in its design.Finally, one might suggest the application of an extra ML algorithm for the task, which is indeed a very interesting direction to proceed but is in objection to the general claim of the current study to carry out predictions without using complex deep learning algorithms for the task.
Therefore, the Weighted Prediction Score (WPS) is introduced as the optimization parameter for the identification of the best predictions.For the calculation of the WPS, both the statistical similarity and deterministic alignment of the predictions are taken into account.Figure 4 shows the division of the time coefficients into positive (red), average (green), and negative (blue) values.Values in the range of [−σ/4, +σ/4] are considered average, and above and below are positive and negative values up to ±2σ.N is going to be the percentage of the time steps in which the predictions are aligned with the ground truth in terms of being in the range of these negative, average, and positive values.By defining the sign function as: ( and assuming the predicted values from the ESN for mode i as P = {p i1 , p i2 , . . ., p in } and the actual values as A = {a i1 , a i2 , . . ., a in }, then the N is calculated as: where 1{} is an indicator function that returns 1 if the condition inside is true (i.e. the sign of prediction matches the sign of the actual value) and 0 otherwise, and n is the total number of predictions.Once the N is calculated for each Mode i, then the WPS is determined as follows: Here, σ a i and σ a ′ i stand for the standard deviation of the Time Coefficient of the POD Modes (TCPM) and their temporal derivative of the ground truth, and σ p i and σ p ′ i are for the predicted counterparts, and E i is the energy percentage of the mode i.With this measure, while N i contributes to the deterministic alignment of the predictions with ground truth, the rest of the equation contributes to the statistical similarity.Finally, the multiplication by E i assures the proportional contribution of the WPS of the individual modes in the WPS of the entire prediction.Thus, when WPS of different predictions are compared, the one with the maximum value will be chosen as the best prediction.

Bluff body flow
Due to the proximity of the cylinders, the flow sees the array as a single bluff body, so one large vortex with a relatively low shedding frequency dominates the flow behind the array, as can be seen in figure 5.This flow shows the fastest convergence of the cumulative energy within the first POD modes (see figure 3).This slow and very unsteady dynamics is represented in the TCPM of the 1, 3, 5, 7, 9th modes in figure 5(a).The TCPM of the modes with lower energies, namely 20, 40, 60, 80, 100th, are shown in the appendix in figure 18.The gray background shows the part of the TCPM that are used for the reservoir training.In the figure, only the even modes are shown since most of the initial modes are associated with vortex shedding and thus accompanied by couples with 90-degree phase shifts.In addition, the spatial distributions of the modes are shown in figure 17 in the appendix for a clearer understanding of their role in flow dynamics.The red curves show the prediction results for the set with the best WPS values.While the slow and unsteady dynamics of the flow make accurate deterministic predictions of the TCPM quite challenging, the general dynamics of the oscillations are imitated to some extent in the predictions.This has been better represented in figure 15 in the appendix, where the dominant frequencies of the TCPM for the prediction and the ground truth are compared.Although these frequencies are not identical, they are from the same range.
Figure 5(b) shows the number of best predictions (with the highest WPS) for each of the 24 random realizations in relation to the SR, LR, and INS.The highlighted squares show the hyperparameter set with the best average prediction for all the random seeds, which is, in this case, SR = 0.1, LR = 0.3, and INS = 0.1.In figure 5(b) left, it is clear that lower LR and SR values result in better prediction.For INS vs LR in figure 5(b) middle, there is a large empty space for LR > 0.3 and INS < 7. Thus, this region is not suitable for predicting the dynamics of the flow behind 2V2H.Lastly, there is a quite linear relationship between INS and SR on the right side, such that as SR increases, INS should also increase for better predictions.
Figure 5(c) shows the reconstructed V y /V ∞ field for the best prediction vs the ground truth for the 350th and 700th time steps.The reason for showing only the vertical velocity field is that it represents much more important/characteristic information about the flow compared to horizontal or total velocity.This is because most oscillations occur vertically in the flow.From a deterministic point of view, the prediction possesses a degree of misalignment with the flow.However, their general shape and structure are well aligned.Considering the slow unsteady dynamics of the flow, as well as the fact that the forecasts are given in time steps far ahead in the predictions, they align reasonably well with the ground truth.A more detailed analysis in this regard will be followed in the discussion section.

Transient flow
Figure 6 shows the prediction results for the 4V2H array.This array represents the transient flow, characterized by a very asymmetric, unsteady flow where four individual vortex streets behind the array are observed.The top cylinder in the first column has a vortex street with a 24 Hz shedding frequency, while the three cylinders in the second column have vortex streets with frequencies of 13, 18, and 6 Hz from top to bottom, respectively.This large variation in the frequencies of the vortex streets, along with their very asymmetric structure (see figure 6(c)) provides very rich and chaotic dynamics for the ESN to encounter.
Figure 6(b) shows the predicted TCPM.Again, higher modes are shown in the appendix in figure 19.All predictions are very close to the ground truth.In other words, seeing the real oscillations in the gray part during the training, one can not choose the real values in the prediction phase if the prediction and ground truth were represented with the same color.From this perspective, mode 7 is particularly interesting.Even though it has a very different oscillatory nature compared to the other four, it has been well reconstructed by the ESN without being affected by the other modes' dynamics.However, deterministic alignment of the predictions with the ground truth occurs in different periods for each TCPM.For instance, while a 1 has deterministic alignment for the first 200 time steps, for a 3 and a 5 , the best alignment occurs between the 300th and 500th time steps.The dominant frequencies in figure 15 are relatively close and sometimes even identical for the first ten modes only modes 3, 4, and 5 show some deviations.From figure 17 it is apparent that these three modes represent features of the vortex street behind the central cylinder in the second column.In order to further investigate this difference, the frequency spectrum of the prediction and ground truth for mode 3 is shown in figure 16 exemplarily.Apparently, there are two peaks for the ground truth in the spectrum, one near 6 Hz and another which is stronger around 18 Hz.However, in the predictions the magnitude of the first peak has increased and the magnitude of the second peak has declined, so that two peaks with almost identical power are present.Therefore, one can conclude that unlike the first impression from figure 15 for dominant frequencies, actually, the frequency domain of the prediction is not that much different from the ground truth.
Next is the hyperparameters of the best predictions within each random seed that are shown in figure 6(b).The highlighted squares show the hyperparameter set with the best prediction for all the random seeds, which is, in this case, SR = 0.9, LR = 0.8, and INS = 15.Compared to the 2V2H array, the main difference is the predominant preference for high INS values for the best predictions.Then, for LR and SR, it is clear that low values are no longer preferred, and moderate to high values have, in general, better predictions.Furthermore, there is also a slight linear relation between SR and LR, meaning that as SR increases, it is necessary to increase the LR values as well in order to produce better predictions.These differences are all due to the significantly different dynamics of the flow.
Finally, figure 6(c) shows the reconstructed V y /V ∞ in comparison to the ground truth for the 350th and 700th time steps in the prediction phase.In general, the main features of the velocity field are well-preserved.However, there are slight shifts in the domain and phase of the values in the predictions.In conclusion, the ESN is capable of preserving the very diverse dynamics of the particular flow in 4V2H array even after 700 time steps.This suggests that the quality of the predictions does not degrade over time and is stable.

Co-shedding flow
The 6V2H array is the final example that is taken for the predictions in this study.For this case, due to the larger vertical distance between the cylinders, the upper and lower cylinders the first column, and the three cylinders in the second column, each produces its own individual vortex street behind the array.However, while the three vortex streets of the second column all have a shedding frequency of 15 Hz, the two vortex streets of the first column in the upper and lower edges have a higher frequency of 18 Hz.From the spatial distribution of the POD modes in the appendix in figure 17, it is apparent that other than the frequency, the synchronization of the vortex streets with each other is changing in time as well.The first and second modes show the time instants when the three vortex streets of the cylinders in the second column have phase synchronization, while the third and fourth modes show the time steps when the central cylinder in the second column oscillates with a 90-degree phase shift compared to the two neighboring vortex streets.This is very much clear in the ground truth TCPM of the a 1 and a 3 in figure 7(a), where the amplitude of oscillations in a 1 is only large when the amplitude is smaller for the a 3 and vise versa.One can see this in the predicted TCPM of the aforementioned modes as well.indicates that the ESN is not only able to reconstruct similar predictions with the ground truth but also grasps the relationship between the oscillations of the different modes and imitates it in its own predictions.
Figure 7(b) shows the hyperparameters of the best predictions for each random seed.The best prediction is for the set of hyperparameters with n = 5, INS = 7, SR = 0.8, and LR = 0.7, which is shown by the highlighted square.The general trend of the best hyperparameters is similar to the case of the 4V2H array.The faster dynamics of the 6V2H array seem to push more toward the higher INS values as the favorable choice.Moreover, higher LR and SR values become more predominant among the best predictions as well.Finally, the reconstructed V y /V ∞ by the ESN, along with their ground truth counterparts for the 350th and 700th time steps, are shown in figure 7(c).The predictions seem to have slightly exaggerated values in general.However, the general features of the flow dynamics are well-preserved, and there is no significant deviation from the ground truth long after predictions have started.

Discussion
Figure 8 shows the average WPS with respect to the hyperparameters for the 2V2H, 4V2H, and 6V2H arrays at the top, middle, and bottom, respectively.It should be noted that this is unlike the previous hyperparameter figures, where only the best predictions of random seed were shown.Generally, shifting between arrays results in clear differences in the average WPS distribution due to variations in the dynamics of the flow.Evidently, high LR values are not favorable for the slow dynamics of the 2V2H array at all.Moreover, compared to the other two arrays, 2V2H will not favor high INS values for the ESN as well.Again, this is related to the slow dynamics of its flow, which prevents the input signals from having a strong effect on the dynamics of the reservoir in the ESN for accurate prediction fidelity.Moving to the 4V2H and 6V2H arrays, it is evident that as the speed of the flow dynamics increases, higher INS, LR, and SR values result in better predictions.High INS values are, in particular, essential to prevent the oscillations from degrading as one moves further in the prediction phase.Whereas, high LR values are essential to keep up with the fast oscillations in the signals.
Figure 9 shows the average V y /V ∞ fields of the ground truth and predictions for the three arrays.Clearly, behind the array at the center for 2V2H, there are some exaggerated up and downward-directed flows in the predictions.This is due to the very unsteady and also slow dynamics of the flow for this case, which prevents the network from acquiring enough samples to produce predictions that align more accurately with the baseline average of the oscillations.In contrast, for the 4V2H and 6V2H arrays, the average fields align well with the ground truth even when they are significant asymmetric structures in the case of 4V2H.
While the average fields are very important in the assessment of prediction quality, the oscillations around these average values in the velocity field are as much important as well.In this regard, figure 10 shows the standard deviation of the values in the V y /V ∞ fields for the ground truth and prediction of all three arrays.Evidently, the predictions for all three arrays produce quite statistically well-aligned oscillations in the field, and standard deviation values are approximated well.This is due to the fact that for all three arrays, the ESN can produce oscillations that neither degrade nor over-exaggerate for the entire span of the prediction.Figure 11 shows the PDF of the V y /V ∞ fields in the predictions and ground truth for the three arrays.Apparently, there is a general agreement between the PDF values of the prediction and ground truth for all three arrays in the moderate ranges of V y /V ∞ .However, the ESN tends to produce slightly higher probabilities for more exaggerated velocity values in the field.With careful observation, this can be seen in the average fields in figure 9 as well.
Finally, in order to have a solid idea about the quality of the predictions from a deterministic point of view, figure 12 shows the correlation coefficient between the ground truth and predicted V y /V ∞ fields during the entire prediction phase for all three arrays.Clearly, the 2V2H array has higher correlation values in general, but this is due to the fact that in this array, a large portion of the flow field behind is simply part of the free stream and not included in the oscillations.Besides that, all three arrays have an interestingly oscillatory nature for their correlation curves.One can even suggest that the period of these oscillations tends to remain quite constant for the entire 700 time steps, especially for the 6V2H array.This shows that there might be a physical reason for the period of oscillations in the flow.One can find a clear similarity between the period of these oscillations and the period of large-scale vortex shedding in the case of single bluff body flow for the 2V2H array.This might suggest that there exists a background weak slow dynamics in the 4V2H and 6V2H arrays like there is in the 2V2H array as well, which somehow affects the oscillations of the individual vortex streets in the flow field.

Conclusion
This study attempted to predict the turbulent flow behind cylinder arrays by the application of ESNs.For this purpose, experimental data of the flow behind three arrays of seven cylinders are chosen.Due to the difference in the distance between the cylinders, these three arrays represent three different flow regimes of single bluff body, transient, and co-shedding flow, respectively.This provided the opportunity to examine the predictive potential of the ESN in different flow dynamics, which are all, in essence, vortex shedding behind bluff bodies.However, they have very different dynamics and degrees of complexity.The ESN algorithm was chosen for the current study due to its ability to handle time series and its simple structure, which supports the theory that a reservoir with random connections between neurons can learn actual physical information in chaotic dynamics.The flow data was reduced by POD and then the Time Coefficients of the POD Modes (TCPM)s were used for prediction by the ESN.The ESN consisted of 3000 neurons, and its main three hyperparameters, LR, SR, and INS, were varied in order to achieve the best performance.WPS of the predicted TCPM, which is a combination of deterministic and statistical convergence of the predictions to the ground truth, was chosen as the optimization parameter.
The best predictions were found to be for INS = 0.1, SR = 0.1, and LR = 0.3 for the bluff body wake in the 2V2H array, while the transient flow in the 4V2H array had the best outcome for INS = 15, SR = 0.9, and LR = 0.8, and finally for the co-shedding regime of the 6V2H array INS = 7, SR = 0.8, and LR = 0.7 had the best prediction.In general, for the slow and unsteady dynamics of the 2V2H, low LR values and relatively lower INS values had the best predictions.However, as the dynamics of the flow speed up in 4V2H and 6V2H arrays, better predictions tend to happen for higher LR and INS values.This was due to the necessity of faster reservoir update and greater input effect in order to cope with such fast and chaotic dynamics.
For all three arrays, the predictions were very close models of the oscillations in the ground truth, but the deterministic exact alignment was only partially possible.Moreover, the ESN also successfully imitated the interconnection between the oscillations of different modes.It was also interesting that the ESN was capable of stable prediction that could long continue without degradation or exaggeration in their values.The reconstructed velocity fields aligned well with the ground truth in terms of the average values, standard deviations, and PDF estimates.However, the deterministic alignment of the reconstructed velocity fields was periodically oscillating over time.
Last but not least, the results of the current study can be continued in several directions in the future.Given the quality of predictions by the ESN in the current study for such chaotic dynamics, ESNs and, more generally, reservoir computing can receive more attention from the community in comparison to much more complex algorithms that might not necessarily result in better outcomes.Moreover, the conclusions of the study regarding the best hyperparameters for different flow dynamics can provide a good foothold for designing more advanced algorithms that can overcome the shortcomings of reservoir computing.Finally, it could be very interesting to try to predict turbulent flow with a kind of forever-evolving dynamics that never repeats itself in the course of time.A very good example, perhaps, is the turbulent Rayleigh-Bénard convection, with its extremely complex and chaotic dynamics [36,37].

Figure 1 .
Figure 1.The water channel along with the cylinders and the field of view for the PIV measurements.The cylinders are arranged in nine different sets in total, with V = [2, 4, 6] mm and H = [2, 4, 8] mm, showing qualitatively different wake flow characteristics [29].

Figure 2 .
Figure 2. A schematic sketch of an Echo State Network (ESN).

Figure 3 .
Figure 3. Sum of kinetic energy percentage of the POD modes for 2V2H, 4V2H, and 6V2H arrays.

Figure 4 .
Figure 4. Schematic explanation of the division of the time coefficients into three main parts (positive, average, and negative) in Weighted Prediction Score (WPS) calculation.

Figure 5 .
Figure 5. Results of predictions for 2V2H array.

Figure 6 .
Figure 6.Results of predictions for 4V2H array.

Figure 7 .
Figure 7. Results of predictions for 6V2H array.

Figure 8 .
Figure 8.Average WPS values for the ESN predictions with respect to the hyperparameters for 2V2H, 4V2H, and 6V2H arrays at the top, middle, and bottom, respectively.

Figure 9 .
Figure 9. Average vertical velocity fields for the ground truth (top) and best prediction (bottom) for 2V2H, 4V2H, and 6V2H arrays at the left, middle, and right, respectively.

Figure 10 .
Figure 10.Standard deviation of the vertical velocity fields for the ground truth (top) and best prediction (bottom) for 2V2H, 4V2H, and 6V2H arrays at the left, middle, and right, respectively.

Figure 11 .
Figure 11.Probability density function of the best prediction vs the respective ground truth for 2V2H, 4V2H, and 6V2H arrays at the left, middle, and right, respectively.

Figure 12 .
Figure 12.The temporal fluctuation of the correlation values between the best prediction set of each array and their respective ground truth for the entire 700 prediction time steps.

Figure 16 .
Figure 16.The frequency spectrum of the ground truth and prediction for 3 of 4V2H array.

Figure 17 .
Figure 17.The first ten spatial POD modes for 2V2H, 4V2H, and 6V2H arrays at the left, middle, and right, respectively.

Figure 18 .
Figure 18.Ground truth and prediction of the temporal coefficients of modes 20, 40, 60, 80, 100 for 2V2H array.The predictions are for the network with the optimized set of hyperparameters with RS = 19, INS = 0.1, SR = 0.1, and LR = 0.3.The first half of the data with the gray background shows the time steps that are used for the training of the ESN.

Figure 19 .
Figure 19.Ground truth and prediction of the temporal coefficients of modes 20, 40, 60, 80, 100 for 4V2H array.The predictions are for the network with the optimized set of hyperparameters with RS = 21, INS = 15, SR = 0.9, and LR = 0.8.The first half of the data with the gray background shows the time steps that are used for the training of the ESN.

Figure 20 .
Figure 20.Ground truth and prediction of the temporal coefficients of modes 20, 40, 60, 80, 100 for 6V2H array.The predictions are for the network with the optimized set of hyperparameters with RS = 5, INS = 7, SR = 0.8, and LR = 0.7.The first half of the data with the gray background shows the time steps that are used for the training of the ESN.