Graebel W.P. Advanced Fluid Mechanics
Подождите немного. Документ загружается.
254 Computational Methods—Ordinary Differential Equations
CVS(2*J,K)=CV1(J)
CVS(2*J+1,K)=CV2(J)
15 END DO
16 END DO
J=MMAX
WRITE(*,*)'------------------------------------------'
WRITE(*,*) J
WRITE(*,*) " ALPHA Q a(J) b(J)"
WRITE(*,*)'------------------------------------------'
DO 17 K=1,28
Q=K
A(2*J-1,K)=CVE(J+1,K)
A(2*J,K)=CVS(J+1,K)
WRITE(*,50)Q,CVE(J+1,K),CVS(J+1,K)
17 END DO
19 END DO
DO 22 L=1,28
IC(L)=L
IB(2*L-1)=L*L
IB(2*L)=L*L
22 END DO
J=0
IXB=0
C WRITE ROWS
WRITE(*,46)J,IXB,IB(1),IB(2),IB(3),IB(4),IB(5),IB(6),IB(7),
& IB(8),IB(9),IB(10),IB(11),IB(12),IB(13),IB(14),IB(15)
WRITE(7,46)J,IXB,IB(1),IB(2),IB(3),IB(4),IB(5),IB(6),IB(7),
& IB(8),IB(9),IB(10),IB(11),IB(12),IB(13),IB(14),IB(15)
WRITE(*,46)J,IXB,IC(1),IC(2),IC(3),IC(4),IC(5),IC(6),IC(7),
& IC(8),IC(9),IC(10),IC(11),IC(12),IC(13),IC(14),IC(15)
WRITE(7,46)J,IXB,IC(1),IC(2),IC(3),IC(4),IC(5),IC(6),IC(7),
& IC(8),IC(9),IC(10),IC(11),IC(12),IC(13),IC(14),IC(15)
DO 24 J=1,14
WRITE(*,45)J,IB(J),A(J,1),A(J,2),A(J,3),A(J,4),A(J,5),A(J,6),
&
A(J,7),A(J,8),A(J,9),A(J,10),A(J,11),A(J,12),A(J,13),A(J,14),
&
A(J,15),A(J,16),A(J,17),A(J,18),A(J,19),A(J,20),A(J,21),
&
A(J,22),A(J,23),A(J,24),A(J,25),A(J,26),A(J,27),A(J,28)
WRITE(7,45)J,IB(J),A(J,1),A(J,2),A(J,3),A(J,4),A(J,5),A(J,6),
&
A(J,7),A(J,8),A(J,9),A(J,10),A(J,11),A(J,12),A(J,13),A(J,14),
&
A(J,15),A(J,16),A(J,17),A(J,18),A(J,19),A(J,20),A(J,21),
&
. A(J,22),A(J,23),A(J,24),A(J,25),A(J,26),A(J,27),A(J,28)
24 END DO
C WRITE COLUMNS
DO 25 J=1,28
WRITE(*,30)J,IB(J),A(1,J),A(2,J),A(3,J),A(4,J),A(5,J),A(6,J),
&
A(7,J),A(8,J),A(9,J),A(10,J),A(11,J),A(12,J)
WRITE(7,30)J,IB(J),A(1,J),A(2,J),A(3,J),A(4,J),A(5,J),A(6,J),
&
A(7,J),A(8,J),A(9,J),A(10,J),A(11,J),A(12,J)
25 END DO
CLOSE(7,STATUS='KEEP')
30 FORMAT(2(I10,","),11(F10.5,","),F10.5)
45 FORMAT(2(I10,","),27(F10.5,","),F10.5)
46 FORMAT(16(I10,","),I10)
50 FORMAT((F6.1,","),F10.4,",",F10.4)
END
DO 16 K=1,28
Q=K
CALL CVA1(3,MMAX,Q,CV1)
CALL CVA1(4,MMAX,Q,CV2)
DO 15 J=1,MMAX/2+1
Program 11.1.2—(Continued)
11.1 Introduction 255
C===================================================================
C
SUBROUTINE CVA1(KD,M,Q,CV)
C
C ============================================================
C Purpose: Compute a sequence of characteristic values of
C Mathieu functions
C Input : M --- Maximum order of Mathieu functions
C q --- Parameter of Mathieu functions
C KD --- Case code
C KD=1 for cem(x,q) ( m = 0,2,4,úúú )
C KD=2 for cem(x,q) ( m = 1,3,5,úúú )
C KD=3 for sem(x,q) ( m = 1,3,5,úúú )
C KD=4 for sem(x,q) ( m = 2,4,6,úúú )
C Output: CV(I) --- Characteristic values; I = 1,2,3,...
C For KD=1,CV(1), CV(2), CV(3),..., correspond to
C the characteristic values of cem for m = 0,2,4,...
C For KD=2,CV(1), CV(2), CV(3),..., correspond to
C the characteristic values of cem for m = 1,3,5,...
C For KD=3,CV(1), CV(2), CV(3),..., correspond to
C the characteristic values of sem for m = 1,3,5,...
C For KD=4,CV(1), CV(2), CV(3),..., correspond to
C the characteristic values of sem for m = 0,2,4,...
C ============================================================
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION G(200),H(200),D(500),E(500),F(500),CV(200)
EPS=1.0D-14
ICM=INT(M/2)+1
IF (KD.EQ.4) ICM=M/2
IF (Q.EQ.0.0D0) THEN
IF (KD.EQ.1) THEN
DO 10 IC=1,ICM
10 CV(IC)=4.0D0*(IC-1.0D0)**2
ELSE IF (KD.NE.4) THEN
DO 15 IC=1,ICM
15 CV(IC)=(2.0D0*IC-1.0D0)**2
ELSE
DO 20 IC=1,ICM
20 CV(IC)=4.0D0*IC*IC
ENDIF
ELSE
NM=INT(10+1.5*M+0.5*Q)
E(1)=0.0D0
F(1)=0.0D0
IF (KD.EQ.1) THEN
D(1)=0.0D0
DO 25 I=2,NM
D(I)=4.0D0*(I-1.0D0)**2
E(I)=Q
25 F(I)=Q*Q
E(2)=DSQRT(2.0D0)*Q
F(2)=2.0D0*Q*Q
ELSE IF (KD.NE.4) THEN
D(1)=1.0D0+(-1)**KD*Q
DO 30 I=2,NM
D(I)=(2.0D0*I-1.0D0)**2
E(I)=Q
30 F(I)=Q*Q
Program 11.1.2—(Continued)
256 Computational Methods—Ordinary Differential Equations
ELSE
D(1)=4.0D0
DO 35 I=2,NM
D(I)=4.0D0*I*I
E(I)=Q
35 F(I)=Q*Q
ENDIF
XA=D(NM)+DABS(E(NM))
XB=D(NM)-DABS(E(NM))
NM1=NM-1
DO 40 I=1,NM1
T=DABS(E(I))+DABS(E(I+1))
T1=D(I)+T
IF (XA.LT.T1) XA=T1
T1=D(I)-T
IF (T1.LT.XB) XB=T1
40 CONTINUE
DO 45 I=1,ICM
G(I)=XA
45 H(I)=XB
DO 75 K=1,ICM
DO 50 K1=K,ICM
IF (G(K1).LT.G(K)) THEN
G(K)=G(K1)
GO TO 55
ENDIF
50 CONTINUE
55 IF (K.NE.1.AND.H(K).LT.H(K-1)) H(K)=H(K-1)
60 X1=(G(K)+H(K))/2.0D0
CV(K)=X1
IF (DABS((G(K)-H(K))/X1).LT.EPS) GO TO 70
J=0
S=1.0D0
DO 65 I=1,NM
IF (S.EQ.0.0D0) S=S+1.0D-30
T=F(I)/S
S=D(I)-T-X1
IF (S.LT.0.0) J=J+1
65 CONTINUE
IF (J.LT.K) THEN
H(K)=X1
ELSE
G(K)=X1
IF (J.GE.ICM) THEN
G(ICM)=X1
ELSE
IF (H(J+1).LT.X1) H(J+1)=X1
IF (X1.LT.G(J)) G(J)=X1
ENDIF
ENDIF
GO TO 60
70 CV(K)=X1
75 CONTINUE
ENDIF
RETURN
END
Program 11.1.2—(Continued)
a search of the World Wide Web provided the needed information in the form of a
downloadable subroutine.
The next program (11.1.3) uses FORTRAN’s capability to deal with complex
numbers, avoiding tedious computations to arrive at the real and imaginary portions of
11.1 Introduction 257
PROGRAM Joukowsky
C=========================================================
C GENERATES A JOUKOWSKY AIRFOIL
C=========================================================
C w = zU* + Ua^2/z + i*GAMMA/2*pi ln z/a
C
C B = GAMMA/4*pi*a*|U| = sin(alpha + sin-1 YCP/a)
C
C dw/dz = U* - U*a^2/z^2 + iGAMMA/2*pi*z
C
C ZP = Z' ZPP = Z''
C
C Z' =Zc' + Z Z'' = Z' + b^2/Z'
C=========================================================
COMPLEX AI,DWDZ,U,USTAR,W,Z,ZCP,ZP,ZPP,ZSP,ZPSP,ZPPSP,ZTE,ZPTE,ZPPTE
PI=3.1415927
A=1.
NPOINT=180
NP=360/NPOINT
AI=CMPLX(0,1.)
DTHETA=2.*PI/NPOINT
OPEN(1, FILE="A:JOUKOWSKY.DAT", STATUS='REPLACE')
WRITE(*,*)
WRITE(*,*) "In this program a and Uinfinity are taken as unity."
WRITE(*,*)
1 WRITE(*,*) "Input Zc':"
WRITE(*,*)
WRITE(*,*) " Thickness increases with Xc'."
WRITE(*,*) " Camber increases with Yc'."
WRITE(*,*) " Xc' = Yc' = 0 stops the program.'"
WRITE(*,*)
WRITE(*,*) " Xc' first. It must lie in -1< Xc' <0. "
READ(*,*) XCP
WRITE(*,*) " Now input Yc'. It must lie in -1< Yc' <1. "
READ (*,*) YCP
WRITE(*,*)
WRITE(*,*) " Enter the angle of attack <degrees>. "
READ(*,*) ALPHA
B=XCP+SQRT(1.-YCP*YCP)
B2=B**2
BETA=ALPHA+ASIN(YCP/A)
GAMMA=4.0*PI*A*SIN(BETA)
R=XCP*XCP+YCP*YCP
SING=(B+XCP)**2+YCP**2-A**2
IF (R.EQ.0.) THEN
STOP
ELSE IF (R.GT.1.0) THEN
WRITE(*,*) " The point Zc' must lie within the unit circle.
&Try again."
GOTO 1
ELSE IF (SING.GE.0) THEN
WRITE(*,*) " The point Z'=-b must lie within the unit circle.
& You will have a singularity in the flow field. Try again."
GOTO 1
ELSE IF (B.LE.0.) THEN
WRITE(*,*) "The parameter b <= 0, try again."
GOTO 1
END IF
Program 11.1.3—Joukowsky Airfoil (program by the author)
258 Computational Methods—Ordinary Differential Equations
WRITE(1,5)
DO I=1, NPOINT+1
THETA=(I-1)*DTHETA
Z=A*CMPLX(COS(THETA),SIN(THETA))
ZP=Z+ZCP
ZPP=ZP+B2/ZP
CALL SPEED(B2,GAMMA,PI,Z,ZP,U,USTAR,UM)
CP=1.-0.5*UM*UM
XPP=REAL(ZPP)
YPP=AIMAG(ZPP)
J=(I-1)*NP
WRITE(*,*) J,XPP,YPP,UM,CP
WRITE(1,100) J,XPP,YPP,UM,CP
END DO
WRITE(*,*)
WRITE(*,*) "Note: The speed calculations are not accurate in "
WRITE(*,*) " the vicinity of the trailing edge. Use"
WRITE(*,*) " l'Hospital's rule here."
WRITE(*,*)
WRITE(*,*) "Press <ENTER> to continue."
READ(*,*)
ZTE=CMPLX(A*COS(ALPHA-BETA),A*SIN(ALPHA-BETA))
ZPTE=ZTE+ZCP
ZPPTE=ZPTE+B2/ZPTE
XTE=REAL(ZPPTE)
YTE=AIMAG(ZPPTE)
ZSP=CMPLX(-A*COS(ALPHA+BETA),-A*SIN(ALPHA+BETA))
ZPSP=ZSP+ZCP
ZPPSP=ZPSP+B2/ZPSP
XSP=REAL(ZPPSP)
YSP=AIMAG(ZPPSP)
WRITE(*,*)
WRITE(*,*) " Trailing edge: ",XTE,YTE
WRITE(*,*) " Stagnation point: ",XSP,YSP
WRITE(*,*) " Xc' , Yc': ",XCP,YCP
WRITE(*,*) " b, beta : ",B,BETA
WRITE(*,*) " Gamma: ",GAMMA
WRITE(1,10) XCP,YCP
WRITE(1,20) B,BETA
WRITE(1,101) XTE,YTE
WRITE(1,102) XSP,YSP
WRITE(1,103) GAMMA
WRITE(*,*)
CLOSE(1)
5 FORMAT("Angle",',',"x''",',',"y''",',',"Speed",',',"Cp")
10 FORMAT(" Xc' = ",F10.7,',',", Yc' = ",F10.7)
20 FORMAT(" b = ",F10.7,',',", beta = ",F10.7," radians")
100 FORMAT (I4,2(',',F10.7),',',F15.7,','F15.7)
101 FORMAT(" Trailing edge : ",2F10.7)
102 FORMAT(" Stagnation point: ",2F10.7)
103 FORMAT(" Gamma: ",F15.7)
END PROGRAM
ALP=ALPHA*PI/180.
U=CMPLX(COS(ALP),SIN(ALP))
USTAR=CMPLX(COS(ALP),-SIN(ALP))
ZCP=CMPLX(XCP,YCP)
Program 11.1.3—(Continued)
11.1 Introduction 259
C===================================================
C SUBROUTINE SPEED
C===================================================
SUBROUTINE SPEED(B2,GAMMA,PI,Z,ZP,U,USTAR,UM)
COMPLEX AI,DWDZ,Z,ZP,U,USTAR
AI=CMPLX(0,1.)
DWDZ=USTAR-U/(Z**2)+AI*GAMMA/(2.*PI*Z)
DWDZ=DWDZ/(1.-B2/(ZP*ZP))
UM=ABS(DWDZ)
RETURN
END SUBROUTINE
Program 11.1.3—(Continued)
complicated complex expressions. The user inputs the three parameters involved in the
Joukowsky transformations described in Chapter 3, as well as the angle of attack of the
airfoil. The program then generates the shape of the airfoil, as well as the velocity and
pressure distributions on its surface.
Program 11.1.4 is an example of the generation of a potential flow about a half-
body consisting of a source in a uniform stream. After specifying source strength and
PROGRAM RANKINE ! RANKINE.FOR
!------------------------------------------------------------------------
! Generates a Rankine half-body in two dimensions
! Potential consists of a source in a uniform stream
! This program is easily generated to the axisymmetric
! half-body with a change in the potential
! from a 2-dimensional source to a three-dimensional source.
! More general shapes can be generated by the addition of more sources
! and/or sinks
! Remember that to generate a closed body the sum of the source
! and sink strengths must be zero.
!------------------------------------------------------------------------
! Data generated here was used in generating data for Figure 2.2.9
! Values used:
! A=m/2*PI*U XSP = Stag point @ -A Half-width = WH = pi*A
! The data in File 1 is easily inserted in any spreadsheet.
! This was made possible by including the commas in FORMAT statement 200.
!------------------------------------------------------------------------
! Configure maximum data size for storage.
DIMENSION S(8),X0(26),Y
0(26),XXS(26),YYS(26),XX(26),YY(11,71)
! Open file for data storage
OPEN(1, FILE='A:RANKINE.DAT', STATUS='UNKNOWN')
A=.125
PI=3.14159
! Generates PSI = 0 streamline on x axis, infinity to source
DO I=1,25
X0(I)=-6.+(12./24.)*(I-1)
Y0(I)=0.
END DO
! Generates PSI = 0 streamline off of x axis, starts at stagnation point
CALL BODYSHAPE(A,PI,XXS,YYS)
! Generates other streamlines
XSTART=-6.
YSTART=0.
Program 11.1.4—Potential flow past a half-body (program by the author)
260 Computational Methods—Ordinary Differential Equations
YS=0.
DX=0.5
DO M=1,25 ! X LOOP
XX(M)=XSTART+(M-1)*DX
X=XX(M)
DO N=1,8 ! PSI LOOP, PSI incremented by 0.2
PSI=0.2*N
CALL SEARCH(PSI,A,X,Y,PSIS,ATD,IS,PI)
YY(N,M)=Y
FER=PSIS-PSI
WRITE(*,100) IS,PSI,FER,X,Y,ATD
END DO
WRITE(*,101)
END DO
WRITE(*,*) " IS PSI PSIS X Y"
DO I=1,25
SM=-YYS(I)
DO J=1,8
S(J)=-YY(J,I)
END DO
WRITE(1,200) X0(I),Y0(I),XXS(I),YYS(I),SM,XX(I),YY(1,I),
& YY(2,I),YY(3,I),YY(4,I),YY(5,I),YY(6,I),YY(7,I),YY(8,I),
& S(1),S(2),S(3),S(4),S(5),S(6),S(7),S(8)
END DO
CLOSE (1)
END PROGRAM
!=====================================================================
SUBROUTINE BODYSHAPE(A,PI,XXS,YYS)
!---------------------------------------------------------------------
! Generates body shape for a 2-D Rankine half-body
! Generates PSI = 0 streamline off of x axis, starts at stagnation point
! Fixes Y, searches for X
!---------------------------------------------------------------------
DIMENSION XXS(26),YYS(26)
WH=A*PI
XSP=-A
DY=WH/25
XXS(1)=XSP
YYS(1)=0.
II=1
XOLD=XSP
Y=DY
X=XOLD-.01
FOLD=Y+A*(ATAN(-X/Y)-.5*PI)
WRITE(*,*) A,WH
! Y loop covers a range of Y's, staring with 0.04*halfwidth
! and continuing to 0.96* half-width.
DO I=1,24
Y=I*DY
X=XOLD
IG=0
WRITE(*,*) "============================="
WRITE(*,*) X,Y,FOLD
! X loop covers a range of X. Since Y is fixed, a starting point is used (XOLD)
! and then incremented a small amount. The search procedure is much as used
! in the subroutine SEARCH, except the roles of X and Y are reversed.
! Comments made in the SEARCH subroutine are true here also regarding ATAN.
100 FORMAT(I6,F10.4,F10.6,3F10.4)
101 FORMAT("============================================")
200 FORMAT(21(F8.4,','),F8.4)
Program 11.1.4—(Continued)
11.1 Introduction 261
!----------------------------------------------------------------------
SUBROUTINE SEARCH(PSI,A,X,YS,PSIS,ATD,IS,PI)
! Generates streamlines
! Fixed X, searches for Y
! said for that).
! The search procedure used is not very elegant, but it works (a lot to be
! A desired value of PSI is inputted, along with a value of x.
! A value of y is guessed at, then a value for PSI based on that is computed.
! The first guess for y is always taken above the desired streamline so it
! gives a value of PSI which is too large.
! The y value is decremented a certain amount dy and the computation is
! repeated, until the error changes sign.
! Then x is incremented and the process repeated.
! Two considerations arise in the computation:
! 1. At and near the stagnation point the slope of the streamline is
! of dy (0.001) are used.
! steep, thus small dy (0.0001) must be used. Later on larger values
! 2. Arctangents are multivalued, and FORTRAN assumes that you want
! values between - pi/2 and + pi/2. Thus in using the ATAN function
! care has to be taken.
!----------------------------------------------------------------------
XXS(25)=6.
YYS(25)=WH
WRITE(*,*) " X Y FOLD"
WRITE(*,*) " I J X Y F ATD"
100 FORMAT(F10.5,','F10.5)
200 FORMAT(2I5,3F10.5,F12.5)
RETURN
END SUBROUTINE
DO J=1,1500
DX=.0001
IF(X.GT.0.234) DX=.0009
IF(I.EQ.24) DX=0.0015
X=DX*(J-1)+XOLD
AT=ATAN(-X/Y)
F=Y+A*(AT-.5*PI)
! Test for convergence.
FTEST=FOLD*F
ATD=AT*180./PI
IF((FTEST.LE.0.).AND.(IG.EQ.0)) THEN
WRITE(*,200)I,J,X,Y,F,ATD
II=I+1
XXS(II)=X
YYS(II)=Y
XS=X
YS=Y
IG=1
ENDIF
FOLD=F
END DO
XOLD=XS-.01
FOLD=YS+A*(ATAN(-(XS-.1)/YS)-.5*PI)
END DO
WH=A*PI
XSP=-A
DY=0.0001
IF(X.GE.-.6) DY=.001
Program 11.1.4—(Continued)
262 Computational Methods—Ordinary Differential Equations
FOLD=Y+A*(ATAN(-X/(Y-DY))-.5*PI)
IG=0
Y=PSI-.015
! This is the Y loop, where Y is continually incremented.
DO I=1,1500
Y=Y+DY
PSIOLD=PSI-.05
AT=ATAN(-X/Y)
F=Y+A*(AT-.5*PI)
FTEST=F-PSI
! Test to see if the desired streamline has been crossed.
! If it has, linear interpolation is used to improve result.
IF((FTEST.GE.0.).AND.(IG.EQ.0)) THEN
IS=I
FF=(F-PSIOLD)/(PSI-PSIOLD)
PSIS=F
YS=Y
IG=1
ATD=AT*180./PI
ENDIF
PSIOLD=F
END DO
RETURN
END SUBROUTINE
Program 11.1.4—(Continued)
uniform stream velocity, the body shape and surrounding streamlines are found. The
computation is complicated by the fact that on a streamline it is not possible to solve
directly for y as a function of x. An iteration procedure must be used. The appropriate
source strength on the panel is computed so as to make the velocity tangent to the body.
11.2 Numerical Calculus
Many schemes are available for computing derivatives and integrals of functions,
depending on the degree of accuracy needed. Many methods are based on the Taylor
series expansion of a function—that is, for yx =fx,on
yx =
n=0
1
n!
f
n
x
0
x −x
0
n
(11.2.1)
Here f
n
x
0
stands for the nth order derivative of f evaluated at x
0
.Ify is known at
neighboring points and knowledge of its derivative is required, then to lowest order in
the point spacing, write yx +x = yx +f
xx +Ox
2
, giving
f
x =
yx +x −yx
x
+Ox (11.2.2)
If greater accuracy is needed, start with
yx +x = yx +f
1
xx +
1
2
f
2
xx
2
+Ox
3
(11.2.3a)
yx −x = yx −f
1
xx +
1
2
f
2
xx
2
+Ox
3
(11.2.3b)
11.2 Numerical Calculus 263
Subtracting one of the two equations from the other gives
f
1
x =
yx +x −yx −x
2x
+Ox
2
(11.2.4)
Further accuracy can be found by including more terms and more points in the Taylor
expansion. Some examples follow:
Central differences
dfx
dx
≈
1
zx
fx +x −fx −x
+O
x
2
Forward differences
dfx
dx
≈
1
x
fx +x −fx
+O
x
(11.2.5)
Backward differences
dfx
dx
≈
1
x
fx −fx −x
+O
x
Newton’s method is frequently used for finding the roots of an equation of the form
fx = 0. (Other methods can be found in Press, Flannery, Teukolsky, and Vetterling
(1986), and in Carnahan, Luther, and Wilkes (1969).) Newton suggested making an
initial guess and then repeatedly improving it according to the formula
x
n+1
=x
n
−
fx
n
f
x
n
(11.2.6)
the prime denoting the derivative of f . This is merely a rearrangement of equation
(11.2.2). Providing that the first guess is reasonably close to the root sought, the iteration
proceeds quickly to the desired root.
Newton’s rule can be used for the computation of the stability curves for
Couette-Rayleigh instability. First, to simplify the calculations, introduce =
−
c
2
0
−a
2
/a
2
. Then
c
2
n
=a
2
1−
cos
2n
3
+i sin
2n
3
n=0 1 2
Since we expect that >1,
c
0
=ai
√
−1,
c
1
=a
1+
1
2
−i
√
3
2
=a
1+ +
2
e
−i/2
c
2
=a
1+
1
2
+i
√
3
2
=a
1+ +
2
e
i/2
with = tan
−1
√
3
2 +
To start the computation, use the known minimum critical values for even distur-
bances to compute the Rayleigh number, and then move first to slightly smaller values
of the wave number, understanding that the corresponding Rayleigh number is always
slightly larger than the previous value. Use Newton’s method repeatedly to guess at the
Raleigh number for this new wave number. (The derivative is computed using equation
(11.2.2).) Repeat the procedure for wave numbers larger than critical, and finally the