Программные продукты и системы

При проектировании и программной реализации моделей сложных систем разработчику иногда приходится сталкиваться с задачами моделирования движения потоков различных жидкостей или газов в трубах, их распределения, взаимодействия и т.д. Примерами таких систем могут служить дренажные системы, системы вентиляции и кондиционирования помещений, газо- и водоснабжения, пожаротушения, подачи топлива в двигатель, подачи масла в различные агрегаты и др. При этом сами они могут быть лишь составной частью более крупных систем или взаимодействовать друг с другом. Для некоторых из систем известны относительно простые модели, другие же, наоборот, сложны настолько, что их адекватное моделирование в реальном времени возможно лишь с привлечением суперкомпьютеров [1–4].

Затруднения в моделировании подобных систем могут быть связаны не только с гро- моздкостью (числом составляющих их элементов и связей между ними), но и со сложностями компьютерной реализации известных для них математических моделей и их построением. 

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

В работе затронуты вопросы практического моделирования переноса жидкости и газа в системе, состоящей из большого количества емкостей, соединенных между собой трубами различных длин и диаметров. При этом допускается, что емкости могут иметь разный объем, содержать как жидкость, так и газ. При моделировании учитываются высота их расположения и присоединения труб, а также давление внутри емкостей. Цель исследования – получить модель, позволяющую вычислять в реаль- ном времени состояния этих емкостей (то есть содержащуюся в них массу жидкости и газа), а также величины потоков жидкости и газа в трубах. Для этого пожертвуем точностью модели, но предусмотрим возможность настройки желаемой точности вычислений. Кроме того, будут использованы технические приемы и эвристические соображения, позволяющие модели кратковременно оказываться в состояниях, невозможных с физической точки зрения. Результатом такого подхода будет возможность применения обычных компьютеров, а не суперкомпьютеров.

Предлагаемая модель (точнее, комплекс технических приемов) не претендует на то, чтобы удовлетворить все возможные ожидания относительно точности и быстроты вычислений. Тем не менее, она может использоваться для решения некоторых задач, где требуется не столько высокая точность реагирования модели на изменение ее состояния, сколько сам факт более-менее адекватного реагирования в реальном времени. Отметим, что предлагаемое описание модели позволяет понять, где и как можно произвести настройку модели таким образом, чтобы сместить баланс между точностью и скоростью в нужную сторону.

Модель была реализована авторами в виде компьютерной программы и протестирована для систем с очень разными свойствами, поэтому там, где это уместно, будут даны комментарии, основанные на результатах тестирования.

Система понятий и обозначений

Опишем понятия, которые будут использоваться в работе.

Емкость. Характеризуется своими геометрическими размерами, из которых будут важны объем V и высота h. Координаты расположения емкости в пространстве соответствуют координатам ее геометрического центра.

Труба. Характеризуется длиной l и площадью сечения S.

Узел. Используется для присоединения к нему концов нескольких труб. Характеризуется координатами расположения.

Жидкость. Характеризуется плотностью r и коэффициентом динамической вязкости h.

Газ. Характеризуется молярной массой M и коэффициентом динамической вязкости h.

Каждая труба соединяет два элемента – две емкости, два узла или же емкость и узел. Коор- динаты присоединения трубы к емкости необходимо задать, а координаты ее присоединения к узлу совпадают с координатами узла. Движение вещества по трубе характеризуется скоростью потока G. Для определенности считаем, что положительное значение G всегда соответствует некоторому фиксированному направлению в трубе, это направление называем положительным. Примем, что в случае, когда труба соединяет узел и емкость, положительным считается направление от узла к емкости.

Емкость содержит некоторое количество жидкости и газа. В каждом случае это количество характеризуется массой m. Во внутренних точках емкостей, труб и узлов возникает давление P, которое также будет учитываться в модели.

При описании модели к указанным обозначениям для величин будем присоединять нижние индексы, значение которых станет понятным из контекста. В случае описания итеративных процессов номер итерации будет указываться в круглых скобках в качестве верхнего индекса вычисляемой величины. Так, например, обозначение  следует понимать как величину потока жидкости в трубе номер i, вычисленную на k-й итерации.

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

Считаем, что единицы измерения для всех используемых в работе величин соответствуют СИ (ГОСТ 8.417–2002).

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

Цель – построить компьютерную модель, позволяющую наблюдать динамическое развитие системы емкостей, узлов и труб с газом и жидкостью в них, вычисляя новые состояния системы через некоторое время Dt на основе уже имеющихся. Поскольку собственно приращение времени Dt можно задавать в цикле, построение такой модели сводится, по сути, к описанию выполнения одного шага этого цикла.

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

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

Выбор математической модели

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

Уравнения первого вида описывают изменение потока в трубе:

,                                (1)

где DP – перепад давления на концах трубы; m – коэффициент трения, зависящий от диаметра трубы D, ее длины l, площади сечения S, коэффициента потерь на трение по длине l (зависящего от h, скорости потока, особенностей трубы) и плотности вещества r,

Поясним, откуда берется (1). Величина потока G вычисляется по формуле G = rSv (v – скорость потока), а импульс p – по формуле p = mv, то есть p = lG. В результате находим, что dp/dt = l×dG/dt. Но, согласно второму закону Ньютона, dp/dt = S(DP - Pп), где Pп – потери давления, связанные с трением, возникающим при движении вещества по трубе. Приравнивая правые части двух последних формул, получаем, что

Величину Pп можно найти по формуле Дарси–Вейсбаха:

Учитывая, что v = G/(rS), получаем абсолютное значение Pп, равное mG2 при указанном выше m. Наконец, учитывая направление по- тока, приходим к (1).

Уравнения второго вида описывают изменение массы вещества в емкостях:

                                                 (2)

Суммирование величин потоков в формуле (2) ведется по всем трубам, присоединенным к рассматриваемой емкости.

Уравнения третьего вида описывают условие, позволяющее вычислять давление в узлах:

                                                           (3)

Суммирование в формуле (3) ведется по всем трубам, присоединенным к рассматриваемому узлу.

По уравнениям вида (1)–(3) несложно построить вычислительную модель, заменив в них производные по времени  и  (то есть dG/dt и dm/dt) на отношения DG/Dt и Dm/Dt соответственно, а затем организовать вычисления при малых значениях Dt. Сами вычисления можно организовать на основе какого-либо из многих известных методов: явного или неявного метода Ньютона, метода Рунге–Кутты или других [5].

О решении задачи

Формально задача решена. Но на практике возникают некоторые трудности в зависимости от выбранного метода.

·     Даже при довольно небольшом значении Dt может оказаться, что модель приходит в состояние, невозможное с физической точки зрения. Например, в ряде случаев, согласно результатам вычислений, емкости должны принять жидкости больше, чем могут вместить, или же отдать больше, чем содержат. Как правило, подобные трудности преодолеваются путем дробления величины шага, однако это не всегда возможно, к тому же такой подход сам порождает новые трудности.

·     Существуют ограничения на дробление величины шага. С одной стороны, слишком большое значение Dt приводит к неадекватному результату, а с другой – достаточно малые значения Dt (обеспечивающие нужную точность вычислений) приводят к тому, что на вычисление состояния модели через интервал Dt требуется время, на несколько порядков превышающее Dt. Варьирование же значения Dt не дает существенного преимущества, так как, если это сделать в последний момент, модель может не удержаться в границах допустимых значений величин. Постоянные откаты на несколько шагов назад с последующим дроблением величины шага и пересчетом тоже не приводят к успеху, если нужно получать результаты в реальном времени.

·        Существует проблема, связанная со сходимостью используемых численных методов. Практика показала, что неудачно выбранный метод или величина шага могут привести к таким неприятным эффектам, как, например, «болтанка».

·     Пока компьютерная модель вычисляет состояния системы только с газом или только с жидкостью, упомянутые выше трудности проявляются не так остро и решаются относительно несложно. При одновременном же наличии и газа, и жидкости сложность вычислений резко возрастает. Это обусловлено тем, что изменение массы вещества в каждой из сред влияет на состояние другой, что приводит к трудности вычисления нового состояния модели в целом, а иногда и к невозможности это сделать: модель «раскачивается», и равновесное состояние не достигается.

·     Стоит упомянуть о проблеме, связанной с погрешностями вычислений. При вычислении нового состояния модели через время Dt может произойти несогласование значений. Например, суммарное значение массы вещества, поступившего в емкость по трубам согласно вычисленным значениям потоков в них, может не соответствовать вычисленному значению массы этого вещества в емкости; найденные значения потоков могут не быть точными решениями уравнений вида (3) и т.п. Опасность в том, что погрешности вычислений могут накапливаться, а это может привести к нарушению условия сохранения суммарного значения массы вещества в системе. Конечно, устранить подобные эффекты несложно, но все же они требуют отдельного внимания.

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

Два этапа вычислений

В целом выполнение шага моделирования по времени представим в виде итеративного процесса. На практике выяснилось, что этот шаг в общем случае может потребовать боль- ших вычислительных затрат, связанных с взаи- модействием уравнений, задающих модель. Чтобы снизить сложность вычислений в рамках этого процесса, каждую итерацию разобьем на два этапа: на одном этапе значения потоков в трубах будут вычисляться на основе состояний емкостей, на другом – состояний узлов. Чтобы это стало возможным, добавим дополнительное требование к модели: каждая труба соединяет узел и емкость. Другими словами, запрещается включать в систему трубы, соединяющие два узла или две емкости непосредственно.

Заметим, что это ограничение несущественно для соединения «емкость–емкость», так как такая труба моделируется двумя трубами и узлом между ними. Что касается соединения «узел–узел», то от него приходится отказаться; тем не менее, практика показала, что такое соединение во многих случаях может быть промоделировано двумя трубами и емкостью малого объема между ними.

Заметим также, что, поскольку каждая из труб теперь соединяет узел и емкость, для нее однозначно определено положительное на­правление (напомним: от узла к емкости), и тем самым вопрос о необходимости уточнять это направление отпадает.

Коллизии и их учет в модели

Рассмотрим ситуации, приводящие к возникновению таких состояний модели, которые невозможны с точки зрения физики описываемого процесса (например, отрицательная масса, бесконечное большое давление, переполнение объема емкости и т.п.). Конечно, лучший вариант – не допускать их возникновения, но, как уже отмечалось выше, для этого требуются непозволительно большие вычислительные ресурсы и не имеющий границ по точности вычислений компьютер.

Введем понятие коллизии.

Коллизией в емкости считаем ситуацию, когда в результате вычислений получена недопустимая масса вещества (отрицательная или превышающая вместимость емкости).

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

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

Корректирующее давление

Коллизии в трубах учесть несложно, чего нельзя сказать о коллизиях в емкостях. Можно было бы дробить шаг Dt, но, как отмечалось выше, это приводит к неприемлемо большим вычислительным затратам.

Эффективным инструментом для работы с коллизиями в емкостях оказался параметр, который был назван корректирующим давлением. Это дополнительное давление, которое способно остановить существующий поток.

Идея заключается в следующем: при возникновении коллизии в емкости вычисляется такое значение давления в ней, которое при неизменности значений остальных параметров не позволило бы этой коллизии возникнуть. Разность между найденным значением и вычисленным ранее давлением в этой емкости и считается корректирующим давлением. Добавление корректирующего давления (оно может быть как положительным, так и отрицательным) к найденному ранее меняет состояние модели, в результате чего либо соответствующая коллизия не возникает, либо ее «разрушающий» эффект значительно снижается.

Отметим, что корректирующее давление меняет вычисленное состояние модели и при этом нарушается соответствие между состоянием компьютерной модели и состоянием описываемой с ее помощью системы. Но это на- рушение кратковременное, зато сохраняется стабильность модели без привлечения существенных дополнительных вычислительных ресурсов.

Буфер в емкости

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

Для устранения возникшей трудности было принято решение позволить быть и тому, и другому, но лишь кратковременно и под строгим контролем. Технически это выражается введением параметра, положительное значение которого соответствует лишней массе веще- ства в емкости, а отрицательное – массе, изъятой из емкости сверх возможного.

Использование буфера в паре с корректирующим давлением оказалось довольно эффективным. Единственное, величина буфера была ограничена разумными значениями (10 % максимально возможной массы).

Сглаживание результатов вычислений

При компьютерной реализации модели, помимо описанных технических приемов, потребовались и некоторые другие, в частности, сглаживание результатов вычислений.

Иногда вычисления шага приводят к слишком резкому изменению значений потоков, что может вызвать «раскачивание» модели. Сглаживание состоит в том, что берется средневзвешенное значение между предыдущим и новым найденным значениями потока в трубе (например, среднее арифметическое), в результате чего достигается «плавность» изменений.

Дробление шага

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

Дробление шага происходит, когда при вычислении следующего состояния модели значение хотя бы одного из параметров выходит за допустимые границы (например, превышает по модулю возможности компьютерного представления чисел с плавающей точкой) и абсолютное значение буфера в емкости превышает некоторые разумные (заранее определенные) границы.

Дробление шага осуществляется рекурсивно в границах текущего шага. При этом шаг моделирования, на котором вычисляется состояние модели через время Dt, разбивается на два последовательных шага, на каждом из которых определяется состояние модели через время Dt/2.

Практика показала следующее: при начальном значении Dt = 50 мс дробление шага возникает не так часто и вложенность соответствующей рекурсии редко превышает величину 3. Тем не менее, при описанном подходе совсем без дробления шага обойтись не удалось.

Итеративные вычисления в узлах внутри шага моделирования

Опишем схемы итерационного вычисления нового состояния модели через время Dt неявным методом.

В результате уравнение (1) преобразуется к следующей разностной схеме:

                 (1а)

Величины G0 и GDt – это значения потока G (газа или жидкости) в некоторой трубе на начало и конец шага моделирования соответственно.

Во избежание нагромождения индексов будем опускать в формулах индекс Dt и вместо GDt писать просто G. Применяя к (1а) метод Ньютона, получаем формулу (1б), позволяющую находить G итеративно с некоторой заданной точностью. Значения величин G(0) и DP(0) определим позже.

На k‑й итерации значение DP для некоторой трубы вычисляется по формуле

                      (4)

где Pуз. – давление в узле, к которому присоединена труба; P – давление в емкости, к которой присоединена труба; HS – гидростатическое давление, создаваемое столбом жидкости. Верхний индекс в обозначении давления в емкости введен, чтобы подчеркнуть неизменность этой величины при соответствующих итеративных вычислениях.

Гидростатическое давление HS в (4) вычисляется по-разному для жидкости и газа.

Для жидкости гидростатическое давление HSж в трубе вычисляется по формуле

                        (5а)

где g – ускорение свободного падения; Dh1 –разность высот между точкой подсоединения трубы к емкости и точкой расположения узла; Dh2 – высота столба жидкости от точки подсоединения трубы к емкости до поверхности жидкости, причем Dh2 учитывается, если поверхность жидкости находится выше точки присоединения трубы к емкости (в этом случае b2 = 1), и не учитывается, если это не так (в этом случае b2 = 0).

Для газа гидростатическое давление HSг в трубе вычисляется по похожей формуле:

                      (5б)

в формуле используются те же величины, что и в (5а), при этом Dh1 учитывается, если на начало вычисления (по Dt) в трубе была жидкость (в этом случае b1 = 1), и не учитывается, если это не так (в этом случае b1 = 0).

Пусть i – идентификатор трубы (например, ее номер в системе). Тогда, учитывая формулы (4), (5а) и (5б), из (1б) получаем формулы (1в) и (1г), описывающие итеративный процесс вычисления в трубе i потоков жидкости и газа соответственно.

В формуле (4) присутствует величина Pуз., соответствующая давлению в узле. Значение Pуз. может быть найдено из уравнения (3), задающего условие баланса потоков. Поскольку для формул (1в) и (1г) соответствующее значение требуется вычислить как для жидкости, так и для газа, получаем две внешне похожие формулы (3а) и (3б), описывающие необходимые вычисления. В обеих формулах суммирование ведется по всем трубам, присоединенным к узлу, для которого вычисляется давление.

Чтобы прийти к уравнениям (3а) и (3б), достаточно в уравнение (3) подставить выражения для потоков (1в) и (1г) соответственно, а затем выразить требуемые величины.

Если на некотором шаге n достигнута нужная точность, то положим G* = G(n), P* = P(n), причем при любых возможных сочетаниях нижних индексов. Эти величины можно рассматривать как приближенные значения искомых величин, но они, вообще говоря, могут быть неокончательными и могут уточняться; как именно – будет описано далее.

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

Итак, моделирование уравнений (1) и (3) получено, обратимся к уравнению (2). В результате применения неявного метода уравнение (2) преобразуется к следующей разностной схеме:

                                    (2а)

Величины m0 и mDt – это значения массы m вещества (газа или жидкости) в некоторой емкости на начало и конец шага моделирования (по Dt) соответственно. Суммирование потоков в (2а) ведется по всем трубам, соединенным с этой емкостью. Далее будем опускать в формулах индекс Dt и вместо mDt писать просто m. Используя неявный метод Ньютона, получим формулы итеративного вычисления как массы в емкости, так и потоков в трубах, соединенных с этой емкостью.

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

                                  (2б)

где значения 𝑋(𝑘) и 𝑌(𝑘) определяются, соответственно, по формулам (6а) и (6б). Как при этом определяется m(0), уточним позже (и для жидкости, и для газа). Суммирование в (6а) и (6б) ведется по всем трубам, присоединенным к данной емкости.

Чтобы найти массу жидкости в емкости, требуется решить уравнение

                                  (2в)

где коэффициенты A, B и C находятся по формулам (7а), а присутствующая в формулах для вычисления B и C величина Z – по формуле (7б). В качестве значения m(k+1) нужно взять тот корень уравнения (2в), который ближе к m(k).

Поступление и расход вещества в емкости (то есть потоки в трубах) находим согласно формуле (1б), которая теперь вместо формул (1в) и (1г) преобразуется в формулы (1д) и (1е). Обращаем внимание, что в формулах (6а), (7б), (1д) и (1е) Piг и Piж обозначают давление газа и жидкости в узле, к которому присоединена i-я труба, а не в емкости, как в случае формул (1в) и (1г). Давление жидкости Pе.ж и давление газа Pе.г в емкости, присутствующие в (1д) и (1е), находятся по формуле

где R – универсальная газовая постоянная; Tг – температура газа (в градусах Кельвина).

Если на n-й итерации достигнута нужная точность, то полагаем m* = m(n), G* = G(n), P* = P(n), причем при любых сочетаниях нижних индексов для всех перечисленных величин. Эти величины можно рассматривать как приближенные значения искомых величин, но они еще не являются окончательными.

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

Итерации итераций

Покажем, как именно протекают два этапа вычислений.

Повторяется следующий процесс вычисления потоков вещества в системе:

-     итеративный процесс в узлах;

-     итеративный процесс в емкостях.

При первом итеративном процессе в узлах полагаем G(0) = G0, P(0) = P0 (при всех возможных сочетаниях нижних индексов), а при последующих – G(0) = G*, P(0) = P*, где значения G* и P* были получены на предыдущем итера- тивном процессе в емкостях (а для некоторых значений давления – в узлах).

При первом итеративном процессе в емкостях полагаем m(0) = m0, G(0) = G*, P(0) = P* (при всех возможных сочетаниях нижних индексов), а при последующих – m(0) = m*, G(0) = G*, P(0) = P*, где значения m*, G* и P* были получены ранее (на предыдущих итеративных процессах в узлах или в емкостях).

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

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

Моделирование буфера

Моделирование буфера для массы вещества в емкости является наиболее простым из упомянутых технических приемов.

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

Эти величины отвечают за закон сохранения массы вещества. Например, если по каким-то причинам в емкость нужно добавить жидкости больше, чем она может вместить, то излишки будут отправлены в bufж (но при первой же возможности они должны быть оттуда изъяты). При изъятии жидкости из емкости значение bufж может оказаться отрицательным. Величина bufг используется аналогично.

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

Если на (k+1)-й итерации при вычислении массы жидкости mж в емкости было получено отрицательное значение, то можно воспользоваться буфером. Но в физической реальности все происходит иначе, поэтому:

-     полагаем

-     вычисляем корректирующее давление  по формуле (8а);

-     далее в вычислениях вместо найденного для этой емкости давления  используем величину .

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

Условием того, что этой коллизии не стало, является следующее: Pкор.ж ³ 0.

Если на (k+1)-й итерации при вычислении массы жидкости mж в емкости получено значение, превышающее возможность емкости вместить такую массу, то делаем следующее:

-     полагаем ;

-     вычисляем корректирующее давление  по формуле (8б);

-     далее в вычислениях вместо найденного для этой емкости давления  используем величину .

Присутствующая в формуле (8б) величина  вычисляется по формуле

                       (9а)

В формуле (9а) величина P[max] соответствует максимальному давлению в емкости, которое допустимо в исходной моделируемой системе. Эта величина может быть как единой для всех емкостей системы, так и определяться для каждой емкости отдельно.

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

Условием того, что этой коллизии не стало, является следующее: Pкор.ж £ 0.

Коллизии и корректирующее давление можно устранять сразу же, как только были выполнены соответствующие условия. Но более практичным оказался другой подход: коллизии сохраняются до конца шага вычисления (по Dt), а в начале следующего шага безусловно отменяются (при необходимости они будут обнаружены и учтены снова).

Хотя корректирующее давление вычисляется для удержания значения массы жидкости в емкости в допустимых границах, вычисленная масса жидкости все-таки может выйти за эти границы из-за взаимодействий различных условий. В таком случае нужно прибегнуть к использованию буфера bufж. Если при этом значение, попавшее в bufж, превышает по модулю установленные разумные границы, то имеет смысл прибегнуть к дроблению шага, то есть разбить процесс вычисления состояния модели через время Dt на два последовательных шага, на каждом из которых вычисляется состояние модели через время Dt/2.

В случае газа корректирующее давление вычисляется и используется аналогично: по формулам (8в) и (8г). При этом величина  вычисляется по формуле, которую несложно получить из (9а), заменив в ней  на , а также  на :

                    (9б)

Кроме того, как и в случае с жидкостью, могут потребоваться использование величины bufг и дробление шага вычисления.

Коллизии в трубах и меры преодоления

Пусть на (k+1)-й итерации в некоторой трубе при вычислении потока вещества G сложилась следующая ситуация: был вычислен ненулевой поток, но по каким-либо причинам реализовать этот поток невозможно (например, из-за отсутствия нужного вещества в емкости в точке присоединения трубы к ней). Тогда поток этого вещества в трубе полагаем нулевым, то есть в зависимости от вещества  или

Возникшая коллизия и нулевой поток сохраняются до конца шага вычисления (по Dt). В начале следующего шага все коллизии в трубах отменяются без каких-либо условий и вычисления продолжаются в штатном режиме (при этом сохраняется возможность повторного возникновения какой-либо коллизии).

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

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

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

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

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

Этап перевода модели в новое состояние

В результате выполнения компьютерной программой итерационных вычислений будет найдено состояние системы, в которое она должна перейти через время Dt. Теперь нужно учесть, что значения всех параметров будут найдены лишь приближенно, в частности, при подстановке найденных значений потоков в уравнения вида (3), вообще говоря, ноль в левой части (3) может и не получиться, то есть может нарушиться сохранение значения суммарной массы вещества в модели.

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

Например, можно поступить следующим образом:

-     для каждого узла найти суммарное значение входящих и выходящих потоков;

-     найти модуль разности между этими значениями и его процент от наибольшего из них;

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

Получившиеся значения потоков вещества не приведут к изменению значения его суммарной массы. Теперь для получения нового состояния системы остается изменить массу жидкости и газа в каждой емкости в соответствии с найденными (скорректированными) значениями потоков.

Итоговые значения для m, G и P можно использовать в качестве начальных для вычисления следующего шага (по Dt).

Отметим, что при этом, возможно, придется задействовать величины bufж и bufг, но изменения их значений будут сравнимы с точностью вычислений.

Еще одно замечание. При описанном подходе значение суммарной массы вещества не будет меняться независимо от точности вычислений, поэтому точность можно уменьшить. Конечно, без необходимости делать это не стоит, но если, например, итеративные вычисления происходят слишком долго, а новое состояние модели нужно получить в реальном времени, можно пожертвовать точностью в пользу своевременности. Пользоваться такой возможностью или нет, зависит от требований, предъявляемых к модели.

Дополнительные возможности построенной модели

Итак, построенная модель позволяет вычислять состояния систем, в которых жидкость и газ могут перемещаться между емкостями по трубам (причем при приемлемых затратах времени, которыми можно еще и управлять). Но в реальных системах их состояние меняется за счет различного рода воздействий, таких как подача газа или жидкости под давлением, открытие или перекрытие труб и т.п.

Покажем, как подобные элементы управления можно встроить в построенную модель.

Клапан. Моделируется как труба, которая имеет два состояния, одно из которых соответ- ствует открытому клапану, другое – закрытому.

Невозвратный клапан. Работает в трех режимах: «открыт», «закрыт» и «невозврат». Первые два моделируются так же, как и в случае соответствующих состояний обычного клапана. В режиме «невозврат» оценивается перепад давления в начале шага, и в случае превышения некоторой заранее заданной уставки клапан открыт, иначе он закрыт.

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

Описанные элементы – клапан, невозвратный клапан и центробежный насос – в компьютерной модели могут быть представлены в виде трубы с дополнительными возможностями, что несложно сделать в объектно-ориентированных языках программирования за счет наследования (при условии, что уже имеется базовый класс, моделирующий трубу). Покажем возможность моделирования и более сложных объектов.

Для этой цели введем еще один тип элементов, который назовем устройством для перемещения вещества, или перемещателем. Перемещатель соединяет две емкости и может находиться в двух состояниях – «включен» и «выключен». В состоянии «включен» он может на шаге переместить определенное количество вещества из первой емкости во вторую, вычисленное согласно некоторому алгоритму (сам алгоритм не является предопределенным, и его можно уточнять в зависимости от задачи).

Редуктор. Этот элемент соединяет два узла и состоит из нескольких элементов: двух клапанов v1 и v2, двух емкостей t1 и t2, а также перемещателя sm. Они связаны между собой в последовательность v1 ↔ t1 → sm → t2 ↔ v2, где стрелками обозначены возможные направления перемещения вещества. С узлами соединяются свободные концы клапанов v1 и v2. Редуктор должен поддерживать в емкости t2 некоторое постоянное давление p0, если это возможно.

Алгоритм работы перемещателя sm в этом случае следующий.

Если давление в t1 превышает p0, а в t2 не превышает p0, то он переносит из t1 в t2 такую порцию вещества, которая обеспечивала бы в t2 давление p0. Если давление в t2 больше, чем p0, то переноса не происходит. Если давление и в t1, и в t2 меньше, чем p0, то редуктор ведет себя как труба. Клапаны v1 и v2 могут использоваться для внешнего управления редуктором; если этого не требуется, то их можно заменить трубами.

Поршневой насос. Этот элемент устроен аналогично редуктору, только вместо клапанов используются невозвратные клапаны (получаем схему v1 → t1 → sm → t2 → v2) и перемещатель работает по другому алгоритму: переносит порцию вещества в соответствии с заданной производительностью насоса и при условии, что значение напора, создаваемого насосом, достаточно для преодоления перепада давления между емкостями t2 и t1.

Отметим также, что при необходимости в систему можно добавить новые параметры, а некоторые параметры, присутствующие в приведенном описании модели в виде констант, можно заменить функциями. Одной из таких констант является коэффициент трения m, меняющийся динамически. Другая константа – это плотность жидкости. Известно, что плотность жидкости зависит от температуры, поэтому можно добавить параметр, отвечающий за температуру жидкости, а плотность r вычислять как функцию, зависящую от этого параметра. Можно изменить описания свойств веществ, позволив, например, жидкости сжиматься при увеличении давления (и тогда нужно добавить в описание и зависимость ее плотности от давления). Кроме того, можно учесть возможности, связанные с перемещением всей системы в пространстве (например, если она является частью машины, корабля, самолета и др.), введя функции, отвечающие за изменение координат, наклонов, вектора гравитации и т.п. – все это несложно встраивается в модель без существенного изменения лежащих в ее основе формул.

Ускорение вычислений

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

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

Тем не менее, на вычисление через 50 мс новых состояний некоторых систем (состоящих примерно из двух тысяч емкостей разных объемов, каждая из которых в среднем соединена с тремя-четырьмя узлами, содержит как жидкость, так и газ, имеет давление, заметно отличающееся от давления в соседних емкостях) требовалось в среднем около 100–300 мс при частоте процессора 3 ГГц и точности вычислений порядка 10–5. В каких-то сложных ситуациях на выполнение такого шага требо- валось несколько секунд, а в некоторых – порядка 10 мс. Число итераций могло превосходить тысячу, и его приходилось ограничивать.

Приведенные цифры показывают, что в случае вычислений в реальном времени имеются сложности, и решение вопроса об ускорении вычислений актуально. На практике значительное ускорение получилось за счет распараллеливания вычислений на процессоре, и в результате удалось добиться вычислений новых состояний описанной системы в реальном времени. В большинстве случаев для тех же данных на вычисление шага в 50 мс требовалось порядка 30 мс.

Отметим, что добиться ускорения вычислений можно попытаться с помощью другого численного метода. Авторы изначально ориентировались на неявный метод Ньютона как достаточно быстрый и надежный. Небольшой опыт применения в решении этой задачи других методов показал, что они заметно уступают ему по скорости и точности.

Тестирование

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

-     емкость от 1 л до 100 000 м3;

-     диаметр труб от 1 мм до 1 м;

-     количество емкостей от 2 до 2 000;

-     количество узлов от 1 до 2 000;

-     давление в емкости до 500 атм.;

-     разница давлений на концах трубы от 0 до 500 атм.;

-     шаг моделирования до 50 мс.

Практика показала, что увеличение длины шага моделирования нежелательно, хотя во многих случаях допустимо. Отметим также, что моделировались системы различных конфигураций, в том числе намеренно сложно устроенные. Например, моделирование систем с газом и жидкостью, в которых последовательно соединены трубами емкости небольшого объема (до 10 л) со значительным перепадом давления на концах получившейся цепи, показало, что модель вычисляет поток жидкости в этой цепи, величина которого почти не отличается от величины потока, идущего по трубе той же длины и того же диаметра, но без промежуточных емкостей; при этом число емкостей при тестировании достигало 100, а величина потока на шаге моделирования могла превосходить объем одной такой емкости на порядок.

Возможные перспективы развития модели

За границами рассмотрения остались следующие моменты.

Во-первых, для ускорения вычислений возможны параллельные вычисления с использованием ресурсов видеокарты. Если такое ускорение удастся получить, то, например, в качестве исходной величины Dt можно было бы взять меньшую, повысить точность вычислений, масштабировать время и т.п. Отметим, что для этого, видимо, придется решить ряд задач, связанных с реорганизацией структуры данных и вычислительной части модели таким образом, чтобы минимизировать обмен данными с видеокартой. Возможна также такая модернизация структуры модели, которая позволит использовать распределенные вычисления в компьютерной сети, суперкомпьютеры [6] или квантовые вычисления [7, 8]. У авторов имеются некоторые наработки в этом направле- нии, но их обсуждение выходит за рамки рассматриваемой задачи.

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

В-третьих, в некоторых системах возможно изменение агрегатного состояния вещества: жидкость может испариться, а газ – конденсироваться или раствориться в жидкости. Кроме того, некоторые вещества способны вступать в химические реакции друг с другом, образуя соединения, которые тоже могут находиться в разных агрегатных состояниях, причем сами реакции могут протекать с выделением или поглощением тепла, и т.п. Известно также, что при смешивании некоторых веществ объем получившейся смеси может отличаться от суммы объемов исходных веществ. Решение задач моделирования подобных процессов станет возможным, видимо, после решения задач, связанных с учетом состава вещества.

Заключение

Предложенная модель основана на системе эвристических соображений. И хотя она вполне работоспособна, для очень больших и сложно устроенных систем в случае недостаточно мелких значений шага Dt не стоит ожидать от нее результатов, очень близко соответствующих действительности, а в случае слишком мелких значений – быстроты и работы в реальном времени. Плюсом является то, что баланс между скоростью вычислений и их точностью (и адекватностью) можно регулировать.

Практика показала, что эта модель ведет себя вполне адекватно в отношении систем, в которых не возникает нагрузок на элементы, приводящих к опасности резкого (неплавного) переполнения или опустошения емкости. При наличии лишь одного компонента из двух рас- смотренных – либо жидкости, либо газа – мо- дель работает почти без задействования эври- стических модулей.

Авторы благодарны кандидату технических наук Александру Васильевичу Барулину за постановку задачи, лежащей в основе предложенной модели, и за возможность проведения исследований, связанных с ее решением, в рамках работы по одному из проектов, а также коллективу отдела моделирования сложных объектов и разработки тренажерных систем НИИ «Центрпрограммсистем» за обсуждения, связанные с этой задачей.

Литература

1.    Селезнев В.Е., Алешин В.В., Прялов С.Н. Математическое моделирование магистральных трубопроводных систем: дополнительные главы. М.: МАКС Пресс, 2009. 356 с.

2.    Селезнев В.Е., Алешин В.В., Прялов С.Н. Основы численного моделирования магистральных трубопроводов. М.: МАКС Пресс, 2009. 436 с.

3.    Селезнев В.Е., Прялов С.Н. Методы построения моделей течений в магистральных трубопроводах. М.: Едиториал УРСС, 2012. 560 с.

4.    Селезнев В.Е., Прялов С.Н. Численное моделирование течений в магистральных системах. М.: Едиториал УРСС, 2014. 800 с.

5.    Вержбицкий В.М. Основы численных методов. М., 2002. 840 с.

6.    Малявко А.А., Менжулин С.А. Суперкомпьютеры и системы. Мультипроцессоры. Новосибирск: Изд-во НГТУ, 2017. 64 с.

7.    Биллиг Ю. Квантовые вычисления. Начальный курс. М.: ИНТУИТ, 2020. 128 с.

8.    Китаев А., Шень А., Вялый М. Классические и квантовые вычисления. М., 1999. 192 с.

References

1. Seleznev V.E., Aleshin V.V., Pryalov S.N. Mathematical Modeling Trunk Pipeline Systems: Additional Chapters. Moscow, MAKS Press, 2009, 356 p. (in Russ.).
2. Seleznev V.E., Aleshin V.V., Pryalov S.N. Basis of Computational Modeling of Trunk Pipelines. Moscow, MAKS Press, 2009, 436 p. (in Russ.).
3. Seleznev V.E., Pryalov S.N. Methods of Construction of Flows in Trunk Pipelines. Moscow, Editorial URSS Publ., 2012, 560 p. (in Russ.).
4. Seleznev V.E., Pryalov S.N. Computational Modeling of Flows in Trunk Pipeline Systems. Moscow, Editorial URSS Publ., 2014, 800 p. (in Russ.).
5. Verzhbitsky V.M. Basis of Computational Methods. Moscow, 2002, 840 p. (in Russ.).
6. Malyavko A.A., Menzhulin S.A. Supercomputers and Systems. Multiprocessors. Novosibirsk, NGTU Publ., 2017, 64 p. (in Russ.).
7. Billig Yu. Quantum Computing for High School Students. 2018. 134 p.
8. Kitaev A., Shen A., Vyaliy M. Classical and Quantum Computing. Moscow, MCNMO, 1999, 192 p. (in Russ.).