Numerical Simulations for the Prediction of Wave Forces on Underwater Vehicle using 3D Panel Method Code

: This study describes a method useful for predicting the hydrodynamic forces on fully appended underwater vehicle when operating beneath the free surface waves. The hydrodynamic forces acting on the underwater vehicle due to waves are sum of the incident wave force, diffraction force and radiation forces. Ansys AQWA which is based on Panel Method is selected to compute wave forces both in frequency domain and time domain. Series of numerical simulations are performed considering the 2nd order Stokes wave theory and head seas conditions. The water depth is assumed infinite hence the interference effect on vertical boundaries and horizontal bottom can be neglected. The simulation procedure requires generating panels on the vehicle's surface, defining the mass properties and providing the desired the sea wave condition. In frequency domain analyses, wave forces and moments are determined at different submergence depths and wave frequencies. These analyses reveal that magnitude of surge force, heave force and pitch moment are the most dominant as compared to other forces and moments. Numerical results in the time domain are also collected in order to investigate the effect of various factors such as wave amplitude, frequency, vehicle depth, orientation, etc., on these forces and moment. The results are then used to derive analytical formulations for wave forces and moment by curve fitting approach. The accuracy of the formulations is ensured by verifying the results with real time simulations using two test conditions in which the different values of parameters are selected. The proposed method brings great convenience to assess the maneuverability characteristics of underwater vehicles operating near sea waves.


INTRODUCTION
From oceanographic surveys to naval operations, the importance of Autonomous Underwater Vehicle (AUV) is increasing with the development of more advanced processing capabilities and robotic technology.AUVs are now considered to be suitable for roles and missions in all water depths, ranging from near-shore to deep oceans (Curtin et al., 2005).While operating in deep underwater, the predominant source of disturbance to the AUV may be the ocean currents.However, when the AUV is sailing near the sea surface in a mission, the effects due to surface waves become the most important (Jensen, 1995;Soylemez, 1996;Xia and Krokstad, 2001).The pressure field of the waves make AUV to experience additional forces, which cause bumps and rocking, thus affecting its maneuverability.Since AUVs are meant to be able to explore a new territory without human control and come back to the position at which it was launched into the water.Therefore, a capable control system that can anticipate and adapt to its surrounding is crucial which in turns necessitate the accurate estimation wave forces acting on AUV.
A number of notable works have been accomplished over last few decades on predicting the forces and moments on surface and submerged platforms due to waves.Wiley (1994) developed a simple model for predicting wave disturbances using linear wave theory, slender body theory and linear time invariant systems theory for the stationary AUV.The wave model is validated by conducting experiments on a stationary, slender body underwater vehicle in inertia dominated wave and force regime, in head and beam seas.Strip theory presented in "Ship Motions and Sea Loads" by Salvesen et al. (1970) has been developed over the past 50 years and used for predicting the unsteady sea-keeping forces of ships.With some modifications and refinements (Rybka, 2005;Milgram, 2007), strip theory is also applied for determining the response of streamlined AUV subjected to waves (Du et al., 2008;Fang et al., 2006).The theory is based on the assumptions of potential flow, slender body and small amplitude motions.While computationally efficient, strip theory cannot give accurate hull pressure predictions and has deteriorating accuracy as slenderness of the vessel decreases.
The pioneered work by Hess and Smith (1967) in the development of Boundary Element Methods (BEM) in aeronautics, laid down the foundation for most of the later work for BEM applications to ship flows.With the application of just a source distribution, the body surface is divided into N flat quadrilateral over which the source strength is assumed to be constant.The boundary conditions are satisfied at each quadrilateral centre (control point) which resulted in N linear equation for unknown source strength.By calculating the source strength, the velocities and pressure at each control point can be easily determined.The flat quadrilaterals are called "panels" and now this method is also known as "Panel Method".With the arrival of fast computers, fully three-dimensional commercial codes based on Panel Method have been developed and extensively used for predicting the response of surface and submerged platforms with reasonable accuracy (Sabra, 2003).The present study utilizes the 3D panel method code, Ansys AQWA, to predict the wave forces and moments experienced by an axisymmetric small AUV having control fins in cruciform stern configuration.ANSYS AQWA software is an engineering analysis suite of tools for the investigation of the effects of wave, wind and current on floating and fixed marine structures (Cheetham et al., 2007;Colby et al., 2011).Series of numerical analyses are performed both in frequency and time domain considering the 2 nd order Stokes wave theory and head seas conditions.The variation of the wave forces and moments with changes in wave amplitude, frequency, vehicle depth and orientation are thoroughly investigated.The results are then used to derive analytical formulations by MATLAB curve fitting function for estimating the wave forces and moments.

MATHEMATICAL DESCRIPTION
Problem definition: As illustrated in Fig. 1, Consider an underwater vehicle at depth h traveling at a constant forward velocity U in response to incoming regular waves with small amplitude ܽ as compared to its wave length ߣ.The fluid domain is bounded by the free surface and the body boundary denoted respectively as ‫ܨ‬ and S. To formulate the problem, two orthogonal and right-handed coordinate systems are introduced.The first is the space fixed system XYZ with its origin O lies in undisturbed free surface plane XY, X-axis in forward velocity direction and Z-axis vertically upwards against the gravitational acceleration g.The space fixed coordinate system facilitates the description of boundary conditions on the free surface.In order to describe rigid body motion, body fixed coordinate system xyz is used such that its origin G lies at body's centre of Gravity (CG).The axes through CG will initially be parallel to the space fixed axes.The body motions in six degree of freedom are denoted as ‫ݔ‬ (1) where, u,v and w are the components of fluid velocity vector ܸ ሬԦ .If the fluid motion is irrotational, then the problem may be formulated in terms of "Potential Flow Theory".This means that the fluid velocity vector ܸ ሬԦ can be represented by the gradient of a scalar function, i.e., velocity potential Φ as: The velocity potential Φ is related to flow velocity components as: (3) Thus, the continuity equation can be expressed in Laplace equation as: Assuming the small oscillatory motion of the vehicle due to incoming waves, the velocity potential can be decomposed into two parts as: The first part ߶ ത is the steady flow associated with the forward motion of the body in calm water and second part ߶ ෨ represents the unsteady flow due to the disturbances of the incident waves, diffraction waves and radiation waves caused by the presence and motions of the body in the steady field.
The steady potential ߶ ௦ can be written as: where, the term ‫ݔܷ−‬ accounts for the uniform current and ߶ ௦ is the steady disturbance due to presence of body.
The unsteady potential ߶ ෨ can be decomposed into three separate components as incident potential ߶ ூ , diffraction potential ߶ and radiation potential ߶ ோ .These unsteady components are harmonic in time with frequency of encounter ߱ determined as: where, ݇ = 2ߨ ߣ ൗ = ߱ ଶ /݃ is the deep water approximation of wave number ߱ is the frequency of incident waves, ߣ is the wave length and ߯ is the heading angle between the vessel and the wave propagation direction (߯ = 0 corresponds to following waves while ߯=ߨ corresponds to head waves).The unsteady potential ߶ ෨ can be written as: The radiation potential߶ ோ is further divided into six oscillatory motion components as: where ߶ the radiation wave potential in j th mode of motion, having amplitude ‫ݔ‬ .

Boundary conditions:
According to equations 4 and 5, the individual potentials should also satisfy the Laplace equation in the fluid domain, therefore: Due to elliptic in nature, Laplace equation has various solutions.In order to obtain the exact solution of a given problem, some prescribed boundary conditions are required to be defined in the fluid domain.Since the unsteady motions are supposed to be small and wave amplitude is also small as compared to wave length, linearized theory is applied for present study.Based on the theory, the boundary conditions for the potentials are as following: • Body boundary condition: For zero flux across the solid boundary of the body S having specified motion, the steady potential satisfies ∇߶ ത .݊ ሬԦ = 0 which implies that: where, ݊ ሬԦ = ݊ ଵ , ݊ ଶ , ݊ ଷ is the normal vector pointing out of the fluid domain.This condition is known as the Neumann condition.The body boundary conditions for unsteady potentials can be expressed as: The boundary condition given by Eq. ( 14) is identical to the condition derived by Timman and Newman (1962).In the equation, the n-term represents the normal velocity of the body due to oscillatory motion, while the second, the m-term, represents the change in the local steady field due to the motion of the body.
• Linearized free surface condition: For steady potential, if ܷ/ඥ݃‫ܮ‬ is small, then free surface boundary can be approximated as: The unsteady potential must satisfy the combined kinematic and dynamic free surface conditions: • Sea bed boundary condition: For infinite depth of water, the sea bed boundary condition is expressed as: • Radiation condition: The radiation condition states that disturbance generated by the body vanish when the distance r from the body tends to infinity: Incident wave potential: The potential for the undisturbed incident wave field at a point (X, Y and Z) in the fluid domain is known.For 2 nd order of Stokes wave as incident wave, then incident wave potential is: ሺ19ሻ

Diffraction and radiation wave potential:
The radiation or diffraction velocity potential located inside the fluid domain can be represented as: The source strength on the body surface is solved by satisfying the following: equation: For vehicle motion in infinite fluid, the surface integral involves only body boundary with: where, ‫ݎ‬ ொ represents the distance between the source point ܲ and field point Q.Thus, we get: By the application of panel method introduced by Hess and Smith (1967), the solution of three dimensional velocity potentials for the underwater vehicle with specified motion can be obtained: Wave forces and moment: Having solved the Laplace equation for Φ, the hydrodynamic pressure distribution ‫‬ on the body with forward speed ܷ can be obtained from the Bernoulli's equation as: From the pressure distribution, the various fluid forces may be calculated by integrating the pressure over the body surface.The expressions for force and moment on submerged body at any instant in time are as follows:

NUMERICAL APPROACH
The first step towards the numerical simulations in AQWA is defining the geometry of the test platform.This is done by creating full scale 3D AUV model having 1850 mm length and 200mm diameter in UniGraphics software.Figure 2 shows the AUV model which is equipped with four identical control fins located at stern in cruciform arrangement.
Once the 3D geometry of AUV is modeled, the next step is importing it to Ansys APDL where the surface mesh is generated using Shell 63 quadrilateral elements.The meshed model is then translated into AQWA using APDL macro.Figure 3 shows the meshed model of AUV with approximately 2000 panels.The function of macro is to map surface mesh elements to panels (or facets) and create an AQWA input file that contains information regarding distribution of panels, body-mass properties, wave conditions, etc.Finally, AQWA simulations are carried out in AQWA-LINE and AQWA-NAUT modules which output hydrodynamic response of the AUV in specified wave conditions.

FREQUENCY DOMAIN ANALYSES
Computations of linearized wave loading on AUV are accomplished by performing frequency domain analyses in AQWA-LINE module.Series of numerical simulations are carried out to investigate the effect of wave forces and moments on AUV at various depths and frequencies in head seas.In all simulations, deep sea conditions are considered with 25 different frequencies ranging from 0.05~5 rad/sec at 1~15 m submergence depths.particle's displacement decreases exponentially with the water depth.It can also be observed that the peak position of amplitude of the wave forces and moments are transferred from higher to lower frequencies as the depth of the vehicle increases.This can be explained by understanding the fact that the waves with higher frequency decrease faster than those of lower frequency with increased in the water depth.Further examining the results reveal that magnitudes of sway force, roll moment and yaw moment are relatively much smaller than surge force, heave force and pitch moment and hence will be ignored in subsequent simulations.

TIME DOMAIN SIMULATIONS
Time domain simulations are performed in the AQWA-NAUT program.The program requires a full hydrostatic and hydrodynamic description of AUV to simulate the response in regular or irregular waves.The description is transferred directly from the output results of the AQWA-LINE analyses done in previous section.Figure 10 and 11 present the response of the wave forces and moment when the AUV operating in a regular wave field at 4m depth with 5m/s velocity.The incident wave is modeled using 2 nd order Stokes theory with frequency 1 rad/s and amplitude 0.3 m.Using MATLAB curve fitting tool, the surge force, ‫ܨ‬ ௫ can be

SURGE FORCE FORMULATION
From the curve fitting function, the surge wave force, ‫ܨ‬ ௫ acting on AUV at any instant of time t can be written as: where, ‫ܣ‬ represents the amplitude of the force and ߜ is the initial phase.If the AUV is moving with constant forward velocity, ܷ, then the above relation becomes: where, ܿ = ߱/݇ is the wave velocity such that: In order to obtain the mathematical model for surge force, ‫ܨ‬ ௫ numerical simulations are performed in AQWA-NAUT to analyze the results by varying following influence factors.
Wave amplitude, : In Fig. 12 the wave force amplitude, ‫ܣ‬ is plotted against different values of wave amplitudes.The plot shows that force amplitude increases with an increase in the wave amplitude and their relationship is linear.Owing to this linear relation between ‫ܣ‬ and ܽ, the next simulations are performed with unit wave amplitude.

Depth of submergence, :
From the above relations, it is found that constants appearing in exponents are the values of wave numbers for corresponding frequencies i.e., ݇ = 0.0255,0.1019,1.723.

Wave frequency, ࣓:
To find the relation between the force amplitude and wave frequency, results are investigated for ߱ = 0.3~1.6ሺ‫ݏ/݀ܽݎ‬ሻwhile the vehicle is at depth, ℎ = 4݉.Figure 14 shows that an increase in wave force amplitude, ‫/ܣ‬ሺܽ݁ ି ሻ with increase in wave frequency, ߱.Using least square fitting, following relation can be obtained:   Figure 15 shows the influence on the force amplitude for various value of the pitch angle i.e., ሺߙ = ±12ሻ with same test conditions as in the side slip angle.
It can be observed that with increment in the pitch angle the wave force curve follows a parabolic trend.Due to small ‫ܮ‬ ߣ ⁄ ratio, the difference in the results for positive and negative pitch angle is insignificant.On curve fitting following equation is obtained.

Initial phase, :
A careful observation of foregoing simulation results reveals that the initial phase δ is only the function of pitch angle, α.Therefore, influence on phase, δ is investigated by varying vehicle's pitch angle, α, while keeping other parameters constant as presented in the Fig. 16.It is found that δ varies linearly with α, such that: Consolidating the relations, the analytical formulation obtained for surge wave force is: Among which ߈ = ݁ ି ߱ ଵ.ଽଶ ‫ߚݏܿ‬ and ‫ܥ‬ ଵ ሺߙሻ = 89.22ߙଶ + 3.291ߙ + 56.48.

HEAVE FORCE AND PITCH MOMENT FORMULATION
As mentioned earlier that the heave force, ‫ܨ‬ ௭ and pitch moment, ‫ܯ‬ ௬ can be represented by 1st order Fourier function as: In the above equation, the unknowns ܽ , ܽ ଵ and b ଵ are determined by adopting the same procedure as explained for the surge force formulation.
In case of the heave force, ‫ܨ‬ ௭ , the unknowns obtained are following: a = 5.529a ଶ.ଵସ ω ସ.ଵ e ିଶ୩୦ + 1.55 Similarly for the pitch moment, ‫ܯ‬ ௬ the unknowns are:

VERIFICATION TESTS
In order to verify the accuracy of analytical formulations, two test conditions are considered by taking different values of the parameters as shown in Table 2.

CONCLUSIONS
The study is focused on a comprehensive numerical study for computing the wave forces and moments acting on an AUV during its operation close to the sea surface.The mathematical formulation regarding the wave-body interaction is presented, which is based on Panel Method.The capabilities of Panel Method based code, Ansys AQWA, are utilized to perform the numerical study.The code consists of several modules, which can perform a variety of simulations both in the frequency domain and time domain thus providing detailed results necessary to assess the effect of wave loading on the underwater vehicle.Several wave force calculations are carried out taking into account the 2nd order Stokes waves and heads seas condition.Firstly, the numerical analyses are performed in the frequency domain with intention to investigate the effect of varying the submergence depths and wave frequencies on the magnitudes of wave forces and moments.The results indicate that as the vehicle's depth increases the effect of wave loading decreases until it becomes almost insignificant at 15m depth.Moreover, the magnitudes of sway force, roll moment and yaw moment are much smaller when compared with surge force, heave force and pitch moment.Numerical simulations in the time domain are performed to examine the individual influence of each parameter on wave loads.This is done by varying one parameter and keeping the other constant.The parameters examined are the wave amplitude, frequency, vehicle depth, orientation and initial phase.Based on the results gathered from the simulations and applying curve fitting functions, analytical formulations representing surge force, heave fore and pitch moment are derived.The accuracy of the formulations is judged through comparison with the real time simulations.This is accomplished by implementing two test conditions in which the different values of parameters are selected.The results show very good agreement with mean square error not more than 2%.Therefore, it can be concluded that, the proposed method provides an efficient means of calculating hydrodynamic forces due to wave and can be used to evaluate the sea keeping and the maneuvering characteristics of underwater vehicles during initial design stages.

Fig. 1 :
Fig. 1: Sketch of the problem geometry where j = 1, 2, 3, 4, 5 and 6 refer to surge, sway, heave, roll, pitch and yaw motions respectively.Governing equations: Assuming that the fluid is absolutely homogeneous and incompressible, the fluid density ߩ remains constant throughout the fluid domain and the continuity equation reduces to: ப୳ ப୶ + ப୴ ப୷ + ப୵ ப = 0(1) is the location of source on the surface S, ‫,ݔ‪ሺ‬ܩ‬ ‫,ݕ‬ ‫,ݖ‬ ܲሻ is the Green function which describes the flow at (x ,y and z) location caused by source ߪሺܲሻ of unit strength at ܲ.The Green function satisfies the Laplace equation and the boundary conditions everywhere except normal velocity ‫ݒ‬ boundary condition on body surface:

Fig. 12 :
Fig. 12: Force amplitude versus wave amplitude at various depths represented by a sine function while heave force, ‫ܨ‬ ௭ and pitch moment, ‫ܯ‬ ௬ by 1st order Fourier function.The unknown coefficients (variables) involved in these functions can be obtained by observing the results of various factors influencing the wave forces and moment.The next section describes the method of deriving the analytical expression of wave forces and moment.
Figure13shows the force amplitude, ‫ܣ‬ at submergence depth, ℎ ranging from 1~15 m for three different ω values.

Fig. 13 :
Fig. 13: Variation of force amplitude with depth at various frequencies

Fig. 17 :
Fig. 17: Comparison of analytical versus AQWA results for surge force

Table 1 :
Force amplitude at various side slip angles Table1shows the magnitude of wave force amplitude at various side slip angles, ߚ when operating depth, h is 2 m.Due to the symmetry of the geometry, the effect of varying only the positive values of β are investigated.Using curve fitting the relationship obtained is:

Table 2 :
Parameters used in the verification tests

Table 3 :
Comparison of analytical and AWQA results SSE It is clear from the figures that analytical formulations accurately predict the wave forces and moment.The accuracy of the analytical formulations can also be observed from Table3, which gives the comparison in terms of the Sum of Square Error (SSE) for the two test conditions.The results show the mean square error for surge force and pitch moment is less than 1% while in case of heave force it is not more than 2%.