A meta-model for added resistance in waves

In this paper, we perform a comprehensive study of various semi-empirical methods using publicly accessible experimental data on added resistance in waves with different ship types and conditions. Based on the analysis results, a new method (So-called ‘‘Combined method’’) is proposed, combining two existing methods, which are available in arbitrary wave headings. The results from the two methods are combined smoothly using a tangent hyperbolic function according to wavelengths and wave headings. The coefficients constituting the function are tuned to minimize mean squared error between predictions and model experiments. Finally, the new Combined method is verified by full-scale measurements of a general cargo ship and a container ship, and it seems to give good agreement with measurements in all analysis areas, compared to existing semi-empirical methods. Especially, it showed better performance in estimating added wave resistance at high waves, resonance frequencies, arbitrary waves, and low speeds.


Introduction
An operating ship experiences additional resistance due to the surrounding weather conditions, resulting in speed reduction and increased fuel consumption, which can be directly related to greenhouse gas emissions.Traditionally, this fact has been of great interest to ship designers and operators from the perspective of speed/power performance.Moreover, with the increasing interest in atmospheric environmental issues recently, IMO has set the EEDI to limit greenhouse gas emissions.In this regard, it is even more necessary to estimate the added resistance of a ship in an accurate and efficient way for the initial design and the management of operations of a ship.
Many theoretical methods have been developed for calculating the added resistance of ships.Havelock (1942) proposed a method of calculating added resistance by integrating longitudinal pressure on the wetted surface of a ship, and Boese (1970) developed a near-field direct pressure integration method using strip theory.Maruo (1957) first introduced the far-field method based on momentum conservation, and it was expanded in later studies (Joosen, 1966;Maruo, 1960Maruo, , 1963)).Radiated energy approach based on Maruo's far-field method was introduced by Gerritsma and Beukelman (1972), and Salvesen (1978) achieved satisfactory results by applying it to the motion of the ship obtained from the strip theory.Faltinsen (1980) presented an asymptotic formula, assuming the added resistance of wall-sided hull forms in short waves.
Panel methods based on potential theory for computing added wave resistance has been extensively studied by many authors (Joncquez, 2009;Kim and Kim, 2011;Seo et al., 2013;Söding et al., 2014;Lee et al., 2021).However, since most approaches were limited to linear theory, it was generally difficult to accurately calculate the nonlinear effect.There were also non-linear panel methods, but they had problems with stability and robustness, and long computational time.Meanwhile, along with the improvement of computational power, Computational Fluid Dynamics (CFD) method based on Reynolds-Averaged Navier-Stokes (RANS) has been widely applied (Orihara and Miyata, 2003;Guo et al., 2012;Sadat-Hosseini et al., 2013;Simonsen et al., 2014;Sigmund and El Moctar, 2018;Lee et al., 2021;T. Kim et al., 2021).It had the advantage of being able to consider nonlinear effects and showed good results overall.However, the output results from the 3D panel method and the RANS equations solver vary depending on the calculation grid and large computational time is required, which leads to a struggle in terms of practicability (Shigunov et al., 2018).Another problem is that they require detailed hull shapes to predict added resistance, which in some cases could serve as an important constraint.
Alternatively, simplified methods based on theory and experimental results have been developed, which could easily estimate added resistance with only a few ship parameters compared to the methods covered earlier.The semi-empirical formula for the added resistance due to wave reflection was first proposed by Fujii (1975), and later further tuned based on more experimental data by Takahashi (1988) and Tsujimoto et al. (2008).In parallel with these studies, for the ship motion-induced added resistance, Jinkine and Ferdinande (1974) https://doi.org/10.1016/j.oceaneng.2022 Ship's design speed waves for sea trial conditions (Boom et al., 2013).The STAWAVE-1 method assumes that the wave reflection contribution dominates the added resistance.From this approach, a practical equation that simplifies the reflection-induced added resistance in irregular waves by approximating the waterline geometry on the bow section and the beam of the ship was presented.Contrary to this, STAWAVE-2 method considers both reflection and radiation contribution in estimating the transfer function for the added wave resistance.Liu and Papanikolaou (2016) originally proposed a statistical method of combining Faltinsen (1980) and Jinkine and Ferdinande (1974).In subsequent studies (Liu andPapanikolaou, 2019, 2020), they introduced wave heading-based trigonometric functions to their previous equation and expanded it to enable calculation for small draft, ballast conditions, and arbitrary waves by regression analysis based on extensive experimental data.Lang and Mao (2020) proposed an added wave resistance model for head seas based on Tsujimoto et al. (2008) and Jinkine and Ferdinande (1974).It was influenced by the formulas presented in Liu and Papanikolaou's paper (Liu and Papanikolaou, 2016;IMO, 2016), and some of its calculations were modified using their experimental datasets.The proposed method was further updated to allow the calculation of the peak position in arbitrary waves by introducing an encountered frequency correction factor (Lang and Mao, 2021).
There are also some simple equations that can directly calculate the added resistance in irregular waves such as Kreitner's method (Kreitner, 1939;ITTC, 2005) and Shopera (Papanikolaou et al., 2015).Most simplified methods have been developed to estimate the resistance of a ship operating in head seas.Although studies for added resistance in arbitrary wave headings have been continuously conducted in recent years, comparative analysis and insight into these methods are still insufficient.
To this regard, we perform a comprehensive study of various semi-empirical methods using publicly accessible experimental data on added resistance in waves with different ship types and conditions, which is presented in Section 2. Thereafter, in Section 3, a new method is proposed combining two methods available in arbitrary wave headings from the studies of Lang and Mao (2021) and Liu and Papanikolaou (2020), which has been elaborately verified in previous research.The results from these two methods are smoothly connected to wave conditions and wave heading through a combining function.This method is developed for the calculation of added resistance for large fleets of ships, so that robustness, computational efficiency, and applicability to a range of different ship types are priorities.The coefficients of this function are tuned using extensive model test data, and their values are presented according to the ship type.Section 4 verifies the performance of the corresponding method using full-scale measurements of two different ships.Conclusions of the study are addressed in Section 5.

Description of the experimental data
The results of various publicly accessible experimental data were obtained to consider the general applicability of added resistance in waves to the fleet level, and the performances of the existing semi-empirical models were compared and analyzed.Fig. 1 shows the distribution of main dimensionless parameters of the ships and experimental conditions used in the study.The whisker in the figure indicates the maximum and the minimum range.The box plot shows 25% and 75% quantiles and the circular marker in it represents the median.The data set consists of a total of 2559 samples of approximately 49 ships and 255 different experimental cases.Most of the experiments were conducted at design loading conditions, without trim.More detailed information on the data set is shown in Tables A.1-A.6 in Appendix A. Liu and Papanikolaou (2016) introduced the length of entrance (  ) parameter in Faltinsen's asymptotic approach to reflect the hull form influence on the component of added resistance due to diffraction effect.Thereafter, corresponding parameters were used in Liu andPapanikolaou (2019, 2020), and Lang andMao (2020, 2021).  is defined as the horizontal distance from the point where the length of the waterline surface reaches 99% of the breadth to Forepeak (Conversely, Length of run (  ) is the horizontal distance from the point where the length of the waterline surface reaches 99% of the breadth to the endpoint of the waterline), as shown in Fig. 2.These are necessary factors for calculating the entrance angle used in the wave reflection contribution to the added resistance.However,   and   cannot be accurately estimated without detailed hull shape information of the ship.Some authors (Liu et al., 2016;Lang andMao, 2020, 2021) also listed such values of several ships in their papers.

Parameter estimation for 𝐿 𝐸 and 𝐿 𝑅
Fig. 3 shows   and   according to the length (  ) of the ships under the design loading conditions we have secured.Overall,   and   increase in proportion to   .  decreases as block coefficient (  ) increases, but in the case of   , the trend according to the   is not clear.In Fig. 4, dimensionless   (  ∕  ) and dimensionless   (  ∕  ) are plotted as a functions of Block coefficient (  ), and a linear regression line for each ship type is also plotted.As   increases, the dimensionless   decreases, and there was a slight difference in the slope and intercept values of the regression line depending on the ship type.On the other hand, the dimensionless   according to   shows a clear difference in trend according to ship type.The dimensionless   values of the tanker, liquefied gas carrier, bulk carrier, and general cargo ship, with relatively high values of   , tend to decrease as   increases, whereas for relatively slender hull types such as ro-ro/ferry and container ship dimensionless   values rather increase.This interpretation is roughly in accordance with what could be expected from knowledge of ship design principles.
Tables 1-2 show the regression equations of dimensionless   and   for each ship type estimated from Fig. 4. If the detailed hull shape or the   and   values of the ship were obtainable from the public source, they were used.Otherwise,   and   values were estimated using the proposed regression equations.It is important to note that one should be careful using regression equations as an alternative to estimating   and   of the ship, in case the input parameters are outside the range listed in Tables 1-2 or if the ship has a specific hull shape such as a bulbous bow or transom stern.In such cases, there may be gaps between the actual value and the estimated value.

Comparison of semi-empirical methods in regular waves
In this section, a comparative analysis of several semi-empirical methods for added resistance in regular waves is presented, where the methods from Boom et al. (2013), Lang and Mao (2021), and Liu and Papanikolaou (2020) are denoted as ''STA2'', ''CTH'', and ''L&P'', respectively.STA2 is also compared as a representation of the method that uses only simple ship dimensions although it is applicable only to head waves.Practically, the greater added wave resistance experienced by ships is of main interest.However, when evaluating the degree of error of the model as a residual, the greater the added resistance of the ship, the greater the residual between the predicted value and the experimental value.Therefore, mean squared error () as defined in Eq. (1), which can give more weight to a larger error by squaring the residual, is used as an evaluation metric.Here, the measurements from model experiments and estimations for the added wave resistance coefficient are compared.As can be found from the figures in the appendices (Fig. B2(a), Fig. B3(c)), there are differences among experimental results for the same ship in the same wave condition, which may be due to the experiment being carried out in different water tank environments.However, the influence of some experimental samples with relative deviations was mitigated by using as many samples as possible.The mean resistance increase of a ship in waves is influenced by many factors related to hull shape, ship operating conditions, wave characteristics, etc.In relation to the nondimensional transfer function of mean resistance increase in regular waves, it is mainly dependent on wave frequency, wave direction, and ship speed, as shown in Eq. ( 17).Therefore, as shown in Figs.5-6, the errors of each method were analyzed by classifying them into Froude number (  ), wavelengths (∕), and ship types according to the wave heading ().For the convenience of analysis, the entire wave heading area in this study is classified into three groups: head seas (180-135 degrees), beam seas (135-45 degrees), and following seas (45-0 degrees).For the detailed formula and application of each method, refer to the original documents (Boom et al., 2013;Lang andMao, 2020, 2021; Liu and  (1) where   represents the nondimensional transfer function of added wave resistance in regular waves from model test, Ĉ is the estimated value from the semi-empirical method such as STA2, CTH, L&P, and  is the number of samples for the model test belonging to the corresponding classification.  is the transfer function of added wave resistance,  is the density of water, and  is the gravity acceleration,   is wave amplitude,  is breadth of the ship, and  is the ship length.
Comparing the overall  for each method in head waves, L&P method had the smallest error.However, CTH was noticeably well estimated in the high-speed region of Froude number more than 0.25 and the short waves of less than 0.3 relative wavelengths.In particular, outstanding performance at short waves appears to be the influence of wavelength correction coefficients used in CTH method, which has been adjusted to capture an increase in resistance in the short wave region.The overall  of STA2 is about 3.5, which is a relatively larger error than the other two methods, and  comparisons as functions of   and ∕ also show no better performance than them.This is because, in contrast to the other two approaches, STA2 mainly assumes a general sea trial and does not employ information pertaining to the detailed hull shape when estimating added wave resistance.
In beam seas, overall, L&P had lower  than CTH, but similarly to the tendency in head waves, CTH had advantages in high Froude number and short wavelength ranges.In following waves, the absolute peak of added resistance was smaller than that of beam or head sea, and accordingly, the  was relatively small.Likewise, L&P showed better performance than CTH in most classifications, but CTH was better in high-speed region.
In terms of error comparison by ship type, the errors of L&P in tanker, liquefied gas, bulk, and container ships were relatively small, while CTH method prevailed in general cargo, ro-ro/ferry ships under head sea conditions.In beam and following wave conditions, except for some cases, L&P had a slightly smaller error than CTH.The difference in error between these two methods is fundamentally due to the introduction of different factors to implement the shape of added wave resistance of a ship.While CTH method modifies the peak position using an encountered frequency correction factor, and the maximum value is derived by using the amplitude adjustment factor, L&P method uses a wave heading-based trigonometric function to estimate the location of resonance frequency and maximum resistance in arbitrary waves.For a more detailed analysis, added resistance according to the wavelength of several cases are plotted in Figs.B1-B3 in Appendix B.
When the prediction results of STA2 against the experimental data in head waves are used as a benchmark, CTH and L&P methods provide significantly smaller s.In addition, since the estimations from the two methods agree fairly well with the experimental data in beam seas and following seas, both are considered to be applicable for estimation of added resistance in the environment of arbitrary waves experienced by ships at sea.
Although L&P method showed a slightly smaller  overall compared to CTH, it was not clearly a better method because they showed different performances depending on experimental conditions such as wavelength, wave heading, and ship type.In addition, it is clear that there is still much research that has to be done in this field as the experiments in beam and following sea have relatively greater uncertainty and the amount of data is limited compared to that of head sea.
Therefore, through the analysis of these existing methods, this study sought to develop a model that can ensure good overall performance at arbitrary waves without deviating significantly from model experiment data depending on ship type and various conditions.Here, we intend to apply a method that can reduce errors by properly combining the results of CTH and L&P, which is explained in detail in the next section.

Procedure for developing a combined method
As seen in the previous results, CTH and L&P methods performed relatively better in almost all comparison cases than STA2, and above all, they had the advantage of estimating results for arbitrary wave headings.Comparing CTH and L&P, the difference in performance according to ship type, wavelength, and wave heading was significant.Therefore, this study attempted to develop a new model capable of improving overall performance based on the CTH and L&P methods.Here, a meta-modeling technique was used, which is to create a new model by combining existing models (It will be denoted as a ''Combined method'' from here).The combined method basically combines the nondimensional added wave resistance coefficients estimated from CTH and L&P to minimize errors with the experimental data.Due to the lack of available model experiment data under irregular wave conditions, it was difficult to develop the model in accordance with various ship conditions.
As a method of blending the two results, a concept similar to Rfunction used in several previous papers was introduced.Fujii (1975) proposed a method of estimating added resistance due to wave reflection by applying the reflection coefficient (R-function) derived by Ursell (1947) to Havelock's formula, where the R-function was initially designed to extend the effect of wave reflection to relatively long waves.Later, this coefficient was further developed and modified by many other researchers to elaborately address the drift force due to the diffraction effect (Kuroda et al., 2008;Liu, 2020;Mourkogiannis and Liu, 2021).On the other hand, Guo and Steen (2011) adopted a method of multiplying R and 1-R by wave reflection term and ship motion term for the entire wavelength, respectively, to gain the advantage that their contribution to the added wave resistance is smoothly transited from short waves to long waves.Recently, Yang et al. (2018) adopted a more simple and practical tangent hyperbolic function as a blending function instead of a R-function composed of Bessel functions.
In this study, R-function is introduced to combine different theoretical calculations as in Guo and Steen (2011) and Yang et al. (2018), and the tangent hyperbolic function is used because of its simplicity which can smoothly connect the results of the two formulas by adjusting a few coefficients.The transition range of the R-function is extended not only to the wave frequency but also to the wave direction to enable the estimation of added resistance in arbitrary waves.Consequently, the added resistance in arbitrary wave headings can be estimated as described in Eq. ( 3).The output of the tangent hyperbolic function is between 0 and 1, which is used as the weight of the two methods for the final result.In addition, added wave resistance estimation is performed by classifying it for each ship type to reflect different characteristics caused by the hull shape.
where  ℎ ,   , and    represent the added wave resistance in head waves, beam waves, and following waves, respectively. () is a function that enables combining the various added wave resistance according to the wave headings as follows: Here, the coefficient  adjusts the slope of the tangent hyperbolic function and it is divided into  1 and  2 according to  as shown in Eq. ( 6).The coefficient  sets the intermediate position for combining the two results and is divided into 135, 45 degrees according to  value, as seen in Eq. ( 7).   and  & represent added resistance coefficients estimated from CTH and L&P.By multiplying 1- and , which are outputs of tangent hyperbolic function, by    and  & , respectively, the two results according to ∕ are smoothly connected.
The values    and  & are combined through the function (∕) given in Eq. ( 8), from which  ℎ (  or    ) can be estimated.
Here,  is a slope coefficient such as , which is divided into  1 ,  2 , and  3 as shown in Eq. ( 9).In other words,  serves to determine the slope of the function when combining the results according to ∕ using the combining function, and  is used to combine the results according to .Coefficient  represents the center position such as , which is divided into  1 ,  2 , and  3 in Eq. ( 10).

Table 3
The coefficients of combining function according to the ship type.The values in parentheses represent the interquartile ranges of 1,000 bootstrap samples for the coefficients.All the coefficients in Eqs. ( 6), (9), and (10) were tuned to minimize the error between the model test data and the estimated value from Eq. (3), and in the process, 10 cross-validations with 1000 bootstrap samplings were performed.As a result, the coefficients that provided the smallest errors were obtained through 10 cross-validations per bootstrapping, and 1000 sets of coefficients were finally obtained through 1000 bootstrap sampling.To avoid the influence of some coefficient estimations that deviate extremely from other values, the median value for 1000 bootstrap samples was adopted as the final value of the coefficient in the equation.Table 3 lists the finally obtained coefficient values for each ship type.

Results of a combined method
In this section, we show how the results of the new method calculated by substituting the coefficients of Table 3 into Eq.( 3) are actually applied and how they differ from the CTH and L&P methods.Figs.9-11 represents the  trend of the added wave resistance predictions in head waves, beam waves, and following waves by each method according to the wavelength.As can be seen from the figures, since the coefficients of the combining function are adjusted to minimize the  with model experimental data, the results of the Combined method generally tended to follow the method that provided lower  for each interval section.
In Figs.12-13, the results of the Combined method have been added to the previously covered ship cases.According to Fig. 12(a) and Fig. 12(b), the Combined method followed CTH in short waves and L&P method in long waves, and these results were well matched with the actual experimental results, and it was opposite in Fig. 12(c).The reason why different blending trends are shown here is that the coefficients of the combined function are applied differently for each ship type.Moreover, since the Combined method was tuned based on two semi-empirical methods, it had the advantage of being smoothly  connected without deviating significantly from the predicted values of the two methods.

Analysis of the combined method by experimental data in regular waves
In addition to , several error metrics with correlation coefficient are used to analyze the quality of the predicted value of the Combined method.The Pearson's correlation coefficient (), Mean Absolute Error (), and Root Mean Squared Error () are defined as in Eqs. ( 11)-( 13).
where   is the true value obtained from the experiment, ŷ is the predicted value from the semi-empirical method.Fig. 14 shows the addition of the results of the Combined method to the  comparison bar charts in Section 2.3.It can be seen that the Combined method shows a significantly smaller error compared to STA2 in head waves, and overall  is reduced compared to CTH and L&P for all wave headings.It is confirmed that the error of the Combined method is reduced compared to other methods in each section divided according to the Froude number, wavelength, and ship type.
The scatter plots between all experimental data of added resistance in regular waves and the predicted values of STA2, CTH, L&P, and Combined methods are shown in Fig. 16.The best match line of the predictions and experimental values and 30% deviation line from it     As mentioned earlier, STA2 is less accurate in head waves than CTH, L&P, and Combined methods.In particular, as many of the predicted values from STA2 for experimental measurements are located in the lower right side of the figure beyond the 30% deviation line, it is likely to underestimate when added resistance is large.This trend is also in line with the large error of STA2 at high speed and resonance positions in Fig. 14.
On the other hand, the predicted values of CTH and L&P are evenly distributed on both sides of the best match line, and the correlation coefficients are 0.82, 0.84 at head seas, and 0.68, 0.75 at beam seas, showing good correlation for both methods.In following waves, the correlation coefficients are relatively lower at 0.33, 0.59, and many predicted values are observed far outside of the 30% deviation line.This might be due to the fact that the experimental data in following sea is limited, highly uncertain, and that many values are close to zero.Most of the predicted samples from the Combined method are evenly distributed around the best match line and remain within the 30% deviation line.Compared to the predictions of CTH and L&P, the correlation coefficients of Combined method in all wave directions are higher by 0.62-0.88,and the MAEs are smaller by 0.53-0.96(see Fig. 15).

Details of full-scale measurements and weather data
In this section, a comparison between the in-service data collected from the two ships and the wave resistance estimated from the Combined method under the corresponding conditions is performed.To this end, not only waves but also the ship resistance factors in calm water, wind, and fouling and roughness conditions, which mainly account for the total resistance of the ship, are obtained through the empirical methods presented in Sections 4.3-4.5.Moreover, the errors between the estimated values of the added wave resistance and the extracted values from the in-service data are compared in Section 4.6.
Table 5 lists the main characteristics of Ship A and Ship B used in the study, and Fig. 17 shows the trajectories of the two ships for the data recording duration.The data of Ship A and Ship B were recorded continuously for 26 months and 2 months from the various fitted sensors and data acquisition systems, and the average values were stored every 15 min and every 1 min, respectively.The collected inservice data of Ship A includes 26 variables and Ship B includes 392 variables.The composition of the data from the two ships is slightly different, but the following variables were commonly used to estimate the resistance components of the ship.Navigation (GPS position, gyro heading, COG heading); Propulsion system (shaft Power, shaft rpm, shaft torque, M/E load); Operating condition (draft, trim, GPS speed, Log speed), etc.The information about the measurement methods of the ship parameters used in this study is shown in Table 6.Additionally, it was possible to extract the data of ships in sea passage operation not at berth or maneuvering by obtaining information on the voyage schedule or navigation state.
In order to calculate the added resistance in wind and waves, information on the surrounding environment the ship experiences during its voyage is required.In this study, weather information such as u and v-components of wind speed, mean wave direction, wave period, and significant wave height was obtained from the re-analysis dataset ERA5 of European Centre for Medium-Range Weather Forecasts (ECMWF), which is a global prediction model and is widely known as one of the most reliable models simulating actual sea weather conditions, instead of onboard measurement (Haiden et al., 2018).There were wind speed and direction data obtained from anemometers installed on ships, but as a result of comparison with wind data from ECMWF, it was found that some parts of the longitudinal wind speed measured from the corresponding ship are changing signs or directions without any probable cause.In addition, there was no data related to waves that could be obtained from the ship.
The horizontal resolution of the dataset is provided based on a grid of 0.25 • × 0.25 • for atmosphere and 0.5 • × 0.5 • for ocean waves, respectively, and the temporal resolution of it is hourly.The weather data at the closest position and time grid can be obtained by matching each data sample of the ship with data from the ECMWF.Through the sequential interpolations on weather data according to the location and timestamp of the ship, the actual environment encountered by the ship can be extracted.It is used as an input value for calculating added resistance, and Figs.18-19 display distributions of wind speed and significant wave height during the data collection periods.

Data pre-processing
Raw measurement data includes all types of navigation status, such as accelerating, decelerating, maneuvering, and even in port.To perform accurate speed-power performance analysis, it is necessary to extract the data sections in which the ship operates steadily.In this study, the steady-state detection algorithm proposed by Dalheim and Steen (2020) was applied.This algorithm identifies a change point among the samples by using a sliding window and its corresponding -value of the local slope.The slope of the fitted regression line from the regression analysis is used to check the unsteady state.In addition, sections with propeller speed below a certain limit were considered to be in the state of maneuvering.The voyage classification of data samples to which steady-state detection and filtering are applied is shown as an example in Fig. 20.As a result, Figs. 21-22 show the histograms of speed through water, propeller speed, mean draft, and engine power for the pre-processed data of Ship A and Ship B.

Estimation of ship resistance factor by empirical approaches
In general, the resistance components that account for most of the total resistance of a ship are calm water resistance, added resistance due to wind, and added resistance in waves.The effects of drifting and rudder on added resistance are neglected in this study.To extract the contribution of the added wave resistance to the power demand from in-service data, the resistance in calm water conditions Y. Kim et al.  and added resistance due to wind, were considered as expressed in Eqs. ( 14)-( 15).
where   is the total resistance of a ship,   is calm water resistance,   is the added resistance due to wind,   is the overall efficiency, and   is engine brake power.
Here, the frictional coefficient was obtained from ITTC-1957 correlation line, and the residual resistance coefficient was calculated using Hollenbach's method (1998) since it was found that Hollenbach's method fits relatively well for the given ship.It was developed based on regression analysis of 433 relatively modern ship models, requiring basic ship design parameters such as length, breadth, draft, displacement, wetted surface area, block coefficient, and propeller diameter for the calculation.According to the original paper, Hollenbach's method had a relatively lower standard deviation of the error in resistance against its validation test cases compared to Holtrop-Mennen (1984), Guldhammer (1974), Lap-Keller (1973), andSeries-60 (1972).
Due to changes in hull roughness caused by marine fouling, differences in powering performance may occur during the operational period, and such differences should be additionally corrected as a roughness allowance.The roughness allowance (  ) can be calculated as the difference between the total resistance coefficient of in-service data (  , ) filtered by wind speeds less than 5.5 m/s, where calm water resistance dominates, and the total resistance coefficient in the same conditions obtained from empirical methods (  , ) as shown in Eq. ( 16) (Gupta et al., 2021).  , represents the total resistance coefficient due to calm water resistance as well as the added resistance due to the wind and waves and fouling, while   , denotes the total resistance coefficient due to calm water resistance, added resistances due to wind and waves without taking fouling into account.Therefore, it is assumed that only fouling contribution is left after the subtraction.Furthermore, trends of roughness allowance were observed over cumulative static time between specific hull cleaning events, and a fitted trend line is used as a corrected roughness allowance (  ), as shown in Fig. 23(c).The overall operating speed range of the two ships used in this study was about Froude number 0.09 to 0.18, the impact of the speed of the wave-making resistance coefficient or viscous resistance coefficient within this range is not that significant, thus the corresponding resistance coefficients, which are less than Froude number 0.2 is assumed to be constant.
There were a total of four propeller and hull cleaning events of Ship A in 2.5 years, and the roughness allowance was estimated considering this.In case of the Ship B, there was no specific data related to fouling and roughness such as cleaning event history.Since the berthing time of the container ship was relatively short, and data of one voyage (approximately 2 months) was used, it was assumed that there would have been no significant change in the hull roughness during the period.As a result, the typical average hull roughness of an operating ship, 150 μm, was applied for the estimation of the roughness allowance according to MARINTEK's formula (Minsaas, 1982;Steen and Aarsnes, 2014).
The wind resistance coefficient was estimated from the regression formula developed by Fujiwara (2006), and the resistance increase due to relative wind was calculated according to the method recommended by ISO 15016 (2015).Wind affected areas of the hull were calculated based on the depth and draft of Ship A and B. Since the wind-affected area above the water surface of Ship A is mainly an accommodation area, the upper structure according to the draft of the ship was estimated using a detailed hull shape.Meanwhile, the container ship has not only accommodation areas but also cargoes on the deck, and since there was no detailed information on cargo volume and arrangement, the parameter estimation method from Kitamura et al. (2017) was used for the above-water structure area for Ship B.
Finally, the relevant resistance components of each ship obtained through the previous estimation process are shown in Figs.23-24.In the case of Ship B, the figure related to the roughness allowance is not included as it was assumed to be constant.The detailed process to estimate added wave resistance is covered in the next section.

Wave spectrum and response amplitude operator in regular waves
The mean resistance increase in short-crested irregular waves (  ) can be calculated as a linear superposition of the transfer function of added resistance in regular waves (  ) and directional wave spectrum (), as expressed in Eq. ( 17).The transfer functions of added resistance in regular waves according to different speeds for the subject ships estimated from the Combined method are shown in the following Figs.25-26.
The directional spectrum is not measured in this study, standard frequency spectrum () with the angular distribution function () is considered as in Eq. ( 18).

𝑆(𝜔) =
5 exp   where  is the circular wave frequency,   is the significant wave height, and   is the mean wave period.
For the angular distribution function for the wind waves, the cosinepower type is applied such as in Eq. ( 22), and the spreading parameter is set to 1 (ITTC, 2017) where  is the primary wave direction,  is a directional spreading parameter, and  is a Gamma function.

Comparison of semi-empirical methods in irregular waves
Theoretical estimations for added resistance in irregular waves of Ship A and Ship B by various methods are plotted in Fig. 28 and Fig. 29, respectively, in which sensitivity analysis according to significant wave height (  ), mean wave period (  ), wave heading (), and ship speed ( ) is performed.The reference conditions for the comparative case of Ship A are   = 1 m,   = 10 s,  = 180, and  = 15.5 knots, and Ship B are   = 1 m,   = 12 s,  = 180, and  = 24.7 knots.It is assumed that the remaining conditions are constant while analyzing variation for each parameter.Basically, the added resistance in irregular waves is calculated through the process of Section 4.4, and the RAOs of the two ships under the reference conditions are shown in Fig. 27.For the comparison in irregular waves, not only STA2, CTH, and L&P but also Shopera and Kreitner's methods that can directly calculate the added resistance in irregular waves are added.The CTH, L&P, and Combined methods take into account the angular distribution function in estimating added resistance in irregular waves as shown in Eq. ( 17), while for Shopera and Kreitner the angular distribution function is not applied, since these methods work directly with irregular waves, so a normal directional spreading is presumably included in the methods.Since STA2 is applicable to the mean resistance increase in long crested irregular head waves according to the original intention, the angular distribution function in Eq. ( 18) was not used.
According to the study of Lang and Mao (2021), it was found that the added resistance rose more drastically as the significant wave height increased compared to the linear superposition.In their work, a wave height correction factor (   = 3.5 √   ) was established that can be used to account for the effects of higher resistance brought on by a large vessel's motion as well as decreased propulsion efficiency in rough seas.This correction is reflected by multiplying the added resistance due to waves (  ) by the wave height correction factor (   ), and Its effectiveness in minimizing errors between full-scale measurements and estimated values from the semi-empirical method has been shown.In this study, since the Combined method basically implemented the CTH method in its original form, the wave height correction factor was applied to the CTH method also when it was used as part of the Combined method.

Wave resistance over the significant wave height
As the significant wave height increased, added resistance due to waves tended to increase very steeply.In particular, the difference in results between methods at a height of 2 m or more was noticeable.For the considered ships A and B, the results of the Combined, CTH methods, and STA2 were the largest.In both cases, the Shopera and L&P methods provided similar results, while the Kreitner provided smaller values compared to the other methods.
Wave resistance over the mean wave period STA2, CTH, L&P, and Combined method formed a peak at around 8-9 s for Ship A and 11-12 s for Ship B, and considering the actual length of the ships, it seemed to match with the resonance frequency positions properly.In addition, weak local crests were formed at around 4-6 s and 5-7 s, respectively, which was interpreted as reflecting the increase in added resistance at a short wavelength.The Combined method seemed to properly reflect the characteristics of CTH and L&P methods according to the mean wave period.Shopera and Kreitner methods showed constant results over the mean wave period, as expected.Since they focus on estimating maximum added wave resistance and do not include the mean wave period in the formula, Shopera provided similar results as the maximum values seen from CTH, L&P, Combined methods, and STA2, while Kreitner had values less than that.

Wave resistance over the wave heading
Since Shopera, Kreitner, and STA2 only consider added wave resistance in head waves, it was assumed that STA2 and Kreitner provided the same values at wave headings from 135 to 180 degrees, and from 150 degrees to 180 degrees for Shopera.Moreover, their results in beam and following waves were set to zero.According to the estimated results from CTH, L&P, and Combined method, the maximum added resistance occurred in head waves, and the magnitude became smaller as it went to the stern direction.It can be seen that the added resistance in the range between 180 degrees and 135 degrees decreased slightly, then     in the range between 135 degrees and 0 degrees dropped sharply.The Combined method seemed to follow the results of the two methods, which were properly weighted for the wave headings.In head waves, STA2, Shopera provided similar resistances to CTH, L&P, and Combined method, and Kreitner tended to underestimate significantly compared to other methods.

Wave resistance over the speed of the ship
Most methods showed a trend of almost linearly increasing added resistance in proportion to speed.CTH showed a drop in resistance at a certain speed, which was found to be due to the influence of the  2 coefficient used in wave motion-induced added resistance.Since the calculation of  2 in the CTH method is different for Froude number smaller and larger than 0.12, where block coefficient (  ) and pitch gyration (  ) were introduced in the  2 equation only at   = 0.12 or higher, it turned out that some combinations of   and   result in a discontinuous  2 as function of   .This trend can be seen from the comparison results of the speed correction factor between the two methods according to Froude number in Fig. 30.Meanwhile, the Kreitner method estimated a constant value according to the ship speed, which was larger or smaller than the other methods depending on the speed.The Combined method mainly followed the CTH method in the case of the general cargo ship and the L&P method for the container ship.As the Combined method integrates the results of the two methods, the discontinuity of added resistance as a function of   occurring in the CTH method may be visible in some cases.

Observations from comparison of the combined method with full-scale measurements
In this section, the results of the combined method were compared with using full-scale measurement data of ships A and B. As shown in Eq. ( 15), the brake power of the ship is estimated by considering the ship resistance factors, propulsive efficiency, and speed.It is compared with the main engine power measured from the shaft horsepower meter on-board.
Figs. 32 and 33 show the percentage of absolute error between the measured power (  ) and the estimated power (  ) for each parameter (  ,   , ,  ), and it is calculated according to Eq. ( 23).Fig. 31 shows an example of the confidence interval, mean line, and histogram for the error of the data samples.The shadowed range in the figure is the confidence interval of 95% mean for the samples, which represents the range of values that there is a 95% probability that the mean value of the samples falls within.The collected data of each parameter were divided at regular intervals to obtain the error defined in Eq. ( 23) for the samples for each section.In addition, the confidence interval and mean line representing each section were estimated, and they were connected.Here, ''no correction'' means that added resistance due to waves is not included in the total estimated power.

𝐸𝑟𝑟𝑜𝑟 𝑜𝑓 𝑝𝑜𝑤𝑒𝑟 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑖𝑜𝑛
Error trend over the significant wave height From Fig. 32(a) and Fig. 33(a), it can be seen that for a significant wave height of 1 m or less, there is little difference in error between the methods, including the ''no correction''.This is as expected since the added resistance is a small fraction of the total for small waves.At a significant wave height of 1 m or more, the difference between ''no correction'' and the other correction methods is clearly visible.Since there were not many data samples in the range between 3 m to 4 m, uncertainty was included in the error trend, but the relative performance difference of each method can be distinguished.It can be seen that the Combined method gives the smallest error compared to the other methods over the entire   range, with values in the range 7%-20%, while CTH and L&P show an error of about 10%-25%.The errors of Kreitner, Shopera, and STA2 methods show much larger errors than these.

Error trend over the mean wave period
The error plotted as a function of the mean wave period in Fig. 32(b) and Fig. 33(b) ranges from 8 to 26% over the entire range, which is somewhat smaller than for the other parameters.Compared with ''no correction'', the effect of wave correction can be identified for periods larger than 5 s for Ship A and 7 s for Ship B. It can also be seen that the errors of CTH, L&P, and Combined methods are significantly smaller in the vicinity of the peak compared to other methods.

Error trend over the wave heading
Since Shopera, Kreitner, and STA2 are applicable only to head waves, they have the same error as ''no correction'' in following and beam waves.The errors of CTH, L&P, and Combined methods differ significantly in the range between 60 to 135 degrees compared to other methods.The Combined method gives the smallest error over most of the heading range.However, the increase of propulsion power due to waves of Ship A and Ship B in the range of 0-30 degrees and 0-45 degrees, respectively, seems to be almost insignificant.

Error trend over the speed of the ship
Looking at the error trends over the speed in Fig. 32(d) and Fig. 33(d), the relative effect of wave correction on the propulsion power generally decreases as the speed of the ship increases, which indicates that added resistance increase less rapidly with speed than the calm water resistance.The Combined method shows a similar error trend as CTH for the general cargo ship and as L&P for the container ship and provides generally good performance for all speeds.
Fig. 34(a) shows  of power prediction by each method for the entire in-service data, and Fig. 34(b) represents the relative , of which all the  results are normalized based on the  of Combined method to identify the relative error degree of each method.The error metrics used in Fig. 34 are defined in Eqs. ( 24)-( 25).Since Fig. 28 and Fig. 29 show theoretical results of added resistances assuming specific conditions, and Figs.32 and 33 represent the error of the predicted value for the actual data in all operating conditions, the results of those figures might give slightly different levels of agreement    between the different methods.

𝑅𝑀𝑆𝐸 𝑜𝑓 𝑝𝑜𝑤𝑒𝑟 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑖𝑜𝑛
where   is the  of the selected empirical method, and   is the  of the Combined method.
It is confirmed that the correction of the propulsion power due to waves is large and strongly varying as function of  ,   , , and   within the range of the collected data used in the study.Referring to Fig. 34, Relative s on the two ships show a similar trend depending on the overall method, although there is a difference in the value.Compared to the  of the Combined method, Ship A showed about 30%-45% larger errors for Kreitner, Shopera, and STA2 and about 8%-9% for CTH and L&P, while Ship B had 14%-19% and 1%-5% larger errors.Although we did not list the results here, applying the wave height correction factor also to the L&P method was effective in reducing errors against in-service data that we collected.However, we decided not to include the wave height correction to the L&P method, since it is not part of that method as it is originally published.One should take caution when it is intended to be used in general as it is still a preliminary concept.
In the process of obtaining the added wave resistance of the ship from the in-service data, many uncertainties may be included, such as errors included in the data and estimation of various resistance components.It should be noted that while the Combined method does not include any contribution from steering and yawing to the added resistance, such effects must be expected to be present, to some extent, in the in-service data.However, the following interpretation was obtained in common through the observations of Figs.32-34.The methods only valid for head waves such as Shopera, Kreitner, and STA2 has the effect of reducing errors in the estimation of the added resistance compared to the case of ''no correction''.More complicated methods such as CTH and L&P, which reflect hull shape information and can be used in any wave headings, can significantly reduce errors.The Combined method generally has smaller errors than other semiempirical methods described here and shows noticeable performance in estimating added wave resistance at high wave height, resonance frequency, arbitrary wave headings, and relatively low ship speeds.

Conclusions
Estimating the added resistance in arbitrary waves of ships in a proper way has always been a challenging task.In this study, several Y. al. semi-empirical methods were compared using abundant public experimental data with various types, and a new meta-model was proposed combining existing semi-empirical methods.This method is developed for the calculation of added resistance for large fleets of ships, so that robustness, computational efficiency, and applicability to a range of different ship types are priorities.
From the thorough investigation against experimental data, CTH and L&P methods were chosen as a basis for the new model due to high accuracy and the availability against arbitrary wave headings.The two methods have been combined smoothly using a tangent hyperbolic function according to wavelengths and wave headings.The coefficients constituting the combining function were tuned in the direction of minimizing  between model experiments and provided for each ship type.The Combined method showed improved results without significantly deviating from the prediction range of existing methods.It has been compared with full-scale measurements of a general cargo and a container ship.For the two vessels, the errors of Kreitner, Shopera, and STA2 were larger with about 14%-45%, and CTH and L&P with 1%-9%, compared to the  of the Combined method.In particular, it was found that the estimation of added resistance in arbitrary waves was more effective in simulating the environment experienced by ships at sea than the methods considering only head seas.It also showed good performance in estimating added wave resistance in the range of high wave height, resonance frequency, arbitrary wave headings, and low ship speed.As can be seen from the comparison of models using experimental test data and full-scale measurements, the Combined method showed good overall performance in various environments.The findings suggest that the new method can be widely applied for any purposes requiring the speed-power performance under the influence of waves such as reference at the initial design stage, speed corrections in sea trials, or overall performance evaluation of a fleet, where detailed hull shape information and advanced tools are not available.
Some of the tank tests in regular waves showed that the interval between wave frequencies was too sparse, that the experiment was not sufficiently conducted in some conditions such as short waves, or that collected samples were scattered even within similar frequency ranges.Due to these problems, it was difficult to estimate added resistance in irregular waves using them.If more experimental data in irregular waves are obtained, detailed comparisons between various estimation methods will be possible, and with more experimental data for various conditions, we expect that the reliability of the model can be improved.Furthermore, since the Combined method is combining two existing methods, it is likely that the Combined method has some of the same   based on 90 degrees.In both methods, the maximum added resistance decreased as the wave heading decreased, but as explained earlier, the peak position tended to be somewhat different.

Comparison in short waves
For added resistance in short waves, wave diffraction due to bow reflection dominates, while the effect by wave radiation is almost insignificant.As can be seen from Fig. B2, STA2 has little curvature in the short waves of wavelength less than 0.3.As Yang et al. (2018) pointed out in their paper, since the reflection coefficient of STA2 becomes unity in the short wave region resulting in constant resistance coefficient, STA2 did not properly estimate an increase of added resistance in the corresponding wavelength range in our case studies.On the other hand, L&P and CTH seems to reasonably estimate the tail shape for short waves.,The CTH method was most consistent with experimental results in short waves.

Comparison at high Froude number
As illustrated in Fig. B3, STA2 was less accurate in estimating the added wave resistance of a high Froude number than other methods because the maximum resistance at the resonance frequency was underestimated or the positions of the resonance frequency did not match.Meanwhile, in accordance with the bar chart (Fig. 5), the added resistance value estimated from L&P had greater  under highspeed operating conditions than that of CTH.As a result of a closer look at the model experiments conducted at Froude number more than 0.25, it generally consisted of datasets of ships that operate in high-speed, such as container ship, general cargo, and Ro-Ro/Ferry.In the case of a container ship at a high Froude number, the two methods provided almost similar results, and even though not described here, there was also no significant difference in the case of other container ships.From such findings, it was determined that this error was not caused by the high Froude number but rather by certain ship types such as general cargo and ro-ro/ferry.

Fig. 1 .
Fig. 1.The composition of the experimental database according to various ship parameters.

Fig. 2 .
Fig. 2. Definition of length of entrance and run.

Fig. 3 .Fig. 4 .
Fig. 3. Scatter plots of (a) length of entrance and (b) length of run against ship length with the color bar showing the block coefficient.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Comparison of  results for added wave resistance methods for all ship types, Froude numbers, and wavelengths.The bar graph in the upper row shows head waves, the middle represents beam waves, and the lower is following waves.

Fig. 6 .
Fig. 6.Comparison of  results for added wave resistance methods according to the ship type.

Fig. 7
Fig. 7 displays the combining function and its coefficients in this study as examples.Fig. 7(a) presents the combining function value  ()

Fig. 7 .
Fig. 7. Coefficients and combining function values according to (a) wave headings, (b) wavelengths proposed in the study.

Fig. 8 .
Fig. 8.The results of (∕) according to the (a) slope coefficient and (b) center position coefficient.Coefficient  also has the same trend as the results shown in Fig. 8(a) as it is a slope coefficient.

Fig. 9 .
Fig. 9. Comparison of  according to wavelengths in head waves.(a) Tanker, (b) Liquefied gas, (c) Bulk carrier, (d) General cargo, (e) Container, (f) Ro-Ro/Ferry.The circle marker stands for the mean of the s of the samples in the corresponding wavelength interval, and the error bar represents the standard error of .

Fig. 14 .
Fig. 14.  results of Combined method for added wave resistance according to all cases, Froude numbers, and wavelengths.The bar graph in the upper row shows head waves, the middle represents beam waves, and the lower is following waves.

Fig. 15 .
Fig. 15. results of Combined method for added wave resistance according to the ship type.

Fig. 16 .
Fig. 16.The scatter plots of the predicted   from (a) STA2, (b) CTH, (c) L&P, (d) Combined method against experimental data.The solid line represents the best match line, and dashed line is 30% deviation line from the best match line.

Fig. 17 .
Fig. 17.Operational routes of Ship A and B for the data recording duration.

Fig. 18 .Fig. 19 .
Fig. 18.The distribution plots of (a) wind, (b) waves encountered by Ship A. Figure (a) shows the degree of wind occurrence according to the true wind direction in %.The color sector shows the true wind speed in m/s. Figure (b) represents the degree of waves occurrence according to the wave direction in %.The color sector shows the significant wave height in meters.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 20 .
Fig. 20.Classification of voyage status according to propeller speed.It shows an example of one short voyage data.

Fig. 22 .
Fig. 21.Histograms of (a) speed, (b) main engine rpm, (c) mean draft, and the (d) engine power of Ship A for the data recording duration.

Fig. 23 .
Fig. 23.Ship resistance components for Ship A: (a) calm water resistance, (b) wind resistance coefficient, (c) roughness allowance.The resistance factors in Figures (a) and (b) are calculated based on the design loading condition.The blue shaded part in (c) is propeller cleaning event, and the solid line shows the trend line of mean roughness allowance according to the ship static time.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig.
Fig. Ship resistance components for Ship B: (a) calm water resistance, (b) wind resistance coefficient.

Fig. 25 .
Fig. 25.Comparison of added resistance coefficient in regular waves of Ship A at (a)  =   , (b)  = 0.75  , (c)  = 0.5  computed using Combined method.  represents the design speed, and the results of the figures are based on the design loading condition of the ship.

Fig. 28 .
Fig. 28.Added resistance of Ship A in irregular waves against (a) significant wave height, (b) mean wave period, (c) wave heading, (d) vessel speed according to various estimation methods.The  axis represents added resistance in irregular waves, and the section indicated by the hatched line is the range beyond the restriction of the corresponding method (Reference conditions:   = 1 m,   = 10 s,  = 180,  = 15.5 knots).

Fig. 31 .
Fig. 31.An example showing the power prediction errors of collected samples with the histogram and confidence interval.It corresponds to the case of the ''no correction'' according to wave headings of Ship A (Refer to Fig.32(c)).The solid line indicates the mean line of the power prediction error, shadowed area represents its confidence interval of 95% mean, and marker is the collected data.

Fig. 32 .
Fig. 32.Power error of Ship A against (a) significant wave height, (b) mean wave period, (c) wave heading, (d) vessel speed according to various estimation methods.The  axis represents the absolute errors between the measurements and the estimation as a percentage.

Fig. 33 .
Fig. 33.Power prediction error of Ship B against (a) significant wave height, (b) mean wave period, (c) wave heading, (d) vessel speed according to various estimation methods.

Fig. 34 .
Fig. 34.Error analysis against full-scale measurements for Ship A.Figure (a) and Figure (b) show  and relative .
Fig. 34.Error analysis against full-scale measurements for Ship A.Figure (a) and Figure (b) show  and relative .
.112749 Received 27 October 2021; Received in revised form 24 September 2022; Accepted 26 September 2022 developed a formula that simplifies the resistance in long waves based on the experimental data of fast cargo ships.Two simple methods have been developed by STA-JIP to correct the added resistance in

Table 1
Regression equations for the dimensionless   according to the ship type. and  in equation represent   and   ∕  .

Table 2
Regression equations for the dimensionless   according to the ship type. and  in equation represent   and   ∕  .
, and    according to wave heading, and Fig.7(b)shows the combining function value (∕) of   and   according to wavelength.As illustrated in Fig.8, coefficients  and  affect the slope of the combining function value, and coefficient  adjusts where the combining weight is half.As the absolute value of  coefficient increases, the slope of the combining function increases ( has the same trend as ), and as the value of  increases, the center point moves in the direction where ∕ increases.

Table 4
Summary of the components of the correlation coefficients and statistical values for the predicted results in Fig.16.
Wang et al. (2021)ther in the figure.In addition, various error metrics such as ,  with  obtained from the corresponding cases are presented in Table4.For the comparison using scatterplots and the table, refer to those presented byWang et al. (2021)in a benchmark study organized by ITTC.

Table 5
Main dimensions and information of ships used in the study.

Table
Measurement methods of the ship parameters used in the study.

Table A .1
Experimental study of added resistance in arbitrary waves for tanker.Experimental study of added resistance in arbitrary waves for liquefied gas carrier.