dim1sun.f

对工业生产过程结晶过程的一个仿真程序软件包

         Y(7)=seed_moment1(I)
         Y(8)=seed_moment2(I)
         Y(9)=seed_moment3(I)
* Solving the moment equations
        CALL lsodes ( FCN,NEQ,Y,T,T+1,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)
*         CALL RNUN(1,noise_conc)
         concentration_measured(I+1)=concentration(I+1)+
     &             betaconc*noise_trans(I)
         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))
*         CALL RNUN(1,noise_trans)
         transmittance(I+1)=DEXP(-ka/2.0D0*cell_length/10.0D0*
     &        moment2(I+1)*(densitys*(1.0D-4)**2.D0))+
     &	      betatrans*noise_conc(I)
         time(I+1)=T
*	 WRITE(98,*) concentration_measured(I+1)
*         WRITE(99,*) transmittance(I+1)
	variance_trans=variance_trans+
     &		(betatrans*noise_trans(I))**2.0D0
        variance_conc=variance_conc+
     &		(betaconc*noise_conc(I))**2.0D0

	WRITE(58,1) concentration_measured(I+1), transmittance(I+1)

 1200 CONTINUE
      print*,'variance_trans = ',variance_trans/160.0D0
      print*,'variance_conc = ', variance_conc/160.0D0
      print*,'sigma_trans = ',DSQRT(variance_trans/160.0D0)
      print*,'sigma_conc = ', DSQRT(variance_conc/160.0D0)


1      FORMAT(2(F13.6,3x))
2      FORMAT(I1)
  
      weight_mean_size=moment4(NN)/moment3(NN)
      cov=DSQRT(moment2(NN)*moment0(NN)/(moment1(NN))**2.0D0-1.0D0)
      mass_ratio=(moment3(NN)-seed_moment3(NN))/
     &     seed_moment3(NN)

      PRINT*,'The weight mean size (in microns) is ',
     &     weight_mean_size
      PRINT*,'The coefficient of variance is ', cov
      PRINT*,'The nucleation to seed mass ratio is ', mass_ratio

      DO 1300 I=1,NN
	real_time(I)=sngl(time(I))
	real_temperature(I)=sngl(Temp(time(I)))
	real_concentration(I)=sngl(concentration_measured(I))
	real_moment0(I)=sngl(moment0(I))
	real_moment1(I)=sngl(moment1(I))*1E-4
	real_moment2(I)=sngl(moment2(I))*(1E-4)**2
	real_moment3(I)=sngl(moment3(I))*(1E-4)**3
	real_transmittance(I)=sngl(transmittance(I))
	real_satn(I)=sngl(relsatn(I))
1300  CONTINUE

*      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=2
      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, kg, g
      REAL*8 r0, alpha, mu00, UA, Msolv
      REAL*8 Csat, birth, growth, Temp
*      COMMON /GROWTH_DATA/kg, g
*      COMMON /BIRTH_DATA/kb, b
      COMMON /EXP_DATA/r0, alpha, mu00, UA, Msolv

*Moments 
*     d(mu0(t))/dt, # of particles/g solvent/minute
      Yprime(1)=birth(Y(6),Csat(Temp(T)),Y(4))
*     d(mu1(t))/dt, microns/g solvent/minute
      Yprime(2)=growth(Y(6),Csat(Temp(T)))*Y(1)+
     &     birth(Y(6),Csat(Temp(T)),Y(4))*r0
*     d(mu2(t))/dt, microns^2/g solvent/minute
      Yprime(3)=2*growth(Y(6),Csat(Temp(T)))*Y(2)+
     &     birth(Y(6),Csat(Temp(T)),Y(4))*r0**2
*     d(mu3(t))/dt, microns^3/g solvent/minute
      Yprime(4)=3*growth(Y(6),Csat(Temp(T)))*Y(3)+
     &     birth(Y(6),Csat(Temp(T)),Y(4))*r0**3
*     d(mu4(t))/dt, microns^4/g solvent/minute
      Yprime(5)=4*growth(Y(6),Csat(Temp(T)))*Y(4)+
     &     birth(Y(6),Csat(Temp(T)),Y(4))*r0**4
*Mass balance
*     d(conc(t))/dt, g solute/g solvent/minute
      Yprime(6)=-alpha*(3.0D0*growth(Y(6),Csat(Temp(T)))*Y(3)+
     &     birth(Y(6),Csat(Temp(T)),Y(4))*r0**3)
*Moments for seed crystals only
*     d(mu'1(t)/dt), microns/g solvent/minute
      Yprime(7)=growth(Y(6),Csat(Temp(T)))*mu00
*     d(mu'2(t)/dt), microns^2/g solvent/minute
      Yprime(8)=2*growth(Y(6),Csat(Temp(T)))*Y(7)
*     d(mu'2]3(t)/dt), microns^3/g solvent/minute
      Yprime(9)=3*growth(Y(6),Csat(Temp(T)))*Y(8)

      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

*       Temp=(28.0D0-32.0D0)*(1.0D0-DEXP(-time/310.0D0))+32.0D0

*      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
*      Temp=32.0D0-0.0625D0*time
 
      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 growth(conc, concs)

*     growth rate for the simulation of a cooling batch
*     crystallizer
*     
*     arguments:  conc - solute concentration
*                 concs - saturation concentration
*     non-argument input: kg, g kinetic rate parameters
*     output: growth - growth rate

      REAL*8 conc, concs, kg, g
      COMMON /GROWTH_DATA/kg, g

      growth=kg*((conc-concs)/concs)**g

      RETURN
      END
************************************************************

      REAL*8 FUNCTION birth(conc, concs,m3)

*     birthth rate for the simulation of a cooling batch
*     crystallizer
*
*     arguments:  conc - solute concentration
*                 concs - saturation concentration
*                 m3 - 3rd moment
*     non-argument input: kb, b kinetic rate parameters
*     output: birth - birth rate

      REAL*8 conc, concs, m3, kb, b
      COMMON /BIRTH_DATA/kb, b

      birth=kb*((conc-concs)/concs)**b*m3*(1.0D-4)**3

      RETURN
      END