*2.2. Flow Solver*

The turbulent flow is solved with a three dimensional LES code, dubbed as virtual flow simulator (VFS-Wind) [28], in which the governing equations are the filtered Navier–Stokes equations for incompressible flows, which read in the compact tensor form as (*i*, *j*, *k*, *l* = 1, 2, 3):

$$J\frac{\partial L^i}{\partial \xi^i} = 0,\tag{5}$$

$$\frac{1}{J}\frac{\partial \mathcal{U}^{j}}{\partial t} = \frac{\xi^{i}\_{l}}{J} \left( -\frac{\partial}{\partial \xi^{j}} \left( \mathcal{U}^{j} u\_{l} \right) + \frac{\mu}{\rho} \frac{\partial}{\partial \xi^{j}} \left( \frac{\mathcal{g}^{jk}}{J} \frac{\partial u\_{l}}{\partial \xi^{k}} \right) - \frac{1}{\rho} \frac{\partial}{\partial \xi^{j}} \left( \frac{\mathcal{z}^{j}\_{l} p}{J} \right) - \frac{1}{\rho} \frac{\partial \pi\_{l j}}{\partial \xi^{j}} + f\_{l} \right), \tag{6}$$

where *ξ<sup>i</sup>* is the curvilinear coordinates, *ξ<sup>i</sup> <sup>l</sup>* = *∂ξi*/*∂xl* is the transformation metrics with *xl* the Cartesian coordinates, *J* denotes the Jacobian of the geometric transformation, *U<sup>i</sup>* = *ξi l* /*J ul* is the contravariant volume flux with *ul* the velocity in Cartesian coordinates, *μ* is the dynamic viscosity, *ρ* is the density, *gjk* = *ξ j l ξk <sup>l</sup>* is the components of the contravariant metric tensor, *p* is the pressure, *fl* is the body forces exerted by the actuator models, and *τij* is the sub-grid stress (SGS) resulted from the filtering operation and is modeled with the Smagorinsky SGS model [29] as follows,

$$
\pi\_{i\bar{j}} - \frac{1}{3}\pi\_{kk}\delta\_{i\bar{j}} = -\mu\_t \overline{S\_{i\bar{j}\bar{\prime}}}\tag{7}
$$

where *μ<sup>t</sup>* is the eddy viscosity and *Sij* is the large-scale strain-rate tensor with (·) denoting the grid filtering operator. The eddy viscosity is computed by

$$
\mu\_t = \mathbb{C}\_s \Delta^2 |\overleftarrow{\mathbb{S}}| \,\tag{8}
$$

where <sup>Δ</sup> is the filter width, |*S*| = (2*SijSij*)1/2 is the magnitude of the strain-rate tensor with *Cs* the Smagorinsky constant computed via the dynamic procedure of [30].

A second-order accurate central differencing scheme is used for space discretization. The time integration uses the fractional step method [31]. The momentum equation is solved with a matrix-free Newton–Krylov method [32] . The pressure Poisson equation is solved with a Generalized Minimal Residual (GMRES) method with an algebraic multi-grid acceleration [33].
