The Simulation of the Nonlinear Consolidation of Porous Media
- Authors: Sheshenin S.V1, Artamonova N.B1
- Affiliations:
- Lomonosov Moscow State University
- Issue: No 1 (2022)
- Pages: 167-176
- Section: ARTICLES
- URL: https://ered.pstu.ru/index.php/mechanics/article/view/2017
- DOI: https://doi.org/10.15593/perm.mech/2022.1.13
- Cite item
Abstract
In this paper, the general formulation of the problem of coupled deformation of a porous deformable medium with a fluid flowing through the pores is formulated, mathematically investigated and numerically implemented within the framework of physical and geometric nonlinearity. We present the formulation of the problem in velocities of solid phase displacements and the rate of pore pressure change in differential and variational forms. A phenomenological approach was used to formulate the mechanical model. The equations of the coupled consolidation model were derived from the general conservation laws of continuum mechanics using spatial averaging over a representative volume element. The consolidation model took into account the change in the porosity and permeability of the medium during deformation. The equations of filtration and porosity change, originally presented in Euler approach, were reformulated in Lagrangian coordinates of the solid phase using the relative fluid velocity according to ALE (Arbitrary Lagrangian - Eulerian) approach. The Gâteaux differentiation technique was used to linearize the variational equilibrium equations. For spatial discretization of the saddle system of equations, the finite element method (FEM) was used: quadratic serendipity elements for approximating the equilibrium equations and Brick type elements for approximating the filtration equation. To solve the system of equilibrium and filtration equations, a generalization of the implicit scheme with internal iterations at each time step by the Uzawa method was used. The results of numerical simulation of elastoplastic deformation of a water-saturated soil under load with fluid outflow are presented. To simulate the constitutive relations of elastoplastic deformation of soil under short-term loads, a generalization of S.S. Grigoryan's model to large deformations is proposed. The calculations were carried out in our own program code. The developed consolidation model can be used to simulate the formation of tracking ruts and unevenness of natural roads, as well as to calculate the uneven settlement of engineering structures.
Full Text
Введение Консолидация - это процесс уплотнения пористого водонасыщенного грунта под нагрузкой при возможности оттока жидкости. Грунты являются естественным природным пористым материалом, который состоит из твердого скелета с порами, заполненными жидкостью. Твердая и жидкая фазы грунта взаимодействуют друг с другом. С одной стороны, увеличение порового давления приводит к росту сжимающих напряжений в скелете грунта и к деформированию скелета. С другой стороны, сжатие скелета грунта может вызвать уменьшение пористости и, как следствие, увеличение давления флюида. Эти связанные явления совместного упругого деформирования грунта и фильтрации жидкости были впервые описаны в теории линейной пороупругости М. Био [1], в которой соединены уравнение равновесия для твердого скелета и закон линейной фильтрации Дарси. Однако процесс деформирования флюидонасыщенных дисперсных грунтов под нагрузкой, как правило, нелинейный и сопровождается пластическими деформациями. В песчаных и глинистых грунтах под действием нагрузки деформации развиваются во времени с гидродинамическим запаздыванием и зависят от скорости оттока жидкости из грунта. Реакция скелета грунта на изменение порового давления может оказаться критичной, например, при ухудшении ровности дорожного покрытия и образовании колеи или при неравномерной осадке инженерных сооружений. Следовательно, разработка строгой и эффективной нелинейной связанной модели консолидации грунта, учитывающей упругопластическое поведение материала, сопровождающееся изменением пористости и проницаемости [2], весьма актуальна и имеет большое практическое значение. Большинство разработанных моделей консолидации базируются на линейной теории упругости или упругопластичности для малых деформаций [1; 3-5]. В настоящее время в литературе представлены обобщения классических уравнений консолидации с бесконечно малых на конечные деформации, например в работах [6-9]. Обобщения постановки классических уравнений консолидации на большие деформации в основном базируются на использовании определяющих соотношений в скоростях [10]. Альтернативная формулировка для упругопластичности с конечными деформациями базируется на мультипликативном разложении градиента деформации [6-9]. Задачи совместного деформирования грунта и движения жидкости могут решаться в несвязанной постановке (упругий режим фильтрации) и в связанной постановке. При решении задач об откачке нефти из глубоких скважин упругий режим фильтрации в отдельных случаях бывает вполне достаточен. Задачи о деформировании водонасыщенного грунта под нагрузкой (задачи консолидации), как правило, решаются в связанной постановке. Наиболее часто для решения подобных задач используется метод конечных элементов (МКЭ) [6-9, 11-13]. В последнее время взаимодействие между течением жидкости и скелетом грунта моделируется с использованием различных связанных схем. Уравнение фильтрации и уравнение равновесия среды могут решаться на каждом шаге по времени в виде общей системы [14; 15]. В этих работах метод решения называется полностью связанным. В другом варианте метода решения уравнение фильтрации и уравнение равновесия связываются последовательно в итерационном процессе типа метода Узавы. Сначала решается или задача фильтрации, или механическая задача, а потом решается другая задача, используя результат решения первой [12-17]. Если заданы достаточно жесткие критерии сходимости, итерационно связанное решение может быть таким же точным, как и полностью связанное, но количество требуемых итераций может быть чрезмерно большим [5]. Итерационно связанные последовательные методы могут использовать независимые решатели механики сплошных сред и средства моделирования пластового потока, развитые в нефтяной промышленности [17]. Численное решение уравнений трехмерной консолидации Био по-прежнему остается сложной задачей, несмотря на интенсивные исследования в этой области [5]. Во-первых, полностью связанный подход приводит к большим алгебраическим системам, которые могут быть довольно сложными для решения. При этом в реальных трехмерных приложениях количество неизвестных может достигать нескольких сотен тысяч. Во-вторых, матрица коэффициентов, полученная в результате численной дискретизации, может быть плохо обусловлена. Следовательно, для работы с полностью связанной схемой необходимы продвинутые решатели [18; 19]. При использовании итерационно связанной схемы сходимость итерационного процесса, зависящая от механических свойств жидкости и грунта, может быть медленной [16]. Третья проблема связана с различными формами неустойчивости при численном решении, которые проявляются в виде ложных осцилляций порового давления. Это может происходить, если сжимаемость жидкости достаточно мала по сравнению с объемной сжимаемостью скелета грунта [5]. Как результат, на начальном этапе скелет грунта ведет себя как несжимаемая среда. Тогда вначале скорость изменения объема скелета грунта близка нулю, что приводит к патологии конечных элементов (КЭ), известной как «запирание» (locking), которое обычно вызывает ложные осцилляции порового давления [5]. Справиться с этой проблемой можно, используя аппроксимации разного порядка для перемещений и порового давления. Известно, что поровое давление следует аппроксимировать линейными конечными элементами, а перемещения - квадратичными КЭ (см., например, [20]). В нашей работе формулируется физически и геометрически нелинейная связанная модель консолидации в скоростях, для решения которой используется итерационный метод Узавы. Для пространственной дискретизации седловой системы уравнений используются квадратичные серендиповы элементы для аппроксимации приращений перемещений и элементы трилинейного типа для аппроксимации изменения давления жидкости. Известно, что при решении связанных динамических задач консолидации различные исследователи часто используют метод ALE (Arbitrary Lagrangian - Eulerian) для исправления сильно искаженной лагранжевой сетки [21; 22], что приводит к перестроению сетки и, соответственно, к некоторой потере точности решения. Следует отметить, что в нашей работе метод ALE используется только для объединения Эйлерова и Лагранжева подходов и не связан с перестроением лагранжевой сетки. 1. Математическая постановка нелинейной связанной задачи консолидации и метод решения Дифференциальная постановка связанной задачи консолидации «в скоростях» включает три уравнения - равновесия, фильтрации и пористости [23-25]: (1) Система дифференциальных уравнений (1) сформулирована относительно дифференциалов вектора перемещений давления жидкости и пористости . Здесь - лагранжевы координаты твердого скелета, и - касательные модули соответственно в текущей и начальной конфигурациях, F - деформационный градиент, , - второй тензор эффективных напряжений Пиола - Кирхгофа, E - тензор деформаций Грина - Лагранжа, f - массовая сила, - эффективные напряжения Коши (осредненные напряжения при нулевом поровом давлении), - средняя плотность пористого водонасыщенного грунта, - тензор коэффициентов фильтрации, - параметр Био, - удельный вес жидкости, - сжимаемость жидкости, - пористость. Уравнения системы (1) выведены из закона фильтрации Дарси, уравнения состояния жидкости, уравнений неразрывности и уравнения равновесия [23; 25]. При выводе были сделаны следующие предположения и допущения. Процесс деформирования грунта изотермический. Поры грунта полностью заполнены флюидом. Жидкость ньютоновская и однородная, изменение плотности жидкости при изменении порового давления подчиняется закону баротропии . Удельный вес жидкости принимается постоянным при конкретной температуре и минерализации. Пористость зависит от перемещений в грунте: коэффициенты фильтрации зависят от пористости: В уравнениях неразрывности материал скелета грунта принимается несжимаемым по сравнению с сжимаемостью жидкости, так как плотность минеральных зерен rs, в отличие от плотности жидкости rf, слабо зависит от изменения давления [26]. Несжимаемость материала скелета не исключает объемную деформацию осредненного каркаса за счет перемещения частиц твердой фазы грунта или переупаковки зерен в дисперсных несвязных грунтах. Отметим, что в уравнении равновесия сжимаемость материала скелета учитывается. Уравнения фильтрации и пористости в (1) были переформулированы в движущейся системе координат, связанной с лагранжевыми координатами x твердого скелета согласно методу ALE (Arbitrary Lagrangian-Eulerian) [27] с использованием относительной скорости движения жидкости (скорости флюида относительно скелета грунта) [23; 25]. В задаче консолидации задавались следующие граничные условия. Часть границы неподвижна и непроницаема, остальная часть границы проницаема, и на ней действует поверхностная нагрузка: , (2) где n - внешняя нормаль к границе области , S - поверхностная сила. Для дальнейшего решения задачи консолидации численными методами была получена вариационная постановка в текущей конфигурации, а уравнение равновесия было линеаризовано с использованием техники дифференцирования по Гато [23; 25]. Системе (1) с граничными условиями (2) соответствует линеаризованная вариационная постановка задачи консолидации [25]: (3) Здесь S - поверхностная сила из граничных условий (2), w, q, h - пробные функции из подпространств Соболева, удовлетворяющих заданным граничным условиям 1-го рода: В качестве определяющих соотношений в (3) могут использоваться соотношения для гиперупругого материала при существовании упругого потенциала, а также теория пластического течения или деформационная теория пластичности. Сначала решаются первые два уравнения (равновесия и фильтрации) системы (3) в предположении постоянной пористости. Для решения используется обобщение неявной схемы, предложенной еще в [28; 29], с внутренними итерациями на каждом шаге по времени по методу Узавы [30]. Для реализации итерационного процесса на каждом шаге по времени применялся метод простой итерации с использованием предобусловливателя для ускорения сходимости. После решения этой задачи вычисляется пористость из третьего уравнения системы (3). Далее посчитанное значение пористости учитывается при решении уравнений равновесия и фильтрации на следующем временном шаге. Дискретизация линеаризованных вариационных уравнений равновесия и фильтрации по пространственным координатам осуществлялась с помощью метода конечных элементов (МКЭ). При дискретизации использовались 20-узловые серендиповы квадратичные элементы (Q2) для аппроксимации приращения перемещений скелета грунта и 8-узловые элементы трилинейного типа (Q1) для аппроксимации приращения давления воды [24; 25]. Такой выбор конечных элементов был обусловлен тем, что система уравнений равновесия и фильтрации является седловой, а при численном решении дискретизированной седловой системы возникает проблема корректности. Достаточным условием корректности этой системы является LBB-условие [31], выполнение которого зависит от выбранных при дискретизации конечных элементов для аппроксимации и . Как показывает опыт моделирования седловых задач [7; 32], аппроксимация квадратичными конечными элементами Q2 перемещений скелета грунта и трилинейными КЭ Q1 давления жидкости обеспечивает корректность численного решения задачи. А выполнение LBB-условия для элементов Q1 - Q1 зависит от соотношения шага по времени и шага сетки h: численное решение седловой задачи корректно при [33; 34]. Уменьшение и, как следствие, нарушение LBB-условия приводит к возникновению ложных осцилляций порового давления, которые появляются на некотором временном шаге и далее усиливаются с каждым шагом по времени, в результате чего решить седловую систему становится невозможно [24; 25]. 2. Пример численного решения задачи об упругопластическом деформировании водонасыщенного грунта Рассмотрим пример численного решения задачи о деформировании водонасыщенного песчаного грунта под действием перфорированного штампа. Для моделирования упругопластического деформирования грунта при кратковременных нагрузках предложено обобщение модели С.С. Григоряна [35] на большие деформации. Для формулировки соотношений теории пластичности использовался тензор Кирхгофа , как предлагается в работах [36-38]. Функция нагружения имеет вид: (4) Здесь - девиатор тензора напряжений, - шаровая часть тензора напряжений. Экспериментальные исследования песчаных грунтов [39] подтвердили справедливость условия пластичности f (4), а также позволили установить, что для песчаных грунтов в реальных условиях при функция F имеет вид где A и M - константы. В работе [35] для сдвиговых свойств выбрана модель идеальной пластичности, но в численных расчетах удобно использовать малое изотропное упрочнение для устойчивости вычислений, как показано в формуле (4). В критерии Кулона - Мора, который лучше всего описывает свойства слабо связных грунтов, параметры A и M из (4) выражаются формулами [40]: где - угол внутреннего трения грунта, c - сцепление. Рис. 1. Зависимость объемных напряжений и деформаций при упругопластическом деформировании грунта Fig. 1. Dependence of volumetric stresses and strains during elastoplastic deformation of the soil Сделаем допущение, что скорость объемной деформации близка к постоянной: Для этой скорости можно принять при нагрузке . Диаграмма для шаровых частей (рис. 1) также соответствует конкретной скорости нагружения. В настоящий момент в расчетах упругопластического деформирования грунта используются определяющие соотношения деформационного типа (5) в рамках предположения, что для деформаций деформационная теория сохраняет приемлемую точность, за исключением областей вблизи угловых точек, где имеют место большие деформации и сложное нагружение. (5) где - функция Ильюшина, - компоненты девиатора накопленной деформации Грина - Лагранжа, - интенсивность накопленной деформации, - мера деформаций Коши - Грина, - касательный модуль на диаграмме - предел текучести при - упругий модуль сдвига, - касательный модуль сдвига, соответствующий малому упрочнению и выбираемый как вычислительный параметр для улучшения сходимости, - тангенс угла внутреннего трения. По деформационной теории решена модельная задача о действии перфорированного штампа на область песчаного грунта кубической формы с ребром 1 м, штамп - квадратный со стороной 0,4 м. Расчеты выполнялись для деформаций 6, 8 и 10 %. На верхней границе области задавалось перемещение штампа равномерно на первых шести, восьми и десяти шагах по времени соответственно. Остальная часть верхней границы свободна от напряжений. Боковые и нижняя границы области непроницаемы и закреплены в направлении, перпендикулярном границе. Вся верхняя граница проницаема, через штамп также может просачиваться жидкость. Общее время - 47,5 с, число шагов по времени - 20, величина шага по времени - 2,375 с. Зависимость коэффициента фильтрации от пористости n (выраженной в долях единиц) задавалась по формуле Козени [41; 42]. Свойства грунта задавались по литературным источникам [43], как для среднезернистого песка: (начальный коэффициент фильтрации до деформирования), (начальная пористость), (рис. 1), (для функции нагружения (4)), Расчеты проводились в собственном программном коде. Результаты расчетов приведены на рис. 2 и 3. На рис. 2 показаны графики зависимости давления воды от времени в двух точках, находящихся на разной глубине - 0,1 м (рис. 2, а) и 0,4 м (рис. 2, б) под серединой штампа при общей деформации 6, 8 и 10 %. На рис. 3 показано распределение давления воды на 10-м шаге по времени, когда только произошла деформация грунта 10 % в результате опускания штампа. На рис. 3, а показано распределение давления во всей области, а на рис. 3, б - в вертикальной плоскости, проходящей через середину штампа. На рис. 4 приведено распределение порового давления в вертикальной плоскости на 11-м (рис. 4, а) и 15-м (рис. 4, б) шагах по времени. Можно заметить, что сразу после опускания штампа непосредственно в области под штампом формируется зона повышенного порового давления (см. рис. 3, б). На следующих шагах по времени (см. рис. 4), по мере просачивания воды через перфорированный штамп, давление воды под штампом постепенно падает, распределение порового давления в исследуемой области становится более равномерным. На 11-м шаге по времени величины порового давления снижаются примерно в 2 раза (см. рис. 4, а), на 15-м шаге по времени - становятся ниже на два порядка в результате оттока жидкости через верхнюю границу модели (см. рис. 4, б). а б Рис. 2. Зависимость давления воды от времени в точке под серединой штампа на глубине: а - 0,1 м; б - 0,4 м Fig. 2. Dependence of water pressure on time at a point under the middle of the stamp at a depth: а - 0.1 m; б - 0.4 m а б Рис. 3. Распределение давления воды на 10-м шаге по времени в области (а), в плоскости (б) (общая деформация 10 %) Fig. 3. Water pressure distribution at the 10th time step in the area (а), in the plane (б) (total deformation 10 %) а б Рис. 4. Распределение давления воды в плоскости на 11-м (а) и 15-м (б) шагах по времени (общая деформация 10 %) Fig. 4. Water pressure distribution in the plane at the 11th time step (а), at the 15th (б) time step (total deformation 10 %) Заключение В работе получена физически и геометрически нелинейная связанная формулировка задачи консолидации в скоростях. Для описания процесса деформирования твердой фазы грунта применялся лагранжев подход, для описания движения жидкости - метод ALE. С помощью метода ALE моделируется относительное движение жидкости через подвижную лагранжеву сетку для твердого каркаса. Для численного решения нелинейной задачи консолидации была получена линеаризованная вариационная постановка. Связывание уравнений деформирования скелета грунта и фильтрации жидкости осуществляется методом Узавы. Решена модельная задача об упругопластическом деформировании пористого водонасыщенного песчаного грунта под действием перфорированного штампа. Для моделирования определяющих соотношений упругопластического деформирования грунта используется обобщение модели С.С. Григоряна и деформационной теории пластичности на конечные деформации. В будущем при решении практических задач предполагается использовать теорию течения.About the authors
S. V Sheshenin
Lomonosov Moscow State University
N. B Artamonova
Lomonosov Moscow State University
References
- Biot M.A. General theory of three-dimensional consolidation // Journal of Applied Physics. - 1941. - Vol. 12, № 2. - P. 155-164. doi: 10.1063/1.1712886
- Finite element analysis of Biot’s consolidation with a coupled nonlinear flow model / Deng Yue-bao, Liu Gan-bin, Zheng Rong-yue, Xie Kang-he // Mathematical Problems in Engineering. - 2016. - Vol. 2016. - P. 1-13. doi: 10.1155/2016/3047213
- Biot M.A. General solutions of the equations of elasticity and consolidation for a porous material // Journal of Applied Mechanics. Trans. ASME. - 1956. - Vol. 23, № 1. - P. 91-96.
- Small J.C., Booker J.R., Davis E.H. Elasto-plastic consolidation of soils // International Journal of Solids and Structures. - 1976. - Vol. 12. - P. 431-448.
- Ferronato M., Castelletto N., Gambolati G. A fully coupled 3-D mixed finite element model of Biot consolidation // Journal of Computational Physics. - 2010. - Vol. 229. - P. 4813-4830. doi: 10.1016/j.jcp.2010.03.018
- Borja R.I., Alarcón E. A mathematical framework for finite strain elastoplastic consolidation. Part 1: Balance laws, variational formulation, and linearization // Computer Methods in Applied Mechanics and Engineering. - 1995. - Vol. 122. - P. 145-171. doi: 10.1016/0045-7825 (94) 00720-8
- Borja R.I., Tamagnini C., Alarcón E. Elastoplastic consolidation at finite strain. Part 2: Finite element implementation and numerical examples // Computer Methods in Applied Mechanics and Engineering. - 1998. - Vol. 159. - P. 103-122. doi: 10.1016/S0045-7825 (98) 80105-9
- Liu Z., Liu R. A fully implicit and consistent finite element framework for modeling reservoir compaction with large deformation and nonlinear flow model. Part I: Theory and formulation // Computational Geosciences. - 2018. - Vol. 22, iss. 3. - P. 623-637. doi: 10.1007/s10596-017-9715-3
- Liu Z., Liu R. A fully implicit and consistent finite element framework for modeling reservoir compaction with large deformation and nonlinear flow model. Part II: verification and numerical example // Computational Geosciences. - 2018. - Vol. 22, iss. 3. - P. 639-656. doi: 10.1007/s10596-017-9716-2
- Carter J.P., Booker J.R., Davis E.H. Finite deformation of an elasto-plastic soil // International Journal for Numerical and Analytical Methods in Geomechanics. - 1977. - Vol. 1. - P. 25-43.
- Mojarad R.S., Settari A. New solution for anisotropic formation damage due to produced water re-injection // Journal of Canadian Petroleum Technology. - 2009. - Vol. 48, iss. 4. - P. 1-7. doi: 10.2118/2007-113
- Settari A., Mourits F.M. A coupled reservoir and geomechanical simulation system // SPE Journal. - 1998. - Vol. 3, iss. 3. - P. 219-226.
- Settari A., Walters D.A. Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction // SPE Journal. - 2001. - Vol. 6, iss. 3. - P. 334-342. doi: 10.2118/74142-PA
- A comparison of techniques for coupling porous flow and geomechanics / R.H. Dean, X. Gai, C.M. Stone, S.E. Minkoff // SPE Journal. - 2006. - Vol. 11, iss. 1. - P. 132-140.
- Accelerating the convergence of coupled geomechanical-reservoir simulations / L. Jeannin, M. Mainguy, R. Masson, S. Vidal-Gilbert // International Journal for Numerical and Analytical Methods in Geomechanics. - 2006. - Vol. 31, iss. 10. - P. 1163-1181. doi: 10.1002/nag.576
- Wheeler M.F., Gai X. Iteratively coupled mixed and Galerkin finite element methods for poro-elasticity // Numerical Methods for Partial Differential Equations. - 2007. - Vol. 23, no. 4. - P. 785-797. doi: 10.1002/num.20258
- Kim J., Tchelepi H.A., Juanes R. Stability, accuracy, and efficiency of sequential methods for coupled flow and geomechanics. - SPE Journal. - 2011. - Vol. 16. - P. 249-262.
- Bergamaschi L., Ferronato M., Gambolati G. Mixed constraint preconditioners for the iterative solution to FE coupled consolidation equations // Journal of Computational Physics. - 2008. - Vol. 227. - P. 9885-9897. doi: 10.1016/j.jcp.2008.08.002
- Ferronato M., Gambolati G., Teatini P. Ill-conditioning of finite element poroelasticity equations // International Journal of Solids and Structures. - 2001. - Vol. 38. - P. 5995-6014. doi: 10.1016/S0020-7683 (00) 00352-8
- Reed M.B. An investigation of numerical errors in the analysis of consolidation by finite elements // International Journal for Numerical and Analytical Methods in Geomechanics. - 1984. - Vol. 8. - P. 243-257.
- Di Y., Sato T.Computational modelling of large deformation of saturated soils using an ALE finite element method // Annuals of Disas. Prev. Res. Inst., Kyoto University. - 2004. - No. 47 C. - P. 1-11.
- El-Amrani M., Seaid M. Eulerian-Lagrangian time-stepping methods for convection-dominated problems // International Journal of Computer Mathematics. - 2008. - Vol. 85. - P. 421-439. doi: 10.1080/00207160701210109
- Артамонова Н.Б., Шешенин С.В. Связанная задача консолидации в нелинейной постановке. Теория и метод решения // Механика композиционных материалов и конструкций. - 2020. - Т. 26, № 1. - С. 122-138. doi: 10.33113/mkmk.ras.2020.26.01.122_138.08
- Шешенин С.В., Артамонова Н.Б. Моделирование нелинейной связанной задачи консолидации // Механика композиционных материалов и конструкций. - 2020. - Т. 26, № 3. - С. 341-361. doi: 10.33113/mkmk.ras.2020.26.03.341_361.04
- Артамонова Н.Б. Численная реализация модели Био при больших деформациях: дис. … канд. физ.-мат. наук. - М.: МГУ, 2020. - 157 с.
- Щелкачев В.Н. Основы и приложения теории неустановившейся фильтрации. Ч. 2. - М.: Нефть и газ, 1995. - 493 с.
- Donea J., Huerta A. Finite element methods for flow problems. - Chichester: John Wiley & Sons Ltd, 2003. - 358 p.
- Киселев Ф.Б., Шешенин С.В. Разностная схема для задачи нестационарной фильтрации в слоистых грунтах // Известия РАН. Механика твердого тела. - 1996. - № 4. - С. 129-135.
- Шешенин С.В., Какушев Э.Р., Артамонова Н.Б. Моделирование нестационарной фильтрации, вызванной разработкой месторождений // Вестник МГУ. Сер. 1. Математика. Механика. - 2011. - № 5. - С. 66-68.
- Быченков Ю.В., Чижонков Е.В. Итерационные методы решения седловых задач. - М.: БИНОМ. Лаборатория знаний, 2014. - 349 с.
- Brezzi F., Fortin M. Mixed and hybrid finite element methods. - New York: Springer-Verlag, 1991. - 223 p.
- Elman H.C., Silvester D.J., Wathen A.J. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Numerical mathematics and scientific computation. - New York: Oxford University Press, 2005. - 400 p.
- Какушев Э.Р., Шешенин С.В., Закалюкина И.М. Итерационные методы решения связанной задачи фильтрации // Вестник МГСУ. - 2012. - № 9. - С. 129-136.
- Какушев Э.Р. Численное решение связанных трехмерных краевых задач упругой пористой среды: дис. … канд. физ.-мат. наук. - М.: МГУ, 2013. - 111 с.
- Григорян С.С. Об основных представлениях динамики грунтов // Прикладная математика и механика. - 1960. - Т. 24, № 6. - С. 1057-1072.
- Simo J.C., Hughes T.J.R.Computational inelasticity. - New York: Springer-Verlag, 1998. - 392 p.
- Bonet J., Wood R.D. Nonlinear continuum mechanics for finite element analysis. - New York: Cambridge University Press, 2008. - 340 p.
- Bonet J., Gil A.J., Wood R.D. Worked examples in nonlinear continuum mechanics for finite element analysis. - New York: Cambridge University Press, 2012. - 138 p.
- Новацкий В.К. Волновые задачи теории пластичности. - М.: Мир, 1978. - 307 с.
- Mróz Z., Zienkiewicz O. Uniform formulation of constitutive equations for clays and sands // Mechanics of Engineering Materials / eds. C.S. Desai, R.H. Gallaher. - New York: Wiley, 1984. - P. 415-449.
- Гидрогеология. Часть 1 / А.Я. Гаев [и др.] - Оренбург: ОГУ, 2010. - 75 с.
- Марышев Б.С. О горизонтальной напорной фильтрации смеси через пористую среду с учетом закупорки // Вестник Пермского университета. - Физика. - 2016. - № 3 (34). - С. 12-21. doi: 10.17072/1994-3598-2016-3-12-21
- Конструирование и расчет нежестких дорожных одежд / под ред. Н.Н. Иванова. - М.: Транспорт, 1973. - 328 с.