- 积分
- 0
- 注册时间
- 2013-9-16
- 仿真币
-
- 最后登录
- 1970-1-1
|
参照关键字手册中自定义状态方程附录,自己编写的JWL状态方程子程序,能够顺利计算。但是为什么炸药还没开始起爆,就存在一个高压区呢?有没有做过类似工作的,能否帮忙看一下子程序?
c example vectorized user implementation of the jwl EOS
c
c*** variables
c lft,llt --- tt,crv,first and last indices into arrays
c iflag ----- =0 calculate bulk modulus
c =1 update pressure and energy
c cb -------- bulk modulus
c pnew ------ new pressure
c hist ------ history variables
c rho0 ------ reference density
c eosp ------ EOS constants
c specen ---- energy/reference volume
c df -------- volume ratio, v/v0 = rho0/rho
c dvol ------ change in volume over time step
c v0 -------- reference volume
c pc -------- pressure cut-off
c dt -------- time step size
c tt -------- current time
c crv ------- curve array
c first ----- logical .true. for tt,crv,first time step
c (for initialization of the history variables)
c
include 'nlqparm'
logical first
c
dimension cb(*),pnew(*),hist(nlq,*),eosp(*),
& specen(*),df(*),dvol(*),pc(*),v0(*)
c
b1=eosp(1)
b2=eosp(2)
r1=eosp(3)
r2=eosp(4)
w1=eosp(5)
bw1dr1=b1*w1/r1
bw2dr2=b2*w1/r2
c
c*** calculate the bulk modulus for the EOS contribution to the sound speed
if (iflag.eq.0) then
do i=lft,llt
r1v=bw1dr1/df(i)
r2v=bw2dr2/df(i)
wdr1v=b1-r1v
wdr2v=b2-r2v
dr1v=w1*specen(i)/df(i)
er1v=exp(-r1*df(i))
er2v=exp(-r2*df(i))
cb(i)=df(i)*df(i)*((wdr1v*r1+r1v/df(i))*er1v+(wdr2v*r2
& +r2v/df(i))*er2v+dr1v/df(i)+specen(i)*w1/(df(i)**2))
if (iburn.eq.3)then
cb(i)=max(1.e-09,cb(i))
endif
enddo
c
c*** update the pressure and internal energy
else
do i=lft,llt
r1v=bw1dr1/df(i)
r2v=bw2dr2/df(i)
wdr1v=b1-r1v
wdr2v=b2-r2v
dr1v=w1*specen(i)/df(i)
er1v=exp(-r1*df(i))
er2v=exp(-r2*df(i))
a=wdr1v*er1v+wdr2v*er2v
b=w1/df(i)
dvov0=0.5*dvol(i)/v0(i)
denom=1.+b*dvov0
pnew(i)=(a+specen(i)*b)/max(1.e-6,denom)
pnew(i)=max(pnew(i),pc(i))
specen(i)=specen(i)-pnew(i)*dvov0
enddo
endif
c
return
end
|
|