Developing a digital twin framework for wind tunnel testing: validation of turbulent inflow and airfoil load applications

Wind energy systems, such as horizontal-axis wind turbines and vertical-axis wind turbines, operate within the turbulent atmospheric boundary layer, where turbulence significantly impacts their efficiency. Therefore, it is crucial to investigate the impact of turbulent inflow on the aerodynamic performance at the rotor blade scale. As field investigations are challenging, in this work, we present a framework where we combine wind tunnel measurements in turbulent flow with a digital twin of the experimental set-up. For this, first, the decay of the turbulent inflow needs to be described and simulated correctly. Here, we use Reynolds-averaged Navier–Stokes (RANS) simulations with k−ω turbulence models, where a suitable turbulence length scale is required as an inlet boundary condition. While the integral length scale is often chosen without a theoretical basis, this study derives that the Taylor micro-scale is the correct choice for simulating turbulence generated by a regular grid: the temporal decay of turbulent kinetic energy (TKE) is shown to depend on the initial value of the Taylor micro-scale by solving the differential equations given by Speziale and Bernard (1992). Further, the spatial decay of TKE and its dependence on the Taylor micro-scale at the inlet boundary are derived. With this theoretical understanding, RANS simulations with k−ω turbulence models are conducted using the Taylor micro-scale and the TKE obtained from grid experiments as the inlet boundary condition. Second, the results are validated with excellent agreement with the TKE evolution downstream of a grid obtained through hot-wire measurements in the wind tunnel. Third, the study further introduces an airfoil in both the experimental and the numerical setting where 3D simulations are performed. A very good match between force coefficients obtained from experiments and the digital twin is found. In conclusion, this study demonstrates that the Taylor micro-scale is the appropriate turbulence length scale to be used as the boundary condition and initial condition to simulate the evolution of TKE for regular-grid-generated turbulent flows. Additionally, the digital twin of the wind tunnel can accurately replicate the force coefficients obtained in the physical wind tunnel.

1 Introduction

Wind energy systems, such as horizontal-axis and vertical-axis wind turbines, operate in a turbulent atmospheric boundary layer, which significantly affects their efficiency. Therefore, it is essential to study the turbulent inflow that they encounter. Important statistical quantities that considerably affect the aerodynamic performance of a rotor blade are the turbulent kinetic energy (TKE) and length scales in the wind. To study their effects, field experiments can be carried out, but they are complex, time-consuming, and costly. Alternatively, experiments can be conducted in a wind tunnel by subjecting a Reynolds-scaled blade section from a real wind turbine blade to turbulent inflow under different inflow conditions, such as homogeneous inflow or gust inflow (see e.g. (Wei et al., 2019; Wester et al., 2022; Mishra et al., 2022; Nietiedt et al., 2022) and (Neunaber and Braud, 2020a)). However, wind turbines have grown to become the largest flexible, rotating machines in the world, with blade lengths approaching now 120 m. The interaction between a highly variable inflow and the unsteady aerodynamics of the moving and deforming blades is pushing the limits of current theory (Veers et al., 2019). At the blade scale, chord-based Reynolds numbers will exceed 15 million, which is unreachable in most available wind tunnels, and reach limits of pressurized facilities that are specifically developed for that purpose (Miller et al., 2019; Brunner et al., 2021). Therefore, having a digital twin model is key to help in the development of very large horizontal-axis wind turbines (HAWTs). Indeed, with these digital twin models, aerodynamic loads from local blade sections, used as input of blade element momentum (BEM) solvers, can be pushed towards very large Reynolds numbers in numerous configurations encountered locally by blade sections (due to elasticity and floating movements). This would also provide the full flow and load description to understand further unsteady aerodynamics of blades in such configurations and thus further improve blade design of very large HAWTs. Among the computational fluid dynamics (CFD) methods available for that purpose, a Reynolds-averaged Navier–Stokes (RANS) formulation will be preferred over solving the Navier–Stokes equations directly with numerical methods using direct numerical simulation (DNS) or large eddy simulation (LES), as it is computationally too expensive for the targeted system. However, many challenges still remain and some of them will be tackled in this paper at low Reynolds number as a first step to be able to provide experimental validation. It concerns the following:

1. the theoretical description of the decay properties of turbulence measures such as the TKE and the dissipation rate in the downstream flow direction

2. the numerical replication of the decay properties of turbulence measures such as the TKE and the dissipation rate in the downstream flow direction

3. the reproduction of the aerodynamic loads acting on an airfoil by means of RANS simulations.

The literature review is therefore split into three sub-sections to give the reader an overview over decaying turbulence, focusing on the known turbulent flow properties behind grids, the great effort made in the literature to replicate these properties in simulated flows, and the impact of turbulence on the aerodynamics of an airfoil.

Brief theoretical description of decaying grid turbulence

The Navier–Stokes equations are nonlinear, and their solutions are non-unique in nature. Therefore, simplifications are often made when describing or simulating turbulent flows, and here, we focus on literature on grid-generated turbulence (GGT). Over the decades, GGT has been investigated intensely, and different works looked into the empirical and analytical description of the evolution of the turbulence decay. According to Batchelor and Townsend (1948), the decay of the turbulence intensity can be described by means of a power law. It has further been discussed that the inlet conditions play a role in the decay of the turbulence (Kurian and Fransson, 2009). The works of Comte-Bellot and Corrsin (1966) show experimentally that the ratio between TKE ($k=\frac{\mathrm{1}}{\mathrm{2}}(\stackrel{\mathrm{‾}}{u_{\mathrm{1}}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{u_{\mathrm{2}}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{u_{\mathrm{3}}^{\mathrm{2}}})$ with the velocity fluctuations u₁, u₂, and u₃ in three axis directions) and dissipation ($\mathit{ϵ}=\mathrm{30}\mathit{\nu }\frac{\stackrel{\mathrm{‾}}{u_{\mathrm{1}}^{\mathrm{2}}}}{{\mathit{\lambda }}^{\mathrm{2}}}$, with the kinematic viscosity ν and the Taylor micro-scale λ) evolves linearly in the downstream direction. Krogstad and Davidson (2009) showed that grid turbulence is Saffman turbulence (Saffman, 1967). They improved the decay exponent of TKE from 1.2, which Saffman gave for perfectly homogeneous, isotropic turbulence, to 1.1 for GGT. Sinhuber et al. (2015) performed experiments with one grid for different grid-mesh-size-based Reynolds numbers (${\mathit{Re}}_{\mathrm{M}}=\frac{UM}{\mathit{\nu }}$, where M is the grid mesh size, and U the mean velocity) and found that the decay exponent of the TKE was equal to 1.18. They also showed that the decay exponent was independent of Re_M. A literature review can be found in Kurian and Fransson (2009). We will use the existing framework that was very briefly summarized above as a starting point for our theoretical framework, and a detailed description of both is given in Sect. 2.

Simulating decaying grid turbulence

Simulating the inflow environment using computational fluid dynamics (CFD) has been the subject of many research works, and different turbulent formulations can be used that are summarized below from the most expensive computational effort to the least one, i.e. RANS models.

Nagata et al. (2008) performed DNS to simulate the turbulent mixing layer for grid-generated turbulence, and they replicated typical GGT including the shear mixing layer. Continuing this work, Suzuki et al. (2010) showed that for a given mesh Reynolds number (Re_M), turbulent mixing is enhanced for fractal grids compared to regular grids. Laizet and Vassilicos (2011) simulated both a regular and different fractal grids using DNS, and they confirmed the characteristic regions of turbulence production and decay downstream of fractal grids. Continuing their work, Laizet et al. (2013) studied inter-scale energy transfer of decaying turbulence for fractal grids. Efforts have also been made to use LES to simulate GGT. For example, Blackmore et al. (2013) developed a grid inlet technique to generate high-intensity turbulence for given length scales in LESs. With this, they successfully imitated the independence of the decay rate of the turbulence intensity (TI) from the mesh-size-based Reynolds number Re_M. Rieth et al. (2014) compared two LES models, namely the sigma model (Nicoud et al., 2011) and the Smagorinsky model (Smagorinsky, 1963), to simulate grid-generated turbulence. They found that the sigma model is a good alternative to the static Smagorinsky model and comparable to the dynamic Smagorinsky model (Germano et al., 1991). Further, Liu et al. (2017) employed 3D LES to simulate GGT, and they found the same trend of absolute values of the mean TI. Djenidi (2006) used the lattice Boltzmann method to simulate GGT, and they simulated the decay power law over a short distance from the grid but found it difficult to find the decay exponent.

Finally, Torrano et al. (2015) have investigated the performance of various two-equation eddy viscosity models for predicting the decay of the turbulent kinetic energy downstream of a regular grid by means of RANS. They calculated the dissipation rate using an integral length scale equal to the grid mesh size. The simulation results were compared with experimental results, and a conclusion of their work was that eddy viscosity models were over-predicting the turbulent kinetic energy in comparison to the experiments. As an alternative to two-equation eddy viscosity models, Reynolds stress transport models (RSTMs) can also be used where transport equations for each of the Reynolds stress terms are solved. For example, Panda et al. (2018) performed RANS simulations using RSTMs. However, because of the closure problem, using RSTMs is quite complex as the nine non-closed components of the nine Reynolds stress transport equations have to be modelled.

Regarding RANS equations, studies focus mostly on improving the eddy viscous model, while not much effort was put on the choice of the length scale that is usually assumed to have an order of magnitude of the integral length scale. Two quantities are generally used at the inlet of RANS simulations: the turbulence intensity and the length scale. Indeed, when using the most popular two-equation eddy viscosity models, namely the k−ϵ and k−ω models (Wilcox, 1988), or a mixture of these two models, the k−ω shear stress transport (SST) Menter 1994 and 2003 models (Menter, 1994; Menter et al., 2003), the following quantities need to be provided at the inlet boundary: k, the turbulent kinetic energy, ϵ, the dissipation rate, and ω, the specific dissipation rate. k is generally given through a turbulence intensity quantity, and ϵ and ω are computed using k and a given characteristic turbulence length scale chosen individually for each flow configuration. For example, in the case of pipe flow, the turbulence length scale is empirically approximated as 0.07 times the diameter of the pipe (Versteeg and Malalasekera, 1995).

We would like to emphasize that a correctly evaluated and properly defined turbulence length scale is a basic requirement to accurately capture the evolution of turbulence. Finally, in principle, if provided with an adequate closure, the RANS model should be able to capture statistical measures of turbulence and their evolution in the flow field. The present work will focus on improving the length scale choice to reproduce the correct evolution of k.

Experimental and numerical aerodynamic research

As stated above, the third challenge of the work presented here is the reproduction of aerodynamic loads on a 2D airfoil section by means of RANS simulations in particular for low Reynolds numbers. There are many experimental works that describe the evolution of the lift and drag coefficients with the variation of the angle of attack (AoA) of the airfoil, and a summary of the design and aerodynamics for wind turbine blades can for example be found in Bak (2022). For this work, the impact of low Reynolds numbers and the impact of turbulence on the performance of 2D blade sections are of interest. Different studies on the influence of the turbulence intensity in the inflow on the aerodynamics of a 2D blade section are summarized in Li and Hearst (2021). There, an NREL S826 profile is exposed to different turbulence levels up to 5.4 % for a chord-based Reynolds number of ${\mathit{Re}}_{\mathrm{c}}=\mathrm{4}\cdot {\mathrm{10}}^{\mathrm{5}}$ defined as ${\mathit{Re}}_{\mathrm{c}}=\frac{c\cdot U_{\mathrm{\infty }}}{\mathit{\nu }}$, where c denotes the chord length, U_∞ denotes the inflow velocity, and ν the kinematic viscosity. They found that the slope of linear part of the lift curve increases with increasing turbulence intensity and that turbulence intensities of up to 1.6 % led to a reduction of the maximum lift, whereas turbulence intensities of 2.1 % and higher led to an increase in lift compared to the reference case. Devinant et al. (2002) and Sicot et al. (2006) studied the influence of turbulence intensities up to 16 % on the aerodynamics of a 2D NACA65(4)421 blade section. The chord-based Reynolds numbers were $\mathrm{1}\cdot {\mathrm{10}}^{\mathrm{5}}\le {\mathit{Re}}_{\mathrm{c}}\le \mathrm{7}\cdot {\mathrm{10}}^{\mathrm{5}}$. Devinant et al. (2002) showed that increasing the turbulence intensity shifts the stall angle towards higher angles of attack. This is attributed to a turbulent boundary layer flow that is known to be less prompt to flow separation, which also displaces the transition towards the trailing edge. Sicot et al. (2006) found that the fluctuations of the surface pressure measurements, characterized by the standard deviation, increased in the separated flow region with increasing turbulence intensity. The average location of the separation line was not affected.

Simulating turbulent flow upstream of the airfoil often requires the use of LES, as demonstrated in the work of Gilling et al. (2009). In this paper, we will show that with appropriate boundary conditions at the inlet, RANS simulations can yield an accurate evolution of turbulence properties.

Finally, fewer investigations of the impact of turbulence have been performed at higher Reynolds numbers due to experimental complexity. The present digital twin model will pave the way for such studies.

Structure of this paper

As stated above, the aim of this work is the creation of a digital twin of a low-Reynolds-number wind tunnel where turbulence is generated with grids. For this, first, a theoretical framework is developed in Sect. 2 to demonstrate that the Taylor micro-scale (λ) is the correct length scale to be used as the turbulence length scale in RANS simulations for homogeneous, isotropic turbulent flows. Additionally, we show that the spatial and temporal decays of turbulent kinetic energy are directly dependent on the Taylor micro-scale, and we derive a relation between the Taylor micro-scale and the downstream position. In Sect. 3, the theoretical framework is then validated using RANS simulations and experiments behind a homogeneous grid; the experimental and numerical set-ups are detailed in Sects. 3.1 and 3.2. A validation and the results are presented in Sect. 3.4. In the next step, the digital twin is used to perform aerodynamic simulations which are compared to experiments. The numerical and experimental set-ups are explained in Sects. 4.2 and 4.1, and the results are compared in Sect. 4.3. Finally, conclusions and perspectives are given in Sect. 5.

2 Theoretical framework

In the following, we will detail the theoretical framework of the decay of homogeneous, isotropic turbulence that will lay the foundation of our digital twin. We will start with some important equations from literature and develop them further.

2.1 Dependence of the temporal evolution of k on λ

Researchers in the physics and mathematics community have studied the case of homogeneous isotropic turbulence in great detail. Batchelor and Townsend (1948) proposed to call the Reynolds number based on the Taylor micro-scale (λ) the “Reynolds number of turbulence” and also suggested that λ is representative of the eddies of large wavenumber, i.e. small eddies, before viscosity becomes relevant. The mathematical definition and the experimental methodology for calculating λ are detailed in Sect. 3.3. In the theory of homogeneous isotropic turbulence, it is necessary to assume the similarity of turbulence at all stages of decay; i.e. changes in the structure of turbulence can be described by two parameters, namely a characteristic length and a characteristic velocity (Stewart and Townsend, 1951). Later on, George (1992) proved that the characteristic length scale of the entire energy spectrum is the Taylor micro-scale. This was confirmed by Speziale and Bernard (1992), who also proposed that the turbulent kinetic energy and the Taylor micro-scale are the appropriate scaling parameters for all scales of motion. They proposed a set of differential equations to calculate the temporal decay of turbulent kinetic energy, k(t), and dissipation rate, ϵ(t) :

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}k\left(t\right)}{\mathrm{d}t}}& {\displaystyle }=\mathit{ϵ}\left(t\right),\text{(2)}& {\displaystyle \frac{\mathrm{d}\mathit{ϵ}\left(t\right)}{\mathrm{d}t}}& {\displaystyle }=-\mathit{\alpha }{\displaystyle \frac{\mathit{ϵ}(t{)}^{\mathrm{2}}}{k\left(t\right)}}.\end{array}$$

Here, α is a constant. Equations (1) and (2) formulate an initial value problem which can easily be solved to give a temporal decay law for initial conditions at time t=0, k(t)=k(0), and ϵ(t)=ϵ(0); cf. (Zhou and Speziale, 1998), which gives

$$\begin{array}{}\text{(3)}& k\left(t\right)=k\left(\mathrm{0}\right){\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{\displaystyle \frac{\mathit{ϵ}\left(\mathrm{0}\right)}{k\left(\mathrm{0}\right)}}t\right)}^{-\mathit{\alpha }}.\end{array}$$

With the equation from Bailly and Comte-Bellot (2015),

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{k\left(\mathrm{0}\right)}{\mathit{ϵ}\left(\mathrm{0}\right)}}={\displaystyle \frac{\left(\mathit{\lambda }\right(\mathrm{0}){)}^{\mathrm{2}}}{\mathrm{20}\mathit{\nu }}},\end{array}$$

we can write Eq. (3) as

$$\begin{array}{}\text{(5)}& k\left(t\right)=k\left(\mathrm{0}\right){\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{\displaystyle \frac{\mathrm{20}\mathit{\nu }}{\left(\mathit{\lambda }\right(\mathrm{0}){)}^{\mathrm{2}}}}t\right)}^{-\mathit{\alpha }}.\end{array}$$

Looking at Eq. (5), we can clearly say that the temporal decay of the turbulent kinetic energy has a direct dependence on the initial Taylor micro-scale λ(0). This is comforted by Mydlarski and Warhaft (1996), who empirically found a relationship between the energy decay exponent and the Taylor Reynolds number R_λ.

2.2 Dependence of the spatial evolution of k on λ

The dependence of the downstream evolution of k on λ is demonstrated here within the framework of homogeneous, isotropic turbulence. For the steady-state case, the transport equation for k can be written as (here, we follow Einstein's summation convention)

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{U}_j{\displaystyle \frac{\partial k}{\partial x_j}}=& {\displaystyle \frac{\partial }{\partial x_j}}\left({\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\stackrel{\mathrm{‾}}{u_{j}p}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{u_i{u}_i{u}_j}-\mathrm{2}\mathit{\nu }\stackrel{\mathrm{‾}}{u_i{s}_{ij}}\right)\\ & -\stackrel{\mathrm{‾}}{u_i{u}_j}{S}_{ij}-\mathit{ϵ}.\end{array}\end{array}$$

U_j and u_j are the mean velocity and the velocity fluctuation, respectively, with $i,j=\mathrm{1},\mathrm{2},\mathrm{3}$. The pressure fluctuations are denoted by p, and ν is the kinematic viscosity of the fluid. s_ij is the strain rate tensor for velocity fluctuations, and S_ij is the strain rate tensor for the mean flow, defined as

$$\begin{array}{}\text{(7)}& {\displaystyle }{S}_{ij}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial U_i}{\partial x_j}}+{\displaystyle \frac{\partial U_j}{\partial x_i}}\right),\text{(8)}& {\displaystyle }{s}_{ij}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right).\end{array}$$

The turbulent kinetic energy is defined as

$$\begin{array}{}\text{(9)}& k={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{u_{\mathrm{1}}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{u_{\mathrm{2}}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{u_{\mathrm{3}}^{\mathrm{2}}}\right).\end{array}$$

It should be noted that for grid-generated turbulence

$$\begin{array}{}\text{(10)}& \stackrel{\mathrm{‾}}{u_{\mathrm{1}}^{\mathrm{2}}}=\mathrm{1.2}\stackrel{\mathrm{‾}}{u_{\mathrm{2}}^{\mathrm{2}}}=\mathrm{1.2}\stackrel{\mathrm{‾}}{u_{\mathrm{3}}^{\mathrm{2}}}=\stackrel{\mathrm{‾}}{u^{\mathrm{2}}}.\end{array}$$

(cf. Comte-Bellot and Corrsin, 1966; Bailly and Comte-Bellot, 2015).

After substituting (10) in (9), we get

$$\begin{array}{}\text{(11)}& k={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\stackrel{\mathrm{‾}}{u^{\mathrm{2}}}.\end{array}$$

If the dominant flow direction is the x₁ direction, then Eq. (6) can be written as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{U}_{\mathrm{1}}{\displaystyle \frac{\partial k}{\partial x_{\mathrm{1}}}}=& \underset{\text{transport}}{\underbrace{{\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\left({\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\stackrel{\mathrm{‾}}{u_{\mathrm{1}}p}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{u_i{u}_i{u}_{\mathrm{1}}}-\mathrm{2}\mathit{\nu }\stackrel{\mathrm{‾}}{u_i{s}_{i\mathrm{1}}}\right)}}\\ & -\underset{\text{production}}{\underbrace{\stackrel{\mathrm{‾}}{u_i{u}_{\mathrm{1}}}{S}_{i\mathrm{1}}}}-\mathit{ϵ}.\end{array}\end{array}$$

If we assume the turbulence to be homogeneous and decaying, then the term representing the transport of k in an inhomogeneous field due to pressure fluctuation, the turbulence itself, and viscous stresses, and the term representing the production of k can be set to zero (Bailly and Comte-Bellot, 2015).

Therefore, Eq. (12) becomes

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial k}{\partial x}}=-\mathit{ϵ}/U.\end{array}$$

Here, we have dropped the indices for the sake of simplification. For homogeneous, isotropic turbulence, the dissipation rate ϵ can be related to both the longitudinal Taylor micro-scale (λ₁) and the transverse Taylor micro-scale (λ₂) as

$$\begin{array}{}\text{(14)}& \mathit{ϵ}=\mathrm{15}\mathit{\nu }{\displaystyle \frac{\stackrel{\mathrm{‾}}{u^{\mathrm{2}}}}{{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}}}=\mathrm{30}\mathit{\nu }{\displaystyle \frac{\stackrel{\mathrm{‾}}{u^{\mathrm{2}}}}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}}.\end{array}$$

By using Eq. (11), we can write Eq. (14) as

$$\begin{array}{}\text{(15)}& \mathit{ϵ}=\mathrm{11.25}\mathit{\nu }{\displaystyle \frac{k}{{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}}}=\mathrm{22.5}\mathit{\nu }{\displaystyle \frac{k}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}}.\end{array}$$

From here onwards, we will only use the relation corresponding to λ₁, and for simplicity, we will drop the subscript “1”. Hence, the relation for the dissipation rate ϵ can be written as

$$\begin{array}{}\text{(16)}& \mathit{ϵ}=\mathrm{22.5}\mathit{\nu }{\displaystyle \frac{k}{{\mathit{\lambda }}^{\mathrm{2}}}}.\end{array}$$

By substituting ϵ from Eq. (16) into Eq. (13), we arrive at

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\partial k}{\partial x}}=-{\displaystyle \frac{\mathrm{22.5}\mathit{\nu }}{U}}{\displaystyle \frac{k}{{\mathit{\lambda }}^{\mathrm{2}}}}.\end{array}$$

Since the evolution of k is only a function of x, for the given case, the partial differential Eq. (17) becomes an ordinary differential equation,

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}k}{\mathrm{d}x}}=-{\displaystyle \frac{\mathrm{22.5}\mathit{\nu }}{U}}{\displaystyle \frac{k}{{\mathit{\lambda }}^{\mathrm{2}}}}.\end{array}$$

Solving the differential Eq. (18) ,

$$\begin{array}{}\text{(19)}& k\left(x\right)=C_{\mathrm{1}}\text{exp}\left(-{\displaystyle \frac{\mathrm{22.5}\mathit{\nu }}{U}}\int {\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}^{\mathrm{2}}}}\mathrm{d}x\right).\end{array}$$

Equation (19) shows that the evolution of k in the downstream direction has a direct dependence on the Taylor micro-scale λ. Here, C₁ is a constant.

To find the final solution, we need to understand the dependence of λ on x. This is performed in the following section.

2.3 Dependence of λ on x

From Eq. (16) we know that ${\mathit{\lambda }}^{\mathrm{2}}\propto k/\mathit{ϵ}$. Experiments have shown that for homogeneous isotropic turbulence, $k/\mathit{ϵ}$ evolves linearly in the downstream direction (e.g. Comte-Bellot and Corrsin, 1966). Thus, we can write

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{k}{\mathit{ϵ}}}={\displaystyle \frac{k_{\text{in}}}{{\mathit{ϵ}}_{\text{in}}}}+{\displaystyle \frac{m}{U}}x,\end{array}$$

where $m/U$ is the slope. Now, using Eqs. (16) and (20), we can write

$$\begin{array}{}\text{(21)}& {\mathit{\lambda }}^{\mathrm{2}}=\mathrm{22.5}\mathit{\nu }\left({\displaystyle \frac{k}{\mathit{ϵ}}}\right)=\mathrm{22.5}\mathit{\nu }\left({\displaystyle \frac{k_{\text{in}}}{{\mathit{ϵ}}_{\text{in}}}}+{\displaystyle \frac{m}{U}}x\right).\end{array}$$

k_in and ϵ_in are the values of the turbulent kinetic energy and the dissipation rate at the starting point x=0 from where the TKE's decay is calculated. Substituting Eq. (21) into Eq. (19), performing integration, and subjecting it to the boundary condition at the starting point, $k(x=\mathrm{0})=k_{\text{in}}$, gives

$$\begin{array}{}\text{(22)}& k\left(x\right)=k_{\text{in}}{\left(\mathrm{1}+m{\displaystyle \frac{{\mathit{ϵ}}_{\text{in}}}{k_{\text{in}}}}{\displaystyle \frac{x}{U}}\right)}^{-\mathrm{1}/m}.\end{array}$$

Using relation (Eq. 16), we can rewrite Eq. (22) as

$$\begin{array}{}\text{(23)}& k\left(x\right)=k_{\text{in}}{\left(\mathrm{1}+m{\displaystyle \frac{\mathrm{22.5}\mathit{\nu }}{{\mathit{\lambda }}_{\text{in}}^{\mathrm{2}}}}{\displaystyle \frac{x}{U}}\right)}^{-\mathrm{1}/m},\end{array}$$

where $\mathit{\lambda }(x=\mathrm{0})={\mathit{\lambda }}_{\text{in}}$. It should be noted that for Eq. (23), the starting point, i.e. x=0, and, thus, the origin of the coordinate system, lies at the point of measurement which is used to give the value of k_in and λ_in. Looking at Eq. (23), we can say that the evolution of k in the downstream direction can be defined for a homogeneous, isotropic flow, provided we have the values of k_in and λ_in at one point.

From the framework of the k−ω series of models for RANS simulations, one can also derive the following evolution equation for k for decaying, homogeneous, isotropic turbulence (Eça et al., 2016),

$$\begin{array}{}\text{(24)}& k\left(x\right)=k_{\text{in}}{\left(\mathrm{1}+{\mathit{\omega }}_{\text{in}}\mathit{\beta }{\displaystyle \frac{x}{U}}\right)}^{-{\mathit{\beta }}^{\ast }/\mathit{\beta }},\end{array}$$

where ω is the so-called turbulence frequency that is related to ϵ and k by

$$\begin{array}{}\text{(25)}& \mathit{\omega }={\displaystyle \frac{\mathit{ϵ}}{{\mathit{\beta }}^{\ast }k}},\end{array}$$

and β and β^∗ are constants with values of 0.0828 and 0.09, respectively (Wilcox, 1988). The value of ω at the inlet boundary is given by ω_in. By substituting Eq. (25) into Eq. (24), we have

$$\begin{array}{}\text{(26)}& k\left(x\right)=k_{\text{in}}{\left(\mathrm{1}+{\displaystyle \frac{\mathit{\beta }}{{\mathit{\beta }}^{\ast }}}{\displaystyle \frac{{\mathit{ϵ}}_{\text{in}}}{k_{\text{in}}}}{\displaystyle \frac{x}{U}}\right)}^{-{\mathit{\beta }}^{\ast }/\mathit{\beta }}.\end{array}$$

Using the relation in Eq. (16), we can rewrite Eq. (27) as follows:

$$\begin{array}{}\text{(27)}& k\left(x\right)=k_{\text{in}}{\left(\mathrm{1}+{\displaystyle \frac{\mathit{\beta }}{{\mathit{\beta }}^{\ast }}}{\displaystyle \frac{\mathrm{22.5}\mathit{\nu }}{{\mathit{\lambda }}_{\text{in}}^{\mathrm{2}}}}{\displaystyle \frac{x}{U}}\right)}^{-{\mathit{\beta }}^{\ast }/\mathit{\beta }}.\end{array}$$

Equations (23) and (27) have a similar form. Even if Eq. (27) comes from the k−ω model, believing in its proven applicability, we can assume that

$$\begin{array}{}\text{(28)}& m={\displaystyle \frac{\mathit{\beta }}{{\mathit{\beta }}^{\ast }}}.\end{array}$$

Therefore, the value of m can be taken as 0.92. It is important to notice that further investigations are required to find the correct value of the parameter m theoretically. The authors wish to emphasize to the readers that, unlike the equations commonly encountered in prior literature, such as those referenced in Comte-Bellot and Corrsin (1966), Kurian and Fransson (2009), Krogstad and Davidson (2009), or Sinhuber et al. (2015), Eq. (23) does not have any fitting parameter, and it is neither an empirical equation nor does it have any virtual origin. The only assumption taken while deriving Eq. (23) is statistical stationarity of fully developed GGT, which is believed to start from $x/M\approx \mathrm{20}$ downstream of a grid (Comte-Bellot and Corrsin, 1966; Bailly and Comte-Bellot, 2015). Upstream of $x/M\approx \mathrm{20}$, one may expect some changes in the form of Eqs. (10), (13), and (20), which are used to derive Eq. (23). However, between $x/M\approx \mathrm{10}$ and $x/M\approx \mathrm{20}$, the turbulent flow can be considered approximately developed (Frisch, 1995). Therefore, for all practical purposes, Eq. (23) is valid downstream of $x/M\approx \mathrm{10}$. For the cases where the assumption of Taylor's hypothesis is valid, i.e. time t equals $x/U$, the spatial evolution equation of k, Eq. (23), becomes equivalent to the temporal decay equation of k (cf. Eq. 5). Hence, measurements of k and λ at a given position in a regular-grid-generated turbulent flow will enable us to obtain the evolution of k when solving the RANS equations. This is essential to reproduce realistic regular-grid-generated turbulent inflow conditions at the inlet boundary of RANS simulations. The experimental and numerical set-ups used to validate the spatial evolution in Eq. (23) are presented in the following section.

3 Digital twin: simulating regular grid inflow in the wind tunnel

In the following, k and λ will be measured from the flow behind a grid in a wind tunnel facility. The objective is to demonstrate that the evolution of k obtained from RANS simulations matches the experiments when using k and λ from measurements at a given location. The wind tunnel facility and the grid set-up are described below together with the RANS simulations performed.

Figure 1 Schematic diagram of the experimental set-up for hot-wire measurements in the wake of a regular grid (side view). The measurement positions are marked.

3.1 Experimental set-up

The experiments were performed in the aerodynamic closed-loop low-Reynolds-number wind tunnel facility at the LHEEA (Laboratoire de recherche en Hydrodynamique, Énergétique et Environnment Atmosphérique) laboratory of CNRS and Centrale Nantes (see Fig. 1). This wind tunnel has a cross section of 0.5 m × 0.5 m and a test section length of 2.3 m with a maximum inflow velocity of 40 m s⁻¹ and a turbulence intensity of less than 0.3 %.

To induce a turbulent inflow, a regular wooden grid with square bars was used. The cross section of the bars used for the frame is 11 mm × 10.5 mm and of the ones used for creating the mesh d = 6 mm × 6 mm. The blockage (b) of the grid is 16 %, and the grid mesh size M is 70 mm × 70 mm.

To measure the downstream evolution of turbulence properties of the inflow, a 1D hot-wire probe of type 55P11 with a wire length of 1.25 mm from Dantec Dynamics was used. It was operated using a DISA55M01 unit. The hot-wire was calibrated in the velocity range 0.5 m s⁻¹ $\le U\le$ 40 m s⁻¹ applying the temperature correction suggested by Hultmark and Smits (2010). During the measurements, the mean velocity was approximately 25 m s⁻¹. The data were recorded at a sampling frequency of 25 kHz for 10 s for calibration and 20 s for the measurements. A hardware low-pass filter with a cut-off frequency of 10 kHz was used. The downstream evolution of k was measured at the centre line of the wind tunnel at nine downstream positions that are indicated in Fig. 1.

3.2 Numerical set-up

The simulations are performed using the finite-volume-based ISIS-CFD incompressible unsteady Reynolds-averaged Navier–Stokes (URANS) solver. This solver, developed by CNRS and Centrale Nantes, also available as a part of the FINE™/Marine computing suite worldwide distributed by Cadence Design Systems, uses an incompressible URANS method. The solver is based on a finite volume method to build the spatial discretization of the transport equations. The unstructured discretization is face-based, which means that cells with an arbitrary number of arbitrarily shaped faces are accepted. A second-order backward difference scheme is used to discretize time. All flow variables are stored at the geometric centre of arbitrarily shaped cells. Volume and surface integrals are evaluated with second-order accurate approximations. As the method is face-based, numerical fluxes are reconstructed on the mesh faces by linear extrapolation of the integrand from the neighbouring cell centres. A centred scheme is used for the diffusion terms, whereas for the convective fluxes, a blended scheme with 80 % central and 20 % upwind is used. In the case of turbulent flows, additional transport equations for the variables in the turbulence model are added. In the following, the transport equations for URANS simulations are presented:

The momentum conservation equation is

$$\begin{array}{}\text{(29)}& \mathit{\rho }{\displaystyle \frac{\partial \mathit{U}}{\partial t}}+\mathit{\rho }(\mathit{U}\cdot \mathbf{\nabla }\mathit{U})=\mathit{\rho }\mathit{g}-\mathbf{\nabla }\stackrel{\mathrm{‾}}{p}+\mathit{\mu }\mathrm{△}\mathit{U}-\mathit{\rho }\mathbf{\nabla }\cdot \stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{R}}.\end{array}$$

The continuity equation for the mean component is

$$\begin{array}{}\text{(30)}& \mathbf{\nabla }\cdot \mathit{U}=\mathrm{0}.\end{array}$$

The continuity equation for fluctuations is

$$\begin{array}{}\text{(31)}& \mathbf{\nabla }\cdot \mathit{u}=\mathrm{0}.\end{array}$$

The Reynolds stress tensor $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{R}}$ present in Eq. (29) is defined as

$$\begin{array}{}\text{(32)}& \stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{R}}={\mathit{\mu }}_{\mathrm{t}}\left(\stackrel{\mathrm{‾}}{\mathit{u}\otimes \mathit{u}}\right),\end{array}$$

where U and u are the mean and fluctuating component of the velocity, and μ_t is the turbulent viscosity. There are many ways to calculate the value of μ_t, and in this paper, a two-equation eddy viscosity method, namely k−ω SST Menter (2003), has been used to do so. In the following, a description of this model is given.

The transport equation for the turbulent kinetic energy k is

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(\mathit{\rho }k\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\mathit{\rho }{U}_{i}k\right)}{\partial x_i}}=& \stackrel{\mathrm{̃}}{P_k}-{\mathit{\beta }}^{\ast }\mathit{\rho }k\mathit{\omega }\\ & +{\displaystyle \frac{\partial }{\partial x_i}}\left[(\mathit{\mu }+{\mathit{\sigma }}_k{\mathit{\mu }}_{\mathrm{t}}){\displaystyle \frac{\partial k}{\partial x_i}}\right].\end{array}\end{array}$$

The transport equation for the turbulence frequency ω is

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(\mathit{\rho }\mathit{\omega }\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\mathit{\rho }{\stackrel{\mathrm{‾}}{u}}_i\mathit{\omega }\right)}{\partial x_i}}=& \mathit{\alpha }\mathit{\rho }{S}^{\mathrm{2}}-\mathit{\beta }\mathit{\rho }{\mathit{\omega }}^{\mathrm{2}}\\ & +{\displaystyle \frac{\partial }{\partial x_i}}\left[(\mathit{\mu }+{\mathit{\sigma }}_{\mathit{\omega }}{\mathit{\mu }}_{\mathrm{t}}){\displaystyle \frac{\partial \mathit{\omega }}{\partial x_i}}\right]\\ & +\mathrm{2}(\mathrm{1}-F_{\mathrm{1}})\mathit{\rho }{\mathit{\sigma }}_{w\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }}}{\displaystyle \frac{\partial k}{\partial x_i}}{\displaystyle \frac{\partial \mathit{\omega }}{\partial x_i}}.\end{array}\end{array}$$

F₁ is the blending function. F₁ goes to zero for the flow away from the surface. Hence the k−ϵ model is applied there, and it goes to one near the wall where the k−ω model is applied. The constant α is computed by $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{1}}{F}_{\mathrm{1}}+{\mathit{\alpha }}_{\mathrm{2}}(\mathrm{1}-F_{\mathrm{1}})$. For this model, the values of constant are ${\mathit{\beta }}^{\ast }=\mathrm{0.09}$, ${\mathit{\alpha }}_{\mathrm{1}}=\mathrm{5}/\mathrm{9}$, ${\mathit{\beta }}_{\mathrm{1}}=\mathrm{3}/\mathrm{40}$, σ_k1=0.85, σ_ω1=0.5, α₂=0.44, β=0.0828, σ_k2=1, and σ_ω2=0.856. Using these constants, both transport equations are solved and the turbulent viscosity is found by Eq. (35).

$$\begin{array}{}\text{(35)}& {\mathit{\mu }}_{\mathrm{t}}={\displaystyle \frac{\mathit{\rho }{a}_{\mathrm{1}}k}{\text{max}(a_{\mathrm{1}}\mathit{\omega },S{F}_{\mathrm{2}})}},\end{array}$$

where S is the measure of the strain rate and F₂ is the second blending function defined by

$$\begin{array}{}\text{(36)}& F_{\mathrm{2}}=\tanh \left[{\left[\text{max}\left({\displaystyle \frac{\mathrm{2}\sqrt{k}}{{\mathit{\beta }}^{\ast }\mathit{\omega }y}},{\displaystyle \frac{\mathrm{500}\mathit{\nu }}{y^{\mathrm{2}}\mathit{\omega }}}\right)\right]}^{\mathrm{2}}\right].\end{array}$$

Figure 2 Simulation domain for decaying turbulence.

Figure 3 Power spectrum obtained from hot-wire measurements of the flow downstream of the regular grid in the frequency domain for the position $x/M=\mathrm{11}$.

Using the above-presented set of equations, 2D simulations are performed for the domain, mimicking the wind tunnel shown in Fig. 2, with a domain length of 3.12 m, a width of 0.5 m (same as the width of the wind tunnel test section), and 8000 cells. Simulations were performed with a higher number of cells (10000 and 12000) as well, but no changes were observed. The top and the bottom walls in the simulations were put to the no-slip condition.

It should be noted that the values put as the boundary condition at the inlet boundary are those obtained at the first point of the measurement, 758 mm downstream of the wind tunnel inlet, which refers to $x/M=\mathrm{11}$ in Figs. 4 and 6.

3.3 Length scales

For homogeneous, isotropic turbulence, both the integral length (L) and λ can be easily estimated using the one-dimensional energy spectrum (Hinze, 1975). L can be obtained by taking the limit of the energy spectrum E(f) in the frequency domain for f→0,

$$\begin{array}{}\text{(37)}& L=\underset{f\to \mathrm{0}}{lim}\left({\displaystyle \frac{E\left(f\right)\cdot U}{\mathrm{4}{\mathit{\sigma }}^{\mathrm{2}}}}\right).\end{array}$$

Here, U is the mean streamwise velocity, u denotes the fluctuations, and σ is the standard deviation of the velocity time series U(t). It should be noticed here that only the frequency range where E(f)≈const just outside of the inertial sub-range is used to determine L.

λ is defined as follows:

$$\begin{array}{}\text{(38)}& \mathit{\lambda }={\left({\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}}{〈{\left(\frac{\partial u}{\partial x}\right)}^{\mathrm{2}}〉}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where $\langle (\frac{\partial u}{\partial x}{)}^{\mathrm{2}}\rangle$ can be determined from the spectrum in the wavenumber (κ) domain, derived through hot-wire measurements using Taylor's hypothesis,

$$\begin{array}{}\text{(39)}& \langle {\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}\rangle =\underset{{\mathit{\kappa }}_{\text{min}}}{\overset{{\mathit{\kappa }}_{\text{max}}}{\int }}{\mathit{\kappa }}^{\mathrm{2}}E\left(\mathit{\kappa }\right)\mathrm{d}\mathit{\kappa },\end{array}$$

where κ_min and κ_max are the wavenumber boundaries of the whole energy spectrum. At the downstream position $x/M=\mathrm{11}$, which represents the first point of measurement in the wind tunnel, we calculated $L/M$, $\mathit{\lambda }/M$, and the Taylor Reynolds number (Re_λ) using hot-wire measurements and Eqs. (37) and (38). An exemplary spectrum is shown in Fig. 3. The calculated values were $L/M=\mathrm{0.39}$, $\mathit{\lambda }/M=\mathrm{0.034}$, and Re_λ=180, respectively.

Figure 4 Comparison of the downstream evolution of $k/\left(\mathrm{0.5}{U}^{\mathrm{2}}\right)$ obtained from theory with experiments and simulations performed with integral length and the Taylor micro-scale as the inlet boundary condition using the k−ω SST Menter 2003 model. (a) Linear and (b) logarithmic.

Figure 5 Comparison of experimental decay of normalized TKE decay obtained from experiments performed in the University of Oldenburg and data from Batchelor and Townsend (1948) (B&T) with Eq. (23). For the application of Eq. (23) to the B&T data, we estimate λ≈ 9 mm from their article.

3.4 Validation

Figure 4 shows the downstream evolution of k obtained theoretically with Eq. (23), experimentally using hot-wire measurements (Sect. 3.1) and from simulation using the k−ω SST Menter 2003 model (see Sect. 3.2) with both L and λ as boundary conditions at the inlet. First, we see that the theoretical equation is validated against experiments. The log–log equivalent of the curve is provided in the inset. Next, simulation results obtained using both the Taylor micro-scale and the integral length as the boundary condition at the inlet are compared with Eq. (23). Here, we can clearly see that the simulation result obtained for the case where the Taylor micro-scale is used as the boundary condition matches very well with the theory. The TKE decay exponent was determined to be 1.087 in Eq. (23) and simulation, while in the experimental data, it was observed to be 1.09. These values are in close proximity to the Saffman decay exponent of 1.1 for grid-generated turbulence, as reported by Krogstad and Davidson (2009). In contrast, when the integral length is used as a boundary condition, the simulation results do not match the theory or the experimental results. The justification for this observation comes directly from Eq. (23) where the derived 1D spatial evolution equation shows a direct dependence of the turbulent kinetic energy on the Taylor micro-scale.

Figure 6 Comparison of the downstream evolution of $k/\left(\mathrm{0.5}{U}^{\mathrm{2}}\right)$ computed using different k−ω models.

To verify the generality of Eq. (23), apart from the validation performed against the TKE decay data obtained in the LHEEA wind tunnel, we also compared it with the results from hot-wire experiments conducted independently in the wind tunnel at the University of Oldenburg and data given in Batchelor and Townsend (1948). The Oldenburg experiments were performed using a passive regular grid with M=115 mm for 33 downstream positions spanning from $x/M\sim \mathrm{8}$ to $x/M\sim \mathrm{170}$ for two inflow speeds: 5 m s⁻¹ and 10 m s⁻¹. Figure 5 shows the log–log plot of the comparison of the TKE decay obtained experimentally with that from Eq. (23). It can clearly be seen that the evolution of TKE given by Eq. (23) matches very well with the experimental data. Note that the deviations are over-accentuated when visualized in the log–log plot. Readers interested in knowing the details of the experiments performed at the University of Oldenburg may refer to Appendix A.

It is also interesting to see how different k−ω models calculate the downstream evolution of k by using λ as the boundary condition. This is shown in Fig. 6. Here, we can see that the results from the k−ω SST Menter 2003, k−ω SST Menter 1994, and k−ω baseline (BSL) Menter overlap with only a little deviation for the oldest k−ω series of the turbulence models, the k−ω Wilcox model.

This emphasizes the very important role of the choice of the turbulent length scale at the simulation domain inlet as compared to the choice of the turbulent model. The following section showcases an instance of the digital twin's functionality, where an airfoil is subjected to testing in both the physical wind tunnel and the digital twin, with a subsequent comparison of the resulting loads.

4 Example of the performance of the digital twin: testing an airfoil section of a wind turbine

The methodology to obtain a digital twin model of the experimented turbulent inflows is now validated. The present section is focusing on modelling the wind tunnel experiments to reproduce loads when a turbulent inflow is set. The same turbulent inflow as described in Sect. 3.1 is used, which induces a homogeneous field with a turbulence intensity of 3 %, measured at the airfoil position before its installation.

We are then introducing an airfoil model described in Sect. 4.1.1, which is equipped with global load sensors (Sect. 4.1.2) and local pressure sensors (Sect. 4.1.3).

Figure 7 Blade section at 82 % of the radius in comparison with a NACA63-3-420 profile with a modified camber of 4 % instead of 2 % (from Neunaber et al., 2022; published under CC BY-NC-ND 4.0, , ).

Figure 8 Computer-aided design (CAD) drawing of the 3D printed 2D blade section mounted in the wind tunnel. The 2D blade section is made of aluminium. It is a hollow model that contains four pressure scanners with 16 channels each, which are connected to integrated pressure taps. The 2D blade section has four removable covers on the suction side to grant access to the pressure scanners.

4.1 Experimental set-up

4.1.1 Airfoil

A 2D blade section was placed 1.79 m downstream of the wind tunnel inlet where the flow field is homogeneous (Neunaber and Braud, 2020b) in the wind tunnel described in Sect. 3.1. The inflow velocity is 25 m s⁻¹. The airfoil shape was derived from scans of a 2 MW wind turbine blade section, at 82 % of its length (Neunaber et al., 2022). It has been scaled down to 1/10th of the original chord length, so that the chord length is c=0.125 m and the chord-based Reynolds number is ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$. The airfoil section closely resembles a NACA63-3-420 profile with a modified camber of 4 % instead of 2 % (see Fig. 7).

Figure 9 Illustration of the deviations (in mm) between the airfoil's original design shape and the shape achieved after the manufacturing process.

The 2D blade section has been designed with multiple sensors to perform a 3D characterization of the wall pressure over the airfoil surface, and future actuators and/or sensors can be implemented on the suction side (see Fig. 8). It was manufactured from aluminium using a 3D metal printer to integrate channels for the pressure measurements. In total, four pressure scanners with 16 channels each are used, and the locations of the pressure taps are given in Sect. 4.1.2. The model is hollow and equipped with four covers on the suction side for access to the sensors. This also allows for the integration of actuators in the future. To perform simulations with the digital twin, it was important to check the shape of the airfoil for deviations and unevenness after manufacturing. Therefore, the down-scaled model has been scanned using the HandyScan 3D 700^MC from CREAFORM, which has an accuracy of 0.03 mm. The scan was then compared to the initially designed shape that is used in simulations (see Fig. 9). At first, the covers of the 2D blade section had steps as high as 0.7 mm, which was not accurate enough to produce experimental results that matched the simulations. These steps were significantly improved manually to an accuracy of less than 0.45 mm, which was sufficient to match simulation results. It should be noted that for this inflow velocity, a TI of at least 3 % was necessary to avoid low Reynolds number effects found in previous investigations (see Mishra et al., 2022).

Table 1 Position of pressure ports.

4.1.2 Local pressure sensors

The blade section is equipped with four differential pressure scanners from EvoScann® (P-Series) with 16 channels each that have a range of ± 50 mbar. The acquisition rate can be up to 1 kHz, but it was limited to 100 Hz for the present measurements as only steady quantities were targeted. These sensors were connected to wall pressure taps through channels integrated in the blade designs and Tygon tubes. Three chord-wise lines of pressure taps were distributed (between the covers): at the mid-span, $z/c=\mathrm{0}$, and at $z/c=\pm \mathrm{1}$. The chord-wise distribution is identical for the three lines, and it is shown in Fig. 10 with exact positions in Table 1. One span-wise line of pressure taps was added at $x/c=\mathrm{0.88}$, with pressure taps spaced by 0.16c from $z/c=-\mathrm{1}$ to $z/c=\mathrm{1}$ and 0.25c otherwise.

Figure 10 Chord-wise distribution of pressure ports around the airfoil – note that the axes are not scaled correctly.

Figure 11 Wind tunnel set-up for load measurements.

4.1.3 Global load sensors

The airfoil was supported on two sides by load cells (cf. Fig. 11). The load cells work by the principle of strain measurement. Two load cells (one on each side) were used to measure the lift force (F_l), and two were used to measure the drag force (F_d). The load measurement system was calibrated in all directions (i.e. +x, −x, +y, and −y) using calibration weights between 500 and 5000 g. The angle of attack was measured using a high-precision voltage-based angle sensor with a resolution of 0.035 ${}^{\circ }{\mathrm{mV}}^{-\mathrm{1}}$. The signal from each load cell was collected at a sampling frequency of f_s = 5000 Hz. The lift coefficient (C_l) and the drag coefficient (C_d) were calculated using Eqs. (40) and (41):

$$\begin{array}{}\text{(40)}& {\displaystyle }{C}_{\mathrm{l}}={\displaystyle \frac{F_{\mathrm{l}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{U}^{\mathrm{2}}}},\text{(41)}& {\displaystyle }{C}_{\mathrm{d}}={\displaystyle \frac{F_{\mathrm{d}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{U}^{\mathrm{2}}}},\end{array}$$

where ρ is the density of air, A = 0.125 m × 0.5 m is the model planform area, U is the inflow speed, and ν is kinematic viscosity.

Figure 12 Views of the mesh: (a) 3D mesh for the simulation of the flow over the airfoil (3 million cells), (b) side view of the mesh, (c) mesh around the airfoil, (d) front view of the mesh, and (e) side view of the mesh at the end of simulation after adaptive grid refinement (16 million cells).

4.2 Numerical set-up for airfoil simulations

Figure 12a, b, and d present the 3D view, front view, and side view of the numerical domain (including mesh), respectively, for an angle of attack of the 2D blade section of α = 0^∘. The mesh consists of 3 million cells, and Fig. 12c displays a close-up of the mesh around the airfoil. In this study, we only simulated the transverse half of the wind tunnel, as shown in Fig. 12a. The shaded area in the figure was not simulated; instead, we applied a mirror boundary condition. The chord length of the airfoil is 0.125 m, which is the same as that used in the experiments. The airfoil is positioned at a downstream distance of 1.035 m from the domain's inlet, as the simulation domain starts at the position of the first measurement point, where the inlet conditions are obtained and from where the turbulence decays. Therefore, the downstream position of the airfoil in the simulation is 0.758 m less than that in the experiments, resulting in a downstream position of 1.035 m.

We have applied a Dirichlet boundary condition at the inlet, and the values are given in Table 2. These values correspond to the values obtained at $x/M=\mathrm{11}$ from the hot-wire measurement (see Sect. 3.4). For pressure, we applied the Neumann boundary condition at the inlet, $\frac{\mathrm{d}p}{\mathrm{d}\mathit{n}}=\mathrm{0}$, where n is the normal vector to the inlet. These same values have been used as the initial conditions as well. In addition, we also use the integral length scale (L=25 mm) to investigate the impact of using the “wrong” length scale at the simulation domain inlet. At the outlet, the velocity is found using Rhie–Chow interpolation. We applied the Dirichlet boundary condition at the outlet for pressure p=p_o, where p_o=0 by default. For TKE and turbulence frequency, we applied the Neumann boundary condition as $\frac{\mathrm{d}k}{\mathrm{d}n}=\mathrm{0}$ and $\frac{\mathrm{d}\mathit{\omega }}{\mathrm{d}n}=\mathrm{0}$, respectively. We have applied a no-slip boundary condition on the airfoil and imposed wall functions

$$\begin{array}{}\text{(42)}& {\displaystyle \frac{\partial U}{\partial y}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{s}}}{\mathit{\kappa }\mathit{\rho }{c}_{\mathit{\mu }}^{\mathrm{1}/\mathrm{4}}\sqrt{k_{\mathrm{w}}}{y}_{\mathrm{w}}}},\end{array}$$

on the top wall (TW), bottom wall (BW), and side wall (SW) to avoid explicitly simulating the boundary layer. Here, U is the velocity, k_w is the TKE at the cell centre of the first cell from the wall, y_w is the perpendicular distance of the cell centre of the first cell from the wall, τ_s is the wall shear stress, κ=0.41, and c_μ=0.09. The y+ values for airfoil, BW, TW, and SW are given in the Table 3.

Table 2 Boundary conditions at the simulation domain inlet. Note that for the inlet length scale L_S, the Taylor micro-scale is used.

Table 3 y⁺ values for the simulation.

We conducted unsteady 3D RANS simulations using the k−ω SST Menter 2003 model (Menter et al., 2003) following the same procedure and using the same equations that we presented in Sect. 3.2. We performed these simulations using a standard AVLSMART numerical scheme. To improve the computational efficiency and to accurately capture the flow downstream of the airfoil, we utilized an in-house adaptive grid refinement (AGR) methodology (Wackers et al., 2012). The total number of cells at the end of the simulation was approximately 16 million (see Fig. 12e).

To obtain well-converged results for each angle of attack (AoA), we conducted 3D simulations using 400 cores on the IDRIS Jeay Supercomputer. This process took almost 100 computing hours for each AoA.

4.3 Comparison between experiments and simulations

In the following sections, we present a comparison between the force coefficients and pressure coefficients obtained from the experiments and simulations.

Figure 13 Comparison of C_l curves obtained from a digital twin (using the Taylor micro-scale (λ) and integral length (L)) and experiments for ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$. Error bars presented in the figure are the standard deviation of the time series data.

Figure 14 Comparison of C_d curves obtained from a digital twin (using Taylor micro-scale, λ, and integral length, L) and experiments for ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$. Error bars presented in the figure are the standard deviation of the time series data.

4.3.1 Comparison of force coefficients obtained from a digital twin and experiments

Figures 13 and 14 present a comparison between the C_l and C_d curves derived from wind tunnel experiments and their digital counterpart for a Reynolds number of 2.0×10⁵. In one case, the inlet length scale for the simulation domain is defined by λ, while in the other case, it is determined by the integral length scale. The experimental tests covered angles of attack (AoAs) ranging from −5 to 24^∘, whereas the simulations were conducted for AoAs between 0 and 20^∘. The figures clearly demonstrate that our digital twin, utilizing λ as the inlet length scale for the simulation domain, successfully replicates the force coefficients obtained in the original wind tunnel experiments under turbulent inflow conditions. In contrast, when the integral length scale is used as the inlet length scale, higher lift values are observed across all AoAs (with the exception of 0^∘), with the disparity increasing as the AoA increases. This increase in lift coefficient values in the simulations, utilizing L as the boundary condition, can be attributed to the greater turbulence intensity experienced by the airfoil in the simulations compared to the experiments (cf. Abbott and Von Doenhoff, 2012). The inadequate representation of the evolution of the turbulent kinetic energy in the simulations subsequently affects the representation of the turbulence intensity, resulting in notable deviations in the lift coefficients (C_l) between the simulations and experimental data.

In the next subsection, a comparison of the C_p obtained from simulations and experiments is made.

Figure 15 Comparative analysis of the C_p curve for AoA = 0^∘ and ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$, considering 3D simulations and experimental data. The comparison is conducted by averaging the results obtained from three span-wise positions (z=0 mm, $z=\pm \mathrm{125}$ mm).

Figure 16 Comparative analysis of the C_p curve for AoA = 4^∘ and ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$, considering 3D simulations and experimental data. The comparison is conducted by averaging the results obtained from three span-wise positions (z=0 mm, $z=\pm \mathrm{125}$ mm).

Figure 17 Comparative analysis of the C_p curve for AoA = 12^∘ and ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$, considering 3D simulations and experimental data. The comparison is conducted by averaging the results obtained from three span-wise positions (z=0 mm, $z=\pm \mathrm{125}$ mm).

4.3.2 Comparison of C_p

Pressure measurements were performed over the surface of the airfoil, both experimentally and in the digital twin, for different angles of attack (AoAs). Specifically, we compared the C_p values at AoAs of 0, 4, and 12^∘, encompassing a broad range of angles. For that, we average over three span-wise positions (z=0 mm, $z=\pm \mathrm{125}$ mm) and plot the average. Figures 15–17 depict the comparison between experimental data and 3D simulation results. In these comparisons, we utilized both the Taylor length scale and the integral length scale as inlet conditions for the simulation domain.

Overall, the simulation results closely align with the experimental results. However, a slight tendency toward higher pressure on the pressure side is evident. The exact cause of this variation remains unknown at present; differences in the extraction of the reference pressure and the dynamic pressure in simulations and experiments might contribute to these variations.

Moreover, the differences in C_p levels between using the Taylor length scale and the integral length scale become more pronounced as the angle of incidence increases. These differences primarily manifest on the suction side of the airfoil, particularly in the region of the leading-edge suction peak.

Figure 18 Comparison of the C_l curve obtained from a 3D digital twin and a 2D digital twin for ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$.

Figure 19 Comparison of the C_d curve obtained from a 3D digital twin and a 2D digital twin for ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$.

4.3.3 Comparison of the performance of 2D digital twin against 3D digital twin

This section provides a comparative study of force coefficients derived from both 2D and 3D digital twin simulations. Consistent initial and boundary conditions were maintained in both models. Figures 18 and 19 visually display the comparison for the lift and drag coefficient, respectively. The 2D outcomes are graphed for angles of attack between -5^∘ and 16^∘, while the 3D results cover 0 to 16^∘.

The findings show a considerable correlation in the C_l values until an AoA of 12^∘. The 2D digital twin exhibits higher C_l values compared to the 3D digital twin for angles of attacks (AoAs) of 14^∘ and above. The C_d values show a close match across all AoAs, with the exception of 16^∘ AoA, at which point the 2D digital twin displays a higher C_d than its 3D counterpart.

The difference in C_l values between the 3D and 2D digital twins can be attributed to the growing importance of 3D effects, for example flow bi-stability, which impacts the position of flow separation at and above 14^∘ AoAs in high Reynolds number experiments at the same airfoil in Neunaber et al. (2022), something which cannot be reproduced in the 2D digital twin. As a result, if the emphasis is purely on force coefficients, the 2D digital twin can be used for low to moderately high AoAs due to its computational efficiency. For higher AoAs, where 3D effects are substantial, the more computationally demanding 3D digital twin is preferable. This strategic blend optimizes the use of the digital twin framework.

5 Conclusions

We have successfully created a digital twin of a wind tunnel with turbulent inflow conditions after creating a theoretical framework for the decay of the TKE, and we successfully expanded this digital twin to perform 3D simulations of the impact of turbulent inflow on a 2D blade section. In the first part, a theoretical proof of the dependence of the downstream evolution of the turbulent kinetic energy on the Taylor micro-scale was developed for regular-grid-generated turbulent flows. It has then been demonstrated that the Taylor micro-scale is the correct turbulent length scale to be used as the boundary condition at the inlet in RANS simulations for the accurate prediction of the downstream evolution of the turbulent kinetic energy. To validate our results, we compared the theoretical development with measurements performed downstream of a regular grid and RANS simulations using the k−ω SST Menter 2003 model. When the Taylor micro-scale (measured experimentally) is used as the inlet condition, both for RANS simulations and for the starting point of the theoretical equations, the turbulent kinetic energy measured experimentally is retrieved. In contrast, when the integral length scale is used as the inlet boundary condition, RANS simulations are far from experiments and theory. Further, we compare the RANS results using several k−ω models with the Taylor micro-scale as the boundary condition, and the results are in good agreement with each other. This work thus demonstrates the validity of closure models in RANS equations to describe homogeneous, isotropic flows, as long as the Taylor micro-scale is used as one inlet boundary condition instead of the integral length scale. Further, our results emphasize fundamental behaviours of grid-generated turbulent flows:

1. The spatial evolution of the turbulent kinetic energy has a fundamental dependence on the Taylor micro-scale (see Eq. 23).

2. For properly capturing the evolution of the turbulent kinetic energy either in space or time, the correct length scale given as the boundary condition is the Taylor micro-scale.

In the third part, we introduced an airfoil into the numerical wind tunnel and conducted 3D simulations at a chord-based Reynolds number of ${\mathit{Re}}_{\mathrm{c}}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{5}}$. Our findings showed a nearly perfect agreement between the force coefficients obtained from experiments in the physical wind tunnel and those obtained from simulations in the numerical wind tunnel when using the Taylor micro-scale (λ) as the simulation domain inlet length scale. This validates the suitability of the numerical wind tunnel for our purposes. In contrast, simulations using the integral length scale at the inlet boundary resulted in significant differences in lift coefficients when compared to the experimental results. This demonstrates that accurately capturing the evolution of the turbulent kinetic energy upstream of the airfoil is crucial for reproducing its aerodynamic behaviour in numerical simulations. The comparison between the chord-wise pressure coefficients obtained from experiments and the digital twin revealed that they were similar on the suction side. However, significant differences in C_p values were observed on the pressure side. As the force coefficients obtained using load cells matched well with those obtained in the digital twin, this suggests that there may be room for improvement in the pressure measurement experiments. Moving forward, we plan to perform simulations at higher Reynolds numbers of approximately O(10⁶) to further validate our numerical wind tunnel.

A comparison between 2D and 3D simulations also indicates that, if the emphasis is purely on force coefficients, the 2D digital twin can be used for low to moderately high AoAs due to its computational efficiency. For higher AoAs, where 3D effects are substantial, the more computationally demanding 3D digital twin is preferable. This strategic blend optimizes the use of the digital twin framework.

The authors advise the simulation community to be aware of the significant impact that the Taylor micro-scale (λ) has on the dissipation rate, which in turn affects the evolution of turbulent kinetic energy (TKE). Although different simulation codes may define length scales differently, they should all produce the same dissipation rate. To ensure consistency, the boundary condition for length scales in simulations should either be the Taylor micro-scale itself or a function of it (depending on the definition used in a particular code). In the authors' study, the Taylor micro-scale is used directly as the boundary condition.

The methodology presented in this paper is tailored to regular-grid-generated turbulent flows within a specific scope. As a next step, we intend to investigate anisotropic flows, particularly shear flows, in order to gain insights into the downstream and transverse evolution of turbulent kinetic energy. Our goal is to extend Eq. (23), currently formulated in 1D, into its 2D counterpart. In future work, we plan to conduct simulations at higher Reynolds numbers, around Re_c∼O(10⁶), to validate our numerical wind tunnel in the context of elevated Reynolds numbers.

Appendix A: Details of experiments performed at the University of Oldenburg

The experiments from the University of Oldenburg were carried out in the large wind tunnel that has an inlet of 3 m × 3 m and a test section length of 30 m. A single hot wire was operated using a StreamLine 9091N0102 frame with a 91C10 CTA (Constant Temperature Anemometry) module. It was sampled at f_s=20 kHz; a hardware low-pass filter with a cut-off frequency of f=10 kHz was set. The experiments were performed using a passive regular grid with M=115 mm, and 33 downstream positions spanning from $x/M\sim \mathrm{8}$ to $x/M\sim \mathrm{170}$ were traversed for two inflow speeds: 5 and 10 m s⁻¹. The experimental data were obtained at the span-wise centre at the height of 8M above the floor of the wind tunnel.