Применение монотонно разностной схемы для численного моделирования волновых процессов

Введение

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

Решение задачи Коши для систем уравнений гиперболического типа с использованием схемы Лакса-Вендроффа или схемы Мак-Кормака , имеющих второй порядок точности, дает осцилляции вблизи разрывов (см.[2,3]).

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

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

 

 

 

 

 

 

 

 

 

 

 

 

 

Рассмотрим линейное дифференциальное уравнение переноса:

        (1)                  ,

   

1. вектор скорости переноса.

Если в уравнении (1) правая часть f=0, то общее решение этого уравнения имеет вид бегущей волны:

        (2)

Для определенности положим  , тогда волна бежит слева направо.

В полуплоскости  введем расчетную сетку с шагами h по оси x и - по оси t. Ниже изображен расчетный шаблон. Здесь точки 1,2,3 находится на старом временном слое, а точка 4-на новом слое по времени.

 

 

 

 

 

рисунок 1

Перепишем уравнение (1) в эквивалентном виде

           (3)

Здесь b – постоянное число. Выбором соответствующего значения b преобразуем левую часть уравнения (3) в полную производную по t вдоль узловой линии.

Введем в рассмотрение функцию

                     (4)

Теперь уравнение (3) запишется в виде

                (5)    

Проинтегрируем уравнение (5) по t вдоль линий 12,24,34. При этом число b возьмем

          

 

соответственно   ,    ,   

Получим три расзностных  уравнения

          (6)

           (7)

     (8)

В формулах (6)-(8) введено число Куранта

Пусть разностную схему I2 составляет уравнение (4), (5) , а разностную схему I3- уравнения (4),(6). Исследуем схемы I2, I3 на устойчивость с помощью спектрального признака. Для этого перепишем их в стандартном виде

 

где

- известные числовые матрицы размерности 22. Для схемы I2

;              ;

  если i-1, i0.

ля схемы I3

;      

 если i-1, i.

 функции оператора  перехода при решении задачи  Коши в бесконечной области  имеют вид

 

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

  

Решая квадратное уравнение  , получим что при k1. Этот вывод справедлив для обеих схем. Таким образом, схемы I2 и I3устойчивы при числе Куранта, меньшем либо равным единице. Это утверждение следует и из необходимого признака устойчивости КФЛ. При измельчении шагов сетки разностное решение сходится к точному со вторым порядком точности в классе решений с ограниченной третьей производной.

Все расчеты задачи с начальными данными проводились до Т=1.0 для двух вариантов шагов сетки: һ=0.1,  и h=0.05,   Параметр k варьировался в интервале

В первом примере начальная  функция имела вид

Она пренадлежит к классу непрерывных функций и имеет  разрывы первой производной в  точках и .

Ниже представлены результаты расчетов по схемам I2, I3 для различных шагов сетки.

 

 

Рассмотрим задачу с двумя  пространственными переменными  в области 

 

        

Скорости переноса по осям a_x, a_y считаем положительными , для простоты, постоянными.

Введем по переменной x сетку а по переменной y сетку . Значения решения в узлах этой сетки обозначим следующим образом:

      

Возьмем шаблон, изображенный жирными линиями и составим на нем схему

 

 

рисунок 2

 

где – шаги по соответствующим направлениям.

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

 

 

Вычисления проводятся послойно. Значение в узле, отмеченном на рисунки 2, выражается по формуле