MES w zagadnieniach nieliniowych

MES w zagadnieniach nieliniowych

Jerzy Pamin e-mail: JPamin@L5.pk.edu.pl

Podziękowania:

A. Wosatko, A. Winnicki

ADINA R&D, Inc.http://www.adina.com ANSYS, Inc. http://www.ansys.com TNO DIANA http://www.tnodiana.com Altair Engineering http://www.comco.com

Tematyka zajęć

Zagadnienia nieliniowe Analiza przyrostowo-iteracyjna Nieliniowość geometryczna Nieliniowość fizyczna Zarysowanie

Katastrofy budowlane Literatura

[1] R. de Borst and L.J. Sluys. Computational Methods in Nonlinear Solid Mechanics. Lecture notes, Delft University of Technology, 1999.

[2] G. Rakowski, Z. Kacprzyk. Metoda elementow skończonych w mechanice kostrukcji. Oficyna Wyd. PW, Warszawa, 2005.

[3] M. Jir´asek and Z.P. Baˇzant. Inelastic Analysis of Structures. J. Wiley &

Sons, Chichester, 2002.

[4] M. Kwasek Advanced static analysis and design of reinforced concrete deep beams. Diploma work, Politechnika Krakowska, 2004.

Metody komputerowe, studia II st.

Źródła nieliniowości

Spowodowane zmianą geometrii ciała (odkształcalnego)

I duże odkształcenia (np.guma, formowanie metali)

I duże przemieszczenia (np. konstrukcje smukłe,cienkościenne)

I kontakt (oddziaływanie stykających się ciał)

I obciążenie śledzące (zależne od deformacji ciała) Spowodowane nieliniowymi związkami konstytutywnymi

I plastyczność (odkształcenia trwałe)

I uszkodzenie (degradacja własności sprężystych)

I zarysowanie (kontynualna reprezentacja rys)

I . . . Uwagi:

I Nie obowiązuje zasada superpozycji.

I Możliwy jest opis ośrodka nieciągłego, w którym części składowe są połączone interfejsami (np. konstrukcje zespolone) lub występują pękniecia (rysy dyskretne). Interfejsy mają zazwyczaj nieliniowe charakterystyki, reprezentując np. tarcie, adhezję, pękanie.

Nieliniowe kontinuum [1,2,3]

Równania równowagi + statyczne warunki brzegowe L^Tσ + b = 0 w V , σν = ˆt na S gdzie:

L – macierz operatorów różniczkowych σ – tensor/wektor uogólnionych naprężeń b – wektor sił masowych

ν

V S ˆ_t

Słaba forma równań równowagi Z

V

δu^T(L^Tσ + b) dV = 0 ∀δu Zasada prac wirtualnych δWint = δWext

Z

V

(Lδu)^Tσ dV = Z

V

δu^Tb dV + Z

S

δu^Tˆt dS

Metody komputerowe, studia II st.

Metoda Galerkina

Przemieszczeniowa wersja MES

u ≈ u^h=

nw

X

i =1

N_i(ξ, η, ζ)u_i= Nu^e

gdzie: N - funkcje kształtu, u^e - wektor stopni swobody elementu, nw - liczba węzłów

Transformacja węzłowych stopni swobody u^e= A^eu^g gdzie: u^g - wektor stopni swobody układu

Słaba forma równań równowagi dla układu zdyskretyzowanego

ne

X

e=1

A^{e T} Z

V^e

B^Tσ dV = f_ext, B = LN

Podejście izoparametryczne, numeryczne całkowanie macierzy ES

Liniowa sprężystość

Prawo Hooke’a

Notacja tensorowa: σ = D^e : , σij = D_ijkl^e kl

Notacja macierzowa:

σ = D^e, σ =









σx

σy

σz

τxy

τyz

τzx







 ,  =









x

y

z

γxy

γyz

γzx







 Izotropia materiału: D^e= D^e(E , ν)

E 1

σ

Liniowe związki kinematyczne 

Notacja tensorowa:  = ¹₂[∇u + (∇u)^T], _ij =¹₂(u_{i ,j}+ u_{j ,i}) Notacja macierzowa:  = Lu

Zatem tensor naprężenia:

σ = D^e = D^eLu = D^eLNu^e = D^eBA^eu^g

Równania równowagi dla układu zdyskretyzowanego

n_e

X

e=1

A^{e T} Z

V^e

B^TD^eB dV A^e u^g = fext, Ku^g = fext

Metody komputerowe, studia II st.

Analiza przyrostowo-iteracyjna

Nieliniowy problem:

fext przykładane w przyrostach t → t + ∆t → σ^t+∆t= σ^t+ ∆σ Równowaga w chwili t + ∆t:

ne

X

e=1

A^{e T} Z

V^e

B^Tσ^t+∆tdV = f_ext^t+∆t

n_e

X

e=1

A^{e T} Z

V^e

B^T∆σ dV = f_ext^t+∆t− f_int^t gdzie: f_int^t =Pne

e=1A^{e T}R

V^eB^Tσ^tdV Linearyzacja lewej strony w chwili czasu t:

∆σ = ∆σ(∆(∆u)) Układ równań dla przyrostu:

K ∆u^g = f_ext^t+∆t− f_int^t

Schemat metody przyrostowo-iteracyjnej

Konieczne poprawki iteracyjne celem osiągnięcia stanu równowagi w chwili t + ∆t → algorytm Newtona-Raphsona

Siły niezrównoważone (residualne): Rj = f_ext^t+∆t− f_int,j^t+∆t→ 0 f

R1

f_int,1

u u^t+∆t

u^t f_ext^t f_ext^t+∆t

∆u₁ du₂

∆f_ext

Kjdu^g= f_ext^t+∆t− f_int,j^t+∆t K- operator styczny Pierwsza iteracja:

∆u^g₁=K⁻¹₀ (f_ext^t+∆t− fint,0^t ) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Poprawki:

du^g_{j +1}= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Kryterium zbieżności:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ

Algorytm zmodyfikowany:

Kj= K0

Metody komputerowe, studia II st.

Schemat metody przyrostowo-iteracyjnej

Konieczne poprawki iteracyjne celem osiągnięcia stanu równowagi w chwili t + ∆t → algorytm Newtona-Raphsona

Siły niezrównoważone (residualne): Rj = f_ext^t+∆t− f_int,j^t+∆t→ 0 f

R1

f_int,1

u u^t+∆t

u^t f_ext^t f_ext^t+∆t

∆u₁ du₂

∆f_ext

Kjdu^g= f_ext^t+∆t− f_int,j^t+∆t K - operator styczny Pierwsza iteracja:

∆u^g₁ = K⁻¹₀ (f_ext^t+∆t− fint,0^t ) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Poprawki:

du^g_{j +1}= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Kryterium zbieżności:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ

Algorytm zmodyfikowany:

Kj= K0

Schemat metody przyrostowo-iteracyjnej

Konieczne poprawki iteracyjne celem osiągnięcia stanu równowagi w chwili t + ∆t → algorytm Newtona-Raphsona

Siły niezrównoważone (residualne): Rj = f_ext^t+∆t− f_int,j^t+∆t→ 0 f

R1

f_int,1

u u^t+∆t

u^t f_ext^t f_ext^t+∆t

∆u₁ du₂

∆f_ext

K ∆u^g= f_ext^t+∆t− fint^t

K - operator styczny Pierwsza iteracja:

∆u^g₁ = K⁻¹₀ (f_ext^t+∆t− fint,0) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Poprawki:

du^g_{j +1}= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Kryterium zbieżności:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ

Algorytm zmodyfikowany:

Kj= K0

Metody komputerowe, studia II st.

Schemat metody przyrostowo-iteracyjnej

Konieczne poprawki iteracyjne celem osiągnięcia stanu równowagi w chwili t + ∆t → algorytm Newtona-Raphsona

Siły niezrównoważone (residualne): Rj = f_ext^t+∆t− f_int,j^t+∆t→ 0 f

u u^t+∆t

u^t f_ext^t f_ext^t+∆t

∆u₁ du₂

∆f_ext

R₂

f_int,2

K ∆u^g= f_ext^t+∆t− fint^t

K - operator styczny Pierwsza iteracja:

∆u^g₁ = K⁻¹₀ (f_ext^t+∆t− fint,0) σ1 → f_int,1^t+∆t6= f_ext^t+∆t Poprawki:

du^g_{j +1}= K⁻¹_j (f_ext^t+∆t−f_int,j^t+∆t) σj +1 → f_{int,j +1}^t+∆t

Kryterium zbieżności:

kf^t+∆t_ext −f^t+∆t int,j k k∆f_extk ¬ δ

Algorytm zmodyfikowany:

Kj= K0

Sposoby przykładania przyrostów

Sterowanie siłą lub przemieszczeniem

Sterowanie parametrem łuku

Metody komputerowe, studia II st.

Nieliniowość geometryczna

X

u

x φ(X, t)

V

V0

S S₀

x1, X1

x2, X2

Konfiguracja początkowa i aktualna Funkcja ruchu: x = φ(X, t)

Wektor przemieszczenia: u(X, t) = x − X

Gradient deformacji (podstawowa miara deformacji): F =^∂φ_∂X = ∇Xx Tensor odkształcenia Greena (jedna z możliwych miar odkształcenia):

E = 1

2(F^TF − I) = 1

2[∇_Xu + (∇_Xu)^T+ (∇_Xu)^T∇Xu]

Nieliniowość geometryczna

Nieliniowe związki kinematyczne, np. ε_x = ε^L_x+ε^N_x =^∂u_∂x +¹₂ ^∂w_∂x
²

∆σ = ∆σ(∆(∆u))

I Równania równowagi opisują równowagę ciała zdeformowanego.

Zasadę prac wirtualnych można zapisać w konfiguracji początkowej lub aktualnej.

I Różnym miarom odkształcenia odpowiadają różne miary naprężenia.

I Małe odkształcenia: E ≈  = ¹₂[∇u + (∇u)^T] < 2%.

I Małe przemieszczenia (uogólnione): V ≈ V0(jeden opis, zasada zesztywnienia).

Metody komputerowe, studia II st.

Nieliniowość geometryczna

Równowaga układu zdyskretyzowanego:

K ∆u^g = f_ext^t+∆t− f_int^t gdzie styczna macierz sztywności:

K = K0+Ku+Kσ

K0- macierz sztywności liniowej

Ku - macierz sztywności przemieszczeniowej

(macierz dyskretnych związków kinematycznych B zależna od przemieszczeń)

Kσ - macierz sztywności naprężeniowej (zależna od naprężeń uogólnionych)

Nieliniowość fizyczna

K ∆u^g = f_ext^t+∆t− f_int^t Linearyzacja lewej strony w chwili czasu t:

∆σ =∆σ(∆(∆u))

∆σ = ^∂σ_∂
^t _∂

∂u

^t

∆u D =^∂σ_∂, L = ^∂_∂u

Dyskretyzacja: ∆u = N∆u^e

Liniowe związki geometryczne → macierz dyskretnych związków kinematycznych B = LN niezależna od przemieszczeń

Styczna macierz sztywności K =

n_e

X

e=1

A^{e T} Z

V^e

B^TDB dV A^e

Metody komputerowe, studia II st.

Uplastycznienie materiału

A B C

przemieszczenie siła

P

A

+

-

σy

σy

σy

σy

σy

σy

+

- -

+ C B

zakres sprężysty

pełne uplastycznienie pełne uplastycznienie zakres sprężysty

Rysy dyskretne lub rozmazane

Energia pękania Gf (zużyta na powstanie jednostki powierzchni rysy)

Metody komputerowe, studia II st.

Symulacja zarysowania żelbetowej tarczy

pakietem ATENA [4]

Katastrofa platformy Sleipner A, Norwegia 1991

I Żelbetowa platforma wirtnicza posadowiona na głebokosci 82 m, podstawa złożona z 24 komór o średnicy 12 m (4 wspierają pomost)

I Przyczyna zatonięcia konstrukcji podstawy podczas operacji posadowienia:

błąd w obliczeniach MES trójnika łączącego komory (niedoszacowanie siły ścinającej o 47%) i niewystarczające zakotwienie zbrojenia w strefie krytycznej

Rysunki z www.ima.umn.edu/∼arnold/disasters/sleipner.html

Metody komputerowe, studia II st.

Airport Paris Charles de Gaulle, Terminal 2E, 2004

I Zespolona przeszklona konstrukcja powłokowa w kształcie rury, swobodnie podparte sklepienie osłabione licznymi otworami

I Zaprojektowany przez architekta Paula Andreu (zaprojektował również terminal 3 w Dubai International Airport, który zawalił się podczas budowy), oddany w roku 2003

I Przyczyna: zbyt mały margines bezpieczeństwa w projekcie,

prawdopodobnie także błedy wykonawcze i/lub niedostatecznie dobry beton

Rysunki zaczerpnięte z www.equipement.gouv.fr

Uwagi końcowe

1. Symulacje komputerowe stwarzają bezcenne możliwości, ale tylko świadomemu użytkownikowi MES.

2. W modelowaniu efektów nieliniowych zaczyna dominować modelowanie 3D.

3. Dla poprawy jakości aproksymacji MES wskazane jest adaptacyjne zagęszczanie siatki na podstawie oszacowania błędu dyskretyzacji.

Metody komputerowe, studia II st.

Adaptacyjne zagęszczenie siatki elementów

Przykład zaczerpnięty ze strony Altair Engineering http://www.comco.com

Generacja siatki elementów

Przykład zaczerpnięty ze strony Altair Engineering http://www.comco.com

Metody komputerowe, studia II st.

Monitoring błędów dyskretyzacji

Adaptacyjne zagęszczenie siatki