Aerodynamic response of micro rotary wings in wind perturbation

3.1. Governing equations

Under the assumption of continuous medium, three dimensional unsteady compressible Reynolds averaged Navier-Stokes (RANS) equations are given by the integral conservation form:

where:

with:

and $\rho$, $p$, $T$, $\tau$ and $k$ are density, pressure, temperature, shear stress and heat conductivity coefficient. $u$, $v$, $w$ denote three Cartesian velocity components of fluid velocity $\mathbf{q}$. $E$ and $H$ express total energy and total enthalpy, and their relation is given by the formula $H=E+p/\rho$. ${\mathbf{I}}_x$, ${\mathbf{I}}_y$, ${\mathbf{I}}_z$ are three unit vectors along three coordinate axes of rectangular coordinate system, and $\partial V$ with unit normal vector $\mathbf{n}=(n_x,n_y,n_z)$ is the boundary of the control volume $V$.

Eqs. (1) are not only suitable for static control volume, but also for arbitrary motion state. However, it is expected that the time derivative can be gotten out of the integral formula so as to conveniently solve equations with the finite volume method. According to the Leibniz theorem [24], there is the formula:

where ${\mathbf{q}}_b$ is the boundary velocity vector of control volume. In other words, ${\mathbf{q}}_b$ is mesh velocity vector. If Eq. (2) are plugged into Eq. (1), we can get the RANS equations with the integral form:

with:

To simulate the flow field of aircraft flight in the wind disturbance environment, it is supposed that perturbation of wind has the temporal synchronization in the computational domain considering the small size of the aircraft and the small wind perturbation, so the field velocity method is adopted. The flow fields can be construed as obtaining an equal and opposite velocity according to the relative motion. Based on the physical principle of the above building mathematical model, the convective terms ${\mathbf{H}}_1$ of the RANS equations should be changed to:

Correspondingly, the governing equations are converted to:

and ${\mathbf{q}}_g$ is the speed vector relative to the atmospheric disturbance. In the paper, the grid velocity for the wind disturbance is only carried out in the two regions of rotor grids. The information is transferred between the rotor grids and intermedial grid, as well as the intermedial grid and background grid via artificial external boundaries and hole boundaries, which embodies the advantage of moving overset girds.

The state equations are needed to guarantee the closure of control equations, which are written by:

where $\gamma$ is the ratio of specific heats. The eddy viscosity is modeled by the zeroth-order algebraic turbulence model provided by Baldwin Lomax [25].

3.2. Conversion of coordinate system

Since the rotary wing motion is broken down by using three levels of overset grids, coordinate transformation is required to carry out among sub domain grids including intermedial grid, background grid and rotor grids during designing program of the unsteady flow field. Three kinds of coordinate systems are involved. One of the related coordinate systems is inertial coordinate system also called absolute coordinate system written as $R_A(O,x,y,z)$. Rotational coordinate system $R_R(O,e_1,y,e_3)$ is obtained after the inertial coordinate rotates ($90°-\psi$) around the $y$ axis, and rotor coordinate system $R_P(O,b_1,b_2,a_3)$ is obtained after the rotational coordinate rotates $(-\theta )$ around the $e_3$ axis. Here $\varphi$ is azimuthal angle.

To define:

where $T_{21}$ is called transformation matrix from coordinate 2 to coordinate 1. Since all of the coordinate systems are orthogonal, the transformation matrices are all orthogonal matrices, that is to say:

Two basic transformation matrices are given by:

where, $\phi =\varphi -90°\text{.}$ Coordinate transformation between overset sub domain grids is implemented by combining the above two basic matrices.

3.3. Discrete method

By using the finite volume method for Eq. (6) on each control volume $(i,j,k)$, ordinary differential equations are obtained as the form:

with $f=i-1/2$, $i+1/2$, $j-1/2$, $j+1/2$, $k-1/2$, $k+1/2$.

For spatial discretization, second-order central differencing scheme is used to solve the viscous fluxes. Roe Riemann solver as a flux difference splitting scheme is adopted to approximate the inviscid fluxes. At the volume interface ($i+1/2,j,k$), Roe scheme is given by the expression [26]:

where the subscript ($L$, $R$) stands for the left and right physical variables across the interface ($i+1/2,j,k$). $\tilde{\mathbf{A}}$ denotes the Roe’s average matrix. It is worth noting that the eigenvalues of the Roe average matrix $\tilde{\mathbf{A}}$ are:

where $\tilde{q}=\tilde{\mathbf{q}}\cdot \mathbf{n}=(\tilde{u},\tilde{v},\tilde{w})\cdot (n_x,n_y,n_z)$, $q_g={\mathbf{q}}_g\cdot \mathbf{n}$, and $q_b={\mathbf{q}}_b\cdot \mathbf{n}$. Entropy fix is required to avoid non-physical solution. We adopt Harten’s entropy fix [27].

To improve accuracy, fifth-order WENO-Z reconstruction [28, 29] is performed for the left and right physical variables in the Roe Riemann solver. In the curvilinear coordinate system with uniform distribution grids, interpolation is conducted dimension by dimension. Left state extrapolated primary variable ${\mathbf{P}}_{i+1/2}^L$ of fifth order accuracy is expressed by:

where ${\mathbf{V}}_k$ ($k=$ 0, 1, 2) is obtained by extrapolating cell average values in the stencil $S_k=(i-k,i-k+1,i-k+2)$:

By nonlinear WENO-Z weight, ${\omega }_k$ can be written as:

where $d_0=$0.3, $d_1=$ 0.6, $d_3=$0.1. The expression of smooth factor ${\beta }_k$ is given by literature [30]. The introduction of $\epsilon >0$ is to avoid the denominator becoming 0. Traditional choice is $\epsilon =10^{-6}$. However, if $\epsilon$ is too large relative to the computational grid, the solution containing discontinuity could be affected. It is proposed that $\epsilon$ is taken as a function of $\mathrm{\Delta }x$ [29, 31], for example $\epsilon =\mathrm{\Delta }{x}^2$. We choose $\epsilon =O\left(\mathrm{\Delta }{x}^3\right)=$10^-18 because the smallest grid magnitude near body surface is 10^-6. Right state value ${\mathbf{P}}_{i+1/2}^R$ can be derived in a similar way on the stencil $S_k=(i-k,i-k+1,i-k+2)$, $k=$1, 2, 3.

The implicit scheme of dual time stepping is used for the time update [21]. One rotation period is divided into 1200 equal physical time stations. The computation of flow fields are iterated 5 steps at every physical time station, which is called pseudo time computational process. Main iteration and sub iteration correspond to the physical time process and the pseudo time process, respectively. The both processes are implicit, so the time update method is called implicit dual time step method.

Specially, to introduce pseudo time term in Eq. (11):

where $\tau$ is pseudo time, and $t$ is physical time. The physical time derivative is approximated by second-order backward difference of three points, and then we can get implicit discrete equations:

where:

stands for sub iterative residual. Therefore, the solutions of governing equations equal to the steady solutions of Eqs. (18) by updating pseudo time $\tau$.

First order implicit backward difference is performed for the pseudo time derivative of Eq. (18):

where $\mathrm{\Delta }{\mathbf{W}}_{i,j,k}={\mathbf{W}}_{i,j,k}^{m+1}-{\mathbf{W}}_{i,j,k}^m$. With linearization for the inviscid flux $\sum _f^{}(\mathbf{H}\cdot \mathbf{S}{)}_f^{m+1}$, it can be derived as following:

with:

The sub iteration is solved by the use of implicit LU-SGS method [32]:

3.4. Initial conditions and boundary conditions

In order to accelerate convergence, pseudo steady is computed until density residual reaches four magnitude or iteration is carried out 1000 steps. The results are regarded as the initial fields.

Boundary conditions consist of body surface boundary condition for the rotary wings, slotted boundary condition for the C-H topological rotor grids, far field boundary condition for the outmost of the background grid, artificial external boundary and hole boundary conditions for transferring information between the rotor grids and intermedial grid, as well as the intermedial grid and background grid. It is a difficulty for the implementation of artificial external boundaries and hole boundaries when using high order schemes. The variables ${\mathbf{P}}_{i+1/2}^L$ and ${\mathbf{P}}_{i-1/2}^R$ of fifth order WENO scheme use the stencil including cell $I_{i-2}$, $I_{i-1}$, $I_i$, $I_{i+1}$, $I_{i+2}$, thus we propose to adopt three layers of cell fingers of hole boundaries and artificial external boundaries so as to the direct implementation of interpolation without reducing order. The information exchange between overset grids is carried out by trilinear interpolation method.

4.1. Flow field simulations without wind perturbation

Caradonna-Tung two-bladed rotor [20] is selected to be the rotary-wing model. The rotary wings have two same blades which are made up of NACA 0012 sections and have a rectangular planform which has an aspect ratio of 6. To verify the codes and the computational overset grids, subsonic flows around hovering rotor are computed based on the available experimental data. The flight states are a pitch angle 8°, a tip Mach number 0.439, and a Reynolds number 1.98×10⁶. As shown in Fig. 2, the computed surface pressure coefficients are compared with the experimental data at the $Z/R=$ 0.68, $Z/R=$0.80, $Z/R=$0.89 and $Z/R=$0.96 spanwise stations. It can be seen that the numerical results are consistent with the experimental data. Then this solver is further developed to compute the flows of micro rotary wings under low Reynolds number with wind perturbation.

Fig. 2Computed and experimental surface pressure coefficient without wind disturbance

4.2. Aerodynamic response of micro rotary wings under wind disturbance

Computation is carried out for the flight conditions of pitch angle 6°, tip Mach number $M_{tip}=\text{0.03}$, and tip Reynolds number $Re=$ 1×10⁵. As shown in Fig. 3, since the most sensitive influence of atmospheric disturbance comes from that of perpendicular to the direction of incoming flow, we study small perturbation with sinusoidal wind that is parallel to the axial direction of the rotary wings.

The function of atmospheric perturbation in axial direction is expressed by:

Then the grid velocity function for the atmospheric disturbance is:

The blade tip velocity can be written as:

where $\omega$ is angular velocity and $R$ is radius rotation.

Fig. 3The sketch of rotary wings suffering disturbance of wind

a) Rotational motion of rotary wings

b) Atmospheric disturbance

Let $f_g=\varpi /2\pi$, $f=\omega /2\pi$. To define ratio of perturbation amplitude and ratio of perturbation frequency:

the grid velocity function can be written as:

In the paper, the aerodynamic response of the micro rotary wings is studied at different disturbance frequencies with same disturbance amplitude, as well as different disturbance amplitudes with same disturbance frequency. Given the basic ratio of perturbation amplitude ${\epsilon }_0=\mathrm{t}\mathrm{a}\mathrm{n}\left(0.35\right)$, the computational cases are expressed as follows: $\xi =$0.5, 1.0, 2.0, as $\epsilon ={\epsilon }_0$; $\epsilon =\text{0.5}{\epsilon }_0$, $\text{1.5}{\epsilon }_0$, $2{\epsilon }_0$, as $\xi =$ 1.0.

Thrust response under different disturbance frequencies is shown in Fig. 4. We can see that curves of thrust coefficient $C_T$ change with the number $N$ of rotary period under different disturbance frequencies based on the same amplitude ${\epsilon }_0=\mathrm{t}\mathrm{a}\mathrm{n}\left(0.35\right)$, except that Fig. 4(a) is the case without wind disturbance. The number of rotation can be given by $N=(\omega /2\pi )t=ft$. And the expression of thrust coefficient is $C_T=T/\left[\frac{\rho }{2}\left({\omega }^2{R}^2\right)\left(\pi R^2\right)\right]$, where $T$ is the total thrust obtained by integral of pressure on the rotor blades [33]. Therefore, both of $N$ and $C_T$ are nondimensional parameters. It can be seen that all of thrust coefficients tend to stability after second rotational period, which reveals that the computational method has high good robustness. The flows are unsteady despite no disturbance as shown in Fig. 4(a). The varying frequencies of thrust agree with the wave frequencies of disturbance wind, however, the varying amplitudes gradually become stronger with the increasing of disturbance frequencies.

When the ratio of perturbation frequency $\eta =\text{1}$, thrust response under different disturbance amplitudes of velocity is shown in Fig. 5. The thrust amplitude increases as disturbance amplitude of velocity increasing. Meanwhile, the varying frequencies of thrust are almost identical. The results are in accord with the actual conditions.

From the computational results, it can be predicted that the influence of wind perturbation could be controlled by adjusting the rotary wing pitch according to the changing frequency and amplitude of wind speed. Maybe, delay effect should be considered so as to nonlinear control algorithm would be adopted through combining CFD technique.

Fig. 4Thrust response under different disturbance frequencies

a)

b)

c)

d)

Fig. 5Thrust response under different disturbance amplitudes of velocity

Furthermore, the computed surface pressure coefficients are illustrated at the spanwise station $r/R=$0.89 in Fig. 6. Whether the pressure response results of different wind disturbance with the azimuthal angle $\varphi =$180° shown in Fig. 6(a)-(b) or different azimuthal angles with the same wind perturbation of $\epsilon ={\epsilon }_0$, $\eta =$ 1.0 shown in Fig. 6(c)-(d) are consistent with the flight situation. And the results conform to the flight aerodynamic characteristics.

Fig. 6Surface pressure coefficients under wind disturbance

a)

b)

c)

d)