二维贴体坐标和常用坐标的转换问题



  • 课题组有一个较老的计算程序,使用的是贴体坐标系。控制方程也需要进行转换在贴体坐标系下相应的形式。如下微信图片_20190311211436.jpg
    里面的alpha和beta等参数定义如下
    微信图片_20190311211421.jpg
    XI和ETA是贴体坐标系的两个方向,我在求以上jacobi,beta,gama,alpha时遇到一些问题,请给位老看掌掌眼。老程序应用的网格中beta等参数没有负的,我的新网格会出现负的。但是从公式来看,确实存在负数的可能性。不得其解。我的网格是
    untitled.png
    参数求解的代码如下:

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !     计算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
    

  • CORE Fluent讲师

    Fortran 代码,令人感动



  • Fortran 代码,令人感动


Log in to reply