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

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2010. Вып. 1

УДК 532.546

Э. Н. Береславский, Л. А. Александрова, Е. В. Пестерев

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ НЕКОТОРЫХ ФИЛЬТРАЦИОННЫХ ТЕЧЕНИЙ В ПОДЗЕМНОЙ ГИДРОМЕХАНИКЕ

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

Построение подземного контура плотины с участками постоянной скорости обтекания. Рассматривается течение под водонепроницаемым подземным контуром заглубленной плотины ABCCiBiAi (рис. 1). Пусть контур основания плотины AAi состоит из двух вертикальных отрезков AB и AiBi одинаковой длины di, среднего горизонтального отрезка CCi и примыкающих к ним дуг кривых BC и BiCi с постоянной величиной скорости их обтекания |u>| = vo. Снизу область течения ограничена криволинейным водоупором FEFi, на котором величина скорости фильтрации также постоянна |u>| = uo(0 < uo < vo). Предполагается, что границы верхнего и нижнего бьефов горизонтальны, грунт однороден и движение подчиняется закону Дарси с известным коэффициентом фильтрации ж = const. Действующий на сооружение напор H, скорость обтекания vo и фильтрационный расход Q считаются заданными.

Береславский Эдуард Наумович — доктор физико-математических наук, профессор, заведующий кафедрой прикладной математики Санкт-Петербургского государственного университета гражданской авиации. Количество опубликованных работ: более 250. Научные направления: конформные отображения, аналитическая теория линейных дифференциальных уравнений класса Фукса, краевые задачи теории аналитических функций, математическое моделирование задач гидро- и аэромеханики. E-mail: eduber@mail.ru.

Александрова Людмила Александровна — аспирант кафедры прикладной математики Санкт-Петербургского государственного университета гражданской авиации. Научный руководитель: проф. Э. Н. Береславский. Количество опубликованных работ: 10. Научные направления: математическое моделирование, разработка прикладных программ. E-mail: mymila@mail.ru.

Пестерев Егор Васильевич — студент 4-го курса кафедры прикладной математики Санкт-Петербургского государственного университета гражданской авиации. Количество опубликованных работ: 5. Научные направления: математическое моделирование, разработка прикладных программ. E-mail: yogurt@live.ru.

© Э. Н. Береславский, Л. А. Александрова, Е. В. Пестерев, 2010

Введем комплексный потенциал движения ш = + гф, где - потенциал скорости,

ф - функция тока (область изменения переменной ш представлена на рис. 2) и комплексную координату г = х + гу, отнесенные соответственно к жЯ и Н, что позволяет при численных расчетах перейти к безразмерным величинам. Тогда V = grad^> = —жgradН, где напор Н = р+ у, р - давление, 7 - удельный вес фильтрующейся жидкости [1]. Задача состоит в определении положения кривых ВС, В\С\ и FE.F1 при краевых условиях

Ai.Fi : у = 0,<^> = —0.5Н, А1В1: х = —I, ф = Q,

С\С: у = —3, ф = Q, АВ: х = 1,ф = Q, (1)

А^: у = 0,^ = 0.5Н, ^Е^1: ф = 0

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

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

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

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

На рис. 3, а изображена область комплексной скорости V, соответствующая краевым условиям (1), которая ограничена дугами концентрических окружностей и отрезками прямых, проходящих через начало координат. Ввиду полной симметрии на плоскостях г, ш и V ограничимся рассмотрением области движения ЛБСЕЕЕ (см. рис. 1) и соответствующих ей одноименных областей на плоскостях ш и V (рис. 2 и 3, а).

Учитывая обилие прямых углов в плоскости V, удобно при конформном отображении в качестве канонической области плоскости г взять прямоугольник [4] (рис. 3, б) О < 11ет < 0.5, 0 < 1тт < 0.5/э, р(к) = К'/К, К' = К(к'), к' = %/1 — А’2,

а

б

Рис. 3. Область комплексной скорости ш (а) и вспомогательной параметрической переменной т (б)

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

т(т) = уо ехр(т — 0.5)пі. (2)

Конформно отобразим прямоугольник вспомогательной переменной т на область комплексного потенциала ш (см. рис. 2). В результате

И-г-ш________Мп(2Кт>А')______ ,,?1 ^

" К (к) Ь ^С8т кл/\ — Х2бп2(2Кт, к)’ 4 (3)

В данной формуле Е (р, т) - эллиптический интеграл первого рода при модуле т = АУ(1 - к'2 а2 [З2) / {1 - к'2а2), А = у/1 -к'2/32, а = вп{2Ка,к'), /3 = вп(2КЬ,к'), вп(р, А:),

сп(р, к) и dn(р, к) - эллиптические функции Якоби (соответственно синус, косинус

и дельта) при модуле к. При этом должно выполняться условие

К'Ы _ 2д

К (т) Я’ ( 1

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

Принимая во внимание соотношения (2) и (3), учитывая, что и> = вш/вг, и поступая аналогично тому, как это сделано в [5, 6], придем к зависимостям

(1ю Мвп(2Кт,к)сп(2Кт,к) (1г Мвп(2Кт, к)сп(2Кт, /г)ехр((0.5 — т)п{)

вт Д(т) 1 вт ^оД(т) ’

(5)

Д(т) = у/[1 — Л28П2(2Кт, к)][а2 + (1 — а2)зп2(2Кт, к)],

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

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

Интегрируя (5) вдоль контура области т, получаем выражения для геометрических и фильтрационных характеристик

0 . 5 0 . 5 0 . 5

§ ХвсвЬ = Д1, § УвсвЬ = Дв, § ФееЛь = 0.5Н, (6)

0 0 0

которые используются для нахождения неизвестных параметров конформного отображения а, в и М. Численным путем определяется монотонность функций, входящих в подынтегральные выражения левых частей уравнений системы (6), и таким образом устанавливается ее однозначная разрешимость. После этого рассчитываются координаты точек подземного контура плотины хвс(Ь), Увс(Ь), а также координаты криволинейной части водоупора хее(Ь), Уее(Ь), 0 ^ Ь ^ 0.5. Здесь Д1 = I — 1\, Дв = в — в\, Хвс, Увс, Фее - выражения правых частей (5) на соответствующих участках контура плоскости т. Полагая в уравнениях для координат Ь = 0.5, находим искомые размеры подземного контура плотины и криволинейного водоупора

= хвс (0.5), = увс (0.5), Ь = I + 12 = хее (0.5), Т = уее (0.5). (7)

Рассмотрим предельные случаи. Если в плоскости течения вертикальный отрезок АВ отсутствует, что соответствует слиянию точек А и В (параметры а = а = 0, = 0),

то, интегрируя уравнения (6), получаем следующие выражения для фильтрационных характеристик:

Л,_Я(1-Л0 Н[Е(У)-Х2К(У) + У}^

ПУ0Л ’ ПУ0Л ’ ПУ0 Л ’

(8)

1 = 11 + Л1 = щт^щун_ц,

ПУ0Л

где Е(А) - полный эллиптический интеграл второго рода при модуле Л = у/1 — /З2. Если в плоскости течения 2 отсутствует горизонтальный отрезок СО, что отвечает слиянию точек С и О (параметры Ь = в = 0, 1\ =0), то интегрирование уравнений (6) приводит к формулам (8) заменой в них параметров а на в, в на I и наоборот. Формулы (8) совпадают с формулами (10.9), (10.13), (10.16), (10.19), (10.22) и (10.24) [1, с. 197-200].

1.70

0.80

-0.10

0.20

1.00

0.50

0.00

IV

6.70

5.30

3.90

0.25 А/ о.ЗО 0.20

V

5.50

4.65

3.8

0.25 д/ 0.30

т

11

0.250

0.313 М 0.375

0.25

0.31

М 0.38

Рис. 4. Зависимости величин йг, (а) и Т, 12 (б) от ад (I), Н (II), Q (III), Д1 (/У) и Дй (У)

На рис. 1 изображена картина течения, рассчитанная при «0 = 1, Н = 2, ^ = 1.14, Д1 = 0.296 и Дв = 0.295 (базовый вариант). Результаты расчетов влияния определяющих физических параметров «0, Н, Q, Д1 и Дй на размеры I, 1\, в, в\, 12 и Т

представлены на рис. 4, I-У в виде зависимостей в\, 12 и 1\, Т от указанных параметров. В каждом из этих рисунков варьируется один из параметров «0, Н, Q, Д1 и Дв, а значения остальных фиксируются базовыми. Анализ графиков позволяет сделать следующие выводы.

Уменьшение скорости обтекания и увеличение действующего на сооружение напора воды приводят к росту всех размеров плотины I, 1\, в и в\, а уменьшение размеров криволинейного водоупора 12 и Т связано с увеличением параметров «0 и Н. Величины 1\ и в\, а следовательно, ширина и толщина плотины могут быть весьма значительными: из рис. 4, I следует, что возрастание скорости в 1.5 раза увеличивает ширину 1\ и толщину соответственно на 329 и 380.4%. Из графиков, приведенных на рис. 4, 111—У, видно, что при фиксированных значениях «0, Н, Q, Д1 и Дв глубина водоупора Т всегда превосходит ширину 12 в среднем на 10-20%.

Графики, приведенные на рис. 4, 1У, У, относящиеся к параметрам Д1 и Дв, отражают закономерность, которая является естественной с физической точки зрения: увеличению разности Д1(Дв) сопутствует убывание (рост) ширины плотины I и рост (убывание) ее толщины в. Так, с увеличением Д1 на 50% ширина 1\ уменьшается в 4.1 раза, толщина становится больше в 110.5 раз, подобное же увеличение параметра Дв приводит к росту ширины 1\ в 5.2 раза и уменьшению толщины в 14.7 раза.

С ростом параметров «0, Н, Д1 и Дв глубина водоупора Т и ширина 12 уменьшаются (хотя и незначительно, в пределах 1.1-1.7 раз) и увеличиваются с возрастанием фильтрационного расхода Q, причем существенно: на 147 и 119% соответственно. При этом величины Т и 12 могут быть весьма значительными и в рассматриваемых случаях превосходить не только параметры 1\ и в\, но и сами размеры плотины I и в соответственно. Так, из рис. 4, II вытекает, что при Н = 1.4 имеем 1\ = 0.285, I = 0.581, 12 = 6.678, значит, 12/1 = 11.5, а из рис. 4, 1У следует, что при Д1 = 0.2 получаем = 0.004, в = 0.299, Т = 6.495, следовательно, Т/в = 21.7. Таким образом, размеры 12 и Т превосходят ширину плотины I и ее толщину в на 1049 и 2072% соответственно.

Об одном случае обтекания шпунта Жуковского. Рассматривается течение жидкости под шпунтом АРЕ, когда на некоторой глубине Т имеется горизонтальный пласт, состоящий из непроницаемого участка ВС и хорошо проницаемого слоя СО, не содержащего напорных грунтовых вод (рис. 5). Грунтовые воды, обтекая шпунт с конечной скоростью Уе на его конце, поднимаются за ним на некоторую высоту РЕ и образуют свободную поверхность ЕО, на которую поступают инфильтрационные воды с интенсивностью е (0 < е < 1). Действующий напор Н, глубина залегания пласта Т, длина шпунта 5, а также скорость обтекания его на конце Уе (0 < Уе < е) считаются заданными. Задача состоит в определении положения кривой депрессии ЕО и, стало быть, размеров в (высота поднятия воды за шпунтом), Ь\ (длина отрезка непроницаемого участка слоя за шпунтом) и Ь2 (проекция свободной поверхности) при следующих краевых условиях:

где Q - искомый фильтрационный расход.

Впервые задача об обтекании шпунта рассматривалась Н. Е. Жуковским в статье «Просачивание воды через плотины» [7], в которой видоизмененный им метод Кирхгофа в теории струй был использован для решения задач фильтрации со свободной поверхностью. Здесь была введена специальная аналитическая функция, впоследствии получившая весьма широкое применение в теории фильтрации, с помощью которой

АВ : у = 0,р = —Н; ВС : у =

ОЕ : р = —у — Т,ф = Q + ех,

Т, ф = 0; СВ : у = —Т, р = 0, АРЕ : х = 0,ф = Q,

(9)

Рис. 5. Картина течения, рассчитанная при VF = 0.3, е = 0.6, T = 5, H = 5, S = 3

Н. Е. Жуковским дано исследование задачи об обтекании шпунта. С тех пор как сама функция, так и шпунт носят имя Жуковского [1]. Работа [7] открыла возможность математического моделирования задач со свободной поверхностью и положила начало исследованиям указанного класса фильтрационных течений [1, 3, 8].

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

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

ги(т) = \ZJtg7TT. (10)

Принимая во внимание соотношение (10) и поступая аналогично [9, 10], решение краевой задачи получим в следующем параметрическом виде:

dw Г sin 7гт dz cos itt

d,r sn(2AY, к)А(т) ’ dr sn(2AY, к)А(т) ’ ^

Рис. 6. Область комплексной скорости ад

здесь Д(т) = ^Ві28п2(2А'т, к) + В2, /Зі = у/1 — /З2, а = агЛі(Ур/л/є)/тт, /3 = вїі(2КЬ, к'), Ь - ордината точки А в плоскости т. Тогда модуль к рассчитывается из уравнения

р = К'/К = 2аі11іл/е/7г. (12)

В данном случае неизвестные параметры отображения Ь и М определяются в результате решения такой системы уравнений:

Ь 0.5

/ГорЗі = Б, / Фоп = Н, (13)

а 0

после чего вычисляются координаты точек свободной поверхности хаб (і) и уав (і), 0 ^ і ^ 0.5. Полагая в уравнениях і = 0.5, установим искомые размеры

Ь = хаб (0.5), З = Т - уав (0.5), (14)

а также

0.5 Ь

Ь1 = Ь2 - І ФаоЗі, = / ФооЗі, (15)

00

где Уор, Фсо, Фао и Фос - выражения правых частей (11) на соответствующих участках контура плоскости т.

На рис. 5 изображена картина течения, рассчитанная при Ур = 0.3, є = 0.6, Т =5, Н = 5, Б = 3 (базовый вариант). Результаты расчетов влияния определяющих физических параметров Ур, Н, Б, Т и є на фильтрационные характеристики представлены на рис. 7, І-ЇУ в виде зависимостей Q, З и Ьі, Ь2 от указанных параметров. Анализ данных графиков позволяет сделать следующие выводы.

1.0

3.5

1.5

-0.5

2.8

III

2.2

1.3

0.4

S 4.5 1.0

IV

3.2

1.1

Li .

d

-1.0 є 0.8 0.4

V

2.30

1.35

0.40

2.8

0.6

Q

Li

S 4.5

———

Xi

e 0.8

Q

Li

3.5

5.0

T 6.5 3.5

5.0

T 6.5

Рис. 7. Зависимости величин L2 и d (а), Q и L1 (б) от Vf (I), H (II), S (III), є (IV) и T (V)

Возрастание скорости Ур и напора Н увеличивает размеры Ь\, Ь2 и расход Q и, наоборот, уменьшает высоту поднятия воды за шпунтом 3. В то же время наблюдается совершенно противоположный характер изменения размеров Ь\, Ь2 и 3 при варьировании параметров Б и е: при снижении интенсивности £ всего в 2 раза размеры Ь\ и Ь2 возрастают на 1744 и 312% соответственно. Наиболее существенное влияние на глубину 3 оказывает напор Н и мощность пласта Т (рис. 7, V). При этом наблюдается линейная зависимость искомых величин от параметра Т.

Литература

1. Полубаринова-Кочина П. Я. Теория движения грунтовых вод. М.: Гостехиздат, 1952. 676 с.; 2-е изд. М.: Наука, 1977. 664 с.

2. Кочина И. Н., Полубаринова-Кочина П. Я. О применении плавных контуров основания гидротехнических сооружений // Прикл. математика и механика. 1952. Т. 16, вып. 1. С. 57—66.

3. Развитие исследований по теории фильтрации в СССР (1917—1967) / авт.-сост.: В. И. Аравин, А. В. Афанасьева, В. Д. Бабушкин и др.; председатель и отв. ред. П. Я. Полубаринова-Кочина. М.: Наука, 1969. 545 с.

4. Береславский Э. Н. О конформном отображении некоторых круговых многоугольников на прямоугольник // Изв. вузов. Математика. 1980. № 5. С. 3—7.

5. Береславский Э. Н. Построение контура постоянной скорости основания гидросооружения при фильтрации двух жидкостей разной плотности // Прикл. математика и механика. 1990. Т. 54, вып. 2. С. 342-346.

6. Береславский Э. Н. Определение подземного контура заглубленного флютбета с участком постоянной скорости при наличии соленых подпорных вод // Прикл. математика и механика. 1998. Т. 62, вып. 1. С. 169-175.

7. Жуковский Н. Е. Просачивание воды через плотины // Н. Е. Жуковский. Собр. соч. М.: Гостехиздат, 1950. Т. 7. С. 297-332.

8. Ведерников В. В. Теория фильтрации и ее применение в области ирригации и дренажа. М.; Л.: Госстройиздат, 1939. 248 с.

9. Береславский Э. Н. Гидродинамическая модель отжима пресными фильтрационными водами покоящихся соленых при обтекании шпунта Жуковского // Докл. АН СССР. 1998. Т. 303, № 4. С. 479482.

10. Береславский Э. Н. К задаче Жуковского об обтекании шпунта // Прикл. математика и механика. 1999. Т. 63, вып. 4. С. 603-610.

Статья рекомендована к печати проф. Л. А. Петросяном.

Статья принята к печати 24 сентября 2009 г.