Определение параметров начала процесса расширения в детандере методом П.Л. Капицы по h-s диаграмме

Постановка задачи

В процессе расчета параметров циклов криогенных систем с детандером часто возникает необходимость решения задачи по определению параметров точки начала процесса расширения по известным параметрам в точке конце расширения - на выходе из детандера. Например, такая задача решается при расчете газового цикла Сименса-Клода в онлайн-калькуляторе на сайте.

В статье рассмотрим два подхода к решению задачи:

• графический метод - альтренатива формуле П.Л. Капицы для определения парамтеров начала процесса расширения в детандере с использованием h-s диаграммы;
• численный метод решения системы уравнений, связывающих соответствующие параметры процесса расширения и характеристику детандера, с использованием библиотеки для расчета термодинамических свойств рабочих агентов.

Для решения задачи потребуются следующие исходные данные:

• h-s диаграмма рабочей среды или ее модель в библитеке для расчета термодинамических свойств;
• температура и давление на выходе из детандера (в конце процесса расширения в детнадере);
• давление на входе в детандер (в начале процесса расширения в детандере);
• изоэнтропный коэффициент полезного действия детендера.

Для большей ясности рассмотрим оба подхода на примере с конкретными численными значениями, представленными в таблице 1.

 Графический метод П.Л. Капицы

Графический метод - аналогия формуле Петра Леонидовича Капицы применяется в случаях, когда расчет ведется для одного или нескольких вариантов исходных параметров криогенной системы, а также в ситуациях, когда нет возможности провести соответствующие вычисления с использованием численного метода. Поскольку метод относится к графическим, то точность расчета прямо определяется точностью и качеством построения соответствующих отрезков и размещения расчетных точек.

Шаг 1. Подготовка h-s диаграммы - выбор рабочей среды

Наиболее простой способ получить h-s диаграмму для принятого рабочего агента - воспользоваться программой RefProp, которую можно бесплатно загрузить в сети Интернет или по ссылке на странице в разделе Работа с библиотекой термодинамических свойств RefProp в Mathcad 15.

После установки программы RefProp следует (см. рисунок 1):

1. На панели управления нажать на кнопку Substance;
2. В выпадающем меню выбрать тип рабочей среды - в данном случае Pure Fluid (чистое вещество);
3. В появившемся окне изучить список доступных сред;
4. Выбрать требуемое вещество - в данном случае Nitrogen (азот);
5. Подтвердить выбор рабочей среды нажатие на кнопу Ok.

Шаг 2. Подготовка h-s диаграммы - построение диаграммы

После выбора рабочей среды следует построить h-s диаграмму, для этого следует (см. рисунок 2):

1. На панели управления нажать на кнопку Plot;
2. В выпадающем меню выбрать тип диаграммы для построения - в данном случае h-s Diagram;
3. В появившемся окне настроить параметры диаграммы - границы диаграммы по этальпии / энтропии и изолинии температур / давления;
4. Подтвердить введенные данные и провести построение диаграммы (см. рисунок 3) нажатием кнопки Ok.

Если получившаяся диаграмма неудобна для построения или содержит нечитаемые надписи, следует повторить шаги 3 и 4.

Шаг 3. Определение параметров процесса расширения

На получившейся диаграмме выполняются следующие построения (см. рисунок 4):

• Определяется положение точки 4 в месте пересечения изотермы (110 K) и изобары (1,05 бар.А).

Если расчет ведется для парожидкостного детандера, то дополнительно используется значение этальпии или энтропии (важно: в той же системе отсчета) в точке 4.

• Через точку 4 проводится горизонталь;
• В произвольной точке b опускается перпендикуляр к ранее проведенной горизонтали до пересечения с изобрарой, проходящей через точку 4 (1,05 бар.А), в точке c;
• Измеряется длина опущенного перпендикуляра [bc] (190 px) в любой системе отсчета (px, энтальпия, см и пр.);
• Рассчитывается длина отрезка [ba] из соотношения:

• Проводится отрезок [ba], перпендикулярный ранее проведенной горизонтали;
• Проводится отрезок [a4];
• В месте пересечения проведенного отрезка [a4] с изобарой (12 бар.А) отмечается точка 3;
• Опускается вертикальный отрезок до пересечения с изобарой в точке 4 (1,05 бар.А), отмечается точка 3s.

Шаг 4. Результаты построений

Температура в точке 3 совпадает с проходящей через эту точку изотермой ~ 182 К.

Температура в точке 3s совпадает с проходящей через эту точку изотермой ~ 90 К.

Численный метод

Численный метод основан на динамическом использование теплофизических параметров рабочей среды, рассчитываемых с применением специализированной библиотеки, в промежуточных точках относительно известных граничных условий. В качестве рабочей применяется динамическая библиотека CoolProp для языка программирования JavaScript. Узнать подробнее о библиотеке можно на официальном сайте.

Суть метода основана на итеративном подборе температуры начала процесса расширения в детандере в некотором диапазоне температур (startT и endT). Для каждой итерации происходит прямой расчет параметров процесса расширения.

В качестве критерия сходимости решения принимается совпадение действительного значения энтальпии в конце процесса расширения и энтальпии, полученной для i-той итерации. Совпадение оценивается методом сравнения модуля разницы с некоторой наперед заданной погрешностью.

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

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

Функция \(h_4(T_i = x)\), определяющая параметры точки начала процесса расширения, в итерационном процессе имеет следующий вид:

/**
 * Возвращает температуру начала процесса расширения по известным параметрам
 *
 * @public
 * @param {number} startPressure - Давление в начале процесса расширения [12 бар]
 * @param {number} endRessure - Давление в конце процесса расширения [1,05 бар]
 * @param {number} endEnthalpy - Энтальпия в конце процесса расширения [112.50 кДж/кг]
 * @param expanderEff - Изоэнтропный КПД детандера [0,76]
 * @param fluid - Название рабочей среды, ["nitrogen"]
 * @param startT - Начальная температура численного поиска решения [Toc = 300 K]
 * @param endT - Конечная температура численного поиска решения [T4 = 110 K]
 * @returns {number} - Температура на входе в детандер, K
 *
 */
calculateExpansionStartTemperature(startPressure, endRessure, endEnthalpy, expanderEff, fluid, startT, endT) {
  return EquationSolver.solve((x) => {
    let h1 = Module.PropsSI('H', 'T', x, 'P', startPressure * 100000 fluid); // бар =  100000 Па
    let s1 = Module.PropsSI('S', 'T', x, 'P', endRessure * 100000, fluid); // бар =  100000 Па
    let h1s = Module.PropsSI('H', 'S', s1, 'P', endRessure * 100000, fluid); // бар =  100000 Па
    let h2 = endEnthalpy * 1000; // кДж/кг =  1000 Дж/кг

    return h1 - h2 - expanderEff * (h1 - h1s);
  }, startT, endT);
}

Для ранее определенных исходных данных таблица с итеративным поиском решения выглядит следующим образом:

// Экспорт класса EquationSolver по умолчанию (default export)
export default class EquationSolver {

  // Статический метод solve для приближённого поиска корня уравнения f(x) = 0
  // Параметры:
  //   f: Function        - функция, для которой ищется корень
  //   leftSide: number   - левая граница начального интервала
  //   rightSide: number  - правая граница начального интервала
  //   accuracy: number   - количество значащих цифр после запятой
  static solve(f: Function, leftSide: number, rightSide: number, accuracy: number = 2): number {

    // Вычисляется допустимая погрешность по длине интервала
    let accuracyDegree = 1 / 10 ** (accuracy + 2);

    // Инициализируется левая и правая границы текущего интервала
    let left = leftSide;
    let right = rightSide;

    // Вычисляется середина интервала
    let middle = (right + left) / 2;

    // Проверяем, имеют ли значения функции на концах начального интервала одинаковый знак.
    // Если произведение положительно (либо оба >0, либо оба <0), то метод половинного деления
    // не гарантирует нахождение корня (либо корней несколько, либо нет вообще).
    if (f(leftSide) * f(rightSide) > 0) {

      // Счётчик шагов для предотвращения бесконечного цикла
      let stepCount = 0;

      // Цикл выполняется, пока ширина интервала больше допустимой погрешности
      for (let i = 1; right - left > accuracyDegree; i++) {

        // Если сделано больше 32 итераций, считается, что корень не найден
        if (stepCount > 32) {
          middle = Infinity; // Устанавливаем middle в Infinity как признак неудачи
          break;             // Прерываем цикл
        }

        // Стандартный шаг метода половинного деления:
        // Если функция меняет знак на отрезке [left, middle], то корень находится слева
        if (f(left) * f(middle) < 0) {
          right = middle;    // Сдвигает правую границу в середину
        } else {
          left = middle;     // Иначе корень справа, сдвигает левую границу
        }

        // Пересчитывается середина нового интервала
        middle = (right + left) / 2;

        // Увеличивается счетчик шагов
        stepCount++;
      }

    } else {
      // Случай, когда знаки на концах разные (произведение <= 0).
      for (let i = 1; right - left > accuracyDegree; i++) {

        // Аналогичный шаг: определяется, в какой половине корень
        if (f(left) * f(middle) < 0) {
          right = middle;
        } else {
          left = middle;
        }

        // Новая середина интервала
        middle = (right + left) / 2;
      }
    }

    // Возвращаем найденное приближение корня (середину последнего интервала).
    // Унарный плюс преобразует строку от toFixed обратно в число.
    return +middle.toFixed(accuracy + 2);
  };
 }