CC BY
23
Об использовании метода граничных интегральных уравнений для определения параметров микронеоднородных сред
Е.Б. Сибиряков
Институт геофизики СО РАН, Новосибирск, 630090, Россия
Работа посвящена развитию метода граничных интегральных уравнений теории упругости. Тензор Грина для полупространства и тензор нагрузок от него записаны не в декартовых координатах, а в обобщенных координатах, связанных с поверхностью тела. При этом тензор имеет достаточно простой вид и удобен для решения задач, так как позволяет сводить краевую задачу не к сингулярным, а к регулярным интегральным уравнениям.
On a boundary integral equation method used to determine parameters of microheterogeneous media
E.B. Sibiryakov
Institute of Geophysics SB RAS, Novosibirsk, 630090, Russia
The paper develops a boundary integral equation method of elasticity theory. Green’s tensor for a half-space and load tensor are written in the generalized coordinates associated with the body surface rather than in the Cartesian coordinates. The tensor has a rather simple form and is convenient for solving problems as it allows reducing a boundary-value problem to regular integral rather than singular equations.
1. Введение
Контрастные микронеоднородные среды, т.е. такие, в которых перепад упругих свойств между матрицей и флюидом составляет многие порядки, требуют для адекватного описания процессов деформирования учета граничных условий на всей поверхности границы раздела скелет - флюид. Попытки игнорировать это обстоятельство и заменить контрастную микронеодно-родную среду некоторым эквивалентным сплошным телом приводят к трудностям теоретического объяснения ряда экспериментальных фактов (например, роста коэффициента Пуассона с ростом дифференциального давления в газонасыщенных средах и падения в насыщенных жидкостью [1], падения и роста коэффициента
Пуассона при компактировании металлических гранул в зависимости от внешнего давления [2]).
Возможность вычислять перемещения (внутри либо на границе) достаточно произвольной замкнутой поверхности в зависимости от заданных на границе нагрузок либо нагрузок по заданным на границе перемещениям (либо в случае смешанной задачи) может быть применена как к вычислению средних упругих характеристик микронеоднородной среды путем замены разностных операторов дифференциальными [3], так и к расчету напряженного состояния нефтегазоносного пласта [4], что позволяет, в частности, определить линии тока флюида.
© Сибиряков Е.Б., 2006
Метод граничных интегральных уравнений предназначен для того, чтобы сводить трехмерную задачу к решению двумерных интегральных уравнений. Его суть состоит в том, что перемещения в любой точке объема (в том числе на поверхности) ищутся в виде свертки тензора Грина для пространства с некоторым потенциалом [5]:
1 г
иі (х) =—}г*(х У)рк(у№у
2п
(1)
где по значку к ведется суммирование, точка х заключается в объеме и может стремиться сколь угодно близко к поверхности изнутри (или извне); точка у бежит по поверхности, в ней же берется элемент поверхности. При этом
г*(х у) = - 1
2ц(Х + 2ц) д— д—
(Х + — -г-+ (Х + 3ц)8ік
дхі дхк
1
г (х, у)' (2)
где г — расстояние между фиксированной точкой и точкой, бегущей по поверхности; 8^. — символ Кронекера; X и ц — упругие константы Ламе. Если точка х стремится к поверхности изнутри, то в ней можно вычислить нагрузки:
дг
Рх (г* )(х, у) = х (г* ) +
дпх
+ М(пх Х Г°х (Гік )) =
1
ц8гк + 3(Х + М)
X + 2 ц
д— д— дхі дхк
д
дпх
+ М
д1
^ (ПхХк )------------------COS (ПхХі)------
дх,- — дхр —
(3)
где пх — нормаль в точке х. В результате, если на поверхности заданы нагрузки, то для нахождения потенциала (с помощью которого можно будет найти перемещения в любой точке объема) получится двухмерное сингулярное интегральное уравнение. Интеграл в этом случае нужно понимать в смысле главного значения (из-за последнего члена в (3)). Численное решение подобного рода уравнений достаточно затруднительно и представляет отдельную проблему. По этой причине метод граничных интегральных уравнений использовался только для решения двухмерных задач, хотя преимущества этого метода для решения краевых задач понимают многие исследователи [6].
Использование матрицы Вейля (антенного потенциала), которая с точностью до множителя совпадает с тензором Грина для полупространства и тензор нагрузок от которой не содержит сингулярных членов, представляется более предпочтительным. Однако при по-
пытке использовать этот метод в качестве теста для решения задачи о всестороннем сжатии шара, получится неверный ответ. Это связано с тем, что тензор нагрузок от фундаментального решения третьего рода вычисляется в декартовых лабораторных координатах. Вычислению упомянутого тензора нагрузок и записи полученных уравнений не в декартовых координатах, а в обобщенных ортогональных координатах, связанных с поверхностью, и посвящена данная работа.
2. Использование фундаментального решения третьего рода и тензора нагрузок от него
Тензор Грина для полупространства (фундаментальное решение третьего рода) записывается в виде
[5, 7]:
М
ік(х У) =
1
(Х + 2ц)Гік(X, у) -1 Zik(X, у)
(4)
где точка у помещается в начале координат, а направление оси Хі совпадает с направлением внутренней нормали к поверхности в точке (у, х) имеет координаты (хр Х2, х3),
Zik =
д 2 2 д 1 2 д 1
дх2 дх,дх2 дх,дх3
д 2 < д ^ д 2 V 2 д
дх,дх2 дх,2 дх| дх2 дх3
д 2 < д 2 < д 2 - д 2 <
дх,дх3 дх2дх3 дх,2 дх|
V = х11п (— + х1) - —. При этом Рх (М к) =
: д*. 72
(5)
(6)
= 3
д— дх,
V 1 У д— д—
дх, дх2
д— д— дх, дх3
д— д— дх, дх2 д— 2
дх-,
д— д—
дх2 дх3
д— д— дх, дх3
д— д— дх2 дх3
(д—_ 72
дх
д
дпх
(7)
и не зависит от упругих модулей.
Однако (7) вычислено в предположении, что «без ограничения общности, пх можно считать параллельным пу» [5, 7]. Вообще говоря, «без ограничения общности», можно считать скорее обратное, что пх параллельно пу только в одной точке (когда х совпадает с у, в этой точке и возникает сингулярность). Теперь вычислим Рх (М л) без упомянутого выше предположения, обозначив за ет1 и ет2 единичные орты в направлении х2 и х3, которые являются касательными в точке у и ортогональны нормали в этой же точке и друг другу.
Х
Х
+
Итак,
Рх (Ма ) = 3
Эх1 дг Эг дх1 дх2 Эг Эг Эх1 Эх3
Эг Эг
Эх1 Эх2
: Эг л2
Эх
Эг Эг
Эх2 Эх3
Эг Эг
Эх1 Эх2 Эг Эг
Эх2 Эх3
: Эг 72
Эх
Э 1+ Эпх г
X + ц
(пх , ет1) (пх , ет2)
Эх1Эх3
0
- (пх , ет2)
Эх1Эх2Эх3
(пх , ет1)
Э\
Эх2Эх3
0
--(пх , ет2)
Э \
Эх3 Эх3
( ) Э у ( )
(«х , ет1^-Г - (пх , ет2)
Эх3
Э\
Э \
Эх1Эх2
--(«х , ет1)
Э3-
Эх1Эх2Эх3
(пх , ет2)~Т - (пх, ет1) Эх2
Э\
Как видно, Рх (Мгк) не содержит сингулярных особенностей, и соответствующие интегральные уравнения будут регулярными.
3. Запись Мй и Р*(МЙ) в ортогональной системе координат, связанной с поверхностью
Используя Мл (4) и Рх (Мгг) (8) можно построить систему регулярных интегральных уравнений, находить потенциалы и решать соответствующие краевые задачи. Однако при данном построении тензоров фундаментальных решений и нагрузок от них следует принять во внимание, что координаты хр х2, х3 не есть декартовы координаты неподвижной системы координат, а заданы в подвижном базисе, при этом х;— проекция радиус-вектора г, проведенного из точки у в точку х, а х2 и х3 — проекции этого же вектора на взаимно ортогональные касательные к поверхности в точке у. Можно было бы поместить начало координат в центр масс тела и, используя матрицы поворота, вычислить как декартовы координаты фиксированной и бегущей точек, так и соответствующие проекции тензоров М* и Рх (Мгк) на оси. При этом декартовы проекции тензоров будут иметь достаточно сложный вид. Более разумным представляется записать их проекции на оси координат, связанные с поверхностью.
Для этого вспомним о физическом смысле тензора Грина для полупространства. М^- — это отклик в направлении i на 8-нагрузку, приложенную в направлении к. Например, М11 — это перемещения в направлении внутренней нормали в бегущей точке, которые вызывает 8-нагрузка, приложенная в этом же направлении в начале координат. Но в каждой фиксированной точке существует своя внутренняя нормаль и свои касательные. Поэтому представляется разумным найти проекции перемещений, вызываемых упомянутой выше нагрузкой на нормаль и касательные в фиксированной точке, т.е.
Эх^Эх3
(пх , ет2)
Э3
Эх2Эх3
Эх^Эх3
- (пх , ет1)
Э\
Эх2Эх2
(8)
Мпхпу = М11(пх, пу ) + М21 (пх, ет1 )+ М31 (пх, ет2 X (9)
Мт10пу = М11(ет10, пу ) + М21 (ет10, ет1) +
+ М31(ет10, ет2),
Мт20пу = М11(ет20, пу ) + М 21 (ет20, ет1) +
+ М31(ет20, ет2),
М яхт1 =М12(пх, пу ) + М 22 (пх, ет1) +
+ М 32 (пх, eт2),
М
т10т1
= М12(ет10, пу ) + М 22 (ет10, ет1) +
+ М32 (ет10, ет2),
Мт20т1 = М12(ет20, пу ) + М22 (ет20, ет1) + + М 32 (ет20, ет2),
Мпхт2 = М13 (пх, пу ) + М23 (пх, ет1) +
+ М33 (пх, ет2 ),
Мт10т2 = М13(ет10, пу ) + М 23 (ет10, ет1) + + М 33 (етю, ех2),
Мт20т2 =М31 (ет20, пу) + М32(ет20, ет1) +
+ М 33 (ет20, ет2),
(10)
(11)
(12)
(13)
(14)
(15)
(16) (17)
при этом пх, ет10, ет20 — единичные взаимно ортогональные векторы в направлении внутренней нормали и касательных направлениях в фиксированной точке. При этом потенциал удовлетворяет следующей системе уравнений:
и„х = (М +М „е р +М „е ^^Му,
ит10 = Т-}(Мт10п Рп +Мт10т1Рт1 + Мт10т2Рт2)^у, (18)
2п у у у
ит20 = (Мт20пуРпу +Мт20т1Рт1 + Мт20т2Рт2)^у •
ц
Рис. 1. Использование сферической системы координат с переменным радиусом
Аналогично вычисляются и проекции тензора нагрузок (8). Полученная система будет простой и наглядной, что можно продемонстрировать на следующем примере.
4. Пример использования метода граничных интегральных уравнений для задачи о всестороннем сжатии шара
Среди тех, кто занимался применением метода граничных интегральных уравнений, существует мнение, что потенциал не может быть константой и что простых примеров, иллюстрирующих этот метод, не существует. Воспользуемся простотой записи уравнений в ортогональной системе координат, связанной с поверхностью для аналитического решения задачи о всестороннем сжатии шара. Пусть на поверхности шара даны нормальные перемещения Un = aR, а касательные перемещения равны нулю. Нужно найти нагрузки на поверхности. Ответ известен: Pn =-а (3Л + 2ц), остальные компоненты вектора нагрузок равны нулю. Для удобства поместим начало координат в фиксированную точку, стремящуюся к поверхности изнутри, а оси координат сориентируем так, чтобы внутренняя нормаль в этой точке имела координаты (0, 0, 1) (рис.1).
Запишем координаты точки у в сферической системе координат с переменным радиусом r и центром в точке х. Обозначим за у полярный, а за ~ — азимутальный углы. При этом
er = ex1 sin у cos ~ + ex2 sin у sin ~ + ex3 cos ~,
еу = ex1 cos у cos ~ + ex2 cos у sin ~ - ex3 sin у, (19)
еф =-ex1 sin~P + ex2 cos
Выразим ex1, ex2, ex3 через единичные орты декартовой системы координат i, j, k и сферические координаты точки х:
где р0,0о — классические сферические координаты фиксированной точки. Расстояние между точками х и у:
r = 2R cos у. (21)
где R — радиус сферы. Внутренняя нормаль в бегущей точке у и элемент поверхности будут равны:
ny = -er cos у - еу sin у, (22)
dSy = 2Rr sin ydydp, (23)
При этом интегрирование по у проводится от 0 до п/2, а по ф — от 0 до 2п. Введение такой системы координат представляется оправданным не только для сферы, но и для любой замкнутой поверхности, так как d» будет всегда пропорционально r, что позволит избавиться от интегрируемой особенности вида 1/r при численном решении уравнений. Поскольку в (4) x1, x2 и x3 — не декартовы координаты, а проекции вектора r на нормаль и две взаимно-ортогональные касательные, выберем направления касательных так, чтобы ет1 в точке х совпадал с ex1, а ет2 — с ex2:
ет1 = er sin у cos ф - еу cos у cos ~ - ep sin ф,
~ ~ ~ ~ (24)
eT2 = er sin у cos ~ - еу cos у cos ~ + ep sin ~ cos ф.
При этом
nx = ex3, ет10 = ex1, ет20 = ex2, (-er, ny) = cos^
(-eT, ет1) = - sin у cos ~,
(-ет, eT2) = - sin у sin
В этом случае тензор М ik принимает очень простой вид:
(25)
M„ =
M22 =
і
2цг
і
2цг
ц 2
і +------+ cos у
А + ц
і+
А + ц
cos у і + cos у
7
+ sin у
2ф
2 ф ц sin ф
cos^ + :
77
А + ц (і + cos у)
/у
M33 =
і , ц cos у і+—---------------------—+
2цг + sin2 у
А + ц і + cos у
: 2 ф . 2- ц cos ф sin 2 ф + -
А + ц (і + cos у)2
(26)
/у
Mu =
і
2цг
sin у cos ф
cos у+
ц
і
77
А + ц і + cos у
ex! і - cos 2 ф0(і + cos 60) - sin ф0 cos ф0 (і + cos 60) sin 60 cos ф0 i
ex2 = - sinф0 cosф0(і + cos 60) і - sin2ф0(і + cos 60) sin 60 sin ф0 j
_ex3 _ - sin 60 cos ф0 - sin 60 sin ф0 - cos 60 k
ц
+
М 21 =
М13 =
М31 =
1
2цг
1
2цг
1
2цг
sin у cos ф
cos у
ц
1
sin у sin ф
cos у +
А + ц 1 + cos у ц 1
sin у sin ф
М32 = М 23 =
1
2цг
sin у sin ф cos ф
cos у
1-
А + ц 1 + cos у
ц1 А + ц 1 + cos у
А + ц (1 + cos у)2
Тензор нагрузок от Мл вычисляется в соответствии с (8) с условием, что: x1 = r cos у, x2 = -r sin у cosф, x3 = -r sin у sin ~,
(nx, ny ) = - cos2^
(ет1, nx) = sin 2у cos ~,
(ny, ет1о) = - sin 2у cos ~,
(ет1, ет1о) = 1 - 2cos2 у cos2 ~, (27)
(ny, ет2о) = - sin 2у sinф,
(ет1, ет20) = -2cos2 у sin ~cos ~,
(nx, ет2) = sin 2у sinф,
(ет2, ет10) = -2cos2 у sinф cos ~,
(ет2, ет20) = 1- 2cos2 у sin2
Теперь воспользуемся (18) и найдем потенциал. Поскольку J Мт10п dSy = 0 и J Мт20п dSy = 0, то решение (18) ищется в виде: Fn = const, Ft1 = Ft2 = °.
М n
1
2цг
ц2
1+ . + cos2 у
1
2цг
sin у cos у
А+ц ц
sin у
А + ц 1 + cos у
cos 2у -sin 2у.
Соответственно
3а(А + ц)
Fn =------------.
(28)
(29)
Получаем, что Рт10 = Рт20 = 0. Для вычисления Рп уч-
д 1 х тем, что Р содержит член-------, а точка х стремится к
дпх г
поверхности изнутри:
Pn =-3а(А + ц) - 3а(А + ц) +ац = -(3А + 2ц)а, (30)
что совпадает с правильным ответом.
5. Комментарий
Если на поверхности заданы нагрузки, то для нахождения потенциала благодаря свойствам оператора д 1
------получим уравнение второго рода, решение ко-
dnx r
торого не представляет трудностей. Если же задан вектор перемещений, то для нахождения потенциала потребуется решать интегральное уравнение первого рода. Для его решения необходимо правильно выбрать нулевое приближение. Представляется физически оправданным выбрать нулевое приближение в виде Fi = = const • Ut и определять константу как среднее значение по поверхности при подстановке в систему (18).
6. Выводы
Таким образом, трудности, препятствующие широкому применению метода граничных интегральных уравнений для решения трехмерных задач, следует считать устраненными. Использование иначе построенного тензора фундаментальных решений и тензора нагрузок от него в системе координат, связанной с поверхностью, дает не только простоту и физическую наглядность, но и существенно упрощает вычислительные процедуры за счет того, что уравнения становятся регулярными.
Литература
1. Carcione J.M., Cavallini F. Poisson ratio at high pore pressure // Geophysical Prospecting. - 2002. - V. 50. - P. 97-106.
2. Carnavas P. C., Page N.W. Elastic properties of compacted metal powders. // Journal of Materials Science. - 1998. - V. 33. - P. 4647-4655.
3. Сибиряков Е.Б. Зависимость между коэффициентом Пуассона и микроструктурой в микронеоднородной среде // Физ. мезомех. -2004. - T. 7. - № 1. - С. 63-68.
4. Сибиряков Б.П., Сибиряков Е.Б., Глебов А.Ф., Нестеров В.Н., Соколов Е.П. Прогноз напряженного состояния месторождений методами многоволновой сейсморазведки как средство интенсификации добычи нефти (на примере верхнеюрского пласта Ю. Аригольского месторождения) // Геология и геофизика. - 2004. -T. 45. - № 4. - С. 709-715.
5. Купрадзе В.Д. Методы потенциала в теории упругости. - M.: Физматгиз, 1963. - 472 с.
6. Джонсон К. Механика контактного взаимодействия. - М.: Мир,
1989. - 510 c.
7. Партон В.З., Перлин П.И. Методы математической теории упругости. - М.: Наука, 1981. - 688 c.
ц
