- Email: [email protected]
Solution of the transonic integrodifferential equation using a decay function
Solution of the transonic integrodif%erential equation using a decay function W. Ogana Department
of Mathematics,
University
of Nairobi,
Nairobi,
Kenya
The two-dimensional, small-disturbance, fransonic integrodifferentiaf equation is solved by using a decay function that is applicable in the entire computational domain rather than on the airfoil surface only. The computational domain is discretized into rectangular elements, and the integrals involving the dependent variable and its derivatives are approximated by boundary element methods in each rectangle. Besides constant elements, use is made of hybrid elements that are based on constant elements in the streamwise direction and variable elements in the transverse direction. The methods are tested for parabolic-arc and NACAO012 airfoils in nonlifting subcritical and supercritical flows. The results compare favorably with finite difference solutions. It is found that even with only one vertical level of rectangular elements, and hence a relatively small number of nodes, the accuracy is still good. Keywords:
transonic flow, boundary elements, decay function
Introduction
Basic equations
Spreiter and Alksne’ developed a method of solving the transonic integral equation that has remained in use for a long time. According to this method, a decay function is introduced so that solution of the integral equation reduces to determination of values on the wing surface only. The method and its variations have been applied to a number of flow phenomena (see, for instance, Refs. 2-4). To extend the domain of application and improve accuracy, Ogana’ derived a twodimensional decay function that can be used in front and aft of the airfoil. Extension to three-dimensional flow soon followed.6 Although direct methods have been developed that determine solutions over the entire computational domain (see, for instance, Refs. 7 and 8), the decay function technique remains attractive for small disturbance problems because it results in a considerable reduction in the number of nodes, storage space, and computation time but still yields fairly good accuracy. In this paper we use recent developments in boundary element methods to solve the two-dimensional transonic integrodifferential equation for nonlifting flow, using the decay function. In addition to the conventional constant elements we derive hybrid elements based on constant elements in the streamwise direction and variable elements in the transverse direction, leading to what we have designated as constant-linear and constant-quadratic elements.
The two-dimensional transonic small-disturbance equation, with inclusion of artificial viscosity, may be written in the form9~i0
Address ematics, Received
30
reprint requests to Dr. Ogana at the Department University of Nairobi, Nairobi, Kenya. 31 January
1989; accepted
Appl. Math. Modelling,
8 September
of Math-
1989
1990, Vol. 14, January
where 4(x, y) is a modified perturbation tential. The term H(x, y) is defined by
velocity po-
H(x, y) = ll2 + 2/_&P
(2)
in which P =
Ax;@ 0
CL=
1 1
- ;u2) if if
ul
(4)
and Ax is the streamwise grid spacing. The velocity components are obtained from U(X,Y) =
$xkY)
06,
Y) = 4J&G
Y)
(5)
Equation (1) is solved subject to the body boundary condition 05x51
otherwise
(6)
where Y,(x) defines the modified airfoil profile. Application of Green’s theorem to equation (1) followed by integration by parts (see, for instance, Refs. 10 and 11) yields an integrodifferential equation that
0 1990 Butterworth
Publishers
Transonic integrodifferential
equation:
W. Ogana
may be written
(7)
where @ is a line integral over the airfoil section and consists of lifting and thickness distribution effects. The kernel is given by
m -
x7
5 - Y) =
-65-x) 45-R - x)2 + (i - Y121
(8)
In the decay function technique we seek solutions at y = 0, so we solve equation (7) in the form 4(x, 0) = @(x, 0) + j-1 K(5 - x, 5)H(& 5) & (9) Since the term H(x, y) is a function of u(x, y) and its derivatives, we must relate u(x, y) to u(x, 0) by what is known as a decay function. Following Ogana,s we let
4&Y) =
a,
- Y: (x)/u
b(x) = 2(X) -
05x5 1 otherwise
TI;[(u - l)U,
in which u and its derivatives
Discretization
of the computational
domain
(10)
1 + a(x)y + b(x)y2
0
-+X
Alrtoll
Figure 1
0)
where u(x) and b(x) are real nonnegative functions. A property of the decay function is that it reduces u(x, y) to u(x, 0) as y += 0 and lets u(x, y) vanish as y + SC. By requiring that equation (10) should satisfy both the irrotationality condition and the transonic partial differential equation in the neighborhood of the x-axis, we find that a(x) =
1,
(11)
Let the node Qj have coordinates (xi, yJ. Application of equation (9) at the node Qj on the x-axis yields $; = @ + $J
Zj
i=
1,2,...,m
(13)
j=l
where
and (14)
+ u3
(12)
are evaluated at y = 0.
Discretization and approximation To solve equation (9) numerically, we divide the computational domain into m horizontal strips parallel to the x-axis with each strip divided into n rectangular elements (Figure 1). Considerable reduction in the total number of rectangular elements, hence nodes, may be achieved by a suitable choice of fewer, and thus larger, elements away from the airfoil. Let M be the total number of rectangular elements denoted A,, A2, . . . , AM starting from the bottom left and proceeding from left to right in each horizontal strip. Then A4 = mn. Let Aj have width ZSj, height 2hj, aspect ratio Lyj = hj6j, and center (Xj, yj). Solutions for the velocity potential are obtained at the nodes Q,, Q2, . . . , Qm located on the x-axis. To implement the numerical scheme, we introduce, off the x-axis, the nodes labelled Qm+ , , Q m+2, f . . 9 QN, from left to right in each row. These nodes are required only for computational convenience because equation (9) does not permit solving for the velocity potential away from the x-axis.
Boundary elements The different boundary element formulas for solving equation (13) depend on how the integral in equation (14) is approximated. Except near shocks and airfoil edges, the velocity component u does not display much variation over short distances in the streamwise direction. We may thus assume that within a given vertical strip, the velocity u and its derivatives are constant in the streamwise direction and take values at the x-coordinate corresponding to the centers of the rectangular elements located within that strip. A local coordinate, C,with origin at (Xj, I;), as shown in Figure 2, can be defined through the transformation I = (i - yj)lhj
(15)
With the application of the local coordinate the righthand side of equation (14) can be integrated with respect to 5 to yield 1
4
zz
I
W,(t)H(X-,
t) dt
(16)
-I
APP~. Math. Modelling,
1990, Vol. 14, January
31
Transonic integrodifferential
equation:
W. Ogana
Substitution
T-7
Ij = h;Hjl + h$Hjz
------
%
(Xj
h$. =
&(t)W,W dt
k= 1,2
(23)
-I
-----_ YI,
is an influence coefficient that defines the interaction between the node Qi and the particular node labelled k in Aj. The influence coefficients may be determined by direct integration. I2 Using equations (13) and (22), we obtain a system identical to that in equation (19) except that N = m(n + I) and the element b, has to be deduced from the influence coefficients.
,I ” Figure 2
(22)
where
t
3
of equation (20) into equation (16) yields
Location of nodes
in which Wti(t) is a logarithmic function that depends on the coordinates (Xj, I$), (Xi, 0) and the geometry Of
Constant-quadratic elements In equation (16) we assume a quadratic variation of H(Xj, t) for t E [ - 1, I]; theni H(Xj, t) = 01Hi, + &Hj, + &Hi,
(24)
Aja
Constant elements
If we assume that H(X,, t) is constant for t E [ - 1, l] and takes the value at the node located in Aj, then we obtain the constant element formulation. If Aj touches the x-axis, the node is situated at the point labelled 1 in Figure 2; but if Aj is away from the x-axis, the node is situated at the center of Aj, namely, at the point labelled 3. Equation (16) reduces to Ij = boHj
where Hjj, Hj2, and Hj3 are the values of H at the nodes labelled 1, 2, and 3, respectively, in Figure 2 and 8,, e2, and 0, are the shape functions
e,(t) = tt(t - 1) e2(t) = ft(t + 1) e3(t) = (1 - t)(l + t)
(25)
(17)
where Hj E H(Qj) and bij =
Wti(t) dt
(18)
-1
Equation (18) can be integrated directly.12 Substitution of equation (17) into equation (13), followed by appropriate evaluation of Hi, leads to the following system of nonlinear equations in the unknowns &, &, . . . , &: +i = fl
+ 5
bcHj
i= 1,2,...,m
(19)
j=l
in which [b,] is an m x N matrix, where N = mn. Since the matrix elements b, do not involve the dependent variable, they can be computed once and stored. Constant-linear elements In equation (16) we assume a linear variation H(Xj, t) for t E [ - 1, 11; thenI H(Xj, t) = 81Hjl + 02Hj2
of (20)
where Hjl and Hj2 are the values of H at the nodes labelled 1 and 2, respectively, in Figure 2 and 8, and 0, are the shape functions e,(t) = $(l - t)
32
Appl.
e*(t) = $(l + t)
Math. Modelling,
(21)
1990, Vol. 14, January
Finik Difforenoe I---D - LY - Q- -Constant-Linoar
Elrmentr
Figure 3 Coefficient of pressure for a parabolic thickness ratio 0.06 at Mm = 0.825
arc airfoil of
Transonic integrodifferential
Substitution
of equation (24) into equation (16) yields
(26)
!j = h$Hjl + h$Hjz + h$Hj, where the influence coefficient
proximated. It is sufficient to solve the system by Jacobi iteration, namely, N
h$ is defined as
($p+
1) =
@ + 2 boHjp’
1
“I=
I &M~fj(t)
h$ =
dt
k = 1,2,3
(27)
-I
and can be evaluated by direct integration.‘* Using equations (13) and (26), we obtain a system identical to that in equation (19) except that N = m(2n + 1) and the element b, has to be deduced from the influence coefficients.
Iteration scheme and shocks It has been established that all formulations lead to the nonlinear system in equation (19), in which the size of the matrix [b,] depends on how equation (14) is ap-
W. Ogana
equation:
p = 0, 1,2, . . .
(28)
I
where the superscript p indicates the iteration level. We now turn our attention to the evaluation of H(@) during the pth iteration. From equations (2) and (3) we observe that H is a function of u and its streamwise derivatives. Hence we should be able to evaluate H(xj, yj) once we know velocities along a row that contains the node Qj. This leads to consideration of how to determine I((x~,Yj). If yi = 0, then at each iteration, U(Xj,0) is determined by differentiating the velocity potential according to equation (5); and if Yj # 0, then u(xj, yj) is determined from equation (10). The streamwise derivative of a function F(x, y) at the node Qj may be approximated by the second-order central difference formula of the type
%
I
I
I
I
I
% -0.5
-
-0-S
I I (
-o-4-
C
-0*3-
-0.4
-0.3
-0.2-0.2
0
-o-i0.1 k
I O0
G I
I
O-l-
I
0.1,
o-2,
G
Finite
Dlfhwonce 4
c 0.3,
-A r
0 Figure 4
-A
Finite
- A- -Constant-Linear d-2
d4
-A
Difference Elements
d6
Coefficient of pressure for NACA0012
at M,
Elements
P
d o’.s
- A- - A -Constant-Linear
1.: = 0.63
A.2
d.4
0:s
o!*
Figure 5 Coefficient of pressure for a parabolic thickness ratio 0.06 at M, = 0.87
Appl.
Math.
Modelling,
1990,
1
arc airfoil of
Vol. 14, January
33
Transonic integrodifferential
equation:
W. Ogana
(aF/a~)~ = (F,+& - F”JAx
or the first-order
backward
(29)
difference
formula
(aF/&+ = (Fj - Fj_ ,)lAx At each iteration, evaluation of
Parabolic arc, Mm = 0.87
(30)
solution of equation (28) requires
ffj = <+:,‘>j+ 2/LjPj
(31)
at node Qj. Before this evaluation is done, a test is carried out to determine whether Qj is in a subsonic or supersonic flow region. If Qj is in a subsonic flow region, the term 6 in #I; is evaluated by equation (29), where we take +jk* = (+j + +je1)/2
(32)
and if Qj is in a supersonic flow region, we use equation (30). Everywhere the derivatives in Pj are approximated by equation (29) to yield Pj = Ax[~ - (+j+, - (bj-,)/2Ax]
x [(+j+l - 2$j + (bj-,)l(Ax)‘l Evaluation
Table 1. Magnitudes of the differences in coefficient of pressure between constant-linear and other elements
of 4: by central differences
(33)
throughout
the
NACA0012, Mm = 0.63
X
&
d cq
dc
d cq
0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
46 99 36 59 65 64 79 75 69 61 52 41 34 30 12 67 55 3 56 46
3 5 2 4 3 3 4 2 2 0 2 5 8 11 5 9 7 6 4 3
12 11 8 7 8 10 12 15 16 18 19 19 20 20 19 18 17 15 12 10
2 2 1 1 2 2 2 2 2 2 2 3 2 2 2 3 2 1 2 1
d, x 10m4 = difference between constant and constant-linear elements results; dcq x 10e4 = difference between constantquadratic and constant-linear elements results.
Table 2. Relative errors in the coefficient of pressure due to the decay function NACAOOI 2, M, = 0.63
Parabolic arc, M, = 0.87 X
-O-6-
. - A- -A
I
34
Appl.
Finite
I
0.4
Elements
1 0.6
Math.
Modelling,
1990,
rcs
r,
rd
rcq
0.025 0.075 0.125 0.175
0 26 0 9
0 0 0 9
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0.275 0.225 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
4 6 3 6 5 5 6 8 11 15 11 85 72 111 0 4
4 0 3 0 2 2 2 4 4 8 4 34 43 56 24 0
0 0 3 0 2 0 2 2 2 0 11 11 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 8
8 0 0 0 0 0 0 0 0 0 0 0 0 0 0
8 0
at M,
domain would yield a weaker shock located upstream of the correct position.
computational
Results and discussion
I 0.8
Coefficient of pressure for NACA0012
rd
r, x 10m4 = relative error due to constant elements; r,, x 10m4 = relative error due to constant-linear elements; rcq x 1O-4 = relative error due to constant-quadratic elements.
Diftrrenccr
- -LS-Constant-Linror
0.2 Figure 6
%
r,
= 0.8
Vol. 14, January
The methods described here were tested on parabolic arc and NACAO012 airfoils in both subcritical and su-
Transonic integrodifferential
percritical flows. Reduction in computing time was achieved by using only one horizontal strip containing 50 rectangular elements, of which 20 were located on the airfoil. The number of nodes was therefore 50,100, and 150 for the constant, constant-linear, and constantquadratic elements, respectively. The use of one horizontal strip was based on an earlier finding by Oganai4 that it yields good accuracy, particularly when applied in conjunction with hybrid elements. Iteration stopped when the relative error between successive values of u was less than 0.01 at all nodes. Figures 3 and 4 show the coefficient of pressure plots for the parabolic arc and NACAO012 airfoils, respectively, in subcritical flow. Figures 5 and 6 show results for supercritical flows. The finite difference results in Figures 3 and 5 were due to Ballhaus et a1.15 while those in Figure 6 were due to Jameson.lh For the current method we graph only the results from constant-linear elements because, as Table 1 shows, the differences between them and the other results are too small to make much impact graphically. However, Table I shows that constant-linear and constant-quadratic element results are closer to one another than they are to constant element results. To determine how significantly the results obtained by use of the decay function deviate from those obtained by direct solution of equation (I), we generated relative errors from the two sets of results as shown in Table 2 for two flow phenomena. Direct solution of equation (1) using rectangular elements identical to those for the current method was carried out by 0gana,14 and the coefficient of pressure values was given correct to four decimal places for the purpose of generating
W. Ogana
ommended in order to reduce the number of nodes and hence computation time.
References 1
2 3 4
5 6
7 8 9
10 11
Table 2.
12
Conclusion
13
Application of the decay function together with hybrid elements yields fairly good results even if the computational domain is discretized into one horizontal row of rectangular elements. Owing to its relatively high accuracy, the constant-quadratic element is recommended for use with one horizontal strip of rectangular elements. However, if the use of more horizontal strips is required, the constant-linear element is rec-
equation:
14 15 16
Spreiter, J. R. and Alksne, A. Y. Theoretical prediction of pressure distribution on nonlifting airfoils at high subsonic speeds. NACA Rept. No. 1217, 1955 Heaslet, M. A. and Spreiter, J. R. Three-dimensional transonic flow theory applied to slender wings and bodies. NACA Rept. No. TN3717, 19.56 Nixon, D. and Hancock, G. J. High subsonic flow past a steady two-dimensional aerofoil. ARC CP Rept. No. 1280, 1974 Takanashi, S. A solution of the three-dimensional transonic integral equation. Proceedings offhe Eleventh Fluid Dynamics Conference, Japan Society of Aero Space Science, 1979, pp. 86-89 Ogana, W. Choosing the decay function in the transonic integral equation. ICTP Internal Rept. No. IC/83/21, Trieste, 1983 Takanashi, S. A numerical solution of the transonic integral equation and its application to three-dimensional transonic wing design. Proceedings of the Second NAL Symposium on Aircraft Computational Aerodynamics, NAL SP-3, Japan, 1984, pp. 280-291 Ogana, W. Solution of transonic flows by an integro-differential equation method. NASA Rept. No. TM 78490, 1978 Nixon, D. Direct numerical solution of the transonic perturbation integral equation for lifting and non-lifting flows. NASA Rept. No. TM 78518, 1978 Jameson, A. Computational fluid dynamics. VKI Lecture Series, von Karman Institute for Fluid Dynamics, Rhode-st-Genese, Belgium, 1976 Ogana, W. Transonic integro-differential and integral equation with artificial viscosity. ICTP Internal Rept. No. IC/87/87, Trieste, 1987 Niyogi, P. Integral Equafion Method in Transonic Flow. Lecture Notes in Physics, Vol. 157. Springer-Verlag. New York, 1982 Ogana, W. Boundary element solution of the transonic integrodifferential equation involving a decay function. ICTP Internal Rept. No. IC/87/84, Trieste,-1987 _ Brebbia. C. A.. Telles. J. C. F. and Wrobel. L. C. Boundarv Element Techniques. Springer-Verlag, Berlin, 1984 ’ Ogana, W. Boundary element solution of the transonic integrodifferential equation. ICTP Internal Rept. No. IC/87/86, Trieste, 1987 Ballhaus, W. F., Jameson, A. and Albert, J. Implicit approximate-factorization schemes for the efficient solution of steady transonic flow problems. NASA Rept. No. TM X-73202, 1977 Jameson, A. Transonic airfoil calculations using the Euler equations. Numerical Methods in Aeronautical Fluid Dynamics. ed. P. L. Roe. Academic Press, London, 1982
Appl. Math. Modelling,
1990, Vol. 14, January
35
of Mathematics,
University
of Nairobi,
Nairobi,
Kenya
The two-dimensional, small-disturbance, fransonic integrodifferentiaf equation is solved by using a decay function that is applicable in the entire computational domain rather than on the airfoil surface only. The computational domain is discretized into rectangular elements, and the integrals involving the dependent variable and its derivatives are approximated by boundary element methods in each rectangle. Besides constant elements, use is made of hybrid elements that are based on constant elements in the streamwise direction and variable elements in the transverse direction. The methods are tested for parabolic-arc and NACAO012 airfoils in nonlifting subcritical and supercritical flows. The results compare favorably with finite difference solutions. It is found that even with only one vertical level of rectangular elements, and hence a relatively small number of nodes, the accuracy is still good. Keywords:
transonic flow, boundary elements, decay function
Introduction
Basic equations
Spreiter and Alksne’ developed a method of solving the transonic integral equation that has remained in use for a long time. According to this method, a decay function is introduced so that solution of the integral equation reduces to determination of values on the wing surface only. The method and its variations have been applied to a number of flow phenomena (see, for instance, Refs. 2-4). To extend the domain of application and improve accuracy, Ogana’ derived a twodimensional decay function that can be used in front and aft of the airfoil. Extension to three-dimensional flow soon followed.6 Although direct methods have been developed that determine solutions over the entire computational domain (see, for instance, Refs. 7 and 8), the decay function technique remains attractive for small disturbance problems because it results in a considerable reduction in the number of nodes, storage space, and computation time but still yields fairly good accuracy. In this paper we use recent developments in boundary element methods to solve the two-dimensional transonic integrodifferential equation for nonlifting flow, using the decay function. In addition to the conventional constant elements we derive hybrid elements based on constant elements in the streamwise direction and variable elements in the transverse direction, leading to what we have designated as constant-linear and constant-quadratic elements.
The two-dimensional transonic small-disturbance equation, with inclusion of artificial viscosity, may be written in the form9~i0
Address ematics, Received
30
reprint requests to Dr. Ogana at the Department University of Nairobi, Nairobi, Kenya. 31 January
1989; accepted
Appl. Math. Modelling,
8 September
of Math-
1989
1990, Vol. 14, January
where 4(x, y) is a modified perturbation tential. The term H(x, y) is defined by
velocity po-
H(x, y) = ll2 + 2/_&P
(2)
in which P =
Ax;@ 0
CL=
1 1
- ;u2) if if
u
(4)
and Ax is the streamwise grid spacing. The velocity components are obtained from U(X,Y) =
$xkY)
06,
Y) = 4J&G
Y)
(5)
Equation (1) is solved subject to the body boundary condition 05x51
otherwise
(6)
where Y,(x) defines the modified airfoil profile. Application of Green’s theorem to equation (1) followed by integration by parts (see, for instance, Refs. 10 and 11) yields an integrodifferential equation that
0 1990 Butterworth
Publishers
Transonic integrodifferential
equation:
W. Ogana
may be written
(7)
where @ is a line integral over the airfoil section and consists of lifting and thickness distribution effects. The kernel is given by
m -
x7
5 - Y) =
-65-x) 45-R - x)2 + (i - Y121
(8)
In the decay function technique we seek solutions at y = 0, so we solve equation (7) in the form 4(x, 0) = @(x, 0) + j-1 K(5 - x, 5)H(& 5) & (9) Since the term H(x, y) is a function of u(x, y) and its derivatives, we must relate u(x, y) to u(x, 0) by what is known as a decay function. Following Ogana,s we let
4&Y) =
a,
- Y: (x)/u
b(x) = 2(X) -
05x5 1 otherwise
TI;[(u - l)U,
in which u and its derivatives
Discretization
of the computational
domain
(10)
1 + a(x)y + b(x)y2
0
-+X
Alrtoll
Figure 1
0)
where u(x) and b(x) are real nonnegative functions. A property of the decay function is that it reduces u(x, y) to u(x, 0) as y += 0 and lets u(x, y) vanish as y + SC. By requiring that equation (10) should satisfy both the irrotationality condition and the transonic partial differential equation in the neighborhood of the x-axis, we find that a(x) =
1,
(11)
Let the node Qj have coordinates (xi, yJ. Application of equation (9) at the node Qj on the x-axis yields $; = @ + $J
Zj
i=
1,2,...,m
(13)
j=l
where
and (14)
+ u3
(12)
are evaluated at y = 0.
Discretization and approximation To solve equation (9) numerically, we divide the computational domain into m horizontal strips parallel to the x-axis with each strip divided into n rectangular elements (Figure 1). Considerable reduction in the total number of rectangular elements, hence nodes, may be achieved by a suitable choice of fewer, and thus larger, elements away from the airfoil. Let M be the total number of rectangular elements denoted A,, A2, . . . , AM starting from the bottom left and proceeding from left to right in each horizontal strip. Then A4 = mn. Let Aj have width ZSj, height 2hj, aspect ratio Lyj = hj6j, and center (Xj, yj). Solutions for the velocity potential are obtained at the nodes Q,, Q2, . . . , Qm located on the x-axis. To implement the numerical scheme, we introduce, off the x-axis, the nodes labelled Qm+ , , Q m+2, f . . 9 QN, from left to right in each row. These nodes are required only for computational convenience because equation (9) does not permit solving for the velocity potential away from the x-axis.
Boundary elements The different boundary element formulas for solving equation (13) depend on how the integral in equation (14) is approximated. Except near shocks and airfoil edges, the velocity component u does not display much variation over short distances in the streamwise direction. We may thus assume that within a given vertical strip, the velocity u and its derivatives are constant in the streamwise direction and take values at the x-coordinate corresponding to the centers of the rectangular elements located within that strip. A local coordinate, C,with origin at (Xj, I;), as shown in Figure 2, can be defined through the transformation I = (i - yj)lhj
(15)
With the application of the local coordinate the righthand side of equation (14) can be integrated with respect to 5 to yield 1
4
zz
I
W,(t)H(X-,
t) dt
(16)
-I
APP~. Math. Modelling,
1990, Vol. 14, January
31
Transonic integrodifferential
equation:
W. Ogana
Substitution
T-7
Ij = h;Hjl + h$Hjz
------
%
(Xj
h$. =
&(t)W,W dt
k= 1,2
(23)
-I
-----_ YI,
is an influence coefficient that defines the interaction between the node Qi and the particular node labelled k in Aj. The influence coefficients may be determined by direct integration. I2 Using equations (13) and (22), we obtain a system identical to that in equation (19) except that N = m(n + I) and the element b, has to be deduced from the influence coefficients.
,I ” Figure 2
(22)
where
t
3
of equation (20) into equation (16) yields
Location of nodes
in which Wti(t) is a logarithmic function that depends on the coordinates (Xj, I$), (Xi, 0) and the geometry Of
Constant-quadratic elements In equation (16) we assume a quadratic variation of H(Xj, t) for t E [ - 1, I]; theni H(Xj, t) = 01Hi, + &Hj, + &Hi,
(24)
Aja
Constant elements
If we assume that H(X,, t) is constant for t E [ - 1, l] and takes the value at the node located in Aj, then we obtain the constant element formulation. If Aj touches the x-axis, the node is situated at the point labelled 1 in Figure 2; but if Aj is away from the x-axis, the node is situated at the center of Aj, namely, at the point labelled 3. Equation (16) reduces to Ij = boHj
where Hjj, Hj2, and Hj3 are the values of H at the nodes labelled 1, 2, and 3, respectively, in Figure 2 and 8,, e2, and 0, are the shape functions
e,(t) = tt(t - 1) e2(t) = ft(t + 1) e3(t) = (1 - t)(l + t)
(25)
(17)
where Hj E H(Qj) and bij =
Wti(t) dt
(18)
-1
Equation (18) can be integrated directly.12 Substitution of equation (17) into equation (13), followed by appropriate evaluation of Hi, leads to the following system of nonlinear equations in the unknowns &, &, . . . , &: +i = fl
+ 5
bcHj
i= 1,2,...,m
(19)
j=l
in which [b,] is an m x N matrix, where N = mn. Since the matrix elements b, do not involve the dependent variable, they can be computed once and stored. Constant-linear elements In equation (16) we assume a linear variation H(Xj, t) for t E [ - 1, 11; thenI H(Xj, t) = 81Hjl + 02Hj2
of (20)
where Hjl and Hj2 are the values of H at the nodes labelled 1 and 2, respectively, in Figure 2 and 8, and 0, are the shape functions e,(t) = $(l - t)
32
Appl.
e*(t) = $(l + t)
Math. Modelling,
(21)
1990, Vol. 14, January
Finik Difforenoe I---D - LY - Q- -Constant-Linoar
Elrmentr
Figure 3 Coefficient of pressure for a parabolic thickness ratio 0.06 at Mm = 0.825
arc airfoil of
Transonic integrodifferential
Substitution
of equation (24) into equation (16) yields
(26)
!j = h$Hjl + h$Hjz + h$Hj, where the influence coefficient
proximated. It is sufficient to solve the system by Jacobi iteration, namely, N
h$ is defined as
($p+
1) =
@ + 2 boHjp’
1
“I=
I &M~fj(t)
h$ =
dt
k = 1,2,3
(27)
-I
and can be evaluated by direct integration.‘* Using equations (13) and (26), we obtain a system identical to that in equation (19) except that N = m(2n + 1) and the element b, has to be deduced from the influence coefficients.
Iteration scheme and shocks It has been established that all formulations lead to the nonlinear system in equation (19), in which the size of the matrix [b,] depends on how equation (14) is ap-
W. Ogana
equation:
p = 0, 1,2, . . .
(28)
I
where the superscript p indicates the iteration level. We now turn our attention to the evaluation of H(@) during the pth iteration. From equations (2) and (3) we observe that H is a function of u and its streamwise derivatives. Hence we should be able to evaluate H(xj, yj) once we know velocities along a row that contains the node Qj. This leads to consideration of how to determine I((x~,Yj). If yi = 0, then at each iteration, U(Xj,0) is determined by differentiating the velocity potential according to equation (5); and if Yj # 0, then u(xj, yj) is determined from equation (10). The streamwise derivative of a function F(x, y) at the node Qj may be approximated by the second-order central difference formula of the type
%
I
I
I
I
I
% -0.5
-
-0-S
I I (
-o-4-
C
-0*3-
-0.4
-0.3
-0.2-0.2
0
-o-i0.1 k
I O0
G I
I
O-l-
I
0.1,
o-2,
G
Finite
Dlfhwonce 4
c 0.3,
-A r
0 Figure 4
-A
Finite
- A- -Constant-Linear d-2
d4
-A
Difference Elements
d6
Coefficient of pressure for NACA0012
at M,
Elements
P
d o’.s
- A- - A -Constant-Linear
1.: = 0.63
A.2
d.4
0:s
o!*
Figure 5 Coefficient of pressure for a parabolic thickness ratio 0.06 at M, = 0.87
Appl.
Math.
Modelling,
1990,
1
arc airfoil of
Vol. 14, January
33
Transonic integrodifferential
equation:
W. Ogana
(aF/a~)~ = (F,+& - F”JAx
or the first-order
backward
(29)
difference
formula
(aF/&+ = (Fj - Fj_ ,)lAx At each iteration, evaluation of
Parabolic arc, Mm = 0.87
(30)
solution of equation (28) requires
ffj = <+:,‘>j+ 2/LjPj
(31)
at node Qj. Before this evaluation is done, a test is carried out to determine whether Qj is in a subsonic or supersonic flow region. If Qj is in a subsonic flow region, the term 6 in #I; is evaluated by equation (29), where we take +jk* = (+j + +je1)/2
(32)
and if Qj is in a supersonic flow region, we use equation (30). Everywhere the derivatives in Pj are approximated by equation (29) to yield Pj = Ax[~ - (+j+, - (bj-,)/2Ax]
x [(+j+l - 2$j + (bj-,)l(Ax)‘l Evaluation
Table 1. Magnitudes of the differences in coefficient of pressure between constant-linear and other elements
of 4: by central differences
(33)
throughout
the
NACA0012, Mm = 0.63
X
&
d cq
dc
d cq
0.025 0.075 0.125 0.175 0.225 0.275 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
46 99 36 59 65 64 79 75 69 61 52 41 34 30 12 67 55 3 56 46
3 5 2 4 3 3 4 2 2 0 2 5 8 11 5 9 7 6 4 3
12 11 8 7 8 10 12 15 16 18 19 19 20 20 19 18 17 15 12 10
2 2 1 1 2 2 2 2 2 2 2 3 2 2 2 3 2 1 2 1
d, x 10m4 = difference between constant and constant-linear elements results; dcq x 10e4 = difference between constantquadratic and constant-linear elements results.
Table 2. Relative errors in the coefficient of pressure due to the decay function NACAOOI 2, M, = 0.63
Parabolic arc, M, = 0.87 X
-O-6-
. - A- -A
I
34
Appl.
Finite
I
0.4
Elements
1 0.6
Math.
Modelling,
1990,
rcs
r,
rd
rcq
0.025 0.075 0.125 0.175
0 26 0 9
0 0 0 9
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0.275 0.225 0.325 0.375 0.425 0.475 0.525 0.575 0.625 0.675 0.725 0.775 0.825 0.875 0.925 0.975
4 6 3 6 5 5 6 8 11 15 11 85 72 111 0 4
4 0 3 0 2 2 2 4 4 8 4 34 43 56 24 0
0 0 3 0 2 0 2 2 2 0 11 11 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 8
8 0 0 0 0 0 0 0 0 0 0 0 0 0 0
8 0
at M,
domain would yield a weaker shock located upstream of the correct position.
computational
Results and discussion
I 0.8
Coefficient of pressure for NACA0012
rd
r, x 10m4 = relative error due to constant elements; r,, x 10m4 = relative error due to constant-linear elements; rcq x 1O-4 = relative error due to constant-quadratic elements.
Diftrrenccr
- -LS-Constant-Linror
0.2 Figure 6
%
r,
= 0.8
Vol. 14, January
The methods described here were tested on parabolic arc and NACAO012 airfoils in both subcritical and su-
Transonic integrodifferential
percritical flows. Reduction in computing time was achieved by using only one horizontal strip containing 50 rectangular elements, of which 20 were located on the airfoil. The number of nodes was therefore 50,100, and 150 for the constant, constant-linear, and constantquadratic elements, respectively. The use of one horizontal strip was based on an earlier finding by Oganai4 that it yields good accuracy, particularly when applied in conjunction with hybrid elements. Iteration stopped when the relative error between successive values of u was less than 0.01 at all nodes. Figures 3 and 4 show the coefficient of pressure plots for the parabolic arc and NACAO012 airfoils, respectively, in subcritical flow. Figures 5 and 6 show results for supercritical flows. The finite difference results in Figures 3 and 5 were due to Ballhaus et a1.15 while those in Figure 6 were due to Jameson.lh For the current method we graph only the results from constant-linear elements because, as Table 1 shows, the differences between them and the other results are too small to make much impact graphically. However, Table I shows that constant-linear and constant-quadratic element results are closer to one another than they are to constant element results. To determine how significantly the results obtained by use of the decay function deviate from those obtained by direct solution of equation (I), we generated relative errors from the two sets of results as shown in Table 2 for two flow phenomena. Direct solution of equation (1) using rectangular elements identical to those for the current method was carried out by 0gana,14 and the coefficient of pressure values was given correct to four decimal places for the purpose of generating
W. Ogana
ommended in order to reduce the number of nodes and hence computation time.
References 1
2 3 4
5 6
7 8 9
10 11
Table 2.
12
Conclusion
13
Application of the decay function together with hybrid elements yields fairly good results even if the computational domain is discretized into one horizontal row of rectangular elements. Owing to its relatively high accuracy, the constant-quadratic element is recommended for use with one horizontal strip of rectangular elements. However, if the use of more horizontal strips is required, the constant-linear element is rec-
equation:
14 15 16
Spreiter, J. R. and Alksne, A. Y. Theoretical prediction of pressure distribution on nonlifting airfoils at high subsonic speeds. NACA Rept. No. 1217, 1955 Heaslet, M. A. and Spreiter, J. R. Three-dimensional transonic flow theory applied to slender wings and bodies. NACA Rept. No. TN3717, 19.56 Nixon, D. and Hancock, G. J. High subsonic flow past a steady two-dimensional aerofoil. ARC CP Rept. No. 1280, 1974 Takanashi, S. A solution of the three-dimensional transonic integral equation. Proceedings offhe Eleventh Fluid Dynamics Conference, Japan Society of Aero Space Science, 1979, pp. 86-89 Ogana, W. Choosing the decay function in the transonic integral equation. ICTP Internal Rept. No. IC/83/21, Trieste, 1983 Takanashi, S. A numerical solution of the transonic integral equation and its application to three-dimensional transonic wing design. Proceedings of the Second NAL Symposium on Aircraft Computational Aerodynamics, NAL SP-3, Japan, 1984, pp. 280-291 Ogana, W. Solution of transonic flows by an integro-differential equation method. NASA Rept. No. TM 78490, 1978 Nixon, D. Direct numerical solution of the transonic perturbation integral equation for lifting and non-lifting flows. NASA Rept. No. TM 78518, 1978 Jameson, A. Computational fluid dynamics. VKI Lecture Series, von Karman Institute for Fluid Dynamics, Rhode-st-Genese, Belgium, 1976 Ogana, W. Transonic integro-differential and integral equation with artificial viscosity. ICTP Internal Rept. No. IC/87/87, Trieste, 1987 Niyogi, P. Integral Equafion Method in Transonic Flow. Lecture Notes in Physics, Vol. 157. Springer-Verlag. New York, 1982 Ogana, W. Boundary element solution of the transonic integrodifferential equation involving a decay function. ICTP Internal Rept. No. IC/87/84, Trieste,-1987 _ Brebbia. C. A.. Telles. J. C. F. and Wrobel. L. C. Boundarv Element Techniques. Springer-Verlag, Berlin, 1984 ’ Ogana, W. Boundary element solution of the transonic integrodifferential equation. ICTP Internal Rept. No. IC/87/86, Trieste, 1987 Ballhaus, W. F., Jameson, A. and Albert, J. Implicit approximate-factorization schemes for the efficient solution of steady transonic flow problems. NASA Rept. No. TM X-73202, 1977 Jameson, A. Transonic airfoil calculations using the Euler equations. Numerical Methods in Aeronautical Fluid Dynamics. ed. P. L. Roe. Academic Press, London, 1982
Appl. Math. Modelling,
1990, Vol. 14, January
35