!=============================================================== ! SUBROUTINE ---- sub_eos_tillotson ---- === ! === ! Subroutine for calculating the pressure === ! and sound velocity from the Tillotson EOS === ! === ! Ver1--2006/07/09 === ! mod 2011/06/23 === !=============================================================== !=============================================================== ! Input Data : imat : Material ID ! rho : Density [kg/m^3] ! e : Energy [J/kg] ! Output Data : p : Pressure [Pa] ! sv : Sound Velocity [m/s] !=============================================================== subroutine sub_eos_tillotson(imat,rho,e,p,sv) use mod_parameter_tillotson_eos implicit none !---- Arguments in This Subroutine ----------------------------- integer imat ! Material ID real*8 rho ! Density [kg/m^3] real*8 e ! Internal Energy [J/kg] real*8 p ! Pressure [Pa] real*8 sv ! Sound Velocity [m/s] !---- Variables Used in This Subroutine ------------------------ real*8 sv2 ! Sound Velocity Squared [m^2/s^2] real*8 p_cold ! Pressure for Cold Expanded Region real*8 p_hot ! Pressure for Hot Expanded Region real*8 sv2_cold ! sv2 for Cold Expanded Region real*8 sv2_hot ! sv2 for Hot Expanded Region real*8 rho0,a,b,ap,bp,e0,alph,beta,eiv,ecv ! Parameter Set for Tillotson EOS !--------------------------------------------------------------- !==== Identify Material ======================================== if (imat .eq. 1) then ! ---- Any Material ---- rho0=any_rho0 a=any_a b=any_b ap=any_ap bp=any_bp e0=any_e0 alph=any_alph beta=any_beta eiv=any_eiv ecv=any_ecv elseif (imat .eq. 2) then ! ---- Murchison ---- rho0=murch_rho0 a=murch_a b=murch_b ap=murch_ap bp=murch_bp e0=murch_e0 alph=murch_alph beta=murch_beta eiv=murch_eiv ecv=murch_ecv elseif (imat .eq. 3) then ! ---- Hard Ryugu ---- rho0=hryugu_rho0 a=hryugu_a b=hryugu_b ap=hryugu_ap bp=hryugu_bp e0=hryugu_e0 alph=hryugu_alph beta=hryugu_beta eiv=hryugu_eiv ecv=hryugu_ecv elseif (imat .eq. 4) then ! ---- Soft Ryugu ---- rho0=sryugu_rho0 a=sryugu_a b=sryugu_b ap=sryugu_ap bp=sryugu_bp e0=sryugu_e0 alph=sryugu_alph beta=sryugu_beta eiv=sryugu_eiv ecv=sryugu_ecv else write(*,*) '*************************************' write(*,*) ' imat is invalid! Program stops!' write(*,*) '*************************************' stop endif !==== Negative Internal Energy ========================================= if (e .le. 0.0d0) then p=0.0d0 sv2=0.0d0 go to 999 end if !==== Compressed Region or Cold Expanded Region ======================== if (rho .ge. rho0 .or. e .le. eiv) then call sub_tillotson_cold(rho,e,rho0,a,b,ap,bp,e0,p,sv2) !ccccc ! if (rho .le. rho0*0.95d0) p=0.0d0 !ccccc go to 999 end if !==== Hot Expanded Region ============================================== if (e .ge. ecv) then call sub_tillotson_hot( & & rho,e,rho0,a,b,ap,e0,alph,beta,p,sv2) go to 999 end if !==== Transit Region =================================================== call sub_tillotson_cold( & & rho,e,rho0,a,b,ap,bp,e0,p_cold,sv2_cold) call sub_tillotson_hot( & & rho,e,rho0,a,b,ap,e0,alph,beta,p_hot,sv2_hot) p=((e-eiv)*p_hot+(ecv-e)*p_cold)/(ecv-eiv) sv2=((e-eiv)*sv2_hot+(ecv-e)*sv2_cold)/(ecv-eiv) !======================================================================= !======================================================================= 999 continue sv2=dmax1(sv2,0.25d0*ap/rho0) sv=dsqrt(sv2) return end !======================================================================= !=== Compressed Region or Cold Expanded Region === !======================================================================= subroutine sub_tillotson_cold( & rho,e,rho0,a,b,ap,bp,e0,p,sv2) use mod_parameter implicit none real*8 rho,e,rho0,a,b,ap,bp,e0,p,sv2 real*8 eta,um real*8 p1,c2,c3,c4,dpde,dpdd eta=rho/rho0 um=eta-1.0d0 p1=e/e0/eta**2+1.0d0 c2=p1**2 c3=p1-1.0d0 c4=b*e*(3.0d0*c3+1.0d0)/c2 dpde=a*rho+b*rho/c2 dpdd=a*e+ap/rho0+2.0d0*bp*um/rho0+c4 p=(a+b/(e/e0/eta/eta+1.0d0))*rho*e+ap*um+bp*um*um sv2=dpdd+dpde*p/rho**2 if (p .lt. 0.0d0 .or. rho .le. D_CUT_OFF*rho0) then p=0.0d0 sv2=0.0d0 endif return end !======================================================================= !=== Hot Expanded Region === !======================================================================= subroutine sub_tillotson_hot( & rho,e,rho0,a,b,ap,e0,alph,beta,p,sv2) implicit none real*8 rho,e,rho0,a,b,ap,e0,alph,beta,p,sv2 real*8 eta,um real*8 p1,p2,p3,p4,vow,p6,p7,c1,c2,c3,c4,c5,c6,dpdd,dpde eta=rho/rho0 um=eta-1.0d0 p1=e/e0/eta**2+1.0d0 p2=a*e*rho p3=b*e*rho/p1 p4=ap*um vow=1.0d0/eta p6=beta*(vow-1.0d0) p6=dmin1(p6,70.0d0) p6=dexp(-p6) p7=alph*(vow-1.0d0)**2 p7=dmin1(p7,70.0d0) p7=dexp(-p7) c1=a*e c2=p1**2 c3=p1-1.0d0 c4=b*e*(3.0d0*c3+1.0d0)/c2 c5=p6*p7*ap c6=2.0d0*alph*(vow-1.0d0) dpdd=c1+p7*c4+p7*p3*rho0*c6/rho**2 & +c5*(1.0d0/rho0+um*rho0/rho/rho*(c6+beta)) dpde=a*rho+p7/c2*b*rho p=p2+(p3+p4*p6)*p7 if (p .lt. 0.0d0) p=0.0d0 sv2=dpdd+dpde*p/rho**2 return end