CC BY
16
УДК 519.63:532.5
А.П. БЛИШУН, М.В. СИДОРОВ, И.Г. ЯЛОВЕГА
ЧИСЛЕННЫЙ АНАЛИЗ СТАЦИОНАРНЫХ ФИЛЬТРАЦИОННЫХ ТЕЧЕНИЙ СО СВОБОДНОЙ ГРАНИЦЕЙ СТРУКТУРНО-ВАРИАЦИОННЫМ МЕТОДОМ
Рассматриваются применение вариационного метода Ритца и структурного метода Я-функций для численного решения задачи расчета фильтрации в области со свободной границей. Уточнение уравнения свободной границы проводится итерационно. Для модельной задачи с известным точным решением проводится вычислительный эксперимент.
Введение
Актуальность исследования. Моделирование и численный анализ течений жидкости в пористой среде имеет большое прикладное значение. Задачи такого класса возникают в инженерном деле (например, при проектировании гидротехнических сооружений), в экологии (рациональное природопользование), при разработке комплекса мер по защите населения в чрезвычайных ситуациях (например, во время паводков) и т.д. Основные сложности, возникающие при численном анализе таких систем, обусловлены разнообразием режимов течений и типов грунта, сложной геометрией области течения и возможным наличием в области фильтрации участков свободной границы. Основные методы моделирования фильтрационных течений рассмотрены в [1 - 5]. Таким образом, разработка новых и усовершенствование существующих методов математического моделирования фильтрационных течений является актуальной прикладной задачей.
Цели и задачи исследования. Целью настоящего исследования является применение для численного анализа фильтрационного течения со свободной границей структурного метода Я-функций и вариационного метода Ритца. Для достижения поставленной цели необходимо решить следующие задачи:
- основываясь на построенной структуре решения краевой задачи, свести ее к вариационной задаче о минимуме функционала энергии;
- разработать алгоритм решения вариационной задачи методом Ритца;
- провести вычислительный эксперимент для модельной задачи с известным точным решением.
1. Постановка задачи
Рассмотрим двумерный случай. Предположим, что оси координат являются главным направлением тензора коэффициента фильтрации. Будем решать смешанную краевую задачу, связанную с определением подземной части гидросооружения:
kll du ]__^
Sx i 5x) Sy
(
k —
k22 -
. ^y.
л
= 0 в Q ,
uIae = ui' uIfc = u3
Su
cN
-c0u
= Уо(х,У).
ABCUA'B'C'
(l) (2) (3)
du , Su , . , Su , . где u - напор; -= k11—cos(n,x) + k22—cos(n,y) - производная по конормали; n -
dN Sx 5y
единичный вектор внешней нормали; u1, u3 - постоянные величины. Часть границы области A'B'C' (рис. 1) - неизвестная линия, которую обозначим через SQ„. В частном случае, когда ABC и A'B'C' - непроницаемые границы, имеем с0 = 0 и у0 = 0 .
Предполагаем, что точки A' и C' неизвестные. Согласно теории на свободном участке границы cQ. задается дополнительное, так называемое «геометрическое условие» (или «кинематическое условие»)
u = S(s) на SQ., (4)
где означает заданный напор на подземном контуре гидросооружения.
A
A'
В
Рис. 1
2. Метод численного анализа
Будем предполагать, что уравнение (1) эллиптическое невырождающееся, т.е. при всех
\2 е Я \{0} в □ выполнено неравенство
к1А2 + к2212 >Ц+12) в П, > 0-
Введем следующие обозначения:
Au = _
Sx
du
dx) 5y i Sy
du
Nu = kn ducos(n,x) + kii ducos(n,y)
dx dy
dQ. = A'B'C' - свободная граница, 5Q1 = AA' U C'C , dQ2 = ABC . Продолжая краевые условия внутрь области Q согласно методу R-функций [6], получаем, что краевую задачу (1) - (4) можно записать в виде
Au = 0 в Q, (5)
u| dQ =^1^ (6)
Nu +Cul8Q2U0Q, =Ню2U0Q , (7)
u = S(x,y) на dQ.. (8)
y
Здесь Y = ЕСу0, c = ECc0, ф= ^ , + Ul„m , фЦ =
ю + ю 1
Íu1 на AA ',
u3 на С' С
а функции ю' (x,y) и
ю "(x,y) строятся с помощью R-функций и удовлетворяют условиям
1) ю' (x, y) = 0 на AA ', ю "(x, y) = 0 на С ' С ;
2)ю' (x,y) > 0 в QU A В ' С 'UС ' Си AB^ ю '' (x,y) > 0 в QU A ВС 'U AA'U ABC;
3)
дю ' дп
дю '
= -1
AA' дп
= -1, п - внешняя нормаль к 5Q .
Если 50, была бы заданной границей, то краевая задача (5) - (7) могла бы быть сведена к нахождению минимума соответствующего квадратичного функционала.
В данном случае мы сразу не можем находить минимум полученного функционала, так
как область о, по которой ведется интегрирование, и часть её границы 50, неизвестны. Поэтому для полного решения задачи необходимо использовать и условие (8).
Для решения поставленной задачи сформируем итерационный процесс следующим образом. Пусть получено к -е приближение к свободной границе 50(к) . Тогда получим конкретную область 0(к) с известной границей. Определим функцию и(к)(х,у) как решение задачи
Аи(к) = 0 в 0(к), (9)
и(к)| =ф|
яо. Ч1эп, ,
(10)
Nu(k) +ou(k)l
=i
1эо2 иэо<к) '1эо2иэо<к) . (11)
Для решения задачи (9) - (11) воспользуемся методом Я-функций. Пусть уравнение 30(к) представлено в виде у = Ь(х), причем у - Ь(х) > 0 в 0(к) и 301 и ЗО2, тогда
ю^у) =
y - h(x)
л/1 + (h ' (x))2
обладает следующими свойствами:
1) ю(к) (x, y) = 0 на dQ(k); 2) a(k)(x,y) > 0 в Q(k) U дЦ U dQ2; 3)
(k) (k)
(k)
дю«
дп
= -1.
SQ<k)
Пусть ю1(х,у) и ю2(х,у) обладают следующими свойствами соответственно:
1) ю1(х,у) = 0 на 501; ю2(х,у) = 0 на 502;
2) ю1(х,у) > 0 в О и д02 и 50(к); и2(х,у) > 0 в О и д01 и 50(к); 5Ю2
дю1
3) —1
дп
= -1;
dQ1
дп
= -1.
В работе [6] показано, что структура решения краевой задачи (9) - (11) имеет вид
u(k) = ф+ю1Ф0
ю1(8Бса + 5ю) + Sea
[-SD^- SD 1(ю1Ф0) + сф + сю1Ф0 -Y + éo Ф1]
, (12)
где 5 = 1 -(1,Vc5)>0, S= J| kn—
dco ^2 Л дю ^
k22
ду
22
v /
_2 _i д^, д^, (k)
ю , D1u =—11 +--12, ю =ю2 ла ю* ,
дx ду
1 = (, I2 ) =
' k11 дю k22 дю л S д^' S ду
Ф0 и Ф1 - неопределенные компоненты структуры. Здесь
ла - знак R-дизъюнкции; u ла v = —1—(u + v - Vu2 + v2 - 2auv).
1 + a
SQ
Возьмем Ф1 = 0 и обозначим Ф0 = ф(к). Перепишем (12) в виде
где
= Ф + -
и(к) = g + В(Ф(к)), 66 -[-80; ф+аф-у]
ю1 (8Б11 са + 5со) +
В(Ф(к)) = с;Ф(к) +-г"6-[-8Б1Т(ш1Ф(к)) + асо1Ф(к)]
ю1(8Б1СЙ + 56) +
В задаче (9) - (11) сделаем замену: и(к) = g + у(к).
Тогда для у(к), в силу построения функции g, получим задачу с однородными краевыми
условиями:
(к)
Лу(к) = -Ag в а
у(к) = 0 1эп, :
Ку(к) + ау(к) (к) = 0
(13)
(14)
(15)
Оператор краевой задачи (13) - (15) действует на области определения Бл = {у е С2(а(к)) пС1(а(к)), Лу е Ь2(а(к)), уЦ = 0, Ку + ау|да^ = 0} с Ь2(а(к)) и является на ней положительно определенным [7]. Энергетическое пространство этого оператора получается пополнением БЛ в норме, порожденной скалярным произведением
[и,у] = Ц
а(к)
, ди ду , ди ду
кц + к22
дх дх ду ду
ахау + | аиуаэ, ||и||л [и,и].
эп2 иап(к)
Тогда в соответствии с теоремой о функционале энергии решение задачи (13) - (15) можно найти как:
у(к) = arg
м \я
уеНл [Д
-2.С.Г
а(к)
ду дх
дх
(ду V
Чду У
д (
к
V
ахау + | ау^ -
ЭП2 исП(;к)
дg
ахау}
Для нахождения минимума этого функционала воспользуемся методом Ритца [7]. При-
т
ближенное решение будем искать в виде ут0 = ^а^^, где {ф1} - координатная последо-
1=1
вательность. Пусть {т1} полна в W2(Q(k)). Тогда в качестве {ф1} возьмем
Ф1 = В(т1) = ю1х1
-[-8Б11 (ю1х1) + аю1х1]
ю1 (8Б11 сс + 56) + 8са
Коэффициенты а1к),..., аПк) найдем из системы линейных алгебраических уравнений
т
Xс(к)[ф1, ф^ = (-Лg, фД ] = 1,2,...,т.
1=1
Найденная функция и^к = g + у^к является приближенным решением задачи (9) - (11) и условию на свободной границе (8) не удовлетворяет.
Условие (8) используем для уточнения уравнения свободной границы. Уравнение свободной границы да(к+1) для следующей (к +1) -й итерации ищем в виде:
у(к+1) = С0к+1) + £ сГк+1)хг.
Г=1
где т - натуральное число; с(к+1) (г = 0,1, 2, ...,п) - неизвестные коэффициенты.
Пусть а(к) и с(к) - абсциссы точек А ' (к) и С ' (к). Тогда положим ^к) = |с(к) - а(к)|.
Возьмем натуральное число N такое, что N > п и на отрезке с концами а(к),с(к) зададим
N — 1 п сетку х1 = а(к) +--d(k), 1 = 1, 2,..., N . Положим у(к+1) = I с(к+1)хГ . Кроме единичной нор— 1 г=0
мали
п = (п1,п2) ведем в рассмотрение конормаль v = (v1,V2) согласно правилу: V =
к11п1 8 ,
к22п2
, где 8 = [(к11, п1)2 + (к22, п2)2]2, п1 = со8(п, х), п2 = со8(п, у). Для вычисления координат точек Р., лежащих на 50(к+1), имеем [6]:
X. = х. + cos(v,x)
ду
— 0(Р.)
У1 = у. +
итк)(р.)—ад)
дОД)
ду
— 0(Р.)
(17)
Для нахождения коэффициентов в уравнении кривой 50(к+1) потребуем, чтобы кривая
каким-то образом определялась через множество точек Р. = (х., у.), 1 = 1,2,...,N. Для этого воспользуемся точечным методом наименьших квадратов: потребуем, чтобы выра-
жение
* = 1
у 1—I сГЧ г(х 1)
имело минимальное значение, где значения х. и у.
берутся из (17). После определения всех параметров из (16) кривая у(к+1) будет известной, значит, имеем известную область О(к+1:1. На следующей (к + 2) -й итерации строится приближенное решение итк+2) для новой области О(к+1:1. Итерация проводится до тех пор, пока " расстояние" между дО(к) и 50(к+1) не станет достаточно малым:
т1ахч/(х(к+1) — х(к))2 — (у(к+1) — у(к))2 , где в> 0 - наперед заданное достаточно малое число.
3. Модельный пример
Вычислительный эксперимент был проведен для задачи
д ( 1 ди ^ д ( 1 ди ^
дх
1 — у2 дх ) ду ^ (1 + х)2 ду
1 ди
/
,,1
АС
1 ди
= 0
в О,
Л 2 я.-с0^п,х) 4Л >2,
1 — у2 дх (1 + х)2 ду
cos(n, х)
1 — у2 дх
Час = 1 — у ,
cos(n,y)
1 а,
(1 + х)2 ду I 3
и1вс = 4.
3 + 3х + 16у2 + 64у2(1 + х)3
cos(n,y)
=0
(18)
(19)
(20)
(21)
(22) 19
вС
Ав
СI 0;
, 2, да1
Л
да.
да 2
Рис. 2
В
Область а представлена на рис. 2. Точку С считаем известной, положение точки в неизвестно, известно лишь, что она находится на прямой Ох , ВС - есть свободная поверхность. Точное решение
этой задачи имеет вид и(х,у) =1——, у =1 >/1 - 3х . В
1 + х 2
качестве начального приближения к уравнению свободной границы возьмем прямую, соединяющую точки
х 1 ^ Г}(0) ( 2 0 и(0)(х) = 5 1
С| 0, 2 1 и Вт I -, 0 1: у = Ь^(х) = --х + -. Построим
решение задачи (18) - (22) в области а(0) = ЛСВ(0) по алгоритму, описанному в разделе 2. Для задачи (18) - (22)
б!0)(х,у) =
- у - 1,25х - 0,5 ^2,5625
с1(х,у) = х> с2(х,у) = у> ®(х,у) = у лс
- у - 1,25х - 0,5 ,/2,5625 :
ф = 1 -у2, а = 0, у = -
3 + 3х + 16у2
- I ,-ту
^9 + 64у2 (1 + х)3
' - у - 1,25х + 0,5
8 =
л/2,5625
( 1 дСо 1 дсо ^
1 дсс
^ ( 1 дсс V
1 - у2 дх У V (1 + х)2 ду
-со
1 - у2 дх . (1 + х)2 ду
Б1и =1
1 ди дСо
1 ди дсо
1 - у2 дх дх (1 + х)2 ду ду
Структура решения краевой задачи (18) - (21) имеет вид
и(0) (х, у) = ф+сФ0 +
ю^Б1 со + 5сс) + 8со
[-8В1ф-8Б1(ю1Ф0) + аф + асо1Ф0 -у + со¥2]
Неопределённую компоненту Ф0 = Ф( ) аппроксимируем выражением вида
т т
Ф(0) «Фт0) = X а1]Т1] = X ачхУ.
1+-=0 1+-=0
Числа а-, 1 + j = 0, 1, ..., т , найдем из условия минимума соответствующего функционала энергии. Координатные функции возьмем в виде
фц(х,у) = с1Т1, -
С® СТЛ1/ \
ю1(8Б1со + 5со) + 8са '
где
5 = 1 -
1 (дсо
( дш^2
1 - у2 V дх У (1 + х)2 V ду
Используя значение и^0) (х, у), уточним уравнение свободной границы с привлечением
условия (22). Будем искать уравнение да,1) в виде
у(1) = с<» + с1(1)х + о(21)х2
Возьмем N = 5 , х1 =
5 -1 10 :
1 = 1,2,..., 5 .
На дП(,0) ^(^х) =
1 1 дю
81-у2 дх
cos(v,x) =
1 1
8О<0)
лежащие на дО(0): у, = -—х, + 1.
* 4 1 2
По формулам (17) находим
1 1 дю
81-у2 дх
(х1.У1)
3
и—"Чвд) - 4
дю
1
(х1.У1)
8(1 + х)2 ду
У1 = У1 +77
Р = (х,,У,) - точкИ,
оО<0)
1 1
дю
8(1 + х)2 дУ
и —0)(х1,У,) - -4
(х1.У1) --у
8
(х1.У1)
Тогда коэффициенты с0), с(), с(2) в (23) найдем из условия
£[У 1 - (с01)+с
с°)х ■
(1)~2 - с2 хх
)]2
■ Ш1П
с, ,с2 ,с3еК
Сходимость с точностью в = 10 3 достигнута на 7-й итерации. Сравнение с точным
решением показало, что и
(7)
- и (7) = 0,12-10-
Нь2(п(7))
Выводы
Впервые на основании структурно-вариационного метода разработан метод численного анализа стационарных фильтрационных течений в областях со свободной границей. Разработанные нами методы могут быть использованы при проектировании сложных гидродинамических сооружений: плотин, дамб и т.п. Этим и определяется научная новизна и практическая значимость полученных результатов.
Список литературы: 1. БуракЯ.Й., Чапля С.Я., Чернуха О.Ю. Континуально-термодинамiчнi моделi мехатки твердих розчитв. К.: Наук. думка, 2006. 272с. 2. ЛяшкоИ.И., СергиенкоН.В.,МистецкийГ.Е., СкопецкийВ.В. Вопросы автоматизации решения задач фильтрации на ЭВМ. К.: Наук. думка, 1977. 288 с. 3. Ляшко Н.И., Великоиваненко Н.М. Численно-аналитическое решение краевых задач теории фильтрации. К.: Наук. думка, 1973. 264 с. 4. ШаманскийВ.Е. Численное решение задач фильтрации грунтовых вод на ЭЦВМ. К.: Наук. думка, 1969. 376 с. 5. Вабищевич П.Н. Численные методы решения задач со свободной границей. М.: Изд-во Моск. ун-та, 1987. 164 с. 6. БлишунА.П., СидоровМ.В., ЯловегаИ.Г. Математическое моделирование стационарных фильтрационных течений со свободной границей методом Я-функций // АСУ и приборы автоматики. 2010. №2 150. С. 18-27. 7. Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970. —12 с.
Поступила в редколлегию 05.06.2010 Блишун Александр Павлович, студент группы ПМм-09-1 фак-та прикладной математики и менеджмента ХНУРЭ. Научные интересы: математическое моделирование, численные методы математической физики. Увлечения и хобби: покер, футбол. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.
Сидоров Максим Викторович, канд. физ.-мат. наук, доц. каф. прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, численные методы, математическая физика, теория Я-функций и её приложения, стохастический анализ и его приложения. Увлечения и хобби: всемирная история, история искусств. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.
Яловега Ирина Георгиевна, канд. техн. наук, ст. преп. каф. прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, анализ динамических систем, стохастический анализ и его приложения. Увлечения и хобби: оригами. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. (057) 7021436.
