Число молей продуктов сгорания, приходящихся на 1 кг топлива, равно

K_CO2 = C/12.01;

K_H2O = H/2*1.008;

K_N2атм = (r_N2)_в/ _вq_т;

K_O2 = (r_O2)_в / _вq_т - .

Суммарное число продуктов сгорания

K_пс = 1/ _вq_т + ; т.к. (r_N2)_в + (r_O2)_в = 1.

Для нормального топлива

K_пс = (0.9280 + q_т) / 2688 q_т (2.3.3)

Объемный состав продуктов сгорания вычисляется по формулам

r_CO2 = C / 12.01*K_пс = 1.9*q_т / (0.9280 + q_т); (2.3.5)

r_H2O = H / 2*1.008*K_пс = 2*q_т / (0.9280 + q_т); (2.3.6)

r_O2 = 1 / K_пс ((r_O2)_в / _в*q_т – C / 12.01 – H / (4*1.008) +O / 32) =

= (0.1944 - 2.902*q_т) / (0.9280 + q_т); (2.3.7)

r_N2атм = (r_N2)_в / ( _в*q_т*K_пс) = 0.7336 / (0.9280 + q_т). (2.3.8)

Кажущийся молекулярный вес продуктов сгорания

_пс = _CO2 r_CO2 + _H2O r_H2O + _O2 r_O2 + _N2атмr_N2 =

= 26.88 (1+q_т) / (0.9280 + q_т). (2.3.9)

Истинная теплоемкость продуктов сгорания

с_pu = [r_CO2( c_p)_CO2 + r_H2O ( c_p)_H2O + r_O2( c_p)_O2 + r_N2( c_p)_N2атм],

или подставляя (2.4.4), (2.4.5), (2.4.6), (2.4.7) и (2.4.8)

с_pu = 1 / (1+q_т){b_O2² ( c_p)_O2 + b_N2( c_p)_N2атм + [b_CO2( c_p)_CO2 +

+ b_H2O( c_p)_H2O + b_O2^”( c_p)_O2 ]q_т}. (2.3.10)

Коэффициенты b приведены в табл. 2.3.1.

Для определения теплоемкостей газов использовались полиномы вида

с_p = _i tⁱ, (2.3.11)

полученные Ривкиным С.Л. [10].

Где t = T/1000 - температура (в К) с масштабным множителем (1/1000).

Коэффициенты полиномов a_i, для используемых газов приведены в табл. 2.3.2. После подстановки (2.3.11) в (2.3.10) получается

с_pu = 1/(1+q_т)[b^’_{O2 iO2} tⁱ + b_{N2 iN2} tⁱ + b_{CO2 iCO2} tⁱ +

+ b_{H2O iH2O} tⁱ + b^’’_{O2 iO2} tⁱ)q_т]

или, изменяя порядок суммирования

с_pu = 1/(1+q_т) ⁱ [a_iO2b^’_O2 + a_iN2 b_N2 + q_т(a_iCO2 b_CO2 + a_iH2O b_H2O + a_iO2 b^’’_O2)]

или, полагая

a_ib = a_iO2 + a_iN2b_N2,

a_iг = a_iCO2 + a_iH2Ob_H2O + a_iO2b^’’_O2, (2.3.12)

с_pu = 1/(1+q_т) ⁱ (a_ib + a_iг*q_т).

Значения коэффициентов a_ib, a_iг приведены в табл. 2.3.2. Размерность истинной теплоемкости с_pu в этом выражении [кДж/кг×К].

Энтальпия в кДж/кг и энтропия в кДж/кг×К определяются в результате интегрирования выражения (2.3.12) по формулам

H = c_pudT = 1000/(1+q_т) ⁱ (a_ib + a_iг*q_т)dt =

= A_i + 1000/(1+q_т) (a_ib + a_iг*q_т) tⁱ⁺¹/(i+1), (2.3.13)

S = c_pu/T dT = 1000/(1+q_т) ^i-1 (a_ib + a_iг*q_т)dt =

= A_s + 1000/(1+q_т) (a_ib + a_iг*q_т) tⁱ/i, (2.3.14)

где условно обозначено t^o = ln t.

Коэффициенты A_i и A_s приведены в табл. 2.3.2.

Таблица 2.3.1

Таблица 2.3.2

A_i = - 1.8741592 + 242.1708460 A_s = +7.94739 + 4.236147

В подпрограмме предусмотрена возможность расчета энтальпии сгоревшего топлива Н_топл. При этом полагается в уравнении (2.3.13) а_ib =0 и q_т = 1.

Универсальная газовая постоянная принята равной

R = 8314.2 Дж кмоль×К

Отсюда

R_пс = R/ _пс = кДж/кг×К (2.3.15)

Формулы (2.3.13), (2.3.14), (2.3.15) обеспечивают прямой расчет свойств, а в модели нужен и обратный, т.е. зависимости

T = f (H),

T = f (S).

Они вычисляются методом последовательных приближений. Для ускорения сходимости полиномы были аппроксимированы параболами вида

t = A_i + B_iH + C_iH², (2.3.16)

t = A_s + B_sS + C_sS². (2.3.17)

коэффициенты которых включены в программу. В процессе последовательных приближений при вычислении температуры по энтропии используются формулы:

t = A_s + B_sS_зад + С_sS²_зад,

S_n = f (t_n) по формуле (2.4.13)

t_n+1 = t_n + (S_зад - S_n) [B_s + C_s (S_зад + S_n)] (2.3.18)

S_n+1 = f(t_n+1) и т.д.,

где n - номер приближения.

Формула (2.3.18) обеспечивает смещение параболы (2.3.17) вдоль оси ординат за счет изменения свободного члена А_s = x. Cмещенная парабола должна проходить через две точки (t_n, S_n) и (t_n+1, S_зад),т.е.

t_n = x + B_sS_n + C_sS²_n ; t_n+1 = x + B_sS_зад + C_sS²_зад.

После исключения "х" получается формула (2.3.18). Температура по энтальпии рассчитывается аналогично. При заданной точности расчета 10^-4 число последовательных приближений обычно равняется 3:6.

В модулях математической модели используются также значения критической температуры Т_кр, которая должна вычисляться в зависимости от температуры или энтальпии торможения. Значение критической температуры в нулевом приближении определяется по приближенной параболе

t_кр A_кр + В_крН^*_зад + С_крН^*2_зад.

Соответствующие ей значения теплоемкости и энтальпии рассчитываются по формулам (2.3.12), (2.3.13). По ним определяется энтальпия торможения.

Н^* = Н + _кр*1000. (2.3.19)

Если она не совпала с заданной, то критическая температура уточняется по формуле, аналогичной (2.4.17)

t_кр2 = t_кр1 + (Н^*_зад-Н^*)[B_кр + С_кр(Н^*_зад + Н^*)]

и повторяется расчет Н^* по формуле (2.3.10) с учетом (2.3.12) и (2.3.13). Если вместо энтальпии задана температура торможения Т^*, то предварительно по формуле (2.3.13) рассчитывается энтальпия Н^*, а затем определяется Т_кр.

Данная программа позволяет за одно обращение к ней определить одну или сразу несколько термодинамических функций. Наименование рассчитываемой функции задается с помощью параметра “PR” в соответствии с нижеследующей табл. 2.3.3.

Таблица 2.3.3

Обращение к подпрограмме имеет вид

CALL FUNKZI (X, F, Q, P, S),

где X - аргумент, в зависимости от которого вычисляется функция;

Q - относительный расход топлива (q_т);

S - количество термодинамических функций, которое требуется рассчитать за одно обращение к модулю;

“P” - массив, заполняемый величинами “PR” в той последовательности, в которой должны рассчитываться требуемые функции;

“F” - массив результатов. Размерность массивов “F”,”P” должна быть не меньше “S”.

Последовательность заполнения массива “F” результатами соответствует последовательности “PR” в массиве “P”.

В подпрограмме предусмотрена возможность, при расчете каждой последующей функции, аргумент заменять результатом расчета предыдущей функции. Для этого к соответствующему значению “PR” в массиве “Р” необходимо прибавлять число “100”. Например, если в массиве S = 3, X = H^*, Q = q_т, а элементы массива “Р” равны: Р₁=1, Р₂=6, Р₃=105, то работа подпрограммы будет происходить в следующей последовательности:

1. Определится R = f(q_т)

2. Определится T^* = f(H^*,q_т)

3. На место аргумента Н^* занесется величина Т^*

4. Определится S^* = f(T^*,q_т).

В результате массив “F” заполнится в следующей последовательности: F₁=R, F₂=T^*, F₃=S^*. На месте аргумента “Х” после работы подпрограммы будет величина Т^*. После расчета функции Т_кр = f(T^*,q_т) на месте аргумента “Х” всегда будет находиться Н^*, соответствующая Т^* и q_т.

При описании алгоритмов модулей узлов, приводимом ниже, в местах, где необходимо определить термодинамические свойства газов, будут вместо расчетных формул приводиться функциональные зависимости типа Т^* = f(H^*, q_т), Н^* = f(Т^*, q_т), S^* = f(T^*, q_т) и т.д., предполагается, что в данном случае этот расчет выполняется при помощи подпрограммы “FUNKZI”.