conf2.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.
*
*       conf2.f
*
*	This program calculates the parameter covariance matrix and 
* 	the variance for each growth and nucleation kinetic 

*     parameter (g1, kg1, g2, kg2, b, and kb).
*
*	This program reads the growth and nucleation 
*     kinetic parameters (g1, kg1, g2, kg2, b, and kb), 
*	concentration, and moments data from data files.
*
*
*       Nset=the number of the data files
*
*       Date:    August 8, 1998
*       Authors: Serena H. Chung and Richard D. Braatz
*                Department of Chemical Engineering
*                University of Illinois at Urbana-Champaign

*	  Modified:March 6, 2000

*	  By:      David L. Ma and Richard D. Braatz

*	  

*	  Modifications were to replace the use of transmittance

*       measurements with moment measurements.
*
*
*



	use MSIMSL
*      PROGRAM MAIN

      INTEGER NSTEP, NEQ, MXPARM, NTHETA, NU,Nsets
      INTEGER npos1, npos2, MSTEP, Cinter, Minter
      PARAMETER (NSTEP=161,NEQ=91,MXPARM=50,NTHETA=6, NU=12)	 
      PARAMETER (Nsets=1, npos1=12, npos2=8, MSTEP=5)
      PARAMETER (Cinter=1, Minter=30)
      INTEGER I, J, K, M, N,Iteration, I1, L, P
      INTEGER Norder, LDA, LDB, IPATH
      PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
	PARAMETER(LDFAC=6, P=6)
	INTEGER IPVT(P)
      REAL*8 T, Y(NEQ)
      REAL*8 F(NU*(NSTEP-1)*Nsets, NTHETA)
      REAL*8 FTVF(NTHETA, NTHETA)
      REAL*8 inv_FTVF(NTHETA,NTHETA)
      REAL*8 FTV(NTHETA,NU*(NSTEP-1)*Nsets)
	REAL*8 DET1, DET2, FAD(LDFAC,LDFAC)
      REAL*8 delt, tfinal, Msolv
      REAL*8 mu00
      REAL*8 cell_length, kv, ka, densityc, densitys
      REAL*8 r0, alpha, g1, kg1, g2, kg2, b, kb
*      REAL*8 weight_mean_size, cov, mass_ratio
      
      REAL*8 moment00(NSTEP), moment10(NSTEP), moment01(NSTEP)
      REAL*8 moment11(NSTEP), moment20(NSTEP), moment02(NSTEP)
      REAL*8 moment21(NSTEP), moment12(NSTEP), moment22(NSTEP)
      REAL*8 moment30(NSTEP), moment31(NSTEP)

	REAL*8 time (NSTEP), concentration(NSTEP) 
      REAL*8 seed_moment10(NSTEP), seed_moment01(NSTEP)
      REAL*8 seed_moment11(NSTEP), seed_moment20(NSTEP)
      REAL*8 seed_moment02(NSTEP), seed_moment21(NSTEP) 
      REAL*8 seed_moment12(NSTEP)
	REAL*8 temperature(NSTEP), relsatn(NSTEP)
      REAL*8 Temp, Csat, detFTVF
      REAL*8 derivtheta(NSTEP,NU,NTHETA)
      
	REAL*8 ConcData(NSTEP-1), Mu00Data(MSTEP)
      REAL*8 Mu10Data(MSTEP), Mu01Data(MSTEP)
      REAL*8 Mu11Data(MSTEP), Mu20Data(MSTEP)
      REAL*8 Mu02Data(MSTEP), Mu21Data(MSTEP)
      REAL*8 Mu12Data(MSTEP), Mu22Data(MSTEP)
	REAL*8 Mu30Data(MSTEP), Mu31Data(MSTEP)
	
      REAL*8 conc_variance, mu00_variance, mu01_variance
      REAL*8 mu10_variance, mu11_variance, mu20_variance
      REAL*8 mu02_variance, mu21_variance, mu12_variance
	REAL*8 mu22_variance, mu30_variance, mu31_variance
	REAL*8 chi_squared
      REAL*8 g1_interval, lnkg1_interval 
	REAL*8 g2_interval, lnkg2_interval
      REAL*8 b_interval, lnkb_interval
*      REAL*8 AA(3,3), BB(3), gamma1(3), lmin, lmax
      REAL*8 mass_seed, mean_seed, width
      CHARACTER*60 file_concen, file_mu


*lsodes' parameters
      
      INTEGER itol, iopt, itask, istate, mf
      INTEGER lrw, liw, iwork(1000)
      REAL*8  rtol, atol, rwork(5800)    

      EXTERNAL FCN, FCNJ

      COMMON /GROWTH_DATA1/kg1, g1
	COMMON /GROWTH_DATA2/kg2, g2
      COMMON /BIRTH_DATA/kb, b 
      COMMON /EXP_DATA/r0, alpha, mu00
      COMMON Iteration


	file_concen='concen_data'

	file_mu='mu_data'

	OPEN(UNIT=90, FILE='param.dat', FORM='FORMATTED',

     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')

	READ(90,*)g1

	READ(90,*)kg1

	READ(90,*)g2
     
      READ(90,*)kg2 
	
      READ(90,*)b

	READ(90,*)kb

	kg1=DEXP(kg1)

	kg2=DEXP(kg2)
     
      kb=DEXP(kb)

	
	print*, "running"
      DO 787 Iteration = 1, Nsets

*Simulation parameters
*     controller time step in minutes
      delt = 1.0D0
*     final time in minutes
      tfinal = DFLOAT(NSTEP)*delt

*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

*     Total mass of seed crystals (grams)	
      mass_seed = 230.0D0
*     Mean size of the seed crystals (microns)
      mean_seed = 196.0D0
*     Percent width of the initial seed distribution
      width = 6.12D0


*Growth and nucleation kinetic parameters 
*     (dimensionaless)
*      g = 0.12889594028434D+01
*      g = 0.13151365743328D+01
*      g = 0.13125608981096D+01
*      g = 0.12912359265617D+01
*      g = 0.12993596644840D+01
*     (mirons/minute)
*      kg = DEXP(0.87458112413202D+01)
*      kg = DEXP(0.88354817511395D+01)
*      kg = DEXP(0.88246805512230D+01)
*      kg = DEXP(0.87309813411979D+01)
*      kg = DEXP(0.87620506647544D+01)
*     (dimensionless)
*      b = 0.17519140794580D+01
*      b = 0.17762965268414D+01
*      b = 0.17704343279926D+01
*      b = 0.17465821622552D+01
*      b = 0.17860382260529D+01

*     (number of particles/cm^3/minute) 
*     (the units have been corrected from that reported in
*     Table 3.1 in Miller)
*      kb = DEXP(0.17048465188505D+02)
*      kb = DEXP(0.17137238240957D+02)
*      kb = DEXP(0.17106281255156D+02)
*      kb = DEXP(0.17005247110939D+02)
*      kb = DEXP(0.17179861650138D+02)

*Initial conditions
*     initial concentration, g/g solvent
      concentration(1) = 0.493D0


*
*       The initial moments were computed using the following
*       population density function:
*
*       f0= -0.00034786x^2 + 0.1363609*x - 13.2743
*	    -0.00034786y^2 + 0.1363609*y - 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=121.575D0
      moment00(1)=mu00
           
*     initial moment10, moment01 micron/g solvent
      moment10(1)=2.383D4
      moment01(1)=moment10(1)
      seed_moment10(1)=moment10(1)
      seed_moment01(1)=moment10(1)
      
*     initia1 moment11, micron^2/g solvent
      moment11(1)=4.67D6
      seed_moment11(1)=moment11(1)
      
*     initial moment20, moment02 micron^2/g solvent
      moment20(1)=4.679D6
      moment02(1)=moment20(1)
      seed_moment20(1)=moment20(1)
      seed_moment02(1)=moment20(1)
      
*     initial moment12, moment21 micron^3/g solvent
      moment12(1)=9.17D8
      moment21(1)=moment12(1)
      seed_moment21(1)=moment12(1)
      seed_moment12(1)=moment12(1)
      
*     initial moment22 micron^4/g solvent
      moment22(1)=1.801D11
      
*     initial moment30 micron^3/g solvent
      moment30(1)=9.203D8

*     initial moment31 micron^4/g solvent
      moment31(1)=1.801D11

*     initial relative supersaturation
      relsatn(1)=(concentration(1)-Csat(Temp(0.0D0)))/
     &     Csat(Temp(0.0D0))



      conc_variance=3.53D-4
      mu00_variance=(0.2D0*moment00(1))**2.0D0
      mu10_variance=(0.2D0*seed_moment10(1))**2.0D0
      mu01_variance=(0.2D0*seed_moment01(1))**2.0D0 
      mu11_variance=(0.2D0*seed_moment11(1))**2.0D0
      mu20_variance=(0.2D0*seed_moment20(1))**2.0D0
	mu02_variance=(0.2D0*seed_moment02(1))**2.0D0
      mu21_variance=(0.2D0*seed_moment21(1))**2.0D0
      mu12_variance=(0.2D0*seed_moment12(1))**2.0D0
      mu22_variance=(0.2D0*moment22(1))**2.0D0
      mu30_variance=(0.2D0*moment30(1))**2.0D0
	mu31_variance=(0.2D0*moment31(1))**2.0D0	


*     Initial conditions for derivatives wrt to theta
      DO 163 I = 1 , NU
         DO 164 J = 1 , NTHETA
            derivtheta(1,I,J)=0.0D0
164      CONTINUE
163   CONTINUE     


*     Read Data in files
      I1=10+Iteration
      DO  J = 1,NSTEP-1
         WRITE(file_concen(npos1:npos1+1), '(I2)') I1 

	   file_concen = file_concen(1:npos1+1) // '.txt'
         OPEN(UNIT=20, FILE=file_concen, FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='OLD')
         READ(20,*) ConcData(J)	
      ENDDO
      CLOSE(UNIT=20)
      DO J=1,MSTEP
         WRITE(file_mu(npos2:npos2+1), '(I2)') I1 

	   file_mu = file_mu(1:npos2+1) // '.txt'
         OPEN(UNIT=20, FILE=file_mu, FORM='FORMATTED',
     &        ACCESS='SEQUENTIAL', STATUS='OLD')
	 READ(20,*) Mu00Data(J), Mu10Data(J), Mu01Data(J),
     &		    Mu11Data(J), Mu20Data(J), Mu02Data(J),
     &		    Mu21Data(J), Mu12Data(J), Mu22Data(J),
     &		    Mu30Data(J), Mu31Data(J)

      ENDDO 
      CLOSE(UNIT=20)


*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=3800
      liw=200
      rtol=1.0d-5
      atol=1.0d-5
      itol=1
******************************************************
*
      time(1)=0.0D0
      T=0.0D0    
      L=1
        
      DO 1200 I=1,(NSTEP-1)
	 temperature(I)=Temp(T)
         Y(1)=moment00(I)
         Y(2)=moment10(I)
         Y(3)=moment01(I)
         Y(4)=moment11(I)
         Y(5)=moment20(I)
         Y(6)=moment02(I)
         Y(7)=moment21(I)
         Y(8)=moment12(I)
         Y(9)=moment22(I)
         Y(10)=moment30(I)
         Y(11)=moment31(I)
         Y(12)=concentration(I)
         Y(13)=seed_moment10(I)
         Y(14)=seed_moment01(I)
         Y(15)=seed_moment11(I)
         Y(16)=seed_moment20(I)
         Y(17)=seed_moment02(I)
         Y(18)=seed_moment21(I)
         Y(19)=seed_moment12(I)

	          K = 20
         DO 263 M = 1 , NU
            DO 264 N = 1 , NTHETA
               Y(K)=derivtheta(I,M,N)
*			 print*, "k=", K
*			 print*, "Y(K)=", Y(K)
	       K = K + 1
264         CONTINUE
263      CONTINUE
*	   pause

*	print*, 'before lsodes'
*	print*, 'istate before', istate
        CALL lsodes ( FCN,NEQ,Y,T,T+Cinter,itol,rtol,
     &          atol,itask,istate,iopt, 
     &           rwork,lrw,iwork,liw, FCNJ, mf )
     
*	print*, 'after lsodes'