% vdwvol.m by: Edward T Peltzer, MBARI % revised: 01 Dec 98 % % CALCULATE THE VOLUME OF A GAS AT A GIVEN T & P % using the Van der Waals Equation of State. % % This equation works best at P = 1 atm or less. % % INPUT: Gas selector: 1 == methane % 2 == carbon dioxide % 3 == nitrogen % 4 == oxygen % 5 == argon % 6 == ethane % 7 == propane % 8 == n-butane % 9 == iso-butane % 10 == n-pentane % 11 == iso-pentane % 12 == n-hexane % % Temperature (T) in degrees C. % Pressure in kPa == P(mbar) / 10. % P(dbar) * 10. % % OUTPUT: Volume of 1 mole of gas [vpm] in liters. % % vpm = vdwvol(n,T,P). function [vpm]=vdwvol(n,T,P) % DEFINE CONSTANTS FOR EQUATION OF STATE A = [2.253 3.592 1.390 1.360 1.345 5.489 8.664 14.47 12.87 19.01 18.05 24.39]; B = [4.278 4.267 3.913 3.183 3.219 6.380 8.445 12.26 11.42 14.60 14.17 17.35]; Bc = B ./ 100; R = 0.082057; % CONVERT INPUT DATA TO APPROPRIATE UNITS Tc = T + 273.15; Pc = P ./ 101.325; % CHECK TO SEE IF PRESSURE 'IN-RANGE' if Pc > 5 disp(' ') disp(' ') disp('Warning: input pressure out-of-range,') disp(' answer may be incorrect!') end % CALCULATE COEFFICIENTS FOR CUBIC EQUATION C2 = 0 - ((Pc .* Bc(n)) + (R .* Tc)); C4 = 0 - A(n) .* Bc(n); % SOLVE CUBIC rts = roots([Pc C2 A(n) C4]); % FIND REAL ROOT == VPM L = length(rts); nreal = 0; for I = 1:L if rts(I) == real(rts(I)) nreal = nreal + 1; idx = I; end end if nreal == 1 vpm = rts(idx); else vdif = abs(rts - idgvol(T,P)); [MinDif,idx] = min(vdif); vpm = rts(idx); end