CFD Simulation of Vortex Induced Vibration for FRP Composite Riser with Different Modeling Methods

Featured Application: Three different risers, i.e., conventional steel riser, composite riser with orthogonal reinforcements and tailored composite riser, are compared for their VIV characteristics. The effects of 2D and 3D models and fluid–structure interaction (FSI) have been considered and two different modeling methods are established to simulate the FRP composite risers. This paper will serve to further understand the dynamic characteristics of the FRP composite risers and improve the utilization of the FRP composite risers in real projects. We believe present paper would appeal to the researchers working on the design and analysis of composite structures and their optimizations. Abstract: Steel risers are widely used in offshore oil and gas industry. However, the production capacity and depths are limited due to their extreme weight and poor fatigue and corrosion resistance. Nowadays, it is conﬁrmed that ﬁber reinforced polymer (FRP) composite risers have apparent advantages over steel risers. However, the study of vortex induced vibration (VIV) for composite risers is rarely involved. Three different risers (one steel riser and two composite risers) were compared for their VIV characteristics. The effects of 2D and 3D models and ﬂuid–structure interaction (FSI) were considered. The models of composite risers are established by effective modulus method (EMM) and layered-structure method (LSM). It is found that 2D model are only suitable for ideal condition, while, for real situation, 3D model with FSI has to be considered. The results show that the displacements of the FRP composite risers are signiﬁcantly larger than those of the steel riser, while the stresses are reversed. In addition, the distributions of the displacements and stresses depend on the geometries, material properties, top-tension force, constraints, etc. In addition, it is obvious that EMM are suitable to study the global working condition while LSM can be utilized to obtain the results in every single composite layer. (Riser 2) and tailored composite riser (Riser 3) are investigated and compared considering the effects of 2D and 3D models, FSI and different modeling methods for the composite risers. The results show that:


Introduction
The exploitation and utilization of marine resources, such as oil and gas, lead to the development of offshore structures. Risers which act as a transport component between the production platforms and sub-sea wellheads are an indispensable part for the offshore exploitation and transportation system. Currently, risers are normally made from high grade steel. Their extreme weight results in the significant increase of the requirement in top tension forces, thereby leading to the limitation of the numbers of risers that can be attached to current platforms or larger platforms to carry the same number of risers with the increase of depths.
It is well acknowledged that fiber reinforced polymer (FRP) composites have great mechanical properties and lower density; therefore, using FRP composites instead of steel can indeed reduce the weight of individual riser, thereby the requirement of top tension forces, which results in lower operational costs [1]. Meanwhile, FRP composites also have better performance in fatigue, corrosion and thermal insulation compared with steel, leading to lower maintenance costs [2]. It should be noted that the design results of a composite riser can be quite different, even for the same operational situation. The design variables of the geometries of a composite riser basically include liner and composite materials, stacking sequences of laminate, fiber orientations, lamina thicknesses, etc. Therefore, more challenges and complexities of design and analysis are produced for the use of composites in risers.
During the past 40 years, several design studies and projects have been conducted to confirm the feasibility of composite risers. The initial attempt to fabricate and test composite risers was accomplished by Institut Francais du Petrole (IFP) and Aerospatiale of France [3] in the 1980s. In the 1990s, tests of mechanical capacities of the composite riser tube were performed by National Institute of Standards and Technology (NIST) and Advanced Technology Programs (ATPs) [4]. Then, more tests, including field tests and optimized designs of composite risers, were conducted by different research institutes and global enterprises [5][6][7][8][9][10][11][12][13]. These prototypes prove that FRP composite risers have significant benefits in weight saving compared to steel risers. In addition, using optimization methodologies in the designs would lead to more benefits for the composite risers with the same load-carrying capacity.
Vortex-induced vibrations (VIV) are motions induced on bodies interacting with an external fluid flow, produced by-or the motion producing-periodical irregularities on this flow. VIV would cause the fatigue damage of risers, leading to the construction and environmental disaster. Therefore, many studies of VIV have been conducted in the past decades. The study of VIV started with experiment in 1960s using wind tunnel [14]; then, many tests were conducted for the VIV of risers [15][16][17][18]. Besides tests and experiments, researchers proposed different empirical models to solve the VIV problem of risers. Initially, Hartlen and Gurrie built the lift-oscillator model [19], and then this model was developed and applied in the VIV study of flexible slender cylinders [20]. More recently, the wake oscillator models were modified and developed by different researchers according to different parameters and situations [21][22][23][24][25][26]. Nowadays, computational fluid dynamics (CFD) has been employed in the VIV simulation and different models and methodologies have been set up. 2D simulation of risers was done by Huang and Wang [27]. Ji and Liu et al. simulated the VIV response of in-line two risers [28,29]; and a top tension riser was used as an example to examine fatigue safety factor of the riser and the effects of the random coefficient [30].
From the previous studies including tests, theoretical analyses and CFD simulations, the VIV of conventional steel riser were understood. However, the study of VIV for composite risers, the late-model risers, are rarely involved. The complexity and design variety of composite risers enhance the difficulty for the study of VIV for composite risers. A global-local analysis methodology was built to investigate the VIV responses of composite riser under semi-empirical fluid load models [31]. Huang and Chen et al. compared the dynamic performance of composite and steel risers using ABQUS and SHEAR7 [32,33]. These studies normally employed empirical models and no fluid-structure interaction (FSI) was included. The design variety for the composite risers' geometries and the VIV response in every lamina of composite riser were not considered either.
In this paper, the VIV responses of three different risers, i.e., conventional steel riser (Riser 1), composite riser with orthogonal reinforcements (Riser 2) and tailored composite riser (Riser 3), are investigated and compared considering the effects of 2D and 3D models and FSI. The simulation of risers' VIV was achieved by ANSYS 17.0 including modules of Geometry, ACP, Transient Structural, Fluid Flow and System Coupling. Effective modulus method (EMM) and layered-structure method (LSM) were built to model composite risers in ACP and Transient Structural and the flow situation is simulated in Fluid Flow (Fluent). System coupling with the data exchange between risers and fluids were utilized to achieve FSI. The VIV responses, such as coefficient of drag (C D ), coefficient of lift (C L ), displacement and stress under VIV, were obtained and compared for all the three risers with different modeling methods.

Riser Geometries
The internal diameter (I.D.) of the riser joints was fixed at 250 mm. The outer diameter (O.D.) varied from riser to riser depending on the wall thickness of segment determined by the design. To satisfy the same functional requirements, the thicknesses for Risers 1, 2 and 3 were 25 mm, 38 mm and 30 mm, respectively. The steel riser was monolithic, while the composite risers consisted of an inner liner, a composite structural body and an external sacrificial layer ( Figure 1). The geometries of these three risers are from the previous studies [6,10]

Materials for the Risers
High grade steel X80 was selected as the material for Riser 1, while AS4/PEEK composites with PEEK liner were chosen for Risers 2 and 3. The properties of these materials are listed in Table 2. The material properties listed in the table were utilized for steel riser directly, as well as composite risers with layered-structure method (LSM). In terms of the composite risers with effective modulus method (EMM), a MATLAB code (Version 9.1 R2016b, MathWorks, Inc., Natick, MA, USA, 2016) of 3D laminate property theory [34,35] was created for calculating composite tube's 3D effective properties using Equations (1)- (6). Table 3 shows the effective properties of FRP composite risers for effective modulus method.
is the principle elasticity matrix of each composite lamina, S is the effective elastic compliance matrix of the composite laminate, C is the laminate stiffness, and S ij is the principal compliance matrix for transversely isotropic composite lamina. Table 3. Effective properties of fiber reinforced polymer (FRP) composite risers (effective modulus method). Subscripts: x, axial direction of the risers; y, hoop direction of the risers; and z, through-thickness direction of the risers.

FE Modeling
The simulation of risers' VIV was achieved by ANSYS17.0 including modules of Geometry, ACP, Transient Structural, Fluid Flow and System Coupling ( Figure 2). The following geometry was utilized to create the basic geometry of flow field: 10D × 20D for 2D simulation and 10D × 20D × 12.5 m for 3D simulation. The basic geometry of the risers were based on previous design results [6,7,10,13].
Flow situation was simulated by Fluid Flow (Fluent) including flow dimension, flow model, flow properties, flow boundary conditions, etc. The surface size of the flow field is defined by 10D × 20D, where D is the outside diameter of the risers. The riser center is 5D from the inlet. One-year winter storm situation in Gulf of Mexico [36], i.e., velocity U = 0.36 m/s and kinematic viscosity coefficient υ = 1.06 × 10 −6 , was used as the design velocity in this paper. Large eddy simulation (LES) which allows the largest and most important scales of the turbulence to be resolved (small eddies are modeled), while greatly reducing the computational cost incurred by the smallest scales, was selected for the CFD simulation. Resolving only the large eddies allows larger time-step and coarser mesh sizes than direct numerical simulation (DNS). However, LES still requires finer meshes than Reynolds Averaged Navier-Stokes (RANS). Meanwhile, LES has to be run for a sufficiently long flow-time to obtain stable statistics of the flow. As a result, the computational cost (memory and CPU time) involved with LES is normally higher than that in RANS calculations. The inlet, outlet, left and right sides and top and bottom sides were defined as boundary conditions of velocity inlet, outflow, symmetry and wall, respectively. More specifically, velocity inlet consisted of 0.36 m/s (all three risers) for velocity magnitude; 3.78% (Riser 1), 3.75% (Riser 2), and 3.77% (Riser 3) for turbulent intensity (Equation (7)); O.D. of each riser for hydraulic diameter; 1.0 was utilized for flow rate weight in outflow boundary; symmetry boundaries were used for the left and right sides to reduce the extent of computational resource; in the wall boundary, stationary wall was used for wall motion and no slip for shear condition. To use LES model, bounded second order implicit transient formulation must be chosen as solution method.
where I is turbulent intensity and Re is Reynolds numbers. ACP (Pre) was utilized to create the lamination of composite risers (Tables 1 and 2) including material properties for liner and composites, liner and composite lamina thicknesses, fiber orientation, numbers of reinforced layers, and laminate stacking sequences (Figure 3) for the layered-structure method. Transient Structural was utilized to simulate the dynamic response of risers including displacement history, global stress distribution, etc. Besides the geometry of the riser, the input parameters also included the gravity, top-tension force (Table 4), buoyancy (based on the geometry of the riser and the density of the fluid), simple support at the top end, fixed support at the bottom end and the fluid solid interface. ACP (Post) was utilized to present the dynamic stress and deformation response for each lamina of the composite risers for the layered-structure method.
System Coupling was utilized to transfer the data between Transient Structural and Fluid Flow to achieve the two-way fluid structure interaction. In this module, analysis type was set as transient and the duration time was set as 200 s (step size is 0.2 s). For the Co-Sim. Sequence, Fluid Flow (Fluent) was 1 and Transient Structural was 2. Table 4 presents the flow situation and Reynolds numbers (Re) and the tension forces for all the risers. Reynolds numbers are calculated by Equation (8), where U is the flow velocity, D is the outside diameter of the riser, and υ is kinematic viscosity coefficient.
To ensure the accuracy of the simulation for vortex shedding phenomenon, the fluid grid for the 5D × 5D area around the riser was refined to take the data exchange between the riser and flow into account. The dynamic mesh was employed in the refined zone. The dynamic mesh zones involved seven faces, i.e. the six surrounding planes of the refined zone (set as stationary) and the contact surface between fluid and riser (set as system coupling). Smoothing and remeshing were utilized as the dynamic mesh methods. For smoothing method, 0.6 and 0.8 were used for spring constant factor and Laplace node relaxation, respectively. Remeshing method included local cell, local face and region face; also, the sizing function is turned on.
In terms of the flow surface (XY-direction), 4478 elements were employed (including the refined area) for both 2D and 3D simulation; in the depth of the flow (Z-direction), 150 elements were applied for 3D simulation. In total, 671,700 elements, which balanced the accuracy and computational cost, were utilized in the flow zone (3D simulation). The number of the elements for the fluid zone was decided based on convergence study, i.e., reducing the number might lead to the loss of accuracy while increasing the number could not change the result obviously.
For the steel riser, the traditional modeling method was utilized with isotropic steel. For the composite risers (Risers 2 and 3), two different methods, effective modulus method and layered-structure method, were established. In effective modulus method, composite riser is regarded as an anisotropic monolithic structure rather than layered-structure using Equations (1)-(6), i.e., only the effective properties for the entire riser in Table 3 is required. In the layered-structure method, the properties for the liner and each composite lamina (Table 2) are required as well as the lamination configuration of the riser (Table 1) to create the real structure of the composite riser. Comparatively speaking, EMM is easier during the modeling process and needs less computational resource to obtain the global results, while LSM can be used to get the response of each composite lamina at the expense of larger computational resource and longer calculation time. The risers are modeled by Solid185; 240 elements in length direction and 20 elements in hoop direction were employed. In total, 4800 elements were utilized in the riser based on convergence study. Compared to the element number in the fluid zone, the number of elements for risers would not lead to significant difference in the requirement of computational resource and calculation time. All the risers were fixed at the bottom and simple supported at the top. Top tension force was also employed at the top end of the risers, which equals to 1.5 times the weight for the steel riser and 2 times the weight for the composite riser (Table 4).
In the two-way FSI, the riser deformation, the flow data variation and the influence between them were considered. Figure 3 shows the FE model of the flow and riser.

Results and Comparison
In this section, VIV characteristics of three different risers are compared. The effects of 2D and 3D models and fluid-structure interaction (FSI) have been considered, as well as the effect of the two different modeling methods for composite risers.

Properties for Different Risers
Based on the risers' geometries, material properties, buoyancy and top tension forces, pre-stress modal analysis is conducted to obtain the natural frequencies and modal shapes for the risers, as shown in Figure 4 and Table 5.
When the vortex shedding frequency becomes close to the natural frequency of the riser, the phenomenon of "lock-in" happens, which results in large and damaging vibrations leading to the fatigue failure of the riser. The vortex shedding frequency can be calculated by Equation (9) which includes Strouhal number S t , the outside diameter of riser D and the flow velocity U. According to the relationship between Re and S t , the Strouhal number S t for all the risers is 0.21.
For a riser, the vortex shedding frequency can also be obtained by a fast Fourier transform (FFT) of the lift coefficient C L ( Figure 5). These results are compared and summarized in Table 5.  In addition, the possibility of a riser's "lock-in" using reduced velocity U r was estimated using Equation (10). When U r ranges 4-8 [5], "lock-in" occurs. The U r of all the risers are listed in Table 5.
From the results listed in Table 5, it can be seen the natural frequencies are far from their vortex shedding frequencies and the reduced velocity U r are not in the range of 4-8. Hence, for all the risers, "lock-in" is unlikely to happen.

Effects of 2D and 3D Models
Figure 6a-c presents the vorticities for Risers 1, 2 and 3 with 2D models while Figure 6d-f illustrates the vorticities for Risers 1, 2 and 3 with 3D models. As can be seen in Figure 6, the vortex shedding modes for all three risers are 2S (two single vortices per cycle). However, the vorticities from the 2D model are more stable as in the ideal condition.  In Figure 7, the value and stability of C L and C D with 2D and 3D models are quite different, although the trend for C L and C D with 2D and 3D simulation are similar. The curves of C L and C D for all the risers with 2D model are more stable and the values are significantly larger than those with 3D models. It also can be seen that the "beat" phenomenon of C L occurs in the 3D simulation, i.e., C L increases and decreases periodically.
In the 2D simulation, the size in depth direction (Z-direction) is omitted; therefore, the effects of depth of the flow field on the VIV are not considered in 2D simulation. However, in any real flow, 3D effects are included and the 3D vorticity profiles of Risers 1, 2 and 3 are presented in Figure 8a-c, respectively. It can be seen in Figure 8 that the VIV of the risers have obvious 3D feature.    Figure 9a-c presents the C L for Risers 1, 2 and 3 with/without FSI while Figure 9d-f illustrates the C D for Risers 1, 2 and 3 with/without FSI.

Effects of the Models with/without FSI
As can be seen in Figure 9, FSI affects all the risers' C D and C L . In general, the FSI has greater effect on C L than C D . From the graph of C D , the curve of C D ascends initially, then descends, after that rises again, and at last tends to be stable. Through the development of vortex shedding, it is found that the reason for the change of the C D curve is that, due to the impingement of the water flow, the riser is subjected to a large resistance in the forward direction of the X-axis for 0.4 s (Figure 10a). After 0.4 s, a pair of symmetrical vortices begins to form, and these symmetrical vortices develop steadily until 8 s (Figure 10b). Because of the existence of the vortices behind the riser, the riser is subjected to a negative force along the X-axis. Hence, the resistance of the riser along the X-axis is reduced and the C D curve decreases. After 8 s, the symmetric tail vortices begin to fall off gradually, and the C D begins to increase until 14 s (Figure 10c). When one of the symmetric vortices completely falls off, the C D returns to a relatively large value.

Results for 3D VIV with FSI Using Different Modeling Methods
The 3D CFD simulation with FSI is adopted for the VIV of all the risers. For Riser 1 (steel riser), monolithic model is utilized and for Risers 2 and 3 (composite risers), two different modeling methods, i.e., effective modulus method (EMM) and layered-structure method (LSM) are employed. The displacement time history of all three types of risers are plotted in Figure 11 for three typical location, which are 1.2 m, 6.2 m and 11.2 m below the water surface. More specifically, Figure 11a     In Figures 11 and 12, it is found that: (1) For the same riser: (i) different modeling methods lead to similar results for riser's displacement time history; (ii) the deformation in the middle part of the riser is larger than that in both ends, and the deformation close to the fixed end is the minimum; (iii) the deformation in the flow-direction is relatively stable at a specific value while the deformation in the cross-flow direction increases and decreases periodically; and (iv) the vibration in cross-flow direction is more obvious than that in flow direction. (2) For the different risers, Riser 1's displacement is less than Riser 2's displacement, which is less than Riser 3's displacement, due to their different bending modulus and tension forces.
The time history of displacement for all the risers in flow direction is quite different from that in cross-flow direction, that is, the flow-direction displacement will maintain at a relatively stable value and vibrate slightly after the initial growth. Figure 13 illustrate the maximum displacement of the risers in flow-direction for all the risers with different modeling methods. Figure 13a shows the maximum flow-direction displacement for Riser 1; Figure 13b,c presents the maximum flow-direction displacement for Riser 2 with EMM and LSM, respectively; and Figure 13d,e illustrates the maximum flow-direction displacement for Riser 3 with EMM and LSM, respectively. The flow-direction displacements are various for different risers: Riser 1's displacement < Riser 2's displacement < Riser 3's displacement. For the same riser, EMM and LSM lead to the similar flow-direction displacement. Figures 14 and 15 show the maximum von Mises stress and the stress distribution diagram for the Risers 1, 2 (EMM) and 3 (EMM).  In Figures 14 and 15, it is obvious that the maximum von Mises stress follows: Riser 1 > Riser 3 > Riser 2. This is not totally dependent on the risers' modulus. Meanwhile, it can be seen that the maximum von Mises stress occurs at the top end of the steel riser (Riser 1) due to large top tension force required to resist the extreme weight of the steel riser; on the contrary, the maximum von Mises stress of the composite risers (Risers 2 and 3) take place in the middle/at the fixed end of the risers due to the large bending moments.
The overall performance and characteristics of composite risers can be obtained using EMM; however, it is impossible to analyze the stress distribution and failure criterion in every composite lamina by EMM. Therefore, LSM is employed to solve this problem. Figure 16 presents the maximum stress distribution including stress in fiber direction (S 1 ), stress in transverse direction (S 2 ) and stress in shear (S 12 ) in every composite lamina for Risers 2 and 3. It is shown in Figure 16 that the stress distribution in each composite lamina in composite Risers 2 and 3 are totally different, since their designed lay-ups are diverse, that is, composite Riser 2 only contains fiber reinforcements in hoop and axial directions while composite Riser 3 includes fiber reinforcements in hoop, axial and inclined directions. Both composite Risers 2 and 3 are dominated by bending deformation under the combination of VIV, buoyancy, top-tension force and gravity. Generally, the stress in each composite lamina is larger in composite Riser 3 than in composite Riser 2, which means tailor-designed composite Riser 3 can take advantage of the fiber strength more efficiently.
Comparing the stress distribution in every composite lamina for composite Risers 2 and 3, it is found that: (1) Stress in fiber direction (S 1 ): S 1 is smaller in hoop reinforced fiber layers than that in axial reinforced fiber layers for composite Riser 2 while S 1 is largest in axial reinforced fiber layers, medium in inclined reinforced fiber layers and smallest in hoop reinforced fiber layers for composite Riser 3. (2) Stress in transverse direction (S 2 ): S 2 is larger in hoop reinforced fiber layers than that in axial reinforced fiber layers for composite Riser 2, while S 2 is largest in hoop reinforced fiber layers, medium in inclined reinforced fiber layers and smallest in axial reinforced fiber layers for composite Riser 3. (3) Stress in shear (S 12 ): Riser 2 only contains orthorhombic fiber reinforcements, and therefore almost no shear stress S 12 is produced; contrarily, composite Riser 3 includes shear stress S 12 in the inclined composite layers.

Conclusions
In this paper, the VIV responses of three different risers, i.e., conventional steel riser (Riser 1), composite riser with orthogonal reinforcements (Riser 2) and tailored composite riser (Riser 3) are investigated and compared considering the effects of 2D and 3D models, FSI and different modeling methods for the composite risers. The results show that: (1) Both 2D and 3D models lead to 2S vortex shedding modes for all three risers and VIV of a riser has obvious 3D characteristics and is affected by the depth of the flow field as well as the effect of FSI. As a result, a 3D model with FSI for the VIV is more accurate and realistic. (2) The amplitude of the cross-flow vibration is much larger than that of the flow-direction vibration, and the cross-flow vibration has a distinct periodicity while the flow-direction vibration tends to be at a relatively stable value. The deformation of the three risers are different due to their geometries, material properties, top-tension force, constraints, etc.