Numerical Study on Wave Radiation by a Barge with Large Amplitudes and Frequencies

: A two-dimensional boundary element method is used to study the hydrodynamics of a single barge with prescribed motions of large amplitudes and high frequencies. The wave radiation problem is solved in the time domain based on the fully nonlinear potential ﬂow theory. For numerical simulations, special treatments like plunging wave cutting and remeshing approaches are presented in detail. The numerical schemes are veriﬁed through comparing with analytical results. Both the generated outgoing wave amplitudes and hydrodynamic coe ﬃ cients can be calculated with su ﬃ cient accuracy. Then, we focus on large heave, sway and roll motions to investigate the nonlinear e ﬀ ects on hydrodynamic forces, respectively. In particular, the heave motion with two frequencies is also simulated to study the interactions between results at di ﬀ erent frequencies. It is interesting to see the sum and di ﬀ erence frequency components and the envelopes in time histories as a result. For forces caused by forced sway or roll motions, there are only even-order harmonics for vertical forces and only odd-order harmonics for horizontal forces. Finally, a single body with combined sway, heave and roll motion is studied to examine the interactions between motion modes.


Introduction
The wave-body interaction problem has long been one of the focuses of hydrodynamics. The study starts from wave interactions with one single body. The understanding of the physical mechanism of hydrodynamic interaction of one body can shed light on the wave interactions with multiple bodies. The study on radiation problem can provide features of hydrodynamic properties of loads acting on the ocean structures and deformations of free surface profile.
The wave radiation problem is related to a rigid body in forced motions in otherwise calm water, which can generate outgoing waves. The oscillating fluid pressure resulting from the fluid motions and radiated waves is the source of hydrodynamic forces. The study of wave radiation problem dates at least back to Ursell [1]. He studied the harmonic heave motion of a horizontal circular cylinder in infinite water depth. The very long cylinder assumption made the problem to be reduced to two-dimensional (2D). Later, Kim [2] proposed a solution of the potential problem associated with the harmonic oscillation of an ellipse (2D) or ellipsoid (3D). Green's function was introduced to represent the potential as the solution of an integral equation, which was obtained numerically. The hydrodynamic coefficients of all six motion modes were provided. Black et al. [3] calculated the radiation problem including heave, sway and roll of a rectangular section (2D) and a vertical circular cylinder (3D) by employing Schwinger's variational formulation. An analytical solution to the heave radiation of a rectangular body was presented by Lee [4], who divided the whole fluid region into three sub-regions and solved the nonhomogeneous boundary value problem separately. Li et al. [5] considered the radiation problem of a circular cylinder submerged below an ice sheet with a crack.
On the free surface f S , expressed as η , there exist both kinematic and dynamic conditions. Mathematically, in the Eulerian form, they can be expressed as And in their Lagrangian form, they have where / / D Dt t φ = ∂ ∂ + ∇ ⋅∇ is the material derivative following a given fluid particle.
On the wetted body surface b S , impermeable condition should be satisfied because no fluid particle can penetrate or leave the rigid body surface. It implies that the normal velocity of the fluid particle equals the normal component of the body surface velocity. Therefore, it yields where / n ∂ ∂ represents the partial differentiation along the normal direction For initial conditions, since the outgoing wave is excited by the body's motion, the water is assumed to be undisturbed initially. Thus, the initial conditions on the free surface have the form Since velocity potential φ exists, the mass continuity equation is reduced to the Laplace equation On the free surface S f , expressed as η, there exist both kinematic and dynamic conditions. Mathematically, in the Eulerian form, they can be expressed as And in their Lagrangian form, they have where D/Dt = ∂/∂t + ∇φ · ∇ is the material derivative following a given fluid particle.
On the wetted body surface S b , impermeable condition should be satisfied because no fluid particle can penetrate or leave the rigid body surface. It implies that the normal velocity of the fluid particle equals the normal component of the body surface velocity. Therefore, it yields where ∂/∂n represents the partial differentiation along the normal direction → n = (n x , n z ) to the body surface, pointing out of the fluid domain. → U = (U, V) denotes the translational velocities of the body motion in the x and z direction, respectively. Ω is the rotational velocity about the body's rotational center, which is positive in the anticlockwise direction. → X = (X, Z) is the position vector in the body-fixed system relative to the rotational center. The impermeable condition also applies to the sea bottom S 0 and the two truncated boundaries S c , which gives ∂φ ∂n = 0. on S 0 , S c For initial conditions, since the outgoing wave is excited by the body's motion, the water is assumed to be undisturbed initially. Thus, the initial conditions on the free surface have the form The above equations completely define the initial-boundary value problem (IBVP) for velocity potential φ. The whole flow field information including velocity and pressure field can be obtained after solving the IBVP.
To satisfy the radiation condition, an artificial damping zone is added to minimize the reflection of the truncated boundary S c . This is achieved by adding a damping term in Equation (3) artificially. So, they become The two parameters α and β are to control the damping strength and length, respectively. λ is the wavelength obtained from the dispersion relation, and ω is the oscillation frequency of the body motion. The performance of the damping zone depends on the combination of α and β. The criterion for appropriate α and β is to make sure that the radiated outgoing waves damp out gradually inside the damping zone and completely vanish at the truncated boundary. A detailed example of the selection of α and β for a single floating body is given by Tanizawa's [20] study.

Hydrodynamic Forces and Moments
The hydrodynamic forces on the body can be obtained by direct integration of the pressure over the instantaneous wetted body surface S b . The pressure in the fluid can be determined through Bernoulli equation The hydrodynamic forces and moments acting on the structures are then obtained through In Equation (9), the time derivative of velocity potential, φ t , is not explicitly given even when φ has been found. The auxiliary function approach [21,22] is adopted to circumvent the need for calculating φ t so that the hydrodynamic forces can be obtained directly.

Numerical Procedures
The initial boundary value problem of velocity potential φ is solved by the boundary element method (BEM), which is derived based on Green's second identity where A p is the solid angle at point P, ∂Ω stands for the closed boundary of the fluid domain (consisting of free surface, wetted body surface, truncated boundaries and sea bottom). Moreover, r pq denotes the distance between two points p and q. Equation (11) is usually referred to as the boundary integral equation, which establishes a relationship between unknown φ at any point and the values of φ and ∂φ/∂n on the boundary. For a 2D problem, the fluid boundary ∂Ω can be discretized into a number of straight line segments. Within each segment, formed by two neighboring nodes, introduce a natural coordinate ξ along it, and postulate φ and ∂φ/∂n as the linear interpolation of the values of the two nodes. Then, rewrite Equation (11) according to the newly discretized boundary, we could obtain N denotes the total number of elements, and (1 + ξ)/2 and (1 − ξ)/2 are shape functions. Finally, applying all the boundary conditions and moving all the unknown terms to the left-hand side and known terms to the right, we can get the matrix equation The subscript S N denotes that the Neumann condition is applied on the corresponding boundaries S b , S c and S 0 . Solving the equation above, the velocity potentials on S N and normal derivatives on the free surface S f at a given time can be obtained. For numerical simulation in the time domain, the BEM should be used together with time stepping method. The BEM deals with the IBVP for velocity potential at a given time and the time stepping method updates the free surface position and the corresponding velocity potential to the next time. A second-order Runge-Kutta time integration is adopted to perform the updating. Let = (x, z, φ) be expressed by a function D /Dt = f (x, z, φ) (the actual form of the derivatives can be found in Equation (7)), then we have where i = (x i , z i , φ i ) are known terms corresponding to time t at the ith step and i+1 are to be determined at time t + ∆t. ∆t is the time step. The fluid velocity on the free surface should be calculated as well after φ and φ n are obtained in order to do the time integration. The calculation of fluid velocity is the same as in Li and Zhang [23]. Particular treatment is needed for intersection points, which belong to both free surface and wetted body surface. What is special here is that there are two different normal directions associated with the same point. Therefore, there are two normal derivatives at each intersection point. A common treatment is to split the normal derivative into two parts in accordance with the element which it belongs to [22]. When dealing with free surface nodes, only the (∂φ/∂n) S f part should be used. The velocity potential is continuous at any point. No special treatment is needed. The velocities at these intersection points are essential because they decide the accuracy of the wave run-up on the body and affect the evaluation of the hydrodynamic loads. Thus, the four-point Lagrangian interpolation, which gives higher accuracy, is used to calculate the velocities at intersection points instead of the linear interpolation. The velocities at the intersection points should be further refined according to the impermeable conditions since these nodes belong to the wetted body surface as well. This treatment not just ensures the higher order accuracy of the velocities, but also maintains smooth deformation of the free surface because the velocities at the intersection points are interpolated through that of the free surface nodes only.

Special Treatments
A thin jet can be formed along the body surface when the body is in large sway or roll motion. When it occurs, the free surface is very close to the body surface. This is likely to cause large numerical errors or even break down the simulation. A simple and common treatment is to cut the thin jet during simulation [24]. For long-time simulations, numerical error will be accumulated step after step. Saw-tooth instability may appear due to the nature of numerical method. To suppress this, the smoothing technique is used. In this paper, two schemes, the five-point-third-order smoothing algorithm [25] and energy smoothing scheme [26], are adopted for different initial meshes.

Plunging Wave Cutting
The plunging wave evolves with time and overturns until hitting the free surface underneath. It is slightly different form the plunging jet in Sun et al. [18] in terms of where the plunging happens. In their paper, the plunging jet occurs close to the body surface, while in this study the plunging wave appears relatively far from the body, usually beyond one wavelength. The simulation will break down once the plunging wave tip touches the free surface underneath. A simple treatment is to cut the plunging wave before it hits the main free surface. The judging condition should be the distance between the lowest point of the plunging wave (Point A in Figure 2 and the free surface underneath. When the distance D jet ≤ d f ∆d (d f is a factor to control the length of the jet to be kept, and ∆d is the reference element size on the free surface), part of the plunging wave should be cut. First, locate the turning point of the plunging wave (Point B in Figure 2) and then find Point C, which is the first node with x value less than that of the turning Point B. Connect points B and C as the new free surface. Obviously, new nodes need to be inserted between B and C.

Special Treatments
A thin jet can be formed along the body surface when the body is in large sway or roll motion. When it occurs, the free surface is very close to the body surface. This is likely to cause large numerical errors or even break down the simulation. A simple and common treatment is to cut the thin jet during simulation [24]. For long-time simulations, numerical error will be accumulated step after step. Saw-tooth instability may appear due to the nature of numerical method. To suppress this, the smoothing technique is used. In this paper, two schemes, the five-point-third-order smoothing algorithm [25] and energy smoothing scheme [26], are adopted for different initial meshes.

Plunging Wave Cutting
The plunging wave evolves with time and overturns until hitting the free surface underneath. It is slightly different form the plunging jet in Sun et al. [18] in terms of where the plunging happens. In their paper, the plunging jet occurs close to the body surface, while in this study the plunging wave appears relatively far from the body, usually beyond one wavelength. The simulation will break down once the plunging wave tip touches the free surface underneath. A simple treatment is to cut the plunging wave before it hits the main free surface. The judging condition should be the distance between the lowest point of the plunging wave (Point A in Figure 2

Remeshing
After some periods of time, the nodes on the boundary may cluster or over stretch, which is a source of numerical inaccuracy and possible instability issue. To avoid over distorted elements, nodes on the boundary should be rearranged every several steps. Check the distance between two adjacent

Remeshing
After some periods of time, the nodes on the boundary may cluster or over stretch, which is a source of numerical inaccuracy and possible instability issue. To avoid over distorted elements, nodes on the boundary should be rearranged every several steps. Check the distance between two adjacent nodes. If they are too close, then delete one of the nodes. Otherwise, add nodes between the two existing nodes. Normally 1.5 ∆l to 2∆l is regarded as too large, and around 0.5 ∆l is regarded as too small. ∆l is the original length of the element. The original nodes on the free surface follow the rule that fine and equal mesh is adopted in the close body region and the size increases gradually until the far-field. The remeshing aims to maintain this node distribution pattern. Considering an original node set X i (i = 1, . . . , N) on the free surface, the node set after remeshing is denoted as X k (k = 1, . . . , N ) with known element size L k (k = 1, . . . , N − 1) and the number of new nodes N is unknown. The remeshing is performed in the following three steps. The first step is to calculate the whole arc length of the boundary to be re-discretized. Since we have known all the original nodes on the free surface and the length of each element L i , the total length of the first j elements S j is calculated through Accordingly, the whole length of the original node set is S N . The second step is to determine the number of new nodes N because the expected nodes' distribution on the original curve is known. If we introduce another symbol S to record the sum of the first j elements of the remeshed curve, we have And the total length of the remeshed node set is S N . Now, N can be determined by the two inequalities The final step is to distribute the new nodes X k on the free surface. The first and last nodes are unchanged during remeshing. So, we have The following part demonstrates how to place nodes X k in between. There are three possibilities, as illustrated in Figure 3. The first one means that only one new node can be located on the original element j − 1. Thus, we have then The second possibility happens when the adjacent original nodes are too close, so no new node can be placed on the elements S k−1 ≤ S j−1 and S k ≤ S j+Q and S k+1 > S j+Q , then Q is the number of original nodes to be skipped. The third possibility occurs when the adjacent nodes are too far away from each other, meaning that remeshing is not performed at every step. Let us say thath new nodes can be placed on the element (h is an integer). That is then nodes. If they are too close, then delete one of the nodes. Otherwise, add nodes between the two existing nodes. Normally 1.5 l Δ to 2 l Δ is regarded as too large, and around 0.5 l Δ is regarded as too small. l Δ is the original length of the element. The original nodes on the free surface follow the rule that fine and equal mesh is adopted in the close body region and the size increases gradually until the far-field. The remeshing aims to maintain this node distribution pattern. Considering an original node set i X ( ) 1, ..., i N = on the free surface, the node set after remeshing is denoted as k X ′ ( ) unknown. The remeshing is performed in the following three steps. The first step is to calculate the whole arc length of the boundary to be re-discretized. Since we have known all the original nodes on the free surface and the length of each element i L , the total length of the first j elements j S is calculated through Accordingly, the whole length of the original node set is N S .
The second step is to determine the number of new nodes ' N because the expected nodes' distribution on the original curve is known. If we introduce another symbol ' S to record the sum of the first j elements of the remeshed curve, we have And the total length of the remeshed node set is N can be determined by the two inequalities ' 1 ' and .
The final step is to distribute the new nodes k X ′ on the free surface. The first and last nodes are unchanged during remeshing. So, we have , .
The following part demonstrates how to place nodes k X ′ in between. There are three possibilities, as illustrated in Figure 3. The first one means that only one new node can be located on the original element 1 j − . Thus, we have then ( ) Special attention should be paid to the second to last node X N −1 because the very last node X N is not calculated through the remeshing formula but remains the same as X N . Thus, the distance between these two nodes may not be desirable. To correct the possible issue of the distance, node X N −1 is adjusted to the mid-point of X N −2 and X N . That is An advantage of this scheme is that it guarantees the quality of the mesh and is quite simple to implement in the programme since no additional equation needs to be solved. This, however, can only be used in 2D problems.

Numerical Results and Discussion
This section firstly provides the convergence study as validation tests. Then, we focus on large respective heave, sway and roll motions of a barge with different frequencies to look at the nonlinear effects due to the increase in amplitude of each degree of freedom. In particular, the heave motion with two frequencies is also simulated to study the interactions between results at different frequencies.
In terms of the hydrodynamic loads, nonlinearity is manifested in the occurrence of the higher order harmonic forces. Finally, a single body with combined sway, heave and roll motion is studied to examine interactions between motions at three degrees of freedom. The hydrodynamic forces and moments of combined motion are compared with those of each single mode.

Convergence Study and Validation
The convergence study is done through the results of heave motion, in which the body is subjected to forced motion with displacement z = A sin(ωt). The initial velocity of the body motion is not zero, which means an abrupt start and probably longer transient period. To allow for a gradual development of the flow field, the velocity is revised by multiplying a modification function F m [27]. In this paper, T m = 2T is adopted to ensure a steady state to be reached soon after the ramp time. T m is the modification time and is related to the period of the body motion T = 2π/ω. The displacement and acceleration of the body motion should be revised accordingly.
The vertical hydrodynamic force F z acting on the body can be obtained via Equation (10). The last term −ρgz in the pressure represents the hydrostatic restoring pressure, which is excluded from the integration of vertical hydrodynamic force. To validate the numerical results, simulations are carried out under the same parameters as in Lee [4] and the results are compared with its linear analytical solutions. The case parameters in Lee [4] are h/D = 3 and B/D = 1. The barge width B = 1 m. Since the analytical solutions are obtained through linear theory, the amplitude A should be small to allow for the linearity assumption. A equals 0.01D in the following simulations. The frequencies of the body motion are set against kh, which is to compare the results in Lee [4] directly. k is the wave number of the generated periodic wave and is linked with the frequency ω through the dispersion relation of finite water depth.
All simulations are carried out in a fluid domain with the body located at the center and with the truncated boundaries placed at x 1 = 5λ + B/2, including the 2λ damping zone. The nodes are initially equally distributed on the wetted surface of the barge with element size ∆d. The free surface part near the body is discretized equally as well with ∆d. Away from the local zone on the free surface, the element size increases gradually until the far end. The biggest element on the free surface should not be larger than λ/20 to ensure that there are at least 40 elements in each side of the damping zone. Larger elements of equal size λ/5 are employed on the truncated boundaries and not changed during mesh convergence study, because they are far away from the barge. The nodes on the seabed right below the structure are placed equally with size 2∆d and nodes on the rest of the seabed with size λ/10 throughout all the simulations. During the simulation, the wetted body surface is remeshed every time step, keeping the segment size as ∆d. The whole free surface is remeshed and smoothed every 5 steps.
The damping zone length is set as 2λ with damping strength α = 0.2 for all the cases. The damping applied is quite effective. It minimizes the reflection and makes the radiated waves steady in time and periodic in time and space.
The convergence tests with reference to element size and time step are conducted based on the case when kh = 2.75, because its frequency lies in the middle of the frequency range. The purpose of this test is to ensure accuracy and to provide reference parameters for the other simulations. The mesh size only refers to the element on the wetted body surface and the local region of the free surface. Three different meshes are tested with ∆t = T/105. The wave runup histories on the right side of the barge for the three meshes, which are given in Figure 4a, show that they are not sensitive to these element sizes. The hydrodynamic force on the barge, shown in Figure 4b, reveals that the coarser mesh ∆d = λ/68 causes a visible difference from the other two meshes. Considering the hydrodynamic forces, the mesh size ∆d = λ/137 is small enough to give convergent results and is chosen to test the temporal convergence. We then rerun the simulations with ∆t = T/52 and ∆t = T/150 and compare the corresponding wave runups and hydrodynamic forces on the barge in Figure 5. We can see that these curves are virtually coincident with each other, which means the larger time step is sufficient to guarantee the temporal convergence.
should not be larger than / 20 λ to ensure that there are at least 40 elements in each side of the damping zone. Larger elements of equal size / 5 λ are employed on the truncated boundaries and not changed during mesh convergence study, because they are far away from the barge. The nodes on the seabed right below the structure are placed equally with size 2 d Δ and nodes on the rest of the seabed with size /10 λ throughout all the simulations. During the simulation, the wetted body surface is remeshed every time step, keeping the segment size as d Δ . The whole free surface is remeshed and smoothed every 5 steps. The damping zone length is set as 2λ with damping strength for all the cases. The damping applied is quite effective. It minimizes the reflection and makes the radiated waves steady in time and periodic in time and space.
The convergence tests with reference to element size and time step are conducted based on the case when 2.75 kh = , because its frequency lies in the middle of the frequency range. The purpose of this test is to ensure accuracy and to provide reference parameters for the other simulations. The mesh size only refers to the element on the wetted body surface and the local region of the free surface. Three different meshes are tested with /105 t T Δ = . The wave runup histories on the right side of the barge for the three meshes, which are given in Figure 4a, show that they are not sensitive to these element sizes. The hydrodynamic force on the barge, shown in Figure 4b, reveals that the coarser mesh / 68 d λ Δ = causes a visible difference from the other two meshes. Considering the hydrodynamic forces, the mesh size is small enough to give convergent results and is chosen to test the temporal convergence. We then rerun the simulations with / 52 t T Δ = and / 150 t T Δ = and compare the corresponding wave runups and hydrodynamic forces on the barge in Figure 5. We can see that these curves are virtually coincident with each other, which means the larger time step is sufficient to guarantee the temporal convergence.   To ensure that the present numerical results are not only convergent but also reliable and accurate, they are compared with the linear analytical solutions presented in Lee [4]. Long-time simulation is conducted to guarantee the accuracy of the added mass and damping coefficients because they are calculated based on the periodic part of the hydrodynamic force history ( ) Given that the hydrodynamic force history becomes periodic at 0 t , the dimensionless added mass coefficient μ and damping coefficient λ are determined through the following equation: where T N is an integer. Free surface elevation is another feature to be considered in radiation To ensure that the present numerical results are not only convergent but also reliable and accurate, they are compared with the linear analytical solutions presented in Lee [4]. Long-time simulation is conducted to guarantee the accuracy of the added mass and damping coefficients because they are calculated based on the periodic part of the hydrodynamic force history F z (t). Given that the hydrodynamic force history becomes periodic at t 0 , the dimensionless added mass coefficient µ and damping coefficient λ are determined through the following equation: where N T is an integer. Free surface elevation is another feature to be considered in radiation problems. It is represented by the amplitude of generated outgoing waves. The generated wave amplitude A w is computed through the surface wave profile in zone λ + B/2 ≤ x ≤ 3λ + B/2, which excludes the local wave close to the body and waves in the damping zone. At each time instant, A w is determined as half of the vertical distance from the crest to trough. The comparison of dimensionless added mass, damping coefficients and generated wave amplitudes is made with analytical solutions and is given in Figure 6. The figure shows an excellent agreement, which indicates that the present numerical model and the code are quite capable of capturing the free surface elevation accurately. Generally speaking, the accuracy of the present numerical results is excellent since the variation tendency and values against kh for all the three physical variables A w , µ, and λ are very well captured compared with the analytical solutions.
accurate, they are compared with the linear analytical solutions presented in Lee [4]. Long-time simulation is conducted to guarantee the accuracy of the added mass and damping coefficients because they are calculated based on the periodic part of the hydrodynamic force history ( ) Given that the hydrodynamic force history becomes periodic at 0 t , the dimensionless added mass coefficient μ and damping coefficient λ are determined through the following equation: where T N is an integer. Free surface elevation is another feature to be considered in radiation problems. It is represented by the amplitude of generated outgoing waves. The generated wave amplitude w A is computed through the surface wave profile in zone excludes the local wave close to the body and waves in the damping zone. At each time instant, w A is determined as half of the vertical distance from the crest to trough. The comparison of dimensionless added mass, damping coefficients and generated wave amplitudes is made with analytical solutions and is given in Figure 6. The figure shows an excellent agreement, which indicates that the present numerical model and the code are quite capable of capturing the free surface elevation accurately. Generally speaking, the accuracy of the present numerical results is excellent since the variation tendency and values against kh for all the three physical variables , , w A μ and λ are very well captured compared with the analytical solutions.

Wave Radiation by a Single Motion Mode
The hydrodynamic properties of a floating barge depend on its body parameter and the motion mode. This section will investigate the nonlinear effects on the hydrodynamics against each motion mode, separately. The mesh used in the following simulations is the same as in the convergence study. The length variables are nondimensionalized by the barge width B. Moreover, the water depth h is assumed to be large compared to the generated wavelength and the body draught. The effect of water depth is ignored and should not be included in the results. The dimensionless parameter kh is therefore replaced by the body motion frequency ω. When analyzing the forces and moments, the contribution of static buoyancy is excluded.

Forced Heave Motion
The results will focus on the nonlinear effects on the hydrodynamic forces, due to the increase in the amplitude of the body motion. The nonlinearity can lead to the existence of higher harmonic force components. The superposition of these higher harmonic components makes the total force history irregular, as illustrated in Figure 7a, which is for D = 0.4 and ω 0 = 1.8 at different oscillation amplitudes. ω 0 = ω 0 B/g denotes the nondimensional frequency for a specific case in the following study. Figure 7b gives their corresponding Fast Fourier Transform (FFT) analysis and shows the amplitudes of the components. It is observed clearly from the figure that the higher harmonic force becomes more pronounced as the oscillation amplitude increases. The first harmonic force, however, still makes up a dominant proportion of the total force. In addition, the phase angle of the total force is not affected by the amplitude A as observed in Figure 7a. force components. The superposition of these higher harmonic components makes the total force history irregular, as illustrated in Figure 7a, which is for denotes the nondimensional frequency for a specific case in the following study. Figure 7b gives their corresponding Fast Fourier Transform (FFT) analysis and shows the amplitudes of the components. It is observed clearly from the figure that the higher harmonic force becomes more pronounced as the oscillation amplitude increases. The first harmonic force, however, still makes up a dominant proportion of the total force. In addition, the phase angle of the total force is not affected by the amplitude A as observed in Figure 7a.  Figure 8a. It can be seen that the larger the amplitude, the bigger the force's peak value. The analysis of their high harmonic components can be done through FFT and is presented in Figure 8b. The second and third harmonic forces grow as the amplitude increases. However, they are still very small compared to first harmonic force even if the amplitude is as large as three quarters of the draught. For heave motion with D = 1, the cases of amplitude A = 0.75 are also simulated. Since the amplitude is very large compared to the draught, the bottom of the body will emerge from water occasionally during simulation if the frequency is also high. This will result in the breakdown of the calculation. Simulations show that the frequency should be lower than ω 0 = 0.8 to keep the body bottom staying in the water. The comparison of vertical forces with different amplitudes when ω 0 = 0.8 is shown in Figure 8a. It can be seen that the larger the amplitude, the bigger the force's peak value. The analysis of their high harmonic components can be done through FFT and is presented in Figure 8b. The second and third harmonic forces grow as the amplitude increases. However, they are still very small compared to first harmonic force even if the amplitude is as large as three quarters of the draught.  The nonlinear features will be more complicated and interesting if the heave motion has multiple frequencies since there will be interactions between components at different frequencies. Let us consider the case with two frequencies. The displacement of the body is expressed as  The nonlinear features will be more complicated and interesting if the heave motion has multiple frequencies since there will be interactions between components at different frequencies. Let us consider the case with two frequencies. The displacement of the body is expressed as z = A 1 sin(ω 1 t) + A 2 sin(ω 2 t). The amplitudes at both frequencies are set as the same with D = 1 and A 1 = A 2 = A = 0.2. The above calculations have shown that the amplitude of body motion cannot be large in order to keep the simulations running. One of the frequencies, denoted as ω 1 , is set as ω 1 = 1.0. The second frequency ω 2 is taken as 0.8, 1.05 and 1.2, respectively. The vertical force histories, given in Figure 9a, show clearly the envelope pattern. The period T in the graph corresponds to the higher one of the two frequencies. The apparent difference lies in the period of the envelope T e . For ω 2 = 0.8 and 1.2, T e is nearly the same because in both cases |ω 2 − ω 1 | = 0.2 and T e depends on their difference. So, when ω 2 = 1.05 is close to ω 1 , the period of envelope becomes very long.
large in order to keep the simulations running. One of the frequencies, denoted as 1 ω′ , is set as  In order to obtain the frequency spectrum of the vertical forces, FFT is performed and given in Figure 9b. The primary frequency components are marked in the figure. As we can see, there exist both sum and difference frequencies in addition to the two excitation frequencies 1 ω′ and 2 ω′ . Due to the existence of difference frequency 2 1 ω ω ′ ′ − , the envelope of vertical force shows a period longer than that of heave motion. These sum and difference frequencies reflect the nonlinear interaction between motions. The spectral analysis of wave runups on the body is also given in Figure 10. The same feature as vertical forces is shown. In order to obtain the frequency spectrum of the vertical forces, FFT is performed and given in Figure 9b. The primary frequency components are marked in the figure. As we can see, there exist both sum and difference frequencies in addition to the two excitation frequencies ω 1 and ω 2 . Due to the existence of difference frequency |ω 2 − ω 1 |, the envelope of vertical force shows a period longer than that of heave motion. These sum and difference frequencies reflect the nonlinear interaction between motions. The spectral analysis of wave runups on the body is also given in Figure 10. The same feature as vertical forces is shown.

Forced Sway Motion
The displacement is prescribed as During simulations, water spray will occur on the free surface when the amplitude of sway motion is large or the frequency of motion is high. Thus, spray cutting is needed, which is done through the method detailed in Section 2. The simulations are conducted with 0.5 A = and 0 0.4 ω′ = , which corresponds to relatively slow motion because of the low frequency. The horizontal and vertical force histories at different draughts are given in Figure  11a,c, respectively. Any differences in these cases are due to the body draughts. It can be seen that the horizontal forces are significantly affected by the draughts. This is because horizontal forces are mainly attributed to the pressure on the side walls, whose magnitudes are directly related to the body draughts. However, the vertical forces are due to the hydrodynamic pressure on the bottom. The FFT analysis of these horizontal and vertical forces is provided in Figure 11b,d, respectively. It is interesting to notice that the leading component of vertical force is the second harmonic force and there are no even order components of horizontal forces. Furthermore, there are large negative mean vertical forces caused by sway motion at this low frequency.
We then consider the cases when the body oscillates very fast horizontally with large frequency Figure 10. Amplitude spectral analysis of wave runups of heave motion with two frequencies.

Forced Sway Motion
The displacement is prescribed as x = A sin(ωt). During simulations, water spray will occur on the free surface when the amplitude of sway motion is large or the frequency of motion is high. Thus, spray cutting is needed, which is done through the method detailed in Section 2. The simulations are conducted with A = 0.5 and ω 0 = 0.4, which corresponds to relatively slow motion because of the low frequency. The horizontal and vertical force histories at different draughts are given in Figure 11a,c, respectively. Any differences in these cases are due to the body draughts. It can be seen that the horizontal forces are significantly affected by the draughts. This is because horizontal forces are mainly attributed to the pressure on the side walls, whose magnitudes are directly related to the body draughts. However, the vertical forces are due to the hydrodynamic pressure on the bottom. The FFT analysis of these horizontal and vertical forces is provided in Figure 11b,d, respectively. It is interesting to notice that the leading component of vertical force is the second harmonic force and there are no even order components of horizontal forces. Furthermore, there are large negative mean vertical forces caused by sway motion at this low frequency.  We then consider the cases when the body oscillates very fast horizontally with large frequency ω 0 = 1.4 and A = 0.5. Their horizontal and vertical force histories are illustrated in Figure 12a,c, respectively. Clear differences can be seen on the peaks and troughs of horizontal forces for different draughts because of the increase in frequency when comparing Figure 11a with Figure 12a. Double peaks appear for deep draught with high frequency. In terms of vertical forces, the high frequency makes them of shallower draught considerably larger than that of deeper draught, while the low frequency allows the draught to have little effect on them. This may be because the hydrodynamic pressure on the bottom is larger closer to the free surface and the high frequency magnifies this difference greatly. The FFT analysis of these forces, which is provided in Figure 12b, gives the frequency components of them. In particular, it is observed in the four spectral graphs, including the above two graphs for ω 0 = 0.4, that there are only even-order harmonic forces for vertical forces and only odd-order harmonic forces for horizontal forces. In fact, this phenomenon is discovered and shown mathematically by Wu [28].
Finally, we will examine the effect of sway motion amplitude with D = 1 and ω 0 = 1.4. The FFT analysis of their horizontal and vertical forces is given in Figure 13a,b, respectively. The amplitude of each order of harmonic force of both horizontal and vertical forces grows with the increase in amplitude as expected. In particular, the fifth harmonic horizontal force and sixth harmonic vertical force become noticeable when the amplitude equals three quarters of the body breadth.  analysis of their horizontal and vertical forces is given in Figure 13a,b, respectively. The amplitude of each order of harmonic force of both horizontal and vertical forces grows with the increase in amplitude as expected. In particular, the fifth harmonic horizontal force and sixth harmonic vertical force become noticeable when the amplitude equals three quarters of the body breadth.

Forced Roll Motion
The roll motion is about the body's longitudinal y -axis and through the gravitational center, which is located at half of the draught below the still water. This means  θ . Therefore, any differences in these cases are due to nonlinear effects, which affect vertical forces the greatest. Additionally, the increase in amplitude will not cause phase shifts of forces. The FFT analysis of forces and moments is provided in Figure 14b,d,f, respectively. According to the conclusions made in Wu [28], there should be only odd-order harmonic components

Forced Roll Motion
The roll motion is about the body's longitudinal y-axis and through the gravitational center, which is located at half of the draught below the still water. This means z g = −D/2. The draught of the body is D = 1. The angular displacement of the body can be expressed as θ = θ 0 sin(ωt). θ 0 is the amplitude of roll motion. Three amplitudes θ 0 = π/15, π/9 and π/6 are considered. For a practical ship, a rolling angle of π/6 is regarded as very large. However, such a scenario still needs to be considered as it poses great danger to ship safety.
During the simulations, a thin jet can be formed along the body surface because of the impact on the free surface. It is cut in this study. The time histories of horizontal and vertical forces and moments are plotted in Figure 14a,c and e for ω 0 = 0.4, respectively. The forces and moments are already normalized by θ 0 . Therefore, any differences in these cases are due to nonlinear effects, which affect vertical forces the greatest. Additionally, the increase in amplitude will not cause phase shifts of forces. The FFT analysis of forces and moments is provided in Figure 14b,d,f, respectively. According to the conclusions made in Wu [28], there should be only odd-order harmonic components existing in horizontal forces and moments, while only even order-harmonic components should be occurring in vertical forces as observed in sway motion cases. The spectral analysis of these forces shows that the present numerical results are consistent with the above statement, especially the horizontal forces. The same observation is made when the frequency is as high as ω 0 = 1.2. Their frequency components are analyzed and shown in Figure 15. One may notice that there are even-order harmonic moments and a small first-order harmonic vertical force when the amplitude of the roll motion is large. This may be because the jet cutting treatment performed damages the mirror images formed about the wetted body surface to some extent.

Wave Radiation by Combined Motion
We have investigated each degree of freedom. The combined motion is now considered. The draught of the body is 1 D = . The displacement of each degree of freedom is specified as before and with the same frequency and no phase difference. The relative phase differences among motion modes are not considered in the present study. The amplitude of both heave and sway motion is The amplitude angle of roll motion is 0 0.2 θ = , or 11.46°. The purpose of this study is to examine the interactions between motion modes. The hydrodynamic forces and moments are compared with that of each single mode.
The frequency of 0 0.8 ω′ = is first simulated. The comparisons of forces and moments are provided in Figure 16. The horizontal force of combined motion is actually a little smaller than that of sway motion only. However, a small second harmonic component appears in the combined motion. As for the vertical force, its mean value and the second harmonic force is excited to a large extent compared to the vertical force in pure heave, while the first harmonic force is unaffected. This can be explained by the conclusions in the previous subsection. Specifically, the sway and roll motion can only trigger the frequency components at  Figure 16e,f. This is because the moment is dependent on both horizontal and vertical forces. The second case is when the frequency of combined motion is set as high as 0 1.4 ω′ = . The motion of the free surface within one period is illustrated in Figure 17. Sprays of the free surface can be seen on both sides of the body. They develop with time and travel in the space. They are cut before touching the main fluid underneath. The hydrodynamic forces and moments are then compared with single mode results in Figure 18. The horizontal force is still smaller in magnitude compared to sway motion force. It, however, consists of more highfrequency components. In comparison, the vertical force and moment are both larger than those of the single mode and contain more high-frequency components, especially the moment. These conclusions are the same as the previous low frequency case.

Wave Radiation by Combined Motion
We have investigated each degree of freedom. The combined motion is now considered. The draught of the body is D = 1. The displacement of each degree of freedom is specified as before and with the same frequency and no phase difference. The relative phase differences among motion modes are not considered in the present study. The amplitude of both heave and sway motion is A = 0.2. The amplitude angle of roll motion is θ 0 = 0.2, or 11.46 • . The purpose of this study is to examine the interactions between motion modes. The hydrodynamic forces and moments are compared with that of each single mode.
The frequency of ω 0 = 0.8 is first simulated. The comparisons of forces and moments are provided in Figure 16. The horizontal force of combined motion is actually a little smaller than that of sway motion only. However, a small second harmonic component appears in the combined motion. As for the vertical force, its mean value and the second harmonic force is excited to a large extent compared to the vertical force in pure heave, while the first harmonic force is unaffected. This can be explained by the conclusions in the previous subsection. Specifically, the sway and roll motion can only trigger the frequency components at 2nω 0 , n = 0, 1, . . . of vertical force. The multiple mode interactions affect the moment the greatest, as shown in Figure 16e,f. This is because the moment is dependent on both horizontal and vertical forces. The second case is when the frequency of combined motion is set as high as ω 0 = 1.4. The motion of the free surface within one period is illustrated in Figure 17. Sprays of the free surface can be seen on both sides of the body. They develop with time and travel in the space. They are cut before touching the main fluid underneath. The hydrodynamic forces and moments are then compared with single mode results in Figure 18. The horizontal force is still smaller in magnitude compared to sway motion force. It, however, consists of more high-frequency components. In comparison, the vertical force and moment are both larger than those of the single mode and contain more high-frequency components, especially the moment. These conclusions are the same as the previous low frequency case.

Concluding Remarks
This paper mainly studies the wave radiation by a single body with large amplitudes and high frequencies through numerical simulations based on the fully nonlinear potential flow theory in the time domain. The associated initial boundary value problem is established in both the Eulerian and Lagrangian frameworks, which is solved by the boundary element method combined with the time stepping scheme. During simulations, some special numerical treatments are also applied to the free surface, such as remeshing, smoothing, thin jet and plunging wave cutting. In particular, new approaches for plunging wave cutting and remeshing are proposed in detail. They have advantages for dealing with meshes in two-dimensional problems. The numerical schemes are verified through comparing with analytical results. Cases of wave radiation by a single barge of different motion amplitudes and frequencies are considered. The findings below are all based on the results obtained from potential flow theory without considering the viscous effect.

Concluding Remarks
This paper mainly studies the wave radiation by a single body with large amplitudes and high frequencies through numerical simulations based on the fully nonlinear potential flow theory in the time domain. The associated initial boundary value problem is established in both the Eulerian and Lagrangian frameworks, which is solved by the boundary element method combined with the time stepping scheme. During simulations, some special numerical treatments are also applied to the free surface, such as remeshing, smoothing, thin jet and plunging wave cutting. In particular, new approaches for plunging wave cutting and remeshing are proposed in detail. They have advantages for dealing with meshes in two-dimensional problems. The numerical schemes are verified through comparing with analytical results. Cases of wave radiation by a single barge of different motion amplitudes and frequencies are considered. The findings below are all based on the results obtained from potential flow theory without considering the viscous effect.
For wave radiation by a single motion mode, we focus on heave, sway and roll motions to look at the higher order harmonics due to the increase in the body oscillation amplitude, respectively. The increase in the body motion amplitude can enhance the nonlinear effects of the radiation problem. For heave motion with two frequencies, it is interesting to see the sum and difference frequency components and the envelopes in time histories as a result. For forces caused by forced sway or roll motions, there are only even-order harmonics for vertical forces and only odd-order harmonics for horizontal forces. Additionally, the moment due to roll motion has only odd-order harmonics. For a single body with combined sway, heave and roll motions, the hydrodynamic forces and moments are compared with those of each single mode. The horizontal force of combined motion is actually a little smaller than that of sway motion only. However, a small second harmonic component appears in the combined motion. As for the vertical force, its mean value and the second harmonic force is excited to a large extent compared to the vertical force in pure heave, while the first harmonic force is unaffected.