sim2d.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.
*******************************************************************************************
*If you use this program, please cite the following papers:
*
* 1. D. L. Ma, D. Tafti, and R. D. Braatz. Compartmental modeling of multidimensional   
*    crystallization. Special Issue on Crystallization and Interfacial Processes, 
*    Int. J. of Modern Physics B, 2002. in press. 
* 2. D. L. Ma, D. K. Tafti, and R. D. Braatz. High resolution simulation of multidimensional    *    crystallization. Special Issue in Honor of William R. Schowalter, Ind. Eng. Chem. Res., *    2002. in press.
* 3.  D. L. Ma, D. K. Tafti, and R. D. Braatz. Optimal control and simulation of 
*     multidimensional crystallization. Special Issue on Distributed Parameter Systems, 
*     Comp. & Chem. Eng., 2002. in press.
*******************************************************************************************
* sim2d.f
*
*   This program is used to simulate batch crystallizer under perfect mixing assumption.  Crystals have two characteristic growth dimensions. 

*     INPUT:  1. kinetic parameters.  These parameters are set manually.
*             

*     OUTPUT: 1. moments (3 files: mu1name, mu2name, mu3name, prefixed by the method's name)
*             2. time (1 file: timename)
*             3. concentration, saturated concentration, supersatuation (1 file: concname)
*             4. temperature, growth rates, birthrates (1 file: tgbname)
*             5. crystal size distribution (1 file: filename)


*     IMPROTANT PARAMETERS:
*     method: upwind=1, lax-wendroff=2, high resolution=3
*     totTime= total simulation time
*     delT=time step
*     delR1=delR2=steps in R1 and R2 directions
*     maxpoint=size the of the array which shold be larger than ceil(totTime/delT)
*     skipT=every skipT seconds, write the crystal size distribution f to file (number of 
*            CSD files=totTime/delT/skipT)
*     skipT1=every skipT1 seconds, write temperature, concentration, etc to the file
*     maxR1R2= MAX(R1,R2)
*     maxf=mesh size=(maxR1R2+2)**2
*     totMass=mass of the solvent.


      PROGRAM Main
      IMPLICIT NONE   

*     Parameters
      INTEGER maxpoint,totTpts,totR1pts,totR2pts
      INTEGER skipT,skipR1,skipNext,skipR2
      INTEGER maxR1R2,skipT1,method
      REAL*8 delT,delR1,delR2,tUnit
      REAL*8 totTime
*     maxf=maxR1R2*maxR1R2*2
      INTEGER maxf
      PARAMETER ( maxpoint=7202, totTime=7201.0D0)
      PARAMETER (delT=5D0, delR1=1D0, delR2=1D0)
      PARAMETER (skipR1=1, skipR2=1, skipT=1440, skipT1=144)
      PARAMETER (maxR1R2=898, maxf=810000,method=3)

*     Integers
      INTEGER i, j, k, m  
*
*     Crystal size distribution      
      REAL*8 f(-1:maxR1R2,-1:maxR1R2)
      REAL*8 ftemp(-1:maxR1R2,-1:maxR1R2)  

*     Moments
      REAL*8 mu00(maxpoint),mu01(maxpoint)
      REAL*8 mu10(maxpoint),mu11(maxpoint)
      REAL*8 mu20(maxpoint),mu02(maxpoint)
      REAL*8 mu21(maxpoint),mu12(maxpoint)
      REAL*8 mu22(maxpoint),mu30(maxpoint)
      REAL*8 mu31(maxpoint),volume(maxpoint)

*     Concentrations
      REAL*8 Concs(maxpoint),Conc(maxpoint)
      REAL*8 superConc(maxpoint)

*     Growth rate, birth rate and temperature
      REAL*8 time(maxpoint),T(maxpoint)
      REAL*8 G1rate(maxpoint), G2rate(maxpoint)
      REAL*8 Brate(maxpoint)

*     Physcial parameters
      REAL*8 g1,kg1,g2,kg2,b,kb
      REAL*8 CA,CB,R0,kv,rhoc,alpha

*     Constants for high resolution method
      REAL*8 vei,vej,veij
      REAL*8 i1theta,i2theta,j1theta,j2theta
      REAL*8 i1phi,i2phi,j1phi,j2phi
      

*     Mesh boundary
      INTEGER meshi1,meshi2,meshi3
      INTEGER meshj1,meshj2,meshj3
      INTEGER starti1,starti2,starti3
      INTEGER startj1,startj2,startj3
      REAL*8 mratei, mratej
      PARAMETER(starti1=5,starti2=150,starti3=250)
      PARAMETER(startj1=5,startj2=150,startj3=250)
*
*     Functions
      REAL*8 Temp,Csat,growth1,growth2,birth,Fzero 

*     File names
      CHARACTER*9 filename
      INTEGER nfile
*
      DATA mu00/maxpoint*0.0d0/
      DATA mu10/maxpoint*0.0d0/
      DATA mu01/maxpoint*0.0D0/
      DATA mu11/maxpoint*0.0D0/
      DATA mu20/maxpoint*0.0D0/
      DATA mu02/maxpoint*0.0D0/
      DATA mu21/maxpoint*0.0D0/
      DATA mu12/maxpoint*0.0D0/
      DATA mu22/maxpoint*0.0D0/
      DATA mu30/maxpoint*0.0D0/
      DATA mu31/maxpoint*0.0D0/
      DATA volume/maxpoint*0.0D0/
      DATA f/maxf*0.0D0/
      DATA ftemp/maxf*0.0D0/

*     Common blocks
      COMMON /unit/ tUnit
      COMMON /g1data/ g1, kg1
      COMMON /g2data/ g2, kg2
      COMMON /rhoc/rhoc
      COMMON /Csat_DATA/CA, CB
      COMMON /bdata/ kb, b

      PRINT*, "I am running"
      filename='popxx.dat'      
*     Conversion factor: tUnit=60=sec, tUnit=1=min      
      tUnit=60D0
*     Save f for every skipT seconds
      skipNext=skipT
*     Meshboundary
      meshi1=DINT(5.0/delR1)
      meshi2=DINT(120.0D0/delR1)
      meshi3=DINT(280.0D0/delR1)
      mratei=0.0D0
      meshj1=DINT(5.0/delR2)
      meshj2=DINT(120.0D0/delR2)
      meshj3=DINT(280.0D0/delR2)
      mratej=0.0D0
      
*     Growth and nucleation kinetic parameters
      b=.2045051881054735D+01
*     #/micron^3/sec/g-solvent
      kb=DEXP(.1531801881484984D+02)/tunit*1.0D-12
      g1=.1477587696083403D+01
*     micron/sec
      kg1=DEXP(.6596298446535068D+01)/tunit
      g2=.1741171935425851D+01
*     micron/sec
      kg2=DEXP(.8706995338237306D+01)/tunit        

*     Initial size of crystal
      R0=0.0D0       
*     volume shape factor 
      kv=1.0D0
*     density of crystal in kg/m^3 
      rhoc=2338D0
*     conversion factor g/micron^3
      alpha=kv*rhoc*(1.0D-6)**3*1.0D3  
*
*     Total time and mesh discretization pts
      totTpts=int(totTime/delT)
      totR1pts=int(maxR1R2/delR1)
      totR2pts=int(maxR1R2/delR2)

      print*, "totTime=", totTime
      print*, "delT=", delT
      print*, "int(tot/T)=", int(totTime/delT)
      print*, "totTpts=", totTpts
      print*, "totR1pts=", totR1pts
      print*, "totR2pts=", totR2pts

*     Initialization
*     Initial time
      time(1)=0.0D0
*     Concentration g/g solvent
      Conc(1)= 0.307D0
*     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)
*     Initial condition at f(r1,r2,0)
      DO i=1, totR1pts
         DO j=1, totR2pts
            f(i,j)=Fzero(delR1*dble(i), delR2*dble(j))
         ENDDO
      ENDDO      
*     Initial moments
       DO i=1,meshi1
        DO j=1,meshj1  
         mu00(1)=mu00(1)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))/4.0D0
     &           *delR1*delR2
	 mu10(1)=mu10(1)
     &    +(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(1)=mu01(1)
     &   +(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(1)=mu11(1)
     &   +(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(1)=mu20(1)
     &   +(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(1)=mu02(1)
     &   +(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(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
       DO i=meshi2,meshi3
        DO j=meshj2,meshj3  
         mu00(1)=mu00(1)+(f(i,j)+f(i,j+1)+f(i+1,j)+f(i+1,j+1))/4.0D0
     &           *delR1*delR2
	 mu10(1)=mu10(1)
     &    +(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(1)=mu01(1)
     &   +(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(1)=mu11(1)
     &   +(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(1)=mu20(1)
     &   +(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(1)=mu02(1)
     &   +(f(i,j)*dble(j-1)**2
     &   +f(i,j+1)*dble(j)**2