An analysis of parametric instability of risers

An analysis of parametric instability of risers

Fabio S. Prado^*; Fernando Y. Sakamoto; Carlos E. N. Mazzilli

Escola Politécnica da Universidade de São Paulo, Av. Prof. Almeida Prado, trav 2 no. 83, CEP 05508-900, São Paulo, SP, Brasil

ABSTRACT

In this work, three catenary riser models subjected to harmonic oscillations are studied. Two are finite-element models, one studied with Orcaflex, an offshore marine system analysis software, and another one with Abaqus, a generalist structural analysis software. The third model is an analytical reduced-order model that represents only the touch-down zone. The aim of this study is to discuss the feasibility, potentialities and limitations of the analytical model in confrontation with the specialist and the generalist softwares for the analysis of risers, under conditions of parametric excitation and unilateral contact at the seabed.

Keywords: Finite-element modeling, reduced-order modeling, riser, parametric instability, nonlinear dynamics, unilateral contact.

1 INTRODUCTION

Risers are extremely slender pipes that convey oil and gas from the seabed up to the offshore platforms. There are different configurations of risers, among them vertical, catenary, lazy-wave and steep-wave risers, but only catenary risers are studied in this paper.

Several sources of nonlinearities are present in a catenary riser analysis, making it a very complex system to study. Due to its slenderness, most part of the riser has cable behavior, where the equilibrium configuration depends on the tension and the tension depends on the equilibrium configuration. Only in two regions the bending behavior is relevant, at the hang-off and at the touch-down zones. Other important geometric nonlinearities are the unilateral contact and friction between riser and soil, and the hydrodynamic interaction between structure and fluid that generates nonlinear damping and lift forces.

The nature and amount of dynamic loads that the riser is subjected to makes the problem even more complex. One may have combination of platform movements, waves, internal flow, currents acting in different directions, levels and intensities, which in turn trigger VIV's (vortex induced vibrations). In this nonlinear scenario, interesting dynamic phenomena may appear, like parametric resonances, which is the main subject of study in this paper. Unlike the classical resonance, where the excitation frequency is equal or close to a natural frequency of the structure, the parametric resonance occurs when some parameter of the equation of motion varies periodically with time, for instance, the stiffness. In this case, one may have, apart from the 1:1 resonance, other ratios, like 2:1 where the excitation frequency is twice of the natural frequency. Parametric excitation can be better understood considering the so-called Mathieu equation:

This is a linear undamped equation of motion where the stiffness varies with time. For some pairs of the control parameters (δ, ε), where δ depends on the ratio between a natural frequency and the forcing frequency, set as Eq. (2), and depends on the forcing amplitude, the solution of Eq. (1) is seen to grow without bound; for other pairs of the parameters, the solution is seen to be limited. The Strutt's diagram shown in Figure 1 indicates the regions of unbounded solution (hatched areas), when parametric instability is said to occur. Of course, when nonlinearities and dissipative effects are taken into account, the instability zones can change and post-critical steady-states can also appear.

In the riser problem, the time variance of the stiffness is attributed to the axial force variation due to the riser vibration. So, to analyze the problem of a catenary riser subjected to parametric excitation, three models were considered. Two are finite-element models (FEM) and are handled by commercial softwares, Orcaflex 9.5 and Abaqus 6.10. The former is a specialist offshore/marine system analysis software and the latter is a generalist structural analysis software. The third model is an analytical reduced-order model and represents only the touch-down zone. It has only one degree of freedom, whose modal variable is the horizontal displacement of the touch-down point (TDP).

The purpose of this work is to analyze the feasibility, potentialities and limitations of the analytical model in confrontation with the specialist and the generalist softwares. In order to enable this comparison and the better understanding of the phenomena, simplified cases are studied. The models are 2D, subjected only to harmonic heave translation (vertical motion of the platform).

In the field of analytical models of semi-infinite slender beams on elastic foundation with unilateral contact, which is in essence the touch-down zone of a riser problem, advances were obtained in the following mentioned works. In Demeio and Lenci (2007), vertical displacements were imposed at the suspended part of the beam and resonances were seen to exist only for conveniently normalized frequencies smaller than one. In Mazzilli and Lenci (2008), normal vibration modes were obtained using the method of multiple scales. Excellent matching was achieved with resonant frequencies, as presented in Demeio and Lenci (2007). In Mazzilli and Mansur (2011), dynamic axial forces were imposed at the suspended part and parametric resonances were observed. In Mazzilli et al. (2012), the problem presented in Mazzilli and Mansur (2011) was solved by the method of multiple scales. The present work considers displacements and dynamic axial forces simultaneously acting at the suspended part, as in Sakamoto (2013).

2 FINITE-ELEMENT MODELS

The Orcaflex model is taken as the reference one and will not be detailed because it is handled in a straightforward way by this specialist software. The following procedures are for the Abaqus model.

To reach the static equilibrium configuration of the riser, the analysis was split in two steps. The first step, depicted as "1" in Figure 2, consisted in applying the axial thrust together with the submerged weight to the riser defined by its section and material properties, as well as by the elastic contact elements, while it is still laying down on the seabed. This step is necessary for a better convergence of the next steps.

The second step aims at reaching the final configuration of equilibrium (catenary-like shape). For that, displacements are gradually applied to the end opposite to the riser anchor, until the final hang-off coordinates are reached, and then the quasi-static current loads are introduced. This procedure is illustrated as "2" in Figure 2.

After the static equilibrium has been determined, the vibration-mode shapes and frequencies of the riser can be evaluated. Then the excitation frequency is strategically chosen to be twice a certain natural frequency ω₀ and within typical heave frequencies. Next, the heave amplitude is imposed at the hang-off so that no undesirable dynamic compression occurs. The analytical formulation, in its present development, could not properly handle it, anyhow.

Eq. (1) can be generalized for a FE model of an undamped linear system under parametric excitation as shown in Eq. (3), where the geometric stiffness matrix [K_G]cos 2 ω₀t is seen to vary parametrically with time.

3 ANALYTICAL REDUCED-ORDER MODEL

3.1 Problem Formulation and Methodology

As previously mentioned, the analytical model represents only the touch-down zone. Only the vertical dynamics is taken into account, although horizontal motion along the riser can be determined once the displacement of the touch-down point (TDP) is known. A semi-infinite beam is used to model the touch-down zone, under the hypotheses of the Bernoulli-Euler theory. The touch-down zone has two parts, one is suspended and the other rests on a Winkler-type elastic support working only in compression and with stiffness coefficient µ, as depicted in Figure 3. To set the extension of the suspended part taken into account in the analytical model, it is assumed that the bending behavior (stiffness EI) is relevant until a distance 4λ from the TDP, λ being the bending length given by Eq. (4), where T_TDP is the static effective tension at the TDP. There is no analytical proof for this assumption. It is simply based on the observation of responses of numerical and experimental analysis. A study about the TDZ's bending behavior can be found in Pesce (1997). So, the distance 4λ is used to define the position of point "O", where axial forces T(t) vertical displacements w(t) and phase angle between them are assumed to be known. The data at point "O" and the static effective tension at the TDP are provided by a global riser analysis made in a FEM model, using Abaqus or Orcaflex, for example.

The equation of motion that describes the problem is given by Eq. (5), where ρ is the mass per unit length, p is the submerged weight per unit length and H is the Heaviside function. As the position of the TDP x_c(t) changes with time, it is a moving-boundary condition problem.

Considering the relations below, it is possible to re-write Eq. (5) in dimensionless form:

In order to transform the problem into a fixed-boundary-condition one, the following variable transformation is introduced:

The price that is paid with such a transformation is that the equation of motion becomes strongly nonlinear due to spatial and temporal differential-operator transformations. For details, see Sakamoto (2013).

The methodology adopted to solve this complex equation is similar to the one proposed in Mazzilli and Mansur (2011). First, the solution of Eq. (9) is written in the form:

where δ(z, τ) is a perturbation, caused by the dynamic loading, around the 'quasi-static' configuration û(z, τ) that depends on time due to displacements imposed at point "O". This general solution is incorporated in Eq. (9).

The 'quasi-static' solution itself is obtained from Eq. (9) under the assumption of null time-derivatives for u and c. To account for the tension decrease in the suspended part as the TDP is approached, and in the contact zone, due to friction, decay-law functions are assumed and incorporated in Eq. (9) in an "ad hoc" manner.

To transform the continuum model into a reduced-order model (ROM) with only one degree of freedom, whose modal variable is the horizontal displacement of the TDP, a nonlinear Galerkin method is used, so that Eq. (9) is projected on one nonlinear vibration mode obtained, for the sake of simplicity, for the simple-bending case. The procedures to obtain these nonlinear vibration modes can be seen in Sakamoto (2013). To complete the reduction, virtual works are imposed to be the same in the continuum model and the ROM.

To eliminate the spatial dependence of the coefficients of the reduced-order-model equation of motion that are functions of the selected vibration mode, the tension decay-law function and the quasi-static configuration, rotation and curvature, they are numerically integrated along the riser length. The integration of these functions is set from the point "O" up to a point at the supported part that is far enough from the TDP for the functions to take values close to zero. This way, a ROM differential equation is obtained and solved numerically to give the response in the time domain. Obviously, the whole process is not fully analytical. The term "analytical model" was used to distinguish it from the high-hierarchy numerical models.

3.2 Quasi-static Solution

The 'quasi-static' configurations for the suspended part û_e(z, τ) and the seabed-supported part û_d(z, τ) are established for the simple-bending problem, but a correction factor ψ₁ can be used to assure a good matching with the configuration given by the FEM model. The dimensionless 'quasi-static' co-ordinate of the TDP is c₀(τ) in a certain normalized time τ. The dimensionless static co-ordinate of TDP is c₀(0) and c₁ is the 'quasi-static' dimensionless displacement of TDP caused by imposed displacements at point "O".

Likewise, correction factors ψ₂ and ψ₃ are used so that rotation and curvature functions (û',û") are close to those of the FEM models.

Knowing that the modal variable is the horizontal displacement of the TDP, the displacement imposed at point "O" is transferred to the TDP using Eq. (11). For z = -1 and τ = 0, the static co-ordinate of point "O" û₀ is found, as depicted in Figure 4. Then, the co-ordinate û₀ + û_1max, as obtained from the FEM model, is put in place of û_e(-1, τ_max) in Eq. (11), and the new position of TDP c₀(0) + c_1max is found. τ_max is the time in which the maximum co-ordinate is reached. The same procedure is done for û₀ + û_1max and c₀(0) + c_1min. The amplitudes c_1maxand c_1min have different values because Eq. (11) establishes a nonlinear relationship between û_e and c₀. But, for the sake of simplicity and knowing that c_1max and c_1min are small compared to the static position c₀(0), it is reasonable to consider for the quasi-static amplitude at the TDP c₁ motion the average between c_1max and c_1min. This allows considering a simple harmonic quasi-static motion at TDP, corresponding to the imposed vertical displacement at point "O.

As mentioned before, spatial dependence is eliminated by numerical integration of the functions that depend on the coordinates along the riser length. Looking at Eq. (12) it is noticed that is not possible to detach time and space variables, since time-dependent terms are inside the arguments of the sine, cosine and exponential functions. The solution was to expand these functions in power series. But to get a good approximation, this procedure generated super-harmonics up to the 37^th order and hundreds integrals to solve. Fortunately, a sensitivity analysis of the super-harmonic terms showed that they have low importance in the response, because they are associated with the term c₁ⁿ, that gets ever smaller compared to c₀ⁿ, with the increase of n, which is a positive integer. In Sakamoto (2013), a fourth-order model presented the same response of a first-order model. Therefore, only the first-order terms were kept:

The same procedure was adopted for rotation and curvature. The terms û_e0, û_e1, û_d0, û_d1, as well as the terms for rotation and curvature, are huge expressions and are not shown here for brevity, but can be found in Sakamoto (2013). It is seen that 6/c₀(0) is the normalized distance from the TDP where these functions take almost-zero values.

3.3 Tension Decay Law

Two exponential-decay laws are adopted for tension along the riser, one for the suspended part and another for the contact zone. The use of such laws is interesting because exponential functions also appear in the static and modal solutions. Its dimensionless form is written as Eq. (16).

For the suspended part, the decay rate a_e needs to be calibrated andis set so that the tension at point "O" and at the TDP are equal to the tensions given by the FEM model. This assumption is reasonable since it is a local model of the TDZ.

For the contact zone itself, the decay rate a_d is chosen so that the tension takes the value of 1% of the tension at the TDP, at a distance L_k from the TDP, where L_k is determined for a null tension under the Coulomb friction hypothesis (with friction coefficient k).

The dimensionless decay rate is:

3.4 Modal Functions

The modal functions for the case of simple bending, Eq. (21) and Eq. (22), were obtained via the method of multiples scales in Mazzilli and Lenci (2008). Parameters a_e, a_d, C_i, β* are terms that depend on and , which is a dimensionless natural frequency, and are presented in Sakamoto (2013).

3.5 Reduced-Order Model

Incorporating the 'quasi-static' tension and modal functions in Eq. (9), it can be seen that:

To solve this equation and reduce the infinite number of degrees of freedom of the continuum model to a single-degree-of-freedom reduced-order model, the nonlinear Galerkin method is used. As the horizontal displacement of the TDP is chosen to be the modal variable and only one vibration mode is used for projection, the following expressions are used:

The displacement of the TDP c = c₀(0) + U, refers to the overall effect of both the imposed motion at point "O" (as partially considered in Eq. (13)) and the imposed tension variation at point "O". To alleviate the notation, the term c₀(0) is now simply termed c₀:

To complete the reduction procedure, each term of Eq. (25) is considered a modal force F(U, , Ü, τ) and its virtual work is calculated as δw = F(ζδU), where δU is a virtual displacement:

To eliminate the space dependency, as previously mentioned, auxiliary constants I_i are used, defined as integrals evaluated from point "O" till a dimensionless distance L_z, where the space-dependent functions reach values close to zero. Note that no damping has been considered so far in Eq. (26). Equivalent linear viscous damping is now incorporated in an "ad hoc" manner to account for energy dissipation. So, Eq. (26) takes its final form:

where:

Note that Eq. (27) has parameters that vary harmonically with time, as the classic Mathieu equation does. So, parametric resonances can be expected in certain scenarios. The dimensionless natural frequency of the reduced-order model, now considering the geometric stiffness effect, can be calculated:

4 CASE STUDY AND RESULTS

A steel catenary riser subjected only to heave motion is modeled with Abaqus, Orcaflex and ROM. Three excitation periods are studied. First the static equilibrium configuration, rotation, curvature and tension decay are compared. Then, the vibration modes calculated by the Abaqus are presented as well as the displacements, velocities and phase diagram for the TDP. For the ROM case, several simulations can be run, and a color map showing the intensity of the peak-to-peak displacement of the TDP is constructed varying excitation frequencies and dynamic tensions at point "O". Looking at this map, it is possible to see which pairs of frequency and dynamic tension trigger parametric resonance.

4.1 Geometric and mechanical properties of the riser

- External diameter (Dex) = 0,2032 m

- Internal diameter (Din) = 0,1651m

- Young's modulus (E) = 200 GPa

- Axial rigidity (EA) = 2204,18 MN

- Bending stiffness (EI) = 9443,3 kN.m²

- Submerged weight per unit length (p)= 727 N/m

- Submerged mass per unit length (considering additional mass) (ρ)= 141,24 kg/m

- Water depth = 158m

- Soil stiffness (Φ) = 10 kN/m/m²

- Z co-ordinate at the hang-off = 158 m

- X co-ordinate at the hang-off = 199 m

- Length of the riser = 300 m

- Angle at the hang-off = 77,5º

4.2 Static Analysis

Figure 5

The static tension at TDP given by the Abaqus model is 28,47kN. In possession of this result, the bending length can be evaluated by Eq. (4): λ = 18,2m. So, the horizontal length of the suspended part of the ROM is 4λ = 72,8m. The tension at point "O", from the Abaqus model, is T₀ = 77,30kN and its height is 66,74m.

The tension curve of the ROM differs from the curves of Abaqus/Orcaflex, mainly in the supported part, because friction between riser and soil was considered in the ROM, but not in the Abaqus/Orcaflex models. This friction only causes the tension to decay in the ROM. It doesn't imply any additional energy dissipation in the reduced-order model.

One of the boundary conditions imposed to the ROM was that the point "O" is hinged (allows rotation). This is confirmed by the curvature equal to zero at point "O".

4.3 Dynamic Analysis

Figure 6

Three cases of heave excitation are studied here: one for excitation period T = 2,9s and heave amplitude z₁ = 0,5m, another for T = 5,2s and z₁ = 1,0m and one last for T = 10s and z₁= 1,0m. The overall equivalent damping ratio adopted in all cases is x = 10%. After running these cases in Abaqus, the following data at point "O" are available for the ROM:

Table 1

The second nonlinear vibration mode of the TDZ was used for projection. The normalized natural frequency evaluated by Eq. (28) is ω₀ = 0,2935, which corresponds to a natural period of T_n = 5,64s, very close to the second vibration mode given by Abaqus.

In Figure 7, the responses for the case T = 2,9s are presented. The frequency ratio for Abaqus/Orcaflex is Ω/ω_n = 1,99 with respect to the second mode and Ω/ω_n= 0,96 with respect to the fourth mode. For the ROM, the ratio is Ω/ω_n= 1,94. In the three models parametric resonance are detected, as the period of the response tipically doubles to T=5,8s. The displacement of the TDP in the three models showed good agreement. The velocities, yet, presented some discrepancies. Abaqus and Orcaflex showed larger values than the ROM. That is because the displacements of Abaqus/Orcaflex contain higher harmonics, whereas ROM was obtained after projection on a single mode, namely the one under parametric resonance. Classical resonance 1:1 also occurs with the fourth mode and this is captured by Abaqus/Orcaflex, but not by MOR. Yet, when parametric resonance builds up, it prevails over the classical resonance. The displacement curves of Abaqus/Orcaflex also show small peaks with period T = 2,9s. When the displacements are derived with respect to time, sharp peaks of velocities appear. This doesn't occur in the ROM, since only one vibration mode is taken into account and all the energy of the system is directed to this mode. In the phase diagram it is seen that MOR's response does not show the "loop" displayed by Abaqus/Orcaflex, meaning that the response is dominated by the parametric instability (frequency Ω/2) and the forcing frequency (Ω) content is small. It is worthy to note that the forcing frequency is not in the 1:1 resonance with respect to the MOR, since this latter was obtained projecting the response onto the mode of frequency Ω/2. That is why the dynamic amplification associated to the forcing term with frequency Ω/2 is small and not perceptible in the MOR response. Yet, for the FE models, the same forcing term will be in the 1:1 resonance with respect to the mode with frequency Ω, which explains the (small, but not negligible) intermediate peaks in the FE time responses and the "loop" in the corresponding phase trajectories.

In Figure 8, the responses for the case T = 5,2s are presented. The frequency ratio for Abaqus/Orcaflex is Ω/ω_n = 2,00 with respect to the first mode and Ω/ω_n = 1,11 with respect to the second mode. For the ROM, the ratio is Ω/ω_n = 1,08. Classical resonance 1:1 occurs in the three models and the response period is the same as the excitation period. The displacement and velocity of the TDP in the three models showed good agreement. Parametric resonance 2:1 could have occurred, at least in principle, with respect to the first mode for the Abaqus/Orcaflex models, but it didn't happen, probably because the force magnitude was not enough to trigger it. This couldn't be detected in the ROM because projection was made onto the second mode whereas the important mode would be the first one.

In Figure 9, the responses for the case T = 10s are presented. The frequency ratio for Abaqus/Orcaflex is Ω/ω_n = 1,04 with respect to the first mode. For the ROM the ratio is Ω/ω_n = 0,56. Classical resonance 1:1 occurs for the Abaqus/Orcaflex models and the responses are amplified. No resonance occurs in the ROM and the responses are much smaller than presented by Abaqus/Orcaflex because projection was made onto the second mode whereas the important mode would be the first one. In this case, this particular ROM is totally inadequate to analyze the dynamics of the riser, for obvious reasons. Again, Abaqus and Orcaflex had good agreement between them.

Figure 10

In the 'displacement intensity' colored maps of Figure 11, the value showed in the square cells are the peak-to-peak displacements of the TDP for each pair of frequency and dynamic tension ratios. The tension axis has the maximum value equal to the unit because the ROM cannot properly handle dynamic compression. The black dot is where the case T = 2,9s was set. The red dot is where the case T = 5,2s was set. And the blue dot is where the case T = 10s was set.

The first map of Figure 11 is for damping ratio ξ = 10%, which is the damping adopted in the case studies. Looking at this map, it is possible to see which pairs of frequency and tension trigger the parametric resonance 2:1. Other resonance ratios are observed, like 1:1 and 1:0,5. Actually, there are several resonance ratios for parametric resonance, as showed in the Srutt's diagram of Figure 1. Resonances occur when δ = n², where is an integer and is defined as in Eq. (2). Transforming the parameter δ into the parameter Ω/ω_n, resonances are observed for values of Ω/ω_n = 2/n, i.e., Ω/ω_n = 2, Ω/ω_n = 1, Ω/ω_n = 0,667, Ω/ω_n = 0,5, and so on. For the case studies, the analogous to Strutt's diagram can be constructed, if Eq. (27) is linearized as follows and the damping is set equal to zero:

Figure 10 shows two Strutt's diagrams for the undamped case, one on the Ω/ω_n parameter base and another on the δ parameter base. Parametric resonances occur for the same frequency ratios presented in Strutt's diagram of Figure 1. Additionally, the complete nonlinear equation is represented in the same diagram by a curve that displays the limits above which responses could not be found by the ROM because of convergence problems. It is possible to see some correlation with the instability boundaries of the linearized equation.

The second map of Figure 11 is for damping ratio ξ = 5%. It shows how the parametric resonance zones grow as damping decreases. The third map is for damping ratio ξ = 0%. The red zone indicates the cases for which the ROM couldn't find a solution. The same zone is given in the first graph of Figure 10.

5 CONCLUSIONS

The static analysis shows that there is a very good agreement between Abaqus and Orcaflex. Correction factors applied to the static equilibrium configuration functions of the ROM leads to a very good fitting with Abaqus/Orcaflex curves.

In the dynamic analysis, the results showed that the three models could detect the phenomenon of parametric resonance. The responses of Abaqus for all case studies are qualitatively and quantitatively close to Orcaflex responses. This means that the generalist software was able to handle this very complicated problem, under parametric instability and unilateral contact conditions.

ROM responses showed good agreement with those of the Abaqus/Orcaflex for parametric resonance 2:1. Even when the FEM model revealed that the system was simultaneously under parametric and classical resonances, the response of the ROM remained valid, because parametric resonance prevails over the classical resonance. For classic resonance 1:1, ROM still predicts responses very well, provided the mode under resonance has been the one kept to project the response. Again, ROM can catch other parametric resonances, like 0,667:1 and 0,5:1, provided the adequate modes have been used to project the response. As these ratios are very close to each other, in the colored maps, a better discretization is required for better seeing the limits of these resonances.

Although simple, but not simplistic, ROM has shown that its outstanding potentiality is that it can run several simulations easily, generating auxiliary outputs like the displacement intensity maps, which help identifying conditions that can trigger parametric resonances. Analytical models are also useful to understand the overall qualitative behavior of the system, because in a low-dimensional phase space parametric analyses are much easier to be carried out. But to use ROM's, it is important to be aware of their limitations. For instance, the particular ROM studied in this paper is not suitable to properly analyze cases with dynamic compression.

Observing the displacement intensity maps, it is possible to see a clear parametric resonance 2:1 zone. Unlike the resonance 1:1, where the displacements grow gradually with the frequency and tension ratios, the resonance 2:1 shows a sudden change of behavior when the parametric resonance zone is reached, where the displacements grow sharply. For the same tension, the displacements in resonance 2:1 are much larger than the displacements under resonance 1:1. This sensitivity to frequency and tension ratios and the high-amplitude responses show how parametric resonance is dangerous. In the maps it is possible to see, as well, that the damping takes an important role in the resonance analysis. If the damping is increased, the phenomenon of parametric resonance quickly disappears.

Received in 14 Dec 2012

In revised form 21 Jun 2013

• Bathe, Klaus-Jürgen (1996). Finite Element Procedures, Prentice-Hall, Inc, New Jersey, 1037p.
• Demeio, L., Lenci, S. (2007). Forced nonlinear oscillations of semi-infinite cables and beams resting on a unilateral elastic substrate. Nonlinear Dyn 49:203-215.
• Gate, Chichester, West Sussex. 440p.
• Mazzilli, C.E.N, Lenci, S. (2008). Normal vibration modes of a slender beam on elastic foundation with unilateral contact. XXII ICTAM, International Congress of Theoretical and Applied Mechanics, Adelaide, Australia.
• Mazzilli, C.E.N, Mansur, A.L. (2011). Parametric resonance of slender beams on unilateral elastic support subjected to bending and time-varying normal force. 7th ENOC, European Nonlinear Dynamics Conference, Rome, Italy.
• Mazzilli, C.E.N, Mansur, A.L., Monticelli, G., Sakamoto, F. (2012). Parametric instability of slender beams with unilateral Winkler support: an application to riser dynamics. XXIII ICTAM, International Congress of Theoretical and Applied Mechanics, Beijing, China.
• Nayfeh, A.R., Mook, D. T. (1979). Nonlinear oscillations, New York.
• Pesce, C. P. (1997). Mechanics of cables and submerged pipes in catenary. Habilitation Thesis (in Portuguese), University of São Paulo, Brazil.
• Prado, F. S. (2013). Nonlinear modeling of risers by the finite element method, Master's Thesis (in Portuguese), University of Sao Paulo, Brazil.
• Sakamoto, F.Y. (2013). Dynamic modeling of the contact zone between riser and seabed under action of displacement and tension imposed, Master's Thesis (in Portuguese), University of Sao Paulo, Brazil.
• Wriggers, P. (2002). Computational Contact Mechanics. England: John Wiley and sons Ltd, The Atrium, Southern

Publication Dates

• Publication in this collection
  13 Dec 2013
• Date of issue
  May 2014

History

• Received
  14 Dec 2012
• Accepted
  21 June 2013