commit 65c6d9394388d2ff9da5e4575e029f3ef91ee1c4
parent 4cb75df8cdcd34f47801c908914a4f226f04d506
Author: Ivan Gankevich <igankevich@ya.ru>
Date: Mon, 29 May 2017 12:56:50 +0300
Add velocity field computation.
Diffstat:
9 files changed, 1057 insertions(+), 38 deletions(-)
diff --git a/arma.org b/arma.org
@@ -4,7 +4,7 @@
#+LATEX_CLASS: scrartcl
#+LATEX_CLASS_OPTIONS:
#+LATEX_HEADER_EXTRA: \input{preamble}
-#+OPTIONS: H:5 num:0 todo:nil toc:nil
+#+OPTIONS: H:5 todo:nil toc:nil
#+PROPERTY: header-args:R :results graphics :exports results
#+begin_abstract
@@ -13,14 +13,15 @@ software-based ship motion modelling applications. These applications generally
use linear wave theory to generate small-amplitude waves programmatically and
determine impact of external excitations on a ship hull. Using linear wave
theory is feasible for ocean waves, but is not accurate for shallow-water and
-storm waves. In this paper we introduce new autoregressive moving-average (ARMA)
-model to cope with these shortcomings. The new model allows to generate waves of
-arbitrary amplitudes, is accurate for both shallow and deep water, and its
-software implementation shows superior performance by relying on fast Fourier
-transform family of algorithms. Integral characteristics of wavy surface
-produced by ARMA model are verified against the ones of real sea surface.
-Despite all its advantages, ARMA model requires a new method to determine wave
-pressures, an instance of which is included in the paper.
+storm waves. To cope with these shortcomings we introduce autoregressive
+moving-average (ARMA) model, which is widely known in oceanography, but rarely
+used for sea wave modelling. The new model allows to generate waves of arbitrary
+amplitudes, is accurate for both shallow and deep water, and its software
+implementation shows superior performance by relying on fast Fourier transform
+family of algorithms. Integral characteristics of wavy surface produced by ARMA
+model are verified against the ones of real sea surface. Despite all its
+advantages, ARMA model requires a new method to determine wave pressures, an
+instance of which is included in the paper.
#+end_abstract
* Introduction
@@ -29,12 +30,12 @@ description of wind waves, which are the major disturbance that displaces the
vessel from equilibrium. Currently, the most popular models for describing wind
waves, are models based on the linear expansion of a stochastic moving surface
as a system of independent random variables. These include models by St. Denis &
-Pearson cite:st1953motions, Rosenblatt cite:rosenblatt1956random, Sveshnikov
-cite:sveshnikov1959determination, and Longuet---Higgins
-cite:longuet1957statistical. The most popular model is that of
-Longuet---Higgins, which is based on a stochastic approximation of the moving
-wave front as a superposition of elementary harmonic waves with random phases
-\(\epsilon_n\) and random amplitudes \(c_n\):
+Pearson\nbsp{}cite:st1953motions, Rosenblatt\nbsp{}cite:rosenblatt1956random,
+Sveshnikov\nbsp{}cite:sveshnikov1959determination, and
+Longuet---Higgins\nbsp{}cite:longuet1957statistical. The most popular model is
+that of Longuet---Higgins, which is based on a stochastic approximation of the
+moving wave front as a superposition of elementary harmonic waves with random
+phases \(\epsilon_n\) and random amplitudes \(c_n\):
\begin{equation}
\label{eq-lhmodel}
\zeta(x,y,t) = \sum\limits_n c_n \cos(u_n x + v_n y - \omega_n t + \epsilon_n),
@@ -146,7 +147,23 @@ becomes an even more significant issue in a nonlinear computation where the wave
model is even more complex. Thus identifying a significantly less time intensive
method for modelling the ambient ocean-wave environment has the potential for
significantly speeding the total simulation process.
+
* Related work
+** Ocean wave modelling
+One of the first application of ARMA model for oceanological processes modelling
+was done in\nbsp{}cite:rozhkov1990. Another approach to modelling wind waves is
+possible in terms of the representation of a stochastic moving surface as a
+linear transformation of white noise with memory. These methods are one of the
+most popular ways of modelling stationary ergodic Gaussian random processes with
+given correlation characteristics\nbsp{}cite:box1976time. However, these methods
+were not used to simulate wind waves for a long time. The first attempts to
+model two-dimensional disturbances were undertaken in the early 70's
+(cf.\nbsp{}cite:kostecki1972stochastic), and the impetus for this was the
+development of the resonance theory of waves in wind. However, the formal
+mathematical framework was developed in\nbsp{}cite:rozhkov1990,gurgenidze1988.
+They built a one-dimensional model of ocean waves, on the basis of an
+autoregressive-moving average (ARMA) model.
+
In\nbsp{}cite:spanos1982arma ARMA model is used to generate time series spectrum of
which is compatible with Pierson---Moskowitz (PM) approximation of ocean wave
spectrum. The authors carry out experiments for one-dimensional AR, MA and ARMA
@@ -158,6 +175,15 @@ number of coefficients than AR model. In\nbsp{}cite:spanos1996efficient the auth
generalise ARMA model coefficients determination formulae for multi-variate
(vector) case.
+In\nbsp{}cite:fusco2010short AR model is used to predict swell waves to control
+wave-energy converters (WEC) in real-time. In order to make WEC more efficient
+its internal oscillator frequency should match the one of ocean waves. The
+authors treat wave elevation as time series and compare performance of AR model,
+neural networks and cyclical models in forecasting time series future values. AR
+model gives the most accurate prediction of low-frequency swell waves for up to
+two typical wave periods. It is an example of successful application of AR
+process to ocean wave modelling.
+
One thing that distinguishes present work with respect to afore-mentioned ones
is the study of three-dimensional (2D in space and 1D in time) ARMA model, which
is mostly a different problem.
@@ -177,20 +203,76 @@ opportunity to visualise output of the programme that allowed to ensure that
generated surface is compatible with real ocean surface, and is not abstract
multi-dimensional stochastic process that is real only statistically.
+** Pressure field determination formulae
+**** Small amplitude waves theory.
+In\nbsp{}cite:stab2012,degtyarev1998modelling,degtyarev1997analysis the
+authors propose a solution for inverse problem of hydrodynamics of potential
+flow in the framework of small-amplitude wave theory (under assumption that wave
+length is much larger than height: \(\lambda \gg h\)). In that case inverse
+problem is linear and reduces to Laplace equation with mixed boundary
+conditions, and equation of motion is solely used to determine pressures for
+calculated velocity potential derivatives. The assumption of small amplitudes
+means the slow decay of wind wave coherence function, i.e. small change of local
+wave number in time and space compared to the wavy surface elevation (\(z\)
+coordinate). This assumption allows to calculate elevation \(z\) derivative as
+\(\zeta_z=k\zeta\), where \(k\) is wave number. In two-dimensional case the
+solution is written explicitly as
+\begin{align}
+ \left.\frac{\partial\phi}{\partial x}\right|_{x,t}= &
+ -\frac{1}{\sqrt{1+\alpha^{2}}}e^{-I(x)}
+ \int\limits_{0}^x\frac{\partial\dot{\zeta}/\partial
+ z+\alpha\dot{\alpha}}{\sqrt{1+\alpha^{2}}}e^{I(x)}dx,\label{eq-old-sol-2d}\\
+ I(x)= & \int\limits_{0}^x\frac{\partial\alpha/\partial z}{1+\alpha^{2}}dx,\nonumber
+\end{align}
+where \(\alpha\) is wave slope. In three-dimensional case solution is written in
+the form of elliptic partial differential equation (PDE):
+\begin{align*}
+ & \frac{\partial^2 \phi}{\partial x^2} \left( 1 + \alpha_x^2 \right) +
+ \frac{\partial^2 \phi}{\partial y^2} \left( 1 + \alpha_y^2 \right) +
+ 2\alpha_x\alpha_y \frac{\partial^2 \phi}{\partial x \partial y} + \\
+ & \left(
+ \frac{\partial \alpha_x}{\partial z} +
+ \alpha_x \frac{\partial \alpha_x}{\partial x} +
+ \alpha_y \frac{\partial \alpha_x}{\partial y}
+ \right) \frac{\partial \phi}{\partial x} + \\
+ & \left(
+ \frac{\partial \alpha_y}{\partial z} +
+ \alpha_x \frac{\partial \alpha_y}{\partial x} +
+ \alpha_y \frac{\partial \alpha_y}{\partial y}
+ \right) \frac{\partial \phi}{\partial y} + \\
+ & \frac{\partial \dot{\zeta}}{\partial z} +
+ \alpha_x \dot{\alpha_x} + \alpha_y \dot{\alpha_y} = 0.
+\end{align*}
+The authors suggest transforming this equation to finite differences and solve
+it numerically.
+
+As will be shown in [[#sec:compare-formulae]] that eqref:eq-old-sol-2d diverges when
+attempted to calculate velocity field for large-amplitude waves, and this is the
+reason that it can not be used together with ARMA model, that generates
+arbitrary-amplitude waves.
+
+**** Linearisation of boundary condition.
+:PROPERTIES:
+:CUSTOM_ID: linearisation
+:END:
+
+LH model allows to derive an explicit formula for velocity field by linearising
+kinematic boundary condition. Velocity potential formula is written as
+\begin{equation*}
+\phi(x,y,z,t) = \sum_n \frac{c_n g}{\omega_n}
+ e^{\sqrt{u_n^2+v_n^2} z}
+ \sin(u_n x + v_n y - \omega_n t + \epsilon_n).
+\end{equation*}
+This formula is differentiated to obtain velocity potential derivatives, which
+are plugged to dynamic boundary condition to obtain pressures.
+
* Three-dimensional ARMA process as a sea wave simulation model
-Another approach to modelling wind waves is possible in terms of the
-representation of a stochastic moving surface as a linear transformation of
-white noise with memory. These methods are one of the most popular ways of
-modelling stationary ergodic Gaussian random processes with given correlation
-characteristics cite:box1976time. However, these methods have were not used
-to simulate wind waves for a long time. The first attempts to model
-two-dimensional disturbances were undertaken in the early 70's (cf.
-cite:kostecki1972stochastic), and the impetus for this was the development of
-the resonance theory of waves in wind. However, the formal mathematical
-framework was developed in cite:rozhkov1990,gurgenidze1988. They built a
-one-dimensional model of ocean waves \(\zeta(t)\), on the basis of an
-autoregressive-moving average (ARMA) model:
+ARMA ocean simulation model defines ocean wavy surface as three-dimensional (two
+dimensions in space and one in time) autoregressive moving average process:
+every surface point is represented as a weighted sum of previous in time and
+space points plus weighted sum of previous in time and space normally
+distributed random impulses. The governing equation for 3-D ARMA process is
\begin{equation}
\zeta_{\vec i}
=
@@ -213,9 +295,9 @@ some substance in water etc.). Equation parameters are AR and MA process
coefficients and order.
Any ARMA process can be uniquely represented as a process moving average and
-autoregression process of general infinite order cite:gurgenidze1988, and the
-parameters of the spectral representation are defined by the rule of division of
-power series (in a rational factorized form cite:rozhkov1990:
+autoregression process of general infinite order\nbsp{}cite:gurgenidze1988, and
+the parameters of the spectral representation are defined by the rule of
+division of power series (in a rational factorized form\nbsp{}cite:rozhkov1990:
\begin{equation*}
S(\omega) =
\frac{\Delta\sigma^2}{\pi}
@@ -553,10 +635,10 @@ In this work both AR and MA model are verified by comparing probability
distributions of different wave characteristics.
*** Verification of wavy surface integral characteristics
-In\nbsp{}cite:рожков1990вероятностные the authors show that several ocean wave
-characteristics (listed in table\nbsp{}[[tab-weibull-shape]]) have Weibull distribution,
-and wavy surface elevation has Gaussian distribution. In order to verify that
-distributions corresponding to generated realisation are correct,
+In\nbsp{}cite:rozhkov1990 the authors show that several ocean wave
+characteristics (listed in table\nbsp{}[[tab-weibull-shape]]) have Weibull
+distribution, and wavy surface elevation has Gaussian distribution. In order to
+verify that distributions corresponding to generated realisation are correct,
quantile-quantile plots are used (plots where analytic quantile values are used
for \(OX\) axis and estimated quantile values for \(OY\) axis). If the estimated
distribution matches analytic then the graph has the form of the straight line.
@@ -657,11 +739,492 @@ because this is the characteristic of the wavy surface (and corresponding AR or
MA process) and is not affected by the type of waves.
** Discussion
+ARMA model, owing to its non-physical nature, does not have the notion of ocean
+wave; it simulates wavy surface as a whole instead. Motions of individual waves
+and their shape are often rough, and the total number of waves can not be
+determined precisely. However, integral characteristics of wavy surface match
+the ones of real ocean waves.
+
+Theoretically, ocean waves themselves can be chosen as ACFs, the only
+pre-processing step is to make them decay exponentially. This may allow
+to generate waves of arbitrary profiles, and is one of the directions of future
+work.
+
* Determining wave pressures for discretely given wavy surface
+Analytic solutions to boundary problems in classical equations are often used to
+study different properties of the solution, and for that purpose general
+solution formula is too difficult to study, as it contains integrals of unknown
+functions. Fourier method is one of the methods to find analytic solutions to
+PDE. It is based on application of Fourier transform to each part of PDE, which
+reduces the equation to algebraic, and the solution is written as inverse
+Fourier transform of some function (which may contain Fourier transforms of
+other functions). Since, it is not possible to write analytic forms of these
+Fourier transforms in all cases, unique solutions are found and their behaviour
+is studied in different domains instead. At the same time, computing discrete
+Fourier transforms on the computer is possible for any discretely defined
+function and efficient when using FFT algorithms. These algorithms use symmetry
+of complex exponentials to decrease asymptotic complexity from
+\(\mathcal{O}(n^2)\) to \(\mathcal{O}(n\log_{2}n)\). So, even if general solution
+contains Fourier transforms of unknown functions, they still can be computed
+numerically, and FFT family of algorithms makes this approach efficient.
+
+Alternative approach to solve PDE is to reduce it to difference equations, which
+are solved by constructing various numerical schemes. This approach leads to
+approximate solution, and asymptotic complexity of corresponding algorithms is
+comparable to that of FFT. For example, stationary elliptic PDE transforms to
+implicit numerical scheme which is solved by iterative method on each step of
+which a tridiagonal of five-diagonal system of algebraic equations is solved by
+Thomas algorithm. Asymptotic complexity of this approach is
+\(\mathcal{O}({n}{m})\), where \(n\)\nbsp{}--- number of wavy surface grid
+points, \(m\)\nbsp{}--- number of iterations. Despite their wide spread,
+iterative algorithms are inefficient on parallel computer architectures; in
+particular, their mapping to co-processors may involve copying data in and out
+of the co-processor in each iteration, which negatively affects their
+performance. At the same time, high number of Fourier transforms in the solution
+is an advantage, rather than a disadvantage. First, solutions obtained by
+Fourier method are explicit, hence their implementations scales with the large
+number of parallel computer cores. Second, there are implementations of FFT
+optimised for different processor architectures as well as co-processors (GPU,
+MIC) which makes it easy to get high performance on any computing platform.
+These advantages substantiate the choice of Fourier method to obtain explicit
+analytic solution to the problem of determining pressures under wavy ocean
+surface.
+
+The problem of finding pressure field under wavy sea surface represents inverse
+problem of hydrodynamics for incompressible inviscid fluid. System of equations
+for it in general case is written as\nbsp{}cite:kochin1966theoretical
+\begin{align}
+ & \nabla^2\phi = 0,\nonumber\\
+ & \phi_t+\frac{1}{2} |\vec{\upsilon}|^2 + g\zeta=-\frac{p}{\rho}, & \text{на }z=\zeta(x,y,t),\label{eq-problem}\\
+ & D\zeta = \nabla \phi \cdot \vec{n}, & \text{на }z=\zeta(x,y,t),\nonumber
+\end{align}
+where \(\phi\)\nbsp{}--- velocity potential, \(\zeta\)\nbsp{}--- elevation (\(z\) coordinate)
+of wavy surface, \(p\)\nbsp{}--- wave pressure, \(\rho\)\nbsp{}--- fluid density,
+\(\vec{\upsilon}=(\phi_x,\phi_y,\phi_z)\)\nbsp{}--- velocity vector, \(g\)\nbsp{}---
+acceleration of gravity, and \(D\)\nbsp{}--- substantial (Lagrange) derivative. The
+first equation is called continuity (Laplace) equation, the second one is the
+conservation of momentum law (the so called dynamic boundary condition); the
+third one is kinematic boundary condition for free wavy surface, which states
+that rate of change of wavy surface elevation (\(D\zeta\)) equals to the change of
+velocity potential derivative along the wavy surface normal
+(\(\nabla\phi\cdot\vec{n}\)).
+
+Inverse problem of hydrodynamics consists in solving this system of equations
+for \(\phi\). In this formulation dynamic boundary condition becomes explicit
+formula to determine pressure field using velocity potential derivatives
+obtained from the remaining equations. So, from mathematical point of view
+inverse problem of hydrodynamics reduces to Laplace equation with mixed boundary
+condition\nbsp{}--- Robin problem.
+
** Two-dimensional case
+:PROPERTIES:
+:CUSTOM_ID: sec:pressure-2d
+:END:
+**** Formula for infinite depth fluid.
+Two-dimensional Laplace equation with Robin boundary condition is written as
+\begin{align}
+ \label{eq-problem-2d}
+ & \phi_{xx}+\phi_{zz}=0,\\
+ & \zeta_t + \zeta_x\phi_x = \frac{\zeta_x}{\sqrt{1 + \zeta_x^2}} \phi_x - \phi_z, & \text{на }z=\zeta(x,t).\nonumber
+\end{align}
+Use Fourier method to solve this problem. Applying Fourier transform to both
+sides of the equation yields
+\begin{equation*}
+ -4 \pi^2 \left( u^2 + v^2 \right)
+ \FourierY{\phi(x,z)}{u,v} = 0,
+\end{equation*}
+hence \(v = \pm i u\). Hereinafter we use the following symmetric form of Fourier
+transform:
+\begin{equation*}
+ \FourierY{f(x,y)}{u,v} =
+ \iint\limits_{-\infty}^{\phantom{--}\infty}
+ f(x,y)
+ e^{-2\pi i (x u + y v)}
+ dx dy.
+\end{equation*}
+We seek solution in the form of inverse Fourier transform
+\(\phi(x,z)=\InverseFourierY{E(u,v)}{x,z}\). Plugging[fn::\(v={-i}{u}\) is not
+applicable because velocity potential must go to nought when depth goes to
+infinity.] \(v={i}{u}\) into the formula yields
+\begin{equation}
+ \label{eq-guessed-sol-2d}
+ \phi(x,z) = \InverseFourierY{e^{2\pi u z}E(u)}{x}.
+\end{equation}
+In order to make substitution \(z=\zeta(x,t)\) not interfere with Fourier
+transforms, we rewrite eqref:eq-guessed-sol-2d as a convolution:
+\begin{equation*}
+ \phi(x,z)
+ =
+ \Fun{z}
+ \ast
+ \InverseFourierY{E(u)}{x},
+\end{equation*}
+where \(\Fun{z}\)\nbsp{}--- a function, form of which is defined in section
+[[#sec:compute-delta]] and which satisfies equation
+\(\FourierY{\Fun{z}}{u}=e^{2\pi{u}{z}}\). Plugging formula \(\phi\) into the boundary
+condition yields
+\begin{equation*}
+ \zeta_t
+ =
+ \left( i f(x) - 1 \right)
+ \left[
+ \Fun{z}
+ \ast
+ \InverseFourierY{2\pi u E(u)}{x}
+ \right],
+\end{equation*}
+where \(f(x)={\zeta_x}/{\sqrt{1+\zeta_x^2}}-\zeta_x\). Applying Fourier transform
+to both sides of this equation yields formula for coefficients \(E\):
+\begin{equation*}
+ E(u) =
+ \frac{1}{2\pi u}
+ \frac{
+ \FourierY{\zeta_t / \left(i f(x) - 1\right)}{u}
+ }{
+ \FourierY{\Fun{z}}{u}
+ }
+\end{equation*}
+Finally, substituting \(z\) for \(\zeta(x,t)\) and plugging resulting equation into
+eqref:eq-guessed-sol-2d yields formula for \(\phi(x,z)\):
+\begin{equation}
+ \label{eq-solution-2d}
+ \boxed{
+ \phi(x,z)
+ =
+ \InverseFourierY{
+ \frac{e^{2\pi u z}}{2\pi u}
+ \frac{
+ \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
+ }{
+ \FourierY{ \Fun{\zeta(x,t)} }{u}
+ }
+ }{x}.
+ }
+\end{equation}
+
+Multiplier \(e^{2\pi{u}{z}}/(2\pi{u})\) makes graph of a function to which Fourier
+transform of which is applied asymmetric with respect to \(OY\) axis. This makes
+it difficult to apply FFT which expects periodic function with nought on both
+ends of the interval. Using numerical integration instead of FFT is not faster
+than solving the initial system of equations with numerical schemes. This
+problem is alleviated by using formula eqref:eq-solution-2d-full for finite
+depth fluid with wittingly large depth \(h\). This formula is derived in the
+following section.
+
+**** Formula for finite depth fluid.
+On the sea bottom vertical fluid velocity component equals nought: \(\phi_z=0\) on
+\(z=-h\), where \(h\)\nbsp{}--- water depth. In this case equation \(v=-{i}{u}\), which came
+from Laplace equation, can not be neglected, hence the solution is sought in the
+following form:
+\begin{equation}
+ \phi(x,z)
+ =
+ \InverseFourierY{
+ \left( C_1 e^{2\pi u z} + C_2 e^{-2\pi u z} \right)
+ E(u)
+ }{x}.
+ \label{eq-guessed-sol-2d-full}
+\end{equation}
+Plugging \(\phi\) into the boundary condition on the sea bottom yields
+\begin{equation*}
+ C_1 e^{-2\pi u h} - C_2 e^{2\pi u h} = 0,
+\end{equation*}
+hence \(C_1=\frac{1}{2}C{e}^{2\pi{u}{h}}\) and
+\(C_2=-\frac{1}{2}C{e}^{-2\pi{u}{h}}\). Constant \(C\) may take arbitrary value
+here, because after plugging it becomes part of unknown coefficients \(E(u)\).
+Plugging formulae for \(C_1\) and \(C_2\) into eqref:eq-guessed-sol-2d-full yields
+\begin{equation*}
+ \phi(x,z) = \InverseFourierY{ \Sinh{2\pi u (z+h)} E(u) }{x}.
+\end{equation*}
+Plugging \(\phi\) into the boundary condition on the free surface yields
+\begin{equation*}
+ \zeta_t = f(x) \InverseFourierY{ 2\pi i u \Sinh{2\pi u (z+h)} E(u) }{x}
+ - \InverseFourierY{ 2\pi u \SinhX{2\pi u (z+h)} E(u) }{x}.
+\end{equation*}
+Here \(\sinh\) and \(\cosh\) give similar results near free surface, and since this
+is the main area of interest in practical applications, we assume that
+\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\). Performing analogous to the
+previous section transformations yields final formula for \(\phi(x,z)\):
+\begin{equation}
+\boxed{
+ \phi(x,z,t)
+ =
+ \InverseFourierY{
+ \frac{\Sinh{2\pi u (z+h)}}{2\pi u}
+ \frac{
+ \FourierY{ \zeta_t / \left(i f(x) - 1\right) }{u}
+ }{
+ \FourierY{ \FunSecond{\zeta(x,t)} }{u}
+ }
+ }{x},
+}
+ \label{eq-solution-2d-full}
+\end{equation}
+where \(\FunSecond{z}\)\nbsp{}--- a function, form of which is defined in section
+[[#sec:compute-delta]] and which satisfies equation
+\(\FourierY{\FunSecond{z}}{u}=\Sinh{2\pi{u}{z}}\).
+
+**** Reducing to the formulae from linear wave theory.
+Check the validity of derived formulae by substituting \(\zeta(x,t)\) with known
+analytic formula for plain waves. Symbolic computation of Fourier transforms in
+this section were performed in Mathematica\nbsp{}cite:mathematica10. In the framework
+of linear wave theory assume that waves have small amplitude compared to their
+lengths, which allows us to simplify initial system of equations
+eqref:eq-problem-2d to
+\begin{align*}
+ & \phi_{xx}+\phi_{zz}=0,\\
+ & \zeta_t = -\phi_z & \text{на }z=\zeta(x,t),
+\end{align*}
+solution to which is written as
+\begin{equation*}
+ \phi(x,z,t)
+ =
+ -\InverseFourierY{
+ \frac{e^{2\pi u z}}{2\pi u}
+ \FourierY{\zeta_t}{u}
+ }{x}
+ .
+\end{equation*}
+Propagating wave profile is defined as \(\zeta(x,t)=A\cos(2\pi(kx-t))\). Plugging
+this formula into eqref:eq-solution-2d yields
+\(\phi(x,z,t)=-\frac{A}{k}\sin(2\pi(kx-t))\Sinh{2\pi{k}{z}}\). In order to reduce
+it to the formula from linear wave theory, rewrite hyperbolic sine in
+exponential form, discard the term containing \(e^{-2\pi{k}{z}}\) as contradicting
+condition \(\phi\underset{z\rightarrow-\infty}{\longrightarrow}0\). Taking real
+part of the resulting formula yields
+\(\phi(x,z,t)=\frac{A}{k}e^{2\pi{k}{z}}\sin(2\pi(kx-t))\), which corresponds to
+the known formula from linear wave theory. Similarly, under small-amplitude
+waves assumption the formula for finite depth fluid eqref:eq-solution-2d-full is
+reduced to
+\begin{equation*}
+ \phi(x,z,t)
+ =
+ -\InverseFourierY{
+ \frac{\Sinh{2\pi u (z+h)}}{2\pi u \Sinh{2\pi u h}}
+ \FourierY{\zeta_t}{u}
+ }{x}.
+\end{equation*}
+Substituting \(\zeta(x,t)\) with propagating plain wave profile formula yields
+\begin{equation}
+ \label{eq-solution-2d-linear}
+ \phi(x,z,t)=\frac{A}{k}
+ \frac{\Sinh{2 \pi k (z+h)}}{ \Sinh{2 \pi k h} }
+ \sin(2 \pi (k x-t)),
+\end{equation}
+which corresponds to the formula from linear wave theory for finite depth fluid.
+
+Different forms of Laplace equation solutions, in which decaying exponent is
+written with either "+" or "-" signs, may cause incompatibilities between
+formulae from linear wave theory and formulae derived in this work, where
+\(\sinh\) is used instead of \(\cosh\). Equality
+\(\frac{\Sinh{2\pi{k}(z+h)}}{\Sinh{2\pi{k}{h}}}\approx\frac{\sinh(2\pi{k}(z+h))}{\sinh(2\pi{k}{h})}\)
+becomes strict on the free surface, and difference between left-hand and
+right-hand sides increases when approaching sea bottom (for sufficiently large
+depth difference near free surface is negligible). So, for sufficiently large
+depth any function (\(\cosh\) or \(\sinh\)) may be used for velocity potential
+computation near free surface.
+
+Reducing eqref:eq-solution-2d и eqref:eq-solution-2d-full to the known formulae
+from linear wave theory shows, that formula for infinite depth
+eqref:eq-solution-2d is not suitable to compute velocity potentials with Fourier
+method, because it does not have symmetry, which is required for Fourier
+transform. However, formula for finite depth can be used instead by setting \(h\)
+to some characteristic water depth. For standing wave reducing to linear wave
+theory formulae is made under the same assumptions.
+
** Three-dimensional case
-** Evaluation
-** Discussion
+Three-dimensional version of eqref:eq-problem is written as
+\begin{align}
+ \label{eq-problem-3d}
+ & \phi_{xx} + \phi_{yy} + \phi_{zz} = 0,\\
+ & \zeta_t + \zeta_x\phi_x + \zeta_y\phi_y
+ =
+ \frac{\zeta_x}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_x
+ +\frac{\zeta_y}{\SqrtZeta{1 + \zeta_x^2 + \zeta_y^2}} \phi_y
+ - \phi_z, & \text{на }z=\zeta(x,y,t).\nonumber
+\end{align}
+Again, use Fourier method to solve it. Applying Fourier transform to both sides
+of Laplace equation yields
+\begin{equation*}
+ -4 \pi^2 \left( u^2 + v^2 + w^2 \right)
+ \FourierY{\phi(x,y,z)}{u,v,w} = 0,
+\end{equation*}
+hence \(w=\pm{i}\sqrt{u^2+v^2}\). We seek solution in the form of inverse Fourier
+transform \(\phi(x,y,z)=\InverseFourierY{E(u,v,w)}{x,y,z}\). Plugging
+\(w=i\sqrt{u^2+v^2}=i\Kveclen\) into the formula yields
+\begin{equation*}
+ \phi(x,y,z) = \InverseFourierY{
+ \left(
+ C_1 e^{2\pi \Kveclen z}
+ -C_2 e^{-2\pi \Kveclen z}
+ \right)
+ E(u,v)
+ }{x,y}.
+\end{equation*}
+Plugging \(\phi\) into the boundary condition on the sea bottom (analogous to
+two-dimensional case) yields
+\begin{equation}
+ \label{eq-guessed-sol-3d}
+ \phi(x,y,z) = \InverseFourierY{
+ \Sinh{2\pi \Kveclen (z+h)} E(u,v)
+ }{x,y}.
+\end{equation}
+Plugging \(\phi\) into the boundary condition on the free surface yields
+\begin{equation*}
+ \arraycolsep=1.4pt
+ \begin{array}{rl}
+ \zeta_t = & i f_1(x,y) \InverseFourierY{2 \pi u \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
+ + & i f_2(x,y) \InverseFourierY{2 \pi v \Sinh{2\pi \Kveclen (z+h)}E(u,v)}{x,y} \\
+ - & \InverseFourierY{2 \pi \Kveclen \SinhX{2\pi \Kveclen (z+h)}E(u,v)}{x,y}
+ \end{array}
+\end{equation*}
+where \(f_1(x,y)={\zeta_x}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_x\) and
+\(f_2(x,y)={\zeta_y}/{\SqrtZeta{1+\zeta_x^2+\zeta_y^2}}-\zeta_y\).
+
+Like in Section\nbsp{}[[#sec:pressure-2d]] we assume that
+\(\Sinh{2\pi{u}(z+h)}\approx\SinhX{2\pi{u}(z+h)}\) near free surface, but in
+three-dimensional case this is not enough to solve the problem. In order to get
+analytic formula for coefficients \(E\) we need to assume, that all Fourier
+transforms in the equation have radially symmetric kernels, i.e. replace \(u\)
+and \(v\) with \(\Kveclen\). There are two points supporting this assumption.
+First, in numerical implementation integration is done over positive wave
+numbers, so the sign of \(u\) and \(v\) does not affect the solution. Second,
+the rate growth of \(\cosh\) term of the integral kernel is much higher than the
+one of \(u\) or \(\Kveclen\), so the substitution has small effect on the
+magnitude of the solution. Despite these two points, a use of more
+mathematically rigorous approach would be preferable.
+
+Making the replacement, applying Fourier transform to both sides of the equation
+and plugging the result into eqref:eq-guessed-sol-3d yields formula for
+\(\phi\):
+\begin{equation*}
+ \phi(x,y,z,t) = \InverseFourierY{
+ \frac{ \Sinh{\smash{2\pi \Kveclen (z+h)}} }{ 2\pi\Kveclen }
+ \frac{ \FourierY{ \zeta_t / \left( i f_1(x,y) + i f_2(x,y) - 1 \right)}{u,v} }
+ { \FourierY{\mathcal{D}_3\left( x,y,\zeta\left(x,y\right) \right)}{u,v} }
+ }{x,y},
+\end{equation*}
+where \(\FourierY{\mathcal{D}_3\left(x,y,z\right)}{u,v}=\Sinh{\smash{2\pi\Kveclen{}z}}\).
+
+** Evaluation and discussion
+:PROPERTIES:
+:CUSTOM_ID: sec:compare-formulae
+:END:
+
+Comparing obtained generic formulae eqref:eq-solution-2d and
+eqref:eq-solution-2d-full to the known formulae from linear wave theory allows
+to see the difference between velocity fields for both large and small amplitude
+waves. In general analytic formula for velocity potential in not known, even for
+plain waves, so comparison is done numerically. Taking into account conclusions
+of [[#sec:pressure-2d]], only finite depth formulae are compared.
+
+**** The difference with linear wave theory formulae.
+In order to obtain velocity potential fields, ocean wavy surface was generated
+by AR model with varying wave amplitude. In numerical implementation wave
+numbers in Fourier transforms were chosen on the interval from \(0\) to the
+maximal wave number determined numerically from the obtained wavy surface.
+Experiments were conducted for waves of both small and large amplitudes.
+
+The experiment showed that velocity potential fields produced by formula
+eqref:eq-solution-2d-full for finite depth fluid and formula
+eqref:eq-solution-2d-linear from linear wave theory are qualitatively different
+(fig.\nbsp{}[[fig-potential-field-nonlinear]]). First, velocity potential contours
+have sinusoidal shape, which is different from oval shape described by linear
+wave theory. Second, velocity potential decays more rapidly than in linear wave
+theory as getting closer to the bottom, and the region where the majority of
+wave energy is concentrated is closer to the wave crest. Similar numerical
+experiment, in which all terms of eqref:eq-solution-2d-full that are neglected
+in the framework of linear wave theory are eliminated, shows no difference (as
+much as machine precision allows) in resulting velocity potential fields.
+
+#+name: fig-potential-field-nonlinear
+#+header: :width 8 :height 11
+#+begin_src R :file build/plain-wave-velocity-field-comparison.pdf
+source(file.path("R", "velocity-potentials.R"))
+par(pty="s")
+nlevels <- 41
+levels <- pretty(c(-200,200), nlevels)
+palette <- colorRampPalette(c("blue", "lightyellow", "red"))
+col <- palette(nlevels-1)
+
+# linear solver
+par(fig=c(0,0.95,0,0.5),new=TRUE)
+arma.plot_velocity_potential_field(
+ file.path("build", "plain_wave_linear_solver"),
+ levels=levels,
+ col=col
+)
+
+# high-amplitude solver
+par(fig=c(0,0.95,0.5,1),new=TRUE)
+arma.plot_velocity_potential_field(
+ file.path("build", "plain_wave_high_amplitude_solver"),
+ levels=levels,
+ col=col
+)
+
+# legend 1
+par(pty="m",fig=c(0.80,1,0.5,1), new=TRUE)
+arma.plot_velocity_potential_field_legend(
+ levels=levels,
+ col=col
+)
+
+# legend 2
+par(pty="m",fig=c(0.80,1,0,0.5), new=TRUE)
+arma.plot_velocity_potential_field_legend(
+ levels=levels,
+ col=col
+)
+#+end_src
+
+#+caption: Velocity potential field of propagating wave \(\zeta(x,y,t) = \cos(2\pi x - t/2)\). Field produced by formula eqref:eq-solution-2d-full (top) and linear wave theory formula (bottom).
+#+label: fig-potential-field-nonlinear
+#+attr_latex: :width \textwidth
+#+RESULTS: fig-potential-field-nonlinear
+[[file:build/plain-wave-velocity-field-comparison.pdf]]
+
+**** The difference with small-amplitude wave theory.
+The experiment, in which velocity fields produced numerically by different
+formulae were compared, shows that velocity fields produced by formula
+eqref:eq-solution-2d-full and eqref:eq-old-sol-2d correspond to each other for
+small-amplitude waves. Two ocean wavy surface realisations were made by AR
+model: one containing small-amplitude waves, other containing large-amplitude
+waves. Integration in formula eqref:eq-solution-2d-full was done over wave
+numbers range extracted from the generated wavy surface. For small-amplitude
+waves both formulae showed comparable results (the difference in the velocity is
+attributed to the stochastic nature of AR model), whereas for large-amplitude
+waves stable velocity field was produced only by formula
+eqref:eq-solution-2d-full (fig.\nbsp{}[[fig-velocity-field-2d]]). So, generic
+formula eqref:eq-solution-2d-full gives satisfactory results without restriction
+on wave amplitudes.
+
+#+name: fig-velocity-field-2d
+#+header: :width 8 :height 11
+#+begin_src R :file build/large-and-small-amplitude-velocity-field-comparison.pdf
+source(file.path("R", "velocity-potentials.R"))
+linetypes = c("solid", "dashed")
+par(mfrow=c(2, 1))
+arma.plot_velocity(
+ file.path("data", "velocity", "low-amp"),
+ file.path("data", "velocity", "low-amp-0"),
+ linetypes=linetypes,
+ ylim=c(-2,2)
+)
+arma.plot_velocity(
+ file.path("data", "velocity", "high-amp"),
+ file.path("data", "velocity", "high-amp-0"),
+ linetypes=linetypes,
+ ylim=c(-2,2)
+)
+#+end_src
+
+#+label: fig-velocity-field-2d
+#+caption: Comparison of velocity field on the ocean wavy surface obtained by generic formula (\(u_1\)) and formula for small-amplitude waves (\(u_2\)). Velocity field for realisations containing small-amplitude (top) and large-amplitude (bottom) waves.
+#+attr_latex: :width \textwidth
+#+RESULTS: fig-velocity-field-2d
+[[file:build/large-and-small-amplitude-velocity-field-comparison.pdf]]
+
* High-performance software implementation for heterogeneous platforms
** White noise generation
** Wavy surface generation
@@ -669,6 +1232,20 @@ MA process) and is not affected by the type of waves.
** Evaluation
** Discussion
* Conclusion
+Research results allow to conclude that a problem of determining pressures under
+sea surface can be solved analytically without assumptions of linear and
+small-amplitude wave theories. This solution coupled with ARMA ocean wave
+simulation model, capable of generating waves of arbitrary amplitudes, can be
+used to determine the impact of wave oscillations on the dynamic marine object
+in a sea way, and give more precise results than analogous solution for
+small-amplitude waves.
+
+Results of the numerical experiments allow to conclude that wavy surface
+generation as well as pressure computation can be efficiently implemented via
+fast Fourier transforms, and long simulation session can be conducted.
+
+The developed mathematical apparatus and its numerical implementation can become
+a base of virtual testbed for marine objects dynamics studies.
bibliographystyle:plain
bibliography:refs.bib
diff --git a/data/performance/factory-vs-openmp.csv b/data/performance/factory-vs-openmp.csv
@@ -0,0 +1,11 @@
+nt,nx,ny,openmp,factory
+2000,32,32,3.20,2.31
+4000,32,32,6.40,4.46
+6000,32,32,10.75,7.35
+8000,32,32,12.43,8.62
+10000,32,32,15.57,11.38
+12000,32,32,18.76,13.33
+14000,32,32,21.96,15.52
+16000,32,32,25.19,17.80
+18000,32,32,28.34,20.12
+20000,32,32,31.44,22.53
diff --git a/data/performance/overlap-factory.csv b/data/performance/overlap-factory.csv
@@ -0,0 +1,7 @@
+t,mark
+856632019315375,
+856632366634850,G₀
+856632502311413,W₀
+856653615740198,W₁
+856653615746769,G₁
+856653634687539,
diff --git a/data/performance/overlap-openmp.csv b/data/performance/overlap-openmp.csv
@@ -0,0 +1,7 @@
+t,mark
+856601220966013,
+856601399800174,G₀
+856619491863579,G₁
+856619505966724,W₀
+856631949369101,W₁
+856632017128894,
diff --git a/data/velocity/high-amp b/data/velocity/high-amp
@@ -0,0 +1,93 @@
+-1.1036507585724333
+-0.5111703936169758
+-0.008759734302216684
+0.2636739870924831
+0.42546308840010627
+0.640452641258267
+0.752327677572453
+0.8318340102844313
+0.8217250351335857
+0.5910327405679315
+0.4430077368070742
+0.5135105237748081
+0.5338033128399037
+0.5753903945535445
+0.5595580697857758
+0.49004939338109643
+0.4722842523464433
+0.40399302079979593
+0.36367499440116224
+0.4828607782519232
+0.6459732315860856
+0.49534926315134
+0.2555553372546335
+0.13589069344080335
+-0.024144668792737287
+0.09206164699370954
+0.28568992264214194
+0.4038382687667773
+0.21478576582942913
+0.04247449724111433
+-0.08633128874141174
+-0.544261552182005
+-0.3888688936004529
+0.8860288296439196
+0.1290043023201055
+-0.17326707635388175
+0.028928265443251826
+0.20033382071469324
+0.2923401243251371
+0.30532161121353674
+0.0771203741320842
+-1.0589631295020974
+-1.0489903530343474
+-0.3266678398246404
+0.07168396548193363
+0.4694645528632
+0.5593195945783879
+0.2424892779425071
+-0.18637839549681515
+-0.4784329254037924
+-0.5679499346748907
+-0.2911135966909082
+-0.000017996099670610878
+0.4362851346611899
+0.8296934305133177
+0.8508698124871845
+0.22531658610201233
+-0.0719461823062096
+0.27917807537566686
+0.4270971178073046
+0.5295261391618792
+0.6618518388522279
+0.5185922202594925
+0.12169369839531728
+-0.15772392465705304
+-0.44662511832662377
+0.052486433631477626
+0.6067445526637256
+0.9536347814759252
+1.0041676067589396
+0.6170873317118108
+0.1602466222593934
+-0.1208137101787172
+-0.09416729659359119
+0.329469654414881
+0.7847178904759198
+1.0845067788352556
+1.081204414765634
+0.9285402309817897
+0.9036274636241375
+0.9416853928532346
+0.9861932612372499
+1.078129573690889
+1.053629580169305
+1.2261137683966643
+1.2035470644699506
+0.7120873747052325
+0.9295311022106962
+0.55538937909681
+0.0792701520208196
+-0.14774888384946672
+-0.25546185646952546
+0.16840510816125093+
\ No newline at end of file
diff --git a/data/velocity/high-amp-0 b/data/velocity/high-amp-0
@@ -0,0 +1,93 @@
+0.736636
+-0.00113116
+-0.728071
+-4.24231
+-16.1775
+-68.5454
+-40.1521
+-16.7891
+-6.32929
+-7.50372
+-6.3322
+-1.97682
+0.718014
+1.99714
+4.07408
+9.38275
+4.09194
+0.562338
+0.869305
+1.07734
+3.59298
+13.1335
+14.6075
+8.57484
+3.0408
+1.01585
+3.19161
+11.5383
+21.6267
+74.1945
+24.0001
+14.9655
+5.75329
+15.8557
+44.9732
+103.304
+295.12
+363.684
+110.686
+44.9031
+23.0951
+15.8491
+23.9933
+49.2657
+133.528
+65.5187
+21.4721
+6.40511
+1.50812
+0.198471
+0.00627411
+-0.608862
+-1.19253
+1.06313
+0.775713
+-0.235812
+-1.47255
+-4.82519
+-14.7861
+-51.3172
+-158.517
+-31.1255
+-16.7888
+-6.69476
+-29.3019
+-84.2644
+-261.993
+-806.782
+-223.341
+-61.0173
+-17.4345
+-6.82616
+-17.8384
+-56.6058
+-224.854
+-528.483
+-445.435
+-233.421
+-329.792
+-327.045
+-158.185
+-35.9012
+-17.3041
+-15.5498
+-4.70248
+-2.70484
+-7.84366
+-30.4055
+-112.084
+-112.078
+-40.3093
+-9.1616
+-6.12593
diff --git a/data/velocity/low-amp b/data/velocity/low-amp
@@ -0,0 +1,93 @@
+-0.18592339204731947
+0.08589265162801145
+0.5255422147798241
+0.8842451508793726
+0.9390034170536219
+0.6011925446051457
+0.05920428715349717
+-0.33068892584188037
+-0.3604264823658643
+-0.11407859921179717
+0.18446019506394268
+0.32332697324708964
+0.15712744433975537
+-0.16992881547290617
+-0.41718988674286933
+-0.3517626425106538
+-0.051899193539124575
+0.3566189412184856
+0.6262446899151234
+0.6735966604305366
+0.5749548388989701
+0.33792537391335536
+0.08205216103274916
+-0.058150420938388805
+-0.07356047736023474
+0.006749725139703499
+0.17434076968780382
+0.40097947355480396
+0.6759436961161973
+0.8415491711789904
+1.0242453980394344
+1.017384168873244
+0.8569097971666798
+0.6858771132491988
+0.5154074300941885
+0.4034179498800563
+0.32209503649896243
+0.2911320408581537
+0.30307260982309087
+0.3789997315567123
+0.526105727036769
+0.6579665639767911
+0.7343726519068272
+0.6679820590809306
+0.44950667885837664
+0.14482607563095914
+-0.08907848319393948
+-0.1801527388550217
+-0.07518625725174084
+0.017726229024175064
+0.1342101316282809
+0.21682142756675443
+0.26035755349654094
+0.4975339536701552
+0.8452233377285836
+1.2518110934811388
+1.6208383822278216
+1.6165870249478722
+1.3888017540068924
+1.0012642709289072
+0.7288789081680658
+0.6939916490420551
+0.8946663344476296
+1.1782170587224947
+1.2075164408069565
+1.0486821137165703
+0.8227172499238671
+0.7028929877730287
+0.7352959416561036
+0.8994114120761405
+1.0512005324995557
+1.0546298296452143
+0.8347920853806886
+0.5770451508234755
+0.46770410631213954
+0.5257968444744331
+0.6894485082008976
+0.8370807257958566
+0.8063142625631204
+0.6212046022689293
+0.3900543562304422
+0.3099053530529887
+0.4745182078767032
+0.7921187158836802
+1.1136843855486644
+1.3333529971825495
+1.4006504984660888
+1.3458632458805744
+1.2764773279956882
+1.2290685020468015
+1.2123009335338795
+1.1372806565191205
+0.9588315439221115+
\ No newline at end of file
diff --git a/data/velocity/low-amp-0 b/data/velocity/low-amp-0
@@ -0,0 +1,93 @@
+-0.0917072
+0.053641
+0.171918
+0.193579
+0.130987
+0.00479723
+-0.226553
+-0.501343
+-0.788274
+-0.699502
+-0.384301
+-0.220301
+-0.157091
+-0.0589546
+-0.14061
+-0.401631
+-0.614626
+-0.821389
+-0.699904
+-0.424918
+-0.344167
+-0.296599
+-0.431693
+-0.678114
+-1.17528
+-1.5362
+-1.14082
+-0.579577
+-0.175253
+0.246136
+0.572435
+0.752649
+0.340125
+-0.114834
+-0.0754783
+0.111857
+0.321068
+0.468458
+0.495278
+0.281224
+-0.193674
+-0.389608
+-0.58564
+-0.677042
+-0.6127
+-0.545627
+-0.166323
+0.234772
+0.370243
+0.245047
+0.0655986
+0.0063266
+0.250332
+0.70624
+0.93047
+1.03047
+1.11407
+1.08168
+0.709051
+0.413053
+0.51106
+0.626707
+0.796466
+1.16404
+1.47804
+1.49801
+1.32694
+0.979558
+0.554724
+0.622674
+0.954948
+1.21362
+1.48913
+1.59508
+1.43442
+1.01447
+0.764845
+0.749179
+0.950156
+1.1008
+1.12541
+0.877175
+0.624883
+0.768627
+1.10514
+1.44319
+1.66924
+1.64651
+1.41551
+1.13722
+1.13222
+1.38316
+1.68314
diff --git a/refs.bib b/refs.bib
@@ -23,7 +23,7 @@
@article{longuet1957statistical,
title={The statistical analysis of a random, moving surface},
- author={Longuet-Higgins, Michael S},
+ author={Longuet---Higgins, Michael S},
journal={Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences},
volume={249},
number={966},
@@ -137,7 +137,7 @@
}
@book{rozhkov1990,
- author = {Rozhkov, V. A. and Y. A. Trapeznikov},
+ author = {Rozhkov, V. A. and Trapeznikov, Yu. A. },
title = {Probabilistic models of oceanographic processes},
address = {Leningrad},
publisher = {Gidrometeoizdat},
@@ -150,3 +150,45 @@
year={1972},
school={CTO, Gdansk}
}
+
+@article{fusco2010short,
+ title = {Short-term wave forecasting for real-time control of wave energy converters},
+ author = {Fusco, Francesco and Ringwood, John V},
+ journal = {IEEE Transactions on sustainable energy},
+ volume = {1},
+ number = {2},
+ pages = {99--106},
+ year = {2010},
+ publisher = {IEEE}
+}
+
+@InProceedings{ stab2012,
+ author = {A. Degtyarev and I. Gankevich},
+ title = {Evaluation of hydrodynamic pressures for autoregression
+ model of irregular waves},
+ booktitle = {Proceedings of 11\textsuperscript{th} International
+ Conference on Stability of Ships and Ocean Vehicles,
+ Athens},
+ year = 2012,
+ pages = {841--852}
+}
+
+@techreport{degtyarev1997analysis,
+ author={Degtyarev, A. and Boukhanovsky, A.},
+ title={Analysis of Peculiarities of Ship-Environmental Interaction},
+ month={September},
+ year={1997},
+ address={Ship Stability Research Center, Glasgow},
+ number={09-97-1AB-1VA},
+ institution={Strathclyde University}
+}
+
+@InProceedings{degtyarev1998modelling,
+ author={Degtyarev, A. B. and and Podoliakin, A. B.},
+ title={Simulation modelling of ship behaviour in real sea waves (in Russian)},
+ address={Saint-Petersburg},
+ volume={V},
+ pages={416--423},
+ booktitle={Proc. of II int. conf. on shipbuilding (ISC'98)},
+ year={1998}
+}+
\ No newline at end of file
