Wind Turbine Aerodynamics Simulation Using the Spectral/hp Element Framework Nektar++

Abstract

Wind power plays an increasingly vital role in sustainable energy development. However, accurately simulating wind turbine aerodynamics, particularly in offshore wind farms, remains challenging due to complex environmental factors such as the marine atmospheric boundary layer. This study investigates the integration and assessment of the Actuator Line Model (ALM) within the high-order spectral/hp element framework, Nektar++, for wind turbine aerodynamic simulations. The primary objective is to evaluate the implementation and effectiveness of the ALM by analyzing aerodynamic loads, wake behavior, and computational demands. A three-bladed NREL-5MW turbine is modeled using the ALM in Nektar++, with results compared against established computational fluid dynamics (CFD) tools, including SOWFA and AMR-Wind. The findings demonstrate that Nektar++ effectively captures velocity and vorticity fields in the turbine wake while providing aerodynamic load predictions that closely align with finite-volume CFD models. Furthermore, the spectral/hp element framework exhibits favorable scalability and computational efficiency, indicating that Nektar++ is a promising tool for high-fidelity wind turbine and wind farm aerodynamic research.

1. Introduction

Wind power has become a key contributor to meeting the growing demand for sustainable energy. The aerodynamics of wind turbines involve complex interactions between airflow and turbine blades, which are crucial to determining both the performance and lifespan of wind turbines. Accurate simulation of these dynamics is essential for the optimal design and configuration of wind farms, especially in offshore environments, where additional factors such as wave action, marine atmospheric boundary layers, and platform motion further complicate these interactions.

Historically, lower-order numerical methods have been used to model wind turbine aerodynamics. Although useful, these methods often fall short in capturing intricate aspects of blade loading and wake effects, particularly in tightly clustered turbines within wind farms. Recent advances in rotor blade modeling, particularly through the Actuator Line Model (ALM) and Computational Fluid Dynamics (CFD), as introduced by Sorensen et al. [1] and others [2,3,4,5], have improved our understanding of wind turbine aerodynamics. ALM simplifies the representation of rotor blades by modeling them as lines with distributed forces along their radial direction, allowing detailed and physically accurate time-resolved rotor simulations that have been validated against load distribution and power coefficient comparisons [6]. In addition, ALM has demonstrated considerable reductions in computational expenses by eliminating the need for geometrically detailed blade simulations [7,8], yet still maintaining satisfactory accuracy and detailed flow insights. This methodology has been validated under complex inflow conditions like high shear and yaw misalignment [8,9], where it successfully replicates the performance predictions typically reserved for more computationally demanding blade-resolved simulations. This validation underscores the robustness and versatility of the ALM in simulating realistic operating environments for wind turbines, providing a comprehensive and cost-effective tool for wind energy research and development.

High-order numerical methods used in fluid dynamics, such as spectral/$hp$ element [10,11,12,13], represent a significant advancement in dealing with complex geometries and high-resolution grid meshes. The spectral/$hp$ element method combines the high accuracy of spectral methods along with the geometric flexibility of h-type finite element methods. In this method, a polynomial expansion of order p is applied in each element domain within a finite element mesh of size h, offering a powerful tool for modeling complex flow fields and providing high-fidelity results. In addition, the spectral/$hp$ element method allows for the accurate simulation of critical scales of motion in turbulent flows (even with a relatively small number of elements) more effectively than their lower-order counterparts [14]. This ability is crucial in wind turbine aerodynamics, allowing us to more accurately predict fatigue loads, optimize turbine spacing, and maximize energy extraction from wind sources. In addition, it improves overall wind farm efficiency and turbine durability. Furthermore, the development of high-fidelity CFD simulations not only supports fundamental research but also enhances applied studies by streamlining the virtual evaluation of new concepts and design solutions.

Several high-order numerical tools such as NEK5000 (Mathematics and Computer Science Division of Argonne National Laboratory, Lemont, IL, USA) [15], LESGO (Turbulence Research Group at Johns Hopkins University, Baltimore, MD, USA) [16], and Xcompact3d (Department of Aeronautics at Imperial College London, London UK) [17] have been utilized for accurate modeling of wind turbine aerodynamics. The integration of ALM within these high-order numerical tools enables more accurate representations of the wind turbine aerodynamic loads and wake, facilitating well-informed design and optimization strategies for large wind farms.

Previous studies have applied the ALM in spectral element simulations to effectively model the interaction between rotor blades and complex flow fields. Peet et al. [18] demonstrated the implementation and validation of the ALM in a high-order spectral element code, highlighting its efficiency in reducing computational costs, accuracy in predicting wake dynamics, and scalability for large-scale wind farm simulations. Chatterjee et al. [19] implemented the ALM within spectral element simulations to investigate the dynamics of wind turbines in atmospheric boundary layer (ABL) flows, demonstrating the model’s efficiency in accurately capturing turbulent wake structures while maintaining computational efficiency in high Reynolds number scenarios. Further studies by Chatterjee et al. [20] used ALM in large-eddy simulations with a spectral element framework to analyze turbulent dynamics and power generation in wind turbine arrays within the ABL. Comparative studies, such as those by Kleusberg et al. [21], validated the wake structures behind wind turbines modeled with spectral element codes against finite-volume codes at comparable resolutions, highlighting the accuracy and reliability of the spectral element method. Furthermore, Melani et al. [22] emphasized the spectral/$hp$ element method’s capability to deliver accurate blade loadings, even with coarse computational fluid dynamics (CFD) grids, leading to improved normal and tangential load predictions.

The integration of ALM with spectral/$hp$ element methods presents a significant advancement in aerodynamic simulations, achieving a balance between high accuracy and computational efficiency. This approach transcends traditional models, such as the Blade Element Momentum (BEM) method [23] and the vortex methods [24,25,26], providing the necessary fidelity to capture the intricate aerodynamics of modern wind turbines [27]. High-fidelity simulations using this approach also deepen our understanding of the flow physics around turbines, which simpler models often overlook, and are crucial to optimizing performance, predicting rotor torque [28], and assessing environmental interactions [29].

This study focuses on integrating and evaluating an Actuator Line Model (ALM) within the high-order spectral/$hp$ element framework Nektar++ [30,31] for wind turbine aerodynamic simulations. Nektar++ is an open-source, high-performance, and scalable framework for developing solvers for partial differential equations using advanced high-order spectral/$hp$ element methods. It accommodates both classical low polynomial order h-type solvers (where h denotes finite element size) and high p-order piecewise polynomials, offering versatility for various applications. In the context of flow simulations by solving the Navier–Stokes equations, Nektar++ models velocity fields using polynomials of order p, utilizing Lagrange interpolants at Gauss–Lobatto–Legendre (GLL) points, with $N=p+1$ GLL points per spatial direction. Pressure fields are represented using ($p-1$) Gauss–Legendre (GL) points, ensuring accurate and efficient simulations. Notably, the current ALM implementation within Nektar++ has certain limitations; it does not account for dynamic stall effects, which are common in wind turbines operating under variable conditions. Moreover, the blades are treated as rigid, excluding the influence of structural deformations on aerodynamic loading.

The primary objective of this work is to examine the implementation and performance of the ALM within the Nektar++ framework, aiming for high-fidelity wind turbine aerodynamic simulations. The research focuses on assessing the aerodynamic loads, wake characteristics, and computational efficiency of the Nektar++ framework in comparison with established CFD tools like SOWFA and AMR-Wind. The remainder of this paper is organized as follows: Section 2 describes the numerical and modeling methodologies, including the ALM approach for rotor blade representation. Section 3 outlines the computational setup, test cases, and simulation configurations. Section 4 presents the results, focusing on performance metrics, parametric studies, aerodynamic loading, and wake characteristics. Finally, Section 5 offers a detailed discussion of the findings, highlighting the advantages and limitations of the Nektar++ framework and suggesting potential directions for future research.

2. Methodology

2.1. Navier–Stokes Equations

The incompressible Navier–Stokes equations are expressed as

$$\begin{array}{cc}\hfill \nabla ·\mathbf{u}& =0\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+(\mathbf{u}·\nabla )\mathbf{u}& =-{\displaystyle \frac{1}{\rho }}\nabla p+\nu {\nabla }^2\mathbf{u}+{\displaystyle \frac{\mathbf{f}}{\rho }}\hfill \end{array}$$

(2)

where $\mathbf{u}$, t, p, $\rho$, $\nu$, and $\mathbf{f}$ denote the velocity field, time, pressure field, fluid density, kinematic viscosity, and external forces, respectively.

In the spectral/$hp$ element method for solving incompressible Navier–Stokes equations, the solution variables within each element are approximated using a polynomial expansion. The number of modes (n) typically refers to the degree of the polynomial basis functions used within each element of the mesh. Nektar++ utilizes polynomials of order p to model velocity fields. However, for pressure fields, Gauss–Legendre (GL) points of order $p-1$ are used. Furthermore, the number of modes, representing $n=p+1$, corresponds to the polynomial order used in the solution approximation.

In the present study, a large eddy simulation (LES) approach is adopted, in which only the large scales are solved for while the small scales are included by means of a subgrid scale (SGS) model. For continuous Galerkin implicit LES, it is typically performed by applying spectral vanishing viscosity or using a gradient jump penalty; see e.g., the discussion in [32]. However, in this work, the explicit filter of [33] is used, which effectively acts as a sub-grid scale model. The filter essentially removes some of the energy contained in the least resolved modes, typically 5% of the top mode.

2.2. Rotor Blades Modeling

In this study, the turbine blades are represented using the actuator line model (ALM) as shown in Figure 1, which effectively simulates the aerodynamic forces on the blades without explicitly resolving the blade geometry [1].

The ALM is a computational technique that enhances the traditional actuator disk model (ADM) by representing the blades as individual lines rather than simplifying the entire rotor into a single disk representation. By modeling blades as lines of aerodynamic forces, the ALM allows for a detailed analysis of these forces along the blade length while avoiding the computational burden of fully resolving the geometry. In practical application, each blade is discretized into several segments called actuator points (see Figure 1). The aerodynamic forces at each actuator point are computed using the blade element theory utilizing tabulated airfoil data based on local flow conditions. These forces are incorporated into the incompressible Navier–Stokes equations as the body force per unit volume through the source terms ($\mathbf{f}$) as illustrated in Equation (2).

As seen in Figure 1, the local velocity relative to the rotating blade ($U_{rel}$) is calculated using a velocity triangle, where $u_n$ and $u_{\theta }$ represent the normal and circumferential velocity components with respect to the rotor plane, respectively. The tangential velocity is also expressed as

$$u_t=\mathsf{\Omega }r-u_{\theta }$$

(3)

where $\mathsf{\Omega }$, r, and $u_{\theta }$ denote the rotational velocity of the blade, the radius of the rotor section, and the local circumferential velocity, respectively. The local flow angle ($\varphi$) and the angle of attack ($\alpha$) are then given by

$$\begin{array}{cc}\hfill \varphi & =ta{n}^{-1}(u_n/u_t)\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \alpha & =\varphi -({\theta }_p+{\theta }_t)\hfill \end{array}$$

(5)

where ${\theta }_p$ and ${\theta }_t$ are the pitch and twist angles of the blade section, respectively. With the angle of attack and relative velocity known, the two-dimensional lift and drag forces per unit span are calculated as

$$\begin{array}{cc}\hfill {\mathbf{L}}^{\prime }& =\rho c{\mathbf{u}}_{rel}^2\left(C_L{\mathbf{e}}_L\right)/2\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\mathbf{D}}^{\prime }& =\rho c{\mathbf{u}}_{rel}^2\left(C_D{\mathbf{e}}_D\right)/2\hfill \end{array}$$

(7)

Here, c, $Re$, ${\mathbf{e}}_L$, and ${\mathbf{e}}_D$ are the chord length, the Reynolds number, and the unit vectors in the directions of lift and drag, respectively. Moreover, the aerodynamic coefficients, $C_L=C_L(\alpha ,Re)$ and $C_D=C_D(\alpha ,Re)$, are dimensionless parameters used to quantify the lift and drag forces experienced by a blade section, respectively. These coefficients are derived from tabulated airfoil data and are dependent on the angle of attack ($\alpha$) and the Reynolds number ($Re$), both of which are determined based on local flow conditions. To account for three-dimensional flow effects, the standard Glauert correction factor ($f_G$) for tip and root losses [23] is applied using

$$\begin{array}{cc}\hfill f_{tip}\left(r\right)& ={\displaystyle \frac{2}{\pi }}{cos}^{-1}\left(exp\left(-{\displaystyle \frac{b\left(R_{tip}-r\right)}{2rsin\left(\varphi \right)}}\right)\right)\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill f_{root}\left(r\right)& ={\displaystyle \frac{2}{\pi }}{cos}^{-1}\left(exp\left(-{\displaystyle \frac{b\left(r-R_{root}\right)}{2rsin\left(\varphi \right)}}\right)\right)\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill f_G\left(r\right)& =f_{tip}\left(r\right)·f_{root}\left(r\right)\hfill \end{array}$$

(10)

where b, $R_{tip}$, $R_{root}$, r, and $\varphi$ denote the number of blades, rotor tip radius, rotor root (hub) radius, local blade radius, and local flow angle, respectively.

The forces computed at each actuator point along the $i^{th}$ blade are projected into the normal and tangential directions with respect to the rotor plane.

$${\mathbf{f}}_i^{2\mathrm{D}}\left(r\right)=\left(f_G\left(r\right)f_i^n\left(r\right),f_G\left(r\right)f_i^t\left(r\right)\right)$$

(11)

These forces are distributed to the spectral elements through a three-dimensional Gaussian regularization kernel ($\eta$) as

$$\begin{array}{cc}\hfill \mathbf{f}\left(\mathbf{x}\right)& =-\sum _{i=1}^{N_b}{\int }_0^R\eta \left(d\right){\mathbf{f}}_{\mathbf{i}}^{2\mathrm{D}}\left(r\right)dr\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \eta \left(d\right)& ={\displaystyle \frac{1}{{\epsilon }^3{\pi }^{3/2}}}exp\left(-{\displaystyle \frac{|\mathbf{x}-r{\mathbf{e}}_i{|}^2}{{ϵ}^2}}\right)\hfill \end{array}$$

(13)

where $N_b$, R, d, ${\mathbf{e}}_i$, and $ϵ$ denote the number of blades, rotor radius, distance between actuator line coordinates (i.e., $r{\mathbf{e}}_i$) and quadrature points (i.e., $\mathbf{x}$), blade direction unit vector, and the kernel spreading width, respectively. Several studies [2,4,18,34,35,36,37] that employed finite volume-based CFD solvers suggest setting the parameter $ϵ$ within the range of $[{\mathsf{\Delta }}_{grid},4{\mathsf{\Delta }}_{grid}]$. ${\mathsf{\Delta }}_{grid}$ is defined as

$${\mathsf{\Delta }}_{grid}={\left({\mathsf{\Delta }}_x{\mathsf{\Delta }}_y{\mathsf{\Delta }}_z\right)}^{1/3}$$

(14)

representing the average size of the grid cells in the x, y, and z directions. In a spectral/$hp$ element framework, ${\mathsf{\Delta }}_{grid}$ may be defined as the average spacing between quadrature points and can be approximated by

$${\mathsf{\Delta }}_{grid}={\displaystyle \frac{{\left({\mathsf{\Delta }}_x{\mathsf{\Delta }}_y{\mathsf{\Delta }}_z\right)}^{1/3}}{p+1}}$$

(15)

where p is the polynomial order [38].

3. Test Case

In this study, the aerodynamic performance of the three-bladed NREL-5MW wind turbine rotor is analyzed, serving as a widely used and validated benchmark case for wind turbine simulations [39]. The simulations are conducted at a Reynolds number of ${10}^4$, calculated based on the rotor radius ($R{e}_R=1\times {10}^4$). This Reynolds number is significantly lower than that encountered in operational wind turbines. However, previous studies have demonstrated that for aerodynamic simulations of wind turbines, the flow behavior becomes relatively insensitive to Reynolds number effects above a threshold value of approximately $Re=1\times {10}^3-1\times {10}^4$ [35,40,41]. Hence, the selected $R{e}_R=1\times {10}^4$ is sufficient to capture the primary aerodynamic characteristics of turbine blades without significantly compromising the accuracy of predictions related to lift and drag forces, as well as wake dynamics [18].

3.1. Computational Domain

Like all Computational Fluid Dynamics (CFD) simulations, Actuator Line Model (ALM) simulations require selecting appropriate domain and grid sizes to ensure accurate results. These selections must also align with the available hardware capabilities to maintain a balance between simulation fidelity and computational feasibility.

Figure 2 presents a schematic representation of the computational domain and the positioning of the turbine. The dimensions of the domain are specified as $L_x=12R$, $L_y=12R$, and $L_z=12R$ in the streamwise, vertical, and spanwise directions, respectively, where R denotes the radius of the turbine rotor. The distance between the inlet and the turbine is set to $2R$, ensuring that the inlet boundary remains unaffected by the turbine’s influence. The NREL-5MW turbine is characterized by a hub height of 90 m and a rotor radius of 63 m. The turbine hub height is positioned centrally within both the vertical and lateral dimensions of the domain, providing adequate space to accurately capture the flow dynamics surrounding the turbine.

To facilitate comparability with other state-of-the-art simulation tools for verification purposes, the grid is uniformly distributed throughout the computational domain. It consists of 9261 elements arranged in a $21\times 21\times 21$ configuration, covering the domain’s width, height, and length, respectively. The grid employs a Cartesian coordinate system with equidistant spectral elements, supporting high-resolution simulations. Each element is discretized using seventh-order polynomials ($p=7$) with 12 Gauss–Lobatto–Legendre (GLL) points in each spatial direction to reduce aliasing errors. This configuration results in approximately 16 million quadrature points, providing a high-resolution representation of the flow field. Such high-order discretization is essential for accurately capturing complex flow structures, particularly in the near-wake region and around the actuator lines, which are critical areas in wind turbine aerodynamic simulations.

3.2. Actuator Line Model (ALM)

In this study, each actuator line is discretized into 80 nodes along the blade radius, providing fine resolution of the blade forces and their effects on the surrounding flow. The current implementation of ALM does not account for the wind turbine tower or nacelle, focusing only on the blade forces within $0.1<r<1$, where r represents the normalized radial position. The turbine rotates at $9.15$ rpm, corresponding to an optimal tip speed ratio of $7.55$. Furthermore, a time step of $0.01$ seconds ensures that each actuator line advances no more than one grid point per time step, optimizing both accuracy and efficiency. The simulations are conducted until a statistical steady-state condition is achieved, which occurs after 14 rotor revolutions. This duration corresponds to a wake propagation distance of approximately five rotor diameters (630 m) downstream.

3.3. Boundary Condition

At the inlet boundary, a Dirichlet condition is imposed, prescribing a steady and uniform inflow velocity of 8 m/s. The outlet boundary condition is treated as a pressure outlet, allowing the flow to exit the domain freely while maintaining a zero-gradient condition for the streamwise velocity. This boundary is located at a sufficient distance (8R) from the rotor to ensure a minimal influence on the blade loading. Periodic boundary conditions are implemented in the vertical (y) and spanwise (z) directions. In addition, the effects of the nacelle and tower on the flow are neglected.

3.4. Simulation Setup

The Nektar++ framework is configured to define the numerical setup for the simulation. The equation type is specified as UnsteadyNavierStokes, suitable for modeling fluid dynamics problems with time-dependent behavior. The solver type is set to VelocityCorrectionScheme for efficiently solving the Navier–Stokes equations. The global system is solved using the IterativeStaticCond method, while the advection form is configured as Convective, a common choice for fluid flow problems.

For time integration, the IMEXOrder2 scheme is employed, providing a second-order implicit–explicit approach balancing stability and accuracy in time-dependent simulations. Additionally, the LowEnergyBlock pre-conditioner is used to enhance the convergence rate of the iterative solver, as recommended in previous studies [30]. The IterativeSolverTolerance, which defines the relative stopping tolerance for the iterative solver, is set to $1.0\times {10}^{-9}$ for both the velocity and pressure fields.

The Actuator Line Model (ALM) is implemented in Nektar++ by using the framework’s ability to introduce external forces into the governing equations. In Nektar++, external forces are incorporated during the assembly of the forcing terms in the momentum equations. This is achieved by modifying the right-hand side of the discretized Navier–Stokes equations to include the desired external force components. These components are implemented through the Forcing mechanism, facilitating the direct application of forces at specified locations or along specific geometrical features, such as actuator lines.

The Forcing functionality in Nektar++ is part of the SolverUtils module, which provides essential data structures and algorithms for developing solvers or auxiliary functionalities. Using the Forcing approach, the ALM imposes aerodynamic forces along the actuator lines, mimicking the effect of blades in wind turbines or similar systems. These forces are added to the momentum equations during the solver’s assembly process.

4. Results

The results obtained from the simulation of wind turbine aerodynamics using the Nektar++ spectral/$hp$ element framework are presented in this section. Some key performance indicators, such as power and thrust coefficients, along with measures of computational efficiency, such as runtime, are analyzed to assess the effectiveness of the proposed approach. Furthermore, a parametric study is performed to examine the impact of various numerical configurations on the aerodynamic forces acting on the turbine blades.

The ALM implementation in Nektar++ is also compared against established models for wind turbine aerodynamic simulations, including the finite-volume CFD solvers SOWFA [3] and AMR-Wind [42]. To ensure a fair comparison, the simulation setup used for Nektar++ is replicated as closely as possible in SOWFA and AMR-Wind. Key simulation parameters, such as computational domain size, boundary conditions, time step size, fluid properties, and ALM-specific parameters, are aligned across all platforms. Furthermore, to maintain consistency in grid resolution, SOWFA and AMR-Wind are configured with a grid resolution of 256 cells along the x, y, and z directions, respectively, aligned with the grid discretization in Nektar++, which used approximately 16 million quadrature points.

4.1. Performance Metrics of Nektar++

Table 1. Aerodynamic performance metrics and computational time for varying polynomial orders (p) and smearing factors ($ϵ$) in NEKTAR++ simulations.

Item	$\mathit{p}=4$			$\mathit{p}=5$			$\mathit{p}=6$			$\mathit{p}=7$
	$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$		$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$		$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$		$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$
$C_P$ [-]	0.548	0.551		0.551	0.556		0.551	0.557		0.552	0.556
$C_T$ [-]	0.789	0.789		0.791	0.792		0.790	0.792		0.791	0.792
t [CPU-h]	4045	4102		5911	6111		10,293	10,978		16,863	17,369

demonstrates how aerodynamic performance metrics, namely the dimensionless power coefficient ($C_P$) and the thrust coefficient ($C_T$), along with the computational time (t) in CPU-hours, change with varying polynomial orders (p) and smearing factors ($ϵ$) in NEKTAR++ simulations.

As observed, increasing the polynomial order p from 4 to 7 results in the convergence of both $C_P$ and $C_T$, indicating enhanced simulation fidelity with higher polynomial orders. In particular, $C_P$ increases by approximately $0.5\%$ to $1.0\%$, with slightly higher values recorded at $ϵ=8$ than at $ϵ=6$. Likewise, $C_T$ shows a consistent increase of about $0.25\%$ as $ϵ$ increases from 6 to 8. However, it should be noted that computational time increases significantly with increasing p, highlighting the inherent trade-off between achieving higher accuracy and incurring greater computational costs.

4.2. Parametric Study in Nektar++

Figure 3 illustrates the distribution of the normal forces ($F_n$) and tangential ($F_t$) forces per blade span along the normalized radial positions ($\mathrm{r}/\mathrm{R}$) of the rotor blade, where $\mathrm{r}$ represents the radial position and $\mathrm{R}$ denotes the radius of the blade. The plot highlights the effects of varying polynomial orders (p) and the smearing factors ($ϵ$) on the distribution of forces.

As the polynomial orders increase, the normal force profile exhibits convergence, as shown in the inset of Figure 3a, which focuses on the region of the tip of the blade ($r/R=0.78$ to $0.90$). For both $ϵ=6$ and $ϵ=8$, the differences in $F_n$ between $p=4$ and $p=7$ are less than $1\%$, indicating that the contributions of higher modal values offer only marginal improvements in accuracy. Furthermore, increasing $ϵ$ from 6 to 8 slightly increases the peak force and smooths the distribution throughout the radial range, although this change remains marginal.

Figure 3b shows the profiles of tangential force per blade span ($F_t$) along the normalized radial positions. Similarly to the normal force ($F_n$), the tangential force demonstrates convergence with increasing polynomial order (p). The sensitivity of $F_t$ to the smearing factor ($ϵ$) is more pronounced near the tip of the blade ($0.8<r/R<0.9$) and in the root region ($0.2<r/R<0.3$), the regions characterized by sharp gradients, compared with the mid-span of the blade. For $ϵ=6$, increasing p from 4 to 7 results in an approximate increase of $1.2\%$ in the peak $F_t$ within the blade tip region. For $ϵ=8$, the increase is slightly lower, at approximately $1.02\%$.

These observations suggest that both p and $ϵ$ influence the normal and tangential force distributions along the rotor blades. Polynomial order (p) has a more pronounced effect on improving the accuracy and convergence of force profiles, while the smearing factor ($ϵ$) mainly affects the smoothness of the distribution, reducing sharp gradients without significantly altering the overall force magnitudes.

4.3. Aerodynamic Loads

Accurate prediction of aerodynamic forces on wind turbine blades is essential for the design, analysis, and optimization of wind turbines. CFD tools play a critical role in capturing aerodynamic forces with high accuracy, providing a comprehensive understanding of the complex flow around the rotor blades. A comparative analysis of aerodynamic load predictions is presented, highlighting results from Nektar++ alongside those from SOWFA and AMR-Wind.

Figure 4 illustrates the comparison of normal ($F_n$) and tangential ($F_t$) forces per blade span along the normalized blade radius ($r/R$) as predicted by various simulation tools. The results of Nektar++ are displayed for a simulation configuration using a polynomial order of 7 ($p=7$) and a smearing factor of 6 ($ϵ=6$).

As shown, the predicted normal force ($F_n$) from Nektar++ exhibits good agreement with those from SOWFA and AMR-Wind across most of the blade span. A similar trend is observed among the computational tools up to approximately $r/R=0.8$, where the normal force gradually increases with the radial position. Beyond this region, the computational tools begin to diverge slightly, particularly near the blade tip ($r/R\approx 0.9$), where AMR-Wind predicts a higher peak force compared with both Nektar++ and SOWFA. This difference may be linked to variations in numerical schemes and grid resolution, particularly in regions with strong velocity gradients, such as near the blade tip, where aerodynamic complexities are more pronounced. Moreover, Nektar++ predicts higher tangential forces in the mid-span region ($0.3<r/R<0.7$) compared with SOWFA and AMR-Wind. However, near the blade root ($r/R<0.2$) and the tip ($r/R>0.9$), Nektar++ shows good agreement with the SOWFA results. The uneven profile of the tangential force distribution observed in the results may be attributed to the coarse grid resolution employed in this study. Enhancing mesh refinement, along with applying induced velocity corrections [43], could potentially mitigate this discrepancy and improve agreement between different computational tools.

This comparative analysis highlights Nektar++’s capability to deliver high-fidelity aerodynamic simulations for wind turbine analysis, demonstrating performance comparable to well-established finite-volume methods like SOWFA and AMR-Wind, while also offering potential advantages in computational efficiency.

4.4. Wake Characteristics

To evaluate the accuracy of the Actuator Line Model (ALM) implemented in Nektar++, a comparative analysis of the time-averaged velocity and vorticity fields is conducted against those from SOWFA and AMR-Wind. Similarly to aerodynamic loads, the results from Nektar++ are presented for a simulation configuration using a polynomial order of 7 ($p=7$) and a smearing factor of 6 ($ϵ=6$).

Figure 5 displays the time-averaged velocity fields at hub height in the horizontal ($xy$) plane (Figure 5a–c) and the vertical ($yz$) plane at the rotor plane (Figure 5d–f). All three simulations demonstrate a similar turbine wake downstream of the rotor, characterized by a velocity deficit directly behind the rotor and higher velocities at the wake edges. Although slight differences in wake width are observed among the tools, these differences remain within an acceptable range.

To further evaluate Nektar++’s capability in modeling wind turbine rotor aerodynamics, a comparative analysis of the time-averaged vorticity magnitude is conducted. Figure 6 shows the results obtained from Nektar++, SOWFA, and AMR-Wind in both the horizontal ($xy$) plane at hub height and the vertical ($yz$) plane at the rotor.

In the horizontal plane (Figure 6a–c), all three simulation tools exhibit similar qualitative trends, with the highest vorticity magnitude observed near the blade tip and root regions. However, minor variations in intensity and downstream diffusion of vortex structures are apparent. In the vertical plane (Figure 6d–f), distinct vortex structures associated with the blade tip and root are clearly visible for each simulation tool. Despite the qualitative similarities, differences in the spatial extent and intensity of these vortices are noticeable. These discrepancies are likely the result of variations in numerical methods, grid resolutions, or other modeling parameters used by each tool.

4.5. Velocity Deficit

Figure 7 presents the time-averaged streamwise velocity (U) profiles plotted against the normalized vertical coordinate ($z/D$) at four distinct streamwise locations ($x/D$ = 0, 1, 2, and 3). Here, D represents the turbine diameter, while x and z denote the axial and vertical directions, respectively. Each subplot provides a comparison between the results from Nektar++ (with $p=7$ and $ϵ=6$), SOWFA, and AMR-Wind.

Across all streamwise locations, a consistent pattern appears in the wake region, characterized by a velocity deficit, which gradually recovers toward the outer edges, where velocities approach the free-stream value of approximately 8 m/s. The three solvers generally show good agreement with only minor discrepancies. These differences are most pronounced at $x/D$ = 2 and $x/D$ = 3, where Nektar++ slightly underestimates the velocity deficit near the rotor tip compared with SOWFA and AMR-Wind.

4.6. Aerodynamic Performance

Table 2. Comparison of aerodynamic performance metrics and computational time across various computational tools with different fidelity levels.

Item	Nektar++ (p = 5)			Nektar++ (p = 7)			AMR-Wind			SOWFA
	$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$		$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$		$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$		$\mathit{ϵ}=\mathbf{6}$	$\mathit{ϵ}=\mathbf{8}$
$C_P$ [-]	0.551	0.556		0.552	0.556		0.575	0.585		0.527	0.549
$C_T$ [-]	0.791	0.792		0.791	0.792		0.831	0.836		0.789	0.810
t [CPU-h]	6111			16,863			8557			10,384

compares aerodynamic performance metrics, such as the dimensionless power coefficient ($C_P$) and the thrust coefficient ($C_T$), as well as computational demands in CPU-hours, in three different CFD tools: Nektar++ (with polynomial orders $p=5$ and $p=7$), AMR-Wind, and SOWFA. It also illustrates the impact of varying the smearing factor ($ϵ$) on $C_P$ and $C_T$ for each of these tools.

The ALM simulation in Nektar++ demonstrates consistent performance across varying smearing factors ($ϵ$). For $p=5$, $C_P$ increases slightly from $0.551$ to $0.556$ (approximately $0.9\%$), while $C_T$ remains relatively unchanged. At a higher polynomial order of $p=7$, similar trends are observed. AMR-Wind exhibits slightly greater sensitivity to changes in $ϵ$, with $C_P$ increasing from $0.575$ to $0.585$ (about $1.7\%$) and $C_T$ from $0.831$ to $0.836$ (about $0.6\%$). SOWFA displays the highest relative sensitivity, with $C_P$ increasing from $0.527$ to $0.549$ ($\sim 4.2\%$) and $C_T$ from $0.789$ to $0.810$ ($\sim 2.7\%$) as $ϵ$ increases. The agreement in trend and magnitude between Nektar++, AMR-Wind, and SOWFA reinforces Nektar++’s capability as a reliable high-fidelity tool for wind turbine aerodynamic analysis.

In terms of computational cost, significant differences are observed. Nektar++ with $p=7$ has the highest computational demand, requiring 16863 CPU-hours, compared with 8557 CPU-hours for AMR-Wind and 10384 CPU-hours for SOWFA. However, reducing the polynomial order to $p=5$ in Nektar++ lowers the computational cost to 6111 CPU-hours, representing approximately $29\%$ and $41\%$ reductions relative to AMR-Wind and SOWFA, respectively. Despite this reduction, Nektar++ still delivers comparable accuracy. This demonstrates that by carefully selecting the polynomial order, Nektar++ can achieve an optimal balance between computational efficiency and the accuracy of rotor aerodynamic simulations.

5. Discussion

The integration of the Actuator Line Model (ALM) within the Nektar++ spectral/$hp$ element framework reveals substantial potential for high-fidelity simulations in wind turbine aerodynamics. The results of these simulations provide a deeper understanding of both numerical methodologies and wind turbine aerodynamics.

Performance metrics within the Nektar++ framework demonstrate an accurate representation of aerodynamic loads, particularly in comparison to finite-volume computational fluid dynamics (CFD) solvers such as SOWFA and AMR-Wind. Notably, the power coefficient ($C_P$) and thrust coefficient ($C_T$) results from Nektar++ closely align with those from CFD models, falling within acceptable operational ranges for wind turbines. This alignment validates the ALM implementation within the Nektar++ spectral/$hp$ element framework. Additionally, we observe a slightly lower power coefficient from Nektar++ relative to AMR-Wind, a leading finite-volume-based CFD solver for wind turbine simulations. This difference may stem from variations in grid resolution and numerical methodologies. Further refinement in mesh resolution or induced velocity adjustments, as suggested by Dağ et al. [43], may improve the alignment between the tools.

Comparisons of numerical parameters and simulation results are summarized in Table 1. The incremental increases in power and thrust coefficients as the polynomial order increases from $p=4$ to $p=7$ suggest that higher-order polynomial approximations capture more refined flow details. However, the marginal improvements observed beyond $p=6$ indicate a potential convergence point, where further increases in polynomial order may not justify the associated computational expense.

The sensitivity of the simulation results to the smearing factor ($ϵ$) deserves special attention. The moderate variations in power coefficient between $ϵ=6$ and $ϵ=8$ highlight the sensitivity of spectral element ALM implementations to force projection parameters, consistent with findings reported for finite-volume-based CFD solvers [2,4,18,34,35,36,37].

Accurate modeling of wake characteristics is crucial for wind turbine simulations, as these dynamics significantly influence the performance of downstream turbines in a wind farm. Our study demonstrates that Nektar++ generates wake profiles comparable to those of SOWFA and AMR-Wind, particularly regarding velocity deficit and vorticity.

The spectral/$hp$ element method employed in Nektar++ appears to provide a balance between computational efficiency and accuracy. Our findings show that computational time scales with the polynomial orders (p), providing valuable insights into optimal configurations that balance resource utilization with simulation fidelity. For example, with the polynomial order of five ($p=5$), the computational time was reduced by approximately $29\%$ and $41\%$ relative to AMR-Wind and SOWFA, respectively, without a compromise in accuracy. This reduction in computational costs supports high-order spectral methods as a practical option for large-scale simulations, which are essential for studying complex inflow conditions and wake interactions, especially in offshore wind farms.

Although this study demonstrates the effectiveness of the Nektar++ framework for wind turbine aerodynamics, future research should extend these findings by incorporating a wider range of test cases. These may include diverse atmospheric boundary layer conditions and multi-turbine configurations, enabling a more comprehensive evaluation of the framework’s potential for wind energy applications. For example, the current implementation of ALM in Nektar++ does not account for dynamic stall, and including dynamic stall modeling could improve the prediction accuracy for turbines operating in unsteady environments. Moreover, incorporating blade flexibility in ALM simulations would provide a more realistic representation of large wind turbines, aligning simulations more closely with operational conditions. Furthermore, the validation of Nektar++ in this study is conducted through comparisons with CFD tools, SOWFA and AMR-Wind, both of which have been thoroughly validated in previous research. However, future work may also involve comparisons with experimental data and high-fidelity blade-resolved simulations to ensure a more comprehensive assessment of the framework’s accuracy and reliability. In conclusion, the successful implementation of ALM within Nektar++ underscores the potential of high-order methods in wind energy applications, serving as a reliable and efficient alternative to traditional CFD solvers. This is especially relevant as the wind energy industry advances toward larger turbines and more complex wind farm designs, where accurate wake predictions and load calculations are essential.