weiye
/
Droplet_LTORC
Watch 1
Star 0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
Browse Source

liquid propert of C3/C7 implemented

binary_fuel
Weiye Wang 1 year ago
parent
293673d25d
commit
42cdd427df
12 changed files with 375 additions and 48 deletions
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. BIN
      DropletCombustionTest7
  2. +1
    -1
      Makefile
  3. +40
    -8
      UserData.cpp
  4. +18
    -2
      UserData.h
  5. +3
    -0
      macros.h
  6. +51
    -2
      main.cpp
  7. +59
    -0
      parse.cpp
  8. +5
    -0
      parse.h
  9. +0
    -1
      parse.hpp
  10. +2
    -0
      readme.md
  11. +190
    -34
      residue.cpp
  12. +6
    -0
      residue.h

BIN
DropletCombustionTest7 View File


+ 1
- 1
Makefile View File

@@ -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 =DropletCombustionTest6
EXE =DropletCombustionTest7-binary
#DESTDIR =~/bin
DESTDIR = ../example



+ 40
- 8
UserData.cpp View File

@@ -339,20 +339,29 @@ UserData allocateUserData(FILE *input){
}
}

ier=parseNumber<char>(input, "dropletComposition", MAXBUFLEN, data->dropSpec);
if(ier==-1){
printf("Enter composition of droplet!\n");
return(NULL);
}
//ier=parseNumber<char>(input, "dropletComposition", MAXBUFLEN, data->dropSpec);
//if(ier==-1){
// printf("Enter composition of droplet!\n");
// return(NULL);
//}

ier=parseDrop(input,"dropletComposition",data->dropSpec,data->dropMole,MAXBUFLEN);
ier=parseDrop(input,"dropletDensity",data->dropSpec,data->dropDens,MAXBUFLEN);

ier=parseNumber<int> (input, "nThreads" , MAXBUFLEN, &data->nThreads);

ier=parseNumber<int> (input, "nThreads", MAXBUFLEN, &data->nThreads);
if(data->nThreads<0 ){
printf("nThreads must be greater than 0!\n");
return(NULL);
}

data->nr=0;
ier=parseNumber<double>(input, "PCAD", MAXBUFLEN, &data->PCAD);
ier=parseNumber<double>(input,"RGTC", MAXBUFLEN, &data->RGTC);
ier=parseNumber<int>(input,"JJRG", MAXBUFLEN, &data->JJRG);
ier=parseNumber<double>(input,"deltaT", MAXBUFLEN, &data->deltaT);
data->nr=0;
//data->np=data->nr+1;
data->nt=data->nr+1;
data->ny=data->nt+1;
@@ -361,9 +370,28 @@ UserData allocateUserData(FILE *input){
data->nvar=data->nsp+4; //assign no: of variables (R,T,P,Mdot,nsp species)


double MW[2];
for(int i=0;i<=1;i++){
int speciesIndex = data->gas->speciesIndex(data->dropSpec[i]);
double weight = data->gas->molecularWeight(speciesIndex);
MW[i]= weight;
}

double massSum = 0.0;
for(int i=0;i<=1;i++){
massSum = massSum + data->dropMole[i] * MW[i];
}
data->dropRho = 0.0;
for(int i=0;i<=1;i++){
data->dropMassFrac[i] = data->dropMole[i]*MW[i]/massSum;
data->dropRho = data->dropRho + data->dropMassFrac[i]*data->dropDens[i];
}
printf("Density of droplet is: %.3f [kg/m^3]\n",data->dropRho);
//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)*684.00; //TODO: The density of the droplet(n-heptane) 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
data->massDrop=1.0/3.0*pow(data->Rd,3)*data->dropRho;

if(!data->adaptiveGrid){
data->uniformGrid = new double [data->npts];
@@ -454,5 +482,9 @@ void setSaneDefaults(UserData data){
data->flameTime[0]=0.0e0;
data->flameTime[1]=0.0e0;
data->nTimeSteps=0;
data->PCAD=0.75;
data->RGTC=1.0;
data->JJRG=0;
data->deltaT=200.0;
}


+ 18
- 2
UserData.h View File

@@ -5,6 +5,7 @@
#endif

#include "gridRoutines.h"
#include <string>

#ifndef USER_DEF
#define USER_DEF
@@ -16,7 +17,15 @@ typedef struct UserDataTag{
* species.*/
Cantera::Transport* trmix;
/* Droplet species composition */
char dropSpec[MAXBUFLEN];
//char dropSpec[MAXBUFLEN];
char dropSpec[2][10];
/* Droplet species mole fractions */
double dropMole[2];
/* Droplet species density at given initialTemperature*/
double dropDens[2];
/* Droplet species mass fractions*/
double dropMassFrac[2];
double dropRho;
/*Length of the domain (in meters):*/
double domainLength;
/*Initial Droplet Radius (in meters)*/
@@ -54,7 +63,8 @@ typedef struct UserDataTag{
size_t k_OH;
size_t k_HO2;
/*Species index of droplet composition*/
size_t k_drop;
/*Index starts with 1 instead of 0*/
size_t k_drop[2];
/*User-defined mass flux (kg/m^2/s):*/
double mdot;
/*Flag to solve isobaric/isochoric problem;*/
@@ -193,6 +203,12 @@ typedef struct UserDataTag{
//double* MW;
FILE* timescaleOutput;

/*Following parameters are used for REGRID function*/
double PCAD;
double RGTC;
int JJRG;
double deltaT;


} *UserData;
UserData allocateUserData(FILE *input);


+ 3
- 0
macros.h View File

@@ -15,6 +15,9 @@
* macros here. Also, we define macros to ease referencing various variables in
* the sundials nvector.
*/
#define WORK(i) WORK[i-1]
#define XNEW(i) XNEW[i-1]

#define psi(i) psidata[i-1]

#define T(i) ydata[((i-1)*data->nvar)+data->nt]


+ 51
- 2
main.cpp View File

@@ -213,6 +213,7 @@ int main(){
double dxRatio=dx/dxMin;
int move=0;
int kcur=0;
int RGCOUNT=0;
if(data->adaptiveGrid){
dxMin=data->grid->leastMove;
xOld=maxCurvPosition(ydata, data->nt, data->nvar,
@@ -229,8 +230,6 @@ int main(){
/*Floor small value to zero*/
floorSmallValue(data, &y);
/*reset the tolerance after ignition*/
resetTolerance(data,&y,&atolv);

if(data->quasiSteady){
ier = IDASolve(mem, finalTime, &tNow, y, ydot, IDA_ONE_STEP);
@@ -338,6 +337,56 @@ int main(){

}
}
/*reset the tolerance after ignition*/
//resetTolerance(data,&y,&atolv);
/*regrid and update the solution based on R,re-initialize the problem*/
/*For the time being,we only allow TORC to REGRID once for each run*/
if(data->JJRG ==1 && (maxT >= data->initialTemperature+data->deltaT)){
if(RGCOUNT<1){
RGCOUNT = RGCOUNT +1;
REGRID(ydata,ydotdata,data);
initializePsiGrid(ydata,data->uniformGrid,data);
printf("REGRID Complete!Restarting Problem at %15.6e s\n",tNow);
ier = IDAReInit(mem,tNow,y,ydot);
if(check_flag(&ier,"IDAReInit",1)){
freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
return(-1);
}

ier= IDASetInitStep(mem,1e-01*dt);
//ier= IDASetInitStep(mem,0.0);
printf("Reinitialized!Calculating Initial Conditions:\n");
printf("Cross your fingers...\n");
ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e-01*dt);
if(check_flag(&ier, "IDACalcIC", 1)){
ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+1e+01*dt);
}
//Every once in a while, for reasons
//that befuddle this mathematically
//lesser author, IDACalcIC fails. Here,
//I desperately try to make it converge
//again by sampling increasingly larger
//time-steps:
for (int i = 0; i < 10; i++) {
ier = IDACalcIC(mem, IDA_YA_YDP_INIT, tNow+(1e-01+pow(10,i)*dt));
if(ier==0){
break;
}
}
//Failure :( Back to the drawing board:
if(check_flag(&ier, "IDACalcIC", 1)){
freeAtLast(mem,&y,&ydot,&res,&id,&atolv,&constraints,data);
return(-1);
}
printf("Initial (Consistent) Conditions Calculated!\n");
printf("------------------------------------------\n\n");

}
}

/*Floor small value to zero*/
floorSmallValue(data, &y);


+ 59
- 0
parse.cpp View File

@@ -18,3 +18,62 @@ void getFromString (const char* buf, char* n){
sscanf(buf,"%s",n);
printf("%s\n",n);
}

/*Extract droplet species and mole fractions*/
int parseDrop(FILE* input, const char* keyword,char dropSpec[][10],double dropMole[],const size_t bufLen){
char buf[bufLen];
char buf1[bufLen];
char comment[1];
char *ret;
while (fgets(buf,bufLen, input)!=NULL){
comment[0]=buf[0];
if(strncmp(comment,"#",1)==0){
}
else{
strcpy(buf1,buf);
ret=strtok(buf,"=");
//DEBUG
//printf("Current KEYWORD in input: %20s \n",ret);
if(strcmp(ret,keyword)==0){
char* modifiedFuel = NULL;
char* equalSign = strstr(buf1,"=");
if(equalSign!= NULL){
modifiedFuel = new char [strlen(equalSign)+1];
strcpy(modifiedFuel,equalSign+1);
//DEBUG
//printf("modifiedFuel:%20s \n",modifiedFuel);
char* token = strtok(modifiedFuel,",");
int index = 0 ;
char* list[2];
while(token!= NULL){
//DEBUG
//printf("TOKEN :%20s \n",token);
list[index] = token;
token = strtok(NULL,",");
index++;
}
//for (int i=0;i<2;i++){
// printf("%20s",list[i]);
//}
for(int i=0;i<2;i++){
char* name= strtok(list[i],":");
char* value = strtok(NULL,":");
//DEBUG
// printf("Name:%10s,Value:%10s \n",name,value);
strcpy(dropSpec[i],name);
dropMole[i]=std::stod(value);
// printf("In the dropArray,Name:%10s,Value:%.3f\n",dropSpec[i],dropMole[i]);
}
delete[] modifiedFuel;
}
printf("%10s:%10s:%.3f,%10s:%.3f\n",keyword,dropSpec[0],dropMole[0],dropSpec[1],dropMole[1]);
//printf("IF statement is execuated. \n");
rewind(input);
//delete[] modifiedFuel;
return(0);
}
}
}
rewind(input);
return(-1);
}

+ 5
- 0
parse.h View File

@@ -3,6 +3,10 @@
#include <string.h> //for strings
#include <stdio.h> //for printf,scanf
#include <stdlib.h> //for atoi, atof
#include <string>
#include <cstring>
#include <iostream>
#include <stdexcept>
#endif

#ifndef PARSE_DEF
@@ -23,6 +27,7 @@ template<typename T>
int parseArray(FILE* input, const char* keyword, const size_t bufLen,
const size_t arrLen, T y[]);

int parseDrop(FILE* input, const char* keyword,char dropSpec[][10],double dropMole[],const size_t bufLen);

#include "parse.hpp"



+ 0
- 1
parse.hpp View File

@@ -73,4 +73,3 @@ int parseArray(FILE* input, const char* keyword, const size_t bufLen,
rewind(input);
return(-1);
}


+ 2
- 0
readme.md View File

@@ -11,4 +11,6 @@ DropletCombustionTest2: Multiple the l function by 100.0 and remove print timeSc
DropletCombustionTest3: Multiple the l function by 10.0 and print the maxTemp and add criteria for regrid
DropletCombustionTest4: dxRatio to be 1.0e-2 and delta_T to be 150 [K]
DropletCombustionTest5: dxRatio to be 1.0 and print x and xOld temperature, modified the fillGrid function;
DropletCombustionTest7: JJRG has been implemented;



+ 190
- 34
residue.cpp View File

@@ -9,6 +9,119 @@
#include <stdio.h>
#include "timing.hpp"

void REGRID(double* ydata,double* ydotdata,UserData data){
size_t nPts = data->npts;
size_t nvar = data->nvar;
double WORK[nPts],XNEW[nPts],XOLD[nPts];
double P,R,R0,R1,R2,TV1,TV2,XLEN,B1,B2;
double DX,DDX;
double DETA,ETAJ,DEL;
size_t ISTART;
P = data->PCAD;
R = data->RGTC;
R0 =1.0 - P;
R1 =P*R/(R+1.0);
R2 =P - R1 ;
TV1 = 0.0;
TV2 = 0.0;

for(size_t i=2;i<=nPts;i++){
TV1 = TV1 + fabs(T(i)-T(i-1));
}
for(size_t i=2;i<=nPts-1;i++){
TV2 = TV2 + fabs( (T(i+1)-T(i))/(R(i+1)-R(i)) - (T(i)-T(i-1))/(R(i)-R(i-1)) );
}
XLEN = R(nPts) - R(1);
B1 = R1* XLEN/ ( (1.0-P)*TV1);
B2 = R2* XLEN/ ( (1.0-P)*TV2);

//Compute partial sum of weight function
WORK(1) =0.0;
DX = 0.0;
for(size_t i=2;i<=nPts-1;i++){
DX = R(i) - R(i-1) ;
WORK(i) = DX + B1*fabs(T(i)-T(i-1)) + WORK(i-1)+B2*fabs( (T(i+1)-T(i))/(R(i+1)-R(i)) - (T(i)-T(i-1))/(R(i)-R(i-1)) );
}
DDX = R(nPts) - R(nPts-1);
WORK(nPts) = WORK(nPts-1) +DDX ;
for(size_t i=2;i<=nPts;i++){
WORK(i) = WORK(i)/WORK(nPts);
}
XNEW(1)=R(1);
XNEW(nPts)=R(nPts);
ISTART =2;
//double DETA;
DETA = 1.0/(nPts-1) ;

//Fill new grid XNEW[nPts],duplicate from PREMIX REGRID SUBROUTINE
for(size_t j=2;j<=nPts-1;j++){
ETAJ= (j-1)*DETA;

for(size_t i=ISTART;i<=nPts;i++){
if(ETAJ <= WORK(i)){
DEL = (ETAJ-WORK(i-1))/(WORK(i)-WORK(i-1)) ;
XNEW(j)=R(i-1)+(R(i)-R(i-1))*DEL;
break;
}else{
ISTART =i;
}

}
}

//Interpolate solution based on XNEW[]&XOLD[]
for(size_t i=1;i<=nPts;i++){
XOLD[i-1]=R(i);
}

INTERPO(ydata,ydotdata,nvar,nPts,XNEW,XOLD);

}


void INTERPO(double* y,double* ydot,const size_t nvar,size_t nPts,const double XNEW[], const double XOLD[]){
double ytemp[nPts],ydottemp[nPts] ;
gsl_interp_accel* acc;
gsl_spline* spline;
acc = gsl_interp_accel_alloc();
spline = gsl_spline_alloc(gsl_interp_steffen,nPts);
gsl_interp_accel* accdot;
gsl_spline* splinedot;
accdot = gsl_interp_accel_alloc();
splinedot = gsl_spline_alloc(gsl_interp_steffen,nPts);

for(size_t j=0;j<nvar;j++){
for(size_t i=0;i<nPts;i++){
ytemp[i]=y[j+i*nvar];
ydottemp[i]=y[j+i*nvar];
}

gsl_spline_init(spline,XOLD,ytemp,nPts);
gsl_spline_init(splinedot,XOLD,ydottemp,nPts);
for(size_t i=0;i<nPts;i++){
y[j+i*nvar]=gsl_spline_eval(spline,XNEW[i],acc);
ydot[j+i*nvar]=gsl_spline_eval(splinedot,XNEW[i],accdot);
}

}

gsl_interp_accel_free(acc);
gsl_spline_free(spline);

gsl_interp_accel_free(accdot);
gsl_spline_free(splinedot);

}

double maxTemperaturePosition(const double* y,const size_t nt,const size_t nvar,const double *x ,size_t nPts){
double maxT = 0.0e0;
double TempT = 0.0e0;
@@ -342,16 +455,24 @@ size_t specIndex(UserData data,char *specName){

int setAlgebraicVariables(N_Vector* id, UserData data){
double *iddata;
char* specPtr[2];

N_VConst(ONE, *id);
iddata = N_VGetArrayPointer_OpenMP(*id);
data->k_bath=BathGasIndex(data);
data->k_oxidizer=oxidizerIndex(data);
data->k_OH=OHIndex(data);
data->k_HO2=HO2Index(data);
data->k_drop=specIndex(data,data->dropSpec);
//use char* pointer to get the address
for(int i=0;i<2;i++){
specPtr[i] = data->dropSpec[i];
}
data->k_drop[0]=specIndex(data,specPtr[0]);
data->k_drop[1]=specIndex(data,specPtr[1]);

printf("Oxidizer index: %lu\n",data->k_oxidizer);
printf("Bath gas index:%lu\n",data->k_bath);
printf("Droplet species index:%lu\n",data->k_drop);
printf("Droplet species index:%lu,%lu\n",data->k_drop[0],data->k_drop[1]);
for (size_t i = 1; i <=data->npts; i++) {
/*Algebraic variables: indicated by ZERO.*/
Rid(i)=ZERO;
@@ -361,7 +482,8 @@ int setAlgebraicVariables(N_Vector* id, UserData data){
}
Mdotid(1)=ZERO;
Rid(1)=ONE;
Yid(1,data->k_drop)=ZERO;
Yid(1,data->k_drop[0])=ZERO;
Yid(1,data->k_drop[1])=ZERO;
Tid(data->npts)=ZERO;
if(data->constantPressure){
Pid(data->npts)=ONE;
@@ -569,8 +691,9 @@ int setInitialCondition(N_Vector* y,

//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*684.0; //TODO:The density of the droplet(n-heptane) should be a user input
if(data->adaptiveGrid){
//massDrop=1.0/3.0*Rd*Rd*Rd*684.0; //TODO:The density of the droplet(n-heptane) should be a user input
massDrop = 1.0/3.0*Rd*Rd*Rd*data->dropRho;
if(data->adaptiveGrid){
psidata = data->grid->xOld;
}
else{
@@ -903,7 +1026,7 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
size_t npts=data->npts;
size_t nsp=data->nsp;
size_t k_bath = data->k_bath;
size_t k_drop = data->k_drop;
size_t k_drop[2] = {data->k_drop[0],data->k_drop[1]};

ydata = N_VGetArrayPointer_OpenMP(y);
ydotdata= N_VGetArrayPointer_OpenMP(ydot);
@@ -952,7 +1075,7 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
double rhomhalf, rhom, lambdamhalf, YVmhalf[nsp],
rho,
rhophalf, lambdaphalf, YVphalf[nsp],
Cpb, Cvb, Cp[nsp], Cpl, dHvl, rhol, wdot[nsp], enthalpy[nsp], energy[nsp],
Cpb, Cvb, Cp[nsp], Cpl[2], dHvl[2], rhol, wdot[nsp], enthalpy[nsp], energy[nsp],
tranTerm, diffTerm, srcTerm, advTerm,
area,areamhalf,areaphalf,aream,areamhalfsq,areaphalfsq;
@@ -982,24 +1105,46 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
//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*/
/*Update specific heat:Cpl for liquid proprane&n-heptane&dodecane*/
/*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
Cpl[0] = 2.423*T(1)+1661.074; //Unit:J/(kg*K),range:1atm,propane
Cpl[1] = 3.336*T(1)+1289.5; //Unit:J/(kg*K),range:1atm,n-Heptane
//Cpl[1] = 3.815*T(1)+1081.8; //Unit:J/(kg*K),range:1atm,n-Dodecane
double Cp_l = 0.0;
for(size_t i=0;i<=1;i++){
Cp_l=Cp_l+Cpl[i]*data->dropMassFrac[i];
}

//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
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
double Tr1 = T(1)/369.8; //Reduced Temperature;
dHvl[0] = 6.32716e5*exp(-0.0208*Tr1)*pow((1-Tr1),0.3766); //Unit:J/kg,Latent Heat of Vaporizaiton of Liquid n-propane,NIST

double Tr2 = T(1)/540.2; //Reduced Temperature:Droplet Temperature divided by critical temperature of liquid n-heptane, Source:NIST
dHvl[1] = 5.366e5*exp(-0.2831*Tr2) * pow((1-Tr2),0.2831); //Unit:J/kg, latent heat of vaporization of liquid n-heptane, Source: NIST
double dHv_l = 0.0;
for(size_t i=0;i<=1;i++){
dHv_l=dHv_l+dHvl[i]*data->dropMassFrac[i];
}

/*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
/*Following section is related to the property of liquid n-heptane and n-dodecane (mainly density rho)*/
/*Density of these two species should be temperature dependent*/
data->dropDens[0] = -1.046*T(1)+823.794; //Unit:kg/m^3,density of propane @1atm
data->dropDens[1] = -0.858*T(1)+933.854; //Unit:kg/m^3,density of n-Heptane @1atm
//data->dropDens[1] = -0.782*T(1)+979.643; //Unit:kg/m^3,density of n-Dodecane @1atm
//rhol = data->dropRho; //Unit:kg/m^3
rhol = 0.0;
for(size_t i=0;i<=1;i++){
rhol = rhol + data->dropMassFrac[i]*data->dropDens[i];
}
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);
@@ -1038,18 +1183,24 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
/*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))) ; //unit:Pa,FOR N-HEPTANE
double p_i[2]={0.0,0.0} ;
p_i[0] = 1.0e5*pow(10,4.53678-(1149.36/(T(1)+24.906)))*data->dropMassFrac[0] ; //unit:Pa,For N-PROPANE
p_i[1] = 1.0e5*pow(10,4.02832-(1268.636/(T(1)-56.199)))*data->dropMassFrac[1] ; //unit:Pa,FOR N-HEPTANE
//p_i = 1.0e3*pow(2.71828,16.7-(4060.0/(T(1)-37.0))); //Unit:Pa,FOR WATER
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);
//Yres(1,k_drop)=Y(1,k_drop) - p_i * data->gas->molecularWeight(k_drop-1)
// / P(1) / data->gas->meanMolecularWeight();
for(size_t i=0;i<=1;i++){
Yres(1,k_drop[i])=Y(1,k_drop[i])-p_i[i]*data->gas->molecularWeight(k_drop[i]-1)/(P(1)*data->gas->meanMolecularWeight());
}

Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass
sum=sum+Y(1,k_drop[0])+Y(1,k_drop[1]);

//Mdotres(1)=Mdot(1) - (YVmhalf(k_drop)*areamhalf) / (1-Y(1,k_drop)); //Evaporating mass
Mdotres(1)= Mdot(1) - ((YVmhalf(k_drop[0])+YVmhalf(k_drop[1])) *areamhalf) / (1-Y(1,k_drop[0])-Y(1,k_drop[1]) ); //Evaporating mass
for (k = 1; k <=nsp; k++) {
if(k!=k_bath && k!=k_drop){
if(k!=k_bath && k!=k_drop[0] && k!=k_drop[1]){
//TODO: May need to include chemical source term
Yres(1,k)=Y(1,k)*Mdot(1) + YVmhalf(k)*areamhalf;
//Yres(1,k)=YVmhalf(k)*areamhalf;
@@ -1076,14 +1227,19 @@ int residue(double t, N_Vector y, N_Vector ydot, N_Vector res, void *user_data){
if(data->dirichletInner){
//Tres(1)=Tdot(1);
//Tres(1)=T(1)-data->innerTemperature;
Tres(1)=lambdamhalf*rhomhalf*areamhalfsq*(T(2)-T(1))/((psi(2)-psi(1))*mass)/dHvl - YVmhalf(k_drop)*areamhalf/(1-Y(1,k_drop));
/*Following code should be revised in the future*/
//Tres(1)=lambdamhalf*rhomhalf*areamhalfsq*(T(2)-T(1))/((psi(2)-psi(1))*mass)/dHv_l - YVmhalf(k_drop)*areamhalf/(1-Y(1,k_drop));
}else{
//Tres(1)=Tdot(1) -
// (rhomhalf*areamhalfsq*lambdamhalf*(T(2)-T(1))/((psi(2)-psi(1))*mass) - Mdot(1)*dHvl)
// / (massDrop * Cpl);
//Tres(1)=Tdot(1) -
// (areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHvl)
// / (massDrop * Cpl);
Tres(1)=Tdot(1) -
(areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHvl)
/ (massDrop * Cpl);
(areamhalf*lambdamhalf*(T(2)-T(1))/(R(2)-R(1)) - Mdot(1)*dHv_l)
/ (massDrop * Cp_l);



@@ -1720,25 +1876,25 @@ void resetTolerance(UserData data, N_Vector* y,N_Vector* atolv)
ydata = N_VGetArrayPointer_OpenMP(*y);
atolvdata = N_VGetArrayPointer_OpenMP(*atolv);
double relTol,radTol,tempTol,presTol,massFracTol,bathGasTol,mdotTol;
double maxT=0.00, ignT=1800;
double maxT=0.00, deltaT=400.0;
relTol = 1.0e-03 ;
radTol = 1.0e-05 ;
tempTol = 1.0e-01 ;
presTol = 1.0e00 ;
massFracTol = 1.0e-08;
bathGasTol = 1.0e-04 ;
mdotTol = 1.0e-09 ;
tempTol = 1.0e-03 ;
presTol = 1.0e-3 ;
massFracTol = 1.0e-06;
bathGasTol = 1.0e-05 ;
mdotTol = 1.0e-10 ;


/*Get the maximum Tempture*/
/*Get the maximum Temperature*/
for (size_t i =1;i <= data->npts;i++){
if(T(i) > maxT){
maxT = T(i);
}
}

/*reset the tolerance when maxT > ignT*/
if(maxT >= ignT){
/*reset the tolerance when maxT > initialTemperature +deltaT*/
if(maxT >= (data->initialTemperature + deltaT)){
data->relativeTolerance = relTol;

for(size_t i =1; i<=data->npts; i++){


+ 6
- 0
residue.h View File

@@ -7,6 +7,8 @@
#ifndef PRINT_DEF
#define PRINT_DEF
#include <string.h> //for strings
#include <string>
#include <cstring>
#include <stdio.h> //for printf,scanf
#include <stdlib.h> //for atoi, atof
#endif
@@ -19,6 +21,10 @@

#include "UserData.h"

void REGRID(double* ydata,double* ydotdata,UserData data);

void INTERPO(double* y,double* ydot,const size_t nvar,size_t nPts,const double XNEW[], const double XOLD[]);

double maxTemperaturePosition(const double* y,const size_t nt,const size_t nvar,const double* x ,size_t nPts);

double maxTemperature(const double* y,const size_t nt, const size_t nvar, size_t nPts);


Write Preview
Loading…
Cancel
Save