SZYBKA WIELOBIEGUNOWA METODA ELEMENTÓW BRZEGOWYCH W HOMOGENIZACJI NUMERYCZNEJ MATERIAŁÓW POROWATYCH

(1)

SZYBKA WIELOBIEGUNOWA

METODA ELEMENTÓW BRZEGOWYCH W HOMOGENIZACJI NUMERYCZNEJ MATERIAŁÓW POROWATYCH

Jacek Ptaszny

^1a

, Marcin Hatłas

^2b

1Instytut Mechaniki i Inżynierii Obliczeniowej, Politechnika Śląska

2Student II stopnia, Wydział Mechaniczny Technologiczny, Politechnika Śląska

ajacek.ptaszny@polsl.pl, ^bmarcin.hatlas91@gmail.com

Streszczenie

W pracy przedstawiono nową wersję szybkiej wielobiegunowej metody elementów brzegowych służącą do analizy zagadnień przestrzennych liniowej teorii sprężystości, w której zastosowano elementy brzegowe o kwadratowych funkcjach kształtu oraz metodę adaptacyjnego całkowania. Metoda została zastosowana w homogenizacji materiałów porowatych zawierających pustki sferyczne. W tym celu analizowano reprezentatywne elementy obję- tości zawierające dużą liczbę pustek. W efekcie obliczono zastępcze stałe sprężystości materiału. Wyniki porówna- no z dostępnymi modelami analitycznymi otrzymanymi metodą samospójną oraz Mori-Tanaki. Porównanie wskazuje na poprawność opracowanej metody oraz wykonanych obliczeń numerycznych.

Słowa kluczowe: szybka metoda wielobiegunowa, metoda elementów brzegowych, materiały porowate, homogenizacja, sprężystość liniowa

FAST MULTIPOLE BOUNDARY ELEMENT METHOD IN THE NUMERICAL HOMOGENIZATION OF POROUS MATERIALS

Summary

In this work, a new version of the fast multipole boundary element method for three-dimensional linear elasticity problems, with boundary elements with quadratic shape functions and adaptive integration, is presented. The method was applied in the homogenization of porous materials with spherical cavities. Representative volume elements containing a high number of cavities were analysed. As a results, overall elastic constants of the porous materials were calculated. The results were compared to available analytical models obtained by the self- consistent and Mori-Tanaka methods. The comparison confirmed vailidity of the method and computations.

Keywords: fast multipole method, boundary element method, porous materials, homogenization, linear elasticity

1. WSTĘP

Zastosowanie materiałów niejednorodnych wymaga znajomości ich własności zastępczych oraz wytrzymało- ściowych. Własności te mogą być wyznaczone za pomocą badań doświadczalnych, metod analitycznych oraz numerycznych. Zwiększenie możliwości obliczeniowych komputerów powoduje ciągły wzrost znaczenia metod

numerycznych. Metody te zapewniają niski koszt badań w porównaniu z metodami doświadczalnymi oraz brak ograniczeń co do geometrii analizowanych układów.

Wśród metod numerycznych najbardziej popularną jest metoda elementów skończonych (MES). Wynika to w dużej mierze z dostępności komercyjnych programów

(2)

komputerowych. Metoda ta wymaga jednak dyskretyzacji całego obszaru, co wpływa niekorzystnie na rozmiar zbioru danych, które należy przygotować, układu rów- nań, który należy rozwiązać oraz rozmiar zbioru warto- ści wynikowych. Korzystną alternatywą może być zastosowanie metody elementów brzegowych (MEB), która w wielu przypadkach o znaczeniu praktycznym wymaga dyskretyzacji jedynie brzegu analizowanego obszaru [1, 2, 11, 12]. Zaletą MEB jest również większa dokładność rozwiązania zagadnień z dużym spiętrzeniem naprężeń, w porównaniu z MES. Wadą konwencjonalnej MEB jest natomiast struktura macierzy, które są pełne i niesymetryczne. Etap budowy układu równań jest operacją rzędu O(N²), gdzie N jest liczbą stopni swobo- dy. Złożoność powoduje, że metoda jest nieefektywna w rozwiązywaniu dużych układów, tzn. przy N rzędu 10⁴ oraz większym, ze względu na wymaganą pamięć kom- putera oraz czas obliczeń. Wymienione niedogodności mogą być wyeliminowane przez zastosowanie nowych wersji metody. Jedną z nich jest szybka wielobiegunowa MEB (SWMEB). W metodzie tej układ równań jest rozwiązywany iteracyjnie, zaś złożoność operacji związa- nych z obliczaniem całek brzegowych jest rzędu O(N).

Metoda wykorzystuje rozwinięcie całek brzegowych w szereg wielobiegunowy oraz jego transformacje prowa- dzące do zmniejszenia liczby operacji całkowania [8, 12].

Metoda elementów brzegowych była stosowana w analizie płaskich i przestrzennych układów zawierają- cych pustki, pęknięcia, wtrącenia i włókna [3, 4, 7, 9, 10, 12, 13, 15, 16, 17, 18, 19, 20, 21]. W wybranych pracach stosowano elementy brzegowe o stałych funkcjach kształtu z analitycznym całkowaniem [12, 15]. Analiza z wykorzystaniem takich elementów może nie być jednak efektywna w przypadku układów poddanych zgina- niu [14, 22].

W ramach niniejszej pracy opracowano program komputerowy SWMEB do analizy zagadnień przestrzennych liniowej teorii sprężystości, wykorzystujący elementy ośmiowęzłowe o kwadratowych fukcjach kształtu oraz całkowanie adaptacyjne z podziałem na podelementy [1, 6]. Taka metoda całkowania nie była jeszcze stosowana w połączeniu z szybką metodą wielobiegunową.

Nową wersję metody zastosowano w homogenizacji numerycznej materiałów porowatych modelowanych jako reprezentatywne elementy objętości z dużą liczbą pustek sferycznych rozmieszczonych losowo [23]. Liczba stopni swobody analizowanych modeli przekraczała 120 000. Wyniki homogenizacji porównano z wynikami analitycznymi uzyskanymi metodą samospójną oraz Mori-Tanaki [5].

Treść niniejszego artykułu została podzielona na pięć rozdziałów. Rozdział drugi zawiera ogólny opis SWMEB. W rozdziale trzecim podano podstawowe zależności stosowane podczas homogenizacji numerycznej. Rozdział czwarty zawiera przykład homogenizacji

z wykorzystaniem SWMEB. Rozdział piąty zawiera wnioski.

2. SZYBKA WIELOBIEGUNOWA METODA ELEMENTÓW

BRZEGOWYCH

Niniejszy rozdział zawiera jedynie podstawowe in- formacje dotyczące SWMEB ze względu na ograniczenie objętości artykułu. Więcej szczegółów dotyczących metody można znaleźć w literaturze, np. [3, 13].

Rozpatrywane jest jednorodne ciało Ω wykonane z izotropowego materiału liniowosprężystego, o brzegu Γ. Ciało obciążone jest siłami powierzchniowymi tj

(j = 1, 2, 3) na wybranym fragmencie brzegu. Na pozo- stałej części brzegu znane są przemieszczenia. Związek pomiędzy siłami powierzchniowymi oraz przemieszcze- niami uj opisany jest brzegowym równaniem całkowym:

, d

, (1)

gdzie x^’ jest punktem kolokacji, x jest punktem całko- wania, cij jest współczynnikiem zależnym od położenia x’, natomiast Tij, Uij są rozwiązaniami podstawowymi zagadnienia sprężystości [1, 2, 12]. Całki występujące w równaniu (1) nazywane są odpowiednio potencjałem warstwy podwójnej oraz pojedynczej o gęstościach uj

oraz tj.

Nieznane wielkości brzegowe można wyznaczyć, sto- sując aproksymację geometrii układu oraz wielkości brzegowych poprzez podział brzegu na elementy brzegowe, zdefiniowane zbiorem węzłów oraz funkcji kształtu.

W niniejszej pracy stosowano ośmiowęzłowy element o kwadratowych funkcjach kształtu (Serendipa) pokaza- ny na rys. 1 [1].

a) b)

Rys. 1. Ośmiowęzłowy element brzegowy o kwadratowych funkcjach kształtu: a) w globalnym układzie współrzędnych

xi (i = 1, 2, 3), b) w lokalnym układzie współrzędnych ξ, η Gdy punkt kolokacji pokrywa się z punktem całko- wania, całki brzegowe występujące w równaniu (1) są osobliwe. W opisywanym zagadnieniu występują osobli- wości typu 1/r (rozwiązanie podstawowe Uij) oraz 1/r² (rozwiązanie podstawowe Tij), gdzie r = |x-x’|.

W przypadku pierwszego z wymienionych typów osobli- wości stosowano podział elementów na elementy trój- kątne oraz transformację układu współrzędnych prowa- dzącą do regularyzacji całek. W drugim przypadku

(3)

stosowano metodę ruchu ciała sztywnego. Całki nieoso- bliwe obliczano stosując kubatury Gaussa.

Gdy dwa fragmenty brzegu są położone blisko siebie, stosowanie standardowych procedur całkowania nume- rycznego może prowadzić do znaczących błędów całko- wania. W celu zminimalizowania błędu można stosować metody regularyzacji całek lub całkowania adaptacyjnego, polegające na podziale elementów na mniejsze elementy. W niniejszej pracy stosowano procedurę adaptacyjnego całkowania opisaną w książce [1], która nie była jeszcze stosowana w połączeniu z szybką metodą wielo- biegunową. Procedura ta bazuje na obliczaniu odległości punktu kolokacji od elementu brzegowego zawierającego punkty całkowania i określaniu na tej podstawie liczby punktów Gaussa. Jeśli liczba punktów dla jednego kierunku przekracza 4, element jest dzielony w danym kierunku na podelementy. Dla rozpatrywanych typów osobliwości procedura zapewnia całkowanie z błędem względnym nieprzekraczającym 10^-3.

Szybka wielobiegunowa MEB wykorzystuje następu- jące rozwinięcie jądra potencjału warstwy pojedynczej:

, ∑ ∑ _, _! " , ! , (2)

gdzie Rn,m oraz Sn,m są funkcjami kulistymi zależnymi od współrzędnych argumentu funkcji w biegunowym ukła- dzie współrzędnych oraz od stowarzyszonych funkcji Legendre’a [12]. Symbol xM oznacza punkt położony w pobliżu punktów całkowania x. Równanie (2) może być bezpośrednio użyte do rozwinięcia w szereg poten- cjału warstwy pojedynczej, zależnego od 1/r. Rozwinię- cie potencjału warstwy podwójnej, zależnego od 1/r², wymaga obliczenia pochodnych funkcji kulistych wystę- pujących w równaniu (2).

Zastosowanie szeregów umożliwia efektywne oblicze- nie potencjałów wielu punktów całkowania dla wielu punktów kolokacji. Wyrazy szeregu zbudowanego do obliczenia potencjałów mają postać sumy iloczynów tzw.

funkcji wielobiegunowych zależnych od Sn,m(x’-xM) oraz momentów zależnych od Rn,m(x-xM) i od gęstości poten- cjałów (wielkości brzegowych w punktach całkowania).

W ten sposób składnik potencjału pochodzący od wielu punktów całkowania może być zredukowany w pojedyn- czym punkcie xM. Kolejną operacją jest przesunięcie punktu rozwinięcia xM do nowego położenia (tzw. trans- lacja M2M, ang. multipole-to-multipole) oraz sumowanie potencjału pochodzącego od wielu grup elementów brzegowych. Następnie buduje się tzw. szereg lokalny poprzez transformację M2L (ang. multipole-to-local), w punkcie xL położonym w pobliżu punktów kolokacji.

Kolejna transformacja, tzw. L2L (ang. local-to-local), pozwala na przesunięcie punktu xL w nowe położenie i dystrybucję potencjału do wielu obszarów zawierają- cych punkty kolokacji. Wreszcie, stosując szereg wielobiegunowy, można obliczyć potencjał dla wielu punktów kolokacji. Taki sposób obliczania potencjałów pozwala

na znaczne zredukowanie liczby operacji całkowania w stosunku do konwencjonalnej MEB, gdzie dla każdego punktu kolokacji oblicza się całki po powierzchni wszystkich elementów brzegowych. Schemat redukcji i dystrybucji potencjałów w SWMEB został przedsta- wiony na rys. 2. Dla uproszczenia przedstawiono układ płaski, jednak ogólna zasada nie zależy od wymiaru zagadnienia i jest stosowana również w przypadku zagadnień przestrzennych.

Metoda wymaga rekurencyjnego grupowania elemen- tów brzegowych w obszarach, a ich hierarchia jest zapisywana w postaci struktury drzewa. W przypadku zagadnień przestrzennych elementy brzegowe grupowane są w obszarach będących sześcianami.

Rys. 2. Schemat obliczania potencjałów stosowany w SWMEB Równanie (1) jest budowane dla wszystkich węzłów brzegowych jako punktów kolokacji. W wyniku dyskretyzacji geometrii oraz wielkości brzegowych, przy której uwzględnia się przyporządkowanie lokalnych numerów węzłów elementów numerom globalnym (incydencji), uzyskuje się układ równań w postaci macierzowej:

#$%^&'()* ($)*^+,' #-%^&'(.* (-.*^+,'. (3) Macierze [H]^bl i [G]^bl zawierają całki zależne od rozwią- zań podstawowych obliczane w sposób bezpośredni.

Pochodzą one od tzw. obszaru bliskiego punktu kolokacji, gdzie szeregi nie są zbieżne i nie mogą być stosowa- ne. Wektory {Hu}^odl oraz {Gt}^odl zawierają składniki potencjałów obliczone za pomocą szeregów. Wektory {u} oraz {t} zawierają brzegowe przemieszczenia i siły powierzchniowe. W wyniku uwzględnienia warunków brzegowych, polegającym na zgrupowaniu niewiadomych przemieszczeń i sił powierzchniowych po lewej stronie w wektorze {x}, a znanych po prawej stronie w wekto- rze {y}, układ równań przyjmuje postać:

#0%^&'( * (0 *^+,' #1%^&'(2* (12*^+,'. (4) Macierze [A]^bl i [B]^bl oraz wektory {Ax}^odl oraz {By}^odl zawierają odpowiednie elementy macierzy występujących w równaniu (3). Układ równań (4) jest rozwiązywany za pomocą iteracyjnej metody GMRES z poprawą uwarun- kowania układu równań za pomocą bloków diagonalnych macierzy [A]^bl [21].

(4)

3. HOMOGENIZACJA NUMERYCZNA

Zagadnienie numerycznej homogenizacji polega na wyznaczeniu zastępczych stałych materiału niejednorod- nego [23]. Stałe te występują w równaniu konstytutyw- nym opisującym związek pomiędzy naprężeniami i odkształceniami w skali makro:

〈σ 〉 6^∗₇₈〈ε78〉, ;, <, =, > 1, 2, 3. (5) W równaniu tym kolejne symbole oznaczają: 〈σij〉 - średnie naprężenia, 〈εij〉 - średnie odkształcenia oraz C∗78 - tensor zastępczych stałych sprężystych materiału.

Stosując notację Voigta równanie (5) można zapisać w następującej postaci:

CDD E DDF〈σ 〉

〈σGG〉

〈σHH〉

〈σGH〉

〈σH 〉

〈σ_G〉IDDJ DDK

LM MM MM N66GG^∗∗

6HH∗

6GH∗

6H∗

6∗G

6∗GG

6GGGG∗

6HHGG∗

6GHGG∗

6H GG∗

6∗GGG

6∗HH

6GGHH∗

6HHHH∗

6GHHH∗

6H HH∗

6∗GHH

6∗GH

6GGGH∗

6HHGH∗

6GHGH∗

6H GH∗

6∗GGH

6∗H

6GGH∗

6HHH∗

6GHH∗

6H H∗

6∗GH

6∗ G

6GG G∗

6HH G∗

6GH G∗

6H G∗

6∗G GOPPPPPQ CDD E DDF〈ε 〉

〈εGG〉

〈εHH〉 2〈εGH〉 2〈εH 〉 2〈εG〉IDDJ

DDK . (6)

W ogólnym przypadku można wyznaczyć wszystkie elementy C∗78 macierzy sztywności materiału, która będzie dalej oznaczona symbolem c^*, wykonując sześć niezależnych testów.

Wartości uśrednione występujące w równaniach (5) i (6) można obliczyć, stosując równania:

〈R 〉 _|T

U| _U V d , (7)

〈W 〉 _G|T

U| _UX Y Y Zd . (8) Poszczególne symbole oznaczają: Ω0 - obszar reprezentatywnego elementu objętości, Γ0 – brzeg zewnętrzny reprezentatywnego elementu objętości, ti – składowe sił powierzchniowych, ui – składowe przemieszczeń, ni – składowe jednostkowego wektora normalnego do Γ0.

W celu rozwiązania zagadnienia brzegowego w skali mikro należy sformułować na brzegu Γ0 odpowiednie warunki brzegowe. Najczęściej stosowanymi w homogenizacji numerycznej warunkami brzegowymi są warunki w postaci sił powierzchniowych, przemieszczeniowe oraz periodyczne warunki brzegowe. W niniejszej pracy stosowano przemieszczeniowe warunki brzegowe opisane równaniem:

|_U 〈W 〉 ∙ V . (9)

Wymienionych wcześniej sześć testów może odpo- wiadać następującym zadanym postaciom tensora śred- nich odkształceń:

〈W 〉^{\ ]\} ^β 0 0 0 0 0

0 0 0a , ^0 0 0 0 β 0

0 0 0a , ^0 0 0 0 0 0 0 0 βa,

^0 0 0 0 0 β

0 β 0a , ^0 0 β 0 0 0

β 0 0a , ^0 β 0 β 0 0

0 0 0a. (10)

4. PRZYKŁAD HOMOGENIZACJI

Analizowano reprezentatywne elementy objętości materiału porowatego zawierającego losowo rozmieszczone jednakowe pustki sferyczne. Liczba pustek wynosiła 125. Porowatość f modeli zmieniano w zakresie od 0.05 do 0.25 z krokiem 0.05, dostosowując promień pustek. Taki zakres porowatości pozwolił na stosunkowo łatwe i szybkie wygenerowanie losowej geometrii modeli.

Obszar układów był ograniczony sześcianem o długości boku równej 1 mm. Stałe materiału jednolitego wynosiły: moduł Younga E0 = 2⋅10⁵ MPa oraz liczba Poissona ν0 = 0.3. Liczba elementów brzegowych wynosiła 13 536, liczba węzłów 40 860, a liczba stopni swobody modeli 122 580. Wnętrze zdyskretyzowanych modeli materiału o największej i najmniejszej porowatości pokazano na rys. 3. Tolerancja metody GMRES wynosiła 10^-6. Rząd szeregu wielobiegunowego wynosił 12.

a)

b)

Rys. 3. Zdyskretyzowana geometria badanych układów o porowatości: a) f = 0.05, b) f = 0.25

Z uwagi na geometrię pustek oraz ich rozmieszczenie analizowany materiał może być uznany za izotropowy w skali makro. Macierz sztywności w takim przypadku ma postać:

b^∗ LM MM

MNκ^∗ 4/3μ^∗ κ^∗ 2/3μ^∗ κ^∗ 2/3μ^∗

00 0

κ^∗ 2/3μ^∗ κ^∗ 4/3μ^∗ κ^∗ 2/3μ^∗

00 0

κ^∗ 2/3μ^∗ κ^∗ 2/3μ^∗ κ^∗ 4/3μ^∗

00 0

0 0 μ0^∗ 00

0 0 00 μ^∗

0

0 0 00 μ0^∗OPPPPQ

. (11)

(5)

Symbole κ^* oraz µ^* oznaczają odpowiednio zastępcze moduły Helmholtza oraz Kirchhoffa. Obydwie stałe są zawarte w elementach macierzy tworzących kolumny 1÷3. W niniejszym przykładzie stosowano test odpowia- dający składowym odkształceń 〈W 〉^\ przy β = 1 (równa- nia (9) i (10)). Znając pierwszą kolumnę macierzy c^* można wyznaczyć zastępcze stałe sprężystości materiału za pomocą równań:

g^∗ ^hⁱⁱ^∗^jGh_H ^ki^∗ , l^∗ ^hⁱⁱ^∗ _G^h^ki^∗ . (12) Wyniki homogenizacji porównano z rozwiązaniem analitycznym uzyskanym za pomocą uogólnionej metody samospójnej (ang. generalized self consistent method, GSCM) oraz metodą Mori-Tanaki (M-T) opisanymi w pracy [5]. W metodzie GSCM zakłada się, że ośrodek zastępczy, w którym umieszczone jest wtrącenie, odpowiada założeniu Voigta (wersja VGSCM metody) lub Reussa (wersja RGSCM). W ten sposób uzyskuje się granice, w których powinny mieścić się zastępcze stałe materiału kompozytowego ze sferycznymi wtrąceniami, którego szczególnym przypadkiem jest analizowany materiał porowaty. Przypadek wtrąceń o małej sztywności, którego granicznym przypadkiem są pustki, jest bliższy wersji VGSCM [5]. W przypadku modułu Helmholtza obydwa warianty metody prowadzą do równania, które uzyskuje się również metodą Mori- Tanaki (M-T). Porównanie zastępczych stałych znormalizowanych względem modułów materiału bez pustek, odpowiednio κ0 i µ0, przedstawiono w tabelach 1 i 2 oraz na rysunkach 4 i 5.

Tabela 1. Porównanie wyznaczonych wartości znormalizowanego zastępczego modułu Helmholtza κ^*/κ0

f GSCM/M-T SWMEB

0.05 0.877 0.877

0.10 0.774 0.772

0.15 0.683 0.681

0.20 0.604 0.602

0.25 0.533 0.531

Tabela 2. Porównanie wyznaczonych wartości znormalizowanego zastępczego modułu Kirchhoffa µ^*/µ0

f VGSCM RGSCM M-T SWMEB

0.05 0.909 0.899 0.909 0.910

0.10 0.824 0.793 0.825 0.826

0.15 0.745 0.691 0.748 0.748

0.20 0.672 0.595 0.677 0.673

0.25 0.603 0.509 0.611 0.599

Rys. 4. Porównanie znormalizowanego zastępczego modułu Helmholtza wyznaczonego różnymi metodami

Rys. 5. Porównanie znormalizowanego zastępczego modułu Kirchhoffa wyznaczonego różnymi metodami

Porównanie wskazuje na poprawność otrzymanych wyników. Różnica względna pomiędzy wynikami nume- rycznymi oraz modelem VGSCM nie przekracza 1%.

5. WNIOSKI

W pracy zaprezentowano zastosowanie nowej wersji szybkiej wielobiegunowej MEB w numerycznej homogenizacji materiałów zawierających losowo rozmieszczone pustki sferyczne. Wyznaczano zastępcze stałe sprężysto- ści tych materiałów. Wyniki homogenizacji numerycznej są zgodne z odpowiednimi modelami analitycznymi uzyskanymi za pomocą metody samospójnej oraz Mori- Tanaki. Mała różnica względna wyników (poniżej 1%) wskazuje na poprawność opracowanej metody oraz wykonanych obliczeń. Opracowana metoda może być stosowana w analizie struktur porowatych. Możliwe jest również zastosowanie sformułowań pozwalających na modelowanie materiałów kompozytowych zawierających różnego rodzaju wzmocnienie [3] i w rezultacie opraco- wanie efektywnej metody analizy szerokiej grupy mate- riałów niejednorodnych.

Niniejsza praca została zrealizowana częściowo w ramach projektu 10/040/BK_15/0006.

(6)

Literatura

1. Beer, G., Smith I., Duenser C.: The boundary element ethod with pProgramming. Wien: Springer-Verlag, 2008.

2. Brebbia C.A., Dominguez J.: Boundary elements: an introductory course. Southampton: WITPress-Computational Mechanics Publications, 1992.

3. Czyż T., Dziatkiewicz G., Fedeliński P (red.), Górski R., Ptaszny J.: Advanced computer modelling in micromechanics.

Gliwice: Wyd. Pol. Śl., 2013.

4. Chen X.L., Liu Y.J.: An advanced 3D boundary element method for characterizations of composite materials. “Engi- neering Analysis with Boundary Elements” 2005, 29, p. 513-523.

5. Dai L.H., Huang Z.P., Wang R.: Explicit expressions for bounds for the effective moduli of multi-phased composites by the generalized self-consistent method. “Composites Science and Technology” 1999, 59, p. 1691-1699.

6. Eberwien U., Duenser C., Moser W.: Efficient calculation of internal results in 2D elasticity BEM. “Engineering Analy- sis with Boundary Elements” 2005, 29, p. 447-453.

7. Fedeliński P., Górski R., Czyż T., Dziatkiewicz G., Ptaszny J.: Analysis of effective properties of materials by using the boundary element method. “Archives of Mechanics” 2014, 66, p. 19-35.

8. Greengard L., Rokhlin V.: A fast algorithm for particle simulations. “Journal of Computational Physics” 1987, 73, p. 325-348.

9. Huang Q. Z., Xu Z. G., Qiang H. F., Wang G., Zheng X. P.: Boundary element method for solid materials with multi- ple types of inclusions. “Acta Mechanica” 2015, 226, p. 547-570.

10. Lei T., Yao Z., Wang H., Wang P.: A parallel fast multipole BEM and its applications to large-scale analysis of 3-D fiber reinforced composites. “Acta Mechanica Sinica” 2006, 22, p. 225-232.

11. Linkov A. M.: Boundary integral equations in elasticity theory. Dordrecht – Boston – London: Kluwer Academic Publishers, 2002.

12. Liu Y.: Fast Multipole boundary element method: theory and applications in engineering. Cambridge: Cambridge University Press, 2009.

13. Liu Y.J., Chen X.L.: Continuum models of carbon nanotube-based composites using the boundary element method.

“Electronic Journal of Boundary Elements” 2003, 1, p. 316-335.

14. Liu Y.J., Li Y.X.: Slow convergence of the BEM with constant elements in solving beam bending problems. “Engineer- ing Analysis with Boundary Elements” 2014, 39, p. 1-4.

15. Liu Y., Nishimura N., Otani Y.: Large-scale modeling of carbon nanotube composites by a fast multipole boundary element method. “Computational Materials Science” 2005, 34, p. 173-187.

16. Ptaszny J., Dziatkiewicz G., Fedeliński P.: Boundary element method modelling of nanocomposites. “International Journal for Multiscale Computational Engineering” 2014, 12, p. 33-34.

17. Ptaszny J., Fedeliński P.: Fast multipole boundary element method for the analysis of plates with many holes. “Ar- chives of Mechanics” 2007, 59, p. 385-401.

18. Ptaszny J., Fedeliński P.: Numerical homogenization by using the fast multipole boundary element method. “Archives of Civil and Mechanical Engineering “ 2011, 11, p. 181-193.

19. Ptaszny J., Fedeliński P.: Numerical homogenization of polymer/clay nanocomposites by the boundary element method. “Archives of Mechanics” 2011, 63, p. 517-532.

20. Rejwer E., Rybarska-Rusinek L., Linkov A.: The complex variable fast multipole boundary element method for the analysis of strongly inhomogeneous media. “Engineering Analysis with Boundary Elements” 2014, 43, p. 105-116.

21. Wang H., Yao Z., Wang P.: On the preconditioners for fast multipole boundary element methods for 2D multi-domain elastostatics. “Engineering Analysis with Boundary Elements” 2005, 29, p. 673-688.

22. Yao Z., Wang H.: Some benchmark problems and basic ideas on the accuracy of boundary element analysis. “Engineer- ing Analysis with Boundary Elements” 2013, 37, 1674-1692.

23. Zohdi T.I., Wriggers P.: An Introduction to Computational Micromechanics. Berlin: Springer-Verlag, 2008.