Monolitick´e a iteraˇcn´ı ˇreˇsen´ı Biotovy poroelasticity

(1)

Monolitick´ e a iteraˇ cn´ı ˇ reˇ sen´ı Biotovy poroelasticity

Bakal´ aˇ rsk´ a pr´ ace

Studijn´ı program: B3901 – Aplikované vˇedy v inˇzenýrstv´ı Studijn´ı obor: 3901R055 – Aplikované vˇedy v inˇzenýrstv´ı Autor práce: Sabina Bednáˇrová

Vedouc´ı pr´ace: Mgr. Jan Stebel, Ph.D.

(2)

Monolithic and iterative solution of Biot’s poroelasticity

Bachelor thesis

Study programme: B3901 – Applied Sciences in Engineering Study branch: 3901R055 – Applied Sciences in Engineering Author: Sabina Bedn´aˇrov´a

Supervisor: Mgr. Jan Stebel, Ph.D.

(3)

Zadání bakalářské práce

Monolitické a iterační řešení Biotovy poroelasticity

Jméno a příjmení: Sabina Bednářová Osobní číslo: M16000097

Studijní program: B3901 Aplikované vědy v inženýrství Studijní obor: Aplikované vědy v inženýrství

Zadávající katedra: Ústav nových technologií a aplikované informatiky Akademický rok: 2018/2019

Zásady pro vypracování:

1. Seznamte se:

(i) s Biotovou teorií poroelasticity a jejím využitím pro modelování hydrogeologických procesů, (ii) s iteračními schématy a dimenzionální redukcí pro Biotův systém,

(iii) s výpočetní knihovnou FEniCS a základními metodami pro řešení soustav lineárních rovnic.

2. Aplikujte iterační schéma ”fixed-stress splitting” na redukovaný Biotův systém. Implementujte toto schéma do existujícího výpočetního kódu pro redukovaný Biotův systém ve 2D. Proveďte výpočetní analýzu

konvergence.

3. Na základě předchozích zkušeností navrhněte vhodnou volbu relaxačních parametrů pro oblast a puklinu v závislosti na parametrech modelu, zajišťující konvergenci iteračního schématu.

4. Porovnejte výpočetní náročnost monolitického a iteračního řešení redukovaného Biotova systému v závislosti na počtu elementů sítě a fyzikálních parametrech.

(4)

Rozsah grafických prací: dle potřeby Rozsah pracovní zprávy: 30 – 40 stran Forma zpracování práce: tištěná/elektronická

Seznam odborné literatury:

[1] CHENG, A. H.-D. Poroelasticity. Springer, 2016.

[2] BOTH, J.W., BORREGALES, M., NORDBOTTEN, J.M., KUMAR, K., RADU, F.A. Robust fixed stress splitting for Biot’s equations in heterogeneous media, Appl Math Lett 68:101-108, 2017.

[3] ANGOT, P., BOYER, F., HUBERT, F. Asymptotic and numerical modelling of flows in fractured porous media, ESAIM:M2AN 43:239-275, 2009.

[4] TEBBENS, J.D., HNĚTYNKOVÁ, I., PLEŠINGER, M., STRAKOŠ, Z., TICHÝ, P. Analýza metod pro maticové výpočty Základní metody. MatfyzPress, 2012.

[5] FEniCS Project [online, cit. 24. 9. 2018]. Dostupné z https://fenicsproject.org.

Vedoucí práce: Mgr. Jan Stebel, Ph.D.

Ústav nových technologií a aplikované informatiky Datum zadání práce: 18. října 2018

Předpokládaný termín odevzdání: 30. dubna 2019

L. S.

prof. Ing. Zdeněk Plíva, Ph.D.

děkan

Ing. Josef Novák, Ph.D.

vedoucí ústavu

(5)

Prohlášení

Byla jsem seznámena s tím, že na mou bakalářskou práci se plně vzta- huje zákon č. 121/2000 Sb., o právu autorském, zejména § 60 – školní dílo.

Beru na vědomí, že Technická univerzita v Liberci (TUL) nezasahuje do mých autorských práv užitím mé bakalářské práce pro vnitřní potřebu TUL.

Užiji-li bakalářskou práci nebo poskytnu-li licenci k jejímu využití, jsem si vědoma povinnosti informovat o této skutečnosti TUL; v tomto pří- padě má TUL právo ode mne požadovat úhradu nákladů, které vyna- ložila na vytvoření díla, až do jejich skutečné výše.

Bakalářskou práci jsem vypracovala samostatně s použitím uvedené literatury a na základě konzultací s vedoucím mé bakalářské práce a konzultantem.

Současně čestně prohlašuji, že texty tištěné verze práce a elektronické verze práce vložené do IS STAG se shodují.

30. 4. 2019 Sabina Bednářová

(6)

Podˇ ekov´ an´ı

V prvn´ı ˇradˇe bych ráda podˇekovala vedouc´ımu své práce, panu Mgr. Janu Steblovi, Ph.D., za neskuteˇcnou trpˇelivost, ochotu a velké mnoˇzstv´ı ˇcasu, které mi vˇenoval, za pomoc s moj´ı ˇceˇstinou a skvˇelé veden´ı, které mˇe donutilo na práci intenzivnˇe pracovat po celý ˇskoln´ı rok.

Dále bych ráda podˇekovala doc. Ing. Janu ˇSemberovi, Ph.D., za ochotnou diskuzi po vyuˇcovac´ıch hodinách a za dodán´ı kuráˇze ve chv´ıl´ıch, kdy mi bylo sp´ıˇse do breku.

Dˇekuji tak´e svoj´ı rodinˇe, bez jej´ıˇz podpory bych se nikdy takhle daleko nedostala, hlavnˇe sv´ym rodiˇc˚um, ˇze mi ˇsli v tomhle smˇeru pˇr´ıkladem.

(7)

Abstrakt

Tato bakaláˇrská práce se zabývá numerickým ˇreˇsen´ım Biotova systému v prostˇred´ı s puklinami. Je odvozeno iteraˇcn´ı schéma typu fixed-stress splitting a popsána jeho implementace s vyuˇzit´ım knihovny FEniCS. Dále je provedena výpoˇcetn´ı analýza konvergence v závislosti na relaxaˇcn´ıch parametrech a na fyzikáln´ıch parametrech. Na základˇe pˇredchoz´ıch zkuˇsenost´ı je narˇzena volba re- laxaˇcn´ıch parametr˚u pro oblast a puklinu, zajiˇst’uj´ıc´ı konvergenci iteraˇcn´ıho schématu. Nakonec je porovnána výpoˇcetn´ı nároˇcnost monolitického a iteraˇcn´ıho ˇreˇsen´ı redukovaného Biotova systému v závislosti na poˇctu element˚u s´ıtˇe a fyzikáln´ıch parametrech.

Kl´ıˇcová slova: poroelasticita, Biot˚uv redukovaný model, fixed- stress splitting, monolitické ˇreˇsen´ı, iteraˇcn´ı ˇreˇsen´ı

Abstract

This bachelor thesis deals with numerical solution of Biot’s sys- tem in domains containing fractures. An iterative fixed-stress type splitting is derived and its implementation using FEniCS library is described. Further we perform computational convergence analysis according to the dependency on tuning parameters and physical parameters. Based on previous experience, we suggest a convenient choice of tuning parameters for the domain and the fracture, en- suring convergence of iterative scheme and compare computational costs of monolithic and iterative solution depending on the quantity of mesh elements and physical parameters.

Keywords: poroelasticity, Biot’s reduced model, fixed-stress splitting, monolithic solution, iterative solution

(8)

Obsah

1 Uvod´ 8

1.1 Porézn´ı materiály . . . 8 1.2 Modelován´ı hydraulického ˇstˇepen´ı . . . 9 1.3 Biot˚uv model a jeho varianta pro prostˇred´ı s puklinami . . . 10 2 Pˇrehled numerických metod pro klasický Biot˚uv model 12 2.1 Monolitické (implicitn´ı) ˇreˇsen´ı - pˇr´ımé a iteraˇcn´ı metody . . . 12 2.2 Iteraˇcn´ı (sekvenˇcnˇe implicitn´ı) ˇreˇsen´ı . . . 14 3 Numerické ˇreˇsen´ı Biotova modelu pro oblast s puklinou 16 3.1 Monolitické ˇreˇsen´ı. . . 21 3.2 Aplikace fixed stress splitting . . . 23 3.3 Analýza konvergence . . . 28

4 Z´avˇer 39

(9)

1 Uvod ´

V dneˇsn´ı dobˇe, kdy ˇreˇs´ıme spoustu environmentáln´ıch problém˚u, se snaˇz´ıme zapojo- vat v´ıce ekologická ˇreˇsen´ı, jako je napˇr´ıklad obnovitelná energie. Obnovitelná energie je taková, která se v lidském ˇcasovém mˇeˇr´ıtku pˇrirozenˇe obnovuje, na rozd´ıl od neob- novitelné energie, jako jsou tˇreba fosiln´ı paliva, které se v lidském ˇcasovém mˇeˇr´ıtku neobnovuj´ı a tak m˚uˇze doj´ıt k jejich vyˇcerpán´ı, nebo náklady na tˇeˇzen´ı tˇechto paliv budou tak vysoké, ˇze je uˇz nebude ekonomicky v˚ubec výhodné pouˇz´ıvat. Mimo to vˇetˇsina neobnovitelné energie negativnˇe pˇrisp´ıvá ke zmˇenˇe klimatu.

Jedn´ım z typ˚u obnovitelné energie je geotermáln´ı energie, kdy vyuˇz´ıváme tepelný projev zemského jádra. Toto vyuˇzit´ı geotermáln´ı energie je ne vˇzdy obnovitelné - nˇekteré zdroje geotermáln´ı energie lze vyˇcerpat v horizontu des´ıtek let, a velkou nevýhodou je také velká závislost na geografickém um´ıstˇen´ı.

V systému Hot Dry Rock (HDR) jde o vyuˇzit´ı akumulovaného tepla v hlu- binné zónˇe s vysokou propustnost´ı. Takovou zónu je obvykle nutné umˇele vytvoˇrit, napˇr´ıklad s vyuˇzit´ım existuj´ıc´ı s´ıtˇe puklin. Pomoc´ı Biotovy teorie poroelasticity chceme sledovat, studovat proces otev´ırán´ı pro vyuˇzit´ı geotermáln´ı energie.

Biotova poroelasticita je soustava rovnic, která popisuje interakci mezi proudˇen´ım a deformac´ı porézn´ıho média. Z takové interakce m˚uˇzeme zjistit, jak se bude puklina v tomto médiu otev´ırat ˇci zav´ırat, pokud do n´ı budeme vtláˇcet kapalinu pod velmi vysokým tlakem.

C´ılem této bakaláˇrské práce je seznámen´ı se s ˇreˇsen´ım Biotových rovnic poroelasticity a porovnán´ı nároˇcnosti dvou metod ˇreˇsen´ı: monolitické a iteraˇcn´ı. K porovnán´ı tˇechto dvou ˇreˇsen´ı je pouˇzita výpoˇcetn´ı knihovna FeniCS.

1.1 Por´ ezn´ı materi´ aly

Celý náˇs okoln´ı svˇet je obklopen porézn´ımi materiály, nˇekteré z nich jsou pˇr´ırodn´ı, jiné jsou umˇele vyrobené ˇclovˇekem. Existuje jich mnoho v r˚uzných formách - patˇr´ı sem granulované maateriály (napˇr. p´ısek), kde kaˇzdá ˇcástice kˇremenu má jinou

(10)

velikost a ˇcástice nejsou pevnˇe spojené, proto jimi m˚uˇze tekutina volnˇe protékat bez vˇetˇs´ıho odporu - toto jsou takzvané granulované materiály. Naopak napˇr´ıklad p´ıskovec, coˇz jsou ˇcástice p´ısku pevnˇe spojené kalcitem a j´ılem, je porézn´ı pevnou látkou. Pˇrestoˇze jsou ˇcástice p´ıskovce k sobˇe pevnˇe spojeny a tekutina nem˚uˇze skrz p´ıskovec protékat tak volnˇe jako hromadou p´ısku, m˚uˇze být obecnˇe porézn´ı pevná látka dobrým vodivým médiem. V kaˇzdém porézn´ım médiu se pˇrirozenˇe vyskytuj´ı póry, puklinky nebo trhliny, v nichˇz se m˚uˇze kapalina ˇs´ıˇrit. Velikost pór˚u se obecnˇe v kaˇzdém porézn´ım médiu liˇs´ı, v závislosti na tom, jak vznikalo. Napˇr´ıklad vul- kanická hornina má póry o velikosti aˇz jednoho centimetru - pˇri chlazen´ı lávy se vytvoˇrily plynové bublinky z rozpouˇstˇen´ı tˇekavých látek [1]. Dalˇs´ı porézn´ı horninou je napˇr´ıklad ˇzula, která je krystalickou horninou a pro HDR systém je ˇzulový masiv vhodný. Schopnost porézn´ıho materiálu propouˇstˇet vodu se nazývá permeabilita (propustnost).

Vˇetˇsina pˇr´ırodnˇe vyskytuj´ıc´ıch se materiál˚u je porézn´ıch. I lidské tˇelo je porézn´ım materiálem. Mezi umˇele vytvoˇrené porézn´ı materiály patˇr´ı napˇr´ıklad beton (tato látka má relativnˇe malou propustnost) nebo r˚uzné stavaˇrské pˇeny.

1.2 Modelov´ an´ı hydraulick´ eho ˇ stˇ epen´ı

Hydraulické ˇstˇepen´ı je metoda pro tvorbu puklin v horninovém masivu. Tato metoda nen´ı pro náˇs svˇet niˇc´ım novým; prvn´ı aplikace probˇehly jiˇz v roce 1947 [2]. Jde o vytvoˇren´ı hlubinného vrtu, do kterého se poté vtláˇc´ı kapalina pod velmi vysokým tlakem. Pokud tlak kapaliny pˇrekoná pevnost horniny, dojde k jej´ımu rozˇstˇepen´ı a tvoˇren´ı puklin, kterými se poté kapalina dále ˇs´ıˇr´ı (nejˇcastˇeji horizontálnˇe) a t´ım zvˇetˇsuje puklinu [3].

Pˇri tvoˇren´ı puklin mohou nastat r˚uzné pˇr´ıpady - pukliny r˚uzné ˇs´ıˇrky, r˚uzného smˇeru a délky, ale jen urˇcité pukliny jsou pouˇzitelné pro konkrétn´ı pr˚umyslové aplikace. Tvoˇren´ı puklin je moˇzné ˇcásteˇcnˇe ovlivnit, napˇr´ıklad do urˇcitého smˇeru ˇci typu pukliny. Ovlivnˇen´ı tvorby prob´ıhá zmˇenou tlaku, teploty nebo zmˇenou kapaliny, která je pouˇzitá pro tvorbu puklin [1].

Systémy HDR vyˇzaduj´ı, aby horniny byly v hloubce rozpraskány. Tvoˇr´ı se (mi- nimálnˇe) dva vrty, do kterých se vtláˇc´ı kapalina, která vytvoˇr´ı pukliny v okol´ı konce vrtu. Tyto dva vrty jsou k sobˇe dostateˇcnˇe bl´ızko, tak, aby vytvoˇrené pukliny byly propojené a mohly mezi sebou komunikovat. Do jednoho vrtu se ˇcerpá voda, která se ohˇreje od okoln´ı horniny a pak se vrát´ı na povrch druhým vrtem.

Modelován´ı hydraulického ˇstˇepen´ı je ˇzádouc´ı a výhodné v pˇr´ıpadech, kdy chceme

(11)

napˇr´ıklad pˇredej´ıt neˇzádouc´ı seismicitˇe nebo ovlivnit vlastnosti propustné zóny.

1.3 Biot˚ uv model a jeho varianta pro prostˇ red´ı s pukli- nami

Maurice Anthony Biot (1905-1985) byl belgicko-americký aplikovaný fyzik, který se vedle pˇr´ıspˇevk˚u v oblasti termodynamiky, letectv´ı, zemˇetˇresen´ı, geofyziky a elektro- magnetismu, povaˇzuje za zakladatele teorie poroelasticity.

Prvn´ım krokem k modern´ı teorii poroelasticity bylo pˇredstaven´ı obecné troj- rozmˇerné teorie pruˇzné deformace v nasyceném tuhém mediu [4]. Tato teorie vˇsak nebyla jeˇstˇe úplnˇe dokonalá, nebot’ byla omezena na izotropn´ı pˇr´ıpad. V 50. letech minulého stolet´ı M. A. Biot svou teorii rozˇs´ıˇril na anizotropn´ı pˇr´ıpady a odvodil rovnice pro dynamickou odezvu porézn´ıch médi´ı [5], [6]. Pˇrestoˇze tato forma Biotovy teorie jiˇz z´ıskala mnoho úspˇech˚u a vyuˇzit´ı, stále ˇslo pouze o teorii s pˇredpokladem lineárn´ı elasticity. Aˇz o nˇekolik let pozdˇeji, v 70. letech 20. stolet´ı, byla tato teorie znovu upravena a zdokonalena pro obecnˇejˇs´ı pˇr´ıpady, kdy M. A. Biot zahrnul i moˇznost nelineárn´ı elasticity [7].

Základn´ı myˇslenky teorie poroelasticity spoˇc´ıvaj´ı v tom, ˇze tlak kapaliny pˇrisp´ıvá k celkovému namáhán´ı v porézn´ım médiu a samotný tlak kapaliny m˚uˇze napnout porézn´ı medium. V takovémto pórovitém médiu docház´ı k pohybu kapaliny kv˚uli rozd´ıl˚um tlaku kapaliny v pórech, které souvis´ı s mechanickým zat´ıˇzen´ım porézn´ıho média. V d˚usledku toho m˚uˇze doj´ıt k posunu, otv´ırán´ı nebo zav´ırán´ı puklin.

Nelineárn´ı chován´ı je obvykle asociováno se vznikem a ˇs´ıˇren´ım puklin [1]. Pro náˇs pˇr´ıpad uvaˇzujeme médium, které má pevný skelet (dále oznaˇcované jako matrice). Pˇredpokládáme, ˇze kapalina je nestlaˇcitelná a ˇze porézn´ı medium stlaˇcitelné je. Biotova teorie pracuje s pˇredpokladem rozkladu vnitˇrn´ıho napˇet´ı v plnˇe saturo- vaném médiu. Toto napˇet´ı lze rozdˇelit na dvˇe ˇcásti: rovnováha sil v horninˇe a zákon zachován´ı hmotnosti v tekutinˇe.

Rovnov´aha sil v horninˇe popisuje zmˇenu tvaru media (deformace matrice) a je pops´ana touto rovnic´ı:

−∇ · σ + α∇p = f , (1.1)

kde σ pˇredstavuje tenzor napˇet´ı, α je Biot˚uv koeficient, p popisuje tlak v p´orech a f pˇredstavuje objemovou s´ılu.

(12)

Materiálový zákon zahrnutý v Biotovˇe teorii se nazývá Hook˚uv zákon:

σ = 2µ(u) + λ(∇ · u)I. (1.2)

Zde µ, λ jsou Lamého parametry, (u) je tenzor malých deformac´ı, u je vektor posunut´ı a I je jednotková matice. Zákon zachován´ı hmotnosti v tekutinˇe je popsán rovnic´ı:

∂_t(Sp + α∇ · u) − ∇ · q = g, (1.3) kde prvn´ı ˇclen ve tvaru ˇcasové derivace je pˇr´ır˚ustek objemu kapaliny za jednotku ˇcasu.Vektor q se nazývá Darcyho tok nebo Darcyho rychlost a popisuje výtok na jednotku plochy. S znaˇc´ı storativitu (nˇekdy znaˇceno také jako 1/M ) a g pˇredstavuje zdroj objemové s´ıly.

Darcyho zákon pro tok tekutin porézn´ım médiem pˇredstavuje rovnice:

q = −κ∇p, (1.4)

kde κ je hydraulická vodivost media, µ popisuje viskozitu média, ∇p popisuje celkový pokles tlaku.

Rovnice (1.1) - (1.4) jsou známé jako tzv. Biot˚uv systém pro deformovatelné porézn´ı médium. Pro prostˇred´ı obsahuj´ıc´ı pukliny byl pomoc´ı tzv. dimenzionáln´ı redukce odvozen Biot˚uv systém, který bude popsán v kap. 3 a který je hlavn´ım pˇredmˇetem studia této práce.

(13)

2 Pˇ rehled numerick´ ych metod pro klasick´ y Biot˚ uv model

Jako témˇer vˇsechny parciáln´ı diferenciáln´ı rovnice (PDR), i Biot˚uv systém nemá (aˇz na nˇekolik specifických pˇr´ıpad˚u) analytické ˇreˇsen´ı, proto se vˇetˇsinou mus´ı vyˇreˇsit numericky.

Numerické metody, jako napˇr. metoda koneˇcných prvk˚u, aproximuj´ı PDR soustavu algebraických rovnic. Pro numerické ˇreˇsen´ı Biotových rovnic jsou v dneˇsn´ı dobˇe pouˇz´ıvány pˇredevˇs´ım dva pˇr´ıstupy: plnˇe implicitn´ı (monolitický) a sekvenˇcnˇe implicitn´ı (iteraˇcn´ı). Plnˇe implicitn´ı pˇr´ıstup vyˇzaduje ˇreˇsen´ı celé soustavy rovnic souˇcasnˇe, coˇz pˇrináˇs´ı výhodu bezpodm´ıneˇcné stability, ale vyˇzaduje robustn´ı metody a pˇredpomiˇnovaˇce.

Sekvenˇcnˇe-implicitn´ı pˇr´ıstup nab´ız´ı vˇetˇs´ı flexibilitu v návrhu kódu a pouˇzit´ı spe- cializovaných metod pro ˇreˇsen´ı d´ılˇc´ıch problém˚u, neˇz plnˇe implicitn´ı pˇr´ıstup [8].

2.1 Monolitick´ e (implicitn´ı) ˇ reˇ sen´ı - pˇ r´ım´ e a iteraˇ cn´ı metody

Numerické ˇreˇsen´ı rovnic (1.1) - (1.4) pˇredstavuje výzvu, protoˇze rovnice pruˇznosti a proudˇen´ı jsou plnˇe svázány a vzniklá matice má sloˇzitou strukturu. Pro monolitické ˇreˇsené se pouˇz´ıvaj´ı dva základn´ı typy metod - pˇr´ımé a iteraˇcn´ı.

Uvaˇzujeme soustavu

Ax = b (2.1)

s matic´ı A ∈ R^n×n a vektorem pravých stran b ∈ Rⁿ. Nejznámˇejˇs´ı pˇr´ımou metodou je ˇr´ıdká LU dekompozice (Gaussova eliminace) pro ˇreˇsen´ı lineárn´ı soustavy rovnic.

Jej´ı maticov´a reprezentace je

A = LU, (2.2)

kde L a U jsou trojúheln´ıkové (tedy snadno invertovatelné) matice. Tato metoda

(14)

je robustn´ı, velmi jednoduchá a pouˇzitelná obecnˇe pro regulárn´ı matice. Je to do- poruˇcená metoda pro systémy s aˇz nˇekolika tis´ıc´ı neznámými (v praxi i v´ıce, pokud je matice ˇr´ıdká) a proto m˚uˇze být dobrou volbou pro mnoho 2D problém˚u a menˇs´ı 3D problémy. Pro vˇetˇs´ı problémy se ˇr´ıdký LU rozklad zpomaluje a rychle se vyˇcerpá pamˇet’, protoˇze ˇr´ıdká matice se rychle zapln´ı [9]. Cena klasické Gaussovy eliminace a Gaussovy eliminace s pivotac´ı mˇeˇrená poˇctem aritmetických operac´ı pro regulárn´ı matici je stejná. Pˇr´ımý chod Gaussovy eliminace vyˇzaduje pro matici ˇrádu n × n aˇz O(n³), viz tabulku 2.1 [10]. Pro velké problémy jsou dobrou alternativou iteraˇcn´ı metody, které jsou flexibilnˇejˇs´ı a vyˇzaduj´ı o dost ménˇe pamˇeti.

Tabulka 2.1: Poˇcty jednotliv´ych operac´ı pˇr´ım´eho chodu Gaussovy eliminace.

Operace celkem

× ¹₃n³− ¹₂n²+ ¹₆n

− ¹₃n³− ¹₂n²+ ¹₆n : ¹₂n²− ¹₂n

Iteraˇcn´ı ˇreˇsiˇce pouˇz´ıvaj´ı v prvn´ım kroku poˇcáteˇcn´ı odhad ˇreˇsen´ı a kaˇzdá k − tá iterace vycház´ı z minulého spoˇcteného ˇreˇsen´ı. Tento proces se opakuje aˇz do doby, kdy metoda zkonverguje. Zastaven´ı iteraˇcn´ı metody je zajiˇstˇeno takzvanými zastavovac´ımi kritérii, pˇr´ıpadnˇe maximáln´ım moˇzným poˇctem iterac´ı.

V problémech numerické lineárn´ı algebry ˇcasto vznikaj´ı matice, které jsou velké a ˇr´ıdké, a mohou být ˇrádu nˇekolika milión˚u ˇci dokonce miliard. V jiných pˇr´ıpadech nemus´ı být matice pˇr´ıliˇs velká, ale m˚uˇze m´ıt komplikovanou strukturu, m˚uˇze napˇr.

obsahovat bloky vyjádˇrené maticovou funkc´ı, jeˇz je velmi drahé sestrojit. ˇCasto pak jedinou operac´ı, kterou jsme schopni s matic´ı levnˇe a efektivnˇe provést, je násoben´ı matice s vektorem.

Iteraˇcn´ı metody si vysvˇetl´ıme na ˇreˇsen´ı klasických systém˚u lineárn´ıch rovnic, ve tvaru Ax = b, kde A ∈ C^n×n je regulárn´ı matice a b ∈ Cⁿ. C´ılem iteraˇcn´ıch metod je nalezen´ı dobré aproximace x_k pˇresného ˇreˇsen´ı x v co nejmenˇs´ım poˇctu iterac´ı. U klasických iteraˇcn´ıch metod posuzujeme jejich konvergenci v asymptotickém smyslu (rychlost konvergence je kvantifikována limitou pro poˇcet iterac´ı rostouc´ı nade vˇsechny meze, pˇriˇcemˇz konvergence metody nebere v úvahu moˇzný poˇcáteˇcn´ı pˇrechodový jev). Naproti tomu metody Krylovových podprostor˚u dávaj´ı v pˇresné aritmetice pˇresné ˇreˇsen´ı soustavy v nejvýˇse n kroc´ıch. V praxi nás zaj´ımá dosaˇzen´ı dostateˇcnˇe pˇresné aproximace ˇreˇsen´ı v co nejmenˇs´ım poˇctu krok˚u i a pˇrechodový jev nemá u metod Krylovových podprostor˚u smysl [10, s. 168].

(15)

Mezi klasick´e iteraˇcn´ı metody patˇr´ı Jacobiho metoda a Gauss-Seidelova metoda.

Tyto metody jsou velmi jednoduché, avˇsak konverguj´ı pomˇernˇe pomalu a lze je pouˇz´ıt pouze na ostˇre diagonálnˇe dominantn´ı matice. Mezi ˇcasto pouˇz´ıvané iteraˇcn´ı metody patˇr´ı pˇredpom´ınˇené Krylovovy ˇreˇsiˇce. Pˇredpodm´ınˇen´ı soustavy je proces, kdy se snaˇz´ımˇe modifikovat p˚uvodn´ı soustavu rovnic tak, aby metoda pouˇzitá na transformovanou úlohu dosáhla mnohem rychlejˇs´ı konvergence, neˇz pˇri pouˇzit´ı na

´

ulohu p˚uvodn´ı [10]. Mezi obl´ıbené pˇredpomiˇnovaˇce patˇr´ı neúplná LU (ILU) faktorizace, nebo algebraický multigrid. Mezi Krylovovy ˇreˇsiˇce napˇr´ıklad patˇr´ı:

CG - metoda sdruˇzen´ych gradient˚u

BiCGStab - metoda bikonjugovaných gradient˚u GMRES - zobecnˇená minimáln´ı reziduáln´ı metoda

Metoda sdruˇzených gradient˚u (CG) ˇreˇs´ı soustavy lineárn´ıch algebraických rovnic Ax = b s reálnou symetrickou pozitivnˇe definitn´ı matic´ı pˇr´ıpadnˇe s hermitovskou pozitivnˇe definitn´ı matic´ı (zkratka CG pocház´ı z anglického conjugate gradients).

Jde o teoreticky velmi silný nástroj s pamˇet’ovˇe nenároˇcným algoritmem[10, s, 81].

Zobecnˇená metoda minimáln´ıch rezidu´ı (GMRES - z anglického generalized mi- nimal residual method) byla vytvoˇrená pro matice které nejsou pozitivnˇe definitn´ı.

Metoda je matematicky ekvivalentn´ı ke generalizovan´e metodˇe sdruˇzen´ych gradient˚u a k metodˇe ORTHODIR.

Metodu bikonjugovaných gradient˚u lze pouˇz´ıt na obecné nehermitovské matice.

Jej´ı princip je podobn´y jako u CG, ale nen´ı zaruˇcen´a konvergence.

2.2 Iteraˇ cn´ı (sekvenˇ cnˇ e implicitn´ı) ˇ reˇ sen´ı

Iteraˇcn´ı ˇreˇsen´ı je postup, kde je problém proudˇen´ı a mechaniky rozdˇelen. V prvn´ım kroku se vyˇreˇs´ı bud’ problém mechaniky nebo proudˇen´ı, následnˇe je vyˇreˇsen druhý problém jiˇz s pouˇzit´ım nejnovˇejˇs´ıch informac´ıch o ˇreˇsen´ı. V kaˇzdém ˇcasovém kroku se tento postup opakuje aˇz do doby, kdy ˇreˇsen´ı dosáhne poˇzadované pˇresnosti.

Nejznámˇejˇs´ımi iteraˇcn´ımi metodami pouˇz´ıvánými na soustavy rovnic popisuj´ıc´ı propojen´ı proudˇen´ı s mechanikou v porézn´ıch médiech jsou: undrained split, fixed- stress split, drained split a fixed stress strain split. Bylo prokázáno, ˇze posledn´ı dvˇe zm´ınˇené metody maj´ı problémy se stabilitou, narozd´ıl od prvn´ıch dvou metod, které jsou robustn´ı a tedy stabiln´ı.

Metoda undrained split spoˇc´ıvá v uloˇzen´ı konstantn´ı hmotnosti kapaliny a následného poˇc´ıtán´ı tlak˚u p^n+1/2 v poloviˇcn´ım ˇcasovém kroku a poté pⁿ⁺¹ a stanoven´ı p^n+1/2 = pⁿ− α_S¹∇ · (u^n+1/2− uⁿ).

(16)

Naopak v metodˇe fixed-stress splitting uloˇz´ıme konstantn´ı celkové objemové napˇet´ı [11]. Metodu fxed-stress split si pop´ıˇseme pomoc´ı pseudokódu:

1. D´ano: (pⁿ⁻¹_h , uⁿ⁻¹_h ), eps, max_it 2. (p^n,0_h , u^n,0_h ) := (pⁿ⁻¹_h , uⁿ⁻¹_h ) 3. P roi = 1, 2, ..., max_it:

• Pomoc´ı (p^n,i−1_h , u^n,i−1_h ) spoˇcti p^n,i_h

• Pomoc´ı (p^n,i_h u^n,i−1_h ) spoˇcti u^n,i_h

4. Pokud i < max_it, pak (pⁿ_h, uⁿ_h) := (p^n,i_h , u^n,i_h )

kde eps je daná pˇresnost a max_it je maximáln´ı poˇcet iterac´ı, pˇri kterém se zastav´ı ˇreˇsen´ı.

Mezi iteraˇcn´ımi schématy byla tato ˇsiroce pouˇz´ıvaná metoda prokázána jako bezpodm´ıneˇcnˇe stabiln´ı ve smyslu Von Neumannovy analýzy a globálnˇe konvergentn´ı, kdyˇz uvaˇzujeme m´ırnˇe stlaˇcitelný tok v homogenn´ım porézn´ım médiu [8].

(17)

3 Numerick´ e ˇ reˇ sen´ı Biotova modelu pro ob- last s puklinou

V t´eto kapitole se budeme vˇenovat ˇreˇsen´ı Biotova modelu pro oblast s puklinou.

Reˇs´ıme probl´ˇ em vstˇrikov´an´ı tekutiny do pukliny γ ve ˇctvercov´em modelu horniny Ω:

Obr´azek 3.1: Model dom´eny s puklinou.

kde γ := (¹₂, 1) × {¹₂} a Ω := (0, 1) × (0, 1). Následuj´ıc´ı model a jeho diskretizace byly odvozeny v diplomové práci Svˇetlany ˇSauˇsové [2].

Biotova poroelasticita je pops´ana syst´emem rovnic:

∂_t(Sp + α∇ · u) − κ4p = g

−∇ · (2µ[u] + λ(∇ · u)I) + α∇p = f )

na (0,T) × Ω

kde neznámé jsou tlak p a posunut´ı u v oblasti Ω. [u] := (∇u + ∇u^T/2) je symetrická ˇcást gradientu u. Storativita S, Biot˚uv koeficient α, hydraulická vodivost κ, Lamého parametry µ, λ, volumetrický zdroj tekutiny g a s´ıla f jsou dány jako konstantn´ı parametry. Puklina m˚uˇze být prázdná (µ = λ = 0), nebo vyplnˇená elastickým materiálem (µ, λ > 0).

Pro redukovaný Biot˚uv model pˇredpokládáme, ˇze tlouˇst’ka pukliny δ je vzhledem k velikosti domény zanedbatelná. Po integraci Biotova systému pˇres ˇs´ıˇrku pukliny

(18)

z´ısk´ame n´asleduj´ıc´ı rovnice viz [2]:

δ{∂t(Sfpf + αf∇τ · uf) − ∇τ· (κf∇τpf)}

+F⁺+ F⁻− ∂_tG = δg_f na (0,T) × γ,

(3.1)

δ−∇_τ · (2µ_f_τ[u_f] + λ_f(∇_τ · u_f)I) + αf∇_τp_f

+Q⁺ + Q⁻− ∇_τ· R = δff na (0,T) × γ. (3.2) Zde p_f, u_f jsou neznámé v puklinˇe, S_f, α_f, κ_f, g_f, µ_f, λ_f, f_f jsou fyzikáln´ımi konstantami v puklinˇe, ∇_τ a _τ znaˇc´ı teˇcný gradient v γ a jeho symetrickou ˇcást.

V dalˇs´ı ˇcásti budeme pouˇz´ıvat toto oznaˇcen´ı: symbol n^± je jednotkový normálový vektor k γ v horn´ım (+) / doln´ım (-) smˇeru a p^±, u^± znaˇc´ı tlak a posunut´ı ze spodn´ıho p⁺, u⁺ / doln´ıho p⁻, u⁻ smˇeru. Neznámé F^±, Q^± reprezentuj´ı toky a normálová napˇet´ı p˚usob´ıc´ı na obou stranách pukliny a G, R jsou podm´ınky vznikaj´ıc´ı z dimenzionáln´ı redukce:

F^± := F^±(p, pf) = 2

δκf(pf − p^±), (3.3)

G := G(u) = α_f(u⁻− u⁺) · n⁺, (3.4)

Q^± := Q^±(p, u, u_f) = 2

δ(µ_f(u_f − u^±) + (µ_f + λ_f)((u_f − u^±) · n^±)n^±) +µ_f∇_τ(u^±· n^±) + λ_f(∇_fu^±)n^±− α_fp^±n^±,

(3.5)

R := R(u) = −µ^f(u⁻− u⁺) ⊗ n⁺+ λf(u⁻− u⁺) · n⁺I. (3.6) Syst´em (3.1) - (3.2) je spojen s rovnicemi (4) na Ω pomoc´ı podm´ınek na rozhran´ı:

κ∇p^±· n^± = F^±, (2µ[u^±] + λ(∇ · u^±)I)n^±− αp^±n^± = Q^± na (0, T ) × γ.

(3.7) Na ˇspiˇcce pukliny x_tip := [¹₂,¹₂]^T uvaˇzujeme homogenn´ı Neumannovy podm´ınky jak pro proudˇen´ı tak pro mechaniku.

Oznaˇcme si n´asleduj´ıc´ı ˇc´asti hranice: Γ_L := {0}×(0, 1) jako levou stranu, Γ_{T B} :=

(19)

(0, 1)×{0, 1} jako horn´ı a doln´ı stranu, Γ_R := {1}×(0, 1) jako pravou stranu hranice Ω. D´ale necht’ γ_R := {x_R}, kde x_R := [1,¹₂]^T je prav´y konec pukliny.

Rovnice Biotova syst´emu na Ω a na γ jsou doplnˇeny o tyto okrajov´e podm´ınky:

• Tekutina je vstˇrikov´ana do pukliny pod dan´ym tlakem p_in:

p_f(t, x_R) = p_in(t) :=

( sin (^π_Tt), t < T /2,

1, t ≥ T /2 pro t ∈ (0, T ).

• Tlak na lev´e stranˇe je pevnˇe dan´y

p = 0 na (0, T ) × Γ_L. (3.8)

• Deformace na horn´ı a doln´ı stranˇe je rovna nule:

u = 0 na (0, T ) × Γ_{T B}. (3.9)

• Prav´a strana stejnˇe jako vnˇejˇs´ı okraj pukliny se m˚uˇze deformovat pouze ve vertik´aln´ım smˇeru:

u₁ = 0 na (0, T ) × Γ_R, u_{f 1}= 0 na (0, T ) × γ_R. (3.10)

• Ve zbytku hranice uvaˇzujeme homogenn´ı Neumannovy podm´ınky.

V poˇc´ateˇcn´ım ˇcase t = 0 je syst´em v klidu, tedy:

p(0, .) = p_f(0, .) = 0 u(0, .) = u_f(0, .) = 0. (3.11) Pomoc´ı slabé formulace najdeme ˇreˇsen´ı (p, p_f, u, u_f) v odpov´ıdaj´ıc´ıch prostorech funkc´ı uspokojuj´ıc´ı poˇcáteˇcn´ı a okrajové podm´ınky takové, ˇze:

d

dt(S_p + α∇u, q)_Ω+ κ(∇p, ∇q)_Ω− (F⁺, q⁺)_γ− (F⁻, q⁻)_γ = (g, q)_Ω, (3.12)

δ d

dt(S_fp_f + α_f∇_τ · u_f,q_f)_γ+ κ_f(∇_τp_f, ∇_τq_f)_γ

+ (F⁺+ F⁻, q_f)_γ− d

dt(G, q_f)_γ = δ(g_f, q_f)γ,

(3.13)

(20)

2µ([u], [v])_Ω+ λ(∇ · u, ∇ · v)_Ω− α(p, ∇ · v)_Ω

− (Q⁺, v⁺)_γ− (Q⁻, v⁻)_γ = (f, v)_Ω,

(3.14)

δn

2µ_f(_τ[u_f], _τ[v_f])_γ+ λ_f(∇_τ · u_f, ∇_τ · v_f)_γ− α_f(p_f, ∇_τ · v_f)o + (Q⁺+ Q⁻, v_f)_γ+ (R, ∇τv_f)_γ= δ(f_f, v_f)_γ

(3.15)

pro vˇsechny testovac´ı funkce (q, q_f, v, v_f) a vˇsechny ˇcasy v (0, T ). Výraz (f, g)_Ω zde znaˇc´ı integrál souˇcinu f g pˇres Ω. Detaily odvozen´ı slabé formulace a definice prostor˚u lze nalézt v [2].

Pro numerické ˇreˇsen´ı pouˇzijeme nespojité Galerkinovy metody pro diskretizaci rovnic na Ω a metodu koneˇcných prvk˚u pro rovnice na γ. Necht’T_h je triangulace s´ıtˇe Ω, tj mnoˇzina vzájemnˇe disjunktn´ıch trojúheln´ıkových prvk˚u tak, ˇze ∪E∈T_hE = Ω.

Mnoˇzina vˇsech stran element˚u v Th bude oznaˇcena jako Eh a mnoˇzina vnitˇrn´ıch stran jako Eh,int. Pˇredpokl´ad´ame, ˇze Th je kompatibiln´ı s γ, tj je zde podmnoˇzina vnitˇrn´ıch stran Eh,γ ⊂Eh,int tak ˇze ∪_{F ∈}_E_h_,γF = γ.

Pro prostorovou diskretizaci na Ω uvaˇzujeme funkce které jsou po ˇcástech lineárn´ı na elementech z Th, ale obecnˇe nespojité na Ω. V puklinˇe uvaˇzujeme funkce po ˇcástech lineárn´ı na okraj´ıch E_h,γ a spojité na γ. D˚uvodem uˇzit´ı nespojitých pol´ı na Ω je, ˇze softwary zaloˇzené na metodˇe koneˇcných prvk˚u jako napˇr´ıklad FEniCS ne- podporuj´ı vnitˇrn´ı nespojitosti podél puklin. Spojitost na rozhran´ı vnitˇrn´ıch element˚u je vynucena dodateˇcnými podm´ınky ve formulaci.

Pro doˇcasnou diskretizaci bereme diskrétn´ı ˇcasový krok 4t := T /N , kde N je pozitivn´ı celé ˇc´ıslo. V diskrétn´ıch ˇcasech n4t budou funkce oznaˇceny horn´ım indexem ()ⁿ. ˇCasové derivace budou aproximovány implicitn´ı Eulerovou metodou, tj ∂_tf|t=n4t≈ ^fⁿ^−f_4tⁿ⁻¹.

Diskr´etn´ı probl´em je: Najdˇete {(pⁿ_h, pⁿ_{f h}, uⁿ_h, uⁿ_{f h})}^N_n=1 tak, ˇze 1

4t(Spⁿ_h + α∇ · uⁿ_h, q_h)_Ω+ κ(∇pⁿ_h, ∇q_h)_Ω− (F_hⁿ⁺, q⁺_h)_γ− (F_hⁿ⁻, q_h⁻)_γ

− X

F ∈Eh,int

(κ{{∇pⁿ_h}}, [[q_h]]n⁺)_F − (κ[[pⁿ_h]], {{q_h}}n⁺)_F

+ X

F ∈E_h,int

(Γ κ

|F |[[pⁿ_h]], [[q_h]])_F = (g, q_h)_Ω+ 1

4t(Spⁿ⁻¹_h + α∇ · uⁿ⁻¹_h , q_h)_Ω,

(3.16)

(21)

δn 1

4t(S_fpⁿ_{f h}+α_f∇_τ · uⁿ_{f hτ}, q_{f h})_γ+ κ_f(∇_τpⁿ_{f h}, ∇_τq_{f h})_γo + (F_hⁿ⁺+ F_hⁿ⁻, q_f)_γ− 1

4t(Gⁿ_h, q_f)_γ

= γ

(g_f, q_f)_γ+ 1

4t(S_fpⁿ⁻¹_{f h} + α_f∇_τuⁿ⁻¹_{f h} , q_{f h})_γ

− 1

4t(Gⁿ⁻¹_h , q_f)_γ,

(3.17)

2µ([uⁿ_h], [v_h])_Ω+ λ(∇ · uⁿ_h, ∇ · v_h)_Ω− α(pⁿ_h, ∇ · v_h)_Ω

−(Qⁿ⁺_h , v⁺_h)_γ− (Qⁿ⁻_h , v⁻_h)_γ

− X

F ∈Eh,int

(λ{{∇ · uⁿ_h}}, [[vh]]n⁺)F − (λ[[uⁿ_h]], {{∇ · vh}}n⁺)F

+2µ({{[uⁿ_h]}}n⁺, [[v_h]])_F − 2µ([[u_h]], {{[v_h]}}n⁺)_F

!

+ X

F ∈h,int

(Γµ + λ

|F | [[uⁿ_h]], [[v_h]])_F = (f, v_h)_Ω,

(3.18)

δ2µ_f_τ[uⁿ_{f h}], _τ[v_{f h}])_γ+ λ_f(∇_τ · uⁿ_{f h}, ∇_τ · v_{f h})_γ− α_f(pⁿ_{f h}, ∇_τv_{f h})_γ

+(Qⁿ⁺_h + Qⁿ⁻_h , v_{f h})_γ+ (Rⁿh, ∇_τv_{f h})_γ = δ(f_f, v_{f h})_γ, (3.19) pro vˇsechny testovac´ı funkce (q_h, q_{f h}, v_h, v_{f h}) a n = 1, 2, ...N . Zde F_h^n±, Gⁿ_h, Q^n±_h a Rⁿh znamená F^±(pⁿ_h, pⁿ_{f h}), G(uⁿ_h), Q^±(pⁿ_h, uⁿ_h, uⁿ_{f h}) a R(uⁿh). Dále {{f }} := (f⁺ + f⁻)/2, [[f ]] = f⁺− f⁻, a Γ > 0 je uˇzivatelem definovaný parametr (napˇr. Γ = 10³).

Zp˚usob ˇreˇsen´ı byl implementován do jiˇz existuj´ıc´ıho výpoˇcetn´ıho kódu pro re- dukovaný Biot˚uv systém. Abstraktn´ı zápis problému (3.16)-(3.19) vypadá v kódu takto:

a(s, t) = l(t) (3.20)

kde s je vektor ˇreˇsen´ı a t vektor testovac´ıch funkc´ı.

Numerické ˇreˇsen´ı Biotova modelu je implementováno pomoc´ı knihovn FEniCS - programován´ı koneˇcných prvk˚u neboli Finite Element (FE) Programming v jazyce Python. FEniCS je open-source výpoˇcetn´ı platforma zaloˇzená na ˇreˇsen´ı parciáln´ıch diferenciáln´ıch rovnic (Partial Diferential Equations - PDE), která umoˇzˇnuje uˇzivatel˚um pˇreloˇzit matematické modely do efektivn´ıho kódu koneˇcných prvk˚u (FE). Ve výpoˇcetn´ı platformˇe FEniCS je moˇzné programovat v C++ i v Pythonu, ale pro náˇs pˇr´ıpad

(22)

byl pouˇzit jazyk Python.

Vˇsechny v´ypoˇcty byly provedeny na poˇc´ıtaˇci s procesorem Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz 1.80GHz a operaˇcn´ı pamˇet´ı 8GB.

Pˇri výpoˇctech byly pouˇz´ıvány tyto hodnoty fyzikáln´ıch parametr˚u:

Tabulka 3.1: Pouˇzit´e parametry v jednotk´ach SI.

Parametr Hodnota

T 0,01

δ 0,1

k 0,1

k_f 1

S 0,01

Sf 0,1

f (0;0)

f_f (0;0)

µ 1

µf 0,5

λ 1

λf 0,5

α 1

α_f 1

3.1 Monolitick´ e ˇ reˇ sen´ı

Jak jiˇz bylo vysvˇetleno v kapitole 2.1, monolitické ˇreˇsen´ı soustavy znamená ˇreˇsen´ı soustavy v celku. Rovnice jsou zde zapsány ve formách a, l. Ve formˇe a je souˇcet levé strany rovnic (3.16) - (3.17), ve formˇe l je souˇcet pravé strany rovnic (3.16) a (3.17).

Ve výpoˇcetn´ı knihovnˇe FEniCS je výpoˇcet tedy zapsán v této formˇe:

Obrázek 3.2: Monolitické ˇreˇsen´ı zapsané v knihovnˇe FEniCS

kde solver parameters urˇcuje parametry naˇseho ˇreˇsiˇce: nejdˇr´ıve vyb´ıráme samo- statný ˇreˇsiˇc, poté pˇredpomiˇnovaˇc.

(23)

Pro náˇs problém jsme srovnávali ˇcasovou nároˇcnost LU faktorizace (bez

pˇredpodmiˇnovaˇce) a metody bikonjugovaných gradient˚u (pˇredpodmiˇnovaˇc ILU - ne- kompletn´ı LU faktorizace). Dalˇs´ı zkouˇsená metoda, GMRES, na tˇechto rovnic´ıch ne- zkonvergovala. Mˇeˇr´ıme ˇcas pouze samotného výpoˇctu, obˇe metody byly spuˇstˇeny na intervalu deseti ˇcasových krok˚u. Z tˇechto deseti krok˚u byl poté vypoˇc´ıtán pr˚umˇerný ˇcas výpoˇctu a smˇerodatná odchylka. Metody jsme srovnávali pro r˚uzné velikosti s´ıtˇe, kde n_x je poˇcet element˚u s´ıtˇe ve smˇeru x/y.

Tabulka 3.2: Srovnán´ı ˇcasového ˇreˇsen´ı ˇreˇsiˇc˚u s pˇredpodmiˇnovaˇci v závislosti na velikosti s´ıtˇe

Reˇˇ siˇc, pˇredpodmiˇnovaˇc Reˇˇ siˇc, pˇredpodmiˇnovaˇc

LU, ˇz´adn´y BICGSTAB, ILU

nx t pr˚umˇer [s] smˇerodatn´a odchylka t pr˚umˇer [s] smˇerodatn´a odchylka

8 3,44 0,033 3,67 0,183

16 12,99 0,118 13,16 0,348

32 50,11 3,041 54,07 0,541

64 207,57 2,638 228,08 1,357

V tabulce 3.1 je vidˇet, ˇze z tohoto srovnán´ı je rychlejˇs´ı LU metoda. BiCGstab je na tomto modelu pomalejˇs´ı, vzhledem k tomu, ˇze ILU je sice robustn´ı, ale ne ideáln´ı pˇredpodmiˇnovaˇc. Je potˇreba ale myslet na to, ˇze LU faktorizace nen´ı pouˇzitelná na velké soustavy rovnic.

Maximáln´ı velikost s´ıtˇe byla 64 × 64 × 2 element˚u, a ˇcasový nár˚ust ˇreˇsen´ı v závislosti na nx byl u LU faktorizace a metody bikonjugovaných gradient˚u kvadra- tický, m˚uˇzeme tedy oˇcekávat, ˇze LU faktorizace bude tedy i pro vˇetˇs´ı s´ıt’ rychlejˇs´ım ˇreˇsiˇcem, dokud nenastane problém zaplnˇen´ı pamˇeti.

(24)

Graf. 3.1: ˇCasové srovnán´ı dvou monolitických metod

3.2 Aplikace fixed stress splitting

Reˇsen´ı Biotova syst´ˇ emu iteraˇcn´ım zp˚usobem je pˇr´ıstup, který zahrnuje postupné implicitn´ı ˇreˇsen´ı proudˇen´ı a elasticity s vyuˇzit´ım nejnovˇejˇs´ıch informac´ı o ˇreˇsen´ı, kdy iterace prob´ıhá v kaˇzdém ˇcasovém kroku aˇz do konvergence.

Nam´ısto ˇreˇsen´ı systému (3.15) - (3.18) plnˇe propojeným zp˚usobem, je obl´ıbenou alternativou pouˇzit´ı iteraˇcn´ıch metod, které oddˇeluj´ı mechaniku od problému s proudˇen´ım a umoˇzˇnuj´ı efektivn´ı ˇreˇsen´ı a pˇredpom´ınˇen´ı jednotlivých podproblém˚u.

Zde omezujeme naˇse úvahy na ˇsiroce pouˇz´ıvanou metodu fixed stress splitting a pˇrizp˚usob´ıme myˇslenku Mikeliće a Wheelerové [11], která povaˇzuje zachován´ı umˇelé objemové konstanty napˇet´ı za umˇelou [8].

Iteraˇcn´ı schéma definuje sekvenci (p^n,i, u^n,i, p^n,i_f , u^n,i_f , kde i ≥ 0.) Po inicializaci u^n,0_h = uⁿ⁻¹_h , p^n,0_h = pⁿ⁻¹_h je kaˇzdý iterát definován ve dvou kroc´ıch. Prvnˇe je problém proudˇen´ı ˇreˇsen nezávisle. Poté je ˇreˇsen problém mechaniky (posunut´ı) pomoc´ı novˇe spoˇc´ıtaného tlaku v doménˇe a v puklinˇe. Do rovnic byl pˇridán relaxaˇcn´ı parametr β. Formy a a l lze rozdˇelit podle rovnic:

a = a_{f low}+ a_mech (3.21)

l = l_{f low}+ l_mech (3.22)

(25)

V´ypoˇcet tlaku a posunut´ı v abstraktn´ım z´apisu:

af low(sf low, tf low) = lf low(tf low), a_mech(s_mech, t_mech) = l_mech(t_mech),

(3.23)

kde s_{f low} := (p^n,i_h , p^n,i_{f h}) a s_mech := (u^n,i_h , u^n,i_{f h}) Pro fixn´ı n, i ∈ N vypad´a detailn´ı sch´ema:

a_{f low} = 1

4t((S + β)pⁿ_h, q_h)_Ω+ κ(∇pⁿ_h, ∇q_h)_Ω + δn 1

4t((S_f + β_f)pⁿ_{f h}, q_{f h})_γ+ κ_f(∇_τpⁿ_{f h}, ∇_τq_{f h})_γo

− (F_hⁿ⁺, q_h⁺)γ− (F_hⁿ⁻, q⁻_h)γ

− X

F ∈h,int

(κ{{∇pⁿ_h}}, [[q_h]]n⁺)_F − (κ[[pⁿ_h]], {{q_h}}n⁺)_F

!

+ X

F ∈h,int

(Γ κ

|F |[[pⁿ_h]], [[q_h]])_F,

(3.24)

l_{f low}= (g, q_h)_Ω+ 1

4t(Spⁿ⁻¹_h + α∇ · uⁿ⁻¹_h , q_h)_Ω + δ

(

(g_f, q_f)_γ+ 1

4t(S_fpⁿ⁻¹_{f h} + α_f∇_τ · uⁿ⁻¹_{f h} , q_{f h})_γ )

− 1

4t(Gⁿ⁻¹_h , q_{f h})_γ

− 1

4t(α∇ · u^n,i−1_h , q_h)

− δ

1

4tα_f∇_τ · (u^n,i−1_f , q_f)_γ

+ 1

4t(G^n,i−1_h , q_{f h})_γ + 1

4t(βp^n,i−1_h , q_h)_Ω + δ

4t(β_fp^n,i−1_{f h} , q_{f h})_γ,

(3.25)

(26)

qⁿ_mech^± _n = 2

δ(µ_f(uⁿ_{f h}− u^n±_h ) + (µ_f + λ_f)((uⁿ_{f h}− u^n±_h ) · n^±)n^±) + µ_f∇_τ · (u^±_h · n^±)n^±

+ λ_f∇_τ· (u^n±_h · n^±)n^±

(3.26)

amech = 2µ([uⁿ_h], [vh])Ω+ λ(∇ · uⁿ_h, ∇ · vh)Ω

δ

2µf(τ[uⁿ_{f h}], τ[vf h])γ+ λf(∇τ · uⁿ_{f h}, ∇τ · vf h)γ

+ (Rⁿh, ∇_τv_{f h})_γ

+ (qⁿ⁺_mechh, v⁺_h)_γ+ (qⁿ⁻_mechh, v⁻_h)_γ

− X

F ∈h,int

(λ{{∇ · uⁿ_h}}, [[v_h]]n⁺)_F − (λ[[uⁿ_h]], {{∇ · v_h}}n⁺)_F

+ 2µ({{[uⁿ_h]}}n⁺, [[v_h]])_F − 2µ([[uⁿ_h]], {{[v_h]}}n⁺)_F

+ X

F ∈h,int

(Γµ + λ

|F | [[uⁿ_h]], [[vⁿ_h]])_F

(3.27)

l_mech = (f , v_h)_Ω+ δ(f_f, v_{f h})_γ + (αp^n,i−1_h ∇ · v_h)_Ω + (δα_fp^n,i−1_{f h} ∇_τ· v_{f h})_γ + α_f(p^(n,i−1)+_h n⁺, v_{f h}− v⁺) + α_f(p^(n,i−1)−_h n⁻, v_{f h}− v⁻)

(3.28)

Pro rovnice uvaˇzujeme tˇri relaxaˇcn´ı parametry: fyzikálnˇe motivovaný parametr β_cl a parametry β_λ, β_opt, odvozené analýzou Mikélica a Wheelera [11], které jsou platné pro homogenn´ı Lamého parametry. Parametry jsou dány takto:

β_cl = α²

2µ

d + λ, β_λ = α²

2λ, β_opt = α²

2(^2µ_d + λ). (3.29)

β_clf = α_f²

2µ_f df + λf

, β_λf = α²_f

2λ_f, β_{f opt}= α²_f 2(^2µ_d^f

f + λf). (3.30) kde d, d_f je dimenze. Biot˚uv model je 2D model, tedy dimenze dom´eny je d = 2, a dimenze pukliny je d_f = 1.

(27)

Obr´azek 3.3: Iteraˇcn´ı ˇreˇsen´ı zapsan´e v knihovnˇe FEniCS.

Metodu fixed stress splitting m˚uˇzeme zapsat pomoc´ı pseudok´odu takto:

1. D´ano: (pⁿ⁻¹_h , pⁿ⁻¹_{f h} , uⁿ⁻¹_h , uⁿ⁻¹_{f h} ), eps, max_it 2. (p^n,0_h , p^n,0_{f h}, u^n,0_h , u^n,0_{f h}) := (pⁿ⁻¹_h , pⁿ⁻¹_{f h} , uⁿ⁻¹_h , uⁿ⁻¹_{f h} ) 3. P roi = 1, 2, ..., maxit:

• Pomoc´ı ((pⁿ⁻¹_h , pⁿ⁻¹_{f h} , uⁿ⁻¹_h , uⁿ⁻¹_{f h} ) spoˇcti (p^n,i_h , p^n,i_{f h})

• Pomoc´ı (p^n,i_h , p^n,i_{f h}, u^n,i−1_h , u^n,i−1_h ) spoˇcti (u^n,i_h , u^n,i_{f h})

• Pokud |(p^n,i_h , p^n,i_{f h}) − (p^n,i−1_h , p^n,i−1_{f h} )|/|(p^n,i−1_h , p^n,i−1_{f h} )| < eps a |(u^n,i_h , u^n,i_{f h}) − (u^n,i−1_h , u^n,i−1_{f h} )|/|(u^n,i−1_h , u^n,i−1_{f h} )| < eps, pak jdi na 4.

4. Pokud i < max_it, pak (pⁿ_h, pⁿ_{f h}uⁿ_h, uⁿ_{f h}) := (p^n,i_h , p^n,i_{f h}, u^n,i_h , u^n,i_{f h})

Po aplikaci iteraˇcn´ıho schématu jsme ovˇeˇrovali pˇresnost ˇreˇsen´ı. Vzhledem k tomu, ˇze rovnice nemaj´ı analytické ˇreˇsen´ı, chyba se poˇc´ıtá jako absolutn´ı hodnota z rozd´ılu od monolitického ˇreˇsen´ı. Pro tlak je výpoˇcet následuj´ıc´ı:

e_p = |p^ite− p^mono|

|p^mono| (3.31)

kde pîte, p^mono je vektor stupˇn˚u volnosti pro tlak v Ω spoˇctený iteraˇcnˇe, resp. mono- liticky, pˇriˇcemˇz monolitické ˇreˇsen´ı se povaˇzuje za pˇresné, a obdobnˇe pro posunut´ı.

(28)

Tabulka 3.3: Absolutn´ı hodnota relativn´ı chyby iteraˇcn´ıho schématu oproti monoli- tickému ˇreˇsen´ı v koneˇcném ˇcase v závislosti na zvolené pˇresnosti.

|mono − ite|

eps e_p e_pf e_u e_uf

1,00E-05 2,40E-06 8,30E-08 1,00E-06 8,30E-08 1,00E-06 5,60E-07 8,30E-08 3,30E-08 5,00E-09 1,00E-07 4,20E-08 5,50E-09 2,30E-09 2,90E-10 1,00E-08 2,00E-09 2,20E-10 9,80E-11 1,10E-11

Naˇs´ım c´ılem bylo m´ıt tak pˇresné ˇreˇsen´ı, ˇze chyba bude menˇs´ı neˇz 1⁻⁶ v tlaku i posunut´ı (v doménˇe i v puklinˇe). Z tabulky 3.2 je vidˇet, ˇze tento poˇzadavek splˇnuje eps = 10⁻⁶. Nejvˇetˇs´ı chyba se vˇzdy vyskytovala v doménˇe, pˇriˇcemˇz tlak p byl naˇs´ı urˇcuj´ıc´ı hodnotou, podle které jsme pˇresnost ˇreˇsen´ı urˇcovali. V puklinˇe byla vˇzdy odchylka menˇs´ı, neˇz v doménˇe. Z tabulky 3.2 je vidˇet, ˇze pˇresnost ˇreˇsen´ı roste, a tak m˚uˇzeme usoudit, ˇze iteraˇcn´ı schéma je konvergentn´ı. Výpoˇcty pˇresnosti byly provedeny v ParaView, pomoc´ı funkce Calculator. ParaView je open-source aplikace pro analýzu a vizualizaci.

Soubˇeˇznˇe se zjiˇst’ován´ım postaˇcuj´ıc´ı pˇresnosti jsme testovali, který parametr β je pro iteraˇcn´ı schéma ten nejvhodnˇejˇs´ı. Prvn´ım testem bylo srovnán´ı iterac´ı a ˇcasové nároˇcnosti vˇsech tˇr´ı parametr˚u v závislosti na velikosti s´ıtˇe.

Tabulka 3.4: Mˇeˇren´ı poˇctu iterac´ı a ˇcasové nároˇcnosti parametr˚u v závislosti na velikosti s´ıtˇe pˇri eps = 10⁻⁶.

β

λ

β

cl

β

opt

nx n ite n ite n ite

16 37 33 16

32 34 33 16

64 34 33 17

Prvn´ım poznatkem je, ˇze jakýkoliv relaxaˇcn´ı parametr β zajiˇst’uje stabiln´ı poˇcet iterac´ı, nezávisle na velikosti s´ıtˇe. Druhým poznatkem je fakt, ˇze pˇri pouˇzit´ı β_opt klesl poˇcet iterac´ı o polovinu vzhledem k β_cl. Vzhledem k výsledk˚um z tohoto testu jsme usoudili, ˇze pro analýzu konvergence pouˇzijeme relaxaˇcn´ı parametr β_opt.

(29)

3.3 Anal´ yza konvergence

V této kapitole se soustˇred´ıme na vliv r˚uzných parametr˚u výpoˇcetn´ıho modelu na poˇcet iterac´ı.

Naˇs´ım c´ılem je zjistit, jaký násobek relaxaˇcn´ıho parametru βoptje ten nejvhodnˇejˇs´ı, proto jsme provedli test závislosti poˇctu iterac´ı na hodnotách relaxaˇcn´ıho parametru β := c · β_opt pro násobek c ve zvoleném rozsahu: [¹₂;³₅;₁₀⁷;⁴₅; 1;³₂; 2;⁵₂; 3;⁷₂; 4;⁹₂; 5]. Zde mluv´ıme o násobku c obou parametr˚u β a β_f, vzhledem k tomu, ˇze parametry vˇzdy mˇen´ıme zároveˇn a jejich vliv na poˇcet iterac´ı je propojen.

Graf. 3.2: Závislost poˇctu iterac´ı na hodnotˇe beta pˇri nx = 16 a eps = 10⁻⁶. Z grafu 3.2 vid´ıme, ˇze pokud se relaxaˇcn´ı parametr rovná nule, poˇcet iterac´ı se limitnˇe bl´ıˇz´ı k nekoneˇcnu; takovýto kód obvykle nezkonverguje. Nejmenˇs´ı poˇcet iterac´ı β je v hodnotˇe β = 0, 25, coˇz je rovno 1 ∗ beta_opt, dále hodnoty lineárnˇe stoupaj´ı. Vzhledem ke zvoleným parametr˚um v grafu 3.3 splývaj´ı hodnoty β_cl, β_λ.

V dalˇs´ı analýze jsme zkoumali vliv parametr˚u Biotova modelu, které maj´ı pˇr´ımý vliv na parametr β. Relaxaˇcn´ı parametr β závis´ı na tˇechto tˇrech parametrech: α, µ a λ, a ekvivalentnˇe parametr βf závis´ı na αf, µf a λf, je-li α = αf = 0, pak rovnice elasticity a rovnice proudˇen´ı nejsou propojené. V´ıme, ˇze s rostouc´ım α, αf roste i poˇcet iterac´ı. V dalˇs´ıch testech jsme zjiˇst’ovali souvislost mezi tˇemito parametry výpoˇcetn´ıho modelu a relaxaˇcn´ım parametrem. Následuj´ı testy, kdy jsme postupnˇe mˇenili hodnoty tˇechto fyzikáln´ıch parametr˚u: (µ, µ_f), (λ, λ_f), (λ_f, µ_f), (λ, λ_f, µ, µ_f), (κ, κ_f).

Pˇri závislosti na zmˇenˇe druhého Lamého parametru (µ, µ_f) jsme postupnˇe provˇeˇrili dvojice parametr˚u (1; 0, 5), (10; 5), (0, 1; 0, 05) a (1; 1), (1; 0, 1), (1; 0, 01). Prvn´ı dvojice parametr˚u (1; 0, 5) je p˚uvodn´ı sada parametr˚u, tedy graf závislosti poˇctu iterac´ı na parametru β je stejný jako graf 3.2.

(30)

Graf. 3.3: Závislost poˇctu iterac´ı na parametru β pˇri druhém Lamého parametru (µ, µ_f) = (10; 5).

Z grafu 3.3 m˚uˇzeme vidˇet, ˇze pˇri zmˇenˇe velikosti druhého Lamého parametr˚u v doménˇe i v puklinˇe se zmˇen´ı i nejvhodnˇejˇs´ı hodnota relaxaˇcn´ıho parametru β_opt, z p˚uvodn´ıho c = 1 na c = ₁₀⁷ . Pokud se tedy zvˇetˇs´ı drhué Lamého parametry, klesne hodnota optimáln´ıho násobku relaxaˇcn´ıho parametru a poˇcet iterac´ı s velikost´ı hodnoty β prudce lineárnˇe roste. β_λ i β_cl pˇri této sadˇe druhých Lamého parametr˚u poskytuj´ı vˇetˇs´ı poˇcet iterac´ı.

Graf. 3.4: Závislost poˇctu iterac´ı na parametru β pˇri druhém Lamého parametru (µ, µ_f) = (0, 1; 0, 05).

Z grafu 3.4 vid´ıme, ˇze pˇri (µ, µ_f) = (0, 1; 0, 05) je nejvhodnˇejˇs´ım násobkem c = 3/2. Pˇri zmenˇsen´ı druhého Lamého parametru se také parametry β_λ a β_cl pˇribl´ıˇz´ı k optimáln´ı hodnotˇe relaxaˇcn´ıho parametru.

Ve druh´em testu byl konstantn´ı parametr µ a postupnˇe zmenˇsov´an paremetr µf

tak, aby pomˇer tˇechto parametr˚u byl postupnˇe 1,₁₀¹,₁₀₀¹ .