Vortex-induced vibrations of a cylinder in inelastic shear-thinning and shear-thickening fluids

Abstract Vortex-induced vibrations (VIV) of a cylinder in a Newtonian fluid is a model problem in fluid–structure interactions and has been studied extensively. In this work, we study the influence of shear-thinning and shear-thickening fluids on the VIV response of a one-degree-of-freedom flexibly-mounted cylinder. We consider a system with a mass ratio of $m^*=2$ and zero structural damping in shear-thinning and shear-thickening power-law fluids at $Re_0 = 15$ and $Re_0 = 200$, respectively, defined based on the zero-shear-rate viscosity of the fluids. We investigate how the VIV amplitude and frequency, flow forces, and the vorticity contours change as the reduced velocity, $U^*$, and fluid's time constant, $\lambda$, change. When the results are compared based on $Re_0$, shear-thinning fluids enhance the oscillations while shear-thickening fluids suppress them. If, however, we define a characteristic Reynolds number, $Re_{char}$, based on a viscosity evaluated at the characteristic shear rate, $\dot {\gamma } = U/D$, then at a constant $Re_{char}$, the amplitude of response stays very similar for the shear-thinning, shear-thickening and Newtonian fluids. Despite this similarity, the observed far wake is different: shear thinning amplifies the generation of vorticity and reduces the extent of the wake, whereas shear thickening limits the generation of vorticity and extends the wake. Our findings show that the local apparent viscosity observed close to the cylinder placed in shear-thinning or shear-thickening fluids governs the VIV response of the cylinder.


Introduction
When a flexible or flexibly-mounted bluff body is placed in flow, formation of vortices downstream of the bluff body generates unsteady forces on the body, which in turn cause the structure to oscillate.When the structure oscillates, the shedding frequency and oscillation frequency are synchronized.This synchronization between the shedding frequency and oscillation frequency is called lock-in, and the resulting oscillations are called vortex-induced vibration (VIV).VIV has been studied extensively for the cases where a structure is placed in a Newtonian fluid.When a one-degree-of-freedom (1DOF) cylinder is placed in Newtonian flow and it is free to oscillate in a direction perpendicular to the direction of incoming flow (crossflow (CF) VIV), lock-in is observed for a range of reduced velocities, U * (defined as U * = U/f n D, in which U is the velocity of the incoming flow, D is the cylinder's diameter, and f n is the natural frequency of the system in vacuum), with amplitudes of oscillations of up to around one cylinder diameter (Sarpkaya 2004;Williamson & Govardhan 2004).In the lock-in region, a 2S shedding pattern (in which two single vortices are shed in the wake of the cylinder in each cycle of oscillations) or a 2P shedding pattern (in which two pairs of vortices are shed in the wake in each cycle of oscillations) is observed.If the 1DOF cylinder is free to oscillate in the direction of flow (inline (IL) VIV) then oscillations are observed over two ranges of U * values, with amplitudes of around 0.1D.The oscillations in the first range of reduced velocities are with a symmetric shedding of vortices in the wake, and in the second range with an asymmetric wake (Cagney & Balabani 2013a,b;Gurian, Currier & Modarres-Sadeghi 2019).When the cylinder is free to oscillate in both the CF and IL directions (2DOF VIV), lock-in occurs in both directions, and figure-eight trajectories are observed in the response of the cylinder together with several different types of vortex shedding patterns, including 2T shedding (in which two triplets of vortices are shed during each cycle of oscillations) (Dahl, Hover & Triantafyllou 2007;Dahl et al. 2010).VIV has also been studied for non-circular cross-sections, such as a square prism or a triangular prism, in which cases, besides VIV, galloping response has been observed, due to the non-zero mean lift forces that act on the structure (e.g.Nemes et al. 2012;Zhao et al. 2014;Seyed-Aghazadeh, Carlson & Modarres-Sadeghi 2017;Carlson, Currier & Modarres-Sadeghi 2021).Studies on VIV have been extended to cases of flexible continuous structures placed in flow, due to the fact that in real-life applications, VIV is observed in long offshore structures such as risers in oil platforms or mooring lines of floating wind turbines.When a flexible structure undergoes VIV, synchronization can be observed between shedding of vortices and several different modes of the structure at different locations along its length, resulting in monoor multi-modal oscillations of the structure (e.g.Vandiver 1993;Bourguet et al. 2011;Wu, Ge & Hong 2012;Seyed-Aghazadeh, Edraki & Modarres-Sadeghi 2019).
When a fixed cylinder is placed in a shear-thinning or shear-thickening fluid (i.e. a non-Newtonian fluid for which the viscosity, η, varies with shear rate, γ ), there are several differences between the flow behaviour in its wake and that of a cylinder placed in a Newtonian fluid, which result in differences in forces that act on the cylinder in these two cases.The critical Reynolds number (where, Re = ρUD/η, in which ρ is the density of the fluid, U is the incoming flow velocity, D is the cylinder diameter, and η is the dynamic viscosity of the fluid) to observe vortex shedding, which is known to be Re crit = 47 for Newtonian fluids (Mathis, Provansal & Boyer 1984), decreases with shear-thinning effects and increases with shear-thickening effects (Lashgari et al. 2012;Şahin and Atalık 2019).Critical Reynolds numbers as low as Re crit = 3 for shear-thinning fluids and as high as Re crit = 193 for shear-thickening fluids have been reported (Lashgari et al. 2012), where Re is defined based on the zero-shear-rate viscosity, η 0 .In the range of Reynolds numbers where the shedding is observed, the vortex shedding frequency increases with shear-thinning effects and decreases with shear-thickening effects (Bailoor, Seo & Mittal 2019;Şahin and Atalık 2019).Shear-thinning effects reduce the formation length in the wake of a fixed cylinder, and shear-thickening effects increase the formation length (Coelho & Pinho 2003a,b, 2004;Lashgari et al. 2012;Bailoor et al. 2019;Şahin and Atalık 2019).Shear-thinning effects reduce the drag coefficient (at least for low Reynolds numbers Re ≤ 300) (Lashgari et al. 2012;Bailoor et al. 2019;Alam et al. 2021) and intensify the magnitude of the vorticity in the region close to the cylinder, due to the reduced shear stress associated with the shear-thinning effects that occur very close to the cylinder (Lashgari et al. 2012).In these studies on shear-thinning and shear-thickening fluid, Şahin and Atalık (2019) use a power-law fluid model in which viscosity, η, is proportional to the shear rate, γ , to the power n − 1, i.e. η = m γ n−1 , where m is the consistency index, and n is the power-law coefficient.For Newtonian fluids, n = 1 and m = η.In other studies mentioned here, the Carreau model is used in which the power-law regime is preceded by a plateau at the zero-shear-rate viscosity and followed by a plateau at the infinite-shear-rate viscosity.
The wake of a cylinder that is forced to oscillate in a shear-thinning fluid has been investigated very recently as well.Hopkins & de Bruyn (2019) showed that the effective viscosity, determined by averaging the viscosity around the circumference of the cylinder and over one period of oscillations, depends significantly on the driving frequency, and it reaches the minimum at a resonant frequency.Alam et al. (2021) studied the vortex dynamics in the locked-in mode and non-locked-in mode of a transversely oscillating cylinder in shear-thinning fluids at Re = 100, where Re is defined based on the zero-shear-rate viscosity.They showed that the vortex separation length is shorter in the shear-thinning fluids compared with the Newtonian fluid due to reduced viscous diffusion near the vortices in the wake.This length decreases with smaller power-law coefficient, n, and larger Carreau numbers, Cu = λU/D, where λ is a time constant of the fluid.However, the effect of n is smaller at higher Cu.In shear-thinning fluids, the 2S vortex shedding mode is observed during locked-in forced oscillations.The paired counter-rotating vortices (P+S and 2P vortex shedding modes) are observed during non-locked-in forced oscillations.This finding implies that in a self-excited VIV case, only 2S shedding would be observed in the wake.Alam et al. (2021) also show that the vortices in a shear-thinning fluid have a stronger asymmetric pattern with higher levels of unsteadiness while departing from the lock-in mode.
In the present work, we will study VIV of a 1DOF cylinder free to oscillate in the CF direction and placed in inelastic power-law fluids.The goal is to understand the shear-thinning and shear-thickening effects in the absence of elasticity on the self-excited response of a 1DOF cylinder undergoing VIV, as a model problem in fluid-structure interactions (FSI).We will discuss the effect of system parameters -such as the reduced velocity, the time constant and the power-law coefficient of the fluid -on the VIV response of the cylinder.We will also introduce a characteristic Reynolds number, which considers the local effects due to shear-rate-dependent viscosity, and can collapse the observed response of the cylinder placed in shear-thinning, shear-thickening or Newtonian fluid.

Governing equations and numerical methods
We consider two-dimensional, incompressible flow of inelastic shear-thinning and shear-thickening fluids in a domain containing a flexibly-mounted cylinder free to oscillate in the crossflow direction.The fluid flow is governed by the unsteady, incompressible Navier-Stokes (N-S) equations (2.1) 10 1 10 2 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 0.0035 0.006 where  = ηγ and γ = ∇u + ∇u T .Unlike in a Newtonian fluid, here, the viscosity is a function of shear rate.Shear-rate-dependent viscosity has been described using the Carreau model as (Morrison 2001) where η 0 is the zero-shear-rate viscosity, η ∞ is the infinite-shear-rate viscosity, n is the power-law coefficient, which describes the slope of increasing or decreasing viscosity curve, and λ is a time constant of the fluid.The value of λ determines the shear rate at which the transition occurs from zero-shear-rate plateau to power-law regime.With increasing time constant, transition occurs at lower shear rates.In the Carreau model, γ describes the local shear rate based on the second invariant II γ of the strain rate tensor as (2.4) Figure 1 shows the shear-rate-dependent viscosity for several time constants of shear-thinning and shear-thickening fluid used here.The characteristic shear rate, defined as the ratio of the incoming flow velocity and the cylinder diameter, γchar = U/D, is shown in the figure using a vertical line.For all the cases presented here, the maximum shear rate remains in a range such that the viscosity does not reach the infinite-shear-rate viscosity plateau.
The finite volume method is used to discretize the N-S equations.Unsteady N-S equations are solved using a coupled algorithm where the momentum equation and pressure-based continuity equations are solved together.A quadratic upwind interpolation for convective kinematics (QUICK) scheme has been applied to discretize the convective terms in the momentum equation.The least squares cell-based method is used to spatially discretize the gradients in the convection and diffusion terms.Pressure and velocity are stored at the cell centres.Since momentum equations require the value of pressure at the face between two adjacent cells of the unstructured grid, the PRESTO (PREssure STaggering Option) scheme has been used to interpolate pressure at the face.The PRESTO scheme uses discrete continuity balance for a staggered control volume about the face to compute pressure at the face.A second-order implicit scheme has been used for temporal discretization of the transient derivative terms.Convergence tolerances for continuity and both velocity components are set to 10 −6 .Once the fluid equations are solved, the net force acting on the cylinder in the y-direction (the crossflow direction) is used to calculate the displacement of the cylinder using its equation of motion.
The equation of motion for a flexibly-mounted cylinder is obtained using Newton's second law of motion as mÿ + cẏ + ky = 1 2 ρC y DU 2 , (2.5) where in which y is the position of the centre of the cylinder, m is its mass, c is a damping coefficient, k is the spring constant, ρ is the density of the fluid, C y is the force coefficient acting on the cylinder in the y-direction, n is the normal unit vector, and j is the unit vector in the y-direction.
Using the diameter, D, as a length scale and the incoming flow velocity, U, as a velocity scale, the equation of motion (2.5) can be written in a non-dimensional form as where are the mass ratio and the reduced velocity, respectively, where (2.10) In (2.7), we have assumed zero structural damping to promote maximum amplitude of oscillations.The second-order ordinary differential equations are converted into two first-order equations and solved using an iterative solver.Following the velocity of the cylinder, the mesh nodes are moved using the diffusion-based smoothing method, in which we solve the modified Laplace equation (2.11) where u is the point velocity field used to modify the position of mesh nodes, x old and x new are the point positions before and after the mesh motion, respectively, and t is the time step.In the modified Laplace equation, β is a constant or variable diffusion field, chosen to govern the mesh motion.We have defined β based on the boundary distance, as β = 1/l a , where l is the distance of the cell centre from the selected boundary, and parameter a describes how the cylinder's motion diffuses through the surrounding mesh: a = 0 934 A39-5 indicates uniform diffusion, and a = 1 and a = 2 indicate linear and quadratic diffusion, respectively, where mesh nodes close to the cylinder move more than mesh nodes far from the cylinder.In the present work, we have used uniform diffusion.Re-meshing would occur only if the quality of mesh anywhere in the domain deteriorates below provided tolerance values during the smoothing step.Parameters used in the simulations are shown in table 1.
Figure 2 shows the 30D × 16D two-dimensional domain that is meshed using a structured grid.The total number of grid nodes is 40 968.The mesh in the domain comprises two zones: the inner zone and the outer zone.The mesh around the cylinder in the inner zone moves along with it as a rigid body, thus maintaining the quality of the mesh near the cylinder.This arrangement facilitates the deformation of the mesh elements lying in the outer zone due to the cylinder movement.The slip boundary condition (zero velocity gradient) is applied at the top and bottom walls.The flow is uniform and steady at the inlet.At the fluid-solid interface, no slip and no mass flux conditions are applied.The reduced velocity is varied by changing the natural frequency of the system while keeping the mass ratio and the incoming flow velocity constant.This approach of changing the reduced velocity is chosen to keep the Reynolds number Re 0 , defined based on the zero-shear-rate viscosity (Re 0 = ρUD/η 0 ), constant.The value of the Reynolds number has been chosen such that the maximum local Reynolds number in the domain stays within the laminar flow regime.The simulations are run for approximately 100 oscillation cycles, depending on the onset of the steady state.When the transient is passed, the steady-state results are collected for at least 20 oscillation cycles and used for analysis.

Verification
We have compared our results with published data for the case of VIV of a Newtonian fluid.This comparison is summarized in figure 3, where the dimensionless amplitude of the cylinder response, A * = A/D, is plotted versus the reduced velocity, U * .In all cases shown in this plot, the Reynolds number is Re = 150, the mass ratio is m * = 2, and the structural damping is ζ = 0.In all cases, the amplitude of oscillations increases initially, as the reduced velocity is increased, and then with further increase in the reduced velocity, the amplitude starts decreasing monotonically until it goes to zero.The maximum amplitude of oscillations in the present case and in the results by Borazjani & Sotiropoulos (2009) is A * = 0.54, and in the results by Ahn & Kallinderis (2006) is only slightly larger, A * = 0.56.The onset and the width of the lock-in range are also in agreement among all three sets of results.Toward the end of the lock-in range, the amplitudes of the present results are slightly smaller than those from the other two cases.For all cases, the lock-in range extends from U * = 3 to U * = 8, and as soon as the structure starts to oscillate, the shedding frequency follows the oscillation frequency, instead of following the Strouhal law, indicating that lock-in has occurred.Overall, this comparison shows that our results  are in agreement with the results of previous studies on predicting the onset and the width of the lock-in range, as well as the amplitude of oscillations.

Response of a 1DOF cylinder in the flow of shear-thinning fluid
The rheology of shear-thinning fluid used here has been described using a Carreau model as shown in figure 1(a).We have used a zero-shear-rate viscosity of η 0 = 0.056 Pa s, an infinite-shear-viscosity of η ∞ = 0.0035 Pa s and a power-law coefficient of n = 0.36, and we have varied the time constant from λ = 0.1 s to λ = 5 s (Cu ≈ 1 to Cu ≈ 40).In all the results presented in this section, unless explicitly mentioned, a non-zero displacement (0.6D) and a zero velocity are given as the initial conditions of the cylinder.Based on these parameters, the Reynolds number at the inlet of the domain is Re 0 = 15.For a flexibly-mounted cylinder placed in a Newtonian fluid at this Reynolds number, VIV is not expected.For Newtonian fluids, the critical Reynolds number to observe shedding of vortices in the wake of a cylinder is Re crit = 47 (Mathis et al. 1984;Jackson 1987;Dušek, Gal & Fraunié 1994) 4 shows the dimensionless amplitude, A * , and frequency, f * , of the cylinder's displacements (figure 4a,b) and the transverse (CF) force coefficient, C y , and frequency, f * C y , that act on the cylinder (figure 4c,d) as a function of reduced velocity, U * , for time constants varying from λ = 0.15 s to λ = 5 s (Cu = 1.2 to Cu = 40).The amplitude of oscillations is normalized by the cylinder diameter, D, and the frequency is normalized by the natural frequency of the system in vacuum, f n .
For the fluid with the largest time constant presented here, λ = 5 s (Cu = 40), the amplitude response of the cylinder in figure 4(a) resembles a typical VIV response for a Newtonian fluid, albeit at an incoming Reynolds number lower than the minimum Reynolds number for which VIV can be observed for a Newtonian fluid.As the time constant is decreased, however, the amplitude of oscillations, as well as the width of the lock-in range, decreases.At λ = 0.15 s (Cu = 1.2), no oscillation is observed, and VIV is completely suppressed for all reduced velocities.For all cases where oscillations are observed, the oscillation frequency stays constant at a value close to f * = f o /f n = 1 (figure 4b).In this range, the frequency of fluctuating force in the direction of oscillation, f * C y , stays close to f * = 1 as well (figure 4d), indicating that the shedding frequency and the oscillation frequency are synchronized, lock-in is observed, and the observed oscillations are indeed VIV.The magnitude of the CF force coefficients shown in figure 4(c) amplifies when the cylinder starts oscillating for all cases with different time constants.The largest magnitude of the CF force coefficient is observed for the largest time constant tested.
The onset of oscillations for λ = 5 s (Cu = 40) is at U * = 3 (figure 4a), and it is delayed with decreasing time constant.Consequently, the width of the lock-in range decreases in the case of shear-thinning fluid with decreasing time constant.The onset of VIV depends on the Strouhal number associated with the shear-thinning fluid.As shown by Bailoor et al. (2019), the Strouhal number (St = f s D/U) increases for stronger shear-thinning effects (i.e.larger time constant, λ, or smaller power-law coefficient, n).This suggests that as the time constant is decreased in the cases discussed here, the Strouhal number decreases, and as a result the shedding frequency decreases for a constant incoming flow velocity.Then the shedding frequency matches the natural frequency of the system at higher

Distribution of local Reynolds numbers
As shown in figure 1, the time constant of a fluid determines the onset of transition from the zero-shear-rate viscosity plateau to the power-law regime.With increasing λ, shear thinning of the viscosity occurs at lower shear rates.As a result, the fluid around the cylinder has a viscosity that depends on both the local shear rate (defined in equation (2.4)) and the time constant of the fluid.This can be observed in figure 5, where the spatial distribution of the local Reynolds number, defined as Re = ρUD/η, where η is the local viscosity of the fluid calculated from the shear rate at that location, is presented for a constant reduced velocity, and four different values of λ: one corresponding to a case with no oscillations (figure 5a) and others corresponding to cases with oscillations (figure 5b-d).The contours of local Reynolds number shown in figure 5 indicate the variation in the viscosity depending on the shear rate in the domain.The time constant determines how far from the cylinder the shear-thinning effect extends.A smaller time constant corresponds to a shear-thinning effect only in close proximity to the cylinder.From the figure, we observe that the shear-thinning effect extends up to a considerable distance downstream from the cylinder for a larger time constant.The reach of the shear-thinning effect in the domain can be understood from the histograms of local Reynolds number for various time constant values.A fixed rectangular bounding box of size 15D × 7D has been created around the cylinder, and the local Reynolds number at each grid-cell inside the box is calculated.Histograms of the local Reynolds number inside the bounding box for four different time constant values are shown in figure 5.For λ = 0.4 s (Cu = 3), where no oscillation is observed, most of the flow is dominated with Reynolds numbers below Re = 20, because shear-thinning occurs only in the regions very close to the cylinder where the shear rate is maximum.As the time constant is increased and the shear rate needed for shear thinning to occur is decreased, the percentage of the flow at large local Reynolds numbers gradually increases until the flow is sufficiently dominated by inertial effects for vortices to separate and shed from the cylinder and drive VIV.We have used the percentage of grid-cells to plot the histograms of local Reynolds number.Since the grid is concentrated close to the cylinder, the results are biased towards the viscosity close to the cylinder.
Note that the incoming Reynolds number, Re 0 , which is defined based on zero-shear-rate viscosity, is the same for all cases shown in figure 5

The response at constant reduced velocities and for varying time constants
In the previous subsection, we discussed the response of a 1DOF system placed in a shear-thinning fluid as the reduced velocity was varied.In this subsection, we focus on the influence of fluid rheology on the system's response as other system parameters stay constant.To do so, we use the time constant, λ, as an independent variable and investigate how the response of the system changes for two constant reduced velocities of U * = 4 and U * = 6, as two sample cases.The Reynolds number has been kept constant at Re 0 = 15 for all these cases.Figure 6 shows the amplitude, A * , and frequency, f * , of response, the CF force coefficient, C y , the ratio of the third harmonic to the first harmonic of the CF forces, C y,3 /C y,1 , and the phase between the displacement and the CF force, φ, for the two sample reduced velocities as λ is varied from λ = 0 s to λ = 5 s (Cu = 0 to Cu = 40).
As shown in figure 6(a), the onset of VIV is at λ = 0.5 s (Cu ≈ 4, Re char ≈ 34) and λ = 0.2 s (Cu ≈ 2, Re char ≈ 22) for U * = 4 and U * = 6, respectively.The critical value for characteristic Reynolds number is different for different U * .The VIV amplitude increases rapidly before it reaches its peak at λ = 1.5 s (Cu ≈ 12, Re char ≈ 59) and λ = 0.8 s (Cu ≈ 6, Re char ≈ 43) for U * = 4 and U * = 6, respectively.By increasing the time constant further, the amplitude decreases slightly for both reduced velocities.The decrease in amplitude is more noticeable in the case of U * = 6.When the time constant of the shear-thinning fluid approaches the largest values tested (approximately at λ = 20 s, Cu = 158), most of the fluid in the wake of the cylinder is at the infinite-shear-rate viscosity and, as a result, the fluid essentially is no longer shear-thinning and the VIV response of the cylinder approaches Newtonian-like behaviour.This is manifested in the form of a plateau in the amplitude at higher time constants.Figure 6(b) shows the variation of frequency of oscillations with time constant for U * = 4 and U * = 6, where the frequency of oscillations is normalized by the natural frequency of the system in vacuum.The f * behaviour in the case of U * = 6 is very similar to the f * behaviour that is typically observed in a Newtonian case when f * is plotted versus U * .Initially, f * is smaller than1.Then, as λ is increased, f * crosses1.This point of crossing 1 corresponds to a switch in the phase between the displacement and force from φ = 0 • to φ = 180 • , and the appearance of a large contribution of the third harmonic force (at values of around λ = 2.4 s (Cu = 19) in the present case).The slight decrease in the amplitude of oscillations for larger values of λ in this case can also be explained by the sudden change in the phase difference, since the displacement of the cylinder and the CF force are out of phase when λ > 2.4 s (Cu > 19).
In the case of U * = 4, however, f * remains lower than 1 for all values of λ, and the large third harmonic component of the force and the sudden phase shift are not observed in the response.For this reduced velocity, the amplitude of oscillations stays constant for larger λ values, since the CF displacement and the CF force stay in phase.
Figure 7 shows the Lissajous curves for several different values of the time constant.As the time constant increases, the Lissajous curves rotate in the counterclockwise direction.For U * = 4 (figure 7a), the Lissajous curves are located in the first and third quadrants, corresponding to a phase difference of φ = 0 • , which means that the displacement and the force coefficient are positively correlated.For U * = 6 (figure 7b), the Lissajous curves move to the second and fourth quadrants for λ > 2.4 s (Cu > 19).This corresponds to the phase jump to φ = 180 • , which means that the displacement and the force coefficient are negatively correlated.The two lobes at the extreme ends in the Lissajous curves represent the third harmonic contribution of the force.It is clear from figure 7 that the contribution of the third harmonic increases for increasing time constants, when the Lissajous curves stay in the first and third quadrants.In figure 7 second quadrant, after which the forcing and displacement become out of phase, and the contribution of higher harmonics decreases.

Subcritical instability in shear-thinning fluids
The results presented in the previous subsections were all for cases where a non-zero initial displacement was given to the structure.Over a range of λ values, we found that oscillations are not observed if the structural initial conditions stay at zero. Figure 8(a) shows a sample amplitude plot for U * = 4 in which the amplitude of response is plotted versus λ for cases with both zero and non-zero initial conditions.A range is observed in the figure from λ = 0.5 s (Cu ≈ 4, Re char ≈ 34) to λ = 0.7 s (Cu ≈ 6, Re char ≈ 40) for which two stable solutions exist: a zero response, and a non-zero response.This plot exhibits a subcritical instability for the flexibly-mounted structure.Similarly to the subcritical VIV that have been observed in a 1DOF system placed in a Newtonian fluid (Mittal & Singh 2005;Boersma et al. 2021), if the cylinder is not given non-zero initial conditions in this range of λ values, then it remains at its initial equilibrium position and no vortices are shed in its wake (figure 8b).If, however, the cylinder is given an initial disturbance, then its wake becomes unstable, vortices are shed, and the cylinder undergoes VIV (figure 8c).The initial disturbance provides an additional shear rate due to the relative motion of the cylinder, which drives the Reynolds number up in regions close to the cylinder and pushes it past the critical Re that is needed to observe VIV.By analysing the histograms of local Reynolds numbers inside a bounding box of 15D × 7D for two cases where oscillations are observed -i.e.λ = 0.5 s (Cu ≈ 4) with an initial disturbance, and λ = 0.75 s (Cu ≈ 6) without any initial disturbance -we find that the distributions are quite similar, with comparable maximum Re values of 116 and 127, respectively.

Wake for shear-thinning fluids
In this subsection, we discuss the flow pattern in the wake of the cylinder as it undergoes VIV.The Reynolds number is kept constant at Re 0 = 15 for all cases.The vorticity 934 A39-14 contours are plotted for various time constants of the fluid in figure 9, where wake patterns are shown for increasing time constants from top to bottom for U * = 4 (a-e) and U * = 6 ( f -j).The snapshots in each row are at different time constants, but they are chosen such that in each row, the extent of the wake is similar.Vortex shedding is not observed for the first sample cases for each reduced velocity, and the cylinder does not oscillate.For all the other cases the cylinder oscillates, and vortex shedding is observed in the wake.For all these cases, a 2S shedding pattern is observed in the wake, in which one single vortex is shed from each side of the cylinder during each cycle of oscillations (Williamson & Roshko 1988).The vortex shedding frequency and the amplitude of oscillations for U * = 4 are higher than those for U * = 6.The lateral distance between counter-rotating vortices is larger for U * = 4 due to higher amplitude of oscillations.The length of recirculation bubble is larger for U * = 6, because the oscillation frequency is lower at this reduced velocity when compared with U * = 4, and the shear layers are cut at a slower rate by the cylinder and they can stretch themselves farther in the wake of the cylinder.The strength of generated vortices depends on the time constant of the fluid.For small time constants, the vorticity generated is small (since a smaller time constant results in a smaller characteristic Reynolds number and therefore a smaller vorticity) and very localized, and it diffuses faster without making an impact on the energy transferred to the cylinder.Therefore, the VIV amplitude is small for low values of the time constant and it increases with increasing time constant due to the increase in the strength of vortices.As the time constant is increased, the vorticity diffusion decreases due to the reduction in viscous dissipation, and the vortices advect farther until reaching the Newtonian limit.

Response of a 1DOF cylinder in the flow of shear-thickening fluids
In this section, we consider the response of a 1DOF flexibly-mounted cylinder placed in the flow of shear-thickening fluids.The rheology of the shear-thickening fluid has been described using the Carreau model as shown in figure 1(b).We have used a zero-shear-rate

Lock-in for shear-thickening fluids
We have conducted VIV simulations for shear-thickening fluids for time constants varying from λ = 1.5 s (Cu ≈ 10) to λ = 650 s (Cu ≈ 4293).The dimensionless amplitude and frequency of oscillations, as well as the CF force coefficient and frequency, are shown in figure 10.Since the incoming Reynolds number, Re 0 , in this case is larger than the minimum required to observe VIV in a Newtonian fluid, we have added the response of the system in a Newtonian fluid to figure 10 as well (λ = 0, Cu = 0).The response of the system in shear-thickening fluids is qualitatively very similar to the Newtonian response; however, the lock-in range is shifted to the right and the amplitude of oscillations is decreased as the time constant is increased.The shift of the onset of the lock-in range to the right is due to the fact that the Strouhal number decreases with increasing shear-thickening effects, opposite to what we discussed for the shear-thinning fluid in § 4. For shear-thickening fluid, we observe that the largest amplitude of oscillations occurs at the smallest value of the time constant, as opposed to what we observed for the shear-thinning fluid.For all cases where oscillations are observed, the shedding frequency and the oscillation frequency are synchronized and lock-in is observed.The width of the lock-in range and the VIV amplitude decrease with increasing time constant.Oscillations are observed for λ values as high as λ = 600 s (Cu ≈ 3962), although for a very small range of reduced velocities.When the time constant increases, the shear-thickening effect becomes more prominent, and the shedding frequency decreases.Then the shedding frequency reaches the natural frequency of the system at higher dimensional flow velocities, and therefore synchronization (and the onset of lock-in) occurs at a higher reduced velocity, as observed in the plots of figure 10.As the time constant is increased, the dimensionless oscillation frequency, f * , decreases slightly, and while it stays close to one, for larger values of the time constant (i.e.λ ≥ 200 s (Cu ≥ 1321)), the dimensionless frequency never reaches one.This implies that for these larger values of λ, the phase jump from 0 • to 180 • is not observed and the flow forces stay in phase with displacement for all reduced velocities.

Wake for shear-thickening fluids
The flow pattern in the wake of the cylinder is shown by plotting normalized vorticity for selected values of time constants at a constant reduced velocity of U * = 6.Similar to the wake in a shear-thinning case (figure 9), for all cases where oscillations are observed, a 2S shedding pattern is observed in the wake of the cylinder.At the lowest time constant (figure 11e), the wake looks very similar to the wake of a cylinder undergoing VIV in a Newtonian fluid.With increasing time constant (moving up from figure 11e to figure11a), shear thickening causes the viscosity to increase in areas of high shear rates, limiting the magnitude of the vorticity.Thus the maximum vorticity in case of shear-thickening fluids is less than that in the case of shear-thinning fluids, but enough to shed the vortices and cause oscillations.The extent of the wake is not limited by the shear-thickening effect since the diffusion becomes less dominant moving away from the cylinder, and the vortices are swept downstream by advection.
For each shear-thickening case on the left in figure 11, we show a shear-thinning case at approximately the same Re char value on the right to compare the wake in detail.The Reynolds number at the inlet is Re 0 = 200 for the shear-thickening cases and Re 0 = 15 for the shear-thinning cases.In figure 11, the double-headed arrow upstream of the cylinder shows the peak to peak amplitude of the cylinder's oscillations.The comparison shows how the shear-thinning and shear-thickening rheology of the fluid lead to a significantly different wake despite the same Re char and similar amplitudes of oscillations, A * .The sizes of the vortices and the recirculation bubble in shear-thickening cases are larger than those in shear-thinning cases.For the same Re char , the extent of the wake is longer in shear-thickening cases, since the advection becomes dominant moving away from the cylinder.In the shear-thinning cases, the maximum vorticity is larger than that in shear-thickening cases; however, the vorticity is diffused faster in shear-thinning cases, which results in wakes that extend much less than those in the case of shear-thickening fluid.constant characteristic Reynolds number for all cases.For all these cases, the reduced velocity is kept constant at U * = 6.This particular numerical experiment is designed for understanding the effect of the power-law coefficient, n, of the shear-thinning and shear-thickening fluids on the amplitude of oscillations and the wake structure when Re char is kept constant.The histogram and the spatial distribution of the local Reynolds numbers for selected shear-thinning and shear-thickening cases are shown in figure 13.The distribution of local Reynolds number is quite different for shear-thinning cases compared with shear-thickening cases.The shear-thinning cases are dominated by low Reynolds numbers globally as seen in figure 13(a,b), while high Reynolds numbers are observed only in close proximity to the cylinder.This local concentration of high Reynolds numbers is enough to drive VIV with amplitudes of A * = 0.45 and A * = 0.43 for shear-thinning cases of figures 13(a) and 13(b), respectively.In shear-thickening cases (figure 13c,d), the local Reynolds numbers in proximity to the cylinder are much smaller than the incoming Reynolds number of Re 0 = 200; however, the drop of the Reynolds number is not significant enough to suppress VIV.Thus we observe VIV at an amplitude of A * = 0.39 for the shear-thickening cases (figure 13c,d), which is very close to the amplitude of oscillations for the shear-thinning cases.Overall, at the same Re char , when the fluid rheology is changed from shear-thinning (n = 0.1 and n = 0.6) to shear-thickening (n = 1.4 and n = 1.9), the VIV amplitude drops only very slightly from A * = 0.45 to A * = 0.39, and the frequency of oscillations drops from f o /f n = 0.9 to f o /f n = 0.78.The For both shear-thinning and shear-thickening cases undergoing VIV, a 2S shedding pattern is observed in the wake.For shear-thinning fluids, the extent of the wake and the strength of vortices increase gradually with increasing time constant as the advection dominates over diffusion.For small time constants, the vorticity generated is small and very localized, and results in VIV with small amplitudes.The extent of the wake is not limited by the shear-thickening effect since the diffusion becomes less dominant moving away from the cylinder, and advection kicks in.The sizes of the vortices and the recirculation bubble in shear-thickening cases are larger than those in shear-thinning cases.The frequency of oscillations and the shedding frequency gradually decrease when the fluid rheology is shifted from shear-thinning to shear-thickening.In proximity to the cylinder, advection dominates for shear-thinning fluids and diffusion dominates for shear-thickening fluids.Thus the maximum generated vorticity is larger in the shear-thinning cases.The advection-dominated regime away from the cylinder in shear-thickening fluids helps the vortices to move downstream, therefore the extent of the wake is larger in shear-thickening fluids.
When the response of the cylinder is compared using the zero-shear Reynolds number, Re 0 , shear-thinning fluids enhance the oscillations while shear-thickening fluids suppress them.The VIV response is very sensitive to the rheology of the fluid and the shear rate.To consider the effect of fluid rheology and the shear rate in a single parameter, we have defined the characteristic Reynolds number using the viscosity at the characteristic shear rate, U/D.To show how well Re char describes the flow in proximity to the cylinder, we consider shear-thinning and shear-thickening fluids with different combinations of the power-law coefficient, n, and the time constant, λ, such that all of them lead to a constant Re char = 30, and all at a fixed reduced velocity of U * = 6.We show that although the distribution of local Re is quite different for shear-thinning and shear-thickening cases even at a constant Re char , the VIV amplitude is quite similar for all the cases, and Re char can be used to collapse the data for both shear-thinning and shear-thickening cases.We also show that the critical values of Re char for the onset of VIV are Re char,crit ≈ 22 for shear-thinning fluids and Re char,crit ≈ 18 for shear-thickening fluids, both comparable with the critical Reynolds number for the onset of VIV in a Newtonian fluid, i.e.Re crit ≈ 19.

Figure 2 .
Figure 2. Schematic of the domain with mesh and boundary conditions.

Figure 3 .
Figure 3.The dimensionless (a) amplitude, A * , and (b) frequency, f * , of the VIV response for a cylinder free to oscillate in the CF direction and placed in Newtonian flow, found in the present study and the published results of Borazjani & Sotiropoulos (2009) and Ahn & Kallinderis (2006) at Re = 150, m * = 2 and ζ = 0.
Figure 4.The dimensionless (a) oscillation amplitude, A * , and (b) oscillation frequency, f * , as well as (c) the force coefficient in the CF direction, C y , and (d) the force frequency, f * Cy , versus the reduced velocity, U * , for shear-thinning fluids with different time constants, λ.

Figure 5 .
Figure 5. Histograms of the local Reynolds numbers around a cylinder undergoing VIV in a shear-thinning fluid with (a) λ = 0.4 s (Cu = 3), (b) λ = 1.5 s (Cu = 12), (c) λ = 5 s (Cu = 40), and (d) λ = 20 s (Cu = 158), all at U * = 4.The incoming Reynolds number is Re 0 = 15 for all cases.All the cells inside the bounding box (dashed blue rectangle) of size 15D × 7D have been used to create the histograms.A snapshot of the spatial distribution of the local Reynolds numbers is shown for each case.The bin size is 5 in these histograms.local Reynolds number resulting from shear thinning of the fluid close to the cylinder wall.Based on this definition of the characteristic Reynolds number, VIV is not observed for Re char = 30, but is observed for Re char = 59 and larger, indicating that similar to what has been observed for the Newtonian case, a critical characteristic Reynolds number exists for the onset of VIV in shear-thinning fluids.The shedding of vortices at Re 0 = 15 in a shear-thinning fluid is purely a result of the shear-thinning effect.How far these vortices are sustained in the wake of a cylinder, and their strength, depend upon the type of shear-thinning fluid.If the shear-thinning effect for a fluid is not strong enough (i.e. the time constant of the fluid is not large enough), then the vortices are observed only in close proximity to the cylinder.The strength of these vortices is not enough to generate large-amplitude oscillations of the cylinder as observed in a Newtonian fluid.Thus we observe reduction in the amplitude of oscillations as the time constant is decreased in a shear-thinning fluid, as observed in figure4(a).934A39-11

Figure 6 .Figure 7 .
Figure 6.The dimensionless (a) amplitude, A * , and (b) frequency, f * , of oscillations, as well as (c) the CF force coefficient, C y , (d) the ratio of the third harmonic to the first harmonic force in the CF direction, C y,3 /C y,1 , and (e) the phase difference between the CF displacement and the CF force, φ, for U * = 4 and U * = 6 versus the time constant, λ.

Figure 8
Figure 8.(a) Dimensionless amplitudes of the CF response, A * , versus the time constant, λ, for simulations with zero and non-zero initial conditions (IC) at U * = 4. Vorticity fields for λ = 0.6 s (Cu ≈ 5) when (b) zero or (c) non-zero initial displacement is given.The incoming Reynolds number is Re 0 = 15 for all cases.In the figure, ω z is the z-component of the vorticity.

Figure 9 .
Figure 9. Wake patterns for shear-thinning fluids at U * = 4 (a-e) and U * = 6 ( f -j) at different λ and Re char values.For all cases, the snapshot is taken when the cylinder is at the centre and moving up.The incoming Reynolds number is Re 0 = 15 for all cases.

6.Figure 10 .λ
Figure 10.The dimensionless (a) amplitude, A * , and (b) frequency, f * , of oscillations as well as the (c) CF force coefficient, C y , and (d) frequency, f * Cy , versus the reduced velocity, U * , for shear-thickening fluids with different time constants.The power-law coefficient is kept constant at n = 1.2.The incoming Reynolds number is Re 0 = 200 for all cases.

Figure 11 .
Figure 11.Wake patterns for shear-thickening (a-e) and shear-thinning ( f -j) fluids at U * = 6, at different λ and Re char values.The cases for comparison are chosen in such a way that in each row the characteristic Reynolds numbers for the shear-thickening and shear-thinning cases are approximately the same.The incoming Reynolds numbers are Re 0 = 200 and Re 0 = 15 for shear-thickening and shear-thinning cases, respectively.For all cases, the snapshot is taken when the cylinder is at the centre and moving up.

Table 1 .
System parameters used in the simulations.
(Boersma et al. 2021)Kou et al. 2017;Dolci & Carmo 2019;Boersma et al. 2021)ical Reynolds number range, i.e. for Re < 47(Mittal & Singh 2005;Kou et al. 2017;Dolci & Carmo 2019;Boersma et al. 2021).The minimum Reynolds number needed to observe subcritical VIV is Re = 18(Kou et al. 2017), and VIV has been observed experimentally for Reynolds numbers as low as Re = 19(Boersma et al. 2021).As a result, with an incoming Reynolds number of Re 0 = 15 in the present work, any observed oscillations will be purely due to the shear-thinning effects of the fluid.In a Newtonian fluid, when VIV starts at Re = 19, its response amplitude increases with Reynolds number up to Re = 33, after which the response reaches a plateau(Boersma et al. 2021).
4.1.The observed response for different time constantsFigure