CC BY
33
Использование метода граничных интегральных уравнений для определения упругих модулей гранулированных геологических сред
Е.Б. Сибиряков, Е.В. Деев
Институт нефтегазовой геологии и геофизики СО РАН, Новосибирск, 630090, Россия
Работа посвящена развитию метода граничных интегральных уравнений теории упругости и возможностям применения этого метода для нахождения эффективных параметров гранулированных геологических сред. В качестве тензора Грина впервые используется модифицированный тензор фундаментальных решений третьего рода, пригодный для использования не только для выпуклых, но и для вогнутых поверхностей. При этом тензор имеет достаточно простой вид и удобен для решения задач, так как позволяет сводить краевую задачу не к сингулярным, а к регулярным интегральным уравнениям. Полученные зависимости позволяют анализировать возможные сценарии развития деформации (хрупкой или пластической) гранулированной среды в зависимости от типа порового флюида, удельной поверхности порового пространства, пористости и порового давления. Показано, что в гранулированной среде, заполненной газом, разрушение будет хрупким с образованием системы трещин. В случае заполнения пор жидкостью будет происходить флюидизация и разжижение среды.
Boundary integral equation method as applied to the determination of elastic moduli of granular geological media
E.B. Sibiryakov and E.V Deev
Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia
In the paper we develop a boundary integral equation method of elastic theory and consider the opportunities of its application to finding the effective parameters of granular geological media. For the first time a modified third-rank tensor of fundamental solutions, which is applicable to both convex and concave surfaces, is used as the Green tensor. The tensor has a simple form and is convenient for problem solving as it allows reducing a boundary problem to regular, rather than singular, integral equations. Based on the found dependences we can analyze possible evolution scenarios of deformation (brittle or plastic) in the granular medium depending on the pore fluid type, specific surface of the pore space, porosity and pore pressure. It is shown that in a gas-filled granular medium the fracture is brittle and occurs with the formation of a system of cracks. If the pores are filled with fluid, the medium fluidizes and liquefies.
1. Введение
Проблема описания параметров и понимания механизмов деформационного поведения гранулированных (зернистых) сред под воздействием внешних статических и динамических нагрузок имеет в геологии как фундаментальный, так и прикладной аспект. Попадая в сферу хозяйственной деятельности человека, зернистые осадки в терминах инженерной геологии рассматриваются в качестве соответствующих грунтов. Особенности физико-механического поведения при сейсмических толчках таких рыхлых грунтов, содержащих жидкий и газовый (в различных сочетаниях) поровый флюид, проявляются в формировании, наряду с разрывными и сейсмогравитационными структурами, специфического класса поверхностных проявлений: излияний (иногда с фонтанированием) водных масс и осадков различной
зернистости из трещин, песчаных вулканов, грифонов, просадок, оплывневых дислокаций. Подобные структуры согласно обобщениям [1] являются характерным атрибутом 8-9-балльных и более сильных землетрясений. Как правило, они возникают на участках с неглубоким залеганием грунтовых вод и преимущественно характерны для осадков от илистой до песчаной размерности.
Инженерные сооружения, расположенные на влагонасыщенных рыхлых грунтах, испытывают существенное воздействие, приводящее к наклону и опрокидыванию зданий, всплыванию трубопроводов, разжижению намывных насыпей и т.д. Существующие оценки показывают, что на водонасыщенных рыхлых и средней плотности песчаных грунтах можно ожидать дополнительного увеличения интенсивности сейсмического воз-
© Сибиряков Е.Б., Деев Е.В., 2008
действия за счет изменения их физико-механических свойств примерно на 1 балл [2].
В то же время, деформационные структуры, связанные с сейсмическими событиями, регулярно выделяются при изучении осадочных толщ и несут важную информацию о палеогеодинамических параметрах территорий. Начиная с работы [3], они выделяются в особый класс, называемый «сейсмиты». К настоящему моменту сейсмиты выявлены в генетически разнотипных континентальных, прибрежно-морских и морских толщах (например [4-11]). Несмотря на то что деформационные структуры этого вида наиболее хорошо изучены в разрезах кайнозойских осадков, в особенности плей-стоцен-голоценовых, морфологически схожие с ними образования обнаружены и в более древних осадочных толщах, начиная с протерозоя (например [12-16]). Лабораторные эксперименты демонстрируют принципиальную возможность получения морфологически подобных структур при вибрационном воздействии на рыхлые осадки в различных режимах нагружения [17, 18].
Но параметры внутреннего напряженного состояния природных гранулированных геологических сред, определяющие их деформационное поведение, в геологии описываются лишь на качественном уровне. Для объяснения процессов формирования разных по морфологии деформационных структур предложены механизмы гидропластичности, флюидизации и разжижения [19-21], различие между которыми базируется на соотношении внутреннего (порового) и внешнего (литостатического) давлений, действующих на слабосцементированный влагонасыщенный осадок.
Представляется, что подобные механизмы также имеют существенное значение и для понимания процессов структурирования неконсолидированных и слабо-консолидированных зернистых сред в прибрежных зонах под действием волновых колебаний водных масс, в том числе штормового и цунамигенного характера [22], а также для выяснения механического поведения зернистых осадков в процессах их переотложения в водных условиях (обломочные, зернистые, флюидизированные и турбидитные потоки).
В настоящей работе предпринята попытка оценки характера деформационного поведения (механизмов деформации) гранулированной флюидонасыщенной среды в зависимости от упругих модулей скелета и структуры порового пространства.
2. Характеристики исследуемой модели гранулированной среды
Составляющие алевриты и пески частицы представляют собой многогранники, приближающиеся с усилением степени окатанности к сферической форме. Поэтому идеализированная модель структуры их минерального скелета может быть представлена в виде набора шаров.
В природе гранулированные среды характеризуются различной степенью заполнения порового пространства водой (степенью водонасыщения): от единицы, соответствующей состоянию полной влагоемкости, т.е. полному заполнению порового пространства водой, до нуля, когда поровое пространство заполнено газовой фазой, в том числе парами воды. В дальнейшем исследуются эти два крайних состояния. Кроме того, в рассматриваемой модели вся свободная вода в поровом пространстве рассматривается как гравитационная, т.е. не связанная с поверхностью частиц и не удерживаемая капиллярными и менисковыми силами. Это объясняется тем, что жидкая влага в песчаных грунтах находится преимущественно в виде свободной воды, при незначительных количествах физически связанной. По существующим оценкам максимальная молекулярная влагоем-кость песков, т.е. их способность удерживать в своей структуре физически связанную воду, не превышает 2 % [23].
С использованием выбранной модели исследованию подлежали вопросы зависимости упругих модулей среды от порового давления, структуры порового пространства, упругих модулей скелета.
Математически среда смоделирована зернами сферической формы с радиусом равным 1. Координатная сетка разбита по углам 0 и ф на 51x51 равных частей. При этом среднее число контактов равно 8, т.е. имеется 8 равных площадок контакта зерна с другими зернами, центры которых образуют вписанный в шар куб и площадь каждой из которых равна 1/16 части от площади сферы. Оставшаяся часть сферы граничит с флюидом. Осреднение производится путем поворота площадок контакта на определенный телесный угол и вычислением среднего арифметического.
3. Перспективы и трудности в использовании метода граничных интегральных уравнений для вычисления средних упругих модулей среды
Очевидно, что гранулированные геологические среды относятся к разновидности контрастных микронеод-нородных сред, у которых перепад упругих свойств между скелетом и флюидом составляет многие порядки. Количественное и качественное описание воздействия возмущения на подобную среду требует адекватного учета граничных условий процессов деформирования на всей поверхности раздела скелет-флюид. Попытки игнорировать это обстоятельство и заменить контрастную микронеоднородную среду некоторым эквивалентным сплошным телом приводят к трудностям теоретического объяснения ряда экспериментальных фактов (например роста коэффициента Пуассона с увеличением дифференциального давления в газонасыщенных средах и падения в насыщенных жидкостью [24], вариаций коэффициента Пуассона при компактировании ме-
таллических гранул в зависимости от внешнего давления [25]).
Задача вычисления перемещений на достаточно произвольной замкнутой поверхности в зависимости от заданных на границе нагрузок или нагрузок по заданным на границе перемещениям, либо напряженного состояния в случае смешанной задачи, когда на одной части поверхности заданы нагрузки, а на другой — перемещения, может быть применена для расчета средних упругих характеристик гранулированной среды путем замены разностных операторов дифференциальными [26].
Если известны структура минерального скелета (упаковка зерен) и законы взаимодействия на контактах, то путем осреднения сил, действующих на зерно, возможно вычислить средние упругие модули среды. Однако случаев, когда такие законы известны, — единицы (например, закон Г ерца — типичная задача смешанного типа). Более того, если среда является многофазной, то совершенно неясно, как учитывать взаимодействие фаз (главный недостаток теории смесей), даже если закон взаимодействия зерен на контактах известен. Таким образом, для нахождения эффективных упругих модулей необходимо решать упругую статическую задачу.
Метод граничных интегральных уравнений предназначен для того, чтобы сводить трехмерную задачу к решению двухмерных интегральных уравнений. Его суть состоит в том, что перемещения в любой точке объема (в том числе на поверхности) ищутся в виде свертки тензора Грина для пространства с некоторым потенциалом [27]:
Ui(x) = ¿J Tik (Х’ У)Fk (y)dSy
(1)
где по значку к ведется суммирование; точка х заключается в объеме и может стремиться сколь угодно близко к поверхности изнутри (или извне); точка у бежит по поверхности, в ней же берется элемент поверхности. При этом
г**(у) = - 1
2ц(к + 2ц)
Эг дr
(к + ц) —— -— + (А + 3M.)Sft дxi axr.
1
r(x, y)
(2)
где r — расстояние между фиксированной точкой и точкой, бегущей по поверхности; — символ Кроне-кера; X и ц — упругие константы Ламе. Если точка х стремится к поверхности изнутри, то в ней можно вычислить нагрузки:
Px (Гл)(х, у) =
= + Xn0divх (гй) + ц(п0 Xrotх (Tik)) =
дп
0
к + 2ц
д— д—
цй* + 3(к + ц) ^ дг
dxi дхк
д
дп„
1
+
r
+ц
cos(n0 xk)
дxi
1 ( \ д -- cos(no x)— r дxk
(3)
где п0 — единичная нормаль в точке х. В результате, если на поверхности заданы нагрузки, то для нахождения потенциала в трехмерном случае (с помощью которого можно будет найти перемещения в любой точке объема), получится двухмерное сингулярное интегральное уравнение (т.е. уравнение, которое содержит неин-тегрируемую особенность). Интеграл в этом случае нужно понимать в смысле главного значения (из-за последнего члена в (3)). Численное решение подобного рода уравнений достаточно затруднительно и представляет отдельную проблему. По этой причине метод граничных интегральных уравнений использовался только для решения двухмерных задач, хотя преимущества этого метода для решения краевых задач понимают многие исследователи [28].
Использование тензора фундаментальных решений третьего рода, который с точностью до множителя совпадает с тензором Грина для полупространства (если поверхность — плоскость) и тензор нагрузок от которого не содержит сингулярных членов, представляется более предпочтительным [29]. Однако данный тензор плохо применим для решения задач на вогнутых поверхностях (т.е. поверхностях с отрицательной кривизной в некоторых точках), а также на поверхностях многосвязных областей. Оказалось, что эту трудность можно преодолеть. Кроме того, оказалось возможным записать в достаточно простой и доступной форме декартовы компоненты тензора фундаментальных решений третьего рода, что для целого ряда задач является значительно более удобным.
4. Совершенствование тензора фундаментальных решений третьего рода в теории потенциала
Как известно, тензор фундаментальных решений третьего рода записывается в виде [30]:
М-к (x, y) =
1
к + ц
(к + 2ц)Г-к ^, у) - 2 Z-k (Л у)
(4)
где точка у помещается в начало координат; направление оси х1 совпадает с направлением внутренней нормали к поверхности в точке у (соответственно х = (г> п)); точка х имеет координаты (х1, х2, х3);
Zik =
д 2 g д 2 g д 2 g
дx12 дx1дx2 дx1дx3
д 2 g д 2 g д^ д 2 g
дx1дx2 дx12 дx3 дx2дx3
д 2 g д 2 g д 2 g д^
дx1дx3 дx2дx3 дx12 дx2
g = x ln(r + x) - —• При этом
(5)
РХ (Мік) = 3
дх дх дх2
1 У 1 2
дг дг { дг 1 дх1 дх2
дх2
дг дг дг дг
дх1 дх3 дх2 дхз
дг дг дх1 дх3
дг дг дх2 дх3
Г дг 12
дх3
д 1
дпо г
А + Ц
(п0,ет1)
0
д3 g
(п0,ет 2)
д3g (п е )
дх1дх| 0’ Т2 дх1дх2дх3
g ^ д3 g
д3'
дх1дх2
--(п0> ет1)
дх1дх2дх3
(п0,ет1)
дх2дх2
д3 -
- (п0>ет2)
(п0,ет2)Т-^ (п0,ет1^ 2_.
дх2 дх2 дх
, ет1)-
дх^дх3 д3 g
(п0> ет1)
д 3 g
л0> - (п
(п0> ет 2)
дх3 д3 £
0> ет2)
дх2дх3
- (п0> ет1)
_д^
дх2дх2 д 3 g
дх2д х
• (7)
Т.е. в уравнениях (1) и (3) для нахождения потенциала и решения задачи нужно использовать Мл и Рх (М1к) вместо Ггк и Рх (Ггк). Далее Рх (Мк) обозначается как рл. При этом уравнения уже не будут сингулярными (т.е. все особенности — интегрируемые типа &5/г). Нужно учитывать, что х1, х2, х3 — не произвольные координаты, а проекции радиус-вектора, проведенного из бегущей точки в фиксированную, на нормаль и две произвольные взаимно ортогональные касательные (т.е. система координат отнюдь не произвольна, как утверждается в [27]).
Как видно, Р1к не содержит сингулярных особенностей и соответствующие интегральные уравнения будут регулярными. Запишем также Мл:
1
2ц
(8)
1 х2 + Ц 1 х1 х2 , Ц д ^
г г3 А + Ц г г3 ' А + Ц дх1дх2
х1х2 1 Ц д 1 . х2 1 Ц д
г3 А + ц дх1дх2 г ' г3 А + ц дх^
хіх3 Ц д 2g х2 х3 Ц д 2g
г3 А + Ц дх1дх3 г3 А + Ц дх2дх3
х1х3 Ц д2 g
х3 :+^г -
А + ц дх1дх3
Ц д2 g
А + ц дх2дх3
Ц д 2 g
А + Ц дх32
В принципе, если за направление оси х2(ет1) выбрать вектор г - (г, п)п (г — радиус-вектор), а за направление х3(ет2) — п X г, то х3 будет равным нулю и в тензорах
(7) и (8) будет не 8 независимых компонент, а только 5.
Тензор (8) и тензор нагрузок (7) можно использовать лишь в том случае, если скалярное произведение радиус-вектора г (проведенного из бегущей точки в фиксированную) на вектор внешней нормали к поверхности п неотрицательно (х1 > 0), т.е. например, если поверхность всюду выпуклая. В самом деле, если х1 < 0, то в
точке, расстояние до которой г = -х1, появится ничем не обоснованная интегрируемая особенность. Т.е. несмотря на то что формально тензор Мл (8) удовлетворяет системе уравнений упругости и его можно использовать для решения краевых задач, сходимость метода последовательных приближений при решении интегрального уравнения типа Фредгольма второго рода (для нахождения вектора потенциала) будет значительно хуже. Что же делать, если поверхность (или ее часть) вогнутая? Физический смысл функции g таков, что в этом случае нужно заменить х1 на - х1, а затем попытаться подправить решение, если столбцы матрицы не будут удовлетворять системе уравнений упругости.
Дело в том, что столбцы матрицы (5) удовлетворяют системе уравнений упругости не только при g = = х11п(г + х1) - г, но и при g = -х11п(г - х1) - г. То есть фактически найден еще один тензор фундаментальных решений, но уже для вогнутых поверхностей. Значит тензор фундаментальных решений третьего рода Мл
(8) можно использовать для решения краевых задач, как для выпуклых, так и для вогнутых областей, если
И =1 х1|1п(г+1 х1|) - г. (9)
Используя Мк (8) и рк (7) можно строить систему регулярных интегральных уравнений, находить потенциалы Fk и решать соответствующие краевые задачи. Однако при данном построении тензоров фундаментальных решений и нагрузок от них следует принять во внимание, что координаты х1, х2, х3 заданы в подвижном базисе. При этом декартовы проекции тензоров (8), использование которых в целом ряде задач представляется более удобным, казалось бы, будут иметь достаточно сложный вид. Однако их можно записать в достаточно простом и понятном виде.
5. Использование декартовых проекций тензоров М1к и р для построения интегральных уравнений
Главное удобство формы записи тензора (2) в том, что не требуется выбирать единичные касательные ортогональные векторы и вычислять х2 и х3 в каждой точке. Зачастую нагрузки или перемещения заданы в неподвижном базисе и компоненты тензоров ((8), а особенно (7), где необходимо вычислять третьи производные) достаточно громоздки. Попробуем записать их в действительно произвольной системе координат. Для удобства направим радиус-вектор из фиксированной точки в бегущую, а за положительное направление нормали возьмем внешнюю нормаль.
Тензор Мл — это отклик перемещений в направлении i на 8-нагрузку, приложенную в направлении к. Сначала вычислим вектор Мк, то есть отклик в направлениях п, ет1, ет2 на 8-нагрузку, приложенную в направлении к (к — направление декартовой системы координат, п, ет1, ет2 — направления единичной нормали и двух ортогональных друг другу касательных векторов в бегущей точке поверхности):
Мк = п[мп(к,п) + М12(к, еТ1) + М1з(к, ет2)]+ + ет1 [М21(к, п) + М22 (к, ет1) + М23 (к, ет2)]+ + ет2 [М31(к, п) + М32(к, ет1) + М33(к, ет2)], (10)
соответственно 1
2ц[ г г 1 д
-3ГГк (г>по) + Ц д
Мк
3
2(Х + ц) дгк где г- и гк
±П 1п(г ± х1) +
Г - хіП
(11)
■ проекции радиус-вектора на произвольные декартовы оси; производная —— берется при усло-
дп 0 Т й
вии, что —L = 0. Теперь можно найти декартовы ком-
дГк
поненты тензора Рк (7):
р = 3 ГГк _д
рк 3 2 3
г дп
11 Ц д - +—----------:
о |_ г ] Х + цдгк
х[(п0^га4?3 Xп)(еті,і) -- ХІ)] =
Х + цдГк
П (п,по) - по- + (Г п,по)[г,п].
Г + К
г (г+ІХ11)
(12)
где (г,п,п0) — смешанное произведение трех векторов, а [г,п] — векторное произведение. Предлагается в дальнейшем именно (11) называть тензором фундаментальных решений третьего рода, а (12) — тензором нагрузок для него.
Для решения интегральных уравнений с интегрируемой особенностью зачастую оказывается полезным аддитивный способ избавления от особенности, который заключается в том, что из исходного интеграла (1), а также из системы интегральных уравнений вычитается значение вектора-потенциала в фиксированной точке х, а затем прибавляется. Т.е. вместо (1) при вычислениях следует использовать (13):
и (х) = 2^1 Мік (х, у) ^ (у) - Fk (х)>^ +
+ I Мк (х, у^у,
(13)
аналогично и для системы интегральных уравнений. Но для использования аддитивного метода необходимо вычислить интегралы от тензоров (11), (12). В частности, упомянутые интегралы по сфере радиуса R можно вычислить, используя преобразование координат, упомянутое в [30]. А именно, начало координат помещается в фиксированную точку, стремящуюся к поверхности изнутри, а оси координат ориентируются так, чтобы внутренняя нормаль в этой точке имела координаты (0, 0, 1). Запишем координаты точки у в сферической системе координат с переменным радиусом г и центром в точке х, обозначив ^ полярный, а ф азимутальный угол [30]. При этом
ег = е х1 sm ^ cos ф + е х2 sm ^ sm ф + е х3 cos ^, е^ = е х1 cos ^ cos ф + е х 2 cos ^ sin ф - е х 3 sin ^,
еф = -е х1 ^п ф +е х 2С^ ф,
(14)
а ех1, ех2, ех3 выражаются через единичные орты декартовой системы координат ^ j, к и обычные сферические координаты точки х:
е
(N1 * е =
3 * е
1 -cos ф0(1 + cos0о) -sinф0cosф0(1 + cos0о) sin0о cosф0
- sin ф0 cos ф0 (1 + cos 0о) - sin 0о cos ф0
1 - sin ф0(1 + cos0о) - sin 0о sin ф0
sin 0о sin ф0
Рис. 1. Компонента х вектора нагрузки на сфере
Рис. 3. Разность между нулевым и 15 приближением х-компоненты вектора потенциала
В этом случае расстояние между точками х и у будет r = 2^cosy, а внешняя нормаль в точке у и элемент поверхности будут равны соответственно:
ny = er cos y + ey sin y, (16)
d£y = 2 Rr sin y dy dcp. (17)
Заметим, что для сферы (r, n, n0) = 0, (r, n) = = -(r, n o) = Xi = r2 /2 R. Результаты вычислений:
1- J M л (x, y )d 5y =
2 R5;¡,
2 n
^JPik (x, y)dSy 2n
V
2 1 —I--------
ц А + ц
У
(18)
(19)
6. Вычисление упругих модулей гранулированной среды
Естественно, что при решении в качестве теста задачи о всестороннем сжатии шара с использованием тензоров (11), (12) получается правильный ответ. Хотелось бы, однако, исследовать скорость сходимости при наличии не только нормальных, но и касательных нагрузок, а заодно и опробовать методику вычисления упругих модулей. Итак, зададим в среде с X = ц = 1 постоянное напряжение а хх = 1. Вырежем в среде сферу радиуса Я = 1. Зная напряжения в каждой точке границы, можно задать нагрузки. Задав нагрузки, найдем вектор потенциала, используя систему уравнений второго рода:
Pi(x) = Fi(х) - т1 J Pik(X У) Fk(y)d Sy >
2П
(20)
где рі = оікпк — вектор нагрузок. Решать уравнение будем методом последовательных приближений. На рис. 1 изображена х-компонента вектора нагрузки (нулевое приближение для той же компоненты вектора потенциала), на рис. 2 — та же компонента вектора потенциала после 15 приближений, на рис. 3 — разность между ними. Как видно, скорость сходимости достаточно хорошая (т.е. уже нулевое приближение достаточно близко к точному решению). Далее, зная вектор потенциала, вычислим перемещения на границе (13). На рис. 4 изображена та же компонента вектора перемещений. Зная перемещения на поверхности, вычислим средние деформации в объеме с помощью теоремы о градиентах:
I gradФd V = I Фпё£,
V Б
соответственно можно вычислить все средние значения по объему от производных, например
j|^ dV = JO«*dS.
V dx S
(21)
После этого, зная нагрузки на поверхности, найдем средние в объеме напряжения. Зная напряжения и деформации, вычислим упругие модули. Получаем, что при разбиении поверхности сферы по углам 0 и ф на 51x51 частей погрешность вычисления упругих модулей составляет менее 0.3 %. Отсюда можно сделать вывод, что с помощью тензора фундаментальных решений третьего рода можно успешно решать краевые задачи теории упругости.
Рис. 2. Компонента х вектора потенциала после 15 приближений
Рис. 4. Компонента х вектора перемещений
Главное отличие микронеоднородной среды от рассмотренного выше теста будет в том, что часть поверхности сферы будет контактировать с флюидом (т.е. нормальная компонента вектора нагрузок должна будет совпадать с давлением во флюиде), а оставшаяся часть будет находиться в контакте со скелетом. При этом следует принять во внимание, что в изолированных порах давление может быть каким угодно и совершенно не зависеть от напряженного состояния скелета, но если поры сообщаются, то связь между давлением в скелете и флюиде становится хотя и сложной, но вполне определенной. Эта связь определяется, прежде всего, удельной поверхностью, а также величиной и знаком кривизны поверхности раздела скелет-флюид. Трещиноватые и зернистые среды обладают разной структурой порового пространства. Если для зернистых сред характерны постоянный (отрицательный) знак кривизны порового пространства, а также близкие значения средней и гауссовой кривизн, то для трещиноватых коллекторов характерно резкое изменение кривизны от нулевых значений до исключительно высоких (как положительных, так и отрицательных). Кроме того, средние и гауссовы кривизны в трещиноватых телах могут существенно отличаться друг от друга. Даже при элементарно простом среднем напряженном состоянии, которое определяется весом вышележащих горных пород и коэффициентом бокового распора (или отношением у,/Ур), соотношение давлений в скелете и флюиде существенным образом зависит от удельной поверхности, а также от знака и величин средней и гауссовой кривизн.
В данной работе ограничимся моделированием зернистых сред (для моделирования трещиноватых сред необходимо использовать многосвязную область). Будем считать, что элементарная представительная микроструктура подвержена действию таких нагрузок, которые вызывают одинаковые деформации твердой фазы и флюида на границах микроструктуры. При этом, естественно, соответствующие нагрузки оказываются различными, что соответствует простому напряженному состоянию плоскопараллельной среды.
Предпринятые ранее попытки оценки упругих модулей [31], основанные на методах механики композитов и теории осреднения дифференциальных операторов для контрастных микронеоднородных сред, не учитывают межфазного взаимодействия, которое заключается в том, что скелет и флюид не могут свободно проникать друг в друга. Объем, равный интегралу |UndS по поверхности раздела скелет-флюид, деленный на объем порового пространства, как раз и является мерой меж-фазного взаимодействия, количественно показывающей разность деформаций в скелете и флюиде при одинаковой деформации, приложенной к внешней поверхности. С помощью подобного подхода можно решать две задачи. Во-первых, можно рассматривать зависи-
мость эффективных упругих модулей от структуры по-рового пространства, рассматривая поровое давление как независимый параметр. Во-вторых, рассматривая поровое давление как функцию от напряженного состояния среды (т.е. при отсутствии внешней нагрузки поровое давление равно нулю), можно, определив деформацию в скелете и флюиде, определить зависимость порового давления от структуры порового пространства и упругих модулей среды.
Будем моделировать среду следующим образом. Центральный шар радиуса 1 окружен идентичными шарами таким образом, что центры масс соседних шаров совпадают с точками пересечения ребер некоторого куба. Этот куб и возьмем в качестве представительного объема, так как этот объем в среде будет повторяться. Несмотря на то что моделируемая среда будет анизотропной, вычислить эффективные упругие модули изотропной среды можно с помощью осреднения по телесному углу поворота куба. Таким образом, в представительном объеме будет один целый шар и восемь октантов.
Пористость и удельная поверхность в таком объеме будут связаны. Например, в случае точечных контактов объем твердой фазы будет равен 8п/ 3, а объем среды — (4Д/3)3. Таким образом, пористость будет равна ~0.32. В случае если контакты не точечные, пористость будет существенно меньше. Например, если выбрать площадки контакта так, чтобы их суммарная площадь поверхности была равна половине площади поверхности сферы (расстояние между центрами масс соседних зерен станет 2 - 0.25 = 1.75), объем среды будет равен ((4 - 0.5)Д/3)3, а объем твердой фазы — 8п/ 3 -я/ 4(1 -1/24). Пористость соответственно равна только 0.08. Если приложить к поверхности данной среды деформацию ехх = 1, то поровое пространство и скелет будут деформироваться по-разному за счет того, что на площадках контакта скелет-флюид будут другие граничные условия. Если вычислить компоненты деформации скелета и флюида, а также средние значения напряжений в объеме, то можно будет вычислить и упругие модули среды.
Средние напряжения в объеме вычисляются через нагрузки на поверхности представительного объема. Средние по объему деформации вычисляются через деформации на поверхности, а также через средние деформации скелета и флюида. Средние параметры деформации зерна можно вычислить с помощью процедуры, изложенной выше, проводя вычисления для отдельного зерна. Если заменить плоские площадки контакта на выпуклые сегменты, то по принципу Сен-Венана нагрузки на сегментах будут такими же, что и на плоских площадках контакта. Таким образом, решение задачи о нахождении упругих модулей сводится к решению задачи теории упругости на отдельном зерне. На рис. 5 изображена х-компонента вектора перемещений на единичной сфере с 8 площадками контакта скелет-скелет,
Рис. 5. Компонента х вектора перемещений на единичной сфере с 8 площадками контакта скелет-скелет, деформация на контакте скелет-скелет е хх = 1, поровое давление р0 = 0
площадь каждой из которых равна 1/16 части площади поверхности сферы, деформация на контакте скелет-скелет ехх = 1, поровое давление р0 = 0. На рис. 6, 7 изображена зависимость эффективного ц и X + 2ц от порового давления при тех же условиях, на рис. 8 — зависимость отношения скоростей поперечных и продольных волн у = У;/Ур от порового давления. Отрицательное значение порового давления может иметь место при наличии капиллярных сил. Видно, что с ростом порового давления растет как X, так и ц. Однако параметр у уменьшается с ростом порового давления. Следует отметить, что параметр у зависит, вообще говоря, не от типа флюида, а от давления во флюиде.
При наличии сильной внешней нагрузки, вызывающей разрушение, данный подход, естественно, не применим для точных расчетов разрушения, однако результаты, полученные выше, позволяют сделать качественный вывод о характере разрушения. Поскольку нагрузка
Рис. 6. Зависимость ц среды от порового давления
Рис. 7. Зависимость X + 2ц среды от порового давления
Рис. 8. Зависимость отношения у= Ур среды от порового дав-
нарастает за какой-то конечный промежуток времени, на начальном этапе разрушения уравнения упругости достаточно адекватно описывают процесс. Если параметр у уменьшается с ростом нагрузки, то будут происходить флюидизация и разжижение среды. Как правило, вода в порах находится под существенно большим давлением, чем воздух. Рост внешней нагрузки приводит к росту параметра |ип&Б, что еще более усиливает поровое давление. Если же поровое пространство заполнено воздухом, то поровое давление под влиянием этого фактора расти практически не будет (так как воздух значительно более сжимаем, чем вода). Соответственно характер разрушения, по-видимому, будет хрупким, с образованием системы трещин.
7. Выводы
Таким образом, трудности, препятствующие широкому применению метода граничных интегральных уравнений для решения трехмерных задач, следует считать устраненными. Впервые тензор фундаментальных решений третьего рода записан в неподвижной (декартовой) системе координат. Кроме того, этот тензор можно использовать для решения упругих задач в многосвязных областях, а также на телах, ограниченных вогнутыми поверхностями. Представляется разумным использовать его в дальнейшем для решения задач на многосвязной области, что позволит моделировать трещиноватые среды. Следует отметить высокую скорость сходимости метода последовательных приближений для системы (20), а также хорошую обусловленность матрицы системы.
Показана возможность использования метода граничных интегральных уравнений для построения эффективных параметров многофазных микронеоднород-ных сред, учитывающих межфазное взаимодействие. Приведены примеры расчета упругих модулей для зернистой среды.
Получено, что отношение у = У^/Ур падает с ростом порового давления.
Получила обоснование гипотеза о том, что два основных типа разрушений (образование системы трещин или флюидизация и разжижение среды), наблюдаемых в сейсмитах, связаны, главным образом, с типом флюида в поровом пространстве, а также с параметром меж-фазного взаимодействия | UndS.
В ближайшем будущем планируется использовать данный подход для нахождения эффективных упругих модулей трещиноватых сред.
Работа выполнена при поддержке РФФИ (грант № 06-05-64772-а), а также междисциплинарного проекта №18 СО РАН.
Литература
1. Попова Е.В. Остаточные деформации грунтов при землетрясениях
Ч. I // Вопросы инженерной сейсмологии. Вып. 16. - М.: Наука, 1974. - С. 209-251; Ч. II. - Вып. 17. - 1975. - С. 150-205; Ч. III. -Вып. 18. - 1976. - С. 155-193; Ч. IV. - Вып. 19. - 1978.- С. 171190; Ч. V. - Вып. 20. - 1980. - С. 190-203.
2. Ершов И.А., Попова Е.В. О влиянии обводненности грунтов на интенсивность сейсмического воздействия // Вопросы инженерной сейсмологии. Вып. 19. - М.: Наука, 1978. - С. 117-139.
3. Seilacher A. Fault-graded bends interpreted as seismites // Sedi-mentology. - 1969. - V. 13. - P. 155-159.
4. Hempton M.R., Dewey J.F. Earthquake-induced deformational structures in young lacustrine sediments, East Anatolian Fault, southeast Turkey // Tectonophysics. - 1983. - V. 98. - Nos. 3-4. - T7-T14.
5. Plaziat J.-C., PurserB.H., PhilobbosE. Seismic deformation structures
(seismites) in the syn-rift sediments of the NW Red Sea (Egypt) // Bull. Soc. Geol. France. Ser. 8. - 1990. - V. VI. - No. 3. - P. 419-434.
6. Obermeier S.F. Liquefaction evidence for strong earthquakes of Holo-cene and latest Pleistocene ages in the states of Indiana and Illinois, USA // Engineering Geology. - 1998. - V. 50. - P. 227-254.
7. Jones A.P., Omoto K. Towards establishing criteria for identifying trigger mechanisms for soft-sediment deformation: a case study of Late Pleistocene lacustrine sands and clays, Onikobe and Nakayama-daira Basins, northeastern Japan // Sedimentology. - 2000. - V. 47. -No. 6. - P. 1211-1226.
8. Moretti M. Soft-sediment deformation structures interpreted as seismites
in middle-late Pleistocene aeolian deposits (Apulian foreland, southern Italy) // Sedimentary Geology. - 2000. - V. 135. - P. 167-179.
9. Becker A., Ferry M., Monecke K. et al. Multiarchive paleoseismic record of late Pleistocene and Holocene strong earthquakes in Switzerland // Tectonophysics. - 2005. - V. 400. - P. 153-177.
10. González de Vallejo L.I., Tsigé M., Cabrera L. Paleoliquefaction features on Tenerife (Canary Islands) in Holocene sand deposits // Engineering Geology. - 2005. - V. 76. - Nos. 3^. - P. 179-190.
11. Деев Е.В., Гибшер А.С., Чигвинцева Л.А. и др. Микросейсмодислокации (сейсмиты) в плейстоценовых осадках Горного Алтая // ДАН. - 2005. - Т. 403. - № 1. - С. 71-74.
12. Guiraud M., Plaziat J.-C. Seismites in the fluviatile Bima
sandstones: identification of paleoseisms and discussion of their magnitudes in a Cretaceous synsedimentary strike-slip basin (Upper Benue, Nigeria)
// Tectonophysics. - 1993. - V. 225. - No. 4. - P. 493-522.
13. Rossetti D.F. Soft-sediment deformation structures in late Albian to Cenomanian deposits, Sao Luis Basin, northern Brazil: evidence for paleoseismicity // Sedimentology. - 1999. - V. 46. - No. 6. - P. 10651081.
14. Netoff D. Seismogenically induced fluidization of Jurassic erg sands, south-central Utah // Sedimentology. - 2002. - V 49. - No. 1. - P. 6580.
15. JewellH.E., Ettensohn F.R. An ancient seismite response to Taconian far-field forces: The Cane Run Bed, Upper Ordovician (Trenton) Lexington Limestone, central Kentucky (USA) // J. Geodynamics. -2004. - V. 37. - P. 487-511.
16. Mazumder R., van Loon A.J., Arima M. Soft-sediment deformation structures in the Earth’s oldest seismites // Sedimentary Geology. -2006. - V. 186. - Nos. 3-5. - P. 19-26.
17. Owen G. Experimental soft-sediment deformation: structures formed by the liquefaction of unconsolidated sands and some ancient examples // Sedimentology. - 1996. - V. 43. - No. 2. - P. 279-293.
18. Moretti M., Alfaro P., Caselles O., Canas J.A. Modelling seismites with a digital shaking table // Tectonophysics. - 1999. - V. 304. -No. 4. - P. 369-383.
19. Lowe D.R. Water escape structures in coarse-grained sediments // Sedimentology. - 1975. - V. 22. - No. 2. - P. 157-204.
20. Allen J.R.L. Sedimentary structures, their character and physical basis. - Amsterdam-New York: Elsevier, 1983. - 663 p.
21. Owen G. Deformation processes in unconsolidated sands // Deformation of sediments and sedimentary rocks / Ed. by M.E. Jones and R.M.F. Preston. - Oxford: Blackwell Scientific Publications, 1987. -P. 11-24.
22. Alfaro P., Delgado J., Estevez A. et al. Liquefaction and fluidization structures in Messinian storm deposits (Bajo Segura Basin, Betic Cordillera, southern Spain) // Int. J. Earth Sci. (Geol Rundsch). -2002.- V. 91. - No. 3. - P. 505-513.
23. ДенисовН.Я. Инженерная геология. - М.: Госстройиздат, 1960. -404 с.
24. Carcione J.M., Cavallini F. Poisson’s ratio at high pore pressure // Geophysical Prospecting. - 2002. - V. 50. - No. 1. - P. 97-106.
25. Carnavas PC., Page N.W. Elastic properties of compacted metal powders // J. Mater. Sci. - 1998. - V. 33. - No. 18. - P. 4647-4655.
26. Сибиряков Е.Б. Зависимость между коэффициентом Пуассона и микроструктурой в микронеоднородной среде // Физ. мезомех. -2004. - Т. 7. - № 1. - С. 63-68.
27. Купрадзе В.Д. Методы потенциала в теории упругости. - М.: Физматгиз, 1963. - 472 с.
28. Джонсон К. Механика контактного взаимодействия. - М.: Мир, 1989. - 510 с.
29. Партон В.З., Перлин П.И. Методы математической теории упругости. - М.: Наука, 1981. - 688 с.
30. Сибиряков Е.Б. Об использовании метода граничных интегральных уравнений для определения параметров микронеоднородных сред // Физ. мезомех. - 2006. - Т. 9. - № 1. - С. 97-101.
31. Шермергор Т.Д. Теория упругости микронеоднородных сред. -М.: Наука, 1977. - 400 с.
Поступила в редакцию 29.10.2007 г.
