#!/usr/bin/gnuplot # hugonita3_tillotson.plt # Hugoniota na grafu (v,p) pro Tillotsonovu stavovou rovnici. # Miroslav Broz (miroslav.broz@email.cz), May 6th 2009 # PROBLEM! Dekompresni adiabata vede na velke finalni objemy!!! # parametry pro vapenec, viz Kenkmann et al. (2003): rho_0 = 2700. # kg/m^3 A = 40.e9 # Pa B = 67.e9 # Pa C = 0. # Pa a = 0.5 b = 0.6 # <- POZOR! TENTO ČLEN JE PRO JEDNODUCHOST ZCELA VYNECHÁN!!! E_0 = 10.e6 # J/kg # hugonita eta(v,v0) = v0/v mu(v,v0) = eta(v,v0)-1. hugoniota(v,v0,p0) = ( p0 * ((2./a+1.)*v0 - v) + 2./a*v * (A*mu(v,v0) + B*mu(v,v0)**2)) / ( (2./a+1.)*v - v0) GPa = 1.e9 # Pa p_atm = 1.e5 # Pa v2 = 1./rho_0 # m^3/kg p2 = 1.*p_atm v1 = 1./3.*v2 # <- PRO v1 = 1/4 v2 JE KOMPRESE PRILIS VELKA (tlak je > 1000 GPa) p1 = hugoniota(v1,v2,p2) print "v2 = ", v2, " m^3/kg" print "p2 = ", p2, " Pa = ", p2/p_atm, " p_atm" print "v1 = ", v1, " m^3/kg" print "p1 = ", p1, " Pa = ", p1/GPa, " GPa = ", p1/p_atm, " p_atm" D = sqrt(v2**2 * (p1-p2)/(v2-v1)) D_u = sqrt(v1**2 * (p1-p2)/(v2-v1)) u = D - D_u print "\nD = ", D, " m/s" print "u = ", u, " m/s" print "D-u = ", D_u, " m/s" x1_ = v2 y1_ = p2 x2_ = v1 y2_ = p1 # Rayleighova linie rayleigh(x) = y1_ + (x-x1_)/(x2_-x1_) * (y2_-y1_) rayleigh2(x) = x > v1 ? (x < v2 ? rayleigh(x) : NaN) : NaN # adiabata #gamma = 1.4 # <- POZOR! TOTO BY MĚLO ODPOVÍDAT EXPANZI PEVNÉ LÁTKY!!! (nikoli idealniho plynu) #alpha = 1./(gamma-1.) #alpha = 5./2. alpha = 3./2. gamma = (alpha+1.)/alpha print "gamma = ", gamma print "alpha = ", alpha # p v^gamma = konst. # p1 v1^gamma = p2 v2^gamma # p1 = p2 (v2/v1)^gamma adiabata(v1,v2,p2) = p2 * (v2/v1)**gamma # hustota po adiabaticke dekompresi # p3 v3^gamma = p1 v1^gamma p3 = p2 v3 = (p1/p3)**(1./gamma) * v1 adiabata2(x,v2,p2) = x > v1 ? (x < v3 ? adiabata(x,v2,p2) : NaN) : NaN print "\nv3 = ", v3, " m^3/kg" print "p3 = ", p3, " Pa = ", p3/GPa, " p_atm" ######################################################################## set xl "{/Helvetica-Oblique V} = 1/{/Symbol-Oblique r} / 10^{{/Symbol -}3} m^3/kg" set yl "{/Helvetica-Oblique p} / GPa" offset +2,0 set xr [0:1.65*v2/1.e-3] y1 = 0. y2 = 1000. # GPa set yr [y1:y2] set mytics 2 set samples 1000 set nokey set style line 1 lt 1 lw 5 # hugoniota set style line 2 lt 7 lw 0.5 # adiabata set style line 20 lt 7 lw 3 set style line 3 lt 0 lw 1 set style line 4 lt 7 lw 1 pt 1 ps 2 set style line 5 lt 7 lw 0.5 # rayleigh set style line 50 lt 7 lw 3 set style line 10 lt 9 set encoding iso_8859_2 #set label "1 " at 0,1 right set label "{/Helvetica-Bold 1} " at v1/1.e-3,p1/GPa-0.04*(y2-y1) right set label " {/Helvetica-Bold 2}" at v2/1.e-3,p2/GPa+0.04*(y2-y1) left set label "{/Helvetica-Bold 3}" at v3/1.e-3,p3/GPa+0.05*(y2-y1) center #tmp = 0.50; set label "{/=14 hugoniota} " at tmp, hugoniota(tmp,v2,p2)/GPa tc lt 1 right #tmp = 0.25; set label "{/=14 adiabata} " at tmp, adiabata(tmp,v2,p2)/GPa tc lt 9 right #tmp = 1.25; set label "{/=14 dekompresní adiabata} " at tmp, adiabata(tmp,v1,p1)/GPa+0.10*(y2-y1) tc lt 7 center #tmp = 0.55; set label " {/=14 Rayleighova linie}" at tmp, rayleigh(tmp)/GPa left x3=0.55*v2; y3=rayleigh(x3)/GPa x4=0.54*v2; y4=rayleigh(x4)/GPa set arrow from x3/1.e-3,y3 to x4/1.e-3,y4 size 0.02,15 lw 3 lt 7 front x3=0.60*v2; y3=adiabata(x3,v1,p1)/GPa x4=0.61*v2; y4=adiabata(x4,v1,p1)/GPa set arrow from x3/1.e-3,y3 to x4/1.e-3,y4 size 0.02,15 lw 3 lt 7 front x3=1.20*v2; y3=adiabata(x3,v1,p1)/GPa x4=1.21*v2; y4=adiabata(x4,v1,p1)/GPa set arrow from x3/1.e-3,y3 to x4/1.e-3,y4 size 0.02,15 lw 3 lt 7 front set lmargin 6.9 set rmargin 2.2 set bmargin 3.1 set tmargin 0.5 set term table set out "hugoniota_tillotson.tmp1" p rayleigh2(x*1.e-3)/GPa set out "hugoniota_tillotson.tmp2" p adiabata2(x*1.e-3,v1,p1)/GPa set term wxt set out p \ "hugoniota_tillotson.tmp3" w filledcurve ls 10,\ hugoniota(x*1.e-3,v2,p2)/GPa ls 1,\ adiabata(x*1.e-3,v2,p2)/GPa ls 3,\ adiabata(x*1.e-3,v1,p1)/GPa ls 2,\ adiabata2(x*1.e-3,v1,p1)/GPa ls 20,\ rayleigh(x*1.e-3)/GPa ls 5,\ rayleigh2(x*1.e-3)/GPa ls 50,\ p2/GPa w l ls 3,\ "