Modified the heat capacity of n-heptane

Makefile

@@ -7,7 +7,7 @@
 compiler	=g++
 CANTERA_DIR	=/opt/scientific/cantera-2.4_gnu_blas
 IDA_DIR		=/opt/scientific/sundials-3.1.1_intel_mkl
-EXE		=DropletCombustion2
+EXE		=DropletCombustion3
 #DESTDIR		=~/bin
 DESTDIR		= ../example
 

UserData.cpp

@@ -358,7 +358,8 @@ UserData allocateUserData(FILE *input){
 	data->nvar=data->nsp+4;		 //assign no: of variables (R,T,P,Mdot,nsp species)
 
 	//Mass of droplet
-	data->massDrop=1.0/3.0*pow(data->Rd,3)*997.0; //TODO: The density of the droplet should be a user input
+	//data->massDrop=1.0/3.0*pow(data->Rd,3)*997.0; //TODO: The density of the droplet should be a user input
+    data->massDrop=1.0/3.0*pow(data->Rd,3)*684.00; //TODO: The density of the droplet(n-heptane) should be a user input
 
 	if(!data->adaptiveGrid){
 		data->uniformGrid = new double [data->npts];

main.cpp

@@ -322,6 +322,7 @@ int main(){
 			       "R=%ld ,"
 			       "o=%d ,",netf, ncfn, njevals, nrevals, kcur);
 			printf("Time=%15.6e s,",tNow);
+            printf("dt=%15.6e s,",dt);
 			printf("frac: %15.6e\n",dxRatio);
 			count++;
 			data->nTimeSteps=count;

residue.cpp

@@ -444,8 +444,8 @@ int setInitialCondition(N_Vector* y,
 	massDrop = data->massDrop;
 
 	//Mass of droplet
-	massDrop=1.0/3.0*Rd*Rd*Rd*997.0; //TODO: The density of the droplet should be a user input
-
+	//massDrop=1.0/3.0*Rd*Rd*Rd*997.0; //TODO: The density of the droplet should be a user input
+    massDrop=1.0/3.0*Rd*Rd*Rd*684.0; //TODO:The density of the droplet(n-heptane) should be a user input 
 	if(data->adaptiveGrid){
   		psidata  = data->grid->xOld;
 	}
@@ -820,15 +820,31 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
     Cvb=data->gas->cv_mass();       	//J/kg/K
 	//TODO: Create user input model for these. They should not be constant
 	//Cpl=4.182e03;						//J/kg/K for liquid water
-	Cpl=  5.67508633e-07*pow(T(1),4) - 7.78060597e-04*pow(T(1),3) + 4.08310544e-01*pow(T(1),2)  //J/kg/K for liquid water
-		- 9.62429538e+01*T(1)        + 1.27131046e+04;
-	
+	//Cpl=  5.67508633e-07*pow(T(1),4) - 7.78060597e-04*pow(T(1),3) + 4.08310544e-01*pow(T(1),2)  //J/kg/K for liquid water  
+    //     - 9.62429538e+01*T(1)        + 1.27131046e+04;
+    
+    /*Update specific heat:Cpl for liquid n-heptane*/
+   /*Based on the existiong data,a linear expression is used*/
+   // Cpl= 12.10476*T(1) - 746.60143; // Unit:J/(kg*K), Need to be improved later
+    Cpl = 2242.871642; //Unit:J/(kg*K),range:25~520K,Ginnings&Furukawa,NIST
+
 	//dHvl=2.260e6;						//J/kg	heat of vaporization of water
-	double Tr = T(1)/647.1;				//Reduced Temperature: Droplet Temperature divided by critical temperature of water
-	dHvl= 5.2053e07*pow(1-Tr,0.3199 - 0.212*Tr + 0.25795*Tr*Tr)/18.01;	//J/kg	latent heat of vaporization of water
-	rhol = 997.0;	
-	massDrop=1.0/3.0*R(1)*R(1)*R(1)*997.0; //TODO: The density of the droplet should be a user input
-	aream= calc_area(R(1),&m);	    
+	//double Tr = T(1)/647.1;				//Reduced Temperature: Droplet Temperature divided by critical temperature of water
+	//dHvl= 5.2053e07*pow(1-Tr,0.3199 - 0.212*Tr + 0.25795*Tr*Tr)/18.01;	//J/kg	latent heat of vaporization of water
+	
+    double Tr = T(1)/540.2;            //Reduced Temperature:Droplet Temperature divided by critical temperature of liquid n-heptane, Source:NIST
+    dHvl = 5.366e5*exp(-0.2831*Tr) * pow((1-Tr),0.2831);           //Unit:J/kg, latent heat of vaporization of liquid n-heptane, Source: NIST
+    
+    /*Following section is related to the property of water*/
+    //rhol = 997.0;	
+	//massDrop=1.0/3.0*R(1)*R(1)*R(1)*997.0; //TODO: The density of the droplet should be a user input
+    
+    /*Following section is related to the property of liquid n-heptane (mainly density rho)*/
+    rhol = 684.0;                          //Unit:kg/m^3
+    massDrop = 1.0/3.0*R(1)*R(1)*R(1)*rhol;  //Unit:kg, this is the mass of liquid droplet  
+
+    aream= calc_area(R(1),&m);	    
+    
 
 	/*******************************************************************/
 	/*Calculate values at j=2's m and mhalf*****************************/
@@ -855,12 +871,20 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
 	sum=ZERO;
 	
 	//TODO: Antoine Parameters should be part of user input, not hardcoded
-	Yres(1,k_drop)=Y(1,k_drop) - 1.0e5*pow(10,4.6543-(1435.264/(T(1)-64.848)))*data->gas->molecularWeight(k_drop-1) 
-				    / P(1) / data->gas->meanMolecularWeight();
+	//Yres(1,k_drop)=Y(1,k_drop) - 1.0e5*pow(10,4.6543-(1435.264/(T(1)-64.848)))*data->gas->molecularWeight(k_drop-1) 
+	//			    / P(1) / data->gas->meanMolecularWeight();
 	//Yres(1,k_drop)=Y(1,k_drop) - 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0)))*data->gas->molecularWeight(k_drop-1) 
 	//			    / P(1) / data->gas->meanMolecularWeight();
 	//Yres(1,k_drop)=Y(1,k_drop)-2.566785e-02;
-	sum=sum+Y(1,k_drop);
+	
+    /*Following using the Antoine Parameters to calculate the pressure of volatile component(n-heptane)*/
+    //Antoine parameter from NIST website
+    double p_i=0.0 ; 
+    p_i = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199))) ;
+    Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1) 
+				    / P(1) / data->gas->meanMolecularWeight();
+	
+    sum=sum+Y(1,k_drop);
 
 	Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass