CC BY
17
УДК 517.54
Численное решение линейной двумерной динамической задачи для пористых сред
Константин Э. Сорокин* Холматжон Х. Имомназаров^
Институт вычислительной математики и математической геофизики СО РАН, Лаврентьева 6, Новосибирск, 630090,
Россия
Получена 18.12.2009, окончательный вариант 25.02.2010, принята к печати 10.04.2010 В работе численно решена линейная двумерная динамическая задача для пористых сред, насыщенных жидкостью. За основу взята линеаризованная модель В.Н.Доровского, в которой среда описывается тремя упругими параметрами. Для решения задачи использовался явный разностный метод второго порядка точности. Также использован метод ГЫЬ (полностью согласованного поглощающего слоя). Проведена серия тестовых численных расчетов для пробной модели среды.
Ключевые слова: пористая среда, метод конечных разностей, шахматная сетка, метод ГЫЬ.
Постановка задачи
В этом разделе представлено краткое описание и постановка задачи распространения сейсмических волн в пористой среде [1]. Пористая среда, состоящая из заполненной вязкой жидкостью упругодеформируемой матрицы, является реалистической моделью. Данная модель позволяет объяснить эффекты сейсмических исследований свойств горных пород при наличии поровой жидкости. Для численного моделирования, как правило, используется модель Френкеля-Био [2,3]. Особенность распространения волнового поля в пористых средах — наличие дополнительной второй продольной волны. Скорости распространения сейсмических волн в теории Френкеля-Био есть функции четырех упругих параметров. В 1989 году В.Н. Доровским [4] была построена нелинейная математическая модель для пористых сред. В этой модели, так же как в теории Френкеля-Био, есть три типа звуковых колебаний: поперечное и два продольных. Но при этом среда описывается тремя упругими параметрами [5,6]. Эти упругие параметры взаимно однозначно выражаются тремя скоростями упругих колебаний. Указанное обстоятельство важно для численного моделирования распространения упругих волн в пористых средах, когда заданными параметрами являются значения скоростей акустических волн, пористости, физических плотностей матрицы и насыщающей жидкости.
Перейдем к математической формулировке модели. Пусть полуплоскость у > 0 заполнена пористой средой и насыщенной жидкостью. Тогда распространение сейсмических волн в данной среде при отсутствии потери энергии описывается следующей начально-краевой задачей [5-7]:
t [email protected] © Siberian Federal University. All rights reserved
+ (дж 0"жж + ду ^жу ) + дЖр -- ^Ж7
Р0,я Р0
д(Му + (дж^жу + ду ^уу ) + дур — 7
Р0,я Ро
д4«ж + джр ^ж,
Р0
д^у + — ду р — ,
Р0
д(0жж + 2РджМж + (_К — 7^Р)(джМж + дуМу)-------~К(дж«ж + ду«у) — О,
,р0,^ 2 | л \ р0,в
----К — 7ГР)(джМж + дуМу)-----------
Р0 3 Р0
д^жу + Р(ду Мж + джМу )
д*°уу + 2Рдуиу + ( Р-^ К — — р)(джМж + дуМу ) — Р-^К(дж«ж + ду«у ) — 0?
Р0 3 Р0
д4р - (К - «Р0Р0,я)(дж Мж + ду Му) + «Р0Р0,г(дж«ж + ду «у) — О,
Мж|4=0 — Му |4=0 — «ж|4=0 — «у |4=0 — ажж|4=0 — ^жу |4=0 — ауу |4=0 — р|4=0 — 0? ауу + р|у=0 — ажу |у=0 — “Р00- р|у=0 — 07 где м — (мж,му) иг/ — («ж,«у) — векторы скорости упругого пористого тела с парциальной плотностью Р0,л и жидкости с парциальной плотностью р0,; соответственно; р — поровое давление; <гжж, <гжу, <гуу — компоненты тензора напряжений; — (^ж,^у) — вектор массовых сил; Р0 — Р0,г + Р0,«, Р0,« — р0,3 (1 - 4>), Р0,г — Ро,;^>, р0,8 и р°,; — физические плотности упругого пористого тела и жидкости соответственно; — пористость; К — А + 2р/3, А > О, р > 0 — коэффициенты Ламе; а — Р0«з + К/р^; р^аз > 0 — модуль объемного сжатия жидкой компоненты гетерофазной среды. Упругие модули К, р, аз выражаются через скорость распространения поперечной волны ея и две скорости продольных волн ср1, ср2 следующими формулами [8,9]:
Р — Р0,вс:;,
К — р0 р0,« /2 , с2 8 Р0,; с2 //2 с2 ) 64 Р0,;Р0,« с4
К — У Р7 Гр1 + СР2 - 3 с* - V (СР1 - СР2) - У “7Т\
1 / с2 , с2 8 Р0,« с2 /(с2 с2 ) 64 Р0,;Р0,« с4
—Р0 Гр1 + Ср2 - 3 77^ - V(Ср1 - Ср2) - ^
PML-модель
Рассматривается РМЬ-модель для линейной двумерной динамической задачи. Одним из основных принципов РМЬ-модели является соединение уравнения в левом полупространстве с уравнением в правом таким образом, чтобы граница этих полупространств не порождала искусственных отражений, а волна в правом полупространстве уменьшалась экспоненциально. Соответственно, вводим функцию й(ж), равную нулю в левом полупространстве и положительную в правом. Эта функция будет играть роль поглощающего множителя. Уравнения РМЬ-модели можно записать в следующем виде:
(д + ^(х))мж дж^жж джр + ^ж, ^ж дустжу + ^ж?
Р0,я Р0 Р0,я
(д4 + ^(х))иу дж^жу + ; д£иУ ду^уу дур + 7
Р0,5 Р0,5 Р0
(д4 + .(ж))^ = -2рджМж - ( К - -р) джмж + Кд^; 54а|ж =
Ро 3 ) ро
_ А Ро,^ 2 \ Р0,8
= - ( ~К - 3 Р ) дУ“у + “Р0 КдУ«у 7
(д4 + .(ж))ст^ = -рджМу; д^у = -рдумж,
(д4 + й(ж))а^„ = - ( ^К - 2р ) дхмж + Р^Кдж«ж; д4а|у =
уу V Р0 “ 3-; — ■ Р0
= -2рду иу - ^ Р-^ К - - Р^ ду иу + Р°- Кду «у ,
V Р0 3 / Р0
(д4 + .(ж))«^ = —1 дхр + ; д4«| = ^ж,
Р0
(д4 + <ж))^ = ; д4«| = -Р0дур + ^у,
(д4 + .(ж))р^ = (К - аР0Р0,8)джМх - аР0Р0,гдх«ж; д^11 = (К - аР0Р0,8)дуМу - аР0Р0,гду«у.
Дискретная РМЬ-модель
Рассмотрим дискретную РМЬ-модель для приведенной выше задачи распространения. Используем конечно-разностный метод. Схема строится на шахматной сетке [10]. Скорости мх и «х считают в точках (г,]), (ж* = г/, у^- = ^’/г.) координатной сетки, му и «у — в точках
+ -,^ + -) , ^хх, ауу, Р — в точках и аЖу — в ^г, + -
записана следующим образом (ввиду громоздкости ниже приведена разностная аппроксимация только уравнения для компоненты скорости мх, остальные разностные уравнения строятся по аналогии [10]):
(“х)^ = (“Х)^ + (иХ)ПП^,
^)£1 - (их)ПП^ , .х (их)ПП+1 + «^
+ гч
Д*
(^хх)”++ 2 - (^х)”^ (р)”++ 2 - (р)”^ ,
* 2 * 2 * 2 * 2_Л + (^ )"+ 2
Р0,я/ Р0/
(■»;?1 - <°хй + .у мс1+«й = ^>А- <^»>;& ,„+2
Д* + * 2 р0,8/ +( х)^' ,
где используются обозначения гк = г + к и ^ = ] + к.
Система уравнений согласована с моделью РМЬ в углах, где поглощение идет в обоих направлениях. Для вычисления решения внутри РМЬ-слоев в направлении ж мы используем ту же самую систему уравнений, где берём .у = 0, и для РМЬ-слоев в направлении у мы берём = 0 (рис. 1).
Рис. 1. PML-слои: в углах dx > 0 и dy > 0, в слоях в направлении x(y) считаем dy = 0(dx = 0)
Численные результаты
В следующем примере мы моделируем распространение сейсмических волновых полей в бесконечной двумерной среде. Представлена однородная, изотропная упруго-пористая среда, заполненная вязкой жидкостью. Материал, заполняющий нашу область, будем характеризовать скоростями продольных и поперечных волн, соответственно cpi = 2000 m/s, cp2 = 450 m/s, cs = 1400 m/s. Определим ограниченную область D квадратной формы со стороной с = 30 m, в которой заданы нулевые начальные условия, и на каждой из четырех границ есть PML-слой. Зададим сетку с количеством узлов по каждой координате N = 400. Временной шаг возьмем равным At = 0, 9h/ (y/2Cp7), где h = c/N. В наших численных экспериментах взрыв будет происходить в точке (xo, уо) = (7, 5m, 7, 5m) и задаваться функцией:
Fx(t,x,y) = f (t) д6 ( Xo) 6 (у - yo) ,Fy(t,x,y) = f (t)6 (x - xo) d6 ^ Уо),
f (t) =1 -2n2f2(t - to)e-n fo (t-to) , t < 2to,
t > 2t
o-
Здесь *0 = 1//0, /0 = е8/(/Ж;) — средняя частота, Ж; = 20 — коэффициент, отражающий число точек на длину волны. Поглощающий коэффициент .(ж) моделируется следующим образом:
~ (ж ^ ~ ( 1 И 3ср1
.(ж) = МV •<!° = 1оня) 1Т'
где £ = 10/ — толщина нашего слоя, а ^0 — функция коэффициента отражения Я = 0,001.
Результаты численных расчетов волновых полей на заданной модели среды представлены на рис. 2-4. На этих изображенях видно образование дополнительной продольной волны в пористой среде и отсутствие отражений на границе исследуемой области вычислений.
(а) (Ь)
Рис. 2. Поля скорости (а): ||и|| = ^+ иуу; (Ь): ||«У = + «2 при * = 2, 76 шэ
(а) (Ь)
Рис. 3. Поле скорости (а): ||и|| = + иуу; (Ь): || = + «2 при * = 8, 35 шэ
(а) (Ь)
Рис. 4. Поле скорости (а): ||и|| = ^+ иуу; (Ь): ||«|| = у^«2 + «2 при * = 16, 7 шэ
Список литературы
[1] Х.Х. Имомназаров, А.А. Михайлов, Использование спектрального метода Лагер-ра для решения линейной двумерной динамической задачи для пористых сред, Сибирский журнал индустриальной математики, 11(2008), №3(35), 86-95.
[2] Я.И. Френкель, К теории сейсмических и сейсмоэлектрических явлений во влажной почве, Изв. АН СССР, сер. География и геофизика, 8(1944), №4, 133-146.
[3] M.A. Biot, Theory of propagation of elastic waves in fluid-saturated porous solid. I. Low-frequency range, J. Acoustic. Soc. America, 28(1956), 168-178.
[4] В.Н. Доровский, Континуальная теория фильтрации, Геология и геофизика, (1989), №7, 39-45.
[5] В.Н. Доровский, Ю.В. Перепечко, Е.И. Роменский, Волновые процессы в насыщенных пористых упругодеформируемых средах, Физика горения и взрыва., (1993), №1, 100-111.
[6] A.M. Blokhin, V.N. Dorovsky, Mathematical Modelling in the Theory of Multivelocity Continuum, N. Y., Nova Science, 1995.
[7] Kh.Kh. Imomnazarov, A mathematical model for the movement of a conducting liquid through a conducting porous medium. I. Excitation of oscillation of the magnetic field by the surface rayleigh wave, Math. Comput. Modelling., 24(1996), №1, 79-84.
[8] Х.Х. Имомназаров, Несколько замечаний о системе уравнений Био, Докл. РАН, 373(2000), №4, 536-537.
[9] Kh.Kh. Imomnazarov, Some remarks on the Biot system of equations describing wave propagation in a porous medium, Appl. Math. Lett., 13(2000), №3, 33-35.
[10] Fr. Collino, Ch. Tsogka, Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Rapport de recherche, №3471, Aout 1998.
Numerical Solving of the Liner Two-dimensional Dynamic Problem in Liquid-filled Porous Media
Konstantin E. Sorokin Kholmatjon Kh. Imomnazarov
In this paper the linear two-dimensional dynamic problem in a liquid-filled porous media is numerically solved. As the foundation we use the V.N. Dorovsky linearized model in which the media is described by three elastic parameters. To solve this problem we use a finite-differences scheme of the second order of accuracy. The PML model is also used in this paper. We present several numerical results for a test media which show the efficiency of this model.
Keywords: porous media, finite-difference schem, staggered grid, PML-model.
