O pewnej interpretacji metod energetycznych dla problemów filtracji ustalonej

I STOSOWAN A 1, 19 (1981)

O PEWN EJ IN TERPRETACJI M ETOD  EN ERG ETYCZN YCH  DLA PROBLEM ÓW F I LTR AC JI U STALON EJ

B O G D AN  W O S I E W I C Z   ( P O Z N A Ń ) 1. Sformułowanie problemu I zależ noś ci podstawowe

W pracy wykaż emy, że poszukiwanie rozwią zania problem ów filtracji ustalon ej m e-todami energetycznymi przez minimalizację  funkcjonał u energii jest równ oważ ne z m i-nimalizacją  metodą  najmniejszych kwadratów waż onych bł ę dów wyznaczenia skł adowych wektora prę dkoś ci filtrują cej cieczy.

Zał oż ymy, że filtracja odbywa się  w trójwymiarowym anizotropowym obszarze fil-tracji  F a osie ukł adu współ rzę dnych pokrywają  się  z gł ównymi osiami an izotropii oś rodka (rys. 1). Przy takich zał oż eniach problem filtracji sprowadza się  d o wyznaczenia funkcji

Rys. 1, Oznaczenia (piaski obszar filtracji dla uproszczenia rysunku).

h(x, y, z), opisują cej rozkł ad cis'nienia piezometrycznego, speł niają cej w obszarze filtracji

równanie róż niczkowe , , . d / , dh\  d / , Bh\  d / . Bh\  , _ n . s „ U) - a-  \k x- ir- 1 +  — \k >- T -  + - T" \k *- x-   + 2 =  0. (x,y,z)eV, dx \  8x I oy \  dy I Bz \  dzj a n a brzegu S mieszane i niejednorodne warunki brzegowe (2) h(xyz) = g, (x,y,z)eS1, oraz

W zależ noś ciach powyż szych przez kx, ky, k. oznaczono współ czynniki filtracji odpo-wiednio wzdł uż osi x, y, z, przez u, v, w skł adowe wektora prę dkoś ci wzdł uż tych osi, przez Q, g, q odpowiednio wydatek ź ródeł wewnę trznych, zadane na brzegu St ciś nienie piezometryczne oraz wydatek przez brzeg S2 a przez a, /?, y ką ty jakie tworzy normalna do brzegu S z osiami współ rzę dnych.

Zwią zek pomię dzy skł adowymi wektora prę dkoś ci a ciś nieniem piezometrycznym dan y jest za pomocą  prawa D arcy'ego

oh u =  — kx- z~, 8x , 8h (4) v=- ky- ^ , . dh

w = —k.

_Bz ' Alternatywnie zagadnienie sprowadza się  do znalezienia takiej funkcji h (x, y, z), która minimalizuje funkcjonał  energii i speł nia zasadniczy dla tego problemu warunek brzegowy (2). W metodach energetycznych stosowanych do zagadnień filtracji ustalonej poszukuje się  przybliż onego rozkł adu ciś nienia piezometrycznego h w postaci

(6) h(x,y,z)

gdzie liczby Ht{i = 1,2, ...,ń ) są  szukanymi współ czynnikami (w metodzie elementów skoń czonych są  t o wprost wartoś ci ciś nień piezometrycznych w wyróż nionych punktach obszaru filtracji) a znane funkcje N- ,(x, y, z) stanowią  zbiór funkcji bazowych (układ współ rzę dnych) w rozpatrywanej przestrzeni funkcyjnej. Po wyznaczeniu współ czynników

Hi przez minimalizację  funkcjonał u (5) moż emy także wyznaczyć przybliż one wartoś ci

skł adowych wektora prę dkoś ci u, v, w n a podstawie prawa D arcy'ego.

2. Funkcjonał  bł ę du wyznaczenia skł adowych wektora prę dkoś ci i jego minimalizacja Jako miarę  bł ę du przybliż onego rozwią zania problemu filtracji uznać moż na nastę -pują ce wyraż enie

(7) e = - r-  (u- u)2  + ~j—(v-  i) 2  + — (w -  ii')2 . ICX Ky Kg

Zależ ność (7) jest waż onym bł ę dem kwadratowym wyznaczania skł adowych wektora prę dkoś ci. Jako wagi przyję to odwrotnoś ci współ czynników filtracji. Wartość bł ę du e jest oczywiś cie funkcją  poł oż enia a także zależy od sposobu aproksymacji skł adowych

Jako najlepsze rozwią zanie, w sensie metody najmniejszych kwadratów, m oż na uzn ać ie, które cał ce takie, które cał ce (8) E = nadaje najmniejszą  wartoś ć.

/ /  /

nadaje najmniejszą  wartoś ć.

Wykaż emy teraz, że wyznaczenie liczb Hx przez minimalizację  funkcjonał u (8) jest równoważ ne minimalizacji funkcjonał u (5).

Wyraż ając u, v, w poprzez funkcję  h(x, y, z) za pomocą  formuł  (4) moż emy funkcjonał (8) zapisać w sposób nastę pują cy

/

Ostatnia cał ka w powyż szej zależ noś ci ma dla każ dego konkretnego zadania stał ą - wartość (choć oczywiś cie inną  dla róż nych zadań ), w ż aden bowiem sposób n ie zależy

ona od wyboru przybliż enia funkcji h(x, y, z) za pomocą  wyraż enia (6). Wartość tej cał ki oznaczymy przez Cx.

Do ostatniego skł adnika w pierwszej cał ce zastosujemy teraz twierdzenie G aussa-Ostrogradzkiego w postaci

C C C (

 8

~

h d

~

h d

~

h

\j- ,  C f C r (

d u (10)  H_ _ + a _ _ +  w \ dV+ /j  _ _ J J J \  8x By 8z j J JJ \  dx

_ +

\ dy dz

= I

Wyraż enie w nawiasie drugiej cał ki obję toś ciowej jest dywergencją  wektora prę dkoś ci. Jest zatem równe wydatkowi ź ródeł wewnę trznych Q w obszarze filtracji. Wyraż enie w nawiasie w cał ce po powierzchni S przedstawia natomiast jednostkowy przepł yw q przez brzeg obszaru filtracji. Wartość tego wydatku przez czę ść S2 brzegu jest okreś lona na podstawie warunku brzegowego (3). N a czę ś ci Si brzegu wartość rzeczywistego wydat-ku CJ nie jest natomiast znana. Pamię tać jednak należ y, że rozwią zania zagadnienia brze-gowego (1), (2) i (3) a także problemu wariacyjnego (5) i (2) poszukiwać należy wś ród funkcji h(x, y, z) speł niają cych zasadniczy dla tego problemu warunek brzegowy (2) [1].

Na brzegu St musi być zatem speł niona zależ ność

0 0 h(x> y, z) =  h(x, y, z) =  g(x, y, z) (x, y, z)eS1.

Tym samym wartość tej cał ki po brzegu Sx jest dla konkretnego zadania stał ą , nieza-leż ną  od wyboru funkcji h(x, y, z). Oznaczymy jej wartość przez  C2.

Podstawiają c do funkcjonał u (9) wyraż enia wynikają ce z twierdzenia G aussa- Ostro-gradzkiego (10) i wykorzystują c wnioski z przeprowadzonych powyż ej rozważ ań uzy-skamy

Porównują c funkcjonał y (5) i (12) widać, że róż nią  się  one tylko o stał y skł adnik, jest bowiem :

(13)

Poszukują c przybliż onego rozwią zania problemu filtracji poprzez minimalizację funkcjonał u (5) lub funkcjonał u (12) uzyskać musimy identyczne wyniki. Tym samym teza, że metody energetyczne dla problemów filtracji ustalonej traktować moż na jako procedurę  minimalizacji metodą  najmniejszych kwadratów waż onych btę dów wyznacza-nia skł adowych wektora prę dkoś ci filtrują cej cieczy został a wykazana.

W identyczny sposób interpretować moż na metody energetyczne dla innych proble-mów fizycznych opisanych równaniem (1) i warunkami brzegowymi (2) i (3). Wystarczy tylko podać fizyczną  Interpretację  funkcjonał u (8) i wyraż eń w formule (10). D la proce-sów fizycznych w których konieczne jest uwzglę dnienie warunku brzegowego trzeciego rodzaju

(14) ucosa + vcosfi + wcosy+ah =  q, (x,y,z) e S3,

gdzie iloczyn ah uwzglę dnia takie zjawiska jak na przykł ad straty konwekcyjne i radia-cyjne, funkcjonał  (8) musi zostać rozbudowany o skł adnik

(15) 2Jf da(h- h)

2

dS.

S3

Skł adnik ten moż na interpretować jako waż ony bł ą d kwadratowy wyznaczenia funkcji /; n a czę ś ci S3 brzegu.

3. Wyprowadzanie funkcjonał u bł ę du wyznaczenia skł adowych wektora prę dkoś ci

W podobny sposób moż na, jak się  wydaje, zinterpretować metody energetyczne dla innych zagadnień brzegowych. Wystarczy jedynie odgadną ć postać takiego funkcjonał u metody najmniejszych kwadratów, którego minimalizacja jest równoznaczna z minima-lizacją  funkcjonał u energii dla danego zagadnienia brzegowego i nadać mu odpowiednią interpretację  fizyczną . Tą  drogą  w pracy [2] zinterpretowano na przykł ad metodę  elemen-tów skoń czonych w zagadnieniach teorii konstrukcji jako procedurę  minimalizacji wa-ż onych bł ę dów wyznaczenia skł adowych dewiatora i aksjatora naprę  minimalizacji wa-ż eń.

Jednakże w pewnych przypadkach odpowiednią  postać funkcjonał u bł ę du moż na, jak się  okazuje wydedukować wprost z postaci równania róż niczkowego problemu i jego interpretacji fizycznej. Rozpatrzmy równanie róż niczkowe

(16) An =  /

gdzie A jest operatorem róż niczkowym problemu  a / i u są  daną  i szukaną  funkcją  w pew- nej przestrzeni funkcyjnej. W przypadku gdy operator A jest operatorem dodatnio okre-ś lonym m oż na wprowadzić operator B taki, że A =  B2

 i dowieść [1] istnienie ograniczo-nego operatora B "1  do niego odwrotnego. W tym wypadku rozwią zanie równania (16) jest równoznaczne rozwią zaniu równania (17) Bu =  B- 1 / .

Ponadto funkcjonał

y metody energetycznej dla równania (16) i najmniejszych kwadra-tów dla równania (17) róż ni

ą  się  wtedy tylko o stał y skł

adnik [1]. Z tego wynika, że w przy-padku dodatnio okreś loneg

o operatora A zastosowanie metody energetycznej do proble-mu opisanego równaniem (16) prowadzi do tego samego wyniku co zastosowanie metody

najmniejszych kwadratów do równania (17).

Problemem pozostaje oczywiś cie konstrukcja równania (17). Nie wydaje się  bowiem

moż liw

e podanie efektywnej drogi znalezienia operatora B w każ dy

m przypadku. P

o-nadto operator B"

1

 jest przecież poszukiwanym rozwią

zaniem problemu, zachodzi bo-wiem B

- 1

Bw =  u. Tym samym konstrukcja równania (17) jest w zasadzie równoznaczna

znalezieniu poszukiwanego rozwią zania.

Tym niemniej jeż el

i uda się  rozł oż y

ć operator równania róż niczkow

e o (16) na ilo-czyn dwóch identycznych operatorów B moż na także zinterpretować prawą  stronę

 rów-nania (17). Tym samym uzyskać moż na funkcjonał  metody najmniejszych kwadratów

dla równania (17), funkcjonał  o okreś lonej interpretacji fizycznej.

Idą c tą  drogą  wykaż emy

, że wprowadzony a priori funkcjonał  (8) jest w analizowanym

problemie identyczny z funkcjonał em metody najmniejszych kwadratów dla równania

(17).

W przypadku problemów filtracji ustalonej operator A z równania (1) jest dodatnio

okreś lony i daje się  przedstawić w postaci iloczynu dwóch identycznych operatorów

macierzowych wzglę dem siebie transponowanych. Istotnie, moż

na sprawdzić bezpo-ś rednim rachunkiem, że

(18) A =  B

r

B

gdzie macierzowy operator B ma postać

(19)

B =  KD  =

- \ 'k

x

 0 0

0 _ j/ fe, 0

0 0 -  l/j

xd

Ty

d dz

Tym samym zależ noś

ć (17) dla problemu opisanego równaniem (1) moż

na zapisać na-stę pują co

,/ r

 dh

\

(20) KD / J

 =

y x

 dx

nr

 8h

- yk

• ' 8y

-  8h

Prawa strona powyż sze

j równoś ci nie jest dana w jawnej postaci. Aby jednak równość

była speł niona winna być ona równa stronie lewej. Biorą c pod uwagę  prawo D arcy'ego

(4) widzimy, że lewa strona zależ noś c

i (20) przedstawia podzielone odpowiednio przez

VK, Vk

y

, \ / k

z

 skł adowe u, v, w wektora prę dkoś ci. Moż na wię c napisać, że

(21)

dh_ dx ~dy dh dz

fk

_x V

\ / k

y

w

P o d st a wi a ją c  t e r a z w m iejsc e h wa r t o ść p r zyb liż o ną h  r ó wn a n i e (21) n ie bę d zie oczy-wiś c ie sp e ł n i o n e .  J a k o n a jle p sze p r zybliż e n ie  m o ż na u wa ż ać t o ,  k t ó r e m in im a lizu je funk-c jo n a ł   m e t o d y n a jm n ie jszyje funk-c h  k wa d r a t ó w d la t e go  r ó wn a n i a

-  2

(22)

£ =

)/ kx

w

- KD/ 1

dV.

Wykonując dział ania w funkcjonale (22) z wykorzystaniem formuł  (4) uzyskamy

u

 lT

-  dh\

2

 I v

 ,T-

 dk\

2

 . /  w

(23) E =

ih

(w —w)

w)

2

Otrzymaliś my zatem wyjś ciow

y funkcjonał  (8). Innymi sł

owy wprowadzony poprzed-nio funkcjonał  (8) jako cał ka z waż onych bł ę dów kwadratowych wyznaczenia skł

ado-wych wektora prę dkoś c

i jest identyczny z funkcjonał

em metody najmniejszych kwadra-tów dla równania (17). Funkcjonał  ten istotnie róż ni się od funkcjonał

u metody energe-tycznej dla równania (1) tylko stał ą co wykazano poprzednio.

Zastosowane tutaj postę powani

e może być, jak się wydaje uż

yte do interpretacji po-dobnych zależ noś c

i dla innych procesów fizycznych.

4. Metoda elementów skoń czonych

W ś r ód  m e t o d e n e r ge t yc z n yc h st o so wa n yc h d o r o z wią z ywa n ia z a ga d n i e ń filt racji usta-l o n e j n a jwię ksze  u z n a n i e z d o b ył a o st a t n i o  m e t o d a e racji usta-l e m e n t ó w sk o ń c z o n ych ( M E S ) [3]. M o ż na  ją o c zywiś c ie t a k ż e  t r a k t o wa ć  j a k o p r o c e d u r ę m in im a liz a c ji  b ł ę du wyzn aczan ia s k ł a d o wy c h  w e k t o r a p r ę d k o ś ci w sen sie n a jm n ie jszyc h  k wa d r a t ó w.  M a  t o o t yle waż ne z n a c z e n i e  p r a k t yc z n e , że z je d n e j st r o n y p o z wa la  n a in t e r p r e t a c ję n i e k t ó r yc h wł aś ciwoś ci r o z wi ą z ań u z ysk i wa n yc h za je j  p o m o c ą a z d r u gie j st r o n y u m o ż l i wia o gó ln e wn ioskowa-n i e o z b i e ż ioskowa-n o ś ci  m e t o d y o r a z  wa r u n k a c h i k r yt e r i a c h t ej zb ie ż n o ś c i.  D l a p r z yk ł a d u :

(i)  R o z wi ą z a n ia p r z yb l i ż o ne m in im a lizu ją ce ś r e d ni  b ł ą d  k wa d r a t o wy m ają  ja k wia-d o m o [2] t e n k wia-d e n c ję  k wia-d o o sc yla c ji  w o k ó ł  wa r t o ś ci  k wia-d o k ł a k wia-d n y c h a  z a t e m p r ę k wia-d k o ś ci obliczone m e t o d ą  e l e m e n t ó w  s k o ń c z o n y ch b ę dą t a k ż e wyk a z ywa ć  t a k ą n a t u r ę . I st o t n i e , t a k

a oscy-lacja skł adowych wektora prę dkoś ci jest powszechnie obserwowana w M E S, szczególnie w pobliżu osobliwoś ci w rozkł adzie prę dkoś ci. Z tego też powodu ja ko ostateczn e wartoś ci przyjmowane są  zwykle wartoś ci uś rednione z róż n ych elementów.

D la ilustracji tego faktu w tabeli 1 zestawiono u/ k* i vjky obliczone w wę ź le 23 dla obszaru filtracji i podział u n a elementy przedstawionego na rys. 2. Obliczenia wyko n an o za pomocą  program u zreferowanego w pracy [4].

Tabela 1 Prę dkoś ci obliczone z elementu 27 28 29 36 37 38 Wartoś ci ś rednie Wartoś ci dokł adne

;( kx - 0 , 5 2 7 - 0 , 3 3 8 - 0 , 3 5 2 - 0 , 5 2 7 - 0 , 5 3 2 - 0 , 3 5 2 - 0 , 4 3 8 - 0 , 4 4 3 V kf - 0, 223 - 0, 412 - 0, 412 - 0,587 - 0,587 - 0,767 - 0,498 - 0, 483

Rys. 2. Obszar Filtracji, warunki brzegowe i podział  na elementy dla zadania przykł adowego (a) oraz nu-meracja wę złów i elementów W otoczeniu wę zł a 23.

(ii) Bł ą d bezwzglę dny wielkoś ci uzyskanych przez minimalizację  bł ę du kwadratowego może być tego samego rzę du dla wszelkich wartoś ci (duż ych i mał ych) obliczonych tym sposobem.

Tym samym należy liczyć się  z wię kszymi bł ę dami wzglę dnymi wyznaczenia skł ado-wych wektora prę dkoś ci dla obszarów gdzie prę dkoś ci t e są  mał e. W tabeli 2 zestawion o przykł adowo wartoś ci u/ kx i v/ ky obliczone dla wę zł ów leż ą cych wewną trz obszaru fil-tracji z rys. 2 wzdł uż prostej x =  0,25. Obliczone bł ę dy wzglę dne wskazują  istotn ie n aj-wię ksze wartoś ci dla prę dkoś ci mał ych co do moduł u. Jest to szczególnie widoczn e w tym przypadku dla skł adowej v.

Tabela 2 Wę zeł 8 13 18 23 28 33 38 M E S - 0, 074 - 0, 165 - 0, 280 - 0, 438 - 0, 663 - 0, 988 - 1, 464 H k_x D okł adnie - 0,077 - 0, 167 - 0, 283 - 0, 443 - 0, 672 - 1, 006 - 1, 497 Błą d wzglę dny w % 3,9 1,2 1.1 1,1 1,3 1,8 2,2 M E S - 0, 215 - 0, 264 - 0, 354 - 0, 498 - 0, 717 - 1, 046 - 1, 534 V k_y D okł adnie - 0, 207 - 0, 255 - 0, 342 - 0, 483 - 0, 699 - 1, 024 - 1, 509 Błą d wzglę dny w % 3,9 3,5 3,5 3,1 2,6 2,1 1,7

(iii) Zbież noś

ć metody elementów skoń czonych

, traktowanej jako procedura mini-malizacji bł ę du kwadratowego wynika z poniż szeg

o twierdzenia [5]:

Jeż eli

: 1) cią g elementów BN,, (n =  1, 2, ...) jest zupeł ny w rozpatrywanej przestrzeni

funkcyjnej, do której należą  zarówno poszukiwania funkcja u jak i znana funkcja/

z rów-nania (16) 2) równanie (17) ma rozwią

zanie oraz 3) istnieje ograniczony operator od-wrotny B "

1

 wówczas ukł ad równań uzyskiwany metodą  najmniejszych kwadratów ma

jedno i tylko jedno rozwią zanie, które przy n - > oo zapewnia zbież noś

ć u„ do u jak i Bu„

d o B-

1

/ .

Zupeł ność elementów ~BN

n

 oznacza tutaj fakt, że funkcje bazowe w MES muszą  być

tak dobrane aby ze wzrostem n dowolnie dokł adnie przybliż y

ć rzeczywisty rozkł ad prę

d-koś ci. N atomiast zbież noś

ć u„ ~> u oraz B«„ do B~

1

/  oznacza w przypadku problemów

filtracji ustalonej zbież noś

ć zarówno funkcji h (x, y, z), opisują cej rozkł ad ciś

nienia piezo-metrycznego, jak i skł adowych wektora prę dkoś ci filtrują cej cieczy do wartoś

ci dokład-nych. Pamię tać jedynie należ y

, że poję cie zbież noś c

i odnosi się  tutaj do zbież noś c

i według

normy w analizowanej przestrzeni funkcyjnej.

Literatura cytowana w tekś cie

1. C . T . M H XJI H H ; BapuaifuoHnue Memodu e MarneMamutecKoii c/ iuMKe, H3fl. H ayua, MocKBa, 1970. 2. L. R. H ERRMAN N ; Interpretation of finite element procedure as stress error minimization procedure,

J. Eng. M ech. D iv. 98, 1972, 1330- 1336.

3. O. C. ZlENKiEwrcz; Metoda elementów skoń czonych, ARKAD Y, Warszawa 1972.

• 4. B. WOSIEWICZ, S. FRĄ CKOWIAK; Program obliczeń pł askiej filtracji ustalonej na EMC ODRA 1204, Wiad. M el. i Łą k., 1/ 1979, s. 25- 27.

- 5. S. G . M ICH LIN , C. L. SMOLICKI; Metody przybliż one rozwią zywania równań róż niczkowych i cał

P e 3 w M c

K H H TEPnPETAlJiH H  3H EPrETH ^ECKH X METOJIOB PEIHEHKLS

3ARAVL  CTAUHOHAPHOfl 4>HJIbTPAUHH

B KaiecTBe MepBi OU IH SKH npn6nH > Kennoro pemenH H  3ap,a^ui cpujibTpannu B Bapuaiwoiuioft ( p o p -jwyjmpoBKe (5) H  (2) H JIH 3i<BHBajieHTH0fi eft KpaeBoii 3afla^ra (1), (2) H  (3) n p n n srro B3BeiueHHyio KBa-flpanWH yio HeBJi3i<y onpeflejienH H  cocTaBJiaiomnx seKTopa CKopocTH  cJMJibTpyiomeHCH  HCHflKOCTH  ( 7) .

J(0Ka3aH0, MTO MHHHMH3aqHH TaKOH  HeBH3KH  n o BCeS oSjiaCTH  dpHJlbTpamiH  (8) npHBOflHT I< MHHH-MH38UHH (JiyHKmioHajia flHCCHnaUHH 3H eprnH  (5).

flaH H bift noAxofl K OHepreTHnecKHM MeTOflaM no3BanJieT onpeflejiHTŁ Hei- coTopwe CBoficTBa pemeH H it, aTaiffl<e, CXOAHMOCTB caMoro MeTo«a.  O H npoaH ajiH aiipoBaH  fljia iwerofla KOHeHHtix

S u m m a r y

AN  IN TERPRETATION  OF VARIATION AL METH OD S FOR STEAD Y SEEPAG E PROBLEM S Analysed boundary value problem for steady seepage is described by partial differential equation (1) with boundary conditions (2) and (3) or by energy functional (5) with boundary conditions (2). The weigh-ted squared error of velocity components (7) is assumed as the error measure of an approximate solution to the problem. It is proved that the minimization of the weighted squared error over whole flow domain (8) is equivalent to the corresponding variational formulation (5). Thus, variational methods for steady seepage problems may be interpreted as the least squared error procedure where the measure of error is the weighted function of the errors in the velocity components. The formulation of variational methods as a minimum least squared error procedure is susceptible to interpretation of some properties of the ob-tained solutions and affords a method of convergence. These aspects have been discussed for the finite elements method.

IN STYTU T BU D OWN ICTWA WOD N O- M E LI OR AC YJN E G O POZN AŃ

Praca został a zł oż ona w Redakcji dnia 13 listopada 1979 roku.