1.1 The Laplace Domain

The Laplace equation is derived from the equation of continuity provided that (1) the fluid is incompressible and (2) inviscid, and (3) the flow is irrotational. Incompressibility refers to the property of the fluid to maintain its density constant within an infinitesimal volume that moves with flow velocity. An inviscid fluid, on the other hand, is a fluid with zero viscosity, namely the term (property) of the fluid that quantifies the resistance of the fluid to its gradual deformation by shear stresses. Finally, irrotational flow assumes that individual parcels of infinitesimal volume of an incompressible fluid cannot be caused to rotate. In other words, the fluid is assumed to be frictionless with no shear stresses applied on the mutual surfaces of adjacent parcels of the fluid.

Properties (1)–(3) are precisely the assumptions taken to formulate large-scale flows of liquid, namely water flows in vast volumes, just like in the open ocean, or in other words, in a marine environment. The Laplace domain has numerous applications in various disciplines of modern science, such as electromagnetics, acoustics, optics, electrostatics, imaging, hydrodynamics, and others. In hydrodynamics, however, the situation is complicated to a major extent by the existence of the free surface and subsequently by the associated boundary conditions that must be satisfied. In fact, the free surface, namely the boundary surface where the water meets the air, is what makes hydrodynamics difficult and challenging, but at the same time, a fascinating discipline.

The Laplace equation in hydrodynamics is associated with the so-called velocity potential, a scalar function, the gradient of which provides the velocity field as well as the pressure at any point of the volume of reference. Having defined the velocity field by the gradient of the velocity potential, the Laplace equation immediately follows from the mass conservation law (or the equation of continuity). Knowing the velocity potential, one may derive all dynamic parameters associated with the concerned hydrodynamic problem, namely, velocities, pressures, and accordingly hydrodynamic loads. All these, however, require the solution of relevant boundary value problems in which the boundary conditions on the free surface are of paramount importance for free-surface flows.

Describing the flow of a fluid through the Laplace equation is admittedly a simplification that leads inevitably to the assumption of an ideal fluid. This simplification, however, is a valid approach for large-scale flows in a marine environment, as has been proven beyond any doubt over the past decades. The reason is simple: in the open ocean the fluid (liquid) is water, and water is practically incompressible, nearly inviscid, while with a large degree of accuracy the flow is irrotational. In addition, use of potential theory (even nonlinear potential theory) allows the employment of analytical approaches that accordingly enable the derivation of robust analytical (or semianalytical) closed-form solutions. That advantage is provided by the separable solutions of the Laplace equation using the method of the separation of variables in various coordinate systems (Moon and Spencer, Reference Moon and Spencer1971), yielding expressions in the form of eigenfunction expansions.

It should be acknowledged, however, that all Newtonian fluids, such as water, are literally viscous. Therefore, for the completeness of the presentation, we derive in the following sections the more general formulations of the flow for Newtonian fluids, namely the Euler and the Navier–Stokes equations and eventually we focus on the flow of the ideal fluid. The analysis starts with the transport theorem, which is discussed in the section immediately following.

1.2 The Transport Theorem

We consider the following volume integral:

$I\left(t\right)=\underset{V\left(t\right)}{\overset{}{\int }}f\left(x,t\right)dV$(1.1)

where $f\left(x,t\right)$ is an arbitrary differentiable scalar function of position $x$ and time $t$. The integral is taken over the given volume $V\left(t\right)$, which changes with time. The surface $S$ that contains $V\left(t\right)$ is also a function of time and we assume that it has a normal velocity equal to $U_n$. Let us further take the same volume integral at the time instant $t+\Delta t.$ The change in integral $\Delta I$ is thus written as

$\Delta I=I\left(t+\Delta t\right)-I\left(t\right)=\underset{V\left(t+\Delta t\right)}{\overset{}{\int }}f\left(x,t+\Delta t\right)dV-\underset{V\left(t\right)}{\overset{}{\int }}f\left(x,t\right)dV$(1.2)

Also, we expand $f\left(x,t+\Delta t\right)$ in a Taylor series in $t$ to obtain

$f\left(x,t+\Delta t\right)=f\left(x,t\right)+\Delta t\frac{\partial f\left(x,t\right)}{\partial t}+\frac{1}{2}{\left(\Delta t\right)}^2\frac{{\partial }^{2}f\left(x,t\right)}{\partial t^2}+\dots$(1.3)

In the following, the function $f\left(x,t\right)$ is denoted simply by $f$. Retaining terms up to $O\left[{\left(\Delta t\right)}^2\right]$, (1.2) becomes

$\Delta I=\underset{V\left(t+\Delta t\right)}{\overset{}{\int }}\left(f+\Delta t\frac{\partial f}{\partial t}\right)dV-\underset{V\left(t\right)}{\overset{}{\int }}fdV$(1.4)

In order to simplify further the change of the volume integral $\Delta I$, we assume that $V\left(t+\Delta t\right)$ is composed of the volume $V\left(t\right)$ and an infinitesimal volume $\Delta V$ that is contained between the adjacent surfaces $S\left(t+\Delta t\right)$ and $S\left(t\right)$, i.e., $V\left(t+\Delta t\right)=V\left(t\right)+\Delta V\left(t\right)$. Hence, (1.4) is expanded in the following manner:

$\Delta I=\underset{V\left(t\right)+\Delta V\left(t\right)}{\overset{}{\int }}\left(f+\Delta t\frac{\partial f}{\partial t}\right)dV-\underset{V\left(t\right)}{\overset{}{\int }}fdV=\underset{\Delta V\left(t\right)}{\overset{}{\int }}fdV+\Delta t\underset{V\left(t\right)}{\overset{}{\int }}\frac{\partial f}{\partial t}dV+O\left[{\left(\Delta t\right)}^2\right]$(1.5)

and it is understood that the $O\left[{\left(\Delta t\right)}^2\right]$ terms also contain terms of $O\left[\Delta t\Delta V\left(t\right)\right]$.

To evaluate the integral over $\Delta V\left(t\right)$ we need to consider the infinitesimal distance between $S\left(t+\Delta t\right)$ and $S\left(t\right)$ that is equal to the thickness of the small region $\Delta V\left(t\right)$. This thickness is the normal component of the distance traveled by $S\left(t\right)$ in time $\Delta t$, and is equal to the product $U_n\Delta t$. Hence, it is immediately deduced that the first integral in (1.5) is proportional to $\Delta t$. In addition, given the infinitesimal thickness of volume $\Delta V\left(t\right)$, we may assume that the integrand $f$ is constant across that thin region in the direction normal to $S\left(t\right)$. Integrating in this direction only, we obtain

$\Delta I=\underset{S\left(t\right)}{\overset{}{\int }}\left(U_n\Delta t\right)fdS+\Delta t\underset{V\left(t\right)}{\overset{}{\int }}\frac{\partial f}{\partial t}dV+O\left[{\left(\Delta t\right)}^2\right]$(1.6)

Taking the limit of $\Delta I/\Delta t$ as $\Delta t\to 0$ yields the differential form of the transport theorem, namely

$\frac{dI}{dt}=\underset{S\left(t\right)}{\overset{}{\int }}f{U}_{n}dS+\underset{V\left(t\right)}{\overset{}{\int }}\frac{\partial f}{\partial t}dV$(1.7)

The surface integral in (1.7) represents the transport of the quantity $f$ out of the volume $V\left(t\right)$ as a result of the motion of the boundary.

1.3 Shear Stresses in Fluid Particles: The Eulerian Approach

Let us now assume a flow field in which the velocity of the fluid particles is defined by the velocity vector $U=\left(u,v,w\right)=\left(u_1,u_2,u_3\right)$, where $u$, $v$, and $w$ (or $u_1,u_2,u_3$, respectively) are the velocity components in the three directions of the Cartesian frame of reference $\left(x,y,z\right)=(x_1,x_2,x_3)$. The density of the fluid is assumed to be constant and equal to $\rho$, while the fluid particles are subjected to the external force vector $F$that is composed of the gravitational force and the surface stresses, or force per unit area, which act on adjacent surfaces of the fluid (Figure 1.1).

Figure 1.1 Stresses acting on a cubical fluid particle.

We note that each surface stress should be defined by both its direction and the orientation of the surface of the fluid particle on which it is acting. For the cubical parcel of Figure 1.1 in particular, a total of $3\times 3=9$ stress components must be defined. The origin of the Cartesian system was chosen to coincide with one of the corners of the cube only for display purposes. The reasoning behind the tensor notation is that ${\tau }_{ij}$ acts in the direction $i$ and is applied on the surface of constant $j$, with $i,j=x,y,z.$ Clearly, ${\tau }_{ij}$, $i\ne j$ denote actual shear stresses while ${\tau }_{ii}$ are normal stresses (often denoted by ${\sigma }_i$).

The definition of the stress components based on the cubical configuration of the fluid particles is in fact too restrictive. A more general configuration arises by assuming an infinitesimal volume of fluid (a fluid particle) in the shape of a tetrahedron with three orthogonal faces normal to the Cartesian coordinates as shown in Figure 1.2. The fourth face of the tetrahedron has an arbitrary oblique orientation. It is assumed the tetrahedron is sufficiently small and that the stresses are effectively constant along each face, and since the volume will be negligible compared with the surface area, the surface forces will dominate the body forces. To achieve an equilibrium state, the forces exerted on the four faces of the tetrahedron, which are induced by the surface stresses, must balance. By analogy with the stress theory in structural mechanics, the equilibrium state is secured by Cauchy’s law according to which there exists a Cauchy stress tensor $\tau$ that maps the normal to a surface to the traction vector acting on that surface. In particular, if $n={\left(n_x,n_y,n_z\right)}^T={\left(n_1,n_2,n_3\right)}^T$ denotes the unit normal vector on the oblique face (see Figure 1.2), Cauchy’s law implies that

$T=\cdot n$(1.8)

or, in Cartesian representation,

$T_x={\tau }_{xx}{n}_x+{\tau }_{xy}{n}_y+{\tau }_{xz}{n}_z$(1.9)

$T_y={\tau }_{yx}{n}_x+{\tau }_{yy}{n}_y+{\tau }_{yz}{n}_z$

$T_z={\tau }_{zx}{n}_x+{\tau }_{zy}{n}_y+{\tau }_{zz}{n}_z$

Figure 1.2 Stresses acting on the surfaces of a volume element in the shape of a tetrahedron with three orthogonal faces.

The stress tensor $\tau$, which is also denoted by ${\tau }_{ij}$, is written as

${\tau }_{ij}=\left[\begin{array}{ccc}{\tau }_{xx}& {\tau }_{xy}& {\tau }_{xz}\\ {\tau }_{yx}& {\tau }_{yy}& {\tau }_{yz}\\ {\tau }_{zx}& {\tau }_{zy}& {\tau }_{zz}\end{array}\right]=\left[\begin{array}{ccc}{\tau }_{11}& {\tau }_{12}& {\tau }_{13}\\ {\tau }_{21}& {\tau }_{22}& {\tau }_{23}\\ {\tau }_{31}& {\tau }_{32}& {\tau }_{33}\end{array}\right]=\left[\begin{array}{ccc}{\sigma }_1& {\tau }_{12}& {\tau }_{13}\\ {\tau }_{21}& {\sigma }_2& {\tau }_{23}\\ {\tau }_{31}& {\tau }_{32}& {\sigma }_3\end{array}\right]$(1.10)

Accordingly, the stress component in the $i$th direction, on a surface element with unit normal $n$, is given by ${\sum }_{j=1}^3{\tau }_{ij}{n}_j$, which for convenience will be denoted by the repeating indices product${\tau }_{ij}{n}_j,i=1,2,3.$

The stress tensor ${\tau }_{ij}$ is symmetric, i.e., ${\tau }_{ij}={\tau }_{ji}$. To demonstrate that statement, we assume that the cube shown in Figure 1.1 is infinitesimal, which allows assuming that the surface moments dominate the body moments. In addition, changes in the magnitude of the stresses may be neglected across the cube and a positive value of ${\tau }_{zx}$ (Figure 1.1) will cause a counterclockwise moment about the centroid. The moment on the opposite face will have the same sign, since the normal vector has the opposite sense and hence a positive ${\tau }_{zx}$ acts downward. Accordingly, the counterclockwise moment is balanced by a clockwise moment that is developed only by the stress ${\tau }_{xz}$ that acts in the $x$-direction on the top and the bottom faces of the cube. Hence, it follows that ${\tau }_{zx}={\tau }_{xz}$ and accordingly ${\tau }_{xy}={\tau }_{yx}$ and ${\tau }_{yz}={\tau }_{zy}$.

1.4 Mass Conservation and Momentum Conservation

Let us now assume that $V\left(t\right)$, defined originally in Section 1.2, is a material volume of fluid that is assembled by a group of material fluid particles. Accordingly, the total mass of the fluid in $V\left(t\right)$ is given by the volume integral of density $\rho$ in $V\left(t\right)$, while the mass conservation law requires that the total mass must be constant and independent of time $t$. Hence

$\frac{d}{dt}\underset{V\left(t\right)}{\overset{}{\int }}\rho dV=0$(1.11)

Further, the momentum conservation law is written as

$\frac{d}{dt}\underset{V\left(t\right)}{\overset{}{\int }}\rho u_{i}dV=\underset{S\left(t\right)}{\overset{}{\int }}{\tau }_{ij}{n}_j\text{}dS+\underset{V\left(t\right)}{\overset{}{\int }}{F}_{i}dV$(1.12)

which implies that the sum of all forces exerted on a volume of fluid must be equal to the rate of change of its momentum with respect to the Newtonian frame of reference. The momentum rate of change on the left-hand side of (1.12) refers to the $i$th component of the velocity, the surface integral on the right-hand side denotes the sum of the stress forces in the $i$th direction, while the last term is the total external force exerted on the concerned volume of fluid, again in the $i$th direction. Clearly, $F_i$ is the $i$th component of the aforementioned vector $F$.

The momentum conservation law can be effectively transformed in order to be expressed by volume integrals only. This can be achieved using the Gauss theorem (also known as the divergence theorem). In particular, for a vector $W$ that is continuous and differentiable in the volume $V\left(t\right)$ and the unit normal $n$ that is fixed on the surface $S$ surrounding $V\left(t\right)$ and pointing to the exterior of $V\left(t\right)$, the Gauss theorem reads in vectorial and tensorial forms:

$\underset{V\left(t\right)}{\overset{}{\int }}\nabla \cdot WdV=\underset{S\left(t\right)}{\overset{}{\int }}W\cdot ndS,\underset{V\left(t\right)}{\overset{}{\int }}\frac{\partial W_i}{\partial x_i}dV=\underset{S\left(t\right)}{\overset{}{\int }}{W}_i{n}_{i}dS$(1.13)

Hence, using (1.13) the momentum conservation law (1.12) becomes

$\frac{d}{dt}\underset{V\left(t\right)}{\overset{}{\int }}\rho u_{i}dV=\underset{V\left(t\right)}{\overset{}{\int }}\left(\frac{\partial {\tau }_{ij}}{\partial x_j}+F_i\right)dV$(1.14)

1.5 The Equation of Continuity

Let us now consider the special case in which $V\left(t\right)$ is a material volume always composed of the same fluid particles, while the surface $S\left(t\right)$ in which the volume is contained moves with the same normal velocity as the fluid and $U_n=U\cdot n=u_i{n}_i$. Accordingly, using (1.1), (1.7), and the Gauss theorem (1.13) it follows that

$\frac{d}{dt}\underset{V\left(t\right)}{\overset{}{\int }}fdV=\underset{V\left(t\right)}{\overset{}{\int }}\frac{\partial f}{\partial t}dV+\underset{S\left(t\right)}{\overset{}{\int }}f{u}_i{n}_{i}dV=\underset{V\left(t\right)}{\overset{}{\int }}\left[\frac{\partial f}{\partial t}+\frac{\partial \left(f{u}_i\right)}{\partial x_i}\right]dV$(1.15)

Returning to the mass conservation (1.11), (1.15) for $f=\rho$ yields

$\frac{d}{dt}\underset{V\left(t\right)}{\overset{}{\int }}\rho dV=\underset{V\left(t\right)}{\overset{}{\int }}\left[\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}\right]dV=0$(1.16)

Clearly, (1.16) should hold for an arbitrary group of fluid particles and for any instant of time $t$. Therefore the integrant itself should be zero, resulting in

$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0$(1.17)

Hence for an incompressible fluid and constant density, (1.17) yields readily the so-called equation of continuity, which reads

$\frac{\partial u_i}{\partial x_i}=0,\nabla \cdot U=0$(1.18)

1.6 Euler Equations

If the transport theorem (1.7) is applied to the momentum conservation law (1.12), namely using $\rho u_i$ instead of $\rho$, it follows that

$\underset{V\left(t\right)}{\overset{}{\int }}\left[\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}\right]dV=\underset{V\left(t\right)}{\overset{}{\int }}\left(\frac{\partial {\tau }_{ij}}{\partial x_j}+F_i\right)dV$(1.19)

or

$\underset{V\left(t\right)}{\overset{}{\int }}\left[\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}-\frac{\partial {\tau }_{ij}}{\partial x_j}-F_i\right]dV=0$(1.20)

Again, for an arbitrary group of fluid particles and for any instant of time, (1.20) is satisfied if the integrand is zero. Thus

$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=\frac{\partial {\tau }_{ij}}{\partial x_j}+F_i$(1.21)

Finally, for a constant density $\rho$ and using the equation of continuity (1.18), (1.21) yields the Euler equations, i.e.,

$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=\frac{1}{\rho }\frac{\partial {\tau }_{ij}}{\partial x_j}+\frac{1}{\rho }{F}_i$(1.22)

1.7 Stress Relations in a Newtonian Fluid

Inevitably, the stress tensor and the kinematical properties of the fluid are tightly correlated. Analogous correlations hold for the stress–strain relations in solid mechanics. Assuming that the fluid is at rest, then no shear stresses are developed between the fluid particles and only a normal pressure component (pressure stress) will exist within the fluid. With reference to Figure 1.2, it is deduced that in order to balance the forces acting across a small tetrahedron, the normal pressure must be isotropic. Given that no viscous shear was considered, the stress tensor will be composed only of the pressure stress, i.e.,

${\tau }_{ij}=-p{\delta }_{ij}$(1.23)

By hypothesis, there are no viscous forces if the fluid moves as a rigid mass without deformation, or with the velocity field

$U=A+B\times r$(1.24)

where $A$ and $B$ are constant vectors equal to the translational and rotational velocities and $r$ is the position vector from the origin of the rotation. Viscous stresses will occur when the fluid velocity differs from (1.24) with relative motion between adjacent fluid particles. The simplest example is a uniform shear flow. Here, and for more general velocity fields, the fundamental assumption of a Newtonian fluid is that the stress tensor is a linear function of the nine gradients $\partial u_k/\partial x_l$. This ensures vanishing of the viscous stress components for uniform translation of the fluid. The rotational term in (1.24) will be stress-free provided the gradients occur only in the form of sums $\partial u_k/\partial x_l+\partial u_l/\partial x_k$. For an isotropic fluid, the values of the stress components must be independent of the choice of coordinates, and for flow in one plane there can be no shear stress in the direction normal to this plane. The most general linear function of the velocity gradients, consistent with these conditions and the requirement that ${\tau }_{ij}$ is symmetric, is of the form

${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),i\ne j$(1.25)

where $\mu$ is the viscous shear coefficient, also known as the coefficient of viscosity.

In order for ${\tau }_{ij}$ be a tensor, the only possible addition to (1.25) for $i=j$ is a second constant times the divergence $\partial u_i/\partial x_i$, which vanishes for an incompressible fluid. Thus, the total stress tensor for an incompressible fluid is given by

${\tau }_{ij}=-p{\delta }_{ij}+\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$(1.26)

or in Cartesian representation

$\left[{\tau }_{ij}\right]=\left[\begin{array}{ccc}-p& 0& 0\\ 0& -p& 0\\ 0& 0& -p\end{array}\right]+\mu \left[\begin{array}{ccc}2\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}& \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\\ \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}& 2\frac{\partial v}{\partial y}& \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\\ \frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}& \frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}& 2\frac{\partial w}{\partial z}\end{array}\right]$(1.27)

The diagonal of the first matrix in (1.27) is occupied by the normal pressure stress. The second matrix is the viscous stress tensor and is proportional to the viscous shear coefficient $\mu$. The diagonal elements of the second matrix are associated with the elongations of fluid elements. The off-diagonal elements are associated with shear deformations.

1.8 The Navier–Stokes Equations

Having derived the Euler equations and expressed the stress–strain relations via (1.25) the Navier–Stokes equations are obtained by direct substitution of (1.25) into (1.22). The Navier–Stokes equations express the conservation of momentum of a Newtonian fluid. Given that the equation of continuity requires that $\frac{{\partial }^2{u}_j}{\partial x_j\partial x_j}=\frac{\partial }{\partial x_i}\frac{\partial u_j}{\partial x_j}=0$, the derivatives of the stress tensor will read

$\frac{\partial {\tau }_{ij}}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\mu \frac{\partial }{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)=-\frac{\partial p}{\partial x_i}+\mu \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}$(1.28)

Hence, the Navier–Stokes equations are written as

$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\nu \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}+\frac{1}{\rho }{F}_i$(1.29)

where $\nu =\mu /\rho$ is the coefficient of the kinematic viscosity. The Navier–Stokes equations are often written in a vectorial form according to

$\frac{\partial U}{\partial t}+\left(U\cdot \nabla \right)U=-\frac{1}{\rho }\nabla p+\nu {\nabla }^{2}U+\frac{1}{\rho }F$(1.30)

Finally, the Cartesian representation of the Navier–Stokes equations is

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu {\nabla }^{2}u+\frac{1}{\rho }{F}_x$(1.31)

$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\nu {\nabla }^{2}v+\frac{1}{\rho }{F}_y$(1.32)

$\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\nu {\nabla }^{2}w+\frac{1}{\rho }{F}_z$(1.33)

The motion of a viscous fluid with constant density and a Newtonian stress–strain relation is fully described by the system of partial differential equations (1.31)–(1.33) and the equation of continuity (1.18). Most common fluids, including water and air, appear to comply with the requirement of a Newtonian stress–strain relation. Exceptions occur for nonisotropic fluids, in which case the stress–strain relation is said to be “non-Newtonian.” Even in the former case, however, the solution of the Navier–Stokes equations raises major difficulties. Major challenges originate from the fact that those equations literally form a coupled nonlinear system that has been solved analytically only for some very simple geometrical configurations, principally those in which the nonlinear convective acceleration terms $\left(U\cdot \nabla \right)U$ can be assumed to vanish (Newman, Reference Newman1977).

1.9 Inviscid, Incompressible Fluid and Irrotational Flow: The Velocity Potential

1.9.1 Inviscid, Incompressible

Those assumptions lead to the concept of an ideal fluid and eventually to potential flow. Incompressibility is secured by the constant density of the fluid $\rho =const$. In addition, for a nonviscous fluid, the Navier–Stokes equations are reduced to

$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{1}{\rho }{F}_i$(1.34)

Further, we assume that the external force $F_i$ is represented by only the gravitational force $\rho g$ ($g$ is the gravitational acceleration), which is directed vertically downward along the $z\equiv x_3$ coordinate. Letting the $z\equiv x_3$ coordinate pointing upward, it follows that $F_i=\left(0,0,-\rho g\right)$. Hence (1.34) may be rewritten as

$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial }{\partial x_i}\left(p+\rho g{x}_3\right)$(1.35)

1.9.2 Irrotational Flow: The Velocity Potential

The correlation of the velocity potential with the flow field is realized through its gradient, which provides the velocity vector of the fluid particles. The gain from using the velocity potential is that instead of three components (namely the velocities in all directions of the three-dimensional space) only one is needed, which apparently depends on the three spatial coordinates and the time. Further, although the velocities have an obvious physical meaning, the potential is somehow a vague term that is materialized physically through its correlation with the velocities. Nevertheless, it can be proved that the potential is indeed an effective representation of the velocity field. To this end, only the assumption of the irrotational flow is required.

In particular let us consider the definite integral

$\Phi \left(x,t\right)=\underset{x_0}{\overset{x}{\int }}{u}_{i}d{x}_i$(1.36)

where the lower limit is some arbitrary constant position $x_0$ and the upper limit is the point $x=(x_1,x_2,x_3)$. This integral is independent of the particular path of integration between the points $x_0$ and $x$, since the difference in value of any two integrals, between the same two points, is equal to circulation around a closed path from $x_0$ to $x$, along one path and back to $x_0$ along the other path, which is equal to zero if the fluid motion is irrotational. Thus, the integration in (1.36) can be performed along any desired path. If we choose a path that approaches the point $x$, along a straight line parallel to the $x_1$-axis, then along the final portion of the path of integration $u_{i}d{x}_i=u_{1}d{x}_1$, so that

$\frac{\partial \Phi }{\partial x_1}=\frac{\partial }{\partial x_1}\underset{x_0}{\overset{x}{\int }}{u}_{1}d{x}_1=u_1$(1.37)

The remaining portion of the integral being a constant, does not contribute to the derivative. Applying a similar argument for the other two coordinates, we have in general that

$u_i=\frac{\partial \Phi }{\partial x_i},U=\nabla \Phi$(1.38)

The existence of the velocity potential therefore is univocally justified by the assumption of an irrotational flow.

Substituting (1.38) into the equation of continuity (1.18) immediately yields

$\frac{{\partial }^2\Phi }{\partial x_i\partial x_i}=0,{\nabla }^2\Phi =\frac{{\partial }^2\Phi }{\partial x^2}+\frac{{\partial }^2\Phi }{\partial y^2}+\frac{{\partial }^2\Phi }{\partial z^2}=0$(1.39)

Equation (1.39) is the Laplace equation given in tensorial and Cartesian representations and ${\nabla }^2={\partial }^2/\partial x^2+{\partial }^2/\partial y^2+{\partial }^2/\partial z^2$ is the Laplace operator.

1.9.3 The Unsteady Bernoulli Equation

The Laplace equation is the governing equation of the flow field under the basic assumptions made in this section. It will also be referred to as the field equation, which we recall is the outcome of the equation of continuity. To complete the problem’s setup within the realm of potential theory we must exploit (1.35), which eventually provides the pressure field that can be also univocally obtained by the potential function only. We are particularly interested for the unsteady case, including the time derivative term of (1.35). The steady case is a special case in which this term is absent. Using (1.38), (1.35) yields

$\frac{\partial }{\partial t}\frac{\partial \Phi }{\partial x_i}+\frac{\partial \Phi }{\partial x_j}\frac{\partial }{\partial x_j}\frac{\partial \Phi }{\partial x_i}=-\frac{1}{\rho }\frac{\partial }{\partial x_i}\left(p+\rho g{x}_3\right)$(1.40)

It can be easily shown that

$\frac{1}{2}\frac{\partial }{\partial x_i}\frac{\partial \Phi }{\partial x_j}\frac{\partial \Phi }{\partial x_j}=\frac{\partial \Phi }{\partial x_j}\frac{\partial }{\partial x_j}\frac{\partial \Phi }{\partial x_i}$(1.41)

Substituting (1.41) into (1.40) and rearranging terms yields

$\frac{\partial }{\partial x_i}\left[\frac{\partial \Phi }{\partial t}+\frac{1}{2}\frac{\partial \Phi }{\partial x_j}\frac{\partial \Phi }{\partial x_j}+\frac{1}{\rho }\left(p+\rho g{x}_3\right)\right]=0$(1.42)

Finally, integrating with respect to $x_i$ for $i=1,2,3$ one gets

$\frac{\partial \Phi }{\partial t}+\frac{1}{2}\frac{\partial \Phi }{\partial x_j}\frac{\partial \Phi }{\partial x_j}=-\frac{1}{\rho }\left(p+\rho g{x}_3\right)+C\left(t\right)$(1.43)

The integration with respect to $x_i$, $i=1,2,3$ will produce an integration factor that is a function of time $C\left(t\right)$. This can be incorporated into the velocity potential, writing

${\Phi }^{\prime }=\Phi -\underset{0}{\overset{t}{\int }}C\left(\xi \right)d\xi$(1.44)

Finally, using (1.44) into (1.43) yields the celebrated unsteady Bernoulli equation

$\frac{\partial {\Phi }^{\prime }}{\partial t}+\frac{1}{2}\frac{\partial {\Phi }^{\prime }}{\partial x_j}\frac{\partial {\Phi }^{\prime }}{\partial x_j}=-\frac{1}{\rho }\left(p+\rho g{x}_3\right)$(1.45)

The prime will be dropped in the sequel. In the Cartesian frame of reference, (1.45) is written as

$\frac{\partial \Phi }{\partial t}+\frac{1}{2}{\left|\nabla \Phi \right|}^2+\frac{p}{\rho }+gz=0$(1.46)

Hence, the flow field for incompressible and inviscid fluid and irrotational flow is fully described by potential flow theory via the Laplace equation (1.39) and the unsteady Bernoulli equation (1.46). Those equations describe explicitly the flow in a boundless medium. Practical hydrodynamic problems, however, involve always boundary conditions that are represented by mathematical constraints, i.e., differential or algebraic equations that must be satisfied together with the Laplace and the Bernoulli equations. Even the conditions at infinity can be regarded as boundary conditions. Thus, in hydrodynamics we are invited to consider boundary value problems formed by several equations accompanying the fundamental equations (1.39) and (1.46).

1.10 Free-Surface Flow in the Laplace Domain

The hydrodynamical applications associated with the marine environment are in most of the cases fluid–structure interaction problems. In relevant situations, the factor that makes the difference is the existence of the free surface, i.e., the mutual surface that separates water and air. The effect of the free surface is of paramount importance in wave–structure interaction problems and also in hydrodynamical problems without incident waves. As an example, which will be analyzed extensively in Chapter 5, we refer to the case of a solid moving under an undisturbed free surface with no waves, namely without periodic or nonperiodic deformation of the free surface. This is the so-called “wave resistance” problem. Clearly, the motion of the solid affects the kinematical condition of the liquid up to the free surface especially in cases where the solid is in close proximity with the free surface. It is evident therefore that the free surface will affect the flow field, i.e., the velocities of the liquid particles and the pressure field exerted on the solid.

There are two main challenges associated with the existence of the free surface: (1) it imposes constraints, or in other words boundary conditions, which are in fact nonlinear and (2) the boundary where those conditions should be satisfied is time dependent and is one of the unknowns to be determined. Both issues are discussed in the present section.

In marine hydrodynamics the fluid is liquid. The domain of interest is the material volume of liquid below the free surface and the alterations of its kinematical condition(s) due to the kinematics of both the liquid (e.g., water waves, current) and the body(ies) that exist within the liquid. We start the analysis with the most typical hydrodynamical problem of a body with arbitrary geometry that floats on the free surface of a liquid field that extends to infinity both in $x$ and $y$ directions as shown in Figure 1.3 (in several cases it is most convenient to assume that the surface at infinity $S_{\infty }$ is a cylindrical one, situated at $r\to \infty$, where $r$ denotes the radial coordinate of a polar frame of reference). The free surface itself is a function of time and apparently of the horizontal plane coordinates $(x,y)$. Therefore a fixed and explicitly defined plane of reference is required to define the free-surface elevation $H$ as a function of $(x,y,t)$. Given the fact that in the general case the topography of the bottom, which is also a fixed surface, can be arbitrary, the only alternative is to assume as a plane of reference the undisturbed free surface and define it on $z=0$. It is evident that the schematic in Figure 1.3 assumes a plane horizontal bottom that is located at a distance $h$ below the undisturbed free surface. The depth $h$ can be finite or infinite.

The control volume $\Omega$ in Figure 1.3 is bounded from above by the free surface, coined $FS$, represented by the function $H(x,y,t)$, from below by the bottom surface, coined $S_{BT}$, and by the impermeable surface of the body $S_0$. The basic condition is to assume infinite extend for $(x,y)\in (-\infty ,+\infty )$, while one can add more complications, e.g., by restricting the extent of the control volume, which is achieved by introducing additional boundary conditions, or putting more solids into the liquid field, floating, fully immersed or bottom seated.

Figure 1.3 A solid that floats in a liquid field of infinite extent.

According to the Laplace domain approach, namely inviscid, incompressible liquid and irrotational flow, the velocity field, described by the velocity components of the liquid particles $(u,v,w)$, in the three directions $(x,y,z)$ respectively, is expressed in terms of the scalar velocity potential $\Phi \left(x,y,z,t\right)$ as

$\left(u,v,w\right)=\nabla \Phi \left(x,y,z,t\right)=\frac{\partial u}{\partial x}\overrightarrow{i}+\frac{\partial v}{\partial y}\overrightarrow{j}+\frac{\partial w}{\partial z}\overrightarrow{k}$(1.47)

where $\left(\overrightarrow{i},\overrightarrow{j},\overrightarrow{k}\right)$ are the unit vectors in $(x,y,z)$.

The basic hydrodynamic boundary value problem is thus defined by the following set of equations:

${\nabla }^2\Phi =\frac{{\partial }^2\Phi }{\partial x^2}+\frac{{\partial }^2\Phi }{\partial y^2}+\frac{{\partial }^2\Phi }{\partial z^2}=0,(x,y,z)\in \Omega$(1.48)

$\nabla \Phi \cdot n=v\cdot n+\cdot \left(x\times n\right),(x,y,z)\in S_0$(1.49)

$\nabla \Phi \cdot n=0,\left(x,y,z\right)\in S_{BT}$(1.50)

$\frac{{\partial }^2\Phi }{\partial t^2}+g\frac{\partial \Phi }{\partial z}+2\nabla \Phi \cdot \nabla \left(\frac{\partial \Phi }{\partial t}\right)+\frac{1}{2}\nabla \Phi \cdot \nabla \left(\nabla \Phi \cdot \nabla \Phi \right)=0,z={\rm H}(x,y,t)$(1.51)

${\rm H}=-\frac{1}{g}\left(\frac{\partial \Phi }{\partial t}+\frac{1}{2}\nabla \Phi \cdot \nabla \Phi \right),z={\rm H}(x,y,t)$(1.52)

Equation (1.48) is the Laplace equation, (1.49) is the boundary condition that must be satisfied on the surface of the solid, (1.50) is the kinematic boundary condition on the bottom, (1.51) denotes the combined dynamic and kinematic boundary condition of the free surface, and finally (1.52) provides the exact position of the free surface. In (1.50) $x$ is the position vector of a point on $S_0$ measured from the center of rotation, while $v=(v_1,v_2,v_3)$ and $=(v_4,v_5,v_6)$ are the translational and rotational vectors of the solid in surge, sway, heave, roll, pitch, and yaw respectively. Clearly $n$ represents the unit normal vector upon the associated impermeable surface that is pointing out of the fluid domain.

More detailed discussion is required for the free-surface condition (1.51). The above boundary value problem should be completed by an appropriate radiation condition at infinity, i.e., as $r=\sqrt{x^2+y^2}\to \infty$.

1.11 Free-Surface Kinematics

The conditions associated with the presence of the free surface are the so-called dynamic and kinematic boundary conditions. Omitting the effect of the surface tension, the former condition implies that the pressure on the free surface should be continuous. The pressure at any point in the liquid is obtained through the Bernoulli equation (1.46) in which $p$ will denote the pressure of the liquid, relative to the atmospheric pressure. As the density of the air is very small compared to the density of the liquid, the motion of the air can be neglected, implying no dynamic fluctuations of the pressure on the free surface. Therefore, $p$ can be taken equal to zero and accordingly the dynamic condition on the free surface, evaluated on $z=H(x,y,t)$, becomes

$\frac{\partial \Phi }{\partial t}+\frac{1}{2}{\left|\nabla \Phi \right|}^2+g{\rm H}=0,z={\rm H}(x,y,t)$(1.53)

Also, continuity of the interface (the surface between the air and the liquid) requires that the liquid particles on the interface remain on it during the motion of the free surface. That is expressed mathematically by

$\frac{D({\rm H}-z)}{Dt}=0,\frac{\partial {\rm H}}{\partial t}+\nabla \Phi \cdot \nabla {\rm H}=\frac{\partial \Phi }{\partial z},z={\rm H}(x,y,t)$(1.54)

where $D/Dt$ denotes the material derivative. Next, we take the gradient and the time derivative of (1.53) to yield

$\nabla \left(\frac{\partial \Phi }{\partial t}\right)+\frac{1}{2}\nabla \cdot \left(\nabla \Phi \cdot \nabla \Phi \right)+g\nabla H=0,z=H(x,y,t)$(1.55)

and

$\frac{{\partial }^2\Phi }{\partial t^2}+\nabla \Phi \cdot \nabla \left(\frac{\partial \Phi }{\partial t}\right)+g\frac{\partial H}{\partial t}=0,z=H(x,y,t)$(1.56)

Substituting $\nabla H$ from (1.55) into (1.54) and the resulting expression into (1.56) establishes (1.51). The exact position of the free surface is obtained directly from (1.53), which yields (1.52). Equation (1.51) is the fundamental, combined dynamic and kinematic free surface boundary condition. The fact that it has been expressed solely in terms of the velocity potential allows some simplification of the intended hydrodynamic problems, which that way, could be considered only in terms of the potential.

1.12 The Taylor Expansion of the Free Surface

The main problem associated with the free surface is that its conditions should be valid in a time-varying boundary, namely on $z=H(x,y,t)$. To overcome this difficulty, Stokes suggested the expansion of the boundary conditions (assuming small displacements of the free surface) around a mean position that coincides with the undisturbed free surface at $z=0$. Therefore, any function $f\left(z\right)$ evaluated on $z=H(x,y,t)$ is written as

${f\left(z\right)|}_{z=H}={f\left(z\right)|}_{z=0}+H{\frac{\partial f\left(z\right)}{\partial z}|}_{z=0}+\frac{1}{2}{H}^2{\frac{{\partial }^{2}f\left(z\right)}{\partial z^2}|}_{z=0}+O\left(H^3\right)$(1.57)

Accordingly, (1.52) gives for the free surface

$\begin{array}{c}H=-\frac{1}{g}[\frac{\partial \Phi }{\partial t}+\frac{1}{2}\nabla \Phi \cdot \nabla \Phi -\frac{1}{g}\frac{\partial \Phi }{\partial t}\frac{{\partial }^2\Phi }{\partial t\partial z}-\frac{1}{2g}\frac{{\partial }^2\Phi }{\partial t\partial z}\nabla \Phi \cdot \nabla \Phi -\frac{1}{\text{g}}\frac{\partial \Phi }{\partial t}\nabla \Phi \cdot \nabla \left(\frac{\partial \Phi }{\partial z}\right)\\ +\frac{1}{2g^2}\frac{\partial \Phi }{\partial t}\frac{\partial \Phi }{\partial t}\frac{{\partial }^3\Phi }{\partial t\partial z^2}]+O\left({\Phi }^4\right)\end{array}$(1.58)

while the boundary condition (1.51) becomes

$\begin{array}{c}\frac{{\partial }^2\Phi }{\partial t^2}+g\frac{\partial \Phi }{\partial z}=-2\nabla \Phi \cdot \nabla \left(\frac{\partial \Phi }{\partial t}\right)-\frac{1}{2}\nabla \Phi \cdot \nabla \left(\nabla \Phi \cdot \nabla \Phi \right)+\frac{1}{g}\frac{\partial \Phi }{\partial t}\left(\frac{{\partial }^3\Phi }{\partial t^2\partial z}+g\frac{{\partial }^2\Phi }{\partial z^2}\right)\\ +\frac{2}{g}\frac{\partial \Phi }{\partial t}\left[\nabla \left(\frac{\partial \Phi }{\partial z}\right)\cdot \nabla \left(\frac{\partial \Phi }{\partial t}\right)+\nabla \Phi \cdot \nabla \left(\frac{{\partial }^2\Phi }{\partial t\partial z}\right)\right]\\ -\frac{1}{g}\left(\frac{1}{g}\frac{\partial \Phi }{\partial t}\frac{{\partial }^2\Phi }{\partial t\partial z}-\frac{1}{2}\nabla \Phi \cdot \nabla \Phi \right)\left(\frac{{\partial }^3\Phi }{\partial t^2\partial z}+g\frac{{\partial }^2\Phi }{\partial z^2}\right)\\ -\frac{1}{2g^2}{\left(\frac{\partial \Phi }{\partial t}\right)}^2\left(\frac{{\partial }^4\Phi }{\partial t^2\partial z^2}+g\frac{{\partial }^3\Phi }{\partial z^3}\right)+O\left({\Phi }^4\right)\end{array}$(1.59)

Both (1.58) and (1.59) have been expressed on the known boundary of the undisturbed free surface on $z=0$.

1.13 Expansion in Perturbations

Although the Taylor series expansion tackles the problem of the time-varying boundary in which the free-surface conditions should apply, it is unable to simplify these conditions, which remain strongly nonlinear. Admittedly, processing of nonlinear systems is a very challenging task. Consideration of (1.59) in its complete form is feasible only in the time domain. The frequency domain approximation requires further simplifications, which in hydrodynamics (and other scientific disciplines) are employed using expansions in series of perturbations (e.g., Kevorkian and Cole, Reference Kevorkian and Cole1981). The expansion of the velocity potential (and the other hydrodynamic components) in perturbation series was first employed by Stokes and accordingly the perturbation scheme is widely refereed as Stokes perturbations corresponding to the sequence of Stokes waves. Use of expansions in perturbations requires a small parameter (the scaling factor) that is the problem’s property and preferentially known. In free-surface hydrodynamics, this role is assumed by the wave steepness $ϵ=k_{0}A\ll 1$, where $k_0$ denotes the wavenumber (being equal to $2\pi /\lambda$) of the regular waves propagating in a liquid field of specified water depth ($\lambda$ is the wave length) and $A$ denotes the linear amplitude of the regular waves. Thus, using the scaling factor $ϵ$, the velocity potential is expressed as an infinite sum according to

$\Phi (x,y,z,t)=\sum _{n=1}^{\infty }{ϵ}^n{\Phi }^{\left(n\right)}(x,y,z,t)$(1.60)

There is a direct correlation between the wavenumber $k_0$ and the (circular) frequency $\omega$ of the propagating waves expressed through the so-called “dispersion” relation, which will be discussed in the sequel. This correlation involves the water depth as well, which modifies the wave length for the same wave frequency. More details will be given in the sequel. Therefore, the actual reasoning behind the expansion (1.60) in terms of the wave steepness is that the total velocity potential involves terms (powers of $ϵ$) depending on the pairs $(A,\omega )$, $(A^2,2\omega$), $(A^3,3\omega )$, and so on.

Further, assuming that at the leading (first-) order, the potential varies periodically with circular frequency $\omega$, we can write

$ϵ{\Phi }^{\left(1\right)}(x,y,z,t)=\text{Re}\left[{\varphi }^{\left(1\right)}\left(x,y,z\right)e^{-i\omega t}\right]$(1.61)

where $\text{Re}$ denotes the real part of a complex argument. The nonlinear products of the inhomogeneous term in (1.59) suggest taking the following forms for the second- and the third-order potentials (phenomena higher than third-order will not be considered):

${ϵ}^2{\Phi }^{\left(2\right)}(x,y,z,t)={\overline{\varphi }}^{\left(2\right)}(x,y,z)+\text{Re}\left[{\varphi }^{\left(2\right)}\left(x,y,z\right)e^{-2i\omega t}\right]$(1.62)

${ϵ}^3{\Phi }^{\left(3\right)}(x,y,z,t)=\text{Re}\left[{\overline{\varphi }}^{\left(3\right)}\left(x,y,z\right)e^{-i\omega t}\right]+\text{Re}\left[{\varphi }^{\left(3\right)}\left(x,y,z\right)e^{-3i\omega t}\right]$(1.63)

Working the same way, the free-surface elevation is expanded in a series of perturbations according to

$\begin{array}{c}H\left(x,y,t\right)=\sum _{n=1}^{\infty }{ϵ}^n{H}^{\left(n\right)}\left(x,y,t\right)=\text{Re}\left[{\eta }^{\left(1\right)}\left(x,y\right)e^{-i\omega t}\right]+{\overline{\eta }}^{\left(2\right)}\left(x,y\right)\\ +\text{Re}\left[{\eta }^{\left(2\right)}\left(x,y\right)e^{-2i\omega t}\right]+\text{Re}\left[{\overline{\eta }}^{\left(3\right)}\left(x,y\right)e^{-i\omega t}\right]+\text{Re}\left[{\eta }^{\left(3\right)}\left(x,y\right)e^{-3i\omega t}\right]+O\left({ϵ}^4\right)\end{array}$(1.64)

Clearly, the second-order (double frequency) problem involves a steady component, while the third-order problem embraces also a term that varies with the leading order frequency $\omega$.

The next step is to introduce (1.61)–(1.63) into the free-surface boundary condition (1.59), and (1.64) into the free-surface elevation (1.58). Doing so and equating like powers of $ϵ$ (up to ${ϵ}^3$) will result in the following boundary conditions and free-surface elevations at the various orders.

• Leading (first-) order problem $O\left(ϵ\right)$

$-K{\varphi }^{\left(1\right)}+\frac{\partial {\varphi }^{\left(1\right)}}{\partial z}=0$(1.65)

${\eta }^{\left(1\right)}=\frac{i\omega }{g}{\varphi }^{\left(1\right)}$(1.66)

• Second-order problem $O\left({ϵ}^2\right)$

$-4K{\varphi }^{\left(2\right)}+\frac{\partial {\varphi }^{\left(2\right)}}{\partial z}=\frac{i\omega }{g}\left[\nabla {\varphi }^{\left(1\right)}\cdot \nabla {\varphi }^{\left(1\right)}-\frac{1}{2}{\varphi }^{\left(1\right)}\left(\frac{{\partial }^2{\varphi }^{\left(1\right)}}{\partial z^2}-K\frac{\partial {\varphi }^{\left(1\right)}}{\partial z}\right)\right]$(1.67)

${\eta }^{\left(2\right)}=\frac{2i\omega }{g}{\varphi }^{\left(2\right)}-\frac{1}{4g}\nabla {\varphi }^{\left(1\right)}\cdot \nabla {\varphi }^{\left(1\right)}-\frac{K^2}{2g}{\varphi }^{\left(1\right)}{\varphi }^{\left(1\right)}$(1.68)

${\overline{\eta }}^{\left(2\right)}=-\frac{1}{4g}\nabla {\varphi }^{\left(1\right)}\cdot \nabla {\varphi }^{\left(1\right)*}+\frac{K^2}{2g}{\varphi }^{\left(1\right)}{\varphi }^{\left(1\right)*}$(1.69)

where the asterisk denotes the equivalent complex conjugate (cc) of the associate component.

• Third-order problem $O\left({ϵ}^3\right)$