Binary  Flutter  Solution  for  Fluid  Power.               Please view in Internet Explorer to see figures



S.P. Farthing, Applied Mathematician,Wing’d Pump Associates   975 Tuam Rd. North Saanich  B.C. V8L 5P2 Canada

Journal of Aerospace Engineering  2018 Vol. 31, Issue 3




The stability of a foil with its ¼ chord center of pressure trailing a pitch axis sprung in heave is solved algebraically to help design a fluttering windmill and perhaps watermill. Its flutter mode and frequency/windspeed do not depend upon its total mass or spring rate.  All contours of this ‘reduced’ frequency in the pitch inertia & imbalance plane pass through a nexus whose total inertia and imbalance are as if just the virtual mass were at the ¾ chord point, with a mode of feathering in the apparent wind at this aerodynamic center. The high frequency flutter amplitude ratio is symmetric in pitch inertia about the nexus.  Similarly from the second factor in its pitch damping, each contour passes through another nearby simple point as if twice its Theodorsen factor times the virtual mass were a ¼ chord divided by this factor behind the ¼ chord. So twice the virtual mass at midchord gives the zero frequency inertia and imbalance “midpost  furthest away from the nexus.  Small trail makes the imbalance greater at the midpost than the nexus so as to slope the zero frequency line downward. Then the imbalance required for quasi-steady  flutter decreases with pitch inertia, even below nil beyond the nexus. Trail also bends the gate of simple points to pass some low frequency contours very slightly below the midpost to locally lower the flutter boundary.  For an oscillating windmill the net virtual mass reaction stiffens heave, opposed by the circulatory lift in flutter since its pitch and heave are necessarily partly in phase. Such new results, and a water flutter demonstration show a practical semi-rotary water blade would need a geared-up pitch flywheel for sufficient inertia to flutter well.  Whereas a wing is so much heavier-than-air it has enough structural pitch inertia to flutter and so pump easily.


Keywords:  mass ratio, stability contours, Theodorsen function



1. Introduction


 It was recognised early in the study of aircraft flutter that its spontaneous phased oscillations of the two ‘binary’ degrees of freedom were being powered by the airstream.  Duncan (1948) even built a heaving ‘engine’ that articulated a balanced foil to pitch and heave (or “plunge”) 90º out of phase to pedantically show this (and nothing more).  In fact to safely tap the highly variable power of ambient flows requires exploiting both free amplitudes of flutter (Farthing 2012).  Our FlutterWing’dPump (FWP) originated in 1978 after the 1976 BBC broadcast of Pocklington School Young Scientists’  fluttering lab models promised  better wind waterpumping  than rotary multiblade windpumps,  especially for developing countries. 

        A fluttering windmill must be designed to have a powerful instability to large amplitude over a wide range of moderate windspeeds whereas  aircraft flutter simply mustn’t start before the never-exceed speed.  Previously unknown was 3D flutter restabilisation in high winds, a key advantage of a fluttering windmill hypothesised  and verified on models  in 1978, then  numerically in 1980, by a full-scale  FWP in storms in 1990 and finally algebraically. (Farthing 2013).

          As early aircraft increased in speed,  flutter critical onset speeds were surpassed with many crashes from the destructive oscillation of control surfaces unless aerodynamically and mass balanced, or of wings unless torsionally stiff.   Yet flutter has scarcely been a problem in marine hydrofoils.  The ratio of foil mass to the virtual mass m of the circumscribing fluid cylinder is much lower because of the 700 times higher density.  Solid steel (Fe) hydrofoils would be very understressed but still below unit mass ratio.  Flutter calculations for hydrofoils by eminent aeroelasticians (Ashley et al 1959) reinforced their experimental evidence of a lower flutter limit of roughly equal real and virtual mass.  Their blades twisted elastically behind the ¼ chord center of pressure so divergence frustrated analytical solution and full understanding of the stability in water against flutter.

     So here the simplest non-divergent pure flutter linear model will be solved to clarify the mass ratio effect. To further motivate the algebra to come,  the next section develops the conceptual niche of flutter pumping vs. heave engines vs. wind and water turbines.


Figure 1 Flutter Wing Pump schematic


2. Flutter & Oscillating Capture of Wind and Waterflow Kinetic Energy


The ease of bird flight, and the speed and leaping of fish, show a remarkable efficiency to oscillating propulsion,  confirmed metabolically.  But for windmills, such an efficiency of power / drag/  flow speed V  bears only on total kinetic capture by an array,  and on capture per structural cost of the drag (Farthing 2007),  most significant in water channels where the stream is not semi-infinite, and very high forces produce the power at very low speeds.  For the higher V but lower power density wind,  the variation of linear oscillation power as amplitude squared means sweeping the maximum area for the area of the wing is more cost effective. Bio- mimicking the wing of the fastest flapping bird was worse than ill-founded because the swift’s sweptback lunate planform has a destructive dynamic divergence rather than a safe cessation of flutter at high airspeeds (Farthing  2013 ) . So actually bird wing flapping is the ultimate in active flutter suppression and pitch control

       As in Fig 1, a semi-rotary 3D oscillation of a wing from vertical about a low streamwise axis close to the ground avoids, like the Vertical xxis wind turbine (Vawt),  the high tower mounting of the rotors, gearboxes and generators of the standard Horizontal Axis (Hawt).   Inverted in waterflow,  these alternatives would keep the bearings and power conversion out of the silty,  if not salty,  bio-corroding water difficult for man to access.  Ideally the Vawt’s power peak,  narrower than the Hawt’s,  might capture a more acceptable fraction of the tidal power than of the wider wind spectrum (Farthing 2009).  But the fouling of a hydro-Vawt’s blades and struts raises the foil drag D to which its power is highly sensitive , and the straight blades of a (grid) synchronous Vawt would vibrate badly from cyclic stall in peak tides. Whereas a semirotary blade projects from just a bit of axle below a floating base,  does not demand  low D and minimises stall (Farthing 2008).  The recent proliferation (Young 2014) of theoretical water heave engines only oscillate in 2D pseudo-heave with swing arms, bearings, roller chains etc often in the water.  Pure heave needs linear bearings, impractical outside or underwater.  Free-in-pitch flutter would be bidirectional for the tides whereas underwater articulations (and Hawts) need to shift their pitch 180º or yaw.  The increase of the wind with height lessens the semi-rotary blade nominal angle of attack decrease with radius, whereas the weaker tidal V drop with depth below the surface worsens it a bit.

          Any reciprocating machine of frequency w is stress limited in its design flowspeed Vd. A wing length R scales geometrically for a given flow’s peak ‘live’ useful wing loading  w so the structural mass per unit wing area should be proportional as DwR /s  where D/s  is the weight/strength  ratio of material density D to fatigue limit stress s. Then with the blade’s speed wR scaling kinematically as design flowspeed Vd for the optimum tip speed ratio or reduced frequency k,  its inertial acceleration w2R multiplies the mass loading for a parasitic inertial loading DVd2/s times the design useful live loading  w. So independent of fluid density there is a scaling maximum for  DVd 2/s.  This DVd 2/s is just the scale of inertial self-loading Dw2R2 to the limit stress s.  It indicates a serious limit on Vd for oscillating windmills (and orthinopters).

        In contrast in rotation, the main inertial stress is centrifugal, or just benignly tensile in Hawt blades.  So high rpm lowers the structural cost per unit power and also the gear up to a generator to make high flowspeed sites the most economic for wind turbines.  In high wind Vawts centrifugal dominance enables a catenary blade with gradual stall and fewer struts to drag.

        A slow oscillating prime mover for light flows could avoid a high ratio gearbox by reciprocating a pump, whereas Hawt’s are very poor and Vawt’s worse at cranking piston pumps,  especially single-acting. Their ideal torque,  varying as windspeed squared,  severely mismatches the mainly angular torque variation of cranking a single-acting pump with constant head.  The difficulties are most intractable in the more widely and rapidly variable wind with the pump constrained deep down a well.  (Dixon 1979) showed a rotary fanmill’s useful pumped work is only 1/10 of the ideal capture of the annual wind energy flux through such a swept area.  Its 20 or more blades maximise the starting torque,  but still not enough to turn over its pump crank in the most efficient wind for its stroke.  Wind tunnel flow visualisation and a flow solution showed for the first time (Farthing 2011) how all of the kinetic energy of the reaction flow to the torque is lost in the wake by centrifuging.  Whereas an oscillating actuator approximates a contra-rotating windmill (Farthing 2010) reducing this loss.  The net ten fold was the primary motivation to develop a pump oscillated by the wind. But many oscillation mechanisms would be worse than Hawts in self-starting against the fixed pump head or handling the mismatch of flowpower as the cube  of  flowspeed V  vs. piston  pump power as just the frequency w. With fixed head, and articulated pump stroke,  their frequency w  V 3 to absorb the best fraction of the cubic windpower but then inertial reactions on the articulation mechanism would increase as V 6!

     Instead at high Vd, flutter’s inherent dynamic balance of inertia, imbalances and flow at near-constant frequency but variable amplitude can be non-linearly converted into a highly variable pump stroke as in Fig 1. (Farthing 2014) composed full amplitude equations of motion from the 2D fluid (Kochin 1964) and 3D rigid body dynamics. Flutter’s non-linearities prove favourable to significant economic extraction of moderate wind power for pumping by the “FlutterWing’d Pump”.

    The wind doesn’t have any merciful ratio of maximum to design average speed like a tidal flow.  But it was speculated that as a resonance of a fixed roll frequency with pitch frequency as windspeed ,  the flutter of  3D semi-rotary roll (Fig 1) might cease above an upper critical ‘cutout’ windspeed.  3D unswept wings of low imbalance indeed feather quickly and stably to the true wind in a storm.  As the cutout windspeed is approached, the pitch amplitude decreases, containing the power and especially the downwind thrust;  whereas to vary their  pitch amplitude immensely complicates articulations. 

     In very high winds 2D feathering to the heave apparent wind gives zero effective pitch stiffness and slow flutter/ divergence. But the basic linear 2D pitch and heave flutter analysis below does explain the ready start of the FlutterWing Pump in light winds yet the stability of heavy ferrocement or even solid steel blades in waterflow.



Figure 2 : Section of Symmetric Airfoil Heaving across Wind  and Pitching ahead of quarter chord

As In Fig 2, consider a unit length section of an infinite (vertical) symmetric airfoil of chord c, whose ¼ chord center of pressure (c.p.) point trails by a distance ec the pivot point free to pitch at any angle g to the stream velocity V.   This trail  ec≥0 eliminates a divergence complicating typical wing bending-torsion flutter.  For the fluttermill of Fig1 small e is structurally desirable so the bottom of the fabric on frame wing can pitch around a tubular steel  axle cantilevered above the pendulum.


Pitch elastic stiffness is not needed on the FWP and would interfere with bidirectionality of a tidal flutter mill. It is also absent in cantilevered ‘spade’ rudders on boats,  if not on all-moving aircraft rudders . But the FWP pivot point is sprung to heave cross-stream with coordinate h.  Use the coordinate vector Q= (g, h/c) and k=wc/V.  Let kn =wnc/V  where wn is the natural roll frequency in heave of only the virtual mass in stationary fluid.  The small amplitude equations of neutral stability oscillation will be non-dimensionalised below as 


                                                                                                                        (–k2A+ikB+C+ kn2 E) Q=0            (1)


with the real matrices,  A inertia, B aerodynamic damping, C aerodynamic stiffness, and E elastic stiffness to be elaborated below. This symbolic solution will identify exactly what (complicated) matrix determinants will be needed.  If the cross determinant is


                                                                         [A,B]=A11B22+A22B11-A12B21-A21B12   then  |A|=½[A,A]               (2)


Then the nil determinant in (1) for  neutral oscillatory stability in powers of k expands to


      k4|A|   -ik3 [A,B]  - k2(|B|+[A,C]+[A,E] kn2)   + ik ( [B,C] +[B,E] kn2)    + ( |C| + [C,E] kn2  + |E| kn4)=0   (3)


where all the crossed out terms will be shown to vanish in this simple case.  For example |E| =0 is because the only elastic spring is E22 in heave. Equating imaginary odd power parts                                                                                

                                                                                       k2[A,B]= kn2[B,E]                                                        (4) 

 (or instead  a flutter/divergence k= kn=0 at infinite windspeed).  Equating even power real terms multiplied by [B,E] to eliminate kn2 gives  

                                                          {|A| [B,E] - [a,E][A,B]} k4  +  { [C,E][A,B]  -  (|B|+[A,C]) [B,E]}k2 = 0 (5)


Either  k= kn = 0 again or the neutral stability criterion is …    

                                                          k2 = [A,B][C,E]- ( |B|+[A,C]) [B,E]      /        [A,B][a,E]- [B,E]|A|   

This ratio of the difference of 4 way products has hitherto defeated meaningful analytic treatment of binary flutter.  Here each product has an E22 factor  so  k the conventional positive root is independent of this heave spring


                                                       k2 = [A,B]C11 - ( |B|+[A,C]) B11      /        [A,B]a11 -B11|A|                         (6)


The naive Routh criteria for stability of determinants not depending upon k is lesser k2 than (6) are unstable; greater, stable.




(6) can be simplified by showing  the flutter contours and mode depend on just inertia and imbalance, and evaluating each cross-determinant in its most convenient moment axis in terms of the key distance factors in the pitch damping.

       For a thin airfoil, potential theory gives a virtual ‘added’ mass of m= ¼prc2 as in the circumscribing cylinder of fluid (of density r) centered at midchord  for virtual imbalance ¼m about the ¼ chord center of pressure.  But the virtual intrinsic pitch inertia about midchord  is mc2 /32,  a factor of 4 less the real inertia of a solid (ice) cylinder’s mc2 /8 . So the virtual inertia about ¼ chord is  3mc2 /32,  just half  of a solid cylinder. These ¼m and 3mc2 /32  will prove insufficient for flutter which needs significant real additions.  

       Use m to non-dimensionalise the total of virtual added to real mass, as pm for p>1,  total pitch axis  (dynamic) imbalance as  mxc,  and total inertia as mjc2. The traditional mass ratio is  p-1. Then |A|= pjx2 is  p /mc2 times the minimum inertia about the center of mass > mc2/32.  So |A|= pjx2>p/32 (equal when the real mass is concentrated at midchord). 

         Let the spring  restoring force be S per unit heave.  More real mass at the axis changes p by Dp but not j and x.  For flutter’s sinusoidal h at frequency w,  the increased inertial force hm Dp may be matched solely by an increased spring force hDS with no other terms changed.  But as above k2 is unchanged by negating just the DE22  proportional to DS  alone so that leaves k2 unchanged by just the Dp, so p must cancel out of (6).  Since the mode or amplitude ratio and phase shift can be determined from the pitch equation ahead (12) which does not involve S or p, then the flutter mode and k are the same for all mass distributions with the same pitch imbalance and inertia about the physical axis, regardless of the total mass p and spring constant S. Exactly the same argument applies to the 3D semi-rotary windmill  that the roll inertia and stiffness do not enter the k equation. Here to give w and kn too, the determinants are given for general p which indeed cancels out of the (19) for k2. 

      The heave real inertia and spring difference is balanced by the virtual mass reaction and the circulatory lift L acting at the ¼ chord point

                                                                             Sh +(p-1)mh”  =  m { -h”+gV +gxc } +L =p0mh”       (7)   


  -gV  comes from a flow Coriolis acceleration in the frame of the foil and negates the virtual roll inertia h  in the feathering mode g=h’/V.  That makes the heave force of the fluid virtual mass m negative for quasi-steady windmilling flutter (vs. propelling) nominal normal velocities dominated by gV over h’.  But the lift L can also have a component in phase with h.   For a given j and x and so oscillation amplitude ratio,  write the sum of fluid  terms as p0mh”  where if  p0>0 the net fluid force is stiffening in heave.    Then

       S=mwn2=(p- p0-1)mw2   or  k2n=(p-p0-1)k2 (8)


     The pitch moment balance per unit span about the pitch axis is                  m {gj c2 +gVcq -hxc} = -ec L  (9)

  where  q=e+½  or the distance in chords from axis to the ¾ chord 

      Fig 2 shows the nominal apparent wind N at angle j at the ¾ chord  aerodynamic center

                                                                                                      j=g-h¾’/V  if  h¾=h-cqg  so  j=g-h’/V+ cqg‘/V (10)


cqg/V is sometimes called the effective camber b  as the lift of an arced foil subtending 4b in steady flow is proportional to  j=g +b  as the angle of attack of the nominal apparent wind N at ¾ chord.  All nominally,

  U=Vj= Vg -h’+ cqg  is the normal flow at ¾ chord, pcU  the circulation and the lift is L=4mVU/c  =4mV2j/c 

        When  Q are small and sinusoidal ,  the unsteady wake corrects these nominal values  by the Fig 3 wake factor function of reduced frequency                        T(k) =K1(ik/2)/{K0(ik/2)+K1(ik/2)} =F-iG

 in terms of modified Bessel functions of the second kind with Fat large k and F(0)=1.  G is small but has infinite slope  at G=k=0, and decays slowly at large k. and so will be carefully added in a further paper (Farthing 2017).  For straightforward  real determinants and initial simplicity here ignore -iG which only gives F  a phase lag of at most .27 or 15° around k=.44. This is the sole and temporary approximation T@F,  apart from linearity.


Figure 3  Theodorsen Function Real and Negative Imaginary Parts

          So extend h/c to the complex domain as (Implied Imaginary Im part of eiwt ) H and  likewise

g=(Implied Imaginary Im part of eiwt ) G where  G=g0eiq  where g0  is the amplitude of g  and q its phase lead over H .  Theodorsen thus extended L=(Implied Im part of eiwt  ) 4mV2T {G+ikqG-ikH }/c 


Then the roll equation (7) divided by mV2eiwt  /c  extends to         G{k2x-ik(1+4Tq)-4T}+H{ k2n- pk2+4ikT}=0   (11)

In the very high wind limit of kn and k and L going to zero,  quasi-steady flutter/divergence mode is  G=ikH=iwh/V  i.e. feathering in the apparent wind j=g -h’/V=0 (the same at all chord positions in view of the low k).   Angular rolling instead of linear heaving cannot have exactly this neutral stability at k= kn =0, and can actually be stable in high winds if  [B,C] is positive rather than nil here.[Farthing 2013]

    Henceforth for T=F, the net pitch damping  is mgVcq+4mVFecqg=4m FqycV gwith y=e+¼/F,  so the pitch moment equation (8) non-dimensionalised by  mcV2 eiwt   extends to    

                                                             G{-k2j +4ikqyF+4eF}+H{k2x-4ikeF}=0                                                    (12)

     Aerodynamic balance e=0  nulls the leading pitch stiffness  so the dominant pitch balance is the aerodynamic damping equalling the cross-forcing by roll via dynamic imbalance or 4ikqyGF+ k2xH=0,  so in the feathering limit G=ikH ,  x= 4qyF=½ for neutral quasi-steady stability.  This net  is twice the virtual  x for an imbalance ratio of 1. At such phase shift q=p/2,  the net energy-conserving dynamic imbalance is conversely damping the roll. In general,  from (11) and (12) 


[A,E]=A11=j    [B,E]=4Fqy=(e+½)(eF+4)   [C,E]=4eF   [A,C]=4F(ep-x)   [A,B]=4F{pqy-(q+y)x+j} (13)


(4) gives the flutter frequencyw  only involving F through y    w2=S [B,E] / m[A,B]=Sqy / {pqy-(q+y)x+j}  (14)


 (4) &(14) are also k2{pqy-(q+y)x+j} =qyk2n   Substituting (8) k2n=(p-p0-1)k2   cancels p for  (p0+1)qy+j=(q+y)x   (15)


 Adding -ec times the lift Eqn(7) to the moment Eqn (9) eliminates ecL  to give the ¼ chord c.p. moment as m  times

                                                                 g”{jc2-ec2x}-hc{x-ep}+cgV/2+ ecwn2h=0    (16)

    This exchanges (complex)  L and its  g, g and h terms for a cross stiffness ecSh= ecwn2h.  The lack of wind pitch stiffness in this c.p. moment (13) means its first pitch ‘inertial’ term (only 0 at j=ex)   must be balanced by some of the cross inertial fling or stiffness by heave being partly in phase with pitch.  So flutter’s phase lead of pitch over heave is ≤90º, actually significantly less at large amplitude [Farthing 2014].  This is what makes the leading g term in the circulatory L partly in phase with h,  so opposing the negative net heave virtual mass effect in (7).

     The eliminated eL is often misidentified as the circulatory moment.  However the remaining ½mcgV  in (14) from (9) is explicit damping and so must be circulatory as potential flow is perfect (without dissipation).  In (9)

mgVc(e +½) arises in the subtraction of a circulatory broadsiding moment in the ¾ chord apparent wind (at angle j) from the inertial broadsiding  moment in the pitch axis apparent wind and so generalises as  mg(Vcosg+h’sing)q  [Farthing 2014].   Non-dimensionalising the linear ¼ chord moment (15) by  mcV2eiwt   

                                                                                                        G{-k2(j-ex)+½ik}+H{k2 (x-ep)+ekn2} =0  (17)


So B,C]=|C|=0         [A,B]= ½p+(1+4Fq)(ep-x)+4F(j-e x) |B|=2F       |B|+[A,C]=2F(1+2ep-2x)  (18)

         As a recap below are the non-dimensional  matrices  with one of the first two rows used at a time. The second ¼ chord center of circulatory lift row was just used

         A= B= C=  E=         


Whereas (13) |A|=pjx2, [A,B]=4F{pqy-(q+y)x+j}=4F{(j-qx)(j-yx)+(pjx2)qy} /j ,  [A,E]=j , [B,E]=4Fqy,[C,E]=4eF 

were the general p results with the first row moments about the physical pitch axis .   For the numerator of (6)  use the top row form of [A,B], A11,B11,C11 but the 2nd&3rd row form of [A,B] for the denominator, and p=0 to avoid  going through the cancellation of the p terms, to get the remarkably simple


                                                  k2 (j-qx)(j-yx)  =  4eF(j-ex)+½x-2Fqy           (19)


This may be the first time the high extent of algebraic cancellation in the difference of products of determinants in (6) has been demonstrated. At large inertia  j , k2=4eF/j  or just the then slow pitch-only aerodynamic reduced frequency.  For low k  For low k, top of Eq. (6) ≈ 0, so from (15)  p0» (x-q)/e, and left of Eq. (19)≈0»4e (j-ex)+.5x/F-2qy

or 0»4e (j-ex)+½/F-2qy so F1-¼pk   gives  je+{⅛ -e2}x q(e+¼)+pk(x-q)/32@ 0,  lines rotated as  ek.

From (4), (13) and  (19)   kn2= k2[A,B]/[B,E]= {4eT(j-ex)+ ½x -2Fqy}[1/qy+ (pjx2)/ (j-qx)(j-yx) ] / j  (20)


giving the flutter V2 = Sc2/  mkn2 where S is just a factor and the p influence is just through |A|= pjx2>0.  {}>0 above the general k  line {}=0 where k2 and V2>0 for  j<yx or j>qx.




Consider the inertia and imbalance for just the virtual mass concentrated and relocated to the ¾ chord point. Then   sddxn=jn/xn=q=e+½  vanishing the numerator as well as the denominator of (19) to solve it for all k and F.  Substituting in the physical axis pitch moment (12) the mode is    H=G(q-i/k) or  h’3/4=gV so  j=0=Fj =L= 0, whatever the complex F.  Thus all the exact (complex Theodorsen)  k contours in the j, x plane must also go through this  “nexus.  From (14) its p0=0.  As  k­, this mode HGq effects pure pitch about ¾ chord which emerges out of the nexus along its ray (from the origin) j=qx  from (19) and (12).  The aerodynamic pitch damping moment about the ¾ chord indeed then vanishes,  so this effective mode is energetically ie stability neutral (whilst not divergent like the unitary pitch about a ¾ chord physical axis).

  (19) will also be satisfied when j=yx instead of j=qx intersects the RHS=0 line at  xg=2Fy=2Fe+ ½ and jg=2Fy2=xgy.  Similarly from (12) its mode is    H=G(y-i/k) or  h’y=gV  (for local apparent wind angle jy=0).  Here the larger imbalance (e>0) and smaller inertia are most simply as if from a point mass of 2mF at  y=e+¼/F.  From (14)  p0=2F-1.  As k­& F↓½, this 2mF, jg & xg fuse into the ¾ chord  nexus. The locus of (jg,xg) or the “gate” all the contours cross,  is for e=0 just the straight k =0 line at x=½ between the (F=½) nexus jn=¼  and F=1 ”midpost.” jg=⅛.  The gate is (xg-½)jg=ex2g, curved by e,   below the k=0 line.

            To summarise, the  nexus and gate  and their modes are as if

 point m    at Nexus    q=e      xn= q     jn=qxn    2(jn-exn)=q    p0=0   j¾=0   H=G(q-i/k)  

point 2mF at gate    y=e+¼/F   xg=2Fy    jg=yxg 2(jg-exg)=y  p0=2F-1  jy= 0   H=G(y-i/k)  (xg-½)jg=ex2g

k2   (21)

   Recasting (19) in terms of these roots as        k2 (j-qx)(j-yx) (¼/F -½)= (j-jn)(x-xg)-(j-jg)(x-xn)   (21)

Since the virtual mass is not a point but distributed with the intrinsic virtual inertia of mc4/32, then to realise the midpost at k=e=0, needs a real point mass of 2m midway between the midpoint and the axis or ⅛c from the axis for a p ≥ 3.  To realise the nexus at e=0 needs less point real mass of  .4m,  but further back for more inertia at ⅝c for a p ≥ 1.4.

           From  (p0+1)qy+j=(q+y)x   (14)   p0=0 when qy-(q+y)x+j=0  as at the nexus and close to its  (F=½) line

 q2-2qx+j=0 which radiates from the  nexus with half the slope of the nexus ray.  To the left at p0>0 the net fluid heave stiffening comes from the dominant negative virtual inertia at high k for the derivative of positive windmill midchord upwash U½. To the right at p0<0 the net fluid heave inertia force comes from the negative stiffness of the circulatory lift from positive ¾ chord upwash U3/4  at  low phase shifts and  k.  p0=-1 on the almost parallel line through the origin.  This effective inertia helps contain the high wind frequency of the 3D semi-rotary windpump of Fig 1 vs. nonlinear stiffening by the pump. (Farthing 2014)

            In Eq. (19) k=0 gives the quasi-steady F=1 flutter boundary line  through the nexus and midpost as

2e(j-ex)+¼x=qy=2ej+2x(⅛-e2) with small e slope about -8e. Comparing these posts at F=1  xm=2e+½=e+xn  and  jmxm2 =(2e+½)(e+¼) <jn=xn2= (e+½)2  until e2=⅛ or e=.35 or the axis 10% ahead of the wing when jm = jn  =.73=⅜+√⅛,  the k=0 boundary line is vertical, and the gate still curved .

       For a given e and k, and therefore F(k), (19) is a (hyperbolic) quadratic in x and  j (tending to parabolic as k­∞); so for instance the two x roots can be plotted vs  j as k contours. Values of  phase shift and amplitude ratio are  tabulated at these x and j so that their own contours can be interpolated. Merged contour plots directly from (18) assuming the limiting F, agree generally quite well. Note there is no guarantee that the same  x and j won’t be generated at different k with different modes so there k contours cross and  mode contours can’t be easily interpolated.

              Figure 4 plots the neutral k Theodorsen contours in the e=0 limit of  aerodynamic balance. On the anticipated x  magnified plotting does show all the contours from the right do cross each other at the nexus at xn  jn= ¼ and then turn  up before x=.4992 and  pass through the xg=½ gate before the midpost  at jm=⅛.  Thus there is a micro zone of recrossing and so double k below the gate. 

                For F=½,  k√2=√(x-½)/│jx│ the square root of the vertical height of the x=½ line over the horizontal distance from the k=∞ x=2j  ray and the k contours are parabolas.  Whilst for F=1 the horzontal distance becomes the geometric mean or root of the product of the distances from the midpost and  nexus rays and the contours hyperbolic.                                              Figure 4 Stability Contours for  Aerodynamically Balanced Foil e=0                                             

            There is no sizeable quasisteady zone to the left of the midpost above the lowest j=3/32 from the virtual inertia about e=0.  In the upper left high k flutter zones the black contours of minimum p of  p>x2/(j-1/32)  show the large implicit p unattainable for a hydrofoil where the real mass p-1 is almost as small  as x and j. 


Figure 5 Heave to Pitch Amplitude Ratio r vs Imbalance and Inertia Factors at e=0

From the pitch equation (12)  H/G=(ji/k)/x  for all real F at  e=0. The first term makes the amplitude ratio r contours just rays around the  nexus ray.  In general r2=│H/G2=( j2-¼k-2)/x2. Using the flutter solution (15),  k-2 also is quadratic in  j  (for fixed F) and so then r is  symmetric in  j at fixed x, eg at F=½, r2=│H/G2 =
¼ x}/ x(2x-1)  which is symmetric about the vertical jthrough the  nexus.  Thus r=½ not only on the pure ¾ chord pitch  k=∞ ray  but also its reflection to finite varying high k.  Regarded  as a quadratic of x as well as j the r2 contours would be hyperbolas asymptoting to the k=∞ line and this reflection.  For F=1, r2={2j(j-3/8)+
x}/ x(2x-1)  which is symmetric about the j= 3/16  vertical midway between the midpost and   nexus and infinite along the xk=0 line.  Figure 5 for F(k)  indeed has nearly symmetric contours of  interpolated amplitude ratio  r whose axis moves from  j=.23  to jas the nexus is approached from above with increasing k,  effectively globally curving the local mirror axis.

      Tanq of the phase lead q,  or the slope in the complex plane of G/H is 1/ 2kj so along a finite k contour, the phase lead is lower with the greater j.  Hence the phase contours are straight(er) crossing the curved k contours downwards in x  to the right of the nexus ray and upwards to the left, as seen in Figure 6.  For F=½,  Tanq = (1-x/ 2j) /√2x-1 where the numerator is the relative difference of  ½  the ray slope x/j  from its identity along the k=∞ line where q has a minimum of zero with the pure pitch about the ¾ chord.   Likewise for F=1,  Tan2q =(1-x/ 2j) (1-x/ 4j) / 2x-1 where the numerator is the root mean of ½ slope differences from the ray through the midpost and  nexus.  As anticipated from (15) q ½p.        Figure 6  Pitch Phase Lead in radians vs Imbalance and Inertia Factors at e=0                 

        The linear aerodynamic power is the circulatory lift heave power - the pitch damping power loss.  With  n=cg0/h0=│G/H│=1/r                           Power/ .5rV 3g0 (h0 +cg0) =     (22)

This power per pitch amplitude  per naïve swept area is optimum at .27,  little changed using T@F henceforth except for the subtraction of  T’s .27 phase lag to reduce the phase shift to q =1.35 (from  1.62) at k=.44,│F│=.73,  r=.88 so H/G= .193+e-.86i  With the power removed in heave, the pitch equation (18) still applies

 H/G=(j-1.19i)/x  requiring  at e=.02 , x= 1.95  j=.6,  more than double the e=.02  nexus jn=.27.  Then the contours show a high onset k=2.2  (indicating a good low starting windspeed). Whereas assuming non-linearly g =+/- 90º,  the optimum k = ⅔ for a tangent blade windmill actuator at q=½p.    [Farthing 2014]. 

       Now fully introduce e>0 with 4eF(j-ex) the new term in (19) for k2with e>0. To first order this term tends to increase the instability as ej. Naively ej makes the pitch motion more naturally harmonic to aid flutter oscillation. The quasi-steady wedge does cross into negative imbalance at the j-intercept  2ej =qy.  Since pitch only is naively critically damped when  ej=F(qy)2 the pitch-only motion at the gate is overdamped by q2/ 2e@ ⅛/e  but at the intercept underdamped by 2qy ≥ ¼  (until critical at nearly vertical interception at  qy=½ or  e=.34 ie the pitch axis 9%  ahead of the foil).                  Figure 7 Flutter Reduced Frequency k  vs Imbalance and Inertia Factors at e=

      In Figure 7 for e=⅛ as, k =0 has down slope about 8e =-1 (intersecting  the previous e= k=0 at j= x=½) with the  nexus √j= x=⅝.  With these the critical k and instability are  increased all the way up the right hand side of  the new nexus ray j>⅝ x for the same j and x. But for x and j on the other left top side above the  nexus ray the k are conversely decreased by its rotation to lesser slope.   Now there is a weak optimum j at the smaller x that give the most unstable highest critical k. From the design for flutter viewpoint,  the reduction of x is welcome (but not so far as structurally unnatural, negative imbalances).  Such j >> jn is desirable to oscillate with low p in the quasisteady zone to the right,  permitting large heave which non-linearly increases the  power (Farthing 2014).  Perturbing for the extra 3D semirotary determinants, flutter can safely cutout when x<xm=2y its midpost value i.e. for the same lower right large j [Farthing 2013] 

           The microzone of double k is now extended below and beyond the midpost. Such contours below the midpost  would indicate  flutter first appears at the lowest contour and then increasing x spreads the unstable range to higher and lower k (until lower k=0 at the q.s. line)  They can be below the midpost because the gate  jg=ex2g /(xg-½)  is curved downwards  and almost identical to the contour km  passing through the midpost at xm= ¾ jm=9/32.  The lower k contours inside of it below the nexus pass leftwards through it and below the midpost whilst higher k contours that emerge from the nexus very slightly below it soon cross it upwards. Very magnified plotting confirms the gate to be very close to the midpost km contour but not identical, because (19) has extra  j2, x and constant terms and very obliquely crosses the gate upwards. In general (19) shows this e>0 k>0 contour through the midpost to have F= k2(e+¼)2/2e which is plotted on Fig 3 for e=⅛ to get km =1.03 and for e=¼ where it has a maximum of km =1.08.  Thus the k@ ½ stability contour is about the most  below the midpost and k= 0 line.     

              From (12) for e>0 at k↓0 H/G­-i∞.  At both F= ½ and F=1 the combination with (19) for r2 can be perturbed for small e  tilting  the original verticals of symmetry  slightly to slope 1/e,  confirmed by exact plots for these limiting F.    In Figure 8 of amplitude ratio r at e=⅛ with F(k), the high k contours are still close to symmetric about a such a sloped line which a priori must go through the nexus.

 Figure 8 Heave to Pitch Amplitude Ratio r  vs Imbalance and Inertia Factors at e=   

The righthand contours basically reflect too, but the close and even crossing k and F confuse numerical contour interpolation for r to the left of the nexus.  With pitch amplitude limited to 90º,  substantial swept area more than four times c requires large r>1.3 , or low x as possible at large j .   In Figure 9 the contours of phase angle are also rotated with the k=0 and k=∞ lines.                Figure 9 Contours of Pitch Lead of Heave in radians at e=   

      Mass balancing without aerodynamic balancing can even be insufficient to prevent flutter at high inertias.  At the intercept x=0  2ej =qy,   the self-wind velocity from the in-phase part of large heave counters the small damping of the small pitch i.e.  erCosq= qy.  The heave damping by its full selfwind is balanced by the quasi-steady lift of the quadature lead component of pitch (and  the Coriolos and effective camber lift  of the velocity of the component of pitch in phase with heave).  From such balances in (11) and (12) with x=0,  q and r2 can be eliminated to yield  this case of (19):  k2j2/2F=2ej-qy.  Feeding k2 into the pitch equation (12),  Cotq=2kj.  As the flutter boundary k=0  q=p/2   is approached along x=0 kr tends up to 1 or feathering.

             At e=√=.35,  k =0  reaches vertical so then unnatural negative imbalance x becomes desirable but is very difficult to obtain  and the required real j contribution becomes absolutely large.  Ultimately in the  e ­∞ limit the  k=0  line tends to x=-j/e which is finally satisfied by the infinite virtual x and j  for  flutter at p=1 without real mass.

             This detailled analysis is a foundation for a further submission reaching exactness by including -iG(k) the small complex part of the Theodorsen function.  There just the extra  G but non-p terms need adding.


6. Flutter  Demonstration  and  Design  in  Water


        In water the x and  j of a thick NACA 0016 of solid steel Fe at e=0  are stable below the gate  in Fig 4,  but for  e=⅛ in Fig 7 are close to the middle of the gate.

          A trailing arm to achieve the desirable  j to the right of the  nexus would be massive.  Instead extending a long pitch axle above the roll axis out of the water to (lowbacklash) gear the pitch up in air by ratio n to a flywheel with inertia mK2 would give sufficient increase Dj=n2K2,  by the square of the stepup ratio,  just like the kinetic energy storage.

       Unless the pitch axis heaves perfectly vertically,  the weight of a dense blade must also be considered.  For instance even a horizontal ferrocement hydrofoil would be severely weight imbalanced in pitch and vertical heave.  A semirotary upright wing in air (Fig 1) has to be over-counterweighted in roll but its tailheaviness forces pitch by roll [Farthing 2013] as well as by roll acceleration.  Given the low design  water current speeds, Vd2 <<gR this 3D gravity forcing effect is dominant in water even at model scale.  Conversely the weight of a blade suspended in water is its own roll spring but must be nose over-counterweighted (again better in air) for net noseheavy gravity forcing

         A very long such pitch air arm allowed a model test of the amount of extra inertia Dj required [Farthing 2015].  Fortunately due to the phase shift q of flutter being less than 90º, the water contact of such a long forward arm is limited up to medium amplitude.  Even at very low uniform e = .02, flutter only appeared above the  nexus inertia jn and improved to  k=.52  at  j= 3/8  vs.  jn = .27    This video might be the first ever public demonstration of  free binary flutter of a hydrofoil.

        A floating  tidal device would need to counter the downstream moment of its underwater flow extractor. Fins angled upwards at the anchored ends would so lift whichever end is upstream and depress the  opposite ‘stern’.  To cancel the pump reactions,  two blades could counteroscillate (assuming sealife can avoid the shearing action.)    The difficult contralinking  pitch shouldn’t be necessary.  But a pitch starting track as in Fig 1 would  definitely be needed  for both ebb and floweach sign of the tide  not being as shifty or gusty as a light wind.  Other design goals besides the essential stepup of each pitch to a flywheel,  are all bearings above or at worst at the water surface  and easy winching the blades against each other to the surface for defouling. To realise any such  Fluttering (Hydro)Foil would require very heavy safety-concious engineering and construction at a shipyard, whereas the density change  of 700 makes the Fluttering Windpump construction lighter than light aircraft (Farthing 2014)


7. Conclusions


All flutter contours radiate from the nexus of total imbalance and inertia about the physical free-to-pitch axis as if from the virtual mass at the ¾ chord aerodynamic center. They also pass through a gate curve from the nexus to a midpost corresponding to twice the virtual mass concentrated at midchord at k=0.

      Both are a substantial increase of imbalance  moment x and especially inertia j about the physical axis above the inherent virtual  x and j. This is how the fluid density (and speed) and material density implicitly affect flutter,  not explicitly by the mass ratio p.  Any wing is sufficiently heavier than air but in a watermill obtaining enough total pitch inertia would require a flywheel geared up from pitch.  Changing from uniform heave to the more practical semi-rotation about the stream axis allows adding gravity forcing and gives a high flow cutout ,  making the FlutterWing’d  pump safe in storms.

      The support of Gifford and Partners of Southampton, The Hamilton (Ontario) Foundation and the Science Council of BC are gratefully acknowledged.





r Fluid density

s   Material Fatigue limit stress

g  pitch angle

g0  pitch angle amplitude

j  Nominal Angle of attack of ¾ chord point

q phase lead of pitch ahead of heave,   

G  complex amplitude of pitch

w circular frequency in radians of phase/unit time

wn  natural no-fluid frequency in heave when V=0

A non-dimensional inertia matrix

B non-dimensional aerodynamic damping matrix  .

C non-dimensional  aerodynamic stiffness matrix

E  non-dimensional elastic stiffness matrix 

c  chord of foil

D Drag

D structural Material density

ec  trail of  quarter chord behind pitch axis in chords

p total mass/virtual

g acceleration due to gravity

G the neglected imaginary part of the Theodorsen function

h   heave of pitch axis

 h3/4    heave of the ¾ chord point

H non-dimensional  h/c complex extended

i  square root of -1

j  pitch inertia/mc2

k  reduced frequency based on chord  wc/V  

km reduced  frequency k >0 of contour passing through midpost

kn  reduced  natural frequency wnc/V

L circulatory lift, L its magnitude

m virtual mass of foil/unit length

n  Pitch amplitude/ H  or 1/r

p ratio of total mass to virtual

N nominal apparent wind at ¾ chord point magnitude N

Q=(h, g)  coordinate vector of heave and pitch  displacements

 q distance of  foil ¾ chord from axis in chords

r   H to pitch amplitude G   ratio or 1/n

S heave stiffness =force / h

F  complex Theodorsen function of k

U  upwash or normal component of the apparent wind N

 V  flowspeed

Vd design flowspeed

w  live flow wing loading

 x  pitch imbalance/ mc 

y distance of  foil midnexus  from axis in chords

denotes time t derivative of  eg h’=dh/dt



Arnold, L. (2001) Extraction of energy from flowing fluids US Patent  6273680 B1


Ashley, H.A., Dugundji, J.,  Henry C.J. (1959 ) “Aeroelastic Stability of Lifting Surfaces in High-Density Fluids”  Journal of Ship Research vol 2 no.4,  10-28


BBCTV 1  18:55 Tuesdays Young Scientists of 1976 Pocklington School  “Oscillating windmills”

Winners: Heat 2 April 20 & Final   May 4


Dixon, J.C.  (1979) “Load Matching Effects on Wind Energy Converter Performance”   Future Energy  Concepts ,  IEE Conference Publication 171,London,  418-421


Duncan, J. W.. “The fundamentals of flutter” (1948) ,  (UK) Ae.Res.Co. R&M  No. 2417, Nov., 1-36


Farthing, S. P. (2007)“Optimal Robust and Benign Horizontal and Vertical Axis Wind Turbines”,  Journal of Power and Energy   Vol 221 No. 7,  971-979.


Farthing, S. P. (2008) “FWP Leading edge smoke”


Farthing, S. P.  (2009)“Vertical axis wind turbine induced velocity vector theory “

Journal of Power and Energy Vol 223  No.2     p104-114


Farthing, S. P. (2010) “Robustly Optimal Contra-Rotating  Hawt” Wind. Engineering Vol 34 No.6, 733-742


Farthing, S.P. (2011)“Analysis of Torque Reacting flow through Starting Windmill”   Appendix toBlasius Boundary Layer Perturbation for Optimal Wind Rotor” Wind Engineering  Vol 35 No.5. 625-634 


Farthing, S.P.  (2012) Wing’d Pumps” [access 5 Sept 2017]


Farthing, S.P. (2013)  Binary Flutter as an Oscillating Windmill – Scaling &Linear Analysis” Wind Engineering  Vol 37 No.5. 483-499


 Farthing, S.P. (2014)“ The Flutterwing WindPumps: Design, NonLinearities, & Measurements” Wind Engineering  2014 Vol 38 No.2 , 217-231


Farthing, S.P. (2015) “Binary Flutter in Water


Farthing, S.P. (2017)”Exact Flutter for the complex Theodorsen function”manuscript submitted to


Kochin, N.E.,  Kibel I.A., & Roze, N. V.  (1964) Theoretical hydromechanics, 1964 by. Translated from the fifth Russian ed. by D. Boyanovitch. Edited by J. R. M. Radok. Publisher: New York, Interscience Publishers


Young,J.,  Lai,J.C.S., Platzer,M.F. (2014)  A review of progress and challenges in flapping foil power generation Progress in Aerospace Sciences  Volume 67, Pages 1-28