2.1. Theoretical derivation of PG-IMCM

2.1.1. Total momentum deficit under pressure gradient

Here, we consider a wind farm with $N$ wind turbines $(WT_1,\ldots,WT_i,\ldots,WT_N)$ under pressure gradient with the base flow $U_b(x)$, as depicted in figure 1. The wind turbines are assumed to have the same rotor diameter $D$. The position of $WT_i$ is labelled as $(x_i,y_i,z_i)$, where $x$, $y$ and $z$ are the streamwise, spanwise and vertical coordinates, respectively. Wind turbines are ordered with respect to their streamwise positions to ensure $x_i\geq x_{i-1}$. The RANS streamwise momentum equation for wind-farm flows can be written as follows:

(2.1)\begin{equation} U_w\frac{\partial U_w}{\partial x}+V\frac{\partial U_w}{\partial y}+W\frac{\partial U_w}{\partial z}={-}\frac{1}{\rho}\frac{\partial P}{\partial x}-\frac{\partial \overline{u^\prime u^\prime}}{\partial x}-\frac{\partial \overline{u^\prime v^\prime}}{\partial y}-\frac{\partial \overline{u^\prime w^\prime}}{\partial z}+\frac{f_t}{\rho}. \end{equation}

Here, $U_w$, $V$ and $W$ are the streamwise, spanwise and vertical velocity of the wind-farm flow, respectively, $\overline {u^\prime u^\prime }$, $\overline {u^\prime v^\prime }$ and $\overline {u^\prime w^\prime }$ are the Reynolds stress components, $P$ is the time-averaged static pressure, $\rho$ is the air density and $f_t$ is the body force caused by wind turbines, which can be expressed as (Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021):

(2.2)\begin{equation} \frac{f_t}{\rho}=\sum_{i=1}^N-\frac{T_i}{\rho\pi R^2}\delta(x-x_i)\mathcal{H}(R^2-[(y-y_i)^2+(z-z_i)^2]), \end{equation}

where $T_i$ is the thrust of $WT_i$, $R$ is the rotor radius, $\delta (x)$ is the Dirac delta function and $\mathcal {H}(x)$ is the Heaviside step function. The viscous term is neglected due to the high Reynolds number in the atmospheric boundary layer. It is noteworthy that we only consider the background pressure variations induced by complex terrain and neglect the sudden pressure changes in the near-wake region of the wind rotor (Dar & Porté-Agel Reference Dar and Porté-Agel2022). Therefore, the streamwise pressure gradient $\partial P/\partial x$ and the streamwise variation of the background velocity $U_b(x)$ satisfy the Bernoulli equation under the considered conditions, i.e. $1/\rho \partial P/\partial x=-({\rm d}U_b/{{\rm d}\kern0.7pt x})U_b$. The governing equation for the wind-farm flow velocity deficit $U_s(x,y,z)=U_b(x)-U_w(x,y,z)$ can be derived by subtracting (2.1) from $U_w \partial U_b/\partial x$, with the assumption of $\partial U_b/\partial y=\partial U_b/\partial z=0$:

(2.3)\begin{equation} U_w\frac{\partial U_s}{\partial x}+V\frac{\partial U_s}{\partial y}+W\frac{\partial U_s}{\partial z}={-}\frac{{\rm d} U_b}{{\rm d}\kern0.7pt x}U_s+\frac{\partial \overline{u^\prime u^\prime}}{\partial x}+\frac{\partial \overline{u^\prime v^\prime}}{\partial y}+\frac{\partial \overline{u^\prime w^\prime}}{\partial z}-\frac{f_t}{\rho}. \end{equation}

In addition, considering the continuity equation,

(2.4)\begin{equation} \frac{\partial U_w}{\partial x}+\frac{\partial V}{\partial y}+\frac{\partial W}{\partial z}=0, \end{equation}

we can obtain the following:

(2.5)\begin{equation} U_w\frac{\partial U_s}{\partial x}+V\frac{\partial U_s}{\partial y}+W\frac{\partial U_s}{\partial z}=\frac{\partial U_w U_s}{\partial x}+\frac{\partial VU_s}{\partial y}+\frac{\partial WU_s}{\partial z}. \end{equation}

Substituting (2.5) into (2.3) yields

(2.6)\begin{equation} \frac{\partial U_wU_s}{\partial x}+\frac{\partial VU_s}{\partial y}+\frac{\partial WU_s}{\partial z}={-}\frac{{\rm d} U_b}{{\rm d}\kern0.7pt x}U_s+\frac{\partial \overline{u^\prime u^\prime}}{\partial x}+\frac{\partial \overline{u^\prime v^\prime}}{\partial y}+\frac{\partial \overline{u^\prime w^\prime}}{\partial z}-\frac{f_t}{\rho}. \end{equation}

Then, we can take the $y$–$z$ plane integral from $x_0$ (a location upstream of the wind farm) to $x$ (an arbitrary location within or downstream of the wind farm). The upper and lower bounds of the $y$- and $z$-integration are selected far from the wind farm to ensure $U_s$, $\overline {u^\prime v^\prime }$ and $\overline {u^\prime w^\prime }$ all vanish (in a given $x$-normal plane). Therefore, the integrals of the terms inclusive of the $y$- or $z$-derivative in (2.6) are zero. In addition, the integral of normal Reynolds stress component $\iint \overline {u^\prime u^\prime }\,{\rm d}\kern 0.05em y\,{\rm d}z|_{x_0}^x$ is negligible, compared with other terms according to the budget comparison of Bastankhah et al. (Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021). Therefore, the final expression for the total momentum deficit under pressure gradient is

(2.7)\begin{align} &\iint U_w(x,y,z)[U_b(x)-U_w(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z\nonumber\\ &\quad =\sum_{i\in B}\frac{T_i}{\rho}-\iint \underbrace{\int_{x_0}^{x}\frac{{\rm d}U_b(x^\prime)}{{{\rm d}\kern0.7pt x}^\prime}U_s(x^\prime, y,z)\,{{\rm d}\kern0.7pt x}^\prime}_{{\mathfrak P}(x,y,z)} \,{{\rm d}\kern0.05em y}\,{\rm d}z. \end{align}

Here, $B=\{i\,|\,x_i< x\}$ and ${\mathfrak P}(x,y,z)$ is the pressure correction term induced by the wind-farm velocity deficit. The equation shows that the total momentum deficit for a wind farm under pressure gradient consists of two parts: the total thrust of wind turbines and the integral of ${\mathfrak P}(x,y,z)$, which is related to the pressure gradient.

Figure 1. Schematic of a wind farm with an arbitrary layout consisting of $N$ wind turbines under FPG. (Dark blue represents the high-pressure region, whereas light blue represents the low-pressure region.)

2.1.2. General expression for the wind-turbine thrust under pressure gradient

To proceed with our derivation, we introduce some new velocity variables. To begin with, we give the rules of nomenclature as follows: $u^i$ represents the velocity variables of $WT_i$, whereas $U^i$ represents the velocity variables associated with $i$ wind turbines; The subscripts $s$, $w$ and $b$ represent the velocity deficit, the wake velocity and the base flow, respectively. By definition, $U_b^i=U_b$. Subsequently, we choose two control volumes, including $i$ and $i-1$ wind turbines, whose upper bound of the streamwise coordinate is larger than $x_i$. The only difference is that the latter excludes $WT_i$, as depicted in figure 2. In this condition, if we apply (2.7) to these two control volumes and subtract one from the other, we can obtain the expression for the thrust of $WT_i$ as follows:

(2.8)\begin{align} \frac{T_i}{\rho}&=\iint U_w^i(x,y,z)[U_b(x)-U_w^i(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z\nonumber\\ &\quad -\iint U_w^{i-1}(x,y,z)[U_b(x)-U_w^{i-1}(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z \nonumber\\ &\quad +\iint \int_{x_i-\delta x}^{x}\frac{{\rm d}U_b(x^\prime)}{{{\rm d}\kern0.7pt x}^\prime}[U_w^{i-1}(x^\prime, y,z)-U_w^i(x^\prime, y,z)]\,{{\rm d}\kern0.7pt x}^\prime\, {{\rm d}\kern0.05em y} \,{\rm d}z, \end{align}

where $\delta _x$ represents an infinitely small distance, $U_w^i$ and $U_w^{i-1}$ are the streamwise velocity of the wind-farm flow with $i$ wind turbines and $i-1$ wind turbines, respectively. Here $U_w^{i-1}-U_w^i$ is essentially the wake velocity deficit $u_s^i$ induced by $WT_i$, which can be modelled by the wake model under pressure gradient. As has been demonstrated by Bastankhah et al. (Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021), the first and second terms on the right-hand side of (2.8) can be simplified as follows:

(2.9)\begin{align} &\iint U_w^i[U_b(x)-U_w^i]\,{\rm d}\kern0.05em y\,{\rm d}z-\iint U_w^{i-1}[U_b(x)-U_w^{i-1}]\,{\rm d}\kern0.05em y\,{\rm d}z\nonumber\\ &=\iint U_w^i(U_w^{i-1}-U_w^i)\,{\rm d}\kern0.05em y\,{\rm d}z-\iint (U_w^{i-1}-U_w^i)[U_b(x)-U_w^{i-1}]\,{\rm d}\kern0.05em y\,{\rm d}z \nonumber\\ &\leq \iint U_w^i(U_w^{i-1}-U_w^i)\,{\rm d}\kern0.05em y\,{\rm d}z=\iint (U_w^{i-1}-u_s^i)u_s^i\,{\rm d}\kern0.05em y\,{\rm d}z. \end{align}

Generally speaking, $U_w^{i-1}(x,y,z)$ is a three-dimensional spatial-varying velocity distribution, which is influenced by the cumulative effect of $i-1$ wind-turbine wakes and is essentially the base flow for $WT_i$. To further simplify our derivation, we assume the use of $u_b^i(x)$ to replace $U_w^{i-1}$. It should be noted that $u_b^i(x)$ may differ significantly from $U_b(x)$ depending on the wind-farm layout, as it can be influenced by upstream wind-turbine wakes. As a result, the modified momentum deficit can be approximated as follows:

(2.10)\begin{equation} \iint (U_w^{i-1}-u_s^i)u_s^i\,{\rm d}\kern0.05em y\,{\rm d}z \approx \iint [u_b^i(x)-u_s^i]u_s^i\,{\rm d}\kern0.05em y\,{\rm d}z, \end{equation}

where $u_b^i(x)-u_s^i=u_w^i$ is the wake velocity of $WT_i$ accounting for the base flow variation.

Figure 2. Schematic of the control volume for (a) $i$ wind turbines and (b) $i-1$ wind turbines under FPG. (b) is the same as (a) in the absence of $WT_i$.

In our study, $u_b^i(x)$ is modelled as the average wind speed in the rotor area behind $WT_i$, which is sampled in $U_w^{i-1}$ to represent different wake interference conditions, such as full wake and partial wake. To save computational costs, we approximate $u_b^i(x)$ using the average of $N_q=16$ velocity sampling lines downstream of $WT_i$, whose spanwise and vertical coordinates ( $y_{i,q}$ and $z_{i,q}$) are the same as the quadrature points evenly divided over the rotor area of $WT_i$, as shown in figure 3. The coordinates of the quadrature points can be found in Allaerts & Meyers (Reference Allaerts and Meyers2019), which is adopted to calculate the rotor equivalent wind speed. The expression for $u_b^i(x)$ is defined as follows:

(2.11)\begin{equation} u_b^i(x)=\frac{1}{N_q}\sum_{q=1}^{N_q}U_w^{i-1}(x,y_{i,q},z_{i,q}) =U_b(x)-\frac{1}{N_q}\sum_{q=1}^{N_q}U_s^{i-1}(x,y_{i,q},z_{i,q}). \end{equation}

Figure 3. Schematic of (a) quadrature points at the rotor plane of $WT_i$ and (b) velocity sampling lines of $WT_i\ (i\geq 2)$.

For $i=1$, $u_b^1(x)=U_b(x)$, whereas for $i\geq 2$, $u_b^i(x)$ can be obtained through a recursive procedure, which is detailed in Appendix A. It should be noted that this modelling for $u_b^i(x)$ is different from Dar & Porté-Agel (Reference Dar and Porté-Agel2024), who use the wind speed in the wake centreline behind $WT_i$ as $u_b^i(x)$, which may be only applicable for the aligned wind farm layout. In contrast, our proposed modelling for $u_b^i(x)$ can be applicable for arbitrary wind farm layouts.

Finally, the thrust of $WT_i$ can be expressed as

(2.12)\begin{align} \frac{T_i}{\rho}&=\iint u_w^i(x,y,z)[u_b^i(x)-u_w^i(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z\nonumber\\ &\quad +\iint \underbrace{\int_{x_i-\delta x}^{x}\frac{{\rm d}U_b(x^\prime)}{{{\rm d}\kern0.7pt x}^\prime}u_s^i(x^\prime, y,z)\,{{\rm d}\kern0.7pt x}^\prime}_{{\mathfrak p}^i(x,y,z)} {{\rm d}\kern0.05em y} \,{\rm d}z. \end{align}

Here, ${\mathfrak p}^i(x,y,z)$ is the pressure correction term induced by individual velocity deficit.

2.1.3. Linearised expression for the wind-turbine thrust under pressure gradient

It has been shown that the wake velocity deficit $u_s^i$ under pressure gradient has a Gaussian distribution, but with different normalised maximum velocity deficit $C^i(x)$ and wake width $\sigma ^i(x)$ from its counterpart under ZPG (Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a). Therefore, $u_s^i$ and $u_w^i$ for $WT_i$ can be expressed as:

(2.13)$$\begin{gather} u_s^i(x,y,z)=u_b^i(x)C^i(x)\exp\left(-\frac{(y-y_i)^2+(z-z_i)^2}{2\sigma^{i}(x)^2}\right), \end{gather}$$

(2.14)$$\begin{gather}u_w^i(x,y,z)=u_b^i(x)\left[1-C^i(x)\exp\left(-\frac{(y-y_i)^2+(z-z_i)^2}{2\sigma^{i}(x)^2}\right)\right]. \end{gather}$$

We give $C^i(x)$ and $\sigma ^i(x)$ in § 2.2, and for now they are assumed to be known. The first term on the right-hand side of (2.12) can be linearised by introducing a convection velocity $u_c^i(x)$ (Zong & Porté-Agel Reference Zong and Porté-Agel2020), and it can be given analytically as

(2.15)\begin{equation} u_c^i(x)=\frac{\displaystyle\iint u_w^i(x,y,z)[u_b^i(x)-u_w^i(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z}{\displaystyle\iint [u_b^i(x)-u_w^i(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z}=u_b^i(x)\left(1-\frac{C^i(x)}{2}\right). \end{equation}

Therefore, the expression for the thrust of $WT_i$ under pressure gradient can be rewritten as

(2.16)\begin{align} \frac{T_i}{\rho}=\underbrace{u_c^i(x)\iint [u_b^i(x)-u_w^i(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z}_{T_{MD}^i(x)}+\underbrace{\iint \int_{x_i-\delta x}^{x}\frac{{\rm d}U_b(x^\prime)}{{{\rm d}\kern0.7pt x}^\prime}u_s^i(x^\prime, y,z)\,{{\rm d}\kern0.7pt x}^\prime \,{{\rm d} y} \,{\rm d}z}_{T_{\mathfrak p}^i(x)}. \end{align}

This equation indicates that for $WT_i$ operating under pressure gradient, its thrust can be decomposed into two parts: the linearised momentum deficit flux $T_{MD}^i$ and the integral of ${\mathfrak p}^i(x,y,z)$ (i.e. $T_{{\mathfrak p}}^i$), which is related to pressure gradient. In the derivation of (2.12) and (2.16), we have made some assumptions, such as the utilisation of the modified momentum deficit and the substitution of $U_w^{i-1}$ by $u_b^i(x)$. However, these simplifications only lead to negligible errors based on the comparison between the modelled wind turbine thrust and LES results shown in § 4, indicating that these assumptions are reasonable.

2.1.4. PG-IMCM superposition method

Substituting (2.16) into (2.7) yields

(2.17)\begin{align} &\iint U_w(x,y,z)[U_b(x)-U_w(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z \nonumber\\ &\quad= \sum_{i\in B}u_c^i(x)\iint[u_b^i(x)-u_w^i(x,y,z)]\,{\rm d}\kern0.05em y\,{\rm d}z \nonumber\\ &\qquad +\sum_{i\in B}\iint \underbrace{\int_{x_i-\delta x}^{x}\frac{{\rm d}U_b(x^\prime)}{{{\rm d}\kern0.7pt x}^\prime}u_s^i(x^\prime, y,z)\,{{\rm d}\kern0.7pt x}^\prime}_{{\mathfrak p}^i(x,y,z)} \,{\rm d} y\,{\rm d}z \nonumber\\ &\qquad -\iint \underbrace{\int_{x_0}^{x}\frac{{\rm d}U_b(x^\prime)}{{{\rm d}\kern0.7pt x}^\prime}U_s(x^\prime, y,z)\,{{\rm d}\kern0.7pt x}^\prime}_{{\mathfrak P}(x,y,z)} \,{{\rm d}\kern0.05em y} \,{\rm d}z. \end{align}

Analogous to (2.15), we can define the combined convective velocity $U_c(x)$ to linearise the above equation:

(2.18)\begin{equation} U_c(x)=\frac{\displaystyle\iint U_w(x,y,z) U_s(x,y,z)\,{\rm d}\kern0.05em y\,{\rm d}z}{\displaystyle\iint U_s(x,y,z)\,{\rm d}\kern0.05em y\,{\rm d}z}. \end{equation}

Finally, we can obtain the analytical expression for wind-farm velocity deficit $U_s$ under pressure gradient:

(2.19)\begin{equation} \text{PG-IMCM:}\quad U_s(x,y,z)=\underbrace{\sum_{i\in B}\frac{u_c^i(x)}{U_c(x)}u_s^i(x,y,z)}_{U_{s,LWS}}+\underbrace{\sum_{i\in B}\frac{ {\mathfrak p}^i(x,y,z)}{U_c(x)}}_{U_{s,{\mathfrak p}}}- \underbrace{\frac{{\mathfrak P}(x,y,z)}{U_c(x)}}_{U_{s,{\mathfrak P}}}. \end{equation}

It can be seen that the wind-farm flow velocity deficit under pressure gradient consists of three terms. The first term is the linear-weighted sum of individual wake velocity deficits; the second term is the sum of ${\mathfrak p}^i$ divided by $U_c(x)$, whereas the third term is ${\mathfrak P}$ divided by $U_c(x)$. Here, we denote them as $U_{s,LWS}$, $U_{s,{\mathfrak p}}$ and $U_{s,{\mathfrak P}}$, respectively. Since the derivation of this superposition method is motivated by the IMCM method and is applicable to different pressure gradient conditions, we denote it as PG-IMCM here.

Equation (2.19) will have the exact solution of $U_s=u_s^1$ when $N=1$, which is essentially the single-wake model under pressure gradient. In addition, (2.19) will be the same as IMCM (1.6) under the ZPG condition, i.e. ${\rm d}U_b(x)/{{\rm d}\kern0.7pt x}=0$. Therefore, our proposed model can be viewed as a generalisation of the IMCM method under pressure gradient.

Similar to the IMCM method, this model also requires an iterative method (Zong & Porté-Agel Reference Zong and Porté-Agel2020). First, we can give an initial guess of $U_c$ based on the maximum value of $u_c^i$ at different streamwise locations, and $U_s$ can be calculated using (2.19). Then $U_s$ can be substituted into (2.18) to update $U_c$. Finally, we can repeat the procedures until a certain criterion is reached. It should be noted that we compute the integral in the second and third terms by setting the upper limit of $x^\prime$ as $x-\Delta x$, where $\Delta x$ is the grid size for numerical integration. In addition, the lower limit of $x^\prime$ in ${\mathfrak p}^i$ is $x_i-\Delta x$. This approximation will only have negligible error for small values of $\Delta x$. A step-by-step procedure of the PG-IMCM model is given in Appendix A.

2.2. Wind-turbine wake models under pressure gradient

The wake model under pressure gradient has been given as (2.13) and (2.14) in § 2.1.3, but with the maximum normalised velocity deficit $C^i(x)$ and wake width $\sigma ^i(x)$ left to be determined. Here $C^i(x)$ can be solved through an ODE or its asymptotic solution (Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a). In this study, we use the asymptotic solution of the ODE to calculate $C^i(x)$ for its low computational cost and ease of use. It has been shown that, compared with the ODE, the asymptotic solution only has some errors in the near-wake region, whereas it can capture the influence of pressure gradient on the far-wake recovery very well (Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a). We also compared the ODE solution with its asymptotic solution under different atmospheric states in Appendix B.1, and it has been shown that the mean absolute percentage error between these two solutions is within 5 %, supporting that our choice of the asymptotic solution has a comparable prediction accuracy. The wake width $\sigma ^i(x)$ under pressure gradient is determined through an invariant ratio $\lambda _0^i(x)$, which is defined based on the maximum normalised velocity deficit $C_0^i(x)$ and wake width $\sigma _0^i(x)$ under the ZPG condition. The maximum normalised velocity deficit $C^i(x)$ and wake width $\sigma ^i(x)$ of $WT_i$ under pressure gradient can be written as follows (Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018b):

(2.20a,b)\begin{equation} C^i(x)=C_0^i(x)\left(\frac{u_0^i}{u_b^i(x)}\right)^{5/3},\quad \sigma^i(x)=\frac{u_b^i(x)}{\lambda_0^i(x)}C^i(x), \end{equation}

where $u_0^i$ is $u_b^i(x_i)$. In addition, the invariant ratio $\lambda _0^i(x)$ under pressure gradient is given as

(2.21)\begin{equation} \frac{u_0^iC_0^i(x)}{\sigma_0^i(x)}=\lambda_0^i(x)=\frac{u_b^i(x)C^i(x)}{\sigma^i(x)}. \end{equation}

To calculate $C^i(x)$ and $\sigma ^i(x)$, $C_0^i(x)$ and $\lambda _0^i(x)$ must be given at first. Here, we use the well-known Gaussian wake model (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2014) to calculate $C_0^i(x)$ as follows:

(2.22)\begin{equation} C_0^i(x)=\left[1-\sqrt{1-\frac{C_T^i}{8(\sigma_0^i(x)/D)^2}}\right]\mathcal{H}(x-x_i), \end{equation}

where $C_T^i$ is the thrust coefficient of $WT_i$ and the Heaviside step function $\mathcal {H}(x)$ is used here to ensure that wind-turbine wakes only influence the downstream region. To partially characterise the influence of the wind-rotor-induced pressure gradient in the near-wake region, the thrust coefficient can be rewritten in the form of the error function $C_T^i(x)=C_T^i[1+\text {erf}((x-x_i)/D)]/2(x_i\leq x< x_i+2D)$, according to Shapiro, Gayme & Meneveau (Reference Shapiro, Gayme and Meneveau2018) and Zong & Porté-Agel (Reference Zong and Porté-Agel2020). $\sigma _0^i (x)$ is the standard deviation of the Gaussian wake profile of $WT_i$, and its empirical expression is written as follows (Zong & Porté-Agel Reference Zong and Porté-Agel2020):

(2.23)\begin{equation} \frac{\sigma_0^i(x)}{D}=0.35+k_w^i\ln\left[1+\exp\left(\frac{x-x_i-x_{th}^i}{D}\right)\right], \end{equation}

where $k_w^i$ is the wake growth rate of $WT_i$, and it can be determined based on its inflow turbulence intensity $I_u^i$ through the relationship $k_w^i=0.38I_u^i+0.004$ (Niayifar & Porté-Agel Reference Niayifar and Porté-Agel2016). It should be noted that this expression is only applicable to the condition of $0.06< I_u^i<0.15$, $C_T^i\approx 0.8$. However, there is a big chance that $I_u^i>0.15$ on the leeward side. We use $k_w^i=0.26I_u^i$ under this condition (Teng & Markfort Reference Teng and Markfort2020). This empirical expression has been validated by Vahidi & Porté-Agel (Reference Vahidi and Porté-Agel2022) that it has a higher accuracy over $k_w^i=0.38I_u^i+0.004$ when the inflow turbulence intensity is high. In other words, we use a piecewise function to determine $k_w^i$.

The near-wake length $x_{th}^i$ is modelled according to Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2016) and Carbajo Fuertes, Markfort & Porté-Agel (Reference Carbajo Fuertes, Markfort and Porté-Agel2018) under the ZPG condition. Although it has been shown that pressure gradient influences the near-wake length to some extent, the sensitivity analysis presented in Appendix B.2 shows that using the model under the ZPG condition only leads to negligible differences in the calculation of $x_{th}^i$ in comparison with the model proposed by Dar et al. (Reference Dar, Gertler and Porté-Agel2023). Hence, we use the model under the ZPG condition for its low computational cost and simplicity of implementation.

The inflow turbulence intensity $I_u^i$ is modelled as $\sqrt {I_b^2+ (\Delta I_u^i)^2 }$, where $I_b$ is the ambient streamwise turbulence intensity and $\Delta I_u^i$ is the added streamwise turbulence intensity at the rotor plane for $WT_i$, which may be induced by multiple wind-turbine wakes similar to the velocity deficit. Here we define the added turbulence intensity induced by $WT_j$ at the rotor plane of $WT_i$ as $\delta I_u^{j\rightarrow i}$, which can be calculated by the three-dimensional Gaussian model proposed by Ishihara & Qian (Reference Ishihara and Qian2018). To calculate $\Delta I_u^i$ in the merged wake region, we resort to the maximum superposition method expressed as $\Delta I_u^i=\max _{j=1}^{i-1}\delta I_u^{j\rightarrow i}$ (Niayifar & Porté-Agel Reference Niayifar and Porté-Agel2016), for more details about the superposition of the turbulence intensity, please refer to Li et al. (Reference Li, Wang, Ge, Huang, Li and Liu2023). It is noteworthy that the added turbulence intensity model is under the ZPG condition, which is also adopted by Dar & Porté-Agel (Reference Dar and Porté-Agel2024). Although pressure gradient may influence the production of turbulence intensity, developing such a model needs further research.

After obtaining the maximum velocity deficit $C_0^i(x)$ and wake width $\sigma _0^i (x)$ under the ZPG condition, $\lambda _0^i(x)$ can be obtained by (2.21). Substituting $\lambda _0^i(x)$ into (2.20), it is easy to know that

(2.24)\begin{equation} \sigma^i(x)=\sigma_0^i(x)\left(\frac{u_0^i}{u_b^i(x)}\right)^{2/3}. \end{equation}

As shown by (2.20) and (2.24), $C^i(x)$ and $\sigma ^i(x)$ can be expressed directly by its counterpart under the ZPG condition and the base flow variation. Under the FPG condition, $u_0^i$ will be less than $u_b^i(x)$, so both $C^i(x)$ and $\sigma ^i(x)$ will lower than its counterpart under the ZPG condition, indicating a faster wake recovery, and vice versa under the APG condition. By substituting these two variables into (2.13) and (2.14), the wake velocity deficit and wake velocity under pressure gradient can be obtained.

2.3. A simplified form of PG-IMCM

Based on the previous two subsections, we have presented the complete form and implementation details of the PG-IMCM method, which is summarised in Appendix A. What makes the PG-IMCM method difficult to implement and time-consuming to calculate is that the last two terms in (2.19) are integral quantities. Therefore, it is of practical interest for us to investigate the magnitude of these terms and to see whether they can be simplified. For instance, if the sum of the last two terms is small enough compared with $U_{s,LWS}$, they can be neglected, resulting in a much simpler and computationally effective form. However, it is difficult for us to prove the assumption theoretically. Therefore, we first make a simple order-of-magnitude analysis of the sum of the last two terms in (2.19), which can be rewritten as

(2.25)\begin{equation} \Delta U_p=U_{s,{\mathfrak p}}-U_{s,{\mathfrak P}}=\frac{\sum_{i\in B} {\mathfrak p}^i-{\mathfrak P}}{U_c(x)}=\frac{1}{U_c(x)}\int_{x_0}^{x}\frac{{\rm d}U_b(x^\prime)}{{{\rm d}\kern0.7pt x}^\prime}\left(\sum_{i\in B} u_s^i-U_s\right) {{\rm d}\kern0.7pt x}^\prime, \end{equation}

where $\sum _{i\in B} u_s^i$ represents the sum of individual velocity deficits using the local linear superposition method. In fact, the term in parentheses is the difference between the local linear superposition method and the PG-IMCM method, which is denoted as $\Delta U_s$. By definition, it is easy to know that the combined convective velocity $U_c(x)$ has the same order as the base flow $U_b(x)$. If we assume that the base flow has a linear speed-up ratio $\gamma$, i.e. ${\rm d}U_b/{{\rm d}\kern0.7pt x}=\gamma$, we can obtain

(2.26)\begin{equation} \Delta U_p\sim\frac{\gamma (x-x_0)}{U_b(x)}\underbrace{\frac{1}{x-x_0}\int_{x_0}^x \Delta U_s \,{{\rm d}\kern0.7pt x}^\prime}_{\overline{\Delta U_s}}, \end{equation}

where $\overline {\Delta U_s}$ is the streamwise-averaged $\Delta U_s$. In general, the wind speed variation $\gamma (x-x_0)$ within the wind farm caused by pressure gradient has the same order as the base flow $U_b(x)$. Therefore, $\Delta U_p\sim \overline {\Delta U_s}$. Although it has been shown that $\overline {\Delta U_s}\approx 0$ under the ZPG condition (Zong & Porté-Agel Reference Zong and Porté-Agel2020), there is no generalised magnitude estimation applicable to FPG and APG conditions. Consequently, it is infeasible for us to analyse the magnitude of $\Delta U_p$ directly. Moreover, it is incomplete that we only know the magnitude of $\Delta U_p$ since whether we can neglect it depends on its relative magnitude to $U_{s,LWS}$. To this end, we further resort to a model sensitivity analysis under different wind-farm configurations and atmospheric states, which are typical of wind farm operational conditions under pressure gradient, as presented in table 1.

Table 1. Case set-up for the sensitivity study of the PG-IMCM method.

Specifically, we consider the aligned wind farm layout and all wind turbines have the same thrust coefficient $C_T=0.8$. The base flow is modelled as a piecewise function:

(2.27)\begin{equation} U_b(x)=\begin{cases} U_{b0}[1+c(x-x_1)/D],\quad c\geq0,\\ U_{b0}[5+c(x-x_1)/D],\quad c<0, \end{cases} \end{equation}

where $x_1$ is the streamwise location of the first wind turbine, $U_{b0}$ is the inflow velocity at $x_1$ and $c=\gamma D/U_{b0}$ is the normalised speed-up ratio. When $c<0$, the inflow velocity is set to $5U_{b0}$ to prevent $U_{b}(x)<0$ at farther downstream, which is unphysical. By systematically varying the number of wind turbines $N$, streamwise spacing $S_x$, normalised speed-up ratio $c$ and ambient turbulence intensity $I_b$, the mean relative percentage of (2.25) to $U_{s,LWS}$ is quantified, which is defined as the streamwise-averaged absolute ratio of (2.25) to $U_{s,LWS}$ at the wake centreline for each case:

(2.28)\begin{equation} MRP_p=\frac{100\,\%}{x_N+S_x-x_1}\int_{x_1}^{x_N+S_x}\frac{|\Delta U_p|}{U_{s,LWS}}\,{{\rm d}\kern0.7pt x}, \end{equation}

where $x_N$ is the streamwise location of the last wind turbine.

Overall, the averaged $\overline {MRP}_p$ for all the cases is only about 0.4 %, indicating that neglect of the last two terms in the complete form will only lead to negligible error. Moreover, this implies that the local linear superposition method under pressure gradient has a very similar performance to the PG-IMCM method, as discussed later in § 5.2. Based on this analysis, we can provide a simplified form of PG-IMCM as follows:

(2.29)\begin{equation} \text{PG-IMCM-S:}\quad U_s(x,y,z)\approx\sum_{i\in B}\frac{u_c^i(x)}{U_c(x)}u_s^i(x,y,z). \end{equation}

As a simplified form of PG-IMCM, we denote it as PG-IMCM-S, and its calculation procedure is very similar to the PG-IMCM method, which is also given in Appendix A. It should be noted that PG-IMCM-S has the same form as Zong & Porté-Agel (Reference Zong and Porté-Agel2020), but the single turbine wake quantities $u_c^i(x)$ and $u_s^i(x)$ are redefined to account for the influence of pressure gradient, as detailed in § 2.1.3. For the first time, our derivation and analysis reveal the similarities and differences between the momentum-conserving wake superposition methods under different pressure gradient conditions and form the basis for a unified framework. On the other hand, based on the form of PG-IMCM-S, it is demonstrated that wake superposition under pressure gradient does not need to explicitly account for the effect of pressure gradient, whereas it is necessary to account for the influence of pressure gradient on background wind speed and wake recovery when calculating individual velocity deficits. These implications pave the way for extending the empirical superposition principles under pressure gradient in § 5.

Although PG-IMCM and PG-IMCM-S have very close prediction results, their physical meanings differ significantly. PG-IMCM is the direct consequence by linearising the integral momentum-conserving equation, which has a solid theoretical foundation. Specifically, the integral of the product of the wind-farm convection velocity with the sum of the first two terms on the right-hand side of (2.19) is the total thrust of the wind farm, whereas the other term represents the influence of background pressure gradient. In contrast, PG-IMCM-S has limited physical insights and is only an indirect result of the momentum conservation equation. One advantage of the PG-IMCM-S method is that it is easy to implement and computationally effective.

The inputs of the model consist of environment variables and wind turbine variables. The former includes the base flow $U_b(x)$ and ambient streamwise turbulence intensity $I_b$, whereas the latter consists of the rotor centre coordinates $(x_i,y_i,z_i)$, thrust coefficient $C_T^i$ and rotor diameter $D$ of $WT_i$. The output of the model is the wind-farm velocity deficit $U_s(x,y,z)$. In this study, $U_b$ is extracted directly from the main computational domain without wind turbines. For the wind farm layout optimisation in complex terrain, the background flow field can be modelled by using some computational fluid dynamics (CFD) tools, such as WAsP and WindSim, for a number of inflow wind direction sectors (Brogna et al. Reference Brogna, Feng, Sørensen, Shen and Porté-Agel2020). In addition, it can be modelled through some empirical functions (Sun, Yang & Gao Reference Sun, Yang and Gao2023) or linearised perturbation equations (Hunt, Leibovich & Richards Reference Hunt, Leibovich and Richards1988; Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018b) under some simple terrains (such as ramps and hills with gentle slopes).

3.1. LES methodology

Here, we utilise the pseudo-spectral solver LESGO for numerical simulation, which has been validated by wind tunnel experiments (e.g. Stevens, Martínez-Tossas & Meneveau Reference Stevens, Martínez-Tossas and Meneveau2018; Liu & Stevens Reference Liu and Stevens2020b; Zhang et al. Reference Zhang, Huang, Bitsuamlak and Cao2022b) and widely adopted in the studies on wind turbines and wind-farm flows (Du et al. Reference Du, Ge, Zeng, Cui and Liu2021; Ge et al. Reference Ge, Yang, Zhang and Zuo2021; Ma et al. Reference Ma, Ge, Wu, Du and Liu2021; Zhang et al. Reference Zhang, Ge, Liu and Yang2021a; Du, Ge & Liu Reference Du, Ge and Liu2022, among others), atmospheric boundary layer over complex terrain (Liu & Stevens Reference Liu and Stevens2020a,Reference Liu and Stevensb; Zhang et al. Reference Zhang, Huang, Bitsuamlak and Cao2022b) as well as urban canopy flows (e.g. Ge et al. Reference Ge, Zhang, Meng and Ma2020; Fan et al. Reference Fan, Ge, Tan and Li2021; Ge et al. Reference Ge, Yang, Zhang and Zuo2021; Zhang et al. Reference Zhang, Ge, Liu and Yang2021a,Reference Zhang, Yang, Du and Geb, Reference Zhang, Du, Ge and Zuo2022a). LES solves the spatially filtered unsteady incompressible Navier–Stokes equations, where the continuity and momentum equations are written as

(3.1)$$\begin{gather} \frac{\partial \tilde{u}_{i}}{\partial x_{i}} = 0, \end{gather}$$

(3.2)$$\begin{gather}\frac{\partial \tilde{u}_{i}}{\partial t}+\tilde{u}_{j}\frac{\partial \tilde{u}_{i}}{\partial x_{j}} ={-}\frac{\partial \tilde{p}^{*}}{\partial x_{i}}-\frac{\partial \tau^{d}_{ij}}{\partial x_{j}}-\frac{f}{\rho}. \end{gather}$$

Here, $\tilde {{\cdot }}$ is the filter at the spatial scale $\Delta$, $t$ is time, $\tilde {u}_i$ is the filtered velocity in the $i$-direction (with $i=1$, 2 and 3 corresponding to the streamwise ($x$), spanwise ( $y$) and vertical ($z$) directions, respectively), $\tilde {p}^*=\tilde {p}/\rho +1/3\tau _{kk}$ is the modified pressure and $\tau _{ij}^d$ is the deviatoric part of the sub-grid scale stress tensor, which is modelled by the Smagorinsky model with wall damping function (Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2005). Since the Reynolds number based on the rated wind speed ($\sim$10 m s$^{-1}$) and rotor diameter ($\sim$100 m) is ${\sim }10^7$, the effect of the viscous term in the momentum equation is neglected. Here $f$ is the external force, and it can originate from different sources as follows:

(3.3)\begin{equation} f=f_{p}+f_{t}-f_{IB}, \end{equation}

where $f_{p}$ is the constant pressure gradient that drives the flow, $f_{t}$ represents the actuator force induced by the wind turbine and $f_{IB}$ is the body force that needs to be enforced to ensure the grid velocity is zero in the solid region.

Here, we adopt the wall-modelled IBM developed by Liu & Stevens (Reference Liu and Stevens2020b), which is a simplified version of the three-step level-set method proposed by Chester, Meneveau & Parlange (Reference Chester, Meneveau and Parlange2007). This method introduces a signed distance function $\phi (\boldsymbol{x})$, and the computational domain is divided into solid region, fluid region and fluid–solid interface region based on $\phi (\boldsymbol{x})$ at first. Then, different treatments are used for different regions; that is, the velocity at the grid of the solid region is enforced to zero, and the sub-grid scale stress of the fluid–solid interface region is calculated using the wall model, whereas no further treatment is done to the grid in the fluid region. For more details about the wall-modelled IBM, please refer to Liu & Stevens (Reference Liu and Stevens2020b).

Wind turbines are modelled using the filtered actuator disk model (Shapiro, Gayme & Meneveau Reference Shapiro, Gayme and Meneveau2019). Based on the one-dimensional momentum theory, the total thrust $F_t$ can be expressed as

(3.4)\begin{equation} F_t=\frac{1}{2}\rho AC_TU_{\infty}^2=\frac{1}{2}\rho AC_T\frac{U_d^2}{(1-a)^2}, \end{equation}

where $U_{\infty }$ is the freestream inflow wind speed, $A$ is the swept area of the wind rotor, $\rho$ is the air density, $C_T$ is the thrust coefficient, $a$ is the axial induction factor and $U_d$ is the rotor-averaged velocity. Finally, the total thrust $F_t$ is distributed to the computation grids to obtain $f_{t}$ according to the Gaussian-filtered indicator function $\Re (\boldsymbol {x})$. For more details about the determination of $\Re (\boldsymbol {x})$, please refer to Shapiro et al. (Reference Shapiro, Gayme and Meneveau2019).

Uniform grids are used in the horizontal direction, and a staggered Cartesian grid is used in the vertical direction. The pseudo-spectral method is used for the horizontal spatial discretisation. In the vertical direction, the second-order central finite difference method is adopted. For time advancement, the code uses the second-order accurate Adams–Bashforth method.

Note that the filter symbol is omitted hereafter for convenience. Although the pressure simulated in this study includes the trace of the sub-grid scale stress tensor, it is expected only to cause a small difference in the momentum balance and the Bernoulli equation (Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021).

3.2. Case set-up

In general, convex (concave) topography will induce FPG (APG). To validate the performance of the proposed WFFAMs under pressure gradient, we set up two types of cases. One type is the ideal convex and concave channels consisting of symmetric upper and lower walls with respect to the centre plane of the channel (Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a), which can be used to simulate FPG and APG conditions, respectively. The surface equation $Z_l(x,y)$ of the lower wall has the following general cosine form:

(3.5)\begin{equation} Z_l(x,y)=\begin{cases} \displaystyle\frac{1}{2}H\left[1+\cos\left(\frac{\rm \pi}{2L_u}(x-x_{ct})\right)\right], & -2L_u\leq x-x_{ct}<0,\\ \displaystyle\frac{1}{2}H\left[1+\cos\left(\frac{\rm \pi}{2L_d}(x-x_{ct})\right)\right], & 0\leq x-x_{ct}<2L_d,\\ 0, & \text{otherwise}, \end{cases} \end{equation}

where $x_{ct}$ is the position of the channel throat, $H$ is the height difference between the throat position and the horizontal ground, and $L_u$ and $L_d$ are the upstream and downstream half-lengths of the cosine functions, respectively. The corresponding expression for the upper wall is $Z_u (x,y)=10D-Z_l(x,y)$. The detailed parameter settings are listed in table 2. In this condition, wind turbines are placed at the centre of the $y$–$z$ plane of the computational domain, as shown in figure 4.

Table 2. Overview of the ideal convex and concave channel cases. (CVW and CCW represent the convex and concave walls, respectively, and WT represents there are wind turbines in this case; $x_T$ is the streamwise location of the rotor centre.)

Figure 4. Schematics of the computational domain for (a) FPG-CVW-WT and (b) APG-CCW-WT.

Another type of case consists of two linear ramps and one platform at the bottom (Dar et al. Reference Dar, Gertler and Porté-Agel2023; Dar & Porté-Agel Reference Dar and Porté-Agel2024), as shown in figure 5. Wind turbines are placed at the same distance over the ground, and their rotors are perpendicular to the ramp. The hub height is the same as the rotor diameter. To facilitate the comparison between the proposed model and the ramp case, we define a local ramp coordinate following the slope of the ramp, whose origin is located in the rotor centre of the first wind turbine in the wind farm array, as shown in figure 5. This case setting is closer to the realistic engineering situation than channel cases. The analytical expression of the ramp height $Z_r(x,y)$ is written as follows:

(3.6)\begin{align} Z_r(x,y)=\begin{cases} H_r\left(\dfrac{x-x_w}{L_1}\right), & x_w\leq x< x_w+L_1,\\ H_r, & x_w+L_1\leq x< x_w+L_1+L_2,\\ H_r\left(1-\dfrac{x-x_w-L_1-L_2}{L_3}\right), & x_w+L_1+L_2\leq x< x_w+L_1+L_2+L_3,\\ 0, & \text{otherwise}, \end{cases} \end{align}

where $x_{w}$ is the starting location of the windward ramp, $H_r$ is the height of the platform, $L_2$ is the length of the platform, and $L_1$ and $L_3$ are the projected lengths of the windward and leeward ramps, respectively. The detailed parameter settings are given in table 3.

Figure 5. Schematics of the computational domain for (a) FPG-RAMP-WT and (b) APG-RAMP-WT.

Table 3. Overview of the ramp cases.

It should be noted that although we made no assumptions about the wind-farm layout in the model derivation, we only simulated the aligned layout in the validation cases, which serves as a basis for verifying the wake superposition method (Zong & Porté-Agel Reference Zong and Porté-Agel2020; Lanzilao & Meyers Reference Lanzilao and Meyers2022; Dar & Porté-Agel Reference Dar and Porté-Agel2024). In this case, the wake generated by different wind turbines will be fully merged, and the wake superposition is very complicated, providing a fair condition to test the performance of different superposition methods.

The main computational domain for the FPG-CVW-WT, FPG-RAMP-WT and APG-RAMP-WT cases is set to $80D\times 12D\times 10D$. However, due to the setting of concave walls, the main computational domain for the APG-CCW-WT case is expanded to $80D\times 12D\times 17D$. The number of grids in both computational domains is $320\times 120\times 144$, which can ensure the rotor diameter is covered with about 10 points in both the spanwise and vertical directions. Such a grid resolution is fine enough for the actuator disk model to predict the velocity deficit in turbine wakes with reasonable accuracy (Wu & Porté-Agel Reference Wu and Porté-Agel2011; Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a,Reference Shamsoddin and Porté-Agelb; Lin & Porté-Agel Reference Lin and Porté-Agel2022). The thrust coefficient of all wind turbines is set to 0.8. In addition to the above four cases with wind turbines, four cases without wind turbines under the same condition are simulated in tables 2 and 3 to obtain the base flow, which is used as inputs to WFFAMs and to calculate the normalised velocity deficit. To overcome the streamwise periodicity and ensure a fully developed inflow to the main computational domain, we use the concurrent precursor method proposed by Stevens, Graham & Meneveau (Reference Stevens, Graham and Meneveau2014), in which the flow field in the fringe region at the end of the main domain is blended with the ideal straight channel flow for the curved channel cases and half-channel flow for the ramp cases, whose domain size and grid size are the same as the main computational domain. The computational accuracy of the LES framework in these channel and half-channel flows has been tested in several validation studies (e.g. Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a,Reference Shamsoddin and Porté-Agelb; Liu & Stevens Reference Liu and Stevens2020b; Zhang et al. Reference Zhang, Huang, Bitsuamlak and Cao2022b). The surface roughness of the ground and topography are both $0.001D$, a typical value of onshore topography.

4.1. Validation of the momentum-conserving method under FPG

In this subsection, we validate the performance of the momentum-conserving method under FPG. Under FPG condition, the pressure decreases streamwise, and the background wind velocity increases, which is beneficial to the wake recovery (Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a).

4.1.1. Validation of results based on the FPG-CVW-WT case

Figure 6(a) shows the streamwise variations of the normalised spanwise-averaged velocity and pressure at hub height for the FPG-CVW case. It can be seen that the pressure decreases as the streamwise distance increases, and the wind speed increases gradually, indicating that we have successfully simulated the ideal FPG case. The mean base flow gradient $\overline {{\rm d}U_b/{{\rm d}\kern0.7pt x}}$ from $16D$ to $50D$ is about $0.0224U_{b0}/D$. In addition, figure 6(b) shows the streamwise variations of the normalised spanwise-averaged kinetic energy gradient and pressure gradient at hub height for the FPG-CVW case. We can see that these two curves nearly coincide, indicating that the pressure and wind speed satisfy the Bernoulli equation $\langle P_{bg}\rangle _y/\rho +1/2 \langle U_{bg} \rangle _y^2=\text {const.}$. Therefore, this case meets the applicability of the wake model and PG-IMCM superposition method under pressure gradient, which is necessary for a fair comparison.

Figure 6. Streamwise variations of the normalised spanwise-averaged (a) velocity and pressure and (b) kinetic energy gradient and pressure gradient at hub height for the FPG-CVW case. (Vertical black dashed lines represent the locations of wind turbines.)

Figure 7 presents the contours of the normalised velocity deficit at the $x$–$z$ rotor symmetry plane. It can be seen that the velocity deficit downstream of the first wind turbine is the smallest, and its recovery is the slowest, whereas the wake velocity deficit downstream of the last wind turbine is the largest. In addition, the second and third wind turbines have similar wake velocity deficit and recovery rates.

Figure 7. Contours of the normalised velocity deficit at the $x$–$z$ rotor symmetry plane for the FPG-CVW-WT case.

Whether the maximum velocity deficit at hub height $U_{s,max}$ in the merged wake region can be accurately predicted is one of the crucial metrics for the assessment of superposition methods. Figure 8 shows the comparison of $U_{s,max}/U_{b0}$ obtained by PG-IMCM and LES for the FPG-CVW-WT case. In addition to the total velocity deficit, its decomposition terms (i.e. $U_{s,LWS}$, $U_{s,{\mathfrak p}}$ and $-U_{s,{\mathfrak P}}$) are also presented in the figure. As expected, the sum of the pressure correction term induced by individual velocity deficits $\sum _{i\in B} {\mathfrak p}^i$ and the pressure correction term induced by total velocity deficit $-{\mathfrak P}$ accumulate with the downstream distance, as they are integral quantities. When divided by the combined velocity $U_c$, they become $U_{s,{\mathfrak p}}$ and $-U_{s,{\mathfrak P}}$, whose variation is very similar to the sum of four stair-step functions located at the wind-turbine locations. As expected, they have opposite symbols but very close magnitudes, and the mean absolute value of their sum is only about $0.008U_{b0}$. Therefore, the sum of these two terms can be neglected, and the prediction result of the PG-IMCM method can be viewed as the linear-weighted sum of individual velocity deficits. In addition, we also present the prediction results of PG-IMCM-S, which is virtually identical with PG-IMCM, further validating our simplification based on the sensitivity study. Overall, PG-IMCM and its simplified form have a good prediction over $U_{s,max}$, only with a slight overestimation in the merged wake region.

Figure 8. Model prediction of $U_{s,max}/U_{b0}$ in comparison with LES for the FPG-CVW-WT case. (Vertical black dashed lines represent the locations of wind turbines.)

Furthermore, we are also concerned about the prediction of velocity deficit profiles by WFFAMs in practical applications since it is closely related to the power calculation of wind turbines. Figure 9 compares the spanwise velocity deficit profiles at hub height obtained by PG-IMCM, PG-IMCM-S and LES downstream of each wind turbine for the FPG-CVW-WT case. Similar to figure 8, a term-by-term comparison is also presented. At first, we can see that $U_{s,{\mathfrak p}}$ and $-U_{s,{\mathfrak P}}$ have a near Gaussian distribution, and their sum can be neglected. The prediction of PG-IMCM and its simplified form is nearly collapsed at each radial location, similar to the rotor centre. Overall, the PG-IMCM method and its simplified form have a good prediction accuracy over the profiles, with a slight overestimation in the rotor centre region.

Figure 9. Model prediction of the spanwise velocity-deficit profiles at hub height downstream of each wind turbine in comparison with LES for the FPG-CVW-WT case.

LES results can not only serve as a baseline for the model validation but also can be used to examine some assumptions in the derivation of the PG-IMCM method. Here, we compare the wind turbine thrust modelling (2.16) based on the single-wake quantities with LES results, as shown in figure 10. It can be seen that the wind thrust modelling can represent the variation and magnitude of the total thrust in the control volume $B$ very well. The sum of momentum deficits of stand-alone turbine wakes $\sum _{i\in B} T_{MD}^i$ under FPG condition is less than the total thrust, whereas the sum of pressure correction terms $\sum _{i\in B} T_{{\mathfrak p}}^i$ makes up for this deficit, finally resulting in a good prediction of the total thrust, which is an essential part of the PG-IMCM method. It should be noted that the slight underestimation in the near-wake region is related to the modified thrust coefficient within the downstream $2D$ region.

Figure 10. Comparison of the modelled wind-turbine thrust with LES for the FPG-CVW-WT case. (Vertical black dashed lines represent the locations of wind turbines.)

4.1.2. Validation of results based on the FPG-RAMP-WT case

The FPG-CVW-WT case is too ideal, leading to minimal wind shear and ambient streamwise turbulence intensity (about 2.3 %) in the rotor region, which does not represent the realistic atmospheric surface boundary layer with high wind shear and high turbulence intensity. Therefore, to further validate the applicability of the proposed momentum-conserving method under FPG condition, we set up the FPG-RAMP-WT case, as depicted in figure 5(a), which is closer to the realistic operation environment of the wind turbine. In this case, streamwise turbulence intensity at hub height is about 6 %. Because of the ramp, the flow has a streamwise component $U_x$ and a vertical component $U_z$. Therefore, the velocity $U$ used in the analysis of the FPG-RAMP-WT and APG-RAMP-WT cases is $\sqrt {U_x^2+U_z^2}$.

Figure 11 shows the contours of the spanwise-averaged wind field and pressure for the FPG-RAMP case. It can be observed that the flow up the windward ramp evolves downstream approximately parallel to the surface and gradually accelerates before approaching the platform. The corresponding pressure gradually decreases, indicating that the FPG condition is successfully simulated through this set-up. The base flow has a nearly constant speed-up ratio above the windward ramp, and its mean value $\overline {{\rm d}U_b/{\rm d}{\kern0.8pt}x_r}$ from $0D$ to $26D$ is about $0.0215U_{b0}/D$. We also examined the kinetic energy gradient and pressure gradient at hub height for the FPG-RAMP case, and they still approximately conform to the Bernoulli equation. Here, we do not show them for brevity.

Figure 11. Contours of the spanwise-averaged (a) wind field and (b) pressure for the FPG-RAMP case. (The white dashed line represents hub height, and black solid lines represent streamlines.)

Figure 12 shows the contours of the normalised velocity deficit and streamwise turbulence intensity at the $x$–$z$ rotor symmetry plane for the FPG-RAMP-WT case. We can see that the velocity deficit computed by LES has an obvious accumulation pattern in this case; that is, the wake velocity deficit of the first wind turbine is the smallest, whereas the wake velocity deficit downstream of the third wind turbine is the largest. As for the streamwise turbulence intensity, its magnitude in the wake of the first wind turbine is low, whereas its magnitude downstream of the second and third wind turbines is high, and their spatial distribution is quite similar. This finding is consistent with Li et al. (Reference Li, Wang, Ge, Huang, Li and Liu2023), who have studied the spatial distribution characteristics of the added streamwise turbulence intensity in the merged wake region of an aligned wind turbine array under the ZPG condition, indicating that the adoption of the maximum superposition method for the wake-added turbulence is reasonable.

Figure 12. Contours of the (a) normalised velocity deficit and (b) streamwise turbulence intensity at the $x$–$z$ rotor symmetry plane for the FPG-RAMP-WT case. (The white dashed line represents hub height.)

Figures 13 and 14 show the comparison of $U_{s,max}/U_{b0}$ and the spanwise velocity deficit profiles at hub height obtained by the momentum-conserving WFFAMs and LES for the FPG-RAMP-WT case, respectively. The term-by-term comparison results are very similar to the FPG-CCW-WT case, thus not repeated here. Overall, the momentum-conserving methods show a remarkably high prediction accuracy in this case, further validating the accuracy of the proposed superposition method. A comparison of the modelled turbine thrust with LES results is further shown in figure 15. Similar to the results presented in figure 10, we can observe that the modelled turbine thrust agrees well with LES results, further confirming the validity of (2.16).

Figure 13. Same as figure 8, but for the FPG-RAMP-WT case.

Figure 14. Same as figure 9, but for the FPG-RAMP-WT case.

Figure 15. Same as figure 10, but for the FPG-RAMP-WT case.

The term-by-term analysis of the proposed momentum-conserving method in the above two cases with different ambient turbulence intensity under FPG conditions confirms the validity of the linearised expression of the turbine thrust derived in § 2.1.3 and shows that our proposed model can predict the maximum velocity deficit and spanwise velocity deficit profiles at hub height very well.

Interestingly, the proposed methods have higher accuracy in the realistic ramp case than the ideal concave channel case. This may be attributed to two reasons. On the one hand, the ambient streamwise turbulence intensity is only 2.3 % in the ideal concave channel case, resulting in a much longer near-wake region. The velocity deficit profiles of the first turbine are closer to super-Gaussian, as depicted in figure 9(a). When we use the Gaussian profile, it will result in a more slender spanwise profile. On the other hand, the ideal curved channel cases may be out of the application range of some empirical relation of wake quantities, especially the wake growth rate. This is because this expression is obtained by curve-fitting based on the LES results and field measurements in the atmospheric boundary layer, which is characterised by the high wind shear and ambient turbulence intensity.

4.2. Validation of the momentum-conserving method under APG

In this subsection, we focus on the validation of the momentum-conserving method under APG. In contrast to the FPG condition, the pressure increases streamwise and the corresponding wind velocity decreases under the APG condition. Compared with the ZPG condition, the APG condition will make the wake recovery slower (Shamsoddin & Porté-Agel Reference Shamsoddin and Porté-Agel2018a).

4.2.1. Validation of results based on the APG-CCW-WT case

The base flow has a nearly constant gradient in the region of interest, and the mean base flow gradient $\overline {{\rm d}U_b/{{\rm d}\kern0.7pt x}}$ from $16D$ to $50D$ is about $-0.0062U_{b0}/D$ for the APG-CCW case. The ambient turbulence intensity is about 2.5 %. In addition, we compare the normalised spanwise-averaged kinetic energy gradient and pressure gradient at hub height for the APG-CCW case, and they almost collapse, indicating that the case set-up satisfies the applicability of the wake model and the PG-IMCM superposition method. For brevity, we do not show them here. Figure 16 shows the contours of the normalised velocity deficit at the $x$–$z$ rotor symmetry plane for the APG-CCW-WT case. It can be seen that the wake recovery is the slowest downstream of the first wind turbine. The wake velocity deficit downstream of the second wind turbine is the largest, and then it recovers very quickly. The wake velocity deficits and their recovery downstream of the third and fourth wind turbines are relatively similar. In this case, the spatial distribution of the merged wake velocity deficit is significantly different from that under ZPG and FPG conditions. The proposed momentum-conserving method cannot accurately predict the variations in the maximum velocity deficit in the merged wake region in this case, as shown in figure 17.

Figure 16. Same as figure 7, but for the APG-CCW-WT case.

Figure 17. Model prediction of $U_{s,max}/U_{b0}$ in comparison with LES for the APG-CCW-WT case.

One possible explanation for this phenomenon is the combined effects of the low ambient turbulence intensity and APG, which stimulates the wake meandering behind the second wind turbine. Observations from various studies, such as Gupta & Wan (Reference Gupta and Wan2019), Li, Dong & Yang (Reference Li, Dong and Yang2022), Gambuzza & Ganapathisubramani (Reference Gambuzza and Ganapathisubramani2023) and Messmer, Hölling & Peinke (Reference Messmer, Hölling and Peinke2024), have supported the mechanism that wake meandering occurs in the far wake due to the amplification of upstream disturbances by shear flow instabilities, especially when the inflow's turbulence intensity is low. Moreover, the optimal perturbations are those large-scale (the same order as the turbine rotor) and low-frequency coherent eddies (Mao & Sørensen Reference Mao and Sørensen2018). This is the situation for the second wind turbine, which is immersed in the wake of the first wind turbine. In this condition, the small perturbations in the far-wake region of the first turbine will be amplified significantly by the complex nonlinear dynamics, leading to enhanced wake recovery. As a result, the wake downstream of the second wind turbine recovers quickly and is in an unsteady physical process, which does not satisfy the applicability of existing wake models and their empirical expressions. A similar phenomenon has also been observed in the study on the influence of swells on the wake evolution of aligned wind turbines (Yang et al. Reference Yang, Ge, Abkar and Yang2022).

Although the same situation is also satisfied in the FPG-CVW-WT case, our model only has a slight underestimation of wake recovery, which may be related to the potential influence of pressure gradient on the wake meandering. In other words, the FPG condition will suppress the selective amplification, making the empirical expression of wake growth rate still applicable only with a minor error, whereas the APG condition will enhance this process, resulting in a high prediction error.

It should be noted that our explanation is a qualitative conjecture based on the simulation results. A quantitative analysis on the influence of pressure gradient on wake meandering, especially the selective amplification mechanism, is beyond the scope of the study, which is an interesting topic for future research.

Therefore, the misestimation of PG-IMCM in this case can be attributed to the invalidity of the stand-alone turbine wake model rather than the superposition method. No further comparison based on this case is given hereafter for brevity.

4.2.2. Validation of results based on the APG-RAMP-WT case

To solve the problems mentioned previously, we further set up a leeward ramp to simulate the APG condition with a high ambient turbulence intensity, as depicted in figure 5(b). The prediction results of the momentum-conserving method in this case are further compared and analysed.

At first, figure 18 shows the contours of the spanwise-averaged wind field and pressure for the APG-RAMP case. We can observe that the flow upper the leeward ramp evolves downstream nearly parallel to the surface and gradually decelerates before approaching the ground, and the corresponding pressure gradually increases. It indicates that the expected APG condition is successfully simulated. The base flow has a nearly constant gradient above the leeward ramp, and its mean value $\overline {{\rm d}U_b/{\rm d}{\kern0.8pt}x_r}$ from $0D$ to $26D$ is about $-0.0128U_{b0}/D$. Furthermore, we compare the kinetic energy gradient and pressure gradient at hub height for the APG-RAMP case. Similar to the FPG-RAMP case, they approximately satisfy the Bernoulli equation in the region of interest and meet the prerequisite of the wake model and the PG-IMCM superposition model under pressure gradient. For brevity, we do not show them here.

Figure 18. Same as figure 11, but for the APG-RAMP case.

The contours of the normalised velocity deficit and streamwise turbulence intensity at the $x$–$z$ rotor symmetry plane for the APG-RAMP-WT case are depicted in figure 19. In the upper region of the leeward ramp, the downslope flow affects the wind turbine, and the streamwise turbulence intensity at the location of $WT_i$ is about 16 %. As a result, wind-turbine wakes recover very quickly. Figures 20 and 21 show the comparison of $U_{s,max}/U_{b0}$ and the spanwise velocity deficit profiles at hub height obtained by the momentum-conserving WFFAMs and LES for the APG-RAMP-WT case, respectively. Due to ${\rm d}U_b/{\rm d}{\kern0.8pt}x_r<0$ under APG condition, $U_{s,{\mathfrak p}}$ is less than zero, whereas $-U_{s,{\mathfrak P}}$ is larger than zero, which is in contrast to the case under the FPG condition. These two terms have a relatively close absolute value, with the mean absolute value of their sum only about $0.004U_{b0}$. Hence, neglecting them only leads to negligible differences throughout the region of interest. The momentum-conserving method performs well in this case, validating the applicability of PG-IMCM and its simplified form under APG condition.

Figure 19. Same as figure 12, but for the APG-RAMP-WT case.

Figure 20. Same as figure 8, but for the APG-RAMP-WT case.

Figure 21. Same as figure 9, but for the APG-RAMP-WT case.

Figure 22 compares the cumulative turbine thrust modelled by (2.16) with LES results in this case. Different from the FPG condition, $\sum _{i\in B}T_{MD}^i$ under APG condition is larger than the thrust, whereas the pressure correction term $\sum _{i\in B} T_{\mathfrak p}^i$ cuts down this excess, finally leading to a satisfactory prediction of the cumulative turbine thrust.

Figure 22. Same as figure 10, but for the APG-RAMP-WT case.

In summary, the proposed model has a satisfactory performance in this APG condition. It is worth mentioning that the APG condition usually occurs in the high-turbulence leeward region in the atmospheric boundary layer, which is far from the ideal case studied in § 4.2.1. Therefore, although the PG-IMCM model has a poor performance in the ideal APG-CCW-WT case, which is related to the inapplicability of empirical expression for individual velocity deficits, it will not influence its engineering applications in realistic APG conditions similar to the APG-RAMP-WT case, in which it still has a high prediction accuracy.

5.1. WFFAMs under pressure gradient based on empirical superposition principles

Specifically, we can develop five empirical WFFAMs based on the empirical superposition methods (i.e. global linear, global square, local linear, local square and wind product) introduced in § 1.2 and the single-wake model under pressure gradient given in § 2.2. In addition, we have shown that these empirical WFFAMs may have a good performance provided that they can account for the base flow variations based on our derivation, which is also supported by Dar & Porté-Agel (Reference Dar and Porté-Agel2024). Therefore, their forms have been changed as follows:

(5.1)\begin{align} &\text{PG-GL:}\quad U_s(x,y,z)=\sum_i(U_b(x)-u_w^i(x,y,z)), \end{align}

(5.2)\begin{align} &\text{PG-GS:}\quad U_s(x,y,z)=\sqrt{\sum_i(U_b(x)-u_w^i(x,y,z))^2}, \end{align}

(5.3)\begin{align} &\text{PG-LL:}\quad U_s(x,y,z)=\sum_i( u_b^i(x)-u_w^i(x,y,z)), \end{align}

(5.4)\begin{align} &\text{PG-LS:}\quad U_s(x,y,z)=\sqrt{\sum_i( u_b^i(x)-u_w^i(x,y,z))^2}, \end{align}

(5.5)\begin{align} &\text{PG-WP:}\quad U_s(x,y,z)=U_b(x)-U_b(x)\prod_i\left(\frac{u_w^i(x,y,z)}{u_b^i(x)}\right). \end{align}

Compared with their counterparts under the ZPG condition, there are three main differences: (i) $U_\infty$ is replaced by $U_b(x)$; (ii) $u_0^i$ is replaced by $u_b^i(x)$, which is the base flow for $WT_i$ defined in § 2.1.2; (iii) the calculation of $u_w^i(x,y,z)$ will account for the influence of pressure gradient. To be consistent with the different physical meanings between global and linear under the ZPG condition, the base flow and turbulence intensity used to calculate $u_w^i$ are different. In the global superposition methods (i.e. PG-GL and PG-GS), $u_w^i$ is calculated based on $U_b(x)$ and ambient turbulence intensity $I_b$, whereas $u_w^i$ is calculated based on $u_b^i(x)$ and inflow turbulence intensity $I_u^i$ in the local superposition methods (i.e. PG-LL and PG-LS) and PG-WP.

In total, we have developed five empirical WFFAMs under pressure gradient, and their classification and details are given in table 4. They have the same inputs and output as the momentum-conserving model. The step-by-step procedures of these empirical WFFAMs are detailed in Appendix C.

Table 4. Classification and details of the proposed empirical WFFAMs under pressure gradient.

5.2. Validation results

In this subsection, we inter-compare the prediction results of these five empirical WFFAMs in the FPG-CVW-WT, FPG-RAMP-WT and APG-RAMP-WT cases. Figure 23 shows the comparison of the $U_{s,max}/U_{b0}$ obtained by the empirical WFFAMs and PG-IMCM with LES for the three validation cases. We can see that the prediction of empirical WFFAMs under pressure gradient have many similarities with their counterparts under the ZPG condition. Overall, the PG-GL model will significantly overestimate the wind farm velocity deficit $U_s$ in the merged wake region. This is likely due to the fact that the background velocity, which does not account for the wake effect of upstream wind turbines, is used as the base flow, leading to the overestimation of individual wake velocity deficits. Moreover, when they are summed directly as the PG-GL model, it will overestimate $U_s$. If they are superimposed as the PG-GS model, it seems that only the largest individual wake velocity deficit will be left, and this may alleviate the overestimation of $U_s$. However, the PG-GS method still fails to have a satisfactory prediction, especially when the ambient turbulence intensity is low to moderate. If we use the base flow updated by accounting for upstream turbine wakes, the individual velocity deficits may be predicted accurately. However, when they are superimposed as the PG-LS model, it may suffer from the same problems as the PG-GS model, leading to the underestimation of the combined velocity deficit. Overall, the PG-LL model has a satisfactory performance. These findings are consistent with the inter-comparison results under the ZPG condition (Niayifar & Porté-Agel Reference Niayifar and Porté-Agel2016; Zong & Porté-Agel Reference Zong and Porté-Agel2020; Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021). In addition, the prediction results of the PG-WP model are almost identical to those of the PG-LL model, which is also observed by Lanzilao & Meyers (Reference Lanzilao and Meyers2022) under the ZPG condition.

Figure 23. Comparison of $U_{s,max}/U_{b0}$ obtained by different WFFAMs with LES for the cases: (a) FPG-CVW-WT, (b) FPG-RAMP-WT and (c) APG-RAMP-WT. Vertical black dashed lines represent the locations of wind turbines.

Figure 24. Comparison of (a) the spanwise velocity deficit profiles at hub height and (b) vertical velocity deficit profiles at the $x$–$z$ rotor symmetry plane obtained by different WFFAMs with LES downstream of each wind turbine for the FPG-CVW-WT case.

Figure 25. Comparison of (a) the spanwise velocity deficit profiles at hub height and (b) vertical velocity deficit profiles at the $x_r$–$z_r$ rotor symmetry plane obtained by different WFFAMs with LES downstream of each wind turbine for the FPG-RAMP-WT case. Black dot-dashed lines represent the hub height.

Figure 26. Same as figure 25, but for the APG-RAMP-WT case.

Interestingly, the performance of PG-LL and PG-WP is very close to the momentum-conserving method, with only discernable differences in the near-wake region of the merged wake, which has also been observed in Zong & Porté-Agel (Reference Zong and Porté-Agel2020). This is related to the gradual variation of the mean convection velocity $u_c$ and the sudden change in the combined convection velocity $U_c$ in the near-wake region of downstream turbines, leading to the weights for individual velocity deficits in these regions deviating from 1, which is exactly the weight of the PG-LL model. In contrast, the weight for individual velocity deficits in the PG-IMCM method in the far-wake region of downstream turbines will approach 1, resulting in a nearly collapsed prediction. A direct comparison between PG-WP and PG-IMCM is more intricate, as they have quite different forms. However, we believe that its high prediction accuracy lies in its intrinsic ability to update the base flow for downstream wind turbines, which is supported by Dar & Porté-Agel (Reference Dar and Porté-Agel2024).

Figures 24–26 further show comparisons of the spanwise and vertical velocity deficit profiles obtained by different WFFAMs against LES results for FPG-CVW-WT, FPG-RAMP-WT and APG-RAMP-WT case, respectively. Overall, the PG-IMCM, PG-LL and PG-WP models have the best performance. As expected, the spanwise and vertical profiles are almost identical in the convex channel case since it is free from the influence of the wall. In contrast, affected by the windward or leeward ramp, the wake velocity deficit has asymmetric characteristics. As a result, the spanwise and vertical wake widths have some differences, especially in the APG-RAMP-WT case. In addition, the wake centre position is slightly shifted downwards with respect to the hub height in the FPG-RAMP-WT case, as shown in figure 12(a), whereas the wake centre shifts upwards in the APG-RAMP-WT case, as depicted in figure 19(a). These pose significant challenges to accurately predicting the spatial distribution of wake velocity deficit and are difficult to model. Nevertheless, the adoption of the axisymmetric wake model and the same spanwise and vertical wake growth rates in the WFFAMs still have satisfactory performance, as shown in figures 25 and 26(b). The difference in the wake width prediction of the APG-RAMP-WT case can be alleviated to some extent by lowering the spanwise wake growth rate slightly and increasing the vertical wake growth rate accordingly to keep the equivalent wake width unchanged. However, developing a general empirical relationship of wake growth rates in different directions is beyond the scope of this study, and it needs to be better understood and addressed in future research.

Overall, the inter-comparison results of the empirical WFFAMs under pressure gradient share many similarities with the inter-comparison results under the ZPG condition. Among the five considered empirical models, the PG-LL and PG-WP models have a satisfactory and consistent performance with the PG-IMCM model in all validation cases, especially in realistic ramp cases. The secret behind their good performance can be attributed to the ability to properly update the base flow for the downstream wind turbines as well as reasonable superposition principles.