1 Metoda elementów brzegowychLu=f plus warunki brzegowe możliwe podejścia: 1) Metoda różnic skończonych 2) Metoda elementów skończonych obszar podzielony na elementy, w każdym z elementów rozwiązanie u w postaci wielomianu interpolacyjnego silna postać równania + ilorazy różnicowe 3) Metoda elementów brzegowych tylko brzeg dzielony na elementy w których zadajemy u a liczymy du/dn (lub na odwrót)

2 resztę łatwo oszacować:Formalizm ważonych reszt: (uniwersalny dla metod przybliżonych) u = rozwiązanie dokładne (r) rozwiązanie przybliżone jego błąd (reszta) resztę łatwo oszacować: chcielibyśmy e(r)=0, ale zadowolimy się: w(r) =fcja wagowa słaba forma równania różniczkowego obustronnie mnożymy przez funkcję wagową, całkujemy obustronnie i żądamy znikania błędy w sensie „słabym” (błąd ortogonalny do funkcji wagowej )

3 1. Metoda różnic skończonych w formalizmie reszt ważonych ( odpowiedność dla słabej formy równania) 1) przestrzeń dzielimy na identyczne elementy, 2) w każdym elemencie rozwiązanie przybliżone rozwijamy w bazie funkcji kształtu (interpolacja Lagrange’a na węzłach narożnych) 3. Funkcje wagowe: (tyle funkcji ile węzłów) każda waga: delty Diraca scentrowana w jednym z węzłów. 1,2,3+wartość główna pochodnej na granicy elementów = równania MRS

4 Metoda różnic skończonych:continuum opisane na dyskretnej siatce + O(dx2,dy2) dx metoda najprostsza w zastosowaniu (duży komfort jeśli rachunek pójdzie tą metodą), lecz: podaje rozwiązanie tylko w dyskretnych punktach wolno zbieżna: wymaga wielkich siatek (często zbyt wielkich) słabo sprawdza się lokalne zagęszczanie siatki (błąd ilorazu różnicowego skacze o jeden rząd przy nierównej siatce) dy

5 2. Metoda elementów skończonych1) elementy np. czworokątne, rozłożone jak nam pasuje 2) w każdym elemencie rozwiązanie przybliżone rozwijamy w bazie funkcji kształtu 3. Funkcje wagowe: fk (metoda Galerkina): równania na ck Błąd ortogonalny do bazy. Metoda równoważna wariacyjnej.

6 3. Metoda elementów brzegowych1) funkcja wagowa = rozwiązanie fundamentalne (swobodna funkcja Greena związana z operatorem L ) = pozostanie nam równanie całkowe wprowadzimy elementy brzegowe, oraz funkcje kształtu, które opiszą wartość u oraz jej pochodnej normalnej policzymy „drugi warunek brzegowy” wartość u w dowolnym punkcie dana będzie przez całkę konturową wzdłuż brzegu

7 Lu=f plus wb. niejednorodne równanie różniczkowe:rozwiązanie fundamentalne (swobodna funkcja Greena) niejednorodne równanie różniczkowe: Lu=f plus wb. (*) rozwiązanie fundamentalne x=zmienna niezależna, a = źródło punktowego zaburzenia na rozwiązanie fundamentalne* nie nakłada się żadnych warunków brzegowych, poza (*swobodną funkcję Greena – w przeciwieństwie do zwykłej funkcji Greena)

8  rozw. fundamentalne będzie określone z dokładnościądo rozwiązania równania jednorodnego Lu=0 funkcja Greena: jedno z rozwiązań fundamentalnych: ma spełniać zadane warunki brzegowe w MEB swobodna funkcja Greena: dowolna, ogólna byle dana prostym wzorem, warunków brzegowych spełniać nie musi (metoda ma działać dla dowolnych warunków brzegowych)

9 = H(xi-x) dalej pracujemy z L=2 rozwiązanie fundamentalne 1D: y x=xiu*(x;xi) x Ciągła funkcja nieciągła pochodna

10 3D 2D 3D 2D

11 rozwiązanie fundamentalne 2D:współrzędne cylindryczne względem źródła ze względu na symetrię poza źródłem: równanie Laplace’a z radialną symetrią poza źródłem równanie fundamentalne spełnione, czy również w r=0 ?

12 czy D ri=0 całkujemy obustronnie, ma wyjść: całka powierzchniowa po D:

13 ri zamiast po D łatwiej po kole otaczającym źródło(można, bo poza nim: ) D n całkujemy obustronnie ri całka powierzchniowa po D: + całkowanie przez części 2D a=1, b=

14 e e n ri =0 =-1 Dla A=-1/2p: laplasjan z spełnia definicje delty Diraca B pozostaje dowolne, zazwyczaj bierze się B=0

15 2D: Podobnie można sprawdzić, że w 3D:

16 niejednorodne równanie różniczkowe rozwiązane przy pomocy swobodnejfunkcji Greena, przykład + wb u(0)=u(1)=0 rozwiązanie: u(x)=x(1-x2)/6 metoda ważonych reszt całkujemy przez części [u’w]’=u’’w+u’w’  u’’w = [u’w]’-u’w’ drugi raz przez części [uw’]’=uw’’+u’w’  u’w’= [uw’]’-uw’’

17 weźmy: ! wymagana znajomość zarówno u jak i jego pochodnej (normalnej) na brzegu !

18 1D: W W 2D: G G odpowiednik równania całkowego,które będziemy rozwiązywać w 2D odpowiednik całki objętościowej odpowiednik całki powierzchniowej w 2D dostaniemy 1D: W pochodna normalna: + d / dx z prawej części przedziału - d / dx z lewej części przedziału W 2D: G G

19 xi z przedziału (0,1), własności delty + warunki brzegowe na u(0)=u(1)=0:można rozwiązać analitycznie

20 ... analitycznie: xi z przedziału (0,1), x=0 x=1 x=1 x=1 x=1

21 bc : jeszcze raz u(1)=0 uwaga: wyznaczyliśmy po drodze w. Neumannax=1 x=1 x=1 bc : jeszcze raz u(1)=0 uwaga: wyznaczyliśmy po drodze w. Neumanna (i tak będziemy robić w MEB)

22 Rozwiązanie fundamentalne a funkcja Greena[z warunkami brzegowymi] + wb u(0)=u(1)=0 plus g(x=0)=g(x=1)=0 doprowadziło nas do równania całkowego: fcja Greena: konstruowana pod równanie niejednorodne odpowiedź na źródło punktowe+ warunki brzegowe narzucone na równanie niejednorodne, które chcemy rozwiązać

23 Funkcja Greena z warunkami brzegowymiplus WB np. g(x=0)=g(x=1)=0 własności g? 1) spełnia określone warunki brzegowe na brzegu pochodna w źródle: 2) skok pochodnej w źródle (mieliśmy tak również w rozw. fund.) 3) ma spełniać równanie jednorodne poza źródłem

24 Axi =B(xi-1) B = - xi A = (1-xi) B-A=-14) zażądajmy ciągłości fcji Greena w xi (różnić się będzie od swobodnej funkcji Greena o rozwiązanie równania Laplace’a, a to – ciągłe) Axi =B(xi-1) plus skok pochodnej: B = - xi A = (1-xi) B-A=-1 xi 1

25 funkcja Greena (z warunkami brzegowymi)swobodna funkcja Greena z u* do g: dodane rozwiązanie Laplace’a -xix y x=xi u*(x;xi) xi g x 1 własność ogólna: g(a;b)=g(b;a) [wpływ źródła w a na punkt b taki jak źródła w b na punkt a]

26 - Jaki zysk z g? + wb u(0)=u(1)=0 weźmy: w=g(x;xi)plus WB np. g(x=0)=g(x=1)=0 zamiast: ogólnie Lu=f Rozwiązanie równania różniczkowego: Niejednorodność scałkowana z funkcją Green’a -

27 znalezienie analitycznej formy Green z warunkiem brzegowym Greena swobodna znacznie prostsze dlaczego w BEM stosuje się trudniejszy przepis? znalezienie analitycznej formy funkcji Greena g z dowolnymi warunkami brzegowymi dla dowolnie skomplikowanego brzegu jest trudne Dlatego zamiast g: swobodna funkcja Greena: prosty wzór, zależny tylko od typu równania

28 Równanie różniczkowe do całkowegosłaba forma równania z fcją Greena jako wagą G n

29 BEM odpowiada metodzie reszt ważonych z wagą w postaci rozwiązania fundamentalnego:Całkowanie przez części TGG:

30 - z oznaczeniami: -

31 - pochodzi z: Podstawowe równanie MEB:reszty ważone z wagą w postaci rozwiązania fundamentalnego druga tożsamość Greena aby przerzucić operator różniczkowy z nieznanego rozwiązania do rozwiązania fundametalnego własności delty Diraca Zamieniliśmy cząstkowe równanie różniczkowe na równanie całkowe (w sposób ścisły): do tej pory : brak przybliżeń

32 Jakie ci ? ci=1 ci=0 ci zależy : gdzie umieścimy źródło ri ri wewnątrzna zewnątrz ci=1 ri ci=0 W W

33 W na brzegu ri G dla metody elementów brzegowych: akurat sytuacja z brzegiem jest najważniejsza ri na brzegu: ci=? brzeg ma prawo mieć kanty na gładkiej części brzegu: w powiększeniu fragment gładkiej części brzegu jest odcinkiem prostej (delta Diraca przecięta na pół) ri w ci=½ wnętrze

34 ci=w / 2p 1D: delta jest „funkcją” parzystącałkujemy po prostych przechodzących przez ri i obracamy je o kąt p: kręcimy całka po kierunku niebieskim: zero [zawsze poza obszarem W] ri całka po kierunku zielonym: ½ u(ri) w w kąt wewnętrzny przyczynek od „kierunków zielonych”: wnętrze: W ci=w / 2p

35 podstawowe równanie MEB- Jeśli znamy u oraz q znajdziemy rozwiązanie w każdym punkcie wewnątrz W! Ale warunki brzegowe zadaje się albo na u albo na q, skąd drugi warunek brzegowy?

36 - Metoda elementów brzegowych: G7 brzeg dzielony na elementy W G1 G G2nie zawiera niewiadomych nie wnosi wiele do idei, Powiedzmy więc że r=0 - G7 brzeg dzielony na elementy W G1 G G2 G6 na każdym elemencie brzegowym zadany warunek brzegowy, „drugi” wyliczany wartości u oraz q opisane przy pomocy interpolacji G3 G5 G4 wielomian węzłowy Lagrange’a

37 w jaki sposób wyliczany drugi warunek brzegowy?N węzłów, N niewiadomych punkt obserwacji ri ustawiany kolejno w każdym węźle brzegowym. każde wstawienie da nam jedno równanie = układ równań liniowych

38 w elemencie Gk: problem oryginalny warunek na u zadany x w węzłachwstawimy dokładny warunek brzegowy wyliczymy q(k,l)

39 Problem modelowy: stagnation flow2u=0, u=x2-y2 grad (u) = (2x,-2y)

40 w każdym elemencie: opisujemy u oraz q na brzegu w bazie funkcji kształtu, weźmy najprostsze (liniowe): x u1 u2 Uwaga: na ogół tylko w węzłach udaje się wprowadzić dokładne warunki brzegowe

41 Wstawimy warunek Dirichleta i będziemy starali się odzyskać warunek Neumanna Dirichlet dokładny przybliżony Dirichlet (wynik interpolacji dokładnego warunku brzegowego w bazie liniowych funkcji kształtu) (x,y)=(1,1) -2= u/n 2= u/n 2= u/n -2= u/n (x,y)=(-1,-1) u/n =n • gradu (u) Problem modelowy: stagnation flow 2u=0, u=x2-y2 grad (u) = (2x,-2y)

42 macierze wpływu F1 F2 pierwszy (lewy) węzełelementu n-tego : on numerze globalnym lg(n,1) [buchalteria węzłów i składanie macierzy dokładnie tak jak w MES]

43 macierze wpływu numeruje lokalizację delty D. (źródła, fcji Greena)numeruje element numeruje fcje kształtu (lokalna numeracja węzłów) podobnie: wtedy:

44 co można przedstawić w postaci Hu=Gqmacierze H i G trzeba poskładać [jak w MES] Węzeł 1 = lokalny 1 elementu 1 oraz lokalny 2 elementu 12

45 Składanie macierzy wpływuzaczepmy źródło w i=1 [pierwszy wiersz] Hu=Gq pierwszy wiersz macierzy H: u1 u2 ... = q1 q2 ... w j-tym wierszu c wchodzi do diagonali [do kolumny j] pierwszy wiersz macierzy G:

46 liczymy lokalne macierze wpływucałkować będziemy po elemencie odniesienia ze zmienną x przestrzeń n wnętrze

47 xa=xi całkować będziemy po elemencie odniesienia ze zmienną xmożliwa próba scałkowana 1/x od 0 xa=xi

48 Zagrożenie jest pozorne: na prostym elemencie, całkować będziemy po elemencie odniesienia ze zmienną x n Zagrożenie jest pozorne: na prostym elemencie, w którego węźle zlokalizowana jest funkcja Greena jej gradient jest prostopadły do wersora normalnego do powierzchni.

49 jeśli element nie będzie prosty: grad u będzie lokalnie prostopadły do n , co wystarczy aby całka skończona w praktyce: ta całka jest najtrudniejsza do obliczenia (najwolniej zbieżna) [można uniknąć jej liczenia korzystając z reguły sumowania podanej poniżej]

50 Hii=ci W problemie laboratoryjnym całki macierzy wpływu hsą łatwe do policzenia dajmy osobliwość do węzła i=2 (składamy drugi wiersz układu równań) na tym fragmencie brzegu: elementy h są zero bo gradient z rozwiązania fundamentalnego jest do brzegu styczny Na pozostałej części brzegu elementy h nie znikają, ale są łatwe do policzenia bo 1/(r-r2) zbytnio się nie zmieniać nie będzie na długości jednego elementu brzegowego Hii=ci wniosek:

51 odpowiednia chwila by zobaczyć, że: macierze wpływu MEB są niesymetryczne, [w MES dla równania Poissona macierze sztywności były symetryczne] pierwszy wiersz H ze źródła wstawionego do węzła 1 zero (węzeł 2 należy do elementu 1 i 2) drugi wiersz ze źródła wstawionego do węzła 2 węzeł 2 nie należy do elementu 12 więc zero nie wyjdzie ten wyraz jest zero węzeł 2 sąsiedni z 1

52 kwadratury Gaussa świetnie się sprawdzajądruga lokalna macierz wpływu: gdy osobliwość znajduje się w jednym z węzłów: trzeba ją scałkować [ta osobliwość całkowalną jest] ln(x) kwadratury Gaussa świetnie się sprawdzają

53 przez Gaussa 32 punktowego dla a=1Całkowanie kwadraturą Gaussa funkcji z logarytmiczną osobliwością punkty wezwane przez Gaussa 32 punktowego dla a=1 wynik num: zamiast: dla a=1

54 niebieski: Gauss numeryczny 32 punktyczerwony: analityczny czarny: błąd [analityczny-numeryczny]

55 Składanie globalnych macierzy wpływuLokalizacja osobliwości rozwiązania fundamentalnego w węźle i-tym produkuje i-ty wiersz macierzy H oraz G zaczynamy od Hii=ci Liczymy macierze lokalne, a następnie składamy globalne: lg(k,l) – globalny numer węzła l-tego w k-tym elemencie brzegowym

56 wyniki dla du/dn na brzegu: (wstawiamy D, liczymy N) 2.22 -2.22 Po złożeniu mamy: Hu=Gq traktujemy u jako dane, q jako szukane dokładne: -2= u/n 2= u/n 2= u/n -2= u/n

57 wynik: zerowa u/n na kantachMEB: przewiduje jedną wartość pochodnej normalnej w węźle jaka ma być pochodna normalna na kancie? jeśli będziemy udawać, że kant jest zaokrąglony -2 2 grad (u) = (-2x,-2y) zero akurat jest wynikiem poprawnym dla takiego węzła.... ale problem jest kłopotliwy numerycznie

58 4 elementy brzegowe na bok kwadratu (równe odstępy):u/n plus plus minus

59 9 elementów brzegowych na bok kwadratu (równe odstępy):0 na kancie ma uzasadnienie, jako wartość główna pochodnej ale powoduje problemy dla interpolacji rozwiązania na pionowych i poziomych końcach. metoda wolałaby dwa parametry niezależne. górny brzeg: u/n

60 pomysł na rozwiązanie problemu-2 -2 2 2 osaczyć kant nowymi wierzchołkami

61 równo rozłożone wierzchołki1 wierzchołek extra o 0.01 od narożnika u/n x dla y=+1 (wzdłuż górnego brzegu)

62 drugi wierzchołek extra o 0.01 od poprzedniegou/n x dla y=+1 (wzdłuż górnego brzegu)

63 trzeci wierzchołek extra o 0.02 od poprzedniegou/n x dla y=+1 (wzdłuż górnego brzegu) wniosek: problem można zawęzić do okolic brzegu, gdzie istotnie występuje

64 Po wyznaczeniu brakującego warunku brzegowego Wyznaczmy u wewnątrz obszaru całkowania ze wzoru ... 1 trzeba macierze wpływu przeliczyć dla każdej (nowej) lokalizacji punktu obserwacji rie Uwaga: dla danego i całek do wyliczenia jest 2 (tyle funkcji kształtu) razy 2 (h i g) razy liczba elementów brzegowych Uwaga: mamy całkowitą swobodę w wyborze siatki punktów, w których znajdziemy ui. Liczymy tylko tam, gdzie rozwiązanie nas interesuje.

65 wyniki na osi x=0 czarne= numeryczne czerwone –y2 9 równych elementów Na 1 bok kwadratu 3 równe elementy na 1 bok kwadratu u u y (dla x=0) y (dla x=0)

66 3 równe elementy na bok kwadratu wzdłuż: x=0.9 czarne= numeryczne czerwone .92 –y2 9 równych elementów na bok kwadratu (węzeł .77 potem 1) wb. D warunek Dirichleta wciąż nie całkiem spełniony [kwadratowe funkcje kształtu rozwiązałyby problem – po prawej]

67 Hu=Gq u=x2-y2 Dx=Cy Inne warunki brzegowe:Dirichlet na niezerowej części brzegu Neumann na reszcie Dirichlet w jednym punkcie u=x2-y2 Hu=Gq Składamy jak zawsze Reorganizujemy macierze przekładając wiersze z jednej do drugiej Dx=Cy układ równań: macierz D, w której pierwszy wiersz bierzemy z H, resztę z G w C: na odwrót y=(u1,q2,q3,...) x=(q1,u2,u3,...) musimy zostawić jedno u na brzegu zostawiłem u(-1,1)=0 Wyniki dla u Wzdłuż górnego brzegu Rozbieżność w wartościach u w węźle

68 Hu=Gq Hu=Gq Hu+H[c]=Gq j Reguła sumjeśli u spełnia równanie Laplace’a , u+c również, ponadto H i G bez zmian dla każdego i Hu=Gq Hu+H[c]=Gq j wektor o stałych wyrazach użyteczność: dla krzywoliniowych elementów najtrudniejsza jest policzenie diagonalnych wkładów do H - gdzie całka jest najwolniej zbieżna (dla krzywoliniowych elementów) Ale: do jej oszacowania można użyć pozostałych elementów, które szybciej zbieżne] j, co więcej: nie musimy martwić się o wyznaczenie ci

69 - druga reguła sum (równanie jednorodne)2) podobnie jeśli rozwiązanie fundamentalne to również +c macierz suma funkcji kształtu dla każdego elementu = 1 stąd po złożeniu każdy elementy macierzy G rośnie o c Hu=Gq Hu=(G+[c])q H nie ulega zmianie, bo w liczeniu elementów h wykorzystywana jest tylko pochodna normalna należy czytać: zdyskretyzowana forma prawa Gaussa strumień zero, jeśli źródeł brak użyteczność drugiej reguły sum w powyższej formie ograniczona do równania jednorodnego -

70 Wersja metody elementów brzegowych dla równania niejednorodnego:- - całka objętościowa Hu=Gq+n -

71 Całkowanie niejednorodności z rozwiązaniem fundamentalnym (osobliwością):- ri przy rozwiązywaniu równań dla „drugiego warunku brzegowego” nie ma problemu, gdy gęstość ładunku jest skończona (całkujemy wtedy logarytm * skończona funkcja – jak elementy macierzy g) r(r) problem z osobliwością rozwiązania fundamentalnego wewnątrz W pojawia się gdy poszukujemy rozwiązania wewnątrz obszaru całkowania [współrzędne kartezjańskie słabo się do całkowania wtedy nadają]

72 ... gdy poszukujemy rozwiązaniawewnątrz obszaru całkowania W z objętości całkowania wycinamy kulę (koło) o zadanym promieniu i środku w osobliwości ri całkujemy we współrzędnych biegunowych względem osobliwości całka po kole: jakobian likwiduje osobliwość

73 r ln (r) ln (r) jakobian likwiduje osobliwość w źródle podobnie 3D:r2 1/r jakobian fcja Greena całka po kole: gdy źródło na brzegu:całkowanie po kącie węższe, lub R(f)

74 MES porównanie MEB Elementy generowane tylko na brzegu G(wymiar o jeden niżej). Ale: niejednorodność należy wycałkować po całym obszarze jej występowania. [co może być łatwe, jeśli np. źródła punktowe] Elementy generowane na całym obszarze całkowania W równania różniczkowego W G 1) Najpierw rozwiązujemy równanie na „drugi wb” na brzegu G Otrzymujemy rozwiązanie na całym W 2) Następnie liczymy rozwiązanie tylko w punktach wnętrza które nas interesują Rozwiązywane jest równanie przybliżone (słaba forma równania z wagami = funkcjami bazowymi). Jeśli dyskretyzacja brzegu dokładna: również warunek brzegowy wstawimy dokładny (ale w ograniczonej bazie). Rozwiązywane jest równanie dokładne (słaba forma równania z funkcją Greena jako wagą = jest dokładnym odpowiednikiem oryginalnego równania). Ale: warunki brzegowe są wstawiane w sposób przybliżony.

75 MES MEB generowane są macierze duże rzadkie i symetryczne(operator samosprzężony) macierze gęste i niesymetryczne, ale mniejsze o tym, która metoda produkuje łatwiejsze w eksploatacji macierze decyduje stosunek objętości do powierzchni, gdy duży MEB preferowana. MEB stosowalna dla nieskończonych obszarów elementy macierzowe całkowalne łatwo: (macierz sztywności = analitycznie lub dokładnie numerycznie [całki z wielomianu]) wektor sztywności: nieco trudniej funkcje podcałkowe są osobliwe całkowanie numeryczne jest trudne stosowalna do tych problemów liniowych, dla których potrafimy wskazać rozwiązanie fundamentalne stosowalność: bardzo szeroka, z problemami nieliniowymi włącznie

76 Przykład: zastosowanie swobodnej fcji Greenaw rachunkach analitycznych, problem warunków brzegowych O delcie Diraca i zupełności fcji własnych operatorów samosprzężonych + wb w x=a,b zał: fn unormowane dla dowolnej f z

77 Przykład: zastosowanie swobodnej fcji Greenaw rachunkach analitycznych, problem warunków brzegowych O delcie Diraca i zupełności fcji własnych operatorów samosprzężonych + wb w x=a,b zał: fn unormowane dla dowolnej f z w szczególności

78 Przykład: zastosowanie swobodnej fcji Greenaw rachunkach analitycznych u=0 1 u=sin(2py) u=0 0 1 u=0 interesuje nas tylko rozwiązanie wewnątrz obszaru całkowania więc ci=1 rozwiązać można separując zmienne: poszukamy rozwiązania fundamentalnego, które znika na brzegach pudła Rozwiązania równania Laplace’a znikające na brzegach:

79 u=0 1 spróbujmy skonstruować fcję Greena dla Zerowych warunków brzegowych: u=sin(2py) u=0 0 1 u=0 fcje własne operatora samosprzężonego: baza ortogonalna

80 u=0 1 u=sin(2py) u=0 (minus bo x=0 to brzeg lewy) 0 1 u=0

81 dokładny y=p/4 suma 20 wyrazów

82 dokładny y=p/4 suma 200 wyrazów

83 dokładny y=p/4 suma 200 wyrazów

84 dokładny y=p/4 suma 600 wyrazów

85 Szereg jest wolno zbieżnysuma 600 wyrazów czasem wygodniej jest pracować z kompaktową formą swobodnej funkcji Greena oraz z pełnym równaniem całkowym.