compart.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.
*******************************************************************************************
* compart.f
*     This program is used to simulate batch crystallization w/o perfect or imperfect 
*     mixning.  Use mpi to run this program.  If the perfect mixing is valid, set the number
*     of process to 1; otherwise, set the number of the process equal to the number of the 
*     compartment.  Crystals have two characteristic growth dimensions.

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

*     OUTPUT: 1. moments (3 files: mu1name, mu2name, mu3name)
*             2. time (1 file: timename)
*             3. concentration, saturated concentration, supersatuation (1 file: concname)
*             4. temperature, growth rates, birthrates (1 file: tgbname)
*             5. CSD


*     IMPROTANT PARAMETERS:
*     method: upwind=1, lax-wendroff=2, high resolution=3(only used for perfect mixing case)
*     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 CSD to file
*     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.
*
*     following parameters are for mpi
*     maxf1=sub-array of array f to be passed to upper and lower neighbors, 
*           with size =maxR1R2*numcols
*     numcols=number of columns
*     iter=maxR1R2/numcols

      program main
      IMPLICIT NONE      
      include 'mpif.h'
   
*     Variables for mpi
      INTEGER i,ierr,com1d,myrank
      INTEGER numprocs,nbrtop,nbrbottom
      INTEGER status_array(MPI_STATUS_SIZE,4),req(4)
      LOGICAL isp, reorder

*     Parameters
      INTEGER totTpts,totR1pts,totR2pts,maxpoint
      INTEGER skipT,skipR1,skipR2,skipNext,skipT1
      INTEGER maxR1R2,maxf,numcols,maxf1,iter,method
      REAL*8 delT,delR1,delR2,tUnit,totTime
      REAL*8 totMass,mass_flow,alpha1,alpha2
      PARAMETER (method=3)
      PARAMETER (maxpoint=7202, totTime=7200.0D0)
      PARAMETER (delT=5.0D0, delR1=1.0D0, delR2=1.0D0)
      PARAMETER (skipR1=1, skipR2=1, skipT=1440, skipT1=288)
      PARAMETER (maxR1R2=998,iter=40,maxf=(maxR1R2+2)*(maxR1R2+2))
      PARAMETER (maxf1=25000,numcols=25)
      PARAMETER (alpha1=0.1D0, alpha2=0.1D0) 
      PARAMETER (mass_flow=10.0D0, totMass=2.0D3)

*     Integers 
      INTEGER j,k,m

*     Crystal size distribution      
      REAL*8 f(-1:maxR1R2,-1:maxR1R2)
      REAL*8 ftemp1(-1:maxR1R2,-1:maxR1R2)  
      REAL*8 ftemp2(-1:maxR1R2,-1:maxR1R2)
      REAL*8 ftemp3(-1:maxR1R2,-1:maxR1R2)
      REAL*8 ftemp4(-1:maxR1R2,-1:maxR1R2)
      REAL*8 ftemp5(-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)
      REAL*8 superConc(maxpoint), conc(maxpoint)
      REAL*8  t_conc,b_conc

*     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 b,kb,g1,kg1,g2,kg2
      REAL*8 R0, kv, rhoc, alpha,CA,CB

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

*     Initial and final temperatures (C)
      REAL*8 init_temp,final_temp

*     total simulation time (second)
      REAL*8 tot_time

*     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=120,starti3=240)
      PARAMETER(startj1=5,startj2=120,startj3=240)
      INTEGER tmpij
*
*     Functions
      REAL*8 Temp, Csat, growth1, growth2, birth, Fzero 

*     Compartment
      REAL*8 Wu, Wd, tau, mass_per_comp

*     file names and units
      CHARACTER*10 fname,mu1name,mu2name,mu3name,tgbname
      CHARACTER*11 timename,concname
      CHARACTER*10 logname
      INTEGER nfile, nUnit

*     Initialization
      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 ftemp1/maxf*0.0D0/
      DATA ftemp2/maxf*0.0D0/
      DATA ftemp3/maxf*0.0D0/
      DATA ftemp4/maxf*0.0D0/
      DATA ftemp5/maxf*0.0D0/

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

*     Initialize topologies
      isp=.false.
      reorder=.true.
      call MPI_INIT(ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr)
      call MPI_CART_CREATE(MPI_COMM_WORLD,1,numprocs,isp,reorder,
     &                     com1d,ierr)
      call MPI_COMM_RANK(com1d,myrank,ierr)
      call MPI_CART_SHIFT(com1d,0,1,nbrbottom,nbrtop,ierr)

*     Initialize the file names
      IF(numprocs.GT.1)THEN
         fname='fxx_xx.dat'
         mu1name='mu1_xx.dat'
         mu2name='mu2_xx.dat'
         mu3name='mu3_xx.dat'
         timename='time_xx.dat'
         concname='conc_xx.dat'
         tgbname='tgb_xx.dat'
         logname='log_xx.dat'
      ELSE
         IF(method.EQ.1)THEN
            fname='fxx_xx.dat'
            mu1name='upwmu1.dat'
            mu2name='upwmu2.dat'
            mu3name='upwmu3.dat'
            timename='upwtime.dat'
            concname='upwconc.dat'
            tgbname='upwtgb.dat'
            logname='log_xx.dat'
         ELSE IF(method.EQ.2)THEN
            fname='fxx_xx.dat'
            mu1name='laxmu1.dat'
            mu2name='laxmu2.dat'
            mu3name='laxmu3.dat'
            timename='laxtime.dat'
            concname='laxconc.dat'
            tgbname='laxtgb.dat'
            logname='log_xx.dat'
          ELSE
            fname='fxx_xx.dat'
            mu1name='higmu1.dat'
            mu2name='higmu2.dat'
            mu3name='higmu3.dat'
            timename='higtime.dat'
            concname='higconc.dat'
            tgbname='higtgb.dat'
            logname='log_xx.dat'
         END IF
      END IF
*     Conversion factor: tUnit=60=sec, tUnit=1=min      
      tUnit=60D0
*     Save CSD for every skipT seconds
      skipNext=skipT
*     Mesh boundary
      meshi1=DINT(DBLE(starti1)/delR1)
      meshi2=DINT(DBLE(starti2)/delR1)
      meshi3=DINT(DBLE(starti3)/delR1)
      meshj1=DINT(DBLE(startj1)/delR2)
      meshj2=DINT(DBLE(startj2)/delR2)
      meshj3=DINT(DBLE(startj3)/delR2)

      
*     Growth and nucleation kinetic parameters (from Rudi)
      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 and final temperature (C)
      init_temp=33.0D0
      final_temp=28.0D0
*     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
*     mass for each compartment 
      mass_per_comp=totMass/dble(numprocs)
*     tau=F/V
      tau=mass_flow/mass_per_comp
*     Total time and mesh discretization pts
      totTpts=INT(totTime/delT)
      totR1pts=INT(DBLE(maxR1R2)/delR1) 
      totR2pts=INT(DBLE(MAXR1R2)/delR2)
      tot_time=totTime
*     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 f(r1,r2,0)      
      IF(numprocs.GT.1) THEN
         IF(myrank.EQ.0) THEN
            DO i=1, totR1pts
               DO j=1, totR2pts
                  f(i,j)=Fzero(delR1*dble(i), delR2*dble(j))*
     &                 DBLE(numprocs)
               ENDDO
            ENDDO
         END IF
      ELSE
         DO i=1, totR1pts
            DO j=1, totR2pts
               f(i,j)=Fzero(delR1*dble(i), delR2*dble(j))
            ENDDO
         ENDDO         
      END IF
*     r2=j=length, r1=i=width
*     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)