Dynamic stall on airfoil is of great importance in engineering applications. In the present work, Fourier neural operator (FNO) is applied to predict flow fields during the dynamic stall process of the NACA0012 airfoil. Two cases with different angles of attack are simulated by Reynolds averaged numerical simulation with the Spalart–Allmaras (SA) model at $ R e = 4 \xd7 10 4$. A prediction model is directly constructed between the flow fields at several previous time nodes and that at the future time node by FNO. The prediction of sequence flow fields based on the iterative prediction strategy is achieved for the dynamic stall. The results show that FNO can achieve a fast and accurate prediction of streamwise velocity, normal velocity, pressure, and vorticity for both cases. The dynamics of vortices around the airfoil is analyzed to demonstrate the prediction accuracy of FNO. In addition, FNOs with different configurations are tested to achieve a lower error and a shorter training time-consuming.

## I. INTRODUCTION

Dynamic stall on airfoils is an important issue in unsteady aerodynamics, which involves complex flow phenomena including flow separation, forming, and shedding of vortices. It occurs on the airfoil undergoing large-amplitude pitching motions with the change of angle of attack that leads to dramatic changes in aerodynamic loads.^{1} Due to its complexity and importance, a fast and accurate prediction of dynamic stall flow fields is required for the design of rotorcraft forward-flight, and wind energy. In the past few years, many experiments and simulations have been studied to analyze the effect of dynamic stall.

Until now, many experiments have been carried out with pitching airfoils. Piziali^{2} conducted a comprehensive experimental investigation of the pressure distribution over the NACA0015 airfoil undergoing pitching motions. Baik *et al.*^{3} conducted an experimental study of the SD7003 airfoil and a flat plate model on a nominally 2D wing undergoing pitching and plunging motion. The phase-averaged 2D instantaneous flow field was captured for selected phases of the motion using the particle image velocimetry (PIV). Favier *et al.*^{4,5} preliminarily investigated the unsteady flow features of large-amplitude fluctuations of both velocity and incidence on the NACA0012 airfoil by an experimental method. Greenblatt *et al.*^{6} examined NACA0012 and NACA0018 pitching airfoils in an oscillatory relative free-stream in a water tunnel with an electric motion rig to oscillate the test article and in a wind tunnel with moving louvers that varied flow speed in the test section. Medina *et al.*^{7} further investigated the high-amplitude combined pitching/oscillatory motions of the NACA0012 airfoil to explain trends in the lift coefficient history through flow visualization and surface pressure measurement.

In the past few decades, computational fluid dynamics (CFD) has been widely applied to many engineering problems. However, high-fidelity numerical simulation is expensive for the dynamic stall problem due to the brand range of simulation time and space scale. For this reason, most of the investigations about dynamic stall simulation were completed using Reynolds averaged numerical simulation (RANS) instead of large eddy simulation (LES) and direct numerical simulation (DNS). The flow around the NACA0015 airfoil at static and dynamic stall using RANS and detached eddy simulation (DES) was simulated by Szydlowski and Costes.^{8} Liu *et al.*^{9} proposed a new delayed detached eddy simulation method with an adaptive coefficient (DDES-AC) for separation and dynamic stall flow past the NACA0015 airfoil. DDES-AC improved the prediction accuracy and agreed excellently with the experimental data than the original DDES and RANS. Rosti *et al.*^{10} used the DNS method to simulate the flow around the NACA0020 airfoil at $ R e = 2 \xd7 10 4$ undergoing a ramp-up motion. Visbal and Garmann^{11} analyzed the unsteady separation and dynamic stall vortex formation over a constant-rate pitching NACA0012 airfoil. The freestream Mach number was *M* = 0.1, and the Reynolds number was $ R e = 2 \xd7 10 5$. Gharali and Johnson^{12} used the SST $ \kappa \u2212 \omega $ model to investigate the effects of horizontal oscillations of the freestream velocity superimposed on a pitch oscillating NACA0012 airfoil for $ R e \u2248 10 5$. Wang and Zhao^{13} accomplished the numerical investigations of the unsteady aerodynamic characteristics of the SC1095 airfoil under conditions of coupled freestream velocity/pitching oscillation using the SA model.

In recent years, many machine learning methods have been widely used in fluid dynamic problems,^{14,15} including turbulent models,^{16–18} hypersonic flow transition,^{19,20} and flow field prediction.^{21–23} Wang *et al.*^{24} proposed a symbiosis method through the experimental data and numerical simulation to predict the lift and moment coefficients of dynamic stall under different balanced parameters. Tracey *et al.*^{25} presented a new method for developing the functional form of the SA model by a machine learning method. A trained neural network successfully reproduced the SA model in a wide variety of flow conditions from 2D flat plate boundary layers to 3D transonic wings. Singh *et al.*^{16} embedded a data-driven framework comprising full-field inversion and machine learning in a traditional RANS solver to improve the applicability of the SA model to strong adverse pressure gradients in the flow past airfoils. Volpiani *et al.*^{26} demonstrated an off-line method for parameter inference and machine learning for turbulence modeling. A neural network was trained to reconstruct the volume force correction to the SA model from variational data assimilation for RANS computations. Bin *et al.*^{27} used experimental and simulation data to recalibrate the SA model by Bayesian optimization and neural networks in free-shear flows, in viscous sublayers, in buffer layers in the presence of adverse pressure gradients. Jäckel^{28} leveraged the field inversion and machine learning approach to obtain a closed-form correction for the SA model to improve predictions of separated flows. Moreover, the additional physical knowledge has been incorporated into deep learning methods to solve differential equations.^{29–32} Raissi *et al.*^{33} employed physical informed neural networks (PINNs) to solve nonlinear partial differential equations. Jin *et al.*^{34} developed the Navier–Stokes flow nets (NSFnets) to simulate incompressible laminar and turbulent flows. Eivazi *et al.*^{35} used PINNs to simulate incompressible turbulent flows by solving the Reynolds averaged Navier–Stokes equations. Mi *et al.*^{36} proposed a Cascade-Net model to predict the fluctuating velocity in the near wall wake of a circular cylinder.

Most neural network models are good at learning mappings between finite-dimensional Euclidean spaces. The generalization of these models is limited by different parameters, initial conditions, or boundary conditions of differential equations.^{32,37–40} Li *et al.*^{41} proposed an FNO framework that can directly learn the mapping function between infinite dimensional spaces from a series of input–output pairs. It is up to three orders faster than traditional computational fluid methods and can model turbulent flows with zero-shot super-resolution successfully. Moreover, the FNO framework for general geometries has been proposed and applied to flow around an airfoil, pipe flow, elasticity, and plasticity problems.^{42} Deng *et al.*^{43} used Unet and FNO to establish the relationship between flow at previous time steps and that at future time steps. A laminar flow around a cylinder at *Re* = 200 and a transonic buffet of a supercritical airfoil at *Ma* = 0.73 were used to demonstrate the effectiveness and robustness of Unet and FNO. Li *et al.*^{44} developed an FNO approach to LES of three-dimensional turbulence. Filtered DNS flow fields of isotropic turbulence at different times were used for training FNO models. The results showed FNO model can perform better than the dynamic Smagorinsky model and the dynamic mixed model in the prediction of velocity spectrum, instantaneous flow structures, and PDFs of velocity and vorticity increments. Peng *et al.*^{45} proposed an attention-enhanced neural network method to predict the nonequilibrium feature of turbulence. A linear attention coupled FNO was further developed to simulate 3D isotropic turbulence and free shear turbulence.^{46} Wen *et al.*^{47} presented the U-net enhanced FNO (UFNO) model to solve complex CO_{2}–water multiphase flow problems. You *et al.*^{48} developed the implicit FNO (IFNO) model to model the increment between layers as an integral operator to capture the long-range dependencies in the feature space. Li *et al.*^{49} proposed an implicit U-Net enhanced FNO (IUFNO) model for the prediction on the long term large-scale dynamics of turbulence. Compared to recurrent neural network (RNN) and long short-term-memory (LSTM) methods, FNO has an advantage in efficiently learning the multiscale dynamics and spatially nonlocal interactions of complex flow structures based on Fourier transform.

This research aims to realize the fast flow field prediction of dynamic stall on a pitching airfoil using Fourier neural operator. The rest of the paper is organized as follows: Sec. II introduces the method of FNO for establishing the mapping between the flow fields at several previous time nodes and that at the future time node. In Sec. III, a mesh of the NACA0012 airfoil is chosen with the balance between computational accuracy and neural network training efficiency. Two numerical unsteady flow datasets with different configurations are used to test the performance of FNO. In Sec. IV, the results show that FNO can achieve a fast and accurate prediction of dynamic stall under different AOAs. The prediction of sequence flow fields is achieved based on the iterative prediction strategy. A deeper analysis of the parameters of FNO is conducted to investigate the performance comprehensively. Conclusions are finally drawn in Sec. V.

## II. FOURIER NEURAL OPERATOR

^{41}Let $ G * : A \u2192 U$ be a non-linear map where $ A = A ( D ; R d a )$ and $ U = U ( D ; R d u )$ are separable Banach spaces of function taking values in $ R d a$ and $ R d u$. Here, $ D \u2282 R d$ is a bounded, open set.

*R*is a real number space, and $ R d a$ and $ R d u$ are the value sets of input

*a*(

*x*) and output

*u*(

*x*).

*d*,

_{a}*d*, and $ d \phi $ are the dimensions of

_{u}*a*(

*x*),

*u*(

*x*), and $ \phi t ( x )$. We aim to learn the non-linear map $ G *$ by constructing a parametric map $ G : A \xd7 \Theta \u2192 U$. FNO needs to find the optimal parameter $ \theta \u2020 \u2208 \Theta $ by data-driven methods. The input

*a*(

*x*) is first transformed to a higher dimensional representation $ \phi 0 = P ( a ( x ) )$ by a shallow fully connected layer

*P*. $ \phi 0$ is updated by Eq. (1) $ \phi t \u2192 \phi t + 1$,

^{41}

Fourier integral operator $ K : A \xd7 \Theta K \u2192 L ( u ( D ; R d \phi ) , u ( D ; R d \phi ) )$ maps to the bounded linear operator on $ u ( D ; R d \phi )$ and is parameterized by $ \varphi \u2208 \Theta K$. *W* is a linear transformation, *σ* is the nonlinear activation function, and $ D \u2208 R d$ is a bounded, open set. Equation (1) shows the iterative architecture, $ \phi 0 \u2192 \phi 1 \u2192 \u2026 \phi M$. Finally, the output *u*(*x*) is projected from $ \phi M$ by applying the local transformation $ u ( x ) = Q ( \phi M )$, where $ Q : R d \phi \u2192 R d u$ is a shallow fully connected layer. Here, $ R d \phi $ and $ R d u$ are the real number sets of $ \phi M ( x )$ and *u*(*x*).

^{41}

*F*is the Fourier transform of a function, and $ F \u2212 1$ is its inverse transform. $ R \varphi $ is the Fourier transform of a periodic function

*κ*parameterized by $ \varphi \u2208 \Theta K$.

*D*with $ n \u2208 N$ points, $ \phi t \u2208 R n \xd7 d \phi $. $ F ( \phi t ) \u2208 C k max \xd7 d \phi $ is obtained by truncating the higher modes. Here,

*C*denotes the complex space. Multiplication by the weight tensor $ R \varphi \u2208 C k max \xd7 d \phi \xd7 d \phi $ gives rise to the following equation:

*F*. For $ f \u2208 R n \xd7 d \phi , l = 1 , 2 , \u2026 , d \phi , k = ( k 1 , k 2 , \u2026 , k d ) \u2208 Z s 1 \xd7 Z s 2 \xd7 \cdots \xd7 Z s d$, and $ x = ( x 1 , x 2 , \u2026 , x d ) \u2208 D$, the FFT $ F \u0302$ and its inverse $ F \u0302 \u2212 1$ are defined as

^{42}Let

*D*be the physical domain, and

_{a}*D*be the computational domain. Let $ x \u2208 D a$ and $ \xi \u2208 D c$ be the corresponding mesh points. Denote the input meshes as $ T i = x ( i ) \u2282 D a$ and the computational meshes as $ T c = \xi ( i ) \u2282 D c$. For simplicity, we assume the output mesh is the same as the input mesh. The points in the computational mesh are transformed to the physical mesh by a coordinate transformation $ \varphi a$,

^{c}^{42}

*T*is non-uniform. The input mesh can be viewed as a multi-dimensional array $ T i [ i 1 , \u2026 , i d ]$, and then, its indexing induces a canonical coordinate map,

^{i}^{42}

Especially, since $ ( i 1 / s 1 , \u2026 , i d / s d )$ is a uniform mesh, such FFT can be directly used in this case.

## III. COMPUTATIONAL SETUP

*α*, of the NACA0012 airfoil is changed in the form of sin function as Eq. (6),

^{12}

*α*

_{0}, $ \Delta \alpha $, and

*f*represent the mean AOA, pitch oscillation amplitude, and pitch oscillation frequency, respectively. The oscillation system has the reduced frequency

*k*= 0.4, where

^{50}

^{,}Figure 1 shows the global mesh and local mesh of the NACA0012 airfoil flow field. To examine the independence of the results from the grid resolution, three meshes with about $ 2 \xd7 10 4$ (mesh 1), $ 4 \xd7 10 4$ (mesh 2), and $ 8 \xd7 10 4$ (mesh 3) mesh points are tested to choose the proper airfoil mesh for FNO training by RANS with the SA model. Table I shows the details of three meshes. The abbreviations PS, SS, and NY, respectively, correspond to the pressure side, the suction side, and the direction normal to the wall. Δn represents the normal spacing for the first cell away from the airfoil surface. Figure 2 shows the pressure coefficient (Cp) of the NACA0012 airfoil in AOA 9.25° and 12°. Cp is computed as

*ρ*is the density of the fluid, and

*p*and $ p \u221e$ are the local pressure and static pressure in the freestream. The DNS results are from Rodríguez

*et al.,*

^{51}and the results of meshes 1–3 are calculated from the RANS solution with the SA model. Comparisons between the RANS model results and DNS results show that for both AOA 9.25° and 12°, the RANS results with different mesh points on the windward side of the NACA0012 airfoil are in good agreement with DNS results; however, the RANS results with different mesh points on the leeward side are in bad agreement with DNS results. With the increase in mesh points, the RANS results in the leeward side cannot be better agreement with DNS results. The reason why RANS results are in bad agreement with DNS results is the lack of the SA model in simulating the flow with separation flow and adverse pressure gradient.

^{8}With the balance of computational accuracy and the training efficiency of the FNO model, mesh 2 with about $ 4 \xd7 10 4$ mesh points is chosen for dynamic stall simulation with the dimensionless wall normal spacing

*y*

^{+}of the first row of cells around the airfoil less than 1.

*y*

^{+}is the distance from the wall measured in viscous lengths denoted by Pope,

^{52}

^{,}

^{41}The training data are obtained from RANS with grid convergence, so the grid resolution is sufficient for FNO to predict flow fields of RANS.

Mesh . | Cell number on PS . | Cell number on SS . | Cell number in wake . | NY . | Total cells . | Δn/c . |
---|---|---|---|---|---|---|

Mesh 1 | 100 | 100 | 101 | 51 | 401 × 51 | 3.7 $ \xd7 10 \u2212 4$ |

Mesh 2 | 100 | 100 | 101 | 101 | 401 × 101 | 3.7 $ \xd7 10 \u2212 4$ |

Mesh 3 | 200 | 200 | 201 | 101 | 801 × 101 | 3.7 $ \xd7 10 \u2212 4$ |

Mesh . | Cell number on PS . | Cell number on SS . | Cell number in wake . | NY . | Total cells . | Δn/c . |
---|---|---|---|---|---|---|

Mesh 1 | 100 | 100 | 101 | 51 | 401 × 51 | 3.7 $ \xd7 10 \u2212 4$ |

Mesh 2 | 100 | 100 | 101 | 101 | 401 × 101 | 3.7 $ \xd7 10 \u2212 4$ |

Mesh 3 | 200 | 200 | 201 | 101 | 801 × 101 | 3.7 $ \xd7 10 \u2212 4$ |

All simulations are performed with the extensively validated Navier–Stokes solver OpenFOAM.^{53} The pimpleFoam algorithm, which is the combination of Pressure-Implicit with Splitting of Operators (PISO) algorithm and Semi-implicit method for pressure-linked equations (SIMPLE) in OpenFOAM, is applied for time marching.^{54} The second-order upwind scheme is applied for the spatial discretization of the convection term. The one equation Spalart–Allmaras (SA) turbulence model^{55} is used to simulate the viscous stresses of airflow around the airfoil. The detailed formulation of the SA model is shown in the Appendix. To resolve the motion of the mesh around the NACA0012 airfoil, a moving mesh has been used. The moving mesh in OpenFOAM is managed through the dynamicMotionSolverFvMesh solver. The oscillatingRotatingMotion function is used to capture the movement of the airfoil. The oscillatingRotatingMotion function creates a non-uniform rotational motion. The motion varies sinusoidally, creating an oscillating pattern. It is useful to model oscillating wheels, pendulums, or other rotating objects that follow a periodic motion.

The Reynolds averaged numerical simulation of dynamic stall of the NACA0012 airfoil is performed with the grid resolution of $ [ 410 \xd7 100 ]$ with the dynamic mesh technique. In Fig. 1(a), the left boundary is set as inlet flow, which is parallel with the x-axis, and the others are set as outlet flow for the case of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$. The left and down boundaries are set as inlet flow, which is in the 10° AOA to the airfoil, and the others are set as outlet flow for the case of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$. The locations of the boundaries are fixed at 10*c* from the airfoil. The no slip boundary condition is used on the surface of the airfoil. The normalized time step $ \Delta t$ is set as $ 1 \xd7 10 \u2212 5$ with the mean Courant number less than 1. The snapshots of the numerical solution are taken every 500 steps as a time node, and there are 200 time nodes in a dynamic stall period. Seven pitch oscillation periods have been simulated, and the first period is abandoned to eliminate the effect of the initial condition with 1200 time nodes (six periods) being saved totally. Each time node contains $ [ 410 \xd7 100 ]$ mesh points with the data of streamwise velocity, normal velocity, and pressure. The mesh points are rearranged according to Eq. (7) and transferred onto a uniform mesh. The flow field data with a tensor size of $ [ 1200 \xd7 410 \xd7 100 \xd7 3 ]$ can be obtained and serve as training and testing datasets with 800 samples (4 periods) used as the training dataset, 200 samples (1 period) used as the validation dataset, and 200 samples (1 period) used as the testing dataset.

## IV. RESULTS AND DISSCUSION

Figure 3 shows the schematic of FNO. In this deep learning model, the flow field data of the previous five time nodes $ [ U 1 , U 2 , U 3 , U 4 , U 5 ]$ are taken as input and the flow field data of the sixth time node $ [ U 6 ]$ are taken as output. Four Fourier layers with 12 modes for each layer are chosen. When the numbers of Fourier layers and modes are set to large values, it will lead to too many model parameters and the computational cost is too large. When the numbers of Fourier layers and modes are small, the model accuracy will be insufficient.^{42} The GELU function^{56} and adaptive moment estimation (Adam)^{57} are taken as the activation function and optimizer, respectively. FNO is trained for 500 epochs on an RTX 3080 GPU with a learning rate of 0.001^{42} and batch size of 10. We have tested different activation functions, optimizers, and learning rates. The performance is similar. The memory of GPU limits the selection of batch size. In our work, an RTX 3080 GPU is used for training the FNO model with 8 GB memory and the largest batch size that we can choose is 10. For a smaller batch size, the model can be converged faster, but the training time become longer. For a larger batch size, the training time become shorter, but the model needs more memory and epochs to achieve the same accuracy.

*l*

_{2}-loss is applied as the loss function for training FNO, as follows:

*u*,

*v*, and

*p*represent the streamwise velocity, normal velocity, and pressure predicted by FNO, and $ u \u0302 , \u2009 v \u0302 , \u2009 and \u2009 p \u0302$ represent the streamwise velocity, normal velocity, and pressure simulated by CFD.

### A. Prediction of streamwise velocity, normal velocity, and pressure

The streamwise velocity, normal velocity, and pressure of the previous five time nodes are set as the input of FNO, and the flow data of the sixth time node are set as the output of FNO.

Figure 4 shows the training loss and testing loss of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ and $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$. The red and blue curves represent the training and testing loss curves of FNO, respectively. As shown in Fig. 4, with the increase in epochs, the training and testing losses decrease rapidly for the first 200 epochs and gradually from 200 to 300 epochs. After 300 epochs, the training and testing losses of the low AOA case are close to the $ 10 \u2212 3$ level, and the training and testing losses of the high AOA case converge at the $ 10 \u2212 2$ level. In the high AOA case, the testing loss is larger than training loss, mainly due to the fact that the amount of training data is few for FNO to achieve more accurate prediction. As shown in Fig. 28–31, with the increase in the amount of data and parameter for FNO, more accurate prediction can be achieved with less training and testing losses. Both of training and testing losses of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ are much lower than those of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$. That is because the flow fields of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ with separation flow and adverse pressure gradient are much more complex than those of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$. Therefore, it is more difficult for FNO to predict the flow fields of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ than $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$.

*l*

_{1}-loss is defined as

*l*

_{1}-loss between the true data and the predicted data is in the range of the $ 10 \u2212 2$ level, which is low enough to predict the dynamic stall process of airfoils. Moreover, the prediction errors are concentrated on the wake flow of the airfoil, indicating that it is difficult for FNO to predict the vortex dynamics of wake flow.

*v*and

_{rms}*p*represent the RMSE of the normal velocity and pressure, respectively. The results indicate that when AOA increases, all RMSEs increase obviously.

_{rms}Figure 9 shows the temporal evolution curves of the pressure coefficient at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in the cycle of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ obtained by CFD simulation and FNO. The results indicate that the pressure coefficients predicted by FNO are of a great consistency with the results of CFD simulation.

Figure 10 shows the temporal evolution curves of the velocity magnitude near the airfoil surface at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in the cycle of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ obtained by CFD simulation and FNO. The data of velocity magnitude are obtained at the tenth mesh points from the airfoil surface. The results show that the velocity magnitudes predicted by FNO are of a great consistency with the results of CFD simulation.

Figures 11–13 show the contour plots of normalized streamwise velocity, normal velocity, and pressure of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in a cycle obtained by the CFD solver, FNO, and the absolute error. As the results show, FNO can give an accurate prediction for the dynamic stall with higher AOA, which contains separation flow and adverse pressure gradient. Compared with prediction results in Figs. 5–7, the error of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ is higher and concentrated on the area of the leeward side and wake of airfoil. The flow in the leeward side of the airfoil with the separation and adverse pressure gradient is much more complex than the windward side, which is almost the attached flow, as shown in Fig. 19. Therefore, the FNO model can achieve a highly accurate prediction for attached flow or flow without complex flow structure, but make less accurate prediction for separation flow or flow with complex flow structure.

Figure 14 shows the temporal evolution curves of root mean squared error (RMSE) of streamwise velocity, normal velocity, and pressure in the cycle of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ calculated by FNO. The results indicate that when the airfoil is in the downstroke, RMSEs reach the maximum and RMSE of pressure is much larger than others due to the larger absolute value.

Figure 15 shows the temporal evolution curves of the pressure coefficient at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in the cycle of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ obtained by CFD simulation and FNO. The results indicate that the pressure coefficients predicted by FNO are of a great consistency with the results of CFD simulation.

Figure 16 shows the temporal evolution curves of the velocity magnitude near the airfoil surface at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in the cycle of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ obtained by CFD simulation and FNO. The data of velocity magnitude are obtained at the tenth mesh point from the airfoil surface. The results show that the velocity magnitudes predicted by FNO are of a great consistency with the results of CFD simulation.

Figure 17 shows the lift coefficients $ ( C L )$ of the NACA0012 airfoil for two cases obtained by CFD simulation and FNO. The lift coefficients are the average of five dynamic stall periods. The results indicate that the lift coefficients predicted by FNO are of a great consistency with the results of CFD simulation.

Because FNO is a purely data-driven method without explicitly embedding any physical constraints, the flow fields obtained by FNO might not strictly satisfy the Navier–Stokers equations and the corresponding physical constraints. However, the FNO model is capable of learning the Navier–Stokers equations from flow data and can approximately satisfy conservation laws.^{49}

### B. Prediction of vorticity

The vorticity of the previous five time nodes is set as the input of FNO, and the vorticity of the sixth time node is set as the output of FNO. Figure 18 shows the contour plots of vorticity of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in a cycle obtained by the CFD solver, FNO, and the absolute error. The prediction results are also the same accurate as Figs. 5–7, which indicate that FNO can achieve an accurate prediction for vorticity. As shown in Fig. 18, the flow around the airfoil is almost attached without any separation or adverse pressure gradient. Moreover, the *l*_{1}-loss is still concentrated on the wake flow of the airfoil where exists a more complex flow structure, similar to that in Figs. 5–7.

Figure 19 shows the contour plots of vorticity of $ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ at different time nodes in a cycle obtained by the CFD solver, FNO, and the absolute error. Similar to the results in Figs. 11–13, the error is still concentrated on the leeward side and wake flow of the airfoil. As shown in Fig. 19, when AOA is low at 0T, laminar flow is attached to the airfoil except in the small trailing edge region, which was also seen in experiments.^{58} As AOA increases, a leading edge vortex (LEV) is forming at 0.1T, which causes a sudden rise in aerodynamic loads.^{58} The prediction error increases around the airfoil with the vortex structure becoming complex. With the growth of LEV, a vortex forming from the shear layer is shedding from the suction surface at 0.2T. At 0.3–0.4T, the first LEV is moving toward the trailing edge and leaving the suction surface, the second LEV is forming and moving toward the trailing edge. Based on the vorticity field, the dynamic stall occurs when the outer surface of the LEV meets the trailing edge.^{12} At the same time, a trailing edge vortex (TEV), which is counterclockwise rotating, is growing and sheds to the wake. More LEVs are growing and mixing with TEV at 0.5–0.7T. There is the largest prediction error in the mixing period at 0.5–0.7T. Finally, there are many small vortices shedding from the suction surface and reattached flow at 0.8–1.0T.

During the process of dynamic stall, the prediction vorticity error becomes larger when vortices are forming, growing, and shedding from the suction surface. This phenomenon is similar to the behavior of the streamwise velocity, normal velocity, and pressure errors. Table II shows comparison of accuracy and time consumption of vorticity between FNO, U-net, and UFNO. The results show that the UFNO model can give a more accurate prediction of flow fields than FNO and U-net with more training time consumption. Although the training time consumption of U-net is similar with FNO, the training loss and testing loss of U-net for both cases are the largest among these models. Figure 20 shows the contour plots of vorticity of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in a cycle obtained by the CFD solver, UFNO, and the absolute error. The results show that errors in the wake flow of the airfoil caused by UFNO at all time nodes are much less than FNO in Fig. 18. The prediction of vorticity can be improved by a more advanced model, such as UFNO. The results indicate that it is difficult for FNO to predict the forming, growing, shedding, and mixing of vortices in the separation flow and adverse pressure gradient. For the attached flow, FNO can achieve an accurate prediction in the cycle of the dynamic stall.

### C. Sequence flow fields prediction

The above-mentioned results demonstrate that FNO can predict the transient flow field with the inputs from the ground truth. In addition, the process of predicting sequence flow fields based on the iterative prediction strategy is achieved. In this process, the outputs of FNO are set as the inputs to perform transient prediction iteratively.

Figures 21–23 show the contour plots of normalized streamwise velocity, normal velocity, and pressure of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in a cycle obtained by the CFD solver, FNO, and the absolute error based on the iterative prediction strategy. In this case, the input and output of FNO are the streamwise velocity, normal velocity, and pressure data separately. The separated training of flow physical quantities can reduce the training error. The flow field data in the first 10 time nodes are set as the input with the flow field data in the next 20 time nodes as the output. FNO needs only ten iterations to achieve a whole dynamic stall prediction in this research. The results show FNO can achieve a whole dynamic stall process prediction iteratively in general. However, there is an error accumulation process in the iterations, which also occurs in Ref. 43. FNO can predict 200 flow fields (a whole dynamic stall period) within 15 s with a single GPU, but it costs about 30 h for a traditional CFD solver with four CPU cores. FNO is much more efficient than CFD simulation for predicting flow fields.

Figure 24 shows the temporal evolution curves of root mean squared error (RMSE) of streamwise velocity, normal velocity, and pressure in the cycle of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ calculated by FNO based on the iterative prediction strategy in three cycles of flow. The results indicate that with the increase in time step, RMSEs of streamwise velocity and pressure increase gradually, which show the error accumulation process, but RMSE of normal velocity maintains a lower level. The reason may be that the variation of normal velocity flow field is relatively simple and the absolute value of normal velocity is much lower than others.

Figure 25 shows the temporal evolution curves of the wall pressure coefficient at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0 T) in the cycle of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ obtained by CFD simulation and FNO based on the iterative prediction strategy. The results indicate that the pressure coefficients predicted by FNO are in good agreement with the results of CFD simulation.

Figure 26 shows the temporal evolution curves of the velocity magnitude near the airfoil surface at different time nodes (0.2T, 0.4T, 0.6T, 0.8T, and 1.0T) in the cycle of $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ obtained by CFD simulation and FNO based on the iterative prediction strategy. The data of velocity magnitude are obtained at the tenth mesh point from the airfoil surface. The results show that the velocity magnitudes predicted by FNO are of a great consistency with the results of CFD simulation.

Figure 27 shows the lift coefficients of the NACA0012 airfoil for $ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ obtained by CFD simulation and FNO based on the iterative prediction strategy. Because the predicted results obtained by FNO are only the lift coefficients of one dynamic stall period, there is a data oscillation between the CFD and FNO results. The results show that the lift coefficient predicted by FNO is in a good agreement with the results of CFD simulation.

### D. Parameters analysis

Figure 28 shows the comparison of accuracy and training time consumption with different cases using different numbers of Fourier layers. With the increase in the number of Fourier layers, the training loss and testing loss become lower, and the time for training becomes longer. When the number of Fourier layers increases from 2 to 4, both losses decrease gradually; when it increases from 4 to 5, the losses decrease a little. The results show that when the Fourier layer is 4, the training loss and testing loss seem to be very close to the minimum. Five Fourier layers improve a little accuracy than four Fourier layers, but more time is consumed. Thus, the four Fourier layers framework is suitable for the prediction of dynamic stall of the airfoil.

Figure 29 shows the comparison of accuracy and training time consumption using different training periods. The results show that as the training periods increase, the training loss and testing loss decrease obviously, but the training time also increases linearly. A more accurate predicted flow field can be obtained by increasing the training periods.

Figure 30 shows the comparison of accuracy and training time consumption using different sampling intervals. $ \Delta N$ represents the normal sampling interval. With the increase in sampling intervals, the training and testing losses increase obviously, and the training time decreases. Increasing the sampling interval can improve the training efficiency but get a prediction result with a larger error.

Figure 31 shows the comparison of accuracy and training time consumption using different modes. The results show that as the number of modes increases, the training loss and testing loss decrease obviously, but the training time also increases linearly. For 12 and 20 modes, the FNO model can achieve a high prediction accuracy. However, the FNO model with 20 modes needs more training time, so the FNO model with 12 modes is chosen for the flow field prediction.

Table III shows the comparison of accuracy and training time consumption between FNO, U-net, and UFNO. The results show that the UFNO model can give a more accurate prediction of flow fields than FNO and U-net with more training time consumption. Although the training time consumption of U-net is similar to FNO, the training loss and testing loss of U-net for both cases are the largest among these models.

. | Model . | Train loss . | Test loss . | Time for training . |
---|---|---|---|---|

$ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.001 78 | 0.001 80 | 3326 s |

U-net | 0.042 10 | 0.022 10 | 3546 s | |

UFNO | 0.001 14 | 0.000 90 | 4843 s | |

$ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.011 03 | 0.025 71 | 3260 s |

U-net | 0.069 40 | 0.041 50 | 3668 s | |

UFNO | 0.007 80 | 0.014 64 | 4756 s |

. | Model . | Train loss . | Test loss . | Time for training . |
---|---|---|---|---|

$ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.001 78 | 0.001 80 | 3326 s |

U-net | 0.042 10 | 0.022 10 | 3546 s | |

UFNO | 0.001 14 | 0.000 90 | 4843 s | |

$ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.011 03 | 0.025 71 | 3260 s |

U-net | 0.069 40 | 0.041 50 | 3668 s | |

UFNO | 0.007 80 | 0.014 64 | 4756 s |

## V. CONCLUSION

In this study, Fourier neural operator is applied to predict the flow field of the NACA0012 airfoil in the dynamic stall process. Three meshes are used to choose the proper mesh with the balance between computational accuracy and neural network training efficiency. Two cases with different angles of attack are simulated by RANS with the SA model at $ R e = 4 \xd7 10 4$. A prediction model is directly constructed between the flow fields at the previous time nodes and that at the future time node by FNO. Four Fourier layer neural network is chosen to predict the flow field with better accuracy and less training time. FNO can achieve a fast and accurate prediction of streamwise velocity, normal velocity, and pressure for both cases. Compared with the absolute error of the low AOA case concentrating on the wake flow of the airfoil, there is a higher error in the high AOA case concentrating on the leeward side and wake flow. It is more difficult for FNO to predict the forming, growing, shedding, and mixing of vortices in the separation flow and adverse pressure gradient. For the attached flow, FNO can achieve an accurate prediction in the cycle of the dynamic stall. The prediction of sequence flow fields based on the iterative prediction strategy is achieved for dynamic stall. FNO can predict the whole dynamic stall process in general, but there is an error accumulation process in the iterative prediction. FNOs with different configurations are tested to demonstrate the effects of Fourier layers, training periods, sampling intervals, and modes on prediction accuracy and time consumption.

## ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (NSFC Grant Nos. 91752202 and 12172161).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Deying Meng:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Yiding Zhu:** Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Writing – review & editing (equal). **Jianchun Wang:** Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). **Yipeng Shi:** Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

### APPENDIX: SPALART–ALLMARAS MODEL

*ν*as follows:

_{t}^{55}

*P*and

*D*are the production and destruction terms of $ \nu \u0303$ given as

*f*is defined as follows:

_{w}The model constants are $ c b 1 = 0.1355 , \u2009 \sigma = 2 / 3 , \u2009 c b 2 = 0.622 , \kappa = 0.41 , \u2009 c w 1 = c b 1 / \kappa 2 + ( 1 + c b 2 ) / \sigma , \u2009 c w 2 = 0.622 , \u2009 c w 3 = 2.0$, and $ c \upsilon 1 = 7.1$. Freestream boundary is set to fully turbulent with $ \nu \u0303 / \nu \u221e = 4$, and $ \nu \u0303$ is set to zero at no-slip walls.

. | Model . | Train loss . | Test loss . | Time for training . |
---|---|---|---|---|

$ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.002 60 | 0.003 06 | 3064 s |

U-net | 0.082 53 | 0.087 46 | 3372 s | |

UFNO | 0.002 26 | 0.002 02 | 4581 s | |

$ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.016 33 | 0.041 33 | 3032 s |

U-net | 0.176 03 | 0.334 78 | 3486 s | |

UFNO | 0.014 72 | 0.027 27 | 4519 s |

. | Model . | Train loss . | Test loss . | Time for training . |
---|---|---|---|---|

$ \alpha = 0 \xb0 + 5 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.002 60 | 0.003 06 | 3064 s |

U-net | 0.082 53 | 0.087 46 | 3372 s | |

UFNO | 0.002 26 | 0.002 02 | 4581 s | |

$ \alpha = 10 \xb0 + 10 \xb0 \u2009 \u2009 sin ( 2 \pi f t )$ | FNO | 0.016 33 | 0.041 33 | 3032 s |

U-net | 0.176 03 | 0.334 78 | 3486 s | |

UFNO | 0.014 72 | 0.027 27 | 4519 s |

## REFERENCES

^{4}and 6 × 10

^{4}

*Ansys Fluent Theory Guide.*

*Turbulent Flows*

*Mathematics, Numerics, Derivations and OpenFOAM®*