Vortex Loops Based Method for Subsonic Aerodynamic Loads Calculation

The new modification of vortex method is considered. It is based on closed vortex loops that can be considered as the separate vortex “superelement”. It is shown that such approach has number of advantages for vortex wakes evolution simulation. The loops generation algorithm permits to simulate the separation zone due to the selforganization of the loops. In this case, there is no need to determine the location of the vortex sheet separation zones. The developed software package can be used in the calculation of aerodynamic loads on unmanned aerial vehicles, as well as in solving the aeroelasticity problems of structural elements of these aircraft.


Introduction
Vortex elements methods are highly efficient for calculation of unsteady aerodynamic loads acting on aircraft moving at low subsonic speed, because they require significantly less computational resources in comparison with grid-based methods [1,2].In these methods Lagrangian particles for vorticity transportation are used.Such particles are called vortex elements.
There are number of known models of vortex elements for flow simulation around 3D bodies: closed vortex framework, vortex filament, vorton, vortex dipole, vortex fragmenton, etc.Each of them has some advantages and restrictions.In "classical" vortex element methods, for example, in the discrete vortex method [3], vorticity is concentrated in vortex framework segments and it is absent outside them in the flow.However, this method requires number of empirical models to determine the location of the vortex sheet separation lines.In case of vortex methods with separated vortex particles (vortons, vortex blobs, etc.) flow separation zones are formed "naturally" due to vorticity flux approachvortex elements are generated on the whole body surface and these elements self-organization in the flow [4,5].The main part of the vorticity is concentrated in neighbourhood of the vortex elements themselves, however, there is distributed non-zero "additional" vorticity in the flow domain, according to the Helmholtz theorems.Its intensity vanishes far from vortex elements [6].This additional vorticity may cause significant errors at computing aerodynamic loads, acting the body.
Some interesting approach has been proposed in computer graphics for 3D smoke dynamics simulation, which implies vorticity flux simulation through vortex loops generation on the body surface [7].The application of this approach to calculation of unsteady aerodynamic loads acting on bluff bodies showed that the obtained results are in good agreement with known experimental data [8].
The aim of this paper is application of vortex element method modification based on closed vortex loops to calculation of aerodynamic loads, acting on wing of the finite span and to simulation of vortex wake behind wingfuselage layout.

Description of Vortex Loops Based Method
The incompressible 3D flow around fixed rigid body is considered.The flow of the media, which has small viscosity and constant density, described by the continuity equation and momentum equation (Navier -Stokes equation) with no-slip boundary condition on the body surface and perturbation decay boundary condition on infinity.Initial conditions correspond to the circulation-free flow.
In the developed modification of vortex method the vortex wake is considered to consist of K closed vortex filaments -vortex loops of the same circulation Γ.The loop with index k is simulated by the polygonal vortex filament with N k vertices.These vertices are considered as the Lagrangian markers r k,i .We assume the loop legs between two neighboring vertices to be rectilinear line segments, which are described by vectors Δr k,i =r k,i+1 -r k,i .
They induce the velocities Γ ν k,i (r), which can be calculated analytically as vortex fragmentons influences, calculated by using Biot -Savart law and regularized by introducing of small smoothing radius ε [6] and finally give contribution to the velocity field V Ω , generated by vortex wake: Markers move in the flow along the streamlines, i.e. their motion is described by ODE: (1) Numerical integration of system (1) is carried out by using numerical method of the first order of accuracy (explicit Euler method) with constant step Δt.Initial conditions for markers positions in (1) are parameters of the loops at the time of their generation on the body surface.
It is well-known from physical point of view, that in order to take into account the presence of the body in the flow, it is possible to replace it with the vortex sheet of unknown intensity, placed on the body surface, which generates the velocity field V γ (r).From mathematical point of view, the velocity V γ (r) potential can be expressed through unknown double layer potential density g(r).This expression coincides with the Biot -Savart law for incompressible flows [3].
For a body having a closed surface, the calculation of the no-sleep boundary condition subject to creation of a vorticity on the whole streamlined surface can be reduced to the no-through condition on body surface r B with unit normal vector n B [9] V n = VB r B = 0, For boundary condition (2) satisfaction linear system should be solved for unknown circulations g j of vortex fragmenton frameworks, placed on M body surface panels.For solving of (3) the regularization variable and additional condition is used: algebraic sum of vortex frameworks circulations should be equal to zero.The summary velocity field V(r) is the superposition of the velocity field, generated by surface panels V γ (r), incident flow velocity V ∞ and velocity field, generated by vorticity inside the flow domain V Ω (r).Finally, the velocity of the i-th marker of the k-th loop has the form Vortex loops generation algorithm consists of three main operations.Firstly, maximal and minimal values of the double layer intensity for all the panels vertices on the body surface g max =max(g j ), and g min =min(g j ), j=1,…,M should be found.Secondary shuold be construct the level lines of the potential density, which correspond to potential values, which differ from one to another on value Γ.These level lines determine the initial shape of the vortex loops.For the given value of the potential density there can be one or more closed level lines, which correspond to separate vortex loops.The number of level lines values is determined as the integer part of expression Then the generated vortex loops become part of vortex wake.The vortex loops, generated on the surface, should be shifted from the surface in normal direction on small constant distance Δ.It should be noted, that the contribution of the vortex loop to the velocity field V Ω is not considered at the time step immediately after its generation, and for velocity field influence of the corresponding panels on the body surface (which give contribution to V γ ) is taken into account.
When vortex elements move in the flow, some marker positions can intersect the body surface, mainly due to numerical errors in velocity field reconstruction and vorticity motion equations integration.In such cases, the loop legs, which intersect the body surface, is replaced with some other loop legs, which lay on the body surface and provide the shortest way.The Dijkstra's algorithm is used for such procedure [10].
At every time step several procedures for smoothing of loops geometry, loop segments length alignment and loops reconnection are used [11].
It is necessary to control the angle values between the neighboring line segments in vortex loops: If for some vortex loop vertex Ψ k,i < φ, where φ is a given constant, then vertex position should be corrected.
After angle correction, when the vortex loops are smooth enough, it is necessary to "rediscretize" them in order to provide nearly the same length h for all the segments.In this procedure cubic spline interpolation is used in order to reconstruct smooth shape of the vortex filaments.
Vortex loops reconnection is carried out according to the following algorithm.For each marker r k,i the nearest marker r l,j is searched ( j=1,…,N l , l=1,…,K ), for which the conditions are satisfied: where μ>0 is some given value.
From all the marker pairs, for which conditions (4) are satisfied, that pair is selected, for which | r l,j -r k,i | has minimal value, and the reconnected is performed for the coppesponding legs.In this procedure either two vortex loops are formed from one, or, vice verca, two loops are merged into one.Then for the newly formed vortex loop (or loops) the reconnection algorithm is repeated.
Pressure field at every time step is being calculated from the positions of loops by using the analog of the Cauchy -Lagrange integral [12].
where A i is the solid angle under which the panel from the point of pressure calculation is visible.Aerodynamic loads can be found by integration of the pressure distribution over the surface of the body.

Simulations Results
To implement the vortex element method modification, the C++ program was developed, which uses MPI technology for parallel computations.Calculations were performed on a cluster with 16 cores (were used 3 GHz Intel I7 processors).The execution time of one calculation did not exceed 1000 seconds.
For calculation of unsteady aerodynamic loads acting on wing of the finite span the model of wing with foil NACA-2217, chord b=0,1 and wingspan elongation 5,0 is considered.The numbers of panels on the streamlined surface are M=5212.Other parameters were chosen as: the smoothing radius of the vortex element ε=0,002, time integration step Δt=0,00219; the constant circulation of the vortex loop Γ=0.00526; nominal length of loop segment h=0,001; loop legs smoothing angle threshold φ =160; loops shifting distance Δ=0,0001; loops reconnection distance μ=0,04; upstream velocity V ∞ =1,0.
The shape of a vortex loops in the wake behind the wing after calculating of 500 time steps is shown in Figure 1.Thin lines denotes the vortex loops.Clear structure of the Prandtl's horseshoe vortex is seen.Note, that this result is obtained without any additional hypotheses about flow separation line position or any equal suggestions.The comparison of calculated aerodynamic loads coefficients with the experimental data is shown on the Figure 2. Points denotes the calculated results and lines denotes the results of experiments [13].
The error in determining coefficients of Cy and Mz does not exceed 8%, which is permissible for engineering calculations.The relative error of Cx is great due to the small value of the drag of the airfoil, but the absolute value of the error is also acceptable for engineering calculations.The reason for the such errors is probably the roughness of the profile model used for the calculations.Here, additional methodological studies are required.For simulation of vortex wake behind wing-fuselage layout the wing model from previous case was supplemented by a blunt cylinder with diameter D=0,2 and elongation L/D=12,75.The numbers of panels on the streamlined surface are M~5500.
The change in the structure of the vortex wake in the flow behind the layout at different angles of attack is shown in Figures 3-5.Each figure consist isometric and top views of wake obtained after 500 time steps.
The figures show the interaction of vortex wakes behind the fuselage and the wing with increasing angle of attack, when, the effects of vortex flux from a smooth streamlined surface begin to appear more and more.

Conclusion
The application of meshfree pure lagrangian vortex element method modification based on closed vortex loops to calculation of aerodynamic loads, acting on wing of the finite span allow to obtain results acceptable for engineering calculations with minimal calculation costs.
The simulation of vortex wake behind wingfuselage layout is shown that developed algorithm can simulate vortex wake evolution when flowing around a complex body geometry without any additional hypotheses about flow separation line position or any equal suggestions.
The developed software package can be used in the calculation of aerodynamic loads on unmanned aerial vehicles, as well as in solving the aeroelasticity problems of structural elements of these aircraft.

Figure 3 .
Figure 3. Vortex wake behind the wing-fuselage layout for angle of attack =10.

Figure 4 .
Figure 4. Vortex wake behind the wing-fuselage layout for angle of attack =20.

Figure 5 .
Figure 5. Vortex wake behind the wing-fuselage layout for angle of attack =30.