Введение

Метод решетчатых уравнений Больцмана (LBM) широко используется при моделировании процессов течения жидкости и воздушных потоков, а гибкость задания граничных условий позволяет учитывать сложную структуру среды.

При исследовании сред со сложной геометрической структурой современное применение метода охватывает такие проблемы как моделирование реактивных транспортных процессов в нефтегазовых резервуарах [1]; расчет потока жидкости и теплопередачи в средах со сложной геометрией, например, металлической пены [2]; моделирование взаимодействия волн и пористых структур [3] и др.

Наряду с методом решетчатых уравнений Больцмана для моделирования фильтрации потока жидкости через пористую среду также используются уравнения Навье-Стокса (NSE) и модели поровой сети (PNM). Уравнения Навье-Стокса решаются на основе метода конечных элементов, конечных объемов или конечно-разностных схем.

В [4] предлагается алгоритм, использующий модифицированные уравнения Навье-Стокса для моделирования фильтрации жидкости в анизотропной пористой среде, в [5] проводится исследование зависимостей результатов, полученных решением уравнений Навье-Стокса и закона Дарси для пористой среды.

Модели поровой сети [6] представляют реальное пространство множеством сферических пор, соединенных цилиндрическими капиллярами. Гидродинамическая модель основывается на законе сохранения массы и на известных характеристиках течения Пуазёйля в цилиндрических капиллярах. К PNM, например, относятся алгоритм преобразования медиальных осей, максимального шара и сегментации по водоразделам.

Каждая из перечисленных групп моделей имеет свои преимущества и недостатки. Уравнение Навье-Стокса всегда дает решение для потока жидкости, находящегося в стационарном режиме, в то же время в LBM поток переходит в стационарный режим движения жидкости после большого числа итераций (10 000…20 000, [6]), поэтому LBM может работать медленнее, чем NSE.

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

Также отметим, что в сравнении с PNM LBM зачастую не способен проанализировать всю структуру порового объекта, а алгоритмы PNM не имеют подобных ограничений.

Методология

Общие сведения о методе решетчатых уравнений Больцмана. LBM представляет собой дискретный аналог уравнения Больцмана [7]. В соответствии с LBM молекулы могут находиться только в узлах решетки, причем число направлений из одного узла в другой ограничено. Обычно решетка разбивается на клетки, имеющие форму квадратов или кубов.

Будем использовать модель D2Q9, которая имеет 2 измерения и 9 возможных направлений скорости. На рис. 1 изображено расположение решетки в декартовой системе координат со скоростями c_i, где i = 0,1…8 индекс направления и c₀ = 0 соответствует неподвижной молекуле.

Рис. 1. Модель D2Q9

Fig. 1. Model D2Q9

Функции распределения f(x, u, t) описывают состояние системы каждого узла решетки. Функция распределения определяет [6], какая часть частиц расположена в окрестности точки x(x, y) в момент времени t с координатами от x до x+Δx, от y до y+Δy и диапазоном скоростей от uux,uy  до uux+Δux,uy+Δuy .

Все методы LBM основаны на принципе разделения этапов столкновения молекул (collision step) и потоковой передаче (streaming step).

На первом этапе [8] происходит релаксация функции fi (x, t) в направлении равновесной функции распределения fieq (x, u, t).

Во время потоковой передачи каждое распределение частиц перемещается в соответствии со скоростями c_i в соседние узлы.

Объединив этапы столкновения и распространения молекул получим формулу (1):

fix+ci∆t, t+∆t=fix,t+Ωi,                                               1

где Ωi  – оператор столкновений.

Равновесная функция распределения [7] определяется по формуле (2):

fieqx=wiρx1+3ci∙ucs2+92(ci∙u)2cs4-32u2cs2,                         2

где веса w_i = 4/9 для частиц находящихся в покое, 1/9 при i = 1, 2, 3, 4 и 1/36 при i = 5, 6, 7, 8; c_s – скорость звука, принимается равной 1/3 ; ρ(x) – плотность жидкости в точке x; u – скорость потока жидкости в точке x. 

Макроскопические параметры скорости и плотности жидкости могут быть получены вычислением нулевого и первого момента функции распределения, то есть:

ρ=ifi,ρu=icifi.                                                        3

Схемы столкновений с двумя параметрами релаксации (TRT). Самой широко применяемой схемой столкновений является оператор BGK в совокупности со стандартными граничными условиями упругого столкновения. Из-за того, что оператор использует лишь один параметр релаксации, [9] появляется вычислительная нестабильность и наблюдается наличие зависимости расположения границ от вязкости жидкости.

С другой стороны, в [9] было показано, что необходимо как минимум два параметра релаксации, чтобы предотвратить наличие нелинейных ошибок, обусловленных вязкостью, при применении граничных условий типа упругого столкновения.

В связи с этим, в основе разрабатываемой модели будет лежать оператор столкновений с двумя параметрами релаксации (TRT-оператор). TRT-оператор имеет [10] параметр ω⁺, который задает вязкость жидкости, и свободный параметр ω^-. Стабильность и точность оператора обусловлена так называемым «магическим параметром» Λ :

Λ=1ω+∆t-121ω-∆t-12.                                                    4

Существуют  известные  значения Λ , которые обеспечивают определенную степень точности,  например, Λ  = 1/12  отвечает   третьему  порядку  точности,  а Λ  =1/6 – четвертому, а Λ  = 1/4 зачастую стабилизирует вычисления.

Модель TRT-оператора основана [10] на декомпозиции ожидаемого фазового пространства молекул на симметричные и ассиметричные части. В методе решетчатых уравнений Больцмана для каждой скорости c_i существует ci=-ci  и таким образом множество скоростей c_i обладает симметрией. По аналогии можно переписать формулы для расчета плотностей распределений:

fi+=fi+fi2,                   fi-=fi-fi2,

fieq+=fieq+fieq2,                   fieq-=fieq-fieq2.

Тогда стандартная TRT модель принимает вид:

fi*=fi-ω+Δtfi+-fieq+-ω-Δtfi--fieq-,                                (5)

fix+ciΔt,t+Δt=fi*x,t.

И кинематическая вязкость жидкости:

ν=cs21ω+Δt-12.                                                          (6)

Граничные условия. Будем считать, что исследуемая область жидкости имеет прямоугольную форму. Тогда на левой и правой границе применим периодические граничные условия, в верхней части зададим вектор скорости, направленный в область жидкости (Zou and He boundary conditions, например, описаны в [11]), в нижней условия нулевого градиента по функции распределения.

Внутри исследуемой области в качестве граничных условий используется видоизмененный half-way bounce-back (описание классического алгоритма приведено, например, в [10]).

Введение внешних сил в уравнение Больцмана. Учтем воздействие внешних сил на молекулы жидкости, добавив элемент ΔtFi  в исходное уравнение (1) [12]:

fix+ciΔt,t+Δt=fix,t+Ωi+ΔtFi.                                      (7)

Чтобы определить значение Fi  рассмотрим схему определения внешних сил, основанную на LGA ( Lattice gas automata – метод, предшествующий LBE):

Fi=ωiF∙cics2,                                                                       (8)

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

iFi=0,               iFici=F.                                                 (9)

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

Для удобства определим [13] переменную gi=ΔfiΔt , задающую изменение функции распределения вследствие влияния внешней силы и рассчитываемую по формуле:

gi=ωiρgciycs2                                                                  (10)

где g  – ускорение (обусловливающее фиктивную силу тяжести), имитирующая действие градиента давления в направлении оси OY, а действующая сила представлена выражением F=ρg .

Внесение изменений в этап распространения частиц. Этап потоковой передачи зачастую включает в себя граничные условия вида упругого столкновения (bounce-back).

Будем применять half-way bounce-back, учитывая частичную непроницаемость узлов.

Введем матрицу ε , содержащую информацию о проницаемости каждого узла решетки, в таком случае εx∈[0, 1]  – где x указывает позицию узла; εx=0,  если узел абсолютно жидкий и εx=1,  если жидкость не может оказаться в точке x (например, в случае если узел полностью состоит из частицы грунта).

Объединим формулы упругого столкновения с непроницаемым узлом и распространения частиц в соседние узлы.

Взаимодействие узлов охарактеризуем с позиции функции вероятности, имеющей направление «в» узел, т.е. с позиции узла с координатами x+ciΔt , а не x.

Пусть в момент времени t частицы перемещаются из узла xj  в x­­­i  по направлению k (рис. 2).

Рис. 2. Момент времени t

Fig. 2. Moment of time t

Тогда во время t+Δt  в узел xj  вернется εxi∙fxjk  часть частиц, двигающихся в направлении k , с другой стороны в момент t+Δt  в xj  из x­­­i  переместится 1-εxj∙fxik  частиц (рис. 3, 4).

Рис. 3. Момент времени t

Fig. 3. Moment of time t

Рис. 4. Момент времени  t+Δ t

Fig. 4. Moment of time t+Δt

Таким образом, формула fix+ciΔt,t+Δt=fi*x,t  видоизменяется следующим образом:

fix+ciΔt,t+Δt=εx∙fi*x+ciΔt,t+1-εx+ciΔt∙fi*x,t.                  (11)

Сравнение с аналогичными уравнениями. Пористую структуру среды можно описать двумя способами:

1. Ввести промежуточный шаг porous medium step [7]:

Обозначим fi** , рассчитанной для этапа столкновения молекул

fi**x,t+Δt=fi*x,t+Ωi*.

Тогда porous medium step принимает вид:

fix,t+Δt=fi**x,t+Δt+nfi**x+ciΔt,t+Δt-fi**x,t+Δt,

или сгруппируем fi**x,t+Δt:

fix,t+Δt=fi**x,t+Δt∙1-n +n∙fi**x+ciΔt,t+Δt,             (12)

где макропараметр n – коэффициент пористости среды.

2. Считать каждый узел решетки частично непроницаемым, как например, в граничных условиях вида partially saturated bounce-back [8].

Выражение (11) схоже по структуре с (12), однако коэффициенты εx  учитываются для каждого узла, что характерно для partially saturated bounce-back. В то же время коэффициенты εx  не зависят от параметра релаксации и в уравнении не используется дополнительный оператор столкновения. Кроме того, в (11) применяется TRT-оператор, в то время как в partially saturated bounce-back BGK-оператор.

Результаты

Проверим, выполняются ли законы сохранения массы и импульса для выведенного уравнения.

Чтобы удостовериться в их справедливости проверим эквивалентность формул для нахождения плотности жидкости и скорости потока [12, 13]:

ρ=ifi=ifieq     и      ρu=icifi=icifieq.                                 (13)

Тогда достаточно убедиться, что:

ifi-fieq=0      и      icifi-fieq=0.                                      (14)

Выполнение законов будем проверять эмпирическим [14] путем на 8000 итераций:

Найдем нормы Фробениуса для матриц ifi-fieq  и ici(fi-fieq)  на каждой из итераций (рис. 5 – 7).

Из рисунков следует, что нормы матриц не превышают значений 8e-15 для любой итерации, т.е. соблюдаются законы сохранения массы и импульса.

Рис. 5. График зависимости iciy(fi-fieq)F  от номера итерации

Fig. 5. Graph of dependence iciy(fi-fieq)F  from the iteration number

Рис. 6. График зависимости icix(fi-fieq)F  от номера итерации

Fig. 6. Graph of dependence icix(fi-fieq)F   from the iteration number

Рис. 7. График зависимости ifi-fieqF  от номера итерации

Fig. 7. Graph of dependence ifi-fieqF   from the iteration number

Покажем, что ux,n-ux,n-1F→0,   uy,n-uy,n-1F→0,  ρn-ρn-1F→0   при n→∞  (рис. 8 – 10).

Рис. 8. График зависимости uy,n-uy,n-1F  от номера итерации n

Fig. 8. Graph of the dependence of uy,n-uy,n-1F   on the iteration number n

Рис. 9. График зависимости ux,n-ux,n-1F  от номера итерации n