init0.for

给予DSMC的二维CFD模拟FORTRAN代码

*   INIT0.FOR
*
      SUBROUTINE INIT0
*
        include 'dsmc0.h'
*
*--set constants
*
      PI=3.141592654
      SPI=SQRT(PI)
      BOLTZ=1.3806D-23
*
      CALL DATA0
      NPR=0
*
*--set information on the cross-species collisions
*
      IF (MNSP.EQ.1) ISPD=0
      DO 100 N=1,MNSP
        DO 50 M=1,MNSP
          IF ((ISPD.EQ.0).OR.(N.EQ.M)) THEN
            SPM(1,N,M)=0.25*PI*(SP(1,N)+SP(1,M))**2
*--the collision cross section is assumed to be given by eqn (1.35)
            SPM(2,N,M)=0.5*(SP(2,N)+SP(2,M))
            SPM(3,N,M)=0.5*(SP(3,N)+SP(3,M))
            SPM(4,N,M)=0.5*(SP(4,N)+SP(4,M))
*--mean values are used for ISPD=0
          ELSE
            SPM(1,N,M)=PI*SPM(1,N,M)**2
*--the cross-collision diameter is converted to the cross-section
          END IF
          SPM(5,N,M)=(SP(5,N)/(SP(5,N)+SP(5,M)))*SP(5,M)
*--the reduced mass is defined in eqn (2.7)
          SPM(6,N,M)=GAM(2.5-SPM(3,N,M))
50      CONTINUE
100   CONTINUE
*
*--initialise variables
*
      TIME=0.
      NM=0
*
      CG(1,1)=XF
      CW=(XR-XF)/MNC
      DO 200 M=1,MNC
        IF (M.GT.1) CG(1,M)=CG(2,M-1)
        CG(2,M)=CG(1,M)+CW
        CG(3,M)=CW
        CC(M)=CW
        DO 150 L=1,MNSG
          DO 120 K=1,MNSG
            CCG(2,M,L,K)=RF(0)
            CCG(1,M,L,K)=SPM(1,1,1)*300.*SQRT(FTMP/300.)
120       CONTINUE
150     CONTINUE
*--the maximum value of the (rel. speed)*(cross-section) is set to a
*--reasonable, but low, initial value and will be increased as necessary
200   CONTINUE
*
*--set sub-cells
*
      DO 300 N=1,MNC
        DO 250 M=1,NSC
          L=(N-1)*NSC+M
          ISC(L)=N
250     CONTINUE
300   CONTINUE
*
*--generate initial gas in equilibrium at temperature FTMP
*
      DO 400 L=1,MNSP
        REM=0
        VMP=SQRT(2.*BOLTZ*FTMP/SP(5,L))
*--VMP is the most probable speed in species L, see eqns (4.1) and (4.7)
        DO 350 N=1,MNC
          A=FND*CG(3,N)*FSP(L)/FNUM+REM
*--A is the number of simulated molecules of species L in cell N to
*--simulate the required concentrations at a total number density of FND
            VYLO = (VFY(2)-VFY(1))*(CG(1,N)+CG(2,N))/2.0/(XR-XF)
C--  VYLO is the local velocity at t=0
          IF (N.LT.MNC) THEN
            MM=A
            REM=(A-MM)
*--the remainder REM is carried forward to the next cell
          ELSE
            MM=NINT(A)
          END IF
          DO 320 M=1,MM
            IF (NM.LE.MNM) THEN
*--round-off error could have taken NM to MNM+1
              NM=NM+1
              IPS(NM)=L
              PP(1,NM)=CG(1,N)+RF(0)*(CG(2,N)-CG(1,N))
	        IF (PP(1,NM).LT.(XF+SIGEMA))PP(1,NM)=XF+SIGEMA
	        IF (PP(1,NM).GT.(XR-SIGEMA))PP(1,NM)=XR-SIGEMA
			PP(2,NM)=0.0D0
	        PP(3,NM)=0.0D0
              IPL(NM)=(PP(1,NM)-CG(1,N))*(NSC-.001)/CG(3,N)+1+NSC*(N-1)
	IF(IPL(NM).GT.MNSC) IPL(NM)=MNSC
	IF(IPL(NM).LT.1)    IPL(NM)=1
*--species, position, and sub-cell number have been set
              DO 305 K=1,3
                CALL RVELC(PV(K,NM),A,VMP)
305           CONTINUE
                  PV(2,NM) = PV(2,NM) + VYLO
*--velocity components have been set
            END IF
320       CONTINUE
350     CONTINUE
400   CONTINUE
      WRITE (*,99001) NM
99001 FORMAT (' ',I6,' MOLECULES')
*
      RETURN
      END