Semi-analytical stochastic analysis of the generalized van der Pol system

The paper is motivated by a series of wind tunnel experiments, which deal with aeroelastic Single Degree of Freedom (SDOF) and Two Degrees of Freedom (TDOF) section models. Most of them can be mathematically expressed by van der Pol-Duffing type equations or their combination. Excitation due to aeroelastic forces consists mostly of a deterministic periodic part and random components, both of them are applied as additive processes. The lock-in state represents an auto-synchronization of the vortex shedding and basic eigen-frequency of the system. This problem seems to be very polymorphous and, therefore, several isolated regimes have been outlined together with their characterization. Parameter setting with solely random excitation is further investigated in the paper. The strategy of stochastic averaging is then employed to formulate normal form of stochastic system for partial amplitudes of harmonic approximates of the response. The random part of excitation is considered as a Gaussian process with significantly variable spectral density. Hence, a conventional way of investigation based on an idea of white noise excitation is no more applicable. Therefore, the general formulation of diffuse and drift coefficients should be used to construct the relevant Fokker-Planck equation (FPE). Semi-analytical solution of FPE is deduced in the exponential form by means of a probability potential. It is later used for stochastic stability investigation together with consideration about the stationary probability distribution existence. Open problems and further research steps are outlined. c © 2018 University of West Bohemia. All rights reserved.


Introduction
The paper is motivated by a series of wind tunnel experiments investigating aeroelastic Single Degree of Freedom (SDOF) and Two Degrees of Freedom (TDOF) section models of various shape and aeroelastic properties.It reveals that most of them can be theoretically modeled by van der Pol-Duffing type equations or their combination adjusting degree of individual nonlinear terms or their coefficients.Many experimental analyzes have been published proving this hypothesis, see for instance papers [1,3].It should be emphasized that this character of the system response is very stable and is obvious in linear as well as in weakly nonlinear domain when post-critical effect emerge.Moreover, many special effects identified by experimental way evoke properties recognized in the pure theory of differential equations.Authors take part in this research as well, see for instance [8] or [11], following both theoretical and experimental way of investigation.
In general, the commonly adopted model is of TDOF type, see Fig. 1.It seems to be meaningful from theoretical point of view to investigate individual variants of linear and nonlinear equations/systems representing large variety of cases which can be encountered.Nevertheless, the universal general TDOF model is complicated and does not enable to get through many mechanisms and phenomena regarding stochastic stability on the level of an (semi) analytical approach and only a numerical simulation of particular settings is possible.However, many particular settings exhibit one dominant component only and, therefore, the second component can be more or less neglected.For precise analysis of individual alternatives with respect to bifurcation in a critical and early post-critical state, see [6,11].In such a case one speaks about the simplified flutter, whatever is the dominating component (pitching or heaving).Both alternatives have been observed in a wind tunnel.In theoretical treatment can be the model formulated as an SDOF system.
Earlier papers, e.g., [12] and others, show a variability of the ratio ϕ(t)/u(t), provided the system finds on the aeroelastic stability limit and the Hopf bifurcation decides about the subsequent post-critical process.In dependence from cross-section shape, elements of the linear part of the stiffness matrix and damping parameters, both components of the response can be commeasurable.Nevertheless, it seems with respect to results of the careful parametric analysis a dominance of one component, either ϕ(t) or u(t), is more probable.For this reason it is worthy to start with an SDOF version in order to verify applicability and effectiveness of mathematical procedures, which are based on Markov processes approach and Fokker-Planck equation analysis.Moreover the SDOF formulation leads to more transparent results and enable to explain explicitly special effects which occur during the process of the stability loss.
Excitation due to aeroelastic forces is in principle a combination consisting of a deterministic periodic part (vortex shedding) and of random components represented by stochastic processes of additive and multiplicative types.Which one is dominant depends on the stream velocity V and, consequently, on aeroelastic processes ruling around the investigated profile.Basically three regimes can be encountered: (i) V 0 < V < V crit : vortex shedding is dominant; despite this rather deterministic process is accompanied by random excitation processes, their influence is very small and excitation and response can be taken to be deterministic with the known frequency, V 0 is a bottom limit of the vortex shedding regime.(ii) V = V crit ± δ, where δ is a small velocity detuning from V crit : influence of the vortex shedding and random perturbations are co-measurable and quasi-periodic effects represent the response.(iii) V > V crit : vortex shedding nearly disappeared and random effects are dominating.
For this reason a series of consecutive papers dealing with various types of special and more general mathematical models considered as a theoretical support of aeroelastic experiments in the wind tunnel are planed.As the first step authors attempt to develop an authentic theoretical counterpart characterizing response of a nonlinear SDOF system associated with aeroelastic model, which is investigated in conditions of the homogeneous stream with the velocity slowly sweeping up and down outside the resonance domain.In this regime the random excitation component of the additive type is dominant.This regime corresponds to super-critical state, see the variant (iii) above, typical for high Reynolds number being characteristic for high velocities of the stream.It is completely different from the lock-in state typical for the variant (ii) ruling in the resonance zone where the vortex shedding with a small detuning δ from the system eigen-frequency and relevant random component makes a combination of deterministic and stochastic excitation producing a number of special effects.Let us take a note, that following experimental results the exciting random processes are different from a white noise and its spectral density cannot be considered as a constant.A sample of the relevant spectral density measured in a wind channel is shown in Fig. 2, where we can see a significant variability of the spectral density in the frequency domain.
In the same time we can regard these processes to be ergodic and correlation stationary.
So that, we concentrate in this paper to the case (iii).In the first step the basic mathematical model will be formulated.Later stationary response using stochastic averaging will be investigated together with some parametric analyzes.

Character of basic mathematical model
Let us consider the van der Pol with extended nonlinear damping part with random excitation of the additive type: where: −4 ]; basically it holds ν > 0, while η should be considered positive or negative, and ϑ varies in an interval (−∞, ϑ max ), where ω 0 -eigen-frequency of the adjoint linear SDOF system, ω -excitation frequency of the vortex shedding [s −1 ], P ω 2 -amplitude of the harmonic excitation force [m•s −2 ], ξ(t) -broadband weakly stationary Gaussian random process When applied in aeroelasticity of the SDOF or TDOF systems, we should be aware that coefficients η, ν, ϑ are functions of the stream velocity V .Regarding ν, ϑ, they vanish for V = 0, as they follow from aeroelastic processes only.Coefficient η consists of two components: (i) elastic part η e < 0 corresponding to conventional damping ratio of a linear SDOF system and (ii) aeroelastic part η a ≥ 0, which is zero for V = 0 and rises monotonously for increasing V .Therefore, the total value of η crosses zero for a certain value of V z and becomes positive.Starting this point, the trivial solution of the homogeneous equation is no more stable and nonlinear part of the damping should stabilize the response in the first limit cycle which is stable and makes an attractor, if ϑ < ϑ max .Consequently, coefficients η, ν, ϑ and their relation are responsible for the response portrait and solution stability, for more details, see, e.g., [13,17].The ratio of these parameters decides about existence of respective Limit Cycles (LC).It can be shown that in aeroelasticity of systems modeled by (1), see, e.g., [9], can exist one (stable) or two (one stable and one unstable) limit cycles, see Fig. 3.If η < 0, the homogeneous system is stable in the trivial state and no LC emerge.Therefore, we concentrate to case for η > 0, which indicates that linear part of the damping becomes negative and only its nonlinear part could stabilize the response on the level of the stable LC, if any.The particular shape of the stabilized state depends on ϑ, see Fig. 3: (i) −∞ < ϑ < 0 -only one stable LC exists, the system cannot be destabilized and moves within the attractor area around this LC; (ii) ϑ = 0 -one stable LC, although the second unstable LC also exists, the infinite volume of energy would be needed to reach it; (iii) 0 < ϑ < ϑ maxone stable (inner) and one unstable (outer) LC exists.The reaching the outer unstable LC means final orbital stability loss.Beyond this limit the system response is no more periodic and its amplitude rises to infinity; (iv) ϑ = ϑ max -both LC coincides into one double LC which is stable on internal side and unstable (repulsive) on outer side; (v) ϑ > ϑ max -the system is absolutely unstable and absolute value of its response is monotonously rising beyond all limits.

Deterministic and random parts of the excitation
Concerning the right side of (1), the vortex shedding is the origin of additive harmonic excitation.First of all we look through the vortex shedding frequency as a function of the stream velocity in the wind tunnel.Observing Fig. 4, the lock-in zone is obvious, for details see [1].From the dynamical viewpoint the lock-in zone represents the resonance zone where the vortex shedding frequency is nearly synchronized with the eigen-frequency of the system ω 0 even if a certain detuning is visible.The detuning can be influenced by means of conditions in which the system is working, e.g., roughness of the surface, damping ratio, etc.Therefore, the mechanism in a particular case should be carefully studied in order to adopt an appropriate measure for suppression of this undesired phenomenon.Fig. 4. Vortex shedding frequency as a function of the stream velocity, see [1] These considerations lead to the finding that in the resonance zone two components of the total response are to be expected: (i) self-excited vibration following from the fact that the trivial solution of the homogeneous variant of (1) already has fallen into unstable state and (ii) forced vibration due to deterministic part harmonic excitation on the right side.The resulting response exhibits a character of quasi-periodic beatings.This phenomenon can be found in theoretical studies of the van der Pol systems, see, e.g., [8], but it has been proved by experimental measurements as well, e.g., [3].
Outside the resonance zone the vortex shedding frequency is given by aerodynamic properties of the model and the stream velocity.In such a frequency domain no self-exciting response component can emerge and hence the response has a character of solely forced vibration, as it has been pointed out in Section 1.In this zone random excitation components are dominant.This particular setting will be treated in this paper.

Semi-analytical solution form -Ito system
It is understood that the terms (η − νu 2 + ϑu 4 ) • u and P ω 2 are of a small order ε, and h • ξ(t) is of order ε 1/2 , without being so indicated in (1).Respecting these conditions we are entitled to apply solution methods supposing a slow variability of the response partial amplitudes following this expression: with auxiliary condition: where argument of the "slow" time is omitted: Substituting assumptions (3) into the system (1) and employing the condition (4), one obtains after several modifications the normal form of the differential system for amplitudes a c , a s , which are slowly variable in time: When detuning 0 < δ < δ r is small, so 2δ = (ω 2 0 − ω 2 )/ω ≈ |ω 0 − ω| is small as well.The same holds regarding right sides of equation ( 5) and, therefore, a c , a s are slowly variable amplitudes.The input process can be considered of the Markov type and, therefore, we are entitled to rewrite the problem in a form of Ito stochastic differential system and moreover the stochastic averaging operation can be applied.This operation eliminates the fast varying parts in (5).As the result we obtain simplified Ito differential system having the general form, see, e.g., [2,5,14]: where D c , D s are averaged deterministic parts of the right-hand sides of (5).Because random processes do not act in parametric position, respective Wong-Zakai parts vanish, see, e.g., [7,16].
As the deterministic parts of the right-hand sides of ( 5) are smooth enough, D c , D s equal to their averaged values without any complement.B w (t) is a unit Wiener processes and σ cc , σ ss are expressions, which can be written in the form, see, e.g., [2,5,14]: where R(τ ) is the auto-correlation function of the process ξ(t).If ξ(t) degenerates to a Gaussian white noise, R(τ ) = dir(τ ) and σ 2 ii (t) = 2πK ξξ g 2 ii (t), where K ξξ is the stochastic intensity of the process ξ(t).Coefficients g ii are coefficients in random parts of the right-hand sides of (5).In particular: If the process ξ(t) is not the Gaussian white noise, as it corresponds with experiments, then the integral (7) should be evaluated with respect to (8): where Φ ξξ (ω) is the spectral density of the process ξ(t) at frequency ω.Subjecting now the system of equations ( 5) to operation of averaging with respect to above considerations, one obtains after a couple of adaptations the Ito system as follows: where B c (t) is the Wiener process related with input excitation ξ(t).

Fokker-Planck equation for the response partial amplitudes
Equations (10) indicate that, given time, the unknown a c , a s will tend to stationary random processes if the averaged system is stable.Then the reduced FPE for the stationary cross probability density function (PDF) p(a c , a s ) (left side of the FPE is put to zero) can be written in the form: FPE (11) should be solved subjected to the boundary conditions: We should be aware that even when the averaged random processes a c , a s reach the stationary state, the original processes a c , a s may remain non-stationary.
Let us add to the FPE a neutral member: D∂ 2 p/∂a c ∂a s − D∂ 2 p/∂a c ∂a s .Putting that into (11) one obtains after a mild modification: where D is an arbitrary constant.It can be shown, see [15], that equation ( 13) can be satisfied, if it holds separately: The above equations can be considered as the system of two linear algebraic equations for ∂p/∂a c and ∂p/∂a s .Hence it holds: where it has been denoted: S = Φ ξξ (ω)/2ω.The solution of ( 11) can be expressed by means of stationary potential: where C is the normalization constant.The potential Φ(a c , a s ) exists provided that the following condition is fulfilled: Upon substituting ( 15) into ( 17) we obtain: Let us fix the free parameter D to eliminate the first part of the left side of equation ( 18): Anyway, for arbitrary a c , a s , equation ( 18) can be satisfied only if it holds: D = 0 which implicates (i) S = 0 or absence of random excitation, or (ii) δ = 0 or coincidence ω = ω 0 and vanishing amplitude of harmonic excitation part.The latter case means that random excitation only is acting.Let us have a look at the case of purely random excitation.As it has been stressed earlier, the excitation process is in general non-Gaussian and, therefore, the spectral density Φ(ω) is not constant.Nevertheless, we still consider the excitation for a broadband process.Its influence in FPE discussed (13) manifests through the fixed value S = Φ(ω 0 )/ω 0 , that is as the spectral density value in the basic eigen-frequency of the system ω 0 .
Returning back to equations ( 15), they can be rewritten with respect to ( 17), ( 18) and ( 19) in the form: Each of equations ( 20) can be integrated following variable a c or a s and to result an arbitrary function of a s or a c can be added, respectively.Finally, we obtain the same result of both integrations providing the stationary potential in the form: It is obvious that the potential function is rotation symmetric.Therefore, the PDF following equation ( 16) depends on the absolute amplitude A 2 = a 2 c + a 2 s only and remains independent on the phase shift ϕ = a c /a s .In other words, the probability distribution in circumferential direction is uniform.The transformation into the polar coordinates and subsequent integration of ( 16) following circumferential direction leads to PDF in the radial direction A ≥ 0: We shall examine the above formula.Its basic character follows from the relation of damping parameters η, ν, ϑ.Let us determine the zero points of the potential Φ: where A 2,4,6 ≤ 0 were omitted as A ≥ 0 only is considered.As amplitude A 2 is a real value, A 3,5 following (23) is meaningful only if holds: where the ν 2 = 32 3 ηϑ means a certain limit case of the double root A 2 3 = A 2 5 , beyond which real roots A 3,5 do not exist.Positions of potential extreme values are: provided it is valid: Several plots of the potential Φ(A) and corresponding PDF p(A) are presented in Figs. 5 and 6 in Section 3.3.

System parameter ϑ
Let us discuss two cases, in particular when ϑ is negative and when ϑ is positive.Provided the parameter ϑ ≤ 0, conditions (24) and ( 26) are satisfied implicitly.The term with A 4 has a stabilizing character, unstable limit cycle vanishes and all trajectories are approaching the inner limit cycle which is stable.In such a case the function ( 22) is integrable and the normalizing constant C can be evaluated, even if it cannot be evaluated exactly and should be evaluated numerically.
We pay attention to some plots exhibited in Fig. 5.They represent potential Φ(a c , a s ) of the PDF with respect to equation ( 21) for settings which are indicated in individual pictures and caption.The bold black curve in every picture means the case for ϑ = 0.It separates cases providing stationary cases (ϑ ≤ 0 and cases leading to collapse of the system ϑ > 0. It is obvious that curves for ϑ ≤ 0 are exclusively in upper half-plane starting a certain zero crossing point and getting the negative potential Φ(a c , a s ) for every higher values A. It enables integrability of Φ(a c , a s ) on the whole space and guarantees existence of the normalization constant C. Unlike from that ϑ positive leads every time to positive potential Φ(a c , a s ) since some starting point and, therefore, to exploding increase of the PDF.Roughly speaking, it is apparent that curves for negative ϑ starting a certain value of A always rise to +∞ as A 6 .Therefore, PDF following (22) rises beyond all limits and normalization constant C does not exist.Provided ϑ is positive, the potential Φ(a c , a s ) becomes negative since the last zero point.Consequently, the exponential in (22) approaches asymptotically zero for A → ∞ and is integrable.Hence, the stationary PDF exists and can be physically interpreted.
Take a notice that ν having the stabilizing character (being always positive) has a diminishing influence on the PDF.The potential takes up even more and more smaller area of negative values and leads to more shallow PDF.Rising ν has stabilizing character and reduces the response variance of the system.We can observe this process also in Fig. 6, where the system response PDF is exhibited for several corresponding parameter settings as in Fig. 5.
As a verification of results obtained in this paper using a semi-analytical procedure can serve the analysis performed numerically by means of the Finite Element Method (FEM), see [4] or [10] and other resources.Although qualitative match is obvious some differences in details can be observed.They follow from a limited extent of trigonometric series (3).They imply slightly more complicated form of the limit cycle and, therefore, not constant crest of the PDF which we can see in Fig. 6.
Let us pay attention now to case with the positive value of the parameter ϑ.This setting changes character of further investigation substantially.Provided the excitation is deterministic (e.g., harmonic) the outer limit cycle represents a limit beyond which the response rises to infinity and the solution as the system response still enable a meaningful interpretation from the viewpoint of the stability analysis.Unlike the deterministic case an existence of the outer unstable (repulsive) limit cycle leads in principle to state (even if with a very small probability) when the stochastic stability limit is reached in final time and is broken through.The stationary or asymptotically stable solution of FPE does not exist any more.It means that the operation of the stochastic averaging which leads to the system (10) is no more applicable.From the viewpoint of the function (22) the PDF is not integrable with respect to amplitude A and the normalization constant does not exist.Therefore, if the excitation process admits infinite values, the solution cannot be considered as a PDF either constant (independent on time) or variable with respect to the slow time having a form of a function either asymptotically converging to a certain integrable limit or to a quasi-periodic function.It is obvious that the PDF of excitation noise ξ(t) has a crucial character.Basically, the excitation PDF could be defined with limited values preventing that the response never reach the outer limit cycle.However, such a formulation is beyond the definition of basic stochastic system leading to the FPE (11) supposing Gaussian character of the excitation.
This facts can be demonstrated by pictures plotted in the right column of Fig. 6.It is obvious that PDF in the area nearby the origin behaves similarly for both negative and positive ϑ if the absolute value of this parameter is identical (let us remember again that PDF is plotted in non-normalized state).Nevertheless since a certain distance from the origin, PDF starts increase dramatically, following the exponential law.
These findings force us to go back to reconsider results obtained in wind tunnels.It is known that model tested in a wind tunnel can be brought to collapse, which corresponds with the model investigated in this paper.Therefore, the task becomes a character of the problem of the first excursion beyond a defined limit.On the other hand the time of the first excursion beyond a critical limit is very long and in practice is hardly to be reached.Another possibility would still be to abandon the perfectly Gaussian excitation character and admit the PDF of the exciting process to be given by a so-called limited distribution of probability.In other words, the PDF of excitation is, for instance, rectangular with limitations before the position of the second (unstable) limit cycle.Of course, this approach would require a completely different solution procedure than that used in this paper.Authors will pay attention to these approaches in the next paper.

Conclusions
The problem of a mathematical model characterizing the generalized van der Pol system under combined harmonic and random excitations has been discussed and classified into a couple of categories.Namely TDOF and SDOF systems have been mentioned as an integral framework.At the moment the SDOF system has been treated in general and in details regarding the response probability density.This basis will be used from a methodological point of view as a background of an extensive study of the general TDOF system with combined deterministic and random excitation.This strategy has been inspired by results of TDOF stability study, which indicates possible dominance of single response component (heaving or pitching) in the moment of stability loss in the Hopf bifurcation point.Regarding the random part of excitation a distinctly variable spectral density in frequency has been considered.Commonly adopted suggestion of its white noise character has been abandoned due to experimental results in a wind tunnel measurements.While the deterministic task admits to work with both stable and unstable limit cycles, the random broadband excitation would lead to unstable collapsing solution.Therefore, the cases with an unstable limit cycle have been avoided and they will be studied separately as a problem of the first excursion of the reliability domain.A particular solution has been obtained using the stochastic averaging, which provides the solution of exponential PDF distribution with the stationary homogeneous probability potential.A short comparison with numerical solution of the relevant Fokker-Planck equation has been done.Applicability and shortcomings of approaches used are commented.A few hints for engineering applications in a design practice are given.Open problems are indicated together with enumeration of variants planed to be investigated in the near future.

Fig. 2 .
Fig. 2. Spectral density of the in wind additive excitation due to pressure fluctuations