714 lines
21 KiB
C
714 lines
21 KiB
C
/*************************************************************************
|
||
* 文件: therm.c
|
||
* 描述: 包含仪表对象的热力学函数
|
||
* 热容、焓、熵、声速
|
||
* 包含以下函数:
|
||
* ThermInit(), ThermDestroy(), Run(), coth(), CpiMolar(), Ho(), So(),
|
||
* CprCvrHS(), HS_Mode(), H(), S()
|
||
* 版本: ver 1.7 2002.11.17
|
||
* 作者: W.B. Peterson
|
||
* 修订:
|
||
* 版权: (c) 2002 美国天然气协会
|
||
**************************************************************************/
|
||
|
||
#include "NGCal.h"
|
||
#include "Detail.h"
|
||
#include "Therm.h"
|
||
#include <math.h>
|
||
#include <stdbool.h> // 添加bool类型支持
|
||
|
||
// 添加函数前置声明解决隐式声明问题
|
||
static void CprCvrHS(AGA10STRUCT *ptAGA10, Detail *ptD);
|
||
static double H(AGA10STRUCT *ptAGA10, Detail *ptD);
|
||
static double S(AGA10STRUCT *ptAGA10, Detail *ptD);
|
||
static void HS_Mode(AGA10STRUCT *ptAGA10, Detail *ptD, double H_target, double S_target, bool bGuess);
|
||
|
||
double CalTH = 4.1840;
|
||
|
||
// 高斯积分的根和权重
|
||
double GK_root[5] = {
|
||
0.14887433898163121088,
|
||
0.43339539412924719080,
|
||
0.67940956829902440263,
|
||
0.86506336668898451073,
|
||
0.97390652851717172008
|
||
};
|
||
|
||
double GK_weight[5] = {
|
||
0.29552422471475286217,
|
||
0.26926671930999634918,
|
||
0.21908636251598204295,
|
||
0.14945134915058059038,
|
||
0.066671344308688137179
|
||
};
|
||
|
||
// 设置积分点数
|
||
int GK_points = 5;
|
||
|
||
// 理想气体热容、焓和熵的方程常数
|
||
double ThermConstants[NUMBEROFCOMPONENTS][11] = {
|
||
{-29776.4, 7.95454, 43.9417, 1037.09, 1.56373, 813.205, -24.9027, 1019.98,-10.1601, 1070.14,-20.0615},
|
||
{-3495.34, 6.95587, 0.272892, 662.738,-0.291318,-680.562, 1.78980, 1740.06, 0.0, 100.0, 4.49823},
|
||
{ 20.7307, 6.96237, 2.68645, 500.371,-2.56429,-530.443, 3.91921, 500.198, 2.13290, 2197.22, 5.81381},
|
||
{-37524.4, 7.98139, 24.3668, 752.320, 3.53990, 272.846, 8.44724, 1020.13,-13.2732, 869.510,-22.4010},
|
||
{-56072.1, 8.14319, 37.0629, 735.402, 9.38159, 247.190, 13.4556, 1454.78,-11.7342, 984.518,-24.0426},
|
||
{-13773.1, 7.97183, 6.27078, 2572.63, 2.05010, 1156.72, 0.0, 100.0, 0.0, 100.0, -3.24989},
|
||
{-10085.4, 7.94680,-0.08380, 433.801, 2.85539, 843.792, 6.31595, 1481.43,-2.88457, 1102.23,-0.51551},
|
||
{-5565.60, 6.66789, 2.33458, 2584.98, .749019, 559.656, 0.0, 100.0, 0.0, 100.0, -7.94821},
|
||
{-2753.49, 6.95854, 2.02441, 1541.22, .096774, 3674.81, 0.0, 100.0, 0.0, 100.0, 6.23387},
|
||
{-3497.45, 6.96302, 2.40013, 2522.05, 2.21752, 1154.15, 0.0, 100.0, 0.0, 100.0, 9.19749},
|
||
{-72387.0, 17.8143, 58.2062, 1787.39, 40.7621, 808.645, 0.0, 100.0, 0.0, 100.0, -44.1341},
|
||
{-72674.8, 18.6383, 57.4178, 1792.73, 38.6599, 814.151, 0.0, 100.0, 0.0, 100.0, -46.1938},
|
||
{-91505.5, 21.3861, 74.3410, 1701.58, 47.0587, 775.899, 0.0, 100.0, 0.0, 100.0, -60.2474},
|
||
{-83845.2, 22.5012, 69.5789, 1719.58, 46.2164, 802.174, 0.0, 100.0, 0.0, 100.0, -62.2197},
|
||
{-94982.5, 26.6225, 80.3819, 1718.49, 55.6598, 802.069, 0.0, 100.0, 0.0, 100.0, -77.5366},
|
||
{-103353., 30.4029, 90.6941, 1669.32, 63.2028, 786.001, 0.0, 100.0, 0.0, 100.0, -92.0164},
|
||
{-109674., 34.0847, 100.253, 1611.55, 69.7675, 768.847, 0.0, 100.0, 0.0, 100.0, -106.149},
|
||
{-122599., 38.5014, 111.446, 1646.48, 80.5015, 781.588, 0.0, 100.0, 0.0, 100.0, -122.444},
|
||
{-133564., 42.7143, 122.173, 1654.85, 90.2255, 785.564, 0.0, 100.0, 0.0, 100.0, -138.006},
|
||
{ 0.0, 4.9680, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0 },
|
||
{ 0.0, 4.9680, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0 }
|
||
};
|
||
|
||
// 全局变量
|
||
static double dSi = 0.0;
|
||
static double dTold = 0.0;
|
||
static double dMrxold = 0.0;
|
||
|
||
/**************************************************************************
|
||
* 函数: ThermInit
|
||
* 参数: void
|
||
* 返回: void
|
||
* 目的: 默认构造函数
|
||
* 修订:
|
||
**************************************************************************/
|
||
void ThermInit(void)
|
||
{
|
||
// 初始化3个历史敏感变量
|
||
dSi = 0.0;
|
||
dTold = 0.0;
|
||
dMrxold = 0.0;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: ThermDestroy
|
||
* 参数: void
|
||
* 返回: void
|
||
* 目的: 默认析构函数
|
||
* 修订:
|
||
**************************************************************************/
|
||
void ThermDestroy(void)
|
||
{
|
||
// 无动态资源需要释放
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: coth
|
||
* 参数: double x
|
||
* 返回: double
|
||
* 目的: 计算双曲余切函数; 用于Ho计算
|
||
* 修订:
|
||
* 备注: 不是Therm对象类成员，只是本文件的工具函数。
|
||
* C语言本身不支持双曲余切函数
|
||
**************************************************************************/
|
||
double coth(double x)
|
||
{
|
||
return cosh(x)/sinh(x);
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: Run
|
||
* 参数: AGA10STRUCT *ptAGA10, Detail *ptD
|
||
* 返回: void
|
||
* 目的: 总体执行控制; SOS和k的顶层数学计算
|
||
* 修订:
|
||
**************************************************************************/
|
||
void Run(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
{
|
||
// 局部变量
|
||
double c, x, y, z;
|
||
|
||
// 首先运行AGA 8(1994)详细方法中的基本函数集
|
||
Detail_Run(ptD, ptAGA10);
|
||
|
||
// 找到Z对D的一阶偏导数
|
||
Detail_dZdD(ptD, ptAGA10->dDf);
|
||
|
||
// 找到真实气体的cv, cp, 比焓和比熵
|
||
CprCvrHS(ptAGA10, ptD);
|
||
|
||
// 真实气体比热比
|
||
ptAGA10->dk = ptAGA10->dCp / ptAGA10->dCv;
|
||
|
||
// 分三步计算c，为了清晰和便于调试
|
||
x = ptAGA10->dk * RGAS * 1000.0 * ptAGA10->dTf;
|
||
y = ptAGA10->dMrx;
|
||
z = ptAGA10->dZf + ptAGA10->dDf * ptD->ddZdD;
|
||
|
||
// 计算c，即SOS^2
|
||
c = (x / y) * z;
|
||
|
||
// 声速
|
||
ptAGA10->dSOS = sqrt(c);
|
||
|
||
// 计算真实气体等熵指数
|
||
// 使用功能上等同于方程3.2的表达式
|
||
ptAGA10->dKappa = (c * ptAGA10->dRhof) / ptAGA10->dPf;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: CpiMolar
|
||
* 参数: AGA10STRUCT *ptAGA10
|
||
* 返回: double
|
||
* 目的: 计算理想气体定压摩尔热容(J/mol-K)，应用Aly, Lee, McFall的方程
|
||
* 备注: 为了连续性，保留了原始常数和方程。
|
||
* 主计算完成后应用从热化学卡路里(th)到焦耳的转换。
|
||
* 修订:
|
||
**************************************************************************/
|
||
double CpiMolar(AGA10STRUCT *ptAGA10)
|
||
{
|
||
double Cp = 0.0;
|
||
double Cpx;
|
||
double DT, FT, HT, JT;
|
||
double Dx, Fx, Hx, Jx;
|
||
double T;
|
||
int i;
|
||
|
||
// 为提高可读性，使用中间变量T
|
||
T = ptAGA10->dTf;
|
||
|
||
// 计算每个组分的比热容
|
||
for (i = 0; i < NUMBEROFCOMPONENTS; i++)
|
||
{
|
||
// 跳过浓度为零的组分
|
||
if (ptAGA10->adMixture[i] <= 0.0) continue;
|
||
|
||
// 初始化组分比热为零
|
||
Cpx = 0.0;
|
||
|
||
// 计算组分中间项
|
||
DT = ThermConstants[i][coefD] / T;
|
||
FT = ThermConstants[i][coefF] / T;
|
||
HT = ThermConstants[i][coefH] / T;
|
||
JT = ThermConstants[i][coefJ] / T;
|
||
|
||
// 使用中间项避免冗余计算
|
||
Dx = DT/sinh(DT);
|
||
Fx = FT/cosh(FT);
|
||
Hx = HT/sinh(HT);
|
||
Jx = JT/cosh(JT);
|
||
|
||
Cpx += ThermConstants[i][coefB];
|
||
Cpx += ThermConstants[i][coefC] * Dx * Dx;
|
||
Cpx += ThermConstants[i][coefE] * Fx * Fx;
|
||
Cpx += ThermConstants[i][coefG] * Hx * Hx;
|
||
Cpx += ThermConstants[i][coefI] * Jx * Jx;
|
||
|
||
// 使用当前摩尔分数加权贡献
|
||
Cpx *= ptAGA10->adMixture[i];
|
||
|
||
// 将此贡献加到总和中
|
||
Cp += Cpx;
|
||
}
|
||
|
||
// 从cal(th)/mol-K转换为J/mol-K
|
||
Cp *= CalTH;
|
||
|
||
return Cp;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: Ho
|
||
* 参数: AGA10STRUCT *ptAGA10
|
||
* 返回: double
|
||
* 目的: 计算理想气体比焓(J/kg)
|
||
* 备注: 为了连续性，保留了原始常数和方程。
|
||
* 主计算完成后应用从热化学卡路里(th)到焦耳的转换。
|
||
* 修订:
|
||
**************************************************************************/
|
||
double Ho(AGA10STRUCT *ptAGA10)
|
||
{
|
||
double H = 0.0;
|
||
double Hx;
|
||
double DT, FT, HT, JT;
|
||
double cothDT, tanhFT, cothHT, tanhJT;
|
||
double T;
|
||
int i;
|
||
|
||
// 为提高可读性，使用中间变量T
|
||
T = ptAGA10->dTf;
|
||
|
||
for (i = 0; i < NUMBEROFCOMPONENTS; i++)
|
||
{
|
||
// 跳过浓度为零的组分
|
||
if (ptAGA10->adMixture[i] <= 0.0) continue;
|
||
|
||
Hx = 0.0;
|
||
|
||
// 计算组分中间项
|
||
DT = ThermConstants[i][coefD] / T;
|
||
FT = ThermConstants[i][coefF] / T;
|
||
HT = ThermConstants[i][coefH] / T;
|
||
JT = ThermConstants[i][coefJ] / T;
|
||
|
||
cothDT = coth(DT);
|
||
tanhFT = tanh(FT);
|
||
cothHT = coth(HT);
|
||
tanhJT = tanh(JT);
|
||
|
||
Hx += ThermConstants[i][coefA];
|
||
Hx += ThermConstants[i][coefB] * T;
|
||
Hx += ThermConstants[i][coefC] * ThermConstants[i][coefD] * cothDT;
|
||
Hx -= ThermConstants[i][coefE] * ThermConstants[i][coefF] * tanhFT;
|
||
Hx += ThermConstants[i][coefG] * ThermConstants[i][coefH] * cothHT;
|
||
Hx -= ThermConstants[i][coefI] * ThermConstants[i][coefJ] * tanhJT;
|
||
|
||
// 使用当前摩尔分数加权贡献
|
||
Hx *= ptAGA10->adMixture[i];
|
||
|
||
// 将此贡献加到总和中
|
||
H += Hx;
|
||
}
|
||
|
||
// 从cal(th)/g-mol转换为kJ/kg-mol
|
||
H *= CalTH;
|
||
|
||
// 从kJ/kg-mol转换为J/kg
|
||
H /= ptAGA10->dMrx;
|
||
|
||
// 返回J/kg
|
||
return H * 1.e3;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: So
|
||
* 参数: AGA10STRUCT *ptAGA10
|
||
* 返回: double
|
||
* 目的: 计算理想气体比熵(J/kg-K)
|
||
* 备注: 为了连续性，保留了原始常数和方程。
|
||
* 主计算完成后应用从热化学卡路里(th)到焦耳的转换。
|
||
* 修订:
|
||
**************************************************************************/
|
||
double So(AGA10STRUCT *ptAGA10)
|
||
{
|
||
double S = 0.0;
|
||
double Sx;
|
||
double DT, FT, HT, JT;
|
||
double cothDT, tanhFT, cothHT, tanhJT;
|
||
double sinhDT, coshFT, sinhHT, coshJT;
|
||
double T;
|
||
int i;
|
||
|
||
// 为提高可读性，使用中间变量T
|
||
T = ptAGA10->dTf;
|
||
|
||
for (i = 0; i < NUMBEROFCOMPONENTS; i++)
|
||
{
|
||
// 跳过浓度为零的组分
|
||
if (ptAGA10->adMixture[i] <= 0.0) continue;
|
||
|
||
Sx = 0.0;
|
||
|
||
// 计算组分中间项
|
||
DT = ThermConstants[i][coefD] / T;
|
||
FT = ThermConstants[i][coefF] / T;
|
||
HT = ThermConstants[i][coefH] / T;
|
||
JT = ThermConstants[i][coefJ] / T;
|
||
|
||
cothDT = coth(DT);
|
||
tanhFT = tanh(FT);
|
||
cothHT = coth(HT);
|
||
tanhJT = tanh(JT);
|
||
|
||
sinhDT = sinh(DT);
|
||
coshFT = cosh(FT);
|
||
sinhHT = sinh(HT);
|
||
coshJT = cosh(JT);
|
||
|
||
Sx += ThermConstants[i][coefK];
|
||
Sx += ThermConstants[i][coefB] * log(T);
|
||
Sx += ThermConstants[i][coefC] * (DT * cothDT - log(sinhDT));
|
||
Sx -= ThermConstants[i][coefE] * (FT * tanhFT - log(coshFT));
|
||
Sx += ThermConstants[i][coefG] * (HT * cothHT - log(sinhHT));
|
||
Sx -= ThermConstants[i][coefI] * (JT * tanhJT - log(coshJT));
|
||
|
||
// 使用当前摩尔分数加权贡献
|
||
Sx *= ptAGA10->adMixture[i];
|
||
|
||
// 将此贡献加到总和中
|
||
S += Sx;
|
||
}
|
||
|
||
// 从cal(th)/mol-K转换为kJ/kg mol-K
|
||
S *= CalTH;
|
||
|
||
// 从kJ/kg mol-K转换为kJ/kg-K
|
||
S /= ptAGA10->dMrx;
|
||
|
||
// 返回J/kg-K
|
||
return S * 1.e3;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: CprCvrHS
|
||
* 参数: AGA10STRUCT *ptAGA10, Detail *ptD
|
||
* 返回: void
|
||
* 目的: 高效地分组计算Cp, Cv, H和S
|
||
* 修订:
|
||
**************************************************************************/
|
||
void CprCvrHS(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
{
|
||
double Cvinc, Cvr, Cpr;
|
||
double Hinc;
|
||
double Sinc;
|
||
double Smixing;
|
||
double Cp, Si;
|
||
double a, b, x;
|
||
int i;
|
||
|
||
// 初始化积分为零
|
||
Cvinc = 0.0;
|
||
Hinc = 0.0;
|
||
Sinc = 0.0;
|
||
|
||
// 初始化混合熵
|
||
Smixing = 0.0;
|
||
|
||
// 找到理想气体Cp
|
||
Cp = CpiMolar(ptAGA10);
|
||
|
||
// 找到理想气体焓
|
||
ptAGA10->dHo = Ho(ptAGA10);
|
||
|
||
// 找到理想气体熵
|
||
Si = So(ptAGA10);
|
||
|
||
// 计算理想气体定压比热容(J/kgK)
|
||
ptAGA10->dCpi = (Cp * 1000.0) / ptAGA10->dMrx;
|
||
|
||
// 从D=0到D=D积分偏导数，应用Gauss-Kronrod求积法
|
||
for (i = 0; i < GK_points; i++)
|
||
{
|
||
// 在+横坐标处设置计算点
|
||
x = ptAGA10->dDf * (1.0 + GK_root[i]) / 2.0;
|
||
|
||
// 在D处获取Z
|
||
Detail_zdetail(ptD, x);
|
||
Detail_dZdT(ptD, x);
|
||
Detail_d2ZdT2(ptD, x);
|
||
|
||
// 在+横坐标处收集贡献，应用权重因子
|
||
Hinc += GK_weight[i] * ptD->ddZdT / x;
|
||
Cvinc += GK_weight[i] * (2.0 * ptD->ddZdT + ptAGA10->dTf * ptD->dd2ZdT2) / x;
|
||
Sinc += GK_weight[i] * (ptD->dZ + ptAGA10->dTf * ptD->ddZdT - 1.0) / x;
|
||
|
||
// 在-横坐标处设置计算点
|
||
x = ptAGA10->dDf * (1.0 - GK_root[i]) / 2.0;
|
||
|
||
// 在D处获取Z
|
||
Detail_zdetail(ptD, x);
|
||
// 计算Z对T的一阶和二阶偏导数
|
||
Detail_dZdT(ptD, x);
|
||
Detail_d2ZdT2(ptD, x);
|
||
|
||
// 在-横坐标处收集贡献，应用权重因子
|
||
Hinc += GK_weight[i] * ptD->ddZdT / x;
|
||
Cvinc += GK_weight[i] * (2.0 * ptD->ddZdT + ptAGA10->dTf * ptD->dd2ZdT2) / x;
|
||
Sinc += GK_weight[i] * (ptD->dZ + ptAGA10->dTf * ptD->ddZdT - 1.0) / x;
|
||
}
|
||
|
||
// 将Z和偏导数恢复到全摩尔密度
|
||
Detail_zdetail(ptD, ptAGA10->dDf);
|
||
Detail_dZdT(ptD, ptAGA10->dDf);
|
||
Detail_d2ZdT2(ptD, ptAGA10->dDf);
|
||
|
||
// 完成Cv摩尔计算
|
||
Cvr = Cp - RGAS * (1.0 + ptAGA10->dTf * Cvinc * 0.5 * ptAGA10->dDf);
|
||
|
||
// Cp的中间值，包含2个偏导数
|
||
a = (ptAGA10->dZf + ptAGA10->dTf * ptD->ddZdT);
|
||
b = (ptAGA10->dZf + ptAGA10->dDf * ptD->ddZdD);
|
||
|
||
// 完成摩尔Cp计算的方程，取消适当的项
|
||
Cpr = Cvr + RGAS * ((a * a)/b);
|
||
|
||
// 从摩尔基准转换为质量基准
|
||
Cpr /= ptAGA10->dMrx;
|
||
Cvr /= ptAGA10->dMrx;
|
||
|
||
// 写入数据结构
|
||
ptAGA10->dCv = Cvr * 1000.0; // 从joules/kgK转换为kilojoules/kgK
|
||
ptAGA10->dCp = Cpr * 1000.0;
|
||
|
||
// 计算比焓
|
||
ptAGA10->dH = ptAGA10->dHo + 1000.0 * RGAS * ptAGA10->dTf *
|
||
(ptAGA10->dZf - 1.0 - ptAGA10->dTf * Hinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx;
|
||
|
||
// 计算混合熵
|
||
for (i = 0; i < NUMBEROFCOMPONENTS; i++)
|
||
{
|
||
if (ptAGA10->adMixture[i] > 0.0) // 避免log(0)
|
||
Smixing += ptAGA10->adMixture[i] * log(ptAGA10->adMixture[i]);
|
||
}
|
||
Smixing *= RGAS;
|
||
|
||
// 计算比熵
|
||
ptAGA10->dS = Si - Smixing - 1000.0 * RGAS *
|
||
(log(ptAGA10->dPf/101325.0) - log(ptAGA10->dZf) + Sinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: HS_Mode
|
||
* 参数: AGA10STRUCT *ptAGA10, Detail *ptD, double H_target, double S_target, bool bGuess
|
||
* 返回: void
|
||
* 目的: 从已知焓和熵计算压力和温度，有或没有先验估计。
|
||
* 此函数在计算C*时有作用。
|
||
* 解决方案基于双重嵌套的试错算法和牛顿法。
|
||
*
|
||
* 为说明目的，此示例支持两种方法。
|
||
* 如果您在不知道P和T的情况下开始，将输入参数bGuess设置为false，
|
||
* 从而指定保守的搜索方法。
|
||
* 但是，如果您有猜测P和T的基础(例如临界流喷嘴的集气室条件)，
|
||
* 通过AGA10STRUCT设置P和T，并设置bGuess = true。
|
||
* 初始猜测允许搜索函数更积极且通常更快。
|
||
* 修订:
|
||
**************************************************************************/
|
||
void HS_Mode(AGA10STRUCT *ptAGA10, Detail *ptD, double H_target, double S_target, bool bGuess)
|
||
{
|
||
double s0, s1, s2, t0, t1, t2, tmin, tmax;
|
||
double h0, h1, h2, p0, p1, p2, px, pmin, pmax;
|
||
double delta1, delta2;
|
||
double tolerance = 0.001; // 收敛容差(用于H和S搜索)
|
||
int i, j;
|
||
|
||
// s0和h0是我们的真实气体参考点
|
||
s0 = S_target;
|
||
h0 = H_target;
|
||
|
||
// 调用函数指定搜索参数是通过ptAGA10提供还是未知
|
||
if (bGuess)
|
||
{
|
||
t1 = ptAGA10->dTf;
|
||
px = ptAGA10->dPf;
|
||
pmax = px * 2.0;
|
||
pmin = px * 0.1;
|
||
tmax = t1 * 1.5;
|
||
tmin = t1 * 0.67;
|
||
}
|
||
else // 使用任意通用限制
|
||
{
|
||
t1 = 273.15;
|
||
px = 1013250.0; // 10个大气压
|
||
pmax = P_MAX;
|
||
pmin = 10000.0; // 10 kPa
|
||
tmax = T_MAX;
|
||
tmin = T_MIN;
|
||
}
|
||
|
||
// 设置温度差
|
||
t2 = t1 + 10.0;
|
||
|
||
// 开始双重试错，搜索T和P
|
||
// 用初始猜测运行计算
|
||
Detail_Run(ptD, ptAGA10);
|
||
|
||
// h1是给定h和h@Tf, Pf之间的差
|
||
h1 = H(ptAGA10, ptD) - h0;
|
||
|
||
// 外循环: 搜索满足恒定焓的t2
|
||
for (i = 0; i < MAX_NUM_OF_ITERATIONS; i++)
|
||
{
|
||
ptAGA10->dTf = t2;
|
||
p1 = px; // 重置一个括号
|
||
p2 = px * 0.1; // 将另一个括号设置为上括号的0.1倍
|
||
ptAGA10->dPf = p1;
|
||
Detail_Run(ptD, ptAGA10);
|
||
s1 = S(ptAGA10, ptD) - s0;
|
||
|
||
// 内循环: 搜索满足恒定熵的p2
|
||
for (j = 0; j < MAX_NUM_OF_ITERATIONS; j++)
|
||
{
|
||
ptAGA10->dPf = p2;
|
||
Detail_Run(ptD, ptAGA10);
|
||
s2 = S(ptAGA10, ptD) - s0;
|
||
|
||
// 计算我们的比例变化
|
||
delta2 = fabs(s1 - s2) / s0;
|
||
|
||
// 足够接近?
|
||
if (delta2 < tolerance) break;
|
||
|
||
// 将我们的估计修正为p2
|
||
p0 = p2;
|
||
p2 = (p1 * s2 - p2 * s1) / (s2 - s1);
|
||
|
||
// 检查负压并钳制到pmin以确保安全
|
||
if (p2 <= pmin)
|
||
{
|
||
p2 = pmin;
|
||
}
|
||
|
||
// 检查是否创建了不切实际的压力
|
||
if (p2 >= pmax) p2 = pmax;
|
||
|
||
// 交换值
|
||
p1 = p0;
|
||
s1 = s2;
|
||
}
|
||
|
||
// 检查是否未能收敛
|
||
if (j >= MAX_NUM_OF_ITERATIONS)
|
||
ptAGA10->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED;
|
||
|
||
// 在猜测的P和当前迭代T处计算焓
|
||
h2 = H(ptAGA10, ptD) - h0;
|
||
|
||
// 计算我们的比例变化
|
||
delta1 = fabs(h1 - h2) / h0;
|
||
|
||
// 足够接近?
|
||
if (delta1 < tolerance && i > 0) break;
|
||
|
||
// 将我们的估计修正为t2
|
||
t0 = t2;
|
||
t2 = (t1 * h2 - t2 * h1) / (h2 - h1);
|
||
|
||
// 检查是否创建了不切实际的温度
|
||
if (t2 >= tmax) t2 = tmax;
|
||
|
||
// 如有必要，修正t2
|
||
if (t2 <= tmin)
|
||
{
|
||
t2 = t0 + 10.0;
|
||
ptAGA10->dTf = t2;
|
||
Detail_Run(ptD, ptAGA10);
|
||
h2 = H(ptAGA10, ptD) - h0;
|
||
}
|
||
|
||
t1 = t0;
|
||
h1 = h2;
|
||
}
|
||
|
||
// 检查是否未能收敛
|
||
if (i >= MAX_NUM_OF_ITERATIONS)
|
||
ptAGA10->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: H
|
||
* 参数: AGA10STRUCT *ptAGA10, Detail *ptD
|
||
* 返回: double
|
||
* 目的: 真实气体比焓
|
||
* 修订:
|
||
**************************************************************************/
|
||
double H(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
{
|
||
double Hinc;
|
||
double x;
|
||
int i;
|
||
|
||
// 初始化积分
|
||
Hinc = 0.0;
|
||
|
||
// 找到理想气体焓
|
||
ptAGA10->dHo = Ho(ptAGA10);
|
||
|
||
// 从D=0到D=D积分偏导数，应用Gauss-Kronrod求积法
|
||
for (i = 0; i < GK_points; i++)
|
||
{
|
||
// 计算Z对T的一阶和二阶偏导数
|
||
x = ptAGA10->dDf * (1.0 + GK_root[i]) / 2.0;
|
||
Detail_zdetail(ptD, x);
|
||
Detail_dZdT(ptD, x);
|
||
Detail_d2ZdT2(ptD, x);
|
||
|
||
Hinc += GK_weight[i] * ptD->ddZdT / x;
|
||
if (i == 10) break;
|
||
|
||
x = ptAGA10->dDf * (1.0 - GK_root[i]) / 2.0;
|
||
Detail_zdetail(ptD, x);
|
||
Detail_dZdT(ptD, x);
|
||
Detail_d2ZdT2(ptD, x);
|
||
|
||
Hinc += GK_weight[i] * ptD->ddZdT / x;
|
||
}
|
||
|
||
Detail_zdetail(ptD, ptAGA10->dDf);
|
||
Detail_dZdT(ptD, ptAGA10->dDf);
|
||
Detail_d2ZdT2(ptD, ptAGA10->dDf);
|
||
|
||
// 计算比焓
|
||
ptAGA10->dH = ptAGA10->dHo + 1000.0 * RGAS * ptAGA10->dTf *
|
||
(ptAGA10->dZf - 1.0 - ptAGA10->dTf * Hinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx;
|
||
|
||
return ptAGA10->dH;
|
||
}
|
||
|
||
/**************************************************************************
|
||
* 函数: S
|
||
* 参数: AGA10STRUCT *ptAGA10, Detail *ptD
|
||
* 返回: double
|
||
* 目的: 真实气体比熵
|
||
* 修订:
|
||
**************************************************************************/
|
||
double S(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
{
|
||
double Sinc;
|
||
double Smixing;
|
||
double x;
|
||
int i;
|
||
|
||
// 初始化积分
|
||
Sinc = 0.0;
|
||
|
||
// 初始化混合熵
|
||
Smixing = 0.0;
|
||
|
||
// 从D=0到D=D积分偏导数，应用Gauss-Kronrod求积法
|
||
for (i = 0; i < GK_points; i++)
|
||
{
|
||
// 计算Z对T的一阶和二阶偏导数
|
||
x = ptAGA10->dDf * (1.0 + GK_root[i]) / 2.0;
|
||
Detail_zdetail(ptD, x);
|
||
Detail_dZdT(ptD, x);
|
||
Detail_d2ZdT2(ptD, x);
|
||
|
||
Sinc += GK_weight[i] * (ptD->dZ + ptAGA10->dTf * ptD->ddZdT - 1.0) / x;
|
||
if (i == 10) break;
|
||
|
||
x = ptAGA10->dDf * (1.0 - GK_root[i]) / 2.0;
|
||
Detail_zdetail(ptD, x);
|
||
Detail_dZdT(ptD, x);
|
||
Detail_d2ZdT2(ptD, x);
|
||
|
||
Sinc += GK_weight[i] * (ptD->dZ + ptAGA10->dTf * ptD->ddZdT - 1.0) / x;
|
||
}
|
||
|
||
// 重置Z和偏导数dZdT和d2ZdT2
|
||
Detail_zdetail(ptD, ptAGA10->dDf);
|
||
Detail_dZdT(ptD, ptAGA10->dDf);
|
||
Detail_d2ZdT2(ptD, ptAGA10->dDf);
|
||
|
||
// 找到理想气体熵，但仅在温度或组成发生变化时
|
||
if (ptAGA10->dTf != dTold || ptAGA10->dMrx != dMrxold)
|
||
{
|
||
dSi = So(ptAGA10);
|
||
dTold = ptAGA10->dTf;
|
||
dMrxold = ptAGA10->dMrx;
|
||
}
|
||
|
||
// 计算混合熵
|
||
for (i = 0; i < NUMBEROFCOMPONENTS; i++)
|
||
{
|
||
if (ptAGA10->adMixture[i] > 0.0) // 避免log(0)
|
||
Smixing += ptAGA10->adMixture[i] * log(ptAGA10->adMixture[i]);
|
||
}
|
||
Smixing *= RGAS;
|
||
|
||
// 计算比熵
|
||
ptAGA10->dS = dSi - Smixing - 1000.0 * RGAS *
|
||
(log(ptAGA10->dPf/101325.0) - log(ptAGA10->dZf) + Sinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx;
|
||
|
||
return ptAGA10->dS;
|
||
}
|