Numerical investigation of regenerative wind farms featuring enhanced vertical energy entrainment

Numerical simulations of wind farms consisting of innovative wind energy harvesting systems are conducted. The novel wind harvesting system is designed to generate strong lift (vertical force) with lifting devices. It is demonstrated that the trailing vortices generated by these lifting devices can substantially enhance wake recovery rates by altering the vertical entrainment process. Specifically, the wake recovery of the novel systems is based on vertical advection processes instead of turbulent mixing. Additionally, the novel wind energy harvesting systems are hypothesized to be feasible without requiring significant technological advancements, as they could be implemented as multi-rotor systems with lifting devices (MRSLs), where the lifting devices consist of large airfoil structures. Wind farms with these novel wind harvesting systems, namely MRSLs, are termed regenerative wind farms, inspired by the concept that the upstream MRSLs actively entrain energy for the downstream ones. With the concept of regenerative wind farming, much higher wind farm capacity factors are anticipated. Specifically, the simulation results indicate that wind farm efficiencies can be nearly doubled by replacing traditional wind turbines with MRSLs under the tested conditions, and this disruptive advancement can potentially lead to a profound reduction in the cost of future renewable energy.

1 Introduction

In the wind energy industry, wind turbines are often arranged in clusters, leveraging closer spacings for economic and operational benefits (Meyers and Meneveau, 2012; Sørensen and Larsen, 2021). These clusters are known as wind farms. However, densely packed wind turbines result in annual energy production (AEP) losses due to the turbine–turbine wake interactions. The more tightly packed the turbines, the more pronounced the negative impact on AEP (Stevens et al., 2016). These losses are substantial, with reported AEP reductions ranging from 10 % to 25 % for large-scale offshore wind farms such as Horns Rev 1 and Nysted (Barthelmie et al., 2009, 2010). Moreover, predictions indicate that AEP losses due to wakes could reach more than 60 % for wind farms on a very large scale (infinite wind farm) with spacings similar to the typical ones (e.g., 7D in streamwise and 5D lateral directions, where D is the rotor diameter) (Dupont et al., 2018; Calaf et al., 2010). Note that the above-mentioned considerations are for conventional wind farms that consist of three-bladed horizontal-axis wind turbines (HAWTs), the prevailing concept in today's commercial wind farms (Manwell et al., 2010).

The AEP drop mentioned in the previous paragraph is attributed to the fact that the kinetic energy carried by the incoming wind is depleted by upstream turbines, and the energy replenishing rates cannot sustain the downstream turbines to extract as much energy as those in the first row of a wind farm (Porté-Agel et al., 2020). Note that the energy is mainly replenished by entraining from above the wind farms. This is due to the fact that wind farms are built close to the ground or sea surface, and they extend in both streamwise and lateral directions. However, without significant mean vertical flow in conventional wind farms, the primary source of vertical energy (momentum) entrainment is through the turbulent mixing process, relying on the Reynolds stress terms (Calaf et al., 2010; VerHulst and Meneveau, 2015). Typical rates of vertical energy entrainment are about 1 to 2 W m⁻² for conventional wind farms with HAWTs (Dupont et al., 2018; VerHulst and Meneveau, 2015) (estimated based on an infinite-wind-farm scenario), which is significantly lower than 7 W m⁻², a typical installation capacity (these estimations are based on a typical wind farm; i.e., the ranges of streamwise and lateral spacings are around 7D and 5D, freestream wind speed is estimated at 10 m s⁻¹, and the power coefficient of the turbines is 0.54) (Barthelmie et al., 2009; Bosch et al., 2019). This indicates that the efficiencies of large wind farms with conventional designs are limited by the low vertical entrainment rates.

Figure 1 Conceptualizing the proposed design of the innovative wind energy harvesting systems, namely multi-rotor systems with lifting devices (MRSLs). (a) MRSL consisting of vertical-axis wind turbines (VAWTs) with airfoils/wings that lift the wake upward. (b) MRSL consisting of horizontal-axis wind turbines (HAWTs) with airfoils/wings that lift the wake downward.

To overcome the aforementioned limitations of the current conventional wind farms, we adopt the strategy of introducing lifting devices onto the wind energy harvesting systems. These lifting devices can induce strong vertical flows, leading to a significant vertical advection process and thus enhancing vertical energy entrainment. To the authors' best knowledge, this concept was first studied by Bader et al. (2018), where they carried out numerical analysis of HAWTs coupled with lifting devices close to their rotors in various configurations. Their promising results showed that the power performance of the downstream turbines was substantially improved with the implementation of the lifting devices. However, they did not propose a way to install the lifting devices, as they were suspended without support in the computational domain.

Very recently, Broertjes et al. (2024) and Avila Correia Martins et al. (2025) have also studied the concept of introducing lifting devices onto wind harvesting devices both experimentally and numerically. These two studies were based on the idea proposed by Ferreira et al. (2024). Unlike Bader et al. (2018), an innovative design, the multi-rotor system with lifting devices (MRSLs), was proposed. The system comprises several sub-rotors, each in the form of vertical- or horizontal-axis wind turbines (VAWTs or HAWTs). The proposed design, illustrated in Fig. 1, highlights the designated positions to mount the lifting devices, where the airfoils/wings themselves serve as structural components. This MRSL design has been realized at a wind tunnel scale by Broertjes et al. (2024), featuring a system with 16 VAWTs and two wings. The results of Broertjes et al. (2024) and Avila Correia Martins et al. (2025) showed that, due to the strong vertical flow induced by the lifting devices, the wake recovery rate of MRSLs can reach more than 90 % at a distance of 5D downstream (based on available power, which is ∝u³), whereas a typical HAWT achieves less than 40 % at a similar distance (Li et al., 2024a), indicating a significant enhancement in wake recovery. Additionally, it should be noted that although the concept of MRSLs came out very recently, the implementation of this design may not require major technological breakthroughs, as the technology for multi-rotor systems already exists (Jamieson and Branney, 2012; Watson et al., 2019).

Building on the work of Broertjes et al. (2024) and Avila Correia Martins et al. (2025), this study further investigated the aerodynamics of wind farms consisting of MRSLs using a numerical method. These wind farms are termed regenerative wind farms by Ferreira et al. (2024). The name reflects the idea that upstream MRSLs actively entrain energy for the downstream ones. At this point, it is suggested that the proposed MRSLs and the concept of regenerative wind farms could be a groundbreaking concept for the wind energy industry. This concept has the potential to revolutionize wind energy by fundamentally altering the process of vertical energy entrainment. Unlike conventional wind farms, regenerative wind farms replenish flow energy vertically through the mean components of the flow rather than relying on Reynolds stress terms, which is likely to significantly elevate their wind farm efficiency. If successfully implemented, this approach promises not only significant economic advantages but also a reduction in the space required to generate the same power output compared to conventional wind farms. Achieving these goals could enhance the benefits of wind energy while minimizing its environmental and spatial impacts, marking a transformative advancement in renewable energy. To validate the groundbreaking potential of MRSLs in transforming the vertical entrainment process, this study conducts a comprehensive numerical analysis of regenerative wind farms, setting the stage for a significant leap forward in wind farm efficiency.

Figure 2 Isosurfaces of the magnitude of streamwise vorticity ($\left|{\mathit{\omega }}_x\right|$, silver) together with the contour plots of the vertical velocity fields w. The arrows on the contours depict the direction of the in-plane velocity, and note that the lengths of the arrows are scaled by the square root of the in-plane velocity's norms. MRSLs are represented with red and blue surfaces, where the blue surfaces are thrusting devices, and the red surfaces are the lifting devices. The plots are based on the solutions of the up-washing (a) and down-washing (b) cases in Table 3. The MRSLs depicted are the ones in the fourth row, and the contours are plotted at $x/D=\mathrm{22.0}$.

Figure 3 A sketch of the simplified vortex/circulation system of a wing is presented. The vortex system is depicted as a horseshoe vortex, indicated by the blue line with angles, with the direction of the circulation Γ shown by an arrow. The solid lines in the horseshoe vortex represent the trailing vortices, while the dashed line represents the bound vortex. Note that, with the orientation of the wing in this figure, the vertical component of the induced flow w_i right behind the wing is upward.

2 Working principles and specifications of multi-rotor systems with lifting devices

2.1 Working principles of MRSLs

As mentioned in the previous section, one of the key reasons that conventional wind farms suffer from slow wake recovery rates is the absence of vertical flow. Regenerative wind farms counter this shortcoming by introducing vertical advection through the placement of lifting devices onto MRSLs. How this concept works is depicted by the vertical velocity fields w inside regenerative wind farms presented in Fig. 2, where the active exchange of flow between the upper and lower layers can be observed. This concept is inspired by the flow field induced by a wing described by the classic lifting-line theory (Anderson, 2011). As depicted in Fig. 3, the vorticity/circulation system of a wing can be simplified as a horseshoe vortex. The horseshoe vortex consists of two trailing vortices and a bound vortex. Due to the induction field of this vortex system, particularly from its trailing vortices, the induced flow u_i behind the wing has a non-zero vertical component (perpendicular to both the freestream and the spanwise directions), resulting in w_i≠0. Additionally, both the strength and the direction of w_i are affected by the wing's configuration. The strength of w_i is governed by the lift per span of the wing, with a higher lift generating a stronger circulation Γ and thus a larger w_i, as explained by the Kutta–Joukowski theorem and Helmholtz's theorem (Anderson, 2011). The direction of w_i can be altered by flipping the wing, that is, swapping the locations of the pressure side and the suction side. Moreover, stacking multiple wings vertically can further amplify w_i. Thus, in this work, MRSLs are equipped with several wings, referred to as the lifting devices, to increase the magnitude of w_i. By arranging these lifting devices as shown in Fig. 1, the flows at different altitudes behind MRSLs are exchanged vertically due to the non-zero w_i. It is this non-zero vertical flow induced by the lifting devices that fundamentally changes the mechanism of vertical energy entrainment within regenerative wind farms (Ferreira et al., 2024).

Based on the configurations of the lifting devices/wings of MRSLs, the lift exerted by MRSLs can be both upward and downward (note that the lift forces experienced by the flow and MRSL are in opposite directions). In this work, the configuration that exerts upward lift onto the flow is termed up-washing (UW), while the one that exerts downward lift is termed down-washing (DW). With the contours of vertical velocity w together with the isosurfaces of the streamwise vorticity magnitudes $\left|{\mathit{\omega }}_x\right|$ in Fig. 2, it can be seen that the flow at lower altitudes is channeled upward, while the flow at higher altitudes is brought downward for both UW and DW, enhancing the vertical exchange process. Note that the isosurfaces of ω_x represent the trailing vortices. The working principle of the lifting devices here is akin to the vortex generators on the wings of modern aircraft and the blades of contemporary wind turbines but on the scale of wind farms, which is much larger (Ferreira et al., 2024).

2.2 Specifications of MRSLs

In this work, the shape of the frontal area of an MRSL is set as a square (as shown in Fig. 1) with a side length D of 300 m, where the height of the rotor center z_rc is 186 m, corresponding to a clearance of 36 m. The lifting devices of MRSLs consist of four straight wings without any twist. These wings are placed at 100 %, 75 %, 50 %, and 25 % of the MRSL's height, as depicted in Figs. 1 and 4. Table 1 lists the key parameters of the MRSL used in this work. Note that the MRSL in Fig. 4 degenerates from Fig. 1, where the sub-rotors are represented with a single actuator disk (blue surface), and the lifting devices/wings are represented with four actuator lines (red surfaces). This simplification enables more efficient numerical modeling (Sorensen and Shen, 2002; Mikkelsen, 2004), and a detailed parameterization is provided in Sect. 3.4.

Figure 4 Illustration of how an MRSL is modeled with an actuator disk (blue) and four actuator lines (red) in the simulation. Note that the sketch showed here is a degenerated form of the MRSL presented in Fig. 1, and this simplification is aimed at facilitating the computational efficiency.

(Selig et al., 1995)

Table 1 Specifications of an MRSL modeled in the current work. D, c, and z_rc are the MRSL's side length, chord length, and height of the rotor center, respectively. The designed C_P (power coefficient) is estimated based on classic actuator disk theory (Manwell et al., 2010) with C_T=0.7. The designed T^R and P^R are the designed thrust and power of an entire MRSL estimated based on C_T=0.7 and C_P=0.54 with u_ref=10 m s⁻¹.

The thrust force exerted by an MRSL is calculated based on the sampled local velocity, where the thrust coefficient C_T is set to 0.7. According to classic actuator disk theory (Manwell et al., 2010), C_T=0.7 gives a power coefficient C_P of 0.54 (see Sect. 3.4 for more explanations). Note that C_P=0.54 is around the design values for modern large-scale wind turbines (Bak et al., 2013; Gaertner et al., 2020). The lifting devices of an MRSL consist of four straight wings with a constant profile, a constant twist angle, a constant chord length c, and a span of D. The chord length of the wings is set to $c=D/\mathrm{8}$, and the airfoil data used are the S1223 airfoil (Selig et al., 1995). The S1223 airfoil is chosen as it is one of the most representative profiles capable of achieving a high lift coefficient (Selig and Guglielmo, 1997). Additionally, the camber and thickness of the S1223 are relatively moderate, potentially making it more practical for real-world implementation. However, it is important to note that the specific airfoil profile is not critical to the implementation of the MRSL. The purpose of the MRSL's wings is to generate strong trailing vortices that enhance wake mixing and thus facilitate wake recovery. The airfoil coordinate and the lift–drag polar (calculated with the chord-based Reynolds number Re_c being 2×10⁷ using XFOIL version 6.99 (Drela, 1989)) for the S1223 airfoil are plotted in Fig. 5.

Figure 5 (a) X−Y plot of the cross-section of airfoil S1223 (Selig et al., 1995). (b) The lift–drag polar of airfoil S1223 obtained with XFOIL version 6.99 (Drela, 1989) with ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{7}}$.

3 Methodology

3.1 Numerical setup and computational domain

Numerical simulations of this work are conducted with OpenFOAM v2106 (OpenCFD Ltd., 2021), an open-source finite-volume-based CFD (computation fluid dynamics) solver. The flow is treated as incompressible and Newtonian (ρ=1.225 kg m⁻³ and $\mathit{\nu }=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹), and neither thermal effects nor Coriolis force is considered. The Reynolds-averaged Navier–Stokes (RANS) approach is employed. While higher-fidelity models such as large-eddy simulation (LES) are available for wind energy applications, RANS is selected for its lower computational demands (Thé and Yu, 2017), making it more suitable for rapid testing of new concepts and allowing for a broader parametric study. For the turbulence closure, the k–ω SST model (Menter, 1994) is chosen as it is the most widely used turbulence model in wind energy applications (Thé and Yu, 2017). A brief overview of the k–ω SST model is given in Appendix A, where the key governing equations are written. Additionally, a sensitivity test on the turbulence model is conducted in Appendix B, showing that the choice of turbulence model has a limited impact on the conclusions drawn from this work.

The spatial discretization schemes used are linear upwind (Gauss linearUpwind) for divergence and second-order central differencing (Gauss linear with limiter) for gradient and Laplacian. The pressure–velocity system is solved using the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm.

Figure 6 Diagram depicting the computational domain and the layout of a regenerative wind farm. The inflow comes from the bottom left to the top right. The wind farm consists of 15 multi-rotor systems with lifting devices (MRSLs), having a layout of five rows and three columns. The deep-blue surfaces represent the rotor part of MRSLs (thrusting devices), while the deep-red surfaces indicate the lifting devices of MRSLs. The semi-transparent volumes annotated with letters are the control volumes used in Sect. 4.5.

3.2 Computational domain and boundary conditions

The computational domain for the simulations is illustrated in Fig. 6. A Cartesian coordinate system is employed, with positive x pointing downstream and positive z pointing upward. The mesh is generated using the blockMesh application, consisting of uniformly sized cubic cells with a grid size of $\mathrm{\Delta }=D/\mathrm{25}$ in all three directions. The dimensions of the computational domain are $\mathrm{42}D\times \mathrm{21}D\times \mathrm{10}D$ in the x, y, and z (streamwise, lateral, vertical) directions, respectively, comprising approximately 137.8 M cells. Additionally, a grid independence test is carried out in Appendix C, confirming that a grid size of $\mathrm{\Delta }=D/\mathrm{25}$ is adequate for this study.

Table 2 Boundary conditions used for the simulation cases of regenerative wind farms immersed in ABL. atmBoundaryLayer and WallFunction are abbreviated as ABL and WF. For instance, ABLInletVelocity stands for atmBoundaryLayerInletVelocity.

A built-in library of OpenFOAM v2106, atmosphericModels (Richards and Hoxey, 1993; Hargreaves and Wright, 2007), is used to model the atmospheric boundary layer (ABL). The inlet profiles for the (mean) streamwise velocity u and turbulence kinetic energy k are given in Eqs. (1) and (2), respectively (see Fig. 12 for the generated freestream profile). z₀ is the surface roughness length, which is set to 10⁻⁴ m, a typical value for an offshore environment (Manwell et al., 2010). u_ref is the reference velocity at the height of the rotor center z_rc, which is set to 10 m s⁻¹. C₁ and C₂ are the two coefficients that are set to 0.814 and 1.0 in order to make the turbulence intensity TI =8 % at z=z_rc, where TI is defined as $\sqrt{\mathrm{2}k/\mathrm{3}}/\left|{\mathit{u}}_{\mathrm{ref}}\right|$. This value corresponds to the turbulence intensity observed in typical offshore environments (Hansen et al., 2012). To make this work more comprehensive, a sensitivity study on inflow turbulence intensity is performed in Appendix D, and it is demonstrated that variation in inflow TI does not overthrow the conclusion drawn later in this paper. The boundary conditions used in this work are listed in Table 2. For more detailed specifications of the boundary conditions used, readers are referred to the OpenFOAM v2106 documentation (OpenCFD Ltd., 2021). Additionally, note that any constants and model coefficients that are not explicitly mentioned are set to their default values, e.g., κ=0.41 and C_μ=0.09.

$$\begin{array}{}\text{(1)}& {\displaystyle }u={\displaystyle \frac{u^{*}}{\mathit{\kappa }}}\ln \left({\displaystyle \frac{z+z_{\mathrm{0}}}{z_{\mathrm{0}}}}\right),u^{*}={\displaystyle \frac{u_{\mathrm{ref}}\mathit{\kappa }}{\ln \left(\frac{z_{\mathrm{rc}}+z_{\mathrm{0}}}{z_{\mathrm{0}}}\right)}}\text{(2)}& {\displaystyle }k={\displaystyle \frac{{\left(u^{*}\right)}^{\mathrm{2}}}{\sqrt{C_{\mathit{\mu }}}}}\sqrt{C_{\mathrm{1}}\ln \left({\displaystyle \frac{z+z_{\mathrm{0}}}{z_{\mathrm{0}}}}\right)+C_{\mathrm{2}}}\end{array}$$

3.3 Wind farm layout

All the simulations in this work share the same wind farm layout, which consists of five rows and three columns. The MRSLs in each column are fully aligned with the direction of the freestream. The mid-column of the wind farm is placed at the centerline of the computational domain, and the first row is located 6 D from the inlet. The lateral distance between any two columns is 5 D, and the streamwise distance between the rows is 6 D. The origin of the coordinate system is set at the first row of the mid-column, as indicated in Fig. 6.

3.4 Modeling multi-rotor systems with lifting devices

The multi-rotor systems with lifting devices (MRSLs) introduced previously are parameterized using a square actuator disk (called disk for historical reasons) together with four actuator lines, as mentioned in Sect. 2. With actuator techniques, the effects of MRSL geometry are replaced by body force fields (term f_body in Eq. A2). This allows avoiding the exceptionally high computational cost required to resolve the boundary layer around the complex geometry (Sorensen and Shen, 2002; Mikkelsen, 2004). These actuator methods are realized in OpenFOAM using a customized library building upon actuationDiskSource (a built-in library of OpenFOAM v2106) and turbinesFoam (Bachant et al., 2019), and we call it flyingActuationDiskSource (Li et al., 2024b).

It should be noted that while the actuator techniques enable efficient simulations, they under-represent certain aerodynamic effects of the MRSL. For instance, the supporting structures of the MRSL are not modeled, the rotational effects of the sub-rotors are not captured, and the inter-spacings between its sub-rotors are not accounted for. Nevertheless, this work focuses on demonstrating the proof of concept for the aerodynamic capabilities of MRSLs within a regenerative wind farm, and detailed investigations of these secondary aerodynamic effects are left for future work. That said, it is worth noting that Broertjes et al. (2024) have already conducted experiments with an isolated MRSL equipped with rotating sub-rotors and supporting frames, demonstrating that the effectiveness of the lifting devices/wings is not significantly affected by these secondary effects, suggesting omitting them has limited impacts.

All the sub-rotors of an MRSL (thrusting devices) are modeled as a single actuator disk, and the actuator disk has 25 by 25 actuator elements that are situated on the same streamwise plane. These actuator elements have the same inter-distance in both lateral and vertical directions. The rotors of each MRSL have non-uniform loading based on the velocities sampled at each actuator element. T^ele and $C_{\mathrm{T}}^{\mathrm{ele}}$ in Eq. (3) denote the thrust force exerted by the actuator element and the corresponding thrust coefficient, respectively. $u_{\mathrm{in}}^{\mathrm{ele}}$ is the undisturbed inflow velocity seen by the actuator element. For all simulation cases in this work, the element-based area A^ele is $D^{\mathrm{2}}/\mathrm{625}$. The $C_{\mathrm{T}}^{\mathrm{ele}}$ targeted for each element for all MRSLs is set to 0.70. However, because the undisturbed inflow velocity perceived by an actuator element ($u_{\mathrm{in}}^{\mathrm{ele}}$) can vary when simulating wind farms and there is no universal method to define where to measure $u_{\mathrm{in}}^{\mathrm{ele}}$, estimating the value of T^ele for an actuator element directly based on $C_{\mathrm{T}}^{\mathrm{ele}}$ using Eq. (3) is challenging. To overcome this challenge, T^ele of this work is estimated based on the locally sampled velocity $u_{\mathrm{ls}}^{\mathrm{ele}}$ and the corrected thrust coefficient $C_{\mathrm{T}}^{*,\mathrm{ele}}$, as expressed in Eq. (4). Note that $u_{\mathrm{ls}}^{\mathrm{ele}}$ is the velocity sampled exactly where the actuator element is situated. Unlike $u_{\mathrm{in}}^{\mathrm{ele}}$, the sampling position of $u_{\mathrm{ls}}^{\mathrm{ele}}$ does not have ambiguity. $C_{\mathrm{T}}^{*,\mathrm{ele}}$ and $C_{\mathrm{T}}^{\mathrm{ele}}$ are linked through the classic actuator disk (one-dimensional momentum) theory (Manwell et al., 2010), which states $C_{\mathrm{T}}^{\mathrm{ele}}$ can be expressed as Eq. (5) based on the axial induction factor a^ele. After dividing/rearranging Eqs. (3) and (4) and applying the classic actuator disk theory (Eq. 5), the expression of $C_{\mathrm{T}}^{*,\mathrm{ele}}$ is obtained with $C_{\mathrm{T}}^{\mathrm{ele}}$ and a^ele, as written in Eq. (6). This method had been successfully implemented by Calaf et al. (2010). Through Eq. (5), it can be calculated that $C_{\mathrm{T}}^{\mathrm{ele}}=\mathrm{0.7}$ infers a^ele=0.23, which leads to $C_{\mathrm{T}}^{*,\mathrm{ele}}=\mathrm{1.17}$.

$$\begin{array}{}\text{(3)}& {\displaystyle }{T}^{\mathrm{ele}}=\mathrm{0.5}\mathit{\rho }(u_{\mathrm{in}}^{\mathrm{ele}}{)}^{\mathrm{2}}{A}^{\mathrm{ele}}{C}_{\mathrm{T}}^{\mathrm{ele}}\text{(4)}& {\displaystyle }{T}^{\mathrm{ele}}=\mathrm{0.5}\mathit{\rho }(u_{\mathrm{ls}}^{\mathrm{ele}}{)}^{\mathrm{2}}{A}^{\mathrm{ele}}{C}_{\mathrm{T}}^{*,\mathrm{ele}}\text{(5)}& {\displaystyle }{C}_{\mathrm{T}}^{\mathrm{ele}}\simeq \mathrm{4}{a}^{\mathrm{ele}}(\mathrm{1}-a^{\mathrm{ele}}),a^{\mathrm{ele}}\stackrel{\Delta }{=}\mathrm{1}-{\displaystyle \frac{u_{\mathrm{ls}}^{\mathrm{ele}}}{u_{\mathrm{in}}^{\mathrm{ele}}}}\text{(6)}& {\displaystyle }{C}_{\mathrm{T}}^{*,\mathrm{ele}}=C_{\mathrm{T}}^{\mathrm{ele}}{\left({\displaystyle \frac{u_{\mathrm{in}}^{\mathrm{ele}}}{u_{\mathrm{ls}}^{\mathrm{ele}}}}\right)}^{\mathrm{2}}\simeq {\displaystyle \frac{C_{\mathrm{T}}^{\mathrm{ele}}}{(\mathrm{1}-a^{\mathrm{ele}}{)}^{\mathrm{2}}}}\end{array}$$

After obtaining the value of T^ele through Eq. (4), the force is projected onto the CFD grid with Eq. (7), where f^ele is the force vector exerted by the actuator element, and ${\mathit{f}}_{\mathrm{body}}^{\mathrm{ele}}\left(\mathit{x}\right)$ is the body force field on the CFD grid projected by f^ele at position x. ξ^ele denotes the position vector of the actuator element. The projection is done by the Gaussian normalization kernel, and it is introduced to improve the robustness of the numerical modeling (Sorensen and Shen, 2002; Mikkelsen, 2004), where ε is called the smearing factor. For the actuator elements of the MRSL's rotors, ${\mathit{f}}^{\mathrm{ele}}=-T^{\mathrm{ele}}{\widehat{\mathit{e}}}_x$ is assigned, and its smearing factor, denoted as ε^R, is set to 1.0Δ, as it is commonly used for actuator disks (Mikkelsen, 2004; Wu and Porté-Agel, 2011). The thrust and power of the rotor (T^R and P^R) are calculated after projecting the body force fields using Eq. (8). Here, i and j represent the indices for positions and actuator elements, respectively. Note that Δ³ is the volume of a cell.

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}& {\mathit{f}}_{\mathrm{body}}^{\mathrm{ele}}\left(\mathit{x}\right)={\mathit{f}}^{\mathrm{ele}}{\mathit{\eta }}_{\mathit{\epsilon }}(\parallel \mathit{x}-{\mathit{\xi }}^{\mathrm{ele}})\parallel )\\ & {\mathit{\eta }}_{\mathit{\epsilon }}\left(d\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathrm{3}}{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}}}\exp \left[-{\left({\displaystyle \frac{d}{\mathit{\epsilon }}}\right)}^{\mathrm{2}}\right]\end{array}\text{(8)}& {\displaystyle }\begin{array}{rl}& T^R=\sum _i\sum _j{f}_{j,\mathrm{body}}^{\mathrm{ele}}\left({\mathit{x}}_i\right){\mathrm{\Delta }}^{\mathrm{3}}\\ & P^R=-\sum _i\sum _{j}u\left({\mathit{x}}_i\right)f_{j,\mathrm{body}}^{\mathrm{ele}}\left({\mathit{x}}_i\right){\mathrm{\Delta }}^{\mathrm{3}}\end{array}\end{array}$$

As mentioned in Sect. 2, the lifting devices of MRSLs are parameterized with four actuator lines, with each having 25 equally spaced actuator elements lining up in the lateral direction, and these actuator elements are in the same plane as those of the rotors. The forces to project are calculated based on the blade element approach, where f^AL is calculated based on the velocity sampled and the airfoil polar, as written in Eq. (9). u^AL is the sampled flow velocity for an actuator element of an actuator line. ${\mathit{f}}_l^{\mathrm{AL}}$, ${\mathit{f}}_d^{\mathrm{AL}}$, C_l, and C_d are the lift and drag forces and their corresponding coefficients. In this work, C_l and C_d are based on the polar data of the S1223 airfoil plotted in Fig. 5. Δ^AL is the span length to which the actuator element corresponds. In this work, ${\mathrm{\Delta }}^{\mathrm{AL}}=D/\mathrm{25}$. ${\widehat{\mathit{e}}}_s$, ${\widehat{\mathit{e}}}_l$, and ${\widehat{\mathit{e}}}_d$ are the unit vectors in the spanwise, lift, and drag directions, respectively. Note that ${\widehat{\mathit{e}}}_s\Vert \pm {\widehat{\mathit{e}}}_y$ (depending on the lifting direction) and ${\widehat{\mathit{e}}}_l\Vert ({\mathit{u}}^{\mathrm{AL}}\times {\widehat{\mathit{e}}}_s)$. In this work, u^AL is obtained by averaging the 20 velocity samples sampled on a circular path with the actuator element at the center (line averaging). The sampling points are equidistant, and the normal direction of the enclosed surface is parallel to the spanwise direction. The radius of the circle is set to $r^{\mathrm{AL}}=\mathrm{3}\mathrm{\Delta }\simeq c$. Single-point sampling is avoided to achieve better robustness (Melani et al., 2021). Note that since ${\mathit{u}}^{\mathrm{AL}}⟂{\widehat{\mathit{e}}}_l$ and the wings are stationary, the lift forces of the wings do not do any work on the flow.

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathit{f}}^{\mathrm{AL}}& =\left({\mathit{f}}_l^{\mathrm{AL}},{\mathit{f}}_d^{\mathrm{AL}}\right)\\ & =\mathrm{0.5}\mathit{\rho }{\left(u^{\mathrm{AL}}\right)}^{\mathrm{2}}c{\mathrm{\Delta }}^{\mathrm{AL}}\left(C_l\right(\mathit{\alpha }){\widehat{\mathit{e}}}_l,C_d(\mathit{\alpha }\left){\widehat{\mathit{e}}}_d\right)\\ & =f_z^{\mathrm{AL}}{\widehat{\mathit{e}}}_z+f_x^{\mathrm{AL}}{\widehat{\mathit{e}}}_x\end{array}\end{array}$$

A Gaussian normalization kernel (see Eq. 7, where f^ele is replaced with f^AL) is used again to project the forces of the lifting devices onto the CFD grid. For the smearing factor ε, instead of assigning a single value, the values of ε for the actuator lines (denoted as ε^W) are calculated based on the relative wing position as described in Eq. (10) ($r/D=\mathrm{0.0}$ corresponds to the middle of the wing). This approach was introduced by Jha et al. (2014). Compared with the experimental results of the load of a finite wing, it has been shown that using this distribution of ε^W outperformed the case using a single value for ε^W (Jha et al., 2014; Jha and Schmitz, 2018). For the current work, n_max is assigned as 3.0, and the distribution of ε^W along the wing used in this work is plotted in Fig. 7.

$$\begin{array}{}\text{(10)}& {\mathit{\epsilon }}^W=n_{\max }\mathrm{\Delta }\sqrt{\mathrm{1}-{\left({\displaystyle \frac{\mathrm{2}r}{D}}\right)}^{\mathrm{2}}},-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\le {\displaystyle \frac{r}{D}}\le {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\end{array}$$

Using a similar method to obtain T^R and P^R (Eq. 8), the total lift L^W (vertical force) and the total induced drag $D_{\mathrm{ind}}^W$ (streamwise force) of the four wings of an MRSL are obtained through Eq. (11). It should be noted that the directions of L^W and ${\mathit{D}}_{\mathrm{ind}}^W$ are based on the global coordinate system, which is different from the directions of ${\mathit{f}}_l^{\mathrm{AL}}$ and ${\mathit{f}}_d^{\mathrm{AL}}$, which are based on the local flow direction seen by the airfoil (u^AL).

Figure 7 The distribution of ε^W across the wings used in this work. $\mathrm{\Delta }=D/\mathrm{25}$ is the grid size of the mesh used in the cases listed in Table 3.

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& L^W=\sum _i\sum _j{\mathit{f}}_{j,\mathrm{body}}^{\mathrm{AL}}\left({\mathit{x}}_{\mathit{i}}\right)\cdot {\widehat{\mathit{e}}}_z{\mathrm{\Delta }}^{\mathrm{3}}\\ & D_{\mathrm{ind}}^W=\sum _i\sum _j{\mathit{f}}_{j,\mathrm{body}}^{\mathrm{AL}}\left({\mathit{x}}_i\right)\cdot {\widehat{\mathit{e}}}_x{\mathrm{\Delta }}^{\mathrm{3}}\end{array}\end{array}$$

To adjust the magnitudes of the lift force exerted by the lifting devices, the wings of MRSLs are pitched in the simulations by varying their pitch angles θ_p. Specifically, θ_p of each wing is adjusted so that the angle of attack α at the midpoint of the wing corresponds to a specified lift coefficient C_l. That is, C_l,mid, the lift coefficient at the midpoint of a wing, is tuned to a specified value. Since the inflow conditions for each of the MRSLs' wings differ, θ_p varies for each wing. The adjustment of θ_p is programmed and carried out automatically during the simulations. Note that each wing is pitched as a whole and has a constant twist angle along its entire span. For a demonstration, see Appendix E, where profiles of α along the wings are presented. In this work, for the MRSLs in the cases equipped with lifting devices, all their wings are pitched to make $C_{l,\mathrm{mid}}=\mathrm{2.5}$, except for the cases in Appendix F.

3.5 Test matrix

The main context of this study includes three simulations. The three cases are without lifting (WL), up-washing (UW), and down-washing (DW), as listed in Table 3. In the WL case, MRSLs are not equipped with lifting devices, serving as the reference case. In the UW case, the lifting devices on MRSLs exert upward vertical force onto the flow, and one of the immediate effects is that the wakes right behind the MRSLs are directed upward. Contrary, in the DW case, the lifting devices are designed to exert downward vertical force (the orientation of the wings is flipped compared to the UW category), sending the wakes right behind the MRSLs downward.

For the lifting devices of the MRSLs in the UW and DW cases, their wings are pitched during the simulations to make C_l,mid for each wing equal to 2.5 (C_l,mid is the lift coefficient at the mid-span of a wing; see the end of Sect. 3.4). Note that based on some rough estimations using the specifications provided in Table 1, $C_{l,\mathrm{mid}}=\mathrm{2.5}$ allows an MRSL to generate a vertical force that is of a similar magnitude to the thrust force of its rotors.

Table 3 Test matrix of the simulation cases considered in the main context.

4 Results and discussions

4.1 Forces exerted by MRSLs

The thrust of the MRSL's rotors together with the lift (vertical force) and the induced drag (streamwise force) of the MRSL's lifting devices/wings are plotted in Fig. 8. The three cases listed in Table 3 are displayed. For the loading profiles of the MRSL's wings, see Appendix E.

${\widehat{T}}^R$, ${\widehat{L}}^W$, and ${\widehat{D}}_{\mathrm{ind}}^W$ in Eq. (12) and Fig. 8 are the normalized thrust T^R, lift L^W, and induced drag $D_{\mathrm{ind}}^W$ of the MRSL (Eqs. 8 and 11), respectively. These forces are normalized against T^R measured at the first row in the mid-column of the WL case in Table 3, denoted as $T^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$, which is 3.87 MN. This value is very close to the designed value of 3.86 MN (based on letting C_T=0.7 and a reference velocity of u_ref=10 m s⁻¹), validating the actuator model described in Sect. 3.4.

$$\begin{array}{}\text{(12)}& {\widehat{T}}^R\stackrel{\Delta }{=}{\displaystyle \frac{T^R}{T^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}}},{\widehat{L}}^W\stackrel{\Delta }{=}{\displaystyle \frac{\left|L^W\right|}{T^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}}},{\widehat{D}}_{\mathrm{ind}}^W\stackrel{\Delta }{=}{\displaystyle \frac{\left|D_{\mathrm{ind}}^W\right|}{T^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}}}\end{array}$$

Operator 〈⋅〉 in this work indicates row averaging. For example, $\langle {\widehat{T}}^R\rangle$ in Fig. 8 denotes the row-averaged normalized rotor thrust. The results show that the value differences between the middle and side columns are at most 1 % for T^R, L^W, and $D_{\mathrm{ind}}^W$ for the three cases in Table 3.

As shown in the left and middle panels of Fig. 8, as designed, $\langle {\widehat{L}}^W\rangle$ for the MRSLs in the two cases with lifting devices is similar to $\langle {\widehat{T}}^R\rangle$, while the WL case has zero lift. Additionally, for both UW and DW, it can be observed that their $\langle {\widehat{T}}^R\rangle$ is much higher than that of WL from the second row onward, despite the fact that the lifting devices also introduce significant $\langle {\widehat{D}}_{\mathrm{ind}}^W\rangle$, as shown in the right panel of Fig. 8. Specifically, it is found that the thrust for the two cases with lifting devices decreases only slightly from the first to the second row, with $\langle {\widehat{T}}^R\rangle$ remaining above 80 %, and the decreasing trend ceases from the third row onward. In contrast, for the case without lifting devices, $\langle {\widehat{T}}^R\rangle$ drops significantly from the first to the second row, falling below 60 %, and continues to decrease row by row. By the third row, $\langle {\widehat{T}}^R\rangle$ for the two cases with lifting devices is more than double compared to the WL case. Additionally, the fact that the forces for the MRSLs in the UW and DW cases remain relatively stable from the third to the fifth row suggests that these values would likely be sustainable if the regenerative wind farms had more rows. Furthermore, higher values of $\langle {\widehat{T}}^R\rangle$ suggest that the streamwise velocity experienced by an MRSL at a given row is much higher when lifting devices are equipped. This is further confirmed by the plots and contours presented in later sections (Sect. 4.3).

Figure 8 The normalized row-averaged thrust of the MRSL's rotor (${\widehat{T}}^R$, a) together with the vertical (${\widehat{L}}^W$, b) and streamwise (${\widehat{D}}_{\mathrm{ind}}^W$, c) force components of the MRSL's lifting devices. The normalization is done by dividing the reference rotor thrust, which is based on the MRSL at the first row in the mid-column of the WL case. The legends correspond to the case names introduced in Table 3.

4.2 Power harvested by MRSLs

Figure 9 presents the normalized row-averaged power $\langle {\widehat{P}}^R\rangle$ harvested by the rotors of MRSLs for the three cases listed in Table 3. These values are plotted alongside those predicted by the Frandsen wake model (Frandsen et al., 2006) (see Appendix G). As in the previous subsection, the rotor power P^R is normalized based on the MRSL located at the first row in the mid-column of the WL case. The reference power, denoted as $P^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$, is 30.1 MW. This value corresponds to a power density of 11.1 W m⁻². This power density is calculated by dividing $P^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$ by the footprint area of an MRSL, which in this study is 6 D×5 D.

Figure 9 The normalized row-averaged rotor power of the MRSL's rotor (${\widehat{P}}^R$). The normalization is done by dividing the values by the rotor power of the MRSL situated at the first row in the mid-column of the WL case. Calculation of the use of the Frandsen wake model is detailed in Appendix G.

Very good agreement was found between the CFD results for the without-lifting (WL) case and the predictions of the Frandsen wake model, which supports the validity of the numerical framework used in this work.

As expected, ${\widehat{P}}^R$ of the three cases is highly correlated with ${\widehat{T}}^R$ (as indicated by Eq. 8), with the cases having lifting devices also exhibiting higher ${\widehat{P}}^R$. However, in terms of the magnitudes, the relative differences in ${\widehat{P}}^R$ between cases with and without lifting devices are greater than those in ${\widehat{T}}^R$, since ${\widehat{P}}^R$ is proportional to the cube of the sampled velocity, while ${\widehat{T}}^R$ is proportional to the square of it (Eqs. 4 and 8).

Examining the values of $\langle {\widehat{P}}^R\rangle$ for the first row of the three cases in Fig. 9, it is observed that $\langle {\widehat{P}}^R\rangle$ for the UW case is higher than that of the WL case. In contrast, the DW case exhibits the opposite behavior. This can be attributed to the wings (lifting devices) of the MRSLs acting as diffuser-like devices. A straightforward explanation is that the bound circulations of the wings (Anderson, 2011) either accelerate or decelerate the flow velocity crossing the rotor (thrusting devices) of an MRSL, depending on the configuration of the lifting devices. Previous studies have reported similar phenomena with comparable configurations (Bader et al., 2018). Although this effect influences the power output of MRSLs, it is overshadowed by the effects of the enhanced wake recoveries due to the lifting devices. Therefore, it is not discussed or quantified in the rest of this work.

When comparing the power output row by row across the entire regenerative wind farm, it is found that the cases with lifting devices have significantly higher values for $\langle {\widehat{P}}^R\rangle$ compared to the case without at and after the second row. Specifically, $\langle {\widehat{P}}^R\rangle$ for the UW and DW cases at the second row is more than double that of the WL case. Remarkably, for the third to the fifth rows, $\langle {\widehat{P}}^R\rangle$ for the two cases with lifting devices is more than triple compared to that without. Furthermore, despite the relatively small spacing (around 5.3 D_cir, considering the shape effects of the rotor; see Appendix G), $\langle {\widehat{P}}^R\rangle$ for the cases with lifting devices remains at least 80 % of the reference power up to the fifth row. This significantly outperforms conventional wind turbines (i.e., HAWT), which typically maintain around 40 % to 60 % when the inflow is aligned with the wind farm layout and when the streamwise spacing is 5 D_cir or 7 D_cir, respectively (Barthelmie et al., 2010; Li et al., 2024a; Wu and Porté-Agel, 2015). These power output results underscore the profound potential of the concept of regenerative wind farms, supporting the current proposal.

The overall performance of the regenerative wind farms is evaluated based on power density, which serves as a measure of the efficiency of the regenerative wind farms. Table 4 lists the relative power densities of the regenerative wind farms, with 100 % corresponding to 11.1 W m⁻², which is the power density of the MRSL at the first row in the mid-column in the WL case.

Similarly, as has been seen in the plot of $\langle {\widehat{P}}^R\rangle$ (Fig. 9), the result of the WL case has very good agreement with the prediction given by the Frandsen wake model in Table 4. By comparing the other values, it is evident that the two cases with lifting devices (cases UW and DW) have power densities that are nearly double the power density of the WL case, increasing from approximately 40 % to around 80 %. In other words, the power losses due to wake interactions among the regenerative wind farms are reduced from roughly 60 % to about 20 % by introducing lifting devices. These results demonstrate the capabilities of MRSLs and the tremendous potential of regenerative wind farms in achieving significantly higher wind farm efficiencies than conventional wind farms.

Table 4 The relative power densities of the regenerative wind farms of the three cases in Table 3 and the values predicted by the Frandsen wake model (see Appendix G). Here, 100 % corresponds to 11.1 W m⁻², which is the power density of the MRSL at the first row in the mid-column in the WL case.

4.3 Flow field characterization

4.3.1 Three-dimensional flow structures

Figure 10 illustrates the three-dimensional flow structures of the simulated wind farms based on streamwise velocity. All three cases in Table 3 are depicted. The plots cover the mid-column of the regenerative wind farms, with the positions of the MRSLs represented by deep-blue surfaces for the rotors and deep-red surfaces for the lifting devices/wings. The low-speed wakes are depicted by light-blue isosurfaces, corresponding to $u/u_{\mathrm{ref}}=\mathrm{0.65}$. Additionally, several x planes color-coded by streamwise velocity u are displayed, and the directions of in-plane velocity are indicated by arrows.

In the plot for the WL case, it is evident that the MRSLs after the second row are generally immersed in the wakes of the upstream ones, resulting in significantly lower inflow velocities compared to the first row. Additionally, based on the arrows in the plot, it can be seen that vertical velocity is generally absent, making its wake recovery rates slow. Consequently, as shown in Fig. 9, the power outputs of the MRSLs after the second row are much lower compared to those in the first row for the WL case.

Figure 10 Three-dimensional flow structures of the regenerative wind farms around their mid-columns. The without-lifting, up-washing, and down-washing cases are plotted at the top, middle, and bottom, respectively. MRSLs are represented by surfaces in deep blue and deep red, which indicate their rotors and wings. The isosurfaces in light blue depict the wakes of the MRSLs, corresponding to where $u/u_{\mathrm{ref}}=\mathrm{0.65}$. Additionally, sections with contours of streamwise velocity in the x planes are plotted, with arrows indicating the directions of the in-plane velocity. Note that the arrows' lengths are scaled by the square root of the in-plane velocity's norm. The frontal projections of the MRSLs are illustrated with light-green squares.

In the UW case, the wakes of the MRSLs are significantly steered upward, where the cores of the wakes (indicated by the light-blue surfaces) are mostly redirected away from the frontal areas of the MRSLs. This results in much higher P^R for the downstream MRSLs compared to the WL case (see Fig. 9). Furthermore, it is observed that the wakes' positions are further elevated as the flow progresses deeper into the regenerative wind farm, indicating that the effects of UW accumulate progressively across rows. Additionally, arrows on the slices of the velocity contour reveal pairs of counter-rotating vortices (CRVs) formed by the trailing vortices released by the lifting devices (these CRVs can be seen clearer in Fig. 15 with the arrows). These CRVs lift the exhaust wakes upward and spread them laterally, simultaneously bringing down fresh, clean flows from above, thereby replenishing the lower layers, where MRSLs are situated, with higher-energy flows. These CRVs enhance the vertical energy entrainment process by promoting mixing in the vertical direction. See Sect. 4.4 and 4.5 for further discussions on CRVs and the vertical energy entrainment process.

In the DW case, the wakes of the upstream MRSLs are also steered away from the frontal areas of the downstream MRSLs, reducing the wake losses experienced by the downstream units. However, the presence of the ground makes the dynamics of the DW case quite different from those of the UW case. In the DW scenario, the wakes are initially directed downward. Then, they are quickly forced to spread sideways as the ground prevents further downward penetration. As the wakes accumulate on the sides as they go deeper into the regenerative wind farm, they eventually start to move upward. Like in the UW case, CRVs are also present in the DW case but rotate in the opposite direction. In this configuration, the CRVs bring fresh, clean flows down from above at the centerlines of the MRSLs while steering the exhausted wakes down and sidewards.

It is important to note that the purpose of the lifting devices is not limited to steering the wakes vertically. In fact, the primary goal of the lifting devices is to introduce a vertical advection process that enhances vertical mixing, as stronger vertical mixing leads to stronger vertical energy entrainment. A key aspect of Fig. 10 is that the blueish areas in the streamwise velocity contours (areas where u<u_ref) for the two cases with lifting devices are significantly larger than those without. This indicates that the mixing process of the cases with lifting devices is more pronounced, and therefore the effectiveness of lifting devices is demonstrated. However, although the significant potential of lifting devices is presented, their effectiveness in regenerative wind farms with different layouts and sizes remains uncertain, necessitating further investigation in future studies.

Figure 11 Contours of streamwise velocity u of the regenerative wind farms with the without-lifting, up-washing, and down-washing cases in Table 3. The slices are cut at $y/D=\mathrm{0.0}$ and u_ref=10 m s⁻¹. The contours are superimposed with the streamlines based on the in-plane velocity (u and w). The thick black lines represent the positions of the MRSLs.

4.3.2 Streamlines

The contours of streamwise velocity u, superimposed with streamlines for the three cases in Table 3, are shown in Fig. 11. These contours are based on the data on the slices at $y/D=\mathrm{0.0}$, corresponding to the middle of the mid-column. In the WL case, no significant vertical flows are indicated by the streamlines, suggesting that the vertical advection process is generally absent. In contrast, for the cases with lifting devices (UW and DW), the streamlines show steep slopes right behind the MRSLs, indicating strong vertical advection and significant vertical mixing.

Additionally, it can be observed that the thickness of the wakes (the blueish area) in the WL case remains nearly constant after the second row (around 1.5 D). In the UW case, the wake thickness progressively increases as it moves deeper into the regenerative wind farm, growing from around 1.0 D to 3.5 D. On the other hand, in the DW case, the wake thickness decreases with each subsequent row of MRSLs, dropping from 1.0 D to 0.5 D. However, it should be noted that the surfaces in Fig. 11 are confined to $y/D=\mathrm{0.0}D$. If the surfaces were shifted along the y direction, it would be evident that the wakes in both the UW and the DW cases penetrate higher than those of the WL case (the thickness of wakes for DW can reach around 2.1 D), as can be confirmed with the contours of u displayed in Fig. 10. This again demonstrates that lifting devices enhance vertical mixing within regenerative wind farms.

4.3.3 Lateral-averaged streamwise velocity profiles

This subsection explores the lateral-averaged velocity profiles in the regenerative wind farms. Two lateral-averaging ranges are considered, which are $-\mathrm{0.5}\le y/D<\mathrm{0.5}$ and $-\mathrm{2.5}\le y/D<\mathrm{2.5}$, and their lateral-averaged velocities are denoted as 〈u〉_±0.5 D and 〈u〉_±2.5 D, respectively. Note that 〈u〉_±0.5 D averages over the frontal area of MRSLs situated in the mid-column, while 〈u〉_±2.5 D averages over the entire mid-column. Figure 12 presents the vertical profiles of 〈u〉_±0.5 D (left) and 〈u〉_±2.5 D (right) at $x/D=\mathrm{10.0}$, 16.0, and 22.0, which are located 4 D downstream from the second, third, and fourth rows of the regenerative wind farms. These positions are selected since they are completely within the regenerative wind farms, where the influences of the wind farm boundaries are relatively small. Additionally, the distance of 4 D is far enough from the upstream MRSLs, while the induction effects of the MRSLs in the next row are minimal.

In the plot of 〈u〉_±0.5 D profiles (Fig. 12a), it is evident that both case UW and case DW exhibit significantly larger values for 〈u〉_±0.5 D around the heights of the MRSLs compared to the WL case, and the differences are greater with larger $x/D$. This has already been reflected in the values of power output reported in Fig. 9. Additionally, the shapes of the velocity profiles differ significantly between the two lifting configurations. In the DW case, the profiles closely resemble the freestream profiles, suggesting that the flow's mean kinetic energy (MKE) is replenished. Upon closer inspection between $\mathrm{1.5}<z/D<\mathrm{3.0}$, the 〈u〉_±0.5 D profiles for the DW case are slightly higher than the freestream profiles. This is related to the strong downward vertical velocities around $y/D=\mathrm{0.0}$, which entrain higher streamwise velocity from the upper layers to the lower ones. In contrast, for the UW case, the 〈u〉_±0.5 D profiles decrease with z from z≃0.2 D to z≃1.5 D, which are atypical velocity profiles for standard atmospheric boundary layers. These shapes indicate that the wakes of the MRSLs are channeled upward in case UW and that the MKE entrainment is primarily from the sides of MRSLs at the lower layers.

Figure 12 Lateral-averaged streamwise velocity profiles. The velocity profiles are sampled at 4 D after the second, third, and fourth rows of the wind farm ($x/D=\mathrm{10.0}$, 16.0, and 22.0). The lateral averaging range for (a) is $-\mathrm{0.5}D\le y<\mathrm{0.5}D$, which covers the frontal area of the MRSL of the mid-column. For (b), the lateral averaging range is $-\mathrm{2.5}D\le y<\mathrm{2.5}D$, covering the entire mid-column. The freestream profiles are based on the results of the without-lifting case at $x/D=-\mathrm{2.0}$. Without lifting, up-washing, and down-washing correspond to the case names in Table 3.

For the profiles of 〈u〉_±2.5 D, case DW notably underperforms compared to case WL around the height of the MRSLs ($\mathrm{0.12}<z/D<\mathrm{1.12}$). This is mainly because the MRSLs of case DW have extracted more power from the flow at these positions, and the induced drag from the lifting devices also negatively impacts 〈u〉_±2.5 D. In contrast, case UW still significantly outperforms case WL around the height of the MRSLs, even though it also extracts more energy and introduces induced drag as case DW. This difference is because case UW ejects most of its exhausted wakes upward, while the wakes in case DW are mostly trapped at lower altitudes.

For both 〈u〉_±0.5 D and 〈u〉_±2.5 D in Fig. 12, it is evident that at higher altitudes (larger $z/D$), case UW has more pronounced effects on altering the velocity profiles compared to case DW. Moreover, for UW, it is found that the maximum deficits of both 〈u〉_±0.5 D and 〈u〉_±2.5 D appear to be higher for larger $x/D$. This observation aligns with the circulation-based analysis carried out in the later section (Sect. 4.4), where it is found that the positions of CRVs ($z_{{\mathrm{\Gamma }}_x}$) for case UW progressively rise as the flow moves deeper into the regenerative wind farm, while this is not the case for DW. Additionally, this observation further suggests that the UW configuration may have inherent advantages in enhancing vertical entrainment, as it can extend its effects to higher layers of the ABL compared to the DW configuration.

Figure 13 Contours of streamwise vorticity ω_x on x planes at different x positions. The x positions are indicated at the top of each column. These x positions are the same as those in Fig. 10, except $x/D=-\mathrm{2.0}$ is excluded. Cases WL, UW, and DW in Table 3 are plotted in the top, middle, and bottom rows, respectively. The frontal projections of the MRSLs are illustrated with light-green squares. The vortical structures enclosed by the dashed magenta lines are used to calculate the streamwise circulation-related quantities analyzed in Sect. 4.4. The arrows indicate the direction of the in-plane velocity. Their lengths are scaled by the square root of the in-plane velocity's norms, and their absolute scales are indicated at the bottom right of the figure.

4.3.4 Vorticity fields

The fields of streamwise vorticity ω_x on the selected x planes (the same as those in Fig. 10 except that $x/D=-\mathrm{2.0}$ is dropped) for cases WL, UW, and DW are plotted in Fig. 13. In this figure, large-scale vortical structures appear in the wakes of both the UW and the DW cases, which we call counter-rotating vortices (CRVs). It is these CRVs that promote the vertical energy entrainment process for cases UW and DW, and it is evident that CRVs are absent in the wake of case WL. Additionally, the plots show that CRVs increase in size as they progress deeper into the regenerative wind farms for the cases with lifting devices. As described earlier, these CRVs originate from the trailing vortices released by the lifting devices of the MRSLs (see Fig. 2 for three-dimensional representations of CRVs with isosurfaces of $\left|{\mathit{\omega }}_x\right|$). Thus, as the downstream rows release their trailing vortices, the existing CRVs are strengthened, as can be assessed qualitatively in Fig. 13. Furthermore, visual inspection reveals that the centers of the CRVs rise as the flow passes through more rows of MRSLs in case UW, while in case DW, the centers of the CRVs are observed to be pushed primarily sideways. These observations highlight that the dynamics of CRVs depend on the lifting configurations of the MRSLs, which is further discussed in Sect. 4.4.

Figure 14 Contour plots of lateral velocity v and vertical velocity w for the WL, UW, and DW cases in Table 3 on x planes at different x positions. The x positions are indicated at the top of each column. The two quantities are plotted side by side in each panel, where the left is v and the right is w. The arrows in the plots represent the in-plane velocities. The arrows indicate the direction of the in-plane velocity. Their lengths are scaled by the square root of the in-plane velocity's norms, and their absolute scales are indicated at the bottom right of the figure.

4.3.5 Lateral and vertical velocities

The contours of lateral and vertical velocities (v and w) at selected x planes (the same as those for ω_x in Fig. 13) are presented in Fig. 14. In each panel, the left side displays the v fields, while the right side shows the w fields. Additionally, to supplement these results, contours of w for UW and DW at $x/D=\mathrm{22.0}$ for the entire wind farm cross-section are provided in Fig. 2.

Contours of v and w are analyzed as they are directly linked to the advection processes that enhance energy entrainment. As shown in Fig. 13, the magnitudes of both v and w increase with $x/D$ for UW and DW, while they remain minimal for WL across all considered x planes. Notably, both v and w exceed 20 % of u_ref after the second row of the wind farm, demonstrating the strong influence of CRVs and significant advection processes.

Furthermore, the regions of high w values for both cases with lifting devices extend above the MRSLs' height, indicating enhanced vertical mixing, which leads to vertical energy entrainment. This effect is particularly pronounced in UW, where regions with $\left|w\right|$ greater than 5 % of u_ref reach up to $z/D=\mathrm{2.5}$.

4.4 Quantification of counter-rotating vortices

Utilizing circulation, this section assesses the counter-rotating vortices (CRVs) identified in Figs. 10 and 13 in a quantitative manner. Based on the fields of ω_x, streamwise circulations Γ_x of all three cases in Table 3 are calculated to represent the strength of the CRVs, which are presented in Fig. 15. The values of Γ_x in Fig. 15 are obtained using Eq. (13). Stokes' theorem is applied in Eq. (13), with C being the contour bounding the surface S.

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathrm{\Gamma }}_x& =\underset{C}{\oint }(u,v,w)\cdot \mathrm{d}\mathit{l}=\underset{S}{\int }\mathrm{\nabla }\times (u,v,w)\cdot \mathrm{d}\mathit{A}\\ & =\underset{S}{\int }({\mathit{\omega }}_x,{\mathit{\omega }}_y,{\mathit{\omega }}_z)\cdot \mathrm{d}\mathit{A}=\underset{S}{\int }{\mathit{\omega }}_x\mathrm{d}A\mathrm{for}\mathrm{d}\mathit{A}\Vert {\widehat{\mathit{e}}}_x\end{array}\end{array}$$

In this work, the strength of CRVs is defined by the magnitudes of Γ_x (denoted as $\left|{\mathrm{\Gamma }}_x\right|$) calculated using the rightmost of Eq. (13). Moreover, the positions of the CRVs are defined based on the center of gravity (CoG) of the vortical structure (Saffman, 1995), which is calculated through Eq. (14), where $z_{{\mathrm{\Gamma }}_x}$ and $y_{{\mathrm{\Gamma }}_x}$ are defined as the z and y positions of CRVs, respectively. Note that when calculating the Γ_x-related quantities in the current work, only the regions within $\mathrm{0.0}<y/D<\mathrm{2.5}$ and $\mathrm{0.0}<z/D<\mathrm{5.0}$ are considered. Furthermore, only the ω_x with the prevailing sign in that region is considered. Examples are illustrated in Fig. 13, where Γ_x as well as $z_{{\mathrm{\Gamma }}_x}$ and $y_{{\mathrm{\Gamma }}_x}$ are calculated based on the regions enclosed by the dashed magenta lines.

$$\begin{array}{}\text{(14)}& z_{{\mathrm{\Gamma }}_x}\stackrel{\Delta }{=}{\displaystyle \frac{{\int }_{S}z{\mathit{\omega }}_x\mathrm{d}A}{{\int }_S{\mathit{\omega }}_x\mathrm{d}A}},y_{{\mathrm{\Gamma }}_x}\stackrel{\Delta }{=}{\displaystyle \frac{{\int }_{S}y{\mathit{\omega }}_x\mathrm{d}A}{{\int }_S{\mathit{\omega }}_x\mathrm{d}A}}\end{array}$$

Figure 15 Plots of $\left|{\mathrm{\Gamma }}_x\right|$ (strength of CRVs, a), $z_{{\mathrm{\Gamma }}_x}$ (heights/z positions of the cores of CRVs, b), and $y_{{\mathrm{\Gamma }}_x}$ (lateral positions/y positions of the cores of CRVs, c). They are calculated based on Eqs. (13) and (14). The examples of the areas considered to calculate these quantities are enclosed by the dashed magenta lines in Fig. 13. See the text for more details. Note that the x ticks are made the same as the x positions of the MRSLs placed in different rows.

4.4.1 Strength of CRVs

As indicated by the values of $\left|{\mathrm{\Gamma }}_x\right|$ on the left of Fig. 15, for both the UW and the DW cases, the strength of CRVs gradually decreases with larger x before reaching the MRSL of the next row (e.g., $\left|{\mathrm{\Gamma }}_x\right|$ drops in the region of $\mathrm{7.0}\le x/D\le \mathrm{11.0}$), indicating that CRVs dissipate as they are convected downstream without further perturbation. However, the plot also reveals that the strength of the CRVs grows stronger as more rows of MRSLs are passed, surpassing the maximum values observed in previous rows. This indicates that the CRVs released by the MRSLs of different rows accumulate, reinforcing their strength row by row. Stronger CRVs result in stronger vertical flows, making the vertical advection process more significant. Furthermore, based on $\left|{\mathrm{\Gamma }}_x\right|$ in Fig. 15, it appears that the strength of the CRVs has not yet reached its maximum or asymptotic value at the fifth row, suggesting that the strength of the CRVs may continue to accumulate if the regenerative wind farms have more rows of MRSLs.

4.4.2 Positions of CRVs

The positions of the CRVs' cores are quantified using $z_{{\mathrm{\Gamma }}_x}$ and $y_{{\mathrm{\Gamma }}_x}$ introduced earlier. They are plotted in the middle and right of Fig. 15. $z_{{\mathrm{\Gamma }}_x}$ and $y_{{\mathrm{\Gamma }}_x}$ indicate the vertical positions (heights) and the lateral positions of the CRVs' cores, respectively.

The self-propelling property of CRVs can be observed by checking the values of $z_{{\mathrm{\Gamma }}_x}$. Specifically, in case UW, the CRVs' cores rise progressively as they move deeper into the regenerative wind farm. Conversely, in case DW, the heights of the CRVs' cores gradually decrease starting from the first row. However, the positions of $z_{{\mathrm{\Gamma }}_x}$ reach a minimum around the third row of the regenerative wind farm in case DW, after which they begin to rise. This is primarily due to the presence of the ground and the induced flow from the MRSLs in the side columns (see Fig. 2 for a contour illustrating all three columns).

Similarly to $z_{{\mathrm{\Gamma }}_x}$, the y positions (lateral positions) of CRVs, $y_{{\mathrm{\Gamma }}_x}$, also depend on the lifting configurations. In case UW, $y_{{\mathrm{\Gamma }}_x}$ consistently remains around $y/D=\mathrm{0.5}$ from the first row to the fifth row. Conversely, in case DW, $y_{{\mathrm{\Gamma }}_x}$ shifts increasingly outward as the flow travels in the positive x direction. This outward shift of the CRVs is due to the boundary condition imposed by the ground, which can be deduced through the method of image vortices (Saffman, 1995). The presence of the ground also influences the locations of $y_{{\mathrm{\Gamma }}_x}$ in case UW, causing them to tend toward $y/D=\mathrm{0.0}$ between any two consecutive rows. However, since $z_{{\mathrm{\Gamma }}_x}$ for case UW is located much higher than that in the DW case, the effects of the ground are much less significant.

An important aspect to mention is that when the positions of the CRVs' cores are higher (larger $z_{{\mathrm{\Gamma }}_x}$), they may be more capable of mixing the exhausted wakes with the fresh freestream, which could be more beneficial for the vertical entrainment process. Therefore, even when the vertical forcing of UW and DW is very similar, case UW may offer inherent advantages in this regard.

4.5 Control volume analysis on energy budget

In this subsection, the energy transport process of the simulated regenerative wind farms is explored using the control volume approach. The calculations are based on Eq. (H4), derived and explained in Appendix H. Six terms are considered, namely mean kinetic energy (MKE) advection, MKE diffusion plus pressure work, TKE advection plus diffusion, power extraction (by the MRSL's rotor), TKE dissipation, and residuals. Five control volumes (CVs) are examined, labeled from A to E (see Fig. 6). Each CV encloses an MRSL in the mid-column, covering a range from its 4 D upstream to its 2 D downstream. This range is designed to assess the energy sources and sinks of the MRSL in a specific CV. The vertical and lateral ranges are $z_{\mathrm{rc}}-\mathrm{0.5}D\le z\le z_{\mathrm{rc}}+\mathrm{0.5}D$ and $-\mathrm{2.5}D\le y\le \mathrm{2.5}D$, encompassing the entire mid-column. The results of the control volume analysis for the three cases listed in Table 3 are presented in Fig. 16.

For case WL (left of Fig. 16), it can be observed that the energy within the CV is primarily supplied by the advection of MKE and the (modeled) turbulent shear stress (MKE diffusion plus pressure work), with both terms having similar magnitudes. Notably, neither of these energy sources alone could supply the power extracted by the MRSL rotors. At first glance, it may appear that the advection of MKE continues to supply energy to the CVs in case WL. However, this is because the wakes of the MRSLs accumulate row by row, continuously depleting the MKE transported in the streamwise direction (freestream MKE), while the contribution from vertical advection is almost negligible. This interpretation is supported by the flow fields shown in the top panels of Figs. 10 and 11.

For the UW and DW cases (middle and right of Fig. 16), unlike the WL case, the contribution of MKE advection is much greater than the work done by turbulent shear stress, with MKE advection alone being sufficient to support the energy extraction by the MRSLs after the second row. Furthermore, the energy supplied to CVs by MKE advection in these two cases is primarily due to the vertical advection process, driven by the strong vertical velocity component (see Fig. 11). This vertical energy entrainment process differs significantly from conventional wind farms, which mostly rely on Reynolds shear stress (turbulent shear stress) (Porté-Agel et al., 2020; Calaf et al., 2010; VerHulst and Meneveau, 2014). Due to the indirect nature of energy entrainment by Reynolds shear stress, its magnitudes are naturally less than those of energy entrainment by advection. The latter directly injects higher energy flows into the control volumes, while the former relies on a secondary process involving Reynolds shear stress and the strength of the shear layer.

By closely inspecting Fig. 16, it can be observed that the MKE advection term is higher for case UW compared to case DW. However, the MKE diffusion plus pressure work term (representing the work done by turbulent shear stress) contributes negatively in case UW after the second row. This phenomenon can be explained by the lateral-averaged velocity profiles shown in Fig. 12, where the shear profiles for case UW are inverted (〈u〉_±2.5 D decreases with increasing z) in the range of $\mathrm{0.2}<z/D<\mathrm{1.5}D$, causing the shear to impact energy entrainment negatively. This highlights one of the key differences in aerodynamics between UW and DW.

Figure 16 Normalized energy transport rates of the terms in Eq. (H4) based on control volumes. The normalization is done by dividing the values by the rotor power of the MRSL situated at the first row in the mid-column of case WL, which is $P^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$. Cases WL, UW, and DW are plotted in the left, middle, and right, respectively. The letters at the abscissa refer to the indices of the control volumes (see Fig. 6), where volume A is the most upstream one, and volume E is the most downstream one. Each volume encloses an MRSL of the mid-column, covering its 4 D upstream to its 2 D downstream.

5 Conclusions and outlooks

This study conducted numerical investigations of regenerative wind farms. Regenerative wind farming is a newly proposed wind farm concept that consists of innovative wind harvesting systems, which are multi-rotor systems with lifting devices (MRSLs; see Fig. 1). In these regenerative wind farms, wake recoveries of MRSLs were engineered to be much faster compared to conventional horizontal-axis wind turbines (HAWTs), significantly reducing the power losses due to wake interactions. These enhanced wake recoveries were achieved by altering the vertical entrainment processes. Instead of turbulent mixing, the entrainment processes were facilitated by the vertical flows induced by trailing vortices generated by the lifting devices of MRSLs (see Fig. 2). To gain a comprehensive understanding of how the MRSLs' lifting devices affect the entrainment processes, cases with different configurations for lifting devices were tested. Additionally, detailed parametric and sensitivity studies on varying lifting magnitudes, inflow turbulence intensities, turbulence models, and grid sizes are performed in the appendixes.

Our results showed that, as the magnitudes of the vertical force were similar to the thrust of MRSLs, the power outputs of MRSLs with the lifting devices could be more than tripled compared to those without after the third row of the regenerative wind farms (see Fig. 9), diminishing the wake losses from around 75 % to 25 %. This significant increase in power output highlights the great potential of the regenerative wind farm. Specifically, to deliver the same amount of power, regenerative wind farms would require only half the land area compared to conventional wind farms with HAWTs because they have much higher wind farm efficiencies (see Table 4). This land use reduction could lower the overall cost of wind energy, making renewable energy more affordable.

Further examinations of how regenerative wind farms could achieve significantly higher power output were conducted by analyzing the flow fields using both qualitative and quantitative methods. Two-dimensional contour plots and three-dimensional isosurfaces illustrated that the low-velocity wakes of MRSLs were guided vertically upward, while high-velocity fresh flows were directed downward, replenishing the available power for MRSLs located further downstream. Circulation-based analysis revealed that the strength of counter-rotating vortices (CRVs), which are the trailing vortices generated by the MRSLs' lifting devices, accumulated progressively as the flow moved deeper into the regenerative wind farms. These CRVs are responsible for inducing the vertical advection process, with stronger CRVs leading to stronger vertical entrainment processes. Energy budget analysis based on control volumes indicated that wind farms with MRSLs equipped with lifting devices underwent a much stronger energy recovery than those with multi-rotor systems lacking such devices. Moreover, the analysis confirmed that the primary contributor to wake recovery in cases with lifting devices was the vertical advection process, contrasting with conventional wind farms, where wake recovery predominantly relies on turbulent shear (Calaf et al., 2010; VerHulst and Meneveau, 2015). These analyses thoroughly investigated the underlying physics of how regenerative wind farms can achieve significantly higher power outputs.

The results and analysis from this study suggest that the concept of regenerative wind farms could potentially lead to wind farms with much higher farm efficiencies than their conventional counterparts. A series of future research efforts is recommended to fully understand the physics and potential of regenerative wind farms and MRSLs, aiding in their realization. Several key aspects related to aerodynamics are outlined below. First, conducting simulations using higher-fidelity models, such as large-eddy simulations, in a time-resolved manner would enable better resolving the aerodynamics within the regenerative wind farms. Moreover, modeling MRSLs with greater detail is desirable, as the effects of their detailed geometry and the rotational of the sub-rotors are omitted in the current work. Additionally, investigating whether the stability properties of atmospheric boundary layers influence the dynamics of CRVs would be of significant interest. Furthermore, exploring how the inflow directions, the layouts of regenerative wind farms, and the sizes of the regenerative wind farms impact the farm efficiency is also critical. In addition, investigating the effectiveness of placing large-wing structures (i.e., MRSLs without rotors) between the turbines of existing wind farms is also of interest. This approach explores the possibility of transforming existing wind farms into regenerative wind farms. Moreover, experimental studies on regenerative wind farms and developing MRSL prototypes should be considered a top priority, as the ultimate goal is to transform this innovative concept into a real-world application. Certainly, there are numerous other practical challenges beyond aerodynamics, such as the structural integrity of MRSLs, the control strategies and mechanisms of MRSLs (such as yaw control), and the economic feasibility of regenerative wind farms, among others. These aspects are also critical, and addressing them adequately will be necessary to bring the concept of regenerative wind farms to a commercial stage.

Appendix A: A brief overview of the RANS k–ω SST model

The key governing equations of the RANS k–ω SST model are written in Eqs. (A1) to (A3), which are the equations for continuity, transport of momentum, and transport of modeled turbulence kinetic energy (denoted as TKE or k). In these equations, u_i, p, k, ω, ρ, ν, f_body,i, S_ij, τ_ij, and ν_T denote the ith component of velocity, static pressure, turbulence kinetic energy, turbulence specific dissipation, fluid density, kinematic (molecular) viscosity, the body forces applied on the flow, the shear strain tensor, the Reynolds stress tensor, and eddy viscosity. Note that all quantities just mentioned are time averaged. The definition of S_ij and the modeling of τ_ij are written in Eq. (A4). For brevity, certain equations related to the k–ω SST model, such as the transport equation of ω and the calculation of ν_T, have been omitted. Readers are referred to the OpenFOAM v2106 documentation (OpenCFD Ltd., 2021) for further details. Moreover, in this work, all model coefficients of the k–ω SST model are set to the default values provided by OpenFOAM v2106 (OpenCFD Ltd., 2021), e.g., ${\mathit{\beta }}^{*}=\mathrm{0.09}$.

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial u_i}{\partial x_i}}=\mathrm{0}\text{(A2)}& {\displaystyle }{u}_j{\displaystyle \frac{\partial u_i}{\partial x_j}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial p}{\partial x_i}}+{\displaystyle \frac{\partial }{\partial x_j}}(\mathrm{2}\mathit{\nu }{S}_{ij}+{\mathit{\tau }}_{ij})+{\displaystyle \frac{f_{\mathrm{body},i}}{\mathit{\rho }}}\text{(A3)}& {\displaystyle }{u}_j{\displaystyle \frac{\partial k}{\partial x_j}}={\mathit{\tau }}_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}-{\mathit{\beta }}^{*}\mathit{\omega }k+{\displaystyle \frac{\partial }{\partial x_j}}\left[(\mathit{\nu }+{\mathit{\nu }}_{\mathrm{T}}){\displaystyle \frac{\partial k}{\partial x_j}}\right]\text{(A4)}& {\displaystyle }{S}_{ij}\stackrel{\Delta }{=}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right),{\mathit{\tau }}_{ij}=\mathrm{2}{\mathit{\nu }}_{\mathrm{T}}{S}_{ij}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}k{\mathit{\delta }}_{ij}\end{array}$$

Appendix B: Sensitivity test of turbulence models

For the steady RANS simulations, all the fluctuating properties are modeled through the turbulence model, and they are modeled differently depending on the model chosen. Thus, model-related uncertainties arise, which may affect the conclusions obtained by analyzing CFD results. To ensure that the conclusions obtained in this work are robust and independent of the chosen turbulence model, a few simulations are conducted with several mostly used turbulence models. In addition to the already-used k–ω SST model (Menter, 1994), configurations of cases WL, UW, and DW in Table 3 are tested with the realizable k–ε model (Shih et al., 1995) and RNG k–ε model (Yakhot et al., 1992). These are three of the most popular turbulence models for wind-energy-related applications, and note that there is currently no unified standard for the optimum RANS turbulence model (Thé and Yu, 2017; Eidi et al., 2021).

Figure B1 Simulations cases with different RANS turbulence models. The configurations used are those of the cases listed in Table 3. $\langle {\widehat{P}}^R\rangle$, $\langle {\widehat{T}}^R\rangle$, and $\langle {\widehat{L}}^W\rangle$ are the normalized row-averaged power, thrust, and vertical force, respectively. The normalization factors are $P^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$ and $T^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$ of the cases using the k–ω SST model, where their values are the same as those used in Figs. 8 and 9.

Figure B2 Contour plots of (modeled) eddy viscosity ν_T and streamwise velocity u for the cases with different RANS turbulence models. The two quantities are plotted side by side in each panel, where the left is ν_T and the right is u. The top, middle, and bottom rows are the solution given by the simulations using the k–ω SST, realizable k–ε, and RNG k–ε models, respectively. The leftmost column presents the inflow conditions at $x/D=-\mathrm{2.0}$ of configuration WL for reference. The rest of the columns plot the quantities at $x/D=\mathrm{4.0}$, with the cases from left to right being WL, UW, and DW. The arrows in the plots represent the in-plane velocities. The lengths of the arrows are scaled by their norms, and their absolute scales are indicated at the bottom right of the figure.

The results of MRSLs' outputted power, thrust, and lift for cases with different turbulence models are presented in Fig. B1. Except for changing the turbulence model, all other parameters remain the same as cases WL, UW, and DW listed in Table 3 (inlet conditions of ε are set to atmBoundaryLayerInletEpsilon). It can be seen that similar results are yielded by the RNG k–ε model when compared to the k–ω SST model. As for the results with the realizable k–ε model, the results significantly deviate from the other two, and it shows that the effectiveness of the wings of the MRSL is less astonishing. This deviation is likely due to the more diffusive behavior of the realizable k–ε model in the current application.

As shown in Fig. B2, the ν_T values modeled by the realizable k–ε model are higher than those of the other two turbulence models across all three configurations. Higher ν_T leads to faster vortex dissipation and increased (modeled) turbulent shear stress (term MKE diffusion in Eq. H4). Faster vortex dissipation reduces the effects of upwash and downwash (see the arrows in Fig. B2), thereby slowing the wake recovery rates for configurations UW and DW. Conversely, higher (modeled) turbulent shear stress (MKE diffusion) accelerates wake recovery for configuration WL, as this is its primary mechanism for wake recovery (see Fig. 16).

Although the power predicted by different turbulence models shows considerable variation, it is evident that significant improvements are achieved by introducing lifting devices when comparing the power outputs of the cases with and without lifting devices across all three turbulence models. Thus, with the results, it can be concluded that even though the selection of the RANS turbulence model may influence the results in the sense of their absolute values, it does not significantly impact the main conclusion of this work, which is that the installation of the lifting devices (wings) can dramatically improve the power performance of the downstream wind harvesting systems (MRSLs) in regenerative wind farms.

Appendix C: Grid independence test

A grid independence test is carried out to ensure that the discretization error of the CFD simulations does not affect the conclusions drawn. The three cases in Table 3, WL, UW, and DW, are tested with three grid sizes Δ. The three tested grid sizes are $\mathrm{\Delta }=D/\mathrm{20}$, $\mathrm{\Delta }=D/\mathrm{25}$, and $\mathrm{\Delta }=D/\mathrm{30}$, and they are labeled as Coarse, Medium, and Fine, and each of them results in a mesh that has 73.9 M, 137.8 M, and 249.5 M cells, respectively. Note that the cases in Table 3 use mesh Medium. Also note that, except for adjustment of the grid sizes, all other parameters are kept the same as the cases in Table 3, including the spacings of the actuator elements for wings and the absolute values of the smearing factors (ε^R and ε^W).

The results of the grid independence test are presented in Fig. C1, where 〈ΔP^R〉, 〈ΔT^R〉, and 〈ΔL^W〉 are the relative deviations of 〈P^R〉, 〈T^R〉, and $\langle \left|L^W\right|\rangle$ from their reference cases, respectively. The reference cases are the cases that used a medium mesh. The definition of 〈ΔP^R〉 is given in Eq. (C1), and 〈ΔT^R〉 and 〈ΔL^W〉 are derived in the same way. It can be seen that 〈ΔT^R〉 and 〈ΔL^W〉 of the cases with coarse and fine meshes both fall in the ranges of ±2 % for the first row and ±4 % for all rows. This suggests that the impacts of grid sizes considered are minimal on the conclusions drawn, as the deviations due to different Δ are at least an order smaller than the differences in T^R between the cases with and without lifting devices (∼50 % to 100 %). For the values of 〈ΔP^R〉 of the WL cases, although they can be up to 6 % for the third, fourth, and fifth rows, their absolute values of the reference 〈P^R〉 are relatively small compared to the upstream rows. Due to their relatively small values compared to those of the upstream rows, their relative deviations are more susceptible to the variations arising from the upstream rows. With these results, it can be concluded that the medium mesh ($\mathrm{\Delta }=D/\mathrm{25}$) is sufficient for the application used in this work.

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& \langle \mathrm{\Delta }{P}^R\rangle \mathrm{of}i\mathrm{th}\mathrm{row}\\ & \stackrel{\Delta }{=}{\displaystyle \frac{\langle P^R\rangle \mathrm{of}i\mathrm{th}\mathrm{row}-\langle P^R\rangle \mathrm{of}i\mathrm{th}\mathrm{row}\mathrm{with}\mathrm{mesh}\mathrm{medium}}{\langle P^R\rangle \mathrm{of}i\mathrm{th}\mathrm{row}\mathrm{with}\mathrm{mesh}\mathrm{medium}}}\end{array}\end{array}$$

Figure C1 Grid independence test with the three cases listed in Table 3, which are the WL, UW, and DW cases. 〈ΔP^R〉, 〈ΔT^R〉, and 〈ΔL^W〉 are the relative deviations of 〈P^R〉, 〈T^R〉, and $\langle \left|L^W\right|\rangle$ from the cases with a medium mesh, where their definitions are in Eq. (C1).

Appendix D: Sensitivity test of inflow turbulence intensity

Previous studies have shown that inflow turbulence intensity (TI) significantly influences wind turbine wake aerodynamics, particularly by disrupting large-scale coherent structures in the wake and reducing their coherence (Li et al., 2024a). In this work, the concept of regenerative wind farms relies heavily on the swirling motions induced by the trailing vortices released from MRSLs, which are also considered coherent structures. Therefore, a parametric study on inflow TI is particularly important and worth conducting. In this appendix, additional simulations are performed with inflow TI values of 5 % and 14 %. Except for the changes in inflow turbulence intensity, all other parameters are kept identical to cases WL, UW, and DW listed in Table 3. To achieve TI levels of 5 % and 14 %, the model coefficient C₁ mentioned in Sect. 3.2 is set to 0.065 and 8.214, respectively.

The results for the power, thrust, and lift output of MRSLs under different inflow TIs are presented in Fig. D1, while the flow fields are analyzed in Fig. D2. As expected, Fig. D1 shows that the effectiveness of the lifting devices diminishes as TI increases. Higher TI causes the trailing vortices in UW and DW to dissipate more rapidly, while it enhances the wake recovery rates in WL by increasing the diffusion of MKE (modeled turbulent mixing). This trend is confirmed by examining the flow fields in Fig. D2, which demonstrate that the upwash and downwash effects in UW and DW become less pronounced with higher TI.

Although the improvement in power outputs is less dramatic under higher TI, this appendix shows that the benefits of the lifting devices remain significant even at TI =14 %, which is around or beyond the upper limit of typical offshore conditions (Hansen et al., 2012). Specifically, MRSLs in UW and DW still achieve over 50 % higher power output compared to WL after the third row when subjected to an inflow TI of 14 %.

Figure D1 Simulation cases that are subjected to different inflow turbulence intensities. The configurations used are those of the cases listed in Table 3. $\langle {\widehat{P}}^R\rangle$, $\langle {\widehat{T}}^R\rangle$, and $\langle {\widehat{L}}^W\rangle$ are the normalized row-averaged power, thrust, and vertical force, respectively. The normalization factors are $P^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$ and $T^R{|}_{\mathrm{1}\mathrm{st},\mathrm{mid}}^{\mathrm{WL}}$ of the cases subjected to TI =8 %, where their values are the same as those in Figs. 8 and 9.

Figure D2 Contour plots of turbulent intensity (TI) and streamwise velocity u for the cases subjected to different inflow TIs. The two quantities are plotted side by side in each panel, where the left is TI and the right is u. Note that TI is defined as $\sqrt{\mathrm{2}k/\mathrm{3}}/u_{\mathrm{ref}}$, where k is the turbulent kinetic energy modeled by the k–ω SST model. The top, middle, and bottom rows are the solution given by the simulations subjected to inflow TI being 5 %, 8 %, and 14 %, respectively. The leftmost column presents the inflow conditions at $x/D=-\mathrm{2.0}$ of configuration WL for reference. The rest of the columns plot the quantities at $x/D=\mathrm{4.0}$, where cases WL, UW, and DW are indicated from left to right. The arrows in the plots represent the in-plane velocities. The lengths of the arrows are scaled by their norms, and their absolute scales are indicated at the bottom right of the figure.

Appendix E: Forcing distribution of MRSLs' wings

The angle of attack and load distributions of the MRSLs' wings for cases UW and DW in Table 3 are presented in Figs. E1 and E2, respectively. Note that the presented angle of attack α, streamwise forcing $f_x^{\mathrm{AL}}$, and vertical forcing $f_z^{\mathrm{AL}}$ are sampled from the mid-column of the wind farms. Definitions of α, $f_x^{\mathrm{AL}}$, and $f_z^{\mathrm{AL}}$ are in Eqs. (E1) and (E2), where f^AL is the force exerted by an actuator element of the wings. Focusing on the force exerted by the first row of the MRSL, tip losses can be identified with the $f_z^{\mathrm{AL}}$ profiles. Additionally, with the $f_x^{\mathrm{AL}}$ profiles, it can be seen that induced drags are mainly concentrated around the tips. These results comply with classical aerodynamic theories (Anderson, 2011), suggesting that the wings' loading predicted by the actuator lines used in this work is reasonable. Some peculiar shapes appear for the loading from the second row onward. This is due to the wakes and vertical flows introduced by upstream MRSLs that complicate the inflow of these wings. Furthermore, in these two figures, it can be seen that α in the middle of the wings is 12.5°, and this α corresponds to C_l=2.5 according to the airfoil polar data of the S1223 airfoil (Fig. 5), confirming $C_{l,\mathrm{mid}}=\mathrm{2.5}$ holds for these two cases.

$$\begin{array}{}\text{(E1)}& {\displaystyle }\mathit{\alpha }\stackrel{\Delta }{=}\left\{\begin{array}{ll}-w/u& \mathrm{for}\mathrm{the}\text{up-washing}\mathrm{cases}\\ \phantom{-}w/u& \mathrm{for}\mathrm{the}\text{down-washing}\mathrm{cases}\end{array}\right.\text{(E2)}& {\displaystyle }{\mathit{f}}^{\mathrm{AL}}=f_x^{\mathrm{AL}}{\widehat{\mathit{e}}}_x+f_z^{\mathrm{AL}}{\widehat{\mathit{e}}}_z\end{array}$$

Figure E1 The profiles of the angle of attack α (top row), vertical forcing $f_z^{\mathrm{AL}}$ (middle row), and streamwise forcing $f_x^{\mathrm{AL}}$ (bottom) for the wings of the MRSLs situated in the mid-column for case UW.

Figure E2 The profiles of the angle of attack α (top row), vertical forcing $f_z^{\mathrm{AL}}$ (middle row), and streamwise forcing $f_x^{\mathrm{AL}}$ (bottom) for the wings of the MRSLs situated in the mid-column for case DW.

Appendix F: Testing MRSLs with different lift magnitudes

To further understand how the magnitudes of MRSLs' lift affect the performance of regenerative wind farms, several auxiliary cases are performed. The cases tested are listed in Table F1. Three different vertical forcing (lifting) magnitudes are tested for each direction of the lift. The lift magnitudes are adjusted by changing C_l,mid (the lift coefficient at the mid-span of the wing) by pitching the wings (see the end of Sect. 3.4). Both directions of lift are tested with C_l,mid being 0.5, 1.5, and 2.5. Note that the cases marked with an asterisk are those listed in Table 3, where WL, U2_5, and D2_5 correspond to WL, UW, and DW, respectively.

Similarly to Fig. 8, Fig. F1 presents the normalized row-averaged thrust, lift, and induced drag ($\langle {\widehat{T}}^R\rangle$, $\langle {\widehat{L}}^W\rangle$, and $\langle {\widehat{D}}_{\mathrm{ind}}^W\rangle$) of the MRSLs for the seven cases listed in Table F1. As designed, $\langle {\widehat{L}}^W\rangle$ increases with higher C_l,mid values. Additionally, regardless of the lift direction, $\langle {\widehat{T}}^R\rangle$ after the second row is higher for all the cases with lifting devices compared to the case without. Moreover, from the second row onward, higher C_l,mid corresponds to higher $\langle {\widehat{T}}^R\rangle$ for both directions of the lift, despite the fact that larger C_l,mid also results in larger $\langle {\widehat{D}}_{\mathrm{ind}}^W\rangle$, as shown in the right panel of Fig. F1.

Table F1 The test matrix of the auxiliary cases. These cases have different values of C_l,mid, which affects the magnitudes of the MRSLs' vertical force (lift). The values of C_l,mid represent C_l at the mid-span of the MRSLs' wings of a case (see Sect. 3.4). Note that the cases marked with an asterisk are the same as those listed in Table 3, where WL, U2_5, and D2_5 correspond to the WL, UW, and DW cases, respectively.

Figure F2 summarizes the power performance of the regenerative wind farms for the cases listed in Table F1. It can be seen that the normalized row-averaged power output of the MRSLs (${\widehat{P}}^R$) progressively increases with higher values of C_l,mid, indicating that the performance of the regenerative wind farms is positively correlated with the lift magnitudes of the MRSLs within the tested range.

Figure F1 The normalized row-averaged thrust of the MRSL's rotor (${\widehat{T}}^R$, a) together with the vertical (${\widehat{L}}^W$, b) and streamwise (${\widehat{D}}_{\mathrm{ind}}^W$, c) force components of the MRSL's lifting devices. The normalization is done by dividing the reference rotor thrust, which is based on the MRSL at the first row in the mid-column of the WL case. The legends correspond to the case number introduced in Table F1.

Figure F2 The normalized row-averaged rotor power of the MRSL's rotor (${\widehat{P}}^R$) for the auxiliary cases in Table F1. The normalization is done by dividing the values by the rotor power of the MRSL situated at the first row of the mid-column of the WL case.

Appendix G: Frandsen wake model

The analytical wake model used in Sect. 4.2 is the well-known Frandsen model, which is proposed by Frandsen et al. (2006) (region I is used as the wakes of different columns do not merge). It is derived from momentum analysis over a control volume covering one or multiple wind turbines aligned in the streamwise direction. The inputs of the Frandsen model include the wind turbine diameter D_F, the thrust coefficient C_T, and the streamwise spacings between the turbines (when there is more than one row of turbines), and the outputs are the wake velocity u_F and wake diameter D_w, which vary along the streamwise direction. The equations for the Frandsen model are briefly written in Eqs. (G1) and (G2), where α_F is the wake expansion factor that is decided empirically. In the current work, α_F=0.0629 is used, based on the CFD results using large-eddy simulations reported by Andersen et al. (2014).

Note that the Frandsen model was developed mainly for horizontal-axis wind turbines, and thus it is not immediately suitable for the current work, as the MRSL has a square frontal area instead of a circular frontal area. Therefore, a correction is needed. In this work, $D_{\mathrm{F}}=D_{\mathrm{cir}}=\mathrm{2}D/\sqrt{\mathit{\pi }}$ is used, making $\mathrm{0.25}\mathit{\pi }{D}_{\mathrm{F}}^{\mathrm{2}}=D^{\mathrm{2}}$ (D is the side length of MRSL). That is, for a circular disk with D_F being the diameter, its swept area will be equal to the one for a square MRSL used in this work.

Since the current work is only interested in the velocity at x positions where there are MRSLs, for simplicity, only the velocity at these positions is calculated. u_F,nth is used to denote the inflow velocity seen by the nth row of the MRSL predicted by the Frandsen model, and x_nth is the streamwise position of the nth row. D_w,nth denotes the wake diameter at x=x_nth. For clarity, u_F of the first and second rows is explicitly written in Eq. (G3). u_F after the third row is calculated through a recursive method using Eq. (G4). After obtaining the values of u_F for all the rows, the relation of Eq. (G5) is utilized to obtain $P_{\mathrm{F}}^R$ (the power output of the MRSL's rotors predicted by the Frandsen model) based on the one-dimensional momentum theory (Manwell et al., 2010). Note that the values for C_T and C_P are 0.7 and 0.54, as mentioned in Sect. 2, and the corresponding power of the first row is 29.9 MW.

$$\begin{array}{}\text{(G1)}& {\displaystyle }{D}_{\mathrm{w}}\left(x\right)=D_{\mathrm{F}}{\left(\mathit{\beta }+{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{F}}x}{D_{\mathrm{F}}}}\right)}^{\mathrm{1}/\mathrm{2}},\mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T}}}}{\sqrt{\mathrm{1}-C_{\mathrm{T}}}}}\text{(G2)}& {\displaystyle }{u}_{\mathrm{F}}=u_{\mathrm{ref}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\mathrm{1}-\mathrm{2}{\left({\displaystyle \frac{D_{\mathrm{F}}}{D_{\mathrm{w}}}}\right)}^{\mathrm{2}}{C}_{\mathrm{T}}}\right)\text{(G3)}& {\displaystyle }\begin{array}{rl}& u_{\mathrm{F},\mathrm{1}\mathrm{st}}=u_{\mathrm{ref}}\\ & u_{\mathrm{F},\mathrm{2}\mathrm{nd}}=u_{\mathrm{ref}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\mathrm{1}-\mathrm{2}{\left({\displaystyle \frac{D_{\mathrm{F}}}{D_{\mathrm{w},\mathrm{2}\mathrm{nd}}}}\right)}^{\mathrm{2}}{C}_{\mathrm{T}}}\right)\end{array}\text{(G4)}& {\displaystyle }\begin{array}{rl}{u}_{\mathrm{F},n\mathrm{th}}& =u_{\mathrm{ref}}-\left[{\left({\displaystyle \frac{D_{\mathrm{w},(n-\mathrm{1})\mathrm{th}}}{D_{\mathrm{w},n\mathrm{th}}}}\right)}^{\mathrm{2}}\left(u_{\mathrm{ref}}-u_{\mathrm{F},(n-\mathrm{1})\mathrm{th}}\right)\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{D_{\mathrm{F}}}{D_{\mathrm{w},n\mathrm{th}}}}\right)}^{\mathrm{2}}{C}_{\mathrm{T}}{u}_{\mathrm{F},(n-\mathrm{1})\mathrm{th}}\right],\mathrm{for}n\ge \mathrm{3}\mathrm{and}x=x_{n\mathrm{th}}\end{array}\text{(G5)}& {\displaystyle }{P}_{\mathrm{F},n\mathrm{th}}^R=\mathrm{0.5}\mathit{\rho }{\left(u_{\mathrm{F},n\mathrm{th}}\right)}^{\mathrm{3}}{D}^{\mathrm{2}}{C}_P\propto {\left(u_{\mathrm{F},n\mathrm{th}}\right)}^{\mathrm{3}}\end{array}$$

Appendix H: Control volume analysis of flow energies

In addition to evaluating the performance of MRSLs based on their power outputs, analysis of this work is also conducted using the terms of the energy transport equations based on the control volume approach. This analysis aims to distinguish the primary source terms for wake recoveries. In this appendix, the control volume approach used to analyze the budget of flow energy is detailed.

We start with the two transport equations for mean kinetic energy (MKE; also denoted as K) and turbulence kinetic energy (TKE; also denoted as k) provided in Eqs. (H1) and (H2), where the physical meaning of each term is labeled. The definition of MKE is given in Eq. (H3). The transport equations for MKE (Eq. H1) and TKE (Eq. H2) are derived by multiplying ρu_i and ρ by Eqs. (A2) and (A3), respectively. Some rearrangements are then performed using the chain rule and the continuity equation (Eq. A1).

$$\begin{array}{}\text{(H1)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{advection}\mathrm{of}\mathrm{MKE}}{\underbrace{{\displaystyle \frac{\partial u_{j}K}{\partial x_j}}}}& =\underset{\mathrm{work}\mathrm{done}\mathrm{by}\mathrm{surface}\mathrm{forces}}{\underbrace{-{\displaystyle \frac{\partial u_{j}p}{\partial x_j}}+\mathit{\rho }{\displaystyle \frac{\partial }{\partial x_j}}\left[u_i\right(\mathrm{2}\mathit{\nu }{S}_{ij}+{\mathit{\tau }}_{ij}\left)\right]}}\\ & \underset{\mathrm{viscous}\mathrm{dissipation}}{\underbrace{-\mathrm{2}\mathit{\rho }\mathit{\nu }{S}_{ij}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}\right)}}\underset{\mathrm{transfer}\mathrm{to}\mathrm{TKE}}{\underbrace{-\mathit{\rho }{\mathit{\tau }}_{ij}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}\right)}}\\ & \underset{\mathrm{work}\mathrm{done}\mathrm{by}\mathrm{body}\mathrm{forces}}{\underbrace{+\mathit{\rho }{u}_i{f}_{\mathrm{body},i}}}\end{array}\text{(H2)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{advection}\mathrm{of}k}{\underbrace{{\displaystyle \frac{\partial u_j\mathit{\rho }k}{\partial x_j}}}}& =\underset{\mathrm{transfer}\mathrm{from}\mathrm{MKE}}{\underbrace{\mathit{\rho }{\mathit{\tau }}_{ij}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}\right)}}\underset{\mathrm{dissipation}\mathrm{of}k}{\underbrace{-\mathit{\rho }{\mathit{\beta }}^{*}\mathit{\omega }k}}\\ & \underset{\mathrm{diffusion}\mathrm{of}k}{\underbrace{+{\displaystyle \frac{\partial }{\partial x_j}}\left[\mathit{\rho }(\mathit{\nu }+{\mathit{\nu }}_{\mathrm{T}}){\displaystyle \frac{\partial k}{\partial x_j}}\right]}}\end{array}\text{(H3)}& {\displaystyle }K\stackrel{\Delta }{=}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{u}_i{u}_i\end{array}$$

The transport equation for the total energy (resolved plus modeled) in differential form can be obtained by adding the two energy equations (Eqs. H1 and H2). This equation is integrated over a control volume (CV) to examine the energy balance, resulting in Eq. (H4). The divergence theorem is applied, with CS denoting the control surface bounding the CV. It is worth noting that the term MKE diffusion plus pressure work essentially represents the work done by surface forces on the control volumes. Due to the high Reynolds number in this study (e.g., ${\mathit{Re}}_{\mathrm{c}}=u_{\mathrm{ref}}c/\mathit{\nu }>{\mathrm{10}}^{\mathrm{7}}$ and ${\mathit{Re}}_D=u_{\mathrm{ref}}D/\mathit{\nu }>{\mathrm{10}}^{\mathrm{8}}$), the primary contributor of MKE diffusion plus pressure work is the turbulent shear stress, which is modeled through the Reynolds shear stress τ_ij (see Eq. A4). That is, the effects of MKE diffusion in this work can be understood as the effects of the modeled turbulent mixing. Additionally, the signs for each term on the left-hand side of the equation are rearranged so that positive values correspond to energy gains for a CV and vice versa for the terms on the right-hand side. The viscous dissipation term is omitted because ν≪ν_T due to the high Reynolds number. Furthermore, a residual term ℛ is introduced to account for discrepancies, including the viscous dissipation term, losses due to the parasitic drag of the wings, errors from discretization, interpolation errors, and other factors.

$$\begin{array}{}\text{(H4)}& \begin{array}{rl}\mathit{\rho }\underset{\mathrm{CS}}{\oint }& \left(\underset{\mathrm{MKE}\mathrm{advection}}{\underbrace{-u_{j}K}}\underset{\mathrm{MKE}\mathrm{diffusion}\mathrm{plus}\mathrm{pressure}\mathrm{work}}{\underbrace{+\left[u_i(\mathrm{2}\mathit{\nu }{S}_{ij}+{\mathit{\tau }}_{ij})-{\displaystyle \frac{u_{j}p}{\mathit{\rho }}}\right]}}\right.\\ & \left.\underset{\mathrm{TKE}\mathrm{advection}\mathrm{plus}\mathrm{diffusion}}{\underbrace{-u_{j}k+\left[(\mathit{\nu }+{\mathit{\nu }}_{\mathrm{T}}){\displaystyle \frac{\partial k}{\partial x_j}}\right]}}\right)\mathrm{d}{S}_j\\ & =\underset{\mathrm{CV}}{\int }\left(\underset{\mathrm{Power}\mathrm{extraction}}{\underbrace{-u_i{f}_{\mathrm{body},i}^R}}\underset{\mathrm{TKE}\mathrm{dissipation}}{\underbrace{+\mathit{\rho }{\mathit{\beta }}^{*}\mathit{\omega }k}}\right)\mathrm{d}V+\underset{\mathrm{Residuals}}{\underbrace{\mathcal{R}}}\end{array}\end{array}$$