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

ИСПОЛЬЗОВАНИЕ РАСПРЕДЕЛЕННЫХ СИСТЕМ ПРИ РЕШЕНИИ ЗАДАЧИ МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ВРЕДНЫХ ВЫБРОСОВ В АТМОСФЕРЕ

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

Ключевые слова: имитационное моделирование, распределенные системы, параллельные вычисления.

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

Необходимо рассчитать в трехмерном пространстве распространение выбросов в пространстве с течением времени.

— направление ветра по розе ветров;

— данные для расчета концентрации (количество вредного вещества, выбрасываемого в атмосферу и т.д.);

Величина концентрации во всех точках пространства в определенные моменты времени.

Концентрация по физическому смыслу является неотрицательной величиной, целесообразно использовать так называемые монотонные схемы, позволяющие получать неотрицательные решения. Воспользуемся методом расщепления по физическим процессам и на каждом малом интервале времени [tj, tj+1] длиной Δt [1].

Перенос примеси по траекториям:

где w_g описывает гравитационное оседание субстанций, u, v и w – скорости ветра по осям x, y, z соответственно. φ – искомая концентрация.

Этот этап является одним из основных в процессе переноса[2].

Для простоты изложения рассмотрим основной элемент схемы на примере одномерного уравнения:

Схема Куранта — Изаксона — Риса. Обобщение схем КИР на квазилинейный случай (при использовании дивергентной формы записи уравнений) очевидно:

где u – искомая функция концентрации, ƒ – функция описывающая источник дыма (его мощность выброса и т.д.), τ – шаг по времени, h – шаг в пространстве.

Рисунок 1 — Распространение выброса

Схема устойчива при выполнении условия Куранта:

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

Существует большое количество полуэмпирических и эмпирических формул, предложенных разными авторами для определения ΔH. По проведенным оценкам (Бем, 1971) наиболее приемлемо ?H определяется для нейтральной стратификации атмосферы по формулам Пристли, Спэра, Берлянда, Дановича–Зайгеля. В связи со сказанным для определения начального подъема ΔH газовой струи была выбрана следующая полуэмпирическая формула (Берлянд, 1975):

где T_в, T_a – соответственно температуры выбросов газа и окружающего воздуха по абсолютной шкале, w – начальная скорость выброса газов, R – радиус устья трубы, g – ускорение свободного падения, u – скорость ветра на высоте флюгера.

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

Учёт ландшафта создаётся вследствие введения такой системы координат:

где H – верхняя граница расчёта, ƒ(x,y) – функция высот ландшафта.

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

Проанализируем алгоритм для применения на распределенных системах.

Выполнять пока выполняется условие Куранта или загрязнение не достигло какой-либо точки:

1. Инициализация матрицы концентраций в пространстве φ[. ]

2. Расчет новой матрицы концентраций φ’[. ] для следующего момента времени t+Δt

3. Визуализация результата.

4. Присвоение φ[. ]=φ’[. ].

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

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

Рисунок 2 — Условная схема распределения вычислений по потокам

Рисунок 3 — Архитектура системы

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

1. В.А. Алексеев, А.В. Арефьев. Адаптивный экологический мониторинг окружающей среды // Экология и промышленность России. 2003. № 10. С. 11-13.

2. Алоян А.Е. Динамика и кинетика газовых примесей и аэрозолей в атмосфере. М.: ИВМ РАН, 2002. 201 с.

3. Корнеев В.Д. Параллельное программирование в MPI. Новосибирск: Изд-во СО РАН, 2000. 213 с.

Применение разностных схем расчёта гиперболических уравнений к решению некоторых задач конвективно-диффузионного переноса

Московский физико-технический институт

ПРИМЕНЕНИЕ РАЗНОСТНЫХ СХЕМ РАСЧЁТА

К РЕШЕНИЮ НЕКОТОРЫХ ЗАДАЧ

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

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

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

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

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

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

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

1) высокой экономичностью и простотой, присущими явным схемам;

2) устойчивостью и монотонностью при произвольном шаге интегрирования, которые характерны для неявных схем;

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

Идея решения задачи

Предлагаемый подход состоит в применении к решению уравнений вида

(1)

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

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

Аппроксимация уравнения переноса

Уравнение типа (1) для переноса вещества пуазейлевским течением жидкости в длинной узкой трубке было впервые выведено Тейлором ещё в 1953 г. [1] , но для ряда задач требуется его обобщить. Пусть профиль скорости в сосуде радиуса R описывается степенным законом , где u – среднемассовая скорость потока. Кроме того, предположим, что концентрация вещества также зависит от расстояния до оси сосуда по степенному закону , где – средняя концентрация.

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

При этих предположениях нетрудно получить дифференциальное уравнение переноса

(1/)

и одновременно его дискретное выражение

, (2)

аппроксимирующее (1/) при достаточно малых приращениях по времени T (по индексу n ) и по координате H ( по индексу m ) . Здесь

– коэффициент обгона потока веществом,

– число Куранта,

– характерное время диффузии вдоль радиуса сосуда ( – коэффициент этой диффузии),

– расстояние, на котором происходит существенное перемешивание вещества в радиальном направлении.

Если Т считать шагом по времени, а Н – шагом по пространству, то дискретное уравнение переноса (2) есть схема Куранта–Изаксона–Риса, используемая для решения гиперболического уравнения

, (3)

то есть уравнения (1/) без диффузии.

Второе дифференциальное приближение этой схемы

позволяет оценить её точность и по отношению к параболическому уравнению ( 1/ ) . Так, при (при стандартном параболическом профиле потока) аппроксимация имеет 3-й порядок, а при всех остальных отличных от 1 числах Куранта (то есть при существенной диффузии) – 2-й порядок.

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

, (4)

то есть уравнением (2) при произвольном шаге интегрирования , связанном с шагом сетки соотношением :

.

Условие 2-го порядка аппроксимации при любых можно обеспечить выбором шага по координате. Условие же 3-го порядка аппроксимации выполняется только в двух случаях – при , т. е. при , а также (приближённо) при , когда и уравнение становится почти гиперболическим.

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

Обобщение на многомерные и ветвящиеся сетки

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

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

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

При моделировании переноса в крупных сосудах становится несправедливым предположение о том, что радиальное выравнивание концентрации вещества в сосудах определяется поперечной диффузией : . С ростом радиуса R сосуда и линейной скорости потока u над диффузионным начинает преобладать конвективный механизм радиального перемешивания вещества. Если допустить, что радиальная конвекция обусловлена преимущественно разветвлениями сосудов, то характерное время этого процесса можно оценить как время прохождения вещества по одному сосуду ( L – длина сосуда).

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

,
где , , , , , , .

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

Результаты и их сравнение с другими методами

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

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

На всех рисунках представлена динамика распространения вещества вдоль последовательности артерий, “вырезанной” из ветвящегося сосудистого дерева человека. Геометрия сосудов реальная, модельным является лишь закон изменения граничных условий.

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

Время расчёта переноса во всех сосудах человека на 600 с