c =====================================================
subroutine combl()
c =====================================================
c
c Create and initialize application specific common-blocks.
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
implicit double precision (a-h,o-z)
c
include "ck.i"
include "cuser.i"
parameter( lin = 5, lrin = 13, lwa = 50 )
dimension Xk(leNsp), Ywk(leNsp), Xkp(leNsp)
double precision Lind, Lpw
c
call cheminit()
c
open(unit=lin,status='old',form='formatted',file='init.dat')
read (lin, *) T0, p0
read (lin, *) sloc, udet, Nstat
read (lin, *) (Xk(k),k=1,Nsp)
read (lin, *) Tp, pp
read (lin, *) x1p, x2p
read (lin, *) (Xkp(k),k=1,Nsp)
read (lin, *) NCal,xe
close (lin)
c
c # Convert from SI-units into Chemkin units
p0 = p0 * 1.d1
udet = udet * 1.d2
c
u0 = 0.d0
c # Use NStat.eq.1 for simulation in Galilean Coordinates
if (Nstat.eq.1) u0 = u0 - udet
c
c # data in unburned state:
sum = 0.d0
do k = 1, Nsp
sum = sum + Xk(k)
enddo
do k=1,Nsp
Xk(k) = Xk(k)/sum
enddo
W0 = 0.d0
do k = 1,Nsp
W0 = W0 + Xk(k)*Wk(k)
enddo
do k = 1,Nsp
Yk(k) = Xk(k)*Wk(k)/W0
enddo
rho0 = (p0*W0)/(T0*RU)
rhou0 = rho0 * u0
h0 = avgtabip(T0,Yk,hms,Nsp)
E0 = rho0 * (0.5d0*u0**2 + h0) - p0
c
c # data in pocket:
sum = 0.d0
do k = 1, Nsp
sum = sum + Xkp(k)
enddo
do k = 1,Nsp
Xkp(k) = Xkp(k)/sum
enddo
Wp = 0.d0
do k = 1,Nsp
Wp = Wp + Xkp(k)*Wk(k)
enddo
do k = 1,Nsp
Ykp(k) = Xkp(k)*Wk(k)/Wp
enddo
rhop = (pp*Wp)/(Tp*RU)
c
call RankineH(Nsp,Yk,udet,rhov,uv,pv,Tv,ifail)
write (6,400) rhov*1.d3,pv*1.d-1,uv*1.d-2,Tv
400 format(' von Neumann-point: rho=',f10.8,' p=',f10.1,
& ' u=',f10.3,' T=',f10.5)
c
call IndLength(Lind, Lpw, xe)
write (6,420) Lind, Lpw, Lind/Lpw
420 format(' Induction length=',f10.8, ' Pulse width=', f10.8,
& ' Ratio=', f10.8)
return
end
c
c # For given mass fractions Y and detonation velocity udet
c # find rho,u,p,T such that these values satisfy the
c # Rankine-Hugoniot conditions with Yk,rho0,u0,p0,T0. See
c # R. Deiterding, Parallel adaptive simulation of
c # multi-dimensional detonation structures, page 43.
c =====================================================
subroutine RankineH(n,Y,D,rho,u,p,T,ifcnfail)
c =====================================================
implicit double precision (a-h,o-z)
c
include "ck.i"
include "cuser.i"
common /Rankine/ Drh, Yrh(leNsp), Wrh
dimension Y(n)
external hcurve
c
Wrh = 0.d0
do k=1,Nsp
Yrh(k) = Y(k)
Wrh = Wrh + Y(k)/Wk(k)
enddo
Wrh = 1.d0/Wrh
Drh = D
c
b = 1.d-16
c = 1.01d0
re = 1.d-12
ae = 1.d-08
call dzeroin(hcurve,b,c,re,ae,iflag)
eta = b
if (iflag.ne.1.and.iflag.ne.4) then
ifcnfail = 1
return
endif
c
b = 1.d-12
c = eta-ae
call dzeroin(hcurve,b,c,re,ae,iflag)
if (iflag.eq.1.or.iflag.eq.4) eta = dmin1(eta,b)
c
c # Eqs. (3.34) and (3.35)
rho = rho0 / eta
u = D*eta
p = p0 + rho0*D*D*(1.d0-eta)
T = (p*Wrh)/(rho*RU)
c
ifcnfail = 0
c
return
end
c
c # Eq. (3.36) has to be solved implicitly.
c =====================================================
double precision function hcurve(eta)
c =====================================================
implicit double precision (a-h,o-z)
c
include "ck.i"
include "cuser.i"
common /Rankine/ Drh, Yrh(leNsp), Wrh
c
h = 0.d0
T = Wrh/RU*((T0*RU/W0+Drh*Drh)*eta - (Drh*eta)**2)
itemp = int(T*TABFAC)
if (itemp.ge.TabS.and.itemp+1.le.TabE)
& h = avgtabip(T,Yrh,hms,Nsp)
hcurve = h - h0 + 0.5d0*Drh*Drh * (eta*eta-1.d0)
c
return
end
c
c =====================================================
c # In order to calculate the induction length
c # now compute the ZND structure first
c ### Xind=sloc-x(MAX dT/dx) ###
c # Calculate detonation profile by integrating
c # production rates in space. Detonation is assumed
c # to be traveling from left to right.
c =====================================================
subroutine IndLength(Lind, Lpw, xe)
c =====================================================
implicit double precision (a-h,o-z)
c
common /grkstat/ jsuc,jfail,jf,ja
include "ck.i"
include "cuser.i"
c
parameter ( lengrkaux = 1 )
dimension grkaux(lengrkaux), scale(LeNsp), qint(LeNsp)
external fcn, Jfcn, odeout
c
parameter ( lengN = 5000 )
dimension dTdx(LengN)
double precision Lind, Lpw, Lind_1, Lind_2
c
dx=xe/lengN
c
tol = 1.d-05
sclimit = 1.d-03
uround = 5.d-15
tst = dx
ts = 0.d0
te = 0.d0
tout = 0.d0
ichan = 0
ijac = 0
iscale = 0
jsuc = 0
jfail = 0
jf = 0
ja = 0
Nint = Nsp
hmax = dx
c
do k=1, Nsp
qint(k) = Yk(k)
enddo
c # Start from the von Neumann states
Told =Tv
c # Fill current grid from right to left and integrate
c # production function by going from cell to cell.
do 150 i=1,LengN
c
ts = (i-1)*dx
te = ts+dx
c
call grk4a(Nint,fcn,Jfcn,ts,qint,te,tol,
& hmax,tst,lengrkaux,grkaux,ichan,ijac,odeout,
& iscale,scale,lengrkw,grkwork,lenint,igrkwork)
call RankineH(Nsp,qint,udet,rho,u,p,T,ifail)
dTdx(i)=(T-Told)/dx
Told=T
150 continue
c # Find the location of the Maxmim Temperature Gradient
dTdx_max=0.0
do 160 i=1,LengN
if (dTdx(i) .ge. dTdx_max) then
dTdx_max=dTdx(i)
Lind=(i-0.5d0)*dx
k_max=i
endif
160 continue
c # Find the location of the half Maxmim Temperature Gradient
do 170 i=1,k_max
if (dTdx(i) .le. 0.5d0*dTdx_max) Lind_1=(i-0.5d0)*dx
170 continue
do 180 i=k_max, LengN
if (dTdx(i) .ge. 0.5d0*dTdx_max) Lind_2=(i-0.5d0)*dx
180 continue
Lpw=Lind_2-lind_1
return
end