CC BY
39
УДК 519.63
М. В. Васильева, А. Е. Колесов, А. А. Тырылгин
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЗАДАЧИ ВЯЗКОПОРОУПРУГОСТИ ГРУНТОВ
Знание процесса консолидации глинистых грунтов необходимо для правильного прогноза скорости осадок сооружений. Наиболее интенсивно вторичная консолидация отмечается в глинистых, суглинистых и песчаных отложениях. В таких грунтах уплотнение происходит за счет ползучести твердых частиц скелета и зависит от длительности действия нагрузки. В работе рассматривается численное решение задачи вторичной консолидации. Рассматривается вязкоупругая формулировка Кельвина-Фойгта для метода конечных элементов. Базовая система уравнений включает в себя нестационарные уравнения для перемещений и давления в пористой среде. Вычислительный алгоритм основан на конечно-элементной аппроксимации по пространственным переменным и неявной аппроксимации по времени. Представлены результаты численного решения модельных задач для двух случаев граничных условий (стационарная и нестационарная нагрузка). Проведено численное сравнение моделей для давления и перемещения при различных значениях параметра вторичной консолидации. Результаты представлены для двумерной и трехмерной формулировок.
Ключевые слова: задача пороупругости, ползучесть, вторичная консолидация, вязкопороупру-гость, метод конечных элементов, математическое моделирование, вязкоупругая формулировка Кельвина-Фойгта, пористая среда.
M. V. Vasilyeva., A. E. Kolesov, A. A. Tyrylgin
Numerical Simulation of the Problem of Visco-Poroelasticity Soils
Knowledge of the process of consolidation of clay soils is necessary for the correct prediction of the rate of sedimentation of structures. The most intensive secondary consolidation is noted in clayey, loamy and sandy sediments. In such soils, condensation occurs due to the creep of the solid particles of the skeleton and
ВАСИЛьЕВА Мария Васильевна - доцент-исследователь научно-исследовательской кафедры вычислительных технологий ИМИ СВФУ им. М.К. Аммосова.
E-mail: vasilyevadotmdotv@gmail.com
VASILYEVA Maria Vasilyevna - Associate Professor-Researcher of the Research Department of Computational Technologies, Institute of Mathematics and Informatics, M.K. Ammosov North-Eastern Federal University.
КОЛЕСОВ Александр Егорович - доцент-исследователь научно-исследовательской кафедры вычислительных технологий ИМИ СВФУ им. М.К. Аммосова.
E-mail: ae.kolesov@s-vfu.ru
KOLESOV Alexandr Еgorovich - Associate Professor-Researcher of the Research Department of Computational Technologies, Institute of Mathematics and Informatics, M.K. Ammosov North-Eastern Federal University.
ТЫРЫЛГИН Алексей Афанасьевич - магистрант научно-исследовательской кафедры вычислительных технологий ИМИ СВФУ им. М.К. Аммосова.
E-mail: koc9tk@mail.ru
TYRYLGIN Alexey Afanas'evich - Master's Student of the Research Department of Computational Technologies, Institute of Mathematics and Informatics, M.K. Ammosov North-Eastern Federal University.
depends on the duration of the action of the load. In this paper considers the numerical solution of secondary consolidation. We consider viscoelastic Kelvin-Voigt formulation for the finite element method. The basic system of equation includes nonstationary equations displacements and pressure in a porous medium. The Computational algorithm is based on finite element approximation in space and an implicit difference scheme for time approximation. We present the results of numerical solution of model problems for two cases of boundary conditions (stationary and non-stationary load). We perform a numerical comparison of the models for pressure and displacement for different values of the secondary consolidation parameter. The results are presented for two-dimensional and three-dimensional formulations.
Keywords: problem of poroelasticity, creep, secondary consolidation, visco-poroelasticity, finite element method, mathematical simulation, viscoelastic formulation of Kelvin, porous medium.
Введение
В работе рассмотрена модель расчета деформаций оснований, сложенных вязкоупру-гими грунтами, к которым можно отнести, например, глинистые грунты. Поведение таких грунтов зависит от длительности действия нагрузки и обуславливается ползучестью скелета грунта (реологическими свойствами). В водонасыщенных грунтах при приложении нагрузки начинается процесс фильтрационной консолидации, т. е. по мере отжатия из грунта воды поровое давление падает, и часть нагрузки переходит на скелет грунта. При этом в грунте происходят процессы как первичной, так и вторичной консолидации.
Присутствие процесса вторичной консолидации может быть объяснено вязкоупругим поведением скелета грунта [1].
Процесс фильтрационной консолидации был впервые исследован Терзахи [2] для одномерного случая, обобщенная феноменологическая модель была предложена и проанализирована Био [3-6].
Базовые математические модели включают уравнения Ламе для перемещений среды, а также закон сохранения массы и закон Дарси для фильтрующейся жидкости [7-10]. Наиболее важной особенностью математических моделей пороупругости является то, что уравнения системы связаны между собой.
В качестве модели вторичной консолидации мы использовали модель Кельвина-Фойгта, в которой имеется дополнительное слагаемое, отвечающее на реологические свойства грунта [11].
Работа состоит из четырех частей. В первой части рассмотрена постановка задачи вторичной консолидации. Во второй части приводится конечно-элементная аппроксимация уравнений. Третья часть включает в себя три раздела: в первом разделе представлено численное решение задачи вязкопороупругости с приложением постоянного граничного условия; во втором разделе рассмотрена задача с граничным условием, зависящим от времени; в третьем разделе представлены расчеты для трехмерной постановки. В конце представлено заключение.
Постановка задачи
Рассматривается система уравнений, описывающая процесс вторичной консолидации в вязкоупругой пористой среде. Система уравнений в частных производных включает в себя уравнение для перемещения и давления
diva + agradp = fu (x,t), (1)
д
a—(divu)-div(kgradp) = fp (x,t)x efi,0< t<T, (2)
где u - вектор перемещения, p - поровое давление жидкости, а - коэффициент Био, s - тензор напряжений, k - проницаемость среды и f - внутренние источники (I = p,u).
Классической моделью, используемой для описания вязкоупругой среды, является модель Кельвина-Фойгта, в которой используются следующие соотношения для тензора деформаций а [12-13]:
е=ее = е\ а=ае +ст\ (3)
где индексы V и е используются для обозначения вязкоупругой и упругой компонентов тензоров, соответственно
г (и) =1 (grad и + grad ит).
Упругое напряжение может быть написано, с точки зрения компонентов тензора деформации, следующим образом:
= СуеЫ = СуеЫ . (4)
Аналогичным образом нижеследующее дифференциальное соотношение дает вязкие компоненты напряжения:
V 1т - V 1т -
=1у£1т =1у£1т . (5)
где ё — производная по времени от тензора деформации.
В уравнении (3) коэффициенты С1™ и п'™ определяются следующим образом:
С; =Щ81т + ^(8и8]т +8т8]1), (6)
П = УхЩ81т + 81т8]1), (7)
где i и] обозначают индексы в Декартовых координатах (/', ] = 1,2 для двумерного случая), Ух и Ур - параметры вторичной консолидации, 1 и т - коэффициенты Ламе
з _ уЕ _ Е
(1 + у)(1 - 2у)' 2 (1 +у)' (8)
Здесь Е и V, соответственно, являются модулем Юнга и коэффициентом Пуассона [14].
Задачи для системы уравнений (1)-(2) рассматриваются в ограниченной области О на границе, которой задаются следующие граничные условия. Для перемещений на границе = Гв + Г, положим
и = 0, хеГд, -ап = 0, хеГы, (9)
где п - единичная нормаль к границе. Аналогично для давления задаются граничные условия первого и второго рода на части границ области:
dp / \
р = 0, х еГ2 иГ3 иГ4^ = р1 (х), х еГ1. (10)
dn
На начальный момент времени ^ = 0) задается условие
е(х,0 ) = 0, р (х,0 ) = р0, х еО (11)
для деформаций и давления.
Аппроксимация по пространству
Для аппроксимации рассматриваемой системы уравнений для давления и перемещений используется метод конечных-элементов. Определим гильбертово пространство Ь2(П), в котором скалярное произведение и норма определяются в виде
Г II II 1
(и, V) = \и (х ^ (х X, р| = (и, и )2 .
п
Для векторных величин используем [Ь2(Л)]а, где ё = 2,3 - размерность пространства. Определим следующие подпространства скалярных и векторных функций:
V = { е Ь1 (П):q(х) = р,х е дП}, Ж = {у е ] (П): V(х) = 0,х е Гв},
V = { е (О): q (х) = 0, х едО}, Ж = Ж.
Запишем вариационную формулировку для уравнений (1)-(2) с начальными граничными условиями (9)-(11): найти и е Ж, р е Vтакие, что
г (", V) + а (и, V) + g (р, V ) = I" (у), V еЖ', (12)
d(и, q) + Ь(р, q) = 1Р ^), q eV', (13)
где билинейные и линейные формы задаются следующим образом:
м
дг
, ч еда4 (и (у ) , ч (• , , , ч -(и, у ) = J-——— ах, а (и, у ) = ]а (и )е(у )dx,
о
§(р,V) = аjgrad pvdx, d(и, q) = а|ди qdx,
о
дг
Ь (р, q) = grad p,grad q )),
о
Г () = ¡Л^х 1Р () = \/р^х +1 рИ^.
О П дП
Билинейные формы г(у), а(у) и Ь(у) симметричны и положительно определены, а ё (у, у) = - g (у, у) для V е V, V е Ж
Определим конечно-мерное пространство для скалярных и векторных функций Vh с V Жь с Ж и запишем конечно-элементную аппроксимацию следующим образом: найти и, е Ж,, р, е V, такие, что
п п7 Г п п 7
г (и,, у) + а (и,, у) + g (р,, у) = Iй (у), V ^, (14)
1 (щ,q) + ь(,q) = 1Р q . (15)
Отметим, что для аппроксимации по времени мы будем использовать неявную разностную схему с шагом , и будем использовать линейные базисные функции как для давления, так и для перемещений.
Вычислительная реализация
В данной части работы проведем численное моделирование задачи вторичной консолидации для задачи пороупругости. Расчеты будем проводить на неструктурированной расчетной сетке, которая представлена на рис. 1. Треугольная сетка содержит 492 вершин, 910 ячеек и 1401 граней.
Численное решение представлено для следующих граничных условий
Рис. 1. Расчетная сетка
их = 0, = 0, х е Г3
иу = 0, ох = 0, х е Г4
^х = = 0 Х еГ2 иГ1>
для уравнения вязкоупругости и для уравнения зададим р = 0 на границах Г2 иГ3 и Г4. На границе Г1 рассмотрим два случая:
• с постоянной нагрузкой по времени р = 100 и;
• с нагрузкой, зависящей от времени р = 100sin(t)
Для расчетов использовались обезразмеренные коэффициенты А = ц = к = а = 1 с областью моделирования [0,1] х [0,1] пористого грунта. Расчет проведем Т = 0.5 с шагом по времени т = 0.01 и Т = п с шагом по времени т = 0.02п постоянной нагрузки и зависящей от времени нагрузки соответственно.
Численное решение с постоянной нагрузкой
Рассмотрим влияние параметра вторичной консолидации у = ух = у, где у = 0.01 и 0. Здесь у = 0 соответствует модели пороупругого грунта (без учета вязкости). На рис. 2 показаны распределения давлений и перемещений по X и У направлениям в последний момент времени для коэффициента у = 0.1.
(МП 0№
......
0.0940
Рис. 2. Распределение давления, перемещения по X и У направлениям на последний момент времени для коэффициента у = 0.1 и (слева направо). Модель с постоянными по времени граничными условиями
Рис. 3. Сравнение моделей пороупругости и вторичной консолидации (вязкопороупругость) с коэффициентом у = 0.1 и у = 0.01. Относительная Ь2 норма разности решений в процентах. Давление и перемещения по X и У направл ениям (слева направо). Модель с постоянными по времени граничными условиями
Рис. 4. Одномерные графики среза вдоль вертикальной серединной линии для перемещения по Х(слева) и Y(справа) направлениям в моменты времени 0.1, 0.3 и 0.5. Модель с постоянными по времени граничными условиями
Проведем численное сравнение моделей с пороупругостью и пороупругостью со вторичной консолидацией. На рис. 3 показано сравнение разности решений пороупругой и вяз-копороупругой моделей по времени. Для сравнения вычислялась относительная Ь2 норма разницы решений для давления и перемещений по X и У направлениям
\\и - и.
2
1п( ие - и) ) dx
[ и2dx J п е
где и и и решения для упругой и вязкоупругой моделей соответственно.
Одномерные графики среза вдоль вертикальной серединной линии для перемещения по X, У направлениям отображены на рис. 4. Графики иллюстрируют различие решений для
Рис. 5. Распределение давления, перемещения по X и У в момент времени 0.8п (слева направо). Модель с граничными условиями, зависящими от времени
Рис. 6. Сравнение моделей пороупругости и вторичной консолидации (вязкопороупругость) с коэффициентом у = 0.1 и у = 0.01. Относительная Ь2 норма разности решений в процентах. Давление и перемещение по X и У. направл ениям (слева направо). Модель с граничными условиями, зависящими от времени
пороупругой и вязкопороупругой среды при различных значениях коэффициента вторичной консолидации у = 0.01, различия практически незаметны на срезах перемещений (рис. 4). Отметим, что при этом на рис. 3 разница решений для у = 0.01 в начальные моменты времени составляет около 10 % для перемещений.
Численное решение с нагрузкой, зависящей от времени
Далее проведем численное моделирование с использованием граничного условия, зависящего от времени. На рис. 5 показаны распределения перемещений по X, Y и давления в момент времени 0.8п. Проведем численное сравнение между моделями пороупругости и пороупругости со вторичной консолидацией (вязкопороупругость). На рис. 6 представлены результаты сравнения разности решений пороупругой и вязкопороупругой моделей по времени при различных значениях коэффициента у = 0.1 и 0.01. Одномерные графики среза вдоль вертикальной серединной линии для перемещения по X, Y направлениям отображены на рис. 7.
Рисунки иллюстрируют различие решений в различные моменты времени. Анализируя рис. 6, можно сделать вывод, что при нагрузке, зависящей от времени, наибольшее различие решений для перемещений возникает на начальные и конечные моменты времени. При этом различие только в начальный момент времени несущественно, поскольку разница составляет менее 10 %. Разница решений для перемещений является достаточно большой и составляет около 50 % при большом значении у и падает до 10 % при уменьшении у.
Рис. 7. Одномерные графики среза вдоль вертикальной серединной линии для перемещения по Х(слева) и Y(справа) направлениям в моменты времени 0,5п, 0.8п и п. Модель с граничными условиями, зависящими от времени
Рис. 8. Геометрическая область в тетраэдральной расчетной сетке Численные результаты для трехмерной задачи
Проведем численное моделирование рассмотренной задачи вязкопороупругости в трехмерной постановке. Моделирование будем проводить в области [0,1] х [0,1] х [0,1] метра при постоянных граничных условиях для давления р = 100 и у = 0.01. Остальные граничные условия были выбраны аналогично двумерной постановке. Расчет проведем при Т = 0.1 с шагом по времени т = 0.01. Расчетная тетраэдральная сетка содержит 1265 вершин, 6472 ячеек и 13416 граней с размером области нагрузки 0.2 х 0.2 (рис. 8).
Рис. 9. Распределение давления и перемещения по Х, Y и Z (слева направо) в конечный момент времени
Численные результаты задачи вязкопороупругости в трехмерной постановке отображены на рис. 9 и иллюстрируют решение задачи на последний момент времени. Время счета составило 77.71 секунды, а количество итераций - 60 с использованием итерационного метода TFQMR с предобуславливателем SOR.
Для построения геометрической области используется программа Gmsh [15]. Визуализацию полученных результатов осуществляем с использованием программы Paraview и Gnuplot [16-17]. Вычислительная реализация построена с использованием библиотеки FeniCS [18].
Заключение
В данной работе мы рассмотрели задачу вязкопороупругости в ограниченной области
При построении математической модели для задачи вторичной консолидации была использована вязкоупругая формулировка Кельвина-Фойгта. Для численного решения возникающей связанной системы уравнений для давления и перемещений построена аппроксимация задачи посредством метода конечных элементов с использованием линейных базисных функций. Проведено численное исследование задачи при различных граничных условиях и построены графики численного сравнения моделей для давления и перемещения при различных коэффициентах вторичной консолидации у . Моделирование проведено как для двумерного, так и для трехмерного случаев.
Отметим, что расчеты проводились для обезразмеренных величин, и в будущем планируется провести исследование моделей при реальных параметрах вязкоупругой пористой среды. Также планируется рассмотреть случай неоднородных коэффициентов упругости, в частности, рассмотреть неоднородности, связанные со слоистым строением грунтов, а также при наличии сети трещин. Для трехмерной задачи планируется рассмотреть параллельные методы решений на вычислительном кластере.
Работа выполнена при поддержке проекта РФФИ №17-01-00732a и мега-гранта Правительства РФ №14. Y26.31.0013.
Л и т е р а т у р а
1. Gaspar F. J. et al. A stabilized method for a secondary consolidation Biot's model // Numerical Methods for Partial Differential Equations. - 2008. - Т. 24, №. 1. - С. 60-78.
2. Terzaghi K. Theory of consolidation. - John Wiley & Sons, Inc., 1943. - С. 265-296.
3. Biot M. A. General theory of three-dimensional consolidation // Journal of applied physics. - 1941. - Т. 12, №. 2. - С. 155-164.
4. Meirmanov A. Mathematical models for poroelastic flows. - Atlantis Press, 2014.
5. Nikolaevskii V. N. Mechanics of porous and fissured media // Nedra Moscow. - 1984.
6. Kolesov A. E., Vabishchevich P. N., Vasilyeva M. V. Splitting schemes for poroelasticity and thermo-elasticity problems // Computers & Mathematics with Applications. - 2014. - Т. 67, №. 12. - С. 2185-2198.
7. Sivtsev P. V., Vabishchevich P. N., Vasilyeva M. V. Numerical simulation of thermoelasticity problems on high performance computing systems / International Conference on Finite Difference Methods. - Springer
International Publishing, 2014. - C. 364-370.
8. Kolesov A. E., Vabishchevich P. N., Vasilyeva M. V. Splitting schemes for poroelasticity and thermo-elasticity problems // Computers & Mathematics with Applications. - 2014. - T. 67, №. 12. - C. 2185-2198.
9. Brown D. L., Vasilyeva M. A generalized multiscale finite element method for poroelasticity problems I: linear problems // Journal of Computational and Applied Mathematics. - 2016. - T. 294. - C. 372-388.
10. Brown D. L., Vasilyeva M. A generalized multiscale finite element method for poroelasticity problems II: Nonlinear coupling // Journal of Computational and Applied Mathematics. - 2016. - T. 297. - C. 132-146.
11. Mesquita A. D., Coda H. B. Alternative Kelvin viscoelastic procedure for finite elements // Applied Mathematical Modelling. - 2002. - T. 26, №. 4. - C. 501-516.
12. Lemaitre J., Chaboche J. L. Mechanics of solid materials. - Cambridge university press, 1994.
13. Rapp P., Kurzyka J., Szostak W. The creep and relaxation in sandwich panels with the viscoelastic cores // Light-weight Steel and Aluminium Structures, Elsevier, Amsterdam. - 1999. - C. 197-204.
14. Salencon J. Handbook of Continuum Mechanics: General Concepts Thermoelasticity. - Springer Science & Business Media, 2012.
15. npoipaMMHtm naKeT GMSH (http://geuz.org/gmsh/).
16. npoipaMMHtm naKeT Paraview (http://www.paraview.org).
17. npoipaMMHtm naKeT Gnuplot (http://www.gnuplot.info).
18. BLmHcnHTenLHLiH naKeT FEniCS (http://fenicsproject.org).
R e f e r e n c e s
1. Gaspar F. J. et al. A stabilized method for a secondary consolidation Biot's model // Numerical Methods for Partial Differential Equations. - 2008. - T. 24, №. 1. - S. 60-78.
2. Terzaghi K. Theory of consolidation. - John Wiley & Sons, Inc., 1943. - S. 265-296.
3. Biot M. A. General theory of three-dimensional consolidation // Journal of applied physics. - 1941. - T. 12, №. 2. - S. 155-164.
4. Meirmanov A. Mathematical models for poroelastic flows. - Atlantis Press, 2014.
5. Nikolaevskii V. N. Mechanics of porous and fissured media // Nedra Moscow. - 1984.
6. Kolesov A. E., Vabishchevich P. N., Vasilyeva M. V. Splitting schemes for poroelasticity and thermoelasticity problems // Computers & Mathematics with Applications. - 2014. - T. 67, №. 12. - S. 2185-2198.
7. Sivtsev P. V., Vabishchevich P. N., Vasilyeva M. V. Numerical simulation of thermoelasticity problems on high performance computing systems / International Conference on Finite Difference Methods. - Springer International Publishing, 2014. - S. 364-370.
8. Kolesov A. E., Vabishchevich P. N., Vasilyeva M. V. Splitting schemes for poroelasticity and thermoelasticity problems // Computers & Mathematics with Applications. - 2014. - T. 67, №. 12. - S. 2185-2198.
9. Brown D. L., Vasilyeva M. A generalized multiscale finite element method for poroelasticity problems I: linear problems // Journal of Computational and Applied Mathematics. - 2016. - T. 294. - S. 372-388.
10. Brown D. L., Vasilyeva M. A generalized multiscale finite element method for poroelasticity problems II: Nonlinear coupling // Journal of Computational and Applied Mathematics. - 2016. - T. 297. - S. 132-146.
11. Mesquita A. D., Coda H. B. Alternative Kelvin viscoelastic procedure for finite elements // Applied Mathematical Modelling. - 2002. - T. 26, №. 4. - S. 501-516.
12. Lemaitre J., Chaboche J. L. Mechanics of solid materials. - Cambridge university press, 1994.
13. Rapp P., Kurzyka J., Szostak W. The creep and relaxation in sandwich panels with the viscoelastic cores // Light-weight Steel and Aluminium Structures, Elsevier, Amsterdam. - 1999. - S. 197-204.
14. Salencon J. Handbook of Continuum Mechanics: General Concepts Thermoelasticity. - Springer Science & Business Media, 2012.
15. Programmnyi paket GMSH (http://geuz.org/gmsh/).
16. Programmnyi paket Paraview (http://www.paraview.org).
17. Programmnyi paket Gnuplot (http://www.gnuplot.info).
18. Vychislitel'nyi paket FEniCS (http://fenicsproject.org).
