CC BY
12
УДК 544.03+519.63
МЕТОД ВЕЙВЛЕТ-ОСРЕДНЕНИЯ ДЛЯ ПРОГНОЗИРОВАНИЯ ЭФФЕКТИВНОГО КОЭФФИЦИЕНТА ПРОНИЦАЕМОСТИ
КОПЫСОВ СП., САГДЕЕВА Ю.А.
Институт прикладной механики УрО РАН, 426067, г. Ижевск, ул. Т.Барамзиной, 34
АННОТАЦИЯ. Предлагается метод прогнозирования эффективных характеристик и анализа осредненных решений уравнений для неоднородных сред с известными структурой и свойствами составляющих компонент на основе вейвлет-преобразования и метода конечных элементов для задачи проницаемости.
КЛЮЧЕВЫЕ СЛОВА: эффективные характеристики, вейвлет-осреднение, коэффициент проницаемости, метод конечных элементов.
ВВЕДЕНИЕ
Для многих физических процессов в неоднородных средах (теплопроводность, электропроводность, фильтрация жидкостей и газов) характерна математически эквивалентная задача определения макроскопической или, так называемой, эффективной проводимости системы на основании информации о структуре локального поля. Определение эффективных параметров различных полей - задачи математически полностью или частично эквивалентные, многие результаты универсально могут использоваться при получении характеристик разной физической природы. Наиболее интересен и практически важен стохастический вариант этой задачи, то есть тот случай, когда локальное поле может трактоваться как случайное.
Рассмотрим двумерную задачу фильтрации жидкости, описываемую дифференциальным уравнением
-Уг (К (х, у)Уи ) + ( = 0 (1)
с быстро осциллирующим коэффициентом проницаемости К(х,у). Непосредственное численное решение уравнений вида (1) (например, методом конечных разностей или методом конечных элементов) требует значительных вычислительных затрат, поскольку предполагает использование расчетной сетки очень малого шага. Поэтому проводят осреднение соответствующего дифференциального оператора в (1).
Целью работы является моделирование осредненных коэффициентов уравнения (1) так, чтобы решения осредненного уравнения (1) были близки к решениям исходных уравнений; такие коэффициенты называются эффективными. Рассматриваются два аспекта -получение осредненных решений уравнения (1) и вычисление эффективного коэффициента К. Эффективный коэффициент проницаемости в условиях выполнения закона Дарси определялся на основе вейвлет-осреднения и решения методом конечных элементов уравнения (1).
В настоящее время существует несколько подходов для описания и анализа эффективных характеристик неоднородных сред. Такими подходами являются, например, метод асимптотического осреднения - этот метод фактически можно отнести к двухмасштабным методам [1], [2], [3] и др., различные микромеханические и смесевые модели (оценки Хашина-Штрикмана, Фойгта-Рейсса [4]). Эти методы не всегда приемлемы в силу существования ограничений по их применению (периодичность среды, двухфазность и некоторые другие). Например, следует иметь в виду, что существующие в настоящее время многочисленные "вилки", правила смесей и инженерные расчетные формулы дают зачастую очень грубое приближение к эффективным характеристикам.
В данной работе предлагается метод осреднения, основанный на многомасштабном анализе с вейвлет-проекцией и аппроксимацией дискретного оператора. В отличии от других
МЕТОД ВЕИВЛЕТ-ОСРЕДНЕНИЯ ДЛЯ ПРОГНОЗИРОВАНИЯ ЭФФЕКТИВНОГО _КОЭФФИЦИЕНТА ПРОНИЦАЕМОСТИ_
методов осреднения, например [1], где учитываются только эффекты самого малого (локального) масштаба и рассматривается случай только периодической структуры материала, в данной работе исследуется многомасштабное вейвлет-осреднение решений уравнений (1) для неоднородных сред.
МНОГОМАСШТАБНЫЙ АНАЛИЗ
Рассмотрим пространство £2(М). Будем называть вейвлетами функции, образующие базис пространства £2 (М) и получаемые сдвигом и сжатием функции у( х) ¥] к (х) = 2; 2 у(2; х - к), у, п е Ъ , где Ъ - пространство целых чисел, а масштабирующими -функции, получаемые сдвигом и сжатием функции <(х) : <р.к (х) = 2; /2<(2;х - к) . В качестве
вейвлет-базиса в работе был выбран базис Хаара, поскольку он обладает свойствами ортогональности и симметрии:
| 1, х е[0,1/ 2) |1, х е [0,1]
[-1, х е [1/2ДУ (О, х й [0,1].
Иерархические свойства вейвлетов лежат в основе многомасштабной декомпозиции функции. Многомасштабный анализ пространства (М) есть последовательность пространств V = к (х)}, имеющих следующие свойства: вложенность
... ^ ^ V ^ ... ^ V ^ ... ^ (М) , и инвариантность относительно масштаба и сдвига (иерархичность) /(х) еК, о /(2х - к) е . Благодаря этим свойствам, пространство определяется как прямая сумма пространства V и дополнения к нему = V Ф W;., W;. = s^an{у;k (х)}, у е 2 . Произвольная функция /(х) е £2(М) аппроксимируется последовательностью функций / (х) е V :
/ (х) = Е(/ (x), <,к (х))<,к (х) = ЕсМ (х)<,к (х) = Есу-1,к<-1,к (х) + 2Й'-1,к¥у-1,к(х).
Функция из V представляется в виде совокупности грубого приближения (первая сумма) и уточняющих добавок (вторая сумма).
Пусть вектор с;. имеет размер 2;. Этот вектор представим как последовательность коэффициентов разложения в многомасштабном анализе некоторой функции / :
2 '
Л = Ёсм (х)<,к (хХ су = у Су ,2 у).
к=1
Коэффициенты разложения в пространствах V, и W; связаны между собой: су-1 = у, = . Вектор су-1 является проекцией вектора су на пространство , то есть его огрубленным представлением. Вектор соответствует уточняющим
коэффициентам. Матрица W; = (р Qj У является матрицей вейвлет-преобразования. Вейвлет-преобразование Хаара является ортогональным, то есть W;-1 = W;г.
Рассмотрим ограниченный линейный оператор : V ^ V. Введем
многомасштабный анализ на некоторой ограниченной области, тогда оператору будет соответствовать матрица конечной длины. Опишем общую процедуру осреднения операторного уравнения £;х = / . Разложение V = Ф позволяет разделить оператор
на четыре части и записать
(а8 в, у d^ (df л
х = f СОЛ
С т - . (2)
С с т <
V 81 А 'х ) V" г )
где А5] : Ж- ^ Ж-, В8} : У-_х ^ , С8] : Ж-_х ^ , Т^ : ^ V-, а dx, df е , 'х, ¿у е V— - ортогональные проекции х и f на подпространства Ж-1 и У;-1. Проекция 'х является компонентой грубого масштаба решения х, dx - мелкомасштабной компонентой. Формально исключая dx из (2) подстановкой dx = A-l (df - В8 'х), получим
(Т5] -С8л-в81 к = -с8А-йг. (3)
Оператор Я8 = Т8 - С8 А-1 В8 является дополнением Шура.
Решение 'х уравнения (3) равно в точности Рх, где Р}-1 - оператор проектирования на У-1, ах - решение системы. Заметим, что данное осредненное уравнение отличается от уравнения, задаваемого оператором проекции Т8 ,
Т8^х = Sf ,
то есть при осреднении учитывается информация о деталях, тогда как при проектировании она теряется.
ВЕЙВЛЕТ-ОСРЕДНЕНИЕ НА ОСНОВЕ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ
Дискретизируем уравнение (1) методом конечных элементов на сетке с помощью линейных элементов. Предполагается, что заданы краевые условия Дирихле или Неймана. В результате применения МКЭ получим систему линейных алгебраических уравнений
Аи = Ь, А = (а1к), Ь = (Ь), и = (иг), ¡, к = 1... 2;, (4)
где матрица А - симметрична и положительно определена.
Введем расширение вейвлет-представления с одномерного случая на двумерный в пространстве Ь2 (М). Базисные функции образуются за счет всевозможных комбинаций тензорных произведений базисных функций для одномерного случая. В этом случае получаются четыре порождающих функций - одна масштабирующая функция Ф (х,у) = ф(х)ф(у) и три вейвлета у/{1 }(х,у) = ф(х)щ(у), у/{11 )(х,у) = ¥(х)ф(у), у/(ш)(х, у) = у/(х)щ(у). Все остальные функции определяются соотношениями
^(х, у) = 2>(/)(2; х - к,2; у - /), у) = 2>(//)(2^ х - к,2 у - /),
^»(х, у) = 2; У 1 )(2;х - к2 у - /), Аналогично одномерному случаю, введем последовательность вложенных пространств Vп и пространства №п:
Vп = 'Рап{Фп,к1 ®Фп,к^ к Е ZК
№п = храп&п,к3 ® Фп,к4, Фп,к5 ® ¥пм6, ¥пМ1 ® ¥пк8, к е 1}.
Как видно dim(Wп) = 3 х dim(Vп ).Пространство Vп = Уп ® Уп отвечает за осредненные
величины, пространство №п = (Жп ® Жп) Ф (Жп ® Уп) Ф (Уп ® Жп) содержит информацию о
взаимосвязи двух направлений и уточняющую информацию.
В двумерном случае вейвлет-преобразование формируется из четырех операторов проекций
( р л
№ =
п
£ 2
^ 3)
МЕТОД ВЕИВЛЕТ-ОСРЕДНЕНИЯ ДЛЯ ПРОГНОЗИРОВАНИЯ ЭФФЕКТИВНОГО _КОЭФФИЦИЕНТА ПРОНИЦАЕМОСТИ_
Р: V, ®V, ^ ®К^, б! : V, ®V, ^ ®Ж^,
б2 :V, ®V, ^Жп-, ®Кп-,, бз:V, ®V, ^Ж,- ® Операторы Р,б,,б2,б3 состоят из двух блоков - один действующий по координате X , второй по координате у : Р = РхРу , б 1 = , б 2 = , бз = бхб^ •
Преобразование ^; является ортогональным, так как. свойства базисных функций сохраняются. Применяя преобразование ^; к системе (5), получим (индексы ] и ] -1 опущены)
Ш = Ш,
Г р ^ ( - Л м Г Ъ-^
б ! б 2 А( Р гб Гб Гб Г) й й м 2 = Ъй Ъ2й
1б з У Vм з У VЪй У
Вводя следующие обозначения
б, А С б,
м„ =
б, А бз б 2 Аб3
гЛ
б з АбГ б з Аб2 б з А б г
М12 =
у
Г б, АР[ ^ б 2 АР Г б з АР Г
Г м-^
м, =
Ъ =
V з У
Г Ъ,й ^
V з у
м2 = мс, Ъ2 = Ъс и с учетом того, что М21 = М1г2 перепишем систему в виде
ГМ,2 Ум, ^ Гъ, ^
Vм 21
112
М22 УуМ2 У
V Ъ2 У
(5)
Размер полученной системы совпадает с размером исходной системы. Выразив м, из первого уравнения и подставив его во второе, получим систему
£м2 = 8, (6) где £ - дополнение Шура: £ = М22 - М^М./М^, 8 = Ъ2 - М^М^Ъ,. Разрешив (6) получаем искомое осредненное решение м2.
Данный подход близок к процедуре исключения неизвестных в блочном методе Гаусса. Отличие состоит в том, что процедура вейвлет-преобразования требует изменение базиса, которое выполняется каждый раз перед шагом редукции. Таким образом, неизвестные в приведенной системе не являются простым подмножеством неизвестных исходной системы. Матрица £ системы (6) получается из конечно-элементной матрицы жесткости и описывает осредненные свойства материала, но непосредственно использовать эту матрицу в качестве грубого оператора на сетке не удается. Алгоритм вейвлет-осреднения, примененный к системе (4), позволяет получить осредненные поля решения дифференциального уравнения (1).
МЕТОД ВЫЧИСЛЕНИЯ ЭФФЕКТИВНЫХ ХАРАКТЕРИСТИК
Для вычисления осредненных характеристик используем осредненное решение системы (6) по следующему алгоритму. Опишем его на примере определения свойств проницаемости среды.
Эффективный коэффициент проницаемости (двумерный случай). Эффективный коэффициент проницаемости в условиях выполнения закона Дарси определялся на основе вейвлет-осреднения и решения методом конечных элементов двумерной задачи фильтрации жидкости, описываемой уравнением
-Уг (К (х, у)Ур) = 0, (7)
в насыщенной пористои неоднородной среде со следующими граничными условиями: на одной границе области задано условие равенства нулю давления р | х=0 = 0, на
противоположной границе задана компонента
К (х. У) £
ох
х
вектора скорости фильтрации
= у0 , на оставшейся границе задано условие непроницаемости границы (поток
х=1
= 0.
у=0, У=1
через границу равен нулю) —
ду
Применим к задаче (7) метод конечных элементов, а к системе, полученной МКЭ,
применим двумерное вейвлет-преобразование Хаара. Обозначим через р*" - осредненное с
помощью вейвлет-преобразования поле давления в узле сетки (xi, у). Эффективный
коэффициент проницаемости рассчитывался как среднее значение по всем узлам ( N - число узлов)
= 1У У кх, кх =- ^.
По аналогии с коэффициентом проницаемости можно определить эффективный коэффициенты теплопроводности и упругие константы неоднородной среды.
РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
Задача 1. Рассмотрим пример осреднения тензора проницаемости в двумерной задаче фильтрации (7) для случайного начального распределения проницаемости с граничными условиями
Н=0 = 0,
К (х, У) %
ох (8)
=40, °Р
ду
=0.
У=0, ,у=1
Тензор проницаемости имеет вид К(х) = d(х)
= 1п(а)
Г1 0^ 0 1
, где величина d(х) имеет
V" V
следующий закон распределения: d(х) = г- 1п(а), г - случайно распределенная величина на интервале (0,1) и величина d (х) предполагается кусочно-постоянной на конечном элементе. Пример данного распределения при а = 2 показан на рис. 1.
х=1
Рис. 1. Случайное поле проницаемости
уп
У
Рис.2. Исследуемые формы сечения волокна
МЕТОД ВЕЙВЛЕТ-ОСРЕДНЕНИЯ ДЛЯ ПРОГНОЗИРОВАНИЯ ЭФФЕКТИВНОГО _КОЭФФИЦИЕНТА ПРОНИЦАЕМОСТИ_
Для распределения такого вида известно (см. [5]), что значение эффективного коэффициента проницаемости определяется как геометрическое среднее значений распределения. Результаты моделирования для значений а = 2, а = 5 и а = 10 на сетке 128 х 128 узлов приведены в табл. 1. В первом-третьем столбцах даны статистические данные о начальном распределении проницаемости - максимальное значение Ктах, минимальное
значение Ктт, среднеквадратическое отклонение а. В четвертом-седьмом столбцах показаны соответственно арифметическое среднее Ка, гармоническое среднее К, геометрическое среднее К8 и значение К№, полученное с помощью вейвлет-осреднения.
Таблица 1
Сравнение эффективного коэффициента проницаемости для случайного распределения
а К тах К тт а Ка Кн К8 К№
2 255,87 1,00 4,45 2,57 1,6 2 2,07
5 6,06е+6 1,00 67936,3 992,47 2,59 4,95 5,13
10 2,5е+10 1,00 2,01е+8 2,21е+6 3,29 9,86 9,5
Очевидно, результаты моделирования близки к геометрическим средним -
I К8 — К№ I
- равны 3,5 %, 3,6 %, 3,7 % для разных а.
относительная погрешность
К8
Арифметическое и гармоническое средние дают настолько широкий интервал оценок, что их использование не имеет смысла. Заметим, что аналитические оценки Хашина-Штрикмана [4] для эффективных свойств материала, а также асимптотический метод осреднения [1] не применимы в данном случае. Другие примеры использования вейвлет-преобразования можно найти в [6].
Задача 2. Рассматривается задача фильтрации (7) с граничными условиями (8). Предполагается, что исследуемая среда состоит из двух изотропных сред (назовем их матрицей и включением), тензор проницаемости которых имеет вид (см. рис. 2)
Г10^,, Г1 0^
К (х) = а
V0 1J
Ух е У(т), К(х) = а
01
Ух е У
(V)
Объемная доля включения составляет см = 0,25 .
Таблица 2
Сравнение аппроксимаций эффективного коэффициента проницаемости
а К№ Ка MG-осреднение Ренорма-лизация К? - К? К? - К?
0,01 0,0168 0,017 0,0180 0,0159 0,013-0,258 0,0164-0,1526
0,1 0,1522 0,152 0,16 0,147 0,129-0,325 0,1514-0,2394
1 1 1 1 1 1 1-1
10 6,545 6,55 6,70 6,17 3,08-7,75 4,176-6,604
100 59,22 59,20 61,74 54,25 3,88-75,25 6,55-60,64
1000 585,31 579,32 611,74 534,25 3,99-750,25 6,95-600,64
Результаты вейвлет-осреднения для коэффициента проницаемости для разных значений константы а сведены в табл. 2. Анализ результатов проводился для значений проницаемости материала У (т) (первая колонка) как больших, так и меньших проницаемости включения. Второй столбец таблицы соответствует результатам вейвлет-осреднения. Результаты из столбцов три-пять взяты из работы [7] и соответствуют различным методам
осреднения - асимптотическому методу, методу MG-осреднения и ренормализации. В столбце шесть приведены оценки Фойгта-Рейсса Kfr, Kf (эти оценки являются
соответственно арифметическим и гармоническим средними значениями), которые задают широкий диапазон допустимых значений. В последнем столбце приведены оценки Хашина-Штрикмана (K - Khus) для трансверсально изотропного двухфазного материала. И хотя рассматриваемая структура материала не является трансверсально изотропной, в силу симметрии задачи и структуры начального коэффициента K отличие в значении осредненного коэффициента по направлению оси, сонаправленной со стороной квадрата, от значения в случае, когда свойства исследуются вдоль диагонали квадрата, будет не очень большим. Поэтому оценки Хашина-Штрикмана также полезны при сравнении. Вейвлет-осреднение и асимптотическое осреднение удовлетворяют оценкам Хашина-Штрикмана для всех а .
ЗАКЛЮЧЕНИЕ
В работе проведено построение и исследование метода осреднения для получения эффективного коэффициента проницаемости неоднородной среды, основанного на вейвлет-осреднении дифференциального оператора. Приведены примеры, показывающие результаты вейвлет-осреднения для сред регулярной и случайной структуры. Достоверность полученных результатов подтверждается сравнением с известными экспериментальными данными и другими методиками получения эффективных коэффициентов.
Работа выполнена при финансовой поддержке РФФИ проект № 10-01-96039-р_урал. СПИСОК ЛИТЕРАТУРЫ
1. Бахвалов Н.С., Панасенко Г.П. Осреднение процессов в периодических средах. М. : Наука, 1984, 352 с.
2. Жиков В.В., Козлов С.М., Олейник О.А. Усреднение дифференциальных операторов. М. : Физматлит, 1993, 464 с.
3. Санчес-Паленсия Э. Неоднородные среды и теория колебаний. М. : Мир, 1984. 472 с.
4. Hashin Z., Shtrikman S. A variational approach to the theory of the elastic behaviour of multiphase materials // Journal of the mechanics and physics of solids. 1963. V.11. P.127-140.
5. Козлов С.М. Осреднение случайных операторов // Математический сборник. 1979. Т.109. С.188-202.
6. Копысов С.П., Сагдеева. Ю.А. Двумерное численное вейвлет-осреднение для получения эффективных характеристик композиционных материалов // Математическое моделирование. 2009. Т.21. С.65-78.
7. Knapek S.. Upscaling techniques based on subspace corrections and coarse-grid approximations // Situ. 1998. V.22, №1. P.35-58.
WAVELET-HOMOGENIZATION METHOD FOR FORECASTING EFFECTIVE PERMEABILITY COEFFICIENT
Kopysov S.P., Sagdeeva Yu.A.
Institute of Applied Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. A new method to forecast effective characteristics and analyze averaged solutions of the partial differential equations for heterogeneous media with given structure and phases is proposed. The method is based on the wavelet transformation and finite element method.
KEYWORDS: wavelet homogenization, effective characteristics, permeability coefficient, finite element method.
Копысов Сергей Петрович, доктор физико-математических наук, старший научный сотрудник, ведущий научный сотрудник ИПМ УрО РАН, тел. (3412)214583, е-таИ: [email protected]
Сагдеева Юлия Альбертовна, кандидат физико-математических наук, старший научный сотрудник ИПМ УрО РАН, тел. (3412)214583, е-таИ: [email protected]
