Параллельная реализация алгоритма для расчета двумерных нестационарных течений газа через пористые объекты с источниками тепловыделения Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 16, № 2, 2011

Параллельная реализация алгоритма для расчета

двумерных нестационарных течений газа через пористые объекты с источниками тепловыделения*

H.A. Луценко, Т.П. Мирошниченко, Д. С. Одякова, Д. И. Харитонов Институт автоматики и процессов управления ДВО РАН, Владивосток, Россия

e-mail: [email protected]

Разработан параллельный вычислительный алгоритм для исследования двумерных (плоских и осесимметричных) нестационарных течений газа через пористые объекты сложной конфигурации с источниками энерговыделения при известном перепаде давления газа на входе в пористый объект и на выходе из него в поле силы тяжести. Показано, что применение объектно-ориентированного подхода при распараллеливании позволяет разработать программный комплекс и автоматизировать задание геометрии пористого объекта без существенного изменения кода последовательной программы. На многопроцессорном вычислительном комплексе МВС-15К проведены расчеты двумерных нестационарных режимов газового охлаждения пористых тепловыделяющих объектов со ступенчато изменяющейся верхней частью сложной конфигурации, в которых возникают застойные зоны; проанализирована эффективность предложенного параллельного алгоритма.

Ключевые слова: параллельные алгоритмы, численное моделирование, пористые среды, газовое охлаждение, тепловыделение.

Введение

Многие природные и техногенные катастрофы приводят к возникновению очагов тепловыделения. Последние часто появляются в пористых средах, к которым с точки зрения механики можно отнести завалы разрушенных зданий, грунты, горные породы и т. д. Учитывая высокую потенциальную опасность таких объектов, необходима эффективная организация локализации и ликвидации самих очагов энерговыделения. Пример аварийного энергоблока Чернобыльской АЭС показал, что единственно приемлемым средством противодействия катастрофе может оказаться охлаждение подобных объектов проходящим через них газом (воздухом) [1]. Следует заметить, что изучение рассматриваемых течений газа в пористых средах может также быть полезным при разработке новых технологий и устройств в энергетике, металлургии, химической и добывающей промышленности.

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

*Работа выполнена при финансовой поддержке Президента РФ (грант № МК-2198.2011.1), РФФИ (гранты № 10-01-90004-Бел-а и 09-01-12056-офи-м), ДВО РАН (проекты № 09-1-09МПУ-04, 09-1-П2-05 и 11-Ш-В-01И-011).

Отличительной особенностью рассматриваемого в настоящей работе класса задач является открытость саморазогревающейся пористой среды в атмосферу снизу и сверху. Таким образом, граничные условия значительно отличаются от условий в классических задачах механики многокомпонентных сред, что затрудняет применение широко известных конечно-разностных схем. Поэтому для моделирования нестационарных одномерных течений газа через пористые объекты с выделением тепла в работах [2, 3] был предложен эффективный оригинальный численный метод, основанный на комбинации явных и неявных конечно-разностных схем. Этот метод был применен при исследовании различных нестационарных одномерных режимов охлаждения пористых тепловыделяющих объектов [2-4], В частности, в [3] обнаружено чрезвычайно опасное явление неограниченного разогрева пористого тепловыделяющего объекта при докритических краевых условиях, приводящее к его разрушению при условиях, удовлетворяющих критерию стационарности,

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

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

1. Математическая модель

Рассмотрим неподвижный пористый объект высотой Н, ограниченный с боков непроницаемыми нетеплопроводными стенками, открытый сверху и снизу (рис, 1), Пусть в результате протекания некоторых процессов в твердой фазе происходит выделение тепла, В нижнюю часть объекта под давлением подается холодный газ, который движется

Ф©ФФ©©©©©>© ФФФФФФФФФФ ФФФФФФ©©©© ФФФФФФФФФФ ФФФФФФФФФФ ФФФФФФФФФ©

ооооооооо© ффффффффф®

фффффффффф фффффффффф

ФФФе®оо©&Ф ©©©ФФФФФФ®

ГАЗ

Рис. 1. Пористый тепловыделяющий объект

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

Модель строится в рамках предположения двух взаимодействующих взаимопроникающих континуумов |8| и включает в себя уравнения энергии, движения, неразрывности и состояния дня каждого компонента (твердого и газообразного). Интенсивность межфазпого теплообмена предполагается пропорциональной разности фазовых температур в рассматриваемой точке среды, тепловыделение прямо пропорционально концентрации реагирующего вещества С, скорость убывания которой прямо пропорциональна самой концентрации. В уравнениях энергии компонентов учитывается не только теплопроводность, но и работа внутренних сил в газе. Дня описания динамики газа используется уравнение сохранения импульса для пористых сред, которое является корректнее классического уравнения Дарси и может применяться в большом диапазоне чисел Рейно.ньдса. Будем считать, что объем и масса конденсированной фазы изменяются незначительно и этими изменениями можно пренебречь, поэтому, так как твердая среда неподвижна и однородна, уравнения неразрывности и движения дня нее вырождаются. Учитывая данные работ |3, 4|, динамическая вязкость газа предполагается зависящей от температуры но формуле Сазерленда, а теплопроводность газа — пропорциональной его вязкости. Полагая справедливость уравнения состояния совершенного газа, запишем систему уравнений, моделирующую нестационарные течения газа через пористые тепловыделяющие объекты:

8Т

(1 - а) рсСс— = -а (Т - Тд) + (1 - а) С + (1 - а) АДТ,

аРср ( ^Г + <Уг ■ V) Тд) = а (Г - Тд) + а + (у8 • + а2^ + «V (АдЧТд) ,

\ р

— + • V) vg) I = -§га&р - pg - а-

3р ^ Т1'5 сР

+ У ■ (аР^ё) = 0. V = -~77Г 1 \ = -5-Ц,

Р (1 + X (1 - а)) ( ^Г + ' ) = -gradР - pg -

31 ' ' — ^ - - -1 с52 + Т/ Рг

8С

р = рКГд, — = -к2С. (1)

Здесь а — пористость; сс — теплоемкость конденсирован ной фазы; ср — теплоемкость газа при постоянном давлении; с^ и с^ — константы в формуле Сазерленда; § — ускорение силы тяжести; к\ — коэффициент проницаемости конденспрованной фазы; к2 — коэффициент, определяющий уменьшение тепловыделения; р — давление газа; Рг — число Прандтля для газа; _ константа, определяющая интенсивность тепловыделения; Н — газовая постоянная; £ — время; Т — температура конденсированной фазы; Тд — температура газа; х — эйлерова координата; — скорость газа; а — константа, определяющая интенсивность межфазового теплообмена; А — теплопроводность конденсированной фазы; Ад — теплопроводноеть газа; ^ — динамическая вязкость газа; р — плотность газа; рс — плотность копденсированной фазы; х ~ коэффициент присоединенной массы, учитывающий инерционное взаимодействие фаз при их ускоренном относительном движении.

Введем обозначение: и = avg — скорость фильтрации газа. Из последнего уравнения системы (1) найдем выражение для концентрации реагирующего вещества: С = ехр (—к2£). Заметим, что поскольку теплопроводность газа мала, то слагаемым, его содержащим, можно пренебречь. При исследовании плоских двумерных течений в декартовой системе координат условимся обозначать индексом 1 — горизонтальную координату, а индексом 2 — вертикальную. При наличии осевой симметрии удобно перейти к цилиндрическим координатам, тогда все параметры состояния трехмерной модели будут зависеть только от двух пространственных координат, В этом случае условимся обозначать индексом 1 — радиальную координату, а индексом 2 — аксиальную (осевую), На входе в пористый объект известны температура газа и давление, на выходе известно давление (истечение газа происходит в открытое пространство). Известны также условия теплообмена на входе в пористый объект, на выходе из него и на ограничивающих непроницаемых стенах. Отличительной особенностью модели является то, что расход и скорость фильтрации газа на входе в пористый объект не известны и должны определяться при решении задачи. Таким образом, краевые условия для системы (1) имеют следующий вид:

Здесь С — поверхность боковых степ, п — нормаль к поверхности боковых стен, в _ коэффициент теплоотдачи.

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

2. Вычислительный алгоритм однопроцессорной реализации

Система уравнений (1), моделирующая двумерное нестационарное движение газа через пористый тепловыделяющий объект, является нелинейной смешанной гиперболически-параболической системой уравнений и в общем случае не может быть решена аналити-

дТ

(2)

чески. Граничные условия в задачах данного класса существенно отличаются от условий в классических задачах механики многокомпонентных сред, что затрудняет применение здесь широко известных конечно-разностных схем. Поэтому для исследования двумерных нестационарных режимов охлаждения пористого объекта с выделением тепла при известном перепаде давлений на входе в объект и на выходе из него разработан оригинальный численный метод, основанный на комбинации явных и неявных конечно-разностных схем [5, 7], Уравнения энергии и сохранения импульса преобразовываются в явные конечно-разностные уравнения, из которых определяются температура газа, температура твердой среды и компоненты скорости фильтрации газа. Уравнение неразрывности преобразовывается в неявное конечно-разностное уравнение, из которого с учетом уравнения состояния совершенного газа находится давление газа. Плотность газа вычисляется из уравнения состояния совершенного газа.

Последовательная программа, производящая расчет плоских течений газа через пористые тепловыделяющие объекты, реализована следующим образом. Входные и выходные данные хранятся в шести трехмерных статических массивах, каждый из которых отвечает за значения одной из следующих величин: температуры газа, давления газа, плотности газа, температуры твердой фазы, вертикальной (или аксиальной) и горизонтальной (или радиальной) скорости фильтрации. Первые две координаты массивов отвечают за размерность сетки вычислений, а третья — за количество временных слоев, необходимых для организации вычислений, В подготовительной фазе вычисляются константы и считываются из соответствующих файлов начальные значения массивов, Далее следует вычислительная фаза, реализующая цикл по времени, в котором значения текущего временного слоя рассчитываются с использованием значений предыдущего слоя. Алгоритм нахождения данных величин на каждом временном слое следующий: фиксируем индекс "горизонтали"; последовательно меняя индекс "вертикали", решаем четыре явных конечно-разностных уравнения системы (из которых определяем температуру газа, температуру твердой среды и компоненты скорости фильтрации газа); затем решаем методом прогонки неявное конечно-разностное уравнение (из которого определяем давление газа) и далее тривиально находим плотность газа; после этого меняем индекс "горизонтали" и повторяем процесс. Следует заметить, что для различных горизонтальных координат алгоритмы прогонки могут различаться в зависимости от проницаемости соответствующей верхней части пористого объекта: если верхняя часть объекта открыта в атмосферу, то применяется обычная прогонка (далее — "прямая" прогонка), при закрытой верхней части прогонку необходимо тривиально изменить ("обратная" прогонка). После расчета внутренних значений всех параметров данного временного слоя вычисляются краевые условия. Затем происходит переход либо к следующему временному слою, либо к заключительной фазе, в которой реализуется запись результатов вычислений в семь файлов: значения параметров каждого из шести массивов записываются в соответствующий им файл, в седьмой файл размещается служебная информация вычислений.

Обратим внимание на то, что общая сложность вычислительного алгоритма оценивается как О ' гДе ^ ~~ шаг 110 пространству, при этом шаг по времени т = гЬ2. Однако эмпирическим путем было выяснено, что с уменьшением шага по пространству Ь, обычно г можно увеличивать, при этом сложность алгоритма достигает О ( —

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

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

3. Технология распараллеливания

Программная реализация параллельного алгоритма готовилась для кластера МВС-15К, состоящего из узлов с двумя процессорами PowerPC, объединенных коммуникационной средой Mvrinet, Особенности вычислительного алгоритма (в частности, "прямая" и "обратная" прогонки, сложная вложенность циклов) не позволяют эффективно применить автоматические и полуавтоматические средства распараллеливания. Анализ наиболее распространенных методов декомпозиции задачи на независимые компоненты [9] показал, что методы разделения индексов и смены циклов также не подходят для распараллеливания рассматриваемого вычислительного алгоритма из-за наличия в нем прогонок и сложных областей вычислений. Неэффективна в этом случае и функциональная декомпозиция, так как по сравнению с размерностью сетки вычислений и количеством итераций время, затраченное на вычисление значений параметров в одной ячейке сетки, невелико. Наиболее подходящим для данной задачи является метод геометрической декомпозиции, широко применяемый для вычислений на регулярных структурах. Используемый вычислительный алгоритм позволяет произвести геометрическую декомпозицию по горизонтальной (или радиальной) координате, декомпозиции же по вертикальной (осевой) координате препятствует наличие прогонок. Параллельная версия вычислительного алгоритма была построена "вручную" на основе исходной последовательной программы с использованием коммуникационной библиотеки MPI, так как особенности среды выполнения расчетов делают нецелесообразным применение технологии ОрепМР,1 При организации межпроцессорного взаимодействия использовались блокирующие процедуры MPI (MPI_Send и MPI_Recv).

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

1 Следует отметить, что технология MPI позволяет запускать программы не только на кластере, но и на SMP-архитектуре.

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

При преобразовании исходной программы к объектно-ориентированному представлению были выполнены следующие модификации:

— функции, отвечающие за вычисления, и функции, выполняющие чтение и запись данных, сгруппированы по классам, соответствующим различным типам вычислений ("ступенька" с открытым или закрытым верхом, краевые зоны);

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

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

Распараллеливание вычислений в программе происходит при помощи геометрической декомпозиции по горизонтальной координате т. Для обеспечения равномерной загрузки вычислительных узлов предварительно решается задача нахождения областей, сопоставимых друг с другом по вычислительной нагрузке. Для этого вычислительная область разбивается по координате т € [0..ттах) на к частей следующим образом: строится возрастающая последовательность из к + 1 натуральных чисел {т^ | 0 ^ ] ^ к, т0 = 0,тк = ттах} такая, что значениями ко ординат т € [mj-1,mj) определяется область вычислений для ^'-го узла. Каждый вычислительный процесс отвечает за расчет своей области, определяемой диапазоном значений координаты ш, обмениваясь на каждой итерации граничными значениями только с ближайшими соседними процессами, Характер вычислений в программе таков, что вычислительная нагрузка практически пропорциональна площади области, поэтому при построении последовательности для оценки нагрузки использовалась зависимость Б площади рассчитываемой области от координаты т. Учитывая, что вычислительная нагрузка ^'-го узла пропорциональна Б (mj) — Б (mj-1), а в идеале площадь каждой расчетной области должна равнятьеяБпееа = Б(ттах)/к, простой итерационный алгоритм позволяет вычислять

Фаза параллельного счета

Подготовительная фаза

Пересылка Г—Л /—Л /—Л /-—\ Пересылка

начальных -сю- • сх> вычисленных

значений ¿Обмен граничными ж \ значениями т значений

Пересылка Пересылка

начальных -СЮ- • ск> вычисленных

значений значений

Пересылка Пересылка

начальных <х> • сх> вычисленных

значений ¿Обмен граничными | \ значениями т /"—Л —N /"—Л /"—Л значений

Пересылка Пересылка

начальных -СЮ- • оо вычисленных

значений значений

Заключительная фаза

Рис. 2. Схема распределения нагрузки в параллельной программе

значения последовательности ш, зная значения с ] = 1 до ] = к из условия 5(ш, - 1) - 5(Ш-1) < Бпееа ^ 5(ш,) - 5(ш^).

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

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

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

4. Тестирование параллельной программы

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

I, ю3с ¿Л04С

для сеток с шагом Л = 1/40 (а) и Л = 1/100 (б)

была протестирована на расчете плоского и осесимметричного течения газа через ступенчато сужающиеся пористые тепловыделяющие объекты, исследованные подробно в [5, 7], Все параллельные вычисления производились на комплексе МВС-15К ЦКП ДВО РАН "Дальневосточный вычислительный ресурс". Сравнительные расчеты независимо от количества используемых узлов дали совпадающие с последовательной программой результаты. Графики зависимости времени работы параллельной программы от количества процессоров приведены на рис, 3,

Теоретическая функция времени работы программы от количества процессоров вы-

С С

ражается формулой Т = + — (sign(iV — 2) + 1), где N — количество процессоров,

N V к к

к — шаг по пространству, константа С отвечает за сложность расчета одной ячейки, константа С2 — за производительность коммуникационной среды. При N = 2 происходит односторонний обмен между узлами, при большем количестве узлов каждый узел производит обмен с двумя соседними за исключением крайних. Графики на рис, 3 совпадают с теоретически определенной формулой, что указывает на качественное распараллеливание программы2 и позволяет определить оптимальное количество узлов для рассматриваемой задачи. Из графиков видно, что оптимальное количество процессоров, после которого наращивание количества вычислительных узлов не имеет смысла, увеличивается с измельчением расчетной сетки. Для относительно грубой сетки с шагом по пространству к = 1/40 целесообразно использовать четыре процесса, а для шага к = 1/100 — соответственно восемь,

5. Газовое охлаждение пористого тепловыделяющего объекта с застойными зонами

С помощью разработанного параллельного программного комплекса исследован процесс течения газа через плоский пористый тепловыделяющий объект со ступенчато изменяющейся верхней частью, схематично изображенный на рис, 4, Для такого объекта

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

mm. m>mm\

©з© G©e©c@

» Si—ni

от »с®сф§ J1 ü

эоэ see®©©©©©©©

©со c©coc©c©gc©

©с© _ сосососоооэ mm я тёфффшш®

осо&ш®£&тсоосо г©эе@шв5>шоэшэ

аоэ@©€^с©©€)Соа©э

Рис. 4. Пористый тепловыделяющий объект

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

о кг 0 Дж о Дж о Дж

Рс = 2.2 • 10 —, сс = 9.2 • 10 ——, а = Ю3 ^ , сг, = 103^—, м3 кг - К м3 - К - с кг • К

с81 = 1.458 • 10"6-св2 = 110.4 К, к\ = Ю-8 м2, С}0 = 5 • 103 Дж

м • с

- л/к' " ' м3-с

Дж n 0 м 7 1 a ín Дж

Л = 1.2 ——-—, (/ = 9.8—, к'2 = Ю-' —, /3=10

м • К • с' ' с2' с' м2 -К ■ с'

Решение системы (1) (см. разде.:: 1) проводилось в обезразмереппом виде, при этом безразмерные переменные вводились следующим образом: X = x/H, í = t/t*, u = u/u*, P = v/v*-, P = P/P*-> T = T/T*, Tg = Tg/T*, где H = Юм, t* = 1c, u* = 1 м/с, T* = 300 K, p* = 105 Па, p* = 1.2 кг/M3, Далее везде тильда опущена, так как использованы только безразмерные переменные.

В задаче положено, что давление па входе в объект быстро возрастет до 1.5 и после этого останется постоянным. Тогда краевые условия получаются из (2) при p0 = 1.5, Tgo = 1, ph = 1.

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

На рис. 5 приведен пример распределения искомых величин но пористому объекту через t = 432 000 (что соответствует 120 ч размерного времени) после начала процесса. Дня лучшей визуализации здесь и далее графики давления и плотности газа развернуты на 180°. Как видно из рисунка, пористый объект эффективно охлаждается проходящим через пего газом везде, кроме двух застойных зон, расположенных но краям

Рис. 5. Распределение температуры твердой среды, температуры газа, давления и плотности газа, вертикальной и горизонтальной скорости фильтрации газа но пористому тепловыделяющему объекту

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

Н

0.8

0.6

0.4

0.2

О

Рис. 6. Поле скорости фильтрации в пористом тепловыделяющем объекте

из исс намного хуже, и следовательно, температура здесь растет быстрее. Максимальный нагрев наблюдается в самой верхней части лечзой застойной зоны вблизи ее правой стенки. Оказывается, что в застойных зонах более сильный нагрев происходит возле вертикальных стенок, которые располагаются ближе к той части объекта, из которых газ выходит в атмосферу. В правой застойной зоне это явление наблюдается еще более наглядно — в нравом верхнем углу зоны температура существенно меньше, чем в ловом. Последнее можно объяснить следующим образом: в застойных зонах, закрытых с трех сторон непроницаемыми стенками, возникают завихрения движений газа, при этом нагревающийся из-за теплообмена с твердой средой газ при своем продвижении отводит тепло все хуже, в силу чего температура становится выше у вертикальной стенки, расположенной ближе к выходу из объекта. Направления движения газа в пористом объекте наглядно демонстрирует рис. 6, где показано поло скорости фильтрации. Потоки газа, поднимающиеся вдоль крайних вертикальных боковых стен, завихряются в нижних частях застойных зон и но мере продвижения к центральной сквозной части объекта значительно ускоряются. Это приводит к возмущениям потоков газа, проходящих но средней части объекта, причем правая менее глубокая застойная зона оказывает большое влияние на движение газа в сквозной части.

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

Заключение

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

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

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

Список литературы

[1] Маслов В.П., Мясников В.П., Данилов В.Г. Математическое моделирование аварийного блока Чернобыльской АЭС. М.: Наука, 1987. 144 с.

[2] луценко H.A. Нестационарные режимы охлаждения пористого тепловыделяющего элемента // Матем. моделирование. 2005. Т. 17, № 3. С. 120-128.

[3] Левин В.А., Луценко H.A. Возникновение неустойчивых режимов охлаждения пористого тепловыделяющего элемента при докритических краевых условиях // Горение и плазмохимия. 2005. Т. 3, № 2. С. 81-90.

[4] Левин В.А., Луценко H.A. Течение газа через пористую тепловыделяющую среду при учете температурной зависимости вязкости газа // Инж.-физ. журн. 2006. Т. 79, № 1. С. 35-40.

[5] Левин В.А., Луценко H.A. Численное моделирование двумерных нестационарных течений газа через пористые тепловыделяющие элементы // Вычисл. технологии. 2006. Т. 11, № 6. С. 44-58.

[6] Левин В.А., Луценко H.A. Движение газа через пористые объекты с неравномерным локальным распределением источников тепловыделения // Теплофизика и аэромеханика. 2008. Т. 15, № 3. С. 407-417.

[7] Левин В.А., Луценко H.A. Нестационарные течения газа через осесимметричные пористые тепловыделяющие объекты // Матем. моделирование. 2010. Т. 22, № 3. С. 26-44.

[8] Нигматулин Р.И. Основы механики гетерогенных сред. М.: Наука, 1978. 336 с.

[9] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. СПб.: БХБ-Петербург, 2002. 608 с.

Поступила в редакцию 7 июля 2010 г., с доработки — 10 сентября 2010 г.