Numerical analysis of poroviscoelastic prismatic solids and halfspaces dynamics via boundary element method
- Authors: Ipatov AA1, Belov AA1, Litvinchuk SY.1
- Affiliations:
- Research Institute for Mechanics, National Research Lobachevsky State University of Nizhny Novgorod
- Issue: No 4 (2016)
- Pages: 248-262
- Section: ARTICLES
- URL: https://ered.pstu.ru/index.php/mechanics/article/view/206
- DOI: https://doi.org/10.15593/perm.mech/2016.4.14
- Cite item
Abstract
Dynamic behavior of poroelastic and poroviscoelastic solids is considered. Poroviscoelastic formulation is based on Biot’s model of fully saturated poroelastic media. The elastic-viscoelastic correspondence principle is applied to describe viscoelastic properties of elastic skeleton. Viscoelastic constitutive equations are introduced. Classical viscoelastic models are used, such as Kelvin-Voigt, Standard linear solid and model with weakly singular kernel of Abel type. Differential equation system of full Biot’s model in Laplace transform and formulas for elastic modules are given. Original problem’s solution is built in Laplace transform and numerical inversion is used to obtain the solution in time domain. Direct boundary integral equation (BIE) system is introduced. Regularized BIE system is considered. Mixed boundary element discretization is introduced to obtain discrete analogues. Gaussian quadrature and hierarchic integrating algorithm are used for integration over the boundary elements. Numerical inversion of Laplace transform is done by means of modified Durbin’s algorithm with a variable integrating step. The described numerical scheme is verified by a comparison with analytical solution in a one-dimensional case. Isotropic poroviscoelastic solids and halfspaces are considered. Results of numerical experiments are presented. Problems of axial force acting on the end of prismatic solid and vertical force acting on a halfspace are solved. Viscous parameter influence on dynamic responses of displacements and pore pressure are studied. Surface waves on poroviscoelastic halfspace are modelled with the help of boundary element method.
Full Text
Введение В развитии механики пористых материалов заинтересованы специалисты химических, нефтехимических отраслей, а также специалисты по механике грунтов и биомеханике. Для широкого класса насыщенных материалов, таких как водонасыщенные грунты, насыщенные нефтепродуктами скальные породы и другие образования с воздушными порами, классическая теория упругости не может быть применена для корректного их описания. Возникла необходимость развития новой теории, и началом стали работы Френкеля Дж. [1] и M. Biot [2-4]. Модель Био является расширением классической теории упругости на случай двухфазной среды, состоящей из упругого скелета и жидкого или газообразного наполнителя. В последующие годы теория получила распространение для анизотропного случая, а также для динамики. Из современных работ по данной тематике стоит отметить работы M. Schanz [5], R. de Boer [6] и др. [7-9]. Некоторые аналитические решения для насыщенных пористых материалов в одномерном случае были получены в работах [10-12]. Классическая теория вязкоупругости представлена, к примеру, у Р. Кристенсена [13]. Теория пористых материалов, обладающих вязкоупругими свойствами, впервые была предложена M. Biot [4]. В дальнейшем на эту тему были опубликованы работы по квазистатике R.K. Wilson и E.C. Aifantis [14] и динамике I. Vgenopoulou и D.E. Beskos [15]. Одномерное аналитическое решение для стержня в случае динамической поровязкоупругости представлено M. Schanz и Cheng [16]. При проведении мониторинга приповерхностных зон земной поверхности требуются эффективные средства, методы и модели. В настоящее время преобладают работы по изучению колебаний в пористых средах с применением метода нормальных мод, лучевого метода, а также метода контурных интегралов. Все возрастающие требования к строгости подходов значительно усложняют схемы решения задач. Метод граничных элементов (МГЭ) как универсальный численно-аналитический подход позволяет для достаточно сложных математических моделей получить высокоточные результаты компьютерного моделирования. Преимущество метода проявляется наиболее существенно в случае бесконечных и полубесконечных тел. В настоящее время сформировались два основных направления развития МГЭ, применяемых при исследовании волновых задач: использование интегрального преобразования и построение шаговых процедур. Первая граничная интегральная формулировка для упругодинамики была опубликована T.A. Cruse и F.J. Rizzo [17, 18]. Эта формулировка применялась в сочетании с преобразованием Лапласа. Первой работой по применению ГИУ непосредственно во временной области является работа W.J. Mansur и C.A. Brebbia [19]. Детальный обзор по начальному этапу упругодинамических аспектов граничных элементов и их применению содержится в [20, 21]. Так как в общем случае отсутствуют фундаментальные решения для пороупругости, то в явном времени соответствующие МГЭ-формулировки не могут быть обобщены. Однако возможности методов ГИУ и МГЭ позволяют успешно моделировать динамику пороупругих тел [22-25]. Приближенный МГЭ-подход для динамических задач - это МГЭ с двойным применением теоремы взаимности [26] и метод, опирающийся на регулярные ГИУ первого рода [27-30]. Данная работа посвящена развитию прямого подхода метода граничных элементов для решения нестационарных динамических трехмерных задач изотропной поровязкоупругости со смешанными краевыми условиями. Метод граничных элементов для решения таких задач недостаточно развит. Он позволяет решать динамические задачи для трехмерных для изотропных упругих и вязкоупругих тел достаточно сложных конфигураций, условий нагружения и закрепления. В методе граничных элементов ключевую роль играют фундаментальные и сингулярные решения, выражения которых в случае пороупругости известны только в пространстве преобразований Лапласа. В работе рассматриваются начально-краевые задачи пороупругих и поровязкоупругих изотропных призматических тел, которые решаются с помощью МГЭ с использованием интегрального преобразования Лапласа. Решение в явном времени получается путем применения модификации численного алгоритма обращения преобразования Лапласа, предложенного в работе Ф. Дурбина [31]. 1. Математическая модель Рассмотрим однородное тело в трехмерном евклидовом пространстве с декартовой системой координат , - граница тела. В теории Био рассматривается полностью насыщенный материал. Пористость обозначается как где - объем взаимосвязанных пор в образце объемом V. Закрытые поры рассматриваются как часть твердого тела. Насыщение считается полным, т.е., , где - объем твердого тела. Уравнения движения пороупругой деформируемой среды в области Ω имеют вид Введены объемные силы скелета и наполнителя , соответствующие плотности обозначены как и . - плотность объемной силы. Эти уравнения дополняются физическим соотношением, геометрическими соотношениями и динамическим законом Дарси где - константы упругости; - компоненты тензора деформации; - компоненты тензора напряжения; - компоненты плотностей объемной силы; - вектор перемещения скелета; - вектор перемещения фильтрации (просачивания); - плотность присоединенной массы; k - проницаемость. После формального применения преобразования Лапласа система дифференциальных уравнений теории Био в обобщенных перемещениях в преобразованиях Лапласа (параметр преобразования s) для перемещений и порового давления принимает вид [5] , , , где - константы упругости; - пористость; k - проницаемость; - коэффициент эффективных напряжений Био; - плотности скелета, присоединенной массы и жидкой среды; - объемные силы. Рассмотрим следующие типы граничных условий: , , , где и - части границы , где заданы соответствующие обобщенные перемещения и обобщенные поверхностные усилия. Вектор перемещений во внутренних точках связан с перемещениями на границе и поверхностными усилиями следующим образом: где и - компоненты тензоров фундаментальных и сингулярных решений. Аналитическое решение одномерной задачи о действии осевой силы для перемещений и порового давления имеет вид [5] Поровязкоупругое решение рассчитывается из пороупругого решения путем пересчета модулей упругости и G в области, преобразованной по Лапласу, и получении функций и . Функции и для различных моделей вязкоупругости имеют следующий вид: модель Кельвина-Фойгта , модель стандартного вязкоупругого тела , , ; модель со слабосингулярным ядром , , , где - величина, обратная характерному времени ползучести в случае модели Кельвина-Фойгта и характерному времени релаксации в случае модели стандартного вязкоупругого тела; и h - параметры ядра ползучести слабосингулярной степенной модели вязкоупругости. Численное обращение преобразования Лапласа. Рассмотрим функцию действительного переменного t. Тогда прямое и обратное преобразования Лапласа определяются формулами где s - комплексный параметр преобразования; - произвольное вещественное число, выбранное таким образом, что все особые точки функции лежат правее прямой . Пусть , тогда обратное преобразование Лапласа по методу Дурбина [31] запишется в виде , Опишем модификацию метода Дурбина на основе формул трапеций с переменным шагом для всей подынтегральной функции. Разбивая промежуток [0, R] на n частей, получаем следующие аппроксимации: , На каждом отрезке аппроксимируем отдельно реальную мнимую части функции : где , , , . Тогда использование линейной интерполянты дает следующие формулы: , С учетом того, что , и , в итоге получаем формулу . Формулы имеют такой же порядок точности, как и формула трапеций численного интегрирования, а именно ). Использование варианта метода Дурбина с переменным шагом позволяет уменьшить время, необходимое для решения задачи, не теряя при этом в точности полученного решения. При решении конкретной задачи рассматривается конечный спектр частот . 2. Гранично-элементный подход Численная схема основана на прямом подходе с использованием формулы Грина-Бетти-Сомильяны. Чтобы ввести ГЭ-дискретизацию, рассматривается регуляризованное уравнение где x - точка коллокации; y - точка наблюдения; и - обобщенные перемещения и поверхностные усилия, соответственно в случае изотропной поровязкоупругости; - соответствующие фундаментальные и сингулярные решения; - статическая часть сингулярного решения; коэффициент равен 1 в случае конечной области и равен -1 в случае бесконечной области . Будем аппроксимировать границу области совокупностью четырехугольных и треугольных восьмиузловых биквадратичных элементов. При этом треугольные элементы рассматриваются как вырожденные четырехугольные элементы. Для аппроксимации обобщенных граничных перемещений применим билинейные элементы, а для аппроксимации обобщенных поверхностных сил - постоянные элементы. Такая согласованность аппроксимаций границы области, граничных перемещений и поверхностных сил выбрана из тех соображений, что напряжения определяются через производные от перемещений, а перемещения зависят не только от координат точки, но и от конфигурации границы в окрестности этой точки. Применяются квадратуры Гаусса и иерархический алгоритм для поэлементного численного интегрирования. В качестве проекционного метода применяется метод коллокации. Для метода коллокации выберем множество узлов, совпадающее с множеством узлов аппроксимации исходных граничных функций. В итоге сформируются системы линейных алгебраических уравнений (СЛАУ) [30]. 3. Верификация методики Рассматривается задача о поровязкоупругом стержне, на который действует сила в виде функции Хэвисайда. Параметры материала: , - коэффициент Пуассона. Задача об одномерном пороупругом стержне имеет аналитическое решение [5]. Получены трехмерное численное и одномерное аналитическое поровязкоупругие решения с учетом принципа соответствия и при одинаковых параметрах модели вязкоупругости. На рис. 1 приведено сравнение аналитического и численного поровязкоупругих решений для перемещений, когда вязкоупругие свойства скелета описываются моделью со слабосингулярным ядром. Рис. 1. Сравнение численного и аналитического поровязкоупругих решений в случае модели со слабосингулярным ядром (h = 1 и ) Fig. 1. Numerical and analytical solutions comparison in case of a model with a weakly singular kernel (h = 1 and ) Одномерное аналитическое решение сравнивается с трехмерным решением при условии, что параметры материала и граничные условия подобраны специальным образом: коэффициент Пуассона равен нулю для того, чтобы провести это сравнение, поскольку аналитическое решение для пороупругости существует только в одномерном случае. Максимальная относительная погрешность по равномерной норме не превосходит 1 % на рассматриваемом временном интервале. 4. Численные результаты 1. Рассматривается задача о поровязкоупругой призматической консоли, жестко закрепленной с одного торца, и действующей осевой силе в виде функции Хэвисайда - с другого (рис. 2). За основу для поровязкоупругого материала взят песчаник Berea. Пороупругие параметры песчаника Berea: , , Для получения численного решения в точке контакта B вводится гранично-элементная сетка из 504 элементов (рис. 3). Длина консоли 3 м. Рис. 2. Постановка задачи Рис. 3. Гранично-элементная дискретизация Fig. 2. Problem statement Fig. 3. Boundary-element discretization Исследуются перемещения вдоль оси действия силы в точке B (см. рис. 2). На рис. 4 представлено сравнение трехмерных поровязкоупругих решений и пороупругого решения в случае модели со слабосингулярным ядром при различной вязкости материала. Рис. 4. Перемещения в точке B в случае модели со слабосингулярным ядром при различных параметрах Fig. 4. Displacements in point B in case of a model with a weakly singular kernel with various model parameters 2. Рассматривается задача о консоли длиной 9 м с одним жестко защемленным торцом и свободным другим. На свободный торец в начальный момент времени действует осевая сила в виде функции Хэвисайда . В качестве пороупругого материала используется песчаник Berea. Исследуются перемещения и поровое давление в точке A, удаленной на 1,5 м от точки приложения (рис. 5). Вводится гранично-элементная дискретизация (рис. 6). Рис. 5. Постановка задачи Рис. 6. Гранично-элементная дискретизация Fig. 5. Problem statement Fig. 6. Boundary-element discretization Получены поровязкоупругие решения для моделей Кельвина-Фойгта, стандартного вязкоупругого тела. На рис. 7-9 приведены графики численных решений для перемещений (в направлении действия силы) и поровых p давлений в точке А, рассмотрены различные значения параметров модели, также на графиках приведено пороупругое решение. По итогам экспериментов, результаты которых приведены на рис. 8, можно заключить, что численно продемонстрирован эффект перестройки волновых полей внутренних перемещений, когда свойства вязкоупругого материала модели стандартного вязкоупругого тела изменялись с мгновенных модулей на длительные. В откликах перемещений изменялись (увеличивались) амплитуда и период искомой функции. Эффект перестройки граничных полей численно описан ранее в [23, 30]. Рис. 7. Перемещения (а) и поровое давление p (б) в случае модели Кельфина-Фойгта Fig. 7. Displacements (a) and pore pressure p (b) in case of Kelvin-Voigt model Рис. 8. Перемещения (а) и поровое давление p (б) в случае модели стандартного вязкоупругого тела ( ) Fig. 8. Displacements (a) and pore pressure p (b) in case of a standard viscoelastic model ( ) Рис. 9. Перемещения (а) и поровое давление p (б) в случае модели со слабосингулярным ядром Fig. 9. Displacements (a) and pore pressure p (b) in case of a model with a weakly singular kernel 3. Рассмотрим задачу о действии вертикальной силы на поверхность составного пороупругого полупространства (рис. 10). Дневная поверхность полупространства свободная и проницаемая: поровое давление p = 0 и поверхностные силы кроме центрального участка площадью 1 м2, где , Н/м2 (см. рис. 9). Исследовалось поведение тела в точке дневной поверхности P(10,0,0), удаленной от источника силы на 10 м. Пороупругий материал - песчаник Berea: , , , , , , , Моделируется поверхностная волна Релея. Чтобы получить численные решения для вертикальных перемещений вводится ГЭ-дискретизация (рис. 11). Рис. 10. Постановка задачи Рис. 11. Гранично-элементная дискретизация Fig. 10. Problem statement Fig. 11. Boundary-element discretization ▬ пороупругое решение поровязкоупругое решение, полученное по модели стандартного вязкоупругого тела: ▬ ▬ ▬ ▬ Рис. 12. Вертикальные перемещения в случае модели стандартного вязкоупругого тела Fig. 12. Vertical displacements in case of a standard viscoelastic model ▬ пороупругое решение поровязкоупругое решение, полученное по модели со слабосингулярным ядром: ▬ , h = 1 ▬ , h = 1 ▬ , h = 1 ▬ , h = 1 Рис. 13. Вертикальные перемещения в случае модели со слабосингулярным ядром Fig. 13. Vertical displacements in case of a model with a weakly singular kernel ▬ пороупругое решение поровязкоупругое решение, полученное по модели со слабосингулярным ядром: ▬ , h = 1 ▬ , h = 2 ▬ , h = 5 Рис. 14. Вертикальные перемещения в случае модели со слабосингулярным ядром Fig. 14. Vertical displacements in case of a model with a weakly singular kernel На рис. 12-14 приведено сравнение пороупругих и поровязкоупругих решений для различных моделей с различными параметрами вязкоупругого материала. На рис. 12 в случае модели стандартного вязкоупругого тела наблюдается уменьшение амплитуды волны Релея при увеличении параметра вязкости материала. Аналогичная ситуация в случае модели со слабосингулярным ядром (см. рис. 13-14); наблюдается уменьшение амплитуды или же «сглаживание» волны Релея при изменении параметров модели. Однако нетипичное поведения наблюдается в случае близости к 1, данный случай требует отдельного изучения в дальнейшем. Заключение С использованием гранично-элементного подхода получены решения следующих краевых задач: о действии скачка силы на торец призматического пороупругого/поровязкоупругого тела, о действии силы в виде функции Хэвисайда на поверхность пороупругого/поровязкоупругого полупространства. Верификация представленного подхода была проведена путем сравнения с аналитическим решением. Установлено, что параметры вязкости оказывают существенное влияние на характер распределения волновых процессов. Проведено сравнение поровязкоупругих решений для ряда моделей при различных значениях параметров моделей. Численно продемонстрирован эффект перестройки волновых полей внутренних перемещений, когда свойства вязкоупругого материала модели стандартного вязкоупругого тела изменялись с мгновенных модулей на длительные. Продемонстрировано влияние вязкости на волну Релея на поровязкоупругом полупространстве.About the authors
A A Ipatov
Research Institute for Mechanics, National Research Lobachevsky State University of Nizhny Novgorod
A A Belov
Research Institute for Mechanics, National Research Lobachevsky State University of Nizhny Novgorod
S Yu Litvinchuk
Research Institute for Mechanics, National Research Lobachevsky State University of Nizhny Novgorod
References
- Frenkel J. On the theory of seismic and seismoelectric phenomena in a moist soil // Journal of Physics, - 1944. - Vol. 8. - P. 230-241.
- Biot M.A. General theory of three-dimensional consolidation // J. Appl. Phys. - 1941. - Vol. 12 (2). - P. 155-164.
- Biot M.A. Theory of deformation of a porous viscoelastic anisotropic solid // J. Appl. Phys. - 1956. - Vol. 27 (5). - P. 459-467.
- Biot M.A. Theory of elasticity and consolidation for a porous anisotropic solid // J. Appl. Phys. - 1955. - Vol. 26 (2). - P. 182-185.
- Schanz M. Wave Propagation in Viscoelastic and Poroelastic Continua. Berlin Springer, 2001. - 170 p.
- Schanz M. Poroelastodynamics: Linear models, analytical solutions, and numerical methods // Appl. Mech. Rev. - 2009. - Vol. 62 (3), - 15 p. doi: 10.1115/1.3090831
- Boer R. de. Highlights in the Historical Development of the Porous Media Theory: Toward a Consistent Macroscopic Theory' // Appl. Mech. Rev. - 1996. - Vol. 49(4). - P. 201-262.
- Nikolaevskiy V.N. Biot-Frenkel Poromechanics in Russia (Review) // J. Eng. Mech. - 2005. - Vol. 131 (9). - P. 888-897.
- Garg S.K., Nayfeh A.H., Good A.J. Compressional waves in fluid-saturated elastic porous media // J Appl Phys. - 1974. - Vol. 45. - P. 1968-1974.
- Simon B.R., Zienkiewicz O.C., Paul D.K. An analytical solution for the transient response of saturated porous elastic solids // Int J Num Anal Meth Geomech. - 1984. - Vol. 8. - P. 381-398.
- Schanz M., Cheng A. H.-D.,Transient Wave Propagation in a One-Dimensional Poroelastic Column // Acta Mech. - 2000. - Vol. 145. - P. 1-8.
- Boer R. de., Wolfqang E., Liu Z.F. One-dimensional transient wave propagation in fluid-saturated incompressible porous media // Arch. Appl. Mech. - 1993. - Vol. 63. - P. 59-72.
- Кристенсен Р. Введение в теорию вязкоупругости. - М.: Мир, 1974. - 338 с.
- Wilson R.K., Aifantis E.C. On the Theory of Consolidation With Double Porosity // Int. J. Eng. Sci. - 1982. - Vol. 20. - P. 1009-1035.
- Vgenopoulou I. Beskos D.E., Dynamic Behavior of Saturated Poroviscoelastic Media // Acta Mech. - 1992. - Vol. 95. - P. 185-195.
- Schanz M., Cheng A. H.-D. Dynamic Analysis of a One-Dimensional Poroviscoelastic Column // J. Appl. Mech. - 2001. - Vol. 68. - P. 192-198. doi: 10.1115/1.1349416#
- Cruse T.A., Rizzo F.J. A direct formulation and numerical solution of the general transient elastodynamic problem. Part I // J. Math. Anal., Applic. - 1968. - Vol. 22. - P. 244-259.
- Cruse T.A., Rizzo F.J. A direct formulation and numerical solution of the general transient elastodynamic problem. Part II // J. Math. Anal, Applic. - 1968. - Vol. 22 (2). - P. 341-355.
- Mansur W.J., Brebbia C.A. Transient Elastodynamics Using a Time-Stepping Technique // In Boundary Elements. - Berlin: Springer-Verla, 1983. - P. 677-698.
- Beskos D.E. Boundary Element Methods in Dynamic Analysis // Appl. Mech. Rev. - 1987. - Vol. 40 (1). - P. 1-23.
- Beskos D.E. Boundary element methods in dynamic analysis: Part II 1986-1996 // Appl. Mech. Rev. - 1997. - Vol. 50(3). - P. 149-197.
- Игумнов Л.А., Петров А.Н., Ипатов А.А. Сравнение численного построения оригинала решения одномерной пороупругой задачи на основе шагового метода и метода Дурбина // Проблемы прочности и пластичности: межвуз. сб. - Н. Новгород: Изд-во Нижегород. гос. ун-та, - 2013. - Вып. 75 (4). - С. 273-279.
- Игумнов Л.А., Ипатов А.А., Сабаева Т.А. Влияние вязкости на динамический отклик в вязкоупругих и поровязкоупругих телах // Проблемы прочности и пластичности: межвуз. сб. - Н. Новгород: Изд-во Нижегород. гос. ун-та. - 2014. - Вып. 76 (2). - С. 106-113.
- Igumnov L.A., Litvinchuk S.Y., Petrov A.N., Belov A.A. Boundary-element modeling of 3-D poroelastic half-space dynamics // Advanced Materials Research - 2014. - Vol. 1040. - P. 881-885.
- Numerical-analytic investigation of the dynamics of viscoelastic and porous elastic bodies / L.A. Igumnov, A.V. Amenitskii, A.A. Belov, S.Y. Litvinchuk, A.N. Petrov // Journal of Applied Mechanics and Technical Physics. - 2014. - Vol. 55 (1). - P. 89-94.
- Nardini D., Brebbia C.A. A New Approach to Free Vibration Analysis Using Boundary Elements // In Boundary Element Methods / ed. C.A. Brebbia. - Berlin: Springer-Verlag, - 1982. - P. 312-326.
- Бабешко В.А., Бабешко О.М., Евдокимова О.В. Об интегральном и дифференциальном методах факторизации // Докл. РАН. - 2006. - Т. 410, № 2. - С. 168-172.
- Ватульян А.О. О граничных интегральных уравнениях I рода в динамических задачах анизотропной теории упругости // Докл. РАН. - 1993. - Т. 333, № 3. - С. 312-314.
- Игумнов Л.А Граничные интегральные уравнения трехмерных задач на плоских волнах // Докл. РАН. - 2006. - Т. 409, № 5. - С. 1-3.
- Баженов В.Г., Игумнов Л.А. Методы граничных интегральных уравнений и граничных элементов в решении задач трехмерной динамической теории упругости с сопряженными полями. - М.: Физматлит, 2008. - 352 с.
- Durbin F. Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate's method // Computer Journal. - 1974. - Vol. 17. - No. 4. - P. 371-376.