I - Wake of Wind Turbines

In this part of the project, two kinds of simulation will be successively conducted in StarCCM+ :

 

 

       

a) Simulation using an Overset Mesh

In this part, the aim is to run a simulation which makes it possible to see what a wind turbine wake actually looks like.

 


The principle

 

The principle of an overset mesh is really simple : the mesh is divided in two parts :

  • a background one for the entire domain. It needs to be refined enough where the turbine is located.
  • a foreground one for the domain close to the turbine. It will be rotating. Obviously, the size of the cells needs to be of the same order as the refined background one.

These two parts share have a small zone in common where the two mesh are superposing :

source : doc StarCCM+

 

More details about this method can be found in the greatly explained StarCCM+ documentation.

Geometry - Region - Physics Continuum


Geometry

 

As in the simulation with virtual disks, the wind turbine studied is the Repower MM70 :

 

As it is difficult to create a turbine geometry from scratch in StarCCM+, a possible solution is to import an *.stl file created with a CAD software. It is easy to find one on the internet. The one used here is taken from the website GrabCAD.

 

Firstly, the *.stl geometry file is imported : Geometry $\rightarrow$ right click on 3D-CAD models $\rightarrow$ New $\rightarrow$ right click on 3D-CAD model 1 $\rightarrow$ import $\rightarrow$ CAD models.

Next, the idea is to modify the bodies to be as close as possible to the geometry of the MM70.

Here, mainly Scaling is used. Finally :

After this step being done, parts need to be created : here, to simplify, only blades and the nacelle will be represented.

Geometry $\rightarrow$  3D-CAD models $\rightarrow$ 3D-CAD model 1 $\rightarrow$ bodies $\rightarrow$ right click on the body $\rightarrow$ New geometry part.

Now, let us create the domain : right click on Part $\rightarrow$ New Shape Part $\rightarrow$ Block. Be careful to create a domain big enough to avoid problems of limit conditions having an impact on the solution. Split the domain's surface Domaine $\rightarrow$ Surfaces $\rightarrow$ right click on Surface $\rightarrow$ Split by angles. Rename the surfaces created (inlet, outlet, side, floor...).

The next step is to create the moving part : it is defined as

To do that, create a cylinder which include the geometry (Part $\rightarrow$ New Shape Part $\rightarrow$ Cylinder). Then, select both part (the geometry part and the cylinder) $\rightarrow$ Boolean $\rightarrow$ Substract parts. Do the same as before to split the surfaces.


The region

Select both the domain and Substract $\rightarrow$ Assign Parts to region :

Rename the regions : domain and moving_region.

Now, create an interface : select both region $\rightarrow$ right click $\rightarrow$ Create interface $\rightarrow$ Overset Mesh. Check that the interface node created has the moving region as region-0 and the domain as region-1.

Set the boundary conditions : for the moving_region, everything has a type "wall" excepted the cylinder surface whose type is "Overset mesh". For the domain :


The Physics Continuum

To create a physics continuum right click on Continua $\rightarrow$ New $\rightarrow$ Physics Continuum. Now select the models : right click on models $\rightarrow$ Select Models. Select :

It is now generally necessary to add a source term : "This source term is required for such cases because the overset mesh method is not conservative and comes with an interpolation error that prevents residuals from converging. STAR-CCM+ automatically calculates the source term." (starCCM+ documentation).
To activate it : right click on Models $\rightarrow$ Edit $\rightarrow$ Segregated Solver

It is now possible to put initial conditions for the velocity (10 m/s here) and turbulence intensity (0.12) :

It is the good moment to set up the boundary conditions (0.12 as turbulence intensity and 10m/s as entering velocity):

 

 

The next step is the creation of the mesh

Mesh

 


The Mesh

The principle is quite simple : the two regions defined before are meshed separately : the software will automatically deactivate the cell of the background mesh which are not useful. Indeed, the Navier Stokes equations won't be solved in the background cells which are close to the blades geometry (in the area close to the blade geometry, the equations will only be solved in the moving mesh).

First, create two mesh continuum : do twice the procedure right click on Continua $\rightarrow$ New $\rightarrow$ Mesh Continuum. Rename one as "moving_mesh" and the other as "background mesh".

On the region domain : region $\rightarrow$ right click on domain $\rightarrow$ Edit $\rightarrow$ Mesh continuum $\rightarrow$ "background mesh". Same procedure for the "moving_region" domain $\rightarrow$ Mesh continuum $\rightarrow$ "moving_mesh".

For both meshes, it is now necessary to select the models to use : righ click on models $\rightarrow$ Polyhedral Mesher and Surface Remesher for the background mesh and Polyhedral Mesher, Surface Remesher, Prism Layer mesher and Surface wrapper for the moving_mesh.

The values used here for the moving mesh are :

 

For the background mesh :

It is important that the cell close to the geometry in the background mesh have the same order of magnitude in term of size. To achieve that, let us create a Volume shape : Tools $\rightarrow$ right click on Volume shapes $\rightarrow$ New Shape $\rightarrow$ Block. This block must be at lease the size of the moving region.

Now it is possible to add a volumetric control :

right click on Volumetric Controls $\rightarrow$ New $\rightarrow$ Select the shape just created and use Polyhedral Mesher and Surface Remesher with a relative size of 2 (which gives 1m, the base size of the moving mesh).

From this point, one can create the mesh by clicking successively on .


Make the Mesh Moving

To make the mesh moving, Tools $\rightarrow$ right click on Motions $\rightarrow$ New $\rightarrow$ Rotation. Click on rotation $\rightarrow$

Now that the rotation is created, it is necessary to let the software now that we want to assign this Motion the moving_region :

here,

Right click $\rightarrow$ Motion $\rightarrow$ Rotation


Visualize the mesh

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.

Right click on Derived part $\rightarrow$ New part $\rightarrow$ Threshold $\rightarrow$

Right click on Scene $\rightarrow$ New Scene $\rightarrow$ Mesh in the onglet Plot/scene ():

 

Finally :

Finally, it is possible to set up the simulation.

Simulation - Vizualisation

 


Simulation parameters

In the node Implicit Unsteady

set up the Time-step. Here, 0.05s is used.

In the node Stopping Criteria the user has to define the conditions which will make the iterations stop : the logical rules are

  • OR : if the condition is satisfied, it does the action.
  • AND : if the condition is satified, it still iterates until all AND conditions are satisfied or a OR condition is satisfied.

Here, two condition are enabled :

In conclusion, the software will do 20 inner iterations for every 0.05s time-step until it reachs 100s.

 

Now, the idea is to save simulation results, let's say every 0.1s to make an animation for example. Right click on Solution Histories $\rightarrow$ New $\rightarrow$ "load". In fact, the results wanted will be solved in the "load.simh" file.

Right click on load $\rightarrow$

That means : when the simulation reachs 70s (the transcient part is not wanted here) the scalar and vector functions chosen (Pressure, TVR and velocity) in the moving_region and domaine are going to be saved in "load.simh" every Delta time =0.1s until the end of the simulation (100s).

Now that everything is set up, one can initialize and launch the simulation . It could be a good idea to launch starccm+ in background with a batch java function and the "nohup" shell command.

 


Create a Visualization

Now that the simulation is done, let us visualize the solution :

To do so, create a scalar scene : right click on Scene $\rightarrow$ New $\rightarrow$ Scalar. 

With this scalar node "pressure", the pressure on the wind turbine is seen. Now, to see the wake on this same scene, the Turbulent viscosity ratio is used. Two choice are possible :

  • add a TVR isosurface
  • add a resampled volume (which is the choice done here)

The second solution will be explained here :

right click on Derived Part $\rightarrow$ Resampled Volume. The idea is to place this volume where the wake should actually be : after the wind turbine.

Now in the same scene as before, in the scene/plot onglet () Right click on Displayer $\rightarrow$ New Displayer $\rightarrow$ Scalar. Here, as a part select the resampled Volume and as a scalar the TVR. In the scalar field properties, disable the Auto Range.

Add a solution time annotation to the scene by draging Solution Time to the scene : Tools $\rightarrow$ Annotation $\rightarrow$ Solution Time. (This way, it is possible to add logos to the scene as well. Create a new 2D image annotation).

Now, it is easily noticeable that only the solution at the last time is seen. To have access to the solutions stored in "load.simh" : Solution Histories $\rightarrow$ right click on load $\rightarrow$ Create recorded Solution View. In the node Solution View, drag "load" to the scene.Here

it is possible to choose the State visualized.

StarCCM+ makes it possible to create a video *.avi of the solution :

Solution Views $\rightarrow$ Load $\rightarrow$ Animation $\rightarrow$ Solution time $\rightarrow$ Solution Time animation. In the Solution Time Animation Properties, choose the length of the video, the starting time and the ending time :

In order to record the animation, click on record (the right button here ).

Then be careful to define the same animation length as before :

 

Eventually :

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.