This benchmark is testing the extended version of a classical heat conduction equation termed ‘heat conduction equation with phase change’ (with a slight abuse of notations, we also call it simply ‘T+freezing’ equation). The initialboundary value problem (IBVP) for this equation models such processes as ice formation and ice melting in watersaturated porous medium. Since the equation is strongly nonlinear in the temperature variable $T$ to be solved for and contains multiple parameters which may affect accuracy of finite element discretization, a carefully designed model and code verification must be performed.
Below, we model heat transfer process – focusing specifically on ice formation – in a cylindrical soil specimen around a borehole heat exchanger (BHE) which contains a refrigerant of subzero temperature. This (negative) temperature is used to prescribe a Dirichlet boundary condition on the specimen boundary adjacent to the BHE, what triggers cooling and consequent freezing of watersaturated soil whose initial temperature is positive.
Simulations are performed using both our OpenGeoSys platform and the FreeFem++ open source finite element code (in the following, simply termed OGS and FF++, respectively), thus enabling crossverification of the numerical codes.
The detailed IBVP problem description for the T+freezing equation, geometric setup, material and model parameters used in the implementation can be found in this document this PDF document. The figures below are taken from this documentation and serve for illustrative purposes to give a hint about the modeled process and simulations outcome.
Note that $(r,z)\in S$ are denoted as $(x,y)$ which are assumed to be unrelated to the coordinate notations in the original 3d formulation. 2. The initial condition for $T$ in $S$ is assumed to be a positive function which decays linearly from surface to bottom. For modeling the (timedependent) boundary conditions on $\Gamma_D$ of $S$, it is assumed that within the first $\widehat{t}$ hours, the temperature on $\Gamma_D$ drops continuously from the initial state to the values prescribed by some continuous piecewise linear function of $y$ and such that at the last depth segment it becomes negative. (The latter mimics the impact of the BHE refrigerant with subzero temperature.) The figure sketches the situation:
Temperature is given in degrees Celsius. For $t>\widehat{t}$, the prescribed temperature on $\Gamma_D$ and, thus, the heat conduction in the modelled case is triggered by a significant difference between the temperature on $\Gamma_D$ and the initial one within $S$. 3. The results of modelling are depicted in the following two figures, where we plot the temperature distribution in the soil block after 720 hours (30 days) of cooling, and also compare the outcomes of the two corresponding packages:
Temperature is given in kelvins. The color legend of $T$ in the corresponding ParaView plots is tuned such that the amount of ice formed around BHEs can be identified. As expected, ice formation occurs in the vicinity of $\Gamma_D$, more specifically, near the segment of $\Gamma_D$ in which the negative temperature has been prescribed. In the rest of the domain, temperature distribution remains almost identical to the initial state, as could also be expected. 4. Finally, the corresponding results from the previous figure are plotted over the three different directed lines within the domain $S$:
Here, origin of the horizontal axis on the right plot corresponds to line’s origin. For the selected lines, the compared data seems identical pointwise, thus supporting the quantitative similarity of the OGS and FF++ results observed earlier.
In the corresponding OGS project file m16m15projectB.prj
the time discretization is different for the “real case study” whose results are presented in the documentation and for the ctest
case, and must be altered manually.
This article was written by Tymofiy Gerasimov. If you are missing something or you find an error please let us know.
Generated with Hugo 0.122.0
in CI job 449919

Last revision: April 23, 2024
Commit: [PL/LD] Use generic cell average output 3557e29
 Edit this page on