📄 dim2.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.
*
* dim2.f
*
* This program simulates the operation of an industrial
* crystallizer with two characteristic dimensions. The
* solubility, heat of crystallization, nucleation kinetics,
* and growth kinetics for one dimension are from Miller's
* 1993 thesis in chemical engineering from U Texas at Austin.
* The growth kinetics for the other direction were selected
* for illustrating two-dimensional crystal growth.
*
* The program reads the temperature profile from the
* FORTRAN function Temp. The output of the
* program are plots of the moments, the solute concentration,
* the relative supersaturation, temperature, transmittance,
* and values for several product quality variables
* appropriate for two-dimensional crystal growth.
*
* Date: January 27, 1998
* Authors: Serena H. Chung and Richard D. Braatz
* Modified: February 19, 2000
* By: David L. Ma and Richard D. Braatz
* Modified: September 27, 2000
* By: Chris Lentz, Rudiyanto Gunawan, and Richard D. Braatz
* Department of Chemical Engineering
* University of Illinois at Urbana-Champaign
*
* Output files:
* mu1.dat: contains moments
* moment00(I), moment10(I), moment01(I), moment11(I)
* mu2.dat: contains moments
* moment20(I), moment02(I), moment21(I), moment12(I)
* mu3.dat: contains moments
* moment22(I), moment30(I), moment31(I)
*
* obj.dat: contains objectives defined in the code
* time(I), avg_r1(I), avg_r2(I), cv_r1(I),
* cv_r2(I), area_ratio(I), tot_surf(I),
* avg_vol(I), mass_ratio(I), tot_mass(I),
* wt_size_r1(I), wt_size_r2(I)
* seed.dat: seed moments
* time(I), seed_moment10(I),
* seed_moment01(I),seed_moment11(I),
* seed_moment20(I),seed_moment02(I),
* seed_moment21(I), seed_moment12(I)
*
* conc.dat: contains
* temperature(I),concentration(I), relsatn(I)
*
* time.dat: records time
*
PROGRAM MAIN
IMPLICIT NONE
*Parameters
INTEGER NN, NEQ, MXPARM
PARAMETER (NN=161, NEQ=19, MXPARM=50)
INTEGER flag, I
REAL*8 T, Y(NEQ)
* REAL*8 PARAM(MXPARM)
REAL*8 delt
REAL*8 mu00
REAL*8 cell_length, Msolv, kv, ka, UA, densityc, densitys
REAL*8 r0, alpha, g1, kg1, g2, kg2, b, kb
*Moments
REAL*8 moment00(NN), moment10(NN), moment01(NN)
REAL*8 moment11(NN), moment20(NN), moment02(NN)
REAL*8 moment21(NN), moment12(NN), moment22(NN)
REAL*8 moment30(NN), moment31(NN)
REAL*8 seed_moment10(NN), seed_moment01(NN)
REAL*8 seed_moment11(NN), seed_moment20(NN)
REAL*8 seed_moment02(NN), seed_moment21(NN)
REAL*8 seed_moment12(NN)
*Objectives
REAL*8 avg_r1(NN), avg_r2(NN), area_ratio(NN), cv_r1(NN)
REAL*8 mass_ratio(NN), avg_vol(NN), tot_surf(NN),cv_r2(NN)
REAL*8 wt_size_r1(NN), wt_size_r2(NN), tot_mass(NN)
*Other variables
REAL*8 time (NN), concentration(NN)
REAL*8 concentration_measured(NN)
REAL*8 temperature(NN), relsatn(NN)
REAL*8 betatrans, betaconc
* REAL*8 transmittance(NN)
* REAL real_time(NN), real_temperature(NN)
* REAL real_moment0(NN), real_moment1(NN), real_moment2(NN)
* REAL real_moment3(NN), real_transmittance(NN), real_satn(NN)
* REAL*8 real_concentration(NN)
*functions
REAL*8 Temp, Csat
*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/kg1, g1
COMMON /GROWTH_DATA2/kg2, g2
COMMON /BIRTH_DATA/kb, b
COMMON /EXP_DATA/r0, alpha, mu00, UA, Msolv
OPEN(UNIT=10, FILE='mu1.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=11, FILE='mu2.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=12, FILE='mu3.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=15, FILE='obj.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=20, FILE='seed.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=25, FILE='conc.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=30, FILE='time.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
*Simulation parameters
*********************************************************
mf=222
itask=1
istate =1
iopt=1
rwork(6)=0.001D0
rwork(5)=0.001D0
iwork(6)=1000000
lrw=900
liw=90
rtol=1.0d-5
atol=1.0d-5
itol=1
*Constant data
**************************************************
* 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
* -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))
* 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)
g1=1.32D0
* (dimensionless)
kg1=DEXP(8.849D0)
* (microns/minute)
g2=1.32D0
* (dimensionless)
kg2=DEXP(8.849D0)
* (microns/minute)
b=1.78D0
* (dimensionless)
kb=DEXP(17.142D0)
* (number of particles/cm^3/minutes)
* (the units have been corrected from that reported in
* Table 3.1 in Miller)
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*, 'The results of the time domain simulation are:'
PRINT*, ' '
*Objectives
avg_r1(1) = moment10(1)/moment00(1)
avg_r2(1) = moment01(1)/moment00(1)
area_ratio(1) = 2*moment02(1)/moment20(1)
cv_r1(1) = DSQRT(moment20(1)*moment00(1)/
& (moment10(1))**2.0D0-1.0D0)
cv_r2(1) = DSQRT(moment02(1)*moment00(1)/
& (moment01(1))**2.0D0-1.0D0)
mass_ratio(1)=(moment21(1)-seed_moment21(1))/
& seed_moment21(1)
avg_vol(1) = moment21(1)/moment00(1)
tot_surf(1) = 4*moment02(1)+2*moment20(1)
tot_mass(1) = alpha*moment21(1)
wt_size_r1(1) = moment31(1)/moment21(1)
wt_size_r2(1) = moment22(1)/moment21(1)
*Simulation starts
time(1)=0.0D0
T=time(1)
DO 1200 I=1,(NN-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)
CALL lsodes ( FCN,NEQ,Y,T,T+delt,itol,rtol,
& atol,itask,istate,iopt,
& rwork,lrw,iwork,liw,FCNJ, mf )
moment00(I+1)=Y(1)
moment10(I+1)=Y(2)
moment01(I+1)=Y(3)
moment11(I+1)=Y(4)
moment20(I+1)=Y(5)
moment02(I+1)=Y(6)
moment21(I+1)=Y(7)
moment12(I+1)=Y(8)
moment22(I+1)=Y(9)
moment30(I+1)=Y(10)
moment31(I+1)=Y(11)
concentration(I+1)=Y(12)
* concentration_measured(I+1)=concentration(I+1)+
* & betaconc*DRNNOF()
seed_moment10(I+1)=Y(13)
seed_moment01(I+1)=Y(14)⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -