The Full Non‑linear Vortex Tube‑Vorton Method: the pre‑stall condition

The present hybrid vortex tube‑vorton method is based entirely on the Full Multi‑wake Vortex Lattice Method (FMVLM) concepts, which means detaching vorticity with pre‑ cise vortex strength and orientation along all separation lines between each discretized element of a shell‑body, including all external edges. Since the classic Vortex Particle Method (VPM) is unstable by itself because it does not conserve the total amount of circulation as time evolves (Kelvin’s circulation theorem), an isolated Vortex (regular‑ ized) Filament Method (VFM) approach is implemented to obtain advection of vorti‑ city, while the induced velocity field is obtained through its corresponding full vorton cloud. Further, a novel vortex squeezing/stretching scheme for such a vortex cylinder‑ sphere approach is proposed based on variation in time for vortex volumes in order to precisely (zero residual) conserve both circulation and vorticity at each time step (for each detached vortex element), while the viscous effect can be accounted for via the Core Spreading Method


Introduction
Eulerian or mesh-based computational fluid dynamics is still the main description used nowadays to approximate the Navier-Stokes (N-S) equations in order to solve a wide range of problems related to fluid-flow interacting with solid walls (e.g., Finite Volume Method).However, in most practical cases, obtaining a solution for a well-resolved simulation can consume large computational resources, mainly due to the fact that such descriptions need to discretize the entire domain, increasing the number of elements as higher precision is required.On the other hand, the Lagrangian or meshless approach has not been completely developed yet, but it is an efficient way to solve dynamics (in a general sense, not only for flows/fluids) precisely because it focuses on discretizing only some (and sufficient) moving elements of interest instead of the entire domain [1][2][3].Such a fact is the reason to obtain, by far (up to double digits), lower computational times through this approach for an equivalent simulation.
Vortex Methods (VMs) are a family of Lagrangian-based numerical methods developed to approximate the N-S equations in their velocity-vorticity form [4,5], but as has been demonstrated in this research, they must be fed by a precise spatially distributed amount of detached vorticity shed from the body's surface in order to solve precisely.Here is the importance of a full detached vorticity model: it provides the correct amount (and orientation) of vorticity to the surroundings (without forcing an attached/viscous flow assumption through a partially misinterpreted Potential Flow Theory), which should be the base model to solve fluid dynamics from a purely meshless approach.Aspects such as vortex strain (stretching), viscous and turbulence effects can also be accounted for through VMs.More complex aspects such as compressibility and thermal effects can be solved in the near future as the development of this new approach evolves.
In the bidimensional case, both the unsteady Discrete Vortex Method (DVM) [6] and the Full Vortex Cloud Method (VCM) [1] (Fig. 1) have demonstrated themselves to be excellent computational tools to help understand the fundamentals of the behavior of the wake dynamics and their consequences in the hydrodynamic load calculations.In the DVM, the detached vorticity occurs only from the trailing edge of the vortex sheet, which represents a small-thickness airfoil (zero-thickness plate) under Thin Airfoil Theory's [7] assumption.On the other hand, in the VCM, the detached vorticity is shed downstream from the entire surface (e.g., a bidimensional cylinder), obtaining satisfactory results even for bluff-body configurations (under a fully detached fluid-flow condition).By analogy to the tridimensional case, the main hypothesis raised in [8] is that vorticity must be detached from the entire surface and not only from some prescribed separation lines as in [9] in order to obtain a precise or even a practical solution (including drag force).Such an approach is followed in the current research presented in this manuscript.
For all the above, VMs and their intrinsically meshless description constitute the most natural way to solve flow and fluid dynamics since their unsteady evolution follows simple and clear rules through well-known equations of motion, as each element of fluid can be represented as a material element that moves and deforms under the action of the remaining ones, including the material belonging to any solid object that naturally disturbs the free-stream.For the same reason, flow (and fluid) detachment should be precisely solved in a straightforward manner, without proposing additional corrections or assumptions as by simply adjusting a boundary layer model or using ad-hoc turbulence models ( κ − ε , κ − ω , etc.) specially developed for solving from an Eulerian description, which invariably, are fed by experimental coefficients, maintaining their semi-empirical Fig. 1 Viscous fluid past a bidimensional circular cylinder through the full vortex cloud method (VCM) with random walk method for viscous diffusion calculation [1] nature (e.g., for the Reynolds-Averaged N-S (RANS) equations) and avoiding being applied to general cases; in the end, all mesh-based approaches (including Large Eddy and Direct Numerical simulations; LES and DNS, respectively) try to approximate fluid dynamics through fixed points in space (including the popular Lattice Boltzmann Method, which is 'meshfree' but full-domain meshgrid; LBM), which in the first instance, and without a deeper analysis or reasoning, seem to be counterintuitive with the etymological definition of "dynamics" that necessarily implies 'movement' .In short, mesh-based approaches seem to complicate the pure fluid dynamics problem by using additional numerical schemes, techniques, methods, corrections, assumptions, and other convenient tricks, for instance, to control both numerical dissipation and numerical viscosity to force convergence and to try to close the N-S equations artificially, which, having said the above, probably do not have an exact solution for general cases from an Eulerian description.Thus, as an alternative, the base method for precisely solving fluid dynamics from a fully meshless approach through VMs is proposed in this manuscript since it is based on purely Lagrangian mechanics and conservation laws.
Since the classic (meshless) Vortex Particle Method (VPM) [3,4] and its reformulated non-turbulent version (rVPM) [10,11] do not conserve the total amount of circulation as their unsteady solutions evolve, both of them tend to numerically diverge as the accumulation error related to the vortex stretching term (based on the temporal variation of circulation) increases.Despite the fact that the rVPM includes new schemes such as vortex volume temporal variation (even for its inviscid version to best conserve circulation) and a numerically precise circulation-vorticity vectorial alignment (divergence relaxation based on [12]), none of these improvements by themselves help to solve the critical instability issue, maintaining its natural unstable characteristic.However, such a development must be considered the basis of the current research since it gives rise to and helps to obtain a deep understanding of general vortex methods in order to explore new proposals and even includes some implementations, such as a LES-based turbulence model, which, according to published results, leads to stabilizing the numerical solution for specific kinds of simulations (e.g., rotor aerodynamics [13]).
In order to solve the implicit numerical limitations described above through the analytical circulation-based vortex stretching calculation, a finite difference vorticity-based scheme is proposed in the current research through the so-called "The Full Non-linear Vortex Tube-Vorton Method" (FTVM).Such an approach ensures that both circulation and vorticity are solved in a precise way since it allows for the modification of the volume of each wake's vortex element at each iteration, obtaining a perfectly numerically stable VM even when it is implemented for low aspect ratio (LAR) configurations at high angles of attack (AoAs).It also allows for the diffusion of viscosity by increasing the vorton's volume (vortex blob) through the Core Spreading Method (CSM).
Unlike both VPM and rVPM, which require additional calculations to approximate a divergence-free vorticity field, the FTVM is based on the isolated Vortex Filament Method (VFM), which offers the great advantage of maintaining both divergence-free velocity and vorticity fields due to the interconnection of all the wake vorticity elements through a grid along the entire simulation in absence of collisions between such vortex elements and the shell-body.If collisions occur as is naturally expected as the unsteady solution evolves, such a divergence-free grid is slightly modified since its interconnection is not enforced but natural, due to the fact that each tube (regularized filament) is isolated from the remaining ones and independently moves in space as it has its own nodes; thus, the current selected scheme ("scheme 2") allows naturally to consider compressibility effects (even if they are minimal as in this case) through a non-divergence-free grid, and it could be applied for future fully compressible flow applications.As far as the calculation of the induced velocity field is concerned, each vortex tube is transformed into its spherical equivalent vorton (hybrid tube-vorton), which offers the advantage of calculating its influence on its surroundings in a spherical manner, not straightly, which is essential to obtaining the results achieved so far.

Methods
The evolution of the wake vorticity vector is obtained by the well-known incompressible N-S equations in their velocity-vorticity form: where the first right-hand side (RHS) term represents the tilting-stretching/squeezing of vorticity (thereinafter vortex stretching) which involves both vorticity and velocity vectors, − → ω and − → u , respectively.In the remaining RHS term, the diffusion of vorticity due to fluid viscosity ( ν ; nu) effect (thereinafter viscous diffusion) is represented.A succinct and clear derivation of such governing equations can be found in [10].
Next, the main contribution in this manuscript is described: a precise way to solve vortex stretching, including viscous diffusion, through a novel growing/shrinking volume of vorticity concept.

A precise vortex stretching and viscous diffusion calculations
In this research, the vortex stretching term is solved through the variation of vorticity [11] (not in circulation as in the classic VPM or rVPM) due to the increasing or decreasing of the vortex core radius, by allowing its precise volume variation at each time step ( t ).Optionally, the viscous diffusion effect can be accounted for via the CSM [14,15].
By definition, the magnitude of circulation (and vectorial circulation) is directly related to the volume of vorticity through: According to Helmholtz's law, the 'circulation' must be conserved between two consecutive time steps, in consequence, along time for each detached vortex element (tube or vorton): According to Eq. ( 4), if vorticity increases due to the elongation of the vortex tube after its advection, its volume must proportionally decrease, and vice versa: Then, the new volume after variation of vorticity is obtained by: (1) ( However, by maintaining a constant volume dV dt = 0 for a vortex tube (cylinder; c) through: as has been proposed in [14], a priori, the conservation of circulation between two consecutive time steps is being violated.According to this, in [11], it is analytically demonstrated that maintaining a constant volume of vorticity, at least in the scope of the classic VPM, such an issue occurs.
The classical mechanism of vortex stretching (by maintaining a constant volume) in a vortex tube is exemplified in Fig. 2, where the vortex core radius ( σ ) must be recalculated after its strain.For the vortex stretching case, the swirl velocity increases, while for the squeezing one, it decreases.

Proposed methodology: dV dt = 0
After advection (through some trajectory integration scheme; see Fig. 3), at t + t , the tilting vector for a vortex tube (regularized vortex filament) is defined as: Where at previous instant t , it was: The time derivative (where capital D denotes the vectorial case) of such vector is then defined by: (5) At this point, it is necessary to calculate the core radius of the vortex tube before its advection (at time t): based on a volume of its corresponding vorton (sphere; s) which is already known from a previous step; it means that at each instant, the vortex tube's (cylinder; c) volume can be represented by a vorton (sphere; s) with the same volume (see Fig. 4): Thus, the time derivative for the vorticity vector [14,16] is defined through: And then, the tube's new vorticity vector at t + t is (for stretching or squeezing, respectively): And according to Eq. ( 5), the vortex tube's new volume is (see Fig. 5): In consequence, the vortex tube's new core radius is: (10) Fig. 3 Advection step for a single isolated vortex tube for one iteration; subindices 0 and 1 denote the previous (at t) and current positions (at t + Δt), respectively Then, the time derivative (first order approximation) for the vortex tube's core radius is obtained by: At this point (optional), the viscous diffusion via the CSM [15] can be added to Eq. ( 16): (15)  where k = 1 is for the Gaussian error function (erf ) or k = 2 is for the 2nd order Gauss- ian regularization function (see Table 1) and ν is the fluid viscosity.Then, the total (invis- cid plus viscous) vortex tube's new core radius is obtained by: Thus, the vortex tube's new volume is obtained, and can be represented by an equivalent sphere with the same volume (only at the same instant; see Fig. 6 for the inviscid and viscous cases, respectively): And the vorton's (called blob in the CSM's scope) new core radius can be obtained through: Then, the vorticity variation due to a growing volume (by viscous diffusion) is recalculated as (for inviscid case V c,t+ t = V total c,t+ t ): (18) High order algebraic And finally, the vorticity vector is modified by (for inviscid case At this point, the vorticity variation through both inviscid and viscous effects can be accounted precisely (zero residual) through the increasing or decreasing of the vorton/blob's volume, ensuring both the conservation of circulation and vorticity for each vortex element between two consecutive time steps and, in consequence, along the entire simulation.A numerical example for two consecutive time steps for a vortex element is shown in Appendix 1 in order to verify this dynamically pure methodology.Further, numerical results between the current volume variation scheme and maintaining constant volumes for all the wake's vortex elements, obtained through the current numerical implementation, are shown in Appendix 2.

Advection of vorticity
After programming the classic VPM [17][18][19][20][21][22][23] and exploring its reformulated non-turbulent version (rVPM) [11] (including variable volumes of vorticity, improvement between circulation and vorticity vector alignment, etc.), some implementations in the current hybrid vortex method (FTVM) are compared against these ones in order to have a reference.In the current numerical implementation, advection of vorticity (including tilting of the vortex tube) is solved through a second-order Adams-Bashforth (2AB) scheme (except at first iteration where a first-order approximation is used): which is applied to both end points (points a 0 and b 0 in Fig. 7) of each vortex (regularized) filament under the influence of the total induced velocity field caused by all the vortex elements (from wake and plate), plus the free-stream velocity.
Figure 7 shows advection of vorticity during the first time step for a single vortex tubevorton, which is located at a prescribed distance ( ε ) over the plate.At t = 0 , the total induced velocity by all the vortex elements is applied on the vortex tube's end points (a 0 and b 0 ); then, at t + t , each point freely moves to its new position (a 1 and b 1 ), by changing its shape (according to the procedure described in subsection 2.1) (including its volume variation between both instants).The same scheme is followed by the remaining detached wake vortex elements at each time step.A way to numerically verify this scheme is to obtain a perfectly interconnected grid along the time (with overlaid nodes for the surrounding wake elements) in absence of collisions of vorticity with the shell-body (see Fig. 26 in Appendix 3).
Note: at this point, in the current computational implementation, the increasing or decreasing of vortex volumes are not restricted (without resetting), except to avoid division by zero for extremely low vorticity values (lower than 1.E-8).(22

Induced velocity field calculation
The total induced velocity field is obtained through the sum of velocities induced by the wake plus the shell-body vortex elements and the free-stream velocity: where for the first two terms, the induced velocity is calculated through a Gaussian error (regularization) function ( g σ ).For both wake ( w ) and shell-body ( b ) calculations, the induced velocity is obtained by: where the tridimensional velocity kernel is defined by: here, − → r is the position vector of the point of interest over which the induced velocity is calculated, while − → r v is the center point's position vector of the corresponding vorton.A 2nd order Gaussian function and a high order algebraic function are also available for testing proposals; both of them show higher scattering compared to the Gaussian error function.
Where ρ = �( − → r − − → r v )�/σ = r/σ is a relation between the linear distance to the point of interest and the vortex core radius.(

Shell-body's induced velocity
The hydrodynamic effect caused by the shell-body to the downstream detached vorticity is considered by transforming the bounded vortex rings, based on vortex segments, into multi-vorton rings as in [22].This scheme is also applied during the influence coefficient matrix (ICM) assembly, and even to obtain the straight wake (or steady-state, if it is assumed that exists) solution, required as an input for the unsteady first iteration.Such a conversion is performed also to the wake's vortex rings, transforming them into multivortons in order to maintain consistency during assembly of the system of equations.From Fig. 8, the spacing vector of the vortons, and thus, the distance ( � − → L v � ) between them, is obtained by: where the default nascent vortex core radius ( σ 0 ) is an input parameter, selected based on the surface mesh discretization in order to maintain an overlapping equal to one (λ = 1) between the closest bounded vortons (before converting each bounded vortex filament/tube to multi-vortons; see Fig. 9a).
Then, the corresponding circulation vector (in some references denoted as − → α ) for each child vorton, which must be recalculated at each iteration, is obtained from the corresponding scalar circulation ( Ŵ ) along its corresponding vortex segment as: In order to maintain a constant vorticity along the original/father regularized vortex filament (tube, represented through a vorton of the same volume) and the new child tubes/vortons after splitting, the volumes of vorticity must decrease according to a new vortex core radius calculation: This last ensures that the total volume (sum of the child vortons' volumes) equals the original (unsplitted) vortex element.Such a calculation is performed only once, and then, the volumes of bounded vortons remain constant along the entire simulation since it depends only on discretization, recalculating only their corresponding circulation strengths at each time step according to its previous temporal solution.
Finally, the default nascent vorton (perpendicular) distance from surface is selected equal to σ 0 ( ε = σ 0 ) to avoid physical crossing/penetration (at t = 0 ) through the plate.Both the longitudinal and lateral positions of each nascent vorton (over the plate) correspond to the midpoint of each discretized bounded vortex tube (see Fig. 4); note that the current scheme only allows detaching one vorton (with its corresponding circulation and volume) per segment/edge at each time step.

Additional considerations
As said before, since the present development is not a VPM approach as [3,4,[9][10][11][17][18][19][20][21][22][23][24], no analytical equations are implemented for vortex stretching through circulation variation ( dŴ/dt = 0 ) and divergence-free vorticity field calculations.Hence, Pedrizzet- ti's divergence relaxation [12] is not implemented because in the VFM the vorticity field is naturally divergence-free, since, as said before, ideally (without collisions with the surface), all the vortex elements remain interconnected through their nodes along the entire simulation, maintaining a well-structured grid.Such fact allows that circulation and vorticity vectors remain aligned between them, avoiding the use of a relaxation scheme.
Related to this last, in the current numerical implementation there are 3 available schemes to treat the crossing/penetration of vorticity through the shell-body (simplified algorithm for flat plates).All of them account for the vortex stretching and viscous diffusion precisely.Next, they are described briefly: 1. Unphysical divergence-free: It allows that vortons cross the shell-body, perfectly maintaining both velocity and vorticity divergence-free fields and, thus, a well-structured wake grid.However, since such behavior is unphysical, it has been discarded to be used during the validation phase (only useful during the verification phase to ensure a perfectly balanced induced velocity field).2. Natural non-divergence-free: It avoids that vortons cross the shell-body, slightly distorting the original divergence-free grid, by relocating them at a prescribed distance (ε) over the plate.This scheme has been selected for the validation phase of the current implementation since it represents the most natural behavior from an isolated VFM viewpoint (this scheme should allow the solution of compressible flow simulations).3. Enforced divergence-free: It avoids that vortons cross the shell-body and automatically corrects the positions of their two corresponding vortex tube nodes and the ones of their surrounding elements, enforcing to maintain the divergence-free grid.
Since it seems that such behavior cannot be physically justified in the understanding that, ideally, all fluids are compressible (strictly, zero compressibility does not naturally exist, even in liquids), this scheme has been discarded for validating the present method.
A numerical comparison between all three previous schemes is shown in Appendix 3.

Aspects to be discussed
One of the main issues in precisely solving flow and fluid dynamics through VM is the Lagrangian grid distortion (LGD), which could be understood as the agglomeration and dispersion of vortex elements in some regions as time evolves [2].Such an issue is even present in bidimensional cases.In the scope of general VMs, some techniques of spatial adaptation have been developed, such as rezoning, remeshing (also known as redistribution), or splitting methods [23] (all are briefly described in [24]), which imply an additional increase in computational cost and code complexity.However, through the current full detached vorticity approach, such an issue should be minimized as the trajectory integration order increases.It means that a 3AB scheme should offer better results than a 2AB one, in the understanding that the full vorton cloud populates, more or less, the surrounding space equally, according to a nascent precise distribution of vorticity, detaching and naturally evolving in time.Note that for a simplified detached vorticity model (e.g., detaching vorticity only from the external edges of a plate), the LGD is under-resolved due to the fact that the distance between most of the detached spherical vortex elements violates by far the ideal overlapping condition (λ = 1), which determines the optimal level of separation between them.On the other hand, such a condition should be improved through a full multiwake approach.Of course, additional research (e.g., parameter variation analysis), even from a theoretical viewpoint, must be performed in this aspect in order to know the precise numerical implications of a tridimensional full vorton cloud (with precise vortex strengths and orientations) in the LGD.However, based on the present numerical results, such an issue seems to be solved, if not completely (which is practically impossible with any numerically discretized method by using a constant time step), at least with sufficient precision for practical purposes.
Another aspect to consider is the lack of a turbulence effect based on Large Eddy Simulation (LES) that has already been satisfactorily implemented in a VM [11].Despite the fact that in the present numerical implementation it has not been included yet, the obtained results show that it doesn't seem to be necessary, at least within the pre-stall range.The justification is that the mixing flow coming from the plate's lateral edges (wingtip vortices) is practically absent, even for the lowest aspect ratio ( AR = 0.5 ) con- figuration analyzed (see Fig. 16).A future improvement to the current numerical method could include such an implementation through a Lagrangian-based LES, among other things, for analyzing a more complex flow condition, such as the post-stall one.

Calculation of hydrodynamic force and pitching moment coefficients
The shell-body's force coefficients, and then the hydrodynamic ones, by projecting to the corresponding wind axes, are calculated through a well-known equation, where the total force on a single discretized element (panel) is obtained through their steady and unsteady force contributions: For the steady part, such a value is obtained the same as in the steady FMVLM [8], that is, by applying the Kutta-Joukowski (KJ) calculation on the leading edge (LE) vortex (LEV) segment of each panel, where ρ ∞ is the flow density, Ŵ and d − → l are the circulation strength and length, respectively, of the particular vortex segment where the force's vector acts, and − → U is the local (total) induced velocity at the vortex segment's midpoint.On the other hand, the unsteady force contribution is obtained by the second term of Eq. (30), where ∂Ŵ/∂t is the time derivative directly obtained by differentiating between the current and previous panel's circulations (same as in the standard unsteady VLM [6]), dA is the panel's area, and − → n the panel's normal vector; additional calculations to obtain the hydrodynamic force coefficients are omitted due to a lack of space in this manuscript (obtained as the projection of the resultant force vector, as said before).
Similarly, the pitching moment coefficient calculation is obtained through its both steady and unsteady contributions: where L is the distance between the perpendicular force application point (panel's lead- ing edge midpoint) and the reference axis of rotation.

Results
Same as in the Unsteady Full (Multi-wake) VLM [25], the hydrodynamic coefficients' calculation (via direct KJ) is highly dependent on the time step selected (low t gives lower values for CL and vice versa).Such an issue seems to be related to the lack of (or excess of ) overlapping between the vortons/blobs [24].However, in the current research it is demonstrated that a t corresponding to a Courant-Friedrichs-Lewy number equal to one ( CFL = 1 ) gives satisfactory results for 3 configurations along the entire AoA (30) range analyzed (within the pre-stall range); a rigorous parameter variation analysis must be performed in order to numerically justify such a selected value (or a similar one).
It is important to keep in mind that as numerical results are obtained for infinitely thin plates (shell-bodies), the effect of the thickness on the hydrodynamic coefficients is neglected; in general terms, as the physical thickness decreases, the CL should increase.On the other hand, as the thickness increases, the CD should increase too, due to the fact that the frontal area increases.Please note that such inferences are made from a purely inviscid viewpoint and are based on a hypothetical reshaping of the plate.However, for the pitching CM, such a consequence is less intuitive to know a priori due to the perpendicular force that generates it, which also depends on the AoA.By the way, for the low AoA range, the numerical viscous contribution should be considered in order to correctly compare with physical experimental data.

Verification and mesh independence analysis
The main input parameters used are shown in Table 2.The remaining ones are selected to be unitary: flow density ( ρ ), free-stream velocity magnitude ( q ∞ ), wake length ( ) and first wake row length factor ( ϕ).
Parameter is related to the total wake length, useful to obtain the steady-state solution (if it is assumed that exists) in order to feed the first iteration by using a large value (e.g., = 40 ), while ϕ helps to determine the first-row vortex elements' positions after each advection step (through a first-order integration scheme for the first wake row), immediately after their creation and locating them over the plate at a prescribed distance ε .Higher values for ϕ throw far downstream the first-row wake vortex elements, deforming sharply the near wake pattern.Due to such parameter is directly related to the time step ( t ), a unitary value ( ϕ = 1 ) can be considered as the nominal one which is selected during the verification and validation phases.

Impulsively started quadrangular flat plate (inviscid; variable volumes)
As the original FMVLM [8] was originally developed to solve a fully detached flow condition (e.g., from α ∼ 12 deg. in AR = 1 configuration), the mesh independence analy- sis is done over such an AoA value (compared to the straight wake one at α = 15 and α = 25 deg.; see Fig. 10); since this last case is out of range for the straight wake assump- tion, a difference between both steady and unsteady numerical results is present.The main input parameters used during this stage are shown in Table 2.
Despite the fact that the middle discretization (10 × 10) offers acceptable results, even for the higher AoA analyzed, the finer discretization (16 × 16) is selected to validate the present method since such a discretization is a common fraction of the base validation (32 × 32) for the steady FMVLM [8], which due to limited computational resources for an unsteady simulation cannot be performed this time, and according to the obtained results, it doesn't seem to be necessary for a low or medium-fidelity simulation.

Inviscid (constant and variable volumes) results
After the verification phase, it has been demonstrated that, at least in the current vortex tube/vorton method (with vortex stretching by varying volumes for each vortex element at each time step), such results are practically the same as those obtained through a pure advection scheme (maintaining constant volumes) as both the number of elements and time steps increase, where the relative error (for all hydrodynamic coefficients) converges to zero.Such behavior allows for obtaining the expected values with less computational effort.This finding could be extended further to more complex flow patterns, such as those expected for the post-stall condition.The corresponding results are shown in Appendix 2. For this subsection, the main parameters are the same as shown in Table 2 for the 16 × 16 discretization case.The pitching moment coefficient (CM) is calculated based on one-quarter chord length for both quadrangular and rectangular plates, while for the remaining plate (swept-back), it is obtained around its apex axis.The remaining parameters are considered to be unitary.
The main input parameters for the three configurations analyzed in this subsection are shown in Table 3.

Validation: CL, CD and CM evolution ( AR = 1 , 16 × 16)
A series of simulations are performed in order to obtain the three main hydrodynamic coefficients from an impulsively started condition ( = 1 ) for a quadrangular flat plate.A comparison with the straight wakes' solution is shown for both lower AoAs (straight

Validation: reached CL, CD and CM ( AR = 1 , 16 × 16)
Despite the slight discrepancies between the experimental data series as the AoA increases, the numerical results match quite well along the entire AoA range analyzed (up to the stall experimental point, located between 35-40 deg.).Note that the current model is only capable of solving within the pre-stall condition, and for this reason, the numerical results after the stall point should be properly interpreted.A further development to solve massively detached flow conditions (including the turbulence effect) should offer comparable results for the post-stall range (after LEV breakdown).
As can be seen from Fig. 12, there are some discrepancies in the averaged CM values for both the lower AoA cases (5 and 15 deg.) compared to the experimental data, showing numerically higher (negative) values.Such a difference could be attributed to the viscous effect on the pitching moment calculation, which is neglected through the present inviscid validation.Further precise viscous simulations could help to improve these results.
In Fig. 13, it is possible to see the well-defined structure of both lateral vortices.Note that the strict symmetry of the wake pattern is due to the fact that there is no numerical randomness as in the bidimensional Full VCM, present in the Random Walk Method (RWM) to account for viscous diffusion.

Validation: CL, CD and CM evolution ( AR = 0.5 , 16 × 8)
A more complex flow pattern case is now analyzed through the inviscid FVTM.The obtained results for the three main hydrodynamic coefficients are shown in Fig. 14.
As in the previous AR case, all hydrodynamic coefficients behave in a very similar way (smoothly converge to a steady value as time evolves).However, in these cases, some discrepancies exist between the averaged and steady-state values.As said before, the time step ( t ) is selected based on the surface discretization size in order to maintain consistency with the remaining analyzed cases.A priori, such a difference seems to be related to wake rollup; however, further investigation should be done in order to clarify such behavior.

Validation: reached CL, CD and CM (
Despite the fact that reliable experimental data for this particular configuration is limited, CD numerical values match excellently almost the entire AoA range, while for CL and CM there are some discrepancies (Fig. 15), probably related to the zero-thickness assumption (higher CL and CM values are obtained).More precise numerical simulations as well as more reliable experimental data are necessary in order to validate this configuration with greater confidence.
Note that, even for this near-stall condition (Fig. 16), the structures of the lateral vortices are well conserved (no mixing flow is present); thus, turbulence effects could be neglected.A future turbulent simulation should confirm this assumption.

Validation: CL, CD and CM evolution (arrow AR = 1 , 16 × 16)
A more complex flow pattern case (45 deg.swept-back flat plate) is solved through the inviscid FTVM.Comparison with straight wakes and variable volumes (at α = 25 and α = 30 deg.) is also shown.In order to maintain congruence with both previous ana- lyzed configurations, the time step was selected based on the distortion of the bounded  elements (panels with a rhomboidal shape) compared to the quadrangular ones previously analyzed.Note that the factor used is 2/3, which is directly related to the total chordwise size of the plates (1 for both quadrangular and rectangular, and 1.5 for sweptback; 1/1.5 = 2/3).
As in both previous analyzed cases, all hydrodynamic coefficients tend to converge relatively smoothly to a well-defined steady-state value as the number of iterations increases (Fig. 17).However, for the variable volume scheme at α = 30 deg. (dashed blue curve), such a solution diverges abruptly after reaching stable behavior; such an issue seems to be related to numerics (treatment of extremely low vorticity values near zero or, even, another cause); a further investigation could help to solve this issue.Note that for the α = 25 deg.case, such behavior does not appear, as for both previous rectangular configurations analyzed, since such a variable volume scheme is completely stable, even for high AoA cases ( α = 40 deg.).

Validation: reached CL, CD and CM (arrow
Leaving aside the certainty of the experimental data since they are limited, from Fig. 18, it is clear that the present method solves precisely for some ranges despite the fact that a deeper analysis must be performed by changing some input parameters since the complexity of the flow increases as the AoA increases; for instance, the same t value could not be the optimal one for the entire AoA range.Probably, a variable t (self-adjusting at each time step) could be the optimal way to obtain the most precise possible solution.
From Fig. 19, as expected, the vortex strength of the wingtips is lower than for both the quadrangular and rectangular cases previously analyzed; however, a more complex flow pattern is obtained close to the rear part of the plate, which could cause some instabilities in the solution if they are not treated properly.

Viscous (growing volumes) results
This subsection is left pending due to a lack of computational resources since higher discretization, thus lower ε is necessary in order to reduce the boundary layer thickness and capture precisely within the viscous regime (at the low AoA range).

Conclusions
The present Full Non-linear Vortex Tube-Vorton Method (FTVM) represents the unsteady flow (and fluid) solution past a shell-body along the pre-stall range, despite the non-intuitively imposed LEV (counter-clockwise rotation; attached LEV) and detached vorticity from the entire surface, an inviscid mechanism that has recently been proved theoretically in [33].
Since the present method is quite different from the classic VPM approach (and its reformulated version, the rVPM [11]), ideally, the FTVM does not require reconstructing any divergence-free vorticity field from a divergence-free velocity one since the circulation and vorticity vectors remain aligned within the scope of the isolated VFM.
As has been demonstrated in this research, maintaining a constant volume of vorticity is a practical way to solve the inviscid flow dynamics, since such results converge to the variable volume ones as both the number of discretized elements and iterations increase.This last avoids the use of an additional scheme (or convenient numerical manipulations) as 'resetting' for the volumes of vorticity to restrict their overgrowing or collapsing into a singular point (and then overshooting their vorticity values) as the simulation evolves.
Despite the fact that a precise way to account for viscous diffusion has been numerically justified through the Core Spreading Method [15], a more robust analysis is needed in order to validate the viscous version of the FTVM.A priori, better resolution simulations (higher discretization, thus, higher computational resources) are needed to precisely account for vorticity interaction near the solid walls; in the limit, such a simulation should be equivalent to a Lagrangian-based Direct Numerical Simulation if the scale resolution is fine enough.
Although the FTVM has been validated for relatively simple configurations and flow patterns, it should work the same as for more complex geometries and orientations (the steady base model has already been satisfactorily tested for sideslip conditions), including curved surfaces (by considering a discretization error due to the displaced bounded vortex ring respect to the physical panel, inherent to VLM).However, an exact solution for curved surfaces should also be possible to implement.Furthermore, such a full detached vorticity approach could be extended for general panel methods (aerodynamic bodies and internal flows).Its extension to the post-stall condition could be possible by including an intuitive LEV (clockwise; detached LEV), which should explain the instantaneous drop in lift (after the stall point) due to the plate's LE becoming an edge with a typical Kutta condition (thus, it does not contribute to the force calculation; an advanced algorithm could automatically detect the inversion of the LEV, depending on the instantaneous flow conditions).However, turbulence effects would be accounted for in order to obtain comparable results with the reference experimental data.

Verification of Eq. (4):
For the second time step (vortex squeezing):    Appendix 2 A series of simulations have been performed in order to compare the numerical behavior between two vortex strain schemes: maintaining a constant [14] and variable volumes of vorticity (described in subsection 2.1) for all the wake elements.
Results for CL, CD and CM are compared between a pure advection calculation (constant volumes) and the novel implementation for vortex squeezing/stretching (variable volumes) (Fig. 20).
Results from Fig. 21 show that even for a high AoA condition (high flow perturbation, thus large wake grid deformation), the relative error for all compared hydrodynamic coefficients is minimal (less than 0.6% for a sufficiently long simulation with relatively few discretization elements); however, for this case, the relative error does not show a trend to converge to zero but remains constant.Case: AR = 1 (10 × 10) at α = 5 • with t = 0.1 , σ 0 = ε = 0.0707 , ν = 0 , ρ = 1 , and q ∞ = 1.
As expected, increasing the number of elements (and testing for a lower AoA to reach the steady-state soon) seems to alleviate the convergence to zero issue present in the previous analyzed case.
Note that all corresponding results between constant and variable volumes are practically overlaid throughout the entire simulation.Thus, it is clear that by increasing the number of emitted vortex elements (40 instead of 12 at each time step), the relative error for all hydrodynamic coefficients tends to zero as the number of iterations increases (Fig. 22).
Even for a higher AoA condition ( α = 40 deg.) which can be considered the practical limit for this configuration through the present method, the maximum relative error remains around 0.05% for all coefficients at the end of simulation (after 100 iterations) (Figs.23 and 24).26 Vortex filaments' wake for a divergence-free grid without any near-wall treatment (the vortex elements cross the shell-body freely; scheme 1) Fig. 27 Vortex filaments' wake for a non-divergence-free grid with a near-wall treatment only for the crossing vortons (they are relocated over the plate at a prescribed distance; scheme 2) Fig. 28 Vortex filaments' wake for an enforced divergence-free grid, with a near-wall treatment for the crossing vortons and their neighboring ones (their connecting nodes are moved to enforce the interconnection of the vortex elements; scheme 3)

Fig. 2
Fig. 2 Schematic of stretching (elongation) due to pure inviscid advection of a vortex tube between two consecutive time steps, from a to b; squeezing (contraction) can be schematized from b to a

2nd order Gaussian 1 − e −ρ 3 √ 2 − ρ 2 π e − ρ 2 2 Fig. 6
Fig. 6 New vorton's volume after a pure inviscid advection step (left) and by viscous diffusion through the CSM (right) by precisely conserving both vorticity and circulation; notice that in the second case the vorton's volume is slightly larger than in the first one

Fig. 7
Fig. 7 Example of pure advection of vorticity for a leading edge's vortex element (2 × 2 discretization) between two instants: previous (gray) and current (blue)

Fig. 8
Fig. 8 Representation of a vortex filament/tube into vortons to ensure an equivalent induced velocity field

Fig. 9
Fig. 9 Conversion from a single vorton per filament to b multiple ones (3 by default), ensuring comparable and consistent results with the original (vortex filament based) version

Fig. 10
Fig. 10 Unsteady CL and CD values vs dimensionless time for different discretizations and angles of attack, compared to straight wake values (variable vortex volume scheme) swept-back ( AR = 1) 16 × 16 0.0417 0.0442 dashed line; within the straight wake assumption).The dashed curve (at α = 40 deg.) corresponds to variable volume (precise vortex stretching calculation) results.From Fig. 11, note that variable volume results (due to strict conservation of both vorticity and circulation) practically converge to the steady constant volume value, and only slight different values are present during the first iterations (along the unsteady phase).

Fig. 11 Fig. 12 Fig. 13
Fig. 11 Unsteady CL, CD and CM vs dimensionless time for an AR = 1 flat plate at different AoAs

Fig. 14
Fig. 14 Unsteady CL, CD and CM vs dimensionless time for an AR = 0.5 flat plate at different AoAs

Fig. 17
Fig. 17 Unsteady CL, CD and CM vs dimensionless time for a 45 deg.swept-back AR = 1 flat plate at different AoAs

Fig. 18 Fig. 19
Fig.18 Averaged lift, drag and pitching moment coefficients vs AoA for a 45 deg.swept-back AR = 1 flat plate compared to experimental data[32] (not available for CD)

Fig.
Fig.26 Vortex filaments' wake for a divergence-free grid without any near-wall treatment (the vortex elements cross the shell-body freely; scheme 1)

Table 2
Main input parameters for mesh independence analysis

Table 3
Main input parameters for the validation analysis