commit 7c184cec46392f521e67086de1c2597691ae6a08
parent 5129cf5a12e337aa1a5f3088edd1d768ac4c7eb8
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Mon, 16 Mar 2020 12:17:52 +0300
Copy everything from math repo.
Diffstat:
4 files changed, 582 insertions(+), 8 deletions(-)
diff --git a/main.tex b/main.tex
@@ -1,10 +1,18 @@
\documentclass[runningheads]{llncs}
\usepackage{amsmath}
+\usepackage{amssymb}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage{url}
+\newcommand{\Jacobian}{\mathbb{J}}
+\newcommand{\Real}[1]{\operatorname{Re}#1}
+\newcommand{\Imag}[1]{\operatorname{Im}#1}
+\newcommand{\VectorL}[1]{\left[\begin{array}{l}#1\end{array}\right]}
+\newcommand{\VectorR}[1]{\left[\begin{array}{r}#1\end{array}\right]}
+\newcommand{\DerivativeT}[1]{\frac{\partial{#1}}{\partial{}t}}
+
\begin{document}
\title{TODO\thanks{Supported by Saint Petersburg State University (grants
@@ -25,14 +33,14 @@
\institute{Saint Petersburg State University\\
7-9 Universitetskaya Emb., St Petersburg 199034, Russia\\
-\email{st049350@student.spbu.ru},\\
-\email{a.degtyarev@spbu.ru},\\
-\email{st047824@student.spbu.ru},\\
-\email{i.gankevich@spbu.ru},\\
-\email{st047437@student.spbu.ru},\\
-\email{st016177@student.spbu.ru},\\
-\email{v.khramusin@spbu.ru}
-\url{https://spbu.ru/}}
+ \email{st049350@student.spbu.ru},\\
+ \email{a.degtyarev@spbu.ru},\\
+ \email{st047824@student.spbu.ru},\\
+ \email{i.gankevich@spbu.ru},\\
+ \email{st047437@student.spbu.ru},\\
+ \email{st016177@student.spbu.ru},\\
+ \email{v.khramusin@spbu.ru}\\
+ \url{https://spbu.ru/}}
\maketitle
@@ -48,7 +56,325 @@ TODO
\end{abstract}
\section{Introduction}
+
+There are two mathematical models that describe rigid body motion and fluid
+particle motion: equations of translational and angular motion of the rigid
+body (Newton's Second Law) and Gerstner equations for ocean waves (which are
+solutions to linearised equations of motions for fluid particles). Usually, we
+use these models independently to generate incident ocean waves and then
+compute body motions caused by these waves. To measure the effect of still
+fluid on an oscillating rigid body (radiation forces) and the effect of fluid
+particles hitting the body (diffraction forces), we use added masses and
+damping coefficients~--- simplified formulae derived for small-amplitude
+oscillatory motion.
+
+But, what if we want to simulate large-amplitude rigid body motion with greater
+accuracy? There are two possible ways. First, we may use numerical methods such
+as Reynolds-Averaged Navier Stokes (RANS) method~\cite{wackers2011rans}. This
+method is accurate, can be used for viscous fluid, but not the most
+computationally efficient. Second, we may solve Gerstner equations with
+appropriate boundary condition and use the solution to compute both rigid body
+and fluid particle motion around it. This paper explores this second option.
+
+%and compares the results to the results of the method of added masses and numerical solver.
+
+
+
\section{Methods}
+
+\subsection{Equations of motions with a moving surface boundary}
+
+An oscillating rigid body that floats in the water and experiences incident
+waves both reflects existing waves and generates new waves:
+\begin{itemize}
+\item fluid particles hit the body, causing it to move, and then reflect from it;
+\item moving body hits fluid particles and makes them move.
+\end{itemize}
+Both wave reflection and generation have the same nature~--- they are
+caused by the collision of the particles and the body~--- hence we
+describe them by the same set of formulae. Hereinafter we borrow the
+mathematical notation for Lagrangian description of the flow
+from~\cite{nouguier2015}.
+
+In Lagrangian description of flow instantaneous particle coordinates
+\(\vec{R}=(x,y,z)\) depend on particle positions at rest (independent initial
+coordinates) \(\vec{\zeta}=(\alpha,\beta,\delta)\) and time \(t\),
+i.e.~\(\vec{R}=\vec{R}(\alpha,\beta,\delta,t)\). Using this notation
+equation of motion (conservation of momentum) is written as
+\begin{equation}
+\vec{R}_{tt} + g \hat{z} + \frac{1}{\rho} \nabla_{\vec{R}} p = 0,
+\label{eq-momentum}
+\end{equation}
+where \(\vec{R}\)~--- particle coordinates, \(\hat{z}\)~---
+unit vector in the direction of positive \(z\), \(p\)~--- pressure,
+\(\rho\)~--- fluid density, \(g\)~--- gravitational acceleration.
+Continuity equation (conservation of mass) is written as
+\begin{equation*}
+\left|\Jacobian\right| = 1, \qquad
+\frac{\partial}{\partial t}\left|\Jacobian\right| = 0, \qquad
+\Jacobian = \left[
+\begin{array}{lll}
+x_\alpha & y_\alpha & z_\alpha \\
+x_\beta & y_\beta & z_\beta \\
+x_\delta & y_\delta & z_\delta \\
+\end{array}
+\right].
+\end{equation*}
+Multiplying both sides of \eqref{eq-momentum} by \(\Jacobian\) and noting that
+\(\nabla_{\vec{\zeta}}p =\Jacobian\nabla_{\vec{R}}p\) gives
+\begin{equation*}
+\Jacobian\vec{R}_{tt} +
+g\vec\nabla\left( \vec{R} \cdot \hat{z} \right) +
+\frac{1}{\rho}\vec\nabla p
+\end{equation*}
+Following \cite{nouguier2015} we seek solution to this equation in the form
+of a simultaneous perturbation expansion for position, pressure, and the
+vorticity function:
+\begin{equation*}
+\begin{aligned}
+\vec{R} &= \vec{R}_0 + \vec{R}_1 + \vec{R}_2 + ... \\
+p &= p_a - \rho g \delta + p_1 + p_2 + ...
+\end{aligned}
+\end{equation*}
+Zeroth order terms are related to particles positions at rest:
+\begin{equation*}
+\begin{aligned}
+\vec{R}_0 &= \vec{\zeta} \\
+p_0 &= p_a - \rho g \delta
+\end{aligned}
+\end{equation*}
+First-order terms are solutions to linearised equation of motion and equation of
+continuity:
+\begin{equation*}
+\begin{aligned}
+& \vec{R}_{1tt} + g\vec\nabla\left( \vec{R}_1 \cdot \hat{z} \right) +
+\frac{1}{\rho}\vec\nabla p_1 = 0 \\
+& \vec\nabla \cdot \vec{R}_1 = 0
+\end{aligned}
+\end{equation*}
+We seek solutions of the form \(\vec{R}_1=\nabla w\) to make the flow
+irrotational. Plugging this form into the equations gives
+\begin{equation}
+\label{eq-mass-momentum}
+\begin{aligned}
+& \vec\nabla\left( w_{tt} + g w_\delta + p_1/\rho \right) = 0 \\
+& \Delta w = 0
+\end{aligned}
+\end{equation}
+The first equation denotes conservation of momentum (Newton's second law,
+equation of motion) and the second equation denotes conservation of mass
+(equation of continuity).
+
+When we have no boundary condition we seek solutions of the form
+\begin{equation}
+\label{eq-inverse-fourier}
+w\left(\alpha,\beta,\delta\right) =
+\Real{
+f(u,v)\exp\left(iu\alpha+iv\beta+k\delta-i\omega{}t\right)
+},
+\end{equation}
+where \(u\) and \(v\) are wave numbers. We plug \eqref{eq-inverse-fourier} into
+continuity equation \eqref{eq-mass-momentum} where \(p_1\) is constant and get
+\(k=\sqrt{u^2+v^2}\). That means that expression~\eqref{eq-inverse-fourier}
+is the solution to this equation when \(k\) is wave vector magnitude,
+i.e.~\(w\) decays exponentially with increasing water depth multiplied by
+wave vector magnitude.
+
+Then we plug \eqref{eq-inverse-fourier} into equation of motion
+\eqref{eq-mass-momentum} and get \(\omega^2=gk\), which is dispersion relation
+from classic linear wave theory. That means that the expression is the solution
+to this equation when angular frequency depends on the wave number,
+i.e.~waves of different lengths have different phase velocities.
+
+Before solving this system of equations for an arbitrary moving surface
+boundary, we consider particular cases to substantiate the choice of the form of
+the solution.
+
+We use \(\Real{\exp\left(iu\alpha+iv\beta+k\delta-i\omega{}t\right)}\) to
+describe fluid particle potential. Here \(u\) and \(v\) are wave numbers,
+\(\omega\) is angular frequency, and \(k\) is wave vector. This notation makes
+formulae short and is equivalent to the description that uses traditional
+harmonic functions. This notation allows for easy transition to irregular
+waves via Fourier transforms which are essential for fast computations. Such
+solutions will be studied in future work.
+
+\subsection{Stationary surface boundary}
+
+In this section we explore solutions stationary surface boundary in a form of
+infinite plain surface. On such a boundary the projection of particle velocity
+to the surface normal is nought. We write boundary conditions and
+corresponding solutions for different orientations of this boundary and then
+generalise these solutions to a parametric surface.
+
+\subsubsection{Infinite wall}
+
+On a vertical surface the boundary condition is written as
+\begin{equation*}
+\frac{d}{d t} \vec\nabla w \cdot \vec{n} =
+\frac{d}{d t} \frac{\partial}{\partial\alpha} w = 0;
+\qquad
+\alpha = \alpha_0;
+\qquad
+\vec{n} = \VectorR{1\\0\\0}.
+\end{equation*}
+Here we consider only \(\alpha\) coordinate, the derivations for \(\beta\) are
+similar. The potential of incident fluid particle has the form
+\begin{equation*}
+w\left(\alpha,\beta,\delta,t\right) =
+\exp\left(k\delta - i\omega t\right)
+\exp\left(iu\alpha + iv\beta\right).
+\end{equation*}
+Velocity vector of this particle is
+\begin{equation*}
+\begin{aligned}
+&
+\frac{d}{dt} \vec\nabla w =
+i\omega \left( \vec{d}_k+i\vec{d}_i \right)
+\exp\left(k\delta - i\omega t\right)
+\exp\left(iu\alpha + iv\beta\right);
+\\
+&
+\vec{d}_k = \VectorR{0\\0\\k};
+\qquad
+\vec{d}_i = \VectorR{u\\v\\0}.
+\end{aligned}
+\end{equation*}
+where \(\vec{d}_i\) is horizontal incident wave direction and \(\vec{d}_k\) is
+a vector that contains amplitude damping coefficient. (We use the vector
+instead of the scalar to shorten mathematical notation, otherwise we would have write a
+separate formula for vertical coordinate.)
+%In addition to this, it makes the solution for finite depth look the same
+%way as all other solutions.
+The law of reflection
+states that the angle of incidence equals the angle of reflection. Then
+the direction of reflected wave is%
+\footnote{Initially, we included the third
+component of incident wave direction making the
+vector complex-valued, however, the solution blew up as a result of mixing real
+and imaginary parts in dot products involving complex-valued vectors. The
+problem was solved by reflecting in two dimensions which is intuitive for
+ocean waves, but not for particles.}
+\begin{equation*}
+\vec{d}_r =
+\vec{d}_i-\vec{d}_s =
+\vec{d}_i-2\vec{n}\left(\vec{d}_i\cdot\vec{n}\right) =
+\VectorR{-u\\v\\0}.
+\end{equation*}
+We seek solutions of the form
+\begin{equation}
+\label{eq-w-alpha-0}
+w\left(\alpha,\beta,\delta,t\right) =
+\left[ C_1 \exp\left(iu\alpha\right) + C_2 \exp\left(-iu\alpha\right) \right]
+\exp\left(k\delta - i\omega t\right) \exp\left(iv\beta\right).
+\end{equation}
+We plug this expression into the boundary condition and get
+\begin{equation*}
+C_1\exp\left(i u\alpha_0 \right) -
+C_2\exp\left(-i u\alpha_0 \right) = 0,
+\end{equation*}
+hence \(C_1=C\exp(-iu\alpha_0)\) and \(C_2=-C\exp(iu\alpha_0)\).
+Constant \(C\) may take arbitrary values, here we set it to 1.
+Plugging \(C_1\) and \(C_2\) into~\eqref{eq-w-alpha-0} gives the final solution
+\begin{equation*}
+w\left(\alpha,\beta,\delta,t\right) =
+\cosh\left(iu\left(\alpha_0-\alpha\right)\right)
+\exp\left(k\delta - i\omega t\right)
+\exp\left(iv\beta\right).
+\end{equation*}
+There are two exponents in this solution with the opposite signs before
+horizontal coordinate \(\alpha\). These exponents denote
+incident and reflected wave respectively.
+The amplitude of the reflected wave does not decay as we go farther
+from the boundary, but decay only when we go deeper in the ocean. This behaviour
+corresponds to the real-world ocean waves.
+
+\subsubsection{Infinite plate}
+
+On a horizontal surface the boundary condition is written as
+\begin{equation*}
+\frac{d}{d t} \vec\nabla w \cdot \vec{n} =
+\frac{d}{d t} \frac{\partial}{\partial\delta} w = 0;
+\qquad
+\delta = \delta_0;
+\qquad
+\vec{n} = \VectorR{0\\0\\1}.
+\end{equation*}
+Analogously to wave direction we write vector form of the incident particle
+trajectory radius damping coefficient as \((0,0,k)\), hence vector form of the
+reflected coefficient is \((0,0,-k)\). We seek solutions of the
+form
+\begin{equation}
+\label{eq-w-delta-0}
+w\left(\alpha,\beta,\delta,t\right) =
+\left[ C_1\exp\left(k\delta\right) + C_2\exp\left(-k\delta\right) \right]
+\exp\left(- i\omega t\right) \exp\left(iu\alpha + iv\beta\right).
+\end{equation}
+We plug this expression into the boundary condition and get
+\begin{equation*}
+C_1\exp\left(k\delta_0\right) - C_2\exp\left(-k\delta_0\right) = 0.
+\end{equation*}
+Hence \(C_1=C\exp\left(-k\delta_0\right)\) and
+\(C_2=C\exp\left(k\delta_0\right)\). Constant \(C\) may take arbitrary values,
+here we set it to \(1/2\). Plugging \(C_1\) and \(C_2\)
+into~\eqref{eq-w-delta-0} gives the final solution
+\begin{equation*}
+w\left(\alpha,\beta,\delta,t\right) =
+\cosh\left(k\left(\delta-\delta_0\right)\right)
+\exp\left(-i\omega t\right)
+\exp\left(iu\alpha + iv\beta\right).
+\end{equation*}
+There are two exponents in this solution with opposite signs before
+vertical coordinate \(\delta\). These exponents make the radius of
+the particle trajectories decay exponentially while approaching the boundary
+\(\delta_0\) (i.e.~with increasing water depth). This is known
+solution from linear wave theory.
+
+\subsubsection{Infinite panel}
+
+On an arbitrary aligned infinite surface the boundary condition is written as
+\begin{equation*}
+\frac{d}{d t} \vec\nabla w \cdot \vec{n} = 0;
+\qquad
+\vec{n} \cdot \left(\vec\zeta - \vec\zeta_0\right) = 0,
+\end{equation*}
+where \(\vec{\zeta_0}\) is the point on the boundary plane and the third
+component of the normal vector is nought: \(\vec{n}=(n_1,n_2,0)\),
+\(\left|\vec{n}\right|=1\).
+The direction of incident wave is \(\vec{d}_i=(u,v,0)\) and the
+direction of reflected wave is \(\vec{d}_r\). We seek solutions of the form
+\begin{equation}
+\label{eq-w-all-0}
+\begin{aligned}
+w\left(\alpha,\beta,\delta,t\right) = &
+C_1\exp\left( \left(i\vec{d}_i+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega_1
+t\right) \\ + &
+C_2\exp\left( \left(i\vec{d}_r+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega_2 t\right).
+\end{aligned}
+\end{equation}
+We plug this expression into the boundary condition and get
+\begin{equation*}
+\left( i\vec{d_i}\cdot\vec{n} \right) C_1
++
+\left( i\vec{d_r}\cdot\vec{n} \right) C_2
+\exp\left( -i\vec{d}_s\cdot\vec\zeta_0 \right) = 0.
+\end{equation*}
+Here we substitute \(\vec{d}_r\cdot\vec{n}\) with \(-\vec{d}_i\cdot\vec{n}\)
+which is derived from the formula for \(\vec{d}_r\)
+(see~sec.~\ref{sec-formulae}).
+Hence, the boundary condition reduces to
+\begin{equation*}
+C_1 - C_2 \exp\left( -i\vec{d}_s\cdot\vec\zeta_0 \right) = 0.
+\end{equation*}
+Hence
+\(C_1=\frac{1}{2}\exp\left(-\frac{1}{2}i\vec{d}_s\cdot\vec\zeta_0\right)\)
+and \(C_2=\frac{1}{2}\exp\left(\frac{1}{2}i\vec{d}_s\cdot\vec\zeta_0\right)\).
+This solution reduces to the solution for the wall when \(\vec{n}=(0,0,1)\).
+
+%\input{stationary-surface.tex}
+%\input{progressively-moving-surface.tex}
+
+
+\subsection{OpenCL implementation}
Virtual testbed is a program for personal computers.
Its main feature is to perform all calculations in real time,
diff --git a/progressively-moving-surface.tex b/progressively-moving-surface.tex
@@ -0,0 +1,156 @@
+\subsubsection{Progressively moving surface}
+
+When a particle touches a parametric surface given by
+\begin{equation*}
+\vec{S}=\vec{S}(a,b,t);
+\qquad
+a,b\in{}A=[0,1];
+\qquad
+\vec{n}=\frac{\partial\vec S}{\partial a} \times \frac{\partial\vec S}{\partial b}
+\end{equation*}
+where \(a\) and \(b\) are parameters, boundary condition is written as
+\begin{equation*}
+\frac{d}{d t} \vec\nabla w \cdot\vec{n} = \frac{d}{d t} \vec S \cdot\vec{n};
+\qquad
+\vec\zeta = \vec S
+.
+\end{equation*}
+We seek solutions of the form
+\begin{equation*}
+\begin{aligned}
+& w\left(\alpha,\beta,\delta,t\right)
+= C_1\exp\left( \left(i\vec{d}_i+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega t\right) + \\
+&
+\iint\limits_{a,b\,\in{}A}
+C_2\exp\left( \left(i\vec{d}_r+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega t\right)
+da\,db\,+ \\
+&
+\iint\limits_{a,b\,\in{}A}
+C_3\exp\left( \left(i\vec{d}_i+i\vec{d}_\upsilon+\vec{d}_{k_3}\right)\cdot\vec\zeta - i\omega_3 t\right)
+da\,db,
+\end{aligned}
+\end{equation*}
+where \(\vec{d}_\upsilon=\vec{n}\left(\DerivativeT{\vec{S}}\cdot\vec{n}\right)\)
+is a projection of surface velocity to the surface normal.
+Here \(C_2\), \(C_3\) and \(\vec{n}\) depend on both \(a\) and \(b\),
+i.e.~we have one reflection for each point of the surface.
+We plug this expression into the boundary condition and get
+\begin{equation*}
+\begin{aligned}
+&
+C_1 f_1
+\exp\left( \left(i\vec{d}_i+\vec{d}_{k}\right)\cdot\vec{S}-i\omega{}t \right) + \\
+&
+\iint\limits_{a,b\,\in{}A}
+C_2 f_2
+\exp\left( \left(i\vec{d}_r+\vec{d}_{k}\right) \cdot \vec{S} -i\omega{}t \right)
+da\,db\,
++ \\
+&
+\iint\limits_{a,b\,\in{}A}
+C_3 f_3
+\exp\left(
+\left(i\vec{d}_i+i\vec{d}_\upsilon+\vec{d}_{k_3}\right)\cdot\vec{S}-i\omega_3{}t\right)
+da\,db
+= \frac{\partial\vec{S}}{\partial t} \cdot \vec{n},
+\end{aligned}
+\end{equation*}
+where
+\begin{equation*}
+\begin{aligned}
+&
+f_1 = \omega \left(\vec{d}_i\cdot\vec{n}\right);
+\qquad
+f_2 = -\omega \left(\vec{d}_i\cdot\vec{n}\right); \\
+&
+f_3 = \left(
+ i\DerivativeT{\vec{d}_\upsilon} +
+ \left(\vec{d}_i + \vec{d}_\upsilon\right) \circ
+ \left(\omega_3-\DerivativeT{\vec{d}_\upsilon}\right)
+\right)\cdot\vec{n}.
+\end{aligned}
+\end{equation*}
+Here \(\circ\) denotes Hadamard (element-wise) vector product.
+Take double integral over \(a\) and \(b\) of each side of the equation,
+assume that expressions under integral sign are equal (i.e.~remove
+integral sign from the equation), and make substitutions for derivatives that
+were made in the previous section. Then, the equation reduces to
+the form, similar to moving and rotating panel, but with different terms for the
+boundary:
+\begin{equation*}
+\begin{aligned}
+&
+C_1 f_1
+\exp\left( \left(i\vec{d}_i+\vec{d}_{k}\right)\cdot\vec{S}\right) +
+S C_2 f_2
+\exp\left( \left(i\vec{d}_r+\vec{d}_{k}\right)\cdot\vec{S}\right)
++ \\
+&
+S C_3 f_3
+\exp\left(
+\left(i\vec{d}_i+i\vec{d}_\upsilon+\vec{d}_{k_3}\right)\cdot\vec{S}-i\omega_3{}t+i\omega{}t\right)
+=
+\left(\DerivativeT{\vec{S}}\cdot\vec{n}\right)\exp\left(i\omega{}t\right),
+\end{aligned}
+\end{equation*}
+Here \(S\) without vector arrow denotes surface area. Choose \(C_1\) and \(C_2\)
+so that the first and the second term cancel each other out, and calculate
+\(C_3\) from the resulting equation:
+\begin{equation*}
+\begin{aligned}
+&
+C_1 = \frac{1}{f_1}
+\exp\left(
+ -\left(i\vec{d}_i+\vec{d}_{k}\right)\cdot\vec{S}
+\right); \\
+&
+C_2 = \frac{1}{S f_2}
+\exp\left(
+ -\left(i\vec{d}_r+\vec{d}_{k}\right)\cdot\vec{S}
+\right); \\
+&
+C_3 = \frac{\DerivativeT{\vec{S}}\cdot\vec{n}}{ S f_3 }
+\exp\left(-
+\left(i\vec{d}_i+i\vec{d}_\upsilon+\vec{d}_{k_3}\right)\cdot\vec{S}+i\omega_3{}t\right)
+\end{aligned}
+\end{equation*}
+
+Plug the solution into continuity equation and get
+\begin{equation*}
+\begin{aligned}
+&
+- \left(u - n_1 \left(\vec{d}_i\cdot\vec{n}\right)\right)^2
+- \left(v - n_2 \left(\vec{d}_i\cdot\vec{n}\right)\right)^2
++ k^2 = 0 \\
+&
+- \left(u + n_1 \left(\vec{d}_\upsilon\cdot\vec{n}\right)\right)^2
+- \left(v + n_2 \left(\vec{d}_\upsilon\cdot\vec{n}\right)\right)^2
++ k_3^2 = 0 \\
+&
+k = \pm\sqrt{u^2+v^2} \\
+&
+k_3 = \pm\sqrt{\smash[b]{u^2+v^2 +
+\left(\vec{d}_\upsilon\cdot\vec{n}\right) \left(
+\vec{d}_\upsilon\cdot\vec{n} + \vec{d}_i\cdot\vec{n} \right)}}
+\end{aligned}
+\end{equation*}
+This equation makes \(k_3\) depend on time, which does not cause trouble,
+because time derivative is always multiplied by \(n\) and the third component
+of \(n\) is always nought.
+
+When \(\vec{n}\) depends on time, i.e.~the surface rotates, the
+coefficients we write coefficients as
+\begin{equation*}
+\begin{aligned}
+&
+f_1 = \left(\omega\vec{d}_i\right)\cdot\vec{n}
+\\
+&
+f_2 =
+\left(i\vec{d}_r\circ\left(\DerivativeT{\vec{d}_s}-i\omega\right)\right)\cdot\vec{n};
+\\
+&
+f_3 =
+\end{aligned}
+\end{equation*}
+
diff --git a/references.bib b/references.bib
@@ -0,0 +1,30 @@
+
+@Article{ wackers2011rans,
+ author = {Wackers, J. and Koren, B. and Raven, H. C. and van der
+ Ploeg, A. and Starke, A. R. and Deng, G. B. and Queutey, P.
+ and Visonneau, M. and Hino, T. and Ohashi, K.},
+ title = {Free-Surface Viscous Flow Solution Methods for Ship
+ Hydrodynamics},
+ journal = {Archives of Computational Methods in Engineering},
+ year = {2011},
+ month = {Mar},
+ day = {01},
+ volume = {18},
+ number = {1},
+ pages = {1--41},
+ issn = {1886-1784},
+ doi = {10.1007/s11831-011-9059-4}
+}
+
+@Article{ nouguier2015,
+ title = {Second-order Lagrangian description of tri-dimensional
+ gravity wave interactions},
+ volume = {772},
+ doi = {10.1017/jfm.2015.179},
+ journal = {Journal of Fluid Mechanics},
+ publisher = {Cambridge University Press},
+ author = {Nouguier, Frédéric and Chapron, Bertrand and Guérin,
+ Charles-Antoine},
+ year = {2015},
+ pages = {165--196}
+}
diff --git a/stationary-surface.tex b/stationary-surface.tex
@@ -0,0 +1,62 @@
+\subsubsection{Stationary surface}
+
+When a particle touches a surface given by
+\begin{equation*}
+\vec{S}=\vec{S}(a,b,t);
+\qquad
+a,b\in{}A=[0,1];
+\qquad
+\vec{n}=\frac{\partial\vec S}{\partial a} \times \frac{\partial\vec S}{\partial b}
+\end{equation*}
+where \(a\) and \(b\) are parameters, boundary condition is written as
+\begin{equation*}
+\frac{d}{d t} \vec\nabla w \cdot\vec{n} = 0;
+\qquad
+\vec\zeta = \vec S
+.
+\end{equation*}
+The third component of the normal vector is nought: \(\vec{n}=(n_1,n_2,0)\),
+\(\left|\vec{n}\right|=1\).
+The direction of incident wave is \(\vec{d}_i=(u,v,0)\) and the
+direction of reflected wave is \(\vec{d}_r\) from
+equation~\eqref{eq-dr}. We seek solutions of the form
+\begin{equation}
+\label{eq-w-all-0}
+\begin{aligned}
+w\left(\alpha,\beta,\delta,t\right) =
+\, &
+C_1\exp\left( \left(i\vec{d}_i+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega
+t\right) \\ +
+\, &
+C_2\exp\left( \left(i\vec{d}_r+\vec{d}_{k}\right)\cdot\vec\zeta - i\omega t\right).
+\end{aligned}
+\end{equation}
+Here \(k\) is damping coefficient and \(\omega\) is wave velocity.
+The coefficient is calculated from continuity equation.
+We plug this expression into the boundary condition and get
+\begin{equation*}
+\left(\vec{d_i}\cdot\vec{n} \right) C_1
++
+\left(\vec{d_r}\cdot\vec{n} \right) C_2
+\exp\left( -i\vec{d}_s\cdot\vec\zeta_0 \right) = 0.
+\end{equation*}
+Hence
+\begin{equation*}
+C_1 =
+\frac{1}{2}
+\exp\left( -\frac{1}{2}i\vec{d}_s\cdot\zeta_0 \right)
+\qquad
+C_2 =
+\frac{1}{2}
+\exp\left( \frac{1}{2}i\vec{d}_s\cdot\zeta_0 \right)
+\end{equation*}
+Here we substituted \(\vec{d}_r\cdot\vec{n}\) with \(-\vec{d}_i\cdot\vec{n}\)
+(see~sec.~\ref{sec:org653b1bb}).
+Now plug the solution into continuity equation to get damping coefficients
+and velocity:
+\begin{equation*}
+k = \pm\sqrt{u^2+v^2};
+\qquad
+\omega^2 = g k.
+\end{equation*}
+
