Evaluation of Actuator Disk Model Relative to Actuator Surface Model for Predicting Utility-Scale Wind Turbine Wakes

: The Actuator Disk (AD) model is widely used in Large-Eddy Simulations (LES) to simulate wind turbine wakes because of its computing efﬁciency. The capability of the AD model in predicting time-average quantities of wind tunnel-scale turbines has been assessed extensively in the literature. However, its capability in predicting wakes of utility-scale wind turbines especially for the coherent ﬂow structures is not clear yet. In this work, we take the time-averaged statistics and Dynamic Mode Decomposition (DMD) modes computed from a well-validated Actuator Surface (AS) model as references to evaluate the capability of the AD model in predicting the wake of a 2.5 MW utility-scale wind turbine for uniform inﬂow and fully developed turbulent inﬂow conditions. For the uniform inﬂow cases, the predictions from the AD model are signiﬁcantly different from those from the AS model for the time-averaged velocity, and the turbulence kinetic energy until nine rotor diameters ( D ) downstream of the turbine. For the turbulent inﬂow cases, on the other hand, the differences in the time-averaged quantities predicted by the AS and AD models are not signiﬁcant especially at far wake locations. As for DMD modes, signiﬁcant differences are observed in terms of dominant frequencies and DMD patterns for both inﬂows. Moreover, the effects of incoming large eddies, bluff body shear layer instability, and hub vortexes on the coherent ﬂow structures are discussed in this paper.


Introduction
Nowadays, large wind farms are constructed to respond to the increasing demand of renewable energy. In these wind farms, turbines are installed in cluster to meet the geographical restriction and to reduce the cable and maintenance cost. A turbine may influence its downwind neighbors significantly with the wake effect, leading to a loss of the power production and an increase of the unsteady load on the structure [1]. Therefore, a need to better understand the wake behavior and its influence on the downwind turbines arises.
Understanding turbine wakes in a wind farm is challenging because of its multi-scale nature. For example, the boundary layer on a wind turbine's blade has a thickness of the centimeters, which is orders of magnitudes smaller than the diameter of the rotor (≈100 m) and the thickness of the Atmospheric Boundary Layer (ABL; ≈1000 m) [2]. Among others, the difficulty in accurately modeling the flow around the blade of a real wind turbine blade arises both in wind tunnel experiments (due to scale effect [3]) and in numerical simulations (due to the resolution requirement).
Facing this challenge, it is often assumed that the actual geometry of the rotor is of less importance [4] if only the far wake and its influence on downwind turbines are of interest, so that the wind turbine can be approximated by equivalent models. In experiments, the simplest model is the porous disk model, which is broadly used [5,6]. Validations in wind tunnels have demonstrated that porous discs can provide time-averaged wake properties with satisfactory accuracy in the region further than 3.5 rotor diameters downstream, especially when the turbulence of the ABL is concerned [7][8][9].
In numerical simulations, a series of blade approximated actuator type models representing the turbine blades using equivalent distributed forces have been proposed. The way in which these forces are calculated and distributed distinguishes different models. The simplest is the Actuator Disk (AD) model, which is the numerical equivalence of a porous disk. The thrust on the disk is calculated with one-dimensional momentum theory and is usually distributed uniformly over the rotor swept area with the rotation effects neglected. The Actuator Line (AL) method was proposed to take into account the effects of individual rotating blades [10]. The AL models a wind turbine blade by a rotating line with lift and drag forces determined from tabulated geometric and aerodynamic data of airfoil. To better take into account the geometrical effects of wind turbine blades, the Actuator Surface (AS) method has been proposed, which models a blade as a two dimensional surface with zero thickness [11,12]. Because of its simplicity and computational efficiency, the AD has been widely used in turbine wake simulations especially in farm-scale simulations [13][14][15][16]. The capability of the actuator disk model in predicting turbine wakes, especially in the far wake region, has been widely validated in the literature [14,17,18]. However, besides the thrust, experiments revealed that the rotor's rotation also influences both the power output and the wake characteristics significantly [19] and including these rotational effects in the actuator disk model can improve the model's accuracy [20]. In [21], the authors showed that the actuator disk model can reasonably predict the mean velocity profiles but underpredict the turbulence kinetic energy (TKE) for the wake of an axial-flow hydrokinetic turbine. It is noticed that most of the validation studies were focused on time-averaged quantities without probing into the dynamic behavior of turbine wakes, e.g., coherent flow structures and the wake meandering, for which the dataset is difficult to obtain from utility-scale wind turbines [22]. Furthermore, inconsistent results were observed in wind tunnel experiments on the dynamic behavior of turbine wakes when different turbine models were used. For instance, regarding the origin of wake meandering, Medici et al. [3,23] found the wake meandering was related to the bluff body vortex shedding in the experiment of a small scale wind turbine, whereas Espana et al. [6] claimed that the meandering was attributed to the inflow large eddies by carrying out an experiment by representing the turbine with a porous disk. These wind tunnel measurements already make it questionable whether the AD (or the porous disk) model can predict correctly the dynamics of small scale wind turbine wakes. Less is known when applying such a model to utility-scale wind turbines.
To this end, the present study employs simulation results from the well-validated AS model proposed in [11] to examine the capability of the AD model in predicting the dynamic behavior of a utility-scale wind turbine under uniform and fully developed turbulent inflow conditions. Large-Eddy Simulation (LES) is employed for turbulent flow simulations. For both models, exactly the same computational setup is employed. We first compare the time-averaged quantities and then employ the dynamic mode decomposition (DMD) to facilitate the comparison of the most dominant dynamic flow structures and the frequency spectra between the AD and AS models.
The remainder of this paper is structured as follows. Section 2 presents the theory of the AD and the AS models together with a brief description of the LES solver and the DMD method. In Section 3 the simulation setup is provided. Section 4 depicts the simulation results and the DMD analysis in both uniform and ABL conditions. A discussion is provided in Section 5 before the final conclusion in Section 6.

Numerical Methods
This section describes the concept of the AD and the AS wind turbine model, the LES flow solver employed in this study, and the dynamic modal decomposition (DMD) employed to analyze the dynamic flow structure from the simulations.

Actuator Disk Model
The AD model neglects the geometry detail of individual wind turbine blades. It represents the rotating blades as a fixed 2D porous disk exerting a uniform thrust on the flow, which is the numerical reflection of the perforated disk model usually used in wind tunnel experiments. Neither the rotational effect nor the non-uniform force distribution are considered in the model employed in this work.
The axial thrust force f T per unit area is uniformly distributed over the entire rotor disk surface A and is expressed with the thrust coefficient C T and the inflow velocity V ref : where ρ is the density of air. The reference velocity V ref is defined to be equal to the freestream velocity in uniform inflow condition. In turbulent inflow simulations, the present work approximately calculates V ref by averaging the velocity on a disk of the rotor's size at one diameter in front of the real turbine. The trust coefficient C T remains to be determined according to the turbine operation state.
In this work, C T is set to be equal to that of the AS simulations to ensure a fair comparison between the two models.

Actuator Surface Model
The AS model represents the geometry of an individual wind turbine blade with a simplified two dimensional surfaces of zero thickness, which is formed by chords at different radial locations [11,12]. In the actuator model employed in this work, the aerodynamic forces on the surface vary with the radial position and are determined by the tabulated airfoil data in the same ways as the AL model as follows: and where L and D are the lift and drag force per unit length, ρ is the density of air, c is the chord length, V ref is the flow velocity relative to the rotating blade, e L and e L are unit directional vectors for lift and drag forces. C L and C D are the lift and the drag coefficients defined in 2D airfoil tables as a function of Reynolds number and the angle of attack. Corrections including the 3D stall delay model of Du and Selig [24] and the tip loss correction of Shen et al. [25,26] are applied. After calculating L and D, the force f per unit area on the surface model is calculated by: The reacting forces exerting by the blade on the air are then distributed to the background Eulerian grid points with a smoothed discrete delta function [27].

Flow Solver
The turbulent flow is solved with a three dimensional LES code, dubbed as virtual flow simulator (VFS-Wind) [28], in which the governing equations are the filtered Navier-Stokes equations for incompressible flows, which read in the compact tensor form as (i, j, k, l = 1, 2, 3): where ξ i is the curvilinear coordinates, ξ i l = ∂ξ i /∂x l is the transformation metrics with x l the Cartesian coordinates, J denotes the Jacobian of the geometric transformation, U i = ξ i l /J u l is the contravariant volume flux with u l the velocity in Cartesian coordinates, µ is the dynamic viscosity, ρ is the density, g jk = ξ j l ξ k l is the components of the contravariant metric tensor, p is the pressure, f l is the body forces exerted by the actuator models, and τ ij is the sub-grid stress (SGS) resulted from the filtering operation and is modeled with the Smagorinsky SGS model [29] as follows, where µ t is the eddy viscosity and S ij is the large-scale strain-rate tensor with (·) denoting the grid filtering operator. The eddy viscosity is computed by where ∆ is the filter width, |S| = (2S ij S ij ) 1/2 is the magnitude of the strain-rate tensor with C s the Smagorinsky constant computed via the dynamic procedure of [30]. A second-order accurate central differencing scheme is used for space discretization. The time integration uses the fractional step method [31]. The momentum equation is solved with a matrix-free Newton-Krylov method [32] . The pressure Poisson equation is solved with a Generalized Minimal Residual (GMRES) method with an algebraic multi-grid acceleration [33].

Dynamic Mode Decomposition
DMD is an equation-free, data-driven method for data analysis and behavior prediction of complex dynamical systems. It was first proposed by Schmid [34] to analyze the high-dimensional fluid dynamics data by decomposing it into coherent spatial structures that oscillate at distinct frequencies.
Thanks to its ability both to analyze and to predict, it has gained successes not only in fluid mechanics, but also in other fields, including video processing and finance, where high-dimensional complex dynamic systems are involved [35].
For fluid mechanics applications, the input of DMD analysis is a sequence of snapshots of the flow field. The snapshot x i is a column vector of dimension n containing all the interested variables at measure points in the flow field at time t = t i . The snapshots are taken at a fixed time interval ∆t. With total m snapshots, the snapshot matrix of the dataset can be written as: where the sub and sup indexes indicate the starting and the ending index of snapshots. The DMD approximates the dynamic system by a time-independent linear matrix A representing the temporal evolution from one snapshot to the next, as follows, In practice, the spatial dimension n is usually very large, which makes the direct solution of A n×n difficult. Instead, the reduced-order Proper Orthogonal Decomposition (POD) projection matrix A r×r with r n is used. r is the number of the most energetic modes kept in the singular value decomposition of X m−1 1 , defined as follows, where the left singular vectors in U are POD modes, the right singular vectors in V are time coefficient of these modes, and Σ is a diagonal matrix containing the first r largest singular values. U and V are orthonormal.Ã is defined asÃ The next step is to calculate the eigenvalues and eigenvectors ofÃ, where Λ is a diagonal matrix containing r eigenvalues and W contains r eigenvectors of unit length.
The snapshots x k = Ax k−1 = A k−1 x 1 can be written as where Φ = UW contains the DMD modes (φ i ) and b = W −1 U T x 1 contains the amplitudes (b i ) of these modes in the first snapshot. Each mode φ i has a corresponding eigenvalue λ i in Λ. The oscillating frequency ( f i ) of mode φ i is equal to and the growth rate g i is equal to For a stable oscillating system, all the growing rate should be equal to 0. A positive growing rate (g i > 0) indicates that the mode φ i has an amplitude b i increasing with time; a negative growing rate (g i < 0), in contrast, indicates a mode damping with time.
With the amplitudes b i , the most energetic DMD modes in the first snapshot can be easily identified and remain the same for other snapshots when the system is stable. In contrast, when growing and damping modes are concerned, the modes should be ordered by the time-averaged amplitudes b i . To calculate the time-averaged amplitudes, the eigenvalue-weighted method proposed by Kou and Zhang [36] is selected for its simplicity and computational efficiency among others [37]. The eigenvalue-weighted amplitudes b i are calculated as,

Simulation Setup
The AD is compared against the AS models by simulating a three-blade Clipper Liberty 2.5 MW wind turbine located at the EOLOS wind energy research field station in University of Minnesota, USA.
The rotor diameter is D = 96 m, the hub height is z hub = 80 m, and the nacelle has a near cuboidal shape with dimensions of 5.3 m × 4.7 m × 5.5 m. The tower as a cylindrical form with a diameter of 3.0 m at the top and 4.1 m at the bottom, respectively. The readers can find more information about this wind turbine in previous works [22,38,39]. Because of proprietary issues, the details of the blade geometry cannot be released in this paper. Interested readers can contact the EOLOS wind energy consortium (Email: eolos@umn.edu, Address: St. Anthony Falls Laboratory, 2 Third Avenue SE, Minneapolis, MN 55414, USA) at the University of Minnesota for these details.
The capability of the employed AS model has been evaluated for different aspects previously. In [40], the employed method was validated against wind tunnel measurements for the time-averaged flow quantities in the wake, such as the velocity deficit and turbulent intensity. Moreover, validations using the same Clipper wind turbine have shown that the AS model is able to predict accurately the near-wake vortex structures as compared with the field measurement using snow-based super-large-scale particle image velocimetry (SLPIV) [41]. Due to these previous validations, in this work, we consider the AS simulation results as references for evaluating the AD model.
Here we present the computational setup for the simulations carried out in this work. In both AD and AS cases, the size of the computational domain is set as L x × L y × L z = 14D × 7D × 10D, where x, y, z represent the stream-wise, the span-wise, and the vertical directions, respectively. The domain is discretized with a Cartesian grid of N x × N y × N z = 281 × 281 × 143. The grid size is uniform in the x, y directions with ∆x = D/20 and ∆y = D/40. In z direction, the mesh is uniform with ∆z = D/40 in z ∈ (0, 2D) region to resolve the wind turbine wake and the interaction with the ground, and is gradually stretched to the top of the computational domain. Figure 1 shows the disk and the surface discretized with unstructured triangular surface mesh. Please note that a nacelle model [11] is employed in both AD and AS cases, otherwise there will be a non-realistic jet flow behind the empty rotor center for the AS method. Furthermore, the vortex shedding from the nacelle plays an important role in the wake evolution because of its interaction with the root vortex and the tip shear layer [21]. Although it would be ideal to take the tower into account to have a complete representation of a realistic wind turbine. However, including the tower (diameter = 3.0 m at the top) gives rise to numerical difficulties since the present study (and most numerical studies on the wind turbine's wake as well) employs a grid which is too coarse to resolve the flow details around the tower, and thus complicates the comparison between the AD and AS models. Furthermore, it was shown in [42] that the effect of tower is limited to the near wake region. For these reasons, the tower was not considered and we focus on the differences caused by the two rotor models.
In the AS simulations, the turbine rotates at a fixed tip speed ratio (TSR = ΩR/U = 8, where Ω is the rotor rotational speed, R is the rotor radius and U is instantaneous streamwise velocity averaged over a disk of radius R located 1D upstream of the turbine). The thrust is recorded at each time step and then averaged to calculate the thrust coefficient C T for the AD model. In the AD simulations, the thrust coefficient, which is computed from the corresponding AS simulations, is employed to compute the thrust on disk using Equation (1).
The simulations are conducted with two inflow conditions, i.e., a uniform and a fully developed turbulent inflow. In both cases, the streamwise velocity at the rotor's hub height is U hub = 9 m/s. The Reynolds number based on D and U hub is equal to Re = DU hub /ν = 5.7 × 10 7 . For the turbulent inflow case, the turbulence density is σ u /U hub = 0.08 at the hub height. The flow at inlet boundary is computed from a precursory LES with a larger computational domain of L x × L y × L z = 62D × 46D × 10D to capture large scale eddies in the incoming flow. In this inflow generation approach, the velocity fields on a y-z plane are first saved for each time step in the precursory simulation and then applied at the inlet of the turbine simulation. If the mesh and the size of time step employed in the precursory simulation are different from those in the wind turbine simulations, linear interpolations in both space and time are carried out to obtain the inflow velocity for the turbine simulations. Periodical boundary condition is applied in the horizontal directions. The upper boundary condition is the free slip. The wall model based on the logarithmic law is applied on the ground (the roughness length is z 0 = 5 × 10 −3 m for the present cases). For the uniform inflow cases, the boundary condition on the lateral walls is the free slip. An extra inflow only simulation with an empty computational domain is carried out with the same setups of the turbulent inflow case to help identify the contribution of turbulent large eddies.

Results
In this section, the simulation results are presented. The uniform inflow cases are presented in the first place and are followed by the turbulent inflow simulations. For each inflow condition, we compare instantaneous flow fields, time-averaged flow fields and DMD modes from the AD simulations with those from the AS simulations. In this section, u, v, w denote the instantaneous flow velocity in the streamwise, spanwise and vertical direction, respectively, with U, V, W for the time-averaged values. Figure 2 depicts the simulated instantaneous velocity fields behind the AD and the AS wind turbine models on the z = z hub plane. For the streamwise velocity contour in Figure 2a,b, it is found that both wake boundaries are first stable in a small distance behind the turbine and then show fluctuations in the far wake. Inside the wake away from the nacelle, the velocity deficit behind the AD model is more evenly distributed along the radial direction than the AS model, since the tip and the root losses and the radial variation of blade sections are not considered in the AD model. In the hub region, the wake of the nacelle is observed in both models. However, a jet which encompasses the nacelle's wake appears uniquely behind the AS model. Figure 2c,d show the spanwise velocity field. Fluctuations appear behind both wind turbine models. In the near wake region, a Kármán vortex street pattern is observed in the centerline of the near wake due to the nacelle. A significant discrepancy emerges in 2 < x/D < 3 region, where the AD model's result shows a regularly oscillating pattern on the wake boundary as observed for both streamwise and spanwise velocities (Figure 2a,c). It is also noticed that the spanwise velocity fluctuations in this region are stronger near the wake boundary than those in the wake center. In contrast, the wake behind the AS model does not have this oscillatory boundary, and the vortex shed from the nacelle grows gradually and starts to influence the wake boundary at x ≈ 4D. In the far wake (x > 6D), the spanwise velocity behind the AS model seems to be more energetic than that of the AD model. Quantitative comparisons shall be conducted in the next section to confirm this observation.

Time-Averaged Flow Field
In Figure 3 we compare the time-averaged flow field computed from the AD and the AS simulations. As seen in Figure 3a, the mean velocity profiles of both models show an overall agreement in the far wake (x = 9D). However, immediately behind the wind turbine, the two velocity profiles differ significantly. The profile of AD is almost uniform (except for the region near the nacelle). In contrast, the profile of AS shows a clear radial variation, which is remarked by a weaker deficit behind the nacelle due to the root loss (y = 0) and a smoother transition on the wake boundary due to the tip loss (y = ±0.5D) at x = 1D. In the near wake, the thickness of the shear layer on the wake boundary is smaller for AD, but it grows faster than that in the AS simulation. At x = 5D, it is obvious that this transitional region is thicker for the AD than the AS. This faster growing of wake boundary thickness denotes a quicker recovery and expansion of the wake of the AD. By comparison, this transitional region has no remarkable development in the result of AS until x = 7D and expands faster from 7D to 9D. In Figure 3b,c, the turbulence kinetic energy (TKE; k) and the primary Reynolds stress (<u v >) both indicate that the AD model shows stronger turbulent effects on the wake boundary in the x < 5D region and is surpassed by the AS in the far wake (x > 7D). At x = 9D the velocity profile of the two models are in reasonable agreement, while the wake computed by the AS model contains more turbulence kinetic energy and larger Reynolds stress. This result confirms the observation from the instantaneous flow field in Figure 2.

DMD Analysis
DMD is conducted to analyze the dynamic coherent structures in the wakes. The velocity field on the horizontal plane at the hub height (z = z hub ) is analyzed. Figure 4 depicts the DMD amplitude spectra and the most dominants modes.
In Figure 4a,e, only the modes with Strouhal number less than three (St = f D/U hub < 3) are plotted although the entire spectra are much longer, because the most energetic modes are within this low frequency range. First, these amplitude spectra reveal a significant difference in the energy distribution over the frequencies of the AD and the AS models' wake. The wake from the AS simulations contains generally more energy in the low frequency range below St < 1.5 and shows a trend of concentration around St ≈ 0.7 (marked with φ 1 in Figure 4e). However, such a concentration trend does not appear in the spectrum in the AD's case, whose DMD modes are of almost the same amplitudes in 0 < St < 2 region, except for a distinct peak at St ≈ 1.8 (marked with φ 1 in Figure 4a). Secondly, when comparing the amplitudes between the two cases, it is found that modes of the AS have slightly larger amplitudes than that of the AD, especially in the low frequency region, showing the wake of AS contains more energetic low frequency oscillations. Thirdly, no energy concentration around the vortex shedding frequency of the bluff body (St ≈ 0.168) is observed in both spectra. Figure 4b-d,f-h show the three most energetic modes of the AD and the AS cases, respectively. Overall, the spatial scale of the oscillation patterns enlarges as the Strouhal number decreases. However, the results from the AD and the AS models have very different dominant modes. Figure 4b shows the mode of the largest amplitude of the AD case. A spatial energy concentration on the wake boundary around 2D < x < 4D is found, which is in agreement with the instantaneous flow field (Figure 2c). Interestingly, no apparent source of disturbance can be traced in the upstream. It suggests that this mode should perhaps be related to the instability of the thin shear layer on the wake boundary that amplifies tiny disturbances in the flow field. It is noticed that this mode dominates for 2D < x < 4D, but becomes weak in the far field. In contrast, the other two modes of the AD (Figure 4c,d) both stem from the nacelle and dominate the far wake. As to DMD modes in the AS case (Figure 4f-h) all the three modes stem from the upstream nacelle and develop until the far wake. They differ from each other by the frequency and the wave length as shown by the spanwise velocity field. No AS mode similar to φ 1 of the AD is found after checking all the modes of the AS (including those not shown in the figure). This observation suggests that for the uniform inflow condition, the wake computed from the AS model is strongly affected by the vortex shedding behind the nacelle, while the AD model predicts a unique mode related to the shear layer instability near the wake boundary, which may not exist in real wind turbine wakes.

Turbulent Inflow
This section presents the simulation results from the AD and AS cases under the fully developed turbulent inflow condition. Figure 5 depicts instantaneous velocity fields computed from (a) the AD and (b) the AS models and (c) the inflow case on the z = z hub plane at the same instant. As seen, the wake shapes and the spanwise velocity contours from the AD and AS cases show a reasonably good agreement, except for some differences, such as the jet flow behind the hub of the AS model. Compared with the uniform inflow condition, strong effects of inflow turbulence is observed. For instance, a much stronger spanwise velocity is observed for the turbulent inflow (Figure 5d,e) when compared with that from the uniform inflow (Figure 2c,d). When comparing with the only inflow case (Figure 5f), we can observe an interesting phenomenon, that the patterns of the spanwise velocity computed from the AD and AS model, are very similar to that from the only inflow case. Moreover, it is noticed that both AD and AS models amplify the spanwise velocity fluctuations as the velocity magnitude is larger in the wake than in the inflow. Apart from this, no more information can be extracted with confidence from this instantaneous flow field yet.  Figure 6 shows the profiles of time-averaged flow quantities. Before the following detailed analysis, it is clear that the results from the AD model are in an overall good agreement in the far wake (x > 7D) for turbulent inflow condition, in contrast to Figure 3 of the uniform inflow condition. For the streamwise velocity, differences between the AD and the AS models are only observed in the near wake region (x < 3D) and become insignificant at far wake locations as shown in Figure 6a. The wakes computed by both models slightly skew towards the −y direction due to the insufficient simulation time, which causes the inflow not perfectly symmetric with respect to y = 0. The spanwise profiles of TKE are shown in Figure 6b. As seen, the TKE is concentrated near the hub (y = 0) and the wake boundary (y ± 0.5D) at x = 1D with larger TKE near the hub and near the boundary for the AS and AD models, respectively. The region of high TKE increases and expands as the wake travels from (3D < x < 5D) and the TKE computed by the AS model develops faster and surpasses that computed by the AD after (x > 3D). At further turbine downwind locations, the two models show consistent results. Similar trends are observed for the Reynold's stress <u v > as shown in Figure 6c. In Figure 6d,e, the vertical distribution of streamwise velocity and the turbulence intensity also confirms the above conclusion that the two models agree well in the far wake and the differences only manifest in the near wake region.

Time-Averaged Flow Field
To this extend, the AD model can reasonably predict the time-averaged flow quantities in the far wake for turbulent inflow conditions.  Figure 7 depicts the DMD amplitude spectra and the most dominants modes. Figure 7a shows the spectrum computed by the AD model with three distinct peaks remarked. As seen, the Strouhal number of the first peak (φ 1 ) is 0.17, which falls in the range of the bluff body vortex shedding frequencies. This dominant mode, illustrated by Figure 7b, shows that the transversal velocity in the wake alters the direction regularly in a way resembling to the wake behind a circular cylinder [43]. The other two modes of the AS illustrated in Figure 7c,d are at higher frequency. The spatial scale of the oscillation decreases as the Strouhal number increases. These two modes have larger amplitudes on the wake boundary than in the wake center. For the second mode φ 2 , stronger velocity amplitude is observed in the far wake, whereas φ 3 shows a local concentration of energy in 2D < x < 6D.

DMD Analysis
On the other hand, the spectrum of the AS case ( Figure 7e) shows a dominant peak close to the bluff body vortex shedding frequency at St = 0.23, which is the second larger peak in the DMD spectrum with the largest mode at St = 0.08. The first mode of the AS case (wavelength ≈ 9D) is significantly larger as shown in Figure 7f. The second mode from the AS case (Figure 7g), on the other hand, is close to the first mode from the AD case, although the wavelength of φ 2 from the AS case is smaller. Moreover, it is observed that the spectrum from the AS case has larger energy at the low frequency range than that of the AD. The two modes of higher frequency (φ 2 and φ 3 ) computed from the AS model are both stronger in the far wake than in the near wake with both larger spatial scales when compared with the AD case.
It is also worth noting that the influence of the nacelle on the far wake is less evident compared to the uniform inflow case, that two (the AS case) or three (the AD case) dominant DMD modes are obviously related to the nacelle for the uniform inflow cases, while only the third dominant mode from the AS case ( Figure 7h) seems to be magnified by the nacelle for the turbulent inflow cases. Detailed space-time correlation study [44] is needed to further examine the effect of nacelle on wake dynamics under turbulent inflows.

Discussion
The previous section has shown the differences between the wakes computed from the AD and AS models for both uniform and turbulent inflow conditions. For the uniform inflow cases, external perturbations are excluded and the wakes computed by both models develop on their own properties. The differences are easily identified already on the instantaneous and the time-averaged flow field: the instability near the wake boundary develops earlier for the AD model and results in a faster growth of turbulence near the wake boundary and quicker wake expansion and recovery in the near to intermediate turbine downwind locations. This phenomenon can be explained with the classical linear stability theory of shear flow [45]. According to this theory, the optimal spatial scale of perturbation to destabilize a shear flow is proportional to the transitional thickness so the sharp transition between the freestream flow and the wake of the AD case (as the tip loss effect is neglected) is more sensitive to small perturbations. This sharp transition across the wake boundary and the resultant stronger mixing and expansion in the near wake of the AD model are in accordance with previous wind tunnel experiments of Lignarolo et al. [8,46]. In the DMD analysis, we have compared the two models with the energy spectra and the dominant modes. For both models, a clear influence of the nacelle on the far wake is shown by the DMD modes, which is in agreement with previous numerical studies [21]. However, obvious disparities are also found that the spectrum of the AS case is concentrated in a lower frequency range, whereas that from the AD case has a distinct high frequency dominant mode corresponding to the shear layer instability. Moreover, the wake computed from the AS model has larger coherent structures and oscillates at a lower frequency than that from the AD model, which may also due to the differences of the shear layer near the wake boundary.
For turbulent inflow cases, the velocity deficit recovers faster than the uniform case due to the inflow turbulence, which is in agreement with previous studies [16,47]. In our test cases, the time-averaged streamwise velocity, TKE, and primary Reynold's stresses computed by the AD model reasonably agree with the AS model starting from x = 7D. This generally good agreement of the AD and the AS in turbulent inflow condition is in accordance with previous studies [2,17,48]. Furthermore, the DMD analysis shows inflow turbulence shifted the energy spectrum significantly to lower frequency ranges for both models. These findings are in accordance with previous studies where the wake is analyzed using Fourier Transform [6,49]. Compared with the Fourier Transform, the DMD analysis provides additional insightful information about the spatial scale of the coherent structures. These coherent structures reveal that for the turbulent inflow condition, the wakes are dominated by DMD modes of larger scale coherent structures related to the inflow eddies and some DMD modes are enhanced by the hub vortex. Moreover, a mode similar to the bluff-body vortex shedding at Strouhal number St = 0.17 is found to be dominant uniquely in the wake behind the AD model, which, on the other hand, happens at St = 0.23 for the AS model. The most dominant mode of the AS case appears at lower frequency of St = 0.08 and has larger scale coherent flow structures, which shall be related to the passive advection of the wake by the large scale inflow eddies. The origin of these large-scale motions of turbine wakes is often attributed to two different mechanisms, namely the bluff body shear layer instability [23] and the inflow large eddies [6], which convect turbine wake as passive scalars. Recent field measurements [50] and computational studies [49] suggested the co-existence of these two mechanisms. Furthermore, the hub vortex behind the nacelle is shown to have a significant impact on the start and enhancement of wake meandering [21,51]. A recent review on the meandering of turbine wakes can be found in [52]. In this work we observed a complex interaction between the turbine and the inflow eddies: (i) the modes at low frequencies are less affected by the turbine (Figure 7f); (ii) the DMD modes close to the bluff body vortex shedding frequency seem to be enhanced (Figure 7b,g); (iii) some DMD modes seem to be amplified by the hub vortex behind the nacelle (Figure 7c,h); (iv) the instability within the shear layer also seems to be a key factor for some modes (Figure 7d). Although the present work still can not provide a direct answer to the origin of the wake meandering, it suggests that the dynamic structures are different in wakes computed from the AD and the AS models. Due to this difference, the AD model should be used with more attention when the wake dynamics are of interest, e.g., to study the wake meandering, because the AD model can lead to different wake meandering patterns.

Conclusions
In this work, we evaluate the capability of an actuator disk model in predicting the wake dynamics of a utility-scale wind turbine by comparing its results with those from an actuator surface model. In the AD model, the same thrust coefficient as that in the AS model is employed. A nacelle model is incorporated into both models. Turbulent flows are simulated using LES with the same mesh and time step for both AD and AS cases. Two inflow conditions are considered, i.e., a uniform inflow and a fully developed turbulent inflow. The wakes computed using the AD model is compared with that from the AS model via the time-averaged field and the DMD analysis. It is found that time-averaged velocity and turbulence kinetic energy computed by the AD model are significantly different from those computed by the AS model until nine turbine rotor diameters downstream for the uniform inflow condition; for fully developed turbulent inflow, the differences between the two models are less significant and agree with each other from seven turbine rotor diameters downstream. The DMD analysis of the uniform inflow cases shows that the vortex shed behind the nacelle triggers the shear layer instabilities on the wake boundary behind both models but of different spatial scales. With a thinner shear layer, the wake predicted by the AD model contains smaller spatial scale oscillations at higher frequency. For the fully developed turbulent cases, the DMD analysis shows that the spectra of both models shift to a lower frequency range and the coherent structures also increase in size. The DMD analysis also reveals significant differences between the two models in the far wake: a bluff-body vortex shedding pattern at St = 0.17 appears uniquely in the wake of the AD model as the most dominant DMD mode, whereas the wake computed by the AS model has the most dominant DMD mode of lower energy and at lower frequency St = 0.08 which is related to the passive transport by the inflow turbulence large eddies, and with the second dominant mode at a frequency St = 0.23 close to the bluff body vortex shedding frequency. It is concluded that the dynamic coherent structures in the wake predicted by the AD model are significantly different from those predicted by the AS models and shall be used with more attention when the dynamics of the wake are of interest. In the present work, the thrust coefficient employed in the AD model is the same as that computed by the AS model. However, since the blade rotation is not modeled in the current AD model, the power coefficient from the AD model is not exactly the same as that from the AS model (the C P from the AD is approximately 5% to 10% higher). In the current AD model, the thrust coefficient and the power coefficient cannot be specified at the same time. Further studies on how the differences in the power coefficient affect wake evolutions will be carried out using more advanced AD models considering the effect of blade rotation (e.g., [20]).
Author Contributions: Conceptualization, investigation, writing-original draft preparation, Z.L.; conceptualization, methodology, software and writing-review and editing, X.Y. All authors have read and agreed to the published version of the manuscript.
Funding: This work is partially supported by NSFC Basic Science Center Program for "Multiscale Problems in Nonlinear Mechanics" (NO. 11988102).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

ABL
Atmospheric Boundary Layer AD Actuator Disk AL Actuator Line AS Actuator Surface DMD Dynamic Mode Decomposition St Strouhal Number