找回密码
 注册
Simdroid-非首页
查看: 666|回复: 4

[二次开发及其理论] JWL 状态方程源代码

[复制链接]
发表于 2014-4-27 15:17:39 | 显示全部楼层 |阅读模式 来自 北京
    参照关键字手册中自定义状态方程附录,自己编写的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
发表于 2016-10-14 02:15:40 | 显示全部楼层 来自 广东广州
Simdroid开发平台
是fortran语言吗
回复 不支持

使用道具 举报

发表于 2017-11-13 21:24:22 | 显示全部楼层 来自 江苏南京
请问楼主,开发状态方程的入口在哪里?我查了查说状态方程是ueos,这是哪个文件里面?是和umat一个文件吗?我怎么没有找到?
回复 不支持

使用道具 举报

发表于 2022-3-7 10:17:21 | 显示全部楼层 来自 湖南湘潭
dabaiyang-m 发表于 2017-11-13 21:24
请问楼主,开发状态方程的入口在哪里?我查了查说状态方程是ueos,这是哪个文件里面?是和umat一个文件吗? ...

自定义状态方程在dyn21b.F中定义
回复 不支持

使用道具 举报

发表于 2023-9-20 11:16:34 | 显示全部楼层 来自 中国
你这个求体积模量CD的方程不对吧
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 15:31 , Processed in 0.035006 second(s), 11 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表