Two-Way Coupling Simulation of Fluid-Multibody Dynamics for Estimating Power Generation Performance of Point Absorber Wave Energy Converters

: The objective of the present study is to develop and validate a two-way coupling simulation method between viscous ﬂuid and multibody dynamics to estimate the power generation performance of point absorber wave energy converters. For numerical analysis of ﬂuid dynamics, an enhanced density correction model was proposed to improve the accuracy and stability of the pressure calculation in DualSPHysics, an open-source code based on smoothed particle hydrodynamics (SPH). Through 2D hydrostatic and wave generation simulations, it was seen that the relative error in the average pressure was reduced from 11.81% to 1.64%. In addition, an interaction interface was developed to enable coupling simulation with RecurDyn, a commercial software for the simulation of multibody dynamics. Simulations were performed for a 3D single-body cylinder with a simple shape and a two-body ﬂoating wave energy converter (WEC) in regular waves, varying the linear damping coeﬃcient of the power take - oﬀ (PTO) system to verify the prop osed coupling simulation method for ﬂuid -multibody dynamics. The results were benchmarked against experimental data, revealing a relative error of 1.05% with the experimental results when employing a high damping coeﬃcient for the PTO system. Furthermore, to improve the eﬃciency of the two - body WEC, two design modiﬁcations were suggested; their impact on power generation performance improvement was examined. The developed method is anticipated to contribute to research aiming to enhance the power gene ration eﬃciency of various wave power devices with multiple elements and joints, pending further validation and reﬁnement of the simulation approach.


Introduction
According to the US Environmental Protection Agency (EPA) report [1], environmental problems such as global warming, extreme weather events, and air pollution caused by emissions of greenhouse gases (e.g., carbon dioxide, methane, etc.), primarily generated from the combustion of fossil fuels, have aggravated to a serious level in recent years.As a countermeasure against this problem, environmental regulations (e.g., IMO 2020, Paris Agreement) have been reinforced at a global level.According to the data presented in the Statistical Review of World Energy by British Petroleum [2], the demand for ecofriendly renewable energy (solar, wind, and wave energy) has sharply increased as alternative energy resources.Ocean waves are estimated to generate vast amounts energy, equivalent to 100 to 400% of the world's total electrical energy demand (see Figure 1).Wave energy converters (WEC) are devices that absorb wave energy from the sea and convert it into electrical energy.The oil crisis in the 1970s has resulted in an increase in the preference for using alternative energy sources instead of oil, and research related to alternative energy and WEC has gained momentum [4].Thereafter, the design and research on various types of WEC have shown some progress; however, compared to other types of renewable energy, the number of cases that have reached the commercialization stage is very small [5,6].As outlined in Figure 2, WECs can be largely classified into oscillating water columns, oscillating body systems, and overtopping devices, based on the working principle [6].Among the various types of WECs, oscillating body systems are designed on the principle of converting kinetic energy from bodies that move sensitively in response to sea surface motions caused by waves into electrical energy.Considering the wave energy directly interacts with the float, the energy extraction efficiency is high; thus, these systems account for a significant proportion of the listed types of WECs.Additionally, within the subtypes of oscillating body systems, the point absorber operates through the motion of a floating body moved by sea waves to generate power.In many simulation studies, the focus is solely on the heave motion for simplification purposes.However, some researchers have incorporated multiple degrees of freedom into the modeling of the floating body [7,8].
With the point absorber type of WEC, multiple structures of the same size can be connected for use, or a hybrid system can be developed via combined installation with other renewable energy devices, which increases the energy generation efficiency per unit area and improves the economic value.Considering these advantages, active research on point observers has been conducted [9][10][11].
In addition, various numerical studies [12,13] have been conducted for the efficient design of WEC considering a significantly complex problem requiring multiphysics simulation (hydrodynamics, multibody dynamics, power take-off system analysis, control, etc.) encompassing different fields.Among these various analysis elements, it is essential to develop a technique for multibody interaction analysis that consists of flow and various constraints based on hydrodynamics.For numerical studies on the hydrodynamic performance of WEC, the potential theory or boundary element method (BEM) has been mainly employed from the early stage to the present [14][15][16].The wave energy converter simulator (WEC-Sim), jointly developed through a research project between the US Sandia National Laboratory (SNL) and the National Renewable Energy Laboratory (NREL), is the most representative and comprehensive simulation program for numerical analysis of WEC based on linear analysis.
However, according to Li and Yu [17], although the above-mentioned methods based on linear analysis have advantages in terms of calculation speed in a comparative analysis between hydrodynamic simulation methods for WEC, they are limited when the accurate effect of viscosity or modeling of nonlinear interactions between the waves and WEC are considered in the simulation.Therefore, recent studies have conducted numerical analysis based on the Navier-Stokes (NS) equation, a suitable approach for addressing these problems.Numerical analysis methods based on NS equations can be largely classified into grid-based methods using the Eulerian approach and particle-based methods using the Lagrangian approach.In general, grid-based analysis methods have been widely applied to conventional WEC modeling.Additionally, simulations utilizing OpenFOAM, an opensource grid-based software, have been performed to model the interaction between various structures and waves [18][19][20].
However, these methods involve complexities, such as algorithm construction for grid generation and movement, which are difficult to conduct for problems with motions with large displacement [21].Moreover, additional equations need to be solved, such as the volume-of-fluid (VOF) method [22] or the Level-Set method [23], for implementing extremely nonlinear free surface motions such as wave breaking [24].Therefore, recent studies have reported on the analysis of WEC performance using particle-based CFD simulation to replace grid-based methods [25,26].
Representative particle-based simulation methods include smoothed particle hydrodynamics (SPH) [27] and moving particle semi-implicit (MPS) [28] methods.In this study, DualSPHysics v4.4.12 [29], an open-source code based on the Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH) method, which has the merit of fast computation speed, is used because in the calculation process, the pressure is solved based on the equation of state using explicit methods.Particle-based simulation methods have a problem of pressure fluctuations during flow analysis due to Lagrangian characteristics from the early stage of the research, with various studies conducted to address this problem.In particular, with DualSPHysics v4.4.12, the δ-SPH method (Molteni and Colagrossi [30]), a density correction method that stabilizes the pressure by adding a density diffusion term to the temporal change in density to stabilize pressure, was introduced to improve the problem of pressure fluctuation.Antuono et al. [31] proposed an improved method to calculate a more accurate diffusion term by excluding the first-order normalized density gradient when calculating the term and adding the diffusive effect to density change.Unlike these methods that add a diffusion term to the temporal change in density, various existing studies discuss a different approach using the Shepard density filter, a method of improving density distribution by filtering the spatial distribution of density [32][33][34].In a recent study, Ren et al. [34] proposed a density filtering method, wherein the conventional Shepard density filter was used for fluid particles and wall particles, a sum with an appropriate proportion between the average density considering the difference in hydrostatic pressure between the density information of the particles and the surrounding fluid particles.However, these improvement methods still have minor problems, such as the discontinuity in the pressure of the particles located on the wall of the object and the pressure of the fluid particles near the wall.
In the present study, to analyze the performance of WEC considering the interactions of fluid-multibody dynamics, the final aim of the research, we developed a two-way coupling [35] simulation method between the enhanced DualSPHysics v4.4.12 with improved stability of pressure calculations and RecurDyn [36], a commercial program for multibody dynamics simulation capable of implementing various constraint conditions.Simulations were performed with a 3D simple-shaped one-body cylinder and a two-body floating WEC according to changes in the power take-off (PTO) damping coefficient to validate the developed simulation system.The simulation results were comparatively analyzed with the experimental data presented by Zang et al. [37] and Yu et al. [38].Finally, two ideas on shape change were proposed to improve the efficiency of the two-body WEC: using the developed simulation method, the improvement in the efficiency of power generation performance according to the shape change was calculated, with the effect of shape change comparatively analyzed.In this study, a High-Performance Computer (HPC) equipped with an Intel i5-4690 CPU @ 3.5 GHz (4 cores), 32 GB of RAM, and a GeForce RTX 3080 GPU was utilized for all simulations.

Governing Equations
In this study, the governing equations of fluid dynamics are continuity Equation (1) and momentum Equation (2); they are presented as follows [29]: where , , and  denote density, time, and velocity, respectively;  and  are tensors indicating directions in the space,  represents total stress, and  indicates external force.
where , , and  denote pressure, mass, and gravitational acceleration, respectively;   ��� = 0.5(  +   ),  2 = 0.01ℎ 2 , and  represent the kernel function the values of which are given according to the distance between particles; furthermore, in this study, the quintic kernel [40] is used as the kernel function.Π  denotes the viscous term with artificial viscosity, as proposed by Monaghan [39].In addition, for additional pressure stabilization, the  -SPH method of Molteni and Colagrossi [30] was applied to the continuity equation, as shown in Equation (3), and coefficient  was set as 0.1.In this case, applying the -SPH method, the type of fluid particles can be classified into inner fluid (IF) and near-wall fluid (NWF) particles based on the distance from the wall (  ), as shown in Figure 3. Wall (W) particles were used to represent the boundary conditions, resulting in three types of fluid particles in total, and the density correction model was applied to these particles.In the case of the default setting in conventional DualSPHysics v4.4.12, the -SPH method is only applied to IF particles.Finally, the pressure was calculated using the  obtained from Equation (3) through the equation of state as in Equation ( 5) below: where  = 7,  =  0 2  0 / ,  0 = 1000  −3 indicates the reference density, and  0 = ( 0 ) = �(/)�  0 denotes the speed of sound in the reference density.

Enhancement of DualSPHysics
In DualSPHysics v4.4.12, the density correction model [30] for stabilizing the density field is applied by categorizing it into three parts: inner fluid (IF), near-wall fluid (NWF), and wall (N) particles, as shown in Figure 3; that is, the model is applied to IF particles but not for NWF and W particles.Additional details on the default settings of DualSPHysics can be found in Crespo et al. [29].
However, instability is still present in the pressure calculation of DualSPHysics v4.4.12, along with the problem of non-physical formation of a relative gap between the fluid and the wall particles close to the free surface.Therefore, to address this problem, we propose an improved density correction model, i.e., the enhanced -SPH method [31], which is still employed for IF particles.However, for NWF and  particles, considering the buffered density filter that was only applied to the W particles by Ren et al. [34], a hydrostatic Shepard filter that considers the difference in hydrostatic pressure is newly proposed and applied, as given in Equation ( 6): The newly proposed hydrostatic Shepard filter enhances filtering accuracy by accounting for differences in hydrostatic pressure compared to the traditional Shepard filter.It offers the benefit of eliminating the need for the artificial coefficient χ, which is required in the buffered density filter in Equation ( 7) [28]: The implementation of the newly proposed density correction model enables more accurate calculation of wall pressure compared to DualSPHysics v4.4.12.It is explained in Section 3.1.1.However, in cases involving high-speed flow, it is necessary to employ a very small computational timestep.This approach prevents excessive particle movement, ensuring that the pressure is sufficient to push fluid particles out without allowing them pass through the wall particles.As this is another major problem among the existing particle methods, the dynamic stabilization (DS) model [41] was employed in this study with a slight modification to prevent such problems.In Tsuruta et al. [41], to maintain a uniform particle distribution in the flow field, if two particles overlapped by a certain amount, the DS model was applied to shift the particles from each other and to all particles.However, in this study, the DS model was applied only to the fluid particles near the wall to prevent escape.As shown in Figure 4, the repulsive force,    , calculated using Equation ( 8) in Tsuruta et al. [37] is applied to the fluid particle along the direction of the relative position vector.For the coefficient (α  ) which determines the permissible rate of penetration, α  = 0.9 was used in Tsuruta et al. [41], but α  = 0.8 was used to minimize the effect on the internal fluid flow in this study.In this study, the ratio of the timestep to the Courant number was increased from 0.1 to 0.2.Despite this, the application of Equation ( 6) enhanced the accuracy and stability of density calculations, which in turn affect pressure, thereby offsetting potential issues.While the enhanced methods in this study might lead to a marginally higher computational load compared to traditional DualSPHysics approaches, these challenges can be efficiently overcome by leveraging high-performance computing hardware.

𝐷𝐷𝜈𝜈 ⃗ 𝐷𝐷𝐷𝐷
where   denotes the volume of particle i,  ⃗  is the relative position vector between two particles,   = α    and   is the particle diameter.

Multibody Dynamics (RecurDyn)
In multibody dynamics (MBD), it is assumed that the object or system to be analyzed is a rigid body without deformation, and the motion of the rigid bodies is tracked over time.Considering the constraints, the equation according to the position, velocity, and acceleration can be expressed as Equation (9) as follows [42]: where  denotes the generalized coordinates of the system; , , Φ indicate the independent coordinate, the dependent coordinate, and the kinematic constraint; Φ  represents the Jacobian matrix obtained by partial differentiation of the constraint in the number of coordinates.
By introducing the Lagrange multiplier theory to the constraint equation composed of algebraic equations, the differential-algebraic equation (DAE) that combines the equations of motion composed of differential equations can be expressed as Equation (10) below, serving as the governing equation for rigid body motion with constraints included [43].
where ,  and  denote the Lagrange multiplier vector, mass matrix, and a force vector that includes the external force. is the term obtained by second-order differentiation of the constraint equation, Φ, and can be represented as the third term on the right side of Equation (9).In this study, a multibody dynamics simulation was performed using RecurDyn V9R3 [36], a commercial program that can set and implement various forces and joint constraint conditions.The details of the numerical method can be found at http://www.functionbay.com(accessed on 15 Feb 2024), the website of FunctionBay.

Two-Way Coupling Simulation between Fluid and Multibody Dynamics
First, to define the interaction forces between the fluid and multibody for 2-way coupling simulation, it is necessary to consider the contact force between the body elements and fluid.In the present study, the dynamic boundary condition (DBC) suggested by Crespo et al. [43] was employed.The DBC calculates W particles in the same governing equation as IF particles and applies a repulsive force by increasing the pressure as the NWF particles approach the W particles within an effective radius.However, unlike fluid particles, each W particle does not move individually, whereas all W particles belonging to the corresponding body element are stationary or moved by a determined motion function.
As depicted in Figure 5, to calculate the force and torque (  ) acting on the body by applying the DBC, the repulsive force acting on each W particle is integrated with respect to the center of gravity; the MBD analysis is performed by including the external force on the right-hand side of Equation (10).Moreover, the coupling simulation is conducted by adding the repulsive force (  ) acting on each fluid particle to momentum Equation ( 4) of the fluid.In this manner, the fluid and multibody are time-integrated by exchanging the interaction forces at every time step after being independently performed through each analysis technique.For additional details, please refer to Yun et al. [35], wherein the MPS method was used instead of DualSPHysics to solve fluid motion.
In the 2-way coupling simulation system, a dynamic link library (DLL) file for the enhanced DualSPHysics solver was created, as mentioned above; the interaction force for coupling analysis was exchanged and time-integrated until the termination time under the control of RecurDyn.

D hydrostatic Problem
To verify the enhanced hybrid density correction model proposed in this study, simulations were performed on the 2D hydrostatic problem, which has been mainly used by a number of researchers to verify the improvement of pressure calculation in the existing research based on the particle method, and on the problem related to 2D wave generation, which is a simulation of hydrodynamic problems necessary for WEC analysis.
For the 2D hydrostatic problem, as shown in Figure 6, a rectangular tank was filled with fresh water at a depth and width of 0.4 m and left as is for 20 s to perform the simulation.The total calculation time was 20 s, and the pressure time series was monitored at the tank bottom point P1.
First, a convergence test of the particle size was performed according to the conditions outlined in Table 1.The analytic solution and time series according to particle size change are shown in Figure 7a, and the average dimensionless pressure at the time of pressure stabilization (15-20 s) after the initial pressure fluctuation is shown with the number of particles used in Figure 7b.The trend of convergence can be seen with a decrease in the particle size, and the grid convergence index (GCI) [44] based on three cases of small particle sizes (0.01, 0.005, and 0.0025 m) was calculated as 0.23%.Consequently, considering the calculation efficiency based on GCI, 0.0025 m was confirmed as the most optimal size for this hydrostatic problem.Next, for the case with converged particle size, the result of comparing the pressure distribution according to the depth of the fluid particles at 20 s with the default Du-alSPHysics is shown in Figure 8.In particular, it can be qualitatively confirmed from the figure that the pressure stability was improved near the free surface and wall boundary conditions by the enhanced density correction model proposed in the present study, and a result closer to the analytic solution was obtained.In addition, the pressure time-series at the P1 point of this study were compared with the results of the conventional Du-alSPHysics with the same particle size as shown in Figure 9. Compared to analytical solution for the average pressure at the time of pressure stabilization (15-20 s) after the initial pressure fluctuation, the relative error was 1.64% when the method from the present study was used, whereas with the conventional DualSPHysics, the relative error was 11.81%, thus quantitatively verifying that the accuracy of the pressure calculation was significantly improved.Therefore, through the 2D hydrostatic pressure test, it can be observed that the stability and accuracy of the pressure calculation can be improved when the enhanced method proposed in this study is applied.

D Wave Generation Problem
According to the manual of DualSPHysics [29] (https://github.com/DualSPHysics/DualSPHysics/wiki/3.-SPH-formulation#38-boundary-conditions,accessed on 15 Feb 2024), the application of dynamic boundary conditions (DBC) [42] introduces a limitation, which includes formation of unphysically large boundary layers between the fluid and the wall particles near the wall boundary.In particular, the layers have a more significant impact in hydrodynamic problems.Although these unphysical boundary layers can be reduced by decreasing the size of the particles, in this case, the computation efficiency may also be decreased.Therefore, as shown in Figure 10, under conditions similar to those of the numerical wave tank of Ren et al. [34], a 2D wave generation test, a hydrodynamic problem, was performed to investigate whether the problem of unphysically large boundary layers can be improved using the enhanced density correction model.However, the reflected wave at the right end of the numerical wave tank was suppressed using the sponge layer function provided by DualSPHysics.Based on the results of the 2D particle size convergence test described in the previous section, a converged particle size of 0.0025 m was used for the simulation.In particular, from comparing the pressure contours near the wavemaker boundary and the bottom surface at t = 7.00 s, as shown in Figure 11a, unphysically large boundary layers were observed along with the occurrence of discontinuity and instability of the pressure field between W and NWF particles in the conventional DualSPHysics.However, after applying the proposed method with enhanced density correction, as shown in Figure 11b, it can be qualitatively confirmed that the pressure field was significantly stabilized, and the unphysical boundary layers almost disappeared.Accordingly, through the 2D hydrostatic and wave generation tests, as described above, the results confirmed that with the application of the method developed in the present study, the stability and accuracy of the pressure calculation can be improved, and the unphysically large boundary layers near the wall boundary, which is a limitation of the conventional DualSPHysics, may also be improved.

Definition of Problem and Grid Convergence Tests
As a preliminary step prior to the application of the proposed method to the WEC of a complex shape, the simulation of the heave motion of a simple 3D one-body cylinder in regular waves was performed with and without the PTO system and according to the changes in the damping coefficient to verify the proposed method of coupling simulation of fluid-multibody dynamics, as shown in Figure 12.In conditions similar to those of the experiments performed in Zang et al. [37], 1-DOF simulation was performed under regular waves (wave height: 0.16 m, period: 1.5 s, wavelength: 3.395 m).Before the main simulation, convergence tests of the particle size for the 3D wave generation simulation were also performed, as shown in Table 2.The trend of convergence can be seen with a decrease in particle size, and the grid convergence index (GCI) [44] based on three cases of small particle sizes (0.08, 0.04, and 0.02 m) was calculated as 1.15%.Finally, considering the calculation efficiency based on GCI, 0.02 m (1/8 size of wave height) was confirmed as the optimal size for the 3D wave generation problem.A numerical wave tank simulation was performed to generate the target waves without a cylinder.From the pressure contour snapshot shown in Figure 13, the pressure calculation can be observed to perform stably without instability in the pressure field, even with the application of the coupling simulation.The wave elevation time series was measured at the position where the cylinder was to be placed, and the results were compared with the analytical solution, as shown in Figure 14.Comparing the average values of wave elevation, it can be seen that relatively accurate wave generation is possible through the coupling simulation of the fluid-multibody dynamics proposed in the present study with a relative error of 0.204%.However, the average water level is approximately 1% higher than the theoretical value.This indicates that the difference is associated with the result in which the water level in the hydrostatic problem was approximately 1% higher than that in the analytical solution after the initial flow field was stabilized, as described in the hydrostatic test in Section 3.1.Next, under the verified wave generation conditions, the cylinder was placed without a PTO system, and modeling was performed with the constraint condition that only allows the heave motion, as shown in Figure 15. Figure 16 shows the results of the time series of the heave velocity of the cylinder in comparison to the experimental results.The experimental results show that the peak values of heave velocity are relatively non-uniform; however, the coupling simulation results in this study showed a large relative error of 8.71% compared to the experimental result for the average peaks over three cycles.One of the major causes of overestimation in the coupling simulation is that the actual experiment consisted of a device for fixing the cylinder so as to only have the heave motion.Furthermore, this simulation did not consider the consequent mechanical damping (MD), but the resistance caused by this fixing experimental device was modeled as velocity damping.Therefore, to estimate mechanical damping, several velocity damping coefficient tests based on the heave motion of the cylinder were performed, and linear interpolation was performed on the results.The results showed that, in the experiment [37], the velocity damping coefficient from mechanical damping was close to approximately 50 Ns/m, and when this value was introduced in the simulation, the relative error with the experimental result was reduced to 1.59%.Next, as shown in Figure 17, a simulation was performed using the PTO system.Velocity damping by the PTO system that was measured in the same experiment [37] was added to mechanical damping, as described above, and two PTO damping coefficients (240 and 1100 Ns/m) were applied.It can be seen that both the results shown in Figure 18a with a small damping coefficient and Figure 18b with a large damping coefficient were smaller than the experimental results, with relative errors of 16.97% and 8.77%, respectively.
As described above, through a 3D cylinder test with a PTO system, the coupling simulation of fluid-multibody dynamics was performed with stable calculation results, and the trend of decreasing heave motion with increasing PTO damping coefficient was properly reproduced in the simulation.Moreover, the relative error was significantly influenced by the Power Take-Off (PTO) damping coefficient.Among the various factors contributing to this issue, the inconsistency of heave velocity in the experimental data, as depicted in Figure 16, stands out.Future simulation studies should focus on enhancing the accuracy of the coupling simulation relative to the PTO damping coefficient and conduct further verifications.

Definition of Problem
To verify the proposed method of coupling simulation of fluid-multibody dynamics for a WEC with a complex shape, as shown in Figure 19, a coupling simulation was performed for 3-DOF motions (surge, heave, and pitch) under regular waves for a 1/33-scale reference Model 3 (RM3), an experimental model of a two-body point absorber WEC of Yu et al. [38] released by the US National Renewable Energy Laboratory (NREL) for WEC research considering a PTO system.As shown in Figure 20, a numerical wave-making tank was constructed considering the computational efficiency under conditions similar to the modeling conditions for RM3.In this case, as shown in Figure 21, RM3 was moored at four points at 90° intervals by the spring force on both sidewalls on the horizontal plane near the stationary water surface.First, a preliminary simulation was performed to generate accurate target regular waves (wave height: 0.091 m, period: 1.393 s, wavelength: 3.018 m) at the point where RM3 was placed.For particle size, 0.011 m (approximately 155 million particles), which is approximately 1/8 of the target wave height, was used.

Simulation Results
Figure 22 shows a snapshot of the pressure field owing to the wave generation in the numerical wave tank, and it can be confirmed that the simulation can be conducted with stable calculation results.In addition, as shown in Figure 23, measuring the wave elevation time series at the position where RM3 is to be placed and comparing the average wave elevation with that of the analytic solution based on the second-order Stokes theory, it can be seen that the target wave elevation is accurately generated within approximately 1% of the relative error.Next, RM3 was placed based on accurate wave generation, and as shown in Figure 24, multibody modeling for RecurDyn was performed to simulate conditions similar to the experimental conditions (four-point mooring, PTO system, relative heave motion constraint between two floats, and 3-DOF constraint of RM3).In this case, since the thickness of the plate (0.003 m) is thinner than particle size (0.011 m), there is a possibility that the buoyancy force in the simulation may be larger than the actual force.Therefore, to address this problem, the simulation was performed by adding the mass of the increased volume (approximately 26.95%) to the total mass (spar/plate + 12.3 kg) due to the use of a larger particle size compared to the actual plate.Table 3 lists the properties of RM3 used in the simulation.
The velocity-based damping force, which is a constraint condition for modeling the PTO system, can be expressed by Equation ( 11): Based on this equation, the instantaneous absorbed power with time can be calculated, as shown in Equation ( 12), and the final averaged absorbed power is calculated using Equation (13).
() =   () , () = c   , () 2 (12) where   () is the constraint force between the two floats by the PTO system according to time,   is the damping coefficient of the PTO system,  , () is the relative heave velocity between the two floats with time,   () is the instantaneous absorbed power with time,   is the average absorbed power,  0 is the start of a new period, and  is the period.Based on the multibody dynamics modeling described above, simulations were performed to verify the proposed method for the three different PTO damping coefficients measured in the experiment, as listed in Table 4. From the snapshot of the pressure field in the central section of the numerical wave tank, a stable calculation of the pressure distribution around the WEC can be confirmed, as shown in the simulation (Figure 25). Figure 26a shows the velocity time series for the relative heave motion between the float and spar/plate.Figure 26b shows the time series of the instantaneous absorbed power calculated according to Equation (12).Similar to the results of single-body cylinder simulation in Section 3.2, the relative heave velocity decreased with an increase in the PTO damping coefficient.
Finally, during the period of approximately 16~17 s when the pattern of the motions is considered to converge after terminating the simulation, the average absorbed power was calculated according to Equation (13), and the value was compared with the experimental result shown in Figure 27.The calculated values of the averaged absorbed power were 1.657, 1.673, and 1.394 W, in the order from the lowest PTO damping coefficient; moreover, compared with the experiment, the relative errors were observed to be 30.78,22.48, and 3.42%.The cause of the relative error in the experimental results can be classified in two ways.First, as shown in Figure 28, in the process of representing a plate thinner than the particle size or the slope of a float by a regular arrangement of particles, distortion may have occurred compared to the actual shape, which may have affected the response characteristics of the object motions in waves.In addition, because the thin plate was expressed with relatively large particle sizes and the mass was increased to adjust for the resulting buoyancy force, this had the largest impact on the relative error with the experimental results in the case with the largest relative heave motion because of the low PTO damping coefficient.To solve this problem, it is necessary to further investigate the particle arrangement and related boundary conditions for arbitrary complex shapes in future research.In this regard, several recent studies [45,46] have been reported by researchers employing the particle method.Second, as shown in Figure 29, when the wave elevation time series is compared at a 1 m point in front of the RM3 structure depending on the presence or absence of RM3, a WEC structure, it is observed that the elevation of the incoming waves decreases by approximately 3.3% owing to the waves diffracted by the WEC structure.Compared with the experimental result, because the energy of the waves also varies with the change in the elevation of the incoming waves, there may be relative errors in the prediction of the average absorbed power of the WEC for the target wave elevation.Therefore, to resolve this problem, it is thought that an active wave absorption technique [47], which can reduce the effect of diffracted or reflected waves by actively controlling the wavemaker motions according to the measured wave elevation data in real time, needs to be introduced in the 3D numerical wave tank method used in this study.

Application on Shape Change for Improving Power Generation Efficiency of Two-Body WEC
To further enhance power generation performance, the shape of the model was partially improved based on the RM3 model described in Section 3.3, as shown in Figure 30.In this case, to keep the draft the same as before, the mass was increased based on the buoyancy force by the volume below the water plane of the added shape (Model A type: float + 3.6 kg, Model B type: spar/plate + 3.5 kg).In Model A, a sloped structure was added to the vertical side of the float shape to increase efficiency.That is, in the case of the vertical side before, when the water level rose, the water plane did not change and a constant buoyancy was applied; however, in the case of the sloped structure, the water plane gradually increased; thus, an additional buoyancy force was applied to increase the relative motion in the vertical upward direction.In the B-type model, a thin sloped plate was attached to the lower part of the spar/plate, such that when the pitch motion of the RM3 spar/plate was decreased, the interference was reduced during the relative heave motion so that the efficiency could be increased.
Nevertheless, the two designs proposed in this study entail increased material costs.Furthermore, the ramped edge profile in Figure 30b induces nonlinear motion in the WEC, complicating its control.Future studies should aim for a detailed comparison of both designs with equivalent material costs [48].
Simulations were performed by applying these two ideas of shape improvement.For the PTO damping coefficient, a value of 1306.475Ns/m was applied, which had the smallest relative error with the experimental result.
As shown in Figure 31a, comparing the time series of relative heave velocity, it can be seen that both types of shape change improved motion performance of the model compared to the existing shape of the model.For quantitative verification, the average absorbed power was calculated from the time series of the instantaneous absorbed power, as shown in Figure 31b.Compared to the average absorbed power at 1.394 W of the existing model, the value of the Model A type was 2.478 W (+77.75%) and that of the Model B type was 1.585 W (+13.74%), showing an increase in both cases.Therefore, the efficiency for both shape changes was confirmed to be improved.In addition, further improvement in efficiency is expected when parameter optimization for the two new shapes proposed in the present study is performed in follow-up research.

Results
In the present study, a coupling simulation method was developed between Du-alSPHysics, an open-source-based particle method code for CFD simulation, and Re-curDyn, a commercial program for elastic multibody dynamics, to analyze the power generation performance of the point absorber WEC in waves.The conclusions obtained from this study are summarized as follows: For flow analysis, an enhanced density correction model was applied to address the limitations of conventional DualSPHysics in relation to the discontinuity and instability of the pressure field near the wall or object.Through 2D hydrostatic and wave generation simulations, a more continuous and stable pressure distribution was obtained near the free surface and wall boundary.In addition, the average value of the pressure measured at the bottom was reduced to 1.64% in the developed model compared to the relative error of 11.81% for conventional DualSPHysics, thus confirming that the accuracy of the pressure calculation was also improved.
A coupling simulation system between the enhanced flow analysis solver and a commercial program for multibody dynamics, RecurDyn, was developed.Simulations were performed on the 1-DOF heave motion of a 3D one-body cylinder and the 3-DOF motions of the commercial program under regular waves to verify the developed coupling simulation.In the simulation of the heave motion of a 3D one-body cylinder without a PTO system, as a result of estimating the mechanical damping coefficient (50 Ns/m) of the device for limiting the motion to heave motion in the actual experiment and applying the estimated coefficient, the relative error with the experimental result was 1.05%.In addition, in the simulation of the 3D two-body point absorber WEC (RM3), which considers 3-DOF motion under regular waves for three different PTO damping coefficients, the average absorbed power was estimated and compared with that of the experiment.The larger the coefficient, the smaller the relative errors compared with the experiment, with the errors in the range of 30.78%~3.42%.
Finally, based on the RM3 model, a point absorber WEC, two ideas of shape change for efficiency improvement were applied, and a simulation was performed accordingly.As a result of calculating the average absorbed power for Model A (shape change in the float) and Model B (shape change in the plate) and comparing the result with that of the existing shape, the average absorbed power increased by 77.75% and 13.74%, respectively, confirming that the efficiency was improved for both ideas of the shape change.

Discussion
The coupling simulation technique developed in this study holds promise for application in the power generation performance analysis of various floating WECs.It proves to be a valuable tool in addressing engineering challenges that involve the complicated interaction of viscous flow with the free-surface and the dynamic coupling of multibody systems featuring diverse elements and joints.
Successful validation will enable virtual performance comparisons and optimization of WEC designs by considering viscosity effects under realistic sea conditions.It is expected that this approach has the potential to significantly reduce costs associated with experimental testing during the design phase of WEC development.
Future research should focus on validating simulations of irregular wave generation, building on the work of Liu and Frigaard [49], to faithfully reproduce irregular waves observed in real sea.This should include particle size convergence tests over a broad spectrum of regular wave scenarios.

Figure 2 .
Figure 2. Classification of WECs based on a working principle modified from [6].

Figure 8 .
Figure 8. Pressure of fluid particles according to the depth at 20 s.

Figure 9 .
Figure 9. Time-history of pressure at the P1 point.

Figure 13 .
Figure 13.Snapshot of 3D wave generation for the cylinder problem (pressure contour).

Figure 14 .
Figure 14.Comparison of wave elevation time-series with the analytic solution.

Figure 15 .
Figure 15.Snapshot of 3D wave generation for the cylinder problem without PTO (pressure contour).

Figure 17 .
Figure 17.Snapshot of 3D wave generation for the cylinder problem with PTO (pressure contour).

Figure 20 .
Figure 20.Schematic of a 3D numerical wave tank with an RM3 model.

Figure 21 .
Figure 21.Constraint for modeling mooring lines by linear spring force.

Figure 22 .
Figure 22.Snapshot of 3D wave generation for the RM3 problem (pressure contour).

Figure 24 .
Figure 24.RecurDyn modeling constraint of multibody dynamics of the RM3 problem with PTO.

Figure 28 .
Figure 28.Difference between CAD file (stl) and particle distribution.

Figure 29 .
Figure 29.Comparison of wave elevation between W and W/O in front of RM3.

Figure 30 .
Figure 30.Modification of RM3 for improving power generation efficiency: (a) Modified float (Model A type); (b) Modified spar/plate (Model B type).

Figure 31 .
Figure 31.Comparison of heave velocity and instantaneous absorbed power with Model A and B types: (a) Relative heave velocity; (b) Instantaneous absorbed power of RM3.

Table 2 .
Convergence tests regarding the particle size for the 3D cylinder problem.Viscosity model Artificial viscosity ( = 0.01) Fluid density (kg/m 3 ) 1000 3.2.2.Heave Motion of Single-Body Cylinder with and without PTO in Regular Waves

Table 3 .
Properties of RM3 for co-simulation.

Table 4 .
Simulation condition for the RM3 problem.