Дипломная (вкр): Моделирование миграции и осаждения примеси на основе конвекционно-диффузионной модели

(3.1)

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

С учетом изложенного рассмотрим задачу Коши

    (3.2)

где  в D при t=0

Рассмотрим простейшие методы аппроксимации задачи (3.2) по времени, полагая, что не зависит от времени. Одна из них - явная схема первого порядка аппроксимации на сетке

    (3.3)

где

Неявная схема первого порядка аппроксимации имеет вид

          (3.4)

где

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

Схема Кранка - Николсона имеет следующий вид [6]:

      (3.5)

где

Эта схема аппроксимирует исходную задачу со вторым порядком по времени.

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

Рассмотрим методы расщепления, применяемые при решении двумерных задач, описывающих процессы различной природы. Рассматриваемое эволюционное уравнение имеет вид (3.1) в . Оператор  не зависит от времени и представим в виде  при условии, что .

Будем полагать, что эта задача уже редуцирована к разностному виду.

3.1 Метод стабилизации

Рассмотрим так называемый метод стабилизации. Для этого рассмотрим разностную схему решения (3.1) в предположении f=0:

      (3.6)

где

Эта схема аппроксимирует исходную задачу со вторым порядком аппроксимации по .

С помощью ряда преобразований приводится к виду

          (3.7)

где .

Отсюда видно, что (3.7) при достаточной гладкости решения совпадает со схемой Кранка - Николсона, т.е. имеет второй порядок аппроксимации по .

Эта разностная схема допускает удобную компьютерную реализацию.

В случае неоднородной задачи

,    (3.8)

где  при t=0.

В этом случае задача запишется следующим образом:

      (3.9)

где ,.

При условии  схема будет обладать вторым порядком аппроксимации по .

3.2 Метод покомпонентного расщепления

В случае зависимости операторов  от времени применяется метод покомпонентного расщепления.

Пусть в (3.1) оператор  не зависит от времени и представим в виде  при условии, что  и . Рассмотрим аппроксимацию этих матриц на интервале  в форме

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

Построим систему разностных уравнений, состоящую из последовательности простейших схем Кранка - Николсона

           (3.10)

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

   (3.12)

Если , , то при достаточной гладкости элементов этих матриц и решения  задачи (3.1) разностная схема (3.10) абсолютно устойчива и аппроксимирует исходное уравнение со вторым порядком по  в случае, если  и  коммутируют, т.е. , и с первым, если не коммутируют.

Рассмотрим двуцикличный метод покомпонентного расщепления. Будем аппроксимировать операторы  и  на интервале . Положим

.

Построим следующие две системы разностных уравнений

         (3.13)

         (3.14)

Имеем

,    (3.15)

Двуцикличный метод абсолютно устойчив, а схема (3.15) аппроксимирует исходное уравнение (3.1) со вторым порядком по

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

Рассмотрим систему разностных уравнений вида (3.13), (3.14), записанных в более удобной форме

      (3.16)

где . Разрешая эти уравнений относительно , получим

       (3.17)

           (3.18)

(3.19)

С помощью разложения по степеням малого параметра  придем к соотношению

   (3.20)

которое преобразуем к виду

    (3.21)

Исключим , используя разложение решения в ряд Тейлора в окрестности точки . С точностью до  будем иметь

    (3.22)

Производную  исключим с помощью соотношения

(3.23)

Подставим (3.23) в (3.22). Тогда

          (3.24)

          (3.25)

Подставим соотношение (3.25) в (3.24). В результате будем иметь

          (3.26)

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

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

4. Описание программной реализации решения двумерной задачи

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

- параметры области решения задачи;

-           компоненты вектора скорости воздушных масс в направлении осей  соответственно ;

-           коэффициенты диффузии в направлении осей  соответственно ;

- мощность источника примеси;

- величина, характеризующая взаимодействие примесей с подстилающей поверхностью (л);

- координаты источников, в которых проводятся наблюдения за распространением примеси;

- количество итераций по времени исследования;

-           количество шагов по временной переменной;

- количество шагов по пространственным переменным.

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

4.1 Выбор среды реализации

Для написания программы была использована среда Compaq Visual Fortran 6.5. Этот выбор обусловлен тем, что Fortran - один из самых мощных языков программирования, позволяющий работать с типами данных повышенной точности, что очень важно при выполнении математических расчётов. Fortran изначально был создан для научных и численных расчетов и всё его последующее развитие ориентировано прежде всего на подобные приложения.

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

4.2 Описание программы

Все входные параметры находятся в отдельном текстовом файле in.txt, их можно изменять непосредственно в этом файле. Ввод данных осуществляется под управлением именованного списка.

Имя списка ввода в данной программе: input. Таким образом, оператор NAMELIST объявления именованного списка ввода есть в разделе объявлений программной единицы и имеет вид: namelist /input/ L, h, tau, T, n, u, w, mu, nu, M, constX, constZ, lambda где после /input/ идёт перечисление заранее определённых в программе в принадлежности какому-либо типу переменных. При вводе именованного списка оператор ввода ищет в файле in.txt начало списка, которое имеет вид: $input. Перечень принадлежащих именованному списку данных завершается знаком доллара ($). Имена переменных во входном файле при использовании именованного списка ввода должны совпадать с соответствующими именами списка переменных оператора NAMELIST.

Написанная программа реализует схему расщепления по физическим процессам.

Рассмотрим первый этап, который соответствует переносу. Разобьем весь промежуток [0, T] на элементарные. Переносу соответствует следующий оператор:

.

Соответствующий разностный аналог имеют вид

где  - шаги сетки вдоль соответствующих осей.

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

.

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

где

За  при этом принимается решение задачи на предыдущем этапе (когда рассматривался только перенос).

В результате получается схема второго порядка аппроксимации по  и по .

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

,

,

,