Simplified vortex methods to model wake vortex roll-up in real-time simulations for fuel-saving formation flight

One possibility for reducing fuel consumption is to fly in the upwind field of the wake vortex generated by an aircraft that is flying ahead. Migratory birds use this principle. Manually flying an aircraft at the point of optimal fuel reduction is not suited for routine flight operations as the pilot workload is excessively high. Hence, an autopilot function has to carry out this task. For designing the autopilot, a flight mechanical simulation with a wake vortex velocity model is required that has the ability to calculate the vortex-induced velocity fields. This paper contributes to the choice of a real-time simulation method for modelling vortex-induced velocities that the wake vortex of a leading aircraft generates and that the trailing aircraft shall use during fuel-saving formation flight.Two different wake vortex velocity models are introduced and compared during steady, horizontal flight. One model is based on the Lifting Line Method (LLM) and the other on the unsteady Vortex Lattice Method (VLM). Both models are able to calculate the wake vortex roll-up phase for arbitrary lift distributions, whereas the commonly used Single Horseshoe Vortex Model (SHVM) ignores the near-field roll up. The differences in the induced upwind distribution and vortex filament position are analysed for coarse spatial and temporal discretisation that the real-time constraint requires. Despite the more stringent simplifications of LLM, both methods yield similar filament positions and similar velocity fields for the same discretisation of the lifting surfaces. Finally, the influence of the discretisation parameters is discussed and parameter values are recommended for using VLM and LLM in real-time flight simulations.


Introduction
Aviation has to noticeably reduce its environmental impact. Innovative airframe and engine concepts as well as sustainable aviation fuels shall contribute [1]. Also aircraft operations shall be part of the solution [2]. Even though not explicitly mentioned in Refs. [1,2], it is long known that formation flying techniques, inspired by migratory birds, have the potential to significantly save energy. Transport aircraft can use this technique as well. The Applied Vehicle Technology (AVT) Panel of the NATO Science and Technology Organization (STO) has published a comprehensive overview on the state of the art of formation flight and the related research activities including an extensive reference compilation [3]. The overview addresses the history of energysaving formation flight, the physical principles, related technologies and systems, as well as operational aspects.
In a formation, the following aircraft can reduce its thrust and hence its fuel consumption, when flying in the upwind field of a wake vortex that a leading aircraft generates. So, the trailing aircraft takes advantage of wake energy that otherwise is left unused. The position for maximum fuel savings in a vortex-induced velocity field is called sweetspot. For flying in the sweet-spot, the trailing aircraft has to control its position relative to the leader. During long haul flights, where the sweet-spot position has to be maintained over many hours, formation flight needs to be automated. This automation requires new autopilot functions. For the development of such autopilot functions, a flight mechanical model with a realistic wake vortex velocity model is needed to simulate the vortex-induced forces and moments. Vortex methods, as comprehensively described in Ref. [4], are used for generation of vortex velocity fields. The research at the TU Berlin aims at increasing the realism of formation flight simulations by means of high-fidelity, real-time methods. The objective of this paper is to evaluate the fidelity of realtime methods for the calculation of the vortex velocity field.
Wieselberger [6] was the first who described the fundamental physics of formation flight by the vortex pair behind a lifting surface. He calculated the upwind field by Biot Savart's Law. Schlichting [7] explained the resulting energy savings by applying Aerodynamic Wing Theory. Hummel [8] extended this approach to inhomogeneous formations. In 1986, Beukenberg and Hummel demonstrated the theoretically explained energy savings by flight tests for the first time [9,10]. Further flight tests with different aircraft types confirmed fuel savings up to 18%, for example in Ref. [11].
A simple vortex model that is commonly used for simulations of formation flight and wake vortex encounters is the Single Horseshoe Vortex Model (SHVM). It approximates the wake vortex system by a pair of two straight vortices, see Refs. [5] and [12]. This vortex model represents a rolled-up vortex. It is valid from 15 up to 150 wing spans behind an aircraft. Hence, it is not well suited to model the near field of the wake that extends up to 15 wing spans behind the vortex generating aircraft. As formation flight can be in within the near field, modelling the vortex roll-up increases realism.
Sarpkaya [13] comprehensively reviews the computational methods that are developed in fluid mechanics to simulate and describe the characteristics of three-dimensional, unsteady vortical flows. The application of those methods to real-time flight simulator investigations is not straight forward, as the real-time requirements on the computational performance significantly differ from accuracy requirements for fluid dynamic applications.
Fuel-saving formation flights are assumed to be carried out when atmospheric turbulence is so low that decay and deformation of the wake can be neglected for the envisioned separation distances. However, to investigate cockpit procedures and pilot workload in real-time flight simulators, vortex models are necessary for distances, where the vortex roll-up is not completed (10 wing spans are 0.5 NM for an 80 m wing span aircraft like Airbus A380). Whereas the computational-cheap SHVM can be used for formations with separations above 15 wing spans, more accurate methods are needed for closer distances.
For the investigation of vortex effects on a trailing aircraft, very high accurate methods, like Large Eddy Simulations (LES), have been applied, e.g. by Bieniek et al. [14]. There, the vortex flow field was pre-computed for a certain distance behind the generating aircraft and stored in a "box", from which the wind field can be used. If LES is not available, LLM can be used in a similar manner-with some loss of realism. However, such simulation techniques require steady flight conditions. If the formation is manoeuvring, e.g. changing altitude or course, a pragmatic solution to stay with the "box" technique by adapting the box location and orientation in space was proposed by Kaden [15]. A superior approach would be to use methods, like LLM or VLM that can online address the impact of the lift changes of the vortex-generating aircraft and the impact on the wing. However, in real-time flight simulations, available CPU performance constrains the discretisation granularity. The motivation for this paper was to prepare the online approach by comparing LLM and VLM vortex models using very coarse discretisation.
Research on vortex methods for wake vortex roll-up simulations dates back to the 1930s [16]. Numerical problems resulting from this approach and methods permitting their resolution have coined a research area, see Refs. [17][18][19]. In addition to summing up the state of the art, Devoria and Mohseni [19] describe the possibility of chaotic vortex trajectories in roll-up simulations. As a solution to this problem, extensions to the roll-up methods are recommended. The extensions can, for example, rediscretise the vortex elements, see Ref. [20], or introduce specific procedures to regulate the distortion of the Langrangian methods, see Refs. [21,22]. This paper compares two well-known methods for the calculation of the vortex velocity field: the Lifting Line Method (LLM) and an unsteady Vortex Lattice Method (VLM). Extensions of the LLM or VLM, as proposed by Refs. [20][21][22], are not considered here. The focus is on the comparison of the unaugmented LLM and VLM calculation schemes and their results. Only the core radius function is introduced as smoothing parameter. The results of the LLM and VLM formulations are compared to the SHVM.
LLM can be used for all distances behind aircraft. Kaden implemented the LLM code that is used at TU Berlin and in his work. He compared it to the simple SHVM, see Refs. [5,15]. A set of parameters for the LLM was optimised to match the vortex-induced velocities in the far field, where the roll-up is complete, with SHVM. The optimised parameters comprise the number of vortex filaments, the integration step size Δt , the core radius r c and its regularisation function. The calculation consists of two main steps: first, an offline computation of the vortex particle field, and second, the subsequent calculation of the vortex-induced velocities at the follower's position. Using simplifying assumptions, the LLM's second step, important for the formation flight simulations, has real-time capability and is therefore currently used at TU Berlin.
Also at TU Berlin, Loftfield developed an unsteady Vortex Lattice Method (VLM) for flight mechanical investigations at separated air flow as part of the project MoSS [23]. The unsteady VLM enables the calculation of wake vortex roll-up behind manoeuvring aircraft. Less simplifications in the VLM yield more realism compared to the LLM. But, the VLM is computational expensive and the present formulation is not suited for real-time application. Here, the unsteady VLM is used. It is applied to a steady, horizontal cruise flight scenario. A one kilometre long wake vortex wind field is generated by an offline simulation. For the same scenario, the wind field is calculated by the LLM. This is achieved by flying the aircraft through planes and capturing the wake properties at their locations, see Fig. 1. Comparing the results of both methods, the effect of the simplifications made in the LLM is investigated.
For the comparison, the filament positions of the wake vortices are analysed. The computed wind field data, resulting from all filament positions and circulations, are used to identify the wake vortex axes. The axes represent the wake vortex position. This position has to be known for fuel-efficient formation flight. In addition to the positionrelated comparison, the upwind areas in the wake vortices are compared. The vertical wind velocity maximum, hence the maximum tangential velocity near the vortex core, are related to the vortex-induced rolling moments, as Bieniek has shown in Ref. [24]. Because the induced rolling moments are important for the simulation of the following aircraft, the vertical wind velocity maxima are used as a substitute in our comparison.
The paper starts with the description of both wake vortex methods LLM and VLM and the resulting models in Sect. 2 addressing assumptions, differences, advantages and disadvantages. To enable a comparison, Sect. 3 explains the methodology to set up the models for a steady, horizontal flight scenario. Section 4 explains how the wake vortex axes are identified. The identification method is used in Sect. 5 for the comparison of VLM and LLM results. Finally, an analysis of the parameter set for this comparison is conducted in Sect. 6.

Vortex filament methods
This section describes the methodologies behind LLM, VLM and SHVM. SHVM is the reference as it is commonly used in wake vortex simulations. The assumptions made by the methods, the differences, advantages and disadvantages of the LLM and VLM are explained.
LLM and VLM are Lagrangian methods, as described in Refs. [4,25,26], tracing particles in a velocity field. The particles are the edges of vortex filaments, which are shed from the aircraft's lifting surfaces. Using the Biot-Savart law, see Sect. 2.4, each filament induces a velocity on every particle and therefore on every filament. The aircraft's lifting surfaces are modelled by vortex filaments, too. These bound filaments also induce velocities, generating the flow field around the lifting surfaces and influencing the wake vortex roll-up. The mutual induced velocity of all vortex filaments leads to a roll-up of the free filaments in the wake, Fig. 1 Concept of the wake vortex roll-up simulation using grid planes according to Ref. [5] generating a rotational velocity field behind the left and right wing of a leading aircraft. The induced velocity from the opposing side induces the typical descent of the wake vortex. In addition to calculating the induced velocity by the bound and free filaments at the vortex filaments' edges, the induced velocity can be evaluated at arbitrary points to compute the induced velocity fields.
Evaluation planes are placed behind an aircraft to calculate and visualise the velocities induced by the vortex sheet in grids. The evaluation planes can be arbitrarily placed. They consist of a grid of particles. At the particle positions, the local velocity is evaluated using the same operations as for the free vortex filaments in the wake roll-up calculation. The induced 3D velocity field is calculated from multiple evaluation planes stacked in longitudinal x a -direction.
Often, an elliptical lift distribution on the wing and tailplane is assumed for LLMs, see Refs. [5,12,27]. In this paper, the regional airliner VFW 614 that is modelled at TU Berlin's research flight simulator is used as example. Detailed geometric and aerodynamic data are available, such as the exact geometry, wing twist and airfoil profiles. Using those data, the elliptical distribution that was used in previous projects is replaced by a lift distribution that is calculated by the VLM. Thus, the wind field computation starts with the VLM calculating the circulation at the lifting surfaces (at the bound filaments). The following steps for the wind field generation are equivalent in the VLM and LLM. These steps are: computing the wake vortex roll-up, calculating the induced velocities at all particles in the evaluation planes and saving the 3D wind field data.

Single horseshoe vortex model
The Single Horseshoe Vortex Model (SHVM) satisfyingly represents the wake vortex velocity field after roll-up and it is commonly utilised for that purpose, see Ref. [5]. Here, it is used as reference for the LLM and VLM.
The SHVM consists of a bounded vortex at the wing (its influence is commonly ignored) and two straight, free counter-rotating vortices that extend behind the generating aircraft to infinity. The circulation of the two straight vortex lines is determined by the Kutta-Joukowski law using the actual parameters of the leading aircraft. For a certain vortex age, vortex decay can be considered by reducing the circulation. In In the following, the global core radius of the SHVM has to be distinguished from the local core radius used to compute the velocities induced by the discrete vortex filaments in LLM and VLM. Figure 2 shows the resulting induced vertical velocities for the elliptical circulation distribution and for the distribution that is calculated by the VLM (for discretisation 1, see Sect. 3). In the VLM calculation, which is used here, the values for aircraft mass, altitude and the airspeed are the same as in [5] and yield Γ = 137.78 m 2 ∕s and a vortex spacing b 0 = 13.88 m . The global core radius r c is the same.

Lifting line method
The LLM calculates the roll-up in several steps that are described in [5]. The wake vortex roll-up is calculated using multiple calculation planes for the wake vortex roll-up calculation, whereas the evaluation planes describe the 3D wind field and define the vortex-induced wind field for an aircraft flying through the planes. Both types of planes are stacked in negative x a -direction behind the aircraft in flight direction, see Fig. 1. The LLM represents lifting surfaces by lifting lines. These lines hold a certain circulation distribution Γ(y) . The local circulation generally increases from a lifting surface's tip to its root. Following the Kutta-Joukowski theorem, the lift is proportional to the circulation Γ , air density and airspeed V [28] The axis of the lifting lines are bound to the aircraft's quarter chord of the aerodynamic surfaces. Here, the bound vortex filaments originate. The continuous circulation distribution Γ(y) is replaced by the sum of the discrete circulations ΔΓ i of the filaments with stepwise constant strength, see Fig. 3. Each discrete step in the bound circulation distribution produces a free vortex filament of the circulation ΔΓ i , see [29]. In other words, the discrete changes in the lifting The origin of the aerodynamic coordinate system (index a) lies in the wing's aerodynamic mean chord. The calculation planes span in y a -and z a -direction. The roll-up of the free vortex filaments is calculated for each plane. The time step size Δt is defined by the distance between those planes, Δt = Δx a ∕V . One can either assume the aircraft to travel through the planes or assume that the planes are following the aircraft in flight. For the previously introduced evaluation planes, both points of view are equivalent as the resultant wake age at each evaluation plane is the same both ways when flying with a constant speed V. The calculation planes are iteratively solved in the LLM to calculate the wake shape. So, visualising the calculation planes behind the aircraft and extending the wake shape in each time step to the next gate is preferred. Another calculation plane is added each time step to extend the wake shape. The iterative calculation starts from the position of the lifting lines and employs the forward Euler integration method. In each calculation plane, the vortex filaments are assumed to be straight and aligned with the x a -axis. By this approximation, the induced velocity in x a -direction vanishes and the roll-up is simplified, see Ref. [5]. So, the LLM only uses the preceding calculation plane that contains the last position of the free filaments to compute the next position of the free filaments. The influence of the bound filaments is considered, but it vanishes downstream. By using this approximation instead of calculating the induced velocities from streaklines, the algorithm time complexity of the calculation decreases.
As a result of this scheme, the runtime of the LLM scales linear with the number of computation planes. Thus, it scales linear with the wake vortex length. The runtime scales quadratically with the number of vortex filaments, as all vortex filaments interact. However, a disadvantage of LLM is, that the lifting line's circulation distribution has to be prescribed, considering the necessary total lift of the aircraft. This can either be done by assuming an elliptical distribution or, as in this paper, by calculating the distribution for example with the VLM.

Vortex lattice method
The unsteady Vortex Lattice Method (VLM) uses a lattice of vortex filaments to represent the lifting surfaces and the wake, as Fig. 4 shows. Each lattice consists of individual vortex rings, each vortex ring comprises four vortex filaments.
Panels model the aircraft's lifting surface geometry, as described in Sect. 3. Based on the panel position and the convention to place the vortex rings' leading segments at the panels' quarter chord [28], the vortex lattice geometry is derived accordingly, see Fig. 4. Each vortex ring has the same oriented circulation Γ ij (circulation of panel i at section j) for all four filaments of the ring, see Ref. [28]. The vortex rings placed on the panel geometry are bound to that geometry. Their circulation Γ ij is calculated so a no-flowthrough condition is reached at the collocation points. The collocation points are placed at the centre of the vortex rings, therefore at three quarter of the panels chord. The condition that no flow passes through surfaces causes a tangential airflow over the lifting surfaces.
The free vortex rings representing the wake are generated at the lifting surfaces' trailing edges and move with their individual local velocity. Obeying the Helmholtz theorems and the conservation of circulation [30], the circulation of the free vortex rings is bound to the filaments Vortex lattice formed by four filaments per vortex ring, adapted from [28] and [23] forming the rings and the free rings circulation is constant over time. In each VLM calculation time step, the induced velocities at the free vortex ring corner points are evaluated and the points are moved according to the induced velocity and airspeed. This procedure yields the wake vortex roll-up in the VLM. The circulation Γ ij of the vortex filaments at the trailing edges of wing and HTP is related to the lift of the respective section, as can be seen in Fig. 4. As each vortex ring circulation is oriented clockwise, the resultant circulation at the intersection of rings must be calculated. This resultant circulation is used to calculate the lift and pitching moment in the VLM. As a consequence of this principle, the circulation of resultant vortex filaments with spanwise orientation adds up downstream. The bound circulation at the trailing edge, Γ 3,1 in Fig. 4, is related to the circulation and the lift of wing Sect. 1.
The trailing-edge Kutta condition Γ TE = 0 , is satisfied by generating a free wake vortex ring with the same circulation, for example Γ 3,1 = Γ w,1,1 in Fig. 4, so they are cancelling each other out [28].
The same superposition principle for the vortex ring circulation applies in the other direction, too. This results in chordwise oriented filaments reflecting the spanwise change in circulation. Particularly, the circulation of those filaments in the first wake row defines the spanwise change in total circulation. This relation is in analogy to the circulation of free filaments in the LLM. The spanwise distance between filaments at the trailing edge is Δy and the separation Δx of the free vortex rings is defined by the time step size Δt , The VLM's roll-up calculation is build upon less restrictive assumptions than the LLM's. For instance, free vortex filaments are not coupled to predefined computation planes. The wake vortex roll-up is calculated by evaluating the induced velocity from all vortex filaments at all vortex ring corner points (wake particles). This results in a quadratic time complexity for the VLM algorithm with respect to the number of bound vortex filaments at the wing and tailplane in spanwise direction and to the number of the free vortex filaments that depends on the simulated vortex lengths. Compared to the VLM, LLM has linear algorithm time complexity with respect to the wake length. Besides the quadratic scaling with wake length, the VLM evaluates all vortex rings while calculating each wake particle's movement in contrast to LLM that evaluates only the latest evaluation plane and the bound vortices. Consequently, computations with the VLM are significantly slower. Another disadvantage is that the simulated wake shape has to be longer than the length for which the wake vortex induced velocities are calculated. Otherwise, a starting vortex might be included in the results that leads to induced velocities not resembling the wake vortex velocity field.
The more complex VLM calculation has advantages: the method derives the circulation distribution from the geometry and actual flight condition as it simulates the tangential flow condition at the lifting surfaces. So, the roll-up calculation can be unsteady because the trailingedge condition is evaluated at each time step. In this way it is possible to include the effect of changing lift distributions during manoeuvres in the circulation of the wake vortex filaments and eventually in the wake vortex roll-up.

Induced velocities
In both methods, the induced velocities are calculated using the Biot-Savart law. The calculation is based on the vortex filament circulations and the distance to a point P, defined by the position vector P at which the induced velocity ind shall be determined as Fig. 5 shows. The velocity induced by a segment j of the filament i at point P is with the distance vector and the circulation Γ i of the vortex filament d .
As stated, for the LLM each vortex filament d is a part of the lifting lines or a line normal to the active calculation plane. In the VLM, each vortex ring comprises four filaments d and all rings are recalculated in every time step. The different application of the Biot-Savart law in both methods results in more calculations per time step for the VLM. As the goal is to compare the methods, both methods use the same function to regularise the induced velocity in the vicinity of a vortex filament, the low-order algebraic regularisation. A review summary of 2D and 3D regularisations for the Biot-Savart equation is given by Ref. [31]. The transcriptions of Eq. (2.2) used in this paper for computation that include regularisation in the Biot-Savart formula are given in [5] for the LLM and for the VLM in [28]. Regularisation functions introduce a viscid core of radius r c around the vortex filaments. From a vortex filament to a distance of r c , the induced velocities increase. At distances larger than r c , the induced velocities decrease and approach the velocities induced by a potential vortex. The choice of r c influences the result of the wake roll-up calculations and thus the resulting vortex-induced flow field [5]. The smaller r c is chosen, the nearer to the vortex axis the transition to the velocity of the potential vortex happens, thus generating higher induced peak velocities.
The wake vortex effect on the trailing aircraft is calculated in two steps: (1) the calculation of the vortex particle field during vortex roll-up, which is currently done offline and (2) the subsequent calculation of the vortex-induced velocities in real-time. The maximum upwind velocity and an upwind area after roll-up are essential for formation flight simulations. The comparison requires that LLM and VLM yield an upwind area outside the core and after rollup that is as similar as possible to the SHVM-despite the relatively coarse gridding that has to be used to comply with computing-time constraints. To preserve the relationship between the wake vortex roll-up and the vortexinduced velocities, the same local core radius must be used for both calculation steps. Kaden investigated the impact of the core size on the roll-up for core radii between 0.015 b and 0.035 b in Ref. [5] and selected r c = 0.025 b for his investigations. Here, r c = 0.02 b is used, see Sect. 3.

Geometry and parameter definition
The circulation distribution in the VLM results from the lifting surfaces geometry and the no-flow-through condition, as described in Sect. 2.3. The lifting surface geometry is modelled by panels and the aircraft is trimmed for steady horizontal flight as defined in Table 2.
The geometry and profile data of the wing, horizontal and vertical tailplane of the VFW 614 are extracted from the documentation of this short haul airliner. The lifting surfaces geometry is modelled in OpenVSP [32] and exported to Matlab. Three models with different discretisation are generated, D1, D2, and D3, see Table 1. D1, D2 and D3 are used in the VLM computations. D1 with the highest resolution is shown in Fig. 6-for the Open-VSP representation-and in Fig. 7 for the Matlab model used for the VLM computations. Additionally to D1, D2  and D3, a finer discretisation for the LLM (DF) has been generated by interpolation of D1. As Fig. 7 shows, the fuselage is neglected in the VLM calculations. Since the flight condition is a steady and horizontal flight without sideslip or rudder input and therefore symmetrical, no forces act on the vertical stabiliser. So, the circulation of the vertical stabiliser is negligible. That is no serious restriction, as the rudder is seldomly used in cruise in airline operations. Therefore, it is not relevant for the investigated formation flight scenario.
In Ref. [5], Kaden initialised the LLM with an elliptical lift distribution. He used a circulation of Γ 0,W = 104.94 m 2 ∕s for the wing and Γ 0,HTP = 22.63 m 2 ∕s for the horizontal tailplane (HTP). The aircraft's centre of gravity was in an aft position to achieve the same lift distribution on wing and tail plane as in Kaden's study [5]. Then, the HTP with zero elevator deflection was used to trim the pitching moment.
This trim condition is used for all three discretisation cases that are applied to compare VLM and LLM. The process starts with loading the geometries D1 to D3 into the VLM computer program and calculating the roll-up with the parameters of Table 2.
The VLM program calculates the circulation of the vortex filaments at the trailing edges. This circulation corresponds to the total lift distribution over the lifting surfaces. Figure 8 shows the circulation distribution for the three discretisation cases. The geometry yields a non-elliptical lift distribution.
The root circulation at the HTP is close to the given reference values [5] for the elliptical loading in this flight condition, whereas the wing root circulation is higher as a result of the non-elliptical loading.
The differences between the SHVM resulting from the actual circulation distribution on the wing and HTP (D1) and the classical SHVM with an elliptical distribution and only considering the wing is shown in Fig. 2.
The SHVM from the actual circulation distribution induces higher velocities and a smaller vortex spacing due to higher wing root circulation and the positive lift on the HTP. Whereas the VLM computes the spanwise circulation distribution from the panel geometry, the circulation distribution along the lifting lines has to be predefined for the LLM. To achieve comparable conditions for the VLM and LLM for the three discretisation cases, the LLM is initialised with the VLM's circulation distribution of the corresponding discretisation. Thus, the VLM defines the spanwise distribution of the LLM's circulation.
The chordwise positioning of the lifting lines is investigated as follows. For three potential positions, LLM and VLM results are compared with the wake vortex filament positions as criterion. Position 1 is a straight line parallel to the y-axis at x = 1∕4 as shown in Fig. 3. It approximates the lifting surfaces with a line while neglecting the sweep. Position 2 also is at the 1/4 line but in addition the sweep is modelled. Position 3 of the lifting line is at the swept lifting surfaces' trailing edges. The last option results in the best fit of the VLM's and LLM's wake filament positions. For that reason it is selected to chordwise position the lifting lines. Thus, by adapting the lifting lines position and using the VLM's circulation distribution, the parameters of the LLM are tuned for the comparison of the roll-up.
Another criterion for the parameter selection is to produce wind velocities in the far field comparable to that of the SHVM with the actual circulation distribution in Fig. 2. References [22] and [33] state that the cores of the vortex filaments should be overlapping throughout the simulation. The overlap ratio K y is the ratio between vortex spacing and core radius. It can be approximated by K y = b∕(r c ⋅ N) , where N stands for the number of spanwise vortex segments along the wing. A core radius was selected, for which all discretisation cases satisfied the criteria of maximum upwind velocity and upwind area best, yielding r c = 0.02 b = 0.43 m . Only D1 and DF have  overlap ( K y < 1 ). In x-direction, D1, D2 and D3 do not achieve an overlap, as K x = Δx∕r c = 2.8 m∕0.43 m = 6.5 . This unfavourable ratio is used under the assumption that the impact on the accuracy is low, if steady or almost steady flight conditions are investigated. For the reference case DF ( Δx = 0.1 m , Δt = 0.0007 s ), the overlap ratio K x = Δx∕r c = 0.1 m∕0.43 m = 0.23 is adequate. The escape of filaments from the spiral can be avoided by reducing the step size and thus the distance each filament travels in between time steps. Reducing the time step size also increases the vortex core overlap.

Vortex identification method
Introducing regularisation to the potential vortex solution in the Biot-Savart equation leads to rotation in the velocity field, see Sect. 2.5. This rotation is bound to the proximity of the vortex cores. Overlapping cores result in common rotational centre [34]. As the cores of the neighbouring filaments are overlapping during wake vortex roll-up, a common rotational centre evolves during vortex roll-up. The rotation in the merging cores allows the identification of the wake vortex axes. Given a monotonously increasing circulation from the tips to the middle of the lifting surfaces, all free vortices on the same side of the aircraft rotate in the same direction. The rolled-up wake vortex core positions for the left and right can thus be identified by calculating the position of two (positive and negative) extrema in the vorticity, see also Ref. [35]. The vorticity is calculated as the curl of the v a -and w a -velocity field in the evaluation planes, neglecting a small u-component in the VLM. The determination of the cores of the rolled-up wake vortices over all evaluation planes yields the trajectory of the wake vortex axes. Furthermore, the wake vortex span can be calculated. Gerz et al. [35] determine the distance between the rolled-up cores by For an elliptical lift distribution, s v equals ∕4 . The root circulation Γ 0 at the middle of a lifting surface is calculated from the mass m to be lifted and the gravitational acceleration g. The root circulation Γ 0 is equal to the lifting surfaces contribution to the total circulation of the wake vortex. Equations (4.1) and (4.2) determine how lift, rolled-up vortex core span and the total circulation interact.

Comparison of filament position and vortex axes position
The differences between the VLM and LLM results are compared by evaluating the position of all free vortex filaments. For each time step, the y-and z-position of the i-th filament in the LLM is compared to the i-th filament of the VLM. This procedure is applied to each lifting surface. As the flight condition is symmetric, only comparing one side of the wake shape (characterised by the filament positions) is sufficient. Figure 9 illustrates this comparison for D2, ten spans downstream-equivalent to about 1.5 s behind the leading aircraft.
On the left side of Fig. 9, only the filament positions from the LLM are visualised and on the right those from the VLM.
The mean value and the standard deviation of the difference in the vortex filament position are illustrated in Figs. 10, 11 and 12. For n given filaments per lifting surface side, the mean value of the spanwise filament position difference Δ between LLM and VLM at a filament age t is The standard deviation of the spanwise positions is calculated by The same formulae apply for the analysis of the z-locations.
For the finest resolution of the aircraft's geometry, D1, the standard deviation increases for filaments older than one second as can be seen in Fig. 10. For D2 and D3, this phenomenon starts later and is less pronounced, see Figs. 11 and 12. In the mean deviations , a large increase cannot be observed. The mean deviations are starting at values up to 0.4 m for the vertical distance of the HTP filaments but are becoming small at large distances behind the aircraft ( < 0.1 m) respectively for older filaments. This leads to the conclusion that the mean filament position does not greatly differ for VLM and LLM calculations.
The standard deviation in y-and z-direction increases with time, as differences in filament positions cause different induced velocities and hence different filament positions in the next time step. Also, finer discretisations cause higher -values, as smaller distances between the filaments lead to higher induced velocities at the neighbouring filaments. But, as Sect. 6 shows, finer discretisations do not lead to significant differences in the upwind profiles. So the standard deviation of the vortex filaments is not a good metric for comparing the fidelity of the methods.
Next, the positions of the wake vortex cores are analysed, defined by the identified wake axis trajectory. The axes are determined with the vortex identification method that is explained in Sect. 4. The comparison is carried out for D1, that resulted in the highest standard deviation for the individual filaments. Figure 13 shows how well the wake trajectories of both methods are matching.
The cores originate from the wingtips at y a = ±10.75 m and form the vortex axes with the typical descent.
Also, the expected convergence to a spanwise distance b 0 is visible. Table 3 lists the y a -positions of the vortex axes at the end of the simulated wind field that is one kilometer behind the aircraft, equivalent to a vortex age and travel time of the aircraft of 7.14 s.
Additionally, the expected value of b 0 ∕2 is calculated by Eq.  Table 3, D3 does not. In Fig. 14, markers are placed on the wake vortex axes one second apart. After 15 b, equivalent to 2.3 s , the roll-up is assumed to be completed. The axes are close to the expected b 0 of 6.94 m at this time. However, they differ after the 2.3 s , rotating about a centre that lies less than 1 m left of the 6.94 m line.

Comparison of induced wind fields
Now, the calculated vertical wind fields of both methods and the vertical wind field of the SHVM are compared. The comparison is carried out at two positions, x a = 10 b and x a = 1 km behind the aircraft. The first position reflects the minimal distance between aircraft in the recent formation flight simulations at TU Berlin [15]. The second distance marks the end of the computed wind field. The position at 1 km is farther downstream than 40 spans, which is a typical distance in formation flight simulations in Ref. [15].
The wind fields are evaluated by an analysis of the upwind difference in the z a -y a -plane and of the vertical velocity at the local height z a of the identified wake vortex axes.
Because the wind field is symmetrical to the x a -z a -plane, only the starboard wind field ( y a > 0 ) is evaluated. Figure 15 shows the upwind at a distance of 10 b behind the leading aircraft, corresponding to a separation of 1.5 s between leader and follower. Figure 16 compares the upwind velocities at the distance of 1 km that corresponds to a flight time of 7.14 s. The right half of the leading aircraft illustrates the scale in Figs. 15 and 16.
The wake vortex parameters of the identified core location and the maximum vertical wind velocity calculated by LLM and VLM are similar as can be seen by the zero crossings of the vertical wind velocity and the adjacent maxima on the right side of Figs. 15 and 16.
The results differ for the three discretisation cases. This means, the calculations are sensitive regarding the choice of spanwise discretisation. The chosen parameter set, see Table 2, is suboptimal as-against expectation and based on a coincidentally favourable position of the strongest wingtip vortices-only the coarsest discretisation, D3, achieves an upwind distribution comparable to the SHVM reference, and only for x a = 10 b , see Fig. 15. The reference distribution in Fig. 2 has no local minima and maxima. Such a smooth profile is not generated with the used parameter set.
The contour plots on the left of Figs. 15 and 16 show that the local differences between VLM and LLM in the y a -z a plane increase with growing distance behind the leader aircraft, in other words for a higher vortex age. The deviations also become more pronounced for finer discretisations, for example D1 in Fig. 16. The deviation contours have their peak values where the wake vortex axes are located, as a comparison with the wake trajectories in Figs. 13 and 14 shows. In Fig. 16, the maximal deviations are vertically located approximately at z a = 12 m for all discretisation cases. This location corresponds to the lateral and vertical wake position in Fig. 14 at the end of the computed field, where the vortex age is t = 7.14 s.
Despite not precisely fitting the SHVM's vertical wind distribution, the curves on the right side of Figs. 15 and 16 exhibit a smooth positive upwind distribution outside the vortex cores. This area is important for formation flight simulations. LLM and VLM calculate similar upwind distributions. So, the maximal vertical velocities that are related to the rolling moment and related to the energy usable in formation flight are also similar. Furthermore, it is expected that the similarity between LLM and VLM also holds for distances larger than 1 km.

Discussion of discretisation parameters
Section 5 compared LLM and VLM for steady horizontal flight. This flight condition is relevant for energy-saving formation flight. At large, VLM and LLM yield similar results for the three discretisation cases (D1, D2 and D3). As the VLM is computationally runtime-expensive, the discretisation parameters are set to relatively coarse values. The values are significantly coarser than in Kaden's LLM study [5].
At x a = 1 km , the vortex has rolled up. This should be visible in the induced vertical velocity. However, none of the three discretisation cases has reached the reference velocities of the SHVM. That is an effect of the discretisation. Although the vortices are rolled-up, that means the vortex filaments have approached to a common axis, the convergence is not sufficient to attain the expected core radius with the corresponding peak velocities. The persistent helical deformation causes the difference in the vortex axis location. Figure 17 shows the expected induced vertical velocity of a rolled-up vortex computed with SHVM (low-order algebraic regularisation) and with the LLM for D1 and the finer discretisation case DF from Sect. 3. The number of bound vortices in DF (128 filaments for the half wing and 64 for the half HTP) is four times higher than for D1 and the time step size Δt is roughly 28 times smaller. Reference [5] used the same discretisation. The induced wind velocities are only computed with LLM as the VLM results are similar. The comparison shows that the discretisation causes differences. The differences are larger for D2 and D3-not shown here. This is not surprising as the overlap ratio in spanwise direction K y is poor for D2 and D3. Therefore, D3 and D2 are not further investigated. But the D1 results deviate as well. Obviously, Δy and Δt are still too large. However, the computer performance requirements for real-time computation prohibit a finer discretisation. This raises the question: are the selected parameters for D1 ( Δy , Δx resp. Δt , r c and the regularisation function) adequate for the task?
To evaluate the suitability of D1 for the formation flight application and to analyse the impact of the discretisation step size, DF is used as reference. In order to assess the impact of r c and the regularisation function, an additional discretisation case was generated with the same Δy and Δt as DF but with r c = 0.015 b and the Gauss regularisation. It is labelled DFmod. Figure 18 compares the induced vertical velocity in spanwise direction for D1, DF, DFmod and the SHVM 10 wing spans behind ( x a = 10 b ). The vertical wind velocities outside the cores ( y a > 10 m ), which are relevant for formation flight, are similar in all cases. Therefore, a finer discretisation than in D1 would not significantly increase the quality of the velocity model. Surprisingly, this can be achieved with a coarse Δt.
The vortex core locations of D1 and DF lie close together at y a = 8 m . They differ by approximately Δy a = 1 m from the theoretical values that the SHVM uses. The vortex rollup is not completed. Within the core ( 7.5 m < y a < 8.5 m ),   16 Upwind comparison for distance x a = 1 km. View from behind the leader aircraft, its contour is overlaid in grey the induced wind velocities are smoother for DF than for D1. However, that effect seems to be less important for formation flight simulations than the difference in amplitudes (D1: 6.8 m∕s and −5.5 m∕s , SHVM: 13 m∕s and −9.9 m∕s ). The difference can be reduced by selecting a regularisation function that provides higher maximum tangential velocities. Using a Gaussian instead of the low-order algebraic regularisation function in the LLM increases the maxima by approximately 20%, which is close to the theoretical value of the SHVM. As the following aircraft in formation flight shall avoid the region near the core, where the induced rolling moment is excessive, the size of the maxima is only relevant for failure case investigations, when this region is unintentionally penetrated. But, a Gaussian regularisation, as well as a high-order algebraic regularisation, can improve the convergence with increasing number of vortex filaments compared to the low-order algebraic regularisation, see Ref. [31]. If required, the LLM results can be further tuned to better fit a reference vertical velocity function. Kaden [5] studied for discretisation DF how variations of the core radius r c between 0.015 and 0.035 b impact the induced vertical velocities.
Heintsch [36] analysed the influence of the step size on vortex filament location and on the induced velocities for a forward Euler method that is also used for the numerical integration in this study. He states, that beginning at a step size above Δt = 0.01 s deviations from the typical and expected spiral form of the vortex roll-up appear. This renders the used step size Δt = 0.02 s for D1, D2 and D3 too large. For DF, Δt = 0.0007 s is appropriate. Heintsch also analysed the impact of the selected integration method by comparing the Euler integration method to second-order Adam-Bashford and fourth-order Runge-Kutta for identical time step sizes. As the results of all methods were similar, Heintsch favoured the Euler method because it is the simplest one. We used the Euler method in this work without investigating more sophisticated integration methods.

Conclusions
This paper contributes to the choice of a real-time simulation method for modelling vortex-induced velocities that are generated by the wake vortex of a leading aircraft and that the trailing aircraft shall use for energy saving in formation flight. Two methods are compared, LLM and VLM. Both are more sophisticated than the simple SHVM that is often used to represent the vortex mid field (> 15 wing spans) but that is not well suited to model the vortex roll-up in the near field.
The comparison is focused on characteristics that are important for fuel-saving formation flight. Hence, the location of the vortex axes and the vortex-induced upwind velocities are compared. They influence the sweet-spot position and the potential energy savings.
The LLM requires the circulation distribution over the wing for the definition of the lifting lines. If it is not available, an elliptical distribution is often used. The VLM calculates the circulation distribution from the lifting surface geometry and the airflow condition. So, it is possible to adapt it to different wing geometries and varying flight conditions. For the comparison, the circulation distribution was computed with the VLM. This circulation distribution significantly differs from the often-assumed elliptical distribution that only considers the main wing. Thus, calculation of the actual distribution increases the realism, especially if the HTP contributes to the overall lift. The circulation distribution calculated by VLM was also used by the LLM in order to eliminate the influence of different distributions.  Both methods use parameters (core radius, discretisation of wing and tailplane, time step size) that are tuned to achieve a velocity field after roll-up that resembles the velocity field of a SHVM. The selected discretisation parameters and the time step size are coarse to allow simulation with both methods on a standard computer. Despite using coarser parameters than in Refs. [18] and [36], the upwind field in the area outside of the cores, which is important for formation flight simulation, has a smooth shape. The limits for the discretisation which still allow achieving the required accuracy for formation flight simulation are strongly coupled to the time step size and the core radius value. Discretisation D1 with 32 filaments on the half of the wing leads to a comparable updraft distribution as the four times higher discretisation DF that additionally uses a 28 times finer integration step size Δt.
The analysis leads to the recommendation to reduce the time step size from 0.02 s below 0.01 s . A selected discretisation can be optimised to better fit a reference velocity profile. This fitting can be achieved by changing the regularisation function and the local vortex core size.
Vortex roll-up computations with LLM and VLM show similar results for the wake vortex axes and for the maximum vertical wind velocity. Deviations are locally observed for the position of the free vortex filaments and for the wind velocities. They increase with finer discretisation as the core radius and the integration time step are kept constant. The mean deviations for the filaments position are very small ( < 0.4 m), standard deviations reach 2 m and deviations in local vertical velocities near the wake cores reach 2 m/s. However, the distance differences are small in relation to the relative distance between the aircraft and the velocity variations will be smoothened as the trailing aircraft flies parallel to the vortex axes. Thus, it can be assumed that the simplifications for the LLM do not significantly impact the accuracy of the velocity field simulation for steady flightwhen compared with the more accurate VLM results.
As the LLM is less computationally expensive than the VLM and delivers similar results for straight trajectories, TU Berlin uses the LLM for vortex field calculations in real-time simulations where the leading aircraft is in a steady cruise flight condition. Even real-time simulations of formation flight manoeuvres have been performed with LLM. Wake vortex velocity fields are simulated along curved trajectories using coordinate transformations as explained in Ref. [15]. However, only VLM can directly model the wake vortex roll-up in manoeuvres as stated in [37]. For real-time simulations of the roll-up, the calculation runtime needs to be significantly reduced.
Future work should consider methods as they are described in [19] and [20], for example to use multiple levels of vortex filament discretisation downstream and thus reducing the number of filaments, to reduce the computation time for both, VLM and LLM, and to achieve finer discretisation for a given computational performance.