# A Highly Accurate Approach for Aeroelastic System with Hysteresis Nonlinearity

Authors: C. C. Cui; S. X. Xie; X. C. Huang; J. K. Liu; Y. M. Chen
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

We propose an accurate approach, based on the precise integration method, to solve the aeroelastic system of an airfoil with a pitch hysteresis. A major procedure for achieving high precision is to design a predictor-corrector algorithm. This algorithm enables accurate determination of switching points resulting from the hysteresis. Numerical examples show that the results obtained by the presented method are in excellent agreement with exact solutions. In addition, the high accuracy can be maintained as the time step increases in a reasonable range. It is also found that the Runge-Kutta method may sometimes provide quite different and even fallacious results, though the step length is much less than that adopted in the presented method. With such high computational accuracy, the presented method could be applicable in dynamical systems with hysteresis nonlinearities.


---

## Body

## 1. 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 long-term 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–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–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–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–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.


## 2. Areoelastic Model with Hysteresis

Figure1 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 ahb 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.Figure 1
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].The coupled motions of the airfoil in incompressible unsteady flow can be modeled as follows [17]:
(1)c0ξ¨+c1α¨+c2+2ζξω¯U∗ξ̇+c3α̇+c4ξ+c5α+c6w1+c7w2+c8w3+c9w4+ω¯U∗2Gξ=ft,d0ξ¨+d1α¨+d2ξ̇+d3+2ζα1U∗α̇+d4ξ+d5α+d6w1+d7w2+d8w3+d9w4+1U∗2Mα=gt,where ξ=h/b is the nondimensional plunge displacement and the dot denotes the differentiation with respect to nondimensional time τ=Vt1/b with t1 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, wi’s, are given as
(2)w1=∫0te−ε1t−σασdσ,w2=∫0te−ε2t−σασdσ,w3=∫0te−ε1t−σξσdσ,w4=∫0te−ε2t−σξσdσ,with the constants ε1=0.0455andε2=0.3. The coefficients c0,…,c9andd0,…,d9 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=x1,x2,…,x8T with x1=α,x2=α̇,x3=ξ,x4=ξ̇,x5=w1,x6=w2,x7=w3,andx8=w4, the coupled state space system can be described by a set of eight first-order ordinary differential equations as follows:
(3)Ẋ=AX+BX.The expressions forA and BX 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 Figure2 is a simplified one, which was 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 Mx1 with x1=α as
(4)Mx1=x1−αf+M0,x1<αf,ẋ1>0,x1+αf−M0,x1>−αf,ẋ1<0,M0+Mfx1−αf,αf≤x1≤αf+δ,ẋ1>0,Mfx1+αf−M0,−αf−δ≤x1≤−αf,ẋ1<0,x1+M0−αf−δ1−Mf,x1>αf+δ,ẋ1>0,x1−M0+αf+δ1−Mf,x1<−αf−δ,ẋ1<0,where the preload M0, the freeplay δ, the beginning of the freeplay αf, and Mf 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.Figure 2
A hysteresis nonlinearity and illustration of the predictor-corrector algorithm for finding exact switching points.


## 3. 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.
(5)Ẋ=AiX+Bi,where AiandBii=1,…,6 are the respective coefficient matrices for subsystems. Note that the inhomogeneous terms Bi in (5) is different from BX in (3). In fact, Bi is a constant vector for each linear subsystem. The function BX in (3) is a matrix associated with the nonlinear force in pitch direction denoted as Mx1. Mx1 is divided into six portions, and in each portion, the nonlinear function is linear and only related to x1. The linear parts related to x1 in BX 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=Bi and Z=X,YT, (5) is further transformed into
(6)Ż=CiZin which
(7)Ci=AiI00,where I is an identity matrix.For each subsystem, the analytical solution of (6) can be expressed as
(8)Z=expCiτZ0,with Z0 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=expCiγ according to the procedures introduced by Zhong [25], we can obtain the state vectors by the following recurrence scheme
(9)Zn+1=TZn,where Zn represents the solution at the nth time point, that is, τ=nγ.If two consecutive vibration states,Zn and Zn+1, are located in the same subsystem, Zn+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 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 Figure2, provided that the vibration velocity is unidirectional. That means that Zn+1 passes one of the switching points such as αf,αf+δ,−αf,−αf−δ. Introduce μ to represent one of switching values αf,αf+δ,−αf,−αf−δ, which is located between the two consecutive states (Znand Zn+1). We define a coefficient as
(10)λ=μ−αnαn+1−αnto reduce the time as λγ for Zn reaching the solution at the switching point. Here, αn corresponds to the counterpart in Zn.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
(11)λ=−αnαn+1−αn,where αn denotes the angular velocity at Zn. 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(12)Zn+1=TλZn.In Figure3, 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 Zn+1) will be updated because the solution Zn+1 is corrected as Zn+1=TλZn. Generally, Zn+1=TλZn will approach the real vibration state at the switching point more closely than Zn+1=TZn. Therefore, the updated solution Zn+1=TλZn 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 Zn+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.Figure 3
Sketch of the predictor-corrector algorithm.


## 4. Numerical Examples

The system parameters are chosen as μ = 100,ah = −0.5, xα = 0.25, rα = 0.5, ζξ = ζα = 0, ω¯=0.2, M0 = Mf = 0.5°, δ = 0.1°, and αf=M0−δ/21−Mf.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.The critical flutter speed is obtained asU0 = 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 Z0=0.1,1,1,0,0,0,0,0T as an illustrative example in the following discussions.In Figure4, 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.Figure 4
Logarithmic errors of modified Euler, Newmark, or PIM results with respect to the exact solutions of the aeroelastic system atU/U0 = 0.9 in pitch direction.The bifurcation charts shown in Figure5 are obtained by the PIM and RK results, respectively, with the same time step Δt=0.1 and 0.8 < U/U0 < 0.95. As a whole, when U/U0 < 0.84 and 0.85 < U/U0 < 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/U0 = 0.84 or 0.862. The amplitudes increase as the flow speed increases beyond U/U0 > 0.86. Interestingly, there are two remarkable discrepancies between the two results as 0.82 < U/U0 < 0.84 and 0.85 < U/U0 < 0.86 or so, respectively, boxed by red dashed lines in Figure 5. When U/U0 > 0.86, both methods provide 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 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.Figure6 shows the time histories in pitch and plunge with U/U0 = 0.9, respectively. Figure 7 shows phase planes for LCOs with U/U0 = 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.Figure 6
Time histories of the aeroelastic system atU/U0 = 0.9 for (a) pitch and (b) plunge directions, respectively.
(a)
(b)Figure 7
Phase planes of LCOs for (a) pitch and (b) plunge directions, respectively, atU/U0 = 0.87.
(a)
(b)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. Figure8 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.Figure 8
Logarithmic errors of RK or PIM result with respect to the exact solutions of the aeroelastic system atU/U0 = 0.87 in pitch direction.More importantly, it should be pointed out the RK solutions are sometimes not reliable. As marked by red dashed lines asU/U0 < 0.86, in some cases the PIM provides LCOs, while the RK solutions converges to equilibriums. When U/U0 < 0.84, both the PIM and RK results imply that the system will not give rise to LCOs. As U/U0 = 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/U0 = 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.Figure 9
Time histories of the aeroelastic system atU/U0 = 0.8 (a) for pitch and (b) for plunge direction, respectively.
(a)
(b)Figure 10
Time histories of the aeroelastic system atU/U0 = 0.85 in pitch direction provided by PIM, RK, and exact solutions, respectively.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/U0 = 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.Figure 11
Time histories of the aeroelastic system atU/U0 = 0.1 obtained by the exact solutions, PIM with the time step as 0.1, and RK with the time step as 0.1, 0.01, or 0.001, respectively.


## 5. 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 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.



---

*Source: 1459830-2017-11-15.xml*