Acta Sci. Pol. Architectura 15 (1) 2016, 15–26
www.acta.media.pl
LINIOWA ANALIZA NUMERYCZNA ZACHOWANIA SIĉ GRUNTU POD FUNDAMENTEM BEZPOĝREDNIM PRZY WYKORZYSTANIU CZTEROWĉZàOWEGO ELEMENTU SKOēCZONEGO Z WYGàADZONYM POLEM NAPRĉĩEē
Tomasz Pasik
Wydziaá Budownictwa i InĪynierii ĝrodowiska
Szkoáa Gáówna Gospodarstwa Wiejskiego w Warszawie
Streszczenie. W pracy wykonano analizĊ zachowania siĊ gruntu pod fundamentem bezpo- Ğrednim z wykorzystaniem metody elementów skoĔczonych. Kod Ĩródáowy programu zo- staá napisany w jĊzyku Mathematica. Grunt zamodelowano jako materiaá liniowo sprĊĪysty.
Obliczenia zostaáy przeprowadzone dla trzech rodzajów stopy fundamentowej: podatnej gáadkiej, sztywnej szorstkiej i sztywnej gáadkiej. Wyniki przeprowadzonych analiz przy- równano z wynikami obliczeĔ programu Midas GTS NX.
Sáowa kluczowe: metoda elementów skoĔczonych, materiaá liniowo sprĊĪysty, páaski stan odksztaácenia, fundament bezpoĞredni
WSTĉP
Jednym z najwaĪniejszych etapów procesu projektowego w budownictwie jest anali- za statyki konstrukcji. W tym celu najczĊĞciej stosowane są programy komputerowe.
Obecnie powszechnie do problemów inĪynierii geotechnicznej wykorzystywane są narzĊdzia, które wykorzystują metodĊ elementów skoĔczonych (MES), na przykáad Mi- das GTS NX [manual.midasuser.com], GEO5, ZSOIL, PLAXIS, ABAQUS, gdzie zamo- delowaü moĪna zachowanie siĊ konstrukcji (tj. skarpy, nasypu, fundamentu bezpoĞred- niego i poĞredniego, konstrukcji oporowej czy tunelu).
Obsáuga tych programów ze wzglĊdu na wygodĊ sprowadza siĊ do intuicyjnego wprowadzania geometrii, obciąĪeĔ, typu konstrukcji oraz parametrów materiaáu. Na- stĊpnie wyniki są czĊsto interpretowane bez znajomoĞci algorytmów i podstawowych Adres do korespondencji – Corresponding author: Tomasz Pasik, Szkoáa Gáówna Gospodarstwa Wiejskiego, Wydziaá Budownictwa i InĪynierii ĝrodowiska, Katedra GeoinĪynierii,
ul. Nowoursynowska 159, 02-766 Warszawa, e-mail: tomasz_pasik@sggw.pl
© Copyright by Wydawnictwo SGGW, Warszawa 2016
zasad rządzących metodą elementów skoĔczonych. Takie podejĞcie doprowadziü moĪe do rozwiązaĔ nieoptymalnych, a nawet báĊdnych. Skąpe w informacjĊ na temat zasto- sowanych procedur obliczeniowych podrĊczniki uĪytkownika doáączane do programów nie dają moĪliwoĞci gáĊbszego zrozumienie zastosowanych rozwiązaĔ. Pomocne w zro- zumieniu podstaw algorytmu MES są na przykáad opracowania Zienkiewicza i innych [2013], Bhatti [2005 i 2006] oraz Bathe [2014]. Dodatkowo dla moĪliwoĞci praktycznych üwiczeĔ i korzystania z wyĪej wymienionych opracowaĔ niezbĊdne są umiejĊtnoĞci pro- gramowania w dowolnym jĊzyku (np. Mathematica [Wellin 2013], Matlab, Fortran, Py- thon czy C++), wspierane podstawami matematyki wyĪszej, w szczególnoĞci z zakresu algebry liniowej [khanacademy.org].
W wyĪej wymienionych narzĊdziach do analiz MES dostĊp do kodów Ĩródáowych jest utrudniony, a nawet niemoĪliwy, co stanowi barierĊ w testowaniu wáasnych rozwią- zaĔ wykraczających poza standard, który oferują. Buduje to takĪe wraĪenie tzw. czarnej skrzynki, gdzie wprowadzane są dane, i bez rozumienia zasad i procedur obliczeniowych otrzymywane są wyniki.
Przedmiotem niniejszej publikacji jest prezentacja wáasnego programu liniowej ana- lizy MES dla wzorcowego przykáadu inĪynierii geotechnicznej, tj. wyznaczenie rozkáadu naprĊĪeĔ w oĞrodku gruntowym pod wpáywem obciąĪenia fundamentem bezpoĞrednim.
Zaprezentowany program bĊdzie dobrą bazą do rozwijania o kolejne elementy: po- czątkowy stan naprĊĪenia, modele konstytutywne, zagadnienie kontaktu i wiele innych.
MODEL WZORCOWY
Obliczenia przeprowadzono dla páaskiego stanu odksztaácenia (2D). Na rysunku 1 przedstawiono podziaá oĞrodka gruntowego na elementy skoĔczone w dwóch wariantach wraz z numeracją tych elementów i numeracją wĊzáów. Rozpatrywany model ma wy-
Rys. 1. Siatka podziaáu oĞrodka gruntowego na elementy skoĔczone w programie Midas GTS NX: a – wariant 1, b – wariant 2
Fig. 1. FEM Mesh for calculated model of soil halfspace in Midas GTS NX: a – option 1, b – option 2
a b
miary 20 m szerokoĞci i 20 m wysokoĞci. Lewa krawĊdĨ modelu jest jednoczeĞnie osią symetrii ukáadu obciąĪenie – grunt.
Analizy zostaáy przeprowadzone dla 3 rodzajów obciąĪenia na dáugoĞci 6 m (rys. 1):
fundament podatny i gáadki (zamodelowany jako obciąĪenie ciągáe 100 kPa), fundament idealnie sztywny i gáadki (zamodelowany jako wymuszone przemieszcze- nie pionowe o –0,01 m bez blokady przemieszczeĔ poziomych),
fundament idealnie sztywny i szorstki (zamodelowany jako wymuszone przemiesz- czenie pionowe o –0,01 m z blokadą przemieszczeĔ poziomych).
Wszystkie obliczenia przeprowadzono we wáasnym programie napisanym w jĊzy- ku Mathematica [Wellin 2013] (zaáącznik). Dodatkowo otrzymane wyniki przyrównano z wynikami z programu do analiz MES Midas GTS NX.
WyĪej wymienione analizy dla siatki z rysunku 1a, lecz dla elementów 8-wĊzáowych, zostaáy przeprowadzone w publikacji Pottsa i Zdravkovica [1999], co stanowi dodatkowe Ĩródáo porównawcze dla przeprowadzonych obliczeĔ.
PrzyjĊto parametry materiaáu: moduá Younga E = 10 MPa i liczbĊ Poissona v = 0,4.
ZASADA DZIAàANIA PROGRAMU – METODA OBLICZEē
Na rysunku 2 przedstawiono schemat blokowy programu, który uáatwia pracĊ nad ko- dem Ĩródáowym. Stanowi on takĪe gotowe Ĩródáo wiedzy inĪynierskiej dla programisty.
Realizacja kaĪdego bloku skutkuje wynikiem, który jest Ĩródáem danych dla nastĊpnego bloku.
– – –
Rys. 2. Schemat blokowy programu liniowej analizy metody elementów skoĔczonych Fig. 2. Flowchart of linear Finite Element Method
W zaprezentowanym programie siatkĊ MES zaimportowano z programu Midas GTS NX, poprzez import tablicy z wspóárzĊdnymi wĊzáów i tablicy z numerami wĊzáów two- rzących element. Wykorzystano elementy czterowĊzáowe czworoboczne (rys. 3), dla któ- rych funkcje ksztaátu (równanie 1), wyprowadziü moĪna z interpolacji Lagrange’a.
W celu formowania siatki MES w dowolny sposób zdecydowano siĊ na elementy izoparametryczne. DziĊki procedurze mapowania elementu skoĔczonego [Bhatti 2005]
(rys. 4a) do postaci bezwymiarowej (rys. 4c), korzystaü moĪna z elementów o dowolnych wymiarach we wspóárzĊdnych globalnych. PamiĊtaü naleĪy tylko, aby funkcja wypro- wadzona z wyznacznika macierzy Jacobiego przyrównana do 0 nie przecinaáa obsza- ru elementu we wspóárzĊdnych bezwymiarowych (rys. 4b). Taki zabieg sprowadza caáą transformacjĊ ukáadu do Jacobianu.
Na rysunku 4a, c umieszczono takĪe punkty, które sobie odpowiadają w obu ukáadach wspóárzĊdnych. Wystarczy skorzystaü z zaleĪnoĞci:
ܰଵሺߦǡ ߟሻ ൌͳ
Ͷሺͳ െ ߟ െ ߦ ߟߦሻ
ܰଶሺߦǡ ߟሻ ൌͳ
Ͷሺͳ െ ߟ ߦ െ ߟߦሻ
ܰଷሺߦǡ ߟሻ ൌͳ
Ͷሺͳ ߟ ߦ ߟߦሻ
ܰସሺߦǡ ߟሻ ൌͳ
Ͷሺͳ ߟ െ ߦ െ ߟߦሻ
(1)
Rys. 3. Element czworoboczny 4-wĊzáowy i funkcje ksztaátu Fig. 3. 4-noded quadrilateral element and shape functions
Rys. 4. Procedura zmiany wspóárzĊdnych 4-wĊzáowego, czworobocznego elementu skoĔczone- go: a – element we wspóárzĊdnych globalnych, b – linia zerowego Jakobianu, c – element we wspóárzĊdnych bezwymiarowych lokalnych
Fig. 4. Mapping procedure for 4-node quadrilateral element: a – element in global coordinate system, b – zero Jakobian determinant line, c – element in local dimensionless coordinate system
a b c
x(ȟ, Ș) = NTxgc – wspóárzĊdne x w zaleĪnoĞci od ȟ, Ș (rys. 4a, c) y(ȟ, Ș) = NTygc – wspóárzĊdne y w zaleĪnoĞci od ȟ, Ș (rys. 4a, c)
xgc = [xgc1, xgc2, xgc3, xgc4]T – wektor wspóárzĊdnych x wĊzáów (rys. 4a) (2) ygc = [ygc1, ygc2, ygc3, ygc4]T – wektor wspóárzĊdnych y wĊzáów (rys. 4a)
NT = [N1, N2, N3, N4] – wektor funkcji ksztaátu
Zagadnienia inĪynierii geotechnicznej modeluje siĊ w 2D jako páaski stan odksztaáce- nia. Na rysunku 5 zobrazowano zaáoĪenia páaskiego stanu odksztaácenia.
ZaleĪnoĞci na rysunku 2 stanowią podstawĊ liniowego algorytmu MES wyczerpu- jąco przedstawioną w publikacjach: Bhatti [2005], Zienkiewicza i innych [2013], Bathe [2014] i wielu innych.
W praktycznych przypadkach celem przeprowadzanych analiz jest wyznaczenie roz- káadu wybranego skáadnika tensora naprĊĪeĔ z zaleĪnoĞci ı = CBTd.
Dla elementów czterowĊzáowych czworokątnych najbardziej odpowiednim miejscem wewnątrz elementu do obliczania wartoĞci naprĊĪeĔ jest poáoĪenie ȟ = Ș = 0 [Bhatti 2006]. W celu wygenerowania gáadkich wykresów konturowych zastosowano „algorytm wygáadzający wykresy konturowe”, interpolując wartoĞci naprĊĪeĔ z początku ukáadu wspóárzĊdnych elementu do wĊzáów (równania 3 i 4).
Dla elementu:
W(e)s = V(e) (3)
gdzie: ( )e T det( )1 11 1 T ,
A
W ³³N NdA ³ ³ N N J dȟdȘ
1 1 ( )
det( )1 1 ,
e e e
A
v s ³³NdA s ³ ³ N J dȟdȘ
Rys. 5. Warunki páaskiego stanu odksztaácenia dla przykáadu waáu przeciwpowodziowego pod- danego parciu wody i redukcja tensora naprĊĪenia do páaskiego stanu
Fig. 5. Plane strain condition: dike subjected to water pressure and reduction of stress tensor to 2D condition
x xy xz
yx y yz
zx zy z
ı IJ IJ
S IJ ı IJ
IJ IJ ı
ª º
« »
« »
« »
« »
¬ ¼
Páaski Stan Odksztaácenia
İz = 0 Ȗyz = 0 Ȗyz = 0
x xy ,
yx y z
ı IJ IJ ı ı
ª º
« »
« »
¬ ¼
s – wektor wyinterpolowanych wartoĞci w wĊzáach elementu,
se – wyliczona wartoĞü np. naprĊĪenia w początku ukáadu wspóárzĊdnych ȟ = Ș = 0.
Dla caáego ukáadu procedura áączenia elementów jest analogiczna do procedury budo- wy globalnej macierzy sztywnoĞci (kG) i wektora przemieszczeĔ wĊzáowych (rG):
WGS = VG =>S = WG–1VG (4)
gdzie: S – wektor naprĊĪeĔ dla wszystkich wĊzáów w rozpatrywanym modelu.
Na rysunku 6 przedstawiono wykresy konturowe naprĊĪeĔ stycznych (IJxy) wykreĞlone dla róĪnych punktów. Dla rekomendowanej lokalizacji (ȟ = Ș = 0) brak jest informacji na temat znacznego obszaru oĞrodka gruntowego (rys. 6a). Wykres (rys. 6b) jest obarczony báĊdem, po- niewaĪ punktami obliczeĔ byáy wĊzáy z lewej strony elementów (wĊzáy 1 i 4 – rys. 3). Rysunek 6c przedstawia wykres konturowy z zastosowaniem „algorytmu wygáadzającego”. „Algorytm wygáadzający wykresy konturowe” (równania 3 i 4) umoĪliwia porównywanie wykresów konturowych z wykresami generowanymi w programie Midas GTS NX (rys. 8 i 9).
WYNIKI OBLICZEē
W celu sprawdzenia poprawnoĞci programu wykonano analizĊ stopy fundamento- wej w 3 wariantach. Na rysunku 7 przedstawiono przemieszczenia powierzchni terenu.
Otrzymano wyniki identyczne z wynikami analiz przeprowadzanych przez Pottsa i Zdrav- kovica [1999]. Dodatkowo zaobserwowano typ przemieszczeĔ gruntu w zaleĪnoĞci od rodzaju stopy, gdzie decydujący wpáyw ma sztywnoĞü stopy.
Na rysunku 8 i 9 przedstawiono wyniki odpowiednio naprĊĪeĔ stycznych i najwiĊk- szych naprĊĪeĔ gáównych dla kaĪdej ze stóp. Dodatkowo analizy zostaáy przeprowadzo- ne w programie Midas GTS NX dla róĪnych rodzajów siatki oraz elementów. Otrzymano wyniki zbliĪone z wynikami z programu Midas GTS NX dla róĪnych siatek, co Ğwiadczy o poprawnoĞci napisanego programu. Na rysunkach 8 i 9 gáówne róĪnice wynikają ze sztywnoĞci stopy fundamentowej. RóĪnice wynikające z kontaktu miĊdzy stopą funda- mentową a gruntem dla przeprowadzonych analiz są nieznaczne.
Rys. 6. Wykresy konturowe IJxy w oĞrodku gruntowym dla róĪnych punktów (siatka MES 1):
a – dla punktów ȟ = Ș = 0, kontury {5, 10, 15, 20} kPa, b – dla wĊzáów z lewej strony elementów 1 i 4, kontury {5, 10, 15, 20, 25} kPa, c – dla wĊzáów wyinterpolowanych z punktów ȟ = Ș = 0, kontury {5, 10, 15, 20, 25} kPa
Fig. 6. Contour plots of IJxy in soil medium for different data point (FEM mesh 1): a – for location points, contours {5, 10, 15, 20} kPa, b – for nodes located on left side of element 1 and 4, con- tours {5, 10, 15, 20, 25} kPa, c – for nodes using algorithm of smoothed stress ¿ eld, contours {5, 10, 15, 20, 25} kPa
a b c
Rys. 7. Przemieszczenia powierzchni terenu: a – obliczenia wáasne dla elementu 4-wĊzáowego, b – obliczenia Pottsa i Zdravkovica [1999] dla elementów 8-wĊzáowych
Fig. 7. Ground surface settlements: a – own calculation for 4 node elements, b – calculation for 8-node elements [Potts and Zdrawkovic 1999]
a b
Program wáasny – elementy czterowĊzáowe, siatka 1
Midas GTS NX – elementy czterowĊzáowe, siatka 1
Midas GTS NX – elementy czterowĊzáowe, siatka 2
a b c
Midas GTS NX – elementy oĞmiowĊzáowe, siatka 1
Rys. 8. Wykresy konturowe naprĊĪeĔ stycznych (IJxy) w oĞrodku gruntowy[kPa]: a – stopa gáadka podatna, kontury {5, 10, 15, 20, 25} kPa, b – stopa gáadka sztywna, kontury {–1, 0, 1, 2, 3, 4, 5} kPa, c – stopa szorstka sztywna, kontury {0, 1, 2, 3, 4, 5, 6} kPa
Fig. 8. Contour plots of tangential stresses (IJxy) in soil medium [kPa]: a – smooth À exible fo- oting, contours {5, 10, 15, 20, 25} kPa, b – smooth rigid footing, contours {–1, 0, 1, 2, 3, 4, 5} kPa, c – rough rigid footing, contours {0, 1, 2, 3, 4, 5, 6} kPa
Program wáasny – elementy czterowĊzáowe, siatka 1
Midas GTS NX – elementy czterowĊzáowe, siatka 1
a b c
Midas GTS NX – elementy czterowĊzáowe, siatka 2
PODSUMOWANIE
Analiza numeryczna uwzglĊdnia wszystkie cztery teoretyczne wymagania: równowa- gi, ciągáoĞci, relacji konstytutywnych i warunków brzegowych [Potts 2003] (rys. 10a), co jest cenne w przypadku, gdy wymagane są informacje dotyczące na przykáad przemiesz- czenia, statecznoĞci oraz wpáywu na pobliskie konstrukcje. Przewaga analiz MES nad innymi metodami analizy konstrukcji (np. metodami moduáu podatnoĞci podáoĪa [Pasik i Koda 2013, Pasik i in. 2015], polega na moĪliwoĞci opisu oĞrodka gruntowego (sta- nu naprĊĪenia i odksztaácenia), konstrukcji oraz zagadnieĔ kontaktu miĊdzy konstrukcją a gruntem.
AnalizĊ moĪna rozszerzyü o zastosowanie bardziej odpowiedniego modelu konsty- tutywnego [Borja 2013], opisu kontaktu miĊdzy fundamentem a gruntem oraz moĪna Rys. 9. Pola wektorowe i kontury najwiĊkszych naprĊĪeĔ gáównych (ı1): a – stopa gáadka po-
datna, kontury {30, 60, 90} kPa, b – stopa gáadka sztywna, kontury {4, 6, 8, 10, 12} kPa, c – stopa szorstka sztywna, kontury {4, 6, 8, 10, 12} kPa
Fig. 9. Vector ¿ elds and contour plots of major principal stress: a – smooth À exible footing, con- tours {30, 60, 90} kPa, b – smooth rigid footing, contours {4, 6, 8, 10, 12} kPa, c – rought rigid footing, contours {4, 6, 8, 10, 12} kPa
Midas GTS NX – elementy oĞmiowĊzáowe, siatka 1
Rys. 10. ZaleĪnoĞci i procedura analizy liniowej MES: a – podstawowe zaleĪnoĞci i równania mechaniki, b – ogólna procedura otrzymywania wielkoĞci w analizie liniowej MES Fig. 10. Dependencies and procedure of linear FEM analysis: a – basic concept and equations in
mechanics, b – overal procedure to get values in FEM rule
a b
uwzglĊdniü początkowy stan naprĊĪenia. Wszystkie wymienione rozszerzenia przybliĪa- ją do opisu rzeczywistoĞci, lecz znacząco rozszerzają prace a nawet zmieniają procedurĊ algorytmu.
PODZIĉKOWANIA
Praca zostaáa zrealizowana dziĊki wsparciu zadania badawczego w ramach wewnĊtrz- nego trybu konkursowego sáuĪącego rozwojowi máodych naukowców i uczestników stu- diów doktoranckich pt. „Projektowanie geotechniczne z zastosowaniem metody elementów skoĔczonych” ze Ğrodków Ministerstwa Nauki i Szkolnictwa WyĪszego.
PIĝMIENNICTWO
Bathe, K.J. (2014). Finite Element Procedures. K.J. Bathe, Watertown.
Bhatti, M.A. (2005). Fundamental Finite Element Analysis and Applications with Matheamatica and Matlab Computations. John Wiley & Sons, Hoboken.
Bhatti, M.A. (2006). Advanced Topics in Finite Element Analysis of Structures with Matheamatica and Matlab Computations. John Wiley & Sons, Hoboken.
Borja, R.I. (2013). Plasticity Modeling and Computation. Springer, Heidelberg.
Logan, D.L. (2012). A First Course in the Finite Element Method. SI, Cengage Learning, Stamford, USA.
Pasik, T., Koda, E. (2013). Analiza siá wewnĊtrznych i przemieszczeĔ rozpieranej Ğciany szczelino- wej. Acta Sci. Pol. Architectura, 12 (4), 121–133.
Pasik, T., Chalecki, M., Koda, E. (2015). Analysis of embedded retaining wall using the subgrade reaction method. Studia Geotechnica et Mechanica, 37 (1), 59–73.
Potts, D.M, (2003). Numerical analysis: a virtual dream or practical reality. Geotechnique, 53, 6, 535–573.
Potts, D.M., Zdravkovic, L., (1999). Finite element analysis in geotechnical engineering: Theory.
Thomas Telford, London.
Wellin, P. (2013). Programming with Mathematica. An Introduction. Cambridge University Press, New York.
Zienkiewicz, O.C., Taylor, R.L., Zhu, J.Z. (2013). The Finite Element Method: Its Basis and Fun- damentals. Butterworth – Heinemann, Oxford.
ZAàĄCZNIK – KOD ħRÓDàOWY ALGORYTMU WYGàADZAJĄCEGO W JĉZYKU MATHEMATICA
Tabela 1. Oznaczenia de¿ nicji wykorzystanych w algorytmie Table 1. Main designation symbols in algorithm.
Symbol Opis
ș Tablica z kątami odpowiadającymi obrotowi kierunków gáównych tensora naprĊĪenia VW[.] Funkcja zwracająca macierz W oraz wektor V dla elementu na podstawie wspóárzĊdnych
i wartoĞci obliczonej wielkoĞci w początku ukáadu wspóárzĊdnych elementu.
WG Tworzy macierz globalną dla caáego ukáadu VG Tworzy wektor globalny dla caáego ukáadu
șval Wektor kątów obrotu kierunków gáównych dla kaĪdego wĊzáa ukáadu
Algorytm wygáadzający IJxy
LINEAR NUMERICAL ANALYSIS OF SOIL BEHAVIOUR UNDER SHALLOW FOUNDATION USING FOUR NODE QUADRILATERAL FINITE ELEMENT WITH SMOOTHED STRESS FIELD
Abstract. The paper present analysis of soil behavior under shallow foundation using ¿ nite element procedure. The source code have been written in Mathematica language. The soil is modeled as linear elastic material. Computations have been carried out for three types of foot: Smooth À exible, rough rigid and smooth rigid. The results of the performed analysis were compared to the results of calculations by the computer software Midas GTS NX.
Key words: ¿ nite element algorithm, elastic solid, plane strain conditions, shallow founda- tion
Zaakceptowano do druku – Accepted for print: 21.03.2016
Cytowanie: Pasik, T. (2016). Liniowa analiza numeryczna zachowania siĊ gruntu pod fundamen- tem bezpoĞrednim przy wykorzystaniu czterowĊzáowego elementu skoĔczonego z wygáadzonym polem naprĊĪeĔ. Acta Sci. Pol. Architectura 15 (1), 15–26.
