A CFD-Based Surrogate for Duct-Jet AUV Maneuvering

1. Introduction

Autonomous Underwater Vehicles (AUVs) have emerged as one of the most actively studied technologies in marine engineering, driven by their versatility across scientific, industrial, and defense applications. Unlike traditional manned or remotely operated systems, AUVs offer persistent and cost-effective operation in hazardous and remote environments, reducing risks to human operators while expanding the scope of underwater exploration. Their capability to operate independently makes them particularly attractive for missions such as seabed mapping, environmental monitoring, infrastructure inspection, and naval reconnaissance [1,2,3].

Recent years have witnessed rapid advancements in the enabling technologies that underpin AUV design and operation. Improvements in energy storage, including high-density lithium-ion batteries and fuel-cell systems, have significantly extended mission endurance. Progress in propulsion system design, aided by computational fluid dynamics (CFD) and reduced-order modeling, has enhanced maneuverability and efficiency in complex hydrodynamic environments. At the same time, miniaturization of high-resolution sensors, such as sonar, optical imaging devices, and chemical probes, has broadened the range of measurable ocean parameters, enabling high-fidelity mapping of physical, chemical, and biological processes [4,5,6].

A particularly transformative development lies in autonomy and onboard intelligence. Advances in machine learning, adaptive control, and real-time data processing now allow AUVs to make mission-critical decisions under uncertainty, adapt their trajectories to environmental variability, and coordinate as swarms for distributed sensing. These capabilities have expanded the operational horizon of AUVs from coastal monitoring to deep-ocean exploration and under-ice missions, where human access is limited or impossible [7,8,9,10,11,12].

Beyond technological progress, global challenges have further motivated AUV research. Climate change, biodiversity conservation, and pollution mitigation require sustained, large-scale ocean observation with spatial and temporal resolutions unachievable by conventional ship-based surveys. Similarly, offshore industries and defense operations demand reliable, long-duration platforms for inspection, surveillance, and asset protection. As a result, AUVs are now regarded as indispensable tools for addressing both fundamental scientific questions and pressing societal needs [13,14,15].

In this context, the development of robust motion prediction models and surrogate frameworks is critical for ensuring the safe, efficient, and autonomous operation of AUVs. High-fidelity CFD simulations provide detailed insights into the complex flow physics of propulsion and maneuvering, however, their computational cost limits real-time applicability. Reduced-order models, surrogate approaches, and hybrid strategies have therefore become central research directions, bridging the gap between accuracy and computational efficiency. The present study contributes to this ongoing effort by developing and validating a surrogate motion prediction model for a duct–jet propelled AUV, with a focus on balancing fidelity and practicality for free-running simulations [16,17,18].

The AUV considered in this study is equipped with a ducted propeller that enhances propulsive efficiency and suppresses cavitation, making it well suited for low-speed, high-thrust operations. Its principal particulars are presented in Section 2. The CFD framework, including the computational setup, test matrix, and governing equations, is described in Section 3. Section 4 presents the methodology for constructing the surrogate model. Results of the CFD simulations are presented in Section 5, while Section 6 discusses the validation of the surrogate model. Finally, conclusions are summarized in Section 7.

2. AUV

2.1. Main Particulars

Table 1 presents the principal specifications of the 2 m-class AUV which is investigated in this study. The mass and center of gravity (COG) of the AUV are obtained from CFD hydrostatic simulations without the use of design variables. The COG is defined as the distance measured from the foremost point of the AUV. The radii of gyration and moments of inertia (MOI) are estimated based on design variables, with the COG as the reference point. The overall configuration of the AUV is shown in Figure 1(a). The holes distributed along the hull, as shown in Figure 1, are designed as mounting slots for vertical and horizontal thrusters. In this study, however, the thrusters themselves are not modeled in either the CFD simulations or the surrogate model, and only the holes remain in the geometry.

The rudders are designed such that the spans of the port and starboard rudders are larger than those of the vertical rudders, in order to enhance controllability in the horizontal plane (Figure 1(b)). Since this study focuses on horizontal-plane maneuvering, only the vertical rudder span is listed in the table. It is also noted that the rudders have relatively small spans compared to the AUV length for its capability to be launched from a tube.

The propeller duct is supported by four upstream and thirteen downstream struts, arranged in an axisymmetric manner (Figure 1(c)). The propeller consists of five axisymmetric blades.

2.2. Coordinate System

This study employs two coordinate systems with parallel axes but different origins in both the CFD and surrogate formulations: the Earth-fixed frame and the body-fixed frame. The origin of the Earth-fixed frame is located at the geometric center of the AUV’s nose, while the origin of the body-fixed frame coincides with the center of gravity. Since the present study focuses on 3 degrees of freedom (DOF) motion, the Earth-fixed frame is used to compute translational displacements in the body-fixed frame—namely surge (x), sway (y), and yaw (ψ)—and to measure the longitudinal offset of the COG from the AUV’s nose. The axes are defined such that the +x axis points toward the tail, the +y axis points to starboard, and the +z axis points upward.

3. CFD

3.1. URANS

In this study, the incompressible URANS equations are numerically solved using unstructured grid–based finite volume method (FVM). The commercial CFD solver STAR-CCM+ V2302 (Siemens Digital Industries Software, Plano, TX, USA) is employed. Equation (1) presents the governing equations of incompressible URANS:

$${\displaystyle \frac{\partial \overline{\mathit{v}}}{\partial t}}+\nabla \cdot \left[\left(\overline{\mathit{v}}-{\mathit{v}}_g\right)\otimes \overline{\mathit{v}}\right]=-\nabla \overline{p}/\rho +\left(\nu +{\nu }_t\right){\nabla }^2\overline{\mathit{v}}$$

(1)

where $\mathit{v}$, $p$, $\nu$, and ${\nu }_t$ denotes the flow velocity, piezometric pressure, kinematic viscosity, and eddy viscosity, respectively, and upper bar indicates Reynolds-averaged quantities. When the grid at certain coordinate system moves with both translational and rotational speed, with respect to other coordinate system, the grid velocity (${\mathit{v}}_g$) is defined as Equation (2):

$${\mathit{v}}_g={\mathit{v}}_{g,t}+{\mathit{\omega }}_g\times \mathit{r}$$

(2)

where ${\mathit{v}}_{g,t}$ is the translational velocity vector of the grid, ${\mathit{\omega }}_g$ is the rotational velocity vector of the grid, and $\mathit{r}$ is the position vector of a mesh vertex.

The current simulations employ Menter’s SST $k-\omega$ model [19], a commonly used anisotropic two-equation turbulence model based on the eddy-viscosity formulation for turbulence closure.

3.2. Numerical Schemes

For temporal discretization, the second-order implicit backward differentiation that incorporates five time levels is used. Convective fluxes of momentum equations and transport equations for turbulence variables are computed with the second-order upwind scheme. Diffusive fluxes are dealt with the central differencing scheme formulated with non-orthogonality adjustment terms, namely the cross-diffusion terms. Gradients are computed with the Gauss linear method with Venkatakrishnan limiter. For pressure-velocity coupling, a semi-implicit method for pressure linked equations (SIMPLE) algorithm is used. The implemented blended wall-functions are used as near-wall treatments for velocities, $k$, and ${\omega }_t$.

3.3. Equations of Motion

In predicting the motion of a moving body, the key quantities obtained from the fluid field are the forces and moments acting on the body. The shear stresses and pressure forces derived from the discretization of Equation (1) are ultimately computed in the Earth-fixed coordinate system. On the other hand, when the body accelerates in a non-inertial frame, the governing equations of motion can be formulated based on Newton’s second law with respect to the body-fitted coordinate system, as shown in Equation (3):

$$m\left[\dot{u}-vr+wq-x_{rg}\left(q^2+r^2\right)+y_{rg}\left(pq-\dot{r}\right)+z_{rg}\left(pr+\dot{q}\right)\right]=X$$

(3a)

$$m\left[\dot{v}-wp+ur-y_{rg}\left(r^2+p^2\right)+z_{rg}\left(qr-\dot{p}\right)+x_{rg}\left(qp+\dot{r}\right)\right]=Y$$

(3b)

$$m\left[\dot{w}-uq+vp-z_{rg}\left(p^2+q^2\right)+x_{rg}\left(rp-\dot{q}\right)+y_{rg}\left(rq+\dot{p}\right)\right]=Z$$

(3c)

$$I_x\dot{p}+\left(I_z-I_y\right)qr+m\left[y_{rg}\left(\dot{w}-uq+vp\right)-z_{rg}\left(\dot{v}-wp+ur\right)\right]=K$$

(3d)

$$I_y\dot{q}+\left(I_x-I_z\right)rp+m\left[z_{rg}\left(\dot{u}-vr+wq\right)-x_{rg}\left(\dot{w}-uq+vp\right)\right]=M$$

(3e)

$$I_z\dot{r}+\left(I_y-I_x\right)pq+m\left[x_{rg}\left(\dot{v}-wp+ur\right)-y_{rg}\left(\dot{u}-vr+wq\right)\right]=N$$

(3f)

where $u, v,w,p,q,r$ and $\dot{u},\dot{v},\dot{w},\dot{p},\dot{q},\dot{r}$ are 6DoF velocity and acceleration components, $X,Y,Z,K,M,N$ are 6DoF forces and moments components, $I_x, I_y,I_z$ are principal components of moment of inertia tensor of a moving body defined in the body-fixed coordinate. The reference point of moments of inertia is the center of rotation, which is assumed the same as center of gravity in this study. $x_{rg},y_{rg},z_{rg}$ are the components of distance vector pointing center of gravity from center of rotation, thus are all zero values in the current study. Due to the relative reference between the inertial and non-inertial frames, apparent inertial effects such as the Coriolis forces ($m\left(p,q,r\right)\times \left(u,v,w\right))$ arise in the non-inertial coordinate system.

Through Euler angles, kinematic variables obtained in each coordinate system can be transformed between the inertial and non-inertial frames. However, since the moments of inertia remain constant in the body-fixed frame, and in the surrogate model all quantities, except for trajectory prediction, are considered in the body-fixed frame, the formulation of the equations of motion in the body-fixed coordinate system is more straightforward in both conceptual and computation aspect.

3.4. Verification Method

This study verifies grid resolution for both propeller characteristics and hull resistance. Instead of evaluating uncertainties through Richardson extrapolation, a more practical approach is adopted by constructing three grid systems with different resolutions and examining the relative differences in the solutions with respect to grid convergence. Richardson extrapolation often exhibits somewhat oscillatory convergence behavior. Such oscillations are considered not to originate from the overset methodology itself, but rather from the fact that the employed grid sizes are already close to the asymptotic solution.

3.5. Grid System

Since the AUV is a submerged vehicle without a free surface, the background grid is configured to move together with the body. Local mesh refinement is applied in regions where strong pressure gradients are expected, such as the nose and stern with rudders, where the cell size is reduced to one-half or one-quarter of the cell size located at mid-body. Around the propeller, an even finer refinement of up to one-eighth of the mid-body size is employed (Figure 2). In the mid-body region, the longitudinal grid size is approximately 0.0141 m in the G2 mesh system. The background grid is designed to gradually expand away from the mid-body, reaching up to eight times the mid-body cell size. For the boundary layer, a near-wall spacing of approximately $6.0\times {10}^{-4} m$ is applied around the hull including appendages, ensuring $y^{+}<30$in the G2 case.

The four rudders, each moving with 1 DOF, and the propeller are allowed to move relative to the hull. To achieve this, grid separation between the hull and appendage meshes is required. An overset grid approach is used such that the rudders and propeller are embedded within the hull and duct grid systems, while the duct is unified with the hull via struts. Boundary values between overset regions are interpolated linearly (Figure 3). Overset interfaces are also applied between rudder regions due to their mutual interference.

As aforementioned, the entire grid system is prepared at three different resolutions. Starting from the baseline G1 mesh, coarser meshes G2 and G3 are generated by increasing the characteristic grid size by a factor of $\sqrt{2}$ successively. Two grid-triplet sets are constructed: one for the propeller open-water (POW) test and another for the remaining simulations including hull and rudders. The corresponding numbers of grid points are listed in Table 2 and Table 3, respectively

3.6. Domain, Boundary Conditions and Initial Conditions

Figure 4 shows the computational domain, and Table 4 summarizes the boundary conditions. For the POW test, only the duct and propeller are included, as illustrated in Figure 2(d), and the domain size is scaled by $L/\left(10D_P\right)$. In general, background grid boundaries impose streamwise ($U$) and transverse ($V$) velocities, and the flow field inside the domain is initialized with the same velocities. In simulations where the body undergoes prescribed displacements, namely forced oscillation, pure sway, pure yaw, and zigzag tests, zero inflow conditions ($U=V=0$) are applied at the boundaries, and the initial condition is also set to zero. In these cases, the forces acting on the hull are computed in terms of the relative velocity obtained from the equations of motion.

3.7. Body-Force Propulsion

In this study, the POW simulations are first conducted using a discretized propeller (DP), and the obtained POW characteristics, i.e., thrust coefficient ($K_T$) and torque coefficient ($K_Q$) are subsequently employed to establish a body-force (BF) propulsion model for the free-running simulations. Since the propeller located inside the duct induces a local acceleration of the flow near the duct inlet, the averaged axial propeller inflow speeds at the duct inlet ($u_A$) are measured at a cross-section (Figure 5) in the discretized propeller simulations and utilized for formulating the advance coefficient at the duct inlet ($J_A=u_A/\left(nd\right)$) rather than utilizing the advance coefficient based on ship speed ($J=u/\left(nd\right))$. Based on these measurements, an additional correlation (Equation 5) between the $J$ and $J_A$ is derived, ensuring that the thrust and torque coefficients could be properly adjusted to act with respect to the hull-relative velocity.

$$J_A=J_{1}J+J_0$$

(5a)

$$K_T=K_{T0}+K_{T1}{J}_A+K_{T2}{J}_A^2$$

(5b)

$$K_Q=K_{Q0}+K_{Q1}{J}_A+K_{Q2}{J}_A^2$$

(5c)

3.8. Controller

In the free-running simulations, a controller is required to compute and adjust the vertical rudder angle ($\delta$) at each time step based on the yaw motion ($\psi$) of the AUV. In this study, the yaw angle reported at every time step is linked to a JavaScript-based macro, which calculates the corresponding rudder angle. The computed rudder angle is then fed back into STAR-CCM+ through a field function to rotate the rudder mesh accordingly.

4. Surrogate

The purpose of the surrogate model is to retain the accuracy of CFD while significantly reducing the computational cost. To achieve this, it is advantageous for the surrogate to adopt modules as similar as possible to those used in CFD. The surrogate is therefore designed to employ the same coordinate system and input variables, as well as the same equation of motion and rudder controller. The only distinction lies in the estimation of forces and moments.

While the CFD solver requires only the inputs listed in Table 1 and the discretized geometry for meshing, the surrogate further modularizes the force and moment estimation. Accordingly, in addition to the values in Table 1, additional inputs specific to force and moment estimation are required. These inputs are detailed in the following section, where the governing expressions for forces and moments are presented.

4.1. Algorithm Flowchart

Figure 6 presents the algorithmic flowchart of the surrogate model used in this study. The surrogate performs all calculations at certain time step using the AUV velocity variables updated from the accelerations of the previous time step through first-order Euler time integration; at the first time step, the initialized velocity values are used. In the force and moment calculation, the propeller loads are computed prior to those of the hull and rudders, since a portion of the propeller thrust contributes to the hull through the thrust-deduction effect.

Although the general form of the equations of motion is identical to that in CFD, the present surrogate employs a nonlinear iterative scheme, in which the surge, sway, and yaw motions are converged sequentially. The convergence tolerance is defined such that the average residual falls below 10⁻¹⁰.

4.2. Effective Flow Speed

The effective flow speed for scaling maneuver coefficients and kinematic properties is defined as the follows:

$$U_e=\sqrt{u^2+v^2}$$

(6)

4.3. Hull Incident Angle

The incident angle of the hull for estimating incident angle of rudder and propeller is defined as follows:

$$\beta ={\tan }^{-1}\left(-v/u\right)$$

(7)

4.4. Forces and Moments – Propulsion

A propulsion model based on the POW curves is established by adopting the propeller force and moment module from the standard MMG formulation [20]. In addition, the correlation between the AUV-relative velocity and the duct inflow velocity obtained from the CFD body-force propulsion approach (Equation 5) is incorporated once the advance coefficient is obtained via Equation (8c).

$${\beta }_P=\beta +\left(x_P/L\right)\left(\psi L/U_e\right)$$

(8a)

$${\omega }_P={\omega }_0-exp\left(-25{\beta }_P^2\right)$$

(8b)

$$J=u\left(1-{\omega }_P\right)/\left(n{D}_P\right)$$

(8c)

where ${\beta }_P$, ${\omega }_P$, ${\omega }_0$, $n$ are incident angle towards the propeller, wake fraction of the propeller inflow in the oblique flow, wake fraction of the propeller inflow in the straight flow, and propeller rotational speed, respectively. The thrust coefficient obtained by Equation (5) with Equation (8c) is used to estimate dimensional thrust of the propeller as follows:

$$X_P=\rho n^2{D}_P^4{K}_T$$

(9)

4.5. Forces and Moments – Hull and Rudder

The hull and rudder force and moment module employs a whole-ship approach [21] and the angle of attack of the rudder during the exposure of the hull in the oblique flow is estimated as follows.

$${\beta }_R=\beta +\left(x_R/L\right)\left(\psi L/U_e\right)$$

(10a)

$$\alpha =\delta -\tan \left({\beta }_R\right)$$

(10b)

The resistance is formulated to be the function of surge speed as follows.

$$R=R_0+R_{1}u+R_2{u}^2$$

(11)

The final form of forces and moment of hull and rudder comprises the effect of sway velocity ($v)$, sway acceleration ($\dot{v}$), yaw rate ($r$), yaw acceleration ($\dot{r}$) and $\delta$ in a Taylor-expansion form.

$$X_H=F\left(-R+X_{vv}{v^{\prime }}^2+X_{vvvv}{v^{\prime }}^4+X_{vvvvvv}{v^{\prime }}^6+X_{\delta }\left|\alpha \right|+X_{\delta \delta }{\alpha }^2\right)+t_p{X}_P$$

(12a)

$$Y_H=F\left(Y_v{v}^{\prime }+Y_{v\left|v\right|}{v}^{\prime }\left|v^{\prime }\right|+Y_{vvv}{v^{\prime }}^3+Y_{\delta }\alpha +Y_{\delta \left|\delta \right|}\alpha \left|\alpha \right|+Y_{\delta \delta \delta }{\alpha }^3+Y_{\dot{v}}{\dot{v}}^{\prime }+Y_r{r}^{\prime }+Y_{r\left|r\right|}{r}^{\prime }\left|r^{\prime }\right|+Y_{\dot{r}}{\dot{r}}^{\prime }\right)$$

(12b)

$$N_H=FL\left(N_v{v}^{\prime }+N_{v\left|v\right|}{v}^{\prime }\left|v^{\prime }\right|+N_{vvv}{v^{\prime }}^3+N_{\delta }\alpha +N_{\delta \left|\delta \right|}\alpha \left|\alpha \right|+N_{\delta \delta \delta }{\alpha }^3+N_{\dot{v}}{\dot{v}}^{\prime }+N_r{r}^{\prime }+N_{r\left|r\right|}{r}^{\prime }\left|r^{\prime }\right|+N_{\dot{r}}{\dot{r}}^{\prime }\right)$$

(12c)

where the kinematic variables are expressed in non-dimensional form ($v^{\prime }=v/U_e$, $r^{\prime }=rL/U_e, {\dot{v}}^{\prime }=\dot{v}L/U_e^2, {\dot{r}}^{\prime }=\dot{r}{L}^2/U_e^2)$ with the non-dimensional representation of forces and moments, which are scaled by $F=\left(0.5\right)\rho LT{U}_e^2$ and $FL$, respectively. $t_p$ denotes thrust deduction factor.

4.6. Equations of Motion

The 3DOF equations of motion, including the added masses ($m_x$, $m_y$) and the added moment of inertia ($J_z$​), can be simplified as follows when the rotational center coincides with the center of gravity by only considering the diagonal added mass matrix [22]. Once the external forces and moment are obtained, the kinematic variables are sequentially updated through nonlinear iteration until convergence is achieved, yielding the final accelerations.

$$\left(m+m_x\right)\dot{u}-\left(m+m_y\right)vr=X_H+X_P$$

(13a)

$$\left(m+m_y\right)\dot{v}+\left(m+m_x\right)ur=Y_H$$

(13b)

$$\left(I_z+J_z\right)\dot{r}=N_H$$

(13c)

5. CFD Results

Table 5 summarizes the types of CFD simulations conducted for constructing the surrogate model, and includes the final zigzag simulation for the validation after the surrogate construction. For evaluating the accuracy of hull and propeller force and moment simulations, the grid triplet is employed, while the medium grid is used for the other simulations. In the POW test with the discretized propeller, a time step of $2.5\times {10}^{-4} s$ is applied, whereas all other simulations, including those with the surrogate model, use a time step of $0.01 s$. The basic composition of the test matrix follows the procedure described in [23].

5.1. Hydrostatic

A hydrostatic simulation is carried out with the AUV fixed and the relative inflow velocity set to zero in order to determine the mass of the fully-attached AUV and the longitudinal center of gravity ($x_G$). Although mass and $x_G$ are often provided from experiments, in this study they are derived directly from the simulated geometry to ensure an internal balance between buoyancy and weight as well as a zero trim moment. Since most simulations are performed under captive conditions and the free-running cases are limited to 3DOF motion, the $x_G$ does not affect the vehicle’s attitude. However, as it serves as the effective center of rotation and the reference for moments, the $x_G$ obtained from the hydrostatic analysis is consistently employed in the subsequent simulations to improve the prediction of yaw moments. The $x_G$ is determined iteratively within the hydrostatic simulation by adjusting its position until the trim moment converges to zero. The final values of mass and $x_G$ are presented in Table 1.

5.2. Propeller Open-Water Test

As noted earlier, the hydrodynamic properties of the propeller are obtained using a discretized propeller under straight-ahead conditions (Figure 7). Figure 8 and Table 6, Table 7 and Table 8 present the thrust, torque, and propulsive efficiency computed from the grid triplet. The results show that the solutions converge toward nearly asymptotic values as the grid resolution increases. Among the three grids, the G2 case, which provides a median value for computational cost and accuracy, is considered suitable for subsequent body-force modeling.

Since the duct-jet propulsive system includes the duct itself, the convergence of the duct resistance is also examined (Figure 8(d)). It is confirmed that the resistance on the duct achieves an accuracy comparable to that of the G1 when computed with the G2.

As shown in Table 6, Table 7 and Table 8, the differences between the solution values ($s_1$, $s_2$, $s_3$) obtained from G1, G2, and G3 exhibit oscillatory convergence, which is reflected in the ratios $R={\epsilon }_{12}/{\epsilon }_{23}$​. Consequently, Richardson extrapolation does not provide reliable uncertainty estimates in many cases.

In the duct-jet propulsive system, deriving the regression formula for POW between the conventional advance coefficient $J$ and the propeller hydrodynamic properties without additionally measuring the flow velocity at the duct inlet leads to difficulties in body-force implementation. This is because the propeller inflow plane, which is typically prescribed at the propeller inlet, exhibits a flow velocity that increases significantly compared to the hull-relative velocity when a duct is present.

To address this, an additional regression relation is introduced between the hull-relative advance coefficient $J$ and the duct-inlet advance coefficient $J_A$​. Accordingly, the regression formula for the propeller hydrodynamic properties is expressed in terms of $J_A$​ rather than $J$. The computation of $J_A$​ is performed at the inflow plane shown previously in Figure 5, and the relation between $J$ and $J_A$​ is presented in Figure 9 and Table 9. In Equation (5a), the values of $J_1$​ and $J_0$​ are set to 0.5123 and 0.5096, respectively. The resulting variations in the POW coefficients are summarized in Table 10 and Table 11. The effect of grid resolution on the computed $J_A$​ values is found to be negligible.

The body-force modeling is carried out using the POW regression formula based on $J_A$​, and validation against the discretized propeller is performed through POW simulations. Figure 10 presents a comparison of the hydrodynamic properties. For thrust, the body-force model shows an overestimation in the region of $J<0.5$, although the overall accuracy remains reasonable. In the case of torque, the body-force model generally exhibits an overestimation. This behavior appears unrelated to the accuracy of $J_A$​ measured for the discretized and body-force propellers, as shown in Figure 11(d). Instead, it is considered to result from the interaction between the duct and the modified pressure field induced by the body-force propeller.

Since this issue is also related to grid resolution, the present study proceeds without further debugging to maintain a practical approach. Ultimately, torque does not influence the 3DOF simulations, and therefore no additional corrective measures are taken. Body-force modeling is not only important for predicting thrust but also for reproducing the flow field around the hull with sufficient fidelity, particularly because the rudder and propulsion system are usually placed in close proximity. As confirmed at representative advance coefficients, the body-force model predicts a velocity distribution similar to that of the discretized propeller (Figure 12).

5.3. Resistance Test

The resistance test (Figure 13) focuses primarily on constructing Equation (11) and performing CFD data verification for AUV hull and rudders using the grid triplet.

In this setup, the duct is included, while the propeller is excluded to enable replacement with the body-force model. The simulation results are presented in Figure 14. In addition, the measured frictional resistance during straight-ahead motion of the AUV shows close agreement with the ITTC 1957 regression formula (frictional resistance coefficient, $C_F=0.075/{\left(\log \left(Rn\right)-2\right)}^2)$ commonly used for ships. This suggests that if pressure-field regression is carried out for an AUV digital twin, the frictional resistance can be modeled as a major feature through a simple regression expression.

The coefficients applied in Equation (11) are listed in Table 12. Furthermore, the grid dependency of the axial flow speed at the duct inlet, expressed as $J_A$​, is examined (Figure 15). Qualitatively, it is confirmed that the G2 resolution yields results nearly identical to those of the G1 grid.

5.4. Self-Propulsion Test

In the self-propulsion test (Figure 16), the body-force propulsion is implemented and the simulation is carried out in 1DOF. At the target Froude number, the propeller rotational speed $n$, the value of ${\omega }_0$ used in Equation (8b), and the value of $t_P$​ used in Equation (12a) are prescribed. As shown in Figure 17(a), three cases using different $n$ are tested, and the final $n$ is determined at the point where the longitudinal force coefficient $X^{\prime }$ becomes zero. The propeller RPS obtained from CFD is directly applied in subsequent free-running simulations.

The values of ${\omega }_0$​ and $t_P$​, however, are determined differently. First, the surrogate expressions required for the self-propulsion test (Equations (5), (8), (9), (11), and (12a)) are constructed. Then, by executing the surrogate and adjusting ${\omega }_0$and $t_P$, values of $K_T$​, $X^{\prime }$, and ${X_P}^{\prime }$ comparable to those from CFD are identified. The final adjusted ${\omega }_0$ and $t_P$ values are 0.233 and 0.13, respectively.

5.5. Static-Drift and Control-Fin Tests

In the static-drift test (Figure 18(a)), the maneuvering coefficients associated with lateral velocity $v$ are identified, while in the control-fin test (Figure 18(b)), the maneuvering coefficients associated with the vertical rudder angle $\delta$ are determined using the least-squares method. The reconstructed forces and moments based on the fitted coefficients are presented in Figure 19. The coefficients applied in Surrogate Equation (12) are summarized in Table 14 and Table 15.

Table 13. Coefficients of regression curves achieved from static-drift simulations.

Coefficient	Value	Coefficient	Value	Coefficient	Value
$X_{vv}$	0.1574	$Y_v$	-0.1686	$N_v$	0.0331
$X_{vvvv}$	-1.9794	$Y_{v\left|v\right|}$	0.2658	$N_{v\left|v\right|}$	0.0616
$X_{vvvvvv}$	7.7387	$Y_{vvv}$	-2.2175	$N_{vvv}$	-0.4145

Table 13. Coefficients of regression curves achieved from static-drift simulations.

Coefficient	Value	Coefficient	Value	Coefficient	Value
$X_{vv}$	0.1574	$Y_v$	-0.1686	$N_v$	0.0331
$X_{vvvv}$	-1.9794	$Y_{v\left|v\right|}$	0.2658	$N_{v\left|v\right|}$	0.0616
$X_{vvvvvv}$	7.7387	$Y_{vvv}$	-2.2175	$N_{vvv}$	-0.4145

Table 14. Coefficients of regression curves achieved from control-fin simulations.

Coefficient	Value	Coefficient	Value	Coefficient	Value
$X_{\left|\delta \right|}$	-0.0018	$Y_{\delta }$	-0.0094	$N_{\delta }$	-0.0023
$X_{\delta \delta }$	0.0294	$Y_{\delta \left|\delta \right|}$	-0.0245	$N_{\delta \left|\delta \right|}$	-0.0093
		$Y_{\delta \delta \delta }$	-0.0082	$N_{\delta \delta \delta }$	-0.0069

Table 14. Coefficients of regression curves achieved from control-fin simulations.

Coefficient	Value	Coefficient	Value	Coefficient	Value
$X_{\left|\delta \right|}$	-0.0018	$Y_{\delta }$	-0.0094	$N_{\delta }$	-0.0023
$X_{\delta \delta }$	0.0294	$Y_{\delta \left|\delta \right|}$	-0.0245	$N_{\delta \left|\delta \right|}$	-0.0093
		$Y_{\delta \delta \delta }$	-0.0082	$N_{\delta \delta \delta }$	-0.0069

Table 15. Coefficients of regression curves achieved from pure-sway simulations.

Coefficient	Value	Coefficient	Value
$Y_{\dot{v}}$	-0.129	$N_{\dot{v}}$	-0.0088

Table 15. Coefficients of regression curves achieved from pure-sway simulations.

Coefficient	Value	Coefficient	Value
$Y_{\dot{v}}$	-0.129	$N_{\dot{v}}$	-0.0088

5.6. Pure-Sway and Pure-Yaw Tests

In the pure-sway test, the maneuvering coefficients associated with lateral acceleration $\dot{v}$ are obtained, while in the pure-yaw test, the coefficients associated with yaw rate $r$ and yaw acceleration $\dot{r}$ are identified. As introduced in the test matrix (Table 5), the pure-sway test is conducted by imposing a sinusoidal sway motion at the target Froude number, and the pure-yaw test is performed by superimposing a sinusoidal yaw motion (Figure 20(a) and (d)). The resulting forces and moments during pure-sway are presented in Figure 20(b) and (c), while those during pure-yaw are shown in Figure 20(e) and (f).

In the pure-sway test, the velocity-dependent forces and moments ($Y_{VEL}^{\prime }$, $N_{VEL}^{\prime }$​) are subtracted from the measured total force and moment ($Y^{\prime }$, $N^{\prime })$, respectively, using the surrogate regression expressions obtained from the static-drift test, and the acceleration-dependent terms ($Y_{ACC}^{\prime }$, $N_{ACC}^{\prime }$​) are identified through the least-squares method (Table 15). In the pure-yaw test, both the velocity-dependent and acceleration-dependent forces and moments are simultaneously determined using the least-squares method (Table 16).

5.7. Forced-Oscillation

In the forced-oscillation test, the AUV is subjected to a sinusoidal sway motion at $Fn=0$​ to determine the added mass coefficient $m_y$​ (Figure 21). To obtain the $m_y$, the simulation is briefly modeled assuming that the measured force is the sum of force affected by damping and inertia, which are assumed to be proportional to velocity and acceleration as follows.

$$m_y\ddot{y}+b\dot{y}=Y$$

(14)

where $b$ denotes damping coefficient. Since the acceleration term exhibits sine function behavior, the side force is convoluted with the same sine function sharing the same frequency to obtain the $m_y$. The coefficient $m_x$​ is found to be negligible in the results and is therefore not reported. Since added mass does not vary significantly, the period and amplitude of the sinusoidal input are selected to represent the typical operating range of the AUV. The non-dimensionalized added moment of inertia in z-direction 0.009 from [X] is directly used for the estimation of dimensional value ($J_z$​). A more accurate estimation of the AUV’s $J_z$ would require a forced-oscillation test in yaw.

6. Validation of Surrogate

Both CFD and the surrogate model are used to predict the results of the 3DOF zigzag maneuver, and the outcomes are compared to validate the surrogate. The validation is carried out with respect to the kinematic variables and $J_A$​ (Figure 22), as well as the forces and moments (Figure 23).

As shown in Figure 22, the surrogate continues to reproduce the overall motion period with reasonable accuracy even beyond the second overshoot. Although the velocity components are occasionally overpredicted, this does not alter the oscillation period. The surrogate accurately captures the decreasing trend of $J_A$​ immediately prior to the onset of turning. While $\beta$ tends to be overpredicted at both low and high values, the level of overestimation does not result in an unrealistic trajectory caused by excessive propulsion.

The comparison of forces and moments in Figure 23 further supports these observations. The propeller thrust follows the same trend as $J_A$​. The side force $Y^{\prime }$ and yaw moment $N^{\prime }$ rise sharply just before turning, with the yaw moment showing a suppressed response in CFD but a pronounced pulse in the surrogate. The surge force $X^{\prime }$ remains consistently smaller than one-tenth of $Y^{\prime }$, and the surrogate successfully reproduces the key trends of the turning behavior observed in CFD.

7. Conclusion

In this study, a surrogate model for AUV maneuvering is established using CFD, and its accuracy is assessed. The required CFD simulations are minimized, while the surrogate is analytically constructed through additional kinematic relations and a limited set of empirical parameters. The free-running predictions obtained from the surrogate demonstrate reasonable accuracy. It is anticipated that this framework can serve as a backbone for future developments, such as more advanced rudder models, propeller models incorporating side-force effects, and full 6DOF extensions.

In the free-running simulations, it is observed that the rudders of the current AUV design are relatively small, resulting in limited lifting force and confining zigzag motions to small values of drift angle $\beta$. This finding implies the need for further surrogate validation at larger $\beta$, particularly when side thruster models are incorporated. At the preliminary design stage, CFD proves to be valuable not only for predicting the flow field, but also for reproducing diverse scenarios through integration with controllers and propeller models. This reaffirms the advantage of predicting physical data directly from geometric information.

Nevertheless, the intrinsic limitation of high computational cost remains. Therefore, the present approach is expected to serve as a reduced-order model for hydrodynamic force prediction, particularly when coupled with control models or employed as training data for multivariate system identification.