sim2d.f

对工业生产过程结晶过程的一个仿真程序软件包

     &   +f(i+1,j)*dble(j-1)**2
     &   +f(i+1,j+1)*dble(j)**2)/4.0D0
     &   *delR1*delR2**3
         mu21(1)=mu21(1)
     &   +(f(i,j)*dble(i-1)**2*dble(j-1)
     &   +f(i,j+1)*dble(i-1)**2*dble(j)
     &   +f(i+1,j)*dble(i)**2*dble(j-1)
     &   +f(i+1,j+1)*dble(i)**2*dble(j))/4.0D0
     &   *delR1**3*delR2**2
	 mu12(1)=mu12(1)
     &   +(f(i,j)*dble(i-1)*dble(j-1)**2
     &   +f(i,j+1)*dble(i-1)*dble(j)**2
     &   +f(i+1,j)*dble(i)*dble(j-1)**2
     &   +f(i+1,j+1)*dble(i)*dble(j)**2)/4.0D0
     &   *delR1**2*delR2**3
     	 mu22(1)=mu22(1)
     &   +(f(i,j)*dble(i-1)**2*dble(j-1)**2
     &   +f(i,j+1)*dble(i-1)**2*dble(j)**2
     &   +f(i+1,j)*dble(i)**2*dble(j-1)**2
     &   +f(i+1,j+1)*dble(i)**2*dble(j)**2)/4.0D0
     &   *delR1**3*delR2**3
         mu30(1)=mu30(1)
     &   +(f(i,j)*dble(i-1)**3
     &   +f(i,j+1)*dble(i-1)**3
     &   +f(i+1,j)*dble(i)**3
     &   +f(i+1,j+1)*dble(i)**3)/4.0D0
     &   *delR1**4*delR2      
         mu31(1)=mu31(1)
     &   +(f(i,j)*dble(i-1)**3*dble(j-1)
     &   +f(i,j+1)*dble(i-1)**3*dble(j)
     &   +f(i+1,j)*dble(i)**3*dble(j-1)
     &   +f(i+1,j+1)*dble(i)**3*dble(j))/4.0D0
     &   *delR1**4*delR2**2
        ENDDO
       ENDDO
      volume(1)=mu21(1)-2.0D0/3.0D0*mu30(1)

*     Current nucleation rate
      Brate(1)=birth(Conc(1),Concs(1),
     &     volume(1))
*     Current growth rate
      G1rate(1)=growth1(Conc(1), Concs(1))
      G2rate(1)=growth2(Conc(1), Concs(1))

*     Mesh boundary
      mratei=mratei+G1rate(1)*delT
      mratej=mratej+G2rate(2)*delT
      meshi1=starti1+DNINT(mratei)
      meshi2=starti2+DNINT(mratei)
      meshi3=starti3+DNINT(mratei)
      meshj1=startj1+DNINT(mratej)
      meshj2=startj2+DNINT(mratej)
      meshj3=startj3+DNINT(mratej)

*     write initial f(i,j) to the file
      nfile=10
      WRITE(filename(4:5), '(I2)')nfile
      OPEN(UNIT=20, FILE=filename, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')     
      DO i=1, totR1pts, skipR1
        DO j=1, totR2pts, skipR2
          WRITE(20,30) i, j, time(1), f(i,j)
        ENDDO
      ENDDO
      CLOSE(UNIT=20)
      m=1      
  
*     Simulation  starts    
      DO k=2,totTpts
	print*, "k=", k
*       Current time
        time(k)=time(k-1)+delT
*       Current temperature
        T(k)=temp(time(k))
*       Current saturated Conc
        Concs(k)=Csat(T(k))
*       Current concentration
        conc(k)=conc(k-1)-dble(2)*alpha*delT*G1rate(k-1)
     &           *(mu11(k-1)-mu20(k-1))
     &           -alpha*delT*G2rate(k-1)*mu20(k-1)
*       Current super concentration
        superConc(k)=(Conc(k)-Concs(k))/Concs(k)        
*       Compute population density f   
        vei=delT*G1rate(k-1)/delR1
        vej=delT*G2rate(k-1)/delR2  
*       Compute f in j-direction              
        DO i=1, meshi1
	  DO j=1, meshj1              
             IF((f(i,j+1)-f(i,j)).EQ.0.0D0)THEN
                j2theta=1.0D10*(f(i,j)-f(i,j-1))
             ELSE
                j2theta=(f(i,j)-f(i,j-1))/(f(i,j+1)-f(i,j)) 
             ENDIF
             IF(((f(i,j)-f(i,j-1)).EQ.0.0D0))then
                j1theta=1.0D10*(f(i,j-1)-f(i,j-2))
             ELSE
                j1theta=(f(i,j-1)-f(i,j-2))/(f(i,j)-f(i,j-1))
             ENDIF
             j1phi=(DABS(j1theta) +j1theta)/(1.0D0+DABS(j1theta))
             j2phi=(DABS(j2theta) +j2theta)/(1.0D0+DABS(j2theta))                
             IF (method.EQ.1) THEN
*               Upwind 
                ftemp(i,j)=f(i,j)-vej*(f(i,j)-f(i,j-1)) 
             ELSE IF (method.EQ.2) then
*               Lax-Wendroff
                ftemp(i,j)=f(i,j)-0.5D0*vej*(f(i,j+1)-f(i,j-1))+
     &               vej*vej*0.5D0*(f(i,j+1)-2.0D0*f(i,j)+f(i,j-1))
             ELSE
*               High resolution             
                ftemp(i,j)=f(i,j)-(vej-0.5D0*vej*(1.0D0-vej)*j1phi)
     &               *(f(i,j)-f(i,j-1))-0.5D0*vej*(1.0D0-vej)*j2phi
     &               *(f(i,j+1)-f(i,j)) 
             ENDIF
	  ENDDO
        ENDDO
        DO i=meshi2, meshi3-1
	  DO j=meshj2, meshj3-1              
             IF((f(i,j+1)-f(i,j)).EQ.0.0D0)THEN
                j2theta=1.0D10*(f(i,j)-f(i,j-1))
             ELSE
                j2theta=(f(i,j)-f(i,j-1))/(f(i,j+1)-f(i,j)) 
             ENDIF
             IF(((f(i,j)-f(i,j-1)).EQ.0.0D0))THEN
                j1theta=1.0D10*(f(i,j-1)-f(i,j-2))
             ELSE
                j1theta=(f(i,j-1)-f(i,j-2))/(f(i,j)-f(i,j-1))
             ENDIF
             j1phi=(DABS(j1theta) +j1theta)/(1.0D0+DABS(j1theta))
             j2phi=(DABS(j2theta) +j2theta)/(1.0D0+DABS(j2theta))
                
             IF (method.EQ.1) THEN
*               Upwind 
                ftemp(i,j)=f(i,j)-vej*(f(i,j)-f(i,j-1)) 
             ELSE IF (method.EQ.2) THEN
*               Lax-Wendroff
                ftemp(i,j)=f(i,j)-0.5D0*vej*(f(i,j+1)-f(i,j-1))+
     &               vej*vej*0.5D0*(f(i,j+1)-2.0D0*f(i,j)+f(i,j-1))
             ELSE
*              High resolution             
                ftemp(i,j)=f(i,j)-(vej-0.5D0*vej*(1.0D0-vej)*j1phi)
     &               *(f(i,j)-f(i,j-1))-0.5D0*vej*(1.0D0-vej)*j2phi
     &               *(f(i,j+1)-f(i,j)) 
            ENDIF
	  ENDDO
        ENDDO
*     Compute f in i-direction
        DO j=1, meshj1
	  DO i=1, meshi1
           IF((ftemp(i+1,j)-ftemp(i,j)).EQ.0.0D0)THEN
              i2theta=1.0D10*(ftemp(i,j)-ftemp(i-1,j))
           ELSE
              i2theta=(ftemp(i,j)-ftemp(i-1,j))/
     &                (ftemp(i+1,j)-ftemp(i,j))
           ENDIF
           IF(((ftemp(i,j)-ftemp(i-1,j)).EQ.0.0D0))THEN
              i1theta=1.0D10*(ftemp(i-1,j)-ftemp(i-2,j))
           ELSE
              i1theta=(ftemp(i-1,j)-ftemp(i-2,j))/
     &                (ftemp(i,j)-ftemp(i-1,j))
           ENDIF
           i1phi=(DABS(i1theta) +i1theta)/(1.0D0+DABS(i1theta))
           i2phi=(DABS(i2theta) +i2theta)/(1.0D0+DABS(i2theta))
           IF(method.EQ.1) THEN
*             Upwind
              f(i,j)=ftemp(i,j)-vei*(ftemp(i,j)-ftemp(i-1,j))
           ELSE IF(method.EQ.2) THEN 
*             Lax-Wendroff
              f(i,j)=ftemp(i,j)-0.5D0*vei*(ftemp(i+1,j)-ftemp(i-1,j))+
     &       vei*vei*0.5D0*(ftemp(i+1,j)-2.0D0*ftemp(i,j)+ftemp(i-1,j))
            ELSE
*             High resolution
               f(i,j)=ftemp(i,j)-(vei-0.5D0*vei*(1.0D0-vei)*i1phi)
     &         *(ftemp(i,j)-ftemp(i-1,j))-0.5D0*vei*(1.0D0-vei)*i2phi
     &         *(ftemp(i+1,j)-ftemp(i,j))
            ENDIF
	  ENDDO
        ENDDO
        DO j=meshj2, meshj3-1
	  DO i=meshi2, meshi3-1
           IF((ftemp(i+1,j)-ftemp(i,j)).EQ.0.0D0)THEN
              i2theta=1.0D10*(ftemp(i,j)-ftemp(i-1,j))
           ELSE
              i2theta=(ftemp(i,j)-ftemp(i-1,j))/
     &                (ftemp(i+1,j)-ftemp(i,j))
           ENDIF
           IF(((ftemp(i,j)-ftemp(i-1,j)).EQ.0.0D0))THEN
              i1theta=1.0D10*(ftemp(i-1,j)-ftemp(i-2,j))
           ELSE
              i1theta=(ftemp(i-1,j)-ftemp(i-2,j))/
     &                (ftemp(i,j)-ftemp(i-1,j))
           ENDIF
           i1phi=(DABS(i1theta) +i1theta)/(1.0D0+DABS(i1theta))
           i2phi=(DABS(i2theta) +i2theta)/(1.0D0+DABS(i2theta))
           IF(method.EQ.1) THEN
*             Upwind
              f(i,j)=ftemp(i,j)-vei*(ftemp(i,j)-ftemp(i-1,j))
           ELSE IF(method.EQ.2) then 
*             Lax-Wendroff
              f(i,j)=ftemp(i,j)-0.5D0*vei*(ftemp(i+1,j)-ftemp(i-1,j))+
     &        vei*vei*0.5D0*(ftemp(i+1,j)-2.0D0*ftemp(i,j)+ftemp(i-1,j))
           ELSE
*             High resolution
              f(i,j)=ftemp(i,j)-(vei-0.5D0*vei*(1.0D0-vei)*i1phi)
     &        *(ftemp(i,j)-ftemp(i-1,j))-0.5D0*vei*(1.0D0-vei)*i2phi
     &        *(ftemp(i+1,j)-ftemp(i,j))
           ENDIF
	  ENDDO
        ENDDO
*       Add nucleation 
        f(1,1)=f(1,1)+Brate(k-1)*delT/delR1/delR2

*	Current moments
       DO i=1,meshi1
        DO j=1,meshj1  
         mu00(k)=mu00(k)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))/4.0D0
     &           *delR1*delR2
	 mu10(k)=mu10(k)
     &    +(f(i,j)*dble(i-1)
     &    +f(i,j+1)*dble(i-1)
     &    +f(i+1,j)*dble(i)
     &    +f(i+1,j+1)*dble(i))/4.0D0
     &    *delR1**2*delR2
	 mu01(k)=mu01(k)
     &   +(f(i,j)*dble(j-1)
     &   +f(i,j+1)*dble(j)
     &   +f(i+1,j)*dble(j-1)
     &   +f(i+1,j+1)*dble(j))/4.0D0
     &   *delR1*delR2**2
	 mu11(k)=mu11(k)
     &   +(f(i,j)*dble(i-1)*dble(j-1)
     &   +f(i,j+1)*dble(i-1)*dble(j)
     &   +f(i+1,j)*dble(i)*dble(j-1)
     &   +f(i+1,j+1)*dble(i)*dble(j))/4.0D0
     &   *delR1**2*delR2**2
         mu20(k)=mu20(k)
     &   +(f(i,j)*dble(i-1)**2
     &   +f(i,j+1)*dble(i-1)**2
     &   +f(i+1,j)*dble(i)**2
     &   +f(i+1,j+1)*dble(i)**2)/4.0D0
     &   *delR1**3*delR2
   	 mu02(k)=mu02(k)
     &   +(f(i,j)*dble(j-1)**2
     &   +f(i,j+1)*dble(j)**2
     &   +f(i+1,j)*dble(j-1)**2
     &   +f(i+1,j+1)*dble(j)**2)/4.0D0
     &   *delR1*delR2**3
         mu21(k)=mu21(k)
     &   +(f(i,j)*dble(i-1)**2*dble(j-1)
     &   +f(i,j+1)*dble(i-1)**2*dble(j)
     &   +f(i+1,j)*dble(i)**2*dble(j-1)
     &   +f(i+1,j+1)*dble(i)**2*dble(j))/4.0D0
     &   *delR1**3*delR2**2
	 mu12(k)=mu12(k)
     &   +(f(i,j)*dble(i-1)*dble(j-1)**2
     &   +f(i,j+1)*dble(i-1)*dble(j)**2
     &   +f(i+1,j)*dble(i)*dble(j-1)**2
     &   +f(i+1,j+1)*dble(i)*dble(j)**2)/4.0D0
     &   *delR1**2*delR2**3
     	 mu22(k)=mu22(k)
     &   +(f(i,j)*dble(i-1)**2*dble(j-1)**2
     &   +f(i,j+1)*dble(i-1)**2*dble(j)**2
     &   +f(i+1,j)*dble(i)**2*dble(j-1)**2
     &   +f(i+1,j+1)*dble(i)**2*dble(j)**2)/4.0D0
     &   *delR1**3*delR2**3
         mu30(k)=mu30(k)
     &   +(f(i,j)*dble(i-1)**3
     &   +f(i,j+1)*dble(i-1)**3
     &   +f(i+1,j)*dble(i)**3
     &   +f(i+1,j+1)*dble(i)**3)/4.0D0
     &   *delR1**4*delR2      
         mu31(k)=mu31(k)
     &   +(f(i,j)*dble(i-1)**3*dble(j-1)
     &   +f(i,j+1)*dble(i-1)**3*dble(j)
     &   +f(i+1,j)*dble(i)**3*dble(j-1)
     &   +f(i+1,j+1)*dble(i)**3*dble(j))/4.0D0
     &   *delR1**4*delR2**2
        ENDDO
       ENDDO
       DO i=meshi2,meshi3-1
        DO j=meshi2,meshj3-1  
         mu00(k)=mu00(k)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))/4.0D0
     &           *delR1*delR2
	 mu10(k)=mu10(k)
     &    +(f(i,j)*dble(i-1)
     &    +f(i,j+1)*dble(i-1)
     &    +f(i+1,j)*dble(i)
     &    +f(i+1,j+1)*dble(i))/4.0D0
     &    *delR1**2*delR2
	 mu01(k)=mu01(k)
     &   +(f(i,j)*dble(j-1)
     &   +f(i,j+1)*dble(j)
     &   +f(i+1,j)*dble(j-1)
     &   +f(i+1,j+1)*dble(j))/4.0D0
     &   *delR1*delR2**2
	 mu11(k)=mu11(k)
     &   +(f(i,j)*dble(i-1)*dble(j-1)
     &   +f(i,j+1)*dble(i-1)*dble(j)
     &   +f(i+1,j)*dble(i)*dble(j-1)
     &   +f(i+1,j+1)*dble(i)*dble(j))/4.0D0
     &   *delR1**2*delR2**2
         mu20(k)=mu20(k)
     &   +(f(i,j)*dble(i-1)**2
     &   +f(i,j+1)*dble(i-1)**2
     &   +f(i+1,j)*dble(i)**2
     &   +f(i+1,j+1)*dble(i)**2)/4.0D0
     &   *delR1**3*delR2
   	 mu02(k)=mu02(k)
     &   +(f(i,j)*dble(j-1)**2
     &   +f(i,j+1)*dble(j)**2
     &   +f(i+1,j)*dble(j-1)**2
     &   +f(i+1,j+1)*dble(j)**2)/4.0D0
     &   *delR1*delR2**3
         mu21(k)=mu21(k)
     &   +(f(i,j)*dble(i-1)**2*dble(j-1)
     &   +f(i,j+1)*dble(i-1)**2*dble(j)
     &   +f(i+1,j)*dble(i)**2*dble(j-1)
     &   +f(i+1,j+1)*dble(i)**2*dble(j))/4.0D0
     &   *delR1**3*delR2**2
	 mu12(k)=mu12(k)
     &   +(f(i,j)*dble(i-1)*dble(j-1)**2
     &   +f(i,j+1)*dble(i-1)*dble(j)**2
     &   +f(i+1,j)*dble(i)*dble(j-1)**2
     &   +f(i+1,j+1)*dble(i)*dble(j)**2)/4.0D0