Numerical Computation of Restoring Time and Prediction of Self-Righting Process

: Excellent self-righting performance is important to guarantee the normal navigation of Unmanned Surface Vehicles (USVs) after overturning, and the restoring time is an important index in design requirements. Traditionally, the static stability method and experiments on full-scale vehicles were used to analyze the large-angle stability of the USV. However, when it comes to the analysis of self-righting performance, the traditional static stability method will cause improper integration, and experiments are not convenient. To solve these problems, an improved static stability method was proposed, and a whole self-righting process simulation of a physical model was finished. The numerical simulation method was used to predict the self-righting process of a USV under four working conditions. Firstly, a midpoint average method based on the static stability theory was adopted to compute the static restoring time, and the results were compared with the results of the references, which verifies the effectiveness of the midpoint average method. Also, the midpoint average method is convenient because it only needs the restoring arm curve, the width and the gravity center height. Then, a numerical simulation of a physical model in static water was finished, and an experiment for a physical model in a towing tank was conducted. Comparing the restoring time of the midpoint average method, the numerical simulation and the experiment, the results show that the numerical simulation has high accuracy. Moreover, the numerical simulation was used to predict the self-righting process and analyze the self-righting performance of a USV under four working conditions.


Introduction
Under real sea conditions, the Unmanned Surface Vehicle (USV) may overturn when executing tasks.If the USV does not have the self-righting ability after overturning, it will lose maneuverability and can only wait for rescuing, which will cause problems when the USV sails out.Therefore, a USV with the self-righting ability is favored and researchers have conducted broad research on it.
At first, researchers mainly analyzed the self-righting performance via experiments; it is intuitive to improve the performance according to experimental results.Hudson et al. evaluated the design requirements of the hull form, the machinery installation and the construction materials [1].The results show that composite materials which were carefully designed could reduce high hydrodynamic influence.Shepard et al. introduced a highspeed and heavy-weather motor lifeboat, which was designed to have the self-righting ability and operate in a wave which is up to 20 feet in height [2].Renilson et al. conducted an experiment on a 1/12.5 scale model and developed two different experimental procedures which were for a single breaking wave and steep irregular waves, and then they tested four working conditions to explore the limit of stability [3].Kim et al. researched materials to prevent the overturn of a rigid inflatable boat [4].The results show that Thermoplastic Polyurethane tubes were better than single Hypalon tubes and that high buoyant materials arranged in empty spaces could increase the spare buoyancy of the ship.Fu researched the configuration of a self-righting lifeboat named Huaying 396 [5].Grenestedt et al. introduced a high-performance USV with applications for surveillance and monitoring tasks, and then they demonstrated an effective trajectory following a riverine environment at peak speeds of 30 mph [6].Akyildiz et al. used a watertight superstructure to provide enough buoyancy to realize the self-righting function [7].The results show the advantage of this structure by comparing it with airbags and self-operating ballast water replacement systems.Relying on experiments, researchers mainly tested self-righting performance under uncomplicated working conditions.
However, the model must be manufactured via complicated steps and the experiment can only research a few working conditions, which costs a lot and is not convenient, especially when it comes to large-scale ships.Furthermore, it is hard to measure some data such as restoring arms, which are important for self-righting analysis in the experiment.With the development of numerical methods and lots of commercial software, researchers began to use these new tools to analyze the self-righting performance.
Normally, the simulation of self-righting ability analyzes the static stability based on the theory of ship hydrostatics and stability.In this analysis method, the static stability curve should be computed first, and researchers should observe whether the restoring arm is always positive when the roll angle is between 0 • and 180 • .In static stability theory, the restoring arm is considered to be the criterium to analyze the self-righting performance.In order to guarantee enough static stability during the self-righting process, the restoring arm should be positive at all roll angles.Yim et al. proposed a stability analysis method under severe sea conditions and built roll-heave models with a high-degree polynomial approximation of restoring forces and moments [8].Bellec et al. applied artificial neural network technology to the parametric rolling prediction on head seas [9].Storey et al. developed a MATLAB simulator to model the vessel's tilting motion on the surface of the sea [10].Jing et al. researched the self-righting performance of a 10 m rescue boat with different trims, gravity centers and water displacements through NAPA [11].Yao et al. established a nonlinear roll equation of a rotary-molded boat based on nonlinear feature analysis of restoring torques and damping torques [12].Then, they dissected its dynamic characteristic of rolling motion via the simulation.Bai et al. computed the restoring time and rolling cycles by solving the motion equations based on the transformed variable method [13].The comparisons of computation results and the model experiment results of a float boat showed the reliability of this method.Tavakoli et al. determined the roll restoring moment by implementing an asymmetric parameter in previous lift equations and proposed an empirical method for computing restoring arms [14].Zaman et al. analyzed the effect of gravity point on the engine room longitudinally, vertically and transversely on the simulation's ability to self-right [15].Lin et al. proposed a new self-righting design method based on the parametric design of decks for an unmanned speedboat [16].Deng et al. simulated the hydrofoil rotation around the center of a wave-powered boat and the passive rotation around the center of the hydrofoil [17].Guan et al. proposed a design method of self-righting decks, and they analyzed stability based on the combination of the particle swarm optimization and the sequential quadratic programming [18].Trimulyono et al. examined the design of an anti-capsize ship by improving the self-righting moment with different heights of deck house and designed a ship whose restoring arm is positive when the roll angle is between 0 • and 180 • [19].With the improvement of algorithms and software, the efficiency of self-righting designs has greatly increased compared with traditional experiments.
As the development of self-righting designs progresses, the period of restoring time has been restricted.Restoring time has been added to the criteria of self-righting performance analysis by more and more researchers [13].According to the design requirement, restoring time refers to the time from the moment when the USV begins the self-righting motion to the moment when the USV returns to the upright position for the first time.
Xu simulated the whole self-righting process of an electric propulsion unmanned ship and computed restoring time [20].Nevertheless, the acquisition of restoring time is mainly based on static stability theory, in which the damping factor and the wave factor are not considered or are simplified via estimation.The above static stability method, which only exists according to restoring arm curves, cannot meet present design requirements which contain restoring time.
Based on the above studies, it is of great importance to improve the traditional static stability method and adopt a new numerical simulation method to predict the self-righting process of the USV.Before the simulation, a midpoint average method based on the static stability theory was used to compute the static restoring time.Even though the midpoint average method is not accurate enough under real sea conditions, its computational process is convenient because it only needs the restoring arm curve, the width of the USV and the gravity center height of the USV.Rigorous verification studies, including the Grid Convergence Index (GCI) analysis and the comparison with the experiment, were carried out to validate the simulation.A USV under four working conditions was used in CFD prediction.Comparisons between the self-righting process on four working conditions were performed to analyze the self-righting performance.
The remainder of this paper is divided into four parts.In Section 2, the disadvantage of the traditional transformed variable method is analyzed, and a new midpoint average method is introduced and validated.In Section 3, the detailed procedure of the simulation is introduced, and validation studies are conducted on the simulation.Then, the results of the simulation are compared with those of the midpoint average method and those of the experiment.In Section 4, the simulation results of a USV under four working conditions are introduced and the restoring time results are analyzed.In Section 5, there is the summary and the reflection for this paper.

Transformed Variable Method
According to the nonlinear roll theory, there are four moments during the self-righting process: the inertia moment, the damping moment, the restoring moment and the wavedisturbing moment.The equation of the nonlinear large-angle roll motion is [21]: where ϕ is the roll angle, • ϕ is the angular velocity, ϕ is the angular acceleration, I ′ xx is the total moment of inertia, 2N is the coefficient of damping moment, D is the weight of displacement, l(ϕ) is the restoring arm, h is the initial stability height, α m0 is the amplitude of the roll angle of the effective wave and ω is the angular frequency of the wave.
Usually, the damping moment term is measured in the experiment, the restoring moment term is computed via fitting the restoring arm curve and the wave-disturbing moment term is computed from the approximate estimation in engineering.
Due to no higher requirement for accuracy, traditionally, the restoring time of the USV in static water is used to represent self-righting performance under real sea conditions.When it comes to static water, the wave-disturbing moment term does not exist.Because damping moment has a slight influence on restoring time, it can be ignored for convenience.Equation (1) can be changed to: The total moment of inertia is itself composed of the moment of inertia I xx and the added moment of inertia J xx .The computation of J xx is complicated in the CFD method.So, the total moment of inertia I ′ xx is computed using an empirical equation: where B is the width of the ship, g is the gravitational acceleration and z g is the height of the center of gravity (measured from the bottom of the ship).
According to the nonlinear large-angle roll theory, Equation ( 2) is changed to: where ϕ m is the maximum roll angle, ϕ m ∈ (0 0 , 180 0 ), θ is a temporal parameter in the integration which represents the roll angle and l d (ϕ) is the dynamic stability arm, l d (ϕ) = ϕ 0 l(θ)dθ.The time from ϕ m to 0 • is: The improper integration in Equation ( 6) is changed to: where . The monotonicity of ϕ is consistent with that of ξ, which means that when ξ approaches 0 • , ϕ will approach 0 • ; when ξ approaches 90 • , ϕ will approach ϕ m .
However, when it comes to the self-righting motion, the maximum angle ϕ m is 180 • .
The restoring arm at 180 • is zero, which means lim This integration is improper and will cause a large error.

Midpoint Average Method
To avoid the improper integration, a new midpoint average method (M-A method) to compute the static restoring time was proposed.The whole self-righting angle is averagely discretized into N intervals and N + 1 nodes.In the nth interval, the variable angular accelerated motion is approximately replaced by a uniformly angular accelerated motion.The angular acceleration is the arithmetic average value of the angular acceleration in two endpoints.Equation ( 2) is discretized: where ϕ n is the average angular acceleration in the nth interval, ϕ n are angular accelerations at the n−1th node and the nth node and l(ϕ n−1 ) and l(ϕ n ) are restoring arms at the n − 1th node and the nth node.
The iterative equation of the angular velocity at the n − 1th node and at the nth node is: where ω 0 is the initial angular velocity while ω n−1 and ω n are angular velocities at the n−1th node and the nth node.
Total static restoring time is the sum of time in every interval.
The overall solution flowchart of the new midpoint average method is presented in Figure 1.

Validation and Analysis
The static restoring time of a trial boat in the transformed variable method and in the experiment were used to verify the accuracy of the midpoint average method.According to the previous study [13], there were two conditions in the experiment and four working conditions in the transformed variable method.The only difference between these conditions was the gravity center height.Restoring arm curves on different conditions were shown in the previous study.Based on these curves, the static restoring time in the midpoint average method was computed via the procedure shown in Figure 1.In this procedure, roll angle nodes and their corresponding restoring arms were straightly measured from restoring arm curves in the previous study.
Based on data from restoring arm curves in the previous study, the results of the new midpoint average method were compared with the available experimental results [13] and transformed variable method results [13].In these methods, gravity center height is measured from the bo om baseline of the boat.According to Equations ( 1) and ( 2), if the damping factor and wave-disturbing factor are ignored, the restoring time of the static condition

Validation and Analysis
The static restoring time of a trial boat in the transformed variable method and in the experiment were used to verify the accuracy of the midpoint average method.According to the previous study [13], there were two conditions in the experiment and four working conditions in the transformed variable method.The only difference between these conditions was the gravity center height.Restoring arm curves on different conditions were shown in the previous study.Based on these curves, the static restoring time in the midpoint average method was computed via the procedure shown in Figure 1.In this procedure, roll angle nodes and their corresponding restoring arms were straightly measured from restoring arm curves in the previous study.
Based on data from restoring arm curves in the previous study, the results of the new midpoint average method were compared with the available experimental results [13] and transformed variable method results [13].In these methods, gravity center height is measured from the bottom baseline of the boat.According to Equations ( 1) and ( 2), if the damping factor and wave-disturbing factor are ignored, the restoring time of the static condition should be shorter than that of the real sea condition because of larger angular acceleration.The results in Figure 2 and Table 1 show that the absolute value of comparison error of the new midpoint average method was lower than that of the transformed variable method, which indicates that the validation of the new midpoint average method was successful.
to the previous study [13], there were two conditions in the experiment and four working conditions in the transformed variable method.The only difference between these conditions was the gravity center height.Restoring arm curves on different conditions were shown in the previous study.Based on these curves, the static restoring time in the midpoint average method was computed via the procedure shown in Figure 1.In this procedure, roll angle nodes and their corresponding restoring arms were straightly measured from restoring arm curves in the previous study.
Based on data from restoring arm curves in the previous study, the results of the new midpoint average method were compared with the available experimental results [13] and transformed variable method results [13].In these methods, gravity center height is measured from the bo om baseline of the boat.According to Equations ( 1) and ( 2), if the damping factor and wave-disturbing factor are ignored, the restoring time of the static condition should be shorter than that of the real sea condition because of larger angular acceleration.The results in Figure 2 and Table 1 show that the absolute value of comparison error of the new midpoint average method was lower than that of the transformed variable method, which indicates that the validation of the new midpoint average method was successful.

Numerical Simulation
Different to the large-angle stability, there are more obvious damping factors and wave factors in the self-righting process because of the larger angular velocity and the larger change in water surface shape.The numerical simulation is necessary to take these two factors in account.

Geometry Model
To verify the validation of the simulation in static water, a physical model was used to compute the restoring time.Figure 3 shows the geometry model, and Table 2 shows the main parameters.The overall appearance of the physical model was a ladder shape, which is convenient to control its geometrical parameter.The origin point of the global coordinate was set in the center of the bottom of the trapezoid.The top length of the trapezoid was set as Lpp (length between perpendiculars).

Numerical Simulation
Different to the large-angle stability, there are more obvious damping factors and wave factors in the self-righting process because of the larger angular velocity and the larger change in water surface shape.The numerical simulation is necessary to take these two factors in account.

Geometry Model
To verify the validation of the simulation in static water, a physical model was used to compute the restoring time.Figure 3 shows the geometry model, and Table 2 shows the main parameters.The overall appearance of the physical model was a ladder shape, which is convenient to control its geometrical parameter.The origin point of the global coordinate was set in the center of the bo om of the trapezoid.The top length of the trapezoid was set as Lpp (length between perpendiculars).3 show the results of the interval quantity convergence analysis.The error between the result of 180 intervals and the result of 200 intervals is less than 1.3%, which means the accuracy of 180 intervals is enough.According to Equation ( 16) in the midpoint average method, the restoring time is 2.22 s. Figure 6 shows the result of the midpoint average method in the time domain.This restoring time result is to be the reference in the simulation.3 show the results of the interval quantity convergence analysis.The error between the result of 180 intervals and the result of 200 intervals is less than 1.3%, which means the accuracy of 180 intervals is enough.According to Equation ( 16) in the midpoint average method, the restoring time is 2.22 s. Figure 6 shows the result of the midpoint average method in the time domain.This restoring time result is to be the reference in the simulation.3 show the results of the interval quantity convergence analysis.The error between the result of 180 intervals and the result of 200 intervals is less than 1.3%, which means the accuracy of 180 intervals is enough.According to Equation ( 16) in the midpoint average method, the restoring time is 2.22 s. Figure 6 shows the result of the midpoint average method in the time domain.This restoring time result is to be the reference in the simulation.

Simulation Procedure
The self-righting motion of the physical model is part of a roll decay motion, and the roll decay simulation was conducted using commercial CFD software Siemens Star-CCM+ version 17.04.The processor was Intel (R) Core (TM) i5-8250U (1.60 Ghz).

Computational Domain and Mesh Generation
The computational domain is shown in Figure 7, which extends over the ranges of −2Lpp < x < 2Lpp, −Lpp < y < Lpp, −Lpp < z < 2Lpp.The origin point of the global coordinate was set in the center of the bottom of the trapezoid.The initial state of the physical model was inverted and completely immersed under the water surface to simulate extreme conditions after capsizing.Figure 8 shows two refinement regions: the overlap refinement region and the free surface refinement region.The overlap refinement is employed between the surface of the physical model and the background to ensure that the transition of these two regions is smooth.The size of the box for the overlap refinement is −0.

Computational Domain and Mesh Generation
The computational domain is shown in Figure 7, which extends over the ranges of −2Lpp < x < 2Lpp, −Lpp < y < Lpp, −Lpp < z < 2Lpp.The origin point of the global coordinate was set in the center of the bo om of the trapezoid.The initial state of the physical model was inverted and completely immersed under the water surface to simulate extreme conditions after capsizing.Figure 8 shows two refinement regions: the overlap refinement region and the free surface refinement region.The overlap refinement is employed between the surface of the physical model and the background to ensure that the transition of these two regions is smooth.The size of the box for the overlap refinement is −0.3 m < x < 0.3 m, −0.2 m < y < 0.2 m, −0.26 m < z < 0.09 m, and this size covers the physical model height.For the free surface refinement, the size is −1.25 m < x < 1.25 m, −1.1 m < y < 1.1 m, −0.02 m < z < 0.02 m.The free surface refinement region is located near the water surface and used to distinguish the water phase and the air phase with high precision.Two mesh tools called the surface remesher and trimmed cell mesh are used to generate the mesh, which is a mixed mesh mainly composed of structured mesh.The prism layer around the physical model surface is divided into five layers of prismatic cells.The free surface refinement region includes 50 grid points in the range of the wavelength [22].For Grid Convergence Index (GCI) analysis, three types of mesh with different accuracies are generated, which are called the coarse mesh, the medium mesh and the fine mesh.The fine mesh has a total of 7,551,460 cells, and it is coarsened with the refinement ratio 2 to generate the medium mesh.The medium mesh has a total of 2,965,373 cells, and it is coarsened with the refinement ratio 2 to generate the coarse mesh.The coarse mesh has a total of 1,091,170 cells.Two mesh tools called the surface remesher and trimmed cell mesh are used to generate the mesh, which is a mixed mesh mainly composed of structured mesh.The prism layer around the physical model surface is divided into five layers of prismatic cells.The free surface refinement region includes 50 grid points in the range of the wavelength [22].For Grid Convergence Index (GCI) analysis, three types of mesh with different accuracies are generated, which are called the coarse mesh, the medium mesh and the fine mesh.The fine mesh has a total of 7,551,460 cells, and it is coarsened with the refinement ratio √ 2 to generate the medium mesh.The medium mesh has a total of 2,965,373 cells, and it is coarsened with the refinement ratio √ 2 to generate the coarse mesh.The coarse mesh has a total of 1,091,170 cells.

Physical Conditions
As for Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations, the Finite Volume Method (FVM) was used to discretize their integral form.In convective terms, a second-order convection scheme was used.The time-domain solution used a first-order temporal discretization for the convergence rate.The turbulence effect in boundary layers was modeled using the Shear-Stress Transport (SST) method [23], which uses a k − ω model near the wall and a k − ε model in the far field.
The VOF (Volume of Fluid) method was used to capture the free surface of the water and air [24].In the VOF method, two phases of the fluid (water and air) are defined by assigning a scalar value of 0 to the air and 1 to the water in each cell, while the value of the interface between the two types of fluid (air and water) is 0.5.A derived surface is set as the interface to simulate the water surface by using the VOF model.Considering the disturbance of the physical model on the water surface in the self-righting process, the flat wave is set in the VOF model.
The back side of the computational domain surface shown in Figure 9 is set as the pressure-outlet face, on which the pressure is defined by the hydrostatic pressure of fluid.The other faces shown in Figure 9 are set as the velocity-inlet face, on which the velocity is defined by the VOF fluid.
As for Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations, the Finite Volume Method (FVM) was used to discretize their integral form.In convective terms, a second-order convection scheme was used.The time-domain solution used a first-order temporal discretization for the convergence rate.The turbulence effect in boundary layers was modeled using the Shear-Stress Transport (SST) method [23], which uses a k   model near the wall and a k   model in the far field.The VOF (Volume of Fluid) method was used to capture the free surface of the water and air [24].In the VOF method, two phases of the fluid (water and air) are defined by assigning a scalar value of 0 to the air and 1 to the water in each cell, while the value of the interface between the two types of fluid (air and water) is 0.5.A derived surface is set as the interface to simulate the water surface by using the VOF model.Considering the disturbance of the physical model on the water surface in the self-righting process, the flat wave is set in the VOF model.
The back side of the computational domain surface shown in Figure 9 is set as the pressure-outlet face, on which the pressure is defined by the hydrostatic pressure of fluid.The other faces shown in Figure 9 are set as the velocity-inlet face, on which the velocity is defined by the VOF fluid.The origin point is set at the gravity center of the model.The self-righting motion can be simplified as the combination of translation along the Z-direction and rotation around the X-direction.The dynamic fluid body interaction (DFBI) module was used to capture these two motions.This module obtains a new position in each time step by computing exciting forces and moments in governing equations.The local coordinate system moves and rotates together with the model all the time during the entire motion process.

Time Step
The period of the roll decay motion after returning to the upright position is approximated using the natural roll period.Based on Equation (4), the natural roll period is as follows [21]: It is recommended that at least 100 time-steps per period are used for the periodic phenomena [25].According to Equation (13) and Table 2, the natural roll period is 1.03 s, and the maximum time step based on the ITTC recommendation is 0.01 s.For time step convergence analysis, three types of time step are generated with the division ratio 2 for coarse (0.01 s), medium (0.005 s) and fine (0.0025 s).The total simulation time is 10 s, which contains at least four roll decay periods.The results of the angular velocity in upright position simulated via coarse, medium and fine time steps are summarized in Table 4. Considering computational accuracy and efficiency, the medium time step (0.005 s) was chosen.

Grid Convergence Index Analysis
Based on the Grid Convergence Index (GCI) analysis method [26][27][28][29], the refinement was to choose the proper quantity of mesh cells, which is shown as follows: where F s is a factor of safety, ε is the relative difference, f 1 is the interested physical value in fine mesh, f 2 is the interested physical value in coarse interval, r is the mesh refinement factor, h 1 is the spacing value in coarse mesh, h 2 is the spacing value in fine mesh and p is the accuracy order of numerical scheme.In this Grid Convergence Index analysis, F s = 1.25 for a three (or more)-grid convergence study.The roll amplitudes of roll decay simulations for the fine, medium and coarse mesh configurations are shown in Figure 10.Fourier series (FS) analysis was used to capture the characteristics in the frequency domain for GCI analysis: where ξ i is the nth-order harmonic amplitude and φ i is the corresponding phase.
right position, which may lead to errors in the FS analysis.As such, the roll amplitude that is used for FS analysis is taken from the roll angle of 0°. Figure 10 shows the time series of roll amplitude for FS analysis, in which the references of time history align.Figure 11 shows the results of frequency domain after FFT, with the peak of each curve being the roll amplitude occurring at the natural frequency.
(a) (b)  GCI analysis shown in Equations ( 23)-( 27) was used as a verification study to quantify the numerical uncertainties of the simulation.This method was based on a convergence verification procedure [30].The GCI result of the fine mesh configuration shown in Table 5 is less than 4%, which indicates that the numerical error in this simulation is small.The verification of fine mesh was successful, and the restoring time in the simulation was 2.575 s.The Fast Fourier transform (FFT) was used to process the time histories and obtain the values, as follows: where T is the period of the time history, ξ 0 is the averaged value of the time history and ξ i is the linear term defined as a fundamental component.
In time history curves, there is a large difference from the initial position to the upright position, which may lead to errors in the FS analysis.As such, the roll amplitude that is used for FS analysis is taken from the roll angle of 0 • .Figure 10 shows the time series of roll amplitude for FS analysis, in which the references of time history align.Figure 11 shows the results of frequency domain after FFT, with the peak of each curve being the roll amplitude occurring at the natural frequency.where T is the period of the time history, 0  is the averaged value of the time history and i  is the linear term defined as a fundamental component.
In time history curves, there is a large difference from the initial position to the upright position, which may lead to errors in the FS analysis.As such, the roll amplitude that is used for FS analysis is taken from the roll angle of 0°. Figure 10 shows the time series of roll amplitude for FS analysis, in which the references of time history align.Figure 11 shows the results of frequency domain after FFT, with the peak of each curve being the roll amplitude occurring at the natural frequency.GCI analysis shown in Equations ( 23)-( 27) was used as a verification study to quantify the numerical uncertainties of the simulation.This method was based on a convergence verification procedure [30].The GCI result of the fine mesh configuration shown in Table 5 is less than 4%, which indicates that the numerical error in this simulation is small.The verification of fine mesh was successful, and the restoring time in the simulation was 2.575 s.GCI analysis shown in Equations ( 23)-( 27) was used as a verification study to quantify the numerical uncertainties of the simulation.This method was based on a convergence verification procedure [30].The GCI result of the fine mesh configuration shown in Table 5 is less than 4%, which indicates that the numerical error in this simulation is small.The verification of fine mesh was successful, and the restoring time in the simulation was 2.575 s.
where S 1 , S 2 and S 3 are the maximum amplitude in frequency domain roll responses obtained from the fine, medium and coarse configurations, R is the convergence ratio, p is the order of accuracy, r 21 is the refinement ratio between the medium configuration and the fine configuration and r 32 is the refinement ratio between the coarse configuration and the medium configuration.

Experiment and Validation
The self-righting experiment was performed to measure the restoring time under real conditions and validate the simulation result.The physical model is shown in Figure 12 and its main parameters are listed in Table 6.The main material of the model was wood because it is normal for boat models and it has enough intensity.There was a rectangular hole through the center of the model.The iron block was used to increase the draft in order to imitate the high-loaded working conditions.There was an iron sheet fixed in the bottom of the model.The iron sheet was used to decrease the gravity center height of the model in order to guarantee enough stability.To prevent water from penetrating the wood, there was a layer of wax wrapping the surface of the whole model.
where 1 S , 2 S and 3 S are the maximum amplitude in frequency domain roll responses obtained from the fine, medium and coarse configurations, R is the convergence ratio, p is the order of accuracy, 21 r is the refinement ratio between the medium configuration and the fine configuration and 32 r is the refinement ratio between the coarse configuration and the medium configuration.

Experiment and Validation
The self-righting experiment was performed to measure the restoring time under real conditions and validate the simulation result.The physical model is shown in Figure 12 and its main parameters are listed in Table 6.The main material of the model was wood because it is normal for boat models and it has enough intensity.There was a rectangular hole through the center of the model.The iron block was used to increase the draft in order to imitate the high-loaded working conditions.There was an iron sheet fixed in the bo om of the model.The iron sheet was used to decrease the gravity center height of the model in order to guarantee enough stability.To prevent water from penetrating the wood, there was a layer of wax wrapping the surface of the whole model.A high-speed camera was used to record the short-time process; the performance parameters of it are listed in Table 7.The self-righting process of the experiment is shown in Figure 13.The restoring time in the experiment was 2.46 s, which is shown in Figure 14.Compared with the restoring time in the midpoint average method, the restoring time of the simulation and the restoring time of the experiment were longer.Damping factor and wave factor are ignored in the midpoint average method, so when they are at the same position, the total angular acceleration in the midpoint average method is larger than that  The restoring time in the experiment was 2.46 s, which is shown in Figure 14.Compared with the restoring time in the midpoint average method, the restoring time of the simulation and the restoring time of the experiment were longer.Damping factor and wave factor are ignored in the midpoint average method, so when they are at the same position, the total angular acceleration in the midpoint average method is larger than that in the simulation and in the experiment.Larger angular acceleration at the same roll angle means less time, so the static restoring time is less than the restoring times of the simulation and the experiment.Compared to the results of the midpoint average method, the simulation results had higher accuracy.The error of restoring time between the midpoint average method and the experiment was less than 10%.The error of restoring time between the simulation and the experiment was less than 5%.The simulation can be used to predict the self-righting process of the USV.

Prediction of Self-Righting Performance
This USV was designed to execute rescue tasks and should have safe navigation ability under severe sea conditions.Parameters of the USV are listed in Table 8, and a simplified model is shown in Figure 15.There are four working conditions when the USV executes tasks: they are leaving port with the design load (condition 1), reaching port with the design load (condition 2), leaving port with maximum load (condition 3) and reaching port with maximum load (condition 4).The differences in the four working conditions can be simplified as differences in displacement and in gravity center position.The reference point is set at the backend of the bo om of the USV.The parameters of the four working conditions are listed in Table 9.The numerical simulation method above was used to predict the self-righting performance of the USV under four working conditions.

Prediction of Self-Righting Performance
This USV was designed to execute rescue tasks and should have safe navigation ability under severe sea conditions.Parameters of the USV are listed in Table 8, and a simplified model is shown in Figure 15.There are four working conditions when the USV executes tasks: they are leaving port with the design load (condition 1), reaching port with the design load (condition 2), leaving port with maximum load (condition 3) and reaching port with maximum load (condition 4).The differences in the four working conditions can be simplified as differences in displacement and in gravity center position.The reference point is set at the backend of the bottom of the USV.The parameters of the four working conditions are listed in Table 9.The numerical simulation method above was used to predict the self-righting performance of the USV under four working conditions.

Prediction of Self-Righting Performance
This USV was designed to execute rescue tasks and should have safe navigation ability under severe sea conditions.Parameters of the USV are listed in Table 8, and a simplified model is shown in Figure 15.There are four working conditions when the USV executes tasks: they are leaving port with the design load (condition 1), reaching port with the design load (condition 2), leaving port with maximum load (condition 3) and reaching port with maximum load (condition 4).The differences in the four working conditions can be simplified as differences in displacement and in gravity center position.The reference point is set at the backend of the bo om of the USV.The parameters of the four working conditions are listed in Table 9.The numerical simulation method above was used to predict the self-righting performance of the USV under four working conditions.The setting of the computational domain is the same as the introduction in Section 3. The length of it is 45 m, the width of it is 21.2 m and the height of it is 30 m, which is shown in Figure 16.The size of the box for the overlap refinement is −4 m < x < 11 m, −4 m < y < 4 m, −5 m < z < 4 m, which is shown in Figure 17.Similar to the simulation procedure in Section 3, the total cell number of the computational domain is 4,268,280 cells and the time step is 0.005 s.The displacement in condition 1 was smaller than that of condition 3, and the gravity center height of condition 1 was higher than that of condition 3. The displacement in condition 2 was smaller than that of condition 4, and the gravity center height of condition 2 was higher than that of condition 4.There were two factors which were different at the same time, and the influences of the two factors on self-righting performance were opposite.According to the restoring arm curves in Figure 18, the restoring arm is less than zero between 126.4° and 146° for condition 2, and the restoring arm is less than zero between 130° and 142.8° for condition 2. The restoring arm curves mean that, based on the static stability theory, the USV should not have the self-righting ability under condition 2 and condition 4.However, the simulation results shown in Figures 19-22 prove the self-righting ability under four working conditions.The change of water surface shape makes the   The displacement in condition 1 was smaller than that of condition 3, and the gravity center height of condition 1 was higher than that of condition 3. The displacement in condition 2 was smaller than that of condition 4, and the gravity center height of condition 2 was higher than that of condition 4.There were two factors which were different at the same time, and the influences of the two factors on self-righting performance were opposite.According to the restoring arm curves in Figure 18, the restoring arm is less than zero between 126.4° and 146° for condition 2, and the restoring arm is less than zero between 130° and 142.8° for condition 2. The restoring arm curves mean that, based on the static stability theory, the USV should not have the self-righting ability under condition 2 and condition 4.However, the simulation results shown in Figures 19-22 prove the self-right- Head wave simulations were carried out to represent a severe sea condition in a wave with wave amplitude A L pp = 0.0183 and wavelength λ L pp = 1, which follows INSEAN captive tests [31].
The displacement in condition 1 was smaller than that of condition 3, and the gravity center height of condition 1 was higher than that of condition 3. The displacement in condition 2 was smaller than that of condition 4, and the gravity center height of condition 2 was higher than that of condition 4.There were two factors which were different at the same time, and the influences of the two factors on self-righting performance were opposite.According to the restoring arm curves in Figure 18, the restoring arm is less than zero between 126.4 • and 146 • for condition 2, and the restoring arm is less than zero between 130 • and 142.8 • for condition 2. The restoring arm curves mean that, based on the static stability theory, the USV should not have the self-righting ability under condition 2 and condition 4.However, the simulation results shown in Figures 19-22 prove the self-righting ability under four working conditions.The change of water surface shape makes the restoring arm in the simulation different from that of the static stability method, which will influence the restoring time.In terms of the restoring time shown in Table 10, the self-righting performance sorted from best to worse is as follows: condition 4, condition 3, condition 1 and condition 2.

Conclusions
In our study, the self-righting process and the restoring time of a physical model were computed by using the midpoint average method, the numerical simulation and the experimental method.Then, the reliability of the simulation was fully verified by comparing the CFD results with the results of the improved static method and the experimental method.Eventually, the self-righting performance of a USV under four working conditions was evaluated and analyzed via the simulation method.It was demonstrated that: (1) The new midpoint average method has computational stability in the integration and is able to estimate the restoring time more conveniently.
(2) Compared with the experiment, the numerical simulation is accurate enough in the self-righting process.This numerical simulation method can be used to compute the restoring time precisely and to predict the self-righting process.
(3) The simulation results of a USV under four working conditions show that under condition 3 (leaving port with maximum load) and under condition 4 (reaching port with maximum load), the USV can return to upright position in 11 s.Under condition 1 (leaving port with the design load) and under condition 2 (reaching port with the design load), the USV returns to the upright position after more than 11 s.The self-righting performance with maximum load is better than that with the design load.
The new midpoint average method and the numerical simulation will serve as tools for the prediction of self-righting performance under different working conditions.However, there are also external forces such as winds and rain under real sea conditions which may greatly influence the self-righting performance of the USV.Future studies will add external forces to the CFD prediction.

Figure 1 .
Figure 1.Solution flowchart for the midpoint average method.

Figure 1 .
Figure 1.Solution flowchart for the midpoint average method.

Figure 2 .
Figure 2. Restoring time results of the experiment [13], T-V method [13] and M-A method.Figure 2. Restoring time results of the experiment [13], T-V method [13] and M-A method.

Figure 2 .
Figure 2. Restoring time results of the experiment [13], T-V method [13] and M-A method.Figure 2. Restoring time results of the experiment [13], T-V method [13] and M-A method.

Table 1 .
Results of validation compared to Bai et al. [13].Gravity center height (m

3. 2 .
Static Stability ResultThe restoring arm curve was computed in a commercial software, Maxsurf version 20.1.0.59.In Maxsurf, the reference point is in the center of the model bottom and the model is divided into 200 stations along the longitudinal direction.The result of restoring arm curve is shown in Figure4.When the roll angle was from 0 • to 180 • , the restoring arm was always positive.When the roll angle is 102.7 • , the restoring arm reaches the largest value, which is 17.5 mm.From the curve, it is clearly demonstrated that the stability range of the model is 180 • , which means the physical model has the ability to return back to the upright position.

Figure 4 .
Figure 4. Restoring arm curve of physical model.

Figure 5
Figure 5 and Table3show the results of the interval quantity convergence analysis.The error between the result of 180 intervals and the result of 200 intervals is less than 1.3%, which means the accuracy of 180 intervals is enough.According to Equation (16) in the midpoint average method, the restoring time is 2.22 s.Figure6shows the result of the midpoint average method in the time domain.This restoring time result is to be the reference in the simulation.

Figure 4 .
Figure 4. Restoring arm curve of physical model.

Figure 5
Figure 5 and Table3show the results of the interval quantity convergence analysis.The error between the result of 180 intervals and the result of 200 intervals is less than 1.3%, which means the accuracy of 180 intervals is enough.According to Equation (16) in the midpoint average method, the restoring time is 2.22 s.Figure6shows the result of the midpoint average method in the time domain.This restoring time result is to be the reference in the simulation.

Figure 4 .
Figure 4. Restoring arm curve of physical model.

Figure 5
Figure 5 and Table3show the results of the interval quantity convergence analysis.The error between the result of 180 intervals and the result of 200 intervals is less than 1.3%, which means the accuracy of 180 intervals is enough.According to Equation (16) in the midpoint average method, the restoring time is 2.22 s.Figure6shows the result of the midpoint average method in the time domain.This restoring time result is to be the reference in the simulation.

Figure 6 .
Figure 6.Restoring arm curve of physical model.

Figure 6 .
Figure 6.Results of the midpoint average method in time domain.
3 m < x < 0.3 m, −0.2 m < y < 0.2 m, −0.26 m < z < 0.09 m, and this size covers the physical model height.For the free surface refinement, the size is −1.25 m < x < 1.25 m, −1.1 m < y < 1.1 m, −0.02 m < z < 0.02 m.The free surface refinement region is located near the water surface and used to distinguish the water phase and the air phase with high precision.J. Mar.Sci.Eng.2024, 12, x FOR PEER REVIEW 9 of 21 The self-righting motion of the physical model is part of a roll decay motion, and the roll decay simulation was conducted using commercial CFD software Siemens Star-CCM+ version 17.04.The processor was Intel (R) Core (TM) i5-8250U (1.60 Ghz).

Figure 8 .
Figure 8. Two refinements and the fine mesh configuration.

Figure 8 .
Figure 8. Two refinements and the fine mesh configuration of the physical model.

Figure 10 .
Figure 10.Comparison of (a) whole time histories and (b) processed time histories in fine, medium and coarse mesh.

Figure 11 .
Figure 11.Roll amplitude in frequency domain from Fourier series analysis.

Figure 10 .
Figure 10.Comparison of (a) whole time histories and (b) processed time histories in fine, medium and coarse mesh.

Figure 10 .
Figure 10.Comparison of (a) whole time histories and (b) processed time histories in fine, medium and coarse mesh.

Figure 11 .
Figure 11.Roll amplitude in frequency domain from Fourier series analysis.

Figure 11 .
Figure 11.Roll amplitude in frequency domain from Fourier series analysis.

Figure 14 .
Figure 14.Restoring time of the midpoint average method, the simulation and the experiment.

Figure 14 .
Figure 14.Restoring time of the midpoint average method, the simulation and the experiment.

J 21 Figure 14 .
Figure 14.Restoring time of the midpoint average method, the simulation and the experiment.

J
. Mar. Sci.Eng.2024, 12, x FOR PEER REVIEW 17 of 21The se ing of the computational domain is the same as the introduction in Section 3. The length of it is 45 m, the width of it is 21.2 m and the height of it is 30 m, which is shown in Figure16.The size of the box for the overlap refinement is −4 m < x < 11 m, −4 m < y < 4 m, −5 m < z < 4 m, which is shown in Figure17.Similar to the simulation procedure in Section 3, the total cell number of the computational domain is 4,268,280 cells and the time step is 0.005 s.Head wave simulations were carried out to represent a severe sea condition in a wave with wave amplitude

Figure 17 .
Figure 17.Two refinements and the fine mesh configuration.

Figure 17 .
Figure 17.Two refinements and the fine mesh configuration.

Figure 17 .
Figure 17.Two refinements and the fine mesh configuration of the boat.

Figure 18 .
Figure 18.(a) Restoring arm curves for condition 1 and condition 2. (b) Restoring arm curves for condition 3 and condition 4.

Table 2 .
Main parameters of model.
Top length of trapezoid 48.2 cm Bottom length of trapezoid 32.24 cm Width of trapezoid 20 cm

Table 2 .
Main parameters of model.

Table 3 .
Results of interval quantity convergence.

Table 3 .
Results of interval quantity convergence.

Table 4 .
Results of time step convergence.

Table 6 .
Main parameters of model in the experiment.
Top length of trapezoid 48.2 cm Bottom length of trapezoid 32.24 cm Width of trapezoid 20 cm

Table 7 .
Parameters of camera.

Table 8 .
Main parameters of the USV.

Table 8 .
Main parameters of the USV.

Table 8 .
Main parameters of the USV.

Table 9 .
Parameters of working conditions.

Table 10 .
Restoring time for four working conditions in the simulation.

Table 10 .
Restoring time for four working conditions in the simulation.

Table 10 .
Restoring time for four working conditions in the simulation.