dim1.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.
*       dim1.f
*
*       This program simulates the operation of an industrial
*       crystallizer scaled-up from the experimental batch cooling
*       crystallizer described in S. M. Miller's Ph.D. thesis
*       published at the University of Texas at Austin in 1993.
*	  The KNO3 crystals have one characteristic size.

*
*       The program reads the temperature profile from the
*       FORTRAN function Temp.  The output of the
*       program are the moments, the solute concentration, the
*       relative supersaturation, temperature, transmittance,
*       and values for three common final crystal properties of
*       interest: (i) weight mean size, (ii) ratio of nucleated
*       crystal mass to seed crystal mass, and (iii) coefficient
*       of variation.
*
*	Output files:
*

*	moments.dat contains:
*		time, moment0, moment1, moment2, and moment3
*	seed.dat contains:
*		time, seedmoment1,seedmoment2,seedmoment3
*	cond.dat contains:
*		time, concentration, temperature, relsatn, transmittance 	
*	obj.dat contains the final objective values:
*		weight_mean_size, cov, mass_ratio
*
*
*	Main parameter inputs are the growth and nucleation 
*       kinetic parameters (g, kg, b, and kb)
*
*	Parameter NN:  the number of discretizations
*	Parameter NEQ: the number of moment equations
*
*
*       Date:    January 27, 1998
*       Authors: Serena H. Chung and Richard D. Braatz
*                Department of Chemical Engineering
*                University of Illinois at Urbana-Champaign

*

*	  Date modified:  February 19, 2000

*       Modified by:    Richard D. Braatz and David L. Ma

*					  Department of Chemical Engineering

*					  University of Illinois at Urbana-Champaign
*
*
*
*	USE MSIMSL
*      PROGRAM MAIN
      INTEGER NN, NEQ, MXPARM
      PARAMETER (NN=161, NEQ=9, MXPARM=50)
      INTEGER IDO, flag, kfinal, I

*	INTEGER IUNIT
      REAL*8 A(NEQ,NEQ), TOL, T, TEND, Y(NEQ)
      REAL*8 PARAM(MXPARM)
      REAL*8 delt, tfinal
      REAL*8 mu00
      REAL*8 cell_length, Msolv, kv, ka, UA, densityc, densitys
      REAL*8 r0, alpha, g, kg, b, kb
      REAL*8 weight_mean_size, cov, mass_ratio
      REAL*8 moment0(NN), moment1(NN), moment2(NN)
      REAL*8 moment3(NN), moment4(NN)
      REAL*8 time (NN), concentration(NN), seed_moment1(NN)
      REAL*8 concentration_measured(NN)
      REAL*8 seed_moment2(NN), seed_moment3(NN), transmittance(NN)
      REAL*8 temperature(NN), relsatn(NN)
      REAL*8 Temp, Csat
      REAL*8 betatrans, betaconc
*      REAL real_time(NN), real_temperature(NN), real_concentration(NN)
*      REAL real_moment0(NN), real_moment1(NN), real_moment2(NN)
*      REAL real_moment3(NN), real_transmittance(NN), real_satn(NN)
	

*lsodes' parameters

      

      INTEGER itol, iopt, itask, istate, mf

      INTEGER lrw, liw, iwork(30)

      REAL*8  rtol, atol, rwork(900)


      EXTERNAL FCN, FCNJ
      COMMON /DATA1/flag
      COMMON /GROWTH_DATA/kg, g
      COMMON /BIRTH_DATA/kb, b 
      COMMON /EXP_DATA/r0, alpha, mu00, UA, Msolv







*Simulation parameters

*********************************************************



      mf=222

      itask=1

      istate =1

      iopt=1

	rwork(5) = 1.1D-14

	rwork(6) = 1.0D-1

	rwork(7) = 1.0D-14

	iwork(6) = 100000

      lrw=900

      liw=90

      rtol=1.0d-4

      atol=1.0d-6

      itol=1


*Simulation parameters
*     controller time step in minutes
      delt = 1.0D0
*     total time step
      kfinal=NN
*     final time in minutes
      tfinal = DFLOAT(NN)*delt
*     noise for transmittance measurement
      betatrans = 0.0003D0
      betatrans = 0.009D0
*     noise for concentration measurement
      betaconc = 0.0001D0
      betatrans = 0.005D0

*Initial conditions
*     initial concentration, g/g solvent
      concentration(1) = 0.493D0
      concentration_measured(1)=concentration(1)
*     
*
*       The initial moments were computed using the following
*       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-212 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.
*
*     initial zeroth moment, number of particle/g solvent
      mu00=1.8996D0
      moment0(1)=mu00
*     initial first moment, micron/g solvent
      moment1(1)=372.322D0
      seed_moment1(1)=moment1(1)
*     initia1 second moment, micron^2/g solvent
      moment2(1)=73072.2D0
      seed_moment2(1)=moment2(1)
*     initial third moment, micron^3/g solvent
      moment3(1)=1.43602D7
      seed_moment3(1)=moment3(1)
*     initial fourth moment, micron^4/g solvent
      moment4(1)=2.82582D9

*     initial relative supersaturation
      relsatn(1)=(concentration(1)-Csat(Temp(0.0D0)))/
     &     Csat(Temp(0.0D0))
*     initial transmittance measurement
      transmittance(1)=DEXP(-ka/2D0*cell_length/10D0*moment2(1)*
     &     (densitys*(1D-4)**2))
  

*Parameters for experimental set-up
*     cell length for spectrophotometer in millimeter
*     This was modified from that in Miller's thesis because
*     his value (2.0) did not agree with his simulation results
      cell_length=1.77D0
*     mass of solvent in grams, converted from 2000 gallons
      Msolv=7.57D6
*     volume shape factor (Appendix C in Miller)
      kv=1.0D0
*     area shape factor (Appendix C in Miller)
      ka=6.0D0
*     heat transfer coefficient multiplied by surface area
*     in calorie/minute/degree C 
*     density of solvent in g/cm^3 (solvent is water)
      densitys=1.0D0
*     density of crystal in g/cm^3 (Appendix C in Miller)
      densityc=2.11D0
*     seed size at nucleation
      r0=0.0D0
*     crystal density*volume shape factor,
*     in gram crystal/micron^3/particle
*     (alpha*L^3=mass of particle)
      alpha=kv*densityc*(1.0D-4)**3

*Growth and nucleation kinetic parameters (Table 4.6 in Miller)
*     (dimensionaless)
      g=1.32D0
 
*     (mirons/minute)
      kg=DEXP(8.749D0)

*     (dimensionless)
      b=1.78D0

*     (number of particles/cm^3/minute) 
*     (the units have been corrected from that reported in
*     Table 3.1 in Miller)
      kb=DEXP(17.742D0)


	 PRINT*, '   '
       PRINT*, 'Which temperature cooling profile would you like?'

	 PRINT *, '  '
       PRINT*, ' - enter 1 for a constant cooling rate'
       PRINT*, ' - enter 2 for natural cooling'
       PRINT*, ' - enter 3 for T(t) specified in function Temp'
       PRINT*, '  '
       READ*, flag
       IF ((flag.NE.1).AND.(flag.NE.2).AND.(flag.NE.3)) THEN
          PRINT*, 'The temperature profile defaulted to a'
          PRINT*, 'constant cooling rate.'

	    PRINT*, '    '
          flag=1
       ENDIF

	 PRINT*, '  '

	 PRINT*, 'Results of time domain simulation:'

	 PRINT*, '  '

*     Parameters for solving the ODE
      TOL=1.00D-2
      PARAM(4)=1.0D6
*     Using Gear's BDF method for integration
      PARAM(12)=2
*     Eror Norm;  Absolute error - max(abs(ei))
      PARAM(10)=1
*     Nonlinear solve method indicator: chord method with a
*     divided-difference Jacobian
      PARAM(13)=2
*     Initial time step Hinit
      PARAM(1)=1.0D-9
*     Maximum time step Hmax
*      PARAM(3)=1.0D-3
      PARAM(3)=1.0D-4
*     Minimum time step Hmin
      PARAM(2)=0.0D0
*     Intrp1
      PARAM(7)=0.0D0
*     Intrp2
      PARAM(8)=0.0D0
      time(1)=0.0D0
	temperature(1)=Temp(time(1))

	OPEN(UNIT=31, FILE='dim1con', STATUS='UNKNOWN')

	OPEN(UNIT=32, FILE='dim1mom', STATUS='UNKNOWN')
      DO 1200 I=1,(NN-1)
         IDO=1
         T=time(I)
         TEND=T+delt
         Y(1)=moment0(I)
         Y(2)=moment1(I)
         Y(3)=moment2(I)
         Y(4)=moment3(I)
         Y(5)=moment4(I)
         Y(6)=concentration(I)
         Y(7)=seed_moment1(I)
         Y(8)=seed_moment2(I)
         Y(9)=seed_moment3(I)


         CALL lsodes ( FCN,NEQ,Y,T,TEND,itol,rtol,

     &          atol,itask,istate,iopt, 

     &           rwork,lrw,iwork,liw,FCNJ, mf )


         moment0(I+1)=Y(1)
         moment1(I+1)=Y(2)
         moment2(I+1)=Y(3)
         moment3(I+1)=Y(4)
         moment4(I+1)=Y(5)
         concentration(I+1)=Y(6)
         concentration_measured(I+1)=concentration(I+1)
         seed_moment1(I+1)=Y(7)
         seed_moment2(I+1)=Y(8)
         seed_moment3(I+1)=Y(9)
         relsatn(I+1)=(concentration(I+1)-Csat(Temp(T)))/
     &        Csat(Temp(T))
	   temperature(I)=Temp(T)
         transmittance(I+1)=DEXP(-ka/2.0D0*cell_length/10.0D0*
     &        moment2(I+1)*(densitys*(1.0D-4)**2.D0))
         time(I+1)=T