Generalized thick strip modelling for vortex-induced vibration of long ﬂexible cylinders

We propose a generalized strip modelling method that is computationally eﬃcient for the VIV prediction of long ﬂexible cylinders in three-dimensional incompressible ﬂow. In order to overcome the shortcomings of conventional strip-theory-based 2D models, the ﬂuid domain is divided into “thick” strips, which are suﬃciently thick to locally resolve the small scale turbulence eﬀects and three dimensionality of the ﬂow around the cylinder. An attractive feature of the model is that we independently construct a three-dimensional scale resolving model for individual strips, which have local spanwise scale along the cylinder’s axial direction and are only coupled through the structural model of the cylinder. Therefore, this approach is able to cover the full spectrum for fully resolved 3D modelling to 2D strip theory. The connection between these strips is achieved through the calculation of a tensioned beam equation, which is used to represent the dynamics of the ﬂexible body. In the limit, however, a single “thick” strip would ﬁll the full 3D domain. A parallel Fourier spectral/ hp element method is employed to solve the 3D ﬂow dynamics in the strip-domain, and then the VIV response prediction is achieved through the strip-structure interactions. Numerical tests on both laminar and turbulent ﬂows as well as the comparison against the fully resolved DNS are presented to demonstrate the applicability of this approach.


Introduction
Vortex-induced vibration (VIV) of long flexible cylinders is widely encountered in various engineering fields. For instance, string-like marine risers exposed to strong currents may experience VIV with large amplitudes, as a consequence of significant interactions of vortex shedding and structural dynamics. Since such vibration, particularly in deep water, has the potential to cause severe fatigue failure, how to accurately predict response features of flexible bodies and understand associated basic flow dynamics has drawn sustainable attention in recent decades (Williamson and Govardhan, 2008;Wu et al., 2012).
Two different ways, respectively, of semi-empirical model and computational fluid dynamic (CFD) simulation, are available in the literature for predicting VIV of long flexible cylinders, according to the method of acquiring hydrodynamic forces exerted on the slender bodies. In the semi-empirical method, such as the wake-oscillator model (Facchinetti et al., 2004;Gabbai and Benaroya, 2005;Chaplin et al., 2005), the principle features of vortex shedding are usually described by a dynamic system, which is derived based on experimental results rather than physical laws. As a consequence, the semi-empirical models are unable to capture the intrinsic physical mechanisms involved and may lack predictable capability to analyze new configurations or scenarios. In CFD simulation, the Navier-Stokes equations are employed to solve the flow dynamics in the physical system and hence compute the fluid forces acting on the bodies, which further couple with structural dynamic model to simulate the interactions between moving bodies and ambient time-dependent flow fields.
With regards to the CFD modelling, most of the full 3D simulation of long flexible cylinders reported in the literature have been carried out for the flow regime where turbulent effects were not significant. The main reason is that full resolutions of cylinder flow with large aspect ratio at high Re are extremely time consuming with current computer power. In particular, Karniadakis (1996, 1997) presented numerical results on VIV of flexible cable at Reynolds numbers of 100 and 200. Later on, Evangelinos and Karniadakis (1999) and Evangelinos et al. (2000) predicted dynamic responses of a flexible cylinder at Re = 1000 using the spectral/hp element-Fourier DNS method. A short cylinder was examined in their simulations under the consideration of computational feasibility and physical reality. A numerical simulation of a much longer flexible cylinder with the same range of Reynolds number was reported in Lucor et al. (2001). The flexible body had an aspect ratio greater than 500 and only larger flow structures were captured due to the very low resolution along the spanwise direction. Bourguet et al. (2011Bourguet et al. ( , 2013 conducted DNS in simulations of VIV of long flexible cylinders within oncoming sheared flows with the maximum Reynolds number up to 1100. The flexible structure was modeled as an extensible tensioned beam with an aspect ratio of 200. Holmes et al. (2006) conducted a numerical study of cylinder VIV with an aspect ratio over 1400 in fully 3-D turbulent flows. In their simulations, they made significant compromises between mesh resolution and computational effort. Even so, the size of the unstructured meshes used to discretize the entire domain still reached 10 million. Constantinides and Oakley (2009) performed similar CFD simulations on a cylinder with an aspect ratio above 4000.
As a more computationally efficient model, strip theory-based VIV modelling technique has also been proposed previously to predict VIV for higher Reynolds number flows (Huang et al., 2011;Sun et al., 2012;Graham, 2001, 2004a,b). In the strip theorybased model, the fluid flow solution is obtained on a series of 2D computational planes (also referred to as "strips") along the cylinder's axis. These strips are coupled with each other through the structural dynamic model of the cylinder, and then the VIV response prediction is achieved by the strip-structure interaction.
Through this model, a large 3D problem can be reduced to a series of smaller 2D problems and decreases the computational cost to a more manageable level. However, an inherent disadvantage of such modelling is that it is unable to resolve the turbulence effects and much of the three-dimensionality of the flow around the body. Motivated by this problem, in this paper we have developed a novel strip model in the context of high-order spectral/hp element framework−Nektar++ (Cantwell et al., 2015). A distinct feature of this approach is that we have constructed a three-dimensional DNS model with local spanwise scale for each strip, and thereby addressed the main objection to the traditional 2D strip model. The outline of this paper is as follows. Section 2 presents the proposed strip model in detail and section 3 shows the numerical method used to solve this fluid-structure interaction problem. In section 4 we present some results for laminar and turbulent flow regimes, including a direct comparison with a fully resolved case, respectively, and in section 5 we discuss on the scale of the width of strips. Some conclusions are drawn in section 6.

"Thick" strip based VIV model
We consider a three-dimensional incompressible flow past a long slender and flexible cylinder with an aspect ratio of L c /D 1, where L c is the spanwise length of the cylinder and D is the cross-sectional diameter. Further, we assume global periodicity along the cylinder's length on both the flow and structure variables. In the 2D strip theory, it is assumed that the flow is locally two-dimensional without spanwise correlation, which allows the problem to be split into various 2D planes. A consequence of 2D strip solution under this assumption is that it is unable to reflect the influence of wake turbulence on the structural dynamics. In order to overcome this problem, a spanwise scale is locally allocated to each of the strips in the current approach, so that the spanwise velocity correlation is reconstructed locally in the flow field. We propose this as a natural extension to the standard strip theory since the largest spanwise wavelength of a turbulent wake (which is of the order of the diameter) is typically much shorter than the excited vibration mode typically observed in a high-aspect ratio cylinder undergoing vortex wave interaction (of the order of the spanwise length). A sketch of the proposed approach is plotted in Fig. 1. In particular, this model lets the fluid domain be divided into N strips with non-dimensional thickness of L z /D evenly distributed along the spanwise (z ) direction. The gap between the neighboring strips, represented by L g , satisfies the relation L c = L z + (N − 1)(L z + L g ). Since the strip in this model has finite scale in the z -direction, we have named it "thick" strip to distinguish it from the conventional 2D strip approximations. Next, the flow dynamics are modeled by a series of incompressible Navier-Stokes equations in an inertial coordinate system x (t) = (x (t), y (t), z ). Then, a sequence of strip-domains is given by and the governing equation is written over a general local domain Ω n in non-dimensional form as where the vector u n = (u n , v n , w n ) denotes the fluid velocity and p n is the fluid's dynamic pressure for the n th strip, and ν is the kinematic viscosity associated with Reynolds number, defined as Re = U ∞ D/ν. U ∞ is the incoming velocity and ρ f is the density of the fluid. The governing equations (2) and (3) are supplemented by boundary conditions of either Dirichlet or Neumann type. For the Dirichlet type, the velocity is imposed in general as where g n is the prescribed boundary value, while for the Neumann type neither velocity nor pressure is known.  A linearized tensioned beam model is employed to govern the dynamics of the flexible structure and can be expressed by the following equation where ρ c is the structural mass per unit length, c is the structural damping per unit length. T is the cylinder tension and EI is the flexural rigidity. F(z , t) denotes the vector of hydrodynamic force (per unit length) with two components of (F D , F L ), respectively of the inline (IL) and crossflow (CF) directions, exerted on the structure's wall. The structural displacement ξ(z , t) has two components which will be denoted as ξ and η, respectively, corresponding to the IL and CF directions. In the following implementation we discretize the spanwise direction of the structural model with a Fourier expansion, so that the partial differential equations reduce to a set of ordinary differential equations, which are further integrated using the Newmark-β method.
How to construct the connection between the discrete strips is an essential aspect of this prediction model, and, typically, we represent it implicitly through the structural dynamic model. As regards to the dynamic continuity condition, the instantaneous hydrodynamic force per unit span length distribution along the n th strip is expressed as where n is the outward-pointing unit vector normal to the surface of the cross-section of the cylinder, represented by Γ .
Under the scenario of N = 1 and L z = L c , the fully resolved case is recovered from the strip modelling, since in this limit a single thick strip would fill the full 3D domain. For the situation of N > 1 with L z > 0 and L g > 0, we need to build the coupling between the strips by filling the gaps between them. To implement this, we need sampling F n (z , t) at equispaced points {z k } n : where K is the number of points (uniformly distributed) over each strip. Further, a set of points is obtained from union operation for all strips as {z nk } := N −1 n=0 {z k } n , which is in general non-uniformly distributed along the spanwise direction. Next, we define the discrete Fourier transform (DFT) as Here, the sampling points in the Fourier domain are equispaced and given by {ζ A reconstruction in the physical domain is then implemented from the following inverse discrete Fourier transform (IDFT), which is defined as The set of points, which is evenly sampling the physical domain after an inverse Fourier transformation, is It is important to remark here that a standard fast Fourier transform (FFT) algorithm, which allows the calculation to be computed in O(N · Klog(N · K)) operations rather than O((N · K) 2 ), is no longer applicable for the computation of the DFT defined by (7) in the situation of K > 1, due to the nonequispaced sampling sequence {z nk }. Over the last decade, a number of numerically efficient algorithms have been developed to overcome this limitation and are often referred to as nonuniform FFTs (NUFFTs). We refer the interested reader to (Greengard and Lee, 2004;Fessler and Bradley, 2003) for the implementation details of the NUFFTs.
Nevertheless, a standard FFT algorithm is still applicable for the case of K = 1, that is, when a single point is used to describe the kinematics of each strip. Under this situation, we take the homogeneous average of the hydrodynamic forces over the strip-domain, such that Equivalently, the above equation can be treated as the zeroth-order wave component of forces for each strip.
As a consequence, the sampling in the DFT then recovers to be a equispaced sequence in the physical domain.
However, such treatment requires an assumption of λ c L z /D, where λ c is the minimum wavelength of the oscillation of the cylinder, since this approximation will be eventually missing some information related to higher order derivatives with respect to z -direction.
In the current work, we consider only the reconstruction based on a standard DFT transform in uniform spaced sampling with K = 1. Because in this simplification the homogeneous periodicity is satisfied in the spanwise direction to the local flow within individual strips, therefore it enables the use of spectral element-Fourier methods to the flow simulations in the strip-domains independently.
The kinematic continuity for the n th strip is stated simply as where ξ n represents the displacement vector of the n th strip. In order to associate to the approximation of (9), the above condition is reduced to The simplified condition of (11) then indicates that the discrete strips with this condition would move as a solid region and so only allow the IL and CF motions.
The coupled system of incompressible flow and moving bodies are solved in the Nektar++ spectral/hp element framework (Cantwell et al., 2015;Web, 2014). A body-fitted coordinate formulation is used to accommodate the moving boundary, therefore, the moving body's influence upon the surrounding flows is embedded in adding acceleration term to the Navier-Stokes equations with non-inertial transformation. More details will be shown in Section 3.

Coordinate Mapping
Solving the flow around flexible cylinders involves fluid-structure interaction, which requires specific treatment of moving boundary conditions. A relatively simple approach that can be used to tackle this situation is setting the coordinate system to be fixed in the moving body (Newman and Karniadakis, 1996;Li et al., 2002). One significant advantage is that it avoids the dynamic re-meshing of the computational domain, but it comes at the cost of increasing the complexity of the governing equations. The coordinate transformation is required to map the flow fields from the inertial domain Ω n (t) to the body-fitted domain Ω n , which is given by We define the mapping with respect to the cylinder's motion as According to this linear mapping, the field variables of velocity u n = (u n , v n , w n ) and p n in the body-fitted system are easily obtained as As the coordinate transformation is applied to the governing equations (2) and (3), the Navier-Stokes and continuity equations, then become with the Dirichlet boundary condition of equation (4) in the body-fitted system transformed as In particular, on the wall of the body, it becomes u n = 0. A high order outflow boundary condition is applied to supplement the governing eqs. (15) and (16), as where n is the outward-pointing unit vector normal to the outlet, and |u n | denotes the magnitude of velocity u n . The smoothing function S 0 (n · u n ) is chosen by Dong et al. (2014) as where U 0 is the characteristic velocity scale, and δ > 0 is a chosen non-dimensional positive constant that is sufficiently small.
The additional forcing term f n = (f nx , f ny , f nz ) appears in the transformed momentum equation (15) to take account of the acceleration of the transformed axes in each strip. The three components of this term are given as Again, if we consider only the zeroth-order component, it reduces to

Fourier Spectral/hp Element Method
We have assumed that the flow variables for each strip are homogeneous in the spanwise direction with a periodic length (i.e., strip thickness) L z . Based on this assumption, a Fourier expansion is introduced in the z-direction to the flow variables u n and p n as with β = 2π/L z . m is a Fourier mode index and M the number of modes in the Fourier expansion. After applying the Fourier transformation to eqs. (15) and (16), we obtain a set of uncoupled two-dimensional equations for each mode, where N(u n ) m represents the Fourier mode of the convective term andf nm is the corresponding Fourier mode of the forcing term. In order to avoid formation of convolution sums, the Fourier transformation of the nonlinear terms is applied in physical space. In eqs. (24) and (25), we defined the differential operators as The stiffly stable high-order splitting scheme proposed by Karniadakis et al. (1991), also referred to as a velocity correction scheme (Guermond and Shen, 2003), is employed for time integration of the eqs. (24) and (25), and it hence consists of the following four steps: (1) Compute a first intermediate velocity field by calculating explicitly the advection term and forcing term in the momentum equation, which are extrapolated from the previous time steps k−q (q = 0, 1, ..., J −1).
(2) Solve a Poisson equation to obtain the pressure solution at the time step k + 1 with consistent Neumann boundary condition prescribed as and a high order Dirichlet boundary condition on the outlet boundarŷ (3) Calculate a second intermediate velocity field by including the pressure gradient term with a high order Neumann boundary condition on the outflow boundary Here J denotes the time integration order, and α q , β q and γ 0 are the corresponding coefficients of the J-th order backward differentiation formulas. In this work, we set J = 2 to obtain second-order accuracy for time integration of the fluid solver.
Next, for application of spectral/hp element methods to the above split scheme, we formulate the weak forms of eqs. (27) Similarly, we define a velocity test function ϕ from the same space such that ϕ ∈ H 1 (Ω 2d ) and ϕ| ∂Ω 2d Taking the same procedure as that applied for the Poisson equation, we can obtain the weak form of the Helmholtz equation for the velocity as follows, Finally, the flow domain is partitioned into a mesh containing elements, such that Ω 2d = Ω 2d e and Ω 2d e1 Ω 2d e2 = ∂Ω 2d e1e2 is the empty set or an interface of two adjacent elements. Spatial discretization proceeds by the approximation of unknown flow variables by polynomial expansions element by element. Here we choose the test functions to be the same as the trial functions, so the flow solution can be represented aŝ where P is the order of polynomial,û δ j andp δ j denotes model coefficients of velocity and pressure. Substituting the polynomial approximation into the elemental-wise weak form of the governing equations, and assembling the local contribution to the solution through the assembly operator, we can form a global solution that can be solved by both direct or iterative solvers (Karniadakis and Sherwin, 2013).

Fluid-Structure Coupling
Similar to the Fourier transformation applied to the governing equations of fluid flow, a periodic assumption is imposed for the cylinder's motion with respect to the full-length scale L c . The motion variables and forces then are expressed as a complex Fourier series, In the following, the superscript and the second subscript s are omitted for simplicity of the notation. The partial differential equation of the tensioned beam model then decouples into a set of ordinary differential equations (38), which are easily solved by a second-order Newmark-β method, A partitioned approach is adopted to solve the fluid-structure coupling system, which requires the exchange of transmission conditions on the interface between the fluid solver and structure solver. A simple and efficient strategy of loose coupling is applied for the time integration of the coupled system. Sequential solution of the fluid and structural fields in a staggered manner causes numerical instability in this scheme due to the so called artificial added mass effect, and that becomes significant particularly at low ratios of structure density to fluid density (Förster et al., 2007). In order to address this problem, Yu et al. (2013) and Baek and Karniadakis (2012) developed a generalized fictitious method to stabilize the partitioned scheme by introducing either fictitious pressure to the fluid solver or fictitious mass and damping to the structure solver.
In the current work, we introduce fictitious mass and damping to the tensioned beam model. Considering the fictitious acceleration and damping terms, which are added to both sides of eq. (38), it becomes where ρ a and c a are fictitious density and damping coefficient, respectively. The term in the square bracket in eq. (39) is approximated explicitly from the values computed in the previous time steps. A combination of fictitious mass and damping (Baek and Karniadakis, 2012) of (ρ a /ρ c , c a /c) = (3.0, 0.1) gives a stable solution for the mass ratios considered in this work.

Hybrid Parallelisation
A hybrid parallelisation approach (Bolis, 2013), in which FFT transposition (modal parallelisation) and mesh decomposition (elemental modal parallelisation) are used concurrently, is adopted for the current model.
The Cartesian virtual topology (rows and columns) shown in Fig. 2 facilitates the conceptual handling of the two parallel techniques. Practically, we split the N strips across a set of P Z processors used for the FFT parallelisation, and then the FFT transposition is implemented by re-splitting the N 2D planes across a subset of P Z /N processors for each strip independently. As can be observed in the graphical illustration of the communicator, the mesh partitions are distributed across columns, requiring a specific communicator for each row. The decomposition of the 2D plane is accomplished on the root process just once, i.e. P0 in Fig. 2. Subsequently, each partition will be sent to the appropriate process. For example the first partition in the strip 0, the yellow one, will be sent to P5, and in the strip 1, P10 will be send to P15. The second partition, the blue one, will be sent to P1, P6 for the strip 0, and P11, P16 for the strip 1 and so on. It also appears from this particular example that P1 is not required to communicate with P2, since the respective partitions do not share DOFs (actually in each row the blue process does not need to communicate with the green one). The composition of all the partitions on a row produces a full 2D plane, but each row refers to a different plane (or set of consecutive planes in the most general case). Therefore operations across planes require another set of communicators, i.e. the column communicators. Before and after an FFT, data needs to be reordered.

Numerical Tests
We now consider some computational validations of the "thick" strip theory by first considering a lower Reynolds number (Re = 100) problem as a fully coupled, single domain problem as well as a higher Reynolds number flow at Re = 3900 where we discuss how the "thick" strip modelling works for the turbulent flow.
In the low Reynolds number test the "thick" strips become equivalent to the classical "thin" strips model since the underlying flow is two-dimensional. At Re = 3900 the flow is certainly three-dimensional and so the "thin" strip model cannot capture these wake features. In the following study we will consider the cases without the addition of a turbulence model, (although as discussed below we will add some high frequency stabilization). We do recognize that in previous "thin" strip models ( At Re = 100, a 3D fully resolved simulation is conducted firstly to verify the code implementation against the previous numerical results in Newman and Karniadakis (1996) and to be used to compare against the results of strip theory-based modelling. In the laminar flow simulation, the vibration of the cylinder is constrained to occur only in the CF direction. The cylinder has a length ratio of L c /D = 4π with the tension parameter set to T = 8.82, choosing the same value reported in Newman and Karniadakis (1996).
The structural damping and flexural rigidity are omitted to simplify the analysis. Two different mass ratios, m * , defined as the ratio of the cylinder's mass per unit length to the fluid density (i.e., m * = ρ c /ρ f D 2 with D = 1), respectively with the values of 1 and 2 are considered.
In the turbulent flow simulation at Re = 3900, the cylinder has an aspect ratio of L c /D = 32π(≈ 100.53).
It is pinned at both ends and free to vibrate in both the CF and IL directions. The mass ratio considered is m * = 1.17, a typical value of that can be encountered in the context of ocean engineering. The tension and the bending stiffness are set equal to 600 and 20.8, respectively. The structural damping coefficient is set to be c = 0.022.

Computational domain and conditions
A C-type domain is used in the computation that has its semicircle upstream boundary with a diameter of 20D, its lateral boundaries located at 10D, and extending over 12.5D downstream. On the inlet and lateral boundaries, a Dirichlet velocity condition is specified as (u n , v n , w n ) ∞ = (U ∞ , 0, 0) with U ∞ = 1, while on the outlet a high-order Neumann boundary condition is determined from eq. (34). The 2D spectral element mesh employed for the computation is also shown in Fig For the multiple strip modelling at Re = 100, the strips with thickness ratio of L z /D = π/8 are evenly decomposed along the z-direction, and for each strip two planes (one pair of complex Fourier modes) are used for the local flow resolution.
For the flow at Re = 3900, we carried out a series of simulation with different strip thickness, which has aggregate impact on the performance of the modelling. 16 strips are distributed along the cylinder's span as shown in Fig. 3b. The spanwise thickness of the strip is chosen to be L z /D = {π/16, π/8, π/4, π/2, π, 2π}.
For the case with L z /D = π, a total of 48 planes ( simulation becomes fully saturated. Fig. 5 shows the displacement evolution for different mass ratios after the VIV response turns into traveling waves. A similar observation was reported by Newman and Karniadakis (1996).

2D strip modelling at Re = 100
The transition from a standing wave pattern to a traveling wave provides the opportunity to compare the strip theory during these two stages. Therefore in this test, we apply the "thick" strip theory to the low Reynolds number problems. Since the three-dimensional wake structures induced by this problem have a larger wavelength than the local "thick" strips (see later Fig. 8), we are essentially solving 2D strips similar to the classical "thin" strip theory. As seen in the temporal evolution of VIV responses plotted in Fig. 6, the model reproduces the standing wave that is also observed in the full 3D model. However, the wave pattern in the drag and lift coefficients are somewhat different from their fully resolved counterparts. This is because, the communication between the flow dynamics in different strips is achieved, not directly, but through the structure dynamics. This then accounts for the principal difference between the full 3D and strip theory-based models. The predicted vibration amplitudes are 0.643 (m * = 1) and 0.589 (m * = 2), respectively, which are approximately lower by 11% and 15% with respect to those in full 3D resolutions. There appears to be a relatively larger difference in prediction of the hydrodynamic forces. Taking lift coefficient as an example, its magnitudes for the cases of m * = 1 and 2 are 1.415 and 0.282, respectively, whereas the full 3D simulations obtained the values of lift coefficient as 1.279 and 0.482, showing differences of 9.6% and 41.5% respectively.
When m * = 1, the cylinder's response converges to a traveling wave pattern that is consistent with the full resolution results. In contrast, at m * = 2, the mode may converge to different states, depending on the initial conditions. As shown in Fig. 7, if the simulation runs from a standing wave perturbation with  The wake features of the 2D strip model at m * = 1 are illustrated in Fig. 8. The spanwise vorticity contours show a consistent convergence of the flow fields with the increase of the number of strip. The standing wave response of the cylinder generates an interwoven structure of vortex shedding in the wake, which is consistent with the wake observed in the full 3D simulations (Newman and Karniadakis, 1996).
Similar to the full 3D results again, the traveling wave responses generate a coherent wake structure with oblique shedding of spanwise vorticity. Therefore, it can be observed that although the wake produced within strip-domain is two dimensional, the intrinsic three dimensionality of the wake structure as a whole Overall many of the full 3D features are captured by what is still essentially a classical 2D "thin" strip theory.

"Thick" strip modelling at Re = 3900
We now turn our attention to a "thick" strip model at a Reynolds number where the wake dynamics vary within a single "thick" strip. Therefore in this section, we investigate the "thick" strip modelling of the flow at Re = 3900. When using a limited resolution it is necessary to stabilize the under-resolved simulations and so we apply a spectral vanishing viscosity (SVV) technique (Tadmor, 1989;Karamanos and Karniadakis, 2000;Kirby and Sherwin, 2006) by introducing a term as an artificial dissipation mechanism to operate only on the highest resolved frequencies of wave space solution. This has the form of a convolution: where is the viscosity amplitude and Q Dim (Dim = 2) is the viscosity kernel. In the framework of Nektar++, this term is incorporated into the velocity splitting scheme of the incompressible Navier-Stokes equations.
In Figs. 9a and 9b, we show the time evolution of the displacement responses predicted by the "thick" strip models with L z /D = π. It is obviously illustrated that a standing wave in second harmonic mode is excited in the CF vibrations over the length of the cylinder. We plot in Figs. 9c and 9d the power spectral densities as a function of spanwise location for the displacement responses shown in Figs. 9a and 9b.
The analysis is performed on the time series after the initial transient effect has died out. From this figure we observe that the frequency profile along the spanwise also illustrates a standing wave response. When compared with the PSD of lift coefficient, the CF vibration is found to be synchronized with the excited lift frequency. On the other hand, the spectrum of the oscillation component of the IL displacement shows three peaks that evenly distributed along the length at a frequency typically twice of that of CF response. The spanwise distribution of a representative XY-trajectory over a time period covering one full vortex shedding cycle is revealed in Fig. 10. A figure of '8' pattern is observed for the two sides of the flexible cylinder, except for the trajectory pattern around the node of the standing wave. with definite nodes and antinodes for the CF displacement of the riser in uniform flow. In particular, the CF response is excited at a single frequency and in a single mode shape. This feature is highly consistent with the current numerical results. In terms of IL response, appreciable contributions at both twice the CF response frequency and at the CF response frequency itself are revealed at some spanwise locations in the contour plot of the power spectral density.
The contour corresponding to the time series of hydrodynamic force distributions along the cylinder and their PSD analysis are presented in Fig. 11. A more detailed discussion on the temporal and frequency characteristics will be provided in the later comparison with the fully resolved case. Finally, the evolution of fully developed wake flows is shown in Fig. 12, which gives a fuller picture of the instantaneous wake features in the discrete strips. Further, the strip coupling is evident in Fig. 12, since the distribution of vortex shedding along the spanwise length shows a second harmonic wave, which is consistent with that observed in the VIV responses in Fig. 9.

Fully resolved simulation at Re = 3900
In this subsection, we present the results of the fully resolved simulation (N = 1 and L z /D = L c /D = 32π) at the same Reynolds number of 3900. A total of 1024 planes (512 complex modes) are employed over the whole length of the cylinder. The simulation was started from a state of fully developed turbulent flow around a stationary cylinder, and approximately 10 cycles of vortex shedding has been developed in the wake. A total of 16 sine modes were taken into account for the cylinder's responses and a cubic spline was used for the interpolation in the mapping from the structure to flow domains.
The displacement contours of both the temporal and frequency responses are plotted in Figs. 13. As it can be now clearly seen, the response characteristics in both the time and frequency domains in the fully resolved simulation are in close agreement to those predicted by the strip' modelling, see Fig. 9. The calculated spatio-temporal root mean square (ξ RM S ) and the spatio-temporal standard deviation (σ 2 ξ ) of the IL displacement are ξ RM S = 1.762 (1.868 for the strip modelling with 5.6% error) and σ 2 ξ =0.143 (0.149 for the strip modelling with 4% error).
We now turn our attention to the contours of the hydrodynamic force coefficients and the associated PSDs are shown in Fig. 14. The spatio-temporal average of the full simulation drag force isC D = 2.266, while for process. We however note that the range of the frequencies are difference in the drag PSD plots (Fig. 14c) of the full simulation and the thick strip theory (Fig. 11c). Whilst the fundamental frequencies are similar there is notably more energy in the fourth harmonic in the full simulation (Fig. 14c) versus a number of high frequency content in the thick strip theory simulation (Fig. 11c). We note that these harmonics are present in the lift forces data but since the lift force is dominated by the vortex shedding frequency they are not evident in these plots (note the different scales of the lift and drag PSD). There are a number of potential sources for this discrepancy. First we recall that we are using a zeroth-order approximation of both the hydrodynamic forces and kinematic response on each strip. Second, we do not allow for far wake interaction in the thick strip modelling. Indeed consideration of Fig. 15 (figure of full simulation wake) as compared to Fig. 12 we observe that in the thick strip model more coherent far wake vortices that are most likely associated with the high frequencies observed in the drag coefficient of the strip model.
To summarize we observe good agreement in the displacement characteristics arising from ability to capture near wake force characteristics with sufficient fidelity, particularly in the lift. However there are discrepancies in the far field force characteristics which is not surprising since far wake interactions are restricted by the strips thickness. Nevertheless the interaction of the far wake has a relatively small effect on the dynamics of the cylinder motion which is the primary motivation behind using the thick strip method.
Whilst we accept that far field interactions will be restricted by the strip thicknesses we also acknowledge the strips need to be sufficiently thick to capture the near wake dynamics. In the following section we therefore explore the dependence in more detail.

Thickness of the strips
The dependence of the VIV responses on the strip thickness is dictated by the behavior of the wake dynamics in each strip and the corresponding forces these wakes produced. In Fig. 16, we illustrate the effect of strip thickness on the wake structures through the comparison of the spanwise velocity contour of the 4 th strip-wake. This highlights that the strip thickness is critical to appropriately capture the three dimensionality of wake behavior. For the 'thinner' strips (e.g. L z /D = π/16, π/8), the wake dynamics are similar to those of a 2D strip without a turbulence model. Although the onset of three dimensionality was already triggered in the wake, it is not 'thick' enough to cover the most important correlation length. This is consistent to the width being smaller than that of the Mode B instability (λ zB /D ≈ 0.825) (Barkley and Henderson, 1996;Williamson, 1996;Barkley et al., 2000). It can also be illustrated from this figure that the 'thicker' strips, with L z /D = π/2, π and 2π, are able to capture the most energized spanwise turbulent scales in the strip-wake. The ability to capture the near wake three-dimensional dynamics leads to a dissipation of higher frequency components of hydrodynamic forces as displayed in the PSD of lift forces as a function of spanwise length shown in Fig. 17. This figure also illustrates the sensitivity of the strip thickness since, in the strips with L z /D = π/16, π/8 and π/4, the wake behaves in a manner more similar to that of a 2D strip theory, where the oscillation energy is scattered over a broader band of frequencies. As the three dimensionality is fully developed and Mode B is fully captured for the cases with L z /D = π/2, π and 2π, the energy is concentrated at the vortex shedding frequency. We also noticed that the energy of lift oscillation becomes statistically saturated as the strips fully captured Mode B. As we have noted earlier, 2D strip models would typically incorporate a turbulence model which would likely suppress the vortex pairing phenomena observed in this model at shorter wavelengths. However, such models can be quite diffusive and they formally make a number of assumptions about turbulence, which are not valid in the near wake region of a cylinder.
For strip-structure coupling, we consider only the mean mode (corresponding to K = 1) inside each strip analogous to the classical 2D strip theory. Therefore each strip can be considered as a spanwise rigid block (able to move inline and transversely) that can capture the anisotropic wake breakdown but cannot respond to other three dimensional features of the structure's movement or indeed capture a linear variation of the incoming flow due to our assumption of spanwise periodic strips. Although this goes beyond the scope of this work, one might consider exploring more intricate deformations of the span even under the locally periodic strips we have adopted. In addition, fully three-dimensional blocks could be incorporated into the model but appropriate boundary conditions, such as the use of symmetry conditions, would also need to be investigated.  A generalized strip model has been presented for prediction of VIV responses of long flexible cylinders, which is also highly parallelizable due to only requiring weak coupling between each strip from the structural model. The flows in this model are split into a sequence of discrete strips, which have spanwise scale to locally resolve the three dimensional wakes that have a substantial influence on the VIV responses in turbulent flows.

Conclusions
This removes the need to rely on 2D turbulence modelling of bluff body flows, thereby reducing flow modelling inaccuracies.
Numerical tests were performed to evaluate the performance of this approach for both laminar and turbulent flow regimes. The working principle of the 2D strip model is investigated for laminar flow, since in this case the 'thick' strip model is essentially equivalent to 2D strip theory. It is shown that the 2D model can accurately predict the main features of the VIV responses of flexible bodies, in particular of traveling wave responses, whereas the performance for standing waves is worse since the coupling between the wake modes has been lost.
Tests on the turbulent flows were conducted on a flexible cylinder with an aspect ratio of L c /D = 32π at Re = 3900. A full simulation (N = 1 and L z = L c ) resolved using 1024 planes were also carried out to evaluate the performance for thick strip modelling. Consistent agreement in terms of VIV responses was observed with RMS and standard deviation of the inline response being within 5.5%. The lift force characteristics also demonstrated close agreement as did the mean drag. However the resolution of the unsteady drag charactersitics was more demanding and highlight a discrepancies in the higher frequencies components. Potential sources of this discrepancy could be related to the zeroth-order approximation of the wave components in each strip as well as the restriction on far wake interaction in this model. Finally, the parametric analysis suggests that the strip thickness is critical for this model and the width of strips needs to be larger than the wavelength corresponds to Mode B instability.