Full Text

Введение Пористые материалы моделируются упругим скелетом, или матрицей, с неким распределением закрытых и открытых пор, содержащих жидкий или газообразный наполнитель. Примерами пористых материалов могут служить водонасыщенные грунты, горные породы, биологические ткани, пенометаллы и т.д. Механика пористых сред имеет большое значение во многих областях науки и техники, таких как геотехника, геомеханика, инженерная геология, биомеханика, машиностроение и материаловедение. В сейсмостойком строительстве важную роль играет прогнозирование потери прочности грунта в результате скачка порового давления жидкости, вызванного сейсмическим воздействием. Началом исследований волновых процессов в насыщенных пористых средах послужила работа Я.И. Френкеля [1]. Л.Я. Косачевский [2] показал, что теория М. Био [3,4], как и подход Я.И. Френкеля, опирается на те же соотношения между напряжениями и деформациями, но отличается большей общностью. Естественным продолжением этих работ стали попытки различных авторов обобщить теорию на случай пористых сред, насыщенных двумя несмешивающимися текучими наполнителями. Принципиальным отличием в данном случае является необходимость учитывать поверхностные натяжения и капиллярные эффекты, возникающие в результате взаимодействия наполнителей. Brutsaert [5] первым расширил теорию Био на случай трехфазной среды, однако исключив из рассмотрения капиллярные эффекты. Общая система уравнений многофазной пористой среды была предложена Николаевским [6]. Работы Garg and Nayfeh [7], Santos и др. [8, 9], Tuncay и Corapcioglu [10, 11], Wei и Muraleetharan [12], Lo и др. [13], Lu и др. [14] по теории пористых сред с двумя жидкостями принадлежат концептуально к общему подходу, отличаясь в деталях реализации. Метод в целом состоял в разработке соотношений между напряжениями и деформациями для составляющих фаз с последующим использованием их в уравнении баланса в рамках стандартной механики сплошных сред. Berryman [15] вывел уравнения для частично насыщенных сред с использованием предположения о незначительном изменении капиллярного давления при прохождении через них акустического сигнала. Smeulders и др. [16] изучали распространение волн сжатия в пористой среде, заполненной жидкостью с небольшой объемной долей газа. В работах Thomas [17, 18] предложена конечно-элементная модель для анализа волновых полей напряжений, перемещений и порового давления в ненасыщенном морском дне. Мейснер и Беккер [19] установили влияние капиллярных эффектов на динамическое поведение частично насыщенного песка. Обзор доступных аналитических решений задач о распространении волн в насыщенных пористых средах можно найти в статье Schanz [20]. Однако для ненасыщенных моделей известно немного существующих аналитических решений. Можно отметить решение, полученное в работе Li и Schanz [21], для одномерного частично насыщенного пороупругого стержня. Поскольку получение аналитического решения часто возможно лишь в некоторых специальных случаях и при определенных видах граничных условий, решение задач динамики частично насыщенных пористых сред требует привлечения развитых численных методов, таких как МКЭ и МГЭ. В случаях бесконечных и полубесконечных тел и сред метод МГЭ является наиболее предпочтительным универсальным подходом к исследованию волн. Для рассмотрения задачи о распространении волны необходимы динамические формулировки исходной системы дифференциальных уравнений и соответствующие ей ГИУ, фундаментальные и сингулярные решения. В работах Gatmiri и Jabbari [22, 23] получены фундаментальные решения для статических и квазистатических задач теории ненасыщенных грунтов в двух- и трехмерной постановке. Соответствующие термо-упругие фундаментальные решения двух и трехмерных статических и квазистатических задач для частично насыщенных сред представлены в работах Jabbari и Gatmiri [24], Gatmiri и др. [25], Maghoul и др. [26]. Динамические формулировки ГИУ и фундаментальные решения для двухмерных задач представлены Maghoul и др. в [27], для трехмерных Ashayeri и др. в [28], Li в [29]. Авторами рассматривается трехкомпонентная пороупругая среда Био с пятью базовыми функциями описания волнового процесса - перемещениями скелета и поровым давлением жидкости и газа. ГИУ получены на основе метода взвешенных невязок, фундаментальные решения представлены в терминах преобразования Лапласа по переменной времени, при построении шаговой гранично-элементной схемы используется метод квадратур свёрток. В настоящей работе представлена гранично-элементная методика на основе шагового метода численного обращения интегрального преобразования Лапласа для решения начально-краевых трехмерных задач динамики частично насыщенных пороупругих тел. Шаговый метод численного обращения преобразования Лапласа опирается на квадратурные формулы для интеграла, получаемого из теоремы операционного исчисления об интегрировании оригинала. Методика исследований основана на граничных интегральных уравнениях прямого подхода трехмерной пороупругости в преобразованиях по Лапласу, на компьютерном моделировании искомых решений методом граничных элементов в сочетании с методом коллокации, локальной поэлементной аппроксимацией на основе согласованной модели интерполирования Гольдштейна. Математическая модель Рассматривается частично насыщенный пористый материал, который представляет собой упругий скелет с порами, заполненными жидкостью и газом, при этом исключается испарение жидкости и растворение газа [29]. Кроме того, предполагается, что все три фазы являются сжимаемыми, имеют постоянную температуру, а любые изменения температуры не учитываются. Для частично насыщенного пористого материала пористость обозначается как , (1) где - объем взаимосвязанных пор в образце объема . Закрытые поры рассматриваются как часть упругого скелета. Степень насыщения материала жидкостью и газом определяется как отношение занимаемого наполнителем объема и соответственно к общему объему пор, т.е. . (2) Существенным отличием частично насыщенных пористых сред от частично ненасыщенных является необходимость учитывать взаимодействие наполнителей. В результате искривления поверхности жидкости, вызванного действием сил поверхностного натяжения на границе раздела сред, возникает разность давлений в жидкости и газе , которая называется капиллярным давлением . Капиллярное давление может быть представлено как функция степени насыщения . (3) Величина определяет давление газа, необходимое для вытеснения жидкости из пор. Коэффициент , характеризующий распределение размера зерен скелета, обычно принимает значения от 0,2 до 3. - эффективная водонасыщенность определяется следующим образом [29]: (4) где - остаточная водонасыщенность; - остаточная газонасыщенность. Определяющие соотношения для суммарного напряжения вслед за [30] могут быть записаны в виде (5) Здесь обозначено эффективное напряжение и - коэффициент эффективных напряжений, определяющий зависимость сжимаемости упругого скелета от его объемного модуля и объемного модуля его зерен . Скелет считается изотропным и однородным с двумя константами упругого материала: объемный модуль и модуль сдвига . Исходя из этих предположений получена формулировка для напряжения: (6) Кроме того, в допущении о малых градиентах деформации полагаем, что компоненты тензора деформаций твердого тела и смещения связаны соотношением (7) Уравнения баланса массы для твердого тела и наполнителей имеют вид , (8) , (9) (10) где , , - плотности твердотельного скелета, жидкости и газа соответственно, и - компоненты смещений наполнителей, относительных по отношению к скелету. Формулировка предполагает, что растворения газообразного наполнителя в жидкости не происходит, а градиенты плотностей, насыщенности и пористости пренебрежимо малы. Для трехфазного материала можно записать , (11) , (12) . (13) Уравнение (11) записано для скелета, уравнения (12), (13) - для жидкости и газа соответственно. Введена объемная сила и усредненная плотность . Перенос наполнителей определяется законом Дарси. Соответствующие фазовые проницаемости определяются как и , где и - относительные фазовые проницаемости; - внутренняя проницаемость материала; и - вязкости наполнителей. Относительная фазовая проницаемость определяется как отношение фазовой проницаемости к абсолютной. Используем вслед за [31] следующие экспериментально установленные зависимости для относительных проницаемостей в случае газоводяной смеси в порах: , . (14) Система уравнений (8)-(13) полностью описывает поведение частично насыщенного пороупругого континуума. Выберем в качестве базовых функций давление в порах , и перемещение скелета . Для дальнейшего рассмотрения применяется интегральное преобразование Лапласа. Пусть - функция по , = 0 для . Функция преобразования по Лапласу и ее обращение определено следующим образом: , (15) (формула Римана-Меллина), (16) где s - комплексный параметр преобразования; a - вещественное число большее, чем вещественные части всех особенностей . Применяя преобразование Лапласа к уравнениям (12), (13) и исключая относительные смещения и , получаем , (17) где , - потоки наполнителей. В (17) введены следующие обозначения: . (18) Подстановка (17) в уравнения (8)-(11) позволяет получить формулировку , (19) , (20) , (21) где , , , . Уравнения (19)-(21) могут быть записаны в операторной форме: (22) с несамосопряженным оператором , , , , , , , , , , . Применение интегрального преобразования позволяет строить метод граничных интегральных уравнений. 2. Граничные интегральные уравнения и фундаментальные решения Фундаментальное решение находится из следующего операторного соотношения: , (23) где - единичная матрица; - дельта-функция Дирака; . Итоговая матрица фундаментальных решений (24) содержит девять компонент - компоненты фундаментальных решений перемещений и давлений , , а также шесть спаренных компонент. Все компоненты матрицы приведены в [29]. Гранично-интегральное уравнение для динамической пороупругости в области Лапласа можно получить на основе традиционных формулировок: теоремы взаимности или метода взвешенных невязок. Пороупругий оператор является несамосопряженным, а значит, для этих методов требуются разные фундаментальные решения, но в конечном счете оба метода дают одно и то же интегральное уравнение. Интегральное уравнение динамической пороупругости, построенное на основе метода взвешенных невязок, имеет вид , (25) где - матрица сингулярных решений; - поверхностная сила; - поток жидкости; - поток газа; - вектор внешней нормали. Построенные граничные интегральные уравнения являются сингулярными и существуют в смысле Коши. Для нахождения решения граничных интегральных уравнений требуется соответствующая методика их дискретного представления. В работе используется гранично-элементный подход. 3. Гранично-элементная (ГЭ) дискретизация Чтобы ввести ГЭ-дискретизацию, рассмотрим регуляризованное уравнение [32]: где - матрицы фундаментальных и сингулярных решений; - матрица компонент, содержащих особенности; - вектор обобщенных усилий; - вектор обобщенных перемещений. Базовый процесс ГЭ-дискретизации состоит в разбиении поверхности на граничных элементов совокупностью четырехугольных и треугольных восьмиузловых биквадратичных элементов. При этом треугольные элементы рассматриваются как вырожденные четырехугольные элементы, каждый из которых отображается на контрольный элемент (каждый - это либо квадрат , либо треугольник ). Неизвестные граничные поля определяются через узловые значения и в интерполяционных узлах. При этом для расчетного значения параметра s будем иметь следующие выражения граничных перемещений и поверхностных сил внутри элемента Sk: Здесь - функции формы для линейного четырехугольного элемента. Для получения дискретного аналога ГИУ применим метод коллокации. В качестве узлов коллокации будем выбирать узлы аппроксимации исходных граничных функций. В итоге формируется система линейных алгебраических уравнений . Уравнения записаны соответственно в узлах аппроксимации обобщенных граничных перемещений и обобщенных поверхностных сил. , 4. Шаговая схема численного обращения преобразования Лапласа Рассмотрим метод, опирающийся на теорему об интегрировании оригинала, - шаговый метод численного обращения преобразования Лапласа. Рассмотрим следующий интеграл: . (26) Интеграл (26) порождает задачу Коши для обыкновенного дифференциального уравнения Интеграл (26) заменяется квадратурной суммой, весовые множители которой определяются с помощью изображения по Лапласу и линейного многошагового метода [33-36]. Дальнейшее изложение идет с учетом результатов этих работ. Традиционный шаговый метод интегрирования оригинала состоит в том, что интеграл (26) считается по следующему соотношению: (27) где . Единственная аппроксимация, которая используется при выводе формулы (27), заключается в применении линейного многошагового метода для решения задачи Коши дифференциального уравнения первого порядка. Все остальные вычисления проводятся непосредственно. Многошаговый метод должен быть порядка точности , являясь строго ноль-устойчивым или A-устойчивым, а функция ограничена в правой полуплоскости относительно прямой , то есть при Если функция аналитична и ограничена в области , где , критерий устойчивости может быть ослаблен до A(a)-устойчивости. Выражаясь в терминах характеристической функции : функция не должна содержать нулей и полюсов внутри замкнутого единичного круга за исключением однократного нуля при где при , где при К соответствующим примерам многошаговых методов относятся методы дифференцирования назад порядка : для A-устойчивого метода дифференцирования назад второго порядка можем записать следующую характеристическую функцию: Следующая аппроксимация, которая используется в методе, - численное интегрирование, при получении значений аппроксимируем с равными шагами : . (28) При условии, что функция в уравнении (28) вычисляется с некоторой погрешностью выбор и допускает погрешность вычисления порядка . 5. Численные эксперименты Рассматривается задача о действии торцевой силы на призматическое пороупругое тело длиной 3 м, высотой и шириной 1 м, с жестко закрепленным концом. Краевая задача представлена на рис. 1. Отклики перемещения и поверхностные усилия , вызванные силой = , наблюдаются соответственно на нагруженном и на закрепленном торцах. Параметры частично насыщенного пористого материала соответствуют параметрам песчаника: = , = , = , = , = , = , = , = , = , = , = , = , = 0,9, = 0, = 1, = 1,5. На временном интервале от 0 с до 0,035 с шаг по времени выбирается из условия L = N, параметр R принят равным 0,997. Задача решается в трёхмерной постановке. Заданные граничные условия позволяют сравнивать решение исходной задачи с решением соответствующей одномерной задачи, представленным в [21]. При решении использовались три гранично-элементные сетки с различной степенью пространственной дискретизации (рис. 2). Наличие двух плоскостей симметрии дает возможность рассматривать только часть сетки, при этом четверть сетки a содержит 14 элементов, сетки б - 56 элементов, сетки в - 126 элементов. Рис. 1. Геометрическая модель и граничные условия задачи о частично насыщенной пороупругой консоли Fig. 1. Geometry and boundary conditions of a partially saturated poroelastic column а б в Рис. 2. Гранично-элементные сетки: a - 56 элементов; б - 224 элемента; в - 504 элемента Fig. 2. Boundary-element meshes: a - 56 elements; b - 224 elements; c - 504 elements Рис. 3. Гранично-элементные решения, полученные на трех различных сетках: a - перемещения; б - усилия Fig. 3. Boundary-element solutions on the three different meshes for: a - displacements; b - tractions Рис. 4. Гранично-элементные решения, полученные при трех различных значениях параметра шагового метода: a - перемещения; б - усилия Fig. 4. Boundary-element solutions for the three different values of the time-step method parameter for: a - displacements; b - tractions Результаты исследования сходимости по расчетной сетке показывают, что с ростом числа элементов гранично-элементные решения для перемещения (рис. 3, а) и напряжения (рис. 3, б) принимают значения, близкие к рассчитанным по аналитическим формулам. Варьирование параметром шаговой схемы приводит к тому, что при снижении его до значения 0,00002185c в расчетах на сетке c наблюдается незначительное уточнение амплитуды решения (рис. 4, а, б), однако при дальнейшем его уменьшении решение расходится. Заключение Дано описание математической модели частично насыщенной пористой среды, основанной на модели Био, с пятью базовыми функциями. В изображениях по Лапласу приведена математическая модель краевой задачи трехмерной динамической теории пороупругости. Представлен шаговый метод численного обращения преобразования Лапласа. Дано описание методики гранично-элементного моделирования. Решена задача о действии нагрузки на торец призматического частично насыщенного пороупругого тела. Полученные гранично-элементные и рассчитанные аналитически значения искомых величин графически неразличимы.

About the authors

L A Igumnov

Research Institute of Mechanics, National Research Lobachevsky State University of Nizhny Novgorod

Email: igumnov@mech.unn.ru

A N Petrov

Don State Technical University

Email: andrey.petrov@mech.unn.ru

References

Statistics

Views

Abstract - 301

PDF (Russian) - 82