2.1 Skin design

The role of the wing skin is to maintain the shape of the wing, so that it has good aerodynamic characteristics. After the skin is subjected to the aerodynamic action, the force is transferred to the connected fuselage wing skeleton. The force is complex, and the skin is directly in contact with the outside world. Therefore, the skin material is not only required to have high strength and plasticity, but also required to have a smooth surface and high corrosion resistance. When designing the deformed wing, it is necessary to select the appropriate layered high elastic polymer material to realise the flexible change of skin in any direction and at any angle; that is, the flexible skin can change synchronously with the shape of the wing. The skin structure is shown in Fig. 2.

Figure 2. Diagram of carbon fiber composite.

In order to predict the mechanical properties of flexible skin, a carbon fiber composite representative volume element (RVE) model with 70A polyurethane rubber as matrix and T300-3k carbon fiber bundle as reinforcement was constructed in this paper. The RVE model of carbon fiber composite material is shown in Fig. 3. A unidirectional carbon fiber with a diameter of 7 × 10^-3mm is laid along the X-axis in the centre. The material properties of the T300-3k carbon fiber bundle are shown in Table 1, and the property parameters of the 70A polyurethane rubber are shown in Table 2, where E, $\mu $ and G respectively represent the elastic modulus, principal Poisson’s ratio and shear modulus of the material.

Figure 3. Finite element model of 1/8 microstructure of composite materials.

Table 1. Elastic constants of T300-3k carbon fiber bundles

Table 2. Calibration results of 70A polyurethane parameters

In this paper, the strain range of polyurethane rubber skin is less than 60%, so the Mooney-Rivlin model was used to describe the tensile mechanical behaviour of polyurethane rubber, and the finite element model was used to calculate and analyse the tensile mechanical behaviour of polyurethane rubber. The expression of the Mooney-Rivlin model is as follows:

(1) \begin{align}W = {C_{10}}\left( {{{\overline I }_1} - 3} \right) + {C_{01}}\left( {{{\overline I }_2} - 3} \right) + {C_{11}}\left( {{{\overline I }_1} - 3} \right)\left( {{{\overline I }_2} - 3} \right) + \frac{1}{D}{\left( {J - 1} \right)^2}\end{align}

In the formula, C ₁₀, C ₀₁ and C ₁₁ are Mooney constants, which are usually obtained by fitting experimental data. ${\overline I _1}$ and ${\overline I _2}$ are strain invariants represented by main elongation ratio, D is incompressible parameter, J describes volume deformation, and J = 1 represents volume incompressibility.

The model describes the variation of strain energy of rubber under different tensile states. Parameters C ₁₀, C ₀₁ and C ₁₁ are closely related to rubber hardness. Therefore, the values of C ₁₀, C ₀₁ and C ₁₁ of rubber with different hardness should be calibrated according to the stress-strain experimental data. Parameter D is usually taken as a minimal positive value (1.e-6 in this paper) to ensure the computability of the model.

The finite element numerical model of polyurethane tensile was established according to uniaxial tensile experiment, and the C ₁₀ and C ₀₁ parameters of the Mooney-Rivlin model were calibrated by numerical simulation. The uniaxial tensile model is shown in Fig. 4.

Figure 4. Uniaxial tensile model of polyurethane rubber.

Because polyurethane rubber is a super-elastic material, based on the Mooney-Rivlin model, the material properties of polyurethane rubber in ANSYS only need to input the parameter values of C ₁₀, C ₀₁, C ₁₁ and D1. Therefore, the polyurethane rubber should be calibrated to correspond to C ₁₀, C ₀₁ and C ₁₁ according to the specific hardness. The mathematical optimisation model of parameter solving is as follows:

(2) \begin{align} & Min\;f(X) = \left| {\sigma _{\max }^i(X) - \sigma _0^i} \right|\nonumber\\ & s.t.\;{\varepsilon ^i} = \varepsilon _0^i\nonumber\\ & \sigma _0^i = \phi \left(\varepsilon _0^i\right)\nonumber\\ & K(X)U = F\nonumber\\ & X = [{{\rm{C}}_{10}}\quad {{\rm{C}}_{01}}\quad {{\rm{C}}_{11}}]^{\rm{T}}\end{align}

where, the objective function $f(X)$ is the difference between the predicted stress and the expected stress, and the expected stress is derived from the stress-strain curve obtained by the experimental test, expressed by $\sigma _0^i = \phi \left(\varepsilon _0^i\right)$ . Under a given strain, the predicted stress is driven to approximate the expected stress by adjusting the design variable.

The adaptive single objective optimisation method of the direct optimisation module in Ansys Workbench was used to solve the optimisation model. Set the parameter optimisation conditions as follows: minimise the maximum difference between the predicted value and the expected value, and the convergence target precision is 0.001.

Uniaxial tensile experiments were carried out on the basis of the Instron microtensile testing machine (as shown in Fig. 5) and the customised polyurethane film (as shown in Fig. 6).

Figure 5. Instron micro tensile testing machine.

Figure 6. Custom polyurethane film.

With a SHAO hardness of 70A polyurethane rubber, its density is 1,250kg/m [Reference Hajarian3]. The uniaxial tensile stress-strain curve obtained by the experiment is shown in Fig. 7.

According to the stress-strain curve shown in Fig. 7, Ansys Workbench platform was used to optimise and calibrate the corresponding values of C ₁₀, C ₀₁ and C ₁₁, as shown in Table 2.

According to the relationship between the engineering elastic constants and the flexibility matrix, the elastic modulus, principal Poisson’s ratio and shear modulus can be obtained, and the engineering elastic constants of the polyurethane-based carbon fiber composites can be obtained as shown in Table 3.

Table 3. Engineering constants of polyurethane based carbon fiber composites

Figure 7. Stress-strain curves of a 70A polyurethane sample.

2.2 Wing rib design

Wing chord bending is completed under the action of aerodynamic load and push rod driving force. Therefore, wing chord bending can be simplified as in-plane deformation of wing rib. The leading edge and trailing edge of the wing are provided with corrugated ribs, and the linear bending angle changes of the leading edge and trailing edge are realised by the tension and compression of the corrugated ribs driven by a motor. The mechanism is shown in Fig. 8.

Figure 8. Plane shape of wing ribs.

After the structure of the wing fin is determined, the optimal thickness of the corrugated structure is determined by analysing the influence of different ripple thicknesses on the maximum deformation and stress of the leading edge structure when the driving force is 40N, 50N and 60N.

As can be seen from Fig. 9, when the driving force is 60 and 50N, the maximum deflection angle of the leading edge structure reaches 60mm, but the maximum stress is large, and it is easy to fail under low-cycle fatigue load. When the driving force is 40N and the thickness of the ripple is 3mm, the maximum deformation of the leading edge is about 35mm and the maximum stress of the structure is small, which meets the design requirements.

Figure 9. Effect of corrugated thickness on structural stress and strain.

Based on the local bending characteristics of the corrugated structure, the bending deformation of gradient corrugated ribs at different positions of the wing is optimised. Finally, the reinforced fin, common fin, beam and girder are assembled into the aircraft wing, as shown in Fig. 10.

Figure 10. Aircraft wing.

The rib plate structure is made of aluminum alloy material through CNC machining, wire cutting and other processes. The aluminum alloy is 7075 aluminum alloy, and the performance is shown in Table 4.

Table 4. 7075 aluminum alloy material properties

In the service state of the deformed wing, the surface aerodynamic pressure changes dynamically with the difference of flight speed and flight altitude, and the aerodynamic pressure acting on the aerofoil also changes dynamically. The distribution vector diagram is shown in Fig. 11. According to the deformation area of the front and rear edges of the deformed wing and the vector diagram of aerodynamic pressure distribution, the deformation area of the front and rear edges is in the area of low aerodynamic pressure and changing direction, which can be simplified as the aerodynamic load acting on the outer profile of the rib plate.

Figure 11. Aerodynamic pressure distribution vector diagram of the wing.

A corrugated rib plate was established, and fixed supports were applied to the middle part of the structure, and driving loads were applied to the front and rear edges, as shown in Fig. 12. Finite element analysis software ANSYS was used to analyse the structure.

Figure 12. Structure drawing of deformed wing fin.

The simulation results show that the maximum deformation of the trailing edge structure is 76.59mm, corresponding to the deflection angle of 10.84°, and the maximum deformation of the leading-edge structure is 33.6mm, corresponding to the deflection angle of 9.54°. The maximum stress of the structure is 195.81MPa, which is lower than the tensile yield limit.

2.3 Drive device design

In the distributed driven deformation structure, the deformation of the deformation structure is shared by multiple distributed drivers. Its advantages are as follows: the size of the driver is reduced, and multiple drivers share the external load, which can reduce the requirements for the stiffness of the deformed structure and reduce the weight of the deformed structure. Its deformation form is more flexible and robust. When part of the actuator fails, the aircraft can still ensure the controllability during flight. The bending degree of each wing rib can be adjusted according to the actual flight conditions and the flight environment monitored by sensors, so that the wing performance is better.

The three-dimensional model of the electric push rod is shown in Fig. 13. The electric push rod is driven by a DC motor with good speed regulation performance and large starting torque, and the reducer in the gear box is used to increase the torque output. The push rod uses the positive and negative rotation of the motor to complete the action of pushing and pulling the rod. It has the absolute self-locking function, which can realise the push and pulling of one side of the lead rod separately, and can also make both sides of the lead rod work simultaneously. The stroke position of push rod is determined by photoelectric encoder.

Figure 13. Schematic diagram of driving push rod.

The force exerted by the push rod driving structure is even the force exerted by the deformation of the corrugated structure, which is similar to the spring structure. Hooke’s Law $F = kx$ can be used to calculate the size of the driving force.

The elastic force when the corrugated structure extends (or shortens) unit length, namely the structural stiffness k, is:

(3) \begin{align}k = \frac{l}{s}\frac{{E{\delta ^3}}}{{12\left( {1 - {v^2}} \right)}}\end{align}

where, l is the half-wave length, s is the half-wave arc length, E is the elastic modulus, $\delta $ is the thickness of the corrugated structure, and v is Poisson’s ratio.

Where,

(4) \begin{align}s = \int_0^l {\sqrt {1 + \frac{{{\pi ^2}{f^2}}}{{{l^2}}}{{\cos }^2}\left(\frac{\pi }{l}x\right)} } dx\end{align}

where, f is the amplitude. The schematic diagram of corrugated structure is shown in Fig. 14.

Figure 14. Corrugated structure.

The corrugated structure parameters designed in this paper are shown in Table 5.

Table 5. Parameter of corrugated structure

By calculation, the structural stiffness of corrugated structure $k = 495.58N/m$ . According to the structural deformation and deformation data before and after wing deformation in 2.2, it can be seen that in the limit state, the change of leading edge length is 17.73mm, the force required to drive leading edge deformation is 8.79N, and the torque required is 1.32Nm. The trailing edge length change is 49.20mm, the force required to drive the wing trailing edge deformation is 24.38N, and the torque required is 6.58Nm. According to the selection of suitable motor, the push rod can achieve the deformation rate of seconds.

The simulation model of the push-rod was established in ADAMS. Set the displacement and speed of the push rod in one direction, as shown in Fig. 15. The displacement and speed of the push rod at the other end are equal in magnitude and opposite in direction. At first, the speed of the putter increases rapidly from 0 to 50mm/s after accelerating, and the putter keeps a uniform motion. At this time, the putter extends at a uniform speed. When the movement reaches 0.8s, the speed of the push rod decreases to 0 at the instant, and it can be seen from the displacement curve that the push rod extends to the maximum. After another 0.2s stop, the push rod will shorten the movement to the initial state stop according to the same motion principle.

Figure 15. Displacement and velocity diagram of push rod in unidirectional motion.

Due to the low motion speed of the mechanism, the elasticity of the component has little influence on the displacement of the mechanism, which can be basically ignored. According to the speed curve of the push rod, the influence of member flexibility on the speed mainly focuses on the moment when the push rod stops completely. Before the putter stops, there will be a swing of about 5mm/s in the opposite direction.

3.1 Fluid-structure interaction analysis

Two methods of fluid-structure deinteraction calculation and fluid-structure interaction calculation were used respectively to obtain the wing variant characteristics under typical working conditions, including structural displacement response, variant aerodynamic force and transient flow field pressure distribution. By comparing the results of the two methods, the influence of fluid-structure interaction effect on the wing variant characteristics is analysed.

The flow of fluid-structure interaction calculation is shown in Fig. 16.

Figure 16. Fluid-structure interaction calculation process.

Under the conditions of inlet flow with Mach number of 0.4, Re number of 1×10⁷ and angle-of-attack of 3.19°, unsteady aerodynamic coefficients and structural displacement response of the leading and trailing edges of the aerofoil are obtained by fluid-structure deinteraction and fluid-structure interaction respectively. The direction of the deflection driving force is the aerofoil chord direction (x direction), and the driving force instruction signal is shown in Fig. 17.

Figure 17. Driving force command signal.

Figure 18 compares the surface deformation responses obtained by the two calculation methods. It can be seen that the calculation deviation of the final leading edge deflection angle obtained by the two methods is relatively small, about 1%. However, the deviation of the trailing edge of the final state is slightly larger, about 6.5%. There is a certain deviation in the aerodynamic response obtained by the two calculation methods. The deviation of the lift coefficient in the final state is about 3.7%, that of the drag coefficient in the final state is about 4.0%, and that of the torque coefficient in the final state is about 3.6%, as shown in Fig. 19(a)-(c). The pressure distribution of the final flow field obtained by the two calculation methods also has a certain deviation, as shown in Fig. 19(d).

Figure 18. Structural deformation response.

Figure 19. Pneumatic response.

The above results show that with the increase of Mach number, the influence of aerodynamic interaction on the variation characteristics gradually increases. Therefore, in order to be closer to the reality, the fluid-structure interaction calculation method is used in this paper to analyse the variant characteristics of aerofoil and wing.

For incompressible fluids, fluid density is constant, and the mass conservation equation can be written as:

(5) \begin{align}\frac{{\partial {u_i}}}{{\partial {x_i}}} = 0\end{align}

where: i is the subscript, indicating the direction of the coordinate axis; u is for fluid velocity.

For turbulence, Reynolds time-mean equation can be expressed as:

(6) \begin{align}\frac{{\partial \overline {{u_i}} }}{{\partial t}} + \frac{{\partial \overline {{u_i}} \overline {{u_j}} }}{{\partial {x_j}}} = {f_i} - \frac{1}{\rho }\frac{{\partial \bar p}}{{\partial {x_i}}} + \frac{1}{\rho }\frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial \overline {{u_i}} }}{{\partial {x_j}}} - \rho \overline {{{u^{\prime}}_i}{{u^{\prime}}_j}} } \right)\end{align}

where: j and i are both subscripts; $\bar u$ is time average flow rate; $u^{\prime}$ is pulsating velocity; t is for time; f stands for mass force; $\rho $ is fluid density; p is for pressure; $\mu $ is the dynamic viscosity of the fluid. Compared with the N-S equation, the Reynolds time mean equation has an extra Reynolds stress $ - \rho \overline {{{u^{\prime}}_i}{{u^{\prime}}_j}} $ , which is a symmetric second order tensor and related to the velocity gradient.

The flexible aerofoil with variable bending adopts homogeneous and isotropic structure, and the following equation can be used:

(7) \begin{align}{\rho _s}\frac{{{D^2}{\bf{u}}}}{{D{t^2}}} - \nabla \cdot {\bf{P}} = 0\end{align}

(8) \begin{align}{\bf{P}} = {\bf{F}} \cdot {\bf{S}}\end{align}

(9) \begin{align}{\bf{F}} = {\bf{I}} + \nabla {\bf{u}}\end{align}

(10) \begin{align}{\bf{S}} = 2\mu {\bf{E}} + \lambda tr({\bf{E}}){\bf{I}}\end{align}

(11) \begin{align}{\bf{E}} = \frac{1}{2}\left(\nabla {\bf{u}} + \nabla {{\bf{u}}^T} + \nabla {{\bf{u}}^T}\nabla {\bf{u}}\right)\end{align}

where: ${\bf{u}}$ is the displacement field; ${\rho _s}$ is the density of the structure. ${\bf{P}}$ is the Piola-Kirchhoff first stress tensor; ${\bf{S}}$ for the Piola-Kirchhoff second stress tensor; ${\bf{F}}$ stands for deformation gradient; ${\bf{E}}$ stands for Green strain tensor; ${\bf{I}}$ represents the characteristic tensor; $\mu $ and $\lambda $ are Lamet parameters.

On the interaction surface, the fluid model and the structure model must meet the kinematic and dynamic conditions for interaction. The kinematic condition means that the velocity and displacement on the interaction surface must be continuous, which can be described by the following expressions:

(12) \begin{align}{{\bf{V}}_f} = {{\bf{V}}_s}\end{align}

(13) \begin{align}{{\bf{u}}_f} = {{\bf{u}}_s}\end{align}

where, ${{\bf{V}}_f}$ , ${{\bf{u}}_f}$ , ${{\bf{V}}_s}$ and ${{\bf{u}}_s}$ are the velocity and displacement of liquid and the velocity and displacement of structure, respectively.

The dynamic condition refers to that the forces on the interaction surface must be in equilibrium, which can be described by the following expressions:

(14) \begin{equation}{\bf{n}} \cdot {{\boldsymbol{\sigma }}_f} = {\bf{n}} \cdot {{\boldsymbol{\sigma }}_s}\end{equation}

where: ${\bf{n}}$ is the normal direction; ${{\boldsymbol{\sigma }}_f}$ and ${{\boldsymbol{\sigma }}_s}$ are the stress of liquid and structure, respectively.

The fluid side and the structure side are coupled by exchanging variables on the interaction surface. At each time step, the fluid part provides pressure and viscous force to the structure, and the structure part provides displacement increment and velocity to the fluid.

3.2 Comparison of aerodynamic characteristics of flexible aerofoils with different leading edge and trailing edge angles

C-type calculation area was adopted, the left half of the calculation domain was semicircle, and the tail of the aerofoil overlapped with the centre of the semicircle. The C-shaped grid division of the whole calculation domain was shown in Fig. 20, and the left end was set as the velocity inflow boundary condition. The upper and lower are set as symmetric boundary conditions: ${{\partial u} / {\partial y}} = v = 0$ ; The right exit is set as Newman boundary condition: ${{\partial u} / {\partial x}} = {{\partial v} / {\partial x}} = 0 $ ; The aerofoil is set as the fluid-structure interaction boundary condition.

Figure 20. Aerofoil initial state aerodynamic mesh.

The grid structure is C-type. 401 grid points are uniformly arranged on the aerofoil surface, 101 grid points are arranged in the normal direction, and 51 grid points are arranged in the wake area. Elliptic partial differential equation is used to smooth the meshes and ensure the good orthogonality of the meshes near the object surface.

Nine groups of aerodynamic profiles were selected with leading edge declinations of 0°, 5° and 10° and trailing edge declinations of 0°, 2° and 6°, and aerodynamic characteristics were analysed under subsonic conditions. The incoming conditions of subsonic flow field are: α = 3.19°, Ma = 0.4, Re = 1 × 10⁷. Among them, α, Ma and Re are the angle-of-attack, Mach number and Reynolds number of the air flow. Figure 21 shows the calculation grid of each typical deflection profile.

Figure 21. Each typical deflection profile is calculated on the grid.

Figure 22 shows the transient pressure distribution under different front and rear edge angles. It can be seen that with the increase of the leading edge declination angle, the low pressure area of the upper wing surface expands and the position of the minimum pressure point moves backward. The stagnation point of the lower wing surface moves forward slightly with the increase of the declination angle of the trailing edge. The combined effect of the two causes the lift coefficient of the aerofoil to increase. From the above analysis results, it can be seen that the deformed wing designed in this paper has better aerodynamic performance than the traditional fixed rudder aerofoil.

Figure 22. Transient pressure distribution for each typical deflection profile.

As shown in Fig. 23, the pressure coefficient distribution of surface under the subsonic velocity shows, it can be seen that the aerodynamic performance can be improved by proper deflection of front and rear edges under the condition of subsonic velocity. The abscissa in the figure represents the relative chord length of each aerofoil.

Figure 23. Subsonic surface pressure coefficient.

3.3 Comparison of aerodynamic characteristics of flexible aerofoil with variable bending under different push rod action schemes

The transient pressure distributions under four typical leading edge and trailing edge deflection timing schemes are calculated under subsonic flow field conditions. Scheme 1 is to deflect the leading edge and trailing edge of the push rod at the same time with 0.2s of actuating time. Scheme 2 is to deflect the leading edge and trailing edge at the same time with 0.04s of actuating time. Scheme 3 is to deflect the trailing edge after the leading edge is completed, and the total time of deformation operation is 0.04s. In Scheme 4, the leading edge was deflected after the trailing edge was completed, and the total time of deformation operation was 0.04 s. By comparing the unsteady aerodynamic responses of the four schemes, the influence of different pushrod action schemes on the unsteady aerodynamic characteristics and unsteady delay characteristics of the variant aerodynamics is analysed.

Figure 24 shows the transient pressure distribution of Scheme 1, where Tu represents the operating time and Ts represents the moment when the flow field reaches the steady state. It can be seen that with the increase of the deflection angle of the front and rear edges, the position of the minimum pressure point moves back to the leading edge deformation operating zone. With the further increase of the deflection angle of the front and rear edges, the area of the low pressure area of the upper wing surface increases significantly, resulting in the increase of the aerofoil lift coefficient. It is noted that the deflection deformation of Scheme 1 has little effect on the pressure distribution of the lower wing surface. Figure 25 shows the transient pressure distribution of Scheme 2, and the variation rule of unsteady pressure distribution with structural deformation is similar to Scheme 1.

Figure 24. Transient pressure distribution (simultaneous action, 0.2s).

Figure 25. Transient pressure distribution (simultaneous action, 0.04s).

Figure 26 shows the transient pressure distribution of Scheme 3. In the first stage, with the increase of the leading edge declination angle, the position of the minimum pressure point moves backward, resulting in the increase of the pressure of the upper wing surface and the decrease of the aerofoil lift. In the second stage, with the increase of the trailing edge declination angle, the area of the low pressure area of the upper wing surface increased significantly, resulting in the increase of the lift coefficient.

Figure 26. Transient pressure distribution (first leading edge next trailing edge, 0.04s).

Figure 27 shows the transient pressure distribution of Scheme 4. In the first stage, with the increase of the trailing edge declination angle, the area of the low pressure area on the upper wing surface increases. In addition, the area of the high pressure area near the stagnation point of the lower wing surface increases, and the combined effect of the two increases the lift coefficient of the aerofoil. In the second stage, with the increase of the leading edge declination angle, the position of the lowest pressure point on the upper wing surface moves back, and the area of the low pressure area increases significantly, resulting in the increase of the lift coefficient.

Figure 27. Transient pressure distribution (first trailing edge next leading edge, 0.04s).

Figure 28 compares the structural deformation and aerodynamic response of the four variant schemes. It can be seen that the unsteady lift coefficient of Scheme 3 gradually decreases at the leading edge deflection stage, but gradually increases at the trailing edge deflection stage. In Scheme 4, the unsteady lift coefficient ends at the trailing edge deflection and local oscillation occurs at the beginning of leading edge deflection, while in Scheme 1 and 2, the unsteady lift coefficient is smoother. In Scheme 3, the unsteady drag coefficient decreases gradually at the leading edge deflection stage and increases gradually at the trailing edge deflection stage. In addition, the peak drag coefficient of Schemes 2, 3 and 4 is higher than that of Scheme 1. The unsteady moment coefficients of Scheme 3 and 4 show obvious local oscillation in the transition stage of the leading and trailing edge deflection, while the unsteady moment coefficients of Schemes 1 and 2 are relatively smooth. In addition, the aerodynamic delay effect of Schemes 2, 3 and 4 is higher than that of Scheme 1.

Figure 28. Structural deformation and aerodynamic response.

In conclusion, the unsteady aerodynamic force in the scheme of simultaneous front and back edge operation is relatively smooth, while the unsteady aerodynamic force in the scheme of sequential front and back edge operation at two stages occurs local oscillation in the transition process of two-stage deflection. At the same time, the rapid deformation rate will also lead to the aerodynamic oscillation of the variant and increase the aerodynamic delay effect.

4.1 Dry wind tunnel test platform

The dry wind tunnel test system is composed of excitation force source system, multi-input and multi-output force controller and real-time unsteady aerodynamic software. The concept of dry wind tunnel test system is shown in Fig. 29. The main equipment includes centralised force loading system, measurement sensor system, data acquisition and control system and test equipment:

1. (1) The concentrated force loading system consists of three groups of 100N shaker systems, each group of shaker system includes one shaker and one power amplifier.

2. (2) The measuring and sensing system includes a set of laser displacement sensor system, six acceleration sensors, three force sensors and three four-channel signal adjusters. The laser displacement sensor system consists of six laser displacement sensors and one displacement sensor controller group.

3. (3) The data acquisition and control system is the core of the excitation source test control system. As the lower computer of the excitation source test control system, the system, after receiving the command issued by the upper computer, converts it into the corresponding time sequence signal according to the command to directly control the corresponding test equipment, and collects the state data of the equipment, converts it into digital signals and then feeds back to the upper computer.

4. (4) The test kit includes one test piece and one test bed base.

Figure 29. Dry wind tunnel test system platform.

Working principle:

To solve the reduced order model of unsteady aerodynamic force, the model can calculate the unsteady aerodynamic force quickly by real-time calculation and feedback. The aerodynamic de-order model can be combined with the excitation force source system. In the excitation force source system, the aerodynamic de-order model calculates the unsteady aerodynamic force by reading the signals of the displacement sensor and the acceleration sensor, and then loads it to the shaker to generate aerodynamic feedback force during the ground vibration test. The dry wind tunnel test system can replace the traditional wind tunnel with a reduced order model of unsteady aerodynamic forces that can be calculated, and a real aircraft structure mounted on the ground vibration test system can replace the reduced scale model.

Working process:

1. (1) The shaker was placed on the base of the test bed, and the main beam of the wing with variable bending was locked and fixed by the fixture. The shaker, force sensor, acceleration sensor and laser displacement sensor were installed in conjunction with the wing with variable bending.

2. (2) After the input signal data is transmitted down through the upper computer, the signal data is converted into excitation signal, which is sent out through the lower computer, and loaded to the shaker through the gain of the power amplifier, so that it can be loaded according to the calculated force load.

3. (3) The response signals of force sensor, acceleration sensor and laser displacement sensor, such as the force of the wing with variable bending, surface acceleration of the test part and surface displacement of the test part, are respectively obtained. After processing by the displacement sensing controller and the four-channel signal adjustability instrument, the signals are sent to the lower machine for collection and processing.

4. (4) The response signal is processed and converted into data through the lower computer to the upper computer for display processing, data analysis, data output storage and other operations.

4.1 Aerodynamic load reduction model

The unsteady panel method can handle complex structures, provide convenient and accurate unsteady aerodynamic prediction, and generate the aerodynamic influence coefficient (AIC) matrix that directly relates downwash to unsteady pressure coefficient:

(15) \begin{align}\left\{ {\Delta {C_p}} \right\} = \left[ {AIC} \right]\left\{ W \right\}\end{align}

where $\Delta {C_p}$ is unsteady pressure jumps, $W$ is the downwash caused by structural oscillations.

The doublet lattice method (DLM) in the unsteady aerodynamic panel method can eliminate the wake modeling in the velocity potential method, and easily generate an AIC matrix, which relates the structural deformation to the aerodynamic force. DLM solves the so-called time-linearised small disturbance (TLSD) equation:

(16) \begin{align}\frac{\partial }{{\partial x}}({\sigma _v}) = {\phi _{xx}} + {\phi _{yy}} + {\phi _{zz}} - 2\frac{{ikV_\infty ^2}}{{{\beta ^2}}}{\phi _x} + \frac{{{k^2}V_\infty ^2}}{{{\beta ^2}}}\phi \end{align}

where $\phi $ is the perturbed unsteady velocity potential, k is the reduced frequency, ${V_\infty }$ is airspeed.

Equation (16) is obtained by linearising the structural amplitude with the nonlinear transonic small disturbance equation, assuming that the structural amplitude is small. The solution of Equation (16) varies linearly with the amplitude of structural oscillation, but includes the nonlinear transonic shock effect embedded in the steady background flow. Let the supersonic perturbed potential be represented by $\phi $ at the control point located either on the wing or the body. Thus, $\phi $ can be expressed as:

(17) \begin{align}\phi ({x_0},{y_0},{z_0}) = \sum\limits_{i = 1}^{na} {{{\bf{A}}_i}{\bf{y}}(j - i)} + \sum\limits_{i = 0}^{nb - 1} {{{\bf{B}}_i}{\bf{u}}(j - i)} \end{align}

In the formula, ${\bf{y}}$ and ${\bf{u}}$ respectively represent the column vector composed of $ny$ outputs and $nu$ inputs, $na$ and $nb$ respectively represent the delay order of output and input, and matrix ${\bf{A}}$ and ${\bf{B}}$ represent the matrix to be identified. For aerodynamic identification, ${\bf{u}}$ is the modal displacement and ${\bf{y}}$ is the calculated modal aerodynamic coefficient.

The unknown matrices ${\bf{A}}$ and ${\bf{B}}$ in the system model were estimated by the least squares (LS) method. In order to facilitate the stability analysis of aeroelasticity, the difference model of Equation (17) is transformed into a state-space model. Define the state vector ${{\bf{x}}_a}(j)$ :

(18) \begin{align}{{\bf{x}}_a}(j) = {\left[ {{\bf{y}}(j - 1), \cdots ,{\bf{y}}(j - na),{\bf{u}}(j - 1), \cdots ,{\bf{u}}(j - nb + 1)} \right]^T}\end{align}

Then the aerodynamic equation of state and output equation in discrete space can be written as:

(19) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{{\bf{x}}_a}(j + 1) = {{{\tilde{\textbf{A}}}}_a}{{\bf{x}}_a}(j) + {{{\tilde{\textbf{B}}}}_a}{\bf{\xi }}(j)}\\[5pt] {{{\bf{f}}_a}(j) = {{{\tilde{\textbf{C}}}}_a}{{\bf{x}}_a}(j) + {{{\tilde{\textbf{D}}}}_a}{\bf{\xi }}(j) + {{\bf{f}}_a}}\end{array}} \right.\end{align}

where:

(20) \begin{align}{{\tilde{\textbf{A}}}_a} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {{{\bf{A}}_1}} & {{{\bf{A}}_2}} & \cdots & {{{\bf{A}}_{na - 1}}} & {{{\bf{A}}_{na}}} & {{{\bf{B}}_1}} & {{{\bf{B}}_2}} & \cdots & {{{\bf{B}}_{nb - 2}}} & {{{\bf{B}}_{nb - 1}}}\\[4pt] {\bf{I}} & 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & 0 & 0\\[4pt] 0 & {\bf{I}} & \cdots & 0 & 0 & 0 & 0 & \cdots & 0 & 0\\[4pt] \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\[4pt] 0 & 0 & \cdots & {\bf{I}} & 0 & 0 & 0 & \cdots & 0 & 0\\[4pt] 0 & 0 & \cdots & 0 & {\bf{I}} & 0 & 0 & \cdots & 0 & 0\\[4pt] 0 & 0 & \cdots & 0 & 0 & {\bf{I}} & 0 & \cdots & 0 & 0\\[4pt] 0 & 0 & \cdots & 0 & 0 & 0 & {\bf{I}} & \cdots & 0 & 0\\[4pt] \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\[4pt] 0 & 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & {\bf{I}} & 0\end{array}} \right]\end{align}

(21) \begin{align}{{\tilde{\textbf{B}}}_a} = {\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {{{\bf{B}}_0}} & 0 & 0 & \cdots & 0 & {\bf{I}} & 0 & 0 & 0 & \cdots \end{array}} \right]^T}\end{align}

(22) \begin{align}{{\tilde{\textbf{C}}}_a} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {{{\bf{A}}_1}} & {{{\bf{A}}_2}} & \cdots & {{{\bf{A}}_{na - 1}}} & {{{\bf{A}}_{na}}} & {{{\bf{B}}_1}} & {{{\bf{B}}_2}} & \cdots & {{{\bf{B}}_{nb - 2}}} & {{{\bf{B}}_{nb - 1}}}\end{array}} \right]\end{align}

(23) \begin{align}{{\tilde{\textbf{D}}}_a} = \left[ {{{\bf{B}}_0}} \right]\end{align}

where, ${\bf{\xi }}$ is the generalised displacement of the structure, ${{\bf{f}}_a}$ is the generalised aerodynamic coefficient.

In dry wind tunnel test platform setup, the number of sensors and shakers is generally very small. As a result, the matrix size in the formula is small, which means that the formula can be very efficient in calculating the aerodynamic forces that generate the shaker position in real time.

4.2 Test results and discussion

Conventional wind tunnels are used to measure wing lift coefficient, drag coefficient and stall angle-of-attack. The experiment with wind tunnels is based on the law of similarity of motion. When the critical angle-of-attack is exceeded, the boundary layer on the upper surface of the aerofoil will be seriously separated, and the lift force will drop sharply and the aerofoil cannot maintain normal flight. This phenomenon is called stall. Reduced frequency is an important factor affecting aerofoil dynamic stall, and its expression is $k = \dfrac{{\omega c}}{{2u}}$ , where $\omega $ is the oscillation angular velocity, c is the chord length, and u is the incoming flow velocity. In this paper, three reduced frequencies were set respectively for two kinds of aerofoils before and after deformation, and the dry wind tunnel platform vibration test with different angles of attack was carried out. The stored data were processed and calculated to analyse the lift coefficient and drag coefficient of the two kinds of aerofoils at different reduced frequencies.

The frequency of dry wind tunnel chatter is less than that of real wind tunnel vibration. Figure 30 shows the variation curve of lift drag coefficient of the deformed wing at different shrinking frequencies. The variation curve of the lift drag coefficient of the deformed aerofoil is consistent with that of the deformed aerofoil as a whole except for the local obvious fluctuation. When k = 0.08, the lift coefficient keeps a low state in the process of increasing the angle-of-attack, and keeps a high state in the process of decreasing the angle-of-attack, and the drag coefficient does not change significantly in the whole process. By comparing the variation curves of the lift drag coefficient of the aerofoil before and after deformation at the same reduction frequency, it can be found that there is no significant difference between them when k is small. When k = 0.08, the lift coefficient of the deformed aerofoil is slightly higher than that of the deformed front aerofoil, the drag coefficient is lower than that of the deformed front aerofoil, and the maximum value of the drag coefficient is obviously lower than that of the deformed front aerofoil. In conclusion, at high reduced frequency, the aerodynamic characteristics of the deformed aerofoil are significantly better than those of the deformed front aerofoil, featuring high lift-drag ratio and low drag.

Figure 30. Variation curve of lift drag coefficient of deformed wing.