Convergence, Node Collocation and the Disc-Edge Singularity of a Vortex-Ring/Vortex-Cylinder Free-Wake Model for the Uniformly Loaded Actuator Disc

1. Introduction

Rotors, and wind-turbine rotors in particular, have been modelled by an axisymmetric actuator disc (AD) for more than 150 years [1,2]. In its classical form the model is evaluated by an integral momentum balance in which only the axial velocity component appears [3]. It is easily shown, however, that such a description is not consistent with the differential equations of inviscid flow, Euler’s equations,

$$\begin{array}{cc}\hfill \nabla ·\mathbf{u}& =0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \rho {\displaystyle \frac{D\mathbf{u}}{Dt}}& =-\nabla p+\mathbf{f},\hfill \end{array}$$

(2)

because the radial velocity component that necessarily accompanies a contracting or expanding slipstream is discarded. Among the ways of restoring consistency without solving (1)–(2) numerically, vortex-ring (free-wake) methods [4,5] are arguably the simplest: the wake boundary is represented by a free vortex sheet whose position and strength are obtained from two physical conditions, and the velocity field follows from the Biot–Savart law in closed (elliptic-integral) form.

Starting from ye’s 1990 paper [6], this approach was developed further and applied by van Kuik and Lignarolo [7] and by Bontempo and Manna [8,9]; Conway [10] generalised the underlying analytical framework to non-uniform heavy loading. Despite this long history, three questions of practical relevance for anybody implementing such a model have, to the best of the author’s knowledge, not been answered in the open literature:

1.

What do the usual convergence measures actually measure? Free-wake AD codes are commonly stopped either on a residual of the wake (sheet) equations or on the deviation of the computed power coefficient $c_P$ from its momentum-theory value. Are these two criteria equivalent, and does either of them certify a converged solution?

2.

Is the free-wake fixed-point iteration stable, and at what price? A cylindrical vortex sheet is Kelvin–Helmholtz (KH) unstable, so any explicit update of the sheet position must be damped. Which damping keeps the scheme stable without corrupting the solution near the disc edge?

3.

What does a converged free-wake solution say about the disc-edge singularity? Two mutually exclusive edge behaviours have been proposed in the literature: a self-similar (Prandtl/Kaden) spiral with $\gamma \propto s^{-1/2}$ [5,11] and a bounded sheet strength with logarithmic velocity singularities [12].

The present paper answers these questions with the help of an independent, open-source FORTRAN 90/95 [13] implementation of the vortex-ring/vortex-cylinder model for a uniformly (constantly) loaded actuator disc, together with a 1:1 Python replica of the same algorithm that allows systematic parameter studies. The main findings are:

• The sheet residual and the momentum-theory deviation are not equivalent (Section 3). The latter possesses a discretisation floor and is not even monotone during the iteration; using it as a stopping criterion can flatter the apparent accuracy by more than an order of magnitude.
• The KH saw-tooth mode is controlled by a damping factor $\lambda \propto \sqrt{z}$, which is stable for all ring numbers tested and equilibrates the edge region about two orders of magnitude faster than the $\lambda \propto z$ damping used previously (Section 3.2).
• With midpoint collocation, $c_P$ retains a small residual fluctuation band that shrinks with the discretisation but not with the number of iterations; it is a property of the unresolved edge region, not an iteration deficiency (Section 3.3). Replacing midpoint collocation by node collocation with a spacing-proportional vortex kernel (Section 2.4) removes the band entirely and yields a unique fixed point with machine-level residuals.
• The converged solutions of this model class show a bounded edge strength, with a fitted exponent $a=0.00\pm 0.04$ in $\gamma \propto s^a$ — different from both proposals cited above. We show why a single-valued, midpoint- or kernel-regularised sheet discretisation is structurally unable to represent a rolled-up spiral, and why this does not affect the computed $c_P$ (Section 5).

The code, its Python replica, the configuration files of all cases discussed here and the GNUPLOT/Python post-processing scripts are available at https://github.com/Schaffarczyk/RingsCylinderActuatorDisk/.

2. Model and Numerical Method

2.1. Model Description

Axisymmetric coordinates $(z,r)$ are used with the velocity $\mathbf{u}=(u_z,u_r)$, non-dimensionalised with the undisturbed wind speed, $u_z{\left(r\right)}_{z\to -\infty }=1$, and with the disc radius, so that the disc occupies $z=0$, $0\le r<1$. Its edge ${\mathbf{r}}_E=(0,1)$ may carry a singularity, i.e. a point at which $\mathbf{u}$ is not defined; this is the subject of Section 5.

The wake boundary is modelled by a piecewise linear vortex sheet which starts at ${\mathbf{r}}_{edge}=(0,1)$ and ends at ${\mathbf{r}}_{Cyl}=(z_{SIVC},r_{SIVC})$, where a semi-infinite vortex cylinder (SIVC) with tangential vorticity is appended. For $0<z<z_{SIVC}$, N vortex rings of strength $\gamma$ are placed with non-uniform (cosine-clustered) spacing at positions $(zs{l}_i,rs{l}_i)$,

$$\begin{array}{cc}\hfill z_i& =z_{SIVC}\left(1-cos{\phi }_i\right),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\phi }_i& ={\displaystyle \frac{\pi }{2}}·{\displaystyle \frac{i}{N}},1\le i\le N,\hfill \end{array}$$

(4)

so that the resolution is highest at the disc edge, where the solution varies most rapidly. This part of the approach is identical to that of [8,9].

The single physical parameter of the model is the thrust coefficient $c_T$, with $c_T>0$ for the propeller case (induced velocities $>0$) and $c_T<0$ for the turbine case; momentum theory is valid for $-1\le c_T<\infty$. Several exact relations follow from momentum theory and are used below both as initial conditions and as reference values. The asymptotic radial extension of the slipstream is

$$r_{z\to \infty }^{sl}\equiv r_{SIVC}=\sqrt{{\displaystyle \frac{1+\sqrt{1+c_T}}{2\sqrt{1+c_T}}}},$$

(5)

giving $r_{\infty }=0.9239$ for $c_T=1$ and $r_{\infty }=\sqrt{2}=1.4142$ for $c_T=-8/9$ (the Betz case). The value of the total (i.e. inflow included) stream function on the slipstream far downstream is

$${\psi }_{wake}={\displaystyle \frac{1}{4}}\left(1+\sqrt{1+c_T}\right),$$

(6)

i.e. ${\psi }_{wake}=0.603553$ for $c_T=1$ and ${\psi }_{wake}=1/3$ for $c_T=-8/9$. (If only the induced part of the stream function is used, the corresponding value is $c_T/\left(4\sqrt{1+c_T}\right)$; the code works with the total stream function, Equation (6).) The model parameters are summarised in Table 1.

2.2. Induced Velocities

An accurate evaluation of the induced axial and radial velocities is the key to any arbitrarily accurate — that is, successful — implementation. The computation is split into three parts:

1.

the velocities induced by the $N-1$ remaining vortex rings on the ring under consideration, $0\le z\le z_{SIVC}$;

2.

the self-induced velocity of the ring under consideration. Bontempo’s expressions [9] (Equations (6) and (7) therein) are used, which, remarkably, do not contain any core parameter:

$$\begin{array}{cc}\hfill v_{z,m\to m}& ={\gamma }_m\left\{-{\displaystyle \frac{{\beta }_{m+1}-{\beta }_m}{8\pi }}cos{\beta }_m-{\displaystyle \frac{d{s}_m}{4\pi r_m}}\left[log\left({\displaystyle \frac{8\pi r_m}{d{s}_m}}\right)-0.25\right]\right\},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill v_{r,m\to m}& ={\gamma }_m\left\{sin{\beta }_m-{\displaystyle \frac{{\beta }_{m+1}-{\beta }_m}{8\pi }}\right\},\hfill \end{array}$$

(8)

with ${\beta }_m$ the local slope of the sheet at ring m, $r_m$ its radial position and $d{s}_m$ the distance between the two neighbouring rings. For the self-induction of vortex rings in general see, e.g., [15];

3.

the velocities induced by the semi-infinite vortex cylinder at $z_{SIVC}$, for which the closed-form expressions of [4] (chapter 36) are used.

To avoid singular velocities on the sheet itself, all velocities are evaluated at the arithmetic midpoints of two neighbouring (singularity-carrying) sheet points. No artificial vortex core is introduced. This midpoint collocation is the classical choice; Section 2.4 discusses a structurally different alternative.

As an independent check against implementation errors, three different sets of equations for the induced velocities can be selected at run time: (i) Equations (1)–(7) of [8,9], (ii) Equation (D1) of [5], and (iii) the expressions of chapters 33 and 36 of [4]. All three agree to the accuracy of the elliptic-integral evaluation. For consistency of nomenclature, the formulation of Branlard [4] is used throughout. The complete and incomplete elliptic integrals (chapter C.4 of [4]) are evaluated with the open-source routine RFL of Carlson and Notis [16]; the implementation was verified against Mathematica to 10 significant digits.

The accuracy of the induced velocities was verified separately for the ring part ($0<z<z_{SIVC}$) and for the cylinder part ($z\ge z_{SIVC}$), and against the data of van Kuik [14].

2.3. Free-Wake (Slipstream) Iteration

The shape of the slipstream is determined by the induced velocities and vice versa, through two conditions [5,8,9]:

• the sheet must be force free, implemented as $\gamma =-c_T/\left(2u_s\right)$ with $u_s^2={(1+u_z)}^2+u_r^2$;
• the sheet must be aligned with the flow, implemented as ${\psi }_w\left(z\right)={\psi }_w(\infty )={\psi }_{wake}$, Equation (6).

An initial guess for the slipstream shape reduces the number of iterations. Three options are implemented: a cylinder of radius $r_{SIVC}$; the prescribed exponential shape

$$r_{sl}\left(z\right)=r_{SIVC}+\left(r_0-r_{SIVC}\right)exp(-az);$$

(9)

and a prescribed algebraic shape deduced from the axial velocity on the axis of a vortex cylinder, attributed by van Kuik [5] to Snel and Schepers,

$$u_z(r=0,z)={\displaystyle \frac{{\gamma }_{\infty }}{2}}\left(1+{\displaystyle \frac{z}{\sqrt{z^2+r_{SIVC}^2}}}\right).$$

(10)

In Equation (9), $r_{SIVC}$ is given by Equation (5) while the decay parameter a (not to be confused with the axial induction factor) has to be estimated, $a\approx 1\dots 2$. Figure 1 shows a two-parameter fit of the $c_T=1$ slipstream data of van Kuik [14], who normalises with $r_{SIVC,vK}:=1$; the fit gives $r_0=1.03498$, $r_{SIVC}=1.00006$ (as it should) and $a=1.01452$.

Several update schemes for the sheet position have been implemented, selected by the parameter iter-sc:

(1)

Equations (13) and (14) of [9];

(2)

shifting the sheet points normal to the current sheet direction;

(3)

van Kuik’s [5] Equation (D.12),

$$\partial r\left(z\right)={\displaystyle \frac{\partial \psi }{r{u}_z}},\partial \psi =\psi (\mathrm{sheet},z)-{\psi }_{wake},$$

(11)

damped with $\lambda \left(z\right)=min(z,0.1)·\mathit{under}$;

(7)

as (3), but with the damping factor $\lambda \left(z\right)=min\left(0.1,\mathit{under}·\sqrt{z}\right)$;

(8)

as (7), plus Anderson acceleration (window 4) of the resulting fixed-point iteration.

Schemes 7 and 8 are the recommended ones; among the older schemes, No. 3 was preferred. All schemes may be under-relaxed, with $\mathit{under}=0.001\dots 0.01$ for scheme 1, $<0.01$ for scheme 2, $0.1\dots 0.9$ for scheme 3 and $0.05\dots 0.1$ for schemes 7 and 8, where $\mathit{under}$ acts as the prefactor of $\sqrt{z}$.

The iteration is terminated when (i) the number of iterations exceeds maxiter, (ii) the sheet residual satisfies $\mathrm{errsl}\le \mathtt{epsSL}$, or (iii) (schemes 7 and 8) $c_P$ has stopped moving, $|\Delta c_P|<{10}^{-6}$ over the last 100 iterations (Cauchy criterion). In earlier versions of the code the iteration was stopped when $\mathrm{err}{c}_P:=|c_P-c_{P,mom}|$ fell below a prescribed tolerance. As shown in Section 3, $\mathrm{err}{c}_P$ is not a valid iteration criterion.

2.4. Node Collocation with a Spacing-Proportional Vortex Kernel

The analysis of Section 3 and Section 5 identifies a structural conflict of the midpoint-collocation scheme: stability requires a position-dependent damping, the edge region remains under-determined, and $c_P$ retains a fluctuation band. A fourth version of the code therefore changes the discretisation itself, following van Kuik and Lignarolo [7] (Appendices A and B):

1.

N vortex rings with circulations ${\Gamma }_i$ are collocated at the nodes $(z_i,r_i)$; ring 1 is pinned at the disc edge $(0,1)$ and exempted from the boundary conditions, in line with the discussion of the leading edge in [7] (Section 4.3 therein).

2.

A vortex kernel regularises all evaluations at distances $\rho <\frac{1}{2}\beta d{s}_i$ using smoothed (Marshall-type) expressions. Deviating from the fixed $\delta =0.002$ of [7], the kernel radius is taken proportional to the local spacing, ${\delta }_i=\beta d{s}_i$ with $\beta =e^{-1/2}\approx 0.607$, which matches the analytic contribution of the local sheet strip. A fixed $\delta$ biases $c_P$ by $\sim {10}^{-2}$ as soon as $ds\gg \delta$, i.e. as soon as the spacing varies appreciably along the sheet, as it does with the clustering of Equations (3)–(4).

3.

The update is single-phase with uniform damping $d=0.05$,

$${\Gamma }_i\leftarrow {\Gamma }_i+d\left({\Gamma }_i^{*}-{\Gamma }_i\right),{\Gamma }_i^{*}=-{\displaystyle \frac{c_T}{2v_i}}d{s}_i,r_i\leftarrow r_i+d{\displaystyle \frac{{\psi }_{wake}-{\psi }_i}{r_i\left(1+v_{z,i}\right)}}.$$

(12)

The kernel removes the saw-tooth mode (Section 3.2), so no z-dependent damping is needed. The normal-velocity (“fine-tuning”) phase of [7] (Appendix B.3) was tested and found to be neither stable nor necessary in this implementation: after the $\psi$-phase, $v_n\approx 3·{10}^{-2}$ at the edge and $\approx 2·{10}^{-3}$ on the sheet; it is retained only as a diagnostic (see Section 5.4).

2.5. Implementation Notes

The code is written in FORTRAN 90/95 [13] and deliberately kept simple: a minimal set of language features is used and no algorithmic “tricks” are applied, so that the source remains readable and maintainable. It consists of three files only: mem.f (global variables), VortexCode.f (the slipstream iteration loop) and SubA.f (the function library, including the elliptic integrals). Compilation uses GNU Fortran [17] with the mandatory options -fno-automatic (local variables in subroutines and functions) and -fdefault-real-8 (64-bit reals, $\approx 15$–16 decimal digits). The latter must be applied to all source files including mem.f; otherwise the global arrays of the module remain 32-bit while the rest of the code works in double precision — an error that is easily overlooked and that limits the attainable residual to $\sim {10}^{-6}$.

The power coefficient is obtained from the axial velocity at the disc following Equations (15)–(20) of [9]. The radial integration of $u_z(z=0,r)$ from $r=0$ to $r=1$ uses Simpson’s rule; since this is second-order accurate, at least ${10}^3$ points are required for an accuracy of ${10}^{-6}$, and 3000 points are used throughout. As expected for a scheme in which every ring interacts with every other ring, the CPU time scales as $\sim N^2$ (Figure 2).

All results reported below were, in addition, reproduced with a 1:1 Python replica of the code (numpy/scipy, identical Carlson elliptic integrals), which made the systematic parameter studies of Section 3 ($N=300\dots 4000$, up to $15000$ iterations) practicable.

3. Convergence Measures, Stability and the Power Coefficient

A numerical scheme deserves to be called meaningful only if it satisfies the elementary requirement

$$\forall ϵ>0\exists n\in \mathbb{N}:\mathrm{after}n\mathrm{iterations}\mathrm{the}\mathrm{result}\mathrm{is}\mathrm{accurate}\mathrm{to}ϵ.$$

(13)

Whether a given free-wake code satisfies Equation (13) depends on which quantity “the result” refers to. This section shows that the two quantities commonly monitored behave very differently in this respect.

3.1. What the Two Convergence Measures Measure

The code monitors

$$\begin{array}{cc}\hfill \mathrm{errsl}& ={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\left|\psi (z_i,r_i)-{\psi }_{wake}\right|d{s}_i,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \mathrm{err}{c}_P& =\left|{\displaystyle \frac{c_P-c_{P,mom}}{c_{P,mom}}}\right|.\hfill \end{array}$$

(15)

These measure two fundamentally different things. $\mathrm{errsl}$ is the residual of the discrete system: it can, in principle, be driven to zero by iteration alone. $\mathrm{err}{c}_P$ is the distance between the discrete solution and the continuum momentum-theory limit. Even at the exact discrete fixed point ($\mathrm{errsl}=0$), $\mathrm{err}{c}_P$ retains a discretisation floor made up of (i) a contribution $\sim 1/N^2$ from the ring discretisation (Table 2), (ii) a truncation error from the finite $z_{SIVC}$, and (iii) the contribution of the disc-edge region, which the panels cannot resolve (Section 5). Consequently $\mathrm{errsl}<{10}^{-8}$ neither implies nor can imply $\mathrm{err}{c}_P<{10}^{-8}$.

Three further mechanisms decouple the two measures:

1.

The $ds$-weighting makes $\mathrm{errsl}$ blind at the disc edge. With the cosine clustering of Equations (3)–(4), $ds\approx z_{SIVC}{\pi }^2/\left(4N^2\right)$ at the first ring — precisely where $|\Delta \psi |$ is largest (edge singularity) and where $c_P$ is most sensitive, because the weight $2rdr$ of the power integration concentrates near $r=1$. The code therefore also reports the unweighted $max|\Delta \psi |$.

2.

A damping factor $\lambda \propto z$ freezes the edge. The relaxation factor $\lambda \left(z\right)=min(z,0.1)$ of scheme 3 vanishes at the disc, so the equilibration time of the first ring is $\sim 1/z_1\approx {10}^4$ iterations for $N=600$, $z_{SIVC}=25$: $\mathrm{errsl}$ reaches its floor after $\approx 200$ iterations while $c_P$ keeps creeping for thousands of iterations. A reference run over $15000$ iterations never reached a fixed point; $c_P$ wandered by $\pm 2·{10}^{-3}$ around $c_{P,mom}$ at $\mathrm{errsl}\approx {10}^{-9}$.

3.

$\mathrm{err}{c}_P$ is not monotone. During this creep $c_P$ crosses $c_{P,mom}$. Values as small as $\mathrm{err}{c}_P\sim {10}^{-6}$ occur at these zero crossings without indicating convergence at all.

Item 3 has a direct practical consequence: accuracy tables obtained by stopping the iteration on $\mathrm{err}{c}_P$ — such as Table 2, which is representative of what is found in the literature — must be read with care, since they may report the position of a zero crossing rather than a converged result.

Two implementation errors were found and corrected in the course of this analysis and are reported here because they are easy to reproduce in any comparable code: the segment length $ds$ entering Equation (14) was taken from a stale variable, which made $\mathrm{errsl}$ about 50 times too large, and the under-relaxation factor was applied twice, giving an effective relaxation $min(z,0.1)·{\mathit{under}}^2$.

3.2. Stability: The Saw-Tooth (Kelvin–Helmholtz) Mode

A cylindrical vortex sheet is KH unstable; the discrete manifestation is a saw-tooth mode on the grid scale. Systematic tests with the Python replica ($N=600$) gave the following picture:

• undamped Newton steps, Equation (11) with $\lambda =1$: the saw-tooth amplitude grows from ${10}^{-2}$ to $O\left(1\right)$ and $c_P$ diverges;
• constant damping: marginally stable at $\lambda \approx 0.003$ for $N=600$ and unstable for larger $\lambda$; the threshold decreases with N;
• $\lambda \propto z$ (scheme 3): stable, but the edge is frozen (see Section 3.1);
• $\lambda =min(0.1,\mathit{under}\sqrt{z})$ (scheme 7): stable for $N=600$ and $N=1200$, with an edge equilibration about 100 times faster than scheme 3.

Purely local alternatives — per-station Newton iterations with a geometry-consistent finite-difference Jacobian, Gauss–Seidel sweeps, and naive streamline tracing — were also tested. They either diverge or converge to a shifted fixed point: the ring-to-ring coupling near the edge is too strong for any purely local update. This is a genuine property of the problem, not of a particular implementation, and explains why free-wake AD codes need global (or strongly damped) updates.

3.3. The $c_P$ Band and Its Reduction by Richardson Extrapolation

With scheme 7, $\mathrm{err}{c}_P$ drops to $\approx {10}^{-4}$ within $\approx 350$ iterations ($N=600$, $z_{SIVC}=25$) and then oscillates slowly around $c_{P,mom}$; with the Anderson-accelerated scheme 8, $\mathrm{errsl}$ reaches its floor three to four times faster still. The remaining oscillation band of $c_P$, $\pm (3\dots 5)·{10}^{-4}$ for $N=600$, $z_{SIVC}=25$, is not an iteration deficiency: at $\mathrm{errsl}\approx {10}^{-8}$ the discrete equations leave the unresolved edge region — and with it $c_P$ — undetermined within this band. Different (all stable) update schemes converge to slightly different pseudo-solutions inside the same band. The band shrinks with increasing N and $z_{SIVC}$, but not with more iterations.

The code therefore reports, in addition to the last iterate, an exponentially weighted moving average of $c_P$ (time constant $\approx 50$ iterations) and the fluctuation band $\pm (max-min)/2$ over the last 100 iterations. We recommend quoting $c_P$ as “average ± band” and reducing the band by an N-study with Richardson extrapolation.

Table 3 and Figure 3 show such a study for the Betz case ($c_T=-8/9$, exact momentum-theory value $c_{P,mom}=-16/27=-0.592593$) computed with scheme 7 at $z_{SIVC}=50$ and 1100–1200 iterations per run. All extrapolation variants land at $c_{P,\infty }=-0.59250\dots -0.59253$, i.e. within $(0.6\dots 1.1)·{10}^{-4}$ of $-16/27$. The observed convergence order is $p\approx 1.2$ — edge-dominated, and clearly not the smooth $1/N^2$ of the ring discretisation away from the edge. The recommended citation for this protocol is $c_P=-0.59251\pm 0.00007$.

A complementary study at fixed $N=2000$ with $z_{SIVC}=50$, 100 and 200 gave the perhaps counter-intuitive result that a larger $z_{SIVC}$ degrades the result: the bands grow from $6.8·{10}^{-5}$ to $2.0·{10}^{-4}$ and $3.9·{10}^{-4}$, and the means move away from $-16/27$. The reason is that at fixed N all cell sizes scale $\propto z_{SIVC}$ — the first (edge) cell even $\propto z_{SIVC}/N^2$ — so the loss of resolution outweighs the smaller truncation error of the cylinder, which is already subdominant at $z_{SIVC}=50$. The practical rule is therefore: use $z_{SIVC}=50$ and invest the available effort in N; a double extrapolation in $(N,z_{SIVC})$ is unnecessary.

3.4. Node Collocation: A Unique Fixed Point

The node-collocation discretisation of Section 2.4 changes this picture qualitatively. Table 4 lists the converged results for the Betz case: the residual $max|\Delta \psi |$ reaches machine level ($\sim {10}^{-12}$), there is no $c_P$ band and no drift, i.e. the iteration converges to a unique fixed point. The error in $c_P$ saturates between $N=1600$ and $N=3200$: the fixed points converge (Aitken limit) to $c_{P,lim}=-0.592545$, leaving an N-independent model bias of $+4.7·{10}^{-5}$ with respect to momentum theory, while the pre-saturation discretisation error decays with an observed order $p\approx 1$ (Figure 4). Already at $N=1600$ the error is $+4.9·{10}^{-5}$, i.e. better than the midpoint scheme at $N=4000$.

Separating the origin of the remaining bias — truncation of the vortex tube at $L_{tube}=30R_{wake}$ versus the kernel model at the tube junction and at the edge — requires runs with $N\propto L_{tube}$; at fixed N, enlarging $L_{tube}$ degrades the result because all cells coarsen, the same resolution/truncation confounding as in the $z_{SIVC}$-study above. In contrast to midpoint collocation, the remaining $c_P$ error is a smooth, pure discretisation error without band or drift, so that Richardson extrapolation is clean (Figure 4). The edge singularity itself remains kernel-smoothed; resolving the conjectured $1/\sqrt{s}$ zone (Section 5) would still require an analytic edge cell or a multivalued sheet parametrisation.

4. Results for the Two Reference Cases

4.1. Propeller, $c_T=1$

The propeller case with $c_T=1$ is only moderately non-linear: the slipstream contracts to $r_{\infty }=0.9239$, i.e. to about 0.93 of the disc radius, and $c_P=1.207$, to be compared with the linear approximation $c_{P,lin}=1+c_T/4=1.25$ of [18] (Equation (9.107) therein). Figure 5, Figure 6 and Figure 7 compare the computed axial velocity and vortex strength along the slipstream, and the streamline pattern, with the data of van Kuik [14]. The agreement is good over the whole range except for the immediate vicinity of the disc edge, $z\lesssim {10}^{-3}$, where the two codes deviate — exactly the region discussed in Section 5. The computed edge streamline coincides with ${\psi }_{wake}=0.60355339$ from Equation (6), which is an independent check of the converged solution.

4.2. Betz Turbine, $c_T=-8/9$

For the Betz case the slipstream expands to $r_{\infty }=\sqrt{2}$ and ${\psi }_{wake}=1/3$. Figure 8, Figure 9, Figure 10 and Figure 11 show the axial and radial induced velocities, the vortex strength and the streamline pattern. Both $v_z$ and $v_r$ tend to finite values as $z\to 0^{+}$ in these regularised solutions — a point taken up in Section 5. The converged power coefficient is $c_P=-0.59251\pm 0.00007$ (Section 3.3) with midpoint collocation and $c_P=-0.5925454$ with node collocation at $N=3200$ (Section 3.4), against the momentum-theory value $-16/27=-0.592593$.

5. The Disc-Edge Singularity

Along a force-free vortex sheet, $\gamma u_s=-c_T/2$ is constant. Van Kuik [5,11] proposes that at the disc edge the wake boundary rolls up into a self-similar (Prandtl/Kaden-type) spiral whose unrolled strength diverges as

$$\gamma \left(s\right)\propto s^{-1/2},s\to 0,$$

(16)

with s the arc length measured from the edge. Is Equation (16) an established result or a hypothesis? A review of the literature together with the numerical evidence of the present model gives the following picture.

5.1. What Is Established

(a) The linear (frozen-wake) model is singular at the rim. For a semi-infinite vortex cylinder of constant strength, the induced velocity field is singular at the rim: $v_z$ diverges logarithmically and $\partial v_r/\partial s\to \infty$. This is analytic — it is visible in the very Branlard expressions [4] used in this code — and Chattot [12] derives from the Biot–Savart law that the linearised slip line behaves as $r-R\propto -xlnx$ near the edge.

(b) The non-linear force-free problem is singular at the leading edge of the sheet. All potential-flow computations that let $\gamma$ adjust to the force-free condition find the sheet strength growing (resolution-limited) towards the leading edge. Van Kuik [19] states that “in all flow cases the distributions have a singularity at the leading edge”, with a load-dependent strength; see also [7]. The qualitative statement “the disc edge is singular” is therefore established.

5.2. The $1/\sqrt{s}$ Spiral: Strongly Supported, Not Proven

The analytical basis of Equation (16) is the local non-linear analysis of Schmidt and Sparenberg [20] for a large constant normal load. Van Kuik revisited the edge in [11] and constructed the flow induced by, and a steady solution for, Prandtl’s self-similar vortex-sheet spirals [21,22]. Equation (16) is exactly the Kaden roll-up law, which is rigorously established for the unsteady two-dimensional similarity problem (see, e.g., [15]). Independent support comes from Spalart [23], who concludes — explicitly labelling it a conjecture — that “the cross-section of the vortex sheet probably begins with a ${45}^{\circ }$ spiral”, supported by a viscous numerical solution. In contrast, Chattot [12] obtains, from a linearised analysis with a prescribed constant sheet strength, logarithmic velocity singularities and a bounded $\gamma$; the disagreement is traceable to that assumption, since it is precisely the force-free coupling $\gamma =-c_T/\left(2u_s\right)$ that drives $\gamma \to \infty$ where $u_s\to 0$. No free-wake computation has resolved the spiral, and no experiment resolves the exponent. Thus the existence of the singularity is established, whereas the specific form (16) is a well-founded, independently supported, but unproven hypothesis.

5.3. Numerical Evidence from the Present Model

The converged solutions of the present model show a bounded edge strength. A log–log fit $\gamma \propto s^a$ over the first rings of converged double-precision solutions (scheme 7) gives the exponents of Table 5: $a=0.00\pm 0.04$ over two decades in s, with ${\gamma }_{edge}\approx 0.55<{\gamma }_{\infty }=2/3$ and, at matched s, N-independent values. This is consistent with the logarithmic decrease seen in Figure 10.

This result does not contradict the spiral picture. A discretisation of the present class is structurally unable to represent it, for three independent reasons:

1.

the single-valued parametrisation $r_{sl}\left(z\right)$ of the sheet excludes any multivalued (rolled-up) geometry;

2.

the self-induction term, Equations (7)–(8) (or the vortex kernel of Section 2.4), regularises the Biot–Savart kernel at the panel scale $ds$;

3.

the force-free update evaluated with this regularised velocity keeps $\gamma$ bounded, because $u_s$ cannot reach zero.

The code therefore converges to the regularised solution of the no-roll-up problem, which is why the present findings differ from van Kuik’s $s^{-1/2}$ and from Chattot’s constant-$\gamma$ picture. Practically, the singularity is integrable ($\int \gamma ds<\infty$), so $c_P$ converges unaffected; its footprint in the computation is the saturation of $max|\Delta \psi |$ at the first ring(s) and the $c_P$ band of Section 3.3. Any attempt to settle the exponent numerically will require either an analytic edge element or a multivalued sheet parametrisation — neither of which is contained in any of the free-wake AD codes known to the author.

5.4. A Note on the “Fine-Tuning” of van Kuik and Lignarolo

The convergence scheme of [7] (Appendix B.3) switches, after convergence of the $\psi$-condition with damping 0.05, to the condition $v_n=0$, applied as $\Delta r_i=-v_{n,i}d$ with $d=0.0025$ (“fine tuning”). Transferring this second phase to a midpoint-collocation code fails structurally: after $\psi$-convergence, $v_n$ at the first $\approx 5$ rings is $O\left(1\right)$ — 0.9 at the first ring — which is a discretisation artefact of midpoint collocation and not a usable residual; stepping on it destroys the solution for either sign of the update. In the node-collocation discretisation of Section 2.4 the artefact disappears ($v_n\approx 3·{10}^{-2}$ at the edge), but the $v_n$-phase is then unnecessary, because the $\psi$-iteration already converges to machine accuracy. It is worth noting that even in [7] the fine-tuning does not resolve the edge: a 40% violation of $\Delta p_{sheet}=0$ is reported at the first ring, and the authors conclude that the boundary conditions hold for a vortex sheet but not for its leading edge.

6. Discussion

Six model and discretisation parameters were found to influence the accuracy of the computed solution: (i) the number of rings N; (ii) the axial location $z_{SIVC}$ of the semi-infinite cylinder; (iii) the number of integration points at the disc; (iv) the model used for the self-induction; (v) the initial and converged shape of the slipstream; and (vi) the number and distribution of vortex rings close to the disc edge. This confirms van Kuik’s finding [5] that these parameters are interrelated rather than independent. Nevertheless, increasing N always increases the accuracy, whereas increasing $z_{SIVC}$ at fixed N does not (Section 3.3).

The results suggest the following practical recommendations for free-wake actuator-disc computations:

1.

Do not stop the iteration on the momentum-theory deviation. Use the sheet residual, monitor the unweighted $max|\Delta \psi |$ in addition to the $ds$-weighted mean, and treat $\mathrm{err}{c}_P$ as a verification quantity, never as a convergence criterion.

2.

Quote $c_P$ as “average ± band” when using midpoint collocation, and reduce the band by an N-study with Richardson extrapolation rather than by more iterations.

3.

Use $z_{SIVC}\approx 50$ and invest the remaining computational effort in N.

4.

Prefer node collocation with a spacing-proportional kernel if a unique fixed point is required, e.g. for gradient-based optimisation or for embedding the model in a larger solver. A kernel radius proportional to the local spacing is essential on clustered grids.

5.

Do not expect the edge exponent from this class of model. A bounded $\gamma$ at the edge is a property of the regularisation, not a physical result.

Returning to Equation (13): for the sheet residual, the improved iteration demonstrably satisfies the requirement, and node collocation does so with machine-level residuals. For $\mathrm{err}{c}_P$ the requirement can hold only in the combined limit (iterations $\to \infty$, $N\to \infty$, $z_{SIVC}\to \infty$), because at fixed discretisation $\mathrm{err}{c}_P$ has a floor that no iteration can undercut. This distinction is elementary, but it is not usually made explicit in the actuator-disc literature, where agreement with momentum theory to a given number of digits is often presented as evidence of numerical convergence.

The limitations of the present work should be stated clearly. The model is restricted to a uniformly loaded, axisymmetric, steady, inviscid disc without swirl; the generalisation to non-uniform loading exists analytically [10] but has not been implemented here. The remaining model bias of the node-collocation scheme ($4.7·{10}^{-5}$ in $c_P$) has not been decomposed into its tube-truncation and kernel contributions, which would require runs with $N\propto L_{tube}$. Finally, neither discretisation resolves the edge region itself, so the results of Section 5 bound what such models can say about the singularity rather than settling the question.

7. Conclusions

An independent, open-source FORTRAN 90/95 implementation of the vortex-ring/vortex-cylinder model for a uniformly loaded actuator disc, together with a 1:1 Python replica, was used to examine the numerical behaviour of this model class for the propeller case $c_T=1$ and the Betz case $c_T=-8/9$. The conclusions are:

1.

The sheet residual $\mathrm{errsl}$ and the momentum-theory deviation $\mathrm{err}{c}_P$ measure different things. $\mathrm{errsl}$ can be driven below ${10}^{-8}$ while $\mathrm{err}{c}_P$ stagnates at a discretisation floor; $\mathrm{err}{c}_P$ is not monotone and is therefore unsuitable as a stopping criterion.

2.

A damping factor $\propto \sqrt{z}$ (optionally with Anderson acceleration) stabilises the Kelvin–Helmholtz saw-tooth mode and reaches the residual floor within a few hundred iterations, where the previously used damping $\propto z$ required ${10}^4$ iterations or never settled.

3.

With midpoint collocation, $c_P$ retains a fluctuation band ($\pm 4·{10}^{-4}$ at $N=600$, $z_{SIVC}=25$) that reflects the under-determined edge region. An N-study with Richardson extrapolation gives $c_P=-0.59251\pm 0.00007$ for the Betz case, i.e. within ${10}^{-4}$ of $-16/27$.

4.

Node collocation with a spacing-proportional vortex kernel and a pinned first ring converges with uniform damping to a unique fixed point (residuals $\sim {10}^{-12}$, no band, no drift) and reaches $c_P-(-16/27)=+4.7·{10}^{-5}$, a residual model bias that is independent of N.

5.

The converged solutions exhibit a bounded edge strength with fitted exponent $a=0.00\pm 0.04$, differing from both the $\gamma \propto s^{-1/2}$ spiral and the constant-$\gamma$ picture. This is the expected behaviour of a single-valued, regularised sheet discretisation and does not decide between the two proposals; the existence of the edge singularity is established, its specific form remains an open question.

Future work should address the decomposition of the remaining model bias, the extension to non-uniformly loaded discs, and — for the edge question — an analytic or multivalued edge element capable of representing a rolled-up sheet.