System and method for measuring wind flow turbulence in complex terrain with lidar

Abstract

A wind flow sensing system for determining turbulence of wind flow at different elevation sets above terrain is provided. The wind flow sensing system includes an input interface configured to receive, for a set of time steps, a set of measurements of radial velocity at a site line point above terrain for each elevation, and a processor configured to estimate a velocity field for each elevation based on data assimilation of the velocity field, estimate an unbiased horizontal velocity at each elevation for each time step, determine a mean of the unbiased horizontal velocity for a time period including the set of time steps, and determine turbulence based on the unbiased horizontal velocity for each time step and the mean of the unbiased horizontal velocity.

Description

TECHNICAL FIELD

[0001] The present invention relates generally to remote sensing, and more particularly, to a wind flow sensing system and method for measuring wind flow turbulence in complex terrain using LiDAR. BACKGROUND

[0002] Atmospheric turbulence is a measure of small-scale fluctuations in wind speed that affects many fields, including air quality, aviation, and numerical weather prediction. In particular, turbulence is an important parameter in the wind energy industry, where high-resolution measurements are often required in remote locations. Typically, turbulence is estimated from anemometer measurements on a meteorological tower. However, anemometer measurements are limited by tower height and can suffer from overspeed and slow response time issues, which can result in inaccurate average wind speed and turbulence measurements.

[0003] To address these issues with anemometer measurements, remote sensing devices such as SODAR (sound detection and ranging) and LiDAR (light detection and ranging) that can make remote measurements are used for average wind speed and turbulence measurements. While anemometers measure wind speed at a small point in space, remote sensing devices report average wind speed from a probe volume and measure at a lower frequency than instruments installed on towers, such as anemometers. These differences in spatial and temporal resolution result in differences in turbulence measurements between anemometers and remote sensing devices. Turbulence motions can range from milliseconds to hours and from centimeters to kilometers, while LiDAR can only measure turbulence motions with a time scale of a few seconds and a spatial scale of a few tens of meters. In addition to the differences in spatial and temporal sampling, the scanning strategy used by remote sensing devices can also introduce errors in different turbulence components.

[0004] LiDAR uses Doppler Beam Swinging (DBS) or Velocity-Azimuth Display (VAD) techniques to capture wind speed measurements. Using LiDAR DBS and VAD scans, the variances of the u, v, and w velocity components are not directly measured, rather, the DBS and VAD techniques combine radial velocity measurements from different points to calculate instantaneous values of the velocity components. The time series of u, v, and w are used to calculate the velocity variances, thereby implicitly assuming that the instantaneous velocity values are constant over the scan circle. In turbulence, this assumption is not true even if the mean flow is uniform over the scan circle. Thus, the DBS and VAD methods of calculating variances have the flaw of variance contamination errors.

[0005] Therefore, there remains a need for a system and method suitable for turbulence measurements. SUMMARY

[0006] One purpose of some embodiments is to provide a wind flow sensing system and a wind flow sensing method for determining the turbulence of a wind flow at a set of different altitudes above a terrain from a set of measurements of radial velocities at each altitude of the set of different altitudes. In one embodiment, the determination of the turbulence includes determining a turbulence intensity. In another embodiment, the determination of the turbulence includes determining a turbulence kinetic energy. A purpose of some embodiments also includes estimating a velocity field at each altitude based on data assimilation. Furthermore, one purpose of embodiments is to estimate an unbiased horizontal velocity at each time step at each altitude.

[0007] In some embodiments, remote sensing instruments, such as LiDARs, are used to measure characteristics of the wind in the atmosphere. Characteristics of the wind include wind speed (horizontal velocity and vertical velocity), turbulence, wind direction, etc. LiDARs measure radial velocities of a line-of-sight (LOS) point of the wind at each altitude for a set of time steps. However, the amount of turbulence related to the horizontal velocity is a parameter of interest.

[0008] To this end, some embodiments aim at determining the amount of turbulence related to the horizontal velocity of the wind flow at each altitude. Some embodiments are based on the realization that, with a geometric relationship and assuming that the wind velocity is uniform on each plane inside the cone of measurements, a horizontal velocity estimate at a time step can be determined from the measurements of radial velocities corresponding to the time step. This horizontal velocity estimate is performed for different time instants. Furthermore, the horizontal velocity is averaged over a certain time period (e.g., 10 minutes). The horizontal velocity and the average of the horizontal velocity are different at each time instant. The square of the difference between the horizontal velocity and the average is called the standard deviation.

[0009] According to some embodiments, the standard deviation of the average of the horizontal velocity defines a turbulence intensity (TI). The turbulence intensity is defined over a certain time period (e.g., 10 minutes). To this end, some embodiments are based on the realization that the turbulence intensity is a function of the shape of the signal and the corresponding values, and, therefore, the TI depends on the instantaneous value of the horizontal velocity. However, for wind flows over complex terrain, such as hills or near large buildings or other urban structures, the uniform velocity assumption considered for the horizontal velocity estimate is not valid. Some embodiments are based on the realization that, for complex terrain, the uniform velocity assumption leads to a bias in the horizontal velocity estimate. This bias in the horizontal velocity estimate leads to a biased horizontal velocity, which introduces a bias in the standard deviation. This bias is due to the variations of the vertical velocity in the vertical direction. According to one embodiment, the bias is given by the height times the horizontal gradient of the vertical velocity, i.e.,

[0010] Some embodiments are based on the realization that the horizontal derivative of the vertical velocity can be used as a correction to the biased horizontal velocity to remove the bias. The horizontal derivative of the vertical velocity is determined by estimating a velocity field. The velocity field is estimated at each altitude based on a data assimilation that is used to fit the measured radial velocity. The velocity field is estimated for a set of time steps. In one embodiment, the data assimilation involves determining boundary conditions for the inlet velocity field and performing a computational fluid dynamics (CFD) simulation by solving the Navier-Stokes equations that define the wind flow using the boundary conditions. Further, the boundary conditions are updated and the simulation is repeated until a termination condition is satisfied.

[0011] However, in the CFD simulation, the boundary conditions are unknown and these boundary conditions are iteratively determined until the operating parameters produce the measured radial velocity. Thus, the data assimilation with CFD is very time consuming. Further, the data assimilation with CFD is cumbersome because the CFD simulation is an optimization process based on the solution of the Navier-Stokes equations. Further, the CFD simulation becomes complex for the wind flow over complex terrain.

[0012] To this end, in another embodiment, the data assimilation involves approximating the non-convex shape of the terrain using a set of convex shapes and determining the boundary conditions for the inlet velocity field. Subsequently, an analytical solution to the Laplace equation that defines the wind flow is derived using the boundary conditions. The boundary conditions are updated and the simulation is repeated until a termination condition is satisfied. Since the data assimilation involves an algebraic solution to the Laplace equation rather than an iterative optimization of the Navier-Stokes equation, the efficiency of the operation is improved. Further, the Laplace equation is easier to solve as compared to the Navier-Stokes equation. However, this approximation to the terrain degrades the simulation of the velocity field. Some embodiments are based on the realization that while such an approximation can not be sufficiently accurate for the determination of the velocity field, such an approximation can be sufficiently accurate for the determination of the horizontal derivative of the vertical velocity that is used as a correction to remove the bias in the estimated unbiased horizontal velocity.

[0013] Such bias removal is performed on the horizontal velocity at each time instance to obtain an unbiased horizontal velocity at the respective time instance. In other words, instantaneous bias removal is performed to obtain an unbiased horizontal velocity at each time instance. Further, an average of the unbiased horizontal velocities for the time period is determined. The standard deviation is determined as the square of the difference between the unbiased horizontal velocities and the average of the unbiased horizontal velocities. Since the standard deviation is determined based on the unbiased horizontal velocities, the bias present in the standard deviation due to the biased horizontal velocities is removed. This then improves the accuracy of the standard deviation. According to some embodiments, at each altitude, the turbulence is determined based on the unbiased horizontal velocities at each time step and the average of the unbiased horizontal velocities. In one embodiment, the turbulence includes a turbulence intensity determined based on a ratio of a root mean square value of the turbulence velocity fluctuations to the average of the unbiased horizontal velocities. In another embodiment, the turbulence includes a turbulent kinetic energy determined based on half the sum of the root mean squares of the turbulence velocity fluctuations. Since the turbulence is determined using the unbiased horizontal velocities and the average of the unbiased horizontal velocities, the accuracy of the determined turbulence is significantly improved.

[0014] Some embodiments are based on the realization that a relationship between the standard deviations of different points on the scan circle at a certain altitude can be established by formulating a function, i.e., an autocorrelation function. The autocorrelation function is also referred to as a correction function. Thus, the autocorrelation function is related to the standard deviations of the points and can be used to measure the standard deviation of the estimated horizontal velocities to produce turbulence. Some embodiments are based on the realization that the correction function can be used to correct the unbiased horizontal velocities before estimating the turbulence. The correction function is trained to reduce the difference between ground truth and the determined unbiased horizontal velocities. The ground truth can correspond to the measurements from the cups.

[0015] Accordingly, one embodiment discloses a wind flow sensing system for determining a turbulence of a wind flow at a set of different altitudes above a terrain from a set of measurements of radial velocities at each altitude of the set of different altitudes, the wind flow sensing system comprising: an input interface configured to receive, for a set of time steps, a set of measurements of radial velocities at a ground track point above the terrain for each altitude; a processor configured to: estimate a velocity field for each altitude based on a data assimilation of the velocity field above the terrain, the data assimilation for fitting the measurements of radial velocities, wherein the velocity field is estimated for the set of time steps; estimate an unbiased horizontal velocity at each altitude for each time step as a horizontal projection of a respective radial velocity corrected with a respective horizontal derivative of a vertical velocity of an estimated velocity field determined for the respective time step; determine a mean of the unbiased horizontal velocities at each altitude for a time period comprising the set of time steps; and determine the turbulence at each altitude based on the unbiased horizontal velocities of each time step and the mean of the unbiased horizontal velocities; and an output interface configured to present the turbulence at each altitude.

[0016] Accordingly, another embodiment discloses a wind flow sensing method for determining a turbulence of a wind flow at a set of different altitudes above a terrain from a set of measurements of radial velocities at each altitude of the set of different altitudes, wherein the wind flow sensing method uses a processor coupled with stored instructions implementing the wind flow sensing method, wherein the instructions, when executed by the processor, carry out steps of the wind flow sensing method, the steps comprising: receiving, for a set of time steps, a set of measurements of radial velocities at a ground track point above the terrain for each altitude; estimating a velocity field for each altitude based on a data assimilation of the velocity field above the terrain, the data assimilation for fitting the measurements of radial velocities, wherein the velocity field is estimated for the set of time steps; estimating an unbiased horizontal velocity at each altitude for each time step as a horizontal projection of a respective radial velocity corrected with a respective horizontal derivative of a vertical velocity of an estimated velocity field determined for the respective time step; determining a mean of the unbiased horizontal velocities at each altitude for a time period comprising the set of time steps; and determining the turbulence at each altitude based on the unbiased horizontal velocities of each time step and the mean of the unbiased horizontal velocities; and outputting the turbulence at each altitude.

[0017] Embodiments of the present disclosure will be further explained with reference to the drawings. The depicted figures are not necessarily to scale and generally emphasize principles of the embodiments of the present disclosure. BRIEF DESCRIPTION OF DRAWINGS

[0018] [ Figure 1 ]

[0019] Figure 1 A schematic overview showing the principles used by some embodiments for fast wind flow measurement in complex terrain.

[0020] [ Figure 2 ]

[0021] Figure 2 A block diagram of a wind flow sensing system for determining wind flow according to some embodiments is shown.

[0022] [ Figure 3A ]

[0023] Figure 3A A schematic diagram of an exemplary remote sensing instrument configured to measure radial velocity of wind flow according to some embodiments is shown.

[0024] [ Figure 3B ]

[0025] Figure 3B A geometric schematic diagram showing radial velocity measured by some embodiments along a surface of a cone and along a centerline of a cone at particular altitudes is shown.

[0026] [ Figure 3C ]

[0027] Figure 3C A remote sensing schematic of wind over complex terrain used by some embodiments is shown.

[0028] [ Figure 4 ]

[0029] Figure 4 A schematic diagram showing exemplary wind flow parameters used by some embodiments to estimate a velocity field of a wind flow is shown.

[0030] [ Figure 5A ]

[0031] Figure 5A A block diagram of a computational fluid dynamics (CFD) framework for resolving a wind flow according to some embodiments is shown.

[0032] [ Figure 5B ]

[0033] Figure 5B A block diagram of a method for determining an unbiased velocity field according to one embodiment is shown.

[0034] [ Figure 6A ]

[0035] Figure 6A ​​​​​​​​​​​​​​​​​​A block diagram showing a framework for obtaining horizontal gradients of vertical velocity based on CFD simulation, according to one embodiment, is shown.

[0036] [ Figure 6B ]

[0037] Figure 6B An embodiment of a mesh determined by some embodiments is shown.

[0038] [ Figure 7 ]

[0039] Figure 7 A block diagram showing a method for selecting operating parameters, according to some embodiments, is shown.

[0040] [ Figure 8 ]

[0041] Figure 8 A flowchart showing a method for determining a current value of an operating parameter, according to some embodiments, is shown.

[0042] [ Figure 9 ]

[0043] Figure 9 A process for assigning different weights to different terms in a cost function, according to some embodiments, is shown.

[0044] [ Figure 10 ]

[0045] Figure 10 A schematic diagram showing an implementation of a direct accompanying loop (DAL) for determining operating parameters and CFD simulation results in an iterative manner, according to one embodiment, is shown.

[0046] [ Figure 11 ]

[0047] Figure 11 An embodiment of various data points on a single plane related to wind flow sensing, according to some embodiments, is shown.

[0048] [ Figure 12 ]

[0049] Figure 12 A schematic diagram showing a method for determining a horizontal velocity of a wind flow, according to one embodiment, is shown.

[0050] [ Figure 13 ]

[0051] Figure 13 A schematic diagram showing a CFD simulation, according to one embodiment, is shown.

[0052] [ Figure 14 ]

[0053] Figure 14 A map of the pressure and velocity field of the wind flow around a cylinder is shown, according to one embodiment.

[0054] [ Figure 15 ]

[0055] Figure 15 A geometry of a terrain flow model is shown, according to one embodiment.

[0056] [ Figure 16 ]

[0057] Figure 16 A combination of uniform flow and source flow is shown, according to one embodiment.

[0058] [ Figure 17A ]

[0059] Figure 17A A combination of uniform flow and dipole flow for determining fluid flow around a cylinder is shown, according to some embodiments.

[0060] [ Figure 17B ]

[0061] Figure 17B A sink flow and a source flow for obtaining a dipole flow of equal strength Λ, according to one embodiment, are shown.

[0062] [ Figure 18 ]

[0063] Figure 18 An exemplary map between a cylinder of radius b and a terrain is shown, according to some embodiments.

[0064] [ Figure 19 ]

[0065] Figure 19 A superposition of a set of cylinders for mapping with a terrain, according to some embodiments, is shown.

[0066] [ Figure 20 ]

[0067] Figure 20 A schematic diagram of constructing and evaluating a cost function including both LOS measurements and LOS according to Laplace superposition, according to some embodiments, is shown.

[0068] [ Figure 21 ]

[0069] Figure 21 A block diagram for an implementation of DAL to determine cylinder radius and upstream velocity, according to some embodiments, is shown.

[0070] [ Figure 22 ]

[0071] Figure 22 A schematic illustration of estimating the horizontal gradient of the vertical velocity is shown, in accordance with some embodiments.

[0072] [ Figure 23A ]

[0073] Figure 23A Collectively, a schematic overview of the principles used by some embodiments for turbulence measurements of wind flow in complex terrain is shown.

[0074] [ Figure 23B ]

[0075] Figure 23B Collectively, a schematic overview of the principles used by some embodiments for turbulence measurements of wind flow in complex terrain is shown.

[0076] [ Figure 24A ]

[0077] Figure 24A A schematic illustration of using the autocorrelation function to correct for standard deviation, in accordance with one embodiment, is shown.

[0078] [ Figure 24B ]

[0079] Figure 24B A schematic illustration of calculating values of the autocorrelation functions p u , p v , and p w based on comparison with anemometer data, in accordance with some embodiments, is shown.

[0080] [ Figure 24C ]

[0081] Figure 24C A schematic illustration of calculating values of the autocorrelation functions p u , p v , and p w based on comparison with high-fidelity CFD simulations, in accordance with some embodiments, is shown.

[0082] [ Figure 25 ]

[0083] Figure 25 A block diagram of a wind flow sensing system for determining turbulence of wind flow, in accordance with some embodiments, is shown.

[0084] [ Figure 26 ]

[0085] Figure 26 A schematic illustration of a wind turbine including a controller in communication with a system employing the principles of some embodiments is shown. DETAILED DESCRIPTION

[0086] In the following description, for purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the present disclosure. It will be apparent, however, to one skilled in the art that the present disclosure can be practiced without these specific details. In other instances, devices and methods are shown in block diagram form in order to avoid obscuring the present disclosure.

[0087] As used in the specification and claims, the terms "for example," "e.g.," and "such as," and the verbs "comprising," "having," "including," and their other verb forms, when used to describe the disclosure, are each meant to encompass the other mutatis mutandis. The term "based on" means at least partially based on. Furthermore, it is to be understood that the phraseology and terminology used herein is for the purpose of description and should not be regarded as limiting. Any use of class names in the present description is merely for convenience and has no legally interpretable effect.

[0088] Figure 1 A schematic overview of the principles used by some embodiments for fast wind flow measurement in complex terrain is shown. Remote sensing instruments, such as LiDARs, are used to measure a subset of the characteristics of the wind in the atmosphere. Different characteristics of the wind include wind speed (both horizontal and vertical rates), turbulence, wind direction, etc. In the method 100, the LiDAR measures the radial velocity of the wind 102 in the line-of-sight (LOS) direction. However, the horizontal velocity vector is a parameter of interest.

[0089] To this end, some embodiments are based on reconstructing the wind from the measured radial velocity 102 using geometric relationships 104. In other words, the horizontal velocity 106 is obtained by a horizontal projection of the measured radial velocity 102. In practice, such a projection is not accurate because the radial velocity is measured differently for different altitudes, and even for the same altitude, the five different radial velocities measured by the LiDAR have different values. In addition, the corresponding horizontal projection of the radial velocity does not take into account the rate variation in the vertical direction. Furthermore, such a projection is not valid for wind flow over complex terrain, such as near hills or large buildings or other urban structures.

[0090] To this end, some embodiments are based on the consideration of the target of the rate variation in the vertical direction, i.e., the vertical variation. Some embodiments are based on the recognition that the horizontal derivative of the vertical velocity can be used as a correction to the horizontal projection of the measured radial velocity to account for the vertical variation. In such embodiments, first, the vertical velocity is determined by data assimilation, which simulates a velocity field of the wind flow to find a velocity field that fits the measured values of the radial velocity 102. This simulation by considering the "closeness" of the simulated radial velocity to the measured values is referred to as data assimilation. In some embodiments, the data assimilation is implemented using computational fluid dynamics (CFD). Some embodiments are based on the recognition that both the horizontal velocity and the vertical velocity can be determined with the simulated velocity field. Moreover, the vertical velocity of the velocity field is used to estimate the corresponding horizontal derivative, which is then used as a correction to the horizontal projection of the measured radial velocity. As a result of this correction, the accuracy of the horizontal velocity estimation is improved.

[0091] However, in data assimilation, the operating parameters, such as the boundary or atmospheric conditions, are unknown and are determined iteratively until the operating parameters produce the measured radial velocity. Therefore, data assimilation with CFD is very time-consuming. Moreover, data assimilation with CFD is cumbersome because CFD simulation is an optimization process based on the solution of the Navier-Stokes equations. Furthermore, CFD simulation becomes complex for wind flow over complex terrain.

[0092] To this end, some embodiments are based on the representation or approximation 110 of the complex terrain with convex shapes, e.g., cylinders. In some embodiments, the complex terrain is approximated with an equivalent cylinder. In some other embodiments, the terrain is approximated with multiple convex shapes. Such representation simplifies the simulation of the velocity field. Moreover, the wind flow around such cylinders is approximated with potential flow 112. Potential flow 112 involves an algebraic solution of the Laplace equation, rather than an iterative optimization of the Navier-Stokes equations, thus improving the computational efficiency. Moreover, the Laplace equation 112 is easier to solve compared to the Navier-Stokes equations, and in the case of simple shapes, e.g., cylinders, there exist closed mathematical form analytical solutions. However, such approximation 110 degrades the velocity field simulation. Some embodiments are based on the recognition that while such approximation 110 can not be sufficiently accurate for the determination of the velocity field, the approximation 110 can be sufficiently accurate for determining the horizontal derivative of the vertical velocity as a correction to improve the horizontal projection of the measured radial velocity 108. Therefore, the accuracy of the horizontal velocity estimation 114 is significantly improved with minimal degradation of the velocity field.

[0093] To this end, some embodiments are based on the recognition that data assimilation is implemented based on fitting the measurements of radial velocity 108 with a convex-shaped topographic approximation 110 to estimate the velocity field. Further, the horizontal velocity is estimated 114 as the horizontal projection of the corresponding radial velocity corrected with the corresponding horizontal derivative of the vertical velocity of the estimated velocity field. Additionally or alternatively, embodiments based on this conception can perform wind reconstruction and / or compute the horizontal velocity on-line (i.e., in real-time).

[0094] Figure 2 A block diagram of a wind flow sensing system 200 for determining wind flow in accordance with some embodiments is shown. The wind flow sensing system 200 includes an input interface 202 for receiving a set of measurements 218 of radial velocity in the site line direction for each elevation. In some embodiments, the measurements 218 are measured by a remote sensing instrument (e.g., a ground-based LiDAR) over a cone. The wind flow sensing system 200 can have interfaces that connect the system 200 with other systems and devices. For example, a network interface controller (NIC) 214 is adapted to connect the wind flow sensing system 200 to a network 216 via a bus 212 that connects the wind flow sensing system 200 with a remote sensing instrument configured to measure radial velocity of wind flow. With the network 216, whether wireless or wired, the wind flow sensing system 200 receives the set of measurements 218 of radial velocity in the site line direction for each elevation.

[0095] Further, in some embodiments, with the network 216, the measurements 218 can be downloaded and stored within a storage system 236 for further processing. Additionally or alternatively, in some embodiments, the wind flow sensing system 200 includes a human-machine interface 230 that connects the processor 204 to a keyboard 232 and a pointing device 234, which can include a mouse, trackball, touchpad, joystick, pointing stick, stylus, or touch screen, among others.

[0096] The wind flow sensing system 200 includes a processor 204 configured to execute stored instructions and a memory 206 that stores instructions executable by the processor. The processor 204 can be a single core processor, multi-core processor, computing cluster, or any number of other configurations. The memory 206 can include random access memory (RAM), read only memory (ROM), flash memory, or any other suitable memory systems. The processor 204 is connected to one or more input and output interfaces and / or devices via a bus 212.

[0097] According to some embodiments, the instructions stored in the memory 206 implement a method for determining a wind flow rate field for a set of different altitudes from a set of measurements of radial velocity at each altitude. To this end, the storage device 236 can be adapted to store different modules storing executable instructions for the processor 204. The storage device 236 stores a CFD simulation module 208 configured to estimate a rate field for each altitude by simulating a computational fluid dynamics (CFD) of the wind flow with current values of the operating parameters. The storage device 236 also stores a CFD operating parameters module 210 configured to determine operating parameter values that reduce a cost function and a horizontal derivative module 236 configured to determine a horizontal derivative of the vertical velocity of the rate field. In addition, the storage device 236 stores a rate field module 238 configured to determine a rate field including a horizontal velocity using the horizontal derivative of the vertical velocity and the measurements of the radial velocity. Furthermore, the storage device 236 stores a Laplace simulation module 240 configured to approximate a shape of the terrain with one or more of a set of convex shapes to fit the measurements of the radial velocity. The Laplace simulation module 240 is configured to solve a plurality of Laplace equations defining the wind flow dynamics for a particular value of the inlet rate field and a radius of the convex shape to approximate the shape of the terrain. The storage device 236 can be implemented using a hard disk, an optical drive, a thumb drive, an array of drives, or any combination thereof.

[0098] The wind flow sensing system 200 includes an output interface 224 for presenting the estimated horizontal velocity for each altitude. In addition, the wind flow sensing system 200 can be linked by means of the bus 212 to a display interface 220 adapted to connect the wind flow sensing system 200 to a display device 222, which can be a computer monitor, a camera, a television, a projector, or a mobile device, among others. In addition, the wind flow sensing system 200 includes a control interface 226 configured to submit the estimated horizontal velocity for each altitude to a controller 228 integrated with a machine, such as a wind turbine. The controller 228 is configured to operate the machine based on the estimated horizontal velocity for each altitude. In some embodiments, the output interface 224 is configured to submit the estimated horizontal velocity for each altitude to the controller 228.

[0099] Figure 3A A schematic diagram of an exemplary remote sensing instrument configured to measure the radial velocity of the wind flow according to some embodiments is shown. The LiDAR 300 is configured to measure the radial velocity 218 of the wind flow at different altitudes. Different embodiments use different remote sensing instruments. Embodiments of these instruments include radar, LiDAR, and SODAR. For the sake of clarity, the present disclosure uses the LiDAR 300 as an exemplary remote sensing instrument.

[0100] The radial velocity of an object relative to a given point is the rate of change of the distance between the object and that point. That is, the radial velocity is the component of the object's velocity in the direction pointing towards the radius connecting the object and the point. In the case of atmospheric measurements, the point is the location on Earth of a remote sensing instrument (such as radar, LiDAR, and SODAR), and the radial velocity represents the speed at which an object moves away from or towards the receiving instrument (LiDAR device 300). This measured radial velocity is also known as the line-of-sight (LOS) velocity.

[0101] Remote sensing instruments determine the flow of fluid (e.g., air) in a volume of interest by describing the velocity field of the airflow. For example, LiDAR 300 includes: a laser 302 or acoustic transmitter and receiver in which a returned signal 306 is spectrally analyzed; a computer 304 for further calculations; and a navigator for aiming the transmitter and / or receiver at a target in space at a considerable distance from the transmitter and receiver. The receiver detects the returned signal 306 scattered along a measurement axis due to contaminants present between the remote sensing system and the target. The laser propagates along a conical surface 308 formed by possible aiming directions. The radial velocity of particles in the volume of interest 310 at the target is inferred from the frequency shift caused by the Doppler effect due to specific air contaminants.

[0102] Figure 3B A geometrical schematic diagram of radial velocities measured at specific altitudes along cone surface 308 and cone centerline 312, by several implementations, is shown. LiDAR measurements provide the radial (line-of-sight) velocity component of the wind, which is difficult to determine precisely due to the so-called "Cyclops" dilemma. This phenomenon refers to the inability to accurately reconstruct an arbitrary three-dimensional velocity field using a single LOS measurement. The radial velocity 314 along a beam shows the projection of the velocity vector 316, with the LiDAR located at position 318 in Cartesian coordinate systems 320, 322, and 324.

[0103] Here, θ320 is the horizontal wind direction measured clockwise from 326 North, ψ328 is the elevation angle of the beam, and (u,v,w) are the x 322, y 320, and z 324 components of the wind speed V at each point in space.

[0104] Horizontal rate v at each altitude h Limited to:

[0105]

[0106] The radial rate (also called the LOS rate) is defined at each altitude as follows:

[0107] v R= u sin 0 sin y + v cos 0 sin y + w cos y Equation 2

[0108] Figure 3C A remote sensing schematic of the wind over a complex terrain 330 is shown for use in some embodiments. A LiDAR 300 disposed at a point 332 (e.g., a hilltop) makes a series of line-of-sight measurements over a cone, including measurements 334, 336, 340, 342 along the surface of the cone and a measurement along the centerline 338. These measurements are made for different elevations 344 shown as different planes. In this way, for each elevation 344, the measurements over the cone are measurements over a circle, including a plurality of measurements of radial velocity in different angular directions measured at different line-of-sight points over the circumference and one measurement of radial velocity in the vertical direction measured at the center of the circle. The line-of-sight measurements correspond to line-of-sight velocities.

[0109] One embodiment aims to determine the horizontal velocity v h of the wind flow at each elevation. Given these measurements, an estimate of the horizontal velocity v R can be determined from the measurements of the radial velocity v h using geometric relationships and assuming that the wind velocity is uniform over each plane. Here V L = (u L , v l , w L ) is the estimated velocity of the uniform wind flow.

[0110] For example, the following equations can be derived for the estimated velocity from the radial velocities V1, V2, V3, V4, V5 corresponding to the beams pointing north, east, south, west, and centerline.

[0111]

[0112] Some embodiments are based on the recognition that for complex terrain (such as terrain 330), the uniform velocity assumption leads to a bias in the LiDAR estimate of the horizontal velocity. The main error is due to the variation in the vertical velocity w in the vertical direction (e.g., along the hill). For this reason, some embodiments are based on the recognition that the uniform velocity assumption can be corrected using the horizontal derivative of the vertical velocity when sensing the wind flow over complex terrain.

[0113] Figure 4 A schematic of exemplary parameters of the wind flow for use in some embodiments to estimate the wind flow velocity field is shown. Some embodiments are based on the recognition that the uniform velocity assumption is incorrect when sensing the wind flow over complex terrain, but can be corrected using the horizontal derivative of the vertical velocity. Figure 4A two-dimensional plot of the wind flow is shown when the LiDAR 300 is placed near the top of a hill (e.g., at position 332). The horizontal derivative of the vertical velocity can show the change in direction and / or magnitude of the vertical velocity 400 at a given altitude. In this embodiment, the horizontal derivative of the vertical velocity shows that the vertical velocity increases along one slope of the hill until the hilltop at point 402, and decreases along the other slope of the hill starting at point 404. Additionally, the precision of the wind flow sensing can be improved using the first derivative showing a linear change in the vertical velocity, as the uniform velocity assumption leads to a first order term of the sensing velocity field error.

[0114] For any point of altitude z above the device 300, its error or bias can be written as a first order:

[0115]

[0116]

[0117] Thus, the bias due to the uniform assumption is proportional to i) the altitude z above the device 300; ii) the horizontal gradients of the vertical velocity dw / dx and dw / dy. This error is not a function of the elevation angle ψ, reducing this angle does not reduce the bias of the horizontal velocity. Some embodiments are based on the recognition that it can not be possible to obtain estimates of dw / dx and dw / dy based on radial velocity measurements alone. The resulting equations are indeterminate due to the symmetry of the scanning beam.

[0118] The incompressibility of a flow refers to a fluid whose density is constant within the material parcel, which is an infinitesimal volume moving at the flow rate. This physical principle is based on the conservation of mass. Some embodiments are based on the recognition that the first order error caused by the uniform velocity assumption is incompressible. In other words, it can be shown that the bias term consisting of the product of the altitude and the horizontal gradient of the vertical velocity is a conservation of mass. This means that, for a wind flow over a complex terrain, enforcing the incompressibility condition on the fluid volume within the field of interest does not correct the first order error caused by the uniform flow assumption.

[0119] Computational fluid dynamics (CFD) is a branch of fluid mechanics that uses numerical analysis and data structures to solve and analyze problems involving fluid flow. Computers are used to perform the calculations needed to simulate the interaction of liquids and gases with surfaces defined by boundary conditions. Some embodiments are based on the general understanding that CFD can be used to estimate the velocity field of the wind from measurements of the wind on the cone sensed by the LiDAR. However, the operating parameters (such as boundary conditions) of the wind flow over a complex terrain are generally unknown, and approximations of these operating parameters can undesirably reduce the precision of the wind flow sensing.

[0120] Some implementations are based on the recognition that while the CFD approximation may not be precise enough for determining the velocity field, the average of the horizontal derivatives of the vertical velocity reconstruction at a given altitude can be sufficiently accurate. This average can then be used to correct for biases arising from the uniform velocity assumption. To this end, some implementations use the CFD approximation to determine the horizontal derivative of the vertical velocity and combine it with a radial velocity measurement of the wind flow at the desired altitude to determine the velocity field at that altitude. In this way, target accuracy can be achieved using the radial velocity measurements for velocity field sensing.

[0121] Figure 5A A block diagram of a computational fluid dynamics (CFD) framework for analytical airflow used in some implementations is shown, with the aim of obtaining accurate measurements of the horizontal gradient of the vertical velocity at each altitude of interest. Using line-of-sight measurements 500, a first approximation of the velocity field 502 is obtained from CFD simulations. For example, the CFD simulation can be performed by... Figure 2 The CFD simulation module 208 shown is used for this purpose. In some cases, CFD simulations require operating parameters, such as boundary and atmospheric conditions. These operating parameters are often unknown. Therefore, some implementations determine operating parameters that reduce the difference between the estimated radial velocity and the measured radial velocity. For example, this estimation can be performed by… Figure 2 The CFD operation parameter module 210 shown is used.

[0122] Although the velocity field in the first approximation provided by CFD is imprecise for the desired purpose, an estimate 504 of the horizontal gradient of the vertical velocity can be extracted with the required accuracy. This extraction can be performed by module 236. The CFD simulation yields the velocity field at discrete points of the grid. Using this velocity field, the x and y derivatives at each discrete point are calculated using the finite difference method. Then, by averaging the derivatives in the x and y directions on the corresponding planes, the x and y horizontal derivatives of the vertical velocity are extracted for each plane. A single value. Then, the horizontal gradient of this vertical rate, together with the geometric relationship between the line-of-sight rate and the wind rate, is used to correct the biased horizontal rate component u based on the uniformity assumption using equations (3a) and (3b). L and v L This estimation can be performed by module 238.

[0123] Figure 5B A block diagram is shown of a method for determining an unbiased velocity field according to one embodiment. To determine a second approximation of the velocity field, this embodiment determines a biased velocity field 508 under the assumption of a uniform velocity field at each altitude, using the horizontal derivative of the vertical velocity at the corresponding altitude. Removing the bias of the uniform velocity assumption of the biased velocity field at each elevation 510.

[0124] For example, equations 3a and / or 3b are used to remove the bias term by subtracting the bias term and from the biased velocity field u L , v L to obtain the unbiased velocity field (u, v).

[0125] Figure 6A A block diagram showing the framework of CFD simulation based implementation of some embodiments for obtaining the horizontal gradient of the vertical velocity 608. This implementation proceeds with a pre-processing step to define the geometry and physical boundaries of the CFD simulation. For example, some implementations use computer aided design (CAD) to define the simulation range. The volume occupied by the fluid (wind) is divided into discrete cells (mesh). GPS is used to extract the geographical location of the terrain 600. This location is compared with available datasets stored in the device memory to generate terrain data. Terrain data can be collected using various resources such as Google or NASA databases. Furthermore, the optimal radius is selected to construct the mesh 610.

[0126] Figure 6B An example of a mesh 610 determined by some embodiments is shown. In various embodiments, the mesh 610 can be uniform or non-uniform, regular or irregular, composed of hexahedral, tetrahedral, prismatic, pyramidal, or a combination of polyhedral elements. The optimal mesh size and number are selected so as to capture important terrain structures in the mesh based on the wind direction. The mesh is generated based on the radius selected on the terrain. Furthermore, the resolution of the mesh is manually set.

[0127] During pre-processing, the values of the operating parameters 604 are also specified. In some embodiments, the operating parameters specify the fluid behavior and properties of all boundary surfaces of the fluid domain. Boundary conditions (inlet velocity) for the field (velocity, pressure) specify the value of the function itself, or the value of the normal derivative of the function, or the form of a curve or surface that gives the value of the normal derivative and the value of the variable itself, or the relationship of the value of the function and the derivative of the function in a given region. The boundary condition at a solid surface defined by the terrain involves the velocity of the fluid, which can be set to zero. The inlet velocity is decided based on the wind direction and velocity, which has a logarithmic curve with respect to height on flat terrain.

[0128] Some embodiments perform a CFD simulation by solving 606 one of the variables of the Navier-Stokes equations for the wind flow with the current value of the operating parameter. For example, CFD solves the Navier-Stokes equations with mass and energy conservation. This set of equations has been shown to represent the mechanical behavior of any Newtonian fluid (e.g., air) and is used for the simulation of atmospheric flows. The discretization of the Navier-Stokes equations is a reformulation of the equations that makes them applicable to computational fluid dynamics. The numerical method can be finite volume, finite element, or finite difference methods as well as all spectral or spectral element methods.

[0129] The control equations (Navier-Stokes equations) are shown below:

[0130]

[0131]

[0132] is the divergence operator. is the gradient operator, is the Laplace operator. Equation 4 can also be extended to the transient case where the variation of the velocity and pressure with time is taken into account.

[0133] Some embodiments express equations 4a and 4b as N(p, V) = 0, where the inlet velocity and direction are expressed by V in , Θ in ; p: air pressure [pa] or [atm], p: air density [kg / m 3 ], v: kinematic viscosity [m 2 / s]. After the CFD simulation, some embodiments extract 608 the horizontal gradient of the vertical velocity.

[0134] Figure 7 A block diagram of a method for selecting operating parameters according to some embodiments is shown. For example, some embodiments select 706 operating parameters 700 based on the sensitivity 702 of the horizontal derivative of the vertical velocity (HDVV) to the variation of the operating parameter value. In one embodiment, among the set of purpose-based operating parameters approximated during the CFD simulation, the operating parameters with a sensitivity higher than a threshold 704 are selected. In this way, some embodiments adapt the unknown operating parameters of the CFD to the purpose of the approximation. Such an operating parameter adaptation can reduce the computational burden without reducing the accuracy of the CFD approximation of the quantities of interest. For example, some embodiments select operating parameters such as terrain roughness, inlet mean velocity, inlet turbulence intensity, and atmospheric stability conditions.

[0135] In some embodiments, the operating parameters include inlet boundary conditions (velocity, direction), surface roughness, and atmospheric stability. In one embodiment, the operating parameters are selected to be inlet boundary conditions (velocity, direction), surface roughness, inlet turbulent kinetic energy, and dissipation. Values for such operating parameters cannot be obtained directly from LiDAR measurements.

[0136]

[0137]

[0138]

[0139]

[0140] C μ is a constant in the k-ε turbulence model,

[0141] κ von Karman constant,

[0142] V * friction velocity [m / s],

[0143] V ref is a reference velocity [m / s] selected at a reference location, which can be arbitrary,

[0144] z ref is a reference elevation [m].

[0145] z0surface roughness.

[0146] Turbulent kinetic energy is the kinetic energy per unit mass of turbulent fluctuations. Turbulent dissipation ε is the rate at which turbulent kinetic energy is converted to heat internal energy.

[0147] Some embodiments are based on the recognition that, in some cases, operating parameters for simulating CFD are unknown. For example, for the case described above, V ref , z ref , z0are unknown operating parameters, and remote sensing measurements do not directly provide these values.

[0148] Figure 8 A flowchart showing a method for determining current values of operating parameters according to some embodiments is shown. In particular, some embodiments determine 802 operating parameters that minimize the error between measurements 800 of radial velocity at a set of site line points and estimates 804 of radial velocity at the same set of site line points by CFD with current values of the operating parameters.

[0149] Some implementations are based on the recognition that, when the CFD is used to extract the horizontal derivative of the vertical velocity, a particular cost function 806 is minimized to obtain the operational parameter estimate. In particular, some implementations are based on the recognition that the horizontal derivative of the vertical velocity depends on the altitude having different effects on the velocity field. To this end, the cost function 806 includes a weighted combination of errors. Each error corresponds to an altitude and includes a difference between a measured velocity at a site line point of the corresponding altitude and a simulated velocity at the site line point of the corresponding altitude as simulated by the CFD using current values of the operational parameters. Moreover, the weights of at least some of the errors are different. For example, the errors include a first error corresponding to a first altitude and a second error corresponding to a second altitude, where the weight of the first error in the weighted combination of errors is different from the weight of the second error in the weighted combination of errors.

[0150] Figure 9 A process of assigning different weights to different terms in the cost function is shown in accordance with some implementations. In some embodiments, the cost function 806 returns a number that represents how well the CFD simulation (the line-of-sight velocity according to the CFD simulation) 902 matches the LiDAR data (the line-of-sight velocity from the LiDAR measurements) 900 along the line of sight of different beams at different altitudes. To this end, in order to determine the horizontal derivative of the vertical velocity by the CFD, the cost function considers different altitudes in different ways, e.g., with different weights 930. For example, in some implementations, the cost function includes a weighted combination of errors that represent the accuracy of the CFD for different altitudes.

[0151] In one implementation, the cost function is:

[0152]

[0153] vi is the line-of-sight velocity at the location of point i, v R,i is the line-of-sight velocity at the location of point i, v R,CFD is the radial velocity according to the CFD simulation operation at the location of point i, w i is a weighting coefficient. The error of each term is proportional to the difference between the measured and the CFD radial velocity. To give more weight to the estimate of the vertical velocity gradient for higher altitudes, some implementations set the weighting coefficient w i to be proportional to the altitude (i.e., the height above the device location). For example, v R,CFD is a set of radial velocities obtained according to the CFD simulation of the wind flow to produce a first approximation of the velocity field that reduces the cost function of the weighted combination of errors given in equation (6).

[0154] vi is the line-of-sight velocity at the location of point i, v RThe set represents the measurements of the radial (or line-of-sight) velocity given by the remote sensing instrument of the wind flow. The error of this value is very small and is used as the true value of the wind in the direction of the beam. Each term in equation (6), denoted by i, corresponds to the error induced by one elevation and includes the measured velocity v R the difference between the simulated velocity at the site line point for the corresponding elevation by the CFD simulation. The weight of each error in the weighted combination of errors is an increasing function of the corresponding elevation value.

[0155] Some embodiments are based on the recognition that unknown values of operational parameters can be estimated using a Direct-Adjoint Loop (DAL) based CFD framework. This framework enables simultaneous correction of unknown parameters serving a common purpose by minimizing a cost function (i.e., the error of the estimated line-of-sight data and its gradient between the forward CFD simulation and the available LiDAR measurements) and then solving the sensitivity (or adjoint-CFD) equations in an iterative manner. The sensitivity of the parameters serving a common purpose indicates the direction of convergence of the DAL based CFD framework. The simultaneous correction reduces the computational time to update multiple operational parameters.

[0156] Some embodiments represent the set of operational parameters to be estimated as (ξ1, ξ2,... ξ n ). Then, the sensitivity of the cost function J with respect to any operational parameter ξ i can be represented as

[0157]

[0158] Figure 10 A schematic diagram showing the implementation of DAL for determining the results of the operational parameters and the CFD simulation in an iterative manner is shown. This embodiment estimates the most likely values of the operational parameters by evaluating the CFD simulation. DAL is an optimization method that solves the CFD equation 1002 and the adjoint (or sensitivity) equation 1004 in an iterative 1014 manner to obtain the sensitivity 1006 of the unknown operational parameters with respect to the current estimate of the operational parameters. DAL is initialized 1000 with a guess or initial estimate of the operational parameters. For example, the inlet velocity is estimated using the Bernoulli equation and the angle is estimated using a uniform assumption. After each iteration, the estimate of the current value of the operational parameters is updated 1008 using a conjugate gradient descent to update in the direction of the maximum reduction of the sensitivity of the cost function. To do this, the CFD simulation is performed multiple times, i.e., once per iteration, and if the change in the estimate from the last iteration is below a threshold, the DAL method is considered to be converged 1010. The DAL method is obtained by formulating the Lagrangian equation:

[0159] L = J + ∫ Ω (p aV a )N(p,V)dΩ Equation 8

[0160] Since N(p,V) = 0 in the equation (Navier-Stokes equation), the equations L and J are equal when the values of p and V are exact. Considering the variation of ξ i , the variation of L can be expressed as:

[0161]

[0162] To determine the term , the adjoint variable is chosen to satisfy:

[0163]

[0164] Thus, the DAL method involves new variables (V a ,p a ), which represent the adjoint velocity and pressure, respectively, such that is computable.

[0165] In one embodiment, the unknown parameters are chosen to be V in ,Θ in , i.e. the inlet velocity and the inlet angle. Thus, the problem of finding V in ,Θ in that minimizes J is converted to the problem of finding V in ,Θ in that minimizes the augmented objective function L. For example, to determine δJ / δV in and δJ / δΘ in , the DAL method can be used by setting ξ i = V in or ξ i = Θ in

[0166] The adjoint equations in step 1004 are given by:

[0167]

[0168]

[0169] The operator corresponds to the transpose of the gradient of the velocity vector.

[0170] The adjoint variable can be used to determine the sensitivity of the cost function to any operational parameter 1006

[0171]

[0172] For example, equation 11 can be taken with respect to the cost function with respect to the inlet velocity V​in The sensitivity of A is written as:

[0173]

[0174] A in : Entrance area of the operating domain Ω [m 2 ]

[0175] n: Unit normal vector of A [m in ] 2

[0176] By using the gradient descent algorithm, the estimate of the operating parameter ξ i can be updated 1008 to:

[0177]

[0178] λ is a positive number representing the step size, which can be chosen using some standard algorithm. Using the DAL method, only one solution of equations (4) and (10) is solved at each iteration regardless of the number of unknown parameters, thus reducing the computational cost and making the solution of the optimization problem feasible. This is an advantage of the companion method compared to the method of determining the sensitivity of the cost function by directly measuring the perturbation of the cost function. After the DAL converges to produce the current value of the external operating parameter 1012, some embodiments extract the quantity of interest (i.e., the vertical rate gradient) to correct the bias error of the wind speed reconstruction on complex terrain using LiDAR Line of Sight (LOS) on the measurement cone.

[0179] Figure 11 An example of various data points on a single plane related to wind flow sensing is shown in accordance with some embodiments. In this example, the points on the circle 1100 are points of the measured radial rate 1102. In some embodiments, the rate field at each altitude includes rate values of the wind inside and outside of the cone (i.e., the circle 1100).

[0180] Additionally or alternatively, in some embodiments, the horizontal derivative of the vertical rate at each altitude defines the gradient of the vertical rate at the center of the cone that defines the measurement of the LiDAR at the respective altitude. For example, some embodiments average the rate and / or the gradient at each altitude to produce the center 1104 of the cone and the circle 1100. In these embodiments, the second approximation of the rate field obtained via the geometric relationship and the removal of the bias using the horizontal gradient of the rate provides a single value of the rate field on each plane (or each altitude). In this manner, the unbiased rate values at 1106 and 1108 are considered to be equal to this single value.

[0181] ​To this end, in one embodiment, the second approximation of the velocity field includes a single value of the velocity field per elevation. Furthermore, this embodiment converts the single value to a dense grid of non-constant values of the velocity field per elevation by enforcing the incompressibility of the wind flow and regularization consistent with the measured values of the radial velocity at each elevation. After such conversion, the horizontal velocities of the points inside and outside the cone (e.g., points 1102, 1104, 1106, and 1108) can have different values.

[0182] Figure 12 A schematic diagram for determining the horizontal velocities of the wind flow is shown, according to one embodiment. This embodiment starts with the same unbiased values for all points on a single plane 1200, enforces the incompressibility of the air to reduce the error caused by the sparsity of the LOS measurements 1202 and the second order error caused by the uniformity assumption to produce a dense grid of velocity values 1204. In various embodiments, the density of the grid of points is a user-specified value. Notably, in this embodiment, the incompressibility of the air is used to correct the second order error, as opposed to enforcing the incompressibility to correct the first order error caused by the uniform velocity assumption.

[0183] In some embodiments, the dense grid of non-constant values of the velocity field is determined using a Direct-Adjoint-Loop based algorithm. This algorithm first interpolates the unbiased velocity values at all discrete points on the grid in each plane 1200. The DAL problem is formulated to enforce the incompressibility in the volume occupied by the fluid 1202 while minimizing a cost function that has two terms: one is the difference between the final velocity field and the initial velocity field at the measured discrete points; the other is a regularization term to increase the smoothness of the velocity field. The adjoint equation for the adjoint variable λ is given by:

[0184]

[0185] where U k is the velocity field at the kth iteration of the DAL loop. At the end of each iteration, the update is made as follows:

[0186]

[0187] The algorithm terminates when convergence is reached.

[0188] When solving the Navier-Stokes equations, the computational cost depends on the velocity and viscosity of the fluid. For atmospheric flows, the computational cost is very large because the wind velocity is high and the viscosity of air is small. This leads to so-called high Reynolds number flows, for which the unstable inertial forces in the flow are significantly larger than the stabilizing viscous forces. In order to fully resolve the fluid dynamics problem and avoid numerical instabilities, all spatial scales of the turbulence are resolved in the computational grid, from the smallest dissipative scale (Kolmogorov scale) up to the integral scale, which is proportional to the domain scale associated with the motions containing most of the kinetic energy.

[0189] Large Eddy Simulation (LES) is a popular CFD technique for solving the governing equations of fluid mechanics. The implication of Kolmogorov’s self-similarity theory is that the large eddies of the flow depend on the geometry, while the smaller scales are more universal. This feature allows the large eddies to be explicitly solved in the computation and the small eddies to be implicitly considered by using a sub-grid scale model (SGS model). CFD simulations using the LES approach can simulate the flow field with high fidelity, but the computational cost is very expensive.

[0190] Some embodiments are based on the recognition that rather than using a high-fidelity CFD solution for every new set of measurement data (e.g., for every new wind direction and / or new terrain), it is preferable to modify a low-fidelity model to learn the internal model parameters required for the desired accuracy in the results.

[0191] Figure 13 A schematic diagram of a CFD simulation is shown, according to one embodiment. A low-fidelity CFD simulation approximates the small scale terms in the flow by means of a model that depends on some internal parameters. For this purpose, this embodiment uses a feature vector that includes the horizontal derivatives of the vertical velocity to apply a Field Inversion and Machine Learning (FIML) method 1302 to learn the dependence of the internal parameters of the low-fidelity model 1303 on the flow features in a high-fidelity LES simulation 1300.

[0192] A low-fidelity CFD model is used, such as the Reynolds-averaged Navier-Stokes equations (or RANS equations). Such a model is a time-averaged equation of motion for a fluid flow. In addition, this model includes internal parameters to approximate terms that cannot be analytically resolved due to the low-fidelity. The correct values of the internal parameters are problem-specific, so in order to make the RANS almost as accurate as the LES, a FIML framework is employed. The significant advantage of using RANS in combination with FIML is that the computational cost of the CFD simulation at high Reynolds numbers is reduced by several orders of magnitude compared to the high-fidelity LES simulation, while maintaining the required accuracy. Once the internal model parameters of the low-fidelity model are fixed offline (in advance), a RANS-based CFD simulation can be performed if the operating parameters are known.

[0193] In this way, in some embodiments, the CFD simulation of the wind flow is performed by solving the Reynolds-averaged Navier-Stokes (RANS) equations, while the internal operating parameters of the RANS equations are determined using Field Inversion with Eigenvectors and Machine Learning (FIML) that utilizes eigenvectors that include the horizontal derivative of the vertical velocity at each elevation.

[0194] In case of using uniform design to compute the horizontal velocity, the relative error of LiDAR to cup anemometer is about 8%, while in case of using CFD and DAL to find the most feasible operating parameters (focusing on the inlet velocity and wind direction), the relative error of LiDAR to cup anemometer is about 1%. Moreover, some embodiments enforce the incompressibility assumption to reconstruct the dense field inside and outside the conical region.

[0195] To this end, it can be appreciated that data assimilation with CFD is very time consuming, because in the CFD simulation, the operating parameters are unknown, and the operating parameters are determined iteratively until the operating parameters produce the measured radial velocity. Moreover, data assimilation with CFD is cumbersome, because the CFD simulation is an optimization process based on the solution of the Navier-Stokes equations. Furthermore, the CFD simulation becomes complex for wind flow over complex terrain.

[0196] Approximation of data assimilation

[0197] To this end, some embodiments are based on representing or approximating the complex terrain with a convex shape (e.g., a cylinder). Some embodiments are based on the appreciation that the wind flow around a cylinder is similar in nature to the wind flow over complex terrain.

[0198] Figure 14 A mapping 1400 of the pressure 1404 and velocity field 1406 of the wind flow around a cylinder 1402 is shown, in accordance with one embodiment. On the surface of the cylinder 1402, the pressure varies, in particular, the pressure is maximum at 1408 and minimum at 1410. Due to this reason, the variation of the wind velocity results in a non-uniform velocity field 1406. In other words, there is a variation / gradient of the vertical velocity in the horizontal direction dw / dx, which can be positive or negative (dw / dx > 0 or dw / dx < 0) based on the location.

[0199] In some embodiments, the wind flow around such convex shape is approximated with potential flow. Potential flow involves the Laplace equation as the governing equation, instead of the Navier-Stokes equation as the governing equation. Thus, the Navier-Stokes equation is replaced by the Laplace equation:

[0200]

[0201] For potential flow, in the Laplace equation, the velocity is expressed in terms of a velocity potential The wind flow is incompressible, so This makes it possible to

[0202] Since the potential flow involves an algebraic solution of the Laplace equation, a simple shape (e.g., a cylinder) is chosen, rather than an iterative optimization of the Navier-Stokes equation, thus improving the efficiency of the operation. Even the iterative solution of the Laplace equation is much cheaper in terms of operation than the iterative solution of the Navier-Stokes equation.

[0203] Figure 15 A geometry of a topographic flow model is shown in accordance with one embodiment. A uniform horizontal wind flow 1500 with a velocity of U in the +ζ direction is incident on a cylinder 1502. The wind velocity U can also correspond to an upstream velocity. The wind flow can also be referred to as a fluid flow. The cylinder 1502 has a radius of "a" 1504 centered at the origin in a rectangular Cartesian coordinate system (ζ, y, η). The potential flow solution yields a stream function

[0204]

[0205] In addition, the streamlines are defined by the following equation:

[0206]

[0207] The horizontal velocity component in Cartesian coordinates is given by the following equation:

[0208]

[0209]

[0210] Some embodiments are based on the recognition that, in a fluid flow, the flow field can be defined by a velocity potential or a stream function that satisfies the Laplace equation. Since the Laplace equation is linear, various solutions can be added to obtain a desired solution. For example, for a linear partial differential equation such as the Laplace equation, the solution for various boundary conditions is the sum of the individual boundary conditions. In a flow field, a streamline can be thought of as a solid boundary because no flow crosses it. In addition, the conditions along the solid boundary and the streamline are the same. Thus, a combination of a velocity potential and a stream function for a basic potential flow results in the shape of a particular body that can be interpreted as a fluid flow around the body. The method of solving such potential flow problems is referred to as superposition.

[0211] To this end, some embodiments are based on the recognition that a potential flow around a cylinder can be determined by a combination of a velocity potential and a stream function for a basic potential flow. The basic potential flows include a uniform flow, a source / sink flow, a dipole flow, etc.

[0212] Figure 16A combination of uniform flow 1600 and source flow 1602 is shown according to one embodiment. The uniform flow 1600 is a uniform flow with a velocity V ∞ The uniform flow 1600 can be defined with a potential function as follows:

[0213]

[0214] The resulting stream function for the uniform flow 1600 is given as:

[0215] Ψ uniform = V ∞ r sin(θ) Equation 15

[0216] In a fluid flow, all streamlines are straight lines converging or diverging from a central point O 1604, and if the flow is converging towards the central point O 1604, it is called a sink flow. Conversely, if the flow is diverging from the central point 1604, the flow is called a source flow 1602. The velocity field resulting from the above flow includes a radial component Vr, which is inversely proportional to the distance from the point O. The potential flow of the source flow 1602 is given by the stream function as:

[0217]

[0218]

[0219] where Λ is the source strength, is the volumetric flow rate from the source, and

[0220] r is the distance from O. A positive Λ value refers to a source flow 1602, while a negative Λ value refers to a sink flow.

[0221] In addition, the source flow 1602 of strength Λ is superimposed with the uniform flow 1600 to produce a combined flow 1606. The resulting stream function from this can be given as

[0222]

[0223] The streamlines of the combined flow 1606 result in fluid flow over a semi-infinite body / shape and are obtained as follows:

[0224]

[0225] The velocity field is obtained from the stream function by differencing in polar coordinates, i.e.:

[0226]

[0227] In the combined flow 1606, the rate due to the source flow 1602 cancels the rate of the uniform flow 1600 at some point, and the flow is stagnant at that point. Some embodiments are based on the realization that the stream line 1608 contains the stagnation point at "B" and separates the flow from the uniform flow 1600 from the flow from the point 1604. The fluid flow outside the stream line 1608 is from the uniform flow 1600, while the fluid flow inside the stream line 1608 is from the source flow 1602. In fluid flow, the rate at the surface of a body is tangential to the body. For this reason, some embodiments are based on the realization that any stream line of the combined flow 1606 can be replaced by a solid surface of the same shape. Thus, in terms of the uniform flow 1600, if the stream line 1608 is replaced by a solid, the flow is not distorted. The stream line 1608 extends downstream to infinity, forming a semi-infinite body, referred to as the Rankine half-body 1610.

[0228] Thus, it can be appreciated that the flow over the semi-infinite body can be determined from the combination of the uniform flow 1600 and the source flow 1602. Some embodiments are based on the realization that a model of the fluid flow around a cylinder can be obtained by the combination of a uniform flow and a dipole flow.

[0229] Figure 17A A uniform flow 1700 and a dipole flow 1702 are shown in combination for determining the fluid flow around a cylinder, according to some embodiments. According to one embodiment, the dipole flow 1702 is obtained by superimposing a source flow and a sink flow of equal strength.

[0230] Figure 17B A source flow 1708 and a sink flow 1710 of equal strength Λ for obtaining the dipole flow 1702 are shown, according to one embodiment. The source flow 1708 emanates from a point 1712, and the sink flow converges toward a point 1714. The source flow 1708 and the sink flow 1710 are separated by a distance of 2d 1716. As the distance between the source-sink pair (1708 and 1710) approaches zero, the dipole flow 1702 is formed. The velocity potential and the stream function of the dipole flow 1702 are given by the following equations

[0231]

[0232]

[0233] The dipole flow 1702 is superimposed with the uniform flow 1700 to obtain a combined flow 1704. The stream function of the combined flow 1704 is given by the following equation:

[0234]

[0235] In addition, the velocity field is obtained from the following equation

[0236]

[0237] In some embodiments, the velocity components in equation (22) are assigned to zero and r and θ are solved simultaneously to locate stagnation points. In the combined flow 1704, the stagnation points are located at (r, θ) = (R, 0) and (R, π), denoted by points A and B, respectively. The streamlines equations passing through the stagnation points A and B are given as:

[0238]

[0239] Equation (23) is satisfied by r = R for all values of θ. Since R is a constant, equation 23 can be interpreted as the equation of a circle centered at the origin with radius R. For all values of R, θ = 0 and π are satisfied. For this reason, the horizontal axis 1718 passing through points A and B, extending infinitely far upstream and downstream, is part of the stagnation streamline.

[0240] In the combined flow 1704, the dividing streamlines are the circle 1706 of radius R. Different values of R can be obtained by varying the uniform velocity and / or the dipole strength. The flow inside the circle 1706 is generated by the dipole flow 1702, while the flow outside the circle 1706 comes from the uniform flow 1700. Therefore, the flow inside the circle can be replaced by a solid (cylinder) without distorting the flow outside the circle 1706. Thus, the fluid flow over a cylinder of radius R can be emulated by adding a uniform flow 1700 with velocity V ∞ and a dipole flow 1702 with strength Λ, and the relationship of R with V ∞ and Λ is

[0241]

[0242] In addition, multiple elementary potential flows (such as uniform flow, source / sink flow, dipole flow, etc.) can be used to approximate fluid flow over complex shapes.

[0243] Some embodiments are based on the goal of determining a mapping between a cylinder and a complex shape (e.g., a complex terrain). In some embodiments, such a mapping can be determined by using a conformal mapping that includes an analytic mapping of complex numbers. In a conformal mapping, a transformation function is used to transform a function of complex numbers from one coordinate system to another coordinate system. In other embodiments, the mapping between a cylinder and a complex shape is determined based on a machine learning approach.

[0244] For example, one technique involves "training" a machine learning program on various shapes and complex terrain representative of typical sites with converged CFD data. On the order of hundreds of such simulations can be required for the training. Once the program is trained, a process (e.g., using Gaussian process regression or deep learning techniques) is used to infer the velocity and pressure as well as the horizontal gradient of the vertical velocity for a new complex terrain shape based on all the previous complex terrain shapes. In the next step, for each shape, the equivalent radii of a single or multiple cylinders can be determined using, for example, the DAL method to solve the inverse problem. The solution of this optimization problem can be used for training. Once a new complex terrain is encountered, the trained regression equation can be used to determine the equivalent radii.

[0245] Figure 18 An exemplary mapping between a cylinder of radius b and terrain is shown, according to some embodiments. One embodiment of such a conformal mapping is the Joukowski airfoil, which refers to the solution of potential flow through a series of airfoil shapes. The process involves finding a mapping that transforms a cylinder into an airfoil shape. Such a mapping is called a conformal mapping. In mathematics, a conformal mapping is a function that locally preserves angles but not necessarily lengths. A transformation is conformal as long as the Jacobian determinant at each point is a positive scalar times a rotation matrix (orthogonal to 1 determinant). Some embodiments restrict conformal to include orientation reversing mappings whose Jacobian determinant can be written as any scalar times any orthogonal matrix.

[0246] Some embodiments are based on the recognition that a set of convex shapes (cylinders) can be superimposed for mapping or approximating complex terrain.

[0247] Figure 19 A superposition of a set of cylinders 1900, 1902, 1904, and 1906 for mapping with terrain 1908 is shown, according to some embodiments. The radii of the cylinders 1900, 1902, 1904, and 1906 are r1, r2, r3, and r4, respectively. The above cylinders are superimposed to form a superimposed shape 1912 or a distribution of cylinders. The radii of the above cylinders are determined such that the streamwise velocity of Laplace flow of these cylinders is closest to the measured value.

[0248] The radii r1, r2, r3, and r4 of the cylinders are unknown. In addition, the velocity U of the wind 1910 is unknown. In some embodiments, the velocity U can correspond to an upstream velocity. Some embodiments are based on the recognition that the radii of the cylinders N is the number of cylinders, and the velocity U can be determined by the adjoint method. Refer to Figure 21 The determination of the cylinder radii and the velocity U is explained in detail.

[0249] Figure 20A schematic diagram showing the construction and evaluation of a cost function according to some embodiments, which includes both the LOS measurements and the simulated LOS according to Laplace superposition. With the line-of-sight rates according to LiDAR measurements 2000 and according to Laplace superposition 2002, the cost function evaluation yields a numerical value that represents how well the line-of-sight rates according to Laplace superposition 2002 match the line-of-sight rates according to LiDAR measurements 2000 along the line of sight of different beams at various altitudes. A weighting factor 2004 is chosen such that each altitude correction term is proportional to the amount of deviation. One embodiment is to use the altitude above the LiDAR as such a weighting, since the deviation becomes higher at higher altitudes. The cost function is evaluated using a multi-cylinder light-weight method that includes an analytical solution. As a result, the computation time is greatly reduced since there is no numerical solution to the partial differential equation (PDE) or difference equation, which then allows for online reconstruction of the wind or online estimation of the horizontal velocity at each altitude.

[0250] Figure 21 A block diagram showing the implementation of DAL to determine the cylinder radii and upstream velocity U according to some embodiments is shown. This implementation estimates the most likely distribution of cylinders and the upstream velocity by minimizing a cost function. The DAL method is initialized 2100, where the initial estimates of the cylinder radii and the upstream velocity are made. Here, DAL is an optimization method that includes obtaining analytical solutions of the potential flow 2102 and the adjoint (or sensitivity) equation 2104 of the cylinder distribution in an iterative 2114 manner. This optimization provides the sensitivity 2106 of the cost function with respect to the unknowns (i.e., the radii of the cylinders) and the upstream velocity, which is the current estimate of the unknowns. After each iteration, the estimates 2108 of the current values of the cylinder radii and the upstream velocity are updated using a conjugate gradient descent. According to some embodiments, updating the current estimates of the unknowns using a conjugate gradient descent involves an update in the direction of the maximum reduction of the sensitivity of the cost function.

[0251] In addition, a convergence criterion 2110 is checked. One embodiment of such a convergence criterion is the variation of the cost function between successive iterations. Another embodiment is that the change in the estimates according to the previous iteration is below a threshold. If the convergence criterion is not met, the next iteration is initiated, where the analytical solution of the potential flow of the cylinder distribution is determined 2102. In the case where the convergence criterion is met, the final estimates of the cylinder radii and the upstream velocity are obtained 2112.

[0252] Figure 22A schematic diagram showing estimation of horizontal gradient of vertical velocity according to some embodiments is shown. Line-of-sight velocities 2200 are obtained from LiDAR measurements. Furthermore, a most likely distribution of cylinders and upstream velocities is determined 2202 using the line-of-sight velocities obtained from the LiDAR measurements. In some embodiments, the most likely distribution of cylinders and upstream velocities is determined by minimizing a cost function 2202. For example, the minimization of the cost function can be performed iteratively based on the sensitivity of the cost function. Based on the determined distribution of cylinders and upstream velocities, the horizontal gradient of vertical velocity can be estimated 2204.

[0253] Turbulence of the wind flow

[0254] Figure 23A and Figure 23B Collectively, some embodiments show a schematic overview of the principle of some embodiments for turbulence measurement of the wind flow in complex terrain. Remote sensing instruments, such as LiDARs, are used to measure the radial velocity of the wind in line-of-sight (LOS) points for each elevation at a set of time steps 2300. However, the amount of turbulence related to the horizontal velocity is the parameter of interest.

[0255] To this end, some embodiments aim at determining the amount of turbulence related to the horizontal velocity of the wind flow for each elevation. Some embodiments are based on the recognition that, with the geometric relations and assuming that the wind velocity is uniform on each plane inside the cone of measurements, an estimate of the horizontal velocity at a time step can be determined from the measurements of radial velocities 2300 corresponding to the time step. Figure 23B The signal 2312 in FIG. 23 represents an exemplary plot of such horizontal velocity estimates at different instants in time. Furthermore, the horizontal velocity is averaged over a certain time period, e.g., 10 minutes. Figure 23B The horizontal line 2314 shown in FIG. 23 represents the average of the horizontal velocity over this time period. The horizontal velocity 2312 and the average of the horizontal velocity 2314 are different at each instant in time. The square of the difference between the horizontal velocity 2312 and the average 2314 is called the standard deviation. According to some embodiments, the standard deviation of the average of the horizontal velocity defines the turbulence intensity (TI). The turbulence intensity is defined over a certain time period, e.g., 10 minutes. To this end, some embodiments are based on the recognition that the turbulence intensity is a function of the shape and the corresponding values of the signal 2312, and thus the turbulence intensity depends on the instantaneous values of the horizontal velocity.

[0256] However, for wind flow over complex terrain, such as near hills or large buildings or other urban structures, the uniform velocity assumption considered for the horizontal velocity estimation is not valid. Some embodiments are based on the recognition that for complex terrain, the uniform velocity assumption leads to a bias in the horizontal velocity estimation. This bias in the horizontal velocity estimation leads to a biased horizontal velocity, which introduces a bias into the standard deviation. This bias is due to the variation of the vertical velocity in the vertical direction. According to one embodiment, the bias is given by the height times the horizontal gradient of the vertical velocity, i.e.

[0257] Some embodiments are based on the recognition that the horizontal derivative of the vertical velocity can be used as a correction to the biased horizontal velocity to remove the bias 2302. The horizontal derivative of the vertical velocity is determined by estimating the velocity field. The velocity field at each time is estimated 2304 based on data assimilation that is used to fit the measurements of the radial velocity 2300. According to some embodiments, the data assimilation is performed by computational fluid dynamics (CFD) that simulates the wind flow. Reference is made to Figures 5A to 13 The data assimilation is explained in more detail. According to other embodiments, the data assimilation is performed by using an analytical fluid mechanics approximation that uses a potential flow approximation, which is explained with reference to Figures 14 to 22 in more detail.

[0258] This bias removal is performed for the horizontal velocity at each time 2312 to obtain the unbiased horizontal velocity at the respective time 2316. In other words, instantaneous bias removal is performed to obtain the unbiased horizontal velocity at each time. In addition, the average of the unbiased horizontal velocities 2316 for the time period is determined 2306. The average of the unbiased horizontal velocities is represented by the horizontal line 2318 in Figure 23B The standard deviation is determined 2308 as the square of the difference between the unbiased horizontal velocities 2316 and the average of the unbiased horizontal velocities 2318. Since the standard deviation is determined 2308 based on the unbiased horizontal velocities, the bias that was present in the standard deviation due to the biased horizontal velocities is removed. This then improves the accuracy of the standard deviation.

[0259] According to some embodiments, the turbulent quantities (e.g., the turbulence intensity (TI) and the turbulent kinetic energy (TKE)) are determined based on the unbiased horizontal velocities at each time step and the average of the unbiased horizontal velocities.

[0260] The component u of the wind velocity can be given as:

[0261]

[0262] where, is the average velocity for the time period T, i.e. In one embodiment, T = 10 minutes, which results in an average rate of 10 minutes. u' is the fluctuating rate. Here, u is the instantaneous rate. The root mean square value is:

[0263]

[0264] For other rate components (e.g., v) and / or horizontal rate v h Equation 23 can be similarly written.

[0265] The turbulence intensity (TI) is determined 2310 from the ratio of the root mean square value of the fluctuating turbulent rate to the average unbiased horizontal rate. TI is given by the following equation

[0266]

[0267] The turbulent kinetic energy (TKE) is the standard deviation of the components of the wind rate and is given by the following equation:

[0268] TKE = u rms + v rms

[0269] Since the turbulent mass is determined using the unbiased horizontal rate obtained by removing the instantaneous bias and the average of the unbiased horizontal rate, rather than using the biased horizontal rate, the accuracy of the turbulent mass is significantly improved.

[0270] Some embodiments are based on the recognition that by formulating a function (i.e., an autocorrelation function), a relationship between the standard deviations of different points on a scan circle at a certain height can be established. The autocorrelation function is also referred to as a correction function. Thus, the autocorrelation function is related to the standard deviations of the points and can be used to measure the standard deviation of the estimated horizontal rate to produce turbulence. Some embodiments are based on the recognition that the correction function can be used to correct the unbiased horizontal rate before estimating the turbulence.

[0271] Figure 24A A schematic diagram showing the correction of the standard deviation using the autocorrelation function is shown. Some embodiments are based on the recognition that the correction function can be used to correct the unbiased horizontal rate before estimating the turbulence. The correction function is trained offline (i.e., in advance). According to one embodiment, the correction function is applied to the determined turbulence to obtain the actual turbulence. The terms "correction function" and "autocorrelation function" can be used interchangeably and have the same meaning.

[0272] According to one embodiment, the actual variances of u and v can be given as:

[0273]

[0274]

[0275] where pu , p v , p w is an autocorrelation function (ACF), and and are obtained from LiDAR measurements. Equations (25) and (26) can be referred to as correction equations.

[0276] Some embodiments are based on the recognition that autocorrelation functions are machine learned based on physics knowledge and can be used to correct the standard deviation of the horizontal velocity with anemometer measurements. Some embodiments are based on the recognition that anemometer (e.g., sonic anemometer) and LiDAR measurements can be used to train autocorrelation functions. and are determined by anemometer (ground truth), and in addition, and determined values are compared to LiDAR measurements to estimate the ACF.

[0277] Figure 24B A plot showing the values of autocorrelation functions p u , p v , and p w calculated based on comparison with cup anemometer data according to some embodiments. Sonic anemometers are used to simulate the measurement technique used by LiDAR 2406. For example, at each measurement height, two sonic anemometers are placed on opposite booms about 11.5 meters apart. The sonic data is projected to the direction of the LiDAR beam position. The projected data from the southern anemometer is time advanced by 2 seconds to simulate the time it takes for the LiDAR beam to move from one side of the scan circle to the other. Based on the time shifted and projected data (i.e., sonic data 2408), the values of autocorrelation functions p u , p v , and p w are calculated 2410.

[0278] According to one embodiment, the average values of p u , p v , and p w are calculated from the sonic data. For example, in one case, under unstable conditions, the average values of p u , p v , and p w may correspond to 0.96, 0.81, and 0.66, respectively. Under stable conditions, the average values of p u , p v , and p wThe average values ​​can be found at 0.95, 0.71, and 0.69. These values ​​indicate that the u, v, and w wind components vary significantly in both space and time. In particular, the value of w suggests that, due to the smaller scale of turbulent motion in the vertical direction, the value of w decorrelates more quickly than the values ​​of u and v.

[0279] In addition, ρ was calculated based on acoustic data. u ρ v , and ρ w The average value is used in conjunction with equations (25) and (26) to correct the variance of LiDAR, where the correction equation is in The value is considered to be the rate variance measured by the LiDAR vertical beam. Under steady-state conditions, when When the value is small, variance correction will not significantly change the variance value, but under unstable conditions, it will... and The estimated value was reduced by more than 20%, thus making the estimated value closer to the value measured by the anemometer.

[0280] Other implementations are based on the understanding that the value of the autocorrelation function can be calculated using the least squares method. Such implementations produce values ​​ρ similar to those calculated from acoustic data. u and ρ v The value, and produces a ρ much lower than the acoustic value. w The value of .

[0281] According to some implementations, the ACF can be trained taking into account the actual terrain. Furthermore, in another implementation, the ACF is trained for terrain approximated by a set of convex shapes. In yet another implementation, the ACF is trained taking into account the actual terrain, and the ACF can be applied during online estimation using terrain approximated by a set of convex shapes to obtain the actual turbulence.

[0282] Figure 24C This illustrates how, according to some implementation methods, the autocorrelation function ρ is calculated based on a comparison with high-fidelity CFD simulations. u ρ v , and ρ wof the values of p u , p v , and p w . In other words, a least squares problem is solved to evaluate the autocorrelation function values.

[0283] Alternatively, such values can be determined by mapping anemometer data or high-fidelity simulation data to the standard deviation of the horizontal velocity extracted from the LOS data and fitting the autocorrelation functions in equations 25 and 26.

[0284] Figure 25 A block diagram of a wind flow sensing system 2500 for determining the turbulence of a wind flow according to some embodiments is shown. The wind flow sensing system 2500 includes an input interface 2502 to receive a set of measurement values 2518 of radial velocity in the site line direction for an elevation for a set of time steps. In some embodiments, the measurement values 2518 on the cone are measured by a remote sensing instrument, such as a ground-based LiDAR. The wind flow sensing system 2500 can have a number of interfaces that connect the system 2500 with other systems and devices. For example, a network interface controller (NIC) 2514 is adapted to connect the wind flow sensing system 2500 to a network 2516 via a bus 2512 that connects the wind flow sensing system 2500 with a remote sensing instrument configured to measure the radial velocity of the wind flow for each time step. With the network 2516, whether wired or wireless, the wind flow sensing system 2500 receives a set of measurement values 2518 of radial velocity at the site line direction for each elevation for a set of time steps.

[0285] Further, in some implementations, measurements 2518 can be downloaded and stored within storage system 2536 by way of network 2516 for further processing. Additionally or alternatively, in some implementations, wind flow sensing system 2500 includes a human machine interface 2530 that connects processor 2504 to a keyboard 2532 and a pointing device 2534, which can include a mouse, trackball, touchpad, joystick, pointing stick, stylus, touch screen, or the like.

[0286] Wind flow sensing system 2500 includes a processor 2504 configured to execute stored instructions and a memory 2506 that stores instructions executable by the processor. Processor 2504 can be a single core processor, multi-core processor, computing cluster, or any number of other configurations. Memory 2506 can include random access memory (RAM), read only memory (ROM), flash memory, or any other suitable memory systems. Processor 2504 is connected to one or more input and output interfaces and / or devices by way of bus 2512.

[0287] According to some implementations, the instructions stored in memory 2506 implement a method for determining, for a set of time steps, a velocity field of a wind flow at a set of different altitudes from a set of measurements 2518 of radial velocities at each altitude of the set of different altitudes. To this end, storage 2536 can be adapted to store different modules that store executable instructions for processor 2504. Storage 2536 stores CFD simulation module 208, CFD operating parameters module 210, horizontal derivative module 236, velocity field module 238, and Laplace simulation module 240. Figure 2 These modules are explained in the description of FIG. 2. Additionally, storage 2536 stores a correction function 2538 that is trained to reduce a difference between ground truth and determined unbiased horizontal velocities to correct the unbiased horizontal velocities. Storage 2536 can be implemented using a hard disk, optical drive, thumb drive, array of drives, or any combination thereof.

[0288] Processor 2504 is configured to estimate, for each time step and each altitude, an unbiased horizontal velocity as a horizontal projection of a respective radial velocity corrected with a respective horizontal derivative of a vertical velocity of an estimated velocity field determined for the respective time step. Additionally, processor 2504 is configured to determine, at each altitude, an average of the unbiased horizontal velocities for a time period that includes the set of time steps, and then determine a turbulence based on the unbiased horizontal velocities for each time step and the average of the unbiased horizontal velocities. According to one implementation, the time period is a multiple of 10 minutes and the difference between time steps is a multiple of one second.

[0289] The wind flow sensing system 2500 includes an output interface 2524 for presenting the turbulence at each altitude. In addition, the wind flow sensing system 2500 can be linked by the bus 2512 to a display interface 2520 adapted to connect the wind flow sensing system 2500 to a display device 2522, which can be a computer display, a camera, a television, a projector, or a mobile device, among others. In addition, the wind flow sensing system 2500 includes a control interface 2526 configured to submit the estimated wind flow at each altitude to a controller 2528 integrated with a machine, such as a wind turbine. According to one embodiment, the controller 2528 is configured to operate the machine based on the estimated wind flow at each altitude.

[0290] Figure 26 A schematic view of a wind turbine 2602 is shown, which includes a controller 2606 in communication with the system 2500 that employs principles of some embodiments. The wind turbine 2602 on the complex terrain 2604 is integrated with the system 2500. The wind turbine 2602 can be equipped with a transceiver, giving the controller 2606 of the wind turbine the ability to communicate with the system 2500 by way of a wired or wireless communication channel. For example, by way of the transceiver, the controller 2606 receives estimates of wind parameters from the system 2500.

[0291] The LiDAR 300 measures the LOS velocity of the wind 2600 flowing through the complex terrain 2604 and the blades of the wind turbine 2602 for each altitude for a set of time steps. The system 2500 receives the LiDAR measurements (as described in Figure 25 Based on the received LiDAR measurements, the system 2500 estimates the unbiased horizontal velocity of the wind flow 2600 at each altitude for each time step. In addition, the system 2500 estimates the turbulent quantities (such as the turbulence intensity and the turbulent kinetic energy) of the wind flow 2600. The system 2500 submits the estimated horizontal velocity and the turbulent quantities to the controller 2606. The controller 2606 generates control inputs based on the estimated horizontal velocity and the turbulent quantities for controlling the wind turbine 2602.

[0292] For example, a horizontal velocity greater than a threshold and turbulence creates irregular wind loads on the wind turbine 2602. Operating the wind turbine 2602 in such conditions affects both energy production and the structure of the blades of the wind turbine 2602. In the case of an estimated horizontal velocity or turbulence greater than a threshold, the controller 2606 generates control inputs that cause the wind turbine 2602 to stop or brake. In this way, the adverse effects of irregular wind on the wind turbine 2602 are prevented. Further, the controller 2606 can actuate the wind turbine 2602 based on the estimated horizontal velocity, rather than arbitrarily actuating the wind turbine. Alternatively, the controller 2606 can actuate the wind turbine 2602 according to both the estimated horizontal velocity and turbulence. Thereby, energy production of the wind turbine 2602 is improved and the wind loads experienced by the wind turbine 2602 are reduced, which in turn prolongs the life of the wind turbine 2602.

[0293] The following description provides example implementations only, and is not intended to limit the scope, applicability or configuration of the disclosure. Rather, the following description of the example implementations will provide those skilled in the art with enabling information to make and use one or more implementations. Various changes can be made in the function and arrangement of elements without departing from the spirit and scope of the subject matter disclosed in the appended claims.

[0294] Specific details are given in the following description to provide a thorough understanding of implementations. However, implementations can be practiced without these specific details. For example, the systems, processes, and other elements in the disclosed subject matter can be shown as components in block diagram form in order not to obscure the implementations in unnecessary detail. In other instances, well-known processes, structures, and techniques can be shown without detailed

[0295] Also, various implementations can be described as a process that is depicted as a flow diagram, a flowchart, a data flow diagram, a structure diagram, or a block diagram. Although a flow diagram can describe operations as a sequential process, many of the operations can be performed in parallel, or concurrently, or in an order of different than that shown. In addition, the order of the operations can be re-arranged. A process can correspond in part to a function, procedure, subroutine, subroutine, or the like. When a process corresponds to a function, its termination can correspond to a return of the function to the calling function or the main function.

[0296] Furthermore, embodiments of the disclosed subject matter can be implemented, at least in part, either manually or automatically. Manual or automatic implementations can be executed or at least assisted with the use of machines, hardware, software, firmware, middleware, microcode, hardware description languages, or any combination thereof. When implemented in software, firmware, middleware or microcode, the program code or code segments to perform the necessary tasks can be stored in a machine readable medium. Processors can perform the necessary tasks.

[0297] The various methods or processes outlined herein can be coded as software that is executable on one or more processors that employ any one of a variety of operating systems or platforms. Additionally, such software can be written using any of a number of suitable programming languages and / or programming or scripting tools, and also can be compiled as executable machine language code or intermediate code that is executed on a framework or virtual machine. Typically, the functionality of the program modules can be combined or distributed as desired in various embodiments.

[0298] Embodiments of the present disclosure can be embodied as a method, of which an example has been provided. The acts performed as part of the method can be ordered in any suitable way. Accordingly, embodiments can be constructed in which acts are performed in an order different than illustrated, which can include performing some acts simultaneously, even though shown as performed sequentially in illustrative embodiments. Further, the use of ordinal terms such as "first," "second," etc., in the claims to modify a claim element does not in and of itself require that the element needs to be the first or second, etc., in an object, nor does it require any priority as to when that element needs to be performed, but rather merely distinguishes that element from another, as used in the specification.

[0299] While the present disclosure has been described with reference to certain preferred embodiments thereof, a variety of other variations and modifications are possible. Accordingly, the scope of the present disclosure encompasses all such variations and modifications.

Claims

1. A wind flow sensing system, the wind flow sensing system being used to determine the turbulent flow rate of wind at each elevation set based on a set of measured radial velocities at each elevation set along the surface of a cone and along the centerline of the cone, the wind flow sensing system comprising: An input interface configured to receive a set of measurements of radial velocity at a field line-of-sight point above the terrain at each elevation for a set of time steps. Processor, the processor being configured to: Based on data assimilation of the velocity field above the terrain, a velocity field is estimated for each elevation, the data assimilation being used to fit measurements of the radial velocity, wherein the velocity field is estimated for the set of time steps; and, The turbulent flow rate at each altitude is determined based on the unbiased horizontal velocity and the average of the unbiased horizontal velocities at each time step; and An output interface configured to present the turbulent flow rate; The processor is further configured to: The unbiased horizontal velocity at each altitude for each time step is estimated as the horizontal projection of the corresponding radial velocity, and the unbiased horizontal velocity is corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field determined for the corresponding time step. Determine the average value of the unbiased horizontal rate at each altitude for a time period including the set of time steps; The turbulent flow rate at each altitude is determined based on the unbiased horizontal velocity and the average of the unbiased horizontal velocities at each time step; and The output interface is configured to display the turbulent flow rate at each altitude. The data assimilation is performed by using an analytical fluid dynamics approximation based on the potential flow approximation.

2. The airflow sensing system according to claim 1, wherein, The processor is configured to: The non-convex shapes of the terrain are approximated using a set of convex shapes; Determine the boundary conditions for the inlet velocity field; An analytical solution to the Laplace equation constraining the airflow is derived using the aforementioned boundary conditions. and Update the boundary conditions and repeat the simulation until the termination condition is met.

3. The airflow sensing system according to claim 1, wherein, The turbulent flow rate includes the turbulence intensity determined based on the ratio of the root mean square value of the turbulence rate fluctuation to the average unbiased horizontal rate.

4. The airflow sensing system according to claim 3, wherein, The turbulent flow rate includes the turbulent kinetic energy determined by half the sum of the root mean square values ​​of the turbulent rate fluctuations.

5. The airflow sensing system according to claim 1, wherein, The time period is a multiple of 10 minutes, and the difference in time steps is a multiple of 1 second.

6. The airflow sensing system according to claim 1, wherein the airflow sensing system further comprises: A memory configured to store a trained correction function to correct the unbiased horizontal rate, wherein the processor applies the correction function to correct the unbiased horizontal rate before estimating the turbulent flow rate.

7. The airflow sensing system according to claim 6, wherein, The correction function is trained to reduce the difference between the ground reality and the unbiased horizontal velocity, which is determined based on the analytical solution of the Laplace equation for wind flow over approximate terrain.

8. The wind flow sensing system of claim 1, further comprising a controller for a wind turbine to control the wind turbine based on one or more estimated unbiased horizontal velocities or turbulent flow rates at each altitude, wherein, The wind turbine is operatively connected to the airflow sensing system.

9. A wind flow sensing method, wherein the wind flow sensing method is used to determine the turbulent flow rate of wind at each elevation set based on a set of measured radial velocities at each elevation along the surface of a cone and along the centerline of the cone, wherein, The airflow sensing method uses a processor coupled to stored instructions for implementing the airflow sensing method, wherein the instructions, when executed by the processor, implement steps of the airflow sensing method, the steps including: For each time step set, a set of measurements of radial velocity at a field line-of-sight point above the terrain at each elevation level is received. Based on data assimilation of the velocity field above the terrain, a velocity field at each elevation is estimated, the data assimilation being used to fit the measured values ​​of the radial velocity, wherein the velocity field is estimated for the set of time steps; The turbulent flow rate at each altitude is determined based on the unbiased horizontal velocity and the average of the unbiased horizontal velocities at each time step; and Output the turbulent flow rate; The steps of the airflow sensing method further include: The unbiased horizontal velocity at each altitude for each time step is estimated as the horizontal projection of the corresponding radial velocity, and the unbiased horizontal velocity is corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field determined for the corresponding time step. Determine the average value of the unbiased horizontal rate at each altitude for a time period including the set of time steps; The turbulent flow rate at each altitude is determined based on the unbiased horizontal velocity and the average of the unbiased horizontal velocities at each time step; and Turbulent flow rates at each altitude are output; The data assimilation is performed by using an analytical fluid dynamics approximation based on the potential flow approximation.

10. The airflow sensing method according to claim 9, wherein, The turbulent flow rate includes the turbulence intensity determined based on the ratio of the root mean square value of the turbulence rate fluctuation to the average unbiased horizontal rate.

11. The airflow sensing method according to claim 10, wherein, The turbulent flow rate includes the turbulent kinetic energy determined by half the sum of the root mean square values ​​of the turbulent rate fluctuations.

Images