//===================================================================
// Thermodynamic package in java
// File Name    : pipe.java
// Author       : Turhan Coban
// EGE University, Engineering Dept. Mechanical Eng. Div.
// email : turhan.coban@mail.ege.edu.tr

// Description  : This file contains the pipe class
//                class pipe calculates pressure drops
//                and internal heat transfer of pipes
// DATA FILE DEFINATION
// pipe datas are written in the data file "Pipe.txt"
// if gas data is not given in the data file, it can be curve fitted
// and added to the data file. Each data has the following form :
//------------------
// pipeName
// shape
// gasName
// d1 d2 (or only d1 depends on shape)
// mass_flow_rate
//-------------------
// a constructor is also available to define this data as passing values
//============================================================
//  VARIABLE IDENTIFICATION
// all the variables that type is not defined is a double variable
// pipeName : String class variable for name
// shape    : String class variable for shape of the cross-section
//            options :
//                    silindir
//                    dikdörtgen
//                    silindirler arası
//                    silindir düzlem
//                    silindir dilimi
//                    dik üçgen
// gasName : String class variable for name of the gas mixtures
// g      : gmix class variable
// mass_flow_rate : mass flow rate of gas mixture flowing
//                  through pipe, kg/s
// eps     : pipe equivalent roughness m
//          -------------------------------------------
//          value of equivalent roughness for some pipes:
//          Pipe             equivalent rougness m
//          -------------    ---------------------
//          sheet steel new      0.05e-3 
//          stainless steel new  0.002e-3
//          commercial steel     0.046e-3
//          Riveted steel        3e-3
//          Rusted steel         2e-3
//          Concrete smooth      0.04e-3
//          Concrete rough       2e-3
//          Wood stave           0.5e-3
//          Cast iron            0.26e-3
//          aphalted cast iron   0.12e-3
//          Galvanised iron      0.15e-3
//          Commercial steel     0.045e-3
//          brass                0.002e-3
//          wrought iron         0.045e-3
//          drawn tubing         0.0015e-3
//          Plastic drawn tubing 0.0015e-3
//          Glass                0.0 (smooth)
//          -------------------------------------------
// P       : Pressure bar
// l       : length m
// dh      : hydraulic diameter m
// A       : pipe area m^2
// KL      : local pressure drop coefficient
//
//============================================================
import java.io.*;
import javax.swing.*;
class pipe3{
protected double d1,d2;  // m pipe dimensions
	                     // actual meaning of d1 and d2
	                     // can be change according to shape
public String  pipeName;
public String  shape;    //pipe shape "silindir","dikdörtgen"...
public String  gasName;  //name of a gas mixture
public liquid_CS g;
public double mass_flow_rate;  // kg/s
public double eps ;            // m
public double P[]  ;             // Pressure bar
public double l ;              // length
public double dh;              // hydraulic diameter
public double A;               // pipe area m^2
public double V[];               //pipe velocity, m/s
public double KL;               //additional local pressure drop
public double z[];

// pipe functions
// =============================================
 
 public static double f(double Re,double eod)
 {
 //Goudar equation 
 double a=2/Math.log(10);
 double b=eod/3.7;
 double d=Math.log(10.0)/5.02*Re;
 double s=b*d+Math.log(d);
 double q=Math.pow(s,(s/(s+1)));
 double g=b*d+Math.log(d/q);
 double z=Math.log(q/g);
 double dLA=(g/(g+1))*z;
 double dCFA=dLA*(1+z/2.0/((g+1)*(g+1)+(z/3)*(2.0*g-1)));
 double f1=a*(Math.log(d/q)+dCFA);
 f1=1.0/(f1*f1);
 //System.out.println("Goudar equation f1="+f1);
 return f1;
 }
public static double fx(double X,double Re,double eod)
{
// colebrook equation to solve
// friction factor for turbulent flow 2000 < Re
double xx=2.0*Math.log(eod/3.7+2.51/Re*X)/Math.log(10.0)+X;
return xx;
}

public  static double dfx(double X,double Re,double eod)
{
//derivative of colebrook equation
double xx;
xx = 1+2.0/(eod/3.7+2.51/Re*X)/Math.log(10.0)*2.51/Re;
return  xx;
}

public  static double fcol(double Re,double eod)
{
//solution of the colebrook equation
// by using newton method
double fi=f(Re,eod);
double x=1.0/ Math.pow(fi,0.5);
int nmax=50;
double tolerance=1.0e-10;
for(int i=0;i<nmax;i++)
{
double fx1=fx(x,Re,eod);
x-=fx1/dfx(x,Re,eod);
if(Math.abs(fx1)<tolerance) return 1.0/(x*x);
}
return 1.0/(x*x);
}

//==============================================
// implementation of constructors functions
// for class pipe

public pipe3(String  name,String liquidname,String shape1,double le,double  d1i,double d2i,double epsi,double P1,double P2,double z1,double z2,double V1,double KL1)
{
	V=new double[2];
	V[0]=V1;
	P=new double[2];
	P[0]=P1;
	P[1]=P2;
	z=new double[2];
	z[0]=z1;
	z[1]=z2;
	pipeName=name;
	gasName=liquidname;
	shape=shape1;
	d1=d1i;
	d2=d2i;
	eps=epsi;
l=le;
//beginning of try block
int ierror=1;
//************

g=new liquid_CS(gasName);
if(g==null)
{
  System.out.println("gmix is not found\n");
}
double pi=4.0*Math.atan(1.0);
if(shape.equals("silindir"))
{
A=pi*d1*d1/4;dh=d1;
}
else if(shape.equals("dikdörtgen"))
{
A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape.equals("silindirler arası"))
{
dh=Math.abs(d2-d1);
A=pi*dh*dh/4.0;
}
else if(shape.equals("silindir düzlem"))
{
dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));
A=pi*dh*dh/4.0;
}
else if(shape.equals("silindir dilimi"))
{
dh=d2/(1+d2)*d1;
A=pi*dh*dh/4.0;
}
else if(shape.equals("dik üçgen"))
{
dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
}
//end of try block

}

void change_d1(double dim1)
{
d1=dim1;
double pi=4.0*Math.atan(1.0);
if(shape=="silindir")
{ A=pi*d1*d1/4;dh=d1;}
else if(shape=="dikdörtgen")
{
 A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="silindirler arası")
{dh=Math.abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape.equals("silindir düzlem"))
{dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="silindir dilimi")
{dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="dik üçgen")
{
dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
}

void  change_d2(double dim2)
{
d2=dim2;
double pi=4.0*Math.atan(1.0);
if(shape=="silindir")
{ d2=0;A=pi*d1*d1/4;dh=d1;}
else if(shape.equals("dikdörtgen"))
{
 A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape.equals("silindirler arası"))
{dh=Math.abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape.equals("silindir düzlem"))
{dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape.equals("silindir dilimi"))
{dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape.equals("dik üçgen"))
{
dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
}

void  change_gas(String str1)
{       g=new liquid_CS(str1);}


public double  get_d1()
//getting the first dimension
{
return d1;
}

public double  get_d2()
//getting the first dimension
{
return d2;
}
public double  fV(double V2,double t)
{   double g=9.81;//m/s^2
    double eod=eps/dh;
    //System.out.println("Z0="+z[0]+"z1="+z[1]);
	return (P[0]-P[1])/2 + g*(z[0]-z[1])+0.5*V[0]*V[0]
	-(0.5*V2*V2+g*(KL+f(Re(V2,t),eod,t))*l/dh*V2*V2/(2*g));

}

public double  Re(double V2,double t)
{
// reynold's number
double T=t+273.0;
double ro=g.rd.ro_l(T);
double vis=g.rd.viscosity_l(T);
double res=(dh*V2*ro/vis);
return res;
}
public double  Re(double t)
{
// reynold's number
double T=t+273.0;
double ro=g.rd.ro_l(T);
double vis=g.rd.viscosity_l(T);
double v=velocity(t);
double res=(dh*v*ro/vis);
//System.out.println("v="+v+"Re="+res);
return res;
}

public double  f(double Re,double eod,double t)
{
// this function calculated friction coefficient
// for both laminar and turbulent flow for various
// cross sections
double Rey=Re; //Reynold's number
double C=64;
if(Rey<2300)
  {
  String c         = "silindir";
  String r         = "dikdörtgen";
  String co_an     = "iki silindir arası";
  String cir_seg   = "silidir düzlem kesiti arası";
  String cir_sec   = "daire dilimi";
  String rig_tri;
  rig_tri=new String("dik üçgen");

  if(shape.equals(c))
     //silindir cross section, laminar flow
     { C=64.0; }
  else if(shape.equals(r))
     {
     // Laminar dikdörtgen channel friction factor
       double x=d1/d2;
       if(x>1) x=1/x;
       double a1[];
       a1=new double[6];
       a1[0]=    .959923253176985E+02;
       a1[1]=   -.130364643234979E+03;
       a1[2]=    .188053177394400E+03;
       a1[3]=  -.166159826967448E+03 ;
       a1[4]=    .952739160606451E+02;
       a1[5]=  -.258949934514987E+02 ;
       C=a1[5];
       for(int i=4;i>=0;i--) { C=C*x+a1[i];}
     }

  else if(shape.equals(co_an))
     {
      //concentric annulus laminar friction factor
      double x=d2/d2;
      if(x>1) x=1/x;
      double a1[];
      a1=new double[5];
      a1[0]=     .717068909524062E+02;
      a1[1]=     .932037998655195E+03;
      a1[2]=    -.947586946520906E+04;
      a1[3]=     .204346570871727E+05;
      a1[4]=    -.118665325115712E+05;
      C=a1[4];
      for(int i=4;i>=0;i--) {C=C*x+a1[i];}
     }
  else if(shape.equals(cir_seg))
     {
     //silindir segment laminar friction factor
      double a1[];
      a1=new double[5];
      a1[0]=     .639999999999997E+02;
      a1[1]=    -.162338014624008E+01;
      a1[2]=     .162620455430327E+01;
      a1[3]=   -.827251633703104E+00 ;
      a1[4]=     .138590724447624E+00;
      C=a1[4];
      for(int i=3;i>=0;i--) {C=C*d2+a1[i];}
     }
  else if(shape.equals(cir_sec))
     {
     //silindir section laminar friction factor
      double a1[];
      a1=new double[4];
      a1[0]=     .480000000000000E+02;
      a1[1]=     .227909935649730E+02;
      a1[2]=    -.134959915493607E+02;
      a1[3]=     .325095868081599E+01;
      C=a1[3];
      for(int i=2;i>=0;i--) {C=C*d2+a1[i];}
     }
  else if(shape.equals(rig_tri))
     {
     //right triangle laminar friction factor
     double a1[];
     a1=new double[5];
      a1[0]=     .480001792918036E+02;
      a1[1]=     .127207559700323E+02;
      a1[2]=    -.108590518340935E+02;
      a1[4]=     .158423382029386E+01;
      a1[5]=     .115425033696128E+01;
      C=a1[5];
      for(int i=4;i>=0;i--) {C=C*d2+a1[i];}
     }
  return C/Rey;
  }
else
  {
  //Goudar equation
  double fc=f(Rey,eod);
  return fc;
  }
}
public double Vmax(double t)
{
double g=9.806;//m/s^2
    double eod=eps/dh;
	double vv=(P[0]-P[1])/2.0 + g*(z[0]-z[1])+0.5*V[0]*V[0];
    vv=Math.sqrt(2*vv);
    return vv;
}

public double velocity(double t)
{
double a=0.001;
double b=Vmax(t);
double test;
double p=0;
double es,ea;
double f1,f2,f3,fp;
int maxit=500,iter=0;
double tol=1.0e-15;
es=0.0000001;
double x1=a;f1=fV(x1,t);
double x2=b;f2=fV(x2,t);
double x3=(x1+x2)/2.0;f3=fV(x3,t);
if(f1==0) return x1;
else if(f2==0) return x2;
else if(f3==0) return x3; 
if(f1*f2>0) System.out.println("No root is existed in the given region");
p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
fp=fV(p,t);
ea=Math.abs(x3-p);
while((ea>es)&&(iter<maxit))
{	if(Math.abs(f3)<tol) return x3; 
    if((p<x1) && (p>x2)) 
    {p=(x1+x2)/2.0;    
     if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
     else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    }
    else
    {        
    if(p>x3) {x1=x3;f1=f3;x3=p;f3=fp;}
    else if(p<x3) {x2=x3;f2=f3;x3=p;f3=fp;}
    p=-(f2*f3*x1*(f2-f3)+f3*f1*x2*(f3-f1)+f1*f2*x3*(f1-f2))/((f1-f2)*(f2-f3)*(f3-f1));
    fp=fV(p,t);
    ea=Math.abs(x3-p);
    }
    iter++;
} 
if(iter>=maxit) JOptionPane.showMessageDialog(null,"Warning : MAximum number of iteration is exceeded \n"+
   " result may not be valid","MAXIMUM ITERATION NUMBER WARNING",JOptionPane.WARNING_MESSAGE);
V[1]=p;
return p;
}

public double  lentry(double t)
{
//hydraulic entry length
double Rey=Re(t);
double le;
if(Rey<2000) le=0.06*dh*Rey;
else le=4.4*dh* Math.pow(Rey,(1.0/6.0));
return le;
}

public double  dP(double t)
{
// pipe pressure drop for the given length plus
// local KL factors
  double T=t+273.0;
  double V=velocity(t);
  double Re=Re(t);
  double ro=g.rd.ro_l(T);
  double eod=eps/dh;
  double ff=f(Re,eod,t);
  return (KL+ff*l/dh)*ro*V*V/2.0*1e-5;
}



public double  dPl_e(double t,double P)
{
	  double T=t+273.0;
 double ro=g.rd.ro_l(T);	
double Q=mass_flow_rate/ro;
double a=d1/2.0;
double b=d2/2.0;
double x=0;
double pi=4.0*Math.atan(1.0);
for(int i=1;i<8000;i+=2)
  {
  x+=Math.tanh(pi*b*i/d1)/(double)(i*i*i*i*i);
  }
x=1.0-192.0*a/(pi*pi*pi*pi*pi*b)*x;
double vis=g.rd.viscosity_l(T);
double xx=b*a*a*a/(6.0*vis)*x;
xx=Q/xx*l*1e-5;
return xx;
}

void  change_length(double l1)
{ l=l1;       }

void  change_mass_flow_rate(double mfr)
{ mass_flow_rate=mfr; }



void  change_KL(double klin)
//change local pressure drop coefficient
{KL=klin; }

void  change_e(double e)
{ //change surface roughness coefficient
eps=e; 
}

public double  hc(double t)
{
//convective heat transfer coefficient
  double T=t+273.0;
double Rey=Re(t);
double Nu=0;
if(Rey<2000)
  {
  if(shape.equals("silindir"))
     { Nu=3.66; }
  else if(shape.equals("dikdörtgen"))
     {
     // Laminar dikdörtgen channel friction factor
       double x=d1/d2;
       if(x>1) x=1/x;
       double a1[];
       a1=new double[7];
       a1[0]=    .754000000000580E+01;
       a1[1]=   -.193789080028494E+02;
       a1[2]=    .333861212302836E+02;
       a1[3]=  -.170496495070800E+02 ;
       a1[4]=   -.304226226912837E+02;
       a1[5]=    .495297442580787E+02;
       a1[6]=   -.206246852871550E+02;
       Nu=a1[6];
      for(int i=5;i>=0;i--) {Nu=Nu*x+a1[i];}
     }
  else if(shape.equals("silindirler arası"))
     {
      double x=d2/d2;
      if(x>1) x=1/x;
      double a1[];
      a1=new double[4];
      a1[0]=     .401698044199492E+01;
      a1[1]=     .921803331830132E+00;
      a1[2]=    -.307845192109755E+00;
      a1[3]=     .229092480046921E+00;
       Nu=a1[3];
      for(int i=2;i>=0;i--) {Nu=Nu*x+a1[i];}
     }
  else if(shape.equals("silindir düzlem"))
     {
     Nu=3.66;
     }
  else if(shape.equals("silindir dilimi"))
     {
     Nu=3.66;
     }
  else if(shape.equals("dik üçgen"))
     {
     Nu=2.47;
     }
  }
else
  {
  double fr  = f(V[1],t);
  Rey = Re(t);
  double Prt = g.rd.Pr_l(T);
  // Gnielinski equation
  // Gnielinski, V, Int. Chem Eng., 16, 359,1976
  // Valid for 2000 < Re <5e5, 0.5 < Pr < 2000
  Nu=fr/8.0*(Rey-1000.0)*Prt/
  (1+12.7* Math.sqrt(fr/8.0)*( Math.pow(Prt,(2.0/3.0))-1.0));
  }
  double kt  =  g.rd.k_l(T);
//  double Prt = g.Pr(t);
  double result;
  result=Nu*kt/dh;
  return result;
}

public double  Nu(double t)
//Nusselt's Number
{
	double T=t+273.0;
return hc(t)*dh/g.rd.k_l(T);
}
public String toString1(double t)
{
String s="";	
String ss[][]=toString(t);
for(int i=0;i<ss.length;i++)
{for(int j=0;j<ss[0].length;j++)
   {s+=ss[i][j];}
}
return s;
}


public String[][] toString(double t)
{    double T=t+273.0;
    String s1[][]=new String[18][3];
    s1[0][0]="Basınç düşüşü ";
    s1[0][1]=""+dP(t);
    s1[0][2]="bar          "+"\n";
    s1[1][0]="f sürtünme katsayısı ";
    double Rey=Re(t);
    double eod=eps/dh;
    s1[1][1]=""+f(Rey,eod);
    s1[1][2]="             "+"\n ";
    s1[2][0]="Akış hızı giriş ";
    s1[2][1]=""+V[0];
    s1[2][2]="m/s"+"\n";   
    s1[3][0]="Akış hızı çıkış ";
    s1[3][1]=""+V[1];
    s1[3][2]="m/s"+"\n";
    s1[4][0]="Akış hızı maksimum çıkış (f=0) ";
    s1[4][1]=""+Vmax(t);
    s1[4][2]="m/s"+"\n"; 
    s1[5][0]="Basınç giriş ";
    s1[5][1]=""+P[0];
    s1[5][2]="Pa "+"\n";   
    s1[6][0]="Basınç çıkış ";
    s1[6][1]=""+P[1];
    s1[6][2]="Pa "+"\n";
    s1[7][0]="yükseklik giriş ";
    s1[7][1]=" "+z[0];
    s1[7][2]=" m "+"\n";   
    s1[8][0]="yükseklik çıkış ";
    s1[8][1]=" "+z[1];
    s1[8][2]=" M "+"\n";       
    s1[9][0]="Reynold sayısı Re ";    
    s1[9][1]=""+Rey;
    s1[9][2]="             "+"\n";
    s1[10][0]="Prandtl sayısı Pr ";
    s1[10][1]=""+g.rd.Pr_l(T);
    s1[10][2]="             "+"\n";
    s1[11][0]="Nusselt sayısı Nu ";
    s1[11][1]=""+Nu(t);
    s1[11][2]="             "+"\n";
    s1[12][0]="Hidrolik giriş gelişme boyu ";
    s1[12][1]=""+lentry(t);
    s1[12][2]="m            "+"\n";
    s1[13][0]=" "+'\u03BC'+" , sıvı viskozitesi ";
    s1[13][1]=""+g.rd.viscosity_l(T);
    s1[13][2]=" Ns/m^2        "+"\n";
    s1[14][0]=" k , sıvı ısıl iletkenliği ";
    s1[14][1]=""+g.rd.k_l(T);
    s1[14][2]=" W/mK        "+"\n";
    s1[15][0]=" Cp, sıvı özgül ısısı ";
    s1[15][1]=""+g.rd.Cp_l(T);
    s1[15][2]=" KJ/kgK        "+"\n";
    s1[16][0]=" "+'\u03C1'+"_l, sıvı yoğunluğu ";
    s1[16][1]=""+g.rd.ro_l(T);
    s1[16][2]=" kg/m^3        "+"\n";
    s1[17][0]="sıvı ismi ";
    s1[17][1]=g.liquidName;
    s1[17][2]="            "+"\n";
    return s1;
}

} //end of class pipe
