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

S.P. Farthing, Applied Mathematician,Wing’d Pump Associates www.econologica.org spfd@cantab.net 975 Tuam Rd. North Saanich B.C. V8L 5P2 Canada

Journal of Aerospace Engineering 2018 Vol. 31, Issue 3

Abstract-

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 V_d. 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 V_d for the optimum tip speed ratio or reduced frequency k, its inertial acceleration w²R multiplies the mass loading for a parasitic inertial loading DV_d²/s times the design useful live loading w. So independent of fluid density there is a scaling maximum for DV_d ²/s. This DV_d ²/s is just the scale of inertial self-loading Dw²R² to the limit stress s. It indicates a serious limit on V_d 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 ³ to absorb the best fraction of the cubic windpower but then inertial reactions on the articulation mechanism would increase as V ⁶_!

Instead at high V_d, 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.

3. FLUTTER MODEL & SYMBOLIC DETERMINANTAL SOLUTION

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 k_n =w_nc/V where w_n 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

(–k²A+ikB+C+ k_n² 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]=A₁₁B₂₂+A₂₂B₁₁-A₁₂B₂₁-A₂₁B₁₂ then |A|=½[A,A] (2)

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

k⁴|A| -ik³ [A,B] - k²(|B|+[A,C]+[A,E] k_n²) + ik ( [B,C] +[B,E] k_n²) + ( |C| + [C,E] k_n² + |E| k_n⁴)=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 E₂₂ in heave. Equating imaginary odd power parts

k²[A,B]= k_n²[B,E] (4)

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

{|A| [B,E] - [a,E][A,B]} k⁴ + { [C,E][A,B] - (|B|+[A,C]) [B,E]}k² = 0 (5)

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

k² = [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 E₂₂ factor so k the conventional positive root is independent of this heave spring

k² = [A,B]C₁₁ - ( |B|+[A,C]) B₁₁ / [A,B]a₁₁ -B₁₁|A| (6)

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

4. EVALUATION OF THE CROSS-DETERMINANTS

(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= ¼prc² 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 mc² /32, a factor of 4 less the real inertia of a solid (ice) cylinder’s mc² /8 . So the virtual inertia about ¼ chord is 3mc² /32, just half of a solid cylinder. These ¼m and 3mc² /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 mjc². The traditional mass ratio is p-1. Then |A|= pj–x² is p /mc² times the minimum inertia about the center of mass > mc²/32. So |A|= pj–x²>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 h”m Dp may be matched solely by an increased spring force hDS with no other terms changed. But as above k² is unchanged by negating just the DE₂₂ proportional to DS alone so that leaves k² 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 k_n too, the determinants are given for general p which indeed cancels out of the (19) for k².

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”+g‘V +g”xc } +L =p₀mh” (7)

-g‘V 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 p₀mh” where if p₀>0 the net fluid force is stiffening in heave. Then

S=mw_n²=(p- p₀-1)mw² or k²_n=(p-p₀-1)k² (8)

The pitch moment balance per unit span about the pitch axis is m {g”j c² +g‘Vcq -h”xc} = -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 =4mV²j/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) =K₁(ik/2)/{K₀(ik/2)+K₁(ik/2)} =F-iG

in terms of modified Bessel functions of the second kind with F =½ at 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 e^{iwt )} H and likewise

g=(Implied Imaginary Im part of e^iwt ) G where G=g₀e^iq where g₀ is the amplitude of g and q its phase lead over H . Theodorsen thus extended L=(Implied Im part of e^iwt ) 4mV²T {G+ikqG-ikH }/c

Then the roll equation (7) divided by mV²e^{iwt /}c extends to G{k²x-ik(1+4Tq)-4T}+H{ k²_n- pk²+4ikT}=0 (11)

In the very high wind limit of k_n 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= k_n =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 mg‘Vcq+4mVFecqg‘=4m FqycV g‘ with y=e+¼/F, so the pitch moment equation (8) non-dimensionalised by mcV² e^iwt extends to

G{-k²j +4ikqyF+4eF}+H{k²x-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+ k²xH=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]=A₁₁=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 w²=S [B,E] / m[A,B]=Sqy / {pqy-(q+y)x+j} (14)

(4) &(14) are also k²{pqy-(q+y)x+j} =qyk²_n Substituting (8) k²_n=(p-p₀-1)k² cancels p for (p₀+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”{jc²-ec²x}-h”c{x-ep}+cg‘V/2+ ecw_n²h=0 (16)

This exchanges (complex) L and its g, g’ and h’ terms for a cross stiffness ecSh= ecw_n²h. 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 ½mcg‘V in (14) from (9) is explicit damping and so must be circulatory as potential flow is perfect (without dissipation). In (9)

mg‘Vc(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 mcV²e^iwt

G{-k²(j-ex)+½ik}+H{k² (x-ep)+ek_n²} =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|=pj–x², [A,B]=4F{pqy-(q+y)x+j}=4F{(j-qx)(j-yx)+(pj–x²)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], A₁₁,B₁₁,C₁₁ but the 2^nd&3^rd 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

k² (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 , k²=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) p₀» (x-q)/e, and left of Eq. (19)≈0»4e (j-ex)+.5x/F-2qy

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

From (4), (13) and (19) k_n²= k²[A,B]/[B,E]= {4eT(j-ex)+ ½x -2Fqy}[1/qy+ (pj–x²)/ (j-qx)(j-yx) ] / j (20)

giving the flutter V² = Sc²/ mk_n² where S is just a factor and the p influence is just through |A|= pj–x²>0. {}>0 above the general k line {}=0 where k² and V²>0 for j<yx or j>qx.

5. FLUTTER Contours AND MODES

Consider the inertia and imbalance for just the virtual mass concentrated and relocated to the ¾ chord point. Then sddx_n=j_n/x_n=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 p₀=0. As k∞, this mode H→Gq 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 x_g=2Fy=2Fe+ ½ and j_g=2Fy²=x_gy. Similarly from (12) its mode is H=G(y-i/k) or h’_y=gV (for local apparent wind angle j_y=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) p₀=2F-1. As k∞ & F↓½, this 2mF, j_g & x_g fuse into the ¾ chord nexus. The locus of (j_g,xg) or the “gate” all the contours cross, is for e=0 just the straight k =0 line at x=½ between the (F=½) nexus j_n=¼ and F=1 ”midpost.” j_g=⅛. The gate is (x_g-½)j_g=ex²_g, 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+½ x_n= q j_n=qx_n 2(j_n-ex_n)=q p₀=0 j_¾=0 H=G(q-i/k)

point 2mF at gate y=e+¼/F x_g=2Fy j_g=yx_g 2(j_g-ex_g)=y p₀=2F-1 j_y= 0 H=G(y-i/k) (x_g-½)j_g=ex²_g

k² (21)

Recasting (19) in terms of these roots as k² (j-qx)(j-yx) (¼/F -½)= (j-j_n)(x-x_g)-(j-j_g)(x-x_n) (21)

Since the virtual mass is not a point but distributed with the intrinsic virtual inertia of mc⁴/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 (p₀+1)qy+j=(q+y)x (14) p₀=0 when qy-(q+y)x+j=0 as at the nexus and close to its (F=½) line

q²-2qx+j=0 which radiates from the nexus with half the slope of the nexus ray. To the left at p₀>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 p₀<0 the net fluid heave inertia force comes from the negative stiffness of the circulatory lift from positive ¾ chord upwash U_{3/4 at} low phase shifts and k. p₀=-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(⅛-e²) with small e slope about -8e. Comparing these posts at F=1 x_m=2e+½=e+x_n and j_m=½x_m² =(2e+½)(e+¼) <j_n=x_n²= (e+½)² until e²=⅛ or e=.35 or the axis 10% ahead of the wing when j_m = j_n =.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 x_n=½ j_n= ¼ and then turn up before x=.4992 and pass through the x_g=½ gate before the midpost at j_m=⅛. Thus there is a micro zone of recrossing and so double k below the gate.

For F=½, k√2=√(x-½)/│j-½x│ 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>x²/(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=(j-½i/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 r²=│H/G│²=( j²-¼k^-2)/x². 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=½, r²=│H/G│² =
{2j(j-½)+¼ x}/ x(2x-1) which is symmetric about the vertical j=¼ through 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 r² contours would be hyperbolas asymptoting to the k=∞ line and this reflection. For F=1, r²={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 x=½ k=0 line. Figure 5 for F(k) indeed has nearly symmetric contours of interpolated amplitude ratio r whose axis moves from j=.23 to j=¼ as 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, Tan²q =(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=cg₀/h₀=│G/H│=1/r Power/ .5rV ³g₀ (h₀ +cg₀) = (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 j_n=.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 k²with 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)² the pitch-only motion at the gate is overdamped by q²/ 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 >> j_n 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<x_m=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 j_g=ex²_g /(x_g-½) is curved downwards and almost identical to the contour k_m passing through the midpost at x_m= ¾ j_m=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 k_m contour but not identical, because (19) has extra j², 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= k²(e+¼)²/2e which is plotted on Fig 3 for e=⅛ to get k_m =1.03 and for e=¼ where it has a maximum of k_m =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 r² 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 r² can be eliminated to yield this case of (19): k²j²/2F=2ej-qy. Feeding k² 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 mK² would give sufficient increase Dj=n²K², 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, V_d² <<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 j_n and improved to k=.52 at j= 3/8 vs. j_n = .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.

Nomenclature

r Fluid density

s Material Fatigue limit stress

g pitch angle

g₀ 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

w_n 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

h_3/4 heave of the ¾ chord point

H non-dimensional h/c complex extended

i square root of -1

j pitch inertia/mc²

k reduced frequency based on chord wc/V

k_m reduced frequency k >0 of contour passing through midpost

k_n reduced natural frequency w_nc/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

V_d 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

References

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” https://www.youtube.com/watch?v=16kB6p-kcC0

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 to “Blasius Boundary Layer Perturbation for Optimal Wind Rotor” Wind Engineering Vol 35 No.5. 625-634

Farthing, S.P. (2012) “ Wing’d Pumps” https://www.youtube.com/watch?v=6zIj7LCtX0U [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) https://www.youtube.com/watch?v=NDm78DOcEOM “Binary Flutter in Water” https://www.youtube.com/watch?v=NDm78DOcEOM

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