Implementation of the Methodology of Coupled Modelling of Aerodynamics and Strength to Determine the Onset Speed of the Bending-Torsional Flutter of the Wing Console

. Time-domain motion simulations using the averaged Navier-Stokes Euler and Reynolds equations were performed for an AGARD 445.6 wing in transonic flow. These simulations were compared with experimental results as well as with the solution in the frequency domain using a method based on the fast Fourier transform. The model matches in frequencies and shapes with experiment.


Introduction
Flutter is a phenomenon of dynamic aeroelasticity caused by the interaction of inertial, elastic and aerodynamic forces.The oscillation mode of the structure determines the possibility of a flutter and depends on many parameters, including the flight speed.In this paper, using numerical modelling methods, the interaction of several simultaneously resonating oscillation modes at different velocity pressures and Mach numbers is considered on the example of the classical experiment Agard 445.6.

Computational Dynamics Modelling
In the experiment [1], a wing console model was used with an elongation λ = 1.6525, a narrowing = 0.6576 and a sweep along the quarter chord 1/ 4 X =45º.The profile used is NACA 65A004.The model is made of laminated red wood.To reduce stiffness, holes were drilled in the model along the normal to the plane of the chord.The holes obtained as a result of drilling were filled with foam to preserve the theoretical contour.The model was loaded with the help of an incoming air flow.During the experiment, the Reynolds number varied from Re = 0.5 million to Re = 6.7 million.The wing console was fixed to the wall of the wind tunnel without imitation of the fuselage in such a way that the root part of the console was immersed in the boundary layer, in order to minimize the effect of the joint of the model with the wall on the results of the experiment [2].The model was equipped with strain gauges mounted externally at the root of the wing console and oriented in such a way as to distinguish the bending and torsion of the wing.The strain gauge signals were recorded on a multichannel oscilloscope and displayed on a cathode ray oscilloscope to help determine the flutter onset.The displacements of the wing console were additionally analysed using video recordings with a frequency of 128 frames per second.
To construct a finite element model, the original geometric characteristics were preserved, and the equality of stiffness values was ensured by the reduced properties of the material [3] To simulate the constraints, an encastre was used on the side of the wing root chord [4].The model consists of two finite element grids: a structural-power mesh, necessary for the perception of external loads [5], and an aerodynamic mesh, necessary for obtaining external aerodynamic loads and transferring them to a structural-power mesh [6].The aerodynamic mesh used first-order quadrangular elements with an average face size of 2.5 mm.The structural and power mesh used hexagonal elements of the first order with an average face size of 5 mm.
At the initial stage, the validation of the finite element model of the wing console was carried out through modal analysis.The calculation methodology implemented in the CAE program is based on the following relation [7]: where i q -displacement components, As can be seen from the above relation, the modal analysis does not take into account the parameters of the external influence.There is also some overestimation of the values of natural frequencies caused by the fact that when constructing a finite element model, all the restrictions of displacement are considered ideal, which does not fully correspond to the conditions for fixing the physical model.The values of the natural frequencies of the wing console oscillations obtained using the finite element model exceed the values of the natural frequencies obtained experimentally by no more than 5%.The calculation results and the error are presented in Table 1.Next, a CFD calculation was carried out.The determination of the pressure values on the surface of the model is carried out by numerical simulation [8].The calculations were performed by the computational gas dynamics method -the solution of the Reynoldsaveraged Navier-Stokes equations [9,10,11]: where , i j u -velocity components,  -density,  -viscosity, i f -a force per unit of volume such as gravitational force, e -is specific internal energy, h -enthalpy, j q -heat flux vector.
Based on the initial geometric characteristics of the experimental setup, a CAD model is constructed, shown in Figure 1.The calculated area is a hemisphere.Sym (the plane of symmetry of the hemisphere) and Freestream (the outer boundaries of the computational domain) are used as boundary conditions.The surface of the experimental setup is a Wall boundary condition (a wall without slippage).The external boundaries of the computational domain, the boundaries of the computational domain at the surface of the model and the resolution of the wall area are shown in Figures 2, 3   The experimental model does not have drainage holes and pressure sensors, so the measurement of the aerodynamic forces acting on the model is also not provided.Thus, it is not possible to validate the CFD solver for this experiment.Figure 4 shows the nature of the distribution of the pressure coefficient over the surface, obtained by numerical simulation (M = 1.141, angle of attack = 0º), in comparison with the results of the authors [12,13,14].
Visualization of the pressure distribution over the surface of the model is shown in Figure 5.  Next, the loads from the CFD solver are transferred to the finite element model and an implicit dynamic analysis is started.The equation for this case is given below [7]: where  -viscous damping.
In contrast to the equation used in the modal analysis, the above ratio takes into account the parameters of the external influence in the form of damping [15].The calculation duration is 0.75 s.The time integrator parameter alpha = 0.At the beginning of the step, the wing is twisted by a moment of 10000 N, acting relative to the normal to the plane of the end rib.The increment value is fixed and is 0.005 s.The load and boundary conditions applied to the model are shown in Figure 6.

Conclusion
As you can see, with an increase in the Mach number, several frequencies begin to affect the nature of the movement due to the interaction of several oscillation modes.The flutter goes from bending to bending-torsional.
Based on the results of the fast Fourier transform, it can be traced that the value of the main peak frequency increases with the growth of the Mach number.
As a result of the work, it can be concluded that with an increase in the Mach number, the number and values of frequencies that determine the oscillation mode of the console

Fig. 4 .
Fig. 4. Distribution of the pressure coefficient in the 26% span section.

Fig. 5 .
Fig. 5. Distribution of the pressure coefficient over the surface.

Fig. 6 .
Fig. 6.Load and boundary conditions in implicit dynamic analysis.

Figures 7 -
Figures 7-10 below show the dependences of sensor displacements on time at different Mach numbers.Further, to calculate the frequencies of wing oscillation during flight, calculations were carried out in the frequency domain.As a result of the work, peak frequency values were allocated for each of the 4 Mach numbers.

Fig. 7 .
Fig. 7.The dependence of sensor 1 movement on time and the results of the short-term Fourier transforms at М=0.3.

Fig. 8 .
Fig. 8.The dependence of sensor 1 movement on time and the results of the short-term Fourier transforms at М=0.5.

Fig. 9 .
Fig. 9.The dependence of sensor 1 movement on time and the results of the short-term Fourier transforms at М= 0.7.

Fig. 10 .
Fig. 10.The dependence of sensor 1 movement on time and the results of the short-term Fourier transforms at М=1.1.