Numerical modelling velocity of thermal energy migration in the vertical heat ground exchanger within unsaturated soils

Wprowadzenie

Autorzy Ngo i Lai z Uniwersytetu w Północnej Dakocie opublikowali w pracy [12] wyniki symulacji numerycznej, której celem było oszacowanie wielkości ciepła, które przejmowane jest na powierzchni rozdzielającej powietrze z gruntem, w którym znajduje się wymiennik ciepła. Dodatkowym mechanizmem, oprócz klasycznego przewodzenia ciepła w ziemi, były ruchy dyfuzyjne spowodowane migracją wody, którą badany grunt był nasączony. Autorzy podają graficznie zależności pomiędzy ilością odebranego ciepła, a wielkością liczby Rayleigh’a opisującą mechanizm przenoszenia ciepła w opływającym płynie oraz liczby Darcy’ego opisującej przepuszczalność ośrodka porowatego czyli jego zdolność do przemieszczania się we wnętrzu danego płynu.

W pracy innych autorów: Koohi-Feyegh i Rosen [10] rozważano problem odwrotny, a mianowicie. Wpływ zainstalowanego gruntowego wymiennika ciepła na wzrost średniej temperatury gruntu w trakcie wieloletniej pracy wymiennika. Wpływ ten, będzie tym większy, im zwiększy się zagęszczenie wymienników ciepła na danym terenie oraz ich liczba. W przedstawionych symulacjach uwzględniono cykliczne sezonowanie wymiennika czyli naprzemienne dostarczanie i odbieranie ciepła z wymiennika.

W pracy autorów: Ocłoń P, Cisek P, Pilarczyk M i Taler D, rozważano rozkład pól temperatury uwzględniając różnorodność składu gleby, głębokość umiejscowienia źródła oraz miejscowe lub częściowe nasycenie gleby wodą. Jako źródło ciepła przyjęto kabel elektryczny przesyłający energię elektryczną w pobliżu jednej z elektrowni. Do obliczeń wykorzystano kod własny oraz wykonano weryfikacje na wynikach z kodu komercyjnego [13].

W kolejnej pracy [14] opisano działanie gruntowego wymiennika ciepła posadowionego na gruntach z występującą wilgocią. Autorzy podali trzy najważniejsze czynniki mające bezpośredni wpływ na wysokość uzyskiwanej temperatury oraz wydajność wymiennika. Są nimi: powierzchnia wymiany ciepła w rurze, tworzenie się miejsc o niskiej dyfuzyjności cieplnej – kawern powietrznych i nieciągłości gruntu oraz migracje ciepła spowodowane występowaniem cieków wodnych w miejscu posadowienia wymiennika.

Autorzy Linlin Z, Lei Z oraz Liu Y w pracy [11] podjęli się wykonania symulacji numerycznych gruntowego wymiennika ciepła posadowionego w pobliżu występowania cieków wodnych. Do symulacji wykorzystano własny kod obliczeniowy oraz dostępny na rynku kod komercyjny. Wykazano, że występowanie wód gruntowych oraz ich migracja mogą mieć wpływ na wielkość uzyskiwanej temperatury w wymienniku, zwłaszcza gdy pracuje on w systemie przerywanym – załadowanie, odczekanie, odbiór. Ustalono również, że nie tylko prędkość przepływu wód gruntowych, ale także właściwości gleby mają znaczący wpływ na sprawność pojedynczego otworu oraz grupy.

Innymi autorami podejmującymi temat przepływu ciepła w gruntach niewysyconych są Wang Zh, Ma F i Wang Xi [19]. Zwracają uwagę na fakt braku pogłębionych badań nad wpływem regeneracji wymiennika oraz wpływu migracji wód w glebach nienasyconych. Do rozwiązania postawionego problemu posłużono się własnym programem numerycznym, którego wyniki porównano z tożsamymi uzyskanymi z kodu komercyjnego. Wykazano, że uzyskane wyniki dobrze pokrywają się z danymi uzyskanymi przez innych autorów.

Autor niniejszego artykułu zgadzając się z wnioskami z wspominanej pracy uważa tym samym, że również istotne jest rozwijanie programów monitorująco-diagnozujących, zdolnych do przewidywania gotowości urządzenia do efektywnej pracy.

Autorzy: Guan Yanling, Zhao Xiaoli, i Wang Guan w pracy [6] wykorzystali własny program napisany w języku Fortran do analizy pól temperatury w gruntowym wymienniku ciepła, w którym występuje migracja wilgoci. W pracy wykonano analizę wpływu prędkości przepływu czynnika w otworze oraz szybkości przesiąkania wilgoci na wydajność cieplną pojedynczego otworu. Wykazano, że wartości oraz pola temperatur różnią się w zależności od rozważanego przypadku.

Autorzy: Kačur Jozef oraz Mihala Patric i Tóth Michal z uniwersytetu w Heidelbergu zaprezentowali wycinkowe efekty swojej pracy [9]. Omówiono w niej numeryczny model transportu masy i ciepła w nasyconych i nienasyconych ośrodkach porowatych. Migracja wód gruntowych napędzona została w wyniku działania sił grawitacji. Zastosowany model w uwagi na cząstkowy charakter prezentowanych wyników nie obejmował transportu oparów oraz przemian fazowych wody przesączającej grunt. Dla weryfikacji modelu wsparto się symulacjami numerycznymi. Autor publikacji, również w przyszłości zamierza rozwinąć swój model numeryczny oraz wykonać aplikację do programu monitorująco-diagnostycznego.

Ciekawe obserwacje opublikowało dwóch naukowców: Abubakar K. Sani oraz Rao M. Singh, z Uniwersytetu Surrey w Guildford w UK. Zajmując się badaniem stosu odwiertów tworzących gruntowy wymiennik ciepła oraz jego odpowiedzi na różne warianty pracy typu regeneracja/ odbiór zauważono, że regenerowanie gruntu zbyt silnym strumieniem cieplnym w glebach zawilgoconych prowadzi do wysychania gruntu, a tym samym do pogorszenia właściwości przewodzenia ciepła. Sytuacja ta, odbija się w późniejszym czasie bezpośrednio na zmniejszeniu wydajności cieplnej gruntu [15].

Jako uzupełnienie, przytoczyć można wyniki z pracy [4,17]. Zauważono w niej, że obszar wymiany ciepła dla pojedynczego otworu nie przekracza średnicy 1m. Po dłuższym czasie pola temperatury pomiędzy rurociągiem wlotowym a wylotowym rozmywają się, tworząc jeden wspólny obszar. Na polach odległych od U-rurki różnice pomiędzy temperaturą gruntu, a pochodzącą z otworu rozmywały się, a tym samym granica pomiędzy nimi, nie była jednoznaczna do wyznaczenia. Na bazie tych informacji obrano wielkość rozpatrywanego obszaru wokoło otworów.

Numeryczne symulacje CFD wykorzystano, także w pracach [2,7,19-23]. Symulacje te, pozwalały na oszacowanie odpowiedzi gruntowego wymiennika ciepła na krótko lub długoterminowe dostarczanie lub pobieranie ciepła. W pracy [20] dla procesu dostarczania ciepła (regeneracji) wykorzystywano energię pochodzącą z kolektora słonecznego. Największą sprawność takiego układu uzyskiwano dla pracy w reżimie 24 godzin, przy udziale ciepła dostarczanego w wielkości 42-52 %.

Potrzebę stworzenia programu projektowego dla gruntowych wymienników ciepła dostrzegają autorzy w pracy [3,16]. Program ten, umożliwia wyliczenie odpowiedzi gruntu na obciążenie termiczne w sposób długo, jak i krótkotrwały. Dokładność szacunków wykonanych przez zaproponowany model obliczeniowy sięgała od 99 do 92% dla dłuższego okresu eksploatacji.

Model obliczeniowy

Na rysunku poniżej (Rys.1) pokazano wycinek gruntu, który przynależy do rozpatrywanego otworu. Na rysunku zaznaczono poszczególne warstwy ziemi. Od góry znajduje się warstwa żwiru i piasków, we wnętrzu których nie stwierdzono występowania cieków wodnych.

Poniżej, znajduje się warstwa żwiru z piaskiem, w której występuje nasączenie ziemi wodą gruntową. Woda ta, wykazuje tendencje do przemieszczania się zgodnie z kinematyką opisaną przez równanie dyfuzji. W warstwie wodonośnej zaznaczono płaszczyznę – Ω – na której wykonano symulacje numeryczne. Szczegółowy opis zaprezentowano na rysunku następnym – Rys.2.

W centralnym miejscu rysunku widoczne są otwory: wlotowy i wylotowy z pojedynczego rurociągu. Całość otoczona jest wylewką bentonitową wypełniającą odwiert. Warstwa bentonitu zamodelowana została, jako ciało stałe, przez które nie wnika przesączająca się w gruncie woda. Na ściankach odwiertów przyłożony został warunek brzegowy II rodzaju opisujący konwekcyjną wymianę ciepła pomiędzy płynącym w otworze glikogenem, a ciałem stałym otworu. Wartości współczynnika przewodzenia zostały wyliczone oraz opublikowane w pracy [17].

Pozostały obszar obliczeniowy zamodelowany został, jako ośrodek porowaty o porowatości średniej, odpowiadającej gruntom typu piasek i żwir [1,18]. Wielkości poszczególnych współczynników zebrano w tabeli. 1. Ze względu na niewielkie liczby Reynoldsa występujące w ośrodku porowatym, przyjęto przepływ przezeń, jako laminarny.

Od strony strzałek przyjęto intuicyjnie wlot strumienia wody, natomiast po przeciwnej stronie jej ujście. Punktami charakterystycznymi zaznaczono miejsca, z których odczytano wyniki zaprezentowane w artykule. Pierwsze dwa punkty znajdują się za otworem akumulatora na kierunku przesączania się wody gruntowej. Trzeci, umiejscowiony jest pomiędzy otworem wlotowym i wylotowym, natomiast czwarty, znajduje się bezpośrednio przy ściance oddzielającej materiał bentonitu wypełniający odwiert, a żwirem wypełniającym objętość akumulatora.

Warunkiem początkowym był rozkład pola temperatury uzyskany po 8 godzinach ładowania akumulatora płynem o temperaturze 25oC przy warunkach brzegowych na ściance, identycznych jak w przypadku odbioru ciepła.

Fig.1
Geological structure of resolved partial soil
Fig.2
Scheme of resolve model with markers characteristic points
Tab.1
Solutions and boundary conditions

Wyniki symulacji numerycznych

Na rysunku poniżej – Rys.3 – zaprezentowano przebiegi w czasie temperatur odczytanych w charakterystycznych punktach pokazanych wcześniej na Rys. 2. Punkty te, umieszczono tak, aby odpowiadały drodze, jaką przebywa woda przesączająca się przez ośrodek porowaty, w tym przypadku – grunt. Linią ciągłą zaznaczono przebieg temperatur nazwany referencyjnym. Przebieg referencyjny uzyskany został dla procesu odbierania ciepła z akumulatora, przy zastosowaniu nasączenia gruntu wodą, lecz bez występowania ruchów dyfuzyjnych.

Przypadek badany, zawierał skład gruntu oraz nasączenie identyczne z rozważanym wcześniej, z tym, że dodatkowo dołączony został mechanizm dyfuzyjnego ruchu wody w kierunku zaznaczonym na Rys. 1 i 2. W obydwu przypadkach, symulacje wykonywano dla warunków początkowych zawierających pola temperatury uzyskane przy braku migracji płynu, po 8 h ładowania obszaru akumulatora temperaturą cieczy wynoszącą 26oC.

W punkcie pierwszym P1 jak i drugim P2, dla przypadku rozpatrywanego istnieje wyraźnie przejście fali wyższej temperatury oddalonej od siebie o pewien godzinowy odstęp czasu. Po tym czasie występuje efekt doliny z temperaturą niższą od tej, która uzyskiwana była dla przypadku nominalnego; „referencyjnego”. Wzrost temperatury tłumaczyć należy przejściem fali płynu, który na wstępie posiadał temperaturę wyższą z racji zalegania w okolicach bliższych otworowi, zatem o temperaturach bliższych tej zasilającej akumulator.

Punkt trzeci P3 umiejscowiony bezpośrednio pomiędzy otworem wlotowym, a wylotowym opisuje przebieg temperatury znacznie niżej od tej uzyskiwanej nominalnie, a wynikającej wyłącznie z mechanizmu dyfuzyjnego przewodzenia ciepła opisywanego strumieniem fourierowskim. Owe przechłodzenie bezpośrednio tłumaczone jest ochładzaniem obszaru otworu, wypełnionego bentonitem przez opływający wokoło ciek wodny.

Najwydatniej sytuacja ta, zarysowana jest dla przypadku czwartego oznaczonego jako P4. Tutaj, powolne obniżanie temperatury spowodowane przewodzeniem ciepła, zastąpione jest gwałtownym spadkiem. Punkt czwarty znajduje się bezpośrednio za otworem, zatem na linii prądu wody opływającej otwór. W tym przypadku spadek temperatury za otworem można określić mianem gwałtownego.

Należy przypuszczać również, że dalsze magazynowanie ciepła i jego dostępność ograniczą się wyłącznie do obszaru występowania bentonitu, silnie ochładzanego od strony zewnętrznej ścianki.

Na rysunkach kolejnych – Rys. 4 i 5 – pokazane zostały pola temperatury odczytane w wybranych krokach czasowych i uzyskane dla przypadku, w którym grunt był nasycony wodą, ale nie występowały ruchy zgromadzonej wody oraz dla przypadku, w którym ruch taki wystąpił.

Dla przypadku pierwszego, pobór ciepła odbywał się wyłącznie za pośrednictwem otworu wlotowego i wylotowego w wymienniku. W każdym z zaprezentowanych kroków czasowych izotermy zmniejszają się proporcjonalnie do środka, czyli kierunku poboru ciepła. Czas trwania rozładunku ciepła jest długi i wynosi 25 000 s (7 h).

W przypadku drugim, pobór odbywa się w dwojaki sposób, mianowicie: do środka w kierunku otworów wymiennika oraz w kierunku zgodnym z dyfuzyjnym ruchem płynu w ośrodku. Jak można zaobserwować wpływ dyfuzyjnego strumienia ciepła związanego z migracją wody gruntowej jest czynnikiem dominującym. Czas opróżnienia akumulatora dla warunków brzegowych przyjętych w symulacji wyniósł połowę czasu opróżniania w przypadku wcześniejszym uznanym, jako referencyjny – 10 000 s (3 h).

Wykresy kolejne – Rys. 6 – prezentują relatywną szybkość zmian wartości poszczególnych strumieni (entalpii, entropii i mocy cieplnej) odczytanych w charakterystycznych punktach modelu. Wielkość relatywna odnosi się do faktu, odjęcia od krzywej szybkości zmian strumienia wyliczonego dla przypadku z dyfuzyjnym ruchem wody, krzywej uzyskanej w modelu referencyjnym.

Na wykresie pierwszym, krzywe dla punktów P1 i P2 wykazują, skokowy przyrost wartości w odległościach odpowiadających szybkości przesuwania się fali ciepła. Efekt doliny tłumaczony jest poborem ciepła przez otwory wymiennika. Zasięg ich odziaływania obejmuje do punktu P1.

Największą szybkość odbioru odczytano dla punktu P4, który najintensywniej obmywany był przez strumień wody przesączający się przez grunt. Obserwację tę, potwierdza wykres następny przedstawiający szybkość zmian strumienia entropii. Największą szybkość przyrostu entropii odczytano dla P4. Po przejściu fali ciepła, szybkość ta stabilizuje się osiągając wartość ustaloną w chwili rozładowania akumulatora ciepła. Największą szybkość spadku entropii odnotowano dla P2, którego nie obejmował obszar początkowych pól temperatury, a w trakcie przejścia fali ciepła wraz z wodą gruntową, sukcesywnie był ogrzewany.

Interesujące wartości przedstawia wykres kolejny. Na trzecim wykresie – Rys.6 – zaprezentowano szybkość przyrostu i spadku mocy cieplnej w trakcie pobierania ciepła. Przypomnieć należy, również, że jest to wykres prędkości względnych, czyli takich, od których odjęta została wartość uzyskana dla rozładowania referencyjnego, a więc bez procesów dyfuzyjnego ruchu wód gruntowych. Największe szybkości przyrostu mocy cieplnej odczytano dla punktu P4 oraz P1. Pierwszy ze wspomnianych tak znaczny przyrost mocy zawdzięcza silnemu obmywaniu przez strumień wody przesączający się przez grunt. Odnotowana wielkość ma charakter piku, zatem objętościowe zagęszczenie mocy w tym miejscu będzie miało charakter chwilowy. Największą intensywność mocy źródła uzyska się w początkowym czasie jego eksploatacji.

Fig. 3
Changing of temperature readied in characteristic points from model for reference case and with migration of ground water
Fig. 4
Fields of temperature write in characteristic points from references case
Fig. 5
Fields of temperature write in characteristic points from case with unsaturated soil
Fig 6.
Relative change of velocity volume flux readied in characteristic points

Wnioski

W pracy przeanalizowano zachowanie się podłoża gruntowego z zamontowanym wymiennikiem ciepła. Rozważono dwa typy gruntów. Pierwszy był nasycony wodą, lecz bez ruchów dyfuzyjnych, drugi zaś posiadał narzucony wektor przemieszczania się wody. Wykonano symulacje numeryczne dla obydwu przypadków.

Zauważono, że przepływająca w gruncie woda pobiera pewną cześć ciepła. Proces ten może być korzystny, gdy z dalszych obszarów ciepło dostarczane jest do miejsca pracy otworu gruntowego lub niekorzystny, w przypadku jego odpływania.

Dla rozpatrywanego przypadku wyliczono przyrosty poszczególnych strumieni opisujących przepływ ciepła w złożu określając szybkość zmian poszczególnych strumieni. Zanotowano, że największe przyrosty i spadki w równaniach bilansu poszczególnych strumieni występują w pierwszym etapie eksploatacji złoża gruntowego. Zaobserwowane wyniki wykazały zgodność z obserwacjami opublikowanymi w czołowych czasopismach naukowych.

L I T E R AT U R A

[1] Ahmed T, Reservoir Engineering Handbook, Elsevier Inc. Oxford 2010.

[2] Cao Shi-Jie, Kong X-Ri, Investigation on thermal performance of steel heat exchanger for ground source heat pump systems using fullscale experiments and numerical simulation, Applied Thermal Engineering 115 (2017), 91-98.

[3] Cullin J.R, Spitler J.D, A computationally efficient hybrid time step methodology for simulation of ground heat exchangers, Geothermics 40 (2011), 144-156.

[4] Dai L.H, Shang Y, Li X.L, Analysis on the transient heat transfer process inside and outside the borehole for a vertical U-tube ground heat exchanger under short-term heat storage, Renewable Energy 87 (2016), 1121-1129.

[5] Darkawa J, Su W, Chow D.H.C, Heat dissipation effect on a borehole heat exchanger coupled with a head pump, Applied Thermal Engineering, 60 (2013), 243-241.

[6] Guan Y, Zhao Xi, Wang G, 3D dynamic numerical programing and calculation of vertical buried tube heat exchanger performance of ground – source heat pumps under coupled heat transfer inside and outside of tube, Energy and Buildings 139 (2017) 186–196.

[7] Hein P, Kolditz O, Görke U.J, A numerical study on the sustainability and efficiency of borehole heat exchanger coupled ground source heat pump systems, Applied Thermal Engineering 100 (2016), 421-422.

[8] Jahangir M.H, Sarrafha H, Kasaeian A, Numerical modelling of energy transfer in underground borehole heat exchanger within unsaturated soil, Applied Thermal Engineering 132 (2018), 697-707.

[9] Kačur J, Mihala P, Tóth M, Numerical modelling of heat exchange and unsaturated-saturated flow in porous media, Computers and Mathematics with Applications (2018 ) (in press) https://doi.org/10.1016/j.camwa. 2018.06.009.

[10] Koohi-Fayegh S, Rosen M, Long-term study of vertical ground heat exchangers with varying seasonal heat fluxes, Geothermics 75 (2018) 15–25.

[11] Linlin Z, Lei Zh, Liu Y, Songtao Hu, Analyses on soil temperature responses to intermittent heat rejection from BHEs in soils with groundwater advection, Energy and Buildings 107 (2015) 355–365.

[12] Ngo C.C, Lai F.C, Heat transfer analysis of soil heating systems, International Journal of Heat and Mass Transfer 52 (2009) 6021–6027.

[13] Ocłoń P, Cisek P, Pilarczyk M, Taler D, Numerical simulation of heat dissipation processes in underground power cable system situated in thermal backfill and buried in a multi-layered soil, Energy Conversion and Management 95 (2015) 352–370.

[14] Platts A.B, Cameron D.A, Ward J, Improving the performance of Ground Coupled Heat Exchangers in unsaturated soils, Energy and Buildings 104 (2015) 323–335.

[15] Sani A. K, Singh R. M, Response of unsaturated soils to heating of geothermal energy pile , Renewable Energy, https://doi. org/10.1016/j.renene.2018.11.032

[16] Sławiński D, Badania numeryczne i eksperymentalne dla potrzeb systemu nadzorującego prace wentylatorów w wymiennikach ciepła, Instal 2 (2017), 27-31.

[17] Sławiński D, termiczna adaptacja pionowego gruntowego wymiennika ciepła w wyniku oddziaływania cyklicznych obciążeń cieplnych, Instal 7-8 (2018)

[18] Vafai K, Handbook of porous media, Taylor & Francis, London 2005.

[19] Wang Zh, Wang F, Ma Zh, Reserch of heat and moisture transfer influence on the characteristics of the ground heat pump exchangers in unsaturated soil, Energy and Buildings 130 (2016) 140–149.

[20] Wołoszyn J, Gołaś A, Modelling of a borehole heat exchanger using a finite element with multiple degrees of freedom, Geothermics 47 (2013), 13-26.

[21] Yang H, Cui P, Fang F, Vertical-borehole ground- coupled heat pumps: A review of models and systems, Applied Energy 87 (2010), 16-27.

[22] Zhang Ch, Guo Z, Liu Y, A review on thermal response test of ground-coupled heat pump systems, Renewable and Sustainable Energy Reviews 40 (2014), 851-867.

[23] Zhang Ch, Wang Y, Liu Y, Computational methods for ground thermal response of multiple borehole heat exchangers: A review, Renewable Energy 127 (2018), 461-473.