Wind energy: from downscaled forecasts to mills production - Part 1
CC Attribution - NonCommercial - NoDerivatives 2.0 Generic:
You are free to use, copy, distribute and transmit the work or content in unchanged form for any legal and non-commercial purpose as long as the work is attributed to the author in the manner specified by the author or licensor.
Rousseau, Antoine
In this talk, I will introduce the stochastic downscaling method (SDM) that borrows techniques from small scale turbulence (S.B. Pope) for the simulation of wind flows thanks to hybrid methods (deterministic-stochastic). I will present the downscaling method used to refine a wind forecast at a sufficiently small scale, and the way wind turbines are implemented in the model. Comparisons with traditional numerical methods (LES) and validation w.r.t. experimental data will also be provided.
 Thank you Emmanuel. I would like to thank the organizers for giving me the opportunity to give this lecture. And I would like to thank SMF and SMAE to let us organize this SEMRAX which is a very nice feature of the mathematical community I guess in France. So I will talk about wind energy in my first talk and I will give a second one related to what we have heard this week regarding optimization for bioreactors. I'm talking about the talk of Amira and the one of Simon yesterday. So it will be related to what you heard this week. So both of my talks will be approximately 30 minutes. Of course don't hesitate to ask questions within the talks. It's quite informal. So this first one is a work done mainly with Mireille Bosse. Mireille is a research director at the Sofia anti-police. She was actually my postdoc advisor 10 years ago at Sofia. And then I moved to Grenoble to do some fluid mechanics. But I keep a strong collaboration with Mireille. And we work with Philippe Drobinski who is a geophysicist at the EPSL and Ecole Polytechnique at CNRS. And Christiane Paris who is an engineer at INREA Chile where we work also for this wind energy production theme. And most of the work that will be presented here has been funded by INREA of course. And Adam. So I wanted to let you know. So SDM holds for stochastic downscaling method. So each time you see SDM that's what it means. So the challenges are to provide some small scale simulations for the wind at the scale of a wind farm for example. And unfortunately for the moment the major scale simulations, the meteorological tools that we have today are not able to provide a wind simulation at a finer mesh than maybe one kilometer in the horizontal scale and a few tens of meters in the vertical. So if you want to consider an area of a few kilometers square for a wind farm, you have to downscale your model to have finer informations. So for that of course in fluid mechanics and geophysics there are a lot of tools based on grid refinement for example. So I'm not at all going to talk about that. What I want to introduce today is a quite new method based on stochastic algorithms and particles that will use the mesoscale provided by some deterministic simulations. It can be also provided by data or whatever you want. And the stochastic downscaling method will launch a small scale simulation on the region of interest such as a wind farm for example. So SDM and wind pose which is the software related to this method work on a smaller domain and a finer mesh. With boundary conditions as I told you such as mesoscale data provided by the meteor of France or WARF, the traditional meteorological solver used in the US and in Europe. And it can be also by reconstructed boundary velocities from measurements, data, whatever you want. So the principle is the following is that you throw a large number of particles inside your domain and you run a Lagrange model that I will present afterwards for which we use a fluid particles discretization. And these fluid particles follow some stochastic differential equations that I will introduce. So the main point of that kind of simulation is that thanks to the particular discretization that we use with fluid particles afterwards we can of course average the amount of fluid which the particle properties to obtain some mean velocities at some locations of the fluid. But we can also obtain turbulence and second order and higher momentums of the velocity for example which means that together with the mean information that is required to obtain so for example the wind speed you can also have some error bars and confidence intervals. That's very important because the stochastic nature of the model allows to obtain both the mean quantities and the second moment. So again the computational domain is one or several cells of a large scale simulation. We want to downscale again. Boundary condition I told you about. And then inside our domain we use the Slagengine model which has no stability condition. Actually I will come back to that. And the statistical advantages that I just told you about which is that to obtain the error bars you don't need to run the code many times with a different boundary condition, different initial conditions, etc. and do some sensibility analysis. You can do this directly thanks to the stochastic nature of the model. So SDM is somehow a new method. Actually it has been developed in the 80s by Stephen Pope, a physician in the US. He only used this for reactive fluids at very small scales. And the main point of our work was to use that kind of method but at a much larger scale for geophysical fluids. And of course it raises new problems in modelling numerical discretization, well-posedness of stochastic systems, etc. Okay. So the outline is the following. I will do a brief introduction on turbulent fluids. Then I will talk about the SDM model. And then I will talk about validation because we have new results regarding this. So the Reynolds decomposition, maybe some of you know about that. So you just write your velocity, pressure and whatever parameter function that you have as the sum of its mean and its turbulent part or its deviation to the mean. And this deviation depends on the alia omega. And of course you recover the mean part thanks to an averaging and integral over dp which is the measure. And you plug this into the Navier-Stokes equation or Euler equations as usual. And if you take the brackets of this, of course, you have the linear parts provide the bracket terms like this one. And for the nonlinear parts, you have some kind of problem which means that if you call this bracket of Ui, if you call it Vi, for example, you don't have an equation for Vi. Because in the equation of Vi, you have Vi here, you have Vi here, you have P here, okay. But here you have quantities that depend on Ui, on Uij, okay, here. Thanks to the nonlinear term. So you cannot close the system. That's the traditional problem with this decomposition. So you can imagine that you can try to obtain the equation for the second moment. So you write the equation of the second order, you take the brackets of this. So if you call it Vi, again, you can try to have an equation for Vi. But if you do this because of the nonlinear terms, you obtain a third order terms, et cetera. So each time you want to have an equation for the nth moment of the velocity, you have the momentum number n plus 1 which comes into the equation and so on. So you cannot close the system. So you have, of course, different ways to close that. So the literature is tremendous on that point and, of course, the main thing is to give a parametrization of that kind of term. And if you do it in the first equation, you call it closure of order 1, order 2 if you take two equations, et cetera. So that's the way Rantz, which holds for Reynolds average Navier-Stokes equations, works. So in the Lagrangian approach, we will describe the fluid thanks to particles with a state vector which is position, velocity, and some other physical quantities such as temperature or salinity if you are in the ocean, humidity if you are in the atmosphere, et cetera. But I will focus on X and U. And this is typical of Landgeuvent type models. And the associated system of SdEs, of course, that we are going to write, have to be consistent in a way with Euler or Navier-Stokes. If you want to downscale Navier-Stokes or Euler in a small box, the model that you run in the small box has to be consistent in a way with the large-scale model. So this is how it works. So of course, dxt equals ut dt. That means that the derivative of the position is the velocity. No surprise with this. And this equation is given for the velocity of the particles. I will forget the third one because I won't use it afterwards. So here you will cover the pressure gradients. And some terms which I will talk about later on. And please notice that there is a Brownian motion w at the end of this equation that makes the system stochastic. So how do you compute brackets of U and brackets of Ui, Uj? So again, as I told you, we are going to run numerically some particles. And at the end, what you want to know is what is the wind velocity at this point. For example, if you put a wind farm and a windmill at some precise location, you want to have the wind velocity at this location. For that, we will average the velocities of the particles that are located around this point. And this will provide the wind velocity at this point. So we go from the Lagrangian formulation to an Eulerian quantity, thanks to Monte Carlo averages, actually. And this is what is written here, actually. The Eulerian density of the particles is computed thanks to an average of the Lagrangian one. In other words, and maybe you can just focus on this one. If you want to know the average of U at some point X, you just multiply Vi by the Lagrangian density, and you divide by the average, which exactly is the formula for conditional expectation in probability. So this means this is the expectation of UiT, knowing that XT equals X. That means that I take all the particles that have a position XT equals X, and I compute their average velocity. And this gives me exactly the Eulerian velocity at position X. Again, remember that a particle carries both its location and its velocity. So when you compute a conditional expectation knowing that X of the particle equals to an X that you chose in the domain, then you are considering the computation of Monte Carlo average of particles around the position X. Of course, this is from the continuous point of view, and XT equals X in the numerical context has to be defined with a nearest gate point or some approximation. Okay, I don't want to go too much into details with the Fokker-Planck equation, but we can discuss about that later on. From the density of the equation, from the equations I showed you here, you can compute the Fokker-Planck equation associated to the density fL, which is a PDE. And if you take this PDE and you multiply it by dv d phi and you integrate, then you recover exactly the mass conservation equation. If you do the same thing, but you multiply by v dv d phi, then you recover the first momentum equation of your Euler or Navier-Stokes equation, Euler here. And if you multiply by v i, v j, et cetera, you obtain exactly the other equations that I showed a few slides ago on the second order. So in summary, and here I use quotes of course, but it's just to say that if you take brackets of our Lagrangian model with the particles, then you recover your Euler model that I showed at the beginning. So this is just to say that the Lagrangian model that I showed at the beginning was compatible in a way with the Navier-Stokes or Euler equations that were running at the large scale. And of course, this has been studied really seriously by Mireille and a former PhD student Jean-François Javier with very nice results of the existence and uniqueness of solutions. Okay. So now the specific SDM model, I just wrote more precisely the terms that were in red a few minutes ago. So this is what is running our code. So it's, okay, again, it's an equation very simple on the position. If you derive the position, you get the velocity. And for the velocity, we have this one where k is the turbulent kinetic energy. So it's actually the sum of the second order, second momentum of the velocity. And epsilon is the production of kinetic energy, and it has to be defined, and I will come to it afterwards. And from the numerical point of view, you will see that we use a kind of prediction-correction scheme as a chorion and temam. And so you recover the pressure thanks to Poisson problem, classically, as we do in fluid dynamics. So it means that the numerical scheme is done in the Lagrangian part for the particles. And when you correct the pressure, it's done in the Eulerian system. So you have to go back and forth from the Lagrange to the Euler coordinates every time in the code. But I will come back to it. So to this, you have to add some boundary conditions. And in the SD configuration, that's how it works. So the x is actually just the velocity that is imposed at the boundary of our domain. It's boundary condition. And that means that when the particle reaches the boundary, then these two terms make the particle get the correct velocity imposed by the large-scale simulation. This is written here. This ensures that the conditional expectation of view with respect to x t equals x, where x is on the boundary, is exactly vx. So just a few words on the boundary layer theory. So we are interested in a domain which is approximately 1,000 meters high. It corresponds to the size of the atmospheric boundary layer. And inside this, you have a sublayer, which is called the surface layer, where important things occur. And inside this, you have the roughness layer, which is, again, a smaller layer. But actually, this one is not taken into account in the... Well, it's taken into account, but it's not simulated by the code. It's just parameterized thanks to a roughness parameter. I will come back to it. So in our code, we simulate this. But with the parametrization of the blue part, which is very difficult to simulate at most scales, of course. So on the ground, the goal is to account for the log low in the surface layer, in the green one. The log low is, I guess, lots of you know about this. But the thing is that you want to reach a vertical profile in z of this kind with u star, which is given by the second momentum, and k is a constant. And z0, again, is the roughness length, which is a parameter which provides information regarding the configuration in the sublayer. So for example, on the ocean, z0 is almost zero. But if you have buildings or things like that, z0 is a few tens of meters or something like that. So it's a typical length that is used in the model. But it's a parameter. And in the whole domain, so on the whole atmospheric boundary layer, we use this formula for the pseudo dissipation. Remember that epsilon was coming into the equations a few minutes ago. And the LM is the mixing length, so it explains how the turbulence is created on the ground and in the layer. Again I don't want to talk too much about this because it could take a few hours, actually. But I have lots of documents regarding this. So if you have questions and queries about this, I can provide lots of information. So again, k, you can compute it thanks to your particles with this formula. It's again a conditional expectation. And epsilon is parametrized thanks to this formula. And then we can go to the numerics. So again, this is our domain, and as I told you, we drop a lot of particles inside the domain. We have an external force, an external velocity which is provided by the large-scale simulation. And if you want to compute f of u, whatever f you have, you do just a Monte Carlo average of the particle field, okay? So you compute f of u, k, k is an index for the particle. For every particle, you compute f of u, k, and you just make sure that this particle k is located in the close to the position x, okay? So it means that if x here is away from the position of x, this capital X, k, then this is equal to zero, okay? And on the contrary, if it's located in the control cell around x, then this is one, and this is what is written here, okay? So it's just an indicator of whether your particle is close to the position you're looking at or not. And if yes, okay, this is one, this is the sum of the j corresponding to the particles at this position, and that's it. So it's just an average, f of u is just an average of the f of u, k, where k is the index of other particle fields. So this is one time step of the algorithm. So you start by computing a new position, sorry, thanks to the former velocity, and you compute thanks to a Scheme for stochastic equations, these terms, you forget the pressure term for the moment. Then you make a correction of the particles because we rely on the hypothesis that the density is constant here, so which is a very strong hypothesis. We agree that it's not correct at this scale. I mean, in a few centimeters or within a meter of vertical size, it's reasonable, but over tens of meters it's not, so that's something that has to be improved. But for the moment, we rely on this approximation. So if we want to have a constant density, that means that from the particular point of view, we have to have a uniform repetition of the particles in the domain. And since we've just moved the particles, thanks to this equation, we have to do some optimal transport problem, actually, to remove the cloud of particles from the position it has to a new one, which is a uniform distribution. And why is this an optimal transport problem? Is that you don't want to move the particle too far from the place there because what we are doing here is not physical, okay? So you don't want to take some physical information from one part of the domain to the other part. You want to stay close to the place where you were. So that's why it's an optimal transport problem. And this step has been done in the Sembrax 2007. The code has been written and the method has been validated, et cetera. And the last step is to correct the particle velocity this time, not the particle position, but velocity, thanks to the classical Poisson equation for the pressure. And here, what I want to stress on is that it's an Eulerian equation, okay? So after step one and step two, you compute from the Lagrangian information on the particles, you compute the Eulerian field. And this is only an equation on the Eulerian field, okay? So it's pretty fast, actually, because you can have many, many, many particles. But here in the box, this kind of thing can be done very easily thanks to fast Poisson solvers such as fish pack or that kind of thing. And of course, classically, you correct the particle velocity from the particle field. Okay. So validation now. We just, we compare the simulations with the LES method that is used by, classically, let's say by geophysicists and in particular, in particular, Philippe Drobinski. So the computational domain is three kilometers, one kilometer large in horizontal dimensions and almost 800 meters in the vertical. So it's approximately the size of the atmospheric boundary layer. Okay, and what you can see here is three couples of curves. So these are the variances of U, V and W. And each time you have computed with LES or with ARM model. And okay, they look very similar even if in the close to the ground, okay, sorry, this is the altitude and this is the variances in Y. So close to the ground, close to X equals zero, which is Z equals zero actually, we have a few discrepancies but afterwards it's pretty good. And here it's the same with the, with the covariances. Yes, covariances. So U, W, V, W, et cetera. Okay. So it can compare pretty well to LES and here are the main momentum. So that's the U velocity, U in the west-east direction and this is north-east. And okay, you can see that in the 100 first meters, you have this log low that I talked to you about. And again, this is a comparison between LES and the SDM, which is at several times and you know, even after a long time of simulation, it's really good. Again we have some discrepancies on the floor for the 4V. Okay and now towards wind farm simulation. So this was validated with no wind farm inside the domain. And Christian Paris put some actor disc model into the numerical code to simulate the presence of windmills. Okay, so I can show you what typical production of a kinetic, turbulent kinetic energy beyond the mill. So you have a constant velocity in U coming from the left and the mill turns and this is the production of a turbulent kinetic energy behind the mill. Okay. This is done with Numesis, the numerical platform of VINREA on anti-police. And finally, we just validated this model thanks to realistic wind data provided by Fernando Portear-Hel in a PFL. So he made up a small wind tunnel, one meter high, in which he implemented some windmills and many measures. He has many data on this. And okay, this is the configuration. You have a mill at maybe one meter from the entrance. You have a constant wind coming from the left and you want to see the profile of the velocity everywhere. So this, as you can see, the wall moves, it's just to see in 3D the velocity field from the numerical point of view. Okay, so that's just the configuration. And what you can see here is the velocity computed at several spots. So this is X. You have your mill here. The wind is coming this way. And the first simulation is computed here. So this is the vertical profile of the u-velocity just behind the mill. So as you can see, there is this kind of belly here, which is due to the fact that the mill pumps some energy from the wind to produce some electricity. Okay, of course. So this is why the u-velocity diminishes just after the mill. So in blue, you have the SDM simulation. In green, you can barely see it, but it's the measurements from Porte-Aguel. And in red, it's the Jensen wake effect model, which is classically used by wind solvers for EDF and the GDF suede, et cetera, in their... So what is interesting here is that SDM does better. And the more you go on the right, the further you go from the mill. And at the end, you almost don't see anymore the impact of the mill on the vertical profile. And if I drew a last plot, you would have that kind of thing, which would be the profile without any impact of the mill. So that was the validation we did with Fernando Porte-Aguel. And sorry. And I think that's it for this presentation. So again, a few conclusions and perspectives. Lagrangian formulation. So this is new for downscaling techniques at that scale. Theoretical problems in probability theory rates. So a few achievements have been done by Mireille Bossin and Jean-François Javier, as I told you. So the link with meteorology has been done thanks to a calibration and comparison, both with experimental data and LES methods that are classically used in that field. And as I told you at the very beginning, our numerical simulations also provide some variability on the fields, which is interesting. Of course, we have a lot of improvements to do. Non-flat domains, it's almost done. Actually, we are dealing with non-flat domains, which is not very easy, but it's almost done. But the non-neutral case means row non-constant. So it would imply probably an equation on the temperature, because you cannot take rid of the temperature evolution at that scale. In the wind tunnel, the reason why our simulations are quite good in comparison with the data is that it's one meter high, the temperature is constant, and everything is in a lab, clean, et cetera. But we don't want to test SDM now on a wind farm simulation. We are sure that the prediction would not be correct, because we are missing some physics still. Okay. And of course, we need to do some optimization and parallelization, and I wanted to end with this because one of the key arguments of using a stochastic and particular method is that it can be parallelized very easily, well, more easily than with the classical deterministic tools, because there is no interaction between the particles themselves. They actually interact through the Eulerian field. Okay. In the equation, you have an equation on the particle, and you have a link with the Eulerian field. So of course, they interact with each other from the continuous point of view, but in the numerical simulation, they interact through the Eulerian field. So if at every time step you compute the Eulerian field, then you can completely parallelize the evolution on the particles, which we think is interesting, but still we have to do it. Okay. And you can browse the SDM.inria.fr. We have some more explanations and simulations. I'm done for this first talk. Thank you for your attention. Okay. So for this first talk, we have to do some questions. Maybe we can take time for a few questions. Please. Thank you. So you said that you will talk about the fact that your method do not have to respect the CFL condition, but in fact you didn't say no. Yes, yes, you're right. Thank you. Thank you for that. Yes. Actually, there is a hidden CFL condition. Well, on stochastic particular algorithms, there is no. But here, when we do the optimal transport thing, for example, we don't want to move the particles too far from their previous position for physical reasons. And that is actually a hidden CFL because it means that you don't want the particle to cross too many cells in one time step. So it's exactly a CFL condition, actually. So it's not, well, it's not, there is no theory about this, but in numerical cases, we could absolutely observe that we could not take any time step. There was a condition on the time step, of course. And linked to the, of course, to the side of the mesh. Okay. Is there another question? I have a question. You said that your method would not work on a real wind farm. So for the time being, what is used to build or to make choice concerning wind farms? Yes. For the moment, it's, you have softwares like, again, what is used in the industry. It's mostly interpolation on the large scale. So of course, they use some meteorological data provided by WAF or whatever weather forecasting system. But at the end, when they want some information locally, it mostly works with interpolations. So there is no physics. What I mean is there is no physics in the downscaling systems that are used today in that kind of... But how do they take into account the presence of mill and so... Presence of? Of mill. For that, they use the same kind of thing with the actuator disk. It's quite, well, it's not difficult. I mean, it's a force, mechanic force that you add at the right-hand side of the equations. And it's not a big deal. I mean, the mill. But the most important point is to have the velocity at some specific locations, getting rid of the mills for the moment. Even getting rid of the mills, the way they downscale the information is done thanks to interpolations, which is completely false. But don't you think that your method, even at its level of development, wouldn't give better result than pure interpolations? I don't know. I don't know because the fact that the temperature is not taken into account is really a big deal, I think. Because the wind, the thermal effects near the ground are really, really important. So I wouldn't, I don't want to rank too bad methods. But these effects are taken into account in the interpolation method? Yes, yes, yes. Okay. Yes, because they interpolate the velocity and the temperature and everything provided by the large scale. So in a way, they are taken into account. Yes. Okay. So another comment or question. So thank you again. Thank you. So. Thank you.
 Thank you Emmanuel.  I would like to thank the organizers for giving me the opportunity to give this lecture.  And I would like to thank SMF and SMAE to let us organize this SEMRAX which is a very  nice feature of the mathematical community I guess in France.  So I will talk about wind energy in my first talk and I will give a second one related  to what we have heard this week regarding optimization for bioreactors.  I'm talking about the talk of Amira and the one of Simon yesterday.  So it will be related to what you heard this week.  So both of my talks will be approximately 30 minutes.  Of course don't hesitate to ask questions within the talks.  It's quite informal.  So this first one is a work done mainly with Mireille Bosse.  Mireille is a research director at the Sofia anti-police.  She was actually my postdoc advisor 10 years ago at Sofia.  And then I moved to Grenoble to do some fluid mechanics.  But I keep a strong collaboration with Mireille.  And we work with Philippe Drobinski who is a geophysicist at the EPSL and Ecole Polytechnique  at CNRS.  And Christiane Paris who is an engineer at INREA Chile where we work also for this wind  energy production theme.  And most of the work that will be presented here has been funded by INREA of course.  And Adam.  So I wanted to let you know.  So SDM holds for stochastic downscaling method.  So each time you see SDM that's what it means.  So the challenges are to provide some small scale simulations for the wind at the scale  of a wind farm for example.  And unfortunately for the moment the major scale simulations, the meteorological tools  that we have today are not able to provide a wind simulation at a finer mesh than maybe  one kilometer in the horizontal scale and a few tens of meters in the vertical.  So if you want to consider an area of a few kilometers square for a wind farm, you have  to downscale your model to have finer informations.  So for that of course in fluid mechanics and geophysics there are a lot of tools based  on grid refinement for example.  So I'm not at all going to talk about that.  What I want to introduce today is a quite new method based on stochastic algorithms  and particles that will use the mesoscale provided by some deterministic simulations.  It can be also provided by data or whatever you want.  And the stochastic downscaling method will launch a small scale simulation on the region  of interest such as a wind farm for example.  So SDM and wind pose which is the software related to this method work on a smaller domain  and a finer mesh.  With boundary conditions as I told you such as mesoscale data provided by the meteor of  France or WARF, the traditional meteorological solver used in the US and in Europe.  And it can be also by reconstructed boundary velocities from measurements, data, whatever  you want.  So the principle is the following is that you throw a large number of particles inside  your domain and you run a Lagrange model that I will present afterwards for which we use  a fluid particles discretization.  And these fluid particles follow some stochastic differential equations that I will introduce.  So the main point of that kind of simulation is that thanks to the particular discretization  that we use with fluid particles afterwards we can of course average the amount of fluid  which the particle properties to obtain some mean velocities at some locations of the fluid.  But we can also obtain turbulence and second order and higher momentums of the velocity  for example which means that together with the mean information that is required to obtain  so for example the wind speed you can also have some error bars and confidence intervals.  That's very important because the stochastic nature of the model allows to obtain both  the mean quantities and the second moment.  So again the computational domain is one or several cells of a large scale simulation.  We want to downscale again.  Boundary condition I told you about.  And then inside our domain we use the Slagengine model which has no stability condition.  Actually I will come back to that.  And the statistical advantages that I just told you about which is that to obtain the  error bars you don't need to run the code many times with a different boundary condition,  different initial conditions, etc. and do some sensibility analysis.  You can do this directly thanks to the stochastic nature of the model.  So SDM is somehow a new method.  Actually it has been developed in the 80s by Stephen Pope, a physician in the US.  He only used this for reactive fluids at very small scales.  And the main point of our work was to use that kind of method but at a much larger scale  for geophysical fluids.  And of course it raises new problems in modelling numerical discretization, well-posedness of  stochastic systems, etc.  Okay.  So the outline is the following.  I will do a brief introduction on turbulent fluids.  Then I will talk about the SDM model.  And then I will talk about validation because we have new results regarding this.  So the Reynolds decomposition, maybe some of you know about that.  So you just write your velocity, pressure and whatever parameter function that you have  as the sum of its mean and its turbulent part or its deviation to the mean.  And this deviation depends on the alia omega.  And of course you recover the mean part thanks to an averaging and integral over dp which  is the measure.  And you plug this into the Navier-Stokes equation or Euler equations as usual.  And if you take the brackets of this, of course, you have the linear parts provide the bracket  terms like this one.  And for the nonlinear parts, you have some kind of problem which means that if you call  this bracket of Ui, if you call it Vi, for example, you don't have an equation for Vi.  Because in the equation of Vi, you have Vi here, you have Vi here, you have P here, okay.  But here you have quantities that depend on Ui, on Uij, okay, here.  Thanks to the nonlinear term.  So you cannot close the system.  That's the traditional problem with this decomposition.  So you can imagine that you can try to obtain the equation for the second moment.  So you write the equation of the second order, you take the brackets of this.  So if you call it Vi, again, you can try to have an equation for Vi.  But if you do this because of the nonlinear terms, you obtain a third order terms, et  cetera.  So each time you want to have an equation for the nth moment of the velocity, you have  the momentum number n plus 1 which comes into the equation and so on.  So you cannot close the system.  So you have, of course, different ways to close that.  So the literature is tremendous on that point and, of course, the main thing is to give  a parametrization of that kind of term.  And if you do it in the first equation, you call it closure of order 1, order 2 if you  take two equations, et cetera.  So that's the way Rantz, which holds for Reynolds average Navier-Stokes equations,  works.  So in the Lagrangian approach, we will describe the fluid thanks to particles with a state  vector which is position, velocity, and some other physical quantities such as temperature  or salinity if you are in the ocean, humidity if you are in the atmosphere, et cetera.  But I will focus on X and U. And this is typical of Landgeuvent type models.  And the associated system of SdEs, of course, that we are going to write, have to be consistent  in a way with Euler or Navier-Stokes.  If you want to downscale Navier-Stokes or Euler in a small box, the model that you  run in the small box has to be consistent in a way with the large-scale model.  So this is how it works.  So of course, dxt equals ut dt.  That means that the derivative of the position is the velocity.  No surprise with this.  And this equation is given for the velocity of the particles.  I will forget the third one because I won't use it afterwards.  So here you will cover the pressure gradients.  And some terms which I will talk about later on.  And please notice that there is a Brownian motion w at the end of this equation that  makes the system stochastic.  So how do you compute brackets of U and brackets of Ui, Uj?  So again, as I told you, we are going to run numerically some particles.  And at the end, what you want to know is what is the wind velocity at this point.  For example, if you put a wind farm and a windmill at some precise location, you want  to have the wind velocity at this location.  For that, we will average the velocities of the particles that are located around this  point.  And this will provide the wind velocity at this point.  So we go from the Lagrangian formulation to an Eulerian quantity, thanks to Monte Carlo  averages, actually.  And this is what is written here, actually.  The Eulerian density of the particles is computed thanks to an average of the Lagrangian one.  In other words, and maybe you can just focus on this one.  If you want to know the average of U at some point X, you just multiply Vi by the Lagrangian  density, and you divide by the average, which exactly is the formula for conditional expectation  in probability.  So this means this is the expectation of UiT, knowing that XT equals X.  That means that I take all the particles that have a position XT equals X, and I compute  their average velocity.  And this gives me exactly the Eulerian velocity at position X.  Again, remember that a particle carries both its location and its velocity.  So when you compute a conditional expectation knowing that X of the particle equals to an  X that you chose in the domain, then you are considering the computation of Monte Carlo  average of particles around the position X.  Of course, this is from the continuous point of view, and XT equals X in the numerical context  has to be defined with a nearest gate point or some approximation.  Okay, I don't want to go too much into details with the Fokker-Planck equation, but we can  discuss about that later on.  From the density of the equation, from the equations I showed you here, you can compute  the Fokker-Planck equation associated to the density fL, which is a PDE.  And if you take this PDE and you multiply it by dv d phi and you integrate, then you  recover exactly the mass conservation equation.  If you do the same thing, but you multiply by v dv d phi, then you recover the first  momentum equation of your Euler or Navier-Stokes equation, Euler here.  And if you multiply by v i, v j, et cetera, you obtain exactly the other equations that  I showed a few slides ago on the second order.  So in summary, and here I use quotes of course, but it's just to say that if you take brackets  of our Lagrangian model with the particles, then you recover your Euler model that I showed  at the beginning.  So this is just to say that the Lagrangian model that I showed at the beginning was  compatible in a way with the Navier-Stokes or Euler equations that were running at the  large scale.  And of course, this has been studied really seriously by Mireille and a former PhD student  Jean-François Javier with very nice results of the existence and uniqueness of solutions.  Okay.  So now the specific SDM model, I just wrote more precisely the terms that were in red  a few minutes ago.  So this is what is running our code.  So it's, okay, again, it's an equation very simple on the position.  If you derive the position, you get the velocity.  And for the velocity, we have this one where k is the turbulent kinetic energy.  So it's actually the sum of the second order, second momentum of the velocity.  And epsilon is the production of kinetic energy, and it has to be defined, and I will come to  it afterwards.  And from the numerical point of view, you will see that we use a kind of prediction-correction  scheme as a chorion and temam.  And so you recover the pressure thanks to Poisson problem, classically, as we do in fluid dynamics.  So it means that the numerical scheme is done in the Lagrangian part for the particles.  And when you correct the pressure, it's done in the Eulerian system.  So you have to go back and forth from the Lagrange to the Euler coordinates every time  in the code.  But I will come back to it.  So to this, you have to add some boundary conditions.  And in the SD configuration, that's how it works.  So the x is actually just the velocity that is imposed at the boundary of our domain.  It's boundary condition.  And that means that when the particle reaches the boundary, then these two terms make the  particle get the correct velocity imposed by the large-scale simulation.  This is written here.  This ensures that the conditional expectation of view with respect to x t equals x, where  x is on the boundary, is exactly vx.  So just a few words on the boundary layer theory.  So we are interested in a domain which is approximately 1,000 meters high.  It corresponds to the size of the atmospheric boundary layer.  And inside this, you have a sublayer, which is called the surface layer, where important  things occur.  And inside this, you have the roughness layer, which is, again, a smaller layer.  But actually, this one is not taken into account in the...  Well, it's taken into account, but it's not simulated by the code.  It's just parameterized thanks to a roughness parameter.  I will come back to it.  So in our code, we simulate this.  But with the parametrization of the blue part, which is very difficult to simulate at most  scales, of course.  So on the ground, the goal is to account for the log low in the surface layer, in the green  one.  The log low is, I guess, lots of you know about this.  But the thing is that you want to reach a vertical profile in z of this kind with u star,  which is given by the second momentum, and k is a constant.  And z0, again, is the roughness length, which is a parameter which provides information  regarding the configuration in the sublayer.  So for example, on the ocean, z0 is almost zero.  But if you have buildings or things like that, z0 is a few tens of meters or something  like that.  So it's a typical length that is used in the model.  But it's a parameter.  And in the whole domain, so on the whole atmospheric boundary layer, we use this formula for the  pseudo dissipation.  Remember that epsilon was coming into the equations a few minutes ago.  And the LM is the mixing length, so it explains how the turbulence is created on the ground  and in the layer.  Again I don't want to talk too much about this because it could take a few hours, actually.  But I have lots of documents regarding this.  So if you have questions and queries about this, I can provide lots of information.  So again, k, you can compute it thanks to your particles with this formula.  It's again a conditional expectation.  And epsilon is parametrized thanks to this formula.  And then we can go to the numerics.  So again, this is our domain, and as I told you, we drop a lot of particles inside the  domain.  We have an external force, an external velocity which is provided by the large-scale simulation.  And if you want to compute f of u, whatever f you have, you do just a Monte Carlo average  of the particle field, okay?  So you compute f of u, k, k is an index for the particle.  For every particle, you compute f of u, k, and you just make sure that this particle  k is located in the close to the position x, okay?  So it means that if x here is away from the position of x, this capital X, k, then this  is equal to zero, okay?  And on the contrary, if it's located in the control cell around x, then this is one, and  this is what is written here, okay?  So it's just an indicator of whether your particle is close to the position you're looking  at or not.  And if yes, okay, this is one, this is the sum of the j corresponding to the particles  at this position, and that's it.  So it's just an average, f of u is just an average of the f of u, k, where k is the index  of other particle fields.  So this is one time step of the algorithm.  So you start by computing a new position, sorry, thanks to the former velocity, and you compute  thanks to a Scheme for stochastic equations, these terms, you forget the pressure term  for the moment.  Then you make a correction of the particles because we rely on the hypothesis that the  density is constant here, so which is a very strong hypothesis.  We agree that it's not correct at this scale.  I mean, in a few centimeters or within a meter of vertical size, it's reasonable, but over  tens of meters it's not, so that's something that has to be improved.  But for the moment, we rely on this approximation.  So if we want to have a constant density, that means that from the particular point of  view, we have to have a uniform repetition of the particles in the domain.  And since we've just moved the particles, thanks to this equation, we have to do some  optimal transport problem, actually, to remove the cloud of particles from the position it  has to a new one, which is a uniform distribution.  And why is this an optimal transport problem?  Is that you don't want to move the particle too far from the place there because what  we are doing here is not physical, okay?  So you don't want to take some physical information from one part of the domain to the other part.  You want to stay close to the place where you were.  So that's why it's an optimal transport problem.  And this step has been done in the Sembrax 2007.  The code has been written and the method has been validated, et cetera.  And the last step is to correct the particle velocity this time, not the particle position,  but velocity, thanks to the classical Poisson equation for the pressure.  And here, what I want to stress on is that it's an Eulerian equation, okay?  So after step one and step two, you compute from the Lagrangian information on the particles,  you compute the Eulerian field.  And this is only an equation on the Eulerian field, okay?  So it's pretty fast, actually, because you can have many, many, many particles.  But here in the box, this kind of thing can be done very easily thanks to fast Poisson  solvers such as fish pack or that kind of thing.  And of course, classically, you correct the particle velocity from the particle field.  Okay.  So validation now.  We just, we compare the simulations with the LES method that is used by, classically,  let's say by geophysicists and in particular, in particular, Philippe Drobinski.  So the computational domain is three kilometers, one kilometer large in horizontal dimensions  and almost 800 meters in the vertical.  So it's approximately the size of the atmospheric boundary layer.  Okay, and what you can see here is three couples of curves.  So these are the variances of U, V and W. And each time you have computed with LES or  with ARM model.  And okay, they look very similar even if in the close to the ground, okay, sorry, this  is the altitude and this is the variances in Y.  So close to the ground, close to X equals zero, which is Z equals zero actually, we  have a few discrepancies but afterwards it's pretty good.  And here it's the same with the, with the covariances.  Yes, covariances.  So U, W, V, W, et cetera.  Okay.  So it can compare pretty well to LES and here are the main momentum.  So that's the U velocity, U in the west-east direction and this is north-east.  And okay, you can see that in the 100 first meters, you have this log low that I talked  to you about.  And again, this is a comparison between LES and the SDM, which is at several times and  you know, even after a long time of simulation, it's really good.  Again we have some discrepancies on the floor for the 4V.  Okay and now towards wind farm simulation.  So this was validated with no wind farm inside the domain.  And Christian Paris put some actor disc model into the numerical code to simulate the presence  of windmills.  Okay, so I can show you what typical production of a kinetic, turbulent kinetic energy beyond  the mill.  So you have a constant velocity in U coming from the left and the mill turns and this  is the production of a turbulent kinetic energy behind the mill.  Okay.  This is done with Numesis, the numerical platform of VINREA on anti-police.  And finally, we just validated this model thanks to realistic wind data provided by Fernando  Portear-Hel in a PFL.  So he made up a small wind tunnel, one meter high, in which he implemented some windmills  and many measures.  He has many data on this.  And okay, this is the configuration.  You have a mill at maybe one meter from the entrance.  You have a constant wind coming from the left and you want to see the profile of the velocity  everywhere.  So this, as you can see, the wall moves, it's just to see in 3D the velocity field from  the numerical point of view.  Okay, so that's just the configuration.  And what you can see here is the velocity computed at several spots.  So this is X.  You have your mill here.  The wind is coming this way.  And the first simulation is computed here.  So this is the vertical profile of the u-velocity just behind the mill.  So as you can see, there is this kind of belly here, which is due to the fact that the mill  pumps some energy from the wind to produce some electricity.  Okay, of course.  So this is why the u-velocity diminishes just after the mill.  So in blue, you have the SDM simulation.  In green, you can barely see it, but it's the measurements from Porte-Aguel.  And in red, it's the Jensen wake effect model, which is classically used by wind solvers  for EDF and the GDF suede, et cetera, in their...  So what is interesting here is that SDM does better.  And the more you go on the right, the further you go from the mill.  And at the end, you almost don't see anymore the impact of the mill on the vertical profile.  And if I drew a last plot, you would have that kind of thing, which would be the profile  without any impact of the mill.  So that was the validation we did with Fernando Porte-Aguel.  And sorry.  And I think that's it for this presentation.  So again, a few conclusions and perspectives.  Lagrangian formulation.  So this is new for downscaling techniques at that scale.  Theoretical problems in probability theory rates.  So a few achievements have been done by Mireille Bossin and Jean-François Javier, as I told  you.  So the link with meteorology has been done thanks to a calibration and comparison, both  with experimental data and LES methods that are classically used in that field.  And as I told you at the very beginning, our numerical simulations also provide some variability  on the fields, which is interesting.  Of course, we have a lot of improvements to do.  Non-flat domains, it's almost done.  Actually, we are dealing with non-flat domains, which is not very easy, but it's almost done.  But the non-neutral case means row non-constant.  So it would imply probably an equation on the temperature, because you cannot take  rid of the temperature evolution at that scale.  In the wind tunnel, the reason why our simulations are quite good in comparison with the data  is that it's one meter high, the temperature is constant, and everything is in a lab, clean,  et cetera.  But we don't want to test SDM now on a wind farm simulation.  We are sure that the prediction would not be correct, because we are missing some physics  still.  Okay.  And of course, we need to do some optimization and parallelization, and I wanted to end with  this because one of the key arguments of using a stochastic and particular method is  that it can be parallelized very easily, well, more easily than with the classical deterministic  tools, because there is no interaction between the particles themselves.  They actually interact through the Eulerian field.  Okay.  In the equation, you have an equation on the particle, and you have a link with the Eulerian  field.  So of course, they interact with each other from the continuous point of view, but in  the numerical simulation, they interact through the Eulerian field.  So if at every time step you compute the Eulerian field, then you can completely parallelize  the evolution on the particles, which we think is interesting, but still we have to do it.  Okay.  And you can browse the SDM.inria.fr.  We have some more explanations and simulations.  I'm done for this first talk.  Thank you for your attention.  Okay.  So for this first talk, we have to do some  questions.  Maybe we can take time for a few questions.  Please.  Thank you.  So you said that you will talk about the fact that your method do not have to respect  the CFL condition, but in fact you didn't say no.  Yes, yes, you're right.  Thank you.  Thank you for that.  Yes.  Actually, there is a hidden CFL condition.  Well, on stochastic particular algorithms, there is no.  But here, when we do the optimal transport thing, for example, we don't want to move  the particles too far from their previous position for physical reasons.  And that is actually a hidden CFL because it means that you don't want the particle  to cross too many cells in one time step.  So it's exactly a CFL condition, actually.  So it's not, well, it's not, there is no theory about this, but in numerical cases,  we could absolutely observe that we could not take any time step.  There was a condition on the time step, of course.  And linked to the, of course, to the side of the mesh.  Okay.  Is there another question?  I have a question.  You said that your method would not work on a real wind farm.  So for the time being, what is used to build or to make choice concerning wind farms?  Yes.  For the moment, it's, you have softwares like, again, what is used in the industry.  It's mostly interpolation on the large scale.  So of course, they use some meteorological data provided by WAF or whatever weather forecasting  system.  But at the end, when they want some information locally, it mostly works with interpolations.  So there is no physics.  What I mean is there is no physics in the downscaling systems that are used today in  that kind of...  But how do they take into account the presence of mill and so...  Presence of?  Of mill.  For that, they use the same kind of thing with the actuator disk.  It's quite, well, it's not difficult.  I mean, it's a force, mechanic force that you add at the right-hand side of the equations.  And it's not a big deal.  I mean, the mill.  But the most important point is to have the velocity at some specific locations, getting  rid of the mills for the moment.  Even getting rid of the mills, the way they downscale the information is done thanks to  interpolations, which is completely false.  But don't you think that your method, even at its level of development, wouldn't give  better result than pure interpolations?  I don't know.  I don't know because the fact that the temperature is not taken into account is really a big  deal, I think.  Because the wind, the thermal effects near the ground are really, really important.  So I wouldn't, I don't want to rank too bad methods.  But these effects are taken into account in the interpolation method?  Yes, yes, yes.  Okay.  Yes, because they interpolate the velocity and the temperature and everything provided  by the large scale.  So in a way, they are taken into account.  Yes.  Okay.  So another comment or question.  So thank you again.  Thank you.  So.  Thank you.
2015-08-13_Rousseau_part1_0_47.png 2015-08-13_Rousseau_part1_48_95.png 2015-08-13_Rousseau_part1_96_155.png 2015-08-13_Rousseau_part1_156_203.png 2015-08-13_Rousseau_part1_204_275.png 2015-08-13_Rousseau_part1_276_395.png 2015-08-13_Rousseau_part1_396_479.png 2015-08-13_Rousseau_part1_480_1847.png 2015-08-13_Rousseau_part1_1848_1895.png 2015-08-13_Rousseau_part1_1896_1943.png 2015-08-13_Rousseau_part1_1944_3407.png 2015-08-13_Rousseau_part1_3408_6635.png 2015-08-13_Rousseau_part1_6636_9119.png 2015-08-13_Rousseau_part1_9120_11183.png 2015-08-13_Rousseau_part1_11184_11339.png 2015-08-13_Rousseau_part1_11340_11387.png 2015-08-13_Rousseau_part1_11388_11447.png 2015-08-13_Rousseau_part1_11448_11543.png 2015-08-13_Rousseau_part1_11544_11591.png 2015-08-13_Rousseau_part1_11592_12899.png 2015-08-13_Rousseau_part1_12900_14039.png 2015-08-13_Rousseau_part1_14040_14123.png 2015-08-13_Rousseau_part1_14124_15983.png 2015-08-13_Rousseau_part1_15984_17327.png 2015-08-13_Rousseau_part1_17328_18443.png 2015-08-13_Rousseau_part1_18444_22319.png 2015-08-13_Rousseau_part1_22320_22487.png 2015-08-13_Rousseau_part1_22488_22547.png 2015-08-13_Rousseau_part1_22548_23039.png 2015-08-13_Rousseau_part1_23040_24779.png 2015-08-13_Rousseau_part1_24780_24839.png 2015-08-13_Rousseau_part1_24840_28247.png 2015-08-13_Rousseau_part1_28248_29699.png 2015-08-13_Rousseau_part1_29700_32699.png 2015-08-13_Rousseau_part1_32700_33131.png 2015-08-13_Rousseau_part1_33132_39959.png 2015-08-13_Rousseau_part1_39960_42107.png 2015-08-13_Rousseau_part1_42108_43319.png 2015-08-13_Rousseau_part1_43320_43391.png 2015-08-13_Rousseau_part1_43392_43439.png 2015-08-13_Rousseau_part1_43440_43487.png 2015-08-13_Rousseau_part1_43488_43571.png 2015-08-13_Rousseau_part1_43572_43619.png 2015-08-13_Rousseau_part1_43620_43967.png 2015-08-13_Rousseau_part1_43968_44687.png 2015-08-13_Rousseau_part1_44688_45479.png 2015-08-13_Rousseau_part1_45480_46187.png 2015-08-13_Rousseau_part1_46188_48839.png 2015-08-13_Rousseau_part1_48840_48899.png 2015-08-13_Rousseau_part1_48900_49007.png 2015-08-13_Rousseau_part1_49008_60179.png 2015-08-13_Rousseau_part1_60180_60299.png 2015-08-13_Rousseau_part1_60300_60359.png 2015-08-13_Rousseau_part1_60360_60407.png 2015-08-13_Rousseau_part1_60408_60490.png