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.
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.
2. Model Outline
2.1. Postulate
The present work hinges on the following postulate. The force exerted by the surrounding fluid on a section of the slender structure is completely determined by the recent histories at that section of the velocities of the structure and of the undisturbed fluid. Several points in this sentence are worthy of discussion.
The “force” includes the components usually distributed into added mass, excitation forces, drag, lift, and so forth.
That the force “at a section of the slender structure” is determined by the history “at that section” implies a “strip theory” in which it is excluded that motions of the structure at a point A cause disturbances in the fluid that affect the force at point B away from A. In other words, it is assumed that there is no significant transmission of information in the axial direction within the water (as opposed to within the slender structure). This would be proved wrong if it turned out that unstable phenomena, like boundary layer shedding, although transmitting little energy along the structure, transmit information that steers how local hydrodynamic energy is channeled at a given point along the structure.
That the force should be “completely determined” implies that the behavior of the structure is deterministic. This does not contradict the observation of hysteretic response of short cylinders mounted on elastic support. Uniqueness of forces for a given position does not imply uniqueness of static equilibrium. Neither does “completely determined” contradict the observation of irregular and unpredictable responses to VIV: nonlinear dynamic systems can have a chaotic behavior. Still, complete determinism is provably wrong, since a short vertical cylinder dragged at uniform speed through water will experience oscillating lift forces. At any given moment, there is nothing in the history of (constant) velocity that allows to predict whether the lift is left or right. So the present work is based on the bet that ignoring such “bifurcations” still leaves us with a useful model.
“Recent” can be defined as anything between the present time and a few times , where the value of still is an object of debate. is likely to be case dependent. The current will transport (convect) away vortices so that they quickly loose significance. The time should then be of the order of where is the cross section diameter and the current velocity. In contrast, if the cylinder is oscillating in still water, it will be traveling in its own wake, and should be related to the rate of diffusion and/or viscous dissipation of vortices, which is likely to result in much higher values of . Tests on periodic forced motion of short cylinders sometimes show a slow drift of the forces (over as many as ten periods). In contrast, force decay tests for a cylinder stopped after oscillations at zero mean velocity, point towards a fraction of a period. In the present work, the idea is to choose an upper bound for , after adequate scaling (cf. Section 3.2).
The “velocity histories” are what count. Accelerations would not do because for example, zero acceleration can correspond to different speeds and hence different forces. On the other hand, the force on a cylinder will not be affected by a uniform translation of its whole trajectory, so a history of positions contains irrelevant information.
In the remainder of this text, the word “trajectory” will be given a very specific meaning. The trajectory is defined as the recent history of the velocity vector of the cylinder relative to the undisturbed surrounding fluid.
2.2. Restrictions
In the present phase of research, the following restrictions are introduced, in order to achieve some simplification of the task. The outer cross section of the slender structure is assumed perfectly circular and smooth. The surrounding fluid is assumed to be infinite, excluding the presence of sea floor, free surface, or neighboring risers. Only fluid flows perpendicular to the cylinder at any point are considered.
2.3. Input and Output
As stated earlier, the VIV model being developed here replaces the Morison model for wave-induced loads. The VIV model is called at each step and iteration, and at each Gauss point or node of each element.
The model is to receive as input: (1)the diameter of the cylinder, the viscosity, and density of the surrounding fluid, (2)the instantaneous velocity of the cross section relative to the undisturbed fluid, (3)the instantaneous velocity and acceleration of the local undisturbed fluid, in a Galilean reference system.
The model uses velocity information stored from previous steps. On this basis, the model produces as output:(1)the vector of hydrodynamic forces per unit length, acting on the cylinder, (2)the matrix containing the derivative of the above with respect to instantaneous values of the cylinder’s behavior.
Gauss integration is then used to compute a consistent load vector and partial derivative matrices (damping, stiffness, and mass) for each element. Note that these element matrices are likely to vary significantly over each VIV oscillation “period”—in contrast to added mass or damping matrices, deemed to be constant over a long time in semiempirical VIV models. The connection of the force model to the finite element analysis is discussed in Section 7.
2.4. Algorithmic Steps
Only the local VIV model is described here, not the whole FE analysis. (1)The relative velocity of the cylinder relative to water (thereafter: “velocity”) is computed. (2)The velocity is scaled (Reynolds scaling) to that the cylinder diameter is the unit of distance (Section 3.2). (3)The trajectory (again: the recent histories of both and components of velocity) is compressed into a small number of “Laguerre coefficients”. This compression is such that it provides detailed information over the recent past and increasingly coarse information for the more distant past (Section 4). (4)The Laguerre coefficients are used to enter an interpolation function (a feed-forward neural network with some specifically tailored properties) which returns and components of hydrodynamic force (Section 5). The fitting of the interpolation function is discussed in Section 8.1. (5)The force is scaled back to the relevant diameter (Section 3.2). (6)The Froude-Krylov forces, which depend on the acceleration of the undisturbed flow, are added (Section 3.1)
The identification of nonlinear systems using a bank of orthogonal filters (including Laguerre filter) to generate multiple signals from a single one, and then using the multiple signals to enter a nonlinear, memory-less function, was introduced by Wiener [19] (Wiener-Laguerre filtering). In the present work, a base of Laguerre polynomials is used, in contrast to Laguerre functions introduced by Wiener. While Wiener apparently did not use neural networks as nonlinear functions (but e.g., Hermite polynomials), neural networks in Wiener models have been studied for some time [20–22]. In the present work, Laguerre filtering is presented without making use of the vocabulary of cybernetics. In particular, the -transform is not introduced here.
3. Exploiting Similitudes
3.1. Froude-Krylov Forces
This section gives the justification for point 6 of Section 2.4. If the undisturbed fluid in which the cylinder is plunged is accelerating (because of surface waves, e.g.,), then it is natural to introduce two reference systems: is a Galilean reference system, for example, fixed relative to the sea floor and is an accelerated reference system, locally following the undisturbed flow. Transforming the equations of equilibrium from (in which we carry out FEM analysis) to (for which we have experimental data, in water that is not accelerated) requires the addition of inertia forces.
The inertial forces create a uniform pressure gradient that was not present in the laboratory test. The effect of a pressure gradient on a submerged body is variously referred to as “Archimedes forces” when the pressure gradient results from the acceleration of gravity, or as “Froude-Krylov forces” when the pressure gradient is due to fluid acceleration in for example, surface waves. As familiar, the integral of the pressure over the wet surface is transformed into a volume integral [23].
It is assumed that this pressure gradient does not affect the turbulent flow, so that the pressure gradient can simply be added to the pressures resulting from turbulence. This seems reasonable enough for incompressible flows, and indeed when it comes to Archimedes forces, the submerged weight of a cylinder is routinely subtracted to laboratory measurements and the relevant correction added again in FEM analysis—even though the Archimedes forces in the laboratory do not necessarily scale with those in the analysis. To the author’s knowledge, there is no experimental indication that a horizontal and vertical cylinder, all other conditions being equal, experience different forces.
To conclude, the hydrodynamic force acting on the cylinder at a given instant is the sum of two terms: (1)a force that is a function of only the cylinder diameter and the recent history of the velocity of the cylinder relative to the undisturbed, steady water flow, (2)Froude-Krylov forces.
All computations in Sections 4 and 5 deal only with the first of the above two terms.
3.2. Scaling
This section details how points 2 and 5 of Section 2.4 are implemented. In order to reduce the amount of experimental data necessary to create the interpolation function used in point 4, one must take advantage of scale similarities. To that effect, all data used to either train or query the database is scaled. Correspondingly, all forces returned by the database are scaled back. The present model uses a scaling that is quite different from the scaling typical used by simplified harmonic models (Section 1.3): reduced amplitudes and frequencies are not used. Instead velocities histories are scaled in a manner familiar from the Reynolds number.
VIV forces are assumed to be uniquely defined by fluid density , kinematic viscosity , cylinder diameter , and the motion history. Hence, in order to create a database that is to be entered with scaled velocities, we wish all experimental data to be scaled to fixed reference values , , and . The choice of , , and is arbitrary, and in this work, all are set to the value 1.
By expressing the units of these quantities, one gets three equations on , , and , which are the scaling factors for the basic units of distance, time, and mass. Solving the system yields 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: The choice [m] hence implies that scaled displacements can be considered to have “1 diameter” as unit. Similarly, the choices [m] and [m2/s] together imply that scaled velocities are expressed as Reynolds numbers since the scaled velocity is calculated as where is the velocity. The Reynolds number is usually computed using some velocity that is characteristic of the system under study. In VIV science, the undisturbed velocity of the current is used. By contrast, in this work, instantaneous local values of the relative velocity vector are multiplied by . The scaled velocities thus obtained are a generalization of the traditional use of Reynolds number: considering an immobile cylinder in a current, the norm of its scaled relative velocity vector is equal to the traditional Reynolds number. To prevent confusion of the present usage of Reynolds number with the more particular classical one, yet emphasize the relation between both, the expression “ilr-Reynolds” (for “instant, local, relative Reynolds”) will be used in this document.
Since scaling is applied consistently to all derived quantities, all nondimensional numbers based on combinations of distance, time and mass (including Reynolds and Froude numbers) are conserved. However, any dimensional quantity with units different from those of , , and is scaled to values that depend of , , and . In particular, (3.3) shows that all accelerations, including the acceleration of gravity , are scaled with a factor proportional to . So while the scaling used here may conserve Froude's number, it does not allow to build a database of forces related to surface wave effects, because the database does not refer to a uniform value .
4. Characterization of Trajectory
4.1. Foreword
This section details how point 3 in Section 2.4 is to be implemented. The objective is, for any given point in time, to distill a “summary” of the recent history of the velocity of the cylinder relative to the surrounding fluid (trajectory). Note that the history of each component of the velocity vector is treated separately in this section and that the procedure is applied to the scaled trajectory.
The trajectory is approximated as a linear combination of some adequate family of functions, and the coefficients in this linear combination are the summary. The family of functions that is used here is the series of Laguerre polynomials (Section 4.2). It is shown in Section 4.3 that if the “Laguerre coefficients” of the linear combination are obtained by integrating the product of the trajectory by adequate “Laguerre analysis functions”, then the difference between the approximating linear combination and the real trajectory is small in the recent past and larger in the further past. This justifies the choice of Laguerre polynomials: they allow to summarize the trajectory in a way that represents recent velocities very precisely, and older velocities in a coarser manner. It is assumed that this corresponds to the information needed to obtain a good estimate of the hydrodynamic force.
Computing the integral of the product of Laguerre analysis functions and trajectory takes time. Luckily, one can show (Section 4.5) that the Laguerre coefficients are the solution of a differential equation driven by the instant value of the velocity. To obtain results that are independent of step size, this differential equation must be carefully discretised in time (Section 4.6) when summarizing experimental data.
4.2. Definitions
The Laguerre polynomial (Figure 1, top) of degree can be defined by its Rodrigues formula [24] Laguerre polynomials verify the orthonormality property 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) where the interpretation of has been discussed in Section 2.1. Functions will now be noted as vectors (in a Hilbert space), marked with overline symbols. An indexed family of functions will be noted as a matrix (symbols with double overline) and so will a linear operator (a distribution of two variables). We introduce the symmetric positive definite operator and a dot product in a suitable space of real valued functions with the canonical norm Further we introduce the base Equation (4.2) can be rewritten in matrix notation as where is the identity matrix. It is useful to introduce the weighted norm or -norm Note that since is of dimension , is of the same dimension as . So when taking as a scaled velocity, is an ilr-Reynolds number.
4.3. Analysis and Synthesis
For a history of either the or component of the velocity, we seek the vector of “Laguerre coefficients” with which to combine the columns , that minimize the weighted error defined as 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: which implies with The “Laguerre analysis functions” (Figure 1, middle) are by definition equal to The Laguerre analysis functions must not be confused with the Laguerre functions (note the factor 2 in (4.18)). Incidentally, Wiener-Laguerre models use Laguerre filters whose impulse response is Laguerre functions (not analysis function), so the present approach is slightly different from the classical Wiener-Laguerre model. A justification for the present choice will appear in Section 6.1.
4.4. Convergence
Laguerre functions, which can be defined as or in matrix notation as 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 This can be used to obtain a result on the convergence of series of Laguerre polynomials. We introduce the change of variables so that We hence have convergence in terms of the quality of approximation that we are seeking, with emphasis on the recent past. Further, on any finite (or “compact”) interval, convergence in the -norm is equivalent to convergence almost everywhere. So under some conditions of continuity on , the series of Laguerre polynomials obtained using as analysis functions converges almost everywhere towards in any finite interval.
Figure 2 illustrates how Laguerre coefficients indeed provide a “summary” of the trajectory
4.5. Differential Equation for Laguerre Coefficients
In the finite element analysis, we need to update the Laguerre coefficients at each iteration of each time step, for every Gauss point of every node of the system. The explicit calculation of (4.14) for every update is hence a CPU-time critical operation, taking in the order of floating point operations (flops), where is the number of Laguerre polynomial used and the number of time steps that the analysis functions take to decay to a negligible value. Further, for each Gauss point, velocity values need to be stored, a severe memory requirement.
In the present section and the next, it is shown how the computation of (4.14) can be carried out by a recursive operation requiring no other storage than that of the Laguerre coefficients and the last velocity values, and taking in the order of flops, which is advantageous because . In this section, it is shown that verifies a differential equation driven by the history of the velocity component. In Section 4.6, this differential equation is solved time-step by time-step in a recursive update.
Equation (4.14) can be rewritten without matrix notation and differentiated Multiplying by and integrating by parts yields A property of Laguerre polynomials is [24] where is a generalized Laguerre polynomial. Hence we can write which is of the form with Equation (4.27) shows that at any time , the rate of the Laguerre coefficients is fully defined by the Laguerre coefficients and the velocity signal.
4.6. Recursive Filter
The discrete integration of (4.27) must be done carefully, for two reasons. Firstly, it is important to obtain Laguerre coefficients that are independent of the sampling rate used (as long as the sampling rate is “adequate”). This is because the experimental data on which the VIV model is based may come from experiments which, after scaling, may have different sampling rates. Further, the numerical analysis in which the VIV model is used may use yet another time step. The choice of time step or sampling rate must not affect the way a trajectory is characterized by Laguerre coefficient.
The second reason for care in discrete integration is that we wish to be able to create synthesized signals of good quality. Synthesized signals are neither used in the numerical process of creating a force interpolation function (Section 5) or in the FEM use of the VIV model. However, visualization is essential to the process of research, both for fault diagnosis and quality control, and to communicate an understanding of the method.
This discrete integration is only used in the analysis of experimental data, to provide an input to the training of the “rotatron” (Section 5.6). In dynamic analysis, the integration of (4.27) is done by means of the Newmark- method, as detailed in Section 7.
Assume that velocity is sampled at regular intervals We seek the values of the Laguerre coefficients at the same intervals The vector (the list of the coefficients for all Laguerre polynomials, taken at step ) must not be confused with scalar (the coefficient for the Laguerre polynomial of degree ). We choose such that , and we approximate by a function that is linear over the interval . Equation (4.27) becomes with This new differential equation can be solved exactly: we seek a solution of the form over the interval. Here stands for a matrix exponential. Replacing this expression into (4.31), noting that and identifying the constant and linear terms and enforcing the initial value leads to Replacing these expressions in (4.33) at , a tedious but straightforward computation yields the recursive filter with
5. Force Interpolation
5.1. Foreword
This section details the implementation of point 4 in Section 2.4. This section presents an interpolation function which, given the Laguerre coefficients, predicts the present value of the force vector. Polynomials were considered initially, but it soon became clear that feed-forward “neural networks” provide a better class of functions to work with. The reason for this is that the number of polynomial coefficients of degree for a polynomial of variables is , and high values of must be expected to be necessary. By contrast, in a neural network, nonlinearity is introduced by “sigmoid” or “threshold” functions, and the coefficients are used to specify in which direction nonlinearity applies. Further, polynomials are infamous for their propensity to oscillate.
The “rotatron” presented here is based on the “perceptron” [26, 27], a well-studied architecture of neural network which provides a flexible tool for the interpolation of scalar-valued functions of a vector (Section 5.2). The rotatron takes advantage of certain symmetry properties of the physics at hand (Section 5.3).
In Section 7, the rotatron is used to predict scaled forces based on the Laguerre coefficients for scaled trajectories.
5.2. Perceptron
The perceptron [26, 27] is a simple feed-forward neural network, consisting of 3 layers. The input layer has neurons where is the number of Laguerre coefficients for each velocity component and the factor 2 comes from the need to analyze in-line and cross-flow speed histories together. The values of the input layer neurons are set to the Laguerre coefficients for both velocity components. The second layer has nhid neurons, whose values are an affine function of the values of the first layer, passed through a sigmoid function like 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, , and are the “weights” or interpolation coefficients, that must be adjusted to fit the perceptron to interpolate some given data. are Laguerre coefficients and are predicted force components. is the index of force direction ( versus ), the index of neuron in the hidden layer, the index of velocity direction, and the index of Laguerre coefficient.
Each output of the perceptron can be seen as a function, which is a sum of sigmoid steps in directions defined by .
5.3. Symmetries
The relation between trajectories (in the sense of the history of the velocity of the cylinder relative to the water) and forces can reasonably be assumed to exhibit several symmetries (Figure 3). Rotational symmetry: if a trajectory can be deduced from the other by a rotation around the origin, then the resulting forces are also deduced from each other by the same rotation. Mirror symmetry: if a trajectory can be deduced from the other by a mirroring around a line crossing the origin, then the resulting forces are also deduced from each other by the same mirroring.
Rotational symmetry and mirror symmetry together imply directionality: if a trajectory is within a line crossing the origin, then the resulting forces are within the same line. In particular, zero velocities must imply zero forces.
The symmetries imply that once experimental data for a trajectory has been obtained, there is no need to acquire data for rotated or mirrored trajectories. However, if one was training a perceptron to interpolate the data, the training set would need to include trajectories and their rotates and mirrors, with the correspondingly rotated and mirrored forces. This would increase memory and CPU usage during training, but also during use of the trained perceptron, because the perceptron will need a larger number of hidden layers to interpolate the training data.
Another approach is hence used in the present work: the classic perceptron is replaced by a “rotatron” (Section 5.5). It is designed so that, whatever the values of the weight coefficient, a rotation or mirroring of the input trajectory results in the same rotation or mirroring of the output force vector.
5.4. Index Notations
In the present work, index notations inspired from tensor analysis are used. However, the present setting differs from tensor analysis in at least three ways.
Firstly, we assume that we are only operating in Euclidean spaces (and not in more general Riemannian manifolds) so that orthogonal bases can be used. This makes it unnecessary to distinguish between co- and contravariant bases and coordinates. Hence, only lowered indexes appear in the present work. Incidentally, it was here assumed that the state of the model is a point in a vector space, which is not true when finite rotations are present and Riemannian geometry should be introduced instead.
Secondly, in tensor notations, each index spans the dimension of the manifold. In an expression like , the indexes range from 1 to 3. Following Einstein's convention, indexes and are summed over, and the relation is valid for any combination of and . The fact that the equation is valid at each point within a solid is implicit in the notation. In the present work, we prepare for the manipulations of arrays in a computer, involving operations that are repeated, for example, for various locations along a riser. If indexes , and were introduced to note the position to which the various tensors refer, one would tend to write , which violates Einstein's convention, because no summation (or rather: no integral) is implied over the positions.
Thirdly, we introduce nonlinear functions. These functions can combine the values of the coordinates for some indexes (these will be noted in brackets) and operate in parallel on the coordinates for other indexes (not brackets). For example, by definition of the notation is equal to , all values of are present in one evaluation of the square root, and the square root is evaluated for each value of .
Appendix A gives a detailed description of the conventions used.
5.5. Rotatron
A modified interpolation function (which will be referred to as “rotatron” in this text) is defined as with In the above, indexes and refer to direction, index to the Laguerre polynomial, and to the hidden layer. , , and are tunable parameters. are Laguerre coefficients, given as input to the “rotatron”. Note that (4.14), (5.3), and (5.8) operate linearly, identically, and independently of the terms related to the and directions, while (5.5) involves a unit vector multiplied by a nonlinear function of its norm.
The rotatron can be shown to enforce the symmetries discussed in Section 5.3: if is multiplied by a 2 by 2 matrix representing a rotation or a mirroring, then is multiplied by the same matrix. The Euclidean norm is hence unchanged. The term in (5.5) is multiplied by , hence so are and : a rotation or mirroring of the Laguerre coefficient results in the same transformation of the forces.
Figure 4 illustrates the flow of information, from right to left, from two vectors containing the histories of the velocity components, to Laguerre coefficient, that are then processed in the rotatron.
The nonlinear function appearing in (5.5) is a sigmoid, whose abruptness is parametrized by (Figure 5). The sigmoid is shown in Figure 5 for various values of the parameter .
5.6. Training
“Training” of a neural network refers to finding weight coefficients , , and such that for any training point number , consisting of Laguerre coefficients and two force components , the outputs computed by the neural network are close to .
5.6.1. Regularization
A common problem when training neural networks is “overspecialization” [28]. In this situation, the neural network predicts the training outputs with high accuracy but behaves wildly between the training points. In contrast, what is implicitly sought is a smooth response of the network to the input, even if this means an imperfect fit to the training data.
Many strategies are described in the literature to address this problem. One of them, which is adopted here, is regularization [28]: the value of the weight parameters , , and are chosen by minimizing the cost function The “regularization coefficient” is an arbitrary input to the training algorithm. High values of favor smoothness of the response of the neural network against precision in reproducing the training set.
Arguments of symmetry by permutation of the numbering of the first index of in (5.6) show that the cost function has multiple minima. Further, there are probably local minima higher than the lowest maxima.
5.6.2. Conjugate Gradient Optimization
is a function of a large number of weight coefficients, and hence it is not practical to compute the Hessian of , because the Hessian is a full matrix. It also proves to be very costly to even compute an approximation to it as done in the Levenberg-Marquardt algorithm [29, 30]. On the other hand, the Nelder-Mead “downhill simplex” algorithm [31], which uses only the values of , proved very slow in this case. Hence a search method is chosen, that determines the search direction from the gradient of [32]. This is a conjugate gradient method, in which the step length is found by deriving the gradient in the direction of the search. In this method, the positive definiteness of the (implicit) Hessian is forced by adding a scaled identity matrix to it, a technique known as “trust region”.
The conjugate gradient method proved far more efficient than the Levenberg-Marquardt and Nelder-Mead methods for the present task.
The weight coefficients are set to random values at the start of the conjugate gradient iterations.
6. Metric
6.1. Euclidean Metric and Distance
In order to describe the available data, it is useful to define a distance between trajectories. This will allow to determine to what extend the set of available data “fills” the set of all possible trajectories, or to detect zones of transition from one hydrodynamic behavior to the other. Finally, this will help detecting contradictions in the available data, arising from a variety of sources, including hidden experimental variables, measurement uncertainties, or inadequate modeling in inverse methods and not least, the natural variability of VIV forces.
The and components of a trajectory are described by a pair of functions: We can define a scalar product between trajectories, that captures any recent differences: By replacing , , , and by their expression in terms of Laguerre polynomials and their respective Laguerre coefficients , , , and , one finds that with The distance is defined from the scalar product in the usual manner: In other words, neighboring vectors of Laguerre coefficients describe trajectories that are similar in the recent past. This is illustrated by taking random samples of Laguerre coefficients around a given value obtained from data analysis and plotting the synthesized trajectories (Figure 6). Given that trajectories are compared using an integral weighted with an exponential (4.2). This important property can be seen as the justification of the choice of the Laguerre coefficients to “summarize” trajectories.
6.2. Rotatron Distance
The above does not account for rotational and mirror symmetries. We seek a distance for which the distance of a trajectory to its transforms by rotation or mirroring is zero. Another distance is hence introduced: where is the set of all rotations of the trajectories around the origin and the set of all mirroring of trajectories around a line passing by the origin. Note that no norm or scalar product associated to the distance is presented here (The vector-space of trajectories, divided by the group of rotations and mirrorings, is not a vector space).
Because is related to by the same linear relation that relates to , linear combinations of and (including rotation and mirroring) are related to the same linear combinations of and . By expressing the distance as a function of the angle of the rotation , and then differentiating it with respect to , it can be shown that the value of that minimizes is where is the angle of a vector with the -axis. Similarly, it can be shown that the mirroring that minimizes is the composition of a rotation of angle by a swap of the sign of the -coordinates. Equations (6.8) and (6.9) allow to compute (6.7).
Figure 7 shows a trajectory and the trajectories within a small database that have the smallest distance to it, measured using .
6.3. Existence of Functional Relation
6.3.1. Introduction
From experiments, one can obtain databases of velocity histories, and the corresponding hydrodynamic forces. The velocity histories can be “summarized” into Laguerre coefficients. An open question at this stage is whether there is actually a functional relationship between the Laguerre coefficients and the forces, that the rotatron could be used to approximate. For example, if the decay time was chosen too small, then the same velocity in the recent past (represented by the Laguerre coefficient) will result in different forces, depending on the velocity in the more remote past. In that case, there is no functional relation. The functional relation could also be lost if too few Laguerre coefficients are used, or if the force in the experiments is affected by incoming turbulence. In the absence of a functional relation, attempting to train the rotatron will not give useful results, in the same way as no curve drawn through a cloud of points on a piece of paper will provide a useful model.
The approach in this section hinges on two ideas. Firstly, the higher the dimension of a spaceis, the faster the volume of a ball increases with its radius. Hence, the statistic of the distances between points in a cloud of data can be used to define a dimension (Section 6.3.2). For example, a set of points on a piece of paper will be shown to follow a curve if the number of neighbors increases linearly with the radius. Secondly, if Laguerre coefficients predict forces , then the dimension of the cloud of experimental values must be the same as the dimension of the cloud of points (Section 6.3.3).
6.3.2. Minkowski-Bouligand Dimension
Considering a relatively uniform cloud of points, the number of points in a ball is proportional to the radius of the sphere to the power of , where is the dimension of the space in which the cloud is defined. For example, using Laguerre coefficients to describe both components of a trajectory, the number of points in the ball would be .
Conversely, one can define the Minkowski-Bouligand dimension (often referred to as the fractal dimension [33]) of a set (in particular, of a “database” of Laguerre coefficients) by counting the number of pairs of points in a set, which have a distance smaller than : In practice, one should either smooth or compute the derivative by finite differences over a large enough interval. Note that the fractal dimension is a function of the scale .
6.3.3. Functional Relation
Imagine that we have a series of data-points , and we are investigating whether can be predicted using and . Let us imagine that the fractal dimension of the set of pairs is 2 (the set of fills the plane). If the fractal dimension of the set of is equal to 2, then the set of is within a surface, and can be predicted using and . If the fractal dimension of the set of is equal to 3, the data forms a cloud, and and are not sufficient to predict , therefore other hidden variables must be at play. These concepts are now applied to the study of the experimental database.
Figure 8 shows the cumulative distribution of the distances between trajectories (“-distance”, black curve) computed using (6.7). is seen to depend on the scale: from afar ([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 ( [ilr Re]), one can discern a cloud of dimension 4.76. At [ilr-Re], the slope decreases to about , and it is believed that this is the dimension of the data set for a given point along the riser. At a small scale ( [ilr Re]), the dimension increases again, possibly due to noise in the data. or weaknesses in the Laguerre approximation.
The red curves in Figure 8 are computed by adding the sum of squares of the differences between force components (suitably scaled) to the squares of the distances between trajectories, and then extracting the square root (“-distance”). The four red curves, from left to right, are drawn using the original force data, to which Gaussian noise of standard deviation 0, , , and [N/m], respectively, has been added (The standard deviation of the original force is about [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 , and , [N/m], and already at [N/m] the curve is distinct from the one based on the original data.
6.3.4. Discussion
The present study suggests that there is indeed a functional relation to be seen in the data set used here. The Laguerre coefficients can be used to predict the forces, with a precision of about [N/m]. The conclusion must be treated carefully, however. The dimensional analysis provides a necessary condition, not a sufficient one: it does not exclude, for example, that there exists a neighborhood of points in which two distinct values of appear.
7. Dynamic Analysis
7.1. Foreword
Once it is possible to predict hydrodynamic forces on a cross section for a given velocity history, the next development is to include the force thus predicted in a dynamic time domain simulation. Because the VIV forces introduce severe nonlinearities, a naive connection (where the forces are just added to the right-hand side of the system) might lead to slow convergence or to divergence of the Newton-Raphson iterations used at each time step. To obtain a proper formulation, it is necessary to jointly treat the system of differential equations composed of the state equations of the structure and the differential equations (4.27) followed by the Laguerre coefficient. However, in doing so, for each displacement degree of freedom, Laguerre coefficients are added, and it is crucial for efficiency to eliminate them before solving a large linear system of equations.
To this effect, in this section, the following sequence of transformations is applied to the differential equations. (1)The differential equations are first set in incremental form (Section 7.3). (2)Time discretisation by the Newmark- method is introduced (Section 7.4). (3)The Laguerre coefficients are condensed out of the system of equations (Section 7.5). (4)Finite element interpolation is introduced (space discretisation), using Gauss quadrature (Section 7.6).
This particular sequence leads to a VIV model that is implemented at the Gauss point level and can easily be introduced in a general purpose FEM software with standard, displacement-based beam, or cable elements. Another sequence, 1, 4, 2, 3, can be used to obtain either a hybrid element, or alternatively, a mixed element which would require a specialized solver for optimal efficiency. These alternatives are more difficult to integrate into existing software working with displacement-based elements and are not discussed here.
7.2. Differential Equations
The dynamic differential equation of a 3D beam subjected to VIV loads can be formalized as where Newton's “dot” notation for a time derivative stands for a derivation with respect to unscaled time , as opposed to scaled time , and with [23, 34] The four terms in the above Morison's equation are the linear drag, the quadratic drag, the sum of diffraction and added mass forces, and the Froude-Krylov forces. The fourth term introduces the correction discussed in Section 3.1.
If , , or are set to values different from zero, then it is necessary to subtract the corresponding values from the forces used to train the rotatron. Experience shows that using , , and 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) so that The indexes and span pairs of directions, orthogonal to the cylinder. Indexes and stand for positions along the cylinder and span a continuous set of values (coordinates along the cylinder). Indexes and refer to the Laguerre coefficients of various degrees. Forces at location only depend on the Laguerre coefficients for the same location. At that location, the force component in direction depends on the Laguerre coefficients of all degrees for both directions . is the fluid density and is the acceleration of the undisturbed fluid. stands for the Froude-Krylov forces. Diffraction forces are present in the laboratory tests and hence accounted for by .
7.3. Incremental Form
The incremental form of (7.1) and (7.4) is with The expression for is presented in Appendix B.
7.4. Time Discretisation
Newmark-β is a method geared towards second order differential equations. Equation (7.6), however, is only of the first order. The reason for using Newmark-β here is that the present VIV model is to be integrated into the model of a larger structure, and the differential system for that structure is of the second order. Preparing a first order equation for a second order solver opens two options: we can treat (7.6) as being of the second order in , but with the coefficient of being zero. Alternatively, we can introduce the antiderivative of , and treat (7.6) as being of the second order in , but with the coefficient of being zero. The latter option was chosen, based on the weak justification that this treats and both as first derivatives, which seems natural considering (4.14).
Applying Newmark- to (7.5) and (7.6) in this way yields with For refinement iterations, , , , and are set to zero. Typically, , . The step refers to unscaled time.
As usual in the Newmark-β method, the increments for the time derivatives are found from the increment as
7.5. Condensation
The time discrete equations can be rewritten in a compact form: with One can then condense out of the above system of equations: Equation (7.22) is forced into “Newmark form” as with and both depend on , , and : the symbol was chosen to indicate that the matrix is handled by the Newmark- solver in the same way as a stiffness. However, this term cannot be interpreted physically as a stiffness.
7.6. Spacial Discretisation
The consistent discretisation by Galerkin finite elements of (7.23) leads to and are typically computed by Gauss quadrature. Note that no space derivative is present in , so no partial integration or Gauss quadrature with curvature shape function appears. One can hence simplify the expression of the element matrix to which means “same quadrature as for a mass matrix”.
7.7. Implementation
In nonlinear FEM code, incremental matrices and vectors are computed by Gauss quadrature. The Gauss quadrature involves shape functions, tensors that are local, continuous versions of the stiffness, damping and mass matrices, and the force imbalance vector. For example, for the drag damping of a beam element, the tensor relates a vector whose components are increments in velocities in three directions, to another vector whose components are increments in forces per unit length in three directions.
Within an iteration, the linear solver provides incremental nodal positions, velocities and accelerations for the model. These are disassembled and provided to the elements. The elements compute positions, velocities, and accelerations (and more) in a corotated reference system at Gauss points. The resulting values are handed to the VIV-Gauss point procedure.
The axial velocities are discarded. The procedure scales the provided values using (3.2) and (3.3).
Having stored the previous approximation of the scaled position, the procedure determines the position increment , and then uses (7.21) to obtain the Laguerre coefficient increment . From there, (7.12) and (7.13) are used to compute and . The values of , , and are updated from previously stored values. is then used to evaluate and its derivative with respect to . These are scaled back, and Froude-Krylov forces are added, leading to and .
The above matrix and vector are padded with zeros to indicate zero force in the axial direction and zero torque.
The condensation of a larger system of time-discretised equations introduces some inelegant features compared to standard dynamic FEM: the VIV-Gauss point must be provided with , , and and a flag showing whether a call is made at a step or within a refinement iteration.
Note that when the VIV model is integrated in the dynamic FEM computation, in the way detailed in this section, the recursive Laguerre filter presented in Section 4.6 is not used. In this work, the filter was used only to process velocities time series and product Laguerre coefficients for the training of the rotatron model. The filter provides a very efficient way to process whole time series.
8. Results
8.1. Training
The Norwegian Deepwater Program was a research effort in which reduced scale tests were carried out on long, flexible riser models, subject to uniform or sheared current [35]. The displacement histories thus acquired at 19 points along the riser model were later prescribed on short stiff cylinders, and the hydrodynamic forces acting on the cylinders directly measured [36, 37]. Some data from [36, 37] is used in this work. It consists of the displacements at 19 points along the NDP riser model, for 3 current profiles (Table 1), so a total of 57 short cylinder test runs. For each of the 57 runs, 100 instants are randomly selected, yielding a training set of the rotatron with 5700 “points”. Each “point” consists of two sets of Laguerre coefficients and the two components of the corresponding force (Figure 11).
The rotatron was trained using Laguerre polynomials, 200 neurons in the hidden layer, and 50 to 1000 iterations of the conjugate gradient optimization algorithm. Other settings have been studied. No precision improvement was obtained from increasing the number of neurons or from increasing the number of iterations. This may suggest that the optimization algorithm converges on a local minima. If that is the case, improvement would require the use of an optimization algorithm better suited to finding absolute minima in a “jungle” of local ones.
Figure 9 shows how the rotatron predicts the forces for the trajectories in the above-mentioned 57 runs of short cylinder tests. The predicted forces are compared to the forces acquired experimentally. The model's ability to predict these forces seems to be good, although we lack a good criteria to judge that yet.
Figure 10 provides a visualization of the different steps of the modelization process, and is hence a useful diagnostic tool. It shows stippled black line: the trajectory for which a force prediction is wanted,smooth black line: the Laguerre approximation to the above trajectory, used to enter the rotatron,stippled black arrow: the measured force for the above trajectory,smooth black arrow: the predicted force for the above trajectory,green: neighboring (in the sense of the rotatron distance, Equation (6.7)) trajectories from experiments, used in the training set (and corresponding Laguerre approximation, experimentally measured force and predicted force),red: same as the above after rotation and/or mirroring.
Figure 10 gives an indication of the quality of the Laguerre approximation, the adequacy of the training set for the trajectory at hand, the presence of contradictions in the training set near the trajectory at hand, the quality of the fit of the rotatron to the training data, and finally the quality of the interpolation between training points.
8.2. Dynamic Analysis of a Flexible Riser
VIV depends not only on current velocities, but also on the type of slender system they act upon. Tension, stiffness, damping, length, and boundary conditions affect the vibration and hence the velocity trajectories that appear in the vibration. Hence the database used in Section 8.1 to train the rotatron is specialized, not only to a few Reynolds numbers for the current velocity, but also to some extend to the particular model used in the NDP program. No study was carried out in this work to determine how much change in the structure this specialized rotatron can accommodate.
Hence, in order to test the performance of the present VIV model within a dynamic analysis, the simplest case was considered: the riser model used in the NDP testing program ([35], characteristics in Table 2) was modeled. The present method was used to analyze the test condition TN2370 (). Numerical results were compared with those obtained experimentally on the NDP model. A laptop, using one core of a dual core processor, took 20 to 50 seconds to compute one second of riser response.
In Figures 12 and 13, the horizontal axis is NDP-laboratory time and the vertical axis is the unscaled length along the riser. The upper subplot shows the response in line (IL) with the flow, the lower subplot the cross-flow (CF) response. The color codes the displacements, with the same color scale used in Figures 12 to 14.
Figures 12 and 13 are for test TN2370 (). The dynamic simulations captures the frequency doubling between CF and IL, as well as the instationary nature of the vibrations. Frequency and amplitude are adequately captured. The 6th mode's dominance of the CF vibration is correctly captured. The dynamic analysis assumes constant tension. By contrast, some small tension modulations seem to displace the position of the lower vibration node in the test. In test 2370, there is a marked tendency for CF vibrations to propagate downwards. In the analysis results, CF waves are of a more static character. The IL vibration in the analysis occurs at a higher mode (11) than in the test (9). This is best seen by counting red dots along diagonals, for example, starting from coordinate 38 m and time 17.9 s in Figure 12. As a consequence, and since, for a given propagation celerity, wavelength and period are related, the phase drift between IL and CF is of opposite sign in the analysis and the test.
It proved impossible to reproduce tests 2340 and 2030 in the same manner. The analysis quickly ends with IL and CF vibrations occurring at the same frequency, and in Figure 14, with in-line motions dominating.
9. Discussion
9.1. Force Prediction and Reynolds Interpolation
Figure 9 shows that given the velocities in the training data as input, the model allows to reproduce the forces in the training data, based on the velocity in the training data.
The same exercise is carried out with trajectories from another test, N2430 (not to be confused with TN2340), which was carried out in shear current ( to ). In comparison, the highest current velocity appearing in the training set is . In N2430, the forces are fairly well predicted at the lower current velocities (Figure 15), while the predictions are very poor at higher velocities.
This illustrates that the present model provides no mechanism to use forces from a test at a given Reynolds number, to predict forces in a situation with another Reynolds number (due to higher current velocity e.g.,). In contrast, a simplified harmonic method like VIVANA enters a “table” with values of reduced amplitude and frequency. The tables are assumed to represent VIV behavior over a relevant range of Reynolds numbers. Because the transformation from added mass and excitation coefficients into forces makes use of the current speed (which is proportional to the Reynolds number), this effectively provides VIVANA with a mechanism for interpolation over Reynolds number. The present model, on the other hand, uses a compressed form (Laguerre coefficient) of recent velocity history, scaled to Reynolds-like numbers. The scaling and data compression scheme used in the present model does not separate the classical Reynolds number as a separate variable. Instead, it is embedded in the compressed data as the average of the recent scaled in line velocity already provided to the rotatron.
It remains an open question whether one should at all try to add such a Reynolds interpolation scheme in models of the present type. In a range of Reynolds number values in which the wake flow remains similar, it is reasonable to assume that viscous forces increase linearly with velocity and inertial forces increase quadratically. For harmonic motions, the component of the force in phase with acceleration can be attributed to an inertial force, and the component in phase with velocity, to a viscous force. For more general motions, “phase” becomes a problematic concept, making it difficult to separate the force into viscous and inertial components that one can scale each according to their own law.
9.2. Dynamic Analysis
The comparison of Figures 12 and 13 is encouraging: when simulating exactly the system from which trajectories were acquired experimentally in the first place, the simulation comes strikingly close to the experimental results. It is understood that a simulation is good if it captures the statistical properties of the real response. With VIV being chaotic, there is no hope to reproduce exactly any given realization of the response.
Figure 14 shows the result of a simulation of a case which is also directly represented in the training data. The simulation fails, in the sense that IL and CF vibrations are simulated to occur at the same frequency. Two explanations are proposed.
The quality of the force prediction by the rotatron, illustrated in Figure 9, was declared “satisfactory”, but this constitutes only an observation that the fitting procedure is operating. On what criteria should one judge that the fit is adequate? When training was carried out, the average norm of the difference between the training force and the predicted force was found to be about 30% of the average norm. For a control group composed of data points from the same experimental database, but not used in training the perceptron, the same ratio was about 40%. These numbers are high and can be reduced to some extend by increasing the number of hidden layers and of training iterations, but this was not found to yield better simulations. Further specializing the perceptron (by training it with data from test 2030) did not lead to successful simulations. One study which might help to understand the observations would be to find the corrective forces that need to be added to the forces predicted by the model, to force the simulation to track the motions observed during test 2030. This can be achieved using dynamic inverse FEM analysis [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 of
One issue that was explored was the adequacy of : the rotatron was trained with . This corresponds roughly to 1/4 of a cross-flow oscillation period in reduced scale for test 2370, and to a smaller fraction for other tests at lower current velocity. It could be argued that this fraction becoming too small could be the cause of the failure of analyses at lower velocities (Figure 14). The rotatron was hence trained again using a suitably increased value of . This was done twice, once with the full training set and increasing the number of Laguerre polynomials to , and once with only the experimental data from lower currents and an unchanged number of polynomials (). Neither rotatrons allowed to perform a successful simulation for the lower current velocities.
9.4. Tension
In Figure 12, one can observe modulations of the positions of the vibration nodes (in particular on the cross-flow graph). One possible explanation for this would be that in the physical system, the tension is modulated by the vibration. This effect, if present, is not captured by the numerical model, which assumes constant tension. The method presented here is designed for use in a nonlinear analysis. However, in this research, time was saved by using a simpler linear structural model.
10. Conclusions
A model for the prediction of VIV forces given the history of velocity of a cylindrical cross section relative to the undisturbed fluid has been developed. The model is closely related to Wiener-Laguerre filters: the recent history of velocity is represented by the coefficients of a Laguerre polynomial series. These coefficients are then used to enter a memory-less nonlinear interpolation function, in this case, a custom made neural network in which some relevant symmetry properties were “hard-wired”. The neural network was trained by using forces and displacements obtained in irregular forced motion tests on a short cylinder.
The proposed model operates in the time domain, making it well suited for integration into fully nonlinear analyses with unsteady currents. It further deals with in-line and cross-flow vibrations as one inseparable issue, which is arguably a necessity to improve on existing VIV models.
The model could provide a “good” reproduction of the forces in the training test, as well as a “good” prediction of forces for “comparable” trajectories. Where the model was queried with trajectories very different from those present in the training data, the model gave very poor results—as should be expected. Due to the limited amount of experimental data available, the present model remains quite specialized to a limited number of situations. What remains unknown at this stage is the size of the training set, and of the neural network model, necessary to create a model with some ambition of generality.
In some dynamic analyses of laboratory tests (NDP TN2030 and TN2430) with a long flexible riser, the numerical solution fell into an unphysical mode of vibration with the same frequency for in-line and cross flow vibration. On the other hand, for another case (TN2470), some unusually fine details were captured by the numerical model. From this, it is concluded that the concept has merit and deserves to be pursued.
At the same time, the present paper leaves several important questions unanswered.
The proposed method assumes that hydrodynamic forces can be predicted by the recent history of velocity of the pipe relative to the surrounding fluid. Fractal dimensional analysis in Section 6.3.2 suggests that this is, to some extent, the case. On the other hand, experimental evidence clearly shows that if a cylinder is forced twice to follow the same trajectory, the histories of hydrodynamic forces will not be exactly the same. Will adding a random force process to the output of the model be of any value in reproducing the dynamics behavior of a pipeline? If so, how should such a random process be studied and modeled?
What is the information content of the relation between velocities and forces? As the size of the experimental database is increased, will data points confirm each other and be interpolated with a relatively simple model (a rotatron with few tunable parameters), or will a landscape of ever increasing complexity appear, defeating the present approach?
The rotatron was found to be able to predict forces “reasonably” well, yet, when included in a dynamic model, the ability to reproduce the dynamics of the system was limited. Can the measure of force agreement, used in training the rotatron be improved? Would it be useful, as suggested to the author by Dr. Halvor Lie to emphasize a good prediction of the power of the force? More generally, how does one ensure that an approximate model reasonably reproduces the statistics of the response of the original nonlinear dynamics system?
Appendices
A. Conventions for Indexed Notations
The following conventions for indexed notations are used in this text. (1)By default, where an index appears more than once in a combination of products and/or divisions, a summation over the index is implied. If that index has a continuous range, then the “sum” is an integration over the range. Points 2, 3, and 4 specify exceptions to this rule. (2)Point 1 notwithstanding, if an equation is preceded by the symbol , followed by a list of indexes, then the listed indexes are not summed over. (3)Point 1 notwithstanding, if within an equation, there is a combination of products and/or divisions within which an index appears only once, then no sum over that index is carried out in the whole equation. (In any situation other than a simple term in the left-hand side, readability should be improved by using the symbol .)(4) Point 1 notwithstanding, if an index appears within an input to a function, and the output of the function is multiplied or divided by one or several terms that have the same index, then no sum within the input to the function is carried out on that index. (5)If an index of an argument to a function is within brackets, then the whole range of index values is used as input to one function evaluation. For example, refers to the evaluation at multiple locations () of a vector-valued () function of a vector (). (6)When the output of the function is shorthanded without explicitly writing its input, then the indexes of the input that are not within brackets are added to the indexes of the function. For example, can be shorthanded . (7)Derivatives of a function are noted with only the bracketed indexes of the input appearing under the fraction: . To refer to the value of that derivative for input , one writes .
B. Rotatron Gradients
The derivative of the force predicted by the rotatron, with respect to the Laguerre coefficients is needed in Section 7. With references to (5.3) to (5.8) that describe the rotatron, we can write with Here index ranges over two values (for two directions orthogonal to the cylinder), and is the other direction than .
The gradients of the rotatron with respect to its coefficients are also needed in order to compute the gradient of the target function with respect to the parameters , , and : with
C. Inverse of
The inverse of (Equation (7.19)), where is of the form with can be verified to be lower triangular banded, with terms on diagonal equal to
Acknowledgments
The author is indebted to Ida Aglen, Celeste Barnardo, Trygve Kristiansen, Carl Martin Larsen, Halvor Lie, Elizabeth Passano, Thomas Sauder, Wu Jie, Yin Decao and an anonymous reviewer for inspiring discussions on the subject of VIV, that either led to or helped the present research and writing. Thanks are extended to CeSOS (Centre of Excellence for Ships and Offshore Structures) at NTNU (Norwegian University of Science and Technology) and to MARINTEK for sponsoring this research. Very special thanks to Professor Carl Martin Larsen of CeSOS who, faced with a barrage of dubious ideas from the author, responded by making this research possible.