Gradbeni vestnik • letnik 65 • maj 2016 109 METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ•Teja Melink, Jože Korelc METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ STOCHASTIC FINITE ELEMENT METHOD IN STRUCTURE MODELLING asist. dr. Teja Melink, univ. dipl. inž. grad. teja.melink@guest.arnes.si Na Pristavi 51, 5290 Šempeter pri Gorici prof. dr. Jože Korelc, univ. dipl. inž. grad. joze.korelc@fgg.uni-lj.si Fakulteta za gradbeništvo in geodezijo, Jamova 2, 1000 Ljubljana Znanstveni članek UDK UDK 519.21:624.07 Povzetek l V članku je predstavljen stohastični pristop k reševanju mehanskih problemov. Stohastični pristop privzame, da spremenljivke nimajo določene vrednosti, temveč z neko verjetnostjo zavzamejo naključno vrednost, vrednosti pa se lahko naključno spreminjajo tudi preko prostorske domene problema. V okviru raziskovalnega dela smo razvili pristop za modeliranje in analizo stohastičnih procesov v okolju za metodo končnih elementov, ki se na ta način razširi v metodo stohastičnih končnih elementov. Posledično ima razvita metoda prednosti, ki odlikujejo metodo končnih elementov, kot sta na primer možnost uporabe za poljubne geometrijsko in materialno nelinearne probleme ter eno- staven opis robnih pogojev problema. Glavna cilja pri razvoju metode sta bila zanesljivost in numerična učinkovitost. Računski časi razvitega pristopa za reševanje stohastičnih problemov so primerljivi z običajnim determinističnim pristopom k analizi konstrukcij, zaradi česar je prikazan stohastični pristop uporaben tudi za praktično projektiranje. Ključne besede: stohastični pristop, verjetnost, metoda stohastičnih končnih elementov, perturbacijska metoda, občutljivostna analiza Summary l The paper presents the stochastic approach in mechanical problems. In stochastic approach, the uncertainties of parameters involved in mechanical problems are taken into account. The developed stochastic approach is compatible with finite ele- ment method, thus allowing application to arbitrary geometrical and material nonlinear mechanical problems, simple manipulation with arbitrary input data and mesh topology. Further on, it is numerically efficient. Two main objectives of our work were the numerical stability and numerical efficiency of the methods. The calculation time of the presented stochastic approach is comparable to that of the deterministic approach to mechanical problems with the use of standard finite element method. Key words: stochastic approach, uncertainty, stochastic finite element method, perturba- tion method, sensitivity analysis. 1•UVOD V inženirskih problemih imamo opravka z naključnostjo pri skoraj vseh parametrih pro- blema. Naključno se spreminjajo lastnosti ma- teriala, obtežba, geometrijske nepopolnosti, korozija, 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. Zaradi tega je stohastični pristop k mode- liranju konstrukcij, ki upošteva, da vhodni para- metri problema nimajo točno določenih vred- nosti, temveč z neko verjetnostjo zavzamejo naključno vrednost na omejenem ali neomeje- nem območju, bolj realen v primerjavi s splošno uveljavljenim determinističnim pristopom, pri Gradbeni vestnik • letnik 65 • maj 2016110 Teja Melink, Jože Korelc•METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ katerem ima vsak vhodni parameter točno določeno (determinirano) vrednost. Kljub temu je stohastični pristop le redko uporabljen v praksi, saj upoštevanje naključnosti pojavov bistveno poveča numerično zahtevnost in kom- pleksnost problema. Zaradi tega se v inženirski praksi večinoma uporablja enostavnejši, deterministični pristop, vse možne scenarije zaradi stohastične narave pojavov pa se poskuša zajeti z upoštevanjem ekstremov ali srednjih vrednosti stohastičnih para- metrov. Vendar so študije z bolj realističnim upoštevanjem stohastičnosti pokazale, da deterministični pristop kljub uporabi ekstrem- nih vrednosti parametrov, ki so na videz na varni strani, ne zajame nujno najslabšega možnega scenarija, ki se lahko zgodi zaradi naključne kombinacije posameznih para- metrov. Poleg tega deterministični pristop ne vodi do optimalnih dimenzij inženirskih kon- strukcij. Zaradi tega je dolgoročno korak v smeri stohastičnega modeliranja neizogiben. V raziskovalnem delu smo se zato posvetili razvoju avtomatizacije stohastičnega pristopa k modeliranju konstrukcij. Glavni cilji pri tem so bili: • računski čas naj bo primerljiv deter- minističnemu pristopu, • metoda naj bo robustna in • sledi naj standardni tehnologiji končnih ele- mentov. 1.�korak:����������������������������������� 2.�korak:�������������������� ������ ����������� ADB�FORMUL�RANA OB�UTLJ��OSTNA� ANAL�ZA R�zv���v�kt�r�����z�v��v�T�y��r��v��vr�t���kr��� �r���k�v���h�vr�����t�. R�zv����t�h��t�������������v�K�rh����-L�èv�v��vr�t�.� S�MBOLN�� PR�STOP METODA�KON�N��� D�FERENC PERTURBAC�JSKA�METODAS�MULAC�JE MONTE�CARLO� P����b����t�v����r����z�c��� �t�h��t�������������z� ���r�b������r�t�r��� ��k������h��t�v��.� �zr�����v�kt�r�����z�v��z�� v��k��r����z�c���� �t�h��t������������. 1.�korak:����������������������������������� 2.�korak:�������������������� ������ ����������� ADB�FORMUL�RANA OB�UTLJ��OSTNA� ANAL�ZA R�zv���v�kt�r�����z�v��v�T�y��r��v��vr�t���kr��� �r���k�v���h�vr�����t�. R�zv����t�h��t�������������v�K�rh����-L�èv�v��vr�t�.� S�MBOLN�� PR�STOP METODA�KON�N��� D�FERENC PERTURBAC�JSKA�METODAMONTE�CARLO� S�MULAC�JE P����b����t�v����r����z�c��� �t�h��t�������������z� ���r�b������r�t�r��� ��k������h��t�v��.� �zr�����v�kt�r�����z�v��z�� v��k��r����z�c���� �t�h��t������������. Slika 1•Procedura predlagane metode stohastičnih končnih elementov Omenjene lastnosti so namreč ključne, da je stohastičen pristop uporaben tudi za praktično projektiranje konstrukcij. Pri razvoju formula- cije stohastičnega pristopa smo iskali metode in numerične postopke, ki so kompatibilni z metodo končnih elementov, saj je le-ta najširše uporabljena metoda v inženirski praksi in ima zato takšen pristop možnost širše uporabe. Na ta način se metoda končnih elementov razširi v metodo stohastičnih končnih elementov. Procedura metode stohastičnih končnih ele- mentov je prikazana na sliki 1 in je v grobem sestavljena iz dveh korakov: 1. reprezentacije stohastičnega polja in 2. izračuna statistike odziva konstrukcije. V nadaljevanju je predstavljena razvita avto- matizacija stohastičnega pristopa za vsak korak posebej. Pristop in postopki so pred- stavljeni na kratko, za podrobnejši opis glej [Melink, 2014b]. Vrednost spremenljivk, ki se stohastično spreminjajo po prostoru, se lahko opiše s stohastičnim poljem. V raziskavi smo se omejili na stohastična polja mehanskih spremenljivk z Gaussovo robno porazdelitvijo, vendar se predstavljen način obravnave stohastičnega polja lahko uporabi tudi za druge porazde- litve ali spremenljivke, le z nekaj dodatnimi računskimi operacijami. Stohastično polje je določeno s pričakovano vrednostjo in kovariančno funkcijo [Ghanem, 2003]. Kovariančna funkcija določa korelacijo vrednosti, tj., koliko so vrednosti naključne spremenljivke v posameznih točkah v pro- storu med seboj povezane. Velika korelacija polja pomeni, da se vrednosti naključne spremenljivke preko prostora zelo malo spre- minjajo, medtem ko pri majhni korelaciji vred- nosti naključne spremenljivke preko domene močneje variirajo. Za boljšo predstavo, kaj korelacija pomeni, so na sliki 2 prikazane 2•REPREZENTACIJA STOHASTIČNEGA POLJA štiri naključno izbrane realizacije dvodimen- zionalnega stohastičnega polja z različnimi korelacijami. Koreliranost stohastičnega polja na tej sliki od primera (a) proti primeru (d) narašča. Primer (a) prikazuje nizko kore- lirano stohastično polje: vrednosti naključne spremenljivke v tem primeru precej variirajo preko prostorske domene. Nasprotni ek- strem je neskončno korelirano stohastično polje, ki ga prikazuje primer (d). Vrednost stohastične spremenljivke je v tem primeru še vedno naključna (in v splošnem različna v vsaki realizaciji), vendar konstantna preko celotne prostorske domene. Vmesni sliki (b) in (c) pa prikazujeta srednje in visoko kore- lirano stohastično polje. V praksi predstavlja močno korelirano polje npr. elastični modul pri posameznem odlitku betona, srednje ko- relirana polja srečamo npr. pri geometrijski nepopolnosti konstrukcij ali pojavu korozije, nizko korelirana polja pa odgovarjajo spre- membam, ki so vezane na kristalno strukturo materialov. Za opis koreliranosti spremenljivk mehanskih problemov se običajno uporablja eksponentna kovariančna funkcija v obliki (1) kjer sta x1 in x2 vektorja prostorskih koordi- nat poljubnih dveh točk na stohastičnem polju, σ2 je varianca stohastičnega polja in lc korelacijska dolžina, ki določa koreliranost stohastičnega polja. Eksponentna kovariančna funkcija in vpliv različnih korelacijskih dolžin na obliko kovariančne funkcije prikazuje slika 3. Ker je stohastično polje zvezna veličina, ga je treba za numerično obravnavo diskretizirati. Za opis in diskretizacijo stohastičnega polja smo izbrali Karhunen-Loèvovo (K-L) dekompozicijo, ki stohastično polje razvije v vrsto in predstav- lja optimalno redukcijo v smislu potrebnega števila členov vrste [Ghanem, 2003]. Po- leg tega je K-L dekompozicija primerna tako ( ) 2 1 2 1 2 σ − − = clC , e x x x x , � � Gradbeni vestnik • letnik 65 • maj 2016 111 METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ•Teja Melink, Jože Korelc za stohastične procese, ki se spreminjajo preko prostorske domene problema, kot za stohastične procese, ki se spreminjajo skozi čas, ter za poljubno korelirane stohastične procese. Stohastično polje označimo z w(x,θ ), x je krajevni vektor, definiran preko območja pro- storske domene, in q je dogodek iz prostora slučajnih dogodkov. Z razvojem stohastičnega polja v K-L dekompozicijo se stohastično polje zapiše kot vsoto pričakovane vrednosti tega polja w (x) in neskončne vrste, ki se jo v numeričnih simulacijah aproksimira z M prvimi členi (a) (b) (c) (d) Slika 2•Stohastična polja z različno korelacijo: (a) nizko, (b) srednje, (c) visoko in (d) neskončno Slika 2•korelirano stohastično polje Slika 3•Eksponentna kovariančna funkcija za tri različne korelacijske dolžine ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 θ λ ξ θ λ ξ θ ∞ = = = + ≈ +∑ ∑ M k k k k k k k k w , w f w fx x x x x , � � enačbe 2. vrste. Ker za splošen primer analitična rešitev te enačbe ne obstaja, smo za numerično reševanje Fredholmove enačbe uporabili Galerkinovo metodo, ki reševanje te enačbe prevede na reševanje posplošenega problema lastnih vrednosti. Za ta korak smo razvili stohastične končne elemente in kjer so ξk (θ) normirane nekorelirane slučajne spremenljivke (ki so v primeru Gaussovega slučajnega polja Gaussove spremenljivke z ničelno pričakovano vrednostjo in enot- sko varianco), λk in ƒk (x) so lastne vred- nosti in lastne funkcije, ki jih dobimo kot rešitev homogene Fredholmove integralske (2) Gradbeni vestnik • letnik 65 • maj 2016112 Teja Melink, Jože Korelc•METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ 3•IZRAČUN STATISTIKE ODZIVA KONSTRUKCIJE Odziv konstrukcije označimo z z(ξk) in pred- stavlja poljubno količino, na primer pomik v določeni točki, maksimalno napetost, plastični del deformacij v določeni točki, mejno nosil- nost itd. Odziv z je odvisen od slučajnih spremenljivk ξk. V okviru izračuna statistike odziva konstrukcije so nas zanimale predvsem statistike nižjega reda, in sicer pričakovana vrednost in standardna deviacija odziva. Za izračun odziva stohastičnih problemov smo uporabili perturbacijsko metodo in simulacije Monte Carlo. Metoda Monte Carlo je bila v raziskavah izbrana predvsem z namenom preveriti rezul- tate razvite formulacije perturbacijske metode. Metoda Monte Carlo je zaradi svoje splošnosti primerna za poljubne stohastične procese, njena implementacija pa je dokaj trivialna. Pri tej metodi se vrednost stohastičnih spremen- ljivk ξk v izrazu (2) 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 problem. Izvede se želeno število simulacij in v vsaki simulaciji beleži odziv konstrukcije. Statistiko odziva se nato določi s standard- nimi statističnimi metodami. Posledica tega je, da je natančnost rešitve odvisna od števila simulacij, kar omejuje praktično uporabnost metode Monte Carlo predvsem pri nelinearnih problemih večjih dimenzij ali ko nas zanimajo repne vrednosti odziva (npr. pri analizi zanes- ljivosti). Posvetili smo se predvsem razvoju numerično učinkovite formulacije perturbacijske metode poljubnega reda, ki za določitev statistike odziva zahteva samo eno direktno simulacijo odziva konstrukcije. Posledično je računski čas, potreben za izvedbo stohastične ana- lize, neprimerno krajši v primerjavi z drugimi metodami za izračun statistik. Perturbacijska metoda odziv konstrukcije aproksimira z raz- vojem odziva v Taylorjevo vrsto kjer je si red razvoja v vrsto za i -ti pa- rameter in 0ξi pričakovana vrednost i -tega slučajnega parametra. V primeru Gausso- vega polja je 0ξi = 0. Red perturbacijske me- tode pomeni red Taylorjeve vrste, v katero se razvije odziv. Koeficienti Taylorjevo vrste v enačbi (3) so odvodi funkcije odziva Ko je razvoj odziva v Taylorjevo vrsto izračunan, sledi izračun statistik odziva Za rešitev tega problema smo kombinirali numerični in simbolni pristop. V uporablje- nem simbolno-numeričnem pristopu smo kombinirali: program za simbolno računanje Mathematica, avtomatski generator pro- gramske kode AceGen [Korelc, 2009b] in okolje za numerično modeliranje po metodi končnih elementov AceFEM [Korelc, 2009b]. V AceGenu so vgrajene različne tehnike za optimiranja programske kode in teh- nika avtomatskega odvajanja, ki omogoča izpeljavo poljubnih odvodov v zaključeni ob- liki. Dodatno smo formulacijo občutljivostne analize prvega reda z uporabo avtomat- skega odvajanja [Korelc, 2009a] razširili za občutljivostno analizo (poljubnega) višjega reda. Z razvito formulacijo perturbacijske metode je za izračun odziva v izrazu (3) in statistik odziva potrebna ena sama izvedba mehanske analize. Rezultat vseh naštetih uporabljenih tehnik in formulacij je numerično zelo učinkovit postopek, posledično je možna izvedba perturbacijske metode višjega reda za neprimerno višje število uporabljenih členov K-L dekompozicije, kot je možno z drugimi metodami za izračun odvodov (npr. metodo končnih diferenc). V zvezi z uporabo perturbacijske metode je treba opozoriti še na njeno omejitev, in sicer je metoda primerna za opis odzivov konstrukcij, ki se v okolici pričakovane vred- nosti odziva obnašajo zvezno. V mehaniki konstrukcij to pomeni, da je uporaba per- turbacijske metode omejena na primere, pri katerih se porušni mehanizem v odvisnosti od slučajnih parametrov ne spreminja. ( ) ( ) ( ) ( ) 1 2 1 2 1 2 1 2 0 1 2 ... 0 0 0 1 2 1 1 2 2 0 0 0 ,1 2 1 2 1,..., 1, ,..., ... ... , ! !... ! ... M M M M M i i s s s l l l l l l M M Ml l l l l l M M i M zz l l l ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + + + = = = = = ∂ = ⋅ − − − ∂ ∂ ∂ ∑∑ ∑ � � 1 2 1 2 0 ... ,1 2 1,..., ... M M i i l l l l l l M i M z ξ ξξ ξ ξ + + + = = ∂ ∂ ∂ ∂ Te se izračuna z enostavnimi matematičnimi enačbami (glej npr. [Sudret, 2002]). Natančnost perturbacijske metode je odvisna od nelinearnosti odziva in reda perturbacij- ske metode. V splošnem velja: višji ko je red perturbacijske metode, bolj natančen je opis odziva v njegovi širši okolici. Kljub temu je perturbacijska metoda višjega reda v litera- turi zelo redko uporabljena. Vzrok je pomanj- kanje numerično učinkovitih in natančnih algoritmov za izračun odvodov višjega reda funkcije odziva, ki so potrebni za izvedbo perturbacijske metode višjega reda. Aproksi- macija odvodov s končnimi diferencami se izkaže v primeru večjega števila stohastičnih spremenljivk za drago metodo, v primeru višjih odvodov pa tudi za nezanesljivo me- todo. Alternativa je analitična občutljivostna analiza, pri kateri poljubne odvode dobimo na podlagi ene same direktne simulacije odziva pri pričakovani vrednosti stohastičnih parametrov [Keulen, 2005]. Dobljeni odvodi so natančni v okviru računske natančnosti računalnika. Vendar pa taka analiza zahteva izpeljavo višjih odvodov ravnotežnih enačb problema v zaključeni obliki, kar se izkaže v primeru splošnih materialno in geometrij- sko nelinearnih problemov za zelo zahteven problem. ( ) ( ) ( ) ( ) 1 2 1 2 1 2 1 2 0 1 2 ... 0 0 0 1 2 1 1 2 2 0 0 0 ,1 2 1 2 1,..., 1, ,..., ... ... , ! !... ! ... M M M M M i i s s s l l l l l l M M Ml l l l l l M M i M zz l l l ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ + + + = = = = = ∂ = ⋅ − − − ∂ ∂ ∂ ∑∑ ∑ � � 1 2 1 2 0 ... ,1 2 1,..., ... M M i i l l l l l l M i M z ξ ξξ ξ ξ + + + = = ∂ ∂ ∂ ∂ . (3) predlagali nekaj izboljšav (npr. modifikacijo kovariančne funkcije), ki občutno skrajšajo računski čas in hkrati izboljšajo robust- nost metode. Nekonsistentno modificirana kovariančna funkcija namreč lahko privede do numerične nestabilnosti Galerkinove me- tode [Melink, 2014a]. Ko izračunamo deterministične člene λk in ƒk (x) v izrazu (2) (slučajne spremenljivke ξk ostanejo neznanke problema), je prvi ko- rak metode stohastičnih končnih elementov zaključen. Gradbeni vestnik • letnik 65 • maj 2016 113 METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ•Teja Melink, Jože Korelc Natančnost in numerično učinkovitost raz- vite avtomatizacije metode stohastičnih končnih elementov smo preverili na različnih računskih primerih, ki smo jih izbrali tako, da pokrivajo vse osnovne primere, pomembne za inženirsko prakso. V računskih primerih so bili zajeti različna mejna stanja in porušni meha- nizmi, in sicer statična analiza pri dani obtežbi, mejno stanje materiala in mejno stanje stabil- nosti. V nadaljevanju bomo zaradi omejenosti prostora predstavili le nekaj primerov uporabe predstavljene metode stohastičnih končnih elementov. 4.1 Prikaz vpliva upoštevanja korelacije 4.1 stohastičnega polja korozije Korozija je kompleksen pojav razpadanja ko- vine, ki so mu bile v preteklosti namenjene že 4•PRIKAZ UPORABE RAZVITEGA PRISTOPA K STOHASTIČNI ANALIZI 4•KONSTRUKCIJ hwtw tf hf PP xz z y Slika 4•Primer tlačenega prostoležečega nosilca Slika 5•Histogram kritične uklonske sile za nosilec dolžine (a) 10 m, (b) 3 m in (c) 5 m (a) (b) (C) Gradbeni vestnik • letnik 65 • maj 2016114 Teja Melink, Jože Korelc•METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ številne študije (npr. [Valor, 2007], [Tarant- seva, 2010]). Zaključki vseh obravnavanih študij in raziskav so, da je korozija primerna za opis s stohastičnim poljem. Kljub temu v literaturi skorajda ni zaslediti študije, ki bi jo obravnavala na ta način. V izbranem računskem primeru smo pokazali, da sta raz- vita avtomatizacija za opis stohastičnih polj in izračun statistike odziva z uporabo metode končnih elementov primerna tudi za opis procesa korozije ter da upoštevanje korelacije stohastičnega polja (tj., da se sprememba debeline pasnic in stojine zaradi procesa korozije preko nosilca naključno spreminja) pomembno vpliva na odziv konstrukcije. Obravnavan je jekleni prostoležeči nosi- lec s prečnim prerezom I-profila z dimen- zijami stojine hw /tw = 276/7 mm in pasnice hf /tf = 300/12 mm, kot prikazuje slika 4. Za preprečitev lokalnih deformacij na mestu vnosa koncentrirane tlačne sile je nosilec na obeh straneh ojačan s čelno pločevino. Dolžine nosilcev so 10 m (v rezultatih oznaka TP1), 5 m (TP2) in 3 m (TP3). Nosilec je obre- menjen s tlačno osno silo P. Pasnici in stojino nosilca smo modelirali z 9-vozliščnimi lupinastimi hiperelastičnimi končnimi elementi. Za opis procesa korozije smo privzeli model jamičaste korozije. Proces korozije smo opisali s srednje koreliranim stohastičnim poljem spremembe debeline pločevine (korelacijska dolžina lc je 300 mm). Stohastično polje spremembe debeline ima pričakovano vrednost µ = -1 mm in stan- dardno deviacijo σ = 0,3 mm. Upoštevanih je 40 členov K-L vrste. Poleg primerov z upoštevanjem korelacije stohastičnega pol- ja je bil za vsak testni primer narejen še primer z enakomerno stanjšano debelino, kjer se je sprememba debeline v vsaki simu- laciji izračunala po Gaussovi porazdelitvi, kar pomeni, da je stohastično polje neskončno korelirano. Ti testni primeri imajo poleg oznake TP še zvezdico (TP1*, TP2*, TP3*). Za vsak testni primer je bilo izvedenih 10.000 simulacij Monte Carlo, v katerih se je z linearizirano sta- bilnostno analizo izračunala kritična uklonska sila Pcrit. V preglednici 1 je primerjava dobljenih rezul- tatov. Izkaže se, da v primeru 10 m nosilca kritična sila povzroči globalni uklon, pri 3 m nosilcu pa pride v vseh primerih do lokalnega izbočenja pločevin, zato sta histograma Pcrit za ti dve dolžini nosilca simetrični (slika 5 (a) in (b)). V obeh primerih se razpršenost rezultatov za primer z upoštevanjem korelacije zelo razlikuje od primera, ko se privzame, da se debelina pasnic in stojine nosilca stanjša konstantno po celotnem nosilcu. Vpliv korela- cije stohastičnega polja v primeru globalnega uklona nima bistvenega vpliva na pričakovano vrednost kritične sile, ta v obeh primerih znaša 1024 kN, izkaže pa se, da ima v primeru, ko je kritično lokalno izbočenje pločevin, korelacija pomemben vpliv na pričakovano kritično silo, ki v tem primeru znaša 4185 kN, medtem ko znaša pri simulacijah z enakomerno stanjšano debelino 4315 kN. V primeru 5 m nosilca se izkaže, da v posameznih realizacijah stohastičnega polja nastopi globalna ali lokalna izguba stabil- nosti. V simulacijah Monte Carlo, v katerih pride na neki dovolj veliki površini stojine ali pasnice nosilca do večjega stanjšanja debe- line, je kritična sila, ki je potrebna za lokalno izbočenje pločevine, manjša od kritične sile, ki bi povzročila globalni uklon nosilca. Histo- gram Pcrit je zaradi tega nesimetričen (slika 5 (c)). Asimetričnost je v primeru upoštevanja korelacije stohastičnega polja -2,01 kN in je v absolutnem smislu znatno večja kot v primeru modela nosilca z enakomerno stanjšanimi pločevinami, kjer znaša -0,90 kN. Temu primer- no je pričakovana vrednost Pcrit manjša v prvem primeru, in sicer 4032 kN, medtem ko je za TP2* 4041 kN. Za primerjavo: kritična sila 5 m nosilca, ki se mu debelina stojine in pasnic enakomerno stanjša za pričakovano vrednost, je 4051 kN in povzroči globalni uklon nosilca. Prikazani računski primer nam pokaže, da neupoštevanje korelacije stohastičnega polja lahko vodi do precej zgrešene ocene dejanske- ga odziva konstrukcije. Pri oceni spodnjih mejnih vrednosti kritične sile je ta ocena v vseh primerih konservativna, medtem ko je ocena pričakovane vrednosti ob neupoštevanju korelacije lahko podcenjena ali precenjena, saj s tako poeno- stavljenim računskim modelom ni mogoče za- jeti vseh možnih odzivov konstrukcije. 4.2 Nelinearna stohastična analiza 4.2 korodiranega tlačeno obremenjenega 4.2 nosilca Tudi v tem primeru smo obravnavali korodirani tlačno obremenjeni jekleni nosilec, vendar smo tokrat upoštevali tudi začetno geome- trijsko nepopolnost nosilca ter odziv izračunali z nelinearno elastoplastično analizo. Nosilec ima prečni prerez I-profila IPE300. Dolžina nosilca je 5,25 m. Primerjalna vitkost nosilca, ki je enakomerno korodiran, je 2,0, zato je ne glede na stopnjo korodiranosti pričakovani porušni mehanizem globalni uklon. Pasnici in stojino nosilca smo modelirali z lupinastimi 9-vozliščnimi elastoplastičnimi končnimi ele- menti, ki upoštevajo velike deformacije. Proces korozije smo modelirali s stohastičnim poljem spremembe debeline lupinastih končnih elementov, v K-L dekompoziciji smo ohranili 66 členov. V okviru perturbacijske metode drugega reda to pomeni, da je treba izračunati 2277 odvodov prvega in drugega reda mejnega obtežnega faktorja po slučajnih parametrih. Nosilec je bil na obeh koncih obremenjen s tlačno osno silo P. Zanimala nas je nosilnost konstrukcije v odvisnosti od slučajnih spremenljivk, s katerimi smo opisali proces korozije. Kritično silo, pri kateri se kon- strukcija ukloni, smo označili s Pcrit. ADB formulirano perturbacijsko metodo drugega reda smo primerjali z metodo Mon- te Carlo, pri kateri je bilo izvedenih 5000 simulacij. Na sliki 6 sta prikazana naključna simulacija procesa korozije in deformiran (za boljšo predstavo so deformacije na sliki 10-kratno povečane) nosilec v trenutku, ko je bila dosežena mejna obtežba za dane vhodne podatke in prikazano porazdelitev korozije. TP1 TP1* TP2 TP2* TP3 TP3* Minimalna vrednost 996 929 3480 3258 3497 3232 Maksimalna vrednost 1119 1134 4183 4494 5783 5497 Pričakovana vrednost 1024 1024 4032 4041 4185 4315 Stand. deviacija 7 28 71 132 176 309 0,05 kvantil 1013 977 3885 3802 3893 3817 Preglednica 1•Karakteristične vrednosti simulacij kritične sile za testne primere (v kN) Gradbeni vestnik • letnik 65 • maj 2016 115 METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ•Teja Melink, Jože Korelc V preglednici 2 je prikazna primerjava rezul- tatov in računskih časov ADB formulirane perturbacijske metode in simulacij Monte Carlo. Rezultati obeh analiz se ujemajo. Primerjava računskih časov analiz pokaže, da je ADB formulirana perturbacijska metoda za izračun rezultatov potrebovala približno 25-krat manj časa v primerjavi s simulacijami Monte Carlo. 4.3 Geometrijska nepopolnost sinusoidnega 4.3 panela Tretji računski primer obravnava upo- gib obojestransko vpetega sinusoidnega strešnega panela s stohastično spremenljivo amplitudo valov. Panel je na zgornji in spodnji strani obdan z jekleno pločevino, na sredini pa zapolnjen z izolacijsko peno. Geometrij- sko idealen panel ima vzdolž abscisne osi obliko sinusoidne funkcije. Oblika in dimen- zije panela in vsi parametri, ki so potrebni za analizo, so prikazani na sliki 7. Zanima nas vertikalni pomik v sredini panela v odvis- nosti od stohastičnega polja geometrijske nepopolnosti. Za modeliranje panela smo izbrali dvodimen- zionalne 4-vozliščne hiperelastične končne elemente. Izolacijska pena v sredini panela je diskretizirana z mrežo 120 × 4 končnih ele- mentov, za vsako izmed pločevin na obodu panela pa je uporabljena diskretizacijska mreža 120 × 2 končnih elementov. Geometrijska nepopolnost panela je mo- delirana z enodimenzionalnim stohastičnim poljem amplitude valov hval (definiranim vzdolž osi x). Privzeta je relativno visoka ko- reliranost stohastičnega polja (lc = 400 cm), zato za ustrezno natančnost izračuna za- dostujejo 4 členi K-L vrste. Odziv panela je izračunan s perturbacijsko metodo četrtega reda. Rezultat perturbacijske metode četrtega reda je razvoj odziva panela (vm) v vrsto po enačbi (3), v kateri se upošteva le člene, katerih vsota potenc slučajnih spre- menljivk (l1 + l2 + l3 + l4) je manjša ali enaka 4. Skupno je treba izračunati 69 odvodov (∂l1+ l2+ l3+ l4vm / ∂ξl11 ∂ξl22 ∂ξl33 ∂ξl44). Odvode smo izračunali z ADB formulirano občutljivostno analizo in standardno metodo končnih dife- renc. Rezultati obeh metod so se popol- noma ujemali. Primerjava računskih časov, potrebnih za izvedbo celotne stohastične analize z različnimi metodami, je prikazana v preglednici 3 in je močno v korist razviti ADB formulirani perturbacijski metodi. Za primer- javo: računski čas za izvedbo ene direktne analize problema je 0,11 s. Vsi izračuni so bili opravljeni v računalniku z Intelovim proce- Slika 6•Naključno izbrana realizacija stohastičnega polja korozije in deformiran nosilec v trenutku, Slika 6•ko je dosežen mejni obtežni faktor (10-kratna povečava deformacij) Perturbacijska metoda (ADB-formulacija) Monte Carlo (5000 simulacij) γmax 359.16 359.14 σγmax 1.38 1.39 t 1 h, 40 min. 41 h, 53 min. Preglednica 2•Pričakovana vrednost (γmax), standardna deviacija (σγmax) in računski čas (t ), potreben Preglednica 2•za izvedbo celotne analize Slika 7•Sinusoidna panelna plošča sorjem i7-950 (3.07 GHz) z 12 GB RAM-a, 4 jedri in 8 nitmi. Naredili smo tudi primerjavo ocene odziva vm z uporabo perturbacijske metode in simulaci- jami Monte Carlo, in sicer dvema serijama z različnim številom simulacij, 1000 in 50.000. Čeprav deluje 1000 simulacij na prvi pogled dokaj veliko število ponovitev, se je izkazalo, da je bilo za konvergenco metode potrebno večje število simulacij, drugače je bil nabor naključno izbranih slučajnih spremenljivk premajhen, da bi predstavljal statistično tipično množico. Rezultati so prikazani v preglednici 3. Gradbeni vestnik • letnik 65 • maj 2016116 Teja Melink, Jože Korelc•METODA STOHASTIČNIH KONČNIH ELEMENTOV V MODELIRANJU KONSTRUKCIJ Odziv vm je nelinearen, zato se na prika- zanem primeru lepo vidi tudi vpliv reda perturbacijske metode na natančnost rezul- tatov. Perturbacijska metoda drugega in višjih redov se z dejanskimi simulacijami veliko bolje ujemajo v primerjavi s perturba- cijsko metodo prvega reda, kar je prikazano na sliki 8. Perturbacijska metoda Metoda Monte Carlo ADB formulirana občutljivostna analiza Metoda končnih diferenc 103 sim. (5 ·104 sim.) νm -4,126 -4,128 -4,126 σνm 0,0586 0,0571 0,0588 Skupni čas 5,27 s 173,28 s 112 s 5604 s Cilj razvoja avtomatizacije metode stoha- stičnih končnih elementov je bil, da bi v analizi mehanskih problemov z upoštevanjem naključnosti vhodnih parametrov dosegli računske čase, primerljive determinističnemu pristopu, ter da bi bila metoda numerično sta- bilna in zanesljiva. Z uporabo kombiniranega Preglednica 3•Primerjava statistik odziva (pričakovane vrednosti νm ter standardne deviacije σνm, Preglednica 3•oboje v cm) in primerjava računskih časov različnih metod Slika 8•Vertikalni pomik na sredini panela (νm), izračunan s perturbacijsko metodo različnih redov Slika 8•in simulacijami Monte Carlo 5•SKLEPI simbolnega in numeričnega pristopa smo ta cilj dosegli. Metode, ki so uporabljene v postopku, smo izbirali tako, da so kom- patibilne z metodo končnih elementov, ki je v inženirski praksi najpogosteje uporabljena metoda za mehansko analizo konstrukcij. Poleg tega smo metode izbrali tako, da so robustne in omogočajo tudi njihovo razširitev na druge tipe stohastičnih problemov (npr. na stohastična polja z Negaussovimi robnimi porazdelitvami, na stohastične procese, kjer se naključnost spreminja s časom, itd.) Zato predstavlja prikazani stohastični pristop k reševanju mehanskih problemov orodje za poglobljene študije mehanskih problemov z upoštevanjem njihove stohastične narave ter potencial za uporabo v splošni inženirski praksi in za nadaljnji razvoj Evrokodov. Ghanem, R., Spanos, P. D., Stochastic Finite Element: A Spectral Approach, Dover, New York, 2003. Keulen, F., Haftka, R. T., Kim, N. H., Review of options for structural design sensitivity analysis, Part 1: Linear systems, Comput. Methods Appl. Mech. Engrg. 194, 3213–3243, 2005. Korelc, J., Automation of primal and sensitivity analysis of transient coupled problems, Comput. Mech, 44, 631–649, 2009a. Korelc, J., AceGen and AceFEM User Manual, www.fgg.uni-lj.si/symech/, 2009b. Melink,T., Korelc, J., Stability of Karhunen-Loève expansion for the simulation of Gaussian stochastic fields using Galerkin scheme, Probabilistic Engineering Mechanics, 37, str. 7–15, 2014a. Melink, T., Metoda stohastičnih končnih elementov v modeliranju konstrukcij, doktorska disertacija, Ljubljana, Univerza v Ljubljani, 191 str., 2014b. Mathematica 9.0, Wolfram Research Inc., povzeto po: http://www.wolfram.com, 2012. Sudret, B., Comparison of the Spectral Stochastic Finite Element method with the Perturbation method for second moment analysis, Proc. 1st ASRANET Int. Coll., 2002. Tarantseva, K.R., Models and Methods of Forecasting Pitting Corrosion, Protection of Metals and Physical Chemistry of Surfaces, Vol. 46, 139–147, 2010. Valor, A., Caleyo, F., Alfonso, L., Rivas, D., Hallen, J. M., Stochastic modeling of pitting corrosion: a new model for initiation and growth of multiple corrosion pits, Corrosion Science, 49, p. 559–579, 2007. 6•LITERATURA