compart.f

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

     &                *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  
           volume(k)=mu21(k)-2.0D0/3.0D0*mu30(k)         
        END IF
*       Current growth rate
        G1rate(k)=growth1(conc(k), Concs(k))
        G2rate(k)=growth2(conc(k), Concs(k))
*       Current nucleation rate
        Brate(k)=birth(conc(k),Concs(k), volume(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)
        IF(meshi3.GT.maxR1R2) THEN
           meshi3=maxR1R2
           WRITE(200+myrank,100)'meshi3 is larger than maxR1R2'
        END IF        
        IF(meshj3.GT.maxR1R2)THEN
           meshj3=maxR1R2
           WRITE(200+myrank,100)'meshj3 is larger than maxR1R2'
        END IF
*       Write f to file
        IF (K .EQ. skipNext) THEN 
           nfile=nfile+1
           nUnit=nfile+myrank
           WRITE(fname(2:3), '(I2)')nfile
           WRITE(fname(5:6), '(I2)')myrank+10
           OPEN(UNIT=nUnit, FILE=fname, 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(nUnit,30) i,j,time(k), f(i,j)
                 ELSE
                    WRITE(nUnit,30) i,j,time(k), 0.0D0
                 END IF
              ENDDO
           ENDDO
          CLOSE(UNIT=nUnit)
          m=m+1
          skipNext=m*skipT
       END IF
      ENDDO
      IF(numprocs.GT.1)THEN
         WRITE(mu1name(5:6), '(I2)')myrank+10
         WRITE(mu2name(5:6), '(I2)')myrank+10
         WRITE(mu3name(5:6), '(I2)')myrank+10
         WRITE(timename(6:7), '(I2)')myrank+10
         WRITE(concname(6:7), '(I2)')myrank+10
         WRITE(tgbname(5:6), '(I2)')myrank+10
      END IF
      OPEN(UNIT=myrank*10+21, FILE=mu1name, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=myrank*10+22, FILE=mu2name, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=myrank*10+23, FILE=mu3name, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=myrank*10+24, FILE=timename, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=myrank*10+25, FILE=concname, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=myrank*10+26, FILE=tgbname, FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      DO k=1,totTpts,skipT1 
         WRITE(myrank*10+24,10) time(k)
         WRITE(myrank*10+25,20) Concs(k), Conc(k), superConc(k)
         WRITE(myrank*10+26,25) time(k), T(k), G1rate(k), 
     &        G2rate(k), Brate(k)
         WRITE(myrank*10+21,15) mu00(k), mu10(k), mu01(k), mu11(k)
         WRITE(myrank*10+22,15) mu20(k), mu02(k), mu21(k), mu12(k)
         WRITE(myrank*10+23,15) mu22(k), mu30(k), mu31(k),volume(k)
      ENDDO
      IF(myrank.EQ.numprocs-1) THEN
         OPEN(UNIT=27, FILE='const.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN') 
         OPEN(UNIT=1000, FILE='domain.dat', FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='UNKNOWN') 
         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
         WRITE(27, 11) numprocs
         CLOSE (UNIT=27)
         CLOSE(UNIT=1000)
      END IF
      
1000  CLOSE (UNIT=myrank*10+21)
      CLOSE (UNIT=myrank*10+22)
      CLOSE (UNIT=myrank*10+23)
      CLOSE (UNIT=myrank*10+24)
      CLOSE (UNIT=myrank*10+25)
      CLOSE (UNIT=myrank*10+26)
      CLOSE (unit=200+myrank)
 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))
 21   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(6(I4,1x))
 70   FORMAT(E13.6,1x,6(I5,1x))
 80   FORMAT(E13.6,1x,4(I5,1x),E13.6)
 90   FORMAT(A,I4,A,I4,A,I4,A,I4,A,I4,A,I4)
 100  FORMAT(A)
 110  FORMAT(A,1x,I6)
 120  FORMAT(A,1X,E13.6)
 130  FORMAT(A,1X,I6,1x,I6)
 140  FORMAT(A,1X,I6,1x,E13.6)
 150  FORMAT(A,1X,E13.6,1X,I6,1X,I6)
 160  FORMAT(A,1X,E13.6,1X,E13.6)
 170  FORMAT(E23.16)
 180  FORMAT(2(I9,1X))
      CALL MPI_FINALIZE(ierr)
      stop 
      end
*************************************************************************
      REAL*8 FUNCTION Temp(time)

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


      REAL*8 time,tot_time,init_temp,final_temp
      COMMON /simTIME/tot_time
      COMMON /temperature/init_temp,final_temp


      Temp=(final_temp-init_temp)/tot_time*time+init_temp


      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
*     from rudy g/g-solvent
      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,rhoc,Temp,g1,kg1 
      COMMON /g1data/g1,kg1
      COMMON /rhoc/rhoc

      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,Temp, rhoc,g2,kg2
      COMMON /g2data/g2, kg2 
      COMMON /rhoc/rhoc

      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
*                  volume - total volume
*     non-argument input: kb, b kinetic rate parameters
*     output: birth - birth rate

      REAL*8 conc, concs, volume, kb, b
      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)*
     &         .3987662063472828D+01
      ELSE
          Fzero=0.0D0  
      END IF
      
      RETURN
      END