Математическое моделирование течений неоднородной несжимаемой жидкости
На правах рукописи
Худобина Юлия Петровна МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ НЕОДНОРОДНОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ 01.02.05 – Механика жидкости, газа и плазмы
Автореферат диссертации на соискание ученой степени кандидата физико-математических наук
Томск - 2009
Работа выполнена на кафедре динамики полета ГОУ ВПО «Томский государст венный университет» и НИИ прикладной математики и механики ТГУ
Научный консультант: кандидат физико-математических наук, ведущий научный сотрудник НИИПММ при ТГУ Либин Эдуард Ефимович
Официальные оппоненты: доктор физико-математических наук, профессор Бубенчиков Алексей Михайлович доктор физико-математических наук, профессор Троицкий Олег Юрьевич
Ведущая организация: Институт гидродинамики им. М.А. Лаврентьева СО РАН (г. Новосибирск)
Защита состоится 25 декабря 2009 г. в 14 час. 30 мин. на заседании диссер тационного совета Д 212.267.13 при ГОУ ВПО «Томский государственный университет» по адресу: 634050, г. Томск, пр. Ленина, 36, НИИ прикладной ма тематики и механики, ауд. 239.
С диссертацией можно ознакомиться в Научной библиотеке ГОУ ВПО «Томский государственный университет» по адресу: г. Томск, пр. Ленина, 34а.
Автореферат разослан 20 ноября 2009 г.
Ученый секретарь диссертационного совета Д 212.267. доктор технических наук Ю.Ф. Христенко
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Актуальность темы: Практическое приложение гидродинамика неодно родной несжимаемой жидкости имеет ко многим геофизическим процессам, протекающим в атмосфере и в океане, и поэтому она интересна для специали стов по геофизике, метеорологии и океанологии. Однако в теоретическом плане гидродинамика неоднородной жидкости изучена еще явно недостаточно. При чиной этого является то, что уравнение неразрывности и условие несжимаемо сти для неоднородной жидкости не совпадают друг с другом и записываются как независимые уравнения. В результате получается, что гидродинамика неод нородной жидкости имеет дело с принципиально вихревыми течениями, кото рые с большим трудом поддаются математическому анализу. Неоднородные жидкости никогда не находятся в состоянии покоя, так как малейшее наруше ние локальной плотности приводит к тому, что отдельные участки жидкости всплывают под действием архимедовой силы, а другие опускаются. Само по себе наличие силы тяжести и неоднородного поля плотности является доста точной причиной для поддержания различного рода специфических движений жидкости. Особый интерес вызывают такие явления как внутренние волны, всплывание и опускание термиков, конвекция, перемешивание в зоне скачков плотности, и др. Степень изученности этих проблем различна. Так теория внутренних волн интенсивно разрабатывалась в последние десятилетия, и в ней получены важные результаты. В работе авторов: А.В. Аксенов, В.А. Городцов, И.В. Стурова «Моделирование обтекания цилиндра стратифицированной иде альной несжимаемой жидкости», (1986, препринт № 282) было получено реше ние задачи в аналитической форме. В подобных задачах, связанных с решением нелинейных дифференциальных уравнений в частных производных, до сих пор остаются значительные затруднения.
Таким образом, своевременными и актуальными являются исследования по построению математических моделей, позволяющих численными методами изучать свойства течений неоднородных жидкостей.
Целью диссертационной работы является:
1. Теоретическое изучение стационарных течений неоднородной жидко сти, таких как обтекание устойчиво стратифицированной жидкостью диполя, источника, и других особенностей. Получение аналитического решения в част ных случаях.
2. Разработка новой математической модели, позволяющей достаточно простыми средствами осуществлять расчет различных нестационарных плаву чих эффектов, проявляющихся в неоднородной жидкости, таких как: всплыва ние плоских и осесимметричных термиков, явление тейлоровской неустойчиво сти, распространение волн вдоль границы раздела двух жидкостей с разной плотностью. Предлагаемая математическая модель позволяет приближенно ис следовать также и поведение волн на свободной поверхности тяжелой жидко сти.
3. Распространение полученных результатов на движение неоднородных жидкостей, подчиняющихся закону Дарси в пористых средах.
Научная новизна работы 1. Впервые получено решение в аналитической форме для стационарного течения неоднородной несжимаемой жидкости в частном случае, когда плот ность жидкости зависит от значений функции тока по показательному закону.
2. Впервые получена формула, представляющая собой закон генерации вихрей на линии, разделяющей несмешивающиеся жидкости с разной плотно стью.
3. На основе этого соотношения построена новая вихревая математическая модель, позволяющая вычислять движение границы раздела двух сред. Несмот ря на такое упрощение, состоящее в рассмотрении только двухкомпонентной жидкости, предлагаемая модель получается достаточно содержательной, чтобы с ее помощью можно было вычислять различные физические эффекты, рас смотренные в диссертации.
Практическая значимость. Создан комплекс компьютерных программ, позволяющих проводить вычисление движения границ разделяющих разнород ные жидкости. На примере решения задачи о нефтяной скважине показано, что разработанная математическая модель применима и к случаю нестационарного движения неоднородных грунтовых вод. При дальнейшем естественном обоб щении предлагаемого метода круг решаемых им задач может быть расширен.
Достоверность полученных результатов обеспечивается строгостью мате матических постановок и подтверждается сопоставлением с эксперименталь ными данными и численными исследованиями других авторов.
Основные положения, выносимые на защиту:
1. Аналитическое решение задачи об установившемся течении неоднород ной несжимаемой жидкости.
2. Новая математическая вихревая модель, позволяющая вычислять дви жение границы между двумя жидкостями с различной плотностью.
3. Численный метод, в основе которого лежит простой закон генерации вихрей, возникающих вдоль границы раздела.
4. Распространение полученных результатов на случай движения неодно родных жидкостей в пористой среде.
Апробация работы. Основные положения и результаты диссертации док ладывались и обсуждались на XI Всероссийской научно-технической конфе ренции под редакцией Шрагера Э. Р. (Томск, 2005). На V Всероссийской кон ференции “Фундаментальные и прикладные проблемы современной механики” (Томск, 2006). Н III Всероссийской конференции молодых ученых “Физика и химия высокоэнергетических систем” (Томск, 2007). На VI Всероссийской конференции “Фундаментальные и прикладные проблемы современной меха ники”, посвященной 130-летию Томского Государственного университета и 40 летию НИИ ПММ ТГУ (Томск, 2008).
Публикации. Основные результаты диссертации представлены в трудах вышеперечисленных конференций, а также в журнале из списка ВАК “Известия вузов. Физика”.
Объем и структура работы. Диссертация состоит из введения, 2-х глав, заключения и списка цитированной литературы. Объем диссертации – страницы, в том числе 124 страницы основного текста, 58 рисунков, 6 фотогра фий. Список литературы – 9 страниц и содержит 100 наименований.
КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ
Во введении сформулированы цель и задачи диссертационного исследо вания, обоснована актуальность, сформулированы выносимые на защиту ос новные положения, кратко излагается содержание диссертационной работы по главам.
Первая глава посвящена обзору состояния вопроса и решению задачи об установившемся течении неоднородной несжимаемой и невязкой жидкости.
Как известно, в двумерном (плоско-параллельном) случае установившееся (не зависящее от времени) течение неоднородной жидкости описывается системой нелинейных дифференциальных уравнений в частных производных:
u u 1 p +v = u x x y v v 1 p u +v = g y x y, (1) u v + = x y +v = u x y где: u, v –горизонтальная и вертикальная скорости: p – давление;
– плот ность жидкости. Если ввести функцию тока ( x, y ), такую, что скорости через нее выражаются по формулам u = y, v = x, то условие несжимаемости жидкости ux + v y = 0 будет выполняться автоматически. Уравнение сплошности u x + v y = 0, выраженное через функцию тока, означает равенство нулю функ ционального определителя D(, ) D( x, y ) = 0. Отсюда следует, что = ( ), и линии тока в установившемся движении одновременно являются и линиями по стоянной плотности. Подстановка этих результатов в первое и второе уравне ния из системы (1), и исключение из них давления, приводит к дифференциаль ному уравнению, которому должна удовлетворять функция тока ( x, y ) :
1 d 1 1 = F ( ). (2) + gy + + d 2 x 2 y Уравнение (2) впервые было получено в работах Дюбреиль-Жакотен (1932) в связи с ее исследованиями о волнах на поверхности неоднородной тяжелой жидкости. Оно содержит две произвольные функции: ( ) и F ( ). Первая из них задает распределение плотности в жидкости, а вторая – распределение вих рей в ней. В безразмерных переменных, когда выбирается масштаб R0 по про странственным переменным и масштаб скорости потока U 0, уравнение (2) со храняет свой вид. Только вместо ускорения силы тяжести g в него будет вхо дить квадрат безразмерного числа Фруда, определяемого как: Fr = gR0 U 0.
В диссертации показывается, что вместо функции тока ( x, y ) удобнее на ходить вспомогательную функцию Q( x, y ), которая связана с функцией соот ношением ( ) Q ( x, y ) = d. (3) Тогда, из (2) и (3) следует, что функция Q должна удовлетворять диффе ренциальному уравнению:
d F ( ). (4) Q + 2 Fr 2 y = 0 d Уравнение (4) является более простым по сравнения с уравнением (2), так как оно не содержит квадратов первых производных от искомого решения. В него нужно подставить зависимость = (Q), найденную из уравнения (3). В диссертации особое внимание уделяется случаю, когда F ( ) = 0 и ( ) = 0 exp(2 ), т. е. когда плотность жидкости распределена по экспонен циальному закону с коэффициентом неоднородности имеющим значение. В этом случае функции Q и связаны между собой соотношениями:
1 e 1, (5) Q= =, ln 1 Q и уравнение (4) становится линейным, т. е.:
2Q 2Q + 2 Fr 2 y Q = Fr 2 y. (6) + x 2 y Частное решение неоднородного уравнения (6) можно взять здесь не зави сящим от переменной x, и оно тогда выражается через функции Эйри. Решение однородного уравнения (6) можно представить в виде, содержащим произ вольную (порождающую) функцию Q0 ( x, y ). Предположим, что если = 0, то решение Q0 ( x, y ), удовлетворяет уравнению Лапласа, и описывает некоторое те чение однородной жидкости. Тогда, применяя преобразование Лапласа (по на правлению x) к уравнению (6), получим операционное соотношение:
e Q( x, y)dx = e s 2 + 2 Fr 2 y x sx Q0 ( x, y )dx 0 Таким образом, связь между функциями Q и Q0 будет такой же, как и связь между оригиналами при замене в изображениях переменной s на величину s 2 + 2 Fr 2 y. Оператор такого преобразования известен из теорем об обобщен ной свертке Эфроса, и в окончательном виде решение линейного уравнения (6) может быть записано в форме:
) ( Q ( x, y ) = Q0 (0, y ) + x J 0 b 1 2 Q0 ( x, y ) d + ) ).
( ( (7) 1 Bi (0) Ai y Fr Ai (0) Bi y Fr 3 2 2 2 + Bi (0) Ai ( 0 ) Ai (0) Bi ( 0 ) Здесь: b = x Fr y, Q0 ( x, y ) = Q0 x, и добавлено частное решение неод нородного уравнения (6), выражающееся через функции Эйри. Входными па раметрами при расчетах по формуле (7), кроме порождающей функции Q0 ( x, y ), является величина, характеризующая степень неоднородности жидкости, и число Фруда Fr. Вид решения (7) показывает, что функция Q( x, y ) зависит от параметра, и от величины произведения Fr, но соотношение (5), выра жающее искомую функцию тока ( x, y ) через Q, содержит только параметр.
Решение в форме (7) имеет достаточно универсальный вид, так как оно по зволяет проследить поведение линий тока ( ( x, y) = const ) при изменении пока зателя неоднородности. В предельном случае, когда 0, уравнения (7) и (5), очевидно, дают: ( x, y ) = Q0 ( x, y ). В диссертации приводится достаточно простая программа для ЭВМ, которая заполняет сеточную прямоугольную об ласть (x,y) значениями функции ( xn, ym ), вычисленными по формулам (7), (5), и рисует ее изолинии. Приводятся соответствующие примеры расчетов для раз личного вида порождающих функций Q0 ( x, y ).
1. Диполь в потоке неоднородной жидкости.
Если в равномерное течение со скоростью U 0 помещен диполь, то порож дающая функция (в безразмерных переменных), и ее частная производная по переменной x, будут иметь вид:
Q0 2 y ( x c).
Q0 ( x, y ) = y 1 =, ( x c) + y x ( x c)2 + y Характер линий тока для этого случая представлен на рисунке 1. Из него видно, что в потоке, в отличие от однородной жидкости, образуются нелиней ные волны вниз по течению.
Рисунок с похожим поведением линий тока, вычисленный с применением разностных методов, встречается в обзорной монографии Дж. Тернера “Эф фекты плавучести в жидкостях”(1977).
Рисунок 1 – Линии тока в неоднородной жидкости, обтекающей диполь:
Fr=3;
=0, 0.1, 0.25.
2. Кусочно-экспоненциальное распределение плотности.
Решение в виде уравнений (5), (7) будет справедливым и для случая, когда показатель неоднородности ( ) принимает постоянные, но различные значе ния на отдельных участках изменения величины. Пусть, например, в задаче об обтекании диполя, на промежутке = 0 ±, параметр принимает значе ние равное 0, и значение равное нулю за пределами этого промежутка. Тогда мы имеем течение двухслойной устойчиво стратифицированной жидкости, в котором слои разделены участком с шириной 2. На рисунке 2 представлен ре зультат расчета, когда 0 = 2, и = 0.1. Линия тока соответствующая значению 0 = 2 (граница раздела) обозначена здесь более темной линией.
Рисунок 2 – Граница раздела между однородными жидкостями со скачком плотности, обтекающими диполь.
Из этого расчета видно, что граница раздела является сложной линией, и сопровождается гидравлическими прыжками. Такую картину можно предста вить себе как начало процесса перемешивания. На удалении от границы раздела линии тока являются плавными, такими как при течении однородной жидкости.
3. Источник в экспоненциально стратифицированной жидкости.
Одним из простейших течений в гидродинамике несжимаемой жидкости является течение от источника, помещенного в начало координат. Функция то ка такого течения, и ее производная по переменной x имеют вид Q y y.
Q0 ( x, y ) = arctg, = x x + y x Линиями тока для такой порождающей функции являются прямые линии, ис ходящие из начала координат. Значения функции Q0 ( x, y ) изменяются в первом квадранте от нуля, на горизонтальной стороне угла, до / 2 на его вертикаль ной стороне. Если плотность неоднородной жидкости убывает по показатель ному закону: = 0 exp(2 ), то справедливо решение в виде формул (7) и (5).
Пример расчета по этим формулам приводится на рисунке 3.
Рисунок 3 – Источник помещен в начало координат. Влияние увеличения числа Фруда на изменение характера линий тока:
Fr 2 = 0.001, 0.1, 0.25, 1, 2;
= 0.05.
Из этого примера видно, что под влиянием силы тяжести и градиента плот ности линии тока (траектории частиц жидкости) имеют как восходящие, так и падающие участки. Течение может сопровождаться гидравлическими скачками в виде ‘языков’ на большом удалении от источника. Число Фруда определяется в данном примере по прежней формуле Fr = gR0 U 0, но под величиной U 0, ха рактеризующей мощность источника, подразумевается радиальная скорость на круге с радиусом равным R0. На рисунке 3 показано, как влияет величина числа Фруда на характер линий тока. При малом числе Фруда, когда Fr 0, источник преодолевает влияние силы тяжести, и картина течения незначительно отлича ется от аналогичного течения однородной жидкости. Наоборот, если преобла дает сила тяжести, то появляется искривленное и расслоенное течение в виде гидравлических скачков.
Численное исследование, проведенное в первой главе диссертации, указы вает на существенное различие между течениями однородной и неоднородной жидкостями. По-видимому, экспоненциальное распределение плотности явля ется единственным частным случаем, когда нелинейное дифференциальное уравнение (2) для функции тока можно свести к линейному уравнению, и даже получить его решение в замкнутой форме.
Во второй главе изучаются отдельные задачи, относящиеся к неустано вившимся движениям неоднородной жидкости. В ней показывается, что на ли нии раздела двух или нескольких жидкостей с разной плотностью, находящих ся в поле силы тяжести, образуется вихревой слой, который свободно движется по своему кинематическому закону, и является возбудителем движения всей массы жидкости. Закон изменения во времени интенсивности этого вихревого слоя можно найти из гидродинамического уравнения движения, записанного в системе координат движущейся вместе с частицами жидкости:
p dV g. (8) = dt Умножая (8) на, и вычисляя циркуляцию скорости по замкнутому кон туру, охватывающему какую-либо точку линии раздела (рисунок 4) получим d dt ( V, ) ds = ( g, ) ds (9).
Здесь – вектор каса тельной (рисунок 4), а ds – линейный элемент дуги контура интегриро вания.
Равенство (9) запи Рисунок 4 – К определению соотношения шем с учетом разных для производной. Интегрирование по значений плотности контуру, охватывающему участок границы жидкости по обе сторо раздела.
ны границы раздела, т.е.
d d (V, ) d s + (V, ) d s = dt dt (10) 2 2 (g, ) d s (g, ) d s = 2 2 Для циркуляции вектора скорости V и для циркуляции постоянного векто ра g имеют место очевидные равенства (V, ) d s (V, ) d s = =, (11) 2 (g, ) d s (g, ) d s (g, L ).
= = 2 Вектор L здесь представляет собой хорду, проведенную между точками пересечения контура интегрирования и линии раздела. Равенство (10), следова тельно, можно представить в следующем виде (1 + 2 ) = ( 1 2 ) (L, g ). (12) Отсюда получаем, что d (g, L ). (13) =2 1 1 + dt При стягивании контура интегрирования в точку, мы получаем основное уравнение, определяющее генерацию вихрей на линии раздела плотностей d d (g ). (14) =2 1, + ds dt 1 Здесь вектор имеет смысл касательной к линии раздела плотностей, а ds – элемент длины дуги вдоль этой линии. Правило обхода, при интегрирова нии по контуру, применяется обычное, такое, чтобы внутренняя область оста валась с левой стороны. Полученная формула (14) означает, что полная произ водная по времени для погонной интенсивности вихрей, образующихся на гра нице раздела, пропорциональна скачку плотности и проекции ускорения силы тяжести на касательную к линии раздела.
В двумерной плоскопараллельной задаче линию вихревого слоя можно представлять последовательностью точечных дискретных вихрей с комплекс ными координатами Z n и соответствующими циркуляциями n, изменение ко торых во времени определяется из соотношения (14). Скорость отдельных вих ревых точек в слое тогда будет равна m (15) = * V.
2 i n Zn Z m m Алгоритм численного моделирования движения такой вихревой пелены, разделяющей несмешивающиеся жидкости с разной плотностью, следователь но, сводится к решению задачи Коши для системы обыкновенных дифференци альных уравнений. Если для решения применять схему Эйлера, то координаты точечных вихрей и их циркуляции находятся по рекуррентным формулам Z nk + 1 = Z nk + t V nk, k + 1 = n + t n.
k k (16) n Здесь величины n и V n вычисляются соответственно из соотношений вида (14) и (15). Практически, для вычисления скорости перемещения сразу всех точек вихревой линии по формуле (15), удобно организовать специальную матрицу Q, такую, чтобы выполнялось равенство: V = Q. Для вычисления векторов касательной и нормали к линии раздела плотностей можно создать матрицу Н, такую, чтобы касательная во всех точках контура определялась из соотношения: = H Z, где Z – вектор, состоящий из списка комплексных ко ординат точек контура. Матрица H обычно имеет простой вид и получается на основе метода центральных разностей. Для замкнутого контура она должна быть циклической.
В безразмерных переменных уравнения (14) и (15) сохраняют свой вид.
При этом за единицу расстояния R принимается характерный размер области.
Безразмерными единицами для скорости и циркуляции тогда будут величины V 0 = gR, 0 = R gR ;
и безразмерное время имеет вид = g R t.
Типичными задачами, которые можно решать с помощью описанной выше вихревой модели движения границы раздела, являются задачи о всплывании или опускании термиков, а также задача о прорыве легкой жидкости из-под тя желой (тейлоровская неустойчивость). Примеры такого расчета, полученные в диссертации, показаны на рисунках 5, 6.
Рисунок 5 – Влияние плавучести на подъем более легкого жидкого тела (а), и на опускание более тяжелого тела (b).
На рисунке 5 (а), показаны кадры всплывания легкого жидкого круга еди ничного радиуса, которые были вычисленные через равные промежутки без размерного времени = t 2 g / r. Расчет проводился с шагом = 0.25, и отно шение плотностей было равным 2 1 = 5. Аналогичный расчет опускания бо лее тяжелого жидкого тела (с обратным отношением плотностей) показан на рисунке 5 (b). Здесь учитывалось влияние дна путем добавления зеркально от раженной вихревой системы.
Если рассматривать незамкнутую границу раздела двух жидкостей, то ее поведение, при аналогичном численном моделировании, проясняет характер неустойчивости по Тейлору. На рисунке 6 показан пример такого расчета.
Рисунок 6 – Изменение границы раздела, когда слой более легкой жидкости в начальном положении находится ниже тяжелой, здесь отноше ние плотностей s = 5.
Начальное положение границы раздела здесь было почти горизонтальным и представлялось уравнением:
y = 1 + 0.1 e x. (17) Такое начальное положение выглядит на рисунке 6 как незначительное возвышение, разделяющее легкую и более тяжелую жидкости. Пунктирной ли нией на рисунке 6 показано промежуточное положение границы, когда безраз мерное время имеет численное значение равное 0.5. Само по себе явление не устойчивости по Тейлору означает, что положение системы, когда более легкая жидкость находится ниже тяжелой, является неустойчивым, и происходит са мопроизвольное нарастание отклонений от невозмущенного состояния. Приме няемая в диссертации вихревая модель показывает, что при прорыве легкой жидкости из-под тяжелой жидкости образуется «грибовидное облако», граница которого имеет характерную спиральную структуру. В случае, когда начальное положение границы между легкой и тяжелой жидкостью имеет несколько сла бых бугорков, все они всплывают в виде характерных отдельных «термиков», как показано на рисунке 7.
Рисунок 7 – Эффект неустойчивости по Тейлору в случае нескольких слабых возвышений на начальной границе раздела между легкой и тяже лой жидкостью. Отношение плотностей s = 1 2 = 5.
Подобные расчеты могут рассматриваться как моделирование возникнове ния естественных конвективных течений в задачах теплофизики.
Если жидкость устойчиво стратифицирована по плотности, когда более тяжелая жидкость находится ниже легкой, то расчет по вихревой модели нико гда не приводил к образованию спиральной закрутки линии раздела. Вдоль нее распространяется обычная нелинейная бегущая волна, сопровождаемая шлей фом из-за влияния дисперсии (рисунок 8).
Рисунок 8 – Распад начального возвышения в виде гауссовского импульса единичной высоты в устойчиво стратифицированной жидкости.
Здесь приведен один из кадров расчета движения границы раздела перво начально покоящейся устойчиво стратифицированной жидкости. Начальное возвышение здесь описывалось импульсом вида (17), но с единичной высотой и с более мелким слоем нижней, тяжелой жидкости. Собственно говоря, расчет, показанный на рисунке 8 – это просто повторение расчета, показанного на ри сунке 6, но при обратном отношении плотностей. В одном случае (s1) имеет место неустойчивость по Тейлору, а в другом, когда s1, по поверхности разде ла распространяется волна, похожая на волны в жидкости со свободной по верхностью.
Кроме вычисления положения границы раздела можно также графически отображать и мгновенные линии тока для отдельных моментов времени. На ри сунке 9 показаны изолинии функции тока для вихревой системы всплывающего легкого «термика».
При таком расчете вместо последовательности точеч ных вихрей применяется по следовательность прямоли нейных отрезков с непрерыв ным распределением вихрей на ней. Такие отрезки заме няют точечные вихри, и ре зультирующая функция тока суммируется от них обычным образом. Характер линий то ка имеет различный вид внутри и снаружи замкнутой области, такой же, как в элек тродинамических задачах с неоднородными средами.
Рисунок 9 – Поле течения (траектории час тиц жидкости) снаружи и внутри всплываю щего «термика».
Далее, в подразделе 2.6 второй главы рассматриваются осесимметричные задачи о «термиках» и тейлоровской неустойчивости плоской границы. Урав нение генерации вихрей на границе раздела остается таким же, как и для плос копараллельного случая (14), но компоненты скорости от вихревых колец, ко торыми покрывается граница между жидкостями с разной плотностью, вычис ляются с помощью полных эллиптических интегралов первого и второго рода:
r 2 ( z z ) r 1 E (k ), (18) u= K (k ) + 2 r1 r22 z z r 2 + ( z z )2 + r E (k ), (19) v= K (k ) 2 rr1 r2 4rr.
где: r12 = ( z z ) + ( r + r ) ;
r22 = ( z z ) + ( r r ) ;
2 2 2 k2 = r Центр вихревого кольца расположен на расстоянии z от начала коорди нат, а само кольцо имеет радиус r, циркуляция кольца равна. Другой отли чительной особенностью вихревых колец является то, что они имеют собствен ную скорость перемещения в осевом направлении. В покоящейся жидкости вихревое кольцо перемещается, не меняя своего радиуса, в отрицательном на правлении оси Oz (если циркуляция положительна) с собственной скоростью U, которая равна 8r. (20) U = ln 4 r Толщина вихревого шнура обычно выбирается из предположения, что ско рость жидкости в центре кольца совпадает с его собственной скоростью, что дает r 86. На рисунке 10 представлен пример расчета движения осесиммет ричных термиков.
Рисунок 10 – Фазы подъема (слева) и опускания (справа) осесимметричного термика.
Первоначальная граница раздела представляла собой покоящийся в жидко сти шар единичного радиуса, слегка приподнятого над горизонтальным экра ном. Графики положения границы раздела получены через одинаковые проме жутки безразмерного времени = 0, 1, 2, 3,..... Начальная интенсивность вихрей на начальной единичной окружности полагалась равной нулю, то есть движе ние начиналось из состояния покоя.
Задача о развитии тейлоровской неустойчивости в осесимметричном слу чае рассматривается в диссертации на примере, когда начальная граница между слоем легкой и тяжелой жидкостями слабо искривлена и описывается уравне нием: z = 1 + 0.1 e r. Начало движения такой границы показано на рисунке 11.
Рисунок 11 – Начальная фаза прорыва легкой жидкости сквозь тяжелую жидкость. = 0, 1, 2, 3, 4.
Сначала деформация границы раздела происходит медленно, но при продолже нии счета (рисунок 12) образуется последовательность всплывающих вихревых сгустков.
Рисунок 12 – Распределение вихревых центров на границе раздела в мо мент безразмерного времени = 7.
Последний пример показывает, что по сравнению с аналогичным плоско параллельным расчетом (рисунок 6) здесь проявляется динамическое взаимо действие вихревых колец, в котором особую роль играет собственная осевая скорость их перемещения. Вихревые центры распределяются вдоль границы разделяющей жидкости неравномерно, что характерно и для плоскопараллель ных задач.
Подраздел 2.7 второй главы диссертации посвящен задаче о поведении волн вблизи наклонного берега. Ее аналитическое решение известно только в приближении линейной теории волн малой амплитуды. Вблизи берега проис ходит разрушение волн, сопровождаемое высокими всплесками. Если свобод ную поверхность рассматривать как линию, разделяющую жидкости с разной плотностью, то задачу о накате волн на берег можно численно исследовать с помощью предлагаемой в диссертации вихревой математической модели. В ней большие всплески воды вдоль линии прибоя объясняются влиянием зеркально отраженной от линии берега вихревой системы. Для точки находящейся вблизи берега зеркально отраженная от него точка находится на малом расстоянии. Та кая пара близких вращающихся в разные стороны вихрей может только сколь зить вдоль линии наклонного берега, причем с довольно высокой скоростью.
Это и является причиной высокого всплеска и разрушения волн вблизи берега.
На рисунке 13 показан характерный вид мгновенных линий тока вместе с их аналитическим продолжением в пространство ниже береговой линии. Расчет здесь проводился с применением вихревой модели движения границы раздела.
Рисунок 13 – Вид мгновенных линий тока (траекторий частиц жидкости) при набегании волн на берег с углом наклона = 15.
В подразделе 2.8 второй главы изучается более сложная задача о подъеме или опускании термиков, которые находятся в неоднородном потоке жидкости, обтекающем эллиптическое препятствие. Комплексный потенциал течения во круг эллипса с полуосями a и b, и скоростью потока на бесконечности U, равен:
az b z 2 (a 2 b 2 ). (21) W0 = U a b В отличие от ранее решаемых задач здесь к ускорению силы тяжести до бавляется конвективное ускорение от заданного неоднородного потока. Для плоскопараллельной задачи его можно выразить через комплексный потенциал W0, в виде * dW0 d 2W A = ( V0, ) V0 =. (22) dz dz В безразмерных переменных в качестве единицы расстояния выбирается характерный размер термика R0, единицей скорости является величина U, а единица циркуляции вихрей равна 0 = R0 gR0. Тогда безразмерное время бу дет определяться как = t g R0. Система уравнений, определяющая движение термика в неоднородном потоке, записывается в виде:
d dz ( ) = Fr V0 + Q(z ). (23) =2 2 i + Fr 2 A, ds;
d 2 + 1 d Число Фруда (Fr) здесь является отношением Fr = U gR0, и единичный вектор i направлен вертикально вверх. Первое из уравнений (23) представляет собой закон генерации вихрей на границе термика, а второе описывает их кинемати ческое движение. Примеры расчета по формулам (23) показаны на рисунках и 15. Течение несущего потока, обтекающего эллипс, здесь направлено слева направо.
Рисунок 14 – Последовательность положений границы легкого термика в потоке, обтекающем эллипс. Fr=1.5, = 0,1,…, 7.
Рисунок 15 – Последовательность положений границы термика, который тяжелее окружающей жидкости. Fr=1.5, = 0,1,…, 7.
На рисунках 16 и 17 показано как изменяются мгновенные линии тока в присутствии всплывающего и опускающегося термиков.
Рисунок 16 – Характер поведения мгновенных линий тока вокруг всплы вающего термика.
Рисунок 17 – Влияние опускающегося термика на ход линий тока.
На этих рисунках видно характерное преломление линий тока при пересе чении границы термиков. Внутри термиков линии тока почти параллельны, что свидетельствует об их поступательном перемещении.
В последнем параграфе второй главы диссертации (подраздел 2.9) числен ным методом исследуется вопрос о применимости вихревой математической модели движения границы, разделяющей жидкости с разной плотностью, к за дачам фильтрации в пористой среде. Если в грунтовой воде имеется включение с другой плотностью, то оно должно всплывать или опускаться под влиянием силы тяжести. Поэтому можно ставить задачу о термиках в грунтовой воде и о других плавучих эффектах. Фильтрационное течение подчиняется закону Дар си, который означает, что скорость фильтрации пропорциональна градиенту давления:
p i. (24) V = k g Здесь коэффициент k характеризует фильтрационную способность данного грунта, и имеет размерность скорости. Единичный вектор i направлен по вер тикали вверх и учитывает влияние силы тяжести. К уравнению (24) добавляется еще условие несжимаемости divV = 0. Тогда существует потенциал скорости, который удовлетворяет уравнению Лапласа. На границе, разделяющей фильт рующиеся жидкости, терпит разрыв касательная составляющая скорости фильтрации. Такую границу, как и в обычной несжимаемой жидкости, можно заменить подвижным вихревым слоем, под действием которого и происходит движение всей массы жидкости. Умножая (24) на плотность, и интегрируя вдоль контура, охватывающего линию раздела, найдем, что условие генерации вихрей имеет вид:
d ( i, ). (25) = 2k 1 + ds Основное отличие этого условия от уравнения (14), полученного для обычной жидкости, заключается в том, что в него входит сама погонная интен сивность вихрей вместо ее полной производной по времени. Этим обстоятель ством определяется специфика поведения термиков в грунтовой воде. Скорость продвижения точек границы термика, как и раньше, определяется кинематиче ским операторным уравнением dz = Q. (26) dt Система обыкновенных дифференциальных уравнений (25) и (26) описы вает движение границы термика в пористой среде. Для ее приведения к безраз мерному виду выберем в качестве единицы расстояния характерный размер термика R0, а за единицу скорости примем величину k. Тогда масштабом цир куляции скорости будет служить величина 0 = R0 k, и безразмерное время опре деляется как R kt =. (27) t=, R0 k В таких безразмерных переменных система уравнений (25) и (26) записы вается в виде 1 dz ( i, ) ds. (28) = Q, = d 1 + Из нее видно, что уравнения задачи о движении термиков в пористой среде полностью нормализованы и не содержат безразмерных констант, характери зующих физическое явление. Влияние коэффициента фильтрации k проявляет ся только в выборе масштаба времени, согласно с формулами (27). Коэффици ент фильтрации k обычно бывает малым по величине. Так, например, для раз личного вида песков он изменяется в пределах (1 0.005 ) см с, а для глины его значения имеют в сотни раз меньшую величину. На приведенных ниже рисун ках представлены результаты расчетов по формулам (28). На рисунке 18 пока заны последовательные положения всплывающего и опускающегося термиков в пористой среде через промежутки времени, указанные на графике.
Рисунок 18 – Процессы подъема легкого и опускания тяжелого термиков в пористой среде. Отношение плотностей: 2 1 = 2, и 2 1 = 1 2.
На рисунке 19 изображен характер мгновенных линий тока, образующихся вокруг грунтовых термиков.
Рисунок 19 – Характер линий тока вблизи поднимающегося и опускающе гося термиков.
Поведение термиков в пористой среде существенно отличается от их пове дения в обычной жидкости. Здесь уже не образуется грибовидное облако, и граница термика не закручивается в спираль. При подъеме в грунте термик почти не меняет свою форму, в то время как опускающийся термик изменяет ее лишь из-за влияния горизонтального водоупора.
По-другому проявляет себя и тейлоровская неустойчивость в пористых средах. На рисунке 20 показаны последовательные положения линии раздела при прорыве более легкой воды через тяжелую. Незначительное возмущение горизонтального уровня довольно быстро приводит к образованию постоянно растущих языков и впадин.
Рисунок 20 – Проявление неустойчивости по Тейлору в пористой среде. Поло жение границы раздела показаны через промежутки безразмерного времени:
= 0, 4, 8, 12, 16, 18.
Если грунтовая жидкость является устойчиво стратифицированной, т.е. ко гда более тяжелая жидкость находится ниже чем легкая, то устойчивой будет и линия их раздела. На рисунке 21 показан результат расчета такой модельной задачи. Там, где в обычной жидкости распространяется волна (рисунок 8), в по ристой среде происходит не колебательный, а монотонный процесс растекания и выравнивания бугров над уровнем грунтовых вод. Задача о выравнивании де прессионной кривой, уравнение которой в начальный момент времени задано, относится к простым задачам об определении не известной свободной границы.
В книге П.Я. Полубариновой-Кочиной «Теория грунтовых вод»,1977 приводит ся аналитическое решение подобной задачи, полученное для случая бесконечно глубокого слоя грунтовой воды. Приведенный в этой книге график растекания бугра над уровнем грунтовых вод в основном совпадает с рисунком 21.
Рисунок 21 – Растекание бугров над уровнем грунтовых вод в случае ус тойчивой стратификации.
В этом же подразделе второй главы, посвященном фльтрации жидкостей в пористых средах, рассматривается также и известная задача о «нефтяной сква жине», которую ранее изучали П.П. Куфарев и П.Я. Полубаринова-Кочина. В классической постановке она выглядит как течение от источника (стока) с ну левым потенциалом на подвижной границе нефтяной залежи. С точки зрения развиваемой в диссертации вихревой модели задачу о «нефтяной скважине» следует представлять как задачу о движении двухслойной жидкости. Сила тя жести в этой задаче никакой роли не играет, так как она перпендикулярна плос кости движения. Но существует поле конвективных ускорений от источника, которое и заставляет завихриваться линию раздела между фильтрующимися жидкостями с разной плотностью. Если источник расположен в начале коорди нат, и его мощность равна J, то комплексный потенциал (для плоской задачи) имеет вид J ln( z ). (29) W ( z) = Скорость и конвективное ускорение от источника, очевидно, будут равны:
* * dW d 2W J2 z dW Jz = 2 4. (30) z = = A= ;
dz 2 z 4 z dz dz Так как конвективное ускорение A пропорционально квадрату мощности источника, то оно будет одинаковым как для источника, так и для стока, и на правлено всегда к центру источника.
Если считать, что расстояния z являются безразмерными величинами, от несенными к характерному размеру нефтяной залежи R0, то уравнение движе ния (закон Дарси) будет иметь вид p z + 4, (31) V = k g z где k — коэффициент фильтрации, и g = J 4 2 R03. Из уравнения (31) можно получить, чему равна циркуляция вихрей в точках границы раздела плотностей:
1 z d. (32) = 2kR0 2, 2 + 1 z ds Скорость течения жидкости равна сумме скоростей от источника и от вих ревой системы линии раздела. Следовательно, кинематическое уравнение за пишется в виде dz J z Q. (33) = + dt 2 R0 z 2 R Уравнения (32) и (33) можно нормализовать таким образом, чтобы они не содержали величин J и k. Для этого за единицу циркуляции скорости нужно принять величину kR0, а за единицу расстояния – величину R0 = J k. Тогда нор мализованные уравнения задачи примут следующий вид:
2 z d dz z = ± 2 + Q. (34) = 2 1,, 1 + 2 z 4 d ds z Знак «+» здесь выбирается для источника, а знак «–» выбирается для стока.
Безразмерное время связано с реальным временем t соотношением:
2 J k =. (35) t= t, 2 J k Уравнения (34) позволяют исследовать различные режимы движения неф тяной залежи. Кроме режима источника или стока, они позволяют рассматри вать принципиально различные случаи, когда вблизи источника находится лег кая либо тяжелая жидкость. В рамках традиционной постановки задачи о «неф тяной скважине» последние такие случаи разделить невозможно.
Расчеты эволюции контура залежи по нормализованным уравнениям (34) показаны на рисунках 22, 23. Результат зависит от смещения центра скважины относительно начальной залежи и от отношения плотностей. На рисунке 22 по казаны положения во времени линии раздела для случая, когда к скважине примыкает более легкая жидкость (нефть). Вода прорывается к скважине гораз до раньше, чем выкачивается нефть. Это вполне соответствует решению клас сической «задачи о нефтяной скважине». Рисунок 23 отображает процесс выка чивания воды из-под нефти. Вихри на границе раздела здесь вращаются в дру гую сторону и граница эволюционирует иначе: вода из-под нефти выкачивается полностью.
Рисунок 22 – Эволюция линии раз- Рисунок 23 – Положение линий дела для отношения плотностей раздела при большем смещении 1 2 = 0.8 скважины относительно центра за лежи. Отношение плотностей 2 1 = 0.8 ;
= 0.05.
В заключении, в диссертации приводятся основные выводы по работе, к ним относятся:
1. Установившееся течение в неоднородной жидкости возможно только в случае, когда плотность зависит лишь от функции тока. Если, например, плот ность зависит от декартовых координат (x,y), то такое течение может быть только неустановившимся. Это объясняет тот факт, что неоднородная атмо сфера не может находиться в состоянии покоя.
2. В отличие от течения однородной жидкости, линии тока в неоднородной жидкости получаются волнистыми, и внутренние волны распространяются на большие расстояния.
3. Моделирование неустановившихся течений путем замены границы раз дела плотности вихревым слоем позволяет решать многие задачи, недоступ ные при использовании иных способов. В частности в задачах, связанных с тейлоровской неустойчивостью, естественным образом моделируется процесс закручивания границы раздела в спиральную линию.
4. Тот факт, что устойчиво стратифицированная жидкость в численных экс периментах не закручивается в спираль, говорит о том, что вихревая модель правильно описывает свойства свободной поверхности. Следовательно, эту математическую модель можно применять и для изучения волн на поверхно сти воды.
5. Предлагаемая в диссертации вихревая модель оказывается пригодной и для течений жидкости в пористой среде, для которой справедлив закон Дарси.
Список опубликованных работ по теме диссертации 1. Худобина Ю.П. Линии тока в установившихся течениях неодно родной тяжелой несжимаемой жидкости // Известия вузов. Физика. – 2007.
– № 9/2. – С. 294–298.
2. Либин Э.Е., Худобина Ю.П. Эволюция границ раздела плотности в неоднородной несжимаемой жидкости // Известия вузов. Физика. – 2007. – № 9/2. – С. 291–293.
3. Акимова Ю.П. Задача о волнах, набегающих на наклонный берег // Ис следования по баллистике и смежным вопросам механики : сб. статей / под ред.
И.Б. Богоряда. – Томск : Изд-во Том. ун-та, 1999. – Вып. 3 – С. 3.
4. Худобина Ю.П. Кинематика волн у наклонного берега // Физика и хи мия высокоэнергетических систем : доклады XI всерос. научно-технич. конфе ренции / под ред. Э.Р. Шрагера. – Томск : Изд-во Том. ун-та, 2005. – С. 67–68.
5. Худобина Ю.П. Диполь в потоке неоднородной несжимаемой жидкости // Фундаментальные и прикладные проблемы современной механики : материа лы конференции, 3–5 окт. 2006 г. – Томск : Изд-во Том. ун-та, 2006. – С. 181 182.
6. Худобина Ю.П. Всплывание легкого жидкого тела // Фундаментальные и прикладные проблемы современной механики : материалы конференции, 3– окт. 2006 г. – Томск : Изд-во Том. ун-та, 2006. – С. 183.
7. Худобина Ю.П. Моделирование плавучих эффектов в жидкости // Фи зика и химия высокоэнергетических систем : сб. материалов III Всерос. конфе ренции молодых ученых, 24–27 апреля 2007 г., Томск. – Томск : ТМЛ-Пресс, 2007. – С. 232–234.
8. Худобина Ю.П., Либин Э.Е. Плоские гидродинамические течения, по рождаемые вихревыми слоями // Фундаментальные и прикладные проблемы современной механики : сборник материалов конференции, 30 сент. – 2 окт.
2008 г. / Том. гос. ун-т. – Томск, 2008. – С. 502–503.
9. Худобина Ю.П. Движение осесимметричного термика // Фундамен тальные и прикладные проблемы современной механики : : сборник материалов конференции, 30 сент. – 2 окт. 2008 г. / Том. гос. ун-т. – Томск, 2008. – С. 500– 501.
Тираж 100 экз.
Отпечатано в ООО «Позитив-НБ» 634050 г. Томск, пр. Ленина 34а