课题组有一个较老的计算程序,使用的是贴体坐标系。控制方程也需要进行转换在贴体坐标系下相应的形式。如下
里面的alpha和beta等参数定义如下
XI和ETA是贴体坐标系的两个方向,我在求以上jacobi,beta,gama,alpha时遇到一些问题,请给位老看掌掌眼。老程序应用的网格中beta等参数没有负的,我的新网格会出现负的。但是从公式来看,确实存在负数的可能性。不得其解。我的网格是
参数求解的代码如下:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 计算JCB,ALPHA,BETA,GAMA,XXI,XETA,YXI,YETA等等
!
! 在计算的过程中,边界采用前差,内点采用中心差分
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PROGRAM MAIN
IMPLICIT NONE
INTEGER I,J,INDCOS
INTEGER NIC,NJC
INTEGER,PARAMETER :: IT=206,JT=92
INTEGER LIMIT(IT,JT),KONES(IT,JT)
REAL X(IT,JT),Y(IT,JT),JCB(IT,JT),ALPHA(IT,JT),BETA(IT,JT)
REAL GAMA(IT,JT),XXI(IT,JT),XETA(IT,JT),YXI(IT,JT),YETA(IT,JT)
!!!!!!!!!!!!网格文件的打开!!!!!!!!!!!!!
OPEN(10,FILE='ZhongXinJieMianErWei.DAT',STATUS='OLD')
READ(10,*)NIC,NJC
DO J=1,NJC
DO I=1,NIC
READ(10,*)X(I,J),Y(I,J),LIMIT(I,J),KONES(I,J)
END DO
END DO
CLOSE(10)
!!!!!!!!!!!!!!JCB等的计算!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!拐角点(1,1)
XXI(1,1)=X(2,1)-X(1,1)
XETA(1,1)=X(1,2)-X(1,1)
YXI(1,1)=Y(2,1)-Y(1,1)
YETA(1,1)=Y(1,2)-Y(1,1)
!!!!!拐角点(1,NJC)
XXI(1,NJC)=X(2,NJC)-X(1,NJC)
XETA(1,NJC)=X(1,NJC)-X(1,NJC-1)
YXI(1,NJC)=Y(2,NJC)-Y(1,NJC)
YETA(1,NJC)=Y(1,NJC)-Y(1,NJC-1)
!!!!!拐角点(NIC,NJC)
XXI(NIC,NJC)=X(NIC,NJC)-X(NIC-1,NJC)
XETA(NIC,NJC)=X(NIC,NJC)-X(NIC,NJC-1)
YXI(NIC,NJC)=Y(NIC,NJC)-Y(NIC-1,NJC)
YETA(NIC,NJC)=Y(NIC,NJC)-Y(NIC,NJC-1)
!!!!!拐角点(NIC,1)
XXI(NIC,1)=X(NIC,1)-X(NIC-1,1)
XETA(NIC,1)=X(NIC,2)-X(NIC,1)
YXI(NIC,1)=Y(NIC,1)-Y(NIC-1,1)
YETA(NIC,1)=Y(NIC,2)-Y(NIC,1)
!!!!!左边界点
DO J=2,NJC-1
XXI(1,J)=X(2,J)-X(1,J)
YXI(1,J)=Y(2,J)-Y(1,J)
XETA(1,J)=(X(1,J+1)-X(1,J-1))/2.0
YETA(1,J)=(Y(1,J+1)-Y(1,J-1))/2.0
END DO
!!!!!上边界点
DO I=2,NIC-1
XXI(I,NJC)=(X(I+1,NJC)-X(I-1,NJC))/2.0
YXI(I,NJC)=(Y(I+1,NJC)-Y(I-1,NJC))/2.0
XETA(I,NJC)=X(I,NJC)-X(I,NJC-1)
YETA(I,NJC)=Y(I,NJC)-Y(I,NJC-1)
END DO
!!!!!下边界点
DO I=2,NIC-1
XXI(I,1)=(X(I+1,1)-X(I-1,1))/2.0
YXI(I,1)=(Y(I+1,1)-Y(I-1,1))/2.0
XETA(I,1)=X(I,2)-X(I,1)
YETA(I,1)=Y(I,2)-Y(I,1)
END DO
!!!!!右边界点
DO J=2,NJC-1
XXI(NIC,J)=X(NIC,J)-X(NIC-1,J)
YXI(NIC,J)=Y(NIC,J)-Y(NIC-1,J)
XETA(NIC,J)=(X(NIC,J+1)-X(NIC,J-1))/2.0
YETA(NIC,J)=(Y(NIC,J+1)-Y(NIC,J-1))/2.0
END DO
!!!!!内部点
DO J=2,NJC-1
DO I=2,NIC-1
XXI(I,J)=(X(I+1,J)-X(I-1,J))/2.0
XETA(I,J)=(X(I,J+1)-X(I,J-1))/2.0
YXI(I,J)=(Y(I+1,J)-Y(I-1,J))/2.0
YETA(I,J)=(Y(I,J+1)-Y(I,J-1))/2.0
END DO
END DO
!!!!!!几何参数的计算
DO J=1,NJC
DO I=1,NIC
ALPHA(I,J)=XETA(I,J)*XETA(I,J)+YETA(I,J)*YETA(I,J)
GAMA(I,J)=XXI(I,J)*XXI(I,J)+YXI(I,J)*YXI(I,J)
BETA(I,J)=XXI(I,J)*XETA(I,J)+YXI(I,J)*YETA(I,J)
JCB(I,J)=XXI(I,J)*YETA(I,J)-XETA(I,J)*YXI(I,J)
END DO
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!结果的输出!!!!!!!!!!!!!!!!!!!!!!
OPEN(11,FILE='GRIDJCB.DAT',STATUS='REPLACE')
DO J=1,NJC
DO I=1,NIC
WRITE(11,*)JCB(I,J),ALPHA(I,J),BETA(I,J),GAMA(I,J),XXI(I,J),&
XETA(I,J),YXI(I,J),YETA(I,J)
END DO
END DO
CLOSE(11)
!!!!!!!!!!!!!!!!!!!!!!!!结果检查!!!!!!!!!!!!!!!!!!!!!!!!!
OPEN(12,FILE='CheckJCB.DAT',STATUS='REPLACE')
DO J=1,NJC
DO I=1,NIC
IF(JCB(I,J).LT.0)WRITE(12,*)"JCB",I,J,JCB(I,J)
IF(ALPHA(I,J).LT.0)WRITE(12,*)"ALPHA",I,J,ALPHA(I,J)
IF(BETA(I,J).LT.0)WRITE(12,*)"BETA",I,J,BETA(I,J)
IF(GAMA(I,J).LT.0)WRITE(12,*)"GAMA",I,J,GAMA(I,J)
IF(XXI(I,J).LT.0)WRITE(12,*)"XXI",I,J,XXI(I,J)
IF(XETA(I,J).LT.0)WRITE(12,*)"XETA",I,J,XETA(I,J)
IF(YXI(I,J).LT.0)WRITE(12,*)"YXI",I,J,YXI(I,J)
IF(YETA(I,J).LT.0)WRITE(12,*)"YETA",I,J,YETA(I,J)
END DO
END DO
CLOSE(12)
END PROGRAM