- Email: [email protected]
Numerical analysis on flutter of Busemann-type supersonic biplane airfoil
Journal of Fluids and Structures 92 (2020) 102788
Contents lists available at ScienceDirect
Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs
Numerical analysis on flutter of Busemann-type supersonic biplane airfoil ∗
Hao Zhou, Gang Wang , Zhikan Liu School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, PR China School of Aeronautics, National Key Laboratory of Aerodynamic Design and Research, China
article
info
Article history: Received 31 July 2018 Received in revised form 9 June 2019 Accepted 29 October 2019 Available online xxxx Keywords: Fluid–structure interaction Busemann supersonic biplane Aeroelastic stabilities Navier–Stokes equations Flutter
a b s t r a c t The Busemann-type supersonic biplane can effectively reduce the wave drag through shock interference effect between airfoils. However, considering the elastic property of the wing structure, the vibration of the wings can cause the shock oscillation between the biplane, which may result in relative aeroelastic problems of the wing. In this research, fluid–structure interaction characteristics of the Busemann-type supersonic biplane at its design condition have been studied. A theoretical two-dimensional structure model has been established to consider the main elastic characteristics of the wing structure. Coupled with unsteady Navier–Stokes equations, the fluid–structure dynamic system of the supersonic biplane is studied through the two-way computational fluid dynamics/computational structural dynamics (CFD/CSD) coupling method. The biplane system has been simulated at its design Mach number with different nondimensional velocities. Different initial disturbance has been applied to excite the system and the effects of the position of the mass center on the system’s aeroelastic stability is also discussed. The results reveal that the stability of the airfoil in supersonic biplane system is decreased compared with that of the airfoil isolated in supersonic flow and such stability reduction effect should be given due attention in practical design. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction In 1935, Busemann (1935) proposed a supersonic biplane configuration to eliminates the wave drag by making use of the shock waves interference effect between airfoils. Later, Licher (1955) made an asymmetric design of the supersonic biplane to further increase its lift-drag ratio. However, due to many reasons, this concept was not actually applied to aircraft design at that time. In 21st century, the pursuit of supersonic commercial aircraft has led designers to refocus on this concept. Matsushima et al. (2006) and Kusunose et al. (2006) studied the aerodynamic performance of supersonic biplane through CFD method and discussed the aerodynamics efficiency decrease in off-design region. Choking and hysteresis phenomena are captured in their simulations and these phenomena were confirmed in experiments (Kuratani et al., 2007; Kuratan et al., 2009). Many efforts were made to counter against the choking problem (Yamashita et al., 2013, 2010) and to improve the off-design aerodynamic efficiency of the supersonic biplane (Hu et al., 2012). Kusunose et al. (2011) reviewed the development of the supersonic biplane, and further discussed the prospects as well as obstacles in the designing of the supersonic biplane. In recent years, Yamazaki and Kusunose (2014) and Yamazaki and Kusunose (2016) ∗ Corresponding author at: School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, PR China. E-mail addresses: [email protected] (H. Zhou), [email protected] (G. Wang), [email protected] (Z. Liu). https://doi.org/10.1016/j.jfluidstructs.2019.102788 0889-9746/© 2019 Elsevier Ltd. All rights reserved.
2
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
studied a biplane-wing/twin-body-fuselage configuration for advanced supersonic transport and evaluated its sonic boom performance. Ban et al. (2017) optimized this configuration through Kringing model and genetic algorithm. In the previous work, most researchers concentrate on the steady aerodynamic performance analysis of the rigid biplane configurations without considering the elastic property of wing structures, and very few of them studied the unsteady aerodynamic characteristics of the supersonic biplane. One relative discussion is given by Zhao et al. (2015) where the unsteady aerodynamic characteristic of the pitched supersonic biplane has been investigated and complex unsteady flow behaviors between the biplane have been revealed. When considering the elastic characteristics of the wing structure in unsteady flow, the vibration of wing structure causes the change of biplane configuration, which leads to the deviation of the designed shock interference and shock wave oscillation between the biplane. The shock wave oscillation can generate severe unsteady aerodynamic loads on the biplane and lead to complex aeroelastic problem. When the dynamic pressure is large enough, the wings in the supersonic biplane system may flutter, which even causes the failure of the wing structure. Such dangerous aeroelastic responses must be avoided in real supersonic aircraft design process. However, so far, the studies of the fluid–structure interaction properties of the supersonic biplane configuration are very limited according to published literatures and this will be the main topic of this paper. The most primary feature of the supersonic biplane system is the complex shock-boundary layer interactions and shock reflections between airfoils of the biplane. For the general supersonic aeroelastic analysis many simplified unsteady aerodynamic models can be used for solving unsteady aerodynamic loads, such as the piston theories (Zhang et al., 2009; Dowell and Bliss, 2013; Meijer and Dala, 2015) and the potential flow theories (Tseng and Merino, 1976; Yang and Sung, 1977). However, these efficient methods are not suited for analyzing the unsteady flow which involves complex shock wave motions and reflections. On the other hand, the unsteady shock motion and shock-boundary layer interactions are quite common in transonic flow (Iovnovich and Raveh, 2012; Fukushima and Kawai, 2018) and some supersonic internal flow (Nishizawa et al., 2012), and mostly accompanied with complex fluid–structure interaction phenomena (Raveh and Dowell, 2014; Poirel et al., 2018; Shinde et al., 2018). These studies involving unsteady shock boundary-layer interactions usually require solving the unsteady Navier–Stokes equations. While most researches about the supersonic biplane are based on inviscid flow, the unsteady shock interaction between the supersonic biplane should be simulated in viscous flow to involve the mechanism of more complex phenomena. The purpose of this paper is to provide a preliminary analysis of the fluid–structure interaction properties of baseline Busemann supersonic biplane and to identify the mechanisms of this aeroelastic system. In this paper, a theoretical 2D supersonic biplane aeroelastic model is established to consider the main elastic characteristics of the wings in biplane system. A loosely two-way CFD/CSD coupling procedure is used for the fluid–structure interaction simulation. The fully turbulence unsteady Reynolds average Navier–Stokes (RANS) equations are solved to include the shock boundary-layer interactions between the airfoils of the biplane. Different types of initial disturbance are applied to excite the system and stability analysis is applied to the airfoil in supersonic biplane system to reveal the effects of the shock interactions between the biplane on the wing’s aeroelastic stability. The effects of the position of the mass center on the system’s aeroelastic stability are also discussed. 2. Flow equations In this paper, we use an in-house CFD solver HUNS3D (Liu et al., 2018) to solve the fully turbulent flow. Spalart– Allmaras (SA) turbulence model is used in this study. The fluid governing equations in the Arbitrary Lagrangian Eulerian (ALE) form for a bounded control volume Ω with boundary ∂ Ω , are expressed in integral form as below:
∂ ∂t
∫∫∫ Ω
U dΩ +
∫∫ ∂Ω
[
F (U , V grid ) − G(U ) · ndS =
]
∫∫∫ Ω
H dΩ
(1)
where U = {ρ, ρ u, ρv, ρw, ρ E , ρ˜ ν}T and F , G are respectively represent convective and viscous fluxes terms. H is the source term, and V grid denotes the grid velocity. The equations are discretized on unstructured hybrid meshes with the cell-centered finite volume method, hence, for each mesh cells Ωi , Eq. (1) can be rewritten as: d(Ωi U i ) dt
+
∑ [
F (U im , V grid ) − G(U im ) S im = Ωi H i
]
(2)
m∈N(i)
The convective fluxes are discretized using Roe’s (1981) method, and the viscous fluxes term is obtained by reconstructed central scheme. For upwind schemes the second order accuracy is achieved via reconstructing the solution following Frink’s interpolation method (Frink, 2008). When the discrete fluxes of each integration surface are obtained, the semi-discretized form of the governing equations are then integrated in time, Eq. (2) is rewritten as: d(Ωi U i )
+ R i (U i ) = 0 (3) dt Dual time stepping method has been applied for solving Eq. (3). The time derivative is discretized by a second order scheme given by: 3(Ωi U in+1 ) − 4(Ωi U ni ) + (Ωi U ni −1 ) 2∆t
+ R i (U ni +1 ) = 0
(4)
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
3
By introducing a pseudo-time veritable t ∗ into Eq. (4), it can be expressed as Eq. (5), which could be solved through many steady flow accelerating methods implemented in HUNS3D (Wang et al., 2012; Liu et al., 2018). dU ni +1
+
dt ∗
1
Ωin+1
R ∗i (U ni +1 ) = 0
(5)
Since the RBF mesh deformation technique (Wang et al., 2018) has been applied to update the grids every step of the simulation, V grid must be calculated in line with the time integration scheme by Eq. (6) to satisfy the Geometric Conservation Law (GCL) where the Vi is control volume, first introduced by Thomas and Lombard (2008). The detailed technique for solving Eq. (6) can be found in Blazek (2005).
∑
3Ωin+1 − 4Ωin + Ωin−1
m
2∆t
(V grid · nS)nm+1 =
(6)
3. Structural dynamics model In the field of aeroelasticity, the classical two-degree-of-freedom (2DOF) model (Isogai, 1979) is widely used in many aeroelastic studies. In this model, the primary elastic characteristics of wing are modeled as pitch and plunge of the airfoil. Although, as early as in the 1950s, some literatures (Biot, 1956; Glasgow and Cronin, 1960) has pointed out that the chordwise mode may be important for the low aspect-ratio thin wings in supersonic aeroelastic analysis, there are still many studies about supersonic aeroelasticity based on 2DOF model after that, such as supersonic compressor rotor flutter (Bendiksen, 1986), supersonic flutter active control and suppression (Marzocca et al., 2002; Rao et al., 2006) as well as some studies about nonlinear flutter phenomena in supersonic flow (Firouz-Abadi et al., 2013; Zhou et al., 2013; Marzocca et al., 2005). In spite of its simplicity, this simplified model should give acceptable predictions of the main dynamic behavior of some structures in supersonic flow to some extent. Therefore, concerning our main thrust and intentions are to explore the basic aeroelastic properties of this novel configuration, it is appropriate to start with the 2DOF model for the preliminary analysis within certain assumptions. We focus on the most primary elastic properties of wing structure. Here, a typical wing section is extracted from the three-dimensional (3D) Busemann wing model, as shown in Fig. 1(a). Since the purpose in this paper is to investigate the theoretical aeroelastic property and mechanism of this configuration at certain Mach region, the structure supporting connection of the two wings at wing root is neglected for simplification and no linkages exist between the upper and lower wings. The 3D effect is not considered at this stage. It should be noticed that when the aspect ratio is low enough, the 3D effects should be considered, and the participation of chordwise mode needs to be studied in future research for practical design. The two-dimensional (2D) model has more theoretical meaning for exploring the basic aeroelastic properties and mechanism. With all these assumptions and intentions, each airfoil of the wing section can be considered as a typical 2DOF airfoil structure model in theoretical aeroelasticity studies and the upper airfoil and lower airfoil interacts only through aerodynamic forces. The structural parameters were originally defined by Isogai (1979). A sketch of this model is shown in Fig. 1(b), where h denotes the displacement away from balance position in z direction and α denotes the angular displacement. E and G represent the elastic axis and mass center respectively. b denotes the half chord length of airfoil and ab denotes the distance between midpoint of the chord and elastic axis. When a is negative, the elastic axis is at the front of the midpoint of the chord to simulate the sweep effect. When a is 0, there is no sweepback for the wing. The dimensional dynamic equations of the system can be easily obtained through the balance of force in z direction and the moment balance around elastic axis, as shown in Eq. (7). The Sα and Iα are defined as mxα b and mrα2 b2 respectively in unit span of the airfoil and m is the mass of airfoil. The Kh and Kα are the elastic coefficient of the springs. The L and Mea are the aerodynamic lift and moment around the elastic axis.
⎧ ∂ 2h ∂ 2α ⎪ ⎪ ⎨ m 2 + Sα 2 + Kh h = −L ∂t ∂t 2 2 ⎪ ∂ h ∂ α ⎪ ⎩Sα + Iα 2 + Kα α = Mea 2 ∂t ∂t
(7)
The equations can be nondimensionalized by the method mentioned in Alonso and Jameson (1994) and Liu et al. (2001), and the time in structure system is also nondimensionalized by multiplying a reference frequency ωα which is also defined as the natural frequency of rotational degree of freedom. The final nondimensional form of Eq. (7) can be written in matrix form as: M
∂ 2q + Kq = F ∂τs2
(8)
where:
[ M =
1
xα
xα
rα2
]
,K =
[
(ωh /ωα ) 0
1 0 −Cl h/b , F = V ∗2 ,q = 2Cm α rα2 π
]
[
]
[
] (9)
4
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 1. 2D airfoil structural model.
V∗ =
V
m
√ ,µ = ωα b µ πρ∞ b2
(10)
The ωh and ωα are the natural frequency of pitch and plunge defined by ωα2 = Kα /Iα , ωh2 = Kh /m. µ is defined as mass ratio. V ∗ is the nondimensional velocity also known as flutter index. Since the time in structure dynamic equation has been nondimensionalized by ωα , the structure dynamic time step ∆τs must match with the flow nondimensional time step ∆τ through Eq. (11) when solving these equations in time domain where Ma is the Mach number.
∆τs =
2Ma · ∆τ V∗
√
(11)
µ
For the biplane system, shown in Fig. 2, two structurally unconnected 2D airfoil model can be simply posed together as follow:
⎤ ∂ 2 qup [ ⎢ ∂τs2 ⎥ ⎥ + K up ·⎢ ⎣ ∂ 2q ⎦ 0 low ∂τs2 ⎡
[
M up 0
0 M low
]
0 K low
] [ ] [ ] q F up · up = qlow
F low
(12)
Noted as: Ms
∂ 2 qs + K s qs = F s ∂τs2
Which is the governing equation of the supersonic biplane FSI system.
(13)
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
5
Fig. 2. Busemann supersonic biplane. Table 1 Structure parameters of NACA64A010 at transonic region.
µ
ωh /ωα
rα2
xα
a
α0 (rad)
60
1
3.48
1.8
−2.0
0
Table 2 Structure parameters of NACA64A010 at supersonic region.
µ
ωh /ωα
rα2
xα
a
α0 (rad)
75
0.5
0.75
0.25
0
0.0873
4. FSI Coupling method By introducing state variables Q = [qTup ,qTlow ,∂ qTup /∂τs ,∂ qTlow /∂τs ] = [qTs ,∂ qTs /∂τs ], Eq. (13) can be written in state space form: ∂Q = A · Q + B · F s (Q , τs ) (14)
∂τs
Where
[ A=
O 1 −M − s Ks
I 1 −M − s Gs
]
[ B=
O 1 M− s
] (15)
A 4th order semi-implicit linear multi-step scheme (Zhang et al., 2007) has been used for solving Eq. (14). The structural part is solved implicitly, the general aerodynamic part is solved explicitly. Its accuracy and stability are proved better than classic Adams predictor–corrector procedure in Zhang et al. (2007). The specific marching steps are as follow:
⎧ F sn = F s (Q n , τsn ) ⎪ ⎪ ⎪ ⎪ f n = A · Q n + B · F sn ⎪ ⎪ ⎪ ⎪ ) ∆τs ( ⎨ Q n+1 = Q n + 55f n − 59f n−1 + 37f n−2 − 9f n−3 24 ⎪ ⎪ ⎪ f n+1 = A · Q n+1 + B · F sn+1 ⎪ ⎪ ⎪ ⎪ ) ∆τs ( ⎪ ⎩Q 9f n+1 + 19f n − 5f n−1 + f n−2 n+1 = Q n +
(16)
24 When Q are solved at each time step, the new wall boundary condition can be calculated according to. Then the mesh is updated by the mesh deformation methods (Wang et al., 2018) for the flow field computation at the next time step. 5. FSI system verification To validate the FSI procedure in this research, a standard transonic 2D flutter case of NACA64A010 (Isogai, 1979) has been studied. The structure parameters are listed in Table 1. This case is studied extensively for Euler equation while the results for Reynolds-Averaged N–S solution are much fewer. Some of the results can be found in Bohbot et al. (2001), Prananta et al. (1998) and Marti and Liu (2017). The N–S flutter boundary for Mach number 0.75–0.95 calculated in this paper is shown Fig. 3, which agrees well with others’ calculations. The multiple equilibrium points discussed in Marti and Liu (2017) are also found at Mach 0.875, as shown in Fig. 4. All these results prove that the FSI procedure implemented in this paper is correct. Since the study in present paper is in supersonic regime, the supersonic unsteady aerodynamic force also need to be validated. However, very few supersonic calculations of airfoil flutter are carried out in viscous flow directly. Here, an
6
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 3. Flutter boundary of NACA64A010 airfoil model.
Fig. 4. Heaving and pitching responses for RANS results of NACA64A010 at Ma = 0.875, Re = 1.0e7.
inviscid 2-D case of NACA64010 airfoil flutter in supersonic Mach region (Zhang et al., 2009) and a standard 3-D flutter case about AGARD445.6 wing (Yates, 1988) in viscous flow, which includes a supersonic Mach condition, is used to validate the calculation in supersonic Mach region. In the 2-D case we compare our N–S results with the inviscid results of Zhang et al. (2009) in Fig. 5(a). The structure parameters are listed in Table 2. The results agree well with each other since the viscous effect in this Mach region is not very obvious, this verifies the aeroelastic model in supersonic Mach region. In the 3-D case, the same unsteady CFD solver is employed and the CFD codes are loosely coupled with an open source CSD solver CALCULIX (Dhondt, 2018). The detailed coupling procedure is listed in Appendix and the date transformation method can refer to Mian et al. (2014). The flutter boundaries obtained by presented codes as well as experiment results and others’ calculations (Chen and Liu, 1985; Lee-Rausch and Batina, 1993; Cavagna et al., 2007; Yang and Mavriplis, 2007; Liu et al., 2001; Im et al., 2012) are shown in Fig. 5(b). The flutter boundary at supersonic Mach condition presented in this paper shows better agreement with the experiment data compared with most of others’ simulation results. These cases should validate the unsteady supersonic aerodynamic force in this study. More supersonic unsteady aerodynamic validation cases can be found in our previous work (Wang et al., 2015). 6. Flutter analysis of supersonic biplane 6.1. Structure parameters and computing mesh In this section, the developed FSI analysis system is applied to simulate the fluid–structure interactions of Busemann supersonic biplane section model. The thickness-chord ratio (t/c) of the Busemann biplane is selected as 0.10, with the gap
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
7
Fig. 5. Validation cases involving supersonic Mach region. Table 3 Structure parameters of Busemann biplane (each airfoil).
µ
ωh /ωα
rα2
xα
a
α0 (rad)
75/2
0.5
0.75
0.25 and 0
0
0
Table 4 Grid convergence. Cell number
CL_low
CM_low
CL_up
CM_up
43 011 57 055 102 535
−1.7111e−1 −1.7110e−1 −1.7110e−1
−2.8229e−3 −2.8210e−3 −2.8205e−3
1.7109e−1 1.7108e−1 1.7107e−1
2.8280e−3 2.8255e−3 2.8248e−3
of half chord length (z/c = 0.5). Minimum drag flow condition was obtained at Mach 1.7 of zero angle of attack (Kusunose et al., 2011), which will be the free-steam condition of this study. The Reynolds number are set to 10 million. The structure parameters of a 2D supersonic flutter case for NACA64A010 airfoil (Zhang et al., 2009) are used as a reference for each airfoil of the biplane. The structure parameters of the supersonic biplane are listed in Table 3. With a = 0, the elastic axis is at the midpoint of the chord, hence the swept effect is not considered in this paper. To consider the effect of sweepback further research work is needed. Noting that, the mass ratio is reduced by half for each airfoil to make sure the structure mass remains the same compared with traditional supersonic single wing. In the real engineering practice, it is much easier to adjust the mass center than the elastic axis, then two different mass center position (xα = 0.25 and xα = 0) are tested to explore its effect to system stability. An O-type hybrid grid shown in Fig. 6 is used for the simulation. First layer height of the grid in the normal direction corresponds to maximal y+ <= 1.0. The convergence of the grid has been studied in Table 4. Considering the computing cost and the force coefficients convergence, the medium size grid is chosen for the simulation. 6.2. Time domain FSI analysis The typical method to determine the stability of an aeroelasticity system is applying a small disturbance to the balanced system and monitoring the response in time domain. However, the aerodynamic loads on each airfoil of the biplane are not symmetrical, making the airfoils unbalanced at zero angle of attack. The static aeroelasticity problem is not considered in this study, however it is an important research content, here the system is balanced at zero angle of attack by external forces to exclude the static aeroelastic effects. A steady simulation has been performed at first to obtain the initial axial force and momentum around the elastic axis of each airfoil of the biplane. When starting the unsteady simulation, the opposite external forces and moments are applied to the upper and the lower airfoil respectively to cancel its initial steady aerodynamic loads, thereby ensuring that the system’s static aeroelasticity balance is excluded in this research. For the unsteady simulation in this paper, about 80 to 100 time-steps per flutter period gives good convergence of the time domain response. 6.2.1. Time domain analysis for xα = 0.25 This parameter is chosen to make the flutter easier to be excited and observed according to the experience of traditional airfoil flutter. As for stability, the key problem is to determine the critical V* of the biplane system at given condition. Two
8
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 6. Typical grid used for FSI analysis.
tentative V* are tested to give the boundary of V*. The time domain responses of two different V* (0.5, 0.7) with the same initial condition E0 = [0 0 0 0 0 0.01 0 0.01] are shown in Fig. 7. For this initial condition, the angle velocity disturbance has the same magnitude and direction. It can be seen that as V* grows, the system goes from stable to divergent, which means that the critical V* is between 0.5 and 0.7. For the stable response solution at V* = 0.5 in Fig. 7(a), (b) the phase of pitch and plunge response for upper and lower airfoil remains almost the same at entire time history and the energy of structure system is dissipated into the flow every cycle of response. In the unstable response history in Fig. 7(c), (d), some kind of phase transition occurred. When the system becomes fully divergent, the upper and lower airfoil has the nearly opposite phase and the structure system absorbs energy from the flow causing the h/b increases to an unacceptable level. The cases of same velocity index are also tested at different initial condition E0 = [0 0 0 0 0 0.01 0 −0.01], shown in Fig. 8. For this initial condition, the angular velocity disturbance has the same magnitude but in opposite direction, which makes the time response has opposite phase at the beginning. In Fig. 8(a), the similar phase transition occurs in the response of α . This implies that when the system is stable the response of α may tends to be in phase, and when the system is unstable the response of seems tends to be in opposite phase. The phase transition is not found in the calculated time series of h/b partly due to h/b is not excited in the initial condition and the limited simulation time. For the results of V* = 0.7 shown in Fig. 8(c), (d), the applied initial condition accelerated the divergence of the biplane system which also exemplifies our point. To determine the critical speed more precisely, the V* section of (0.5, 0.7) needs to be divided into more couple of small sections. The final critical V* is found near 0.61, The time domain solution for V* = 0.65, 0.625 and 0.61 with different initial condition are presented in Figs. 9–14 to indicate the critical state of the system. For the V* = 0.65 shown in Fig. 9, the amplitude of decreases gradually at first and then increase to divergence accompanied by a very clear phase transition phenomenon. The amplitude of h/b increase at entire time history. When applying the antiphase disturbance, shown in Fig. 10, both α and h/b increase to divergence directly. Noticing that when the opposite phase disturbance is applied, the α and h/b have the same trend of convergence or divergence, and when the in-phase disturbance is applied, some component of the system may have a convergent trend in the first period of time but will become divergent in the end due to the phase transition process mentioned previously. These phenomena imply that the initial condition does not affect the system’ stability however antiphase disturbance is recommended for time saving when to determine the system’s stability boundaries. In Fig. 11(a), the for V* = 0.625 has obvious convergence pattern in the simulated time period, but the system is not stable since the time response of h/b shown in Fig. 11(b) is divergent. The system’s stability can be seen more clearly when the antiphase disturbance is applied, shown in Fig. 12, which is near the critical state but slightly divergent. It should be noticed that the phenomenon of beats can be clearly observed in Figs. 11 and 13. The beats phenomenon usually occurs when the flutter is about to happen where the frequency of two mode is close to each other. The beats phenomenon proved that these two conditions are close to the critical flutter point. Figs. 13 and 14 actually draws the lower boundary of the critical V*. The more precise critical V* can be obtained by the damping interpolation method mentioned in Marti and Liu (2017). Limited to our purpose here, the actual flutter boundaries are not precisely calculated in this paper and the complete flutter boundaries of the supersonic biplane require a lot of computations which will be carried out in our future work. The main difference between supersonic biplane and its isolated airfoils is the shock interactions and reflections between the biplane, as shown in Fig. 15, which may make the aeroelastic system more complicated than the isolated airfoil. Here, the same stability analysis is also applied to the isolated upper airfoil of the supersonic biplane to identify the effects of the shock interactions and reflections between the biplane to the airfoil’s aeroelastic stability. Since the structure
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
9
Fig. 7. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0.25, E0 = [000000.0100.01].
is symmetrical, the upper airfoil is selected. For the upper airfoil of the supersonic biplane isolated in the supersonic flow, the critical V* is about 0.83, shown in Fig. 16, nearly 27% higher than that of supersonic biplane at same flow condition, indicating that the aeroelasticity stabilities of the airfoil in supersonic biplane system been greatly decreased due to the shock interactions and reflections between the airfoils. 6.2.2. Time domain analysis for xα = 0 For xα = 0, the mass center is on the elastic axis, and M becomes a diagonal matrix. The pitch and plunge are no longer coupled in structure dynamics. According to the discussion in the previous section, it is easier to determine the system’s stability by antiphase disturbance, hence all the analysis about xα = 0 is carried out under antiphase disturbance (E0 = [0 0 0 0 0 0.01 0 −0.01]). Starting the simulation from V* = 0.61, the time domain results are shown in Fig. 17 which is quite different from Fig. 14. In Fig. 17, the time domain response can be divided into two sections, the first section is before τ = 200 where the α and h/b respond around zero. In the first section, the magnitude of α increases gradually and then reaches a nearly fixed value for a short period of time. After that, the α decreases a little and then increases rapidly, goes into the next section and the response of h/b is similarly. In first section, the α , h/b of each airfoil of the biplane responses in the opposite phase respectively, the system is actually divergent. In the second section, the α and h/b deviate from zero and respond around a new position in quite different form. The α in the second section is divergent while the h/b is convergent and each airfoil of the biplane responses in the same phase. It implies that the initial equilibrium point is no longer stable and the flow pattern transformation occurs. Although the response in second section may actually not happen in real structure due the large bending deformation, these high amplitude oscillations are shown to give a more complete and intuitive feeling of the whole divergent process, and the response of ideal structure model is an interesting topic which may help to understand the mechanism of such aeroelastic system when the shock interactions are involved.
10
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 8. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0.25, E0 = [0 0 0 0 0 0.01 0 −0.01].
Fig. 9. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.65.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
11
Fig. 10. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 −0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.65.
Fig. 11. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.625.
Fig. 12. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 −0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.625.
12
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 13. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.61.
Fig. 14. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 −0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.61.
Fig. 15. Flow pattern of supersonic biplane and its isolated airfoil at Ma1.7, Re10e7.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
13
Fig. 16. Heaving and pitching for RANS results of isolated upper airfoil of supersonic biplane at critical condition V* = 0.83, Ma = 1.7, Re = 1.0e7, xα = 0.25, E0 = [0 0 0 0.01].
Fig. 17. Heaving and pitching for RANs results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0.0 for different initial condition, V* = 0.61, E0 = [0 0 0 0 0 0.01 0 −0.01].
Here we extract a cycle of flow field form these two sections respectively. For the first section the snapshot is extracted from τ = 135 and for the second section, it is extracted from τ = 284. The results are shown in Figs. 18 and 19. In Fig. 18, the flow pattern is not stable, it changes as the airfoils oscillate. It can be seen clearly that a normal shock occurs near the throat of the biplane and moves forward against the flow. As the normal shock moves forward it becomes stronger. When it reaches somewhere ahead of the throat but still holds in the channel constructed by the biplane, it moves backward and becomes weaker and weaker until disappear and then another cycle starts. The severe shock oscillation causes strong pressure fluctuation which aggravated the structure instability and the unstable flow pattern finally turns into a choked flow. Fig. 19 shows this choked flow pattern for the biplane, which is much more stable than that of last section since the flow is nearly not change in its whole oscillation cycle. However, due to the strong shock ahead of the airfoils the shock interactions between the drag reduction effects of supersonic biplane is no longer exist. By decreasing V*, the final critical V* is found between 0.25 and 0.30. The time domain results of V* = 0.25 and V* = 0.3 as well as another divergent state V* = 0.4 are shown in Figs. 20–22. In these time domain results, the upper and lower airfoil are not responding symmetrically, the mechanism of such phenomena needs further research. Noticing that the critical V* for xα = 0.25 is about 0.61, the results here indicating that the stability of the system is greatly reduced as the mass center moving to the elastic axis. This is inconsistent with the general experience of traditional airfoil flutter that the critical flutter index usually increases when reducing the distance between mass center and elastic axis of the wing section. The same simulations are also performed on the isolated upper airfoil of the biplane and its critical V* is found near 1.4, shown in Fig. 23, which is higher than that of xα = 0.25 (shown in Fig. 16) and agree with our general experience
14
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 18. Flow field of Busemann biplane in an oscillation cycle at τ = 135 ∼ 142, Ma = 1.7, Re = 1.0e7, xα = 0, V* = 0.61.
on wing flutter. What is more, it shows that the shock interaction between the biplane destabilizes the system much severely when the mass center is moved to the elastic axis.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 19. Flow field of Busemann biplane in an oscillation cycle at τ = 284 ∼ 291, Ma = 1.7, Re = 1.0e7, xα = 0, V* = 0.61.
15
16
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 20. Heaving and pitching for RANs results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0 for different initial condition, V* = 0.25, E0 = [0 0 0 0 0 0.01 0 −0.01].
Fig. 21. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0 for different initial condition, V* = 0.30, E0 = [0 0 0 0 0 0.01 0 −0.01].
Fig. 22. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0 for different initial condition, V* = 0.40, E0 = [0 0 0 0 0 0.01 0 −0.01].
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
17
Fig. 23. Heaving and pitching for RANS results of isolated upper airfoil of supersonic biplane Ma = 1.7, Re = 1.0e7, xα = 0, E0 = [0 0 0 0.01], V* = 1.4.
Fig. 24. HUNS3D/Calculix coupling procedure.
7. Conclusions In this paper, the fluid–structure interaction characteristics of the Busemann supersonic biplane has been studied. A theoretical two-dimensional biplane structure model is established to consider the main elastic property of the structure. The structure dynamic equations coupled with fully turbulent URANS equations are solved through a loosely coupled CFD/CSD method. At the designed flight condition of the supersonic biplane, the time domain responses of the structure at different flutter indexes with different initial disturbance are studied. The shock interaction between the biplane is critical for drag reduction, however, the time domain simulations reveal that the aeroelastic stability of airfoil in the supersonic biplane system can be greatly decreased compared with that of airfoil isolated in the supersonic flow due to the shock interactions and reflections between the biplane. When the disturbance is applied, the shock wave deviates from its original interference point and reflected between the biplane. The reflection point moves as the airfoils vibrate, resulting in severe change in pressure distribution on the airfoils. Such unsteady load is much severe than that on the isolated supersonic airfoil since the shock position of the isolated airfoil in unsteady supersonic flow is nearly unchanged, which should be the main reason for the decrease of the system stability. Different initial disturbance is applied to excite the dynamic system and we find that the initial disturbance affects the
18
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
system’s time domain response history but does not have much effects on the final divergence pattern. For the divergent state, α and h/b of the upper and lower airfoil tends to respond in opposite phase and the anti-phase disturbance can accelerate the divergent process of the system. The position of mass center has great effects on the system’s aeroelastic stability. The critical flutter index of the airfoil in supersonic biplane system is greatly decreased when the mass center is moved forward to its elastic axis. This result contradicts the general engineering experience that moving the mass center forward can increase the critical flutter speed of the wing, and our calculations on the isolated upper airfoil agrees with this experience. The contradictory reveals the effects of the unsteady shock wave on the airfoil’s aeroelastic stability and implies that some kind of different flutter mechanism may exist in this aeroelastic system, which requires further study. The results in the presented work show the theoretical aeroelastic properties in 2D model of supersonic biplane system, some large deformation in divergent state may actually not happen in real structure due to the structural limits, these results are shown to provide a more complete and intuitive feeling of the whole divergent process. Limited by the current research intention and content, the 3D effect as well as the chordwise mode is not considered at this stage and it is important to reinvestigate and confirm these theoretical properties through 3D cases in future work once the real structure parameters are obtained from structure engineering. Acknowledgments This study was sponsored by the seed Foundation of Innovation and Creation for Graduate Students in Northwestern Polytechnical University, PR China (No. ZZ2019001). And was also co-supported by the National Natural Science Foundation of China (No. 11772265) and the ‘‘111’’ Project of China (No. B17037). Appendix. Coupling procedure of HUNS3D/Calculix codes The coupling procedure is performed in a staggered form (see Fig. 24). The structure deformation is sent from Calculix to HUNS3D and the pressure distribution along the fluid–structure interface is sent back to Calculix at each time step. Since the coupling procedure is in a loosely coupled way, the time steps should be small enough to achieve both the stability and convergence condition. References Alonso, J., Jameson, A., 1994. Fully-implicit time-marching aeroelastic solutions. In: 32nd Aerospace Sciences Meeting and Exhibit. American Institute of Aeronautics and Astronautics. Ban, N., Yamazaki, W., Kusunose, K., 2017. Low-Boom/Low-Drag Design optimization of innovative supersonic transport configuration. J. Aircr. 55 (3), 1–11. Bendiksen, O., 1986. Role of shocks in transonic/supersonic compressor rotor flutter. AIAA J. 24 (7), 1179–1186. Biot, M.A., 1956. The divergence of supersonic wings including chordwise bending. J. Aeronaut. Sci. 23 (3), 237–251. Blazek, J., 2005. Computational Fluid Dynamics: Principles and Applications, second ed. Bohbot, J., Garnier, J., Toumit, S., Darracq, D., 2001. Computation of the flutter boundary of an airfoil with a parallel Navier–Stokes solver. In: 39th Aerospace Sciences Meeting and Exhibit. American Institute of Aeronautics and Astronautics. Busemann, A., 1935. Aerodynamic Lift at Supersonic Speeds, twelveth ed. Cavagna, L., Quaranta, G., Mantegazza, P., 2007. Application of Navier–Stokes simulations for aeroelastic stability assessment in transonic regime. Comput. Struct. 85 (11–14), 818–832. Chen, P.-C., Liu, D., 1985. A harmonic gradient method for unsteady supersonic flow calculations. J. Aircr. 22 (5), 371–379. Dhondt, G., 2018. Calculix. [Online]. Available: http://www.calculix.de/. Dowell, E.H., Bliss, D.B., 2013. New look at unsteady supersonic potential flow aerodynamics and piston theory. AIAA J. 51 (9), 2278–2281. Firouz-Abadi, R.D., Alavi, S.M., Salarieh, H., 2013. Analysis of non-linear aeroelastic response of a supersonic thick fin with plunging, pinching and flapping free-plays. J. Fluids Struct. 40, 163–184. Frink, N.T., 2008. Upwind scheme for solving the Euler equations on unstructured tetrahedral meshes. AIAA J. 30 (1), 70–77. Fukushima, Y., Kawai, S., 2018. Wall-modeled large-eddy simulation of transonic airfoil buffet at high Reynolds Number. AIAA J. 56 (6), 2372–2388. Glasgow, J.C., Cronin, J.L., 1960. Chordwise-bending effects on wing-vibration modes. J. Aerosp. Sci. 27 (12), 962–963. Hu, R., Jameson, A., Wang, Q., 2012. Adjoint-based aerodynamic optimization of supersonic biplane airfoils. J. Aircr. 49 (3), 802–814. Im, H., Chen, X., Zha, G., 2012. Prediction of a supersonic wing flutter boundary using a high Fidelity detached Eddy simulation. In: 50th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition. American Institute of Aeronautics and Astronautics. Iovnovich, M., Raveh, D.E., 2012. Reynolds-averaged Navier–Stokes study of the Shock-Buffet instability mechanism. AIAA J. 50 (4), 880–890. Isogai, K., 1979. On the transonic-dip mechanism of flutter of a swept back wing. AIAA J. 17 (7), 793–795. Kuratan, N., Ozak, S., Obayash, S., Ogawa, H., Matsuno, Toshihiro, Kawazoe, Takashi, 2009. Experimental and computational studies of low-speed aerodynamic performance and flow characteristics around a supersonic biplane. Trans. Japan Soc. Aeronaut. Space Sci. 89–97. Kuratani, N., Ogawa, T., Yamashita, H., Yonezawa, M., Obayashi, S., 2007. experimental and computational fluid dynamics around supersonic biplane for sonic-boom reduction. In: 13th AIAA/CEAS Aeroacoustics Conference (28th AIAA Aeroacoustics Conference), Aeroacoustics Conferences. Kusunose, K., Matsushima, K., Goto, Y., Yamashita, H., T. B. T.-A. A. S. M., Nakano, E., 2006. A fundamental study for the development of boomless supersonic transport aircraft, In: 44th AIAA Aerospace Sciences Meeting and Exhibit. Kusunose, K., Matsushima, K., Maruyama, D., 2011. Supersonic biplane A review. Prog. Aerosp. Sci. 47 (1), 53–87, Elsevier Ltd. Lee-Rausch, E.M., Batina, J., 1993. Calculation of AGARD wing 4456 flutter using Navier–Stokes aerodynamics. In: 11th Applied Aerodynamics Conference. American Institute of Aeronautics and Astronautics. Licher, R.M., 1955. Optimum Two-Dimensional Multiplanes in Supersonic Flow. Report No.SM-18688, Douglass Aircraft Co.. Liu, F., Cai, J., Zhu, Y., Tsai, H.M., Wong, A.S.F., 2001. Calculation of wing flutter by a coupled fluid-structure method. J. Aircr. 38 (2), 334–342. Liu, Y., Wang, G., Ye, Z., 2018. Dynamic mode extrapolation to improve the efficiency of dual time stepping method. J. Comput. 352, 190–212.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
19
Marti, F., Liu, F., 2017. Multiple equilibrium points of airfoil flutter in viscous flow. In: In 55th AIAA Aerospace Sciences Meeting. American Institute of Aeronautics and Astronautics. Marzocca, P., Librescu, L., Kim, D.-H., Lee, I., 2005. Supersonic flutter and LCO of airfoils via CFD/analytical combined approach. In: 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. American Institute of Aeronautics and Astronautics. Marzocca, P., Librescu, L., Silva, W.A., 2002. Flutter, postflutter, and control of a supersonic wing section. J. Guid. Control. Dyn. 25 (5), 962–970. Matsushima, K., Kusunose, K., Maruyama, D., Matsuzawa, T., 2006. Numerical design and assessment of a biplane as future supersonic transport. In: 25th International Congress of the Aeronautical Sciences. Meijer, M.-C., Dala, L., 2015. Generalized formulation and review of piston theory for airfoils. AIAA J. 54 (1), 17–27. Mian, H.H., Wang, G., Ye, Z.Y., 2014. Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach. J. Fluids Struct. 49, 186–201. Nishizawa, U., Kameda, M., Watanabe, Y., Yamamoto, S., 2012. Computational simulation of shock oscillation around a supersonic air-intake, In: Collection of Technical Papers - 36th AIAA Fluid Dynamics Conference, Vol. 1, pp. 406–415. Poirel, D., Goyaniuk, L., Benaissa, A., 2018. Frequency lock-in in pitch–heave stall flutter. J. Fluids Struct. 79, 14–25. Prananta, B.B., Hounjet, M.H.L., Zwaan, R.J., 1998. Two-dimensional transonic aeroelastic analysis using thin-layer navier–stokes method. J. Fluids Struct. 12 (6), 655–676. Rao, V.M., Behal, A., Marzocca, P., Rubillo, C.M., 2006. Adaptive aeroelastic vibration suppression of a supersonic airfoil with flap. Aerosp. Sci. Technol. 10 (4), 309–315. Raveh, D.E., Dowell, E.H., 2014. Aeroelastic responses of elastically suspended airfoil systems in transonic buffeting flows. AIAA J. 52 (5), 926–934. Roe, P.L., 1981. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys. 43 (2), 357–372. Shinde, V., McNamara, J.J., Gaitonde, D.V., Barnes, C.J., Visbal, M.R., 2018. Panel flutter induced by transitional shock wave boundary layer interaction. In: 2018 Fluid Dynamics Conference. Thomas, P.D., Lombard, C.K., 2008. Geometric conservation law and its application to flow computations on moving grids. AIAA J. 17 (10), 1030–1037. Tseng, K., Merino, L., 1976. Fully unsteady subsonic and supersonic potential aerodynamics of complex aircraft configurations for flutter applications. In: 17th Structures, Structural Dynamics, and Materials Conference. American Institute of Aeronautics and Astronautics. Wang, G., Chen, X., Liu, Z., 2018. Mesh deformation on 3D complex configurations using multistep radial basis functions interpolation. Chin. J. Aeronaut. 31 (4), 660–671. Wang, G., Jiang, Y., Ye, Z., 2012. An improved LU-SGS implicit scheme for high Reynolds number flow computations on hybrid unstructured mesh. Chin. J. Aeronaut. 25 (1), 33–41. Wang, G., Zeng, Z., Suo, Q., 2015. Trajectory simulation of a spinning projectile based on variable step size CFD/RBD method. In: AIAA Atmospheric Flight Mechanics Conference. American Institute of Aeronautics and Astronautics. Yamashita, H., Obayashi, S., Kusunose, K., 2010. Reduction of drag penalty by means of plain flaps in the boomless busemann biplane. Int. J. Emerg. Multidiscip. Fluid Sci. 1 (2), 141–164. Yamashita, H., Yonezawa, M., Obayashi, S., Kusunose, K., 2013. A study of busemann-type biplane for avoiding choked flow. In: Collection of Technical Papers - 45th AIAA Aerospace Sciences Meeting, Vol .12, pp. 8514–8533. Yamazaki, W., Kusunose, K., 2014. Biplane-wing/Twin-Body-Fuselage configuration for innovative supersonic transport. J. Aircr. 51 (6), 1942–1952. Yamazaki, W., Kusunose, K., 2016. Aerodynamic/Sonic Boom Performance evaluation of innovative supersonic transport configurations. J. Aircr. 53 (4), 942–950. Yang, Z., Mavriplis, D.J., 2007. Higher-order time integration schemes for aeroelastic applications on unstructured meshes. AIAA J. 45 (1), 138–150. Yang, T.Y., Sung, S.H., 1977. Finite-element panel flutter in three-dimensional supersonic unsteady potential flow. AIAA J. 15 (12), 1677–1683. Yates, E.C., 1988. AGARD Standard Aeroelastic Configurations for Dynamic Response I-Wing4456. AGARD Report 765, North Atlantic Treaty Organization, Group for Aerospace Research and Development. Zhang, W., Jiang, Y., Ye, Z., 2007. Two better loosely coupled solution algorithms of CFD based aeroelastic simulation. Eng. Appl. Comput. Fluid Mech. 1 (4), 253–262. Zhang, W.-W., Ye, Z.-Y., Zhang, C.-A., Liu, F., 2009. Supersonic flutter analysis based on a local piston theory. AIAA J. 47 (10), 2321–2328. Zhao, C.X., Hua, R.H., Ye, Z.Y., Jiang, Y.W., 2015. Unsteady aerodynamic characteristics of the pitched supersonic biplane. Appl. Mech. Mater. 798, 523–530. Zhou, L., Chen, Y., Chen, F., 2013. Chaotic motions of a two-dimensional airfoil with cubic nonlinearity in supersonic flow. Aerosp. Sci. Technol. 25 (1), 138–144.
Contents lists available at ScienceDirect
Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs
Numerical analysis on flutter of Busemann-type supersonic biplane airfoil ∗
Hao Zhou, Gang Wang , Zhikan Liu School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, PR China School of Aeronautics, National Key Laboratory of Aerodynamic Design and Research, China
article
info
Article history: Received 31 July 2018 Received in revised form 9 June 2019 Accepted 29 October 2019 Available online xxxx Keywords: Fluid–structure interaction Busemann supersonic biplane Aeroelastic stabilities Navier–Stokes equations Flutter
a b s t r a c t The Busemann-type supersonic biplane can effectively reduce the wave drag through shock interference effect between airfoils. However, considering the elastic property of the wing structure, the vibration of the wings can cause the shock oscillation between the biplane, which may result in relative aeroelastic problems of the wing. In this research, fluid–structure interaction characteristics of the Busemann-type supersonic biplane at its design condition have been studied. A theoretical two-dimensional structure model has been established to consider the main elastic characteristics of the wing structure. Coupled with unsteady Navier–Stokes equations, the fluid–structure dynamic system of the supersonic biplane is studied through the two-way computational fluid dynamics/computational structural dynamics (CFD/CSD) coupling method. The biplane system has been simulated at its design Mach number with different nondimensional velocities. Different initial disturbance has been applied to excite the system and the effects of the position of the mass center on the system’s aeroelastic stability is also discussed. The results reveal that the stability of the airfoil in supersonic biplane system is decreased compared with that of the airfoil isolated in supersonic flow and such stability reduction effect should be given due attention in practical design. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction In 1935, Busemann (1935) proposed a supersonic biplane configuration to eliminates the wave drag by making use of the shock waves interference effect between airfoils. Later, Licher (1955) made an asymmetric design of the supersonic biplane to further increase its lift-drag ratio. However, due to many reasons, this concept was not actually applied to aircraft design at that time. In 21st century, the pursuit of supersonic commercial aircraft has led designers to refocus on this concept. Matsushima et al. (2006) and Kusunose et al. (2006) studied the aerodynamic performance of supersonic biplane through CFD method and discussed the aerodynamics efficiency decrease in off-design region. Choking and hysteresis phenomena are captured in their simulations and these phenomena were confirmed in experiments (Kuratani et al., 2007; Kuratan et al., 2009). Many efforts were made to counter against the choking problem (Yamashita et al., 2013, 2010) and to improve the off-design aerodynamic efficiency of the supersonic biplane (Hu et al., 2012). Kusunose et al. (2011) reviewed the development of the supersonic biplane, and further discussed the prospects as well as obstacles in the designing of the supersonic biplane. In recent years, Yamazaki and Kusunose (2014) and Yamazaki and Kusunose (2016) ∗ Corresponding author at: School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, PR China. E-mail addresses: [email protected] (H. Zhou), [email protected] (G. Wang), [email protected] (Z. Liu). https://doi.org/10.1016/j.jfluidstructs.2019.102788 0889-9746/© 2019 Elsevier Ltd. All rights reserved.
2
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
studied a biplane-wing/twin-body-fuselage configuration for advanced supersonic transport and evaluated its sonic boom performance. Ban et al. (2017) optimized this configuration through Kringing model and genetic algorithm. In the previous work, most researchers concentrate on the steady aerodynamic performance analysis of the rigid biplane configurations without considering the elastic property of wing structures, and very few of them studied the unsteady aerodynamic characteristics of the supersonic biplane. One relative discussion is given by Zhao et al. (2015) where the unsteady aerodynamic characteristic of the pitched supersonic biplane has been investigated and complex unsteady flow behaviors between the biplane have been revealed. When considering the elastic characteristics of the wing structure in unsteady flow, the vibration of wing structure causes the change of biplane configuration, which leads to the deviation of the designed shock interference and shock wave oscillation between the biplane. The shock wave oscillation can generate severe unsteady aerodynamic loads on the biplane and lead to complex aeroelastic problem. When the dynamic pressure is large enough, the wings in the supersonic biplane system may flutter, which even causes the failure of the wing structure. Such dangerous aeroelastic responses must be avoided in real supersonic aircraft design process. However, so far, the studies of the fluid–structure interaction properties of the supersonic biplane configuration are very limited according to published literatures and this will be the main topic of this paper. The most primary feature of the supersonic biplane system is the complex shock-boundary layer interactions and shock reflections between airfoils of the biplane. For the general supersonic aeroelastic analysis many simplified unsteady aerodynamic models can be used for solving unsteady aerodynamic loads, such as the piston theories (Zhang et al., 2009; Dowell and Bliss, 2013; Meijer and Dala, 2015) and the potential flow theories (Tseng and Merino, 1976; Yang and Sung, 1977). However, these efficient methods are not suited for analyzing the unsteady flow which involves complex shock wave motions and reflections. On the other hand, the unsteady shock motion and shock-boundary layer interactions are quite common in transonic flow (Iovnovich and Raveh, 2012; Fukushima and Kawai, 2018) and some supersonic internal flow (Nishizawa et al., 2012), and mostly accompanied with complex fluid–structure interaction phenomena (Raveh and Dowell, 2014; Poirel et al., 2018; Shinde et al., 2018). These studies involving unsteady shock boundary-layer interactions usually require solving the unsteady Navier–Stokes equations. While most researches about the supersonic biplane are based on inviscid flow, the unsteady shock interaction between the supersonic biplane should be simulated in viscous flow to involve the mechanism of more complex phenomena. The purpose of this paper is to provide a preliminary analysis of the fluid–structure interaction properties of baseline Busemann supersonic biplane and to identify the mechanisms of this aeroelastic system. In this paper, a theoretical 2D supersonic biplane aeroelastic model is established to consider the main elastic characteristics of the wings in biplane system. A loosely two-way CFD/CSD coupling procedure is used for the fluid–structure interaction simulation. The fully turbulence unsteady Reynolds average Navier–Stokes (RANS) equations are solved to include the shock boundary-layer interactions between the airfoils of the biplane. Different types of initial disturbance are applied to excite the system and stability analysis is applied to the airfoil in supersonic biplane system to reveal the effects of the shock interactions between the biplane on the wing’s aeroelastic stability. The effects of the position of the mass center on the system’s aeroelastic stability are also discussed. 2. Flow equations In this paper, we use an in-house CFD solver HUNS3D (Liu et al., 2018) to solve the fully turbulent flow. Spalart– Allmaras (SA) turbulence model is used in this study. The fluid governing equations in the Arbitrary Lagrangian Eulerian (ALE) form for a bounded control volume Ω with boundary ∂ Ω , are expressed in integral form as below:
∂ ∂t
∫∫∫ Ω
U dΩ +
∫∫ ∂Ω
[
F (U , V grid ) − G(U ) · ndS =
]
∫∫∫ Ω
H dΩ
(1)
where U = {ρ, ρ u, ρv, ρw, ρ E , ρ˜ ν}T and F , G are respectively represent convective and viscous fluxes terms. H is the source term, and V grid denotes the grid velocity. The equations are discretized on unstructured hybrid meshes with the cell-centered finite volume method, hence, for each mesh cells Ωi , Eq. (1) can be rewritten as: d(Ωi U i ) dt
+
∑ [
F (U im , V grid ) − G(U im ) S im = Ωi H i
]
(2)
m∈N(i)
The convective fluxes are discretized using Roe’s (1981) method, and the viscous fluxes term is obtained by reconstructed central scheme. For upwind schemes the second order accuracy is achieved via reconstructing the solution following Frink’s interpolation method (Frink, 2008). When the discrete fluxes of each integration surface are obtained, the semi-discretized form of the governing equations are then integrated in time, Eq. (2) is rewritten as: d(Ωi U i )
+ R i (U i ) = 0 (3) dt Dual time stepping method has been applied for solving Eq. (3). The time derivative is discretized by a second order scheme given by: 3(Ωi U in+1 ) − 4(Ωi U ni ) + (Ωi U ni −1 ) 2∆t
+ R i (U ni +1 ) = 0
(4)
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
3
By introducing a pseudo-time veritable t ∗ into Eq. (4), it can be expressed as Eq. (5), which could be solved through many steady flow accelerating methods implemented in HUNS3D (Wang et al., 2012; Liu et al., 2018). dU ni +1
+
dt ∗
1
Ωin+1
R ∗i (U ni +1 ) = 0
(5)
Since the RBF mesh deformation technique (Wang et al., 2018) has been applied to update the grids every step of the simulation, V grid must be calculated in line with the time integration scheme by Eq. (6) to satisfy the Geometric Conservation Law (GCL) where the Vi is control volume, first introduced by Thomas and Lombard (2008). The detailed technique for solving Eq. (6) can be found in Blazek (2005).
∑
3Ωin+1 − 4Ωin + Ωin−1
m
2∆t
(V grid · nS)nm+1 =
(6)
3. Structural dynamics model In the field of aeroelasticity, the classical two-degree-of-freedom (2DOF) model (Isogai, 1979) is widely used in many aeroelastic studies. In this model, the primary elastic characteristics of wing are modeled as pitch and plunge of the airfoil. Although, as early as in the 1950s, some literatures (Biot, 1956; Glasgow and Cronin, 1960) has pointed out that the chordwise mode may be important for the low aspect-ratio thin wings in supersonic aeroelastic analysis, there are still many studies about supersonic aeroelasticity based on 2DOF model after that, such as supersonic compressor rotor flutter (Bendiksen, 1986), supersonic flutter active control and suppression (Marzocca et al., 2002; Rao et al., 2006) as well as some studies about nonlinear flutter phenomena in supersonic flow (Firouz-Abadi et al., 2013; Zhou et al., 2013; Marzocca et al., 2005). In spite of its simplicity, this simplified model should give acceptable predictions of the main dynamic behavior of some structures in supersonic flow to some extent. Therefore, concerning our main thrust and intentions are to explore the basic aeroelastic properties of this novel configuration, it is appropriate to start with the 2DOF model for the preliminary analysis within certain assumptions. We focus on the most primary elastic properties of wing structure. Here, a typical wing section is extracted from the three-dimensional (3D) Busemann wing model, as shown in Fig. 1(a). Since the purpose in this paper is to investigate the theoretical aeroelastic property and mechanism of this configuration at certain Mach region, the structure supporting connection of the two wings at wing root is neglected for simplification and no linkages exist between the upper and lower wings. The 3D effect is not considered at this stage. It should be noticed that when the aspect ratio is low enough, the 3D effects should be considered, and the participation of chordwise mode needs to be studied in future research for practical design. The two-dimensional (2D) model has more theoretical meaning for exploring the basic aeroelastic properties and mechanism. With all these assumptions and intentions, each airfoil of the wing section can be considered as a typical 2DOF airfoil structure model in theoretical aeroelasticity studies and the upper airfoil and lower airfoil interacts only through aerodynamic forces. The structural parameters were originally defined by Isogai (1979). A sketch of this model is shown in Fig. 1(b), where h denotes the displacement away from balance position in z direction and α denotes the angular displacement. E and G represent the elastic axis and mass center respectively. b denotes the half chord length of airfoil and ab denotes the distance between midpoint of the chord and elastic axis. When a is negative, the elastic axis is at the front of the midpoint of the chord to simulate the sweep effect. When a is 0, there is no sweepback for the wing. The dimensional dynamic equations of the system can be easily obtained through the balance of force in z direction and the moment balance around elastic axis, as shown in Eq. (7). The Sα and Iα are defined as mxα b and mrα2 b2 respectively in unit span of the airfoil and m is the mass of airfoil. The Kh and Kα are the elastic coefficient of the springs. The L and Mea are the aerodynamic lift and moment around the elastic axis.
⎧ ∂ 2h ∂ 2α ⎪ ⎪ ⎨ m 2 + Sα 2 + Kh h = −L ∂t ∂t 2 2 ⎪ ∂ h ∂ α ⎪ ⎩Sα + Iα 2 + Kα α = Mea 2 ∂t ∂t
(7)
The equations can be nondimensionalized by the method mentioned in Alonso and Jameson (1994) and Liu et al. (2001), and the time in structure system is also nondimensionalized by multiplying a reference frequency ωα which is also defined as the natural frequency of rotational degree of freedom. The final nondimensional form of Eq. (7) can be written in matrix form as: M
∂ 2q + Kq = F ∂τs2
(8)
where:
[ M =
1
xα
xα
rα2
]
,K =
[
(ωh /ωα ) 0
1 0 −Cl h/b , F = V ∗2 ,q = 2Cm α rα2 π
]
[
]
[
] (9)
4
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 1. 2D airfoil structural model.
V∗ =
V
m
√ ,µ = ωα b µ πρ∞ b2
(10)
The ωh and ωα are the natural frequency of pitch and plunge defined by ωα2 = Kα /Iα , ωh2 = Kh /m. µ is defined as mass ratio. V ∗ is the nondimensional velocity also known as flutter index. Since the time in structure dynamic equation has been nondimensionalized by ωα , the structure dynamic time step ∆τs must match with the flow nondimensional time step ∆τ through Eq. (11) when solving these equations in time domain where Ma is the Mach number.
∆τs =
2Ma · ∆τ V∗
√
(11)
µ
For the biplane system, shown in Fig. 2, two structurally unconnected 2D airfoil model can be simply posed together as follow:
⎤ ∂ 2 qup [ ⎢ ∂τs2 ⎥ ⎥ + K up ·⎢ ⎣ ∂ 2q ⎦ 0 low ∂τs2 ⎡
[
M up 0
0 M low
]
0 K low
] [ ] [ ] q F up · up = qlow
F low
(12)
Noted as: Ms
∂ 2 qs + K s qs = F s ∂τs2
Which is the governing equation of the supersonic biplane FSI system.
(13)
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
5
Fig. 2. Busemann supersonic biplane. Table 1 Structure parameters of NACA64A010 at transonic region.
µ
ωh /ωα
rα2
xα
a
α0 (rad)
60
1
3.48
1.8
−2.0
0
Table 2 Structure parameters of NACA64A010 at supersonic region.
µ
ωh /ωα
rα2
xα
a
α0 (rad)
75
0.5
0.75
0.25
0
0.0873
4. FSI Coupling method By introducing state variables Q = [qTup ,qTlow ,∂ qTup /∂τs ,∂ qTlow /∂τs ] = [qTs ,∂ qTs /∂τs ], Eq. (13) can be written in state space form: ∂Q = A · Q + B · F s (Q , τs ) (14)
∂τs
Where
[ A=
O 1 −M − s Ks
I 1 −M − s Gs
]
[ B=
O 1 M− s
] (15)
A 4th order semi-implicit linear multi-step scheme (Zhang et al., 2007) has been used for solving Eq. (14). The structural part is solved implicitly, the general aerodynamic part is solved explicitly. Its accuracy and stability are proved better than classic Adams predictor–corrector procedure in Zhang et al. (2007). The specific marching steps are as follow:
⎧ F sn = F s (Q n , τsn ) ⎪ ⎪ ⎪ ⎪ f n = A · Q n + B · F sn ⎪ ⎪ ⎪ ⎪ ) ∆τs ( ⎨ Q n+1 = Q n + 55f n − 59f n−1 + 37f n−2 − 9f n−3 24 ⎪ ⎪ ⎪ f n+1 = A · Q n+1 + B · F sn+1 ⎪ ⎪ ⎪ ⎪ ) ∆τs ( ⎪ ⎩Q 9f n+1 + 19f n − 5f n−1 + f n−2 n+1 = Q n +
(16)
24 When Q are solved at each time step, the new wall boundary condition can be calculated according to. Then the mesh is updated by the mesh deformation methods (Wang et al., 2018) for the flow field computation at the next time step. 5. FSI system verification To validate the FSI procedure in this research, a standard transonic 2D flutter case of NACA64A010 (Isogai, 1979) has been studied. The structure parameters are listed in Table 1. This case is studied extensively for Euler equation while the results for Reynolds-Averaged N–S solution are much fewer. Some of the results can be found in Bohbot et al. (2001), Prananta et al. (1998) and Marti and Liu (2017). The N–S flutter boundary for Mach number 0.75–0.95 calculated in this paper is shown Fig. 3, which agrees well with others’ calculations. The multiple equilibrium points discussed in Marti and Liu (2017) are also found at Mach 0.875, as shown in Fig. 4. All these results prove that the FSI procedure implemented in this paper is correct. Since the study in present paper is in supersonic regime, the supersonic unsteady aerodynamic force also need to be validated. However, very few supersonic calculations of airfoil flutter are carried out in viscous flow directly. Here, an
6
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 3. Flutter boundary of NACA64A010 airfoil model.
Fig. 4. Heaving and pitching responses for RANS results of NACA64A010 at Ma = 0.875, Re = 1.0e7.
inviscid 2-D case of NACA64010 airfoil flutter in supersonic Mach region (Zhang et al., 2009) and a standard 3-D flutter case about AGARD445.6 wing (Yates, 1988) in viscous flow, which includes a supersonic Mach condition, is used to validate the calculation in supersonic Mach region. In the 2-D case we compare our N–S results with the inviscid results of Zhang et al. (2009) in Fig. 5(a). The structure parameters are listed in Table 2. The results agree well with each other since the viscous effect in this Mach region is not very obvious, this verifies the aeroelastic model in supersonic Mach region. In the 3-D case, the same unsteady CFD solver is employed and the CFD codes are loosely coupled with an open source CSD solver CALCULIX (Dhondt, 2018). The detailed coupling procedure is listed in Appendix and the date transformation method can refer to Mian et al. (2014). The flutter boundaries obtained by presented codes as well as experiment results and others’ calculations (Chen and Liu, 1985; Lee-Rausch and Batina, 1993; Cavagna et al., 2007; Yang and Mavriplis, 2007; Liu et al., 2001; Im et al., 2012) are shown in Fig. 5(b). The flutter boundary at supersonic Mach condition presented in this paper shows better agreement with the experiment data compared with most of others’ simulation results. These cases should validate the unsteady supersonic aerodynamic force in this study. More supersonic unsteady aerodynamic validation cases can be found in our previous work (Wang et al., 2015). 6. Flutter analysis of supersonic biplane 6.1. Structure parameters and computing mesh In this section, the developed FSI analysis system is applied to simulate the fluid–structure interactions of Busemann supersonic biplane section model. The thickness-chord ratio (t/c) of the Busemann biplane is selected as 0.10, with the gap
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
7
Fig. 5. Validation cases involving supersonic Mach region. Table 3 Structure parameters of Busemann biplane (each airfoil).
µ
ωh /ωα
rα2
xα
a
α0 (rad)
75/2
0.5
0.75
0.25 and 0
0
0
Table 4 Grid convergence. Cell number
CL_low
CM_low
CL_up
CM_up
43 011 57 055 102 535
−1.7111e−1 −1.7110e−1 −1.7110e−1
−2.8229e−3 −2.8210e−3 −2.8205e−3
1.7109e−1 1.7108e−1 1.7107e−1
2.8280e−3 2.8255e−3 2.8248e−3
of half chord length (z/c = 0.5). Minimum drag flow condition was obtained at Mach 1.7 of zero angle of attack (Kusunose et al., 2011), which will be the free-steam condition of this study. The Reynolds number are set to 10 million. The structure parameters of a 2D supersonic flutter case for NACA64A010 airfoil (Zhang et al., 2009) are used as a reference for each airfoil of the biplane. The structure parameters of the supersonic biplane are listed in Table 3. With a = 0, the elastic axis is at the midpoint of the chord, hence the swept effect is not considered in this paper. To consider the effect of sweepback further research work is needed. Noting that, the mass ratio is reduced by half for each airfoil to make sure the structure mass remains the same compared with traditional supersonic single wing. In the real engineering practice, it is much easier to adjust the mass center than the elastic axis, then two different mass center position (xα = 0.25 and xα = 0) are tested to explore its effect to system stability. An O-type hybrid grid shown in Fig. 6 is used for the simulation. First layer height of the grid in the normal direction corresponds to maximal y+ <= 1.0. The convergence of the grid has been studied in Table 4. Considering the computing cost and the force coefficients convergence, the medium size grid is chosen for the simulation. 6.2. Time domain FSI analysis The typical method to determine the stability of an aeroelasticity system is applying a small disturbance to the balanced system and monitoring the response in time domain. However, the aerodynamic loads on each airfoil of the biplane are not symmetrical, making the airfoils unbalanced at zero angle of attack. The static aeroelasticity problem is not considered in this study, however it is an important research content, here the system is balanced at zero angle of attack by external forces to exclude the static aeroelastic effects. A steady simulation has been performed at first to obtain the initial axial force and momentum around the elastic axis of each airfoil of the biplane. When starting the unsteady simulation, the opposite external forces and moments are applied to the upper and the lower airfoil respectively to cancel its initial steady aerodynamic loads, thereby ensuring that the system’s static aeroelasticity balance is excluded in this research. For the unsteady simulation in this paper, about 80 to 100 time-steps per flutter period gives good convergence of the time domain response. 6.2.1. Time domain analysis for xα = 0.25 This parameter is chosen to make the flutter easier to be excited and observed according to the experience of traditional airfoil flutter. As for stability, the key problem is to determine the critical V* of the biplane system at given condition. Two
8
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 6. Typical grid used for FSI analysis.
tentative V* are tested to give the boundary of V*. The time domain responses of two different V* (0.5, 0.7) with the same initial condition E0 = [0 0 0 0 0 0.01 0 0.01] are shown in Fig. 7. For this initial condition, the angle velocity disturbance has the same magnitude and direction. It can be seen that as V* grows, the system goes from stable to divergent, which means that the critical V* is between 0.5 and 0.7. For the stable response solution at V* = 0.5 in Fig. 7(a), (b) the phase of pitch and plunge response for upper and lower airfoil remains almost the same at entire time history and the energy of structure system is dissipated into the flow every cycle of response. In the unstable response history in Fig. 7(c), (d), some kind of phase transition occurred. When the system becomes fully divergent, the upper and lower airfoil has the nearly opposite phase and the structure system absorbs energy from the flow causing the h/b increases to an unacceptable level. The cases of same velocity index are also tested at different initial condition E0 = [0 0 0 0 0 0.01 0 −0.01], shown in Fig. 8. For this initial condition, the angular velocity disturbance has the same magnitude but in opposite direction, which makes the time response has opposite phase at the beginning. In Fig. 8(a), the similar phase transition occurs in the response of α . This implies that when the system is stable the response of α may tends to be in phase, and when the system is unstable the response of seems tends to be in opposite phase. The phase transition is not found in the calculated time series of h/b partly due to h/b is not excited in the initial condition and the limited simulation time. For the results of V* = 0.7 shown in Fig. 8(c), (d), the applied initial condition accelerated the divergence of the biplane system which also exemplifies our point. To determine the critical speed more precisely, the V* section of (0.5, 0.7) needs to be divided into more couple of small sections. The final critical V* is found near 0.61, The time domain solution for V* = 0.65, 0.625 and 0.61 with different initial condition are presented in Figs. 9–14 to indicate the critical state of the system. For the V* = 0.65 shown in Fig. 9, the amplitude of decreases gradually at first and then increase to divergence accompanied by a very clear phase transition phenomenon. The amplitude of h/b increase at entire time history. When applying the antiphase disturbance, shown in Fig. 10, both α and h/b increase to divergence directly. Noticing that when the opposite phase disturbance is applied, the α and h/b have the same trend of convergence or divergence, and when the in-phase disturbance is applied, some component of the system may have a convergent trend in the first period of time but will become divergent in the end due to the phase transition process mentioned previously. These phenomena imply that the initial condition does not affect the system’ stability however antiphase disturbance is recommended for time saving when to determine the system’s stability boundaries. In Fig. 11(a), the for V* = 0.625 has obvious convergence pattern in the simulated time period, but the system is not stable since the time response of h/b shown in Fig. 11(b) is divergent. The system’s stability can be seen more clearly when the antiphase disturbance is applied, shown in Fig. 12, which is near the critical state but slightly divergent. It should be noticed that the phenomenon of beats can be clearly observed in Figs. 11 and 13. The beats phenomenon usually occurs when the flutter is about to happen where the frequency of two mode is close to each other. The beats phenomenon proved that these two conditions are close to the critical flutter point. Figs. 13 and 14 actually draws the lower boundary of the critical V*. The more precise critical V* can be obtained by the damping interpolation method mentioned in Marti and Liu (2017). Limited to our purpose here, the actual flutter boundaries are not precisely calculated in this paper and the complete flutter boundaries of the supersonic biplane require a lot of computations which will be carried out in our future work. The main difference between supersonic biplane and its isolated airfoils is the shock interactions and reflections between the biplane, as shown in Fig. 15, which may make the aeroelastic system more complicated than the isolated airfoil. Here, the same stability analysis is also applied to the isolated upper airfoil of the supersonic biplane to identify the effects of the shock interactions and reflections between the biplane to the airfoil’s aeroelastic stability. Since the structure
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
9
Fig. 7. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0.25, E0 = [000000.0100.01].
is symmetrical, the upper airfoil is selected. For the upper airfoil of the supersonic biplane isolated in the supersonic flow, the critical V* is about 0.83, shown in Fig. 16, nearly 27% higher than that of supersonic biplane at same flow condition, indicating that the aeroelasticity stabilities of the airfoil in supersonic biplane system been greatly decreased due to the shock interactions and reflections between the airfoils. 6.2.2. Time domain analysis for xα = 0 For xα = 0, the mass center is on the elastic axis, and M becomes a diagonal matrix. The pitch and plunge are no longer coupled in structure dynamics. According to the discussion in the previous section, it is easier to determine the system’s stability by antiphase disturbance, hence all the analysis about xα = 0 is carried out under antiphase disturbance (E0 = [0 0 0 0 0 0.01 0 −0.01]). Starting the simulation from V* = 0.61, the time domain results are shown in Fig. 17 which is quite different from Fig. 14. In Fig. 17, the time domain response can be divided into two sections, the first section is before τ = 200 where the α and h/b respond around zero. In the first section, the magnitude of α increases gradually and then reaches a nearly fixed value for a short period of time. After that, the α decreases a little and then increases rapidly, goes into the next section and the response of h/b is similarly. In first section, the α , h/b of each airfoil of the biplane responses in the opposite phase respectively, the system is actually divergent. In the second section, the α and h/b deviate from zero and respond around a new position in quite different form. The α in the second section is divergent while the h/b is convergent and each airfoil of the biplane responses in the same phase. It implies that the initial equilibrium point is no longer stable and the flow pattern transformation occurs. Although the response in second section may actually not happen in real structure due the large bending deformation, these high amplitude oscillations are shown to give a more complete and intuitive feeling of the whole divergent process, and the response of ideal structure model is an interesting topic which may help to understand the mechanism of such aeroelastic system when the shock interactions are involved.
10
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 8. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0.25, E0 = [0 0 0 0 0 0.01 0 −0.01].
Fig. 9. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.65.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
11
Fig. 10. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 −0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.65.
Fig. 11. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.625.
Fig. 12. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 −0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.625.
12
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 13. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.61.
Fig. 14. RANS results of Busemann biplane for in-phase disturbance, E0 = [0 0 0 0 0 0.01 0 −0.01], Ma = 1.7, Re = 1.0e7, xα = 0.25, V* = 0.61.
Fig. 15. Flow pattern of supersonic biplane and its isolated airfoil at Ma1.7, Re10e7.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
13
Fig. 16. Heaving and pitching for RANS results of isolated upper airfoil of supersonic biplane at critical condition V* = 0.83, Ma = 1.7, Re = 1.0e7, xα = 0.25, E0 = [0 0 0 0.01].
Fig. 17. Heaving and pitching for RANs results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0.0 for different initial condition, V* = 0.61, E0 = [0 0 0 0 0 0.01 0 −0.01].
Here we extract a cycle of flow field form these two sections respectively. For the first section the snapshot is extracted from τ = 135 and for the second section, it is extracted from τ = 284. The results are shown in Figs. 18 and 19. In Fig. 18, the flow pattern is not stable, it changes as the airfoils oscillate. It can be seen clearly that a normal shock occurs near the throat of the biplane and moves forward against the flow. As the normal shock moves forward it becomes stronger. When it reaches somewhere ahead of the throat but still holds in the channel constructed by the biplane, it moves backward and becomes weaker and weaker until disappear and then another cycle starts. The severe shock oscillation causes strong pressure fluctuation which aggravated the structure instability and the unstable flow pattern finally turns into a choked flow. Fig. 19 shows this choked flow pattern for the biplane, which is much more stable than that of last section since the flow is nearly not change in its whole oscillation cycle. However, due to the strong shock ahead of the airfoils the shock interactions between the drag reduction effects of supersonic biplane is no longer exist. By decreasing V*, the final critical V* is found between 0.25 and 0.30. The time domain results of V* = 0.25 and V* = 0.3 as well as another divergent state V* = 0.4 are shown in Figs. 20–22. In these time domain results, the upper and lower airfoil are not responding symmetrically, the mechanism of such phenomena needs further research. Noticing that the critical V* for xα = 0.25 is about 0.61, the results here indicating that the stability of the system is greatly reduced as the mass center moving to the elastic axis. This is inconsistent with the general experience of traditional airfoil flutter that the critical flutter index usually increases when reducing the distance between mass center and elastic axis of the wing section. The same simulations are also performed on the isolated upper airfoil of the biplane and its critical V* is found near 1.4, shown in Fig. 23, which is higher than that of xα = 0.25 (shown in Fig. 16) and agree with our general experience
14
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 18. Flow field of Busemann biplane in an oscillation cycle at τ = 135 ∼ 142, Ma = 1.7, Re = 1.0e7, xα = 0, V* = 0.61.
on wing flutter. What is more, it shows that the shock interaction between the biplane destabilizes the system much severely when the mass center is moved to the elastic axis.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 19. Flow field of Busemann biplane in an oscillation cycle at τ = 284 ∼ 291, Ma = 1.7, Re = 1.0e7, xα = 0, V* = 0.61.
15
16
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
Fig. 20. Heaving and pitching for RANs results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0 for different initial condition, V* = 0.25, E0 = [0 0 0 0 0 0.01 0 −0.01].
Fig. 21. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0 for different initial condition, V* = 0.30, E0 = [0 0 0 0 0 0.01 0 −0.01].
Fig. 22. Heaving and pitching for RANS results of Busemann biplane at Ma = 1.7, Re = 1.0e7, xα = 0 for different initial condition, V* = 0.40, E0 = [0 0 0 0 0 0.01 0 −0.01].
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
17
Fig. 23. Heaving and pitching for RANS results of isolated upper airfoil of supersonic biplane Ma = 1.7, Re = 1.0e7, xα = 0, E0 = [0 0 0 0.01], V* = 1.4.
Fig. 24. HUNS3D/Calculix coupling procedure.
7. Conclusions In this paper, the fluid–structure interaction characteristics of the Busemann supersonic biplane has been studied. A theoretical two-dimensional biplane structure model is established to consider the main elastic property of the structure. The structure dynamic equations coupled with fully turbulent URANS equations are solved through a loosely coupled CFD/CSD method. At the designed flight condition of the supersonic biplane, the time domain responses of the structure at different flutter indexes with different initial disturbance are studied. The shock interaction between the biplane is critical for drag reduction, however, the time domain simulations reveal that the aeroelastic stability of airfoil in the supersonic biplane system can be greatly decreased compared with that of airfoil isolated in the supersonic flow due to the shock interactions and reflections between the biplane. When the disturbance is applied, the shock wave deviates from its original interference point and reflected between the biplane. The reflection point moves as the airfoils vibrate, resulting in severe change in pressure distribution on the airfoils. Such unsteady load is much severe than that on the isolated supersonic airfoil since the shock position of the isolated airfoil in unsteady supersonic flow is nearly unchanged, which should be the main reason for the decrease of the system stability. Different initial disturbance is applied to excite the dynamic system and we find that the initial disturbance affects the
18
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
system’s time domain response history but does not have much effects on the final divergence pattern. For the divergent state, α and h/b of the upper and lower airfoil tends to respond in opposite phase and the anti-phase disturbance can accelerate the divergent process of the system. The position of mass center has great effects on the system’s aeroelastic stability. The critical flutter index of the airfoil in supersonic biplane system is greatly decreased when the mass center is moved forward to its elastic axis. This result contradicts the general engineering experience that moving the mass center forward can increase the critical flutter speed of the wing, and our calculations on the isolated upper airfoil agrees with this experience. The contradictory reveals the effects of the unsteady shock wave on the airfoil’s aeroelastic stability and implies that some kind of different flutter mechanism may exist in this aeroelastic system, which requires further study. The results in the presented work show the theoretical aeroelastic properties in 2D model of supersonic biplane system, some large deformation in divergent state may actually not happen in real structure due to the structural limits, these results are shown to provide a more complete and intuitive feeling of the whole divergent process. Limited by the current research intention and content, the 3D effect as well as the chordwise mode is not considered at this stage and it is important to reinvestigate and confirm these theoretical properties through 3D cases in future work once the real structure parameters are obtained from structure engineering. Acknowledgments This study was sponsored by the seed Foundation of Innovation and Creation for Graduate Students in Northwestern Polytechnical University, PR China (No. ZZ2019001). And was also co-supported by the National Natural Science Foundation of China (No. 11772265) and the ‘‘111’’ Project of China (No. B17037). Appendix. Coupling procedure of HUNS3D/Calculix codes The coupling procedure is performed in a staggered form (see Fig. 24). The structure deformation is sent from Calculix to HUNS3D and the pressure distribution along the fluid–structure interface is sent back to Calculix at each time step. Since the coupling procedure is in a loosely coupled way, the time steps should be small enough to achieve both the stability and convergence condition. References Alonso, J., Jameson, A., 1994. Fully-implicit time-marching aeroelastic solutions. In: 32nd Aerospace Sciences Meeting and Exhibit. American Institute of Aeronautics and Astronautics. Ban, N., Yamazaki, W., Kusunose, K., 2017. Low-Boom/Low-Drag Design optimization of innovative supersonic transport configuration. J. Aircr. 55 (3), 1–11. Bendiksen, O., 1986. Role of shocks in transonic/supersonic compressor rotor flutter. AIAA J. 24 (7), 1179–1186. Biot, M.A., 1956. The divergence of supersonic wings including chordwise bending. J. Aeronaut. Sci. 23 (3), 237–251. Blazek, J., 2005. Computational Fluid Dynamics: Principles and Applications, second ed. Bohbot, J., Garnier, J., Toumit, S., Darracq, D., 2001. Computation of the flutter boundary of an airfoil with a parallel Navier–Stokes solver. In: 39th Aerospace Sciences Meeting and Exhibit. American Institute of Aeronautics and Astronautics. Busemann, A., 1935. Aerodynamic Lift at Supersonic Speeds, twelveth ed. Cavagna, L., Quaranta, G., Mantegazza, P., 2007. Application of Navier–Stokes simulations for aeroelastic stability assessment in transonic regime. Comput. Struct. 85 (11–14), 818–832. Chen, P.-C., Liu, D., 1985. A harmonic gradient method for unsteady supersonic flow calculations. J. Aircr. 22 (5), 371–379. Dhondt, G., 2018. Calculix. [Online]. Available: http://www.calculix.de/. Dowell, E.H., Bliss, D.B., 2013. New look at unsteady supersonic potential flow aerodynamics and piston theory. AIAA J. 51 (9), 2278–2281. Firouz-Abadi, R.D., Alavi, S.M., Salarieh, H., 2013. Analysis of non-linear aeroelastic response of a supersonic thick fin with plunging, pinching and flapping free-plays. J. Fluids Struct. 40, 163–184. Frink, N.T., 2008. Upwind scheme for solving the Euler equations on unstructured tetrahedral meshes. AIAA J. 30 (1), 70–77. Fukushima, Y., Kawai, S., 2018. Wall-modeled large-eddy simulation of transonic airfoil buffet at high Reynolds Number. AIAA J. 56 (6), 2372–2388. Glasgow, J.C., Cronin, J.L., 1960. Chordwise-bending effects on wing-vibration modes. J. Aerosp. Sci. 27 (12), 962–963. Hu, R., Jameson, A., Wang, Q., 2012. Adjoint-based aerodynamic optimization of supersonic biplane airfoils. J. Aircr. 49 (3), 802–814. Im, H., Chen, X., Zha, G., 2012. Prediction of a supersonic wing flutter boundary using a high Fidelity detached Eddy simulation. In: 50th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition. American Institute of Aeronautics and Astronautics. Iovnovich, M., Raveh, D.E., 2012. Reynolds-averaged Navier–Stokes study of the Shock-Buffet instability mechanism. AIAA J. 50 (4), 880–890. Isogai, K., 1979. On the transonic-dip mechanism of flutter of a swept back wing. AIAA J. 17 (7), 793–795. Kuratan, N., Ozak, S., Obayash, S., Ogawa, H., Matsuno, Toshihiro, Kawazoe, Takashi, 2009. Experimental and computational studies of low-speed aerodynamic performance and flow characteristics around a supersonic biplane. Trans. Japan Soc. Aeronaut. Space Sci. 89–97. Kuratani, N., Ogawa, T., Yamashita, H., Yonezawa, M., Obayashi, S., 2007. experimental and computational fluid dynamics around supersonic biplane for sonic-boom reduction. In: 13th AIAA/CEAS Aeroacoustics Conference (28th AIAA Aeroacoustics Conference), Aeroacoustics Conferences. Kusunose, K., Matsushima, K., Goto, Y., Yamashita, H., T. B. T.-A. A. S. M., Nakano, E., 2006. A fundamental study for the development of boomless supersonic transport aircraft, In: 44th AIAA Aerospace Sciences Meeting and Exhibit. Kusunose, K., Matsushima, K., Maruyama, D., 2011. Supersonic biplane A review. Prog. Aerosp. Sci. 47 (1), 53–87, Elsevier Ltd. Lee-Rausch, E.M., Batina, J., 1993. Calculation of AGARD wing 4456 flutter using Navier–Stokes aerodynamics. In: 11th Applied Aerodynamics Conference. American Institute of Aeronautics and Astronautics. Licher, R.M., 1955. Optimum Two-Dimensional Multiplanes in Supersonic Flow. Report No.SM-18688, Douglass Aircraft Co.. Liu, F., Cai, J., Zhu, Y., Tsai, H.M., Wong, A.S.F., 2001. Calculation of wing flutter by a coupled fluid-structure method. J. Aircr. 38 (2), 334–342. Liu, Y., Wang, G., Ye, Z., 2018. Dynamic mode extrapolation to improve the efficiency of dual time stepping method. J. Comput. 352, 190–212.
H. Zhou, G. Wang and Z. Liu / Journal of Fluids and Structures 92 (2020) 102788
19
Marti, F., Liu, F., 2017. Multiple equilibrium points of airfoil flutter in viscous flow. In: In 55th AIAA Aerospace Sciences Meeting. American Institute of Aeronautics and Astronautics. Marzocca, P., Librescu, L., Kim, D.-H., Lee, I., 2005. Supersonic flutter and LCO of airfoils via CFD/analytical combined approach. In: 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. American Institute of Aeronautics and Astronautics. Marzocca, P., Librescu, L., Silva, W.A., 2002. Flutter, postflutter, and control of a supersonic wing section. J. Guid. Control. Dyn. 25 (5), 962–970. Matsushima, K., Kusunose, K., Maruyama, D., Matsuzawa, T., 2006. Numerical design and assessment of a biplane as future supersonic transport. In: 25th International Congress of the Aeronautical Sciences. Meijer, M.-C., Dala, L., 2015. Generalized formulation and review of piston theory for airfoils. AIAA J. 54 (1), 17–27. Mian, H.H., Wang, G., Ye, Z.Y., 2014. Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach. J. Fluids Struct. 49, 186–201. Nishizawa, U., Kameda, M., Watanabe, Y., Yamamoto, S., 2012. Computational simulation of shock oscillation around a supersonic air-intake, In: Collection of Technical Papers - 36th AIAA Fluid Dynamics Conference, Vol. 1, pp. 406–415. Poirel, D., Goyaniuk, L., Benaissa, A., 2018. Frequency lock-in in pitch–heave stall flutter. J. Fluids Struct. 79, 14–25. Prananta, B.B., Hounjet, M.H.L., Zwaan, R.J., 1998. Two-dimensional transonic aeroelastic analysis using thin-layer navier–stokes method. J. Fluids Struct. 12 (6), 655–676. Rao, V.M., Behal, A., Marzocca, P., Rubillo, C.M., 2006. Adaptive aeroelastic vibration suppression of a supersonic airfoil with flap. Aerosp. Sci. Technol. 10 (4), 309–315. Raveh, D.E., Dowell, E.H., 2014. Aeroelastic responses of elastically suspended airfoil systems in transonic buffeting flows. AIAA J. 52 (5), 926–934. Roe, P.L., 1981. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys. 43 (2), 357–372. Shinde, V., McNamara, J.J., Gaitonde, D.V., Barnes, C.J., Visbal, M.R., 2018. Panel flutter induced by transitional shock wave boundary layer interaction. In: 2018 Fluid Dynamics Conference. Thomas, P.D., Lombard, C.K., 2008. Geometric conservation law and its application to flow computations on moving grids. AIAA J. 17 (10), 1030–1037. Tseng, K., Merino, L., 1976. Fully unsteady subsonic and supersonic potential aerodynamics of complex aircraft configurations for flutter applications. In: 17th Structures, Structural Dynamics, and Materials Conference. American Institute of Aeronautics and Astronautics. Wang, G., Chen, X., Liu, Z., 2018. Mesh deformation on 3D complex configurations using multistep radial basis functions interpolation. Chin. J. Aeronaut. 31 (4), 660–671. Wang, G., Jiang, Y., Ye, Z., 2012. An improved LU-SGS implicit scheme for high Reynolds number flow computations on hybrid unstructured mesh. Chin. J. Aeronaut. 25 (1), 33–41. Wang, G., Zeng, Z., Suo, Q., 2015. Trajectory simulation of a spinning projectile based on variable step size CFD/RBD method. In: AIAA Atmospheric Flight Mechanics Conference. American Institute of Aeronautics and Astronautics. Yamashita, H., Obayashi, S., Kusunose, K., 2010. Reduction of drag penalty by means of plain flaps in the boomless busemann biplane. Int. J. Emerg. Multidiscip. Fluid Sci. 1 (2), 141–164. Yamashita, H., Yonezawa, M., Obayashi, S., Kusunose, K., 2013. A study of busemann-type biplane for avoiding choked flow. In: Collection of Technical Papers - 45th AIAA Aerospace Sciences Meeting, Vol .12, pp. 8514–8533. Yamazaki, W., Kusunose, K., 2014. Biplane-wing/Twin-Body-Fuselage configuration for innovative supersonic transport. J. Aircr. 51 (6), 1942–1952. Yamazaki, W., Kusunose, K., 2016. Aerodynamic/Sonic Boom Performance evaluation of innovative supersonic transport configurations. J. Aircr. 53 (4), 942–950. Yang, Z., Mavriplis, D.J., 2007. Higher-order time integration schemes for aeroelastic applications on unstructured meshes. AIAA J. 45 (1), 138–150. Yang, T.Y., Sung, S.H., 1977. Finite-element panel flutter in three-dimensional supersonic unsteady potential flow. AIAA J. 15 (12), 1677–1683. Yates, E.C., 1988. AGARD Standard Aeroelastic Configurations for Dynamic Response I-Wing4456. AGARD Report 765, North Atlantic Treaty Organization, Group for Aerospace Research and Development. Zhang, W., Jiang, Y., Ye, Z., 2007. Two better loosely coupled solution algorithms of CFD based aeroelastic simulation. Eng. Appl. Comput. Fluid Mech. 1 (4), 253–262. Zhang, W.-W., Ye, Z.-Y., Zhang, C.-A., Liu, F., 2009. Supersonic flutter analysis based on a local piston theory. AIAA J. 47 (10), 2321–2328. Zhao, C.X., Hua, R.H., Ye, Z.Y., Jiang, Y.W., 2015. Unsteady aerodynamic characteristics of the pitched supersonic biplane. Appl. Mech. Mater. 798, 523–530. Zhou, L., Chen, Y., Chen, F., 2013. Chaotic motions of a two-dimensional airfoil with cubic nonlinearity in supersonic flow. Aerosp. Sci. Technol. 25 (1), 138–144.