Category UNSTEADY INCOMPRESSIBLE POTENTIAL FLOW

As was mentioned in the introduction to this chapter, the solution of the complete flowfield requires the solution of the outer potential flow and the inner viscous boundary layer (near the solid surface). Typical information that can be obtained by the solution of the attached boundary layer problem is

1. Displacement effects due to the slower velocity inside the boundary layer.

2. Surface skin friction, so that the contribution of friction to the drag force can be estimated.

3. Indications about the tendency of the flow to detach (or separate). If the boundary layer does separate then the boundary layer solution beyond the separation point cannot be calculated (without using more complete viscous flow equations or interaction techniques).

In this section these topics will be discussed very briefly, too, and a more comprehensive description of them can be found in books dealing with viscous flows and boundary layers (e. g., Schlichting,1 6 AGARD CP-291,14 2 etc.).

The conditions that the wake will move with the local streamlines (and carry no loads) were introduced as early as Section 4.7 for thin lifting surfaces and later in Section 9.3 when discussing the wake model for panel methods. From the steady-state flow point of view, the shape of the wake is not known, and the process of finding the proper wake shape (wake rollup) is often denoted as a “slight nonlinearity” in the solution process. Typical remedies to this problem are:

Prescribe wake shape. This is done in Chapters 4, 8, and 11 for the lifting line and lifting surface type of solution (where the wake is placed on the z = 0 plane). A more refined alternative of this option is to prescribe the wake shape based on flow visualizations. This approach is very helpful when analyzing multielement wings where, for example, in the case of a two-element airfoil the wake of the main airfoil is very close to the trailing-edge flap upper surface.

Wake relaxation. This is a process used by several steady-state numerical solutions, and to demonstrate the principle of this method let us follow the approach used in the code VSAERO.92,12 11 The initial wake geometry is specified by the programmer (usually as a planar wake extending backward from the trailing edge) and then several wake grid planes (normal to the free stream) are established, as shown in Fig. 14.1. For the first iteration, the flowfield over the wing and the initial wake shape are calculated using the method described in Section 12.5.

For the second iteration the velocity induced by the wing and wake (n, v, w)h at each of the wake points (formed by the intersection between the

Trailing vc

wake grid planes and the wake lines) is calculated. Next, the wake points are moved with the local induced velocity (see Fig. 14.1) by

(Ax, Ay, Az), = (u, v, w), At

(Some methods, for simplicity, will move the wake in the wake grid plane only. If this is done in the free-stream coordinate system, then the wake grid lies in the x = const, plane and only (Ay, Az)t are required.) In Eq. (14.1) At is an artificial time parameter and its value can be approximated as

(14.2)

where Ax/ is the distance of the wake grid plane from the trailing edge (or between the wake grid planes) and К has values between 0.5 and 5. Once all the wake points have been moved, due to the local induced velocity, the flow is computed with the new wake geometry and the second iteration cycle (or wake relaxation iteration) has been completed.

These wake relaxation iterations can be continued until convergence is obtained or when sufficient wake rollup has been achieved (decision made by the programmer). Since there is always a risk that by too many iterations the wake can reach levels of nonphysical rollup (where the sum of the iteration time steps E At is much larger than AxJQJ) it is recommended to limit the number of wake rollup iterations to less than three. Results of such a procedure (after two wake iterations) are presented in Fig. 14.2. Here the VSAERO9 21211 code was used and the interaction between a close-coupled wing-canard configuration was calculated. Figure 14.3 shows the effect of the

01______ і_______ I_______ 1_______ I———- 1———- 1———-

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

2 у lb

FIGURE 14.3

Effect of canard position on wing’s spanwise loading. From Katz, J., “Evaluation of an Aerodynamic Load Prediction Method on a STOL Fighter Configuration”, AIAA Papier No. 86-0590, 1986. Reprinted with permission. Copyright AIAA.

canard on the wing’s spanwise loading. Note the noticeable effect of the canard and its wake, which induces a downwash at the wing root area and thereby reduces its lift there. The proper placing of the wake in such cases of closely spaced lifting surfaces is critical and estimation of the wake motion is important for the solution.

Time-stepping method. This approach was demonstrated in Chapter 13 (Sections 13.8.2 and 13.10) and in principle is similar to the wake relaxation method, but now the time step is directly related to the motion. (Therefore, the apparent “slight nonlinearity” does not exist.) From the point of view of computations, the number of wake points increases with time and, for example, for N wake lines during К time steps approximately NK/2 wake point velocity computations are required. When using the wake relaxation approach, even for the first iteration, all wake grids are used and therefore NK such velocity calculations are required. Thus, even for steady-state flows, this time-stepping wake rollup method may require less computational effort.

As an example for this wake rollup calculation consider the rollup of a single horseshoe vortex. In this case, the wing bound vortex is modeled by a

FIGURE 14.4

Instability of a pair of trailing vortices, and comparison between calculated and observed vortex formations. Figure 14.4a shows the wake behind the airplane after its passage and Fig. 14.46 depicts the Crow instability which is shown later at a distance of about 80 wing spans. More details about such calculations can be found in Rossow, V. J., Journal of Aircraft, vol. 24, no. 7, 1987, pp. 433-440. Photo from Ref. 14.1. Reprinted with permission of AIAA and Meteorology Research, Inc. Photo originally appeared in Smith, T. B. and Beemer, К. M., “Contrail Studies of Jet Aircraft," MRI Report, April 1959.

single vortex line that sheds two wake line segments of length Q „ At at its tips during each time step. As this shedding process continues, the two long trailing vortices are formed, but because of the instability of these two vortex lines a sinusoidal pattern will develop. This instability was first analyzed by Crow14 1 who presented the photographs appearing in Fig. 14.4 in Ref. 14.1. The numerical solution presented in Fig. 14.4a is obtained by using only one panel in the method described in Section 13.12, and the above instability is also visible in the computation. Calculations such as this one and that of Fig. 13.29 indicate that this approach for calculating the rollup of vortex sheets yields satisfactory results (at least when modeling trailing wakes behind wings).

As a closing remark to most of the wake rollup modeling efforts, we must emphasize that the velocity induced by a vortex point or line is singular (see, e. g., Fig. 3.8a). Therefore, an artificial vortex core (or cut-off distance) must be defined for the purpose of numerical solutions. One possibility is to define the self-induced influence as zero within this radius; however, in some methods a solid-body rotation model is used within this core (which is very similar to Fig. 2.11 with є being the core size).

Toward the end of Chapter 1 (Section 1.8) it is postulated that many flowfields of interest to the low-speed fluid dynamicist lie in the range of high Reynolds number. Consequently, for attached flowfields, the fluid is divided into two regions: (a) the thin inner boundary layer, and (b) the mainly inviscid irrotational outer flow. Chapters 2-13 are entirely devoted to the solution of the inviscid outer flow problem, which indeed is capable of estimating the resulting pressure distribution and lift due to the shape of the given solid boundaries. For the solution of the complete flowfield, however, viscous effects must be considered too, which for the attached flows will provide information such as the displacement thickness and the skin friction on the solid surface—or the drag force component due to this surface friction. Also, more advanced viscous methods should be capable of indicating whether the flow will have a tendency to detach (e. g., predicting location of separation points, or lines). The objective of this chapter is to provide a brief survey of some frequently occurring low-speed flowfields and to help the student to place in perspective the relative role of the potential-flow methods (presented in this book) and of the viscous effects in order to comprehend the complete real flowfield environment. Additionally, several simplified enhancements to the potential-flow model that account for some viscous effects will be surveyed.

The modifications presented in this chapter will begin with methods of

calculating the wake rollup, which from the classical potential-flow solution point of view was denoted as a “slight nonlinear effect.” The rest of the presented improvements (or modificatons) deal with efforts to include the effects of viscosity and some of them are logical extensions to the potential – flow model. Some others (e. g., modeling of two-dimensional flow separation) will clearly fall into the “daring and imaginative” category and their impor­tance is more in providing some explanation of the fluid-mechanical phenom­ena rather than being in such a stage that they can predict unknown flowfields.

In the following discussion, for the sake of simplicity, mainly the lifting characteristics of the experimental observations and the resulting flow models are presented. In a limited number of cases the drag force also is discussed but important effects such as side forces, moments and possible crosscoupling of the aerodynamic loads is omitted in favor of brevity. Therefore, the treatment of the various topics in this chapter is by no means complete or comprehensive and the reader is encouraged to further investigate any of the following topics in the referenced literature.

Principles of converting a steady-state, potential-flow solution into a time – dependent mode were summarized in Section 13.6, and were demonstrated in the sections that followed. The complexity of these examples in terms of geometry increased gradually and in the previous section the three-dimensional thin lifting-surface problem was illustrated. Although this method was capable of estimating the fluid dynamic lift, the calculation of the drag was indirect and inefficient from the computational point of view. Therefore, a similar conversion of a three-dimensional panel method into the unsteady mode can provide, first of all, the capability of treating thick and complex body shapes, and in addition the fluid dynamic loads will be obtained by a direct integration

of the pressure coefficients. Since the pressure coefficient is obtained by a local differentiation of the velocity potential (and not by summing the influence of all the panels) this approach yields an improved numerical efficiency. Also, the drag force is obtained as a component of the pressure coefficient integration and there is no need for a complicated estimation of the leading-edge suction force.

The following example is based on the conversion of a steady-state panel method based on constant-strength source and doublet elements12 11 (described in Section 12.5), which resulted in the time-dependent version,12 13 presented in this section. Also, it is recommended for the reader to be familiar with Sections 12.5, and 13.12 since some details mentioned in these sections are described here only briefly.

The method of the conversion is described schematically in Fig. 13.25, and the potential-flow solution will be included in a time-stepping loop that will start at t = 0. During each of the following time steps the strength of the latest wake row is computed by using the Kutta condition, and the previously shed wake vortex strengths will remain unchanged. Thus, at each time step, for N panels N equations will result with N unknown doublet strengths. If the geometry of the body does not change with time then the matrix is inverted only once. In a case when the body geometry does change (e. g., when a propeller rotates relative to a wing) the influence coefficients and matrix inversion are calculated at each time step. The description of the method, based on the eight step procedure, is then:

Choice of singularity element. The basic panel element used in this method has a constant-strength source and/or doublet, and the surface is also planar (but the doublet panels that are equivalent to a vortex ring can be twisted). Following the formulation of Section 9.4, the Dirichlet boundary condition on a thick body (e. g., Eq. (13.18)) can be reduced to the following form (see Eq. (12.29)):

N Nw N

2 Ckvk + 2 c, fi, + 2 Bkak = 0 (13.156)

* = 1 /=1 k = 1

which condition must hold at any moment t. This equation will be evaluated for each collocation point inside the body and the influence coefficients Ck, C, of the body and wake doublets, respectively, and Bk of the sources are calculated by the formulas of Section 10.4. (In this example only the Dirichlet boundary condition is described but with a similar treatment the Neumann condition can be applied to part or all of the panels.)

Kinematics. Let us establish an inertial frame of reference X, Y, Z, as shown in Fig. 13.39, such that this frame of reference is stationary while the airplane is moving to the left of the page. The flight path of the origin and the orientation of the x, y, z system are assumed to be known and the boundary

condition (Eq. (13.13a)) on the solid surface becomes

ЭФ

-г-= (V0 + vrel + ft X r) • n (13.157)

an

The kinematic velocity components at each point in the body frame due to the motion (V0 + vte, + ft X r) are given as [t/(/), V(t), W(t)] by Eq. (13.141). If the combined source/doublet method is used (see Section 13.2) then the Dirichlet boundary condition requires that the source strength is given by Eq.

(13.19) :

о = —n • (V0 + vrel + ft X r) (13.19)

Discretization and grid generation. In this phase the geometry of the body is divided into surface panel elements (see for example Fig. 12.22) and the panel comer points, collocation points (usually slightly inside the body) and the outward normal vectors n* are identified. Also, the counter к for each panel is assigned and a typical example of generating a wing grid and the unfolded patch are shown in Fig. 12.23.

The wake shedding procedure is described schematically by Fig. 13.40. A typical trailing-edge segment is shown with momentary upper fiu and lower ju, doublet strengths. The Kutta condition requires that the vorticity at the trailing edge stays zero:

Hw, = (j*u – M,), (13.158)

Thus, the strength of the latest wake panel /uWt is directly related to the wing’s (or body’s) unknown doublets. Note that the spanwise segment (parallel to the

trailing edge) of the latest wake panel (which is actually equivalent to a vortex ring) is placed in the interval covered by the trailing edge during the latest time step (of length Q ■ At). Usually it must be placed closer to the trailing edge within 0.2-0.3 of the above distance (see discussion about this topic at the beginning of Section 13.8.2). During the second time step the wing trailing edge has advanced and a new wake panel row can be created using the new aft points of the trailing edge. The previous (/ – Af)th wake row will remain in its previous location (as observed in the inertial frame) so that a continuous wake sheet is formed. The wake corner points will then be moved with the local velocity, in the wake rollup calculation phase. Once the wake panel is shed, its strength is unchanged (recall the Helmholtz theorems in Section 2.9), and the wake vortex carries no aerodynamic loads (and therefore moves with the local velocity). Thus, the strengths of all the previous wake panels are known from previous time steps. This shedding procedure is repeated at each time step and a row of new trailing-edge wake vortex rings are created (wake shedding procedure).

Influence coefficients. In order to specify the time-dependent boundary condition the kinematic conditions need to be known (from Eq. (13.141)) and a time-stepping loop (shown in Fig. 13.25) is initiated with I, as the time step counter:

Let us assume that at t = 0 the two coordinate systems x, y, z and X, Y, Z in Fig. 13.39 coincided and the wing was at rest. The calculation is initiated at t = At and the wake at this moment consists of one wake panel row (the wake panel row adjacent to the trailing edge in Fig. 13.40). The Dirichlet boundary condition (Eq. (13.156)) when specified, for example, at the ith panel’s collocation point (inside the body) is influenced by all the N body and Nw wake panels and will have the form:

N Nw N

2 Cl7tju* + 2 C„ni + 2 Bik°k = 0 (13.156a)

*=i /=і л=1

But the strength of all the wake panels is related to the unknown doublet values of the trailing-edge upper and lower panels, via the Kutta condition (Eq. 13.158). Therefore, by resubstituting the trailing-edge condition (see also a similar explanation in Section 12.5), this boundary condition can be reduced to include only the body’s unknown doublets and for the first time step it becomes:

2 AuJtk + 2 B№ak = 0 t = At (13.159)

k=1

where Aik = Cik if no wake is shed from this panel and Axk = Cik ± CH if it is shedding a wake panel.

During the subsequent time steps wake panels will be shed, but as noted, their strength is known from the previous computations. Thus, Eq. (13.159) is valid only for the first time step, and for t> At the influence of these wake doublets nw (excluding the latest row) must be included in the boundary conditions. So for all the other time steps Eq. (13.156a) will have the form

N Mw N

2 Aikfik + 2 Can, + 2 Влок = 0 t> At (13.160)

*=i /=i *=i

Note that now the wake counter Mw does not include the latest wake row.

Establish RHS vector. Since the source value is set by the value of the local kinematic velocity (Eq. (13.19)), the second and third terms in Eq. (13.160) are known at each time step and, therefore, can be transferred to the right-hand side of the equation. The RHS vector is then defined as

(Again, note that fi, w and ok are known.) In the case when the body geometry is not changing with time the bkt coefficients are calculated only once, but the ckt coefficients of the wake must be recomputed at each moment because of the wake’s time dependent rollup.

This matrix has a nonzero diagonal (akk = when the panel is not at the

trailing edge) and has a stable numerical solution.

The results of this matrix equation can be summarized in indicial form (for each collocation point к) as

If the shape of the body remains unchanged then the matrix inversion occurs only once. For time steps larger than 1 the calculation is reduced to

^ = 2>*/RHS, (13.164)

/=i

where akll are the coefficients of the inverted matrix. In situations when a large number of panels are used (more than 2000) then from the computational point of view it is often more economical to iterate for a new instantaneous solution of Eq. (13.163), at each time step, than to store the large inverted matrix a*/ in the memory.

and in the normal direction

where /, m are the local tangential coordinates (see Fig. 12.25). For example, the perturbation velocity component in the / direction can be formulated (e. g., by using central differences) as

Яі = ^гЛРі+і~ Pi-i) (13.166)

2 Л/

In most cases the panels do not have equal sizes and instead of this simple formula, a more elaborate differentiation can be used. The total velocity at collocation point к is the sum of the kinematic velocity plus the perturbation velocity:

Q* = [1/(0, V(0, wm • (/, rn, n)k + (q„ qm, qn)k (13.167)

where lk, mk, nk are the local panel coordinate directions (shown in Fig. 12.25) and of course the normal velocity component for a solid surface is zero. The pressure coefficient can now be computed for each panel using Eq. (13.28):

Here Q and p are the local fluid velocity and pressure values, ЭФ/dt = Эц/dt (since Ф, = 0) and pref is the far-field reference pressure and uref can be taken as the magnitude of the kinematic velocity as appears in Eq. (13.8):

Vref = -[Vo + ft X r]

or as the translation velocity of the origin V0. For nonlifting bodies the use of Eq. (13.28a) instead of Eq (13.168) is recommended when the body rotation axis is parallel to the direction of motion. (In the case of more complex motion the use of the pressure equation and selection of vref should be investigated more carefully.)

The contribution of an element with an area of A5* to the aerodynamic loads AF* is then

AF* = – CPkaPv2Kf)k AS* n* (13.170)

In many situations, off-body velocity field information is required too. This type of calculation can be done by using the velocity influence formulas of Chapter 10 (and the singularity distribution strengths of о and of p are known at this point).

Vortex wake rollup. Since the wake is force-free, each wake panel (or wake vortex ring) must move with the local stream velocity (Eq. (13.21a)). The local velocity is a result of the kinematic motion and the velocity components induced by the wake and body and are usually measured in the inertial frame of reference X, Y, Z, at each panel’s corner points. This velocity can be calculated (using the velocity influence formulas of Section 10.4.1 for the doublet and of Section 10.4.2 for the source panels) since the strength of all the singularity elements in the field is known at this point of the calculation.

To achieve the wake rollup, at each time step, the induced velocity (и, v, w)i at each wake panel corner point l is calculated and then the vortex elements are moved by

(Ajc, Ay, Az), = (и, v, w)t At (13.171)

Summary

The time-stepping solution is best described by the block diagram in Fig.

13.25. For cases with fixed geometry (e. g., a maneuvering airplane) the geometrical information, such as panel corner points, collocation points, and normal vectors, must be calculated first. Then the time-stepping loop begins and based on the motion kinematics the geometry of the wake panel row adjacent to the trailing edge is established. Once the geometry of the trailing edge area is known, the influence coefficients ak, of Eq. (13.62) can be calculated. Also, using the same kinematic velocity information (e. g. Eq.

(13.141) ) the body’s source strength (Eq. (13.19)) and the RHS vector of Eq. (13.161) are obtained. Next, the unknown doublet distribution is obtained and the surface velocity components and pressures are calculated. Prior to advancing to the next time step, the wake rollup procedure is performed and then the time is increased by At, the body is moved along the flight path, and the next time step is treated in a similar manner.

Some examples of using the unsteady, constant-strength singularity element based panel method of Ref. 12.13 are presented in the following paragraphs.

Example 1. Large-amplitude pitch oscillation of a NACA 0012 airfoil. The

previous examples on the pitch oscillations of an airfoil were obtained by thin airfoil methods which do not provide the detailed pressure distribution on the surface. In this case the computation’s are based on a thick-airfoil model and the two-dimensional results were obtained by using a large aspect ratio (Ж = 1000) rectangular wing. The lift and pitching-moment histograms, during a fairly large-amplitude pitch oscillation cycle, of this NACA 0012 two-dimensional airfoil are presented in Fig. 13.41. Comparison is made with experimental results of Ref. 13.14 for oscillations about the airfoil’s quarter-chord. The computations are reasonably close to the experimental values of the lift coefficient through the cycle. During the pitchdown motion, however, a limited flow separation reduces the lift of the airfoil in the experimental data. The shape of the pitching-moment loop is close to the experimental result with a small clockwise rotation. This is a result of the inaccuracy of computing the airfoil’s center of pressure, since only nine chordwise panels were used.

This example indicates, too, that if the flow stays attached over the airfoil then the Kutta condition based load calculation is applicable to engineering analysis even for these large trailing edge displacements.

Example 2. Sudden acceleration of an airplane configuration. The transient load on a thin airfoil that was suddenly set into motion was first reported during

the 1920s13 2 and only recently with the use of panel methods could this type of analysis be applied to more realistic airplane configurations. Such computation for a complex aircraft shape is presented in Fig. 13.42, and the panel grid consists of 706 panels per side of the model. The transient lift growth of this wing/canard combination differs somewhat from the monotonic lift increase of a single lifting surface as presented in Fig. 13.37. At the first moment the lift of the wing and canard grow at about the same rate, with the lift of the wing being slightly lower because of the canard-induced down wash. Then the wing’s lift increases beyond its steady-state value, since the canard wake has not yet reached the wing. At about Q^t/c — 1.0 the canard wake reaches the wing and its influence begins to reduce the wing’s lift. This behavior results in lift-overshoot, as shown in the figure.

Example 3. Helicopter rotor. The flexibility of this method can be demon­strated by rotating a pair of high aspect ratio, untwisted wings around the 2 axis, to simulate rotor aerodynamics. The trailing-edge vortices behind this two-bladed rotor, which was impulsively set into motion, are presented in Fig. 13.43. Similar information on wake trajectory and rollup, for more complex rotorcraft geom­etries and motions (including forward flight), can easily be calculated by this technique. The spanwise lift distribution on one rotor blade of Fig. 13.43, after one-quarter revolution (Aip = 90°), is presented in Fig. 13.44. The rotor for this example is untwisted and has a collective pitch angle of aR = 8°, to duplicate the geometry of the rotor tested by Caradonna and Tung.1315 The large difference between this spanwise loading (АЦі = 90°) and the experimental loading measured in Ref. 13.15, for a hovering rotor, are due to the undeveloped wake. This solution can be considerably improved by allowing about eight revolutions of the rotor, so that the wake induced flow will develop. This spiral vortex wake – induced downwash did reduce the spanwise lift distribution on the wake to values that are close to those measured by Caradonna and Tung,1315 as shown in Fig. 13.44 (by the “steady hover” line). The corresponding chordwise pressures, for three blade stations, are presented in Fig. 13.45. The computed pressures fall close to the measurements of Ref. 13.15 and the small deviations could be a result of the sparse panel grid used or could be due to experimental errors.

Increasing the complexity of the motion is fairly simple and the forward flight of this rotor with a generic body is shown in Fig. 13.46

Example 4. Coning motion of a generic airplane. The coning motion is described schematically in Fig. 13.47 for the generic airplane geometry that was modeled by 718 panels per side. In principle the x coordinate of the body system translates forward at a constant speed Q„ and the model angle of attack is set within this frame of reference. The rotation <j> is performed about the x axis, as shown in the figure. Computed and experimental normal force C2, side-force Cy, and rolling moment C, are presented in Fig. 13.48a and b. The aircraft model was rotated about its center of gravity at a rate of up to <ub/2I7 = 0.04. This rate is fairly low, but representative of possible aircraft flight conditions and was selected to match the experiments of Ref. 13.16. The normal force is not affected by this low rotation rate and both experimental and computed lines are close to being horizontal. For higher angles of attack, the computational results are lower than the experimental data due to the vortex lift of the strakes. The side-force, in this type of motion, is influenced by the side-slip of the vertical and horizontal tail surfaces. Consequently, the computed values of Cy, for the above angle-of-attack range, are close to the experimental data.

The computed rolling moment of the configuration Cy at a = 0 (Fig. 13.48c) is much larger than shown by the experiment. However, most important is that the trend of the curve slope (which is really the roll damping) becoming negative at the larger angles of attack is captured by the computation. This slope is also a function of the distance between the wing’s center of pressure and the rotation axis, and the error in computing this distance is probably the reason for the larger (computed) rolling moments.

PROBLEMS

13.1. Consider a two-dimensional version of the relative motion described in Fig. 13.1, between a body-fixed frame of reference (дг, z) and an inertial frame (X, Z) such that

(X0, Zu) = (-Vj, – Wjt) в = sin (Ot and

(Xn, Z„) = (-fA, — WQ
6 = to cos cot

(a) Use the chain rule to evaluate the derivatives Э/ЭХ, Э/ 3Z, and Э/dt in terms of the body coordinates.

(b) Using your results from (a) transform the Bernoulli equation

p„- p _ (УФ)2 ЭФ p 2 dt

into the (x, z) frame of reference.

13.2. The two-dimensional flat plate, shown in Fig. 13.49, is initially at rest and at t = 0+ it moves suddenly forward at a constant speed {/». Obtain the time – dependent circulation T(f) and lift L(t) of the flat plate using two chordwise lumped-vortex elements (select the vortex and collocation points as suggested in Section 11.1.1) with a discrete-vortex model for the wake and present your results graphically (as in Fig. 13.8).

(a) Study the effect of time step U„ At/c in the range U„ At/c = 0.02-0.2. Note that a smaller time step simulates a faster acceleration to the terminal speed U, o, and therefore has a physical effect on the results (in addition to the numerical effect).

(b) Study the effect of wake vortex positioning by placing the latest trailing-edge vortex at the beginning, center, and end of the interval covered by the trailing edge during the latest time step (see Fig. 13.49). Compare your results with the more accurate calculations in Fig. 13.34 (for Ж = <»).

13.3. Use the flat-plate model of the previous example to study the constant acceleration of a flat plate. Assume that the forward speed is U(t) = at and the time step is a At1 lie = 0.1. Calculate the time-dependent circulation F(t) and lift L(f) of the flat plate for several values of the acceleration a and present your results graphically (as in Fig. 13.8). For simplicity, place the latest vortex shed from the trailing edge at one-third of the distance covered by the trailing edge during that time step.

13.4. Convert any of the two-dimensional panel codes of Chapter 11 (e. g., a constant-strength doublet method) to the unsteady mode and validate it by calculating the lift and circulation after a sudden acceleration. (This problem requires a larger effort and can be given as a final project.)

The method of transforming steady-state flow-based numerical solutions into the time-dependent mode is described schematically in Fig. 13.25 and in the introduction to Section 13.10. In this section the same approach is applied to the three-dimensional thin lifting-surface problem. In this example, as in the case of the lifting surface of Section 12.3, the wing’s bound circulation and the

vortex wake will be modeled by vortex ring elements (see Fig. 13.30). The main advantage of this element is in the simple programming effort that it requires (although its computational efficiency can be further improved). Also, in this numerical example, the boundary conditions are specified on the actual wing surface, which can have camber and various planform shapes.

The solution is again based on the time-stepping technique, and at the beginning of the motion only the wing bound vortex rings exist (upper part of Fig. 13.30). Note that the closing segment of the trailing edge vortex elements in Fig. 13.30 will represent the starting vortex. Consequently, during the first time step there will be no wake panels and if the wing is represented by К unknown vortex rings (K = 8 in Fig. 13.30), then by specifying the zero normal flow boundary condition on the К collocation points a solution at t = At is possible. During the second time step, the wing is moved along its flight path and each trailing edge vortex panel sheds a wake panel with a vortex strength equal to its circulation in the previous time step (lower part of Fig. 13.30). With this model, therefore, we do not have to add an additional equation to enforce the Kelvin condition since the vortex ring formulation inherently fulfills this condition. During this second time step, there will be only one row of wake vortices, but with a known strength. Therefore, the wing К bound vortices can be calculated for this time step, too, by specifying the boundary condition on the same К collocation points. This time-stepping methodology can then be continued for any type of flight path and at each time step the vortex wake corner points can be moved by the local velocity, so that wake rollup can be simulated. One of the advantages of this wake model is in its simplicity, and its being equal in its formulation to the wing’s bound vortex rings (or constant-strength doublet panels); consequently, this wake model can be used for more advanced panel methods (and is actually used by the code PMARC9 6,9 7,12 13).

It is recommended that prior to reading this section the reader should be familiar with the steady-state solution presented in Section 12.3. The presenta­tion of the unsteady version of this method then proceeds with the same sequence (steps) of the previous sections. So, in general, the wing is divided into panel elements containing vortex ring singularities as shown in Fig. 13.30 and the solution procedure is as follows:

Choice of singularity element. The method by which the thin-wing planform is divided into panels is similar to what is described in Section 12.3 and shown schematically in Figs. 13.30 and 13.31. The leading segment of the vortex ring is placed on the panel’s quarter-chord line and the collocation point is at the center of the three-quarter chord line. The normal vector n is defined at this point, too, which falls at the center of the vortex ring. A positive Г is defined here according to the right-hand rotation rule, as shown in the figure.

j

From the numerical point of view these vortex ring elements are stored in rectangular patches (arrays) with i, j, indexing as shown by Fig. 13.31 (see also Fig. 12.10). The velocity induced at an arbitrary point P (де, у, z), by a typical vortex ring (see Fig. 13.32) at location ij can be computed by applying the vortex line routine VORTXL (Eq. (10.116)) to the ring’s four segments:

(и,, vu w,) = VORTXL(x, y, z, xu, уіф zQ, xiJ+l, yiJ+l, z, v+1, Г,.,)

The velocity induced by the four vortex segments is then

(u, v, w) = (uit vu w^) + (u2, v2, w2) + (m3, v3, w3) + (m4, v4, w4) (13.133)

It is convenient to include these computations in a subroutine (see Eqs.

(10.117) and (12.17)) such that

(и, v, tv) = VORING(x, y, z, i, j, Г)

Note that in this formulation it is assumed that by specifying the i, j counters the (jc, y, z) coordinates of this panel are automatically identified (see Fig. 12.10).

The use of this subroutine can considerably shorten the programming effort; however, for the inner vortex rings the velocity induced by the vortex segment is computed twice. For the sake of simplicity this routine will be used for this problem, but more advanced programming can easily correct this loss of computational effort.

It is recommended at this point, too, to calculate the velocity induced by the trailing vortex segments only (the vortex lines parallel to the free stream,

as in Fig. 12.5). This information is needed for the induced drag computations and if done at this point will only slightly increase the computational effort. The influence of the trailing segments is obtained by simply omitting the (ut, Vi, Wj) + (и3) v3, и’з) part from Eq. (13.133):

(u, v, w)* = (u2, v2, w2) + (u4, ^4, w„) (13.135)

and it is assumed that (и, v, w)* is automatically obtained as a byproduct of subroutine VORING.

Kinematics. Let us establish an inertial frame of reference X, Y, Z, as shown in Fig. 13.1, such that this frame of reference is stationary while the wing is moving to the left of the page. The flight path of the origin and the orientation of the x, y, z system is assumed to be known and is prescribed as:

X0 = X0(t) Yn=Y0(t) Zu = Z0(t) (13.136)

ф = ф(і) в = 6(t) ф = p(t) (13.136a)

and the momentary velocity of the origin and the frame’s rotations about the axes are

*о = ЛГо(0 Z0 = Z0(t) (13.137)

Ф = Ф( t) e = d(t) Ф = Ф(1) (13.137a)

For example, if the wing moves to the left (parallel to the X axis) at a constant speed of f/„ and has a constant sideslip (parallel to the Y axis) of 14, then

and

It is useful to establish a transformation between the two coordinate systems (Eq. (13.7a)) that depends on the translation of the origin and the orientation of the body frame of reference:

/X-Хо

=/(0, в, V)(

Z-Z0

This transformation, without the translation part, can be used also for the velocity transformation. To define such a three-dimensional transformation uniquely, the order of rotation must be specified. For example, if we adopt the order of rotation such that first we rotate about the z axis, then about the у

axis, and finally about the x axis then the transformation becomes:

Rotation by ip (sideslip)

Rotation by в (angle of attack)

Rotation by ф (roil angle)

where U3, КЗ, W3 are the velocity component, observed in the x, y, z frame due to the translation of the origin.

The time-dependent kinematic velocity components U(l), K(f), W(t), in the x, y, z frame are then a combination of the translation velocity and the rotation of the body frame of reference:

(13.141)

Since the instantaneous rotation and translation rates are known, these kinematic terms are known, too, at each time step.

Discretization and grid generation. The method by which the thin wing planform is divided into elements is the same as presented in Section 12.3 and is shown in Figs. 13.30 and 13.31. Some typical panel elements are also shown in Figs. 12.8-12.10. Also, if only the wing’s semispan is modeled then the mirror image method must be used to account for the other semispan. The leading segment of the vortex ring is placed on the panel’s quarter-chord line and the collocation point is at the center of the three-quarter chord line (see Fig. 13.31). The lifting surface is usually given by z = r)(x, y) and is divided into N span wise and M chord wise panels. Using a procedure such as shown in Fig. 12.13 allows the scanning and calculation of geometrical information such as the panel area Sljt normal vector n,; and the coordinates of the collocation points. A simple and fairly general method for evaluating the normal vector is shown in Fig. 12.11. The panel opposite corner points define two vectors A
and B, and their vector product will point in the direction of n

A XB П-|АХВ|

A positive Г for a vortex ring is defined here using the right-hand rotation convention, as shown in Fig. 13.31. Also, for increased surface curvature the above-described vortex rings will not be placed exactly on the lifting surface, and a finer grid needs to be used, or the wing surface can be redefined accordingly. By placing the leading segment of the vortex ring at the quarter-chord line of the panel the two-dimensional Kutta condition is satisfied along the chord (recall the lumped-vortex element).

At this point also the wake shedding procedure must be addressed. Consider a typical trailing-edge vortex ring placed on the last panel row (as shown in Fig. 13.33). The trailing segment (parallel to the trailing edge) is placed in the interval covered by the trailing edge during the latest time step (of length Q ■ At). Usually it must be placed closer to the trailing edge within

0. 2-0.3 of the above distance (see discussion about this topic at the beginning of Section 13.8.2). The wake vortex ring corner points must be created at each time step, such that at the first time step only the two aft points of the vortex ring are created. Therefore, during the first time step there are no free wake elements and the trailing vortex segment of the trailing edge vortex ring represents the starting vortex. During the second time step the wing trailing edge has advanced and a wake vortex ring can be created using the new aft points of the trailing edge vortex ring and the two points where these points were during the previous time step (see Fig. 13.30). This shedding procedure is repeated at each time step and a set of new trailing-edge wake vortex rings are created (wake shedding procedure).

The strength of the most recently shed wake vortex ring (Г^( in Fig. 13.33) is set equal to the strength of the shedding vortex ГТЕ.(_Л( (placed at the trailing edge) in the previous time step (as if it was shed from the trailing edge

FIGURE 13.33

Nomenclature for the wake shed­ding procedure at a typical trailing edge panel.

and left to flow with the local velocity):

ІЧ = Гт. е.,„лі (13.142)

Once the wake vortex is shed, its strength is unchanged (recall the Helmholtz theorems in Section 2.9), and the wake vortex carries no aerodynamic loads (and therefore moves with the local velocity).

This is the unsteady equivalent of the Kutta condition. For the steady-state flow conditions, all wake panels shed from a particular trailing edge panel will have the same vortex strength, which is equal to the strength of the shedding panel. Thus, the spanwise-oriented vortex lines of the adjacent vortex rings will cancel each other and only a horseshoe-like vortex will remain.

Influence coefficients. At this point, the zero normal flow across the solid – surface boundary condition is implemented. In order to specify this condition, the kinematic conditions must be known and the time-stepping loop (shown in Fig. 13.25) is initiated. Let us select again I, as the time-step counter, so that the momentary time is

t = /, • Af

Let us assume that at t = 0 the two coordinate systems x, y, z and X, Y, Z coincided and the wing was at rest. The calculation is initiated at t = At and the wake at this moment consists of the vortex line created by the trailing segments of the trailing edge vortex rings (Fig. 13.30). The location of the trailing edge at t = 0 and at t = Af, which is needed for specifying the wake panels’ corner points, is obtained by using coordinate transformations such as Eq. (13.140).

The normal velocity component at each point on the camberline is a combination of the self-induced velocity, the kinematic velocity, and the wake-induced velocity. The self-induced part can be represented by a combination of influence coefficients, as in the steady-state flow case. If the shape of the wing rj(x, у ) remains constant with time then these coefficients will be evaluated only once. The normal velocity component due to the motion of the wing is known from the kinematic equations (Eqs. (13.13a), or (13.141)) and will be transferred to the right-hand side (RHS*) of the equation. The strength of the other wake vortices is known from the previous time steps and the wake-induced normal velocity on each panel will be transferred to the right-hand side too.

Let us establish a collocation point scanning procedure (similar to that of Section 12.3) that takes the first chordwise row where і = 1 and scans spanwise with j = l—*N and so on (recall that we have M chordwise panels—see Fig. 12.10). This procedure can be described by two “DO loops” shown in Fig.

12.13. As the panel scanning begins, a sequential counter assigns a value К to each panel (the sequence of К is shown in Fig. 12.14) which will have values from 1 —* M x N.

Once the collocation point scanning has started, К — 1 (which is point

(і = 1, j = 1) on Fig. 12.12. The velocity induced by the first vortex ring is then
(u, v, wOu = VORING(jc, y, z, i = , j = ,Tl = 1.0)

Note that a unit strength vortex is used for evaluating the influence coefficient

«и = («, v, tv),! • nj (13.143)

To scan all the vortex rings influencing this point, an inner scanning loop is needed with the counter L – 1 —>N X M (see Fig. 12.13). Thus, at this point, the К counter is at point 1, and the L counter will scan all the vortex rings on the wing surface, and all the influence coefficients aXL are computed (also, in Eq. (13.143) the ( )ц index means К = 1, L = 1):

aiL = (u, v, w)1L – Пі

and for the K, Lth panel

aKL = («, v, w)KL • nK (13.143a)

As mentioned before, parallel to the computation of the aKL coefficients, the normal velocity component induced by the streamwise segments of the wing vortex rings can also be computed. These bKL coefficients, which will be used for the induced drag calculation, are calculated by using Eq. (13.135)

bKL = (u, v, w)*KL • nK (13.144)

and it is assumed that these coefficients are a byproduct of the aKL calculations and do not require additional computational effort.

This procedure continues until all the collocation points have been scanned and a Fortran example for this influence coefficient calculation is presented in Fig. 12.13 of Chapter 12.

Establish RHS vector. Specifying the zero normal velocity boundary condition on the surface (Q„K = 0) at an arbitrary collocation point К results in

QnK = «апГ| + ак2Г2 + акзГ3 + • • • + aKmTm

+ [f/(f) + Uw> Y(0 + vw> lY(t) + • и* = 0

where [U(t), V(t), W(r)k are the time-dependent kinematic velocity com­ponents due to the motion of the wing (Eq. (13.141)), (uw, vw, ww)K are the velocity components induced by the wake vortices, and m = M x N. The wake influence can be calculated using Eq. (13.134) since the location of all vortex points is known (including the wake vortex points). Since these terms are known at each time step, they can be transferred to the right-hand side of the equation. Consequently, the right-hand side (RFIS) is defined as

RHS* = -[l/(f) + uw, V{t) + vw, W(t) + wwK • nK (13.145)

Solve linear set of equations. Once the computations of the influence coefficients and the right-hand side vector are completed, the zero normal flow

boundary condition on all the wing’s collocation points will result in the following set of algebraic equations:

(Recall that К is the vertical and L is the horizontal matrix counter and the order of this matrix is m = M X N.)

The results of this matrix equation can be summarized in indicial form (for each collocation point К) as

m

2>KLrL = RHS* (13.146)

L= 1

If the shape of the wing remains unchanged then the matrix inversion occurs only once. For time steps larger than 1 the calculation is reduced to

Г* = 2 e«.RHSt (13.147)

L=1

where a^L. are the coefficients of the inverted matrix.

Computation of velocity components, pressures and loads. For the pressure distribution calculations the local circulation is needed, which for the leading – edge panel is equal to Гц but for all the elements behind it is equal to the difference Г^ — Г,_іг The fluid dynamic loads then can be computed by using the Bernoulli equation (Eq. (13.24)) and the pressure difference is given by Eq. (13.123):

Ар = A – P. – p[(f)_ – (f)(+ (^)_ – (™) J

t, and t are the panel tangential vectors in the і and j directions (of course, these vectors are different for each panel and the ij subscript from t,-.. is dropped for the sake of simplicity).

The velocity-potential time derivative is obtained by using the definition ф*=± y/2 dl and by integrating from the leading edge. Since for this vortex

ring model ДФ = Г,

(13.149)

Substituting these terms into the pressure difference equation results in Ару = pj[I/(f) + uw, V(t) + %, W(t) + wwtj • t, ‘ ^ –

Г _ j-*. ^ і

+ [U(t) + uw, V(t) + vw, W(t) + Hv],; • x, + J( r./j (13.150)

The contribution of this panel to the loads, resolved along the three body axes, is then

AF=-(ApAS)0n& (13.151)

The total forces and moments are then obtained by adding the contribution of the individual panels.

The total force obtained by this pressure difference integration will have some of the thin lifting-surface problems since it does not account for the leading-edge suction force. In general, the lifting properties of the wing will be estimated adequately by this method but the induced drag will be overesti­mated. Also, in the case of an arbitrary motion, the definition of lift and drag is more difficult and even the definition of a reference velocity (e. g., free-stream velocity) is not always simple. For example, presenting the pressure coefficient data on a helicopter blade in forward flight can be based on the local blade velocity or on the helicopter flight speed. So for the simplicity of this discussion on the induced drag, we shall limit the motion of the lifting surface such that it moves forward along a straight line without sideslip (but the forward speed may vary). The induced drag is then the force component parallel to the flight direction and each panel contribution is

A A, = p[(»,nd + wv),;,(r, y – Г,-еу) А Ьи + -| Г Ц А Бц sin or,,] (13.152)

and if the panel is at the leading edge then

(13.152a)

where ay; is the panel’s angle of attack relative to the free-stream direction. Also, the first term here is due to the downwash induced by the wing’s streamwise vortex lines ^ and due to the wake ww and the second term is
due to the fluid acceleration. The induced downwash wind at each collocation point i, j is computed by summing up the velocity induced by all the trailing segments of the wing bound vortices. This can be done during the phase of the influence coefficient computation (Eq. (13.144)) by using the VORING routine with the influence of the spanwise vortex segments turned off. This procedure can be summarized by the following matrix formulations where all the bKL and the are known:

where again m = N x M.

The main difficulty in the induced drag calculation for a general motion lies in the identification of the force component which will be designated as drag. Once this problem is resolved then the above method can be extended to more complex wing motions (and then angles such as <Хц in Eq. (13.152) must be defined).

Vortex wake rollup Since the vortex wake is force-free, each vortex must move with the local stream velocity (Eq. (13.21a)). The local velocity is a result of the velocity components induced by the wake and wing, and is usually measured in the inertial frame of reference X, Y, Z, at each vortex ring corner point. To achieve the vortex wake rollup, at each time step the induced velocity (и, v, w), is calculated and then the vortex elements are moved by

(Ax, Ay, Az)t = (u, v, w). At (13.153)

The velocity induced at each wake vortex point is a combination of the wing and wake influence and can be obtained by using the same influence routine (Eq. (13.134)):

m

(и, V, w)i = 2 VORING (xt, yh zt, i, j, TK)

K=l

Nw

+ 2 VORING (X/, У/, Zi, iw, jw, rW) (13.154)

*=i

and there are m wing panels and Nw wake panels.

In the case of a strong wake rollup the size of the wake vortex ring can increase (or be stretched) and if a vortex line segment length increases its strength must be reduced (from the angular momentum point of view). For the methods presented in this section it is assumed that this stretching is small and therefore is not accounted for.

Summary

The solution procedure is described schematically by the flowchart of Fig.

13.25. In principle, at each time step the motion kinematics is calculated (Eqs.

(13.141) ), the location of the latest wake vortex ring is established, and the RHS, vector is calculated. The influence coefficients appearing in Eq. (13.146) are calculated only during the first time step and the matrix is inverted. At later time steps, the wing vortex distribution can be calculated by the momentary RHS, vector, using Eq. (13.147). Once the vortex distribution is obtained the pressures and loads are calculated, using Eq. (13.151). At the end of each time step, the wake vortex ring corner point locations are updated due to the velocity induced by the flow field.

A student program based on this algorithm is given in Appendix D, Program No. 14.

that the * coordinate is parallel to the motion and the kinematic velocity components of Eq. (13.141) become [U(t), 0, 0]. The angle of attack effect is taken care of by pitching the wing in the body frame of reference and for the planar wing then all the normal vectors will be n = (sin a, 0, cos a). Consequently, the RHS vector of Eq. (13.145) becomes

RHS* = —{[£/(/) + uw sin a + ww cos a}K (13.155)

and here the wake influence will change with time. The rest of the time-stepping solution is as described previously in this section. For the numerical investigation the wing is divided into 4 chordwise and 13 spanwise equally spaced panels, and the time step is L4 At/c = 1/16.

Following the results of Ref. 13.12 the transient lift coefficient variation with time for rectangular wings with various aspect ratios is presented in Fig. 13.34. The duration of the first time step actually represents the time of the acceleration during which the ЭФ/dt term is large. Immediately after the wing has reached its steady-state speed (t/„) the lift drops due to the influence of the starting vortex and most of the lift is a result of the ЭФ/dt term (due to the

{/,/

c

change in the downwash of the starting vortex—see Fig. 13.8, too). Also, the initial lift loss and the length of the transient seem to decrease with a reduction in the wing aspect ratio (due to the presence of the trailing vortex wake).

The transient drag coefficient variation with time for the same rectangular wings is presented in Fig. 13.35. Recall that this is the inviscid (induced) drag and it is zero for the two-dimensional wing (/R = ^-). Consequently, the larger aspect ratio wings will experience the largest increase in the drag due to the downwash of the starting vortices. The length of the transient is similar to the results of the previous figure, that is, a smaller aspect ratio wing will reach steady state in a shorter distance (in chord lengths).

The above drag calculation results allow us to investigate the components of Eq. (13.152). For example, Fig. 13.36 depicts the drag C0l due to induced downwash (first term in Eq. (13.152)) and due to the fluid acceleration term C„, (which is the second term in Eq. (13.152)) for a rectangular wing with an aspect ratio of 8. At the beginning of the motion, most of the drag is due to the ЭФ/dt term, but later the steady-state induced-drag portion develops to its full value.

The effect of wing aspect ratio on the nondimensional transient lift of

uncambered, rectangular wings that were suddenly set into a constant-speed forward flight is shown in Fig. 13.37. This figure is very useful for validating a new calculation scheme, and is sensitive to the spacing of the latest wake vortex from the trailing edge (actually this is one method to establish the distance of the trailing edge vortex behind the trailing edge for a given time step). For comparison, the results of Wagner13 2 for the two-dimensional case are presented as well (in his case the acceleration time is zero and the lift at t = 0+ is infinite). It is clear from this figure, too, that the length of the transient and the loss of initial lift decreases with decreasing wing aspect ratio. The difference between the computed curve and the classical results of Wagner can be contributed to the finite acceleration rate during the first time step. The effect of this finite acceleration is to increase the lift sharply during the acceleration and to increase it moderately later (this effect of finite acceleration is discussed in Ref. 13.13).

Example 2. Heaving oscillations of a rectangular wing. As a final example this method is used to simulate the heaving oscillations of a rectangular wing near the ground. The boundary conditions for this case were established exactly as in the example of Section 13.9.1. The ground effect is obtained by the mirror image method and the results1312 for an aspect ratio of 4, planar rectangular wing, are presented in Fig. 13.38. The upper portion shows the effect of frequency on the lift without the presence of the ground and the loads increase with increased frequency. The lower portion of the figure depicts the loads for the same motion, but with the ground effect. This case was generated to study the loads on the front wing of a race car due to the heaving oscillation of the body and the data

indicates that the ground effect does magnify the amplitude of the aerodynamic loads.

The potential flow examples, as presented in Chapters 3 and 4, indicate that the solution for lifting flows is not unique for a given set of boundary conditions. This difficulty was resolved by requiring that the flow leave smoothly at the trailing edge of two-dimensional airfoils thereby fixing the amount of circulation generated by the airfoil. The above two-dimensional Kutta condition was almost automatically extended to the three-dimensional steady-state case and in this chapter was used for unsteady flows as well. Although from the mathematical point of view a condition to fix the amount of circulation is required, it is not obvious that this condition is the best candidate. However, prior to arriving at any conclusion in this regard, let us use the method of this section to study the wake rollup behind a thin airfoil undergoing a small-amplitude heaving oscillation.

FIGURE 13.28

Variation of the vertical displacement and normal force coefficient during one heaving cycle. From Katz, J. and Weihs, D., “Behavior of Vortex Wakes from Oscillatory Airfoils”, Journal of Aircraft, Vol. 15, No. 12, 1978. Reprinted with permission. Copyright AIAA.

Consider the small-amplitude heaving oscillation of the flat plate shown in the inset to Fig. 13.28. Assume that the motion of the origin of the x, z coordinates is given by

Xa =~UJ X0 (f) = – t/„

Zq = ~ho sin (Ot Z0(0 = —h0a) cos Ш
0=0 0=0

The time-dependent kinematic velocity components of Eq. (13.116) then become

Since the normal vector to the flat plate is n = (0,1) the right-hand side (downwash) vector of Eq. (13.117) becomes

RHS, = ~(h0(o sin Ш + ww)

The wake-induced downwash is obtained by using Eq. (13.113) and then at each moment Eq. (13.118) is solved for the airfoil vortex distribution. The

pressure difference is then obtained by using Eq. (13.127) and the wake rollup is obtained by moving the wake vortices with the local induced velocity (Eq. (13.132)). The rest of the details are as presented in the previous section and results of the wake rollup computation for the flat plate oscillating at various frequencies is shown in Fig. 13.29. The comparison in this figure indicates that up to a high reduced frequency of (oc/2Q„ = 8.5 the calculated wake shape is similar to the results of flow visualizations. Since the wake shape is a direct result of the airfoil’s circulation history and the calculated wake shape is similar to the experimentally observed shape, it is safe to assume that the calculated airfoil’s lift history is similar to the experimental one (which was not measured in this case). As an example, the vertical load Cz on the airfoil during one cycle is presented in Fig. 13.28 next to the motion history (note the phase shift due to the wake influence).

Now that we have generated a “good example” in favor of the unsteady Kutta condition, let us investigate some possible parameters affecting its validity. It is clear that conditions such as very high oscillation frequency, large amplitudes, and large angles of attack will cause some trailing edge separation. Such local flow separation automatically violates the Kutta condition, but in practice may not have a noticeable effect on the lift. On the other hand it may cause a lag in the aerodynamic loads. Experimental investigations of the unsteady Kutta condition13 7-1310 usually indicate that the streamlines do not leave parallel to the trailing edge at reduced frequencies of (Oc/2Q„>0.6, but the lift and pressure distributions are not affected in a visible manner even at higher frequencies. These experiments were based on small-amplitude oscilla­tions of airfoils where the trailing-edge vertical displacement was small.

So, based on the indirect results of Fig. 13.29 and some cited references we can try to establish some guidelines for the boundaries for the validity of the unsteady Kutta condition. First and most important, large angles of attack where trailing edge separation begins to develop must be avoided. Also, it is clear that as the reduced frequency increases the “allowed” trailing edge displacement amplitude (e. g., h0 in the previous example) must be smaller. So, for example, with h0 = 0.1c and with reduced frequencies of up to (od2Qx = 1.0 calculations based on the Kutta condition may provide reasonable load calculations. The vertical kinematic velocity of the trailing edge (e. g., hJU,„ in the previous example) is an important parameter, too, and in the case of the highest frequency oscillation in Fig. 13.29 it has a value of about 0.35. So in addition to the previously mentioned limits on the reduced frequency and trailing-edge amplitude, if we limit ourselves to trailing-edge vertical displace­ment velocity of hJU. « 1, then for practical purposes we can assume that the unsteady Kutta condition is valid. Also, we must remember that characteristic airplane maneuvers will fall in a category where the reduced frequency is far less than 1, and therefore the use of the Kutta condition is justified in most cases. For a rapidly pitching helicopter rotor in forward flight this may not be the case!

The above discussion was aimed primarily at the lift calculation; however, the lag (due to viscous effects) in the adjustment of the flow at the trailing edge may cause some lags in the aerodynamic loads, and this effect is still not explored sufficiently.

As was mentioned earlier, with only a few minor modifications, steady-state solution techniques can be updated to treat unsteady flows. As a first numerical example, the discrete-vortex model of thin lifting airfoils (Section 11.1.1) will be modified and the three areas of the program that will be affected are:

1. The normal velocity component on the solid boundary should include the components of the unsteady motion as well (as in Eq. (13.13a))—this is a minor and local modification to the steady-state program.

2. Similar corrections due to the unsteady motion should be included in the pressure calculations (using the modified Bernoulli equation, Eq. (13.23))— again this is a limited local modification.

3. An unsteady-wake model has to be established (e. g., the time-stepping discrete-vortex model as presented in Section 13.8.2). Note that such a wake model can be added on, in a simple manner, to the steady-state solvers.

(In this example the discrete vortex model for the thin lifting airfoil of Section

11.1.1 will be transformed to the unsteady mode. But any of the methods presented in Chapter 11 can be modified easily and can be given as a student project.)

The mechanics of such an upgrade are demonstrated by the generic flow chart of Fig. 13.25. Comparing this diagram with the steady-state diagram of Fig. 9.15 reveals that, first, a time-stepping loop has to be established (only one programming statement). Then a new element appears (the flight path information block), which has all the kinematic information on the body’s or

wing’s motion. The rest of the program will follow the methodology of Chapter 11, and the only additional block is the “wake rollup” block that will perform the wake rollup at each time step. Consequently, we shall follow the same sequence used in the previous numerical chapters.

Choice of singularity element. For this discrete-vortex method the lumped – vortex element is selected and its influence is given in Section 11.1.1:

(13.102)

where

rj = (x – Xj)2 + (z – z,)2

Thus, the velocity at an arbitrary point (x, z) due to a vortex element of circulation Гу located at (x,, zj) is given by this equation. This can be included in a subroutine, which was defined by Eq. (11.2):

(и, и») = VOR2D (Гу,, x, z, xy, z,) (11.2)

Using this lumped-vortex element, the airfoil will be represented by a set of discrete vortices placed on the camberline, as shown in Fig. 13.26. If the airfoil’s circulation changes with time, then vortex wake elements are shed at the trailing edge and the wake will be modeled by using the same discrete vortex model (Fig. 13.26).

Kinematics. Let us establish an inertial frame of reference X, Z, shown in Fig.

13.27, such that this frame of reference is stationary while the airfoil is moving to the left of the page. Next, the airfoil’s camberline is placed in a moving frame of reference x, z with the leading edge at the origin. The flight path of the origin and the orientation of the x, z system is prescribed as

*o = *o(0 Z0 = Z0(r) 0 = 6(0 (13.103)

and the instantaneous velocity of the origin and the frame’s rotation 0 about

FIGURE 13.26

x Discrete vortex model tor the un­steady thin airfoil problem (shown during the first time step, t = At)

the у axis are

X0*=X0(t) Z0 = Z0(f) 0=0(f) (13.104)

For example, if the airfoil moves to the left at a constant speed of t/„ and sinks at a speed of W„ then

X0=-U„t Y0(f) = -£4

Z0 = – w„t Z0(t) = – W. (13.105)

0 = 0, 0 = 0

Or if the airfoil flies at a constant forward speed £4, and performs pitch oscillations at a frequency со about the у axis, then

X0=-UJ, X0 (t)=-U„

Zo = 0, Z0(f) = о (13.106)

0 = sin cof, 0=co cos (at

It is useful to establish a transformation between the two coordinate systems shown in the figure such that

and similarly the transformed velocity components are

(X_/ cos 0(f) sin 0(f) /i Z/ v—sin 0(f) cos 0(f) Ai/

The inverse transformation is also useful, and the velocity components U„ W, observed in the x, z frame due to the translation of the origin are

/t/Л _/cos 0(f) – sin 0(f)/-Aro sin0(f) cos 0(f) A—Z0/

Discretization and grid generation. At this phase the thin-airfoil camberline (Fig. 13.26) is divided into N subpanels, which may be equal in length. The N vortex points (xj, Zj) will be placed at the quarter-chord of each planar panel and the zero normal flow boundary condition can be fulfilled on the camberline at the three-quarter chord point of each panel. These N collocation points (*,, Zj) and the corresponding N normal vectors n, along with the vortex points can be computed numerically or supplied as an input file. The normal n,- pointing outward at each of these points is found in the x-z frame from the surface shape tj(x), as shown in Fig. 13.27:

where the angle a, is shown in Fig. 13.26 for panel number 4. Similarly the tangential vector t, is

tі = (cos <*,, —sin a,) (13.111)

Since the lumped vortex element inherently fulfills the Kutta condition for each panel, no additional specification of this condition is required.

If the airfoil’s geometry does not change with time, then the calculation of the vortex points and collocation points can be done before the beginning of the time loop, as shown in Fig. 13.25 (where the box “definition of geometry” is placed outside of the time stepping loop).

Influence coefficients. At this point, the zero normal flow across the solid – surface boundary condition is implemented. In order to specify this condition, the kinematic conditions need to be known and the time-stepping loop (shown in Fig. 13.25) is initiated. Let us select I, as the time step counter, so that the momentary time is

t = • At

For simplicity, we assume that at t = 0 the two coordinate systems x, z and X, Z coincided and the airfoil was at rest (hence there was no wake). Consequently, the calculation begins at t = At (Fig. 13.26) and the wake at this moment consists of a single vortex Г^, which is placed along the path of the trailing edge (see Fig. 13.16, too). The location of the trailing edge at t = 0 and at t = At is obtained by using the transformations of Eqs. (13.107). The wake vortex is then placed usually at 0.2-0.3 of the distance covered by the trailing edge during the latest time step (see the discussion on this subject in the beginning of Section 13.8.2)

In general, the normal velocity component at each point on the camberline is a combination of the self-induced velocity, the kinematic velocity, and the wake-induced velocity. The self-induced part can be represented by a combination of influence coefficients, as in the steady-state flow case. If the shape of the airfoil rj(x) remains constant with time then these coefficients will be evaluated only once. The normal velocity component due to the motion of the airfoil is known from the kinematic equations (Eqs. (13.13a)) and will be transferred to the right-hand side of the equation. The velocity induced by the most recent wake vortex is unknown and will be resolved by adding an additional equation (the Kelvin condition). The strength of the other wake vortices is known from the previous time steps (for the general case when /, > 1 but since at t = At only one wake vortex is present, the rest of the wake contribution for the first time step is zero) and their effect on the normal velocity will be transferred to the right-hand side, as well.

To formulate the momentary boundary condition, let us use Eq. (13.13a), and for simplicity we allow only one component of the relative velocity vrei = (0, drf/dt) limited to small amplitudes, within the coordinate system x-z. Also, it is convenient to divide the perturbation potential into an

airfoil potential Фв and a wake potential Ф№ and both parts of the velocity potential will be modeled by discrete vortex elements of circulation Г. Consequently, the boundary condition of Eq. (13.13a) becomes

(УФВ + УФ*, – V0 – vreI – П X r) • n = 0 (13.112)

To establish the self-induced portion of the normal velocity (УФВ • n in Eq. (13.112)), at each collocation point, consider the velocity induced by the airfoil’s Tyth element at the first collocation point (in order to get the influence due to a unit strength Г, assume Г, = 1):

(u, w)v = VOR2D (Г, = 1 ,xuzu Xj, zj) (13.113)

The influence coefficient a(> is defined as the velocity component induced by the airfoil’s unit strength Г, element normal to the surface (at collocation point /). Consequently, the contribution of a unit strength singularity element j, at collocation point 1 is

au = (u, w)u – nj (13.114)

The induced normal velocity component qnl, at collocation point No. 1, due to all N vortex elements and the latest wake vortex is therefore

q„i — ЯцГі + а12Г2 + аі3Г3 + • • • + а^уГ/у + aiWrWi

Note that the strength of the airfoil vortices T; and of the latest wake vortex TW( is unknown at this point (see also Fig. 13.27 where T^, = rw„).

Fulfilling the boundary condition on the surface requires that at each collocation point the normal velocity component will vanish (Eq. (13.112)). Specifying this condition for the /th collocation point yields: а,]Г, + а,2Г2 + а,3Г3 + • • • + a, NrN + diw^w, + [C(t) + uw, W(f) + hv], • n, =0

(13.115)

where [U{t), W{t)j are the kinematic velocity components due to the motion of the airfoil and (uw, ww)t are the velocity components induced by the wake vortices (except the latest wake vortex—shown in Fig. 13.27). The wake influence can be calculated using Eq. (13.102) since the location of all wake vortex points is known. The time-dependent kinematic velocity components U(t), W(f), (see Eq. (13.49), too) are calculated with the help of Eq. (13.109):

/c°s 0(0 – sin 0(f)/-*o [ g

W(t)/ Vsin0(r) cos 0(f) A-Z0/

Since these terms are known at each time step, they can be transferred to the right-hand side of the equation. Consequently, the right-hand side (RHS) is defined as

Specifying the boundary condition for each of the collocation points results in the following set of algebraic equations:

RHSj ‘

rhs2

RHS*

m-A/) where the last equation represents the Kelvin condition,

Г(г) – Г(t – At) + TWi = 0 (13.119)

and the instantaneous airfoil circulation is the sum of all the airfoil’s vortices

Г(0 = £г,. (13.120)

This influences coefficient calculation procedure can be accomplished by using two “DO loops” where the outer loop scans the collocation points (xt, z,) and the inner scans the vortices Г,. If the airfoil camberline geometry does not vary with time [r/(jc, t) = r/(jc)], the influence coefficient calculation needs to be carried out only once at the beginning of the computations.

Establish RHS vector. The right-hand side vector, the normal component of the kinematic velocity and the wake induced velocity can be computed within the outer loop of the previously described “DO loops” by using Eq. (13.117).

Solve linear set of equations. The results of the previous calculations (shown by Eq. (13.118)) can be summarized in indicial form (for each collocation point i) as

N + l

2«,yr, = RHS, (13.121)

y=1

Again, if the shape of the airfoil remains unchanged then the matrix inversion occurs only once. For time steps larger than 1, the calculation is reduced to

N+l

r,= 2a,7’RHS, (13.122)

/ = 1

where ajj’ are the coefficients of the inverted matrix.

Computation of velocity components, pressures and loads. The resulting pressures and loads can be computed by using the Bernoulli equation

(Eq. (13.24)) near the panel surface.

Prcf-P _ Q2 ^ref ЭФ

p 2 2 dt

The pressure difference between the camberline upper and lower surfaces is then

(13.123)

and the tangential velocity component Q, is found from

ЗФ

a = [U(t) + uw, W(t) + иv], • T; ± ■— (13.124)

‘ otj

where the ± sign stands for above and under the surface, respectively. The tangential derivative of the thin airfoil potential can be approximated as

±эф,±Х„±Л

Эт) 2 2 A/у

and here A/, is the )th panel length. (Note that here Э/Эту is used for tangential derivative and d/dt for a derivative with respect to time.)

The velocity-potential time derivative is obtained by using the definition

Ф* = ±Jo(y/2)d/:

so that the local potential is the sum of the vortices from the leading edge to the )th vortex point. After substituting these terms into Eq. (13.123), the pressure difference between the airfoil’s upper and lower surfaces becomes

ДPi = p{V{i) + »w, W(t) + ww], • Ту~j + ~ 2) rfc) (13-12?)

For example, in the case of a translatory motion parallel to the X axis, as in Eq. (13.105), and for a flat thin airfoil placed on the x axis the normal and tangential vectors become

n, = (0, 1) ty = (1, 0)

The upper and lower tangential velocity components are then

Q, = u„±

and the pressure difference becomes

The total lift and moment are obtained by integrating the pressure difference along the chordline:

N

L’ = Fz = 2 Apy A/; cos ос, (13.128)

i=і

N

M0 = – 2 Apy cos о, Д/, дсу (13.129)

y=i

The drag of the two-dimensional airfoil during the unsteady motion is due to the induced angle caused by the wake and due to the added-mass effect (caused by the relative fluid acceleration, see also discussion leading to Eq.

(13.37) )

D’ = jt р("*Г, Г* Д4 sin a*) (13.130)

y=i 4 c,**=i >

Here the first term is due to the wake-induced downwash ww which in the lumped vortex element case is evaluated at the panel’s three-quarter chord point. The second term is due to the fluid acceleration (second term in Eq. (13.127)) and its center of pressure is assumed to act at the panel center.

Vortex wake rollup. Since the vortex wake is force-free, each vortex must move with the local stream velocity (Eq. (13.21a)). The local velocity is a result of the velocity components induced by the wake and airfoil, and is usually measured in the inertial frame of reference X, Z. To achieve the vortex wake rollup, at each time step the induced velocity (u, w), at each vortex wake point is calculated and then the vortex elements are moved by

(Ax, Az)i = (и, w), At (13.131)

The velocity induced at each wake vortex point is a combination of the airfoil Гу and wake Г* vortices and can be obtained by using the same influence routine (Eq. (11.2)):

N Nw

(u, iv), = 2 VOR2D (Гу, Ху., ZwXj, Zj) + 5) VOR2D (ГWk, xWi, zw., xWk, zWk) y=i *=i

(13.132)

In this simple scheme the velocity components of the current or previous (or a combination of) time step were used. But more refined techniques can be applied here that may improve the wake shape near the trailing edge.

Summary

The solution procedure is described schematically by the flow chart of Fig.

13.25. In principle at each time step the motion kinematics is calculated (Eqs. (13.103-104)), the location of the latest wake vortex is established and the RHS, vector is calculated. Then, during the first time step the influence coefficients appearing in Eq. (13.118) are calculated and the matrix equation is solved. At later time steps, the airfoil vortex distribution can be calculated by the momentary RHS; vector, using Eq. (13.122). Once the vortex distribution is obtained the pressures and loads are calculated, using Eqs. (13.128-1301. To conclude the time step, the wake vortex locations are updated due to the velocity induced by the flow field using Eq. (13.132).

As an example, consider the sudden acceleration of a flat plate (discussed in Section 13.7) placed along the x axis. The angle of attack a can be obtained by rotating the plate frame of reference by в = a relative to the direction of Q,_. For this case at t > 0 the velocity of the origin is (X0, Z0) = (-<2=°, 0) and Eq. (13.116) results in the following free-stream components

Since the normal vector to the flat plate is n = (0,1) the right-hand side (downwash) vector of Eq. (13.117) becomes

RHS, = —(Qoo cos a + uw, Qx sin a + ww) • (0, 1) = -((?» sin a + vvw)

The wake-induced downwash is obtained by using Eq. (13.113) and then at each moment Eq. (13.118) is solved for the airfoil vortex distribution. The pressure difference is then obtained by using Eq. (13.127). Results of this computation for the case of the sudden acceleration of a flat plate are presented in Figs. 13.8-13.10 along with the results obtained with the one-element lumped vortex method.

Program No. 13 in Appendix D includes most of the elements of this method except the matrix inversion phase and may be useful in developing a computer program based on this method. The matrix inversion is included, though, in Program No. 14 which is a more complicated three-dimensional model.

Solution of the vortex distribution for any given time t, at each x = const, station, is reduced now to solving this equation (Eq. (13.82)) for
y(_y) = y(yt t). Because of the similarity between this case and the steady-state slender-wing case we know that the spanwise circulation (load) distribution is elliptic, as in Eq. (8.74),

Г(у) = Гтах[і – ] (13.83)

and again Г(у) = Г(у, t). The spanwise vorticity distribution (shown in Fig.

8.18)

is obtained by differentiating with respect to у (as in Eq. (8.41) or (8.75)):

Substituting this into the integral equation, Eq. (13.82), results in

But this integral has already been evaluated in Chapter 8 (see Eq. (8.77)) and resulted in a constant: —лЬ(х)/2. Substituting this result into Eq. (13.85) yields

1 4Гтах Г-л;Ь(дО] _пдт] ■ dy

2л b(x)2 L 2 J Эх X dt

and after rearranging the terms:

(13.86)

which shows that the spanwise induced downwash due to an elliptic lift distribution is constant and independent of y. The value of Tmax (at each x station and time t) is easily evaluated now and is

Гmax=-b(x)W(x, t) (13.87)

Recalling that the velocity potential can be defined by a path of integration along the local у axis (for an x = const, section),

where the integration starts at the left leading edge of the x = const, station and the integration path is above (0+) or under (0—) the wing.

By substituting y(y) and Гтах into Eqs. (13.79)-(13.81), the crossflow

potential and its derivatives are obtained:

(13.88)

“(*> 3;, 0±, t) = ~(x, y, 0±, t) = |w(jc, ‘)д/[^г] (13.89)

This differentiation can be executed only if the wing planform shape b(x) and chordwise downwash W(x, t) are known. The other derivatives of the velocity potential are

and based on Eq. (13.77), the downwash on the wing is

Once the velocity field is obtained the pressure distribution on the wing can be calculated by using the Bernoulli equation. If the reduced frequencies and the accelerations in the z direction are small, then the pressure is given by Eq. (13.66a):

where the first term is similar to the steady-state (circulatory) term, and the second term is a result of the change of the flow with time. The pressure difference across the thin wing is then

Г ЭФ ЭФ) Э д

Ар =р(х, у, 0-) – р(х, у, 0+) = 2p[t/— + -^-J = pU— АФ + р — ДФ

(13.92)

(13.93)

The longitudinal wing loading is obtained by integrating the spanwise pressure

difference and by recalling the result of Eq. (8.88) that

(13.95)

This is the unsteady version of the slender wing lift (see Section 8.2.3) and a similar formulation was derived by Lighthill13 6 when studying the swimming of slender fish in small-amplitude motion.

In regard to the momentary drag force, recall that the x axis of the coordinate system used for this problem remains parallel to the flight path (Fig. 13.21) and the normal component of the force was designated as lift. Similarly it is possible to define the axial component of the force as drag. Since at each moment for a given set of boundary conditions the potential flow problem is independent of time (excluding the wake influence—which is neglected in this example) the drag component due to the circulatory part of the force can be approximated by the steady-state results of Eq. (8.95). Consequently, owing to the leading-edge suction, the drag due to the circulatory lift (the first term of Eq. (13.95)) is half of the projected pressure difference component:

npU Э

2~~4~~dx

while the drag due to the fluid acceleration (second term in Eq. (13.95)) is not reduced by the leading-edge suction. Thus, the instantaneous drag force becomes

Example 1. Heaving oscillations of a slender delta wing. As one of the simplest examples let us consider the small-amplitude heaving oscillations of a slender delta wing. The x, z coordinate system is selected such that it moves to the left of the page at a constant velocity U(t) = t/„. The wing remains parallel to the jt axis, but oscillates up and down at a frequency to and amplitude h0 (see Fig. 13.22). The small displacement of the wing relative to the x axis is then

i](x, t) = h0 sin Ш

FIGURE 13.22

Heaving oscillations of a wing with an amplitude of hn.

Then with v, ei = (0,0, dr}/dt) the time-dependent downwash W(x, t) of Eq. (13.55) becomes

IV(де, f) = — = h0w cos tot at

The longitudinal loading dL/dx is obtained from Eq. (13.95):

— = — л?^°° — 1/|0о)Ь(х)2 cos tot] —[кпшЬ(х)г cos tot] dx 4 dx 4 at

= – —[b(x)2]h0(o cos tot^bixfh^2 sin cot

For a flat triangular delta wing with a chord c and trailing edge span of bT E the span b(x) becomes

x

b(x) — b-r. E. ‘

c

and the longitudinal lift distribution is

dL дгр1У»/>|.Е. , , , лгр(ЬТЕх)2, 2 .

— —————- r – xhnto cos tot + ——— 5-— h0<o sin tot

dx 2 c2 4 c

The lift of the wing is obtained by integrating dLjdx along the chord:

L= [ ~dx = ~~^b2 Eh0(o cos ш + ch0<a2 sin ш (13.97)

Jo dx 4 12

The pitching moment about the apex (де = 0) is

M0 = — f ^X dx = – br. E ch0(o cos cot+ ^§b2T Bc2h0<o2 sin cot (13.98) J0 dx 6 10

As we can see the loads on the wing are created by two terms that have a phase shift between them. The lift of the wing, for example can be divided such that

npU„ 2 xpU„,2 .

—— b^E. hoW cos tot =———– — by E. T)

The time-dependent vertical displacement of the wing r(i) and these two terms of the lift (Lj, L2) are shown schematically in Fig. 13.23. The first term L, resembles the steady-state (circulatory) term and the lift is a result of the instantaneous effective angle of attack. This lags in phase with the motion such that when the wing moves downward it creates lift L, and vice versa. The second term L2 is a result of the wing vertical acceleration (added mass) and is in phase with the motion.

Example 2. Pitch oscillations of a slender delta wing. Another simple example is the small-amplitude pitching oscillations of a slender delta wing. The origin of the x, z coordinate system is now moving to the left of the page at a constant velocity U(t) = U„. The wing is pitching about the point xcg with a frequency (i) and a small amplitude в0 (shown in Fig. 13.24):

в = в0 sin o>t

The wing’s chordline rj(x, t) is given then by

tj(x, /) = — (x – jcC|?) tan 0 ~ —(x – xcg)6 = -(дг – xcg)e0 sin cot

By substituting the derivatives of r into Eq. (13.50) the time-dependent downwash W(x, t) is obtained:

W(x, t) = – f/,0o sin cot – в0со(х – xcg) cos cot

Note that the same result can be obtained from Eq. (13.55) by placing the chordline on the x axis and pitching the x-z plane with 0. The longitudinal loading dL/dx is obtained from Eq. (13.95):

dx 4 dx

For a flat triangular delta wing with a chord c and trailing edge span of bT E. the local span b(x) becomes

and the longitudinal lift distribution is

+ [U~.e<fi>x2 cos cot – в0(о2(х3 – x2xcg) sin to/]

The lift of the wing is obtained by integrating dL/dx along the chord:

The pitching moment about the apex (x = 0) is

As we can see in this case, too, the loads on the wing can be divided into three terms that have a phase shift between them. The three components can be

rewritten as:

ЭМо npU„bn /Зс4 2с3 лрЬтЕ r, c4 npU„L2 ( 2 2c дв ~ 4 7~ T~TXcs)~T~U~4~~ 4 Ьт EC 3 Xcg)

(13.101)

As the simplest example for the conversion of a three-dimensional wing theory to the time-dependent mode, consider the planar motion of a slender wing in the x, z plane.13 5 (This is a three degrees of freedom motion with X, Z, and в as shown in Fig. 13.21.) Since for a slender wing the longitudinal dimension is much larger than the other two dimensions (x »y, z) we can assume that the derivatives are inversely affected such that:

э_ д_ d_

Эх ^ dy ’ dz

As in the case of the steady-state flow over slender wings and bodies, substitution of this condition into the continuity equation (Eq. (13.12)) allows us to neglect the first term, compared to the other derivatives:

This suggests that the crossflow effect is dominant, and for any x = const, station, a local two-dimensional solution is sufficient. An interesting aspect of this simplification is that the wake influence is negligible, too, as long as the longitudinal time variations (e. g., wing’s forward acceleration) are small.

The slender, thin lifting surface with a chord length of c is shown schematically in Fig. 13.21. At t = 0 the wing is at rest in the inertial system X, Y, Z, and at f > 0 it moves along a time-dependent curved path, S (for this particular case S is assumed to be two-dimensional). For convenience, the coordinates x, z are selected such that the origin О is placed on the path 5, and the x coordinate is always tangent to the path. The wing shape (camberline) is given in this coordinate system by r/(x, t), which is considered to be small (jj/c«1, since small-disturbance motion is assumed). Also, the normal component of the kinematic velocity is small (e. g., ec/U(t) « 1).

The time-dependent version of the boundary condition requiring no normal flow across the surface (at any x = const, station) for this case is given by Eq. (13.50):

(13.77)

where Ф is the wing’s perturbation potential and the subscript В is not used for the slender wing case (since the wake effect is neglected).

and it is clear that y, U, rf, and в are functions of time. Comparing this form of the boundary condition with the formulation for high-aspect-ratio wings (Eq. (8.11)) clearly indicates that due to the slender-wing assumption the effect of the vortex lines parallel to the у axis (including the time-dependent portion of the wake) were neglected.

The fluid dynamic pressures and loads generated by the airfoil can be calculated by using the unsteady Bernoulli equation:

By recalling Eqs. (13.45)-(13.46), the terms inside the second parentheses become – (V0 + ft X r) = [17(f) – 0(0*1, 0, 0(0*]

If the reduced frequency (вс/U) is small and the point of interest is on the airfoil then the pressure equation becomes

n — n Дф ЭФ ЭФ

and in this case ЭФ/dz has the same value above and under the chordline and does not contribute to the pressure difference. The pressure equation thus becomes

p ox dt

where the first term is similar to the steady-state term, but for the time – dependent case also the change in the potential contributes to the pressures (due to the acceleration of the fluid).

The pressure difference across the airfoil Ap (positive Ap is in the +z direction) is then

Г ЭФ ЭФП г э э 1 Ар=Р’-Ри = 2р[Щі) — + — ] = р|и(0-ЛФ+-Дф] (13.67)

where ДФ(х, г) = Ф(х, 0+, г) – Ф(х, 0—, г) = /5 у(*0. 0 dx0 = Г(дс, г) and there­fore the pressure difference in terms of the airfoil chordwise circulation у becomes

Ap = plJ(t)Y(x, 0 + J-J y(*o.0«ko] (13.67a)

The normal force on the thin airfoil is then

L’ = Fz = j Apdx = j р[и(і)у(х, t) + pjT(x, t)] dx

= pU(t)T(t) + p f f Г(х, 0 dx Jo dt

where the first term is due to the instantaneous circulation (and similar to the steady-state circulatory term) and the second term includes the contribution of the time dependency.

To evaluate the time derivative of the velocity potential in terms of the coefficients A„ (appearing in Eq. (13.59)), recall that Ф = J q • d/ (or ДФ = Sydl)

д Э Г д Ґ c

— ДФ(х, t) = — j y(xo,

For a given airfoil geometry, the mean camberline t](x, t) is a known function and the coefficients A0(t), Ax(t), A2{t), …. can be computed by Eqs. (13.61) and (13.62) (assuming that the wake influence is known). The pressure difference of Eq. (13.67) can be evaluated since all terms in this equation depend on the coefficients An(t). The force in the z direction is then

L’ = FZ= Apdx = 2p[ (е/2(0|А)(0 1 + C-°^+ X /f„(0sin(/tfl)1 Jo Jo 1 L sin v n=i J

+ Bo(& + sin i?) + sin 20^

і B„[^” ~ Ц* – ™(M 4-1) £si" ***

“2 L 2(n — 1) 2(n + 1) JJ 2

These integrals are similar to those treated in Section 5.3 and after their

evaluation we get

L'(t) = pcj— 80 + —81 + ~ 82 + 7tU2A0 + — U2A^ (13.71)

In terms of the i4„’s (using Eq. (13.70)), we get

L’(t) = ;rpc{[t/2A0 + — (IMo)] + ‘у + (UA0 + ljt (UA*

(13.71a)

and it is clear that the velocity and the coefficients are a function of time (e. g., U = 1/(0, Ая-Ая(0).

The pitching moment about the airfoil’s leading edge is M0(t) = – jCbpxdx = – j P^WyАФ + |дф]*<*х

= -2p f f І/2(г)[ло(0 + 2 A„(0 sin (n0)l

Jo l L sin v ns) J

/0 1

+ 8O(0 + sin 0) + ~ 4 s^n

■A rsin(n-l)0 sin(n +1)011 c c. .

and after an evaluation of these integrals we get

With the use of Eq. (13.71) for L’, the moment about x = 0 becomes,

Щ0 = “PC2 | [‘у (a> + Al “ у) + ^ B° + Bl + 4 B2 ~ Jg 83] (13.72) and in terms of the A„’s,

Щ0 = – pc2f [y A, + II(£M0) + у Л, + f |(^0

и2 с д с d 1

“ +–(ВД-—-(iM3)

4 8 3r 2′ 32 ЭГ 3/J

Example 1. Small-amplitude oscillations of a thin airfoil. One of the simplest and yet important examples is the small-amplitude unsteady motion of a flat plate airfoil, which was analyzed by Theodoresen13 3 and by von Karman and Sears.13 4 For this case let us assume that the (x, z) frame of reference in Fig. 13.14 moves to the left of the page at a constant speed U(t) =U = const, in an otherwise stationary fluid. Also, the (x, z) frame does not rotate for this example

(0 = 0 = 0) and the small-amplitude unsteady motion will be introduced through the vrcl term (or the dr}/dt term in Eq. (13.55)) in the boundary conditions.

The time-dependent chordline position can be represented by a vertical displacement h(t) (positive in the z direction) and by an instantaneous angle of attack a(t) (Fig. 13.18). The chordline shape is then

t) = h — a(x — a)

where a is the pitching axis location. For simplicity, first, we shall assume that the pitching axis is at the origin (a = 0, and h is the vertical displacement of the leading edge) and r then becomes

tf = h — ax

and the vertical displacement is small (e. g., r « c). The derivatives of t] are

in

dx

where the dot denotes a time derivative. Substituting this into the downwash W(x, t) term of Eq. (13.55) we get

ЭФуу

W(x, t) = – Ua + h — ax —-—

Since the wake effect is a function of the motion history let us concentrate first on the loads due to the motion only. This portion of the downwash, W*(x, t), is then

• . c c

W*(x, t) = – Ua + h — da = —Ua + h — – dr + – dr cos &

and here x was replaced by the trigonometric variable d, using Eq. (13.58). Substituting this term of the downwash into Eqs. (13.61) and (13.62) provides the

values of the An coefficients:

ac ‘2U

The circulation due to the downwash W* can be obtained by recalling the results of Eq. (5.58):

г*(0 = I У(*. 0 dx = jicu(a0 + y)

and after substitution of the A„ coefficients the circulation becomes

Г*(г) = nc(Ua – h + led)

The lift per unit span is then obtained from Eq. (13.71a):

and in terms of the displacement A and the angle of attack a, we get

L* = npUc(^Ua – h + ^ cdj + -^pc2^ (Ua – h) + ^ dj

In the derivation of this expression the downwash of the unsteady wake was not included. Theodorsen13 3 and von Karman and Sears134 showed that for a small-amplitude oscillatory motion the final result will include similar terms and the effect of the wake is to reduce the lift due to the first term in L* by a factor of C(k), which is called the lift deficiency factor. Now, if we consider the harmonic heave and pitch oscillations such that

h=ha sin (ot a = a0 sin Ш

then the lift per unit span becomes

L’ = KpUcC(k)i^Ua — h + ^cd^ + —(ua — h + ^oc^ (13.73a)

In the case when the pitch axis is moved to a location a (and also h is measured at this point), as shown in Fig. 13.18, then the lift per unit span will have a form similar to the result of Ref. 13.3:

L’ = jzpUcC(k^Ua — h + ^-^cdj +—~ j^(f/d — A) + -^dj (13.73b)

This can be rewritten as

L’ = L[ + L’2

where L[ is similar to the circulatory lift term in a steady motion and L is the lift
due to the acceleration, called the added mass. The lift deficiency factor C(k) is plotted in Fig. 13.19 versus the reduced frequency k, which is defined similarly to the nondimensional number of Eq. (1.52):

we 2U

As Fig. 13.19 indicates, the wake has a delaying effect on the circulatory part of the lift such that

L’x(t) = L sin (cot — w)

and m represents the time shift effect of the wake (note that m changes with the reduced frequency as shown in Fig. 13.19).

After a similar treatment of the moment about the leading edge we get

M0=-^[-^h+^a + ^c2a + UC(k)^-h + Ua+jdr]} (13.74a)

Again, when the pitch axis is moved to point a (Fig. 13.18) then the pitching moment about this point is

4a 4a

+ ——- jd

c c /

– (j – l) t/C(*)[-A + Ua + c( – *)<*]} (13.74A)

The most important portion of Theodorsen’s analysis is that the basic nature of the unsteady effect can be briefly summarized by the C(k) diagrams in Fig. 13.19. As the reduced frequency к increases, the magnitude of the pUT term in the lift is reduced. Additionally, the lift lag initially increases with the reduced frequency, but for к > 0.4 a gradual decrease in the phase shift is shown.

The above model for the small-amplitude oscillation of a thin airfoil is useful in estimating the unsteady loads in cases such as wing flutter or propulsion. The propulsion effect due to the heaving oscillations of a flat plate is shown schematically in Fig. 13.20. Recall that due to the leading edge suction the circulatory part of the lift (pUT) is normal to the instantaneous motion path, which clearly results in a propulsive (forward-pointing) com­ponent. If the heaving motion is relatively slow then the second term in Eq.

(13.73) is relatively small, too, and high propulsive efficiencies131 can be obtained.