MODELOWANIE INŻYNIERSKIE ISNN 1896-771X 32, s. 521-528, Gliwice 2006
EFEKTYWNOŚĆ SZACOWANIA BŁĘDÓW A POSTERIORI METODĄ WYRÓWNOWAŻONYCH RESIDUÓW ELEMENTOWYCH
W ADAPTACYJNEJ ANALIZIE PŁYT I POWŁOK
GRZEGORZ ZBOIŃSKI
Instytut Maszyn Przepływowych PAN w Gdańsku, Uniwersytet Warminsko-Mazurski w Olsztynie
Streszczenie. Prezentowane badania dotyczą globalnej (w strukturze) i lokalnej (w elemencie skończonym) efektywności szacowania błędów w adaptacyjnej analizie struktur cienkosciennych. Interesuje nas residualna (resztkowa) metoda szacowania błędów a posteriori zwana techniką wyrównoważania residuów.
Porównano wyniki dwóch testów adaptacyjnych. Uzyskane w obu testach współczynniki efektywności ilustrują problem lokalnego pogorszenia jakości estymacji. Przedstawiono zródła takiego pogorszenia oraz ewentualne sposoby zaradzenia mu w oparciu o wyrównoważanie wyższego stopnia.
1. WSTĘP
Prezentowane badania dotyczą globalnej i lokalnej (odnoszących się odpowiednio do struktury i do elementu skończonego) efektywności szacowania błędów w adaptacyjnej metodzie elementów skończonych. Efektywność ta jest mierzona za pomocą odpowiednich współczynników efektywności, określających stosunek oszacowanej wartości błędu do jego wartości dokładnej.
W niniejszej pracy interesują nas residualne (resztkowe) metody szacowania błędów a posteriori, a w szczególności tzw. technika wyrównoważania residuów [1]. Metoda ta została zmodyfikowana przez nas [2] w ten sposób, że błąd jest określony jako różnica pomiędzy energiami potencjalnymi odpowiadającymi rozwiązaniu ścisłemu i rozwiązaniu numerycznemu (a nie jako ekwiwalentna różnica energii odkształcenia), a wielkościami poszukiwanymi w problemach lokalnych elementów jest nie samo oszacowanie błędu, ale oszacowanie rozwiązania dokładnego. Tak zmodyfikowana metoda została następnie przystosowana przez nas [2, 3, 4, 5] do analizy płyt i powłok modelowanych zarówno elementami bryłowymi, powłokowymi hierarchicznymi, jak i powłokowymi pierwszego rzędu. Warto tutaj zaznaczyć, że oszacowania błędów aproksymacji obliczone za pomocą podejścia oryginalnego [1] i zmodyfikowanego [2] są identyczne. Podejście zmodyfikowane ma jednak tę przewagę, iż można je wykorzystać także do szacowania błędu całkowitego i błędu modelowania, czego nie zapewnia podejście oryginalne. Przypomnijmy tutaj, że błąd całkowity stanowi sumę błędu aproksymacji i błędu modelowania wynikającego z zastosowanie modelu mechanicznego uproszczonego względem modelu ścisłego, za jaki uznawać tutaj będziemy model trójwymiarowej teorii sprężystości.
Zastosowana metoda szacowania składa się z kilku etapów [5], wśród których najważniejszy polega na rozwiązywaniu problemów lokalnych w elementach, w oparciu o wyznaczone wcześniej liniowe funkcje rozkładu naprężeń międzyelementowych, wchodzące do definicji wyrównoważonych residuów elementowych. Jak to już wspomniano, uzyskane rozwiązania lokalne w elementach reprezentują oszacowane wartości rozwiązania dokładnego.
Ich porównanie z globalnym rozwiązaniem numerycznym, uzyskanym metodą elementów skończonych, pozwala na wyznaczenie elementowych indykatorów błędu oraz globalnego oszacowania błędu.
W naszych badaniach interesować nas będzie pogorszenie jakości procesu estymacji błędu w strukturach cienkościennych, wynikające z zastosowania hierarchicznych modeli powłokowych lub modeli trójwymiarowych, i obserwowane w stosunku do analiz w oparciu o modele powłokowe pierwszego rzędu.
2. METODA SZACOWANIA BŁĘDÓW I JEJ EFEKTYWNOŚĆ 2.1. Metoda wyrównoważania residuów elementowych
Jak to już zaznaczono we wstępie metoda wyrównoważania residuów została przez nas przedstawiona w pracy [5]. Tam też znaleźć można inne źródła literaturowe opisujące tę metodę w szczegółach. Tutaj ograniczymy się do przywołania podstawowej zależności charakteryzującej metodę:
(1)
Powyżej e oznacza numer elementu skończonego,
Π to energia potencjalna w elemencie, e B e to jego podwojona energia odkształcenia,
e
V i
e
S to jego objętość i brzeg, podczas gdy P i Q oznaczają odpowiednio obciążoną i utwierdzoną część brzegu całego ciała. Z kolei e a reprezentuje wektor błędu aproksymacji, na którym zdefiniowana jest globalna norma błędu, a
) ( M q
ηa to globalny estymator tegoż błędu. Ostatni składnik drugiej linii powyższej zależności definiuje tzw. residuum elementowe, w którym uq( M) to poszukiwane rozwiązanie dokładne, a
hp M q( ),
u to rozwiązanie numeryczne.
W pracy [5] znaleźć można definicje energii potencjalnej i energii odkształcenia odpowiadające zależności (1). Tam też znajdują się analogiczne zależności dla przypadków błędu całkowitego i błędu modelowania.
( )
( ) 22 ) (
), ( ), ( )
( ) (
) (
\
), ( )
( ),
( ), ( )
(
), ( )
2 ( ), (
] ) ,
( ) , ( [
] ) (
) ( 2
) ,
( ) ( 2 [
) (
2 ) ( 2 )
(
M q a e
M q a e e
hp M q hp M e q M q M e q
e Q
P S
hp M e q T M q e
hp M q hp M e q M e q
hp M q M
q V
hp M q a
B B
dS B
e e
e
η η
Π
Π Π
=
=
−
=
+
−
−
≤
+
−
=
∑
∑
∑ ∫
∪
∑
u u
u u
u r u
u u
u
u u
u e
2.2. Liniowe funkcje rozkładu naprężeń międzyelementowych
W zależności (1) występuje składnik odkreślony przez nas jako residuum elementowe. Jego drugi czynnik reprezentuje wektor tzw. wyrównoważonych naprężeń międzyelementowych.
Wektor taki definiuje się dla elementu e w oparciu o jego tensor naprężeń σe oraz analogiczny tensor σf w elemencie f posiadającym z nim wspólną ścianę, przy czym wykorzystujemy tutaj także wektor normalny do wspólnej ściany oznaczony symbolem ne . Iloczyny odpowiednich tensorów i wektora normalnego do brzegu prowadzą do wektorów naprężeń re i rf . W definicji wektora wyrównoważonych naprężeń międzyelementowych występują także wektorowe funkcje rozkładu tych naprężeń αe i αf , przy czym obowiązuje: αf = 1−αe . Definicja rozważanego wektora jest następująca:
(2)
Wektorowa liniowa funkcja rozkładu zdefiniowana jest na wspólnej ścianie sąsiednich elementów w oparciu o funkcje kształtu χ na ścianie elementu oraz współczynniki rozkładu i
α , wyznaczone w węzłach narożnych i ściany w oparciu o procedurę przedstawioną w [5]: i
(3)
gdzie w oznacza liczbę węzłów wierzchołkowych na wspólnej ścianie elementów e i f. Sv
2.3. Współczynniki efektywności
Miarą efektywności oszacowań błędu aproksymacji, uzyskanych w oparciu o metodę wyrównoważania residuów elementowych, będą zdefiniowane w sposób standardowy współczynniki efektywności, globalny θaq( M) i elementowe qa( M)
θe , zdefiniowane jako stosunek oszacowanej wartości błędu do jego wartości dokładnej
(4)
W przedstawionych poniżej testach numerycznych wartości błędu przyjęte za dokładne wynikać będą z porównania badanego rozwiązania numerycznego i najlepszego rozwiązania numerycznego możliwego do uzyskania (maksymalne zagęszczenie siatki i maksymalny stopień wielomianów aproksymujących).
)]
( ) (
)[
( )]
( ) (
[
) (
) (
) (
) (
) (
), ( ),
( ),
( ),
( 2 1
), ( ),
( ),
( ),
( ),
(
hp M f q hp M e q hp e
M f q hp M e q
hp M f q hp e M e q e f hp M f q e hp e M e q hp f
M e q
u r u
α r u
r u
r
u αr u
αr n σ u
α n σ u
α u
r
−
−
− +
=
+
= +
=
1
∑=
=
v
wS
i i i e
1
α
α χ
∑
≠
=
=
∑ e
M q a M e
q a
V hp M a q e
M q a e M
q a e
V hp M q a
M q M a
q a
e e
e
) ( )
(
), (
) ( )
( ),
( ) ( )
( ,
) (
, ) (
θ η θ
η θ θ
u u e
e
3. TESTY ILUSTRUJĄCE ZJAWISKO POGORSZENIA EFEKTYWNOŚCI 3.1. Płyta modelowana w oparciu o teorię pierwszego rzędu
Pierwszy z testów dotyczy symetrycznej ćwiartki utwierdzonej płyty kwadratowej obciążonej równomiernie siłami powierzchniowymi. Przyjęty do obliczeń model mechaniczny odpowiada teorii pierwszego rzędu Reissnera-Mindlina. Rys. 1 przedstawia stan naprężeń efektywnych uzyskanych dla siatki zadaptowanej automatycznie (przez program komputerowy). Przypomnijmy tutaj, że stosowana przez nas procedura adaptacyjna składa się z czterech kroków, odpowiadających sieci początkowej, zmodyfikowanej, pośredniej i końcowej. Procedura ta opisana jest w pracy [2].
Podstawową cechą uzyskanego rozwiązania jest jego gładki charakter, przejawiający się brakiem istotnych koncentracji naprężeń, co widać na rysunku. Globalna oszacowana wartość błędu dla tego przykładu wynosi 0,4%, a uzyskany globalny współczynnik efektywności jest równy 1,658 (patrz rys.2).
Rys.1. Stan naprężeń efektywnych w symetrycznej ćwiartce płyty
Rys.2. Współczynniki efektywności w symetrycznej ćwiartce płyty
3.2. Powłoka o złożonym modelu mechanicznym
Drugi test dotyczy symetrycznej powłoki półwalcowej obciążonej równomiernymi pionowymi siłami powierzchniowymi. Brzeg prostoliniowy powłoki jest utwierdzony, a brzeg zakrzywiony pozostaje swobodny. Powłoka ma charakter struktury złożonej ze względu na niejednorodny model mechaniczny, przejawiający się występowaniem w jednej sieci, elementów powłokowych pierwszego rzędu (środek powłoki), elementów powłokowych hierarchicznych i elementów bryłowych (brzegi: utwierdzony i swobodny powłoki), a także łączących je elementów przejściowych.
Rys. 3 przedstawia stan naprężeń efektywnych uzyskanych dla sieci końcowej wynikającej z czterokrokowego procesu adaptacji. Widać wyraźnie koncentracje naprężeń wzdłuż brzegu krzywoliniowego, wynikające z wystąpienia tzw. efektu brzegowego (inaczej, warstwy brzegowej), polegającego na wykładniczej zmianie naprężeń w kierunku normalnym do tego brzegu. Dodatkowo, w jednym z widocznych naroży powłoki zaobserwować można wystąpienie osobliwości rozwiązania, wiążącej się z gwałtownym (dążącym teoretycznie do nieskończoności) wzrostem naprężeń.
Uzyskany dla niniejszego przykładu rzeczywisty poziom błędu globalnego w całej strukturze wynosi 5,4%, a odpowiadający mu globalny współczynnik efektywności metody wyrównoważonych residuów elementowych jest równy 7,298 (patrz rys. 4). Warto tutaj zwrócić uwagę na fakt, iż lokalnie, w elemencie skończonym modelującym naroże z osobliwością, rzeczywisty błąd jest równy 11%, a lokalny współczynnik efektywności przekrocza 100, co oznacza, iż błąd w tym elemencie jest przeszacowany ponad stukrotnie.
Rys.3. Stan naprężeń efektywnych w symetrycznej ćwiartce powłoki
Porównanie wyników pierwszego i drugiego testu wskazuje, że metoda wyrównoważania residuów elementowych, oparta na liniowych funkcjach rozkładu, nie radzi sobie z przypadkami, kiedy rozwiązanie nie ma charakteru gładkiego.
Rys.4. Współczynniki efektywności w symetrycznej ćwiartce powłoki
4. PROPOZYCJA POPRAWY EFEKTYWNOŚCI I JEJ PRZESŁANKI
W celu poprawienia efektywności szacowania błędu aproksymacji autor proponuje zastosowanie wyrównoważania wyższego stopnia, które wykorzystywać będzie wielomianowe funkcje rozkładu naprężeń międzyelementowych. Takie podejście wymaga wyznaczenia węzłowych współczynników rozkładu nie tylko w węzłach narożnych, jak to ma miejsce w przypadku liniowych funkcji rozkładu, ale także w węzłach wyższego stopnia, tzn.
na krawędziach i w środku ścian. Sytuację tę zilustrowano na rys. 3, dla przypadku dwuwymiarowego. Zauważmy, iż na rysunku tym dowolnie zmienną liczbę węzłów krawędziowych oznaczono za pomocą kresek. W przypadku zaproponowanego podejścia, modyfikacji ulec musi także zależność (3). Przyjmie ona teraz postać następującą:
(5)
gdzie w oznacza dodatkową liczbę węzłów wyższego stopnia na ścianie elementu. Sh
Niniejsza propozycja wykorzystania funkcji rozkładu bierze za podstawę następujące przesłanki wynikające z definicji (patrz [1, 2]) współczynników rozkładu naprężeń między- elementowych oraz przeprowadzonych testów pilotażowych. Otóż, liniowe współczynniki rozkładu są tak zdefiniowane, że błąd jaki powodują w przypadku zastosowania ich do elementów liniowych (pierwszego stopnia) jest równy zeru. Z tego powodu w problemach lokalnych definiuje się elementy skończone jako o jeden stopień wyższe niż w problemie globalnym, na podstawie którego te współczynniki są wyznaczane. Autor przeprowadził
∑+
= =
h S v
S w
w i
i i e
1
α
α χ
wstępny test dla przypadku współczynników drugiego stopnia zastosowanych do szacowania błędów aproksymacji w elementach drugiego stopnia. Także i w tym przypadku błędy te okazały się być równe zeru. Świadczy to o poprawności naszej propozycji. Z kolei dokładniejsza, niż w przypadku liniowym, aproksymacja funkcji rozkładu stanowi przesłankę do oczekiwania wyższej efektywności szacowania błędów.
Rys.3. Węzły współczynników rozkładu w przypadku liniowym i wyższego stopnia
3. WNIOSKI
Przyczyną pogorszenia jakości szacowania błędu aproksymacji w problemach struktur cienkościennych, w przypadku zastosowania hierarchicznych modeli powłokowych zamiast modelu powłokowego pierwszego rzędu, jest różny charakter odpowiadających im rozwiązań.
W drugim przypadku bowiem mamy do czynienia z rozwiązaniem gładkim, podczas gdy w pierwszym wystąpić mogą spiętrzenia naprężeń w postaci tzw. warstwy brzegowej i osobliwości wierzchołkowych lub krawędziowych. W wypadku ich pojawienia się dochodzi do lokalnego pogorszenia efektywności zastosowanej residualnej metody szacowania błędu, co wpływa także na pogorszenie globalnej efektywności takiego oszacowania.
Jako jedną z możliwych metod poprawienia efektywności szacowania błędu aproksymacji, w sytuacji wystąpienia osobliwości lub efektu brzegowego, zaproponowano zastosowanie wyrównoważania wyższego stopnia z funkcjami rozkładu naprężeń międzyelementowych w postaci wielomianów, w miejsce wykorzystywanego dotychczas przez autora i poprzedników wyrównoważania pierwszego stopnia z liniowymi funkcjami rozkładu. Przesłanki takiego podejścia wynikają z rozważań teoretycznych oraz ze wstępnych testów numerycznych, których wyniki omówiono w punkcie poprzednim.
LITERATURA
1. Ainsworth, M., Oden, J. T., Wu, W., A posteriori error estimation for hp approximation in elastostatics. Appl. Numer. Math., 14, 1994, 23-55.
2. Zboiński G.: Modelowanie hierarchiczne i metoda elementów skończonych do adaptacyjnej analizy struktur złożonych (rozprawa habilitacyjna). Zesz. Nauk. IMP PAN.
Studia i Materiały, 520/1479/2001. Gdańsk: Wydawnictwo IMP PAN.
3. Zboiński G.: A posteriori error estimation for hp-approximation of the 3D-based first order shell model. Part I. Theoretical aspects. Applied Mathematics Informatics and Mechanics, 8, 1, 2003, 104-125.
4. Zboiński G.: A posteriori error estimation for hp-approximation of the 3D-based first order shell model. Part II. Implementation aspects. Applied Mathematics, Informatics and Mechanics, 8, 2, 2003, 59-83.
5. Zboiński G.: Adaptacyjna analiza struktur złożonych. Szacowanie błędów a posteriori, Zesz. Nauk. Katedry Mechaniki Stosowanej, 23, 2004, 501-506.
EFFECTIVITY OF A POSTERIORI ERROR ESTIMATION WITH ELEMENT RESIDUAL EQULIBRATION METHOD
IN ADAPTIVE ANALYSIS OF PLATES AND SHELLS
Summary. The presented research deals with the global (within the structure) and local (within a finite element) effectivity of a posteriori error estimation in adaptive analysis of thin-walled structures. We are interested in the Element Residual Method (ERM) of a posterirori error estimation with equilibration of the residuals.The results of two adaptive tests have been compared. The effectivity indices from both examples reflect the problem of local worsening of the quality of error estimation. The source of such local worsening has been indicated and the method of overcoming it, based on higher order equilibration, has been suggested.
