//====================================================== // Thermodynamic Package in java // Class Gas Properties of perfect gases // Dr. Turhan Coban // Ege University, School of Engineering, // Department of Mechanical Engineering // Bornova İZMİR - TURKEY // // email : turhan.coban //====================================================== // File Name : Gas.java // This file contains the Gas class // this class sets basic properties of perfect gases // required data is read from Gas.dat. // ===================================================== // Description : This file contains the gas class // class gas calculates thermophysical properties of // perfect gasses // following properties can be calculated // T() : Temperature degree K // h(T) : enthalpy KJ/kmol // hf : formation enthalpy KJ/kg // ht(T) : total enthalpy KJ/kg (h+hf) // M : molar mass kg/kmol // HT(T) : total enthalpy KJ : M*ht(T) // P() : presuure bar // s(T,P) : entropy KJ/kmol K // Cp(T) : specific heat at constant pressure KJ/kmol K // Cv(T) : specific heat at constant volume KJ/kg K // gamma(T): adiabatic constant Cp/Cv // u(T) : Internal energy KJ/kmol // c(T) : speed of sound m/s // vis(T) : viscosity Ns/m^2 // k(T) : thermal conductivity KJ/kg K // DATA FILE DEFINATION // gas datas are written in the data file "Gas.dat" // if gas data is not given in the data file, it can be curve fitted // and added to the data file. Additional curve fitting programs supplied // in Numerical Analysis package. Each data has the following form : //------------------ // gasName // n_equation M h0 hf sf // xa[0] xb[0] xc[0] xd[0] tl[0] th[0] // ......... // xa[n_equation-1] xb[n_equation-1]......th[n_equation-1] // n_vis // xvis[0] // ......... // xvis[n_vis-1] // n_k // xk[0] // ......... // xk[n_k-1] //------------------- // unit of the xa : kcal/kmol // note : if any curvefitting applied for a new gas temperature values // in K nad enthalpy values in the unit of Kcal.kmol should be supply // for the Cp curve fitting //============================================================ // VARIABLE IDENTIFICATION // all the variables that type is not defined is a double variable // PROTECTED VARIABLES : // xa,xb,xc,xd ,tl,th : double pointers. This values used to calculate // specific heat at constant pressure from the following equation : // Cp(T) = xa[i]+xb[i]*1e-3*T+xc[i]*1.0e5/T^2+xd[i]*1e-6*T^3 // where tl[i] <= T <= th[i] // n_equation : number of equations (xa,xb,xc,xd,tl,th) for a gas // xvis : real pointers to define viscosity according to formula : // vis(T)=sum(xvis(i)*T^i) , for(i=0;i=0;i--) { visg=visg*T+xvis[i]; } visg*=1.0e-7; } else visg=0; if(!SI) visg/=1.488; return visg; } //=================================================================== public double k(double T) { // thermal conductivity of the gas if(!SI) T/=1.8; double kg; if(n_k!=0.0) { int nk=n_k-1; kg=xk[nk]; for(int i=n_k-2;i>=0;i--) { kg+=kg*T+xk[i]; } kg*=1.0e-3; } else kg=0; if(!SI) kg/=1.731; return kg; } //=================================================================== public double h(double T) { if(!SI) T/=1.8; // enthalpy KJ/kmol //integration of function dh=Cp(T)*dT double hh = h0; for(int i=0;ith[i]) && (i== (n_equation-1) ) ) || ((T tl[i])) { hh+= xa[i]*(T- tl[i]) + xb[i]*1.0e-3/2.0*(T*T-tl[i]*tl[i]) - xc[i]*1e5*(1/T-1/tl[i]) + xd[i]*1e-6*(T*T*T-tl[i]*tl[i]*tl[i])/3.0; } else if(T>th[i]) { hh+= xa[i]*(th[i]- tl[i]) + xb[i]*1.0e-3/2.0*(th[i]*th[i] - tl[i]*tl[i]) - xc[i]*1e5*(1/th[i]-1/tl[i]) + xd[i]*1e-6*(th[i]*th[i]*th[i]-tl[i]*tl[i]*tl[i])/3.0; } } if(!SI) hh*=0.42992; if(!mole) hh/=M; return hh; } //=================================================================== public double ht( double t) { double hf1=hf; double h1,href; double Tref=298.0; if(!SI) {Tref=536.67;hf1=hf*0.42992;} h1=h(t); href=h(Tref); double ht= h(t)-h(Tref)+hf1; //System.out.println("h1="+h1+"href="+href+"hf1="+hf1+"ht="+ht); if(!mole) ht/=M; return ht; } public double s0t( double t) { double sf1=sf; double Tref=298.0; if(!SI) {Tref=536.67;sf1*=0.238846;} double st=s(t)-s(Tref)+sf; if(!mole) st/=M; return st; } public double st( double t,double p) { double sf1=sf; double s1,sref; double Tref=298.0,Pref=1.0; if(!SI) {Tref=536.67;sf1*=0.238846;Pref=14.503684;} s1=s(t,p); sref=s(Tref,Pref); double sst=s1-sref+sf1; if(!mole) sst/=M; //System.out.println("s1="+s1+"sref="+sref+"sf1="+sf1+"sst="+sst+"P="+p+"Pref="+Pref); return sst; } //=================================================================== public double H(double t) { return h(t)*N; } //=================================================================== public double HT(double t) { return ht(t)*N; } //=================================================================== public double u(double T) { // internal energy KJ/kmol // Integration of function du = Cv(T)*dT double ht=h(T); double R=8.3145; if(!SI) {R=1.986;} double at= R*T; if(!mole) at/=M; return ht-at; } //=================================================================== public double v(double T,double P) { // specific volume of the gas m^3/kmol if(!SI) {T/=1.8;P/=14.503684;} double vt= 8314.5*T/(P*1e5); if(!mole) vt/=M; if(!SI) vt*=16.016949; return vt; } public double v(double T) { if(!SI) {T/=1.8;} double P=1.0; // specific volume of the gas m^3/kmol double vt=8314.5*T/(P*1e5); if(!mole) vt/=M; if(!SI) vt*=16.016949; return vt; } //=================================================================== public double c(double t) { double g=gamma(t); if(!SI) t/=1.8; // speed of sound m/s double ct=Math.sqrt(8314.5/M*t*g); if(!SI) ct/=0.3048; return ct; } //=================================================================== public double si(double T, double P) { // entropy KJ/kmol K // integration of function // ds = Cp(T) * dt/T - R dP/P // 0 point at 298 K // si=s(t,P)-s0(298) if(!SI) {T/=1.8;P/=14.503684;} double sot=s0; if(!mole) sot/=M; if(!SI) sot*=0.2388444444; return s(T,P)-sot; } public double s(double T, double P) { //entropy KJ/kmol K //integration of function // ds = Cp(T) * dt/T - R dP/P // 0 point at 0 K (so is added up at 298 K) if(!SI) {T/=1.8;P/=14.503684;} double ss=s0; for(int i=0;i th[i] && i==n_equation - 1 ) ||( T < tl[i] && i==0 )) { ss+=xa[i]*Math.log(T/tl[i]) + xb[i]*1.0E-3*(T-tl[i]) - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i])) + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]); } else if((T <= th[i]) && (T > tl[i])) { ss+=xa[i]*Math.log(T/tl[i]) + xb[i]*1.0E-3*(T-tl[i]) - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i])) + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]); } else if( T > th[i] ) { ss+=xa[i]*Math.log(th[i]/tl[i]) + xb[i]*1.0E-3*(th[i]-tl[i]) - xc[i]*1e5/2.0*(1.0/(th[i]*th[i]) - 1.0/(tl[i]*tl[i])) + xd[i]*1e-6/2.0*(th[i]*th[i] - tl[i]*tl[i]); } } double sst= (ss-8.3145*Math.log(P)); if(!mole) sst/=M; if(!SI) sst*=0.238846; //System.out.println("sst="+sst+"s0="+s0+"sf="+sf); return sst; } public double S(double T, double P) { return s(T,P)*N; } public double Si(double T, double P) { return si(T,P)*N; } public double s(double T) { //entropy KJ/kmol K //integration of function // ds = Cp(T) * dt/T - R dP/P double P=1; if(!SI) {T/=1.8;} double ss=s0; for(int i=0;i th[i] && i==n_equation - 1 ) ||( T < tl[i] && i==0 )) { ss+=xa[i]*Math.log(T/tl[i]) + xb[i]*1.0E-3*(T-tl[i]) - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i])) + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]); } else if((T <= th[i]) && (T > tl[i])) { ss+=xa[i]*Math.log(T/tl[i]) + xb[i]*1.0E-3*(T-tl[i]) - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i])) + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]); } else if( T > th[i] ) { ss+=xa[i]*Math.log(th[i]/tl[i]) + xb[i]*1.0E-3*(th[i]-tl[i]) - xc[i]*1e5/2.0*(1.0/(th[i]*th[i]) - 1.0/(tl[i]*tl[i])) + xd[i]*1e-6/2.0*(th[i]*th[i] - tl[i]*tl[i]); } } if(!mole) ss/=M; if(!SI) ss*=0.238846; return ss; } public double si(double T) { return s(T)*N; } //=================================================================== public double s0(double T) { return s(T); } public double pr(double T) { double Tref, Rref; if(!SI) {Tref=491.67;Rref=1.986;} else {Tref=273.15;Rref=8.3145;} if(!mole) Rref/=M; return Math.exp((s(T)-s(Tref))/Rref); } public double vr(double T) { double Rref; if(!SI) {Rref=1.986;} else {Rref=8.3145;} return (Rref/M)*T/pr(T)*10; } //=================================================================== public double g(double T,double P) { return h(T)-T*s(T,P); } public double gt(double T,double P) { return ht(T)-T*st(T,P); } public double gt(double T) { double Pref; if(!SI) {Pref=14.503684;} else {Pref=1.0;} return ht(T)-T*st(T,Pref); } public double g(double T) { double Pref; if(!SI) {Pref=14.503684;} else {Pref=1.0;} return h(T)-T*s(T,Pref); } public double G(double T,double P) { return g(T,P)*N; } public double G(double T) { return g(T)*N; } public double GT(double T) { return gt(T)*N; } public double GT(double T,double P) { return gt(T,P)*N; } //=================================================================== public double g0(double T) { return h(T)-T*s0(T); } //=================================================================== public double Cp(double T) { //specific heat at constant pressure KJ/kmol K if(!SI) {T/=1.8;} double cp=0.0; for (int i=0;i th[i] && i==n_equation - 1 ) ||( T < tl[i] && i==0 )) { cp=xa[i]+xb[i]*1.0e-3*T+xc[i]*1.0e5/T/T+xd[i]*1.0e-6*T*T; break; } else if((T <= th[i]) && (T >= tl[i]) ) { cp=xa[i]+xb[i]*1.0e-3*T+xc[i]*1.0e5/T/T+xd[i]*1.0e-6*T*T; break; } } if(!mole) cp/=M; if(!SI) cp*=0.2388444444; return cp; } //=================================================================== public double Cv(double T) { //specific heat at constant volume KJ/kmol K double Rref; if(!SI) {Rref=1.986;} else {Rref=8.3145;} if(!mole) Rref/=M; double cv; cv=Cp(T) - Rref; return cv; } //=================================================================== public double gamma(double T) { //adiabatic constant return Cp(T)/Cv(T); } //=================================================================== public double T( char name,double y0,double p) { // name can have values h : for enthalpy // u : for internal energy // s : for entropy // v : specific volume // yo : the value of the variable given by variable name double t; if(!SI) {t=540;} else {t=300;} if(name=='v') { double R; if(!SI) {R=1545/144;} //R else {R=8.3145e3/1e5;} //K if(!mole) R/=M; t= p*y0/R; } else { double dt=0; int nmax=400; double tolerance=1.0e-8; for(int i=0;i