ZASTOSOWANIE FUNKCJI RADIALNYCH W ANALIZIE STATYCZNEJ BELKI

ZASTOSOWANIE FUNKCJI RADIALNYCH W ANALIZIE STATYCZNEJ BELKI

Leszek Majkut

^1a

1AGH Akademia Górniczo – Hutnicza, Katedra Mechaniki i Wibroakustyki

amajkut@agh.edu.pl

Streszczenie

W pracy opisano bezsiatkową metodę kolokacyjną Kansy i jej zastosowanie w analizie statycznej belki. W anali- zie zastosowano funkcję wielokwadratową, opisano metodę doboru wartości parametru kształtu oraz podano współczynniki korygujące wyniki analizy w przypadku wymuszenia punktowego. Wszystkie wyniki porównano z wynikami analitycznymi z zastosowaniem siedmiu różnych błędów względnych pozwalających ocenić jakość aproksymacji.

Słowa kluczowe: metody bezsiatkowe, linia ugięcia, aproksymacja

STATIC ANALYSIS OF THE BEAM LIKE STRUCTURES WITH RADIAL BASED FUNCTION

Summary

The work concerns the static analysis of the Bernoulli-Euler beam with the Radial Based Functions. The Kansa collocation method was used for determination deflection, slope, bending moment and shear force of the beam. All results were compared with analytical ones. All results indicate that using of multiquadric (MQ) RBF provide a results with very high accuracy in comparison to analytical results in static analysis of beam-like structural components.

Keywords: meshfree methods, beam deflection, approximation

1. WSTĘP

Modelowanie matematyczne układów mechanicznych prowadzi zazwyczaj do układu równań różniczkowych zwyczajnych lub cząstkowych, który łącznie z warunka- mi brzegowymi stanowi problem początkowo-brzegowy.

Tradycyjne metody numeryczne rozwiązania takiego problemu wymagają podziału na elementy: całego anali- zowanego obszaru, w metodach takich jak metoda różnic skończonych (MRS), metoda elementów skończonych (MES) czy metoda objętości skończonych (MOS), lub jedynie brzegu obszaru w metodach takich jak metoda elementów brzegowych (MEB). Olbrzymi nakład pracy włożony w rozwój MES uczynił z niej najbardziej popu- larną, a nawet podstawową, metodę przybliżonego rozwiązania problemów początkowo-brzegowych.

Największy problem z zastosowaniem metod wyma- gających siatki elementów to definicja właśnie siatki.

Generacja modelu składającego się z kilku czy kilkunastu tysięcy elementów o różnych kształtach, wielkościach

i orientacji nie jest łatwa i jest czasochłonna. Ze względu na obniżające się koszty obliczeń komputerowych gene- racja siatki staje się głównym składnikiem kosztów analizy numerycznej układów mechanicznych.

Metody przybliżonego rozwiązania problemu począt- kowo-brzegowego niewymagające siatki elementów to tzw. metody bezsiatkowe. Metody te podzielić można na dwa typy: metody obszarowe, m.in. bezelementowa metoda Galerkina [1], metoda interpolacji punktowej [12], metoda funkcji Rvacheva [13] i metody brzegowe, do których należą m.in. metoda węzłów brzegowych [8], metoda parametrycznych układów równań całkowych (PURC) [16]. W obu typach metod punkty kolokacyjne umieszczone są w całym analizowanym obszarze lub jedynie na jego brzegu. Większość wymienionych metod jest bezsiatkowa jedynie w sensie interpolacji rozwiąza- nia, niezbędne są natomiast „komórki” wykorzystywane

do całkowania słabego sformułowania wariacyjnego problemu.

Innym podejściem poszukiwania przybliżonego roz- wiązania problemu początkowo-brzegowego jest metoda Kansy [6], bazująca na sformułowaniu silnym i radial- nych funkcjach bazowych (RFB). Metoda ta nie wymaga zastosowania jakichkolwiek elementów czy komórek, dlatego też zwana jest metodą „prawdziwie” bezsiatko- wą.

W ciągu ostatnich 20 lat kolokacyjna metoda Kansy znalazła zastosowanie do rozwiązania wielu różnych problemów mechaniki, m. in. analizy dwu- i trójfazowych modeli tkanki [4], symulacji prądów i pływów wód płytkich [5], problemów wymiany ciepła [15], analizy równań Naviera-Stokesa [2], analizy pola elektromagne- tycznego [11], problemów związanych z optymalizacją [9] i wielu innych.

Sukces metod związanych z zastosowaniem RFB wią- że się z ich podstawową własnością polegającą na trans- formacji problemów wielowymiarowych do problemów jednowymiarowych. W niniejszej pracy zastosowano kolokacyjną metodę Kansy wraz z wielokwadratową (ang. multiquadric) radialną funkcją bazową do analizy statycznej belki. W celu wykazania poprawności i do- kładności proponowanej metody obliczeń wszystkie wyniki analizy przybliżonej porównano z wynikami uzyskanymi metodą analityczną.

Przy stosowaniu funkcji wielokwadratowej bardzo istotne jest wyznaczenie tzw. parametru kształtu, od którego zależy kształt funkcji bazowej. W pracy zapro- ponowano metodę wyboru „optymalnego” parametru kształtu, która związana jest z analizą aproksymacyjną funkcji wymuszenia.

Problem ugięcia belki był już przedmiotem analizy np. [10], natomiast według wiedzy autora niniejsza praca jest jedyną, w której analizowane jest wymuszenie siłą skupioną.

2. RADIALNE FUNKCJE BAZOWE

Funkcja radialna to każda funkcja jednej zmiennej postaci:

( ) (

^j

)

j r =φ x x

φ − (1)

gdzie: x−x_j jest Euklidesową odległością pomiędzy punktami x (współrzędna niezależna funkcji radialnej (1)) i xj. Punkt xj jest nazywany jest centrum funkcji radialnej (1). Zmieniając położenia centrów, otrzymuje się rodzinę funkcji, która tworzy bazę stosowaną do interpolacji lub aproksymacji dowolnej funkcji.

Każdą z funkcji bazowych postaci (1) zaliczyć można do jednej z dwóch kategorii: funkcje o zwartym nośniku, tj. funkcje, które są różne od zera jedynie w sferze

o promieniu r (najczęściej r = 1) lub funkcje o nośniku nieograniczonym (r → ∞∞∞∞).

Najczęściej stosowanymi funkcjami z pierwszej kate- gorii są, zamieszczone w tabeli 1, funkcje Wendlanda [14].

Tab. 1. Bazowe funkcje radialne o zwartym nośniku wymiar przestrzeni definicja funkcji

d = 1

( )

r+

= φ(r) 1−

( ) (

¹

^r

^{3 3r}

⁺

¹

)

=

φ(r) − ₊

( )

¹

^r

⁵

(

^8r2

⁺

^5r

⁺

¹

)

=

φ(r) − ₊

d = 2,3

( )

¹

r

+²

=

φ(r) −

( ) ( ^{1 r}

⁴

^{4r + 1} )

=

φ(r) −

₊

( )

¹

^r

⁶

(

³⁵

^r

²

⁺

¹⁸

^{r +}

³

)

=

φ(r) − ₊

gdzie:

( ) ( )





 − ∈

− 0, dla 1

0,1 dla

1 1

>

r r ,

= r r

=

φ(r)

n n +

W tabeli 2 zawarto przykłady funkcji o nośniku nie- ograniczonym.

Tab. 2.

Bazowe funkcje radialne o nośniku nieograniczonym definicja funkcji nazwa funkcji

r r

)=

φ( liniowa

) 3

(

r

=

r

φ sześcienna

2 2

+ c r

=

φ (r)

wielokwadratowa

r r

= (r)

²ln

φ cienkiej płyty

2 2

1 c + r

=

φ(r) wielokwadratowa od-

wrotna

cr2

e

= (r)

⁻

φ

^Gaussa

Ze względu na dużą popularność w zastosowaniach oraz dobre własności aproksymujące w pracy zastosowa- no funkcję wielokwadratową, która w przypadku 1D ma postać:

( ) ^{r =} ( ^{x x}

j

)

²

^{+ c}

²

j −

φ (2)

Jak wspomniano wcześniej, kształt funkcji (2) zależy od parametru kształtu c. Wraz ze wzrostem jego warto- ści funkcja (2) staje się coraz bardziej płaska, przez to mało wrażliwa na zmiany odległości pomiędzy punktem x i centrum funkcji radialnej xj. Problem doboru wartości parametru c jest jak dotąd nierozwiązany.

3. METODA KOLOKACYJNA KANSY

Analizowany problem początkowo brzegowy opisany jest równaniem postaci:

( ) ^{x , x}

^Ω

f

=

Lu

∈ (3)

wraz z warunkami brzegowymi postaci:

( ) ^{x , x}

^Γ

g

=

Bu

∈ (4)

gdzie: L jest liniowym operatorem różniczkowym, B jest operatorem opisującym warunki brzegowe, Ω to analizo- wany obszar, Γ to brzeg tego obszaru.

Idea metody Kansy polega na aproksymacji rozwią- zania problemu początkowo-brzegowego (3) i (4) za pomocą sumy szeregu rodziny funkcji radialnych tj.:

^{u =}

^α j

( ) ^r

N

j=

jφ

∑

1

ˆ ₍₅₎

Współczynniki αj wyznacza się w procedurze koloka- cji. W tym celu należy wybrać zbiór No punktów {x1, x2,

…, xNo} należących do obszaru Ω, w których żąda się, by przybliżone rozwiązanie (5) spełniało równanie (3):

( )

^{( )}

) ˆ (

1

i i j N

j=

j

i α

L r = f x

x f

= u

L

^⇔

∑

φ ⁽⁶⁾

Podobnie należy wybrać zbiór Nb punktów {xNo+1, xNo+2, …, xNo+Nb} na brzegu Γ analizowanego obszaru.

W tych punktach muszą zostać spełnione równania warunków brzegowych (4):

( )

^{( )}

1

i i j N

j=

j

B r = g x

α φ

∑

⁽⁷⁾

Równania (6) i (7) stanowią liniowy układ równań, który zapisać można w postaci macierzowej:

f

Aα= (8)

Z układu (8) wyznaczyć można poszukiwane współ- czynniki α^j. W przypadku, gdy suma liczb wybranych punktów obszaru Ω tj. No i liczby punktów brzegowych Nb jest równa liczbie punktów centralnych N (No+ Nb=N) układ (8) rozwiązać można, stosując eliminację Gaussa, w przypadku gdy No+ Nb > N układ (8) jest układem nadokreślonym i należy poszukiwać rozwiązania metodą najmniejszych kwadratów, np. stosując rozkład SVD.

4. ANALIZA STATYCZNA BELKI

W pracy przez analizę statyczną rozumiane jest wy- znaczenie linii ugięcia, kątów obrotu przekrojów (stycz- nej do linii ugięcia), przebiegu momentu gnącego i siły

tnącej belki opisanej modelem Bernoullego-Eulera. Ze względu na to, że celem pracy jest wykazanie poprawno- ści i dokładności metody bezsiatkowej, wszystkie opisane powyżej wielkości zostały wyznaczone w 200 punktach belki metodą przybliżoną i metodą analityczną.

Linia ugięcia w(x) belki opisana jest równaniem:

) ) (

(

4 0 4

x f dx

x w

EI d

= (9)

gdzie: EI to sztywność na zginanie, f0(x) to obciążenie zewnętrzne na jednostkę długości belki

Funkcja w(x) z równania (9) musi dodatkowo speł- niać warunki brzegowe. W pracy analizowano typowe warunki brzegowe belki, tj. utwierdzenie (w = 0 i w' = 0), swobodne podparcie (w = 0 i w'' = 0) i wolny koniec (w'' = 0 i w''' = 0).

4.1 ROZWIĄZANIE ANALITYCZNE

Funkcje linii ugięcia, stycznej, momentu gnącego i siły poprzecznej są opisane przez:

(

^{4( ) 1 3/6 2 2/2 3 4}

)

/ 1 )

(

x EI f x C x C x C x C

w

= + + + +

- linia ugięcia

(

2 2 3

)

1

3( ) /2

/ 1 ) (

'

x EI f x C x C x C

w

= + + + - styczna

2 1 2( ) ) (

''

x f x C x C w

EI

= + + - moment gnący

1 1( ) ) ( ''

'

x f x C

w

EI

= + - siła tnąca

w równaniach ^{fi− x =}

∫

^{l fi x}

0

1( ) ( ), i = 1,2,3,4. Stałe C1-C4

zależą od warunków brzegowych.

4.2 ROZWIĄZANIE PRZYBLIŻONE

Zgodnie z metodą Kansy równanie linii ugięcia belki (9) może być zapisane w postaci:

( )

0^{( )}

2 2 4

4

1

i x

x j

N

j=

j

x x c f x

dx

α

d EI

i

= +

− =

∑

⁽¹⁰⁾

gdzie: xi to współrzędne punktów kolokacyjnych i = 1,2,..., No, xj to współrzędne centrów wielkokwadra- towych funkcji radialnych xj = 1,2,...N.

We wszystkich analizowanych w pracy przypadkach przyjęto jedną wartość parametru kształtu dla wszyst- kich funkcji bazowych, tj. cj = c. Punkty kolokacyjne {xi}, w każdym z analizowanych przypadków były równomiernie rozłożone wzdłuż długości belki.

Dla każdego i równanie (10) opisuje i-ty wiersz ma- cierzy A z równania (8), 4 ostatnie wiersze tej macierzy zależą od równań opisujących warunki brzegowe.

Po wyznaczeniu współczynników α^j, możliwe jest wy- znaczenie linii ugięcia, stycznej, momentu gnącego i siły poprzecznej, które opisane są odpowiednio funkcjami:

( )

^{2 2}

1

)

(

x x x c

w

_j

N

j

j − +

=

∑

=

α - linia ugięcia

( )

^{2 2}

1

) (

'

x x c

dx x d

w

_j

N

j

j − +

=

∑

=

α ^{- styczna}

( )

^{2 2}

2 2

1

) ( '

'

x x c

dx x d

w

EI

_j

N

j

j − +

=

∑

=

α - moment gnący

( )

^{2 2}

3 3

1

) ( ''

'

x x c

dx x d

w

EI

_j

N

j

j − +

=

∑

=

α - siła tnąca

W przypadku wymuszenia punktowego każda z po- wyższych funkcji musi zostać pomnożona przez współ- czynnik korekcyjny, równy 1/d, gdzie d jest średnią odległością pomiędzy punktami kolokacyjnymi.

Dobór „optymalnej” wartości parametru kształtu jest, jak dotąd, nierozwiązanym problemem. Szczegółowy przegląd różnych strategii pozwalających na wyznaczenie tego parametru znaleźć można w pracy [3].

W niniejszej pracy autor proponuje inną metodę opartą na wyznaczeniu aproksymacji funkcji wymuszenia f0(x) z równania (10), wraz z uwzględnieniem warunków brzegowych, dla różnych c. Za wartość „optymalną”, w proponowanej metodzie przyjmuje się tę wartość parametru c, dla której błąd zdefiniowany równaniem (11) osiąga minimum.

( )

∫ ∑

= = 









 − + −

=

l

x N

j

i j

j x x c f x EI dx

dx RMS d

0

2

1

2 2 4

4

/ )

α ( ^{(11) (11)}

Na rys. 1 pokazano przykład zmian błędu (11) w funkcji parametru kształtu

Rys. 1. Zmiana wartości błędu aproksymacji (11) w funk- cji parametru kształtu c wyznaczone dla belki swobodnie

podpartej z obciążeniem opisanym funkcją liniową Górnym ograniczeniem wartości c jest wartość wskaźnika uwarunkowania macierzy aproksymacji.

Standardowe biblioteki algebry liniowej prowadzą do niestabilności rozwiązania, gdy wskaźnik uwarunkowania osiąga wartość rzędu 10¹⁷ [3]. Zmiany wartości wskaźnika

uwarunkowania macierzy aproksymacyjnej w funkcji parametru kształtu pokazano na rys. 2

Rys. 2. Zmiana wartości wskaźnika uwarunkowania w funkcji parametru kształtu c wyznaczone dla belki swo- bodnie podpartej z obciążeniem opisanym funkcją liniową

4.3 PORÓWNANIE WYNIKÓW ANALIZY

W celu oceny „jakości” wyników aproksymacji autor zdefiniował następujące błędy względne (subskrypt a oznacza wyniki analizy analitycznej, subskrypt p – analizy przybliżonej, każdą z funkcji wyznaczono w = 200 punktów):

% max 100

max max

1 − ⋅

=

a p a

w w w

E

,

% max 100

max max

2 ⋅

′′

− ′′

= ′′

a p a

w w w

E

,

Błędy E1 i E2 związane są z maksymalnymi warto- ściami odpowiednio strzałki ugięcia i maksymalnego momentu gnącego.

( )

n x w x w E

n

i

i p i

∑

a

−

−

= ¹

2

3

) ( ) (

,

( )

( )

% 100

0 2 0

2

4 ⋅

−

=

∫

∫

dx w

dx w w

E l

a l

p a

,

( ) ( )

( )

^100%

0

2 2 0

2 2

5 ⋅

+ ′

− ′ + ′

−

=

∫

∫

dx w w

dx w w w w

E l

a a l

p a p a

,

0.2 0.25 0.3 0.35 0.4 0.45 0.5

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

parametr kształtu

RMS

0.2 0.25 0.3 0.35 0.4 0.45 0.5

0 2 4 6 8 10 12 14x 10¹⁷

parametr kształtu

wskaźnik uwarunkowania

( )

^100%

0

2 2 2

0 1

6 ⋅

+ ′′

+ ′

=

∫

∫

dx w w w

dx L

E

_l

a a a

l

,

( ) (

²

) (

²

)

²

1 wa wp wa wp wa wp

L = − + ′ − ′ + ′′− ′′

( )

^100%

0

2 2 2 2

0 2

7 ⋅

+ ′′′

+ ′′

+ ′

=

∫

∫

dx w w w w

dx L

E

_l

a a a a

l

gdzie:

( ) (

²

) (

²

) (

²

)

²

2

w

a

w

p

w

a

w

p

w

a

w

p

w

a

w

p

L

= − + ′ − ′ + ′′− ′′ + ′′′− ′′′

Błędy E3-E7 pozwalają na określenie jakości aprok- symacji w porównaniu do wyników analizy analitycznej.

Analizy błędów przeprowadzono dla belki o czterech różnych warunkach brzegowych:

• swobodnie podparta - swobodnie podparta (s-s),

• utwierdzona - swobodnie podparta (u-s),

• utwierdzona – wolny koniec (u-w),

• utwierdzona – utwierdzona (u-u)

i dla sześciu różnych wymuszeń:

• stałe

f

01=

q

,

• liniowe

f

02 =

( q

2 −

q

1

)

/

l

⋅

x

+

q

1,

• sinusoidalne

^f

03 =

^e

⋅^sin

(

π

^x

^/

^l )

,

• paraboliczne

f

04 =

a

(

x

−

l

/2)²+

bx

+

q

,

• pierwiastkowe

f

₀₅=

d x

+

q

,

• skupione

f

₀₆ =

P

δ(

x

,

l

/2).

W obliczeniach przyjęto dane: EI = 1000 N/m², l = 2.0 m, q = q1 = -100 N/m, e = -100 N/m, a = -50 N/m³, b = 25 N/m², d = -75 N/m^3/2, P = -100 N.

Na rys.3 pokazano przykładowe wyniki analizy sta- tycznej belki u-w obciążonej funkcją

f

₀₃. Linią ciągłą oznaczono wyniki analizy przybliżonej znakami x – wyniki analizy analitycznej.

Wyniki analizy – wartości błędów E1 - E7, wyznaczo- ne dla wszystkich przypadków warunków brzegowych i funkcji opisujących obciążenie zamieszczono w tab. 3.

Rys. 3. Przykładowe wyniki analizy statycznej belki wsporniko- wej, linią ciągłą oznaczono wyniki analizy przybliżonej, x –

analiza analityczna

Dodatkowo, w celu oceny jakości aproksymacji wy- znaczono współczynniki korelacji Pearsona oddzielnie dla funkcji ugięcia (K1), stycznej (K2), momentu gnącego (K3) i siły poprzecznej (K4). We wszystkich przypadkach warunków brzegowych i obciążeń ciągłych, wszystkie (każdy z 80) współczynniki korelacji są równe 1.0.

W przypadkach obciążenia punktowego współczynniki K1, K2 i K3 są równe 1.0. Współczynniki K4 (pomiędzy funkcjami opisującymi funkcje siły poprzecznej) mają wartości 0.997-0.999. Przyczyną tego nie jest błąd w aproksymacji, a efekt Gibbsa – siła tnąca opisana jest funkcją nieciągłą.

0 0.5 1 1.5 2

-0.15 -0.1 -0.05 0 0.05

ugięcie

x[m]

0 0.5 1 1.5 2

-0.08 -0.06 -0.04 -0.02 0

styczna

x[m]

0 0.5 1 1.5 2

-150 -100 -50 0 50

moment gnący

x[m]

0 0.5 1 1.5 2

-50 0 50 100 150

siła poprzeczna

x[m]

Tab. 3. Błędy analizy przybliżonej belki warunki

brzegowe

wymu-

szenie E1 [%] E2 [%] E3 E4 [%] E5 [%] E6 [%] E7 [ %]

p-p

f01 0.0012 0.0049 0.0004 0.0031 0.0026 0.0024 0.0037 f02 0.0031 0.0016 0.0018 0.0061 0.0045 0.0035 0.0041 f03 0.0018 0.0005 0.0004 0.0032 0.0022 0.0019 0.0044 f04 0.0016 0.0025 0.0005 0.0040 0.0033 0.0030 0.0040 f05 0.0014 0.0025 0.0008 0.0034 0.0027 0.0026 0.0050 f06 0.1141 0.7707 0.0139 0.1194 0.1210 0.1840 1.8502

u-p

f01 0.0030 0.0042 0.0006 0.0109 0.0089 0.0064 0.0148 f02 0.0008 0.0041 0.0008 0.0079 0.0077 0.0057 0.0160 f03 0.0023 0.0045 0.0003 0.0066 0.0058 0.0048 0.0135 f04 0.0133 0.0065 0.0011 0.0247 0.0155 0.0103 0.0156 f05 0.0078 0.0051 0.0006 0.0143 0.0092 0.0065 0.0100 f06 0.0416 0.1904 0.0022 0.0456 0.0539 0.2358 0.3733

u-w

f01 0.0160 0.0001 0.0007 0.0251 0.0131 0.0049 0.0138 f02 0.0176 0.0039 0.0016 0.0303 0.0172 0.0081 0.0143 f03 0.0118 0.0005 0.0004 0.0167 0.0093 0.0040 0.0121 f04 0.0223 0.0025 0.0008 0.0391 0.0215 0.0089 0.0148 f05 0.0206 0.0014 0.0015 0.0334 0.0171 0.0066 0.0130 f06 0.4688 0.3631 0.1943 0.4286 0.4627 0.4559 1.9055

u-u

f01 0.0757 0.0571 0.0730 0.0729 0.0750 0.0748 0.0680 f02 0.0780 0.0624 0.1852 0.0752 0.0774 0.0772 0.0720 f03 0.0649 0.0454 0.0369 0.0615 0.0640 0.0634 0.0554 f04 0.0769 0.0571 0.0607 0.0738 0.0762 0.0760 0.0680 f05 0.0761 0.0590 0.0392 0.0734 0.0754 0.0751 0.0691 f06 0.0015 0.2061 0.0061 0.0031 0.0071 0.2524 0.9122

5. PODSUMOWANIE I WNIOSKI

W pracy opisano metodę kolokacyjną Kansy oraz jej zastosowanie do analizy statycznej belki. W celu oceny poprawności i dokładności tej metody analizy porównano przebiegi ugięcia, pochodnej, momentu gnącego i siły tnącej z przebiegami wyznaczonymi metodą analityczną.

Analizę taką przeprowadzono dla belki o czterech

różnych warunkach brzegowych i przy sześciu różnych wymuszeniach.

Wyniki analizy (wartości siedmiu różnych błędów względnych) zamieszone w tab.3 wskazują jednoznacznie na bardzo dużą dokładność analizy z zastosowaniem tej metody bezsiatkowej. Wyniki analizy wskazują również na poprawność zaproponowanej metody doboru parame- tru kształtu oraz współczynnika korygującego w przy- padku wymuszenia punktowego.

Literatura

1. Belystcho T., Lu Y., Gu L.: Element free Galerkin methods. “International Journal for Numerical Methods in Engineering” 1994, 37, p. 229-256.

2. Chinchapatnam, P.P., Djidjeli, K., Nair, P.B.: Radial basis function meshless method for the steady incompressible Navier–Stokes equations. “International Journal for Numerical Methods in Engineering” 2007, 84, p. 1509-1526.

3. Fasshauer, G.E., Zhang, J.G.: On choosing ,,optimal'' shape parameter for RBF approximation. “Numerical Algorithms” 2007, 45, p. 346-368.

4. Hon Y.C., Lu M.W., Xue W.M., Zhou X.:A new formulation and computation of the triphasic model for mechano-electrochemical mixtures. “Computational Mechanics” 1999, 24, p. 155-165.

5. Hon Y.C., Cheung K.F., Mao X.Z., Kansa E.J.: Multiquadric solution for shallow water equations. “ASCE Journal of Hydraulic Engineering” 1999, 125, p. 524-533.

6. Kansa E.J.: Multiquadric-a scattered data approximation scheme with applications to computational fluid dynamics. “Computers & Mathematics with Applications” 1990, 19, p. 147-161.

7. Mui-Duy N., Tanner R.I.: Computing non-Newtonian fluid flow with radial basis function networks.

“International Journal for Numerical Methods in Fluids” 2005, 48, p. 1309-1336.

8. Mukherjee Y.X., Mukherjee S.: The boundary node method for potential problems. “International Journal for Numerical Methods in Engineering” 1997, 40, p. 797-815.

9. Pearson, J.W..: A radial basis function method for solving PDE-constrained optimization problems. “Numerical Algorithms” December 2012 DOI 10.1007/s11075-012-9675-6 (Article not assigned to an issue - Online First).

10. Tiago C.M., Leitao V.M.A.: Application of radial basis functions to linear and nonlinear structural analysis problems. “Computers and Mathematics with Applications” 2006, 51, p.1311-1334.

11. Vu, P., Fasshauer, G.E.: Application of two radial basis function-pseudospectral meshfree methods to three- dimensional electromagnetic problems. “IET Science, Measurement & Technology” 2011, 5, p. 206-210.

12. Wang J.G., Liu G.R.: A point interpolation meshless method based on radial basis functions. Int. Journal for Numerical Methods in Engineering 54, 2002, pp. 1623-1648

13. Wawrzynek A., Detka M., Cichoń, Cz.: Zastosowanie metody R-funkcji do wyznaczania współczynnika przejmowania ciepła. „Modelowanie Inżynierskie” 2012, nr 43, s. 255-263.

14. Wendland, H.: Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. “Advances in Computational Mathematics” 1995, 4, p. 389-396.

15. Zerroukat M., Power H., Chen C.S.: A numerical method for heat transfer problem using collocation and radial basis functions. “International Journal for Numerical Methods in Engineering” 1998, 42, p. 1263-1278

16. Zieniuk E., Sawicki D.: Metoda PURC w analizie nieustalonego pola temperatury w obszarach płaskich.

„Modelowanie Inżynierskie” 2012, nr 44, s. 285-292