谁有fortran语言的坐标转换程序?

lihong

我现在准备作一个水翼绕流计算的程序,急用fortran语言的坐标转换程序,从直角坐标转换至一般曲线坐标。若哪位有这东东，万望帮我一把。多谢了！

gongjianwei0

唉 高手都去哪了？？

legolasGG

For generalized curvilinear grid, please try Fengyan Shi';s coastal Grid
  http://www.coastal.udel.edu/~fyshi
while, it';s a matlab code,
  you may also go to Steven Owen';s mesh survey and see if u can find something:
   http://www.andrew.cmu.edu/user/sowen/softsurv.html

caiqd

我没有明白你的问题，是要生成一个贴体网格吗？如果是，应该有许多软件可以做到。

houston

在此转一下别人的贴子
      PROGRAM MAIN
      DIMENSION X(509) , Y(509) , R(509) , THETA(509),
     1 R1(509), THETA1(509)
      DIMENSION XG(509,509),YG(509,509)
      COMPLEX Z1,Z2,Z3,Z4,Z5,Z6,Z7
C
C ALGEBRAIC O- GRID GENERATOR BASED ON
C KARMAN-TREFFZ MAPPING
C
C ******************************************************
C
C EXPALNATION OF INPUT QUANTITIES
C
C FSYM : FLAG THAT INDICATES IF AIRFOIL IS SYMMETRIC
C FSYM > 0 MEANS AIRFOIL IS SYMMETRIC
C FSYM = 0 MEANS CAMBERED AIRFOIL
C NU NO. OF DATA POINTS WHERE UPPER SURFACE
C ORDINATES ARE SPECIFIED
C NL NO. OF DATA POINTS WHERE LOWER SURFACE
C ORDINATES ARE SPECIFIED
C X,Y INPUT X,Y COORDINATES OF THE AIRFOIL
C
C XN,YN COORDINATES OF A SINGULAR POINT INSIDE THE
C AIRFOIL NEEDED BY THE KARMAN-TREFFZ MAPPING
C
C IMAX TOTAL NUMBER OF I- NODES ON THE BODY-FITTED
C GRID. MUST BE LESS THAN 81
C I = IMAX-1 IS UPPER-SURFACE TRAILING EDGE
C I = 2 IS LOWER-SURFACE TRAILING EDGE
C
C JMAX TOTAL NUMBER OF J- NODES ON THE BODY FITTTED
C GRID MUST BE LESS THAN 31
C NOTE; J = 2 IS AIRFOIL SURFACE
C
C ************************************************************
      READ (5,1)
    1 FORMAT(1X)
      READ (5,1)
      READ (5,*) FSYM ,NU,NL , IMAX , JMAX
      IF(IMAX.GT.500) STOP
      IF(JMAX.GT.500) STOP
      IF(FSYM.GT.0) NL = NU
      N = NU + NL - 1
C
C READ UPPER SURFACE ORDINATES
C
      READ (5,1)
      DO 1000 I = NL , N
1000 READ (5,*) X(I),Y(I)
C
      IF(FSYM.GT.0.0) THEN
C
C SYMMETRIC AIRFOIL... FLIP SIGN TO GET
C LOWER SURFACE COORDINATES
C
      DO 2000 I = NL , 1 , -1
      X(I) = X(N+1-I)
      Y(I) = - Y(N+1-I)
2000 CONTINUE
      ELSE
C
C CAMBERED AIRFOIL... READ LOWER SURFACE DATA TOO
C
      READ (5,1)
      DO 2010 I = NL , 1 , -1
2010 READ (5,*) X(I) , Y(I)
      ENDIF
      READ (5,1)
C
C READ LEADING ENDGE SINGULAR POINT
C
      READ (5,*) XN,YN
      SLOP1 = ATAN2((Y(1)-Y(2)),(X(1)-X(2)))
      SLOP2 = - ATAN2((Y(N)-Y(N-1)),(X(N)-X(N-1)))
      EPS = ABS(SLOP1 + SLOP2)
      PI = 3.1415926535
      POWER = 2.0 - EPS / PI
      XT = 0.5 * (X(1) + X(N))
      YT = 0.5 * (Y(1) + Y(N))
      Z1 = CMPLX(XN,YN)
      Z2 = CMPLX(XT,YT)
      Z7 = 0.5 * (Z1+Z2)
      U = 1.0
      V = 0.0
      ANGL = 8.0 * ATAN(1.0)
      POWER = 1.0 / POWER
      DO 3000 I = 1 , N
      Z3 = CMPLX(X(I),Y(I))
      Z4 = (Z3 - Z2) / (Z3-Z1)
      Z4 = CLOG(Z4) * CMPLX(POWER,0.0)
      Z4 = CEXP(Z4)
      Z5 = CMPLX(1.0,0.0)
      Z4 = (Z5+Z4) * Z7 / (Z5-Z4)
      X4 = REAL(Z4)
      Y4 = AIMAG(Z4)
      R(I) = SQRT(X4*X4+Y4*Y4)
      A = X4
      B = Y4
      ANGL = ANGL + ATAN2((U*B-V*A),(U*A+V*B))
      U = A
      V = B
C 计算得到机翼表面离散点对应的复平面半径和幅角
      WRITE (6,*) ';I= ';,I,'; MAG=';,R(I),';THETA=';,ANGL
      THETA(I) = ANGL
3000 CONTINUE
C
C LINEARLY INTERPOLATE R VS. THETA TO GET R AT UNIFORMLY
C SPACED THETA LOCATIONS
C
      DTHETA = 8. * ATAN(1.0) / (IMAX-3)
      DO 3001 I = 3 , IMAX - 2
      THETA1(I) = (I-2) * DTHETA
C 我在这里加入了尾缘附近的处理，直接让离散数据点角度范围外的
C 点的半径等于最近的离散点的半径，不知道合适与否
      IF( THETA1(I) .LT. THETA(N) ) THEN
        R1(I) = R(N)
        WRITE (6,*) "<< ",I ,THETA1(I) , R1(I)
        GOTO 3001
      ELSE IF( THETA1(I) .GT. THETA(1) ) THEN
    R1(I) = R(1)
        WRITE (6,*) ">> ", I ,THETA1(I) , R1(I)
        GOTO 3001
      ENDIF
C 以上是加入的代码
      
      DO 3010 L = 1 , N-1
      IF((THETA(L) - THETA1(I)) *
     1 (THETA(L+1) - THETA1(I)).LE.0.0) THEN
          DRDT = (R(L+1) - R(L)) / (THETA(L+1)-THETA(L))
          R1(I) = R(L) + DRDT * (THETA1(I) - THETA(L))
          WRITE (6,*) "** ",I ,THETA1(I) , R1(I)
          GOTO 3011
      END IF
3010 CONTINUE
      WRITE (6,*) ';INTERPOLATION FAILED....';
      STOP
3011 CONTINUE
3001 CONTINUE
      R1(2) = 0.5 * (R(1) + R(N))
      R1(IMAX-1) = R1(2)
      R1(IMAX) = R1(3)
      R1(1) = R1(IMAX-2)
      THETA1(2) = 0.0
      THETA1(IMAX-1) = 8. * ATAN(1.0)
      THETA1(1) = - DTHETA
      THETA1(IMAX)= THETA1(IMAX-1) + DTHETA
C
C GENERATE GRID
C
      DR = 0.99 / (JMAX-2)
      DO 4000 J = 2 , JMAX
      ETA = (J-2) * DR
      ETA1 = 1. + (5.*ETA)**2.7
C ETA1 = 1. / (1. - ETA)
      DO 4000 I = 1 , IMAX
      X4 = R1(I) * ETA1 * COS(THETA1(I))
      Y4 = R1(I) * ETA1 * SIN(THETA1(I))
      Z4 = CMPLX(X4,Y4)
      Z4 = (Z4-Z7) / (Z4+Z7)
      Z4 = CMPLX(1.0/POWER,0.0) * CLOG(Z4)
      Z4 = CEXP(Z4)
      Z4 = (Z2 - Z4 * Z1) / (CMPLX(1.0,0.0) - Z4)
      X4 = REAL(Z4)
      Y4 = AIMAG(Z4)
      IF(J.EQ.2) WRITE (6,*) I , X4 , Y4
      XG(I,J) = X4
      YG(I,J) = Y4
4000 CONTINUE
C
C EXTRAPOLATE GRID TO THE INTERIOR OF THE AIRFOIL
C TO OBTAIN COORDINATES AT J = 1
C
      DO 4001 I = 1 , IMAX
      XG(I,1) = 2. * XG(I,2) - XG(I,3)
4001 YG(I,1) = 2. * YG(I,2) - YG(I,3)
      OPEN(UNIT=9,FILE=';OG.OUT';,STATUS=';UNKNOWN';)
c WRITE (9,*) IMAX,JMAX
      WRITE(9,*) ';variables="x","y"';
      WRITE(9,*) ';Zone T="O-grid" I=';,IMAX-2,';,';,';J=';,JMAX-3, ';f=point';
      DO 911 J=2,JMAX
      IF( (J .NE. 3) .AND. (J .NE. 4) ) THEN
           WRITE (9,4002) (XG(I,J),YG(I,J),I=2,IMAX-1)
      ENDIF
911 CONTINUE
      STOP
4002 FORMAT(2F10.4)
      END

iamfirst

推荐一本书，陈景仁的《湍流模型及有限分析法》，在坐标变换那一章附有一个坐标转换的源程序，介绍得非常详细，可以去参考一下！

wangye03219

非常感谢 我也正在做这样的程序