Unsteady forces and flows in low Reynolds number hovering flight:two-dimensional computations vs robotic wing experiments

The demonstrated importance of unsteady effects in insect flight has prompted recent development of better experimental and computational tools to investigate unsteady forces and vortical flows around a flapping wing. In particular, dynamically scaled robots of both the hawkmoth(Ellington et al., 1996) and fruitfly (Dickinson et al.,1999) have been developed to perform controlled measurements of flows and forces. In parallel, researchers have developed direct numerical simulations to probe the detailed vortex dynamics and unsteady forces in flapping flight (Liu et al.,1998; Wang, 2000a,b; Sun and Tang, 2002; Ramamurti and Sandberg,2002).

Given the complexity of modeling fluid flows in three dimensions(Liu et al., 1998; Sun and Tang, 2002; Ramamurti and Sandberg, 2002),it would be desirable to determine if simpler models provide results that are consistent with those generated in experiments. Here, we compare two-dimensional (2D) computations of hovering flight against robotic wing experiments. 2D computations are appealing partly because of their relative simplicity and efficiency. Obviously, 2D computations cannot predict three-dimensional (3D) effects; on the other hand, it is almost impossible to attach the significance of 3D effects without knowing what happens in 2D flow. Therefore, in addition to being relevant to cases where the flow is approximately 2D, as with large wing aspect ratio, when compared with 3D experiments or computations, 2D computations can offer useful insight into the relative significance of 3D effects, as we will discuss at the end of the paper.

Comparing computations and experiments is delicate, partly because it is almost impossible to match any two setups exactly, and partly because it is tempting to present results that compare well, thereby biasing the interpretations. Therefore, it is essential to test the methods in qualitatively different flows generated by different wing kinematics.

In addition to comparing the experimental and computational forces, we also evaluate the relative importance of unsteady effects. These include wing acceleration, both in translation and rotation, and interactions between the wing the existing flow. Most recent work using a robotic fruit fly focused on kinematics based on tethered flight measurements. These kinematics have relatively constant translational velocity in the mid-stroke and large accelerations and sharp rotations at the end of strokes. In these strokes,pronounced peaks appear near the end of each stroke. These peaks were attributed to either wing rotation and wake capture(Dickinson et al., 1999; Birch and Dickinson, 2003), or rotational and translational acceleration(Sun and Tang, 2002; Sane and Dickinson, 2002). The sinusoidally varying strokes studied here offer a set of kinematics where the relative contribution of some of the dynamic effects can be theoretically estimated. For example, we estimate the relative contribution that wing rotation and acceleration make to the quasi-steady forces. We also estimate the wing inertia relative to the fluid force, as well as the non-inertial forces due to wing acceleration relative to the pressure forces associated with vorticity flux. Given that the free flight kinematics of fruit flies appear to be more sinusoidal than those derived from tethered flight(Fry et al., 2003), our results, though using an idealized kinematics, may nevertheless relate to the forces generated by the real flies.

We obtained both force and flow data using a dynamically scaled robotic fly, as described in prior studies(Dickinson and Götz, 1993; Dickinson et al., 1999; Birch and Dickinson, 2003). At the base of one arm was attached a 2D force sensor that measured forces parallel and perpendicular to the wing surface. The force sensor at the wing base was designed to be insensitive to the force moments; thus the force distribution on the wing did not influence overall force magnitude. Lift and drag forces were then calculated from the perpendicular shear forces measured by the sensor. The wing consisted of a 2.25 mm thick piece of Plexiglas, cut to the planform of a Drosophila wing. When attached to the force sensor, total wing length was 0.25 m. The wing and arm apparatus were lowered into a 1 m×1 m×2 m Plexiglas tank filled with 1.8 m³ of mineral oil (density 0.88×10³ kg m^–3;kinematic viscosity=115 cSt). By changing flapping frequency we could operate the robot at Reynolds numbers Re between 50 and 200. The design of the robotic arms permitted the wing to move with three degrees of freedom(measured as stroke amplitude, stroke deviation from horizontal, and wing rotation) and a custom program written in Matlab (Mathworks v.5.3) allowed us to program arbitrary kinematic patterns. To simplify the comparisons with the 2D computational fluid dynamics model, we created a simple back-and-forth wing beat pattern with no stroke plane deviation.

We used digital particle image velocimetry (DPIV) to measure the flow structure in a 841 cm² area centered on the wing. The oil was seeded with air forced through a ceramic water filter stone, creating a dense bubble field. After the larger bubbles rose to the surface, the remaining bubbles, although slightly positively buoyant, were effectively stationary for the duration of each exposure pair. A commercial software package controlling a dual Nd-YAG laser system (Insight v.3.2, TSI Inc., St Paul, MN, USA) created two identically positioned light sheets approximately 2.5 mm thick separated by 2000 μs. These light sheets were parallel to the wing chord and positioned at 0.65R, where R is the wing span, and timed to fire when the wing chord was directly in front of the high-speed video camera placed perpendicular to the laser sheet. We chose 0.65R as our point of measurement because in a prior DPIV study in which the wake was viewed from the rear, 0.65R was the position in which the circulation was the greatest (Birch and Dickinson,2003). We captured one image per stroke from a 29 cm×29 cm area centered on the wing during each of the four strokes. After saving the captured images, the trigger for the laser was advanced and the starting position of the wing was adjusted to line up with the camera at the appropriate time before starting the next trial(Birch and Dickinson, 2003). We repeated this procedure until we had divided each stroke into 10 equally spaced intervals. In this way, we quantified the fluid flow from the perspective of the wing through four downstroke/upstroke cycles, although this paper will only report on the wake dynamics during the fourth stroke. For each image-pair captured, a 2-frame cross-correlation of pixel intensity peaks with 50% overlap of 64 pixel×64 pixel interrogation areas resulted in 900 velocity vectors/image. Vector validation resulted in the removal of only 2 of 9000 vector values; these were filled by interpolation of a mean value from a 3 × 3 nearest neighbor matrix. A custom program written in Matlab(v.5.0) calculated vorticity, ω=Δ×u, from velocity fields smoothed using using a least-squares finite difference scheme.

The computation models a thin wing element of elliptic cross section under the same kinematics as performed in the experiments. The computation of flows around this hovering wing employs a fourth-order finite difference scheme of Navier–Stokes equation in vorticity-stream function formulation(E and Liu, 1996). The scheme is implemented in the elliptic coordinates with appropriate boundary conditions to account for the wing motion (Wang, 2000a,b). See also Russell and Wang(2003) for an alternative method employing Cartitian grids appropriate for multiple wings.

The Navier–Stokes equation is solved in a frame fixed to the wing. In the 2D vorticity stream function formulation, the non-inertial frame introduces only one extra term, the rotational acceleration of the wing. Other non-inertial terms can be expressed as a gradient of a potential function. Thus they can be absorbed into the pressure term. The curl of the gradient of pressure is zero. The body motion is reflected in the far field boundary conditions, and the no-slip boundary condition at the wing is enforced explicitly through the vorticity and stream function boundary conditions. More specifically, on the wing, we set Ψ=c, where c is a constant, to satisfy the no-penetration boundary condition, and ∂Ψ/∂n=0 to satisfy the no-slip boundary condition. At far field, we setΔ×Ψ=–(U₀+r×Ω₀),where U₀ and Ω₀ are the translational and rotational velocity of the wing, respectively, r is the position relative to the wing center, and ω=0. The exact boundary condition onΨ can be recovered by solving the Poisson equation twice(Wang, 1999). For this computation, the far field boundary condition is correct to the dipole order.

The translational motion of the wing is completely specified by two dimensionless parameters, the Reynolds number, Re=U_maxc/ν=πfA₀c/ν,and A₀/c, where U_max is the maximum wing velocity, and c the chord. In the subsequent studies, we fix f but vary A₀/c and study its effect on the flow. For clarity, we will report the value of A₀/c directly instead of Re. A₀/c varies from 2.8 to 4.8, with resulting Re from 75 to 115, appropriate for fruitflies. Other parametersα ₀, β and f are fixed to be π/2, π/4 and 0.25 Hz, respectively.

A main variable of interest in this study is the phase delay between rotation and translation, φ, which was shown to be a sensitive parameter in force generation (Dickinson et al.,1999; Wang,2000b). Three cases, φ=π/4, 0 and –π/4,corresponding to the advanced, symmetrical and delayed rotation(Dickinson et al., 1999), will be studied for each A₀/c.

We normalize the computational and experimental forces by the maxima of their corresponding quasi-steady forces, as described in the next section. Our choice for the normalization is dimensionally the same as the conventional choice, Gρu²_rmsc, but has the advantage of normalizing away features specific to the wing, such as its thickness and geometry. This is because that force dependence of the wing geometry is sometimes relatively simple. For example, the force coefficients of ellipses of different thickness were shown to have almost the same functional dependence on the angle of attack but different magnitude (see fig. 6 in Wang, 2000a). The experimental force coefficients of the robotic fly wing also show little dependence on the wing planform. If the dependence on the geometry in the steady and unsteady forces is similar, then their ratio does not depend sensitively on the geometry of the wing. This would allow us to compare wings of different cross sections and planforms. Although comparing the force coefficients appears to be a natural thing to do, one must be cautious when comparing 2D and 3D force coefficients. The lift coefficient of 1 has different meanings in experiments and computations, unless the sectional lift coefficient in a 3D wing is a constant. Strictly speaking, the numbers should only be compared within the computations or the experiments.

Before discussing the unsteady forces, we first describe the calculation of the quasi-steady forces based on both the translational and rotational velocity. Because the wing operates at a large range of angles of attack, from 0° to 135°, the Kutta–Joukowski lift, which works for attached flow associated with small angles of attack, is clearly inapplicable. Instead,we determine the quasi-steady coefficients empirically, using both the robofly experiments and computation of a steady translating wing at a fixed angle of attack. The lift and drag coefficients, defined with respect to the far field flow, are measured at a time when the forces reach a temporary plateau after the initial transients (see for example, fig. 2 in Dickinson et al., 1999; fig. 5 in Wang, 2000a). Forces at all angles are measured at a fixed time, t=2 in dimensionless time scale.

Note that Equation 18includes both the pressure lift of a translating and rotating wing in the absence of circulation (Magnus force) and the lift due to circulation given by the Kutta's condition. There is no a priori reason as to which of the quasi-steady forms should fit in the unsteady case better, since both use assumptions that are invalid in unsteady cases. Recently, revised quasi-steady models have been proposed to fit these forces(Sane and Dickinson, 2002). For the purpose of this paper, we simply compare L_R, L_T and D_T with the unsteady forces. It turns out that for the prescribed motions here, L_Rdeviates substantially from the unsteady forces, while L_Tapproximates the unsteady forces reasonably well. Therefore, in the subsequent discussions, we will use L_T as an estimate for the quasi-steady forces.

Twelve kinematics are analyzed here: four stroke amplitudes, 60°,80°, 100° and 120°, at three phase delays, φ=π/4, 0 and–π/4. The range of amplitudes corresponds to the end-to-end amplitude:chord ratios of 2.8, 3.6, 4.2 and 4.8 at 70% span. The computational parameters are chosen to match both the Reynolds number Re and the amplitude:chord ratio, A₀/c. In particular, A₀ is estimated to be the projection of the arc length at 70% span onto a 2D plane. The location at 70% span was chosen because the sectional circulation is near maximal there (J. M. Birch, W. Dickson and M. H. Dickinson, manuscript submitted for publication) and the interference from the trailing edge vortex is small. Once we fix the computational units for time and length, i.e. T=1 and c=2.0, the viscosity is also fixed in computational units.

Figs 2, 3, 4 summarize the results for variation in φ. Each figure shows the comparison of experimental and computational force coefficients over the first four cycles. In addition,instantaneous force vectors are shown in superposition with the traveling wing during the second stroke.

The wing motion in these cases differs in the angle of attack at the end of stroke. The angles of attack are π/4, π/2 and 3π/4, respectively, in the advanced (φ=π/4), symmetrical (φ=π/2) and delayed(φ=–π/4) rotation cases, as shown in the force vector plots. The delayed rotation case is unusual from the point of view of operating an airfoil. After each reversal, the wing has angles of attack greater thanπ/2, which leads to a downward lift (see Fig. 4). However, insects or bio-mimetic devices may use such a stroke to reduce the force on one wing, and thus generate a torque to turn. In addition, when the wing is moving at an angle greater than π/2, the flow separates quickly, which is qualitatively different from the other two cases. Thus it also provides a good case for testing computations and experiments in different scenarios.

In all cases, the forces quickly settle into an almost periodic state after two strokes. The computational drag follows the experimental drag closely in all three cases. Lift agrees well in the first two cases (Figs 2, 3), but shows a clear phase delay in the case of delayed rotation (Fig. 4). Notice that the shift occurs only after the first stroke. The averaged experimental lift and drag coefficients are (0.93, 1.28), (0.86,1.34) and (0.38, 1.10), for the advanced, symmetrical and delayed rotation,respectively. The averaged computational lift and drag coefficients are (1.10,1.36), (0.82, 1.44) and (0.19, 1.21), respectively, for the corresponding cases. We will return to the presence and absence of the phase shift in lift in these three different cases when we discuss the 3D effects.

The averaged force coefficients depend weakly on stroke amplitude, as shown in Fig. 5. In the case ofφ=π/4, the average experimental lift coefficients are 0.93, 0.99, 0.95 and 0.93 at A₀/c=2.8, 3.6, 4.2 and 4.8,respectively. The corresponding computational lift coefficients are 1.07, 1.0,0.9 and 0.9. The drag coefficients are 1.28, 1.19, 1.12 and 0.93 in experiments, and 1.36, 1.34, 1.24 and 1.16 in computation. This indicates that the total force scales roughly with

The average lift depends sensitively on φ, as emphasized before(Dickinson et al., 1999; Wang, 2000b). For example, in the case of A₀/c=2.8, the averaged values for experimental lift, C_L are 0.93, 0.86 and 0.38, for φvalues of 0.25π, 0 and –0.25π, respectively. The comparable quasi-steady lift coefficients are 0.75, 0.95 and 0.75. This dependence onφ can be understood intuitively, based on two facts. First, the deviation between the unsteady forces and quasi-steady forces occurs mostly after the flip of each stroke. Second, the instantaneous forces in all these cases are typically normal to the wing, as indicated in the force vector plots in Figs 2, 3, 4, and as discussed in Materials and methods. Therefore, the contribution to lift and drag can be correlated with the instantaneous orientation of the wing at the end of each stroke: π/4, π/2 and 3π/2, for φ=π/4, φ=0 andφ=–π/4, respectively. One expects an increase in both lift and drag when φ=π/4, a decrease of lift and increase of drag whenφ=–π/4, and relatively small change in lift, but a large increase in drag for φ=0. These indeed are consistent with both the experimental and computational forces.

As a further comparison between computation and experiments, we show side by side the snapshots of the vorticity field near the wing from experimental DPIV measurements and computational results(Fig. 6). Ten frames are shown for each period. The colors indicate the strength of the vorticity field. In Fig. 6, columns A and C are computed vorticity, and B and D are experimental vorticity in a 2D slice. The simulations appear to capture the major features of vortex dynamics through a complete stroke cycle. Notice that the fluid momentum is directed downward by pairs of vortices, similar to those shown in asymmetric strokes that model dragonfly wing kinematics (Wang, 2002b). In 3D, the pairs of vortices can be cross-sections of a donut-shaped vortex ring. The structure of the downward momentum jet, characterized by the averaged velocity field over a cycle, is examined elsewhere for both the symmetric and asymmetric strokes strokes (Z. J. Wang, manuscript submitted for publication). Also notice that even the kinematics of left and right strokes are identical, but the flow field differs slightly. This can be seen by comparing columns A and C in Fig. 6. The wing positions are mirror images, but the flows deviate slightly from the mirror symmetry. Such a deviation may be inconsequential in terms of average lift, but worth keeping in mind when interpreting the precise time course of forces.

The differences between the unsteady and quasi-steady forces have been analyzed extensively for fruitfly kinematics, based on results from tethered animals, with relatively constant translational velocity in the mid-stroke and large acceleration and sharp rotations at the end of strokes. The unsteady effects were dominant near the wing reversal, where they contribute to the rotational effect and the wake capture(Dickinson et al., 1999). The discrepancy between experimental measures of forces and flows(Birch and Dickinson, 2003) and a CFD model of nearly identical conditions(Sun and Tang, 2002) raises a debate about the physical basis of these unsteady effects. Here we do not attempt to resolve these discrepancies, but probe the presence of unsteady effects in a different set of kinematics. In Fig. 7, we compare the unsteady forces with the steady state forces based on the translational velocity. In the case of advanced rotation, the unsteady effects can contribute an additional 50% to the total lift, and in the case of delayed rotation, they can reduce the total lift by a factor of 2–3. It is also clear that the quasi-steady forces based on translational velocity alone do not predict the time-dependent forces during a sinusoidal flapping.

To gauge the relative importance of different unsteady effects, we examine the difference between unsteady forces and estimates based on the translational velocity, as shown in Fig. 7D–F. Ideally we wish to decompose, if possible, the unsteady force approximately into a sum of separate terms, which can be related to wing acceleration, the coupling between rotation and translation,wing–wake interaction, etc. However, in the absence of quantitative prediction of these effects, we can only offer a plausible decomposition by correlating the force peaks with the time course of translational and rotational velocity.

The coupling between translation and rotation can be modeled by C_rotΩ(t)u(t), a form predicted by classical theory of a translating and pitching motion(Munk, 1925), and tested in robotic fruitfly experiments (Sane and Dickinson, 2002), where C_rot is assumed to be a constant that depends on the center of rotation. Part of C_rotΩ(t)u(t) is the Magnus force caused by the pressure difference due to velocity difference,given by Bernoulli's law, and another part is due to additional circulation caused by the rotational motion to satisfy the Kutta condition(Munk, 1925). In the three kinematic patterns studied here, the peaks (labeled `r' in Fig. 7) associated with rotation, are picked by matching (with a small shift) the force curve to the maximum of Ω(t)u(t). The positions of these measured force peaks vary in the three cases, in accordance with the shift of the peak positions of u(t)Ω(t). The variation occurs at the same time scale as u(t)Ω(t).

Other unsteady effects occur near the wing reversal (labeled `u' in Fig. 7). The position of these effects occur roughly at the same time in all three kinematics. These force peaks do not follow the trace of du(t)/dt, thus do not behave as the classical added mass. These force peaks are likely related to the unsteady growth of vorticity and wake–wake interaction, which do not have simple analytic expressions in general. Regardless of its physical basis, the most substantial contribution of this unsteady effect is on drag(Fig. 7B,C,E,F).

We also note that the peaks alternate in size from stroke to stroke in the experimental lift, most obvious in Fig. 7A,D. One possible explanation is that this asymmetry reflects the mechanical artifact due to gear backlash. However a small degree of asymmetry is also observed in the computational data, e.g. the vorticity field as shown in Fig. 6 and the forces in Figs 2, 3, 4. At Re=0, the left and right strokes, which are mirror images about the vertical axis, would generate forces that have the same symmetry; i.e. the lift from left and right strokes are identical and the drag forces are of equal size but in the opposite direction. Here, the Reynolds number is finite and sufficiently large that the force can depend on the history of flow. One possibility for breaking the symmetry in forces is by the initial condition.

While identification of the above features is useful when dissecting the unsteady force traces, it is also relevant to determining their net contributions. The force directly proportional to the linear acceleration can have sharp peaks, but it has a zero contribution in reciprocal motions. The pressure force of a rotating and translating plate, approximated by u(t)Ω(t), has a non-zero net contribution in the cases where φ↑0, since[sin(2πft)cos(2πft+φ)]∼sinφ. The unsteady vortex force due to wing acceleration has eluded simple analytical expressions, except for power-law start up flow, where both the added mass term and the vortex force are calculated analytically and numerically(Pullin and Wang, 2003). The unsteady forces contribute to both lift and drag, both predicted in theory(Pullin and Wang, 2003) and seen here in Fig. 7.

Fig. 8 illustrates the contribution of F_ω (broken line) to the total force(solid line). In kinematics where there is a fast acceleration at the end of the stroke, the force will have sharp peaks at the end of the stroke, but the net contribution is zero, as discussed before.

Given that flow around real and model insect wings exhibits 3D effects such as spanwise flow, it is tempting to conclude that 2D computations have little to offer. Here we see that the success and failure of a 2D model in capturing the forces in 3D experiments can provide important insights. In both the advanced and symmetrical rotation cases, the 2D forces are very similar to the 3D forces. Why do they agree? In addition, what do we learn from comparisons between the 2D computation and the 3D experiments?

First, a notable difference between the experimental and computational forces is seen in the delayed rotation, where there is a clear phase shift between the computed and measured lift. In the canonical example of flow past a 2D cylinder and a 3D sphere, the forces during von Karman vortex shedding also show a phase shift, which was argued to be a generic feature between 2D and 3D flows (Mittal and Balachandar,1995). In view of this, the absence of the phase shift in the advanced and symmetrical cases are particularly interesting. The difference in flow structure in the three cases may be worth further investigation.

Second, these results are relevant to recent discussions about the role of 3D effects on delayed stall. Insects are known to flap their wings at angles of attack much higher, around 35°(Ellington, 1984), than the stall angle of a conventional airfoil, about 12°. As suggested by Ellington et al. (1996), the pressure gradient from root to tip within the vortex core might drive spanwise flow that stabilizes the leading edge vortex by convecting away the vorticity. The spanwise flow was indeed seen by smoke visualization in the robotic hawkmoth experiment, where Re≈5000(Ellington et al., 1996; Willmott et al., 1997). This proposed mechanism is thought to be analogous to that occurring on delta wing aircraft, in which spanwise flow through the vortex core maintains the stability. But as discussed previously(Ellington et al., 1996), the exact conditions for establishing spanwise flow in leading-edge vortex for rotary wings are not completely understood. For example, a helicopter rotor also experiences a pressure gradient, centrifugal and Coriolis forces, but no large-scale spanwise flow is observed(Harris, 1966). Recent smoke visualization of free-flying butterflies also did not observe substantial spanwise flow, but reported high variability of 3D flow patterns(Srygley and Thomas, 2003). DPIV images of flow field in a robotic fruitfly experiment, where Re≈150, showed no substantial spanwise flow inside the core of leading edge vortex, but instead indicted substantial spanwise flow behind the leading edge vortex, which connects to the tip vortex(Birch and Dickinson, 2001). Strictly speaking, there is no contradiction among these experiments regarding the spanwise flow. It is likely that the spanwise flow within the vortex core occurs only at sufficiently large Reynolds number as in the case of hawkmoth,but not at low Reynolds number, as in the case of fruitflies. The details of the spanwise flow can also depend on the wing shape and kinematics. The high variability of the 3D flow patterns shown by these different experiments,however, makes it difficult to conclude that spanwise flow is crucial for generating sufficient lift by a hovering insect.

An alternative explanation for why the conventional stall does not seem to affect a flapping insect wing relates to the time scale governing the flow separation that leads to stall. For example, a 2D translating wing at an angle of attack of 40° and Re=1000, does not show a drop in time-dependent force until the chord travels for about 4 chords, after which the forces become oscillatory due a von Karman shedding(fig.∼6 in Wang, 2000a). Therefore, in theory there is no need for additional mechanisms to stabilize the leading edge vortex if the wing travels less than about 4 chords. The early data compiled by Weis-Fogh (1973)showed that ratios of stroke-arc to wing-chord of different species during hovering, including bats (Plecotus auritus), birds (hummingbirds),butterflies and moths (Lepidoptera), wasp and bees (Hymenoptera) and flies(Diptera), have values less than 4. The beetles (Coleoptera) have values between 5 and 6. A main difference between a 3D revolving wing and a 2D translating wing, as noted in recent literature, is that a revolving wing does not appear to shed its leading edge vortex after a distance much longer than the stroke length of a typical insect(Dickinson et al., 1999; Usherwood and Ellington,2002). No doubt such a difference would affect the force coefficients observed in 2D and 3D in the steady state. On the other hand, the difference in terms of vortex shedding may not affect the transient values. It is worth re-examining the results of 3D experiments on a flapping wing(fig. 2D in Dickinson et al., 1999), which show that while the 2D steady state lift coefficients underpredict substantially their the 3D counterparts, the 2D transient values follow closely the 3D coefficients, up to an angle of attack of about 72°. The 3D steady force is slightly lower than the unsteady 2D counterpart, due the well-known downwash due to tip vortices. Similarly, in the cases studied here,the chord is moving between 3 to 5 chords, and the leading edge vortex does not appear to separate during each stroke in the cases of advanced and symmetrical rotation, as indicated by the absence of phase shift between the 2D and 3D forces. In these cases, the 2D forces are good approximations of 3D experiments.

Z.J.W. acknowledges support from AFOSR, NSF Early Career Award and the ONR Young Investigator Program. M.H.D. thanks the Packard Foundation and the National Science Foundation (IBN-0217229) for support of this work.