Commit 00df5bf6 authored by Lluis Jofre Cruanyes's avatar Lluis Jofre Cruanyes
Browse files

Entropy calculation added

parent 0dea4277

src/ThermodynamicModel.cpp

+93 −0
Original line number Diff line number Diff line
@@ -104,6 +104,19 @@ double IdealGasModel::calculateInternalEnergyFromPressureTemperatureDensity(cons 


};



double IdealGasModel::calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho) {



    /// Ideal-gas model:

    /// s = c_v*ln( P/rho^gamma ) is specific entropy



    double c_v = R_specific/( gamma - 1.0 );



    double s = c_v*log( P/pow( rho, gamma ) );



    return( s );



};



void IdealGasModel::calculatePressureTemperatureFromDensityInternalEnergy(double &P, double &T, const double &rho, const double &e) {



    /// Ideal-gas model:
@@ -312,6 +325,19 @@ double StiffenedGasModel::calculateInternalEnergyFromPressureTemperatureDensity( 


};



double StiffenedGasModel::calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho) {



    /// Stiffened-gas model:

    /// s = c_v*ln( ((P+P_inf)/rho^gamma ) is specific entropy



    double c_v = R_specific/( gamma - 1.0 );



    double s = c_v*log( ( P+P_inf )/pow( rho, gamma ) );



    return( s );



};



void StiffenedGasModel::calculatePressureTemperatureFromDensityInternalEnergy(double &P, double &T, const double &rho, const double &e) {



    /// Stiffened-gas model:
@@ -590,6 +616,20 @@ double PengRobinsonModel::calculateInternalEnergyFromPressureTemperatureDensity( 


};



double PengRobinsonModel::calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho) {



    /// Peng-Robinson model:

    /// D. Y. Peng, D. B. Robinson.

    /// A new two-constant equation of state.

    /// Industrial and Engineering Chemistry: Fundamentals, 15, 59-64, 1976.



    double bar_v = molecular_weight/rho;

    double s     = ( 1.0/molecular_weight )*this->calculateMolarEntropyFromPressureTemperatureMolarVolume( P, T, bar_v );



    return( s );



};



void PengRobinsonModel::calculatePressureTemperatureFromDensityInternalEnergy(double &P, double &T, const double &rho, const double &e) {



    //double P_norm   = this->critical_pressure;	/// Set pressure normalization factor
@@ -735,6 +775,14 @@ double PengRobinsonModel::calculateMolarInternalEnergyFromPressureTemperatureMol 


};



double PengRobinsonModel::calculateMolarEntropyFromPressureTemperatureMolarVolume(const double &P, const double &T, const double &bar_v) {



    double bar_s = this->calculateMolarStdEntropyFromNASApolynomials( T ) + this->calculateDepartureFunctionMolarEntropy( P, T, bar_v );



    return( bar_s );



};



double PengRobinsonModel::calculate_eos_a(const double &T) {



    /// Peng-Robinson model:
@@ -872,6 +920,34 @@ double PengRobinsonModel::calculateMolarStdEnthalpyFromNASApolynomials(const dou 


};



double PengRobinsonModel::calculateMolarStdEntropyFromNASApolynomials(const double &T) {



    double std_bar_s = 0.0;



    if( (T >= 200.0 ) && ( T < 1000.0 ) ) {

	std_bar_s = R_universal*( NASA_coefficients[7]*log(T) + NASA_coefficients[8]*T + NASA_coefficients[9]*pow( T, 2.0 )/2.0 + NASA_coefficients[10]*pow( T, 3.0 )/3.0 + NASA_coefficients[11]*pow( T, 4.0 )/4.0 + NASA_coefficients[13] );

    } else if( ( T >= 1000.0 ) && ( T < 6000.0 ) ) {

	std_bar_s = R_universal*( NASA_coefficients[0]*log(T) + NASA_coefficients[1]*T + NASA_coefficients[2]*pow( T, 2.0 )/2.0 + NASA_coefficients[3]*pow( T, 3.0 )/3.0 + NASA_coefficients[4]*pow( T, 4.0 )/5.0 + NASA_coefficients[6] );

    } else if( T < 200.0 ) {

	// Assume linear interpolation from T = 200 K 

        double T_min = 200.0;

	    

	double std_bar_s_min   = R_universal*( NASA_coefficients[7]*log(T_min) + NASA_coefficients[8]*T_min + NASA_coefficients[9]*pow( T_min, 2.0 )/2.0 + NASA_coefficients[10]*pow( T_min, 3.0 )/3.0 + NASA_coefficients[11]*pow( T_min, 4.0 )/4.0 + NASA_coefficients[13] );

	double std_bar_s_slope = R_universal*( NASA_coefficients[7]/T_min + NASA_coefficients[8] + NASA_coefficients[9]*T_min + NASA_coefficients[10]*pow( T_min, 2.0 ) + NASA_coefficients[11]*pow( T_min, 3.0 ) );



	std_bar_s = std_bar_s_min + std_bar_s_slope*( T - T_min );

    } else {

#if _ACTIVATE_COUT_

	cout << endl << "NASA 7-coefficient polynomials for std bar s. T = " << T << " is above 6000 K." << endl << endl;

        MPI_Abort( MPI_COMM_WORLD, 1 );										

#endif

        int* trash = nullptr; *trash = 42;	// Causes segmentation fault

    }



    return( std_bar_s );



};



double PengRobinsonModel::calculateDepartureFunctionMolarCp(const double &P, const double &T, const double &bar_v) {



    /// Peng-Robinson model:
@@ -928,6 +1004,23 @@ double PengRobinsonModel::calculateDepartureFunctionMolarEnthalpy(const double & 
  

};



double PengRobinsonModel::calculateDepartureFunctionMolarEntropy(const double &P, const double &T, const double &bar_v) {



    /// Peng-Robinson model:

    /// D. Y. Peng, D. B. Robinson.

    /// A new two-constant equation of state.

    /// Industrial and Engineering Chemistry: Fundamentals, 15, 59-64, 1976.



    double Z                      = this->calculate_Z( P, T, bar_v );

    double A                      = this->calculate_A( P, T );

    double B                      = this->calculate_B( P, T );



    double Delta_bar_s = R_universal*( log( Z - B ) + ( A/( 2.0*sqrt( 2.0 )*B ) )*( eos_kappa*sqrt( T/critical_temperature )/( 1.0 + eos_kappa*( 1.0 - sqrt( T/critical_temperature ) ) ) )*log( ( Z + ( 1.0 - sqrt( 2.0 ) )*B )/( Z + ( 1.0 + sqrt( 2.0 ) )*B ) ) );



    return( Delta_bar_s );

  

};



double PengRobinsonModel::calculateTemperatureFromPressureMolarVolume(const double &P, const double &bar_v) {



    /// Numerical Recipes in C++, Second Edition.

src/ThermodynamicModel.hpp

+25 −0
Original line number Diff line number Diff line
@@ -64,6 +64,11 @@ class BaseThermodynamicModel { 
	#pragma acc routine

        double calculateInternalEnergyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho) { return 0.0; };	// ... modified for OpenACC



        /// Calculate entropy from pressure, temperature and density

        //virtual double calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho) = 0;

	#pragma acc routine

        double calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho) { return 0.0; };	// ... modified for OpenACC



        /// Calculate pressure and temperature from density and internal energy

        //virtual void calculatePressureTemperatureFromDensityInternalEnergy(double &P, double &T, const double &rho, const double &e) {};

	#pragma acc routine
@@ -157,6 +162,10 @@ class IdealGasModel : public BaseThermodynamicModel { 
	#pragma acc routine

        double calculateInternalEnergyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho);



        /// Calculate entropy from pressure, temperature and density

	#pragma acc routine

        double calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho);



        /// Calculate pressure and temperature from density and internal energy

	#pragma acc routine

        void calculatePressureTemperatureFromDensityInternalEnergy(double &P, double &T, const double &rho, const double &e);
@@ -236,6 +245,10 @@ class StiffenedGasModel : public BaseThermodynamicModel { 
	#pragma acc routine

        double calculateInternalEnergyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho);



        /// Calculate entropy from pressure, temperature and density

	#pragma acc routine

        double calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho);



        /// Calculate pressure and temperature from density and internal energy

	#pragma acc routine

        void calculatePressureTemperatureFromDensityInternalEnergy(double &P, double &T, const double &rho, const double &e);
@@ -322,6 +335,10 @@ class PengRobinsonModel : public BaseThermodynamicModel { 
	#pragma acc routine

        double calculateInternalEnergyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho);



        /// Calculate entropy from pressure, temperature and density

	#pragma acc routine

        double calculateEntropyFromPressureTemperatureDensity(const double &P, const double &T, const double &rho);



        /// Calculate pressure and temperature from density and internal energy

	#pragma acc routine

        void calculatePressureTemperatureFromDensityInternalEnergy(double &P, double &T, const double &rho, const double &e);
@@ -358,6 +375,10 @@ class PengRobinsonModel : public BaseThermodynamicModel { 
	#pragma acc routine

        double calculateMolarInternalEnergyFromPressureTemperatureMolarVolume(const double &P, const double &T, const double &bar_v);



        /// Calculate molar entropy from pressure, temperature and molar volume

	#pragma acc routine

        double calculateMolarEntropyFromPressureTemperatureMolarVolume(const double &P, const double &T, const double &bar_v);



        /// Calculate attractive-forces a coefficient

	#pragma acc routine

        double calculate_eos_a(const double &T);
@@ -389,6 +410,8 @@ class PengRobinsonModel : public BaseThermodynamicModel { 
        double calculateMolarStdCpFromNASApolynomials(const double &T);

	#pragma acc routine

        double calculateMolarStdEnthalpyFromNASApolynomials(const double &T);

	#pragma acc routine

        double calculateMolarStdEntropyFromNASApolynomials(const double &T);



        /// Calculate high-pressure departure functions

	#pragma acc routine
@@ -397,6 +420,8 @@ class PengRobinsonModel : public BaseThermodynamicModel { 
        double calculateDepartureFunctionMolarCv(const double &P, const double &T, const double &bar_v);

	#pragma acc routine

        double calculateDepartureFunctionMolarEnthalpy(const double &P, const double &T, const double &bar_v);

	#pragma acc routine

        double calculateDepartureFunctionMolarEntropy(const double &P, const double &T, const double &bar_v);



        /// Calculate temperature from pressure and molar volume

	#pragma acc routine

stuff/rhea_thermodynamics_transport_coefficients.py

+84 −14
Original line number Diff line number Diff line
@@ -39,6 +39,11 @@ class BaseThermodynamicModel: 
    return 0.0





  def calculateEntropyFromPressureTemperatureDensity(self, P, T, rho):



    return 0.0





  def calculatePressureTemperatureFromDensityInternalEnergy(self, P, T, rho, e):



    return 0.0
@@ -120,15 +125,26 @@ class IdealGasModel(BaseThermodynamicModel): 


  def calculateInternalEnergyFromPressureTemperatureDensity(self, P, T, rho):



    #Equation of Specific heat at constant volume

    # Specific heat at constant volume

    c_v = self.R_specific/(self.gamma - 1.0)



    #Equation of Internal Energy 

    # Specific internal energy 

    e = c_v*T



    return e





  def calculateEntropyFromPressureTemperatureDensity(self, P, T, rho):



    # Specific heat at constant volume

    c_v = self.R_specific/(self.gamma - 1.0)



    # Specific entropy 

    s = c_v*np.log( P/( rho**self.gamma ) );



    return s





  def calculatePressureTemperatureFromDensityInternalEnergy(self, P, T, rho, e):



    c_v = self.R_specific/(self.gamma - 1.0)
@@ -315,6 +331,14 @@ class PengRobinsonModel(BaseThermodynamicModel): 
    return e





  def calculateEntropyFromPressureTemperatureDensity(self, P, T, rho):



    bar_v = self.molecular_weight/rho

    s = (1.0/self.molecular_weight)*self.calculateMolarEntropyFromPressureTemperatureMolarVolume(P, T, bar_v)



    return s





  def calculatePressureTemperatureFromDensityInternalEnergy(self, P, T, rho, e):



    # Calculate molar volume
@@ -448,6 +472,13 @@ class PengRobinsonModel(BaseThermodynamicModel): 
    return bar_e

   



  def calculateMolarEntropyFromPressureTemperatureMolarVolume(self,P, T, bar_v):



    bar_s = self.calculateMolarStdEntropyFromNASApolynomials(T) + self.calculateDepartureFunctionMolarEntropy(P, T, bar_v)

    

    return bar_s





  def calculate_eos_a( self, T ):

     

     eos_a = (0.457*((self.R_universal*self.critical_temperature)**2)/(self.critical_pressure))*(1+self.eos_kappa*(1-np.sqrt(T/self.critical_temperature)))**2
@@ -517,47 +548,63 @@ class PengRobinsonModel(BaseThermodynamicModel): 




  def calculateMolarStdCpFromNASApolynomials(self, T):



    std_bar_c_p = 0.0

    if 200.0 <= T < 1000.0:

        std_bar_c_p = self.R_universal*(self.NASA_coefficients[7] + self.NASA_coefficients[8]*T + self.NASA_coefficients[9]*T**2.0 + self.NASA_coefficients[10]*T**3.0 + self.NASA_coefficients[11]* T**4.0)

        std_bar_c_p = self.R_universal*( self.NASA_coefficients[7] + self.NASA_coefficients[8]*T + self.NASA_coefficients[9]*(T**2.0) + self.NASA_coefficients[10]*(T**3.0) + self.NASA_coefficients[11]*(T**4.0) )

    elif 1000.0 <= T < 6000.0:

        std_bar_c_p = self.R_universal*(self.NASA_coefficients[0] + self.NASA_coefficients[1]*T + self.NASA_coefficients[2]*T**2.0 + self.NASA_coefficients[3]*T**3.0 + self.NASA_coefficients[4]* T**4.0)

        std_bar_c_p = self.R_universal*( self.NASA_coefficients[0] + self.NASA_coefficients[1]*T + self.NASA_coefficients[2]*(T**2.0) + self.NASA_coefficients[3]*(T**3.0) + self.NASA_coefficients[4]*(T**4.0) )

    elif T < 200:

        # Assume constant temperature below T = 200 K	    

        T_min = 200.0	    

        std_bar_c_p = self.R_universal*(self.NASA_coefficients[7] + self.NASA_coefficients[8]*T_min + self.NASA_coefficients[9]*T_min**2.0 + self.NASA_coefficients[10]*T_min**3.0 + self.NASA_coefficients[11]*T_min**4.0)

        std_bar_c_p = self.R_universal*( self.NASA_coefficients[7] + self.NASA_coefficients[8]*T_min + self.NASA_coefficients[9]*(T_min**2.0) + self.NASA_coefficients[10]*(T_min**3.0) + self.NASA_coefficients[11]*(T_min**4.0) )

    else:

        print(f"\nNASA 7-coefficient polynomials for std bar c_p. T = {T} is above 6000 K.\n\n")

        exit()



    return std_bar_c_p

     



  def calculateMolarStdEnthalpyFromNASApolynomials(self, T):

    std_bar_h = 0.0



    std_bar_h = 0.0

    if T >= 200.0 and T < 1000.0:

      #std_bar_h = R_universal*T*( NASA_coefficients[7] + NASA_coefficients[8]*T/2.0 + NASA_coefficients[9]*pow( T, 2.0 )/3.0 + NASA_coefficients[10]*pow( T, 3.0 )/4.0 + NASA_coefficients[11]*pow( T, 4.0 )/5.0 + NASA_coefficients[12]/T) - R_universal*NASA_coefficients[14];

      std_bar_h = self.R_universal*T*( self.NASA_coefficients[7] + self.NASA_coefficients[8]*T/2.0 + self.NASA_coefficients[9]*(T**2.0)/3.0 + self.NASA_coefficients[10]*(T**3.0)/4.0 + self.NASA_coefficients[11]*(T**4.0)/5.0 + self.NASA_coefficients[12]/T )

    elif T >= 1000.0 and T < 6000.0:

      #//std_bar_h = R_universal*T*( NASA_coefficients[0] + NASA_coefficients[1]*T/2.0 + NASA_coefficients[2]*pow( T, 2.0 )/3.0 + NASA_coefficients[3]*pow( T, 3.0 )/4.0 + NASA_coefficients[4]*pow( T, 4.0 )/5.0 + NASA_coefficients[5]/T ) - R_universal*NASA_coefficients[14];

      std_bar_h = self.R_universal*T*( self.NASA_coefficients[0] + self.NASA_coefficients[1]*T/2.0 + self.NASA_coefficients[2]*(T**2.0)/3.0 + self.NASA_coefficients[3]*(T**3.0)/4.0 + self.NASA_coefficients[4]*(T**4.0)/5.0 + self.NASA_coefficients[5]/T )

    elif T < 200.0:

      T_min = 200.0



      #std_bar_h_min   = R_universal*T_min*( NASA_coefficients[7] + NASA_coefficients[8]*T_min/2.0 + NASA_coefficients[9]*pow( T_min, 2.0 )/3.0 + NASA_coefficients[10]*pow( T_min, 3.0 )/4.0 + NASA_coefficients[11]*pow( T_min, 4.0 )/5.0 + NASA_coefficients[12]/T_min ) - R_universal*NASA_coefficients[14];

      std_bar_h_min   = self.R_universal*T_min*( self.NASA_coefficients[7] + self.NASA_coefficients[8]*T_min/2.0 + self.NASA_coefficients[9]*(T_min**2.0)/3.0 + self.NASA_coefficients[10]*(T_min**3.0)/4.0 + self.NASA_coefficients[11]*(T_min**4.0)/5.0 + self.NASA_coefficients[12]/T_min )

      std_bar_h_slope = self.R_universal*( self.NASA_coefficients[7] + self.NASA_coefficients[8]*T_min + self.NASA_coefficients[9]*(T_min**2.0) + self.NASA_coefficients[10]*(T_min**3.0) + self.NASA_coefficients[11]*T_min**4.0)

      std_bar_h_slope = self.R_universal*( self.NASA_coefficients[7] + self.NASA_coefficients[8]*T_min + self.NASA_coefficients[9]*(T_min**2.0) + self.NASA_coefficients[10]*(T_min**3.0) + self.NASA_coefficients[11]*(T_min**4.0) )

      std_bar_h = std_bar_h_min + std_bar_h_slope*(T-T_min)

      

    else:

       

       print(f"\nNASA 7-coefficient polynomials for std bar c_p. T = {T} is above 6000 K.\n\n")

       print(f"\nNASA 7-coefficient polynomials for std bar h. T = {T} is above 6000 K.\n\n")

       exit()



    return std_bar_h





  def calculateMolarStdEntropyFromNASApolynomials(self, T):



    std_bar_s = 0.0

    if T >= 200.0 and T < 1000.0:

      std_bar_s = self.R_universal*( self.NASA_coefficients[7]*np.log(T) + self.NASA_coefficients[8]*T + self.NASA_coefficients[9]*(T**2.0)/2.0 + self.NASA_coefficients[10]*(T**3.0)/3.0 + self.NASA_coefficients[11]*(T**4.0)/4.0 + self.NASA_coefficients[13] )

    elif T >= 1000.0 and T < 6000.0:

      std_bar_s = self.R_universal*( self.NASA_coefficients[0]*np.log(T) + self.NASA_coefficients[1]*T + self.NASA_coefficients[2]*(T**2.0)/2.0 + self.NASA_coefficients[3]*(T**3.0)/3.0 + self.NASA_coefficients[4]*(T**4.0)/4.0 + self.NASA_coefficients[6] )

    elif T < 200.0:

      T_min = 200.0

      std_bar_s_min   = self.R_universal*( self.NASA_coefficients[7]*np.log(T_min) + self.NASA_coefficients[8]*T_min + self.NASA_coefficients[9]*(T_min**2.0)/2.0 + self.NASA_coefficients[10]*(T_min**3.0)/3.0 + self.NASA_coefficients[11]*(T_min**4.0)/4.0 + self.NASA_coefficients[13] )

      std_bar_s_slope = self.R_universal*( self.NASA_coefficients[7]/T_min + self.NASA_coefficients[8] + self.NASA_coefficients[9]*T_min + self.NASA_coefficients[10]*(T_min**2.0) + self.NASA_coefficients[11]*(T_min**3.0) )

      std_bar_s = std_bar_s_min + std_bar_s_slope*(T-T_min)

    else:

       print(f"\nNASA 7-coefficient polynomials for std bar s. T = {T} is above 6000 K.\n\n")

       exit()



    return std_bar_s





  def calculateDepartureFunctionMolarCp(self, P, T, bar_v):



    # Peng-Robinson model:

    # D.Y. Peng, D. B. Robinson 

    # A new two-constants equation of state
@@ -578,6 +625,11 @@ class PengRobinsonModel(BaseThermodynamicModel): 


  def calculateDepartureFunctionMolarCv(self, P, T, bar_v):



    # Peng-Robinson model 

    # D. Y. Peng, D. B. Robinson 

    # A new two-constant equations of State

    # Industrial and engineering Chemistry: Fundamental, 15, 59-64 , 1976.



    eos_a_second_derivative = self.calculate_eos_a_second_derivative( T )

    Z                       = self.calculate_Z(P, T, bar_v)

    B                       = self.calculate_B(P, T)
@@ -588,10 +640,12 @@ class PengRobinsonModel(BaseThermodynamicModel): 




  def calculateDepartureFunctionMolarEnthalpy(self, P, T, bar_v):



    # Peng-Robinson model 

    # D. Y. Peng, D. B. Robinson 

    # A new two-constant equations of State

    # Industrial and engineering Chemistry: Fundamental, 15, 59-64 , 1976.



    eos_a                  = self.calculate_eos_a( T )

    eos_a_first_derivative = self.calculate_eos_a_first_derivative( T )

    Z                      = self.calculate_Z(P, T, bar_v)
@@ -602,6 +656,22 @@ class PengRobinsonModel(BaseThermodynamicModel): 
    return Delta_bar_h

 



  def calculateDepartureFunctionMolarEntropy(self, P, T, bar_v):



    # Peng-Robinson model 

    # D. Y. Peng, D. B. Robinson 

    # A new two-constant equations of State

    # Industrial and engineering Chemistry: Fundamental, 15, 59-64 , 1976.



    Z = self.calculate_Z(P, T, bar_v)

    A = self.calculate_A(P, T)

    B = self.calculate_B(P, T)



    Delta_bar_s = self.R_universal*( np.log(Z - B) + ( A/( 2.0*np.sqrt( 2.0 )*B ) )*( self.eos_kappa*np.sqrt( T/self.critical_temperature )/( 1.0 + self.eos_kappa*( 1.0 - np.sqrt( T/self.critical_temperature ) ) ) )*np.log((Z + (1.0 - np.sqrt(2.0))*B)/(Z + (1.0 + np.sqrt(2.0))*B)) )

    

    return Delta_bar_s





  def calculateTemperatureFromPressureMolarVolume(self, P, bar_v):



    # Initial temperature guess using ideal-gas model