MODELOWANIE PROPAGACJI FAL W PŁYTACH KOMPOZYTOWYCH

MODELOWANIE INŻYNIERSKIE ISNN 1896-771X 32, s. 323-330, Gliwice 2006

MODELOWANIE PROPAGACJI FAL W PŁYTACH KOMPOZYTOWYCH

PAWEŁ KUDELA

Instytut Maszyn Przepływowych PAN w Gdańsku

WIESŁAW OSTACHOWICZ

Instytut Maszyn Przepływowych PAN w Gdańsku Wydział Nawigacyjny, Akademia Morska w Gdyni

Streszczenie. W pracy przedstawiano wyniki symulacji propagacji poprzecznej fali sprężystej w płycie kompozytowej. Zagadnienie to rozwiązywane jest za pomocą zaproponowanych w pracy modeli zbudowanych na podstawie metody spektralnych elementów skończonych. Badano wpływ kąta ułożenia włókien wzmacniających oraz ich procentowej zawartości na propagację fali. Zwrócono uwagę na możliwości wykorzystania zmian w propagującej się fali do wykrywania uszkodzeń w płytach.

1. WSTĘP

W ostatnich dwóch dekadach obserwuje się stały wzrost zainteresowania zastosowaniem materiałów kompozytowych w różnych elementach konstrukcyjnych. Spowodowane to jest faktem, iż materiały kompozytowe cechują znakomite własności mechaniczne. Podobnie jak w przypadku materiałów izotropowych, elementy konstrukcyjne wykonane z materiałów kompozytowych narażone są na różnego rodzaju uszkodzenia (delaminacja, pęknięcie włókien, etc.). Ważnym zagadnieniem jest relatywnie szybkie i trafne wykrycie ewentualnych defektów.

W tym celu wykorzystuje się zmiany różnorodnych wielkości fizycznych określających stan konstrukcji. Najnowsze metody wykrywania uszkodzeń wykorzystują zmiany w propagującej się fali sprężystej.

Aby poprawnie zamodelować zachowanie propagującej się fali sprężystej, istotne jest zbudowanie dokładnego modelu numerycznego. Z uwagi na stosunkowe wysokie częstotliwości propagujących się fal oraz ich na duże prędkości propagacji istotną trudność stanowi dyskretyzacja modelu numerycznego. Oczywistym jest, że siatka klasycznych elementów skończonych musi być w takim przypadku bardzo gęsta, a czas obliczeń bardzo długi. Alternatywą jest zastosowanie metody spektralnych elementów skończonych [1–3], która, dzięki odpowiedniemu doborowi funkcji bazowych i punktów całkowania numerycznego, pozwala rozprzęgnąć równania i znacznie skrócić czas symulacji.

Obliczenia numeryczne przeprowadzono dla różnych parametrów materiału kompozytowego. Jako materiał przyjęto żywicę epoksydową wzmacnianą włóknami szklanymi lub grafitowymi. W pierwszej kolejności analizowano płytę jednowarstwową. Z obliczeń

numerycznych wynika, że prędkość oraz kierunek propagacji fali są funkcjami orientacji włókien oraz ich względnej zawartości. W przypadku wielowarstwowych płyt kompozytowych zachowanie propagującej fali jest bardziej skomplikowane niż dla jednowarstwowej płyty kompozytowej. W płycie wielowarstwowej kształt propagującej fali stanowi wynik superpozycji, rozumianej w sensie homogenizacji fal wytworzonych przez poszczególne warstwy kompozytu.

2. PŁYTOWY SPEKTRALNY ELEMENT SKOŃCZONY 2.1. Zdefiniowanie węzłów elementu

Węzły płytowego spektralnego elementu skończonego zdefiniowane są w lokalnym układzie współrzędnych elementu ξη jako pierwiastki następującego wyrażenia:





′ =

−

′ =

−

0 ) ( ) 1 (

0 ) ( ) 1 (

2 2

η η

ξ ξ

N N

P

P (1)

gdzie ξ, η∈[–1, 1], P^’N oznacza pierwszą pochodną wielomianu Legendre’a rzędu N. W ten sposób węzły elementu mogą zostać określone w lokalnym układzie współrzędnych elementu.

W bieżącym sformułowaniu wybrano wielomian Legendre’a 5 rzędu. Stąd można otrzymać 36 węzłów w układzie współrzędnych ξη jako:

} 1 , ,

, ,

, 1 { ,

6 ..., , 2 , 1 , ), , (

7 3

2 3 1 7 3

2 3 1 7 3

2 3 1 7 3

2 3

1 + − − − +

−

−

∈

=

n m

n

m m n

η ξ

η ξ

(2)

Na rysunku 1 można zauważyć, że wynikiem takiej definicji węzłów elementu jest ich nierównomierny rozkład wewnątrz elementu.

η

ξ

η=-1 η=1

ξ=-1 ξ=1

Rys. 1. Rozkład węzłów w lokalnym układzie współrzędnych elementu

2.2. Funkcje kształtu w elemencie

Na podstawie tak zdefiniowanych węzłów można zbudować zbiór funkcji kształtu aproksymujących poprzeczne i kątowe przemieszczenia wewnątrz elementu. W tym celu zastosowano dyskretnie ortogonalną aproksymację Lagrange’a opartą na węzłach elementu.

2.3. Macierze sztywności i mas

Kompozytowy element płytowy opracowano na podstawie teorii Mindlina [4, 5].

Umożliwia to dokładne modelowanie zjawiska propagacji fali (w przeciwieństwie do klasycznej teorii płyt laminowanych). Obliczenie charakterystycznych macierzy sztywności [K]

i bezwładności [M] skończonego płytowego elementu spektralnego przebiega w sposób identyczny z klasycznym sformułowaniem metody elementów skończonych. Proces ten może być opisany w postaci związków:

∫ ∫

∑∑

= ⁶

1 6

1 1

1 1

1[ ] [ ]det[ ] [ ] [ ]det[ ]

] [

m n

mn mn

T mn n m

T N J d d w w N N J

N

M ρ η ξ ρ (3)

∫ ∫

∑∑

+

− =

= ⁶

1 6

1 1

1 1

1[ ] [ ][ ]det[ ] [ ] [ ][ ]det[ ]

] [

m n

mn mn

T mn n m

T D B J d d w w B D B J

B

K η ξ (4)

gdzie [N] jest macierzą funkcji kształtu, [J] jest macierzą Jakobianów, [D] oznacza macierz stałych sprężystości materiału kompozytowego, natomiast w i _m w są wagami Gaussa–_n Lobatto, które obliczono w węzłach elementu m i n za pomocą wzoru [6]:

6 , , 1 , ), ( 15

1

, 5

, = m n= K

w P

n m n

m ξ . (5)

Z uwagi na ortogonalność funkcji kształtu i zastosowanie reguły całkowania Gaussa–Lobatto macierz bezwładności elementu [M] jest diagonalna. Własność ta wynika bezpośrednio z definicji. W przypadku techniki skupiania mas można również otrzymać diagonalną macierz bezwładności. Wówczas jednak w rozwiązaniu równań ruchu pojawiają znaczące błędy, [7].

2.4. Rozwiązanie równań ruchu

W przypadku braku tłumienia równania ruchu można zapisać w postaci macierzowej:

t t

t K q F

q

M]{ } [ ]{ } { }

[ && + = (6)

gdzie symbol t oznacza czas {F} jest wektorem sił wzbudzających, symbol {q} oraz {q&& } oznaczają wektory węzłowych przemieszczeń oraz przyspieszeń. Ponieważ globalna macierz bezwładności [M] w równaniu ruchu przyjmuje postać diagonalną, równanie (6) można dalej uprościć i rozwiązać, stosując bardzo szybką i efektywną metodę całkowania opartą na schemacie różnic centralnych [8]. Ponadto macierz sztywności jest stosunkowo rzadka (zawiera około 35% elementów niezerowych), a więc wymagane są mniejsze zasoby pamięci niż w przypadku klasycznej metody elementów skończonych.

3. OBLICZENIA NUMERYCZNE

3.1. Geometria płyty i dane materiałowe

Do obliczeń przyjęto płytę kompozytową o wymiarach: długość (100 cm), szerokość (100 cm). Przyjęto założenie, że płyta została wykonana z żywicy epoksydowej wzmacnianej włóknami szklanymi lub grafitowymi. Własności materiałowe zostały zestawione w tabeli 1.

Przeprowadzono badania płyt kompozytowych jednowarstwowych o całkowitej grubości równej 1 cm. W dalszej kolejności badano płyty dziesięciowarstwowe o tej samej grubości.

Tabela 1. Własności materiałowe Żywica Włókna szklane Włókna grafitowe

Moduł Younga [GPa] 3.43 66.5 275.6

Współczynnik Poissona 0.35 0.23 0.20

Gęstość [kg/m³] 1200 2250 1900

3.2. Wpływ objętościowej zawartości włókien wzmacniających na propagację fali

Prędkość grupowa cg fali poprzecznej jest funkcją względnej objętościowej zawartości włókien wzmacniających, kierunku propagacji fali oraz częstotliwości sygnału wzbudzenia [9].

Prędkość ta może zostać obliczona analitycznie, w wyniku czego otrzymujemy profile prędkości pokazane na rysunkach 2a i 2b. Założono, że kąt ułożenia włókien jest stały i wynosi 0˚. W tym też kierunku wartości prędkości grupowych są największe.

Rys. 2. Profile prędkości grupowych w jednowarstwowej płycie kompozytowej wzmacnianej włóknami a) szklanymi, b) grafitowymi, dla zawartości włókien 0, 20, 40, 60, 80, 100%

Wraz ze wzrostem zawartości włókien wzmacniających prędkość fali w kompozycie wzrasta, przy jednoczesnej zmianie kształtu profilu prędkości. Okręgi na rysunkach 2a i 2b odnoszą się do przypadków izotropowych, kiedy prędkość fali w każdym kierunku jest jednakowa.

a) b)

Zależnie od zastosowanego materiału włókien kształty profili prędkości różnią się. Z uwagi na znacznie wyższy moduł Younga włókien grafitowych od modułu Younga włókien szklanych prędkość fali na rysunku 2b jest wyższa, a kształty profili prędkości są mniej gładkie.

3.3. Wpływ kąta ułożenia włókien na propagację fali

Prędkość grupowa fali zależy również od kąta ułożenia włókien wzmacniających.

Przeprowadzono teoretyczne i numeryczne obliczenia dla płyty kompozytowej o stałej zawartości włókien wzmacniających. Wyniki wskazują, że kształt frontu propagującej fali jest zachowany, podczas gdy eliptyczne wydłużenie obrócone jest zgodnie z zadaną orientacją włókien wzmacniających.

3.4. Propagacja fali w płycie wielowarstwowej

W tym przypadku zachowanie propagującej fali jest bardziej skomplikowane niż obserwowane w przypadku pojedynczej warstwy kompozytu. Kształty profili prędkości grupowych zależą od objętościowej zawartości włókien wzmacniających, jak i od ich orientacji w poszczególnych warstwach. Można założyć, że propagująca fala stanowi wynik superpozycji (rozumianej w sensie homogenizacji) fal z poszczególnych warstw kompozytu.

Rysunek 3 przedstawia fragmenty symulacji numerycznych dla dwóch przykładowych płyt kompozytowych. W obu przypadkach założono, że płyta składa się z 10 warstw, zawartość włókien wzmacniających w każdej z nich jest jednakowa i wynosi 25%. W lewej kolumnie symulacja dotyczy płyty wykonanej z żywicy epoksydowej wzmacnianej włóknami grafitowymi, a w kolumnie po prawej stronie płyty wykonanej z żywicy epoksydowej wzmacnianej włóknami szklanymi. Różne układy orientacji włókien w poszczególnych warstwach wpływają na kształt propagującej się fali. Znajomość rozkładu prędkości w zależności od kierunku propagacji jest niezwykle istotna z uwagi na zastosowanie do detekcji uszkodzeń.

Na szczególną uwagę zasługuje dobór odpowiedniej częstotliwości sygnału wzbudzającego.

Wiąże się to z efektem dyspersji, czyli zależności prędkości fali od częstotliwości. Rysunek 4 przedstawia przykładowe krzywe dyspersji dla płyty kompozytowej o układzie warstw θ=[+60^o/–60^o]5 i grubości każdej z warstw 1 mm. Dyspersja powoduje deformacje sygnału wymuszającego, co z kolei pociąga za sobą błędy w szacowaniu prędkości fali. Tradycyjnie do detekcji uszkodzeń stosowane są sygnały wymuszające w postaci sygnału typu sinusoidalnego, modulowanego za pomocą okna (hanning, trójkąt), które tworzą tzw. paczki [10]. Efekt dyspersji można zminimalizować poprzez stosowanie tego typu sygnałów wejściowych o wąskim zakresie częstotliwości i przez skupienie energii wejściowej w punkcie na krzywej dyspersji, gdzie dyspersja jest niska, [11]. Na rysunku 4 miejsce takie stanowi zakres częstotliwości od około 15 kHz do około 50 kHz, gdzie krzywa jest najbardziej płaska.

Jednocześnie w zakresie tym propaguje tylko jeden mod, dzięki temu łatwiej analizować sygnał wyjściowy.

w. grafitowe, θ=[+60^o/–60^o]5, vol=25% w. szklane, θ=[0^o/90^o]5, vol=25%

t=0.24 ms t=0.24 ms

t=0.36 ms t=0.36 ms

t=0.48 ms t=0.48 ms

Rys. 3. Symulacje numeryczne propagującej fali poprzecznej w wielowarstwowych płytach kompozytowych

0 25 50 75 100 125 150 175 200 0

0.2 0.4 0.6 0.8 1 1.2 1.4

Częstotliwość [kHz]

Prędkość grupowa [km/s]

Rys. 4. Przykładowe krzywe dyspersji dla wielowarstwowej płyty kompozytowej

4. WNIOSKI

Metoda spektralnych elementów skończonych stanowi efektywne i dokładne narzędzie do modelowania zjawiska propagacji fal w płytach kompozytowych. Przeprowadzone symulacje numeryczne pozwalają obserwować zachowanie propagującej fali i dostarczają cennych informacji, szczególnie w kontekście detekcji uszkodzeń.

W płytach kompozytowych fala propaguje w każdym kierunku z inną prędkością. Dlatego też w rzeczywistych systemach detekcji uszkodzeń bazujących na czujnikach piezoelektrycznych preferowane będą układy koncentryczne (np. czujniki rozmieszczone w układzie zegarowym). Wtedy istnieje możliwość zmierzenia odpowiednich czasów propagacji fali i uwzględnienia kształtu profilu prędkości (w pewnym przybliżeniu).

Pokazano, że na przebieg propagującej fali wpływają również czynniki takie jak:

objętościowa zawartość włókien wzmacniających, rodzaj zastosowanych włókien, układ warstw. Szczególną uwagę należy zwrócić na dobór odpowiedniej częstotliwości sygnału wymuszającego, w celu zminimalizowania efektu dyspersji.

PODZIĘKOWANIA

Autorzy niniejszej pracy pragną wyrazić podziękowania za wsparcie finansowe prowadzonych przez nich badań w ramach projektu ARTIMA (Aircraft Reliability Through Intelligent Materials Application – numer referencyjny 502725), prowadzonego w ramach 6 Programu Ramowego Unii Europejskiej, jak również pragną podziękować Ministerstwu Edukacji i Nauki za dodatkowe wsparcie finansowe tego projektu w postaci projektu SPUB– ARTIMA).

LITERATURA

1. Patera A.T.: A spectral element method for fluid dynamics: Laminar flow in a channel expansion. Journal of Computational Physics, 54, 1984, p. 468–488.

2. Boyd J.P.: Cheybyshev and Fourier: Spectral Methods. Springer 1989.

3. Pozrikidis C.: Introduction to Finite and Spectral Element Methods using MATLAB®.

Chapman & Hall/CRC, 2005.

4. Ochoa O.O., Reddy J.N.: Finite element analysis of composite laminates. Kluwer Academic Publishers, 1992.

5. Vinson J.R., Sierakowski R.L.: Behavior of structures composed of composite materials.

Martinus–Nijhoff Inc, 1989.

6. http://mathworld.wolfram.com/LobattoQuadrature.html

7. Dauksher W., Emery A.F.: Accuracy in modeling the acoustic wave equation with Chebyshev spectral finite elements. Finite Elements in Analysis and Design, 26, 1997, p. 115–128.

8. Kleiber M.: Incremental finite element modelling in non–linear solid mechanics. J. Wiley &

Sons, New York, 1989.

9. Liu G.R., Xi Z.C.: Elastic Waves in Anisotropic Laminates, CRC Press, 2002.

10. Doyle J.F.: Wave Propagation in Structures. Springer–Verlag, 1997.

11. Wilcox P., Lowe M., Cawley P.: The effect of dispersion on long–range inspection using ultrasonic guided waves. NDT&E Int., 34, 2001, p. 1–9.

WAVE PROPAGATION MODELLING IN COMPOSITE PLATES

Summary. This paper presents results of numerical simulation of the propagation of a transverse elastic wave in a composite plate. The problem is solved by the use of proposed algorithm formulated on the base of the Spectral Finite Element Method. The influence of orientations of reinforcing fibres and their volume fraction on wave propagation has been investigated. It has been paying attention to application of changes in propagating waves on damage detection in structures.