Przedmiotem zainteresowań projektu GeoModel jest wykorzystanie niskotemperaturowych zasobów geotermalnych. Z tego powodu zakłada się, że przepływ płynu otworem geotermalnym jest przepływem jednofazowym, a dokładnie przepływem płynu w ciekłym stanie skupienia.
Górotwór traktujemy jako ośrodek ciągły, izotropowy, bez wewnętrznych źródeł ciepła, o stałych własnościach termofizycznych. Zakładamy, że w górotworze nie ma przepływu płynu – jedynym zachodzącym procesem jest wymiana ciepła, która zachodzi w warunkach nieustalonych. Do opisu zmian temperatury w czasie w górotworze sąsiadującym z otworem, wykorzystano schemat jawny metody różnic skończonych (MRS).
Opis geometrii modelowanego obszaru
Przyjęte założenia pozwalają rozważać proces we współrzędnych walcowych (Rysunek 1).
Rysunek 1. Strefa sąsiedztwa otworu we współrzędnych walcowych, schemat dyskretyzacji górotworu i lokacji otworu w przestrzeni
Bazując na poczynionych założeniach, zmiana temperatury w czasie opisana będzie równaniem:
Równanie 1
Wymiana ciepła w górotworze
Zgodnie z przyjętymi założeniami, opisu nieustalonej wymiany ciepła w górotworze można dokonać posługując się metodą różnic skończonych, dokonując uprzednio dyskretyzacji przestrzeni górotworu (Rysunek 1). Przy wykorzystaniu do opisu pochodnych względem współrzędnych przestrzennych ilorazów różnicowych symetrycznych o odpowiednim rzędzie, postać różnicową równania 1 można zapisać następująco (wyprowadzenie zależności i więcej informacji można znaleźć w literaturze (Wiśniewski i Wiśniewski, 1997)):
Równanie 2
Niewiadomą w powyższym równaniu jest ti,j,k+1 – temperatura w węźle i,j w kolejnym kroku czasowym k+1:
Równanie 3
Warunki początkowe
Ustalenie warunków początkowych jest jednym z warunków jednoznaczności rozwiązania równania Fouriera, opisanego we współrzędnych walcowych przez równanie 1. Warunki początkowe mogą oznaczać warunki, w których otwór zaczyna być eksploatowany po okresie długotrwałego przestoju. Jeżeli przestój otworu był długotrwały na tyle, by pole temperatury w górotworze powróciło do stanu naturalnego, warunki początkowe mogą odzwierciedlać krzywą profilowania temperaturowego otworu (rozkład temperatury z głębokością) – jeżeli jest ona znana. Jeżeli krzywa profilowania termicznego nie jest znana, to warunki początkowe można opisać zgodnie z regionalnym gradientem geotermicznym. W obu przypadkach temperatura w węzłach siatki obliczeniowej uzależniona jest jedynie od głębokości.
Innym przykładem, w których konieczne będzie ustalenie warunków początkowych to chęć odwzorowania zaburzonego pola temperatury w górotworze, spowodowanego np. wcześniejszą eksploatacją otworu. W tym przypadku ustalenie warunków początkowych wymaga przypisania temperatury wszystkim węzłom w modelowanym obszarze. Najwygodniej jest tego dokonać importując plik wsadowy o odpowiednim formacie zapisu danych. Będzie to ułatwione, jeżeli program, który określa rozkład temperatury w modelowanej przestrzeni będzie zapisywał rozkład temperatury jako jeden z plików wynikowych.
Warunki brzegowe
Warunki brzegowe są następnym z warunków jednoznaczności rozwiązania nieustalonego równania wymiany ciepła w modelowanym obszarze. Zaproponowano zastosowanie warunków brzegowych Dirichleta, tzw. warunków 1-szego rodzaju. Ich zastosowanie przynosi dużą stabilność modelu numerycznego, a implementacja w analizowanym procesie, sprowadza się do ustalenia temperatury w wybranych węzłach modelu (Rysunek 2). Pamiętać należy, że węzły z warunkami brzegowymi 1-szego rodzaju powinny być oddalone od otworu na tyle daleko, by zaburzenie termiczne spowodowane przepływem płynu nie dochodziło do węzłów z ustaloną temperaturą. W praktyce, dla analizowanego procesu oznacza to konieczność ustalenia warunku na pobocznicy walca rzędu kilkudziesięciu metrów (~50 m), w przypadku dna rzędu kilkunastu – do 20 m. Sugerowane rozpiętości modelu oszacowano bazując na (Pająk i Bujakowski, 2000).
Rysunek 2. Sugerowany schemat ustalenia warunków brzegowych
W przypadku strefy przypowierzchniowej dodano dodatkową warstwę elementów, na której ustalono warunek brzegowy 1-szego rodzaju. Zastępuje ona warstwę powietrza atmosferycznego oraz odzwierciedla efekty wymiany ciepła między powierzchniową warstwą gruntu, a powietrzem atmosferycznym. Gęstość i ciepło właściwe materiału zastępującego warstwę powietrza odpowiada parametrom typowym dla powietrza. Natomiast współczynnik przewodzenia ciepła jest wartością zastępczą, znacząco większą niż współczynnik przewodzenia ciepła nieruchomej warstwy powierza, która została ustalona tak, by uwzględnić warunki konwekcyjnej wymiany ciepła, zachodzące w warstwie powietrza, przy powierzchni Ziemi. Zakładając na podstawie danych klimatycznych, że średnia roczna prędkość wiatru w Polsce to ok. 3,24 m/s (Dygulska i Perlańska, 2015), to współczynnik konwekcyjnej wymiany ciepła szacować można jako (Recknagel, Sprenger, Honmann i Schramek, 1994):
Równanie 4
Jeżeli założymy, że wymiana ciepła powietrza atmosferycznego z gruntem powoduje zmianę temperatury powietrza w warstwie do Δz = 2030 cm i wprowadzimy zastępczy współczynnik przewodzenia ciepła λair’ o takiej wartości, że różnica temperatury między temperaturą powierzchni ziemi i powietrzem Δtsa pozostanie niezmieniona:
Równanie 5
Stąd przyjęta na rysunku 2 wartość zastępczego współczynnika przewodzenia ciepła powietrza λair’, uwzględniająca konwekcyjną wymianę ciepła z otoczeniem.
Wymiana ciepła między płynem a ścianą wewnętrzną otworu
Opisu rozkładu temperatury płynu przepływającego otworem dokonać można wykorzystując metodę bilansów strefowych. Zakładając, że w otworze zachowany jest bilans masy (m ̇ w całym wymienniku jest stałe), temperatura w analizowanej strefie wymiennika n (Rysunek 1) zmienia się liniowo od temperatury wejścia do strefy t(n)j-1, do temperatury wyjścia ze strefy t(n)j. Moc wymieniana między płynem płynącym w otworze, a górotworem powoduje zmianę temperatury płynu, zgodnie z równaniem:
Równanie 6
Równanie powyższe można rozwiązać iteracyjnie, dokonując określenia temperatury na wyjściu płynu ze strefy n. Daje to możliwość dokładnego określenia ciepła właściwego płynu w strefie, które znacząco wpływa na efekt jego nagrzania lub ochłodzenia.
Wymiana ciepła zachodząca między płynem a wewnętrzną ścianą rur okładzinowych jest konwekcyjnym wnikaniem ciepła i opisana jest równaniem:
Równanie 7
Określenie wartości temperatury twsf(n)k na powierzchni wewnętrznej rury otworu w strefie n, na wysokości j, a aktualnym czasie (kroku czasowym) k można dokonać korzystając z równania:
Równanie 8
Współczynnik konwekcyjnej wymiany ciepła może zostać określony przy wykorzystaniu liczby Nusselta Nu:
Równanie 9
W przypadku przepływu laminarnego w rurze (wartość liczby Reynoldsa < 2000), wartość współczynnika konwekcyjnej wymiany ciepła można obliczyć z następującej zależności (Wiśniewski i Wiśniewski, 1997):
Równanie 10
W przypadku przepływów turbulentnych (Re>2000), dla wody i innych cieczy o małej lepkości (1<Pr<20) skorzystać można z zależności empirycznej (Wiśniewski i Wiśniewski, 1997):
Równanie 11
W przypadku otworów eksploatujących płyn geotermalny należy się liczyć z tym, że nie będzie to czysta woda, a solanka. Z punktu widzenia modelowania parametrów płynu w otworze geotermalnym, mineralizacja wpływa na ciepło właściwe, gęstość i lepkość płynu. Parametry płynu, uwzględniające mineralizację mogą zostać określone na podstawie danych literaturowych. Równania takie zaproponowano na stronie projektu (IGSMiE PAN, 2024).
Opory przepływu
Rozkład ciśnienia we wnętrzu otworu można określić określając opory przepływu płynu geotermalnego przez każdą strefę n otworu ∆p(n) (Recknagel, Sprenger, Honmann i Schramek, 1994):
Równanie 12
Ciśnienie całkowite wody na końcu n-tej strefy otworu, w chwili k określane jest zależnością:
Równanie 13
Ciśnienie statyczne płynu na wyjściu ze strefy n, w chwili k wynosi ps(n)j,k – od którego zależą parametry płynu, takie jak: lepkość, gęstość i ciepło właściwe. Określa się odejmując od ciśnienia całkowitego panującego na początku strefy pc(n)j-1,k ciśnienie dynamiczne, spadek ciśnienia hydrostatycznego ∆ps(n)k i opory przepływu ∆p(n)k:
Równanie 14
Literatura:
Dygulska, A. i Perlańska, E. (2015). Mapa wietrzności Polski. Słupsk: Akademickie Centrum Czystej Energii (https://kierunkizamawiane.upsl.edu.pl/pliki/czystaenergia/raport2_II.pdf).
IGSMiE PAN. (2024, marzec 15). Geomodel. Pobrano z lokalizacji Parametry termofizyczne wody zasolonej: https://geomodel.pl/parametry-termofizyczne-wody-zasolonej/
Pająk, L. i Bujakowski, W. (2000). Efektywność wykorzystania głębokich odwiertów wiertniczych jako wymienników ciepła wykorzystujących energię geotermiczną. Kraków: Wydawnictwo Instytutu Gospodarki Surowcami Mineralnymi i Energią Polskiej Akademii Nauk.
Recknagel, H., Sprenger, E., Honmann, W. i Schramek, E. R. (1994). Poradnik Ogrzewanie i Klimatyzacja z uwzględnieniem chłodnictwa i zaopatrzenia w ciepłą wodę. Gdańsk: EWFE – Wydanie 1.
Wiśniewski, S. i Wiśniewski, T. S. (1997). Wymiana ciepła. Warszawa: Wydawnictwa Naukowe i Techniczne.
Wspólnie działamy na rzecz Europy zielonej, konkurencyjnej i sprzyjającej integracji społecznej
Projekt GeoModel jest dofinansowany ze środków Funduszu Współpracy Dwustronnej Mechanizmu Finansowego EOG na lata 2014–2021 i Norweskiego Mechanizmu Finansowego 2014-2021 w ramach Programu Środowisko, Energia i Zmiany klimatu.
Operatorzy Programu: Ministerstwo Klimatu i Środowiska