断熱火炎温度を算出する際には、色々な化学種の比熱やエンタルピーの値が欲しくなります。気体の比熱やエンタルピーは、温度の多項式で表現されることが多いです。有名なのは、NASA多項式です。NASA多項式は、比熱を4次方程式で、エンタルピーを5次方程式で表現します。
$$\frac{C_p}{R}=a_1+a_2T +a_3T^2 +a_4T^3 +a_5T^4 $$
$$\frac{H}{RT}=a_1+ \frac{1}{2}T + \frac{1}{3} a_3T^2 + \frac{1}{4} a_4T^3 + \frac{1}{5} a_5T^4 + \frac{a_6}{T}$$
各係数については、Gri-mechの熱物性データ(thermo30.dat)で公開されています。今回は、この熱物性データを読込むプログラムを作成しました。
まずは、thermo30.datを読み込むためのヘッダーです。メンバ関数getThermoDataに化学種の名前を渡すとデータ数の14個のvector<double>が返ってきます。
#ifndef READ_FILE_HPP_INCLUDED
#define READ_FILE_HPP_INCLUDED
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
class ThermoFile {
public:
ThermoFile(const std::string &file_name_) : delim(' ') {
file_name = file_name_;
openFile(file_name);
}
void openFile(const std::string file_name);
void closeFile() { ifs.close(); }
std::vector<std::string> split(const std::string &input);
void DeleteSpace(std::string &buf);
void skip(std::string &str);
void inputData(const std::string &T, std::vector<std::string> &information);
std::vector<double> getThermoData(const std::string &);
private:
std::string file_name;
std::ifstream ifs;
const char delim;
};
void ThermoFile::openFile(const std::string file_name) {
ifs.open(file_name, std::ios::in);
if (ifs.fail()) {
throw std::logic_error("Fail to open " + file_name);
}
}
void ThermoFile::DeleteSpace(std::string &buf) {
size_t pos;
while ((pos = buf.find_first_of(" \t")) != std::string::npos) {
buf.erase(pos, 1);
}
}
std::vector<double> ThermoFile::getThermoData(const std::string &name) {
std::string str;
std::vector<std::string> information;
std::vector<double> data(14);
std::stringstream ss;
while (getline(ifs, str)) {
std::stringstream ss(str);
getline(ss, str, ' ');
if (str == name) {
for (int i = 0; i < 3; i++) {
getline(ifs, str);
DeleteSpace(str);
information.resize(0);
inputData(str, information);
for (long long j = 0, size = information.size(); j < size; j++) {
data.at(5 * i + j) = std::stod(information.at(j));
}
}
break;
}
}
ifs.clear();
ifs.seekg(0, std::ios_base::beg);
return data;
}
void ThermoFile::inputData(const std::string &str,
std::vector<std::string> &information) {
size_t found = 0;
size_t newfound = str.find_first_of("E", found + 1);
while (newfound != std::string::npos) {
if (found == 0) {
information.push_back(str.substr(found, newfound + 4));
} else {
information.push_back(str.substr(found + 4, newfound - found));
}
found = newfound;
newfound = str.find_first_of("E", found + 1);
}
}
#endif
使い方はこんな感じ。CO2の比熱とエンタルピーを求めてみます。
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
#include "readfile.hpp"
static const double UNIVERSAL_GAS_CONSTANT = 1.9872; // unit [cal/K/mol]
double getCp(const double temp, const std::vector<double> &thermalcoefficient) {
double cp = 0.0;
if (temp >= 1000.0) {
for (int i = 0; i < 5; i++) {
cp += thermalcoefficient[i] * pow(temp, i);
}
} else {
for (int i = 7; i < 12; i++) {
cp += thermalcoefficient[i] * pow(temp, i - 7);
}
}
cp *= UNIVERSAL_GAS_CONSTANT;
return cp;
}
double getHp(const double temp, const std::vector<double> &thermalcoefficient) {
double hp = 0.0;
if (temp >= 1000.0) {
for (int i = 0; i < 6; i++) {
if (i != 5) {
hp += thermalcoefficient[i] * pow(temp, i + 1) / (i + 1);
} else if (i == 5) {
hp += thermalcoefficient[5];
}
}
} else {
for (int i = 7; i < 13; i++) {
if (i != 12) {
hp += thermalcoefficient[i] *
pow(temp, i - 6) / (i - 6);
} else if (i == 12) {
hp += thermalcoefficient[i];
}
}
}
hp *= UNIVERSAL_GAS_CONSTANT;
return hp;
}
int main() {
std::string inputthermodata = "./file/thermo30.dat";
ThermoFile thermo(inputthermodata);
std::string name = "CO2";
double temperature;
std::vector<double> a_coef(14);
a_coef = thermo.getThermoData(name);
thermo.closeFile();
std::cout << "cp , hp" << std::endl;
for (temperature = 300; temperature <= 3000; temperature += 100) {
std::cout << getCp(temperature, a_coef) << "," << getHp(temperature, a_coef)
<< std::endl;
}
}上記のプログラムで得られた比熱です。ちゃんと算出できていそう。単位はcal/mol-Kなので注意してください。
下記がエンタルピーです。
これで断熱火炎温度の計算ができそうです。
最後までお読みいただきありがとうございました。
コメント