optconsun.f

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

*      print*,gamma(1), gamma(2), gamma(3)



*     initial zeroth moment, number of particle/g solvent
      mu00 = gamma(1)*(lmax**3-lmin**3)/3.0D0+
     &	     gamma(2)*(lmax**2-lmin**2)/2.0D0+
     &	     gamma(3)*(lmax-lmin)
      moment0(1)= mu00
*     initial first moment, micron/g solvent
      moment1(1) = gamma(1)*(lmax**4-lmin**4)/4.0D0+
     &		   gamma(2)*(lmax**3-lmin**3)/3.0D0+
     &	           gamma(3)*(lmax**2-lmin**2)/2.0D0
      seed_moment1(1)=moment1(1)
*     initia1 second moment, micron^2/g solvent
      moment2(1) = gamma(1)*(lmax**5-lmin**5)/5.0D0+
     &		   gamma(2)*(lmax**4-lmin**4)/4.0D0+
     &	           gamma(3)*(lmax**3-lmin**3)/3.0D0
      seed_moment2(1)=moment2(1)
*     initial third moment, micron^3/g solvent
      moment3(1) = gamma(1)*(lmax**6-lmin**6)/6.0D0+
     &		   gamma(2)*(lmax**5-lmin**5)/5.0D0+
     &	           gamma(3)*(lmax**4-lmin**4)/4.0D0
      seed_moment3(1)=moment3(1)
*	print*,moment3(1)*Msolv*densityc*(1.0e-12)
*	print*,THETA_T(Ntemp1+1)
*     initial fourth moment, micron^4/g solvent
      moment4(1) = gamma(1)*(lmax**7-lmin**7)/7.0D0+
     &		   gamma(2)*(lmax**6-lmin**6)/6.0D0+
     &	           gamma(3)*(lmax**5-lmin**5)/5.0D0

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

*Growth and nucleation kinetic parameters (Table 4.6 in Miller)
*     (dimensionaless)
      g=1.32D0
      g = 0.13136920450278D+01
*     (mirons/minute)
      kg=DEXP(8.849D0)
      kg = DEXP(0.87898608553756D+01)
*     (dimensionless)
      b=1.78D0
      b = 0.18393707105317D+01   
*     (number of particles/cm^3/minute) 
*     (the units have been corrected from that reported in
*     Table 3.1 in Miller)
      kb=DEXP(17.142D0)
      kb = DEXP(0.17382782278677D+02)

*Simulation parameters
*********************************************************
*
      mf=222
      itask=1
      istate =1
      iopt=0
      lrw=3800
      liw=200
      rtol=1.0d-10
      atol=1.0d-9
      itol=1


      time(1)=0.0D0
      T=0.0D0

* Start the loop to solve the moment equations
      DO 1200 I=1,(NN-1)
	 temperature(I)=Temp(T)
         Y(1)=moment0(I)
         Y(2)=moment1(I)
         Y(3)=moment2(I)
         Y(4)=moment3(I)
         Y(5)=moment4(I)
         Y(6)=concentration(I)
         Y(7)=seed_moment1(I)
         Y(8)=seed_moment2(I)
         Y(9)=seed_moment3(I)
         
         CALL lsodes ( MOMENTS,NEQ,Y,T,T+1,itol,rtol,
     &          atol,itask,istate,iopt, 
     &           rwork,lrw,iwork,liw,MOMENTSJ, 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)
         concentration_measured(I+1)=concentration(I+1)
         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))
         transmittance(I+1)=DEXP(-ka/2.0D0*cell_length/10.0D0*
     &        moment2(I+1)*(densitys*(1.0D-4)**2.D0))
         time(I+1)=T

 1200 CONTINUE
  

* weigth_mean_size, Eqn 2.6-5 in Randolph and Larson
      weight_mean_size=moment4(NN)/moment3(NN)
* coefficient of variation for population-size distribution
* Eqn 2.6-6 in Randolph and Larson
      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

* To maximized the weight mean size
        function=-weight_mean_size
* To minimize the coefficient of variance
*        function=cov
* To minimize the ratio of the nucleated mass to seed mass
*        function=mass_ratio

      RETURN
      END

*****************************************************************

      SUBROUTINE cntr(Nvar,jjj,x,gj)
*
*	cntr is the subroutine for the constraints.  The
*	constraints are listed in the following order:   
*	Nonlinear inequality constraints, linear inequality 
*	constraints, nonlinear equality constraints, and linear 
*	equality constraints.
*
*     input:   Nvar - number of parameters
*	       jjj  - indicates the jjj_th constraint
*	       x    - Nvar-dimensional vector of parameters
*	       
*     output:  gj   - the jjj_th constraint
*

      INTEGER Nvar,jjj, Ntemp2, I
      PARAMETER (Ntemp2 = 8)
      REAL*8 x(*),gj, Coeff(Ntemp2+3), umin, umax, Temp, Cf_max
      REAL*8 concentration(161), function
      COMMON Coeff
      COMMON concentration

      DO 628 I = 1, Nvar
         Coeff(I)=x(I)
628   CONTINUE

*     Maximum temperature allowed
      umax=32.3D0
*     Minimum temperature allowed
      umin=22.0D0

* Cf_max is the maximum final time concentration that is required to
* satisfy the minimum yield constraint.  It is set equal to the
* the final time concentration for the linear profile temperature
* temperature = 32.0-0.0625*time (the linear profile with the
* steepest descend without violating the minimum temperature
* constraint).
         CALL FCN(Nvar,jjj,x,function)
*	print*,concentration(161)

      Cf_max=0.342D0
      IF(jjj.EQ.1) THEN
         CALL FCN(Nvar,jjj,x,function)
*	print*,concentration(161)
         gj=concentration(161)-Cf_max
      ELSE 
         gj=umin-Temp(160.0D0)
      ENDIF

*	IF (jjj.EQ.9) print*,gj

      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, Ntemp3
      PARAMETER(Ntemp3 = 8)
      REAL*8 time, interval, sum
      REAL*8 Coeff(Ntemp3+3)
      COMMON Coeff

      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
      
****************************************************************
      SUBROUTINE MOMENTS(NEQ,T,Y,Yprime)
*
*	Subroutine MOMENTS lists the moment equations and
*	the mass balance equations.  The subroutine
*	is called by lsode subruoution, which
*	solves the ODE's
*
      INTEGER NEQ
      REAL*8 T, Y(NEQ), Yprime(NEQ)
      REAL*8 r0, alpha, mu00, UA, Msolv
      REAL*8 Csat, birth, growth, Temp
      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 MOMENTSJ(NEQ,T,Y,DYPDY)
      INTEGER NEQ
      REAL*8 T,Y(NEQ),DYPDY(NEQ,*)


      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 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.0D0

      RETURN
      END