Generalized Dynamical Model for the Drag Coefficient of Fixed and Moving Cylindrical Bodies

1. Introduction

The flow over a cylindrical body is commonly encountered in different applications, such as offshore oil risers [1,2], wind turbine towers [3,4,5], solar updraft systems [6,7], overhead power distribution poles [8,9,10], chimneys of power plant and combustion processes [11,12,13], above-ground piping network of gases or liquids [14,15,16], and various structural elements in the buildings sector or the energy sector [17,18,19]. Therefore, theoretical investigations of such a problem have been conducted by several researchers [20,21,22]. Reduced (simplified or low-dimensional) modeling [23,24] allows for faster and less-demanding analysis of this problem, and thus can facilitate further control and design computations [25].

There are several reduced-order models (ROMs) [26,27] for the lift and drag coefficients for a fixed cylinder. These lift and drag coefficients are nondimensionalized forms of the two orthogonal components of the exerted hydrodynamic force due to the time-dependent variations in the near wake. As in other fluid dynamics or gas dynamics applications [28,29], the lift is the force component acting on the cylinder in the cross-stream direction, while the drag is the force component acting on the cylinder in the free-stream direction [30,31]. For convenience and simplified discussion in the current work, we may use the term “lift” or “lift force” to refer to the lift coefficient, and use the term “drag” or “drag force” to refer to the drag coefficient. This should not pose a fundamental concern since either force coefficient is derived from its relevant force component using the same scaling factor [32,33,34], which is ½ ρ U² D L, where ½ ρ U² is the dynamic pressure, and D L is the projected area of the cylinder perpendicular to the incoming flow [35,36]. While some of the recent phenomenological reduced-order models (i.e., models that explain one or more observed phenomena) and fully analytical reduced-order models (i.e., models whose parameters are fully related to some variables of the flow case of interest and thus can be systematically estimated) have proved to be very accurate and can even model the transient phase, in addition to the steady-state phase of development of these forces very well, they need to be extended to the moving-cylinder case. In fact, some early attempts to do that go back to 1970 [37,38], and there are existing models that represent the induced lift force for a moving cylinder in a free stream direction (in-line motion) [39,40]. Several of these models focus on the case of a cylinder that moves in the cross-stream direction [41,42], and these models couple this motion to the lift oscillator, which is a nonlinear ordinary differential equation (ODE) [43] representing a self-excited dynamical system [44,45,46], without paying attention to the drag [47,48,49,50,51,52,53,54,55].

Other models that consider the drag force as well as the lift impose a specific coupling between this drag force and the lift force, which is purely quadratic. Another example is the work based on modeling the structural response of elastic cable suspensions, where a model of the drag was introduced and the drag was represented as a function of time and the cable’s longitudinal coordinate in the static condition [56]. In that work, the drag was decomposed into the product of a temporal (time-dependent) function and a spatial (space-dependent) function. The time-dependent function is modeled by a van der Pol equation (cubic nonlinear oscillator) whose linear natural frequency is twice the shedding frequency (which was assumed to be equal to the angular natural frequency of the cable). As in drag-modeling studies for a fixed cylinder, that research work used quadratic coupling between the lift (which is modeled by another van der Pol equation whose linear natural frequency is equal to the shedding frequency) and the drag. Therefore, the two-to-one frequency relationship between the lift and the drag (the fundamental frequency of the oscillatory drag signal is twice the fundamental frequency of the oscillatory lift signal) was assumed in that work. While this assumption is usually acceptable for a wide range of recorded data of VIV (vortex-induced vibration), we show here in the current study that this assumption is not accurate for some types of flows, such as when the cylinder is vibrating at a prescribed tilt angle relative to the free stream [57,58]. We denote such cases when the two-to-one frequency relationship between the lift and drag is not satisfied as “non-traditional wakes”; this is in contrast to the case when this relationship is satisfied, which we denote by “traditional wakes”. Also, the discussed semi-empirical drag model included two empirical constants that were selected and were not well related to the data of interest being modeled (the training dataset). The proposed reduced-order drag model here is fully analytical, meaning that we provide closed-form mathematical expressions to relate its parameters to some variables that are extracted from the training data to be modeled such that good agreement between the reduced-order model (ROM) results and these training data is achieved.

A research team also worked on modeling both the lift and drag forces for fixed and moving cylinders with two degrees of freedom [59]. Instead of using a wake oscillator, they followed a force-decomposition approach, where the lift was assumed to be composed of a main frequency component at the shedding frequency and the odd harmonics, whereas the drag was assumed to be composed of a mean value in addition to a main frequency component at twice the shedding frequency and the even harmonics. For a fixed cylinder, the harmonics were neglected, and one-term lift and one-term drag (in addition to the mean-drag value) were used to represent the signals (time series) of these two fluid forces, which are used to excite the cylinder and move it in both the x and y directions. An earlier reduced-order model (ROM) accounted for the third harmonic of the lift when modeling the lift induced on a fixed cylinder. That drag model is consistent with the one used by another research team in terms of using one harmonic only for the drag of a fixed cylinder. For a moving cylinder, the work described claimed that the harmonics can still be dropped if the shedding frequency is away from the natural frequency of the cylinder, but these harmonics should be retained if the two frequencies are close to each other (i.e., at resonance [60,61]). The coefficients of these harmonics can be obtained from any arbitrary available time-history lift and drag data using auto-regressive moving averaging (ARMA) [62,63,64].

We contribute to this class of problems by extending existing drag reduced-order models (ROMs) that are based on a simple quadratic relationship between the lift and drag, to situations where this relationship does not hold. One of these situations is demonstrated here, which occurs when the cylinder vibrates along directions other than the cross-stream. This can be observed in fluid-structure interaction (FSI) [65,66] or vortex-induced vibration (VIV) [67,68,69] problems, where the oscillatory motion of the cylinder (or cylinder-like element) cannot be restrained to the perpendicular cross-stream direction [70,71,72].

Figure 1 illustrates this generalized situation, where the cylinder undergoes an oscillation that is inclined (i.e., lies between the in-line direction and the cross-stream direction).

The change in the lift-drag relationship is accompanied by other changes in their spectral characteristics. However, these changes are significant for the drag, and thus the available models fail in these cases. When the traditional wake regime and its implied simple two-to-one lift-drag relationship and spectra are restored, the proposed reduced-order model (ROM) [73,74,75] automatically converges to a form that is identical to traditional existing drag models (the extra terms we add in the proposed drag ROM automatically vanish).

Aside from the technical aspect, the current study can be viewed an educational asset and a method for teaching-research nexus [76,77,78]. It exposes college students (at both the undergraduate and the postgraduate levels) to an application of nonlinear dynamics, numerical simulations, and fluid mechanics [79,80,81].

2. Flow Analysis

In order to provide a training dataset of fluid force signals with a drag-lift relationship that does not obey the traditional two-to-one frequency relationship (we need such a dataset to develop and test the new universal reduced-order drag model), and compare its performance to existing drag models that are limited to the case of two-to-one frequency relationship between the lift and drag; we generate this training dataset using computational fluid dynamics (CFD) modeling [82,83,84,85]. For this part, we numerically solve (integrate in time and space) the two-dimensional incompressible Navier-Stokes (NS) equations [86,87,88], which are

$$\begin{array}{c}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}=0\\ {\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{\partial u^2}{\partial x}}+{\displaystyle \frac{\partial uv}{\partial y}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}+\nu {\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+\nu {\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}\\ {\displaystyle \frac{\partial v}{\partial t}}+{\displaystyle \frac{\partial uv}{\partial x}}+{\displaystyle \frac{\partial v^2}{\partial y}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial y}}+\nu {\displaystyle \frac{{\partial }^{2}v}{\partial x^2}}+\nu {\displaystyle \frac{{\partial }^{2}v}{\partial y^2}}\end{array}$$

(1)

In the above equation, u and v are the axial (streamwise) and vertical (cross-flow) components of the fluid flow velocity vector in the x and y coordinates, respectively. The symbol t designates the time, ρ is the fluid density, p is the fluid pressure, and ν is the fluid’s kinematic viscosity. Due to the lack of heat sources or sinks [89], heat transfer mechanisms, thermal radiation [90], temperature-dependent fluid properties [91,92], and high-speed compressible-fluid viscous heating [93], the energy equation is not solved [94,95,96]. This implies a reasonable assumption of isothermal (constant-temperature) cold flow. This is considered a reasonable assumption.

Before the numerical solution, these NS equations are nondimensionalized using the free-stream velocity U, the cylinder diameter D, and twice the dynamic pressure ρU². We then discretize and solve the resulting system of partial differential equations (PDEs) in the time domain using direct numerical simulation (DNS) [97,98]. We then compute the unsteady flow field around the vibrating cylinder, including the pressure distribution and the shear-stress distribution on its surface.

Although the lift coefficient (C_L) and the drag coefficient (C_D) for a cylinder should be theoretically obtained through integrating the nondimensional forces over the cylinder’s surface (over circular circumference in the current two-dimensional problem), in computational fluid dynamics (CFD) modeling [99,100], such continuous symbolic integration needs to be discretized and replaced by an quantitative summation process [101]. This transformation is needed because a closed-form analytical solution does not exist for the lift coefficient or the drag coefficient for a cylinder subject to a viscous fluid flow 102,103, although it exists in a very simplified case of potential flows where the viscosity is totally neglected [104], which is referred to as an inviscid flow [105,106], governed by Euler equations [107,108]. Thus, numerical approximation is necessary, and it is achieved by representing the perfectly circular circumference of the cylinder as a series of connected straight-line segments. Each segment is formed by connecting two adjacent grid points located on the cylinder’s surface (circumference). The incremental distances Δx and Δy are the horizontal and vertical projections of the distance between two adjacent grid points on the cylinder surface. Figure 2 illustrates a portion of the cylinder surface between two grid points and shows Δx and Δy.

We numerically integrate the nondimensional pressure coefficient and the nondimensional friction coefficient over the cylinder surface to compute the corresponding nondimensional lift and drag coefficients, which are recorded as two time signals according to

$$\begin{array}{l}{C}_L\left(t\right)=-{\displaystyle \sum _{\mathrm{surface}}\left[C_P\left(t\right)\hspace{0.33em}\Delta \hspace{0.17em}x\right]}+{\displaystyle \sum _{\mathrm{surface}}\left[C_f\left(t\right)\hspace{0.33em}\Delta \hspace{0.17em}y\left(t\right)\right]}\\ C_D\left(t\right)=-{\displaystyle \sum _{\mathrm{surface}}\left[C_P\left(t\right)\hspace{0.33em}\Delta \hspace{0.17em}y\right]}-{\displaystyle \sum _{\mathrm{surface}}\left[C_f\left(t\right)\hspace{0.33em}\Delta \hspace{0.17em}x\left(t\right)\right]}\end{array}$$

(2)

where C_P is the time-dependent pressure coefficient (or nondimensional pressure). It is the gauge pressure (the absolute pressure minus the freestream ambient or atmospheric pressure) divided by the dynamic pressure as a scaling factor [109,110,111]. Similarly, Cf is the time-dependent friction coefficient (or nondimensional shear stress). It is the shear stress divided by the dynamic pressure as a scaling factor.

At the present Reynolds number (Re) of 300, the pressure contribution (the first term on the right-hand side of Eq. 2) to the lift and drag coefficients is more significant than the friction contribution (the second term on the right-hand side of Eq. 2). However, we did not neglect the friction term when computing the lift and drag coefficients. It is useful to mention here that the pressure contribution to the lift and drag coefficients arises from the asymmetric pressure distribution over the cylinder’s surface [112]. Vertical asymmetry (deviation between the pressure distribution over the top half of the cylinder and the bottom half of the cylinder) contributes to the lift coefficient, while horizontal asymmetry (deviation between the pressure distribution over the front, upstream half of the cylinder and the rear, downstream half of the cylinder) contributes to the drag coefficient. This pressure asymmetry, in turn, is primarily caused by the low-pressure wake zone formed downstream of the cylinder due to eddies (vortices) [113]. A vortex has a low-pressure zone at its core [114] due to the high spinning speeds developed there, in accordance with Bernoulli’s principle [115] (energy conservation) [116], where pressure energy is converted into kinetic energy, as well as due to the centrifugal effects that pull the spinning fluid outward and establish a suction power. It is also useful to elaborate on the flow regimes over a stationary cylinder based on the Reynolds number. At extremely low Reynolds numbers (Re < 5), reflected by a very slow incoming flow, the flow can be described as “creeping” [117,118]. For such a limiting case, the flow around the cylinder is nearly double-symmetric (having both front-rear symmetry as well as top-bottom symmetry), and no separation of the boundary layer occurs. The creeping flow is laminar and steady, where no irregularities or unsteadiness are observed. This nearly double-symmetric flow regime and double-symmetric pressure distribution cause a cancellation of the pressure effects, leading to vanishing contributions to either the lift coefficient or the drag coefficient. However, the skin friction effect still results in a non-zero drag coefficient. As the Reynolds number increases up to around 40, the boundary layer separates at the downstream side of the cylinder’s surface, and two small vortices (Fopple vortices) develop [119,120], but they remain confined within the near-wake zone rather than being shed away downstream. Thus, this flow regime features a static (steady) separated near-wake region. As the Reynolds number increases further up to around 200, a laminar periodic vortex street is established in the wake. This oscillatory pattern of the flow around the cylinder in its near wake zone causes the lift and drag coefficients to also exhibit an oscillatory behavior. The laminar vortex street gradually transitions into a turbulent one until a fully turbulent vortex street is reached, approximately at a transition threshold Reynolds number [121,122] or critical Reynolds number [123] of 3 × 10⁵ [124]. Up to this Reynolds number, the boundary layer over the unseparated part of the cylinder’s surface is laminar. Laminar boundary layers exhibit relatively weak resistance to separation compared to turbulent boundary layers. This causes the laminar boundary layer to separate early [125], and this enlarges the width of the low-pressure wake. This magnifies the imbalance in the oscillatory pressure distribution, leading to a large contribution to the drag and the lift coefficients. The portion of the drag caused by the pressure imbalance is called form drag or pressure drag. On the other hand, the skin friction effect of a laminar boundary layer is smaller than the skin friction effect of a turbulent boundary layer. This makes the frictional contribution to the drag small compared to the pressure when the Reynolds number is below the transition threshold. As the Reynolds number increases beyond this transition threshold, the boundary layer undergoes a gradual change from a laminar type to a turbulent type, and eventually becomes a turbulent boundary layer. This retards the boundary layer separation, causing the two separation points to move further downstream along the cylinder’s surface, leading to a narrower zone of low-pressure wake. This suppresses the pressure imbalance, and this weakens the influence of the pressure distribution on the lift and drag coefficients. On the other hand, the larger skin friction influence of the turbulent boundary layer increases its contribution to the lift and drag coefficients.

The cylinder undergoes a sinusoidal vibration with a frequency that is equal to the Strouhal frequency (i.e., Ω = ω_S) [126,127], and with a prescribed vibrational amplitude. The cylinder’s vibrational motion starts from the beginning of the CFD simulation.

The values of C_P and C_f are resolved at the grid points; then, an average value is considered to be applied at the midpoint. These midpoint values are used in Eq. 2 to find C_L and C_D at every time step after the flow field is resolved. The arbitrary Lagrange Euler (ALE) technique was used to implement the cylinder motion and couple it to the flow solver [128,129,130,131]. In this ALE technique, Eq. 1 is written in an inertial coordinate system (having a global fixed center). Initially, the local cylinder center coincides with this global center. The grid deforms every time step, and it is re-generated algebraically to fit the new fluid region.

Although numerical modeling inevitably involves approximations and errors [132,133], a rigorous numerical approach was followed to ensure that the resolved flow features are physical and do not reflect large numerical errors. Several test cases were performed, which include solving the flow field past a fixed cylinder and a moving cylinder under different simulation parameters (e.g., time step, grid resolution, and grid size) so that we get converged solutions that are satisfactorily independent of these preset parameters that are manually selected rather than automatically derived.

In order to validate the adopted numerical simulation procedure, we performed several calculations in order to reproduce other reported benchmarking results [134] from numerical or experimental studies for both fixed and moving cylinders. For example, we ran a test case at Reynolds number 500 for both fixed and moving cylinders and compared the results with a similar two-dimensional direct numerical simulation of an earlier independent study [135], where they used the spectral element method with ninth-order Gauss-Lobatto-Legendre polynomials [136]. For the fixed cylinder, we got fS = 0.225 (compared to 0.228), time-averaged CD = 1.429 (compared to 1.46), and an amplitude of CL = 1.159 (compared to 1.20). Our values for this case also agree with the measured value of fS = 0.22 [137] and the measured value of time-averaged CD = 1.5 [138]. We also solved the case of a moving cylinder in the cross-stream direction at different frequencies and got similar results as reported by our benchmarking study [135]. In Figure 3, we show the agreement of our results with the reported ones in terms of the limit cycle projection onto the CL–η plane at motion frequencies 0.875 fS and 0.975 fS. The lift response is periodic in both cases.

3. Model Analysis

We start this section by illustrating the quadratic coupling between the drag and lift for the fixed cylinder. Figure 4 shows typical power spectra of the lift and drag coefficients for this case, which we obtained from the numerical simulation. The main component of the lift occurs at f_S. We denote this component by a_1L. The main component of the drag occurs at twice this frequency; that is why we denote it by a_2D (the numeral “2” in the subscript means twice the frequency of the fundamental frequency component of the lift). As a consequence of this two-to-one frequency relationship, the projection of the limit cycle of the wake in the lift-drag plane consists of two loops, as shown in Figure 5. This two-loop pattern was also observed at a higher Reynolds number of 500 [139,140].

Based on this two-to-one frequency relationship, proposed drag models for the fixed cylinder use quadratic coupling terms with the lift [141]. The first proposed model for this category was

$$C_D\hspace{0.33em}=\hspace{0.33em}{C}_{D,\hspace{0.33em}\mathrm{mean}}\hspace{0.33em}-\hspace{0.33em}2\hspace{0.33em}{\displaystyle \frac{a_{2D}}{\hspace{0.33em}{\omega }_S\hspace{0.17em}{a}_{1L}^2\hspace{0.33em}}}\hspace{0.33em}{C}_L\hspace{0.17em}{\dot{C}}_L$$

(3)

The rationale for the above ROM drag model is that it has the least number of necessary terms, which is two. So, it was made as simple as possible as an initial attempt to model the lift-drag coupling. The above ROM is based on the idea of splitting the drag coefficient signal into two independent components. First, there is a constant term ($C_{D,\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}})$, which is the mean (the time-averaged value or the DC component) of the drag coefficient signal. It is represented by the first term in the above ROM equation. Then, there is the oscillatory component (AC component) of the drag coefficient signal, which has a zero mean. It is represented by the second term in the above ROM equation. Given that the lift and drag are Cartesian force components of the same hydrodynamic force exerted on the cylinder surface, it was considered reasonable in the literature to relate the oscillatory component of the drag to the lift, which is a pure oscillatory signal (has zero mean). In addition, because a two-to-one frequency relationship was observed in the literature between the oscillatory component of the drag and the oscillatory lift, it was decided to model the oscillatory component of the drag as a quadratic function of the oscillatory lift signal. This quadratic dependence can take six forms, namely ($C_L^2$), ($-C_L^2$), (${\dot{C_L}}^2$), ($-{\dot{C_L}}^2$), ($\dot{C_L}{C}_L$), and ($-\dot{C_L}{C}_L$). Because the phase angle between the oscillatory component of the drag and the oscillatory lift was found in the literature to be near 270°, the last option was considered the most suitable. The multiplier factor of $2{\displaystyle \frac{a_{2D}}{{\omega }_S{a}_{1L}^2}}$ in the second term of the above two-term ROM drag model is added for proper scaling, such that the reconstructed time-domain drag coefficient signal from the ROM equation is close in magnitude to the observed training data for a fixed cylinder.

While the above drag model is able to reproduce the simulated C_D in the time domain to a good extent, it cannot reproduce the phase angle between the drag and the lift (i.e., the phase with which a_2D leads a_1L). Therefore, the ROM represented by Eq. 3 was extended to another amended ROM that has two quadratic terms, as shown in Eq. 4, instead of one term. The weights q₁ and q₂ are related to the phase angle found in the training dataset. This amended ROM is

$$C_D\hspace{0.33em}=\hspace{0.33em}{C}_{D,\hspace{0.33em}\mathrm{mean}}\hspace{0.33em}+2q_1{\displaystyle \frac{a_{2D}}{\hspace{0.33em}\hspace{0.17em}{a}_{1L}^2\hspace{0.33em}}}\left(C_L^2\hspace{0.17em}-\hspace{0.17em}{\left[C_L^2\right]}_{\hspace{0.33em}\mathrm{mean}}\right)+\hspace{0.33em}2q_2\hspace{0.33em}{\displaystyle \frac{a_{2D}}{\hspace{0.33em}{\omega }_S\hspace{0.17em}\hspace{0.17em}{a}_{1L}^2\hspace{0.33em}}}\hspace{0.33em}{C}_L\hspace{0.17em}{\dot{C}}_L$$

(4)

In the current study, we aim to extend this amended ROM further such that it holds valid for fixed and moving cylinders with a more general lift-drag relationship than the commonly assumed quadratic one. In the latter case, the two-to-one frequency relationship is distorted. To identify the behavior of the appropriate coupling in these new situations, we analyzed the simulated lift and drag coefficients for a vibrating cylinder at a vibration amplitude of 7.5% D. This value is relatively small for typical problems of VIV (vortex-induced vibration), but we selected it after conducting an extensive amount of CFD simulations at lower and higher amplitudes. Our amplitude selection criterion is the preference for an amplitude value that manifests a wide variety of lift and drag patterns as the motion inclination angle of the cylinder is varied while keeping the amplitude of this motion unchanged. This is important to achieve the aims of the current study and to provide challenging profiles of these forces (i.e., the lift and drag) and their relationship. Lower amplitudes can result in force patterns with a two-to-one relationship (or close). Existing drag models for this type already perform satisfactorily, and we then become unable to demonstrate the usefulness of the proposed new model. On the other hand, higher oscillation amplitudes may result in the dominance of chaotic or quasi-periodic force patterns [142,143,144]. In these cases, the limit cycle is highly distorted, and we then become unable to demonstrate the relationship between the model parameters and the temporal and spectral properties of the lift and drag. The reader is reminded that this study is not concerned primarily with the CFD simulations and numerically resolving the flow field or the lift and drag forces, although it is a large component of the current work. Rather, the CFD simulation conducted here shows how existing reduced-order drag models need to be extended to accommodate cases when the lift and drag relationship does not follow a commonly perceived type for which these models were derived. Also, the CFD simulation provides training datasets for computing the proposed reduced-order model (ROM) for the drag coefficient. We later show several analysis techniques that lead to the proper modifications of these existing reduced-order drag models, and emphasize the influence of the proposed modifications on the modeled results. The proposed reduced-order drag model here contains some parameters that can then be tuned (based on closed-form rigorous analytical expressions, not by trial and error) to other training datasets generated at higher motion amplitudes than the selected 7.5% D. The Reynolds number is maintained at 300, which is adequate here.

We first look at the projection of the limit cycle onto the C_L-C_D plane and the distortion in the two symmetric loops about the C_L = 0 axis. In Figure 6, the motion axis is rotating toward the cross-stream direction as we proceed from case (a) to case (f). The inclination angles of the oscillation line are from 40° in case (a), 50° in case (b), 60° in case (c), 70° in case (d), 80° in case (e), and 90° (thus cross-stream oscillation) in case (f). We show the change of the lift-drag relationship in the projection of the limit cycle onto the plane of their nondimensional coefficients (C_L and C_D) through this figure. The two-to-one frequency relationship (or the two-loop curve) that existed for the fixed-cylinder case is altered for all cases except for the case of cross-stream motion (the traditional case of one-degree-of-freedom motion). As the motion axis tilts away from the cross-stream direction, the lower loop (negative lift) becomes thinner while the upper loop (positive lift) becomes thicker. Eventually, the lower loop disappears, and we end up with one loop only, as in case (a) within Figure 6. This fact suggests adding a linear coupling of the drag to the lift. In case the two-loop curve is transformed into a multi-loop curve, then a better alternative is to propose a higher-order coupling (e.g., cubic) rather than a lower-order coupling (i.e., linear). This simple tool of graphical analysis for this nonlinear dynamical system revealed a lot of useful details about the dynamics of the problem, and allowed qualitative determination of the appropriate modeling type. However, this elementary analysis should be complemented with a more thorough and quantitative one. For this purpose, we analyze the same lift and drag signals (time series) in the spectral domain and identify the frequency components of each signal. The power spectra of the lift and drag for the cases presented in Figure 6 are shown in Figure 7. Again, case (f) corresponds to a cross-stream oscillatory motion. The spectra of the lift and drag with such a cross-stream motion are qualitatively similar to those obtained for the fixed cylinder (was shown in Figure 4). This confirms what we mentioned earlier about the similarity of the two cases when examined in the lift-drag plane. So, we were able to expose new lift-drag relationship patterns that were not identified before. With a cross-stream motion, the lift is composed of a dominant fundamental lift frequency component (a_1L) and its odd harmonics (such as a_3L), whereas the drag is composed of a mean quantity (C_{D, mean}) plus a dominant drag frequency component (a_2D) and its even harmonics (such as a_4D). So, a_2L and a_1D in this case are very small and can be neglected. As the motion axis tilts away from the cross-stream direction, a_1L increases slowly, and a_3L and a_2D show weak non-monotonic variations. However, a_2L and a_1D increase rapidly. Since a_2L remains about one order of magnitude less than a_1L, the distortion of the periodic lift signal is small. In contrast, a_1D approaches and even exceeds a_2D at some point, which is reflected in the strong deviation in the drag signal from the periodic pattern obtained either with a cross-stream motion or without any motion (a fixed cylinder). Based on the aforementioned findings and analysis, we propose the following extended (we refer to it as “universal”) reduced-order drag model with five terms:

$$\begin{array}{l}{C}_D\hspace{0.33em}=\hspace{0.33em}{C}_{D,\hspace{0.33em}\mathrm{mean}}\hspace{0.33em}+{\lambda }_1{\displaystyle \frac{a_{2D}}{\hspace{0.33em}\hspace{0.17em}{a}_{1L}\hspace{0.33em}}}{C}_L+\hspace{0.33em}{\lambda }_2\hspace{0.33em}{\displaystyle \frac{a_{2D}}{\hspace{0.33em}\Omega \hspace{0.17em}\hspace{0.17em}{a}_{1L}\hspace{0.33em}}}\hspace{0.33em}\hspace{0.17em}{\dot{C}}_L\\ +2q_1{\displaystyle \frac{a_{2D}}{\hspace{0.33em}\hspace{0.17em}{a}_{1L}^2\hspace{0.33em}}}\left(C_L^2\hspace{0.17em}-\hspace{0.17em}{\left[C_L^2\right]}_{\hspace{0.33em}\mathrm{mean}}\right)+\hspace{0.33em}2q_2\hspace{0.33em}{\displaystyle \frac{a_{2D}}{\hspace{0.33em}\Omega \hspace{0.17em}\hspace{0.17em}{a}_{1L}^2\hspace{0.33em}}}\hspace{0.33em}{C}_L\hspace{0.17em}{\dot{C}}_L\end{array}$$

(5)

The quadratic terms having the model coefficients q₁ and q₂ are almost identical to those presented earlier in Eq. 4 for the amended drag ROM, but the Strouhal frequency ω_S in Eq. 4 is replaced here by the sinusoidal vibration frequency Ω (both are equal in the current work). The new proposed linear-terms coefficients λ₁ and λ₂ are determined such that the phase of a_1D relative to a_1L (we denote this phase angle by ψ[a_1D, a_1L]) is matched with the training CFD-based value. So

$${\lambda }_1=\cos \left(\psi \left[a_{1D}/a_{1L}\right]\right); {\lambda }_2=\mathrm{s}\mathrm{i}\mathrm{n}\left(\psi \right[a_{1D}/a_{1L}\left]\right)$$

(6)

Similarly, the existing coefficients q₁ and q₂ for the quadratic terms are determined such that the phase of a_2D relative to a_1L (we denote this phase angle by ψ[a_2D, a_1L]) is matched with the training CFD-based value. So

$$q_1=\cos \left(\psi \left[a_{2D}/a_{1L}\right]\right); q_2=\mathrm{s}\mathrm{i}\mathrm{n}\left(\psi \right[a_{2D}/a_{1L}\left]\right)$$

(7)

In the old reduced-order drag models of Eq. 3 (the primitive two-term drag ROM) and Eq. 4 (the amended three-term drag ROM), the angular shedding frequency ω_S appears because for a fixed cylinder [145] (which was the scope of these models), the frequency of the main lift component takes place at this value. However, in the proposed five-term ROM of Eq. 5, we replaced this angular shedding frequency ω_S with the angular vibration frequency Ω, which can be different from ω_S, although Ω is set to be equal to ω_S here, reflecting a locking phenomenon [146,147]. This is to account for the fact that the shedding frequency (and the frequency of the main lift component) when the cylinder vibrates corresponds to the vibration frequency rather than to the Strouhal frequency that pertains to the fixed cylinder, when the wake is not excited by the mechanical vibration.

4. Results

We solve the flow field around a vibrating cylinder by direct numerical simulation (DNS) [148,149,150] of the two-dimensional Navier-Stokes equations at Re 300, which is low enough to justify the two-dimensionality and the use of DNS without turbulence modeling [151,152]. It should be noted that the previous statement does not imply a connection between the dimensionality of the Navier-Stokes equations and the Reynolds number. For example, turbulence (at high Reynolds number) [153] can be investigated using the two-dimensional Navier-Stokes equations [154,155]. On the other hand, laminar flow dynamics (at low Reynolds numbers) can be studied using the three-dimensional Navier-Stokes equations. Instead, we point out that in a real flow regime at Re 300, the three-dimensionality and turbulence exist, but not at a predominant scale [156,157]. The transition from a two-dimensional to a three-dimensional flow field, and from a laminar to a turbulent flow field, happens gradually as the Reynolds number increases. The selected Reynolds number of 300 is viewed as a suitable trade-off between (1) much lower values in which the flow is trivially steady without interesting dynamics, and (2) much larger values in which the complexities of the real flow are partly missed due to the inability of the two-dimensional laminar version of the Navier-Stokes equations to capture them.

The time-dependent lift and drag coefficients are recorded as time signals during each case of computational fluid dynamics (CFD) simulation, and these time signals are obtained after a stable periodic time response is reached.

Spectral analysis is then applied to these signals to evaluate the significant components and their lead angles (relative to a1L). These evaluated spectral parameters of the simulated lift coefficient are used in Eqs. 5~6 to generate a fitted drag-coefficient signal (i.e., obtained by applying the proposed five-term reduced-order drag model after identifying the numerical values of its five parameters through closed-form analytical expressions [158]), which is then compared with the one obtained from the CFD simulation (the training set). We conducted several cases and got a good agreement in them. Some of these cases are discussed next.

In Figure 8, we compare the two drag-coefficient signals (the CFD-simulated and the reduced-order modeled) for the fixed-cylinder case. In this case, a1L = 0.905, a1D = 0.0016, a2D = 0.077, ψ[a1D, a1L] = 256°, and ψ[a2D, a1L] = 335.2°. Consequently, λ1 = –0.241 and λ2 = –0.970; while q1 = 0.908 and q2 = –0.419. The mean drag coefficient (its DC or zero-frequency component) based on the training CFD signal is 1.376. As shown in the figure, the agreement for the drag coefficient is excellent. We would like to add here that the superharmonics in the lift coefficient beyond a3L are negligible, which allows replacing the lift coefficient signal obtained from the CFD simulation with a third-order approximation of the form:

$$\begin{array}{l}{C}_L\hspace{0.33em}\cong \hspace{0.33em}{a}_{1L}\hspace{0.33em}\cos \left(\Omega \hspace{0.17em}t\right)+a_{2L}\hspace{0.33em}\cos \left(2\hspace{0.17em}\Omega \hspace{0.17em}\hspace{0.17em}t\hspace{0.33em}+\hspace{0.33em}\psi \left[a_{2L}/a_{1L}\right]\right)\\ + a_{3L}\hspace{0.33em}\cos \left(3\hspace{0.17em}\Omega \hspace{0.17em}\hspace{0.17em}t\hspace{0.33em}+\hspace{0.33em}\psi \left[a_{3L}/a_{1L}\right]\right)\end{array}$$

(8)

The definitions of ψ[a_2L, a_1L] and ψ[a_3L, a_1L] are similar to the definition of ψ[a_1D, a_1L]. In Figure 9, we compare this harmonic approximation of the lift coefficient to the simulation results, and also the modeled drag coefficient that corresponds to either of the two lift signals. The difference is very small, in fact, which is again due to the insignificance of the neglected superharmonic components in the lift coefficient. Thus, Eq. 8 can be viewed as an algebraic model for the lift coefficient.

In Figure 10, we compare the two drag-coefficient signals for a case where a_1D is significant and even larger than a_2D. This case corresponds to case (b) in Figure 6~7 (the tilt angle of the mechanical excitation line is 50° from the horizontal free stream), a_1L = 0.902, a_1D = 0.16, a_2D = 0.123, ψ[a_1D, a_1L] = 223°, and ψ[a_2D, a_1L] = 352.9°. Consequently, λ₁ = –0.736 and λ₂ = –0.677; while q₁ = 0.992 and q₂ = –0.124. The mean drag coefficient (its DC or zero-frequency component) based on the training CFD signal is 1.445. In the same figure, we also show the C_L signal (as obtained from the CFD simulation). Again, the fact that the main harmonic of the lift is more than one order of magnitude larger than any of its higher harmonics (superharmonics) minimizes the modulation in the C_L signal and makes it close to a pure-harmonic signal. In the same figure, we demonstrate the significance of the added linear-coupling terms by comparing the simulated drag-coefficient signal with the modeled one, but after dropping the linear terms. The disagreement is remarkable, which is a result of assuming a two-to-one frequency relationship between the lift and drag in this case, and such an assumption is not true in this case.

As the motion axis gets closer to the cross-stream, the amplitude a_1D decreases, whereas a_2D exhibits slight variations. For case (d) in Figure 6~7 (the tilt angle of the mechanical excitation line is 70° from the horizontal free stream), a_1D = 0.0796 and a_2D = 0.124. Although a_2D is larger now than a_1D, the contribution of the latter is still important. In Figure 11, we indicate this by comparing the modeled drag-coefficient signal with and without the linear-coupling terms (i.e., with and without accounting for a_1D) to the simulated signal. Again, the proposed universal drag model performs very well, whereas the original model (without the extra added terms) shows some level of discrepancy that is less than the one in Figure 10 (due to decreased ratio of a_1D/a_2D), but it is still considerable.

The final case is for the cross-stream motion, which is case (f) in Figure 6~7. This case does not represent a challenge to the proposed model since the two-to-one frequency relationship between the drag and lift is restored (as was originally the case for a fixed cylinder, before introducing any motion). This explains the big similarity between the modeled drag signal with and without the linear terms, as indicated in Figure 12.

We would like to add that after running the simulations for several cases of interest (i.e., at different tilt angles of the cylinder’s vibration axis), one can fit the various needed spectral quantities for the proposed model, such as a_1L and a_1D, into suitable functions of the inclination angle, which can be used later to produce the modeled drag signal at any arbitrary inclination angle. This saves time remarkably, knowing that generating the modeled signal takes about three orders of magnitude less time than generating the CFD-based signal. This also allows for analyzing the tilt sensitivity [159,160] of the flow parameters. Similarly, this concept of one-time construction of a lookup table (LUT) or database [161,162] of model parameters can be repeated at different Reynolds numbers of interest or different oscillation amplitudes. After building such a database of parameters, the reduced modeling can be used to conveniently generate interpolated signals within a continuous parameter space.

We also would like to comment on the applicability of the proposed extended drag model in more complex three-dimensional flow regimes (such as in the case of a finite-length cylinder with a high Reynolds number of 10⁶, for example) and whether the two added nonlinear terms in the drag coefficient model are able to explain the drag-lift coupling in such three-dimensional flow dynamics. Due to the three-dimensional nature of such a situation, the signals of either the lift coefficient or the drag coefficient in a real problem become less regular than in the case of simpler two-dimensional flow. However, the three-dimensional flow still shows a coherent vortex street with distinguished frequency components in the lift coefficient [163,164], in a similar manner to the two-dimensional flow problem [165,166]. Therefore, the proposed drag model remains useful in three-dimensional problems.

5. Discussion

Although the current study focuses on reduced-order modeling of the drag coefficient, it can be useful to briefly discuss the reduced-order oscillator model for the lift coefficient also. The current section is devoted to achieving this aim.

Despite the availability of more than one form for modeling the lift coefficient, we consider here one simple form (van der Pol oscillator) [167,168,169] for a fixed cylinder. Such a simple oscillator facilitates the discussion and avoids the extensive details that are involved with a more complicated modeling approach.

A van der Pol oscillator [170,171,172] for the lift coefficient ($C_L$) has the following form:

$$\ddot{C_L}+{\omega }_S C_L^2-\mu \dot{C_L}+\alpha C_L^2 \dot{C_L}=0$$

(9)

where $\ddot{C_L}$ is the second time derivative of the lift coefficient, ${\omega }_S$ is the Strouhal angular frequency, $\mu$ is a coefficient for linear destabilizing damping [173,174], $\dot{C_L}$ is the first time derivative of the lift coefficient, and $\alpha$ is a coefficient for nonlinear (cubic) stabilizing damping [175,176,177].

At a low value of $C_L$ (and thus $C_L^2$), the nonlinear damping term becomes weaker than the linear destabilizing damping term. Therefore, the lift coefficient tends to be amplified. On the other hand, a high value of $C_L$ (and thus $C_L^2$) causes the nonlinear damping term to be stronger than the linear destabilizing damping term. Therefore, the lift coefficient tends to decay. The opposing influences of the two damping terms lead to a balanced response that is a periodic limit cycle [178,179,180], where the time signal of the lift coefficient oscillates with a fixed amplitude.

6. Conclusions

We considered non-traditional wakes formed as a result of applied mechanical oscillatory perturbations, where the wake of a fixed cylinder is excited by introducing a prescribed vibrational motion along an inclined straight line that lies between the cross-stream direction and the stream-wise direction. This revealed new types of the relationship between the induced lift and drag forces acting on the cylinder. The typical quadratic lift-drag coupling alone fails to reproduce this relationship. We first analyzed different lift and drag signals obtained from computational fluid dynamics (CFD) simulations. The analysis techniques were in the time domain (time series analysis), the frequency domain (spectral analysis), and the lift-drag plane (limit cycle projection). Our observation that the two-loop curve of the lift-drag projection of the limit cycle gradually transforms into a single-loop curve suggested that additional linear lift-coupling terms in the drag reduced-order model (ROM) are needed, especially when the previously neglected component of the drag coefficient at the frequency of the main lift component (a_1D) increases. We added two linear terms to an existing three-term drag ROM so that it is extended and can handle non-traditional wake cases where the two-to-one frequency relationship between the drag and lift is broken. We related the model parameters to the spectral variables of the lift and drag signals. We performed several tests of the proposed five-term universal drag model, including the fixed cylinder case and the cross-stream vibration case. These two special cases exhibit a two-to-one frequency relationship between the drag and lift. We also performed tests of the proposed five-term algebraic drag model for cases where the drag has significant components at the motion frequency in addition to the typical component at twice this frequency, so both linear and quadratic couplings become important. All the comparisons showed excellent agreement between the proposed algebraic reduced-order drag model and the time-consuming CFD simulations, which require solving the Navier-Stokes equations that govern the flow field.