A Highly Accurate Approach for Aeroelastic System with Hysteresis Nonlinearity

State Key Laboratory of Complex Electromagnetic Environment Effects on Electronics and Information System, P.O. Box 085 Luoyang, China Department of Mechanics, Sun Yat-sen University, No. 135 Xingang Road West, Guangzhou 510275, China Basic Education Department, Jiangxi Industry Polytechnic College, Qinshanhu Road 1988, Nanchang, China Department of Mechanical Engineering, Ganzhou Institute of Technology, 88 Loutiling Rd, Ganzhou 341000, Jiangxi, China


Introduction
As more and more approaches have been developed successfully, numerical solution techniques play a pivotal role in nonlinear dynamic research.It has still been a tough job, however, to obtain a reliable numerical solution in a longterm duration for nonlinear dynamical systems [1], especially when nonsmooth nonlinearities are included [2,3].Hysteresis is one typical kind of complex nonsmooth nonlinearities, which can result from smart materials and composite structures [4][5][6].This study aims at proposing an effective approach for dynamical systems, illustrated by the aeroelastic systems of airfoils with hysteresis.
Nonlinear aeroelastic problems have received a large amount of research interests.Well known, Dowell et al. have done many significant contributions to nonlinear aeroelastics [7][8][9].Xu et al. have proposed a control strategy appropriate for the suppression of a nonlinear aeroelastic system [10].The nonlinearities of interests to aeroelasticians consist of two categories, that is, the structural and aerodynamics nonlinearities [11].Specifically, the concentrated structural nonlinearities usually refer to cubic stiffnesses, freeplay, and hysteresis.And the freeplay and hysteresis are the typical nonsmooth nonlinearities arising from aeroelastic systems.They can result in complex nonlinear dynamic responses, such as limit cycle oscillations (LCOs), bifurcations, and chaos [12][13][14][15][16].
It is tougher and even impossible to find an analytic solution of a nonlinear aeroelastic system.Alternatively, researchers pay their attention to numerical and/or semianalytical solution approaches.The nonlinear aeroelastic responses are usually generated by numerical solution techniques such as the Runge-Kutta (RK) method and the finite-difference method [9].Semi-analytical solvers are usually based on harmonic balancing [17][18][19], such as the perturbation-based method [12] and the homotopy analysis method [20] to mention a few.
To the best of our knowledge, additional theoretical and/ or computational obstacles will usually be confronted when using semianalytical approaches to solving aeroelastic systems with nonsmooth nonlinearities.For example, Liu and Dowell [21] introduced FFT technique to handle the freeplay nonlinearities when applying the harmonic balance method to solve an airfoil aeroelastic system with a control surface.Liu et al. [22] proposed a new technique to transform time domain solutions into frequency ones so that the incremental harmonic balance method is capable of analyzing the LCOs of an airfoil aeroelastic system with a hysteresis.
Though in principle numerical solution approaches can be directly implemented to nonsmooth aeroelastic systems, troublesome problems regarding reliability and accuracy of the results still have been observed.For instance, Conner et al. [23] found that switching points in piecewise linear systems can affect the computational accuracy significantly and applied Hénon's method to locate the switching points.More recently, Cui et al. [24] employed the precise integration method (PIM) initiated by Zhong [25] to propose a promising solution technique for aeroelastic systems with freeplay.In this method, the switching points are accurately and effectively determined by a predictor-corrector algorithm.Furthermore, Chen and Liu [26] implemented this method to study LCOs and bifurcation of an airfoil with both an external store and a pitch freeplay.
Hysteresis phenomena have been observed in aeroelastic systems of airfoil for decades [11,27].To some extent, hysteresis can be considered as a composite of two directional freeplays, that is, the value of the nonlinearity depends not only on the vibration displacement but also on the velocity.Therefore, it is much tougher to solve hysteresis systems.Maybe for this reason, it is rare to find effective solution approaches to do this job [12,22,28].It is necessary and worthwhile to extending our newly proposed method [24,26] to study the airfoil aeroelastic system with hysteresis.Due to the directional features of the hysteresis, we should split the nonlinear system into subsystems of much higher dimensions.Different from the routine procedures, moreover, a predictor-corrector algorithm has to be constructed to treat with the directional features when determining the switching points.

Areoelastic Model with Hysteresis
Figure 1 shows the physical model of a two-dimensional airfoil.It oscillates in pitch and plunge directions, denoted by α, positive with the nose up, and by h, positive in the downward direction, respectively.The symbol b denotes the semichord length.The elastic axis is located at a distance a h b from the midchord, and the mass center is located at a distance x α b from the elastic axis.Both distances are positive when measured towards the trailing edge of the airfoil.
The coupled motions of the airfoil in incompressible unsteady flow can be modeled as follows [17]: where ξ = h/b is the nondimensional plunge displacement and the dot denotes the differentiation with respect to nondimensional time τ = Vt 1 /b with t 1 as the real time and V as the flow speed.And ω is given by ω = ω ξ /ω α , with the symbols ω ξ and ω α as the natural frequencies of the uncoupled plunging and pitching modes, respectively.The symbol U is a nondimensional flow velocity defined as U = V/ bω α .The damping coefficients are denoted as ζ α and ζ ξ corresponding with the pitch and plunge degrees of freedom, respectively; G ξ and M α represent the nonlinear plunge and pitch stiffnesses, respectively.The other variables, w i 's, are given as with the constants ε 1 = 0 0455 and ε 2 = 0 3.The coefficients c 0 , …, c 9 and d 0 , …, d 9 are given in Liu and Dowell [17].In this paper, we consider the hysteresis nonlinearity in pitch direction and a linear spring in plunge direction.Introducing a variable vector x 7 = w 3 , and x 8 = w 4 , the coupled state space system can be described by a set of eight first-order ordinary differential equations as follows: The expressions for A and B X are given in the Appendix.Given any initial conditions, the state space vector X can be generated by time-stepping integration methods such as the RK method and Newmark method.
A freeplay nonlinearity usually results from backlash in loose and/or worn control surface hinges.When friction and backlash occur simultaneously, the nonlinearity usually exhibits a hysteresis feature.Actually, the considered model shown in Figure 2 is a simplified one, which was 2 International Journal of Aerospace Engineering widely used in aeroelastic analysis for airfoil [11,18,21,26,27].As mentioned above, mathematically, the hysteresis can be treated as a directional composite of two freeplay nonlinearities [11].Figure 2 is the general sketch of a hysteresis stiffness denoted by M x 1 with x 1 = α as where the preload M 0 , the freeplay δ, the beginning of the freeplay α f , and M f are all constants.Graphically, the hysteresis is divided into six portions shown in Figure 2. From portions (1) to (3), the direction coincides with the positive angular velocity in pitch direction, and vice versa for portions (4) to (6).The switching points, separating different portions as subsystems, depend not only on the angular displacement but also on the velocity.

PIM with Predictor-Corrector Algorithm
The PIM is chosen because it can provide highly precise numerical results for linear dynamical systems given by ordinary differential equations [25].The nonsmooth system (3) can be rewritten as a series of linear subsystems.
where A i and B i i = 1, …, 6 are the respective coefficient matrices for subsystems.Note that the inhomogeneous terms B i in ( 5) is different from B X in (3).In fact, B i is a constant vector for each linear subsystem.The function B X in ( 3) is a matrix associated with the nonlinear force in pitch direction denoted as M x 1 .M x 1 is divided into six portions, and in each portion, the nonlinear function is linear and only related to x 1 .The linear parts related to x 1 in B X is put to coefficient matrix A. As a result, the nonsmooth system (3) can be rewritten as a series of linear subsystems, that is, (5).
A simple-yet-efficient dimension expansion method can be employed to avoid computing the reverse of the coefficient matrix.With Y = B i and Z = X, Y T , ( 5) is further transformed into in which where I is an identity matrix.
For each subsystem, the analytical solution of ( 6) can be expressed as with Z 0 as the initial condition.A time step length γ is chosen repeatedly to generate solutions for a given time series, denoted as 0, γ, 2γ, …, Nγ .Calculating the exponential matrix T = exp C i γ according to the procedures introduced by Zhong [25], we can obtain the state vectors by the following recurrence scheme where Z n represents the solution at the nth time point, that is, τ = nγ.
If two consecutive vibration states, Z n and Z n+1 , are located in the same subsystem, Z n+1 can be obtained straightforwardly by implementing (9).As the recurrence computation proceeds, the vibration state sometimes inevitably switches from one subsystem to another, as illustrated in Figure 2. The latter state is in essence not a true solution unless it locates exactly at a switching point.Instead, it 3 International Journal of Aerospace Engineering can only be considered as an approximation.Though the latter state can still provide better approximation as long as the time step γ is refined, the computation cost will dramatically increase.As we know, the most outstanding merit for PIM lies in its high accuracy and efficiency.In order to maintain this merit, a predictor-corrector method has to be proposed to accurately search the vibration response happening exactly at the switching point.
Intuitively, the switches between subsystems can be judged by the fact that the vibration displacement increases (decreases) beyond the corresponding parametric region to one subsystem.For instance, the state vector switches from portion (1) to (2), or (2) to (3), etc., as shown in Figure 2, provided that the vibration velocity is unidirectional.That means that Z n+1 passes one of the switching points such as which is located between the two consecutive states (Z n and Z n+1 ).We define a coefficient as to reduce the time as λγ for Z n reaching the solution at the switching point.Here, α n corresponds to the counterpart in Z n .Different from the case with a freeplay nonlinearity [24,26], the considered aeroelastic system with a hysteresis exhibits directional features.Implicitly, a switch may happen as the direction changes, even though the displacement does not pass any one of the abovementioned switching points.This kind of switch depends on the fact that the angular velocity passes 0 rather than the displacement, for example, from portion (2) to (4) and vice versa.In this case, we define where α n denotes the angular velocity at Z n .According to (11), the coefficient converges to 1 as α n+1 approaches 0, which denotes implicitly a switch point.Once an explicit or implicit switching is revealed, an approximation of the solution at the switching point is given as In Figure 3, we can see sometimes that the vibration state switches from one subsystem to another.Once the vibration state switches from one subsystem to another, we defined a ratio as λ = μ − α n / α n+1 − α n being located between 0 and 1.The vibration state α n+1 will be the real one appearing at the switching point as long as λ = 1.This ratio is in fact the proportional mean of the switching value μ between the vibration states α n , α n+1 .With this coefficient, α n+1 (as the first entry of Z n+1 ) will be updated because the solution Z n+1 is corrected as Z n+1 = T λ Z n .Generally, Z n+1 = T λ Z n will approach the real vibration state at the switching point more closely than Z n+1 = TZ n .Therefore, the updated solution Z n+1 = T λ Z n presents a better approximation for the real vibration state at the switching point.The approximation precision can be improved by further updating the ratio λ according to the updated solution Z n+1 , and vice visa.The predictor-corrector algorithm is performed by repeating (10) (or (11)) and ( 12) until λ approaches 1 close enough.In the following calculation, the tolerance error between λ and 1 is set to be 10 −10 .

Numerical Examples
The system parameters are chosen as μ = 100, a h = −0.5, The results provided by the presented method will be compared with those by the modified Euler method, Newmark method, widely-used fourth-order RK method, and the exact results, respectively.Note that the exact solutions are obtained by solving the subsystems one by one using the theories for ordinary differential equations.The procedures are tedious due to the switches of the vibration state, because the solution at the switching points has to be determined by trial-and-error calculations.Figure 4: Logarithmic errors of modified Euler, Newmark, or PIM results with respect to the exact solutions of the aeroelastic system at U/U 0 = 0.9 in pitch direction.

International Journal of Aerospace Engineering
The critical flutter speed is obtained as U 0 = 6.2851 by linear flutter analysis method.When the speed is beyond the critical flutter speed, the system will lose its stability and all dynamic responses will increase infinitely.The initial condition is chosen as Z 0 = 0 1, 1, 1, 0, 0, 0, 0, 0 T as an illustrative example in the following discussions.
In Figure 4, we compare the numerical results provided by the PIM, modified Euler method, and Newmark method with the exact ones.We can see that when the time step chosen is much smaller, the modified Euler method improves the accuracy about one order of magnitude, while the Newmark method improves about three orders of magnitude at the beginning and produces a big error in the back.We find that the algorithm PIM with a predictor-corrector has high computational accuracy.The RK method that is widely used has better accuracy and stability than the modified Euler method and Newmark method.So in the following examples, we use the RK method compared with the proposed approach to illustrate the effectiveness.
The bifurcation charts shown in Figure 5 are obtained by the PIM and RK results, respectively, with the same time step Δt = 0 1 and 0.8 < U/U 0 < 0.95.As a whole, when U/U 0 < 0.84 and 0.85 < U/U 0 < 0.86, the pitch motions tends to the zero equilibrium even as the flow speed increases.Sudden jumping for the pitch motion is observed at U/U 0 = 0.84 or 0.862.The increase as the flow speed increases beyond U/U 0 > 0.86.Interestingly, there are two remarkable discrepancies between the two results as 0.82 < U/U 0 < 0.84 and 0.85 < U/U 0 < 0.86 or so, respectively, boxed by red dashed lines in Figure 5.When U/U 0 > 0.86, both methods provide  5 International Journal of Aerospace Engineering nearly the same solutions.To further investigate the differences, more refined simulations will be presented in the following discussion as the flow speed near the jump positions.
Figure 6 shows the time histories in pitch and plunge with U/U 0 = 0.9, respectively.Figure 7 shows phase planes for LCOs with U/U 0 = 0.87.Both the time histories and phase planes indicate that the PIM results are in nice agreement with the RK solutions.Note that the time steps are both chosen as 0.1.
In order to check the accuracy of PIM and RK method, we compare the numerical results provided by the PIM and RK with the exact ones.Figure 8 shows the logarithmic error of PIM as well as RK results compared with the exact ones.The PIM solution is much more accurate than the RK result by about eight orders of magnitude with the same time step as 0.01 (or 0.1).Importantly, the high precision can be maintained as the step length is chosen much longer for PIM.When RK is used, the accuracy improves at the earlier stage of simulation, but it remains at the same order in long-term duration.
More importantly, it should be pointed out the RK solutions are sometimes not reliable.As marked by red dashed lines as U/U 0 < 0.86, in some cases the PIM provides LCOs, while the RK solutions converges to equilibriums.When U/U 0 < 0.84, both the PIM and RK results imply that the system will not give rise to LCOs.As U/U 0 = 0.8, Figure 9 shows that the solutions in pitch direction obtained by PIM are well consistent with those by the RK method.However, observable discrepancies appear in the plunge direction.Figure 10 shows the time histories in pitch at U/U 0 = 0.85 obtained by the PIM, RK, and exact, respectively.We can see that the PIM solutions are in excellent agreement with the exact ones, whereas the discrepancy for the RK result grows to an observable magnitude during a relative long-term duration.Compared with PIM, possibly, the error for RK method is caused by its relatively lower accuracy and the approximate determining of the switching points.
Last but not least, the PIM can maintain a high precision as long as the time step is chosen at a reasonable range [25], whereas the RK method could possibly provide fallacious results even though the time step is refined significantly.Figure 11 shows the time histories presented by PIM, RK, and the exact solutions at U/U 0 = 0.1, respectively.The solutions obtained by PIM with the time step 0.1 agrees well with the exact ones.The RK solutions are very different from the PIM and exact solutions, even as the step is refined from 0.1 to 0.01 and 0.001.

Conclusion and Remarks
We have employed the PIM to propose an accurate algorithm for solving dynamical systems with hysteresis nonlinearities.The procedures are illustrated by the aeroelastic system of a two-dimensional airfoil.The solutions obtained by PIM are in excellent agreement with those by exact methods.Importantly, the accuracy can remain to be a high level as long as   International Journal of Aerospace Engineering the step length locates at a reasonable range.In principle, the accuracy can reach any desired magnitude.
Special attention has been paid on the RK method, as it has been implemented widely even considered as a benchmark for comparison.It is found that, however, the RK may sometimes provide erroneous results even if the step is very refined.With such high accuracy, the presented method could be applicable to more aeroelastic systems with nonsmooth nonlinearities, especially those with hysteresis.Moreover, it has the potential to become a benchmark for comparison in verifying some other numerical solution techniques for such systems.
It has been reported a lot on the complexity in numerical simulation of the control problems raising in nonlinear dynamical systems, especially when nonsmooth nonlinearities are included [29,30].For instance, Li et al. [31] employed Hénon's method to locate the switching points in the numerical integration procedures, for closed-loop controlled systems of airfoil with a freeplay.In the near future, we will try to investigate the feasibility and effectiveness of using the presented method in such problems.

Figure 1 :
Figure1: Sketch of an airfoil oscillating in pitch (α) direction with respect to the elastic axis and in plunge (h) direction measured from the mean position[14].

Figure 2 :
Figure 2: A hysteresis nonlinearity and illustration of the predictor-corrector algorithm for finding exact switching points.

Figure 5 :Figure 6 :
Figure 5: Bifurcation charts of the aeroelastic system obtained by RK (a) and PIM (b) with Δt = 0 1, respectively.Major differences are boxed by red dashed lines.

Figure 8 :Figure 7 :Figure 9 :Figure 11 :
Figure8: Logarithmic errors of RK or PIM result with respect to the exact solutions of the aeroelastic system at U/U 0 = 0.87 in pitch direction.

Figure 10 :
Figure10: Time histories of the aeroelastic system at U/U 0 = 0.85 in pitch direction provided by PIM, RK, and exact solutions, respectively.