Модельное исследование и анализ закономерностей формирования солянокупольного бассейна Текст научной статьи по специальности «Физика»

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

Том 2, № 6, 1997

МОДЕЛЬНОЕ ИССЛЕДОВАНИЕ И АНАЛИЗ ЗАКОНОМЕРНОСТЕЙ ФОРМИРОВАНИЯ СОЛЯНОКУПОЛЬНОГО БАССЕЙНА

Е. Б. ОсиповА Институт механики и машиноведения МН—АН РК Алматы, Казахстан

A model research of united mechanism of arising and forming of structures of shaped like a dome as a result of "creeping"motins in the limits of Rayleigh — Taylor's hydrodynamic instability in two layer medium. This medium models an existing formation viz the oversalt sedimentary rocks — a salt layer.

1. Введение

Рассматривается весомое многослойное с выраженной инверсией по плотности и вязкости полупространство, которое описывается уравнениями ползущего течения, или уравнениями движения вязкой несжимаемой жидкости в приближении Стокса. Гидродинамическая неустойчивость, развивающаяся под действием силы тяжести между слоями с инвертированным распределением плотностей, называется неустойчивостью Релея — Тэйлора. Начальная линейная стадия, в течение которой возмущения с малой амплитудой экспоненциально возрастают во времени, достаточно полно изучена [2-8]. Нелинейные эффекты развития неустойчивости Релея — Тэйлора частично исследованы в [9-12]. Линеаризация основных уравнений и граничных условий допустима в случае "ползущего"медленного движения вязкой жидкости и малости амплитуды возмущения по сравнению с длиной волны. По мере нарастания возмущения наблюдается нарушение линеаризованных и действительных для начала развития неустойчивости граничных условий и усиление нелинейных эффектов.

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

Рассмотрим в трехмерной постановке в полупространстве П неустановившееся движение двухслойной сильно вязкой несжимаемой жидкости с общей поверхностью раздела £ (г, ^ = 0. Оба слоя являются несжимаемыми и однородными с динамической вязкостью

и ^2, плотностью р1 и р2 (Р1 > Р2) соответственно. Состояние равновесия с плоской (невозмущенной) поверхностью раздела примем за основное. Сила тяжести направлена от более тяжелого (верхнего) слоя к более легкому (нижнему) слою жидкости.

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

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

Зададим следующие граничные (кинематические и динамические) условия: на верхней поверхности осадочного слоя (1 (г, £) = 0 компоненты вектора скоростей ограничены и отсутствуют напряжения;

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

непрерывность компонент вектора скоростей,

И) И) И)

непрерывность компонент £13, £2з , £зз тензора напряжений с учетом поверхностного натяжения а и кривизны поверхности раздела Я,

изменение наклона поверхности раздела определяется из условия равенства нулю полной лагранжевой производной по времени; на поверхности (г, £) = 0 жидкость прилипает к твердой неподвижной поверхности, скорость обращается в нуль.

Рис. 1.

Расчетная схема модели приведена на рис. 1. Начальные условия определяют при £ = 0 состояние покоя.

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

3. Линейная стадия

Рассмотрим основной процесс в круговой цилиндрической системе координат:

х = р сое у = р вт г = г,

0 < р< то, 0 < ^ < 2п, —то < г < то.

(1)

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

А*(?)2 = д(Л2 — т (2)

получим в виде

V? = Л (?)

п

— ( 4Лел*(^ + А^е-^*) + (сЦЛел(^ + С^е-^

+

+^+1 (А*(^)р) |а*0) (А?^*+А^е-^*) +АО) (в^е^* — В^е-л(^

сов п^ ■ е

„шЬ

V? = { ¿> (?) п

А^е^ + А" е

А(?) (в?ел*0)'г — В^--л*00*

— ( С^е^* + СГе

+

1 е — В2 е " ~ ) + А

«С?) (с(°)ел(^)* + с2О)е-л(')г

вт п^ ■ е

шЬ

V,

(?)

{—А*? ■ ^ (А*(?)р^ е

1 (?)ел*(^)* — +

+ (В^е^ + В^ е

сов п^ ■ е

_шЬ

(3)

В заданной системе координат, подставляя соответствующие выражения для компонент вектора скорости и выражения компонент тензора напряжений в начальные и граничные условия, получим систему однородных линейных уравнений относительно неизвестных коэффициентов решения. Порядок системы зависит от количества рассматриваемых слоев и свойств симметрии решения. Систему решаем матричным способом в системе аналитических вычислений "МАТНЕМАТ1СА".

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

Полученная система имеет ненулевое решение, если определитель, составленный из коэффициентов при неизвестных, равен нулю. Приняв определитель равным нулю, получим характеристическое уравнение для определения функции скорости роста волнового возмущения ш = ^(А*(?), А(О), ,р?, Л?. Из аналитического решения следует, что скорость роста возмущений зависит от неоднородности физико-механических (д? — вязкости, р? —

р

плотности) и геометрических параметров (мощности слоев) рассматриваемой среды. При наличии поверхностного натяжения которое оказывает общее стабилизирующее действие, скорость роста возмущений равна нулю А*. =0 и при А*. = то; при промежуточных значениях из интервала 0 < А*. < то скорость роста принимает конечные значения. Существует хотя бы одно значение Ат., при котором скорость роста достигает максимума. В общем случае начальное возмущение содержит все волновые числа из интервала 0 < А*. < то , но так как для решения применим принцип суперпозиции, то конечная картина определяется максимальной скоростью роста и соответствующим вол-

л *(.)

новым значением Атах, которое определяется из условия

^ 0. (4)

Характеристическое уравнение, соответствующее решению (3), для круговой цилиндрической области, верхняя и нижняя поверхности которой ограничены плоскостями г = Ъ1 и г = Ъ2 = то, поверхность раздела слоев задана уравнением г = 0, может быть представлено полиномом пятого порядка и разрешено явно в системе "МАТНЕМАТ1СА". В этом случае в выражениях (3) А12) = В(2) = С12) = 0. Решение имеем в виде функциональной зависимости скорости роста волнового возмущения от всех физико-механических параметров среды.

Так как значение А*.) определяется доминирующим куполом по условию оптимальности скорости роста возмущений то вторичные купола, развиваясь при оседании края первого купола, должны расти медленнее. Конкретная зависимость определяется для каждого множества физико-механических параметров. Параметр п должен быть положительным целым числом.

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

д^

дА*(.)

д^

г

0,

0. (5)

дп

Полученные выражения компонент скорости, поверхности раздела можем физически интерпретировать как возмущение исследуемого параметра вследствие начального возмущения /0(А*(.)р), на которое накладывается вторичное возмущение вида /га(А*(.)р). Какое именно п определяет фундаментальную моду с максимальной скоростью роста вдоль первого кругового хребта при р = 0 можно выяснить, анализируя систему (5), положительные действительные нули и экстремумы функции Бесселя /га(А*(.)р).

Решение системы (5) проводим числено в системе "МАТНЕМАТ1СА"при п = 0,1, 2, ... , 57. Первый корень р = 0 соответствует куполу в центральной части области, второй корень х2 = А*(:,)р = 3.8317 дает минимум функции /0(А*(:,)р) и определяет радиус первой впадины. Третий корень = А*(:,)р = 7.0156 дает максимум функции /0(А*(^)р).

Расчеты показывают, что при п = 5 (или п = 6), п = 24, п = 27 влияние вторичного возмущения максимально, значит, вокруг первичного купола вероятнее образование конфигурации с пятью или шестью вторичными куполами по первому кругу радиуса р = ж3/А*(.), затем вероятны последующие конфигурации с 24 и 27 куполами.

4. Нелинейная стадия

Линеаризация уравнений Навье — Стокса, приводящая к уравнениям в приближении Сток-са, допустима в случае "ползущего"медленного движения вязкой жидкости и малости амплитуды возмущения по сравнению с длиной волны. По мере нарастания возмущения наблюдается нарушение всех определенных и действительных для начала развития неустойчивости граничных условий. В частности:

соответствующие компоненты напряжения не согласуются с осями заданной системы координат, а варьируются в соответствии с наклоном поверхности раздела;

непрерывность компонент напряжения и скорости следует определять по действительному (возмущенному) £ *(г, ¿) = 0, а не начальному уровню поверхности раздела £ (г, ¿) = 0;

на поздних стадиях развития неустойчивости условие непрерывности скорости вдоль поверхности раздела следует заменить условием скольжения.

Очевидно, что учет всех этих условий значительно усложнит задачу и соответствующий анализ. Следуя предложенному в [3] анализу нелинейной стадии развития неустойчивости, предполагаем, что по мере развития основного процесса:

на возмущенной поверхности раздела слоев среды £ * (г, ¿) = 0 отсутствует явление диффузии, скольжения;

две точки различных слоев, находящиеся на бесконечно малом расстоянии, остаются на этом расстоянии и их скорости отличаются на пренебрежимо малую величину;

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

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

£ [г - С *(г, *)]

I + V-

[г - С*(М)]. (6)

В круговой цилиндрической системе координат, заданной согласно (1), выполняя дифференцирование в выражении (6) и учитывая, что полная лагранжева производная по времени уравнения поверхности раздела слоев равна нулю, получаем соотношение

= дС*(р,М , дс*(р,М , 1 дс*(р,М (7)

^ = ^ + р др + р^ д^ . (7) Непосредственное решение нелинейного уравнения (7) аналитическими методами затруднительно. Следуя приближенному методу решения согласно [3], предполагаем, что:

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

по мере развития неустойчивости и ростом куполов нелинейные члены становятся доминирующими.

Согласно вышеприведенным предположениям на линейной стадии из соотношения (7) имеем:

*

^ = ж

Тогда при г = £ * первая аппроксимация нелинейной части выражения (7) пре

образуется к виду

(8)

ас * (p,v,t) , 1 ас * (p,v,t)

dp

+ - V

P

Jn (Л*p) - [- (AieA*z* + A2e-A*z*) +

z=Z * (p,<p,t) ^ p

+ (CieAZ* + C2e-AH] + Jn+1 (Л*р) [Л* (A^ + A2e-A*z*) + Л (BieAZ* - £2e-^)] }

Л* " "

x < -

и

n

p

Jn (Л* p) - Л* Jn+i (Л* p)

x

[(AieA*z* -A2e-A*z*) + (BieAZ* + £2e-AZ*)] \ cos2 • e2wt+

1

+ - \ Jn (Л*p) — [(AieA*z* + A2e-^^ - (CieAC* + C2e-A^] +

pp

+Jn+i (Л*p) [Л (BieA*z* - B2e-A*z*) + Л* (CieAZ* + C2e-AZ*)]

x -{Л* Jn (Tp) [(Aiev^ - A2e-A^) + (£ieAZ* + £2e-AC*)] } sin2 • e2wt. (9) и

5. Анализ неустойчивости поверхности раздела системы "надсолевые осадочные породы — соляное тело"

Исследование особенностей формирования купольных структур выполнено применительно к данным обстановки Прикаспийской впадины:

вязкость соли д2 = 2.54 • 10i8 пуаз и вязкость верхнего осадочного слоя = 1020 пуаз; плотность соли p2 = 2.16 г/см3 и плотность верхнего осадочного слоя pi = 2.65 г/см3; мощность верхнего осадочного слоя hi = 1.4 • 105 см; ускорение свободного падения g = 981.6 см/с2.

Рис. 2.

Как показывает анализ полученных решений, физико-механическому смыслу задачи подходит одно решение, которое дает непрерывный положительный спектр параметра ы. На рис. 2 показано это решение, определяющее зависимость функции волнового возмущения ы = ы(А*) (1/с) от волнового числа А* = д • 10-6 (1/см).

V

р

X

Решение и исследование системы (4) позволяет определить значение параметра Л*, обеспечивающее максимальную скорость роста возмущений и соответственно доминирующую моду развития неустойчивости: Л* = 0.00007607 1/см, ш = 0.00003052 1/с, Л« = 0.00007607 1/см, Л(2) = 0.00007607000000000016 1/см.

В системе "МАТНЕМАТ1СА" определено в общем виде, что ранг матрицы данной системы линейных однородных уравнений равен 8.

Фундаментальное решение соответствующей системы линейных однородных уравнений определено с точностью до постоянной

= ¿1, 41} = 1.441115142723703 • 1024 •

= -4.413583732956841 • , В2 =0. • , С^ = 1.622970873027692 • 1014 • , С21} = 2.888081685393142 • 1023 • ,

А22) = 0.5516979666196051 • , В22) = 2.871781624383739 • 1015 • С22) = -2.120176249007957 • 1016 • .

Непосредственная подстановка полученных значений в соответствующее выражение нормальной компоненты вектора скорости позволяет составить и решить уравнение относительно параметра по известным средним значениям скорости (10-3 см/год ) процесса в масштабах геологического времени:

Ж = -4.586 • 10

-423

Рис. 3.

Рис. 4.

На рис. 3 для сравнения приведены графики амплитудных значений на линейной (1) и нелинейной (2) стадиях сечения плоскостью х = 5.5 • 105 см поверхности раздела: надсо-левые осадочные породы — соляное тело в области ограниченной 0 < х < 5.5 • 105 см, 0 < у < 5.5 • 105 см, при £ = 1.0022 млн лет и £ = 1.0068 млн лет соответственно. На линейной стадии имеем увеличение амплитудных величин куполов и впадин в пределах, допускаемых линейной теорией. На рис. 4 приведен график равных амплитудных значений на линейной при £ = 1.0022 млн лет и нелинейной при £ = 1.0028 млн лет стадиях сечения плоскостью х = 5.5 • 105 см области Очевидна правильность предположения о пренебрежении нелинейными слагаемыми в уравнении (7) до тех пор, пока мал угол наклона

поверхности раздела в заданной точке. Нелинейные члены становятся доминирующими по мере развития неустойчивости и с ростом куполов.

На рис. 5 приведены графики амплитудных значений на линейной (1) при Ь = 1.0022 млн лет и нелинейной (2) стадиях при Ь = 1.0043 млн лет сечения плоскостью х = 5.5 • 105 см поверхности раздела: надсолевые осадочные породы — соляное тело в области Д^.

Рис. 5. Рис. 6.

На рис. 6 приведены графики амплитудных значений на нелинейной (1, 2, 3) стадии при Ь = 1.0048, 1.0053 и 1.0064 млн лет соответственно сечения плоскостью х = 5.5 • 105 см поверхности раздела: надсолевые осадочные породы — соляное тело в области Д. Усложненный характер взаимодействия функций Бесселя отражает сложный механизм образования и развития межкупольных структур. Очевидно, что Л* — компонента разложения "порождает" и "питает" все более высокие функции Бесселя, т. е. купольные поднятия "питаются" с близлежаших площадей соседних куполов и "подпитка" распределена неравномерно.

Общий анализ соотношения (7) с учетом нелинейной части в виде (9) позволяет выявить некоторые характерные для всего процесса в целом особенности:

наличие выражения возмущенной поверхности раздела в показателе экспоненциальной функции определяет асимметрию роста относительно начальной поверхности;

наличие временного множителя во второй степени показывает, что вклад нелинейных слагаемых на ранних стадиях развития неустойчивости пренебрежимо мал, но по мере развития процесса возрастает и должен преобладать;

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

На рис. 7 приведены фрагменты графиков поверхности раздела: надсолевые осадочные породы — соляное тело на нелинейной стадии при Ь = 1.0058 млн лет областей Д, ограниченной 0 < х < 5.5 • 105 см, 0 < у < 5.5 • 105 см (рис. 7, а); Д2, ограниченной 5.5 • 105 < х < 1.2 • 106 см, 0 < у < 5.5 • 105 см (рис. 7, б); Д3, ограниченной 1.2 • 106 < х < 5.6 • 106 см, 0 < у < 5.5 • 105 см (рис. 7, в); Д4, ограниченной 9.9 • 106 < х < 4.37 • 107 см, 0 < у < 5.5 • 105 см (рис. 7, г).

На основе исследования неустойчивости Релея — Тэйлора поверхности раздела в кусочно-однородных (двухслойных) сильно вязких средах с инвертированным распределением плотности установлено: формирование солянокупольных структур в осадочных бассейнах, возможно, является гидродинамической неустойчивостью Релея — Тэйлора, развивающейся на поверхности раздела системы "надсолевые осадочные породы — соляное тело".

Рис. 7. График поверхности раздела "надсолевые осадочные породы — соляное тело" в области

Di (a), D2 (б), D3 (в), D4 (г).

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

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

[1] БЕйТМЕН Г., ЭРДЕйИ А. Высшие трансцендентные функции: Функции Бесселя, функции параболического цилиндра, ортогональные многочлены. Т. 2, Наука, М., 1974.

[2] Бигкгоф Г. Неустойчивость Гельмгольца и Тейлора. В "Гидродинам. неустойчивость", Мир, М., 1964, 68-94.

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

[4] Layzer D. On the instability of superposed fluids in a gravitational field. Astrophys. J., 122, No. 1, 1955, 1-12.

[5] Lewis D.I. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. II. Proc. R. Soc. London, Ser. A, 202, 1950, 81-96.

[6] Menikoff R., Mjolsness R. C., Sharp D. H., Zemach C. Unstable normal mode for Rayleigh —Taylor instability in viscous fluids. Phys. Fluids, 20, No. 12, 1977, 2000-2004.

[7] Menikoff R., Mjolsness R. C., Sharp D. H., Zemach C., Doyle B. J. Initial value problem for Rayleigh — Taylor instability of viscous fluids. Ibid., 21, No. 10, 1978, 16741680.

[8] Miles J.W. and Dienes J.K. Taylor instability in a viscous liquid. Ibid., 9, No. 12, 1966, 2518-2519.

[9] Nayfen A. H. On the non-linear Lamb —Taylor instability. J. Fluid Mech., 38, p. 3, 1969, 619-631.

[10] ОсиповА Е. Б. Релей-тэйлоровская неустойчивость поверхности раздела сильно вязких сред в круговой цилиндрической области. Деп. в КазгосИНТИ 05.02.97, №7416-Ка97, Алматы, 1997.

[11] Осипова Е.Б. Нелинейные эффекты развития неустойчивости Релея — Тэйлора. Там же, №7417-Ка97, Алматы, 1997.

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

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