A constant phase approach for the frequency response of stochastic linear oscillators

. When studying a mechanical structure, evaluation of its frequency response function (FRF) over a given frequency range is one of the main interests. Computational cost aside, evaluating FRFs presents no methodological dif-ﬁculty in the deterministic case. Doing this when the model includes some uncertain parameters may however be more difﬁcult as multimodality can arise around resonances. Indeed, even for a single degree of freedom system, it can be shown that usual methods of the probabilistic frame such as generalized Polynomial Chaos may fail to properly describe the probability density function of the response amplitude. This study proposes another approach which involves a shift in the usual quantities used to draw FRFs. Instead of computing the stochastic response for a given excitation frequency, this work adopts a constant response phase point of view. For each phase value of the oscillator response, the uncertainty over some parameters is propagated to the corresponding uncertain amplitudes and excitation frequencies. This provides much smoother variations of the involved quantities which are much easier to describe using a simple Polynomial Chaos approach. This also provides a mean to conduct a straight-forward stochastic study of special points such as the maximum of amplitude (which matches a given phase) or the point for which the response is in quadrature with the excitation. Both analytical and numerical results will be exposed for a single degree of freedom oscillator whose squared eigen frequency (or stiffness) follows a uniform law.


INTRODUCTION
When studying a mechanical structure, evaluation of its frequency response function (FRF) over a given frequency range is one of the main interests.Computational cost aside, evaluating FRFs presents no methodological difficulty in the deterministic case.Doing this when the model includes some uncertain parameters may however be more difficult as multimodality can arise around resonances as demonstrated by Pagnacco et al. (2009Pagnacco et al. ( , 2011)).Indeed, even for a single degree of freedom system, it can be shown that usual methods of the probabilistic frame such as generalized Polynomial Chaos (gPC, (Xiu and Karniadakis, 2002)) may fail to properly describe the probability density function (PDF) of the response amplitude (Pagnacco et al., 2013).This latter work shows that more complex methods such as Multi-Element generalized Polynomial Chaos can be used to address this problem increasing the computational cost in return.This study proposes another approach to handle the frequency study of stochastic linear systems.It involves a shift in the usual quantities used to draw FRFs: instead of computing the stochastic response for a given excitation frequency, this work adopts a constant phase point of view.For each phase value of the oscillator response, the uncertainty over some parameters is propagated to the corresponding uncertain amplitudes and excitation frequencies.This work will be illustrated by a simple single degree of freedom (sdof) linear damped oscillator whose eigen frequency follows a uniform law.This system is describe in Sec. 2. Section 3 illustrates multimodality of system response amplitude when the response is sought for a given excitation frequency but variable (free) response phase.Section 4 then develops the proposed approach: the response is sought for a given phase but variable excitation frequency.Section 4.1 provides the equations while Sec.4.2 illustrates the approach on the sdof system.Finally, Sec. 5 proposes a comparison of both methods efficiency -constant excitation frequency and constant phase -when combined to a Polynomial Chaos Expansion.

STOCHASTIC SYSTEM STUDIED
Deterministic single degree of freedom oscillator Let us consider the sdof damped oscillator undergoing a harmonic load depicted in Fig. 1.Its movement is governed by Eq. (1).q + 2ηω 0 q + ω 2 0 q = f 0 cos(ωt) (1 where q is the mass displacement, q and q are its velocity and acceleration respectively, ω 0 is its eigen circular frequency, η is the damping ratio, f 0 is the excitation amplitude and ω is the excitation circular frequency.In the frequency domain, Eq. ( 1) becomes where j 2 = −1 and Q is the complex amplitude of q: Using this complex notations, one can easily write the complex solution Q as a function of the mechanical parameters ω 0 , η and f 0 and the excitation frequency ω: Finally, decomposing the complex quantity Q into its amplitude a and phase ϕ, one gets separate expressions for each component: Amplitude and phase response of the oscillator are displayed on Fig. 2 for f 0 = 1, ω 0 = 2π and η = 0.05.

Stochastic single degree of freedom oscillator
Let us now consider a probability space (Θ, A , P) with Θ the event space, A the σ -algebra on Θ, and P a probability measure.Random variables will be denoted by the capital letter which matches the deterministic variable.Hence if x is a deterministic variable, the associated random variable will be denoted X.Its cumulative distribution function (cdf) and probability density function (pdf) will be denoted P X (x) = P(X ≤ x) and p X (x) = dP X dx respectively.
We assume that ω 2 0 varies and can be modeled using a random variable Ω 2 0 (θ ) : Θ → R which follows a uniform distribution: This may happen when the oscillator stiffness has bounded variations.Numerical values for later numerical applications are: ω 2 0 = (2π) 2 , ∆ω 2 0 = 0.3ω 2 0 , η = 0.05 and f 0 = 1.Diagrams in Fig. 2 then match the mean sdof system response over the frequency range [0; 4π].Four operating points around which the stochastic response will be detailed are marked on Fig. 2 using letters (a) to (d).

PROBLEMS ARISING WHEN CONSIDERING A CONSTANT EXCITATION FREQUENCY
When studying a linear system over a given frequency range, it is natural to observe the variation of a for a set of excitation frequencies ω.That is what is usually done when considering deterministic structures: for several values of ω, amplitude a and phase ϕ are evaluated and plotted on graphics similar to Fig. 2. It then seems natural to use a similar procedure when studying stochastic linear structures over a given excitation frequency range: for several values of ω, a sample of Ω 2 0 realizations is generated and corresponding realizations of A and Φ are evaluated.The problem is that, for some values of ω, the probability density function of A is discontinuous as demonstrated by Pagnacco et al. (2011Pagnacco et al. ( , 2013)).This is illustrated on Fig. 3 and Fig. 4.These Monte Carlo simulations are obtained by considering 501 ω values equally distributed over the range [0; 4π].For each ω value, 20 001 realizations of Ω 2 0 are considered; these realizations are equally distributed over the range I (Ω 2 0 ).For each Ω 2 0 realization, the corresponding value for A is evaluated using Eq. ( 6) and stored.The pdf p A is then evaluated for each excitation frequency ω. Figure 3 displays p A for the whole ω range using colors while Fig. 4 displays p A in a classical way for the four ω values defined by operating points (a-d) marked in Fig. 2.
Panes (a) and (d) in Fig. 4 show smooth pdfs whereas discontinuous pdfs similar to pane (b) curve can be observed for the excitation range ω ∈ [1.66π; 2.28π].The exception in this range is the pdf obtained in pane (c).Detailed explanations for these behaviors can be found in the previously mentioned references.Only the main phenomenon will be outlined here using Fig. 5 which plots a versus ω 2 0 for the same four excitation frequencies values as in Fig. 4. Continuous pdfs of panes (a) and (d) in Fig. 4 are related to a bijective relation between ω 2 0 and a as displayed by panes (a) and (d) in Fig. 5. On the contrary, Fig. 4.(b) shows that different values of ω 2 0 map identical values of a. Being represented twice, theses values for a are linked to suddenly higher pdf values.The exception of Fig. 3.(c) comes from the symmetrical property visible in Fig. 5.(c): each a value has 2 preimages.Figure 6 displays plots similar to Fig. 5 except that they display ϕ instead of a when ω 2 0 varies.These discontinuous and possibly multimodal pdfs are difficult to obtain when applying the widely used Polynomial Chaos method to approximate A and Φ as illustrated further in Sec. 5.More complex methods must then be deployed to handle the problem (Pagnacco et al., 2013).

CONSIDERING A CONSTANT PHASE
To avoid the previously mentioned drawbacks of using a constant excitation frequency method which can be explained by the non-bijective link between the square eigen frequency ω 2 0 and the response amplitude a for a given ω value, let us observe a when the response phase ϕ is kept constant.Equation (4) creates a link between the triplet (a, ϕ, ω); instead of choosing ω and evaluating subsequent a and ϕ values, let us choose a ϕ value and evaluate subsequent a and ω values.From a mechanical point of view, this makes sense: when ω 2 0 varies, responses sharing a same phase ϕ match similar operating points (maximum response amplitude for example).From a mathematical point of view, it is interesting as, for a sdof oscillator, the link between ω 2 0 and a is bijective when ϕ is kept constant as illustrated by Fig. 9.This figure displays graphics equivalent to Fig. 5 (that is a versus ω 2 0 ) but for four given values of ϕ instead of using given values for ω. Figure 10 shows that the link between ω and ω 2 0 for these for given ϕ values is also bijective.The next subsection develops the equations giving the expressions of ω and a for a given ϕ value.An expression for A probability density function in the case when Ω 2 0 follows a uniform distribution is then derived, showing its continuity.The second subsection illustrates this constant phase method and provides graphics equivalent to those in Figs. 3 to 6.

Expressions of ω, a and p A for a given ϕ
For a given ϕ ∈] − π, 0[ value, Eq. ( 7) imposes ω as follows:    Case ϕ = − π 2 : which has two possible solutions: This leads to three cases for the displacement amplitude a formula due to Eq. ( 6): In every case, a can be rewritten where a ϕ is a coefficient depending on ϕ but not on ω 2 0 .Hence, formulas for A and its cdf can be derived: In the case when Ω 2 0 ֒→ U (ω 2 0 − ∆ω 2 0 ; ω 2 0 + ∆ω 2 0 ), one gets: which is smooth, unlike the case when ω is kept constant.

Numerical example
To illustrate the constant phase method, Monte Carlo simulations similar to those of Sec. 3 are carried out.501 φ values equally distributed over the range [−0.98π; −0.02π].For each φ value, 20 001 realizations of Ω 2 0 are considered; these realizations are equally distributed over the range I (Ω 2 0 ).For each Ω 2 0 realization, the corresponding values for Ω and A are evaluated using Eq. ( 12) and Eq. ( 14) or Eq. ( 9) and Eq. ( 15) or Eq. ( 13) and Eq. ( 16) depending on ϕ value.The pdf p A is then evaluated for each phase ϕ. Figure 7 displays p A for the whole ϕ range colors while Fig. 8 displays p A in a classical way for the four ϕ values defined by operating points (a-d) marked in Fig. 2.
All panes (a-d) in Fig. 8 show smooth pdfs while discontinuous pdfs where observed in the constant excitation frequency case, for ω ∈ [1.66π; 2.28π].This is justified mathematical by Eq. ( 21) which proves the smoothness of p A .It can also be understood by considering the bijective link between a and ω 2 0 emphasized by Fig. 9 which was previously mentioned.Finally let us point out that the link between ω and ω 2 0 is bijective too as illustrated by Fig. 10.To figure out the main differences between the spaces involved by each method (constant excitation frequency or constant phase), Fig. 11 displays the variations around operating points (a-d) when ω 2 0 varies over I (Ω 2 0 ) for both methods in classical diagrams: amplitude a and phase ϕ versus excitation frequency ω.

CONSEQUENCES ON A POLYNOMIAL CHAOS STUDY
Previous illustrations are based on Monte Carlo simulations: the direct problem (considering either a constant excitation frequency or a constant phase) is solved for a large sample of Ω 2 0 realizations.This can be afforded here because the system is very small and the cost of the direct evaluation is trivial.However, stochastic systems are often studied using an approximation of the stochastic response (Lucor and Karniadakis, 2004;Finette, 2006;Sarrouy et al., 2013) in order to decrease the computational cost.Among the different methods used to compute such approximations, the Polynomial Chaos expansion (PCE) introduced by Wiener (1938) and recently expanded to generalized Polynomial Chaos (gPC) expansion and Multi-Element generalized Polynomial Chaos (MEgPC) expansion (Xiu andKarniadakis, 2002, 2003;Wan and Karniadakis, 2005) is one of the most famous.
The next subsection provides a brief description of the simple Polynomial Chaos expansion.While the second one compares the results obtained when combining each method (constant excitation frequency and constant phase) with PCE.

Brief summary about Polynomial Chaos expansion
Only the principle is recalled here for a dimension-one stochastic space, that is when only one random variable ξ is used to introduce randomness in the system.The reader is referred to the references cited in above for a complete presentation of PCE.Considering a second-order random process X, the Polynomial Chaos expansion proposes to express it as a function X which is a polynomial series using a set of N orthogonal polynomials denoted ψ n in the variable ξ : where the order N is theoretically infinite for general situations.The deterministic coefficients x n are now used to represent X.They can be evaluated in two ways: using an intrusive method or a non-intrusive one.The intrusive method follows a Galerkin approach: Eq. ( 22) is introduced in the equations governing X and theses equations are projected onto the set of orthogonal polynomials ψ n .The non-intrusive method uses the orthogonality of the polynomials with respect to a scalar product denoted < •, • >: where the numerator is usually evaluated using a quadrature rule.
The main difference between both methods is that the intrusive method provides a set of m×N coupled algebraic equations (where m is the size of the underlying deterministic problem) and often requires a special implementation while the nonintrusive approach determines the set of coefficients x n one after the other in an independent manner and reuses existing codes to evaluate X realizations needed for the quadrature.The choice of the polynomial basis is somehow arbitrary even if some bases are considered as optimal to describe some distributions by some authors, as Xiu and Karniadakis (2002).In the present case, the random input Ω 2 0 follows a uniform distribution which makes the Legendre polynomial basis the most natural choice.The first 6 polynomials are: This set of polynomials is orthogonal with respect to the following scalar product ω/π (rad/s) ω/π (rad/s) ω/π (rad/s) ω/π (rad/s) Figure 11: Deterministic frequency response diagrams for the mean sdof damped oscillator plus variation of (A,C) amplitude and (B,D) phase when ω 2 0 varies around operating points (a-d) using a (A,B) constant frequency or a (C,D) constant phase approach.Grey patches display the of possible amplitudes or phases considering the whole ω 2 0 variation range.
The adequacy of the Legendre polynomial basis and the expansion on a random variable ξ that follows a uniform distribution U (−1; 1) and hence has p ξ (x) = 1 2 as probability density function may become visible if the numerator of Eq. ( 23) is rewritten as follows: where E[X] denotes the expected value of random variable X.
Once PCE coefficients x n are evaluated, there are two ways to post-process them.First, the mean and variance can be directly computed provided ψ 0 (x) = 1: Second, cdf and pdf can be evaluated based on MC simulations.The difference with the usual processing is that X realizations are computed using its PCE (i.e.Eq. ( 22)) rather than solving the direct problem which saves a lot of computational time and resource when the samples are large.

Application of PCE for both methods
Let us define the random variable ξ as follows: Its cdf P ξ and pdf p ξ can be easily established: It follows that ξ has a uniform distribution U (−1; +1).This random variable will serve to develop all the stochastic quantities around the 4 operating points (a-d) defined in Fig. 2 and which correspond to the 4 following triplets: When the constant excitation frequency method is applied, ω is set to ω op x (x ∈ {a, b, c, d}) and A and Φ expansions are evaluated.In the case when the constant phase method is applied, ϕ is set to ϕ op x (x ∈ {a, b, c, d}) and A and Ω expansions are evaluated.A degree 6 expansion is used in every case.The coefficients are evaluated using a non-intrusive method relying on a Gauss-Legendre quadrature with 7 nodes.
Figures 12 to 15 provide comparisons of Monte Carlo simulations and PCE results around operating points (a-d) when the constant excitation frequency method is used.As expected, this method provides correct result for operating points (a), Fig. 12 and (d), Fig. 15: variations of a and ϕ with ω 2 0 is well reproduced by the PCE and so are the respective pdfs.However, PCE does not provide a proper description of A and Φ for operating points (b), Fig. 13 and (c), Fig. 14.As stated in (Pagnacco et al., 2013), increasing the expansion degree would not provide better results.Figures 16 to 19 provide comparisons of Monte Carlo simulations and PCE results around operating points (a-d) when the constant phase method is used.In this case, the variation of a and ω with ω 2 0 is well described by the PCE as well as the corresponding pdf, both being much smoother than when the constant excitation frequency is used.

Comments on the practical use of the constant phase approach
The constant phase approach is obviously useful when one wants to study the variation of the system response in a particular configuration which is characterized by the phase ϕ: variation of the resonance peak (obtained for ϕ = arctan(− 1 − 2η 2 /η)) , variation of the system response when in quadrature with the excitation (ϕ = −π/2), ...However, one frequently wants to check that the system response will not exceed some given values over a range of excitation frequency.In this case, the constant excitation frequency approach seems more adapted but returns erroneous results for some system parameters when using PCE (see Pagnacco et al. (2011Pagnacco et al. ( , 2013))), especially around resonance that is where amplitudes are generally controlled.The constant phase approach can still be used if combined to a little post processing: let us denote [ f min ; f max ] the frequency range of interest.Evaluation of the response of the mean system for f min and f max provides coarse upper and lower bounds for the phase response.As depicted in Fig. 11, it is necessary to enlarge this range to properly cover the whole frequency range [ f min ; f max ].By enlarging this phase range, one gets the systems stochastic response over the desired frequency range.Depending on the desired statistical indicators (confidence interval, quantiles, moments, ...), an adapted post-processing can be implemented.This said, it is nonetheless interesting to keep in mind that being able to describe the variation of the resonance peak both in terms of amplitude and frequency is much more interesting that knowing that the amplitude may stay below a given value over f min and f max and ignoring it will explode for a little lower excitation frequency.

CONCLUSION
An original approach to study the dynamic response of a single degree of freedom system has been developed.This approach proposes to expand the system response on the Polynomial Chaos when imposing the response phase and freeing the system excitation frequency rather than the other way around as it is usually done.The proposed approach was applied to a single dof system whose squared eigen frequency follows a uniform law.This method was proven to provide much better results than the usual approach which suffers from the inability to describe the multimodality of the stochastic response.The numerical application also demonstrated its ability to follow some phase defined points such as the response of the system when in quadrature with the excitation which is usually close to the resonance point and is easy to detect experimentally.This work addressed a single dof system: further work should handle the case of multi-dofs systems.

REFERENCES
Finette, S., 2006, "A stochastic representation of environmental uncertainty and its coupling to acoustic wave propagation in ocean waveguides", Journal of the Acoustical Society of America, 120(5):2567-2579.

Figure 3 :Figure 4 :
Figure 3: Constant excitation frequency study: A pdf (p A ) over a given excitation frequency range.Color scale maps log(p A ). (a-d) cuts refer to operating points defined in Fig. 2.

Figure 7 :Figure 8 :
Figure 7: Constant phase study: A pdf over a given excitation frequency range.Color scale maps log(p A ). (a-d) cuts refer to operating points defined in Fig. 2.

Figure 12 :
Figure 12: Constant excitation frequency study: PC simulations compared to MC simulations for operating point (a).

Figure 13 :
Figure 13: Constant excitation frequency study: PC simulations compared to MC simulations for operating point (b).

Figure 14 :
Figure 14: Constant excitation frequency study: PC simulations compared to MC simulations for operating point (c).

Figure 15 :
Figure 15: Constant excitation frequency study: PC simulations compared to MC simulations for operating point (d).

Figure 16 :
Figure 16: Constant phase study: PC simulations compared to MC simulations for operating point (a).

Figure 17 :
Figure 17: Constant phase study: PC simulations compared to MC simulations for operating point (b).

Figure 18 :
Figure 18: Constant phase study: PC simulations compared to MC simulations for operating point (c).

Figure 19 :
Figure 19: Constant phase study: PC simulations compared to MC simulations for operating point (d).