Slender structures immersed in a cross flow can experience vibrations induced by vortex shedding (VIV), which cause fatigue damage and other problems. VIV models that are used in structural design today tend to assume harmonic oscillations in some way or other. A time domain model would allow to capture the chaotic nature of VIV and to model interactions with other loads and nonlinearities. Such a model was developed in the present work: for each cross section, recent velocity history is compressed using Laguerre polynomials. The compressed information is used to enter an interpolation function to predict the instantaneous force, allowing to step the dynamic analysis. An offshore riser was modeled in this way: some analyses provided an unusually fine level of realism, while in other analyses, the riser fell into an unphysical pattern of vibration. It is concluded that the concept is promising, yet that more work is needed to understand orbit stability and related issues, in order to produce an engineering tool.

1. Introduction

1.1. Relevance

Vortex-induced vibration (VIV) is a vibration of a flexible structure that occurs when a fluid flowing around the structure sheds vortices at near-regular intervals, locked with the structure's own vibration. VIV is a major concern in the offshore oil industry in particular, where marine currents can cause slender structures like pipelines, risers, umbilicals, and cables to vibrate, inducing fatigue damage. While design tools are available, they are still improvable. Today, VIV is still actively studied.

A brief classification of existing VIV models is presented in the following. The classification is biased in the sense that it aims at comparing existing models with the model proposed here. More comprehensive overviews of existing models can be found in [1, 2].

1.2. Detailed Wake Models

In this group of models, the details of the wake flow behind the structure are resolved, to various levels of detail, by using various techniques of computational fluid dynamic. Such models can be coupled to a structural model, which typically uses beam elements. Because the water behaves in a strongly nonlinear fashion, such models operate in the time domain. While all models use some sort of “strip theory”, computing the flow at a limited set of points along the riser, the computation in a given strip may allow 3D turbulence or limit it to 2D, the later being now recognized as unsuitable.

Orcina's vortex tracking method is based on [3]: the vorticity is assumed concentrated in a curve that is convected by itself and the incoming current.

Deepflow [4], the USP code [5], and VIVIC [6] use discrete vortex solutions in a series of planes along the riser.

ACUSOLVE is a general purpose computational fluid dynamic software. It has been used for VIV modeling [7]. Such an approach is computationally intensive.

1.3. Simplified Harmonic Models

The common denominator of the models in this group is that they operate by characterizing the oscillation, at any given point along the riser, by frequency, amplitude, and possibly phase difference between oscillations in two orthogonal directions (in-line and cross-flow). These values are used to enter tables that yield excitation and added mass coefficients. Such tables first appeared, to the author’s knowledge, in [8]. Models differ widely on how they make use of the above coefficients to estimate a solution.

Because they assume harmonic vibration, the models in this class tend to share the same approach to similitude: the reduced frequency is computed as the ratio of the time it takes a particle in the undisturbed flow to travel one cylinder diameter, divided by the oscillation period.

VIVA [9] models the response of the structure as a superposition of amplitude-modulated traveling waves. Again, force coefficients are obtained from tables.

In [10], a time domain solution is used, in which, at any step and point along the cable, the recent computed velocity is approximated by a harmonic function of time. This is used to enter the above-mentioned tables. The instantaneous value of the force is then computed from the hydrodynamic coefficients, allowing to pursue the time domain integration.

SHEAR7 [11] starts from a modal analysis of the structure. Modes are then examined for their susceptibility to lock-in. Tables are used for lift coefficients.

In VIVANA [12, 13], an iteration scheme is used to arrive at a harmonic solution for the whole structure. When relevant, the solution can be a superposition of such oscillation “modes”.

1.4. Simplified Nonharmonic Models

Models in this group forgo a detailed description of the flow in the wake, replacing it by a highly simplified nonlinear model with very few degrees of freedom. The model is repeated at several points along the oscillating structure. The analysis operates in the time domain, and the response of the structure is typically computed using finite element analysis. The challenge in such models is to capture the influence of structural motions on the wake, and of the wake forces on the structure, in a compact model.

Several models make use of simple nonlinear oscillators to represent the self-exciting and self-limiting nature of VIV response: A single degree of freedom van der Pol oscillator has been used in several models [1417]. Orcina also uses a wake oscillator with few degrees of freedom [18].

1.5. Discussion

The drawback of detailed wake models is that for the relevant Reynolds numbers, they are computationally very demanding. Hence they do not really offer a practical option in structural design and design verification, where extensive computations are needed to adequately sample the statistics of currents and other operating conditions that the structure is likely to encounter.

Despite the fact that they currently provide the most used tools in VIV design, simplified harmonic models have several drawbacks. Most of them do not operate in the time domain, which makes it difficult to account for structural nonlinearities. The models characterize the oscillation of a cross section by amplitude and frequency, which is not adequate to describe more general types of motions.

Time domain models based on nonlinear oscillators have been proved to be able to reproduce some aspects of VIV behavior, but so far seem limited in their ability to capture the details of the response in a range of current conditions.

A good time domain, nonharmonic, simplified model, if it existed, would open new possibilities, compared to harmonic models: (1)study of VIV on nonlinear structures, for example, studying the damping effect of seafloor interaction in a steel riser or using a hysteretic cross section model for VIV on flexible pipes, (2)accounting for VIV caused by unsteady water flows, in particular by waves or vessel motions, (3)accounting for the increase in drag at wave frequency due to VIV, (4)accounting for the superposition of wave-frequency and VIV-frequency stresses in fatigue analysis, (5)accounting for the asymmetry of oscillation patterns in the vicinity of, for example, a seafloor.

1.6. Objective

The objective of the work reported here is to demonstrate the viability of a local, deterministic, time-domain force model for VIV on slender bodies with cylindrical cross sections.

The model is to treat in-line and cross flow vibrations jointly.

It is to characterize the recent history of velocity of the cross section relative to the surrounding fluid without making a harmonic assumption. The characterization is to be used to enter a “table”, necessarily more complex than those used under harmonic assumption, to predict the instantaneous value of the hydrodynamic force. The model is to handle external steady or unsteady water currents.

Like many other models discussed above, this force model is to be used at each Gauss point of the dynamic finite element (FE) model of a slender structure. The model is to be used within a time domain analysis (e.g., Newmark-β time integration with Newton-Raphson iteration).

Hence the FE model resembles that commonly used in a slender structure analysis, with degrees of freedom for the structure, and none for the surrounding fluid. In other words, the proposed model takes the place usually held in software by the Morison model for wave induced loads.

2. Model Outline

2.1. Postulate

The present work hinges on the following postulate. The force exerted by the surrounding fluid on a section of the slender structure is completely determined by the recent histories at that section of the velocities of the structure and of the undisturbed fluid. Several points in this sentence are worthy of discussion.

The “force” includes the components usually distributed into added mass, excitation forces, drag, lift, and so forth.

That the force “at a section of the slender structure” is determined by the history “at that section” implies a “strip theory” in which it is excluded that motions of the structure at a point A cause disturbances in the fluid that affect the force at point B away from A. In other words, it is assumed that there is no significant transmission of information in the axial direction within the water (as opposed to within the slender structure). This would be proved wrong if it turned out that unstable phenomena, like boundary layer shedding, although transmitting little energy along the structure, transmit information that steers how local hydrodynamic energy is channeled at a given point along the structure.

That the force should be “completely determined” implies that the behavior of the structure is deterministic. This does not contradict the observation of hysteretic response of short cylinders mounted on elastic support. Uniqueness of forces for a given position does not imply uniqueness of static equilibrium. Neither does “completely determined” contradict the observation of irregular and unpredictable responses to VIV: nonlinear dynamic systems can have a chaotic behavior. Still, complete determinism is provably wrong, since a short vertical cylinder dragged at uniform speed through water will experience oscillating lift forces. At any given moment, there is nothing in the history of (constant) velocity that allows to predict whether the lift is left or right. So the present work is based on the bet that ignoring such “bifurcations” still leaves us with a useful model.

“Recent” can be defined as anything between the present time and a few times 𝑡𝑤, where the value of 𝑡𝑤 still is an object of debate. 𝑡𝑤 is likely to be case dependent. The current will transport (convect) away vortices so that they quickly loose significance. The time 𝑡𝑤 should then be of the order of 𝐷/𝑈 where 𝐷 is the cross section diameter and 𝑈 the current velocity. In contrast, if the cylinder is oscillating in still water, it will be traveling in its own wake, and 𝑡𝑤 should be related to the rate of diffusion and/or viscous dissipation of vortices, which is likely to result in much higher values of 𝑡𝑤. Tests on periodic forced motion of short cylinders sometimes show a slow drift of the forces (over as many as ten periods). In contrast, force decay tests for a cylinder stopped after oscillations at zero mean velocity, point towards a fraction of a period. In the present work, the idea is to choose an upper bound for 𝑡𝑤, after adequate scaling (cf. Section 3.2).

The “velocity histories” are what count. Accelerations would not do because for example, zero acceleration can correspond to different speeds and hence different forces. On the other hand, the force on a cylinder will not be affected by a uniform translation of its whole trajectory, so a history of positions contains irrelevant information.

In the remainder of this text, the word “trajectory” will be given a very specific meaning. The trajectory is defined as the recent history of the velocity vector of the cylinder relative to the undisturbed surrounding fluid.

2.2. Restrictions

In the present phase of research, the following restrictions are introduced, in order to achieve some simplification of the task. The outer cross section of the slender structure is assumed perfectly circular and smooth. The surrounding fluid is assumed to be infinite, excluding the presence of sea floor, free surface, or neighboring risers. Only fluid flows perpendicular to the cylinder at any point are considered.

2.3. Input and Output

As stated earlier, the VIV model being developed here replaces the Morison model for wave-induced loads. The VIV model is called at each step and iteration, and at each Gauss point or node of each element.

The model is to receive as input: (1)the diameter of the cylinder, the viscosity, and density of the surrounding fluid, (2)the instantaneous velocity of the cross section relative to the undisturbed fluid, (3)the instantaneous velocity and acceleration of the local undisturbed fluid, in a Galilean reference system.

The model uses velocity information stored from previous steps. On this basis, the model produces as output:(1)the vector of hydrodynamic forces per unit length, acting on the cylinder, (2)the matrix containing the derivative of the above with respect to instantaneous values of the cylinder’s behavior.

Gauss integration is then used to compute a consistent load vector and partial derivative matrices (damping, stiffness, and mass) for each element. Note that these element matrices are likely to vary significantly over each VIV oscillation “period”—in contrast to added mass or damping matrices, deemed to be constant over a long time in semiempirical VIV models. The connection of the force model to the finite element analysis is discussed in Section 7.

2.4. Algorithmic Steps

Only the local VIV model is described here, not the whole FE analysis. (1)The relative velocity of the cylinder relative to water (thereafter: “velocity”) is computed. (2)The velocity is scaled (Reynolds scaling) to that the cylinder diameter is the unit of distance (Section 3.2). (3)The trajectory (again: the recent histories of both 𝑥 and 𝑦 components of velocity) is compressed into a small number of “Laguerre coefficients”. This compression is such that it provides detailed information over the recent past and increasingly coarse information for the more distant past (Section 4). (4)The Laguerre coefficients are used to enter an interpolation function (a feed-forward neural network with some specifically tailored properties) which returns 𝑥 and 𝑦 components of hydrodynamic force (Section 5). The fitting of the interpolation function is discussed in Section 8.1. (5)The force is scaled back to the relevant diameter (Section 3.2). (6)The Froude-Krylov forces, which depend on the acceleration of the undisturbed flow, are added (Section 3.1)

The identification of nonlinear systems using a bank of orthogonal filters (including Laguerre filter) to generate multiple signals from a single one, and then using the multiple signals to enter a nonlinear, memory-less function, was introduced by Wiener [19] (Wiener-Laguerre filtering). In the present work, a base of Laguerre polynomials is used, in contrast to Laguerre functions introduced by Wiener. While Wiener apparently did not use neural networks as nonlinear functions (but e.g., Hermite polynomials), neural networks in Wiener models have been studied for some time [2022]. In the present work, Laguerre filtering is presented without making use of the vocabulary of cybernetics. In particular, the 𝑧-transform is not introduced here.

3. Exploiting Similitudes

3.1. Froude-Krylov Forces

This section gives the justification for point 6 of Section 2.4. If the undisturbed fluid in which the cylinder is plunged is accelerating (because of surface waves, e.g.,), then it is natural to introduce two reference systems: 𝔊 is a Galilean reference system, for example, fixed relative to the sea floor and 𝔄 is an accelerated reference system, locally following the undisturbed flow. Transforming the equations of equilibrium from 𝔊 (in which we carry out FEM analysis) to 𝔄 (for which we have experimental data, in water that is not accelerated) requires the addition of inertia forces.

The inertial forces create a uniform pressure gradient that was not present in the laboratory test. The effect of a pressure gradient on a submerged body is variously referred to as “Archimedes forces” when the pressure gradient results from the acceleration of gravity, or as “Froude-Krylov forces” when the pressure gradient is due to fluid acceleration in for example, surface waves. As familiar, the integral of the pressure over the wet surface is transformed into a volume integral [23].

It is assumed that this pressure gradient does not affect the turbulent flow, so that the pressure gradient can simply be added to the pressures resulting from turbulence. This seems reasonable enough for incompressible flows, and indeed when it comes to Archimedes forces, the submerged weight of a cylinder is routinely subtracted to laboratory measurements and the relevant correction added again in FEM analysis—even though the Archimedes forces in the laboratory do not necessarily scale with those in the analysis. To the author’s knowledge, there is no experimental indication that a horizontal and vertical cylinder, all other conditions being equal, experience different forces.

To conclude, the hydrodynamic force acting on the cylinder at a given instant is the sum of two terms: (1)a force that is a function of only the cylinder diameter and the recent history of the velocity of the cylinder relative to the undisturbed, steady water flow, (2)Froude-Krylov forces.

All computations in Sections 4 and 5 deal only with the first of the above two terms.

3.2. Scaling

This section details how points 2 and 5 of Section 2.4 are implemented. In order to reduce the amount of experimental data necessary to create the interpolation function used in point 4, one must take advantage of scale similarities. To that effect, all data used to either train or query the database is scaled. Correspondingly, all forces returned by the database are scaled back. The present model uses a scaling that is quite different from the scaling typical used by simplified harmonic models (Section 1.3): reduced amplitudes and frequencies are not used. Instead velocities histories are scaled in a manner familiar from the Reynolds number.

VIV forces are assumed to be uniquely defined by fluid density 𝜌, kinematic viscosity 𝜈, cylinder diameter 𝐷, and the motion history. Hence, in order to create a database that is to be entered with scaled velocities, we wish all experimental data to be scaled to fixed reference values 𝜌𝑜, 𝜈𝑜, and 𝐷𝑜. The choice of 𝜌𝑜, 𝜈𝑜, and 𝐷𝑜 is arbitrary, and in this work, all are set to the value 1.

By expressing the units of these quantities, one gets three equations on 𝜆𝑚, 𝜆𝑠, and 𝜆𝑘𝑔, which are the scaling factors for the basic units of distance, time, and mass. Solving the system yields𝜆𝑚=1𝐷,𝜆𝑠=𝜈𝐷2,𝜆𝑘𝑔=1𝜌𝐷3.(3.1) Once the scaling of basic units is known, the scaling of any derived quantities, for example, velocities, accelerations, and forces per unit length can be expressed:𝜆𝑚𝑠1=𝐷𝜈,𝜆(3.2)𝑚𝑠2=𝐷3𝜈2,𝜆(3.3)𝑁𝑚1=𝐷𝜌𝜈2.(3.4) The choice 𝐷𝑜=1 [m] hence implies that scaled displacements can be considered to have “1 diameter” as unit. Similarly, the choices 𝐷𝑜=1 [m] and 𝜈𝑜=1 [m2/s] together imply that scaled velocities are expressed as Reynolds numbers since the scaled velocity is calculated as 𝐷𝑣/𝜈 where 𝑣 is the velocity. The Reynolds number is usually computed using some velocity that is characteristic of the system under study. In VIV science, the undisturbed velocity of the current is used. By contrast, in this work, instantaneous local values of the relative velocity vector are multiplied by 𝐷/𝜈. The scaled velocities thus obtained are a generalization of the traditional use of Reynolds number: considering an immobile cylinder in a current, the norm of its scaled relative velocity vector is equal to the traditional Reynolds number. To prevent confusion of the present usage of Reynolds number with the more particular classical one, yet emphasize the relation between both, the expression “ilr-Reynolds” (for “instant, local, relative Reynolds”) will be used in this document.

Since scaling is applied consistently to all derived quantities, all nondimensional numbers based on combinations of distance, time and mass (including Reynolds and Froude numbers) are conserved. However, any dimensional quantity with units different from those of 𝜌, 𝜈, and 𝐷 is scaled to values that depend of 𝜌, 𝜈, and 𝐷. In particular, (3.3) shows that all accelerations, including the acceleration of gravity 𝑔, are scaled with a factor proportional to 𝐷3/𝜈2. So while the scaling used here may conserve Froude's number, it does not allow to build a database of forces related to surface wave effects, because the database does not refer to a uniform value 𝑔𝑜.

4. Characterization of Trajectory

4.1. Foreword

This section details how point 3 in Section 2.4 is to be implemented. The objective is, for any given point in time, to distill a “summary” of the recent history of the velocity of the cylinder relative to the surrounding fluid (trajectory). Note that the history of each component of the velocity vector is treated separately in this section and that the procedure is applied to the scaled trajectory.

The trajectory is approximated as a linear combination of some adequate family of functions, and the coefficients in this linear combination are the summary. The family of functions that is used here is the series of Laguerre polynomials (Section 4.2). It is shown in Section 4.3 that if the “Laguerre coefficients” of the linear combination are obtained by integrating the product of the trajectory by adequate “Laguerre analysis functions”, then the difference between the approximating linear combination and the real trajectory is small in the recent past and larger in the further past. This justifies the choice of Laguerre polynomials: they allow to summarize the trajectory in a way that represents recent velocities very precisely, and older velocities in a coarser manner. It is assumed that this corresponds to the information needed to obtain a good estimate of the hydrodynamic force.

Computing the integral of the product of Laguerre analysis functions and trajectory takes time. Luckily, one can show (Section 4.5) that the Laguerre coefficients are the solution of a differential equation driven by the instant value of the velocity. To obtain results that are independent of step size, this differential equation must be carefully discretised in time (Section 4.6) when summarizing experimental data.

4.2. Definitions

The Laguerre polynomial (Figure 1, top) of degree 𝑖1 can be defined by its Rodrigues formula [24]𝑖𝑒(𝑥)𝑥𝑑𝑖!𝑖𝑑𝑥𝑖𝑥𝑖𝑒𝑥.(4.1) Laguerre polynomials verify the orthonormality property0𝑖(𝑥)𝑗(𝑥)𝑒𝑥𝑑𝑥=𝛿𝑖𝑗.(4.2) We seek to describe the recent trajectory with a precision that is good for the immediate past, and decreasing for the further past. To this end, we introduce a weight function which emphasizes “recent past” (Figure 1, bottom)𝑒𝒲(𝑡)𝑡/𝑡𝑤𝑡𝑤,𝑡,(4.3) where the interpretation of 𝑡𝑤 has been discussed in Section 2.1. Functions will now be noted as vectors (in a Hilbert space), marked with overline symbols. An indexed family of functions will be noted as a matrix (symbols with double overline) and so will a linear operator (a distribution of two variables). We introduce the symmetric positive definite operator𝑊𝑡1,𝑡2𝑡=𝛿1,𝑡2𝒲𝑡1(4.4) and a dot product in a suitable space of real valued functions𝑓𝑇𝑔0𝑓(𝑡)𝑔(𝑡)𝑑𝑡(4.5) with the canonical norm|||𝑓|||𝑓𝑇𝑓.(4.6) Further we introduce the base𝐿(𝑡,𝑖)𝑖𝑡𝑡𝑤,𝑡,𝑖{1,,𝑛}.(4.7) Equation (4.2) can be rewritten in matrix notation as𝐼=𝐿𝑇𝑊𝐿,(4.8) where 𝐼 is the 𝑛×𝑛 identity matrix. It is useful to introduce the weighted norm or 𝑤-norm|||𝑓|||𝑤𝑓𝑇𝑊𝑓.(4.9) Note that since 𝒲(𝑡) is of dimension [1/𝑠],  |𝑓|𝑤 is of the same dimension as 𝑓. So when taking 𝑓 as a scaled velocity, |𝑓|𝑤 is an ilr-Reynolds number.

4.3. Analysis and Synthesis

For a history 𝑣(𝑡) of either the 𝑥 or 𝑦 component of the velocity, we seek the vector of “Laguerre coefficients” 𝜏 with which to combine the columns 𝐿, that minimize the weighted error 𝐽 defined as1𝐽=2|||𝑣𝐿𝜏|||2𝑤=12𝑣𝐿𝜏𝑇𝑊𝑣𝐿𝜏.(4.10) A notation borrowed from physics is used here: the dot in the above equations symbolizes a sum, as would appear in a matrix-vector product or the scalar product of two vectors. In this notation, the sum acts on the last index of the left argument and the first index of the right argument. Vector transpositions are hence without effect, but have been added in the text for readers that prefer matrix notations.

To this effect, we require that the derivative be zero:𝜕𝐽𝜕𝜏=𝐿𝑇𝑊𝐿𝜏𝐿𝑇𝑊𝑣,(4.11) which implies𝜏=𝐿𝑇𝑊𝐿1𝐿𝑇𝑊𝑣=(4.12)𝐿𝑇𝑊𝑣,(4.13)𝜏=𝐷𝑇𝑣(4.14) with𝐷𝑊𝐿.(4.15) The “Laguerre analysis functions” 𝐷 (Figure 1, middle) are by definition equal to𝐷(𝑡,𝑖)=𝒟𝑖𝑡𝑡𝑤=𝑖𝑡𝑡𝑤𝑒𝑡/𝑡𝑤𝑡𝑤.(4.16) The Laguerre analysis functions 𝐷 must not be confused with the Laguerre functions (note the factor 2 in (4.18)). Incidentally, Wiener-Laguerre models use Laguerre filters whose impulse response is Laguerre functions (not analysis function), so the present approach is slightly different from the classical Wiener-Laguerre model. A justification for the present choice will appear in Section 6.1.

4.4. Convergence

Laguerre functions, which can be defined as𝐹(𝑡,𝑖)𝑖𝑡𝑡𝑤(4.17)𝑖𝑡𝑡𝑤𝑒𝑡/(2𝑡𝑤)𝑡𝑤(4.18) or in matrix notation as𝐹=𝑊𝐿,(4.19) have been extensively studied. Series of Laguerre functions are known to converge almost everywhere (under some conditions of continuity) [25]. In matrix notation, this result can be stated aslim𝑛||||𝐹𝐹𝑇𝑓𝑓||||=0.(4.20) This can be used to obtain a result on the convergence of series of Laguerre polynomials. We introduce the change of variables𝑓=𝑊𝑔(4.21) so that||||𝐹𝐹𝑇𝑓𝑓||||=|||||𝐹𝐷𝑇𝑔𝑊𝑔|||||=|||||𝑊𝐿𝐷𝑇𝑔𝑔|||||=||||𝐿𝐷𝑇𝑔𝑔||||𝑤.(4.22) We hence have convergence in terms of the quality of approximation that we are seeking, with emphasis on the recent past. Further, on any finite (or “compact”) interval, convergence in the 𝑤-norm is equivalent to convergence almost everywhere. So under some conditions of continuity on 𝑔, the series of Laguerre polynomials obtained using 𝐷 as analysis functions converges almost everywhere towards 𝑔 in any finite interval.

Figure 2 illustrates how Laguerre coefficients indeed provide a “summary” of the trajectory

4.5. Differential Equation for Laguerre Coefficients

In the finite element analysis, we need to update the Laguerre coefficients at each iteration of each time step, for every Gauss point of every node of the system. The explicit calculation of (4.14) for every update is hence a CPU-time critical operation, taking in the order of 𝑛×𝑁 floating point operations (flops), where 𝑛 is the number of Laguerre polynomial used and 𝑁 the number of time steps that the analysis functions take to decay to a negligible value. Further, for each Gauss point, 2𝑁 velocity values need to be stored, a severe memory requirement.

In the present section and the next, it is shown how the computation of (4.14) can be carried out by a recursive operation requiring no other storage than that of the Laguerre coefficients and the last velocity values, and taking in the order of 𝑛×𝑛 flops, which is advantageous because 𝑛𝑁. In this section, it is shown that 𝜏 verifies a differential equation driven by the history 𝑣 of the velocity component. In Section 4.6, this differential equation is solved time-step by time-step in a recursive update.

Equation (4.14) can be rewritten without matrix notation and differentiated𝜕𝜏𝑖=𝜕𝑡0+𝑒𝜃𝑖(𝜃)𝜕𝑣𝜕𝑡𝑡𝑡𝑤𝜃1𝑑𝜃=𝑡𝑤0+𝑒𝜃𝑖(𝜃)𝜕𝑣𝜕𝜃𝑡𝑡𝑤𝜃𝑑𝜃.(4.23) Multiplying by 𝑡𝑤 and integrating by parts yields𝑡𝑤𝜕𝜏𝑖𝑒𝜕𝑡=𝜃𝑖(𝜃)𝑣𝑡𝑡𝑤𝜃+0+𝑒𝜃𝑖(𝜃)+𝑒𝜃𝜕𝜕𝜃𝑖(𝑣𝜃)𝑡𝑡𝑤𝜃𝑑𝜃.(4.24) A property of Laguerre polynomials is [24]𝜕𝜕𝜃𝑖(𝜃)=(1)𝑖1(𝜃)=𝑖1𝑗=0𝑗(𝜃),(4.25) where 𝑖(1)(𝜃) is a generalized Laguerre polynomial. Hence we can write𝑡𝑤𝜕𝜏𝑖𝜕𝑡=𝑖(0)𝑣(𝑡)𝜏𝑖0+𝑒𝜃𝑖1𝑗=1𝑗(𝜃)𝑣𝑡𝑡𝑤𝜃𝑑𝜃=𝑣(𝑡)𝜏𝑖𝑖1𝑗=1𝜏𝑗=𝑣(𝑡)𝑖𝑗=1𝜏𝑗.(4.26) which is of the form𝜕𝜏𝜕𝑡(𝑡)=𝜇𝜏(𝑡)+𝑛𝑣(𝑡)(4.27) with𝜇𝑖𝑗=1𝑡𝑤𝑛𝑗𝑖0𝑗>𝑖,𝑖=1𝑡𝑤.(4.28) Equation (4.27) shows that at any time 𝑡, the rate of the Laguerre coefficients is fully defined by the Laguerre coefficients and the velocity signal.

4.6. Recursive Filter

The discrete integration of (4.27) must be done carefully, for two reasons. Firstly, it is important to obtain Laguerre coefficients that are independent of the sampling rate used (as long as the sampling rate is “adequate”). This is because the experimental data on which the VIV model is based may come from experiments which, after scaling, may have different sampling rates. Further, the numerical analysis in which the VIV model is used may use yet another time step. The choice of time step or sampling rate must not affect the way a trajectory is characterized by Laguerre coefficient.

The second reason for care in discrete integration is that we wish to be able to create synthesized signals 𝐿𝜏 of good quality. Synthesized signals are neither used in the numerical process of creating a force interpolation function (Section 5) or in the FEM use of the VIV model. However, visualization is essential to the process of research, both for fault diagnosis and quality control, and to communicate an understanding of the method.

This discrete integration is only used in the analysis of experimental data, to provide an input to the training of the “rotatron” (Section 5.6). In dynamic analysis, the integration of (4.27) is done by means of the Newmark-𝛽 method, as detailed in Section 7.

Assume that velocity is sampled at regular intervals𝑣𝑗𝑡=𝑣0.+𝑗𝑑𝑡(4.29) We seek the values of the Laguerre coefficients at the same intervals𝜏𝑗=𝜏𝑡0+𝑗𝑑𝑡.(4.30) The vector 𝜏𝑗 (the list of the coefficients for all Laguerre polynomials, taken at step 𝑗) must not be confused with scalar 𝜏𝑖 (the coefficient for the Laguerre polynomial of degree 𝑖). We choose 𝑡0 such that 𝑡0+𝑗𝑑𝑡=0, and we approximate 𝑣 by a function that is linear over the interval [0,𝑑𝑡]. Equation (4.27) becomes𝜕𝜏𝜕𝑡(𝑡)=𝜇𝜏(𝑡)+𝛼+𝛽𝑡(4.31) with𝛼=𝑛𝑣(0),𝛽=𝑛𝑣(𝑑𝑡)𝑣(0).𝑑𝑡(4.32) This new differential equation can be solved exactly: we seek a solution of the form𝜏(𝑡)=exp𝜇𝑡𝑎+𝑏𝑡+𝑐(4.33) over the interval. Here exp(𝜇𝑡) stands for a matrix exponential. Replacing this expression into (4.31), noting that𝜕𝜕𝑡exp=𝜇𝑡𝜇exp,𝜇𝑡exp0=𝐼,(4.34) and identifying the constant and linear terms and enforcing the initial value leads to𝑏=𝜇1𝛽,𝑐=𝜇2𝛽𝜇1𝛼,𝑎=𝜏(0)+𝜇2𝛽+𝜇1𝛼.(4.35) Replacing these expressions in (4.33) at 𝑡=𝑑𝑡, a tedious but straightforward computation yields the recursive filter𝜏𝑗+1=𝑀𝜏𝑗+𝑉1𝑣𝑗+𝑉2𝑣𝑗+1(4.36) with𝑀=exp,𝑚𝑑𝑡𝜇1=𝜇1𝑛,𝜇2=𝜇2𝑛1,𝑑𝑡𝑉1=𝑀𝜇1𝜇2+𝜇2,𝑉2=𝑀𝜇2𝜇1𝜇2.(4.37)

5. Force Interpolation

5.1. Foreword

This section details the implementation of point 4 in Section 2.4. This section presents an interpolation function which, given the Laguerre coefficients, predicts the present value of the force vector. Polynomials were considered initially, but it soon became clear that feed-forward “neural networks” provide a better class of functions to work with. The reason for this is that the number of polynomial coefficients of degree 𝑑 for a polynomial of 𝑛 variables is 𝑛𝑑, and high values of 𝑑 must be expected to be necessary. By contrast, in a neural network, nonlinearity is introduced by “sigmoid” or “threshold” functions, and the coefficients are used to specify in which direction nonlinearity applies. Further, polynomials are infamous for their propensity to oscillate.

The “rotatron” presented here is based on the “perceptron” [26, 27], a well-studied architecture of neural network which provides a flexible tool for the interpolation of scalar-valued functions of a vector (Section 5.2). The rotatron takes advantage of certain symmetry properties of the physics at hand (Section 5.3).

In Section 7, the rotatron is used to predict scaled forces based on the Laguerre coefficients for scaled trajectories.

5.2. Perceptron

The perceptron [26, 27] is a simple feed-forward neural network, consisting of 3 layers. The input layer has 2𝑛 neurons where 𝑛 is the number of Laguerre coefficients for each velocity component and the factor 2 comes from the need to analyze in-line and cross-flow speed histories together. The values of the input layer neurons are set to the Laguerre coefficients for both velocity components. The second layer has nhid neurons, whose values are an affine function of the values of the first layer, passed through a sigmoid function like2𝜎(𝑥)=1𝑒2𝑥.+1(5.1) Finally, the third layer gives the output of the perceptron, and its values are an affine function of the values of the second layer. This can be summarized as𝑓𝑖=𝑀𝑖𝑗𝑁𝜎𝑗𝑘𝑙𝜏𝑘𝑙+𝑉𝑗+𝑈𝑖.(5.2)𝑀𝑖𝑗, 𝑁𝑗𝑘𝑙, 𝑈𝑖 and 𝑉𝑗 are the “weights” or interpolation coefficients, that must be adjusted to fit the perceptron to interpolate some given data. 𝜏𝑘𝑙 are Laguerre coefficients and 𝑓𝑖 are predicted force components. 𝑖 is the index of force direction (𝑥 versus 𝑦), 𝑗 the index of neuron in the hidden layer, 𝑘 the index of velocity direction, and 𝑙 the index of Laguerre coefficient.

Each output of the perceptron can be seen as a function, which is a sum of sigmoid steps in directions defined by 𝑁𝑗𝑘𝑙.

5.3. Symmetries

The relation between trajectories (in the sense of the history of the velocity of the cylinder relative to the water) and forces can reasonably be assumed to exhibit several symmetries (Figure 3). Rotational symmetry: if a trajectory can be deduced from the other by a rotation around the origin, then the resulting forces are also deduced from each other by the same rotation. Mirror symmetry: if a trajectory can be deduced from the other by a mirroring around a line crossing the origin, then the resulting forces are also deduced from each other by the same mirroring.

Rotational symmetry and mirror symmetry together imply directionality: if a trajectory is within a line crossing the origin, then the resulting forces are within the same line. In particular, zero velocities must imply zero forces.

The symmetries imply that once experimental data for a trajectory has been obtained, there is no need to acquire data for rotated or mirrored trajectories. However, if one was training a perceptron to interpolate the data, the training set would need to include trajectories and their rotates and mirrors, with the correspondingly rotated and mirrored forces. This would increase memory and CPU usage during training, but also during use of the trained perceptron, because the perceptron will need a larger number of hidden layers to interpolate the training data.

Another approach is hence used in the present work: the classic perceptron is replaced by a “rotatron” (Section 5.5). It is designed so that, whatever the values of the weight coefficient, a rotation or mirroring of the input trajectory results in the same rotation or mirroring of the output force vector.

5.4. Index Notations

In the present work, index notations inspired from tensor analysis are used. However, the present setting differs from tensor analysis in at least three ways.

Firstly, we assume that we are only operating in Euclidean spaces (and not in more general Riemannian manifolds) so that orthogonal bases can be used. This makes it unnecessary to distinguish between co- and contravariant bases and coordinates. Hence, only lowered indexes appear in the present work. Incidentally, it was here assumed that the state of the model is a point in a vector space, which is not true when finite rotations are present and Riemannian geometry should be introduced instead.

Secondly, in tensor notations, each index spans the dimension of the manifold. In an expression like 𝜎𝑖𝑗=𝐶𝑖𝑗𝑘𝑙𝜀𝑘𝑙, the indexes range from 1 to 3. Following Einstein's convention, indexes 𝑘 and 𝑙 are summed over, and the relation is valid for any combination of 𝑖 and 𝑗. The fact that the equation is valid at each point within a solid is implicit in the notation. In the present work, we prepare for the manipulations of arrays in a computer, involving operations that are repeated, for example, for various locations along a riser. If indexes 𝑥, 𝑦 and 𝑧 were introduced to note the position to which the various tensors refer, one would tend to write 𝜎𝑖𝑗𝑥𝑦𝑧=𝐶𝑖𝑗𝑘𝑙𝑥𝑦𝑧𝜀𝑘𝑙𝑥𝑦𝑧, which violates Einstein's convention, because no summation (or rather: no integral) is implied over the positions.

Thirdly, we introduce nonlinear functions. These functions can combine the values of the coordinates for some indexes (these will be noted in brackets) and operate in parallel on the coordinates for other indexes (not brackets). For example, by definition of the notation |𝑦[𝑗]𝑘| is equal to 𝑦21𝑘+𝑦22𝑘, all values of 𝑗 are present in one evaluation of the square root, and the square root is evaluated for each value of 𝑘.

Appendix A gives a detailed description of the conventions used.

5.5. Rotatron

A modified interpolation function (which will be referred to as “rotatron” in this text) is defined as𝑓𝑖=𝑉𝑘𝜎𝑖𝑘(5.3) with𝜎𝑖𝑘=𝜎𝑖𝑦[𝑗]𝑘=𝑦(5.4)𝑖𝑘||𝑦[𝑗]𝑘||||𝑦1+[𝑗]𝑘||𝛼𝑘,||𝑦(5.5)[𝑗]𝑘||=𝑦21𝑘+𝑦22𝑘,𝛼(5.6)𝑘=1𝑒𝑈𝑘,𝑦(5.7)𝑗𝑘=𝑀𝑘𝑙𝜏𝑗𝑙.(5.8) In the above, indexes 𝑖 and 𝑗 refer to direction, index 𝑙 to the Laguerre polynomial, and 𝑘 to the hidden layer. 𝑉𝑘, 𝑈𝑘, and 𝑀𝑘𝑙 are tunable parameters. 𝜏𝑗𝑙 are Laguerre coefficients, given as input to the “rotatron”. Note that (4.14), (5.3), and (5.8) operate linearly, identically, and independently of the terms related to the 𝑥 and 𝑦 directions, while (5.5) involves a unit vector multiplied by a nonlinear function of its norm.

The rotatron can be shown to enforce the symmetries discussed in Section 5.3: if 𝜏𝑗𝑙 is multiplied by a 2 by 2 matrix 𝑠𝑗𝑗 representing a rotation or a mirroring, then 𝑦𝑗𝑘 is multiplied by the same matrix. The Euclidean norm |𝑦[𝑗]𝑘| is hence unchanged. The term 𝑦𝑖𝑘 in (5.5) is multiplied by 𝑠𝑖𝑗, hence so are 𝜎𝑖𝑘 and 𝑓𝑖: a rotation or mirroring of the Laguerre coefficient results in the same transformation of the forces.

Figure 4 illustrates the flow of information, from right to left, from two vectors containing the histories of the velocity components, to Laguerre coefficient, that are then processed in the rotatron.

The nonlinear function appearing in (5.5) is a sigmoid, whose abruptness is parametrized by 𝑈𝑘 (Figure 5). The sigmoid is shown in Figure 5 for various values of the parameter 𝑈𝑘.

5.6. Training

“Training” of a neural network refers to finding weight coefficients 𝑉𝑘, 𝑀𝑘𝑙, and 𝑈𝑘 such that for any training point number 𝑚, consisting of Laguerre coefficients 𝜏𝑗𝑙𝑚 and two force components 𝑓𝑖𝑚, the outputs 𝑓𝑖𝑚 computed by the neural network are close to 𝑓𝑖𝑚.

5.6.1. Regularization

A common problem when training neural networks is “overspecialization” [28]. In this situation, the neural network predicts the training outputs with high accuracy but behaves wildly between the training points. In contrast, what is implicitly sought is a smooth response of the network to the input, even if this means an imperfect fit to the training data.

Many strategies are described in the literature to address this problem. One of them, which is adopted here, is regularization [28]: the value of the weight parameters 𝑉𝑘, 𝑀𝑘𝑙, and 𝑈𝑘 are chosen by minimizing the cost function𝐽𝑉[𝑘],𝑈[𝑘],𝑀[𝑘],𝑓[𝑖𝑚],𝜏[𝑗𝑙𝑚]=12𝑓𝑖𝑚𝑓𝑖𝜏[𝑗𝑙]𝑚21+𝜌2𝑈2𝑘+𝑉2𝑘+𝑀2𝑘𝑙.(5.9) The “regularization coefficient” 𝜌 is an arbitrary input to the training algorithm. High values of 𝜌 favor smoothness of the response of the neural network against precision in reproducing the training set.

Arguments of symmetry by permutation of the numbering of the first index of 𝑦𝑗𝑘 in (5.6) show that the cost function has multiple minima. Further, there are probably local minima higher than the lowest maxima.

5.6.2. Conjugate Gradient Optimization

𝐽 is a function of a large number of weight coefficients, and hence it is not practical to compute the Hessian of 𝐽, because the Hessian is a full matrix. It also proves to be very costly to even compute an approximation to it as done in the Levenberg-Marquardt algorithm [29, 30]. On the other hand, the Nelder-Mead “downhill simplex” algorithm [31], which uses only the values of 𝐽, proved very slow in this case. Hence a search method is chosen, that determines the search direction from the gradient of 𝐽 [32]. This is a conjugate gradient method, in which the step length is found by deriving the gradient in the direction of the search. In this method, the positive definiteness of the (implicit) Hessian is forced by adding a scaled identity matrix to it, a technique known as “trust region”.

The conjugate gradient method proved far more efficient than the Levenberg-Marquardt and Nelder-Mead methods for the present task.

The weight coefficients are set to random values at the start of the conjugate gradient iterations.

6. Metric

6.1. Euclidean Metric and Distance

In order to describe the available data, it is useful to define a distance between trajectories. This will allow to determine to what extend the set of available data “fills” the set of all possible trajectories, or to detect zones of transition from one hydrodynamic behavior to the other. Finally, this will help detecting contradictions in the available data, arising from a variety of sources, including hidden experimental variables, measurement uncertainties, or inadequate modeling in inverse methods and not least, the natural variability of VIV forces.

The 𝑥 and 𝑦 components of a trajectory are described by a pair of functions:𝑓𝑓𝑥,𝑓𝑦.(6.1) We can define a scalar product between trajectories, that captures any recent differences:𝑓𝑇𝑔𝑓𝑇𝑥𝑊𝑔𝑥+𝑓𝑇𝑦𝑊𝑔𝑦=0𝑒𝑡/𝑡𝑤𝑡𝑤𝑓𝑥(𝑡)𝑔𝑥(𝑡)+𝑓𝑦(𝑡)𝑔𝑦(𝑡)𝑑𝑡.(6.2) By replacing 𝑓𝑥, 𝑓𝑦, 𝑔𝑥, and 𝑔𝑦 by their expression in terms of Laguerre polynomials and their respective Laguerre coefficients 𝜏𝑓𝑥, 𝜏𝑓𝑦, 𝜏𝑔𝑥, and 𝜏𝑔𝑦, one finds that𝑓𝑇𝑔=𝜏𝑇𝑓𝑥𝜏𝑔𝑥+𝜏𝑇𝑓𝑦𝜏𝑔𝑦=𝜏𝑇𝑓𝜏𝑔(6.3) with𝜏𝑓𝜏𝑓𝑥𝜏𝑓𝑦,𝜏𝑔𝜏𝑔𝑥𝜏𝑔𝑦.(6.4) The distance is defined from the scalar product in the usual manner:|||𝑓𝑔|||𝑓𝑔𝑇𝑓𝑔=||,(6.5)𝜏𝑓𝜏𝑔||.(6.6) In other words, neighboring vectors of Laguerre coefficients describe trajectories that are similar in the recent past. This is illustrated by taking random samples of Laguerre coefficients around a given value obtained from data analysis and plotting the synthesized trajectories (Figure 6). Given that trajectories are compared using an integral weighted with an exponential (4.2). This important property can be seen as the justification of the choice of the Laguerre coefficients to “summarize” trajectories.

6.2. Rotatron Distance

The above does not account for rotational and mirror symmetries. We seek a distance for which the distance of a trajectory to its transforms by rotation or mirroring is zero. Another distance is hence introduced:𝑑𝑓,𝑔minmin𝑅|||𝑓𝑅𝑔|||,min𝑆𝔖|||𝑓𝑆𝑔|||,(6.7) where is the set of all rotations of the trajectories around the origin and 𝔖 the set of all mirroring of trajectories around a line passing by the origin. Note that no norm or scalar product associated to the distance 𝑑 is presented here (The vector-space of trajectories, divided by the group of rotations and mirrorings, is not a vector space).

Because 𝑓𝑥 is related to 𝜏𝑓𝑥 by the same linear relation that relates 𝑓𝑦 to 𝜏𝑓𝑦, linear combinations of 𝑓𝑥 and 𝑓𝑦 (including rotation and mirroring) are related to the same linear combinations of 𝜏𝑓𝑥 and 𝜏𝑓𝑦. By expressing the distance |𝑓𝑅(𝑔)| as a function of the angle 𝛼 of the rotation 𝑅, and then differentiating it with respect to 𝛼, it can be shown that the value of 𝛼 that minimizes |𝑓𝑅(𝑔)| is𝛼=arctan𝜏𝑓𝑥𝜏𝑔𝑦𝜏𝑓𝑦𝜏𝑔𝑥,𝜏𝑓𝑦𝜏𝑔𝑦+𝜏𝑓𝑥𝜏𝑔𝑥,(6.8) where arctan(𝑦,𝑥)]𝜋,𝜋] is the angle of a vector [𝑥,𝑦]𝑇 with the 𝑥-axis. Similarly, it can be shown that the mirroring that minimizes |𝑓𝑆(𝑔)| is the composition of a rotation of angle𝛽=arctan𝜏𝑓𝑥𝜏𝑔𝑦+𝜏𝑓𝑦𝜏𝑔𝑥,𝜏𝑓𝑦𝜏𝑔𝑦𝜏𝑓𝑥𝜏𝑔𝑥(6.9) by a swap of the sign of the 𝑥-coordinates. Equations (6.8) and (6.9) allow to compute (6.7).

Figure 7 shows a trajectory and the trajectories within a small database that have the smallest distance to it, measured using 𝑑.

6.3. Existence of Functional Relation
6.3.1. Introduction

From experiments, one can obtain databases of velocity histories, and the corresponding hydrodynamic forces. The velocity histories can be “summarized” into Laguerre coefficients. An open question at this stage is whether there is actually a functional relationship between the Laguerre coefficients and the forces, that the rotatron could be used to approximate. For example, if the decay time 𝑡𝑤 was chosen too small, then the same velocity in the recent past (represented by the Laguerre coefficient) will result in different forces, depending on the velocity in the more remote past. In that case, there is no functional relation. The functional relation could also be lost if too few Laguerre coefficients are used, or if the force in the experiments is affected by incoming turbulence. In the absence of a functional relation, attempting to train the rotatron will not give useful results, in the same way as no curve drawn through a cloud of points on a piece of paper will provide a useful model.

The approach in this section hinges on two ideas. Firstly, the higher the dimension of a spaceis, the faster the volume of a ball increases with its radius. Hence, the statistic of the distances between points in a cloud of data can be used to define a dimension (Section 6.3.2). For example, a set of points on a piece of paper will be shown to follow a curve if the number of neighbors increases linearly with the radius. Secondly, if Laguerre coefficients 𝜏 predict forces 𝑓, then the dimension of the cloud of experimental 𝜏 values must be the same as the dimension of the cloud of [𝜏,𝑓] points (Section 6.3.3).

6.3.2. Minkowski-Bouligand Dimension

Considering a relatively uniform cloud of points, the number 𝑚 of points in a ball is proportional to the radius 𝑟 of the sphere to the power of 𝑝, where 𝑝 is the dimension of the space in which the cloud is defined. For example, using 2×10 Laguerre coefficients to describe both components of a trajectory, the number of points in the ball would be 𝑚𝑟20.

Conversely, one can define the Minkowski-Bouligand dimension (often referred to as the fractal dimension [33]) 𝑝 of a set (in particular, of a “database” of Laguerre coefficients) by counting the number 𝑚(𝑟) of pairs of points in a set, which have a distance smaller than 𝑟:𝑝(𝑟)𝜕log𝑚(𝑟).𝜕log𝑟(6.10) In practice, one should either smooth 𝑚(𝑟) or compute the derivative by finite differences over a large enough interval. Note that the fractal dimension 𝑝 is a function of the scale 𝑟.

6.3.3. Functional Relation

Imagine that we have a series of data-points (𝑥,𝑦,𝑧), and we are investigating whether 𝑧 can be predicted using 𝑥 and 𝑦. Let us imagine that the fractal dimension of the set of (𝑥,𝑦) pairs is 2 (the set of (𝑥,𝑦) fills the plane). If the fractal dimension of the set of (𝑥,𝑦,𝑧) is equal to 2, then the set of (𝑥,𝑦,𝑧) is within a surface, and 𝑧 can be predicted using 𝑥 and 𝑦. If the fractal dimension of the set of (𝑥,𝑦,𝑧) is equal to 3, the data forms a cloud, and 𝑥 and 𝑦 are not sufficient to predict 𝑧, therefore other hidden variables must be at play. These concepts are now applied to the study of the experimental database.

Figure 8 shows the cumulative distribution of the distances between trajectories (“𝜏-distance”, black curve) computed using (6.7). 𝑝 is seen to depend on the scale: from afar (𝑟>2×104[ilr Re]), the slope of the curve is zero, hence the dimension is zero: all the data are lumped into a point. Zooming into the data set (𝑟=3×103 [ilr Re]), one can discern a cloud of dimension 4.76. At 𝑟=1.5×103 [ilr-Re], the slope decreases to about 𝑝=2, and it is believed that this is the dimension of the data set for a given point along the riser. At a small scale (𝑟<1×103 [ilr Re]), the dimension increases again, possibly due to noise in the data. or weaknesses in the Laguerre approximation.

The red curves in Figure 8 are computed by adding the sum of squares of the differences between force components (suitably scaled) to the squares of the distances between trajectories, and then extracting the square root (“𝜏𝑓-distance”). The four red curves, from left to right, are drawn using the original force data, to which Gaussian noise of standard deviation 0, 107, 108, and 109 [N/m], respectively, has been added (The standard deviation of the original force is about 2×108 [N/m]). The first curve is close to the black one, which is strongly suggestive that indeed, there is a functional relation. The two first red curves are indistinguishable, which seems to indicate that we cannot expect to achieve a 10% precision in force predictions. The marked difference with curves 3 and 4 shows, however, that we have assets in hand to predict the force. Similar curves have been produced with added noise of standard deviations 1×107, 2×107 and 10×107, [N/m], and already at 2×107 [N/m] the curve is distinct from the one based on the original data.

6.3.4. Discussion

The present study suggests that there is indeed a functional relation to be seen in the data set used here. The Laguerre coefficients can be used to predict the forces, with a precision of about 107 [N/m]. The conclusion must be treated carefully, however. The dimensional analysis provides a necessary condition, not a sufficient one: it does not exclude, for example, that there exists a neighborhood of points 𝜏 in which two distinct values of 𝑓 appear.

7. Dynamic Analysis

7.1. Foreword

Once it is possible to predict hydrodynamic forces on a cross section for a given velocity history, the next development is to include the force thus predicted in a dynamic time domain simulation. Because the VIV forces introduce severe nonlinearities, a naive connection (where the forces are just added to the right-hand side of the system) might lead to slow convergence or to divergence of the Newton-Raphson iterations used at each time step. To obtain a proper formulation, it is necessary to jointly treat the system of differential equations composed of the state equations of the structure and the differential equations (4.27) followed by the Laguerre coefficient. However, in doing so, for each displacement degree of freedom, 𝑛 Laguerre coefficients are added, and it is crucial for efficiency to eliminate them before solving a large linear system of equations.

To this effect, in this section, the following sequence of transformations is applied to the differential equations. (1)The differential equations are first set in incremental form (Section 7.3). (2)Time discretisation by the Newmark-𝛽 method is introduced (Section 7.4). (3)The Laguerre coefficients are condensed out of the system of equations (Section 7.5). (4)Finite element interpolation is introduced (space discretisation), using Gauss quadrature (Section 7.6).

This particular sequence leads to a VIV model that is implemented at the Gauss point level and can easily be introduced in a general purpose FEM software with standard, displacement-based beam, or cable elements. Another sequence, 1, 4, 2, 3, can be used to obtain either a hybrid element, or alternatively, a mixed element which would require a specialized solver for optimal efficiency. These alternatives are more difficult to integrate into existing software working with displacement-based elements and are not discussed here.

7.2. Differential Equations

The dynamic differential equation of a 3D beam subjected to VIV loads can be formalized as𝑟𝑑𝑖𝑥[𝑏𝑗],̇𝑥[𝑏𝑗],̈𝑥[𝑏𝑗],𝑡=𝜆1𝑁𝑚1𝑓𝑑𝜏[𝑝𝑏]𝑖+𝐸𝑑𝑖,(7.1) where Newton's “dot” notation for a time derivative stands for a derivation with respect to unscaled time 𝑡, as opposed to scaled time 𝑡, and with [23, 34]𝐸𝑑𝑖=𝐶𝐿̇𝑤𝜌𝜈𝑑𝑖̇𝑥𝑑𝑖+𝐶𝑄12𝜌𝐷𝑖||̇𝑤𝑑𝑖̇𝑥𝑑𝑖||̇𝑤𝑑𝑖̇𝑥𝑑𝑖+𝐶𝑀𝜋4𝜌𝐷2𝑖̈𝑤𝑑𝑖̈𝑥𝑑𝑖+𝜋4𝜌𝐷2𝑖̈𝑤𝑑𝑖.(7.2) The four terms in the above Morison's equation are the linear drag, the quadratic drag, the sum of diffraction and added mass forces, and the Froude-Krylov forces. The fourth term introduces the correction discussed in Section 3.1.

If 𝐶𝐿, 𝐶𝑄, or 𝐶𝑀 are set to values different from zero, then it is necessary to subtract the corresponding values from the forces 𝑓𝑑𝑖 used to train the rotatron. Experience shows that using 𝐶𝑀=1, 𝐶𝑄=1, and 𝐶𝐿=0 contributes to the stability of the dynamic analysis.

Equation (4.27) must be scaled to keep only derivatives with respect to unscaled time, for the application of Newmark-𝛽 (Section 7.4)𝜕𝜏𝑙𝑏𝑖𝜕𝑡=𝜇𝑙𝑝𝜏𝑝𝑏𝑖+𝑛𝑙𝜆1𝑚𝑠̇𝑥𝑏𝑖̇𝑤𝑏𝑖(7.3) so that𝜆𝑠1̇𝜏𝑙𝑏𝑖=𝜇𝑙𝑝𝜏𝑝𝑏𝑖+𝑛𝑙𝜆𝑚𝑠1̇𝑥𝑏𝑖̇𝑤𝑏𝑖.(7.4) The indexes 𝑑 and 𝑏 span pairs of directions, orthogonal to the cylinder. Indexes 𝑖 and 𝑗 stand for positions along the cylinder and span a continuous set of values (coordinates along the cylinder). Indexes 𝑙 and 𝑝 refer to the Laguerre coefficients of various degrees. Forces 𝑓𝑑𝑖=𝑓𝑑(𝜏[𝑝𝑏]𝑖) at location 𝑖 only depend on the Laguerre coefficients 𝜏𝑝𝑏𝑖 for the same location. At that location, the force component in direction 𝑑 depends on the Laguerre coefficients of all degrees 𝑏 for both directions 𝑝. 𝜌 is the fluid density and ̈𝑤𝑑𝑖 is the acceleration of the undisturbed fluid. 𝜋/4𝜌𝐷2𝑖̈𝑤𝑑𝑖 stands for the Froude-Krylov forces. Diffraction forces are present in the laboratory tests and hence accounted for by 𝑓𝑑.

7.3. Incremental Form

The incremental form of (7.1) and (7.4) is𝑟𝑑𝑖+𝑘𝑑𝑖𝑏𝑗𝑑𝑥𝑏𝑗+𝑐𝑑𝑖𝑏𝑗𝑑̇𝑥𝑏𝑗+𝑚𝑑𝑖𝑏𝑗𝑑̈𝑥𝑏𝑗=𝜆1𝑁𝑚1𝑓𝑑𝑖+𝑑𝑖𝑝𝑏𝑗𝑑𝜏𝑝𝑏𝑗+𝐸𝑑𝑖,𝜆(7.5)𝑠1̇𝜏𝑙𝑏𝑖+𝑑̇𝜏𝑙𝑏𝑖=𝜇𝑙𝑝𝜏𝑝𝑏𝑖+𝑑𝜏𝑝𝑏𝑖+𝑛𝑙𝜆𝑚𝑠1̇𝑥𝑏𝑖+𝑑̇𝑥𝑏𝑖̇𝑤𝑑𝑖(7.6) with𝑘𝑑𝑖𝑏𝑗=𝜕𝑟𝑑𝑖𝜕𝑥𝑏𝑗,𝑐𝑑𝑖𝑏𝑗=𝜕𝑟𝑑𝑖𝜕̇𝑥𝑏𝑗+𝐶𝑄𝜌𝐷𝑖𝛿𝑖𝑗||̇𝑤[𝑝]𝑖̇𝑥[𝑝]𝑖||𝛿𝑏𝑑+̇𝑤𝑑𝑖̇𝑥𝑑𝑖̇𝑤𝑏𝑖̇𝑥𝑏𝑖||̇𝑤[𝑝]𝑖̇𝑥[𝑝]𝑖||1+𝐶𝐿𝜌𝜈𝛿𝑖𝑗𝛿𝑏𝑑,𝑚𝑑𝑖𝑏𝑗=𝜕𝑟𝑑𝑖𝜕̈𝑥𝑏𝑗+𝐶𝑀𝜋4𝜌𝐷𝑖𝛿𝑖𝑗𝛿𝑏𝑑,𝑑𝑖𝑝𝑏𝑗=𝜆1𝑁𝑚1𝜕𝑓𝑑𝑖𝜕𝜏𝑝𝑏𝛿𝑖𝑗.(7.7) The expression for 𝜕𝑓𝑑𝑖/𝜕𝜏𝑝𝑏 is presented in Appendix B.

7.4. Time Discretisation

Newmark-β is a method geared towards second order differential equations. Equation (7.6), however, is only of the first order. The reason for using Newmark-β here is that the present VIV model is to be integrated into the model of a larger structure, and the differential system for that structure is of the second order. Preparing a first order equation for a second order solver opens two options: we can treat (7.6) as being of the second order in 𝜏𝑙𝑏𝑖, but with the coefficient of ̈𝜏𝑙𝑏𝑖 being zero. Alternatively, we can introduce the antiderivative 𝑇𝑙𝑏𝑖 of 𝜏𝑙𝑏𝑖, and treat (7.6) as being of the second order in 𝑇𝑙𝑏𝑖, but with the coefficient of 𝑇𝑙𝑏𝑖 being zero. The latter option was chosen, based on the weak justification that this treats 𝜏𝑙𝑏𝑖 and ̇𝑥𝑏𝑗 both as first derivatives, which seems natural considering (4.14).

Applying Newmark-𝛽 to (7.5) and (7.6) in this way yields𝑘𝑑,𝑖,𝑑𝑖𝑏𝑗+𝛾𝑐𝛽𝑑𝑡𝑑𝑖𝑏𝑗+1𝛽𝑑𝑡2𝑚𝑑𝑖𝑏𝑗𝑑𝑥𝑏𝑗𝛾𝛽𝑑𝑡𝑑𝑖𝑝𝑏𝑗𝑑𝑇𝑝𝑏𝑗=𝜆1𝑁𝑚1𝑓𝑑𝑖+𝐸𝑑𝑖𝑟𝑑𝑖+𝑐𝑑𝑖𝑏𝑗𝑏𝑥𝑏𝑗+𝑚𝑑𝑖𝑏𝑗𝑎𝑥𝑏𝑗𝑑𝑖𝑝𝑏𝑗𝑏𝜏𝑝𝑏𝑗,𝛾𝑛𝛽𝑑𝑡𝑙𝜆𝑚𝑠1𝑑𝑥𝑏𝑖+1𝛽𝑑𝑡2𝜆𝑠1𝛿𝑙𝑝𝛾𝜇𝛽𝑑𝑡𝑙𝑝𝑑𝑇𝑝𝑏𝑖=𝑛𝑙𝜆𝑚𝑠1̇𝑥𝑏𝑖̇𝑤𝑏𝑖+𝜇𝑙𝑝𝜏𝑝𝑏𝑖𝜆𝑠1̇𝜏𝑙𝑏𝑖𝑛𝑙𝜆𝑚𝑠1𝑏𝑥𝑏𝑖𝜇𝑙𝑝𝑏𝜏𝑝𝑏𝑖+𝑎𝜏𝑙𝑏𝑖(7.8) with𝑎𝑥𝑏𝑗=1𝛽𝑑𝑡̇𝑥𝑏𝑗+12𝛽̈𝑥𝑏𝑗,𝑏𝑥𝑏𝑗=𝛾𝛽̇𝑥𝑏𝑗+𝛾2𝛽1𝑑𝑡̈𝑥𝑏𝑗,𝑎𝜏𝑝𝑏𝑗=1𝜏𝛽𝑑𝑡𝑝𝑏𝑗+12𝛽̇𝜏𝑝𝑏𝑗,𝑏𝜏𝑝𝑏𝑗=𝛾𝛽𝜏𝑝𝑏𝑗+𝛾2𝛽1𝑑𝑡̇𝜏𝑝𝑏𝑗.(7.9) For refinement iterations, 𝑎𝑥𝑏𝑗, 𝑏𝑥𝑏𝑗, 𝑎𝜏𝑝𝑏𝑗, and 𝑏𝜏𝑝𝑏𝑗 are set to zero. Typically, 𝛾=1/2, 𝛽=1/4. The step 𝑑𝑡 refers to unscaled time.

As usual in the Newmark-β method, the increments for the time derivatives are found from the increment as𝑑̇𝑥𝑏𝑗=𝛾𝛽𝑑𝑡𝑑𝑥𝑏𝑗𝑏𝑥𝑏𝑗,(7.10)𝑑̈𝑥𝑏𝑗=1𝛽𝑑𝑡2𝑑𝑥𝑏𝑗𝑎𝑥𝑏𝑗,(7.11)𝑑𝜏𝑝𝑏𝑗=𝛾𝛽𝑑𝑡𝑑𝑇𝑝𝑏𝑗𝑏𝜏𝑝𝑏𝑗,(7.12)𝑑̇𝜏𝑝𝑏𝑗=1𝛽𝑑𝑡2𝑑𝑇𝑝𝑏𝑗𝑎𝜏𝑝𝑏𝑗.(7.13)

7.5. Condensation

The time discrete equations can be rewritten in a compact form:𝑠1𝑑𝑖𝑏𝑗𝑑𝑥𝑏𝑗𝑠2𝑑𝑖𝑝𝑏𝑗𝑑𝑇𝑝𝑏𝑗=𝑠3𝑑𝑖,𝑠4𝑙𝑑𝑥𝑏𝑖+𝑠5𝑙𝑝𝑑𝑇𝑝𝑏𝑖=𝑠6𝑙𝑏𝑖(7.14) with𝑠1𝑑𝑖𝑏𝑗=𝑘𝑑𝑖𝑏𝑗+𝛾𝑐𝛽𝑑𝑡𝑑𝑖𝑏𝑗+1𝛽𝑑𝑡2𝑚𝑑𝑖𝑏𝑗,𝑠(7.15)2𝑑𝑖𝑝𝑏𝑗=𝛾𝛽𝑑𝑡𝑑𝑖𝑝𝑏𝑗,𝑠(7.16)3𝑑𝑖=𝜆1𝑁𝑚1𝑓𝑑𝑖+𝐸𝑑𝑖𝑟𝑑𝑖+𝑐𝑑𝑖𝑏𝑗𝑏𝑥𝑏𝑗+𝑚𝑑𝑖𝑏𝑗𝑎𝑥𝑏𝑗𝑑𝑖𝑝𝑏𝑗𝑏𝜏p𝑏𝑗,𝑠(7.17)4𝑙𝛾=𝑛𝛽𝑑𝑡𝑙𝜆𝑚𝑠1,𝑠(7.18)5𝑙𝑝=1𝛽𝑑𝑡2𝜆𝑠1𝛿𝑙𝑝𝛾𝜇𝛽𝑑𝑡𝑙𝑝,𝑠(7.19)6𝑙𝑏𝑖=𝑛𝑙𝜆𝑚𝑠1̇𝑥𝑏𝑖̇𝑤𝑏𝑖+𝜇𝑙𝑝𝜏𝑝𝑏𝑖𝜆𝑠1̇𝜏𝑙𝑏𝑖𝑛𝑙𝜆𝑚𝑠1𝑏𝑥𝑏𝑖𝜇𝑙𝑝𝑏𝜏𝑝𝑏𝑖+𝜆𝑠1𝑎𝜏𝑙𝑏𝑖.(7.20) One can then condense 𝑑𝑇𝑝𝑏𝑖 out of the above system of equations:𝑑𝑇𝑝𝑏𝑗=𝑠51𝑝𝑙𝑠6𝑙𝑏𝑗𝑠4𝑙𝑑𝑥𝑏𝑗𝑠,(7.21)1𝑑𝑖𝑏𝑗+𝑠2𝑑𝑖𝑝𝑏𝑗𝑠51𝑝𝑙𝑠4𝑙𝑑𝑥𝑏𝑗=𝑠3𝑑𝑖+𝑠2𝑑𝑖𝑝𝑏𝑗𝑠51𝑝𝑙𝑠6𝑙𝑏𝑗.(7.22) Equation (7.22) is forced into “Newmark form” as𝑘𝑑𝑖𝑏𝑗+𝑘𝑑𝑖𝑏𝑗+𝛾𝑐𝛽𝑑𝑡𝑑𝑖𝑏𝑗+1𝛽𝑑𝑡2𝑚𝑑𝑖𝑏𝑗𝑑𝑥𝑏𝑗=𝑓𝑑𝑖𝑟𝑑𝑖+𝑐𝑑𝑖𝑏𝑗𝑏𝑥𝑏𝑗+𝑚𝑑𝑖𝑏𝑗𝑎𝑥𝑏𝑗(7.23) with𝑘𝑑𝑖𝑏𝑗=𝑠2𝑑𝑖𝑝𝑏𝑗𝑠51𝑝𝑙𝑠4𝑙,𝑓𝑑𝑖=𝜆1𝑁𝑚1𝑓𝑑𝑖+𝐸𝑑𝑖𝑑𝑖𝑝𝑏𝑗𝑏𝜏𝑝𝑏𝑗+𝑠2𝑑𝑖𝑝𝑏𝑗𝑠51𝑝𝑙𝑠6𝑙𝑏𝑗(7.24)𝑘𝑑𝑖𝑏𝑗 and 𝑓𝑑𝑖 both depend on 𝑑𝑡, 𝛽, and 𝛾: the symbol 𝑘𝑑𝑖𝑏𝑗 was chosen to indicate that the matrix is handled by the Newmark-𝛽 solver in the same way as a stiffness. However, this term cannot be interpreted physically as a stiffness.

7.6. Spacial Discretisation

The consistent discretisation by Galerkin finite elements of (7.23) leads to𝐾𝑛𝑚=𝑁𝑑𝑖𝑛𝑘𝑑𝑖𝑏𝑗𝑁𝑏𝑗𝑚,𝐹𝑛=𝑁𝑑𝑖𝑛𝑓𝑑𝑖.(7.25)𝐾𝑛𝑚 and 𝐹𝑛 are typically computed by Gauss quadrature. Note that no space derivative is present in 𝑘𝑑𝑖𝑏𝑗, so no partial integration or Gauss quadrature with curvature shape function appears. One can hence simplify the expression of the element matrix to𝐾𝑛𝑚=𝑁𝑑𝑖𝑛𝑘𝑑𝑖𝑏𝑁𝑏𝑖𝑚(7.26) which means “same quadrature as for a mass matrix”.

7.7. Implementation

In nonlinear FEM code, incremental matrices and vectors are computed by Gauss quadrature. The Gauss quadrature involves shape functions, tensors that are local, continuous versions of the stiffness, damping and mass matrices, and the force imbalance vector. For example, for the drag damping of a beam element, the tensor relates a vector whose components are increments in velocities in three directions, to another vector whose components are increments in forces per unit length in three directions.

Within an iteration, the linear solver provides incremental nodal positions, velocities and accelerations for the model. These are disassembled and provided to the elements. The elements compute positions, velocities, and accelerations (and more) in a corotated reference system at Gauss points. The resulting values are handed to the VIV-Gauss point procedure.

The axial velocities are discarded. The procedure scales the provided values using (3.2) and (3.3).

Having stored the previous approximation of the scaled position, the procedure determines the position increment 𝑑𝑥𝑏𝑗, and then uses (7.21) to obtain the Laguerre coefficient increment 𝑑𝑇𝑝𝑏𝑗. From there, (7.12) and (7.13) are used to compute 𝑑𝜏𝑝𝑏𝑗 and 𝑑̇𝜏𝑝𝑏𝑗. The values of 𝑇𝑝𝑏𝑗, 𝜏𝑝𝑏𝑗, and ̇𝜏𝑝𝑏𝑗 are updated from previously stored values. 𝜏𝑝𝑏𝑗 is then used to evaluate 𝑓𝑑𝑖 and its derivative with respect to 𝜏𝑝𝑏𝑗. These are scaled back, and Froude-Krylov forces are added, leading to 𝑘𝑑𝑖𝑏𝑗 and 𝑓𝑑𝑖.

The above matrix and vector are padded with zeros to indicate zero force in the axial direction and zero torque.

The condensation of a larger system of time-discretised equations introduces some inelegant features compared to standard dynamic FEM: the VIV-Gauss point must be provided with 𝛽, 𝛾, and 𝑑𝑡 and a flag showing whether a call is made at a step or within a refinement iteration.

Note that when the VIV model is integrated in the dynamic FEM computation, in the way detailed in this section, the recursive Laguerre filter presented in Section 4.6 is not used. In this work, the filter was used only to process velocities time series and product Laguerre coefficients for the training of the rotatron model. The filter provides a very efficient way to process whole time series.

8. Results

8.1. Training

The Norwegian Deepwater Program was a research effort in which reduced scale tests were carried out on long, flexible riser models, subject to uniform or sheared current [35]. The displacement histories thus acquired at 19 points along the riser model were later prescribed on short stiff cylinders, and the hydrodynamic forces acting on the cylinders directly measured [36, 37]. Some data from [36, 37] is used in this work. It consists of the displacements at 19 points along the NDP riser model, for 3 current profiles (Table 1), so a total of 57 short cylinder test runs. For each of the 57 runs, 100 instants are randomly selected, yielding a training set of the rotatron with 5700 “points”. Each “point” consists of two sets of 𝑛=30 Laguerre coefficients and the two components of the corresponding force (Figure 11).

The rotatron was trained using 𝑛=30 Laguerre polynomials, 200 neurons in the hidden layer, and 50 to 1000 iterations of the conjugate gradient optimization algorithm. Other settings have been studied. No precision improvement was obtained from increasing the number of neurons or from increasing the number of iterations. This may suggest that the optimization algorithm converges on a local minima. If that is the case, improvement would require the use of an optimization algorithm better suited to finding absolute minima in a “jungle” of local ones.

Figure 9 shows how the rotatron predicts the forces for the trajectories in the above-mentioned 57 runs of short cylinder tests. The predicted forces are compared to the forces acquired experimentally. The model's ability to predict these forces seems to be good, although we lack a good criteria to judge that yet.

Figure 10 provides a visualization of the different steps of the modelization process, and is hence a useful diagnostic tool. It shows stippled black line: the trajectory for which a force prediction is wanted,smooth black line: the Laguerre approximation to the above trajectory, used to enter the rotatron,stippled black arrow: the measured force for the above trajectory,smooth black arrow: the predicted force for the above trajectory,green: neighboring (in the sense of the rotatron distance, Equation (6.7)) trajectories from experiments, used in the training set (and corresponding Laguerre approximation, experimentally measured force and predicted force),red: same as the above after rotation and/or mirroring.

Figure 10 gives an indication of the quality of the Laguerre approximation, the adequacy of the training set for the trajectory at hand, the presence of contradictions in the training set near the trajectory at hand, the quality of the fit of the rotatron to the training data, and finally the quality of the interpolation between training points.

8.2. Dynamic Analysis of a Flexible Riser

VIV depends not only on current velocities, but also on the type of slender system they act upon. Tension, stiffness, damping, length, and boundary conditions affect the vibration and hence the velocity trajectories that appear in the vibration. Hence the database used in Section 8.1 to train the rotatron is specialized, not only to a few Reynolds numbers for the current velocity, but also to some extend to the particular model used in the NDP program. No study was carried out in this work to determine how much change in the structure this specialized rotatron can accommodate.

Hence, in order to test the performance of the present VIV model within a dynamic analysis, the simplest case was considered: the riser model used in the NDP testing program ([35], characteristics in Table 2) was modeled. The present method was used to analyze the test condition TN2370 (Re[0,24300]). Numerical results were compared with those obtained experimentally on the NDP model. A laptop, using one core of a dual core processor, took 20 to 50 seconds to compute one second of riser response.

In Figures 12 and 13, the horizontal axis is NDP-laboratory time and the vertical axis is the unscaled length along the riser. The upper subplot shows the response in line (IL) with the flow, the lower subplot the cross-flow (CF) response. The color codes the displacements, with the same color scale used in Figures 12 to 14.

Figures 12 and 13 are for test TN2370 (Re[0,24300]). The dynamic simulations captures the frequency doubling between CF and IL, as well as the instationary nature of the vibrations. Frequency and amplitude are adequately captured. The 6th mode's dominance of the CF vibration is correctly captured. The dynamic analysis assumes constant tension. By contrast, some small tension modulations seem to displace the position of the lower vibration node in the test. In test 2370, there is a marked tendency for CF vibrations to propagate downwards. In the analysis results, CF waves are of a more static character. The IL vibration in the analysis occurs at a higher mode (11) than in the test (9). This is best seen by counting red dots along diagonals, for example, starting from coordinate 38 m and time 17.9 s in Figure 12. As a consequence, and since, for a given propagation celerity, wavelength and period are related, the phase drift between IL and CF is of opposite sign in the analysis and the test.

It proved impossible to reproduce tests 2340 and 2030 in the same manner. The analysis quickly ends with IL and CF vibrations occurring at the same frequency, and in Figure 14, with in-line motions dominating.

9. Discussion

9.1. Force Prediction and Reynolds Interpolation

Figure 9 shows that given the velocities in the training data as input, the model allows to reproduce the forces in the training data, based on the velocity in the training data.

The same exercise is carried out with trajectories from another test, N2430 (not to be confused with TN2340), which was carried out in shear current (Re=0 to Re=40500). In comparison, the highest current velocity appearing in the training set is Re=24300. In N2430, the forces are fairly well predicted at the lower current velocities (Figure 15), while the predictions are very poor at higher velocities.

This illustrates that the present model provides no mechanism to use forces from a test at a given Reynolds number, to predict forces in a situation with another Reynolds number (due to higher current velocity e.g.,). In contrast, a simplified harmonic method like VIVANA enters a “table” with values of reduced amplitude and frequency. The tables are assumed to represent VIV behavior over a relevant range of Reynolds numbers. Because the transformation from added mass and excitation coefficients into forces makes use of the current speed (which is proportional to the Reynolds number), this effectively provides VIVANA with a mechanism for interpolation over Reynolds number. The present model, on the other hand, uses a compressed form (Laguerre coefficient) of recent velocity history, scaled to Reynolds-like numbers. The scaling and data compression scheme used in the present model does not separate the classical Reynolds number as a separate variable. Instead, it is embedded in the compressed data as the average of the recent scaled in line velocity already provided to the rotatron.

It remains an open question whether one should at all try to add such a Reynolds interpolation scheme in models of the present type. In a range of Reynolds number values in which the wake flow remains similar, it is reasonable to assume that viscous forces increase linearly with velocity and inertial forces increase quadratically. For harmonic motions, the component of the force in phase with acceleration can be attributed to an inertial force, and the component in phase with velocity, to a viscous force. For more general motions, “phase” becomes a problematic concept, making it difficult to separate the force into viscous and inertial components that one can scale each according to their own law.

9.2. Dynamic Analysis

The comparison of Figures 12 and 13 is encouraging: when simulating exactly the system from which trajectories were acquired experimentally in the first place, the simulation comes strikingly close to the experimental results. It is understood that a simulation is good if it captures the statistical properties of the real response. With VIV being chaotic, there is no hope to reproduce exactly any given realization of the response.

Figure 14 shows the result of a simulation of a case which is also directly represented in the training data. The simulation fails, in the sense that IL and CF vibrations are simulated to occur at the same frequency. Two explanations are proposed.

The quality of the force prediction by the rotatron, illustrated in Figure 9, was declared “satisfactory”, but this constitutes only an observation that the fitting procedure is operating. On what criteria should one judge that the fit is adequate? When training was carried out, the average norm of the difference between the training force and the predicted force was found to be about 30% of the average norm. For a control group composed of data points from the same experimental database, but not used in training the perceptron, the same ratio was about 40%. These numbers are high and can be reduced to some extend by increasing the number of hidden layers and of training iterations, but this was not found to yield better simulations. Further specializing the perceptron (by training it with data from test 2030) did not lead to successful simulations. One study which might help to understand the observations would be to find the corrective forces that need to be added to the forces predicted by the model, to force the simulation to track the motions observed during test 2030. This can be achieved using dynamic inverse FEM analysis [3840].

Even if the above corrective forces were strictly zero, so that the model was accurately predicting forces for the riser motion from test 2030, this would not be sufficient to ensure that the simulation adequately mimics test 2030. A given trajectory could still have very different stability properties in the physical system and in the simulation. It could be that in the physical system, the trajectories follow the “bottom of a valley” while the model renders it as the “crest of a mountain”. A measure of stability for this is the Lyapunov exponent [41, 42]. Procedures exist to compute Lyapunov exponents from experimental data, and this should be compared to Lyapunov exponents for the simulation.

9.3. Influence of 𝑡𝑤

One issue that was explored was the adequacy of 𝑡𝑤: the rotatron was trained with 𝑡𝑤=5105. This corresponds roughly to 1/4 of a cross-flow oscillation period in reduced scale for test 2370, and to a smaller fraction for other tests at lower current velocity. It could be argued that this fraction becoming too small could be the cause of the failure of analyses at lower velocities (Figure 14). The rotatron was hence trained again using a suitably increased value of 𝑡𝑤. This was done twice, once with the full training set and increasing the number of Laguerre polynomials to 𝑛=60, and once with only the experimental data from lower currents and an unchanged number of polynomials (𝑛=30). Neither rotatrons allowed to perform a successful simulation for the lower current velocities.

9.4. Tension

In Figure 12, one can observe modulations of the positions of the vibration nodes (in particular on the cross-flow graph). One possible explanation for this would be that in the physical system, the tension is modulated by the vibration. This effect, if present, is not captured by the numerical model, which assumes constant tension. The method presented here is designed for use in a nonlinear analysis. However, in this research, time was saved by using a simpler linear structural model.

10. Conclusions

A model for the prediction of VIV forces given the history of velocity of a cylindrical cross section relative to the undisturbed fluid has been developed. The model is closely related to Wiener-Laguerre filters: the recent history of velocity is represented by the coefficients of a Laguerre polynomial series. These coefficients are then used to enter a memory-less nonlinear interpolation function, in this case, a custom made neural network in which some relevant symmetry properties were “hard-wired”. The neural network was trained by using forces and displacements obtained in irregular forced motion tests on a short cylinder.

The proposed model operates in the time domain, making it well suited for integration into fully nonlinear analyses with unsteady currents. It further deals with in-line and cross-flow vibrations as one inseparable issue, which is arguably a necessity to improve on existing VIV models.

The model could provide a “good” reproduction of the forces in the training test, as well as a “good” prediction of forces for “comparable” trajectories. Where the model was queried with trajectories very different from those present in the training data, the model gave very poor results—as should be expected. Due to the limited amount of experimental data available, the present model remains quite specialized to a limited number of situations. What remains unknown at this stage is the size of the training set, and of the neural network model, necessary to create a model with some ambition of generality.

In some dynamic analyses of laboratory tests (NDP TN2030 and TN2430) with a long flexible riser, the numerical solution fell into an unphysical mode of vibration with the same frequency for in-line and cross flow vibration. On the other hand, for another case (TN2470), some unusually fine details were captured by the numerical model. From this, it is concluded that the concept has merit and deserves to be pursued.

At the same time, the present paper leaves several important questions unanswered.

The proposed method assumes that hydrodynamic forces can be predicted by the recent history of velocity of the pipe relative to the surrounding fluid. Fractal dimensional analysis in Section 6.3.2 suggests that this is, to some extent, the case. On the other hand, experimental evidence clearly shows that if a cylinder is forced twice to follow the same trajectory, the histories of hydrodynamic forces will not be exactly the same. Will adding a random force process to the output of the model be of any value in reproducing the dynamics behavior of a pipeline? If so, how should such a random process be studied and modeled?

What is the information content of the relation between velocities and forces? As the size of the experimental database is increased, will data points confirm each other and be interpolated with a relatively simple model (a rotatron with few tunable parameters), or will a landscape of ever increasing complexity appear, defeating the present approach?

The rotatron was found to be able to predict forces “reasonably” well, yet, when included in a dynamic model, the ability to reproduce the dynamics of the system was limited. Can the measure of force agreement, used in training the rotatron be improved? Would it be useful, as suggested to the author by Dr. Halvor Lie to emphasize a good prediction of the power of the force? More generally, how does one ensure that an approximate model reasonably reproduces the statistics of the response of the original nonlinear dynamics system?


A. Conventions for Indexed Notations

The following conventions for indexed notations are used in this text. (1)By default, where an index appears more than once in a combination of products and/or divisions, a summation over the index is implied. If that index has a continuous range, then the “sum” is an integration over the range. Points 2, 3, and 4 specify exceptions to this rule. (2)Point 1 notwithstanding, if an equation is preceded by the symbol , followed by a list of indexes, then the listed indexes are not summed over. (3)Point 1 notwithstanding, if within an equation, there is a combination of products and/or divisions within which an index appears only once, then no sum over that index is carried out in the whole equation. (In any situation other than a simple term in the left-hand side, readability should be improved by using the symbol .)(4) Point 1 notwithstanding, if an index appears within an input to a function, and the output of the function is multiplied or divided by one or several terms that have the same index, then no sum within the input to the function is carried out on that index. (5)If an index of an argument to a function is within brackets, then the whole range of index values is used as input to one function evaluation. For example, 𝜎𝑖(𝑦[𝑗]𝑘) refers to the evaluation at multiple locations (𝑘) of a vector-valued (𝑖) function of a vector (𝑗). (6)When the output of the function is shorthanded without explicitly writing its input, then the indexes of the input that are not within brackets are added to the indexes of the function. For example, 𝜎𝑖(𝑦[𝑗]𝑘) can be shorthanded 𝜎𝑖𝑘. (7)Derivatives of a function are noted with only the bracketed indexes of the input appearing under the fraction: 𝜕𝜎𝑖/𝜕𝑦𝑗. To refer to the value of that derivative for input 𝑘, one writes 𝜕𝜎𝑖𝑘/𝜕𝑦𝑗.

B. Rotatron Gradients

The derivative of the force predicted by the rotatron, with respect to the Laguerre coefficients is needed in Section 7. With references to (5.3) to (5.8) that describe the rotatron, we can write𝜕𝑓𝑖𝑛𝜕̇𝜏𝑗𝑙=𝜕𝑓𝑖𝑛𝜕𝜎𝑘𝜕𝜎𝑖𝑘𝑛𝜕𝑦𝑗𝜕𝑦𝑗𝑘𝑛𝜕̇𝜏𝑙=𝑉𝑘𝜕𝜎𝑖𝑘𝑛𝜕𝑦𝑗𝑀𝑘𝑙(B.1) with𝜕𝜎𝑖𝑘𝑛𝜕𝑦𝑗1=||𝑦[𝑗]𝑘𝑛||3||𝑦[𝑗]𝑘𝑛||𝛼𝑘+12×𝛼𝑘𝑦𝑖𝑘𝑛𝑦𝑗𝑘𝑛||𝑦[𝑗]𝑘𝑛||𝛼𝑘+(1)𝛿𝑖𝑗𝑦¬𝑖𝑘𝑛𝑦¬𝑗𝑘𝑛||𝑦[𝑗]𝑘𝑛||𝛼𝑘.+1(B.2) Here index 𝑖 ranges over two values (for two directions orthogonal to the cylinder), and ¬𝑖 is the other direction than 𝑖.

The gradients of the rotatron with respect to its coefficients are also needed in order to compute the gradient of the target function with respect to the parameters 𝑉𝑘, 𝑈𝑘, and 𝑀𝑘𝑙:𝜕𝑓𝑖𝑛𝜕𝑉𝑙=𝜎𝑖𝑙𝑛,𝜕𝑓𝑖𝑛𝜕𝑀𝑘𝑙=𝜕𝑓𝑖𝑛𝜕𝜎𝑘𝜕𝜎𝑖𝑘𝑛𝜕𝑦𝑗𝜕𝑦𝑗𝑘𝑛𝜕𝑀𝑘𝑙=𝑉𝑘𝜕𝜎𝑖𝑘𝑛𝜕𝑦𝑗̇𝜏𝑗𝑙𝑛,𝜕𝑓𝑖𝑛𝜕𝑈𝑙=𝑉𝑘𝜕𝜎𝑖𝑘𝑛𝜕𝑈𝑙(B.3) with𝜕𝜎𝑖𝑘𝑛𝜕𝑈𝑙=𝛿𝑘𝑙𝑒𝑈𝑙𝑦𝑖𝑘𝑛||𝑦[𝑗]𝑘𝑛||𝛼𝑘1||𝑦log[𝑗]𝑘𝑛||||𝑦[𝑗]𝑘𝑛||𝛼𝑘+12.(B.4)

C. Inverse of 𝑠5

The inverse of 𝑠5 (Equation (7.19)), where 𝑠5 is of the form𝑠5𝑖𝑗=𝛼𝑇𝑖j+𝛽𝛿𝑖𝑗(C.1) with𝑇𝑖𝑗=1𝑗𝑖0𝑗>𝑖,(C.2) can be verified to be lower triangular banded, with terms on diagonal 𝑖 equal to𝑄1=1,𝑄𝛼+𝛽𝑖=𝛼𝛽𝑖2(𝛼+𝛽)𝑖,𝑖{2𝑛}.(C.3)


The author is indebted to Ida Aglen, Celeste Barnardo, Trygve Kristiansen, Carl Martin Larsen, Halvor Lie, Elizabeth Passano, Thomas Sauder, Wu Jie, Yin Decao and an anonymous reviewer for inspiring discussions on the subject of VIV, that either led to or helped the present research and writing. Thanks are extended to CeSOS (Centre of Excellence for Ships and Offshore Structures) at NTNU (Norwegian University of Science and Technology) and to MARINTEK for sponsoring this research. Very special thanks to Professor Carl Martin Larsen of CeSOS who, faced with a barrage of dubious ideas from the author, responded by making this research possible.