Computational Mechanics (1989) 4, 231-243
Computational Mechanics © Springer-Verlag 1989
Boundary integral equation for wave equation with moving boundary and applications to compressible potential aerodynamics of airplanes and helicopters L. Morino Universitfi "La Sapienza", Dipartmento di Meccanica e Aeronautica, 00184 Rome, Italy B.K. Bharadvaj Douglas Aircraft Co., Long Beach, CA 90446, USA M.I. Freedman Boston University, Dept. of Mathematics, Boston/MA 02215, USA K. Tseng United TechnologiesResearch Center, East Hartford/CT 06108, USA
Abstract. This work presents a general boundary-integral-equation methodology for the solution of the wave equation around objects moving in arbitrary motion, with applications to compressible potential aerodynamics of airplanes and helicopters. The paper includes the derivation of the boundary integral equation for the wave equation, in a frame of reference moving in arbitrary motion (in particular, in translation and in rotation). The formulation is then applied to study unsteady potential compressible aerodynamic flows around streamlined bodies, such as airplanes and helicopters. The formulation is given in terms of the velocitypotential, for which an explicit treatment of the wake is required; a discussion of the formulation for the wake transport is included. The advantages of the velocity-potential formulation over the acceleration-potential formulation are discussed. The boundary-elementalgorithm used for the computational implementation is briefly outlined. Validation of the formulation is presented for airplane wings and helicopter rotors in hover. The test cases fall into two categories, prescribed-wake and free-wake analyses. The validation of the prescribed-wake analysis is presented for compressible flows, subsonic for helicopter rotors, transonic for airplanes. The numerical validation of the free-wake analysis of helicopter rotors is presented for incompressible flows. 1 Introduction The objective of this paper is to present a new boundary integral equation for the solution of the wave equation around an object moving in arbitrary motion, with applications to potential aerodynamics. The integral equation is expressed in a frame of reference that also moves in arbitrary motion. The frame of reference is typically chosen so as to minimize the motion of the surface of the body with respect to the frame of reference itself; in particular, it is rigidly connected with the body when this moves in rigid-body motion (in practical applications, the body moves with small oscillations around a rigid-body configuration). The frame of reference which moves with the body will be referred to as the 'body frame of reference', to distinguish it from the frame of reference connected with the medium. The formulation presented here is important in aerodynamics and acoustics, but, in general, is applicable to all the physical phenomena governed by the wave equation in any medium. In view of the fact that, for all the applications in this paper, the medium is air, we will refer to the fixed frame of reference as the 'air frame of reference'. The methodology is applied to the analysis of unsteady compressible potential aerodynamics of lifting objects, such as airplanes, helicopters, and propellers. In this case the problem is complicated by the presence of the wake, i.e., a layer of vortices which, in potential aerodynamics, manifests itself as a discontinuity in the velocity potential. F o r airplanes, the wake geometry is typically assigned ('prescribed-wake analysis'), whereas for helicopters the wake geometry itself is obtained as part of the solution process ('free-wake analysis'). The development of an efficient methodology for the analysis of unsteady viscous compressible flows around complete aircraft configurations
232
Computational Mechanics4 0989)
having arbitrary motion is the long-range objective of the field of computational aerodynamics. Such a development would represent a maj or improvement in the methodology available for aircraft design. The formulation presented in this paper should be viewed as a contribution in that direction. The effects of compressibility are addressed here primarily for the case of inviscid (potential) flows; a viscous-flow correction is used in connection with transonic flows. A brief review of recent developments in potential aerodynamics is presented here. Some related developments in acoustics are also included. The equation governing compressible potential flows is a nonlinear wave equation. A general integral equation for potential flows around objects moving in arbitrary motion (with respect to a frame of reference in uniform translation) was presented in Morino (1974). Numerical implementations of the formulation for isolated wings were presented in Morino and Kuo (1974) and Morino et al. (1975), for subsonic and supersonic, steady and unsteady flows, and in Morino and Tseng (1982) and Tseng (1983, 1984a, 1984b), for transonic flows. The use of a frame of reference in uniform translation is acceptable in airplane aerodynamics. However, in the case of propellers and helicopter rotors, it is more convenient to formulate the integral equation in a frame of reference that moves with the rotor. A general boundary integral equation for the nonlinear wave equation for applications to potential aerodynamics was developed by Morino et al. (1985), who assumed that the frame of reference moves in arbitrary motion; the surface of the body itself is then allowed to have arbitrary motion with respect to the frame of reference. Some mathematical aspects of the formulation are based on the work of Ffowcs Williams and Hawkings (1969), who developed a related integral formulation for acoustics. Their formulation, however, differs from that of Morino et al. (1985) in several respects. First, it is strictly related to the physical problem of acoustic waves and is not directly applicable to the mathematical problem. Also, in their formulation the surface is rigidly connected with the frame of reference (the implications of this assumption are discussed at the end of Appendix A). However, the main difference is that their final result is not an integral equation but an integral representation that, given certain quantities in the near field, allows the evaluation of the pressure distribution in the far field. More recently, boundary integral equations related to that by Ffowcs Williams and Hawkings (1969), have been developed and used in aerodynamics of rotors, for instance by Hanson (1983), Long (1983), Farassat (1984), Tai and Runyan (1984), and Williams and Hwang (1986). A stricly mathematical formulation, very similar to that of Morino et al. (1985), was developed recently by Farassat and Myers (1987), which provides several expressions for the integral representation. It should be noted that the above-referenced works by Hanson (1983), Long (1983), Farassat (1984), Tai and Runyan (1984), and Williams and Hwang (1986), are based on the use of the acceleration potential, which has typically been used in connection with integral-equation formulations for rotor aerodynamics. On the other hand, the velocity potential has typically been used in connection with finite-difference solutions of rotor aerodynamics (for a review, see Caradonna 1985; for recent developments on the free-wake analysis, see Steinhoff and Ramachandran 1986). We believe that the velocity-potential formulation is more advantageous than the accelerationpotential formulation for the following reasons: (i) the availability of exact explicit expressions for the nonlinear terms in the differential equation; (ii) the simplicity of the implementation of the boundary condition on the surface of the blades; and (iii) the explicit treatment of the wake. It should be emphasized that, indeed, the last item is what makes the acceleration potential (which is continuous across the wake) more attractive than the velocity potential for lightly loaded rotors, where helicoidal wakes may be used. However, for helicopter rotors in hover, the loads are strongly affected by the geometry of the wake used in the calculations, and therefore, the use of an accurate wake geometry is essential to the solution of the problem. A methodology for the analysis of the wake motion for helicopter rotors was presented by Morino et al. (1985) and refined by Morino and Bharadvaj (1985), using an approach that is closely related to that introduced by Suciu and Morino (1977) for isolated wings. These papers are limited to incompressible flows. The formulations by Morino et al. (1985) and Morino and Bharadvaj (1985) were unified and further developed by Morino et al. (1987a and 1987b), who presented a general formulation for potential compressible flows, with a frame of reference in arbitrary motion, which includes the analysis of the wake transport. Their integral equation for rotors is simpler to use than that by Morino et al. (1985) and is remarkably similar to that developed by Morino (1974) for airplane aerodynamics.
L. Marina et al.: Boundary integral equation for wave equation with moving boundary
233
This paper presents an outline of the formulation and numerical results which demonstrate its feasibility as a computational methodology.
2 Boundary integral equation for wave equation The derivation of a boundary integral formulation for the nonlinear wave equation which governs the motion of the air around a body that moves in arbitrary motion, is presented here, using a body frame of reference (defined above); details are available in Marina, Freedman and Tseng (1987b). Following Ffowcs Williams and Hawkings (1969), the derivation of the boundary integral formulation is facilitated if the problem is first formulated in the air frame of reference and then the integral equation is expressed in the body frame of reference. We will indicate with 4 the location of a point in the air frame of reference and with ~ the corresponding time (x and t will be used for the body frame of reference). In the air frame of reference, the wave equation is V~4
1 a2~b a~ av2 - Z
~ outside a(t)
(1)
where a indicates the boundary surface and Z represents the nonlinear terms (which in aerodynamic applications are important in the transonic regime). The boundary condition at infinity is q5 = 0. As pointed out in Marina et al. (1987 a), the derivation of the integral formulation is facilitated if the boundary-value problem is replaced with an equivalent infinite-space problem. This may be accomplished by introducing the 'domain function', E(4, z) = 1 outside a, = 0 otherwise
r> 0
(2)
and the function ~ = EqS. Combining with Eq. (1) one obtains 1
a~ ~'c2 - )~
(3)
where )? is defined as 2(~,t)=Ez+V~
E'V¢¢+V~'(OV~E)
a~[a.c
ar + a r k
a~]_]'
(4)
The solution to Eq. (3) is CO
(5)
(4,, r,) = f i l l G (~ - 4,, ~ - %) Z (G ~) d Vdr, --oO
with G being the fundamental solution for the wave equation, -1 a ( 4 - 4,, z - %) = - - g ) ( r 4rc¢
(6)
- z, + o/aco),
where ~o= [4-4,1 and 6 is the Dirac delta function. Using the expression for 2 given by Eq. (4), integrating by parts and noting that V¢ G = - V~, G (where V¢, indicates the gradient with respect to 4,) and aG/ar = - aG/~z,, one obtains co co / ~ , [ aq~ aE E(4*' r*) O (¢*' r*) = ~I~ GV'O ' V~-EdVdr + V¢* " ~ f G ~ Vc E d V d r - a~ O7 a a v a-~ d V d r --COO
--cO
1 a2
~
oo_.
az~fl~ _ G f 3 aE a~
co + ~E IfGzd~Vdr_
--
(7)
This is an integral representation of the solution of Eq. 1. Note that this expression contains two Dirac delta functions: one in the fundamental solution, G, and the other arising from V E [Eq. (2)]. Thus, the four-dimensional integrals may be reduced to two-dimensional ones. However, before
234
C o m p u t a t i o n a l Mechanics 4 (1989)
doing this, it is convenient to express the integrals in the body frame of reference. Let ~ = ~ (x, t), -c = t indicate the transformation relating the two frames of reference and x = x (~, -c), t = ~ indicate the inverse transformation• Note that the Jacobian of the spatial-temporal transformation is equal to one. Furthermore, indicating with V the gradient with respect to x, we have Vf= V J . On the other hand, = 0f +
• V=
(8)
0¢
where vR = ax/Oz -- - a~/at is the velocity of a point ~ of the fixed frame of reference as seen in the moving frame of reference (equal and opposite to the velocity of a point x of the moving frame of reference as seen in the fixed frame of reference). Combining with Eq. (7), introducing the function g(x, x,, t, t,) = t - t, + Q/aoo[where 0 = ]el, with e = ~ (x,t)-
~ (x,,t,)] and noting that for any two functions f(t) and g(t) ~ f(t)3[g(t)] dt = ~ ] g (t)]-lf(ti) (where the ti'S a r e the roots of g = 0, whereas g = ag/Ot), one obtains --O0
-1
~[4@~d
dq~ dE] 1 d ~[-1 dt dt d V - a~ dt-, . 4~
1 ~/[ -1 a~_ . ~?
ti
q~VE] dV
dE] ~[E-I q~dV+ • ti
] 4~Z
dV
(9)
~ ti
where [ ], indicates evaluation at t -- tj and V. indicates the gradient with respect to x. whereas ~ is defined as ~ =p[g[ -- el (1 - MR)I, where M R = e" vR/9 ao~with vR defined above. A more convenient expression for 0 may be obtained by introducing the function ?(x,, x, t,, t) = a~ (t
-
t , ) 2 - ~2
(10)
Note that 7(0 = 0(t < t,) yields the same roots as g(t) = 0. In addition, 1
2aoo
~87 t=t,-o~/%-
1
a2 ( t -
a~
t,) + e"
O~t t = t , - a / % - - o~ 1 - -e. VaZ=
(11)
Equation (9) may be simplified further by taking advantage of the fact that the derivatives of E are equal to zero everywhere except on the surface o-. For simplicity, assume that the surface ~ moves in arbitrarily-prescribed rigid-body motion [the general case is examined in Morino et al. (1985) and is outlined in Appendix B]. Under this assumption, if the moving frame of reference is chosen to be rigidly connected with the surface o-, the function E, defined by Eq. (2), becomes time independent (for t > 0); in this case [Eq. (B.1)]
~ f V E d V = j~ f ndcr
(12)
--oo
Also, assuming f = 0 at t = 0 in order to avoid the initial-condition contribution arising from OE/~t (see Morino 1974 for an analysis of this contribution), one obtains (Eq. B. 1) oo dE
_~f ~t dV= - [,Sofvoda
(13)
where v~ = vR'n is the velocity of the surface 0. Combining Eqs. (9), (12) and (13) and using Eq. (8) to express dq)/dt, one obtains
[-1
84)
8(-1)
- 1 049{Oti
v~)l dcr
L. Morino et al.: Boundary integral equation for wave equation with moving boundary
235
where ~h
~n
+
1 ;~a+-V+VR"V
1 ©~
-
n • V, -
_5- v~
ao~
(15) d
(16)
dt,
Equation (14) is the desired integral representation of q5(x,, t,) in terms of the values of q5 a n d ~(a/~n on a; if x, approaches a point of o-, Eq. (14) yields an integral equation that may be used to solve the problem. It is apparent that if the velocity of the body frame of reference equals zero, Eq. (14) reduces to the well known Kirchhoff formula for the wave equation.
3 A p p l i c a t i o n s to a i r p l a n e s a n d h e l i c o p t e r r o t o r s
In this Section the results of Sect. 2 are used to study the aerodynamic flow field around airplanes and helicopters. The flow field is assumed to be isentropic, nonconducting and initially irrotational. Such a field remains irrotational at all times and therefore the velocity may be expressed in terms of the velocity potential, ~b, as v -- VqS. The differential equation for the potential, ~b, is given by Eq. (l) [for an explicit expression for the nonlinear terms, Z, see for instance Morino (1985); as mentioned above, in aerodynamic applications these terms are important in the transonic regime]. The surface a is assumed impermeable; hence, the boundary condition on a point of a is (v - vb) • n = 0, or ©(o/On = vb • n is the velocity of the point on the surface of the body a. The boundary condition at infinity is q5 = 0. As mentioned above, in the case of lifting bodies, the problem is complicated by the presence of the wake. The treatment of the wake for compressible flows is outlined in Appendix A, with emphasis on the free-wake analysis. The main result of the analysis of Appendix A is that the wake manifests itself as a surface of discontinuity on the potential qS, which is equivalent to a vortex layer. The boundary condition on the wake is that the potential discontinuity, AqS, remains constant following a wake point. In transonic and supersonic flows an additional complication arises because of the presence of shock waves; for a discussion on the analysis of shock waves, in particular on the shock-capturing nature of the algorithm used for transonic computations, the reader is referred to Morino and Tseng (1982) and Tseng (1983, 1984a, 1984b). As mentioned above, in airplane aerodynamics, one may assume that the body frame of reference is in uniform translation, with velocity v b = - U~ i (if the motion of the airplane includes small oscillations about a reference configuration, these are typically accounted for only in the boundary conditions, with the integral equation written in a frame of reference rigidly connected with the reference configuration). Thus, the transformation relating the frames of reference is ~ = x - Uoo ti, z -- t. This yields, for the function 7 defined by Eq. (10), ? = a~ 0 2 - [r 2 + 2 Uoo (x - x,) 0 + U~ 0 2]
(17)
with 0 = t, - t and r = Ix - x, [. For the sake of simplicity, only the case Moo = U®/aoo < 1 is considered here [this includes subsonic and transonic flows; for the supersonic case, see Morino (1974), Morino et al. (1975) and Tseng and Freeman (1985)]. Then the equation ? = 0 (t < t,) yields only one root 1
0-
aoo
B2 [rp + M + (x - x,)]
where/3 = ]/1 - M 2 r++ =
(18)
and
( x - x , ) +2+ p+r+ =
- x+) + + p + ( y - y + ) + ÷
- z,)+
(19)
C o m p u t a t i o n a l M e c h a n i c s 4 (1989)
236
In addition, using Eq. (11), one obtains ~ = D" Next, consider the integral representation for qS, Eq. d (14). The above transformation implies ve = U~ i, dt - St + U~ , and v~ = vb .n = - U®nx. Combining with Eq. (14) one obtains
E(x,)(p(x,,t,) = ~o 4 - D OR
OR ~
Z
+ 4-~rp Ot 8h - 2
a~ / b , - o
dV
d~r (20)
t, --0
with ~/~ h = ~/~n - M ~ nx ~/Sx and 0 and r~ given by Eqs. (18) and (19) respectively, in agreement with the results of Morino (1974). Equation (20) is the desired integral representation (for airplanes) of q5(x,, t,) in terms of the values of ~b and ©O/Onon o; if x, approaches a point of o-, Eq. (20) yields an integral equation that may be used to solve the problem. Note that the volume integral is relevant only for transonic flows. The wake is dealt with as indicated in Appendix A. Next, consider a helicopter rotor in hover (i. e., in uniform rotation). Let us choose the z-axis in the direction of the angular velocity of the rotor, e) = ok. In this case the transformation is given by ~ = U(t) x, where U(t) is a rigid-body rotation matrix given by
U(t) =
coscot
- sincot
sinco t
cosco t
0
0
0
1
(0
0) = U r ( - t)
(21)
It is easy to verify that U(fi + t2) = U(tt)U(t2) and U = eat = ~
a =
co
0
0
0
0
0
1
A~t ", where
(22)
The function 7, defined by Eq. (10), is now given by 7 = a% (t
-
t,) 2 -
[U(t)x
-
U(t,)x,I
2 = a 2 0 2 q-- 2 x T [ U ( O )
-- I]x,
-- r 2
(23)
where 0 = t , - t and r = Ix - x,J. Note that the roots 0i of the equation 7 = 0 are independent of t,; this is to be expected since, from a physical point of view, the formulation for uniform rotation is invariant to time-shifts. Also, according to Eq. (11),
= aoo 0 + xrAU(O)x,/ao~.
(24)
Next, consider the integral representation. Note that in the body frame vR = - a ~ M , where M = e) x x/a~ is the 'Mach vector' of the point x as seen in the air frame of reference. Hence, using Eq. (8), one obtains d/dt = O/Ot- a~M" V. Also vo = ao~M" n. This implies [Eq. (15)]
a/Oh = O/On - M . n M . V.
(25)
Next, note that ~ and 0 are independent of t, and that, for any f u n c t i o n f t h a t is independent of t,, Of/Off = ~f/Oh with [Eq. (16)]
O/Oh = - n" V, + M . n 21//," V,
(26)
(where M, = 09 x x,/aoo defines the Mach vector of the point x,, as seen from the air frame of reference). For the sake of simplicity, assume that the tip Mach number M r is less than one (i. e., that the flow is subsonic or transonic) and that V is completely inside the 'sonic cylinder', defined by IMI < 1; then the equation 7 = 0 (t < t,) yields only one root. Combining the above equations with Eq. (14) one obtains
L. Morino et al.: Boundary integral equation for wave equation with moving boundary
E(x,)O(x,,t,) = ~
-1 4~-~ ~h
-1
E
+ ~f~v
~h
X
- 1 a4 ~ + 4-~-~ ~t ~n + 2
(
dV
aoo
237
t, --0
da (27)
t, --0
with ~/~fi, ~/~h, ~, and 0 defined above. Equation (27) is the desired integral representation (for helicopter rotors in hover) of ~b(x,, t,) in terms of the values of q5 and ~gp/On on o; if x, approaches a point of o-, it yields an integral equation that may be used to solve the problem. Note that the volume integral is relevant only in the tip region and only when the tip region is transonic. The wake is dealt with as indicated an Appendix A. Considerable simplification in the expression for 0 and 0 is obtained if the terms of order 03 in U (0) are neglected (this is acceptable if the tip Mach number M r is relatively small, say M r < 0.6). This yields U(0) = I + A0 + A202/2and (noting that A = - A r a n d A x = o~ x x)
xr(U(O) - I)x,/a~ = r " MAO _ _1M . M,O 2 2
(28)
(where MA = ( M + M,)/2 is the arithmetic average of the Mach vectors M a n d M, at the points x and x, of the rotor frame of reference, as seen in the air frame of reference). Combining with the expression for 7 given above yields 7 = (1 - M . m , ) a~ 02 + 2 r ' M Aa~ 0 - r 2.
(29)
Solving the equation 7 = 0 (t > 0) yields
0 = [-- r. MA + ]/@" MA) 2 + (1 -- M" M,)r2]/aoo (1 - M " M,).
(30)
Also, according to Eq. (11),
O = ~(r-. MA)2+ (1 - - M . M , ) r 2.
(31)
Note the similarity between the above equations for q~, Eq. (27), for 0, Eq. (30), and for d, Eq. (31), and the corresponding ones for airplanes, i.e., Eqs. (20), (18) and (19), respectively. 4 Numerical results
The integral equations [Eq. (20) for airplanes, and Eq. (27) for helicopter rotors] are discretized by dividing the boundary surface into quadrilateral elements and assuming the u n k n o w n ~b to be constant within each element. In other words, the formulation used may be classified as zerothorder direct-approach boundary element method. For Eq. (27), the evaluation of the source and doublet integrals is facilitated by introducing a spacial transformation similar to the Prandtl-Glauert transformation for airplanes and neglecting certain terms in the expression for the doublets (see Morino et al. 1987b for details). As mentioned above, the formulation has been validated for several cases. For the case of airplanes, the formulation is the same as that of Morino (1974), for which extensive results have been presented in the works referenced above. Results for the free-wake analysis of helicopter rotors hovering in incompressible flows are available in Morino et al. (1985). Only new results are presented here, with the emphasis on helicopter rotors in compressible flows (Figs. 1 and 2). In order to put these results into proper perspective, recent results for helicopter rotors (incompressible free-wake analysis, Figs. 3 and 4) and for airplane wings (transonic analysis, Figs. 5 and 6) are also included. The hovering rotor considered in Figs. 1 and 2 is the same as that considered in the experimental and analytical studies performed by Caradonna and Tung (1981): the rotor has a N A C A 0012 section, radius of 45 in, constant chord of 7.5 in, and no twist. The case considered here is that with collective pitch of 8 ° and rotational speed of 1250 rpm, which corresponds to a tip Mach number of M r = 0.439. The wake geometry is prescribed to be close to that measured experimentally by Caradonna and Tung (1981) (see also Fig. 3). The parameters describing the numerical discretization
Computational Mechanics 4 (1989)
238 0.32
0.24 4 3
I0.16
7
2
1
0.08 Fig. 1. Convergence Analy-
0 0
0.2
0.4
0.6
0.8
1.0
r/R
sis: I N 1 = 4, 3 N 1 = 6, 4 N 1 = 7
2 N~ = 5,
0.32
0.24
I0.16
1
',S
/
2
•
4
0.08
3
•
0
J
0.2
2
0.4
0.6 r/R.
=
0.8
1.0
Fig. 2. Effect of the number of spirals; I N s = 3, 2 N s = 4, 3Ns=5, 4Ns=6, • comparison with experiments
are: the number of blade elements in the chordwise direction, N1, the number of blade elements in the radial direction, N2 (equal to the number of wake elements in the radial direction), the number of wake elements per spiral, in the azimuthal direction, N3, and the number of wake spirals, Ns. A convergence analysis (performed using three spirals) is shown in Fig. 1, which depicts the sectional lift coefficient as a function of the nondimensional radius, for N1 = 4, 5, 6, and 7 (with N 2 = 3N1 and N3 = 2N1). The results for N1 = 6 may be considered close to the converged three-spiral solution. The effect of the number of spirals is presented in Fig. 2, which depicts the sectional lift coefficient (equal to the sectional lift divided by the product of the tip dynamic pressure times the chord) as a function of the nondimensional radius, for N~ = 3, 4, 5, and 6 (with N1 = 6, N2 = 18, and N3 = 12). One may notice that the effect of the number of spirals is not strong. The solution for six spirals appears close to being converged. Also shown in Fig. 2 are the experimental results of Caradonna and Tung (1981). The agreement is good. The difference is attributed to the effect of viscosity and to the model used for the wake geometry (see next paragraph). Next consider the free-wake analysis which is performed following the procedure outlined in Appendix A. Here the validation for the free-wake analysis is limited to incompressible flows around helicopter rotors in hover, in which case the geometry of the wake has a strong influence on the load distribution (for airplanes see Suciu and Morino 1985). The free-wake results are presented in Figs. 3 and 4. The discretization corresponds to NI = 3, N2 = 7 and N3 = 12. Five spirals are used:
L. Morino et al.: Boundary integral equation for wave equation with moving boundary
0 ------~* ~.
i
-
i
239
r
0.6
&
-0.2.4
3
A
0.4
I-0.48I
i 0.2
I -0.96! -
1.20 0
[]B
0 i
-0.2
0.2
3
0.4
0.6
0.8
-0.4
1.0
r/R -
0'.2
5
0..40t
0'.4 x/c
o'.6
oL8
1.0
---
2.0. 1.5.
0.32!
~
Inviscid
~iscous
1.0
I 0.24I G
0.16i_
a~
0
'
"
1,
-0.5-4
0.08)
-1.0 [ O r
0
4
0.2
0.4 r/R
0.6 =
0.8
-1.5
1,0 6
o
~, Experiment
0'.2
o2
o'.8
o18
1.o
x/c
Figs.3-6. 3 Wake geometry at ~p = 0 : 1 prescribed-wake, 2 free-wake, • experiments. 4 Comparison. with experiments: 1 prescribed-wake, 2 free-wake, • experiments. 5 Comparison with finite-difference results ( - ) , 1 t / = 0.1, 2 t / = 0.5, 3 q = 0.9. 6 Comparison with experiments: present vs experimental results, 1 upper, 2 lower
the velocities of the wake nodes are evaluated only for the first one [the last four spirals move according to a prescribed velocity obtained from the experimental results; see Morino and Bharadvaj (1985) for details]. Figure 3 depicts the wake geometry as obtained from the free-wake analysis. Also shown is the wake used for the results of Figs. 1 and 2, which connects (in simple fashion) the vortex sheet and the tip vortex available from Caradonna and Tung (1981). It may be seen that the wake obtained from the free-wake analysis is in good agreement with the experimental results, and differs from the prescribed wake primarily in the region that connects the vortex sheet with the tip vortex. The difference in the wake geometry is the only reason for the discrepancy between the two curves in Fig. 4 which depicts the sectional lift coefficient as a function of the nondimensional radius. It is apparent that, as mentioned above, the use of the more rolled-up geometry obtained from the free-wake analysis tends to improve the agreement with the experimental results (increasing the sectional lift coefficient in the tip region and decreasing it near the root). Next, consider the use of the volume integral to account for the effects of the nonlinear terms for transonic flows. As the Mach number increases, the nonlinear terms in the equation become important (in the case of rotors this occurs only in the region near the tip, whereas for airplanes this occurs in the whole region around the aircraft). The validation of the formulation for transonic
240
Computational Mechanics4 (1989)
flows has been performed only for airplane wings. Some results obtained by Tseng (1983b) are presented in Figs. 5 and 6; additional results are available in Tseng (1983 a, 1983 b, 1984). The results shown in Fig. 5 are for a rectangular wing, at M~ = 0.908 and zero angle of attack. The wing has a 6% thick biconvex section, and an aspect ratio of four. The discretization corresponds to NI = 20 and N2 = 10. The pressure distributions at ~/= 2y/b = 0.1, 0.5, 0.9 (where b indicates the span of the wing) are compared with those obtained in Baley and Steger (1972) (at t/= 2y/b = 0.0, 0.5, 1.0) using a finite-difference solution of the same equations. The agreement is very good. Next consider a comparison with experimental results. [In the case of transonic flows the viscous effects are important; therefore, any comparison with experiments is meaningless without a viscous correction. This correction is performed here by using the technique developed by Wai and Yoshihara (1980); see Tseng (1983 b) for details.] The two-dimensional configuration under consideration is that of a RAE 2822 airfoil at an angle of attack c~= 3.19 ° and Mo~ = 0.730. The discretization corresponds to N1 = 30. Figure 6 presents the pressure distribution on the upper and lower surface of the airfoil as obtained by Tseng (1983 b) (with and without viscous correction) along with the experimental results of Cook et al. (1979). The agreement is very good.
5 Concluding remarks A new integral-equation formulation for the wave equation with a boundary moving in arbitrary motion has been presented. The formulation has been applied to the analysis of unsteady threedimensional compressible potential flows about airplanes and helicopters in arbitrary motion [Eq. (20) for airplanes, Eq. (27) for helicopter rotors in hover, Eq. (9) for arbitrary rigid-body motion, and Eq. (B.4) for completely arbitrary motion]. The formulation differs from preceding ones in that it includes, as a coherent package, the following items: (i) arbitrary motion of the frame of reference; (ii) formulation for wake transport; and (iii) nonlinear (transonic) terms of the differential equation. It may be worth noting that the formulation as presented above is directly applicable to aeroacoustics. In this case, once the integral equation has been solved [Eq. (20) for airplanes, Eq. (27) for helicopter rotors in hover, Eq. (9) for arbitrary rigid-body motion, and Eq. (B.4) for completely arbitrary motion], the integral representation may be used to evaluate q5 anywhere in the far field, and hence the pressure, using the Bernoulli theorem. Each of the above items has been validated independently of the other two (using a zeroth-order boundary-element algorithm); more precisely, validations of the formulation have been presented for: (i) prescribed-wake analysis of compressible flows about helicopter rotors in hover; (ii) free-wake analysis of incompressible flows about helicopter rotors in hover; and (iii) prescribed-wake analysis of transonic flows about airplanes wings. It should be noted that, in the free-wake analysis, steady-state flows have been studied using a time-accurate, time-marching technique. The validation includes a study of the numerical convergence as well as a comparison with existing numerical and experimental results. Additional work not discussed here includes recent developments of the formulation for steady and unsteady supersonic flows; numerical results are available in Tseng and Freedman (1985). Additional work includes an extension of the potentialflow formulation to viscous flows based on the use of the Helmholtz scalar/vector-potential decomposition (Morino 1986; Morino and Bharadvaj 1985). A completely independent approach (boundary-integral formulation for viscous flows in primitive variables) is presented in Piva and Morino (1987).
Appendix A For incompressible potential flows, the formulation governing the transport of the wake and the equation for the potential discontinuity A~b may be obtained by using the formulation introduced
L. Morino et al.: Boundary integral equation for wave equation with moving boundary
241
by Suciu and Morino (1977), and refined by Morino, Kaprielian and Sipcic (1985) and Morino and Bharadvaj (1985). In this Appendix, this formulation is extended to compressible potential flows; details are available in Morino et al. (1987b). From a physical point of view, the wake is the region where the vorticity generated at the surface of the body is transported. In potential aerodynamics this region has zero thickness. A detailed analysis of the relationship between the mathematical representation (potential wakes) and physical reality is given in Morino (1986). Here it is sufficient to point out that the proof that an inviscid isentropic initially-irrotational flow remains irrotational at all times fails for the fluid points that come in contact with the surface of the body. These points form a surface called the wake. This fact implies that the flow is not necessarily potential for the points of the wake. Therefore, in deriving the integral representation for q~, Eq. (14), the surface a should surround both the body and the wake. The contribution of the wake may be isolated by taking the limit as the two sides of the surface surrounding the wake approach each other. This yields an integral over the (open) surface of the wake which contains the discontinuities on (1) ~, (2) its tangential derivative, (3) its normal derivative, aO/~n, and (4) ~4/~t. For wakes, the discontinuity on ~(o/~n is identically equal to zero. The potential discontinuity A~b (and hence those on the tangential and time derivatives of A~b)may be obtained by starting from the condition that the pressure discontinuity is equal to zero and using Bernoulli's theorem to express the pressure in terms of the velocity potential ~b. This yields a boundary condition on the wake as
+ vw" V~ ) A(~= D D~ A0 :
0,
(A.1)
where Vw indicates the velocity of a point of the wake, defined as vw = (i~1 -~- v2)/2, with 1 and 2 indicating the two sides of the wake. The above equation implies that A4) remains constant following a point having velocity vw. The value of Aq5 is obtained by using the Joukowski hypothesis and is equal to the value that it had when it left the trailing edge (see Morino 1986 for a clarification of this point). The above results are used to study the motion of the wake as follows: the wake is divided into elements that are convected by the flow; the value of Aq~connected with each panel remains constant in time (and is obtained from the value of the discontinuity on q5 at the trailing edge). The location of each node on the wake is obtained by displacing it, at any time step, by vwA-c,where vwis obtained by taking the discretized form of the gradient of Eq. (20) (for airplanes) or Eq. (27) (for helicopter rotors). It should be pointed out that, in general, the wake cannot be assumed to be fixed in the body frame of reference (although this is true in the cases considered here). Therefore, in general the formulation for a surface moving with respect to the frame of reference (presented in Appendix B) may be required even if the body moves in rigid-body motion (e. g., a rigid helicopter rotor in forward flight).
Appendix B In this Appendix the extension of the integral representation given by Eq. (14) to the case in which the surface moves in arbitrary motion with respect to the body frame of reference (obtained by Morino et al. 1985) is outlined, Let S(x, t) be a continuously differentiable function that is positive outside o-, negative inside o., and zero on o.. Thus, the equation S(x, t) = 0 defines the motion of the surface with respect to the moving frame of reference. Note that, E(x, t) = H(S(x, 0), where H( ) is the Heaviside function, and VE = ~(S)VS, where g ( S ) = dH/dS. Note that, for any smooth function f, f~f6(S)dV = ~ ~
do-
--oO
Hence, for any smooth function f, using the notation [ ]~o= [ ]~=~o,one obtains
(B.1)
242
Computational Mechanics 4 (1989)
[.S~[fVEltodV= [.S~[fVSf(S)]todV= ~j" [fVS]~o6(Sto)dV --o0
--o0
--o0
~o,,, [fVS]to IVStol]dato = ~,o [.fn]to[VS]to]vst{odato = ~,o [fn]toddto, where d 6 to =
Sto = S[x, to(X,X,,t,)]
whereas the 'retarded surface,'
¢7to,
(B.2) is defined by
S,o = 0 and
(IV Sl,o/IVS,ol)d%.
Similarly,
oo fdE
oo
o~ fdS
?,]1
= f~,,, dt toIV~ol dat°=-~'f~,o [fvA°lVSLda'ool=-lI~,lVS,
,~[fv~,ltodd,o,
(B.3)
where v~ = - (dS/dt)/[VS[ is the velocity of the surface a. Combining Eqs. (9), (B.2), and (B.3) one obtains
E(x,,t,)(o(x,,t,)=~f,%[~4_~n]tddt,+V, 1 d
[-1
+a2 ~,~;;~,o ~6q~v~
-1 ] Jtl
de,,+~
~[
~
1;
-1 ] E 4 ~ z dV
o[~d~ (B.4)
ti
This is the desired generalization of Eq. (14). If the surface is rigid, then V, and d/dt, may be moved under the integral sign. Using V, [q~lt~= V, ti(O4/~t)and a similar expression for d[(o]t]dt,, Eq. (B.4) reduces to Eq. (14). It may be noted that Eq. (B.4) is considerably more difficult to use than Eq. (14). Therefore, a convenient frame of reference to be used in connection with a given surface is the one with respect to which the motion of the surface is small, so that Eq. (B.3) may be approximated with Eq. (14). It may be worth noting that, if more than one surface is involved, a different frame of reference may be used for each surface. In particular, for the case of a helicopter rotor-fuselage configuration, it is convenient to use a rotor frame of reference for the rotor, and a fuselage frame of reference for the fuselage. The best frame of reference for the wake depends on the type of flow. For instance, for an isolated rotor in hover, the wake is fixed in the rotor frame of reference and therefore such frame of reference is the most convenient one in this case. On the other hand, for an isolated rotor in forward flight, the motion of the wake in the air frame of reference may be treated as being small, in some cases negligible: in this case the air frame of reference may be the most convenient to express the wake integral.
Acknowlegements This work was partially supported by NASA Contract No. NAS 1-17 317 to ICARUS Inc. and by contracts DAAG-29-83K-0050 and DAAL-03-87-K-0022 from the U. S. Army Research Office to Boston University.
References Baley, F.R.; Steger, J.L. (1972): Relaxation techniques for three-dimensional transonic flows about wings. AIAA Paper 72-189 Caradonna, F.X. (1985): Finite-difference method for the solution of unsteady potential flows. In: Morino, L. (ed.) Computational methods in potential aerodynamics. Southampton/UK: Comput. Mech. Publ. Caradonna, F.X.; Tung, C. (1981): Experimental and analytical studies of a model helicopter rotor in hover. USAAVRADCOM TR-81-A-23
L. Morino et al.: Boundary integral equation for wave equation with moving boundary
243
Cook, P. H.; McDonald, M.A.; Firmin, M. C. P. (1979): Aerofoil RAE 2822 - pressure distribution and boundary layer and wake measureements. Barche, J. (ed.) Experimental data base for computer program assessment. AGARD-Advisory Report No. 138 Farassat, F. (1984): A new aerodynamic integral equation based on an acoustic formula in the time domain. AAIA J. 22, 1337-1340 Farassat, F.; Myers, M.K. (1987): Extension of Kirchhoff's formula to radiation from moving surfaces. NASA Technical Memorandum 89149 Ffowcs Williams, J.E.; Hawkings, D.L. (1969): Sound generated by turbulence and surfaces in arbitrary motion. Philos. Trans. R. Soc. London Ser. A 264, 321-342 Hanson, D.B. (1983): Compressible helicoidal surface theory for propeller aerodynamics and noise. AIAA J. 21,881-889 Long, L.N. (1983): The aerodynamics of propellers and rotors using an acoustic formulation in the time domain. AIAA Paper 1821, AIAA Applied Aerodynamics Conf. Morino, L. (1974): A general theory of unsteady compressible potential aerodynamics. NASA CR-2464 Morino, L. (1985): Foundation of potential flows. In: Morino, L. (ed.) Computational methods in potential aerodynamics. Southampton/UK: Comput. Mech. Publ. Morino, L. (1986): Helmholtz decomposition revisited: vorticity generation and trailing edge condition. 1. Incompressible flows. Comput. Mech. 1, 65-90 Morino, L.; Bharadvaj, B. (1985): Two methods for viscous and inviscid free-wake analysis of helicopter rotors. Boston University, Center for Computational and Applied Dynamics, CCAD-TR-85-02-R Morino, L.; Chen, L.T.; Suciu, E. O. (1975): Steady and oscillatory subsonic and supersonic aerodynamics around complex configurations. AIAA J. 13, 368-374 Morino, L.; Freedman, M.I.; Deutsch, D.J.; Sipcic, S.R. (1985): An integral equation method for compressible potential flows in an arbitrary frame of reference. In: Morino, L. (ed.) Computational methods in potential aerodynamics. Southampton/UK: Comput. Mech. Publ. Morino, L.; Freedman, M.I.; Tseng, K. (1987a): Unsteady three-dimensional compressible potential aerodynamics of helicopter rotors - a boundary-element formulation. American Helicopter Society Specialists Meeting on Aerodynamics and Aerocaustics, Arlington TX, February 1987 Morino, L; Freedman, M. I.; Tseng, K. (1987 b): Unsteady compressible potential aerodynamics of rotors. Boston University, Center for Computational and Applied Dynamics, CCAD-TR-87-02. Morino, L.; Kaprielian, Jr., Z.; Sipcic, S.R. (1985): Free-wake analysis of helicopter rotors. Vertica 9, 127-140. Morino, L.; Kuo, C.C. (1974): Subsonic potential aerodynamics for complex configurations: a general theory. AAIA J. 12, 191-197 Morino, L.; Tseng, K. (1982): Nonlinear Green's function method for unsteady transonic flows. In: Nixon, D. (ed.): Transonic aerodynamics. American Institute of Aeronautics and Astronautics, New York/NY, USA, 565-603 Piva, R.; Morino, L. (1987): Vector Green's function method for unsteady Navier-Stokes equations. Meccanica, 22, 76-85 Steinhoff, J.; Ramachandran, K. (1986): Free-wake analysis of compressible rotor flow fields in hover. Twelfth European Rotorcraft Forum, Garmisch-Partenkirchen, FRG Suciu, E.O.; Morino, L. (1977): Nonlinear steady incompressible lifiting-surface analysis with wake roll-up. AIAA J. 15, 54-58 Tai, H.; Runyan, H.L. (1984): Lifting surface theory for a helicopter rotor in forward flight. NASA TM 86 315 Tseng, K. (1983a): Nonlinear Green's function method for transonic potential flow. Ph.D. Dissertation, Mathematics Department, Boston University Tseng, K. (1983 b): A Green's function method for viscous transonic flows. ICARUS Inc., Rept. 83-02 Tseng, K. (1984): Application of the Green's function method for two- and three-dimensional steady transonic flows. AIAA Paper 84-0425 Tseng, K.; Freedman, M.I. (1985): A first order Green's function method for supersonic oscillatory flows. ICARUS Inc., Rept. 85-01 Wai, J.C.: Yoshihara, H., (1980): Viscous transonic flows over airfoils. In: Reynolds, W.C.; Maccormack, R.W. (eds.) Seventh int. conf. numerical methods in fluid dynamics, Stanford/CA, 1980. Berlin, Heidelberg, New York: Spriltger (Lecture Notes in Physics, vol. 141) Williams, M.H.; Hwang, G. (1986): Three dimensional unsteady aerodynamics and aeroelastic response of advanced turboprop. AIAA Paper 86-0846, AIAA/ASME/ASCE/AHS 27th Structures, Structural Dynamics and Materials Conference, pp. 1%21