2.1 Equivalent aerofoil

The steady lifting-line method uses the aerodynamic data of 2D aerofoil at a series of angles of attack to calculate the aerodynamic forces of 3D wing. Similarly, the key to the improved method in this paper lies in the utilisation of unsteady aerodynamic data of 2D aerofoil.

For the pitching motion, the literature [Reference Parenteau, Plante and Laurendeau19] used the unsteady aerodynamic data of the aerofoil as the input for 3D calculations. In this paper, we consider that every phase of the pitching motion corresponds to an aerofoil, which can be called an equivalent aerofoil. Meanwhile, the angle-of-attack of the equivalent aerofoil is equal to the mean angle-of-attack of the periodic motion.

Figure 3 gives the illustration of equivalent aerofoil. Keeping the same pitching amplitude, we can obtain the aerodynamic curve of the equivalent aerofoil in every phase by calculating the pitch motion under different mean angles of attack. Then, these aerodynamic data can be used as the input for the augmented lifting-line method in every time step.

Figure 2. Schematic of the augmented lifting-line method.

Figure 3. Illustration of equivalent aerofoil.

For example, $\alpha = {\alpha _0} + 4.0 \cdot \sin (\omega t)$ , the mean angle-of-attack ${\alpha _0}$ varies from - 8° to 8°. If the angle interval is 1º, there are 17 calculation cases for 2D CFD simulation, and we can obtain 17 aerodynamic curves. If the phase angle is 90° (t/ T = 0.25), by taking the aerodynamic data at 90° phase angle in these 17 aerodynamic curves, the aerodynamic curve of the equivalent aerofoil at 90° phase angle can be obtained. The other phases are similar, so we can obtain the equivalent aerofoil of all phases in one cycle.

2.2 Unsteady Kutta-Jouowski(K-J) Theorem

The steady lifting- line method uses K-J theorem for the circulation iteration, and finally realises the convergence between the lift of angle-of-attack and the lift of the circulation. For unsteady conditions, the unsteady K-J theorem must be considered. Reference [Reference Katz and Plotkin1] presents the lift formula of 2D aerofoil under unsteady flows as follows:

(7) \begin{align}L'(t) = \rho {U_\infty }\Gamma (t) + \rho \frac{{\partial \Gamma (t)}}{{\partial t}}c\end{align}

In Equation (7), the lift consists of two parts in unsteady condition. The first term is the proportional term, which represents the circulation lift, and the second term is the time derivative term, which represents the influence of the additional mass.

For the jth section on the wing, the time derivative term of the unsteady K-J theorem is discretised by first-order explicit scheme:

(8) \begin{align}{L'_{\!\!j}}(t) = \rho {U_\infty }{\Gamma _j}(t) + \rho \frac{{({\Gamma _j}(t) - {\Gamma _j}(t - \Delta t))}}{{\Delta t}}{c_j}\end{align}

After the derivation, the expression of the circulation in unsteady flows can be obtained as follows:

(9) \begin{align}{\Gamma _j}(t){\rm{ = }}\frac{{{{L'}_{\!\!j}}(t) \cdot \Delta t + \rho {\Gamma _j}(t - \Delta t){c_j}}}{{\rho \left[ {{U_\infty }\Delta t + {c_j}} \right]}}\end{align}

In each time step, once the lift changes, the circulation will be updated.

2.3 Circulation correction factor

For a 2D moving aerofoil, its wake is constructed by the infinite span spanwise vorticity. However, for a moving wing with finite span, its wake contains both spanwise vorticity and streamwise vorticity, and the length of the spanwise vorticity is finite, which leads to different influence for the bound vorticity compared to the 2D one. Thus, in the augmented lifting-line method, we add the circulation correction to consider this difference.

Studies in references [Reference Smyth, Young and Mare16, Reference Ahmadi and Widnall29] have shown that when aspect ratio or reduction frequency is large enough, 2D and 3D results will tend to be consistent. Therefore, the expression of the circulation correction factor in this paper is given as follows:

(10) \begin{align}F = {f_1} \cdot {f_2}\end{align}

The circulation correction factor consists of the product of two factors.

(11) \begin{align}{f_1} = 1.0 - \frac{{2.0}}{\pi }arccos\left( {{e^{\frac{{ - 4\left(\frac{b}{2} - \left| z \right|\right)}}{{A{R^2}\left| z \right|}}}}} \right)\end{align}

$\left| z \right|$ indicates the span position of the section. ${f_2}$ is related to the sign of the circulation difference between adjacent time step.

(12) \begin{align}{f_2}{\rm{ = }}\left\{ {\begin{array}{*{20}{l@{\quad}l}}{\dfrac{2}{{AR}}\dfrac{k}{{{k^2} + 1}}} & {\Gamma (t) - \Gamma (t - \Delta t) \gt 0}\\[12pt] 0 & {\Gamma (t) - \Gamma (t - \Delta t) = 0}\\[9pt] {{e^{ - \left(\frac{{A{R^2}}}{{1000}}\right)}}{e^{ - \left(\frac{{{k^2}}}{3}\right)}}}& {\Gamma (t) - \Gamma (t - \Delta t) \lt 0}\end{array}} \right.\end{align}

The circulation correction factor $F$ is related to the aspect ratio, the reduction frequency and the span position of the section. ${f_1}$ reflects the circulation correction along the spanwise, which is similar to the Prandlt’s tip loss factor [Reference Branlard30] and closer to 1 at wingtip; ${f_2}$ represents the correction of the 2D spanwise vorticity. Through the circulation correction factor, the improved method can consider the difference between 2D and 3D spanwise vorticity to a certain extent.

2.4 Aerodynamic force calculation of 3D wing

The improved method can obtain the aerodynamic force coefficient and induced angle-of-attack of each section on the wing during the calculation process, and the aerodynamic force of 3D wing at each time step can be calculated based on these quantities.

Figure 4 illustrates the downwash of a finite wing. For the $j{\rm{th}}$ section, the span is $\Delta {b_j}$ , and the chord is ${c_j}$ . ${C_l}$ and ${C_d}$ are aerodynamic coefficient correspond to effective angle-of-attack ${\alpha _e}$ , we can obtain aerodynamic forces per unit span,

(13) \begin{align}\begin{array}{*{20}{c}}{{{L'}_{\!\!j}} = \dfrac{1}{2}\rho U_\infty ^2{c_j}{C_l}}\\[14pt]{{{D'}_{\!\!j}} = \dfrac{1}{2}\rho U_\infty ^2{c_j}{C_d}}\end{array}\end{align}

Figure 4. Effect of downwash over a local aerofoil section of a finite wing.

If the induced angle-of-attack is ${\alpha _i}$ , the lift ${L_j}$ and drag ${D_j}$ relative to freestream can be expressed as follows:

(14) \begin{align}\begin{array}{*{20}{c}}{{L_j} = \dfrac{1}{2}\rho U_\infty ^2\Delta {b_j}{c_j}({C_l}\cos {\alpha _i} - {C_d}\sin {\alpha _i})}\\[14pt] {{D_j} = \dfrac{1}{2}\rho U_\infty ^2\Delta {b_j}{c_j}({C_l}\sin {\alpha _i} + {C_d}\cos {\alpha _i})}\end{array}\end{align}

Figure 4 shows the effect of downwash. The total force of the 3D wing can be obtained by summing the aerodynamic force of Equation (14) along the span.

The following will use the improved method to perform a series of calculations, and compare results with the 3D CFD calculations to verify the effectiveness of the improved method.

2.5 CFD method verification

In order to compare with CFD method, this paper first tests the effectiveness of CFD method. The NACA 0012 aerofoil for pitching motion is selected for the verification.

The sliding mesh method in the commercial software Fluent is used to simulate the pitch motion, and the calculation is based on the density-based solver. The velocity inlet condition is adopted for the left, upper and lower boundaries of the computation domain, while the pressure outlet boundary is adopted for the right. The non-slip condition is used on the aerofoil surface.

The computation domain shows in Fig. 5, which is divided into two zones. The outer zone is rectangular with a length of $160c$ and a width of $120c$ . The inner zone is a circular zone with a radius of $2.0c$ , and the centre of the zone is rotation position of the aerofoil, which is $0.37c$ from the leading edge of the aerofoil. During the process of the computation, the inner zone moves with the aerofoil. The inner and outer zones interact with each other by the interpolation on the interface surface.

Figure 5. CFD mesh of NACA0012 aerofoil.

The aerofoil surface is divided into 268 mesh elements, and a value of ${y^ + }$ below 1 is realised to ensure the accuracy of wall flow. The mesh element in the inner zone is 25,000, while the mesh element in the outer zone is 21,000, and the total mesh element is 46,000.

The general form of pitching motion is

(15) \begin{align}\alpha = {\alpha _0} + {\alpha _m}\sin (\omega t)\end{align}

In Equation (15), ${\alpha _0}$ is the mean angle-of-attack, and ${\alpha _m}$ is the amplitude of the oscillation. In this case, ${\alpha _0}{\rm{ = }}{0^{\rm{o}}}$ , ${\alpha _m}{\rm{ = }}{6.74^{\rm{o}}}$ , $k = \dfrac{{\omega c}}{{2{U_\infty }}} = 0.4$ .

The free stream velocity is ${U_\infty } = 34.0294m/s$ , and the air density is $\rho {\rm{ = }}1.225kg/{m^3}$ . The chord length is $c{\rm{ = }}0.43m$ , and dynamic viscosity is $\upsilon {\rm{ = }}1.7894 \times {10^{{\rm{ - }}5}}kg/(m \cdot s)$ , corresponding to Reynolds number $Re = 1.0 \times {10^6}$ .

The incompressible Unsteady Reynolds-averaged Navier-Stokes equations are solved in conjunction with the S-A turbulence model. Every cycle of pitching motion takes 250 time steps, and the number of iterations in each time step is 50 to ensure the convergence. After four cycles, the motion reaches a stable period. Take the lift coefficient in the fifth cycle and compare it with the experiment [Reference Halfman31] as follows:

Figure 6. Comparisons of lift coefficient.

It can be seen in Fig. 6 that the lift coefficient calculated by CFD method is in good agreement with the experiment, indicating that the CFD method used in this paper is reliable.

The second verification example is NACA0015 aerofoil for pitching motion. In this case, ${\alpha _0}{\rm{ = 4}}{\rm{.}}{{\rm{0}}^{\rm{o}}}$ , ${\alpha _m}{\rm{ = 4}}{\rm{.}}{{\rm{3}}^{\rm{o}}}$ , $k = 0.133$ . More flow conditions can be found in reference [Reference Piziali32]. Take the lift and drag coefficients and compare them with the experiment [Reference Piziali32] as follows:

Figure 7. Comparisons of lift and drag coefficient.

It can be seen in Fig. 7 that CFD results is in good agreement with the experiment data, which further validates the numerical method.

In the following, the pitch motion of different wings will be calculated by the CFD method, and the results will be compared with the improved method at the same time.

3.1 Pitching motion for large aspect ratio wing

Firstly, the NACA 0012 wing with AR = 20 is selected to verify the effectiveness of the improved method for the calculation of unsteady aerodynamic force of wing with large aspect ratio.

3.1.1 CFD calculations

In order to perform the improved method, the aerodynamic data of 2D equivalent aerofoil must be obtained. Thus, NACA 0012 aerofoil undergoing pitching motion is calculated by CFD method firstly.

The free stream velocity is ${U_\infty } = 30m/s$ , and the air density is $\rho {\rm{ = }}1.225kg/{m^3}$ . The chord length is $c{\rm{ = 1}}{\rm{.0}}m$ , and dynamic viscosity is $\upsilon {\rm{ = }}1.7894 \times {10^{{\rm{ - }}5}}kg/(m \cdot s)$ , corresponding to Reynolds number $Re = 2.0 \times {10^6}$ . The CFD calculation method remains the same as the 2.6 section; 250 time steps are taken in one cycle, corresponding to 250 equivalent aerofoils.

The amplitude of the oscillation remains ${\alpha _m}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ . Given the reduced frequency $k = 0.3$ , the mean angle-of-attack ${\alpha _0}$ is −8° to 8°, the lift curves corresponding to the equivalent aerofoil under five typical phases in a cycle are presented as follows:

Figure 8. Lift curves of equivalent aerofoil at different phases.

In Fig. 8, one lift curve reflects the lift characteristics of the equivalent aerofoil at a phase angle. For the same angle-of-attack, taking the force coefficients of different lift curves, the lift characteristics of the aerofoil in a motion cycle can be obtained. It can be seen that the lift curve of the equivalent aerofoil is similar to the ordinary aerofoil, which keeps linear at small angles of attack.

For 3D CFD calculation of the wing, the calculation adopts a half-model considering the symmetry and the mesh is showed in Fig. 9. The length of the X and Y directions of the computational domain is consistent with 2D calculation, and the length of the Z direction is $160c$ . The computational domain is also divided into two zones. The inner zone is a cylindrical zone with a mesh of 1.9 million, which can rotate with the wing; the outer zone is a stationary zone with a mesh of 2.6 million and a total mesh for 3D calculation is 4.5 million.

Figure 9. 3D mesh of NACA 0012 straight wing.

Symmetry boundary condition is used for the planes of symmetry in the inner and outer zones, and the other settings are consistent with the 2D CFD calculation. The sliding mesh method is used to calculate the pitching motion again.

3.1.2 Comparisons of different methods

The reduction frequencies are $k = 0.1$ , 0.3 and 1.0, respectively. For pitching motion of ${\alpha _0}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ and ${\alpha _m}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ , the lift and drag coefficients obtained by 2D and 3D CFD methods are compared as follows:

Figure 10. Comparisons of 2D and 3D CFD results for NACA 0012 wing with AR = 20.

In Fig. 10, it can be seen that there exists difference between the aerodynamic curves of the aerofoil and the wing, which comes from the three-dimensional effect. With the increase of reduction frequency, the effect of additional mass is enhanced, and the phase of aerodynamic curve is advanced [Reference Yu and Bernal33] for both 2D and 3D.

The improved method can be called UNLL (unsteady lifting-line). Firstly, the influence of the number of horseshoe vortices on the solution of UNLL is studied. The wing is divided into 10, 30, 50 and 80 elements, respectively. For $k = 0.3$ , the lift and drag coefficients calculated by UNLL with different number of horseshoe vortices are compared as follows:

Figure 11. Comparisons of UNLL results with different number of horseshoe vortices.

It can be seen in Fig. 11 that for N = 10, the lift coefficient calculated by UNLL has some differences with the other three curves. However, for N > 30, UNLL results are basically not affected by the number of horseshoe vortices. Thus, 50 horseshoe vortices along the spanwise is sufficient for UNLL in this paper.

The comparisons between UNLL and 3D CFD results are as follows:

Figure 12. Comparisons of UNLL and 3D CFD results for NACA 0012 wing with AR =  20.

In Fig. 12, the aerodynamic curves calculated by UNLL are in better agreement with 3D CFD results compared with 2D CFD results at every reduced frequency. Thus, the improved method can capture the three-dimensional effect and can be used for the calculation of unsteady aerodynamic force of wings with large aspect ratio.

3.1.3 Comparisons of time cost

For the pitching motion of ${\alpha _0}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ and ${\alpha _m}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ , 2D CFD results with the mean angle-of-attack −2° to 5° are needed in the lifting-line method. The angle interval is taken as ${1^{\rm{o}}}$ , and a total of eight 2D cases should be calculated.

Table 1. The time cost of UNLL and 3D CFD

The comparisons of time cost between UNLL and 3D unsteady CFD calculation are reported in Table 1.

In Table 1, we can be seen that the time cost on the improved method is mainly in 2D URANS simulations, and the cost of the lifting-line is small. Compared with 3D unsteady CFD simulation, the computational efficiency for unsteady aerodynamic simulation is significantly improved at least 12 times by the improved method. Thus, the improved method is suitable for the preliminary design in engineering.

3.2 Pitching motion for wing with medium aspect ratio

Because the improved method is based on the lifting-line method, it’s convenient to predict the aerodynamic force of 3D wings under different aspect ratio. Based on the same 2D URANS data, NACA 0012 wings with aspect ratio AR = 10 and AR = 5 are selected for pitching calculations again.

Figure 13. Comparisons of 2D and 3D CFD results for NACA 0012 wing with AR = 10.

Figure 14. Comparisons of 2D and 3D CFD results for NACA 0012 wing with AR = 5.

Figure 15. Comparisons of UNLL and 3D CFD results for NACA 0012 wing with AR = 10.

Figure 16. Comparisons of UNLL and 3D CFD results for NACA 0012 wing with AR = 5.

The reduction frequencies are $k = 0.1$ , 0.3 and 1.0, respectively. For ${\alpha _0}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ and ${\alpha _m}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ , the calculation results of different methods are compared as follows:

From Figs 13 and 14, it can be seen that the lift curve of 2D aerorfoil is always above the lift curve of 3D wing, but the drag curve is just the opposite. As the aspect ratio decreases, the difference between 2D and 3D becomes larger, indicating that the downwash of the wake is enhanced, and the 3-D effect becomes more and more obvious.

From Figs 15 and 16, the lift curve obtained by the improved method is in good agreement with 3D CFD under different aspect ratio. For the drag curve, although the error at the peak has increased, the drag obtained by the improved method is more consistent with 3D results throughout the cycle compared with 2D CFD. Therefore, the improved lifting-line method in this paper is also suitable for the estimation of unsteady aerodynamic forces of wings with medium aspect ratio.

3.3 Pitching motion for NACA 4412 wing

In order to verify the applicability of the improved method to different type of wings, the NACA4412 wing was selected to perform unsteady CFD calculations, and then the results were compared with the improved method as well.

The CFD calculation method remained the same as the previous section. Firstly, 2D CFD pitching calculations were carried out to obtain aerodynamic data of equivalent aerofoils. Then, NACA4412 wing with AR = 20 was selected for the verification of the improved method.

The reduction frequencies are $k = 0.1$ , 0.3 and 1.0, respectively. For ${\alpha _0}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ and ${\alpha _m}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ , the calculation results of different methods are compared as follows:

Figure 17. Comparisons of 2D and 3D CFD results for NACA 4412 wing with AR = 20.

In Fig. 17, there is still an obvious difference between 2D and 3D results. However, the aerodynamic curves obtained by the improved method in Fig. 18 are in good agreement with the 3D CFD results, indicating that the proposed method is also effective for the wing with camber.

Figure 18. Comparisons of UNLL and 3D CFD results for NACA 4412 wing with AR = 20.

From another point of view, the improved method can also be regarded as a reduced order model. Based on the idea of lifting-line, the mapping relationship between 2D and 3D is constructed, and low-dimensional aerodynamic data is used to predict high-dimensional aerodynamic forces, achieving the purpose of time-saving for unsteady calculations.

3.4 Pitching motion for tandem configuration

The above calculations are all about the single wing, but in engineering, there are often multiple wing aerodynamic problems, such as biplane wing and tandem wing [Reference Hosseini, Tadjfar and Abba34]. For 3D multiple wings CFD calculations, the time cost on mesh generation and computation will also increase, so it is necessary to develop a fast estimation method for these problems. Thus, the tandem configuration is selected in order to verify the effectiveness of the improved method in calculating unsteady aerodynamic forces of multiple wings.

The tandem shows in Fig. 19. The front wing adopts NACA 0012 aerofoil with a geometric installation angle of 4°, and the rear wing adopts NACA 4412 aerofoil with a geometric installation angle of 0°, and the front and rear wing chords are both 1.0 m. In this case, the front wing will perform the pitching motion while the rear wing keep stationary. The pitching centre of the front wing is around its own 0.25c position and the distance between the front and rear wings is D = 4 m. The span of front and rear wing is 10 m, corresponding to AR = 10. 3D mesh is showed in Fig. 20.

Figure 19. Schematic of tandem wing.

Figure 20. 3D mesh of tandem wing.

In the previous sections, we expands the concept of aerofoil and considers that every phase of the pitching motion corresponds to an aerofoil. Furthermore, multiple aerofoil can also be considered as a single equivalent aerorfoil, and we can use 2D CFD results of multiple aerofoils to calculate the unsteady aerodynamic force of 3D multiple wings. In other words, the multiple wings can be considered as a single equivalent wing, which avoids modeling each wing separately and can further improve the calculation efficiency for complex configuration.

Figure 21 shows the equivalent aerofoil for the tandem. For each phase, the aerodynamic force of the equivalent aerofoil is defined as the sum of the front and rear aerofoils, and the aerodynamic force of the equivalent wing is the sum of the front and rear wings. By changing the direction of the free stream flow, the 2D aerodynamic force of the equivalent aerofoil at different angles of attack can be obtained. Then, the aerodynamic force of 3D equivalent wing is calculated based on the improved method. The reference chord length and reference area of the equivalent wing are defined as the chord length and area of a single wing.

Figure 21. Equivalent aerofoil for tandem configuration.

The CFD calculation method remains the same as the previous section. The front wing performs periodic pitching movement with ${\alpha _m}{\rm{ = }}{{\rm{4}}^{\rm{o}}}$ , while the rear wing remains stationary. The reduction frequencies are $k = 0.1$ , 0.3 and 1.0. For 0º angle-of-attack, the total aerodynamic force of the tandem configuration calculated by different methods is compared as follows:

In Fig. 22, it still can be seen obvious difference between 2D and 3D aerodynamic curves of the tandem configuration. However, the aerodynamic curves obtained by the improved method in Fig. 23 are in good agreement with 3D CFD results, indicating that the improved method can also be used to estimate the overall unsteady aerodynamic force of 3D multiple wing configuration.

Figure 22. Comparisons of 2D and 3D CFD results for tandem wing.

Figure 23. Comparisons of UNLL and 3D CFD results for tandem wing.