dim1.f

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

 1200 CONTINUE
	DO I=1,(NN-1)

      WRITE(31,*) concentration(I)

	ENDDO

	DO I=1,(NN-1),30

	WRITE(32,*) moment0(I), moment1(I), moment2(I),

     &            moment3(I), moment4(I)

	ENDDO

	CLOSE(UNIT=31)

	CLOSE(UNIT=32)


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 variation is ', cov
      PRINT*,'The nucleation to seed mass ratio is ', mass_ratio
	PRINT*,'  '


      OPEN(UNIT=20, FILE='moments.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')

      OPEN(UNIT=25, FILE='seed.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')

      OPEN(UNIT=30, FILE='cond.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')

      OPEN(UNIT=35, FILE='obj.dat', FORM='FORMATTED',
     &     ACCESS='SEQUENTIAL', STATUS='UNKNOWN')

      DO I=1,NN-1
	 WRITE(20,1400)time(I),moment0(I), moment1(I), 
     & 		       moment2(I), moment3(I), moment4(I)
	 WRITE(25,1500)time(I), seed_moment1(I), 
     & 		       seed_moment2(I), seed_moment3(I)
	 WRITE(30,1400)time(I), concentration(I), 
     & 	 temperature(I), relsatn(I), transmittance(I) 
       ENDDO
       WRITE(35,1600)weight_mean_size, cov, mass_ratio 

 1400  FORMAT(6(E25.16,1X))
 1500  FORMAT(4(E25.16,1X))
 1600  FORMAT(3(E25.16, 1X))


*the following section should be commented for PCs
*      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
*      error_flag = af_scatter_plot(timax,realtime,realcout,0)
*      call af_set_titles('Outlet Concentration','time','Concentration')
*      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, 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,*)


*the following commands exist only to remove compiler warnings

	T = T

	IF (Y(1).LT.0) THEN

	ENDIF

	IF (DYPDY(1,1).LT.0) THEN

	ENDIF


      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=-(32.0D0-28.0D0)/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=-(32.0D0-28.0D0)/160.0D0*time+32.0D0
      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.9971D0-1.06D0*conc/(1.0D0+conc)

     &     + 1.007D0*(conc/(1.0D0+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.75D0+92.82D0*conc-99.97D0*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
	if ((conc-concs) .gt. 0) then

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

	   elseif ((conc-concs) .eq. 0) then

	   growth = 0.0D0

	   else

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

	  endif


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

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

*     birth 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
	  

	  if ((conc-concs) .gt. 0) then

	     birth = kb*(((conc-concs)/concs)**b)*m3*(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)

	&	         *m3*(1.0D-4)**3.0D0

	  endif


      RETURN
      END