sim1d.f

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

* Copyright c 1998-2002 The Board of Trustees of the University of Illinois
*                 All rights reserved.
* Developed by: Large Scale Systems Research Laboratory
*               Professor Richard Braatz, Director
*               Department of Chemical Engineering
*               University of Illinois
*               http://brahms.scs.uiuc.edu
* 
* Permission hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to
* deal with the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the 
* Software is furnished to do so, subject to the following conditions:
*               1. Redistributions of source code must retain the above copyright
*                  notice, this list of conditions and the following disclaimers.
*               2. Redistributions in binary form must reproduce the above 
*                  copyright notice, this list of conditions and the following 
*                  disclaimers in the documentation and/or other materials 
*                  provided with the distribution.
*               3. Neither the names of Large Scale Research Systems Laboratory,
*                  University of Illinois, nor the names of its contributors may
*                  be used to endorse or promote products derived from this 
*                  Software without specific prior written permission.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
* THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
* OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
********************************************************************************
* sim21.f
*
*   This program is used to simulate batch crystallizer under perfect mixing assumption.
*   Crystals have one characteristic growth dimensions 


      PROGRAM Main
      IMPLICIT NONE      
      INTEGER i,j, maxpoint, totTpts, totR1pts, maxr1
      INTEGER skipT, skipR1, skipNext, k, nfile,maxf
      REAL*8 delT, delR1, totTime, tUnit
      PARAMETER ( maxpoint=4830, totTime=4800.0D0, maxr1=1600)
      PARAMETER (delT=1D0,delR1=1D0,skipR1=1,skipT=4801,maxf=3202)
*     f(0,1) and f(0,2) are imaginary boudary points, used for method 2. 
      REAL*8 f(0:maxr1, 2), Conc(maxpoint)
      REAL*8 moment1(maxpoint),moment2(maxpoint)
      REAL*8 moment0(maxpoint), moment3(maxpoint)
      REAL*8 time(maxpoint), Concs(maxpoint)
      REAL*8 superConc(maxpoint)
      REAL*8 g, kg, b, kb, R, kv, rhoc
      REAL*8 T(maxpoint), Grate(maxpoint), Brate(maxpoint)
      REAL*8 Temp, Csat, growth, birth, Fzero 
      CHARACTER*9 filename
      INTEGER mesh1, mesh2, mesh3
      REAL*8 mrate, theta1, theta2, phi1, phi2, ve
      
      DATA moment0/maxpoint*0.0D0/
      DATA moment1/maxpoint*0.0D0/
      DATA moment2/maxpoint*0.0D0/
      DATA moment3/maxpoint*0.0D0/
      DATA f/maxf*0.0D0/
         
      COMMON g, kg, b, kb
      COMMON /unit/ tUnit
      
      tUnit=60D0
      skipNext=skipT
      mesh1=DINT(10.0D0/delR1)
      mesh2=DINT(90.0D0/delR1)
      mesh3=DINT(320.0D0/delR1)
      mrate=0.0D0
      filename='popxx.dat'


*     Growth and nucleation kinetic parameters (Table 4.6 in Miller)
*     (dimensionaless)
      g=1.32D0
      print*, "g=", g
*     (mirons/sec, or min depending on tUnit)
      kg=DEXP(8.849D0)/tUnit
      print*, "kg=", kg
*     (dimensionless)
      b=1.78D0
      print*, "b=", b
*     (number of particles/micron^3/sec or min depending on tUnit) 
*     (the units have been corrected from that reported in
*     Table 3.1 in Miller)
      kb=DEXP(17.142D0)*(1.0D-04)**3/tUnit 
      print*, "kb=", kb
*     Initial size of crystal
      R=0.0D0       
*     volume shape factor (Appendix C in Miller)
      kv=1.0D0
*     density of crystal in g/cm^3 (Appendix C in Miller)
      rhoc=2.11D0  
*     Total time discretization pts
      totTpts=int(totTime/delT)+1
      totR1pts=int(dble(maxr1)/delR1)+1 

*     Initialization
*     Initial time
      time(1)=0.0D0

*     Concentration g/g solvent
      Conc(1)=0.493D0

*     Initial temperature
      T(1)=temp(time(1))

*     Initial saturated Conc
      Concs(1)=Csat(T(1))

*     Inital super concentration
      superConc(1)=(Conc(1)-Concs(1))/Concs(1)


*     Boundary condition at f(r1,0)
      DO i=1, totR1pts-1
         f(i,1)=Fzero(delR1*dble(i))
      ENDDO

*     Boundary condition at f(maxr1,t)
      f(totR1pts,1)=0.0D0

*     Boundary condition at f(-1,t)
      f(0,1)=0.0D0
      f(0,2)=0.0D0

*     Initial moments
      DO i=1,mesh1
         moment0(1)=moment0(1)+(f(i,1)+f(i+1,1))*0.5D0*delR1
	 moment1(1)=moment1(1)+(f(i,1)*dble(i-1)*delR1+f(i+1,1)*dble(i)
     &              *delR1)*0.5D0*delR1
	 moment2(1)=moment2(1)+(f(i,1)*dble(i-1)**2*delR1**2+f(i+1,1)*
     &              dble(i)**2*delR1**2)*0.5D0*delR1
     	 moment3(1)=moment3(1)+(f(i,1)*dble(i-1)**3*delR1**3+f(i+1,1)*
     &              dble(i)**3*delR1**3)*0.5D0*delR1
      ENDDO
      DO i=mesh2,mesh3
         moment0(1)=moment0(1)+(f(i,1)+f(i+1,1))*0.5D0*delR1
	 moment1(1)=moment1(1)+(f(i,1)*dble(i-1)*delR1+f(i+1,1)*dble(i)
     &              *delR1)*0.5D0*delR1
	 moment2(1)=moment2(1)+(f(i,1)*dble(i-1)**2*delR1**2+f(i+1,1)*
     &              dble(i)**2*delR1**2)*0.5D0*delR1
     	 moment3(1)=moment3(1)+(f(i,1)*dble(i-1)**3*delR1**3+f(i+1,1)*
     &              dble(i)**3*delR1**3)*0.5D0*delR1
      ENDDO


*     Initial nucleation rate
      Brate(1)=birth(Conc(1),Concs(1), moment3(1))

*     Current growth rate
      Grate(1)=growth(Conc(1), Concs(1))

*     Compute mesh boundary
      mrate=mrate+Grate(1)*delT
      mesh1=mesh1+DNINT(mrate)
      mesh2=mesh2+DNINT(mrate)
      mesh3=mesh3+DNINT(mrate)
      
*     Write f to 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
         WRITE(20,30) i, time(1), f(i,1)
      ENDDO
      CLOSE(UNIT=20)
      k=1
      
*     Simulation starts
      print*, "running"
      DO j=2,totTpts
*       Current time
        time(j)=time(j-1)+delT
*       Current temperature
        T(j)=temp(time(j))
*       Current Csat
        Concs(j)=Csat(T(j))
*       Current concentration
        Conc(j)=Conc(j-1)-dble(3)*delT*rhoc*kv*Grate(j-1)
     &          *moment2(j-1)*(1.0D-04)**3+delT*rhoc
     &          *kv*Brate(j-1)*R**3*(1.0D-04)**3
*       Current relative  super Csat
        superConc(j)=(Conc(j)-Concs(j))/Concs(j)
*       Boundary condition at f(0,t) and f(maxr1,t)
        f(totR1pts,2)=0.0D0
*       Current population density f           
        DO i=1, mesh1
           if((f(i+1,1)-f(i,1)).EQ.0.0D0)then
              theta2=0.0D0
           else
              theta2=(f(i,1)-f(i-1,1))/(f(i+1,1)-f(i,1))
           endif
           if(((f(i,1)-f(i-1,1)).EQ.0.0D0))then
              theta1=0.0D0
           else
              theta1=(f(i-1,1)-f(i-2,1))/(f(i,1)-f(i-1,1))
           endif
           phi1=(DABS(theta1) +theta1)/(1.0D0+DABS(theta1))
           phi2=(DABS(theta2) +theta2)/(1.0D0+DABS(theta2))
           ve=delT/delR1*Grate(j-1)
            f(i,2)=f(i,1)-ve*(f(i,1)-f(i-1,1))-ve/2.0D0*
     &      (1.0D0-ve)*((f(i+1,1)-f(i,1))*phi2-
     &                 (f(i,1)-f(i-1,1))*phi1)
        ENDDO
        f(1,2)=f(1,2)+Brate(j-1)*delT/delR1
*        f(1,2)=Brate(j-1)/Grate(j-1)

        DO i=mesh2, mesh3
           if((f(i+1,1)-f(i,1)).EQ.0.0D0)then
              theta2=0.0D0
           else
              theta2=(f(i,1)-f(i-1,1))/(f(i+1,1)-f(i,1))
           endif
           if(((f(i,1)-f(i-1,1)).EQ.0.0D0))then
              theta1=0.0D0
           else
              theta1=(f(i-1,1)-f(i-2,1))/(f(i,1)-f(i-1,1))
           endif
           phi1=(DABS(theta1) +theta1)/(1.0D0+DABS(theta1))
           phi2=(DABS(theta2) +theta2)/(1.0D0+DABS(theta2))
           ve=delT/delR1*Grate(j-1)
           f(i,2)=f(i,1)-ve*(f(i,1)-f(i-1,1))-ve/2.0D0*
     &          (1.0D0-ve)*((f(i+1,1)-f(i,1))*phi2-
     &          (f(i,1)-f(i-1,1))*phi1)
          
        ENDDO
******************************************************************************       
*	Current moments
        DO i=1,mesh1
         moment0(j)=moment0(j)+(f(i,2)+f(i+1,2))*0.5D0*delR1
	 moment1(j)=moment1(j)+(f(i,2)*dble(i-1)*delR1+f(i+1,2)*dble(i)
     &              *delR1)*0.5D0*delR1
	 moment2(j)=moment2(j)+(f(i,2)*dble(i-1)**2*delR1**2+f(i+1,2)*
     &              dble(i)**2*delR1**2)*0.5D0*delR1
     	 moment3(j)=moment3(j)+(f(i,2)*dble(i-1)**3*delR1**3+f(i+1,2)*
     &              dble(i)**3*delR1**3)*0.5D0*delR1
        ENDDO

        DO i=mesh2,mesh3
         moment0(j)=moment0(j)+(f(i,2)+f(i+1,2))*0.5D0*delR1
	 moment1(j)=moment1(j)+(f(i,2)*dble(i-1)*delR1+f(i+1,2)*dble(i)
     &              *delR1)*0.5D0*delR1
	 moment2(j)=moment2(j)+(f(i,2)*dble(i-1)**2*delR1**2+f(i+1,2)*
     &              dble(i)**2*delR1**2)*0.5D0*delR1
     	 moment3(j)=moment3(j)+(f(i,2)*dble(i-1)**3*delR1**3+f(i+1,2)*
     &              dble(i)**3*delR1**3)*0.5D0*delR1
        ENDDO        
*       Current growth rate
        Grate(j)=growth(Conc(j), Concs(j))
*        print*, "growth=", Grate(j)
*       Current nucleation rate
        Brate(j)=birth(Conc(j),Concs(j), moment3(j))
*        print*, "birth=", Brate(j)

*       Compute mesh boundary
        mrate=mrate+Grate(j)*delT
        mesh1=10+DNINT(mrate)
        mesh2=90+DNINT(mrate)
        mesh3=320+DNINT(mrate)       
*       Update f(i,1)
        DO i=1, totR1pts
           f(i,1)=f(i,2)
        ENDDO
*       Write f to file 
        IF ( j .EQ. skipNext) THEN            
          nfile=nfile+1
          WRITE(filename(4:5), '(I2)')nfile
*          WRITE(50,40)time(j), mesh1, mesh2, mesh3, mrate
          OPEN(UNIT=20, FILE=filename, FORM='FORMATTED',
     &    ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
          DO i=1, totR1pts,skipR1
             IF((f(i,1).GE.1.0D-16).OR.(f(i,1).LE.-1.0D-16)) then
                WRITE(20,30) i,time(j), f(i,1)
             ELSE 
                WRITE(20,30) i,time(j), 0.0D0
             ENDIF
          ENDDO
          CLOSE(UNIT=20)
          k=k+1
          skipNext=k*skipT
        ENDIF


      ENDDO

      OPEN(UNIT=21, FILE='moments.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=22, FILE='time.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=23, FILE='conc.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=24, FILE='tgb.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
      OPEN(UNIT=25, FILE='const.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')      
      DO j=1,totTpts,skipT 
        WRITE(22,10) time(j)
        WRITE(23,20) Concs(j), Conc(j), superConc(j)
        WRITE(24,15) time(j), T(j), Grate(j), Brate(j)
        WRITE(21,15) moment0(j), moment1(j), moment2(j), moment3(j)
      ENDDO
      WRITE(25, 11) maxr1
      WRITE(25, 10) totTime
      WRITE(25, 10) delT
      WRITE(25, 10) delR1
      WRITE(25, 11) skipT
      WRITE(25, 11) skipR1
      WRITE(25, 11) totR1pts
      WRITE(25, 11) totTpts

      CLOSE (UNIT=20)
      CLOSE (UNIT=21)
      CLOSE (UNIT=22)
      CLOSE (UNIT=23)
      CLOSE (UNIT=24)
      CLOSE (UNIT=25)
      CLOSE (UNIT=50)


      
 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)) 
 30   FORMAT(I8, 2(E13.6,1X)) 
 40   FORMAT(E13.6,1X, 3(I5,1X),E13.6)
 51   FORMAT(4(E13.6,1X))

      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
      REAL*8 tUnit
      COMMON /unit/ tUnit
      
 
*     constant-rate cooling
*          Temp=-(32D0-28D0)/160D0/tUnit*time+32D0
*     Natural cooling
      Temp=(28.0D0-32.0D0)*(1.0D0-DEXP(-time/310.0D0/tUnit))
     &     +32.0D0 
      RETURN
      END
*************************************************************

      REAL*8 FUNCTION Csat(T)

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

      REAL*8 T

      Csat=0.1286D0+0.00588D0*T+0.0001721D0*T**2

      RETURN
      END
****************************************************************

      REAL*8 FUNCTION growth(conc, concs)

*     growth rate for the simulation of a cooling batch
*     crystallizer
*     
*     arguments:  conc - solute concentration
*                 concs - saturation concentration
*     non-argument input: kg, g kinetic rate parameters
*     output: growth - growth rate  microns/sec

      REAL*8 conc, concs, kg, g, b, kb
      COMMON g, kg, b, kb

      growth=kg*((conc-concs)/concs)**g
      growth=0.2D0
 
      RETURN
      END
************************************************************

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

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

      REAL*8 conc, concs, m3, kb, b, g, kg
      COMMON g, kg, b, kb

*      birth=kb*((conc-concs)/concs)**b*m3
      birth=0.0D0

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

*       population density function:
*
*       f0= -0.00034786x^2 + 0.1363609*x - 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
      IF (r1 .LT. 180.5D0) THEN
          Fzero=0.0D0
      ELSE IF ( r1 .GT. 211.5D0) THEN
          Fzero=0.0D0
      ELSE
          Fzero= -0.00034786D0*r1**2 + 0.1363609D0*r1 - 13.2743D0  
      END IF
      
      RETURN
      END