Содержание к диссертации
Введение
Глава 1. Двумерная модель распространения облака тяжелого газа . 31
1.1. Исходная система уравнений 31
1.2. Получение двумерной модели 31
1.3. Упрощенная модель 40
Глава 2. Численные методы решения уравнения Бюргерса . 43
2.1. О свойствах методов 43
2.1.1. Погрешность аппроксимации 43
2.1.2. Устойчивость 44
2.1.3. Диссипация, дисперсия 48
2.2. Явные схемы для невязкого уравнения Бюргерса 49
2.2.1. Метод Лакса 49
2.2.2. Метод Лакса-Вендроффа 51
2.2.3. Метод Мак-Кормака 53
2.2.4. Сводная таблица. Вывод 54
2.3. Неявные схемы для невязкого уравнения Бюргерса 56
2.3.1. Центрированный по времени неявный метод 56
2.3.2. Неявный метод Эйлера 59
2.3.3. Вывод 61
2.4. Обзор явных схем для вязкого уравнения Бюргерса 61
2.4.1. Метод Браиловской 61
2.4.2. Метод Аллена-Чена 62
2.4.3. Метод Лакса-Вендроффа 63
2.4.4. Метод Мак-Кормака 64
2.4.5. Вывод 65
Глава 3. Схема решения двумерной системы . 66
3.1. Описание метода расщепления по физическим процессам 66
3.2. Этапы решения двумерной системы 69
3.3. Разностный метод решения двумерной задачи 71
Глава 4. Анализ полученных результатов 77
4.1. Растекание жидкости по неровной поверхности 77
4.2. Расползание газа по неровной поверхности 84
4.3. Моделирование аварии 90
Заключение. 95
Литература
- Получение двумерной модели
- Явные схемы для невязкого уравнения Бюргерса
- Этапы решения двумерной системы
- Расползание газа по неровной поверхности
Введение к работе
1. Постановка задачи.
Диссертация посвящена математическому моделированию распространения облаков тяжелых газов и растекания жидкости, с учетом неровностей подстилающей поверхности (рельефа), трения о подстилающую поверхность и препятствия, смешения облака газа с окружающим воздухом при движении облака, а также теплообмена с подстилающей поверхностью и окружающим воздухом.
Учитывая объективные трудности трехмерного моделирования рассматриваемых течений, в работе построена двумерная модель гравитационного растекания жидкости и тяжелого газа в условиях орографической и термической неоднородности. Модель рассматривается для случая, когда облако имеет неоднородное распределение плотности, концентрации и температуры по высоте облака, что позволяет, в частности, естественным образом описать механизм разбавления газа в облаке окружающим воздухом при движении облака.
Двумерная модель построена методом осреднения по высоте исходных трехмерных уравнений газовой динамики. Полученная двумерная система уравнений выглядит следующим образом:
St ас ау 2 ас Sx 10
а де ду 2 ду *ку } ду ш v* vti
apwd ерша spvws _ P _ „ і ~a~ ~~ёГ+~^T " Fdw+FwH + F+ WU Q'
SpS Spa5 dpvS _ n
St Зс ду
SpcS SpucS dpved _ -, _ F ,
SphS dpahS dpvhS _ ,
1zeS0
где
s( „da
pvS— + — pvS
ax.
\ SxJ Sy^ Sy
Аналогично определяется F^
дсу дх ) фу ду j Аналогично определяются Fqc , Fq^ .
F
Полученная двумерная система уравнений решалась методом расщепления по физическим процессам, а возникающие при этом подсистемы решались разностным методом на двумерной прямоугольной сетке с использованием явных разностных схем. В качестве этапов решения выделялись подсистемы, учитывающие рельеф подстилающей поверхности, описывающие процессы адвекции, диффузии, трения, влияния источника.
С помощью этой модели, реализованной на ШМ PC, были промоделированы задачи о растекании жидкости и газа при разрушении резервуара цилиндрической формы на ровной, наклонной поверхностях и поверхности с препятствием; а так же реальная авария под Уфой (03.06.89 г.), связанная с разрывом трубопровода и распространением облака пропан -бутановой смеси на пересеченной местности.
С помощью построенной модели возможно решение широкого круга задач динамики жидкости и газа, экологии и промышленной безопасности, например, изучение природных катастроф, связанных с затоплением местности, технологических аварий и катастроф, связанных с разлитиями жидкостей (нефти, нефтепродуктов и т. д.) и распространением облаков токсических и горючих газов.
2. Обзор литературы. 2.1. Современное состояние вопроса об исследовании механизма рассеяния тяжелых газов.
Значительные количества выбрасываемых в приземную зону пограничного слоя атмосферы химически активных газов и аэрозолей создают во многих промышленных регионах напряженную экологическую ситуацию [9]. Для разработки эффективной стратегии решения задач промышленной безопасности и экологических проблем необходимо проведение большой серии исследований, которая должна включать в себя создание современных методов расчета нагрузки, создаваемой различными выбросами на окружающую среду и человека.
Одним из элементов системы методической поддержки программ в области экологии и промышленной безопасности являются модели распространения тяжелых газов вблизи поверхности земли в условиях, близких к реальным промышленным площадкам и жилым массивам.
Тяжелые газы представляют наибольшую опасность для природы и человека,
поскольку распространяются вдоль поверхности земли. Примеры аварий, произошедших на
химических предприятиях и трубопроводах и повлекших за собой тяжелые последствия :
авария 09.12.70 г. в Порт - Хадсоне (штат Миссури, США) [11]. После разрыва
трубопровода было разрушено 37 сооружений, а размер зоны разрушения составил 8 км
от эпицентра взрыва;
авария 03.06.89 г. под Уфой (Башкирия, СССР) [12]. При разрыве трубопровода вытекло несколько тысяч тонн пропан - бутановой смеси, зона поражения составила 2500000 м2, пострадали 623 человека, погибло 575 человек !!!
20.03.89 г. в Ионаве (Литва, СССР) в результате полного разрушения резервуара произошел выброс 7000 т жидкого аммиака [10]. Распространение паров аммиака зафиксировано на максимальном расстоянии 23 км. В результате аварии пострадали 57 человек, погибли 7 человек.
Эти опасные явления требуют для своего описания подробного исследования механизма процесса и разработки специальных моделей переноса примесей в окружающем пространстве.
В основе упомянутых выше промышленных аварий лежат общие механизмы их возникновения и развития. Так, нефтяные газы и многие токсические вещества, широко используемые в промышленности, например, хлор и аммиак, хранятся в сжиженном виде под давлением. В случае потери герметичности резервуара или трубопровода происходит разлитие и мгновенное испарение жидкости и образование облака топливно-воздушной смеси (облака ТВС) или токсического облака. Дальнейшее протекание таких аварий и их последствия определяются в основном механизмом распространения облака в условиях термической и орографической неоднородности. К числу параметров токсического облака, требуемых для определения зоны поражения и токсических нагрузок на человека, следует отнести конфигурацию облака, его смещение (дрейф) и концентрация газа в облаке.
Проблема экспериментального исследования физического механизма рассеяния тяжелых газов впервые возникла в семидесятые годы в связи с необходимостью практического решения задач промышленной безопасности [9]. В первой серии экспериментов основное внимание уделялось исследованию макроэффектов распространения облаков. Изучались размеры облаков и средние концентрации в них примесей тяжелых газов в зависимости от
мощности и характера выброса, состояния окружающей среды. В небольшом числе работ основное внимание было уделено изучению турбулентных характеристик процесса рассеяния. Измерялись флуктуации скорости газа, концентрации примесей и на этой основе оценивались значения коэффициентов переноса.
Для разрешения многих вопросов относительно механизма процессов переноса вещества и энергии в тяжелых газах были поставлены специальные эксперименты. Наиболее полная серия опытов по рассеянию- тяжелого газа была выполнена в Англии в 1982-1984 гг. на полигоне острова Торни-Айленд [8]. Исследовалось рассеяние больших объемов предварительно подогретой до температуры окружающего воздуха смеси фреона-12 и азота. Контейнер с газом объемом 2000 мЗ представлял собой разборный пластиковый цилиндр высотой 13 м. Его конструкция предусматривала возможность быстрого раскрытия стенок и практически мгновенного истечения смеси. Эксперименты финансировались группой промышленных компании топливно-энергетического направления различных стран мира.
Так как проведение натурных экспериментов в условиях реальных предприятий в большинстве случаев является невозможным, то основным инструментом исследования процесса рассеяния облаков тяжелых газов, испускаемых промышленным объектом, является метод математического моделирования. Он включает в себя анализ физики всех фаз явления, построение на этой основе математической модели общего процесса, разработку численных методов и алгоритма решения задачи, написание компьютерной программы и проведение вычислительных экспериментов. На сегодняшний день задача описания рассеяния облака тяжелого газа в условиях термической и орографической неоднородностей, как отмечается в обзорах [22], [24], является одной из наиболее актуальных в промышленной безопасности.
При численной реализации трехмерных моделей рассеяния тяжелых газов используются конечно-разностные методы, реже - метод конечных элементов [21], [23]. Так, например, в модели обтекания рельефа [18] для решения системы уравнений движения и сохранения массы использовался модифицированный вариант метода крупных частиц, а в полуэмпирической модели турбулентного стратифицированного течения [19] применялся метод дробных шагов с использованием схемы стабилизирующей поправки. В работе [20] трехмерное течение рассчитывалось на основе метода конечных объемов, реализованного при помощи явной схемы Мак - Кормака в консервативной форме. Для численного моделирования пространственных течений горячих масс газа (в качестве математической модели выбрана полная система уравнений Навье - Стокса для сжимаемого теплопроводного газа) использовалась конечно-разностная методика, основанная на явной 3-х шаговой схеме расщепления по физическим процессам [25]. Иногда используется описание диффузии в
переменных Лагранжа наряду с описанием конвекции в переменных Эйлера [23]. В работе [26] на основе явного метода Годунова построена схема для численного моделирования 3-х мерных нестационарных уравнений Навье - Стокса. Конвективный член дискретизируется явно, а вязкий - неявно.
В настоящее время известно достаточно большое количество численных методов решения уравнений Навье - Стокса, описывающих течение несжимаемой вязкой жидкости. Большая часть этих подходов была разработана применительно к системе уравнений относительно функции тока у/ и вихря ш [27]. Общим недостатком этих методов является
использование в том или ином виде граничного условия для вихря на твердой поверхности тела, которое отсутствует в физической постановке задачи. Наличие дополнительного итерационного процесса, связанного с указанным граничным условием для вихря, лимитирует скорость сходимости численных алгоритмов. Следовательно, высказывается предположение о том, что разностная схема, позволяющая рассчитывать течения вязкой несжимаемой жидкости без использования граничных условий для вихря на твердой поверхности, при всех прочих равных условиях обладает большей эффективностью. Очевидная ограниченность методов решения (у/, бт) -системы связана с невозможностью
развития их на случай исследования пространственных потоков вязкой жидкости, а также течений сжимаемого газа, объясняет возросший в последнее время интерес к численному решению уравнений Навье -Стокса, записанных в естественных переменных.
По данным [22] к 1987 г. в мире существовало более сорока специализированных программ расчета аварийного рассеяния в трехмерном пространстве. Основной трудностью при численном моделировании процесса распространения примеси в атмосфере является ограниченная емкость памяти и скорость современных вычислительных машин. Так, например, для расчета реального процесса, протекающего одну секунду, требуется от 10 до 100 сек. вычислительного времени на машине CRAY-1.
Рассмотрим еще один важный момент. Одной из главных задач вычислительных методов расчета выброса в окружающую среду вредных химических отходов предприятий является расчет полей температур и концентраций, вызываемых такими выбросами в атмосфере и поверхностных водах [13]. Распространение загрязняющих веществ в жидких средах определяется двумя процессами: конвективным переносом вследствие осредненного движения среды и диффузией за счет турбулентности. Поэтому математическая модель должна правильно описывать как поле средних скоростей, так и характеристики турбулентной диффузии.
Система точных уравнений, описывающих во времени все детали эволюции поля скорости и скалярных полей в практической задаче не может быть решена с помощью современных вычислительных средств. В настоящее время существует единственный экономически оправданный выход: решать уравнения осредненного движения, которыми определяется распределение осредненных по времени величин [13]. При этом время осреднения должно быть много больше временного масштаба турбулентности, но много меньше временного масштаба осредненного течения (например, суточного цикла в пограничном слое атмосферы). Обычно только средние величины и имеют практический смысл. Уравнения осредненного движения содержат члены турбулентного переноса. Для замыкания системы уравнений эти члены должны быть аппроксимированы с помощью определенной модели турбулентности. В некоторых случаях вполне достаточно лишь приближенного описания характера турбулентности. Так, в задачах о больших водных массах, значения турбулентной вязкости и коэффициент турбулентной диффузии обычно принимают const. Более сложные модели в таких задачах себя не оправдывают из-за значительной неопределенности в задании граничных условий и погрешности в численном решении. Для других задач, в частности, при расчетах полей вблизи источников, требуются уже более точные модели. Сложность модели турбулентности зависит от условий конкретной задачи.
Во многих случаях в потоках со свободными поверхностями средние значения характеристик потока слабо меняются в вертикальном направлении, так что изменение этих величин в горизонтальном направлении можно определить, решая двумерные уравнения среднего течения для осредненных по глубине значений [13]. Уравнения среднего течения для осредненных по глубине величин можно вывести, интегрируя 3-х мерные уравнения по глубине в предположении, что распределение давления является гидростатическим. Такая модель была построена для жидкости, однако газ автором модели не рассматривался.
К недостаткам известных численных моделей рассеяния тяжелых примесей следует отнести излишнюю идеализацию граничных условий на поверхности земли [9]. Обычно она полагалась гладкой, что не соответствует наиболее важным случаям реальной действительности. Это говорит лишь о начальной стадии исследования явления с помощью современных математических и вычислительных средств.
2.2. Общий обзор численных методов решения задач гидромеханики.
В последние годы резко возросла роль и значение нового подхода к решению задач гидродинамики и теплообмена, который получил название вычислительной гидромеханики,
Основной особенностью этого вычислительного (или численного) подхода является то, что уравнения (чаще всего уравнения в частных производных), описывающие интересующий нас физический процесс, решаются численно. В гидромеханике обычно приходится решать нелинейные задачи, так как давление, плотность, температура и скорость должны быть определены из решения нелинейной системы уравнений в частных производных.
2.2.1. Уравнение Бюргерса (невязкое течение).
Рассмотрим одно простое нелинейное уравнение, аналогичное уравнениям гидромеханики [2]. Это уравнение должно включать в себя члены, описывающие те же физические процессы, что и члены, входящие в уравнение гидромеханики, то есть конвективный, диффузионный и нестационарный члены. Такое простое нелинейное уравнение было предложено Бюргерсом. Оно имеет вид
ckt ди д и
-- + ^- = ^-^. (1)
Ф ^ дх1
Первое и второе слагаемые в левой части этого уравнения являются соответственно
нестационарным и конвективным членами, а в правой части стоит вязкий член. Если вязкий
член не равен нулю, то уравнение (1) параболическое; если же он равен нулю, то в уравнении
остаются лишь нестационарный и нелинейный конвективный члены. Такое уравнение
гиперболическое и имеет вид
( <&
Х + «^Г = 0. (2)
ду (к
Его можно рассматривать как модельное для уравнений Эйлера, описывающих движение
идеального газа. Уравнение (2) можно интерпретировать и как нелинейное волновое
уравнение, при этом скорость распространения волны в различных точках будет разной. Так
как скорость распространения возмущений меняется, то характеристики начинают
пересекаться и в решении возникают разрывы, аналогичные ударным волнам в газовой
динамике. Следовательно, рассматриваемое одномерное модельное уравнение позволяет
изучать свойства разрывных решений.
Запишем уравнение (2) в дивергентном виде
ди dF -j
- + - = 0, F = J,2 (3)
В общем случае и неизвестная и, и функция F(u)-векторы.
Физическая интерпретация этого определения [4]; каждая компонента и является величиной, отнесенной к единице объема. Если слой заключен между х^ и х2, то полное
количество этой величины на единицу площади слоя равно интегралу от Хг до Х2от этй компоненты. Но
d *2
— \udx = ~F(ui + F(u)\
dt '>Х=Х2 Х У| Х=Х{
х1 Следовательно, соответствующая компонента F(u)\ есть поток этой величины через
единичную площадку в точке х в положительном направлении оси х. Перепишем уравнение (3) в виде Ш fa
где в общем случае А = А(и)- матрица Якоби cFj і di,-, а в нашем примере A =dF I du. Так
как наше уравнение в частных производных гиперболическое, то все собственные значения матрицы А вещественные.
Нелинейные гиперболические уравнения в частных производных обладают, согласно Лаксу, решениями двух типов [2]. Гладким называется такое решение уравнения (4), когда функция и непрерывна внутри области, а ее производная может иметь разрыв на границе (то есть решение уравнения является непрерывным по Липшицу). Слабым называется решение уравнения (4) гладкое всюду, кроме некоторой поверхности в пространстве (x,t), на которой функция и может иметь разрыв. На величину скачка функции и при переходе через поверхность разрыва накладываются определенные ограничения. Пусть w- произвольная непрерывная векторная функция с непрерывной первой производной, равная нулю вне некоторой ограниченной области, тогда и называется слабым решением уравнения (3), если
JJ (wt u + wx F)dxdt + J w(x ,Q)
причем
всякое непрерывное слабое решение является гладким. Одним из примеров слабого решения
являются ударные волны, возникающие в сверхзвуковых течениях невязкой жидкости.
Найдем необходимые условия существования слабого решения уравнения Бюргерса,
то есть необходимые условия существования решения с разрывом, как изображено на рис. 1.
траектория
движения
разрыва
Рисі. Задача о распространении разрыва для уравнений Бюргерса.
Пусть w(x,t) -произвольная непрерывная функция, имеющая непрерывную первую
производную. Кроме того, пусть она обращается в нуль на границе В области D и вне D (на дополнении к D); D -произвольная прямоугольная область в плоскости (x,t) . Ясно, что
в этом случае
(6)
(7)
Если функции и и F непрерывны и имеют непрерывные первые производные, то уравнения (6) и (7) эквивалентны. Входящий в уравнение (5) второй интеграл в последних двух уравнениях отсутствует, так как функция w на границе равна нулю. Функция u(x,t), удовлетворяющая условию (7) для любой функции w, называется слабым решением невязкого уравнения Бюргерса.
Рассмотрим случай, когда прямоугольная область D в плоскости (x,t) разделена кривой r(x,t) - 0, на которой функция и имеет разрыв (см, рис.2).
w=0
Рисі Схематическое изображение произвольной области с расположенным в ней разрывом.
Предположим, что функция и непрерывна и имеет непрерывные первые производные в подобластях, лежащих слева от z{D{) и справа от r(D2). Используя формулы интегрирования по частям и учитывая, что функция и равна нулю на границе области D и вне D, из (7) получим
jjj ^ + —)w(x,t)dxctt + jj| — + — \v(x,t)dxdt + \І[и\со$щ + [F]co$a2)ds = 0. (8)
Последний интеграл вычисляется вдоль кривой г(х,г) = 0, разделяющей подобласти Z)] и D2. Он появляется при интегрировании по частям вследствие того, что кривая r(x,t)~0 является границей подобластей D\и D2. Квадратными скобками обозначена разность значений заключенной в них величины по- разные стороны от разрьша ("скачок" этой величины при переходе через разрыв), щ и а2 -углы между направлением нормали к кривой r(x,/) = 0 и осями /их соответственно.
Согласно (6), входящие в (8) интегралы по подобластям Z)j и D2 равны нулю, поэтому и последний интеграл равен нулю для любой функции w. Следовательно,
[w]cOS«i + (^0)80 =0. (9)
Последнее соотношение и является условием, которому должно удовлетворять слабое решение и уравнения Бюргерса.
Рассмотрим движущийся разрыв. Пусть начальное распределение и(х,0) имеет вид,
показанный на рис.1, где щ и и2_значения и слева и справа от разрыва. В одномерном случае уравнение поверхности r(x,t) = 0 можно представить в виде /-/j(r)=0. Тогда входящие в (9) направляющие косинусы определяются соотношениями
cos«l=r „цд ' cosa2=-- :~ш
1+>12
1+/J2
(штрихом обозначено дифференцирование по х ). Следовательно,
1/2
l + i2
1 + г'2
И=> или И2-«1=^—^
Окончательно
dx щ+и2 ~dt^ 2
(10)
то есть скорость распространения разрыва равна полусумме скоростей слева и справа от него. Зная, что по обе стороны разрьша скорости постоянны и что сам он движется с
постоянной скоростью і.щ+и2)/2 , легко провести сравнение решений по различным
численным методам расчета течений с разрывами с точным решением.
Волны разрежения встречаются в сверхзвуковых течениях не реже, чем ударные волны. Известно точное решение уравнения Бюргерса, описывающее волну разрежения. Пусть начальное распределение и(х,0) имеет вид, изображенный на рис.3.
u(s,0)
u=l
u=0
Рис- 3. Начальные условия для волны разрежения.
Характеристики уравнения Бюргерса описываются соотношением
dx и
Вид характеристик в плоскости (x,t) показан на рис.4.
(")
Рис. 4. Характеристики для случая центрированной волны разрежения.
В левой полуплоскости характеристики - вертикальные прямые, а справа от характеристики, ограничивающей волны разрежения, они составляют с осью х угол я7 4 рад. Рассматриваемая задача похожа на задачу о распространении центрированной волны разрежения в течении сжимаемой жидкости, В случае уравнения Бюргерса волна разрежения ограничена слева линией х = О, а справа - проходящей через начало координат характеристикой, которая изображена на рис.4 точечной линией. Математически решение задачи о распространении волны разрежения можно записать в виде и = О, х < О,
u = xlt, 0<х
и = 1, x>t. Итак, заданное начальное распределение и приводит к образованию центрированной волны разрежения, ширина которой растет по времени линейно.
Описанные выше две задачи, часто встречающиеся в сверхзвуковых
газодинамических течениях - ударные волны и волны разрежения можно моделировать при помощи уравнения Бюргерса.
2.2.2. Уравнение Бюргерса (вязкое течение).
Полное нелинейное уравнение Бюргерса
(Ги
- + !/-
<% ' дс дх1
является параболическим уравнением в частных производных [1,2]. Оно используется как модельное для уравнений пограничного слоя, "параболизованных" уравнений Навье-Стокса и полных уравнений Навье-Стокса.
Как и в случае невязкого уравнения, для вязкого уравнения Бюргерса существуют точные аналитические решения при некоторых начальных и граничных условиях. Эти решения полезно использовать для сравнения различных разностных схем. Линейное и нелинейное вязкие уравнения Бюргерса можно скомбинировать в обобщенное уравнение
ut+(c+bu)ux = рихх, (13)
где с и Ъ - свободные параметры. При Ь = 0 получаем линейное уравнение Бюргерса, а при с-0 и Ъ~\ - нелинейное уравнение Бюргерса. Если с-У* и * = ~ 1, то обобщенное уравнение Бюргерса имеет точное стационарное решение, то есть lim и(х, і) при t->ao
1+/А-
(14)
с(х-х0)
и- —
которое при р. = Уд показано на рис.5.
0.5
-і 1—&
2 3 Mi
Рис.5. Точное решение уравнения (14)
Следовательно, если начальное распределение и задано соотношением (14), то точное решение не меняется по времени, а остается равным заданному начальному распределению. Если ц = 0, то получаем
гиперболическое уравнение. Как говорилось выше, уравнение (15) допускает разрывы
(слабое решение) [16]. Если эти разрывы удовлетворяют неравенству, называемому условием
энтропии, то слабое решение задачи с начальными условиями для - со < х < со является
единственным. Более того, доказано, что если // —> 0, то решение задачи (1) с некоторыми
начальными условиями стремится к слабому решению уравнения (15) с теми же начальными
условиями. Поэтому, если v мало, то решение уравнения (1) будет иметь тот же самый вид,
что и решение уравнения (15) вне области, примыкающей к границе. Например, в случае
plu сжимаемой вязкой жидкости, когда число Рейнолъдса Re велико (Re = —, где /-
И-
характерный линейный размер в задаче), в основном потоке могут существовать узкие области с очень большими градиентами. Если жидкость невязкая, эти области будут представлять собой ударные волны, то есть разрывы решения. Если же вязкость мала, но не равна 0, то толщина области с большими градиентами настолько мала, что не может быть описана с помощью конечно-разностной сетки и численное решение для вязкой жидкости ведет себя подобно решению для идеальной жидкости. Поэтому численная схема, предназначенная для решения уравнения (1) с малыми значениями /л должна быть приспособлена к расчету разрывных решений.
2.2.3. О свойствах методов.
В этом параграфе изложены основные понятия и методы, используемые при конечно-разностном решении уравнений в частных производных. Основой метода конечных разностей является дискретизация - замена непрерывной области совокупностью изолированных точек (сеткой), причем решение уравнений ищется лишь в этих точках (узлы сетки). Производные аппроксимируются конечными разностями и решение уравнения в частных производных сводится к решению системы алгебраических уравнений. Также рассматривается вопрос о том, сколь точно решение разностных уравнений приближается к решению исходной задачи. Для этого анализируется погрешность аппроксимации, устойчивость и согласованность разностных схем.
А) Погрешность аппроксимации.
Для каждого уравнения в частных производных существует множество его конечно-разностных аналогов, из которых обычно нельзя выбрать наилучший со всех точек зрения. В первую очередь при использовании метода конечных разностей надо стремиться к правильной аппроксимации уравнений поставленной задачи, а во вторую очередь выбрать "наилучшую" схему, то есть оптимизировать ее, учитывая ее точность, удобство программной реализации на ЭВМ и т.д.
Погрешностью аппроксимации называется разность значений частной производной и ее конечно-разностного аналога [2,5,6]. Можно характеризовать погрешность аппроксимации стандартным математическим обозначением порядка малой величины (О). Представление погрешности аппроксимации в виде, например, с(Дг) обозначает, что погрешность аппроксимации по абсолютной величине не превосходит К\Ах\ при Ах —»0 (для достаточно малых Ах), причем К > 0 -вещественная константа. В общем случае выражение /(х) = <\4{х)) означает, что существует такая не зависящая от х константа К,
что |/(jc)| < АТ|(0(х)[ для всех х из области S. В качестве области S часто выбирается область х —> да (достаточно большие х) или, как обычно бывает при использовании разностных методов, область х -> 0 (достаточно малые х ).
Б) Погрешность округления, точность.
Любое численно полученное решение, даже точное аналитическое решение уравнения в частных производных, зависит от ошибок округления, связанных с конечным числом знаков, используемых при арифметических операциях на ЭВМ. Возникающая при этом погрешность называется погрешностью округления. Она может оказать существенное влияние на решение конечно-разностных уравнений, так как получение этого решения обычно связано с выполнением большого числа однотипных арифметических операций. В ряде случаев погрешность округления пропорциональна числу узлов разностной сетки, поэтому измельчение сетки, снижая погрешность аппроксимации, может увеличивать погрешность округления. Следовательно, погрешность полученного на ЭВМ решения уравнения в частных производных равна сумме погрешностей аппроксимации и округления. Точность численного решения уравнения в частных производных определяется погрешностью аппроксимации не только самого уравнения, но и граничных условий.
В) Устойчивость, СХОДИМОСТЬ.
Согласованной называется разностная схема, аппроксимирующая уравнение в частных производных [2]. Условием согласованности разностной схемы является стремление к нулю погрешности аппроксимации при измельчении сетки. Это условие безусловно выполняется , если погрешность аппроксимации убывает при измельчении сетки, то есть если погрешность аппроксимации имеет порядок о(Дг), о(Дг) и т.д. Однако, если порядок погрешности аппроксимации равен, например, o(At I Ах), то схема будет согласованной лишь в том случае, когда измельчение сетки проводится в соответствии с условием Д//Дг->0.
Устойчивой называется разностная схема, если на каждом шаге по времени любая ошибка (погрешность округления, погрешность аппроксимации, просто ошибка) не возрастает при переходе от одного шага к другому [2].
Выполнения условий устойчивости и согласованности достаточно для сходимости разностной схемы. Под сходимостью в данном случае понимается стремление решения конечно-разностного аналога уравнения в частных производных к решению исходного уравнения (для одинаковых начальных и граничных условий) при измельчении сетки.
Для анализа устойчивости конечно-разностной схемы применяют, например, метод Фурье [2, 6, 14]. Метод предложен Дж. Фон Нейманом в Лос-Аламосе в 1944 г. Краткое описание метода впервые появилось в работе Кранка и Николсона (1947) [14]. В этом методе решение модельного уравнения представляется рядом Фурье с конечным числом членов и устойчивость (или неустойчивость) определяется тем, что каждое отдельное колебание затухает (или нарастает). Погрешность округления не будет возрастать на каждом шаге по времени, если коэффициент перехода не превосходит единицы. Выражение для коэффициента перехода можно получить при помощи подстановки разложения погрешности в ряд Фурье в конечно-разностное уравнение. Иллюстрация метода Фурье на примере модельного уравнения представлена в главе 2.
Рассмотрим прием, который позволит сделать некоторые выводы об устойчивости уравнений для сжимаемого газа на основании модельного уравнения (1) [14]. Типичным условием устойчивости для этого уравнения является ограничение на число Куранта:
С = \иу--<\ (16)
Это неравенство следует из рассмотрения невязких членов в уравнении (1), то есть при
,0 = 0. В уравнении (1) при /1 = 0 информация переносится с конвективной скоростью
непрерывной среды. На каждом расчетном шаге по времени возмущение в / -й точке влияет
на новое значение в (i+J) -й точке. Это означает, что за каждый шаг по времени Д/ информация переносится на расстояние Д* и таким образом вычислительная скорость распространения информации равна Дх/Д/. Неравенство (14) означает, что скорость распространения информации в непрерывной среде не должна превышать вычислительную скорость распространения информации, то есть
і і д*
В течениях сжимаемой жидкости члены с grad давления меняют скорость распространения информации в среде; она равна не и, несколько больше. Рихтмайер и Мортон словесно описали путь получения зависимости скорости распространения информации в среде от давления.
Ограничимся элементарным газодинамическим соотношением. Малое возмущение давления распространяется с местной скоростью звука а относительно газа, который сам движется со скоростью и. Возмущения давления распространяются во всех направлениях и необходимо рассматривать только О. > О. Таким образом, скорость распространения
информации в сжимаемой жидкости равна |и| + а. Следовательно, ограничение на число
Куранта запишется в виде:
с - (їм + а) < 1.
Переформулированное в других терминах условие устойчивости КФЛ констатирует, что область влияния конечно - разностных уравнений должна быть по меньшей мере столь же велика, как область влияния исходных ДУ. Физически это означает, что за один шаг по времени звуковая волна не должна проходить расстояние, большее размера ячейки.
Этот метод исследования устойчивости часто приводит к тем же результатам, что и строгий метод исследования устойчивости, и дает по крайней мере необходимое условие устойчивости. Однако тот метод не подходит для неявных схем.
В случае течения сжимаемой жидкости размерность задачи также влияет на условие устойчивости. При Дх = Ду применение одномерного метода одновременно для всех измерений обычно меняет условие устойчивости следующим образом:
С< -- ИЛИ С< —ft,
л/2 л/3
соответственно для 2-х или 3-х пространственных координат. Использование схемы расщепления по времени Марчука с раздельным расчетом вкладов от каждого измерения приводит к одномерному условию устойчивости.
Г) Диссипация, дисперсия.
Коэффициент перехода для некоторой конечно-разностной схемы зависит от шагов сетки и волнового числа. Например, для конечно-разностной схемы Лакса для волнового уравнения коэффициент перехода имеет вид [2]
G = cos-ivsin = |G|e^, где ф - фазовый угол, определяемый соотношением
Im(G) -vsin/ї / а
ф=агс ^іеМ="* tg-^T=агс '*- «т-
Из последнего соотношения ясно, как коэффициент G зависит от числа Куранта v и параметра частоты J3.
Чтобы показать связь коэффициента перехода и дифференциального приближения разностной схемы, запишем волновое уравнение в виде [15]:
щ+сих = ^
<х> ( ~3.п.. д2я+гЛ
Здесь Сіп и С-іп+і" коэффициенты перед членами уравнения с производными четного и
нечетного порядков по пространству.
Импульс или сигнал произвольной формы может быть разложением в интеграл или
—i((t]t~k х\
ряд Фурье представлен как суперпозиция гармонических волн вида е * \ где оз-
частота волны, кт = 2л/к- волновое число, Л- длина волны. Скорость, с которой фаза
волны ктх - Ш перемещается в пространстве, называется фазовой скоростью волны и равна
v = co/km. Зависимость фазовой скорости распространения волны от длины волны (или, что
dv то же, от волнового числа) называется дисперсией волны. Если — > 0, т. е. волны большей
длины распространяются с большей фазовой скоростью, дисперсия называется нормальной,
если — < 0, т. е. волны большей длины распространяются с меньшей фазовой скоростью, -dX
~ ,-r dv _
аномальной. При — = 0 дисперсия отсутствует. В случае дисперсии гармонические
составляющие сигнала смещаются относительно друг друга, в результате чего профиль
сигнала искажается.
Образ Фурье оператора шага разностной схемы равен
G=eAt\ где Ats = -ikmcAt + AtJ^(ikfcn .
я =2
В решении разностной схемы Фурье - компонента решения с длиной волны 2я!кт за время At меняется на величину G, причем изменение фазы решения равно
ф = arg G=arctg^) = Atlms = -ikmcAt + At^-Xf kl"+lCln+l,
а величина затухания в решении схемы за время А* есть
|G| = ехр(Д/ Re s) = ехр (А^(-1)Ч^С2п).
В решении дифференциального уравнения Фурье - компонента е' "х решения с длиной волны 2к1кт за время At меняется на величину Ge. Для определения точного значения коэффициента перехода подставим в волновое уравнение его фундаментальное решение
и~еа ег тХ и найдем, что a~~ikmc. Тогда и = е' №^x с> и, следовательно, коэффициент перехода для точного решения имеет вид
причем изменение фазы решения равно фе = arg Gi = ~kmcAt,
a|Ge| = l.
Дисперсией разностной схемы (фазовой ошибкой), называют величину [2,16]
Аф = ф - фе . Ясно, что
Аф = М^1)Ч^+1С2„+1 . «=1
Диссипацией разностной схемы (относительная погрешность в амплитуде) называют
величину [2,16]
Z = \Ge\-W. В рассматриваемом случае
<х> X = 1-ехр(Д/ Re 5) = 1-ехр (Д/(-1).
я=1
Диссипативная схема неявно вводит в уравнение искусственную вязкость, которую часто называют неявной (схемной) искусственной вязкостью в отличие от явной искусственной вязкости, которая преднамеренно вводится в разностное уравнение. Искусственная вязкость сглаживает решение уравнения, уменьшая градиенты всех параметров независимо от причины возникновения этих градиентов, физической или вычислительной. Такое свойство разностной схемы обусловлено наличием в выражении для погрешности аппроксимации производных четного порядка. Дисперсия разностной схемы связана с производными нечетного порядка в выражении для погрешности аппроксимации. Дисперсия приводит к искажению соотношения фаз различных волн. Совместное воздействие диссипации и дисперсии на решение иногда называют диффузией. Диффузия приводит к растяжению крутых линий раздела, которые могут появляться в расчетной области.
а) Ь) с)
Рис.6. Влияние диссипации и дисперсна, (а) Точное решение, (Ь) Численное решение, полученное в том случае, когда ошибка является в основном дисснпативной (такое решение типично для схем первого порядка точности), (с) Численное решение, полученное в случае, когда ошибка является в основном дисперсионной (такое решение типично для схем второго порядка точности).
На рис. 6, показаны эффекты диссипации и дисперсии на расчет разрыва. Обычно если главный член в выражении для погрешности аппроксимации содержит производную четного порядка, то схема обладает в основном диссипативными свойствами, а если производную нечетного порядка - то дисперсионными.
Г) Консервативность.
Уравнение в частных производных записано в "дивергентной форме", или, что эквивалентно, в "консервативной форме", если коэффициенты при производных являются либо константами, либо функциями, производные которых в уравнение не входят. Обычно уравнения в частных производных, описывающие законы сохранения, записываются в дивергентной форме тогда, когда в них явно входит дивергенция той величины, для которой
этот закон формулируется. Например, дивергентная форма уравнения неразрывности (описывается закон сохранения массы) имеет вид
др дри dpv dpw
at сх ду dz
Последнее соотношение можно записать в векторной форме
-^- + divpV =0.
Уравнения в частных производных описывают физические законы сохранения в точке. Конечно-разностная схема обеспечивает близкую аппроксимацию уравнений в частных производных в небольшой области, содержащей несколько узлов разностной сетки. Те же законы сохранения, из которых выводятся уравнения в частных производных, справедливы для любой конечной области (контрольного объема). Если конечно-разностная схема даст близкую аппроксимацию уравнения в частных производных в окрестности каждого узла разностной сетки, то можно ожидать, что законы сохранения будут приближенно выполняться и для большего контрольного объема, содержащего довольно большое число узлов разностной сетки. Консервативной схемой называется разностная схема, обеспечивающая точное выполнение законов сохранения (исключая погрешности округления) на любой сетке в конечной области, содержащей произвольное число узлов разностной сетки [2]. Законы сохранения для всей сеточной области ("интегральные законы сохранения") для консервативных схем должны быть алгебраическим следствием разностных уравнений [5].
На важность требования консервативности схемы обратили внимание в начале 50-х годов А.Н. Тихонов и А.А. Самарский [3]. Ими был предложен интегро - интерполяционный метод для конструирования консервативных разностных схем. Сущность этого метода состоит в том, что разностные уравнения строятся на основе интегральных соотношений, выражающих законы сохранения для элементарной ячейки сетки. При этом на сетке вводится определенная интерполяция искомого решения и коэффициентов уравнения, изменяя которую можно получать различные разностные схемы.
Главным в определении консервативности разностной схемы является слово "точное" [2]. Любая согласованная разностная схема обеспечивает приближенное выполнение законов сохранения в большой области, но лишь консервативная разностная схема обеспечивает точное выполнение этих законов, вследствие взаимного уничтожения ряда членов уравнения. Проиллюстрируем это на примере решения уравнения неразрывности для установившегося течения: V pV = 0, Предположим, что мы построили конечно-разностную схему, аппроксимирующую это уравнение, и решили разностные уравнения во всей области
течения. Из закона сохранения массы, примененного к любому контрольному объему, который может совпасть со всей областью течения или составлять какую-то ее часть, следует, что суммарный расход газа через границу этого объема равен нулю (сколько вещества втекает в объем, столько и вытекает из него). Формально это можно показать, проинтегрировав записанное в дивергентной форме уравнение по всему объему и воспользовавшись формулой Гаусса-Остроградского
fflv-pV dR = \\pV nds=0 .
R S
Чтобы показать, что для решения уравнения неразрывности используется консервативная разностная схема, следует установить, что решение разностных уравнений удовлетворяет конечно-разностному аналогу этого интегрального тождества. Обычно это проверяется для контрольного объема, совпадающего со всей областью течения. Для этого вычислим интеграл в левой части, просуммировав конечно-разностные аналоги уравнения в частных производных во всех узлах разностной сетки. Если конечно-разностная схема консервативна, то при суммировании сократятся все члены, кроме тех, которые описывают поток массы через границу. Последние следует перегруппировать так, чтобы эти члены совпали с конечно-разностным аналогом интеграла в правой части закона сохранения массы. Для этого примера результат будет проверкой равенства втекающего в объем и вытекающего из него количества вещества. Если используемая конечно-разностная схема не консервативна, то внутри области могут появиться источники и стоки небольшой интенсивности.
Многие расчеты подтверждают тот факт, что применение уравнений в консервативной форме дает более точные результаты при расчете течения со скачками [14]. Это связано с тем, что ошибка аппроксимации конечно - разностных уравнений зависит от величины отброшенных высших производных при разложении в ряды Тейлора.
Существуют схемы, для которых выполнены не только разностные аналоги основных законов сохранения (массы, импульса и полной энергии), как для классических консервативных схем, но также ряд дополнительных сеточных соотношений, необходимость выполнения которых диктуется физическими соображениями. Схемы такого типа были названы полностью консервативными [3]. Полностью консервативные разностные схемы представляют собой сужение класса обычных консервативных схем. Для получения полностью консервативных схем можно использовать те же приемы, что и при построении консервативных схем. При этом дополнительно следует соблюдать некоторое формальное правило отбора: разностная схема должна одновременно аппроксимировать различные виды записи исходной дифференциальной системы уравнений, имеющие непосредственный физический смысл. Другими словами, отдельные разностные уравнения схемы должны быть
сформулированы таким образом, чтобы они допускали преобразования, аналогичные дифференциальному случаю.
Впервые полностью консервативная разностная схема для уравнений газовой динамики в эйлеровых координатах предложена Кузьминым А.В., Макаровым В.Л., Меладзе Г.В. Эта схема является трехслойной. Ими исследовано достаточно общее семейство двуслойных схем, содержащее 11 параметров, и показано, что среди них нет полностью консервативных.
Д) Монотонность.
При использовании разностных схем важно знать заранее характер поведения численного решения вблизи разрыва, в частности наличие или отсутствие осциляций. При построении разностных схем для нелинейных уравнений обычно требуют монотонность [5] этих схем для линейных уравнений. Но расчеты показывают [15], что монотонность в нелинейном случае часто нарушается по причинам, необъяснимым в рамках линейной теории. Описанный в [15] подход к оценке монотонности использует дифференциальное представление разностной схемы.
2.2.4. Обзор методов решения задач гидромеханики.
Часто начало современного численного анализа (вычислительной математики) связывают с появлением знаменитой работы Куранта, Фридрихса и Леви, опубликованной в 1928 г. [2]. В ней обсуждались характеристические свойства уравнений и в общих чертах излагался известный метод характеристик. В этой работе было так же получено и объяснено знаменитое необходимое условие устойчивости Куранта - Фридрихса - Леви, гласящее, что при расчетной сетке, не совпадающей с характеристикой, область зависимости разностных уравнений должна, по крайней мере, включать в себя область зависимости дифференциальных уравнений. Это условие устойчивости КФЛ справедливо для уравнений гидромеханики как в лагранжевых, так и в эйлеровых переменных [3,14]. (Существо подходов Лагранжа и Эйлера хорошо отражает следующее сравнение. Изучение движения воды в реке можно вести, либо плывя на лодке от истоков реки до ее устья и наблюдая за судьбой отдельных частиц жидкости (подход Лагранжа), либо наблюдая за течением с берега в определенных местах (подход Эйлера). Оба эти подхода эквивалентны друг другу в том смысле, что тот и другой дают полную картину движения среды [3].)
В течение 2-й мировой войны и сразу после ее охончания многие исследования были посвящены применению численных методов к решению задач гидромеханики. Именно в эти года проф. фон Нейман создал свой метод анализа устойчивости разностных схем решения
нестационарных задач. Примерно в то же время была опубликована статья Лакса, который одним из первых разработал метод расчета газодинамических течений с ударными волнами, которые являются поверхностями разрыва газодинамических параметров. При этом для расчета ударных волн не требовалось задания каких-то дополнительных условий. Лаксом и Вендроффом (1960 г.) было показано, что разрывные решения могут быть рассчитаны без специального рассмотрения разрывов, если уравнение рассматривается в консервативной форме и если схема консервативна [2]. Опыт показывает [14], что консервативные схемы, вообще говоря, дают более точные результаты. Чен и Аллен показали, что с помощью консервативной схемы получаются существенно более точные результаты для некоторых решений уравнения Бюргерса. Кроме того, дивергентная форма уравнений более осмысленна физически и облегчает постановку граничных условий для течений сжимаемой жидкости.
Утверждается [14], что для двумерных задач эйлеровы методы предпочтительнее, однако при их использовании затрудняется расчет ударных волн. Если размер ячейки сетки не меньше, чем толщина ударной волны, то появляются осциляции, снижающие точность. Эти осциляции на дискретной сетке имеют физический смысл. Кинетическая энергия, высвобождающаяся из-за потери скорости при переходе через ударную волну, превращается во внутреннюю энергию случайных соударений молекул; при расчетах роль молекул играют ячейки конечно-разностной сетки. В эйлеровых переменных невозможно брать особенно мелкую сетку вблизи скачка, так как его расположение заранее нешвестно. В лагранжевых переменных скачок может быть рассчитан на мелкой сетке. Этот подход плодотворен в случае одномерных задач, но трудноосуществим для многомерных задач. В общем случае наиболее эффективным методом расчета скачков является их искусственное размазывание. При этом утрачиваются детали течения внутри скачка, но выполнены законы сохранения при переходе через скачок.
Наиболее обычным подходом к расчету ударных волн на эйлеровой сетке является "размазывание" скачка на несколько ячеек сетки путем явного или неявного введения искусственной вязкости, не оказывающей влияния на решение на некотором расстоянии от ударных волн. В 1950 г. фон Нейман и Рихтмайер предложили схему искусственной вязкости, в которой "коэффициент вязкости" был пропорционален квадрату градиента скорости. Для размазывания скачка вместо явного введения искусственной вязкости можно использовать и неявную вязкость, имеющую место в конечно-разностных аппроксимациях. Это было осуществлено в методе Лакса и других методах. В 1960 г. Лаке и Вендрофф предложили новый метод расчета газодинамических течений с ударными волнами, позволяющий строить разностные схемы второго порядка точности, существенно меньше размазывающие ударные волны, чем использовавшиеся ранее конечно-разностные методы.
Предложенный Мак-Кормаком вариант этого метода является и сейчас одним из наиболее популярных методов расчета течений с ударными волнами [2].
В противоположность подходу с размазыванием скачка на несколько расчетных ячеек можно, наоборот, выделять разрыв. Одной из ранних была работа Гари, в которой этот подход использован для расчета движущихся скачков. В работах Моретти схемы с выделением скачков применены для расчета многомерного сверхзвукового обтекания с ударными волнами тел различной формы. Одна из первых схем с выделением скачка описана Рихтмайером и Мортоном. И сегодня для расчета течений с ударными волнами используются разностные схемы как с выделением, так и с размазыванием скачков.
Для решения уравнения вида
ду дх,
F = uA /2
в 1974 г. Лера и Пейре построили общий вид схемы из класса схем со следующими свойствами [16]:
явная консервативная схема;
двухшаговая, типа предиктор-корректор;
я+1
3) трехточечная двухуровневая (решение щ зависит от трех значений и на уровне
я);
4) второй порядок точности по времени и пространству. Схемы, обладающие всеми этими свойствами, образуют общий класс схем, зависящих от двух параметров а и р. Эти параметры характеризуют положение точки, в которой вычисляется предиктор щ (рис.7).
п+1
0Дх ,
К Ц
0Лк ,
X
п+а
-#
\^Л
t-1
i+1
Рис.7. Положение узлов вычисления предиктора «;
Эти схемы, обозначаемые S %, имеют следующий вид:
' 2аАх
[(а - fi^lb + (20 - \)Ftn + (\-а- ^Д + /} - FM ],
где Ft = Fyut j, a * 0.
Все конечно-разностные схемы, обладающие четырьмя вышеуказанными свойствами, относятся к классу S д -схем. Если же а = /3 = 112, то предиктор вычисляется в центральной
точке сетки: это есть двухшаговый вариант схемы Лакса-Вендроффа, предложенный Рихтмайером (1962). Если а~\, /? = 1/2, получается схема, предложенная Рубиным и
Берстейном (1967). Более общие схемы Syj были рассмотрены Мак-Гиром и Моррисом
(1973).
Схемы, соответствующие а = 1 и /? = 0 или fi = 1, являются двумя вариантами схемы
Мак-Кормака (1969). Если /7 = 0, предиктор щ определяется в точках xt = iAx и
вычисляется с помощью схемы разностей вперед, а щ - с помощью разностей назад. Если /? = 1, щ определяется в точках х,- = (/+1)Дх и вычисляется с помощью схемы разностей
назад, а щ - с помощью разностей вперед.
Класс схем Sq и S" был рассмотрен Уормингом, Катлером и Ломэксом (1973).
Численные эксперименты Лера и Пейре показали, что различные So -схемы дают
неидентичные результаты, в частности, в случае ударных волн. Некоторые из этих схем дают профили ударной волны с сильными осциляциями, другие - более регулярные профили, в некоторых случаях - монотонные. Такие различия в поведении численных решений появляются только в нелинейных случаях, но в линейном случае все схемы идентичны.
Почти все расчеты практически важных многомерных задач о течениях сжимаемой жидкости, опубликованные к 1980 году [14], проводились при помощи явных схем. Построению неявных схем для гиперболических систем уравнений посвящены ранние работы Вендроффа, Анучиной и Гари. В работе Браиловской с соавторами отмечена причина пессимизма относительно использования неявных схем: многие неявные схемы, безусловно устойчивые в применении к модельному уравнению, не являются таковыми в применении к системе уравнений, описывающих течение сжимаемой жидкости. Это подтверждается
опытом Шредера и Томсена; они построили неявную схему для многомерных гиперболических уравнений, которая по-прежнему требовала ограничения (16).
Полежаев рассчитал течения сжимаемого газа при помощи метода чередующихся направлений, при котором устранялось ограничение на шаг по времени, связанное с диффузией, но явная аппроксимация членов с градиентом давления приводила к сохранению условия (16).
Гурли и Митчелл также рассматривали метод чередующихся направлений и для двумерных гиперболических уравнений разработали безусловно устойчивую схему метода чередующихся направлений, основанную на девятиточечном шаблоне, применяемом на обоих слоях по времени. Однако эта схема не была опробована на нелинейных задачах и на реальных газодинамических расчетах. Шварц и Вендрофф при помощи нелинейно неявной схемы, используя итерации на каждом шаге по времени, рассчитали одномерную задачу о распространении ударной волны, но эту схему трудно обобщить на случай двух или трех пространственных переменных [14].
3. Краткое содержание работы.
Получение двумерной модели
Для адекватного учета явления разбавления облака газа окружающим воздухом будем считать, что плотность р, концентрация с, вертикальная скорость w и энтальпия h имеют неоднородное распределение по высоте облака р p(x,y,z,t), c=c(x,y,z,t), w = w(x,y,zj), h= h(x,y,z,t); za(x,y) z H(x,y,t), где z0- рельеф подстилающей поверхности, H -верхняя граница облака.
Будем предполагать, что все частицы облака газа или жидкости по высоте при любых фиксированных (x,y,t) движутся с одинаковыми горизонтальными скоростями, за исключением лишь окрестности нижней границы облака, то есть u(x,y,z,t)&u(x,y,t). v(x,y,z,t)&v(x,y,t). Если на практике для конкретного газа или жидкости и конкретной подстилающей поверхности это условие не выполнено, то в любом таком течении можно выделить слои, где это предположение выполнено, а далее для каждого конкретного слоя провести приведенные ниже рассмотрения и построить общую модель течения ИЗ послойных моделей, "сшитых" с помощью граничных условий.
В окрестности нижней границы течения у поверхности земли скорости резко замедляются и на нижней границе выполнены условия и = vj = 0, их = uv = иг = vx = vy - vz = wx = wy = 0, при любых (xyy,i). Нас будут интересовать осредненные по высоте облака параметры течения j Н{х,у,і) p(x,yj)=y, л jp(x,y,zj)dz, (1.2) где S(x,y,t) = H(xty,t) - zo(x,y) - толщина облака газа или жидкости. Аналогично определяются c(x,yj), h(x,y,t).
По величине этих параметров можно судить, например, о достижении предельно допустимых концентраций ядовитых веществ в облаке или о достижении пределов воспламенения.
Проинтегрируем уравнения системы (1.1) по высоте от нижней границы z0(x,y) до верхней границы Н(x,yj) и получим двумерную (х,у) систему уравнений для определения изменения в горизонтальном направлении средних по высоте характеристик течения. При этом толщину $(х,у, t) будем также рассматривать как одну из искомых функций двумерной модели.
Рассмотрим третье уравнение исходной системы; dt ах су az cz где Зс\ дх) ду\ су) dz\ dz j
Так как газ (жидкость) - молекулярно слабо связанная система, можем расщепить уравнение для w на два уравнения, одно из которых опирывает процесс переноса вертикальной скорости, а второе - изменение давления за счет гравитации:
Последний член описывает трение о землю и препятствия и его можно аппроксимировать формулой [17]: a dz pv (ЇЛО) := „ = ( / № 1 1 = 2 + 2 + 2 (1.11) FuH = PV Сгруппируем все оставшиеся граничные члены следующим образом (М dudH дидН Л z=H dz дх дх dy dy г fu =\pu ач ант (1.12) vu— + V M а дх dy 2=Ц
Член Fuff описывает действие силы сопротивления воздуха вдоль оси х на верхней границе облака и его можно аппроксимировать формулой, аналогичной (1.10): FUII=-4PO№- (1.13) Рассмотрим член /м. Отметим, что при последующем интегрировании всех остальных уравнений системы, получатся подобные выражения / для соответствующих субстанций р с множителем p(pt(p = (utv,w,c,h), поэтому естественно предположить, что вертикальная скорость w на верхней границе облака и изменение профиля границы связаны соотношением FWJJ описывает действие вертикальной компоненты силы сопротивления воздуха (подъемную силу) в результате чего происходит процесс смешения (разбавления) облака окружающим воздухом и его можно аппроксимировать формулой, аналогичной (1.10), (1.13) Fwff=- fi(p\w. (1.18) WZn описывает изменение вертикальной скорости из-за теплообмена с подстилающей поверхностью и его можно аппроксимировать аналогичной формулой FWZu=-iwpW- (1.19)
С учетом сделанных выше предположений, после интегрирования уравнения (1.4) получаем; a dpwS dpuw S dpvwS ck ду FDw+FwH +Fwzu + и zeSn (1.20) где a( „ тдЛ ,_(- s] F =&\fV ck FwH = -%wPo\U\w , K z0 = Cwp№ Аналогично будут осредненьї и уравнения для концентрации с и энтальпии h. При этом граничные члены вида (1.12) /с = 0, //, = 0 в силу (1.14),(1.15), а оставшиеся члены на верхней и нижней границе FcH , F JJ , F , Ff,7 вида (1.17) будут аппроксимированы формулами вида (1.18), (1.19):
Явные схемы для невязкого уравнения Бюргерса
В настоящее время схемы первого порядка точности почти не используются для решения гиперболических уравнений в частных производных. Метод Лакса выбран как типичный метод первого порядка точности для того, чтобы показать, что такие методы позволяют решать нединейные уравнения, но обладают сильными диссипативными свойствами [2,4,14-16].
Для построения разностной схемы воспользуемся дивергентної формой записи уравнения Бюргерса - + - = 0, = /2. а ас
Применим метод Лакса. Для этого выпишем первые два члена ряда Тейлора для функции и в точке (x,t): u(xJ + At) = u(x,t) + Al\ При помощи исходного уравнения заменим производную по времени, тогда запишем u(x,t + At) = u(x,t)-At\—\ +- .
Следуя методу Лакса, для аппроксимации производной яспользуем центральные разности, а первое слагаемое в правой части представим как среднее арифметическое значение в двух соседних узлах. В результате получим "j 2 Ах 2 Для уравнения Бюргерса F = и / 2,
Это явная одношаговая схема первого порядка точности с погрешностью аппроксимации «(ДїДДх) / Af). Условие устойчивости схемы Лакса имеет вид А/_ Дх"тах SI, где мгаах - максимальное собственное значение матрицы А, в случае уравнения Бюргерса состоящей лишь из одного элемента и. Модифицированное уравнение имеет вид
Эта схема не всегда удовлетворяет условию согласованности, так как отношение (Ах) I At может не стремиться к нулю при At и Дх, стремящихся к нулю. Однако, если при стремлении At и Ах к нулю число Куранта сохраняется постоянным, то условие согласованности выполняется. Схема Лакса отличается высоким уровнем диссипации при umaxAt I Ах &\. Причем, как можно видеть из численных расчетов, диссипативные свойства метода проявляются тем сильнее, чем больше отличается u At/Ax от 1. Важным свойством метода Лакса является его монотонность, то есть отсутствие осциляций решения. Схема Лакса очень просто обобщается на случай двух или трех пространственных переменных. Несмотря на свои недостатки, Схема Лакса имеет важное преимущество - простоту. Именно из-за простоты программирования и надежности схема Лакса может применяться на ранних стадиях разработки алгоритма решения задачи с тем, чтобы впоследствии заменить ее более сложной схемой [14].
Метод Лакса-Вендроффа - один из первых конечно-разностных методов второго порядка точности, созданных для решения гиперболических уравнений в частных производных. Для нелинейных уравнений разностную схему можно построить, исходя из разложения в ряд Тейлора [2,14-16]: а) x.t 2 \32) к . ди) ШУ u(x,t + At) = u(x,t) + М—\ + +. XJt
Первую производную по времени можно заменить при помощи исходного уравнения в частных производных. Для замены второй производной запишем исходное уравнение в виде скі dF dt дх Дифференцируя его по времени, получаем &и F &F ±_{ F др- dtdx ckdt dt \ dt. где порядок дифференцирования функции F изменен. Так как F = F(u), то ди (X7 SF ди ди dF dF ди ди dt дх да дх дх St ди dt dt Следовательно, производную dF I St можно заменить по формуле
В случае уравнения Бюргерса матрица Лкоби А состоит лищь из одного элемента и. Если решается система уравнений, то и и F - векторы, а А - матрица. Подставляя найденные производные в разложение функции и в ряд Тейлора, получаем (cF Ш)2 д( 6F\
Для построения схемы Лакса-Вендроффа теперь достаточно вместо производных подставить их центрально-разностные аппроксимации
Матрица Якоби А вычисляется в середине между узлами разностной сетки. В случае уравнения Бюргерса А = и, то Aj+l/2 =(»/ +Uj+\)I% Aj-Vl = («у + г/y-i )/2. Это явная одношаговая схема второго порядка точности с погрешностью аппроксимации ЩЬх) ,(Д0 ). Условие устойчивости имеет вид Дх max 1. Модифицированное уравнение имеет вид щ + иих = (At) (" t x)xx V) Коэффициент перехода для рассматриваемого метода вычисляется по формуле (Ы Л2 М G=l-2—A\ (1 - cos /J)-2i A sin fi. \Ax J Ax На рис.2.3 показаны результаты расчета методом Лакса-Вендроффа той же модельной задачи, что и для метода Лакса. Осциляции решения подчеркивают преимущественно дисперсионные свойства разностной схемы.
Этапы решения двумерной системы
Вернемся к рассмотрению двумерной системы (1.26). Произведем расщепление системы (1.26) по физическим процессам аналогично параграфу 3,1, Опишем поэтапно (и+1)-й шаг численного решения этой системы, считая, что йп ,vn ,рп ,с" ,дп известны с предыдущего шага по времени.
Таким образом, вместо одной сложной двумерной системы (1.26), описывающей несколько различных физических процессов, рассматривается пять небольших подсистем, каждая из которых описывает конкретный физический процесс. Э. Разностный метод решения двумерной задачи.
Для решения каждого из этапов, описанных в предыдущем параграфе использовались явные разностные схемы на двумерной прямоугольной сетке с постоянными шагами по л и по у [2-6,Ю]: i=i-hx, i = \Nx, yj j-hy, j = \Ny. 1 этап. Аппроксимация с порядком Olr,/ . +hy\ в полуцелых точках явной разностной схемой: (-х\а ""+ J+V2J »"/+1/2,; g( ,„2\по / ,о\« ( \0 і (-я\я v iJ+№ -v І./+1/2 gt ,x2\"o J ІХ\я { \0 X; = Ml/ =u", V
Этап. Схема с расщеплением по времени. Решается уравнение переноса по явной схеме типа предиктор-корректор с расщеплением по времени. Эта схема проверена на модельных задачах, полученные результаты сравнивались с известными результатами.
Метод с расщеплением по времени таким образом "расщепляет" оригинальную схему, что решение многомерной задачи сводится к последовательному решению одномерных задач [2]. Благодаря этому условие устойчивости разностной схемы становится менее жестким (менее ограничительным). Другими словами, расщепление по времени позволяет получить решение в каждом направлении с максимально допустимым шагом по времени. Особенно заметно преимущество расщепления в том частном случае, когда максимально допустимые шаги по времени Atx и Aty сильно различаются друг от друга из-за различия шагов Дх и Лу разностной сетки.
Чтобы записать метод расщепления, воспользуемся одномерными разностными операторами Lx[Atx) и іДд/ І. Если оператор Lx[AtxJ применяется к величине U"j, то выражение Погрешность аппроксимации этой разностной схемы 0\т ,hx2 +hy21, где 2 2 2 т =тх +Ту.
В общем случае разностная схема, полученная при применении такой последовательности операторов, удовлетворяет следующим условиям: 1) устойчивости, если для каждого оператора шаг по времени не превосходит максимально допустимый для этого оператора; 2) согласованности, если суммы шагов по времени для каждого оператора совпадают; 3) аппроксимации со вторым порядком точности, если последовательность операторов симметрична.
Разности вперед (В) и назад (Н) можно последовательно чередовать как на шагах предиктор-корректор, так и при аппроксимации производных по двум пространственным
Производные в вязких членах Е и F дискретизируют следующим образом. Имеющиеся в Е производные по координате х аппроксимируются разностями противоположного направления относительно тех, которые используются при аппроксимации Ж / дх, тогда как производные по направлению у аппроксимируются центральными разностями. Аналогично, производные по координате у в F аппроксимируются разностями противоположного направления относительно тех, которые используются при аппроксимации 3F і ду,
Применим схему с расщеплением по времени к уравнению (3.15). Для оператора Lx получим (входными данными для второго этапа служат и и v, вычисленные на первом этапе ,и S и с, взятые с предыдущего шага по времени):
Расползание газа по неровной поверхности
Рассмотрим растекание жидкости по неровной поверхности. Изначально жидкость находится в резервуаре цилиндрической формы, стоящем на поверхности (zo -рельеф поверхности), а затем, после мгновенного раскрытия стенок резервуара, жидкость начинает растекаться по поверхности под действием гравитационных сил. Таким образом, имеем задачу с начальными условиями. Система, описывающая данную задачу, выглядит следующим образом: + .+. Начальные условия вне источника: «о/=0 = » vor=0 = 0, Ч=0 = () Начальные условия на источнике:
Источник импульсный, действует только в начальный момент времени. Площадь источника равна площади основания цилиндрического резервуара. Рассматривалось две формы основания цилиндра, четырехточечная и двенадцатиточечная: При одинаковых начальных данных для обоих источников из расчетов можно было видеть, что при растекании двенадцатиточечного источника жидкость имеет более округлую форму. Это связано с тем, что более подробно задано основание цилиндра, т, е. окружность. Значения const, с которыми производились вычисления; % = hy = 0.82 , р 1000, 6 = 0.05 (толщина фона), d\ 1 (высота цилиндра).
Из численных экспериментов было установлено, что большое влияние на толщину жидкости оказывает коэффициент турбулентной вязкости v. В самых первых расчетах коэффициент v задавался не формулой, a const (0.01 -ь 1). При этом после нескольких шагов вычислений на сетке появлялись осциляции, причем они были тем меньше, чем больше был коэффициент вязкости v. Однако увеличение v приводило к более медленному расползанию жидкости и, следовательно, к долгим временам счета. Поэтому возникла идея менять коэффициент v со временем (чем дальше от начала процесса, тем меньше v). На первых порах это проделывалось вручную, но затем было решено задавать v формулой, вытекающей из известной модели Колмогорова, что погасило осциляции: v = const \U\ д .
Также было выяснено, что увеличение фонового значения S гасило осциляции. Поэтому, пока коэффициент v задавался covst , увеличение фоновых значений 5 в совокупности с увеличением v давало более - менее гладкие результаты. Когда же v стали задавать формулой, уменьшение фоновых значений 8 не приводило к осциляциям, только в начале вычислений значительно уменьшался автоматический шаг по времени из-за большой разницы высот источника и фона.
Расчеты проводились для 4-х основных случаев: 1. Растекание жидкости по ровной поверхности ( ZQ = 0). 2. Растекание жидкости по наклонной поверхности (z0 =а х,а-коэффициент наклона плоскости). Было рассмотрено несколько вариантов в зависимости от наклона плоскости a = (aha2): а) щ = 0.001 + 0.1,а2 =0 (рис. 4.2; щ =0.1); б) щ =а2 = 0.001+0.1. 3. Обтекание жидкостью препятствия. В качестве препятствия бралась пирамидка с основанием "крест" (пять точек). Высота пирамидки варьировалась от У высоты S на источнике до высоты источника. 4. Натекание жидкости в углубление.
Динамику расползания жидкости по ровной поверхности можно видеть на рис.4.1. На рис.4.2 представлено сползание жидкости по наклонной поверхности (жидкость ползет в сторону наклона), а на рис.4.3 - огибание жидкостью препятствия. В качестве препятствия бралась пирамидка с основанием "крест" (пять точек). Как можно видеть из рисунка, жидкость за препятствием сомкнулась, что говорит в пользу используемой численной модели. Подобные расчеты проводились ИАЭ им. Курчатова [8], однако, при расчетах использовалась упрощенная численная модель и жидкость за препятствием не смыкалась. Также в качестве неровной поверхности рассматривалась поверхность с углублением. Постепенно расползаясь, жидкость заполняет углубление (рис.4.4). Как видно из представленных рисунков, предложенная двумерная модель воспроизводит качественные особенности явления.
