1. Introduction

1.1. Background

The NACA airfoils are a series of airfoils proposed by the National Advisory Committee for Aeronautics (NACA) in the U.S., where each airfoil is assigned a four-digit number that describes its geometric properties. The first digit represents the camber in percent of the chord, the second digit describes the distance between the maximum camber and the leading edge of the airfoil, in tenths of the chord, and the last two digits indicate the ratio between the maximum thickness and the chord length, in percentage [1]. The NACA 0012 airfoil is a symmetric airfoil, which has no camber, with a maximum thickness of 12% of the chord. It is a standard airfoil that is significant in aerodynamic research, and is commonly used as a CFD benchmark. Its section is also commonly used on helicopter blades and airplane tails [2].

CFD, or Computational Fluid Dynamics, plays an important role in aerodynamic design, since it can reduce the experimental wind-tunnel tests and hence reduce the research cost [3]. The modern turbulence models, such as k-ω SST, k-ε, LES, LBM, etc., has already been developed to accurately predict aerodynamic coefficients [4].

Beyond research on single airfoil, scientists are now taking accounts of studying tandem-wing configurations. The tandem-wing configurations can enhance the airfoils’ aerodynamic performance over single airfoil configuration [5]. It also offers potential benefits, such as increasing the lift to drag ratio of the downstream wing due to wake interactions. When a higher lift to drag ratio is attained, the aircraft will gain more lift and less drag, and thus reduce the amount of energy required for flights [6, 7].

1.2. Previous studies

Researchers have put huge effort into studying the behavior of NACA 0012 airfoil, using various methods. Khalid (2022) applied CFD on analyzing the relationship between the angle of attack (AoA) and the lift coefficient of one airfoil, and reached a result that the lift coefficient reaches its maximum value at AoA = 14° [6]. In 2018, Li compared the performance of smooth and rough airfoils [8]. Kurtulus in 2019 studied the lift to drag ratio of a pitching airfoil at low Reynolds number, which helped understanding the flapping airfoil vortex formation [9].

Previous studies on tandem airfoil configurations show that axial spacing and phase difference between oscillating airfoils critically influence downstream lift and drag cycles, since upstream wake disrupts or energizes the downstream boundary layer. In 2011, Lee discovered that two in-tandem NACA 0012 airfoils oscillating sinusoidally with a phase difference of 180° will improve the downstream performance significantly, compared to the in-phase motion [7].

However, most studies on one airfoil only consider the flapping motion or the pitching motion individually, instead of a multi-step analysis. Also, the studies on tandem airfoils only consider limited phase differences between the two motions, where the phase differences except 0 or π may have better performance. Moreover, most studies fail to build the connection between the distance between the tandem airfoils, the phase difference between the pitching motion of two airfoils, and the lift to drag ratio, which is important on understanding how the vortices produced by the front airfoil affect the rear airfoil. Finally, a lot of previous research considers the lift force and drag force individually, instead of analyzing the lift to drag ratio. However, simply maximizing the lift force may not reach an optimal configuration that minimizes the energy, since the drag force could still be high.

1.3. Contribution

This research uses COMSOL Multiphysics to conduct 2D CFD simulations on the NACA 0012 airfoils, and aimed to reach an optimal configuration that maximizes the lift to drag ratio of the airfoils. This research first simulates one airfoil that has a linear motion and a rotational motion, to get the optimal amplitude of rotational motion and the optimal phase difference between the two motions. It then simulates two airfoils, placed one behind another, with the motions obtained in the simulation of one airfoil. The study further quantifies how the phase difference in the linear motions of the two airfoils under different distances between the two airfoils affects the maximum lift to drag ratio of two airfoils, and the relative L/D enhancement of the downstream wing (rear airfoil) compared to the upstream wing (front airfoil).

This research contributes to deepening the understanding of wake‐induced aerodynamic interactions and to the study of energy harvesting devices for airfoils. It is also useful in guiding the design new aircrafts, such as the aircrafts with multi‐wing systems.

2. Methodology and results

2.1. Setup

The computational domain of the experiment is a rectangular region with a size of 20 m × 8 m, centered at the origin. The chord length of a standard NACA 0012 airfoil is c = 1.8 m. Laminar flow with a free-stream velocity of U_∞ = 0.01 m/s is used as the flow field, with a free-stream density of ρ_∞ = 1.225 kg/m³, and a dynamic viscosity of μ_∞ = 1.83 × 10^-5 kg/(m·s). The Reynolds number is calculated by

$$$ \begin{array}{c}\text{R}\text{e}\text{=}\frac{{\text{U}}_{\text{∞}} {\text{ρ}}_{\text{∞}}\text{ c}}{{\text{μ}}_{\text{∞}}}\text{=}\text{1204.9}\text{.}\end{array}$$$(1)

The inlet velocity of the rectangular region is set to be U_∞, entering from the left boundary (Figure 1), while the outlet pressure is p₀ = 0 Pa, which is the boundary condition of the right side (Figure 1). The AoA of the airfoil is denoted by α. Moreover, the lift force, L, is calculated by the line integration over the airfoil:

$$$ \begin{array}{c}\text{L}\text{=}\int \left({\text{T}}_{\text{x}}\text{sin}\text{α}\text{-}{\text{T}}_{\text{y}}\text{cos}\text{α}\right)\text{ d}\text{l}\text{.}\end{array}$$$(2)

Similarly, the drag force, D, is calculated as:

$$$ \begin{array}{c}\text{D}\text{=}\int \left(\text{-}{\text{T}}_{\text{x}}\text{cos}\text{α}\text{-}{\text{T}}_{\text{y}}\text{sin}\text{α}\right)\text{d}\text{l}\text{.}\end{array}$$$(3)

Where T_x and T_y are the total stress on the x and y direction, respectively [10].

2.2. Single airfoil

The airfoil is centered at the origin. The mesh grid for computation is chosen to be “Extra Fine” airfoil in the COMSOL. The geometry of the computation domain and the mesh grid are shown in Figure 1.

Figure 1. The geometry (a) and the mesh grid (b)

The airfoil has a linear (flapping) motion and a rotational motion, where linear motion means that the airfoil moves periodically along the y direction, and the rotational motion means that the airfoil rotates periodically along an axis perpendicular to the plane for simulation and locate at the leading edge of the airfoil. The linear motion is described by

$$$ \begin{array}{c}\text{Δ}\text{h}\text{=}\text{0.01}\text{sin}\text{(}\text{ωt}\text{)}\text{,}\end{array}$$$(4)

and the rotational motion can be achieved by changing the angle of attack of the airfoil which is described by

$$$ \begin{array}{c}\text{α}\text{=}\text{15}\text{+}\text{A}\text{s}\text{i}\text{n}\text{(}\text{ω}\text{t}\text{+}\text{φ}\text{)}\text{,}\end{array}$$$(5)

where t represents the time, A is the amplitude of the rotational motion, φ is the phase difference between two motions, and ω in both formulas represent the angular frequency of the periodic motions. This experiment aims to maximize the average lift to drag ratio, L/D, in one period, by varying A and φ, with three angular frequencies, 1, 10, and 100.

The simulations are conducted using COMSOL, with a time-dependent solver. The experiments are conducted within one period of the motion of the airfoil, and the time steps are set to be 0.05/ω seconds. Each pair of parameters will have 125 data points, and their average lift to drag ratio is recorded.

Firstly, fix φ = 0 and change the amplitude A. Choosing A to be the integers between and including 0 and 20, the average lift to drag ratios are reached as results, and are recorded in Table 1 below.

By Table 1, whatever ω is, the maximum lift to drag ratio exists at A = 12.

Now fix A = 12 and set  $$ \text{φ }\text{= }\text{kπ}\text{/4}$$, where k are the integers between (and including) 0 and 8. The experiment results are listed in Table 2, and shown in Figure 3.

By Table 2, the lift to drag ratio is maximized at φ = 3π/4. Hence, the optimal choices of the amplitude and phase difference are A = 12, φ = 3π/4. The flow field of the optimal choices are shown in Figure 4.

Figure 4. The flow field for single airfoil with ω = 1 (a), 10 (b), 100 (c), with A = 12 and φ = 3π/4, at time t = 0s

2.3. Dual airfoil system

This experiment explores two identical airfoils, placed one behind another, with a distance ℓ measured from the front airfoil’s tail to the back airfoil's tail. The distance ℓ was chosen to be 2.8, 6.8, and 11.8 meters. The configurations are shown in Figure 5.

Each of the airfoils has a periodical linear motion and rotational motion. In this case, there is a phase difference ψ between two airfoils’ linear motion. That is, the airfoil in the front has a linear motion described by

$$$ \begin{array}{c}{\text{Δ}\text{h}}_{\text{1}}\text{=}\text{0.01}\text{sin}\text{(}\text{ωt}\text{)}\text{,}\end{array}$$$(6)

while the airfoil at the back has a linear motion described by

$$$ \begin{array}{c}{\text{Δ}\text{h}}_{\text{2}}\text{=}\text{0.01}\text{sin}\text{(}\text{ωt }\text{+ }\text{ψ}\text{)}\text{.}\end{array}$$$(7)

Furthermore, the rotational motion is

$$$ \begin{array}{c}\text{α}\text{=}\text{15}\text{+}\text{12}\text{s}\text{i}\text{n}\text{(}\text{ω}\text{t}\text{+}\text{7}\text{π}\text{/}\text{4}\text{)}\text{.}\end{array}$$$(8)

Here the phase difference between the rotational motion and the linear motion is chosen to be 7π/4, since when φ = 3π/4, the COMSOL cannot generate convergent solutions.

In this experiment, the lift to drag ratio is maximized by varying the phase difference between two airfoils, ψ. The experiment was conducted under different distance between two airfoils, ℓ, and the angular frequency, ω. As a result, a dimensionless parameter is introduced here:

$$$ \begin{array}{c}{\text{t}}^{\text{*}}\text{=}\frac{\mathcal{l}\text{ω}}{{\text{U}}_{\text{∞}}}\text{,}\end{array}$$$(9)

so that the experiment was conducted under different t^*. The values of ω (in s^-1) are chosen to be 10, 50, and 100. The phase different, ψ, was chosen to be 0, π/4, π/2, π, and 3π/2.

All values of t^* are shown in Table 3 below.

In this case the lift to drag ratios of the two airfoils are calculated separately, by linear integration on each airfoil. The time intervals for simulation are the same as in the single airfoil case.

The lift to drag ratio of the front airfoil is shown in Table 4 and Figure 6.

The lift to drag ratio of the back airfoil is shown in Table 5 and Figure 6.

Table 4 and Table 5 show that, in most cases, the lift to drag ratios of the rear airfoil are larger than these of the front airfoil; and that the lift to drag ratios are maximized when ψ = 3π/2. The flow field of several situations are shown in Figure 8 below.

Figure 8. Flow field around the airfoils, at time t = 0

3. Discussion

3.1. Single airfoil result

The simulation of single airfoil shows that under low Reynolds number, a NACA 0012 airfoil with a periodical linear motion and a periodical rotational motion that starts at the angle of attack is 15° has a maximum lift to drag ratio when the amplitude of the rotational motion is A = 12, and the phase difference between the two motions is φ = 7π/4. This result is consistent for all three frequencies (ω = 1, 10, 100) of the linear motion.

Considering the ω = 1 case: by Table 1 shown, when there is no phase difference between the motions, at A = 0, which means only linear motion exists, the lift to drag ratio of the airfoil is L/D = 3.72. While at A = 12, L/D = 5.57. The ratio has increased by 49.62%. For the other two cases, ω = 10 and ω = 100, L/D are increased by 49.71% and 50.31%, respectively. Among the three frequencies, when ω = 100, L/D = 5.61 is the largest value, and it is also increased the most when ω = 100.

Furthermore, Table 2 shows that when ω = 1, the lift to drag ratio is L/D = 6.32 when there is a phase difference of φ = 3π/4 between the two motions. The lift to drag ratio is further increased by 13.46%. At ω = 10 and ω = 100, L/D is further increased by 13.60% and 13.93%, respectively. Again, ω = 100 results in the largest value of L/D among the three frequencies, which is 6.38. The lift to drag ratio is increased the most when ω = 100.

3.2. Dual airfoil result

The dual-airfoil experiment shows that the lift to drag ratio is generally larger in the rear airfoil than in the front airfoil. This is caused by the vortices shed by the front airfoil, which convects downstream with the wake and lift the rear airfoil [11]. By Table 4, when t^* ≤ 11800, the L/D of the front airfoil is maximized when ψ = 0; when t^*=34000, it is maximized when ψ = π/2; for other t^* values, it is maximized when ψ = 3π/2. By Table 5, the rear airfoil has L/D maximized at ψ = 0 when t^* ≤ 6800, at ψ = π when t^* = 68000, and at ψ = 3π/2 for other t^* values. The maximum L/D exists at t^* = 14000, where L/D_front = 6.143 and L/D_rear = 6.213.

Comparing the data recorded in Table 4 and Table 5, most of the L/D values are larger in the rear airfoil than the front airfoil. An unusual observation is that at t^*=34000, the L/D of the front airfoil is maximized at ψ = π/2, which is 6.173, while the maximum L/D of the back airfoil is 6.171. However, Figure 6 shows that all t^* values other than 34000 will result in the L/D to be minimized near ψ = π/2. Hence, this deviation is most probably caused by anomalous spikes in the simulation results, due to the time steps or mesh grid taken.

Moreover, by Table 4 and Table 5, the rear airfoil outperformed the front by 1.40% L/D at t^* = 14000, which is the maximum percentage. At t^* = 6800, the L/D of the rear airfoil is 0.08% less than that of the front airfoil, and at t^* = 68000 and ψ = 3π/2, the L/D of the rear airfoil is only 0.03% more than that of the front. This can be explained by that when t^* ≤ 6800, the vortices are insufficiently developed, and when t^* ≥ 68000, the vortex decays before reaching the rear airfoil.

3.3. Limitations

The experiments use COMSOL as the tool of simulation. “Extra fine” meshes are used in both simulations of single and double airfoils, since this is the only mesh that can make most solutions converge.

Additionally, in the simulation of double airfoils, instead of α = 15 + 12 sin (ωt + 3π/4), the rotational motion is set to be α = 15 + 12 sin (ωt + 7π/4). The reason is that when φ = 3π/4, the COMSOL solver encounters repeated error test failures, which may be because it has reached a singularity. This problem cannot be resolved by changing the mesh of the simulation. However, situations could be different from φ = 3π/4. Thus, further experiments are needed. One may try to refine the mesh grid, or use a different simulation tool to solve the problem.

Moreover, in these simulations, the lift to drag ratios are only calculated within one period. To make the results more reliable, one can repeat the simulation with more periods. Also, the time intervals used in the simulation is 0.05/ω, which means that each simulation will generate only 125 data points. To get a more precise result, the time intervals can be further decreased to 0.01/ω, or less.

Finally, the flow in the rectangular region is the laminar flow, and the simulations are conducted under a low Reynolds number. For more sophisticated situations, turbulent flows and a larger Reynolds number can be used.

4. Conclusion

This study systematically investigates the aerodynamics of flapping NACA 0012 airfoils through integrated CFD simulations at low Reynolds numbers. Key findings demonstrate that, for one airfoil with a flapping motion and a pitching motion around AoA = 15°, the amplitude A = 12 of pitching motion and the phase difference φ = 3π/4 will amplify the lift to drag ratios by up to 71.11%, compared to A = 0 and φ = 0 situations. Tandem airfoil systems achieve optimal configurations when two airfoils pitching with a phase difference of ψ = 3π/2, which enhances the L/D of the rear airfoil by at most 1.40% through leading edge vortex impingement. The dimensionless parameter, t^*, governs efficiency of the vortex, peaking at t^* = 14000.

These results are useful in designing new aircraft, and new energy recovering systems. However, limitations exist in this research, including the insufficient density of the mesh grid chosen, the time interval and number of periods chosen, laminar flow assumptions, etc. Future work may build denser mesh grids, conduct simulation in more time periods, and implement turbulence models for high Reynolds numbers validation.