Эффективные коэффициенты в задаче распространения акустических волн в многомасштабной изотропной среде

Вычислительные технологии

Том 19, № 6, 2014

Эффективные коэффициенты в задаче распространения акустических волн в многомасштабной изотропной среде

Е.П. Курочкинл1, О.Н. Соболева2

В рамках метода подсеточного моделирования получены эффективные коэффициенты для задачи распространения акустических волн в многомасштабной неоднородной среде. Коррелированные поля плотности и модуля упругости моделируются непрерывными мультипликативными каскадами с логарифмически нормальными распределениями вероятностей. Предполагается, что длина волны много больше максимального масштаба неоднородностей среды. Полученные теоретические результаты сравниваются с результатами прямого численного моделирования.

Ключевые слова: распространение волн, эффективные коэффициенты, подсе-точное моделирование, многомасштабные случайные среды.

Введение

Моделирование распространения волн в неоднородной среде — задача, представляющая большой научный и практический интерес. Такие задачи возникают при исследовании подземных ядерных взрывов, землятресений, в океанографии и в медицине [1]. В процессе численного моделирования крупные неоднородные включения в среде учитываются в математической модели непосредственно с помощью граничных условий. При этом учёт всех масштабов вариаций физических параметров требует громадных вычислительных затрат. Чтобы решить данную проблему, при численном моделировании неупорядоченных сред используют эффективные коэффициенты, задачи построения которых для слоистой среды рассматривались в работах [2, 3]. Для построения коэффициентов в периодической детерминированной среде применялись методы гомогенизации, хорошо работающие в искусственных средах. Методы гомогенизации в детерминированных непериодических двухмасштабных средах для двумерной задачи SH-волн, а затем для P- и SV-волн рассматривались в работах [4, 5]. Геофизическая среда, как правило, является непериодической и многомасштабной, к тому же координаты мелкомасштабных неоднородностей точно не известны. Поэтому в данном случае среда описывается статистическими полями. Геофизические неупорядоченные среды хорошо представляются мультипликативными иерархическими каскадами с негауссовым распределением вероятностей и фрактальными моделями [6, 7]. Данный факт позволяет для получения эффективных коэффициентов применять метод подсеточного моделирования. В этом случае строится более простая модель, требующая меньшего количества вычислительных затрат, правильно описывающих поведение решения в крупномас-

1 Институт теплофизики СО РАН, Новосибирск, Россия

2 Институт вычислительной математики и математической геофизики, Новосибирск, Россия

Контактный e-mail: olga@nmsf.sscc.ru

штабном пределе [8]. В настоящей работе с помощью метода подсеточного моделирования получены уравнения для эффективных коэффициентов в волновом уравнении, если модуль упругости аппроксимируется непрерывным мультипликативным каскадом с логарифмически нормальным распределением. Предполагается, что длина волны много больше максимального масштаба неоднородностей среды. Волны возбуждаются источником с доминантной частотой и шириной спектра Полученные теоретические результаты сравниваются с результатами прямого численного моделирования для двумерной среды.

1. Постановка задачи

Распространение акустической волны в неоднородной многомасштабной среде описывается уравнением

. . д2и (х,£) д Л / ч д . Л .

р

где р (х) — плотность, х — вектор пространственных координат, £ — время, Л (х) — модуль объемной упругости, и (х, £) — смещение, ^(х, £) — источник с доминантной частотой ш0 и шириной спектра ш1. По повторяющимся индексам подразумевается суммирование. Для моделирования случайных полей р (х), Л (х) используется подход, подробно описанный в работе [9].

Пусть поле модуля упругости Л (х) известно. Это означает, что выполнено его измерение на некотором масштабе /0 в каждой точке х. Чтобы перейти к более грубой сетке масштабов, недостаточно сгладить Л (х)г по масштабу /, I > /0, так как сглаженное поле не будет адекватно представлять физический процесс, описываемый уравнением (1) на интервале масштабов (/, Ь), где Ь — максимальный масштаб неоднородности среды. Это объясняется тем, что флуктуации упругости Л (х) на интервале масштабов (/0,/) корд

релируют с флуктуациями поля ——и (х, £) и эти корреляции могут быть достаточно

дхз

большими. Для построения модели среды, как и в работе [10], рассматривается безразмерное поле 0, равное отношению полей, полученных сглаживанием Л (х)г по двум различным, близким к /0, масштабам /, /'. Обозначим через Л(х) сглаженное по масштабу I поле Л1о(х). Тогда 0(х,/,/') = Л(х)^/Л(х), I' < /. Случайное поле 0 по сравнению с полями Л(х)г/, Л(х) меняется плавно. Раскладывая поле 0 в ряд относительно I' — I и оставляя только члены первого порядка малости, при I' ^ I получим уравнение

^ =

где <^(х,/') = (д—(х, /', 1'у)/ду) |у=1. Фактически мелкомасштабные флуктуации поля ^ могут наблюдаться только в некотором конечном диапазоне масштабов 10 < I < Ь. Решение уравнения (2) имеет вид

Л^о

\ 1о

где Л0 - константа. Согласно теореме о суммах независимых случайных полей [11], если дисперсия <^

стремится к полю с нормальным распределением вероятностей. Если дисперсия поля <^(х, /) бесконечна и существует невырожденное (не сосредоточенное в одной точке) предельное распределение суммы случайных величин, то это распределение является устойчивым. В данной работе предполагается, что поле <^(х, /) имеет нормальное распределение вероятности и изотропную однородную корреляционную функцию

<р

(угловые скобки означают статистическое усреднение). Из формулы (4) следует, что флуктуации поля ^ в разных масштабах не коррелируют. Это обычное предположение для скейлинговых моделей соответствует тому факту, что статистическая зависимость становится незначительной, если масштабы флуктуаций параметров различны по величине [10]. В случае масштабноинвариантной среды для любого положительного значения К выполняется условие

Ф^(|х - y | , l, l') = Ф^(К |x - y| , Kl, Kl').

Коэффициент плотности р(х) так же, как поле проводимости, моделируется мультипликативным каскадом:

Pío(х) = Ро exP | ^У X(x,li)^

Предполагается, что функция х(х,1) имеет нормальное распределение вероятностей и дельта-коррелированна по логарифму от масштаба /. Исходя из физического смысла крэффициента плотности для любого / должно удовлетворяеться условие <рДх)> = р0. В случае масштабноинвариантной среды это равенство выполняется при

ФХХ = 2 < X > .

Корреляционная функция между полями плотности и модулем упругости определяется корреляцией между полями х(х,/) и <^(х,/):

Фдх

Если среда не масштабноинвариантна, то величины Фдд, Ф^х, Фдх зависят от масштаба /.

2. Подсеточная модель

Функции плотности и модуль упругости р(х) = р(х)1о, Л (х) = Л (х)г разделим на две компоненты относительно масштаба /. Крупномасштабные компоненты Л (х,/), р (х, /) получены статистическим усреднением по всем <^(х, /1) и х(х, /1) для 10 < 1\ < /, / - 10 = где мало. Мелкомасштабные (подсеточные) компоненты равны Л'(х) = Л(х) -Л(х, /), р'(х) = р(х) - р(х, /):

р(х, /) = р0 ехр

— / Х(х.'1) £

ехр

— / х(хЛ) ^

р'

ехр

— I £

ехр

—I ^ £

. ¿о

¿0

-1

<р'(х)> = 0,

Л(х, /) = Л0 ехр

— I ) £

ехр

—/ £

Л'

ехр

— / ^(хЛ) 17"

ехр

— / ^(хЛ) 111

1

<Л'(х)> = 0.

Из формул (8) следует, что с точностью до членов второго порядка малости по 1///

Г 1/ 1 1/1

р(х,/) * 1 — <Х> у + ^ФХХ (0,/) у рг(х),

1/ 1 1/" Л(х, /) * 1 — <^> 1 + (0,/) 1 Лг(х),

1/

<Л' (х) Л' (х')> * Ф^ (|х — х'| , /) Л (х, /) у,

1/

<р' (х) р' (х')> * Фхх (|х — х'| ,/) р (х,/) у,

1/

<р' (х) Л' (х')> * (|х — х'| , /) р (х, /) Л (х, /) у.

Сделаем в уравнении (1) преобразование Фурье по времени:

д ( д \

ш2р (х) и (ш, х) + —— ( Л (х) —— и (ш, хм = —^ (ш, х).

дж,

дж,

10)

Крупномасштабная (надсеточная) компонента смещения и (ш, х, /) получается как усредненное решение уравнения (10), в котором крупномасштабные компоненты Л(х, /), р (х, /) фиксированы, а мелкомасштабные Л'(х), р'(х) — случайные поля. Подсеточная компонента смещения и' (х) = и (ш, х) — и (ш, х,/). Подставим выражения для и (ш, х), Л(х), р (х) в уравнение (10) и усредним по мелкомасштабным компонентам:

дд

ш2р (х, /) и (ш, х, /) + ш2 <р' (х) и' (х)> + —— ( Л (х, /) —— и (ш, х, /) ) +

джА дж, /

д

дж,

(Л' (х) (х^ = —^ ^ х).

11)

, д \

Подсеточные члены <р' (х) и' (х) >, ( Л' (х) ——-и' (х) у в уравнении (11) не известны и без

предварительной оценки не могут быть учтены. Оценки данных выражений определяют подсеточную модель. Эти члены оцениваются с помощью теории возмущений. Вычитая уравнение (11) из уравнения (1) и оставляя только члены первого порядка малости, получим подсеточное уравнение

и2р

дх2 oxj ¿Х

Крупномасштабная переменная смещения и (и, х,/) в правой части (12) считается известной. Используя метод "замораживания коэффициентов", запишем решение уравнения (12) в двухмерной среде:

и'

V j j

р

где г = |х - х'|. Обозначим к = и\ ——г-. Учитывая (13), оценим подсеточные члены

Л

в уравнении (11):

* = и2 <р' (х) и' (х) > =

I Я1(кг) I

4Л

Я° (кг^Л' (х') £¿7и (и, х', /) - и2р' (х') и (и, х', /^ ¿х'^} ,

6'2 = (л' <х> дХт и' (х)> = ( ^ дХт кх/у / Я <кг> (4Л' <х'> ^«<и.х'-/)+

\ V '' ''

+и2р'

Для полей, в которых небольшое изменение масштаба вызывает значительные изменения самого поля

д д

можно вынести из под знака интеграла в (2). Учитывая, что —— Н (кг) = -^-гН (кг),

дxj дxj

получим

* = 41^7) / Я (кг) ¿1 <р' (х) Л' (х,)> ¿х'д|и (и. хО-

V

,4

ги ' Яо1

4Л (х,/) у

V

=— щя! аХт Я0 4 <Л' <х>Л' 1х' 4и <ш-х-/)+

V 3 3

?ш2 Р д

' д -Я0

4Л (х, /) у дХ

т

V

Подставляя три последние формулы из (2) в (15), имеем

й = 4 / Я0 (кг) ^ФХ^ (г, /) 1х' уш2р (х, /) ^и (ш, х', /) —

V 3 3 7^2 Г 1/

— Т Я° (кг) ФХХ (г,/) 1х' уш2р (х,/) и (ш, х,/),

V

7 С д д 1/ д

й = — 4 У ^Я (кг) ^(г, /) 1х'-/"Л (x, /) — и (ш, x, /) +

V т 3 3

7 Г д 1/

+7 -дгЯо1

4 ,/ дхт /

V

Поскольку масштабы неоднородностей много меньше длины волны и, следовательно, корреляционные функции очень быстро затухают вне области Ь3, в формулах (16) можно перейти к бесконечным пределам интегрирования. Переходя в (16) к полярным ко-

2п 2п

ординатам, интегрируя по частям и используя формулы У = 0, У птп3- = ,

о о

где п, = ж,/г, получим

те

7к2п Г 1/

& =--— Яо1 (кг) гФхх (г, /) 1г—ш2р (х, /) и (ш, х, /) ,

0

те

й = — ук2 /" Я1 (кг) гФ^ (г, /) 1гуЛ (х, /) -дги (ш, х, /) —

4 и / дхт

о

1 1/ д

— ^

где Д = 2 — размерность пространства. Интегралы в (17) порядка Ь. Если Ь много меньше длины волны, то выполняется неравенство Ьш2р (х, /) /Л (х, /) ^ 1, где ш1 — ширина спектра источника ^. Поэтому интегральными членами в (17) можно пренебречь. Это означает, что подсеточный член <р' (х) и' (х)> в первом приближении по масштабу сглаживания 1/// не влияет на параметры эффективного уравнения. Потребуем, чтобы эффективный упругий модуль определялся по формуле того же вида, что и Л^0 (х):

Лг

Подставляя полученную оценку подсеточных членов в надсеточное уравнение и используя (2), имеем

ш2р (х, /) и (ш, X, /) + Л - (0,/) У) Л - ^ | + (0,1) у) X

5 / / [ , 5

хЛот^ I ехр I - ^

,

Пренебрегая членами второго порядка малости по У///, получим

Ло 1 =(1 -<*>) у + Ф- (0,/) у ) Ло.

Переходя в этом уравнении к пределу по У/, получим уравнение, которому удовлетворяет эффективный коэффициент Ло1:

^ = 2

В масштабноинвариантной среде коэффициент Ф^ (0, /) не зависит от масштаба / и решение уравнения (20) имеет вид

/ « \ *Г-<*>)

Ло 1 = Ло Л ^ .

Эффективный коэффициент не зависит от ш. В следующем разделе формула (21) проверяется численным моделированием.

3. Численное моделирование

Численный эксперимент проводится в прямоугольной области на сетке 2048 х 1024 точек с шагом Н = 5 м. Источник находится в каждой точке с координатами х1 = 0, х2 = ^Н, временная зависимость (импульс Риккера) с доминирующей частотой 1 Гц

/(¿) = (1 - 2п2 (¿о - ¿)2) ехр (-п2 (¿о - ¿)2) .

Область интегрирования разбивается на три подобласти. В областях 0 < х1 < 900Н, 1924Н < х1 < 2048Н, 0 < х2 < 1024Н среда однородная, ро = 2000кг/м2, скорость распространения волны = 3000 м/с. В области 900Н < х1 < 1924Н, 0 < х2 < 1024Н средняя скорость волны равна V, а коэффициенты Л, р моделируются мультипликативными консервативными каскадами. Интегралы в (3), (5) аппроксимируются конечно-разностными формулами, в которых удобно перейти к логарифму по основанию 2:

0 0 - 12 <р(х,Тк)Аг - £ х(х,Тк)Ат

Л

В формулах (22) / = 2Т, Дт — шаг по т. В расчётах Дт принимается равным 1. В показателе степени формул (22) стоит сумма статистически независимых слоёв, поскольку предполагается статистическая независимость для полей с различными масштабами.

Для расчётов используются три слагаемых е к = -5, -4, -3. Это соответствует тому, что неоднородности имеют три масштаба 1/128,1/64,1/32 от длины волны. Остальные слагаемые полагаются равными нулю. Количество слагаемых выбиралось так, чтобы масштаб Ь самых крупных вариаций позволил заменить приближённо вероятностные средние величины усреднёнными по х2 и выполнялось неравенство Ь^2А/р ^ 1. Масштаб самых мелких вариаций выбирался так, чтобы разностная задача хорошо аппроксимировала уравнение (1). При каждом тк поля х строятся по формулам

/ФГ ФГ

) = УьГ С1 (Х,Тк ) + ,

х(х, Тк) = (гС1 (Х, Тк) + УГ-72С2(Х, тк)) + ФХХ,

где , С2 _ гауссовы поля с единичной дисперсией, нулевым средним и корреляционной функцией

«1 (х,Тк)С1 (У, т)) - (С1 (х,Тк)) (С1 (У, т )) = (С2 (х,Тк)С2(у, т )) - «2 (х,Тк)) «2 (У, Т; )> =

= ехр [-

Случайные поля ^, С2 моделируются с помощью метода, изложенного в [12]. Коэффициенты Ф0 должны выбираться из экспериментальных данных. В нашем расчёте Фо° = ФдХ = 0.2, < х >= 0.1, < ^ >= 0.1, г = 0.9. Уравнение (1) решалось разностной схемой второго порядка точности по времени и пространству, с постоянным шагом и граничными условиями ди/дх2 = 0 на границах 0 < х1 < 2048^, х2 = 0, и х2 = 1024^. Результаты расчётов приведены на рисунках 1, 2. В соответствии с процедурой вывода подсеточных формул для проверки необходимо много раз решить полную задачу

Рис. 1. Мгновенный снимок реализации поля смещения. Плотность и модуль упругости моделируются по формуле (22) для трёх масштабов т = -5, -4, -3, поля плотности и модуль упругости коррелированны, коэффициент г = 0.9

Рис. 2. Мгновенный снимок усреднённого поля смещения. Минимальный масштаб неоднород-ностей 1/128, максимальный — 1/32 от длины волны. 1 — модуль упругости и плотность постоянны и равны Ао, р0; 2 — эффективное смещение; 3 — смещение, усреднённое по 200 реализациям, расчёты проводились с постоянной плотностью р0; 4 — смещение, усреднённое по 200 реализациям, плотность и модуль упругости — коррелированные случайные поля, коэффициент г = 0.9

и выполнить вероятностное усреднение по мелкомасштабным вариациям. В результате получится решение, которое можно сравнивать с решением эффективного уравнения (см. рис. 2, кривая 2). Для сравнения приводится решение уравнения со средними значениями модуля упругости и плотности (кривая 1). Усреднение по ансамблю Гибб-са требует многократного решения полной задачи. В работе используется более экономный вариант проверки. Решение при каждом значении х\ усредняется по отрезку 100Н < х2 < 924Н, а затем проводится доусреднение по ансамблю. Расчёты проведены для 200 реализаций. Кривая 3 на рис. 2 — среднее смещение при постоянной плотности р0, кривая 4 — результат численного моделирования, если плотность и модуль упругости моделируются мультипликативными каскадми. Близость кривых 2, 3 и 4 подтверждает теоретический вывод о том, что плотность в первом приближении по масштабу измерения не влияет на эффективные параметры.

Заключение

В работе получены эффективные коэффициенты для волнового уравнения в случае, если они описываются крайне нерегулярными непрерывными полями, близкими к мульти-фракталам. Если минимальный масштаб 10 устремить к нулю, то поле в (3) стремится к мультифракталу. Поскольку минимальный масштаб остаётся конечным, то какие-

либо сингулярности отсутствуют. При этом канторовы множества не возникают и весь анализ не выходит за пределы аппарата дифференциальных уравнений и теории случайных процессов. Данный подход позволяет избежать проблем с предельным переходом от дискретного поля к непрерывному, часто возникающих при исследовании поведения волн во фрактальных средах. В масштабноинвариантной среде эффективные коэффициенты степенным образом зависят от масштаба сглаживания. Прямое численное моделирование показывает, что представленное эффективное уравнение дает хорошее согласие с усреднённым решением, полученным прямым численным моделированием.

Список литературы

[1] Bleistein N., CohenJ.K., Stockwell J.W. Mathematics of Multidimensional Seismic Imaging, Migration, and Inversion. New York: Springer, 2001.

[2] Backus G. Long-wave elastic anisotropy produced by horizontal layering //J. Geophys. Res. 1962. Vol. 67, No. 11. P. 4427-4440.

[3] Молотков Л.А. Об эквивалентности слоисто-периодических и трансверсально-изотропных сред // Математические вопросы теории распространения волн: Записки научного семинара ЛОМИ. 1979. № 89. С. 219-233.

[4] Guillot L., Capdevile Y., Marigo J.J. 2-D non periodic homogenization for the SH wave equation // Geophys. J. Intern. 2010. Vol. 182. P. 1438-1454.

[5] Capdeville Y., Marigo J.J. A non-periodic two scale asymptotic method to take account of rough topographies for 2-D elastic wave propagation // Ibid. 2013. Vol. 192. P. 163-189.

[6] Sahimi M., Tajer S.E. Self-affine fractal distributions of bulk density, elastic moduli and seismic wave velocities of rock // Phys. Rewiew. E. 2005. Vol. 71. 046301.

[7] Koohi lai Z., Vasheghani Farahani S., Jafari G.R. Non-Gaussianity effect of petrophysical quantities by using q-entropy and multi fractal random walk // Phys. A. Statist. Mech. and Appl. 2013. Vol. 392, iss. 20. P. 5132-5137.

[8] Soboleva O.N., Kurochkina E.P. Effective coefficients of quasi-steady Maxwell's equations with multiscale isotropic random conductivity // Ibid. 2011. Vol. 390, No. 2. P. 231-244.

[9] Кузьмин Г.А., Соболева О.Н. Подсеточное моделирование фильтрации в пористых автомодельных средах // ПМТФ. 2002. Т. 43, № 4. C. 115-126.

[10] Kolmogorov A.N. A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number //J. Fluid Mech. 1962. Vol. 13. P. 82-85.

[11] Гнеденко Б.В., Колмогоров А.Н. Предельные распределения для сумм независимых случайных величин. Л.: Гостехиздат, 1954.

[12] Ogorodnikov V.A., Prigarin S.M. Numerical Modeling of Random Processes and Fields: Algorithms and Applications. Utrecht, The Netherlands, 1996.

Поступила в 'редакцию 18 июня 2014 г., с доработки — 9 июля 2014 г.

Effective coefficients of a multiscale isotropic medium applied for the problem of propagation of acoustic waves

Kurochkina Ekaterina P.2, Soboleva Olga N.1

Wave propagation in heterogeneous media is a fundamental phenomenon of great scientific and practical interest. It is relevant to such important problems as detecting underground nuclear explosions, understanding the scale structure of oil, gas, and geothermal reservoirs. The spatial geometry of small-scale heterogeneities of medium is not exactly known. It has been shown that the irregularity of elastic parameters, density, permeability, porosity, increases as the scale of measurements decreases. It is customary to assume these parameters to be random fields characterized by the joint probability distribution functions. However, it is difficult to measure higher-order statistical moments for the geophysical parameters. At best, only the mean values and correlation functions of the second order are known. Hence, effective solutions cannot be constructed using the conventional perturbation theory. Geophysical parameters, for example, porosity, density, elastic modules can be well approximated by fractals and multiplicative hierarchical cascade models with lognormal distribution.

As the first step towards the goal of finding effective coefficients in the problem of propagation of elastic waves in strongly heterogeneous solids, in this paper we study the propagation of acoustic waves in the same type of media in which local elastic parameters have essentially all variations of scales from a certain interval at each spatial point. The density of a medium and the elastic stiffness are approximated by a multiplicative continuous cascade. We use the subgrid modeling method and obtain effective coefficients for the estimation of the first statistical moment of the displacement in the acoustic equation if wavelength essentially exceeds then a maximum scale of heterogeneity. If a medium is assumed to satisfy the improved Kolmogorov similarity hypothesis, the effective coefficients take especially the simple form. Differential equations for obtaining effective coefficients are also derived for the media that do not satisfy the improved similarity hypothesis. The derived formulas are verified by direct 2D numerical modeling.

Keywords: propagation of waves, acoustic equation, effective coefficients, subgrid modeling, multiscale random media, multiplicative continuous cascade.

Received 18 June 2014 Received in revised form 9 July 2014

1 Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk, 630090, Russia

2Institute of Thermophysics SB RAS, 630090, Novosibirsk, Russia Corresponding author: Soboleva Olga N., e-mail: olga@nmsf.sscc.ru