Neural Network-Based Model Reduction of Hydrodynamics Forces on an Airfoil

: In this paper, an artiﬁcial neural network (ANN)-based reduced order model (ROM) is developed for the hydrodynamics forces on an airfoil immersed in the ﬂow ﬁeld at different angles of attack. The proper orthogonal decomposition (POD) of the ﬂow ﬁeld data is employed to obtain pressure modes and the temporal coefﬁcients. These temporal pressure coefﬁcients are used to train the ANN using data from three different angles of attack. The trained network then takes the value of angle of attack (AOA) and past POD coefﬁcients as an input and predicts the future temporal coefﬁcients. We also decompose the surface pressure modes into lift and drag components. These surface pressure modes are then employed to calculate the pressure component of lift C pL and drag C pD coefﬁcients. The train model is then tested on the in-sample data and out-of-sample data. The results show good agreement with the true numerical data, thus validating the neural network based model.


Introduction
Airfoils are cross sectional curved surfaces which give required ratio of hydrodynamic forces (i.e., lift and drag). Flow separation is an active area of research in airfoils due to stall phenomenon when the angle of attack goes beyond a certain limit. In response, there is vivid decrease in the hydrodynamic forces (i.e., lift and drag) [1,2], so there is a need to control flow separation.
Computational fluid dynamics (CFD) is widely used for numerical solution of complex, nonlinear, and transient flows to study the flow separation phenomenon. To further develop flow control strategies, reduced-order models (ROM) are often developed. Proper orthogonal decomposition (POD) [3] is mostly used to project these transient flows on the low dimensional subspace where these flows behave linearly. POD was first used for reduced-order modeling of turbulence by projecting the Navier-Stokes equations onto the low dimensional space [4]. Different model reduction techniques have been developed in recent years to reduce the dimensionality of dataset, such as POD, dynamic mode decomposition (DMD) [5,6], and its variants [7][8][9], spectral proper orthogonal decomposition (SPOD) [10,11] and machine learning-based methods, such as autoencoders [12], etc. Taira et al. [13,14] did the detail review of such methods and techniques in their paper. POD/Galerkin models can work well but often require careful tuning to give satisfactory results. In addition, POD-based models representing complex (nonlinear) dynamics are known to be inaccurate for long-time integration [15], and different closure model strategies were employed in literature [15][16][17][18][19][20] which were aimed at modeling small-scale dynamics. In addition, if DMD is used for system identification, it results in a linear system [13,14]. For nonlinear systems, the Koopman operator, its eigenfunctions, and corresponding Koopman modes may be computed using a version of DMD [7][8][9]21]. These eigenfunctions can then be used to determine coordinates [22,23], in which the system is linear, and the resulting Koopman modes may be used to understand coherent structures that occur at a particular temporal frequency.
Reduced order models are broadly used for modeling of fluid systems [14] as they reduce the CPU time to achieve the desired results and preserve the physics of interest. There are mainly three approaches which are used to develop reduce order models: (a) The first is Phenomenon-based [24,25]. (b) The second is Direct approach, which is the intrusive way of modeling the fluid systems. This approach is dependent on the source code of the governing equations. POD-Galerkin-based models suffer from instability and efficiency issues [15,18,19,[26][27][28][29][30]. Other intrusive models which show promising results can be seen in references [31][32][33][34]. (c) The third is Non-intrusive ROM which do not require modification of the source code rather these are purely data-driven. Xiao et al. presented such reduced order models for Navier-Stokes equations [30,35]. Machine learning has also been explored to infer POD coefficients in non-intrusive parametric ROMs which can be seen in references [36][37][38][39]. San and Maulik [40], San et al. [41] used artificial neural networks (ANN) to predict the temporal coefficients of POD in non-intrusive way. Wang et al. [42] used long-short term memory networks (LSTM) to predict the future time steps of POD coefficients. ANN and LSTM have also been explored for developing the closure for truncated POD modes [40,43]. Murata et al. [44], Fukami et al. [45] used convolutional neural networks (CNN) for data compression and reconstruction of high dimensional flow fields. These models require a large amount of data for training, which also increases the computational expense [44][45][46][47][48][49].
Nonetheless, most ROM based on POD employs velocity field information for modeling [41,50], whereas pressure specifically on the surface also plays an important role in the accuracy of POD-based ROMs [51]. Many studies have been done on the ROM based on pressure mode decomposition (PMD) that also enables accurate reproduction of lift and drag forces (hydrodynamic forces) [15,52].
To understand the flow physics over an airfoil under unsteady motion, it is better to have a fundamental knowledge of physics under steady conditions. Flows at a low Reynolds number over the static airfoil at different angles of attack can also be used to understand the flapping phenomenon [53]. In this study, we simulate the flow past a NACA-0012 airfoil at Re = 1000 at different angles of attack (AOA = 10 • , 15 • , 20 • ). Since the pressure plays an important role in the development of the hydrodynamic forces [15,52] and the wake of the vortex street behind the structure, we numerically calculate the pressure field in the present study. We then perform the POD analysis and investigate the pressure distribution over NACA0012 airfoil. Using these surface POD modes, we calculate the hydrodynamic forces. In literature, many researchers used "National Advisory Committee for Aeronautics (NACA)" airfoils for investigation of aerodynamic/hydrodynamic forces. For the validation purpose of our DNS technique, we also considered the NACA 0012. However, the proposed neural network model depends on the input data which are the POD temporal coefficients but not on the airfoil shape. The same neural network setting can be employed to the data of other NACA airfoils.
The procedure to perform the pressure mode decomposition (PMD) analysis consists of the following steps: First, the direct numerical simulation (DNS) of the incompressible flow field data over the airfoil is calculated for a long enough duration to mitigate the effects of the transient phase. Second, the snapshots of the pressure field at different time steps are assembled in a matrix that characterizes the coherent structures in the flow field. Third, the finite numbers of dominant modes or eigenfunctions, that optimally spanned the low-dimensional solution space of the pressure field, are extracted on the airfoil surface and decomposed into lift and drag components.

Governing Equations
The flow under consideration is a two-dimensional incompressible unsteady flow over NACA0012 airfoil at different AOA in the flow field. The governing equations include the continuity and momentum (Navier-Stokes) equations. The two-dimensional nondimensional form of these governing equations are shown in Equations (1) and (2), as follows: where the u i , p, and ρ represent the Cartesian components of velocity, pressure, and density of the fluid, respectively. We use the length of the foil, denoted as c, and uniform velocity, labeled as U ∞ , as the reference scales for the length and velocity, respectively, in order to define the Reynolds number (Re = ρU ∞ c/µ), where µ is dynamic viscosity of the fluid.
In this study, we used a non-staggered O-type grid to simulate the flow around the airfoil, as shown in Figure 1 with the representation of inflow and outflow boundaries. All the state variables u,v,p are defined at the cell centers; however, the volume fluxes are defined at the midpoint of the corresponding faces. To handle the curved form of the physical domain, we transform Equations (1) and (2) from the Cartesian coordinate system (x, y) into a curvilinear coordinate system (ξ, η) to accurately handle curved grid lines. Consequently, these equations takes the following forms: where the terms for fluxes and geometric matrices are defined as:

Discretization Scheme
A fractional step method is used by splitting the momentum equations into an advection-diffusion equations (by omitting the pressure term) and the pressure Poisson equation. The pressure Poisson equation is then obtained by taking the divergence of the momentum equations and subjected to satisfy the continuity equation. This splitting method leads to a three-step predictor-corrector method. First, the intermediate or predicted velocity field is computed by solving the advection-diffusion equations without satisfying the continuity equation. Secondly, the pressure Poisson equation is solved to determine the pressure field, which satisfies the continuity equation. Finally, the divergencefree velocity field is computed by correcting the predicted velocity field. Note that we generate the orthogonal grid between the airfoil and circular outer domain of radius 30c. We numerically solved the elliptic equations as used by Ryskin and Leal [54] to compute the grid nodes.
The differential operators involved in the governing equations are discretized using the finite-difference approximation. A second-order accurate centered difference scheme is employed to discretize all the spatial derivatives, except the convective terms. The convective terms are discretized using a variation of QUICK method [55]. The velocity variables (u i ) are calculated at the faces from the nodal values using quadratic upwinding interpolation. The positive and negative fluxes ( ), respectively, are computed, and a generic stencil is used to perform upwinding of the QUICK method. For marching in time, a semi-implicit scheme is used in which the diagonal viscous terms are advanced using the second-order accurate Crank-Nicolson method, and all the other terms are advanced using the second-order accurate Adams-Bashforth method. The Adams-Bashforth scheme is chosen because of its computational efficiency when coupled with the fractional step method. This methodology is similar to that employed by Zang [56]. The only difference is that we employed the line successive over-relaxation (LSOR) scheme to solve the advection-diffusion equation and the bi-conjugate gradient stabilized (BiCGStab) method to solve the pressure Poisson equation.

Grid Generation
A two-dimensional orthogonal grid around the symmetric airfoil is generated using a body-fitted grid generation technique as proposed by Thompson et al. [57]. A finite-difference method is used to solve the system of coupled elliptic partial differential equations [54,[57][58][59][60][61][62]. The governing equations for the generation of the two-dimensional orthogonal body-fitted grid, as used by Ryskin and Leal [54], are given by ∂ ∂ξ f ∂y ∂ξ where x and y are the coordinates in the physical plane, ξ and η are the coordinates in the transformed plane (see Figure 2), and f is the distortion function which is the ratio of scale factors in the η-direction to that of the ξ-direction, i.e., The scale factors are defined as follows: After some manipulations and rearrangements, Equations (7) and (8) yield A central-difference scheme at the interior points is employed to discretize the partial differential equations. All the scalar factors and the distortion function are calculated from the previous approximation of the solution, leading to a system of linear algebraic equations. The tri-diagonal matrix algorithm (TDMA), which is an iterative technique, is used to obtain the solution until the relative error of the distortion function converges to the desired accuracy [63].

Hydrodynamic Forces
The pressure and shear stresses are the two main stimuli through which a fluid exerts forces on a body immersed in it. These forces can be decomposed into lift and drag components. We use the free-stream dynamic pressure ρU 2 ∞ /2 and the length of the foil to define the dimensionless lift and drag coefficients, which are expressed as: where ω z = ∂v ∂x − ∂u ∂y is the spanwise vorticity component. Note that the terms for shear stress in Equation (14) diminish as the Reynolds number increases; therefore, the pressure is considered as the main source in the development of C l and C d . Figure 3 shows the time histories of pressure and viscous shear stresses components of hydrodynamic forces lift and drag at AOA = 20 • . It is noted that approximately 98% of the lift and 88% of the drag are generated because of the pressure forces.

Validation of the Numerical Simulation
To validate the solution methodology, NACA0012 airfoil subjected to the incompressible flow field at different angles of incidence is investigated at Re = 1000. A mesh of size 361 × 257 is used in this study having 257 mesh points on the airfoil surface and 361 mesh points in the radial direction with an outer domain diameter of 50D. We simulate the flow for long enough time to reach the steady-state solution. Figure 4 shows the time-averaged lift coefficient (C l ) distribution versus angle of attack (α) obtained in the present study, along with the results computed at the same Reynolds number by Kurtulus [64] for 0 • < α ≤ 90 • , Liu et al. [65] and Khalid and Akhtar [66] for α < 30 • , Suzuki et al. [67] at α = 15 • , HOA-RAU et al. [68] at α = 20 • , and Mittal and Tezduyar [69] at α = 10 • . The numerical results are in good agreement with the published data, hence validating the present simulations. The time histories of C l and C d at α = 20 • are shown in Figure 5a, and the power spectrum of C l is shown in Figure 5b, along with the results of Khalid and Akhtar [66]. The fast Fourier transformation (FFT) algorithm is used to calculate the amplitude spectrum for the lift coefficient. Figure 5b shows that the lift spectrum has a dominant peak at f s = 0.541, and it contains both the odd and even harmonics, which indicates the existence of a quadratic, as well as a cubic, nonlinearity in the system, as observed by Khalid and Akhtar [66].

Parallel Implementation
CFD of the fluid-structure problems on large mesh sizes generally involves several million degrees of freedom within each time loop to simulate the unsteady flow behavior. The demand for high-order accuracy and dealing with complex geometries increases the computational workload. Therefore, to reduce the overall computational time, we employ parallel computing of the incompressible flow over the airfoils. In the parallel implementation, the problem domain is partitioned to allocate its regions to different processing units to reduce overall execution time. These processing units simultaneously compute the results on their assigned part of the domain. For inter-process communication among the processors, we utilize the Massage Passing Interface (MPI) library. In Figure 6, a computational domain decomposed in 8 processing units is shown. The different colors represent the part of the domain to be assigned to different processing units. For a detailed implementation of the parallel algorithm, the interested readers are referred to Reference [70].

Artificial Neural Network
In the current work, feed forward neural network has been employed to map the past time steps to future time steps. The inputs to the neural network are the POD temporal coefficients of the pressure and velocity fields of airfoils along with the angle of attack. The basic architecture used is shown in Figure 7. The inputs are then passed through two hidden layers, and each layer has 512 neurons followed by relu activation function. For simplicity, only five neurons are shown in Figure 7. Tensorflow [71] is used for the training of neural network, which minimizes the error between predicted and the target quantities of interest (i.e., temporal coefficients) to determine the best fit parameters (weights and biases). For optimization of neural network a stochastic optimization algorithm, Adam [72] was employed. Autokeras [73] was used to find the best possible ANN architecture (number of hidden layers and neurons), which uses the Bayesian Optimization [74]. Here, 370 time steps from each case of AOA = 10 • , 15 • , 20 • have been used to train the network, and 80 time steps were set aside for the testing purpose. Time coefficients for the case of AOA = 12 • are not used for the training, but these are predicted by trained ANN. After training, the network can be used as a function, which takes AOA and past time steps to predict the future time steps, as shown in Equation (16).
where X is the input matrix that contains the temporal coefficients and angles of attack (AOA) values. It has one column per coefficient. The weight matrix W contains all the weights, except for the ones from the bias neuron. The bias vector b contains all the connection weights between the bias neuron and the artificial neurons. It has one bias term per artificial neuron. The function ψ is called the activation function, which is 'relu' in this study.

Proper Orthogonal Decomposition (POD)
POD modes can be computed through either Singular Value Decomposition (SVD) or the method of snapshots [75,76]. In this study, we employ the latter technique since it is computationally inexpensive as compared to SVD for large mesh sizes used in high-fidelity simulations. The set {U n = U(x, t n ), t n = (n − 1) t, n = 1, 2, . . ., S}, where S represents the number of snapshots of time discrete snapshots and is used as the input for POD of the state vector U = (u, v, p) T over one time-period T of vortex shedding. The unknown, U, is decomposed into its mean part,Ū(x), and fluctuatingÚ(x, t) part. Next,Ú(x, t) is expanded in terms of the linear combination of temporal coefficients, which need to be determined, and POD bases as: where M represents the finite number of modes in the expansion, Φ j (x) denotes POD-modes or simple modes of U, and α j (t) are temporal coefficients. Moreover, Φ j (x) = (ψ u j , ψ v j , φ i ) T , where ψ u j and ψ v j are the u and v velocity modes, respectively, and φ i are the pressure modes. In addition, we have α j (t) = (q u j , q v j , a j ), where q u j , q v j , and a j are the temporal coefficients for u, v, and p, respectively.
The decomposition of the pressure field in terms of its meanp(x) and the fluctuatinǵ p(x, t) parts is carried out by: The instantaneous solutions at time t 1 , t 2 , · · · , t s ∈ [0, T] of U are stored in matrix P of size N × S, where S N and N denotes the number of grid points. The objective is to seek a low dimensional basis {Φ 1 , Φ 2 , · · · , Φ M } that satisfies the following optimization problem.
subject to the condition of orthogonality, Φ i , Φ j H = δ ij , where δ ij is the Kronecker delta function, and H is a real Hilbert space, such thatÚ(., t) ∈ H|t ∈ [0, T]. In order to solve Equation (19), we consider the following eigenvalue problem: where v k for k = 1, 2,· · ·, S are eigenvectors, λ 1 ≥ λ 2 ≥ · · · ≥ λ S > 0 are eigenvalues, and P T P = M ∈ R S×S is the correlation matrix that can be defined as: It can be shown [76] that the solution of Equation (19) is given by: For our present work, we utilize 50 snapshots of flow fields captured in one timeperiod (T = 1/ f s ) and assemble them in P. For a pressure field, the correlation matrix M of the snapshots is defined as:

Results & Discussion
Using the pressure POD data, the surface mean pressure (p s ) and the surface pressure modes (ψ s j ) are extracted from the pressure modes ψ j . Symbolically, the sine and cosine components of ψ s k over the airfoil surface are, respectively, denoted by −ψ s k sin Γ and −ψ s k cos Γ, which contribute in the lift and drag coefficients, respectively. Likewise, the distribution of the mean pressurep over the airfoil surface contributes to the mean lift and drag coefficients denoted by −p s sin Γ and −p s cos Γ, respectively.

Hydrodynamics Forces and Wake
In Figure 8a, instantaneous lift coefficient C l is shown at AOA=10 • , 15 • , and 20 • . The lift coefficient at AOA = 20 • is relatively higher in magnitude but oscillating about its nonzero mean value with low frequency as compared to that of the other two cases. The power spectra of C l of the simulated cases are shown in Figure 8b. A snapshot of vorticity pattern in the wake of the inclined airfoils is shown in Figure 9. The shear layer for the case AOA = 10 • is close to the surface, while, for the cases 15 • and AOA = 20 • , the shear layer is clearly separated. importantly, we find that the first four pressure modes constitute over 95% of the total contribution of the system (see Figure 10b). This observation demonstrates the effectiveness of the modal decomposition based on pressure data of flow fields around the airfoil.

PMD Analysis
We present the mean mode (Mode 0), defined as the average of the snapshots, in Figure 11 (left column), along with the datap s on the airfoil surface. The decomposition of p s in its sine and cosine components is also depicted. Here, we represent the upper and lower surfaces of the airfoil by chord length from 1 to 0 and 0 to 1, respectively. For plots in the right column of Figure 11, the 1s represents the trailing edge, whereas 0 represents the leading edge of the airfoil. In each case, we observed higher pressure near the edges of the airfoil and its magnitude increases with AOA. The pressure distribution over the upper surface and near vicinity of the lower leading edge is the main source to the generation of mean lift coefficient, while the pressure at the edges of the airfoil contributes to the mean drag coefficient.

Lift and Drag Models
The ROM for the lift and drag coefficients can be obtained by replacing the pressure term (p) in Equation (14) with the pressure modal expansion given in Equation (18) over the airfoil surface (φ s k ). As mentioned earlier, the viscous contribution is neglected because of its insignificance at high Reynolds numbers. The resulting forms of the lift and drag decomposition models are given as: where a i (t) is the temporal coefficient for both lift and drag models, and M is the number of modes in the expansion. The terms L i and D i are scalars called lift decomposition coefficient (LDC) and drag decomposition coefficient (DDC), respectively, and contribute, respectively, to the lift and drag coefficients as the weight of the temporal coefficient. Thus, the decomposition models are independent of any spatial distribution, unlike in a typical POD-ROM framework. The LDC and DDC are defined as the integral of sine and cosine components of surface pressure modes: Here, the subscripts • and k represent the decomposition ofp s and ψ s k , respectively. The LDC and DDC of mean and the first ten modes of each simulated case are plotted in Figure 13. It can be seen that contribution of the mean mode is significant to the lift and drag, and its significance in terms of magnitude increased with AOA. Figure 13a,b show that mean plus first modes are enough to model the lift and drag coefficients of AOA = 10 • using model (24). In the case of AOA = 15 • , the mean plus the first four modes are sufficient for lift and drag models, while, in the case of AOA = 20 • , as we can infer from Figure 13c,d, the LDC and DDC of first eight modes are significant, which would contribute to lift and drag models.

ANN-Based ROM Using PMD
ANN is trained on the POD coefficients obtained at three AOA = (10 • , 15 • , 20 • ), and out of sample POD coefficients of AOA = 12 • are predicted. The first four pressure POD coefficients are shown in Figure 14, indicating good agreement with the true data. We then perform the PMD analysis and predict the lift coefficients for both training data and out of sample case. Figure 15 shows the predicted and true lift coefficients for all AOA (10 • , 15 • , 20 • , 12 • ).

Conclusions
In this study, we simulated the incompressible flow over a stationary airfoil NACA0012 at different AOA. The main objective of the study was to develop the ROM for hydrodynamic forces acting on the airfoil surface using machine learning tools. Unlike the traditional POD-ROM-based models, we developed the ANN-ROM using POD temporal coefficients and PMD. We took flow field snapshot data at AOA = 10 • , 15 • , 20 • . We then performed POD to obtain optimal modes and temporal coefficients. Then, we trained the neural network using these temporal coefficients. ANN was used as a function which takes past time steps and angle of attack to predict the future time steps. We also performed interpolation of out of sample data at angle of attack 12 • . We used these predicted temporal coefficients to calculate lift coefficients using PMD. We employed pressure in this study because it plays an important role in the accuracy of POD-based ROM, as indicated by Reference [15,51,52]. The results show that the ANN model is more accurate than the traditional POD-ROM in terms of predicting in-sample and out-of-sample POD coefficients. Furthermore, the pressure component of lift is reconstructed using PMD analysis and shows excellent agreement, which is found with the true data.