sim2d.f

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

     &   *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
       volume(k)=mu21(k)-2.0D0/3.0D0*mu30(k)

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

*      Mesh boundary
       mratei=mratei+G1rate(k)*delT
       mratej=mratej+G2rate(k)*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 f(i,j) to the file
       IF ( K .EQ. skipNext) THEN   
          nfile=nfile+1
          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
                IF (f(i,j).GE.1.0D-16.OR.(f(i,j).LE.-1.0D-16)) THEN
                   WRITE(20,30) i,j,time(k), f(i,j)
                ELSE
                   WRITE(20,30) i,j,time(k), 0.0D0
                END IF
             ENDDO
          ENDDO
          CLOSE(UNIT=20)
          m=m+1
          skipNext=m*skipT
       END IF
      ENDDO
      
      IF(method.EQ.1)THEN
         OPEN(UNIT=21, FILE='upwmu1.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=22, FILE='upwmu2.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=23, FILE='upwmu3.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=24, FILE='time.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=25, FILE='upwconc.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=26, FILE='upwtgb.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=27, FILE='const.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')  
      ELSE IF(method.EQ.2) THEN
         OPEN(UNIT=21, FILE='laxmu1.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=22, FILE='laxmu2.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=23, FILE='laxmu3.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=24, FILE='time.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=25, FILE='laxconc.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=26, FILE='laxtgb.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=27, FILE='const.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN') 
      ELSE
         OPEN(UNIT=21, FILE='higmu1.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=22, FILE='higmu2.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=23, FILE='higmu3.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=24, FILE='time.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=25, FILE='higconc.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=26, FILE='higtgb.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
         OPEN(UNIT=27, FILE='const.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN') 
      END IF
      DO k=1,totTpts,skipT1 
        WRITE(24,10) time(k)
        WRITE(25,20) Concs(k), Conc(k), superConc(k)
        WRITE(26,25) time(k), T(k), G1rate(k), G2rate(k), Brate(k)
        WRITE(21,15) mu00(k), mu10(k), mu01(k), mu11(k)
        WRITE(22,15) mu20(k), mu02(k), mu21(k), mu12(k)
        WRITE(23,15) mu22(k), mu30(k), mu31(k), volume(k) 
      ENDDO
      WRITE(27, 60) maxR1R2
      WRITE(27, 10) totTime
      WRITE(27, 10) delT
      WRITE(27, 10) delR1
      WRITE(27, 10) delR2
      WRITE(27, 11) skipT
      WRITE(27, 11) skipR1
      WRITE(27, 11) skipR2
      WRITE(27, 11) totR1pts
      WRITE(27, 11) totR2pts
      WRITE(27, 11) totTpts
      CLOSE (UNIT=20)
      CLOSE (UNIT=21)
      CLOSE (UNIT=22)
      CLOSE (UNIT=23)
      CLOSE (UNIT=24)
      CLOSE (UNIT=25)
      CLOSE (UNIT=26)
      CLOSE (UNIT=27)      
 10   FORMAT(E13.6, 1X)
 11   FORMAT(I9)
 12   FORMAT(2(E13.6, 1X))
 15   FORMAT(4(E13.6,1X))
 20   FORMAT(3(E13.6,1X))
 25   FORMAT(5(E13.6,1X))  
 30   FORMAT(I8, I8, 2(E13.6,1X)) 
 40   FORMAT(E13.6, 1X)
 50   FORMAT(E13.6,1X, 3(I5,1x),E13.6)
 60   FORMAT(I9)

      PRINT*, "Done!"
      STOP
      END

*************************************************************************
      REAL*8 FUNCTION Temp(time)

*     Crystallizer temperature setpoint profile for the
*     simulation of a batch cooling crystallizer
*
*     input:   time - minutes
*     output:  Temp - temperature in  degrees Centigrade


      REAL*8 time

      Temp=(28.0D0-33.0D0)/7200.0*time+33.0D0

     
      RETURN
      END
      
      
**-*-*-*-*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-***

      REAL*8 FUNCTION Csat(Temp)

*     saturation concentration for the simulation of a
*     cooling batch crystallizer (potassium nitrate-water)
*     system, from Appendix C in Miller
*
*     input:  Temp - temperature (30-50 degree Centigrade)
*     output: Csat - saturation concentration 
*                    (g KDP/g water)

      REAL*8 Temp
      REAL*8 CA, CB
      COMMON /Csat_DATA/CA, CB

*      Csat=CA*Temp+CB
      Csat=9.3027D-5*Temp*Temp - 9.7629D-5*Temp + 0.2087D0

      RETURN
      END

**-*-*-*-*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-***

      REAL*8 FUNCTION growth1(conc, concs)

*     growth rate for the simulation of a cooling batch
*     crystallizer
*     
*     arguments:  conc - solute concentration
*                 concs - saturation concentration
*     non-argument input: D1, E1 kinetic rate parameters
*     output: growth - growth rate

      REAL*8 conc,concs,g1,kg1 
      COMMON /g1data/g1, kg1


      IF(conc.GE.concs)THEN
         growth1=kg1*((conc-concs)/concs)**g1
      ELSE
         growth1=0.0D0
      END IF

      RETURN
      END
      
**-*-*-*-*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-***
 
      REAL*8 FUNCTION growth2(conc, concs)

*     growth rate for the simulation of a cooling batch
*     crystallizer
*     
*     arguments:  conc - solute concentration
*                 concs - saturation concentration
*     non-argument input: D2, E2  kinetic rate parameters
*     output: growth - growth rate

      REAL*8 conc,concs,g2,kg2
      COMMON /g2data/g2, kg2


      IF(conc.GE.concs)THEN
         growth2=kg2*((conc-concs)/concs)**g2
      ELSE
         growth2=0.0D0
      END IF

      RETURN
      END

**-*-*-*-*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-*-*-*-*-*_*-*-***


      REAL*8 FUNCTION birth(conc, concs,volume)

*     birthth rate for the simulation of a cooling batch
*     crystallizer
*
*     arguments:  conc - solute concentration
*                 concs - saturation concentration
*                 m21 - 3rd moment
*     non-argument input: kb, b kinetic rate parameters
*     output: birth - birth rate

      REAL*8 conc,concs,kb,b,volume
      COMMON /bdata/kb, b


      IF(conc.GE.concs)THEN
         birth=kb*((conc-concs)/concs)**b*volume
      ELSE
         birth=0.0D0
      END IF

      RETURN
      END

*******************************************************************
      REAL*8 FUNCTION Fzero(r1, r2)
      
*     Initial population density f
*
*     arguments:  r1 - first length of crystal
*                 r2 - second length of crystal
*     output:     Fzero - inital population density for the given r1, r2
*                         # crystals /g solvent/micron**2 

*       population density function:
*
*       f0= -0.00034786r1^2 + 0.1363609*r1 - 13.2743
*           -0.00034786r2^2 + 0.1363609*r2 - 13.2743
*
*       which is in units of number of crystals/g solvent/micron.
*       This function was based on assuming 0.05 g as the mass in
*       1650 grams of solvent and 180.5-211.5 microns as the size range
*       for the seed crystals, and assuming a quadratic distibution
*       function.  The values are valid if mass of seed crystals
*       is scaled up proportionally to the mass of solvent.

      REAL*8 r1, r2
      LOGICAL inRange

      inRange=(r1 .GE.180.5) .AND. (r1 .LE. 211.5) .AND. 
     &        (r2 .GE.180.5) .AND. (r2 .LE. 211.5)
*
      IF (inRange) THEN
          Fzero= -0.00034786*r1**2 + 0.1363609*r1 - 13.2743
     &           -0.00034786*r2**2 + 0.1363609*r2 - 13.2743
      ELSE
          Fzero=0D0  
      END IF
      
      RETURN
      END