This study extends the analysis by Wood et al. (2021) of the velocities induced at the blades by the helical trailing vortices for lifting line (LL) analysis of horizontal-axis turbines and propellers. Wood et al. (2021) considered only identical, straight, equispaced blades that lie in the radial plane of rotation. The vortices were assumed to have constant pitch and radius as this is the only case with closed-form expressions for the induced velocities. The restriction of constant pitch is not significant as Wood and Hammam (2022) showed that the pitch is constant throughout the wake of optimal rotors as it is the ratio of torque to thrust. The restriction of constant radius is possibly more significant, but at least, the constant radius analysis captures the important features of circumferential periodicity for identical, equispaced blades. Some further comments on the effects of wake expansion and radial velocity are made in Subsection 2.1. Wood et al. (2021) surveyed the approximate methods to compute the velocities along a straight LL based on the exact Kawada–Hardin (KH) equations described in Section 2.1. There are, however, situations that involve the induced velocities throughout the rotor plane rather than just along straight LLs. An example is a swept rotor, with the LL curved in the plane of rotation, which has been the subject of a significant amount of recent research, e.g., Fritz et al. (2022), Li et al. (2022b), and Gemaque et al. (2022). The first two studied swept wind turbine blades using different approximations for the induced velocities. Fritz et al. (2022) used a restricted version of the analysis developed here, whereas Li et al. (2022b) used a vortex-cylinder model of the wake. It is suggested that the present analysis is simpler and more general. Gemaque et al. (2022) investigated blade sweep for a hydrokinetic turbine without explicit incorporation of the spatial dependence of the induced velocities. In other words, they used Prandtl’s well-known tip loss factor in their blade element/momentum theory (BEMT) calculations. Bergmann et al. (2021) analyzed propeller performance also using BEMT and Prandtl’s tip loss factor. A possibly more important example is a coned rotor (Li et al., 2022a). Most modern large blades have pre-bend; the tips are upwind of the hub at a low wind turbine. As the wind speed and aerodynamic loads increase, the blades deflect downwind out of the radial plane of rotation due to their inherent flexibility. No study of coned rotors has considered the effect on the induced velocities of any axial displacement between the start of a trailing vortex and the control point where the induced velocity is required.

Figure 1 shows a doubly infinite helical vortex whose velocity field is given by the KH equations. When the angle θ in the figure is zero, the velocity induced along the LL by the singly infinite blue vortex shed at radius t, also along the LL, is one-half the doubly infinite value and so follows immediately from the KH equations. The present analysis accounts for a non-zero θ between the start of each trailing helical vortex and the control point, often the center of a blade element. Note that the terms “trailing” and “singly infinite” are synonymous. It will be shown that non-zero θ requires more than the KH equations. Initially, the blades are assumed to remain in the plane of rotation, and the extension to coned blades will be outlined subsequently as a small extension of the analysis.

The next section describes the KH equations and their two-term approximations for the velocities along the LLs and the new extension to the whole rotor plane. Then, the third-order approximations are derived and compared to those due to Okulov (2004) and Okulov and Sørensen (2020), which they called “remainders.” The testing of the accuracy of the approximations is described in Section 3; the remainders are generally more accurate than the third-order corrections. As mentioned above, a major complication is that the KH equations apply to a doubly infinite helical vortex, whereas the velocities at x = 0 are needed for the trailing blue vortex in Figure 1. When θ = 0, the singly infinite induced velocities are one-half the doubly infinite ones; it is shown at the start of Section 4 that the doubly infinite results for non-zero θ give the sum of the induced velocities at ±θ for a trailing vortex, and a further relation between these velocities is needed. Section 4 continues by using the Biot–Savart law to derive an approximate further relation for the difference in the velocities at ±θ. Section 5 describes numerical tests of the singly infinite equations. Section 6 discusses the results, shows how the analysis can be extended to coned rotors, and contains the conclusions. Appendix A describes the derivation of the sums needed for the θ − dependent approximations, and Appendix B shows the three terms in the approximation integration to zero over the region of planar symmetry as a check on the consistency of the derivation. Appendix C presents the analytic remainders added to the numerical integral of the Biot–Savart law for the singly infinite helices.

Wood et al. (2021) documented the range of pitch values for a number of simulations and experiments on wind and hydrokinetic turbines and propellers, which led them to choose representative values of p = 0.1, 1 and N = 4 for their numerical tests to determine the accuracy of the approximate equations for θ = 0. The accuracy depends on N, but the dependence on p is stronger; decreasing p, which corresponds to increasing the tip speed ratio for turbines, increases the accuracy. Because of the close similarity between S ₁ and S ₃, only the former is considered here for N = 3. S ₁ was calculated directly from Eq. 5 by summing to a maximum value of m = m _max , chosen to give an accuracy of 10^–8 using Eq. 18 of Wood et al. (2021), which assumes θ = 0. Since the Kapteyn series will converge more rapidly for θ ≠ 0, using an upper limit of m _max ∀ θ should be safe. S _1,A is the A − term in the approximation, that is, the first term in (21), and S _1,W includes all three terms.

The results shown in Figure 2 correspond to a tip speed ratio of approximately 7, and t = 1 and r = 0.99 represent the tip vortex and the center of the blade element at the tip, respectively, when the number of elements per blade is around 50. Since S ₁ and S ₃ are even in θ, 200 equispaced points over the positive region of symmetry, θ = [0, π], were computed, but only every fifth point is plotted in the main figure. The inset amplifies the region near θ = 0 for all points in that region. For these conditions, m _max = 59, indicating the considerable computational burden of summing the exact equations. Wrench’s approximation is very accurate; the r.m.s difference between S ₁ and S _1,W over the 200 points calculated for 0 ≤ θ ≤ π/N was almost within the tolerance of calculating S ₁. Using the A − term only gives reasonable accuracy at p = 0.1. For θ = 0, S ₁ = −0.43172, S _1,A = −0.42838, and S _1,W = −0.43172; so using the first term only in the approximation gives a modest error of less than 1%.

It is inferred from these results that the C − term or Okulov’s alternative is not significant at this p. Their behavior is shown in Figure 3. Okulov’s third term, Eq. 24, is more accurate than C _W from (21), but both could be used to determine the accuracy of using only two terms, i.e., Eqs 16, 17. Increasing p to p = 1 increases the time required for the evaluation of the exact solution. For the same N, t, and r, m _max becomes 519, and an accurate summation was not possible. r was decreased to 0.98, for which m _max = 250. The results are shown in Figures 4, 5. At θ = 0, S ₁, S _1,W , and S _1,A were 1.7152, 1.7154, and 1.6933, respectively, all multiplied by 10⁴. The accuracy of the approximations is slightly worse than that at p = 0.1, and Figure 5 shows that the two versions of the C − term have become nearly equal. In particular, it is noted that the A − term is a significantly poorer approximation to the exact solution at its higher p.

Other combinations of N, p, t, and r were investigated but not presented as the results shown are typical. The extension of Wrench’s approximations for θ = 0 to any θ was found to be accurate for a representative range of the parameters, and the accuracy can be improved and/or assessed by the C − term or Okulov’s alternative.

All the equations derived and tested so far apply to a doubly infinite helical vortex. They give twice the velocities induced by a singly infinite (trailing) vortex only when θ = 0. For non-zero θ, it will be shown below that the doubly infinite equations give u(r, θ) + u(r, − θ) and w(r, θ) + w(r, − θ) for a singly infinite vortex, and a further relation between the velocities is needed. Approximate relations for u(r, θ) − u(r, − θ) and w(r, θ) − w(r, − θ) can be derived from the Biot–Savart (BS) law for a line vortex.

Since the BS law provides a different method for determining the induced velocities, its relation to the KH equations, (1)–(4), is now briefly considered. The KH equations were derived from the velocity potential of a flow containing a doubly infinite helical vortex. The relationship with the BS law was analyzed by Morgan and Wrench (1965) and, more accessibly, by Boersma and Yakubovich (1998), whose Eq. 1 establishes the BS law as the “integral representation” of the series involving S ₁ and S ₃; see also Eq. (2.6) of Boersma and Wood (1999).

A straightforward application of the BS law gives u r , θ , w r , θ = Γ 4 π I u , B S , I w , B S = Γ 4 π ∫ 0 ∞ i u r , θ , i w r , θ d 3 d β , (26) where Γ is the vortex strength and β is the angle which is positive in the direction of positive θ. I _u and I _w , called the “influence functions,” depend only on the geometry of the vortex. The limits indicate a singly infinite vortex so that β = [0, ∞] for both θ and −θ. The integrands contain i u r , θ = t 2 − r t ⁡ cos β − θ and i w r , θ = p r − t ⁡ cos β − θ − t β ⁡ sin β − θ , (27)

and the distance d is given by d 2 r , θ = r 2 + t 2 − 2 r t ⁡ cos β − θ + p 2 β 2 . (28)

The relationship between the doubly and singly infinite velocities can be seen from ∫ − ∞ ∞ i u r , θ d 3 r , θ d β = ∫ 0 ∞ i u r , θ d 3 r , θ d β + ∫ − ∞ 0 i u r , θ d 3 r , θ d β = ∫ 0 ∞ i u r , θ d 3 r , θ d β + ∫ 0 ∞ i u r , − θ d 3 r , − θ d β = I u , B S r , θ + I u , B S r , − θ . (29)

It follows from (26) that ∫ p i u r , θ + r i w r , θ d 3 r , θ d β = p β t 2 + r 2 − 2 r t ⁡ cos θ − β + p 2 β 2 . (30) By this extended version of helical symmetry, w(r, θ) can be obtained immediately from u(r, θ), so only the latter will be analyzed. Despite (30), the integrals for a helical vortex do not have closed-form solutions. Thus, the following approximation is made for the difference in the influence functions: Δ I u = I u r , θ − I u r , − θ ≈ Δ I u , 1 = ∫ 0 2 θ i u r , θ r 2 + t 2 − 2 r t ⁡ cos β − θ + z 2 3 / 2 d β , (31) where z = pθ. In words, the difference in the velocities is approximately that induced by a sector of a vortex ring of radius t located a distance pθ behind the control point. The approximation is sketched and further explained in Figure 6. In terms of β, measured from the intersection of the blue vortex and the green plane, the sector starts at (t, 0) and ends at (t, 2θ). The integral in (31) has an analytic form that must be modified to give the correct asymptotic behavior as r ↓ 0 for any t. In that limit, the integrand in (26) becomes independent of θ, and so ΔI _u → 0. The result obtained using Mathematica is Δ I u , 1 = 2 r − t 2 + z 2 r 2 − t 2 + z 2 r + t 2 + z 2 E θ 2 , − 4 r t r − t 2 + z 2 − F θ 2 , − 4 r t r − t 2 + z 2 + 4 r t r 2 − t 2 + z 2 sin θ r 4 + 2 r 2 z 2 − t 2 + t 2 + z 2 2 r 2 + t 2 + z 2 − 2 r t ⁡ cos θ + 2 t 2 θ t 2 + z 2 3 / 2 , (32) where E(. , .) and F(. , .) are the incomplete elliptic integrals described in chapter 19 of DLMF (2023). The extra subscript was added to ΔI _u to denote the use of the full form of the integrals. Two more approximations with different additional subscripts are described below.

As for the doubly infinite results in Section 3, the tests on the accuracy of the approximations are restricted to r < t. A single vortex (N = 1) only will be considered, but some comments on the important cases of N > 1 are made below. The sum of the influence functions is Σ I u = I u r , θ + I u r , − θ ≈ Σ I u , W = t p 2 S 1 , W (33) when the C − terms are ignored on the grounds that they are likely to be comparable to, or smaller than, the errors in using (32) for most turbine applications. The extra subscript in I _u,W indicates that it was determined from Eq. 16. For comparison, the “exact” values of I _u,BS (θ) and I _u,BS (−θ) were determined by the numerical evaluation of the Biot–Savart integral for I _u in (26) using MATLAB's integral function with the default error settings, over β = [0, β _max ], from which ΔI _u,BS and ΣI _u,BS follow. To the latter was added an analytic remainder approximating the integrals over [β _max , ∞], as described in Appendix C. β _max = 1000π was used for all results presented here. By varying β _max , this value was found to result in 6-figure accuracy for the results presented in the next section.

ΣI _u,BS provides a check on the accuracy of ΣI _u,W from (33). ΔI _u,BS will be used to show the errors in ΔI _u,1 from Eq. 32, and further approximations using the simplification of the incomplete elliptic integrals described below. It is emphasized that the KH equations were found to give nearly identical results to those from the BS integration, but they are not considered further as they do not lead to an equation for ΔI _u .

Eq 33 is valid for all N, and it is unlikely that further simplification is possible. Eq. 32, however, has to be summed over N blades spaced 2π/N apart in a BEMT code, which raises three important considerations. First, the incomplete elliptic integrals are only partially periodic, as demonstrated by the small-modulus approximations of formulas (62:9:5) and (62:9:6) of Oldham et al. (2009): F θ 2 , m ∼ θ 2 + m 4 θ − 1 2 sin θ a n d E θ 2 , m ∼ θ 2 − m 4 θ − 1 2 sin θ . (34) Thus, summation over N blades will not allow simple Kawada cancellation. Second, the integrals are computationally expensive to evaluate, so the accuracy of using (34) was tested. Surprisingly, the most accurate overall results came from ignoring the terms in m, leaving only θ/2 for both F and E, which is very easy for summing. Using (34), without any further modification to (32), leads to the second approximation to ΔI _u : Δ I u , 2 = 2 t θ t t 2 + z 2 3 / 2 − r + t r + t 2 + z 2 r − t 2 + z 2 + 4 r t r 2 − t 2 + z 2 sin θ r 4 + 2 r 2 z 2 − t 2 + t 2 + z 2 2 r 2 + t 2 + z 2 − 2 r t ⁡ cos θ , (35) which will be tested for accuracy along with ΔI _u,1 from (32).

The third consideration is that the trigonometric term will simplify by Kawada cancellation when summed over N blades, but the presence of cos(θ) in the denominator makes the result more complicated than that for the KH equations. Summing over N blades is equivalent to using the trapezoidal rule to approximate the integral of the periodic function: ∫ θ θ + 2 π sin β r 2 + t 2 + z 2 − 2 r t ⁡ cos β d β = 1 r 2 + t 2 + z 2 ∫ θ θ + 2 π sin β 1 − a ⁡ cos β d β = 0 , (36) where a = 2rt/(r ² + t ² + z ²) satisfies 0 < a < 1 for a BEMT code. Note that the integral is independent of θ, but generally, the trapezoidal approximation is not, even though Kawada cancellation always applies. In the limit of a = 0, however, the sum over N is zero for any θ. Thus, the sum of the trigonometric term, which is common to ΔI _u,1 and ΔI _u,2, tends to zero as r ↓ 0 for N > 1. When a = 1, the integrand reduces to 2 cos ( β / 2 ) for β < 2π and − 2 cos ( β / 2 ) for β > 2π, again giving a closed-form sum. The general form of the sum is complicated and does not converge monotonically to zero at large N, probably because the complex variable form of the second integrand in (36) is not analytic when a = 1 (Trefethen and Weideman, 2014). When a < 1, the integrand is analytic, and the sums for large N are consistent with the exponential convergence of the trapezoidal rule for analytic functions that are periodic (Trefethen and Weideman, 2014). Nevertheless, the sum for small N when a = 1 appears to allow a simple approximation for a < 1 and the small values of N in a BEMT code: ∑ i = 0 N − 1 sin θ + 2 π i N 1 − cos θ + 2 π i N = 2 ∑ i = 0 N − 1 cos θ 2 + π i N = 2 cos θ 2 − cot π 2 N sin θ 2 , (37) provided Nθ ≤ 2π. This suggests approximating the sum for any a as ∑ i = 0 N − 1 sin θ + 2 π i N 1 − a ⁡ cos θ + 2 π i N ≈ sin θ 1 − a ⁡ cos θ − 2 cot π 2 N sin θ 2 , (38) which is exact for N = 1 as cot(π/2) = 0 and is accurate for other values of N < 10 approximately over a wide range of θ. Eq 37 shows that the sum becomes zero when N = π/θ, which is large for small θ. Thus, the sum of the non-analytic function when a = 1 is algebraic rather than N times the exponential convergence of the trapezoidal rule when a < 1 and N is sufficiently large.

The convergence of the sum in (38) for a < 1 inferred from that of the trapezoidal rule shows that there are no sweep effects on the induced velocities for an actuator disc that is approached as N → ∞, that is, much larger values than considered in the previous paragraph. Periodicity also requires u(π) = u(−π), that is, u(θ) → u(−θ) as θ → π. This requirement is not satisfied by either (32) or (35), but this error is unlikely to be significant for moderate sweep and coning. As shown at the end of the next section, the error can be easily corrected.

Figure 7 plots the influence functions over θ = [0, π] for the conditions of Figure 2, except that now N = 1, and Figure 8 shows the various ΔI _u and ΣI _u . As expected from the previous results, Wrench’s approximation for ΣI _u is very accurate at all θ. Both approximations for I _u (θ) are accurate at small θ when I _u (θ) is large in magnitude. Small θ is likely to occur frequently when the control point, r, is close to the vortex at t. There is, however, a very large difference between I _u (θ) and I _u (−θ), which remains close to zero for all θ. Only the exact influence functions asymptote to equality as θ → π, but the effect of the error in the others is likely to be small. The error is due almost entirely to the error at large θ in ΔI _u,1 and ΔI _u,2.

In Figure 9 for I _u and Figure 10 for ΔI _u and ΣI _u , r has been reduced to 0.1, which places it in the hub region where θ for t = 1 is likely to be larger than for r = 0.99 in a BEMT implementation. Now, I _u (θ) is nearly equal and opposite to I _u (−θ), and the typical magnitude has decreased substantially. Both approximations for I _u differ from the exact solution at θ = π, but the error may not be significant.

It is possible that the errors in ΔI _u are related to placing the vortex sector at the arbitrary location of pθ behind the control point. Changes to the position of the sector were investigated, but no overall improvement in accuracy was found, so all the present results were obtained using z = pθ.

A simple way to restore the periodicity to I _u is to replace θ with sin(θ) in Eq. 35. This gives the third approximation: Δ I u , 3 = 2 t ⁡ sin θ t t 2 + z 2 3 / 2 − r + t r + t 2 + z 2 r − t 2 + z 2 + 2 r r 2 − t 2 + z 2 r 4 + 2 r 2 z 2 − t 2 + t 2 + z 2 2 r 2 + t 2 + z 2 − 2 r t ⁡ cos θ . (39) This approximation is not shown in Figures 7–9 to ensure their clarity. Figures 11, 12 plot the accuracy of I _u,i for i = 1, 2, 3 for t = 1.0, p = 0.1, and r = 0.99 and 0.1, respectively. The “error” is I _u,BS − I _u,i , where ΣI _u,W is used for all i. Overall, the third approximation is better than the second and is sometimes better than the first. As expected from ΣI _u,W being substantially more accurate than ΔI _u,i for any i, the errors in I _u,i (θ) and I _u,i (−θ) are equal in magnitude and opposite in sign.

The main aim of this study was to extend the approximate methods for computing the induced velocities along straight, radial blades of horizontal-axis turbines and propellers, as described by Wood et al. (2021), to the whole rotor plane for blades that are swept. This extension also provides accurate approximations of the induced velocities for a swept or nonswept rotor anywhere in the plane of rotation. The basis for the analysis is the approximations derived by Wrench (1957) for the exact Kawada–Hardin (KH) equations for the velocity field of a doubly infinite vortex of constant pitch and radius. These equations involve infinite sums of products of Bessel functions and their derivatives, which cannot be evaluated quickly and accurately in many cases. The extension of the approximate equations is straightforward because the generalized sums described in Appendix A are not significantly more complicated than those needed to determine the velocities along straight and radial blades. Furthermore, blade element/momentum theory (BEMT) requires the induced velocities in the axial and circumferential directions. For all parts of the present analysis, there is a close connection between them, so only the axial velocity was considered in detail.

Comparison to the exact solutions showed the extended approximate equations are accurate. For some operating conditions, sufficient accuracy may be obtained using the first or A − term. Two different third-order or C − terms were used. The first extended the expansion of the Bessel functions by Wrench (1957) to the third order, and the second was the leading-order approximation to the alternative treatment of Okulov (2004) and Okulov and Sørensen (2020). The two expressions give similar results, but Okulov’s expression is more accurate.

The application of the analysis in a blade element/momentum theory (BEMT) code must consider the differences between the doubly infinite results from the KH equations and the singly infinite vortices trailing from the junctions of the blade elements. For an angle θ between the origin of the trailing vortex and the control point, the KH equations give the sum of the velocities at ±θ. An additional, approximate equation for the difference in the velocities was derived from the Biot–Savart law. This equation contains incomplete elliptic integrals, which are computationally expensive to evaluate, so the first term in their low-modulus expansion was tested. This term was found to lead to accurate values of the induced velocity for small values of θ, but the accuracy reduced as θ approached π, which may well be applicable for the velocity induced in the hub region by the tip vortex.

The most interesting results from the previous section are those in Figures 7, 8 for small |θ|, where u(θ) is much larger in magnitude than u(−θ). This suggests that the induced velocities when the blade is swept forward, whereby the tip is “ahead” of the remainder of the blade, will be significantly different from those due to sweepback. Furthermore, the large magnitudes of u(θ) occur for small θ, suggesting that the LL representation of the blade may need revising to account for finite blade chords.

This analysis was developed for a single, swept blade, but it can be extended to multiple swept blades or to the more common situation of coned blades. The present equations require the control point to have the same axial (windward) location as the start of the vortex, which is not the general situation for a coned rotor. The difference can be accommodated as follows. A vortex that starts behind (downwind of) a control point can be extended to the point’s axial location. θ is now the angle between the extended vortex and the point, and the present analysis is applied to the trailing plus extended vortex. The velocity induced by the extension on its own can be determined using a simple modification to Eqs 32, 35 and then subtracted. A corresponding subtraction of a vortex segmentwould be the basis for computing the induced velocities when the control point lies behind the start of the vortex.

The further development of the present equations for use in a BEMT code for swept and coned rotors is currently underway. The computationally efficient extension of the equations for ΔI _u for an arbitrary number of blades is being developed together with modifications to make ΔI _u (2π/N) = 0 for any N. A major aim is to determine the necessity of the present method of determining the induced velocities in comparison to the much faster use of Prandtl’s tip loss function.