b) Simulation using Virtual Disks


The context

In this part, the idea is to simulate an existing wind farm and evaluate its performances. The farm under study is the following :

 

 

 

 

 

 

 

 

 

 

 

 

It is composed by 5 wind turbines Repower MM70 aligned and separated from a distance of $2.85 D_{blades}$.

Here, the wind is represented by the vectors of the top right of the figure.

The geometrical characteristics of the MM70 are gathered in the figure :

 

In StarCCM+: first, the domain is created, then a first simulation is launch. To finish the results are exploited.


How does this simulation work ?

 

The wind turbines are represented by actuator disks (in blue here):

 

The "virtual disk" is a model in StarCCM+ which simulate the presence of a wind turbine by modifying the Navier-Stokes equations adding a forcing term: in the cells composing the virtual disk, two forces are added to the equation solved by the software : a normal force, and a tangential one.

StarCCM+ takes as inlet :

 

and geometrical consideration of the turbine : the position, the orientation...

 

With that, the software adds two forces to the momentum equation for every cell whose center is inside the virtual disk :

  • a Thrust force $$F^{\perp}_{cell}=T \frac{V_{cell}}{\sum V_{cell}}$$

where T is the real thrust force : $$T=\frac{1}{2} \rho_0 U_0^2 C_t \pi (R_0^2-R_1^2)$$

V the volume of the cell, $U_0$ the velocity impacting the turbine, $C_t$ the thrust coefficient (given by the curve), $\rho_0$ the density of the air and $R_0$ $R_1$ the geometry of the turbine.

  • a Torque $$F^{//}_{cell}=Q \frac{r_{cell}^2 V_{cell}}{\sum r_{cell}^2 V_{cell}}$$

where Q is the real torque of a wind turbine : $$Q=\frac{P}{\Omega}$$

V the volume of the cell, r its distance from the center P the power produced (given by the power curve) and $\Omega$ the rotation rate.

 

Domain - Region - Physics Continuum

 


The Domain

 

First, it is important to create the domain. Its characteristics are gathered here :

In order to do that, first create a new simulation. Then the idea is to create the block:

geometry  $\rightarrow$ parts. Right click New Shape Part $\rightarrow$ Block

 

The domain can be visualized by right clicking on Scene $\rightarrow$ Geometry.

Right click Block $\rightarrow$ rename "Domain"

It is important to split the surface. Indeed, at this point, the surface gather all the boundaries of the domain which would make it impossible to define different limit conditions.

To do that, Parts $\rightarrow$ domain $\rightarrow$ Surfaces. Right click on split by angle. Ok.

Rename every boundaries by right clicking on the surface $\rightarrow$ rename.


The region

Then we have to create the region :

Parts right click on Domain, $\rightarrow$ Assign part to regions.

(be careful to select Create a Boundary for Each Part Surface).

 


Physics continuum

It is now important to create a physics continuum : right  click on Continua $\rightarrow$ New $\rightarrow$ Physics Continuum.

Select the physic required by the simulation. Here :

It is now necessary to enter the characteristics of the wind turbine : StarCCM+ requires the power and thrust curves. So, a file *.csv is created with the wind velocity in the first column, the power in the second and the thrust coefficient in the third.

In StarCCM+ : Tools $\rightarrow$ Table. Right click New $\rightarrow$ File Table.  Choose the *.csv you previously created.

 

Now, one can create the virtual disk in the physics continuum :

right click $\rightarrow$ New right click on virtual disk $\rightarrow$ edit. Method $\rightarrow$ 1D momentum.

The configuration chosen for the first disk is the following :

It is now easier to copy and paste this virtual disk and translate the copy of $2.85 D_{blades}$ by right clicking on Virtual Disk 1 $\rightarrow$ copy $\rightarrow$ paste.

 


Boundary Conditions and Initialization

Now, let us define the boundary conditions (right click on the Boundary $\rightarrow$ Edit) :

One can define a field function as velocity inlet and initialization to represent accurately the atmospheric boundary layer :

Tools $\rightarrow$ right click on field function $\rightarrow$ New $\rightarrow$ Vector. Here, a class 2 wind is used :

Now, this function "in_velocity" can be used as initialized velocity and inlet velocity :

Regions $\rightarrow$ Region $\rightarrow$ Boundaries $\rightarrow$ Domain.inlet $\rightarrow$ Physics values $\rightarrow$ Velocity magnitude (right click $\rightarrow$ edit $\rightarrow$ field function ) and then choose "in_velocity".

The turbulence intensity needs to be 12% on the inlet as well. Physics values $\rightarrow$ turbulence intensity $\rightarrow$ 0.12.

 

The next step is to create the mesh.

Mesh

 


The mesh

Now is the time to create the mesh :

In order to do so, first, create a mesh continuum : right click continua $\rightarrow$ New $\rightarrow$ Mesh continuum. To choose models : right click models $\rightarrow$ Select models. Choose :

Here are the values chosen for the poly mesh :

The parameters are quite easily understandable. It is of utmost importance to refine the mesh in the wake of wind turbines. To achieve that, creating Volumes mesh is necessary :

Tools $\rightarrow$ Volume shapes $\rightarrow$ right click New $\rightarrow$ Cylinder.

 

It is facultative but cones immediately after the virtual disk can be created as well (the same method can be followed). To redefine cell size in those volume shapes, go in Continua $\rightarrow$ Mesh 1 $\rightarrow$ Volumetric controls right click New. In Volume Shape, define Shapes as the Cylinder previously created. Use the parameters :

 

The mesh can now be generated by clicking successively on :

 

 


Creating a visualization of the mesh with virtual disks

To get a 3D view of the mesh (like the first image of this page), a "Threshold" derived part is used. It makes it possible for the user to display cell which meet a specific criteria.

To display the virtual disk, the cells whose centers are inside the disk are displayed. Right click on Derived parts $\rightarrow$ New Part $\rightarrow$ Treshold. Use the following :

 

Note that the disks will only appear after an initialization.

It is possible to do the same to have a cut through the domain. This time, the scalar field will be "Position[X]". It will display the cells whose center are inside the chosen range of X.

Finally, one obtains:

Simulations - Results

 


Simulations

From this point, it is possible to run a simulation. First, initialize then run . After convergence (it might be necessary to use the first order : Continua $\rightarrow$ Physics 1 $\rightarrow$ right click on Models $\rightarrow$ Edit $\rightarrow$ Segregated Solver $\rightarrow$ first order) one wants to actually see the wake :

 


Visualization

To do so, the velocity is seen on a plane cutting the domain.

Right click on Derived part $\rightarrow$ Section $\rightarrow$ Plane.

Then, a scalar scene is created. Right click on Scene $\rightarrow$ New Scene $\rightarrow$ Scalar scene.

Now that the scalar scene is opened, click on Scene/Plot .

Displayers $\rightarrow$ Scalar 1 $\rightarrow$ Parts. Choose the threshold part created for every virtual disks and the Plane section. In Scalar Field, choose velocity magnitude.

Finally :

 

The visualization obtained is :

 


Quantify the results

The next step of the study is to compare the simulation results with experimental data.

Those experimental data are normalized power produced by each turbine(normalized by the power produced by the first turbine) in two cases :

  • every turbine working
  • the 4th turbine stopped

To quantify the simulation results it is necessary to get the velocity entering each wind turbines. To do so, the velocity is averaged on a disk just before the virtual disk :

The method will be quickly explained here :

first create a part :geometry $\rightarrow$ right click on parts $\rightarrow$ New Shape part $\rightarrow$ cylinder. Place it conveniently. As it was done for the domain, split the surface in 3 (right click on surface $\rightarrow$ split by angle). To be faster, copy and paste the cylinder then transform $\rightarrow$ translate to place it just in front the second disk.

However, now the created cylinders are independant from the region. So if the averaged velocity is taken at this point, it will be 0. That is why a derived part is used.

Right click on derived part $\rightarrow$ New part $\rightarrow$ Section $\rightarrow$ Arbitrary section. Use the following parameters :

Now, let us create a report : right click on report $\rightarrow$ New report $\rightarrow$ surface average :

To display the velocities : right click on report $\rightarrow$ Run all reports.

With these velocities, one can interpolate on the MM70 power curve (with the function interp1 of Matlab for example) and obtain the power produced by each turbine. After normalization :

The order of magnitude remains the same, however, the drop in term of power after the first turbine is bigger in the experimental data. Moreover, the simulations do not predict the stabilization of the loss after the second turbine.


Conclusion and possible improvements

In conclusion, this part of the project made it possible to simulate an existing wind farm and compare its production to experimental data. To improve the results, one can do again the simulations with some changes :

  • creating a bigger domain to make sure that the limit conditions do not have an impact on the solution
  • try to use the first order and then use the second order to improve precision.
  • improve the way of getting the velocities with StarCCM+ to be closer to reality.
  • use the V2F turbulence model which theorically predict better a mixing flow.