LINEBURG


<< Пред. стр.

страница 7
(всего 9)

ОГЛАВЛЕНИЕ

След. стр. >>

? ij = ? ? S ( xi ? x L , y i ? y M ) ,
S =1

где N – число точечных источников, а функция

70
? S ( xi ? x L , yi ? y M ) = Q при ( xi = x L , yi = y M ).
На рис. 5 показаны изолинии концентрации аэрозоля для точечного
источника, полученные в аналитической модели для параметров среды
v = 0.1 м/с, u = 0 , ? = 0 , ? = 20 м2/с, интенсивность источника Q = 1 .
На рис. 6 показаны изолинии концентрации для точечного источника
при Q = 1 , ? = 1 м2/с, ? = 0 , v = 10 м/с. Как видно из рис.6 , наличие силь-
ного ветра приводит к сносу выбросов по ветру довольно узкой полосой.
Из приведенных расчетов видно, что аналитическая модель хорошо




Рис. 5. Распределение концентрации аэрозоля, выбрасываемого точечным ис-
точником. Случай близкий к штилевой ситуации




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




Рис. 6. Распределение концентрации аэрозоля в выбрасываемого точечным ис-
точником при скорости ветра 10 м/сек




72
Список литературы
1. Марчук Г.И. Математическое моделирование в проблеме окружающей
среды. М.: Наука. Гл.ред. физ.-мат. лит., 1982. 320 с.
2. Янке Е., Эмдэ Ф., Леш Ф. Специальные функции. М.: Наука, 1968.
344 с.



3.3. Двумерная численная модель


3.3.1. Формулировка стационарной задачи
( x, y )
В плоскости с декартовыми координатами задана
прямоугольная область:
0 ? x ? x0
(1)
0 ? y ? y0

На границе области задана концентрация примеси ? :

? = ?1 ( y )
x =0

? = ? 2 ( y)
x = x0
(2)
? = ? 3 ( x)
y =0

? = ? 4 ( x)
y = y0


В случае, когда граница области расположена достаточно далеко от
источников, можно считать, что примесь на ней отсутствует, то есть
?1 = ? 2 = ? 3 = ? 4 = 0 . Именно последний случай, как правило, мы и будем
рассматривать.
Внутри прямоугольника требуется решить дифференциальное уравне-
ние

? ? ?? ? ? ? ?? ?
? ?
(u? ) + (v? ) + ? k ? ?y ? + ?? = Q( x, y ) ,
? + ?k (3)
?
?x ?y ?x ? ?x ? ?y ? ?



73
где ? – концентрация примеси, (u, v) – компоненты вектора скорости в на-
правлении по x и по y; k ( x, y ) – коэффициент диффузии; ? ( x, y) – коэффи-
циент поглощения; Q(x,y) – функция источника. Функции k ( x, y ) ,
? ( x, y) > 0 , Q(x,y) и числа (u, v) должны быть заданы, ? ( x, y) – искомая
функция.

3.3.2. Общая схема численного решения задач
Одним из распространенных методов решения стационарных задач
является метод установления, в его рамках формулируется нестационарная
задача, решения которой с течением времени приближаются к искомому
решению стационарной задачи. При этом процесс установления численно-
го решения не обязательно соответствует физическому процессу установ-
ления. Более того, стационарное решение достигается быстрее, если этого
соответствия не требовать. Здесь мы не претендуем на описание нестацио-
нарных явлений, все зависимости от времени понимаются как формальный
прием.
При численном решении сложных задач, описывающих совокупность
нескольких физических процессов, эффективным приемом является так
называемое расщепление [1].
Процессы диффузии и поглощения удобно рассматривать совместно.
Обозначим:
? ? ?? ? ? ? ?? ?
A? = ? ? ?y ? + ?? ,
? ? ?k
?k ?
?x ? ?x ? ?y ? ? (4)
?
(u? ) + ? (v? ) .
B? =
?x ?y
Тогда нестационарный аналог (3) запишется в виде:
??
+ A? + B? = Q . (5)
?t



74
Идея метода расщепления по физическим процессам [1] состоит в
том, что шаг по времени ? разбивается на две половины, на каждой из ко-
торых один процесс интегрируется по явной схеме, а другой – по неявной.
Обозначим индексами 0, 1, 2 значения ? в моменты времени t, t + ? / 2 ,
t + ? . Тогда при заданном ? 0 требуется найти ? 2 из уравнений:
?1 ? ? 0
+ A? 0 + B? 1 ? Q = 0 ,
? /2 (6)
?2 ? ?0
+ A? 2 + B? 1 ? Q = 0 .
? /2
Промежуточное значение ?1 определяем независимо из первого урав-
нения, после чего подставляем во второе, из которого определяем ? 2 . И
этот шаг по времени повторяется достаточно долго, пока решение не
установится. В этом случае ? 0 = ?1 = ? 2 , и из (6) получаем:
A? 0 + B? 0 ? Q = 0 ,
то есть на стационарное решение расщепление, во всяком случае, не
влияет.

3.3.3. Аппроксимация
Выразим ? 0 и ? 2 через ?1 из уравнений (6):

? ? ? ?? ? ? ??
? 0 = ?1 ? A ? ??1 + B ?? 1 ? Q ? ,
? 2 ? ?? 2 ? 2?
?1
? ?? ? ? ??
? ?
? 2 = ?1 + 1 ? B ?? 1 ? Q ?
A? ?? 2
? 2? ?? ? 2?

при малых ? , когда ? << 1 A , ? << 1 B , приближенно получаем

? ? ?
? ?
? 0 = ?1 + B ?? 1 ? Q + O (? 2 ) ,
A+
? 2?
2 2

?? ?? ?
? 2 = ?1 ? A ? B ?? 1 + Q + O (? 2 ) .
? 2?
2 2
Поэтому

75
?0 + ?2
= ? 1 + O(? 2 ) .
2
Подставим это выражение в полусумму уравнений (6):
?2 ? ?0 ? + ?0 ? + ?0
? Q = O (? 2 ) ,
+A 2 +B 2
? 2 2
то есть построенный итерационный процесс при ? > 0 действительно
аппроксимирует нестационарное уравнение.



3.3.4. Организация итераций
В итерационном процессе ?1 находим из первого уравнения (6):

?2 ?
? + B ? ? 1 = P1 , (7)
?? ?
?2 ?
где P1 = ? ? ? + A ? ? 0 + Q , ? 0 – уже известно из предыдущего шага.
?? ?
Решив (3), найдем ? 2 из уравнения

?2 ?
A ? ? 2 = P2 ,
?+ (8)
?? ?


?2 ?
где P2 = ? ? ? + A ? ? 1 + Q .
?? ?



3.3.5. Выбор итерационного параметра
При постоянных коэффициентах и равномерной сетке лучшая сходи-
мость к стационарному решению имеет место при
4 ?2
? ? h x0 , (9)
k 1 + ? x (k? )
2 2
0


где h – шаг сетки, x0 – размер области.
В частности, при больших значениях коэффициента поглощения ?

76
4? 4h x
, в обратном случае, когда ? = 0 – ? 0 ? 2 0 .
? ?h
?k
k?
На практике установлено, что выгоднее (численное решение быстрее
сходится к точному) чередовать ? через шаг по времени. В большинстве
моделируемых случаев поглощение малозначительно, поэтому для опреде-
ления ? используется приближение ? = 0 . Значения ? для нечетного и чет-
ного шагов по времени шагов
?1 = ? 0 N,
(10)
?2 =?0 N ,
где N – число делений разностной сетки.



3.3.6. Дискретная модель для диффузии и поглощения
Горизонтальные и вертикальные линии сетки нумеруем целыми ин-
дексами n и m соответственно. Соответственно узлы сетки нумеруем парой
чисел (n,m), горизонтальные отрезки – (n+1/2,m), вертикальные отрезки –
(n,m+1/2).
Определим также два типа ячеек сетки. Ячейка (n+1/2,m+1/2) огра-
ничена ребрами сетки, а ячейка с целыми индексами (n,m), заштрихован-
ная на рис. 9, имеет в качестве узлов центры полуцелых ячеек.
Построим оператор An, аппроксимирующий оператор A в уравнении
(4), в проинтегрированном по ячейке (n,m) виде. Считаем функцию ? оп-

ределенной в узлах сетки. Соответствующие значения обозначаем ? n . По-
m



скольку сетка прямоугольная, расстояние между узлами (n+1,m) и (n,m), то
есть шаг сетки по x, не зависит от y, и шаг сетки по y не зависит от x. По-
этому шаги можно пронумеровать одним полуцелым индексом hn +1 / 2 и

h m +1 / 2 .




77
Рис.7. Элементы конечно-разностной сетки


Производные ? определим на ребрах сетки:
m
? n +1 ? ? n
? ?? ?
m m
=
? ?
?x ? n +1 / 2 hn +1 / 2
?
(11)
m +1 / 2
? n +1 ? ? n
? ?? ?
m m
=
? ?
h m +1 / 2
?x ? n
?
Коэффициент диффузии k изначально определен постоянным в каж-
дой полуцелой ячейке (n+1/2,m+1/2). Сопоставим каждому ребру сетки
среднее по прилегающим ячейкам значение:
h m +1 / 2 kn ++1/ /22 + h m ?1 / 2 kn +?1/ /22
m m
=
m
k 1 1
(12)
n +1 / 2 m +1 / 2 m ?1 / 2
+h
h



78
Тогда поток через правую границу ячейки (n,m), который мы сопос-
тавляем ребру (n+1/2,m), можно вычислить как
1 ? n +1 ? ? n h m +1 / 2 + h m ?1 / 2
m m
= ?k + ?
m m
J (13)
n +1 / 2 n
2 hn +1 / 2 2
Площадь ячейки (n+1/2,m)

( )(
1 m +1 / 2
+ h m ?1 / 2 hn +1 / 2 + hn ?1 / 2 )
Sn =
m
h (14)
4
Для оператора поглощения простейшей аппроксимацией будет
Sn ? n ? n .
mmm
(15)
Объединяя формулы (13,15), получаем

m ?1 / 2 m ?1 / 2 ? n +1 ? ? n
m m
( )
1 m +1 / 2 m +1 / 2
An? n = ? h k n +1 / 2 + h k n +1 / 2
hn +1 / 2
2
? n ?1 ? ? n
m m
( )
1 m +1 / 2 m +1 / 2
k n ?1 / 2 + h m ?1 / 2 k n ??1/ /22
?h m
1
2 hn ?1 / 2
? n +1 ? ? n
m m
( )
1
? hn +1 / 2 k n ++1/ /22 + hn ?1 / 2 k n ?+1/ /22
m m
1 1
h m +1 / 2
2 (16)
? n ?1 ? ? n
m m
( )
1
? hn +1 / 2 k n +?1/ /22 + hn ?1 / 2 k n ??1/ /22
m m
1 1
h m ?1 / 2
2
+ Sn ? n ? n
mmm



Индексом h подчеркиваем, что ? n – это набор чисел ? n , используе-
m



мый для аппроксимации значений ? ( x, y ) в узлах сетки.
При постоянных коэффициентах и шагах сетки получается классиче-
ская пятиточечная разностная аппроксимация оператора Лапласа, допол-
ненная младшим членом при ? ? 0 .



3.3.7. Способ решения дискретных уравнений диффузии
Основные вычислительные трудности при решении уравнения (6) свя-
заны с обращением оператора диффузии с поглощением А (16). Поэтому


79
решать уравнение (6) будем наиболее эффективным из известных для эл-
липтических уравнений многосеточным методом Федоренко [2].
Будем использовать последовательность сеток, шаг которых увеличи-
вается в 2,4,8... раз. На каждой конкретной сетке уравнение решаем мето-
дом Зейделя с параметром релаксации ? an . Оптимальное значение ? an для
уравнения Лапласа:
? an = 1.25
Новое значение решения на k + 1 итерации в i-ом узле ? ik +1

выражается через ? ik следующим образом:

? ik +1 = ? ik ? ? an H ik L?1 ,

<< Пред. стр.

страница 7
(всего 9)

ОГЛАВЛЕНИЕ

След. стр. >>

Copyright © Design by: Sunlight webdesign