optexp2.f

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

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

      DO 203 I = 1 , NU
	DO 204 J = 1 , NU
	   JACOBIAN(I,J)=0.0D0
204     CONTINUE
203   CONTINUE


*The following equations assume r0 = 0 (nucleation crystal size)
      JACOBIAN(1,7)=birth(Y(12),Csat(Temp(T)),Y(7))/Y(7)
      JACOBIAN(1,12)=b*birth2(Y(12),Csat(Temp(T)),Y(7))/
     &			(Y(12)-Csat(Temp(T)))
      JACOBIAN(2,1)=growth1(Y(12),Csat(Temp(T)))
      JACOBIAN(2,12)=g1*growth3(Y(12),Csat(Temp(T)))*Y(1)/
     &			(Y(12)-Csat(Temp(T)))
      JACOBIAN(3,1)=growth2(Y(12),Csat(Temp(T)))
      JACOBIAN(3,12)=g2*growth4(Y(12),Csat(Temp(T)))*Y(1)/
     &			(Y(12)-Csat(Temp(T)))
      JACOBIAN(4,2)=growth2(Y(12),Csat(Temp(T)))
      JACOBIAN(4,3)=growth1(Y(12),Csat(Temp(T)))
      JACOBIAN(4,12)=g1*growth3(Y(12),Csat(Temp(T)))*Y(3)/
     &			(Y(12)-Csat(Temp(T)))+
	&			g2*growth4(Y(12),Csat(Temp(T)))*Y(2)/
     &			(Y(12)-Csat(Temp(T)))
	JACOBIAN(5,2)=2.0D0*growth1(Y(12),Csat(Temp(T)))
      
	JACOBIAN(5,12)=2.0D0*g1*growth3(Y(12),Csat(Temp(T)))*Y(2)/
     &			(Y(12)-Csat(Temp(T)))
	JACOBIAN(6,3)=2.0D0*growth2(Y(12),Csat(Temp(T)))
      
	JACOBIAN(6,12)=2.0D0*g2*growth4(Y(6),Csat(Temp(T)))*Y(3)
     &			/(Y(12)-Csat(Temp(T)))
	JACOBIAN(7,4)=2.0D0*growth1(Y(12),Csat(Temp(T)))
	JACOBIAN(7,5)=growth2(Y(12),Csat(Temp(T)))
	JACOBIAN(7,12)=2.0D0*g1*growth3(Y(12),Csat(Temp(T)))*Y(4)/
     &			(Y(12)-Csat(Temp(T)))+
	&	        g2*growth4(Y(12),Csat(Temp(T)))*Y(5)/
     &			(Y(12)-Csat(Temp(T)))
	JACOBIAN(8,4)=2.0D0*growth2(Y(12),Csat(Temp(T)))
	JACOBIAN(8,6)=growth1(Y(12),Csat(Temp(T)))
	JACOBIAN(8,12)=g1*growth3(Y(12),Csat(Temp(T)))*Y(6)/
     &			(Y(12)-Csat(Temp(T)))+
	&	        2.0D0*g2*growth4(Y(12),Csat(Temp(T)))*Y(4)/
     &			(Y(12)-Csat(Temp(T)))
	JACOBIAN(9,7)=2.0D0*growth2(Y(12),Csat(Temp(T)))
	JACOBIAN(9,8)=2.0D0*growth1(Y(12),Csat(Temp(T)))
	JACOBIAN(9,12)=2.0D0*g1*growth3(Y(12),Csat(Temp(T)))*Y(8)/
     &			(Y(12)-Csat(Temp(T)))+
	&	        2.0D0*g2*growth4(Y(12),Csat(Temp(T)))*Y(7)/
     &			(Y(12)-Csat(Temp(T)))
	JACOBIAN(10,5)=3.0D0*growth1(Y(12),Csat(Temp(T)))
	JACOBIAN(10,12)=3.0D0*g1*growth3(Y(12),Csat(Temp(T)))*Y(5)/
     &			(Y(12)-Csat(Temp(T)))
	JACOBIAN(11,7)=3.0D0*growth1(Y(12),Csat(Temp(T)))
	JACOBIAN(11,10)=5.0D0*growth2(Y(12),Csat(Temp(T)))
	JACOBIAN(11,12)=3.0D0*g1*growth3(Y(12),Csat(Temp(T)))*Y(7)/
     &			(Y(12)-Csat(Temp(T)))+
	&	        g2*growth4(Y(12),Csat(Temp(T)))*Y(10)/
     &			(Y(12)-Csat(Temp(T)))
	JACOBIAN(12,4)=-2.0D0*alpha*growth1(Y(12),Csat(Temp(T)))
	JACOBIAN(12,5)=-1.0D0*alpha*growth2(Y(12),Csat(Temp(T)))
	JACOBIAN(12,12)=-2.0D0*alpha*g1*growth3(Y(12),Csat(Temp(T)))*
	&          Y(4)/(Y(12)-Csat(Temp(T)))-
	&	        alpha*g2*growth4(Y(12),Csat(Temp(T)))*
	&          Y(5)/(Y(12)-Csat(Temp(T)))

      K=20
      DO 213 I =1 , NU
	DO 214 J = 1, NTHETA
	    W(I,J)=Y(K)
            K=K+1
214	CONTINUE
213   CONTINUE

      DO 223 I = 1 , NU
         DO 224 J = 1 , NTHETA
            F_THETA(I,J)=0.0D0
224      CONTINUE
223   CONTINUE

*Again, the following equations assume r0 = 0
      F_THETA(1,5)=birth(Y(12),Csat(Temp(T)),Y(7))*
     &		DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(1,6)=birth(Y(12),Csat(Temp(T)),Y(7))
      F_THETA(2,1)=Y(1)*growth1(Y(12),Csat(Temp(T)))*
     &		DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(2,2)=Y(1)*growth1(Y(12),Csat(Temp(T)))/kg1
      F_THETA(3,3)=Y(1)*growth2(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(3,4)=Y(1)*growth2(Y(12),Csat(Temp(T)))/kg2
      F_THETA(4,1)=Y(3)*growth1(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(4,2)=Y(3)*growth1(Y(12),Csat(Temp(T)))/kg1 
     	F_THETA(4,3)=Y(2)*growth2(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))	
      F_THETA(4,4)=Y(2)*growth2(Y(12),Csat(Temp(T)))/kg2
      F_THETA(5,1)=2.0D0*Y(2)*growth1(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(5,2)=2.0D0*Y(2)*growth1(Y(12),Csat(Temp(T)))/kg1	 
	F_THETA(6,3)=2.0D0*Y(3)*growth2(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))   
	F_THETA(6,4)=2.0D0*Y(3)*growth2(Y(12),Csat(Temp(T)))/kg2
	F_THETA(7,1)=2.0D0*Y(4)*growth1(Y(12),Csat(Temp(T)))*
	&          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(7,2)=2.0D0*Y(4)*growth1(Y(12),Csat(Temp(T)))/kg1
	F_THETA(7,3)=Y(5)*growth2(Y(12),Csat(Temp(T)))*
	&          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(7,4)=Y(5)*growth2(Y(12),Csat(Temp(T)))/kg2     
	F_THETA(8,1)=Y(6)*growth1(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
     	F_THETA(8,2)=Y(6)*growth1(Y(12),Csat(Temp(T)))/kg1
      F_THETA(8,3)=2.0D0*Y(4)*growth2(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(8,4)=2.0D0*Y(4)*growth2(Y(12),Csat(Temp(T)))/kg2
	F_THETA(9,1)=2.0D0*Y(8)*growth1(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
     	F_THETA(9,2)=2.0D0*Y(8)*growth1(Y(12),Csat(Temp(T)))/kg1
      F_THETA(9,3)=2.0D0*Y(7)*growth2(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(9,4)=2.0D0*Y(7)*growth2(Y(12),Csat(Temp(T)))/kg2
	F_THETA(10,1)=3.0D0*Y(5)*growth1(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(10,2)=3.0D0*Y(5)*growth1(Y(12),Csat(Temp(T)))/kg1
	F_THETA(11,1)=3.0D0*Y(7)*growth1(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
     	F_THETA(11,2)=3.0D0*Y(7)*growth1(Y(12),Csat(Temp(T)))/kg1
      F_THETA(11,3)=Y(10)*growth2(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(11,4)=Y(10)*growth2(Y(12),Csat(Temp(T)))/kg2
	F_THETA(12,1)=-2.0D0*alpha*Y(4)*growth1(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
     	F_THETA(12,2)=-2.0D0*alpha*Y(4)*growth1(Y(12),Csat(Temp(T)))/kg1
      F_THETA(12,3)=-alpha*Y(5)*growth2(Y(12),Csat(Temp(T)))*
     &          DLOG(DABS((Y(12)-Csat(Temp(T)))/Csat(Temp(T))))
	F_THETA(12,4)=-alpha*Y(5)*growth2(Y(12),Csat(Temp(T)))



*Sensitivity Equation:   W' = J W + d f/d theta
*where
* 			 d u/ d t = f(t,u;theta)     (solved above)
* 			 W =  d u /d theta           (unknown)
* 			 J = d f / du                (Jacobian matrix)

      DO 101 I = 1,NU
	DO 102 J = 1,NTHETA
           WPRIME(I,J)=F_THETA(I,J)
102     CONTINUE
101   CONTINUE

      DO 521 I = 1,NU
	DO 522 J = 1,NTHETA
	   DO 523 K = 1,NU
             WPRIME(I,J)=WPRIME(I,J)+JACOBIAN(I,K)*W(K,J)
523        CONTINUE
522     CONTINUE
521   CONTINUE

      K=20	
      DO 243 I = 1 , NU
         DO 244 J = 1 , NTHETA
		Yprime(K)=WPRIME(I,J)
		K=K+1
244	 CONTINUE
243   CONTINUE
      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 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_DATA1/kg1, g1

      if ((conc-concs) .gt. 0) then
	   growth1=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 growth3(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_DATA1/kg1, g1

      if ((conc-concs) .gt. 0) then
	     growth3 = kg1*((conc-concs)/concs)**(g1-1.0D0)
	  elseif ((conc-concs) .eq. 0) then
	     growth3 = 0.0D0
	  else
	     growth3 = -1.0D0*kg1*(-1.0D0*(conc-concs)/concs)**(g1-1.0D0)
	  endif


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

      REAL*8 FUNCTION growth4(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
	     growth4 = 5.0D0*kg2*((conc-concs)/concs)**(g2-1.0D0)
	  elseif ((conc-concs) .eq. 0) then
	     growth4 = 0.0D0
	  else
	     growth4 = -5.0D0*kg2*(-1.0D0*(conc-concs)/concs)**(g2-1.0D0)
	  endif


      RETURN
      END

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

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

************************************************************
      REAL*8 FUNCTION birth2(conc, concs,m21)

*     birthth 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
	     birth2 = kb*(((conc-concs)/concs)**(b-1.0D0))
	&			 *m21*(1.0D-4)**3.0D0
	  elseif ((conc-concs) .eq. 0) then
	     birth2 = 0.0D0
	  else
	     birth2 = -1.0D0*kb*((-1.0D0*(conc-concs)/concs)**(b-1.0D0))
	&	         *m21*(1.0D-4)**3.0D0
	  endif

			
      RETURN
      END