# Ainslie Wake Model with Forcing Term

Norbert G W Warncke, Julie Paddeu, Nathan Mutin, Etienne Muller, Felix

Houtin Mongrolle, Pierre Bénard

► To cite this version:

Norbert G W Warncke, Julie Paddeu, Nathan Mutin, Etienne Muller, Felix Houtin Mongrolle, et al. Ainslie Wake Model with Forcing Term. Wake Conference 2025, Jun 2025, Uppsala, Sweden. pp.012016, 10.1088/1742-6596/3016/1/012016. [hal-05491734](https://hal.science/hal-05491734v1)

HAL Id: hal-05491734

https://hal.science/hal-05491734v1

Submitted on 3 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

PAPER • OPEN ACCESS

# Ainslie Wake Model with Forcing Term

To cite this article: Norbert G.W. Warncke et al 2025 *J. Phys.: Conf. Ser.* **3016** 012016

View the [article online](#) for updates and enhancements.

## You may also like

- Modelling higher trapped fields by pulsed field magnetisation of composite bulk MgB₂ superconducting rings
V Cientanni, M D Ainslie, H Fujishiro et al.

- Analysis of critical current inhomogeneity in r-z plane of GdBCO superconducting bulk and simulation of flux jumps during pulsed field magnetization
J T Hu, W J Yang, D F Zhou et al.

- Numerical simulation of flux jump behavior in REBaCuO ring bulks with an inhomogeneous J_c profile during pulsed-field magnetization
Tatsuya Hirano, Hiroyuki Fujishiro, Tomoyuki Naito et al.
# Ainslie Wake Model with Forcing Term

Norbert G.W. Warncke¹, Julie Paddeu¹, Nathan Mutin¹, Etienne Muller², Felix Houtin Mongrolle¹, Pierre Bénard²

¹Siemens Gamesa Renewable Energy B.V., Prinses Beatrixlaan 800, 2595 BN Den Haag, NL

²Univ Rouen Normandie, INSA Rouen Normandie, CNRS, CORIA UMR6614, 675, avenue de l'Université, Saint-Etienne-du-Rouvray, 76801, France

E-mail: norbert.warncke@siemensgamesa.com

**Abstract.** The Ainslie Wake model, as part of the Dynamic Wake Meandering model, is the industry standard wake model for wind farm loads. Its defining set of equations is parabolic; it can be solved efficiently via a marching scheme by omitting the pressure term and using a prescribed eddy viscosity model. However, eliminating the pressure term requires generating an inflow condition for an already expanded wake some distance downstream of the rotor via wake expansion models that do not preserve the intrinsic conservation of mass and momentum for the full computational domain. Consequently, ensuring conservativeness within the Ainslie flow solver alone is insufficient to obtain a conservative wake model.

The aim of this work is to reintroduce a part of the pressure term into the Ainslie wake model that can be precomputed and is compatible with a marching scheme. An efficient closed-form solution exists analogous to the vector potential of a vortex ring, which has been integrated into a Finite Volume solver for the momentum equation. The resulting solver conserves the applied thrust force in the integral momentum deficit to very small relative errors related to the truncation of the domain size. A first comparison to LES reference simulations shows a good qualitative agreement of the obtained velocity fields, though still being limited by the lack of a suitable eddy viscosity model.

## 1. Introduction

The Dynamic Wake Meandering model [1] (DWM) is a certified method for computing farm wake effects in aeroelastic wind turbine loads computations [2], and has found wide adaptation in the industry due to its relatively moderate computational costs. These are a consequence of using a parabolised Navier-Stokes equation for the axial momentum in cylinder coordinates (Ainslie wake model [3]) in combination with an initial wake expansion model introduced to neglect the pressure term. However, the use of an expansion model based on Momentum Theory as an inflow condition bears intrinsic problems with the conservation laws for mass and momentum for the full computational domain.

Momentum theory is based on mass and momentum conservation inside a stream tube. Using the mass balance ($a = \frac{1}{2}(1 - \sqrt{1 - C_t})$, $1 - a = \frac{1}{2}(1 + \sqrt{1 - C_t})$, $1 - 2a = \sqrt{1 - C_t}$)

$$U_{\infty}\pi R_{x \to -\infty}^2 = U_{\infty}\frac{1}{2}(1 + \sqrt{1 - C_t})\pi R_{\text{rotor}}^2 = U_{\infty}\sqrt{1 - C_t}\pi R_{x \to +\infty}^2 \quad (1)$$

with the thrust coefficient $C_t = 4a(1 - a) = \frac{F_{\text{thrust}}}{\frac{1}{2}\rho U_{\infty}^2 \pi R_{\text{rotor}}^2}$ and the axial induction factor $a$, the axial momentum balance over the infinite cylinder domain $(r, \phi, x) \in \mathbb{R}^+ \times [0, 2\pi) \times \mathbb{R}$ ($r$ radial,

<figure><img src="image_3.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_4.png" /><figcaption>(a) Integral momentum deficit relative to the applied thrust force (b) Near wake at $U_\infty = 3.5 \text{ m/s}$ $F_{\text{thrust}}$ over wind speed $U_\infty$ for selected wake expansion models (DWM after expansion [4])</figcaption></figure>

<caption>Figure 1: DWM momentum conservation and near wake compared to FAWM, SG 14-236 turbine</caption>

$\phi$ azimuthal, $x$ - axial/wind direction) is given by

$$F_{\text{thrust}} = 2\pi\rho \int_{0}^{R_{x \to +\infty}} \left[ \underbrace{U_{\infty}^{2}}_{=U^{2}(r, x \to -\infty)} - \underbrace{U_{\infty}^{2}(1-C_{t})}_{=U^{2}(r, x \to +\infty)} \right] r dr + \underbrace{2\pi\rho \int_{R_{x \to +\infty}}^{\infty} U_{\infty}^{2} - U_{\infty}^{2} r dr}_{=0} \\ = 2\pi\rho U_{\infty}^{2} \int_{0}^{R_{x \to +\infty}} C_{t} r dr = C_{t}\rho U_{\infty}^{2}\pi R_{x \to +\infty}^{2} = C_{t}\frac{\rho}{2}U_{\infty}^{2}\left(1+\frac{1}{\sqrt{1-C_{t}}}\right)\pi R_{\text{rotor}}^{2}. \quad (2)$$

The full-domain momentum balance therefore differs from the momentum balance of the stream tube by a factor $1+\frac{1}{\sqrt{1-C_t}}$: for $C_t \ll 1$ it approaches the value 2; at the Betz limit with $C_t = 8/9$ the value 4, and for $C_t \to 1$ it grows unbound and is only limited by a possible truncation that avoids the wake velocities becoming zero or imaginary ($C_t \ge 1$).

Figure 1a demonstrates that this effect is also observed in wake expansion models based on Momentum Theory as used inside DWM: both models based on a single stream tube with constant $a(C_t)$ (Momentum Theory, IEC 61400-1 model) converge to the asymptotic value of 2 for large wind speeds (weak induction), reach a plateau at values ~3 below the rated wind speed, and increase steeply at very low wind speeds. When using an expansion model[4, 5] based on the outputs of a Blade Element Momentum Theory (in this case the in-house BEMT model PentaBEM), more problems can arise if the axial induction factor per stream tube is derived from the BEMT solution directly and not from the total axial force per stream tube (effectively ignoring the tip-loss correction factor), as seen for larger wind speeds. BEMT-based expansion models introduce various correction factors that violate the mass and momentum balance for the stream tubes, but typically lower the overprediction of the momentum deficit somewhat. Notwithstanding, all BEMT-based inflow models significantly overpredict the momentum deficit for the operational range of a wind turbine.

For mass- and momentum-conserving discretisations of the flow solver (Finite Volume and Finite Element methods), the overprediction of the momentum deficit and the non-zero mass deficit are carried through the full computational domain and result in an overprediction of the wake deficits further downstream. It therefore has implications for the full DWM tool chain, and requires additional empirical corrections and the use of overly strong eddy viscosity models to diffuse the over-predicted wake deficits away. For increasing farm sizes, this requires repeated
tuning and recalibration of the model parameters and raises concerns about the general validity of the DWM approach.

The approach presented in this work is an attempt to eliminate BEMT-based inflow conditions by reintroducing a part of the pressure term defined by the time-averaged forcing of the actuator disk into the momentum equation, as illustrated in Figure 1b.

## 2. Problem Statement

Mass conservation for a fluid with constant density $\rho$ is described by the continuity equation

$$ \mathrm{div} \, \vec{u} = 0 \qquad (3) $$

and the conservation of momentum of a stationary flow $\vec{u}(\vec{x})$ subject to the time-independent force density $\vec{f}$ by the Reynolds-Averaged Navier Stokes equation (RANS)

$$ \mathrm{div} (\vec{u}\vec{u}^T) = \vec{u} \cdot \mathrm{grad} \, \vec{u} = -\frac{1}{\rho} \mathrm{grad} \, p + \nu_T \Delta \vec{u} + \frac{1}{\rho} \vec{f}. \qquad (4) $$

To eliminate the dependent variables, a model must be chosen for the (scalar) eddy viscosity

$$ \nu_T(\vec{x}) = \nu_T(\vec{u}(\vec{x})), \qquad (5) $$

and a Poisson equation for the pressure $p$ is obtained by taking the divergence of Equation (4)

$$ \Delta p = \mathrm{div} \, \mathrm{grad} \, p = \underbrace{-\rho \mathrm{div} [\mathrm{div}(\vec{u}\vec{u}^T)] + \rho \mathrm{grad} \, \nu_t \cdot \Delta \vec{u} + \rho \nu_T \Delta(\underbrace{\mathrm{div} \, \vec{u}}_{=0}) + \mathrm{div} \, \vec{f}}_{=\Delta p_u} = \Delta(p_u + p_f). \qquad (6) $$

By the linearity of the Laplace operator $\Delta$, the pressure is given by a superposition of contributions from the velocity field $p_u$ and the forcing $p_f$. Similarly to Equation (4), the pressure Poisson equation is elliptic and requires a full-domain solver with prescribed boundary conditions. The computational cost for solving the system of Equations (3) to (6) are therefore prohibitively high for solvers developed for wind turbine loads computations with a tight budget of $O(10^3)$ seconds on a single CPU core only.

Equation (4) can be made parabolic by eliminating the second axial derivative in the Laplacian. This step essentially neglects diffusion in the predominant direction of advection $u$, which is often less severe than other simplifications made (e.g. axisymmetric flow). Changing Equation (4) to cylinder coordinates $(r, \phi, x)$ and assuming uniformity in the azimuthal direction with $\partial(\dots)/\partial\phi = 0$ yields the axial momentum equation of the Ainslie wake model [3]

$$ \mathrm{div} \left( \begin{pmatrix} u(r, x) \\ v(r, x) \end{pmatrix} u(r, x) \right) = u(r, x) \frac{\partial u(r, x)}{\partial x} + v(r, x) \frac{\partial u(r, x)}{\partial r} = \frac{\nu_t(x)}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u(r, x)}{\partial r} \right), \qquad (4a) $$

which can be computed efficiently by solving a 1D non-linear problem with a purely radial discretisation in $r$ for a finite step size in the downstream direction $x$ (marching scheme).

The force exerted on the air flow by the rotor of a wind turbine can be approximated by an actuator disk force $\vec{f}_{AD}$; the approach presented in the following is to find closed-form solutions for the pressure gradient of the forcing contribution $p_f$ and only neglect the elliptic problem of finding the solution of $p_u$. In extension to the parabolised Ainslie momentum equation (4a), this results in the axial momentum equation (compare [6, Art. 371b], 2D Boundary Layer equations)

$$ u(r, x) \frac{\partial u(r, x)}{\partial x} + v(r, x) \frac{\partial u(r, x)}{\partial r} = \frac{1}{\rho} f_{AD,x} + \frac{\nu_t(x)}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u(r, x)}{\partial r} \right) \qquad (4b) $$

with the forcing given either as the curl of a vector potential $\vec{f}_{AD,1} = \mathrm{curl} \, \vec{A}(r, x)$ or by the gradient of the pressure (scalar potential) $\vec{f}_{AD,2} = \mathrm{grad} \, p_f(r, x)$.
### 3. Closed-form solutions of the forcing term

#### 3.1. Vector potential

A constant axial force density $F_x$ over an infinitesimally thin circular disk with radius $R$ and centered at the origin is given by

$$ \vec{f}_{AD}(r, \phi, x) = (f_r, f_{\phi}, f_x)^T \text{ with } f_r = 0, f_{\phi} = 0 \text{ and } f_x(r, \phi, x) = F_x \delta(x)(1 - H(r - R)), \quad (7) $$

using the Heaviside function $H(r)$ and Dirac delta function $\delta(x)$. Taking the curl of this force density yields

$$ \mathrm{curl} \vec{f} = \left( \frac{1}{r} \frac{\partial f_x}{\partial \phi} - \frac{\partial f_{\phi}}{\partial x}, \frac{\partial f_r}{\partial x} - \frac{\partial f_x}{\partial r}, \frac{1}{r} \frac{\partial(r f_{\phi})}{\partial r} - \frac{1}{r} \frac{\partial f_r}{\partial \phi} \right)^T = (0, F_x \delta(r-R) \delta(x), 0)^T . \quad (8) $$

The azimuthal component $F_x \delta(r - R) \delta(x)$ of curl $\vec{f}$ is mathematically identical to the vorticity of an infinitesimal vortex ring [7, §6][6, Art. 161] or to the current density of a circular loop [8, 5.5], and its solution is given by the (azimuthal component of the) vector potential

$$ A_{\phi}(r, x) = \frac{F_x}{2\pi r} \frac{1}{k} \left[ (2 - k^2)K(k^2) - 2E(k^2) \right] . \quad (9) $$

with the complete elliptic integrals of the first $K(k^2)$ and second kind $E(k^2)$ and the square elliptic modulus $k^2 = \frac{4rR}{(r+R)^2+x^2}$. $K(k^2)$ has a logarithmic singularity for $k^2 \to 1$ ($r \to R$, $x \to 0$); a desingularised square modulus $k_{\epsilon}^2 = \frac{4rR}{(r+R)^2+x^2+\epsilon^2}$ with a small parameter $\epsilon^2$ can be used to avoid this singularity numerically. Equation (9) suffers from floating point cancellation errors for small arguments; for $k^2 < 0.01$ the (factorised) Taylor series expansion

$$ A_{\phi,k^2(r,x)<0.01}(r,x) \approx \frac{F_x}{r} \left( \frac{k^3}{32} + \frac{3k^5}{128} + \frac{75k^7}{4096} + \frac{245k^9}{16384} + \frac{6615k^{11}}{524288} + \frac{22869k^{13}}{2097152} + \mathcal{O}(k^{15}) \right) \quad (10) $$

has lower numerical errors and is also computationally more efficient.

Asymptotically [6, 9], Equation (9) approaches the vector potential of a dipole

$$ A_{\phi, r^2+x^2 \gg R^2}(r, x) \approx \frac{F_x \pi R^2 \sin \phi}{4\pi(r^2+x^2)}, \quad (11) $$

which allows to combine the far-field forcing of multiple co-located vortex rings by the area-weighted sum $\sum_i F_{x,i} \pi R_i^2$. In analogy to vortex rings [9], the consequence of the weak asymptotic convergence is that the integral $\int f_x dV$ is divergent and does not exist for the infinite domain (see also [10]), which requires computing the total forcing for the truncated finite computational domain as detailed in Section 4.3.

#### 3.2. Forcing term via the scalar potential

Taking instead the divergence of the forcing as found in the pressure Poisson equation results in

$$ \Delta p_f = \mathrm{div} \vec{f} = F_x \frac{\partial}{\partial x} \delta(x)(1 - H(r - R)), \quad (12) $$

and a convolution with the fundamental solution of the Laplacian $G(\vec{x}, \vec{x}') = \frac{1}{4\pi|\vec{x} - \vec{x}'|}$ gives

$$ p_f = \underbrace{\left[ \frac{f_x}{4\pi|\vec{x} - \vec{x}'|} \right]}_{=0} \to \frac{1}{4\pi} \int_{\mathbb{R}^3} \frac{f_x \hat{e}_x \cdot (\vec{x} - \vec{x}')}{|\vec{x} - \vec{x}'|^3} dV = \frac{F_x}{4\pi} \int_{\mathrm{AD}} \frac{(\vec{x} - \vec{x}') \cdot \vec{n}}{|\vec{x} - \vec{x}'|^3} dS. \quad (13) $$
The convolution integral obtained after one step of integration-by-parts is well-known in Boundary Element Methods [11] as the induced potential of a uniform dipole or doublet distribution over the disk. Geometrically, it has the interpretation of the solid angle covered by the disk as seen from the target point of induction. Solutions for the solid angle of a circular disk exist [12, 13], albeit more complex ones than the vector potential that require quadrature for the point evaluation instead. Solutions via the scalar potential were therefore not implemented.

### 3.3. The Forced Ainslie Wake Model (FAWM)

The solver is based on a second-order accurate (first-order on non-uniform spacing) Finite Volume Method on a staggered mesh in the $r-x$-plane for the Ainslie transport model (Parabolised Navier Stokes Equation, PNSE). Eliminating the radial velocity component by a radial integral of the Continuity Equation results in a non-linear system of equations solved by Newton-Raphson iteration. The Jacobian has a lower Hessenberg form, resulting in a $O(N_{radial}^2)$ complexity for the matrix inversion and overall a very fast solver (single core wall clock times of approximately 10 s on an AMD Ryzen 5 PRO 3500U using meshes with $2^{20}$ nodes). The solver is complemented by an adaptive sub-stepping (reduction and relaxation of the axial step size) that improves convergence in zones of strong velocity gradients.

The solver implements the IEC 61400-1 prescribed eddy viscosity model [2] as well as a constant prescribed eddy viscosity over the full computational domain. The IEC 61400-1 model was calibrated in [14] for an initial wake expansion model based on Momentum Theory with correction factors [4]. When used with the forcing term, the radial gradients in the diffusion term are much weaker behind the rotor disk, and diffusion is overall strongly underestimated, resulting in locally zero or negative velocities and an invalid state for PNSE solvers. The computations typically fail a short distance behind the actuator disk at almost all ambient wind speeds and turbulence intensities besides for very weakly loaded rotors at wind speeds significantly above the rated wind speed. Therefore, all solver results in this work were obtained with a constant prescribed eddy viscosity for the full domain.

### 3.4. Far-Field Approximation and Projection on Finite Volume Cell

The maximum relative error in the forcing term by combining Equations (9) and (10) is $< 10^{-11}$.
The projection of the forcing term onto a Finite Volume cell is given by

$$ \int_{x_i}^{x_{i+1}} \int_{r_j}^{r_{j+1}} \frac{1}{r} \frac{\partial}{\partial r} (rA_{\phi}(r, x)) r dr dx = r \int_{x_i}^{x_{i+1}} A_{\phi}(r_{j+1}, x) - A_{\phi}(r_j, x) dx. \quad (14) $$

The remaining integral over the axial interval is computed using a Gauss-Legendre quadrature that converges even in close proximity to the logarithmic singularity. An adaptive number of quadrature points was implemented to reduce the order in the far field, given by

$$ N_{\text{quadrature}}(r, x) = \min \left( 63, \left\lfloor 2^{2-\min\left(1, \frac{1}{2}\log_{10}\left(\frac{(r-R)^2+x^2}{R^2}\right)\right)} \right\rfloor \right). \quad (15) $$

Its parameters were empirically found and tested to ensure a maximum relative error of $10^{-12}$.

## 4. Calibration with reference CFD simulations

In order to obtain reference numerical results, CFD simulations have been performed using the incompressible finite-volume solver of the YALES2 library [15, 16]. This library was specifically designed for high-performance LES simulations of multiphysics problems on unstructured meshes. The numerical schemes used in the mentioned solver are 4th-order on regular grids in both space and time. The sub-grid scales tensor is calculated using the Dynamic Smagorinsky model [17].
<table><caption>Table 1: Configurations simulated with the YALES2 library.</caption><thead><tr><th>Config. #</th><th>U<sub>∞</sub> [m/s]</th><th>I<sub>REF</sub></th></tr></thead><tbody><tr><td>1</td><td>4.0</td><td>0.16</td></tr><tr><td>2</td><td>11.4</td><td>0.12</td></tr><tr><td>3</td><td>11.4</td><td>0.16</td></tr><tr><td>4</td><td>19.25</td><td>0.16</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 2: Domain size and cell size $\delta_c$ mapping used in all configurations from Table 1.</figcaption></figure>

### 4.1. Reference simulation setup

The simulated configurations are given in Table 1. The quantity $I_{REF}$ refers to the notation used in the IEC standard 61400-1 to define the turbulence classes. The turbulence intensity (TI) at the domain inlet is thus defined as follows:

$$ \text{TI} = \frac{\sigma_1}{U_{\infty}} \qquad (16) $$

where $\sigma_1$ corresponds to the 90% quantile of the velocity fluctuations amplitude, defined in the IEC standard as:

$$ \sigma_1 = I_{REF} (0.75 U_{\infty} + 5.6) \qquad (17) $$

In each case, the rotor of the DTU-10MW reference wind turbine [18] has been modeled using the Actuator Line Method (ALM) [19] implemented as a module in the YALES2 library [20, 21]. The computational domain and the complete cell size mapping are shown in Figure 2. The total cell count was over 230 millions. The usual guidelines coming with the ALM have been followed carefully; to facilitate, this work does not make use of the classical corrections on the polars or on the aerodynamic forces given by the Blade Element method. The rotor speed and collective pitch angle are set in each case through linear interpolation of the design data corresponding to an infinitely rigid blade. Consistently with the actuator disk formulation used in FAWM, the cone and tilt angles have been set to zero in all simulations, ignoring the actual design values. Second, to obtain axisymmetric inflow conditions, the inlet flow field is obtained by superposing the steady uniform wind speed $U_{\infty}$ and synthetic turbulence [22, 23], set to be homogeneous and isotropic (without vortex stretching ($\Gamma = 0$) and high-frequency correction).

A two-step approach was used to collect the reference time-averaged loads along the blades. First, the inlet turbulence is convected downstream for a period of time $T_1 = T_{ft}$, where $T_{ft}$ is the flow-through time of the injected turbulence fields. This step also allows the wind turbine wake to build up and converge. The second step consists of accumulating, during a period $T_2 = 2 T_{ft}$, the statistics of the blade loads given in the blade coordinate system, at the radial position of each actuator point.
<figure><img src="image_6.png" /><figcaption>Figure 3: Piece-wise constant axial force Figure 4: Systematic and residual solver error in densities ($\rho\pi(r_{i+1}^2 - r_i^2)^{-1} \sum_{j=1,2,3} F_{i,j}$ the momentum deficit integral, comp. Figure 1a</figcaption></figure>

## 4.2. Deducing the Actuator Disk Forcing from the Actuator Line Model

Given the total time-averaged axial force $F_{i,j}$ per span-wise segment $[r_i, r_{i+1}]$ and blade $j$, the (kinematic) strength of the forcing ring located at $r_i$ is obtained by

$$F_{x,i} = \frac{1}{\rho} \left( \frac{\sum_{j=1,2,3} F_{i,j}}{\pi(r_{i+1}^2 - r_i^2)} - \frac{\sum_{j=1,2,3} F_{i-1,j}}{\pi(r_i^2 - r_{i-1}^2)} \right). \quad (18)$$

Figure 3 shows the piece-wise constant force densities ($\rho\pi(r_{i+1}^2 - r_i^2)^{-1} \sum_{j=1,2,3} F_{i,j}$ of the actuator disk as obtained from the YALES2 Actuator Line model.

## 4.3. Conservation of Integral Properties

The chosen numerical scheme fulfils the continuity equation point-wise, and the resulting relative errors in the integral conservation of the mass flux are of the order of $10^{-14}$. Without the forcing term, the solver conserves momentum to relative errors of about $10^{-12}$. With the forcing term, the solver has an unavoidable error

$$\epsilon_{\text{domain}}(x) = \rho_{\text{air}} \left( \int_{0}^{R_{\text{max}}} U(x, r)^2 - U_{\infty}^2 r dr + \int_{X_{\text{min}}}^{x} U(\tilde{x}, R_{\text{max}})V(\tilde{x}, R_{\text{max}})d\tilde{x} \right) - \sum_{i} \sum_{j=1,2,3} F_{i,j} \quad (19)$$

due to the divergent force integral of a forcing ring and the finite size of the computational domain $\Omega = [X_{\text{min}}, X_{\text{max}}] \times [0, R_{\text{max}}]$ that truncates the total forcing. Figure 4 shows $\epsilon_{\text{domain}}(x)$ for $U_{\infty} = 4\frac{m}{s}$ and $\nu_t = 10\frac{m^2}{s}$. This systematic error can be as large as 11% and strongly depends on the aspect ratio of the computational domain; choosing $R_{\text{max}} \gg |X_{\text{max}} - X_{\text{min}}|$ can reduce this error by two orders of magnitude or more. Taking the total forcing over the computational domain into account, the residual solver error in the momentum conservation can be defined as

$$\epsilon_{\text{solver}}(x) = \epsilon_{\text{domain}}(x) - \int_{X_{\text{min}}}^{X_{\text{max}}} \int_{0}^{R_{\text{max}}} f_x(x, r) r dr dx \quad (20)$$

This resulting error in the momentum conservation is about two to three orders of magnitude smaller than $\epsilon_{\text{domain}}$; Figure 4 again shows lower errors for larger radial domain sizes. The error magnitudes obtained here should be compared to the systematic overprediction of the integral momentum deficit in BEMT-based wake expansion models shown in Figure 1a.

The radial mesh has consequently been adapted to a uniform spacing near the axis of symmetry that gradually changes to an exponentially increasing step size in the far field. This allows to maintain a spatial resolution of $\le 1m$ in the zone of interest $0 \le r \le D$ while extending the computational domain to $R_{\text{max}} > 100D$ with about 500 nodes in the radial direction.
<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution of the cumulative errors (MAPE and RRMSE) between FAWM and YALES2 simulations for several configurations from Table 1, along with the corresponding calculation masks $\mathcal{M}$.</figcaption></figure>

### 4.4. Comparison YALES2 - FAWM

For comparing the streamwise velocity fields $u_{Y2}$ and $u_F$ obtained with YALES2 and FAWM in the x-r-plane, two metrics will be used in the following, the Mean Absolute Percentage Error (MAPE) and the Relative Root Mean Square Error (RRMSE):

$$ \text{MAPE}(u_{Y2}, u_F) = \frac{1}{n} \sum_{i}^{n} \left| \frac{u_{Y2,i} - u_{F,i}}{u_{Y2,i}} \right|, \quad \text{RRMSE}(u_{Y2}, u_F) = \frac{\sqrt{\frac{1}{n} \sum_{i}^{n} (u_{Y2,i} - u_{F,i})^2}}{\frac{1}{n} \sum_{i}^{n} |u_{Y2,i}|} \quad (21) $$

Such metrics were calculated on a subpart of the full domain, given by the mask $\mathcal{M}(x, r)$:

$$ \mathcal{M} = (|r| \le D) \& (x \ge -2D) \& (x \le 3D) \& \left( \left| \frac{u_F - U_{\infty}}{U_{\infty}} \right| \ge 0.01 \right) \quad (22) $$

For each configuration in Table 1, several FAWM flow fields were simulated using different values of the global eddy viscosity $\nu_t$ in the range $U_{\infty} \cdot [0.5, 2.0]$. The resulting MAPE and RRMSE are shown in Figure 5 along with illustrations of the corresponding calculation masks $\mathcal{M}(x, r)$. The resulting errors are reasonably low overall, but mostly above 5% in the range below the rated wind speed. The value of $\nu_t$ that minimizes error metrics does not vary linearly with $U_{\infty}$, but shows a significant evolution (see configurations 1 and 3). It should be emphasized again that the results are obtained from a subpart of the domain, including only the upstream induction zone and the near-wake region.

Following this initial calibration of the eddy viscosity, the time-averaged velocity fields can be compared with a focus on the near-wake expansion and velocity deficit profiles, as illustrated in Figure 6. In general, the velocity fields and wake expansion rates are very similar in FAWM and YALES2 results, which is consistent with the approach used to calibrate $\nu_t$ in FAWM. The velocity deficit profiles match rather well up to a downstream distance of 3D for all configurations, and even further downstream for the above-rated wind speed, with a tendency of FAWM to overpredict both the wake velocity deficit and diffusion. For configurations 1 and 3 with higher thrust coefficients, a significant mismatch in velocity deficits can be observed in the far wake region, as FAWM underestimates both wake expansion and recovery.

Finally, a closer look is needed at the impact of inlet turbulence intensities on the YALES2 results (see Figure 7). Relative to the target magnitude $u_{RMS,T}$ of the RMS velocities in configuration 3, the inlet turbulence significantly decayed before reaching the wind turbine. Although still noticeable, the initial difference in turbulence intensities between configurations 2 and 3 is reduced at the rotor position. Ultimately, the resulting time-averaged velocity fields
<figure><img src="image_8.png" /><figcaption>Figure 6: Side by side comparison of the non-dimensional streamwise time-averaged velocity fields for selected configurations of Table 1, including a comparison of the streamwise velocity deficits at several downstream positions. All FAWM setups were using the optimal $v_t$ values among those assessed to minimize the MAPE/RRMSE with respect to the YALES2 results.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Overview of the RMS (left), time-averaged and instantaneous (right) velocity fields obtained with YALES2 in configurations 2 and/or 3 from Table 1.</figcaption></figure>

are almost identical in both cases. This setback will need to be addressed if extended FAWM calibrations with YALES2 results are to be performed.

Using an adaptive Gauss-Legendre quadrature of the forcing term in the far field, single core wall clock times of about 200 s have been achieved on an AMD Ryzen 5 PRO 3500U using 100$D \times 200$D meshes with 2$^{20}$ nodes and 65 forcing rings, with about 85% of the computational time spent on the elliptic integrals in the forcing term. Order reductions of the forcing rings have been implemented and reduce the simulation time further, but the current implementation does not preserve the higher-order terms in the multipole expansion and thereby introduces small additional errors. All presented results were obtained without the order reduction. The multipole expansion and the vectorisation of the forcing term computation remain ongoing work.
**5. Conclusion**

This work presents an attempt to introduce the forcing term of a circular actuator disk into the axial momentum equation. Using a closed-form solution for the vector potential avoids solving an elliptic pressure Poisson equation while retaining compatibility with the parabolised axial transport equation introduced in the Ainslie model. Some means (adaptive quadrature, far-field reduction of the number of vortex rings) were introduced and tested to reduce the computational cost of calculating this forcing term.

The resulting flow fields compare surprisingly well to reference simulations when an optimal global eddy viscosity parameter is chosen; wake expansion and sizes are very similar, but some discrepancies are observed in the near- and far-field due to diffusion being locally over- or under-predicted. The match can be expected to improve once calibrated eddy viscosity models are available; the IEC 61400-1 model was found to be incompatible with this solver. To conserve the total axial force in the momentum deficits, a very large radial domain size must be chosen, otherwise the spatial truncation of the forcing term introduces a systematic error. Furthermore, a weak residual solver error remains in the integral momentum conservation whose origin needs to be clarified, but which is likely related to the omitted terms in the pressure equation. Nonetheless, FAWM achieves several orders of magnitude lower integral conservation errors than the BEMT-based wake expansion models in DWM, and sufficiently short simulation times make the solver usable for wind turbine loads computations.

**Acknowledgement:** The LES results were obtained using the Austral cluster of CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie), www.criann.fr, under the allocation 2012006.

[1] Larsen G, Madsen H A and Bingoel F 2007 *Dynamic wake meandering modeling*

[2] International Electrotechnical Commission et al. 2019 IEC 61400-1 Ed.4 Wind Turbine Generator Systems Part 1, International Standard, International Electrotechnical Commission

[3] Ainslie J F 1988 *Journal of Wind Engineering and Industrial Aerodynamics* **27** 213–224

[4] Madsen H A, Larsen G, Larsen T, Troldborg N and Mikkelsen R 2010 *Journal of Solar Energy Engineering* **132** 041014 (14 pages) ISSN 0199-6231

[5] Keck R E, Veldkamp D, Wedel-Heinen J and Forsberg J 2013 *A consistent turbulence formulation for the dynamic wake meandering model in the atmospheric boundary layer* Ph.D. thesis Denmark

[6] Lamb H 1993 *Hydrodynamics* (Cambridge University Press)

[7] Helmholtz H 1858 *3. Über Integrale der hydrodynamischen Gleichungen, welche den Wirbelbewegungen entsprechen* (Berlin, Boston: De Gruyter) pp 25–55 ISBN 9783112336489

[8] Jackson J D 2021 *Classical electrodynamics* (John Wiley & Sons)

[9] Akhmetov D G 2009 *Vortex rings* (Springer Science & Business Media)

[10] Paddeu J 2024 *Finite Element Method for Wind Turbine Wake Modelling* Master's thesis Delft University of Technology URL http:// resolver.tudelft.nl / uuid:5234e49e-42b0-45a1-bd0c-d5f100a72cb3

[11] Sauter S A, Schwab C, Sauter S A and Schwab C 2011 *Boundary Element Methods* (Springer)

[12] Masket A V 1957 *Review of Scientific Instruments* **28** 191–197

[13] Paxton F 1959 *Review of Scientific Instruments* **30** 254–258

[14] Larsen T, Madsen H A, Larsen G and Hansen K 2013 *Wind Energy* **16** 605–624 ISSN 1095-4244

[15] Moureau V, Domingo P and Vervisch L 2011 *Comptes Rendus Mécanique* **339** 141–148

[16] Moureau V and Lartigue G 2023 *Yales2* public page https://www.coria-cfd.fr/index.php/YALES2

[17] Lilly D K 1992 *Physics of Fluids A: Fluid Dynamics* **4** 633–635

[18] Bak C, Zahle F, Bitsche R, Kim T, Yde A, Henriksen L, Hansen M, Blasques J, Gaunaa M and Natarajan A 2013 *The DTU 10-MW Reference Wind Turbine* https://orbit.dtu.dk/en/publications/the-dtu-10-mw-reference-wind-turbine

[19] Sørensen J N and Shen W Z 2002 *Journal of Fluids Engineering* **124** 393–399

[20] Benard P, Viré A, Moureau V, Lartigue G, Beaudet L, Deglaire P and Bricteux L 2018 *Computers & Fluids* **173** 133–139 ISSN 0045-7930

[21] Houtin-Mongrolle F 2022 *Investigations of yawed offshore wind turbine interactions through aero-servo-elastic Large Eddy Simulations* Ph.D. thesis Normandie Université

[22] Mann J 1994 *Journal of fluid mechanics* **273** 141–168

[23] Mann J 1998 *Probabilistic engineering mechanics* **13** 269–282