computational fluid dynamics
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.
Geological structure of resolved partial soil
Scheme of resolve model with markers characteristic points
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.
Changing of temperature readied in characteristic points from model for reference case and with migration of ground water
Fields of temperature write in characteristic points from references case
Fields of temperature write in characteristic points from case with unsaturated soil
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.
Unstationary thermal analysis of the vertical heat ground exchanger within unsaturated soils
Wprowadzenie
Modelowanie migracji ciepła w gruntowych wymiennikach ciepła jest jednym z najważniejszych umiejętności podczas przewidywania czasu eksploatacji złoża, jego wydajności czy czasu trwania poszczególnych cykli pracy. Aby móc w sposób poprawny wyznaczyć poszczególne składowe strat niezbędny jest dokładny model obliczeniowy. Na jakość modelu obliczeniowego wpływ ma wiele czynników, w tym: czy przedstawia obraz poszukiwanych parametrów w sposób jedno, czy trójwymiarowy, opisuje otaczający grunt, jako ciało stałe, ośrodek porowaty… Wewnątrz gruntu zachodzą przemiany fazowe, istnieje przepływ wód gruntowych lub możliwość dyfuzyjnego przemieszczania się gazów, płynów wieloskładnikowych itd.
Zaproponowany model obliczeniowy oparty jest na metodzie objętości skończonych. Przewiduje trójwymiarowy opis otaczającego gruntu traktując go, jako ośrodek porowaty, wewnątrz którego w sposób dyfuzyjny przemieszczają się gazowe i płynne składniki płynów występujących w gruncie. Istnieje przepływ wód gruntowych, uwarunkowany naturalnymi procesami geologicznymi w rozpatrywanej parceli – Jabłonna koło Warszawy.
Dla dokładniejszego odwzorowania poszarpanego charakteru mierzonej temperatury przez czujniki w otworze roboczym, zaproponowano definicję temperatury całkowitej. Wzbogacona jest o dodatkowy wyraz opisujący dyfuzję zarówno ciepła jak i masy w elemencie podsypki oraz otaczającego gruntu. Dokładniejszy opis przedstawiono w dalszym punkcie. Wyniki wykazały dużo dokładniejsze odwzorowanie pomiarów od klasycznie stosowanych modeli obliczeniowych, co w przyszłości będzie miało znaczący wpływ na dobór urządzenia odbierającego ciepło z gruntu, wydajności danego złoża oraz szybkości degradowania w trakcie eksploatacji [1]
Analizując inne dostępne modele wyróżnić można Ngo i Lai [2]. Dodatkowym mechanizmem, oprócz klasycznego przewodzenia ciepła w ziemi, były ruchy dyfuzyjne spowodowane migracją wody. Autorzy podają graficznie zależności pomiędzy ilością odebranej energii, 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.
Autorzy [3] omówili 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 z uwagi na cząstkowy charakter prezentowanych wyników nie obejmował transportu oparów oraz przemian fazowych wody przesączającej grunt, co miało swoje odniesienie w krzywych weryfikujących model.
W pracy [3] 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. W modelu z uwagi na rozwiązywanie zagadnienia wymiany ciepła wokół przewodu emitującego silne pole elektro-magnetyczne, brakowało wpływu dedykowanego strumienia na odchylanie strugi migrującej wilgoci, czy przepływających wód podskórnych.
Autorzy [4,5,6] 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 przedstawionego problemu posłużono się modelem opisującym przepływ Darcy pomijając efekty zamarzania. Przepływ powietrza modelowany był, jako gaz idealny. Uzyskane wyniki uzyskiwały zbieżność nadającą się do zgrubnej analizy zasobów złoża.
Interesujące wyniki opublikowano w pracy [7]. Rozpatrując model ośrodka porowatego z przepływami wieloskładnikowymi płynu stwierdzono, że wzrost nasycenia gleby, przyjęty czas ładowania GHE, oraz szybkość przepływu czynnika zatłaczającego ciepło wpływa znacząco na zmiany temperatury w glebie. Ponadto zauważono, że regenerowanie akumulatora zbyt silnym strumieniem cieplnym w glebach zawilgoconych prowadzi do wysychania gruntu, a tym samym do pogorszenia właściwości przewodzenia ciepła.
Li i Cleall opublikowali wyniki, w których podkreślają wpływ zwartości wilgoci w glebie oraz wysokości poziomu wody gruntowej na wydajność pracy pompy ciepła, powiązanej GHE. Zaniedbanie zmian wilgoci może prowadzić do niedoszacowania pojemności cieplnej gruntu [8]. Podniesienie poziomu wód gruntowych miało korzystny wpływ na działanie powiązanej z wymiennikiem GSHP. Podkreślają również, że biorąc pod uwagę migrację energii cieplnej wraz z przepływem wód gruntowych ważne jest, aby uwzględniać wysokość poziomu lustra wody. Jego wahania mogą spowodować różnice wskazań temperatury od 3-4 %. Wartości te nie odbiegają od tych, jakie prezentują autorzy w niniejszej pracy.
W pracy [9,10] zaprezentowano model matematyczny opisujący dodatkowo przepływy konwekcyjne powietrza, pary wodnej oraz wód powierzchniowych. Wykazano, że około 83,3% energii magazynowanej jest przechowywane w medium do końca szóstego dnia, a 10% zmagazynowanej energii nie jest wykorzystywane w wyniku konwekcyjnego przepływu ciepła do końca 10 dnia.
Numeryczna analiza wymiany ciepła
Model fizyczny
Model fizyczny przedstawia klasyczny wymiennik typu U-rurka, umieszczony w otworze wywierconym w glebie. Posiada dwa kanały o przekroju kołowym, w których przepływa płyn, ścianki rozdzielające kanał od materiału zasypki oraz zasypkę od otaczającego gruntu. Szczegółowe wymiary pokazano na rys.1.
Materiał zasypki oraz płynu wewnątrz kanałów posiada stałe właściwości fizyczne, natomiast gruntu zmieniają się w zależności od głębokości, na której rozpatrywano wycinek wymiennika. Całkowita wysokość modelu to 80 m, a długość i szerokość to 1 i 0,7 m. Średnica otworu, to 0,18 m.
Wewnątrz gruntu oraz materiału zasypki istnieje możliwość przemieszczania się wody oraz powietrza. Nie istnieje kontakt pomiędzy płynem roboczym, a płynami przesiąkającymi. Płyn roboczy dostarcza ciepło do wymiennika, natomiast grunt i zasypka służą jako magazyn. Przesiąkające płyny odbierają/dostarczają dodatkowo ciepło do magazynu.
Structure and major dimensions of the vertical U-tube GHE model
Model gruntu posiada geometrycznie odwzorowaną porowatość. Polegało to, na zamodelowaniu szkieletu ośrodka porowatego. Stosunek powierzchni kanałów przepływowych do całkowitej powierzchni zajmowanej przez grunt odpowiada rzeczywistej porowatości np. dla żwiru i piasków wynosi 0,3.
W pracy porowatość gruntu odwzorowana została zarówno matematycznie, jak i geometrycznie.
Model matematyczny
Model matematyczny opisuje niestacjonarną wymianę ciepła dla rzeczywistych wymiarów opisujących wycinek gruntowego wymiennika ciepła. Pozwala uwzględnić wymianę ciepła pomiędzy płynem roboczym, a ciałem stałym, przewodzenie ciepła w ciele stałym oraz ośrodku porowatym.
Zawiera równania dyfuzyjne opisujące proces przewodzenia ciepła pomiędzy i-tymi składnikami przesiąkającego płynu, a szkieletem ośrodka porowatego.
Model matematyczny oparto na następujących założeniach:
- Gleba i materiały zasypki posiadają właściwości materiału izotropowego dla niewysyconego ośrodka porowatego.
- Właściwości fizyczne materiału nie zmieniają się w zależności od temperatury.
- Ścianka U-tube jest w dobrym kontakcie z materiałem zasypki, a także materiał zasypki jest w dobrym kontakcie z otaczającym gruntem. Opór cieplny styku poszczególnych powierzchni jest pominięty.
- Efekty związane z generowaniem ciepła w wyniku tarcia czynnika roboczego o powierzchnie ścianek rurociągu zostały pominięte.
- Efekty związane z generowaniem ciepła w wyniku oporów przepływu w ośrodku porowatym i-tego składnika płynu zostały, również pominięte.
- Nie istnieje zjawisko parowania lub zamarzania płynu w glebie. Powietrze w przestrzeni porów traktowane jest, jako gaz idealny.
- Zjawiska chemiczne w glebie wynikające z rozpuszczonych soli zostały pominięte.
Zgodnie z przyjętym modelem obliczeniowym geometrię modelu można podzielić na: płyn roboczy, U-rurkę, materiał zasypki oraz grunt. Przepływ płynu roboczego traktowany będzie jako turbulentny przepływ trójwymiarowy [11]. Przepływ ciepła przez ściankę U-rurki traktowano jako klasyczne przewodzenie. Materiał zasypki oraz gruntu opisano, jako niewysycony ośrodek porowaty. Wymiana ciepła w ośrodku porowatym odbywa się poprzez strumienie dyfuzyjne wewnątrz szkieletu oraz pomiędzy szkieletem, a i-tym składnikiem płynu opływającego ciało stałe.
Model matematyczny poruszającego się płynu
W przeprowadzonych symulacjach numerycznych, w których analizowano zagadnienia niestacjonarnego przewodzenia ciepła pomiędzy płynem, a ciałem stałym, posłużono się równaniami bilansu masy, pędu i energii [11]:
gdzie:
T, υ→, m, p, t↔ , I↔ , q→, b → , są odpowiednio: temperaturą, wektorem prędkości, lepkością molekularną, ciśnieniem, tensorem naprężeń lepkich, tensorem jednostkowym, strumieniem cieplnym Fouriera oraz siłą masową.
Wartość dysponowanej energii e = cpT + 1/2υ2 eqυ podzielono na dwie składowe: energie kinetyczną ruchu molekularnego, reprezentowaną przez pierwszy człon z prawej strony równania oraz energię kinetyczną ruchu makroskopowego zapisaną przez skalarny ekwiwalent wektora prędkości liczony wewnątrz elementarnej objętości:
Wymiana ciepła na ściance rurociągu Rura, w której przepływa płyn roboczy wykonana jest z polietylenu (PE), a posiadając dużą gęstość została potraktowana jako ciało stałe. Wymianę ciepła w ściance rury zamodelowano, jako niestacjonarne przewodzenie ciepła w ciele stałym wykorzystując formułę [11,12]:
Gdzie: ks to współczynnik przewodzenia ciepła w ciele stałym. Całkowita temperatura zastosowana w materiale gruntu oraz zasypki Definicje temperatury całkowitej zapisano, jako sumę molekularnego ruchu cząsteczek wyrażoną przez pierwszy człon z prawej strony oraz wynikający z kinetycznego ruchu makroskopowego zapisanego równaniem (3):
Gdzie: f1, f2, T, υ→, cp to współczynniki korekcyjne dobierane w trakcie symulacji oraz: temperatura, wektor prędkości, pojemność cieplna płynu, odpowiednio
Sposób rozwiązania symulacji numerycznej
Parametry fizyczne
Właściwości materiałowe dla gruntowego wymiennika ciepła typu U-rurka (GHE) wykorzystane w pracy pokazano w tab.1. Zastosowana siatka oraz warunki brzegowe
Siatkę objętości skończonych wykonano w programie Gambit 2.4. Wykonana siatka jest w pełni heksagonalną z ręcznym dopracowaniem miejsc wrażliwych, które mają duży wpływ na szybkość obliczeń oraz uzyskiwane niskie wartości liczb rezydualnych. W miejscach styku płynu roboczego ze ścianką ciała stałego wykonano dodatkowe zagęszczenie siatek, w celu odwzorowania profilu prędkości oraz temperatury. Model obliczeniowy składa się z 410 000 elementów. Szczegółowy widok siatki z zaznaczeniem obszarów z przypisanymi właściwościami materiałowymi oraz przyjętymi warunkami brzegowymi pokazano na rys. 2.
Obliczenia numeryczne przeprowadzono w programie Ansys Fluent 12.1. Rodzaje przyjętych warunków brzegowych oraz miejsca ich przyłożenia pokazano w tab.2.
W miejscu styku płynu roboczego ze ścianką rurociągu, zastosowano warunek brzegowy II rodzaju – Neumana. W miejscu styku nieznana była gęstość strumienia cieplnego.
Strumień cieplny Fouriera dla płynu rozumiany jest zgodnie z (4), jako suma części składowych pochodzących od dyfuzyjnego ruchu ciepła q→ oraz członu turbulentnego q→ t. W przypadku ciała stałego wielkość całkowitego strumienia ciepła przyjmuje wartości odpowiadające jedynie dyfuzyjnemu ruchowi ciepła [25,26].
Do opisania zjawiska turbulentnej wymiany ciepła pomiędzy ciałem stałym, a płynem zastosowano trójwymiarowy model turbulencji k – ∈. Równanie opisujące gęstość strumienia ciepła podzielono na część kondukcyjną i turbulentną [27-29].
Turbulentną liczbę Prandtla zdefiniowano poprzez zamienne wprowadzenie współczynnika lepkości:
Turbulentny współczynnik przewodzenia ciepła w płynie oraz turbulentny współczynnik lepkości zdefiniowano formułami:
Gdzie: k, ∈, jest energią kinetyczną turbulencji oraz intensywnością dyssypacji w 2 równaniowym modelu turbulencji k – ∈. Wielkość współczynnika korekcyjnego cm przyjęto 0,09.
![Tab.1 Physical parameters of the vertical U-tube GHE [18-24]](https://informacjainstal.com.pl/wp-content/uploads/2019/10/nistacjonarna-analiza-wody-tab.1.jpg "nistacjonarna-analiza-wody tab.1")
Physical parameters of the vertical U-tube GHE [18-24]
Grid mech and boundary of the vertical U-tube GHE warunki brzegowe dla pionowego wymiennika typu U-rurka
Weryfikacja eksperymentalna modelu
Analizowane otwory gruntowego wymiennika ciepła wykonane zostały w Centrum Badawczym Konwersja Energii i Źródła Odnawialne PAN w Jabłonnej koło Warszawy. Pełnią rolę dolnego źródła dla pracy pompy ciepła oraz jako miejsce zrzutu nadmiaru ciepła m.in. pochodzącego z kolektorów słonecznych, a gromadzonego w buforze ciepła niskotemperaturowego.
Wymiennik składa się z 52 otworów o głębokości 80 m rozmieszczonych w sposób równomierny na działce o wymiarach 25/26 m. Odległości pomiędzy poszczególnymi otworami wynoszą 2,5 m w obu kierunkach – rys.3.
Projektowana moc cieplna gruntowego wymiennika to 136 240 W, przy mocy pojedynczego otworu 2620 W. Nominalna moc cieplna wysokotemperaturowej pompy ciepła, to 100 000 W. Wymiennik zbudowany jest na 3 warstwach gleby: pierwszą i drugą stanowią żwir i piasek suchy i mokry kolejno. Trzecią warstwę wypełniają iły i muły z wstawkami węgla brunatnego.
Czujniki temperatury umieszczone zostały w otworze o oznaczeniu BHE22, na trzech głębokościach. Czujniki przymocowane zostały opaskami zaciskowymi bezpośrednio do ścianki rurociągu, po stronie zatłaczającej. Odczyt temperatury przekazywany był do komputera i zapisywany automatycznie w 10 minutowych odstępach czasu.
Eksperyment przeprowadzony został w terminie od 2018.05.23 od godz. 9:53 do 2018.05.24 9:45, czyli w okresie wczesnoletnim. Podzielony został na dwa etapy: process ładowania wynoszący 13,8 h oraz swobodnego stygnięcia trwający do zakończenia eksperymentu. Pomiary temperatury wykonywano za pomocą omówionych czujników pokazanych na rys. 3.
Temperatury początkowe dla gruntu przyjmowano, jako pierwsze rejestrowane przez czujnik i tak: na głębokości 25 m – 285,65 K, 50 m – 287,65 K oraz 75 m – 289,45 K. Strumień płynu roboczego wynosił 0,01 kg/s oraz posiadał temperaturę: 291,65 K.
Na rys. 4 pokazano wyniki dla przypadku wykorzystującego klasyczne równanie przewodzenia ciepła. W symulacjach tych, nie istniała migracja wilgoci oraz przepływy wód gruntowych. Model odwzorowuje poprawnie ładowanie, natomiast stygnięcie jedynie w górnych obszarach.
Rysunek kolejny – rys. 6 – przedstawia weryfikację danych eksperymentalnych przeprowadzoną przy pomocy modelu niewysyconego ośrodka porowatego. W modelu tym, uwzględniono możliwość migracji wilgoci oraz przepływy wód gruntowych.
Temperatury początkowe oraz wielkość strumienia płynu roboczego pozostały identyczne jak wyżej. Zadano strumień masowy dla przepływu wód gruntowych: 0,001 kg/s dla wody oraz 0,000 15 kg/s dla gazu. Temperatury obu składników były powiązane z temperaturą początkową obejmującą grunt w oddaleniu od otworu.
Wykonano symulacje przyjmując następujące wartości współczynników: macierz diagonalna dla gruntu: 3∙1010 dla wody oraz, 3∙1011 1/m2 dla gazu. Dla materiału zasypki: 3∙1011 i 3∙1010 1/m2 dla wody i gazu, odpowiednio. Macierz bezwładności dla gruntu: 350 1/m w przypadku wody oraz 350 dla gazu. Dla materiału zasypki, to: 350, 350 1/m dla wody i gazu, odpowiednio.
Model ośrodka porowatego z przepływami wieloskładnikowymi lepiej odwzorowuje swobodne stygnięcie akumulatora ciepła. Dzięki odpowiedniemu doborowi współczynników możliwe jest wyznaczenie dolnej i górnej granicy pola temperatur oraz bifurkacje rejestrowane w pomiarach.
Distribution of the vertical U-tube and the temperature measurement points for the GHE w wymienniku oraz czujników temperatury w otworze GHE
The experiment and simulation value of temperature for vertical U-tube oraz wyniki symulacji numerycznej dla wymiennika typu U-rurka
Experiment and simulation value of total temperature in porous model wartości temperatury całkowitej dla modelu porowatego
Wyniki symulacji oraz dyskusja
Rozkład temperatur w gruntowym wymienniku typu U-rurka
Rysunek 6 przedstawia pole temperatury odczytane w środkowym przekroju rozpatrywanego wycinka gleby na głębokości 75 m.
Pole temperatury stopniowo wzrasta w trakcie ładowania wymiennika, a następnie maleje podczas stygnięcia.
Po czasie 8,5 h ładowania – rys.6.c – obszar o największej temperaturze nie przekracza połowy rozpatrywanego obszaru gleby. Widoczne jest, także wyraźne przesunięcie izoterm w kierunku prawym. Wyraźniejsza migracja pola temperatur widoczna jest w trakcie procesu stygnięcia – obrazy na rys.6. d-f.
Za przemieszczanie się pola temperatur odpowiedzialna jest migracja wód gruntowych zadana w warunkach brzegowych. Przyjęte prędkości odpowiadają pomiarom, jakie wykonano w obiekcie w trakcie wykonywania odwiertów [20].
Wyliczenie wielkości temperatury odbywało się w sposób klasyczny, bez uwzględniania dodatkowych strumieni pochodzących od dyfuzyjnej wymiany masy pomiędzy i-tymi składnikami wody gruntowej, a materiałem gruntu.
The temperature field in the radial direction for vertical U-tube promieniowym dla pionowego wymiennika typu U-rurka
Pole wektorów prędkości w materiale zasypki
Rys. 7 przedstawia prędkości przepływu strumieni dyfuzyjnych odczytanych dla materiału zasypki. Czas rejestracji obejmował całkowity czas prowadzenia eksperymentu. Na poszczególnych wykresach widoczne są silne bifurkacje wartości spowodowane dużą niestacjonarnością tego typu przepływów. Bliżej to zagadnie, pokazano na rysunku 8.
Rys. 8 przedstawia pola prędkości wyliczone w płaszczyźnie promieniowej dla materiału zasypki. Wyniki podano dla głębokości 75 m.
Dla wszystkich zarejestrowanych kroków czasowych, widoczne są siatki połączonych kanałów kapilarnych, w których przemieszcza się wilgoć wraz ze składnikiem gazowym. Widoczna sieć połączeń nie jest stacjonarna, zmieniając swoją konfigurację z każdym następnym krokiem czasowym.
Wytworzona sieć połączeń pozostaje aktywna, pomimo wyłączenia przepływu płynu roboczego w widocznych otworach rurociągu. Drugi człon definiujący temperaturę całkowitą w równaniu (3) posiada duży wpływ na generowaną wartość temperatury całkowitej zwłaszcza w procesie stygnięcia wymiennika.
Rozkład temperatury całkowitej w materiale zasypki
Na rys. 9 pokazano pole całkowitych temperatur (3) odczytane w przekroju promieniowym dla materiału zasypki. Górny rząd odpowiada za ładowanie wymiennika, natomiast dolny za proces swobodnego stygnięcia.
Na rysunku widoczne są kanały o niższej temperaturze powodowane przez wnikanie wilgoci do wnętrza zasypki. Wnikająca wilgoć zanim osiągnie temperaturę materiału zasypki, powoduje lokalne schłodzenie ośrodka.
Dodanie strumieni dyfuzyjnych do całkowitego bilansu temperatury umożliwiło odtworzenie dużych wahań rejestrowanych przez czujniki insitu w trakcie prowadzenia eksperymentu.
Rysunek kolejny – rys. 10 – pokazuje liczby Reynoldsa odczytane na ściance oddzielającej materiał zasypki od otaczającego gruntu. Wyliczone niskie wartości wskazują na dominujący przepływ laminarny. Widoczne bifurkacje wskazują na miejsca, w których wytworzył się kanał kapilarny. Zmienność miejsca, w którym odczytywane są impulsy dla wybranych kroków czasowych potwierdza dużą niestacjonarność procesu.
Straty ciepła spowodowane działaniem strumieni dyfuzyjnych
Tab.3 przedstawia straty, jakie wnoszą strumienie dyfuzyjne. Całkowitą moc cieplną, jaka migruje w sposób dyfuzyjny porównano z mocą cieplna osiąganą przez otwór. Jej wielkość oszacowano na 1÷2% w czasie ładowania oraz 2÷4% w trakcie stygnięcia gruntu.
The Darcy velocity reading in backfill material odczytane dla materiału zasypki
The Darcy velocity field in the backfill material for deep 75 m
Calculation of losses of thermal energy resulting from diffusion streams
The total temperature field in backfill material from deep 75 m dla materiału zasypki na głębokości 75 m
The Reynolds numbers read on the wall separating the backfill from the surrounding soil
Wnioski
W pracy przeanalizowano model wymiennika ciepła posadowionego na gruntach niewysyconych. Zamodelowano geometrycznie i matematycznie porowatość gruntu, a także uwzględniono przepływ wód gruntowych w postaci ciekłej oraz gazowej. Zaproponowano definicję temperatury całkowitej uwzględniającą wpływ dyfuzyjnej wymiany masy na wielkość przechowywanego ciepła. Odwzorowano 24 h prace wymiennika odnosząc się do pomiarów wykonanych na 3 głębokościach w jednym otworze. Analizując uzyskane wyniki postawiono następujące wnioski:
1. Pole temperatur w obszarze zasypki zależy silnie od ilości wilgoci, która przesącza się przez materiał. Zjawisko to, najsilniej widoczne jest w procesie swobodnego stygnięcia gleby.
2. Ilość ciepła, która migruje wraz ze strumieniem dyfuzyjnym można przyjąć na, około 2-4 %.
3. Kanały kapilarne, którymi porusza się wilgoć nie mają charakteru stacjonarnego. Proces ten jest silnie zmienny w czasie wynikiem czego, obserwowane są bifurkacje temperatur rejestrowanych przez czujniki.
4. Występowanie wód gruntowych oraz ich przepływ powoduje przesunięcie pola temperatur w kierunku strugi prądu. Wielkość odchylenia zależna jest od szybkości oraz objętości przepływającego strumienia.
5. Zaproponowany model obliczeniowy dokładniej odwzorowuje wartości mierzone przez czujniki, umożliwiając w przyszłości dokładniejsze oszacowanie wydajności cieplnej oraz degradacji złoża w trakcie eksploatacji.
L I T E R AT U R A
[1] D. Sławiński, Numerical estimation of the thermal power available at a ground heat exchanger built on loamy silt, Acta Energetica 1/34 (2018) 12-17.
[2] C.C. Ngo, F.C. Lai, Heat transfer analysis of soil heating systems, International Journal of Heat and Mass Transfer 52 (2009) 6021–6027.
[3] J. Kačur, P. Mihala, M. Tóth, 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.
[4] P. Ocłoń, P. Cisek, M. Pilarczyk, D. Taler, 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.
[5] Z. Linlin, Z. Lei, Y. Liu, H. Songtao, Analyses on soil temperature responses to intermittent heat rejection from BHEs in soils with groundwater advection, Energy and Buildings 107 (2015) 355–365.
[6] Z. Wang, F. Wang, Z. Ma, Research 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.
[7] Yanling Guan, Xiaoli Zhao, Guanjun Wang, 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.
[8] A.K. Sani, R.M. Singh, Response of unsaturated soils to heating of geothermal energy pile , Renewable Energy, https://doi.org/10.1016/j. renene.2018.11.032
[9] Ch. Lia, P.J. Cleall, J. Mao, J.J Munoz-Criollo, Numerical simulation of ground source heat pump systems considering unsaturated soil properties and groundwater flow, Applied Thermal Engineering 139 (2018) 307–316.
[10] M.H. Jahangir, H. Sarrafha, A. Kasaeian, Numerical modelling of energy transfer in underground borehole heat exchanger within unsaturated soil, Applied Thermal Engineering 132 (2018), 697- 707.
[11] M.H Jahangir, M. Ghazvini, F. Pourfayaz, A numerical study into effects of intermittent pump operation on thermal storage in unsaturated porous media, Applied Thermal Engineering 138 (2018) 110–121.
[12] L.H. Dai, Y. Shang, X.L. Li, 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.
[13] Ansys Fluent, User’s Guide, Ansys Inc. April 2009.
[14] J. Darkwa, W. Su, D.H.C Chow, Heat dissipation effect on a borehole heat exchanger coupled with a head pump, Applied Thermal Engineering, 60 (2013), 243-241.
[15] M. Hruška, Ch. Clauser, R.W. de Doncker, Influence of dry ambient conditions on performance of underground medium-voltage DC cables, Applied Thermal Engineering 149 (2019) 1419–1426.
[16] Karsten Pruess, numerical simulation of multiphase tracer transport in fracture geothermal reservoirs, Geothermics 31 (2002) 475–499.
[17] Yan Gao, Rui Fan, HaiShan Li, Thermal performance improvement of a horizontal groundcaupled heat exchanger by rainwater harvest, Energy and Buildings 110 (2016) 302–313.
[18] K. Vafai, Handbook of porous media, Taylor & Francis, London 2005.
[19] T. Ahmed, Reservoir Engineering Handbook, Elsevier Inc. Oxford 2010.
[20] J. Kondracki, Regional geography of Poland, WNT, 2002.
[21] A. Alrtimi, M. Rouainia, S. Haigh, Thermal conductivity of a sandy soil, Applied Thermal Engineering 106 (2016) 551–560.
[22] Y. Li, Xu Han, Zh. Xiaosong, S. Geng, Ch. Li, Study of performance of borehole heat exchanger considering layered subsurface based on field investigations, Applied Thermal Engineering 126 (2017) 296–304.
[23] Y. Guo, G. Zhang, S. Liu, Investigation on the thermal response of full-scale PHC energy pile and ground temperature in multi-layer strata, Applied Thermal Engineering 143 (2018) 836–848.
[24] B. Bouhacina, R. Saim, H.F. Oztop, Numerical investigation of a novel tube design for the geothermal borehole heat exchanger, Applied Thermal Engineering 79 (2015) 153–162.
[25] Di Qi, Liang Pu, Zhenjun Ma, Lei Xia, Yanzhong Li, Effects of ground heat exchangers with different connection configurations on the heating performance of GSHP systems, Geothermics 80 (2019) 20–30.
[26] Giuseppe Brunetti, Hirotaka Saito, Takeshi Saito, A computationally efficient pseudo-3D model for the numerical analysis of borehole heat exchangers, Applied Energy 208 (2017) 1113-1127.
[27] A.B. Platts, D.A. Cameron, J. Ward, Improving the performance of Ground Coupled Heat Exchangers in unsaturated soils, Energy and Buildings 104 (2015), 323–335.
[28] Thomas Mimouni, Fabrice Dupray, Lyesse Laloui, Estimating the geothermal potential of heat-exchanger anchors on a cut-and-cover tunnel, Geothermics 51 (2014) 380–387.
[29] Rui Fan, Yiqiang Jiang, Yang Yao, A study on the performance of a geothermal heat exchanger under coupled heat conduction and groundwater advection, Energy 32 (2007) 2199–2209.
[30] Xinguo Li, Jun Zhao, Qian Zhou, Inner heat source model with heat and moisture transfer in soil around the underground heat exchanger, Applied Thermal Engineering 25 (2005), 1565–1577.