paramsun.f

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

*	Then
*
*	(lmax^6-lmin^6)/6 a + (lmax^5-lmin^5)/5 b + (lmax^4-lmin^4)/4 c =
*	mass/(mass_solvent*crystal_density)
*
*	lmax^2 a + lmax b + c = 0
*	lmin^2 a + lmax b + c = 0
*
*	Note: In the implementation, the first equation is scaled.
*

      lmin = mean_seed-width*(1.0D-2)*mean_seed
      lmax = mean_seed+width*(1.0D-2)*mean_seed

      AA(1,1) = (lmax**6-lmin**6)/6.0D0/1.0D12 
      AA(1,2) = (lmax**5-lmin**5)/5.0D0/1.0D12
      AA(1,3) = (lmax**4-lmin**4)/4.0D0/1.0D12
      AA(2,1) = lmin**2
      AA(2,2) = lmin
      AA(2,3) = 1.0D0
      AA(3,1) = lmax**2
      AA(3,2) = lmax
      AA(3,3) = 1.0D0
      BB(1)=mass_seed/
     &		(Msolv*densityc)*(1.0D4)**3/1.0D12
      BB(2)=0.0D0
      BB(3)=0.0D0

      CALL DLSARG (Norder, AA, LDA, BB, IPATH, gamma)

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



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


      time(1)=0.0D0
      T=0.0D0
      
      DO 1201 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)
*	Write(888,*)Coeff(1),Coeff(2)

         
         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)
         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(TEND)))/
     &        Csat(Temp(TEND))
         transmittance(i+1)=DEXP(-ka/2D0*cell_length/10D0*
     &        moment2(i+1)*(densitys*(1D-4)**2))
         time(I+1)=T
 1201 CONTINUE

      DO 1251 I = 1, (NN-1)
        Phi = Phi+
     &	       weight_conc*(ConcData(Iteration,I)-
     &		concentration(I+1))**2.0D0+
     &	       weight_trans*(TransData(Iteration,I)-
     &		transmittance(I+1))**2.0D0
	error_conc=error_conc+
     &	   (ConcData(Iteration, I)-concentration(I+1))**2.0D0
	error_trans=error_trans+
     &     (TransData(Iteration,I)-transmittance(I+1))**2.0D0
1251  CONTINUE
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, 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*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 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.0D2

      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