Symulacja przepływu płynu przez przestrzeń porową skał i estymacja

i modelowania opartego na fizycznych podstawach mechaniki płynów

Istnieje kilka metod wyznaczenia przepuszczalności skał. Jedną z nich jest bezpośredni pomiar na próbce skalnej z zastosowaniem przepuszczalnościomierza, którego możliwości pomiarowe mieszczą się w zakresie ułamków milidarcy lub nanodarcy. Kolejną możliwość stanowią empiryczne formuły do wyznaczania wartości przepuszczalności, opracowane na podstawie wyników uzyskanych za pomocą NMR. Jedną z najnowszych metod obliczania przepuszczalności skał jest modelowanie z zastosowaniem komputerowej mechaniki płynów (CFD, ang. Computational Fluid Dynamics). Metoda bazuje na odtworzeniu rzeczywistego układu geometrycznego analizowanej próbki i przeprowadzeniu symulacji przepływu płynu, rozwiązując równania zachowania masy oraz pędu (Arns et al., 2004; Fredrich et al., 2006; Bielecki et al., 2012).

Utworzenie trójwymiarowego modelu geometrycznego przestrzeni porowej, niezbędnego do wykonania symulacji, było możliwe dzięki wynikom komputerowej mikrotomografii rentgenowskiej. Analizę przepływu płynów przez przestrzeń porową przeprowadzono przy użyciu Metody Objętości Skończonych (MOS, FVM, ang. Finite Volume Method). Metoda polega na wygenerowaniu siatki numerycznej w wyniku podziału całego modelu geometrycznego na objętości skończone (objętości kontrolne, komórki kontrolne) i obliczaniu poszukiwanych wartości w każdej komórce siatki. Metoda objętości skończonych jest najbardziej popularną metodą dyskretyzacji równań zachowania masy, pędu (równanie Naviera-Stokesa) i energii w komputerowej mechanice płynów. Wynikami przeprowadzonej symulacji przepływu płynu były wartości prędkości lokalnej oraz ciśnienia w każdej komórce siatki, w zależności od różnicy ciśnienia na wlocie i wylocie analizowanej próbki (warunki brzegowe). Przepuszczalność obliczono na podstawie uzyskanych wyników symulacji (średnia wartość prędkości na wylocie próbki), własności analizowanego płynu i zmodyfikowanego prawa Darcy (Dvorkin et al., 2008; Kaczmarczyk et al., 2011; Krakowska & Madejski, 2013).

Trójwymiarowe obrazy mikrotomograficzne zostały użyte do utworzenia modelu geometrycznego przestrzeni porowej. Obiektem analizy była karbońska próbka piaskowcowa 888. Wybór próbki był podyktowany występowaniem ciągłych połączeń porów i gardzieli porowych w próbce oraz składem mineralnym, w szczególności najniższą zawartością minerałów ilastych. Analizowana próbka składa się w 93% z kwarcu, 2% z dolomitu i anhydrytu oraz około w 1% z minerałów ilastych. Porowatość całkowita, porowatość

efektywna z MP oraz przepuszczalność wynosiły odpowiednio 13,3%, 2,2% and 94,4 mD. Porowatość całkowita z analizy micro-CT wynosiła 12,6%. Figury 3.9.1 i 3.9.2 prezentują trójwymiarowy obraz przestrzeni porowej z pomiaru micro-CT. Ze środka całkowitego obrazu mikrotmograficznego przestrzeni porowej próbki 888 o wymiarach 950x950x400 wokseli został wycięty fragment o wymiarach 100x100x100 wokseli (fig. 3.9.1). W celu oszacowania przepuszczalności wykonano symulację przepływu wody i azotu przez przestrzeń porową próbki.

Pierwszym krokiem było obliczenie liczby Knudsena dla analizowanego przypadku, w celu przyporządkowania odpowiedniego reżimu przepływu. Liczba Knudsena stanowi podstawę oceny prawidłowego podejścia modelowego dla przepływu gazu w ośrodku porowatym. Istnieje kilka podejść stosowanych w komputerowej mechanice płynów, zależnych od liczby Knudsena, ważniejsze zostały zestawione w tabeli 3.9.1 (Civan, 2011).

Tabela 3.9.1. Rodzaje modeli stosowanych w symulacjach przepływu w zależności od liczby Knudsena

Model przepływu

płynu ^{Liczba Knudsena (Kn) Równania w modelu} ośrodek ciągły Kn≤0,001 ^{równanie Boltzmanna, Eulera,}

Naviera-Stokesa bez poślizgu

z poślizgiem 0,001<Kn<0,1 równanie Boltzmanna, Naviera-Stokesa z poślizgiem przejściowy 0,1<Kn<10 równanie Boltzmanna, Burnetta

swobodne molekuły 10≤Kn równanie Boltzmanna

Zakładając wartość średniej drogi swobodnej molekuł azotu w standardowych warunkach temperatury i ciśnienia równą λ = 10*10^-8 m (Couture & Zitoun, 2000) i długość charakterystyczną (średnica kanałów porowych w największym przewężeniu) L z micro-CT otrzymano liczbę Knudsena Kn=0,0008 dla analizowanej próbki. Oszacowana wartość liczby Knudsena pozwoliła na wybór podejścia modelowego z przedziału Kn≤0,001, co oznacza możliwość zastosowania równania Naviera-Stokesa bez poślizgu (brak poślizgu molekuł gazu na ścianach przestrzeni porowej). Komercyjne programy do komputerowej mechaniki płynów (m. in. zastosowany w pracy program Star-CCM+) bazują na rozwiązywaniu równań Naviera-Stokesa. W tym zakresie wartości liczby Knudsena możliwe jest stosowanie również innych równań, tj. równania Eulera i Boltzmanna. Równanie Boltzmanna przedstawia matematyczną interpretację zjawisk przepływowych na poziomie pojedynczych molekuł i stosowane jest głównie w zakresach niepozwalających na wykorzystanie założenia ciągłości ośrodka. W analizowanym w pracy przypadku, ze względu na odpowiednio duże wymiary geometryczne próbki, nie ma konieczności stosowania tego podejścia. Równanie Eulera

opisuje przepływ płynu nielepkiego, tzn. nieuwzględniający tarcia pomiędzy poszczególnymi warstwami płynu oraz pomiędzy płynem a ścianką. To podejście do modelowania przepływów w analizowanym przypadku nie może być stosowane z uwagi na bardzo duży wpływ lepkości płynu na ostateczne wyniki. W przypadku uzyskania wyższych wartości liczby Knudsena należałoby zastosować równanie Naviera-Stokesa z poślizgiem. W przypadku wystąpienia poślizgu molekuł gazu na ściankach porów, nie można traktować płynów jako ośrodków ciągłych. Jedną z konsekwencji tego zjawiska jest istnienie efektu Klinkenberga, który jest związany z poślizgiem molekuł gazu na ścianach porów przy niskim,

średnim ciśnieniu, skutkując zawyżeniem wartości przepuszczalności (Tiab, 2012). W przypadku symulacji przepływu cieczy poślizg na ściankach porów nie występuje, a więc nie ma wpływu na szacowane wartości przepuszczalności.

W celu sprawdzenia, czy symulacja może być przeprowadzona w zakresie przepływu laminarnego (równanie Naviera-Stokesa dla przepływu laminarnego w komputerowej mechanice płynów, reżim prawa Darcy w obliczeniach przepuszczalności) obliczono także liczbę Reynoldsa (Amao, 2007).

Zastosowane podejście do obliczenia przepuszczalności zostało szczegółowo opisane w pracach Narsillio et al. (2009), Zalewska et al. (2010) i Krakowska et al. (2013). Bazując na wynikach symulacji (rzeczywista prędkość przepływającego płynu na wylocie próbki) oraz na wynikach mikrotomografii (przekrój poprzeczny na wylocie próbki) wyznaczono objętościowy strumień przepływu na wylocie z analizowanej próbki dla różnych wartości spadku ciśnienia ∆p. Objętość przepływającego płynu q wzrastała liniowo podczas zwiększania gradientu ciśnienia ∆p, co jest charakterystyczne dla przepływów w zakresie stosowalności równania Darcy’ego. Wartość przepuszczalności została obliczona przy użyciu równania (3.9.1), wykorzystując parametry modelu geometrycznego przestrzeni porowej, tj. długość próbki L i przekrój poprzeczny na wylocie próbki A_p, porowatość próbki Kp ef, obliczoną wartość strumienia objętościowego q oraz wartość lepkości dynamicznej płynu

μ.

1 2 * * * A p L vdA Kpef K A p p

µ

∫

gdzie: K [mD] – przepuszczalność, Kp ef [ułamek] – porowatość efektywna, A_p [m²] – przekrój poprzeczny przez przestrzeń porową na wylocie próbki (Φ=A_p/A), p₁, p₂ [Pa] – ciśnienie na wlocie i wylocie próbki, L [m] – długość próbki, µ [Pa·s] – lepkość dynamiczna

płynu, v [m/s] – prędkość lokalna płynu na wylocie (wynik symulacji), A [m] – przekrój poprzeczny na wylocie próbki.

Przepływ płynu przez przestrzeń porową piaskowca karbońskiego symulowano z użyciem wody i azotu (fig. 3.9.2) w programie Star CCM+ firmy Cd-Adapco. Ustawiono następujące parametry do symulacji przepływu wody: gęstość wody ρ_w=997 kg/m³, lepkość dynamiczna wody µ_w=8,887*10^-4 Pa*s oraz azotu: ρ_a=f(p_a, T_a), gdzie p₂=p_a=100 kPa i T_a=273,15K, ρ_a=1,123 kg/m³ i µ_a=1,663*10^-5 Pa*s. Symulacja przepływu azotu była rozważana bez występowania poślizgu molekuł na ściance porów. Założono wartości ciśnienia na wlocie próbki: p1=60 lub 600 kPa oraz na wylocie: p2=100 kPa, przy długości próbki L=560*10^-6 m.

Symulacja przepływu wody została wykonana dla dwóch typów siatek (liczba komórek kontrolnych w MOS) dla spadku ciśnienia ∆p=60 Pa oraz dla jednego typu siatki, ale różnych spadków ciśnień ∆p=60 i 600 Pa (fig. 3.9.3). Symulacja przepływu azotu była przeprowadzona dla spadku ciśnienia równego ∆p=600 Pa.

Przepuszczalność karbońskiej próbki piaskowcowej została obliczona z użyciem wyników komputerowej mechaniki płynów (Andersson, 2012). Wykonano cztery symulacje dla różnych płynów i parametrów przepływu w celu potwierdzenia wiarygodności wyników i estymacji przepuszczalności. Tabela 3.9.2 przedstawia wyniki obliczeń przepuszczalności dla założonych parametrów w symulacji.

Tabela 3.9.2. Wynik obliczeń przepuszczalności dla założonych parametrów w modelu. Oznaczenia: εr – błąd względny

Lp. Płyn ^{Liczba komórek}

siatki ^{∆p = p1 – p2 K εr} [Pa] [mD] [%] 1 woda 790 495 60 73,3 22 2 222 777 60 74,2 21 3 222 777 600 73,8 21 4 azot 222 777 600 76,3 19

Uzyskane wartości przepuszczalności różnią się w zależności od użytego płynu. Przepuszczalność obliczona na podstawie wyników symulacji przepływu azotu miała najwyższą wartość i najbardziej zbliżoną do wartości uzyskanej z standardowych badań laboratoryjnych (przepuszczalnościomierz). Niewielkie różnice w wartościach można było zaobserwować dla przepuszczalności obliczonej na podstawie wyników symulacji przepływu wody. Najwyższą wartość przepuszczalności, równocześnie najbliższą wartości

przepuszczalności z badań laboratoryjnych, uzyskano dla spadku ciśnienia ∆p=60 Pa i obniżonej ilości objętości skończonych w modelu geometrycznym. Jednakże, wartości przepuszczalności uzyskane z wyników symulacji przepływu wody różnią się o 0,9 mD od uzyskanej z przepuszczalnościomierza, co można uznać za błąd obliczeń przepuszczalności z wyników symulacji.

Jednym z celów symulacji było sprawdzenie wpływu liczby objętości kontrolnych w modelu na wyniki symulacji. W przypadku skał niskoporowatych i niskoprzepuszczalnych dokładność obliczeń przepuszczalności odgrywa bardzo ważną rolę i musi być jak największa. Policzony został błąd względny, jako iloraz różnicy wartości przepuszczalności otrzymanej z badań laboratoryjnych i z wyników symulacji przez otrzymaną z badań laboratoryjnych (wyrażony w procentach), który oscylował wokół wartości 20%. Wyniki uzyskane dla różnych mediów (woda i azot), ciśnień na wlocie (spadkach ciśnienia) i liczb komórek kontrolnych w siatce modelu geometrycznego były zbliżone. Stanowiło to dowód na poprawne przeprowadzenie symulacji i założenie parametrów wejściowych, gdyż wartość przepuszczalności zależna jest od własności geometrycznych przestrzeni porowej.

Fig. 3.9.1. Trójwymiarowy obraz próbki piaskowcowej 888, okno programu Star CCM+ (zmodyfikowane)

Fig. 3.9.2. Trójwymiarowy model geometryczny przestrzeni porowej próbki piaskowca karbońskiego 888 z micro-CT,

okno programu Star CCM+

Fig. 3.9.3. Podział modelu geometrycznego na objętości skończone, okno programu Star CCM+ (zmodyfikowane)

Fig. 3.9.4. Rozkład prędkości wody dla gradientu ciśnienia ∆p=60 Pa, okno programu Star CCM+ (zmodyfikowane)

Fig. 3.9.5. Linie przepływu wody (kolorem oznaczono prędkość) dla gradientu ciśnienia ∆p=60 Pa, okno programu Star CCM+ (zmodyfikowane)

Fig. 3.9.6. Linie przepływu wody (kolorem oznaczone ciśnienie) dla gradientu ciśnienia ∆p=60 Pa, okno programu Star CCM+ (zmodyfikowane)

Figury 3.9.4 i 3.9.5 przedstawiają rozkład prędkości wody w przestrzeni porowej. Najwyższe prędkości wody zauważalne były w centralnej części próbki, w największych przewężeniach (minimalna średnica kanałów porowych). Analiza rozkładu linii przepływu wody informuje o głównych kierunkach przepływu strumienia wody w przestrzeni porowej. Wartości ciśnienia przedstawiono na figurze 3.9.6; osiąga ono najwyższe wartości na wlocie próbki i minimalne na wylocie. Szczegółowa analiza rozkładu pola prędkości i ciśnienia oraz linii przepływu płynu niesie informację na temat wpływu stopnia komplikacji układu geometrycznego przestrzeni porowej i krętości kanałów porowych na przepuszczalność skał.

Obliczenia przepuszczalności absolutnej karbońskiego piaskowca 888 zostały przeprowadzone na podstawie wyników symulacji komputerowej mechaniki płynów. Przepływ wody i azotu przez przestrzeń porową był przedmiotem analizy ilościowej (wielkości liczbowe parametrów) i jakościowej (rozkład parametrów w przestrzeni porowej).

Różnica między wartościami przepuszczalności obliczonymi w wyniku symulacji i otrzymanej z badań przepuszczalnościomierzem wynosi około 20%, definicję błędu względnego podano powyżej. Wartości przepuszczalności, obliczone w wyniku symulacji przepływu wody i azotu, a także różnych spadków ciśnień i liczby komórek kontrolnych w modelu geometrycznym są zbliżone. Porównanie uzyskanych wyników przepuszczalności wskazuje, że symulacje komputerowej mechaniki płynów mogą dostarczyć poprawnych wyników, jeżeli wartość liczby Knudsena przyjmuje wartości niższe lub równe niż 0,001.

Symulacje komputerowe przepływu płynów przez przestrzeń porową nie byłyby możliwe bez narzędzia pozwalającego uzyskać trójwymiarowy obraz przestrzeni porowej skały. Komputerowa mikrotomografia lub nanotomografia jest obecnie najlepszym narzędziem uzyskania obrazów przestrzeni porowej w sposób nieinwazyjny. W przypadku obliczania przepuszczalności przy użyciu metod komputerowej mechaniki płynów i tomografii rentgenowskiej należy brać pod uwagę czułość tych metod na poprawne wyodrębnienie przestrzeni porowej z całości obrazu mikrotomograficznego oraz wybór fragmentu reprezentatywnego dla całej próbki. Połączenie micro-CT i modelowania CFD stanowi przyszłościowe narzędzie do analizy przepływu mediów przez przestrzeń porową skał i obliczeń przepuszczalności.