POROVNÁNÍ VÝSLEDKŮ NUMERICKÉHO MODELOVÁNÍ S KLASICKÝMI TEORIEMI, CHYBY PŘI MODELOVÁNÍ
Eva Hrubešová Lukáš Ďuriš Josef Aldorf Taťána Petrášová Fakulta stavební, VŠB- Technická univerzita Ostrava TunelářskéTunelářské odpoledne , 26.3.2014, odpoledne březen 2014 -Hotel Praha Duo Praha
Osnova přednášky 1. Úvod 2. Přehled metod pro stanovení zatížení výztuže tunelu, porovnání výsledků 3. Formulační chybové aspekty numerických modelů 4. Diskretizační chybové aspekty numerických modelů 5. Numerické chybové aspekty numerických modelů
Tunelářské odpoledne březen 2014 - Praha
Přehled metod pro stanovení zatížení výztuže tunelu Základní odlišnost od klasické stavební mechaniky: zatížení výztuže podzemních konstrukcí není explicitně známo, ale vyplývá až ze vzájemné spolupráce celého systému „ horninový masív- konstrukce“. Základní členění metod pro výpočet zatížení: klasické metody : klenbové metody – zatížení výztuže je dáno tíhou horniny pod klenbou přirozené rovnováhy (Protodjakonov, Mostkov, …) metody zohledňující snížení tíhy sloupce zeminy nad stropem tunelu v důsledku třecích sil od bočního tlaku, vycházejí z rovnováhy sil metody založené na rovnicích teorie mechaniky kontinua ( silové, deformační, analytické, numerické)
Model je vždy více či méně přesná aproximace skutečného stavu. Tunelářské odpoledne březen 2014 - Praha
Klasické metody stanovení zatížení výztuže Základní obecná charakteristika metody silové, nezohledňují deformační chování výztuže, výztuž je pouze pasivním prvkem základní vztahy jsou formulovány pro homogenní , izotropní nezohledňují tvar díla, pouze jeho velikost
Tunelářské odpoledne březen 2014 - Praha
Klenbová Protodjakonova metoda
Polovina šířky klenby: a 1 = a + h (cot gα + cot g (45 + ϕ / 2 )) pro f 〈 40
Svislé zatížení:
a 1 = a pro f 〉 = 40
qv = b γ
ϕ = arctg f
Výška klenby: b=
a1 f
γ- objemová tíha prostředí (výšku nadloží metoda nezohledňuje !) Tunelářské odpoledne březen 2014 - Praha
Metoda dle Bierbäumera Snižuje účinek celé tíhy nadloží G o účinek tření T, které vzniká podél sloupce horniny nad klenbou díla T
G
T
h
N
N
Tíha sloupce nad výrubem:
G = γ . h .b
b
Rovnováha sil:
Aktivní tlak od klínu zeminy
1 2 2 ϕ N = γ .h .tg 45 − 2 2 Tření:
Q = G − 2.T Zatížení stropu:
Q pv = b
T = N .tgϕ Tunelářské odpoledne březen 2014 - Praha
Metody vycházející z teorie mechaniky kontinua
vycházejí ze základních rovnic mechaniky kontinua (diferenciální rovnice rovnováhy, rovnice kompatibility, …)-silové nebo deformační metody analytické (přímé řešení diferenciálních rovnic) nebo numerické (přibližná řešení, převádějí řešení diferenciálních rovnic na řešení soustav lineárních algebraických rovnic) silové metody – nezohledňují spolupráci s výztuží deformační metody - zohledňují spolupráci s výztuží
Tunelářské odpoledne březen 2014 - Praha
Silová metoda vycházející z mechaniky kontinua metoda dle Terzaghiho Využívá Janssenovu rovnici pro zatížení v silu, nezohledňuje spolupráci s výztuží.
Vertikální zatížení výztuže díla:
γb qv = 2 tan ϕ K b
2h − K b tan ϕ b 1 − e
Kb- koeficient bočního tlaku Kolymbas, D. Tunneling and Tunnel Mechanics. Springer 2005
Tunelářské odpoledne březen 2014 - Praha
Deformační metody mechaniky kontinua pro stanovení zatížení zatížení výztuže je dáno průsečíkem Fenner-Pacherovy křivky horniny a pracovně-deformační charakteristiky výztuže – bod rovnováhy systému závislost zatížení na deformacích proběhlých před instalací výztuže
závislost zatížení na tuhosti výztuže
Tunelářské odpoledne březen 2014 - Praha
Analytická deformační metoda dle Bulyčeva řeší rovnice mechaniky kontinua analytickými prostředky – využívá teorie analytických funkcí komplexní proměnné, vztahů Kolosova-Muschelišviliho a komplexních potenciálů metoda vychází z rovnosti posunů na kontaktu horninového prostředí a výztuže řešení předpokládá pružné, izotropní, homogenní horninové prostředí výztuž je uzavřená, má konstantní tloušťku a konstantní materiál po celém obvodu metoda předpokládá střed díla v hloubce rovné minimálně 5-ti násobku poloměru díla (v této hloubce lze již zanedbat vliv okrajové podmínky na povrchu ) – převod úlohy z pružné poloroviny na řešení úlohy v pružné rovině, primární napětí je v modelu konstantní a odpovídá hloubce díla pod povrchem lze řešit i jiná než kruhová podzemní díla, jedinou podmínkou je alespoň jedna osa symetrie relativně složitý matematický aparát, ale konečné řešení je získáno v uzavřeném tvaru, dobře algoritmizovatelné a programovatelné, časově nenáročný výpočet Tunelářské odpoledne březen 2014 - Praha
Numerické metody vycházející z teorie mechaniky kontinua metoda konečných prvků, hraničních prvků, konečných diferencí univerzální metody umožňují zohlednit širokou geometrickou i materiálovou variabilitu, různé konstitutivní modely chování materiálu umožňují implementovat do modelu různé typy výztužních prvků umožňují zohlednit časové fáze výstavby díla Prezentované srovnávací výpočty provedeny softwarem PLAXIS 2D ( předpoklad: výztuž instalována ihned po vyražení, nebyla aplikována β –metoda pro zohlednění časového odstupu mezi vyražením díla a jeho vyztužením) Tunelářské odpoledne březen 2014 - Praha
Řešená úloha pro srovnání výsledků metod pro stanovení zatížení výztuže Tvar díla : kruh o poloměru 5 m Výška nadloží: variantní 5-30 m Tloušťka výztuže: variantně 5 cm, 20 cm Materiál výztuže: Eb=20 000 MPa, µ=0.2, γ=24 kN/m3 Zeminové prostředí: homogenní, izotropní, jílovité prostředí Objemová tíha zemin v okolí tunelu: 20 kN/m3 Modul pružnosti prostředí: 10 MPa Poissonovo číslo: 0.4 Úhel vnitřního tření zeminového prostředí: 25 °
? zatížení stropu díla
Soudržnost: 20 kPa Tunelářské odpoledne březen 2014 - Praha
Srovnání hodnot zatížení výztuže díla ve stropě stanoveného různými metodami
700 zatížení stropu díla (kPa)
600 500 400 300 200
výška klenby dle Protodjakonova
100 0 0
5
10
15
Výška nadloží (m)
20
25
30
35
minim. hloubka platnosti výpočtu dle Bulyčeva
tíha celého nadloží Terzaghi Analytický (Bulyčev), tloušťka = 5cm Plaxis (d=5 cm)
Protodjakonov Bierbaumer Analytický (Bulyčev), tloušťka d = 20cm Plaxis (d=20 cm)
Tunelářské odpoledne březen 2014 - Praha
Srovnání hodnot zatížení výztuže díla ve stropě stanoveného různými metodami (detail pro hloubky 5-20 m) zatížení stropu díla (kPa)
390
290
190
90 5
6
7
8
9
10
11
12
13
14
15
16
17
18
výška nadloží (m)
tíha celého nadloží
Terzaghi
Plaxis (d=5 cm)
Plaxis (d=20 cm)
Tunelářské odpoledne březen 2014 - Praha
Bierbaumer
19
20
zatížení dle přijaté metody /zatížení plného nadloží (%)
Srovnání hodnot zatížení výztuže díla ve stropě stanoveného různými metodami (v procentech plné tíhy nadloží) 110 100 90 80 70 60 50 40 30 20 10 0
87,8
90,5
81,9 75,2 59,67 52,7
43,2
výška klenby dle Protodjakonova
0
5
10
15
Výška nadloží (m)
tíha celého nadloží Terzaghi Analytický (Bulyčev), tloušťka = 5cm Plaxis (d=5 cm)
20
25
30
35
minim. hloubka platnosti výpočtu dle Bulyčeva
Protodjakonov Bierbaumer Analytický (Bulyčev), tloušťka d = 20cm Plaxis (d=20 cm)
Tunelářské odpoledne březen 2014 - Praha
Chybové aspekty numerických modelů
Chyby formulační – zadání geometrie, volba konstitutivních vztahů, materiálových vlastností, okrajových podmínek, volba typu analýzy (lineární, nelineární, odvodněné, resp. neodvodněné podmínky atd.),… Výchozí motto formulace modelu: Make everything as simple as possible, but not simpler. Albert Einstein
Chyby diskretizace – vyplývají z generace sítě a volby typu prvků Chyby numerické –
integrační chyby , chyby zaokrouhlovací, chyby iterační,…
Tunelářské odpoledne březen 2014 - Praha
Formulační chyby numerických modelů – volba dimenze modelu Volba dimenze modelu: 2D x 3D, 2D model – úlohy, v nichž jsou splněny podmínky rovinné deformace (např. liniová díla (tunely apod.) nebo rovinné napjatosti (např. tenké desky) nebo stav rotační symetrie – kruhové základy, piloty apod. ( ! nejen symetrická konstrukce, ale i podloží, včetně hladiny podzemní vody … ) Rovinná deformace:
Rotační symetrie:
3D model – nejsou splněny podmínky pro 2D model např. stav v blízkosti čelby a na čelbě tunelu (i když lze částečně simulovat ve 2D zadáním koeficientu vlivu čelby) Tunelářské odpoledne březen 2014 - Praha
Formulační chyby numerických modelů – volba charakteru prostředí Prostředí: homogenní x nehomogenní x kvazihomogenní izotropní x transversálně izotropní x anizotropní drénované x nedrénované kontinuální x diskontinuitní Drénované x nedrénované prostředí: Drénované – při přitěžovaní resp. odlehčování nevznikají v prostředí změny pórových tlaků (pomalé zatěžování, velmi propustné prostředí (např. štěrky), řešení dlouhodobé stability) Nedrénované - při přitěžovaní resp. odlehčování vznikají v prostředí změny pórových tlaků, (rychlé zatěžování, málo propustné prostředí (např. jíly), řešení krátkodobé stability) Tunelářské odpoledne březen 2014 - Praha
Formulační chyby numerických modelů – volba charakteru prostředí
Kontinuální x diskontinuální model prostředí (např. skalní prostředí s blokovou strukturou, sypký granulární materiál, vliv diskontinuit je primární )
volba numerické metody
numerická metoda modelování kontinua (metoda konečných prvků, metoda hraničních prvků,metoda konečných diferencí) –předpoklad kontinuálního přetváření tělesa, bez možnosti modelování odtržení
numerická metoda modelování diskontinua (metoda oddělených elementůsoftwary UDEC, PFC) – umožňuje modelovat separaci bloků, otevírání resp. uzavírání trhlin …
Tunelářské odpoledne březen 2014 - Praha
Formulační chyby modelů - volba konstitutivního modelu Základní aspekty konstitutivních modelů Základní aspekty chování zemin
Lineárně Mohrpružný Coulomb model
Cam-Clay model
Hypoplastický model(Mašín)
Mezní plocha stavů
N
A
A
A (pouze plocha porušení)
Závislost chování na středním napětí a pórovitosti
N
N
A
A
Nelineární chování zemin
N
N
N
A
Závislost tuhosti na historii zatěžování
N
N
N
A
Závislost tuhosti na úrovni napětí
N
N
A (pouze
A
objem.tuhost)
Tunelářské odpoledne březen 2014 - Praha
Formulační chyby modelů - volba konstitutivního modelu Základní charakteristika nejčastěji používaného konstitutivního modelu Mohr-Coulomb Mohr-Coulomb (pružný-ideálně plastický model) nejčastěji využívaný v geotech. praxi, i když nemusí poskytovat vždy zcela objektivní výsledky pružný-ideálně plastický model 5 základních charakteristik (modul pružnosti, Poissonovo číslo, soudržnost,úhel vnitřního tření, úhel dilatance), nezohledňuje stavový charakter charakteristik– využívá standardní laboratorní výstupy nezohledňuje změnu tuhosti v závislosti na přetvoření, stejný modul pružnosti při zatěžování i odlehčování - mnohdy nad tunelem dochází v modelu při vyražení díla (odlehčení) ke zvedání povrchu přehodnocena nedrénovaná smyková pevnost reálnější výsledky při řešení stabilitních úloh, deformace není schopen modelovat zcela objektivně Tunelářské odpoledne březen 2014 - Praha
Formulační chyby modelů - volba konstitutivního modelu Nicoll Highway Collapse (Singapore, 2004) Jednou z hlavních příčin bylo selhání pažícího systému dimenzovaného na základě nesprávného nastavení Mohr-Coulombova materiálového modelu v programu Plaxis.
Metoda A: nedrénovaná analýza v efektivních napětích, efektivní smykové parametry c, ϕ, přeceněná nedrénovaná smyková pevnost cu Metoda B: nedrénovaná analýza v efektivních napětích, Tunelářské odpoledne březen 2014 - Praha c=cu , ϕ=0
Formulační chyby modelů - volba konstitutivního modelu
Optimalizace: dostupnost vstupních charakteristik x výstižnost chování zeminového prostředí
www.soilmodels.info – otevřená databáze konstitutivních modelů pro numerickou analýzu , interface pro implementaci do softwarů
Tunelářské odpoledne březen 2014 - Praha
Formulační chyby modelů – zadání materiálových charakteristik Chyby: subjektivní (lidský faktor) i objektivní (komplikovaný a proměnlivý materiál) Nejistoty lze snížit především: kvalitním a dostatečným průzkumem, poskytujícím dostatečný počet výsledků laboratorních i polních zkoušek pro stanovení spolehlivých mater. charakteristik zvyšováním odborností pracovníků provádějících průzkum, lab. i polní zkoušky aplikací stochastických metod modelování, metod inverzní analýzy
Tunelářské odpoledne březen 2014 - Praha
Formulační chyby modelů – volba okrajové podmínky Standardní okrajové podmínky statické rovnováhy: musí zabránit rotačnímu i translačnímu posunu celého modelu ( v geotechnických úlohách nejčastěji podmínka tzv. tuhé vany) Okrajové podmínky konsolidační – definují v modelu propustnost či nepropustnost dané hranice vzhledem ke konsolidačním procesům – jednostranná či dvoustranná konsolidace (např. při modelování procesu konsolidace pod násypy budovanými na zvodnělém měkkém jílovitém podloží) – volba determinuje časový průběh sedání a vývoj pórových tlaků v podloží Okrajové podmínky omezující proudění vody Nezadání nebo chybné zadání okrajových podmínek vede k problémům s řešitelností výsledné soustavy rovnic, matice tuhosti není regulární a není zajištěna řešitelnost výsledné soustavy rovnic. Tunelářské odpoledne březen 2014 - Praha
Formulační chyby modelů – volba rozsahu modelu Rozsah modelu by měl být takový, aby okrajové podmínky zadávané na hranicích, neovlivňovaly výpočet v zájmové oblasti, tj. deformační hranice by měly být v místech, ve kterých se již nepředpokládají deformační změny. Testovací úloha vlivu velikosti modelu na výsledky řešení (Plaxis 2D 2010): nevyztužené dílo kruhového příčného průřezu o poloměru r= 5 m výška nadloží: h= 5 m objemová tíha okolní horniny: γ= 20 kN/m3 (homogenní prostředí) modul pružnosti okolního prostředí: E=20 MPa materiálový model: lineárně pružný variantní rozměry modelu: vzdálenost bočních svislých hranic a spodní hranice od středu díla vždy v k-násobcích poloměru díla (k=4,6,8,10,12)
Tunelářské odpoledne březen 2014 - Praha
Formulační chyby modelů – volba rozsahu modelu 60 m=12*r
Schéma parametrické modelové studie (síť identická):
50 m=10*r 30 m=6*r 20 m=4*r
Tunelářské odpoledne březen 2014 - Praha
40 m=8*r
Formulační chyby modelů – volba rozsahu modelu Srovnání svislých posunů v závislosti na rozsahu modelu(strop, počva) 0,14 vertikální posuny (m)
0,12
0,115
0,11 0,102
0,1 0,08 0,06 0,071 0,053
0,04
0,043
0,02
0,035 0,027
0 4
5
6
7 8 9 k-násobek poloměru r
10
maximální svislý posun stropu(15-ti uzlové prvky) maximální zdvih počvy (15-ti uzlové prvky) Tunelářské odpoledne březen 2014 - Praha
11
12
Formulační chyby modelů – volba rozsahu modelu Srovnání maximálních napětí pod počvou v závislosti na rozsahu modelu
max.napětí v počvě (kPa)
520 510 508
500
504
490
495
480
483
470 460 450
457
4
5
6
7 8 9 k-násobek poloměru r
10
max. hlavní napětí pod počvou(15-ti uzlové prvky)
Tunelářské odpoledne březen 2014 - Praha
11
12
Chyby diskretizační – volba typu prvku Základní faktory určující tvar aproximační funkce posunů na prvku : typ prvků počet uzlových bodů
Tunelářské odpoledne březen 2014 - Praha
Chyby diskretizační – volba typu prvku Nejčastěji používaný prvek v rovině : trojúhelník Uzlové body: ve vrcholech trojúhelníka (nejjednodušší prvek v rovině) (3 uzlový prvek) – aproximace posunů na prvku je lineární, není příliš přesný, nevystihuje zejména lokální extrémy deformací ani napětí, ve většině komerčních softwarů se nevyužívá ve vrcholech trojúhelníka a ve středech stran (6-ti uzlový prvek) – aproximace funkce posunů na prvku je polynomem 2. řádu, dostatečná přesnost v případě deformační analýzy, pro stabilitní analýzu nepřesný ve vrcholech trojúhelníka, ve středech stran i uvnitř trojúhelníka (15-ti uzlový prvek) – aproximace funkce posunů na prvku polynomem 3. řádu, doporučuje se především v případě napěťové analýzy (stabilitní úlohy, vyhodnocení čerpání smykové pevnosti apod.) Vyšší počet uzlových bodů umožňuje zpřesnit řešení, avšak představuje zvýšení dimenze soustavy rovnic, vyšší nároky na výpočetní čas, kapacitu operační paměti i disku, … Tunelářské odpoledne březen 2014 - Praha
Chyby diskretizační – volba typu prvku
vertikální posuny (m)
Srovnání svislých posunů pro různé typy prvků(strop, počva) 0,14 0,12 0,1 0,08 0,06 0,071 0,04 0,02 0 4
0,102
0,115
0,11
0,053
5
6
7
0,043
0,035
8 9 k-násobek poloměru r
10
0,027
11
maximální svislý posun stropu(15-ti uzlové prvky)
maximální zdvih počvy (15-ti uzlové prvky)
maximální svislý posun stropu (6-ti uzlové prvky)
maximální zdvih počvy (6-ti uzlové prvky)
12
Posuny pro oba typy trojúhelníkových prvků (6-ti i 15-ti uzlové) jsou posuny identické Tunelářské odpoledne březen 2014 - Praha
Chyby diskretizační – volba typu prvku
max.napětí v počvě (kPa)
Srovnání maximálních napětí pod počvou pro různé typy prvků 540
532 523
520
509 494
500 480
508
504 495
466
483
460 440
457
4
5
6
7 8 9 k-násobek poloměru r
10
11
12
max. hlavní napětí pod počvou(15-ti uzlové prvky) max. hlavní napětí pod počvou (6-ti uzlové prvky)
Maximální napětí kolem díla (v počvě) je pro různé prvky rozdílné Tunelářské odpoledne březen 2014 - Praha
Chyby diskretizační – kvalita sítě Základní chybové faktory sítě Málo hustá síť, hustší síť je nutno zvolit v místech s vyšším gradientem změny( kolem vyraženého tunelu, v okolí paty svahu, v okolí vyhloubené jámy …) – zachycení lokálních extrémů Velké zkosení prvků (ostré úhly) Příliš velký poměr mezi největším a nejmenším rozměrem prvků (tzv. aspect ratio-AR) Příliš velké rozdíly ve velikosti sousedních prvků – optimální je postupná změna velikosti prvků (do 20 %)
Tunelářské odpoledne březen 2014 - Praha
Chyby diskretizační – kvalita sítě Správné tvary AR = cca 1
Nevhodné tvary AR vysoké, ostré úhly
Velikosti sousedních prvků Příliš velký rozdíl
Špatná kvalita sítě způsobuje nepřesné řešení, numerické problémy, výsledná soustava rovnic je tzv. špatně podmíněná – tj. malá změna ve vstupních datech znamená velkou změnu v řešení. Konvergenci úlohy může rovněž narušit kombinace různých prvků v jedné úloze, kdy při spojení mají prvky na společné hraně odlišný počet uzlů. Tunelářské odpoledne březen 2014 - Praha
Chyby numerické Chyby zaokrouhlovací – zejména při aplikaci Gaussovy eliminační metody pro řešení soustavy rovnic dochází k jejich akumulaci
Chyby integrace – chyby spojené s numerickou integrací např. pro stanovení matice tuhosti s využitím určitého počtu Gaussových integračních bodů, čím vyšší počet integračních bodů, tím vyšší přesnost Chyby iteračních metod – při nevhodné volbě počáteční aproximace, iteračního kroku, nastavení přesnosti výpočtu nemusí být splněna podmínka konvergence metody, může docházet k oscilacím iteračního procesu Chyby, k nimž dochází při řešení výsledné soustavy rovnic, často souvisí se špatnou kvalitou sítě popř. špatně zadanými okrajovými podmínkami modelu. Tunelářské odpoledne březen 2014 - Praha
Za výsledek modelování vždy odpovídá realizátor výpočtu, nikoliv aplikovaná metoda ani software !!!
Děkuji za pozornost
Tunelářské odpoledne březen 2014 - Praha