Discrete Vortex Cylinders Method for Calculating the Helicopter Rotor-Induced Velocity

A new vortex model of a helicopter rotor with an infinite number of blades is proposed, based on Shaidakov ’ s linear disk theory for calculating inductive speeds at any point in space in the helicopter area. It is proposed to consider the helicopter rotor and the behind vortex column as a system of discrete vortex cylinders. This allows building a matrix of the influence of the vortex system under consideration on any set of points, for example, the calculated points on the rotor itself, on the tail rotor, etc. The model allows calculating inductive velocities at any point near the helicopter using matrix multiplication operation. It is shown that the classical results for the momentum theory remain constant even in the discrete simulation of the helicopter rotor vortex system. The structure of the air flow behind the rotor and the simulation results obtained by the proposed method is compared with the structure of the tip vortices and the results of the blade vortex theory. In addition, the experimental data were compared with the simulation results to verify the correctness of the model under real operating conditions by the helicopter trimming.


Introduction
The character of the load distribution on the disk rotor vortex theory affects induced velocity. In turn, the inductive flow is the most important factor affecting the value of the inductive losses, as well as forces and moments acting on the helicopter's rotor. Therefore, the efforts of many authors are aimed for creating theories and methods for the simplest way to calculate the induced-velocity field, without which it is impossible to calculate the determination of the aerodynamic characteristics of the rotor.
In conditions of low velocities, induced-velocity field is particularly irregular. This leads to significant changes of aerodynamic forces acting on the blade. The blades begin to oscillate with higher amplitudes, causing significant variable tensions inside blades.

Coordinate systems of the vortex cylinder
The properties of a cylindrical vortex surface are considered in detail by Shaidakov [1]. He analytically investigated the properties of a vortex surface that completely covers a beveled vortex surface. Shaidakov studied a vortex surface that starts from a disk plane, has an arbitrary shape in section, and pointed with its one end to infinity. Vortices on the surface are parallel to the base of the cylinder, which lies in the plane of the beginning of the vortex cylinder. When applied to the rotor disk of a helicopter, it is more appropriate to consider not an arbitrary shape of the vortex surface, but a very specific shape, the cross section of which is shown in Figure 1. The disk plane is filled with closed contours formed on two sides by arcs of circles and on the other two sides by radial segments. The number of closed contours depends on the number of calculated points along the blade and the number of points along the azimuth. However, the shape of closed contours remains the same.
To study the velocity field caused by a discrete vortex cylinder, we will select two typical sections of it. One section is located at the beginning plane of the vortex cylinder (lies in the plane of the screw disk); the other is parallel to it and intersects the cylinder at an infinite distance from the first plane. The sections of the cylinder with these planes are conventionally designated 1-1 and 2-2, respectively.
It is convenient to calculate induced velocities in the coordinate systems Oxyz (Ox 1 y 1 z 1 ) shown in Figure 2. The vortex cylinder is tilted from the axis Oz(Oz 1 ) by the inclination angle of the vortex cylinder δ. It is easy to see that the inducedvelocity calculated point is always located in a plane parallel to the disk plane, and the projection of the vortex cylinders on this plane is a disk with a radius equal to the radius of the rotor. The origin of the coordinate system Oxyz (Ox 1 y 1 z 1 ) is always located in the center of this disk.
The right rectangular coordinate system Oxyz is used to record the components of the induced velocity v x , v y , v z . The right-linked coordinate system Ox 1 y 1 z 1 with guide orts e 1 , e 2 , e 3 is used to calculate the components of the induced velocity v x 1 , v y 1 , v z 1 . The axis Ox 1 and Oz 1 coincides with the axes Ox and Oz. The axis Oy 1 is the axis of the cylinder and is inclined to the axis Ox 1 at an angle δ. In both systems, the base plane Oxz (Ox 1 z 1 ) is parallel to the plane of the vortex cylinder and contains the induced-velocity calculated point A. The left skew coordinate Oryψ system is derived from the Ox 1 y 1 z 1 . The left oblique cylindrical coordinate system Oryψ is obtained by moving the origin to a point A. It is used to bring contour integrals to a form that is convenient for integration.
The last two axis systems are most convenient for deriving equations of the vortex cylinder surface. It is enough to know the distance h along the axis Oy 1 from the base plane of the vortex cylinder and the equation of the projection of the cylinder base on the base plane.
The ratio of induced-velocity's components in systems Oxyz and Ox 1 y 1 z 1 is determined by the following dependencies These component interdependences are true for any vector [3]. In the accepted coordinate system Ox 1 y 1 z 1 for an arbitrary pair of vectors a, b we find a vector product с and a scalar product ab Coordinate system of the vortex cylinder with arbitrary form [2].
The modulus of the vector a is calculated by the formula (Eq. (5)) The projection of the vector c on the axis of the cylinder is denoted c 0 . Then, using (Eq. (6)), we will have

Components of the induced-velocity vector at any point in area around the rotor
Let us consider a certain part of a cylindrical vortex surface, with the beginning at the plane of the disk, bounded on two sides by two generatrices and leaving the other side to infinity (Figure 2). The beginning of the generatrices is denoted by the point's m and n. Select a vortex element ds with circulation dΓ at any point M ξ, η, ζ ð Þ in the vortex surface: where γ is the running circulation in the direction of the cylinder generatrices. We will calculate the induced velocities from this element at the pointA x 1 , 0, z 1 ð Þusing the formula of Biot-Savart where l is connecting the point A to the point M vector; ds is a vortex element represented as a vector. The sign of circulation dΓ (or γ) is determined by the direction of the vector ds. A positive value corresponds to a positive direction in the accepted coordinate system. When calculating the inductive effect from the vortex surface, the contour integral is calculated along the vortex lines. The positive direction of the contour traversal in the right coordinate system corresponds to the right rotation, in the left coordinate system-to the left.
We express the vectors included in the formula (Eq. (8)) in terms of affine coordinates ds ¼ dξ e 1 þ dζ e 3 ; (10) The modulus of the vector l is defined by the formula (Eq. (6)) Calculate the components of the inductive velocity vectors from the vortex path ds. For this purpose, it is necessary to define expressions for projections of the vector product ds Â l on the coordinate axis.
After describing the determinants, we will have the following expressions for induced velocities Let us perform integration of formulas (Eqs. (16)-(18)) along the cylinder creators After the conversion, we will have In formulas (Eq. (22)-(24)) we will denote by Using the following formulas Integrals (Eq. (26)) are table integrals and can be denoted by the following: where Using the already known relationships cos converting the received formulas where Integrals were obtained for the first time by Shaidakov [3] and are referred to in this chapter as Shaidakov's integrals. Now it is possible to get expressions for differentials of induced velocity's components on the axis of a coordinate system Ox 1 y 1 z 1 We get the formula for the component of the induced velocity directed along the vortex cylinder generatrices. We use the expression (Eq. (7)) After simple transformations we get To calculate the velocity from a limited width vortex cylindrical surface, we need to take the contour integral from the projection of this surface on the base plane (the integral along the length of the arc s).

Discrete vortex cylinders method
The method of discrete vortex cylinders is based on the linear disk theory of Shaidakov described above. It allows you to calculate the induced velocity from the rotor at any point in the space around the rotor. Consider a main rotor with an infinite number of blades [2]. Imagine a vortex system that descends from the rotor in the form of a vortex column, starting at the plane of the disk and going to infinity. The vortex column is supported by a circle with a radius equal to the radius of the rotor. The angle of vortex column inclination to the disk plane depends on the helicopter forward flight speed and the thrust of the rotor. It is calculated using the Shaidakov formula for the angle of inclination of the vortex column.
Each partial volume of the vortex column can be considered as an elementary column of dipoles with a constant density of circulation. Alternatively consider it as an elementary vortex cylinder of arbitrary shape with a linear circulation of closed vortices γ along the generatrix.
In case of the beveled cylinder filling with dipoles, to calculate the inductive velocity from the entire vortex column, it is necessary to make integral sums from the n number final volumes at n ! ∞ for the limit case. In this case, the area of the base of the cylinder is divided into n number areas dσ 1 , dσ 2 , … , dσ n . The area of the cylinder base is divided into finite regions, and the entire vortex column is divided into infinite volumes.
When filling a vortex column with vortex cylinders, the area of the base is filled with closed contours of a specific shape (Figure 1), and the column is entirely filled with vortex cylinders of linear circulation γ along the generatrix. The generatrices of the vortex cylinders are parallel to the axis of the vortex column and inclined at an angle δ to the plane of the disk (Figure 3).
We propose to consider a vortex column as a collection of a finite number of vortex cylinders resting on the plane of the rotor disk. The plane of the disk is filled with closed vortex contours on two sides by arcs of circles and on the other two sides by radial segments (Figure 4). When the disk is split in this way the size of the contours will depend on the number of calculated points along the blade radius and along the circumference of the rotor disk. In this case, the induced velocities will be calculated at points located at the vortex cylinder's axis, and outside the vortex column at any point other than the vortex cylinders surface. This avoids computational difficulties when calculating contour integrals on the surface of cylindrical columns. This point is indicated by a letter A [3]. We will also follow to this designation.
To calculate the components of the inductive velocity at point A by the method of discrete vortex cylinders we will use Shaidakov's formulas (Eqs. (31)-(33), (35)), derived for a discrete vortex cylinder (Figure 4).
We assume that the circulation along the contour and along the generatrix of each discrete vortex cylinder is constant. Then in formulas (Eq. (36)) it is possible to take the linear circulation γ as an integral and calculate the induced velocities from the vortex cylinder of the unit circulation, integrating along four segments of the contour (Figure 4).
After the obvious transformations, we will have: where l z ¼ ζ À z 1 ð Þ, l x ¼ ξ À x 1 ð Þ (Figure 2). The derivatives dξ ds , dζ ds are the cosine and sine of the angle of inclination of the arc element ds to the axis Ox 1 . Their values do not depend on the position of the point A, but only on the direction of the vector ds in the base coordinate system (Figure 4). This angle is denoted by To calculate integrals along the contour, we will use integrating matrices [3]. For convenience, we will calculate the closed loop integral as the sum of integrals ð Þ dv over four contours. On segments ab and cd the angle τ is constant along the contour and is equal to the azimuthal angles ψ ab and ψ cd þ π, respectively. On segments bc and da, the angle τ changes along the contour length, equal to the angle of inclination of the tangent in the middle of the arc ds to the axis Ox 1 . It depends only on the midpoint of the arc ds azimuthal position.
With the proposed natural partition of the disk into discrete vortex cylinders, the problem of calculating inductive velocities in the disk plane is reduced to multiplying the matrix of influence on the column of linear circulations γ at the calculated points along the rotor disk.
In hovering mode, the main diagonal (i ¼ j) of the influence matrix A ½ in the disk plane is 0.5, and for the influence matrix in a plane far from the disk plane is 1. Non-diagonal elements of the matrix A ½ are close to zero.

The aerodynamic and inertia loads on the blades
The aerodynamic load on the blades is calculated from the known equations of the aerodynamics of a helicopter rotor is the same for main and tail rotors of the helicopter for the spatial movements relative to the longitudinal and transverse axis of the helicopter We add inertial loads [4] to the distributed aerodynamic forces in the blade cross section As a result, we obtain the blade flap equations. The integration over the length of the blade gives us the equation With a trimmed helicopter flight, the flapping of the blades can be represented as a Fourier series. This uniquely determined by the blade flaps angular velocity dβ b =dt and angular acceleration d 2 β b =dt 2 . The parameters of the helicopter state are set by the main and tail rotor control (θ 0 , θ c1 , θ s1 and θ tp ) and the load computation is reduced to the computation of the Fourier series coefficients β 0 , β c1 , β s1 , … , β cn , β sn of the blade flaps β b . The Fourier coefficients are determined from the equation (Eq. (49)) by Newton's method. The number of coefficients in this case should be equal to the number of rotor azimuth steps.

The fuselage aerodynamic loads
The fuselage aerodynamic loads depend on the angle of attack or the angle of sliding is evaluated by wind tunnel or calculation by CFD-method.
The angle of attack The angle of sliding In this case, the coefficients in equations (Eqs. (50) and (51)) are calculated with CFD-method and corrected by the results of the tests in the wind tunnel.

Trimming model
The helicopter trimming equations derived from the helicopter dynamics equations. Equating to zero the linear and angular accelerations of the helicopter is defined as a data source. As a result of the system of equations (Eq. (52)) solving we obtain the vector X ¼ γ, ϑ, θ 0 , θ c1 , θ s1 , θ tp È É T . This vector is calculated by Newton's method.

Results
To illustrate the comparison of result, we distribute the induced velocity in the plane of the rotor disk and at a far distance from it ( Figure 5). This confirms the results of the Momentum theory: the inductive velocities in the plane of the rotor disk are two times less than the inductive velocities at an infinite distance from it [5].
The tip vortices structure of the rotor shown in Figure 6. Comparison of the structure under the main rotor with the blade theory results shows an adequate behavior of the vortex surface in modeling.
Introducing the air flow configuration, we can see which areas of the helicopter are influenced by the induced flow and used in the analysis of information corresponding loads. Figure 7 shows a comparison of the position of blade vortex theory (disorderly line) and disk theory (black line) for horizontal flight. The tip vortex is shown only from one blade, but the influence of all blades is taken into account.
The results of calculating the normal component of induced velocity were compared with experimental data for forward flight from 75 to 180 km/h. Comparison with the experiment gave good results (Figure 8). Experimental and calculated data have got an adequate correlation. The application of the described method consists in forming a matrix of influence of the main and tail rotor for any group of points around the rotor. For example, a matrix of influence of the rotor on the fuselage, the matrix of influence of the main rotor on the tail rotor, the matrix of influence of the main rotor on the stabilizer, the matrix effect of the tail rotor on the main rotor, etc., the dimensions of the matrices depend on the selected number of design points on the stabilizer, fuselage, and tail rotor. In this case, these matrices can be formed taking into account the mutual influence of the main and tail rotors.  For a single arbitrary point, this will be a matrix-string for the main rotor and a matrix-string for the tail rotor. Multiplying the row by the column of known air circulations on the disk, we calculate the induced velocities at the selected point.
The values of the influence matrix elements depend on only the rotor geometric characteristics and the position of the selected point relative to the rotor and the angle of inclination δ of the vortex cylinder. In the process of calculations, the values of linear circulations and their corresponding induced velocities are refined. Influence matrices must be built in advance for a possible range of angles. During the calculation process, interpolate between the calculated matrices by the necessary values δ.
At different flight speeds, induced velocities have different effects on the tail rotor and stabilizer. In the figures, you can see that induced velocity have a special influence on the low flight speeds mode of V ¼ 60 km/h reach 4 m/s (Figures 7 and 9). This is 24% of the flight speed. Therefore it is very important to know and take into account the field of induced velocity's at low flight speeds.
The horizontal and vertical components of induced velocity at a characteristic point in the tail rotor and stabilizer area significantly depend on the flight speed (Figures 7 and 9).
By drawing the configuration of the air flow, we can see which areas of the helicopter are affected by the induced flow. This allows us to correctly configure the loads calculating program. Figure 10 shows the position of the tip vortex according to the blade theory (disordered line) and the disk theory (black line) for forward flight.   Comparison with the blade theory for calculating the position of the tip vortex shows that the method of discrete vortex cylinders gives satisfactory results. The angle of vortex cylinder δ calculated from the disk theory inclination coincides with the angle δ calculated from the blade theory.
As the graphs show, calculating the position of the tip vortex by the method of discrete vortex cylinders gives satisfactory results. Figure 11 shows the results of comparison of the normal component of the average induced velocity calculated in the cross section along the vane theory [3,6] and the discrete vortex cylinder method with experiment [2]. Error experiment for flight speed μ ¼ 0:095 is AE15% and the angle of attack α ¼ À9:2 ∘ .
From these results, we can conclude that the greatest convergence with experimental method of discrete vortex cylinder has got in a cross section X ¼ 0. In some cases, we have strong disagreement with experiment, but similar to the results obtained by the blade theory.

Conclusion
In this chapter, the authors have developed a method of discrete vortex cylinders based on the Shaydakov's disk vortex theory. The capabilities of the discrete vortex cylinder method are demonstrated using a helicopter balancing program based on a "semi-rigid" model of the main and tail rotors. Data from numerical calculations are proposed with experimental data from actual flights. Simulations have been conducted. We found the following: 1. The complex procedure for calculating inductive velocities at any point in the space around the rotor was reduced to the procedure of multiplying the row by column. In this case, the row of influence of discrete vortex cylinders coefficients is multiplied by the column of the corresponding vortex cylinders circulations.
2. The use of pre-calculated influence matrices makes it much easier to calculate circulations.