MES w 1D

w stronę metody elementów skończonych

Lu=f (na

Ω)

Bu=g (na d

Ω)

Rozwiązanie dokładne (silnej postaci równania) jest „trudne”.

szukamy rozwiązania przybliżonego w bazie funkcji

Działając operatorami L i B na rozwiązanie przybliżone dostajemy funkcje resztkowe (rezydualne) zamiast zera:

zależy nam, aby reszty r i s były jak najmniejsze

dla metody Galerkina bierzemy funkcje bazowe jako wagi: w_j=v_j

c wyznaczamy z ważenia reszty

:

(rozwiązanie w podprzestrzeni wektorowej rozpiętej przez wektory bazy)

Silna forma równania:

Lu=f (równość funkcji w każdym punkcie

obszaru całkowania)

ważone reszty:

słaba forma równania,

Lu=f

chcemy znaleźć taki element przestrzeni żeby:

słaba forma równania

residuum (błąd) r=Lu-f: ortogonalny do każdego wektora bazowego

(r,v

)=0

residuum znajduje się poza przestrzenią generowaną przez wybraną bazę

metoda Galerkina: residuum a przestrzeń wektorowa

rozpięta przez wektory wybranej bazy

metoda Galerkina rzutuje rozwiązanie dokładne na wektory bazy

ilustracja: ue to rozwiązanie dokładne (przekątna sześcianu),

u to rozwiązanie przybliżone R tutaj to ue-u

od (a) do (c) dodajemy elementy bazowe φ1,φ2,φ3.

błąd metody: residuum – jest ortogonalne do podprzestrzeni wyznaczonej przez wektory bazy

metoda Galerkina jest zbieżna: gdy baza obejmie całą przestrzeń – nie ma miejsca na residuum

SY=F stiffness matrix macierz sztywności load vector wektor obciążeń

S

=(Lv

,v

)

F

=(f,v

)

-1.00 -0.50 0.00 0.50 1.00 -0.20 -0.10 0.00 0.10 0.20 n=2 n=4 n=6 rozwiązanie -1.00 -0.50 0.00 0.50 1.00 -0.02 -0.01 0.00 0.01 0.02

błąd ε (nie residuum tylko różnica dokładne – Galerkina):

baza używana poprzednio: SY=F

powyższy przykład: baza wielomianów określonych na całym pudle obliczeniowym. Z wielu powodów jest to zły pomysł.

Wysokie potęgi wielomianów niewygodne w użyciu: całkowanie, efekt Rungego, powód najważniejszy:

macierz S byłaby gęsta, problem nie do rozwiązania przy dużym N.

najprostszy wybór funkcji kształtu(*): baza funkcji odcinkami liniowych zbieżność dostaniemy w przestrzeni funkcji odcinkami liniowych

Galerkin z bazą odcinkami wielomianowych funkcji zdefiniowanych w sposób rozłączny przestrzennie→metoda elementów skończonych

SY=F

Metoda elementów skończonych: funkcje rozłączne tak, żeby S = rzadka

Zobaczymy w działaniu metodę elementów skończonych, ale na razie: bez jej charakterystycznych narzędzi:

bez lokalnych macierzy sztywności związanych z każdym elementem bez ich składania do macierzy globalnej

bez mapowania przestrzeni fizycznej do przestrzeni referencyjnej

będziemy mówili o metodzie z punktu widzenia węzłów:

tak najłatwiej wprowadzić metodę, ale dla 2D i 3D takie podejście okazuje się niepraktyczne

x x_i x_i+1 x_i-1 v_i(x) 1 element K_idługości h_i=x_i-x_i-1

element K_i+1długości

h_i+1 = x_i+1-x_i

węzły

funkcje bazowe i brzeg

Dla (jednorodnych) warunków Dirichleta mamy

y_pierwsze=y_ostatnie =0

W każdym elemencie: mamy 2 funkcje,

każda z innym węzłem związana

fcja kształtu

niezerowe tylko dla

i=j, i=j-1 oraz i=j+1 [bez przekrywania całka znika]

-niech j = i+1

v_i v_i+1

i-1 i i+1 i+2 gdy jedna pochodna

dodatnia druga ujemna

długość elementu o numerze większym z dwóch indeksów S

-Macierz sztywności dla n węzłów wiersz n-1 SY=F F_i=(v_i,f) + warunek y₁=y_n=0 po elemencie K_i po K_i+1 1 2 3 4 h₂ h₃ h₄

Macierz sztywności dla n węzłów wiersz n-1 SY=F F_i=(v_i,f) + warunek y₁=y_n=0 po elemencie K_i po K_i+1

dla równoodległych węzłów S jak macierz metody RS (razy h=dx), ale wektor obciążeń F – nie! w MRS mielibyśmy F_i=f(x_i) dx

1 2 3 4

dla f(x) = - sin (πx)

warunki brzegowe (jednorodne Dirichleta): forma S oraz F₁=F_n=0

ten URL wygląda prawie jak dla MRS...

zobaczmy wyniki

Układ równań z macierzą trójprzekątniową – przypomnienie. Jak rozwiązac?

SY=F

S=LU (LU – trójkątne) (LU)Y=F

UY=x

Lx=F - najpierw rozwiązujemy ten układ

Dekompozycja LU mecierzy trójprzekątniowej

dwuprzekątniowe

bez zmian

-1.00 -0.50 0.00 0.50 1.00 -0.20 -0.10 0.00 0.10 0.20

Wynik: równoodległe węzły

-1.00 -0.50 0.00 0.50 1.00 -0.01 -0.01 0.00 0.01 0.01 Błąd: MES

dokładny

-1.00 -0.50 0.00 0.50 1.00 1.50 -0.20 -0.10 0.00 0.10 0.20

MES (równoodległe węzły) a MRS (węzły w tych samych punktach):

MES dla laplasjanu

bez pochodnej z funkcjami liniowymi: w węzłach wynik dokładny !!! -1.00 -0.50 0.00 0.50 1.00 -0.01 -0.01 0.00 0.01 0.01 błąd MRS i MES

znikanie błędu MES (1D, liniowe f.kształtu) w węzłach zachodzi również dla nierównomiernego rozkładu węzłów:

-1.00 -0.50 0.00 0.50 1.00 -0.10 -0.05 0.00 0.05 0.10 -1.00 -0.50 0.00 0.50 1.00 -0.01 -0.01 0.00 0.01 0.01 błąd:

Dla MRS: dla nierównomiernej siatki musielibyśmy używać

niesymetrycznych ilorazów o

[jak widzieliśmy] niższej dokładności

-1.00 -0.50 0.00 0.50 1.00 -0.01 -0.01 0.00 0.01 0.01 błąd:

Równanie Poissona,

funkcje kształtu liniowe

wynik MES dokładny

w węzłach

MES: produkuje oszacowanie

wyniku również między węzłami

MRS: tylko w węzłach

MRS: wartości w węzłach,

są dokładne

TYLKO

w granicy Δx→0

laboratorium -1.0 -0.5 0.0 0.5 1.0 0.00 0.04 0.08 0.12 0.00 0.20 0.40 0.60 0.80 1.00

rozwiązanie (bardzo) dokładne

MRS: gęsta siatka

ρ

-1.00 -0.50 0.00 0.50 1.00 0.00 0.04 0.08 0.12 0.00 0.20 0.40 0.60 0.80 1.00

ρ

u

rozwiązanie

metodą elementów skończonych

z liniowymi funkcjami kształtu

na dziewięciu węzłach

rozwiązanie MES w tej wersji (liniowe fcje kształtu 1D) jest dokładne w węzłach

-1.00 -0.50 0.00 0.50 1.00 0.00 0.04 0.08 0.12 0.16 0.00 0.20 0.40 0.60 0.80 1.00 dla porównania: MRS na 9 – ciu węzłach MES

uwaga – na rysunku - dla MRS

punkty łączymy liniami tylko dla ilustracji dla MES z liniowymi funkcjami kształtu połączenie punktów: ma znaczenie dosłowne

-1.00 -0.50 0.00 0.50 1.00 0.00 0.04 0.08 0.12 0.00 0.20 0.40 0.60 0.80 1.00

ρ

u

między wynikiem MES a dokładnym

-1.00 -0.50 0.00 0.50 1.00 0.00 0.04 0.08 0.12 0.00 0.20 0.40 0.60 0.80 1.00

ρ

u

rozwiązanie MES w naszej bazie jest odcinkami liniowe

a rozwiązanie dokładne jest liniowe tam gdzie ρ=0 wniosek: tam gdzie ρ=0 wystarczy jeden element!

pomysł: przesunąć wszystkie węzły poza brzegowymi do obszaru gdzie nie znika gęstość ładunku – tam gdzie u zaokrąglone.

Kryterium wyboru węzłów? (bx)

przy okazji dyskusji metod relaksacyjnych dowiedzieliśmy się, że najbliższe prawdzie jest rozwiązanie, które minimalizuje funkcjonał całki działania

wykorzystajmy działanie jako kryterium jakości rozwiązania w metodzie elementów skończonych

i=2,8 x1= -x9 = -1 zacieśniamy węzły wokół x=0

c1=c9=0 (warunki brzegowe)

do oceny jakości wyboru węzłów użyjemy macierzy A i F, które i tak musimy wyznaczyć aby wyliczyć c.

0.00 0.25 0.50 0.75 1.00 bx -0.0120 -0.0118 -0.0116 -0.0114 -0.0112

a

działanie w MES w funkcji parametru bx

numer iteracji MRS

-1.00 -0.50 0.00 0.50 1.00 0.00 0.04 0.08 0.12 0.00 0.20 0.40 0.60 0.80 1.00 rozwiązanie dla optymalnego bx: czerwone: MES

niebieskie: wynik dokładny

Wybór węzłów: przez optymalizację funkcjonału ...

metoda elementów skończonych ma charakter wariacyjny

w ogóle: metoda Galerkina dla dowolnej bazy jest równoważna metodzie Reyleigha-Ritza gdy ta stosowalna

metoda Reyleigha-Ritza:

rozwiązanie w bazie funkcyjnej

dla ustalonych funkcji bazowych (w naszym przykładzie: dla ustalonych węzłów)

c wyznaczone przez warunek minimum a.

(pokazać, że warunek min a produkuje Ac=F)

Wybór węzłów: przez optymalizację funkcjonału ...

metoda elementów skończonych ma charakter wariacyjny

w ogóle: metoda Galerkina dla dowolnej bazy jest równoważna metodzie Reyleigha-Ritza gdy ta stosowalna

metoda Reyleigha-Ritza:

rozwiązanie w bazie funkcyjnej

dla ustalonych funkcji bazowych (w naszym przykładzie: dla ustalonych węzłów)

c wyznaczone przez warunek minimum a.

(pokazać, że warunek min a produkuje Ac=F)

uwaga: w naszym przykładzie : dodatkowo optymalizowaliśmy funkcje bazowe

(położenie węzłów). Zasada najmniejszego działania wykorzystana została więc dwukrotnie.

jeśli tylko znamy funkcjonał dla równania różniczkowego: przyda się do optymalizacji kształtu elementów (2D)

-1.00 -0.50 0.00 0.50 1.00 -0.01 -0.01 0.00 0.01 0.01 błąd:

Równanie Poissona,

funkcje kształtu liniowe

wynik MES dokładny

w węzłach

w MRS wszystko co możemy

otrzymać, to wartości w węzłach,

które są dokładne

TYLKO

w granicy Δx→0 !!!

Wniosek:

cały rachunek

na 2 elementach:

-1

q

+1

f(q)

niezależnie od wyboru

q: f(q) da dokładne

rozwiązanie równania

WRACAMY

-1.00 -0.50 0.00 0.50 1.00 1.50 -0.20 -0.10 0.00 0.10 0.20

q = 0.1

q=-0.6

q możemy ustawić gdziekolwiek, zawsze dostaniemy

rozwiązanie dokładne

wniosek:

wystarczy przeskanować

q przez pudło

-1

q

+1

f(q)

jeden węzeł wewnątrz pudłą

funkcja kształtu

v

= (x+1)/(q+1) dla x<q

v

= (1-x) / (1-q) dla x>q

rozwiązanie przybliżone: v=f(q)v

(x)

Lu=

ρ

F(v)=½(v,Lv)-(

ρ

,v)

funkcjonał:

wyliczymy F jako funkcję f(q)

v

= (x+1)/(q+1) dla x<q

v

= (1-x) / (1-q) dla x>q

v=f(q)v

(x)

=0

=0

policzmy pochodne f(q) po q

spełnia silną formę równania

u’’= -ρ

silna a słaba forma równania: różnica

u(x)’’= -ρ(x)

-1

q

+1

v(x

)

v’(x)

v’’(x)

funkcja v nie spełnia silnej

formy równania różniczkowego,

tylko słabą:

(v,v’’)f(q)=-(v,

ρ

)

druga pochodna potencjału = ładunek

druga pochodna potencjału delta D,

Podobny zabieg dla MRS: 3 węzły.

Metoda różnic skończonych, siatka nierównomierna

Iloraz różnicowy drugiej pochodnej dla nierównej siatki:

Δl Δp + tracimy jeden rząd dokładności w porównaniu z siatką równomierną Problem rozwiązany w metodzie elementów skończonych. Wzór trójpunktowy

W MES: nie ma problemu

Δ

p=1-q

Δ

l=q+1

-1

q

+1

v(x

)

u(q)

rozwiązanie dokładne

liniowe fcje kszałtu a warunki Neumanna

tylko v₁ oraz v₂ wnoszą przyczynek do pochodnej na lewym końcu:

pierwszy wiersz macierzy S

(-1/h₂ 1/h₂ 0 0 0 .... ) prawa strona pierwszy wiersz F : C

wybrane narzędzia MES umożliwiające

jej automatyzację w więcej niż 1D:

1) macierze sztywności pojedynczych elementów

oraz ich

2) składanie do globalnej macierzy sztywności

Przestrzeń referencyjna [odniesienia]

Problem fizycznie zadany jest na siatce [x₁,x₂,x₃,...x_N] Rachunki (całkowanie elementów macierzowych)

dla każdego elementu chcemy przenieść do przedziału (-1,1)

Element K_m=(x_m-1, x_m) → (-1,1) mapowanie z (-1,1) do K_m: x=(x_m+x_m-1)/2+(x_m-x_m-1)/2 ξ , gdzie ξ z przedziału (-1,1) Modelowy operator w 1D x= -1 1 y= -1 y=1

będziemy całkować jego elementy macierzowe w przestrzeni odniesienia

element w przestrzeni fizycznej

element w przestrzeni odniesienia

x(

ξ

=

ξ

skala transformacji m-tego elementu: (czynnik skali, jakobian)

przy transformacji: granice całki

zmieniają się na –1,1, poza tym dx=J_mdξ transformacja pochodnych:

Całkowanie macierzy sztywności w przestrzeni referencyjnej

1D: J nie zależy od ξ w 2D: zobaczymy, że nie zawsze tak jest [gdy element zmienia

swój kształt w mapowaniu. w 1D: odcinek -> odcinek] całkę i pochodne przenosimy do przestrzeni odniesienia:

element macierzowy całkowany w elemencie [fizycznym]

J_m=(x_m-x_m-1)/2

pole elementu fizycznego / pole elementu odniesienia -1

x(

ξ

=

ξ

całkowanie wektora sztywności: całka (f,v_j) transformuje się jak wyraz z a₀.

Całkowanie macierzy sztywności w przestrzeni referencyjnej

odcinkowo liniowe funkcje kształtu w przestrzeni odniesienia

x(

ξ

)

=

(x

+x

)/2+(x

- x

)/2

ξ

v_i(x(

ξ

v_i+1(x(

ξ

dwie funkcje kształtu

x v_i(x) 1 fcja kształtu x_i-1 K_i+1 K_i x_i x_i+1

fizyczna

K

-1 1

v

v

odniesienia

Przykład całkowanie w przestrzeni odniesienia dla bazy odcinkami liniowej (całka po elemencie K_i+1)

v_i(x(

ξ

ξ

J_i+1=(x_i+1-x_i)/2

pole elementu fizycznego / pole elementu odniesienia

tu prim to pochodna po x

a tu po ξ

Macierz sztywności pojedynczego elementu

składanie macierzy globalnej

Zmieniamy punkt widzenia:

(z funkcji kształtu na elementy) x_m-1 x_m

-1 1 _J m=hm/2 um₍

ξ

φ

ξ

φ

ξ

ξ

ξ

φ

φ

Zmieniamy punkt widzenia:

(z funkcji kształtu na elementy) x_m-1 x_m

-1 1 _J

m=hm/2

macierz sztywności elementu m

[wymiar taki jak liczba funkcji kształtu na element] um₍

ξ

φ

ξ

φ

ξ

φ

φ

zależność od m w J

:

ξ

=

ξ

na przekrywających się węzłach:suma

Em₌ Em11 Em12

Em

21 Em22

węzły na granicy elementów

obsługują więcej niż jeden element

Składanie (assembly) globalnej macierzy sztywności

macierz globalna S (rozmiar = liczbie węzłów)

U

U

U

U

1 2 3 4

ξ

φ

ξ

φ

ξ

u

_{= u}

globalna [U] i lokalna [u] numeracja

węzłów

u(x=-1)=0 u(x=1)=0 Przedział (-1,1) Podzielony na 7 elementów (8 węzłów) x_m-1 x_m -1 1 u₁ u₂ J_m=h_m/2 u(

ξ

φ

ξ

φ

ξ

φ

φ

case study

1 2 3 4

h₂ h₃ h₄ Składanie (assembly) macierzy sztywności z całek po elementach

Element 2

Element 3

Forma już znana

dodajemy elementy

z różnych macierzy lokalnych

Wektor obciążeń pojedynczego elementu/składanie globalnego

po elemencie K_i po K_i+1

x_i

x_i-1 x_i+1

U_i-1 U_i U_i+1

druga funkcja elementu i

pętla po elementach m=1,M

pętla po węzłach lokalnych k=1,N pętla po węzłach lokalnych l=1,N

identyfikacja numeru globalnego węzła

i=nr(k,m) j=nr(l,m)

S(i,j)=S(i,j)+E(m,k,l)

Składanie globalnej macierzy sztywności – buchalteria węzłów

m=1 m=2 m=3 m-numeruje elementy g=1 _g=2 g=3 _g=4 g – globalna numeracja węzłów 1 2 1 2 1 2 lokalne numery węzłów

nr (k,m) – przyporządkowanie numeru globalnego

węzłowi o lokalnym numerze k w elemencie m 1 = nr (1,1)

2 = nr (2,1) = nr (1,2) 3 = nr (2,2) = nr (1,3) 4 = nr (2,3)

identycznie składa się macierze dla wyższych funkcji kształtu i w więcej niż 1D

pętla po elementach m=1,M

pętla po węzłach lokalnych k=1,N pętla po węzłach lokalnych l=1,N

identyfikacja numeru globalnego węzła

i=nr(k,m) j=nr(l,m)

S(i,j)=S(i,j)+E(m,k,l)

Składanie globalnej macierzy sztywności – buchalteria węzłów

m=1 m=2 m=3 g=1 _g=2 g=3 _g=4 1 2 1 2 1 2 1 = nr (1,1) 2 = nr (2,1) = nr (1,2) 3 = nr (2,2) = nr (1,3) 4 = nr (2,3)

pętla po elementach m=1,M

pętla po węzłach lokalnych k=1,N

identyfikacja numeru globalnego węzła

i=nr(k,m)

F(i)=F(i)+P(k,m)

o potrzebie używania wyższych funkcji kształtu (i o laboratorium):

z liniowymi funkcjami kształtu: poza węzłami

nie uzyskamy dokładnego rozwiązania tego równania

(nigdy nie uzyskamy rozwiązania silnej postaci równania, druga pochodna wewnątrz elementów jest zawsze równa zeru a ma być równa niejednorodności dla równań elektrostatyki – źródło potencjału, dla równania przew. ciepl. – źródło ciepła)

u(1)=u(-1)=0

0.0

0.2

0.4

0.6

0.8

1.0

bx

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

a

całka działania a rozkład elementów dla funkcji odcinkowo linowych:

dokładne działanie

optymalne rozwiązanie odcinkami liniowe -1.0 -0.5 0.0 0.5 1.0 -0.4 -0.2 0.0 0.2 0.4 u rozwiązanie dokładne i=2,8

Funkcje kształtu wyższych rzędów:

u₁ u2 u4

funkcje kształtu

jeden element, cztery (n) węzły

u₃

ξ

_ξ

2

ξ

ξ

Funkcje kształtu wyższych rzędów:

u₁ u2 u4

funkcje kształtu

jeden element, cztery (n) węzły

u₃

ξ

_ξ

2

ξ

ξ

wielomian stopnia n-1, taki, że

wiemy jak go wskazać:

wielomian węzłowy

Lagrange’a

Funkcje kształtu wyższych rzędów:

u₁ u2 u4

funkcje kształtu

jeden element, cztery (n) węzły

u₃

ξ

_ξ

2

ξ

ξ

wielomian stopnia n-1, taki, że

wiemy jak go wskazać:

wielomian węzłowy

Lagrange’a

funkcje kształtu Lagrange’a: rozwiązanie interpolowane wielomianowo

w każdym z elementów. jedynie ciągłość rozwiązania między

elementami. w przeciwieństwie do problemów z KSN: wartości funkcji

w węzłach

nie są znane

. należy je wyliczyć. istota FEM.

Elementy wyższych rzędów:

Jeden element, trzy funkcje bazowe, 3 parametry węzłowe

Funkcje bazowe : w danym węźle tylko jedna z nich niezerowa

(co min. gwarantuje liniową niezależność funkcji bazowych)

φ

₁

=

ξ

(

ξ

-1)/2

φ

₂

=-(

ξ

-1)(

ξ

+1)

φ

₃

=

ξ

(

ξ

+1)/2

funkcje wierzchołkowe (vertex functions) 1 na krawędziach elementu

funkcja bąbelkowa (bubble function) węzeł wewnątrz elementu

-1 0 1

ξ

-1.00 -0.50 0.00 0.50 1.00 0.00 0.20 0.40 0.60 0.80 1.00

węzły= granice elementów

-0.80 -0.60 -0.40 -0.20 -0.40 0.00 0.40 0.80 1.20

liniowa baza Lagrange’a kwadratowa baza Lagrange’a każda strzałka to węzeł

granice elementów czerwone

φ

₁

=

ξ

(

ξ

-1)/2

φ

₂

=-(

ξ

-1)(

ξ

+1)

φ

₃

=

ξ

(

ξ

+1)/2

liczone numerycznie metodą Gaussa

całki wyliczone analitycznie: ilu punktowego Gaussa

należałoby użyć aby dokładnie

scałkować m.sztywności numerycznie?

przy równym podziale przedziału E takie samo dla każdego elementu

lecz P nie! [inny zakres x(ξ)]

x(ξ )=(x_m+x_m+1)/2+(x_m+1 - x_m)/2 ξ

Macierz sztywności

Liczba wierszy: 2n+1 (n-liczba elementów)

Składanie globalnej macierzy sztywności i wektora obciążeń

dla kwadratowych funkcji Lagrange’a

lokalne

-1.00 -0.50 0.00 0.50 1.00 -0.40 -0.20 0.00 0.20 0.40 Laboratorium: 4 elementy: x x x

2 centralne elementy o długości bx krzyżyki: węzły bąbelkowe

dokładne

MES dla bx=0.5 (równoodległe węzły) widzimy: dokładne

dla węzłów granicznych

Laboratorium: 0.00 0.20 0.40 0.60 0.80 1.00 bx -1.00 -0.80 -0.60 -0.40 -0.20 4 elementy:

dokładne 2 centralne elementy o długości bx krzyżyki: węzły bąbelkowe

x x x x

a

Wyniki dla problemu modelowego

funkcje odcinkami liniowe 9 węzłów (8 elementów)

funkcje kwadratowe 9 węzłów (4 elementy)

rozwiązanie dokładne -1.00 -0.50 0.00 0.50 1.00 -0.20 -0.10 0.00 0.10 0.20

Funkcje liniowe i kwadratowe -1.00 -0.50 0.00 0.50 1.00 -0.01 -0.00 0.00 0.00

0.01 funkcje odcinkami liniowe 9 węzłów (8 elementów)

funkcje kwadratowe 9 węzłów (4 elementy) [ błąd =0 tylko na węzłach granicznych ]

jeden rząd funkcji kształtu więcej:

maksymalny błąd

zmniejszony trzykrotnie

rozmiar problemu liniowego bez zmian, ale S

ma więcej niezerowych elementów

-1.00 -0.50 0.00 0.50 1.00 -0.20 -0.10 0.00 0.10 0.20

dla kwadratowej bazy Lagrange’a pochodna jest nieciągła na granicy elementów

-1.00 -0.50 0.00 0.50 1.00 -0.00 -0.00 0.00 0.00 0.00 nieciągła pochodna