1 SUBROUTINE STRECH(N,P,Q,S)

3 С COMPUTES ONE-DIMENSIONAL STRETCHING FUNCTION,

4 С S = P*ETA + {1.-P)M1.-TANH(QM1.-ETA))/TANH(Q)),

5 С FOR GIVEN CONTROL PARAMETERS, P AND Q. € С

7 DIMENSION S(51)

S AN = N-1

Э DETA l./AN

10 TQI *.1./TANH(Q)

11 С

Хг DOU L 1,N

13 .AL L - 1

14 ETA AL*DETA

15 .DUM s Q*(l. - ETA)

16 DUM * 1. - TANH(DUM)*TQI 17- Sa) P*ETA + (l.-P)*DUM

18 1 CONTINUE

19 RETURN

20 END

Рис. 13.27. Распечатка программы STRECH.

стей (п. 13.3.3) обеспечивает ортогональность сетки вблизи поверхностей Zi и Z4.

Поскольку Zi(r) известно, достаточно провести через каждую точку сетки XS(1,/), YS(1,/) прямую линию, перпендикулярную Zi, до пересечения с Z2.Точками пересечения будутXS(2,У), YS(2,/). Данное построение эквивалентно этапу предиктор, изображенному на рис. 13.17. Поскольку задание узла / определяет значение г на Zi, точки XS(2,/), YS(2, У) имеют то же значение г и ортогональны к Zi в точках XS(1,/), YS(1,/). Расчет ортогональных точек сетки на Z2 осуществляется в подпрограмме SURCH (рис. 13.29). Точки сетки XS(3, У), YS(3,/) на поверхности Z3 строятся так, что они оказываются ортогональными поверхности Z4 в точках XS(4,/), YS(4, У). Расчет этих точек также осуществляется в подпрограмме SURCH.

При заданных координатах на поверхностях от Zi до Z4 алгоритм многих поверхностей в случае четырех поверхностей осуществляется следующим образом [Eiseman, 1979]:

x{J. K)=jl SH(L)XS(L, /), у {J, K)=t SH(L)YS(L, /),

L=l L=\

(13.76)

SH (1) = (1 - 5)2 (1 - a,s), SH (2) = s (1 - sf (a, + 2),

SH(3) = S2(1 - 5)2(2+2), SH(4) = 52(1 - 2(1 ~S)),

a, = 2/(3a - 1), 2 = 2/(2 - 3aJ. (13.78)

1 SUBROUTINE FOIL(INT,T,RAB,X,Y)

3 С SURFACE PROFILE IS NACA-OOT AEROFOIL

4 с IF INT 0, NUMERICALLY INTEGRATE TO OBTAIN SURFACE

5 с COORDINATE

6 С IF INT = 1, INTERPOLATE SURFACE COORDINATES TO OBTAIN

7 С CORRESPONDING (X,Y)

9 DIMENSION A(5),XA(51),RX(51)

10 COMMON RA,XA

11 DATA A/1.4779155,-0.624424,-1.727016,1.384087,-0.489769/

12 PI 3.14159265

13 IF (INT .EQ. DGOTO 2

14 С

15 С NUMERICALLY INTEGRATE TO OBTAIN RA(L) AS A FUNCTION OF XA(L)

16 С

17 RA(1) 0.

18 XA(1> = 0.

19 XA(2) = (A(l)/(1./T - A(2)))**2

20 RA(2) = 0.5*PI*XA(2)

21 DUM = 3.*A(4) + 4.*A(5)*XA(2)

22 DUM = 2.*A(3) + DUM*XA(2)

23 DUH T*(0.5/SQRT(XA(2)) + A(2) + DUM*XA(2))

24 FLP SQRTd, + DUM*DUM)

25 DO 1 L 2,50 26 AL = L

27 LP L + 1

28 XA(LP> 0.02*AL

29 DX = XA(LP) - XA(L)

30 FL = FLP.

31 DUM .3>A(4) + 4.*A(5)*XA{LP)

32 DUM 2.*A(3) + DUM*XA(LP)

33 DUM T*{0.5/SQRT(XA(LP)) + A(2) + DUM*XA(LP))

34 FLP SQRTd. + DUM*DUM)

35 RA(LP> = RA(L) + 0.5MFL + FLP)*DX

36 1 CONTINUE

37 RAB RA{51)

38 RETURN

39 с

40 с INTERPOLATE RA(L) TO OBTAIN X CORR:ESPONDING TO RAB

41 с SUBSEQUENTLY OBTAIN Y FROM ANALYTIC NACA-OOT* PROFILE

42 С

43 2 DO 3 L 2,51

44 IF(RAB .GT. RA(L))GOTO 3

45 LM L - 1

46 X XA(LM) + (XA(L)-XA(LM))*(RAB-RA(LM))/(RA(L)-RA(LM))

47 IF(X .LT. l.OE-06)X l.OE-06

48 DUM A(4) A(5)*X

49 DUM A(3) + DUM*X

50 DUM A(2) + DUM*X

51 Y T*{Ad)*SQRT(X) + DUM*X)

52 RETURN

53 3 CONTINUE

54 VRITE(6,4)RAB,RAd) ,RA(51)

55 4 FORMATC RAB OUTSIDE RANGE,5X, RAB \E10.3,

56 1* RAd)-,E10.3,* RA(51) 4E10.3)

57 RETURN

58 END

SUBROUTINE SURCH(JMAX.S 2.S 3, XS, YS)

GENERATES SURFACES 2 AND 3 TO CREATE ORTHOGONAL ,BOUNDARY GRIDS

DIMENSION XS2(51),YS2{51),XS3(51),YS3(51),XS(4,51),YS(4,51) JMAP JMAX - 1

PRELIMINARY GENERATION OF SURFACES 2 AND 3

DO 1 J = 1,JMAX DXS XS(4,J) - XS(1,J) DYS YS(4,J) - YS(1,J) XS(2,J) = XS(1,J) + S2*DXS YS(2,J) = yS(l,J) + S2*DYS XS(3,J) XS(1,J) + S3*DXS YS(3,J) = YS(1,J) + S3*DYS

1 CONTINUE

PROJECT ORTHOGONALLY FROM SURFACE 1 ONTO SURFACE 2 DO 9 J = 2,JMAP

IF(ABS(XS(1,J+1)-XS(1,J-1)) .GT. 1.0E-06)GOTO 2 EMI 1.0E06*(YS(1,J+1)-YS(1,J-1)) GOTO 3

2 EMI = (YS(1,J+1)-S(1,J-1))/(XS(1,J+1)-XS(1,J-1))

3 IF(ABS(XS(2,J)-XS(2,J-1)) .GT. 1.0E-06)GOTO 4 EM2 = 1.0E+06*(YS(2,J)-YS(2,J-1))

GOTO 5

4 EM2 = (YS(2,J)-YS(2,J-1))/(XS(2,J)-XS(2,J-1))

5 X2 = (EM1*(YS(1,J)-YS(2,J)+EM2*XS{2,J))+XS(1,J))/(1.+EM1*EM2) Y2 = YS(2,J) + EM2*(X2 - XS(2,J))

STJM = SQRT((X2-XS(2,J-1))**2 (Y2-YS (2, J-1)) **2)

SJJM = SQRT((XS(2,J)-XS(2,J-1))**2 + (YS(2,J)-YS(2,J-1))**2)

IF(STJM .LT. SJJM)GOTO 8

IF(ABS(XS(2,J+1)-XS(2,J)) .GT. 1.0E-06)GOTO б EM2 = 1.0E+06*(YS(2,J+1)-YS(2,J)) GOTO 7

6 ЕИ2 = (YS(2,J+1)-YS(2,J))/(XS(2,J1)-XS(2,J))

7 X2 = (£Ml*(YS(l,J)-YS(2,J)+EM2*XS(2,J))+XS(l,J))/(l.+Eni*EM2) Y2 = YS(2,J) + EM2*(X2-XS(2,J))

8 XS2(J) = X2 YS2(J) = Y2

9 CONTINUE

PROJECT ORTHOGONALLY FROM SURFACE 4 ONTO SURFACE 3 DO 17 J 2,JMAP

IF(ABS(XS(4,J+1)-XS(4,J-1)) .GT. 1.0E-06)GOTO 10 EM4 = 1.0E+06*(YS(4,J+1)-YS(4,J-1)) GOTO 11

10 EM4 = (YS(4,J+1) - YS(4,J-1))/(XS(4,J+1) - XS(4,J-1))

11 IF(ABS(XS(3,J) - XS(3,J-1)) .GT. 1.0E-06)GOTO 12 EM3 = 1.0E+06*(YS(3,J)-YS(3,J-1))

GOTO 13