Thompson, Warsi, Mastin - Numerical Grid Generation (523190), страница 43
Текст из файла (страница 43)
Note that (2) and (3) are totally unreasonable.4. Generate an algebraic grid for a circular simply-connected region by (1) unidirectionalinterpolation, (2) tensor product interpolation, and (3) transfinite interpolation. Note that hereonly (3)gives a reasonable grid.
Compare linear and Hermite interpolation for (3).5. Repeat Exercise 4 with a triangular boundary.6. Using the boundary configuration of Exercise 3, but with a hyperbolic tangent pointdistribution on the right-hand boundary of the physical region with smaller spacing at thecenterline than at the top and bottom. Compare algebraic grids generated using (1) linearinterpolation between the inner and outer boundaries, (2) nonlinear interpolation, based onthe hyperbolic tangent, between the inner and outer boundaries, (3) transfinite interpolationwith linear blending functions, and (4) transfinite interpolation using the boundary pointdistribution (in terms of relative arc length) as the blending functions. Note that only (2) and(4) preserve the boundary point distribution in the field.7.
Generate an algebraic grid for a square inside a rectangle using linear interpolationbetween the inner and outer boundaries. Note the propagation of the boundary slopediscontinuities into the field. Generate a grid from an elliptic generation system for the sameboundary point distribution and note the difference.8. Generate an algebraic grid for a square inside a circle using linear interpolation betweenthe inner and outer boundaries. Show that it is possible to position the points on the circlesuch that the grid overlaps the corners of the square. Generate a grid from an ellipticgeneration system for the same boundary point distribution and note the difference.3. Listing of Routine for Computer ExercisesSUBROUTINE INTERPPARAMETER(NI=20,NJ=20,N=5)COMMON/COORD/X(NI,NJ),Y(NI,NJ)COMMON/CONST/CHOICE,IMAX,JMAX,NA,DS1,DS2COMMON/ATTR/IAL(N),IAX(N),IAY(N),JAL(M),JAX(N),JAY(N)COMMON/COEF/AI(N),BI(N),CI(N),DI(N),AJ(N),BJ(N)COMMON/COEF/CJ(N),DJ(N)DIMENSION P(NI,NJ),Q(NI,NJ),XX(0:NI,O:NJ),YY(0:MI,O:NJ)DIMENSION X1(NI,NJ),X2(NI,NJ),Y1(NI,NJ),Y2(NI,NJ)INTEGER CHOICECCCCCCCCCCCCCCCCCCCCCCCBOUNDARY INTERPOLATIONXYIMAXJMAXNADSlDS2ATTRCOEFX ARRAY OF XI-ETA COORDINATEY ARRAY OF XI-ETA COORDINATEMAX.
NUMBER OF GRID IN XI AXISMAX. NUMBER OF GRID IN ETA AXISMAX. NUMBER OF ATTRACTIONSSPECIFIED LENGTH OF INITIAL INTERVALSPECIFIED LENGTH OF FINAL INTERVALARRAY OF ATTRACTION TO LINES/POINTSARRAY OF COEFFICIENT FOR ATTRACTIONCHOICE12345678VERTICAL INTERPOLATIONHORIZONTAL INTERPOLATIONTENSOR PRODUCT INTERPOLATIONTRANSFINITE INTERPOLATIONHERMITE CUBIC INTERPOLATIONHYPERBOLIC TANGENT INTERPOLATIONELLIPTIC GRID GENERATION ( SOR ITERATION )ATTRACTION TO COORDINATESIF(CHOICE.EQ.1)IF(CHOICE.EQ.2)IF(CHOICE.EQ.3)IF(CHOICE.EQ.H)IF(CHOICE.EQ.5)IF(CHOICE.EQ.6)IF(CHOICE.EQ.7)IF(CHOICE.EQ.S)CCCGOGOGOGOGOGOGOGOTOTOTOTOTOTOTOTO100200300400500600700800**** VERTICAL INTERPOLATION ****100C110CCDO 110 I 1,IMAXDO 110 J 1,JMAXRJ1=FLOAT(JMAX-J)/FLOAT(JMAX-1)RJ2=FLOAT(J-1)/FLOAT(JMAX-1)*** ( EQ. 8-1 )X(I,J)=RJ1*X(I,l)+RJ2*X(I,JMAX)Y(I,J) RJl*Y(I.1)+RJ2*Y(I,JMAX)RETURN**** HORIZONTAL INTERPOLATION ****C200C210CCCDO 210 I=l,JMAXDO 210 J=1,IMAXRI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)RI2=FLOAT(I-1)/FLOAT(IMAX-1)*** ( EQ.
8-1 )X(I,J)=RI1*X(1,J)+RI2*X(IMAX,J)Y(I,J)=RI1*Y(1,J)+RI2*Y(IMAX,J)RETURN**** TENSOR PRODUCT INTERPOLATION ****300C310CCCDO 310 I=1,IMAXDO 310 J=l,JMAXRI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)RI2=FLOAT(I-1)/FLOAT(IMAX-1)RJl=FLOAT(JMAX-J)/FLOAT(JMAX-1)RJ2=FLOAT(J-1)/FLOAT(JMAX-1)*** ( EQ. 8-69 )X(I,J)=RI1*RJ1*X(1,1)+RI1*RJ2*X(1,JMAX)*+RI2*RJ1*X(IMAX,1)+RI2*RJ2*X(IMAX,JMAX)Y(I,J)=RI1*RJ1*Y(1,1)+RI1*RJ2*Y(1,JMAX)*+RI2*RJ1*Y(IMAX,1)+RI2*RJ2*Y(IMAX,JMAX)CONTINUERETURN**** TRANSFINITE INTERPOLATION ****400410420C430CCCDO 410 I=l,IMAXDO 410 J=1,JMAXRI1=FLOAT(I-1)/FLOAT(IMAX-1)RI2=FLOAT(IMAX-I)/FLOAT(IMAX-1)X1(I,J)=RI1*X(IMAX,J)+RI2*X(1,J)Yl(I,J)=RI1*Y(IMAX,J)+RI2*Y(1,J)DO 420 I=1,IMAXDO 420 J=1,JMAXRJ1=FLOAT(J-1)/FLOAT(JMAX-1)RJ2=FLOAT(JMAX-J)/FLOAT(JMAX-1)X2(I,J)=RJ1*(X(I,JMAX)-X1(I,JMAX))+RJ2*(X(I,1)-X1(I,1))Y2(I,J)=RJ1*(Y(I,JMAX)-Y1(I,JMAX))+RJ2*(Y(I,l)-Y1(I,1))*** ( EQ.
8-73 )DO 430 I=1,IMAXDO 430 J=1,JMAXX(I,J)=X1(I,J)+X2(I,J)Y(I,J)=Y1(I,J)+Y2(I,J)IF(CHOICE.NE.4) GO TO 740RETURN***500510520HERMITE CUBIC INTERPOLATION (ORTHOGONAL BOUNDARY) ***DO 510 I=1,IMAXDO 510 J=1,JMAXXX(I,J)=X(I,J)YY(I,J)=Y(I,J)DO 520 J=1,JMAXXX(O,J)=XX(IMAX-1,J)YY(O,J)=YY(INAX-1,J)XX(IMAX+1,J)=XX(2,J)YY(IMAX+1,J)=YY(2,J)DO 530 I=1,IMAXDO 530 J=1,JMAXRJJ=FLOAT(J-1)/FLOAT(JMAX-1)*** ( EQ. 8-6 a and b, n=2 )PHI1=(1.+2.*RJJ)*(1.-RJJ)*(1.-RJJ)PHI2=(3.-2.*RJJ)*RJJ*RJJPSI1=(1.-RJJ)*(1.-RJJ)*RJJPSI2=(RJJ-1.)*RJJ*RJJCCCC** CAL.
NORMAL DERIV. **CC530540CCCXXI1=.5*(XX(I+1,1)-XX(I-1,1))XXI2=.5*(XX(I+1,JMAX)-XX(I-1,JMAX))YXI1=.5*(YY(I+1,1)-YY(I-1,1))YXI2=.5*(YY(I+1,JMAX)-YY(I-1,JMAX))UNIT1=SQRT(XXI1*XXI1+YXI1*YXI1)UNIT2=SQRT(XXI2*XXI2+YXI2*YXI2)*** ( EQ. 3-108 )XNl=-YXI1/UNIT1*DSlXN2=-YXI2/UNIT2*DS2YN1=XXI1/UNIT1*DSlYN2=XXI2/UNIT2*DS2*** ( EQ. 8-5 )XX(I,J)=PHI1*XX(I,l)+PHI2*XX(I,JMAX)+PSI1*XN1+PSI2*XN2YY(I,J)=PHI1*YY(I,1)+PHI2*YY(I,JMAX)+PSI1*YN1+PSI2*YN2DO 540 I=1,IMAXDO 540 J=1,JMAXX(I,J)=XX(I,J)Y(I,J)=YY(I,J)RETURN**** HYPERBOLIC TANGENT SPACING INTERPOLATION ****600CC610620630C650660CCCCTOL=1.0E-10*** ( EQ.
8-49, 50 and 51A=SQRT(DS2/DS1)B=1./(FLOAT(JMAX-1)*SQRT(DS1*DS2))*** INITIAL GUESS BY SERIES EXPANSIONDELTA=SQRT(6.*(B-1.))DO 610 IT=1,20RESID=SINH(DELTA)/(DELTA*B)-1.IF(ABS(RESID)LT.TOL) GO TO 630CALL AITKEN(DELTA,RESID,DELTO,RO,RSO)PRINT 620, RESID,DELTA,IT-1FORMAT(//, 5X, ’DELTA IS NOT CONVERGE ?’, 5X, 2E15.5,*5X, I3, //)GO TO 660CONTINUE*** ( EQ. 8-52, 53 and 54 )DO 650 I=1,IMAXDO 650 J=2,JMAX-1RATIO FLOAT(J-1)/FLOAT(JMAX-1)U=.5*(1.+TANH(DELTA"(RATIO-.5))/TANH(.5*DELTA))S=U/(A+(1.-A)*U)X(I,J)=X(I,1)+(X(I,JMAX)-X(I,1))*SY(I,J)=X(I,1)+(Y(I,JMAX)-Y(I,1))*SRETURN**** ELLIPTIC GRID GENERATION ( SOR ITERATION ) ******* CAL. P AND Q ON THE BOUNDARY **700DO 710 I=1,IMAXDO 710 J=1,JMAXC710CCCCXXI=.5*(X(I+1,J)-X(I-1,J))XXIXI=X(I+1,J)-2.*X(I,J)+X(I-1,J)XETA=.5*(X(I,J+1)-X(I,J-1))XETA2=X(I,J+1)-2.*X(I,J)+X(I,J-1)YXI=.5*(X(I+1,J)-Y(I-1,J))XXIXI=Y(I+1,J)-2.*X(I,J)+Y(I-1,J)YETA=.5*(Y(I,J+1)-Y(I,J-1))YETA2-Y(I,J+1)-2.*Y(I,J)+Y(I,J-1)IF(ABS(XETA2).LT.10E-3) XETA2=0.IF(ABS(YETA2).LT.10E-3) YETA2=0.RXI2=XXI*XXI+YXI*YXIRETA2=XETA*XETA+YETA*YETA*** ( EQ.
8-70 )P(I,J)=(XXI*XXIXI+XXI*YXIXI)/RXI2Q(I,J)=(XETA*XETA2+YETA*YETA2)/RETA2CONTINUE**720CCCDO 720 I=1,IMAXDO 720 J=l,JMAXRJ1=FLOAT(JMAX-J)/FLOAT(JMAX-1)RJ2=FLOAT(J-1)/FLOAT(JMAX-1)RI1=FLOAT(IMAX-I)/FLOAT(IMAX-1)RI2=FLOAT(I-1)/FLOAT(IMAX-1)P(I,J)=RJl*P(I,l)+RJ2*P(I,JMAX)Q(I,J)=RI1*Q(1,J)+RI2*Q(IMAX,J)CONTINUE**740CCCCINTERPOLATE P AND Q BETWEEN BOUNDARY **P : VERTICAL, Q : HORIZONTALINITIAL GUESS WITH TRANSFINITE INTERPOLATION **GO TO 400CONTINUE*** ITERATION ( SOR ) ***ITMAX=200TOL=10.E-5W=1.8DO 760 IT=1,ITMAXERRX=0.ERRY=0.DO 750 J=2,JMAX-1DO 750 I=2,IMAX-1XXI=.5*(X(I+1,J)-X(I-1,J))YXI=.5*(Y(I+1,J)-Y(I-1,J))XXIXI=X(I+1,J)+X(I-1,J)YXIXI=Y(I+1,J)+Y(I-1,J)XETA=.5*(X(I,J+1)-X(I,J-1))YETA .5*(Y(I,J+1)-Y(I,J-1))XXIETA=.25*(X(I+1,J+1)-X(I+1,J-1)-X(I-1,J+1)+X(I-1,J-l))YXIETA=.25*(Y(I+1,J+1)-Y(I+1,J-1)-Y(I-1,J+1)+Y(I-1,J-1))XETA2=X(I,J+1)+X(I,J-1)YETA2=Y(I,J+1)+Y(I,J-1)*** ( EQ.
6-18 and 6-20 )G11=XXI*XXI+YXI*YXIG22=XETA*XETA+YETA*YETAG12=XXI*XETA+YXI*YETAXTEMP=.5*(G22*(P(I,J)*XXI+XXIXI)+G11*(Q(I,J)*XETA+XETA2)*-2.*G12*XXIETA)/(G11+G22)750760770780CCCYTENP=.5*(G22*(P(I,J)*YXI+YXIXI)+G11*(Q(I,J)*YETA+YETA2)*-2.*Gf2*YXIETA)/(Gll+G22)XTENP=W*XTEMP+(1.-W)*X(I,J)YTEMP=W*YTEMP+(1.-W)*Y(I,J)ERRX=AMAXO(ERRX,ABS(XTEMP-X(I,J)))ERRY=AMAXO(ERRY,ABS(YTEMP-Y(I,J)))X(I,J)=XTEMPY(I,J)=YTEMPCONTINUEIF(ERRX.LT.TOL.AND.ERRY.LT.TOL) GO TO 78OCONTINUEPRINT 770,ERRX,ERRY,IT-1FORMAT(//, 5X, ’X AND Y ARE NOT CONVERGE ?’, 2E15.5,*5X, I5, //)CONTINUEIF(CHOICE.EQ.8) GO TO 830RETURN**** ATTRACTION TO COORDINATE LINE/POINT ****800810C820830DO 810 I=1,IMAXDO 810 J=1,JMAXP(I,J)=0.Q(I,J)=0.DO 820 NS=1,NADO 820 I=1,IMAXDO 820 J=1,JMAXXL=FLOAT(I-IAL(NS))XI=FLOAT(I-IAX(NS))XJ=FLOAT(J-IAY(NS))YL=FLOAT(J-JAL(NS))YI=FLOAT(I-JAX(NS))YJ=FLOAT(J-JAY(NS))*** ( EQ.
6-30 )P(I,J)=P(I,J)-AI(NS)*(XL/ABS(XL))*EXP(-CI(NS)*ABS(XL))*-BI(NS)*(XI/ABS(XI))*EXP(-DI(NS)*SQRT(XI*XI+XJ*XJ))Q(I,J)=Q(I,J)-AJ(NS)*(YL/ABS(YL))*EXP(-CJ(NS)*ABS(YL))*-BJ(NS)*(YJ/ABS(YJ))*EXP(-DJ(NS)*SQHT(YI*YI+YJ*YJ))CONTINUEGO TO 400CONTINUERETURNEND4. Examples for Computer ExercisesREFERENCES1. Thompson, Joe F., Warsi, Z. U.