modification to empirical region formula

CHANGES.md

@@ -1,6 +1,10 @@
 # Changelog
 
 
+## [0.2.3] - Nov 2021
+
+Minor change in how the empirical region is handled based on suggestion from Kenneth Lønbæk.  In effect it is Buhl(F = 1)*F.  This forces CT -> 0 as F -> 0 at the tip.
+
 ## [0.2.2] - Aug 2021
 
 - Two bug fixes for AlphaReMachAF functions

Project.toml

@@ -1,7 +1,7 @@
 name = "CCBlade"
 uuid = "e1828068-15df-11e9-03e4-ef195ea46fa4"
 authors = ["Andrew Ning <[email protected]>"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 FLOWMath = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3"

src/CCBlade.jl

@@ -271,16 +271,13 @@ function residual(phi, rotor, section, op)
         if k >= -2.0/3  # momentum region
             a = k/(1 - k)
 
-        else  # empirical region
-            g1 = F*(2*k - 1) + 10.0/9
-            g2 = F*(F - 2*k - 4.0/3)
-            g3 = 2*F*(1 - k) - 25.0/9
-
-            if isapprox(g3, 0.0, atol=1e-6)  # avoid singularity
-                a = 1.0/(2.0*sqrt(g2)) - 1
-            else
-                a = (g1 + sqrt(g2)) / g3
-            end
+        else  # empirical region. Not Buhl's correction but instead uses Buhl with F = 1 then multiplied by F.  
+            # (Buhl(F = 1)*F).  The original method does not force CT -> 0 as f->0.  This can be problematic if 
+            # using a or CT directly as a design variable.  Suggestion courtesy of Kenneth Lønbæk.
+            g1 = 2*k + 1.0/9
+            g2 = -2*k - 1.0/3
+            g3 = -2*k - 7.0/9
+            a = (g1 + sqrt(g2)) / g3
         end
 
         u = a * Vx

test/runtests.jl

@@ -137,15 +137,15 @@ theta = atan.(pitch./(2*pi*r))
 chord = 0.10
 
 
-function affunc(alpha, Re, M)
+function affunc2(alpha, Re, M)
 
     cl = 6.2*alpha
     cd = 0.008 - 0.003*cl + 0.01*cl*cl
 
     return cl, cd
 end 
 
-sections = Section.(r, chord, theta, Ref(affunc))
+sections = Section.(r, chord, theta, Ref(affunc2))
 
 
 # --- inflow definitions ---
@@ -203,7 +203,7 @@ eff, CT, CQ = nondim(T, Q, Vinf, Omega, rho, rotor, "propeller")
 
 alpha0 = -3*pi/180
 
-function affunc(alpha, Re, M)
+function affunc3(alpha, Re, M)
 
     cl = 6.2*(alpha - alpha0)
     cd = 0.008 - 0.003*cl + 0.01*cl*cl
@@ -212,7 +212,7 @@ function affunc(alpha, Re, M)
 end 
 
 
-sections = Section.(r, chord, theta, Ref(affunc))
+sections = Section.(r, chord, theta, Ref(affunc3))
 
 Vinf = 5.0
 Omega = RPM * pi/30 
@@ -244,7 +244,7 @@ Q = sum(r.*out.Tp*(r[2]-r[1]))*B
 
 theta = atan.(pitch./(2*pi*r)) .- 3*pi/180
 
-sections = Section.(r, chord, theta, Ref(affunc))
+sections = Section.(r, chord, theta, Ref(affunc3))
 out = solve.(Ref(rotor_no_F), sections, ops)
 T = sum(out.Np*(r[2]-r[1]))*B
 Q = sum(r.*out.Tp*(r[2]-r[1]))*B
@@ -364,8 +364,8 @@ for i = 1:ntsr
     cpvec[i], ctvec[i], _ = nondim(T, Q, Vinf, Omega, rho, rotor, "windturbine")
 end
 
-cpvec_test = [0.02350364982213745, 0.07009444382725231, 0.13891408294580307, 0.21795999362154536, 0.30793657066245234, 0.39220453169574504, 0.43584242088313013, 0.4590001212905916, 0.4695346630018769, 0.4680559467091491, 0.45853627726342805, 0.44346560401350865, 0.4247660950655427, 0.4031517121734798, 0.37850875685858926, 0.3506421939620393, 0.31947670368255254, 0.28492212475211376, 0.24684802187254656, 0.2053516821754716]
-ctvec_test = [0.1234606206655405, 0.19135957414241966, 0.2753772801035205, 0.3637320188268945, 0.46087282130376256, 0.5702751788685602, 0.6510892321304749, 0.7117745993566226, 0.7600809643657047, 0.7994008741637466, 0.8332719078069747, 0.8638031244208144, 0.8924359419919573, 0.9199569340127343, 0.9465496984524919, 0.9723679911393518, 0.9975143642860824, 1.0220510627832238, 1.0458478359218666, 1.0685364647098319]
+cpvec_test = [0.02350364982213779, 0.07009444382724848, 0.1389142676582008, 0.21796006904596332, 0.30793648085907793, 0.39220377428841097, 0.4358126426438294, 0.45890373500510734, 0.4692861367774551, 0.4675297063104674, 0.45763823228659106, 0.442190136700292, 0.4231819091376136, 0.4013611296914438, 0.37661228050658546, 0.3487266071360946, 0.31761016554953614, 0.2831513131286119, 0.24519843834864297, 0.20383577257448457]
+ctvec_test = [0.12346062066554207, 0.19135957414241508, 0.2753773487370145, 0.3637320163247399, 0.4608728263820876, 0.5702746116812709, 0.6510654098500409, 0.7116957736180881, 0.759888802150446, 0.799015527437428, 0.8326366210606186, 0.8629179420352486, 0.8913475399157195, 0.9187322806175379, 0.9452557920123842, 0.971065072969849, 0.9962534829381495, 1.0208678404660274, 1.0447602003736, 1.067546263425622]
 
 for i = 1:ntsr
     @test isapprox(cpvec[i], cpvec_test[i], atol=1e-3)