П.Я. Полубаринова-Кочина - Теория движения грунтовых вод (1132345), страница 88
Текст из файла (страница 88)
е. коэффициент фильтрации не является постоянным), то напор Ь= — — = — +г л рд не будет удовлетворять уравнению Лапласа. Но если уравнения Дарси да да да и= — Ь вЂ”, п= — Ь вЂ”, ти= — Ь— дк' ду' да имеют место, то они показывают, что расход ЬЯ по-прежнему будет выражаться формулой (6.4). Однако условие (6.5) не будет выполняться; вместо него будем иметь соотношение бр тль — = й — = сопз1, аа и где Ь вЂ” осредиенный по вертикали коэффициент фильтрации. й 7. Метод конечных элементов для интегрирования эллиптических уравнений*). В последнее время большое применение находят вариационные методы решения уравнений эллиптического и параболического типа (см., например, Зенкевич и Цанг 1974).
Остановимся здесь на так называемом методе конечных элементов (для уравнения нестационарной фильтрации в главе Х1Х будет рассмотрен вариационно-разностный метод). Пусть дано уравнение — (Р, д ) + д (Рз д ) +та(х, У)=0, (7.1) где Р,(х, у) и рз(х,у) — заданные функции х, у. Требуется найти его решение в области О, ограниченной контуром 7., при условии, что на контуре значения Ь(х, у) известны: Ь(х, у) ~л —— 1(х, у).
(7.2) Будем искать минимум функционала $ $ ( о ~Р' ( д» ) + Рз ( др ) ~ тиЬ Ф с(х ату (7.3) и при условии, что Ь(х, у) удовлетворяет условию (7.2). ') При составлении етого параграфа был использован отчет В. А. Васильева. 596 хстхновившнвся движения гюнтовых вод (гл. хчш Как известно (см., например, Михлин 1970), если интеграл 7(и) = 3 3 г '(х у и д д ) ах ау о минимизируется, то необходимым условием является то, что функция г должна удовлетворять дифференциальному уравнению Эйлера дк ! д (ди/дх) ) ду ! д (ди/ду)) ди в той же области, а и(х,у) удовлетворяет тем же граничным условиям.
При достаточно широких условиях задача (7.!) — (7.2) эквивалентна задаче (7.3) — (7.2). Можно интерпретировать выражение (7.3) для 7 как энергию, запасенную жидкостью в данной области иа фиксированный момент времени. Рис. 334. Разобьем всю область В на треугольники произвольного вида (см. рис. 334) и пронумеруем их от 1 до М, а все вершины — от ! до М. Выделим какой-нибудь треугольник Л(й 1, Й) с вершинами М,(хну,), М;(хну;), Мь(хмуь). Обозначим значения Й(х,у) в узлах через Йь Й;, Йм В пределах каждого треугольника Й(х,у) считается представимой линейной функцией координат: (7.4) Й (х, у) = а, + а,х + азу Следовательно, Й(хь у~) =Й, = а, + а1х, + азуо Й (хй у() = Й( — — а, + а2х( + азуй Й (хм уь) = Й, = а, + а,хь+ азуа. $7] метод конечных элементов 597 Из этой системы уравнений находим а, и подставляем в (7.4).
Получим Ь(х, у) = Ь!А+ Ь!!Ь!+ Майн (7.5) где И,= — (а,+ Ь,х+ с,у) (э=!, 1, Ь). (7.6) 1 При этом а! — — х!уь — хьу1, а, = хьу, — х,уы аь — — х,у, — х!у!, (7.7) Ь =у! уы Ь1= уь — у! Ь~=у! у! (7.8) с; =ха — х1, с! =х! — хы сь=х! — х!, (7.9) ! Ь = —, [(х, — х,) (у! — уь) — (х, — хь) (у, — у!)]. (7.10) Теперь интеграл 7 по треугольной области будет функцией Ь,, Ь1, Ь!д 7 (Ь (х, у)) = 7'(Ь,, Ь1, Ьь).
Чтобы найти минимум 7' по переменным Ь„Ь1, Ьы полагаем дУ' д!' дУ' — — — 0 да! да! даь Получаются основные уравнения метода конечных элементов для треугольной области, которые можно записать в матрич- ной форме: (7.11) Здесь [д] — квадратная матрица (хм+ аи) хц икь [Ы] = х!! (У!! + а!ь) и!ь, (7.12) уы ум (уы + а1ц) ! кц = ча (РАь! + Рзс!с!), 1 йи ча (Р!Ь!Ьь + Р~с!с ) 1 а!х = — (р! Ь!Ьь + рзс!сь) (7.13) (функции р,(х, у) и рз(х, у) в пределах каждой треугольной ячейки считаются постоянными), Далее, имеем (Ь), (бв!3)— 598 гстлновившиеся движения гегнтовых вод 1гл. хчп! матрицы-столбцы: (й) = й! (7.!4) ! 1 — би!! 3 (З 1 зб (7.15) 1 — Йвь з [И] = [А] ' Х [В] Х [А] (7.17) где Ь! сс Ь! с! (7 18) Ь„ сь с! с! «ь [А]-' ' Ь, Ь, Ь, сЬ с с с О О О [В]= о р, о [А] а (7.!9) Для составления системы уравнений для всех неизвестных й, в области Р нужно объединить матрицы проводимости для всех элементов области.
Получается система уравнений вида [Р] Х (г!) = (! !), (7.20) где Ь| Ь2 с1 Ч2 В развернутом виде уравнения (7.11) имеют вид Ии (й! — й!) + Ись (йь — й!) = О! Иь! ("» — "!) + Ис! (и! й!) = 9! (7.16) И!ь (й! — йь) + Ись (й! — йь) = Чь где о! —— Ош!!3 и т. д. Структура этих уравнений позволяет трактовать коэффициенты Ип как «проводимости» некоторых фиктивных трубок, по которым жидкость якобы перемещается от одного узла треугольника к другому. Матрицу (7.12) называют матрас!ей фильтрационных проводимостей треугольного элемента.
Эта матрица симметрична: Ио = Им,... Для программирования ее записывают в виде метОД кОнечных элементов 699 — (3„+ 3,,) И!+ д!48, + а„'й =9„ 3,!й! — (841,. + д4!3) й! + 343йб = 941, ЯА+ 334й! — (И + 334) йб = 93. (7.21) Рассмотрим схему из шести треугольных элементов с семью вершинами (рис. 334,6) и составим для них шесть систем уравнений вида (7.21). Выпишем лишь какие-нибудь две из них: например, для первого элемента с узлами 1, 3 и 4: (в!3+ в14) 1+ 613 3+ в!4 4 41! л31 1 (831+ л34) 3+ л34 4 !13' л41 1+ л43 3 !,л41+ л43) 4 л4' (7.22) для четвертого элемента с узлами 3, б н б: (84 1 646) ь 1 34 13 ! ~У вЂ”,14 л53 3 !л543+ л56) 5+ л55 Б !15' в463 3+ йб5 5 4663+ 6465) б '16' (7.23) Всего получим 18 уравнений с 7 значениями й: 151, ..., йь Их группировка в систему из 7 уравнений производится так.
Три уравнения, правые части которых содержат 411, 413 и дь встре- чающиеся лишь по одному разу (что видно и по ячейкам с но- мерами 1, 6 и 6), выписываются без изменений, так что для первого уравнения объединенной системы 7 Х ~1.Ь.= Ф л 1 будет, согласно первому из уравнений (7.22), !3 М Так же составим вторую и седьмую строки матрицы 6. Далее, заметим, что величина дб (с разными верхними знач- ками) встречается четыре раза: во втором из уравнений (7.22), Составление матрицы (41] является трудоемким, оно производится по тем же правилам, которые разработаны в электротехнике при составлении баланса токов в замкнутых и разветвленных цепях (Бессонов 1967). Поясним на простом примере, как строится система уравнений (7.20).
Предварительно перепишем систему уравнений (7.16), снабдив в ней все д!3 значком ! в знак принадлежности их ячеике с номером 1: ооо УСТАНОВИВШИЕСЯ ДВИЗКЕНИЯ ГРУНТОВЫХ ВОД !ГЛ. ХРЕН в первом из уравнений (7.23) и еще в двух невыписанных уравнениях. Сложим почленно четыре соответствуюгцих уравнения.
В правой части объединенного уравнения получим '~з = Чз+ 032+ уз+ Чз. Суммируя коэффициенты при Ьь Ь,, ..., Ь„найдем 33 (ВЗ! + 334 + 636 + 634 + 633 + 636 + Взг + УЗЗ)! 34 з'34 + » 34' 35 »435 + з 635' 36 » 36 + з'м' 47 Нетрудно заметить правило составления элементов 6м, Для диагональных элементов имеем 66 = — ~' у,'„. 3 ! Недиагональные элементы можно составить, рассматривая на рис. 334 отрезок (1, Ь): если он примыкает только к одной ячейке 1, то 6, =д!3, если к двум ячейкам 1! н 13, то 6,„= д',.„+ д,'*3, если узлы 4 и Ь не лежат на одном отрезке, то 64А = О.
Для нашего примера можно записать систему уравнений (7,20) в виде Матрица 6 симметрична; 646= 634, она имеет ленточную структуру, т. е. 60 = О, как только (1 — 1! ) Ш, где Ш вЂ” ширина ленты матрицы проводимостей. Если граничное условие отличается от (7.22) — обычно оно задается в виде Ь „— соз (и, х) + ЬУ вЂ” соз (и, у) + д (х, у) + ОЬ = О, да да то в выражении (7.3) для 7(Ь) добавляется интеграл по контуру Л, ограничивающему область 6: дЬ 51~+ — ( ~Ь'415 2 ! и ищется минимум такого расширенного функционала.
0„ О 0„ О 022 022 02! 0зг 022 О, О 0 О / 032 032 О О ~ 0 о о о О, ~ О О О, 024 023 044 043 034 0и 0„ 6м 024 О О О О 024 1 О 0»а 0и 0 0 044 О О О, а! Ьг Ьз Ь» «а Ьа аг О! 4»'г 443 9» Оа 444 О! во1 основнык соотпоншния метода эгдх эн На рис, 335 представлены результаты расчета линий равного напора, выполненного О. Зенкевичем с соавторами (см. Зенкевич и 11анг 1974) для обтекания плотины в сильно неоднородной анизотропией среде с криволинейными границами. Главные Рис. 333. оси анизотропии каждого слоя вместе с коэффициентами фильтрации вдоль них указаны на рисунке в виде построенных в одинаковом масштабе взаимно ортогональных векторов.
Б. МЕТОД ЭЛЕКТРОГИДРОДИНАМИЧЕСКИХ АНАЛОГИЯ (ЭГДА) Э 8. Основные соотношения метода ЭГДА. Между установившимся движением грунтовых вод и движением электрического тока существует аналогия, выражающаяся в том, что оба явления описываются с помощью дифференциальных уравнений одного и того же вида с одинаковыми граничными условиями. Эта аналогия была использована Н.
Н. Павловским, который еще в !918 г. применил ее для экспериментального решения ряда задач, встречающихся при изучении движения грунтовых вод (Павловский 1922). Предположим, что у нас имеется тонкий слой жидкого или твердого проводника (теперь получила широкое применение электропроводная бумага), толщина которого г может быть переменной. Дифференциальные уравнения установившегося б02 У!ЛВНОИ1ВШНИСЯ ЛВИЛОЗ1У!Я ГВУНТОВЫХ ВОД (ГЛ. ХУН1 движения электрического тока ду юх == — с —, х дх имеют вид дУ 1 = — С вЂ”, д (г1с) =о, ду (8.1) д (г1х) + (8.2) дх где 1„, 1„ — компоненты плотности тока, У вЂ” электрический по- тенциал или функция напряжения, с — коэффициент электро- проводности, х, у — координаты точек электрической модели. Введем функцию тока Ч' с помощью равенств дЧ' .
дЧ' — = — г1, — = 21' У ду х. (8.3) Сравнивая эти выражения с (8.1), получим дУ 1 дЧ" дк 1 дЧ' дх сг ду ' ду сг дх ' Обратимся теперь к уравнениям теории фильтрации да да и= — Й вЂ”, о= — Й вЂ”, дх ' ду ' — + — =О, дх ду (8.6) (8.4) где Ь вЂ” напор. Введя функцию тока ф, найдем зависимость между Ь и !р: да 1 д!В дл ! дЧ! (8.7) дх й ду ду й дх Сравнение (8.7) и (8.4) показывает, что линиям тока грунтового потока отвечают линии электрического тока, а линиям равного напора — линии равного электрического потенциала. Если смоделировать дно по условию Й = сг и г построить область, заполненную проводником, подобную области грунтового потока (рис.