*2.1. Governing Equations*

In order to represent the payload's posture in the wave, two Cartesian frames are defined, as shown in Figure 1a. The world frame (*ow* − **x***w***y***w***z***w*) defined based on the 3D NWT. *ow* is the midpoint of the inlet. *ow***x***w* is the direction of wave propagation. *ow***z***w* points straight upwards. The body frame is fixed with the payload. As for the body frame (*ob* − **x***b***y***b***z***b*) of the cuboid payload, *ob* is the centroid of the payload, and the three axes follow the directions of the three edges of the cuboid, respectively. The cuboid payload's posture in the wave can be expressed by the three Euler angles θ, φ,ψ (pitch, roll, and yaw angles), which represent the pose relationship between the payload's body frame (*ob* − **x***b***y***b***z***b*) and the world frame (*ow* − **x***w***y***w***z***w*) as shown in Figure 1b.

*J. Mar. Sci. Eng.* **2020**, *8*, 433

Both air and water are assumed to be incompressible laminar fluid. The motion of the fluid continuum is described with the governing equations, i.e., the Navier–Stokes equations and the continuity equation [16],

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{U}) = 0 \tag{1}$$

$$\frac{d\rho\mathbf{U}}{dt} + \nabla \times (\rho\mathbf{U}\mathbf{U}) - \nabla \times (\mu\nabla\mathbf{U}) - \rho\mathbf{g} = -\nabla p - f\_{\sigma} \tag{2}$$

where **U** is the fluid velocity, ρ is the fluid density, *p* is the fluid pressure, μ is the dynamic viscosity, *t* is the time, *g* is the gravity acceleration, and *f*σ is the surface tension. Only the laminar flow is considered in the study.

(**a**)

(**b**)

**Figure 1.** Two different descriptions of frames and the postures of the cuboid payload in the 3D numerical wave tank (NWT). (**a**) The overall description of two frames and a cuboid payload in the 3D NWT; (**b**) description of Euler angles in the top, side, and front view.

#### *2.2. Free Surface Tracking*

The Volume of Fluid (VOF) method is applied for tracking the free surface in OpenFOAM. In the VOF method, a phase function α is defined in each cell, which indicates the quantity of water in the cell. α is 1 if the cell is full of water, and it is 0 in empty cells. On the air-water interface, the value of α is between 0 and 1. The fluid density ρ and the dynamic viscosity μ in each cell are calculated with the equations,

$$\begin{cases} \rho = a\rho\_1 + (1 - \alpha)\rho\_2\\ \mu = a\mu\_1 + (1 - \alpha)\mu\_2 \end{cases} \tag{3}$$

where the subscripts 1 and 2 mean the values of water and air, respectively. The phase function α can be determined by solving an advection equation,

$$\frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha \mathbf{U}) + \nabla \cdot (\alpha (1 - \alpha) \mathbf{U}\_a) = 0 \tag{4}$$

where the last term on the left-hand side is an artificial compression term and **U**α is the relative compression velocity [25].

#### *2.3. Waves2Foam Library and WaveFoam Solver*

The library waves2Foam is used to generate regular waves. The boundary condition and solve procedures are listed below.

### 2.3.1. Waves2Foam Library

The library waves2Foam is a toolbox for generating and absorbing water waves [26]. Waves are generated at the inlet and absorbed at the outlet.

The velocities of regular waves are based on the linear Stokes' wave theory,

$$u(\mathbf{x}, z, t) = \frac{g k A}{\omega} \frac{\cosh k(z + h)}{\cosh kh} \sin \varphi \tag{5}$$

$$w(\mathbf{x}, z, t) = \frac{g k \mathcal{A}}{\omega} \frac{\sinh k(z + h)}{\cos kh} \cos q \tag{6}$$

where *u*(*x*, *z*, *t*) is the horizontal velocity distribution, *A* is the wave amplitude, ω is the wave frequency, ϕ = *kx* − ω*t* and *k* is the wave number, *h* and is the water depth.

The relaxation zone technique is used to absorb waves at the outlet. The relaxation function is

$$a\_R(\chi\_R) = 1 - \frac{\exp(\chi\_R^{3.5}) - 1}{\exp(1) - 1} \text{ for } \chi\_R \in [0:1] \tag{7}$$

It is applied into the relaxation zone as follows,

$$
\lambda = \alpha\_{\mathbb{R}} \lambda\_{\text{computed}} + (1 - \alpha\_{\mathbb{R}}) \lambda\_{\text{target}} \tag{8}
$$

where λ is either **U** or α. The variation of α*<sup>R</sup>* is the same as given in [27], and χ*<sup>R</sup>* represents a certain point in the relaxation zone. The definition of χ*<sup>R</sup>* is such that it is always 1 at the interface between the nonrelaxed part of the computational domain and the relaxation zone.
