📄 param2.f
字号:
mean_seed = 196.0D0
* Percent width of the initial seed distribution
width = 6.12D0
*Growth and nucleation kinetic parameters (Table 4.6 in Miller)
* (dimensionaless)
g1=THETA(1)
* (mirons/minute)
kg1=DEXP(THETA(2))
* (dimensionless)
g2=THETA(3)
* (microns/minute)
kg2=DEXP(THETA(4))
* (dimensionless)
b=THETA(5)
* (number of particles/cm^3/minute)
* (the units have been corrected from that reported in
* Table 3.1 in Miller)
kb=DEXP(THETA(6))
* Optimization parameters:
* sigma_trans=0.009D0
* sigma_conc=0.0005D0
* weight_trans=(1.0D0/sigma_trans)**2.0D0
* weight_conc=(1.0D0/sigma_conc)**2.0D0
*Initial conditions
* initial concentration, g/g solvent
concentration(1) = 0.493D0
* concentration_measured(1)=concentration(1)
*
*
* 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))
*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-8
atol=1.0D-8
itol=1
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 ( MOMENTS,NEQ,Y,T,T+1,itol,rtol,
& atol,itask,istate,iopt,
& rwork,lrw,iwork,liw,MOMENTSJ, 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)
seed_moment11(I+1)=Y(15)
seed_moment20(I+1)=Y(16)
seed_moment02(I+1)=Y(17)
seed_moment21(I+1)=Y(18)
seed_moment12(I+1)=Y(19)
relsatn(I+1)=(concentration(I+1)-Csat(Temp(T)))/
& Csat(Temp(T))
* transmittance(I+1)=DEXP(-ka/2.0D0*cell_length/10.0D0*
* & moment2(I+1)*(densitys*(1.0D-4)**2.D0))+betatrans*DRNNOF()
time(I+1)=T
1200 CONTINUE
K = 1
DO I = 1,(NN-1),Cinterval
Phi = Phi+
& weight_conc*(ConcData(Iteration,K)-
& concentration(I))**2.0D0
K = K+1
ENDDO
K=1
DO I=1,(NN-1),Minterval
Phi=Phi+weight_mu00*(Mu00Data(Iteration,K)-
& moment00(I))**2.0D0+
& weight_mu10*(Mu10Data(Iteration,K)-
& moment10(I))**2.0D0+
& weight_mu01*(Mu01Data(Iteration,K)-
& moment01(I))**2.0D0+
& weight_mu11*(Mu11Data(Iteration,K)-
& moment11(I))**2.0D0+
& weight_mu20*(Mu20Data(Iteration,K)-
& moment20(I))**2.0D0+
& weight_mu02*(Mu02Data(Iteration,K)-
& moment02(I))**2.0D0+
& weight_mu21*(Mu21Data(Iteration,K)-
& moment21(I))**2.0D0+
& weight_mu12*(Mu12Data(Iteration,K)-
& moment12(I))**2.0D0+
& weight_mu22*(Mu22Data(Iteration,K)-
& moment22(I))**2.0D0+
& weight_mu30*(Mu30Data(Iteration,K)-
& moment30(I))**2.0D0+
& weight_mu31*(Mu31Data(Iteration,K)-
& moment31(I))**2.0D0
K = K + 1
ENDDO
777 CONTINUE
RETURN
END
****************************************************************
SUBROUTINE MOMENTS(NEQ,T,Y,Yprime)
INTEGER NEQ
REAL*8 T, Y(NEQ), Yprime(NEQ)
REAL*8 r0, alpha, mu00, UA, Msolv
REAL*8 Csat, birth, growth1, growth2, Temp
COMMON /EXP_DATA/r0, alpha, mu00, UA, Msolv
*Moments
* 00 moment, number of particles/g solvent/minute
Yprime(1)=birth(Y(12),Csat(Temp(T)),Y(7))
* 10 moment, microns/g solvent/minute
Yprime(2)=growth1(Y(12),Csat(Temp(T)))*Y(1)
* 01 moment, microns/g solvent/minute
Yprime(3)=growth2(Y(12),Csat(Temp(T)))*Y(1)
* 11 moment, microns^2/g solvent/minute
Yprime(4)=growth1(Y(12),Csat(Temp(T)))*Y(3)
& + growth2(Y(12),Csat(Temp(T)))*Y(2)
* 20 moment, microns^2/g solvent/minute
Yprime(5)=2*growth1(Y(12),Csat(Temp(T)))*Y(2)
* 02 moment, microns^2/g solvent/minute
Yprime(6)=2*growth2(Y(12),Csat(Temp(T)))*Y(3)
* 21 moment, microns^3/g solvent/minute
Yprime(7)=2*growth1(Y(12),Csat(Temp(T)))*Y(4)+
& growth2(Y(12),Csat(Temp(T)))*Y(5)
* 12 moment, microns^3/g solvent/minute
Yprime(8)=growth1(Y(12),Csat(Temp(T)))*Y(6)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(4)
* 22 moment, microns^4/g solvent/minute
Yprime(9)=2*growth1(Y(12),Csat(Temp(T)))*Y(8)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(7)
* 30 moment, microns^4/g solvent/minute
Yprime(10)=3*growth1(Y(12),Csat(Temp(T)))*Y(5)
* 31 moment, microns^4/g solvent/minute
Yprime(11)=3*growth1(Y(12),Csat(Temp(T)))*Y(7)+
& growth2(Y(12),Csat(Temp(T)))*Y(10)
*Mass balance
* d(conc(t))/dt, g solute/g solvent/minute
Yprime(12)=-alpha*2.0D0*growth1( Y(12), Csat(Temp(T)) )
& *Y(4) - alpha*growth2( Y(12),Csat(Temp(T) ) )*Y(5)
*Moments for seed crystals only
* 10 moment, microns/g solvent/minute
Yprime(13)=growth1(Y(12),Csat(Temp(T)))*mu00
* 01 moment, microns/g solvent/minute
Yprime(14)=growth2(Y(12),Csat(Temp(T)))*mu00
* 11 moment, microns^2/g solvent/minute
Yprime(15)=growth1(Y(12),Csat(Temp(T)))*Y(14)
& + growth2(Y(12),Csat(Temp(T)))*Y(13)
* 20 moment, microns^2/g solvent/minute
Yprime(16)=2*growth1(Y(12),Csat(Temp(T)))*Y(13)
* 02 moment, microns^2/g solvent/minute
Yprime(17)=2*growth2(Y(12),Csat(Temp(T)))*Y(14)
* 21 moment, microns^3/g solvent/minute
Yprime(18)=2*growth1(Y(12),Csat(Temp(T)))*Y(15)+
& growth2(Y(12),Csat(Temp(T)))*Y(16)
* 12 moment, microns^3/g solvent/minute
Yprime(19)=growth1(Y(12),Csat(Temp(T)))*Y(17)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(15)
RETURN
END
****************************************************************
SUBROUTINE MOMENTSJ(NEQ,T,Y,DYPDY)
INTEGER NEQ
REAL*8 T,Y(NEQ),DYPDY(NEQ,*)
RETURN
END
****************************************************************
REAL*8 FUNCTION Temp(time)
* Crystallizer temperature setpoint profile for the
* simulation of a batch cooling crystallizer
*
* input: time - minutes
* output: Temp - temperature in degrees Centigrade
* INTEGER I, J
INTEGER Ntemp3,Iteration
PARAMETER(Ntemp3 = 8)
REAL*8 time
* REAL*8 interval, sum
* REAL*8 Coeff(Ntemp3)
COMMON Iteration
* constant-rate cooling
* Temp=-(32D0-28D0)/160D0*time+32D0
Temp=(28.0D0-32.0D0)*(1.0D0-DEXP(-time/310.0D0))+32.0D0
* IF(Iteration.EQ.1) THEN
* Coeff(1) = -1.0000000000000D-01
* Coeff(2) = 5.0835793601210D-35
* Coeff(3) = -4.8624953603150D-02
* Coeff(4) = -1.0000000000000D-01
* Coeff(5) = -1.0000000000000D-01
* Coeff(6) = -1.0000000000000D-01
* Coeff(7) = -5.1375046396850D-02
* Coeff(8) = -5.2673237458465D-34
* ENDIF
* IF(Iteration.EQ.2) THEN
* Coeff(1) = -1.0000000000000D-01
* Coeff(2) = 0.0D0
* Coeff(3) = -1.0000000000000D-01
* Coeff(4) = -1.0000000000000D-01
* Coeff(5) = -1.0000000000000D-01
* Coeff(6) = -1.0000000000000D-01
* Coeff(7) = 0.0D0
* Coeff(8) = -2.7489751596006D-32
* ENDIF
* interval = DFLOAT(160/Ntemp3)
* sum=0.0D0
* DO 2666 I = 1, (Ntemp3-1)
* IF(time.LE.(DFLOAT(I)*interval)) THEN
* DO 2667 J = 1, (I-1)
* sum=sum+Coeff(J)
*2667 CONTINUE
* Temp=32.0D0+sum*interval+
* & Coeff(I)*(time-DFLOAT(I-1)*interval)
* RETURN
* ENDIF
*2666 CONTINUE
* DO 2668 I = 1, (Ntemp3-1)
* sum=sum+Coeff(I)
*2668 CONTINUE
* Temp=32.0D0+sum*interval+Coeff(Ntemp3)*
* & (time-DFLOAT(Ntemp3-1)*interval)
*
RETURN
END
*************************************************************
REAL*8 FUNCTION Csat(T)
* saturation concentration for the simulation of a
* cooling batch crystallizer (potassium nitrate-water)
* system, from Appendix C in Miller
*
* input: T - temperature (20-40 degree Centigrade)
* output: Csat - saturation concentration
* (g KNO3/g water)
REAL*8 T
Csat=0.1286D0+0.00588D0*T+0.0001721D0*T**2.0D0
RETURN
END
****************************************************************
REAL*8 FUNCTION heat(conc)
* heat capacity for the simulation of a cooling
* batch crystallizer (potassium nitrate-water
* system), from Appendix C in Miller
*
* input: conc - soluction concentration
* output: heat - heat capacity, calories/g solvent/ degree C
REAL*8 conc
heat=0.9971D0-1.06D0*conc/(1.0D0+conc)+1.007D0*(conc/(1.0D0
& +conc))**2.0D0
RETURN
END
*****************************************************************
REAL*8 FUNCTION enthalpy(conc)
* enthalpy of crystallization for the simulation of a
* cooling batch crystallizer (potassium nitrate-water
* system), from Appendix C in Miller
*
* input: conc - solute concentration
* output: enthalpy - enthalpy of crystallization
* (calories/g solvent)
REAL*8 conc
enthalpy=-85.75D0+92.82D0*conc-99.97D0*conc**2.0D0
RETURN
END
************************************************************
REAL*8 FUNCTION growth1(conc, concs)
* growth rate for the simulation of a cooling batch
* crystallizer
*
* arguments: conc - solute concentration
* concs - saturation concentration
* non-argument input: kg1, g1 kinetic rate parameters
* output: growth - growth rate
REAL*8 conc, concs, kg1, g1
COMMON /GROWTH_DATA1/kg1, g1
if ((conc-concs) .gt. 0) then
growth1=kg1*((conc-concs)/concs)**g1
elseif ((conc-concs) .eq. 0) then
growth1 = 0.0D0
else
growth1 = -1.0D0*kg1*(-1.0D0*(conc-concs)/concs)**g1
endif
RETURN
END
************************************************************
REAL*8 FUNCTION growth2(conc, concs)
* growth rate for the simulation of a cooling batch
* crystallizer
*
* arguments: conc - solute concentration
* concs - saturation concentration
* non-argument input: kg2, g2 kinetic rate parameters
* output: growth - growth rate
REAL*8 conc, concs, kg2, g2
COMMON /GROWTH_DATA2/kg2, g2
if ((conc-concs) .gt. 0) then
growth2 = 5.0D0*kg2*((conc-concs)/concs)**g2
elseif ((conc-concs) .eq. 0) then
growth2 = 0.0D0
else
growth2 = -5.0D0*kg2*(-1.0D0*(conc-concs)/concs)**g2
endif
RETURN
END
************************************************************
REAL*8 FUNCTION birth(conc, concs,m21)
* birth rate for the simulation of a cooling batch
* crystallizer
*
* arguments: conc - solute concentration
* concs - saturation concentration
* m21 - moment21
* non-argument input: kb, b kinetic rate parameters
* output: birth - birth rate
REAL*8 conc, concs, m21, kb, b
COMMON /BIRTH_DATA/kb, b
if ((conc-concs) .gt. 0) then
birth = kb*(((conc-concs)/concs)**b)*m21*(1.0D-4)**3.0D0
elseif ((conc-concs) .eq. 0) then
birth = 0.0D0
else
birth = -1.0D0*kb*((-1.0D0*(conc-concs)/concs)**b)
& *m21*(1.0D-4)**3.0D0
endif
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -