Stress wave propagation analysis on vortex-induced vibration of marine risers

To analyze the stress wave propagation associated with the vortex-induced vibration (VIV) of a marine riser, this paper employed a multi-signal complex exponential method. This method is an extension of the classical Prony’s method which decomposes a complicated signal into a number of complex exponential components. Because the proposed method processes multiple signals simultaneously, it can estimate the “global” dominating frequencies (poles) shared by those signals. The complex amplitude (residues) corresponding to the estimated frequencies for those signals is also obtained in the process. As the signals were collected at different locations along the axial direction of a marine riser, the phenomena of the stress wave propagation could be analyzed through the obtained residues of those signals. The Norwegian Deepwater Program (NDP) high mode test data were utilized in the numerical studies, including data sets in both the in-line (IL) and cross-flow (CF) directions. It was found that the most dominant component in the IL direction has its stress wave propagation along the riser being dominated by a standing wave, while that in the CF direction dominated by a traveling wave.


Introduction
Vortex-induced vibration (VIV) of long flexible structures such as marine risers exposed to the marine environment is ubiquitous in the offshore industry. Understanding the vibration behavior of marine risers excited by vortex shedding is of great importance to the offshore structure designs (Trim et al., 2005;Song et al., 2011;Wu et al., 2012). Analyzing the data of field experiment on a long flexible cylinder densely instrumented with fiber optic strain gauges, (Vandiver et al., 2009) first revealed a previously unknown phenomenon: the dominance of traveling waverather than standing wave-on riser's strain/stress response.
A traditional way to analyze the stress wave propagation associated with the fundamental vortex shedding frequency is often based on using a colored image plot (Wu et al., 2012;Vandiver et al., 2009;Ge et al., 2009). The processed data used in the image plot are usually extracted by a two-step signal processing procedure from the raw sensor signals: (1) carrying out the fast Fourier transform (FFT) of the time series at each sensor location one by one, and (2) filtering the signals at the dominant fundamental frequency which usually involves a relatively subjective peak-picking process. From these processed data, the image plot was generated with its horizontal axis being the passage of time, the vertical axis being the non-dimensional axial location on the riser, and the color being the instantaneous amplitude of the filtered strain signals.
This paper proposes a multi-signal complex exponential method to investigate the stress wave propagation. By processing all signals simultaneously, the proposed method-an extension of the classical Prony's method-estimates the fundamental vortex shedding frequency through an eigen analysis on the raw signals. The complex residue associated with the fundamental vortex shedding frequency for each signal is also obtained in the process. Sequentially, the phenomena of stress wave propagation can be studied by ana-lyzing the obtained complex residues. In comparison with the traditional Fourier-based methods which always suffer from the problems related to frequency leakage and resolution, the proposed method completely avoids those problems. Another advantage of the proposed method is that it requires only short duration signals to conduct the stress wave propagation analysis (Hu et al., 2013(Hu et al., , 2014. In the numerical studies, measured acceleration signals in both the in-line (IL) direction and the cross-flow (CF) direction from the Norwegian Deepwater Program (NDP) high mode test will be analyzed (Trim et al., 2005).

Preliminaries
2.1 p-th order constant-coefficient linear differential equation a n A p-th-order homogenous linear ordinary differential equation (ODE) with constant real-valued coefficients can be written as: p X n=0 a n y (n) (t) = 0; for 0 6 t < T; (1) . Substituting into Eq. (1) yields the characteristic equation: p X n =0 a n n = 0: (2) y = e t y 1 = e 1 t ; ¢ ¢ ¢ ; y p = e p t If λ is a root of Eq.
(2), then is a solution of Eq. (1). And if all the p roots λ 1 …λ n of Eq. (2) are different, then the p solutions constitute a basis for Eq. (1). Thus the corresponding general solution of Eq. (1) is a linear combination of those exponentials, i.e., n e n t ; 0 6 t < T: (3) In Eq. (3), when is a real-valued signal, must either be real numbers or occur in complex conjugate pairs. Let, and is defined as the damping factor in s -1 and is the frequency in radians. The coefficients corresponding to complex exponents must also appear in complex conjugate pairs. Let , then is the amplitude and is the sinusoidal initial phase in radians associated with . Taking the Laplace transform of Eq. (3) yields Usually and are also referred to as the residues and poles associated with the function .
2.2 p-th order constant-coefficient linear difference equation y(t) y k´y (t k ) t k = k¢t k = 0; 1; : : : ; N ¡ 1 ¢t In digital signal processing problems, one is almost always interested in a time signal only for the discrete values. Denote , , for , in which =time step increment, and N=the num-ber of time steps. In parallel to the ordinary differential equation shown at Eq. (1), a p-th-order homogenous linear, constant-coefficient difference equation can be written as: In the above equation, if both b 0 and b p are different from zero, the positive integer p is called the order of the equation.
To find the solutions of Eq. (5), one might try, as with the analogous differential equation, the substitution y k = e t k = (e ¢t ) k : However, it is more convenient to write y k = z k ; (7) where . Substituting Eq. (7) into Eq. (5), and dividing out , one can obtain the characteristic equation of the difference equation Eq. (5) as: If the roots of Eq. (8) are distinct, then n z n = e n ¢t n = ln(z n )=¢t As Eq. (9) is the discrete counterpart of Eq. (3), the two equations must have the same residues . Since , .
3 Extended Prony's method for multiple signals y k n n y k k = 0; 1; : : : ; N ¡ 1 Prony's method is a technique for fitting the exponential model Eq. (3) (or Eq. (9)) using equally spaced data . In other words, Prony's method is to determine and based on the set , . While the original Prony's method presented a method of exactly fitting as many purely damped exponentials as needed to fit the available data points, the modern version of Prony's method generalizes to damped sinusoidal models and also makes use of the least squares analysis to approximately fit an exponential model for cases where there are more data points than needed to fit to the assumed number of exponential terms.

y k
The exponential model Eq. (9) for Prony's method readily suggests that the corresponding difference equation model is Eq. (5). In fact, the first step of Prony's method is to use the equally spaced data to estimate the constant coefficients of the linear difference equation. Once the difference equation Eq. (5) is determined, one can proceed the calculation for poles and residues associated with the difference equation. In summary, three sequential steps, using Eqs. (5), (8) and (9), respectively, are implemented in Prony's method.
Step 1: Without losing the generality, let . By repeatedly using Eq. (5), one determines the linear prediction parameters that fit the given discrete signal. In matrix form, it is written as: y 0 y 1 ¢ ¢ ¢ y p¡1 y 1 y 2 ¢ ¢ ¢ y p : : : : : : : : : : : : y p y p+1 : : : z n ;n = 1; : : : ; p Step 2: With the linear prediction coefficients b obtained from Step 1, the roots (denoted as ) of the characteristic equation Eq. (8) are computed. z n (n = 1; : : : ; p) Step 3: Once have been computed, the third step involves the solution of the second set of linear equations to yield the estimates of the exponential amplitude and sinusoidal initial phase. Based on Eq. (9), a matrix form is written as: 2 : : : : : : : : Recall that the roots of Eq.
(2), namely , are related to those of Eq. (8) through . Appropriate linear least squares procedures for the first and third steps of the three-step Prony's method has sometimes been called the extended Prony's method.

Multiple signals decomposition
Prony's method can be further extended to perform the decomposition for multiple signals (Trudnowski et al., 1999). Denote the column vector for the signal data at where is the number of signals. An intended feature is that all signals share the same model-a p-th-order difference equation as Eq. (5). Thus, one can write With this feature, a consistent set of "global" poles (natural frequencies and damping factors) is among all signals.
The extension of Prony's method from the single signal to multiple signals is quite straightforward for the three steps: Step 1: By repeatedly using Eq. (13), one determines the linear prediction parameters (noting ) from , :  : (14) z n ; n = 1; : : : ; p Step 2: With the linear prediction coefficients b obtained from Step 1, the roots (denoted as ) of the characteristic equation Eq. (8) are computed. z n (n = 1; : : : ; p) Step 3: Once have been computed, the third step involves the computation for the residues of each signal. Denote for the n-th residue vector for all signals. A straight forward extension of Eq. (12) leads to: 2 : : : : : : : : Since the second step of Prony's method is an ill-conditioned problem and round-off errors must exist for the linear prediction parameters b computed in the first step, the estimation of in the second step of Prony's method can have significant error.

Numerical implementation
To avoid dealing with an ill-conditioned problem on determining the roots of a high order polynomial encountered in Prony's method, an alternative approach which begins with a first-order matrix difference equation-a state-space model-to replace a high order difference equation used in Prony's method has been advocated (Hu et al., 2013). The state-space model approach also involves three steps. While its third step is identical to that of Prony's method, the first two steps of the proposed method are: (1) obtaining a realization for the state matrix of first-order matrix difference equation from the Hankel matrices composed by the signal data; and (2) computing the eigenvalues of the realization matrix which are the global poles of the signals. One obvious advantage of the state-space model approach is to completely avoid the ill-conditioned problem of solving the zeros of a high order polynomial.

Experimental data studies
The VIV signals investigated are from the Norwegian Deepwater Program (NDP) high mode test, which was a high length-to-diameter ratio (L/D) riser hydrodynamics testing program: a riser model of L =38 m and D = 0.027 m. For each test run, the riser was setup in a way, and towed at a constant speed, to simulate the riser subjected to either uniform or sheared current. Measurements of micro bending strain and acceleration in the in-line (IL) and cross-flow (CF) directions along the riser were taken. The interested reader can refer to Trim et al. (2005) for the details of the test.
In the present numerical study, eight acceleration data sets in the IL direction and eight in the CF direction from Test 2350, for which the sheared current profiles are corresponding to the towing speed 0.7 m/s, will be analyzed. The eight acceleration signals were collected at z/L = 0.1314, 0.2404, 0.3381, 0.4370, 0.5555, 0.6400, 0.7735 and 0.8907, respectively, where the vertical coordinate z is positive upward and with the origin at the bottom of the riser. To show typical measured signals, what displayed in Fig. 1 are the acceleration signals at z/L = 0.8907 in both the IL and CF directions. It seems that the segments between 20 s and 30 s of both signals are stationary.
In the following study, only segments between 20 s and 30 s are utilized to implement the proposed method. The original sampling rate of the measured signals has been 1200 Hz. This sampling rate is unnecessarily high while implementing the proposed method. For reducing the data size in order to save the computational time, the original signals have been decimated with the new sampling rate of 300 Hz. As a result, each 10-second signal to be analyzed contains 3000 data points. Fig. 2 shows the acceleration signals in the IL and CF directions at z/L = 0.8907 between 20 s and 30 s with sampling rate of 300 Hz. Signals in both the IL and CF directions will be investigated, but they will be analyzed separately.

In-line signals analysis
All the eight signals between 20 s and 30 s-each is a segment containing 3000 data points-in the IL direction are processed simultaneously. In the numerical implementation, a first-order matrix difference equation is employed to replace a high order difference equation for avoiding an illconditioned problem on determining the roots of a high order polynomial. The detailed procedure of employing a first-order matrix difference equation has been described by Hu et al. ( 2013Hu et al. ( , 2014. The numerical implementation of the multiple signals decomposition based on a state-space model approach for the present data sets involves the following three steps: (1) obtain a realization for the state matrix of first-order matrix difference equation with the model order , (2) compute 40 eigenvalues , which are 20 complex conjugate pairs, of the realization matrix, and obtain the corresponding poles , and (3) compute the residues associated with the eight signals by Eq. (15).
Although 20 complex exponential components are obtained from the current analysis, only the most dominant component is presented below to show its stress wave propagation phenomenon. For a damped sinusoidal signal, the variance of a signal is a good index to represent the vibration level of the signal because variance reflects energy in physics. In this paper, the most dominant component has been judged based on the largest variance among all components (Hu et al., 2014). Denoting the pole associated with the most dominant component in the IL direction by , one obtains . Letting the corresponding frequency be defined by and damping ratio by , one shows Hz and , respectively. Table 1 listed the corresponding complex residues associated with the eight signals, , . Furthermore, using the polar coordinate for the complex residues , the corresponding amplitudes and phase angles are also listed in Table 1.   where the superscript * indicates the complex conjugate. Plotted in Fig. 3 are the reconstructed signals , , based on the parameters listed in Table 1, compared with the original signals (in dotted line) shown in the background. Overall, the reconstructed dominant components appear to match well with their original signals.
The 8 complex residues in Table 1 normalized by the maximum of are plotted in Fig. 4a, and a linear interpolation method has been applied to generate 800 complex residues at the equally-spaced locations along the riser. Shown in Fig. 4b and Fig. 4c are the corresponding amplitudes and phase angles. Notice that the interpolation of the amplitudes and the interpolation of the phase angles do not follow a linear way while the interpolated amplitudes and phase angles have been computed directly from the interpolated complex residues.
The physics of the stress wave propagations can be interpreted easily from Fig. 4. As shown in Fig. 4a, when the imaginary part of a complex residue vector diminishes, it implies that the vibration mode becomes a real-valued mode. In the physics of wave propagation, a real-valued mode indicates a standing wave. In contrast, if a residue vector also contains a significant imaginary part, it reflects a traveling wave. In Fig. 4b, the amplitudes provides the physical information of vibration level. Fig. 4b suggests that the dominant mode of the riser experienced a stronger vibration at the top portion of the riser. As shown in Fig. 4c, when the phase angles at various locations are mainly either in phase or 180° out of phase, it suggests a standing wave phenomenon. It is recommended that the amplitudes and phase angles be read together. From Fig. 4b and Fig. 4c, one concludes that the top portion of riser vibrates more severely in a standing wave form while the bottom portion vibrates mildly in a less standing wave form. I The above residue vector analysis represents a novel way to interpret the stress wave propagation associated with the dominant VIV component. In contrast, the traditional way to observe the wave propagation phenomenon has been carried out by using an image plot. The image plot is a three-dimensional representation, where the horizontal axis is the passage of time, the vertical axis is the non-dimensional axial location on the riser, and the color is the instantaneous amplitude of signals. For completeness, an image plot based on the proposed method is also conducted below. Firstly, the reconstructed signals associated with the dominant frequency component are obtained based on the dominant pole and the corresponding residues, including those from the linear interpolation. The image plot shown in Fig. 5 has been obtained based on 800 equally-spaced reconstructed signals. It is clearly shown in Fig. 5 that the local maximum and minimum at various locations occur almost at the same time, which indicates that the phase angles of the reconstructed signals are nearly in phase or 180° out of ph- Similar to the procedure and tasks carried out for the signals in the IL direction, all 8 signals between 20 s and 30 s in the CF direction are processed simultaneously in the remaining study. The following presentation will avoid repeating the same description, but highlight the different features observed in the CF direction. Denote the pole of the most dominant component in the CF direction by , and the corresponding frequency and damping ratio . Taking the model order p=40, one obtains , Hz and . The corresponding complex residues , amplitudes and phase angles for , in the CF direction are listed in Table 2.  Table 2, as well as the amplitudes and phase angles . Basically, Fig. 4 and Fig. 7 exhibit very different patterns. As shown in Fig. 7a, both the real and imaginary parts of the normalized complex residue vector are equally significant. In other words, the vibration mode is a truly complex-valued mode which implies a travelling wave in physics. From Fig. 7b, it suggests that the dominant mode in the CF direction has a slightly stronger vibration level near the bottom of the riser. As shown in Fig. 7c, the phase angles vary almost linearly along the riser which also suggests a traveling wave. Note that a phase angle π is equal to a phase angle -π, thus a sudden jump of the phase angle in Fig. 7c from π to -π really means a very small change in phase angle.
The image plot of the dominant component in the CF direction based on the reconstructed signals for seconds is shown in Fig. 8. Fig. 8 clearly reveals the traveling wave phenomenon which is completely different from the characteristics shown in Fig. 5. From the slope of the di-C  , and (c) initial phase angle , associated with . C Fig. 8. Image plot associated for 21≤t≤22.
f C agonal line associated with the peaks or valleys in Fig. 8, the speed of the traveling wave can be roughly estimated. The speed of the traveling wave can also be estimated from the corresponding wave length, which can be computed from the phase angles at two locations in Fig. 7c, together with the already known wave frequency .

Comparison between IL and CF dominant components
From the estimated = 6.7475 and = 3.3716, the ratio is 2.0013. This ratio is in a perfect agreement with the traditional wisdom that the dominant frequency in the IL direction is twice that in the CF direction. Both the estimated damping ratios and are near zero, which indicates that both dominant components are steady without noticeable decay. An interesting observation from this study is that the wave propagations in the IL and CF directions are fundamentally different. While a traveling wave is dominant in the CF direction, a standing wave is dominant in the IL direction.

Concluding remarks
With multiple measured acceleration signals, this paper developed an efficient method, a multi-signal complex exponential method, to analyze the stress wave propagation associated with the vortex-induced vibration of marine risers. The novelty of the multi-signal complex exponential method is to simultaneously process many signals to determine the dominant "global" poles (i.e. dominant frequencies and the associated damping factors) among those signals. In the process, the corresponding residues (i.e. amplitudes and initial phase angles) of each individual signal are also determined. From the extracted residues, the stress wave propagation phenomena could be interpreted easily.
In the numerical studies, the Norwegian Deepwater Program (NDP) high mode test data were utilized, including eight acceleration data sets in the in-line (IL) direction and 8 in the cross-flow (CF) direction. Agreeing perfectly with the traditional wisdom, the numerical studies found that the dominant frequency in the IL direction is twice as large as the dominant frequency in the IL direction. Another interesting observation about the dominant components in the IL and CF directions is that the stress wave propagation along the riser is dominated by the standing wave in the IL direction, but dominated by the traveling wave in the CF direction.