RohanCFD's Blog

Happy Life, Happy Cfd

二维热传导计算程序

只适用于二维情况下的矩形,超级简单的程序,带的参数也少,是自己写的第一个程序~哈哈

初温为0,边界条件自己设定

 

        PROGRAM HEATTRANSFER_2D
        IMPLICIT NONE

        CHARACTER :: FILENAME*30,STEP*6
        INTEGER :: XMAX,YMAX,ZMAX,N,NMAX,STEPD
        INTEGER :: I,J,K,S(4)
        REAL :: DX,DY,DT,DT0,Kr,T0(1:4),A,B,T
        REAL,ALLOCATABLE :: X(:,:,:),Y(:,:,:),Z(:,:,:),TEMP(:,:,:)

        WRITE(*,*) "PLEASE INPUT FILENAME:"
        READ(*,*) FILENAME
        !----------------------------------------------------读取网格文件--------------------------------------------------------------------------
        OPEN(10,FILE=TRIM(FILENAME))       !Read Mesh File
        READ(10,*) XMAX,YMAX               !Get Mesh Size          
        ZMAX=1                             
        ALLOCATE(X(XMAX,YMAX,ZMAX),Y(XMAX,YMAX,ZMAX),Z(XMAX,YMAX,ZMAX))   !生成网格坐标存储空间
        ALLOCATE(TEMP(XMAX,YMAX,ZMAX))                                    !生成温度存储空间
        READ(10,*) (((X(I,J,K),I=1,XMAX),J=1,YMAX),K=1,ZMAX),(((Y(I,J,K),I=1,XMAX),J=1,YMAX),K=1,ZMAX),(((Z(I,J,K),I=1,XMAX),J=1,YMAX),K=1,ZMAX)
        CLOSE(10)
        !-----------------------------------------------------确定步长时间------------------------------------------------------------------------
        DX=(MAXVAL(X)-MINVAL(X))/XMAX
        DY=(MAXVAL(Y)-MINVAL(Y))/YMAX
        
        WRITE(*,*) 'Please input the rate of heat transfer k'
        READ(*,*) Kr
        DT0=(1/2.0/K)/(1/DX**2+1/DY**2)/2.0
		WRITE(*,*) 'Please input the step time, suggest it smaller than ',DT0
		READ(*,*) DT
		
		!------------------------------------------------------确定边界\初值条件-----------------------------------------------------------------------
		TEMP(1:XMAX,1:YMAX,1:ZMAX)=0
		!------------------确定边界条件类型------------------
		DO I=1,4
				WRITE(*,*) "Select Boundary Condition Type For"
				SELECT CASE(I)
				    CASE (1) 
				        WRITE(*,*) "Left side"
				    CASE (2)
				        WRITE(*,*) "Bottom"
				    CASE (3)
				        WRITE(*,*) "Right side"
				    CASE (4) 
				        WRITE(*,*) "Top"
				END SELECT
				WRITE(*,*) "1\Temperature"
				WRITE(*,*) "2\Rate of Temperatur Spread"
				READ(*,*) S(I)
		END DO
		!-----------------确定边界条件数值-------------------
		WRITE(*,*) "Please input the number for boundaries of left, bottom, right and up,use comma to seperate./The number should be wrote in 7 chracter and 3 float number(including comma)"
		READ(*,"(F7.2,','F7.3,',',F7.3,','F7.3)") T0(1:4)
		!----------------------------------------------------------------输入迭代次数\显示代数------------------------------------------------------------
		WRITE(*,*) 'Please input the number of steps'
		READ(*,*) NMAX
		WRITE(*,*) "Please input the number between two steps to show"
		READ(*,*) STEPD
		N=0
		T=0
100     N=N+1
        T=T+DT

		SELECT CASE(S(1))
		    CASE(1) 
		        TEMP(1,1:YMAX,1:ZMAX)=T0(1)
		    CASE(2) 
		        TEMP(1,1:YMAX,1:ZMAX)=(0-T0(1)*2*DX-TEMP(3,1:YMAX,1:ZMAX)+4*TEMP(2,1:YMAX,1:ZMAX))/3
        END SELECT
        SELECT CASE(S(2))
		    CASE(1) 
		        TEMP(1:XMAX,1,1:ZMAX)=T0(2)
		    CASE(2) 
		        TEMP(1:XMAX,1,1:ZMAX)=(0-T0(2)*2*DY-TEMP(1:XMAX,3,1:ZMAX)+4*TEMP(1:XMAX,2,1:ZMAX))/3
        END SELECT
        SELECT CASE(S(3))
		    CASE(1) 
		        TEMP(XMAX,1:YMAX,1:ZMAX)=T0(3)
		    CASE(2) 
		        TEMP(XMAX,1:YMAX,1:ZMAX)=(0-T0(3)*2*DX-TEMP(XMAX-2,1:YMAX,1:ZMAX)+4*TEMP(XMAX-1,1:YMAX,1:ZMAX))/3
        END SELECT
        SELECT CASE(S(4))
		    CASE(1) 
		        TEMP(1:XMAX,YMAX,1:ZMAX)=T0(4)
		    CASE(2) 
		        TEMP(1:XMAX,YMAX,1:ZMAX)=(0-T0(1)*2*DY-TEMP(1:XMAX,YMAX-2,1:ZMAX)+4*TEMP(1:XMAX,YMAX-3,1:ZMAX))/3
        END SELECT

        !-----------------------------------------------------二阶中心差分格式显式推进计算---------------------------------------------------------------------------
       	A=Kr*DT/DX**2
		B=Kr*DT/DY**2
		K=1
		XDO:DO I=2,XMAX-1
			YDO:DO J=2,YMAX-1
				TEMP(I,J,K)=A*(TEMP(I+1,J,K)+TEMP(I-1,J,K))+B*(TEMP(I,J+1,K)+TEMP(I,J-1,K))+(1-2*(A+B))*TEMP(I,J,K)
			END DO YDO
		END DO XDO
        !------------------------------------------------------收敛之判断-----------------------------------------------------------------------------------
        IF (MOD(N,STEPD) == 0) THEN
			WRITE(STEP,"(I6)") N
			OPEN (1,FILE='OUTPUT'//STEP//'_'//TRIM(FILENAME))
			WRITE(1,*) 'VARIABLES=','"X","Y","T"'
			WRITE(1,*) 'ZONE,I=',XMAX,',J=',YMAX,',F=POINT'
			DO I=1,XMAX
				DO J=1,YMAX
					WRITE(1,*) X(I,J,K),Y(I,J,K),TEMP(I,J,K)
				END DO
			END DO
			CLOSE(1)
		ENDIF
		
        WRITE(6,*) 'TIME=',T,',TEMPERATURE=',TEMP(XMAX/2,YMAX/2,K)
		IF(N < NMAX)  GOTO 100
		IF(N == NMAX) GOTO 101

101		WRITE(*,*) "THE COMPUTATION IS OVER!"
        
        WRITE(STEP,'(I6)') NMAX
        OPEN (2,FILE='OUTPUT'//STEP//'_'//TRIM(FILENAME))
        WRITE(2,*) 'VARIABLES=','"X","Y","T"'
        WRITE(2,*) 'ZONE,I=',XMAX,',J=',YMAX,',F=POINT'
        DO I=1,XMAX
            DO J=1,YMAX
            	WRITE(2,*) X(I,J,K),Y(I,J,K),TEMP(I,J,K)
            END DO
        END DO
        CLOSE(2)
        

        STOP
        END PROGRAM

PLot3D with Fortran

 

PLot3D是一种在CFD领域很常见的格式。几乎所有的网格生成软件都支持这种格式的输出,因此搞清楚这种文件的格式的读入和输出对CFD研究人员的工作带来很大的帮助。这里给出的都是读入ASCII格式文件的代码,输出方式类似。

PLot3D的网格文件后缀名为"xyz",可能是二进制(bin)格式或者文本文件格式(ASCII)。

1.单块网格(Single Grid

 读入3维单块网格

PROGRAM PLOT3D_SINGLE 
IMPLICIT NONE
CHARACTER :: BUFFER*30
INTEGER :: I,J,K,IMAX,JMAX,KMAX
REAL,ALLOCATABLE::MESH(:,:,:,:) !MESH存储网格的点坐标

WRITE(*,*)"Please input file name:"
READ(*,*)BUFFER
OPEN(200,FILE=TRIM(buffer))
READ(200,*) IMAX,JMAX,KMAX !获得网格尺寸
ALLOCATE(MESH(3,IMAX,JMAX,KMAX)) !划分网格所需要的内存空间
READ(200,*)      (((MESH(1,I,J,K), I=1,IMAX), J=1,JMAX), K=1,KMAX),&
                       (((MESH(2,I,J,K), I=1,IMAX), J=1,JMAX), K=1,KMAX),&
                       (((MESH(3,I,J,K), I=1,IMAX), J=1,JMAX), K=1,KMAX)
CLOSE(200)
STOP
END PROGRAM

2.多块网格(Multiple Grids

读入3维多块网格(似乎GRIDGEN读入网格时需采用这种格式的文件,你可以将其改成输出程序,所输出的网格就可以用GRIDGEN读入了)

PROGRAM PLOT3D_MULTI
IMPLICIT NONE
CHARACTER :: BUFFER*30
INTEGER :: I,J,K,IMAX,JMAX,KMAX,NMESH,N
INTEGER,ALLOCATABLE::IJK(:,:)   !IJK存储每块网格的尺寸
REAL,ALLOCATABLE::MESH(:,:,:,:,:) !MESH存储网格的点坐标

WRITE(*,*)"Please input file name:"
READ(*,*)BUFFER
OPEN(200,FILE=TRIM(buffer))
READ(200,*)NMESH !读入网格块数
ALLOCATE(IJK(3,NMESH))
READ(200,*) ((IJK(1,N),IJK(2,N),IJK(3,N)), N=1,NMESH) !获得多块网格尺寸
DO N=1,NMESH
  IMAX=MAX(IMAX,IJK(1,N))
  JMAX=MAX(IMAX,IJK(2,N))
  KMAX=MAX(IMAX,IJK(3,N))
END DO
ALLOCATE(MESH(3,IMAX,JMAX,KMAX,NMESH)) !划分网格所需要的内存空间
READ(200,*)    ((((MESH(1,I,J,K,N), I=1,IJK(1,N)), J=1,IJK(2,N)), K=1,IJK(3,N)),&
                      (((MESH(2,I,J,K,N), I=1,IJK(1,N)), J=1,IJK(2,N)), K=1,IJK(3,N)),&
                      (((MESH(3,I,J,K,N), I=1,IJK(1,N)), J=1,IJK(2,N)), K=1,IJK(3,N)),N=1,NMESH)
CLOSE(200)
STOP
END PROGRAM