Windmill Toy

globals [initCount initBlade radius baseRot step-size normSpeed lambda rhoAir areaAir powerT coefP matBlade powerW powerMW]
turtles-own [real-heading areaBlade]

to turtleIt
  clear-turtles
  ;; make turtles, equally spaced in heading
  create-ordered-turtles numBlade [
    ;; this code calculates the area of the blades. if can be easily modified to include any blade
    if shape = "bladea" [set areaBlade 0.5 * 6 * 9 ]
    if shape = "bladeb" [set areaBlade (0.5 * (6 * 9 + 3 * 9)) ] ; slightly off
    if shape = "bladec" [set areaBlade 0.5 * 20 * 9 ]
    ;; move to edge of circle
    fd radius
    ;; turn to face clockwise
    rt 90
    ;; change blade length
    set size sizeB
    ;set color white
    set real-heading heading
    ;turn the turtle in the apparent heading for display only
    facexy 0 0
    set shape bladePick
  ]
end 

to initialize
  clear-ticks
  clear-patches
  ;these two variables control the speed of turtle rotation. This is not physically related to true rotation rate.
  set normSpeed (( windSpeed - 12 ) / 12 )
  set step-size 10 * (1 + normSpeed)
  ;sets the length of the blade
  set radius .5 * sizeB
  ;;background color
  ask patches [set pcolor blue - 2]

  turtleIt

  set initCount count turtles
  set initBlade sizeB
  ;defines swept area of turbine
  set areaAir pi * sizeB ^ 2
  ;defines air density
  set rhoAir 1.225 ; 0.001225
  set baseRot 1.2
  reset-ticks
end 

to go-distance
  power
  ask turtles [
    set heading real-heading
    arc-forward-by-angle step-size
    set real-heading heading
    facexy 0 0
  ]
  refresh
    tick
end 

;;this command calculates the theoretical power of the turbine based on the wind speed and blade pitch:

to power

  let cpnom 0.48
  let genPower (1.5e6 / 0.9)
  ;lambda is the wing-tip speed/windspeed
  let windNom windSpeed / 12
  set lambda baseRot / windSpeed

  let pitchR  pitch * 0.01745329252
  ;;these coefficient definitions come from the reference equations --> they dpend on windspeed, and angle of attack
  ;let lambdaI ((1 / ( lambda + (0.08 * pitch) )) - (0.035 / (pitch ^ 3 + 1))) ^ -1
  let lambdaI (1 / ( lambda + (0.08 * pitchR) ) - (0.035 / (pitchR ^ 3 + 1)) )

  ;this inverses the inverse
  ;let lambdaInv (1 / lambda )
  let c1 0.5176 let c2 116 let c3 0.4 let c4 5 let c5 21 let c6 0.0068
  ;set coefP (c1 * (c2 / lambdaI - ( c3 * pitch ) - c4) * e ^ ((-1 * c5 ) / lambdaI) + (c6 * lambda))
  set coefP (c1 * (c2 / lambdaI - c3 * pitchR - c4) * e ^ ((-1 * c5 ) / lambdaI) + c6 * lambda)
  ;this scales the coefP in accordance to the efficiency changes that the turbine could possibly experience by changing the number of blades. These values max at 3 blades and drop off. they are not yet derived from actual data.
  set coefP (coefP * ((-0.00003 * (numBlade ^ 5) + 0.0006 * (numBlade ^ 4) + 0.0012 * (numBlade ^ 3) - 0.0816 * (numBlade ^ 2) + 0.4151 * (numBlade) + 0.3455) / 0.93011))

  ;this corrects for the material, (non-physical toy scalar)
  ;ifelse material = "balsawood" [ set matBlade 0.8 ]
  ; [ ifelse material = "aluminium" [set matBlade 0.7] [set matBlade 1]]
  ;this is the power calculation, divide by 1e6 to make it in Megawatt
  set powerW (rhoAir * areaAir * .5 * windSpeed ^ 3)
  set powerT coefP * powerW
  set powerMW (powerT / 1000000)

  ;set powerT matBlade * (rhoAir * areaAir) * (kP * coefPnom * (windNom ^ 3)) * (1 / baseRot)
end 

to refresh ;for debugging, strange blade behavior
  ask turtles [set size sizeB set shape bladePick]
  set areaAir pi * sizeB ^ 2
  ;sets the length of the blade
  set radius .5 * sizeB
  set normSpeed (( windSpeed - 12 ) / 12 )
  set step-size 10 * (1 + normSpeed) * coefP
  if numBlade != initCount or initBlade != sizeB [
    turtleIt
    ;go-distance
  set initCount count turtles
  set initBlade sizeB
  ]
end 

;; This is the blade spin procedure.  It moves the blade turtles to the next point
;; on the circle, the given distance along the curve. Code copied from the model library

to arc-forward-by-angle [angle]
  ;; turn to face the next point we're going to
  rt angle / 2
  ;; calculate the distance we'll have to move forward
  ;; in order to stay on the circle. Go there.
  fd 2 * radius * sin (angle / 2)
  ;; turn to face tangent to the circle
  rt angle / 2
end 

;; these are some relevant equations for air foils that were not included, but could be useful for understanding
;air density 1.225 kg/m3
;Coef = 2 * pi * pitch ;;where pitch is the angle of attack in radians
;Lift = 1/2 x coefL x densityWING x windSpeed squared x wing area
;Drag = drag + dragInd = Cd = D / (A * .5 * r * V^2) + dragInd
;dragInd=Cdi = (Cl^2) / (pi * AR * e) ;; for delta wing, AR = wingspan/averageChord ;; use e = 0.7 for approximation
;Power in the Wind = 1/2 x densityAIR x swept area x windSpeed cubed
;Power in turbine = torque * omega = torque * tangentialSpeed / radius ~= (lift - drag)*velocity
;acceleration of turbine = (lift-drag)/(mass)
;probably needs to be fudged so that turbulent air slows down the foil --
;;;Re = (airDensity * V * L) / mu ; mu = 1.784 x 10e-5, in case of trying momentum simulation


;material from public domain and simulink documentation: 
;
; Siegfried Heier, “Grid Integration of Wind Energy Conversion Systems,” John Wiley & Sons Ltd, 1998, ISBN 0-471-97143-X

;This work was created by Jonathan Harrington and is licensed under the Creative Commons Attribution-NonCommercial 3.0 United States License.
;To view a copy of this license, visit http://creativecommons.org/licenses/by-nc/3.0/us/
;or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.