The main equation of the DFSEM (Divergence Free SEM) is the following:

Obsiously the (3) could be replaced by several other functions, since it is just a "shape function". is a independent variable with a zero mean distribution and unit variance that basically represents the eddy intensity. We notice the crossproduct in the (1), which implies a relation between the 3 components of the and , very important in our following considerations about the Reynold Stresses.
From our first implementation of this DFSEM, we notice our incapability to keep under control the calculated Reynolds stress component. So I decided to calculate them from the (1), to check from which variable they depends. Because of the cross product in the (1), I was forced to calculate them component by component. Here I post two of them, as an example:

As we can see from (4) and (5), in both the calculated components there is a dependence from the spatial coordinates (mainly the reciprocal position of the grid point considered and every single eddies position). I stress that in our calculation, .
From all these discussions, we stressed the need to implement a formulation of the DFSEM in whitch we can have a stronger control over the Reynolds stresses output.
First of all we focus our attention on the calculated Reynolds stress in our simulation.
We can easily see from the first graph a qualitative behaviour of the calculated Reynold stresses. The components converge towards a value in all the simulations, instead the mixed component converge towards a value lower than 0. Such a behaviour is expected from the theory, where we saw in equatios (4) and (5). These equation in fact have some different feature:
The first one explain us why the calculations converge towards a positive or negative value, the second one explain the different absolute value (in (4) we have the sum of two squared quantities instead in (5) is just a product).
Unluckily, as I said previously, at the moment we do not have the situation in hand. In fact, as we can see in (1), the Reynold stress are not used in the velocity calculation, and this is the main limit of the DFSEM method at the moment.
Question: how the SEM method considered the Reynolds stress imput?
The SEM method used the Lund coefficient (a Cholesky decomposition of the Reynolds tensor) that modified the "random" calculated velocity field with the purpouse to generate the right turbulent stresses (basically these coefficients strech or press the velocity field in order to get what we want). Unluckily such a method cannot be used in a Divergence Free method, otherwise we lose this feature of the field, making all our modified algorithm useless.
Simulation with the following shape function
Steady value: 0.317  

Simulation with the following shape function
Steady value: 0.195 ( This morning I found a bug in this simulation!! A new one is running right now!! ) 
We try to generate a velocity field that satisfies the 7 followings relations:
The main problem is to make cohabit both the divergence free and the Reynolds stress condition. Further to best fit the turbulence characteristics we have to use all these relations: a model that do not consider the RANS Reynolds stresses loses a very important information we can not afford to lose and beyond the divergence free is the main feature of any turbulent fenomena.
First of all let's have a look at what we can modify. Basically, we can modify only the shape function and the eddy length scale, but none of these parameters allow us to have a complete control over the stresses. At the moment seems the DFSEM as in (1) does not allow us any control. We have then to modify such a method.
I thought to use the old SEM and to make it divergencefree using an appropriate eddyintensity generation. This trik does not give us any good results: in fact if we use appropriate eddyintensity generation, these three components will not be anymore independents with a zero mean distribution and unit variance, and this is one hypothesis we need to keep the Reynold stress under control in this SEM.