📄 dim2.f
字号:
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()
* Objectives
avg_r1(I+1) = moment10(I+1)/moment00(I+1)
avg_r2(I+1) = moment01(I+1)/moment00(I+1)
area_ratio(I+1) = 2*moment02(I+1)/moment20(I+1)
cv_r1(I+1) = DSQRT(moment20(I+1)*moment00(I+1)/
& (moment10(I+1))**2.0D0-1.0D0)
cv_r2(I+1) = DSQRT(moment02(I+1)*moment00(I+1)/
& (moment01(I+1))**2.0D0-1.0D0)
mass_ratio(I+1)=(moment21(I+1)-seed_moment21(I+1))/
& seed_moment21(I+1)
avg_vol(I+1) = moment21(I+1)/moment00(I+1)
tot_surf(I+1) = 4*moment02(I+1)+2*moment20(I+1)
tot_mass(I+1) = alpha*moment21(I+1)
wt_size_r1(I+1) = moment31(I+1)/moment21(I+1)
wt_size_r2(I+1) = moment22(I+1)/moment21(I+1)
time(I+1)=T+delt
1200 CONTINUE
*Print final results
****************************************************************************
PRINT*,'The weight mean size (in microns) of r1 is ',
&wt_size_r1(NN)
*
PRINT*,'The weight mean size (in microns) of r2 is ',
& wt_size_r2(NN)
PRINT*, 'The average r1 is ', avg_r1(NN)
PRINT*, 'The average r2 is ', avg_r2(NN)
PRINT*, 'The area ratio is ', area_ratio(NN)
PRINT*,'The coefficient of variation of r1 is ', cv_r1(NN)
PRINT*,'The coefficient of variation of r2 is ', cv_r2(NN)
PRINT*,'The nucleation to seed mass ratio is ',mass_ratio(NN)
PRINT*, 'The total surface area is ', tot_surf(NN)
PRINT*, 'The total mass is ', tot_mass(NN)
PRINT*, 'The average volume is ', avg_vol(NN)
DO 1300 I=1,NN
WRITE(30,1800) time(I)
WRITE(25,1900) temperature(I),concentration(I), relsatn(I)
WRITE(10,1500) moment00(I), moment10(I),
& moment01(I), moment11(I)
WRITE(11,1500) moment20(I), moment02(I),
& moment21(I), moment12(I)
WRITE(12,1900) moment22(I), moment30(I), moment31(I)
WRITE(20,1600) time(I), seed_moment10(I),
& seed_moment01(I),seed_moment11(I),
& seed_moment20(I),seed_moment02(I),
& seed_moment21(I), seed_moment12(I)
WRITE(15,1700)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)
1300 CONTINUE
1500 FORMAT(4(E13.6, 2X))
1600 FORMAT(11(E13.6, 2X))
1700 FORMAT(12(E13.6, 2X))
1800 FORMAT(E13.6)
1900 FORMAT(3(E13.6))
*the following commented commands only work with exponent graphics
* call af_set_titles('Temperature (degree C) vs Time (min)',
* & 'Time', 'Temperature')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_temperature,1)
* call af_set_titles('Concentration (g/g solvent) vs Time (min)',
* & 'Time', 'Concentration')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_concentration,1)
* call af_set_titles('Transmittance vs Time (min)',
* & 'Time', 'Transmittance')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_transmittance,1)
* call af_set_titles('Relative Supersaturation vs Time (min)',
* & 'Time', 'Rel. Supersaturation')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_satn,1)
* call af_set_titles('Zeroth moment (#/g solvent) vs Time (min)',
* & 'Time', 'Zeroth Moment')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_moment0,1)
* call af_set_titles('First Moment (cm/g solvent) vs Time (min)',
* & 'Time', 'First Moment')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_moment1,1)
* call af_set_titles('Second Moment (cm^2/g solvent) vs Time (min)',
* & 'Time', 'Second Moment')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_moment2,1)
* call af_set_titles('Third Moment (cm^3/g solvent) vs Time (min)',
* & 'Time', 'Third Moment')
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & NN,real_time, real_moment2, 1)
* CALL AF_SET('.1.axis.1.data.2 number x_copy y_copy flag$',
* & NN,real_time, real_moment3, 1)
* iunit = 2
C error_flag = af_scatter_plot(timax,realtime,realcout,0)
C call af_set_titles('Outlet Concentration','time','Concentration')
C call af_plot(iunit, 0)
* CALL AF_SET('.1.axis.1.data.1 number x_copy y_copy flag$',
* & timax, realtime, realcout, 1)
* CALL AF_SET('.1.axis.1.data.2 number x_copy y_copy flag$',
* & timax, realtime, realp, 1)
* CALL AF_SET('.2.axis.1.data.1 number x_copy y_copy flag$',
* & timax, realtime, residuals, 1)
* CALL AF_MULTIPLE_PLOT(2, iUNIT, 0)
* IUNIT=12
* CALL AF_PLOT(IUNIT,0)
STOP
END
**************************************************************************
SUBROUTINE FCN(NEQ,T,Y,Yprime)
INTEGER NEQ
REAL*8 T, Y(NEQ), Yprime(NEQ)
* REAL*8 kb, b, kg1, g1, kg2, g2
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 FCNJ(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 flag
REAL*8 time
COMMON /DATA1/flag
IF(flag.EQ.1) THEN
* constant-rate cooling
Temp=-(32D0-28D0)/160D0*time+32D0
ELSEIF (flag.EQ.2) THEN
* natural cooling
Temp=(28.0D0-32.0D0)*(1.0D0-DEXP(-time/310.0D0))+32.0D0
ELSE
* user-supplied
Temp=-(32D0-28D0)/160D0*time+32D0
ENDIF
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
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.9971-1.06*conc/(1+conc)+1.007*(conc/(1+conc))**2
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.75+92.82*conc-99.97*conc**2
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_DATA/kg1, g1
if ((conc-concs) .gt. 0) then
growth1 = 1.0D0*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 + -