Rozwiązywanie zagadnień nieliniowych

Rozwiązywanie zagadnień nieliniowych

Wykład 4 z MKwIL, kierunek Budownictwo

Jerzy Pamin

Katedra Technologii Informatycznych w Inżynierii Politechnika Krakowska

Podziękowania:

A. Wosatko, A. Winnicki, Ł. Kaczmarczyk ANSYS, Inc. http://www.ansys.com TNO DIANA http://www.tnodiana.com FEAP http://www.ce.berkeley.edu/feap ATENA http://www.cervenka.cz

MKwIL, Budownictwo 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ęknięcia (rysy dyskretne). Interfejsy mają zazwyczaj nieliniowe charakterystyki, reprezentując np. tarcie, adhezję, pękanie.

Nieliniowe kontinuum [1]

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 δW_int = δW_ext

Z

V

(Lδu)^Tσ dV = Z

V

δu^Tb dV + Z

S

δu^Tˆt dS

MKwIL, Budownictwo II st.

Metoda Galerkina

Przemieszczeniowa wersja MES u ≈ u^h =

nw

X

i =1

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

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

Transformacja węzłowych stopni swobody d^e = A^ed gdzie: d - wektor stopni swobody układu

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

n_e

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

σ



loading unloading

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^eLNd^e = D^eBA^ed

Równania równowagi dla układu zdyskretyzowanego

ne

X

e=1

A^{e T} Z

V^e

B^TD^eB dV A^e d = f_ext, Kd = f_ext

MKwIL, Budownictwo II st.

Analiza przyrostowo-iteracyjna

Nieliniowy problem:

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

n_e

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 = Pn_e

e=1A^{e T}R

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

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

K ∆d = 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): R_j = f_ext^t+∆t− f_int,j^t+∆t → 0

f

R₁

f_int,1

u u^t+∆t

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

∆u₁ du₂

∆fext

K ∆d = f_ext^t+∆t − f_int^t K - operator styczny Pierwsza iteracja:

∆d1=K⁻¹₀ (f_ext^t+∆t− f_int^t ) σ1 → f_int,1^t+∆t 6= f_ext^t+∆t Poprawki:

ddj +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_ext^t+∆t−f^t+∆t int,j k k∆f_extk ¬ δ

Algorytm zmodyfikowany:

Kj = K0

MKwIL, Budownictwo 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): R_j = f_ext^t+∆t− f_int,j^t+∆t → 0

f

R₁

f_int,1

u u^t+∆t

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

∆u₁ du2

∆f_ext

K ∆d = f_ext^t+∆t − f_int^t K - operator styczny Pierwsza iteracja:

∆d1 = K⁻¹₀ (f_ext^t+∆t− f_int,0^t ) σ1 → f_int,1^t+∆t 6= f_ext^t+∆t Poprawki:

ddj +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_ext^t+∆t−f_int,j^t+∆tk 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): R_j = f_ext^t+∆t− f_int,j^t+∆t → 0

f

R₁

f_int,1

u u^t+∆t

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

∆u₁ du₂

∆fext

K ∆d = f_ext^t+∆t − f_int^t K - operator styczny Pierwsza iteracja:

∆d1 = K⁻¹₀ (f_ext^t+∆t− fint,0) σ1 → f_int,1^t+∆t 6= f_ext^t+∆t Poprawki:

ddj +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_ext^t+∆t−f^t+∆t int,j k k∆f_extk ¬ δ

Algorytm zmodyfikowany:

Kj = K0

MKwIL, Budownictwo 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): R_j = 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₁ du2

∆f_ext

R₂

f_int,2

K ∆d = f_ext^t+∆t − f_int^t K - operator styczny Pierwsza iteracja:

∆d1 = K⁻¹₀ (f_ext^t+∆t− fint,0) σ1 → f_int,1^t+∆t 6= f_ext^t+∆t Poprawki:

ddj +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_ext^t+∆t−f_int,j^t+∆tk k∆f_extk ¬ δ

Algorytm zmodyfikowany:

Kj = K0

Sposoby przykładania przyrostów

Sterowanie siłą lub przemieszczeniem

Sterowanie parametrem łuku

MKwIL, Budownictwo II st.

Nieliniowość geometryczna

X

u

x φ(X, t)

V

V₀

S S₀

x₁, X₁ x₂, X₂

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).

MKwIL, Budownictwo II st.

Nieliniowość geometryczna

Równowaga układu zdyskretyzowanego:

K ∆d = f_ext^t+∆t − f_int^t

gdzie styczna macierz sztywności:

K = K₀+K_u +K_σ K₀ - macierz sztywności liniowej

K_u - 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)

Teoria umiarkowanie dużych ugięć Karmana [2]

Dopuszczamy ugięcie rzędu grubości Odkształcenia powierzchni środkowej ε_x = ε^L_x + ε^N_x = ^∂u_∂x + ¹₂ ^∂w_∂x
²

ε_y = ε^L_y + ε^N_y = ^∂v_∂y + ¹₂ 

∂w

∂y

² γ_xy = γ_xy^L + γ_xy^N = ^∂u_∂y + ^∂v_∂x + ^∂w_∂x ^∂w_∂y Krzywizny i spaczenie jak w teorii liniowej

κ_x = κ^L_x = −^∂_∂x^2w₂ , κ_y = κ^L_y = −^∂_∂y^2w₂, χ_xy = χ^L_xy = −2_{∂x ∂y}^∂2w

Układ dwu równań teorii Karmana płyt o umiarkowanie dużych ugięciach

∇²∇²F (x , y ) + ^Eh₂ L(w , w ) = 0 D^m∇²∇²w (x , y ) − L(w , F ) − ˆp_z = 0 gdzie:

F (x , y ) - funkcja naprężenia (n_x = F_,yy, n_y = F_,xx, n_xy = −F_,xy) L(a, b) = a_,xxb_,yy − 2a_,xyb_,xy + a_,yyb_,xx

MKwIL, Budownictwo II st.

Teoria umiarkowanie dużych ugięć

Aproksymacja MES

Całkowita energia potencjalna (dodatkowa energia stanu tarczowego) Π˜^m = U^m + ˜Uⁿ − W^m

Dyskretyzacja

uⁿ =

 u(x , y ) v (x , y )



=

 N_u 0 N_v 0

  dⁿ d^m



w^m = [w (x , y )] = 

0 N_w 

 dⁿ d^m



Nieliniowe związki kinematyczne B_{(6×LSSE )} = B^L_{(6×LSSE )}+ B^N_{(6×LSSE )}

B^L_{(6×LSSE )} =

 Bⁿ 0 0 B^m



, B^N_{(6×LSSE )} =

 0 Bⁿ_w 0 0



Teoria umiarkowanie dużych ugięć

Aproksymacja MES

Macierze dyskretnych związków kinematycznych

Bⁿ =





N_u,x N_{v ,y} N_u,y + N_{v ,x}



, B^m =





−N_{w ,xx}

−N_{w ,yy}

−2N_{w ,xy}





Bⁿ_w =





w_,xN_{w ,x} w_,yN_{w ,y} w_,xN_{w ,y} + w_,yN_{w ,x}



 Gradienty ugięcia

g =

 w_,x w_,y



=

 N_{w ,x} N_{w ,y}



d^m = G^md^m

MKwIL, Budownictwo II st.

Teoria umiarkowanie dużych ugięć

Aproksymacja MES

Macierz styczna elementu k^e_T = k^e₀+ k^e_u + k^e_σ

k^e₀ = Z Z

A^e

 B^nTDⁿBⁿ 0 0 B^mTD^mB^m

 dA

k^e_u = Z Z

A^e

 0 B^nTDⁿBⁿ_w B^nT_w DⁿBⁿ 0

 dA

k^e_σ = Z Z

A^e

 0 0

0 G^mTSⁿG^m



dA, Sⁿ =

 n_x n_xy n_xy n_y



Wektor sił węzłowych reprezentujących siły wewnętrzne f_int^e =

Z Z

A^e

 B^nT 0 B^nT_w B^mT

  Sⁿ S^m

 dA

Płyta kwadratowa

Umiarkowanie duże ugięcia

ˆ p_z

x

z, w y

a

a C

Dane:

L = L_x = L_y = 2a = 1.0 m h = 0.002 m

E = 200000 MPa, ν = 0.25

Zależność ugięcia od obciążenia:

0.001 0.003

0.30

0.15

0.002

w_C [m]

ˆ p_z [kPa]

MES

(nieliniowość geometryczna)

w_C = 0.00406^pˆz_D^L4

rozwiązanie liniowe:

MKwIL, Budownictwo II st.

Nieliniowość fizyczna

K ∆d = 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∆d^e

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

Styczna macierz sztywności K =

n_e

XA^{e T} Z

V^e

B^TDB dV A^e

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

MKwIL, Budownictwo II st.

Rysy dyskretne lub rozmazane (DIANA)

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

Symulacja zarysowania żelbetowej tarczy pakietem ATENA [3]

MKwIL, Budownictwo II st.

Nośność tarcz wieloprzęsłowych (DIANA) [4]

Geometrycznie i materiałowo nieliniowa analiza [5]

Schemat strategii obliczeniowej realizowanej na wielu poziomach:

I konstrukcji

I elementu skończonego I warstwy

I punktu

Uwzględnione efekty:

I śledzenie wytężenia w przekroju I zarysowanie betonu

I sprężysto-plastyczne zbrojenie I duże przemieszczenia i ich

gradienty Cele:

I wyznaczenie ewolucji przemieszczeń

I określenie mechanizmu uszkodzenia

I oszacowanie nośności

MKwIL, Budownictwo II st.

Model powłoki żelbetowej

I Zdegenerowany element powłokowy 8-węzłowy (teoria Mindlina-Reissnera)

I Model warstwowy powłoki żelbetowej (5 warstw betonu, 4 warstwy stali reprezentujące 2 siatki zbrojenia)

I Model sprężysty z rozmazanym zarysowaniem dla warstwy betonu (osłabienie betonu, redukcja sztywności na ścinanie)

I Model sprężysto-plastyczny dla warstw stali

Analiza numeryczna powłoki chłodni kominowej [5]

MKwIL, Budownictwo II st.

Analiza numeryczna powłoki chłodni kominowej

Zależności λ − wK otrzymane dwoma pakietami MES przy sterowaniu siłą lub przemieszczeniem dla obciążenia g + λ(w + s)

Obciążenia chłodni:

I ciężar własny g

I wiatr w

I ssanie wewnętrzne s

I obciążenia termiczne

I osiadania podłoża

Analiza numeryczna żelbetowej powłoki

Wyniki analizy numerycznej powłoki z otworem technologicznym

Deformacja Mapa warstwicowa membranowych sił południkowych

MKwIL, Budownictwo II st.

Analiza numeryczna żelbetowej powłoki

Kierunki naprężeń głównych w warstwie zewnętrznej Wizualizacja rozmazanych rys

Uszkodzona chłodnia kominowa [6]

MKwIL, Budownictwo II st.

Analizowane przypadki przy obciążeniu g + λ(w + s)

I chłodnia projektowana

I chłodnia wybudowana ze strefami słabego betonu (f_cm=11 MPa)

I chłodnia z dwoma otworami obwodowymi (długości 25m i 14m)

I chłodnia naprawiona i wzmocniona (5cm zbrojonego torkretu w

Liniowe wyboczenie powłoki chłodni (DIANA) [7]

Obciążenie Projektowana Skonstruowana Uszkodzona Naprawiona

λg 25.52 22.22 18.59 19.11

λ(w + s) 13.15 14.72 6.24 20.71

λ(g + w + s) 11.29 11.39 5.49 20.84

Krytyczne mnożniki obciążenia λ1

Forma wyboczenia dla powłoki uszkodzonej

MKwIL, Budownictwo II st.

Wyniki analizy nieliniowej [6]

Błąd technologiczny nie miał znaczącego wpływu na doraźną nośność, lecz na trwałość konstrukcji (korozję zbrojenia, lokalne przeciążenia betonu).

Predykcja stref zarysowania (DIANA)

Strefy zarysowania wewnętrznej i zewnętrznej warstwy betonu dla λ = 3.2 (DIANA)

MKwIL, Budownictwo II st.

Model chłodni z otworem (DIANA)

Obliczenia nie spełniają kryterium zbieżności już dla λ ≈ 1.0 - potrzebne adaptacyjne zagęszczenie siatki i bardziej stabilny model materiału.

Uwagi końcowe

1. W analizie efektów nieliniowych dominuje modelowanie 3D.

2. Konsystentna linearyzacja równań zapewnia kwadratową zbieżność procedury Newtona-Raphsona.

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

4. W projektowaniu akceptuje się zazwyczaj połączenie liniowo sprężystych obliczeń statycznych celem wyznaczenia naprężeń (sił przekrojowych) z analizą stanów granicznych uwzględniających uplastycznienie lub zarysowanie.

5. W obliczeniach nieliniowych szacuje się mnożnik obciążenia, przy którym następuje uszkodzenie/zniszczenie/utrata stateczności konstrukcji - ma on interpretację globalnego współczynnika

bezpieczeństwa, więc obliczenia powinno się prowadzić dla średnich wartości obciążeń i wytrzymałości.

MKwIL, Budownictwo II st.

Literatura

[1] R. de Borst, M.A. Crisfield, J.J.C. Remmers and C.V. Verhoosel. Non-linear Finite Element Analysis of Solids and Structures. Second Edition, J. Wiley & Sons, Chichester 2012.

[2] M. Radwańska. Ustroje powierzchniowe, podstawy teoretyczne oraz rozwiązania analityczne i numeryczne. Wydawnictwo PK, Kraków, 2009.

[3] M. Kwasek Advanced static analysis and design of reinforced concrete deep beams.

Diploma work, Politechnika Krakowska, 2004.

[4] M. Asin The behaviour of reinforced concrete continuous deep beams. PhD Thesis, Delft University of Technology, 2000.

[5] Z. Waszczyszyn, E. Pabisek, J. Pamin, M. Radwańska. Nonlinear analysis of a RC cooling tower with geometrical imperfections and a technological cut-out. Engineering Structures, 2, 480-489, 2000.

[6] A. Moroński. Analiza zarysowania i utraty stateczności uszkodzonej powłoki żelbetowej chłodni kominowej. Praca dyplomowa, Politechnika Krakowska, Kraków, 1996.

[7] A. Moroński, J. Pamin, M. Płachecki, Z. Waszczyszyn. Fracture and loss of stability of a partly-damaged cooling tower shell. Proc. 2nd Int. DIANA Conf. on Finite Elements in Engineering and Science, Eds M.A.N. Hendriks et al, 107-110, Kluwer Academic Publishers, Dordrecht, 1997.