Модификация схемы «Донор-акцептор» для расчета диффузии завихренности и ее применение в методе «Вихрь в ячейке» Текст научной статьи по специальности «Физика»

УДК 532.527

МОДИФИКАЦИЯ СХЕМЫ «ДОНОР-АКЦЕПТОР» ДЛЯ РАСЧЕТА ДИФФУЗИИ ЗАВИХРЕННОСТИ И ЕЕ ПРИМЕНЕНИЕ В МЕТОДЕ «ВИХРЬ В ЯЧЕЙКЕ»

© 2003 В. В. Никонов, В. Г. Шахов Самарский государственный аэрокосмический университет

Предлагается модернизация схемы расщепления уравнений Навье-Стокса для численного метода "вихрь в ячейке". Модернизация заключается в введении в задачу малого параметра, пропорционального величине вязкости потока, что позволяет находить решение задачи в виде асимптотического ряда. Для тестирования предлагаемой схемы расщепления рассматривается двумерная задача о движении вихря Ламба-Озеена, имеющая аналитическое решение. Сравнение результатов, полученных методом «донор-акцептор» (Д-А) и модернизированным методом (МД-А), с аналитическим решением показало, что модернизированный метод является более предпочтительным при моделировании течений с малой вязкостью. Метод МД-А в отличие от Д-А сохраняет точность в широком диапазоне изменения шага интегрирования по времени.

1. Основные соотношения и алгоритм метода

В рамках метода «вихрь в ячейке» [1] уравнение Навье - Стокса для вязкого двумерного несжимаемого потока в переменных завихренность - скорость имеет вид

¿со 2 — = vV со dt

О)

где со - завихренность, I - время, V - кинематическая вязкость, V2 - оператор Лапласа. Здесь и далее все уравнения записаны в безразмерной форме. В численном методе «вихрь в ячейке» уравнение (1) расщепляется на подзадачи конвекции и диффузии завихренности. Ошибка такой схемы [1] имеет порядок шага по времени Ы. На этапе конвекции используется лагранжев подход к описанию движения жидкости:

г/со ~dt

dx ~dt

= 0

(2)

(3)

5(0 Т72

— -- vV со

dt

(4)

Вектор скорости и , стоящий в правой части (3) обыкновенного дифференциального уравнения (ОДУ), обычно представляется

в виде суммы потенциальной ир и соленои-

дальной (вихревой) и5 его составляющих. Потенциальная составляющая в данном случае равна скорости набегающего потока

1Л р — Ы оо . Соленоидальная составляющая поля скорости определяется с помощью векторного потенциала V)/ ег, где \\) - функция тока, ег - единичный вектор вдоль оси г:

Us = Vxfiy ez) ■

(5)

Предварительно функция тока находится из решения уравнения Пуассона

V2v|/ = -со

(6)

где х - радиус-вектор центра «вихря в ячейке», и - вектор скорости рассматриваемого вихря. На диффузионном этапе решается уравнение в частных производных (ДУЧП), аналогичное уравнению теплопроводности:

с помощью быстрого преобразования Фурье (БПФ). Данный подход, как показано в [2, 3], позволяет существенно экономить машинное время по сравнению с непосредственным определением скорости по закону Био-Сава-ра от всех вихревых частиц. Уравнение (6) решается с помощью подпрограммы библиотеки IMSL языка программирования Fortran 90.

Алгоритм численного метода «вихрь в ячейке» на примере одного расчетного шага имеет следующий вид:

- поле завихренности берется из начальных условий задачи или с предыдущего расчетного шага;

- расчет диффузии завихренности проводится методом «донор-акцептор» Д-А или модифицированным методом Д-А (МД-А);

- выполняется подготовка граничных условий для решения уравнения (6);

- определяется функция тока с помощью БПФ из уравнения Пуассона (6);

- рассчитывается поле скорости с помощью конечно-разностных формул;

- решается система ОДУ (3) методом Эйлера 1 -го порядка для каждого «вихря в ячейке»;

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

- осуществляется переход к следующему расчетному шагу.

2. Подзадача конвекции

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

В противном случае вклад в функцию тока вычисляется от каждого вихря рассматриваемого кластера

2л

(7)

члены, начиная с третьего порядка малости, и суммируя по всем вихрям данного кластера, можно получить выражение для функции тока для всего кластера. Данной точности аппроксимации вполне достаточно [3,4].

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

После расчета поля скоростей течения новые координаты «вихрей в ячейках» получаются численным интегрированием системы ОДУ (3) методом Эйлера. Новое местоположение вихрей не обязательно совпадет с координатами расчетной сетки, поэтому используется процедура перераспределения их интенсивности в ячейки сетки

(и )

Т'(хк,у,)К(х1-хк)\(у]-у1), (8)

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

Тейлора по расстоянию от рассматриваемого вихря до центра кластера [3]. Отбрасывая

где А - интерполяционная функция; Г* - интенсивность перераспределяемого вихря, находящегося в произвольной точке с координатами (хк, у); Г(./} - циркуляция, получаемая вихрем в ячейке (/,_/) от перераспределяемого вихря. Регулярность поля завихренности необходима для подшага диффузии и расчета скорости с использованием быстрого преобразования Фурье. При этом количество вихревых частиц не растет с течением времени, как в бессеточном методе. В работах [5], [6], [7] дается обзор различных интерполяционных формул, применяемых в вихревых методах. В качестве интерполяционных функций А используются В-сплайны различных степеней. Ошибка такой интерполяции уменьшается с ростом степени используемого полинома. В данной работе используется так называемая формула М4' [6], которая является формулой третьего порядка:

АГ'(х) =

5 .2 3 .3

1--X + — X , • .

2 2 <1

(\-Х')(2-Х')2/2, \<Х'<29 (9)

О, х" >2

где х* =|х|//г. В работах [5], [6], [7] было

показано, что схемы Эверетта и М4' имеют наименьшую «численную диффузию» и сохраняют моменты завихренности. Формула (9) вносит наименьшую «численную диффузию» в численный метод [6], [7], и по этой причине она используется.

Для проверки алгоритма подшага конвекции рассматривается задача о диффузии двумерного вихря Ламба-Озеена, которая для вихревого поля имеет аналитическое решение:

- ехр[ -г2/(4\>т) ]

471VI

Т = ( +

£ 4у '

(Ю)

Здесь г = х-хс

, Г - интенсивность

вихря, Хс - радиус-вектор центра и р - «радиус ядра» вихря Ламба-Озеена.

При моделировании процесса чистой конвекции вихрь Ламба-Озеена превращается в экспоненциальный вортон с распределением завихренности, получаемым из (10) при / = 0. При этом его центр за время Г = 2,0 переместится в направлении вектора скорости набегающего потока ик = 1,0 на расстояние

= 2,0 . Распределение завихренности внутри такого вортона при перемещении в системе координат, связанной с его центром, не меняется.

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

1

т.

р р

со (хр,г)-<йа(хрл)

(И)

где /ир - число «вихрей в ячейке», попавших в ядро вихря Ламба-Озеена.

Результаты численных экспериментов, представленные на рис. 1, позволяют отметить особенности работы применяемого метода. Ошибка этапа конвекции складывается из двух частей: ошибки метода численного интегрирования системы (3) и ошибки от перераспределения интенсивности вихрей в ячейки сетки (8). Рассмотрим случай поступательного (без вращения) перемещения вортона (Г0 = 1,0, р = 0,1) за один расчетный шаг А? = 2,0 сразу в конечную точку. В результате

получается погрешность 6р = 4 = 4,6 • 10"" (сетка 400x200) и б"4'= 9,1-Ю-7 (сетка

0,23 0,20 0,18 0,15 0,13 0,10 0,08 0,05 0,03 0,00 1,00Е-06

А

/

* / /

/

/

/

/

//

1,00Е-05 1.00Е-04 1.00Е-03 1.00Е-02

М

Рис. 1. Изменение погрешности решения на этапе конвекции от шага по времени Расчетная область 4x2, конечное время I = 2,0, их = 1,0: □ - сетка 400x200 при Г0 = 3,0, (3 = 0,2; ' Х- - сетка 400x200 и -(> - сетка 800x400 при Г0 = 1,0, (3 = 0,1; -д- - сетка 400x200 и -+- - сетка 800x400 при Г0 = 1,0, (3 = 0,2

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

значение шага по времени А^р', при котором погрешность этапа конвекции метода «вихрь в ячейке» на конечном интервале времени

будет минимальна. Величина А^р' несколько уменьшается при увеличении разницы между максимальными и минимальными значениями завихренности вихревого поля (рис. 1). Так для вортона Г0 = 1,0, р = 0,2 оптимальный шаг Д?,0*" «Ю-4, а для Г0 = 1,0,

Р = 0,1 оптимальный шаг А^1" « 5 -10"5. Для случая Г0 = 1,0, р = 0,2 при Д/ = 10"3 и Д/ = 10"5 получается одинаковый порядок ошибки (рис. 1), хотя время расчета увеличивается в 100 раз. По этой причине на этапе конвекции уменьшение шага по времени дальше

величины А^"' является нецелесообразным. 3. Подзадача диффузии

3.1. Схема метода «донор-акцептор» в двумерной постановке задачи

На диффузионном подшаге метода «вихрь в ячейке» для решения ДУЧП (4) применяется метод «донор-акцептор» [1,5]. Для удобства вместо завихренности «вихря в ячейке» со/; в вихревых методах используют интенсивность вихря (циркуляцию)

Гр((,х) = И2а>р(1,х). Здесь применяется

сквозная нумерация вихрей и индексу р вихря соответствует пара индексов (г,у) для хранения его интенсивности в двумерной матрице. Математическая формулировка метода «донор-акцептор» [5] имеет следующий вид:

ЧеМр

г// + д/; = г//;+ X (гщ0)в'„(х)х

ч

*&„(У)-Тр«Х?9(х)СГт(у)), (12)

где

ег/

г+к/2-2п

Р_ч_

-ег/

р_ч

л/

(13)

Здесь 2-х или ур - координаты центра ячейки вихря р, ег/{г) - интеграл вероятности. Вследствие характера функции распределения Гаусса суммирование в (12) ведется в некоторой окрестности Мр ячейки р [4], которая далее будет называться «диффузионной молекулой». «Диффузионная молекула» Мр представляет собой квадрат со стороной (2п +1)/г, где п - количество соседних ячеек в одном из направлений системы координат. Коэффициенты (13) не зависят от времени и для «диффузионной молекулы» заданного размера вычисляются один раз.

Наилучшие результаты численного моделирования диффузии вихря Ламба-Озеена методом Д-А получаются [5], если шаг по времени выбирается равным

Д/ = кор' А2 / V.

(14)

Численное моделирование задачи при

V > 10"4 показало [5], что величина кор' для оптимального процесса диффузии зависит только от размера п «диффузионной молекулы». Далее для сокращения времени вычислений используется молекула минимального размера п = 1. Величина коэффициента оптимального шага по времени находится в пределах 0,20 <кор' <0,21 для погрешности 18 | □ < 0,02, полученной в конечный момент времени счета г = 3,75 [5]. Здесь 8 - относительное отклонение численного решения задачи от аналитического для поля завихренности в центре вихря Ламба-Озеена (10)

8 = Гсо-соет;/сое1.

(15)

Для дальнейшего повышения точности необходим подбор кор'.

3.2. Модернизированная схема метода «донор-акцептор» (М-ДА)

Применим метод Пуанкаре [8] к уравнению (4), рассматривая его как возмущенное уравнение соответствующего однородного уравнения, которое соответствует невязкой задаче (условие сохранения завихренности). Решение можно записать в виде

о = со 0 + sco, + 0( s2),

если порождающее решение со 0 существует [8]. При этом пренебрегаем членами порядка 0(s2) и выше. Здесь со t - возмущенное решение, s - малый параметр. Полагая для удобства

е = v / s.

(17)

подставляя (16) и расщепляя уравнение (4) по степеням s, получим

ôco0 Эсо, —— = 0. —L = е2Дсо0,

dt

dt

(18)

где е, - некоторая постоянная.

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

Д t = kop'h2/zT

(19)

xG'jj)-rp(t)G'ji)Glp(j)),

(20)

где

Erf

Ч+0-5-О

2^

(16) ~Erf

f \ i -0.5-г

Р ч

(21)

Если шаг Д? задается, например, схемой этапа конвекции, то величины е, и г находятся последовательно из (19) и (17). При этом проверяется, чтобы величина е действительно была малой. В данной работе использовалось ограничение: е < 0,1.

Уравнения численной схемы метода М-ДА с учетом выражений (17) - (19) записываются в следующем виде:

Г// + М) = Гр(1)+Е X

Коэффициенты (21) зависят только от взаимного положения вихрей друг от друга, если шаг по времени выбирается в виде (19) для схемы МД-А или (14) для схемы Д-А.

3.3. Верификация схемы метода М-ДА

В качестве тестовой используется задача о диффузии двумерного вихря Ламба-Озе-ена (10). Результаты моделирования представлены в табл. 1 и на рис. 2, 3, погрешность решения определялась по формуле (15).

Наблюдаемый рост погрешности решения при уменьшении расчетного шага Д/ с 10 ~3 до 10 ^ для случая V = 10 ~5 (табл. 1) объясняется:

а) возрастанием в 10 раз количества расчетных шагов, при этом ошибка на конечном интервале времени складывается из ошибок каждого шага по времени;

б) уменьшением в 10 раз (с 2,2-10 4 до 2,2-10~5) величины вихревой интенсивности, передаваемой ближайшему вихрю в «диффузионной молекуле». При этом изменение вихревой интенсивности на величину меньше, чем 10 "8, не будет фиксироваться из-за ограниченности количества значащих цифр в машинном представлении числа (использовалась одинарная точность, 8 значащих цифр). При очень малых Дг и V величина передаваемой соседнему вихрю вихревой интенсивности в «диффузионной молекуле» может содержать только одну значащую цифру и даже обнуляться. Для обеспечения достаточной точности необходимо иметь запас хотя бы в три ячейки матрицы представления числа в ЭВМ. Чтобы повысить точность метода МД-А при моделировании задач диффузии с V < Ю-5 и Дг < ЮЛ необходимо использовать

Рис. 2. Распределение завихренности при моделировании диффузии вихря Ламба-Озеена методом МД-А

Г0 = 1,0, Р = 0,1, V = Ю-3, / = 2,0, ЬхН =2x2.

Х- М = 103 г< □ - Аг = 104 для сетки 200x200 вихрей;

+ - А1 = 5-10 4 и О - А1 = 10 * для сетки 400x400 вихрей; штриховая линия - аналитическое решение в момент времени 1 = 0, сплошная линия - аналитическое решение в момент времени X = 2,0

СО

Рис. 3. Распределение завихренности при моделировании диффузии вихря Ламба-Озеена методом МД-А

Г0 = 1,0, р = 0,1, V = Ю-4,1 = 2,0, ЬхН =2x2.

X - А1 = Ю-3 и □ - А1 = Ю-4 для сетки 200x200 вихрей;

+ - А1 = 10 3 и О - А( = /0 4 для сетки 400x400 вихрей; штриховая линия - аналитическое решение в момент времени 1 = 0, сплошная линия - аналитическое решение в момент времени г = 2,0

переменные с двойной точностью (16 значащих цифр), а также более точно подобрать коэффициент к°р' (точности в два знака после запятой к°р' = 0,21 может оказаться недостаточно).

Анализируя данные, представленные на рис. 1, 2, 3 и в табл. 1, можно сделать вывод,

что метод МД-А (20) позволяет получать хорошие результаты в широком диапазоне изменения размера ячейки сетки /г, коэффициента кинематической вязкости V и шага по времени Д?. Метод Д-А дает аналогичные результаты только при двух значениях шага А? (14), что не всегда удобно на практике [5].

Таблица 1. Влияние изменения параметров численной схемы метода МД-А на погрешность решения в задаче о диффузии вихря Ламба-Озеена

V ъ Д1 8 5

10"3 0,01 ю-3 4,76-10"2 0,0089

ю-3 0,01 Ю-4 4,76-10"3 0,0087

ю-3 0,005 5-10^ 9,52-Ю-2 0,0080

Ю-3 0,005 ю-4 1,90-10"2 0,0084

ю-4 0,01 ю-3 4,76-10"3 0,0013

ю-4 0,01 ю-4 4,76-10^ 0,0013

ю-4 0,005 ю-3 1,90-10"2 0,0012

ю-4 0,005 ю-4 1,90-10~3 0,0011

10~5 0,01 ю-3 4,76-10^ 0,00014

10~5 0,01 10"4 4,76-10"5 0,00054

ю-5 0,005 ю-3 1,90-10"3 0,00011

ю-5 0,005 ю-4 1,90-Ю"4 0,00068

4. Моделирование конвекции-диффузии вихря Ламба-Озеена

Рассмотрим применение метода «вихрь в ячейке» к задаче конвекции - диффузии вихря Ламба-Озеена (10), (11). В качестве критерия ошибки получаемого решения используется величина 5р (11).

Сравнивая численные результаты (рис. 4, 5 и табл. 2), полученные при использовании схем Д-А и МД-А в методе «вихрь в ячейке», с аналитическим решением задачи, можно сделать следующие выводы.

При уменьшении коэффициента кинематической вязкости V величина изменения

завихренности вследствие диффузии может становиться меньше «численной диффузии», возникающей на подшаге конвекции. Это приводит к необходимости измельчения сетки и уменьшения шага по времени численного моделирования. Как видно из представленных данных, если для численного моделирования движения вихря с Г0 = 1,0, (3 = 0,1 при V = 10 "3 было достаточно сетки 400x200 и шага Дг = 10 "4, то при V = 10 4 для достижения такого же порядка погрешности необходима сетка 800x400 и = 5-10 5. Расчетное время задачи при этом возрастает больше, чем в 8 раз.

V ь Д-А МД-А

дг 5В дг 5В

ю-1 0,04 3,36-10"3 0,0062 3,36-Ю"4 0,0079

ю-2 0,02 8,40-10"3 0,036 8,00-10 0,015

ю-2 0,01 2,10-Ю"3 0,0044 2,10-10^ 0,017

ю-3 0,01 2,10-Ю-2 0,56 5,00-Ю"4 0,024

Ю"4 0,53 ю-4 0,0077

- - 0,014

ю-4 0,01 ю-4 0,035 ю-4 0,024

0,005 10"4 0,034 10"4 0,019

ю-4 0,005 5,00-Ю-5 0,045 5,00-Ю-5 0,0097

ю-4 0,0025 ю-4 0,032 ю-4 0,019

ю-5 0,005 ю-4 0,020 10"4 0,023

Таблица 2. Погрешность численного метода в зависимости от применяемых методов расчета диффузии и изменения параметров моделирования

Рис. 4. Распределение завихренности при моделировании конвекции-диффузии вихря Ламба-Озеена Г0 = 1,0, $=0,1, и= Ю-3, I = 2,0, ЬхН=4x2, сетка 400x200 вихрей.

Х-МД-А А1 = 5,0-10 ', +-МД-А А( = 10',О-МД-А А( = 10'5,

П-Д-А М = 10-4,А-Д-А А1 = 2,1-Ю-2, штриховая линия - аналитическое решение в момент времени I = О, сплошная линия - аналитическое решение в момент времени I = 2,0

Рис. 5. Распределение завихренности при моделировании конвекции-диффузии вихря Ламба-Озеена Г0 = 1,0, Р = 0,1,\ = 10',Х = 2,0, ЬхН=4x2, сетка 800x400 вихрей.

О - МД-А Л1 = 10 ', +-МД-А А1 = 5, О Ю -5,

А-Д-А А( = 10-', П-Д-А М = 5,0-10'5, штриховая линия - аналитическое решение в момент времени / = О, сплошная линия - аналитическое решение в момент времени 1 = 2,0

При стремлении шага по времени к нулю численное решение задачи методом «вихрь в ячейке» со схемой Д-А стремится к начальным условиям задачи (10), при применении МД-А оно стремится к нестационарному решению вязкой задачи (рис. 4, 5). Численная схема метода Д-А накладывает гораз-

до более жесткое требование на расчетную сетку (14), которое сложно выполнить при малых V. По этой причине при значениях коэффициента кинематической вязкости у<10"3 более целесообразно применение схемы МД-А на подшаге диффузии в методе «вихрь в ячейке».

Схема МД-А также сохраняет точность при величине коэффициента кинематической вязкости V = 10"5 на сетке 800x400 с шагом A/ = 10"4 для задачи диффузии (табл. 1). При таких же данных схема Д-А теряет чувствительность, так как коэффициенты (13) пренебрежимо малы. Однако «численная диффузия» на этапе конвекции уже на порядок превышает физическую диффузию (табл. 2), что и объясняет одинаковый порядок ошибки при использовании схем Д-А и МД-А при v=10"5. Такая же ошибка 5р « 0,02 получается и при отсутствии этапа диффузии для рассматриваемого случая (рис. 1). Для повышения точности метода «вихрь в ячейке» при v < 10 "5 необходимо дальнейшее измельчение сетки или использование на этапе конвекции других численных схем, имеющих большую точность.

Список литерату ры

1. Basin М., Kornev N. Beruecksichtigung der R eibung in der Wirbelmethode, ZAMM, 1998, vol. 78, N5. Pp. 335-344.

2. Greengard L., Rokhlin V. A Fast Algorithm for Particle Simulations, J. Comput. Physics,

1987, vol. 73. Pp. 325-348.

3. Kornev N., Leder A., Mazaev K. Comparison of Two Fast Algorithms for the Calculation of Flow Velocities Induced by a Three-Dimensional Vortex Field, Schiffbauforschung, 2001, vol. 40, N 1. Pp. 47-55.

4. Taranov A., Kornev N., Leder A. Development of the С omputational Vortex Method for Calculation of Two-Dimensional Ship Sections with Flow Separation, Schiffbauforschung, 2000, vol. 39, N 2. Pp. 95-105.

5. Nikonov V., Kornev N., Leder A. The Ratio between Spatial and Time Resolutions for the Diffusion Substep in 2D Computational Vortex Methods, Schiffbauforschung, 2002, vol. 41, N 3/4. Pp. 5-12.

6. Cottet G.-H., Koumoutsakos P.D. Vortex Method: Theory and Practice, Cambridge University Press, 2000. - 320 p.

7. Monaghan J.J. Extrapolating B-Splines for Interpolation, J. Comput. Phys., 1985, vol. 60. Pp. 253-262.

8. Моисеев H.H. Асимптотические методы нелинейной механики. М.: Наука, 1981. -400 с.

COMPUTING VORTICITY DIFFUSION WITH A MODIFIED "DONOR-ACCEPTOR" SCHEME AND THE INTEGRATION OF THE SCHEME INTO THE "VORTEX-IN-CELL" METHOD

© 2003 V. V. Nikonov, V. G. Shakhov Samara State Aerospace University

A modernised scheme o f splitting Navier-Stokes equations for the numerical method of "vortex-in-cell" is proposed. The method consists in introducing a minor parameter proportional to kinematic viscosity into the scheme. It allows finding the solution of the problem in the form of asymptotic series. An advection-diffusion problem of the Lamb-Ozeen vortex is considered for the verification of the proposed method. Computation results obtained by "donor-acceptor" (D-A) and MD-A methods are compared with an analytical solution of the advection-dififusion problem. The MD-A method is shown to be more suitable for low viscosity flow simulation. Unlike the D-A method, the MD-A method remains accurate in a wide range of space and time resolution.