2.1. Problem set-ups

We investigate the flow around a prolate spheroid pitching in a uniform upstream flow $U_{\infty }$ . Figure 1(a) sketches the geometric model of the prolate spheroid and its associated Cartesian coordinate systems. The earth-fixed system ( $X$ , $Y$ , $Z$ ) has its $X$ -axis aligned with the streamwise direction. The body-fixed system ( $\hat {X}$ , $\hat {Y}$ , $\hat {Z}$ ) is defined with the $\hat {X}$ -axis along the prolate spheroid’s major axis. The prolate spheroid pitches around the $\hat {Z}$ -axis, which is located at the centre of the prolate spheroid. The origin and $Z$ -axis of both coordinate systems coincide at the initial moment. In the body-fixed coordinate system, the coordinates ( $\hat {x}$ , $\hat {y}$ , $\hat {z}$ ) of the surface of the prolate spheroid (Miao & Xiao Reference Miao and Xiao2021; Hu et al. Reference Hu, Lin, Zhu, Yu, Lin and Li2024) are consistent with

(2.1) \begin{align} \frac {\hat {x}^{2} }{a^{2}} +\frac {\hat {y}^{2} }{b^{2}}+\frac {\hat {z}^{2} }{c^{2}}=1, \end{align}

where $a$ is the semi-major axis, and $b=c$ is the semi-minor axis. Here, the prolate spheroid’s minor-axis diameter is $D=2b$ , and its aspect ratio is defined as $\beta = a/b$ .

Figure 1.

Schematic of the computational set-up for the problem: (a) the two Cartesian coordinate systems; (b) the flow configuration in a three-dimensional rectangular box.

Figure 1. Long description

Panel A: A diagram showing two Cartesian coordinate systems. The first coordinate system is labeled with axes X, Y, and Z, and the second coordinate system is labeled with axes X-hat, Y-hat, and Z-hat. The prolate spheroid is positioned at the origin of both coordinate systems, with its major axis aligned along the X-axis. The angle phi represents the pitch angle between the two coordinate systems. The dimensions of the prolate spheroid are given as a equals 3D and b equals 0.5D, where D is a characteristic length. The free stream velocity U-infinity is directed along the X-axis. Panel B: A three-dimensional rectangular box representing the flow configuration. The box has dimensions 24D by 11D by 40D. The prolate spheroid is placed within the box, and its position and orientation are defined by the overset region. The flow enters the box at the inlet and exits at the outlet. The wake region downstream of the prolate spheroid is also depicted. The axes X, Y, and Z are labeled, and the pitch angle is indicated.

The kinematics of the prolate spheroid are prescribed by changing the inclination angle

(2.2) \begin{align} \phi \left ( t \right ) =\phi _{0} +\phi _{m}\sin \left ( 2\pi f t \right ), \end{align}

where $\phi _{0}$ is the average intersection angle of the $X$ -axis and $\hat {X}$ -axis, and $\phi _{m}$ is the pitching amplitude. The non-dimensional pitching frequency is given by $f=f_{*}D/U_{\infty }$ , where $f_{*}$ denotes the dimensional pitching frequency. Similarly, non-dimensional time is defined as $t=t_{*}U_{\infty }/D$ , with $t_{*}$ representing dimensional time. A complete non-dimensional pitching period is given by $T=1/f$ .

2.2. Numerical method

The unsteady flow around the pitching prolate spheroid is governed by the Navier–Stokes equations for incompressible flows:

(2.3a) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla }\boldsymbol{u} & = - \boldsymbol{\nabla }p+\frac {1}{Re_{D}} \nabla ^{2} \boldsymbol{u}, \\[-12pt] \nonumber \end{align}

(2.3b) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} & = 0, \\[10pt] \nonumber \end{align}

where $\boldsymbol{u} = ( u_{x} ,u_{y} ,u_{z})$ denotes the velocity vector in the Cartesian coordinates, $p$ denotes the pressure, and $ \textit{Re}_{D}=U_{\infty }D/\nu$ denotes the Reynolds number, with $\nu$ the kinematic viscosity. The uniform free-stream velocity ( $U_{\infty }$ ) is arranged at the inlet, and the fixed pressure condition is stipulated at the outlet. The slip-wall condition is operated at the other boundaries, and the no-slip condition is stated on the prolate spheroid.

The governing equations are discretised on an unstructured mesh using the finite volume method (Manik, Dalal & Natarajan Reference Manik, Dalal and Natarajan2018). The time derivative is discretised with a second-order implicit backward scheme with an adaptive time-stepping protocol employing a dynamic global Courant number constraint ( $\text{CFL} \leqslant 0.8$ ). Spatial discretisation adopts the second-order central differencing scheme. For solving the discretised equations, the preconditioned bi-conjugate gradient method with a diagonal-based incomplete LU preconditioner is used for the pressure fields. The reader can refer to Mokalled, Mangani & Darwish (Reference Mokalled, Mangani and Darwish2016) for more detailed descriptions of the numerical implementation.

The three-dimensional view of the pitching prolate spheroid is shown in figure 1(b). The present simulation is performed in a rectangular domain $[ -8 D, 32 D ] \times [ -8 D, 16 D ] \times [ -5.5 D, 5.5 D ]$ . This study concerns a prolate spheroid with a constant aspect ratio $\beta = 6$ and baseline inclination $\phi _{0} = 45^{\circ }$ . Consequently, the computational domain dimensions are identical for all cases; their verification is provided in Appendix A. In this study, we consider a range of Reynolds numbers ( $ \textit{Re}_{D}$ ) $50{-}300$ . Notably, the diameter of the volume-equivalent sphere ( $d =2 \sqrt [3]{ab^{2} }=\sqrt [3]{\beta }\, D$ ) is used as the characteristic length to calculate the Reynolds number in some of the literature (Zhang, Ni & Magnaudet Reference Zhang, Ni and Magnaudet2021; Sanjeevi, Dietiker & Padding Reference Sanjeevi, Dietiker and Padding2022; Wang et al. Reference Wang, Yang, Andersson, Zhu, Wu, Wang and Liu2022; Shi, Zhang & Magnaudet Reference Shi, Zhang and Magnaudet2024), i.e. $ \textit{Re}_{d} = U_{\infty } d / \nu$ . We perform a series of numerical simulations using the OpenFOAM package, with the parameters varied as listed in table 1. The effects of different pitching frequencies on the hydrodynamic and flow-field characteristics of the prolate spheroid are considered. The $45^{\circ }$ initial inclination angle selected in this study serves as a typical threshold for boundary-layer separation in bluff body flows (Lim, Castro & Hoxey Reference Lim, Castro and Hoxey2007; Kumar et al. Reference Kumar, Sourav, Sen and Yadav2018; Zigunov, Sellappan & Alvi Reference Zigunov, Sellappan and Alvi2020). It represents a critical condition where small changes in vortex dynamics strongly influence unsteady loads. This configuration provides an opportunity to probe the nonlinear interplay between added-mass forces and viscous separation, making it an ideal test case for assessing hybrid models. In the low-frequency regime, the combination of a $5^{\circ }$ pitching amplitude and a non-dimensional pitching frequency ( $f=0{-} 0.16$ ) can simulate either the attitude fine-tuning of underwater vehicles or the pulsed pitching kinematics witnessed in biological propulsion systems.

Table 1.

List of cases investigated in this study.

Table 1. Long description

TA table with three columns and nine rows. The columns are labeled ‘Cases’, ‘Re_D’, and ‘f’. The ‘Cases’ column lists specific case numbers or ranges. The ‘Re_D’ column lists Reynolds numbers: 200 for case 1 and 50, 100, 150, 200, 250, 300 for cases 2-43. The ‘f’ column lists frequencies: 0 for case 1 and 0.032, 0.040, 0.053, 0.080, 0.106, 0.128, 0.160 for cases 2-43. The table organizes cases by ranges and provides corresponding Reynolds numbers and frequencies for each case.

Figure 2.

The unstructured overset mesh: three levels of local Cartesian mesh refinement, and a detailed view of the near-wall boundary layer and overset mesh.

An overset mesh approach (Brazell, Sitaraman & Mavriplis Reference Brazell, Sitaraman and Mavriplis2016; Morse & Mahesh Reference Morse and Mahesh2021; Lin et al. Reference Lin, Zhu, Liu and Wang2026) is used to handle the prolate spheroid’s pitching motion. The unstructured overset mesh with local mesh refinement is shown in figure 2. The size of the grid is multiplicatively related to the refinement levels of the grid. A locally refined mesh is employed, with minimum spacing ( $0.05D$ in all three directions) near the body surface (level 3), specifically within the prismatic mesh layers (Roget et al. Reference Roget, Sitaraman, Lakshminarayan and Wissink2020; Ye et al. Reference Ye, Liu, Ni and Chen2025). The dimensionless wall distance $y^{+} \lt 1$ ensures resolution of the viscous sublayer. Coarser spacing ( $0.4D$ in all three directions) is used in the far field (level 0) to reduce computational cost. The overset grid and the nearby background grid have the same mesh density (level 3) to facilitate data exchange at the overset boundaries. The total number of cells is approximately 17 million. By altering the size of the grid, a convergence test has confirmed the validity of the current set-up, as detailed in Appendix A.

2.3. Force and torque coefficients

The drag force ( $F_{x}$ ) acting upon a prolate spheroid is aligned with the direction of the flow and is characterised by the drag coefficient ( $C_{\kern-1pt D}$ ), which is defined as

(2.4) \begin{align} C_{\kern-1pt D} =\frac {F_{x}}{\frac {1}{2} \rho U_{\infty }^{2} \frac {\pi }{4} d^{2} } , \end{align}

where $\rho$ denotes the fluid density. Here, defining the reference area as the cross-sectional area of a volume-equivalent sphere ensures that it remains constant irrespective of the angle of incidence.

The lift force ( $F_{y}$ ) acts in a direction perpendicular to the upstream flow velocity and can be characterised by a lift coefficient ( $C_{L}$ ), which is defined as

(2.5) \begin{align} C_{L} =\frac {F_{y}}{\frac {1}{2} \rho U_{\infty }^{2} \frac {\pi }{4} d^{2} }. \end{align}

The pitching torque around the $Z$ -axis ( $M_{z}$ ) can be characterised by the torque coefficient ( $C_{T}$ ), which is defined as

(2.6) \begin{align} C_{T} =\frac {M_{z}}{\frac {1}{2} \rho U_{\infty }^{2} \frac {\pi }{8} d^{3} }. \end{align}

The time-averaged force and torque coefficients are defined as

(2.7a) \begin{align} \bar {C}_{D} = \frac {1}{\frac {1}{2} \rho U_{\infty }^{2} \frac {\pi }{4} d^{2} } \times\frac {1}{\textit{NT}} \int _{t_{0} }^{t_{0}+\textit{NT}} F_{x} (t) \,\mathrm{d}t, \\[-12pt] \nonumber \end{align}

(2.7b) \begin{align} \bar {C}_{L} = \frac {1}{\frac {1}{2} \rho U_{\infty }^{2} \frac {\pi }{4} d^{2} } \times\frac {1}{\textit{NT}} \int _{t_{0} }^{t_{0}+\textit{NT}} F_{y} (t) \,\mathrm{d}t, \\[-12pt] \nonumber \end{align}

(2.7c) \begin{align} \bar {C}_{T} = \frac {1}{\frac {1}{2} \rho U_{\infty }^{2} \frac {\pi }{8} d^{3} } \times\frac {1}{\textit{NT}} \int _{t_{0} }^{t_{0}+\textit{NT}} M_{z}(t) \,\mathrm{d}t, \\[10pt] \nonumber \end{align}

where $t_0$ denotes the starting instant for averaging, commencing once the flow has attained a statistically periodic state, and $N$ represents the number of pitching periods over which the averaging is performed. In this study, value $N = 5$ is adopted.

4.1. Effect of added mass

In § 3, we explore the time-averaged force coefficients associated with the pitching prolate spheroid. However, the pitching motion exhibits unsteady characteristics that necessitate additional scrutiny. These attributes complicate the mathematical description of the forces and torques acting on the prolate spheroid, as they vary significantly with time and position during each cycle of pitching. The unsteady effect arises from the transient dynamics of vortex shedding and boundary-layer separation, which occur at distinct phases of the pitching motion. To accurately model the prolate spheroid’s response, both viscous effects and added-mass effects on forces and torques must be accounted for: the latter represent inertial loads due to fluid acceleration in the body’s vicinity, even in the absence of viscosity.

Classical potential-flow theory provides analytical estimates of added mass and associated hydrodynamic coefficients for canonical geometries such as spheres and ellipsoids (Kochin et al. Reference Kochin, Kibel, Roze, Boyanovitch and Radok1964). These results establish that unsteady inertial forces originate from the acceleration of the surrounding ideal fluid, and that the added-mass tensor depends solely on body geometry. Building on these classical formulations, we derive the added-mass forces and torques for a pitching prolate spheroid at a large inclination angle.

The dimensionless results for the added-mass forces and torques can be found by referring to the forms of (2.4)–(2.6), respectively:

(4.1) \begin{align} \left . \begin{array}{lll} C_{D\hbar }=8F_{x\hbar }/ \rho \pi U_{\infty }^{2} d^{2} ,\\[6pt] C_{L\hbar }=8F_{y\hbar }/ \rho \pi U_{\infty }^{2} d^{2},\\[6pt] C_{T\hbar }=16M_{z\hbar }/ \rho \pi U_{\infty }^{2} d^{3}, \end{array}\right \} \end{align}

where $F_{x\hbar }$ , $F_{y\hbar }$ and $M_{z\hbar }$ denote the added-mass forces and torques in earth-fixed coordinates, which are obtained via the transformation from the body-fixed coordinates,

(4.2) \begin{align} \left . \begin{array}{lll} F_{x\hbar }=F_{\hat {x}}\cos \phi - F_{\hat {y}}\sin \phi = \left (\lambda _{22}-\lambda _{11}\right ) U_{\infty } \dot {\phi } \sin 2\phi ,\\[6pt] F_{y\hbar }=F_{\hat {x}}\sin \phi + F_{\hat {y}}\cos \phi = \left (\lambda _{11}-\lambda _{22}\right ) U_{\infty } \dot {\phi } \cos 2\phi ,\\[6pt] M_{z\hbar }=M_{\hat {z}}=-\lambda _{66} \ddot {\phi } + \frac {1}{2} \left (\lambda _{11}-\lambda _{22}\right ) U_{\infty }^{2} \sin 2\phi , \end{array}\right \} \end{align}

where $F_{\hat {x}}$ , $F_{\hat {y}}$ and $M_{\hat {z}}$ denote the added-mass forces and torques in body-fixed coordinates, $\lambda _{jk}$ ( $j, k = 1, 2$ ) are the added masses, $\dot {\phi } = \partial \phi / \partial t$ is the rotational angular velocity around the $Z$ -axis, and $\ddot {\phi } = \partial ^{2} \phi / \partial t^{2}$ is the rotational angular acceleration around the $Z$ -axis. The specific derivation process for these quantities can be found in the supplementary material.

Building upon the quasi-steady correlations established by Sanjeevi et al. (Reference Sanjeevi, Dietiker and Padding2022), this study first develops an enhanced hydrodynamic model by incorporating added-mass effects into the base formulation. Figure 9 shows comparative results between the improved model and classical quasi-steady formulations under various pitching frequency conditions when $ \textit{Re}_{D}=200$ . Here, hydrodynamic coefficients across three pitching cycles ( $T$ ) at two distinct non-dimensional frequencies are presented. Critically, the classical quasi-steady model exhibits demonstrable frequency independence, manifested through invariant phase relationships and amplitude of force and torque. Analytical results reveal that the classical quasi-steady model significantly underestimates the amplitude variations of drag, lift and torque coefficients compared to unsteady simulation results, with the most pronounced discrepancy observed in torque coefficient predictions – exhibiting a maximum NRMSE of 35.4 % as quantified in figure 10. When accounting for added-mass effects under unsteady flow conditions, the corrective influence on lift coefficient magnitude remains relatively constrained (relative error improvement rate $\lt 1\,\%$ ). However, notable suppression effects emerge in NRMSE metrics for both drag coefficient and pitching torque coefficient, with maximum reductions approaching 12 % at elevated frequencies. This frequency-dependent divergence originates from the anisotropic characteristics of added mass within relatively high-frequency domains: the phase-lag phenomenon in vortex shedding induced by fluid–structure interaction amplifies dynamic response deviations from quasi-steady solutions for force coefficients. Incorporation of added-mass parameters effectively compensates for such transient inertial effects, thereby attenuating NRMSE values in hydrodynamic load characteristics. In addition, implementation of the added-mass correction enhances temporal congruence with numerical simulations by resolving inertial phase lags inherent in unsteady flows. This is conclusively evidenced in figure 9 at $f=0.160$ , where the corrected lift coefficient exhibits significantly improved synchronisation with numerical results – particularly in peak magnitude alignment and phase correspondence. However, residual discrepancies persist in amplitude predictions, as evidenced by maximum NRMSE values 16.7 %, 24.2 % and 34.6 % for drag, lift and torque coefficients, respectively. This deviation originates from the inherent limitations of the current potential flow-based added-mass formulation, which fails to account for nonlinear coupling between vortical evolution and solid-body kinematics. Subsection 4.2 endeavours to focus on incorporating a viscous correction factor to refine the theoretical framework.

Figure 9.

Comparative time-history analysis of simulation data and correlations of the drag, lift and torque coefficients considering the effect of added mass on the prolate spheroid at various pitching frequencies when $ \textit{Re}_{D}=200$ .

Figure 10.

The NRMSE of the simulated values versus the correlations for drag, lift and torque coefficients considering the effect of added mass at different pitching frequencies when $ \textit{Re}_{D}=50$ and 300.

Figure 10. Long description

Two line graphs compare the normalized root mean square error (NRMSE) of simulated values for drag, lift, and torque coefficients at different pitching frequencies for two Reynolds numbers. Panel A: The left graph shows data for a Reynolds number of 50. The x-axis represents the pitching frequency (f) ranging from 0.02 to 0.17, and the y-axis represents the NRMSE in percent, ranging from -5 to 40 percent. The graph includes multiple data series represented by different symbols and colors: squares for drag coefficients, circles for lift coefficients, and triangles for torque coefficients. Each series is shown with and without added mass effects. Panel B: The right graph shows data for a Reynolds number of 300. The x-axis represents the pitching frequency (f) ranging from 0.02 to 0.17, and the y-axis represents the NRMSE in percent, ranging from -5 to 40 percent. Similar to Panel A, the graph includes multiple data series represented by different symbols and colors: squares for drag coefficients, circles for lift coefficients, and triangles for torque coefficients. Each series is shown with and without added mass effects. Both graphs indicate that the NRMSE varies with pitching frequency and Reynolds number, with different trends observed for drag, lift, and torque coefficients.

Figure 11.

Comparative analysis of angular acceleration phase diagrams of simulation data and correlations of the drag, lift and torque coefficients considering the effect of viscosity on the prolate spheroid at various pitching frequencies when $ \textit{Re}_{D}=200$ .

Figure 12.

Comparative analysis of angular acceleration phase diagrams of simulation data and correlations of the drag, lift and torque coefficients considering the effect of viscosity on the prolate spheroid at various Reynolds numbers when $f=0.080$ .

Figure 12. Long description

The image contains nine panels arranged in a 3x3 grid, each depicting phase diagrams of angular acceleration for different Reynolds numbers. Each panel compares simulation data with correlation data for drag, lift, and torque coefficients. Panel A: The top row shows the drag coefficient (C_Dv) versus angular acceleration (phi dot) for Reynolds numbers 50, 150, and 250. The red circles represent simulation data, and the blue line represents correlation data. Panel B: The middle row shows the lift coefficient (C_Lv) versus angular acceleration (phi dot) for the same Reynolds numbers, with the same color scheme. Panel C: The bottom row shows the torque coefficient (C_Tv) versus angular acceleration (phi dot) for the same Reynolds numbers, again with the same color scheme. Each panel demonstrates how the coefficients vary with angular acceleration and Reynolds number, highlighting the differences between simulation and correlation data.

4.2. Effect of viscosity

Fluid viscosity introduces additional forces and modifies the flow patterns around the prolate spheroid. Viscosity affects the drag, lift and torque experienced by the prolate spheroid, and thus impacts the force coefficients $C_{i}$ , where $i$ is 1 for drag ( $C_{\kern-1pt D}$ ), 2 for lift ( $C_{L}$ ), and 3 for torque ( $C_{T}$ ). The presence of viscosity means that the theoretical predictions based on potential flow conditions will diverge from the actual simulation outcomes, which are conducted in a viscous environment. To bridge this gap and provide a more accurate representation of the physical phenomena, an adjustment term is introduced, denoted $C_{iv}$ . This term is specifically designed to reconcile the discrepancies between the actual coefficients $C_{i}$ obtained from the simulation and the sum of their time-averaged quasi-steady counterparts $\bar {C}_{i, \phi }$ and the unsteady added-mass contributions $C_{i\hbar }$ (see (4.2)). This study investigates the temporal evolution characteristics of force coefficients within the pitching parameter space by quantitatively characterising the combined effects of added mass and viscosity on hydrodynamic behaviour. Consequently, incorporating the viscosity correction term to achieve a more accurate representation, the final correlations of the unsteady force coefficients ( $C_{i}$ ) can be expressed as

(4.3) \begin{align} C_{i} = \underset {\text{quasi-steady term}}{\underbrace {\bar {C}_{i, \phi }}} + \underset {\text{unsteady term} }{\underbrace {C_{i\hbar }}} + \underset {\text{viscosity correction} }{\underbrace {C_{iv}} } . \end{align}

In this subsection, we propose new correlations for the unsteady drag, lift and torque coefficients considering the effect of viscosity on the prolate spheroid across various pitching frequencies and Reynolds numbers, as illustrated in figures 11 and 12. The simulation results show how $C_{Dv}$ , $C_{Lv}$ and $C_{Tv}$ correlate with the angular velocity function $\dot {\phi }$ at different frequencies $f$ and Reynolds numbers ( $ \textit{Re}_{D}$ ). Each row in the figure corresponds to one of these hydrodynamic coefficients, while each column represents a distinct frequency or Reynolds number. The functional form in (4.4) below is motivated by the observed elliptical feature in the phase-space trajectories of the viscous force and torque coefficients. As shown in figures 11 and 12, the hysteresis loops traced by these coefficients over a pitching period strongly resemble closed, approximately elliptical contours – a signature of quadrature-phase coupling between body motion and fluid response. This elliptical behaviour indicates that the viscous force correction acts as a linear damping term within a harmonic oscillator framework, a hallmark of vortex-induced forces under frequency lock-in conditions (Bearman Reference Bearman1984; Williamson & Govardhan Reference Williamson and Govardhan2004). Such behaviour is characteristic of weakly nonlinear, periodically forced systems near resonance. The adopted general elliptic regression model thus captures the amplitude and phase effect of the unsteady viscous loads, providing a compact, physics-informed empirical representation suitable for reduced-order modelling. This phenomenological approach aligns with the idea of the Morison equation decomposition (Sarpkaya Reference Sarpkaya2010), where the unsteady force is partitioned into inertial and viscous components. The phase lag introduced by the trigonometric terms in (4.4) is analogous to the effects captured by unsteady hydrodynamic theories such as Theodorsen’s function (Fung Reference Fung2008; Cordes et al. Reference Cordes, Kampers, Meißner, Tropea, Peinke and Hölling2017), which reflects the influence of the wake vorticity on the force phase.

Based on the distribution characteristics observed in the simulation results, nonlinear regression analyses are conducted using elliptic equations for $C_{Dv}=C_{1v}$ , $C_{Lv}=C_{2v}$ and $C_{Tv}=C_{3v}$ . The formally consistent correlations for $C_{iv}$ ( $i = 1, 2, 3$ ) are given by

(4.4) \begin{align} C_{iv} = - \dot {\phi } \sin \varsigma _{i1} - \mathrm{sgn} \left ( \ddot {\phi } \right ) \frac {\varsigma _{i2}\sqrt {1- \dot {\phi } ^{2}/\phi _{f}^{2}}}{ \cos \varsigma _{i1} }+\varsigma _{i3}, \end{align}

where $\mathrm{sgn} ( \ddot {\phi } )$ denotes the sign function corresponding to the upper and lower halves of the pitching period, $\phi _{f} = 2 \pi f \phi _{m}$ is the amplitude of angular velocity, and $\varsigma _{\textit{ij}}$ ( $i, j = 1, 2, 3$ ) denote the fitting coefficients of the hydrodynamic coefficients taking into account the effect of viscosity. Specifically, $\varsigma _{i1}$ primarily reflects phase shift to adjust the orientation of the $C_{iv}$ – $\dot {\phi }$ phase diagram, $\varsigma _{i2}$ reflects viscous damping magnitude, and $\varsigma _{i3}$ captures baseline offset due to mean shear stresses.

In periodic motion, the formation and shedding of vortices do not occur instantaneously with the object’s position, but lag by a certain phase. The parameter $\varsigma _{i1}$ corresponds to the time delay of vortex-induced forces, analogous to the phase difference between damping forces and displacement in a damped oscillatory system. It accounts for the time delay between the instantaneous angular motion of the prolate spheroid and the formation and shedding of vortical structures, and therefore determines the inclination of the elliptical phase-plane trajectories. For example, (4.4) can simplify to a standard elliptic equation when $\varsigma _{i1}=0$ , namely

(4.5) \begin{align} \frac {\dot {\phi } ^{2}}{\phi _{f}^{2}}+\frac {\left ( C_{iv}- \varsigma _{i3}\right ) ^{2} }{\varsigma _{i2}^{2} }=1. \end{align}

The parameter $\varsigma _{i2}$ reflects the amplitude of the viscous contribution, corresponding to the size of the ellipse’s major axis in the $C_{iv}$ – $\dot {\phi }$ phase diagram, and can be interpreted as a viscous damping coefficient. It characterises the strength of the separated shear layers and hairpin vortices, which directly influence the surface pressure distribution. As the Reynolds number and pitching frequency increase, the parameter $\varsigma _{i2}$ grows, reflecting the intensification of unsteady vortex structures, thus quantifying the effective viscous damping.

The parameter $\varsigma _{i3}$ represents a bias term corresponding to the mean offset of the viscous correction. This arises from the asymmetry of boundary-layer separation between the upward and downward parts of the pitching cycle, and encapsulates the residual mean force produced by viscous shear.

Nonlinear distortions at higher frequencies are empirically captured through the coefficients $\varsigma _{\textit{ij}}$ , a well-established strategy for modelling complex fluid–structure interactions (Holmes Reference Holmes2012). To satisfy the fitting accuracy of nonlinear equations in two variables, the Taylor series rational function is used to construct the regression equation between the coefficients $\varsigma _{\textit{ij}}$ and the pitching frequencies and Reynolds numbers given by

(4.6) \begin{align} \varsigma _{\textit{ij}} \left ( \breve {f} , \breve {Re_{D}} \right ) = \frac { \backepsilon _{ij1} + \backepsilon _{ij2} \breve {f} + \backepsilon _{ij3} \breve {Re_{D}} + \backepsilon _{ij4} \breve {Re_{D}}^{2} + \backepsilon _{ij5} \breve {f} \breve {Re_{D}} } { 1 + \backepsilon _{ij6} \breve {f} + \backepsilon _{ij7} \breve {f}^{2} + \backepsilon _{ij8} \breve {Re_{D}} + \backepsilon _{ij9} \breve {Re_{D}}^{2} + \backepsilon _{ij10} \breve {f} \breve {Re_{D}} }, \end{align}

where $\backepsilon _{ijk}$ ( $i, j = 1, 2, 3$ and $k=1{-}10$ ) denote the fitting coefficients for the correlations of unsteady force and torque coefficients, as shown in table 3. To examine whether a simpler representation is sufficient, we also tested a two-dimensional polynomial with six parameters. As detailed in the supplementary material, this reduced model exhibits markedly larger residuals and phase errors, particularly for $C_{Lv}$ and $C_{Tv}$ at higher Reynolds numbers. By contrast, the 10-parameter form achieves consistently lower Akaike information criteria ( $AICc$ ) and a higher R-squared ( $R^{2}$ ) value, confirming it as the most reliable approximation.

Table 3.

Summary of fitting coefficients for the correlations of unsteady force coefficients using Taylor series rational functions.

Table 3. Long description

The table presents fitting coefficients for the correlations of unsteady force coefficients using Taylor series rational functions. It has 10 columns labeled k = 1 to k = 10 and 9 rows labeled with subscripts. Each cell contains a numerical value representing the fitting coefficients. The table is structured to show the relationship between the coefficients and the pitching frequencies and Reynolds numbers. Row 1: backepsilon 11k, -0.071, 0.297, 0.164, 0.047, 0.767, 3.838, -0.654, 1.793, 42.561, 30.202. Row 2: backepsilon 12k, 0.032, 0.017, 0.047, -0.046, 0.021, -0.537, 0.317, 0.811, -0.942, 0.476. Row 3: backepsilon 13k, 0, 0.006, 0.001, -0.003, 0.047, -1.556, 1.156, 1.384, 1.046, -1.103. Row 4: backepsilon 21k, -0.104, -0.371, -1.765, -4.983, 9.181, 5.726, 39.559, -72.515, 876.769, 148.870. Row 5: backepsilon 22k, 0.023, 0.011, -0.022, 0.019, -0.015, -0.576, 0.962, -0.701, 0.899, -1.228. Row 6: backepsilon 23k, 0, 0.007, 0.001, -0.004, 0.055, -1.265, 0.775, 3.682, 1.110, -1.815. Row 7: backepsilon 31k, 0.827, 2.301, 12.996, 0.776, -6.151, 2.677, -0.214, 20.177, 82.541, -1.642. Row 8: backepsilon 32k, 0.024, 0.271, -0.064, 0.099, -0.381, 2.527, -2.096, -2.488, 2.899, 0.567. Row 9: backepsilon 33k, 0, -0.005, -0.001, 0.002, -0.021, -2.094, 1.723, 1.783, 0.928, -2.208.

The phase diagrams shown in figures 11 and 12 serve to illustrate the periodic structure of the viscous component and its dependence on angular motion, thereby clarifying the physical basis of the proposed correlation. For the $C_{Tv}$ – $\dot {\phi }$ phase diagram, the elliptical form is consistently preserved at different values of $f$ and $ \textit{Re}_{D}$ , indicating strong periodic behaviour. Consequently, the correlation given by (4.3) is closely aligned with the simulation data. The phase diagrams for both drag and lift exhibit predominantly elliptical characteristics at frequencies not exceeding 0.1, as shown in figure 11. This behaviour correlates well with numerical simulations, showing satisfactory agreement with theoretical expectations. As the frequency $f$ approaches 0.16, nonlinear effects emerge, accompanied by phase shift phenomena. Consequently, the phase diagrams develop asymmetry – particularly observable in the lift phase diagram – owing to the growing influence of higher harmonic components and vortex-induced flow separation. This transformation suggests varying flow phenomena at higher frequencies, consistent with the alterations in wake vortex patterns observed at corresponding frequencies in figure 5. Notwithstanding these nonlinear distortions, the fundamental elliptical feature persists as an approximation. This suggests that lower-order elliptical approximations retain utility for preliminary performance estimates even within weakly nonlinear flow conditions. Figure 12 shows the results of the hydrodynamic coefficients under different Reynolds numbers. Crucially, whilst mean values display significant Reynolds number dependency, the topological characteristics of all hydrodynamic coefficients phase diagrams maintain similarity. The lift coefficient manifests subtle deviations attributable to nonlinear effects at increased pitching frequencies. However, the formally consistent correlation (see (4.4)) establishes a robust predictive framework capturing the fundamental lift response. Since the correlations were derived via regression against the same dataset, these figures primarily illustrate fit quality rather than independent predictive validation. Nevertheless, the close agreement across a wide range of frequencies and Reynolds numbers supports the functional form’s ability to capture the dominant unsteady physics. This is substantiated by validation against time histories, as presented in figures 13 and 14, which compare hydrodynamic coefficient evolution before and after introducing the viscous correction term $C_{iv}$ across varying pitching frequencies and Reynolds numbers. Post incorporation of the $C_{iv}$ term, all hydrodynamic coefficients exhibit significantly enhanced congruence with numerical data, showing quantifiable improvements in both amplitude fidelity and phase synchronisation.

Figure 13.

Comparative time-history analysis of simulation data and correlations of the drag, lift, and torque coefficients considering the effect of added mass and viscosity on the prolate spheroid at various pitching frequencies when $ \textit{Re}_{D}=150$ .

Figure 14.

Comparative time-history analysis of simulation data and correlations of the drag, lift, and torque coefficients considering the effect of added mass and viscosity on the prolate spheroid at various Reynolds numbers when $f=0.106$ .

Quantitative validation against unsteady numerical simulations reveals that incorporating viscous dissipation mechanisms achieves significant improvements in force and torque coefficient prediction accuracy, with the NRMSE of torque coefficient showing a notable reduction of approximately 33 %, as shown in figures 10 and 15. These findings corroborate the critical role of viscous effects in unsteady flow regimes. From a quantitative perspective, figure 15 illustrates the NRMSE between the simulated values and the correlations for hydrodynamic coefficients, considering the effect of added mass and viscosity across various pitching frequencies and Reynolds numbers. A lower NRMSE with an average value less than 1.8 % signifies a better fit of the correlation given by (4.3) to the simulation results. Specifically, when $0.04 \leqslant f \leqslant 0.1$ and $50 \leqslant Re_{D} \leqslant 250$ , the maximum NRMSE values for drag, lift and torque coefficients are 0.7 %, 4.5 % and 2.0 %, respectively. For the torque coefficient, the NRMSE remains relatively stable throughout the $f{-}Re_{D}$ range, with a minor increase at lower frequencies and Reynolds numbers. The average NRMSE remains near 0.7 % for drag coefficient, and 1.5 % for torque coefficient. The developed analytical framework shows particular efficacy in modelling ellipsoidal bodies executing low-frequency pitch oscillations, as evidenced by progressive error amplification in drag and lift coefficients at elevated pitching frequencies and Reynolds numbers. The pronounced gradient is observed in lift coefficient profiles, inducing growing deviations from base benchmarks. Such frequency-dependent discrepancies originate from the framework’s inherent limitation in resolving flow separation mechanisms during high-rate pitch manoeuvres (see figure 6). Overall, the lift coefficient exhibits the most pronounced variation with frequency and Reynolds number, indicating complex flow phenomena that are not fully captured by the current correlation. This suggests a need for further refinement of the correlation to achieve accurate predictions at higher frequencies. The NRMSE of nearly 10 % is observed in lift coefficient prediction error within characteristic phase intervals when the dimensionless frequency reaches 0.160. This phenomenon exhibits intrinsic correlation with phase portrait analysis conclusions presented in figure 11, where intensified flow separation effects induce progressive divergence between vortical structure evolution patterns and quasi-steady assumptions. Specifically, this deviation manifests as phase-shifted dynamic vortex formation and discrepancies in vortex intensity magnitude. Such nonlinear unsteady effects are further evidenced in the temporal response characteristics through an amplitude difference phenomenon at lift coefficient extreme points, as illustrated in figure 13. In contrast, drag and torque coefficients show less variability, implying that their respective correlations are more robust across the examined frequency range.

Figure 15.

The NRMSE of the simulated values versus the correlations for drag, lift and torque coefficients considering the effect of added mass and viscosity at different pitching frequencies and Reynolds numbers.

Furthermore, a randomised sampling approach was applied within the Reynolds number–pitching frequency parameter space, generating five additional simulation cases beyond the original dataset. These cases were specifically designed to cover both intermediate and extrapolated regimes, thereby providing a stringent test of the robustness and general applicability of the proposed correlation model. As shown in figure 16, the correlation predictions match the numerical results across all newly simulated conditions. The NRMSE for all force coefficients does not exceed 10 % over the full range of Reynolds numbers and pitching frequencies, confirming that the correlation maintains predictive accuracy even under flow conditions not included in the original training data. It is noted, however, that at relatively high pitching frequencies, the NRMSE of the force coefficients increases, with the lift coefficient exhibiting the most pronounced deviation.

Figure 16.

Validation of the proposed correlations is performed using five supplementary simulations obtained by interpolation and extrapolation in the $ \textit{Re}_{D}{-}f$ parameter space. (a) Locations of the additional validation cases, designed to sample intermediate and extreme regions beyond the original dataset. (b) The NRMSE of the predicted drag ( $C_{\kern-1pt D}$ ), lift ( $C_{L}$ ) and torque ( $C_{T}$ ) coefficients for each case. (c–f) Time-history comparisons between simulation results and correlation predictions for the force coefficients at (c,d) $f=0.145$ , $ \textit{Re}_{D}=275$ , and (e,f) $f=0.028$ , $ \textit{Re}_{D}=175$ .