On a simple oscillator problem describing ice-induced vibrations of an offshore structure

In this paper, a new simple oscillator model is considered describing ice-induced vibrations of upstanding, water-surrounded, and bottom-founded offshore structures. Existing models are extended by taking into account deformations of an ice floe and a moving contact interaction between an ice rod, which is cut out from the floe, and the oscillator which represents the offshore structure. Special attention is paid to a type of ice-induced vibrations of structures, known as frequency lock-in, and characterized by having the dominant frequency of the ice forces near a natural frequency of the structure. A new asymptotical approach is proposed that allows one to include ice floe deformations and to obtain a nonlinear equation for the simple oscillator vibrations. The instability onset, induced by resonance effects for the oscillator and generated by the ice rod structure interaction, is studied in detail.


Introduction
Vertically sided, bottom-founded offshore structures occasionally experience sustained vibrations due to drifting ice sheets crushing against them. Usually, three regimes of interaction are distinguished: intermittent crushing, frequency lock-in, and continuous brittle crushing. In this paper, we introduce a mathematical model for a special type of ice-induced vibrations (IIV) of structures, known as frequency lock-in, and characterized by having the dominant frequency of the ice forces near a natural frequency of the structure. The first model of IIV was proposed in [1], where the ice failure was considered as a sequence of discrete events. This model was extended in [2] and more recently in [3], where randomness of the ice failure is taken into account. Mathematically, these models describe oscillators under an external-time dependent force, which simulates an action of the discrete events of ice failure. These models exhibit a resonance effect as a possible source of large IIVs. Other IIV models treat the ice failure as a continuous process (see, for example [4]) and can be applied for large ice velocities. In this case, the crushing ice force has a relatively low magnitude, and structure oscillations have low magnitudes and high frequencies. As was mentioned in [3], these situations are less important for safety applications since dangerous vibrations start with smaller velocities. Also some other models were proposed in [5][6][7][8][9][10]. Depending on the model, the interaction is described as a function of different sets of parameters related to the ice and the structure. In this paper, we propose a model extending the previous ones, in particular, those suggested in [1,2] and [3,10]. A novelty with respect to the previous investigations is that we study the ice deformations in more detail. We describe deformations of the ice rod taking into account possible viscous ice behavior (in [10] the small ice floe nonlinear deformation was taken into account) and a moving contact between the structure and the ice. For the simple oscillator-ice interactions, we also take into account extrusion effects. As a result, this leads to a difficult problem, which involves partial differential equations (PDEs) for the ice rod and ordinary differential equations (ODEs) for the oscillator. The main difficulty is that these PDE and ODE are coupled through boundary conditions for the ice rod deformations on a contact line between the ice rod and the oscillator. This contact line is unknown. The problem is solved by means of an asymptotic approach and by using a mechanical model for the ice rod deformation during the interaction with the oscillator. This asymptotic approach allows us to find a nonlinear ODE for the oscillator, where ice deformations are included. The resulting equation for the IIV terms involves many parameters, but a crucial parameter is the ice velocity v. The dynamics of such oscillator models can be studied by well-known methods (see, for example [11,12]). The following results are obtained by an asymptotic investigation and by numerical simulations. We show that ice rod deformation patterns arise which are generated by an interaction between the oscillator and the rod. This ice rod deformation is exponentially decreasing along the ice rod length, and moreover, the oscillator is forced to oscillate with a frequency close to the oscillators natural frequency . We obtain a plot describing how the oscillator's amplitude A depends on the ice velocity v. For some parameter choices, this plot shows a significant peak for the amplitude A. For small rod speeds v, we have a small oscillation amplitude A, as well as for large v. The height and width of the peak depend on the system parameters. These results essentially depend on the ice velocity v and are consistent with experimental data (see [2]). We find a new mechanism for the oscillator's behavior which can be described by a resonance between the external load, the oscillator, and the part of the ice rod (boundary layer) inducing an oscillator-ice rod interaction during the rod's motion with speed v. This paper is organized as follows. In Sect. 2 of this paper, the problem is formulated, and in Sects. 3 and 4, the problem is solved approximately by using asymptotic methods. Stability issues and resonances are discussed in Sect. 5, and finally in Sect. 6, some conclusions are drawn.

Statement of the problem
Following [1][2][3]5], we consider an offshore structure as a one-degree-of-freedom oscillator and the ice floe as a system of ice rods, whose properties include local failure (for details, see [5]). In particular, we consider a simple oscillator, which interacts with one of such rods. Considering only one rod follows from the theory which was proposed in [13]. As mentioned in [5], in order to simulate the behavior of the structure during the frequency lock-in regime, the condition was set that the ice should always be in contact with the structure, i.e., no gap should exist between the ice edge and the structure. On the basis of this requirement, the equation describing the simple oscillator dynamics is given by: where q = q(t) is the oscillator displacement, 2 = G M , where is the oscillator frequency, and M , and G are the mass and rigidity of the oscillator respectively, and α > 0 is a positive damping coefficient. For (1) the following initial conditions are used: The term μ in (1) defines a force which is acting on the oscillator due to the ice rod, and has the form where u = u(x, t) is the longitudinal ice rod displacement, E is an ice Young's modulus, F is the ice rod cross-sectional area, and δ 0 is the ice bulk viscosity coefficient. The term u x in the right-hand side of (3) defines the contribution of linear deformations, and the term u xt is the ice viscosity. It should be observed that the boundary conditions at x = q(t) also play an important role in other applications (see in [14][15][16][17][18]) when, for instance, longitudinal oscillations in axially moving cables are considered. The following equation describes the dynamics of the ice rod, which is defined on the semi-infinite domain: where u(x, t) is the unknown ice rod displacement, m is the ice rod mass per unit length, Q is a force per unit length occurring in the rod due to its sidesurface contact with other ice rods in the space around the rod in the ice floe that is considered in (4). The ice floe is drifting along the x-axis. The parameters β, δ 0 , and k 0 are positive. The parameter β is the ice shear viscosity, and the parameter k 0 characterizes the rod compression which is caused by stresses due to the ice rod compression in the transverse direction by other ice rods and thus defines "the load spreading capacity of the foundation" according to [13]. The parameter k 0 = Ehγ 2(1−ν 2 ) , where h is the ice floe thickness, the parameter ν is the Poisson ratio, the parameter γ is a coefficient determining the rate of decrease in the displacement u over the ice rod length and can be found experimentally (see also [13] for further details). In other words, the ice floe behavior can be modeled by a generalized spring and a generalized dashpot as suggested in [19]. The function s(t) describes the shift of the ice rod, and we suppose that s(t) is defined by where v > 0 is the relative ice velocity, and Here, t n are time moments when the ice rod crushes at x = q(t); d n are the lengths of ice blocks that split off, and H (z) stands for the Heaviside step function. The time moments t n are defined by the condition that is, the pressure p in the ice attains a critical level p c . The pressure p can be computed by the relation where t n is the moment of the previous break, and p 0 = (u x mc 2 0 )| t=0 is the initial pressure in the rod, where c 0 = √ r/m is the ice sound velocity. Therefore, the breaks are determined by the relation We introduce the following boundary conditions: The first one is a contact relation between the ice rod and the oscillator, and the second one is a radiation condition at infinity. The initial conditions are given by where φ j (x) are fast decreasing functions in x for x → ∞. For example, one can suppose that Notice that the differential equations, boundary and initial conditions can be transformed to a dimensionless form when we rescale the variables. For the rescaling, the following relations are used: It is obvious that the parameter is small. To simplify notations, we omit now the bars and obtain the final equations: and and the initial conditions are In the next two sections, we will approximate the solution of the initial boundary value problems (14)- (18).

Asymptotic solutions of the oscillator equation
The asymptotic approach to study Eq. (14) is well known; see [12,14,18]. Let τ = t be a slow time. It is assumed that α = ᾱ, where 0 <ᾱ < C, with C a positive constant independent of . Furthermore, it should be observed that the right-hand side of (14) is of order and depends on t and q. For that reason, we seek solutions q of (14) in the form where and where the amplitude A and the phase φ are unknown slowly varying functions in time. From (19), it follows that By substituting these relations into (1), and by taking together terms of order , one obtains the following equation for q 1 : where and Finally, it follows from (23) and (24) that the following system of equations for the amplitude A and the phase φ are obtained: and We investigate this system in the next section. To simplify the formulas, we will use the notation for the averaged quantities.

Assumptions
The aim of this subsection is to express the displacement u(x, t) in q and obtain an equation involving q only. We use the following assumption: i.e., the natural frequency of the oscillator is small with respect to the cutoff frequency of the ice rod. We also assume that all coefficients associated with friction and damping effects are small, i.e., It is useful to introduce two other small parameters To find u, we first define an auxiliary function V (t) as a solution of the following second-order ODE: We seek u in the form u = V (t) +ū, whereū satisfies the following equation, and the boundary and initial conditions: where V = dV /dt. In the next subsections, we determine the auxiliary function V and construct an asymptotic approximation forū.

Computation of V (t)
To simplify computations, let us assume that d n = d.
Then, s(t) is a periodic function with period T = d/v. We use a Fourier expansion for V . Consider the plot of s(t) (see Fig. 1). In this plot, a sawtooth curve is presented. Computing the Fourier coefficients of s, we have where T is the period and where ω = 2π/T = 2πv/d is the frequency of ice floe splitting off. Therefore, where the functionṼ (t) is the homogeneous solution of (31 ) and decreases in t with an exponential rate. The Fourier coefficients in (38) are defined byV n =ŝ The coefficientV 0 = d/2 is the mean shift of the oscillator. Note that the termṼ (t) is exponentially decreasing in t, and therefore, we can remove that term for large times t is small, then V n is large for n ≈ k 0 /ω. This means that resonances occur induced by the periodic breaking-off of ice from the ice rod. These results forV n will be used in the next subsections, and we will see that these resonances can create an onset for instability in the oscillator dynamics.

Asymptotics forū
Let us introduce a new variable: where L denotes the differential operator defined by and g is a functional of u and q defined by with The boundary conditions becomẽ and the initial data take the form Due to our assumptions (28) and (29) on the parameters , β and δ 0 , the term g is small. Thus, we can construct an approximation of the solution to the initial boundary value problems (IBVP) (40)-(45) as follows.
We setũ =ũ (0) +ũ (1) +ũ (2) , whereũ k are solutions of the following IBVPs: Let us findũ 0 . We write this function as a sum where and Problem (57) can be solved by using the Fourier decomposition method, yielding where the coefficients U n are solutions of the following BVPs: whereV n is defined in (39), and the complex frequenciesω n are defined bŷ In order to satisfy the boundary condition at y → +∞, we choose the complex roots such that Reω n < 0. In order to find an asymptotic approximation for U (q) , we use the assumption (28). Due to this assumption, the free oscillations of the ice rod defined by relation (56) have frequencies larger than k 1/2 0 , while q(t) oscillates with the frequency approximately equal to . Therefore, to obtain an approximation for U (q) , we can remove the derivatives with respect to t in the righthand side of eq. (56). That procedure gives the asymptotics where λ is defined by (30), and U (y) satisfies and In (65), as usual the signs of the complex roots are chosen such that the real parts ofω are negative. Then, for 2 k 0 (when the ice rod speed v is small), we havê To find U n (y) on the interval y ∈ (0, ∞), we perform the following standard procedure. Let us replace the conditions (60) by We solve the corresponding BVP problem for each L and then let L → +∞. Then, one obtains U n (y, L) = −V n exp(ω n y) +V n exp(ω n (L − y)).
The second term vanishes as L → +∞ for all y ∈ [0, R], where R is fixed. As a result, we obtain the following solutions of the above IBVPs: wherē These representations are valid for all compact intervals, y ∈ [0, R]. Note that the oscillator and the ice rod interact at the boundary y = 0 ; thus, these formulas allow us to describe this interaction. As a result, we find that The Fourier series in the right-hand side of the (71) converges sinceV n is O(n −2 ). Note that for small velocities v one has This solution describes a localized slowly oscillating oscillator. For large velocities, one can get This solution has two terms. The first term describes a fast exponentially decaying (in y) function generated by the ice-oscillator interaction which oscillates with the frequency . The second term describes a wave slowly decreasing in y = x − q. The wave propagates along the ice rod from the rod edge to x = ∞. The rates of decrease in space of the ice rod structure interaction and the wave are different. The first rate is proportional to √ k 0 and the second one to β/2. Our next step is to find the functionũ (1) , which is the solution of the nonlinear initial boundary value problems (48) and (50). Equation (48) forũ (1) is nonlinear and extremely complicated to solve in general. However, if we take into account that A 1 and 1, i.e., consider small oscillations with a small frequency, then we observe thatũ (1) ,ũ (2) ũ (0) , and thus we can removeũ 1 andũ (2) in g [ũ, q]. Furthermore, to obtain an asymptotic approximation forũ (1) , we can repeat the analysis as given earlier. We again use assumption (28), and therefore we can suppose that in (48) the expression for g [u (0) , q] is dominated by the frequency in the functions q, q t , q tt . Then, according to (41), a natural approximation forũ (1) is: where the functions G i (y, t) are defined by LG i = g i . Using the Fourier decomposition forũ (0) and the relation (42), one obtains nV n exp(ω n y + iωnt).

Deformation at the edge of the rod
In this subsection, it will be shown which ODE the oscillator displacement function q(t) satisfies. At y = 0, the boundary condition (14) becomes forũ(y, t): The relations obtained in the previous subsections show that the deformationũ y at y = 0 has the form u y (0, t) = S(t, q, q t , q tt ), where S(t, q, q t , q tt ) = Reωq + F 0 (t) where and As a result, we obtain and We also need the expressioñ which has the form S(t, q, q t , q tt ) = (G 10 +G 11 q +G 12 q t )q 2 t + (G 20 +G 21 q)q tt q t wherẽ and Using (86), (88), (89), and (90), one obtains the following equation for q(t): where ξ = S(t, q, q t , q tt ) + δ 0 dS dt −S(t, q, q t , q tt ) . (98) In the next section, the ODE for q(t) will be studied by using the averaging method as presented in Section 3 of this paper.

Equations for A and φ
In Sect. 3, it has been shown that the solution q(t) of (97) can be approximated by q 0 = A(τ ) sin( t + φ(τ )), where A(τ ) and φ(τ ) have to satisfy (25) and (26), respectively. From Sect. 3 and from the "Appendix," it follows that A(τ ) has to satisfy: where the coefficients D i are also given in "Appendix." The analytical and numerical analysis of the expressions for D i shows that the D i 's have the following properties. The coefficient D 3 is proportional to the small parameter δ 0 and is negative. For β δ 0 and √ k 0 β, the coefficient D 0 has a sharp maximum when ω ≈ , i.e., when we have a resonance between the external load and the oscillator. The other coefficients D 1 , D 2 show weaker resonances for 2ω ≈ and 3ω ≈ . Moreover, for β δ 0 and √ k 0 β, there exists a weak resonance between the oscillator and the ice floe structure for ω ≈ √ k 0 . Note that if β δ 0 , then this last resonance is much stronger than the resonances between the oscillator and the external load (ice crushing). For the phase φ, one has (see also (26)): where the right-hand side does not depend on the phase φ. Equation (99) can be studied completely in an analytical way. Consider Fig. 2 for the possible plots of the polynomial P(A, v) as a function of A. The polynomial can have either three real roots, two real roots ( of which one consists of two coinciding roots) or a single real root a 1 . In the third case that single root is a global attractor, i.e., A(τ ) → a 1 (v) as τ → +∞. In the first case if we have the roots a 1 (v) < a 2 (v) < a 3 (v), the solution A(τ ) goes either to a 1 or to a 3 as τ → +∞.
If A(0) = 0 and a 2 > 0 then A(τ ) → a 1 . If A(0) = 0 and a 2 < 0, then A(τ ) → a 3 . Note that a 3 has the order δ −1 0 1 for bounded D i and a 1 = O(1). So, when a 2 changes sign, we can observe a transition from small oscillator amplitudes to large ones. Therefore, Eq. (99) can describe dynamical bifurcations as D i (v) changes in v. A geometric model of these bifurcations is shown in Fig. 2.
The amplitude plots are given in Figs. 3 and 4; that is, in these figures, the stable equilibrium amplitudes A of Eq. (99) are given as function of v. The first plot shows the case k 0 2 , and the second one corresponds to a more interesting situation when k 0 is close to 2 . In the second case, we have a more complicated resonance between the external load, the rod, and the oscillatorice structure localized at the edge of the rod, and defined by (71). This resonance, in our opinion, leads to a sharp drop in amplitude as observed in experiments [4].

Discussion and conclusions
In this paper, a simple oscillator model is proposed to describe a continuous contact between the offshore structure (i.e., the oscillator) and the ice floe. It is a Fig. 5 Dependence of the oscillator displacement on time. Comparison of an asymptotic and numerical solutions characteristic property for IIV regimes describing continuous brittle crushing and frequency lock-in and, to a certain degree, intermittent crushing. The analytical methods as proposed in this paper allow analyzing the dynamic behavior of the oscillator. To verify the obtained asymptotic approximations of the solution, it is firstly compared with the numerical solution of the system of equations describing the problem. To find numerical solutions of the equations, the finite differences method is used and a uniform spatial grid is introduced. A monotone finite-difference scheme of the second order is used. The scheme is a generalization of the well-known Crank-Nicolson method. To obtain the solutions, an artificial viscosity is introduced in the finite-difference scheme. The obtained finite-difference equations are solved by a sweep method (the Thomas algorithm). The numerical method has been tested on elasticity model problems and has been compared with the known analytical solutions. Smooth and discontinuous solutions have been considered. This testing proved the correctness of the used finite-difference scheme. The calculations (see its results in Fig. 5) show that the asymptotic (curve 2 ) and numerical solutions (curve 1) are close to each other at times of order 1/ . The model as proposed in this paper has been compared with some current IIV models. A separate problem is to choose the system parameters necessary for the calculations. Such a parameter as bulk viscosity (or a relaxation time value related to it) has been found experimentally in [20][21][22][23].
In the papers [20][21][22][23], the authorsevaluated the ice vis-  [24,25] have computed the ice viscosity characteristic parameters. Therefore, the authors of this paper decided to take the δ 0 values from [20][21][22][23][24] and determined the effect of the viscosity coefficient value variations on the oscillator displacements. For that reason, also some calculations in the lock-in regime have been presented. The results are given in Fig. 6. Curve 1 corresponds to the viscosity coefficient value δ 0 = 3 × 10 9 Pa s; Curve 2 corresponds to the value δ 0 = 10 9 Pa s; and Curve 3 corresponds to the value δ 0 = 12 × 10 9 Pa s. The results presented in Fig. 6 reveal that the variation of δ 0 by a factor 12 does not effect significantly the behavior of the curves. The calculations reveal that we have satisfactory results that agree well with the results obtained by using other IIV models and experiments, when the parameter value δ 0 is taken in agreement with the data from [21,24]. Because of that, we take δ 0 = 3 × 10 9 Pa s and β = 0.75 × 10 9 Pa s in our calculations. Figure 7a compares the oscillator displacement in the lock-in regime with the results obtained in [26] for a 0.023 ms −1 speed. Curve 1 corresponds to the displacements calculated by the model proposed in this paper, Curve 2 is built up by using the displacement values found in the model from [26], and Curve 3 is based on experimental data as given in [26]. The amplitudes of displacements as obtained in this paper are close to the displacements found in [26] and are in a good agreement with the experimental data. In Figure 7b, Curve 1 describes displacements obtained with the help of the model proposed in this paper for the continuous brittle crushing at 0.13 ms −1 ; Curve 2 describes displacements for the same regimes, but obtained by using the model proposed in [26], and Curve 3 describes experimental data from [26].The curves presented in Fig. 7 show that the displacement amplitudes calculated by using the model proposed in this paper are half the value as those in [26] for times t < 39 s but for times t > 39 s the difference in values is not so big. At the same time, the values of the oscillator amplitude calculated with the help of the proposed model in this paper and the experimental one are in a satisfactory agreement. In Figure 8a, the ice force, which is acting on the oscillator in the lock-in regime, is compared with the results obtained in [26] for a 0.023 ms −1 speed. Curve 1 corresponds to the values calculated by the model proposed in this paper, and Curve 2 is based on experimental data as given in [26] for a 0.023 ms −1 speed. In Fig. 8b, Curve 1 describes the ice forces as obtained with the help of the proposed model in this paper for the continuous brittle crushing at 0.13 ms −1 , and Curve 2 describes the experimental data from [26]. In Fig. 9a, Curve 1 depicts the oscillator displacements in the lock-in regime at 0.18 ms −1 as obtained by the model presented in this paper; Curve 2 and Curve 3 give the model results and experimental data from [27], respectively. For the lock-in regime presented in Fig. 9a, all curves are in good agreement. In Fig. 9b, oscillator displacement comparisons are presented in the intermittent crushing regime at a 0.02 ms −1 speed. In this case, the displacement amplitude as calculated by using the proposed model is in good agreement with that presented in [27]. In Fig. 10a, Curve 1 describes the ice forces in the lock-in regime at 0.18 ms −1 speed as obtained by the proposed model in this paper, and (a) (b) Fig. 9 Dependence of the oscillator displacement on time during frequency lock-in at an indention speed of 0, 18 ms −1 and intermittent crushing regime at a velocity of 0.02 ms −1 Curve 2 gives the experimental data as given in [27]. In Fig. 10b, the proposed model results and the experimental data are compared for the intermittent crushing regime at 0,02 msec-1. The amplitudes of the forces as obtained in this paper are close to the values found in [26,27] and are in a good agreement with the experimental data. In Fig. 11a, one can compare the results as obtained by using the model as proposed in this paper and the model as proposed in [3] for the lock-in regime. The figure shows that the displacement amplitudes as obtained by using the model proposed in this paper are less than the ones as obtained in [3]. Figure 11b shows results for the continuous brittle crushing. The amplitudes as obtained by using the proposed model are 20% less than the ones as obtained in [3]. Thus, based on the comparison with the models that were proposed earlier, we can conclude that the model as proposed in this paper can be used for the evaluation of the structure vibration amplitudes in the lock-in regime. For the other vibration regimes, the model results underestimate the values when compared to the results obtained by using other current models. The calculations show that the asymptotic solutions as obtained can be used for the primary evaluation of the behavior of a structure interacting with the ice floe. Of course, the simple model proposed in the paper needs updating and further improvement with factors that have not been taken into account. In this paper, a new model is proposed to describe the interaction between an ice rod and an oscillator. This model takes into account deformations of the ice floe. The model is analytically investigated by asymptotic methods. The model presented in the paper provides an explanation for the initiation of nonlinear dynamic oscillator frequency regimes as a result of its interaction with a moving ice floe. The second possible cause of the nonlinear dynamic frequency regimes is the nonlinear ice behavior (breaking stress/deformation rate), which has not been taken into account by the proposed model, but that can be introduced into an updated model in future research. The model describes the start of the amplitude sharp raise at the lock-in regime and the amplitude sharp drop when the regime is changing to the continuous brittle crushing. In the authors' opinion, the cause of the amplitude sharp rise at the lock-in regime is the oscillator resonance caused by the moving ice floe. The resonance observed depends not only on the oscillator parameters and the ice floe velocity, but also on physical-mechanical ice characteristics, in particular, on relations among such param-eters as a Young's modulus, an ice thickness, a contact area size, dynamic viscosity, ice shear viscosity and a parameter determining an ice rod contraction caused by other rods in the oscillator's area. Apart from the ω = frequency, resonances are possible at 2ω = and 3ω = frequencies that contribute to the amplitude growth. Besides, at certain conditions given in the paper, a resonance between the oscillator and an ice rod part adjoining the contact area is also possible. A complex dynamic regime occurs in this zone, including two processes. The first process is localized frequencies of the rod end at an oscillation frequency. The second process is a slowly dampening wave depending on x-q variables and propagating along the ice rod. Thus, these localized frequencies lead to amplitude disruption and transformation to a continuous brittle crushing. This complicated resonance between the external load, the rod, and the oscillator-ice structure located at the edge of the rod can occur and can lead to sharp drops in the amplitudes as also has been observed in experiments (see [4]). The model used in [24][25][26] represents a structure as continuous and, therefore, is more general. We would like to point out that the model proposed in this paper can be further updated. It might be necessary to consider a continuous structure instead of the oscillator, to take into account water and wind effects, loss of contact between an oscillating structure and ice in intermittent crushing regimes and other factors.
Relations for D 0 , D 1 and D 2 Assuming that δ 0 > 0 is small, to compute D 1 and D 2 , we take into account only the averages, which have order O(1) and do not involve the parameter δ 0 . We observe that the coefficient D 1 is defined by the average and To compute these averages, let us introduce the auxiliary functions J + m,n = We see that these expressions describe resonance effects at nω ≈ m . Then, we find that and D 2 ≈ −0.125 2 n∈Zω nVn (J + 3,n + 3J + 1,n ).
Similarly, for D 0 one has