Numerical simulation of a wind excited quasi-periodic regime of a stochastic van der Pol-type equation

Abstract. The contribution regards a mathematical single-degree-of-freedom model of a slender structure vibrating in an air flow. Based on an experimental investigation, movement of such structures can be expressed by van der PolDuffing-type equations. Several particular configuration parameter settings for a white and non-white Gaussian random excitation together with deterministic harmonic forcing are considered and numerically analysed. The results support recently published analytic formulas.


Introduction
Based on experimental investigations, movement of slender structures vibrating in an air flow can be modelled by van der Pol-Duffing-type equations. Although the vortex shedding, which appears behind the structure, may induce a regular excitation, turbulence generated by an air circumfluence around the profile and its surface irregularities may cause random fluctuations of pressure. This effect particularly influences ice-accreted bridge or mast cables or power lines, [1,2]. A correct mathematical model thus has to consider both a deterministic (harmonic) excitation, caused by the vortex shedding, and a random excitation, produced by surface irregularities and turbulence.
The single-degree-of-freedom (SDOF) system following Fig. 1 is generally used for modelling the heave component of the structure. It is popular because of its simplicity and capability to characterise main dynamic properties of the structure. The basic SDOF model may be complemented by the pitch component, however, in the case of rotationally symmetric bodies is its influence only marginal.
It was shown that the presence of surface irregularities, e.g., ice accretion, introduces frequencies which slightly differ from the vortex-shedding frequency. The corresponding velocity spectral density of the random excitation component have the shape of the von Kármán spectrum [3], see also Fig. 2. Such spectral properties have to be respected also by random processes representing the stochastic part of the excitation. On the other hand, the real measurements confirm that the random component can be predominantly supposed Gaussian.
The recent paper of the authors [4] investigates several particular configuration-parameter settings of the non-linear system for a non-white Gaussian random part of an excitation process and characterizes corresponding response properties. The random part of the response was modelled using the Fokker-Planck equation (FPE). Solution of this equation is always problematic and only in very special cases can be performed analytically. A numerical approximation of its solution can employ, e.g., the method of moments, see [5], finite element  method, [6], or Monte Carlo and other simulation based methods, [7]. The authors in [4] expressed analytically its stationary solution in an exponential form based on a probability potential. Then the response characteristics were harmonically approximated and the corresponding partial amplitudes were determined using the stochastic averaging strategy. It was shown that only the coincidence of the excitation frequency (due to vortex shedding) and the system eigen-frequency can provide a stationary response.
Numerical simulations of the non-linear system from the aforementioned paper are summarised in this contribution. The input parameters were selected to be compatible with parameters which were used for illustration of the theoretically obtained results. The type of results presented in [4] enables a direct comparison only for the case when only the random component of excitation is present. For the case of a combination of harmonic/random excitation is taken into account, the analytic formula was derived to distinguish two distinct phase shifted components of the response. This information is unavailable in the numerical results, thus, in such a case only approximative comparison is presented.
The paper is organized as follows. First, after this introduction, the stochastic generalized van der Pol model is described. The numerical procedure used for solution of the Itô system is described in section 3, approximations of the probability density functions (PDFs) of the response are shown in section 4. Finally, section 5 concludes.

Mathematical model
A behaviour of the aeroelasic model shown in Fig. 1 can be represented by the generalized van der Pol equation. The quadratic or bi-quadratic damping terms may attain low or even negative values of damping and thus they enable to model the concurrent effect of forced and self-excited vibrations which appears when the two vibration frequencies are close to each other and lead to a quasiperiodic response. So that the corresponding stochastic differential equation for the quadratic damping used in [4] can be written as or, alternatively, for the bi-quadratic damping where u(t) denotes displacement [m]; ω 0 is the eigen-frequency of the adjoint The non-linear damping coefficients ν, ϑ follow from aeroelastic processes and depend on the stream velocity V. When V = 0 also the non-linear damping coefficients vanish. Coefficient η comprises an elastic conventional damping ratio, which is always negative (stabilizing), and a positive (destabilizing) aeroelastic part, which increases for increasing stream velocity. Consequently, for a certain stream velocity the total value of η may become positive and self-exited vibrations occurs.
The particular configuration of values η, ν, ϑ determine basic properties of the response. Particularly, the trivial state is unstable for η > 0 and stabilization of the response can occur: one stable limit cycle occurs when ν > 0, ϑ = 0 and two limit cycles (one stable, one unstable) take place if ν > 0, 0 < ϑ < ϑ max . See also [8,9] for detailed discussion. Without significant loss of generality, in the further work it will be assumed that ϑ = 0, i.e., only the quadratic damping is taken into account, Eq. (1). Although the simulation approach is the same regardless the damping complexity, the existence of a single stable limit cycle of Eq. (1) allows the integration procedure to use a simple method with a fixed steplength.
The main results of the paper [4] can be summarized as follows. The stationary response of the generalized van der Pol equation allows for characterization using a mono-harmonic approach where α c = α c (τ), α s = α s (τ) are slowly variable partial amplitudes. The FPE may be set up for the PDF p(a c , a s ) of their stochastically averaged variants, a c = α c , a s = α s . Solution to this FPE can be expressed by means of a stationary potential in an exponential form. When only a random part of excitation is assumed, P = 0, its general solution depends only on the absolute amplitude A 2 = a 2 c + a 2 c and can be written as Here the η, ν are damping coefficients, S = Φ(ω 0 )/(2ω 2 0 ) is the value of the spectral density Φ(ω) of the random excitation in the basic eigen-frequency of the system Eq. (1) and C is a normalization constant.
In the case of combined harmonic and random excitation, the solution of the FPE yields the following formula of the joint PDF of the partial amplitudes It is worthwhile to note that the character of the random excitation enters the theoretical formulas Eqs. (4,5) only through the value S , which corresponds to the value of the corresponding power spectral density at a single frequency. Vice versa, a particular shape of the spectral density function should not influence the resulting PDF of the response.

Numerical simulation
For the sake of simplicity, the numerical solution of Eq. (1), which is in fact a stochastic differential equation, will be performed using a simple digital simulation procedure. This facilitated approach is acceptable because the random part of Eq. (1) is represented only as an additive term. The integration procedure applied in simulation follows the standard fourth-order Runge-Kutta scheme, e.g., [10]. Due to presence of a discrete random process on the right hand side, the fixed steplength δt = 0.01s was used in all simulations. The same value of the random process was used also for the intermediate step (t = t i + δt/2) of the RK4 procedure for simplicity reason. Because the system eigen-frequency was set to ω 0 = 1rad s −1 , the chosen δt is sufficiently small for the assumed frequency content of the response. Moreover, the non-explosive character of the van der Pol equation (1), [11], allows the simple method to remain stable for a long integration time.
Two types of random excitation processes were used for simulation, both properly scaled to give the prescribed value S = 1 at ω = ω 0 : (i) a Gaussian white noise process which was generated using a normally distributed random number generator; (ii) a narrow band Gaussian process in the form of filtered periodic white noise generated with the method proposed by Shinozuka [12], i.e., the demanded signal is assembled of sine functions with random uniformly distributed phase angles. The spectral density of the narrow band process has a bell shape with the maximum at ω = 1 The simulation procedure starts from zero (unstable) initial conditions. The transient part of the response is discarded. Then the amplitudes of the stationary response are collected and analysed in the form of smoothed histograms which approximate the corresponding PDFs.

Effect of the random excitation
A rigorous comparison between an analytic PDF and an approximative histogram, which is based on a finite set of realizations, is always problematic. Positions of maxima of numerically obtained histograms may serve as a simple criterion. In addition to visual similarity of the histograms' shape is this test sufficient for the purpose of the paper. With regard to the analytical PDF given by Eq. (4) is the relation between the maximal probability amplitude  A m and damping parameters η, ν given in the form The analytical PDFs of the total amplitude computed according to Eq. (4) for different values of damping parameters are shown in Fig. 3. The normalization parameter C was determined for each curve, S = 1, ω = 1. The formula Eq. (4) corresponds to a rotationally symmetrical joint probability p(a c , a s ), which fact reflects, in this case, in a symmetry for positive and negative values of A. Four different cases for distinct values of quadratic damping ν are shown in Fig. 3; each one contains colour curves for different values of η = −0.6, −0.4, −0.2, 0, 0.2, 0.4. Only negative values of η are assumed for ν = 0, because the response amplitude diverges for ν = 0 and η > 0. The square plots in Fig. 3 show the dependence of the maximal probability amplitude on the damping coefficient η for individual values of the quadratic damping ν. The black dots in p(A) and A m plots mutually relate to values of η, the dashed lines in A m plot show the course of the analytic relation.
Numerically obtained approximations of PDFs of the response in the case of the Gaussian white noise excitation are shown in Fig. 4. The structure and used colours conform to that of Fig. 3. With an exception of the case ν = 0, the correspondence between analytical and numerical results is noticeable. On the other hand, positions of the maximal probability amplitudes A m fit well the theoretical dependence for ν = 0.33, the accuracy in cases ν = 0.67, 1.0 is acceptable. The inaccuracy of the case ν = 0 may be caused by a poor selectivity of the original system with respect to excitation by a white noise.
Visually better results are presented in Fig. 5 regarding the excitation by a band limited Gaussian white noise. Both the shape of individual histograms and the position of maximal probability amplitudes approximately fit the theoretically expected effects.

Combined random and harmonic excitation
The main advantage of the analytic result given by Eq. (5) is a possibility to distinguish effects of the booth amplitude components a c , a s . Unfortunately, the transformation into the amplitude and phase shift of the response is not straightforward. Identification of matching quantities using numerical procedures is hardly feasible. For this reason, the results presented in this paragraph do not relate directly to those theoretically obtained in [4]. Instead, two particular cases will be presented. Both plots in Fig. 6 show smoothed histograms of the response amplitudes computed for a combined harmonic and random excitation. The Gaussian white noise with spectral density Φ(ω) = 1 was considered in all cases for the random part of the excitation and the harmonic signal with the resonant frequency (ω = ω 0 = 1) for the deterministic part. Unlike plots in Fig. 3-5, only one particular damping configuration η = −0.4, ν = 0.33 is used in the both graphs in Fig. 6 . The left hand plot shows five curves which correspond the PDFs computed for increasing amplitude of the harmonic part, P = 0, . . . , 20, where the blue curve (P = 0) repeats the pure random case from the previous section, whereas the leftmost violet colour curve depicts attributes of the case with P = 20, where the deterministic part of excitation is dominant. The latter case is characteristic by a well localised maximum with a small dispersion. The same tendency was reported in [4] and directly follows from an analysis of the relation in Eq. (5).
The right hand plot in Fig. 6 describes the case where the amplitude of the harmonic part is kept constant (P = 1). The individual curves show the PDF for increasing amplitude (intensity) of the random part of excitation. When h = 0, the blue line ω = 1.6 shows the purely deterministic stationary amplitude. For increasing h the ordinate of the maximal probability increases, the same is also valid for the dispersion.

Conclusions
The approximative numerical stochastic analysis of a non-linear single-degree-of-freedom model was presented as a complement of previously published analytical formulas. A whitenoise and band-limited random excitation was applied to several parameter configurations. A direct comparison of numerically obtained histograms with analytically evaluated probability density functions was possible only for a purely random excitation; the comparison revealed in most cases an acceptable correspondence in both the shape of the PDF and in the ordinates of maximal probabilities. The numerical results were obtained using a simple simulation method. More advanced numerical methods, which are specific for the stochastic differential equations, will be used for real applications in the further work.
The kind support of Czech Science Foundation project No. GA19-21817S and of RVO 68378297 institutional support are gratefully acknowledged.