Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo PODIPLOMSKI ŠTUDIJ GRADBENIŠTVA DOKTORSKI ŠTUDIJ Kandidatka: TEJA MELINK, univ. dipl. inž. grad. METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ Doktorska disertacija štev.: 239 STOCHASTIC FINITE ELEMENT METHOD IN STRUCTURE MODELLING Doctoral thesis No.: 239 Soglasje k temi doktorske disertacije je dala Komisija za doktorski študij UL na 8. redni seji 8. julija 2010. Za mentorja je bil imenovan prof. dr. Jože Korelc. Ljubljana, 16. junij 2014 Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo Komisijo za oceno ustreznosti teme doktorske disertacije v sestavi:  izr. prof. dr. Jože Korelc,  prof. dr. Darko Beg,  prof. dr. Boris Štok, UL FS,  prof. dr. Goran Turk, je imenoval Senat Fakultete za gradbeništvo in geodezijo na 10. redni seji 21. aprila 2010. Poročevalce za oceno doktorske disertacije v sestavi:  prof. dr. Darko Beg,  prof. dr. Boris Štok, UL FS,  prof. dr. Goran Turk, je imenoval Senat Fakultete za gradbeništvo in geodezijo na 3. redni seji 25. septembra 2013. Naknadno je Senat Fakultete za gradbeništvo in geodezijo na 9. redni seji 26. marca 2014, zaradi smrti poročevalca prof. dr. Darka Bega, imenoval dodatna poročevalca:  prof. dr. Tomaža Rodiča, UL NTF,  izr. prof. dr. Matjaža Dolška. Komisijo za zagovor doktorske disertacije v sestavi:  prof. dr. Matjaž Mikoš, dekan UL FGG, predsednik,  prof. dr. Jože Korelc, mentor,  prof. dr. Boris Štok, UL FS,  prof. dr. Goran Turk,  prof. dr. Tomaž Rodič, UL NTF,  izr. prof. dr. Matjaž Dolšek, je imenoval Senat Fakultete za gradbeništvo in geodezijo na 11. seji 28. maja 2014. Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo IZJAVA O AVTORSTVU Podpisana TEJA MELINK, univ. dipl. inž. grad., izjavljam, da sem avtorica doktorske disertacije z naslovom METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ. Izjavljam, da je elektronska različica v vsem enaka tiskani različici. Izjavljam, da dovoljujem objavo elektronske različice v repozitoriju UL FGG. Ljubljana, 16. junij 2014 ……………………………….. (podpis) Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. I Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Popravki Stran z napako Vrstica z napako Namesto Naj bo II Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Ta stran je namenoma prazna. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. III Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Bibliografsko-dokumentacijska stran UDK: 624.07:519.21:519.63:(043) Avtorica: Teja Melink Mentor: prof. dr. Joºe Korelc Naslov: Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij Obseg in oprema: 191 str., 67 sl., 37 pregl., 194 en. Klju£ne besede: metoda stohasti£nih kon£nih elementov, stohasti£no polje, Karhunen- Loèvova dekompozicija, stohasti£ni kon£ni elementi, avtomatsko odvajanje, perturbacijska metoda Izvle£ek: V disertaciji je razvit splo²ni simbolni pristop k avtomatizaciji metode sto- hasti£nih kon£nih elementov v modeliranju konstrukcij. Uporabljen je hibridni simbolno- numeri£ni pristop, v katerem so kombinirani: program za simbolno ra£unanje Mathe- matica, avtomatski generator kode z vgrajeno tehniko avtomatskega odvajanja AceGen, ter okolje za kon£ne elemente AceFEM. Za opis in diskretizacijo stohasti£nega polja je izbrana Karhunen-Loèvova (K-L) dekompozicija. Za izra£un deterministi£nih £lenov K-L dekompozicije je razvita formulacija stohasti£nih kon£nih elementov, ki je primerna in uporabna v poljubnem standardnem okolju za kon£ne elemente, ki omogo£a vklju£evanje uporabni²ko deniranih kon£nih elementov. Posebna pozornost je namenjena poenosta- vljeni eksponentni kovarian£ni funkciji, ki se z namenom pocenitve numeri£nega postopka v primeru srednje koreliranih stohasti£nih polj pogosto uporablja v literaturi. V disertaciji je podan dokaz, da predlagana poenostavljena kovarian£na funkcija ni pozitivno denitna, zaradi £esar je problem numeri£no nestabilen. Predlagani sta dve alternativni modika- ciji kovarian£ne funkcije, ki pocenita numeri£no re²evanje problema in hkrati zagotavljata bolj²o numeri£no stabilnost. Za izra£un odziva stohasti£nih problemov je uporabljena perturbacijska metoda vi²jega reda. V ta namen je razvita na avtomatskem odvajanju bazirana ob£utljivostna analiza (poljubnega) vi²jega reda. Izra£unani odvodi so analiti£ni. Zaradi visoke numeri£ne u£inkovitosti postopka je moºna izvedba perturbacijske metode vi²jega reda za neprimerno vi²je ²tevilo uporabljenih £lenov Karhunen-Loèvove dekompo- zicije, kot je moºno z drugimi metodami za izra£un odvodov. Ra£unski primeri so izbrani tako, da je razvita formulacija metode stohasti£nih kon£nih elementov vericirana za vse osnovne primere, pomembne za inºenirsko prakso, in sicer za razli£ne koli£ine, ki se jih modelira s stohasti£nim poljem ter za razli£ne poru²ne mehanizme, tako za mejno stanje materiala, kot za mejno stanje stabilnosti. IV Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Ta stran je namenoma prazna. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. V Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Bibliographic and documentalistic infor- mation UDC: 624.07:519.21:519.63:(043) Author: Teja Melink Supervisor: Prof. Joºe Korelc, Ph.D. Title: Stochastic nite element method in structure modelling Notes: 191 pp., 67 g., 37 tab., 194 eq. Key words: stochastic nite element method, stochastic eld, Karhunen-Loèvove de- composition, stochastic nite element, automatic dierentiation, perturbation method Abstract: In the thesis, the automation of stochastic nite element method to deal with stochastic elds is developed. The automation is based on hybrid symbolic-numeric approach. In this approach symbolic and algebraic computational system Mathematica, AceGen for the automatic code generation and the general-purpose nite element envi- ronment AceFEM are combined. The discretization of the stochastic eld is done via Karhunen-Loève (K-L) expansion. Stochastic nite elements are developed to calculate the deterministic terms in K-L expansion. Further on, the approximation of exponential covariance function that reduces the cost of numerical calculation of K-L expansion is investigated. The thesis provides a proof that the approximation of covariance function that can often be found in the literature is not positive denite. The consequence is possible numerical instability of the problem expressed as loss of positive deniteness of covariance matrix. The thesis presents two alternative modications of the covariance function that reduce the computation cost and at the same time signicantly improve the numerical stability of the procedure. For the calculation of the stochastic system response, higher-order perturbation method is chosen, based on the sensitivity analysis. For this step, the automatic dierentiation based (ADB) formulation of higher-order sensitivity analysis is developed. The advantages of the developed formulation are calculation of the exact derivatives, high numerical eciency and abstract symbolic description of the nite element code to derive higher-order derivatives. The proposed formulation of sto- chastic nite element method is veried on all basic types of numerical examples that are important in the engineering practice, i.e. for various quantities modelled by stochastic eld and for dierent failure modes: for excessive plastic deformations, as well as loss of structural stability. VI Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Ta stran je namenoma prazna. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. VII Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Zahvala Zahvaljujem se Ministrstvu za visoko ²olstvo, znanost in tehnologijo, ki je moje razisko- valno delo nanciralo v okviru programa mladih raziskovalcev. Zahvaljujem se vam, mentor prof. dr. Joºe Korelc, da ste mi zaupali in mi omogo£ili raziskovalno delo. Hvala vam, da ste me vedno sprejeli z dobro voljo, me usmerjali in mi bili vedno pripravljeni pomagati. Hvala, da ste me poleg raziskovalnega dela vklju£ili tudi v pedago²ki proces in mi zaupali vlogo somentorice pri dveh diplomskih nalogah. šelela bi se zahvaliti ºal preminulemu prof. dr. Darkotu Begu, ki je bil v £asu mojega raziskovalnega dela predstojnik katedre in mi prav tako stal ob strani ter mi pomagal pri ²irjenju znanja tudi na strokovnem podro£ju. Zahvaljujem se vsem sodelavcem na katedri: Niku Kristani£u, Ur²i ’olinc, Blaºu Hudo- bivniku, Franciju Sinurju, Blaºu ƒermelju, Barbari Gorenc, Primoºu Moºetu, Romani Hudin, Petru Skubru in Luki Pavlov£i£u ter ostalim mladim raziskovalcem: Maji, Mojci, Nata²i, Klemnu, Zlatku, Matevºu, Robertu, Jaki, Davidu, Mihi, Mateju, Nani... za dobro delovno vzdu²je, za bodrenje, ko je bilo potrebno, in za pomo£, ko sem jo potrebovala. Posebna zahvala gre tebi mama Milena in tebi tata Janko, ker sta od malega verjela vame, me usmerjala, mi svetovala in mi pri moji poti brezmejno pomagala. Hvala tudi moji ²ir²i druºini, da ste me podpirali in mi pomagali po svojih najbolj²ih mo£eh. Hvala tebi Igor, ki si me ºe kot trener atletike v dolgih letih nau£il, da se nekaj lahko naredi dobro le, £e se dela s srcem in vztrajnostjo. Hvala za vso tvojo podporo in spodbudo tudi na moji raziskovalni poti. Na koncu hvala ²e najinim trem sinovom Jakobu, Vidu in Levu, ki ste v tem obdobju s prihodom na svet dali mojemu ºivljenju in delu dodaten smisel. VIII Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Ta stran je namenoma prazna. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. IX Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. KAZALO VSEBINE 1 Uvod 1 1.1 Predstavitev problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Predstavitev stanja na obravnavanem podro£ju . . . . . . . . . . . . . . . 2 1.2.1 Diskretizacija stohasti£nih polj . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Izra£un statistike odziva konstrukcije . . . . . . . . . . . . . . . . . 5 1.2.3 Ob£utljivostna analiza . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.4 Programska okolja za metodo kon£nih elementov, ki omogo£ajo apli- kacijo stohasti£nega pristopa . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Tema doktorske disertacije . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.4 Zasnova doktorske disertacije . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 Numeri£no re²evanje stohasti£nih problemov v mehaniki kontinuuma 15 2.1 Diskretizacija stohasti£nega polja . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.1 Diskretizacija stohasti£nega polja v Karhunen-Loèvovo vrsto . . . . 18 2.1.2 Galerkinova procedura za numeri£no re²evanje homogene Fredhol- move integralske ena£be drugega reda . . . . . . . . . . . . . . . . . 21 2.2 Diskretni sistemi v mehaniki kontinuuma . . . . . . . . . . . . . . . . . . . 26 2.2.1 Osnovne ena£be v mehaniki kontinuuma . . . . . . . . . . . . . . . 27 2.2.2 Formulacija mehanskega problema z uporabo avtomatskega odvajanja 30 2.2.3 Klasikacija mehanskih problemov in njihova formulacija v metodi kon£nih elementov . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3 Perturbacijska metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3.1 Analiti£na ob£utljivostna analiza . . . . . . . . . . . . . . . . . . . 36 2.3.2 Simbolna ob£utljivostna analiza . . . . . . . . . . . . . . . . . . . . 36 2.3.3 Metoda kon£nih diferenc . . . . . . . . . . . . . . . . . . . . . . . . 39 X Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 2.4 Metoda Monte Carlo za izra£un statistike odziva konstrukcije . . . . . . . . 41 3 Hibridni simbolno-numeri£ni pristop 43 3.1 Uporaba AceGen-a v metodi kon£nih elementov. . . . . . . . . . . . . . . . 44 3.2 Avtomatsko odvajanje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4 Avtomatizacija diskretizacije stohasti£nega polja 51 4.1 Stohasti£ni kon£ni elementi . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Izvedba diskretizacije stohasti£nega polja v okolju za kon£ne elemente Ace- FEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5 Problem stabilnosti Karhunen-Loèvove dekompozicije 63 5.1 Eksponentna kovarian£na funkcija . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Okrnjena eksponentna kovarian£na funkcija . . . . . . . . . . . . . . . . . 65 5.2.1 Dokaz, da okrnjena eksponentna kovarian£na matrika ni pozitivno denitna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.2.2 Priporo£ilo o izbiri minimalne efektivne dolºine, da je kovarian£na matrika pozitivno denitna . . . . . . . . . . . . . . . . . . . . . . . 68 5.3 Modicirana okrnjena eksponentna kovarian£na funkcija . . . . . . . . . . 69 5.4 Gladka kovarian£na funkcija . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.5 Modicirana gladka kovarian£na funkcija . . . . . . . . . . . . . . . . . . . 72 5.6 Vpliv efektivne dolºine na rezultirajo£o Karhunen-Loèvovo vrsto . . . . . . 73 6 Analiti£na ob£utljivostna analiza vi²jega reda 81 6.1 Denicija ob£utljivostnega problema . . . . . . . . . . . . . . . . . . . . . 81 6.2 Ob£utljivostna analiza prvega reda . . . . . . . . . . . . . . . . . . . . . . 82 6.2.1 Avtomatizacija ob£utljivostne analize prvega reda . . . . . . . . . . 84 6.3 Ob£utljivostna analiza drugega reda . . . . . . . . . . . . . . . . . . . . . . 86 6.3.1 Avtomatizacija ob£utljivostne analize drugega reda . . . . . . . . . 90 6.4 Ob£utljivostna analiza vi²jih redov . . . . . . . . . . . . . . . . . . . . . . 93 6.4.1 Avtomatizacija ob£utljivostne analize vi²jega reda . . . . . . . . . . 94 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XI Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 6.5 Avtomatizacija ob£utljivostne analize za £asovno neodvisen nepovezan pro- blem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.6 Elasto-plasti£ni kon£ni element z vgrajeno ob£utljivostno analizo drugega reda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.6.1 Simbolni opis kon£nega elementa za elasto-plasti£no ob£utljivostno analizo drugega reda . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.7 Izvedba ob£utljivostne analize v okolju za kon£ne elemente AceFEM . . . . 111 7 Ra£unski primeri 119 7.1 Nelinearna analiza z upo²tevanjem stohasti£ne za£etne geometrijske nepo- polnosti. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 7.1.1 Za£etna geometrijska nepopolnost . . . . . . . . . . . . . . . . . . . 119 7.1.2 Upogib sinusoidnega panela . . . . . . . . . . . . . . . . . . . . . . 119 7.1.3 Cikli£no obremenjevanje jeklene epruvete . . . . . . . . . . . . . . . 126 7.2 Nelinearna stohasti£na analiza korodiranih konstrukcij . . . . . . . . . . . 139 7.2.1 Proces korozije . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 7.2.2 Tla£ena vrtljivo podprta plo²£a . . . . . . . . . . . . . . . . . . . . 140 7.2.3 Tla£en korodiran nosilec z globalnim uklonom kot merodajnim po- ru²nim mehanizmom . . . . . . . . . . . . . . . . . . . . . . . . . . 150 7.2.4 Tla£en korodiran nosilec s spremenljivim poru²nim mehanizmom . . 158 8 Zaklju£ki 167 9 Povzetek 173 10 Summary 179 XII Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Ta stran je namenoma prazna. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XIII Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. KAZALO SLIK 1.1 Procedura metode stohasti£nih kon£nih elementov. . . . . . . . . . . . . . 13 2.1 Prve ²tiri lastne funkcije za enodimenzionalno stohasti£no polje z ni£elno pri£akovano vrednostjo in enotsko standardno deviacijo ter eksponentno kovarian£no funkcijo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2 Prvih deset lastnih vrednosti za enodimenzionalno stohasti£no polje z ni- £elno pri£akovano vrednostjo in enotsko standardno deviacijo ter eksponen- tno kovarian£no funkcijo s tremi razli£nimi korelacijskimi dolºinami. . . . . 24 2.3 Primerjava realizacij dvodimenzionalnega stohasti£nega polja z ni£elno pri- £akovano vrednostjo, enotsko standardno deviacijo ter eksponentno kova- rian£no funkcijo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4 Procedura primarne analize in ob£utljivostne analize vi²jega reda. . . . . . 33 2.5 Newton-Raphsonova procedura re²evanja, prilagojena za re²evanje razvoja vektorja re²itve v vrsto za £asovno neodvisne mehanske probleme. . . . . . 40 3.1 Simultana stohasti£na poenostavitev izrazov, vgrajena v AceGen. . . . . . 45 3.2 Hibridni simbolno-numeri£ni pristop k avtomatizaciji metode stohasti£nih kon£nih elementov. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.1 Procedura avtomatizacije diskretizacije stohasti£nega polja. . . . . . . . . . 51 4.2 Razdelitev zikalne domene D na nD poddomen. . . . . . . . . . . . . . . 52 4.3 Diskretizacija dvodimenzionalne domene D in formulacija 2ˆ4-vozli²£nih stohasti£nih kon£nih elementov. . . . . . . . . . . . . . . . . . . . . . . . . 53 4.4 Diskretizacija dvodimenzionalne domene D in formulacija 2ˆ4-vozli²£nih stohasti£nih kon£nih elementov na razdalji le, z upo²tevanjem simetrije problema. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.5 Izoparametri£na preslikava iz globalnega v referen£ni koordinatni sistem za primer dvodimenzionalnega, 2ˆ4-vozli²£nega stohasti£nega kon£nega ele- menta. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 XIV Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 4.6 Skica zikalne domene, preko katere je denirano stohasti£no polje za prikaz izvedbe diskretizacije v AceFEM-u. . . . . . . . . . . . . . . . . . . . . . . 58 4.7 Prve ²tiri lastne funkcije za prikazani primer postopka diskretizacije stoha- sti£nega polja v AceFEM-u. . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.1 Eksponentna kovarian£na funkcija CA za tri razli£ne korelacijske dolºine. . 64 5.2 Eksponentna kovarian£na funkcija in njene modikacije (CA - eksponen- tna, CB - okrnjena eksponentna, CC - modicirana okrnjena eksponentna kovarian£na funkcija). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.3 Konvergenca prvih treh negativnih lastnih vrednosti K-L vrste kot funkcije ²tevila kon£nih elementov (nel) vzdolº enotske domene. . . . . . . . . . . . 67 5.4 Minimalna efektivna dolºina za eno-, dvo- in tridimenzionalno domeno, ki je potrebna, da je kovarian£na matrika ob uporabi kovarian£ne funkcije CB pozitivno denitna. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.5 Minimalna efektivna dolºina za eno-, dvo- in tridimenzionalno domeno, ki je potrebna, da je kovarian£na matrika ob uporabi kovarian£ne funkcije CC pozitivno denitna. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.6 Gladka kovarian£na funkcija (CD) in modicirana gladka kovarian£na funk- cija (CE). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.7 Minimalna efektivna dolºina za eno-, dvo- in tridimenzionalno domeno, ki je potrebna, da je kovarian£na matrika ob uporabi kovarian£ne funkcije CE pozitivno denitna. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.8 Primerjava prve lastne funkcije za razli£ne korelacijske in efektivne dolºine. 75 5.9 Primerjava £etrte lastne funkcije za razli£ne korelacijske in efektivne dolºine. 76 5.10 Najmanj²a lastna vrednost za razli£ne kovarian£ne funkcije in lc “ 0.05. . . 79 5.11 Primerjava globalne relativne napake variance stohasti£nega polja za apro- ksimirano kovarian£no funkcijo. . . . . . . . . . . . . . . . . . . . . . . . . 80 6.1 Upogib konzole. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.2 Skica deformirane konzole. . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.3 Ob£utljivost navpi£nega pomika konzole v odvisnosti od dolºine konzole in vi²ine pre£nega prereza konzole: (a) Dv{DL, (b) Dv{Dh in (c) D2v{DLDh.117 7.1 Sinusoidna panelna plo²£a. . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XV Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 7.2 Prve ²tiri lastne funkcije K-L vrste stohasti£nega polja spremembe amplitude.122 7.3 Povpre£na vrednost pomika na sredini panela ¯vm v odvisnosti od ²tevila simulacij Monte Carlo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 7.4 Standardna deviacija pomika na sredini panela σv v odvisnosti od ²tevila m simulacij Monte Carlo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 7.5 Navpi£ni pomik na sredini panela vm, izra£unan s perturbacijsko metodo razli£nih redov in simulacijami Monte Carlo. . . . . . . . . . . . . . . . . . 128 7.7 Diagram cikli£nega obremenjevanja epruvete. . . . . . . . . . . . . . . . . . 128 7.6 Skica modela cikli£no obremenjene jeklene epruvete. . . . . . . . . . . . . . 129 7.8 Izra£unane lastne vrednosti K-L dekompozicije stohasti£nega polja geome- trijske nepopolnosti. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.9 Ob£utljivosti drugega reda D2¯peq , primarna analiza izvedena s paralelizacijo Dξ2 2 na procesorju s 4 jedri in 8 nitmi. . . . . . . . . . . . . . . . . . . . . . . . 132 7.10 Vodoravni pomik (levo) in detajl diagrama vodoravnega pomika v to£ki A (desno) v odvisnosti od slu£ajnega parametra ξ1 v primeru adaptivno izbranih £asovnih korakov. . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.11 Diagram pomika uAptq (a) ter ob£utljivosti: (b) DuA{Dξ1, (c) D2uA{Dξ21 in (d) D2uA{Dξ2Dξ3 v odvisnosti od £asa. . . . . . . . . . . . . . . . . . . 135 7.12 Akumulirane plasti£ne deformacije ¯p v to£ki T in ob£utljivosti ¯p pr- eq eq vega in drugega reda na slu£ajne parametre (a) ¯peqptq, (b) D¯peq{Dξ3, (c) D2¯ p , (d) D¯p eq {Dξ2 3 eq {Dξ2Dξ3 ( ADB formulirana ob£utljivostna analiza, ˝ metoda kon£nih diferenc). . . . . . . . . . . . . . . . . . . . . . . . . . . 136 7.13 Primerjava rezultatov perturbacijske metode drugega reda in 100 simulacij Monte Carlo za slu£ajna parametra ξ1 in ξ25. . . . . . . . . . . . . . . . . . 139 7.14 Skica plo²£e. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.15 Model geometrijske nepopolnosti plo²£e. . . . . . . . . . . . . . . . . . . . 142 7.16 Sledenje obteºne poti preko vodenega pomika. . . . . . . . . . . . . . . . . 144 7.17 Diagram obteºne poti za primer, ko je debelina stanj²ana enakomerno za pri£akovano vrednost. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 7.18 Obteºni faktor v odvisnosti od pomika in slu£ajne spremenljivke ξ2. . . . . 145 7.19 Diagram mejnega obteºnega faktorja v odvisnosti od debeline plo²£e. . . . 147 7.20 Prva lastna funkcija za stohasti£no polje s kovarian£no funkcijo CA. . . . . 148 XVI Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 7.21 ƒetrta lastna funkcija za stohasti£no polje s kovarian£no funkcijo CA. . . . 148 7.22 Naklju£no izbrana realizacija stohasti£nega polja korozije plo²£e. . . . . . . 149 7.23 Deformirana plo²£a v trenutku, ko je doseºen mejni obteºni faktor (7-kratna pove£ava deformacij). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 7.24 Prvih ²est uklonskih oblik nosilca. . . . . . . . . . . . . . . . . . . . . . . . 153 7.25 Model nepopolnosti nosilca. . . . . . . . . . . . . . . . . . . . . . . . . . . 154 7.26 Prva lastna funkcija za stohasti£no polje s kovarian£no funkcijo CA. . . . . 155 7.27 ƒetrta lastna funkcija za stohasti£no polje s kovarian£no funkcijo CA. . . . 155 7.28 Diagram obteºne poti za primer, ko je debelina nosilca enakomerno stanj- ²ana za pri£akovano vrednost. . . . . . . . . . . . . . . . . . . . . . . . . . 156 7.29 Naklju£no izbrana realizacija spremembe debeline nosilca. . . . . . . . . . . 156 7.30 Deformiran nosilec v trenutku denirane γ “ γmax (10-kratna pove£ava deformacij). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 7.31 Histograma kriti£ne uklonske sile za kovarian£ni funkciji CA in CC. . . . . 159 7.32 Mejni obteºni faktor γmaxpξ1q in γmaxpξ5q, izra£unan s perturbacijsko me- todo drugega reda in 100 Monte-Carlo simulacijami. . . . . . . . . . . . . . 159 7.33 Prvih ²est uklonskih oblik nosilca. . . . . . . . . . . . . . . . . . . . . . . . 161 7.34 Model nepopolnosti nosilca. . . . . . . . . . . . . . . . . . . . . . . . . . . 162 7.35 Diagram obteºne poti za primer, ko je debelina nosilca enakomerno stanj- ²ana za pri£akovano vrednost. . . . . . . . . . . . . . . . . . . . . . . . . . 163 7.36 Primer, ko zaradi stohasti£ne porazdelitve korozije nosilec odpove lokalno. 164 7.37 Primer, ko zaradi stohasti£ne porazdelitve korozije nosilec odpove globalno. 164 7.38 Mejni obteºni faktor v odvisnosti od slu£ajnega parametra ξ3, izra£unan z dejanskimi simulacijami. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XVII Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. LIST OF FIGURES 1.1 Stochastic nite element procedure. . . . . . . . . . . . . . . . . . . . . . . 13 2.1 First four eigenfunctions for one-dimensional stochastic eld with zero mean, unit variance and exponential covariance function. . . . . . . . . . . 23 2.2 First ten eigenvalues for one-dimensional stochastic eld with zero mean, unit variance and exponential covariance function for three dierent corre- lation lengths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3 Comparison of realizations of stochastic eld with zero mean, unit variance and exponential covariance function. . . . . . . . . . . . . . . . . . . . . . 25 2.4 Procedure of primal analysis and higher-order sensitivity analysis. . . . . . 33 2.5 NewtonRaphson solution procedure adjusted for power series solution to the system of nonlinear equations. . . . . . . . . . . . . . . . . . . . . . . . 40 3.1 Simultaneous optimization of the expressions that is built in AceGen. . . . 45 3.2 Hybrid symbolic-numeric approach to automation of stochastic nite ele- ment method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.1 Procedure of automation of stochastic eld discretization. . . . . . . . . . . 51 4.2 Discretization of physical domain D into nD subdomains. . . . . . . . . . . 52 4.3 Discretization of twodimensional domain D and formulation of 2ˆ4-noded stochastic nite elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.4 Discretization of two-dimensional domain D and formulation of 2ˆ4-noded stochastic nite elements on le region with consideration of problem sim- metry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.5 Isoparametric mapping from global to reference coordinate system for two- dimensional 2ˆ4-noded stochastic nite element. . . . . . . . . . . . . . . . 54 4.6 Sketch of physical domain on which twodimensional stochastic eld is de- ned for demonstration of AceFEM procedure. . . . . . . . . . . . . . . . . 58 XVIII Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 4.7 First four eigenfunctions of numerical example for demonstration of sto- chastic eld discretization procedure in AceFEM. . . . . . . . . . . . . . . 61 5.1 Exponential covariance function CA for three dierent correlation lengths. . 64 5.2 Exponential covariance function and its modications (CA - exponential, CB - truncated exponential, CC - modication of truncated exponential). . 65 5.3 Convergence of rst three negative eigenvalues of K-L decomposition as a function of the number of nite elements (nel) along unit domain. . . . . . 67 5.4 Minimum eective length to obtain positive eigenvalues in the case of trun- cated exponential covariance function for 1D, 2D and 3D domain. . . . . . 69 5.5 Minimum eective length to obtain positive eigenvalues in the case of co- variance function CC for 1D, 2D and 3D domain. . . . . . . . . . . . . . . 71 5.6 Smooth (CD) and modied smooth (CE) covariance function. . . . . . . . . 72 5.7 Minimum eective length to obtain positive eigenvalues in the case of mo- died smooth covariance function for 1D, 2D and 3D domain. . . . . . . . 74 5.8 Comparison of rst eigenfunction for dierent correlation and eective leng- ths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.9 Comparison of fourth eigenfunction for dierent correlation and eective lengths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.10 Minimum eigenvalue for dierent covariance functions and lc “ 0.05. . . . . 79 5.11 Comparison of global relative error in case of modied covariance function. 80 6.1 Bending of cantilever beam. . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.2 Deformed cantilever beam. . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.3 Sensitivity of vertical displacement of cantilever beam in dependence of beam length and cross section height: (a) Dv{DL, (b) Dv{Dh in (c) D2v{DLDh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.1 Sinusoidal double skin cladding. . . . . . . . . . . . . . . . . . . . . . . . . 121 7.2 First four eigenfunctions of K-L decomposition of stochastic eld of ampli- tude change. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 7.3 Mean deection in the middle of cladding ¯vm in dependence on the number of Monte Carlo simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XIX Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 7.4 Standard deviation of the deection in the middle of cladding σv in de- m pendence on the number of Monte Carlo simulations. . . . . . . . . . . . . 127 7.5 Deection in the middle of cladding, obtained with perturbation method of dierent orders and Monte Carlo simulations. . . . . . . . . . . . . . . . 128 7.7 Diagram of cycling loading. . . . . . . . . . . . . . . . . . . . . . . . . . . 128 7.6 Numerical example of cyclic plasticity of circle bar. . . . . . . . . . . . . . 129 7.8 Calculated eigenvalues of K-L decomposition of stochastic eld of the geo- metric imperfections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.9 Second-order sensitivities D2¯peq , primal analysis performed with parallel Dξ2 2 computing on processor with 4 cores and 8 threads. . . . . . . . . . . . . . 132 7.10 Horizontal displacement (left) and a detail of diagram of horizontal displa- cement (right) in point A in dependence on sensitivity parameter ξ1 in case of adaptive time-stepping. . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.11 Diagram of displacement uAptq (a) and diagrams of sensitivities: (b) DuA{Dξ1, (c) D2uA{Dξ2 and (d) D2u 1 A{Dξ2Dξ3 in dependence of time. . . . . . . . . 135 7.12 Accumulated plastic deformations and rst- and second-order sensitivities of accumulated plastic deformations in point T with respect to random parameters: (a) ¯p , (d) D¯p eq ptq, (b) D¯ peq{Dξ3, (c) D2¯peq{Dξ23 eq {Dξ2Dξ3 ( ADB formulated sensitivity analysis, ˝ nite dierence method). . . . . . . 136 7.13 Comparison of the results obtained with second-order perturbation method and 100 Monte Carlo simulations for parameters ξ1 and ξ25. . . . . . . . . 139 7.14 Sketch of the plate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.15 Geometrical imperfections of the plate. . . . . . . . . . . . . . . . . . . . . 142 7.16 A load-displacement curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 7.17 A load-displacement curve for the case of expected value of plate thickness change. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 7.18 Load factor in dependence on displacement and random variable ξ2. . . . . 145 7.19 Limit load factor in dependence of plate thickness. . . . . . . . . . . . . . . 147 7.20 First eigenfunction of stochastic eld with covariance function CA. . . . . . 148 7.21 Fourth eigenfunction of stochastic eld with covariance function CA. . . . . 148 7.22 A random sample of stochastic eld of corrosion. . . . . . . . . . . . . . . . 149 XX Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 7.23 Deformed plate in the moment of limit load factor (deformations are 7- times magnied). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 7.24 First six buckling modes of the beam. . . . . . . . . . . . . . . . . . . . . . 153 7.25 Imperfections of the beam. . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 7.26 First eigenfunction of stochastic eld with covariance function CA. . . . . . 155 7.27 Fourth eigenfunction of stochastic eld with covariance function CA. . . . . 155 7.28 Load-displacement curve for the beam with uniformly reduced thickness. . 156 7.29 A random sample of beam thickness change. . . . . . . . . . . . . . . . . . 156 7.30 Deformed beam in the moment of γ “ γmax (deformations are 10-times magnied). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 7.31 Histograms of critical buckling load for covariance functions CA and CC. . 159 7.32 Limit load factor γmaxpξ1q and γmaxpξ5q, obtained with second-order per- turbation method and 100 Monte Carlo simulations. . . . . . . . . . . . . . 159 7.33 First six buckling modes of the beam. . . . . . . . . . . . . . . . . . . . . . 161 7.34 Imperfections of the beam. . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 7.35 Load-displacement curve for the beam with uniformly reduced thickness. . 163 7.36 Local buckling of the beam due to stochastically distributed corrosion. . . 164 7.37 Global buckling of the beam due to stochastically distributed corrosion. . . 164 7.38 Limit load factor in dependence on parameter ξ3, obtained with simulations.165 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XXI Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. KAZALO PREGLEDNIC 3.1 Tipi£ne izjeme avtomatskega odvajanja. . . . . . . . . . . . . . . . . . . . 50 4.1 Parametri stohasti£nega polja. . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.1 Ocena parametrov α in β za ovojnico le,min in kovarian£no funkcijo CB. . . 69 5.2 Ocena parametrov α in β za ovojnico le,min in kovarian£no funkcijo CC. . . 70 5.3 Ocena parametrov α in β za ovojnico le,min in kovarian£no funkcijo CE. . . 73 5.4 Relativna napaka prve lastne funkcije v temenu za razli£ne korelacijske in efektivne dolºine. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.5 Relativna napaka prve lastne funkcije na robu domene za razli£ne korela- cijske in efektivne dolºine. . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.6 Povpre£na relativna napaka prve lastne funkcije preko celotne domene za razli£ne korelacijske in efektivne dolºine. . . . . . . . . . . . . . . . . . . . 78 5.7 Povpre£na relativna napaka prvega £lena K-L dekompozicije za razli£ne korelacijske in efektivne dolºine. . . . . . . . . . . . . . . . . . . . . . . . . 78 6.1 Izjeme, ki jih je potrebno podati h klicu avtomatskega odvajanja po para- metru ξs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 1 6.2 Izjeme, ki jih je potrebno podati h klicu avtomatskega odvajanja po para- metru ξs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 2 6.3 Izjeme, ki jih je potrebno podati h klicu avtomatskega odvajanja po para- metru ξs , k ă n. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 k 6.4 Izjeme, ki jih je potrebno podati h klicu avtomatskega odvajanja po para- metru ξs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 n 6.5 Izjeme, ki jih je v primeru nepovezanega in £asovno-neodvisnega problema potrebno podati h klicu avtomatskega odvajanja po parametru ξs . . . . . 100 1 6.6 Izjeme, ki jih je v primeru nepovezanega in £asovno-neodvisnega problema potrebno podati h klicu avtomatskega odvajanja po parametru ξs . . . . . 100 2 XXII Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 6.7 Izjeme, ki jih je v primeru nepovezanega in £asovno-neodvisnega problema potrebno podati h klicu avtomatskega odvajanja po parametru ξs . . . . . 101 n 6.8 Materialni in geometrijski parametri. . . . . . . . . . . . . . . . . . . . . . 112 6.9 Primerjava rezultatov metode kon£nih elementov in Bernoullijeve grede. . . 116 7.1 Prve ²tiri lastne vrednosti K-L vrste stohasti£nega polja spremembe am- plitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 7.2 Ra£unski CPU £as, potreben za izra£un razvoja odziva konstrukcije v vrsto. 124 7.3 Pri£akovana vrednost (¯vm) in standardna deviacija (σv ) navpi£nega po- m mika vm v cm ter skupni CPU £as, potreben za izra£un odziva. . . . . . . . 125 7.4 Materialni parametri epruvete. . . . . . . . . . . . . . . . . . . . . . . . . . 127 7.5 Parametri stohasti£nega polja geometrijske nepopolnosti. . . . . . . . . . . 130 7.6 Primerjava absolutnega ra£unskega £asa za perturbacijsko metodo drugega reda z ADB formulirano ob£utljivostno analizo in metodo kon£nih diferenc (elasto-plasti£en material). . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.7 Primerjava absolutnega ra£unskega £asa za perturbacijsko metodo drugega reda z ADB formulirano ob£utljivostno analizo in metodo kon£nih diferenc hiperelasti£en material). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 7.8 Materialni parametri plo²£e. . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.9 Dimenzije plo²£e. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.10 Parametri stohasti£nega polja korozije. . . . . . . . . . . . . . . . . . . . . 146 7.11 Pri£akovana vrednost (¯γmax) in standardna deviacija (σγ ) mejnega ob- max teºnega faktorja γmax, ²tevilo odvodov, ki jih je potrebno izra£unati pri perturbacijski metodi 2. reda (nodv) in skupni ra£unski CPU £as, potreben za posamezno metodo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 7.12 Materialni parametri nosilca. . . . . . . . . . . . . . . . . . . . . . . . . . . 152 7.13 Dimenzije nosilca. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 7.14 Parametri stohasti£nega polja korozije. . . . . . . . . . . . . . . . . . . . . 152 7.15 Potrebno ²tevilo £lenov K-L vrste za 80 % natan£nost, £as centralne pro- cesne enote (CPU) za sestavo kovarian£ne matrike in gostota kovarian£ne matrike za razli£ne kovarian£ne funkcije. . . . . . . . . . . . . . . . . . . . 154 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XXIII Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 7.16 Pri£akovana vrednost (¯γmax) in standardna deviacija (σγ ) mejnega ob- max teºnega faktorja γmax, ²tevilo odvodov, ki jih je potrebno izra£unati pri perturbacijski metodi (nodv) in skupni CPU £as za razli£ne metode. . . . . 158 7.17 Parametri stohasti£nega polja korozije. . . . . . . . . . . . . . . . . . . . . 160 7.18 Potrebno ²tevilo £lenov K-L vrste za 80 % natan£nost variance, £as cen- tralne procesne enote (CPU) za sestavo kovarian£ne matrike in gostota kovarian£ne matrike za razli£ne kovarian£ne funkcije. . . . . . . . . . . . . 160 7.19 Pri£akovana vrednost (¯γmax) in standardna deviacija (σγ ) mejnega ob- max teºnega faktorja γmax, ²tevilo odvodov, ki jih je potrebno izra£unati pri perturbacijski metodi (nodv) in skupni CPU £as za razli£ne metode. . . . . 163 XXIV Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Ta stran je namenoma prazna. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XXV Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. LIST OF TABLES 3.1 Typical automatic dierentiation exceptions. . . . . . . . . . . . . . . . . . 50 4.1 Parameters of stochastic eld. . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.1 Parameters α and β that dene envelopes for estimation of le,min for co- variance function CB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.2 Parameters α and β that dene envelopes for estimation of le,min for co- variance function CC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3 Parameters α and β that dene envelopes for estimation of le,min for co- variance function CE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.4 Relative error of rst eigenfunction in its turning point for dierent corre- lation and eective lengths. . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.5 Relative error of rst eigenfunction in boundary point of the domain for dierent correlation and eective lengths. . . . . . . . . . . . . . . . . . . . 77 5.6 Average error of rst eigenfunction, calculated over the whole domain for dierent correlation and eective lengths. . . . . . . . . . . . . . . . . . . . 78 5.7 Average relative error of the rst term in K-L decomposition for dierent correlation and eective lengths. . . . . . . . . . . . . . . . . . . . . . . . . 78 6.1 Exceptions that need to be dened for automatic dierentiation with re- spect to parameter ξs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 1 6.2 Exceptions that need to be dened for automatic dierentiation with re- spect to parameter ξs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 2 6.3 Exceptions that need to be dened for automatic dierentiation with re- spect to parameter ξs , k ă n. . . . . . . . . . . . . . . . . . . . . . . . . . 96 k 6.4 Exceptions that need to be dened for automatic dierentiation with re- spect to parameter ξs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 n 6.5 Exceptions that need to be dened for automatic dierentiation with re- spect to parameter ξs in case of uncoupled and time-independent problem. 100 1 XXVI Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 6.6 Exceptions that need to be dened for automatic dierentiation with re- spect to parameter ξs in case of uncoupled and time-independent problem. 100 2 6.7 Exceptions that need to be dened for automatic dierentiation with re- spect to parameter ξs in case of uncoupled and time-independent problem. 101 n 6.8 Material parameters and geometry. . . . . . . . . . . . . . . . . . . . . . . 112 6.9 Comparison of the results of nite element method and Bernoulli beam. . . 116 7.1 First four eigenvalues of stochastic eld of amplitude change. . . . . . . . . 122 7.2 CPU time needed for calculation of the response. . . . . . . . . . . . . . . 124 7.3 Expected value (¯vm) and standard deviation (σv ) of the vertical displa- m cement vm in cm and total CPU time, needed to calculate the response of the cladding. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7.4 Material parameters of the circle bar. . . . . . . . . . . . . . . . . . . . . . 127 7.5 Parameters for stochastic eld of geometric imperfection. . . . . . . . . . . 130 7.6 Comparison of absolute calculation time needed to perform analysis with second-order perturbation method with ADB formulated sensitivity ana- lysis and nite dierence method (elasto-plastic material). . . . . . . . . . 137 7.7 Comparison of absolute calculation time needed to perform analysis with second-order perturbation method with ADB formulated sensitivity ana- lysis and nite dierence method (hyperelastic material). . . . . . . . . . . 138 7.8 Material parameters of plate. . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.9 Plate dimensions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.10 Parameters for the stochastic eld of corrosion. . . . . . . . . . . . . . . . 146 7.11 Expected value (¯γmax) and standard deviation (σγ ) of the ultimate load max factor γmax, number of derivatives that have to be calculated in case of second-order perturbation method (nodv) and total CPU time, needed to perform the analysis with dierent methods. . . . . . . . . . . . . . . . . . 151 7.12 Material parameters of beam. . . . . . . . . . . . . . . . . . . . . . . . . . 152 7.13 Beam dimensions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 7.14 Parameters of the stochastic eld of corrosion. . . . . . . . . . . . . . . . . 152 7.15 Number of retained terms in K-L decomposition to attain 80 % accuracy, central processor unit (CPU) time to assemble covariance matrix and co- variance matrix density for dierent covariance functions. . . . . . . . . . . 154 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. XXVII Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 7.16 Expected value (¯γmax) and standard deviation (σγ ) of the limit load max factor γmax, number of derivatives that have to be calculated in case of perturbation method (nodv) and total CPU time, needed to perform the analysis with dierent methods. . . . . . . . . . . . . . . . . . . . . . . . . 158 7.17 Parameters for the stochastic eld of corrosion. . . . . . . . . . . . . . . . 160 7.18 Number of retained terms in K-L decomposition to attain 80 % accuracy of stochastic eld variance, central processor unit (CPU) time to assemble covariance matrix and covariance matrix density for dierent covariance functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 7.19 Expected value (¯γmax) and standard deviation (σγ ) of the ultimate load max factor γmax, number of derivatives that have to be calculated in case of perturbation method (nodv) and total CPU time, needed to perform the analysis with dierent methods. . . . . . . . . . . . . . . . . . . . . . . . . 163 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 1 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 1 Uvod 1.1 Predstavitev problema V inºenirskih problemih imamo opravka z naklju£nostjo pri skoraj vseh parametrih pro- blema. Naklju£no se spreminjajo lastnosti materiala, obteºba, geometrijske nepopolnosti, proces korozije, poroznost in druge spremenljivke. To pomeni, da konkretna vrednost posameznih parametrov ni dolo£ljiva vnaprej, temve£ je bolj ali manj naklju£na in se obi£ajno naklju£no spreminja tudi preko obmo£ja problema. Skozi zgodovino £love²tva se je razvilo ve£ razli£nih pristopov, kako re²evati inºenirske probleme. Najenostavnej²a metoda, ki se je uporabljala v dolgem obdobju v zgodovini in se na manj razvitih obmo£jih sveta uporablja ²e danes, je metoda posku²anja in napak. Na ta na£in so se konstrukcije spreminjale in izbolj²evale predvsem na podlagi preteklih izku²enj. Z razvojem analiti£nih matemati£nih metod in kasneje z razvojem ra£unalnikov je bil narejen pomemben korak k novim pristopom in izbolj²evanju metod re²evanja inºe- nirskih problemov. Dandanes je, poleg uporabnosti in udobja, ki naj bi ga konstrukcija nudila uporabnikom, glavno vodilo pri na£rtovanju konstrukcij varnost. Zahteva se, da je verjetnost poru²itve manj²a od predpisane. Kolik²na mora biti zanesljivost konstrukcije, je odvisno predvsem od namembnosti konstrukcije, od predvidenega ²tevila uporabni- kov v njej ter od velikosti ²kode, ki bi jo imela poru²itev konstrukcije na okolje (npr. jedrske elektrarne in naftne plo²£adi). V tem okviru sta se razvila dva pristopa: stoha- sti£ni in deterministi£ni. Zahtevnej²i izmed njiju, ki probleme obravnava bolj realno, je stohasti£ni pristop. Ta pristop upo²teva, da vhodni parametri problema nimajo to£no dolo£enih vrednosti, temve£ z neko verjetnostjo zavzamejo naklju£no vrednost na ome- jenem ali neomejenem obmo£ju. Nadalje obstajata v grobem dva na£ina stohasti£nega pristopa: enostavnej²i, ki privzame, da nek parameter zavzame sicer naklju£no vrednost, vendar je le-ta konstantna po celotnem obmo£ju problema. Drugi, zahtevnej²i stohasti£ni pristop, upo²teva, da je vrednost stohasti£nega parametra naklju£na in da se naklju£no spreminja tudi po domeni problema. Glavna teºava stohasti£nega pristopa je, da upo- ²tevanje stohasti£nosti pojavov bistveno pove£a numeri£no zahtevnost in kompleksnost problema. Zaradi tega se v inºenirski praksi ve£inoma uporablja enostavnej²i, determini- sti£ni pristop. V tem primeru ima vsak vhodni parameter to£no dolo£eno (determinirano) vrednost, vse moºne scenarije zaradi stohasti£ne narave pojavov pa se posku²a zajeti z upo²tevanjem ekstremov ali srednjih vrednosti stohasti£nih parametrov. Standardi, ki urejajo podro£je projektiranja gradbenih konstrukcij in k uporabi katerih 2 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. so inºenirji v Sloveniji in ve£ini preostalih evropskih drºav zakonsko zavezani, se imenu- jejo Evrokodi in dovoljujejo uporabo obeh pristopov, stohasti£nega in deterministi£nega. Vendar se pomanjkanje enostavnih in u£inkovitih metod za stohasti£ni pristop kaºe tudi tu, saj predpisani standardi temeljijo predvsem na deterministi£nem pristopu, stohasti£ni pristop pa se uporablja predvsem pri nadaljnjem razvoju Evrokodov, kar je v delu Evro- kodov, ki obravnava osnove projektiranja konstrukcij (SIST EN 1990), tudi eksplicitno zapisano. Primerjave rezultatov, dobljenih s stohasti£nim in deterministi£nim pristopom, pokaºejo, da deterministi£en pristop kljub uporabi ekstremnih vrednosti parametrov, ki so na videz na varni strani, ne zajamejo nujno najslab²ega moºnega scenarija, do katerega lahko pride zaradi naklju£ne kombinacije posameznih parametrov. Poleg tega deterministi£en pristop ne vodi do optimalnih dimenzij inºenirskih konstrukcij. Ob upo²tevanju dejstva, da se zaradi konkuren£nega boja na svetovnem trgu i²£e vedno nove izbolj²ave in moºnosti ni- ºanja stro²kov, je dolgoro£no korak v smeri stohasti£nega modeliranja neizogiben. Prav zaradi tega je stohasti£ni pristop (kot najbolj realen pribliºek opisa dejanskega problema) ºe od nekdaj deleºen velike pozornosti, kljub temu pa je le redko uporabljen v praksi. Kot re£eno, so glavni razlogi, da je stohasti£ni pristop tako redko uporabljen v praksi, nepoznavanje stohasti£nih lastnosti inºenirskih problemov, kompleksnost stohasti£nega pristopa in dejstvo, da uporaba parametrov z verjetnostno porazdelitvijo enormno pove£a velikost in numeri£no zahtevnost problema. Razvoj novih metod je tako pogojen z zmo- gljivostjo ra£unalnikov. V zadnjih dveh desetletjih je, vzporedno z rastjo in napredkom v ra£unalni²ki tehnologiji, raslo tudi zanimanje za stohasti£en pristop in posledi£no je raslo ²tevilo novo razvitih metod. Vendar potreba po enovitem pristopu, ki bi omogo£al sin- tezo modeliranja stohasti£nih pojavov ter mehanske analize inºenirskih problemov, ostaja. Prav tako ²e vedno ostaja potreba po izbolj²avah in razvoju novih metod, s katerimi bi bilo modeliranje stohasti£nosti enostavnej²e in u£inkovitej²e. 1.2 Predstavitev stanja na obravnavanem podro£ju Vrednost spremenljivk, ki se stohasti£no spreminjajo po prostoru in/ali £asu, se lahko opi²e s stohasti£nim poljem (prostor) ali procesom (£as). V doktorski disertaciji se omejimo na probleme, kjer se obravnavane koli£ine stohasti£no spreminjajo po domeni problema, kar opi²emo s stohasti£nim poljem (Phoon et al., 2002, Ghanem in Spanos, 2003, Grigoriu, 2006). Kot je bilo ºe v uvodu omenjeno, se kljub velikemu interesu raziskave na tem podro£ju intenzivneje opravljajo ²ele v zadnjih dveh desetletjih, saj upo²tevanje stoha- sti£nosti bistveno pove£a numeri£no zahtevnost problema. Tendenca ve£ine v zadnjem £asu razvitih metod je v smeri moºnosti uporabe metode kon£nih elementov, ki se na ta na£in raz²iri v metodo stohasti£nih kon£nih elementov. Metoda stohasti£nih kon£nih ele- mentov je v splo²nem sestavljena iz dveh korakov. Prvi korak predstavlja reprezentacijo Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 3 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. in diskretizacijo stohasti£nega polja, drugi korak pa izra£un statistike odziva konstrukcije. Mreºa kon£nih elementov je za posamezen korak lahko razli£na. V naslednjih treh pod- poglavjih je predstavljeno obstoje£e stanje in do danes razvite metode za vsakega izmed obeh korakov ter posebej za perturbacijsko metodo, ki je ena izmed metod za izra£un statistike odziva konstrukcije. Tem trem podpoglavjem sledi ²e podpoglavje, v katerem so predstavljena pomembnej²a do danes razvita programska okolja za metodo kon£nih elementov, ki omogo£ajo aplikacijo stohasti£nega pristopa. 1.2.1 Diskretizacija stohasti£nih polj Stohasti£no polje je za implementacijo v metodo stohasti£nih kon£nih elementov potrebno ustrezno diskretizirati. Diskretizacija stohasti£nega polja pomeni, da se zvezno stoha- sti£no polje aproksimira s kon£nim ²tevilom stohasti£nih spremenljivk. Metode diskreti- zacije na splo²no razdelimo v dve skupini. V prvo skupino lahko razvrstimo metode, po katerih se stohasti£no polje reprezentira tako, da se celotno domeno, na kateri je problem deniran, razdeli na ve£ poddomen ozi- roma kon£ne elemente. Na vsakem kon£nem elementu je stohasti£no polje diskretizirano z eno ali ve£ slu£ajnimi spremenljivkami. ’tevilo slu£ajnih spremenljivk je tako proporci- onalno ²tevilu kon£nih elementov. Koreliranost vrednosti stohasti£nega polja preko roba elementa dolo£a kovarian£na funkcija. V splo²nem je za vi²jo natan£nost potrebnih ve£ kon£nih elementov, vendar gostej²a mreºa pomeni vi²jo koreliranost kon£nih elementov in posledi£no korelacijsko matriko blizu singularnosti, kar lahko povzro£i numeri£ne ne- stabilnosti (Matthies et al., 1997, Liu, 1993, Keese 2003, Stefanou 2009, Schuëller 2006). Poleg tega z nara²£anjem ²tevila kon£nih elementov nara²£a ²tevilo slu£ajnih spremenljivk. ’tevilo slu£ajnih spremenljivk dolo£a dimenzijo stohasti£nega prostora, numeri£na zah- tevnost problema zato na splo²no raste eksponentno s stohasti£no dimenzijo. Re²ljivost je tako pogojena z velikostjo deterministi£nega problema, ºeleno natan£nostjo rezultatov in hitrostjo ra£unalnika. V to skupino lahko uvrstimo: • interpolacijsko metodo ali metodo oblikovnih funkcij (Liu et al. 1986a,b), po kateri se slu£ajno polje interpolira z oblikovnimi funkcijami kon£nih elementov, • metodo sredi²£ne to£ke (ang. midpoint method), pri kateri vrednost slu£ajne spre- menljivke posameznega kon£nega elementa predstavlja vrednost stohasti£nega polja v sredi²£ni to£ki tega elementa (Li in Kiureghian, 1993), • metodo lokalnega povpre£enja (Vanmarcke, 1983), pri kateri se stohasti£ni proces preko kon£nega elementa predstavi s povpre£no vrednostjo procesa na tem elementu ter • metodo uteºnih ostankov (ang. weighted integral method, Noh, 2005). 4 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Korak naprej v smislu re²ljivosti problema so naredile metode, ki stohasti£no polje repre- zentirajo z vrsto. Prednost teh metod je redukcija velikega ²tevila koreliranih stohasti£nih spremenljivk na manj²e ²tevilo nekoreliranih stohasti£nih spremenljivk. V to skupino me- tod spadajo: • Razvoj v Karhunen-Loèvovo vrsto ali Karhunen-Loèvova dekompozicija (Ghanem in Spanos, 2003), ki predstavlja razvoj stohasti£nega procesa v ortogonalno vrsto, v kateri so ortogonalne funkcije izbrane kot lastne funkcije homogene Fredholmove integralske ena£be druge vrste s kovarian£no funkcijo kot jedrom. • Spektralna reprezentacija (predstavljena npr. v Stefanou, 2007) stohasti£no polje razvije v trigonometri£no vrsto z amplitudami in faznim zamikom kot stohasti£nimi spremenljivkami problema. • Razvoj z linearnim optimalnim ocenjevanjem (ang. expansion optimal linear esti- mation method, Li, 1993) stohasti£ni proces izrazi kot linearno funkcijo vektorja vozli²£nih vrednosti v smislu njegove spektralne dekompozicije. U£inkovitost zgoraj na²tetih metod je pogojena z zahtevnostjo ra£unanja deterministi£- nih £lenov v vrsti ter odvisnostjo natan£nosti rezultatov od ²tevila £lenov v vrsti. Glede slednje zahteve je metoda razvoja v Karhunen-Loèvovo (v nadaljevanju K-L) vrsto op- timalna, saj je napaka aproksimacije za izbrano ²tevilo £lenov vrste minimalna v smislu srednje kvadratne napake (ang. mean square error). Zaradi tega v okviru doktorske disertacije za opis stohasti£nega polja uporabljamo K-L vrsto in si je potrebno nekoliko natan£neje pogledati stanje do danes razvitih metodologij za izra£un £lenov K-L vrste. Pri razvoju v K-L vrsto deterministi£ne dele £lenov vrste predstavljajo lastne vrednosti in lastne funkcije, ki so re²itev homogene Fredholmove integralske ena£be druge vrste, pri kateri je jedro kovarian£na funkcija. ’ibka to£ka K-L vrste je prav v teºavnosti re²evanja te ena£be, saj analiti£ne re²itve obstajajo le za dolo£ene tipe kovarian£nih funkcij ter za enostavne in pravilne geometrijske domene problema (za tak enostaven primer so izpeljane analiti£ne re²itve npr. v Ghanem in Spanos, 2003). V splo²nem je to ena£bo potrebno re- ²iti numeri£no. Razvitih je bilo ve£ metod numeri£nega re²evanja, npr. uporaba diskretne kosinusne transformacije (Courmontagne, 1999), re²evanje z uporabo razli£nih ortogonal- nih baz (npr. triangularnih, Babolian et al., 2008, val£kov (ang. wavelet), Phoon et al., 2002, metoda mnogokratnih polov, Frauenfelder et al., 2005, in Galerkinova metoda, Ghanem in Spanos, 2003). Za implementacijo problema v okoljih za metodo kon£nih elementov je od na²tetih najprimernej²a Galerkinova metoda, ki Fredholmovo integral- sko ena£bo prevede na posplo²en problem lastnih vrednosti. Slabost te metode je, da je dobljena kovarian£na matrika polna in zato zelo draga za re²evanje. Da bi pocenili re²eva- nje tega posplo²enega problema lastnih vrednosti, je bila zato s strani nekaterih avtorjev (Phoon et al., 2002, Li et al., 2008, Stefanou et al., 2005, Stefanou, 2009) uporabljena poe- Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 5 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. nostavljena kovarian£na funkcija, s katero je v primeru nizko koreliranih stohasti£nih polj kovarian£na matrika razpr²ena in je zato moºna uporaba u£inkovitej²ih algoritmov. Ven- dar imajo ti poskusi uporabe modicirane kovarian£ne funkcije tudi negativne posledice, saj le-ta lahko privede do numeri£ne nestabilnosti postopka, ki pa je v obstoje£ih £lankih ve£inoma prezrta. Vzroki in posledice numeri£ne nestabilnosti so podrobneje predstavljeni v poglavju 5. 1.2.2 Izra£un statistike odziva konstrukcije Drugi del analize v metodi stohasti£nih kon£nih elementov predstavlja izra£un statistike odziva konstrukcije. Obstoje£e metode za izra£un statistik odziva lahko razdelimo na tri tipe: • simulacije Monte Carlo. Pri tej metodi se vrednost stohasti£nih spremenljivk v vsaki simulaciji dolo£i naklju£no (z generatorji naklju£nih ²tevil, ki so vgrajeni v ve£ino programskih jezikov). S temi vrednostmi se nato re²uje deterministi£ni pro- blem. Izvede se ºeleno ²tevilo simulacij, v vsaki simulaciji se beleºi ºelena koli£ina (npr. pomik, kriti£na obteºba). Statistiko odziva se nato dolo£i s standardnimi statisti£nimi metodami. Posledica tega je, da je natan£nost re²itve odvisna od ²te- vila simulacij, kar omejuje uporabnost te metode predvsem pri problemih ve£jih dimenzij (deterministi£nih in/ali stohasti£nih) ali ko nas zanimajo repne vrednosti odziva (npr. pri analizi zanesljivosti). Kljub temu ostaja metoda Monte Carlo ena najbolj uporabljenih, saj je zelo robustna, enostavna za implementacijo, ra£unski £as pa je zaradi medsebojne neodvisnosti posameznih simulacij moºno skraj²ati tudi s paralelizacijo procesa. • Perturbacijska metoda (Sudret, 2002, Stefanou, 2009, Lin et al., 2008) vektor odziva razvije v Taylorjevo vrsto okrog pri£akovanih vrednosti in omogo£a izra£un statistik drugega reda. Za izra£un posameznih £lenov Taylorjeve vrste je potrebno izra£u- nati odvode odziva konstrukcije po slu£ajnih spremenljivkah. Analiza, s katero se odvode izra£una, se imenuje ob£utljivostna analiza. Obstaja ve£ metod za izvedbo ob£utljivostne analize, ki bodo posebej predstavljene v poglavju 1.2.3. Natan£nost perturbacijske metode je odvisna od zveznosti in nelinearnosti odziva, koecienta variacije stohasti£nega polja in reda, do katerega se Taylorjeva vrsta izra£una. Sla- bost te metode je, da je uporabna le za primere, ko je odziv v odvisnosti od spre- minjajo£ih se vhodnih parametrov zvezen. V mehanskih problemih to pomeni, da lahko metoda dovolj natan£no opi²e le primere, pri katerih se poru²ni mehanizem v odvisnosti od spreminjajo£ih se vhodnih parametrov ne spreminja. Prednost per- turbacijske metode je, da je ne glede na na£in izra£una £lenov Taylorjeve vrste praviloma numeri£no veliko u£inkovitej²a v primerjavi z Monte Carlo metodo. 6 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. • Spektralna metoda stohasti£nih kon£nih elementov odziv konstrukcije aproksimira z ortogonalnimi Hermitovimi polinomi (polinomski kaos, Ghanem in Spanos, 2003, Field in Grigoriu, 2004, Gaignaire et al., 2006). Ob zahtevi, da je napaka zaradi uporabe kon£nega ²tevila polinomov minimalna, se problem pretvori v raz²irjen sis- tem ena£b, katerega re²itev omogo£a neposreden izra£un statistike odziva. Slabost te metode je, da je uporabna (numeri£no u£inkovita) le za linearne probleme. Me- toda je tudi omejena z zmogljivostjo dana²njih ra£unalnikov, saj je raz²irjeni sistem linearnih ena£b neprimerno ve£ji v primerjavi z deterministi£nim sistemom. 1.2.3 Ob£utljivostna analiza Ob£utljivostna analiza je analiza, ki izra£una ob£utljivost (odvod) funkcije odziva glede na izbrani parameter. Glede na red iskanega odvoda je potrebno izvesti ob£utljivostno analizo tega reda. Odvod prvega reda se izra£una z ob£utljivostno analizo prvega reda, odvod drugega reda se izra£una z ob£utljivostno analizo drugega reda ali v splo²nem: odvod i-tega reda se izra£una z ob£utljivostno analizo i-tega reda. V literaturi je zelo redko zaslediti uporabo ob£utljivostne analize vi²jega reda (Bebamzadeh in Haukaas, 2008, RezaieePajand in Kadkhodaye Bahre, 2011, Ozaki, Kimura in Berz 1995). Razlog je predvsem v zahtevnem in zamudnem ro£nem izpeljevanju analiti£nih izrazov za vi²je odvode ali na drugi strani velikem ra£unskem £asu in akumuliranju numeri£nih napak ob uporabi numeri£nih metod. Do danes se je razvilo ve£ metod za izra£un ob£utljivostne analize: • Metoda kon£nih diferenc (glej npr. Haftka in Adelman, 1989) je pribliºna metoda in najenostavnej²a izmed metod za izra£un ob£utljivostne analize. Metoda je iz- peljana iz razvoja funkcije v Taylorjevo vrsto. Lahko so uporabljene backward (izbira diferenc nazaj), forward (izbira diferenc naprej) ali simetri£ne (ang. cen- tral) diference. Prednost metode kon£nih diferenc je, da je robustna in primerna za poljubne probleme. Glavna slabost te metode je, da je zaradi neodstranljivih napak nestabilna. Obstajata namre£ dva vira napak: prvi vir je napaka zaradi zanema- ritve vi²jih £lenov Taylorjeve vrste, iz katere je metoda kon£nih diferenc izpeljana. Drugi vir pa je zaokroºitvena napaka, ki nastane zaradi kon£ne natan£nosti zapisa ²tevil v ra£unalniku (t.i. formata plavajo£e vejice), kar pomeni, da je ra£unalni²ki zapis posameznega realnega ²tevila omejen le na dolo£eno ²tevilo decimalnih mest natan£no. Medtem ko se prva napaka z manj²anjem kon£ne diference manj²a, se druga napaka ve£a. Zato je potrebno biti zelo pazljiv pri izbiri ustrezne diference, pri £emer pri kompleksnej²ih primerih, med katere v tem okviru spada tudi metoda kon£nih elementov, izbira ustrezne diference ni tako enozna£na in je obi£ajno potreb- nih ve£ iteracij, da se doseºe ustrezno natan£nost rezultata. Zato je metoda kon£nih diferenc v primeru velikega ²tevila ob£utljivostnih parametrov in/ali kompleksnih Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 7 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. mehanskih problemov numeri£no zelo draga, pri posebnih, kompleksnej²ih primerih in ob£utljivostih, ki so znatno niºjega velikostnega reda, kot je odziv sam, pa dovolj natan£ne re²itve zaradi neodstranljivih napak sploh ni moºno dobiti. V ta namen je bilo narejenih nekaj izbolj²av, s katerimi je doseºena bolj²a stabilnost metode (npr. Voorhees, Millwater in Bagley 2011, Cho in Kim, 2006) ali izbolj²ana numeri£na u£inkovitost (npr. Press et al. 1988). Uporabe kon£nih diferenc za izra£un vi²jih odvodov v literaturi skorajda ni, saj je problem predvsem v dolgih ra£unskih £asih ter akumuliranju numeri£nih napak, ki so v primeru odvodov vi²jega reda velikokrat neodstranljive (Vorhees, Millwater in Bagley, 2011). • Variacijska metoda ob£utljivostne analize (ang. continuum derivatives). Po tej me- todi se osnovne ena£be mehanskega problema (to so parcialne diferencialne ena£be ali integralska ena£ba, dobljena z uporabo principa o virtualnem delu) najprej od- vaja, sledi diskretizacija ena£b za uporabo v metodi kon£nih elementov (Akbari et. al, 2010). Obi£ajno se uporabi enako diskretizacijo za mehansko in ob£utljivostno analizo, lahko pa sta tudi razli£ni. • Diskretna metoda ob£utjivostne analize (ang. discrete derivatives). Po tej metodi se osnovne ena£be problema najprej diskretizira in ²ele nato odvaja (Keulen et al., 2005, Choi in Kim, 2005, Akbari et. al, 2010). Za slednji, variacijsko in diskretno ob£utljivostno analizo, se lahko uporabi pristop z di- rektnim odvajanjem (ang. direct dierentiation) ali pridruºeni pristop (ang. adjoint approach) (Choi in Twu 1989, Michaleris et al. 1994, Keulen et al. 2005). Oba pristopa vodita do istega rezultata, razlika pa je v numeri£ni u£inkovitosti. V primeru stacio- narnih linearnih sistemov je direktno odvajanje primernej²e, ko imamo majhno ²tevilo ob£utljivostnih parametrov in veliko ²tevilo funkcij odzivov, medtem ko je drugi pristop primernej²i za primer, ko nas zanima manj²e ²tevilo funkcij odzivov in ve£je ²tevilo ob- £utljivostnih parametrov. Pri kompleksnej²ih, £asovno odvisnih povezanih problemih ter iskanih vi²jih odvodih je izra£un, po katerem pristopu je re²evanje problema numeri£no u£inkovitej²e, bolj zapleten, razlika med numeri£no u£inkovitostjo enega ali drugega pri- stopa pa je manj²a (za konkretne izra£une glej npr. Michaleris et al. 1994). Ena£be, dobljene z diskretizacijo v okviru metode kon£nih elementov, se obi£ajno impli- citno re²uje z uporabo Newton-Raphsonove iterativne metode. ƒe se isti pristop uporabi tudi pri ob£utljivostni analizi, se dobi analiti£ne vrednosti odvodov, zato sta variacijska in diskretna ob£utljivostna analiza analiti£ni metodi. Kljub vsemu pa je, ne glede na uporabljeno variacijsko ali diskretno ob£utljivostno analizo, izpeljevanje odvodov osnov- nih ena£b in ro£no vna²anje dobljenih izrazov v kodo kon£nega elementa zelo kompleksno in dolgotrajno delo, ²e posebej v primeru £asovno odvisnih povezanih problemov in pri iskanih vi²jih odvodih. Zato se v praksi pogosto uporablja poenostavitev, ko se odvode na nivoju kon£nih elementov izra£una numeri£no z metodo kon£nih diferenc, nato pa se 8 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. globalni problem re²i analiti£no. Tak pristop k variacijski ali diskretni analizi se imenuje semi-analiti£na ob£utljivostna analiza. Problem semi-analiti£ne ob£utljivostne analize je enak, kot je opisan pri metodi kon£nih diferenc. Numeri£na nestabilnost postopka, po- sebno v primerih iskane ob£utljivosti na oblikovne parametre problema, lahko vodi do velikih napak. Temu problemu je pozornost posvetilo veliko ²tudij, ki so natan£nost semi- analiti£ne ob£utljivostne analize posku²ale izbolj²ati z uporabo diferenc vi²jega reda ali z bolj temeljnimi modikacijami na nivoju kon£nih elementov (Keulen in Boer, 1998, Boer et al., 2002). Napredek v ra£unalni²ki znanosti je ponudil ²e dve moºnosti: odvajanje z uporabo simbol- nega in algebrai£nega sistema kot je npr. Mathematica (Mathematica 9.0) ter avtomatsko odvajanje (ang. automatic dierentiation ali computational dierentiation). Medtem ko odvajanje z uporabo splo²nega algebrai£nega sistema lahko vodi do problema nenadzoro- vane rasti izrazov, pri avtomatskem odvajanju tega problema ni. Avtomatsko odvajanje je zasnovano na dejstvu, da vsak ra£unalni²ki program izvaja zaporedje elementarnih aritmeti£nih operacij (se²tevanje, od²tevanje, mnoºenje ...) in elementarnih funkcij (log, sin, cos ...). Posledi£no se odvodi lahko izra£unajo avtomatsko s ponavljajo£o uporabo veriºnega pravila na te operacije, dobljene vrednosti pa so analiti£ne, to£ne do natan£- nosti ra£unalni²kega zapisa ²tevil (Griewank, 2000). Avtomatsko odvajanje je za izra£un ob£utljivostne analize lahko uporabljeno na globalnem nivoju problema ali pa na nivoju kon£nih elementov. Slabost uporabe avtomatskega odvajanja na globalnem nivoju je, da je aplikacija te tehnike neprakti£na in v primerih velikih komercialnih okolij za kon£ne elemente celo neizvedljiva. Na drugi strani aplikacija avtomatskega odvajanja na nivoju posameznega kon£nega elementa prav tako ni enostavna, saj so lahko prisotne odvisnosti med posameznimi spremenljivkami, ki iz izrazov v programski kodi na nivoju kon£nega elementa niso eksplicitno razvidne. Ta problem je re²il Korelc (2009a), ki je raz²iril teh- niko avtomatskega odvajanja z vpeljavo izjem in na ta na£in omogo£il uporabo avtomat- skega odvajanja tudi na nivoju kon£nih elementov. Formulacija inºenirskih problemov, ki temelji na uporabi raz²irjene tehnike avtomatskega odvajanja, se imenuje Automa- tic Dierentiation Based formulation (ali okraj²ano ADB formulacija) in je podrobneje predstavljena v poglavju 3.2. 1.2.4 Programska okolja za metodo kon£nih elementov, ki omogo£ajo aplika- cijo stohasti£nega pristopa V zadnjih dveh desetletjih se, skladno s pove£evanjem zanimanja za stohasti£ni pristop, po svetu razvija ve£ programskih kod, ki omogo£ajo stohasti£ni pristop. V ve£ini pred- stavljenih programov se poleg moºnosti ocene vpliva spremembe slu£ajnih parametrov na odziv konstrukcije (s £imer se ukvarjamo v disertaciji) veliko pozornost namenja predvsem analizi zanesljivosti konstrukcij, s katero se oceni verjetnost poru²itve konstrukcije. Oceno zanesljivosti se sicer lahko dolo£i tudi z uporabo simulacij Monte Carlo (ki jih omogo£ajo Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 9 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. prakti£no vsi predstavljeni programi) ali perturbacijske metode (oz. drugih podobnih response surface metod, ki odziv razvijejo v vrsto), obstajajo pa tudi aproksimativne metode za izra£un verjetnosti odpovedi konstrukcije. Najbolj znani sta metoda zaneslji- vosti prvega reda (ang. rst-order reliability method ali okraj²ano FORM) in metoda za- nesljivosti drugega reda (ang. second-order reliability method ali okraj²ano SORM), ki jih ima implementiranih ve£ina predstavljenih programov. Ve£ina predstavljenih programov omogo£a tudi izra£un ob£utljivostne analize, pri £emer se ob£utljivosti razen izjemoma izra£unajo z metodo kon£nih diferenc. Pomembnej²e razlike med programi se pokaºejo tudi pri tem, ali omogo£ajo tudi opis stohasti£nega polja ali le slu£ajnih spremenljivk (ki se preko domene problema ne spreminjajo) ter v numeri£ni u£inkovitosti glede na ²tevilo slu£ajnih parametrov, kompleksnost problema, ºeleno natan£nost rezultatov itd. NESSUS je programski paket, ki je bil prvotno zasnovan na Southwest Research in²titutu kot del Nasinega projekta (Thacker et al., 2006). NESSUS omogo£a simulacijo slu£ajnih spremenljivk, in sicer obteºbe, geometrije, materialnih parametrov in drugih, denira- nih s strani uporabnika. V program je vgrajenih ve£ razli£nih probabilisti£nih metod, npr. simulacije Monte Carlo, FORM, SORM, Monte Carlo s prioritetnim vzor£enjem (ang. importance sampling) in druge. NESSUS omogo£a tudi izvedbo ob£utljivostne analize z metodo kon£nih diferenc. Glavna prednost programa NESSUS je, da omogo£a tako izvedbo analize s pripravljenimi analiti£nimi modeli kot tudi povezavo z razli£nimi komercialnimi programskimi okolji (kot npr. ABAQUS, ANSYS, NASTRAN itd.). COSSAN je programski paket, katerega razvoj se je za£el leta 1992 na In²titutu inºenirske mehanike univerze Leopolda Franzensa v Innsbrucku (Avstrija). Avtorja programa sta Schuëller in Pradlwarter (2006). COSSAN sestavljata dve mnoºici komponent in sicer Stand Alone Tool Box ter Third Party Communication Tools. Prvi je zasnovan kot odprto splo²no namensko programsko orodje, ima predpripravljenih ve£ razli£nih modulov, razvr²£enih v posamezne skupine, in omogo£a re²evanje razli£nih mehanskih problemov, kot so optimizacija, modeliranje ²irjenja razpok, naklju£ne vibracije itd. Omogo£a tudi modeliranje stohasti£nih polj, vendar v predstavitvi programa ni podatka o na£inu in uporabljeni metodi za modeliranje le-teh. Komponenta Third Party Communication Tools ponuja moºnost komunikacije z zunanjimi kodami in omogo£a stohasti£no analizo brez spreminjanja izvorne kode deterministi£nih kon£nih elementov. COSSAN omogo£a tudi izvedbo ob£utljivostne analize z uporabo kon£nih diferenc. Na Kalifornijski univerzi so bili v okviru raziskovalnega dela ²tudentov na doktorskem ²tudiju razviti trije programi za probabilisti£no analizo: CalREL, FERUM in OpenSees (Kiureghian, Haukaas in Fujimura, 2006). Vsi trije programi so brez gra£nih vmesni- kov, zna£ilnih za komercialne programe, temve£ so razviti za namen raziskovanja in se kontinuirno razvijajo in dopolnjujejo. Prvi je bil razvit CalREL leta 1980, druga dva programa izhajata iz tega, zato imata ve£inoma iste lastnosti in moºnosti. CalREL je namenjen predvsem analizi zanesljivosti konstrukcij in omogo£a uporabo razli£nih metod, 10 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. npr. simulacije Monte Carlo, FORM, SORM, Monte Carlo s prioritetnim vzor£enjem itd. Pri tem omogo£a modeliranje medsebojno neodvisnih ali dolo£enih odvisnih slu£ajnih spremenljivk z ve£ razli£nimi funkcijami gostote verjetnosti (ki jih CalREL avtomati£no preslika v prostor standardnih normalnih slu£ajnih spremenljivk). Nadalje omogo£a ob- £utljivostno analizo, in sicer avtomati£no s strani programa po metodi kon£nih diferenc, lahko pa izraze za odvode uporabnik v izvorni kodi implementira sam (pristop z direktnim odvajanjem). Najve£ja razlika med vsemi tremi programi so programski jeziki, v kate- rih so kode napisane. CalREL je napisan v programskem jeziku FORTRAN, FERUM v Matlab skriptu (ang. Matlab scripting), OpenSees pa v C++. Prednost FERUM-a pred drugima dvema je tako predvsem v tem, da je napisan v programskem jeziku, ki je za£etnim uporabnikom laºje razumljiv. Slabost je slab²a numeri£na u£inkovitost. Glavna razlika OpenSees-a v primerjavi z drugima dvema pa je predvsem ta, da uporablja objek- tno orientiran pristop k modeliranju mehanskih problemov, kar pomeni, da ima vnaprej predpripravljene standardizirane izmenljive komponente, ki uporabniku olaj²ajo branje in uporabo v primerjavi z akademskimi programi, pri katerih je opis problema stvar razisko- valca in je zato lahko napisan nekonsistentno in nepregledno. Po drugi strani je prednost tega pristopa tudi njegova slabost, saj je zaradi tega bolj omejen in veliko manj eksibilen v primerjavi s sistemom, ki omogo£a popolno svobodo pri opisu problema in izpeljavi podprogramov. Sadcheva in Keane (2007) v svojem £lanku predstavita formulacijo projekcijske sheme z reducirano bazo (ang. reduced basis projection schemes) za re²evanje linearnih stohasti£- nih problemov z razvojem v sistem linearnih ena£b, ki je primerna za implementacijo v programska okolja za kon£ne elemente. V tem pristopu se izhaja iz semi-diskretiziranih parcialnih diferencialnih ena£b mehanskega problema. Re²itev problema se na ta na£in aproksimira z linearno kombinacijo neznanih koecientov in stohasti£nih baznih vektorjev, ki so izbrani z metodo stohasti£nih Krylovih podprostorov ustreznih dimenzij (Nair, 2001). Metoda je uporabna za stohasti£no polje. Za diskretizacijo le-tega avtorja uporabita ra- zvoj v K-L vrsto z uporabo Galerkinove metode, pri £emer ne predstavita podrobnosti modeliranja tega koraka. Omejitev predstavljene formulacije modeliranja stohasti£nih procesov je, da je uporabna le za linearne mehanske probleme. Reh et al. (2006) v £lanku predstavijo analizo probabilisti£nih kon£nih elementov z upo- rabo programa ANSYS. Razviti sta dve komponenti, ki omogo£ata stohasti£ni pristop: ANSYS Probabilistic Design System in ANSYS DesignXplorer. Prva komponenta ima enak gra£ni vmesnik kot ANSYS in je namenjena predvsem zahtevnej²im uporabnikom, medtem ko slednja temelji na objektno orientiranem pristopu in je namenjena predvsem parametri£nim ²tudijam. DesignXplorer omogo£a le modeliranje stohasti£nih spremen- ljivk, medtem ko Probabilistic Design System omogo£a tudi modeliranje stohasti£nega polja, vendar o podrobnostih o na£inu in metodi modeliranja le-tega ni podatka. Prednost ANSYS-a v primerjavi z drugimi programi, ki omogo£ajo stohasti£ni pristop, je moºnost Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 11 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. izra£una odvodov z avtomatskim odvajanjem. ANSYS je po na²ih podatkih edini program s stohasti£nim pristopom, ki omogo£a razvoj odziva v Taylorjevo vrsto vi²jega reda (kar pomeni perturbacijsko metodo vi²jega reda) z uporabo tehnike avtomatskega odvajanja. Avtomatsko odvajanje je aplicirano brez moºnosti uporabe izjem in podajanja podatka o implicitnih odvisnostih, ki iz kode kon£nega elementa niso razvidne, zaradi £esar je iz- vedljivost odvajanja omejena z dostopnostjo programske kode in odvisna od tipa vhodne spremenljivke, katero se odvaja, ter tipa parametra odziva, ki je predmet raziskave. Allaix in Carbone (2012) v svojem £lanku predstavita numeri£no orodje za diskretizacijo stohasti£nega polja v K-L vrsto. Program je napisan v jeziku FORTRAN 90. Za izra- £un lastnih vrednosti in lastnih vektorjev avtorja uporabita Galerkinovo proceduro, pri £emer je izra£un kovarian£ne matrike izveden s kombiniranjem dveh poddomen diskre- tizacijske mreºe zikalne domene na podoben na£in, kot je predstavljen v disertaciji (in smo ga premierno predstavili na konferenci ECCOMAS v Benetkah, Korelc in Melink, 2008). Program se v okviru diskretizacije stohasti£nega polja posve£a predvsem avtoma- tizaciji dolo£itve ustrezno goste mreºe, da je napaka variance zaradi diskretizacije polja ºeleno majhna. Poleg tega avtorja v £lanku predstavita prednosti opisa stohasti£nega odziva konstrukcije z uporabo perturbacijske metode vi²jega reda, vendar ni podatka o uporabljeni metodi za izra£un odvodov. Shang in Yun (2013) predstavita moºnost diskretizacije stohasti£nih polj v K-L vrsto z uporabo splo²nega programskega okolja za metodo kon£nih elementov ABAQUS. ƒlanek je posve£en predvsem predstavitvi uporabe razli£no gostih mreº za diskretizacijo zikalne domene za izra£un lastnih vrednosti in lastnih funkcij K-L vrste z uporabo Galerkinove metode ter za izvedbo izra£una odziva konstrukcije. Nadalje je predstavljena splo²na metoda preslikave in interpolacije izra£unanega stohasti£nega polja na mreºo za izvedbo izra£una odziva konstrukcije, ki je implementirana v predstavljeno programsko kodo. 1.3 Tema doktorske disertacije Glavni cilj doktorske disertacije je bil razviti splo²en simbolni pristop k avtomatizaciji mo- deliranja mo£no in srednje koreliranih stohasti£nih procesov v okolju za kon£ne elemente. Kot re£eno, je stohasti£ni pristop k modeliranju inºenirskih problemov zelo kompleksen za analizo. Zato je bilo glavno vodilo iskanje numeri£nih algoritmov in izbolj²av obsto- je£ih razvitih metod, ki bi bile kompatibilne z obstoje£imi numeri£nimi postopki, ki se uporabljajo v metodi kon£nih elementov. Metoda kon£nih elementov je v inºenirskih pro- blemih zaradi svoje robustnosti namre£ najpogosteje uporabljena metoda po svetu. Temu primerno je preizku²ena, ustrezno prilagojena in potrjeno natan£na za zelo ²irok spekter inºenirskih problemov. Razvoj avtomatizacije modeliranja stohasti£nih procesov, ki je kompatibilna za uporabo in analiziranje v okolju za kon£ne elemente, se zato ponuja kot najbolj²a izbira, vkolikor ºelimo, da je razvita avtomatizacija doprinos k ²ir²i uporabnosti 12 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. v inºenirskih problemih. Razvoj tak²ne avtomatizacije ima posledi£no enake prednosti ro- bustnosti in neomejenosti uporabe kot metoda kon£nih elementov sama, in sicer se jo da aplicirati na splo²ne geometrijsko in materialno nelinearne stati£ne mehanske probleme. Za razvoj avtomatizacije smo uporabili hibridni simbolno-numeri£ni pristop, ki zdruºuje program za simbolno ra£unanje, avtomatski generator kode AceGen (Korelc, 2009b) z vgrajeno tehniko avtomatskega odvajanja in programski jezik za numeri£ne operacije. Prednost tega pristopa je simbolni zapis ena£b kon£nega elementa, ki jih avtomatski generator kode AceGen nato prevede v poljubni izbrani programski jezik. Razviti kon£ni elementi so tako brez spreminjanja simbolnega zapisa uporabni v poljubnem izbranem programskem jeziku (npr. FORTRAN, C, Mathematica script, Matlab script itd.) in poljubnem okolju za metodo kon£nih elementov (npr. ABAQUS, AceFEM, FEAP itd.). V okviru diskretizacije stohasti£nega polja smo se posvetili razvoju formulacije stohasti£- nih kon£nih elementov, ki omogo£ajo operiranje s stohasti£nimi polji za poljubne topo- logije in poljubno geometrijo problema. Cilj je bil razviti formulacijo, ki je robustna in ki jo je moºno uporabiti v poljubnem splo²nem okolju za metodo kon£nih elementov. Za diskretizacijo stohasti£nega polja smo uporabili Karhunen-Loévovo dekompozicijo, saj le- ta predstavlja optimalno redukcijo v smislu potrebnega ²tevila £lenov vrste in je hkrati primerna za poljubno variabilne stohasti£ne procese. Za izra£un £lenov K-L dekompozi- cije smo zaradi njene aplikativnosti na metodo kon£nih elementov uporabili Galerkinovo metodo. Nadalje smo v okviru razvoja v K-L vrsto raziskovali moºnost za numeri£no poce- nitev tega postopka. Znotraj teh raziskav smo posebno pozornost posvetili poenostavljeni kovarian£ni funkciji, ki se z namenom pocenitve re²evanja v primeru srednje koreliranih stohasti£nih polj pogosto uporablja v literaturi. Dokazali smo, da poenostavljena kova- rian£na funkcija (predlagana v Phoon et al., 2002, Li et al., 2008, Stefanou et al., 2005 in Stefanou, 2009) ni pozitivno denitna. Posledi£no je problem numeri£no nestabilen. V doktorski disertaciji predstavimo dve alternativni modikaciji kovarian£ne funkcije, ki pocenita numeri£no re²evanje problema in hkrati zagotavljata bolj²o numeri£no stabilnost. Za realizacije stohasti£nih problemov smo uporabili simulacije Monte Carlo in perturba- cijsko metodo. V okviru slednje smo avtomatizacijo ob£utljivostne analize prvega reda, ki z uporabo avtomatskega odvajanja izra£una analiti£ne vrednosti odvodov, raz²irili za ob£utljivostno analizo (poljubnega) vi²jega reda. Na ta na£in dobimo analiti£ne vredno- sti odvodov vi²jega reda, ki omogo£ajo razvoj odziva v Taylorjevo vrsto vi²jega reda, kar pomeni perturbacijsko metodo vi²jega reda. Razvito ob£utljivostno analizo vi²jega reda, ki temelji na avtomatskem odvajanju, odlikujejo natan£nost, numeri£na u£inkovitost ter zelo kratek in strnjen zapis kode kon£nega elementa. Metoda izra£una analiti£ne vredno- sti odvodov ob visoki numeri£ni u£inkovitosti. Posledi£no je moºno izvesti ob£utljivostno analizo vi²jega reda za znatno ve£je ²tevilo ob£utljivostnih parametrov, kot je bilo izve- dljivo z drugimi metodami v preteklosti. Enako numeri£no u£inkovitost in natan£nost bi bilo moºno dose£i le ²e z ro£nim izpeljevanjem izrazov za izra£un odvodov znotraj Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 13 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. programske kode kon£nega elementa, vendar je ta postopek v primeru kompleksnej²ih problemov in vi²jega reda ob£utljivostne analize tako zamuden, da je vnos celotnih izra- zov v kodo kon£nih elementov prakti£no neizvedljiv in se ga zato ne uporablja, temve£ se za izra£un odvodov uporablja druge metode, ki pa so numeri£no draºje in manj zane- sljive. Celotna shema razvite avtomatizacije metode stohasti£nih kon£nih elementov je prikazana na sliki 1.1. 1. korak: DISKRETIZACIJA STOHASTIČNEGA POLJA Razvoj stohastičnega polja v Karhunen-Loèvovo vrsto. 2. korak: IZRAČUN STATISTIKE ODZIVA KONSTRUKCIJE MONTE CARLO PERTURBACIJSKA METODA SIMULACIJE Poljubno število realizacij ADB FORMULIRANA stohastičnega polja z SIMBOLNI METODA KONČNIH OBČUTLJIVOSTNA uporabo generatorja PRISTOP DIFERENC ANALIZA naključnih števil. Izračun vektorja odziva za vsako realizacijo Razvoj vektorja odziva v Taylorjevo vrsto okrog stohastičnega polja. pričakovanih vrednosti. Slika 1.1: Procedura metode stohasti£nih kon£nih elementov. Figure 1.1: Stochastic nite element procedure. 1.4 Zasnova doktorske disertacije Doktorska disertacija poleg uvodnega vsebuje ²e 8 poglavij: • 2. poglavje: Predstavitev teorije in temeljnih ena£b metod, ki so v disertaciji upo- rabljene za numeri£no re²evanje stohasti£nih problemov. • 3. poglavje: Predstavitev hibridnega simbolno-numeri£nega pristopa, uporabljenega v disertaciji. • 4. poglavje: Opis razvite avtomatizacije postopka diskretizacije stohasti£nega polja v K-L vrsto. 14 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. • 5. poglavje: Predstavitev problema numeri£ne nestabilnosti postopka ob uporabi aproksimirane K-L dekompozicije ter predlogi modikacije aproksimacij za izbolj- ²anje numeri£ne stabilnosti postopka. • 6. poglavje: Prikaz postopka izvedbe ob£utljivostne analize vi²jega reda po metodi direktnega odvajanja ter predstavitev postopka ob uporabi tehnike avtomatskega odvajanja. • 7. poglavje: Predstavitev rezultatov in u£inkovitosti predstavljene avtomatizacije diskretizacije stohasti£nega polja in perturbacijske metode vi²jega reda z uporabo tehnike avtomatskega odvajanja na ra£unskih primerih, v katerih sta stohasti£no modelirani geometrijska nepopolnost in korozija. • 8. poglavje: Zaklju£ki disertacije. • 9. poglavje: Povzetek vsebine doktorske disertacije. • 10. poglavje: Povzetek vsebine doktorske disertacije v angle²£ini. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 15 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 2 Numeri£no re²evanje stohasti£nih pro- blemov v mehaniki kontinuuma V tem poglavju bomo podrobneje predstavili teorijo in temeljne ena£be za vsakega izmed korakov stohasti£nega pristopa, ki smo ga privzeli v tej doktorski disertaciji in je prikazan na sliki 1.1. Predstavljene bodo tudi prednosti in slabosti ter obrazloºitve odlo£itve za vsako izmed izbranih metod. Kot je zapisano ºe v uvodu, nam je bilo temeljno vodilo pri naboru ustreznih metod razvoj splo²nega simbolnega pristopa k avtomatizaciji mode- liranja stohasti£nih procesov v okolju za kon£ne elemente, ki se da aplicirati na splo²ne geometrijsko in materialno nelinearne stati£ne mehanske probleme z mo£no in srednje koreliranimi stohasti£nimi polji. Zaradi tega je bila morda prezrta katera izmed metod, ki za speci£ne probleme daje natan£nej²e rezultate ob morebitni vi²ji numeri£ni u£inko- vitosti. Vendar se ob uporabi oºje omejenih metod izgubi eleganca re²itev, ki so robustne in uporabne neodvisno od specike in topologije problema. Vodilo k izboru primernih metod nam je bilo, da je na podro£ju stohastike pomanjkanje predvsem v smeri uporab- nih aplikacij stohasti£nih metod, ki bi omogo£ale re²evanje problemov na podoben na£in kot se z metodo kon£nih elementov re²uje splo²ne mehanske probleme. 2.1 Diskretizacija stohasti£nega polja Stohasti£no polje je denirano v dveh razli£nih prostorih: zikalnem prostoru ter prostoru slu£ajnih dogodkov. Narava prostora slu£ajnih dogodkov je tak²na, da se v tem prostoru stohasti£nega polja ne da diskretizirati z deljenjem na kon£no ²tevilo poddomen, kot se to naredi pri zikalni domeni v deterministi£ni analizi kon£nih elementov. Zato je pristop k diskretizaciji stohasti£nega polja druga£en kot pri obi£ajni metodi kon£nih elementov. Razvoj teorij, ki so ponudile temeljne koncepte teorije predstavitve stohasti£nih polj, sega v sredino prej²njega stoletja. V tem obdobju se je razvilo ve£ pristopov od katerih jih ve£ina razvije slu£ajno polje v vrsto. Naj bo wpX, θq stohasti£no polje, denirano na d-dimenzionalni zikalni domeni D Ă d R , kjer je X vektor prostorskih koordinat in θ dogodek iz prostora slu£ajnih dogodokov Ω. Stohasti£no polje je dolo£eno s pri£akovano vrednostjo polja ¯ wpXq in kovarian£no funkcijo CpX1, X2q. Kovarian£na funkcija je funkcija, ki dolo£a medsebojno koreliranost vrednosti preko stohasti£nega polja. Stohasti£no polje se lahko razvije v vrsto 16 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. M ÿ wpX, θq “ NkpXqξkpθq “ NpXq ξpθq, (2.1) k“1 kjer so ξk slu£ajne spremenljivke, Nk funkcije (najpogosteje so to oblikovne funkcije kon£nih elementov), ki jih zdruºimo v vektorja NpXq “ tN1pXq, N2pXq, ..., NMpXqu ter ξpθq “ tξ1pθq, ξ2pθq, ..., ξM pθquT . Pri£akovana vrednost in kovarian£na funkcija aproksimiranega polja sta ¯ w T pXq “ NpXq ¯ w C T pX1, X2q “ NpX1q C NpX2q, (2.2) kjer je ¯w “ Epξpθqq vektor pri£akovanih vrednosti in C kovarian£na matrika slu£ajnega vektorja ξpθq. V preteklosti je bilo razvitih ve£ metod, ki se razlikujejo v izbiri funkcij Nk in slu£ajnih spremenljivk ξk v ena£bi (2.1). Nekatere izmed njih so predstavljene v naslednjih vrsticah (povzeto po Keese, 2003). Interpolacijska metoda ali metoda oblikovnih funkcij (Liu et al., 1986a, 1986b) interpolira stohasti£no polje z oblikovnimi funkcijami kon£nih elementov NkpXq v vozli²£ih Xj. Ker so oblikovne interpolacijske funkcije Nk obi£ajno vzete tako, da so v vozli²£u k enake 1, v preostalih vozli²£ih pa 0, iz tega sledi ξk “ ξpXkq. Pri£akovana vrednost in kovarian£na matrika iz izraza (2.2) sta v tem primeru ¯ wk “ ¯ wpXkq Ckj “ CpXk, Xjq. (2.3) Metoda sredi²£ne to£ke je poseben primer metode oblikovnih funkcij, pri kateri so NkpXq odsekoma konstantne in kjer je Xk sredi²£na to£ka kon£nega elementa. Metoda lokalnega povpre£enja (Vanmarcke, 1983) je metoda, pri kateri se prav tako upo- rabi odsekoma konstantne funkcije NkpXq, le da se po tej metodi ξk izbere kot prostorsko povpre£je stohasti£nega polja wpX, θq preko domene k-tega kon£nega elementa (Dk) z volumnom ş |Vk|: ξkpθq “ 1 wpX, θqdX. Pri£akovana vrednost in kovarian£na matrika |Vk| Dk se v tem primeru izra£unata po izrazu 1 ż ¯ wk “ ¯ wpXqdX |Vk| Dk 1 ż ż Ckj “ CpXk, XjqdXkdXj. (2.4) |Vk| |Vj| Dk Dj Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 17 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Zaradi povpre£enja to£ne statistike metoda lokalnega povpre£enja nekoliko podceni vari- abilnost stohasti£nega polja. Kot je bilo predstavljeno ºe v uvodnih poglavjih, je glavna slabost zgoraj na²tetih metod ta, da je ²tevilo slu£ajnih spremenljivk sorazmerno ²tevilu kon£nih elementov, s kate- rimi je diskretizirana zikalna domena problema. Poleg tega so slu£ajne spremenljivke medsebojno korelirane. Z gostenjem mreºe kon£nih elementov se koreliranost slu£ajnih spremenljivk vi²a, kar posledi£no pomeni kovarian£no matriko blizu singularnosti. Ta problem je bil re²en z razvojem metod, ki so stohasti£no polje razvile v ortogonalno vrsto. Metode razvoja v ortogonalno vrsto za funkcije Nk izberejo medsebojno ortogonalne funk- cije. Posledi£no so slu£ajne spremenljivke medsebojno nekorelirane 1 ż ξkpθq “ wpX, θqNkpXqdX. (2.5) ||Nk|| D Pri£akovana vrednost in kovarian£na matrika se v tem primeru izra£unata po izrazih ż ¯ w “ ¯ wpXqNpXqdX D ż ż C T “ NpX1qCpX1, X2qNpX2q dX1dX2. (2.6) D D Glavna omejitev uporabnosti metod ortogonalnega razvoja v vrsto je, da so primerne le za stohasti£na polja z Gaussovo porazdelitvijo. Izmed metod razvoja v ortogonalno vrsto je v smislu srednje kvadratne napake (ang. mean square error) pri omejenem ²tevilu £lenov v vrsti ponudila najbolj²o re²itev Karhunen-Loèvova (K-L) dekompozicija, ki so jo leta 1947 neodvisno razvili Karhunen ter Kac in Siegert ter leta 1948 Loève in je predstavljena v veliko literaturah, npr. Ghanem in Spanos, 2003, Huang, Quek in Phoon, 2001, Stefanou in Papadrakakis, 2007. K-L dekompozicija je spektralna reprezentacija stohasti£nega po- lja, pri kateri so za funkcije Nk izbrane lastne funkcije homogene Fredholmove integralske ena£be drugega reda. Glavna omejitev uporabnosti K-L dekompozicije je, kot pri ostalih metodah razvoja v ortogonalno vrsto, da je primerna le za stohasti£na polja z Gaussovo porazdelitvijo. Za primere drugih, ne-Gaussovih porazdelitev, so bili sicer razviti pristopi, kako aplicirati razvoj teh procesov v K-L vrsto (npr. Schevenels et al., 2004), vendar se bomo v doktorski disertaciji omejili le na Gaussova stohasti£na polja in se bomo opisu teh pristopov zato na tem mestu izognili. Ker je zaradi svojih prednosti K-L vrsta izbrana za diskretizacijo stohasti£nega polja v disertaciji, je podrobneje predstavljena v naslednjem poglavju. 18 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 2.1.1 Diskretizacija stohasti£nega polja v Karhunen-Loèvovo vrsto Stohasti£no polje wpX, θq se lahko razvije v K-L vrsto, ki se jo za potrebe aplikativnosti diskretizacije odreºe in obdrºi le prvih M £lenov vrste (Ghanem in Spanos, 2003) 8 M ÿ a ÿ a wpX, θq “ ¯ wpXq ` λkfkpXqξkpθq « ¯ wpXq ` λkfkpXqξkpθq, (2.7) k“1 k“1 kjer je ¯ wpXq pri£akovana vrednost stohasti£nega polja, λk in fkpXq so lastne vrednosti in lastne funkcije, θ je dogodek iz prostora slu£ajnih dogodokov Ω, ξkpθq pa so normi- rane nekorelirane slu£ajne spremenljivke. V primeru Gaussovega stohasti£nega polja, so ξkpθq spremenljivke z Gaussovo porazdelitvijo z ni£elno pri£akovano vrednostjo ter enotsko varianco ali t. i. standardizirano normalno porazdelitvijo. λk in fkpXq so re²itev homo- gene Fredholmove integralske ena£be drugega reda v kateri je jedro kovarian£na funkcija CpX1, X2q ż CpX1, X2qfkpX1qdX1 “ λkfkpX2q. (2.8) D Integralske ena£be Fredholmovega tipa so podrobneje opisane v knjigi Mikhlina (1957). Ena£ba (2.8) se imenuje homogena, ker je t.i. prosti £len Fredholmove ena£be v tem primeru enak ni£. Kovarian£na funkcija je po deniciji omejena, simetri£na in pozitivno denitna. λk v tej ena£bi, ki odgovarjajo netrivialni re²itvi ena£be, so karakteristi£ne vrednosti ali lastne vrednosti, fk pa so karakteristi£ne funkcije ali lastne funkcije. Lastne funkcije tvorijo kompletno mnoºico ortonormiranih funkcij, tako da velja ż fkpXqflpXqdX “ δkl, (2.9) D kjer je δkl Kroneckerjev simbol. Lastne vrednosti predstavljajo magnitudo lastnih funk- cij. Ker je kovarian£na funkcija pozitivno denitna, so tudi lastne vrednosti pozitivne in £e jih uredimo v padajo£em redu λ1 ě λ2 ě ... ą 0, ima vsak naslednji £len K-L vrste manj²o ali enako magnitudo in zato vpliv vsakega naslednjega £lena na obna²anje stohasti£nega polja v splo²nem pada. Kot je znano iz funkcionalne analize, hitreje kot pada kovarian£na funkcija proti ni£ (kar pomeni, da je stohasti£no polje manj korelirano), ²ir²a je odgovarjajo£a sprektralna gostota, kar pomeni da je za dovolj natan£no pred- stavitev stohasti£nega polja potrebnih ve£ £lenov K-L dekompozicije. Ali obratno, bolj kot je stohasti£no polje korelirano, hitreje padajo vrednosti λi z nara²£ajo£im indeksom i. Lastne vrednosti z vi²jim indeksom imajo zato v tem primeru bistveno manj²e (do nekaj velikostnih redov) vrednosti v primerjavi z lastnimi vrednostmi z za£etnimi indeksi in zato znatno manj²i vpliv na odziv polja. Stohasti£no polje je zato v primeru visoke koreliranosti dovolj natan£no opisano ºe z nekaj za£etnimi £leni K-L vrste. V skrajnem primeru, ko je CpX1, X2q “ σ2, zadostuje samo en £len vrste. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 19 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Iz (2.7) z upo²tevanjem ortogonalnosti sledi spektralna dekompozicija kovarian£ne funkcije 8 ÿ CpX1, X2q “ λkfkpX1qfkpX2q. (2.10) k“1 Kovarian£na funkcija in njene lastnosti bodo podrobneje predstavljene v poglavju 5. V disertaciji se omejimo na stohasti£ne procese drugega reda (ang. second order random process). Stohasti£ni proces drugega reda je proces, za katerega velja, da je pri£akovana vrednost Epw2q ă `8 za vsak slu£ajni dogodek iz prostora slu£ajnih dogodkov. To pomeni, da so pri£akovana vrednost, korelacija in kovarian£na funkcija tega procesa dobro denirane in kon£ne. S to zahtevo ne izgubimo splo²nosti, saj je velika ve£ina zikalno merljivih stohasti£nih procesov, s katerimi se sre£amo pri mehanskih problemih, drugega reda. Za bolj²o predstavo, kak²na je povezava med ena£bama (2.8) in K-L dekompozicijo, si poglejmo naslednjo izpeljavo. Zapi²imo kovarian£no funkcijo za stohasti£no polje wpX, θq: CpX1, X2q “ E pwpX1, θqwpX2, θqq É pwpX1, θqq E pwpX2, θqq (2.11) Ob razvoju stohasti£nega polja v vrsto (ena£ba (2.7)), sledi ˜˜ 8 ¸ ˜ 8 ¸¸ ÿ a ÿ a CpX1, X2q “ E ¯ wpX1q ` λkfkpX1qξkpθq ¯ wpX2q ` λmfmpX2qξmpθq k“1 m“1 ˜ 8 ¸ ˜ 8 ¸ ÿ a ÿ a É ¯ wpX1q ` λkfkpX1qξkpθq E ¯ wpX2q ` λmfmpX2qξmpθq k“1 m“1 ˜ 8 8 ¸ ÿ a ÿ a “ ¯ wpX1q ¯ wpX2q È λkfkpX1qξkpθq λmfmpX2qξmpθq ´ k“1 m“1 ´ ¯ wpX1q ¯ wpX2q 8 8 ÿ ÿ a “ Epξkpθqξmpθqq λkλmfkpX1qfmpX2q. (2.12) k“1 m“1 Ker sta ξkpθq in ξmpθq slu£ajni spremenljivki z Gaussovo porazdelitvijo z ni£elno pri£ako- vano vrednostjo in enotsko varianco, je bilo v izpeljavi upo²tevano Epξkpθqq “ 0. (2.13) Pri£akovana vrednost produkta dveh slu£ajnih spremenljivk s standardizirano normalno 20 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. porazdelitvijo je Epξkpθqξmpθqq “ δkm. (2.14) Iz tega sledi 8 ÿ a CpX1, X2q “ λkλkfkpX1qfkpX2q. (2.15) k“1 Ker so lastne vrednosti pozitivna realna ²tevila, velja ?λkλk “ λk. Sedaj pomnoºimo obe strani ena£be (2.15) z fnpX2q in integriramo po domeni D. ƒe nato upo²tevamo ortogonalnost lastnih funkcij (en. (2.9) ), dobimo ż ż 8 ÿ CpX1, X2q fnpX2qdX2 “ λkfkpX1qfkpX2qfnpX2qdX2 D D k“1 8 ż ÿ “ λkfkpX1q fkpX2qfnpX2qdX2 (2.16) k“1 D “ λnfnpX1q, kar je to£no homogena Fredholmova integralska ena£ba druge vrste. Poglejmo si ²e dokaz, da je posplo²eni koordinatni sistem, ki ga denirajo lastne funkcije, optimalen v smislu srednje kvadratne napake, ki nastopi zaradi zanemaritve vi²jih £lenov K-L dekompozicije. Naj bo hk pXq kompletna ortonormalna mnoºica funkcij. Stohasti£no polje wpX, θq z ni£elno pri£akovano vrednostjo ¯ wpXq “ 0 aproksimiramo s konvergen£no vrsto oblike 8 ÿ wpX, θq “ λkhkpXqξkpθq. (2.17) k“1 Odreºimo vrsto (2.17) po £lenu M. Napaka zaradi zanemaritve £lenov je 8 ÿ M “ λkhkpXqξkpθq. (2.18) k“M `1 Ena£bo (2.17) pomnoºimo s hlpXq, integriramo po domeni D in upo²tevamo ortogonalnost mnoºice funkcij hkpXq. ƒe ena£bo preuredimo tako, da izrazimo ξpθq, sledi 1 ż ξlpθq “ wpX, θqhlpXqdX. (2.19) λl D Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 21 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Sedaj £len ξkpθq v (2.18) zamenjamo z (2.19) ter zapi²emo srednjo kvadratno napako 8 8 ż ż ÿ ÿ 2 h E M “ k pXqhlpXq pwpX1, θq wpX2, θqqhkpX1qhlpX2qdX1dX2 k“M `1 l“M `1 D D 8 8 ż ż ÿ ÿ “ hkpXqhlpXq CpX1, X2qhkpX1qhlpX2qdX1dX2. (2.20) k“M `1 l“M `1 D D V naslednjem koraku integriramo ena£bo (2.20) po domeni D. Ob upo²tevanju ortogo- nalnosti mnoºice hnpXq sledi ż 8 ż ż ÿ 2 dX C M “ pX1, X2qhkpX1qhkpX2qdX1dX2. (2.21) D k“M `1 D D Da bo napaka minimizirana, zahtevamo da je ortogonalna na prostor, ki ga razpenjajo funkcije hnpXq, kar pomeni da bo re²itev minimizirala funkcional 8 ż ż ÿ F phkpXqq “ CpX1, X2qhkpX1qhkpX2qdX1dX2 k“M `1 D D ˆż ˙ ´λk hk pXq hk pXq dX ´ 1 . (2.22) D Ker zahtevamo, da je funkcional minimiziran, ga odvajamo po hipXq in i²£emo re²itev, ko je odvod enak 0 ˆ ˙ BF ph ż ż k pXqq “ CpX1, X2qhipX1qdX1 ´ λihipX2q dX2 “ 0, (2.23) BhipXq D D kar je izpolnjeno v primeru, ko je ż CpX1, X2qhipX1qdX1 “ λihipX2q. (2.24) D 2.1.2 Galerkinova procedura za numeri£no re²evanje homogene Fredholmove integralske ena£be drugega reda Fredholmova integralska ena£ba (2.8) je analiti£no re²ljiva le za dolo£ene kovarian£ne funkcije in za enostavnej²e oblike domen. Izpeljava re²itve za enega takih primerov je podana npr. v Ghanem in Spanos, 2003. V splo²nem je Fredholmovo integralsko ena£bo potrebno re²iti z numeri£nim postopkom. Za aplikacijo metode stohasti£nih kon£nih ele- mentov znotraj standardnega okolja za metodo kon£nih elementov je za numeri£no re²eva- nje Fredholmove ena£be najprimernej²a Galerkinova procedura. Po tej proceduri se vsaka 22 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. izmed lastnih funkcij aproksimira z linearno kombinacijo J oblikovnih funkcij NjpXq J ÿ fkpXq « fjkNjpXq, k “ 1, 2, ..., M (2.25) j“1 kjer so fjk neznani koecienti. Ena£bo (2.25) ustavimo v (2.8) in zamenjamo vrstni red se²tevanja in integriranja J ˆż ˙ ÿ fjk CpX1, X2qNjpX1qdX1 ´ λkNjpX2qdX2 « 0. (2.26) j“1 D Leva in desna stran ena£be (2.26) nista povsem enaki, zaradi aproksimacije lastnih funkcij (en. 2.25)). Z uporabo Galerkinovega pristopa, se zahteva, da je napaka zaradi aproksi- macije otrogonalna na prostor, ki ga razpenjajo oblikovne funkcije. Sledi J ˆż ż ż ˙ ÿ fjk CpX1, X2qNjpX1qNipX2qdX1dX2 ´ λk NjpX2qNipX2qdX2 “ 0(2.27) j“1 D D D Ob uporabi oznak ż ż Cij “ CpX1, X2qNjpX1qNipX2qdX1dX2, D D ż Nij “ NjpX2qNipX2qdX2, (2.28) D Λlk “ δlkλk, ena£bo (2.27) lahko preoblikujemo v matri£no obliko Cf “ ΛNf . (2.29) Ena£ba (2.29) predstavlja posplo²en problem lastnih vrednosti. Re²itev te ena£be so matrika lastnih vektorjev f in lastne vrednosti λk. Lastne funkcije je potrebno dodatno normirati skladno z ena£bo (2.9). Za bolj²o predstavo, kako izgledajo lastne funkcije, so na sliki 2.1 prikazane prve ²tiri la- stne funkcije za primer enodimenzionalnega stohasti£nega polja, deniranega na intervalu r0, 1s, s pri£akovano vrednostjo 0, standardno deviacijo 1 ter eksponentno kovarian£no funkcijo z enotsko korelacijsko dolºino lc }X2´X1} CpX l 1, X2q “ σ2eć “ e´}X2´X1}. (2.30) Kot je razvidno na sliki, se z vsakim vi²jim indeksom pove£a ²tevilo polvalov lastne funkcije. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 23 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. fi 1.0 f1 0.5 f2 x f3 0.2 0.4 0.6 0.8 1.0 f4 - 0.5 - 1.0 Slika 2.1: Prve ²tiri lastne funkcije za enodimenzionalno stohasti£no polje z ni£elno pri£a- kovano vrednostjo in enotsko standardno deviacijo ter eksponentno kovarian£no funkcijo. Figure 2.1: First four eigenfunctions for one-dimensional stochastic eld with zero mean, unit variance and exponential covariance function. Slika 2.2 pa prikazuje prvih 10 lastnih vrednosti λk, k “ 1, ..., 10 za enodimenzionalno stohasti£no polje, denirano na intervalu r0, 1s, s pri£akovano vrednostjo 0, standardno deviacijo σ “ 1 ter eksponentno kovarian£no funkcijo za tri razli£no korelirana stohasti£na polja. Koreliranost stohasti£nega polja je dolo£ena s korelacijsko dolºino lc. Ve£ja kot je korelacijska dolºina, bolj so vrednosti preko polja povezane in, kot je razvidno tudi na sliki 2.2, hitreje padajo lastne vrednosti. Da bi bralec dobil bolj²i ob£utek, kaj pomeni koreliranost stohasti£nega polja, so na sliki 2.3 prikazane ²tiri naklju£no izbrane realizacije dvodimenzionalnega stohasti£nega polja z ni£elno pri£akovano vrednostjo in enotsko standardno deviacijo ter korelacijskimi dol- ºinami lc “ 0.1; 0.5; 2.0 in 8. Kot je razvidno s slike, je v primeru ve£je koreliranosti stohasti£no polje bolj gladko, medtem ko vrednosti stohasti£nega polja pri niºji korelira- nosti opazno bolj variirajo. V nadaljevanju bomo namenili nekaj pozornosti izbiri gostote diskretizacijske mreºe za uporabo Galerkinove procedure ter izbiri potrebnega ²tevila £lenov, ki jih obdrºimo v K-L dekompoziciji. V deterministi£ni metodi kon£nih elementov se velikost kon£nih elemen- tov obi£ajno dolo£a glede na pri£akovani gradient napetosti in geometrijo problema. Pri diskretizaciji stohasti£nega polja pa je gostota mreºe odvisna predvsem od koreliranosti stohasti£nega polja. Kot je razvidno na sliki 2.2, v splo²nem velja: bolj kot je stohasti£no polje korelirano, hitreje padajo lastne vrednosti. Spomnimo se, da lastne vrednosti pred- stavljajo magnitudo oz. uteºi lastnih funkcij. Ker v primeru visoko koreliranega polja vrednosti lastnih vrednosti z vi²anjem indeksa sorazmerno hitro padajo, to pomeni da imajo prvi £leni K-L vrste sorazmerno ve£ji vpliv v primerjavi s £leni z vi²jimi indeksi. V primeru nizko koreliranih stohasti£nih polj pa lastne vrednosti padajo zelo po£asi, kar 24 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Λ i æ 0.8 0.6à æ lc=2.0 à lc=0.5 0.4 ì lc=0.1 à ì 0.2 ì ì æ ì à ì ì à ì æ à ì ì æ ì æ à æ æ à æ à æ à æ à i 2 4 6 8 10 Slika 2.2: Prvih deset lastnih vrednosti za enodimenzionalno stohasti£no polje z ni£elno pri£akovano vrednostjo in enotsko standardno deviacijo ter eksponentno kovarian£no funk- cijo s tremi razli£nimi korelacijskimi dolºinami. Figure 2.2: First ten eigenvalues for one-dimensional stochastic eld with zero mean, unit variance and exponential covariance function for three dierent correlation lengths. pomeni, da imajo tudi £leni K-L vrste z vi²jim indeksom ²e vedno pomemben doprinos k obna²anju polja, vrednosti stohasti£nega polja pa preko domene mo£neje variirajo (slika 2.3 (a)). Zato mora biti v tem primeru mreºa elementov gostej²a, da lahko zajame bistvene zna£ilnosti slu£ajnega polja. Narejenih je bilo ºe ve£ ²tudij, v katerih se je raziskovalo ustrezno gostoto mreºe za opis stohasti£nega polja, iz katerih lahko povle£emo splo²no oceno, da mora biti velikost kon£nih elementov tak²na, da sta vzdolº ene korelacijske dol- ºine vsaj 2 kon£na elementa. Natan£nej²a ocena je, da naj bi bila razdalja med vozli²£i kon£nih elementov manj²a od ene £etrtine vala lastne funkcije z najvi²jim indeksom, pri kateri (glede na izbrano ºeleno natan£nost) odreºemo K-L vrsto (povzeto po Matthies et al., 1997). Natan£nost diskretizacije stohasti£nega polja se lahko oceni s pomo£jo variance stohasti£- nega polja wpX, θq, ki je denirana kot 8 Var ÿ pwpX, θqq “ λif 2ipXq (2.31) i“1 ter variance aproksimiranega stohasti£nega polja, ki je denirana kot řM λ i“1 if 2 i pXq. Lo- kalna napaka diskretizacije je tako ř8 λ řM if 2 λif 2 err i“1 i pXq í“1 i pXq w pXq “ . (2.32) ř8 λ i“1 if 2 i pXq Ocena globalne diskretizacijske napake je formulirana v smislu povpre£ne napake stoha- Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 25 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. (a) (b) (c) (d) Slika 2.3: Primerjava realizacij dvodimenzionalnega stohasti£nega polja z ni£elno pri£a- kovano vrednostjo, enotsko standardno deviacijo ter eksponentno kovarian£no funkcijo za (a) lc “ 0.1, (b) lc “ 0.5, (c) lc “ 2.0 in (d) lc Ñ 8. Figure 2.3: Comparison realizations of stochastic eld with zero mean, unit variance and exponential covariance function with (a) lc “ 0.1, (b) lc “ 0.5, (c) lc “ 2.0 and (d) lc Ñ 8. sti£nega polja 1 ż errw “ errwpXqdX, (2.33) V D kjer je V volumen (oz. v primeru dvodimenzionalnega problema povr²ina, v primeru enodimenzionalnega problema pa dolºina) domene. V praksi je pogosteje uporabljena numeri£no nekoliko cenej²a ocena napake, ki zanemari, da se varianca stohasti£nega polja, diskretiziranega s K-L vrsto, dejansko spreminja preko domene. V tem primeru je VarpwpX, θqq “ σ2, kjer je σ2 varianca stohasti£nega polja. Ocena napake je tako 1 ż σ2 řM řM ´ λif 2 λi err i“1 i pXq i“1 w “ dX “ 1 ´ . (2.34) V σ2 V σ2 D Izraz za desnim ena£ajem sledi iz ena£b (2.9) in (2.31). 26 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Ker je Galerkinova procedura numeri£na metoda re²evanja Fredholmovega integrala, je na tem mestu prav omeniti, da ima kot vse numeri£ne metode tudi ta metoda diskretizacijsko napako. Dobljene lastne vrednosti so manj²e ali enake analiti£nim ter so to£nej²e lastne vrednosti kot lastne funkcije (Ghanem in Spanos, 2003). To pomeni, da aproksimacija K- L vrste, dobljena z Galerkinovo metodo, izgubi del natan£nosti, zaradi katere je K-L vrsta sicer superiorna nad ostalimi razvoji v vrsto. Li in Kiureghian (1993) v £lanku, v katerem predstavita metodo razvoja v vrsto z linearnim optimalnim ocenjevanjem, trdita, da so aproksimacije, ki z vpeljavo oblikovnih funkcij re²ujejo razvoj v vrsto preko diskretnega problema lastnih vrednosti, med katere spada tudi predstavljena Galerkinova procedura za re²evanje Fredholmove ena£be, v smislu napake variance celo slab²e od njune razvite metode, ki uporablja optimalne oblikovne funkcije. Omeniti velja ²e, da je napaka K- L vrste nekoliko ve£ja v robnih obmo£jih domene, kar je zna£ilnost metod spektralne dekompozicije. Prav tako K-L vrsti ni v prid, da je homogenost stohasti£nega polja, opisanega s K-L vrsto vpra²ljiva, saj se varianca preko domene nekoliko spreminja. Kljub navedenim slabostim metode, smo se v doktorski disertaciji odlo£ili za K-L vrsto, saj je v primerjavi z drugimi razvitimi metodami ²e vedno najve£ argumentov v prid K-L vrsti. Najpomembnej²a prednost je, da od vseh razvitih metod K-L vrsta najbolj reducira ²tevilo dimenzij stohasti£nega prostora ter je posplo²eni koordinatni sistem, ki ga denirajo lastne funkcije, optimalen v smislu srednje kvadratne napake. Poleg tega Galerkinova procedura z gostenjem mreºe konvergira k analiti£nim re²itvam, kar z drugimi besedami pomeni, da se ob dovolj gosti mreºi s to metodo ²e vedno lahko poljubno pribliºamo optimalni re²itvi. 2.2 Diskretni sistemi v mehaniki kontinuuma Za celoten vpogled v re²evanje mehanskih problemov s perturbacijsko metodo je prav, da se na za£etku najprej predstavi osnove re²evanja problemov v mehaniki kontinuuma. Teo- rija kontinuuma zanemari, da je snov sestavljena iz delcev (npr. atomov na mikroskopski ravni ali diskontinuitet, vklju£kov, razpok itd. na makroskopski ravni), zato je natan£na le za primere, ko so dimenzije mehanskih problemov veliko ve£je v primerjavi z razda- ljami med temi delci. Mehanika kontinuuma se deli na mehaniko teko£in in mehaniko trdnih teles. V doktorski disertaciji se ukvarjamo s primeri, ki sodijo v mehaniko trdnih teles, vendar se ideja in pristop lahko prenese tudi na mehaniko teko£in. V disertaciji bomo osnove mehanike trdnin predstavili kolikor se da na kratko in zgo²£eno. Za poglo- bljene ²tudije naj se bralec obrne npr. na Wriggers 2008, Zienkiewicz in Taylor (1991) ali Criseld (1996). Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 27 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 2.2.1 Osnovne ena£be v mehaniki kontinuuma Trdno telo je denirano na nekem obmo£ju D, ki ga omejuje rob obmo£ja BD. Na posame- znih delih roba obmo£ja so lahko predpisani pomiki ali obteºba. Zanima nas odziv telesa (pomik, deformacije, napetosti itd.), do katerega pride zaradi delovanja sil in predpisanih pomikov. Obna²anje teles se opi²e s tremi sklopi ena£b: kinemati£nimi, konstitutivnimi in ravnoteºnimi. Ena£be so lahko zapisane glede na za£etno ali trenutno konguracijo teles. V spodnjih vrsticah je ob vsakem sklopu ena£b zapisana mo£na oblika teh ena£b za stati£en problem in hiperelasti£en material glede na za£etno konguracijo. V tem pri- meru imamo ve£ moºnosti opisa mer napetosti, t. i. prve Piola-Kirchoove napetosti ali druge Piola-Kirchoove napetosti, kar vodi do dveh razli£nih formulacij (prva moºnost je prikazana v levem, druga pa v desnem stolpcu): • Kinemati£ne ena£be, s katerimi se opi²e deformacije in gibanje telesa, mere defor- macij in £asovne odvode kinemati£nih koli£in F E “ 1pFT F ´ 1q, 2 kjer je F deformacijski gradient in E Green-Lagrangev deformacijski tenzor. • Konstitutivne ena£be, s katerimi se opi²e materialne lastnosti telesa (npr. odnos napetost-deformacija) P “ BW S “ BW , BF BE kjer je P prvi Piola-Kirchoov napetostni tenzor, S drugi Piola-Kirchoov napeto- stni tenzor in W funkcija, ki denira elasti£no energijo. • Ravnoteºne ena£be, s katerimi se zagotovijo zakoni ravnoteºja (ravnoteºje sil in momentov, ohranjanje mase, prvi zakon termodinamike) DivP ` ρ ¯ ¯ 0b “ 0 DivpF Sq ` ρ0b “ 0, kjer ρ ¯ 0b denira volumsko silo. Tako deniran problem, pri katerem moramo re²iti ena£be ob upo²tevanju danih robnih pogojev (predpisanih pomikov in predpisane obteºbe), se imenuje robni problem. Ne- znanke problema, ki nas zanimajo, so obi£ajno pomiki u. V disertaciji se omejimo na me- hanske probleme, pri katerih je temperatura konstantna in je termodinamsko ravnovesje 28 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. med telesi in okolico avtomati£no izpolnjeno. Za re²itev robnega problema se kinema- ti£ne ena£be in ena£be snovi vnese v ravnoteºne ena£be. Rezultat je v splo²nem sistem diferencialnih ena£b, tak opis problema se imenuje mo£na oblika robnega problema. Re²itev mo£ne oblike robnega problema izpolnjuje pogoj ravnoteºja v vsaki to£ki telesa. Ker so diferencialne ena£be v taki obliki re²ljive le za dolo£ene enostavnej²e oblike telesa in dolo£ene robne pogoje, so se v preteklosti razvile aproksimativne metode, npr. metoda kon£nih diferenc, Ritzova metoda, Galerkinova metoda in metoda kon£nih elementov, s katerimi se robne probleme re²uje aproksimativno z numeri£nimi metodami (Zienkiewicz in Taylor, 1991). Izmed numeri£nih metod je zaradi svoje splo²nosti in numeri£ne stabilnosti v inºenirskih problemih najpogosteje uporabljena metoda kon£nih elementov. Metoda kon£nih elemen- tov prevede robni problem na variacijski problem. Po metodi kon£nih elementov se mo£no obliko ena£b preoblikuje v ²ibko obliko, kar vodi v formulacijo minimiziranja napake, ki nastane zaradi aproksimacije s kon£nimi elementi. Pri formulaciji ²ibke oblike ravnoteºnih ena£b dejanski pomik u, ki predstavlja to£no re²itev ravnoteºnih ena£b, aproksimiramo s pomikom uh. DivPpu ¯ hq ` ρ0b “ Rh, (2.35) kjer rezidual Rh predstavlja napako zaradi aproksimacije pomikov. Z namenom minimizacije napake zaradi aproksimacije nato zahtevamo, da je napaka enaka 0 v ²ibkem smislu in sicer tako, da se rezidual skalarno pomnoºi z uteºno funkcijo η ter njun produkt integrira preko celotne domene. Funkcija η mora izpolnjevati predpisane robne pogoje pomikov in se imenuje virtualni pomik ali testna funkcija. Re²ujemo torej ena£bo ż ż ż R ¯ h ¨ η dV “ 0 ùñ DivPpuhq ¨ η dV ` ρ0b ¨ η dV “ 0. (2.36) D D D Ena£ba (2.36) je izpolnjena tudi za to£en pomik u. ż ż DivP ¨ η dV ` ρ ¯ 0b ¨ η dV “ 0. (2.37) D D Z integracijo per partes prvega £lena ena£be (2.37), z upo²tevanjem divergen£nega izreka in vpeljavo robnega pogoja predpisane obteºbe se ²ibko obliko ravnoteºne ena£be lahko zapi²e kot ż ż ż P ¨ Grad η dV ´ ρ ¯ ¯ 0b ¨ η dV ´ t ¨ η dA “ 0. (2.38) D D BDσ Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 29 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Gradient testne funkcije Gradη se lahko interpretira kot smerni odvod deformacijskega gradienta D F ¨ η, poznan kot variacija δF deformacijskega gradienta. V ena£bi (2.38) se prvi Piola-Kirchoov napetostni tenzor ob upo²tevanju povezave P “ F S lahko zamenja z drugim Piola-Kirchoovim napetostnim tenzorjem, kar vodi v P 1 ¨ Grad η “ S ¨ FT Gradη “ S ¨ pFT Grad η ` GradT η Fq “ S ¨ δE, (2.39) 2 pri £emer je bilo upo²tevano, da je skalarni produkt simetri£nega tenzorja S z antisi- metri£nim delom tenzorja enak ni£. δE predstavlja variacijo Green-Lagrangevega defor- macijskega tenzorja. Ob upo²tevanju ena£be (2.39), se ena£bo (2.38) lahko preoblikuje v ż ż ż S ¨ δE dV ´ ρ ¯ ¯ 0b ¨ η dV ´ t ¨ η dA “ 0. (2.40) D D BDσ Prvi £len v tej ena£bi predstavlja virtualno delo notranjih sil, preostala dva £lena pa virtualno delo zunanjih sil. Ena£ba (2.40) je primerna za vse probleme. V primeru da ima problem potencial, lahko uporabimo enostavnej²i pristop. Pristop se imenuje princip stacionarnega elasti£nega potenciala. Za hiperelasti£ni material obstaja funkcija W , ki denira elasti£no energijo, shranjeno v materialu. Pod predpostavko konservativne obteºbe (to pomeni, da je obteºba neodvisna od poti), lahko opi²emo stati£ni problem s funkcionalom ż ż Π ` pϕq “ W pCpϕqq ´ ρ ¯ ¯ 0b ¨ ϕ˘ dV ´ t ¨ ϕ dA ùñ STAT, (2.41) D BDσ kjer je C desni Cauchy-Greenov deformacijski tenzor. Izmed vseh moºnih deformacij ϕ je deformacija, ki zgotovi pogoje ravnoteºja, tista de- formacija, ki ima za posledico stacionarno vrednost potenciala Π. Stacionarna vrednost Π se lahko izra£una z variacijo Π glede na deformacije. V ta namen se uporabi smerni odvod d ˇ ˇ δΠ “ DΠpϕq ¨ η “ Πpϕ ` αηqˇ , (2.42) dα ˇα“0 kar vodi do ż ˆ ˙ BW ż DΠpϕq ¨ η “ ¨ η ´ ρ ¯ ¯ 0b ¨ η dV ´ t ¨ η dA “ 0. (2.43) D Bϕ BDσ Denicija mehanskega problema z uporabo principa stacionarnega elasti£nega potenciala ima ve£ prednosti. Med drugim vodi (ob uporabi avtomatskega odvajanja) do numeri£no najbolj u£inkovite formulacije kon£nih elementov. 30 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 2.2.2 Formulacija mehanskega problema z uporabo avtomatskega odvajanja Kot je bilo predstavljeno, je ²ibka oblika ravnoteºnih ena£b mehanskega problema oblike ş a ¨ δb dV ` ... “ 0, kjer sta a in b tenzorja poljubnega reda, δb pa je smerni odvod ali D variacija tenzorja b. V naslednjem koraku se ²ibko obliko ena£be diskretizira. Variacijo tenzorja b se izra£una kot Bbppq δbppq “ Dbppqδp “ δp, (2.44) Bp kjer je p “ tp1, p2, ..., pn uT mnoºica neznanih parametrov problema in ps δp “ tδp1, δp2, ..., δpn uT variacije neznanih parametrov problema. Skalarni produkt ps appq ¨ δbppq se izrazi kot ˆ ˙ Bbppq Bbppq appq ¨ δbppq “ appq ¨ δp “ appq ¨ δp. (2.45) Bp Bp Izrazimo produkt a Bbppq ppq ¨ po komponentah Bp ˆ ˙ ˆ ˙ Bbppq Bbppq a T Bbppq ppq ¨ “ appq ¨ “ tr appq , (2.46) Bp Bp Bp m m m kjer tr pomeni sled matrike. Iz tega sledi diskretizirana oblika ²ibke oblike ena£be me- hanskega problema n ż ps ˆż ˙ ÿ Bbppq appq ¨ δbppq dV ` ... “ appq ¨ dV δpm ` ... “ 0. (2.47) D Bp m m “1 D in se lahko formira v sistem nps algebrajskih ena£b ż Bbppq R “ appq ¨ dV ` ... “ 0. (2.48) D Bp V primeru, da je mehanski problem deniran kot minimum potenciala Π ş “ W ppq dV , D kjer je p “ tp1, p2, ..., pn uT mnoºica neznanih parametrov problema, se variacijo Πppq ps izra£una kot BΠppq ż δW ppq δΠppq “ δp “ dV δp “ 0. (2.49) Bp D Bp Ena£ba (2.49) se lahko formira v sistem nelinearnih algebrajskih ena£b oblike ż BW ppq R “ dV “ 0. (2.50) D Bp Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 31 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Rezultirajo£ sistem algebrajskih ena£b imenujemo primarni problem. Ena£be so lahko linearne ali nelinearne, odvisno od narave mehanskega problema. 2.2.3 Klasikacija mehanskih problemov in njihova formulacija v metodi kon£- nih elementov Mehanske probleme se najpogosteje klasicira glede na to ali so £asovno odvisni ali neod- visni ter glede na to, kako so diskretizirane ena£be formirane, v povezane ali nepovezane. V primeru najenostavnej²ih izmed mehanskih problemov, t.j. nepovezanih in £asovno neodvisnih problemov re²ujemo sistem algebrajskih ena£b Rppq “ 0. (2.51) Kot re£eno, se v metodi kon£nih elementov domeno problema razdeli na poddomene ali kon£ne elemente. Nato se na nivoju posameznega kon£nega elementa izra£una prispevek tega elementa Re k globalnemu rezidualu. Prispevke vseh kon£nih elementov se nato s proceduro zdruºevanja kon£nih elementov se²teje v globalno matriko konstrukcije R “ A Re, (2.52) e pri £emer je potrebno upo²tevati tudi kinemati£ne kompatibilnosti med elementi. Ker je rezidual R deniran z integracijo preko celotne domene problema, je tudi prispevek kon£nega elementa Re formiran z integracijo preko domene elementa. Za integriranje se najpogosteje uporabi numeri£na integracija in sicer Gaussova kvadratura. Gaussova kvadratura je numeri£na metoda integracije, pri kateri se integral izra£una tako, da se v dolo£enih to£kah g domene izra£una vrednosti funkcije in jih pomnoºi z ustreznimi uteºmi ter produkte se²teje Ng ÿ Re “ wgRg (2.53) g“1 Ena£be (2.51) so, odvisno od mehanskega problema, lahko linearne ali nelinearne. Ne- linearnost lahko izvira iz razli£nih vzrokov in sicer je nelinearnost lahko posledica npr. geometrije problema ali nelinearnih konstitutivnih ena£b. Metoda, ki je najpogosteje uporabljena za re²evanje nelinearnih sistemov algebrajskih ena£b in ki jo bomo upora- bljali tudi v disertaciji je Newton-Raphsonov iterativni postopek. Newton-Raphsonov postopek namre£ odlikuje kvadrati£na konvergenca v bliºini pravilne re²itve. Do tu je bil v grobem predstavljen postopek re²evanja najenostavnej²ega izmed razredov mehanskih problemov - £asovno neodvisnega in nepovezanega. V doktorski disertaciji bomo predstavili in izpeljali formulacijo perturbacijske metode za najzahtevnej²i razred 32 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. mehanskih problemov: £asovno odvisne povezane probleme, saj je aplikacija na ostale nato enostavna, pri £emer se iz ena£b le £rta odve£ne £lene. V skupino £asovno odvisnih povezanih problemov spadajo npr. elasto-plasti£ni problemi. Za primerjavo, kako izgleda perturbacijska metoda za enostavnej²e primere, bo na koncu podano tudi poglavje, v katerem bojo predstavljeni izrazi za £asovno neodvisne in nepovezane probleme, kot so npr. elasti£ni ali hiperelasti£ni problemi. V primeru £asovno odvisnih povezanih problemov se znotraj Newton-Raphsonove zanke (ki se izvaja na globalnem nivoju problema) uporabi ²e eno Newton-Raphsonovo zanko v kateri se iterativno linearizira ta dodatni sistem nelinearnih ena£b. Povezan problem torej sestavljata odvisen in neodvisen problem, ki sta izraºena z reziduali Qg “ 0, g “ 1, 2, ..., Ng (2.54) R “ 0, (2.55) kjer je R neodvisni rezidual, deniran in formiran na globalnem nivoju problema, Qg pa je g-ti odvisni rezidual, obi£ajno deniran v Gaussovi integracijski to£ki g in re²en na nivoju kon£nega elementa z vgnezdeno Newton-Raphsonovo iteracijsko zanko. Re²itve ena£b (2.54) so odvisni vektorji re²itve hg, re²itev ena£be (2.55) pa neodvisni vektor re²itve p. Pri £asovno odvisnih povezanih problemih, se £asovna domena (psevdo ali realnega £asa) razdeli v nkor £asovnih korakov. Vsakemu £asovnemu koraku pripada neodvisni vektor re²itve p0, p1,...,pn in odvisni vektorji re²itve h ´1, pn,pn`1,..., pnkor g,0, hg,1,...,hg,n´1, hg,n, hg,n , g `1,..., hg,n “ 1, 2, ..., N kor g . V splo²nem je re²itev v trenutnem £asovnem koraku odvisna od re²itev v vseh prej²njih korakih, vendar so v procesu avtomatizacije pomembni le tisti vektorji re²itve, ki se v formulaciji kon£nih elementov pojavijo eksplicitno. V doktorski disertaciji bomo tako privzeli najenostavnej²o moºnost, pri kateri so £asovni odvodi aproksimirani s kon£nimi diferencami prvega reda. V tem primeru je problem eksplicitno odvisen le od vektorjev re²itve v prej²njem (pn in hg,n) in trenutnem £asovnem koraku (pn`1 in hg,n`1, ki jih bomo, z namenom bolj²e preglednosti, pisali brez indeksov: p in hg). V primeru, da je re²itev odvisna od ve£ predhodnih korakov, ostaja proces avtomatizacije enak, le da je v posamezne ena£be potrebno dodati ²e vektorje re²itve iz predhodnih korakov, od katerih je re²itev odvisna. Na podlagi privzete formulacije se povezan £asovno odvisen problem lahko zapi²e kot R pp, h, pn, hnq “ 0, Qg ppe, hg, pe,n, hg,nq “ 0, g “ 1, 2, ..., Ng. (2.56) Za izra£un prispevka Re je potrebno poznati tudi vektorje re²itve povezanega problema hg in hg,n. Zato je v vsaki Newton-Raphsonovi iteraciji potrebno re²iti najprej sisteme ena£b Qg “ 0, ki se jih denira in re²i v vsaki Gaussovi integracijski to£ki kon£nega Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 33 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. elementa. Odvisen rezidual Qg je eksplicitno odvisen le od hg, hg,n in podmnoºic pe in pe,n globalnega neodvisnega vektorja re²itve p. V primeru elasto-plasti£nih problemov je hg vektor spremenljivk stanja za g-to integracijsko to£ko in vsebuje komponente plasti£- nega deformacijskega tenzorja in plasti£nega mnoºitelja, ena£be Qg “ 0 pa predstavljajo nelinearni sistem evolucijskih ena£b. Vektor p predstavlja prostostne stopnje problema, ki so obi£ajno pomiki in zasuki v vozli²£ih kon£nih elementov. Rezultirajo£ povezani sis- tem ena£b (2.56) predstavlja ena£be primarnega problema in se obi£ajno re²i z uporabo Newton-Raphsonovega iterativnega postopka, ki je prikazan v levem okvirju na sliki 2.4. Ko vektorji re²itve p in hg skonvergirajo, kar pomeni da je norma razlike re²itve v trenutni iteraciji od re²itve v prej²nji iteraciji manj²a od izbrane natan£nosti, se primarna analiza zaklju£i. Primarna analiza Občutljivostna analiza prvega reda 0 p : p n Dp Dp K : i I   R  , i  1, 2,..., M D D i i 0 h : h g gn Dh Q g  Dp i 1 g e   : Z  A , j   j h :  A Q g g D p D g  g  1  j g i e i  j 1  j  j h : h  h i  1, 2,..., M , g  1, 2,..., Ng g g g ?  j Ne h   Občutljivostna analiza drugega reda g j : j  1 2 2 g := g + 1 Ja D p D p K : ij II   R  , i, j  1, 2,..., M  D D D D i  i  i K p : R i j i j  i  1  i  i p  p   2 : p 2 h  i := i + 1 D Q g  D p ij 1 : g e  Z  A , g g D D p D D i j e i j ? Ne  i p   Ja i, j  1, 2,..., M , g  1, 2,..., Ng ... n := n + 1 Naslednji časovni korak Občutljivostna analiza n-tega reda Slika 2.4: Procedura primarne analize in ob£utljivostne analize vi²jega reda. Figure 2.4: Procedure of primal analysis and higher-order sensitivity analysis. 2.3 Perturbacijska metoda Perturbacijska metoda za izra£un statistike odziva konstrukcije v mehaniki konstrukcij je direktno povezana z matemati£no formulacijo perturbacije. Bistvo perturbacijske metode 34 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. je, da se izhaja iz enostavnej²ega (v stohasti£nem pristopu to pomeni deterministi£nega oz. primarnega) sistema, za katerega je re²itev enostavno izra£unljiva. Nato se v ta eno- stavnej²i sistem vnese motnjo (perturbacijo). Na ta na£in se zikalne koli£ine (npr. vektor re²itve primarnega problema p, ki v mehaniki trdnin obi£ajno predstavlja prostostne sto- pnje - pomike in zasuke vozli²£ kon£nih elementov) ob predpostavki zveznosti sistema lahko izrazijo kot popravki enostavnega sistema. V perturbacijski metodi se popravke izrazi z metodo aproksimacije z razvojem odziva v Taylorjevo vrsto. Pri perturbacijski metodi nas torej ne zanima le deterministi£na re²itev vektorjev p in hg, temve£ i²£emo razvoj vektorjev re²itve v Taylorjevo vrsto okrog pri£akovanih vrednosti spremenljivk, katerih vpliv na odziv sistema nas zanima. V metodi stohasti£nih kon£nih elementov so to slu£ajne spremenljivke ξ. s1 s2 sM ÿ ÿ ÿ 1 ppξ1, ξ2, ..., ξM q “ . . . ¨ l l 1!l2! . . . lM ! 1“0 l2“0 lM “0 ˇ Bl1`l2`...`lM p ˇ ` 0 ˘l1 ` 0 ˘l2 0 ˘lM ¨ ˇ ξ ξ ξ ξ . . . `ξ ξ , 0 1 ´ 1 2 ´ 2 M ´ M Bξl1 ... ˇ ξi “ ξi, 1 Bξl2 2 BξlM M i “ 1, ..., M (2.57) s1 s2 sM ÿ ÿ ÿ 1 hgpξ1, ξ2, ..., ξM q “ . . . ¨ l l 1!l2! . . . lM ! 1“0 l2“0 lM “0 ˇ Bl1`l2`...`lM hg ˇ ` 0 ˘l1 ` 0 ˘l2 0 ˘lM ¨ ˇ ξ ξ ξ ξ . . . `ξ ξ , 0 1 ´ 1 2 ´ 2 M ´ M Bξl1 ... ˇ ξi “ ξi, 1 Bξl2 2 BξlM M i “ 1, ..., M (2.58) kjer je si red razvoja v vrsto za i-ti parameter in M ²tevilo parametrov, za katere nas zanima odziv konstrukcije oz. sistema v odvisnosti od tega parametra. V okviru me- tode stohasti£nih kon£nih elementov predstavlja M ²tevilo £lenov K-L dekompozicije. V primeru Gaussovega stohasti£nega polja, ki smo ga privzeli v disertaciji, je 0ξ “ 0 in bi se zato iz zgornjih izrazov lahko odstranili, vendar smo jih kljub temu pustili v iz- razih z namenom, da je postopek primeren tudi za druge verjetnostne porazdelitve ali nestohasti£ne spremenljivke. Ena£be (2.57) in (2.58) namre£ veljajo tudi za splo²en, deterministi£en mehanski problem - v tem primeru so ξ1, ξ2, ..., ξM poljubni parametri, za katere nas zanima vpliv spremembe teh parametrov na odziv konstrukcije. Zaradi bolj²e Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 35 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. preglednosti, zapi²imo izraza (2.57) in (2.58) v bolj splo²ni obliki s1 s2 sM ÿ ÿ ÿ p ` 0 ˘l1 ` 0 ˘l2 0 ˘lM pξ1, ξ2, ..., ξM q “ . . . pl ξ ξ ξ ξ . . . `ξ ξ , 1l2...lM 1 ´ 1 2 ´ 2 M ´ M l1“0 l2“0 lM “0 (2.59) s1 s2 sM ÿ ÿ ÿ h ` 0 ˘l1 ` 0 ˘l2 0 ˘lM g pξ1, ξ2, ..., ξM q “ . . . hg,l ξ ξ ξ ξ . . . `ξ ξ , 1l2...lM 1 ´ 1 2 ´ 2 M ´ M l1“0 l2“0 lM “0 (2.60) kjer so pl in h koecienti razvoja v Taylorjevo vrsto. Priporo£ilo o izbiri 1l2...lM g,l1l2...lM ustreznega reda, do katerega se vrsto razvije, je odvisno predvsem od nelinearnosti pro- blema. V splo²nem velja, vi²ji je red Taylorjeve vrste, ve£je je obmo£je okolice odziva, na katerem je odziv dovolj natan£no opisan. Zato je red, do katerega bi bilo potrebno raz- viti Taylorjevo vrsto za izpolnitev ºelene natan£nosti, odvisen od narave in nelinearnosti posameznega problema. Bolj kot je problem nelinearen, ve£ £lenov vrste je potrebnih. Prednosti razvoja odziva v vrsto, izraºeno kot funkcijo slu£ajnih spremenljivk, sta pred- vsem dve. Prva je, da je metoda, ne glede na uporabljeno tehniko izra£una koecientov vrste, znatno cenej²a v primerjavi z metodo Monte Carlo. Druga pomembna prednost pa je, da je rezultat polinom, pri katerem so slu£ajne spremenljivke izraºene simbolno, kar omogo£a enostavno analizo narave polinoma. Glavna slabost perturbacijske metode je, da je primerna le za zvezen odziv sistema in ne zmore zajeti diskretnih dogodkov. Kar pomeni, da £e je npr. v primeru uklonske analize rezultat deterministi£nega sistema globalni uklon konstrukcije, bo perturbacijska metoda opisala odziv konstrukcije za spre- membe slu£ajnih spremenljivk v smislu globalnega uklona, v resnici pa bi pri dolo£enih vrednostih slu£ajnih spremenljivk lahko pri²lo prej do lokalnega uklona kot globalnega. Druga slabost je, da nam tudi ob predpostavki zveznega odziva sistema da odgovor o obna²anju odziva dovolj natan£no le na dolo£enem obmo£ju v okolici to£ke, okrog katere smo razvili Taylorjevo vrsto, medtem ko se iz dobljene vrste ne da z gotovostjo oceniti odziva za vrednosti stohasti£nih spremenljivk, ki so dale£ od pri£akovanih. Vendar, kljub navedenim slabostim metode, ostaja perturbacijska metoda ena bolj priljubljenih in po- gosto uporabljenih metod. Ob upo²tevanju njenih omejitev so rezultati te metode dokaj natan£ni (pri£akovane vrednosti in variance odziva) in numeri£no u£inkoviti. V naslednjih poglavjih so predstavljene tri moºnosti za izra£un koecientov pl : sim- 1l2...lM bolna ob£utljivostna analiza, analiti£na ob£utljivostna analiza in metoda kon£nih diferenc. Simbolna ob£utljivostna analiza spada v kategorijo semi-analiti£nih re²itev primarnega problema (Iokamidis, 1992, Iokamidis, 1993, Zhiming in Jianming, 1998, Pavlovi¢, 2003, Korelc, 2011). Ta metoda uporabi programsko okolje, ki omogo£a simbolno ra£unanje. Na ta na£in se razvoj v vrsto izra£una med samim re²evanjem problema po metodi kon£nih elementov in sicer tekom izvajanja programa znotraj Newton-Raphsonovih iteracij. Zato 36 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. je metoda lahko direktno aplicirana za poljuben red vrste. Analiti£na ob£utljivostna analiza ra£una odvode Bl1`l2`...`lM p in Bl1`l2`...`lM hg (ena£be l l l l l l Bξ 1 2 ... M 1 2 ... M 1 Bξ2 Bξ Bξ Bξ M 1 Bξ2 M (2.57) in (2.58)). Dobljene odvode je nato za izra£un £lenov vrste potrebno pomnoºiti z 1 . Pri tej metodi je, za razliko od simbolne ob£utljivostne analize, kot je prikazano l1!l2!...lM ! v naslednjem poglavju, potrebno predhodno izpeljati ustrezne izraze in jih vnesti v kodo za kon£ne elemente. Metodo odlikuje visoka numeri£na u£inkovitost. Metoda kon£nih diferenc je numeri£na metoda, ki izra£una odvode Bl1`l2`...`lM p in l l l Bξ 1 2 ... M 1 Bξ2 BξM Bl1`l2`...`lM hg in je popolnoma splo²na, vendar nestabilna in zato nezanesljiva metoda. l l l Bξ 1 2 ... M 1 Bξ2 BξM Njena slabost je tudi numeri£na neu£inkovitost. 2.3.1 Analiti£na ob£utljivostna analiza Ob£utljivostna analiza i²£e odvode (poljubnega reda) vektorjev re²itve po izbranih ob- £utljivostnih parametrih ξ “ tξiu, i “ 1, 2, ..., M. Za iskane odvode se uporablja tudi izraz ob£utljivosti. V splo²nem je lahko predmet prou£evanja ob£utljivosti lahko poljuben funkcional odziva, npr. napetosti, deformacije, prostornina, pomiki ipd. Ker so funk- cionali odziva odvisni od vektorjev re²itve p in hg, je zato potrebno najprej izra£unati ob£utljivosti vektorjev re²itve. Za izra£un odvodov vektorjev re²itve je potrebno odvajati osnovne ena£be (2.56) po Bl1`l2`...`lM . Ker ºelimo, da je avtomatizacija ob£utljivostne analize uporabna za po- l l l Bξ 1 2 ... M 1 Bξ2 BξM ljubne primere, privzamemo da ob£utljivostni parametri v ena£bah rezidualov lahko na- stopajo eksplicitno in/ali implicitno, zato osnovne rezidualne ena£be £asovno odvisnega povezanega problema zapi²emo kot R pp pξq , h pξq , pn pξq , hn pξq , ξq “ 0, Qg ppe pξq , hg pξq , pe,n pξq , hg,n pξq , ξq “ 0, g “ 1, 2, ..., Ng. (2.61) Ena£ne nato odvajamo po enem ali ve£ parametrih ξi, ξj,... Po preureditvi dobimo linearen sistem ena£b za izra£un ob£utljivosti vektorja p, ki je analogen lineariziranemu sistemu ena£b za primarni problem. Sledi izra£un ob£utljivosti vektorjev hg, ki se izvede na nivoju posameznega kon£nega elementa. Celoten postopek ob£utljivostne analize vi²jega reda in pristop z avtomatskim odvajanjem, razvit v okviru doktorske disertacije, je predstavljen v poglavju 6. 2.3.2 Simbolna ob£utljivostna analiza Pri simbolni ob£utljivostni analizi se spremenljivke, okrog katerih ºelimo razviti odziv konstrukcije v Taylorjevo vrsto, pusti v simbolni obliki. Zaradi tega je ta na£in re²evanja Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 37 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. stohasti£nih problemov moºen le ob uporabi programskih okolij, ki omogo£ajo simbolno ra£unanje, to je npr. Mathematica. Prednost uporabe simbolnega sistema je, da je za uporabnika vloºek, ki ga mora v primerjavi z opisom problema za primarno analizo dodati v opis problema, da dobi razvoj re²itve v Taylorjevo vrsto, zanemarljiv, dobljena re²itev pa je analiti£na. Glavna slabost simbolne ob£utljivostne analize je, da je uporaba simbol- nega sistema £asovno zelo potratna in obi£ajni mehanski problemi zelo hitro prerastejo zahtevnost, ki bi bila ²e sprejemljiva za izra£un problema v doglednem £asu. Bistvena razlika med analiti£no in simbolno ob£utljivostno analizo je, da se pri prvi odvaja le rezidual za konvergirano re²itev, medtem ko se pri slednji odvaja vse korake Newton- Raphsonove iterativne zanke in je zato potrebno razviti v vrsto tudi tangentno matriko K. Zato se bomo pri aplikaciji simbolne ob£utljivostne analize omejili na najenostav- nej²e, numeri£no najcenej²e mehanske probleme: £asovno-neodvisne in nepovezane. Upo- raba simbolne ob£utljivostne analize nam bo sluºila predvsem za kontrolo izpeljane ADB formulirane ob£utljivostne analize. Pri £asovno neodvisnih in nepovezanih problemih v rezidualni ena£bi nastopa le iskani vektor re²itve p R ppq “ 0. (2.62) V splo²nem je £asovno-neodvisni nepovezani problem lahko nelinearen, bodisi zaradi ge- ometrije problema ali zaradi materialnega modela. Zato se za re²evanje lahko uporabi Newton-Raphsonova metoda. Pri simbolni ob£utljivostni analizi se vektor re²itve v i-ti Newton-Raphsonovi iteraciji ppiq, tangentni operator Kpiq ter rezidual Rpiq izrazi v odvi- snosti od spremenljivk ξ Kpiqpξq ∆ppiqpξq ` Rpiqpξq “ 0 (2.63) BRpiqpξq Kpiqpξq “ (2.64) Bppiqpξq ppì1qpξq “ ppiqpξq ` ∆ppiqpξq, (2.65) kjer je ∆ppiq inkrement vektorja re²itve p v i-ti Newton-Raphsonovi iteraciji. Ker se pri simbolni ob£utljivostni analizi odvaja celotno Newton-Raphsonovo zanko, je potrebno 38 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. razviti v vrsto vse koli£ine, ki v Newton-Raphsonovi zanki nastopajo s1 s2 sM ÿ ÿ ÿ Kpiq ` 0 ˘l1 ` 0 ˘l2 0 ˘lM pξ1, ξ2, ..., ξM q “ . . . Kpiq ξ ξ ξ ξ . . . `ξ ξ , l 1 ´ 1 2 ´ 2 M ´ M 1l2...lM l1“0 l2“0 lM “0 s1 s2 sM ÿ ÿ ÿ Rpiq ` 0 ˘l1 ` 0 ˘l2 0 ˘lM pξ1, ξ2, ..., ξM q “ . . . Rpiq ξ ξ ξ ξ . . . `ξ ξ , l 1 ´ 1 2 ´ 2 M ´ M 1l2...lM l1“0 l2“0 lM “0 s1 s2 sM ÿ ÿ ÿ ∆ppiq ` 0 ˘l1 ` 0 ˘l2 0 ˘lM pξ1, ξ2, ..., ξM q “ . . . ∆ppiq ξ ξ ξ ξ . . . `ξ ξ . l 1 ´ 1 2 ´ 2 M ´ M 1l2...lM l1“0 l2“0 lM “0 (2.66) Newton-Raphsonove iteracije se izvajajo dokler ni doseºena konvergenca vseh £lenov ∆ppiq . l1l2...lM › › ›∆ppiq , l › ă ε › l k “ 1, ..., sk 2. (2.67) 1l2...lM › V Mathematici sta Kpiq g in Rpiq g izraºena z vrsto avtomati£no, £im je kateri od vhodnih parametrov izraºen z vrsto. Koecienti Kpiq in Rpiq so izra£unani v Gaussovih g,l1l2...lM g,l1l2...lM to£kah in se po pravilih Gaussove kvadrature pomnoºijo z ustreznimi uteºmi ter se²tejejo v prispevek kon£nega elementa Kpiq in Rpiq . Prispevki posameznih kon£nih e,l1l2...lM e,l1l2...lM elementov se nato s standardno proceduro formiranja se²tejejo v globalne matrike in vektorje koecientov Kpiq and Rpiq . Ker sta Kpiq in Rpiq v obliki vrste, je tudi l1l2...lM l1l2...lM izra£unani ∆ppiq v obliki vrste. Zaradi uporabe simbolnega algebrajskega sistema za ra£unanje z vrstami je za izvedbo simbolna ob£utljivostna analiza re²evanja linearnega sistema ena£b potrebno veliko ²tevilo operacij, zaradi £esar so ra£unski £asi zelo dolgi. Obstaja tudi druga£en pristop, ki ²tevilo operacij zmanj²a in s tem skraj²a ra£unski £as. Po tem postopku se koeciente ∆ppiq l1l2...lM izra£una z rekurzivno formulo l1 l2 lM ÿ ÿ ÿ Kpiq ∆ppiq ¨ ¨ ¨ Kpiq ∆ppiq , 00...0 l “ Ŕpiq ´ 1l2...lM l1l2...lM l1l2...lM pl1ŕ1qpl2ŕ2q...pl r M ŕM q 1“0 r2“0 rM “0 except:r1“r2“¨¨¨“rM “0 (2.68) lk “ 0, 1, ..., sk. Ob tem velja opozoriti, da imajo koecienti matrik Kpiq enako strukturo razpr²enosti l1l2...lM kot matrika Kpiq ter da je v tem primeru potrebno izra£unati le razcep (ki je draga numeri£na operacija) matrike Kpiq . Iz ena£be (2.68) je tudi razvidno, da koecienti vi²jih 00...0 redov ∆ppiq nimajo vpliva na koeciente niºjih redov. Na podlagi te ugotovitve, se l1l2...lM lahko numeri£no u£inkovitost ²e pove£a, kot je prikazano na sliki 2.5, pri £emer se Newton- Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 39 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Raphsonovo zanko razdeli na dva dela. V prvem delu se i²£e re²itve ni£elnih koecientov p00...0. Ta del Newton-Raphsonove zanke je popolnoma numeri£en in predstavljata re²itev deterministi£nega primarnega problema za pri£akovane vrednosti slu£ajnih parametrov 0ξ. Za re²itev tega dela se torej lahko uporabi povsem numeri£ne programe, kot je npr. CDriver (ki je napisan v programskem jeziku C). Po doseºeni konvergenci Newton-Raphsonovih iteracij za ni£ni £len p00...0 se pri£ne drugi del Newton-Raphsonove zanke. V tem delu se izvaja rekurzivno ra£unanje preostalih koecientov vrste, kot prikazuje slika (2.5). V tem delu je potrebno uporabiti simbolni algebrajski sistem, ki omogo£a simbolno operiranje z vrstami. V ta namen se lahko upo- rabi MDriver, ki je komplementaren modul CDriverju. MDriver je v celoti zapisan v simbolnem jeziku programskega okolja Mathematica in ima enak input, strukturo podat- kov in ukazni jezik kot CDriver. Avtomatski generator kode AceGen izvede avtomatsko generacijo podprogramov za kon£ne elemente za MDriver in CDriver iz istega simbolnega opisa kon£nega elementa (podrobnej²i opis AceGen-a in simbolno-numeri£nega pristopa je v poglavju 3.1). Ko je doseºena konvergenca vseh koecientov, se pri£ne izra£un za naslednji £asovni korak. Lastnost Newton-Raphsonove procedure, da je kvadrati£no konvergentna, velja tudi za koeciente razvoja v vrsto. To pomeni, da ko je enkrat doseºena konvergenca ni£elnega koecienta p00...0, se ²tevilo pravilno izra£unanih koecientov vrste z vsako iteracijo pod- ´ ´ ¯¯ voji. Kar pomeni, da je potrebno najve£ ceiling log 1 řM s iteracij. Funkcija 2 ì“1 i ceiling pxq “ min tm P Z|m ě xu izra£una najmanj²e celo ²tevilo ve£je ali enako vrednosti x. 2.3.3 Metoda kon£nih diferenc Tudi tretja predstavljena metoda, metoda kon£nih diferenc, sluºi v doktorski disertaciji predvsem za kontrolo pravilnosti rezultatov in primerjavo numeri£ne u£inkovitosti. Zato bo v tem poglavju predstavljena metoda kon£nih diferenc kolikor je mogo£e strnjeno. Osnovne ena£be metode kon£nih diferenc so izpeljane iz razvoja funkcije v Taylorjevo vrsto. Prvi in drugi odvodi funkcije ve£ spremenljivk f px, yq se po tej metodi izra£unajo Bf px, yq f px ` h, yq ´ f px ´ h, yq B2f px, yq f px ` h, yq ´ 2f px, yq ` f px ´ h, yq “ , “ , Bx 2h Bx2 h2 B2f px, yq f px ` h1, y ` h2q ´ f px ´ h1, y ` h2q ´ f px ` h1, y ´ h2q ` f px ´ h1, y ´ h2q “ . BxBy 4h1h2 (2.69) Najve£ja slabost metode kon£nih diferenc je, da metoda ni stabilna. Razlog temu je, da ima dva izvora napak: napaka metode zaradi zanemaritve vi²jih £lenov Taylorjeve vrste in zaokroºitvena napaka, ki je posledica ra£unalni²kega zapisa ²tevil na omejeno 40 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 1.del: CDriver, numerično reševanje 2.del: MDriver, simbolno reševanje n : 0  j  i j := 1 p   p : 0 : p 00...0 00...0 00...0,0  j p : , 0 l  1, 2,..., s l l ... l k k 1 2 M i : 0 ( i ) p : p 00...0 00...0, n  j  j  j K p : R  00...0 l l ... l l l ... l 1 2 M 1 2 M l l l 1 2 M   j  j i  i  i    K  K p : R ... p r r ... r    1 2 M  l r l r ... l r 1 1  2 2   M M  00...0 00...0 00...0 r 0 r 0 r 0  1 2 M i   1  i  i :1 p : p  p except r  r ... r 0 1 2 M 00...0 00...0 00...0 ii  j  1  j  j   p : p  p :1 :1 l l ... l l l ... l l l ... l jj 1 2 M 1 2 M 1 2 M nn ? l : 0,1,..., s Ne  k k i p   00...0 Ja ? M ( 1  )    Ja p : i  p 00...0, n 1  00...0 j  ceiling  log 1  s   2 i   i 1   Ne Ja n  nkor Ne Slika 2.5: Newton-Raphsonova procedura re²evanja, prilagojena za re²evanje razvoja vek- torja re²itve v vrsto za £asovno neodvisne mehanske probleme. Figure 2.5: NewtonRaphson solution procedure adjusted for power series solution to the system of nonlinear equations. natan£nost v t.i. formatu plavajo£e vejice. Medtem ko napaka zaradi zanemaritve £lenov pada z manj²anjem h, se zaokroºitvena napaka z manj²anjem h pove£uje. Posledi£no je ustrezen h potrebno izbrati zelo previdno, da se skupna napaka minimizira. Vendar pa izbira optimalnega h, ²e posebej pri kompleksnej²ih funkcijah, ni tako preprosta in jo je nemogo£e dolo£iti vnaprej. Zato je obi£ajno potrebnih nekaj iteracij, da se oceni ustrezen h in doseºe ºeleno natan£nost. ƒe ob tem upo²tevamo, da v metodi kon£nih elementov vsak klic funkcije f px, yq v izrazih (2.69) pomeni izvedbo celotne primarne analize, to pomeni, da je v vsaki iteraciji potrebno izvesti 2 do 4 primarne analize ter to ²tevilo pomnoºiti s ²tevilom iteracij, potrebnih da se oceni ustrezen h. Hitro postane jasno, da je metoda, predvsem v primeru kompleksnih mehanskih problemov, zelo draga. Za pocenitev ra£unskega postopka metode kon£nih diferenc v disertaciji uporabljamo iz- bolj²ano metodo kon£nih diferenc. V tej metodi se uporabi Riddersova implementacija Richardsonove ekstrapolacije h Ñ 0 ter Nevillovega algoritma, po katerem je vsak izra£un kon£ne diference uporabljen hkrati za ekstrapolacijo vi²jega reda in ekstrapolacijo pred- hodno izra£unanih kon£nih diferenc, vendar tokrat z manj²im h. h se v vsaki iteraciji manj²a. Z uporabo ekstrapolacije, ki je neprimerljivo numeri£no u£inkovitej²a od celo- tne primarne analize, se ²tevilo potrebnih iteracij in posledi£no numeri£na u£inkovitost Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 41 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. postopka praviloma izbolj²a. Podrobnosti algoritma so podane v Press et al. (1988-92). 2.4 Metoda Monte Carlo za izra£un statistike odziva konstrukcije Metoda Monte Carlo je metoda, pri kateri se izvede ve£je (poljubno) ²tevilo simulacij. V vsaki simulaciji se skladno z verjetnostno porazdelitvijo naklju£nih spremenljivk gene- rirajo naklju£ne vrednosti teh spremenljivk. Nato se z uporabo teh vrednosti re²i (de- terministi£ni) primarni problem in beleºi ºeleni funkcional odziva (npr. pomik, plasti£ne deformacije, porabljeno energijo itd.). Iz mnoºice re²itev se nato s standardnimi statisti£- nimi metodami dolo£i statistiko odziva. Prednost metode Monte Carlo je v tem, da je zelo splo²na in uporabna ne glede na naravo problema. Tudi natan£nost vseh novo razvitih metod za izra£un verjetnosti se zaradi tega obi£ajno preverja prav z metodo Monte Carlo. Slabost metode je, da rezultat, ne glede na to, koliko simulacij izvedemo, vedno vsebuje nekaj negotovosti, saj zaradi naklju£ne izbire parametrov ne moremo vedeti ali je bila izbrana tipi£na mnoºica moºnih vrednosti slu£ajnih parametrov. V splo²nem velja groba ocena, da je natan£nost metode Monte Carlo za ve£je ²tevilo simulacij proporcionalna ? 1{ n, kjer je n ²tevilo simulacij. Za generiranje naklju£nih ²tevil se obi£ajno uporablja psevdo-naklju£en generator na- klju£nih ²tevil (ang. pseudorandom number generator). Ve£ina programskih jezikov ima generatorje naklju£nih spremenljivk ºe vgrajene. Vrednosti, ki jih dobimo na ta na£in niso povsem naklju£ne, temve£ se po dolo£enem algoritmu iz relativno majhne mnoºice za£etnih vrednosti generira zaporedje ²tevil, ki aproksimirajo lastnosti slu£ajnih spremen- ljivk. Obstajajo tudi algoritmi, pri katerih so generirane spremenljivke ²e bliºje slu£ajnim, vendar so numeri£no bolj zahtevni. Za oceno ustreznosti generatorjev naklju£nih ²tevil obstajajo posebni testi. Eni najbolj znanih so naprimer t.i. diehard testi, ki jih je razvil G. Marsaglia, 1995. V splo²nem so za simulacije Monte Carlo dovolj natan£ni tudi generatorji psevdo-naklju£- nih ²tevil, ki so numeri£no u£inkovitej²i in ²e vedno dovolj dober pribliºek povsem na- klju£nim spremenljivkam. V disertaciji so naklju£ne spremenljivke za simulacije Monte Carlo generirane po metodi extended cellular automata, ki je ena od vgrajenih metod za generiranje psevdo-naklju£nih ²tevil v Mathematici in jo odlikuje visoka kvaliteta ge- neriranih ²tevil. Po tej metodi se generira vektor stanja dolºine 5120, ki ga skladno z deterministi£nim pravilom sestavljajo ni£le in enice. Nov vektor stanja se nato dolo£i na podlagi vrednosti dolo£enih sosednjih celic v prej²njem vektorju stanja. Podmnoºice celic predstavljajo naklju£ne bite, iz katerih se generirajo psevdonaklju£na ²tevila. Na ta na£in obstaja dolo£ena korelacija med trenutno in petimi predhodnimi celicami. V praksi se izkaºe, da se ob sestavljanju psevdonaklju£nih ²tevil na na£in, da se v vektorju stanja preskakuje po 4 celice (kar smo uporabili tudi pri generiranju naklju£nih ²tevil v diserta- ciji), doseºe visoko stopnjo kvalitete naklju£nosti generiranih ²tevil, ki izpolnijo tudi zelo 42 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. stroge teste naklju£nosti. Za metodo Monte Carlo v splo²nem velja: ve£je kot je ²tevilo simulacij, ve£ja bo natan£- nost izra£unane statistike odziva. Kar pomeni, da je metoda numeri£no precej draga, saj n simulacij Monte Carlo pomeni izvedbo n primarnih analiz. Dobra stran pri tem je, da so posamezne simulacije med seboj neodvisne in se ra£unski £as zato lahko skraj²a s paralelnim izvajanjem analiz na ve£ procesorjih. Ob tem je potrebno biti pozoren le na to, kak²en algoritem uporablja program za generiranje naklju£nih ²tevil oz. za generiranje semena, iz katerih se pri£ne generiranje psevdo naklju£nih ²tevil. V primeru, da bi se pri paralelnih izra£unih osnovnih ²tevil, uporabilo isto seme, bi se namre£ v teh izra£unih generirala identi£no ista ²tevila, kar bi posledi£no pomenilo korelacijo med generiranimi ²tevili in rezultat ne bi imel ºelene naklju£nosti. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 43 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 3 Hibridni simbolno-numeri£ni pristop Kot pove ime, simbolno-numeri£ni pristop predstavlja kombinacijo razli£nih pristopov in metod k re²evanju numeri£nih problemov. Cilj kombiniranja razli£nih metod je, da se uporabi prednosti vsake izmed kombiniranih metod, medtem ko se na podro£ju, na katerem je posamezna metoda ²ibka, aplicira drugo metodo. Kon£ni rezultat je zato bolj²i kot £e bi uporabili le eno izmed njih. V uporabljenem simbolno-numeri£nem pristopu v doktorski disertaciji smo kombinirali • program za simbolno ra£unanje Mathematica, • avtomatski generator programske kode z vgrajeno tehniko avtomatskega odvajanja AceGen in • okolje za numeri£no modeliranje po metodi kon£nih elementov AceFEM. Program za simbolno ra£unanje (ang. computer algebra system) omogo£a operiranje z izrazi na simbolni ravni. V disertaciji je izbran program za simbolno ra£unanje Mathe- matica (Mathematica 9.0). Glavna prednost programov za simbolno ra£unanje je, da omogo£ajo simbolno izpeljavo in poenostavitev izrazov. Uporaba programa za simbolno ra£unanje v metodi kon£nih elementov tako omogo£a zapis ena£b in algoritma kon£nega elementa v zelo strnjeni in splo²ni obliki. Prednost tako kratkega in preglednega zapisa ena£b je tudi v tem, da je iskanje morebitnih napak zato veliko laºje. Slabost programov za simbolno ra£unanje je, da so zaradi simbolnega zapisa funkcij izra£uni v teh programih zelo zahtevni, velikost izrazov pa eksponentno nara²£a z nara²£ajo£im ²tevilom matema- ti£nih operacij (kot so npr. odvajanje, razli£ne transformacije) in ²tevilom prostostnih stopenj. Zaradi tega so ra£unski £asi analiz lahko zelo dolgi in je izvedljivost analize mehanskih problemov na simbolni ravni zato omejena le na matemati£no enostavnej²e probleme z majhnim ²tevilom kon£nih elementov. Programi za simbolno ra£unanje sicer omogo£ajo tudi numeri£ne operacije, vendar so kljub temu numeri£no manj u£inkoviti v primerjavi s programskimi jeziki kot sta na primer FORTRAN ali C. Na tem mestu se po- kaºe prednost hibridnega pristopa, ki s prevedbo simbolnega zapisa v poljubni programski jezik omogo£a, da se izkoristi prednosti obeh pristopov - simbolnega in numeri£nega. Povezavo med simbolno izpeljanimi ena£bami in prevedbo ena£b v izbrani programski jezik omogo£a avtomatski generator kode kon£nih elementov AceGen (Korelc, 2009b). AceGen je programski paket, zdruºljiv s programom Mathematica. Uporabnik ena£be zapi²e v simbolni obliki. AceGen ena£be nato prevede v izbran programski jezik. Tu se 44 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. pokaºe ²e ena prednost AceGena, ki omogo£a prevedbo simbolnega zapisa ena£b v mnoge programske jezike (C, C++, FORTRAN, Mathematica, Matlab), ne da bi bilo potrebno v izvornem simbolnem zapisu karkoli spreminjati. Pri tem sta za numeri£no u£inkovitost izpeljane kode pomembni predvsem dve tehniki, ki sta vgrajeni v AceGen: avtomatsko odvajanje in simultana stohasti£na poenostavitev numeri£nih kod. Avtomatsko odvajanje je numeri£no zelo u£inkovita tehnika odvajanja, ki temelji na uporabi veriºnega pravila in je podrobneje predstavljena v poglavju 3.2. Simultana stohasti£na poenostavitev nume- ri£nih kod je alternativa obi£ajnemu pristopu optimiranja izrazov, po katerem se izraze optimira, tik preden se jih prevede v numeri£ni jezik in sicer tako, da se i²£e morebitne enake £lene ali dele £lenov v izrazih. Tak pristop se pri splo²nem nelinearnem mehanskem problemu namre£ izkaºe kot neu£inkovit. S tehniko simultane stohasti£ne poenostavitve izrazov se enakost izrazov i²£e z izvrednote- njem le-teh z uporabo naklju£nih vrednosti posameznih spremenljivk med samo izpeljavo izrazov (Korelc, 1997). Delom izrazov, ki se v razli£nih izrazih ponavljajo, se nato av- tomatsko dodeli nove pomoºne spremenljivke, kot je prikazano na sliki 3.1. Izkaºe se, da je simultana stohasti£na poenostavitev numeri£nih kod u£inkovita in primerna tudi za bolj kompleksne, nelinearne probleme. Po kon£ani optimizaciji izraza ima uporabnik omogo£eno svobodno manipuliranje z optimiranim izrazom, ki ga izpelje AceGen. Re- zultat vseh omenjenih tehnik, vgrajenih v AceGen, je uspe²no manj²anje problema rasti programske kode in izbolj²anje numeri£ne u£inkovitosti izpeljane in prevedene program- ske kode. Shema uporabljenega hibridno simbolno-numeri£nega pristopa k avtomatizaciji stohasti£ne metode kon£nih elementov je predstavljena na sliki 3.2. Podoben pristop, ki je v praksi pogosto uporabljen, je hibridni objektno-orientiran pristop (Eyheramendy in Zimmermann, 1999, Beall in Shephard, 1999, Logg, 2007). 3.1 Uporaba AceGen-a v metodi kon£nih elementov. V splo²nem so inºenirski problemi opisani s skupkom diferencialnih ali integralnih ena£b. Ob predpostavki poljubnih domen je direktna analiti£na re²itev le-teh celo ob enostav- nej²ih konstitutivnih zakonih materialov neizra£unljiva. Zato se za re²evanje ena£b, ki opisujejo inºenirske probleme, obi£ajno uporablja numeri£ne metode. Najpogosteje upo- rabljena metoda v re²evanju inºenirskih problemov je metoda kon£nih elementov. Metoda kon£nih elementov je variacijska metoda, po kateri se zvezni problem diskretizira in nato z ustreznimi matemati£nimi postopki pretvori v sistem nelinearnih algebrajskih ena£b, ki se jih lahko nato re²i s standardnimi metodami, npr. Newton-Raphsonovo metodo. Opis metode kon£nih elementov za analizo mehanskih problemov je bil podan v poglavju 2.2. AceGen je zasnovan tako, da je primeren tako za generiranje poljubnih programskih kod, deniranih s strani uporabnika, kot za generiranje kode za uporabo v metodi kon£nih Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 45 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Originalna matrika: Poenostavljena matrika, izražena z novimi 1 2 E I 6 E I 1 2 E I 6 E I    pomožnimi 3 2 3 2 L L L L spremenljivkami: 6 E I 4 E I 6 E I 2 E I v v  v v  4 5 4 5 2 2 L L L L AceGen v v  v v  5 6 5 7 K  0 K0 1 2 E I 6 E I 1 2 E I 6 E I     v v v v 4 5 4 5 3 2 3 2 L L L L v v  v v 5 7 5 6 6 E I 2 E I 6 E I 4 E I  2 2 L L L L Vektor s 7 novimi pomožnimi spremenljivkami: 1 2 v v 6 v v 4 v v v 1 2 1 2 1 2 6 v  { E , I , L , ,  , , } 3 2 v v v 2 3 3 3 Slika 3.1: Simultana stohasti£na poenostavitev izrazov, vgrajena v AceGen. Figure 3.1: Simultaneous optimization of the expressions that is built in AceGen. elementov. V primeru, da se kodo generira za uporabo v metodi kon£nih elementov, AceGen izpeljano kodo kon£nega elementa prilagodi zahtevam okolja za kon£ne elemente, v katerem bo generirana koda kon£nega elementa uporabljena (npr. AceFEM (Korelc, 2009c), ABAQUS, ELFEN, FEAP). Ob tem ni potrebno v izvorni simbolni kodi kon£nega elementa ni£esar spreminjati. Za re²evanje problemov z metodo kon£nih elementov ima AceGen pripravljene standar- dne podprograme s privzetimi imeni in argumenti. To so npr. podprogram za izra£un tangentne matrike in reziduala, podprogram za postprocesiranje, podprograma za izra£un ob£utljivosti neodvisnega in odvisnih vektorjev re²itve, podprogram za izvedbo poljubne naloge, ki zahteva izvedbo standarne procedure sestavljanja matrik na globalnem nivoju (ang. assembly procedure) itd. V doktorski disertaciji smo za diskretizacijo stohasti£nega polja uporabili zadnjo moºnost. Za izra£un statistike odziva konstrukcije, smo za simula- cije Monte Carlo, kjer je bilo to potrebno, ustrezno raz²irili oz. prilagodili podprogram za izra£un tangentne matrike in reziduala. Za perturbacijsko metodo smo morali ustrezno denirati podprograma, ki sta prilagojena za izra£un ob£utljivostne analize. Ker je AceGen uporabljen za deniranje in generiranje programskih kod, ki so razvite v okviru te disertacije in so deli posameznih kod v kasnej²ih poglavjih tudi predstavljeni, je prav, da se predstavi standardno proceduro za generiranje kode v AceGenu. AceGen je dodatni programski paket Mathematice, zato uporablja jezik Mathematice, raz²irjen z dodatnimi ukazi, ki imajo predpono SMS ter dodatnimi operatorji ($, %, ( in )), 46 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Generator kode AceGen Simbolna izpeljava modela Optimizacija Vpeljava vmesnih Avtomatsko izrazov spemenljivk odvajanje Avtomatsko generiranje kode Okoljski vmesnik · C, C++ · vmesnik za podatke · C# · interpretacija opravil · FORTRAN · splošne numerične rutine · Mathematica script · Matlab script Numerične uporabniške podrutine C/C++/C# Mathematica FORTRAN Matlab AceFEM ELFEN CDriver MDriver ABAQUS FEAP Matlab FE Slika 3.2: Hibridni simbolno-numeri£ni pristop k avtomatizaciji metode stohasti£nih kon£- nih elementov. Figure 3.2: Hybrid symbolic-numeric approach to automation of stochastic nite element method. ki spremenljivki priredijo vrednost in nadome²£ajo standardni operator Mathematice =. Standardna procedura za generiranje kode kon£nega elementa v AceGenu je sestavljena iz naslednjih korakov • SMSInitialize - znotraj tega ukaza se poda ime kode generiranega kon£nega ele- menta, v katerem jeziku naj bo generiran in za katero programsko okolje. • SMSTemplate - tu se poda glavne karakteristike kon£nega elementa (²tevilo vozli²£, topologija, katere vhodne parametre je potrebno podati in njihove privzete vrednosti ipd.) ter pravila za simbolno-numeri£ni vmesnik. • SMSStandardModule - denira za£etek standardnega podprograma s privzetim stan- dardnim imenom in argumenti. Posamezen element ima lahko ve£ razli£nih stan- dardnih podprogramov, ki se jih poda enega za drugim, v poljubnem vrstnem redu. • SMSWrite - ukaz, ki AceGenu sporo£i, da je vhodna datoteka za kon£ni element zaklju£ena in za£ne generiranje datoteke s kodo kon£nega elementa. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 47 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Posamezni deli procedure si morajo slediti v vrstnem redu, kot je predstavljen v zgornjih alinejah. 3.2 Avtomatsko odvajanje Avtomatsko odvajanje je alternativa numeri£nemu in simbolnemu odvajanju. Avtomatsko odvajanje je zasnovano na dejstvu, da vsak ra£unalni²ki program izvaja zaporedje elemen- tarnih aritmeti£nih operacij (se²tevanje, od²tevanje, mnoºenje...) in elementarnih funkcij (log, sin, cos...). Posledi£no se odvodi izra£unajo avtomatsko, s ponavljajo£o uporabo veriºnega pravila na te operacije, dobljene vrednosti pa so analiti£ne, to£ne do natan£no- sti ra£unalni²kega zapisa ²tevil (Griewank, 2000). Da se tehniko avtomatskega odvajanja lahko uporabi za izpeljavo kon£nih elementov, jo je potrebno raz²iriti z dodatnimi ope- ratorji, ki jih imenujemo izjeme avtomatskega odvajanja. Pristop, kjer bomo probleme opisali z uporabo avtomatskega odvajanja bomo imenovali na avtomatskem odvajanju ba- ziran opis problema. Zanj bomo v nadaljevanju uporabljali okraj²avo ADB opis, ki izhaja iz angle²kega izraza automatic dierentiation based. Na tem mestu bomo predstavili kratek povzetek ADB pristopa, podrobnej²i opis je podal Korelc (2009a). Odvajanje poljubne funkcije f po mnoºici medsebojno neodvisnih spremenljivk a z upo- rabo tehnike avtomatskega odvajanja ozna£imo z ˆ δf paq ∇f :“ . (3.1) ˆ δpaq Z operatorjem ˆδp‚q{ˆδp‚q opi²emo, da je za odvajanje funkcije f po spremenljivkah a uporabljena tehnika avtomatskega odvajanja. Obstajata dve moºnosti avtomatskega odvajanja: forward in backward na£in. Prikaz obeh moºnosti avtomatskega odvajanja si poglejmo na naslednjem primeru (povzeto po zapiskih semi-plenarnega predavanja Korelc, 2012). Imamo funkcijo f f “ b c, (3.2) n ÿ b “ a2 in (3.3) i i“1 c “ sinpbq, (3.4) kjer so ai neodvisne spremenljivke. Forward na£in akumulira odvode vmesnih spremen- 48 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. ljivk po neodvisnih spremenljivkah " db * ∇b “ “ t2aiu, i “ 1, 2, ..., n, (3.5) dai " dc * ∇c “ “ tcos(b)∇biu, i “ 1, 2, ..., n, (3.6) dai " df * ∇f “ “ t∇bi c ` b∇ciu, i “ 1, 2, ..., n. (3.7) dai Backward na£in avtomatskega odvajanja propagira pridruºene spremenljivke ¯x “ Bf , ki Bx so odvodi kon£nih vrednosti po vmesnih spremenljivkah ¯ df f “ “ 1, 1, (3.8) df df Bf ¯ c “ “ ¯ f “ b ¯ f , 1, (3.9) dc Bc ¯ df Bf Bc b “ “ ¯ f ` ¯ c “ c ¯ f ` cospbq ¯c, 1, (3.10) db Bb Bb " db * ∇f “ t¯ a ¯ ¯ iu “ b “ 2ai b( , i “ 1, 2, ..., n. (3.11) dai Pri odvajanju N skalarnih funkcij F “ tfi, i “ 1, ..., Nu po M neodvisnih spremenljivkah a “ taj, j “ 1, ..., M u je numeri£na u£inkovitost forward na£ina v splo²nem sorazmerna s ²tevilom neodvisnih spremenljivk M, medtem ko je backward na£in sorazmeren s ²te- vilom funkcij N. V primerih, ko je ²tevilo funkcij majhno, je backward na£in superioren nad forward na£inom. Njegova slabost je le, da potrebuje potencialno ve£ prostora za shranjevanje vmesnih podatkov med izvrednotenjem funkcije, ki je lahko enako ²tevilu iz- vedenih numeri£nih operacij. Za u£inkovito uporabo avtomatskega odvajanja je zaºeleno, da je moºna uporaba obeh na£inov. AceGen omogo£a uporabo obeh na£inov, pri £emer v posameznem klicu procedure avtomatskega odvajanja samodejno izbere na£in odvajanja glede na razmerje ²tevila funkcij in spremenljivk. Kot je bilo ºe predstavljeno, je tehniko avtomatskega odvajanja za uporabo v metodi kon£nih elementov potrebno raz²iriti z dodatnimi operatorji, ki jih imenujemo izjeme avtomatskega odvajanja. V nadaljevanju so predstavljeni formalizmi za vpeljavo razli£nih tipov izjem in razli£nih na£inov vpeljave izjem, lokalno in globalno deniranih (povzeto po Korelc, 2009a). Vpeljavo lokalno denirane izjeme k proceduri avtomatskega odvajanja ozna£imo z ˆ ˇ δf pa, bpaqq ˇ ∇f :“ ˇ , (3.12) ˆ δ ˇ paq Ďpbq D “M paq kjer je f poljubna funkcija, a je mnoºica medsebojno neodvisnih spremenljivk, b je mno- ºica medsebojno neodvisnih vmesnih spremenljivk, ki so del izvrednotenja funkcije f, G je Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 49 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. mnoºica poljubnih funkcij, odvisnih od a, tako da je b :“ Gpaq, in M je poljubna matrika. Izjema Dpbq “ M denira, da je odvod poljubne mnoºice medsebojno neodvisnih vme- Dpaq snih spremenljivk b po neodvisnih spremenljivkah a (ne glede na dejansko medsebojno odvisnost ali neodvisnost teh spremenljivk) enak matriki M. V primeru, da za odvajanje funkcije f po a ni denirana nobena izjema, je rezultat procedure avtomatskega odvajanja enak parcialnemu odvodu f po a ˆ δf paq Bf paq “ . (3.13) ˆ δpaq B paq Lo£imo ²tiri osnovne tipe izjem avtomatskega odvajanja, ki so prikazane v preglednici 3.1. Tip A je predstavljen ºe v ena£bi (3.12). Uporaba tega tipa izjeme avtomatskega odvajanja je lahko numeri£no ugodna tudi v primeru, ko odvisnost b od a dejansko obstaja, vendar se na ta na£in izognemo ra£unanju odvodov b po a preko veriºnega pravila. Tipi izjem B, C in D predstavljajo posebne primere izjeme tipa A. Tip izjeme B se uporabi, ko v algoritmu programske kode sicer obstaja odvisnost b od a, vendar se zahteva, da se med proceduro avtomatskega odvajanja te odvisnosti ne upo²teva. Ta tip izjeme pride v mehanskih problemih v po²tev npr. ko je namesto totalne variacije potrebno izvrednotiti poljubno variacijo obravnavane koli£ine. Tip izjeme C pride v po²tev, ko v algoritmu ne obstaja odvisnost b od a, vendar iz formulacije problema sledi, da obstajajo (obi£ajno implicitne) odvisnosti, ki jih je pri odvajanju potrebno upo²tevati. Izjemo tipa A se lahko tudi posplo²i. Naj bo c mnoºica vmesnih spremenljivk, ki so del izvrednotenja funkcije f, tako da je b :“ Gpcq in H mnoºica poljubnih funkcij, odvisnih od a, tako da je c :“ Hpaq. Na ta na£in se v veriºnem pravilu lahko naredi most, ki premosti vmesno pot med b in c. Veriºno pravilo, ki povezuje c z a se nato propagira avtomatsko s proceduro avtomatskega odvajanja. Izjema avtomatskega odvajanja, ki se ti£e mnoºice vmesnih spremenljivk b je lahko vpe- ljana na dva na£ina: • Lokalna denicija izjeme: izjema (ena ali ve£) je denirana neposredno na mestu, kjer se pokli£e proceduro avtomatskega odvajanja, kot je prikazano v (3.12) in v srednjem stolpcu v preglednici 3.1. • Globalna denicija izjeme: podatek o izjemi avtomatskega odvajanja, ki se ti£e mnoºice vmesnih spremenljivk b se poda v algoritmu na mestu, kjer je mnoºica vmesnih spremenljivk b vpeljana, kar za splo²en primer lahko zapi²emo b :“ Gpaq | Dpbq D “M paq ˆ δf pa, bpaqq ∇f :“ (3.14) ˆ δpaq 50 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. in je za posamezne tipe izjem prikazana v desnem stolpcu v preglednici 3.1. Preglednica 3.1: Tipi£ne izjeme avtomatskega odvajanja. Table 3.1: Typical automatic dierentiation exceptions. tip Lokalna denicija Globalna denicija izjeme izjeme izjeme ˆ ˇ A ∇f δf pa,bpaqq ˇ b :“ Gpaq| D A :“ pbq ˆ δpaq Ďpbq D “M paq D “M paq ˆ ∇f δf pa,bpaqq A :“ ˆ δpaq ˆ ˇ B ∇f δf pa,bpaqq ˇ b :“ Gpaq| D B :“ pbq ˆ δpaq Ďpbq D “0 paq D “0 paq ˆ ∇f δf pa,bpaqq B :“ ˆ δpaq ˆ Č ∇f δf pa,bq ˇ b :“ Gpaq| D C :“ pbq ˆ δpaq Ďpbq D “M paq D “M paq ˆ ∇f δf pa,bpaqq C :“ ˆ δpaq c : ˆ ˇ “ Hpaq D ∇f δf pa,bpcpaqqq Ď :“ ˆ δpaq Ďpbq D “M pcq b :“ Gpcq| Dpbq D “M pcq ˆ ∇f δf pa,bpcpaqqq D :“ ˆ δpaq Izbira med enim in drugim na£inom vpeljave izjem je odvisna od posameznega primera. V splo²nem je bolj prikladna uporaba globalnih denicij izjem, ²e posebej, ko imamo opravka s kompleksnej²imi ra£unskimi modeli, v katerih je med izvajanjem algoritma potrebno ve£krat ra£unati odvode. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 51 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 4 Avtomatizacija diskretizacije stohasti£- nega polja 4.1 Stohasti£ni kon£ni elementi V okviru doktorske disertacije je bila razvita avtomatizacija postopka diskretizacije sto- hasti£nega polja v K-L vrsto, predstavljenega v poglavju 2.1. V ta namen smo razvili stohasti£ne kon£ne elemente. Koda stohasti£nega kon£nega elementa je sestavljena iz treh podprogramov: • podprogram za izra£un prispevka posameznega stohasti£nega kon£nega elementa h globalni matriki oblikovnih funkcij N • podprogram za izra£un prispevka posameznega stohasti£nega kon£nega elementa h globalni kovarian£ni matriki C • podprogram za izra£un integrala produkta dveh lastnih funkcij (leva stran ena£be (2.9)) Procedura avtomatizacije diskretizacije stohasti£nega polja v KL vrsto je prikazana na sliki 4.1. DISKRETIZACIJA STOHASTIČNEGA POLJA Diskretizacija domene D z D poddomenami k Definicija 2× N-vozliščnih stohastičnih končnih elementov Izračun matrik C in N Rešitev posplošenega problema lastnih vrednosti Cf = ΛNf Normiranje lastnih funkcij f w  M X,   w  X     f X k k   k  k 1 Slika 4.1: Procedura avtomatizacije diskretizacije stohasti£nega polja. Figure 4.1: Procedure of automation of stochastic eld discretization. 52 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Glavni problem, ki ga je bilo potrebno re²iti v procesu formulacije stohasti£nih kon£nih elementov je izra£un prispevka posameznega stohasti£nega kon£nega elementa h kova- rian£ni matriki C, saj je v tem izra£unu potrebno poznati razdaljo med posameznimi elementi. Klasi£no formulirani kon£ni elementi informacije o koordinatah preostalih ele- mentov ne vsebujejo, zato bi, £e bi ºeleli re²iti problem na ta na£in, morali izvesti temeljne posege v okolje za kon£ne elemente. Temu smo se ºeleli izogniti in problem re²iti tako, da je re²itev primerna in uporabna v poljubnem standardnem okolju za kon£ne elemente, ki omogo£a odprto kodo kon£nih elementov. Tako se nam je porodila ideja denicije novih kon£nih elementov po naslednjem postopku. Najprej smo zikalno domeno D razdelili na nD poddomen, kot je prikazano na sliki 4.2 nD ď D “ Dk (4.1) k“1 ... DnD ... Dk ... D D D1 D2 D3 ... Slika 4.2: Razdelitev zikalne domene D na nD poddomen. Figure 4.2: Discretization of physical domain D into nD subdomains. V naslednjem koraku se vsaka poddomena Dk, k “ 1, 2, ..., nD kombinira z vsako posame- zno poddomeno Dl P D, l “ k, k`1, ..., nD, vklju£no s samo sabo. Vsaka izmed kombinacij dveh poddomen predstavlja 2ˆN-vozli²£ni stohasti£ni kon£ni element. Primer takega sto- hasti£nega kon£nega elementa za dvodimenzionalno domeno D, diskretizirano z linearnimi ²tirivozli²£nimi elementi je prikazan na sliki 4.3. V primeru, da je polje nizko korelirano in se kovarian£na funkcija poenostavi tako da se upo²teva le do neke omejene razdalje med elementi, ki jo ozna£imo z le (kar je podrobneje razloºeno kasneje v poglavju 5), se vsako izmed poddomen kombinira le s poddomenami, ki so od domene Dk oddaljene manj ali enako od le, kot je prikazano na sliki 4.4. Po obi£ajnem postopku, ki se uporablja v metodi kon£nih elementov, se nato koli£ine oz. prispevki posameznih stohasti£nih kon£nih elementov izra£unajo na referen£nih kon£nih elementih. Referen£na kon£na elementa, na katerih se izra£una prispevke, sta oblike kvadrata z dolºino stranice 2, kot je prikazano na sliki 4.5. Za prera£unavanje koli£in med globalnim in referen£nim koordinatnim sistemom se uporabi izoparametri£na preslikava. Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 53 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 2×4-vozliščni stohastični končni element: 8 7 Dl Dl D 5 6 k 4 3 D Dk 1 2 Slika 4.3: Diskretizacija dvodimenzionalne domene D in formulacija 2ˆ4-vozli²£nih sto- hasti£nih kon£nih elementov. Figure 4.3: Discretization of twodimensional domain D and formulation of 2ˆ4-noded stochastic nite elements. 2×4-vozliščni stohastični končni element: 8 7 Dl Dl Dk 5 6 4 3 leff D Dk 1 2 Slika 4.4: Diskretizacija dvodimenzionalne domene D in formulacija 2ˆ4-vozli²£nih sto- hasti£nih kon£nih elementov na razdalji le, z upo²tevanjem simetrije problema. Figure 4.4: Discretization of two-dimensional domain D and formulation of 2ˆ4-noded stochastic nite elements on le region with consideration of problem simmetry. Ker originalni in referen£ni kon£ni element v splo²nem nista enake oblike in dimenzij, je potrebno koli£ine, izra£unane v referen£nem koordinatnem sistemu, pri prevedbi v globalni koordinatni sistem pomnoºiti z determinanto Jakobijeve matrike preslikave. Prispevek posameznega stohasti£nega kon£nega elementa h globalni kovarian£ni matriki C in matriki oblikovnih funkcij N se znotraj stohasti£nega kon£nega elementa izra£una z uporabo numeri£ne integracije Ng Ng ÿ ÿ Ce w , (4.2) ij “ g whC gh ij g“1 h“1 Ng ÿ N e w , (4.3) ij “ g N g ij g“1 54 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. 8 7 8 7 (-1,1) (1,1) 2 5 Dl 2 (-1,-1) (1,-1) 3 4 6 5 6 4 3 (-1,1) (1,1)  D 1 k 1 Y 2 1 (-1,-1) (1,-1) X 1 2 Slika 4.5: Izoparametri£na preslikava iz globalnega v referen£ni koordinatni sistem za primer dvodimenzionalnega, 2ˆ4-vozli²£nega stohasti£nega kon£nega elementa. Figure 4.5: Isoparametric mapping from global to reference coordinate system for two- dimensional 2ˆ4-noded stochastic nite element. kjer je Ng ²tevilo integracijskih to£k, wg in wh sta uteºi Gaussovih integracijskih to£k g in h ter Cgh in Ng sta prispevka k matrikama C in N, izra£unana v Gaussovih to£kah po ij ij ena£bi Cgh ij “ JgJhσ2C pXg, Xhq Nj pXgq Ni pXhq , N gij “ JgNi pXgq Nj pXgq , (4.4) kjer sta Jg in Jh determinanti Jacobijeve matrike preslikave referen£nih koordinat ξ1, η1 oz. ξ2, η2 v globalni koordinatni sistem. Ra£unanje prispevka Cgh se izvede po vseh ij stohasti£nih kon£nih elementih, medtem ko se izra£un prispevka Ng izvede le, £e sta ij indeksa k-te in l-te domene, ki sta kombinirani v stohasti£nem kon£nem elementu (Dk in Dl), enaka. Razlog je v tem, da je za izra£un matrike N potrebno izra£unati le enojni integral produktov oblikovnih funkcij. Integral je sicer potrebno izvesti preko celotne domene, vendar so oblikovne funkcije v metodi kon£nih elementov denirane tako, da so le lokalno razli£ne od ni£ (npr. oblikovna funkcija Ni je razli£na od ni£ le v poddomenah, ki se stikujejo v vozli²£u i), zato je produkt dveh oblikovnih funkcij NipXgq in NjpXgq razli£en od ni£ le v primeru, da vozli²£i i in j pripadata isti poddomeni. Potrebno je torej narediti zanko, ki gre le po vsaki poddomeni enkrat. Zato se, da se izognemo ve£kratnemu integriranju po isti domeni, izra£un prispevka Ng izvede le, £e sta v stohasti£nem elementu ij kombinirani isti poddomeni (Dk “ Dl). Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 55 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Globalni matriki C in N sta nato formirani s standardno proceduro sestavljanja po kon£- nih elementih ď C “ Ce, (4.5) e ď N “ Ne. (4.6) e Ob tem je potrebno dodati ²e eno opombo. Lastnost, da je matrika C simetri£na, smo izkoristili za izbolj²anje numeri£ne u£inkovitosti, se pravi da se sestavi le zgornja trikotna matrika C in N. Ker so domene Dk in Dl kombinirane po vnaprej dolo£enem zaporedju (k “ 1, 2, ..., nD, l “ k, k ` 1, nD), je bilo v okviru formulacije kode za proceduro sesta- vljanja prispevkov Cij posameznih stohasti£nih elementov v globalno matriko potrebno upo²tevati, da se ob upo²tevanju simetri£nosti matrike sestavlja le zgornji trikotnik ma- trike. Indeks i naredi namre£ zanko po vozli²£ih domene Dk, indeks j pa po vozli²£ih domene Dl. Izra£unani prispevek Ce se v tem primeru namesto po izrazu (4.2) izra£u- ij najo po izrazu N # g Ng ÿ ÿ 2 i “ j ^ k ‰ l Ce w ; a (4.7) ij “ g whaC gh ij “ 1 sicer g“1 h“1 ƒe pj ă iq ^ pk ‰ lq, shrani Ce na mesto Ce ij ji Izra£unani prispevek Cgh je zaradi vnaprej dolo£enega vrstnega reda domen D ij k in Dl potrebno pri²teti k Ce v primeru da je indeks j manj²i od indeksa i in da domeni D ji k in Dl nista isti (k ‰ l), saj se pri simetri£ni tangenti sestavlja le zgornji trikotnik tangentne matrike in se tega prispevka sicer ne bi upo²tevalo. Druga posebnost, ki jo je potrebno upo²tevati ob simetri£ni tangenti, pa je, da v primeru, ko domeni Dk in Dl nista isti (k ‰ l), imata pa skupno vozli²£e (v tem primeru i “ j), je prispevek (zopet zaradi vnaprej dolo£enega vrstnega reda indeksov kombinacij domen Dk in Dl, ker se indeksa ne zamenjata) potrebno upo²tevati dvakrat, saj bi se sicer izgubil del prispevka, ki bi ga imela kombinacija domene Dl ´ Dk (in je zaradi i “ j enak, od tu faktor 2). Ko sta globalni matriki N in C sestavljeni, sledi izra£un lastnih vrednosti in lastnih funk- cij. Ta izra£un je opravljen v Mathematici, ki ima vgrajen algoritem ra£unanja problema lastnih vrednosti in vektorjev in omogo£a tudi nekaj dodatnih opcij za izbolj²anje nume- ri£ne u£inkovitosti in natan£nosti, npr. Arnoldijev algoritem, ki v primeru, da nas zanima le majhno ²tevilo £lenov K-L vrste in je matrika C razpr²ena in pozitivno denitna, re²eva- nje problema lastnih vrednosti ob£utno poceni. Na tem mestu morda ni odve£ opomniti, da so vrednosti dobljenih lastnih vektorjev zaradi narave oblikovnih funkcij, ki imajo vre- dnost 1 v enem vozli²£u in 0 v vseh ostalih, ravno vrednosti (sicer ²e nenormiranih) lastnih 56 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. funkcij v vozli²£ih (ena£ba (2.25)). Izra£unani lastni vektorji oz. lastne funkcije so v splo²nem nenormirane, saj je njihova norma izbrana avtomati£no s strani numeri£ne knjiºnice. Lastne funkcije je zato potrebno normirati skladno z ena£bo (2.9). Tretji podprogram stohasti£nega kon£nega elementa izra£una prispevek posameznih elementov k izra£unu integrala ż Ifk “ fkpXqfkpXqdX. (4.8) D Zaradi ortogonalnosti lastnih funkcij je potrebno kombinirati le produkt lastnih funkcij z istim indeksom. Za integriranje kot obi£ajno v metodi kon£nih elementov uporabimo Gaussovo numeri£no integracijo, po kateri izra£unamo vrednost izraza, ki ga integriramo, v to£no dolo£enih to£kah. Te vrednosti nato pomnoºimo z ustreznimi uteºmi Ng ÿ Ie w , f k “ g I g f k g“1 ˜ J ¸2 ÿ Ig f , (4.9) f k “ Jg jk,nenormNj pXgq j“1 ď Ifk “ Ie . f k e Prispevek Ig se izra£una le v tistih stohasti£nih kon£nih elementih, pri katerih sta kombi- f k nirani dve domeni z istim indeksom (k=l). Razlog je enak kot pri izra£unu matrike N in sicer se na ta na£in izognemo ve£kratnemu integriranju po isti domeni. Skladno z ena£bo (2.9) se lastne funkcije nato normira po izrazu f f k,nenorm k “ . (4.10) aIfk Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 57 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Algoritem 1 ADB formulacija stohasti£nega kon£nega elementa pDk Y Dlq if task = Nij and k “ l then Ne :“ 0 for g :“ 1, ..., ng do Ź Zanka po integracijskih to£kah ´ ˆ ¯ J δXg g :“ det ˆ Ź Jacobijeva determinanta δΞg for i :“ 1, ..., J; j :“ i, ..., J do N g : : ij “ JgNi pXgq Nj pXgq; N e ij “ N e ij ` wg N g ij end for end for end if if task = Cij then Ce :“ 0 for g :“ 1, ..., ng; h :“ 1, ..., ng do Ź Zanki po integracijskih to£kah na Dk in Dl ´ ˆ ¯ ´ ˆ ¯ J δXg δXh g :“ det ; J ˆ h :“ det δΞ ˆ g δΞh for @i, j P r1, 2, ..., Js do Cgh : ij “ JgJhσ2C pXg, Xhq Nj pXgq Ni pXhq Ź Prispevek Gaussove to£ke k Cij a :“ If pi “ j ^ k ă lq, then 2, else 1 If pj ă i ^ k ‰ lq, then Ce : , else Ce : ji “ C e jì wg whC gh ji ij “ Ce ij ` wg whaC gh ij end for end for end if if task = Ifk and k “ l then Ie : f k “ 0 for g :“ 1, ..., ng do ´ ¯ J δXg g :“ det ; Ig :“ f :“ Ie δΞ k .fk Jg ; Ie g f k f k f k ` wg I g f k Ź Prispevek Gaussove to£ke k Ifk end for end if 4.2 Izvedba diskretizacije stohasti£nega polja v okolju za kon£ne elemente Ace- FEM V tem poglavju je prikazan primer vhodne datoteke za izvedbo diskretizacije dvodimen- zionalnega stohasti£nega polja. Na sliki 4.6 je prikazana oblika domene problema, preko katere je denirano stohasti£no polje. Podatki stohasti£nega polja so podani v preglednici 4.1. Fizikalna domena problema je v prvem koraku diskretizirana s 600 poddomenami (15 razdelkov v navpi£ni smeri in 40 razdelkov v vodoravni smeri). Preglednica 4.1: Parametri stohasti£nega polja. Table 4.1: Parameters of stochastic eld. Pri£akovana vrednost ( ¯ w) 0 Standardna deviacija (σw) 1 Korelacijska dolºina (lc) 0.5 Kovarian£na funkcija CC, le “ 3lc 58 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. .57 .53 8 Slika 4.6: Skica zikalne domene, preko katere je denirano stohasti£no polje za prikaz izvedbe diskretizacije v AceFEM-u. Figure 4.6: Sketch of physical domain on which twodimensional stochastic eld is dened for demonstration of AceFEM procedure. Generiranje mreºe pomoºnih generi£nih 4-vozli²£nih Q1 elementov << AceFEM‘; L = 10.; R1 = L  4; R2 = L  2; h = 2.5; m1 = FlattenB MapThreadBTable, :::x + L  4, h  2 + R12 - x2 >, :x + 3 L  4, h  2 - R12 - x2 >>, 88x, -L  4 + L  10, L  4, L  20<, 8x, -L  4 + L  20, L  4 - L  10, L  20<<>F, 1F; m3 = TableB:x + L  2, - h  2 - R22 - x2 >, 8x, - L  2 + L  10, L  2 - L  10, L  20F; 2 mastermesh = 8m3, m2, m1<; SMTInputData@D; SMTAddDomain@"D", "Q1", 8 Σw, "lc*" -> lc, "leff*" ® leff "AceFEM"D; ngh = 5; lgh = ngh + 1 + ngh idata$$@"NoSensParameters"D; leh = lgh * es$$@"id", "NoIntPoints"D; SMSTemplate@"SMSTopology" ® "Q1", "SMSSymmetricTangent" ® False, "SMSNoTimeStorage" ® leh, "SMSPostIterationCall" ® True, "SMSGroupDataNames" -> 8"E -elastic modulus", "Ν -poisson ration", "t -thickness", "Σy -initial yield stress", "Kh -hardening coefficient", "ΣyInf -residual flow stress", "∆ -saturation exponent", "Ρ0 -density", "bX -force per unit mass X", "bY -force per unit mass Y"<, "SMSShapeSensitivity" -> True, "SMSDefaultData" -> 821 000, 0.3, 1, 24, 0, 24, 0, 1, 0, 0<, "SMSIntSwitch" -> Table@0, 8maxSensPar HmaxSensPar + 3L hgD; H * en. H6.21L *L Rgti £ SMSD@Rg, Ξi, "Dependency" -> 88hg, Ξi, Zgi< TrueD; , 8k, 1, 8 hgD; H * en. H6.43L *L Rgti £ SMSD@SMSD@Rg, ΞiD, Ξj, "Dependency" ® 88DhgDΞiΞj@@1DD, Ξj, Zgij< TrueD; , 8k, 1, 8 NumberQD; 8Λ, Μ< £ SMSHookeToLame@Em, ΝD; wgp ¢ SMSReal@es$$@"IntPoints", 4, IgDD; hgP1T hgP4T 0 Cgpi = IdentityMatrix@3D + hgP4T hgP2T 0 ; 0 0 hgP3T Γ = hgP5T; SMSFreeze@be, F.Cgpi.Transpose@FD, "Ignore" ® NumberQ, "Symmetric" -> TrueD; Jbe £ Det@beD; 1 1 1 W = Μ HTr@beD - 3L - Μ Log@JbeD + Λ HJbe - 1 - Log@JbeDL; 2 2 4 SMSFreeze@Τ, Simplify@2 be.SMSD@W, be, "Ignore" ® NumberQ, "Symmetric" ® TrueDD, "Symmetric" ® True, "Ignore" -> NumberQD; 1 s = Τ - IdentityMatrix@3D Tr@ΤD; 3 Σy £ HΣy0 + Kh Γ + HΣyInf - Σy0L H1 - Exp@- ∆ ΓDLL; 3 ΣMises = SMSSqrtB Total@s s, 2DF; 2 F £ ΣMises - Σy; hgnP1T hgnP4T 0 Cgpin = IdentityMatrix@3D + hgnP4T hgnP2T 0 ; 0 0 hgnP3T Γn = hgnP5T; A = SMSD@F, Τ, "Ignore" ® NumberQ, "Symmetric" ® TrueD; M £ - 2 HΓ - ΓnL A; Z £ Simplify@F.Cgpi - SMSMatrixExp@M, "Order" ® 4D.F.CgpinD; Qg £ 8ZP1, 1T, ZP2, 2T, ZP3, 3T, ZP1, 2T, F<; 6.7 Izvedba ob£utljivostne analize v okolju za kon£ne elemente AceFEM V tem poglavju je predstavljena procedura za izvedbo ob£utljivostne analize v okolju za kon£ne elemente AceFEM. Za prikaz procedure je izbran primer elasti£nega upogiba konzole. Skica primera je prikazana na sliki 6.1, materialni in geometrijski podatki pa so 112 Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. podani v preglednici 6.8. Izbran je linearno elasti£en materialni model. I²£e se navpi£ni pomik na prostem robu konzole vA ter ob£utljivost pomika vA na dolºino konzole L ter na vi²ino pre£nega prereza konzole h. Zaradi laºjega razumevanja predstavljene procedure ra£unanja ob£utljivostne analize, parametra L in h nista stohasti£ne narave. y P t A h v x A L Slika 6.1: Upogib konzole. Figure 6.1: Bending of cantilever beam. Preglednica 6.8: Materialni in geometrijski parametri. Table 6.8: Material parameters and geometry. E 1000 GPa ν 0.0 L 1000 mm h 10 mm t 10 mm P 1 kN Melink, T. 2014. Metoda stohasti£nih kon£nih elementov v modeliranju konstrukcij. 113 Doktorska disertacija. Ljubljana, UL FGG, Podiplomski ²tudij gradbeni²tva, Konstrukcijska smer. Denicija domene, robnih pogojev ter ob£utljivostnih parametrov << AceFEM‘; L = 1000; h = 10; t = 10; P = - 1; Em = 1000; Ν = 0.0; nx = 25; ny = 1; SMTInputData@D; SMTAddDomain@"Cantilever", "HookeQ2", 8"E *" ® Em, "Ν *" ® Ν, "t *" ® t 0.<, 8"X" == 0 && "Y" Š 0 &, 2 -> 0.<