Исследование фундаментальной моды неустойчивости Релея - Тэйлора в бассейне эллиптической формы Текст научной статьи по специальности «Физика»

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

Том 4, № 2, 1999

ИССЛЕДОВАНИЕ ФУНДАМЕНТАЛЬНОЙ МОДЫ НЕУСТОЙЧИВОСТИ РЕЛЕЯ-ТЭЙЛОРА В БАССЕЙНЕ ЭЛЛИПТИЧЕСКОЙ ФОРМЫ

Е. Б. ОсиповА Институт механики и машиноведения МН-АН РК Алматы, Казахстан e-mail: imms@kaznet.kz

The fundamental mode of Rayleigh — Taylor hydrodynamic instability in an elliptical basin has been numerically and graphically analyzed for the data of geological situation in the Caspian depression.

Введение

Неустойчивость Релея — Тэйлора представляет собой частный случай общего понятия гидродинамической неустойчивости слоистой среды с выраженной неоднородностью физических свойств, обусловленной градиентами плотности, температуры и давления [13]. Качественное проявление этого вида поверхностной неустойчивости, развивающейся в поле силы тяжести, в физических процессах общепризнанно. Естественное уплотнение осадочных пород приводит в некоторых регионах земной коры к инверсии плотностей, т. е. породы с большей плотностью оказываются налегающими на менее плотные. Уникальным в этом плане является соляно-купольный бассейн Прикаспийской впадины. По гравиметрическим данным, средняя плотность верхних осадочных слоев здесь составляет 2.5 - 2.9 г/см3, соляной толщи — 2.1 - 2.2 г/см3 [14]. На поверхности раздела системы перекрывающие осадочные толщи — соль осадочный чехол Прикаспийской впадины содержит разнообразные соляно-купольные структуры. По гипотезе гравитационной неустойчивости, формирование соляно-купольных структур в седиментационных областях осадочного чехла земной коры является следствием неустойчивости Релея — Тэйлора [2, 10, 15]. При этом деформирование толщ каменной соли и осадочных горных пород в геологическом масштабе времени рассматривается как следствие “ползущих11 движений сильно вязкой среды и описывается уравнениями движения вязкой несжимаемой жидкости в приближении Стокса [3, 5, 12]. В настоящей работе исследовано в линейном приближении развитие неустойчивости в бассейне эллиптической формы для данных геологической обстановки Прикаспийской впадины. Вывод, преобразование основных соотношений и численно-графический анализ характеристического уравнения с применением результантов выполнены в системе MATHEMATICA [1].

© Е. Б. Осипова, 1999.

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

В общем виде расчетная схема и трехмерная постановка задачи исследования закономерностей развития поверхностной неустойчивости Релея — Тэйлора приведены в [3, 11]. Отметим основные положения.

Рассмотрим “ползущее” движение двуслойной сильно вязкой несжимаемой жидкости, имеющей общую поверхность раздела £ (г,¿) = 0 и распространенной по полупространству П. Оба слоя являются несжимаемыми и однородными с динамической вязкостью р\ и р2, плотностью р1 и р2 > р2) соответственно. Сила тяжести направлена от более тяжелого

(верхнего) слоя к более легкому (нижнему) слою жидкости.

Конкретизация физического закона состояния в виде обобщенного закона Ньютона, условие несжимаемости, линеаризация уравнения движения в пространстве, арифмети-зированного эйлеровыми переменными, определяют основное уравнение движения в приближении Стокса.

Основной процесс рассматривается при граничных условиях:

на верхней поверхности г = Н\ осадочного слоя компоненты вектора скоростей ограничены и напряжения отсутствуют;

на поверхности раздела осадочного и соляного слоев £ (г, ¿) = 0 имеют место непре-

И) И) И)

рывность компонент вектора скоростей и непрерывность компонент ¿33 , ¿2з , ¿13 тензора напряжений.

Начальные условия определяют при ¿ = 0 состояние покоя. Основное состояние — равновесие с плоской (невозмущенной) поверхностью раздела. Поверхность раздела приведенной системы динамически неустойчива [13]. Определим в линейном приближении фундаментальную моду возмущенного состояния.

Общее решение для несжимаемых жидкостей можно выразить по теореме Гельмгольца, согласно которой векторное поле скоростей (у^,'уП'),'уР)) — однозначное, непрерывное и обращающееся в нуль на бесконечности, может быть представлено в виде суммы скалярного градиента р и ротора векторного потенциала Ф [8, 9]. Доказана единственность этого разложения. Исходной задачей будет:

V = —Ур + Ух Ф. (1)

Здесь массовые силы потенциальны и выражаются через скалярный потенциал и:

g = -Уи = —к(ди/дг) = — кд. (2)

Подставляя (1) в уравнение Стокса, с учетом условия несжимаемости и потенциальности массовых сил (2) после ряда преобразований и изменения порядка дифференцирования получаем систему, представляющую полное множество уравнений движения [3]:

У2и = 0,

У V = 0,

V У2Ф — (дФ/&) = 0,

(др/эд = р/р + и, (3)

где V — вектор скорости, Р — давление, р — плотность, g — сила тяжести, vр = р — коэффициент динамической вязкости, V — коэффициент кинематической вязкости, ¿ — время, индекс ] (] = 1, 2) означает принадлежность величины к ^-му слою.

Обобщение понятия разделимости на векторные поля и соответствующие решения с учетом свойства соленоидальности допускает в векторном уравнении (3) разделение на три независимых скалярных уравнения даже в криволинейных системах координат. Отделением временной зависимости в виде экспоненциального множителя получаем векторное уравнение Гельмгольца для зависящей от пространственных координат части решения. Разложение имеет вид [8, 9]

Ф = Уф + Ух (а^ф) + (1/к)У х У х (а^х), (4)

где ф, ф, х — собственные функции скалярных уравнений типа

У2ф + к2ф = 0. (5)

Задача состоит в выборе к и соответствующих функций ф, ф, х таким образом, чтобы

удовлетворялись все заданные граничные условия. Множитель 1/к уравнивает размерность функций ф, шф, шх.

2. Общее решение в линейном приближении

Рассмотрим основной процесс в эллиптической цилиндрической системе координат O^Vz [8]:

x = Cch £ cos n, y = Csh £ sin n, z = z,

0 < £ < to, 0 < n < 2n, C > 0 (6)

с метрическими коэффициентами

Hi = H2 = CV(ch 2£ - cos 2n)/2, Я3 = 1.

Решение (5) определяется методом разделения переменных. В эллиптической цилиндрической системе координат имеем решение, частично разделяющее константы разделения, в виде соответствующих рядов по собственным функциям [8, 9]. Рассмотрим частное решение системы (3), удовлетворяющее заданным граничным и начальным условиям, в виде функций Матье первого рода целого порядка. Граничные условия определяют для функции, зависимой от координаты n, свойство периодичности T = п; для функции, зависимой от £, — свойство конечности и непрерывности при £ = 0; для функции, зависимой

от z, — свойство конечности при z ^ то. Частным решением системы (3) для каждого слоя среды с учетом граничных и начальных условий имеем:

[- (A(/¿exp(A*(j)z) + A2mexp(-A*(j)z)) + (A(j)/fc(j))x

m=0

x (cJ£ exp(A(j)z) — C2mm exp(-A(j)z)) Sem(d(j), cos n) [Jem(d(j), ch £ )]^ +

+ ^B(mlexp(A(j)z) + E»exp(—A(j)z)) [Som(d(j), cos n)]'v Jom(d(j), ch£) > exp(wí),

m=1 )

[— (A?„mexp(A*(j)z) + A2m exp(—A*(j)z)) + (A(j)/k(j))x

m=0

х (Cmexp(A(j)z) - Cjexp(—A(j)z)] [Sem(dj), cos n)]l Jem(d(j), ch£) —

— iB(mmexp(A(j)z) + B2mlexp(—A(j)zH Som(d(j), cos n) [Jom(d(j), ch£)]^ > exp(wt),

m= 1

,j)

m=0

х (Cexp(A(j)z) + C^miexp(—A(j)z)) Sem(d(j), cos n)Jem(d(j), ch £) !• exp(wt)

^ I—A*(j) (A(imiexp(A*(j)z) — Ajexp(—A*(j)z)) + (A(j)2 + fc(j)2) /k(j)

х

(7)

где

d(j) = C\/ A(j)2 + k(j)2, k(j) = —upj /pj.

В выражении (7) обозначены функции Матье первого рода целого порядка m (m > 0) Jom(d(j), ch £), Som(d(j), cos n), Jem(d(j), ch £), Sem(d(j), cos n), определенные, согласно [7, 8], в виде соответствующих рядов по тригонометрическим функциям и функциям Бесселя первого рода действительного аргумента целого порядка n.

Подставляя соответствующие выражения компонент скорости в соотношения компонент тензора напряжений и тензора скоростей деформаций, а затем в заданные граничные условия, получаем однородную линейную систему из 9 уравнений относительно неизвестных 4 m, A2m, в 1 m, Bm, с mi c2m(j = i> 2 =в (m=cs=0). Полученная система

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

w(d(j),^j ,pj ,hj) = 0.

(8)

В общем виде имеем нелинейное трансцендентное уравнение, которое можно упростить, не нарушая общности рассуждений, и решить, если принять значение 1

одинаковым для ] слоев среды. После отделения нулевых решений уравнение (8) преобразуется к виду

Ао^ + А !^ + А.2^ + Аэ^ + ^4^ + А5 = 0, (9)

где А = А*(1,ру ,ру, Ну).

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

1тах (наименьший параметрический нуль), которое определяется численно из системы

w(d,pj ,pj ,hj) = ° w(d,pj ,pj ,hj )d = °.

(10)

В случае решения задачи в эллиптической цилиндрической системе координат (6) функция , Ру, Ну) проверяется на максимум через вторые производные по общим прави-

лам.

СО

оо

z

Соответствующее выражение поверхности раздела можем физически интерпретировать как возмущение исследуемого параметра вследствие начального возмущения, на которое накладывается вторичное возмущение. Очевидно, что в общем случае значение двух параметров 1тах и т определяет пространственное расположение куполов и, как следствие, согласно [12], расстояние между смежными гребнями определяется расположением максимумов функции Бесселя Лп(1 • еЬ£). Параметр 1 изменяется непрерывно, п — принимает только целые значения. Но так как в общем случае функция Бесселя Лп(1 • еЬ £) является непрерывной функцией порядка п для п > 1, то рассмотрим п и соответственно т как непрерывные переменные. При этом допущении максимальная скорость как функция непрерывных переменных 1 и т определится из следующей системы:

( ду^ = 0

= ’ (11)

| дуУ) = 0 ^ дт ’

Фундаментальную моду в начальном возмущении с максимальной скоростью роста и соответствующие значения параметров т и п получим, анализируя систему (11), положительные действительные нули и экстремумы функции Бесселя Лп(1 • еЬ £). Для этого подставим в (11) выражение соответствующей компоненты вектора скорости (7) для

элементов ]-го слоя при £ = 0 и исследуем численно поведение точек при п = 0. В раз-

вернутом виде систему (11) не приводим из-за громоздкости. При дифференцировании функций Бесселя по порядку использована формула [12]

П — 1

дЛ„(ж)/дп = (п/2)Уга(ж) + (п!/2) ^(ж/2)к—пп [Л*.(ж)/(&!(& - п))]. (12)

к=0

Дифференцирование функций Бесселя по аргументу выполнено по общим правилам [8].

3. Численно-графический анализ

Решение характеристического уравнения (8) и определение соответствующих параметров систем (10), (11) выполнено для входных данных бассейна Прикаспийской впадины [5, 14]: вязкость соли р2 = 2.54 • 1018 пуаз и верхнего осадочного слоя р1 = 1020 пуаз; плотность соли р2 = 2.16 г/смэ и верхнего осадочного слоя р1 = 2.65 г/смэ; мощность верхнего осадочного слоя Н1 = 1.4 • 105 см; ускорение свободного падения д = 981.6 см/с2.

По “Структурной карте надсолевого комплекса Прикаспийской впадины” (М 1: 1000000, 1980 г., ред. Л. Ф. Волчегурский, О. С. Турков, А. Е. Шлезингер; Мин. геол. СССР, Всесоюзное аэрогеолог. научно-производ. объед. “Аэрогеология”) определено расположение соляно-купольной области:

по долготе 440.0 - 570.9 в.д. = 130.9, по широте 450.8 - 510.8 с.ш. = 60.0.

Аппроксимируем рассматриваемую область эллипсом, отношение малой полуоси которого = Ь к большой = а (С2 = а2 — Ь2) равно 5/6; общая площадь области составляет 600 000 кв. км. Имеем:

а = 478.731 км, Ь = 398.942 км, = 264.628 км.

Соответствие расчетного эллиптического контура и границы Прикаспийской НГП установлено по обзорно-тектонической карте “Прикаспийская впадина и прилегающие районы” (М 1: 1000000, 1989 г., ред. Т. А. Акишев, Э. С. Воцалевский, Ю. С. Кононов, С.К. Курма-нов, О. С. Турков, С. У. Утегалиев, Д. Л. Федоров, Мин. геол. СССР, Всесоюзное аэрогеолог. научно-производ. объед. “Аэрогеология”).

При указанных данных система (10) сводится к двум полиномам пятого порядка относительно переменной ш, коэффициенты которых являются функциями параметра 1. Основное алгебраическое решение состоит в исключении переменной ш через результант. Решение системы многочленов (10) сводится к определению корней одного уравнения относительно параметра 1, так как результант принадлежит кольцу коэффициентов рассматриваемых полиномов и равен нулю тогда и только тогда, когда эти полиномы имеют общий корень [6]. В данном случае результантом И,в2 = Д(1) будет определитель десятого порядка матрицы Сильвестра [4], составленной из коэффициентов полиномов системы

(10). В системе МАТНЕМАТ1СА [1] с помощью Подфункции выполнены определение начальных приближений и численное решение полученного уравнения относительно параметра 1 в формате Е1^И,оо1;-оператора. Фрагмент значений функции Я(1) и параметра 1, по данным Р1о1-функции, приведен в таблице, соответствующий график — на рисунке.

Т а б л и ц а

а ща)

0.001177382812499999 0.001179192708333333 0.001181002604166666 0.001182812499999999 0.001184622395833333 0.001186432291666666 0.001188242187499999 0.001190052083333333 0.001191861979166666 0.001193671874999999 0.001224440104166666 0.00122625 0.001228059895833333 9.05805751050194Е — 17 3.103725686987218Е — 16 1.242275593451Е — 16 —(2.069851730685726Е — 17) 2.814565612540769Е — 16 3.998260308097988Е — 16 4.517947612154636Е — 16 2.742328657565452Е — 16 4.744527734405079Е — 16 6.044940228284266Е — 16 2.498786899069588Е — 16 —(7.056608984262915Е — 17) 1.371107094501782Е — 15

Берем такое расчетное значение параметра 1, которое является корнем функции Д(1), многочленов системы (10), и при этом старшие коэффициенты этих многочленов не равны нулю. Имеем 1 = 0.001182812499998678.

Решение системы (11) проводим с учетом (12) численно при 1 = 0.001182812499998678 и т =1, 2,..., 57 в системе МАТНЕМАТ1СА [11]. Аналогично случаю круговой цилиндрической системы координат первый корень соответствует куполу в центральной части области, второй корень определяет положение первой впадины. При т = 5 (возможно, т = 6)

Рис. 1. График функции Rez = R(d).

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

Vz = dZ */ôi. (13)

Подставляя в выражение (13) vz, согласно (7), при условии Л*(1) = Л*(2) = Л*, Л(1) = Л(2) = Л

и интегрируя его по t, определим уравнение поверхности раздела:

У [-Л* (Aimexp^*z) - А2ТОвхр(-Л*г)) + [(Л2 + fc2)/fc] х

m=0

х (С1твхр(Лг) + C2mexp(-Лг))] Sem(d, cos n)Jem(d, ch £)} exp(wt). (14)

Соотношение (14) при указанных параметрах d и m определяет фундаментальную моду развития неустойчивости Релея — Тэйлора. При этом в начале развития возмущения имеем увеличение амплитудных величин куполов и впадин в пределах, допустимых линейной теорией [13]. Усложненного взаимодействия доминирующих и вторичных мод в плане слияния куполов и впадин не наблюдается [11].

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

[1] АлАДЪЕВ В.З., ШИШАКОВ М. Л. Введение в среду пакета MATHEMATICA 2.2. Инф.-изд. дом “Филинъ”, М., 1997.

[2] ARRHENIUS Sv. Sur Physik der Salzlagerstatten. Meddel. K. Vetenskabsakademiens, Nobelinstitut, 2, No. 20, 1912.

[3] DANES Z. F. Mathematical formulation of salt — dome dynamics. Geophysics, XXIX, No. 3, 1964, 414-424.

[4] ДЭВЕНПОРТ Дж., Сирэ И., Турнье Э. Компьютерная алгебра. Мир, М., 1991.

[5] ЕРЖАНОВ Ж. С., БЕРГМАН Э. Н. Ползучесть соляных пород. Наука, Алма-Ата, 1977.

[6] Курош А. Г. Курс высшей алгебры. Наука, М., 1971.

[7] МАК-ЛАХЛАН Н. В. Теория и приложения функций Матье. ИЛ, М., 1953.

[8] МОРС Ф. М., ФЕШБАХ Г. Методы теоретической физики. Т. 1, ИЛ, М., 1958.

[9] МОРС Ф. М., ФЕШБАХ Г. Методы теоретической физики. Т. 2, ИЛ, М., 1960.

[10] NETTLETON L. L. Fluid mechanics of salt-domes. Bull. Amer. Ass. Petrol. Geol., 18, No. 9, 1934, 1175-1204.

[11] ОСИПОВА Е. Б. Модельное исследование и анализ закономерностей формирования соляно-купольного бассейна. Вычислительные технологии, 2, №6, 1997, 61-70.

[12] SELIG F. A theoretical prediction of salt dome patterns. Geophysics, XXX, No. 4, 1965, 633-643.

[13] TAYLOR G. I. The instability of liquid surfaces when accelerated in a direction

perpendicular to their planes I. Proc. Royal Soc. London, Ser. A., 201, No. 1065, 1950,

192-196.

[14] Тектоническое строение и перспективы нефтегазоносности Прикаспийской впадины. Под ред. Н. П. Казакова, Н. Н. Чарыгина. Гостоптехиздат, М., 1958.

[15] TRUSHEIM F. Mechanism of salt migration in northern Germany. Amer. Assoc. Petrol. Geol. Bull., 44, No. 9, 1960, 1519-1540.

Поступила в редакцию 3 июня 1998 г.