Численный анализ стационарных фильтрационных течений со свободной границей структурно- вариационным методом

УДК 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'

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

где означает заданный напор на подземном контуре гидросооружения.

A

A'

В

Рис. 1

2. Метод численного анализа

Будем предполагать, что уравнение (1) эллиптическое невырождающееся, т.е. при всех

\2 е Я \{0} в  выполнено неравенство

к1А2 + к2212 >Ц+12) в П, > 0-

Введем следующие обозначения:

Au = _

Sx

du

dx) 5y i Sy

du

Nu = kn ducos

dx dy

dQ. = A'B'C' - свободная граница, 5Q1 = AA' U C'C , dQ2 = ABC . Продолжая краевые условия внутрь области Q согласно методу R-функций [6], получаем, что краевую задачу (1) - (4) можно записать в виде

Au = 0 в Q,

u| dQ =^1^

Nu +Cul8Q2U0Q, =Ню2U0Q ,

u = S

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(к) с известной границей. Определим функцию и(к)(х,у) как решение задачи

Аи

яо. Ч1эп, ,

Nu(k) +ou(k)l

=i

1эо2 иэо<к) '1эо2иэо<к) .

Для решения задачи (9) - (11) воспользуемся методом Я-функций. Пусть уравнение 30(к) представлено в виде у = Ь(х), причем у - Ь(х) > 0 в 0(к) и 301 и ЗО2, тогда

ю^у) =

y - h

л/1 + (h ' (x))2

обладает следующими свойствами:

1) ю

дю«

дп

= -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]

,

где 5 = 1 -(1,Vc5)>0, S= J| kn—

dco ^2 Л дю ^

k22

ду

22

v /

_2 _i д^, д^,

ю , 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) - (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

ду

— 0

У1 = у. +

итк)

дОД)

ду

— 0

Для нахождения коэффициентов в уравнении кривой 50(к+1) потребуем, чтобы кривая

каким-то образом определялась через множество точек Р. = (х., у.), 1 = 1,2,...,N. Для этого воспользуемся точечным методом наименьших квадратов: потребуем, чтобы выра-

жение

* = 1

у 1—I сГЧ г

имело минимальное значение, где значения х. и у.

берутся из (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

1 — у2 дх

Час = 1 — у ,

cos

1 а,

(1 + х)2 ду I 3

и1вс = 4.

3 + 3х + 16у2 + 64у2(1 + х)3

cos

=0

(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: у = Ь^(х) = --х + -. Построим

решение задачи

б!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 дх

3

и—"Чвд) - 4

дю

1

8(1 + х)2 ду

У1 = У1 +77

Р = (х,,У,) - точкИ,

оО<0)

1 1

дю

8(1 + х)2 дУ

и —0)(х1,У,) - -4

(х1.У1) --у

8

Тогда коэффициенты с0), с(), с(2) в (23) найдем из условия

£[У 1 - (с01)+с

с°)х

(1)~2 - с2 хх

)]2

Ш1П

с, ,с2 ,с3еК

Сходимость с точностью в = 10 3 достигнута на 7-й итерации. Сравнение с точным

решением показало, что и

- и (7) = 0,12-10-

Нь2

Выводы

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

Список литературы: 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.