Numerical Study on Unsteady Pressure Distribution on Bulk Carrier in Head Waves with Forward Speed

: This study deals with wave-induced unsteady pressure on a ship moving with a constant forward speed in regular head waves. Two different numerical methods are applied to solve wave– ship interaction problems: a Rankine panel method which adopts velocity potential, and a Cartesian-grid method which solves the momentum and mass conservation equations under the assumption of inviscid and incompressible ﬂuids. Before comparing l1ocal pressure distributions, the computational methods are validated for global quantities, such as ship motion responses and added resistance, by comparison with available experimental data. Then, the computational results and experimental data are compared for hydrodynamic pressure, particularly focusing on the magnitude of the ﬁrst-harmonic component in different sections and vertical locations. Furthermore, the Cartesian-grid method is used to simulate the various wave-amplitude conditions, and the characteristics of the zeroth-, ﬁrst-, and second-harmonic components of wave-induced pressure are investigated. The nonlinearity of pressure distribution is observed mostly from the pressure near the still-water-level of the ship bow and the normalized ﬁrst-harmonic component of wave-induced pressure decreases as the wave steepness increases. Lastly, to understand the local characteristics of wave-induced unsteady pressure, the time-averaged added pressure and added local force are analyzed. It is found that the major contribution of the time-averaged added local force that occurs around the ship stem above the design waterline.


Introduction
An estimation of seakeeping performances is one of the important procedures in safe ship design. Thanks to rapid progress in computing power, many numerical methods, such as potential flow solvers and computational fluid dynamics (CFD), have been applied to estimate a ship's behavior in various environmental conditions [1][2][3]. In addition to motion responses, the added resistance of a ship due to waves also becomes an important topic to estimate a ship's performance in waves [4][5][6][7][8]. However, those quantities are global or integrated values, rather than local or primitive physical values. Thus, the validation of local physical variables is essential to establish reliable numerical methods for the prediction of ship performance in waves. Furthermore, based on the observation of local flows, physical understanding of the flow around a ship in waves can be deepened, and numerical methods can be improved.
In an experiment, it is relatively easy to measure hydrodynamic forces acting on a ship by a dynamometer, whereas the measurement of local variables, such as hydrodynamic Cartesian-grid method can be found in Yang and Kim [17], and the numerical procedure is briefly summarized here.
The right-handed coordinate system which is moving with the constant forward speed U of the ship is introduced and the origin is located at the center of gravity of the ship. The x-axis is parallel to the ship's longitudinal direction pointing to the ship stern and the starboard side of the ship is the positive y-axis. The z-axis is normal to the other two axes, that means positive upward. Many wave-ship interaction problems are dominated by the inertia effect. Thus, the fluid viscosity can be ignored, and the following continuity and Euler equations are solved: where → u is the velocity vector, ρ means the fluid density, p is the pressure, and b is the body-force vector. A fractional-step method is applied for velocity and pressure coupling and the finite-volume method with staggered variable allocation is used for the spatial discretization. The pressure Poisson equation is solved using a multi-grid algorithm to update the pressure in the fluid domain: where ∆t means the time step. The hydrostatic pressure was eliminated for the ship part that was initially submerged in the water, when calculating the time history of pressure. Thus, the pressure in the subsequent section includes both linear and nonlinear hydrodynamic pressure, and hydrostatic pressure above the still water level. An immersed boundary method is used and the volume fraction of each phase in each cell is calculated [18]. The tangent of hyperbola for interface capturing (THINC) and weighted line interface calculation (WLIC) methods [19,20] are applied to update the volume-fraction function of the water phase. In this formulation, incompressible fluid is assumed for both air and water phases and the free surface boundary conditions are implicitly satisfied. The volume fraction of a solid body is calculated by the transformation of a signed distance field between the ship surface and a grid point. The slip boundary condition is satisfied by a volume-weighted formula [17,21]: where→ u indicates the updated fluid velocity, → U body is the velocity of ship surface closest to the center of velocity-control volume, α → u 3 is the volume fraction of the ship in the corresponding velocity-control volume, I is the identity matrix, and N is the normal dyad.
To generate the incident waves and eliminate reflecting waves from the ship, forcing terms are added in the momentum conservation equations. The forcing term is calculated using the difference between numerical solutions and Stokes analytic solutions of fluid velocities and wave elevation. Besides forcing terms, the boundary values for the velocity and the volume-fraction function of the liquid phase are fixed as Stokes analytic solutions. To reduce the transient time, the velocity field, that is a sum of uniform velocity and Stokes solution, is assumed as an initial condition.
The verification and validation of the present Cartesian-grid method for wave-body interaction problems have been made in Yang and Kim [17] and Yang et al. [21]. The sensitivity of time steps is less than the grid sensitivity and almost ignorable if the Courant-Friedrichs-Lewy (CFL) condition is satisfied. For the ship motion problems, the initial time step ∆t 0 = T e /400 and the CFL value 0.4 are recommended, where T e means the encounter wave period. Unlike implicit pressure-correction methods, the present Cartesian-grid method does not require so-called 'outer' iteration because the continuity equation is satisfied by solving the pressure Poisson equation at the projection step. In general, the residual of pressure is less than order of 10 −6 after three or five iterations using a multigrid method, while the velocity field is updated using explicit method and thus an iteration scheme is not necessary. A typical calculation using a single processor required about 3.5 h per one encounter period and the simulation was conducted during 20 encounter periods. Grid resolution is the most critical factor to ensure the convergence of the present numerical solutions. Thus, a grid convergence test was conducted and is summarized in the next section based on the 'grid convergence index' concept [22].

Rankine Panel Method
A numerical program based on a time-domain Rankine panel method, called the computer program for linear and nonlinear Wave-Induced loads and SHip motion (WISH), was developed by Seoul National University under the support of several large shipbuilding companies. This program is used for the simulation of nonlinear ship motions in waves, and it has been extended and applied to the seakeeping problems for cruise ships and offshore structures, hydro-elasticity analysis such as springing, ship maneuvering, and so forth [2].
The body-fixed coordinate system is introduced for a ship advancing with a constant forward speed U. Assuming the ideal fluid and irrotational flow, a velocity potential φ can be introduced. To linearize the boundary value problem, the velocity potential and the wave elevation are decomposed as follows: where Φ is the double-body basis potential, subscript I is the component associated with the incident wave, and subscript d is the component associated with the disturbed wave. Then, by assuming that Φ is the order of O(1) and the others: φ I , φ d , ζ I , and ζ d , are the order of O(ε), the linearized boundary value problem can be summarized as follows: ∂φ d ∂n Here, (m 1 , m 2 , , ∂n means the spatial derivative in the normal direction to the ship surface, S B is the mean body surface, and g is the gravitational acceleration. Green's second identity is applied to solve the above linearized boundary value problem. In the Rankine panel method, the boundary surfaces are discretized into a grid of quadrilateral panels, and Rankine sources are distributed on each panel. Furthermore, the variables, such as the velocity potentials on the boundary surfaces, are represented by the bi-quadratic B-spline basis function, using the values in the neighboring panels.
The first-order pressure about the mean body position can be obtained by Bernoulli's equation and the perturbations of variables with respect to the mean body position, as follows: where, the linear displacement → δ is defined as x , and → ξ T and → ξ R indicate linear translational and rotational displacement vectors. An artificial wave absorbing area is set near the truncated far-field boundary to satisfy the radiation condition. In the absorbing area, the damping term is added in the kinematic free-surface boundary condition and the damping strength gradually increases toward the far-field.
The first-order pressure acting on the ship hull is integrated to calculate the resultant forces and moments. The added resistance in waves is calculated by integrating the second-order pressure, which is obtained by perturbations of physical variables, such as the displacement, pressure, wave elevation, and normal vector with respect to the mean body position. The detailed derivation of the added resistance formula in the Rankine panel method can be found in Joncquez [23] and Kim and Kim [6].

Ship Model and Test Case
The RIOS bulk carrier, which was the target ship in an innovative experiment [14] to measure the unsteady pressure distribution on the whole ship surface, is considered in this study, and Table 1 gives its main particulars. The Froude number, Fn = U/(gL) 1/2 , is equal to 0.18. In the experiment [14], a total of 234 pressure sensors were attached on the RIOS bulk carrier. The Fiber Bragg Gratings (FBG) sensor was adapted, which is based on optical-fiber sensing technology. In this study, the wave-induced pressure for some representative positions of the ship surface is compared with the experimental data. The distance between aft perpendicular (AP) and forward perpendicular (FP) is divided into evenly spaced ordinates, which are equal to 0.0 and 10.0 for AP and FP, respectively. In the present comparison, four sections-ordinate (Ord.) = 9.5, 9.0 (near bow), 5.0 (mid-ship), and 0.5 (near stern)-are selected, and the position of pressure measurement in each section is represented as the angle θ from the keel line. That means 0 • indicates the keel-line, and 90 • indicates the still-water-level, as shown in Figure 1.   wave cases are considered in this study, half domain is used in both numerical calculations. It should be noted that the present Rankine panel method is based on the linear assumption, and thus the ship below the still-water-level is considered, while the full ship up to ship freeboard is modeled in the Cartesian-grid method.

Motion Responses and Added Resistance
To check the grid convergence in the Cartesian-grid method, the calculating procedure of grid convergence index (GCI) suggested by Celik et al. [22] was applied for three different grids. The subscripts 1, 2, and 3 indicate the finest grid, medium grid, and coarsest grid, respectively. The average grid spacing h for each grid is defined as follows: Ord.=10

Motion Responses and Added Resistance
To check the grid convergence in the Cartesian-grid method, the calculating procedure of grid convergence index (GCI) suggested by Celik et al. [22] was applied for three different grids. The subscripts 1, 2, and 3 indicate the finest grid, medium grid, and coarsest grid, respectively. The average grid spacing h for each grid is defined as follows: where N indicates the number of elements in the index set I c . (i, j, k) is the index pair of the i-th, j-th, and k-th grid point in the x-, y-, and z-direction, respectively. L is the ship length, B is the ship beam, T is the ship draught, A is the amplitude of incident wave, xc i , yc j , and zc k are the center positions, and ∆x i , ∆y j , and ∆z k are the size of the i-th, j-th, and k-th cell, respectively. The grid spacing outside of the range in Equation (13) remains the same to generate the incident wave with the same grid resolution, while the inside grid is systematically refined, which has the ratio of the average grid spacing larger than the square root of two. The total number of grids varied from about two million to five million and the GCI can be obtained using the following procedure: Here, sgn means sign function, Table 2 summarizes the ratio of average grid spacing and the resultant GCI values for the magnitudes of heave and pitch motions and added resistance in waves. The heave and pitch motion amplitudes are almost identical, even for different grid spacing, whereas the added resistance is quite sensitive to the grid. The uncertainty of the grid in the added resistance using the Cartesian-grid method is about 10-20% for the resonance case, and about 26% and 38% for the short-wave case. The uncertainty of the grid in the smaller wave amplitude case is larger than that in the large-wave amplitude case. To generate an accurate incident wave in the Cartesian-grid method, the aspect ratio of the grid is important, as well as the number of grids in wave height and wavelength. Thus, it becomes difficult to generate an accurate incident wave for small-wave steepness and short wavelength cases, and the inaccuracy of the incident wave is reflected in the grid uncertainty. However, the maximum difference between the computed added resistance and experimental result is 18% for the resonance case (λ/L = 1.25), with H/λ = 1/125, and the corresponding grid uncertainty is 17%, which is an acceptable range of uncertainty.  Figure 3 shows the magnitude of the wave-induced heave and pitch motions of the RIOS bulk carrier. The normalized motion magnitudes are defined as follows: Here, k means the wave number of the incident wave and the horizontal axis in the figure represents the normalized encounter frequency ω e (L/g) 1/2 as well as the normalized wavelength λ/L. The numerical computations using the Cartesian-grid method were conducted with H/λ = 1/80 and 1/40, where H and λ are the wave height and wavelength of the incident wave. In addition, the smaller wave steepness condition (H/λ = 1/125) was simulated for the cases where the wavelength is 1.25 times longer than the ship length. In References [13,14], the wave steepness was varied from H/λ = (1/200 to 1/56). It is difficult using the Cartesian-grid method to generate the incident wave accurately with small-wave steepness in short-wave cases for both experiment and numerical simulation, while in the Rankine panel method, the linearized free-surface and body boundary conditions were applied. The overall magnitude of heave and pitch motions are similar to each other, whereas for the Cartesian-grid method, smaller heave response can be found near the resonance region (λ/L = 1.25). However, if the amplitude of the incident wave decreases, the magnitude of heave motion becomes close to the others. This nonlinearity is known to be responsible for the cross-coupling damping of heave and pitch if the ship has a large flare angle.  To investigate the dependency of wave steepness more clearly, Figure 5 plots the magnitude of heave and pitch motions and the added resistance in waves using the Cartesian-grid method for different wave steepness. Figure 5 also represents the corresponding uncertainty of grid spacing as an error bar. As discussed before, the vertical motions and the added resistance decrease as the wave steepness increases, especially for the resonance case (λ/L = 1.25). If the steepness of incident wave increases, the bow wave becomes easily broken, and thus the variation of wetted area is not proportional to the wave amplitude. Moreover, the nonlinearity of bow geometry is one of the reasons that the added resistance in waves is not proportional to the square of the wave amplitude.  Figure 4 shows the added resistance of the RIOS bulk carrier. The normalized added resistance is defined as follows: where R AW means the dimensional added resistance. Both numerical codes provide similar results to the experiment in the overall wavelength range. Likewise, in the motion responses, the normalized added resistance decreases as the wave steepness increases in the computational results of the Cartesian-grid method. It should be noted that the present Cartesian-grid method solves Euler equations, which means the fluid viscosity is ignored. The ship motion and added resistance in waves are known as inertia-dominant problems and the viscous effects are secondary. But, the viscosity cannot be ignorable depending on the local Keulegan-Carpenter number. This is beyond the scope of the present study and further studies are required using a flow solver of Navier-Stokes equations with a suitable turbulence model.  To investigate the dependency of wave steepness more clearly, Figure 5 plots the magnitude of heave and pitch motions and the added resistance in waves using the Cartesian-grid method for different wave steepness. Figure 5 also represents the corresponding uncertainty of grid spacing as an error bar. As discussed before, the vertical motions  The difference between experimental data (D) and numerical results (S) is defined as follows: Thus, if the value of numerical results is larger than that of experimental data, the difference becomes negative. The difference values for the amplitude of heave and pitch motions and added resistance are summarized in Table 3. The difference of motion amplitude is smaller than that of added resistance in waves. In general, added resistance in waves is a small quantity to measure from the model tests and thus the uncertainty is also higher than that of motion amplitude. Table 3. Difference between experiment and numerical results for amplitude of heave and pitch motion and added resistance. To investigate the dependency of wave steepness more clearly, Figure 5 plots the magnitude of heave and pitch motions and the added resistance in waves using the Cartesian-grid method for different wave steepness. Figure 5 also represents the corresponding uncertainty of grid spacing as an error bar. As discussed before, the vertical motions and the added resistance decrease as the wave steepness increases, especially for the resonance case (λ/L = 1.25). If the steepness of incident wave increases, the bow wave becomes easily broken, and thus the variation of wetted area is not proportional to the wave amplitude. Moreover, the nonlinearity of bow geometry is one of the reasons that the added resistance in waves is not proportional to the square of the wave amplitude.

Distribution of Unsteady Pressure
The first harmonic amplitudes of the wave-induced pressure distribution over t whole ship surface for three different wavelengths are compared in Figure 6. The fi harmonic amplitude of unsteady pressure is normalized with the maximum value of l ear dynamic pressure of the incident wave ρgA. In each figure set, the left column is t side view and the right column is the front view. The upper figure in the left column the result of the Rankine panel method (WISH), while the lower figure is the result of t Cartesian-grid method (SNU-MHL-CFD). For the short-wave case (λ/L = 0.5), the ship m tions can be ignored, and the incident wave is passing through the mid-ship section. Th the magnitude of first harmonic pressure is similar to that of the incident wave. On t other hand, the first harmonic amplitude of unsteady pressure is almost zero at the m ship for the long-wave condition (λ/L = 2.0) because the ship is surf-riding on the incide wave. The amplitude of the first harmonic from the Rankine panel method near the st water-level is larger than that from the Cartesian-grid method, especially in the ship bo region, because of the linear assumption. Near the bulbous bow and the ship stern, a slig difference can be observed, whereas the overall patterns below the still-water-level a similar to each other. In the bulbous bow and the ship tern, the aspect ratio of the pane higher than that in the other parts, and this gives an inaccurate solution.

Distribution of Unsteady Pressure
The first harmonic amplitudes of the wave-induced pressure distribution over the whole ship surface for three different wavelengths are compared in Figure 6. The first harmonic amplitude of unsteady pressure is normalized with the maximum value of linear dynamic pressure of the incident wave ρgA. In each figure set, the left column is the side view and the right column is the front view. The upper figure in the left column is the result of the Rankine panel method (WISH), while the lower figure is the result of the Cartesian-grid method (SNU-MHL-CFD). For the short-wave case (λ/L = 0.5), the ship motions can be ignored, and the incident wave is passing through the mid-ship section. Thus, the magnitude of first harmonic pressure is similar to that of the incident wave. On the other hand, the first harmonic amplitude of unsteady pressure is almost zero at the mid-ship for the long-wave condition (λ/L = 2.0) because the ship is surf-riding on the incident wave. The amplitude of the first harmonic from the Rankine panel method near the still-water-level is larger than that from the Cartesian-grid method, especially in the ship bow region, because of the linear assumption. Near the bulbous bow and the ship stern, a slight difference can be observed, whereas the overall patterns below the still-water-level are similar to each other. In the bulbous bow and the ship tern, the aspect ratio of the panel is higher than that in the other parts, and this gives an inaccurate solution.   Figure 7 shows the time history of unsteady pressure for Ord. = 9.5. The left column is the resonance case (λ/L = 1.25), while the right column is the short-wave condition (λ/L = 0.5). The pressure is again normalized with the maximum value of linear dynamic pressure of the incident wave ρgA. The magnitude of normalized pressure is more than three for the resonance case, while it is approximately 1.5 for the short-wave condition. The Rankine panel method provides harmonic time series, while the results of the Cartesiangrid method show the half-sine shape for θ = 90, 82.5, and 75° in the resonance case. Those positions are near the still-water-level, and thus the measuring positions are regularly exposed to the air, because of large ship motion. On the other hand, the exposed time to the air becomes shorter in the short-wave condition, because in this condition, the ship motion  Figure 7 shows the time history of unsteady pressure for Ord. = 9.5. The left column is the resonance case (λ/L = 1.25), while the right column is the short-wave condition (λ/L = 0.5). The pressure is again normalized with the maximum value of linear dynamic pressure of the incident wave ρgA. The magnitude of normalized pressure is more than three for the resonance case, while it is approximately 1.5 for the short-wave condition. The Rankine panel method provides harmonic time series, while the results of the Cartesiangrid method show the half-sine shape for θ = 90, 82.5, and 75 • in the resonance case. Those positions are near the still-water-level, and thus the measuring positions are regularly exposed to the air, because of large ship motion. On the other hand, the exposed time to the air becomes shorter in the short-wave condition, because in this condition, the ship motion can be ignored. The pressure time histories at the other positions show similar magnitude and oscillation period between the two numerical computations.
In the case of stern flow (Ord. = 0.5), the two numerical computations show different magnitudes of wave-induced pressure for the resonance case, even for deeply submerged positions, as shown in Figure 8. In particular, the pressure magnitudes around the crest are different, whereas those around the trough region show similar values to each other. This implies that there are higher harmonic components in the results of the Cartesian-grid method. The magnitude of the unsteady pressure is very small in the stern region for the short-wave condition because of shadowing effects. It should be noted that the time histories of pressure in the Cartesian-grid method show slightly oscillatory behavior. It is known that an immersed boundary method suffers spurious oscillation in pressure fields, especially when a solid body moves, and a grid point is abruptly changed from one phase to another [24]. In the case of stern flow (Ord. = 0.5), the two numerical computations show different magnitudes of wave-induced pressure for the resonance case, even for deeply submerged positions, as shown in Figure 8. In particular, the pressure magnitudes around the crest are different, whereas those around the trough region show similar values to each other. This implies that there are higher harmonic components in the results of the Cartesiangrid method. The magnitude of the unsteady pressure is very small in the stern region for the short-wave condition because of shadowing effects. It should be noted that the time histories of pressure in the Cartesian-grid method show slightly oscillatory behavior. It is known that an immersed boundary method suffers spurious oscillation in pressure fields,  especially when a solid body moves, and a grid point is abruptly changed from one phase to another [24]. The wave-induced pressure of the Cartesian-grid method and experiment can be decomposed into the Fourier components. Figures 9-11 show the magnitude of the firstharmonic component p (1) for different wavelengths. The incident wave steepness is different in the experiment, and it varies from H/λ = (1/200 to 1/56) (λ/L = 2.0 to 0.5), while H/λ = (1/125 and 1/40) cases are simulated using the Cartesian-grid method. The uncertainty of the grid in the Cartesian-grid method is represented as an error bar. The maximum grid uncertainty is 15% near the still-water-level of the bow region in the short-wave condition. Other positions and conditions show less than 5% grid uncertainty. The wave-induced pressure of the Cartesian-grid method and experiment can be decomposed into the Fourier components. Figures 9-11 show the magnitude of the firstharmonic component p (1) for different wavelengths. The incident wave steepness is different in the experiment, and it varies from H/λ = (1/200 to 1/56) (λ/L = 2.0 to 0.5), while H/λ = (1/125 and 1/40) cases are simulated using the Cartesian-grid method. The uncertainty of the grid in the Cartesian-grid method is represented as an error bar. The maximum grid uncertainty is 15% near the still-water-level of the bow region in the short-wave condition. Other positions and conditions show less than 5% grid uncertainty. theory that the linear dynamic pressure is exponentially decaying in the vertical direction with the factor of kz, where k is the wave number. The discrepancy implies that in these cases, the scattered waves-radiation and/or diffraction waves-are equally important to the incident wave. The first-harmonic component of wave-induced pressure is distributed almost uniformly along the vertical direction in the stern region, and its magnitude is less than or equal to half of the linear dynamic pressure.    The magnitude of the first-harmonic component of wave-induced pressure is higher near the ship bow than for the other parts of the ship, and the maximum value is about 1.5 times the maximum value of linear dynamic pressure of the incident wave ρgA for the short-wave (λ/L = 0.5) and long-wave (λ/L = 2.0) cases. In the resonance condition (λ/L = 1.25), the maximum of the first-harmonic pressure is about four times larger than the linear dynamic pressure, because of bow submergence. Near the still-water-level, pressure sensors are regularly exposed to the air, because of the relative motion between ship and water surface, and at those positions, the pressure becomes zero. Thus, the first-harmonic component of pressure near the still-water-level from the Cartesian-grid method and experiment is smaller than that of the Rankine panel method, in which the pressure was calculated for the mean body position. The magnitude of the first-harmonic components of the wave-induced pressure gradually decreases to zero above the still-water-level in the results of the Cartesian-grid method.
A similar tendency can be found in the results at the mid-ship section for the shortwave case (λ/L = 0.5), whereas the value is very small for the resonance condition (λ/L = 1.25) and the long-wave case (λ/L = 2.0). All those results are different from the linear wave theory that the linear dynamic pressure is exponentially decaying in the vertical direction with the factor of kz, where k is the wave number. The discrepancy implies that in these cases, the scattered waves-radiation and/or diffraction waves-are equally important to the incident wave. The first-harmonic component of wave-induced pressure is distributed almost uniformly along the vertical direction in the stern region, and its magnitude is less than or equal to half of the linear dynamic pressure. Figures 12-14 compare the magnitude of the zeroth-, first-, and second-harmonic components of wave-induced pressure for λ/L = (0.5, 1.25, and 2.0), based on the results of the Cartesian-grid method. The steady pressure p s in calm water condition is subtracted from the zeroth-harmonic component. This quantity is closely related to the added resistance in waves and has been defined as the time-averaged added pressure. In the case of ship bow (Ord. = 9.5), it is a positive value near the still-water-level because of wave inertia and this effect is exponentially decaying in the vertical direction. There exists negative added pressure due to suction pressure, which is corresponding to the velocity square term in Bernoulli's equation. However, if the relative ship motion becomes negligible (λ/L = 0.5 and 2.0), the difference of inertia and suction effects is small, which means the added resistance becomes small. In contrast, the added pressure shows relatively large magnitude in the resonance condition, as shown in Figure 13a, and this implies that the added resistance is larger than that in the short-wave or long-wave conditions. It is also clearly observed that the positive added pressure occurs near the still-water-level, while the negative value exists below the still-water-level. Moreover, the vertical distribution range of positive added pressure increases with increasing wave amplitude.
range of positive added pressure increases with increasing wave amplitude.
The normalized values of first-harmonic component for the short-wave condition are less sensitive to wave steepness than that for the resonance condition. In other words, the first-harmonic component is proportional to the wave amplitude for the short-wave condition. On the contrary, if the wave amplitude increases, the first-harmonic component decreases for the resonance condition, as shown in Figure 13b. It implies that the firstharmonic component is proportional to the lower order, rather than to the first order of the wave amplitude for this condition. The nonlinear effect of the incoming wave amplitude can be clearly seen in the bow region, whereas this effect is diminished near the stern section. In the stern section, the nonlinear effects of the hull geometry and viscosity are more important than the incoming wave itself. The magnitudes of the second-harmonic components in both the bow and stern regions are very small, compared with that of the first-harmonic component in the bow region. In the long-wave condition, the zeroth-and second-harmonic components are close to zero, and the first-harmonic component is slightly reduced as the incident wave steepness increases. However, the degree of change according to the amplitudes of incident wave is smaller than that in the resonance condition because the relative wave height is different.   To investigate the characteristic of added resistance in head waves more clearly, Figure 15 plots the time-averaged added pressure and the added local force in the longitudinal direction obtained using the Cartesian-grid method. In the previous study [17], the added pressure Δpk was defined as follows: The normalized values of first-harmonic component for the short-wave condition are less sensitive to wave steepness than that for the resonance condition. In other words, the first-harmonic component is proportional to the wave amplitude for the short-wave condition. On the contrary, if the wave amplitude increases, the first-harmonic component decreases for the resonance condition, as shown in Figure 13b. It implies that the firstharmonic component is proportional to the lower order, rather than to the first order of the wave amplitude for this condition. The nonlinear effect of the incoming wave amplitude can be clearly seen in the bow region, whereas this effect is diminished near the stern section. In the stern section, the nonlinear effects of the hull geometry and viscosity are more important than the incoming wave itself. The magnitudes of the second-harmonic components in both the bow and stern regions are very small, compared with that of the first-harmonic component in the bow region. In the long-wave condition, the zerothand second-harmonic components are close to zero, and the first-harmonic component is slightly reduced as the incident wave steepness increases. However, the degree of change according to the amplitudes of incident wave is smaller than that in the resonance condition because the relative wave height is different.
To investigate the characteristic of added resistance in head waves more clearly, Figure 15 plots the time-averaged added pressure and the added local force in the longitudinal direction obtained using the Cartesian-grid method. In the previous study [17], the added pressure ∆p k was defined as follows: where p k (t) is the unsteady pressure at the k-th triangular surface mesh and superscript '0 indicates the steady pressure for the calm-water condition, and N face is the total number of triangular surface mesh. where pk(t) is the unsteady pressure at the k-th triangular surface mesh and superscript indicates the steady pressure for the calm-water condition, and Nface is the total number triangular surface mesh.  The added pressure is a meaningful quantity only for a fixed ship or a ship in short waves, where the ship motion can be neglected. In the present study, the added local force ∆ → F k is defined by considering the variation of surface normal vector in time, as follows: where → n k (t) is the surface normal vector of the k-th triangular surface mesh and superscript '0 indicates the mean value of surface normal vector for the calm-water condition. It should be noted that the surface area is assumed to be unity, because in the present Cartesian-grid method, the change of wetted area is automatically reflected in the pressure value. In this study, only the x-component of added local force is considered, because it is directly related to the added resistance in waves.
As the wavelength increases, the ship motion becomes large, and thus the timeaveraged added pressure increases and the considerable added pressure is occurred in the wide range of the bow. If the wavelength is longer than that for the resonance condition (λ/L = 1.25), the ship tends to move up and down following the incoming waves. Thus, the relative wave elevation would decrease and consequently, the added pressure becomes smaller. The added local force shows much more clearly the contribution of the ship bow to the added resistance in waves. Because the variation of the surface normal vector in time is also reflected in the added local force, the contribution from the ship shoulder becomes less significant, whereas the 'hot spot' is observed around the ship stem above the design waterline. This result also indicates why the sharp bow, such as the Ax-bow type, reduces the added resistance even in the resonance condition, compared with the conventional bow shape.

Conclusions
This study compared the computational and experimental results of heave motion, pitch motion, and added resistance of the RIOS bulk carrier in head waves. Moreover, it investigated the distribution of unsteady, wave-induced pressure on the ship surface. The results show that the overall accuracy of both numerical methods is acceptable, and the following conclusions can be drawn:

•
The nonlinearity of pressure distribution was observed mostly from the pressure near the still-water-level of the ship bow, where the measurement position was regularly exposed to the air. Consequently, the pressure time history near the still-water-level showed half-sine shape, and the first-harmonic component was smaller than that from the Rankine panel method, in which the linearized boundary value problem is solved.

•
The time-averaged added pressure, which is defined as the difference between the zeroth-harmonic component of the wave-induced pressure and the steady pressure in calm-water condition, tended to show positive value near the still-water-level of the bow, whereas negative value was observed below the still-water-level. This tendency can be clearly seen for the resonance wave condition in this study, and the vertical distribution range of positive added pressure became wider as the amplitude of incident wave increased.

•
The normalized first-harmonic component of wave-induced pressure decreased as the wave steepness increased, especially in the ship bow section for the resonance wave condition. This nonlinearity indicates that the variation of wetted surface is not proportional to the amplitude of the incident waves, and it also affects the added resistance due to waves for different wave steepness.

•
The added local force, which includes the variation of pressure and surface normal vector in time due to the ship motion, was introduced. The major contribution of the time-averaged added local force that occurs around the ship stem above the design waterline and the distribution of the added local force in the longitudinal direction will provide useful information to understand the added resistance in waves in more detail.