# A Wiener-Laguerre Model of VIV Forces Given Recent Cylinder Velocities

Authors: Philippe Mainçon
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

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.


---

## Body

## 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 [14–17]. 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.


## 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 [14–17]. 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 timestw, where the value of tw still is an object of debate. tw is likely to be case dependent. The current will transport (convect) away vortices so that they quickly loose significance. The time tw should then be of the order of D/U where D is the cross section diameter and U the current velocity. In contrast, if the cylinder is oscillating in still water, it will be traveling in its own wake, and tw should be related to the rate of diffusion and/or viscous dissipation of vortices, which is likely to result in much higher values of tw. 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 tw, 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 therecent 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 Section7.


### 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 (Section3.2).(3)
The trajectory (again: the recent histories of bothx and y 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 returnsx and y 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 (Section3.2).(6)
The Froude-Krylov forces, which depend on the acceleration of the undisturbed flow, are added (Section3.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 functionsintroduced 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 [20–22]. In the present work, Laguerre filtering is presented without making use of the vocabulary of cybernetics. In particular, the z-transform is not introduced here.


## 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 timestw, where the value of tw still is an object of debate. tw is likely to be case dependent. The current will transport (convect) away vortices so that they quickly loose significance. The time tw should then be of the order of D/U where D is the cross section diameter and U the current velocity. In contrast, if the cylinder is oscillating in still water, it will be traveling in its own wake, and tw should be related to the rate of diffusion and/or viscous dissipation of vortices, which is likely to result in much higher values of tw. 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 tw, 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 therecent 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 Section7.


## 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 (Section3.2).(3)
The trajectory (again: the recent histories of bothx and y 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 returnsx and y 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 (Section3.2).(6)
The Froude-Krylov forces, which depend on the acceleration of the undisturbed flow, are added (Section3.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 functionsintroduced 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 [20–22]. In the present work, Laguerre filtering is presented without making use of the vocabulary of cybernetics. In particular, the z-transform is not introduced here.


## 3. Exploiting Similitudes

### 3.1. Froude-Krylov Forces

This section gives the justification for point 6 of Section2.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 Sections4 and 5 deal only with the first of the above two terms.


### 3.2. Scaling

This section details how points 2 and 5 of Section2.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 D, 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 ρo, νo, and Do. The choice of ρo, νo, and Do 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λm, λs, and λkg, which are the scaling factors for the basic units of distance, time, and mass. Solving the system yields(3.1)λm=1D,λs=νD2,λkg=1ρD3.
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:(3.2)λms-1=Dν,(3.3)λms-2=D3ν2,(3.4)λNm-1=Dρν2.
The choice Do=1 [m] hence implies that scaled displacements can be considered to have “1 diameter” as unit. Similarly, the choices Do=1 [m] and νo=1 [m2/s] together imply that scaled velocities are expressed as Reynolds numbers since the scaled velocity is calculated as Dv/ν where v 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 D/ν. 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 D is scaled to values that depend of ρ, ν, and D. In particular, (3.3) shows that all accelerations, including the acceleration of gravity g, are scaled with a factor proportional to D3/ν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 go.


## 3.1. Froude-Krylov Forces

This section gives the justification for point 6 of Section2.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 Sections4 and 5 deal only with the first of the above two terms.


## 3.2. Scaling

This section details how points 2 and 5 of Section2.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 D, 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 ρo, νo, and Do. The choice of ρo, νo, and Do 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λm, λs, and λkg, which are the scaling factors for the basic units of distance, time, and mass. Solving the system yields(3.1)λm=1D,λs=νD2,λkg=1ρD3.
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:(3.2)λms-1=Dν,(3.3)λms-2=D3ν2,(3.4)λNm-1=Dρν2.
The choice Do=1 [m] hence implies that scaled displacements can be considered to have “1 diameter” as unit. Similarly, the choices Do=1 [m] and νo=1 [m2/s] together imply that scaled velocities are expressed as Reynolds numbers since the scaled velocity is calculated as Dv/ν where v 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 D/ν. 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 D is scaled to values that depend of ρ, ν, and D. In particular, (3.3) shows that all accelerations, including the acceleration of gravity g, are scaled with a factor proportional to D3/ν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 go.


## 4. Characterization of Trajectory

### 4.1. Foreword

This section details how point 3 in Section2.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 (Section4.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 (Section4.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 (Figure1, top) of degree i-1 can be defined by its Rodrigues formula [24](4.1)Li(x)≡exi!didxi(xie-x).
Laguerre polynomials verify the orthonormality property(4.2)∫0∞Li(x)Lj(x)e-xdx=δij.
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)W(t)≡et/twtw,t∈R-,
where the interpretation of tw 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(4.4)W¯¯(t1,t2)=δ(t1,t2)W(t1)
and a dot product in a suitable space of real valued functions(4.5)f¯T∘g¯≡∫-∞0f(t)g(t)dt
with the canonical norm(4.6)|f¯|≡f¯T∘f¯.
Further we introduce the base(4.7)L¯¯(t,i)≡Li(-ttw),t∈R-,i∈{1,…,n}.
Equation (4.2) can be rewritten in matrix notation as(4.8)I¯¯=L¯¯T∘W¯¯∘L¯¯,
where I¯¯ is the n×n identity matrix. It is useful to introduce the weighted norm or w-norm(4.9)|f¯|w≡f¯T∘W¯¯∘f¯.
Note that since 𝒲(t) is of dimension [1/s],  |f¯|w is of the same dimension as f¯. So when taking f¯ as a scaled velocity, |f¯|w is an ilr-Reynolds number.Figure 1
Laguerre polynomials (top), Laguerre analysis functions (middle), and weight function (bottom). The convolution of a signal by the Laguerre analysis functions yields Laguerre coefficients (analysis). The linear combinations of the Laguerre polynomials weighted by the coefficients give an approximation of the original signal (synthesis), with a quality that decreases towards the past in a way related to the weight function.


### 4.3. Analysis and Synthesis

For a historyv¯(t) of either the x or y component of the velocity, we seek the vector of “Laguerre coefficients” τ¯ with which to combine the columns L¯¯, that minimize the weighted error J defined as(4.10)J=12|v¯-L¯¯⋅τ¯|w2=12(v¯-L¯¯⋅τ¯)T∘W¯¯∘(v¯-L¯¯⋅τ¯).
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)∂J∂τ¯=L¯¯T∘W¯¯∘L¯¯⋅τ¯-L¯¯T∘W¯¯∘v¯,
which implies(4.12)τ¯=(L¯¯T∘W¯¯∘L¯¯)-1⋅L¯¯T∘W¯¯∘v¯(4.13)=L¯¯T∘W¯¯∘v¯,(4.14)τ¯=D¯¯T∘v¯
with(4.15)D¯¯≡W¯¯∘L¯¯.
The “Laguerre analysis functions” D¯¯ (Figure 1, middle) are by definition equal to(4.16)D¯¯(t,i)=Di(-ttw)=Li(-ttw)et/twtw.
The Laguerre analysis functions D¯¯ 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 analysisfunction), 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)F¯¯(t,i)≡Fi(-ttw)(4.18)≡Li(-ttw)et/(2tw)tw
or in matrix notation as(4.19)F¯¯=W¯¯∘L¯¯,
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 as(4.20)lim⁡n→∞|F¯¯⋅F¯¯T∘f¯-f¯|=0.
This can be used to obtain a result on the convergence of series of Laguerre polynomials. We introduce the change of variables(4.21)f¯=W¯¯∘g¯
so that(4.22)|F¯¯⋅F¯¯T∘f¯-f¯|=|F¯¯⋅D¯¯T∘g¯-W¯¯∘g¯|=|W¯¯∘(L¯¯⋅D¯¯T∘g¯-g¯)|=|L¯¯⋅D¯¯T∘g¯-g¯|w.
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 w-norm is equivalent to convergence almost everywhere. So under some conditions of continuity on g¯, the series of Laguerre polynomials obtained using D¯¯ as analysis functions converges almost everywhere towards g¯ in any finite interval.Figure2 illustrates how Laguerre coefficients indeed provide a “summary” of the trajectoryFigure 2
Example of Laguerre approximation for two components of a velocity history (arbitrary scaling). The red curve is the original cyclic signal. Each black dot marks a present time and the black curves are the corresponding Laguerre approximations for the recent past.


### 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 n×N floating point operations (flops), where n is the number of Laguerre polynomial used and N the number of time steps that the analysis functions take to decay to a negligible value. Further, for each Gauss point, 2N 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 n×n flops, which is advantageous because n≪N. In this section, it is shown that τ¯ verifies a differential equation driven by the history v¯ 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(4.23)∂τi∂t=∫0+∞e-θLi(θ)∂v∂t(t-twθ)dθ=-1tw∫0+∞e-θLi(θ)∂v∂θ(t-twθ)dθ.
Multiplying by tw and integrating by parts yields(4.24)tw∂τi∂t=-[e-θLi(θ)v(t-twθ)]+∫0+∞[-e-θLi(θ)+e-θ∂∂θLi(θ)]v(t-twθ)dθ.
A property of Laguerre polynomials is [24](4.25)∂∂θLi(θ)=-Li-1(1)(θ)=-∑j=0i-1Lj(θ),
where ℒi(1)(θ) is a generalized Laguerre polynomial. Hence we can write(4.26)tw∂τi∂t=Li(0)v(t)-τi-∫0+∞e-θ∑j=1i-1Lj(θ)v(t-twθ)dθ=v(t)-τi-∑j=1i-1τj=v(t)-∑j=1iτj.
which is of the form(4.27)∂τ¯∂t(t)=μ¯¯⋅τ¯(t)+n¯v(t)
with(4.28)μij={-1twj≤i0j>i,ni=1tw.
Equation (4.27) shows that at any time t, 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 signalsL¯¯·τ¯ 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” (Section5.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(4.29)vj=v(t0+jdt).
We seek the values of the Laguerre coefficients at the same intervals(4.30)τ¯j=τ¯(t0+jdt).
The vector τ¯j (the list of the coefficients for all Laguerre polynomials, taken at step j) must not be confused with scalar τi (the coefficient for the Laguerre polynomial of degree i). We choose t0 such that t0+jdt=0, and we approximate v by a function that is linear over the interval [0,dt]. Equation (4.27) becomes(4.31)∂τ¯∂t(t)=μ¯¯⋅τ¯(t)+α¯+β¯t
with(4.32)α¯=n¯v(0),β¯=n¯v(dt)-v(0)dt.
This new differential equation can be solved exactly: we seek a solution of the form(4.33)τ¯(t)=exp⁡(μ¯¯t)⋅a¯+b¯t+c¯
over the interval. Here exp(μ¯¯t) stands for a matrix exponential. Replacing this expression into (4.31), noting that(4.34)∂∂texp⁡(μ¯¯t)=μ¯¯⋅exp⁡(μ¯¯t),exp⁡(0¯¯)=I¯¯,
and identifying the constant and linear terms and enforcing the initial value leads to(4.35)b¯=-μ¯¯-1⋅β¯,c¯=-μ¯¯-2⋅β¯-μ¯¯-1⋅α¯,a¯=τ¯(0)+μ¯¯-2⋅β¯+μ¯¯-1⋅α¯.
Replacing these expressions in (4.33) at t=dt, a tedious but straightforward computation yields the recursive filter(4.36)τ¯j+1=M¯¯⋅τ¯j+V¯1⋅vj+V¯2⋅vj+1
with(4.37)M¯¯=exp⁡(m¯¯dt),μ¯1=μ¯¯-1⋅n¯,μ¯2=μ¯¯-2⋅n¯1dt,V¯1=M¯¯⋅(μ¯1-μ¯2)+μ¯2,V¯2=M¯¯⋅μ¯2-μ¯1-μ¯2.


## 4.1. Foreword

This section details how point 3 in Section2.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 (Section4.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 (Section4.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 (Figure1, top) of degree i-1 can be defined by its Rodrigues formula [24](4.1)Li(x)≡exi!didxi(xie-x).
Laguerre polynomials verify the orthonormality property(4.2)∫0∞Li(x)Lj(x)e-xdx=δij.
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)W(t)≡et/twtw,t∈R-,
where the interpretation of tw 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(4.4)W¯¯(t1,t2)=δ(t1,t2)W(t1)
and a dot product in a suitable space of real valued functions(4.5)f¯T∘g¯≡∫-∞0f(t)g(t)dt
with the canonical norm(4.6)|f¯|≡f¯T∘f¯.
Further we introduce the base(4.7)L¯¯(t,i)≡Li(-ttw),t∈R-,i∈{1,…,n}.
Equation (4.2) can be rewritten in matrix notation as(4.8)I¯¯=L¯¯T∘W¯¯∘L¯¯,
where I¯¯ is the n×n identity matrix. It is useful to introduce the weighted norm or w-norm(4.9)|f¯|w≡f¯T∘W¯¯∘f¯.
Note that since 𝒲(t) is of dimension [1/s],  |f¯|w is of the same dimension as f¯. So when taking f¯ as a scaled velocity, |f¯|w is an ilr-Reynolds number.Figure 1
Laguerre polynomials (top), Laguerre analysis functions (middle), and weight function (bottom). The convolution of a signal by the Laguerre analysis functions yields Laguerre coefficients (analysis). The linear combinations of the Laguerre polynomials weighted by the coefficients give an approximation of the original signal (synthesis), with a quality that decreases towards the past in a way related to the weight function.


## 4.3. Analysis and Synthesis

For a historyv¯(t) of either the x or y component of the velocity, we seek the vector of “Laguerre coefficients” τ¯ with which to combine the columns L¯¯, that minimize the weighted error J defined as(4.10)J=12|v¯-L¯¯⋅τ¯|w2=12(v¯-L¯¯⋅τ¯)T∘W¯¯∘(v¯-L¯¯⋅τ¯).
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)∂J∂τ¯=L¯¯T∘W¯¯∘L¯¯⋅τ¯-L¯¯T∘W¯¯∘v¯,
which implies(4.12)τ¯=(L¯¯T∘W¯¯∘L¯¯)-1⋅L¯¯T∘W¯¯∘v¯(4.13)=L¯¯T∘W¯¯∘v¯,(4.14)τ¯=D¯¯T∘v¯
with(4.15)D¯¯≡W¯¯∘L¯¯.
The “Laguerre analysis functions” D¯¯ (Figure 1, middle) are by definition equal to(4.16)D¯¯(t,i)=Di(-ttw)=Li(-ttw)et/twtw.
The Laguerre analysis functions D¯¯ 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 analysisfunction), 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)F¯¯(t,i)≡Fi(-ttw)(4.18)≡Li(-ttw)et/(2tw)tw
or in matrix notation as(4.19)F¯¯=W¯¯∘L¯¯,
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 as(4.20)lim⁡n→∞|F¯¯⋅F¯¯T∘f¯-f¯|=0.
This can be used to obtain a result on the convergence of series of Laguerre polynomials. We introduce the change of variables(4.21)f¯=W¯¯∘g¯
so that(4.22)|F¯¯⋅F¯¯T∘f¯-f¯|=|F¯¯⋅D¯¯T∘g¯-W¯¯∘g¯|=|W¯¯∘(L¯¯⋅D¯¯T∘g¯-g¯)|=|L¯¯⋅D¯¯T∘g¯-g¯|w.
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 w-norm is equivalent to convergence almost everywhere. So under some conditions of continuity on g¯, the series of Laguerre polynomials obtained using D¯¯ as analysis functions converges almost everywhere towards g¯ in any finite interval.Figure2 illustrates how Laguerre coefficients indeed provide a “summary” of the trajectoryFigure 2
Example of Laguerre approximation for two components of a velocity history (arbitrary scaling). The red curve is the original cyclic signal. Each black dot marks a present time and the black curves are the corresponding Laguerre approximations for the recent past.


## 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 n×N floating point operations (flops), where n is the number of Laguerre polynomial used and N the number of time steps that the analysis functions take to decay to a negligible value. Further, for each Gauss point, 2N 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 n×n flops, which is advantageous because n≪N. In this section, it is shown that τ¯ verifies a differential equation driven by the history v¯ 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(4.23)∂τi∂t=∫0+∞e-θLi(θ)∂v∂t(t-twθ)dθ=-1tw∫0+∞e-θLi(θ)∂v∂θ(t-twθ)dθ.
Multiplying by tw and integrating by parts yields(4.24)tw∂τi∂t=-[e-θLi(θ)v(t-twθ)]+∫0+∞[-e-θLi(θ)+e-θ∂∂θLi(θ)]v(t-twθ)dθ.
A property of Laguerre polynomials is [24](4.25)∂∂θLi(θ)=-Li-1(1)(θ)=-∑j=0i-1Lj(θ),
where ℒi(1)(θ) is a generalized Laguerre polynomial. Hence we can write(4.26)tw∂τi∂t=Li(0)v(t)-τi-∫0+∞e-θ∑j=1i-1Lj(θ)v(t-twθ)dθ=v(t)-τi-∑j=1i-1τj=v(t)-∑j=1iτj.
which is of the form(4.27)∂τ¯∂t(t)=μ¯¯⋅τ¯(t)+n¯v(t)
with(4.28)μij={-1twj≤i0j>i,ni=1tw.
Equation (4.27) shows that at any time t, 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 signalsL¯¯·τ¯ 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” (Section5.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(4.29)vj=v(t0+jdt).
We seek the values of the Laguerre coefficients at the same intervals(4.30)τ¯j=τ¯(t0+jdt).
The vector τ¯j (the list of the coefficients for all Laguerre polynomials, taken at step j) must not be confused with scalar τi (the coefficient for the Laguerre polynomial of degree i). We choose t0 such that t0+jdt=0, and we approximate v by a function that is linear over the interval [0,dt]. Equation (4.27) becomes(4.31)∂τ¯∂t(t)=μ¯¯⋅τ¯(t)+α¯+β¯t
with(4.32)α¯=n¯v(0),β¯=n¯v(dt)-v(0)dt.
This new differential equation can be solved exactly: we seek a solution of the form(4.33)τ¯(t)=exp⁡(μ¯¯t)⋅a¯+b¯t+c¯
over the interval. Here exp(μ¯¯t) stands for a matrix exponential. Replacing this expression into (4.31), noting that(4.34)∂∂texp⁡(μ¯¯t)=μ¯¯⋅exp⁡(μ¯¯t),exp⁡(0¯¯)=I¯¯,
and identifying the constant and linear terms and enforcing the initial value leads to(4.35)b¯=-μ¯¯-1⋅β¯,c¯=-μ¯¯-2⋅β¯-μ¯¯-1⋅α¯,a¯=τ¯(0)+μ¯¯-2⋅β¯+μ¯¯-1⋅α¯.
Replacing these expressions in (4.33) at t=dt, a tedious but straightforward computation yields the recursive filter(4.36)τ¯j+1=M¯¯⋅τ¯j+V¯1⋅vj+V¯2⋅vj+1
with(4.37)M¯¯=exp⁡(m¯¯dt),μ¯1=μ¯¯-1⋅n¯,μ¯2=μ¯¯-2⋅n¯1dt,V¯1=M¯¯⋅(μ¯1-μ¯2)+μ¯2,V¯2=M¯¯⋅μ¯2-μ¯1-μ¯2.


## 5. Force Interpolation

### 5.1. Foreword

This section details the implementation of point 4 in Section2.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 d for a polynomial of n variables is nd, and high values of d 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 Section7, the rotatron is used to predict scaledforces 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 2n neurons where n 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 like(5.1)σ(x)=1-2e2x+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)f̂i=Mij⋅σ(Njkl⋅τkl+Vj)+Ui.Mij, Njkl, Ui and Vj are the “weights” or interpolation coefficients, that must be adjusted to fit the perceptron to interpolate some given data. τkl are Laguerre coefficients and f̂i are predicted force components. i is the index of force direction (x versus y), j the index of neuron in the hidden layer, k the index of velocity direction, and l 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 byNjkl.


### 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 (Figure3). 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.Figure 3
For circular cross sections, it is assumed that if two trajectories that can be deduced from each other by rotation or mirroring; then the corresponding forces are deduced from each other by the same operation.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 trajectoriesand 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” (Section5.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σij=Cijklɛkl, the indexes range from 1 to 3. Following Einstein's convention, indexes k and l are summed over, and the relation is valid for any combination of i and j. 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 x, y and z were introduced to note the position to which the various tensors refer, one would tend to write σijxyz=Cijklxyzɛklxyz, 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|y[j]k| is equal to y1k2+y2k2, all values of j are present in one evaluation of the square root, and the square root is evaluated for each value of k.AppendixA 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)f̂i=Vkσik
with(5.4)σik=σi(y[j]k)(5.5)=yik|y[j]k|(1+|y[j]k|αk),(5.6)|y[j]k|=y1k2+y2k2,(5.7)αk=-1-e-Uk,(5.8)yjk=Mklτjl.
In the above, indexes i and j refer to direction, index l to the Laguerre polynomial, and k to the hidden layer. Vk, Uk, and Mkl are tunable parameters. τjl 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 x and y 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 Section5.3: if τjl is multiplied by a 2 by 2 matrix sjj′ representing a rotation or a mirroring, then yjk is multiplied by the same matrix. The Euclidean norm |y[j]k| is hence unchanged. The term yik in (5.5) is multiplied by sij′, hence so are σik and f̂i: a rotation or mirroring of the Laguerre coefficient results in the same transformation of the forces.Figure4 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.Figure 4
Laguerre analysis and rotatron transform velocity histories into a hydrodynamic force. The matrixD is the discrete form of the Laguerre analysis functions, which appear in (4.14). The dots symbolize a matrix-matrix product.The nonlinear function appearing in (5.5) is a sigmoid, whose abruptness is parametrized by Uk (Figure 5). The sigmoid is shown in Figure 5 for various values of the parameter Uk.Figure 5
Log-logistic sigmoid functions.


### 5.6. Training

“Training” of a neural network refers to finding weight coefficientsVk, Mkl, and Uk such that for any training point number m, consisting of Laguerre coefficients τjlm and two force components fim, the outputs f̂im computed by the neural network are close to fim.


#### 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 Vk, Mkl, and Uk are chosen by minimizing the cost function(5.9)J(V[k],U[k],M[k],f[im],τ[jlm])=12(fim-f̂i(τ[jl]m))2+ρ12(Uk2+Vk2+Mkl2).
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 ofyjk 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

J is a function of a large number of weight coefficients, and hence it is not practical to compute the Hessian of J, 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 J, proved very slow in this case. Hence a search method is chosen, that determines the search direction from the gradient of J [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.


## 5.1. Foreword

This section details the implementation of point 4 in Section2.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 d for a polynomial of n variables is nd, and high values of d 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 Section7, the rotatron is used to predict scaledforces 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 2n neurons where n 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 like(5.1)σ(x)=1-2e2x+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)f̂i=Mij⋅σ(Njkl⋅τkl+Vj)+Ui.Mij, Njkl, Ui and Vj are the “weights” or interpolation coefficients, that must be adjusted to fit the perceptron to interpolate some given data. τkl are Laguerre coefficients and f̂i are predicted force components. i is the index of force direction (x versus y), j the index of neuron in the hidden layer, k the index of velocity direction, and l 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 byNjkl.


## 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 (Figure3). 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.Figure 3
For circular cross sections, it is assumed that if two trajectories that can be deduced from each other by rotation or mirroring; then the corresponding forces are deduced from each other by the same operation.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 trajectoriesand 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” (Section5.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σij=Cijklɛkl, the indexes range from 1 to 3. Following Einstein's convention, indexes k and l are summed over, and the relation is valid for any combination of i and j. 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 x, y and z were introduced to note the position to which the various tensors refer, one would tend to write σijxyz=Cijklxyzɛklxyz, 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|y[j]k| is equal to y1k2+y2k2, all values of j are present in one evaluation of the square root, and the square root is evaluated for each value of k.AppendixA 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)f̂i=Vkσik
with(5.4)σik=σi(y[j]k)(5.5)=yik|y[j]k|(1+|y[j]k|αk),(5.6)|y[j]k|=y1k2+y2k2,(5.7)αk=-1-e-Uk,(5.8)yjk=Mklτjl.
In the above, indexes i and j refer to direction, index l to the Laguerre polynomial, and k to the hidden layer. Vk, Uk, and Mkl are tunable parameters. τjl 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 x and y 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 Section5.3: if τjl is multiplied by a 2 by 2 matrix sjj′ representing a rotation or a mirroring, then yjk is multiplied by the same matrix. The Euclidean norm |y[j]k| is hence unchanged. The term yik in (5.5) is multiplied by sij′, hence so are σik and f̂i: a rotation or mirroring of the Laguerre coefficient results in the same transformation of the forces.Figure4 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.Figure 4
Laguerre analysis and rotatron transform velocity histories into a hydrodynamic force. The matrixD is the discrete form of the Laguerre analysis functions, which appear in (4.14). The dots symbolize a matrix-matrix product.The nonlinear function appearing in (5.5) is a sigmoid, whose abruptness is parametrized by Uk (Figure 5). The sigmoid is shown in Figure 5 for various values of the parameter Uk.Figure 5
Log-logistic sigmoid functions.


## 5.6. Training

“Training” of a neural network refers to finding weight coefficientsVk, Mkl, and Uk such that for any training point number m, consisting of Laguerre coefficients τjlm and two force components fim, the outputs f̂im computed by the neural network are close to fim.


### 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 Vk, Mkl, and Uk are chosen by minimizing the cost function(5.9)J(V[k],U[k],M[k],f[im],τ[jlm])=12(fim-f̂i(τ[jl]m))2+ρ12(Uk2+Vk2+Mkl2).
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 ofyjk 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

J is a function of a large number of weight coefficients, and hence it is not practical to compute the Hessian of J, 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 J, proved very slow in this case. Hence a search method is chosen, that determines the search direction from the gradient of J [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.


## 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 Vk, Mkl, and Uk are chosen by minimizing the cost function(5.9)J(V[k],U[k],M[k],f[im],τ[jlm])=12(fim-f̂i(τ[jl]m))2+ρ12(Uk2+Vk2+Mkl2).
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 ofyjk 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

J is a function of a large number of weight coefficients, and hence it is not practical to compute the Hessian of J, 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 J, proved very slow in this case. Hence a search method is chosen, that determines the search direction from the gradient of J [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.Thex and y components of a trajectory are described by a pair of functions:(6.1)f¯≡(f¯x,f¯y).
We can define a scalar product between trajectories, that captures any recent differences:(6.2)f¯T⋄g¯≡f¯xT∘W¯¯∘g¯x+f¯yT∘W¯¯∘g¯y=∫-∞0et/twtw(fx(t)gx(t)+fy(t)gy(t))dt.
By replacing f¯x, f¯y, g¯x, and g¯y by their expression in terms of Laguerre polynomials and their respective Laguerre coefficients τ¯fx, τ¯fy, τ¯gx, and τ¯gy, one finds that(6.3)f¯T⋄g¯=τ¯fxT⋅τ¯gx+τ¯fyT⋅τ¯gy=τ¯fT⋅τ¯g
with(6.4)τ¯f≡[τ¯fxτ¯fy],τ¯g≡[τ¯gxτ¯gy].
The distance is defined from the scalar product in the usual manner:(6.5)|f¯-g¯|≡(f¯-g¯)T⋄(f¯-g¯),(6.6)=|τ¯f-τ¯g|.
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.Figure 6
A set of neighboring trajectories according to (6.6). Typical distance between trajectories: 103  [ilr Re].


### 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:(6.7)d(f¯,g¯)≡min⁡(min⁡R∈R|f¯-R(g¯)|,min⁡S∈S|f¯-S(g¯)|),
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 d is presented here (The vector-space of trajectories, divided by the group of rotations and mirrorings, is not a vector space).Becausef¯x is related to τ¯fx by the same linear relation that relates f¯y to τ¯fy, linear combinations of f¯x and f¯y (including rotation and mirroring) are related to the same linear combinations of τ¯fx and τ¯fy. By expressing the distance |f¯-R(g¯)| as a function of the angle α of the rotation R, and then differentiating it with respect to α, it can be shown that the value of α that minimizes |f¯-R(g¯)| is(6.8)α=arctan⁡(τ¯fx⋅τ¯gy-τ¯fy⋅τ¯gx,τ¯fy⋅τ¯gy+τ¯fx⋅τ¯gx),
where arctan(y,x)∈]-π,π] is the angle of a vector [x,y]T with the x-axis. Similarly, it can be shown that the mirroring that minimizes |f¯-S(g¯)| is the composition of a rotation of angle(6.9)β=arctan(τ¯fx⋅τ¯gy+τ¯fy⋅τ¯gx,τ¯fy⋅τ¯gy-τ¯fx⋅τ¯gx)
by a swap of the sign of the x-coordinates. Equations (6.8) and (6.9) allow to compute (6.7).Figure7 shows a trajectory and the trajectories within a small database that have the smallest distance to it, measured using d.Figure 7
A trajectory and its neighbors in terms of rotatron distance. Smooth curve: Laguerre approximation of trajectory, stippled arrow: true force, smooth arrow: predicted force. Black is for the trajectory used to enter the model to find the force. Green and red are used for the three closest points in the database, respectively, before and after rotation or mirroring.


### 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 timetw 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 (Section6.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 f¯, then the dimension of the cloud of experimental τ¯ values must be the same as the dimension of the cloud of [τ¯,f¯] points (Section 6.3.3).


#### 6.3.2. Minkowski-Bouligand Dimension

Considering a relatively uniform cloud of points, the numberm of points in a ball is proportional to the radius r of the sphere to the power of p, where p 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 m∝r20.Conversely, one can define the Minkowski-Bouligand dimension (often referred to as the fractal dimension [33]) p of a set (in particular, of a “database” of Laguerre coefficients) by counting the number m(r) of pairs of points in a set, which have a distance smaller than r:(6.10)p(r)≡∂log⁡m(r)∂log⁡r.
In practice, one should either smooth m(r) or compute the derivative by finite differences over a large enough interval. Note that the fractal dimension p is a function of the scale r.


#### 6.3.3. Functional Relation

Imagine that we have a series of data-points(x,y,z), and we are investigating whether z can be predicted using x and y. Let us imagine that the fractal dimension of the set of (x,y) pairs is 2 (the set of (x,y) fills the plane). If the fractal dimension of the set of (x,y,z) is equal to 2, then the set of (x,y,z) is within a surface, and z can be predicted using x and y. If the fractal dimension of the set of (x,y,z) is equal to 3, the data forms a cloud, and x and y are not sufficient to predict z, therefore other hidden variables must be at play. These concepts are now applied to the study of the experimental database.Figure8 shows the cumulative distribution of the distances between trajectories (“τ-distance”, black curve) computed using (6.7). p is seen to depend on the scale: from afar (r>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 (r=3×103 [ilr Re]), one can discern a cloud of dimension 4.76. At r=1.5×103 [ilr-Re], the slope decreases to about p=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 (r<1×103 [ilr Re]), the dimension increases again, possibly due to noise in the data. or weaknesses in the Laguerre approximation.Figure 8
Statistics of the distance between vectors of Laguerre coefficients (“τ distance”, black) and vectors of Laguerre coefficients and forces (“τ-f distance”, red). the four red curves, from left to right correspond to an added random noise to the force with standard deviation 0, 1×107, 2×107, and 10×107 [N/m].The red curves in Figure8 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 (“τf-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 about107 [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 f¯ appear.


## 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.Thex and y components of a trajectory are described by a pair of functions:(6.1)f¯≡(f¯x,f¯y).
We can define a scalar product between trajectories, that captures any recent differences:(6.2)f¯T⋄g¯≡f¯xT∘W¯¯∘g¯x+f¯yT∘W¯¯∘g¯y=∫-∞0et/twtw(fx(t)gx(t)+fy(t)gy(t))dt.
By replacing f¯x, f¯y, g¯x, and g¯y by their expression in terms of Laguerre polynomials and their respective Laguerre coefficients τ¯fx, τ¯fy, τ¯gx, and τ¯gy, one finds that(6.3)f¯T⋄g¯=τ¯fxT⋅τ¯gx+τ¯fyT⋅τ¯gy=τ¯fT⋅τ¯g
with(6.4)τ¯f≡[τ¯fxτ¯fy],τ¯g≡[τ¯gxτ¯gy].
The distance is defined from the scalar product in the usual manner:(6.5)|f¯-g¯|≡(f¯-g¯)T⋄(f¯-g¯),(6.6)=|τ¯f-τ¯g|.
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.Figure 6
A set of neighboring trajectories according to (6.6). Typical distance between trajectories: 103  [ilr Re].


## 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:(6.7)d(f¯,g¯)≡min⁡(min⁡R∈R|f¯-R(g¯)|,min⁡S∈S|f¯-S(g¯)|),
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 d is presented here (The vector-space of trajectories, divided by the group of rotations and mirrorings, is not a vector space).Becausef¯x is related to τ¯fx by the same linear relation that relates f¯y to τ¯fy, linear combinations of f¯x and f¯y (including rotation and mirroring) are related to the same linear combinations of τ¯fx and τ¯fy. By expressing the distance |f¯-R(g¯)| as a function of the angle α of the rotation R, and then differentiating it with respect to α, it can be shown that the value of α that minimizes |f¯-R(g¯)| is(6.8)α=arctan⁡(τ¯fx⋅τ¯gy-τ¯fy⋅τ¯gx,τ¯fy⋅τ¯gy+τ¯fx⋅τ¯gx),
where arctan(y,x)∈]-π,π] is the angle of a vector [x,y]T with the x-axis. Similarly, it can be shown that the mirroring that minimizes |f¯-S(g¯)| is the composition of a rotation of angle(6.9)β=arctan(τ¯fx⋅τ¯gy+τ¯fy⋅τ¯gx,τ¯fy⋅τ¯gy-τ¯fx⋅τ¯gx)
by a swap of the sign of the x-coordinates. Equations (6.8) and (6.9) allow to compute (6.7).Figure7 shows a trajectory and the trajectories within a small database that have the smallest distance to it, measured using d.Figure 7
A trajectory and its neighbors in terms of rotatron distance. Smooth curve: Laguerre approximation of trajectory, stippled arrow: true force, smooth arrow: predicted force. Black is for the trajectory used to enter the model to find the force. Green and red are used for the three closest points in the database, respectively, before and after rotation or mirroring.


## 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 timetw 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 (Section6.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 f¯, then the dimension of the cloud of experimental τ¯ values must be the same as the dimension of the cloud of [τ¯,f¯] points (Section 6.3.3).


### 6.3.2. Minkowski-Bouligand Dimension

Considering a relatively uniform cloud of points, the numberm of points in a ball is proportional to the radius r of the sphere to the power of p, where p 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 m∝r20.Conversely, one can define the Minkowski-Bouligand dimension (often referred to as the fractal dimension [33]) p of a set (in particular, of a “database” of Laguerre coefficients) by counting the number m(r) of pairs of points in a set, which have a distance smaller than r:(6.10)p(r)≡∂log⁡m(r)∂log⁡r.
In practice, one should either smooth m(r) or compute the derivative by finite differences over a large enough interval. Note that the fractal dimension p is a function of the scale r.


### 6.3.3. Functional Relation

Imagine that we have a series of data-points(x,y,z), and we are investigating whether z can be predicted using x and y. Let us imagine that the fractal dimension of the set of (x,y) pairs is 2 (the set of (x,y) fills the plane). If the fractal dimension of the set of (x,y,z) is equal to 2, then the set of (x,y,z) is within a surface, and z can be predicted using x and y. If the fractal dimension of the set of (x,y,z) is equal to 3, the data forms a cloud, and x and y are not sufficient to predict z, therefore other hidden variables must be at play. These concepts are now applied to the study of the experimental database.Figure8 shows the cumulative distribution of the distances between trajectories (“τ-distance”, black curve) computed using (6.7). p is seen to depend on the scale: from afar (r>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 (r=3×103 [ilr Re]), one can discern a cloud of dimension 4.76. At r=1.5×103 [ilr-Re], the slope decreases to about p=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 (r<1×103 [ilr Re]), the dimension increases again, possibly due to noise in the data. or weaknesses in the Laguerre approximation.Figure 8
Statistics of the distance between vectors of Laguerre coefficients (“τ distance”, black) and vectors of Laguerre coefficients and forces (“τ-f distance”, red). the four red curves, from left to right correspond to an added random noise to the force with standard deviation 0, 1×107, 2×107, and 10×107 [N/m].The red curves in Figure8 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 (“τf-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 about107 [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 f¯ appear.


## 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 timetw 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 (Section6.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 f¯, then the dimension of the cloud of experimental τ¯ values must be the same as the dimension of the cloud of [τ¯,f¯] points (Section 6.3.3).


## 6.3.2. Minkowski-Bouligand Dimension

Considering a relatively uniform cloud of points, the numberm of points in a ball is proportional to the radius r of the sphere to the power of p, where p 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 m∝r20.Conversely, one can define the Minkowski-Bouligand dimension (often referred to as the fractal dimension [33]) p of a set (in particular, of a “database” of Laguerre coefficients) by counting the number m(r) of pairs of points in a set, which have a distance smaller than r:(6.10)p(r)≡∂log⁡m(r)∂log⁡r.
In practice, one should either smooth m(r) or compute the derivative by finite differences over a large enough interval. Note that the fractal dimension p is a function of the scale r.


## 6.3.3. Functional Relation

Imagine that we have a series of data-points(x,y,z), and we are investigating whether z can be predicted using x and y. Let us imagine that the fractal dimension of the set of (x,y) pairs is 2 (the set of (x,y) fills the plane). If the fractal dimension of the set of (x,y,z) is equal to 2, then the set of (x,y,z) is within a surface, and z can be predicted using x and y. If the fractal dimension of the set of (x,y,z) is equal to 3, the data forms a cloud, and x and y are not sufficient to predict z, therefore other hidden variables must be at play. These concepts are now applied to the study of the experimental database.Figure8 shows the cumulative distribution of the distances between trajectories (“τ-distance”, black curve) computed using (6.7). p is seen to depend on the scale: from afar (r>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 (r=3×103 [ilr Re]), one can discern a cloud of dimension 4.76. At r=1.5×103 [ilr-Re], the slope decreases to about p=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 (r<1×103 [ilr Re]), the dimension increases again, possibly due to noise in the data. or weaknesses in the Laguerre approximation.Figure 8
Statistics of the distance between vectors of Laguerre coefficients (“τ distance”, black) and vectors of Laguerre coefficients and forces (“τ-f distance”, red). the four red curves, from left to right correspond to an added random noise to the force with standard deviation 0, 1×107, 2×107, and 10×107 [N/m].The red curves in Figure8 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 (“τf-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 about107 [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 f¯ 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, n 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 (Section7.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 (Section7.5).(4)
Finite element interpolation is introduced (space discretisation), using Gauss quadrature (Section7.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(7.1)rdi(x[bj],ẋ[bj],ẍ[bj],t)=λNm-1-1f̂d(τ[pb]i)+Edi,
where Newton's “dot” notation for a time derivative stands for a derivation with respect to unscaled time t, as opposed to scaled time t*, and with [23, 34](7.2)Edi=CLρν(ẇdi-ẋdi)+CQ12ρDi|ẇdi-ẋdi|(ẇdi-ẋdi)+CMπ4ρDi2(ẅdi-ẍdi)+π4ρDi2ẅdi.
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.IfCL, CQ, or CM are set to values different from zero, then it is necessary to subtract the corresponding values from the forces fdi used to train the rotatron. Experience shows that using CM=1, CQ=1, and CL=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)(7.3)∂τlbi∂t*=μlpτpbi+nlλms-1(ẋbi-ẇbi)
so that(7.4)λs-1τ̇lbi=μlpτpbi+nlλms-1(ẋbi-ẇbi).
The indexes d and b span pairs of directions, orthogonal to the cylinder. Indexes i and j stand for positions along the cylinder and span a continuous set of values (coordinates along the cylinder). Indexes l and p refer to the Laguerre coefficients of various degrees. Forces f̂di=f̂d(τ[pb]i) at location i only depend on the Laguerre coefficients τpbi for the same location. At that location, the force component in direction d depends on the Laguerre coefficients of all degrees b for both directions p. ρ is the fluid density and ẅdi is the acceleration of the undisturbed fluid. π/4ρDi2ẅdi stands for the Froude-Krylov forces. Diffraction forces are present in the laboratory tests and hence accounted for by f̂d.


### 7.3. Incremental Form

The incremental form of (7.1) and (7.4) is(7.5)rdi+kdibjdxbj+cdibjdẋbj+mdibjdẍbj=λNm-1-1f̂di+hdipbjdτpbj+Edi,(7.6)λs-1(τ̇lbi+dτ̇lbi)=μlp(τpbi+dτpbi)+nlλms-1(ẋbi+dẋbi-ẇdi)
with(7.7)kdibj=∂rdi∂xbj,cdibj=∂rdi∂ẋbj+CQρDiδij[|ẇ[p]i-ẋ[p]i|δbd+(ẇdi-ẋdi)(ẇbi-ẋbi)|ẇ[p]i-ẋ[p]i|-1]+CLρνδijδbd,mdibj=∂rdi∂ẍbj+CMπ4ρDiδijδbd,hdipbj=λNm-1-1∂f̂di∂τpbδij.
The expression for ∂f̂di/∂τpb 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 τlbi, but with the coefficient of τ̈lbi being zero. Alternatively, we can introduce the antiderivative Tlbi of τlbi, and treat (7.6) as being of the second order in Tlbi, but with the coefficient of Tlbi being zero. The latter option was chosen, based on the weak justification that this treats τlbi and ẋbj both as first derivatives, which seems natural considering (4.14).Applying Newmark-β to (7.5) and (7.6) in this way yields(7.8)∀d,i,[kdibj+γβdtcdibj+1βdt2mdibj]dxbj-γβdthdipbjdTpbj=λNm-1-1f̂di+Edi-rdi+cdibjbbjx+mdibjabjx-hdipbjbpbjτ,-γβdtnlλms-1dxbi+[1βdt2λs-1δlp-γβdtμlp]dTpbi=nlλms-1(ẋbi-ẇbi)+μlpτpbi-λs-1τ̇lbi-nlλms-1bbix-μlpbpbiτ+albiτ
with(7.9)abjx=1βdtẋbj+12βẍbj,bbjx=γβẋbj+(γ2β-1)dtẍbj,apbjτ=1βdtτpbj+12βτ̇pbj,bpbjτ=γβτpbj+(γ2β-1)dtτ̇pbj.
For refinement iterations, abjx, bbjx, apbjτ, and bpbjτ are set to zero. Typically, γ=1/2, β=1/4. The step dt refers to unscaled time.As usual in the Newmark-β method, the increments for the time derivatives are found from the increment as(7.10)dẋbj=γβdtdxbj-bbjx,(7.11)dẍbj=1βdt2dxbj-abjx,(7.12)dτpbj=γβdtdTpbj-bpbjτ,(7.13)dτ̇pbj=1βdt2dTpbj-apbjτ.


### 7.5. Condensation

The time discrete equations can be rewritten in a compact form:(7.14)sdibj1dxbj-sdipbj2dTpbj=sdi3,sl4dxbi+slp5dTpbi=slbi6
with(7.15)sdibj1=kdibj+γβdtcdibj+1βdt2mdibj,(7.16)sdipbj2=γβdthdipbj,(7.17)sdi3=λNm-1-1f̂di+Edi-rdi+cdibjbbjx+mdibjabjx-hdipbjbpbjτ,(7.18)sl4=-γβdtnlλms-1,(7.19)slp5=1βdt2λs-1δlp-γβdtμlp,(7.20)slbi6=nlλms-1(ẋbi-ẇbi)+μlpτpbi-λs-1τ̇lbi-nlλms-1bbix-μlpbpbiτ+λs-1albiτ.
One can then condense dTpbi out of the above system of equations:(7.21)dTpbj=(s5)pl-1(slbj6-sl4dxbj),(7.22)[sdibj1+sdipbj2(s5)pl-1sl4]dxbj=sdi3+sdipbj2(s5)pl-1slbj6.
Equation (7.22) is forced into “Newmark form” as(7.23)[kdibj*+kdibj+γβdtcdibj+1βdt2mdibj]dxbj=f̂di*-rdi+cdibjbbjx+mdibjabjx
with(7.24)kdibj*=sdipbj2(s5)pl-1sl4,f̂di*=λNm-1-1f̂di+Edi-hdipbjbpbjτ+sdipbj2(s5)pl-1slbj6kdibj* and f̂di* both depend on dt, β, and γ: the symbol kdibj* 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)Knm*=Ndinkdibj*Nbjm,F̂n*=Ndinf̂di*.Knm* and F̂n* are typically computed by Gauss quadrature. Note that no space derivative is present in kdibj*, 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)Knm*=Ndinkdib*Nbim
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 incrementdxbj, and then uses (7.21) to obtain the Laguerre coefficient increment dTpbj. From there, (7.12) and (7.13) are used to compute dτpbj and dτ̇pbj. The values of Tpbj, τpbj, and τ̇pbj are updated from previously stored values. τpbj is then used to evaluate f̂di and its derivative with respect to τpbj. These are scaled back, and Froude-Krylov forces are added, leading to kdibj* and f̂di*.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 dt 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 Section4.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.


## 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, n 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 (Section7.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 (Section7.5).(4)
Finite element interpolation is introduced (space discretisation), using Gauss quadrature (Section7.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(7.1)rdi(x[bj],ẋ[bj],ẍ[bj],t)=λNm-1-1f̂d(τ[pb]i)+Edi,
where Newton's “dot” notation for a time derivative stands for a derivation with respect to unscaled time t, as opposed to scaled time t*, and with [23, 34](7.2)Edi=CLρν(ẇdi-ẋdi)+CQ12ρDi|ẇdi-ẋdi|(ẇdi-ẋdi)+CMπ4ρDi2(ẅdi-ẍdi)+π4ρDi2ẅdi.
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.IfCL, CQ, or CM are set to values different from zero, then it is necessary to subtract the corresponding values from the forces fdi used to train the rotatron. Experience shows that using CM=1, CQ=1, and CL=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)(7.3)∂τlbi∂t*=μlpτpbi+nlλms-1(ẋbi-ẇbi)
so that(7.4)λs-1τ̇lbi=μlpτpbi+nlλms-1(ẋbi-ẇbi).
The indexes d and b span pairs of directions, orthogonal to the cylinder. Indexes i and j stand for positions along the cylinder and span a continuous set of values (coordinates along the cylinder). Indexes l and p refer to the Laguerre coefficients of various degrees. Forces f̂di=f̂d(τ[pb]i) at location i only depend on the Laguerre coefficients τpbi for the same location. At that location, the force component in direction d depends on the Laguerre coefficients of all degrees b for both directions p. ρ is the fluid density and ẅdi is the acceleration of the undisturbed fluid. π/4ρDi2ẅdi stands for the Froude-Krylov forces. Diffraction forces are present in the laboratory tests and hence accounted for by f̂d.


## 7.3. Incremental Form

The incremental form of (7.1) and (7.4) is(7.5)rdi+kdibjdxbj+cdibjdẋbj+mdibjdẍbj=λNm-1-1f̂di+hdipbjdτpbj+Edi,(7.6)λs-1(τ̇lbi+dτ̇lbi)=μlp(τpbi+dτpbi)+nlλms-1(ẋbi+dẋbi-ẇdi)
with(7.7)kdibj=∂rdi∂xbj,cdibj=∂rdi∂ẋbj+CQρDiδij[|ẇ[p]i-ẋ[p]i|δbd+(ẇdi-ẋdi)(ẇbi-ẋbi)|ẇ[p]i-ẋ[p]i|-1]+CLρνδijδbd,mdibj=∂rdi∂ẍbj+CMπ4ρDiδijδbd,hdipbj=λNm-1-1∂f̂di∂τpbδij.
The expression for ∂f̂di/∂τpb 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 τlbi, but with the coefficient of τ̈lbi being zero. Alternatively, we can introduce the antiderivative Tlbi of τlbi, and treat (7.6) as being of the second order in Tlbi, but with the coefficient of Tlbi being zero. The latter option was chosen, based on the weak justification that this treats τlbi and ẋbj both as first derivatives, which seems natural considering (4.14).Applying Newmark-β to (7.5) and (7.6) in this way yields(7.8)∀d,i,[kdibj+γβdtcdibj+1βdt2mdibj]dxbj-γβdthdipbjdTpbj=λNm-1-1f̂di+Edi-rdi+cdibjbbjx+mdibjabjx-hdipbjbpbjτ,-γβdtnlλms-1dxbi+[1βdt2λs-1δlp-γβdtμlp]dTpbi=nlλms-1(ẋbi-ẇbi)+μlpτpbi-λs-1τ̇lbi-nlλms-1bbix-μlpbpbiτ+albiτ
with(7.9)abjx=1βdtẋbj+12βẍbj,bbjx=γβẋbj+(γ2β-1)dtẍbj,apbjτ=1βdtτpbj+12βτ̇pbj,bpbjτ=γβτpbj+(γ2β-1)dtτ̇pbj.
For refinement iterations, abjx, bbjx, apbjτ, and bpbjτ are set to zero. Typically, γ=1/2, β=1/4. The step dt refers to unscaled time.As usual in the Newmark-β method, the increments for the time derivatives are found from the increment as(7.10)dẋbj=γβdtdxbj-bbjx,(7.11)dẍbj=1βdt2dxbj-abjx,(7.12)dτpbj=γβdtdTpbj-bpbjτ,(7.13)dτ̇pbj=1βdt2dTpbj-apbjτ.


## 7.5. Condensation

The time discrete equations can be rewritten in a compact form:(7.14)sdibj1dxbj-sdipbj2dTpbj=sdi3,sl4dxbi+slp5dTpbi=slbi6
with(7.15)sdibj1=kdibj+γβdtcdibj+1βdt2mdibj,(7.16)sdipbj2=γβdthdipbj,(7.17)sdi3=λNm-1-1f̂di+Edi-rdi+cdibjbbjx+mdibjabjx-hdipbjbpbjτ,(7.18)sl4=-γβdtnlλms-1,(7.19)slp5=1βdt2λs-1δlp-γβdtμlp,(7.20)slbi6=nlλms-1(ẋbi-ẇbi)+μlpτpbi-λs-1τ̇lbi-nlλms-1bbix-μlpbpbiτ+λs-1albiτ.
One can then condense dTpbi out of the above system of equations:(7.21)dTpbj=(s5)pl-1(slbj6-sl4dxbj),(7.22)[sdibj1+sdipbj2(s5)pl-1sl4]dxbj=sdi3+sdipbj2(s5)pl-1slbj6.
Equation (7.22) is forced into “Newmark form” as(7.23)[kdibj*+kdibj+γβdtcdibj+1βdt2mdibj]dxbj=f̂di*-rdi+cdibjbbjx+mdibjabjx
with(7.24)kdibj*=sdipbj2(s5)pl-1sl4,f̂di*=λNm-1-1f̂di+Edi-hdipbjbpbjτ+sdipbj2(s5)pl-1slbj6kdibj* and f̂di* both depend on dt, β, and γ: the symbol kdibj* 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)Knm*=Ndinkdibj*Nbjm,F̂n*=Ndinf̂di*.Knm* and F̂n* are typically computed by Gauss quadrature. Note that no space derivative is present in kdibj*, 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)Knm*=Ndinkdib*Nbim
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 incrementdxbj, and then uses (7.21) to obtain the Laguerre coefficient increment dTpbj. From there, (7.12) and (7.13) are used to compute dτpbj and dτ̇pbj. The values of Tpbj, τpbj, and τ̇pbj are updated from previously stored values. τpbj is then used to evaluate f̂di and its derivative with respect to τpbj. These are scaled back, and Froude-Krylov forces are added, leading to kdibj* and f̂di*.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 dt 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 Section4.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 n=30 Laguerre coefficients and the two components of the corresponding force (Figure 11).Table 1
Reynolds number in NDP tests.
Test nameReynolds numberCurrentTN203013500uniformTN23400–16200shearTN23700–24300shearThe rotatron was trained usingn=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.Figure9 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 9
Quality of prediction for 5 nodes in test TN 2030. Velocity (left), training force (middle), and predicted force (right). All velocities and all forces presented at the same scales.Figure10 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
Illustration of the approximation process. See Section8.1.Figure 11
Force vector and Laguerre approximation of velocity, for a fraction (1%) of the data used to train the rotatron. The blue cross marks the origin (zero velocity relative to water).Figure10 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 Section8.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.Table 2
Characteristics of the NDP reduced scale riser model.
QuantityValuelength38mouter diameter0.027mEI37.2Nm2EA5.09·105Nmass0.933kg/mtension3000NIn Figures12 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.Figure 12
NDP test results, test 2370 (Re∈[0,24300]). The color coding describes the displacement.Figure 13
Analysis results, test 2370 (Re∈[0,24300]).Figure 14
Analysis results, test 2030 (Re=13500).Figures12 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 Figure14, with in-line motions dominating.


## 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 n=30 Laguerre coefficients and the two components of the corresponding force (Figure 11).Table 1
Reynolds number in NDP tests.
Test nameReynolds numberCurrentTN203013500uniformTN23400–16200shearTN23700–24300shearThe rotatron was trained usingn=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.Figure9 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 9
Quality of prediction for 5 nodes in test TN 2030. Velocity (left), training force (middle), and predicted force (right). All velocities and all forces presented at the same scales.Figure10 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
Illustration of the approximation process. See Section8.1.Figure 11
Force vector and Laguerre approximation of velocity, for a fraction (1%) of the data used to train the rotatron. The blue cross marks the origin (zero velocity relative to water).Figure10 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 Section8.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.Table 2
Characteristics of the NDP reduced scale riser model.
QuantityValuelength38mouter diameter0.027mEI37.2Nm2EA5.09·105Nmass0.933kg/mtension3000NIn Figures12 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.Figure 12
NDP test results, test 2370 (Re∈[0,24300]). The color coding describes the displacement.Figure 13
Analysis results, test 2370 (Re∈[0,24300]).Figure 14
Analysis results, test 2030 (Re=13500).Figures12 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 Figure14, with in-line motions dominating.


## 9. Discussion

### 9.1. Force Prediction and Reynolds Interpolation

Figure9 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.Figure 15
Quality of prediction on a control data set. See for example, color code of Figure9.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 Figures12 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.Figure14 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 Figure9, 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 [38–40].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 oftw

One issue that was explored was the adequacy oftw: the rotatron was trained with tw=5·10-5. 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 tw. This was done twice, once with the full training set and increasing the number of Laguerre polynomials to n=60, and once with only the experimental data from lower currents and an unchanged number of polynomials (n=30). Neither rotatrons allowed to perform a successful simulation for the lower current velocities.


### 9.4. Tension

In Figure12, 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.


## 9.1. Force Prediction and Reynolds Interpolation

Figure9 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.Figure 15
Quality of prediction on a control data set. See for example, color code of Figure9.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 Figures12 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.Figure14 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 Figure9, 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 [38–40].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 oftw

One issue that was explored was the adequacy oftw: the rotatron was trained with tw=5·10-5. 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 tw. This was done twice, once with the full training set and increasing the number of Laguerre polynomials to n=60, and once with only the experimental data from lower currents and an unchanged number of polynomials (n=30). Neither rotatrons allowed to perform a successful simulation for the lower current velocities.


## 9.4. Tension

In Figure12, 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 Section6.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?



---

*Source: 414702-2011-07-10.xml*