CC BY
71
ДОКЛАДЫ АКАДЕМИИ НАУК РЕСПУБЛИКИ ТАДЖИКИСТАН ________________________________________2011, том 54, №9______________________________________
МЕХАНИКА
УДК 532.546.:550.820.7
Н.К.Корсакова, В.И.Пеньковский, М.А.Саттаров* МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПРОЦЕССОВ АБЛЯЦИИ ДИАПИРОВ
Институт гидродинамики им. М.АЛаврентьева СО РАН, Новосибирск,
Институт водных проблем, гидроэнергетики и экологии АН Республики Таджикистан
(Представлено академиком АН Республики Таджикистан З.Д.Усмановым 03.10.2011 г.)
Предлагается математическая модель процесса уноса массы с поверхности непроницаемых массивов каменной соли (например, диапиров), замкнутых включений мерзлоты или айсбергов, обтекаемых потенциальным потоком воды. Рассматриваются случаи двумерного потенциального обтекания, для которых можно выписать точные аналитические решения. В общем случае применяется численный метод конечных элементов. Создана экспериментальная установка и проведены предварительные опыты, результаты которых качественно согласуются с расчётами.
Ключевые слова: каменная соль - диапир - мерзлота - фильтрация - обтекание - модель - эксперимент - расчёт - плотина - гидроузел.
1. Постановка задачи. Рассмотрим плоскую задачу, соответствующую плановому потоку грунтовых вод, обтекающему одиночный диапир, газогидрат, или плоско-вертикальному потоку, обтекающему ленточный массив мерзлоты. Пусть: 1) поток потенциальный, то есть существует аналитическая функция /(г) = р(х, у, + I /(х, у, ^) такая, что
и = др/дх = д// ду, V = др/ду = -д// дх, (1)
где и,V - компоненты скорости, р,/ — сопряжённые гармонические функции, х,у - координаты, ?
- время, г = х + 1у - комплексная переменная;
2) нормальная составляющая скорости перемещения точек контура у пропорциональна проекции ^ скорости потока на касательную к контуру у
Vnc =-яЫ, (2)
где а — постоянная, определяемая экспериментально.
Пусть уравнение Е (х, у; I) = 0 представляет изменяющийся во времени обтекаемый контур у. Представляя его полный дифференциал по времени ?,
дЕ / д^ + / Л = 0
х с у •у с
Адрес для корреспонденции: Саттаров Малик Абдусатторович. 734002, Республика Таджикистан, г.Душанбе, ул.Парвин, 12, Институт водных проблем, гидроэнергетики и экологии АН РТ. E-mail: [email protected]
и обозначая п0 внешнюю по отношению к контуру нормаль, получим
Vnc = Vc К = (FxdXc /dt + Fydyc / dt)/ Fn , Fn = >/ FX + Fy ; Vnc = -Ft / Fn
или, принимая во внимание соотношение (2),
Ft / Fn = a\vc\ = a<p\ • (3)
Отметим, что скорость фильтрации отличается от истинной скорости движения частиц на величину пористости. В процессах оттаивания мерзлоты или растворения массивов в пористых средах этот факт может быть учтён с помощью специально подобранного коэффициента кинетики а.
Нормальная составляющая относительной скорости vrn потока выражается в виде vrn = др/ dn — vnc, а поток массы жидкости qo сквозь контур формулой q0 = р0 (<рп — vnc) . Здесь р0
- плотность жидкости. Приток qx массы жидкости, обусловленный фазовым переходом в процессе оттаивания массива льда (или разложения газогидрата) с массовым содержанием воды р , вычисляется в виде q1 = P1vnc. В случае растворения массива соли: р1 = 0 . В общем случае условие непроницаемости контура запишем так: q0 + q = 0 или
Pn =svnc ,S = 1 — Pl/ Po . (4)
Учитывая соотношения Римана (1) в натуральной системе координат (рп = щ s;ps = —щп, нахождение аналитической функции f (z) и формы контура F(x, y, t) сводится к решению начально-краевой задачи
дщ/ds = еа|др/д^,z = zk е у
Ft /Fn = а\др/&|,F(zk;0) = Fo(zk) (5)
lim df / dz = , z ^го.
Здесь иш обозначает вектор скорости течения на бесконечности. Рассмотрим некоторые тривиальные примеры, допускающие аналитические решения.
2. Обтекание круга циркуляционным потоком жидкости. Пусть F = x2 + у2 — R2 (t) = 0 -
искомое уравнение контура кругового массива (R(0) = 1), обтекаемого циркуляционным потоком с циркуляцией Г = const. Комплексный потенциал потока ищем в виде вихреисточника
—Г + tN N — Г 0 Г N0 i&
f (z) =-ln z, ю = — lnr л-----------------------, щ =— lnr Л-, (z = re ).
2^i 2ж 2ж 2ж 2ж
Поскольку элемент дуги контура ds = —R(t)d0, то для касательной к контуру скорости получаем V = (ps = (P@®s = Г /(2^R(t)). Первое условие задачи (5) дает щ s = —N /(2яК) = аеГ /(2яК)., откуда N = —ааГ . Из второго условия (5) имеем
р / Р = —dR / Ш = аГ/2Ж,
откуда с учётом начального условия следует R(t) = у/1 — аГ/ж . Массив исчезает, когда I = Т = ж /(аГ), а - определятся из опыта.
3. Обтекание бесконечного массива в режиме установления. Пусть существует некоторая форма симметричного относительно оси х массива, образуемого в процессе длительного его обтекания агрессивной жидкостью. Будем искать уравнение контура массива в виде Р(х, у, t) = у — Ас (х — vct) = 0 и рассмотрим область течения, ограниченную линиями у = 0, х е (—*, А) и у > fc (х - vct), х е (А, да) . Поскольку р > 0, то первое условие задачи (5) можно проинтегрировать. Без ограничения общности произвольную константу считаем равной нулю, и для А1 =р1 +1/1 = f (г1) получим краевое условие
/1 = 0:у = у = 0х1 = х — Vkt < 0;/1 =а/р1:у1 = Л(х1). (6)
Таким образом, область комплексного потенциала А (г) составляет внутренность угла (1 — а)ж между лучом р е (—да,0),/ = 0 и прямой линией / = шру, которая наклонена под углом жа = аг^(аа) к полуоси р1 > 0.
Отобразим конформно область потенциала А на верхнюю полуплоскость т] > 0 вспомогательной комплексной переменной £ = £ + 1т так, чтобы точка А = 0 перешла в точку £ = 0, а бесконечно удалённая точка А = да - в бесконечно удалённую точку £ = да. Потенциал w(£) имеет вид
w(£) = в1жа£1-а . (7)
Пусть £ = g(^ + /^2 представляет конформное отображение области движения в плоскости г на верхнюю полуплоскость вспомогательной переменной £ с соответствием точек g(0) = 0, g(да) = да . Тогда уравнению контура Р = 0 в плоскости £ будет соответствовать уравнение Р = 1т g = ^ = 0; т = 0, £ е (0, да). Далее получаем выражения для производных Р( = —Vс 1тdg/Лг1,Рп = |dg/= 0,£ > 0 и для второго условия в задаче (5)
— V Im(dg/ Лг) /|dg/ Щ = a|dg/ Яе dw / Л£; е у,т = 0, £ > 0. (8)
Пусть г = г (£) аналитическая функция, обратная к функции £ = g(). Тогда
dg / dzx = 1/(dz / d£,lm(dg / dzx)/ dg / dz = Im dzl / d£) = — Im(dz / d£ и условие (8) с учетом формулы (7) преобразуется к виду
Im(dz: /d£) = аcosжal(1 — а)£~а /vk,ц = 0,£> 0. (9)
На линии симметрии потока Г] = 0, £ < 0 имеем Im zx = 0 или, дифференцируя, Im(dzx / d£) = 0 . По условиям на действительной оси — да < £ < да производная искомой функции zx (£) легко восстанавливается:
dzx / d£ = (а/ vk )^жа(1 — а)егжа£~а . (10)
Удовлетворяя условию на бесконечности £ ^*;lim(( dw / d£)/(dzx / d£)) = иш , найдем значение скорости продвижения контура vc = иш / е и интеграл соотношения (10) для контура в параметрическом виде
X = x — vkt = £1~а cosжа/ vc, У = У = £1~а sin жа/ vc, £ е (0, да).
Таким образом, контур описывается прямой линией у = а(ех — ujt ), а обтекаемый массив есть клин с углом раствора 2жа = 2arctg(еa).
4. Абляция солевого диапира под основанием ГЭС. Пусть имеется солевой диапир, оголовок которого находится на глубине у °° м под руслом реки. Стационарное положение оголовка диапира можно объяснить естественным ростом последнего с некоторой скоростью Vd м/год. Пусть 21 м -ширина флютбета, перепад напора между верхним и нижним бьефом плотины АН м, естественный гидравлический уклон грунтовых вод под руслом реки i . Исходя из этих данных, можно определить динамику снижения оголовка во времени после возведения плотины. Предварительно рассмотрим режим равновесного с потоком грунтовых вод положения диапира. Пусть ось x направлена вдоль русла реки, а ось у - вертикально вниз. Чаще всего слоисто неоднородная толща пород представляет собой предельно анизотропную среду, а коэффициент фильтрации k( y) кусочно-постоянная функция. В пределах каждого слоя компоненты скорости фильтрации U, v определяются в соответствии с законом Дарси и = к(у)др/ дх, v = к(у)др/ ду и в отсутствии источников (стоков) удовлетворяют условию неразрывности
ди / дх + дv/ду = 0. (11)
Будем исходить из хорошо зарекомендовавшей себя в инженерной практике гидравлической схемы фильтрации воды в слоистом грунте [2]. Исследования земной толщи в районе ГЭС показывают, что проницаемость слоёв с глубиной убывает, следовательно, убывать будет с глубиной и горизонтальная проекция U скорости фильтрации, а вертикальную проекцию скорости v считаем пренебрежимо малой. Полагая в уравнении (11) второе слагаемое равным нулю и интегрируя оставшуюся часть дважды по X, получим выражение (р = C(у)х + <р0 . Здесь р0 - некоторая постоянная, а зависимость С(у) = и(у) может быть только линейной, поскольку в каждом слое грунта для вертикальной скорости фильтрации д2р/ду2 =д2С/ду2 = 0 . Учитывая линейное распределение гори-
зонтальной скорости по глубине в условиях естественного режима фильтрации, для ординаты диапира можно выписать дифференциальное уравнение вида
Фа / dt = аи01—(уа/ Н)]—Vd,
где и 0 и Н - значения скорости фильтрации у поверхности земли и глубина толщи, на которой скорость обращается в нуль, соответственно. В равновесном состоянии скорость абляции диапира равна скорости его роста и правая часть выписанного уравнения равна нулю. Обозначим это значение
у = у0 . Поскольку и0 = ki0, то для кинетического коэффициента получаем его естественную связь со скоростью роста диапира в виде соотношения
ак\> = Vd/(1—у^/ Н) .
После сооружения плотины с перепадом напора АН и длиной флютбета 2l гидравлический уклон i = АН /(2l) >> i0. Как следует из классической формулы для комплексной скорости фильтрации [2]
и — iv = к АН /(ж^1^—ХЛу)', в любой точке с координатами (х, у) под флютбетом плотины убывание скорости с глубиной у стремится к нулю по линейному закону. Таким образом, и в этом случае уравнение для ординаты оголовка диапира будет в некотором приближении справедливым дифференциальным уравнением вида
^у„ / dt = Vd [ft/ 0(1—ул / Н )/(1 — у0/ Н) — 1]. (13)
Начальным данным для этого уравнения будет условие уа (0) = у0. Функция уй экспоненциально убывает с течением времени
Уd (J) = уда — (уда— у0) exp(—т),
где т = Vdt / Н - безразмерное время, а у да - стационарное состояние диапира
уда = н—(ij аН—у0).
В таблице указаны моменты времени растворения диапира до глубины 100 м под флютбетом плотины при различных Н и Vd при ширине флютбета 1560 м, перепаде напора плотины 300 м ,
уклоне реки 0.002.
Таблица
Vj, м/год 51G-2 2.51G-2
H, м 200 300 200 300
t, год 22.S 19.6 11.4 9.S
5. Обтекание округлых массивов. Пусть симметричный относительно осей координат массив имеет первоначально округлую форму, очертания которой в каждый момент времени можно аппроксимировать луночками, ограниченными дугами окружности, радиус которых изменяется с тече-
72S
нием времени в соответствии с протекающими процессами размыва, разложения или растворения. Область движения представляет верхняя полуплоскость с выброшенной круговой луночкой, образующей изменяющийся с течением времени угол жас (t) с осью X. Радиус луночки
R = sin-1(жас (t)) , а её возвышение ус (t) = tg(жас /2) . В точке контура A(—1,0) раздвоения потока флюида и точке схода -6(1,0) скорость обращается в нуль, поэтому эти точки в процессе движения будут неподвижными. Если в начальный момент времени сечение массива представляло собой полный круг, то ас (0) = 1/2, Яс (0) = ус (0) = 1. Конформное отображение области движения в плоскости z на верхнюю полуплоскость вспомогательной переменной £ с соответствием точек z = ±1 <^£ = ±1,z ^ да ^ ^ да осуществляет функция
£ = g(z, t)=С1+N1(z, t))/(1—N1(z, t)),
где N = )(z), D(z) = (z —1)(z +1) . Функция, представляющая обратное отображение, имеет
вид z = z(£;t) = (1 + N2)/(1 — N2), где N2 = D(£).
Отображение области потенциала f (z; t) на верхнюю полуплоскость вспомогательной переменной £ с соответствием точек границы, указанных на рисунке, осуществляет функция w(£) по формуле Кристоффеля-Шварца [1]
£
w(£;t) = M(t) JDа (£)d£ .
Здесь множитель М ^) определяется условием задачи (5) на бесконечности
= 1т( dw/Л£)/(/ Л£), £ ^ да .
Подставив сюда значения функций и выполнив разложение в бесконечно удаленной точке, в пределе получим М^) = 2мш /(1 — ас ^)). Изменяющееся с течением времени возвышение ус ^) круговой
луночки (или угол ас (t)) можно определить из условия (2), записанного для одной точки = /ус, которая отвечает точке вспомогательной плоскости £ = 0. Учитывая симметрию течения относительно оси у , получим уравнение
dyc ж dаc df .. .
v = —- =--------- ------------ = —а Re—(iyc) = —а Re
dt 2 cos (жас /2) dt dz
Г— (0)/ А (0)Л d£ d£
После подстановки значений функций w(£; t), г(£; t) и несложных выкладок получим обыкновенное дифференциальное уравнение
ж — /3 dp /1-1, 2 2
= —жяигг, /л/1 + а е
(1 + cos 3) dt
Решение этого уравнения при начальном условии 3(0) = ж / 2 имеет вид
/>А + а2е2 . (14)
1
Р
I (Р) = J (ж — x)(1 + cos x)—2 dx = 1/6(—1 — 4ln(cos(p/2)) + cos ~2(p/2) + (ж —p)(2 + cos2(p/2)tg (p/2)
0
В общем случае задача (5) может быть решена численно.
Рис. Массив в виде квадрата. Расчеты методом конечных элементов.
На рисунке приведены результаты расчетов на различные моменты условного времени, полученных методом конечных элементов, для массива, первоначальной формой которого был квадрат, описывающий круг единичного радиуса.
Поступило 13.07.2011 г.
ЛИТЕРАТУРА
1. Лаврентьев В.А., Шабат Б.В. Методы теории функций комплексного переменного. - М.: Наука, 1973, 680 с.
2. Полубаринова-Кочина П.Я. Теория движения грунтовых вод. - М.: Наука, 1977, 664 с.
Н.К.Корсакова, В.И.Пенковский, М.А.Саттаров*
МОДЕЛИ МАТЕМАТИКИИ ОБШУШТАИ ^ИНС^ОИ ХОС АЗ ДОХИЛИ ЦИШРИ МАНША
Институти гидродинамикаи Шуъбаи Сибирии Академияи илмхои Россия, *Институти масъала^ои об, гидроэнергетика ва экологияи Академияи илмхои Цум^урии Тоцикистон
Модели математикии чараёни тагйирёбии массаи аз сатхд обногузаронандаи намаксанг (ба мисли диапирх,о), дар дохили кишри манша ва ё айсбергх,о, хднгоми шустани онх,о бо об, пешних,од шудааст. Х,алли х,олатх,ои дученакаи обшуй пешних,од шуда, хдлли амики аналитикии онх,о мавриди назар карор гирифтааст. Дар долати умумй усули элементх,ои чузъй истифода шудааст. Тач^изоти тачрибавй сохта шуда, дар он тачрибах,о гузаронида шуданд, ки натичааш бо натича^ои х,исобй сифатан хдмохднганд.
Калима^ои калиди: намаксанг - диапир - яхбастагй - полоиш - шорогузоиш - модел - озмоиш -уисобкунй - сарбанди об - гидрогиреу.
H.K.Korsakova, V.I.Penbkovskii, M.A.Sattarov*
MATHEMATICAL MODEL OF PROCESSES OF THE DIAPIR ABLATION
Institute of Hydrodynamics, Siberian Branch of the Russian Academy of Science *Institute of Water problems, Hydroenergetics and Ecologic,
Academy of Sciences of the Republic of Tajikistan A mathematical model of mass loss from surfaces of impermeable massifs of rock salt (diapir, for example), closed frozen earth inclusions or icebergs in potential flow are proposed. The two-dimensional potential flow cases, for which one can write some analytical solutions, are considered. Generally, the finite-element method is applied. The experimental arrangement has been constructed. The preliminary experiments have conducted. The results of the preliminary experiments are agreed qualitatively with the numerical results.
Keywords: rock-salt -diapir -frozen ground -filtration -flow round - model -experiment -calculation -dam
- hydro-electric power station.
