Hopf bifurcation for a fractional van der Pol oscillator and applications to aerodynamics: implications in flutter

Flutter is an important instability in aeroelasticity. In this work,we derive a model for this phenomenon which naturally leads to an equation similar to a van der Pol oscillator in which the friction term is given by a fractional derivative. Motivated by these considerations,we study a fractional van der Pol oscillator and show that it exhibits a Hopf bifurcation. The model is based on a one-dimensional reduction where the instabilities associated with flutter are preserved. However, due to the fractional derivative, the bifurcation analysis differs from the standard case. We present both analytical and numerical results and discuss the implications to aerodynamics. Additionally, we contrast our qualitative results with experimental data.


Introduction
Mathematical modeling of memory effects requires non-local in time constitutive equations. Fractional derivatives are thus well suited to capture memory effects in viscoelastic materials, as noted in [1]. Fluid-structure interaction is one scenario where a system acquires effective viscoelastic properties, even though the materials themselves are not viscoelastic, as explained in [2]. A well-established problem in engineering involving mechanical instabilities is that of the aircraft wing in flight. The general theory of aerodynamic instability and the mechanism of flutter is attributed to Theodorsen [3]. In the 1935 technical report, Theodorsen relies on potential flow to model both the effects of translation and rotation of the wing. This line of research focused on studying the behavior of classical flutter for complex aircraft configurations, including wing-fuel tank interactions by Sewall [4] and Reese [5], wing-flap interactions, among others.
There are two kinds of flutter, classical and stall flutter. While the former emerges as flow velocity exceeds a critical value, the latter ensues as a result of leading-edge vortex shedding. Dynamic stall was extensively studied by McCroskey in [6,7]. A first approach to model it using point vortex dynamics was used by Ham in [8,9]. However, these attempts do not account for the fluid-structure interactions involved in flutter, a process which introduces memory effects [6]. More recent approaches to the problem have predicted classical flutter with reasonable accuracy, but stall flutter remains elusive [10]. Hence, taking advantage of a fractional derivative-based model is appropriate since it captures memory effects.

A fractional model for flutter
From experimental results, it is known that classical flutter can only materialize in a system that effectively behaves as if it had at least two degrees of freedom. On the other hand, it is known that this phenomenon can take place in one degree-of-freedom systems if it originates from stall. The degree of freedom corresponds to the vertical bending of the cross section of the wing, as it can be observed in the diagram of the mechanical system proposed (Fig. 1). Therefore, it is desirable to understand how to obtain a simplified mathematical model for these instabilities from the conservation laws of solid and fluid mechanics.
We choose a fractional derivative model to describe these instabilities due to the limitations of classical models. This deduction is not straightforward and frequently little emphasis is put on the details of the procedure. Consequently, we believe it is important to review the relevant steps in order to accurately interpret the physical meaning of the fractional derivative in this context. In order to simplify the argument, we make the following assumptions: (1) the airfoil is thin enough to be approximated by a flat plate of length L and width w; and (2) the motion of the plate is exclusively a vertical translation. That is, the position of the plate is given by where X is the equilibrium position. By integrating the conservation law of momentum for the solid in an Eulerian reference frame, i.e., whereü is the second derivative with respect to time, ρ is the density and σ is the stress tensor. Notice that the density can be interpreted as a piecewise continuous function across ∂Ω. Considering the aircraft wing as a homogeneous linearly elastic body, we arrive atü where ω 2 is the normalized stiffness coefficient, η is the normalized dissipation coefficient and f is the force per unit mass that the fluid is exerting on the structure, which in particular depends on the parameter α. The latter corresponds to the fractional derivative order, which is associated with the degree of viscoelasticity. Thus, f contains both the information of the fluid and the interaction between the two media. Of course, Eq.
(3) can be transformed using Eq. (1) into its familiar form in a Lagrangian reference frame. Next, we proceed to deduce the functional form of f . We manage to do this by integrating the conservation law of momentum for the fluid in an Eulerian reference frame, i.e., and by incorporating the incompressibility condition. In the equations above, σ is the stress tensor which is decomposed in the familiar spherical and deviatoric parts, σ = τ − pI. We use the following fractional, non-local in time, constitutive equation where κ is the normalized volumetric dissipative coefficient, included here just for the sake of completeness, and C 0 D α t is the Caputo derivative where 0 < α ≤ 1, which is consistent with a homogeneous linearly viscoelastic material. On the other hand, κ 0 , η 0 , T α−1 λ and T α−1 μ are factors which are introduced to preserve dimensional homogeneity. Observe that when α = 1, the constitutive relation for a Newtonian fluid is recovered. Equation (5) is invariant under translations (Galilean invariance) and rotations. This is a particular case of the constitutive equations proposed by Rivlin in [11][12][13]. For more information regarding this case, see [14]. The result of the aforementioned procedure leads tö where J is the Jacobian of the transformation. Due to Newton's third law, we know Eq. (7) must be equal to f . Hence, the simplified model of this mechanical system in its most general form isẍ Equation (10) is crucial because integrated forces I 1 and I 2 have different physical origins and, as a matter of fact, only I 2 contemplates the interaction between the flow and the structure. Depending on the functional forms of I 1 and I 2 , bifurcations concomitant with coefficients of the model will be associated with divergence, flutter or a combination of both. More specifically, I 1 is linked to external interactions, while I 2 will determine the abruptness of the self-excitation onset [15].
Equations (5) and (9) are novel, in particular the latter relates the divergence-flutter effects.
which is the standard van der Pol equation. Since it was proposed and studied [16][17][18], it has been used to model robust oscillations in many fields. Flutter on aerohydrofoils is no exception [15,19].
Equation (11) has an equilibrium point at (0, 0) and it is a well-known fact that it exhibits a Hopf bifurcation for μ = 0. Based on the previous section, we study the following fractional van der Pol equation where ω is the frequency, μ is a dissipation parameter and α is the fractional derivative order. When α = 1, we recover the classical van der Pol equation. The idea of this proposal is to substitute the termẋ representing the usual friction or dissipation as a fractional derivative due to viscoelastic effects. We choose the Caputo derivative as previously discussed for the following procedure. Let's define We write Eq. (12) as a first-order system, allowing to set the subsequent two variants: or Both systems (13)-(15) and (13), (16), (17) have a stationary point at X 0 = (0, 0, 0) and possess the same system matrix, which is where ζ =ẑ,z. Their corresponding linearizations around X 0 are , respectively, where It must be emphasized that ω and μ on Eq. (21) become bifurcation parameters.

Stability theorems and results
We will analyze the stability of equation (12) based on the stability of the systems (19) and (20). Moreover, by virtue of the next theorems [20,21], they are the same and can be studied directly from the matrix A.
are rational numbers such that 0 < α i ≤ 1. Let M be the lowest common multiple (LCM) of the denominators of the α i 's. Then for the system  Theorem 2 If we assume that the conditions of Theorem 1 hold except replacing the Caputo derivative and the initial value X 0 = X (0) by the Riemann-Liouville derivative and the initial value X 0 = RL D α t 0 ,t X (t) t=0 , respectively, then the stability result is still available.
If α 1 , . . . , α n are not rational numbers but real numbers between 0 and 1, then we have the following result, which was introduced using the initial-value theorem of Laplace's transform.

Theorem 3 If all the roots of the characteristic equation
det diag s α 1 , s α 2 , . . . , s α n − A = 0 have negative real parts, then the zero solution of the system is asymptotically stable, where α i is real and lies in (0, 1).
Our interest is to study the stability for the systems (19) and (20). Suppose α = p/q, then the denominators of {1, α, 1 − α} are 1, q and a divisor of q, respectively, thus in the theorem M = q. From Theorem 1,we rely on the zeros of the polynomial det diag λ qα 1 , λ qα 2 , . . . , λ qα n − A , for A given by (21) and α =α orᾱ. For both systems (19) and (20), the resulting polynomial equation is or, for Λ = λ q , we have It is worth noticing that there are 2q complex roots for 0 < α ≤ 1.
To find these roots,we can always consider ω ∈ {0, 1}. Otherwise, the scalinĝ μ = μω p/q−2 is used to reduce it to this case. In this manner, from the polynomial p(τ ) ≡ τ 2q −μτ p + 1, and functionĝ (T ) ≡ T 2 −μT α + 1, the zeros of (22) and (23) are obtained, respectively, with the change of variables Moreover, this change of variables does not affect the argument angle of the original roots. Example For α = 0.5, the characteristic polynomial is of 4th order: that can be solved by radicals in the form: where s 1 , s 2 ∈ {−1, 1}. To continue with the stability analysis, we study the real part of Λ = λ q , e.g.,| arg(Λ)| = |q arg(λ)|. Rising (24) to the q = 2 power, we obtain For μ > 0, the values of Λ are given by: s 1 = −1: • if R > 0 then 2μ > A 3 and so Λ 3 and Λ 4 are real numbers with the same sign as: and, in particular, for s 2 = −1, Λ 4 is positive.

Hopf bifurcation
A Hopf bifurcation occurs at a parameter value where a periodic solution arises as the stability of an equilibrium point changes, i.e., the appearance or the disappearance of a periodic orbit for the dynamical system [22]. For Eq. (12), it means that a pair of complex conjugate zeros of the polynomial (22) crosses the critical argument angle Let λ 0 be a root with |q arg(λ 0 )| = π/2 and ω = 0, it implies and as the gcd( p, q) = 1 then λ p 0 cannot be a real number, this leads to the only possibility of μ = 0 and |λ 0 | 2q = ω 2 , so arg(λ 0 ) is not affected by the choice of ω 2 . Conversely for μ = 0, is a root of (12) satisfying |q arg(λ 0 )| = π/2. That is, the only possibility for the existence of a stable critical point at 0 is in the case of μ = 0.
We should notice also that in Eqs. (22)-(23) So, from Bolzano's (intermediate-value) Theorem, for a fixed μ ≥ 1+ω 2 , there would be at least one real positive root, i.e., | arg(Λ 0 )| = | arg(λ 0 )| = 0. From the previous analysis and using the implicit function theorem, the arguments of the root are continuous functions of μ. It is clear that the inequality holds for all μ > 0. Therefore solutions are stable, being μ = 0 the bifurcation point (see Fig. 2).

Numerical solutions and phase planes
The (linearized) fractional differential equations (19)- (20) were solved using the program FDE_PI1_Im for MATLAB created by Garrappa and published in [23]. The latter relies on an implicit method that solves the equation system as shown in (13), (14) and (15), or equivalently (13), (16) and (17). The implicit method requires introducing the Jacobian of the right-hand side of the equations. The solution is then achieved via Newton-Raphson iterations. The following results were obtained using a 14.9467 • initial value and considering the system begins at rest for all cases. Figure 3 shows multiple solutions as functions of time for different values of α. It can be observed that out of the solutions shown in the Figure, the one corresponding to the value α = 0.1 approximates best the case when α is an integer. Notice that as α → 0 the solution will tend to a harmonic oscillator. As α grows, a major damping effect occurs. It is also worth noticing a change of phase between the different solutions.
The phase diagram of the solutions is presented in Fig. 4. The damping effect is observed as well for the different values of α. It can be observed that even when the value is 0.9, the solution tends to zero. In this graph, the change of phase between solutions can be better appreciated.
Given the fact that a fractional differential equation is solved, a corresponding phase diagram for each value of α was obtained. In this case, the fractional derivative is plotted against the state variable, shown in Fig. 5. It is in fact a modified phase  portrait which can be called a fractional phase plane. Since μ and ω are fixed, the present case is in fact a generalization of the classical phase diagram. This implies the angle between this fractional phase diagram and the classical phase diagram is directly related to the degree of viscoelasticity of the system under study. Furthermore, there is a phase difference between the responses for distinct values of α. A potential interesting physical interpretation for this fractional phase plane may be proposed. The nonlinear fractional system given by equations (13), (16), (17) was solved using the program MT_FDE_PI1_Im for MATLAB created by Garrappa and published in [23]. Figure 6 presents areas formed by initial conditions pairs composed by position and velocity that give stable soluitions, i.e., solutions that eventually tend to zero. Each value of α is an area of different color, bounded by a solid curve for better comprehension. Figure 6 is made through the following procedure. First, an array of approximately 270-by-270 points was defined with a given set of initial conditions for the nonlinear equation system given by (13)- (15). By solving the aforementioned equation system   and only retaining those points which led to a bounded evolution in time, the basins of attraction for different values of α were obtained. Several features of the basins of attraction can be observed. First of all, the basin is rotated clockwise as α increases. At the same time, the area covered by the stable region becomes larger. Two sharp edges can be identified on the limits of the basin corresponding to α = 0.1. By introducing (x 0 ,ẋ 0 ) pairs of values in the vicinity of the sharp edges and increasing the computing time of the solution, it is apparent the sharp edges are not an artifact of the numerical method. An analysis of the qualitative behavior of the solution near these edges is suitable. Figure 7 highlights a single solution spiraling toward zero. This implies that the orbit corresponding to a set of initial conditions which result in a bounded evolution of the system can nonetheless exit the basin of attraction at later times. This is to be contrasted with what happens in standard dynamical systems and it is due to the non-local effects of the fractional derivative, i.e., memory. Self-intersections of the orbit can be observed near the attractor owing, once again, to the memory effects that result from the fractional derivative.

Conclusions
The fractional model for flutter proposed in Eq. (12) presents a more realistic approach to study the phenomenon, allowing for the consideration of memory effects due to the effective viscoelastic properties of the system. It has been pointed out in [6] that these systems present history effects (memory). These materialize as hysteresis loops in the aerodynamic forces and moments of the wing. Furthermore, it is emphasized in [6] that these "are not clearly correlatable based on any simple parameter of the mechanics of oscillation." The present approach may provide a solution by introducing an additional quantity, i.e., the fractional viscoelastic parameter α. This parameter may be a physical property of the system. The latter changes the phase of the response even if the other parameters remain fixed and leads to generalization of the classical phase plane. Even when Eqs. (11) and (12) exhibit a Hopf bifurcation, the way in which it happens is different. In particular, the amplitude of the oscillations is, for the fractional case, larger than those of the standard model. The value of α affects the abruptness of the onset of the bifurcation while leaving the critical value of μ for which the onset takes place unchanged. This has important implications in real airfoils, given the fact that the effects of flutter in the standard case are underestimated.
There are several differences of the proposed model with the classical van der Pol. An important one is that there are no stable limit cycles. This is closer to reality, for wings never reach limit cycles, they either tend to an equilibrium or have a breaking point. A second one is that, as opposed to the classical case in which orbits do not intersect, in the fractional setting, this can happen. In the fractional model, solutions can temporarily escape the basin of attraction as a direct consequence of memory effects, in contrast with the standard one. The instance where α = 0, which corresponds to a solid material with no dissipation becomes a limiting case. This is more realistic from the standpoint that all materials have some degree of dissipation.