c
subroutine flow(xi,v,u,p,rho)
implicit double precision (a-h,o-z)
implicit integer (i-n)
c
include "cuser.i"
c
double precision lam
c
c Implementation of Chisnell's approximate solution of Guderley
c Shock implosion problem
c Source; Chisnell, R.F., JFM V354, pp357-375, 1998
c
c Copyright (C) 2003-2007 California Institute of Technology
c Dale Pullin
c
c Purpose: calculation of boundary conditions on
c r = r0 at t > t0
c
s = 2.
c
c s= 2 correponds to cylindrical flow
c Following numerical values are only for this case
c
c use variable xi to estimate v by interpolation
c al = 0.835319
c
v0 = al*0.653563
eta = s*al*gamma/2./(1.-al) -1.
eta = 1./eta
lam = (al/v0-1.)**2
c
q = -al/(v0*(1.-s*lam))
e = eta
d = (s-1.)/(1.-s*lam) - eta
f = al - (1.-lam)/(1.-s*lam)
b = eta+(s-1.)*(2.*lam+gamma-1.)/(1.-s*lam)
c
c Post shock conditions; xi = 1.0
c
c0 = sqrt(gamma)
rhos = (gamma+1.)*ms**2/((gamma-1.)*ms**2 + 2.)
us = -ms*c0*(2.*(ms**2+1.)/((gamma+1.)*ms**2))
ps = p0*(1.+2.*gamma/(gamma+1.)*(ms*ms-1.))
fact = 2.*(gamma-1.)/(gamma+1.)**2*(gamma*ms*ms+1.)
cs = c0*sqrt(fact*(ms*ms-1.)/(ms**2) + 1.0)
temps = temp0*(cs/c0)**2
gs = rhos/rho0
vs = t0*us/r0
zs = cs*cs*t0*t0/(r0*r0)
c
c use variable xi; put y = 1./xi, y=1-> shock. y=0-> infinity
c
c Estimate value of v by interpolation
c
y = 1./xi
n = 10002
dy = 1./n
ky = (1.-y)/dy
del = 1.-y - ky*dy
if(ky.ge.0.and.ky.le.n-3)then
v = var(ky) + del/dy*(var(ky+1)-var(ky))
else
print*,y,dy,del,ky,' y,dy,del,ky'
stop
endif
22 format(2x,i5,3f12.4,' xi,v,vs')
do j = 1,10
fun = ((q + v)/(q + vs))**f
fun = (vs/v)**al*((q + v)/(q + vs))**f - xi
dfun = ((f*v - al*(q + v))*(vs/v)**al*
& ((q + v)/(q + vs))**f)/(v*(q + v))
c print*,j,xi,fun,dfun,' j,xi,fun,dfun'
v = v - fun/dfun
if(abs(fun).le.1.d-10)then
go to 10
endif
enddo
print*,ky,xi,v,' i, xi,v'
stop
10 continue
c print(22),j,xi,v,vs
g = ((al-v)/(al-vs))**eta
g = g*((v+q)/(vs+q))**d
rho = g
u = xi*v/vs
p = (v/vs)**(2.*(1.-al))*((v+q)/(vs+q))**(b+d+2.*f)
c print*, ky,xi,r,v,rho,u,p,' i,xi,r,v,rho,u,p, ratio'
rho = rho*rhos
u = u*us
p = p*ps
c print*, ky,xi,r,v,rho,u,p,' i,xi,r,v,rho,u,p, physical'
100 format(2x,i5,5f12.4,e12.4,f12.4)
return
end
c
subroutine guderley()
implicit double precision (a-h,o-z)
implicit integer (i-n)
c
include "cuser.i"
c
double precision lam
c
c open(unit=8,file='CylV1.2_sh5.0_ratio.dat',status = 'unknown')
c open(unit=10,file='CylV1.2_shM5.0_phy.dat',status = 'unknown')
c
c Implementation of Chisnell's approximate solution of Guderley
c Shock implosion problem
c Source; Chisnell, R.F., JFM V354, pp357-375, 1998
c
c Purpose: to setup xi and v arrays for the implementation
c of a boundary condition at ro = 1.
c The shock is at r =r0 at t = t0, with Mach number ms
c to <0 must be calculated and returned
c
s = 2.
c
c s= 2 correponds to cylindrical flow
c Following numerical values are only for this case
c
al = 0.835319
v0 = al*0.653563
eta = s*al*gamma/2./(1.-al) -1.
eta = 1./eta
lam = (al/v0-1.)**2
c print*,al,' al'
c print*,v0,' v0'
c print*,eta,' eta'
c print*,lam,' lam'
c print*,gamma,' gamma'
c
q = -al/(v0*(1.-s*lam))
e = eta
d = (s-1.)/(1.-s*lam) - eta
f = al - (1.-lam)/(1.-s*lam)
b = eta+(s-1.)*(2.*lam+gamma-1.)/(1.-s*lam)
c print*,q,' q'
c print*,e,' e'
c print*,d,' d'
c print*,f,' f'
c print*,b,' b'
c
c pre-shock conditions
c
c0 = sqrt(gamma)
rhos = (gamma+1.)*ms**2/((gamma-1.)*ms**2 + 2.)
us = -ms*c0*(2.*(ms**2+1.)/((gamma+1.)*ms**2))
ps = p0*(1.+2.*gamma/(gamma+1.)*(ms*ms-1.))
fact = 2.*(gamma-1.)/(gamma+1.)**2*(gamma*ms*ms+1.)
cs = c0*sqrt(fact*(ms*ms-1.)/(ms**2) + 1.0)
temps = temp0*(cs/c0)**2
gs = rhos/rho0
c print*,cs,' cs'
c print*,rhos,' rhos'
c print*,us,' us'
c print*,fact,' fact'
c print*,temps,' temps'
c
c Choose shock radius and time
c
t0 = -al*r0/(ms*c0)
r0dot = -r0/abs(t0)*(-t0/abs(t0))**al
c print*,t0,' t0'
vs = t0*us/r0
zs = cs*cs*t0*t0/(r0*r0)
c print*,vs,' vs'
c
c use variable xi; put y = 1./xi, y=1-> shock. y=0-> infinity
c
c
c use variable xi; put y = 1./xi, y=1-> shock. y=0-> infinity
c
v = vs
c print*,v,' v'
n = 10002
dy = 1./n
i = 0
y = 1.
xi = 1./y
r = 1.
v = vs
rho=1.
u = 1.
p = 1.
rhoinv = 1./rho
yar(0) = y
var(0) = vs
uar(0) = u
par(0) = p
rhar(0) = rho
c
c write(8,100)i,xi,r,v,rho,u,p
rho = rho*rhos
u = u*us
p = p*ps
rhoinv = 1./rho
mach = abs(us)/sqrt(gamma*p/rho)
c write(10,100)i,xi,r,v,rho,u,p,mach
c
do i = 1,n-2
y = 1.- i*dy
xi = 1./y
r = xi*r0
do j = 1,10
c print*,v,vs,' v,vs'
c print*,q,al,' q,al'
fun = ((q + v)/(q + vs))**f
fun = (vs/v)**al*((q + v)/(q + vs))**f - xi
dfun = ((f*v - al*(q + v))*(vs/v)**al*
& ((q + v)/(q + vs))**f)/(v*(q + v))
c print*,j,xi,fun,dfun,' j,xi,fun,dfun'
v = v - fun/dfun
if(abs(fun).le.1.d-10)then
go to 10
endif
enddo
print*,i,xi,v,' i, xi,v'
stop
10 continue
c print*,al,vs,' al,vs'
c print*,al,v,' al,v'
g = ((al-v)/(al-vs))**eta
g = g*((v+q)/(vs+q))**d
rho = g
u = xi*v/vs
p = (v/vs)**(2.*(1.-al))*((v+q)/(vs+q))**(b+d+2.*f)
yar(i) = y
var(i) = v
uar(i) = u
par(i) = p
rhar(i) = rho
c
rhoinv = 1./rho
c write(8,100)i,xi,r,v,rho,u,p
rho = rho*rhos
u = u*us
p = p*ps
rhoinv = 1./rho
mach = abs(us)/sqrt(gamma*p/rho)
c write(10,100)i,xi,r,v,rho,u,p,mach
enddo
100 format(2x,i5,5f12.4,e12.4,f12.4)
return
end