Third-order high-temperature Birch-Murnaghan equation of state

Solve for volume, and for integral of V dP, for solid phases at high pressure and temperature.

#pragma version=1.3

// set bracketing values and tolerance for volume
// absolute brackets will be (1 +/- kBracketFactor) * V0
// should be okay for EOS for solids
static constant kBracketFactor=0.1
static constant kTol=1e-5

// V = volume at P, T
// Tref = reference temperature
// V0 = volume at Tref
// thermal expansion coefficient = c1 + c2*T + c3/T + c4/T^2 + c5/T^0.5
// K0 is the bulk modulus at Tref
// Kprime is the pressure deriviative of bulk modulus
// dKdT is temperature derivative of bulk modulus

// For room temperature EOS, set a, b, c, and dKdT to zero.

// 3rd order high-temperature Birch-Murnaghan equation of state, returns P
threadsafe function HTBM(V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    variable V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT
   
    variable KT=K0+dKdT*(T-Tref)
    variable V0T=V0*exp(c1*(T-Tref)+0.5*c2*(T^2-Tref^2)+c3*ln(T/Tref)-c4*(T^-1-Tref^-1)+2*c5*(T^0.5-Tref^0.5))
    return 3/2*KT*((V0T/V)^(7/3) - (V0T/V)^(5/3))*(1 + 3/4*(Kprime-4)*((V0T/V)^(2/3) - 1))
end

// wrapper function for FindRoots
threadsafe function HTBM_wrapper(wave w, variable V)
    // w contains V0, T, Tref, a, b, c, dKdT, K0, Kprime
    return HTBM(V, w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9], w[10])
end

// solve for V that satisfies P=pressure
threadsafe function solveHTBM(pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT, [limH,limL])
    variable pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT, limL, limH
    limL = ParamIsDefault(limL) ? (1-kBracketFactor)*V0 : limL
    limH = ParamIsDefault(limH) ? (1+kBracketFactor)*V0 : limH
   
    Make /free/D w={temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT}
    FindRoots /H=(limH)/L=(limL) /Q /T=(kTol)/Z=(pressure) HTBM_wrapper, w
    return V_flag==0 ? V_Root : NaN
end

// calculate integral of V dP from P=P0 to P=pressure at T=temperature
threadsafe function integrateVolume(P0, pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    variable P0, pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT
   
    variable volume=solveHTBM(pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    return (pressure-P0)*volume-integrateP(volume, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
end

// isothermal integral of 3rd order Birch-Murnaghan wrt V from V0 to V at T=T
threadsafe function integrateP(V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
    variable V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT
   
    variable KT=K0+dKdT*(T-Tref)
    variable V0T=V0*exp(c1*(T-Tref)+0.5*c2*(T^2-Tref^2)+c3*ln(T/Tref)-c4*(T^-1-Tref^-1)+2*c5*(T^0.5-Tref^0.5))
   
    variable integral = 0
    integral = -(Kprime-4)*V0T^3*V^-2 - (2-3*(Kprime-4))*V0T^(7/3)*V^(-4/3) - (3*(Kprime-4)-4)*V0T^(5/3)*V^(-2/3)
    integral += V0T*((Kprime-4) + (2-3*(Kprime-4)) + (3*(Kprime-4)-4))
    integral *= 9/16*KT
    return integral
end

 

 

 

Well, yes, I should be able to provide plenty of examples.

I posted this after digging out my EOS-solving code to make a calculation that involved iron at high pressure. Of course iron tends to sit in the centre of planets where the VdP contribution to enthalpy is important...

I used values from here:

https://www.sciencedirect.com/science/article/abs/pii/S0012821X1300294X

And a more general collection of reference values can be found here:

https://agupubs.onlinelibrary.wiley.com/doi/book/10.1029/RF002

Incidentally, if you're interested in EOS for real fluids, I have code for plenty of EOS for water, including IAPWS, and some other formulations that are optimised for high pressure. I'd be happy to share.

Forum

Support

Gallery

Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More