//====================================================== // java Termodinamik paketi // Class sogut soğutma akışkanların özellikleri // Dr. Turhan Coban // Ege Üniversitesi Mühendislik fakultesi // Makina Bölümü // //====================================================== // Dosya adı : sogut.java // Bu dosya soğutucu akışkanların termodinamik özelliklerini // Martin-Hou hal denklemiyle hesaplar // // ===================================================== // P : pressure as a function of temperature and specific volume // equation of state // Pi : pressure as a function of temperature and specific volume // also predict saturation presure if the volume in saturation region // dPdv : derivative of P with respect to v // dPdt : derivative of P with respect to t // property : method to calculate several thermodynamic properties // u : internal energy as a function of t and v // h : enthalpy as a function of t and v // s : entropy as a function of t and v // Ps : saturated vapor pressure as a function of T // Ts : saturated vapor temperature as a function of T // vsl_t :saturated liquid specific volume as a function of T // vsg_t :saturated vapor specific volume as a function of T // v_tp : specific volume as a function of T and v // v_th : specific volume as a function of T and h // v_tu : specific volume as a function of T and u // v_tx : specific volume as a function of T and x // v_ts : specific volume as a function of T and s // t_pv : temperature as a function of pressure and specific volume // t_ph : temperature as a function of pressure and enthalpy // t_pu : temperature as a function of pressure and internal energy // t_ps : temperature as a function of pressure and internal entropy // x_tv : quality as a function of temperature and specific volume // toString : result output through a string variable import java.io.*; import java.net.*; import java.util.*; class sogut { boolean SI=true; //switch for SI/EN units public sogutData rd; final String u[]={"SI","EN"}; String unit=u[0]; //=======constroctors======= public sogut() { readsogut("R22");} public sogut(String rName) { readsogut(rName);} public String readGasNames() { String s1=rd.refrigerantName;String s2=readsogut("null");readsogut(s1);return s2;} //property calculation methods=========== public double[] property_EN(String s, double v1, double v2) { if(s.charAt(0)=='t') v1=(v1-32.0)/1.8; else if(s.charAt(0)=='p') v1/=145.04; else if(s.charAt(0)=='v') v1/=16.01877892; if(s.charAt(1)=='t') v2=(v2-32)/1.8; else if(s.charAt(1)=='p') v2/=14.504; else if(s.charAt(1)=='v') v2/=16.01877892; else if(s.charAt(1)=='h') v2/=0.42992; else if(s.charAt(1)=='u') v2/=0.42992; else if(s.charAt(1)=='s') v2/=0.2388444444; double pp[]=property_SI(s,v1,v2); pp[0]*=145.04; pp[1]=pp[1]*1.8+32; pp[2]*=16.01877892; pp[3]*=0.42992; pp[4]*=0.42992; pp[5]*=0.2388444444; pp[7]=1.0/pp[2]; return pp; } public double[] property_SI(String s, double v1, double v2) { //---- cases of String s -------- // tv tp th tu ts tx // pv pt ph pu ps px // vp vt //---- Exit properties --------- // r[0] P Pressure MPa // r[1] t temperature degree C // r[2] v specific volume m^3/kg // r[3] h KJ/kg // r[4] u KJ/kg // r[5] s KJ/kgK // r[6] x kg vapor/kg total phase // r[7] ro density kg/m^3 double[] r=new double[8]; char s1=s.charAt(0); char s2=s.charAt(1); r[6]=x_tv(v1,v2); if(s.equals("tv")) //tv case { r[6]=x_tv(v1,v2); if(r[6]>=0 && r[6]<=1) { r[0]=Ps(v1); double vsl=vsl_t(v1); double vsg=vsg_t(v1); r[3]=h(v1,vsl)*(1-r[6])+ h(v1,vsl)*r[6]; r[4]=u(v1,vsl)*(1-r[6])+ u(v1,vsl)*r[6]; r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsl)*r[6]; } else { r[0]=P(v1,v2); r[3]=h(v1,v2); r[4]=u(v1,v2); r[5]=s(v1,v2); } r[1]=v1; r[2]=v2; } //end of tv case else if(s.equals("pv")) // pv case { r[6]=x_pv(v1,v2); if(r[6]>=0 && r[6]<=1) { r[1]=Ts(v1); double vsl=vsl_t(r[1]); double vsg=vsg_t(r[1]); r[3]=h(r[1],vsl)*(1-r[6])+ h(r[1],vsl)*r[6]; r[4]=u(r[1],vsl)*(1-r[6])+ u(r[1],vsl)*r[6]; r[5]=s(r[1],vsl)*(1-r[6])+ s(r[1],vsl)*r[6]; } else { r[1]=t_pv(v1,v2); r[3]=h(r[1],v2); r[4]=u(r[1],v2); r[5]=s(r[1],v2); } r[0]=v1; r[2]=v2; } //end of pv case else if(s.equals("tp")) //tp case { r[0]=v2; r[1]=v1; r[2]=v_tp(v1,v2); r[3]=h(r[1],r[2]); r[4]=u(r[1],r[2]); r[5]=s(r[1],r[2]); r[6]=x_tv(r[1],r[2]); } //end of tp case else if(s.equals("th")) { double hsl=h(v1,vsl_t(v1)); double hsg=h(v1,vsg_t(v1)); double vth=v_th(v1,v2); r[0]=P(v1,vth); r[1]=v1; r[2]=vth; r[3]=v2; r[4]=u(v1,vth); r[5]=s(v1,vth); r[6]=x_tv(v1,r[2]); if(r[6]>=0 && r[6]<=1) { r[0]=Ps(v1); r[1]=v1; double vsl=vsl_t(v1); double vsg=vsg_t(v1); r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=v2; r[4]=h(v1,vsl)*(1-r[6])+ h(v1,vsg)*r[6]; r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsg)*r[6]; } } //end of th case else if(s.equals("tu")) { double vtu=v_tu(v1,v2); r[0]=P(v1,vtu); r[1]=v1; r[2]=vtu; r[3]=h(v1,vtu); r[4]=v2; r[5]=s(v1,vtu); r[6]=x_tv(v1,r[2]); if(r[6]>=0 && r[6]<=1) { r[0]=Ps(v1); r[1]=v1; double vsl=vsl_t(v1); double vsg=vsg_t(v1); r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=h(v1,vsl)*(1-r[6])+ h(v1,vsg)*r[6]; r[4]=v2; r[5]=s(v1,vsl)*(1-r[6])+ s(v1,vsg)*r[6]; } } // end of tu case else if(s.equals("ts")) { double vts=v_ts(v1,v2); r[0]=P(v1,vts); r[1]=v1; r[2]=v_ts(v1,v2); r[3]=h(v1,vts); r[4]=u(v1,vts); r[5]=v2; r[6]=x_tv(v1,r[2]); if(r[6]>=0 && r[6]<=1) { r[0]=Ps(v1); double vsl=vsl_t(v1); double vsg=vsg_t(v1); r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=h(v1,vsl)*(1-r[6])+ h(v1,vsg)*r[6]; r[4]=u(v1,vsl)*(1-r[6])+ u(v1,vsg)*r[6]; } } //end of ts case else if(s.equals("tx")) { r[6]=v2; r[0]=Ps(v1); r[1]=v1; double hl=hsl_t(v1); double hg=hsg_t(v1); double vsl=vsl_t(v1); double vsg=vsg_t(v1); double ul=hl-r[0]*1e3*vsl; double ug=hg-r[0]*1e3*vsg; r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=hl*(1-r[6])+ hg*r[6]; r[4]=ul*(1-r[6])+ ug*r[6]; r[5]=ssl_t(v1)*(1-r[6])+ ssg_t(v1)*r[6]; } //end of tx case else if(s.equals("vt")) { r=property_SI("tv",v2,v1); } //end of vt case else if(s.equals("vp")) { r=property_SI("pv",v2,v1); } //end of vp case else if(s.equals("pt")) { r=property_SI("tp",v2,v1); } //end of tp case else if(s.equals("ph")) { double ts=Ts(v1); double vsl=vsl_t(ts); double vsg=vsg_t(ts); double hsl=h(ts,vsl); double hsg=h(ts,vsg); double vph=v_ph(v1,v2); r[0]=v1; r[1]=t_ph(v1,v2); r[2]=vph; r[3]=v2; r[4]=u(r[1],vph); r[5]=s(r[1],vph); if( v1 < rd.Pc) { r[6]=(v2-hsl)/(hsg-hsl);} else r[6]=3; if(r[6]>=0 && r[6]<=1) { r[0]=v1; r[1]=ts; r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=v2; r[4]=u(ts,vsl)*(1-r[6])+ u(ts,vsg)*r[6]; r[5]=s(ts,vsl)*(1-r[6])+ s(ts,vsg)*r[6]; } } //end of ph case else if(s.equals("pu")) { double ts=Ts(v1); double vsl=vsl_t(ts); double vsg=vsg_t(ts); double usl=u(ts,vsl); double usg=u(ts,vsg); double vpu=v_pu(v1,v2); r[0]=v1; r[1]=t_pu(v1,v2); r[2]=vpu; r[3]=h(r[1],r[2]); r[4]=v2; r[5]=s(r[1],r[2]); if( v1 < rd.Pc) { r[6]=(v2-usl)/(usg-usl);} else r[6]=3; if(r[6]>=0 && r[6]<=1) { r[0]=v1; r[1]=ts; r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=h(ts,vsl)*(1-r[6])+ h(ts,vsg)*r[6]; r[4]=v2; r[5]=s(ts,vsl)*(1-r[6])+ s(ts,vsg)*r[6]; } } //end of pu case else if(s.equals("ps")) { double ts=Ts(v1); double vsl=vsl_t(ts); double vsg=vsg_t(ts); double ssl=s(ts,vsl); double ssg=s(ts,vsg); double vps=v_ps(v1,v2); r[0]=v1; r[1]=t_ps(v1,v2); r[2]=vps; r[3]=h(r[1],r[2]); r[4]=u(r[1],r[2]); r[5]=v2; if( v1 < rd.Pc) { r[6]=(v2-ssl)/(ssg-ssl);} else r[6]=3; if(r[6]>=0 && r[6]<=1) { r[0]=v1; r[1]=ts; r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=h(ts,vsl)*(1-r[6])+ h(ts,vsg)*r[6]; r[4]=u(ts,vsl)*(1-r[6])+ u(ts,vsg)*r[6]; r[5]=v2; } } //end of ps case else if(s.equals("px")) { r[6]=v2; double ts=Ts(v1); r[0]=v1; r[1]=ts; double hl=hsl_t(ts); double hg=hsg_t(ts); double vsl=vsl_t(ts); double vsg=vsg_t(ts); double ul=hl-r[0]*1e3*vsl; double ug=hg-r[0]*1e3*vsg; r[2]=vsl*(1-r[6])+ vsg*r[6]; r[3]=hl*(1-r[6])+ hg*r[6]; r[4]=ul*(1-r[6])+ ug*r[6]; r[5]=ssl_t(v1)*(1-r[6])+ ssg_t(v1)*r[6]; } //end of px case r[7]=1.0/r[2]; return r; } public double[] property(String s, double v1, double v2) { if(unit.equals(u[0])) return property_SI(s,v1,v2); else return property_EN(s,v1,v2); } // ===== Equation of State P(t,v) ======= public double P(double t, double v) { //basic equation of state double T=t+273.15; double V=v*rd.M; double pi; double gamma=rd.K*T/rd.Tc; double kk=rd.u*V; double eg=Math.exp(-gamma); pi=rd.R*T/(V-rd.b); double mult=(V-rd.b); for(int i=1;i<5;i++) {mult*=(V-rd.b); pi+=(rd.A[i]+rd.B[i]*T+rd.C[i]*eg)/mult;} if(kk < 200) { double euV=Math.exp(rd.u*V); double euV1; if(rd.CP==0) euV1=euV; else euV1=euV*(1+rd.CP*euV); pi+=(rd.A[5]+rd.B[5]*T+rd.C[5]*eg)/euV1; } return pi; } double Pi(double t, double v) { if(vsl_t(t)<=v && v>=vsg_t(t)) {return Ps(t);} return P(t,v); } //==== derivatives of equation of state ==== public double dPdT(double t, double v) { //derivative of equation of state dP/dT double T=t+273.15; double V=v*rd.M; double kk=rd.u*V; double pi; double gamma=rd.K*T/rd.Tc; double ga=rd.K/rd.Tc; double eg=Math.exp(-gamma); pi=rd.R/(V-rd.b); double mult=(V-rd.b); for(int i=1;i<5;i++) {mult*=(V-rd.b); pi+=(rd.B[i]-ga*rd.C[i]*eg)/mult;} if(kk < 200) { double euV=Math.exp(rd.u*V); double euV1; if(rd.CP==0) euV1=euV; else euV1=euV*(1+rd.CP*euV); pi+=(rd.B[5]*T-rd.C[5]*ga*eg)/euV1; } return pi; } public double dPdV1(double t,double x) { // numerically calculated dP/dv double h0=0.0256; int i,m; int n=7; //derivative of a simple function double T[][]; T=new double[n][n]; double h[]; h=new double[n]; //vector h(n,0); for(i=0;i 0.0) System.out.println("xxx Root must be bracketed in ZBRENT"); fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=P(t,b)-Pi; } System.out.println("Maximum number of iterations exceeded in BRENT"); return b; } public double[] v_th_braket(double t,double hi,double x1,double x2) { //returns a bracketed root //if a root is not bracketed, enlarge region outwards each time int NTRY=100; double a[]=new double[2]; double FACTOR=1.343253; int j; double f1,f2; if (x1 == x2) System.out.println("Bad initial range in ZBRAC"); f1=h(t,x1)-hi; f2=h(t,x2)-hi; for (j=1;j<=NTRY;j++) { if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;} if (Math.abs(f1) < Math.abs(f2)) f1=h(t,(x1 += FACTOR*(x1-x2)))-hi; else f2=h(t,(x2 += FACTOR*(x2-x1)))-hi; } {a[0]=x1;a[1]=x2;return a;} } public double v_th(double t,double hi) { //specific volume as a function of temperature and enthalpy double T=t+273.15; double x=vg_app(t)+0.1; double dx=1.5436543624e-3; double x1=x-dx; double x2=x+dx; //find a bracket with the root double aa[]=v_th_braket(t,hi,x1,x2); x1=aa[0]; x2=aa[1]; double ITMAX=100; double EPS=3.0e-8; int iter; double a=x1,b=x2,c,d,e,min1,min2; double tol=1.0e-5; c=0; d=0; e=0; double fa=h(t,a)-hi,fb=h(t,b)-hi,fc,p,q,r,s,tol1,xm; if (fb*fa > 0.0) System.out.println("Root must be bracketed in ZBRENT"); fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=h(t,b)-hi; } System.out.println("Maximum number of iterations exceeded in BRENT"); return b; } public double[] v_tu_braket(double t,double ui,double x1,double x2) { //returns a bracketed root //if a root is not bracketed, enlarge region outwards each time int NTRY=50; double a[]=new double[2]; double FACTOR=1.343253; int j; double f1,f2; if (x1 == x2) System.out.println("Bad initial range in ZBRAC"); f1=u(t,x1)-ui; f2=u(t,x2)-ui; for (j=1;j<=NTRY;j++) { if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;} if (Math.abs(f1) < Math.abs(f2)) f1=u(t,(x1 += FACTOR*(x1-x2)))-ui; else f2=u(t,(x2 += FACTOR*(x2-x1)))-ui; } {a[0]=x1;a[1]=x2;return a;} } public double v_tu(double t,double ui) { //specific volume as a function of temperature and internal energy double T=t+273.15; double x=vg_app(t)+0.1; double dx=1.5436543624e-3; double x1=x-dx; double x2=x+dx; //find a bracket with the root double aa[]=v_tu_braket(t,ui,x1,x2); x1=aa[0]; x2=aa[1]; double ITMAX=100; double EPS=3.0e-8; int iter; double a=x1,b=x2,c,d,e,min1,min2; double tol=1.0e-5; c=0; d=0; e=0; double fa=u(t,a)-ui,fb=u(t,b)-ui,fc,p,q,r,s,tol1,xm; if (fb*fa > 0.0) System.out.println("Root must be bracketed in ZBRENT"); fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=u(t,b)-ui; } System.out.println("Maximum number of iterations exceeded in BRENT"); return b; } public double[] v_ts_braket(double t,double si,double x1,double x2) { //returns a bracketed root //if a root is not bracketed, enlarge region outwards each time int NTRY=50; double a[]=new double[2]; double FACTOR=1.343253; int j; double f1,f2; if (x1 == x2) System.out.println("Bad initial range in ZBRAC"); f1=s(t,x1)-si; f2=s(t,x2)-si; for (j=1;j<=NTRY;j++) { if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;} if (Math.abs(f1) < Math.abs(f2)) f1=s(t,(x1 += FACTOR*(x1-x2)))-si; else f2=s(t,(x2 += FACTOR*(x2-x1)))-si; } {a[0]=x1;a[1]=x2;return a;} } public double v_ts(double t,double si) { //specific volume as a function of temperature and entropy double T=t+273.15; double x=vg_app(t)+0.1; double dx=1.5436543624e-3; double x1=x-dx; double x2=x+dx; //find a bracket with the root double aa[]=v_ts_braket(t,si,x1,x2); x1=aa[0]; x2=aa[1]; double ITMAX=100; double EPS=3.0e-8; int iter; double a=x1,b=x2,c,d,e,min1,min2; double tol=1.0e-5; c=0; d=0; e=0; double fa=s(t,a)-si,fb=s(t,b)-si,fc,p,q,r,s,tol1,xm; if (fb*fa > 0.0) System.out.println("Root must be bracketed"); fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=s(t,b)-si; } System.out.println("Maximum number of iterations exceeded in BRENT"); return b; } public double v_tx(double t,double xi) { //specific volume as a function of temperature and quality return (1-xi)*vsl_t(t)+xi*vsg_t(t); } public double v_ph(double pi,double hi) { //specific volume as a function of pressure and enthalpy return v_tp(t_ph(pi,hi),pi); } public double v_pu(double pi,double ui) { //specific volume as a function of pressure and enthalpy return v_tp(t_pu(pi,ui),pi); } public double v_ps(double pi,double si) { //specific volume as a function of pressure and enthalpy return v_tp(t_ps(pi,si),pi); } // ======= t relations========== public double t_pv(double Pi,double v) { // x in f(x)=0 given also df(x)/dx // Newton - Raphson type method double V=v*rd.M; double x=(V-rd.b)*Pi/rd.R-273.15; System.out.println("starting value ="+x); int nmax=500; double tolerance=1.0e-15; double fx,dfx; for(int i=0;i 0.0) System.out.println("Root must be bracketed"); fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=h(b,v_tp(b,pi))-hi; } System.out.println("Maximum number of iterations exceeded in BRENT"); return b; } public double[] t_pu_braket(double pi,double ui,double x1,double x2) { //returns a bracketed root //if a root is not bracketed, enlarge region outwards each time int NTRY=50; double a[]=new double[2]; double FACTOR=1.343253; int j; double f1,f2; if (x1 == x2) System.out.println("Bad initial range in ZBRAC"); f1=u(x1,v_tp(x1,pi))-ui; f2=u(x2,v_tp(x2,pi))-ui; for (j=1;j<=NTRY;j++) { if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;} if (Math.abs(f1) < Math.abs(f2)) f1=u((x1 += FACTOR*(x1-x2)),v_tp(x1,pi))-ui; else f2=u((x2 += FACTOR*(x2-x1)),v_tp(x2,pi))-ui; } {a[0]=x1;a[1]=x2;return a;} } public double t_pu(double pi,double ui) { //specific volume as a function of pressure and enthalpy double x=Ts(pi); double dx=2.3e-1; double x1=x-dx; double x2=x+dx; //find a bracket with the root double aa[]=t_pu_braket(pi,ui,x1,x2); x1=aa[0]; x2=aa[1]; double ITMAX=100; double EPS=3.0e-8; int iter; double a=x1,b=x2,c,d,e,min1,min2; double tol=1.0e-5; c=0; d=0; e=0; double fa=u(a,v_tp(a,pi))-ui,fb=u(b,v_tp(b,pi))-ui,fc,p,q,r,s,tol1,xm; if (fb*fa > 0.0) System.out.println("Root must be bracketed"); fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=u(b,v_tp(b,pi))-ui; } System.out.println("Maximum number of iterations exceeded in BRENT"); return b; } public double[] t_ps_braket(double pi,double si,double x1,double x2) { //returns a bracketed root //if a root is not bracketed, enlarge region outwards each time int NTRY=50; double a[]=new double[2]; double FACTOR=1.343253; int j; double f1,f2; if (x1 == x2) System.out.println("Bad initial range in ZBRAC"); f1=s(x1,v_tp(x1,pi))-si; f2=s(x2,v_tp(x2,pi))-si; for (j=1;j<=NTRY;j++) { if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;} if (Math.abs(f1) < Math.abs(f2)) f1=s((x1 += FACTOR*(x1-x2)),v_tp(x1,pi))-si; else f2=s((x2 += FACTOR*(x2-x1)),v_tp(x2,pi))-si; } {a[0]=x1;a[1]=x2;return a;} } public double t_ps(double pi,double si) { //specific volume as a function of pressure and enthalpy double x=Ts(pi); double dx=2.3e-1; double x1=x-dx; double x2=x+dx; //find a bracket with the root double aa[]=t_ps_braket(pi,si,x1,x2); x1=aa[0]; x2=aa[1]; double ITMAX=100; double EPS=3.0e-8; int iter; double a=x1,b=x2,c,d,e,min1,min2; double tol=1.0e-5; c=0; d=0; e=0; double fa=s(a,v_tp(a,pi))-si,fb=s(b,v_tp(b,pi))-si,fc,p,q,r,s,tol1,xm; if (fb*fa > 0.0) System.out.println("Root must be bracketed"); fc=fb; for (iter=1;iter<=ITMAX;iter++) { if (fb*fc > 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=s(b,v_tp(b,pi))-si; } System.out.println("Maximum number of iterations exceeded in BRENT"); return b; } // ======= h relations ============= public double h1(double t, double v) { double T=t+273.15; double V=v*rd.M; double hh=0.0; double J=1e3; double gamma=rd.K*T/rd.Tc; double expgamma=Math.exp(-gamma); double ueuv=rd.u*Math.exp(rd.u*V); double mult=1.0; for(int i=1;i<5;i++) { mult*=(V-rd.b); hh+=J/(mult*i)*( rd.A[i] +(1.0+gamma)*expgamma*rd.C[i]); } if(rd.A[5]!=0) hh+=J*rd.A[5]/ueuv+J*P(t,v)*V; hh+=rd.G[0]*T+rd.G[1]*T*T/2.0+rd.G[2]*T*T*T/3.0+rd.G[3]*T*T*T*T/4.0 -rd.G[4]/T+rd.G[5]*T*T*T*T*T/5+rd.X; return hh/rd.M; } public double h(double t, double v) { double h1=h1(t,v); double x=x_tv(t,v); if(x>1) { h1=h1(t,v); } else if(x<0) { double P=P(t,v); h1=hsl_t(t)+vsl_t(t)*(P(t,v)-Ps(t))*1.0e3; } else if(x==0) { h1=hsl_t(t); } else if(x==1) { h1=hsg_t(t); } else if(x>0 || x<1) { h1=hsl_t(t)*(1-x)+hsg_t(t)*x; } return h1; } // ======= u relations ============= public double u(double t, double v) {return h(t,v)-1.0e3*P(t,v)*v;} // ======= s relations ============= public double s1(double t, double v) { double T=t+273.15; double V=v*rd.M; double ss; double J=1e3; double gamma=rd.K*T/rd.Tc; double expgamma=Math.exp(-gamma); double ueuv=rd.u*Math.exp(rd.u*V); double mult=1.0; ss=J*rd.R*Math.log(V-rd.b); for(int i=1;i<5;i++) { mult*=(V-rd.b); ss+=J/(mult*i)*(rd.B[i]-rd.K/rd.Tc*expgamma*rd.C[i]); } if(rd.B[5]!=0) ss+=(-J)*rd.B[5]/ueuv; ss+=rd.G[0]*Math.log(T)+rd.G[1]*T+rd.G[2]*T*T/2.0+ rd.G[3]*T*T*T/3.0+rd.G[4]/(2.0*T*T)+rd.G[5]*T*T*T*T/4.0+rd.Y; return ss/rd.M; } public double s(double t, double v) { double s1=s1(t,v); double x=x_tv(t,v); if(x>1) { s1=s1(t,v); } else if(x<0) { double P=P(t,v); s1=ssl_t(t)+vsl_t(t)*(P(t,v)-Ps(t))*1.0e3; } else if(x==0) { s1=ssl_t(t); } else if(x==1) { s1=ssg_t(t); } else if(x>0 || x<1) { s1=ssl_t(t)*(1-x)+ssg_t(t)*x; } return s1; } // ======= ideal gas relations ============= public double Cv(double t) {//ideal gas heat capacity double T=t+273.15; double cv0=rd.G[0]+rd.G[1]*T+rd.G[2]*T*T+rd.G[3]*T*T*T+rd.G[4]/(T*T)+rd.G[5]*T*T*T*T; return cv0; //KJ/kmol K } public double Cp(double t) {//ideal gas heat capacity double T=t+273.15; return Cv(T)+8.31434; //KJ/kmolK } // saturation state //========================================================= public String phase(double x) { //phase of the steam-water system if( x==3) return "kritik üstü sıvı"; else { if(x<0) return "sıvı"; else if( x==0.0) return "doymuş sıvı"; else if(x==1.0) return "doymuş buhar"; else if( x>0.0 && x<1.0) return "doymuş sıvı-buhar karışımı"; else // (x>1) return "kızgın buhar"; } } public String phase(double t, double v) { //phase of the vapour-liquid system if( Ps(t) > rd.Pc) return "kritik üstü sıvı"; else { double ts=t; double vsl=vsl_t(ts); double vsg=vsg_t(ts); if(vvsl && vvsg) return "kızgın buhar"; } } public double x_tv(double t, double v) { //phase and/or quality of the steam-water system double xx=0; double p=P(t,v); double Ts=Ts(p); if(t> rd.Tc) xx=3; else if(t > Ts) xx=2; else if( t < Ts) xx=-1; else if( t==Ts) { double vsl=vsl_t(t); double vsg=vsg_t(t); if(v==vsl) xx=0; else if(v==vsg) xx=1; else if(v>vsl && v rd.Tc) xx=3; else if(t > Ts) xx=2; else if( t < Ts) xx=-1; else if( t==Ts) { double vsl=vsl_t(t); double vsg=vsg_t(t); if(v==vsl) xx=0; else if(v==vsg) xx=1; else if(v>vsl && v= rd.Tc) return rd.Pc; double pi=0.0; if(t<=rd.th[0]) { pi=rd.aVP[0][0]+rd.aVP[0][1]/T+rd.aVP[0][2]/(T*T)+rd.aVP[0][3]*Math.log10(T)+rd.aVP[0][4]*T +rd.aVP[0][5]*T*T+rd.aVP[0][6]*T*T*T; if(rd.aVP[0][8]!=0) pi+=(rd.aVP[0][7]*(rd.aVP[0][8]-T)/T)*Math.log10((rd.aVP[0][8]-T)*MM); } else if(t>rd.th[0] && t<= rd.th[1]) { pi=rd.aVP[1][0]+rd.aVP[1][1]/T+rd.aVP[1][2]/(T*T)+rd.aVP[1][3]*Math.log10(T)+rd.aVP[1][4]*T+ rd.aVP[1][5]*T*T+rd.aVP[1][6]*T*T*T; if(rd.aVP[1][8]!=0) pi+=(rd.aVP[1][7]*(rd.aVP[1][8]-T)/T)*Math.log10((rd.aVP[1][8]-T)*MM); } else if(t>rd.th[1] && t<= rd.th[2]) { pi=rd.aVP[2][0]+rd.aVP[2][1]/T+rd.aVP[2][2]/(T*T)+rd.aVP[2][3]*Math.log10(T)+rd.aVP[2][4]*T+ rd.aVP[2][5]*T*T+rd.aVP[2][6]*T*T*T; if(rd.aVP[2][8]!=0) pi+=(rd.aVP[2][7]*(rd.aVP[2][8]-T)/T)*Math.log10((rd.aVP[2][8]-T)*MM); } else if(t>rd.th[2] && t<= rd.th[3]) { pi=rd.aVP[3][0]+rd.aVP[3][1]/T+rd.aVP[3][2]/(T*T)+rd.aVP[3][3]*Math.log10(T)+rd.aVP[3][4]*T+ rd.aVP[3][5]*T*T+rd.aVP[3][6]*T*T*T; if(rd.aVP[3][8]!=0) pi+=(rd.aVP[3][7]*(rd.aVP[3][8]-T)/T)*Math.log10((rd.aVP[3][8]-T)*MM); } else if(t>rd.th[3] && t<= rd.th[4]) { pi=rd.aVP[4][0]+rd.aVP[4][1]/T+rd.aVP[4][2]/(T*T)+rd.aVP[4][3]*Math.log10(T)+rd.aVP[4][4]*T+ rd.aVP[4][5]*T*T+rd.aVP[4][6]*T*T*T; if(rd.aVP[4][8]!=0) pi+=(rd.aVP[4][7]*(rd.aVP[4][8]-T)/T)*Math.log10((rd.aVP[4][8]-T)*MM); } else if(t> rd.th[4]) { pi=rd.aVP[5][0]+rd.aVP[5][1]/T+rd.aVP[5][2]/(T*T)+rd.aVP[5][3]*Math.log10(T)+rd.aVP[5][4]*T+ rd.aVP[5][5]*T*T+rd.aVP[5][6]*T*T*T; if(rd.aVP[5][8]!=0) pi+=(rd.aVP[5][7]*(rd.aVP[5][8]-T)/T)*Math.log10((rd.aVP[5][8]-T)*MM); } //System.out.println("pi="+pi); pi=Math.pow(10.0,pi); //System.out.println("pi="+pi); if(rd.refrigerantName=="R123" || rd.refrigerantName=="R134a") pi*=1e-3; return pi; //MPa } public double Ts(double Pi) { // x in f(x)=0 given also df(x)/dx // Newton - Raphson type method if( Pi >= rd.Pc) return rd.Tc; double x=Ts1(Pi); //double x=(rd.Tc-273.15); int nmax=500; double tolerance=1.0e-15; double fx,dfx; for(int i=0;i 0.0) { c=a; fc=fa; e=d=b-a; } if (Math.abs(fc) < Math.abs(fb)) { a=b; b=c; c=a; fa=fb; fb=fc; fc=fa; } tol1=2.0*EPS*Math.abs(b)+0.5*tol; xm=0.5*(c-b); if (Math.abs(xm) <= tol1 || fb == 0.0) return b; if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) { s=fb/fa; if (a == c) { p=2.0*xm*s; q=1.0-s; } else { q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0); } if (p > 0.0) q = -q; p=Math.abs(p); min1=3.0*xm*q-Math.abs(tol1*q); min2=Math.abs(e*q); if (2.0*p < (min1 < min2 ? min1 : min2)) { e=d; d=p/q; } else { d=xm; e=d; } } else { d=xm; e=d; } a=b; fa=fb; if (Math.abs(d) > tol1) b += d; else b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1)); fb=Ps(b)-Pi; } System.out.println("Maximum number of iterations exceeded in algorithm, root might not be valid"); return b; } public double dPsdT(double t) { //derivative of saturation pressure as a function of saturation temperature double MM=rd.aVP[0][9]; double T=t+273.15; double log10=Math.log(10); double pi=0.0; if(t<=rd.th[0]) { pi=(-rd.aVP[0][1]/(T*T)-2.0*rd.aVP[0][2]/(T*T*T)+rd.aVP[0][4]+2.0*rd.aVP[0][5]*T+3.0*rd.aVP[0][6]*T*T)*Math.log(10.0) +rd.aVP[0][3]/T; if(rd.aVP[0][7]!=0.0) { pi+=rd.aVP[0][7]*(-1.0/T-(rd.aVP[0][8]-T)/(T*T))*Math.log((rd.aVP[0][8]-T)*MM)+ (rd.aVP[0][7]*(rd.aVP[0][8]-T)/T)*(-1.0/(rd.aVP[0][8]-T)); } } else if(t>rd.th[0] && t<= rd.th[1]) { pi=(-rd.aVP[1][1]/(T*T)-2.0*rd.aVP[1][2]/(T*T*T)+rd.aVP[1][4]+2.0*rd.aVP[1][5]*T+3.0*rd.aVP[1][6]*T*T)*Math.log(10.0) +rd.aVP[1][3]/T; if(rd.aVP[1][7]!=0.0) { pi+=rd.aVP[1][7]*(-1.0/T-(rd.aVP[1][8]-T)/(T*T))*Math.log((rd.aVP[1][8]-T)*MM)+ (rd.aVP[1][7]*(rd.aVP[1][8]-T)/T)*(-1.0/(rd.aVP[1][8]-T)); } } else if(t>rd.th[1] && t<= rd.th[2]) { pi=(-rd.aVP[2][1]/(T*T)-2.0*rd.aVP[2][2]/(T*T*T)+rd.aVP[2][4]+2.0*rd.aVP[2][5]*T+3.0*rd.aVP[2][6]*T*T)*Math.log(10.0) +rd.aVP[2][3]/T; if(rd.aVP[2][7]!=0.0) { pi+=rd.aVP[2][7]*(-1.0/T-(rd.aVP[2][8]-T)/(T*T))*Math.log((rd.aVP[2][8]-T)*MM)+ (rd.aVP[2][7]*(rd.aVP[2][8]-T)/T)*(-1.0/(rd.aVP[2][8]-T)); } } else if(t>rd.th[2] && t<= rd.th[3]) { pi=(-rd.aVP[3][1]/(T*T)-2.0*rd.aVP[3][2]/(T*T*T)+rd.aVP[3][4]+2.0*rd.aVP[3][5]*T+3.0*rd.aVP[3][6]*T*T)*Math.log(10.0) +rd.aVP[3][3]/T; if(rd.aVP[3][7]!=0.0) { pi+=rd.aVP[3][7]*(-1.0/T-(rd.aVP[3][8]-T)/(T*T))*Math.log((rd.aVP[3][8]-T)*MM)+ (rd.aVP[3][7]*(rd.aVP[3][8]-T)/T)*(-1.0/(rd.aVP[3][8]-T)); } } else if(t>rd.th[3] && t<= rd.th[4]) { pi=(-rd.aVP[4][1]/(T*T)-2.0*rd.aVP[4][2]/(T*T*T)+rd.aVP[4][4]+2.0*rd.aVP[4][5]*T+3.0*rd.aVP[4][6]*T*T)*Math.log(10.0) +rd.aVP[4][3]/T; if(rd.aVP[4][7]!=0.0) { pi+=rd.aVP[4][7]*(-1.0/T-(rd.aVP[4][8]-T)/(T*T))*Math.log((rd.aVP[4][8]-T)*MM)+ (rd.aVP[4][7]*(rd.aVP[4][8]-T)/T)*(-1.0/(rd.aVP[4][8]-T)); } } else if(t> rd.th[4]) { pi=(-rd.aVP[5][1]/(T*T)-2.0*rd.aVP[5][2]/(T*T*T)+rd.aVP[5][4]+2.0*rd.aVP[5][5]*T+3.0*rd.aVP[5][6]*T*T)*Math.log(10.0) +rd.aVP[5][3]/T; if(rd.aVP[5][7]!=0.0) { pi+=rd.aVP[5][7]*(-1.0/T-(rd.aVP[5][8]-T)/(T*T))*Math.log((rd.aVP[5][8]-T)*MM)+ (rd.aVP[5][7]*(rd.aVP[5][8]-T)/T)*(-1.0/(rd.aVP[5][8]-T)); } } pi*=Ps(t); return pi; //MPa } public double rol(double t) {//saturated liquid density double T=t+273.15; if( T >= rd.Tc) return rd.M/rd.Vc; double roli=0; if(rd.refrigerantName.equals("R11") || rd.refrigerantName.equals("R14") || rd.refrigerantName.equals("R22") || rd.refrigerantName.equals("R23") || rd.refrigerantName.equals("R113") || rd.refrigerantName.equals("R115") || rd.refrigerantName.equals("RC318") || rd.refrigerantName.equals("R500") || rd.refrigerantName.equals("R502") || rd.refrigerantName.equals("R123") || rd.refrigerantName.equals("R134a") || rd.refrigerantName.equals("R503")) { roli=rd.E[0]; double to=(rd.Tc-T)/rd.Tc; for(int i=1;i<5;i++) {roli+=rd.E[i]*Math.pow(to,((double)i/3.0));} } else if (rd.refrigerantName.equals("R12") || rd.refrigerantName.equals("R13") || rd.refrigerantName.equals("R13B1")) { double to=(rd.Tc-T); roli=rd.E[0]+rd.E[1]*to+rd.E[2]*Math.sqrt(to)+rd.E[3]*Math.pow(to,(1.0/3.0))+ rd.E[4]*to*to; } else if (rd.refrigerantName.equals("R21")) { roli=rd.E[0]+rd.E[1]*T+rd.E[2]*T*T; } return roli*rd.M; } public double rog_app(double t) { if( (t+273.15) >= rd.Tc) return rd.M/rd.Vc; //saturated vapor density (approximation) double T=t+273.15; double roc=1.0/rd.Vc; double to=(rd.Tc-T)/rd.Tc; double ro=0; if(rd.refrigerantName.equals("R11") || rd.refrigerantName.equals("R14") || rd.refrigerantName.equals("R22") || rd.refrigerantName.equals("R23") || rd.refrigerantName.equals("R113") || rd.refrigerantName.equals("R115") || rd.refrigerantName.equals("RC318") || rd.refrigerantName.equals("R500") || rd.refrigerantName.equals("R502") || rd.refrigerantName.equals("R503") || rd.refrigerantName.equals("R12") || rd.refrigerantName.equals("R13") || rd.refrigerantName.equals("R13B1") || rd.refrigerantName.equals("R21")) { ro=rd.F[24]*Math.log(T/rd.Tc); for(int i=-10;i<14;i++) { ro+=rd.F[i+10]*Math.pow(to,(double)(i/3.0)); } ro=Math.exp(ro)/rd.Vc*rd.M; } else if(rd.refrigerantName.equals("R123") || rd.refrigerantName.equals("R134a")) { ro=rd.F[0]; double car=t; for(int i=1;i<25;i++) { ro+=rd.F[i]*car;car*=t;} ro=Math.pow(10,ro); } return ro; } public double vg_app(double t) { //saturated vapor specific volume (approximation) return 1.0/rog_app(t); } public double rog(double t) {return 1.0/vsg_t(t);} public double vsg_t(double t) { if( (t+273.15) >= rd.Tc) return rd.Vc/rd.M; // x in f(x)=0 given also df(x)/dx // Newton - Raphson type method //return v_tp(t,Ps(t)); double Pi=Ps(t); double T=t+273.15; double x=x=vg_app(t); int nmax=100; double tolerance=1.0e-2; double fx,dfx; for(int i=0;i