# Optimization of a wind farm with StarCCM+

PICOUET Erwan, MOUSSAOUI Hamza

Encadré par : STOUKOV A. (enseignant chercheur IMFT)

LEGENDRE D. (enseignant chercheur IMFT)

IMBERT J. (Sereema consulting)

PINTO B. (Sereema consulting)

The project is divided in two parts with two different objectives :

• simulate an existing wind farm and evaluate its global performances. The results are compared to experimental data. (Erwan is in charge of this part)

• determine the ideal position of a sensor on the nacelle to get an accurate value for the incoming velocity. (Hamza will be in charge)

StarCCM+ is the CFD software used in both cases.

# I - Wake of Wind Turbines

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

• one using the Virtual Disk Method. It will make it possible to simulate an existing wind farm.

• one using an overgrid mesh. It will make it possible to see what a wake looks like.

# 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+

# 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

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 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.

# II - Position of the velocity sensor

Collection of accurate wind speed data is one of the most problematic elements in conducting wind turbine power performance tests. It allows us to send a signal to the control systems in order to stop the turbine when the wind velocity is too hight, in order to prevent material damage.

The owners of wind farms and turbine manufacturers have shown interest in the use of anemometers placed on the platform for the wind speed data collection, rather than using a meteorological tower upstream of the turbine because of the low cost of this solution.

However, the most significant problem with this practice is that the wind flow is disturbed by the rotor and the nacelle and the wind speed measurements gathered by an anemometer located at the rear of the nacelle does not exactly represent the speed of undisturbed wind upstream of the rotor. This problem can be avoided if the measures can be adjusted.

This second part of the project concerns the modeling of flow at the surface of a wind turbine nacelle in order to find the optimal position for an anemometer.

An optimal position means a spot where the flow is the most uniform, and the less turbulent possible. That would make a correlation between the inlet velocity $V_{\infty}$ and the velocity measured by the anemometer $V_{nacelle}$ with the less residual possible.

Key Parameters:

First we have to sort out the parameters influencing the measurement, which are:

• the nacelle geometry,
• the boundary layer around the nacelle,
• the near wake generated by the rotor and the nacelle,
• the atmospheric boundary layer,
• the pitch angle and the position of the anemometer on the nacelle.

The effect of the nacelle geometry will be studied as a first step (a 2D simulation without the nacelle), then a 3D simulation with the presence of the rotor will be carried out:

# a) Simulation of the nacelle (2D)

In this part, we will simulate the wake of the nacelle

We already know that the sensor will be placed on the upper part of the nacelle, so we will focus our study on this upper domaine.

# Sketch of the Domain - Region - Physics Continuum

Sketch

First we create a sketch of the fluid domain.

Right click 3D-CAD Models $\rightarrow$ New

Then we creat a new sketch on the XY plane,

Then we use the entities shown above to creat our sketch.

Then we create a domain with the following dimensions :

$L_H \simeq$ 8 x hight of the nacelle

$L_{upstream} \simeq$  lenght of nacelle

$L_{downstream} \simeq$ 2 x lenght of the nacelle

Then we make an extrude of this sketch, in order to be able to make the simulation :

The domain looks like this:

Then we rename each boundary by right clicking on the surface $\rightarrow$ rename.

The region

Then we create the region:

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

Physics continuum

It is now important to define each kind of Boundary Condition:

Then we have to create a physics continuum by right clicking on Continua $\rightarrow$ New $\rightarrow$ Physics Continuum and select the physic models required by the simulation, here we shose:

# Mesh

Mesh

Now we have to create the mesh,

right click continua $\rightarrow$ New $\rightarrow$ Mesh Continuum, to choose models right click on Models $\rightarrow$ Select Models

The mesh consists of predominantly non-uniform rectangular cells with fine resolution near the nacellein order to capture the boundary around the nacelle. This will allow us to optimize the computing time.

a zoom on the nacelle region shows the refined mesh in this part:

We can also make the mesh finer on a selected block, as the following figure shows:

Finally the mesh is generated by successively clicking on the the following buttons

# Simulations - Results

Results

we follows similar steps to the first part in order to get the results,

We can create a line probe to get the velocity values on this line,

Derived $\rightarrow$ Parts $\rightarrow$ New Part $\rightarrow$ Probe $\rightarrow$ Line...

We plot the streamlines afterwards,

Derived $\rightarrow$ Parts $\rightarrow$ New Part $\rightarrow$ Streamline...

Velocity field is represented in the following figure for an upstream velocity of $10m/s$

A zoom on the nacelle:

We plot here the velocity measured by the sensor $s_1$ as a function of the velocity upstream the rotor:

• Position of sensor $s_1$: $0,5\ m$ above the pod and $0,5\ m$ from the end,
• calibration equation: $y=-0,176+1,044 x$
• the residual norme for this method: $1,6\ .10^{-4}$.

One can notice that the residual is very tiny, because we didn't simulate the rotor which is the most important source of turbulence.
Once we simulate the rotor, one can plot the calibration equation for different points above the nacelle, and find the point with the less residual norme.

# b) Simulation of the rotor (3D)

In this part, we will simulate the wake of the rotor.

# Domain - Region

Domain

First we create the nacelle 3D and a box containing this nacelle, then we substract the nacelle from the box:

Then we define the Bondary Conditions as done in the previous part,

# Mesh - Results

Mesh

We follow exactly the same procedure as in the previous part in order to creat the mesh, and here we create the virtual disk which represent the rotor,

it is explained in the first part how to create a virtual disk,

Results of this 3D simulation couldn't be exploited because of lack of time.

# Conclusion

Conclusion

From this study, we propose the ideal position to place our anemometer:

• in the vertical symmetry plane of the nacelle (by means of geometric symmetry),
• vertically, far from the nacelle's boundary layer,
• not too far from the nacelle (because the wake created by the profiled part of the blades is more turbulent than the wake created by the cylindrical part of it),
• horizontally, far from the blades that create an important wake.