Graebel W.P. Advanced Fluid Mechanics
Подождите немного. Документ загружается.
294 Multidimensional Computational Methods
varied, and the airfoil is highly cambered, two such rows may be needed, one to
take care of each problem individually.
A program using the preceding procedure with source panels only to develop the
body shape is given as Program 12.3.1. Once the panel strengths are known, there is
enough information to find the velocity and pressure at any point in the exterior flow.
PROGRAM PANELBOOK
REAL*8 LAMBDA
DIMENSION ACS(40),ASN(40),BETA(40),RHS(40),THETA(40),XCP(40),
& XLEN(40),XP(40),YCP(40),YP(40),VEL(40)
DIMENSION AA(40,40),AA1(40,40),LAMBDA(40),SUM1(40,40),SUM2(40,40)
PI=3.14159265
WRITE(*,*)"This is a panel program designed to calculate the inviscid flow about"
WRITE(*,*)" an arbitrary body. The body shape is determined by the user. "
WRITE(*,*)"The program first asks for the number of panels to be used."
WRITE(*,*)"Next it asks for the X,Y coordinates of the end points of the panels."
WRITE(*,*)" These are called the nodes of the panel. For the program to run"
WRITE(*,*)" properly the nodal points should be numbered in a clockwise manner"
WRITE(*,*)" as one goes around the body."
WRITE(*,*)"'Finally the free stream velocity is requested."
WRITE(*,*)
WRITE(*,*)"The output will consist ofthe velocities calculated at the control points."
WRITE(*,*)"The control points were chosen so as to be located at the center of the panels."
WRITE(*,*)
WRITE(*,*)"Enter the number of panels you wish to use. Maximum allowable is 40. "
READ(*,*) N
WRITE(*,*)"The coordinates of the end points of the panels must now be entered."
WRITE(*,*)"his program is two dimensional, so only x and y coordinates should be entered."
WRITE(*,*)
WRITE(*,*)"You will be told which node you are entering. Panel number"
WRITE(*,*)"one has nodes one and two, panel two has nodes two and three,"
WRITE(*,*)"and the last panel has nodes n and one."
WRITE(*,*)
WRITE(*,*)"Results are in file A:PANELBOOK.DAT"
WRITE(*,*)
!============================================================
! ENTER NODES
!============================================================
DO I=1,N
WRITE(*,*) "Enter node number ",I
READ(*,*) XP(I),YP(I)
END DO
OPEN(1, FILE='A:PANELBOOK.DAT', STATUS='UNKNOWN')
WRITE(*,*) "Your entered points are listed below:"
DO I=1,N
WRITE(*,100)I,XP(I),I,YP(I)
WRITE(1,100)I,XP(I),I,YP(I)
END DO
WRITE(*,*) " Enter the free steam velocity. "
READ(*,*) U
!============================================================
! CALCULATE CONTROL POINTS
!============================================================
L=N-1
WRITE(1,*) "LOCATION OF CONTROL POINTS:"
DO I=1,L
XCP(I)=(XP(I)+XP(I+1))/2.
YCP(I)=(YP(I)+YP(I+1))/2.
WRITE(1,100)I,XCP(I),I,YCP(I)
END DO
XCP(N)=(XP(N)+XP(1))/2.
YCP(N)=(YP(N)+YP(1))/2.
WRITE(1,110)N,XCP(N),I,YCP(N)
Program 12.3.1—Panel method for flow past a nonlifting body (program by the author)
12.3 Surface Singularities 295
!============================================================
! CALCULATE ANGLES
!============================================================
DO I=1,L
IF(XP(I+1).EQ. XP(I).AND. YP(I+1).LT.YP(I)) THEN
THETA(I)=-PI/2.
ELSE IF(XP(I+1).EQ.XP(I).AND.YP(I+1).GT.YP(I)) THEN
THETA(I)=PI/2.
ELSE
THETA(I)=ATAN((YP(I+1)-YP(I))/(XP(I+1)-XP(I)))
END IF
IF (XP(I+1).LT.XP(I)) THETA(I)=THETA(I)-PI
END DO
IF(XP(1).EQ.XP(N).AND.YP(1).LT.YP(N)) THEN
THETA(N)=-PI/2.
ELSE IF(XP(1).EQ.XP(N).AND.YP(1).GT.YP(N)) THEN
THETA(N)=PI/2.
ELSE
THETA(N)=ATAN((YP(1)-YP(N))/ (XP(1)-XP(N)))
END IF
IF( XP(1) .LT. XP(N)) THETA(N)=THETA(N)-PI
!============================================================
! CALCULATE THE LENGTH AND ANGLES FOR ALL PANELS
!============================================================
DO I=1,L
XLEN(I)=SQRT( (XP(I+1)-XP(I))**2 + ( YP(I+1)-YP(I))**2 )
ASN(I)= SIN( THETA(I) )
ACS(I)=COS( THETA(I))
END DO
XLEN(N)=SQRT( (XP(1)-XP(N))**2 + (YP(1)-YP(N))**2 )
ASN(N)=SIN( THETA(N) )
ACS(N)=COS( THETA(N) )
!=======================================================================
! CALCULATE I(I,J)=INTEGRAL ( D/DN LN(R) DR(J) ) OF THE Ith PANEL
!=======================================================================
DO I=1,N
DO J=1,N
IF(J.EQ.I) THEN
SUM1(I,J)=PI
SUM2(I,J)=0.
ELSE
A=-(XCP(I)-XP(J))*ACS(J)-( YCP(I)-YP(J))*ASN(J)
B=( XCP(I)-XP(J))**2 + ( YCP(I)-YP(J) )**2
C=SIN( THETA(I)-THETA(J) )
D=COS( THETA(I)-THETA(J) )
E=( XCP(I)-XP(J))*ASN(J) - ( YCP(I)-YP(J) )*ACS(J)
F=LOG( 1.0 +XLEN(J)*( XLEN(J)+2.*A)/ B )
G=ATAN( (XLEN(J)+A)/E )-ATAN(A/E)
SUM1(I,J)=C*F/2.-D*G
SUM2(I,J)=-D*F/2.-C*G
END IF
END DO
END DO
!============================================================
! ASSIGN THE MEMBERS OF THE MATRIX EQUATION TO SOLVE FOR
! THE INDIVIDUAL PANEL STRENGTH'S.
!============================================================
DO I=1,N
DO J=1,N
AA(I,J)=SUM1(I,J)
END DO
END DO
DO I=1,N
DO J=1,N
AA1(I,J)=AA(I,J)/AA(1,1)
END DO
END DO
DO I=1,N
BETA(I)=THETA(I)+PI/2.
RHS(I)=-U*COS( BETA(I) )/AA(1,1)
END DO
Program 12.3.1—(Continued)
296 Multidimensional Computational Methods
!============================================================
! SOLVE THE MATRIX SYSTEM FOR THE PANEL STRENGTHS
! USING A GAUSS-SIEDEL ROUTINE
!============================================================
DO I=1,N
LAMBDA(I)=0.
END DO
IMAX=20
DO K=1,IMAX
ICOUNT=1
DO I=1,N
XSTAR=LAMBDA(I)
LAMBDA(I)=RHS(I)
DO J=1,N
IF (J.GT.I) LAMBDA(I)=LAMBDA(I)-AA1(I,J)*LAMBDA(J)
END DO
!============================================================
! TEST FOR CONVERGENCE
!============================================================
IF( DABS(XSTAR-LAMBDA(I)).GT.0.00001) THEN ICOUNT=0
END DO
IF (ICOUNT.EQ.1) GOTO 10
END DO
!============================================================
! CALCULATE PANEL STRENGTHS
!============================================================
10 DO I=1,N
VI=U*ACS(I)
DO J=1,N
VI=VI+LAMBDA(J)*SUM2(I,J)
END DO
VEL(I)=VI
END DO
WRITE(1,120)
DO I=1,N
WRITE(1,130) I,LAMBDA(I)
END DO
WRITE(1,140)
DO I=1,N
WRITE(6,150) I,VEL(I)
END DO
100 FORMAT(1X,'XP(',I2,')=',F10.5,1X,'YP(',I2,')=',F10.5,/)
110 FORMAT(1X,'XCP(',I2,')=',1X,F10.5,2X,'YCP(',I2,')=',F10.5)
120 FORMAT("CALCULATED STRENGTHS OF THE PANELS:")
130 FORMAT(1X,'PANEL',1X,I2,'=',1X,F10.5)
140 FORMAT(//1X,'THE SOLUTION IS GIVEN BELOW'//)
150 FORMAT(1X,'THE VELOCITY AT CONTROL POINT',1X,I3,1X,'IS',F10.5)
CLOSE(1)
END PROGRAM
Program 12.3.1—(Continued)
The use of surface singularity distributions can be continued into three dimensions
with, however, the loss of the use of complex functions. Green’s second identity in
three dimensions is
c
h
f
n
−f
h
n
dS =
V
h
2
f −f
2
hdV (12.3.22)
Let h = , where
2
=0, and f =1/R , where R =
x −x
0
i +y −y
0
j+z −z
0
k
.
Inserting this into equation (12.3.22) gives
=
1
4
S
−
1
R
n
+
n
1
R
dS (12.3.23)
12.4 One-Step Methods 297
The 1/R term represents a surface distribution of sources of strength /n, and the
second term represents a surface doublet distribution of strength . Thus, in principal,
the two- and three-dimensional methodologies are the same.
Singularity distribution techniques have been widely used by naval architects and
aircraft designers. The approach was pioneered by Hess and Smith (1962, 1967) and
widely adapted and used by others. In those industries the primary emphasis is the
external flow past bodies. The methodology chosen can either suppress the “flow”
within the body or allow for that flow. For these applications the internal flow is ignored
as being without physical interest. It is also possible to use these methods for wave
problems where there are free surfaces with or without solid boundaries (Schwartz
(1981), Forbes (1989), Zhou and Graebel (1990)).
A combination of the surface integral method and the FEM is the boundary integral
method (BIM). It uses the idea of surface singularities to solve for the interior flow and
uses paneling as in FEM to solve for the magnitude of the singularities.
Parabolic Partial Differential Equations
12.4 One-Step Methods
One-step methods take the results from a previous computation in the space dimension(s)
and advance it in the time dimension.
12.4.1 Forward Time, Centered Space—Explicit
In one space dimension, the diffusion equation becomes
T
t
=
2
T
x
2
(12.4.1)
Approximating the time derivative by
T
t
≈
T
n+1
i
−T
n
i
t
(12.4.2)
and the space derivative by
2
T
x
2
≈
T
n
i+1
−2T
n
i
+T
n
i−1
x
2
(12.4.3)
equation (12.4.1) can be approximated by
T
n+1
i
−T
n
i
t
≈
T
n
i+1
−2T
n
i
+T
n
i−1
x
2
(12.4.4)
Here, the superscript refers to the time dimension, and the subscript refers to the space
dimension (Figure 12.4.1). Solving for T at the most recent time gives
T
n+1
i
≈T
n
i
+
t
x
2
T
n
i+1
−2T
n
i
+T
n
i−1
(12.4.5)
This is an explicit (i.e., it’s not necessary to solve a system of equations, just solve one
equation at a time) equation for T in terms of its value at three points at the previous
time. It is a stable method as long as t/x
2
< 1/2. This stability criterion is called the
298 Multidimensional Computational Methods
i
1
1
i
+
1
i
–
1
n
–2
1
n
+ 1
Figure 12.4.1 Computational molecule for the forward time centered space explicit method
Courant condition. It guarantees that the solution does not “blow up,” but it does not
tell us anything about the accuracy of the method. Since equations (12.4.2) and (12.4.3)
come about from Taylor series expansions, it can be shown that accuracy is of first
order in t and second order in x. If more terms in the Taylor expansions are taken,
providing we choose t/x
2
< 1/6, the accuracy improves to second order in t and
fourth order in x.
12.4.2 Dufort-Frankel Method—Explicit
Changing the derivatives somewhat in the previous method leads to another explicit
method. Instead of equation (12.4.2), use
T
t
≈
T
n+1
i
−T
n−1
i
2t
(12.4.6)
and instead of equation (12.4.3), use
2
T
x
2
≈
T
n
i+1
−
T
n+1
i
+T
n−1
i
+T
n
i−1
x
2
(12.4.7)
Then, with =t/x
2
, equation (12.4.5) is replaced by
T
n+1
i
1
2
+
≈
T
n
i+1
+T
n
i−1
+
1
2
−
T
n−1
i
(12.4.8)
The computational molecule is shown in Figure 12.4.2. This method is explicit and
unconditionally stable (i.e., stable for all values of ). It does, however, require another
method to start up and generate the first line, and must be small, or the method is not
even first-order accurate.
12.4.3 Crank-Nicholson Method—Implicit
If we stay with equation (12.4.2) for the time derivative but replace equation (12.4.3) by
2
T
x
2
≈
T
n
i+1
−2T
n
i
+T
n
i−1
2x
2
+
T
n+1
i+1
−2T
n+1
i
+T
n+1
i−1
2x
2
(12.4.9)
instead of equation (12.4.5), we have
−T
n+1
i+1
+2
1+
T
n+1
i
−T
n+1
i−1
≈T
n
i+1
+2
1−
T
n
i
+T
n
i−1
(12.4.10)
12.4 One-Step Methods 299
n
+
1
n
n
–
1
i
+
1
ii
–
1
Figure 12.4.2 Computational molecule for the Dufort-Frankel method
where =t/x
2
. This is an implicit method, meaning a system of equations must be
solved simultaneously, but since the system is tridiagonal, the solution method is easy.
Depending on boundary conditions, it is usually unconditionally stable (i.e., stable for
all values of ), but it can be unstable if the boundary conditions include derivatives of
the independent variable. Also, if there are two or more space dimensions, it becomes
block tridiagonal, increasing the difficulty of solution.
12.4.4 Boundary Layer Equations—Crank-Nicholson
We saw earlier that similarity solutions of the boundary layer equations could be
obtained, reducing a partial differential equation to an ordinary equation. This, however,
depended on having a minimum number of dimensions and also special outer-velocity
profiles. With the present-day easy access to personal computers, a numerical program
can easily be written to solve the two-dimensional steady equations with nonsimilar
solutions with little difficulty.
An easy procedure for the two-dimensional boundary layer equations uses the
Crank-Nicolson method. Start with the steady form of the boundary layer equations—
namely,
u
u
x
+v
u
y
=
U
t
+U
U
x
+
2
u
y
2
(12.4.11)
The second-order derivative in y and the first in x tell us that this is a parabolic partial
differential equation. It is convenient to first suppress as much as is easily possible the
growth of the boundary layer thickness in the x direction. Introduce variables suggested
by our Falkner-Skan solution according to
=x/L = y
Ux
=
Ux
f (12.4.12)
Then, since
x
=
1
L
+
2
−
1
x
+
1
U
dU
dx
and
y
=
U
x
300 Multidimensional Computational Methods
we have
u =
y
=
U
x
Ux
f
=U
f
(12.4.13)
v =−
x
=−
1
2
x
U
dU
dx
+
U
x
f +
1
L
Ux
f
+
2
Ux
−
1
x
+
1
U
dU
dx
f
Substituting these into equation (12.4.1), find, after some arranging, that
U
f
dU
dx
f
+
U
L
2
f
+
U
2
−
1
x
+
1
U
dU
dx
2
f
2
−U
2
2
f
2
1
2
1
U
dU
dx
+
1
x
f +
1
L
f
+
2
−
1
x
+
1
U
dU
dx
f
=
U
2
x
3
f
3
+U
dU
dx
For convenience in writing introduce the parameter = x/UdU/dx. Multiply the
preceding equation by x/U
2
and simplify. The result is
f
f
+
2
f
−
2
f
2
+1
2
f +
f
=
3
f
3
+ (12.4.14)
It is seen that equation (12.4.14) is a third-order partial differential equation. In
numerical solutions, it is somewhat easier to deal with first-order partial differential
equations. To get to that point here, make the further change of variables
F =
f
V= v
x
U
where F is then the dimensionless x velocity component and V is a scaled dimensionless
version of the y velocity component. Then, from the second of equation (12.4.13)
find that
V =−
+1
2
f +
f
+
−1
2
f
(12.4.15)
Differentiating this with respect to gives
V
=−
F +
F
+
−1
2
F
(12.4.16)
Our momentum equation (12.4.4) can now be written as
F
F +
F
+
F
V +
2
−1F
=
2
F
2
+ (12.4.17)
To put these equations in numerical form, set up a grid with spacing by ,
with grid points at
=n −1 n = 1 2N
=m −1 m =1 2M
12.4 One-Step Methods 301
Using centered differences and keeping in mind that we want to finish up with a set of
linear algebraic equations, the various derivatives can be replaced by
V
=
1
2
V
m+1n
−V
m+1n−1
+V
mn
−V
mn−1
F =
m+1
4
F
m+1n
+F
m+1n−1
+F
mn
+F
mn−1
F
=
m+05
2
F
m+1n
−F
m+1n−1
+F
mn
−F
mn−1
−1
2
F
=
n−05
m+05
−1
4
F
m+1n
−F
m+1n−1
+F
mn
−F
mn−1
1−F
2
=
m+1
1−F
m+1n
F
mn
(12.4.18)
F
F
=
m+05
F
mn
F
m+1n
−F
mn
V +
−1
2
F
F
=
1
4
V
mn
+
n
m
−1
2
F
mn
F
m+1n
−F
m+1n−1
+F
mn
−F
mn−1
2
F
2
=
1
2
2
F
m+1n+1
−2F
m+1n
+F
m+1n−1
+F
mn+1
−2F
mn
+F
mn−1
Putting these into equation (6.6.5) find that
V
m+1n
=V
m+1n−1
−V
mn
+V
mn−1
+2P
m
F
m+1n
+Q
m
F
m+1n−1
+R
m
F
mn
+S
m
F
mn−1
(12.4.19)
with
P
m
=−
m+1
4
−
m+05
2
−
n−05
4
m+05
−1
Q
m
=−
m+1
4
−
m+05
2
+
n−05
4
m+05
−1
R
m
=−
m+1
4
+
m+05
2
−
n−05
4
m+05
−1
S
m
=−
m+1
4
+
m+05
2
+
n−05
4
m+05
−1
(12.4.20)
For the momentum equation (12.4.14), similar substitutions give
−
m+1
1−F
m+1n
F
mn
+
m+05
F
mn
F
m+1n
−F
mn
+
1
4
V
mn
+
n
m
−1
2
F
mn
F
m+1n
−F
m+1n−1
+F
mn
−F
mn−1
−
1
2
2
F
m+1n+1
−2F
m+1n
+F
m+1n−1
+F
mn+1
−2F
mn
+F
mn−1
(12.4.21)
=0
302 Multidimensional Computational Methods
These two equations can be rearranged in a form suitable for solution as
follows:
A
mn
F
m+1n−1
+B
mn
F
m+1n
+C
mn
F
m+1n+1
=D
mn
(12.4.22)
V
m+1n+1
=V
m+1n
+E
mn
(12.4.23)
with 1 <n<N, and the parameters defined by
A
mn
=−
V
mn
+05
n
m
−1F
mn
4
−
1
2
2
B
mn
=
m+1
F
mn
+
m+05
F
mn
+
1
2
C
mn
=
V
mn
+05
n
m
−1F
mn
4
−
1
2
2
(12.4.24)
D
mn
=
m+1
+
m+05
F
2
mn
−
V
mn
+05
n
m
−1F
mn
F
mn+1
−F
mn−1
4
+
F
mn+1
−2F
mn
+F
mn−1
2
2
E
mn
=−V
mn
+V
mn−1
+2P
m
F
m+1n
+Q
m
F
m+1n−1
+R
m
F
mn
+S
m
F
mn−1
This numerical integration procedure requires knowledge of F and V at the initial
line m =1. These starting conditions usually would come from a Falkner-Skan or other
of the similarity solutions previously discussed. We also must know as input data at
the various values of .
For a given value of m, the A
mn
B
mn
, and C
mn
are all known. Equation (12.4.22)
is then a tridiagonal set of algebraic equations, which is an algebraic system that is
particularly amenable to solution by a number of methods. From the boundary con-
ditions we know that F
m1
= 0V
m1
= 0, and F
mN
= 1, so equation (12.4.13) rep-
resents N −2 equations in terms of the N −2 unknowns F
mn
. Solve for the F
mn
first, and then go to equation (12.4.24) and solve for the V
mn
step by step. Then
increment m by one, and repeat the procedure. As the solution proceeds downstream,
it is a good idea to check the wall velocity gradient as the integration goes for-
ward in . When it becomes zero, flow reversal will occur downstream from this
point, and the boundary layer approximation is no longer valid. This would normally
determine M.
Numerical solution of these boundary layer equations would usually proceed
by taking a relatively coarse grid spacing (small M and N ), solving the alge-
braic system, and then repeating the solution after halving the grid spacing (dou-
bling M and N ). This procedure would repeat until the changes are acceptably
small.
The numerical procedure just outlined is not the only numerical scheme that could
have been used. It has, however, the virtue of leading to a tridiagonal matrix, which
is about the simplest algebraic system to solve. It also is a method well within the
capabilities of a personal computer and provides results that appear to be at least as
good as those found by other numerical methods.
12.4 One-Step Methods 303
12.4.5 Boundary Layer Equation—Hybrid Method
The previous methods can be combined in various ways into a hybrid method for the
boundary layer equations. Write
u
n+1j
−u
nj
+u
n+1j−1
−u
nj−1
2x
+
v
n+1j−v
−v
n+1j−1
y
=0 (12.4.25)
ru
n+1j
+1 −ru
nj
u
n+1j
−u
nj
x
+
rv
n+1j
u
n+1j+1
−u
n+1j−1
2y
+1 −rv
nj
u
nj+1
−u
nj−1
2y
=
rU
n+1
+1 −rU
n
U
n+1
−U
n
x
+
y
2
r
u
n+1j+1
−u
n+1j
−
u
n+1j
−u
n+1j−1
(12.4.26)
+ 1 −r
u
nj+1
−u
nj
−
u
n+1j
−u
n+1j−1
Notice the following:
1. r = 0: The method is explicit, and the error is first order in x and second in y.
For stability,
2x
u
nj
y
2
≤1
v
2
nj
x
u
nj
≤2
2. r = 1/2: This becomes the Crank-Nicholson method. The error is second order
in both x and y. There is no stability constraint.
3. r =1: The method now is fully implicit. The error now is first order in x and
second order in y. There is no stability constraint.
12.4.6 Richardson Extrapolation
From the preceding it is seen that there is always a trade-off in stability, accuracy, and
ease of solution. For the case of nonlinear differential equations, such as the boundary
layer equations, the nonlinearity adds additional complexity. Sometimes, however, it is
possible to make some gains in accuracy without necessarily sacrificing computational
time and complexity.
Suppose that we have a nonlinear first order system of equations
dy
dx
=fx y ya = y
0
(12.4.27)
Let u
j
be a finite difference solution at point x
j
; for example, let u
j+1
=u
j
+hfx
j
u
j
.
Further, let e
j
be the error at point x
j
e
j
=u
j
−yx
j
.Ife
0
=h
0
, where
0
does
not depend on h, then e
j
= h
x
j
+Oh
2
where
x
j
solves
d
dx
=
f
y
x yx
−
1
2
d
2
y
dx
2
a =0. Then u
j
h =yx
j
+hx
j
+Oh
2
. If the problem is solved again
with half the step size, then u
j
1
2
h
= yx
j
+
1
2
hx
j
+Oh
2
. Subtracting these two
results to eliminate the first-order term in h gives
u
j
≡2u
j
1
2
h
−u
j
h
=yx
j
+Oh
2
(12.4.28)