Dynamic Wake Meandering Model

Fig. 7 Exemplary application of the DWM model to the simulation of the Lillgrund Wind Farm in QBlade.

The Dynamic Wake Meandering (DWM) model is based on the premise that the physical mechanisms governing wind turbine wakes can be mathematically decoupled. The model splits the wake dynamics into three distinct components:

1. Wake Deficit Evolution: The quasi-steady axisymmetric expansion and recovery of the velocity deficit in a rotor-bound frame of reference.

2. Wake Meandering: The stochastic, transverse displacement of the wake planes driven by large-scale atmospheric turbulence.

3. Added Turbulence: The micro-scale, high-frequency turbulence mechanically generated by the rotor shear layer.

Initial Deficit and Boundary Condition

The initial velocity deficit immediately behind the rotor is derived from the blade axial induction factor, \(a(r)\). Because the flow expands as it passes through the rotor disc, the control volume radius must be expanded to conserve mass.

QBlade implements several boundary condition models to determine the initial expanded wake radius, \(r_{out}\), and the initial normalized wake velocity, \(U(r)\):

None (Simple Momentum): Assumes no stream-tube expansion.

\[U(r) = 1 - 2a(r)\]

Madsen: Applies an expansion factor \(f_w = 1 - 0.45 a_{ave}^2\) based on the area-averaged induction \(a_{ave}\). The expanded radial coordinates are computed iteratively:

\[r_{out, i+1} = \sqrt{\frac{1-a_i}{1-2a_i}(r_{i+1}^2 - r_i^2) + r_{out, i}^2} \cdot f_w\]

IEC: Calculates the total thrust coefficient \(C_T = 4(1-a_{ave})a_{ave}\) and an expansion parameter \(m = 1/\sqrt{1-C_T}\) to determine the final expanded radius:

\[r_{out} = 2r (1 - 0.45 a_{ave}^2) \sqrt{\frac{1+m}{8}}\]

Keck: Uses an empirical relationship to account for turbulence build-up:

\[r_{out} = r \sqrt{\frac{1-a_{ave}}{1-1.98 a_{ave}}}\]

\[U(r) = 1 - 2.1 a(r)\]

To ensure numerical stability across all models, the local induction factor is strictly capped at \(|a| \le 0.47\).

Wake Deficit Evolvution

Once the initial boundary condition is set, the quasi-steady wake deficit recovers as it propagates downstream. This is solved using the Thin Shear Layer approximation of the Reynolds-Averaged Navier-Stokes (RANS) equations in an axisymmetric coordinate system.

Momentum Conservation

The axial momentum equation describes how the deficit is smoothed out by turbulent mixing (eddy viscosity, \(\nu\)):

\[U \frac{\partial U}{\partial x} + V \frac{\partial U}{\partial r} = \frac{1}{r} \frac{\partial}{\partial r} \left( \nu r \frac{\partial U}{\partial r} \right)\]

This partial differential equation is solved numerically in QBlade using an implicit finite-difference scheme. A transition matrix is constructed mapping the velocity profile at step \(x\) to \(x + dx\) and solved via a Sparse LU decomposition (Eigen::SparseLU).

Continuity (Radial Velocity)

The radial velocity, \(V(r)\), required for the momentum equation, is derived from the conservation of mass:

\[\frac{\partial U}{\partial x} + \frac{1}{r} \frac{\partial (rV)}{\partial r} = 0\]

Integrating this yields the radial velocity at any point:

\[V(r) = -\frac{1}{r} \int_0^r r' \frac{\partial U}{\partial x} dr'\]

Eddy Viscosity Models

The momentum recovery is entirely dependent on the turbulent eddy viscosity, \(\nu\), which is modeled as a linear combination of the ambient atmospheric turbulence and the wake-added turbulence generated by the shear layer:

\[\nu = \nu_{ambient} + \nu_{wake}\]

QBlade includes several closure models for the eddy viscosity (e.g., Madsen, Larsen, IEC, Keck). These models generally scale the viscosity using the ambient turbulence intensity (\(TI\)), the downstream distance (\(x\)), the instantaneous wake width, and the maximum velocity deficit (\(1 - \min(U)\)).

For example, the shear-driven wake viscosity is typically formulated as:

\[\nu_{wake} \propto \text{width} \cdot (1 - \min(U(r)))\]

Wake Meandering and Advection

While the thin shear layer equations govern the shape of the deficit, atmospheric turbulence dictates where that deficit is physically located.

Spatial Averaging

Large atmospheric eddies (larger than the rotor diameter) push the entire wake structure bodily. To isolate these large-scale movements from the background micro-turbulence, the ambient wind field is spatially averaged over a moving polar grid at the wake plane’s location.

The solver computes two distinct spatial averages:

1. Advection Velocity: Averaged over a polar grid of diameter \(C_{advect} \cdot D\). This determines the downstream propagation speed (\(dx/dt\)).

2. Meandering Velocity: Averaged over a polar grid of diameter \(C_{meander} \cdot D\). This determines the transverse displacement in the \(y\) and \(z\) directions.

To filter the ambient wind field smoothly, weighting functions are applied to the polar grid points. QBlade utilizes uniform weighting or Jinc filters, where the JINC2 spatial filter is defined as:

\[w(x) = \frac{J_1(x)}{x} \frac{J_1(x/2)}{x/2}\]

Low-Pass Filtering and Deflection

The properties of the emitted wake planes (thrust, yaw, tilt, and ambient free-stream velocity) are smoothed using an exponential moving average low-pass filter to prevent unphysical high-frequency jitter:

\[x_{lp,t} = x_{lp,t-1} \cdot e^{-2\pi f_c \Delta t} + \left(1-e^{-2\pi f_c \Delta t}\right) \cdot x_t\]

If the rotor is yawed or tilted, a lateral or vertical deflection velocity is superimposed onto the wake plane’s kinematics:

\[\Delta y_{yaw} = f_{yaw} \cdot \theta_{yaw} \cdot \Delta x\]

\[\Delta z_{tilt} = f_{tilt} \cdot \theta_{tilt} \cdot \Delta x\]

Where \(f_{yaw}\) and \(f_{tilt}\) are tuning factors scaling the deflection relative to the downstream travel \(\Delta x\).

Wake Added Turbulence

Because the spatial averaging process removes high-frequency fluctuations, the micro-scale turbulence generated by the turbine’s shear layer must be artificially reintroduced.

A pre-generated, isotropic, unit-variance turbulence wind field (\(\vec{V}_{field}\)) is superimposed onto the flow field during velocity evaluations. The magnitude of this added turbulence is scaled dynamically based on the local velocity deficit and the radial velocity gradient (shear):

\[\vec{V}_{turb} (\vec{x},t) = k_m(x,r) \cdot \vec{V}_{field}(\vec{x},t)\]

The scaling factor, \(k_m\), evaluates the local state of the wake using two tunable constants (\(k_{m1}\) and \(k_{m2}\)):

\[k_m(x,r) = k_{m1} |1 - U(x,r)| + k_{m2} \left|\frac{\partial U(x,r)}{\partial r}\right|\]

This ensures that added turbulence peaks in the highly-sheared annular region of the wake (the blade tips/roots) and smoothly dissipates as the velocity deficit recovers downstream.

Fig. 8 3D Visualization of the wake velocity deficit evaluated with the DWM model.