MODELOWANIE OSOBLIWYCH PÓL NAPRĘŻEŃ W ZAGADNIENIACH MECHANIKI KRUCHEGO PĘKANIA Z WYKORZYSTANIEM ANALITYCZNYCH ELEMENTÓW SKOŃCZONYCH

33, s. 5-10, Gliwice 2007

MODELOWANIE OSOBLIWYCH PÓL NAPRĘŻEŃ W ZAGADNIENIACH MECHANIKI KRUCHEGO PĘKANIA Z WYKORZYSTANIEM

ANALITYCZNYCH ELEMENTÓW SKOŃCZONYCH

ADAM ADAMOWICZ, ANDRZEJ SEWERYN

Katedra Mechaniki i Informatyki Stosowanej, Wydział Mechaniczny, Politechnika Białostocka e-mail: adamow@pb.edu.pl, seweryn@pb.edu.pl

Streszczenie. W pracy przedstawiono metodę elementów analitycznych, służącą do wyznaczania wartości parametrów opisujących osobliwe pola naprężeń w pobliżu wierzchołków ostrych naroży. Wykorzystano ją do wyznaczenia współczynników intensywności naprężeń oraz współczynników stojących przy członach wyższych rzędów rozwinięcia asymptotycznego opisującego pole naprężeń w pobliżu wierzchołka szczeliny. Otrzymane wyniki posłużyły do wyznaczenia warunków krytycznych propagacji szczelin z zastosowaniem nielokalnego naprężeniowego kryterium kruchego pękania.

1. WSTĘP

Identyfikacja pól naprężeń w pobliżu wierzchołka szczeliny jest nieodzownym elementem prognozowania pękania. Sprowadza się ona do wyznaczenia wartości parametrów analitycznych opisujących te pola. Wykorzystywane do tego techniki obliczeniowe podzielić można na trzy grupy [1]. W metodach asymptotycznych poszukiwane parametry analityczne wyznaczane są na podstawie porównania wyników obliczeń, np. metody elementów skończonych z rozkładami teoretycznymi [2]. W metodach energetycznych wartości poszukiwanych parametrów wyznaczane są na podstawie zmian energii potencjalnej układu, wywołanej zmianą wymiaru szczeliny [3], lub wykorzystania twierdzenia o wzajemności prac [4]. Prezentowana w pracy metoda elementów analitycznych [5], obok elementów hybrydowych [6] i metody więzów analitycznych [7], należy do trzeciej grupy - metod bezpośrednich. Poszukiwane parametry analityczne znajdują się w wektorze niewiadomych zagadnienia metody elementów skończonych, obok składowych przemieszczeń węzłów elementów skończonych, i wyznaczane są z układu równań równowagi sił węzłowych.

2. ROZKŁAD PÓL NAPRĘŻEŃ W OTOCZENIU WIERZCHOŁKA SZCZELINY

Opis osobliwego pola naprężeń i przemieszczeń w pobliżu wierzchołka szczeliny, w zagadnieniach liniowej teorii sprężystości, może być przedstawiony w układzie współrzędnych biegunowych (r, ϑ) (rys. 1) w postaci następującego rozwinięcia [8]:

( ) ( )

( ) ( )

I II

I II

1 1

I II

2 2

I II

1 1

I II

2 2

I II C

1 1

,

,

ϑ ϑ

ϑ ϑ

− −

= =

= =

= +

= + +

∑ ∑

∑ ∑

k k

m m

k k

ij k ij k ij

k k

k k

m m

k k

i k i k i i

k k

σ K r f K r f

u K r g K r g u

(1)

gdzie KIn, KIIn – współczynniki przy n-tym członie rozwiązania asymptotycznego, a KI1 i KII1

odpowiadają klasycznym współczynnikom intensywności naprężeń KI i KII:

( )

I II

0 0

i lim 2 _ϑϑ i _ϑ

ϑ σ τ

= → +

 

+ = _r  π + ^r 

K K r (2)

Wyrażenia opisane wzorami (1) zostaną zapisane w postaci macierzowej:

= K

σ fu , u=gu _K (3)

gdzie σ i u są wektorami kolumnami zawierającymi odpowiednio: składowe tensora naprężeń i składowe przemieszczeń, f i g są macierzami funkcji współrzędnych, a uK jest wektorem parametrów analitycznych, opisujących osobliwe pole naprężeń i przemieszczeń:

{

}

T

K = K KI, I2,KI3,...,KIm ,KII,KII2,KII3,...,KIIm ,uCr,uC_ϑ

u (4)

uCr i uCϑ są składowymi przemieszczenia wierzchołka szczeliny.

2l y

x r

Rys.1. Szczelina z biegunowym układem

współrzędnych (r, ϑ) Rys.2. Ciało ze szczeliną oraz obszary ΩA

i ΩO

3. METODA ELEMENTÓW ANALITYCZNYCH

Idea metody elementów analitycznych polega na wykorzystaniu specjalnego elementu skończonego do modelowania obszaru przywierzchołkowego szczeliny lub ostrego naroża.

W wektorze kolumnie parametrów węzłowych q zagadnienia MES wydzielono podwektor uO, zawierający składowe przemieszczeń węzłów standardowych elementów skończonych (obszar ΩO – rys. 2), podwektor uA składowych przemieszczeń węzłów leżących na granicy obszaru analitycznego ΩA oraz uK – wektor poszukiwanych parametrów analitycznych (4):

{ }

T T T T

O , A , K

=

q u u u (5)

Energia odkształcenia sprężystego rozpatrywanego układu przyjmie postać:

{

}

A K

1 , ,

2

   

  

=        U

K K 0 u

u u u K K 0 u

0 0 K u

(6)

gdzie KA jest macierzą sztywności elementu analitycznego ΩA.

Postać macierzy sztywności elementu analitycznego KA wyprowadzona została z wyrażenia na energię odkształcenia sprężystego obszaru analitycznego ΩA przy wykorzystaniu związków

pomiędzy wektorami kolumnami odkształceń ε i naprężeń σ oraz przy zastosowaniu opisu pola naprężeń za pomocą parametrów analitycznych (3):

Α A A

T T 1 T T 1 T

A K K K A K

1 1 1 1

d d d

2_Ω Ω 2 ⁻ 2  ⁻  2

=

∫

∫

∫

Ω Ω

U σ εt σ C σt Ω u f C ft Ω u u K u (7)

gdzie t – grubość elementu, C – macierz stałych sprężystych, f - macierz funkcji współrzędnych (3). Stąd macierz sztywności elementu analitycznego:

A

T 1

A =

∫

Ω

t Ω

K f C f (8)

Ciągłość pola przemieszczeń i naprężeń pomiędzy obszarem analitycznym ΩA a obszarem modelowanym za pomocą standardowych elementów skończonych ΩO zapewniono przez zastosowanie metody więzów analitycznych [5], [7]. Wykorzystano wymuszenie przemieszczeń węzłów granicznych obu obszarów zgodnie z zadanym rozkładem teoretycznym, uzależniając je od poszukiwanych wartości parametrów analitycznych:

O

O A

K K

0 0 0

   

  =  

      

   

   

u I

u Ψ u

u I u

u_A =Ψ u _K (9)

Postać macierzy więzów analitycznych Ψ wynika bezpośrednio z przyjętego opisu osobliwego pola przemieszczeń (1). Wprowadzenie macierzy więzów analitycznych spowoduje modyfikację globalnej macierzy sztywności. Ostatecznie energia odkształcenia sprężystego układu przybierze postać:

{

}

K

1 ,

2

∗ 

=  

U  u

u u K

u (10)

gdzie K^* jest zmodyfikowaną macierzą sztywności układu:

OO OA

*

T T

AO AA A

 

=  + 

K K Ψ

K Ψ K Ψ K Ψ K (11)

4. PROGNOZOWANIE PĘKANIA ELEMENTÓW ZE SZCZELINAMI

W nielokalnym naprężeniowym kryterium pękania, zaproponowanym przez Seweryna i Mroza [9], zakłada się, że inicjacja lub propagacja szczeliny następuje wówczas, gdy uśredniona na odcinku d0 długości strefy pękania funkcja naprężeń normalnych σn i tnących τn

w płaszczyźnie fizycznej osiągnie wartość krytyczną, czyli:

( )

0

( , ) 0 0

max 1 , d 1

ϑ  _σ σ τ 

 =

 

_d ^d

∫

x0 (12)

gdzie, x0 – początek lokalnego układu współrzędnych (r, ϑ) określający miejsce pękania. Do przewidywania pękania materiałów kruchych (takich jak polimetakrylan metylu) z wykorzystaniem kryterium (12), wystarczające jest przyjęcie lokalnej funkcji pękania w postaci warunku naprężeń normalnych [1]:

( )

σ _n = _n

R (13)

Długość strefy pękania d0 wyznacza się z równoważności kryterium Griffitha – Irwina (K_I =K ) oraz nielokalnego kryterium pękania (12) dla przypadku rozrywanej szczeliny [1]: _IC

2 IC 0

C

2 σ

 

=   π  

d K (14)

gdzie σC, KIC – naprężenia niszczące i krytyczna wartość współczynnika intensywności naprężeń.

Nielokalne kryterium kruchego pękania posłużyło do wyznaczenia warunków krytycznych pękania dla zagadnienia tarczy ze szczeliną nachyloną pod różnymi kątami do kierunku działania obciążenia (rRys.3), a wyniki porównano z badaniami doświadczalnymi [10].

0 30 60 90

0.0 0.4 0.8 1.2 1.6 2.0

0 30 60 90

-1.0 -0.5 0.0 0.5 1.0

0 30 60 90

0.0 0.1 0.2 0.3 0.4 m = 1

m = 2 m = 3

2l

Rys.3. Tarcza ze szczeliną ukośną oraz wartości współczynników intensywności naprężeń i współczynników stojących przy członach wyższych rzędów rozwinięcia asymptotycznego

w funkcji kąta pochylenia szczeliny γ dla szczeliny o długości l = 12.7 mm

W tarczy wykonanej z polimetakrylanu metylu, o wymiarach: szerokość 152.4 mm, wysokość 304.8 mm, grubość t = 3.175 mm, przygotowano szczeliny o długościach l = 7.62, 12.7, 17.78, 25.4 mm, o różnym kacie nachylenia γ względem osi próbki. Próbkę poddano rozciąganiu, rejestrując poziom obciążenia, przy którym następuje zniszczenie. Przykładowe wartości współczynników intensywności naprężeń i współczynników stojących przy członach wyższych rzędów rozwiązania asymptotycznego obliczone za pomocą metody elementów analitycznych przy wykorzystaniu różnej liczby członów (m) rozwiązania asymptotycznego (1) , unormowano zgodnie z zależnościami (15) i przedstawiono na rys. 3.

( )

( )

* *

2 2

I_n = I_n σ π ⁻ⁿ, II_n = II_n σ π ⁻ⁿ

K K l K K l (15)

Wyznaczona w doświadczeniu krytyczna wartość współczynnika intensywności naprężeń wyniosła KIC = 1.370 MPa⋅m^½. W obliczeniach przyjęto wartość naprężeń krytycznych σC = 102.8 MPa, co na podstawie zależności (14) dało wartość strefy pękania d0 = 0.113 mm.

Kryterium pękania (12) pozwoliło wyznaczyć kierunki, w których następować będzie propagacja szczeliny. Wyniki obliczeń przeprowadzonych przy uwzględnieniu różnej liczby członów rozwiązania asymptotycznego zestawiono z wynikami doświadczalnymi na rys. 4.

W przypadku, gdy γ → 0°, wartości współczynników intensywności naprężeń powinny dążyć do zera (KI → 0 i KII → 0) oraz stosunek KI/KII → 0, co sugerowałoby przypadek czystego ścinania wzdłużnego (II sposób obciążania) i wartości kąta propagacji szczeliny ϑ0

zawierające się w przedziale -80°÷ -70°. Liczne badania doświadczalne wykazują, że kąt propagacji w takim przypadku zbliżony jest do -90° [10]. Użycie do obliczeń numerycznych

opisu pola naprężeń wykorzystującego nie tylko człony osobliwe rozwiązania asymptotycznego, ale także wyższych rzędów, daje rezultaty zgodne z eksperymentalnymi (linia m = 2, 3 - rys. 4).

0 15 30 45 60 75 90

0 -10 -20 -30 -40 -50 -60 -70 -80 -90

0 15 30 45 60 75 90

0 -10 -20 -30 -40 -50 -60 -70 -80

l = 7.6 mm -90 wyniki obliczeń dane doświadczalne

l = 12.7 mm

wyniki obliczeń dane doświadczalne

a) b)

m = 2 m = 3 m = 1

m = 2 m = 3 m = 1

Rys.4. Kąt propagacji szczeliny ϑ0 w funkcji kąta nachylenia γ szczeliny środkowej o długości 7.6 mm (a) i 12.7 mm (b), dane doświadczalne – [10]

Kryterium (12) pozwoliło wyznaczyć poziom obciążeń krytycznych, przy których następuje pękanie. Wyniki tych obliczeń, razem z rezultatami badań doświadczalnych, przedstawiono na rys. 5. Krytyczne wartości obciążeń osiągają minimum dla kąta pochylenia szczeliny γ = 70°÷60°, co zgodne jest z wynikami badań doświadczalnych. Jest to wyraźnie widoczne jedynie w przypadku obliczeń uwzględniających efekt członów wyższych rzędów rozwiązania asymptotycznego.

0 15 30

0 2 4 6 8 10 12 14 16 18 20

30 45 60 75 90

0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7

0 15 30

0 2 4 6 8 10 12 14 16 18 20

30 45 60 75 90

0.9 1.0 1.1 1.2 1.3 1.4 1.5

m = 2,3 m = 2

m = 3 m = 1

m = 1 m = 1

1.6

a) 1.7 b)

m = 2,3

m = 2 m = 3 m = 1

l = 7.6 mm

wyniki obliczeń dane doświadczalne

l = 12.7 mm

wyniki obliczeń dane doświadczalne

Rys.5. Wartości krytycznych obciążeń w zależności od kąta γ pochylenia szczeliny środkowej o długości l = 7.6 mm (a) i l = 12.7 mm (b), dane doświadczalne – [10]

4. PODSUMOWANIE

Na podstawie wyników przeprowadzonych obliczeń można stwierdzić, że zastosowanie w obliczeniach opisu pola naprężeń i przemieszczeń wykorzystującego jedynie człony osobliwe rozwiązania asymptotycznego powoduje obarczenie otrzymywanych wyników znacznym błędem numerycznym, w wielu przypadkach dyskwalifikującym obliczenia. Wpływ członów wyższych rzędów w asymptotycznym rozwinięciu, opisującym pola naprężeń w pobliżu wierzchołka szczeliny, jest dominujący dla szczelin pochylonych pod małym kątem do kierunku działania obciążenia (γ = 0°÷10°). Uwzględnienie tylko członów osobliwych w obliczeniach (m = 1) powoduje znaczne zawyżanie poziomu obciążeń krytycznych – zwłaszcza dla małych kątów γ.

LITERATURA

1. Seweryn A.: Metody numeryczne w mechanice pękania. Biblioteka Mechaniki Stosowanej, Seria A. Monografie. Warszawa: IPPT PAN, 2003.

2. He W.J., Lin Y., Ding H.J.: A tree-dimensional formula for determining stress intensity factors in finite element analysis of cracked bodies. “Eng. Fract. Mech.”, 1997, 56, s.409- 415.

3. Yang Z.J., Chen J.F., Holt G.D.: Efficient evaluation of stress intensity factors using virtual crack extension technique. “Comput. Struct.”, 2001, 79, s.2705-2715.

4. Sinclair G.B., Okajima M., Griffin J. H.: Path independent integrals for computing stress intensity factors at sharp notches in elastic plates. “Int. J. Numer. Meth. Eng.”, 1984, 20, s.999-1008.

5. Seweryn A., Adamowicz A.: On analytical constraints and elements methods in modeling stresses near the tips of cracks and V-notches. “Material Science”, 2005, 41, 4.

6. Lin K.Y, Tong P.: Singular finite elements for the fracture analysis of V-notched plate.

“Int. J. Num. Meth. Eng.”, 1980, 15, s.1343-1354.

7. Seweryn A.: Modelling of singular stress fields using finite element method. “Int. J. Solids Struct”., 2002, 39, s.4787-4804.

8. Williams M.L.: On the stress distribution at the base of stationary crack. Trans. ASME,” J.

Appl. Mech.”, 1957, 24, s.109-114.

9. Seweryn A., Mróz Z.: A non-local stress failure condition for structural elements under multiaxial loading. “Eng. Fract. Mech.”, 1995, 51, s.955-973.

10. Williams J.G., Ewing P.D.: Fracture under complex stress - the angled crack problem.

„Int. J. Fract.”, 1971, 8, s.441-446.

Praca naukowa finansowana ze środków budżetowych na naukę w latach 2005 - 2007 jako projekt badawczy nr 4 T07A 030 28

MODELLING OF SINGULAR STRESS FIELDS USING FINITE ANALYTICAL ELEMENT METHODS

IN BRITTLE FRACTURE PROBLEMS

Summary. The paper deals with the problems of applications of analytical elements method in modelling of stress fields near crack tips in elastic bodies. The method of analytical elements is applied to find the stress intensity factors and the coefficients of the higher terms of the asymptotic solution in the case of the sheet containing an angled crack. The derived calculations were used to develop critical condition of crack propagation.