Flowfield Estimation in the Wake of a Pitching and ... AWS
Flowfield Estimation in the Wake of a Pitching and Heaving Airfoil Brian T. Hinson and Kristi A. Morgansen
Abstract— Biological systems have shown great abilities for navigating dense and uncertain environments with variable winds or currents. These animals have unique sensors that allow them to detect the local flowfield and adjust their motor controls accordingly. One approach to leveraging biological sensing capabilities for engineered systems is to understand the sensing from a theoretical perspective. In this paper, flowfield observability and estimation is addressed in the framework of detecting vortices in the flow around a flapping airfoil. An unsteady potential flow model is developed, and the flowfield is measured using velocity sensors placed on the surface of the airfoil. The observability characteristics are identified using two methods, which reveal that for a set of surface velocity sensors, nonzero pitch and heave controls are required to achieve full state observability of airfoil velocities, angle-of-attack, and vortex locations and strengths. These interesting results indicate that wing flapping increases environment observability for animals with wing-mounted velocity sensors. The observability results are demonstrated in simulation with an unscented Kalman filter to estimate a discrete number of vortices in the wake of the flapping foil. Although the filter assumes only a small number of vortices to be estimated, the UKF is able to capture the structure of the wake vorticity.
I. I NTRODUCTION Biological systems have a keen ability to navigate dense environments in the presence of temporally and spatially varying flowfields. The hypothesis has been proposed that animals use vortices in the flow to their advantage for propulsion [1], [2], e.g. a hovering hummingbird or a fish flapping its caudal fin. An underlying assumption in this hypothesis is that animals are assumed to be estimating or predicting the local flowfield magnitude, direction, and structure using limited on-board processing. The ability to observe the flowfield dynamics could explain how animals are able to locate obstacles, catch prey, and avoid predators. Obviously, vision can help tremendously with some of these tasks, but research has shown that even blind fish can school due to a set of sensors called the lateral line [3]. These sensors are essentially a set of “hairs” of varying lengths that vibrate in response to disturbances in the surrounding flowfield. Flying animals also have sensors on their wings to observe the local flowfield. Recent experiments with bats has shown that stiff, microscopic hairs on their wings are used to sense vortices at the onset of stall [4]. Bats without these sensory hairs were generally less maneuverable, exhibiting B. T. Hinson and K. A. Morgansen are with the Department of Aeronautics and Astronautics, University of Washington, Box 352400, Seattle, WA, 98195-2400. Contact {bthinson,morgansen}@aa.washington.edu. This material is based upon work supported in part by a National Science Foundation Graduate Research Fellowship under Grant No. DGE-0818124 and in part by ONR grant N000141010952.
faster flight speeds and wider turns. Ruffling of feathers on birds may be functioning similarly during landing as indicators of vortices or flow separation. Clearly, animals are able to sense their environment in a way that engineered systems are far from accomplishing, but how can we leverage these capabilities to enhance our engineering designs? Recent advances in MEMS technology has afforded the design and manufacture of lateral line type sensors that are small enough to be used in practice on engineered systems [5], [6]. However, effective use of this developing technology lags well behind the perception capabilities of biological systems. A promising first step to improve engineering designs is to understand how the flowfield can be estimated from a theoretical perspective, using the inspiration of biological sensors as a guide. Previous work in theoretical vortex flowfield observability has looked at applying dynamic state estimators to a simplified flowfield model. Ide and Ghil [7], [8] developed an extended Kalman filter to estimate the location and strength of both point- and Rankine-type vortices in a two-dimensional flow. The observer development considered both Lagrangian and Eulerian measurements, concluding that Lagrangian measurements are preferred to Eulerian measurements due to loss of tracking when sensors are near the singularity of a point vortex. Tadmor [9] suggested an extended Kalman filter to estimate the dynamics of a rotating vortex pair, using fluid velocity measurements at a fixed point in space. The estimated state was then used in feedback to control the flow. Krener [10] explored the observability characteristics of a system of two point vortices using Lagrangian and Eulerian measurements, showing that the vortex locations and strengths are observable using either one Lagrangian measurement or two Eulerian measurements. Li and Saimek [11] estimated the strength of vortices in the flow around a pitching and heaving airfoil using surface-mounted pressure sensors and an extended Kalman filter, assuming that the vortex locations were known a priori. However, their model did not consider shed vorticity, as is done here. The contributions of this paper are two-fold—first, to analyze the observability characteristics of r point vortices in a flow around an airfoil, and second, to design an observer that is capable of estimating vortex structures in the wake of an oscillating airfoil. The flow dynamics of r point vortices around an airfoil are briefly developed in Section II. In Section III, the observability of vortices is analyzed using two methods. The development of an unscented Kalman filter to estimate vortices in the wake of the foil is presented in Section IV, and simulation results are provided in Section V. Ongoing and future work is discussed in Section VI.
Fig. 1. Geometry of the flowfield estimation problem with p velocity sensors located on the cylinder surface to measure the strength and location of r vortices in the flow.
Fig. 2. Example Joukowski transformation of circle ζ-domain to airfoil z-domain with ζ0 = 0.15 + i0.05.
z = ζ + ζ0 + II. F LOWFIELD M ODEL The flowfield model that will be used here is based on conformal mapping of a circular cylinder to a Joukowski airfoil. Unsteady vortex shedding is implemented through the shedding of discrete vortices at every time step to simultaneously satisfy the Kutta condition and Kelvin’s circulation theorem. To this end, consider a circular cylinder of radius a ∈ R that encompasses the origin of the complex ζ-plane and whose center is at some location, ζ0 = ξ0 + iη0 ∈ C, where ξ0 , η0 ∈ R, as pictured in Fig. 1. The cylinder has p surface-mounted sensors that measure velocity tangent to the cylinder’s surface, which are chosen to emulate the small velocity-sensing hairs on bats or on the lateral line of fish. The estimation problem is to determine the location, ζvj = ξvj + iηvj ∈ C, and strength, κvj ∈ R, of r vortices in the flow using the set of p velocity sensors. The flowfield is assumed to be inviscid, incompressible, and irrotational, which allows it to be modeled using potential flow theory. Under these assumptions, the flowfield can be completely described at each instant by the complex potential, w(ζ, t) = φ(ζ, t) + iψ(ζ, t) ∈ C, where φ(ζ, t), ψ(ζ, t) ∈ R, which is assembled through superposition of elementary flows. The complex potential around a circular cylinder centered at the origin is first derived in the ζ-plane through use of the Milne-Thompson circle theorem [12], and is given by a2 w(ζ, t) = V (t)eiα(t) ζ + V (t)e−iα(t) + iκ0 log(ζ) ζ r X iκj log(ζ) + log(ζ − ζj (t)) +
(1)
j=1
a2 . − log ζ − ¯ ζj (t)
This √ equation defines the flowfield moving with velocity V = u2 + v 2 around a cylinder at an angle α with the ξ-axis and a circulation of strength 2πκ0 at the cylinder’s origin to satisfy the Kutta condition. The corresponding dynamics in the airfoil z-plane are obtained through thepJoukowski transformation, given by (2), where l = −ξ0 + a2 − η02 is the distance between the circle’s center and the intersection with the negative ξ-axis [12]:
l2 . ζ + ζ0
(2)
The circle’s center, ζ0 , and radius, a, define the parameters of the resulting Joukowski foil. The ξ-coordinate of the circle’s center determines the thickness of the profile (ξ0 = 0 yields zero thickness), and the η-coordinate of the circle’s center determines the camber of the profile (η0 = 0 yields a symmetric profile). An example transformation is shown in Fig. 2. In order to facilitate airfoil flapping, vortex shedding is modeled through point vortices shed at discrete time intervals. New vortices are introduced into the flow at a pre-determined location of ζr = −1.5a to simulate vortex shedding from the airfoil trailing edge, similar to [13] (see [14] for a thorough review of discrete vortex shedding methods). The strength of the newly shed vortex, κr , and the airfoil bound circulation, κ0 , are calculated to simultaneously satisfy the Kutta condition (finite velocity at the airfoil trailing edge) and Kelvin’s circulation theorem (total circulation is constant). The Kutta condition is satisfied by forcing a stagnation point at ζ = −a in the ζ-domain: iκ0 dw (−a, t) = V eiα(t) − V e−iα(t) − dζ a r X 1 1 iκj − + + a −a − ζj (t) j=1 # 1 − = 0. 2 −a − ζ¯ja(t)
(3)
Kelvin’s circulation theorem states that the circulation in a potential flow is constant in time. Applying this theorem X to the current flowfield yields the relation that κ0 = − κj . Therefore, the airfoil bound circulation is equal and opposite the circulation in the wake. Substituting this expression for the bound circulation into (3) and solving for the shed vortex strength yields " # r−1 X 1 1 −2V sin α − − 2 −a − ζj −a − aζ¯j j=1 κr = . (4) 1 1 a2 −a−ζr − −a− ζ¯
r
Therefore, the airfoil initially starts with zero circulation and no vortices in the wake, and the circulation changes
as vortices are introduced into the wake. Each vortex in the wake convects according to Routh’s rule [15], which gives the velocity in the z-domain due to all other elements of the flow. If the complex potential due to all sources except the k th vortex is denoted Wk (z, t) ∈ C and the corresponding potential in the ζ-domain is denoted wk (ζ, t), then Routh’s rule states that the position of the k th vortex, zvk = xvk + iyvk ∈ C, where xvk , yvk ∈ R, convects according to the real and imaginary parts of dWk (ζk + ζ0 )2 dwk = dz zk dζ ζk (ζk + ζ0 )2 − l2 −
iκk l2 (ζk + ζ0 )
((ζk + ζ0 )2 − l2 )
(5)
2.
Rewriting the flowfield dynamics in the airfoil body reference frame in control-affine form, one obtains x˙ = f0 (x) +
2 X
fj (x)uj
j=1
0 0 0 x˙ v1 y˙ v1 0 .. .
= x˙ v r y˙ vr 0
0 1
0 0 1 yv1 −xv1 0 .. .
0 0 u1 + u2 0 .. . 0 yvr −xvr 0 0 0 dζ yi = (−u(ζsi ) sin θsi − v(ζsi ) cos θsi ) (ζs ) , dz +
x2 x21 +x22
(6)
where the measurement equation yi is evaluated for i ∈ {1, 2, . . . , p}. The state vector is x ∈ R3+3r = T u v α xv1 yv1 κv1 . . . xvr yvr κvr , and the available controls are the heave acceleration, u1 = v, ˙ ˙ Note that the state vector size and pitch rate, u2 = θ. increases by three every time a new vortex is shed into the flow. The output function, y = h(x) ∈ Rp consists of the tangential velocity magnitude and direction (positive counterclockwise) measured on the surface of the airfoil at each sensor location. In general, the functions x˙ vj and y˙ vj are significantly complex in their structure, however, the linearity of potential flow does provide some relief. Superposition of each elementary flow yields the overall flow dynamics. Thus, addition of a new vortex into the flow (such as those shed from the trailing edge) is additive to the dynamics. Another benefit of the linearity of potential flow arises in the estimation problem. Rather than tracking every vortex in the flow, one can simply estimate the mean vorticity and its location within the wake without losing much information, as will be discussed in Section IV.
III. O BSERVABILITY A NALYSIS Before an observer can be constructed to estimate the local flowfield around the oscillating airfoil, the observability characteristics of system (6) must be determined. The basic idea behind observability is that the initial condition, x(t0 ) = x0 can be determined from the time history of the output, y(t), and the time history of the controls, u(t), for t ∈ [t0 , tf ]. More than one approach can be used to check the observability of the control affine dynamics (6), and two different approaches are used here to thoroughly analyze the system observability. First, the observability problem is addressed from a differential geometry perspective in Section III-A to determine what (if any) controls are useful in improving system observability. Then, the system is linearized about a trajectory in Section III-B to calculate the condition number of the estimation problem. A. Nonlinear Observability from Geometric Methods Observability of the control affine system given by (6) can be analyzed using a differential geometry approach. A brief review of the relevant key elements of nonlinear observability based on differential geometry will be presented here. The reader is referred to [16] for a more complete discussion. First some useful nomenclature is defined for the following analysis. The Lie derivative of the output function, h(x), with respect to a vector field, fi (x) ∈ Rn , where x ∈ Rn , is defined as ∂h fi . (7) Lfi h = ∂x Repeated Lie derivatives are subsequently computed by ∂ h)fi , and mixed Lie derivatives are com(Lfk−1 Lkfi h = ∂x i ∂ puted by Lfi fj h = Lfi Lfj h = ∂x (Lfj h)fi . Now define the observability Lie algebra, O, which is the Lie algebra of the output function, h, the drift vector field, f0 , and the control vector fields, fi (in the work here, i ∈ {1, 2}): O = span{LX1 LX2 . . . LXk h} k = 1, 2, . . .
(8)
where Xi ∈ {f0 , f1 , . . . , fm }, for i ∈ {1, 2, . . . , k}. Note that the terms in O include Lie derivatives with respect to the drift vector field and with respect to the control vector fields. The effect of the control vector fields in the observability test for nonlinear systems is counter-intuitive to what we know from linear systems, where studying the uncontrolled system is sufficient to analyze the observability. For nonlinear systems, however, the actuation and sensing can be coupled, which means that they must be analyzed simultaneously. The Jacobian of the observability Lie algebra can be used to determine local invertability of the relations between the measurements and the states and thus is an indicator of observability in the system. The Jacobian is given by dO = ∂/∂x(O). A condition for local observability at some x0 is that dO have rank n at x0 . Unfortunately, for many systems, and in particular those that are relatively complex or those that have high-order dynamics, the rank of dO may be difficult or impossible to compute analytically in reasonable time. In particular, for the present system, the analytical rank calculation is not practical, even for the case of one vortex
in the flow. To address this issue, a numerical approach is used via a Monte Carlo rank check over a large range of state values. For the case of one vortex in the flowfield and p = 2 sensors on the surface of the airfoil, the Monte Carlo rank check is performed for 100 uniformly distributed random state values and randomly selected sensor locations (excluding the trailing edge). The choice of random sensor locations is simply to avoid selecting sensor locations that may skew the results. Some sensor locations may in fact be better than others, but this is the subject of ongoing work. The numerical observability rank check reveals that h ∂ Lf0 h (9) dO = ∂x Lf1 h Lf2 h
has a rank of six, which is equal to the dimension of the system. This result indicates that the vortex location and strength, as well as the airfoil velocities and angle-of-attack are observable using the array of two sensors. One can consider higher order Lie derivatives with respect to just the drift rather than with respect to the controls. However, these derivatives, while nonzero, do not increase the rank of the Jacobian. Therefore, the Lie derivatives with respect to the controls are necessary for full state observability, and active control is required to achieve this result. The required controls correspond to motion in the f1 and f2 directions, which means that nonzero heave acceleration and pitch rate inputs must be present to have full state observability. The states requiring actuation to be observed can be isolated by incrementally T direct measurements of each state T adding ′ xi ) and noting how the rank of dO (i.e h = h changes. The results show that vortex location and strength are uniquely determined without any control actuation, however, airfoil angle-of-attack and velocity cannot be uniquely determined from velocity measurements on the surface of the airfoil without control input. This result is interesting because it gives insight into how biological systems are able to sense the surrounding environment using sensory hairs on their wings. To extend the above results from two sensors and a single vortex to a more general setting with arbitrary numbers of sensors and vortices, a Monte Carlo observability rank check was performed for a range of number of sensors and vortices in the flow. In this study, state values were bounded to constrain the system to conditions typical of flapping flight, while remaining in a regime where the inviscid flow assumption was valid. The u velocity was constrained in 0 10 the range of , the v velocity was constrained to −5 5 , angle-of-attack was constrained to −π/4 π/4 , vortex position must remain outside the unit circle in the ζdomain and remain within −5 5 , and the vortex strength was allowed to fall in the range of −1 1 . In all cases, a symmetric Joukowski airfoil was used with ζ0 = 0.1. Each combination of number of sensors and number of vortices was run with 100 uniformly distributed random values within the prescribed state bounds. Results of the Monte Carlo
rank check are shown in Table I, which lists the number of sensors required for local observability. From these results, one can infer that at least 3r − 1 sensors are required for local observability of u and v velocities, angle-of-attack, and location and strength of r vortices. TABLE I M ONTE C ARLO OBSERVABILITY RANK ANALYSIS RESULTS # of vortices 1 2 3 4
# of sensors for local observability 2 5 8 11
B. Linearization About a Trajectory The observability Lie algebra yields useful information about what types of controls are needed to gain full state observability. Unfortunately, it does not tell us how wellconditioned the observation problem is. From the results of the Monte Carlo rank check, we see that the problem of estimating r vortices with two sensors becomes more illconditioned as the number of vortices grows, therefore we need a measure of how observable the system is. An alternate method for checking controlled system observability has been developed by Bastista, et. al. [17], which may yield additional information about the system that the observability Lie algebra does not. This method is based upon the linearization of (6) about a trajectory to obtain a linear, timevarying system, ˜˙ (t) = A(t, x0 (t), u0 (t))˜ x x(t) + B(t)˜ u(t) ˜ (t) = C(t, x0 (t))˜ y x(t),
(10)
where x0 (t) and u0 (t) are the nominal trajectory and corresponding controls about which the linearization is performed, ˜ (t) ∈ Rn = x(t)−x0 (t) and u ˜ (t) ∈ Rm = u(t)−u0 (t) are x ˜ (t) ∈ Rp is the deviations from the nominal trajectory, and y the measurement deviation. The observability of linear timevarying systems can be determined by calculating the rank of the observability gramian, P (t0 , tf ) =
Ztf
ΦT (t, t0 )C T (t, t0 )C(t, t0 )Φ(t, t0 )dt, (11)
t0
where Φ(t, t0 ) ∈ Rn×n is the state transition matrix defined ˙ t0 ) = A(t, x0 (t), u0 (t))Φ(t, t0 ). System (10) is by Φ(t, globally observable and system (6) is locally observable about (x0 (t), u0 (t)) if P (t, t0 ) has rank n. The controls are not explicitly present in the computation of P (t, t0 ), as would be expected for a linear system, however, the controls are implicitly present through the selection of the nominal controls, u0 (t), in the linearization. Thus the sensing and actuation are still coupled in this method. Given the results of the Lie algebra observability analysis above for the vortex system, we know that actuation in the f1 and f2 directions is required to obtain full state observability; therefore, a nominal trajectory must be chosen with nonzero u01 (t) and u02 (t). Note that “good” choices of controls is still an open question with this analysis. Many options are
possible (infinitely many actually) that provide the necessary excitation of the sensors, but some will be better than others. For highly nonlinear or complicated systems, however, the state transition matrix may be difficult to compute for a given nominal trajectory, and it cannot be computed in most cases for a generic trajectory. Prior to the developments of Batista, et. al. [17], Krener and Ide [18] pursued a similar route for uncontrolled nonlinear systems. They recognized the difficulty in computing either the observability Lie algebra or the observability gramian for complex systems and introduced a numerical method for computing the observability gramian using an empirical approach. The empirical local observability gramian, as they label it, is given by (∆y(1) )T t f Z (∆y(2) )T 1 (1) P˜ = 2 . . . ∆y(n) dt, (12) ∆y .. 4ǫ . t0
Fig. 3. Local estimation condition number as a function of number of sensors, p, and number of vortices, r.
number levels off as the number of sensors becomes large, indicating that the addition of more sensors does not provide much more information.
(∆y(n) )T
where ∆y(i) = y+i (t) − y−i (t). The output, y±i (t) = h(x±i ), is computed by perturbing the initial condition of the nominal trajectory, x±i (0) = x0 (0)±ǫei , where ei is the single-entry vector with one in the ith row. The simulation is therefore run 2n + 1 times — once to compute the chosen nominal trajectory, x0 (t), and two times for each state direction to compute the plus and minus perturbed trajectories, x±i (t). The empirical local observability gramian can be shown to converge to the observability gramian in the limit as ǫ approaches zero. In practice, a small value of ǫ on the order of 10−3 is sufficient for calculating the empirical local observability gramian. Clearly, the main benefit of this method is the ability to calculate observability for any system with a working simulation, regardless of complexity. A secondary benefit that is useful in the scenario here is the fact that if the states are properly scaled, then the condition number of the observability gramian, denoted the local estimation condition number, is a direct measure of the well-posedness of the estimation problem. Additionally, the inverse of the minimum singular value of P˜ gives a measure of how close the system is to singular, which Krener and Ide denote the local unobservability index. The local estimation condition number is computed for system (6) linearized about a nominal trajectory to determine the well-posedness of the estimation problem as a function of the number of vortices and the number of sensors. The nominal trajectory is chosen to emulate an airfoil flapping sequence, with u01 (t) = 4π sin (2πt − π/2) and u02 (t) = π 2 /3 sin (2πt + π/2). The unsteady potential flow simulation was run for one period of oscillation, and the resulting variation of average local estimation condition number as a function of number of sensors for one, two, three, and four vortices in the flow is shown in Fig. 3. Notice that, in general, the local estimation condition number increases with the number of vortices. The condition number is largest with two sensors for each number of vortices and decreases with increasing number of sensors. However, the condition
IV. O BSERVER S YNTHESIS Due to the highly nonlinear dynamics of (6), an unscented Kalman filter (UKF) was chosen as the best option for estimating vortex position and strength in the wake of the oscillating airfoil. The UKF algorithm used here is presented in [19] in discrete-time form and is applied to a discrete time version of system (6) in MATLAB. The augmented estimation T T ˆk w ˆ kT v ˆ kT , consists of the ˆ ak ∈ R2ˆn+p = x vector, x ˆ k ∈ Rnˆ , the estimated process noise, estimated state, x n ˆ ˆ k ∈ R , and the estimated measurement noise, v ˆ k ∈ Rp at w discrete time step k. Here, n ˆ = 3+3ˆ r, where rˆ is the number of vortices to estimate in the flow. The UKF is based upon an unscented nonlinear transformation, where the propagation of noise through the system is estimated both with the state and by the use of sigma points on either side of the state estimate. The matrix of state and sigma points is a p p ˆ ak 1T − Pka . ˆk x ˆ ak 1T + Pka x (13) χak = x
The augmented covariance matrix Pka ∈ R2ˆn+p×2ˆn+p consists of the state and process and measurement noise variance and covariances. At each time step, the UKF propagates 4ˆ n + 2p + 1 sigma points through the state equations to estimate the state and noise propagation through time. Although the oscillating airfoil sheds a vortex at each time step of the simulation, tracking every vortex in the flow is not practical or necessary. The overall coherent structure of vortices can be captured by estimating a small number of vortices. The method for determining rˆ as a function of time is based upon thepestimated steady-state airfoil bound circulation, κ ˆ0 = ˆ21 + x ˆ22 sin x ˆ3 . If κ ˆ 0 changes by a specified amount, 2a x then a new vortex is added to the estimation state vector, and rˆ is incremented by one. The choice of threshold for change in κ ˆ 0 affects the number of vortices that will be estimated. If the threshold is set very low (e.g. equivalent to one degree change in angle-of-attack), then new vortex estimates will frequently be added to the flow. If the threshold is set very
high (e.g. equivalent to thirty degree change in angle-ofattack), then large intervals of time will pass between vortex estimates. The threshold should be chosen to allow at least two new vortex estimates per cycle in order to capture the wake structure. Although the steady-state airfoil circulation will always over-estimate the true circulation in the unsteady potential flow, a change in κ ˆ 0 directly correlates to a change in the true κ0 . The initial estimate for shed vortex strength will be slightly scaled, but the UKF has been shown here in simulation to quickly converge to the true value of local vorticity. The initial position for the new estimated vortex is predetermined, and the strength is set by ζˆj (0) = −1.5a κ ˆ j (0) = κ ˆ 0 (kj−1 ) − κ ˆ0 (kj ).
Evolution of the wake of an oscillating airfoil, t/T = 4.
(14)
This vortex initialization method is essentially enforcing Kelvin’s circulation theorem by requiring another vortex estimate to be made when the change in circulation is significant, as defined by the κ ˆ 0 threshold. Although this method leads to fewer discrete vortex estimates than the true number of vortices, the estimates will converge to the average local vorticity centroid, capturing the overall wake structure while saving on computation. At the addition of a new vortex estimate into the flow, the estimation state vector and covariance must be augmented. TThe state vector increases its T dimension by three to ˆ k−1 x ˆk = x ˆj (0) yˆj (0) κ ˆ j (0) , where x x ˆj (0)+iˆ yj (0) is the Joukowski transformation of the point ζˆj (0). The augmented covariance matrix is modified by adding six rows and columns into the matrix (three each for the state and process noise covariances) with predetermined variances and zero covariances. The UKF then proceeds with the estimation algorithm using the new state vector and covariance matrix. Note that the number of sigma points to propagate through the state equations increases by twelve for each vortex estimate added. Thus, the number of discrete vortex estimates must be kept to a minimum to allow for practical computation time and effort. V. S IMULATION R ESULTS The potential flow simulation governed by (6) is integrated in time using an Euler integration routine in MATLAB with a time step of ∆t/T = 0.02, where T = 2π/ω is the period of flapping for a frequency of ω. Each simulation is run with cyclic control inputs u1 (t) = a1 sin (ωt + φ1 ) u2 (t) = a2 sin (ωt + φ2 ).
Fig. 4.
(15)
Simulations were first run to validate the behavior of the unsteady potential flow model. The first simulation was run for a cyclic pitch rate input with a2 = π 2 /3 and φ1 = π/2. The velocity u was chosen to produce a reduced frequency of k = ωa/u = 1. Fig. 4 shows the resulting wake formation behind the flapping airfoil after four flapping cycles. Notice how the small shed vortices interact and roll-up to create the often observed vortex roll-up downstream of the airfoil.
Fig. 5. Comparison of unsteady simulation model and quasi-steady theory for a pitchup maneuver, k = 0.25.
A second simulation was run to demonstrate the effect of unsteady fluid flow on the airfoil bound circulation. Quasisteady theory states that for a symmetric airfoil in the limit of small reduced frequency (k ≪ 1), the bound circulation is a function of the instantaneous angle-of-attack, κ0qs = 2aV sin α. However, for an unsteady model, the circulation will always remain bounded by the quasi-steady assumption due to a phase lag and gain reduction that is a function of the reduced frequency. Fig. 5 shows the response of the unsteady simulation model to a thirty degree pitchup maneuver. The pitchrate input was sinusoidal for t/T ≤ 0.25, followed by zero control input for the remainder of the simulation. Notice how the quasi-steady model directly follows the angle-ofattack trajectory, while the unsteady model asymptotically approached the quasi-steady result in steady state. This simple simulation further validates the simulation model and yields confidence in the observability and estimation results. The UKF was run using the output of the simulation to compute the measurement function at each time step and the system (6) dynamics to propagate the sigma points. The estimation routine was run for a simulation with control inputs in both the u1 and u2 channels, with amplitudes of a1 = 4π and a2 = π 2 /3, frequencies of ω = 2π, and a phase shifts of φ1 = −π/2 and φ2 = π/2. Vortex estimates were added to the flow, as described in Section IV, when a change in steady-state circulation was greater than the equivalent of a ten degree change in angle-of-attack. This threshold for circulation change resulted in approximately four shed vortices per period of oscillation. A minimum of two vortices per period is required to capture the wake structure, and four to six vortices have been sufficient in practice. Fig. 6 shows
tion code. This condition number showed that the estimation problem becomes more ill-conditioned with growing number of vortices and gave an indication of how many sensors to use for a given number of vortices. Furthermore, the numerical study indicated that sensor position may have a direct effect on the observability of the system, which is an interesting path for future work. A UKF was designed to estimate the average vorticity in the wake as it propagates downstream. Overall, the UKF performed well due to the estimation of only a discrete number of vortices, rather than adding another vortex to be estimated at every time step. New vortex estimates were initialized into the flow after sufficient change in the airfoil bound circulation was detected. Although this method of determining how many vortices to estimate is somewhat crude, it proved to be successful in simulation. Future work will include estimating the airfoil bound circulation as another parameter in the UKF. Additional ongoing work is directed toward applying these results to biologically observed data and sensing capabilities. R EFERENCES
Fig. 6.
Vortex estimation simulation results, t/T = 1, 1.5, and 2.
the results of the estimation scheme after 1, 1.5, and 2 cycles, where the small filled circles are the shed vortices, and the large unfilled circles are the estimated vortices, each with a radius that is indicative of the vortex strength. Notice that the estimated vortices tend to group around the strongest concentrations of vorticity in the wake. This behavior can be expected since the UKF is estimating the average vorticity and its centroid in the wake. VI. C ONCLUSIONS An unsteady potential flow model was developed to represent the dynamics of N vortices flowing past a heaving and pitching airfoil. Observability analysis of the system dynamics was performed using two methods. The first method is based on the observability Lie algebra, which gives a matrix rank condition to check for observability. The observability Lie algebra is beneficial because it explicitly accounts for the coupling between sensing and actuation inherent in many nonlinear systems and gives an indication of what types of control inputs are required to obtain full state observability. The observability Lie algebra rank check showed that control actuation in the u1 and u2 channels is required to achieve full state observability in the scenario considered here. This result is interesting because it states that the act of wing flapping in fact makes the local flowfield more observable to velocity sensors on the wing. A second observability check was performed to gauge how the observability of the system changes with varying numbers of sensors and vortices. A numerical method from Krener and Ide [18] was used to compute the estimation condition number using only simula-
[1] D. Weihs, “Hydromechanics of fish schooling,” Nature, vol. 241, pp. 290–291, Jan. 1973. [2] B. Ahlborn, D. G. Harper, R. W. Blake, D. Ahlborn, and M. Cam, “Fish without footprints,” J. Theor. Biol., vol. 148, pp. 521–533, 1991. [3] T. J. Pitcher, B. L. Partridge, and C. S. Wardle, “A blind fish can school,” Science, vol. 194, no. 4268, pp. 963–965, Nov. 1976. [4] S. Sterbing-D’Angelo, M. Chadha, C. Chiu, B. Falk, W. Xian, J. Barcelo, J. M. Zook, and C. F. Moss, “Bat wing sensors support flight control,” PNAS, vol. 108, no. 27, pp. 11 291–11 296, July 2011. [5] I. Kao, A. Kumar, and J. Binder, “Smart MEMS flow sensor: Theoretical analysis and experimental characterization,” IEEE Sensors J., vol. 7, no. 5, pp. 713–722, May 2007. [6] Z. Fan, J. Chen, J. Zou, J. Li, C. Liu, and F. Delcomyn, “Development of artificial lateral-line flow sensors,” in Solid-State Sensor, Actuator and Microsystems Workshop, June 2002. [7] K. Ide and M. Ghil, “Extended Kalman filtering for vortex systems. Part I: Methodology and point vortices,” Dynamics of Atmospheres and Oceans, vol. 27, pp. 301–332, 1997. [8] ——, “Extended Kalman filtering for vortex systems. Part II: Rankine vortices and observing-system design,” Dynamics of Atmospheres and Oceans, vol. 27, pp. 333–350, 1997. [9] G. Tadmor, “Observers and feedback control for a rotating vortex pair,” IEEE Trans. Contr. Syst. Technol., vol. 12, no. 1, pp. 36–51, Jan. 2004. [10] A. J. Krener, “Observability of vortex flows,” in Proc. of the 47th CDC, 2008, pp. 3884–3889. [11] P. Y. Li and S. Saimek, “Modeling and estimation of hydrodynamic potentials,” in Proc.of the 38th CDC, 1999, pp. 3253–3258. [12] L. M. Milne-Thomson, Theoretical Aerodynamics, 4th ed. Dover, 1973. [13] J. Cochran, S. D. Kelly, H. Xiong, and M. Krstic, “Source seeking for a Joukowski foil model of fish locomotion,” in Proc. of the ACC, June 2009, pp. 1788–1793. [14] H. Xiong, “Geometric mechanics, ideal hydrodynamics, and the locomotion of planar shape-changing aquatic vehicles,” Ph.D. dissertation, University of Illinois at Urbana-Champaign, 2007. [15] P. K. Newton, The N-Vortex Problem. Springer, 2001. [16] H. Nijmeijer and A. J. van der Schaft, Nonlinear Dynamical Control Systems. Springer, 1990. [17] P. Batista, C. Silvestre, and P. Oliveira, “Single range aided navigation and source localization: Observability and filter design,” Systems & Control Letters, vol. 60, no. 8, pp. 665–673, Aug. 2011. [18] A. J. Krener and K. Ide, “Measures of unobservability,” in Proc. of the 48th CDC, 2009, pp. 6401–6406. [19] J. L. Crassidis and J. L. Junkins, Optimal Estimation of Dynamic Systems. Chapman & Hall, 2004.
Abstract— Biological systems have shown great abilities for navigating dense and uncertain environments with variable winds or currents. These animals have unique sensors that allow them to detect the local flowfield and adjust their motor controls accordingly. One approach to leveraging biological sensing capabilities for engineered systems is to understand the sensing from a theoretical perspective. In this paper, flowfield observability and estimation is addressed in the framework of detecting vortices in the flow around a flapping airfoil. An unsteady potential flow model is developed, and the flowfield is measured using velocity sensors placed on the surface of the airfoil. The observability characteristics are identified using two methods, which reveal that for a set of surface velocity sensors, nonzero pitch and heave controls are required to achieve full state observability of airfoil velocities, angle-of-attack, and vortex locations and strengths. These interesting results indicate that wing flapping increases environment observability for animals with wing-mounted velocity sensors. The observability results are demonstrated in simulation with an unscented Kalman filter to estimate a discrete number of vortices in the wake of the flapping foil. Although the filter assumes only a small number of vortices to be estimated, the UKF is able to capture the structure of the wake vorticity.
I. I NTRODUCTION Biological systems have a keen ability to navigate dense environments in the presence of temporally and spatially varying flowfields. The hypothesis has been proposed that animals use vortices in the flow to their advantage for propulsion [1], [2], e.g. a hovering hummingbird or a fish flapping its caudal fin. An underlying assumption in this hypothesis is that animals are assumed to be estimating or predicting the local flowfield magnitude, direction, and structure using limited on-board processing. The ability to observe the flowfield dynamics could explain how animals are able to locate obstacles, catch prey, and avoid predators. Obviously, vision can help tremendously with some of these tasks, but research has shown that even blind fish can school due to a set of sensors called the lateral line [3]. These sensors are essentially a set of “hairs” of varying lengths that vibrate in response to disturbances in the surrounding flowfield. Flying animals also have sensors on their wings to observe the local flowfield. Recent experiments with bats has shown that stiff, microscopic hairs on their wings are used to sense vortices at the onset of stall [4]. Bats without these sensory hairs were generally less maneuverable, exhibiting B. T. Hinson and K. A. Morgansen are with the Department of Aeronautics and Astronautics, University of Washington, Box 352400, Seattle, WA, 98195-2400. Contact {bthinson,morgansen}@aa.washington.edu. This material is based upon work supported in part by a National Science Foundation Graduate Research Fellowship under Grant No. DGE-0818124 and in part by ONR grant N000141010952.
faster flight speeds and wider turns. Ruffling of feathers on birds may be functioning similarly during landing as indicators of vortices or flow separation. Clearly, animals are able to sense their environment in a way that engineered systems are far from accomplishing, but how can we leverage these capabilities to enhance our engineering designs? Recent advances in MEMS technology has afforded the design and manufacture of lateral line type sensors that are small enough to be used in practice on engineered systems [5], [6]. However, effective use of this developing technology lags well behind the perception capabilities of biological systems. A promising first step to improve engineering designs is to understand how the flowfield can be estimated from a theoretical perspective, using the inspiration of biological sensors as a guide. Previous work in theoretical vortex flowfield observability has looked at applying dynamic state estimators to a simplified flowfield model. Ide and Ghil [7], [8] developed an extended Kalman filter to estimate the location and strength of both point- and Rankine-type vortices in a two-dimensional flow. The observer development considered both Lagrangian and Eulerian measurements, concluding that Lagrangian measurements are preferred to Eulerian measurements due to loss of tracking when sensors are near the singularity of a point vortex. Tadmor [9] suggested an extended Kalman filter to estimate the dynamics of a rotating vortex pair, using fluid velocity measurements at a fixed point in space. The estimated state was then used in feedback to control the flow. Krener [10] explored the observability characteristics of a system of two point vortices using Lagrangian and Eulerian measurements, showing that the vortex locations and strengths are observable using either one Lagrangian measurement or two Eulerian measurements. Li and Saimek [11] estimated the strength of vortices in the flow around a pitching and heaving airfoil using surface-mounted pressure sensors and an extended Kalman filter, assuming that the vortex locations were known a priori. However, their model did not consider shed vorticity, as is done here. The contributions of this paper are two-fold—first, to analyze the observability characteristics of r point vortices in a flow around an airfoil, and second, to design an observer that is capable of estimating vortex structures in the wake of an oscillating airfoil. The flow dynamics of r point vortices around an airfoil are briefly developed in Section II. In Section III, the observability of vortices is analyzed using two methods. The development of an unscented Kalman filter to estimate vortices in the wake of the foil is presented in Section IV, and simulation results are provided in Section V. Ongoing and future work is discussed in Section VI.
Fig. 1. Geometry of the flowfield estimation problem with p velocity sensors located on the cylinder surface to measure the strength and location of r vortices in the flow.
Fig. 2. Example Joukowski transformation of circle ζ-domain to airfoil z-domain with ζ0 = 0.15 + i0.05.
z = ζ + ζ0 + II. F LOWFIELD M ODEL The flowfield model that will be used here is based on conformal mapping of a circular cylinder to a Joukowski airfoil. Unsteady vortex shedding is implemented through the shedding of discrete vortices at every time step to simultaneously satisfy the Kutta condition and Kelvin’s circulation theorem. To this end, consider a circular cylinder of radius a ∈ R that encompasses the origin of the complex ζ-plane and whose center is at some location, ζ0 = ξ0 + iη0 ∈ C, where ξ0 , η0 ∈ R, as pictured in Fig. 1. The cylinder has p surface-mounted sensors that measure velocity tangent to the cylinder’s surface, which are chosen to emulate the small velocity-sensing hairs on bats or on the lateral line of fish. The estimation problem is to determine the location, ζvj = ξvj + iηvj ∈ C, and strength, κvj ∈ R, of r vortices in the flow using the set of p velocity sensors. The flowfield is assumed to be inviscid, incompressible, and irrotational, which allows it to be modeled using potential flow theory. Under these assumptions, the flowfield can be completely described at each instant by the complex potential, w(ζ, t) = φ(ζ, t) + iψ(ζ, t) ∈ C, where φ(ζ, t), ψ(ζ, t) ∈ R, which is assembled through superposition of elementary flows. The complex potential around a circular cylinder centered at the origin is first derived in the ζ-plane through use of the Milne-Thompson circle theorem [12], and is given by a2 w(ζ, t) = V (t)eiα(t) ζ + V (t)e−iα(t) + iκ0 log(ζ) ζ r X iκj log(ζ) + log(ζ − ζj (t)) +
(1)
j=1
a2 . − log ζ − ¯ ζj (t)
This √ equation defines the flowfield moving with velocity V = u2 + v 2 around a cylinder at an angle α with the ξ-axis and a circulation of strength 2πκ0 at the cylinder’s origin to satisfy the Kutta condition. The corresponding dynamics in the airfoil z-plane are obtained through thepJoukowski transformation, given by (2), where l = −ξ0 + a2 − η02 is the distance between the circle’s center and the intersection with the negative ξ-axis [12]:
l2 . ζ + ζ0
(2)
The circle’s center, ζ0 , and radius, a, define the parameters of the resulting Joukowski foil. The ξ-coordinate of the circle’s center determines the thickness of the profile (ξ0 = 0 yields zero thickness), and the η-coordinate of the circle’s center determines the camber of the profile (η0 = 0 yields a symmetric profile). An example transformation is shown in Fig. 2. In order to facilitate airfoil flapping, vortex shedding is modeled through point vortices shed at discrete time intervals. New vortices are introduced into the flow at a pre-determined location of ζr = −1.5a to simulate vortex shedding from the airfoil trailing edge, similar to [13] (see [14] for a thorough review of discrete vortex shedding methods). The strength of the newly shed vortex, κr , and the airfoil bound circulation, κ0 , are calculated to simultaneously satisfy the Kutta condition (finite velocity at the airfoil trailing edge) and Kelvin’s circulation theorem (total circulation is constant). The Kutta condition is satisfied by forcing a stagnation point at ζ = −a in the ζ-domain: iκ0 dw (−a, t) = V eiα(t) − V e−iα(t) − dζ a r X 1 1 iκj − + + a −a − ζj (t) j=1 # 1 − = 0. 2 −a − ζ¯ja(t)
(3)
Kelvin’s circulation theorem states that the circulation in a potential flow is constant in time. Applying this theorem X to the current flowfield yields the relation that κ0 = − κj . Therefore, the airfoil bound circulation is equal and opposite the circulation in the wake. Substituting this expression for the bound circulation into (3) and solving for the shed vortex strength yields " # r−1 X 1 1 −2V sin α − − 2 −a − ζj −a − aζ¯j j=1 κr = . (4) 1 1 a2 −a−ζr − −a− ζ¯
r
Therefore, the airfoil initially starts with zero circulation and no vortices in the wake, and the circulation changes
as vortices are introduced into the wake. Each vortex in the wake convects according to Routh’s rule [15], which gives the velocity in the z-domain due to all other elements of the flow. If the complex potential due to all sources except the k th vortex is denoted Wk (z, t) ∈ C and the corresponding potential in the ζ-domain is denoted wk (ζ, t), then Routh’s rule states that the position of the k th vortex, zvk = xvk + iyvk ∈ C, where xvk , yvk ∈ R, convects according to the real and imaginary parts of dWk (ζk + ζ0 )2 dwk = dz zk dζ ζk (ζk + ζ0 )2 − l2 −
iκk l2 (ζk + ζ0 )
((ζk + ζ0 )2 − l2 )
(5)
2.
Rewriting the flowfield dynamics in the airfoil body reference frame in control-affine form, one obtains x˙ = f0 (x) +
2 X
fj (x)uj
j=1
0 0 0 x˙ v1 y˙ v1 0 .. .
= x˙ v r y˙ vr 0
0 1
0 0 1 yv1 −xv1 0 .. .
0 0 u1 + u2 0 .. . 0 yvr −xvr 0 0 0 dζ yi = (−u(ζsi ) sin θsi − v(ζsi ) cos θsi ) (ζs ) , dz +
x2 x21 +x22
(6)
where the measurement equation yi is evaluated for i ∈ {1, 2, . . . , p}. The state vector is x ∈ R3+3r = T u v α xv1 yv1 κv1 . . . xvr yvr κvr , and the available controls are the heave acceleration, u1 = v, ˙ ˙ Note that the state vector size and pitch rate, u2 = θ. increases by three every time a new vortex is shed into the flow. The output function, y = h(x) ∈ Rp consists of the tangential velocity magnitude and direction (positive counterclockwise) measured on the surface of the airfoil at each sensor location. In general, the functions x˙ vj and y˙ vj are significantly complex in their structure, however, the linearity of potential flow does provide some relief. Superposition of each elementary flow yields the overall flow dynamics. Thus, addition of a new vortex into the flow (such as those shed from the trailing edge) is additive to the dynamics. Another benefit of the linearity of potential flow arises in the estimation problem. Rather than tracking every vortex in the flow, one can simply estimate the mean vorticity and its location within the wake without losing much information, as will be discussed in Section IV.
III. O BSERVABILITY A NALYSIS Before an observer can be constructed to estimate the local flowfield around the oscillating airfoil, the observability characteristics of system (6) must be determined. The basic idea behind observability is that the initial condition, x(t0 ) = x0 can be determined from the time history of the output, y(t), and the time history of the controls, u(t), for t ∈ [t0 , tf ]. More than one approach can be used to check the observability of the control affine dynamics (6), and two different approaches are used here to thoroughly analyze the system observability. First, the observability problem is addressed from a differential geometry perspective in Section III-A to determine what (if any) controls are useful in improving system observability. Then, the system is linearized about a trajectory in Section III-B to calculate the condition number of the estimation problem. A. Nonlinear Observability from Geometric Methods Observability of the control affine system given by (6) can be analyzed using a differential geometry approach. A brief review of the relevant key elements of nonlinear observability based on differential geometry will be presented here. The reader is referred to [16] for a more complete discussion. First some useful nomenclature is defined for the following analysis. The Lie derivative of the output function, h(x), with respect to a vector field, fi (x) ∈ Rn , where x ∈ Rn , is defined as ∂h fi . (7) Lfi h = ∂x Repeated Lie derivatives are subsequently computed by ∂ h)fi , and mixed Lie derivatives are com(Lfk−1 Lkfi h = ∂x i ∂ puted by Lfi fj h = Lfi Lfj h = ∂x (Lfj h)fi . Now define the observability Lie algebra, O, which is the Lie algebra of the output function, h, the drift vector field, f0 , and the control vector fields, fi (in the work here, i ∈ {1, 2}): O = span{LX1 LX2 . . . LXk h} k = 1, 2, . . .
(8)
where Xi ∈ {f0 , f1 , . . . , fm }, for i ∈ {1, 2, . . . , k}. Note that the terms in O include Lie derivatives with respect to the drift vector field and with respect to the control vector fields. The effect of the control vector fields in the observability test for nonlinear systems is counter-intuitive to what we know from linear systems, where studying the uncontrolled system is sufficient to analyze the observability. For nonlinear systems, however, the actuation and sensing can be coupled, which means that they must be analyzed simultaneously. The Jacobian of the observability Lie algebra can be used to determine local invertability of the relations between the measurements and the states and thus is an indicator of observability in the system. The Jacobian is given by dO = ∂/∂x(O). A condition for local observability at some x0 is that dO have rank n at x0 . Unfortunately, for many systems, and in particular those that are relatively complex or those that have high-order dynamics, the rank of dO may be difficult or impossible to compute analytically in reasonable time. In particular, for the present system, the analytical rank calculation is not practical, even for the case of one vortex
in the flow. To address this issue, a numerical approach is used via a Monte Carlo rank check over a large range of state values. For the case of one vortex in the flowfield and p = 2 sensors on the surface of the airfoil, the Monte Carlo rank check is performed for 100 uniformly distributed random state values and randomly selected sensor locations (excluding the trailing edge). The choice of random sensor locations is simply to avoid selecting sensor locations that may skew the results. Some sensor locations may in fact be better than others, but this is the subject of ongoing work. The numerical observability rank check reveals that h ∂ Lf0 h (9) dO = ∂x Lf1 h Lf2 h
has a rank of six, which is equal to the dimension of the system. This result indicates that the vortex location and strength, as well as the airfoil velocities and angle-of-attack are observable using the array of two sensors. One can consider higher order Lie derivatives with respect to just the drift rather than with respect to the controls. However, these derivatives, while nonzero, do not increase the rank of the Jacobian. Therefore, the Lie derivatives with respect to the controls are necessary for full state observability, and active control is required to achieve this result. The required controls correspond to motion in the f1 and f2 directions, which means that nonzero heave acceleration and pitch rate inputs must be present to have full state observability. The states requiring actuation to be observed can be isolated by incrementally T direct measurements of each state T adding ′ xi ) and noting how the rank of dO (i.e h = h changes. The results show that vortex location and strength are uniquely determined without any control actuation, however, airfoil angle-of-attack and velocity cannot be uniquely determined from velocity measurements on the surface of the airfoil without control input. This result is interesting because it gives insight into how biological systems are able to sense the surrounding environment using sensory hairs on their wings. To extend the above results from two sensors and a single vortex to a more general setting with arbitrary numbers of sensors and vortices, a Monte Carlo observability rank check was performed for a range of number of sensors and vortices in the flow. In this study, state values were bounded to constrain the system to conditions typical of flapping flight, while remaining in a regime where the inviscid flow assumption was valid. The u velocity was constrained in 0 10 the range of , the v velocity was constrained to −5 5 , angle-of-attack was constrained to −π/4 π/4 , vortex position must remain outside the unit circle in the ζdomain and remain within −5 5 , and the vortex strength was allowed to fall in the range of −1 1 . In all cases, a symmetric Joukowski airfoil was used with ζ0 = 0.1. Each combination of number of sensors and number of vortices was run with 100 uniformly distributed random values within the prescribed state bounds. Results of the Monte Carlo
rank check are shown in Table I, which lists the number of sensors required for local observability. From these results, one can infer that at least 3r − 1 sensors are required for local observability of u and v velocities, angle-of-attack, and location and strength of r vortices. TABLE I M ONTE C ARLO OBSERVABILITY RANK ANALYSIS RESULTS # of vortices 1 2 3 4
# of sensors for local observability 2 5 8 11
B. Linearization About a Trajectory The observability Lie algebra yields useful information about what types of controls are needed to gain full state observability. Unfortunately, it does not tell us how wellconditioned the observation problem is. From the results of the Monte Carlo rank check, we see that the problem of estimating r vortices with two sensors becomes more illconditioned as the number of vortices grows, therefore we need a measure of how observable the system is. An alternate method for checking controlled system observability has been developed by Bastista, et. al. [17], which may yield additional information about the system that the observability Lie algebra does not. This method is based upon the linearization of (6) about a trajectory to obtain a linear, timevarying system, ˜˙ (t) = A(t, x0 (t), u0 (t))˜ x x(t) + B(t)˜ u(t) ˜ (t) = C(t, x0 (t))˜ y x(t),
(10)
where x0 (t) and u0 (t) are the nominal trajectory and corresponding controls about which the linearization is performed, ˜ (t) ∈ Rn = x(t)−x0 (t) and u ˜ (t) ∈ Rm = u(t)−u0 (t) are x ˜ (t) ∈ Rp is the deviations from the nominal trajectory, and y the measurement deviation. The observability of linear timevarying systems can be determined by calculating the rank of the observability gramian, P (t0 , tf ) =
Ztf
ΦT (t, t0 )C T (t, t0 )C(t, t0 )Φ(t, t0 )dt, (11)
t0
where Φ(t, t0 ) ∈ Rn×n is the state transition matrix defined ˙ t0 ) = A(t, x0 (t), u0 (t))Φ(t, t0 ). System (10) is by Φ(t, globally observable and system (6) is locally observable about (x0 (t), u0 (t)) if P (t, t0 ) has rank n. The controls are not explicitly present in the computation of P (t, t0 ), as would be expected for a linear system, however, the controls are implicitly present through the selection of the nominal controls, u0 (t), in the linearization. Thus the sensing and actuation are still coupled in this method. Given the results of the Lie algebra observability analysis above for the vortex system, we know that actuation in the f1 and f2 directions is required to obtain full state observability; therefore, a nominal trajectory must be chosen with nonzero u01 (t) and u02 (t). Note that “good” choices of controls is still an open question with this analysis. Many options are
possible (infinitely many actually) that provide the necessary excitation of the sensors, but some will be better than others. For highly nonlinear or complicated systems, however, the state transition matrix may be difficult to compute for a given nominal trajectory, and it cannot be computed in most cases for a generic trajectory. Prior to the developments of Batista, et. al. [17], Krener and Ide [18] pursued a similar route for uncontrolled nonlinear systems. They recognized the difficulty in computing either the observability Lie algebra or the observability gramian for complex systems and introduced a numerical method for computing the observability gramian using an empirical approach. The empirical local observability gramian, as they label it, is given by (∆y(1) )T t f Z (∆y(2) )T 1 (1) P˜ = 2 . . . ∆y(n) dt, (12) ∆y .. 4ǫ . t0
Fig. 3. Local estimation condition number as a function of number of sensors, p, and number of vortices, r.
number levels off as the number of sensors becomes large, indicating that the addition of more sensors does not provide much more information.
(∆y(n) )T
where ∆y(i) = y+i (t) − y−i (t). The output, y±i (t) = h(x±i ), is computed by perturbing the initial condition of the nominal trajectory, x±i (0) = x0 (0)±ǫei , where ei is the single-entry vector with one in the ith row. The simulation is therefore run 2n + 1 times — once to compute the chosen nominal trajectory, x0 (t), and two times for each state direction to compute the plus and minus perturbed trajectories, x±i (t). The empirical local observability gramian can be shown to converge to the observability gramian in the limit as ǫ approaches zero. In practice, a small value of ǫ on the order of 10−3 is sufficient for calculating the empirical local observability gramian. Clearly, the main benefit of this method is the ability to calculate observability for any system with a working simulation, regardless of complexity. A secondary benefit that is useful in the scenario here is the fact that if the states are properly scaled, then the condition number of the observability gramian, denoted the local estimation condition number, is a direct measure of the well-posedness of the estimation problem. Additionally, the inverse of the minimum singular value of P˜ gives a measure of how close the system is to singular, which Krener and Ide denote the local unobservability index. The local estimation condition number is computed for system (6) linearized about a nominal trajectory to determine the well-posedness of the estimation problem as a function of the number of vortices and the number of sensors. The nominal trajectory is chosen to emulate an airfoil flapping sequence, with u01 (t) = 4π sin (2πt − π/2) and u02 (t) = π 2 /3 sin (2πt + π/2). The unsteady potential flow simulation was run for one period of oscillation, and the resulting variation of average local estimation condition number as a function of number of sensors for one, two, three, and four vortices in the flow is shown in Fig. 3. Notice that, in general, the local estimation condition number increases with the number of vortices. The condition number is largest with two sensors for each number of vortices and decreases with increasing number of sensors. However, the condition
IV. O BSERVER S YNTHESIS Due to the highly nonlinear dynamics of (6), an unscented Kalman filter (UKF) was chosen as the best option for estimating vortex position and strength in the wake of the oscillating airfoil. The UKF algorithm used here is presented in [19] in discrete-time form and is applied to a discrete time version of system (6) in MATLAB. The augmented estimation T T ˆk w ˆ kT v ˆ kT , consists of the ˆ ak ∈ R2ˆn+p = x vector, x ˆ k ∈ Rnˆ , the estimated process noise, estimated state, x n ˆ ˆ k ∈ R , and the estimated measurement noise, v ˆ k ∈ Rp at w discrete time step k. Here, n ˆ = 3+3ˆ r, where rˆ is the number of vortices to estimate in the flow. The UKF is based upon an unscented nonlinear transformation, where the propagation of noise through the system is estimated both with the state and by the use of sigma points on either side of the state estimate. The matrix of state and sigma points is a p p ˆ ak 1T − Pka . ˆk x ˆ ak 1T + Pka x (13) χak = x
The augmented covariance matrix Pka ∈ R2ˆn+p×2ˆn+p consists of the state and process and measurement noise variance and covariances. At each time step, the UKF propagates 4ˆ n + 2p + 1 sigma points through the state equations to estimate the state and noise propagation through time. Although the oscillating airfoil sheds a vortex at each time step of the simulation, tracking every vortex in the flow is not practical or necessary. The overall coherent structure of vortices can be captured by estimating a small number of vortices. The method for determining rˆ as a function of time is based upon thepestimated steady-state airfoil bound circulation, κ ˆ0 = ˆ21 + x ˆ22 sin x ˆ3 . If κ ˆ 0 changes by a specified amount, 2a x then a new vortex is added to the estimation state vector, and rˆ is incremented by one. The choice of threshold for change in κ ˆ 0 affects the number of vortices that will be estimated. If the threshold is set very low (e.g. equivalent to one degree change in angle-of-attack), then new vortex estimates will frequently be added to the flow. If the threshold is set very
high (e.g. equivalent to thirty degree change in angle-ofattack), then large intervals of time will pass between vortex estimates. The threshold should be chosen to allow at least two new vortex estimates per cycle in order to capture the wake structure. Although the steady-state airfoil circulation will always over-estimate the true circulation in the unsteady potential flow, a change in κ ˆ 0 directly correlates to a change in the true κ0 . The initial estimate for shed vortex strength will be slightly scaled, but the UKF has been shown here in simulation to quickly converge to the true value of local vorticity. The initial position for the new estimated vortex is predetermined, and the strength is set by ζˆj (0) = −1.5a κ ˆ j (0) = κ ˆ 0 (kj−1 ) − κ ˆ0 (kj ).
Evolution of the wake of an oscillating airfoil, t/T = 4.
(14)
This vortex initialization method is essentially enforcing Kelvin’s circulation theorem by requiring another vortex estimate to be made when the change in circulation is significant, as defined by the κ ˆ 0 threshold. Although this method leads to fewer discrete vortex estimates than the true number of vortices, the estimates will converge to the average local vorticity centroid, capturing the overall wake structure while saving on computation. At the addition of a new vortex estimate into the flow, the estimation state vector and covariance must be augmented. TThe state vector increases its T dimension by three to ˆ k−1 x ˆk = x ˆj (0) yˆj (0) κ ˆ j (0) , where x x ˆj (0)+iˆ yj (0) is the Joukowski transformation of the point ζˆj (0). The augmented covariance matrix is modified by adding six rows and columns into the matrix (three each for the state and process noise covariances) with predetermined variances and zero covariances. The UKF then proceeds with the estimation algorithm using the new state vector and covariance matrix. Note that the number of sigma points to propagate through the state equations increases by twelve for each vortex estimate added. Thus, the number of discrete vortex estimates must be kept to a minimum to allow for practical computation time and effort. V. S IMULATION R ESULTS The potential flow simulation governed by (6) is integrated in time using an Euler integration routine in MATLAB with a time step of ∆t/T = 0.02, where T = 2π/ω is the period of flapping for a frequency of ω. Each simulation is run with cyclic control inputs u1 (t) = a1 sin (ωt + φ1 ) u2 (t) = a2 sin (ωt + φ2 ).
Fig. 4.
(15)
Simulations were first run to validate the behavior of the unsteady potential flow model. The first simulation was run for a cyclic pitch rate input with a2 = π 2 /3 and φ1 = π/2. The velocity u was chosen to produce a reduced frequency of k = ωa/u = 1. Fig. 4 shows the resulting wake formation behind the flapping airfoil after four flapping cycles. Notice how the small shed vortices interact and roll-up to create the often observed vortex roll-up downstream of the airfoil.
Fig. 5. Comparison of unsteady simulation model and quasi-steady theory for a pitchup maneuver, k = 0.25.
A second simulation was run to demonstrate the effect of unsteady fluid flow on the airfoil bound circulation. Quasisteady theory states that for a symmetric airfoil in the limit of small reduced frequency (k ≪ 1), the bound circulation is a function of the instantaneous angle-of-attack, κ0qs = 2aV sin α. However, for an unsteady model, the circulation will always remain bounded by the quasi-steady assumption due to a phase lag and gain reduction that is a function of the reduced frequency. Fig. 5 shows the response of the unsteady simulation model to a thirty degree pitchup maneuver. The pitchrate input was sinusoidal for t/T ≤ 0.25, followed by zero control input for the remainder of the simulation. Notice how the quasi-steady model directly follows the angle-ofattack trajectory, while the unsteady model asymptotically approached the quasi-steady result in steady state. This simple simulation further validates the simulation model and yields confidence in the observability and estimation results. The UKF was run using the output of the simulation to compute the measurement function at each time step and the system (6) dynamics to propagate the sigma points. The estimation routine was run for a simulation with control inputs in both the u1 and u2 channels, with amplitudes of a1 = 4π and a2 = π 2 /3, frequencies of ω = 2π, and a phase shifts of φ1 = −π/2 and φ2 = π/2. Vortex estimates were added to the flow, as described in Section IV, when a change in steady-state circulation was greater than the equivalent of a ten degree change in angle-of-attack. This threshold for circulation change resulted in approximately four shed vortices per period of oscillation. A minimum of two vortices per period is required to capture the wake structure, and four to six vortices have been sufficient in practice. Fig. 6 shows
tion code. This condition number showed that the estimation problem becomes more ill-conditioned with growing number of vortices and gave an indication of how many sensors to use for a given number of vortices. Furthermore, the numerical study indicated that sensor position may have a direct effect on the observability of the system, which is an interesting path for future work. A UKF was designed to estimate the average vorticity in the wake as it propagates downstream. Overall, the UKF performed well due to the estimation of only a discrete number of vortices, rather than adding another vortex to be estimated at every time step. New vortex estimates were initialized into the flow after sufficient change in the airfoil bound circulation was detected. Although this method of determining how many vortices to estimate is somewhat crude, it proved to be successful in simulation. Future work will include estimating the airfoil bound circulation as another parameter in the UKF. Additional ongoing work is directed toward applying these results to biologically observed data and sensing capabilities. R EFERENCES
Fig. 6.
Vortex estimation simulation results, t/T = 1, 1.5, and 2.
the results of the estimation scheme after 1, 1.5, and 2 cycles, where the small filled circles are the shed vortices, and the large unfilled circles are the estimated vortices, each with a radius that is indicative of the vortex strength. Notice that the estimated vortices tend to group around the strongest concentrations of vorticity in the wake. This behavior can be expected since the UKF is estimating the average vorticity and its centroid in the wake. VI. C ONCLUSIONS An unsteady potential flow model was developed to represent the dynamics of N vortices flowing past a heaving and pitching airfoil. Observability analysis of the system dynamics was performed using two methods. The first method is based on the observability Lie algebra, which gives a matrix rank condition to check for observability. The observability Lie algebra is beneficial because it explicitly accounts for the coupling between sensing and actuation inherent in many nonlinear systems and gives an indication of what types of control inputs are required to obtain full state observability. The observability Lie algebra rank check showed that control actuation in the u1 and u2 channels is required to achieve full state observability in the scenario considered here. This result is interesting because it states that the act of wing flapping in fact makes the local flowfield more observable to velocity sensors on the wing. A second observability check was performed to gauge how the observability of the system changes with varying numbers of sensors and vortices. A numerical method from Krener and Ide [18] was used to compute the estimation condition number using only simula-
[1] D. Weihs, “Hydromechanics of fish schooling,” Nature, vol. 241, pp. 290–291, Jan. 1973. [2] B. Ahlborn, D. G. Harper, R. W. Blake, D. Ahlborn, and M. Cam, “Fish without footprints,” J. Theor. Biol., vol. 148, pp. 521–533, 1991. [3] T. J. Pitcher, B. L. Partridge, and C. S. Wardle, “A blind fish can school,” Science, vol. 194, no. 4268, pp. 963–965, Nov. 1976. [4] S. Sterbing-D’Angelo, M. Chadha, C. Chiu, B. Falk, W. Xian, J. Barcelo, J. M. Zook, and C. F. Moss, “Bat wing sensors support flight control,” PNAS, vol. 108, no. 27, pp. 11 291–11 296, July 2011. [5] I. Kao, A. Kumar, and J. Binder, “Smart MEMS flow sensor: Theoretical analysis and experimental characterization,” IEEE Sensors J., vol. 7, no. 5, pp. 713–722, May 2007. [6] Z. Fan, J. Chen, J. Zou, J. Li, C. Liu, and F. Delcomyn, “Development of artificial lateral-line flow sensors,” in Solid-State Sensor, Actuator and Microsystems Workshop, June 2002. [7] K. Ide and M. Ghil, “Extended Kalman filtering for vortex systems. Part I: Methodology and point vortices,” Dynamics of Atmospheres and Oceans, vol. 27, pp. 301–332, 1997. [8] ——, “Extended Kalman filtering for vortex systems. Part II: Rankine vortices and observing-system design,” Dynamics of Atmospheres and Oceans, vol. 27, pp. 333–350, 1997. [9] G. Tadmor, “Observers and feedback control for a rotating vortex pair,” IEEE Trans. Contr. Syst. Technol., vol. 12, no. 1, pp. 36–51, Jan. 2004. [10] A. J. Krener, “Observability of vortex flows,” in Proc. of the 47th CDC, 2008, pp. 3884–3889. [11] P. Y. Li and S. Saimek, “Modeling and estimation of hydrodynamic potentials,” in Proc.of the 38th CDC, 1999, pp. 3253–3258. [12] L. M. Milne-Thomson, Theoretical Aerodynamics, 4th ed. Dover, 1973. [13] J. Cochran, S. D. Kelly, H. Xiong, and M. Krstic, “Source seeking for a Joukowski foil model of fish locomotion,” in Proc. of the ACC, June 2009, pp. 1788–1793. [14] H. Xiong, “Geometric mechanics, ideal hydrodynamics, and the locomotion of planar shape-changing aquatic vehicles,” Ph.D. dissertation, University of Illinois at Urbana-Champaign, 2007. [15] P. K. Newton, The N-Vortex Problem. Springer, 2001. [16] H. Nijmeijer and A. J. van der Schaft, Nonlinear Dynamical Control Systems. Springer, 1990. [17] P. Batista, C. Silvestre, and P. Oliveira, “Single range aided navigation and source localization: Observability and filter design,” Systems & Control Letters, vol. 60, no. 8, pp. 665–673, Aug. 2011. [18] A. J. Krener and K. Ide, “Measures of unobservability,” in Proc. of the 48th CDC, 2009, pp. 6401–6406. [19] J. L. Crassidis and J. L. Junkins, Optimal Estimation of Dynamic Systems. Chapman & Hall, 2004.
