Modele matematyczne wielkoskalowych zasobników ciepła typu PTES /Mathematical models of large-scale PTES thermal energy storages – GWiTS – Gaz Woda i Technika Sanitarna – czasopismo Polskiego Zrzeszenia Inżynierów i Techników Sanitarnych

Afiliacja

¹ Wydział Inżynierii Środowiska, Politechnika Warszawska
* Kontakt / Correspondence: hubert.pugacz@pw.edu.pl

DOI: 10.65545/GWITS.2026.06.01, GWiTS 6/2026, czerwiec 2026

Streszczenie:
W artykule omówiono modele matematyczne cieplnych wielkoskalowych zasobników energii typu PTES (Pit Thermal Energy Storage) pracujących w warunkach nieustalonych. Omówiono modele bilansowe, półanalityczne oraz uproszczone modele numeryczne. Przeanalizowano sposób odwzorowania zjawisk przewodzenia, konwekcji, stratyfikacji termicznej oraz wymiany ciepła z gruntem. Wskazano zalety i ograniczenia poszczególnych podejść w kontekście dokładności odwzorowania procesów fizycznych, kosztów obliczeniowych oraz możliwości zastosowania w projektowaniu i optymalizacji systemów magazynowania ciepła.

Słowa kluczowe: modelowanie numeryczne, analiza stanu nieustalonego, zasobniki ciepła, PTES, magazynowanie ciepła

Abstract:
In the paper mathematical models used for the thermal analysis of large-scale Pit Thermal Energy Storage (PTES) systems operating under transient conditions have been described. Energy balance models, semi-analytical models, and simplified numerical models were discussed. The study analyzes the representation of heat conduction, convection, thermal stratification, and heat exchange with the surrounding ground. The advantages and limitations of the individual approaches were identified in terms of the accuracy of physical process representation, computational cost, and applicability in the design and optimization of thermal energy storage systems.

Keywords: numerical modeling, transient-state analysis, thermal energy storage tanks, PTES, heat storage

1. Wstęp

Systemy magazynowania ciepła umożliwiają rozłożone w czasie wykorzystywanie nadprodukowanego ciepła. Powszechność zasobników ciepła wzrasta wraz ze zwiększaniem się zainstalowanych mocy w intermitujących źródłach energii, przy których często występuje niespójność pomiędzy produkcją, a popytem. Magazynowanie energii jest również kluczowym elementem systemu w kontekście optymalizacji kosztów wytwarzania energii, ponieważ umożliwia stabilną pracę urządzeń wytwórczych z najwyższą sprawnością. Dodatkową korzyścią jest zmniejszenie wykorzystania paliw kopalnych. W zaawansowanych technologicznie systemach ciepłowniczych Unii Europejskiej często stosowane są tzw. hybrydowe systemy magazynowania ciepła. Opierają się one na wykorzystaniu kilku magazynów różnego typu w jednym systemie ciepłowniczym. W zależności od produkcji ciepła oraz wahań poboru optymalne mogą być magazyny o różnym czasie magazynowania ciepła.

2. Znaczenie modelowania numerycznego w analizie systemów magazynowania ciepła

Procesy zachodzące w wielkoskalowych zasobnikach ciepła mają charakter silnie nieustalony oraz obejmują jednoczesne występowanie przewodzenia, konwekcji i zjawisk wyporu cieplnego. Dodatkową trudność stanowi obecność stratyfikacji termicznej oraz wymiana ciepła pomiędzy zasobnikiem a otoczeniem. Złożoność tych zjawisk powoduje, że analiza pracy zasobników wymaga zastosowania zaawansowanych metod modelowania numerycznego. Wykorzystanie modelu numerycznego dokładnie opisującego magazyn ciepła pozwoli w procesie projektowania dobrać odpowiednią pojemność, geometrię czy parametry pracy, które umożliwią poprawę sprawności pracy systemu ciepłowniczego. Modele numeryczne są również podstawą do obliczeń optymalizacyjnych, umożliwiają efektywniejsze zarządzanie zasobnikiem. Najważniejszą przewagą nad badaniami eksperymentalnymi jest znacząco niższy koszt wykonania analizy oraz uniwersalność metody dla różnych rozwiązań konstrukcyjnych. W zakresie magazynowania ciepła w zasobnikach typu TES, model numeryczny przedstawiony w [1] umożliwia analizę rozkładu temperatury czynnika wewnątrz zasobnika w funkcji czasu magazynowania. Zastosowanie tego modelu pozwala określać straty ciepła oraz co za tym idzie umożliwia wyznaczanie optymalnych okresów magazynowania ciepła. Obecnie najszybciej rozwijającym się typem zasobników są wielkoskalowe zasobniki ciepła zagłębione w gruncie typu PTES (Pit Thermal Energy Storage). Jest to technologia magazynowania ciepła w zbiorniku wykopanym w gruncie. Najczęstszym medium magazynującym ciepło jest woda, choć zdarzają się również para wodna, żwir. Geometrię zasobników typu PTES uzależnia teren inwestycji oraz poziom zwierciadła wód gruntowych, które zaburzają pracę zasobnika. Celem modeli numerycznych opisujących zasobniki ciepła zagłębione w gruncie typu PTES jest odwzorowanie zjawisk termicznych i hydraulicznych zachodzących w zasobniku, jak i oddziaływania zasobnika na otaczający go grunt. W przypadku tego typu zasobników, szczególnego znaczenia nabierają zagadnienia związane z analizą procesów cieplno-przepływowych zachodzących w warunkach silnie dynamicznych i nieustalonych. Proces ładowania oraz rozładowania zasobnika obejmuje jednoczesne występowanie przewodzenia ciepła, konwekcji naturalnej i wymuszonej, zjawisk wyporu cieplnego, stratyfikacji termicznej oraz wymiany ciepła z otaczającym gruntem. Zjawiska te są wzajemnie sprzężone i zależne od czasu, co znacząco utrudnia ich analityczny opis. Ze względu na wielkoskalowy charakter instalacji PTES, prowadzenie badań eksperymentalnych jest procesem kosztownym oraz długotrwałym. W rezultacie modelowanie numeryczne staje się podstawowym narzędziem wykorzystywanym do projektowania, optymalizacji oraz oceny efektywności sezonowych magazynów ciepła. Rozwój metod numerycznych stosowanych w analizie zasobników PTES doprowadził do powstania szerokiego spektrum modeli różniących się zarówno stopniem odwzorowania procesów fizycznych, jak i wymaganiami obliczeniowymi. W literaturze spotykane są między innymi modele bilansowe, modele warstwowe o skupionych parametrach, modele półanalityczne, a także wysokorozdzielcze modele CFD oparte na rozwiązywaniu równań Naviera-Stokesa oraz równań transportu energii. Modele wysokiej dokładności umożliwiają szczegółową analizę lokalnych zjawisk cieplno-przepływowych, takich jak degradacja stratyfikacji, wpływ dyfuzorów, konwekcja naturalna czy oddziaływanie przepływu wód gruntowych na bilans cieplny zasobnika. Jednocześnie ich zastosowanie wiąże się z bardzo dużym kosztem obliczeniowym oraz koniecznością stosowania zaawansowanych metod dyskretyzacji przestrzennej i czasowej. Z kolei modele uproszczone pozwalają na szybką realizację analiz parametrycznych oraz wieloletnich symulacji pracy systemów energetycznych, jednak kosztem ograniczonej dokładności odwzorowania lokalnych zjawisk fizycznych. W ostatnich latach obserwowany jest intensywny rozwój metod hybrydowych oraz modeli zredukowanego rzędu (Reduced-Order Models), których celem jest osiągnięcie kompromisu pomiędzy dokładnością numeryczną a efektywnością obliczeniową. Szczególne znaczenie mają również techniki współsymulacji (co-simulation), umożliwiające sprzężenie modeli CFD z dynamicznymi modelami systemów energetycznych. Podejście takie pozwala na jednoczesną analizę lokalnych procesów zachodzących w zasobniku oraz wpływu jego pracy na funkcjonowanie całego systemu ciepłowniczego [2–8].

3. Model numeryczny A. Dahash’a

Jednym z podstawowych podejść stosowanych w analizie dynamicznych procesów cieplnych w zasobnikach energii jest modelowanie oparte na bilansie energii [1, 2]. Metody tego typu należą do grupy modeli uproszczonych o skupionych parametrach, w których przestrzeń obliczeniowa dzielona jest na ograniczoną liczbę warstw lub objętości kontrolnych o jednorodnych właściwościach termicznych. Takie podejście umożliwia opis zmian temperatury w funkcji czasu przy relatywnie niewielkich wymaganiach obliczeniowych. Modele bilansowe są wykorzystywane w analizach systemowych oraz symulacjach długookresowych, szczególnie w przypadkach wymagających wielokrotnego wykonywania obliczeń parametrycznych. Ich istotną zaletą pozostaje możliwość stosunkowo prostego odwzorowania procesów ładowania i rozładowania zasobnika przy zachowaniu akceptowalnej dokładności predykcji globalnych parametrów cieplnych. Jednocześnie modele tego typu posiadają ograniczoną zdolność odwzorowania lokalnych zjawisk przepływowych oraz szczegółowej struktury pola temperatury wewnątrz zasobnika. W analizowanym modelu numerycznym [7] przyjęto jednowymiarowy opis procesu magazynowania ciepła z wykorzystaniem warstwowego podziału objętości zasobnika, wymiana energii pomiędzy sąsiednimi warstwami opisywana jest z wykorzystaniem równań bilansu energii. Model uwzględnia podstawowe mechanizmy transportu ciepła, obejmujące przewodzenie pomiędzy warstwami, procesy związane z ładowaniem i rozładowaniem zasobnika oraz straty ciepła do otoczenia. Zastosowanie dyskretyzacji czasowej umożliwia analizę nieustalonego charakteru procesów cieplnych zachodzących w zasobniku oraz wyznaczenie ewolucji pola temperatury podczas pracy instalacji.

Rys. 1. Schemat dwuwymiarowej domeny obliczeniowej zasobnika typu PTES zastosowanej w modelu [2]
Fig. 1. Schematic diagram of the two-dimensional computational domain of the PTES type container used in the model [2]

3.1. Bilans masy

Bilans przepływu masy dla pojedynczego segmentu obliczeniowego można zapisać w postaci:

\mathrm{\Sigma}{\dot{m}}_{in}=\mathrm{\Sigma}{\dot{m}}_{out}=\rho\bullet\ \dot{V} \quad(1)

gdzie:
m_in– masa czynnika wpływająca do i-tego segmentu, kg
m_out – masa czynnika wypływająca z i-tego segmentu, kg
ρ – gęstość czynnika, kg/m³
V – strumień objętościowy czynnika, m³/h

3.2. Bilans energii

Ogólny wzór bilansu energii dla dowolnego segmentu:

\frac{\partial E_i(t)}{\partial t}=\mathrm{\Sigma}{\dot{m}}_i\bullet h_i+\left({\dot{Q}}_y-{\dot{Q}}_{y+dy}\right)-{\dot{Q}}_{str,i} \quad(2)

gdzie:
m_i h_i– strumień entalpii, kW
Q_y – Q_y+dy – przewodzenie ciepła pomiędzy segmentami, kW
Q_str,i – straty ciepła do otoczenia, kW

Po wyprowadzeniu, równanie przyjmuje postać:

(ρA_i c_p) \frac{∂T_i}{∂t}= \frac{ρc_p}{∆y_i} ((\dot{V}_y \bullet T_y )-(\dot{V}_{y+dy} \bullet T_{y+dy} ))+ \frac{ρc_p}{∆y_i} ((\dot{V}_{in,II} \bullet T_{in,II} )-(\dot{V}_{out,II}∙T_i )) +A_i \frac{∂}{∂y_i} (λ_w \frac{∂T_i}{∂y_i})-\dot{q}_{str} \quad(3)

{\dot{q}}_{str}=U_{side}\bullet P_i\bullet(T_i\left(t\right)-T_{g,i}\left(t\right)) \quad(4)

gdzie:
c_p – ciepło właściwe, kJ/kgK
λ_w– współczynnik przewodzenia ciepła, W/mK
U_side – współczynnik przenikania ciepła, W/m2K
P_i – długość strefy przenikania ciepła, m

W wielkoskalowych zasobnikach ciepła zagłębionych w gruncie najczęściej stosuje się układ dwóch dyfuzorów, odpowiadających za doprowadzanie oraz odbiór czynnika roboczego. Przy takim rozwiązaniu można przyjąć, że strumień objętościowy wpływający do zasobnika jest równy strumieniowi odpływającemu (V_in,II = V_out,II) . W okresie magazynowania energii cieplnej, gdy przepływ przez układ dyfuzorów nie występuje (V_in,II = V_out,II = 0), ruch czynnika wewnątrz zasobnika jest znacząco ograniczony. W takich warunkach dominującym mechanizmem transportu energii stają się procesy związane z oddziaływaniem sił wyporu wynikających ze zmian gęstości medium magazynującego ciepło. Istotnym zagadnieniem eksploatacyjnym pozostają straty ciepła występujące w górnej części zasobnika. Intensywna wymiana ciepła przez pokrywę może prowadzić do lokalnego obniżenia temperatury w przypowierzchniowej warstwie wody. W rezultacie może dojść do zaburzenia naturalnej stratyfikacji termicznej oraz powstania zjawiska określanego jako odwrócona termoklina [2-6].

W celu ograniczenia błędów obliczeniowych związanych z odwróconą termokliną w metodzie Dahasha stosuje się poniższy wzór, gdzie zwiększona przewodność cieplna jest uzależniona od liczny Nusselta i Rayleigha:

\lambda_i\bullet \frac{\partial}{\partial y_i} \bullet ( \lambda_w \frac{\partial T_i}{\partial y_i} ) = \left\{ \begin{array}{ll} A_i \bullet \frac{\partial}{\partial y_i} \bullet ( \lambda_w \frac{\partial T_i}{\partial y_i} ), & \dot{V}_w \neq 0, \\[6pt] A_i \bullet \frac{\partial}{\partial y_i} \bullet ( \lambda_{w,\mathrm{enh}} \bullet \frac{\partial T_i}{\partial y_i} ), & \dot{V}_w = 0 \land \dfrac{\partial T_i}{\partial y_i} < 0 \end{array} \right\} \quad(6)

\lambda_{w,enh}=\lambda_w\bullet\ Nu \quad(7)

Liczba Nusselta opisana jest z wykorzystaniem dwóch stałych C oraz k, które są wyznaczane eksperymentalnie.

Nu=C\ \bullet\ {Ra}^k \quad(8)

Ra=\frac{(gβ∙∆T∙y^3)}{να} \quad(9)

gdzie:
g- przyśpieszenie ziemskie, m/s²
β – współczynnik rozszerzalności cieplnej, 1/K
ν – współczynnik lepkości kinematycznej, m²/s
α – współczynnik dyfuzyjności cieplnej, kW/mK
ΔT – różnica temperatur sąsiadujących segmentów, K

Po przekształceniu powyższego, autorzy modelu uzyskali:

\lambda_{w,enh}=C\bullet{(\frac{\partial T_i}{\partial y_i})}^k \quad(10)

gdzie:
C – stała, która łączy różne parametry wymiarowe (np. objętość, wysokość) z właściwościami termofizycznymi (np. gęstość, pojemność cieplna właściwa, współczynnik rozszerzalności cieplnej)
k – dostosowywany w zależności od zastosowania (w opisanym w modelu przypadku k = 0,5)

Wymiana ciepła zachodząca również w gruncie jest opisana następującym równaniem:

{(\rho}_g \bullet c_{p,g}) \bullet \frac{{\partial T}_g(t)}{\partial t}=\nabla\bullet(\lambda_g \bullet\nabla T_g) \quad(11)

gdzie:
ρ_g – gęstość gruntu, kg/m³
c_p,g – ciepło właściwe gruntu przy stałym ciśnieniu, J/(kg⋅K)
T_g – temperatura gruntu, ° C
t – czas, s
λ_g – współczynnik przewodzenia ciepła dla gruntu, W/(m⋅K)

Model opracowany przez zespół A. Dahasha opiera się na jednowymiarowym podziale domeny obliczeniowej na poziome
warstwy o zmiennej szerokości wynikającej z geometrii zasobnika. Wysokość poszczególnych segmentów stanowi parametr dyskretyzacji przestrzennej i wpływa bezpośrednio na rozdzielczość numeryczną rozwiązania. W modelu uwzględniono transport energii pomiędzy sąsiednimi warstwami oraz straty ciepła do otoczenia. Istotnym problemem modelu pozostaje jednak sposób wyznaczania parametrów korelacyjnych wykorzystywanych w równaniu liczby Nusselta. W analizowanych publikacjach nie przedstawiono jednoznacznej procedury określania współczynników C oraz k, mimo że mają one bezpośredni wpływ na intensywność modelowanego transportu ciepła. Dodatkowo wartości tych parametrów różnią się pomiędzy poszczególnymi opracowaniami, co wskazuje na ich zależność od warunków geometrycznych oraz parametrów pracy analizowanego zasobnika. W praktyce oznacza to konieczność indywidualnej kalibracji modelu dla każdej analizowanej konfiguracji. Przyjęcie stałej temperatury w obrębie pojedynczej warstwy stanowi kolejne uproszczenie modelu. Podejście takie ogranicza możliwość odwzorowania lokalnych gradientów temperatury oraz szczegółowej struktury pola cieplnego, szczególnie w rejonie termokliny. Model nie umożliwia również bezpośredniej analizy lokalnych zjawisk przepływowych, ponieważ efekty konwekcyjne uwzględniane są wyłącznie pośrednio poprzez korektę efektywnego współczynnika przewodzenia ciepła. Dokładność uzyskiwanych wyników zależy przede wszystkim od przyjętego podziału warstwowego oraz jakości danych wykorzystywanych podczas kalibracji modelu, co jest wnioskiem identycznym do postawionego w [1] dotyczącego modelowania zasobników typu TES. Pomimo przyjętych uproszczeń model umożliwia efektywną analizę globalnego bilansu energetycznego zasobnika przy relatywnie niskim koszcie obliczeniowym.

4. Model półanalityczny M.Gao

Model półanalityczny stanowi rozwinięcie metod uproszczonego modelowania procesów cieplnych zachodzących w wielkoskalowych zasobnikach typu PTES. Implementacja modelu została wykonana w środowisku TRNSYS z wykorzystaniem metody różnic skończonych do dyskretyzacji równania bilansu energii. Podejście zaproponowane przez autorów opiera się na połączeniu numerycznej dyskretyzacji domeny obliczeniowej z rozwiązaniami analitycznymi wykorzystywanymi do opisu przewodzenia ciepła w gruncie otaczającym zasobnik. Zastosowanie półanalitycznej struktury modelu umożliwia ograniczenie kosztu obliczeniowego przy jednoczesnym zachowaniu zdolności odwzorowania długookresowych procesów wymiany ciepła pomiędzy zasobnikiem a otaczającym środowiskiem gruntowym.

W modelu przyjęto podział zasobnika oraz otaczających warstw gruntu na segmenty obliczeniowe, dla których wyznaczane są zmiany temperatury w kolejnych krokach czasowych. Szczególną uwagę poświęcono odwzorowaniu dyfuzji cieplnej w gruncie, mającej istotny wpływ na długoterminową sprawność magazynowania energii. Do oceny przebiegu nieustalonego transportu ciepła autorzy wykorzystali liczbę Fouriera, stanowiącą bezwymiarowy parametr opisujący relację pomiędzy czasem trwania procesu a zdolnością ośrodka do przewodzenia ciepła. Poprawność działania modelu została zweryfikowana przez autorów z wykorzystaniem danych eksploatacyjnych pochodzących z rzeczywistego zasobnika PTES zlokalizowanego w Dronninglund. Przeprowadzona walidacja umożliwiła ocenę zdolności modelu do odwzorowania zmian temperatury oraz strat cieplnych występujących podczas długookresowej pracy instalacji.

Rys. 2. Podział na regiony wymiany ciepła [7]
Fig. 2. Heat exchange regions [7]

Rys. 3. Oznaczenia węzłów w zależności od pionowego i poziomego położenia względem węzła j [7]
Fig. 3. Node designations depending on the vertical and horizontal position relative to node j [7]

Uproszczenia zastosowane w modelu:

stałe właściwości fizyczne wody oraz gruntu,
w gruncie uwzględnia się wyłącznie przewodzenie ciepła, wpływ przepływu wód gruntowych i naturalnej konwekcji jest pomijany,
pomijane jest konwekcyjne przenoszenie ciepła wewnątrz zasobnika, współczynnik mieszania inwersyjnego uproszczono do stałej,
założono magazyn i obszar gleby jako symetryczne, ma to na celu pominięcie efektu anizotropii,
lepkie rozpraszanie się płynu jest pomijane,
temperatura początkowa węzła jest stałą uśrednioną dla jednego kroku czasowego.

4.1. Wyznaczenie współczynnika przenikania ciepła

Ze względu na brak izolacji pomiędzy zasobnikiem a glebą,
współczynnik przenikania ciepła jest określany na podstawie konwekcji po stronie wody oraz przewodzenia po stronie gleby. Pominięto konwekcję wymuszoną.

U=\frac{A_{mea}}{A_{sim}}\bullet\frac{1}{\frac{1}{h_w}+\frac{\delta HDPE}{\lambda HDPE}+\frac{\delta geo}{\lambda geo}} \quad(12)

gdzie:
U – współczynnik przenikania ciepła, W/m²K
A_mea, A_sim – powierzchnia rzeczywista oraz symulowana, m²
h_w – współczynnik przejmowania ciepła po stronie wody, W/m²K
δHDPE, δgeo – grubość warstw HDPE oraz geowłókniny, m
λHDPE, λgeo – współczynniki przewodzenia cieplnego HDPE i geowłókniny, W/mK

h_w=\ \frac{0.6\lambda_w}{H}\bullet{(Gr\bullet P r)}^\frac{1}{3} \quad(13)

Gr=\ \frac{g\bullet\beta_w\bullet(T_{geo}-T_w)\bullet D^3}{64\bullet\nu_w^2} \quad(14)

gdzie:
h_w – współczynnik przejmowania ciepła po stronie wody, Wm²K
λ_w – współczynniki przewodzenia cieplnego wody, W/mK
H – wysokość ściany bocznej, m
Gr – liczba Grashofa
Pr – liczba Prandtla
g – współczynnik przyśpieszenia ziemskiego, m/s²
B_w – współczynnik rozszerzalności objętościowej, 1/K
T_geo – temperatura geowłókniny, °C
T_w – temperatura wody w zasobniku, °C
D – ekwiwalentna średnica ściany bocznej, m
ν – współczynnik kinematycznej lepkości wody, m²/s

4.2. Wymiana ciepła w regionie magazynowania wody (Region V)

Region V w zasobniku modelowany jest za pomocą jednowymiarowego równania przewodzenia ciepła. Założono, że równanie przewodzenia zostanie zapisane w jawnej centralnej postaci różnicowej jak przedstawiono poniżej. Dodatkowo składnik nieustalony jest pierwszego rzędu, a składnik dyfuzyjny drugiego rzędu.

\frac{\rho_wC_pV_j(T_j^{t+1}-T_j^t)}{\mathrm{\Delta t}}=\ K_{j+1}\left(T_{j+1}^t-T_j^t\right)+K_{j-1}\left(T_{j-1}^t-T_j^t\right)+U_BA_j\left(T_B^t-T_j^t\right)+S_T \quad(15)

\rho_w– gęstość wody, kg/m³
C_p – ciepło właściwe wody, J/kgK
V_j – objętość przypisana do j-tego węzła, m³
T_j^t+1 – temperatura węzła j w chwili t+1, K
Δt – krok czasowy, s
K_j+1 – współczynnik przewodzenia ciepła pomiędzy węzłami j oraz j+1, W/K
U_B – współczynnik konwekcji na powierzchni granicznej, W/m²K
A_j – powierzchnia wymiany ciepła w węźle j, m²
T_B^t – temperatura warstwy granicznej, K
S_T – źródło ciepła (dopływ gorącej wody), W

W modelu przedstawione zostały równania służące do wyznaczenia A_j, K_j+1 oraz S_T.

A_{side,j}=\frac{1}{2} π∙∆r_j∙√(∆y_j^2+ ∆r_j^2 ) \quad(16)

K_{j±1}=\frac{λw(∆r_j+∆r_(j±1)²}{2(∆y_j+∆y_(j±1)}+m_j C_p+m_{mix,j} C_p \quad(17)

S_T=V_{in} ρ_w C_p (T_{in}-T_j) \quad(18)

Zaproponowano przybliżone rozwiązanie analityczne, które oferuje kilka korzyści w porównaniu z metodą numeryczną. Model rozwiązuje problem matematyczny, eliminując konieczność opierania się na niestandardowym sformułowaniu numerycznym. Ponadto rozwiązanie analityczne nie zależy od kroku czasowego, choć wymaga iteracyjnego rozwiązywania w podprogramie w celu rozwiązania równań różniczkowych. Rozwiązanie przyjmuje formę:

\frac{dT}{dt}=aT_j+b \quad(19)

gdzie:

a= \frac{-(K_{j+1}+K_{j-1}+U_B A_j+V_{in} ρ_w C_p)}{V_j ρ_w C_p} \quad(20)

b=\frac{(K_{j+1} T_{j+1}+K_{j-1} T_{j-1}+U_B A_j T_B+V_{in} ρ_w C_p T_{in})}{V_j ρ_w C_p} \quad(21)

T_j^{t+1}=(T_j^t+\frac{b}{a})∙e^a∆t-\frac{b}{a} \quad(22)

Podstawową zaletą rozwiązania analitycznego jest niezależność od kroku czasowego. W przeciwieństwie do standardowych metod numerycznych, które wymagają precyzyjnie dobranego kroku czasowego w celu zachowania stabilności, rozwiązanie analityczne eliminuje ten problem. Niemniej jednak, konieczne jest iteracyjne rozwiązywanie równań różniczkowych w podprogramie, aby zapewnić poprawność wyników. Rozwiązanie tego równania różniczkowego jest bezpośrednie w przypadku, gdy 𝑏 jest stałą. Jeśli jednak 𝑏 zmienia się w czasie, problem można rozwiązać przybliżeniowo, zakładając, że wartość 𝑏 jest stała w danym kroku czasowym i równa jej wartości średniej. Choć założenie stałości 𝑏 w czasie nie jest w pełni z rzeczywistością, stanowi ono aproksymację dla małych kroków czasowych. W badanym przypadku, gdzie maksymalny krok czasowy wynosił 1 godzinę, podejście to zapewniło wystarczającą dokładność w symulacji dużych objętości PTES.

4.3. Wymiana ciepła w obszarze gruntu i izolacji (Region VI)

Grunt i izolacja są modelowane jako dwuwymiarowe siatki obliczeniowe, co oznacza, że wymiana ciepła jest opisana za pomocą dwuwymiarowego równania przewodnictwa cieplnego. Do rozwiązania tego równania zastosowano metodę różnic centralnych, w której i k odpowiadają numerom węzłów w kierunku promieniowym i pionowym.

Podstawowe równanie bilansu energii dla każdego węzła siatki numerycznej przyjmuje postać:

V_{i,k} ρ_g C_g \frac{T_{i,k}^{t+1}-T_{i,k}^t}{∆t}=K_y^{i,k+1} (T_{i,k+1}^t-T_{i,k}^t)+K_y^{i,k-1} (T_{i,k-1}^t-T_{i,k}^t )+K_r^{i+1,k} (T_{i+1,k}^t-T_{i,k}^t )+ K_r^{i-1,k} (T_{i-1,k}^t-T_{i,k}^t )+Q_{i,k} \quad(23)

Wzory na K i Q_i, k różnią się w zależności od regionu wymiany ciepła. Region VI jest fizycznym regionem izolacji, a ogólny opór przewodzenia ciepła jest przedstawiony jako:

K_{i\pm1,k}^r=\frac{( ∆y_{k+1}+∆y_k)\bullet\lambda_{ins/g}}{2lj(\frac{∆r_{i+1}}{∆r_i})} \quad(24)

K_{i\pm1,k}^z=\frac{( ∆r_{i±1}+∆r_i)\bullet\lambda_{ins/g}\bullet∆r_i}{(∆y_{k±1}+{∆y})} \quad(25)

Wymianę ciepła pomiędzy górną granicą izolacji a otoczeniem (Region VI) opisać można jako:

Q_{i,1}=h_{am} \bullet A_{i,1} \bullet (T_{am}-T_{i,1}) \quad(26)

Dodatkowo, proponowane jest rozpatrzenie wymiany ciepła w kierunkach pionowym oraz poziomym względem węzła, w którym wykonywane są obliczenia. Wzór na konwekcję pomiędzy górną powierzchnią zasobnika PTES, a otoczeniem zasobnika ma następującą postać:

Q_{i,NY}=U_{top} \bullet A_{i,NY} \bullet (T_{top}-T_{i,NY}) \quad(27)

gdzie:
U_top – współczynnik przejmowania ciepła dla górnej powierzchni zasobnika, W/(m²⋅K)
A_i,NY – pole powierzchni górnej warstwy zasobnika w komórce siatki i na poziomie NY, m²
T_top – temperatura otoczenia nad zasobnikiem ciepła, °C
T_i,NY – temperatura górnej powierzchni zasobnika w węźle i NY, °C

Wzór na wymianę ciepła pomiędzy ostatnią warstwą izolacji pokrywy zasobnika, na krawędzi zasobnika (skrajne poziome położenie), a otoczeniem zasobnika ma następującą formę:

Q_{NR,k} = \left( \frac{\log\left(\dfrac{\Delta r_{NR} + \Delta r_{NR+1}}{2\Delta r_{NR}}\right)} {2(\Delta y_{k+1} + \Delta y_k)\,\lambda_{\mathrm{ins}}} + \frac{\log\left(\dfrac{2\Delta r_{NR}}{\Delta r_{NR} + \Delta r_{NR+1}}\right)} {2(\Delta y_{k+1} + \Delta y_k)\,\lambda_{s}} \right)^{-1} \cdot A_{NR,k}\,(T_{NR+1,k} – T_{NR,k}) \quad(28)

Powyższy wzór bierze pod uwagę szereg połączonych oporów przewodzenia dla izolacji oraz gleby. Dokładne opisanie strat
ciepła poprzez region VI jest kluczowe dla obliczeń oraz prowadzonych na podstawie modelu symulacji, ponieważ to pokrywą zasobnika tracona jest największa ilość ciepła.

4.4. Wymiana ciepła w gruncie

Zgodnie z rys. 2 i 3 wymiana ciepła w gruncie przebiega w 4 regionach (I, II, III oraz IV).

Region IV – kontakt gleby z powietrzem K_NR,k^r=0; K_R,k^r=0, założono brak przepływu ciepła w kierunku poziomym (promieniowym). K_i,1^y=0, założono brak przepływu ciepła w kierunku pionowym.

-\lambda_s\frac{\partial T}{\partial n}=h_s(T-T_\infty)

z warunków brzegowych Robina wyprowadzono wzór na przenikanie ciepła z gruntu do otoczenia.

Region II – gleba bez bezpośredniego styku z zasobnikiem ciepła. K_NR,k^r=0; K_R,k^r=0, założono brak przepływu ciepła w kierunku poziomym. K_i,1^y=0, założono brak przepływu ciepła w kierunku pionowym. T=T_far z warunków brzegowych Direchleta, temperatura w głębi gruntu, bez punktu styku z zasobnikiem jest stała. Region II jest termicznie stabilny, wahania temperatur są pomijalnie małe.

Region I – grunt, na który oddziałuje dno zasobnika

-\lambda_f\left(\frac{\partial T}{\partial n}\right)_f=\frac{T_f{-T}_s}{\frac{\delta}{k}}=\lambda_s\left(\frac{\partial T}{\partial n}\right)_s

wzór na przewodzenie ciepła z dolnej warstwy zasobnika do gruntu. Warunek brzegowy Robina zapewnia uwzględnienie ciągłości strumienia ciepła.

Region III – grunt przy bocznej ścianie zasobnika. Obliczenie współczynników K_i±1,k^r; K_i±1,k^z odbywa się ze wzorów (24) oraz (25) wykorzystywanych w przypadku regionu VI. Ze względu na kształt regionu zmienia się równanie bilansu energii. Przyjmuje ono poniższą postać:

V_{i,k} \bullet \rho_g \bullet C_g \frac{T_{i,k}^{t+1}-T_{i,k}^t}{2∆t} = K_y^{i,k+1}\left(T_{i,k+1}^t-T_{i,k}^t\right)+K_r^{i+1,k}\ \left(T_{i+1,k}^t-T_{i,k}^t\right)+Q_{i,k}

Zastosowano tak zwaną różnicę w przód. Uwzględniono strumień cieplny zarówno w pionie, jak i poziomie oraz lokalną stratę ciepła. W celu odwzorowania wymiany energii pomiędzy poszczególnymi regionami modelu, w równaniu podstawowym uwzględniono strumienie ciepła opisujące transport energii pomiędzy sąsiadującymi obszarami obliczeniowymi. Autorzy zaprezentowali odpowiednie zależności matematyczne definiujące poszczególne składowe wymiany ciepła.

Strumień ciepła pomiędzy gruntem, a otoczeniem:

Q_{i,1}=h_{am}A_{i,1}\left(T_{am}-T_{i,1}\right),\ dla\ i=(1,R) \quad(29)

Powyższy wzór określa wymianę ciepła pomiędzy warstwą gruntu, a powietrzem atmosferycznym, wykorzystano warunek brzegowy Robina, gdzie:
h_am – współczynnik przejmowania ciepła między glebą a powietrzem, W/(m²⋅K)
A_i,1 – powierzchnia wymiany ciepła dla danego węzła siatki, m²
T_am – temperatura otoczenia, °C
T_i,1 – temperatura wierzchniej warstwy gruntu, °C

Strumień ciepła pomiędzy regionem III, a regionem II: Określono strumień ciepła wymieniany pomiędzy boczną częścią (region III) a strefą dalekiego zasięgu (far field – region II):

Q_{i,y}= \frac{2\lambda_sA_{i,k}}{∆ y_Y+∆ y_{Y-1}} \bullet A_{R,k}\bullet(T_{far}-T_{R,k}),\quad dla \space k=(1,Y) \quad(30)

gdzie:
𝜆𝑠 – współczynnik przewodzenia cieplnego gruntu, W/(m⋅K)
r_far – odległość w osi poziomej, gdzie zakłada się stałą temperaturę T_far, m
Δr_R – odległość pozioma (radialna) od centrum zasobnika do węzła R, m

Wzór uwzględnia cylindryczną geometrię zasobnika oraz logarytmiczny charakter przewodzenia ciepła w takiej geometrii. Strumień ciepła w głąb gruntu:

Q_{i,y}= \frac{2\lambda_sA_{i,k}}{∆ y_Y+∆ y_{Y-1}} \bullet (T_{deep}-T_{i,NY}) \quad(31)

gdzie:
T_deep – temperatura gruntu poniżej zasobnika, zgodnie z warunkiem brzegowym Direchleta, temperatura w głębi gruntu jest stała, °C
T_i,NY – temperatura gruntu przy zasobniku, °C
Δy_Y – wysokość siatki w pionie, m

Strumień ciepła z dna zasobnika do gruntu:

Q_{i,Ny}=U_{bot}A_{i,Ny}(T_{bot}-T_{i,Ny}) \quad(32)

gdzie:
U_bot – współczynnik przejmowania ciepła na dnie zasobnika (obejmuje konwekcję w wodzie, przewodnictwo materiału oddzielającego wode z zasobnika od gruntu, przewodnictwo geowłókniny oraz przewodnictwo gruntu), W/(m²⋅K)
T_bot – temperatura wody na dnie zasobnika, °C
T_i,NY – temperatura gruntu przy zasobniku, °C

Wzór uwzględnia sprzężenie ciecz-ciało stałe. Strumień ciepła ze ścian bocznych do gruntu:

Q_{i,k}=U_{side}A_{i,k}(T_{side,i,k}-T_{i,k}) \quad(33)

gdzie:
U_side – współczynnik przejmowania ciepła dla bocznych ścian zasobnika, W/(m²⋅K)
T_side,i,k – temperatura na granicy zasobnika i gruntu, °C
T_i,k – temperatura gruntu w badanym węźle, °C

Powierzchnia wymiany ciepła A_ik:

A_{i,k}= \frac{1}{4} π \bullet (∆r_{i+1,k}+∆r_{i,k})^2 \quad(34)

Dla zasobnika o innej geometrii należy zastosować wzór na powierzchnię odpowiedniej bryły.

4.5. Bilans energii

W celu przeprowadzenia bilansu energii autorzy opracowali zestaw zależności opisujących wymianę energii pomiędzy magazynem a otaczającym środowiskiem. Zaproponowane równania umożliwiają wyznaczenie energii dostarczanej i odbieranej z zasobnika, określenie strat ciepła oraz analizę zmian poziomu akumulacji energii w czasie. Pierwszym wzorem określa się ilość energii przekazywaną lub odbieraną z magazynu ciepła w czasie jednego kroku czasowego Δτ:

Q_{ch/dis}=\rho_wc_p\Delta\tau(V_{out}T_{out}-V_{in}T_{in}) \quad(35)

gdzie:
𝜌_𝑤 – gęstość wody, kg/m3
c_p – ciepło właściwe wody, J/kgK
Δτ – krok czasowy, s
V_outT_out – objętość i temperatura wody wypływającej z PTES, m³/s, K
V_inT_in – objętość i temperatura wody wpływającej do PTES, m³/s, K

Następnym elementem są straty ciepła z zasobnika zagłębionego w gruncie typu PTES:

Q_{loss}=Q_{loss,top}+Q_{loss,side}+Q_{loss,bot} \quad(36)

gdzie:
Q_loss,top – straty ciepła przez pokrywę zasobnika, W
Q_loss,side – straty ciepła przez powierzchnię boczną zasobnika, W
Q_loss,bot – straty ciepła przez dno zasobnika, W

Energia zmagazynowana, w postaci energii wewnętrznej cieczy magazynującej ciepło (wody), w danym momencie opisana jest wzorem:

Q_{sto}=\sum_{j=1}^{N}{\rho_wc_pV_jT_j} \quad(37)

gdzie:
N – liczba warstw (węzłów) w siatce zasobnika ciepła,
𝜌_𝑤 – gęstość wody, kg/m³
c_p – ciepło właściwe wody, J/kgK
V_j – objętość przypisana do j-tego węzła, m³
T_j – temperatura j-tego węzłą, K

Zmiana energii względem stanu początkowego:

Q_{change}=\sum_{j=1}^{N}{\rho_wc_pV_j{(T}_{j,now}}-T_{j,start}) \quad(38)

gdzie:
T_j,now – aktualna temperatura węzła j, °C
T_j,start – temperatura początkowa węzła j, °C

Opisany model półanalityczny stanowi zaawansowane narzędzie wykorzystywane do symulacji procesów termodynamicznych zachodzących w zasobnikach typu PTES. Wykorzystuje on dwuwymiarową dyskretyzację domeny obliczeniowej w układzie cylindrycznym oraz umożliwia analizę wymiany ciepła zarówno w obrębie zasobnika, jak i otaczającego gruntu. Zaproponowane podejście umożliwia analizę zmian temperatury gruntu otaczającego zasobnika oraz wyliczenie sprawności magazynowania ciepła, która będzie zmieniać się w czasie wraz z wygrzewaniem się otaczającej gleby. Do podstawowych zalet modelu należy wysoka dokładność odwzorowania modelu, zweryfikowana dokładnymi danymi pomiarowymi zasobnika w Dronninglund. Model potwierdza, że zastosowanie uproszczenia w postaci dwuwymiarowej geometrii umożliwia uzyskanie wysokiej dokładności przy zmniejszonych nakładach obliczeniowych. Ograniczeniem modelu pozostaje pominięcie bezpośredniego modelowania procesów konwekcyjnych w objętości magazynowanej wody. Powoduje to ograniczoną możliwość odwzorowania lokalnych zjawisk przepływowych i szczegółowej struktury pola temperatury podczas dynamicznych procesów ładowania i rozładowania zasobnika. Ze względu na skomplikowany charakter modelu, zaleca się stosowanie metody do procesów optymalizacji pracy istniejącego system. Umożliwia to weryfikację oraz walidację modelu dzięki danym pomiarowym. W zakresie projektowania, brak szczegółowych danych wynikających z pomiarów może generować niedokładności w obliczeniach oraz brak możliwości weryfikacji wyników symulacji.

5. Model numeryczny A.Tosatto

Model numeryczny wielkoskalowego zasobnika do sezonowego magazynowania ciepła, opracowany w środowisku COMSOL Multiphysics. Z założenia model ma umożliwić szybką ocenę efektywności zastosowania wielkoskalowego zasobnika ciepła zagłębionego w gruncie. Opracowanie opiera się na uproszczonym modelu uwzględniającym straty ciepła, przewodzenie w gruncie oraz wpływ sił wyporu. Model można zastosować przy obliczaniach dotyczących zasobników zagłębionych w gruncie typu PTES oraz TTES zagłębionych w gruncie.

Rys. 4. Poglądowy model 2-D opisywanego w [8] zasobnika ciepła
Fig. 4. Illustrative 2-D model of the heat accumulator described in [8]

Założenia przedstawionego modelu:

w pełni zakopany cylindryczny zasobnik ciepła, modelowany jako 2D osiowo-symetryczny,
zasobnik podzielony został na równe warstwy wzdłuż osi pionowej, które wykorzystywane są do obliczeń bilansu energii,
grunt wokół zasobnika modelowany w 2D (wykorzystano moduł Heat Transfer in Solids w programie Comsol),
obudowa zasobnika może być modelowana z wykorzystaniem dowolnych materiałów, umożliwia to zarówno symulację dla izolowanych zasobników o betonowych konstrukcjach, jak i zasobników, których ściany boczne wyłożone są materiałem oddzielającym zasobnik od gruntu bez izolacji,
pełnym cyklem pracy zasobnika uznaje się cztery fazy (ładowanie, magazynowanie, rozładowywanie oraz bezczynność). Temperaturę otoczenia można modelować dowolnie.

5.1. Równania opisujące model

Sprawność zasobnika 𝜂_TES:

\eta_{TES}=1-\frac{Q_{loss}}{E_{max}} \quad(39)

gdzie:
Q_loss – straty ciepła, kWh
E_max – zmagazynowana energia cieplna, kWh

Bilans energii dla pojedynczej warstwy:

\frac{\partial E_i}{\partial t}=\dot{m}\bullet{(h}_z-h_{z+dz})+\left({\dot{q}}_z-{\dot{q}}_{z+dz}\right)-{\dot{Q}}_{loss,side} \quad(40)

gdzie:
m (h_z-h_z+dz) – strumień entalpii, W
q_z-q_z+dz – przewodzenie ciepła pomiędzy segmentami, W
Q_loss,side – straty ciepła do otoczenia, W

Wymiana ciepła pomiędzy warstwami:

(\rho A{)c}_p\frac{\partial T_i}{\partial t}=-\rho Avc_p\left(\frac{\partial T}{\partial z}\right)-A\frac{\partial}{\partial z}\left(\lambda\frac{\partial T}{\partial z}\right)-\pi DU(T-T_{ext}) \quad(41)

gdzie:
ρ – gęstość czynnika, kg/m³
A – pole przekroju każdej warstwy, m²
c_p – ciepło właściwe, KJ/kgK
v – prędkość przepływu, m/s
π_D – obwód warstwy, m
U – współczynnik przenikania ciepła, W/m²K
T–T_ext – różnica temperatury pomiędzy wodą a otoczeniem, K

W modelu zastosowano efektywny współczynnik przewodzenia ciepła, którego zadaniem jest pośrednie uwzględnienie wpływu naturalnej konwekcji na transport energii w objętości magazynowanej wody. Wartość współczynnika określana jest jako funkcja pionowego gradientu temperatury, co pozwala na częściowe odwzorowanie intensyfikacji wymiany ciepła wywołanej siłami wyporu. Podejście takie umożliwia uwzględnienie efektów konwekcyjnych bez konieczności rozwiązywania pełnego układu równań przepływu płynu, jednak wymaga stosowania parametrów empirycznych, których procedura wyznaczania nie została szczegółowo przedstawiona przez autorów modelu. Warunkiem jego wykorzystania jest spełnienie warunku 𝜕T/𝜕z < 0 dla wymiany ciepła pomiędzy warstwami. W takim przypadku drugi człon równania przyjmuje postać:

-A\frac{\partial}{\partial z}\left((\lambda+\lambda_{w,enh})\frac{\partial T}{\partial z}\right) \quad(42)

Współczynnik przewodności o zwiększonej wydajności oblicza się na podstawie wzoru:

\lambda_{w,enh}=c\sqrt{{(\frac{\partial T}{\partial z})}^n} \quad(43)

Grunt otaczający zasobnik ciepła modelowany jest z wykorzystaniem Heat Transfer in Solids, pozwalającym na uwzględnienie wymiany ciepła przez przewodzenie. Funkcją wyjściową dla tego zjawiska jest równanie Fouriera:

\rho C_p\frac{dT}{dt}=\nabla\bullet\lambda\nabla T \quad(44)

Powyższy wzór pozwala obliczyć, jak straty ciepła z zasobnika ciepła wpłyną na grunt otaczający.

Przedstawiony model numeryczny stanowi uproszczone narzędzie do analizy stanów nieustalonych procesów cieplnych zachodzących w zasobnikach typu TES/PTES. Model opiera się na sprzężeniu jednowymiarowego opisu rozkładu temperatury w objętości magazynowanej wody z dwuwymiarowym modelem przewodzenia ciepła w gruncie, wykorzystującym założenie symetrii osiowej. Ze względu na umiarkowane wymagania obliczeniowe może być wykorzystywany do analiz parametrycznych oraz oceny wariantów projektowych na etapie koncepcyjnym. Ograniczeniem modelu pozostaje przede wszystkim uproszczony sposób odwzorowania transportu energii wewnątrz objętości magazynowanej wody. Zastosowanie jednowymiarowego opisu rozkładu temperatury nie pozwala na analizę lokalnych niejednorodności pola termicznego ani efektów mieszania występujących w pobliżu dyfuzorów. Model nie uwzględnia również przestrzennego rozkładu prędkości przepływu, przez co nie jest możliwe odwzorowanie lokalnych zjawisk hydrodynamicznych wpływających na degradację stratyfikacji termicznej. W efekcie analiza wyników otrzymanych z opisanego modelu uniemożliwia np. Porównanie wpływu smukłości zasobnika na jego sprawność. Pomimo tego, model pozwala skutecznie analizować warianty koncepcyjne na etapie projektowania systemu magazynowania ciepła.

6. Wnioski

Przeprowadzona analiza modeli numerycznych stosowanych do symulacji zasobników ciepła wykazała, że poszczególne podejścia obliczeniowe różnią się zarówno zakresem odwzorowania procesów fizycznych, jak i stopniem złożoności numerycznej oraz wymaganiami obliczeniowymi. Wybór odpowiedniego modelu powinien być uzależniony przede wszystkim od celu analizy, oczekiwanej dokładności oraz dostępności danych wejściowych. Analizowane modele uproszczone i półanalityczne umożliwiają efektywną symulację długookresowych procesów magazynowania energii. Zastosowanie dyskretyzacji warstwowej oraz uproszczonych metod odwzorowania konwekcji pozwala na analizę globalnego bilansu energii zasobnika, ocenę strat ciepła oraz wyznaczanie zmian temperatury w czasie. Podejście takie jest szczególnie przydatne w analizach parametrycznych oraz na etapie koncepcyjnym projektowania instalacji. Jednocześnie stwierdzono, że uproszczone modele numeryczne posiadają ograniczoną zdolność odwzorowania lokalnych zjawisk cieplno-przepływowych zachodzących wewnątrz zasobnika. W szczególności dotyczy to procesów związanych z naturalną konwekcją, degradacją stratyfikacji termicznej oraz lokalnymi efektami hydrodynamicznymi występującymi w rejonie dyfuzorów. W większości analizowanych modeli efekty te uwzględniane są pośrednio poprzez modyfikację efektywnych parametrów transportu ciepła, bez bezpośredniego rozwiązywania równań ruchu płynu. Istotnym problemem pozostaje również konieczność stosowania parametrów empirycznych wykorzystywanych do kalibracji modeli uproszczonych. Parametry te wykazują zależność od geometrii zasobnika, warunków eksploatacyjnych oraz właściwości materiałowych. Modele wykorzystujące sprzężenie jednowymiarowego opisu zasobnika z dwuwymiarowym modelem przewodzenia ciepła w gruncie pozwalają na bardziej dokładną analizę strat energii do otoczenia oraz długoterminowego oddziaływania cieplnego na grunt. Jednocześnie zachowują one znacznie niższe wymagania obliczeniowe niż pełne modele CFD, co umożliwia wykonywanie wieloletnich symulacji pracy systemu. Na podstawie przeprowadzonej analizy można stwierdzić, że modele uproszczone i półanalityczne stanowią efektywne narzędzie do analiz systemowych, oceny wariantów projektowych oraz szacowania długoterminowej efektywności energetycznej zasobników ciepła. Przeprowadzona analiza potwierdza, że nie istnieje uniwersalny model numeryczny odpowiedni dla wszystkich zastosowań związanych z symulacją zasobników PTES.

7. Literatura

[1] Osiadacz, A. J., Kotyński, Ł., & Kwestarz, M. (2023). A comparative study of one-dimensional models for stratified thermal energy storage. Gaz Woda i Technika Sanitarna, (97)1, 10-16, https://doi.org/10.15199/17.2023.9.2

[2] Dahash, A., Ochs, F., Tosatto, A., & Streicher, W. (2020). Toward efficient numerical modeling and analysis of large-scale thermal energy storage for renewable district heating. Applied Energy, 279, 115840. https://doi.org/10.1016/j.apenergy.2020.115840

[3] Chang, C., Nie, B., Leng, G., Li, C., She, X., Peng, X., & Deng, J. (2017). Influences of the key characteristic parameters on the thermal performance of a water pit seasonal thermal storage. Energy Procedia, 142, 495–500. https://doi.org/10.1016/j.egypro.2017.12.077

[4] Dahash, A., Janetti, M. M. B., & Ochs, F. (2020). Numerical analysis and evaluation of Large-Scale hot water tanks and pits in district heating systems. Building Simulation Conference Proceedings, 16, 1692–1699. https://doi.org/10.268 68/25222708.2019.210566

[5] Dahash, A., Ochs, F., Tosatto, A., & Streicher, W. (2019). Co-simulation of CFD and system models for large-scale thermal energy storage applications. In Proceedings of Building Simulation 2019 Conference. International Building Performance Simulation Association (IBPSA). Rome, Italy

[6] Xiang, Y., Xie, Z., Furbo, S., Wang, D., Gao, M., & Fan, J. (2022). A comprehensive review on pit thermal energy storage: Technical elements, numerical approaches and recent applications. Journal of Energy Storage, 55, 105716. https://doi.org/10.1016/j.est.2022.105716

[7] Gao, M., Shao, S., Xiang, Y., Wang, D., Furbo, S., & Fan, J. (2024). Semi-analytical model of a large-scale water pit heat storage for the long-term thermal applications. Energy and Buildings, 307, 113975. https://doi.org/10.1016/j.enbuild.2024.113975

[8] Tosatto, A., Dahash, A., Ochs, F., & Janetti, M. B. (2019). Development of a numerical model for large-scale seasonal thermal energy storage for district heating systems. In COMSOL Conference (Vol. 2019, pp. 1-7).