commit 2971b9eaae560c2be1fcff0b529afdf41a0902fd
parent c2a3e6ad754c2038babb31aa4e5f6eaf611e8fe9
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Fri, 9 Apr 2021 20:01:14 +0300
Introduction.
Diffstat:
5 files changed, 235 insertions(+), 43 deletions(-)
diff --git a/main.bib b/main.bib
@@ -92,3 +92,71 @@
url = {https://www.sciencedirect.com/science/article/pii/S0196890408000228},
author = {José A. Carta and Penélope Ramírez and Celia Bueno}
}
+
+@InProceedings{ cosgrove2007ultra,
+ title = {Ultra-low-cost logging anemometer for wind power generation
+ feasibility surveys},
+ author = {Cosgrove, Michael and Rhodes, Bruce and Scott, Jonathan B.},
+ year = {2007},
+ booktitle = {Proceedings of the 14\textsuperscript{th} New Zealand
+ Electronics Conference},
+ url = {https://hdl.handle.net/10289/1485},
+ address = {Wellington, New Zealand},
+ pages = {153--158}
+}
+
+@InProceedings{ yakunin20173d,
+ title = {3D Ultrasonic Anemometer with tetrahedral arrangement of
+ sensors},
+ author = {Yakunin, A. G.},
+ booktitle = {Journal of Physics: Conference Series},
+ volume = {881},
+ number = {1},
+ pages = {12--30},
+ year = {2017},
+ organization = {IOP Publishing}
+}
+
+@Article{ arens2020anemometer,
+ title = {Measuring 3D indoor air velocity via an inexpensive low-power
+ ultrasonic anemometer},
+ journal = {Energy and Buildings},
+ volume = {211},
+ pages = {109805},
+ year = {2020},
+ issn = {0378-7788},
+ doi = {10.1016/j.enbuild.2020.109805},
+ url = {https://www.sciencedirect.com/science/article/pii/S0378778819331044},
+ author = {Edward Arens and Ali Ghahramani and Richard Przybyla and
+ Michael Andersen and Syung Min and Therese Peffer and Paul
+ Raftery and Megan Zhu and Vy Luu and Hui Zhang}
+}
+
+@Article{ kaimal1972spectral,
+ title = {Spectral characteristics of surface-layer turbulence},
+ author = {Kaimal, J.C. and Wyngaard, J.C. and Izumi, Y. and Cot{\'e},
+ O.R.},
+ journal = {Quarterly Journal of the Royal Meteorological Society},
+ volume = {98},
+ number = {417},
+ pages = {563--589},
+ year = {1972},
+ publisher = {Wiley Online Library}
+}
+
+@InProceedings{ veers1984modeling,
+ title = {Modeling stochastic wind loads on vertical axis wind
+ turbines},
+ author = {Veers, Paul},
+ booktitle = {25\textsuperscript{th} Structures, Structural Dynamics and
+ Materials Conference},
+ pages = {910},
+ year = {1984}
+}
+
+@TechReport{ veers1988wind,
+ title = {Three-dimensional wind simulation},
+ author = {Veers, Paul S},
+ year = {1988},
+ institution = {Sandia National Labs., Albuquerque, NM (USA)}
+}
diff --git a/main.tex b/main.tex
@@ -49,7 +49,106 @@
\section{Introduction}
-\cite{yim2000}
+Wind simulation in the context of ship motion simulation is the topic where
+multiple mathematical models are possible, and the choice of the model depends
+on the purpose of the model. In the course of the research where we apply
+autoregressive model to ocean wave simulation, we decided to investigate
+whether the same model can be used to simulat wind flow round ship hull.
+
+Wind simulation is studied at different scales and the closest scale for a ship
+in the ocean is wind turbine. One of the model that is used to
+describe air flow around wind turbines~\cite{veers1984modeling} is similar to
+Longuet---Higgins model which is typically used for ocean wave simulations.
+This model uses wind velocity frequency spectrum to determine coefficients
+and use them to generate wind velocity vector components.
+
+Wind velocity distribution is described by Kaimal
+spectrum~\cite{kaimal1972spectral,veers1988wind}
+\begin{equation*}
+S(f) = \frac {c_1 u_{*}^2 z / U(z)} {1 + c_2 \left(f z / U(z)\right)^{5/3} },\quad
+u_{*}=\frac{k U(z)}{\ln{(z/z_0)}},\quad c_1=105, \quad c_2=33,
+\end{equation*}
+Here \(u_{*}\) is shear velocity,
+\(k=0.4\) is von Karman constant, \(z_0\) is surface roughness
+coefficient, \(f\) is frequency,
+\(z\) is height above ground, \(U(z)\) is mean speed at height \(z\).
+
+This spectrum is used to simulate each wind velocity vector component at a
+specified point in space. For each component the same spectrum is used, but the
+coefficients are different~\cite{veers1984modeling}.
+%\begin{align*}
+%& c_1 = 12.3 && c_2 = 192 && \text{for }x,\\
+%& c_1 = 4 && c_2 = 70 && \text{for }y,\\
+%& c_1 = 0.5 && c_2 = 8 && \text{for }z.
+%\end{align*}
+The spectrum describes wind velocity vector in the plane that is perpendicular
+the mean wind direction vector and travels in the same direction
+with mean wind speed. Time series is generated as Fourier series coefficients
+of which are determined from the spectrum and phases are random variables:
+\begin{align*}
+& V(t) = \overline{V} + \sum\limits_{j=1}^{n}
+ \left( A_j \sin\omega_jt + B_j \cos\omega_jt\right), \\
+& A_j = \sqrt{\frac{1}{2} S_j \Delta\omega} \sin\phi_j, \\
+& B_j = \sqrt{\frac{1}{2} S_j \Delta\omega} \cos\phi_j. \\
+\end{align*}
+Here \(S_j\) is spectrum value at frequency \(\omega_j\),
+\(\phi_j\) is random variable which is uniformly distributed in \([0,2\pi]\).
+The result is one-dimensional vector-valued time series, each element of
+which is velocity vector at a specified point in time and space.
+
+In order to simulate wind velocity vector in multiple points in space,
+the authors use the function of coherency~--- the amount of correlation between
+wind speed at two points in space. This function has a form of an exponent and
+depends on frequency~\cite{veers1988wind}:
+\begin{equation*}
+\text{Coh}_{jk}(f) = \exp\left( -\frac {C\Delta r_{jk} f} {U(z)} \right),
+\end{equation*}
+where \(\Delta{}r_{jk}\) is the distance betweein \(i\) and \(j\) points and
+\(C\) is coherency decrement. Time series for each wind velocity vector
+component are generated independenty and after that their spectra are modified
+in accordance with coherence function (the formulae are not presented here).
+
+In order to simulate wind with autoregressive model it is easier to use
+autocovariance function instead of spectra and coherence function. We can obtain
+autocovariance function from spectra as inverse Fourier transform using
+Wiener---Khinchin theorem. The formula that we obtained this way using
+computer algebra programmes is to complex, but can be approximated by
+a decaying exponent:
+\begin{equation*}
+\gamma(t) = \sigma^2 \exp\left( -0.1 \frac{c_2^{3/5}}{c_1} t\right),
+\end{equation*}
+where \(\sigma^2\) is process variance (area under the spectrum).
+
+Unfortunately, this autocovariance function is one-dimensional and there is
+easy way of obtaining three-dimensional analogue from the spectrums. In order
+to solve this problem we looked into datasets of wind speed measurements
+available for the research. However, most of them contain either one or two
+wind velocity vector components (wind speed and direction), they are difficult
+to get and their discretization is very coarse (we found only one paper that
+deals with three components~\cite{yim2000}). For our purposes we need
+descretization of at least one second to simulate gusts and some form of
+turbulence. To summarise, our requirements for the dataset is to provide all
+three wind velocity vector components and has discretization of at least one
+second.
+
+We failed to find such datasets and continued our research into anemometers
+that can generate the required dataset. One of the anemometers that can record
+all three components is ultrasonic
+anemometer~\cite{cosgrove2007ultra,arens2020anemometer,yakunin20173d}. As
+commercial anemometers are too expensive for this research, we built our own
+version from the generally available electric components. However, this version
+failed to capture any meaningful data, and incidentally we decided to build an
+anemometer from load cells and strain gauges. This anemometer is
+straightforward to construct, the electrical components are inexpensive, and it
+is easy to protect them from bad weather. This anemometer is able to record
+all three wind velocity vector components multiple times per second.
+
+In this paper we describe how load cell anemometer is built, then we collect
+dataset of all three wind speed components with one-second resolution, verify
+this anemometer using commercial analogue, verify measurements from the dataset
+using well-established distribution for wind speed and direction, and finally
+estimate autocovariance functions for autoregressive model from our dataset.
+Finally, we present preliminary wind simulation results using this model.
\section{Methods}
@@ -62,19 +161,21 @@ theorem it is easy to compute the function from the spectrum. Unfortunately,
most of the existing wind velocity historical data contains only wind velocity
magnitude and direction. We can use them to reconstruct \(x\) and \(y\)
spectrum, but there is no way to get spectrum for \(z\) coordinate from this
-data. To solve this problem, we decided to build our own three-axis anemometer
-that measures wind velocity for all three axes at one point in space.
+data. Also, data discretisation in historical data is too coarse for our
+purposes (wind simulation for ship motions). To solve these problems, we decided
+to build our own three-axis anemometer that measures wind velocity for all
+three axes at one point in space.
To measure wind velocity we used resistive foil strain gauges mounted on the
arms aligned perpendicular to the axes directions. Inside each arm we placed
aluminum load cell with two strain gauges: one on the bending side and another
-one on the lateral side. Load cell uses Wheatstone bridge to measure the
-resistance of the gauges and compensate for temperature changes. Load cells are
-connected to the circuit that measures the resistance and transmits it to the
-microcontroller in digital format. Microcontroller then records the value for
-each load cell and transmits all of them in textual form to the main computer.
-The main computer then adds a timestamp and saves it to the database.
-3-D anemometer model is show in figure~\ref{fig-anemometer}.
+one on the lateral side. Load cells use Wheatstone bridge to measure the
+resistance of the gauges and are connected to the circuit that measures the
+resistance and transmits it to the microcontroller in digital format.
+Microcontroller then records the value for each load cell and transmits all of
+them in textual form to the main computer. The main computer then adds a
+timestamp and saves it to the database. 3-D anemometer model is shown in
+figure~\ref{fig-anemometer}.
\begin{figure}
\centering
@@ -104,11 +205,11 @@ cell is written as
\vec{F} = p S \vec{n},
\end{equation}
where \(S\) is area of the side of the load cell on which the force is applied
-and \(\vec{n}\) is normal vector. Arms in which the load cell reside have only
-front and back sides with much larger areas than the left and right side, therefore
+and \(\vec{n}\) is normal vector. Arms in which load cells reside have
+front and back side with much larger areas than the left and right side, therefore
we neglect forces on them. Additionally, on the ground \(\rho{}gz\) term vanishes.
-All in all, for all load cells we have a system of three equations
+For all load cells we have a system of three equations
\begin{equation}
\begin{cases}
\upsilon_x^2 \propto F_x \\
@@ -119,7 +220,7 @@ All in all, for all load cells we have a system of three equations
Hence \(\upsilon_{x}=\alpha_{x}\sqrt{F_{x}}\),
\(\upsilon_{y}=\alpha_{y}\sqrt{F_{y}}\),
\(\upsilon_{z}=\alpha_{z}\sqrt{F_{z}}\) where \(\alpha_{x,y,z}\) are constants of
-proportionality. Therefore, to obtain wind velocity to take square root of the
+proportionality. Therefore, to obtain wind velocity we take square root of the
value measured by the load cell and multiply it by some coefficient determined
empirically during anemometer calibration.
@@ -135,33 +236,49 @@ Here \(\upsilon>0\) is scalar wind velocity, \(b>0\) is scale parameter and \(c>
is shape parameter. This function is defined for positive wind velocity, since
scalar wind velocity is a length of wind velocity vector
\(\upsilon=\sqrt{\smash[b]{\upsilon_x^2+\upsilon_y^2+\upsilon_z^2}}\). However, projection of
-wind velocity vector on \(x\), \(y\) or \(z\) axis may be negative. We could simply
-extend Weibull distribution for negative values, but this may lead to multimodal distribution
-because the maximum value may not be located at \(\upsilon=0\). Our solution to
-this problem is to make substitution \(\upsilon_x'=\left|\upsilon_x-\Mode\upsilon_x\right|\),
-where \(\Mode\upsilon_x\) is the mode of wind velocity projection to the corresponding axis,
-and modulus makes the graph simmetric with respect to \(\upsilon_x'=0\).
-The substitution yields
+wind velocity vector on \(x\), \(y\) or \(z\) axis may be negative.
+Our solution to this problem is to use two Weibull distributions: one both positive and
+one for negative projections --- but with different parameters.
\begin{equation}
- f\left(\upsilon_x';b,c\right) = b c
- \exp\left(-\left|b\upsilon_x'-g\left(c\right)\right|^{c}\right),
- \qquad
- g\left(c\right) = \left(\frac{c-1}{c}\right)^{1/c},
-\end{equation}
-which is the formula for the right side of the distribution. In order to get the
-formula for the left side, we have to use an additional set of parameters \(b\), \(c\).
-Then the final formula for both sides is written as
-\begin{equation}
- f\left(\upsilon_x';b_1,c_1,b_2,c_2\right) = \begin{cases}
- b_1 c_1 \exp\left(-\left|b_1\upsilon_x'-g\left(c_1\right)\right|^{c_1}\right)
- \qquad \text{if } \upsilon_x < \Mode\upsilon_x \\
- b_2 c_2 \exp\left(-\left|b_2\upsilon_x'-g\left(c_2\right)\right|^{c_2}\right)
- \qquad \text{if } \upsilon_x \geq \Mode\upsilon_x \\
+ f\left(\upsilon_x;b_1,c_1,b_2,c_2\right) = \begin{cases}
+ b_1 c_1 \left(b_1 \upsilon_x\right)^{c_1-1} \exp\left(-\left(b_1\upsilon_x\right)^{c_1}\right).
+ \qquad \text{if } \upsilon_x < 0 \\
+ b_2 c_2 \left(b_2 \upsilon_x\right)^{c_2-1} \exp\left(-\left(b_2\upsilon_x\right)^{c_2}\right).
+ \qquad \text{if } \upsilon_x \geq 0 \\
\end{cases}\label{eq-velocity-distribution}.
\end{equation}
Here \(b_{1,2}>0\) and \(c_{1,2}>0\) are parameters of the distribution
-that control the scale and the shape.
+that control the scale and the shape, \(\upsilon_x\) is the projection of the velocity
+vector on \(x\) axis. The same formula is used for \(y\) and \(z\) axis.
+
+%We could simply
+%extend Weibull distribution for negative values, but this may lead to multimodal distribution
+%because the maximum value may not be located at \(\upsilon=0\). Our solution to
+%this problem is to make substitution \(\upsilon_x'=\left|\upsilon_x-\Mode\upsilon_x\right|\),
+%where \(\Mode\upsilon_x\) is the mode of wind velocity projection to the corresponding axis,
+%and modulus makes the graph simmetric with respect to \(\upsilon_x'=0\).
+%The substitution yields
+%\begin{equation}
+% f\left(\upsilon_x';b,c\right) = b c
+% \exp\left(-\left|b\upsilon_x'-g\left(c\right)\right|^{c}\right),
+% \qquad
+% g\left(c\right) = \left(\frac{c-1}{c}\right)^{1/c},
+%\end{equation}
+%which is the formula for the right side of the distribution. In order to get the
+%formula for the left side, we have to use an additional set of parameters \(b\), \(c\).
+%Then the final formula for both sides is written as
+%\begin{equation}
+% f\left(\upsilon_x';b_1,c_1,b_2,c_2\right) = \begin{cases}
+% b_1 c_1 \exp\left(-\left|b_1\upsilon_x'-g\left(c_1\right)\right|^{c_1}\right)
+% \qquad \text{if } \upsilon_x < \Mode\upsilon_x \\
+% b_2 c_2 \exp\left(-\left|b_2\upsilon_x'-g\left(c_2\right)\right|^{c_2}\right)
+% \qquad \text{if } \upsilon_x \geq \Mode\upsilon_x \\
+%
+% \end{cases}\label{eq-velocity-distribution}.
+%\end{equation}
+%Here \(b_{1,2}>0\) and \(c_{1,2}>0\) are parameters of the distribution
+%that control the scale and the shape.
\subsection{Three-dimensional ACF of wind velocity}
@@ -264,10 +381,11 @@ figure~\ref{fig-intervals}.
Unique feature of three-axis anemometer is that it measures both velocity of
incident air flow towards the arms and the turbulent flow that forms behind the
-arms. Turbulence is complicated phenomena and there is no simple way of
-excluding it from the statistics, therefore, it is important to choose
-intervals with unimodal distributions (with one mean wind direction) to get
-statistics that is not distorted by the turbulence.
+arms. Turbulent flow velocity distribution mode is often smaller than the mode
+of incident flow. To exclude it from the measurements one can choose the side
+with the largest mode (either positive or negative part
+of~\eqref{eq-velocity-distribution}), but in this paper we left them as is
+for the purpose of the research.
\begin{figure}
\centering
@@ -356,7 +474,7 @@ distributions for each axis are presented in figure~\ref{fig-velocity-distributi
Wind direction was approximated by von Mises distribution using least-squares
fitting. Following the common
practice~\cite{feng2015sectors} we divided direction axis
-into four sectors: from -180\textdegree to -90\textdegree, from from -90\textdegree to
+into sectors: from -180\textdegree to -90\textdegree, from from -90\textdegree to
0\textdegree, from 0 to 90\textdegree and from 90\textdegree to 180\textdegree --- and fitted each sector independently. We chose four sector to have one sector for each side of the anemometer.
The best-fit and worst-fit distributions for each axis
are presented in figure~\ref{fig-direction-distributions}.
@@ -446,8 +564,8 @@ corresponding interval.
In this paper we proposed three-axis anemometer that measures wind speed for
each axis independently. We analysed the data collected by this anemometer and
verified that per-axis wind speeds fit into Weibull distribution with the
-largest RMSE of 14\% and wind directions fit into von Mises distibution with
-the largest RMSE of 15\%. We estimated autocovariance functions for wind speed
+largest RMSE of 6.7\% and wind directions fit into von Mises distibution with
+the largest RMSE of 11\%. We estimated autocovariance functions for wind speed
for each axis of the anemometer and used this approximations to simulate wind
flow in Virtual Testbed. The parameters of these functions allow to control
both wind speed and mean direction. The future work is to construct new
diff --git a/papers/Scottanemometerpaper.pdf b/papers/Scottanemometerpaper.pdf
@@ -0,0 +1 @@
+../.git/annex/objects/j2/wx/SHA256E-s418316--b5a1c28612294eeda7bfc2152c9088ffe0578298e66941c1199a63050cbfc0be.pdf/SHA256E-s418316--b5a1c28612294eeda7bfc2152c9088ffe0578298e66941c1199a63050cbfc0be.pdf+
\ No newline at end of file
diff --git a/papers/Yakunin_2017_J._Phys.__Conf._Ser._881_012030.pdf b/papers/Yakunin_2017_J._Phys.__Conf._Ser._881_012030.pdf
@@ -0,0 +1 @@
+../.git/annex/objects/w3/mK/SHA256E-s1546897--c2b64368f683a922fb640b683009bd13bdce5e70116e214af998f8a30bfd61b0.pdf/SHA256E-s1546897--c2b64368f683a922fb640b683009bd13bdce5e70116e214af998f8a30bfd61b0.pdf+
\ No newline at end of file
diff --git a/papers/qt43c525tg.pdf b/papers/qt43c525tg.pdf
@@ -0,0 +1 @@
+../.git/annex/objects/f4/GK/SHA256E-s1239284--42d7ddf55bbec578c9cfbc5e38d98cedd6be856bcdab492fb0e6768f8cd8db4f.pdf/SHA256E-s1239284--42d7ddf55bbec578c9cfbc5e38d98cedd6be856bcdab492fb0e6768f8cd8db4f.pdf+
\ No newline at end of file
