optexp.f

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

******************************************************
*
      time(1)=0.0D0
      T=0.0D0    
        
      DO 1200 I=1,(NSTEP-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)

         K = 10
         DO 263 M = 1 , NU
            DO 264 N = 1 , NTHETA
               Y(K)=derivtheta(I,M,N)
	       K = K + 1
264         CONTINUE
263      CONTINUE
	 
        CALL lsodes ( MOMENTS,NEQ,Y,T,T+1.0D0,itol,rtol,
     &          atol,itask,istate,iopt, 
     &           rwork,lrw,iwork,liw, MOMENTSJ, mf )

	
	   if(istate.EQ.-1) then
		  istate=1
	   else if (istate.LE.0) then
		  print*, "istate=", istate
	      pause
	   
	   end if
     
         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)=(Y(6)-Csat(Temp(T)))/Csat(Temp(T))
         time(I+1)=T
         K = 10
         DO 363 M = 1 , NU
            DO 364 N = 1 , NTHETA
               derivtheta(I+1,M,N)=Y(K)
	         F(Nu*(I-1)+M,N)=Y(k)
	         K = K + 1
364         CONTINUE
363      CONTINUE
1200  CONTINUE

      DO 1370 I = 1, NTHETA
         DO 1375 J = 1, NU*(NSTEP-1), Nu
            FTV(I,J)=F(J,I)/mu0_variance
            FTV(I,J+1)=F(J+1,I)/mu1_variance
            FTV(I,J+2)=F(J+2,I)/mu2_variance
            FTV(I,J+3)=F(J+3,I)/mu3_variance
            FTV(I,J+4)=F(J+4,I)/mu4_variance
            FTV(I,J+5)=F(J+5,I)/conc_variance
1375     CONTINUE
1370  CONTINUE

      DO 1400 I = 1, NTHETA
        DO 1500 J = 1, NTHETA
           FTVF(I,J) = 0.0D0
           DO 1600 K = 1, Nu*(NSTEP-1)
               FTVF(I,J)=FTVF(I,J)+FTV(I,K)*F(K,J)
1600       CONTINUE
1500    CONTINUE
1400  CONTINUE
	pause
	DO I=1,NTHETA
	 WRITE(*,1650) FTVF(I,1), FTVF(I,2), FTVF(I,3), FTVF(I,4)
	ENDDO
1650	FORMAT(4(E13.6, 1X))
*	PAUSE	 


*Calculate determinant of FTVF, 
*which is symmetric with dimensions NTHETA by NTHETA

      detFTVF =	(FTVF(1,3)*FTVF(2,4)-FTVF(1,4)*FTVF(2,3))**2.D0+
     &		(FTVF(1,3)*FTVF(3,4)-FTVF(3,3)*FTVF(1,4))*
     &			(FTVF(1,4)*FTVF(2,2)-FTVF(1,2)*FTVF(2,4))+
     &		(FTVF(2,3)*FTVF(3,4)-FTVF(3,3)*FTVF(2,4))*
     &			(FTVF(1,1)*FTVF(2,4)-FTVF(1,4)*FTVF(1,2))+
     &		(FTVF(1,3)*FTVF(4,4)-FTVF(3,4)*FTVF(1,4))*
     &			(FTVF(1,2)*FTVF(2,3)-FTVF(1,3)*FTVF(2,2))+
     &		(FTVF(2,3)*FTVF(4,4)-FTVF(3,4)*FTVF(2,4))*
     &			(FTVF(1,2)*FTVF(1,3)-FTVF(1,1)*FTVF(2,3))+
     &		(FTVF(3,3)*FTVF(4,4)-FTVF(3,4)*FTVF(3,4))*
     &			(FTVF(1,1)*FTVF(2,2)-FTVF(1,2)*FTVF(1,2))
	print*, 'detFTVF=', detFTVF
		fj =-DLOG(detFTVF)

*	print*, "fj=", fj

      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
      PARAMETER (Ntemp2 = 8)
      REAL*8 x(*),gj, Coeff(Ntemp2+3), umin, umax, Temp
      COMMON Coeff

      DO 628 I = 1, Nvar
         Coeff(I)=x(I)
628   CONTINUE
*	 print*, "in cntr"

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

      gj=umin-Temp(160.0D0)

*      IF(jjj.LE.Ntemp2) THEN
*         gj=Temp(160.0D0/DFLOAT(Ntemp2)*DFLOAT(jjj))-umax
*      ELSE
*         gj=umin-Temp(160.0D0/DFLOAT(Ntemp2)*DFLOAT(jjj-Ntemp2))
*      ENDIF

*	print*,"out cntr"

      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
	 
*	print*, "in Temp"

      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)
*	 print*, "out Temp"
 
      RETURN
      END
****************************************************************
      SUBROUTINE MOMENTS(NEQ,T,Y,Yprime)
      INTEGER NEQ, I, J, K
      INTEGER NTHETA, NU
      REAL*8 T, Y(NEQ), Yprime(NEQ)
      REAL*8 r0, alpha, mu00, kg, g, kb, b
      REAL*8 Csat, birth, growth, Temp
      REAL*8 JACOBIAN(6,6), W(6,4)
      REAL*8 WPRIME(6,4), F_THETA(6,4) 
      COMMON /EXP_DATA/r0, alpha, mu00
      COMMON /GROWTH_DATA/kg, g
      COMMON /BIRTH_DATA/kb, b
 
      NTHETA=4
      NU=6
*	print*, "in moments"

*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.0D0*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.0D0*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.0D0*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.0D0*growth(Y(6),Csat(Temp(T)))*Y(7)
*     d(mu'2]3(t)/dt), microns^3/g solvent/minute
      Yprime(9)=3.0D0*growth(Y(6),Csat(Temp(T)))*Y(8)

      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,4)=birth(Y(6),Csat(Temp(T)),Y(4))/Y(4)
      JACOBIAN(1,6)=b*birth(Y(6),Csat(Temp(T)),Y(4))/
     &			(Y(6)-Csat(Temp(T)))
      JACOBIAN(2,1)=growth(Y(6),Csat(Temp(T)))
      JACOBIAN(2,6)=g*growth(Y(6),Csat(Temp(T)))*Y(1)/
     &			(Y(6)-Csat(Temp(T)))
      JACOBIAN(3,2)=2.0D0*growth(Y(6),Csat(Temp(T)))
      JACOBIAN(3,6)=2.0D0*g*growth(Y(6),Csat(Temp(T)))*Y(2)/
     &			(Y(6)-Csat(Temp(T)))
      JACOBIAN(4,3)=3*growth(Y(6),Csat(Temp(T)))
      JACOBIAN(4,6)=3.0D0*g*growth(Y(6),Csat(Temp(T)))*Y(3)/
     &			(Y(6)-Csat(Temp(T)))
      JACOBIAN(5,4)=4*growth(Y(6),Csat(Temp(T)))
      JACOBIAN(5,6)=4.0D0*g*growth(Y(6),Csat(Temp(T)))*Y(4)/
     &			(Y(6)-Csat(Temp(T)))
      JACOBIAN(6,3)=-3.0D0*alpha*growth(Y(6),Csat(Temp(T)))
      JACOBIAN(6,6)=-3.0D0*alpha*g*growth(Y(6),Csat(Temp(T)))*Y(3)
     &			/(Y(6)-Csat(Temp(T)))

      K=10
      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,3)=birth(Y(6),Csat(Temp(T)),Y(4))*
     &		DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(1,4)=birth(Y(6), Csat(Temp(T)),Y(4))
      F_THETA(2,1)=Y(1)*growth(Y(6),Csat(Temp(T)))*
     &		DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(2,2)=Y(1)*growth(Y(6),Csat(Temp(T)))
      F_THETA(3,1)=2.0D0*Y(2)*growth(Y(6),Csat(Temp(T)))*
     &          DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(3,2)=2.0D0*Y(2)*growth(Y(6),Csat(Temp(T)))
      F_THETA(4,1)=3.0D0*Y(3)*growth(Y(6),Csat(Temp(T)))*
     &          DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(4,2)=3.0D0*Y(3)*growth(Y(6),Csat(Temp(T)))
      F_THETA(5,1)=4.0D0*Y(4)*growth(Y(6),Csat(Temp(T)))*
     &          DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(5,2)=4.0D0*Y(4)*growth(Y(6),Csat(Temp(T)))
      F_THETA(6,1)=-3.0D0*alpha*Y(3)*growth(Y(6),Csat(Temp(T)))*
     &          DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
      F_THETA(6,2)=-3.0D0*alpha*Y(3)*growth(Y(6),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=10	
      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
*	print*, "in Csat"
      Csat=0.1286D0+0.00588D0*T+0.0001721D0*T**2.0D0
	
*	print*, "Csat=", Csat

      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
	INTEGER fj_flag
      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)

*     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
	INTEGER fj_flag
      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