📄 optexpsun.f
字号:
derivIg(I+1)=-const1*DEXP(-const1*Y(3))*Y(18)
derivIlnkg(I+1)=-const1*DEXP(-const1*Y(3))*Y(19)
derivIb(I+1)=-const1*DEXP(-const1*Y(3))*Y(20)
derivIlnkb(I+1)=-const1*DEXP(-const1*Y(3))*Y(21)
derivconcg(I+1)=Y(30)
derivconclnkg(I+1)=Y(31)
derivconcb(I+1)=Y(32)
derivconclnkb(I+1)=Y(33)
F(2*(I-1)+1,1)=derivIg(I+1)
F(2*(I-1)+1,2)=derivIlnkg(I+1)
F(2*(I-1)+1,3)=derivIb(I+1)
F(2*(I-1)+1,4)=derivIlnkb(I+1)
F(2*(I-1)+2,1)=derivconcg(I+1)
F(2*(I-1)+2,2)=derivconclnkg(I+1)
F(2*(I-1)+2,3)=derivconcb(I+1)
F(2*(I-1)+2,4)=derivconclnkb(I+1)
1200 CONTINUE
* WRITE(*,*) 'done'
* WRITE(*,*) 'done'
1 FORMAT(2(F13.6,3x))
2 FORMAT(I1)
3 FORMAT(2(E16.6,3x),2(F16.6,3x))
DO 1370 I = 1, NTHETA
DO 1375 J = 1, 2*(NSTEP-1), 2
FTV(I,J)=F(J,I)/(sigma_trans**2.D0)
FTV(I,J+1)=F(J+1,I)/(sigma_conc**2.0D0)
1375 CONTINUE
1370 CONTINUE
DO 1400 I = 1, NTHETA
DO 1500 J = 1, NTHETA
FTVF(I,J) = 0.0D0
DO 1600 K = 1, 2*(NSTEP-1)
FTVF(I,J)=FTVF(I,J)+FTV(I,K)*F(K,J)
1600 CONTINUE
1500 CONTINUE
1400 CONTINUE
* weight_mean_size=moment4(NSTEP)/moment3(NSTEP)
* cov=DSQRT(moment2(NSTEP)*moment0(NSTEP)/
* & (moment1(NSTEP))**2.0D0-1.0D0)
* mass_ratio=(moment3(NSTEP)-seed_moment3(NSTEP))/
* & seed_moment3(NSTEP)
* 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
*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*, "Determinant of FTVF = ", detFTVF
fj = -DLOG(detFTVF)
* PRINT*, "fj = ", fj
* WRITE(11,666) fj
666 FORMAT('-DLOG(detFTVF) = ',F13.6)
* DO 1482 I = 1, NSTEP
* IF(Temp(time(I)).GT.umax) THEN
* H=H+Kpenalty*(Temp(time(I))-umax)
* ELSE IF (Temp(time(I)).LT.umin) THEN
* H=H+Kpenalty*(umin-Temp(time(I)))
* ENDIF
*1482 CONTINUE
* DO 1300 I=1,NSTEP
* real_time(I)=sngl(time(I))
* real_temperature(I)=sngl(Temp(time(I)))
* real_concentration(I)=sngl(concentration(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))
* real_derivIg(I)=sngl(derivIg(I))
* real_derivIlnkg(I)=sngl(derivIlnkg(I))
* real_derivIb(I)=sngl(derivIb(I))
* real_derivIlnkb(I)=sngl(derivIlnkb(I))
* real_derivconcg(I)=sngl(derivconcg(I))
* real_derivconclnkg(I)=sngl(derivconclnkg(I))
* real_derivconcb(I)=sngl(derivconcb(I))
* real_derivconclnkb(I)=sngl(derivconclnkb(I))
*1300 CONTINUE
* DO 1850 I = 1 , NSTEP
* DO 1860 J = 1, NU
* DO 1870 K = 1, NTHETA
* real_derivtheta(I,J,K)=sngl(derivtheta(I,J,K))
*1870 CONTINUE
*1860 CONTINUE
*1850 CONTINUE
* PRINT*
* WRITE(*,11)FTVF(1,1), FTVF(1,2),FTVF(1,3),FTVF(1,4)
* WRITE(*,11)FTVF(2,1), FTVF(2,2),FTVF(2,3),FTVF(2,4)
* WRITE(*,11)FTVF(3,1), FTVF(3,2),FTVF(3,3),FTVF(3,4)
* WRITE(*,11)FTVF(4,1), FTVF(4,2),FTVF(4,3),FTVF(4,4)
*11 FORMAT(4(F16.4,2x))
* K=1
* DO 643 I = 1 , 2
* DO 653 J = 1, 2
* PNginverse(K)=FTVF(I,J)
* K=K+1
*653 CONTINUE
*643 CONTINUE
* K=1
* DO 663 I = 3 , 4
* DO 673 J = 3, 4
* PNbinverse(K)=FTVF(I,J)
* K=K+1
*673 CONTINUE
*663 CONTINUE
* CALL DEVCSF(2, PNginverse, 2, EVALg, EVECg, 2)
* CALL DEVCSF(2, PNbinverse, 2, EVALb, EVECb, 2)
* CALL DLINRG(NTHETA,FTVF,NTHETA, inv_FTVF,NTHETA)
* chi_squared=11.1
* g_interval=DSQRT(chi_squared*inv_FTVF(1,1))
* lnkg_interval=DSQRT(chi_squared*inv_FTVF(2,2))
* b_interval=DSQRT(chi_squared*inv_FTVF(3,3))
* lnkb_interval=DSQRT(chi_squared*inv_FTVF(4,4))
* PRINT*
* PRINT*,'g = ', g, '+/- ', g_interval
* PRINT*,'ln kg = ',DLOG(kg),'+/- ', lnkg_interval
* PRINT*,' b = ',b,'+/- ', b_interval
* PRINT*,'ln kb = ',DLOG(kb),'+/- ', lnkb_interval
* DO 762 K = 1, (NSTEP-1)
* DO 763 M = 1 , NU
* WRITE(23,12)derivtheta(K,M,1),derivtheta(K,M,2)/kg,
* & derivtheta(K,M,3), derivtheta(K,M,4)/kb
*763 CONTINUE
* WRITE(23,*)
*762 CONTINUE
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
* 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*,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)
INTEGER NEQ, I, J, K, nr
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)
* real randnum
* external rnun
COMMON /EXP_DATA/r0, alpha, mu00
COMMON /GROWTH_DATA/kg, g
COMMON /BIRTH_DATA/kb, b
NTHETA=4
NU=6
nr = 1
* call rnun(nr,randnum)
* if (randnum .lt. 0.001) then
* write(999,*) ' T = ', T
* write(999,*) ' Y= ', Y
* end if
*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*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
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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -