ZAGADNIENIE DOBREGO POSTAWIENIA PROBLEMU OPISUJĄCEGO SPRĘŻYSTO-KRUCHY PROCES ROZWOJU USZKODZEŃ PRZY OBLICZENIACH METODĄ ELEMENTÓW SKOŃCZONYCH

77

ZAGADNIENIE DOBREGO POSTAWIENIA PROBLEMU OPISUJĄCEGO SPRĘŻYSTO- KRUCHY PROCES ROZWOJU USZKODZEŃ PRZY OBLICZENIACH

METODĄ ELEMENTÓW SKOŃCZONYCH

Piotr Mika

Instytut Technologii Informatycznych w Inżynierii Lądowej, Politechnika Krakowska e-mail: pm@l5.pk.edu.pl

Streszczenie

Jednym z warunków przeprowadzenia poprawnych obliczeń inżynierskich, wykonanych metodą elementów skończonych, jest niezależność otrzymanych wyników od przyjętej gęstości siatki modelującej konstrukcję.

Z matematycznego punktu widzenia błędne obliczenia mogą być spowodowane utratą eliptyczności równania konstytutywnego, co prowadzi do niepoprawnego postawienia problemu. Stosując standardowe procedury metody elementów skończonych, należy uzasadnić poprawne sformułowanie przyjętego modelu konstytutywnego.

W zamieszczonych analizach konstrukcji, wykonanych z materiałów podatnych na rozwój uszkodzeń, odpowiedź ośrodka jest opisana przez równanie konstytutywne materiału liniowo-sprężystego z uszkodzeniami, natomiast bieżący stan uszkodzonego materiału przez symetryczny tensor uszkodzeń rzędu drugiego. Prawo konstytutywne, równanie ewolucji uszkodzeń i kryterium zniszczenia jest wyprowadzone na bazie teorii reprezentacji izotropowych funkcji tensorowych.

ISSUE OF THE WELL-POSEDNESS PROBLEM DESCRIBING ELASTIC BRITTLE DAMAGE EVOLUTION PROCESS

AT THE CALCULATIONS WITH THE FINITE ELEMENT METHOD

Summary

The response of the material is modeled by the constitutive equation of the elastic medium with damages. As a result of the possible loss of the ellipticity for the constitutive equations, which leads to the ill-posedness of the problem, the standard finite element methods typically yield mesh-dependent results. It is necessary to check the well-posedness of the assumed model. A current state of the deteriorated material is described by a symmetric second order damage tensor. The constitutive law, damage evolution equation and failure criterion are formulated by employing the theory of isotropic tensor function representation. Numerical computations are performed with the ABAQUS programme by means of finite element method and user procedure UMAT.

1. WSTĘP

Rozwój metod obliczeniowych oraz możliwość włączenia do programów komputerowych złożonych modeli konstytutywnych pozwala na wykonywanie szeregu analiz, uwzględniających np. efekty lepkosprężyste. Szczególnie ważne jest to dla konstrukcji z materiałów o właściwościach zmiennych w trakcie procesu, gdzie należy określić czas pojawienia się pierwszych makrouszkodzeń czy czas upływający do naruszenia stanu nośności. Tematyka ta jest przedmiotem badań kontynualnej mechaniki uszkodzeń [1], a jej rozwoju dowodzą propozycje nowych modeli konstytutywnych zamieszczanych w literaturze oraz wzbogacanie biblioteki programów komputerowych o nowe sformułowania konstytutywne, [2]. Poprawność sformułowań ma zapewnić spełnienie postulatów termodynamiki materiałów i uwzględnienie wyników badań doświadczalnych. Modele te powinny umożliwiać opis anizotropii rozwoju defektów, jednostronności więzów w obszarach zarysowań czy złożonych, nieproporcjonalnych obciążeń. Ponadto ważne są możliwości efektywnego przeprowadzenia obliczeń numerycznych dla rzeczywistych konstrukcji o skomplikowanych kształtach.

Stosowany w pracy tensorowy model rozwoju uszkodzeń, poprzez wprowadzenie kryterium zniszczenia dla materiału z malejącą wytrzymałością w procesie ewolucji, wprowadza informację związaną z kierunkowością rozprzestrzeniania się zarysowań poprzez kierunki własne tensora uszkodzeń. W ramach proponowanych analiz jest możliwe ustalenie współzależności pomiędzy konfiguracją wymiarów geometrycznych, warunkami brzegowymi, sposobem obciążenia konstrukcji i mechanizmami rozprzestrzeniania uszkodzeń. Poprawność otrzymywanych rezultatów jest warunkowana zachowaniem eliptyczności równań przyjętego modelu konstytutywnego. W celu przeprowadzenia obliczeń numerycznych model konstytutywny włączony został w ramach procedury użytkownika User Material (UMAT) do programu ABAQUS [3].

2. MODEL MATEMATYCZNY I OBLICZENIA

2.1 MODEL MATEMATYCZNY I POPRAWNE POSTAWIENIE PROBLEMU

Model fizyczny został wyprowadzony na podstawie zależności pomiędzy tensorami naprężeń i odkształceń z zastosowaniem rozwinięcia wynikającego z teorii reprezentacji funkcji tensorowych. Zakładając liniowo- tensorową zależność, właściwości fizyczne materiału

z uszkodzeniami są modelowane przez równanie konstytutywne ośrodka liniowo-sprężystego (1), w którym macierz podatności Aijkl (2) jest izotropową funkcją tensorową niezmienników tensora uszkodzeń D, wartości głównej tensora uszkodzeń D1 oraz stałych materiałowych zależnych od temperatury – modułu Younga E i współczynnika Poissona ν _[4]:

(

)

ijkl

ij A D σ

ε = ⋅ (1)

( )

( ) (

)

1 1

jk il jl ik kl

ij ijkl

D δ D δ D δ D Eδ D 1 4

D

δ δ δ 2E δ

ν δ 1 Eδ A ν

+ + + +

+

+ + +

+

−

= (2)

W powyższym równaniu ij oznacza symbol Kroneckera. Zmienne w czasie parametry materiałowe, powodujące zmienną w czasie sztywność konstrukcji, odwzorowują ewolucję uszkodzeń w materiale, opisaną przez dwie równoważne, tensorowe miary uszkodzeń Murakami-Ohno J and D, [3], przy czym wartości główne tensora uszkodzeń J zawarte są w przedziale zamkniętym [0, 1], gdzie wartość 0 oznacza brak uszkodzeń, a 1 oznacza całkowite zniszczenie wiązań materialnych na kierunku prostopadłym do danej wartości głównej.

Wartości główne obu tensorów powiązane są równaniem:

(

)

i

i Ω 1 Ω

D = − i=1,2,3 (3)

Warunki początkowe są sformułowane dla materiału dziewiczego bez uszkodzeń i = Di = 0.

Równanie ewolucji uszkodzeń jest zależne od zmiennych w czasie niezmienników tensora naprężeń T, dewiatora naprężeń S, tensora uszkodzeń D oraz od kombinacji stałych materiałowych v, E i składowej tensora uszkodzeń Ω1, która jest również zmienna w trakcie procesu ewolucji

[ ]

1,2,3

= i ), H(σ σ

D) tr(T trS T 2E tr Ω 2E

ν 1 6E

2ν k 1 Ω

i i

2 2 T 2 2 1 i

t

⋅

 ⋅















 − +

∂ =

(4) Funkcja Heaviside’a H w równaniu (4) eliminuje wzrost uszkodzeń na kierunkach naprężeń ściskających, natomiast k jest dodatkową stałą materiałową.

Kryterium zniszczenia jest postulowane jako trójparametrowa funkcja tensorowa, zależna od krytycznej kombinacji składowych tensorów naprężeń i uszkodzeń

[

]

C ^{2 2 2 T}− ²_u= (5)

gdzie σu jest zależną od temperatury granicą wytrzymałości materiału nieuszkodzonego. Wektor stałych materiałowych C=[C1, C2, C3] definiuje

79 aktualną hiperpowierzchnię uszkodzenia. Stan zniszczenia utożsamia się z chwilą czasu, w której hiperpowierzchnia zniszczenia, na skutek kontrakcji wynikającej z ewolucji uszkodzeń, osiągnie wektor naprężeń, co prowadzi do makrouszkodzenia.

Składowe wektora C wyznacza się z liniowego układu równań algebraicznych, który jest otrzymany przez zastosowanie równania (5) do trzech niezależnych stanów naprężeń – dwóch jednoosiowych rozciągań we wzajemnie prostopadłych kierunkach naprężeń głównych i biosiowego rozciągania w tych kierunkach [4].

I C

W ^T= (6)

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

















−

−

−

−

−

−

−

−

−

=

1 1 2

1 2

1

1 2 1 2 2

1 2 2

1 2

1 2 1 2

1 2

1

Ω Ω 1 2 Ω

31 Ω 2

1 4

Ω ν Ω ν 1 Ω ν 31 Ω 2

ν 1

Ω Ω 1 Ω

31 Ω 2

1

W (7)

Przekształcając równanie (1) do postaci:

•

•

=H:ε

σ (8)

można sformułować konieczny i wystarczający warunek dla poprawnego postawienia problemu [5]:

[

]

det ⋅ ⋅ ≠ ∀ ∈ ³ ≠ (9)

Jeśli warunek (9) nie jest spełniony, równanie traci eliptyczność. Równanie (1) przy uwzględnieniu początkowej współosiowości kierunków własnych tensorów naprężeń i uszkodzeń, po zróżniczkowaniu i uwzględnieniu członów istotnych w dalszych rozważaniach, przyjmuje postać (w zapisie macierzowym):

σ E D I Ω E

ν σ 1 EItr

ε& ν & ¹ &



 



 + +

+

−

= (10)

Dla sformułowania (10) trudno jest analitycznie znaleźć formę odwrotną ( Rσσσσ=σσσσ(εεεε)R R zamiastR Rεεεε=εεεε(σσσ)σ R RRRRdlatego w celu dalszych obliczeń przyjęto to oznaczenie:

1 1D E +Ω E I

ν +

= 1 B

−



 



 (11)

co po przekształceniach prowadzi do następującej formy równania konstytutywnego:



 



 +

⋅

= Itrσ

E ε ν B

σ& & & (12)

Kinematyczne warunki zgodności Maxwella, dla lokalizacji w paśmie o normalnej n, przyjmują postać [5]:

(

)

2 ε 1 ε

ε& =&₁−&₀= o + o (13) 0

n σ n σ 0 n

σ&⋅ = ⇒&₁⋅ −&₀⋅ = (14)

Mnożąc obustronnie równanie (10) przez n, rozpisując warunki równowagi w paśmie i uwzględniając równania (13) i (14), można wyprowadzić warunek na utratę eliptyczności:

(

)

det (15)

następnie po rozpisaniu zależności (15) zgodnie z rachunkiem macierzowym (B^c – macierz kofaktorów):

( ) ( )

n) (B B n detB

n n) (B B tr detB

= Bn°n + B det

C c

⋅

⋅

⋅ +

=

=

⋅

⋅

+ o (16)

W głównych kierunkach macierz B posiada niezerowe współrzędne na kierunkach głównych, co po przekształceniach pozwala wykazać, że dla każdych wartości RΩ i D:

0 n) (B B n oraz 0

detB> ⋅ ^C⋅ ⋅ > (17) ponieważ poszczególne elementy macierzy B są wyrażone przez dodatnie ilorazy parametrów materiałowych - dlai 1,2,3

) D Ω ν (1

E

i 1

+ = +

2.2 PRZYKŁADY NUMERYCZNE

W analizie numerycznej przebadano szczegółowo lokalizację strefy uszkodzeń, która powinna być niezależna od gęstości podziału siatki MES, a otrzymane rozwiązania numeryczne powinny być zbieżne. Dokonano również oceny czasu pierwszych pęknięć w odniesieniu do czasu zniszczenia przekroju nośnego konstrukcji. Zamieszczone wizualizacje umożliwiają przewidywanie możliwych mechanizmów rozwoju frontu uszkodzeń.

Do wykonania poglądowych analiz wybrano przestrzenną konstrukcję (rys. 1) modelowaną kontynualnymi, ośmiowęzłowymi elementami 3D.

Warunki swobodnego podparcia przyjęto wzdłuż krótszych boków, obciążenie równomiernie rozłożone jest przyłożone do górnej powierzchni płaszczyzny w strefie znajdującej się nad otworem. Konstrukcja jest wykonana ze stali ANSI, pracującej w temperaturze 811 ⁰K o następujących parametrach materiałowych: wytrzymałość σu=288 MPa

417 E/σ

E= u= , ν=0.47, k=kσ³_uτ=82.1, gdzie τ=1[h]

jest czasem jednostkowym. Bieżący stan rozwoju uszkodzeń reprezentowany przez poszczególne współrzędne tensora oraz informacje o makrouszkodzeniach są przechowywane w zmiennych użytkownika włączonych do procedury user material (UMAT) programu ABAQUS.

Aby dowieść braku zależności rozwiązań od gęstości podziału siatki MES, przeprowadzono obliczenia dla różnych gęstości siatki, od siatki zgrubnej do bardzo gęstej zilustrowanej na rys. 2.

Obszar wokół otworu w środkowej części został dodatkowo zagęszczony.

Poglądowe mapy konturowe maksymalnej wartości głównej tensora uszkodzeń ΩΩΩΩ, ilustrujące przebieg procesu rozwoju mikro i makrouszkodzeń, w początkowym stadium procesu, zamieszczono na poniższych rysunkach (rys. 3 i rys. 4).

Poszczególne barwy odpowiadają stopniu zaawansowania procesu wzrostu mikrouszkodzeń, który jest proporcjonalny do wartości własnych tensora ΩΩΩΩ Najintensywniej proces rozwija się w elementach wokół otworu, ewoluując w kierunku

płaszczyzny górnej. Rozwój procesu jest również dość intensywny w pewnym oddaleniu od podpór (ok. 20%

odległości od krawędzi konstrukcji do krawędzi otworu), w płaszczyźnie dolnej.

Na przebieg procesu, oprócz gęstości energii sprężystej, zgodnie z równaniem ewolucji (4), mają wpływ główne naprężenia rozciągające, [6], co potwierdza ich rozkład zilustrowany na rys. 5.

Naprężenia te na kierunkach prostopadłych do otworu powodują intensywny rozwój uszkodzeń w tym

Rys. 1. Geometria konstrukcji i warunki podparcia Rys. 2. Siatka MES

Rys. 3. Ewolucja wartości głównej tensora uszkodzeń Ω1 w płaszczyźnie dolnej

Rys. 4. Ewolucja wartości głównej tensora uszkodzeń Ω1

w płaszczyźnie górnej

Rys. 5. Ewolucja naprężeń głównych σ

81 obszarze. Naprężenia ściskające w dolnej warstwie elementów przypodporowych nie powodują rozwoju mikrouszkodzeń.

Wartości główne Ω1 tensora uszkodzeń w chwili krytycznej, w płaszczyźnie górnej konstrukcji, zilustrowano na rys. 6. Wartości krytyczne składowych tensora uszkodzeń, przy których powstają makrouszkodzenia, są poniżej jedności, co jest zgodne z wynikami badań doświadczalnych. Analiza historii ewolucji dowodzi również, że kolejne makropęknięcia w elementach przekroju poprzecznego pojawiają się w krótkim interwale czasowym.

Rys. ilustruje makrouszkodzenia w elementach, zaznaczone czerwonym kolorem, powstałe w trakcie procesu ewolucji, w których spełniono kryterium zniszczenia (5).

W rejonach przypodporowych makrouszkodzenia powstają z powodu przekroczenia granicy wytrzymałości w obszarze naprężeń ściskających – są zaznaczone zielona barwą.

Rys. . ilustruje zmianę ugięć, odniesionych do wartości początkowej, w obszarach o maksymalnych wartościach przemieszczeń w zależności od gęstości podziału siatki MES. Chwila asymptotycznego wzrostu funkcji ugięcia wskazuje na moment krytyczny dla konstrukcji. Porównywalne jakościowo rozwiązania dla wszystkich siatek dowodzą poprawnego postawienia problemu, co potwierdza rozważania teoretyczne. Wykres dowodzi również poprawiającej się zbieżności wyników w miarę zagęszczania siatki MES.

Rys. 6. Ewolucja wartości Ω1 w płaszczyźnie górnej, w końcowym stadium procesu

Rys. 7. Makropęknięcia w elementach konstrukcji

3. PODSUMOWANIE

Zastosowany model konstytutywny umożliwia efektywne wykonywanie obliczeń dla konstrukcji, w których występuje proces rozwoju uszkodzenia.

Włączenie modelu materiałowego do komercyjnego pakietu MES, wyposażonego w pre- i post procesor, ułatwia wykonanie i wizualizację wyników numerycznych. Analizy teoretyczne i przykład numeryczny dowodzą poprawnego postawienia modelu

konstytutywnego. Ponadto jest możliwe, co ważne z inżynierskiego punktu widzenia, określenie czasu pierwszych mikropęknięć i symulowanie procesu rozwoju uszkodzeń, ze wskazaniem kierunku przemieszczania się frontu zniszczenia. Obliczenia również pozwalają określić, czy proces będzie przebiegał w sposób lawinowy po pojawieniu się pierwszych uszkodzeń, a jeśli nie, to określić czas bezpiecznej pracy konstrukcji.

Literatura

1. Krajcinovic D.: Damage mechanics. Elsevier Science, 1996

2. ABAQUS manuals, ver. 6.9. Dassault Systèmes Simulia Corp., Providence, RI, USA.

3. Murakami S., Ohno N.: A continuum theory of creep and creep damage. In: Creep in Structures, Eds.: A. R. S.

Ponter and D. R. Hayhurst. Berlin: Springer, 1981, p. 422-444.

4. Litewka A.: Effective material constants for orthotropically damaged elastic solid. Arch. Mech.1985, 37, p. 631- 642.

5. Benallal A.: A note on ill-posedness for rate-dependent problems and its relation to the rate-independent case.

Comput. Mech., 2008, 42.2, p. 261-269.

6. Mika, P.: Modelowanie konstrukcji powłokowej z uwzględnieniem procesu rozwoju uszkodzeń. Zesz. Pol. Rzesz.

„Budownictwo i Inżynieria Środowiska” 2011, nr 3, z. 58, II, s.405-412.

Badania zrealizowano w ramach projektu „Politechnika XXI wieku - Program rozwojowy Politechniki Krakowskiej;

najwyższej jakości dydaktyka dla przyszłych polskich inżynierów”; współfinansowanego ze środków Unii Europejskiej w ramach Europejskiego Funduszu Społecznego; umowa nr UDA-POKL.04.01.01-00-029/10-00

Proszę cytować ten artykuł jako:

Mika P.: Zagadnienie dobrego postawienia problemu opisującego sprężysto-kruchy proces rozwoju uszkodzeń przy obliczeniach Metodą Elementów Skończonych. „Modelowanie Inżynierskie” 2013, nr 46, t. 15, s. 77 – 82.

Rys. 8. Historia unormowanych ugięć w wybranych elementach dla różnych gęstości siaki MES