MODELOWANIE INŻYNIERSKIE ISSN 1896-771X 35, s. 23-30, Gliwice 2008
MODELOWANIE KOMPUTEROWE PRÓB PĘKANIA PRZY OBCIĄŻENIU DYNAMICZNYM
PIOTR FEDELIŃSKI
Katedra Wytrzymałości Materiałów i Metod Komputerowych Mechaniki, Politechnika Śląska e-mail: piotr.fedelinski@polsl.pl
Streszczenie. W pracy przedstawiono sformułowanie w dziedzinie czasu metody elementów brzegowych (MEB) zastosowane do modelowania prób pękania przy obciążeniu dynamicznym. Opracowano program komputerowy, który wykorzystano do wyznaczenia kształtu wzrastającego pęknięcia i dynamicznych współczynników intensywności naprężeń w specjalnej próbce obciążonej za pomocą dzielonego pręta Hopkinsona. W obliczeniach numerycznych wykorzystano wyznaczone doświadczalnie obciążenie dynamiczne próbki i prędkość wzrostu pęknięcia. Porównano kształty pęknięcia określone metodą numeryczną i doświadczalną.
1. WPROWADZENIE
Celem prób pękania przy obciążeniu dynamicznym jest wyznaczenie odporności materiałów na pękanie, a także określenie zależności od czasu dynamicznych współczynników intensywności naprężeń (DWIN), kierunków, prędkości wzrostu i kształtu pęknięć. Pomiary umożliwiają określenie praw opisujących zależności kierunków i prędkości wzrostu od DWIN i innych parametrów. Ze względu na swoją złożoność badania eksperymentalne modelowane są różnymi metodami komputerowymi: metodą różnic skończonych (MRS), metodą elementów skończonych (MES), metodą elementów brzegowych (MEB) i metodami bezsiatkowymi (MB). W ostatnich latach metody doświadczalne i komputerowe często łączone są ze sobą.
Klasyfikację takich metod hybrydowych podali Nishioka i inni [1] i [2]. Zależne od czasu warunki brzegowe (przemieszczenia i siły powierzchniowe), które są konieczne w obliczeniach numerycznych, wyznacza się metodami doświadczalnymi. Programy komputerowe określają kierunki i prędkości wzrostu pęknięć na podstawie praw pękania. Przyjęte prawa można zweryfikować poprzez porównanie wyznaczonych numerycznie i eksperymentalnie DWIN, kształtów pęknięć i prędkości ich wzrostu.
Najczęściej stosowaną metodą analizy pęknięć jest metoda elementów skończonych. Bui, Maigre i Rittel [3], [4], [5] i [6] wykorzystali MES do analizy specjalnej próbki z pęknięciem o stałej długości, obciążonej dzielonym prętem Hopkinsona. Wyznaczone numerycznie DWIN i prędkości obciążonych krawędzi porównano z wynikami badań doświadczalnych. Weisbrod i Rittel [7] analizowali jednopunktowo zginaną krótką próbkę obciążoną prętem z czujnikami odkształceń. Porównano DWIN wyznaczone za pomocą MES z wynikami doświadczalnymi
do chwili wzrostu pęknięcia. Nishioka i inni [2] analizowali za pomocą MES trójpunktowo zginaną próbkę obciążoną spadającym młotem udarowym. Porównano wyznaczone MES kształty pęknięć i DWIN z wynikami doświadczalnymi dla różnych miejsc uderzenia młota.
Metoda elementów brzegowych jest szczególnie odpowiednia dla analizy wzrostu pęknięć.
Sposób modelowania jest prostszy niż w metodzie różnic skończonych lub elementów skończonych, ponieważ dyskretyzuje się wyłącznie powierzchnie zewnętrzne ciała i powierzchnie pęknięć. Spośród różnych wariantów MEB (Dominguez [8]) do analizy dynamicznie wzrastających pęknięć stosowane jest sformułowanie w dziedzinie czasu.
Sformułowanie metody dla wzrostu pęknięcia ze zmienną prędkością z uwzględnieniem kontaktu powierzchni pęknięcia przedstawili Sellig i Gross [9] i [10]. Metodę zastosowali Sellig, Gross i Pothmann [11] do analizy tej samej próbki, którą wcześniej analizowali Bui, Maigre i Rittel [3]. Uwzględniono dodatkowo dynamiczny wzrost pęknięcia. Porównano wyniki wyznaczone MEB, MES i doświadczalne.
Pierwsze sformułowanie MEB, w którym otrzymano rozwiązanie numeryczne dla dynamiki pęknięć o stałej długości za pomocą MEB i dyskretyzacji wyłącznie brzegów układu przedstawili Fedeliński, Aliabadi i Rooke [12]. Ci sami autorzy w pracy [13] przedstawili po raz pierwszy rozwiązanie numeryczne za pomocą MEB dla wzrastającego pęknięcia, w którym program komputerowy wyznaczał kierunek wzrostu na podstawie kryterium pękania. Różne praktyczne zastosowania metody w mechanice pękania przedstawiono w pracach Fedelińskiego [14] i [15].
Celem pracy jest krótkie omówienie sformułowania dualnego w dziedzinie czasu MEB w dynamice pęknięć. Metodę zastosowano do rozwiązania nowego przykładu numerycznego - specjalnej próbki z pęknięciem obciążonej dzielonym prętem Hopkinsona. Badania doświadczalne i modelowanie komputerowe tej próbki za pomocą rozszerzonej metody elementów skończonych (ang. eXtended Finite Element Method X-FEM) przedstawiono w pracach Gregorie [16] i Combescure [17] i innych. W obliczeniach numerycznych za pomocą MEB wykorzystano wyznaczone doświadczalnie obciążenie dynamiczne próbki i prędkość wzrostu pęknięcia. Wyznaczono zmienność DWIN i odkształcenia próbki w czasie wzrostu pęknięcia. Porównano kształty pęknięcia określone metodą numeryczną i doświadczalną.
2. METODA ELEMENTÓW BRZEGOWYCH W MODELOWANIU WZROSTU PĘKNIĘĆ PRZY OBCIĄŻENIU DYNAMICZNYM
Dla kompletności pracy zostanie krótko omówione sformułowanie metody elementów brzegowych w modelowaniu wzrostu pęknięć przy obciążeniu dynamicznym. Dokładniejszy opis metody znajduje się w pracach Fedelińskiego [13] i [14].
2.1. Brzegowe równania całkowe dla ciała z pęknięciem
Metoda będzie stosowana dla ciał liniowosprężystych, jednorodnych i izotropowych z szybko wzrastającymi pęknięciami. Brzeg ciała, oznaczony przez Γ(t), jest funkcją czasu t z powodu wzrostu pęknięcia. Dla ciała obciążonego tylko siłami powierzchniowymi i dla zerowych warunków początkowych, przemieszczenie punktu x’ może być określone za pomocą następującego brzegowego równania całkowego
0 0
( ') ( ', ) [ ( ', , , ) ( , ) ( )] [ ( ', , , ) ( , ) ( )]
t t
ij j ij j ij j
c x u x t U x t xτ t xτ d x dτ T x t xτ u xτ d x dτ
Γ Γ
=
∫ ∫
Γ −∫ ∫
Γ , (1)gdzie Uij(x’,t,x,τ) i Tij(x’,t,x,τ) są rozwiązaniami fundamentalnymi elastodynamiki, uj(x,τ) i tj(x,τ) są odpowiednio przemieszczeniami brzegowymi i siłami powierzchniowymi, cij(x’) jest stałą, która zależy od położenia punktu kolokacji, x’ jest punktem kolokacji, x jest punktem brzegowym, a t jest czasem, w którym określa się przemieszczenie.
W prezentowanym sformułowaniu metody elementów brzegowych, które nazywa się metodą dualną, stosowane są jednocześnie brzegowe równania całkowe przemieszczeń i sił powierzchniowych dla punktów na powierzchniach pęknięcia. Brzegowe równanie całkowe przemieszczeń dla pokrywających się punktów x’ i x” na przeciwległych gładkich powierzchniach pęknięcia ma formę
0
0
1 1
( ', ) ( ", ) [ ( ', , , ) ( , ) ( )]
2 2
[ ( ', , , ) ( , ) ( )]
t
i i ij j
t
ij j
u x t u x t U x t x t x d x d
T x t x u x d x d
τ τ τ
τ τ τ
Γ
Γ
+ = Γ
− Γ
∫ ∫
∫ ∫
. (2)
Brzegowe równanie całkowe sił powierzchniowych dla tych samych punktów ma postać
0
0
1 1
( ', ) ( ", ) ( '){ [ ( ', , , ) ( , ) ( )]
2 2
[ ( ', , , ) ( , ) ( )] }
t
j j i kij k
t
kij k
t x t t x t n x U x t x t x d x d
T x t x u x d x d
τ τ τ
τ τ τ
Γ
Γ
− = Γ
− Γ
∫ ∫
∫ ∫
, (3)
gdzie Ukij(x’,t,x,τ) i Tkij(x’,t,x,τ) są rozwiązaniami fundamentalnymi elastodynamiki, a ni(x’) jest jednostkowym wektorem normalnym do powierzchni pęknięcia, zwróconym na zewnątrz materiału, w punkcie kolokacji.
2.2. Numeryczna realizacja metody
Modelowanie numeryczne układów obciążonych dynamicznie wymaga dyskretyzacji w przestrzeni i czasie. Brzeg ciała Γ(τ) podzielono na kwadratowe elementy brzegowe, a czas analizy t na N kroków czasowych. Przemieszczenia i obciążenia brzegowe interpolowano elementami ciągłymi na powierzchniach nienależących do pęknięcia, półciągłymi w miejscu połączenia brzegu zewnętrznego z pęknięciami krawędziowymi i prostoliniowymi elementami nieciągłymi na powierzchniach pęknięcia. Geometrię brzegu interpolowano elementami ciągłymi. W każdym kroku czasowym przemieszczenia interpolowano funkcjami liniowymi, a obciążenia funkcjami przedziałami stałymi. W metodzie wykorzystuje się równanie przemieszczeń (1) dla węzłów, które nie należą do powierzchni pęknięcia, równanie przemieszczeń (2) i sił powierzchniowych (3), które stosuje się jednocześnie dla węzłów na powierzchniach pęknięcia. Całki ze względu na czas, dla prostych funkcji interpolujących, całkowano analitycznie. W wyniku dyskretyzacji i całkowania otrzymuje się układ równań algebraicznych dla kroku czasowego N, który może być zapisany w następującej postaci macierzowej
1
1
( )
−
=
= +
∑
N −NN N NN N Nn n Nn n
n
H u G t G t H u , (4)
gdzie un i tn zawierają wartości węzłowe przemieszczeń i sił powierzchniowych w kroku czasowym n, HNn i GNn zależą od rozwiązań fundamentalnych elastodynamiki i funkcji interpolujących. Równanie macierzowe rozwiązuje się krokowo, żeby wyznaczyć nieznane przemieszczenia i obciążenia na brzegu. W każdym kroku czasowym oblicza się dwie nowe macierze HNn i GNn i przechowuje w celu obliczeń w następnych krokach. Dla wzrastającego pęknięcia macierze obliczone w poprzednich krokach czasowych zwiększa się poprzez dodanie nowych podmacierzy, które odpowiadają nowym punktom kolokacji i elementom dodanym w czasie ostatniego wzrostu pęknięcia. Proces rozwiązywania staje się stopniowo coraz powolniejszy na skutek zwiększającej się liczby punktów kolokacji, elementów brzegowych i koniecznych modyfikacji wszystkich macierzy obliczonych w poprzednich krokach czasowych.
2.3. Modelowanie dynamicznego wzrostu pęknięcia
Dynamiczne współczynniki intensywności naprężeń (DWIN) wyznaczono na podstawie względnych przemieszczeń powierzchni pęknięcia
( )
(
22)
21
2 2 2 2
1
1 4
1 4
2 2 u
KI r ∆
− +
= −
β β
β β
β
µ π , (5)
( )
(
22)
12
2 2 2 2
1
1 4
1 4
2 2 u
KII r ∆
− +
= −
β β
β β
π β
µ , (6)
gdzie:
2 1 1= 1−(c c )
β , β2 = 1−(c c2)2 , (7)
gdzie µ jest modułem sprężystości poprzecznej, ∆u1 i ∆u2 są względnymi przemieszczeniami w kierunku stycznym i prostopadłym do powierzchni pęknięcia na przeciwległych powierzchniach pęknięcia, r jest odległością tych punktów od wierzchołka pęknięcia, c jest prędkością wzrostu pęknięcia, c1 i c2 są odpowiednio prędkościami fali podłużnej i poprzecznej. Rozkład naprężeń w pobliżu wierzchołka pęknięcia jest obliczany na podstawie DWIN i prędkości wzrostu pęknięcia. Opracowana metoda może być stosowana do analizy wzrostu pęknięcia ze zmienną prędkością, gdy kształt pęknięcia nie jest wcześniej określony.
Założono, że pęknięcie będzie wzrastało w kierunku określonym przez maksymalne naprężenie obwodowe w otoczeniu wierzchołka pęknięcia. Wzrost pęknięcia modelowano poprzez dodawanie w wierzchołku pęknięcia nowych prostoliniowych elementów brzegowych o długości
t c a= ∆
∆ , (8)
gdzie ∆t jest długością kroku czasowego.
3. PRZYKŁAD NUMERYCZNY
Metodę zastosowano do rozwiązania nowego przykładu numerycznego - specjalnej próbki z pęknięciem obciążonej dzielonym prętem Hopkinsona. Badania doświadczalne
i modelowanie komputerowe tej próbki za pomocą rozszerzonej metody elementów skończonych (ang. eXtended Finite Element Method X-FEM) przedstawiono w pracach Gregorie [16] i Combescure [17] i innych. Dane potrzebne do przeprowadzenia symulacji komputerowej zaczerpnięto z pracy [16]. Specjalna próbka została umieszczona pomiędzy dzielonym prętem Hopkinsona składającym się z pręta 1 i 2 (rys.1). W pręt 1 uderza bijak.
Prędkość bijaka jest mierzona przez czujnik optyczny. Odkształcenia prętów są mierzone przez czujniki tensometryczne 1 i 2 naklejone na pręt 1 i czujnik 3 naklejony na pręt 2. Próbka jest oświetlona i 4 kamery rejestrują stan próbki w czasie próby. Bijak uderza w pręt 1 powodując propagację fali przez pręt 1, próbkę i pręt 2. Wymiary próbki w milimetrach podano na rys. 2.
Przesunięcie pęknięcia w stosunku do osi próbki powoduje jednoczesne rozciąganie i ścinanie wzdłużne pęknięcia przy obustronnym ściskaniu próbki. Próbka wykonana jest z polimetakrylanu metylu o następujących własnościach: gęstość ρ=1180 kg/m3, moduł Younga E=3.3 GPa, współczynnik Poissona ν=0.42 i układ znajduje się w płaskim stanie odkształcenia. Modelowano tylko próbkę, którą dyskretyzowano 70 elementami brzegowymi.
W czasie wzrostu pęknięcia w każdym kroku czasowym dodawano 2 elementy brzegowe w wierzchołku pęknięcia. Końcowa liczba elementów brzegowych wynosiła 130. Krok czasowy był równy ∆t= 10 µs. Lewą i prawą krawędź obciążono siłami równomiernie rozłożonymi t1 i t2, które zarejestrowano doświadczalnie w czasie próby. Zmienność sił wypadkowych na krawędziach F1 i F2 przedstawiono na rys. 3. Pęknięcie było stacjonarne od momentu obciążenia próbki do 200 µs, następnie wzrastało z prędkością 210 m/s, w czasie od 270 do 320 µs nastąpiło zatrzymanie wzrostu, a następnie ponowny wzrost z prędkością 160 m/s. Obliczenia wykonano dla 500 µs.
Rys. 1. Układ pomiarowy
Rys. 2. Wymiary próbki w milimetrach i obciążenie krawędzi
Zmienność w czasie dynamicznych współczynników intensywności naprężeń (DWIN) KI
i KII przedstawiono na rys. 5. Po dotarciu fali podłużnej do wierzchołka pęknięcia następuje wzrost DWIN. W chwili kiedy rozpoczyna się wzrost pęknięcia zmniejszają się DWIN. Kiedy współczynnik KI ma minimalną wartość, tzn. w czasie od 270 do 320 µs, wówczas następuje zatrzymanie wzrostu pęknięcia. Przy ponownym zwiększaniu się KI następuje dalszy wzrost pęknięcia. Współczynnik KII ma bardzo małe wartości w czasie wzrostu pęknięcia. Gregorie i inni [16] w symulacji komputerowej tej próby przyjęli krytyczną wartość DWIN KIC=1.33 MPam1/2. Z rys. 4 wynika, że kiedy następuje zatrzymanie wzrostu pęknięcia współczynnik KI
ma wartość mniejszą od krytycznej.
Rys.3. Zmienność w czasie sił wypadkowych obciążających krawędzie próbki
Rys.4. Zmienność w czasie dynamicznych współczynników intensywności naprężeń Próbkę z końcowym pęknięciem przedstawiono na rys. 5. Kształt pęknięcia wyznaczony numerycznie porównano z kształtem określonym doświadczalnie [16] na rys. 6. Widoczna jest bardzo dobra zgodność wyników
Rys. 5. Próbka z końcowym pęknięciem
Rys.6. Porównanie kształtów pęknięć wyznaczonych numerycznie i doświadczalni
Na rys. 7 przedstawiono odkształconą próbkę w różnych fazach wzrostu pęknięcia: rys. 7a – dla 200 µs, kiedy rozpoczyna się wzrost pęknięcia, rys. 7b – dla 300 µs, w czasie zatrzymania wzrostu pęknięcia, rys. 7c – dla 400 µs, przy ponownym wzroście pęknięcia i rys.
7d – dla 500 µs, dla końca analizy numerycznej.
a) b)
c) d)
Rys. 7. Kształt próbki z pęknięciem w różnych chwilach czasowych:
a) 200 µs, b) 300 µs, c) 400 µs, d) 500 µs 4. PODSUMOWANIE
W pracy przedstawiono zastosowanie sformułowania w dziedzinie czasu metody elementów brzegowych do modelowania próby pękania przy obciążeniu dynamicznym. Metoda umożliwia określenie zmienności w czasie przemieszczeń układu, kształtu wzrastającego pęknięcia i dynamicznych współczynników intensywności naprężeń. Porównanie kształtu pęknięcia wyznaczonego numerycznie i doświadczalnie wskazuje na wysoką dokładność metody elementów brzegowych.
LITERATURA
1. Nishioka T.: Hybrid numerical methods in static and dynamic fracture mechanics. „Optics and Lasers in Engineering” 1999, 32, s. 205-255.
2. Nishioka T., Tokudome H., Kinoshita M.: Dynamic fracture-path prediction in impact fracture phenomena using moving finite element method based on Delaunay automatic mesh generation. „International Journal of Solids and Structures” 2001, 38, s. 5273-5301.
3. Bui H.D., Maigre H., Rittel D.: A new approach to the experimental determination of the dynamic stress intensity factors. „International Journal of Solids and Structures” 1992, 29, s. 2881-2895.
4. Maigre H., Rittel D.: Mixed-mode quantification for dynamic fracture initiation:
Application to the compact compression specimen. „International Journal of Solids and Structures” 1993, 30, s. 3233-3244.
5. Maigre H., Rittel D.: Dynamic fracture detection using the force-displacement reciprocity:
application to the compact compression specimen. „International Journal of Fracture”
1995, 73, s. 67-79.
6. Rittel D., Maigre H.: An investigation of dynamic crack initiation in PMMA. „Mechanics of Materials” 1996, 23, s. 229-239.
7. Weisbrod G., Rittel D.: A method for dynamic fracture toughness determination using short beams. „International Journal of Fracture” 2000, s. 89-103.
8. Dominguez J.: Boundary elements in dynamics. Computational Mechanics Publications, Southampton, 1993.
9. Sellig Th., Gross D.: Analysis of dynamic crack propagation using a time-domain boundary integral equation method. „International Journal of Solids Structures” 1997, 34, s. 2087-2103.
10. Seelig Th., Gross D.: On the stress wave induced curving of fast running cracks – a numerical study by a time-domain boundary element method. „Acta Mechanica” 1999, 132, s. 47-61.
11. Seelig Th., Gross D., Pothmann K.: Numerical simulation of a mixed-mode dynamic fracture experiment. „International Journal of Fracture” 1999, 99, s. 325-338.
12. Fedeliński P., Aliabadi M.H., Rooke D.P.: A single-region time-domain BEM for dynamic crack problems. „International Journal of Solids and Structures” 1995, 32, s. 3555-3571.
13. Fedeliński P., Aliabadi M.H., Rooke D.P.: The time-domain DBEM for rapidly growing cracks. „International Journal for Numerical Methods in Engineering” 1997, 40, s. 1555- 1572.
14. Fedeliński P.: Metoda elementów brzegowych w analizie dynamicznej układów odkształcalnych z pęknięciami. Zeszyty Naukowe Politechniki Śląskiej, Mechanika, 137, Gliwice, 2000.
15. Fedeliński P.: Boundary element method in dynamic analysis of cracks, „Engineering Analysis with Boundary Elements” 2004, 28, s. 1135-1147.
16. Gregorie D., Maigre H., Rethore J., Combescure A.: Dynamic crack propagation under mixed-mode loading – comparison between experiments and X-FEM simulations.
„International Journal of Solids and Structures” 2007, 44, s. 6517-6534.
17. Combescure A., Gravouil A., Gregorie D., Rethore J.: X-FEM a good candidate for energy conservation in simulation of brittle dynamic crack propagation. „Computer Methods in Applied Mechanics and Engineering” 2008, 197, , s. 309-318.
COMPUTER MODELLING OF FRACTURE TESTS UNDER DYNAMIC LOADING
Summary. In this work the time-domain formulation of the boundary element method (BEM) is presented and applied for modelling fracture tests under dynamic loading. A computer code is developed and applied to determine the crack shape and dynamic stress intensity factors in a special specimen loaded by the split Hopkinson bar. In numerical computations the dynamic loading and crack growth velocity determined experimentally are used. Numerical and experimental crack shapes are compared.
