//===================================================================
// 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 :
//                    circular
//                    rectangular
//                    concentric_annulus
//                    circular_segment
//                    circular_sector
//                    right_triangle
// 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
//          -------------    ---------------------
//          Riveted steel    0.9e-3 - 9e-3
//          Concrete         0.3e-3 - 3e-3
//          Wood stave       0.18e-3 - 0.9 e-3
//          Cast iron        0.26e-3
//          Galvanised iron  0.15e-3
//          Commercial steel 0.045e-3
//          wrought iron     0.045e-3
//          drawn tubing     0.0015e-3
//          Plastic,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.*;

class pipe{
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 "circular","rectangular"...
public String  gasName;  //name of a gas mixture
public Gmix 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 KL;               //additional local pressure drop
// pipe functions
//=============================================

public static double tanh(double x)
{
return (Math.exp(x)-Math.exp(-x))/(Math.exp(x)+Math.exp(-x));
}

public static double f(double Re,double eod)
{
// friction factor for turbulent flow 5000<Re<10^8
double f1=Math.log10(eod/3.7+5.74/ Math.pow(Re,0.9));
f1=1.325/(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.log10(eod/3.7+2.51/Re*X)+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 pipe(String  name,double le)
{
//this function reads pipe configuration
//from file Pipe.txt
l=le;
//beginning of try block
try{
   BufferedReader cfin=new BufferedReader(new FileReader("Pipe.txt"));
int ierror=1;
//************
try{
	while(cfin!=null)
	{
	pipeName=Text.readString(cfin);
	if(pipeName.equals(name)) { ierror=0; break;}
    } //end of while
        }  catch(EOFException e_eof)
     {
     System.out.println("error required pipe "+name+" is not found");
     cfin.close();return;
     }
//************
if(ierror==1)
{System.out.println("error required pipe name cannot be found in the data file\n");return;}
shape=Text.readString(cfin);
gasName=Text.readString(cfin);
g=new Gmix(gasName);
if(g==null)
{
  System.out.println("gmix is not found\n");
}
double pi=4.0*Math.atan(1.0);
if(shape.equals("circular"))
{
d1=Text.readDouble(cfin);
d2=0.0;
A=pi*d1*d1/4;dh=d1;
}
else if(shape.equals("rectangular"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape.equals("concentric_annulus"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
dh=Math.abs(d2-d1);
A=pi*dh*dh/4.0;
}
else if(shape.equals("circular_segment"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));
A=pi*dh*dh/4.0;
}
else if(shape.equals("circular_sector"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
dh=d2/(1+d2)*d1;
A=pi*dh*dh/4.0;
}
else if(shape.equals("right_triangle"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);;
dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
}
eps=0.03e-3;
mass_flow_rate=Text.readDouble(cfin);
P=1.0;;
KL=0.0;
//end of try block
}
   catch(IOException e)
   {
   System.err.println("Error Opening File \n"+e.toString());
   System.exit(1);
   }
}

public pipe(String  name)
{
//this function reads pipe configuration
//from file Pipe.txt
l=1.0;
//beginning of try block
try{
      BufferedReader cfin=new BufferedReader(new FileReader("Pipe.txt"));
int ierror=1;
while(cfin!=null)
{
pipeName=Text.readString(cfin);
if(pipeName.equals(name)) {ierror=0;break;}
}
if(ierror==1)
{System.out.println("error required pipe name cannot be found in the data file\n");return;}
shape=Text.readString(cfin);
gasName=Text.readString(cfin);
g=new Gmix(gasName);
if(g==null)
{
  System.out.println("gmix is not found\n");
}
double pi=4.0*Math.atan(1.0);
if(shape.equals("circular"))
{
d1=Text.readDouble(cfin);
d2=0.0;
A=pi*d1*d1/4;dh=d1;
}
else if(shape.equals("rectangular"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape.equals("concentric_annulus"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
dh=Math.abs(d2-d1);
A=pi*dh*dh/4.0;
}
else if(shape.equals("circular_segment"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));
A=pi*dh*dh/4.0;
}
else if(shape.equals("circular_sector"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);
dh=d2/(1+d2)*d1;
A=pi*dh*dh/4.0;
}
else if(shape.equals("right_triangle"))
{
d1=Text.readDouble(cfin);
d2=Text.readDouble(cfin);;
dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
}
eps=0.03e-3;
mass_flow_rate=Text.readDouble(cfin);
P=1.0;;
KL=0.0;
//end of try block
}
   catch(IOException e)
   {
   System.err.println("Error Opening File \n"+e.toString());
   System.exit(1);
   }

}

void change_d1(double dim1)
{
d1=dim1;
double pi=4.0*Math.atan(1.0);
if(shape=="circular")
{ A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
{
 A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
{dh=Math.abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape.equals("circular_segment"))
{dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
{dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
{
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=="circular")
{ d2=0;A=pi*d1*d1/4;dh=d1;}
else if(shape.equals("rectangular"))
{
 A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape.equals("concentric_annulus"))
{dh=Math.abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape.equals("circular_segment"))
{dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape.equals("circular_sector"))
{dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape.equals("right_triangle"))
{
dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
}


 pipe(String  pn,
	  String  sh,
	  String  gn,
	  double dim1,
      double dim2,
	  double  m,double e,
	  double pp,
      double k,
	  double length)
{
// direct assignment constructor
// this constructor can assign all
// the required values directly without
// going to the data file
pipeName=pn;
l=length;
shape=sh;
    try{
       g=new Gmix(gn);
       } catch(IOException ioe) {System.out.println("IOException");}

double pi= 4.0*Math.atan(1.0);
d1=dim1;
d2=dim2;
if(shape.equals("circular"))
//d1 is diameter of the pipe
{ A=pi*d1*d1/4;dh=d1;d2=0;}
else if(shape.equals("rectangular"))
// d1 and d2 are two sides of the rectangle
{ A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape.equals("concentric_annulus"))
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{ dh=Math.abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape.equals("circular_segment"))
// d1 is diameter and d2 is the angle of the circular segment
{ dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape.equals("circular_sector"))
//d1 is the diameter, d2 is the height of the circular section
{ dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape.equals("right_triangle"))
// d1 and d2 two sides of right triangle
{
dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
eps=e;
mass_flow_rate=m;
P=pp;
KL=k;
}

 pipe(String  pn,
	  String  sh,
	  Gmix gmix1,
	  double dim1,
      double dim2,
	  double m,
      double e,
	  double pp,
      double k,
	  double length)
{
// direct assignment constructor
// this constructor can assign all
// the required valued directly without
// going to the data file
pipeName=pn;
shape=sh;
//setup gas name
    try{
       g=new Gmix(gmix1);
       } catch(IOException ioe) {System.out.println("IOException");}

double pi= 4.0*Math.atan(1.0);
d1=dim1;
d2=dim2;
l=length;
if(shape.equals("circular"))
//d1 is diameter of the pipe
{ A=pi*d1*d1/4;dh=d1;d2=0;}
else if(shape.equals("rectangular"))
// d1 and d2 are two sides of the rectangle
{ A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape.equals("concentric_annulus"))
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{ dh=Math.abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape.equals("circular_segment"))
// d1 is diameter and d2 is the angle of the circular segment
{ dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape.equals("circular_sector"))
//d1 is the diameter, d2 is the height of the circular section
{ dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape.equals("right_triangle"))
// d1 and d2 two sides of right triangle
{ dh=2.0*d1*Math.sin(d2)/(1.0+Math.sin(d2)+Math.cos(d2)); A=pi*dh*dh/4.0;}
eps=e;
mass_flow_rate=m;
P=pp;
KL=k;
}

void  change_gas(String str1)
{
 try{
       g=new Gmix(str1);
    }
    catch(IOException ioe) {System.out.println("IOException");}
}


void  change_pipe(
      String pn,
	  String sh,
	  String gn,
	  double dim1,
      double dim2,
	  double e,
      double m,
	  double pp,
      double k)
{
// change assignment
// this constructor can assign all
// the required valued directly without
// going to the data file
pipeName=pn;
shape=sh;
 try{
       g=new Gmix(gn);
       } catch(IOException ioe) {System.out.println("IOException");}
double pi= 4.0*Math.atan(1.0);
d1=dim1;
d2=dim2;
if(shape.equals("circular"))
//d1 is diameter of the pipe
{ A=pi*d1*d1/4;dh=d1;}
else if(shape.equals("rectangular"))
// d1 and d2 are two sides of the rectangle
{ A=d1*d2;dh=2.0*d1*d2/(d1+d2);
}
else if(shape.equals("concentric_annulus"))
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{ dh=Math.abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape.equals("circular_segment"))
// d1 is diameter and d2 is the angle of the circular segment
{ dh=d1*(1+Math.sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape.equals("circular_sector"))
//d1 is the diameter, d2 is the height of the circular section
{ dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape.equals("right_triangle"))
// d1 and d2 two sides of right triangle
{ dh=4.0*d1*d2/(d1+d2+ Math.pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
eps=e;
mass_flow_rate=m;
P=pp;
KL=k;
}


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

public double  get_d2()
//getting the first dimension
{
return d2;
}

public double  velocity(double t)
{
// velocity of the pipe m/x
double d =g.M/g.v(t,P);
return (double)(mass_flow_rate/d/A);
}

public double  Re(double t)
{
// reynold's number
double d=g.M/g.v(t,1.0);
double vis=g.vis(t);
double vv=velocity(t);
double res=(dh*vv*d/vis);
return res;
}

public double  f(double t)
{
// this function calculated friction coefficient
// for both laminar and turnulent flow for various
// cross sections
double Rey=Re(t); //Reynold's number
double C=64;
if(Rey<2000)
  {
  String c         = "circular";
  String r         = "rectangular";
  String co_an     = "concentric_annulus";
  String cir_seg   = "circular_segment";
  String cir_sec   = "circular_sector";
  String rig_tri;
  rig_tri=new String("right_triangle");

  if(shape.equals(c))
     //circular cross section, laminar flow
     { C=64.0; }
  else if(shape.equals(r))
     {
     // Laminar rectangular 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))
     {
     //circular 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))
     {
     //circular 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
  {
  double eod=eps/dh;
  double fc=fcol(Rey,eod);
  return fc;
  }
}

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,double P)
{
// pipe pressure drop for the given length plus
// local KL factors
  double V=velocity(t);
  return (KL+f(t)*l/dh)*(g.M/g.v(t,P))*V*V/2.0*1e-5;
}

public double  dP(double t)
{
double P=1;
// pipe pressure drop for the given length plus
// local KL factors
  double V=velocity(t);
  return (KL+f(t)*l/dh)*(g.M/g.v(t,P))*V*V/2.0*1e-5;
}
public double  dPl_e(double t,double P)
{
double Q=mass_flow_rate*g.v(t,P)/g.M;
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+=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 xx=b*a*a*a/(6.0*g.vis(t))*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_P(double Pin)
//change inlet pressure
{ P=Pin; }

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 Rey=Re(t);
double Nu=0;
if(Rey<2000)
  {
  if(shape.equals("circular"))
     { Nu=3.66; }
  else if(shape.equals("rectangular"))
     {
     // Laminar rectangular 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("concentric_annulus"))
     {
      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("circular_segment"))
     {
     Nu=3.66;
     }
  else if(shape.equals("circular_sector"))
     {
     Nu=3.66;
     }
  else if(shape.equals("right_triangle"))
     {
     Nu=2.47;
     }
  }
else
  {
  double fr  = f(t);
  Rey = Re(t);
  double Prt = g.Prandtl(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.k(t);
//  double Prt = g.Pr(t);
  double result;
  result=Nu*kt/dh;
  return result;
}

public double  Nu(double t)
//Nusselt's Number
{
return hc(t)*dh/g.k(t);
}
} //end of class pipe
