Конечно-элементное моделирование локальных процессов переноса в пористых средах

УДК 532.51+536.24

Ю. И. Димитриенко, А. И. Левина, П. Боженик

КОНЕЧНО-ЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЕ ЛОКАЛЬНЫХ ПРОЦЕССОВ ПЕРЕНОСА В ПОРИСТЫХ СРЕДАХ

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

Процессы, происходящие в материалах на микро-, мезо- или на-ноуровнях, объединяют общим понятием — локальные процессы. Эти процессы во многом определяют макроскопические свойства материалов. Так, газопроницаемость пористых материалов определяется особенностями локального переноса на уровне отдельной поры. В отличие, например, от упругих локальных свойств, хорошо изученных в микромеханике композитов [1-4], локальные процессы переноса значительно менее изучены. В работах [5-7] был предложен метод асимптотического осреднения для исследования газодинамических процессов в периодических пористых средах. Возникающие в этом методе задачи на ячейке периодичности впервые были решены в работах [610] приближенно-аналитическими и конечно-разностными методами. Целью настоящей работы является разработка конечно-элементного метода решения задач на ячейке периодичности и исследование с его помощью особенностей локальных процессов переноса в пористых средах, а также изучение погрешностей численного решения, вызванных особеностями самого метода конечных элементов.

Формулировка локальной задачи. Рассмотрим пористую среду периодической структуры, поры которой заполнены линейно-вязким совершенным газом (рис. 1). Обозначим 10 — линейный размер ячейки периодичности V среды и х0 — характерный размер всей области V среды. Положим, что выполняется соотношение к = 10/х0 ^ 1, где к — малый параметр, и введем локальные и глобальные безразмерные координаты £1 = ХХ/к и XX = хг/х0. Тогда, согласно методу асимптотического осреднения (МАО) [9], движение газа в отдельной поре (в ячейке периодичности (ЯП)) описывается системой уравнений

Навье-Стокса, которая в криволинейных координатах X1, связанных с локальными декартовыми координатами , имеет следующий вид:

V* • и(0) =0, е е ^;

- V*p(1) + ^Аи(0) = Vжp(0), Г е ^д; и(0) = 0, Г е

[ [и(0)] ] = 0, [ [р(1)] ] = 0, <р(1)> =0.

Здесь и(0) — вектор скорости газа в нулевом приближении МАО; р(0) и р(1) — давления в нулевом и первом приближении; ^ — коэффициент вязкости газа; V* f = г1дf/дX1 = е1 дf/д^l — набла-оператор дифференцирования по локальным координатам; Vxf=el дf/дxl — набла-оператор дифференцирования по глобальным координатам; г1 = еР дХ1 /д^ — локальные векторы взаимного базиса; еР — ор-тонормированный декартов базис; ^} = (1/<£д) / fdV*g — оператор

осреднения по ячейке периодичности; = / dV*g — пористость.

V*g — область ячейки периодичности, заполненная газом (пора), а — поверхность раздела газа и твердой фазы в ЯП (твердая стенка). В задачу (1) входят условия периодичности функций на противолежащих границах ЯП, которые обозначаются как ] ] =0.

Задача (1) рассматривается относительно четырех неизвестных функций локальных координат и(0)(£Р) и р(1)(£Р), заданных в V*g, причем вся область ЯП V* в координатах представляет собой куб: -1/2 < < 1/2. Давление р(0)(ХР), зависящее только от глобальных координат, в задаче (1) полагается заданным, поэтому решение задачи (1) параметрически зависит от р(10)(ХР).

Вследствие наличия интегрального условия (р(1)) = 0 задача (1) имеет интегродифференциальный тип, что в совокупности с наличием условий периодичности делает ее непосредственное решение достаточно сложным даже для современных численных методов. Отметим также, что функция р(1) фактически представляет собой пульсацию давления по отношению к нулевому приближению р(0) (согласно МАО, полное давление определено соотношением кр(1) = р — р(0) + о(к), где последнее слагаемое означает величины малые по сравнению с к), поэтому функция р(1) может быть как положительной, так и отрицательной, что также не характерно для классических задач газовой динамики.

Сведение к локальной задаче на 1/4 ЯП. Будем рассматривать далее только двумерную осесимметричную пористую структуру (рис. 1),

Рис. 1. Пористая периодическая структура с ячейкой периодичности

у которой поры образуют продольные каналы вдоль оси 0£3. Такая структура возникает, например, в композиционных однонаправленных армированных материалах, обладающих продольными порами — трещинами. Выберем в качестве криволинейных координат Xг цилиндрические координаты X1 = г, X2 = 9, X3 = г, в этом случае решение локальной задачи (1) зависит только от двух координат — г и г.

В силу линейности локальной задачи (1) ее решение ищем в виде линейной функции входных данных, т.е. градиента Ур(0) = р,3е3, где Р,3 = др/дХ3 (остальные компоненты градиента Ур(0) в силу однока-нальной структуры равны нулю, так как отсутствует сквозное течение по всей пористой структуре по направлениям Ох1 и Ох2):

3

Подставляя выражения (2) в локальную задачу (1), после исключения градиента Ур(0) получаем

V- W = 0;

-VP + ^Д • W = ег, в

W

< р >= 0, [^]]=0, [ [р] ] = 0

— локальную задачу для определения давления Р и вектора скорости W = "гег + "ег, где ег, ег, ев — векторы физического базиса цилиндрической системы координат X\

Система (3) представляет собой стационарную задачу течения некоторой фиктивной линейно-вязкой несжимаемой жидкости. Положим, что структура ЯП имеет зеркальную симметрию относительно плоскостей 0£ 1£2 и 0£2£3, тогда, как было показано в [10], решение Р, W задачи (3) может быть построено с помощью симметричного и антисимметричного продолжений функций Р, , определенных в 1/4

ячейки периодичности (в первом квадранте У = г : 0 ^ £г ^ 1/2}) и являющихся решениями краевой задачи

V- ^ = 0;

Р + = ег, в %;

£: 1 = 1 = 0,

5 д _

г = 0 : = 0, 1 = 0, Р/2 = 0;

г = 0 : \¥т/т = 0, 1 = 0, Р = 0; г = 1/2 : 1 = 1 = 0, Р = 0,

где У^д = У П У^д — область течения, а = У П — твердая стенка в 1/4 ЯП.

Вариационная формулировка локальной задачи. Введем обозначения для тензора напряжений Т и тензора скоростей деформаций D

Т = — РЕ + ^Б, Б = ^5 <8> 1р + V << 1),

2

имеющих в осесимметричном случае следующие ненулевые компоненты Ттт, Тгг, Твв, Ттг и Отт, , , Втг, определяемые как

01РТ ЗШг 1¥г _ 1/ 01РТ

ДГГ = —^-, = —^-, Евв = -, = - --1--~-

ог ох г 2 \ ох ог

Тгг = — Р + Н-Етт, Т22 = — р + , Твв = — Р + Н^вв, Ття = ,

где E — метрический тензор. Тогда уравнение равновесия в задаче (4) можно представить в дивергентном виде V • Т = ег. Если ввести кинематически возможное поле скоростей 8\У и вариацию давления 8р, то, умножая скалярно уравнение равновесия на 8\У, а уравнение несжимаемости в (4) на 8р, а затем интегрируя получившиеся выражения по области р, получаем систему двух вариационных уравнений

J Т • •8БдУд = ^ 5 • 81 вУд — У ег • 81вУд;

П П

/ V • 18рвУд = 0.

Здесь мы использовали формулы преобразования произведения [11] (V • Т) • = V • (Т • — Т • -8^ < ^, а также формулу Гаусса-Остроградского и ввели обозначения для вектора нормали к поверхности £д области Уд С У^д п и вектора усилий Б = п • Т. Заметим, что на всей границе области течения У^д в соответствии

с граничными условиями задачи (4) работа вектора поверхностных усилий обращается в нуль: Б • = 0 при г = 0, г = 0, 1/2. В цилиндрической системе координат система (7) имеет вид

/ (Тгг бБгт + Тгг + 2Trz + Твв ¿Дде =

Vg

¿Wz dVg;

v0

(Drr + Dzz + D00 )£pdVg = 0.

v0

Метод конечных элементов для решения вариационной задачи.

Для решения системы вариационных уравнений (8) применен метод конечных элементов с 6-узловым треугольным конечным элементом (КЭ), который отличается от классического 6-узлового КЭ [12] тем, что имеет 15 степеней свободы: по две компоненты вектора скорости И", в каждом узле и по одному значению давления Р = р в каждой вершине треугольника. Аппроксимация в каждом КЭ по скоростям — квадратичная, а по давлению — линейная:

{W} = [Ф] {q}, p = [Фр}т {y},

2x12 12

где {И} = (И1,И2) = (ИИГ, ИИ г) — координатный столбец скоростей в КЭ; {д} — координатный столбец скоростей в узлах; {у} — координатный столбец давлений в вершинах КЭ.

Матрица [Ф] и столбец {Фр} имеют следующий вид:

[Ф] = [Фав] =

2x12

Ф1 0 Ф2 0 Ф1 0

0 Ф3 0 Ф4 Ф2 0 Ф3 0

0 Ф5 0

Ф4 0 Ф5

Фб 0

0

Фб

{Фр}т = ( L1 L2 L3 ),

где

Ф1 = L1(2L1 - 1); Ф2 = L2(2L2 - 1); Ф3 = ¿3^3 - 1) Ф4 = 4L1L2; Ф5 = 4L2L3; Фб = 4L1L3;

1

L = 2s ^+Ь(г)Г+С(г)^;

a(1) = r(2)z(3) - r(3)z(2); b(1) = z(2) - z(3); c(1) = r(3) - r(2); 2Se = b(1)c(2) - b(2)c(1),

2

3

а г(¿), ^(г) — координаты вершин треугольника. Остальные коэффициенты а(г), Ь(г), С(г) получаются с помощью круговой перестановки индексов; 5е — площадь КЭ; Уд = 2п5е — объем КЭ.

Координатный столбец скоростей деформаций {Я}т = = (^1, Я2, Я3, Я4) = (Ягг, Д^, Бвв, \f2Drz) представим в виде {Я} = = [В]{д}, где

[В] =[{В(1)},..., {В(х)},..., {В(12)}];

4x12

{В(х)}т = (ЬХ1,ЬХ2,ЬХз, (^ + ^)/2);

LX = фах dHa + X V dHa.

Ha dXa HaHß dXe + ав HaH7 dXY.

X = 1... 12; а, в = 1, 2.

Вариации компонент вектора скоростей и тензора скоростей деформаций, а также давления р с помощью (9) принимают вид

¿{Ж} = [Ф]%}; ¿р = {Фр}т ¿{у}; ¿{Я} = [В]%}.

Введем координатный столбец напряжений {Т}т=(Т1, Т2, Т3, Т4) = = (Тгг , Тгг, Т^, -\/2Тгг), который на основании определяющих выражений (6) связан с координатным столбцом скоростей деформаций линейным соотношением

{Т} = - [Фр] {у} {Я},

cp

4 4x3 3 4

где [Фр] = {/0} {Фр}т — матрица, образованная диадным умножением

4x3 4 3

столбцов; {/о}т = (1, 1, 1, 0); {/*}т = (0, 1).

Введем также координатный столбец внешних усилий {5} = = (51, 52) = (вг, ), где вг, — компоненты вектора усилий S на поверхности тела. Тогда, подставив (9) и (14) в систему (8), получим систему линейных алгебраических уравнений (СЛАУ) для каждого отдельного КЭ

[К] М- [Вр] {у} = {/};

12x12 12 12x3 3 12

[Вр]т м = 0,

3x12 12

где обозначены следующие локальные матрицы жесткости и векторы

4

правой части:

[K] = [B]T [B])dVg; [Bp] = [ [B]T №P]dV; f = fs + /

12

Ш = / [Ф ]T {S}d£; |fs} = 1 [Ф]т {Iz }dVg

g

Vg

Численный алгоритм. На основе СЛАУ (15), записанной для одного КЭ, далее строилась СЛАУ для всей области (глобальная СЛАУ). При формировании глобальной СЛАУ суммарные векторы правых частей {/;}, соответствующие вектору поверхностных усилий S, обращаются в нуль, поскольку все внутренние стороны КЭ при суммировании перечисляются дважды, но направление вектора S при этом меняется из-за изменения направления вектора нормали на этих сторонах. На сторонах КЭ, которые принадлежат внешней границе области , работа вектора S, как отмечалось выше, обращается в нуль в силу граничных условий задачи, поэтому в глобальной СЛАУ {/е} = 0. Оставшиеся граничные условия задачи

£^sg : Wr = Wz = 0; r = 0 : Wz = 0; z = 0, 1/2 : Wz = 0, P = 0

присоединяются к глобальной СЛАУ.

Поскольку матрица СЛАУ в соответстии с (15) не является симметричной, то и матрица глобальной СЛАУ также не симметрична, для ее решения применялись прямой метод Гаусса и итерационный метод QMR.

Генерация конечно-элементной (КЭ) треугольной сетки для рассматриваемой области осуществлялась с помощью специализированного программного продукта [11], позволяющего с использованием компьютера задавать двумерные области в интерактивном режиме, а также генерировать нерегулярную КЭ сетку по методу Делоне, исходя из предварительно задаваемых ограничений площади и минимального угла КЭ: Se < Smax, > $min. Кроме того, в целях сопоставления результатов рассматривались регулярные КЭ сетки, построенные методом параметрического преобразования четырехугольника.

Результаты численного моделирования. Рассмотрим некоторые результаты численных расчетов. Вычисления были проведены для пористой одноканальной среды, у которой сквозные каналы пор имеются только в одном направлении (O£3). При расчетах варьировались следующие параметры: 1) тип КЭ сетки — регулярный и нерегулярный; 2) общее число КЭ; 3) неравномерность распределения КЭ; 4) радиус входного канала поры; 5) форма поры.

Рис. 2. Распределение скорости \¥г = ТТ в поре при различных значениях радиуса входного канала а для случая регулярной сетки

На рис. 2 показано поле распределения скорости Ж в поре при различных значениях относительного радиуса входного канала а = 0,1; 0,15; 0,2 для случая регулярной сетки. Параметр Ь канала — максимальный радиус поры — во всех трех случая был одинаков (Ь = 0,4). С увеличением радиуса а распределение скорости Ж по оси 0£3 поры становится более равномерным (уменьшается разность между максимальным и минимальным значениями скорости), при этом максимальное значение Ж, достигаемое во входном сечении, не зависит от а (рис. 3).

mio1

7т

\\

\л\ \ N \

\ \ Ч \ \ \ \ \ л-0,2 / / ' / / /

\ Ч ч 4^-0,15 ..... У / /

\ ч^я-о.?""-—

- 0,5 0 Ъ

Рис. 3. Распределение проекции скорости \¥г вдоль продольной оси канала при различных значениях радиуса входного канала а

При измельчении регулярной сетки КЭ (т.е. при уменьшении максимальной площади КЭ Smax), которое приводило к увеличению общего числа КЭ с 200 до 800 и соответственно к увеличению размерности СЛАУ, было установлено резкое увеличение погрешности вычислений. В результате абсолютные значения скорости \Уг возрастали примерно на три порядка при использовании метода Гаусса, и на два — при использовании метода QMR (рис. 4). Кроме того, искажался сам характер решения — появлялся максимум скорости в зоне максимального сечения поры, противоречащий физическим законам. Лучшие результаты были получены с применением нерегулярных КЭ сеток, вычислительная погрешность в этом случае уменьшалась примерно в три раза при использовании метода Гаусса и в полтора раза при использовании метода QMR (рис. 5). Нефизических максимумов при этом не возникало, а характер распределения скорости 1уг и других функций был такой же, как и для крупных сеток; число КЭ для нерегулярной сетки составляло 576 элементов. Практически полного устранения вычислительной погрешности удалось достичь при использовании неравномерной КЭ сетки (адаптивной сетки), у которой происходит сгущение узлов в зоне максимума решения; в рассматриваемой задаче это зона вблизи оси симметрии канала О£3 (см. рис.5,г). При использовании метода QMR для адаптивной сетки погрешность не превышала 10%. Все результаты сравнивались с решением, полученным конечно-разностным методом функций вихря и тока, описанным в работе [8]. При использовании адаптивной КЭ сетки и метода Гаусса погрешность снова возрастала до 10 раз (хотя характер решения сохранялся правильным и нефизические максиму-

»15.322

<15.322

<14.5158

<1 3.7095

<12.5033

<12.097

<11.2907

<10.4845

<0 67821

<8 87194

<8.06568

<7.25942

<6.45316 П <5.6469

<4 84083

<4 03437

<3 2281 1

<2 42185

<1.61558

<0.809321

>3.05326

<3.05328

<2.39257

<2.73188

<2.57119

<2.4105

<2.24981

<2.08912

<1.92843

<1.76774

<1.80705

<1.44636

<1.28567

<1.12498

<0.964288

<0.803598

<0.642907

<0.482217

<0.321528

<0.180838

Рис. 4. Распределение скорости Т'¥г = ТТ в поре для случая регулярной сетки с применением методов Гаусса

мы не появлялись) из-за большого числа КЭ — для адаптивных сеток оно составляло около 2000.

Для приложений важное значение имеет вопрос о влиянии геометрической формы поры на характер течения в ней газа. В работе [6] была высказана гипотеза о том, что для пор с узкими каналами существуют застойные зоны, в которых скорость движения газа практически нулевая, а основное движение осуществляется по тонким каналам, соединяющим застойные зоны. Там же было приведено экспериментальное подтверждение расчетов коэффициентов газопроницаемости пористых структур, вычисленных на основе этого допущения. Разработанный в настоящей работе метод решения локальных задач газовой динамики на ЯП дал возможность численного исследования влияния геометрии на характер течения газа в порах. На рис. 6

Рис. 5. Распределение скорости \¥х = \¥2 в поре для случая нерегулярной КЭ сетки с применением методов Гаусса

приведены результаты расчетов для продольной скорости \Уг в поре, у которой имеется полость, отделенная от сквозного канала перешейком. Вычисления показали (см. рис.5,г и 6,в), что действительно в таких порах формируется застойная зона, в которой скорость течения практически нулевая, а основное движение газа происходит по сквозному каналу, причем значение продольной скорости \¥г слабо зависит от геометрической формы застойной зоны и определяется только диаметром сквозного канала поры. На рис. 5 показано влияние типа КЭ сетки (адаптивная или равномерная), а также общего числа КЭ на точность расчетов. Как и ранее, использование адаптивной сетки в сочетании с методом QMR приводило к лучшим результатам даже для значительного числа КЭ (около 2000 элементов).

Выводы. В настоящей работе продолжены исследования, начатые в [5-10]. В результате дана вариационная формулировка локальной задачи газовой динамики на ячейке периодичности и на ее основе

Рис. 6. Распределение скорости = УУ в поре с застойной зоной для случая нерегулярной КЭ сетки с применением метода QMR и разным числом КЭ:

а — 576 КЭ, б — 751 КЭ, в — адаптивная сетка с числом КЭ 1837

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

СПИСОК ЛИТЕРАТУРЫ

1. Санчес-Паленсия Э. Теория колебаний и неоднородные среды. - М.: Мир, 1984. -472 с.

2. Бахвалов Н. С., Панасенко Г. П. Осреднение процессов в периодических средах. - М.: Наука, 1984. - 352 c.

3. П о б е д р я Б. Е. Механика композиционных материалов. - М.: Изд-во МГУ, 1984. -336 с.

4. Димитриенко Ю. И., Кашкаров А. И. Конечно-элементный метод расчета эффективных характеристик композитов с периодической структурой // Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". - 2002. - № 2. -С. 95-108.

5. Dimitrienko Y u. I. Heat-mass-transfer and thermal stresses in porous charring materials// Transport in Porous Media. - 1997. - V. 27(2). - P. 143-170.

6. Д и м и т р и е н к о Ю. И. Механика композиционных материалов при высоких температурах. - М.: Машиностроение, 1997. - 356 с.

7. Dimitrienko Y u. I. Dynamic transport phenomena in porous polymer materials under impulse termal effects // Transport in Porous Media. - 1999. - V. 35. - P. 299-326.

8. Димитриенко Ю. И., Глазиков М. Л. Численный расчет проницаемости и процессов фильтрации в пористых средах // Аэрокосмические технологии. - М.: Изд-во МГТУ им. Н.Э. Баумана. - 2002. - С. 132-137.

9. Димитриенко Ю. И., Глазиков М. Л. Моделирование процессов фильтрации в периодических пористых средах // Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". - 2003. - № 1. - С. 59-71.

10. Димитриенко Ю.И., Глазиков М.Л. Р азработка метода асимптотического осреднения для решения задач газовой динамики в пористых средах // Математика в современном мире / Под ред. Ю.А. Дробышева. - Калуга.: Изд-во КГПУ - 2004. - С. 163-177.

11. Димитриенко Ю. И. Тензорное исчисление. - М.: Высш. шк., 2001. -575 c.

12. Z i e n k i e v i c z O. Z., T a y l o r L. R. The finite element method. - Oxford: Batterworth-Heinemann. - 2000. - 663 p.

13. Д и м и т р и е н к о Ю. И., Ц а р е в С. М., Веретенников А. В. Разработка метода конечных элементов для расчета элементов конструкций из несжимаемых материалов с большими деформациями // Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". - 2006. - № 3. - С. 102-116.

Статья поступила в редакцию 8.10.2007

Юрий Иванович Димитриенко родился в 1962 г., окончил в 1984г. МГУ им. М.В.Ломоносова. Д-р физ.-мат. наук, профессор, заведующий кафедрой "Вычислительная математика и математическая физика" МГТУ им. Н.Э.Баумана, действительный член Академии инженерных наук. Автор более 140 научных работ в области вычислительной механики, нелинейного тензорного анализа, термомеханики композитов, математического моделирования в материаловедении.

Yu.I. Dimitrienko (b.1962) graduated from the Lomonosov Moscow State University in 1984. D. Sc. (Phys.-Math.), professor, head of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Full member of the Russian Academy of Engineering Sciences. Author of more than 140 publications in the field of computational mechanics, nonlinear tensor analysis, thermomechanics of composite materials, mathematical simulation in science of materials.

Александра Игоревна Левина родилась в 1980 г., окончила МГТУ им. Н.Э. Баумана в 2004 г. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им.Н.Э.Баумана. Специализируется в области математического моделирования процессов газовой динамики.

A.I. Levina (b. 1980) graduated from the Bauman Moscow State Technical University in 2004. Post-graduate of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Specializes in the field of mathematical simulation of processes of gas dynamics.

Павел Боженик родился в 1985 г., окончил МГТУ им. Н.Э. Баумана в 2007 г. Инженер ОАО "НПО Машиностроения".

P. Bozhenik (b. 1985) graduated from the Bauman Moscow State Technical University in 2007. Engineer of open stock company NPO Mashinostroenie.

Вниманию авторов журнала "Вестник МГТУ им. Н.Э. Баумана"!

В связи с изменением условий распространения журнала "Вестник МГТУ им. Н.Э. Баумана" журналы можно приобрести в киоске издательства за наличный расчет. Иногородние авторы могут оформить заявку на приобретение журнала в отделе распространения Издательства (тел. (499) 263-60-45; факс.(495) 261-45-97; E-mail: press@bmstu.ru).

Ориентировочная стоимость одного экземпляра журнала для авторов статей — 150 руб.

Редакция журнала "Вестник МГТУ им. Н.Э. Баумана"