The Full Multi-wake Vortex Lattice Method: a detached flow model based on Potential Flow Theory

One of the main issues concerning the standard Vortex Lattice Method is its application to partially or fully detached flow conditions, where non-linear aerodynamic characteristics appear as the angle of attack increases and/or the aspect ratio decreases. In order to solve such limitations, a pure numerical approach based entirely on the Vortex Lattice Method concepts has been developed. The so-called steady “Full Multi-wake Vortex Lattice Method” comes from the main hypothesis that each discretized element on the body’s surface detaches their own wakes downstream. The obtained results match for lift, drag and moment coefficients for the entire aspect ratio range configurations (under straight wakes and inviscid assumptions). Future unsteady versions of such a multi-wake approach could improve the current results obtained through Vortex Element Methods (as vortons or isolated vortex filaments).


Introduction
As is well-known, the standard Vortex Lattice Method (single trailing edge wake; VLM) [1,2] has the limitation to be applied successfully only to relatively medium/ high aspect ratio (AR) plates (where wing tip vortex interactions are negligible) and low angles of attack (AoA; α ), where the lift coefficient (CL) physically maintains, more or less, a linear relation as the AoA increases.Although such vorticity-based method has been used extensively for preliminary aircraft design in the last decades, a few improvements have been done in the last years due to the widespread belief that potential flow is limited only to be applied to attached flow conditions, discrediting it to solve more complex flow patterns, as near to stall point where partial flow detachment is present.An extensive review about the VLM history until nowadays is presented in [3], however, references to separated/detached flow research are limited to two old publications, with more than three decades since they were published.The first one is related to practical (not on fundamentals) unsteady aerodynamics [4], while the second to high-speed (transonic and supersonic) aerodynamics [5], in both cases by detaching vorticity only from external edges in delta wings.
Related with flow detachment from external edges, it has been demonstrated that the inclusion of lateral wakes along the side edges (or wing tips) in a rectangular plate allows to simulate the detached flow along such lines and offers better results for higher AoAs and lower AR [6], but maintains the lack of precision for partial or full flow detachment (massively separated) from the surface, as in a bluff shell-body (e.g., round parachute canopy).By following this idea, the inclusion of additional wakes should offer better results for conditions where the flow is considered detached behind the plate.
From the numerical (inviscid potential) point of view, the main hypothesis behind the previous proposal comes from the fact that, ideally, the detached flow must exist at the entire plate, independently of its AoA (due to their sharp edges and zero-viscosity assumption).It means that even for a flat plate at extremely low AoA (e.g., α = 1 deg.), the detached flow is present on the entire plate including their external edges (even along its leading edge; LE).Then, the intensity of such detached wakes must be determined by the instantaneous circulation strength calculated on the surface.For this particular case, the flow detachment could be negligible, but nevertheless exists.On the other hand, from the physical (viscous) point of view, such detached flow assumption should be only valid for the AoA range where viscous effects could be considered as negligible (after viscous-inviscid transition range; moderate-high AoA).
There are a few models and codes already published which try to solve the detached flow condition via panel or vorticity-based methods.One of these is Gersten's non-linear model [7] (based on a theory with the same name), which has the lack that considers a constant circulation distribution only along spanwise and an arbitrary wake alignment ( α/2).As is demonstrated in the validation part of this manuscript, such value must cor- respond to α (the most logical value at infinite downstream after a fully developed wake).The wake alignment has a strong implication for low AR configurations for CL (and thus, for induced drag coefficient; CDi), CD and CM calculation (drag and moment coefficients, respectively).Such issue could be negligible for high AR configurations, according to numerical evidence.
The "nonlinear Vortex Lattice Method" [8] (straight wake version) has the inconvenience that detaches the wakes downstream based on the Lifting Line Theory (LLT), which means from the 1  4 panel chord, not as in VLM (from trailing segments of each bounded vortex ring; thereinafter BVR).This issue is increased due to the fact that the discretization along the chord shown in such publication is extremely low (two or four panels), in consequence, underestimating by far the LE wake detachment.Another lack with this model is that it considers the vortex separation only on lateral edges of the rectangular flat plate, while inner vortices remain embedded.Furthermore, it has the same lack that Gersten's non-linear model (in which it is based), by considering an arbitrary wake alignment with α/2 .Finally, the wake rollup version of this model shows an overes- timated CL from α > 15 deg.(quadrangular flat plate case; AR = 1).
LinAir [9] is a code based on Weissinger's LLT [10] represented by discrete horseshoe vortices, which includes Polhamus' suction analogy [11] to take into account the LE suction effect (Fig. 1).Due to this last approach, it turned out to be uninteresting to consider as a base of the current model development, since it requires calculating suction parameters based on external (experimental) values.
Until now, no numerical vorticity-based model (or method) has been reported yet truly based on VLM concepts, which should offer the advantage that computes the circulation distribution on the entire surface (along chordwise and spanwise directions; fully non-linear in the sense of detached flow), by considering the flow detachment exactly between all discretized bounded elements (including all external edges).The Full Multi-wake Vortex Lattice Method (FMVLM) has been developed to solve such issues.

Problem description
The Potential Flow Theory states that a potential flow must be incompressible (divergence-free): and irrotational (vorticity-free): However, compressible potential flow also can be defined, but not within the classical Potential Flow Theory (PFT).Then, strictly, and, in a broader sense, a potential flow must only fulfill the irrotationality condition (2).
At this point, the term "viscosity" is absent from such definition; thus, the fluid-flow viscosity (ν) value can be considered as zero: Then, an "inviscid potential flow" now is called as "ideal flow" in order to be more specific.However, a viscous-forced potential flow (Kutta-Joukovski "pseudo-real" flow) also can be defined, being this the typical attached flow approach described by the wellaccepted PFT, presented in all references until nowadays.According to the previous proposed idea, a "potential flow" should not be interpreted as a synonym of "ideal flow" as most authors misunderstand.In other words, the ideal flow (particular) is a potential flow (general), but a potential flow (general) is not necessarily the ideal flow (particular).See Fig. 2 for better comprehension.
In order to be clearer, the previous affirmation is explained next, and its discussion constitutes the basis of the current development.It must be remembered that, originally, (1) (3) ν = 0. the PFT was applied by M. Kutta and N. Zhukovski (Joukovski) to justify the generation of total (viscous plus inviscid) lift force on an immersed body (e.g., bidimensional thin flat plate), but not the (total) drag.In any case, only the induced drag is obtained by projection of the lift vector.However, from such approach, the viscosity is being imposed arbitrarily by attaching the flow (vorticity) to the plate, which must be considered as a forced viscous assumption, being a clear contradiction within an inviscid approach: if such flow is inviscid (ideal), then vorticity should be detached from the sharp (zero thickness) LE, independently of the AoA (which clearly does not occur from the Kutta-Joukovski's viewpoint), and finally, obtaining the true ideal flow result (a lower lift caused only by the inviscid contribution, within the low AoA range; viscous regime), which must be independent of the Reynolds number.
According to the previous order of ideas but applied to the well-known d'Alembert's Paradox potential flow problem, where an unphysical result (no drag) is theoretically obtained, a detached potential flow approach could be explored in order to obtain the expected (real) results.It must be considered that, under the current approach (forced attached flow assumption), the streamline patterns around a bidimensional cylinder are hypothetically the same between both inviscid and viscous (infinitely/very large value for viscosity) cases, constituting a clear inconsistency through the commonly accepted theory (see Fig. 3).Thus, such fact does not mean that the current attached flow assumption is wrong by itself, however, it must be treated as a particular case (a forced/imposed viscous version of the potential flow, in fact) of an "extended"/general theory, which should treat all flows as fully detached (true inviscid; the general case), being the fluid viscosity, an external factor to the original definition of a potential flow, as said before, "a" cause to maintain the flow attached to the body's surface; another cause could be due to flow recirculation (with reattachment behind the body), caused by the unsteadydeveloped velocity field, condition that should not be forced a priori through an embedded/attached/simplified vorticity model.
An alternative PFT could be explained through some potential numerical methods.For instance, based on the VLM, it has been demonstrated through the current development that the theoretical aspects behind a new hypothetical PFT remain the same as the original, but only incorporate additional hypotheses for flow detachment.Fortunately, due to the linear superposition principle, the PFT allows to generate more complex flow patterns from simpler ones, which should be the justification to obtain such convenient results through a purely numerical approach.
Continuing with the same order of ideas, through the VLM it had been widely accepted "by default" that flow detachment (and thus, vorticity generation) occurs only from trailing and lateral (wing tips) edges, which seems to be a crude simplification that limits its application range.Such approach avoids that partial or full flow detachment cases can be precisely solved (through spatially precise shedding vorticity), since the current simplified models consider vorticity embedded on the entire plate, being a zero-vorticity assumption just behind it (no flow perturbation), a more than questionable physical behavior within a continuum medium assumption; even the boundary layer must be treated as a kind of detached flow but "damped" due to the fluid viscosity, which exists on the entire surface and not only in some regions (within the viscous regime; at low AoA range).
Physically, the internal vorticity generation (in a rectangular plate case, exemplified for simplification purpose) must be explained as a "vorticity interconnection" due to a transport mechanism (chordwise and spanwise) through the external vorticity generated along their edges; it means, vorticity should not be suddenly cut-off due to the continuity assumption on the shell-body (numerically discretized, in the understanding that conservation of mass is satisfied because the distance between all discretized elements is considered as infinitesimal; in other words, zero-flux across the plate).The same analogy can be applied, for instance, for temperature distribution or structural stresses applied along the plate's edges; such physical properties must be naturally transported to the inner part of the plate, and not confined within certain regions.
Thus, the physical justification to include internal detached vorticity is that, by definition, vorticity is a consequence of shear between layers of materials with different velocities (solid-flow/fluid or flow/fluid-flow/fluid) and can be generated by the interaction of local streams on surface, not only in the surroundings (due to a continuum medium assumption).On the other side, in a recent publication [12], it is shown theoretically that vorticity is generated/created on a vortex sheet (shell-body/plate) due to the difference in flow velocities between their two faces, being a purely inviscid mechanism; the numerical justification to both previous, physical and theoretical approaches, is explained through the current research, presenting both strictly numerical verification and validation cases, for high and low AR flat plates (including a swept-back one), even at sideslip flow conditions.
In the same line, in the 2D case, the Kirchhoff-Rayleigh inviscid separated (perpendicular to plate) flow shows an underestimated drag coefficient.According to the sub-chapter 6.13 Free-streamline theory, and steady jets and cavities in [13], for this case is being solved a cavity flow (with vacuum behind the plate), or in other words, no vorticity is modeled just behind it (Fig. 4); thus, by analogy to 3D, a model based on detaching vorticity only along external separation lines (trailing, laterals, and leading edges of a plate) would also model a three-dimensional cavity case, at least for the post-stall condition.For the pre-stall one, the imposed viscosity assumption should be forcing the flow deflection to round the body, then avoiding a cavity condition, but in any case, the Lagrangian grid distortion (in the scope of a vortex element or particle method) would not be properly solved.

System assembly (thin rectangular flat plate; non-sideslip condition)
As the current FMVLM approach is developed to further be applied to massive-detached flow conditions, a complete arrangement of straight detached wakes which extend downstream is set in order to alleviate the cavity flow issue behind the plate described previously (Fig. 5).
Consequently, the procedure to obtain the FMVLM, based purely on PFT and VLM concepts, is shown below.From the well-known linear system of equations that defines the flow solution in the standard VLM (and panel methods), the instantaneous vortex intensities (or circulation strengths) on the shell-body can be obtained by isolating the corresponding vector, from: Where A is the (general) influence coefficient's matrix (ICM), − → Ŵ is the BVR's (or panel's) circulation strength vector and −−→ RHS is the right-hand side vector (due to the free-stream contribution).Notice that only A is modified during the current model development, by adding the contribution of the additional wakes to the system of equations.

Addition of extra trailing wakes (Multi trailing edges' model; MTE)
The development of this model has been in the most natural way, which means by adding to the standard VLM (single trailing edge wake model) the chordwise vortex wakes (except the plate's LE at this stage, which must be treated in a different manner due to its different nature) called longitudinal or trailing wakes (Fig. 6).It must be remembered that the effect of the trailing edge wake in the standard VLM must be considered by adding the contribution of each wake vortex ring to the particular matrix element of the general (bound plus wake vortex rings; thereinafter WVRs) ICM during the system assembly.By doing this, the Kutta condition is implicitly imposed at the trailing edge of each trailing panel.The same reasoning must be applied for the addition of the remaining trailing wakes, since the Kutta condition is fulfilled automatically.Now, all panels have the trailing edge condition.Because of this, all vortex ring's elements corresponding to the trailing wakes have their own effect on the general ICM.
For instance, for a 2 × 2 flat plate, the general (body plus wake) ICM (based on unitary vortex strengths) yields as: Where the first numerical subscript represents the control point (related to the body panel's numeration) and the second one corresponds to the effect of each bounded (b (5) A = A b11 + A wT 11 A b12 + A wT 12 A b13 + A wT 13 A b14 + A wT 14 A b21 + A wT 21 A b22 + A wT 22 A b23 + A wT 23 A b24 + A wT 24 A b31 + A wT 31 A b32 + A wT 32 A b33 + A wT 33 A b34 + A wT 34 A b41 + A wT 41 A b42 + A wT 42 A b43 + A wT 43 A b44 + A wT 44 .Fig. 6 Body and wake vortex rings for 2 × 2 MTE model subscript) or trailing wake (wT subscript) vortex ring.At this point, the ICM is assembled; as said before, the RHS calculation remains without any modification, compared to the standard VLM procedure.
Note: two additional ways, including by imposing explicitly the Kutta condition, to assembly A were explored during the current development (by taking the difference between upstream and downstream BVRs' circulations for wake strength calculation).However, both additional results were the same between them, but unphysical; for the circulation distribution, in consequence, the lift tends to concentrate towards the trailing edge (instead to LE) and the obtained CL values are out of range (by discarding not well-conditioned matrix issues).

Addition of all lateral wakes (MTE + all lateral)
The next stage in the development of the FMVLM model consists of adding all (internal and external) lateral wakes, by following the same reasoning for the previous group of wakes added, which means to assign a positive unitary circulation value to each WVR's segment, according with the direction criterion (clockwise) of Katz's [2] circulation.Notice that, during the A assembly, each panel must add its trailing (T) edge wake plus their both left (L) and right (R) lateral wakes' influence coefficients.The schematic of this model is shown in Fig. 7.
For the 2 × 2 example, the ICM yields: with Where A Kij represents the sum of all panels' j-wakes with the Kutta (K) condition (trail- ing, left and right; T, L and R) on i-control points.Notice that A wTij comes from the previ- ous model (MTE).

Addition of the plate's LE wake (MTE + all lateral + LE)
The addition of the LE wake detachment improves substantially the CL (also CD and CM) results as can be noticed in the validation part of this manuscript.For this case, unlike for the two previous models, the value for the WVR (and their segments) detached from the LE body panels must be imposed with the opposite sign (their directions are inverted).By following Katz's convention for the direction of the circulation, its value must be −Ŵ (minus one during the system assembly).
Mathematically, such sign inversion in the circulation strength for a LE vortex ring can be explained through the next definition: the WVR's strength ( Ŵ w ) detached between two panels is equal to their upstream ( Ŵ u ) minus downstream ( Ŵ d ) BVRs' circulations [2] (see Fig. 8).( 6) Thus, if upstream BVR is absent ( Ŵ u = 0 ), as in the plate's LE panels, which means that the LE WVR's circulation must be equal to minus its BVR's circulation ( Ŵ w = −1 ).This simple assumption applied to the entire plate's LE allows to take into account the LE flow detachment numerically, as is demonstrated further, improving the numerical results compared to analytical values (and simpler VLM models) for CL and CDi obtained via Jones' linear theory [14,15] for different low AR flat plates.
Then, schematically the imposing of (LE) plate's flow detachment along the LE panels/ BVRs is shown in Fig. 9 ( 2 × 2 example).
By continuing with the 2 × 2 example, the ICM yields: with A LEij being the influence coefficient of the (inverted LE circulation) j-wake (nega- tive sign is implicit) on the i-control point.As can be seen, only LE panels (1 and 2) provide such contribution to the particular matrix's elements (first two columns).Notice that A Kij comes from the previous model (MTE + all lateral).Then, the system assembly for the steady FMVLM (with straight wakes) has been developed, where each BVR corresponds to four WVRs (upstream, downstream and two laterals) or sixteen segments, aligned with the free-stream direction.The upstream WVR must be considered as a shared (merged) wake between upstream and downstream panels, as shown earlier during the development of both simpler models.

Aerodynamic coefficient calculation (Kutta-Joukovski)
After solving the linear system of equations (4) (same as in the VLM) to find the BVR's circulation strengths, the well-known Kutta-Joukovski (KJ) aerodynamic calculation method (see Fig. 10) is applied to find the force vector ( d − → F ) on each bounded vortex segment (on surface): Where ρ ∞ is the flow density, Ŵ and d − → l are the circulation strength and the length, respectively, of the particular vortex segment where force's vector is calculated, and − → U is the local velocity at vortex segment's midpoint.The local velocity at the midpoint of each vortex segment is considered as the sum of the induced velocity from BVR and WVRs plus the free-stream velocity, which means the total flow velocity.
In general, force calculation through KJ does not act on bounded vortex segments which share a KJ-type detached wake (e.g., trailing or lateral wakes for non-sideslip condition), since its circulation strength is nullified ( Ŵ = 0 ) in (11).Such fact avoids that any force can be obtained in a stalled panel configuration (see Fig. 10c), unlike for a single panel lifting surface, where an aerodynamic force (at pre-stall regime; see Fig. 10b) is ( 9) obtained due to the contribution of its LE bounded vortex segment ( Ŵ = 0 ) to force cal- culation.This will be explained in more detail in section 3.2.

Wake vortex rings' circulation strength calculation
In general, all external WVRs (detached from trailing, sides and plate's LE) take its circulation value from the BVR that generates it, in order to impose implicitly the Kutta condition (except LE WVR which inverts its sign) on each separation edge (vortex segment).
As said before, for all internal WVRs (trailing and lateral) the value for its circulation strength must be calculated as the difference in intensity between both (contiguous) BVRs that generate it.In order to generalize this procedure, it is necessary to propose that lower circulation value must be subtracted to the higher.Such operation applies independently of the location of both BVRs.For instance, for an internal trailing WVR (for non-sideslip case), the circulation strength for that WVR is, in most cases (for a positive AoA) the difference between the upstream minus downstream BVRs circulation, as (8) states.For a generalized case (for both internal trailing and lateral wakes), the computational algorithm must compare the two BVRs' circulation strengths that generate such internal WVR to subtract in the correct order.

Bounded vortex segment's vector strength calculation
In general, same as in the standard VLM, the vortex segments corresponding to the panel's trailing edges do not contribute to the aerodynamic load calculation via KJ.This is also true for all edges that have detached wakes, as panel's lateral edges in the FMVLM model (non-sideslip case).The exception is the LE panel; due to the direction of its circulation strength (opposite sign to its BVR circulation) it does not cancel its strength and therefore the Kutta condition cannot be applied along those particular edges.For instance, for a single panel configuration with all wakes detached (see Fig. 10c), this last consideration avoids that no bounded vortex segment contributes to the aerodynamic load calculation, with the force's vector being equal to zero, without providing any lift or drag, an unphysical result for a lifting surface.
Related with the last hypothesis, such LE bounded vortex segment's vector strength cannot be considered as twice ( 2Ŵ ), due to the assumption of a wake discon- tinuity (related to its inverted sign wake) along such edge (see Fig. 11), where strictly only the bounded segment (not the wake, due to its force-free assumption) contributes to aerodynamic coefficients/load calculation.
Schematically, the aerodynamic load calculation via KJ for the FMVLM ( 2 × 2 ) can be represented as follows (Fig. 12), where only the LE bounded vortex segments (with red markers) contribute to aerodynamic coefficients/load calculation (remaining ones, marked in gray, automatically fulfill the Kutta condition).
As can be noticed, for this particular discretization (and non-sideslip condition), Ŵ 1 = Ŵ 2 and Ŵ 3 = Ŵ 4 , thus, the circulation value for both internal-lateral WVRs must be zero, but nevertheless, for higher discretization along the spanwise direction, the internal-lateral WVR exists because there is a difference in circulation between each pair of contiguous BVRs that generate it.At this point, notice that both contiguous internal lateral WVRs have been merged/collapsed into one, obtaining exactly the same numerical result.
Then, for the case of the bounded vortex segments where the aerodynamic coefficient calculation acts (red points), these must take the circulation value that corresponds to their BVR.That is the reason why the aerodynamic load calculation must be performed from the BVR (panel) viewpoint, as if each detached WVR separates (isolates) the surface into pieces (panels).
From the physical point of view, the previous hypothesis has a logical sense.For the FMVLM where all BVRs are isolated, each panel acts like a small flat plate (with sharp edges) with all their four detached wakes (and then, merged with the corresponding ones) and the sum of the contribution of each piece gives the aerodynamic characteristics of the entire geometry.The distance between panels must be considered as infinitesimal between their edges (being the mass flow across the shellbody equal to zero), maintaining the constitution of the entire solid geometry (from Figs. 6, 7, 8, 9, 10, 11 and 12 panels/BVR are shown in an exploded view).

Mesh independence analysis
In order to compare the behavior as number of panels increases between the standard VLM and the FMVLM, a mesh independence analysis for CL and CD to a medium ( AR = 5 ) rectangular flat plate is performed.A low (5 deg.) and a medium (10 deg., which could be considered as the practical limit of the VLM) AoA values are proposed.Unitary values are considered for the chord, magnitude of flow velocity and flow density.The wake length is 40 times the chord [16], which after some tests, continues being valid for the FMVLM.The straight wakes are considered aligned with the AoA (free-stream alignment; α ).The constant discretization (chordwise × spanwise) proposed is: 4 × 8 , 8 × 16 , 16 × 32 , 32 × 64 and 64 × 128 (cosine discretization along one or two directions is avoided due to its high dependence on the vortex core radius [17] selected).

CL convergence
In the next graphs the convergence of CL as the number of panels increases is shown for the standard VLM and the FMVLM.First, the VLM is compared against well-tested Tornado and XFLR5 codes by aligning the wake along the plate/chord (0 deg., according to XFLR5's manual), not along the free-stream direction ( α ).As is demonstrated later, such assumption offers practically the same results (compared with free-stream alignment) for medium or high AR configurations.Thus, in order to generalize the results for any AR configuration, the wake alignment along the free-stream direction continues as originally proposed, which means along the AoA.
As can be noticed from Fig. 13, the FMVLM is slower to reach the convergence than VLM.For this particular case, compared to the number of panels where VLM is considered converged (2048; criterion: less than 1% in difference for the CL), the CL variation for the FMVLM is around 3% (difference with the corresponding finer mesh shown; 8192 panels) for both AoA cases.
Notice that the difference in the converged values must be understood because two different models are compared (linear and non-linear in the sense of flow detachment); thus, there is no need to converge to the VLM's CL value as the number of panels increases.The results for different discretization for the FMVLM are shown in Fig. 14.As mentioned above, an increase in the number of panels for the FMVLM is necessary to obtain comparable results to VLM.
Despite the results are quite similar to the linear VLM model, a CL's non-linear behavior for the FMVLM (slight curvature) appears as the AoA increases.Such characteristic is better observed further for lower AR configurations.

CD convergence
As is well-known, the CDi calculation in panel methods continues being tricky.For that reason, other approaches (e.g., Trefftz plane) have been proposed in the standard VLM in order to obtain a better result for this particular value.However, in the current development, the CD (CDi is not calculated explicitly from the obtained numerical results) is directly obtained as the projection on the wind axis of the total force's vector, calculated via the KJ method on the entire surface.As same as for the CL validation, the results for the CD in the standard VLM have been compared to Tornado and XFLR5 codes.
In the same way as for the CL comparison, the CD is much slower to converge in the FMVLM, being around 5% in difference between the finer mesh and the corresponding one for the VLM where the value is considered practically converged (2048 panels).
As can be seen from Fig. 15, the CD value is more than doubled for the FMVLM compared to the VLM.Apparently, such values correspond to a full non-linear CD, which includes the effect of the detached flow.For that reason, such value cannot be obtained Fig. 14 CL vs α for an AR = 5 rectangular flat plate.Experimental data is obtained from [7] Fig. 15 CD vs number of panels for an AR = 5 flat plate at α = 5 and α = 10 deg through the standard VLM.A brief description of the nature of this non-linear coefficient is exposed in section 2.4 of [18].The validation of this coefficient (also CL and CM) for low AR rectangular flat plates based on small perturbation (linear) theory [14] is shown in the validation part of this manuscript (by using the same computational routine in order to maintain consistency).

Visual results (circulation distribution)
As a result of the postprocessing of the developed open-source code [19], Fig. 16 shows the visualization of circulation strength along each (bounded and detached) vortex filament through the FMVLM for the AR = 5 flat plate at α = 15 deg.case (wake length: 40 times the chord; steady state assumption).
From such visual results, it is clear that concentration of circulation (in consequence, lift force) occurs towards the LE, as physically expected [20].Additional results are shown (see Appendix A) to compare between seven detached flow models (described next) for the α = 10 deg.case.Circulation's ranges are set according to VLM and FMVLM in order to compare adequately between models of the same group (attached and detached vorticity-based).

Quadrangular flat plate (AR=1)
In order to follow the logical evolution from the standard VLM to the FMVLM, seven representative models have been selected for comparison by solving a low AR quadrangular ( AR = 1 ) flat plate (Fig. 17) ( AR = 5 does not show significant difference for CL calculation).As can be noticed, results are totally congruent with expected ones, reaffirming and validating the following methodology.
All models compared are briefly described, and then relevant results are discussed (based on the AR = 1 case): 1. VLM: is the standard VLM with only a trailing edge wake.Its behavior is linear for the CL, and therefore is not capable to solve correctly for low AR configurations.2. VLM plus laterals: is the standard VLM including both lateral external wakes (does not include the LE flow detachment).Its behavior is non-linear for CL and generally For low AoA, CL and CD also fit excellent with Jones' analytical solution (the variation between numerical CD and analytical CDi is less than 1 drag count up to 5 deg.).
Fig. 17 The seven numerical models for comparison; straight wakes (blue lines) extend downstream aligned with the free-stream (physical panels are omitted to avoid line saturation) 7. FMVLM without LE (FMVLM -Leading Edge): is quite similar to the previous model but avoids the wake detachment along the plate's LE (rounded LE).This model is useful to know and compare the LE flow detachment effect on the plate.

Numerical results: CL and CD comparison
The CL analytical value based on Jones' linear theory for low AR configurations ( 0.5 < A < 1 ) is obtained by: Where A is the AR and α is the AoA.In the same way, for the CDi (only valid for very low AoA range due to the linearity of CL), Notice that both the obtained numerical CD and analytical CDi should be calculated in different projections (from total resultant and lift vectors, respectively); however, they must converge for the low AoA range due to the small angle/perturbation assumption.In addition to the results shown in Fig. 18, the complete validation results for CL, CD and CM (for different ARs) are shown at Appendix B, which should be interpreted with sufficient criterion and comprehension of low AR aerodynamics from an ideal (incompressible, irrotational and, among all, inviscid) detached flow perspective.The numerical data corresponding to Fig. 18 is shown in Tables 1 and 2.

Visual results: circulation distribution
From Fig. 19 it is clear that, for both attached (forced viscosity) flow models (VLM plus laterals and OEW), a slight difference exists between them for bounded circulation distribution, while for the detached flow one (FMVLM), such result is quite different, obtaining a similar behavior previously obtained for AR = 5 plate, which means (12) Fig. 18 CL and CD comparison between seven models for an AR = 1 flat plate.Experimental: [21,22] concentrating circulation (thus, lift force) towards LE, as physically expected [14,23].
Notice that such result coincides with the proposed result for a single discretized element/panel, where the resultant aerodynamic force acts independently on each panel's LE, reaffirming the hypothesis that the sum of all individual pieces gives the total aerodynamic characteristics of the plate.

Future improvement: curved surfaces
In order to simplify the current development (see.Fig. 20; physical panels, BVRs and detached wakes are represented in 2D for better comprehension), it was considered that both panels and BVRs coincide each other (see Fig. 20a) in contradistinction to the standard model presented in a well-known state of the art, where these are displaced to one-quarter of the panel length along the chordwise direction (the same idea is applied   to the multi-wake approach; see Fig. 20b and Fig. 20c).The justification to do this simplification comes from the fact that the flow solution is simply displaced on the vertical symmetry plane, being the same in both cases (see Fig. 20a and Fig. 20b).Strictly, such assumption is only valid for flat (not cambered) plates with a constant discretization along the chordwise direction, otherwise, such discretization would not be corresponding exactly to the simulation case that is intended to be solved (see Fig. 20d; here an extra detached wake should be included per pair of BVRs; four instead three).A more detailed explanation with such discretization issues in VLM could be found at comments' section in [24].Despite such simplification, the current development can be applied as well to curved surfaces, by maintaining the hypothesis that each wake detaches exactly between a pair of BVRs (with its corresponding discretization error; a practical way to implement; see Fig. 20e).However, an exact solution could be obtained by detaching two instead one wake per a pair of BVRs along chordwise (see Fig. 20c), increasing the complexity of the model.

Conclusion
During the current numerical model development, it has been demonstrated that both the system assembly and the aerodynamic coefficient calculation for the Full Multiwake Vortex Lattice Method follow a logical sense based entirely on the potential flow and VLM concepts.The developed numerical model avoids the main drawback of the standard VLM and its improvement (VLM plus lateral wakes) by adding both lateral wakes, which means to consider the flow fully attached to the body's surface by imposing viscosity artificially (high Reynolds number assumption).Despite the Full Multi-wake Vortex Lattice Method is computationally more complex and expensive than simpler models, it allows to obtain the correct values for CL, CD and CM (including non-linear effects by flow detachment) not only for high AR but also for low AR configurations.Additional validation cases (not shown due to lack of space) demonstrate that such model modifies the jump in pressure ( p ) distribution and the local lift coefficient (local CL), which causes the drop in the global lift coefficient as the AoA increases (tested for high AR flat plate configuration because reliable data does not exist for low AR).Furthermore, the FMVLM allows obtaining excellent results for CL, not only for non-sideslip but also for sideslip flow conditions (where two leading and two trailing bounded vortex segments exist per panel; see Appendix C).Also, such Fig.20 Bidimensional-simplified schematics for flat and cambered plates from a full multi-wake detached vorticity approach model has been tested not only for quadrangular/rectangular configurations but for a low AR swept-back flat plate, obtaining a satisfactory result.
At this point, the extension to two unsteady versions of this full detached vorticity approach, based on vortex lattice [25] and vorton wakes [26], have been already satisfactorily performed; related future publications will be presented as a continuity of the current research.[17]); vertical dashed lines would be showing the straight wakes limit.Experimental: [28]

Fig. 1
Fig. 1 LinAir's horseshoe vortices arrangement for the first (red lines) and the second (blue lines) row's bounded elements

Fig. 2 Fig. 3
Fig. 2 Potential Flow Theory's flowchart.The new proposed approach (red arrows) strictly fulfills the zero-viscosity definition

Fig. 4 Fig. 5
Fig. 4 Flow past a perpendicular plate (2D) with detachment from their edges and an extended cavity region behind it (black zone)

Fig. 7 Fig. 8
Fig. 7 Body and wake vortex rings for 2 × 2 MTE plus all lateral model

Fig. 19
Fig. 19 Bounded circulation distribution for an AR = 1 flat plate at α = 20 deg.(values are set to visible circulation range)

Fig. 21
Fig. 21 Bounded circulation distribution for an AR = 5 flat plate at α =10 deg.(set to global circulation range)

Table 1
CL vs α for the seven models compared ( AR = 1)

Table 2
CD vs α for the seven models compared ( AR = 1)