Skip to content
Snippets Groups Projects
Unverified Commit d1b043c9 authored by Dmitri Naumov's avatar Dmitri Naumov Committed by GitHub
Browse files

Merge pull request #2084 from HBShaoUFZ/bhe_array_benchmark 

add a new benchmark "BHE Array" to the web content.
parents b2fc6144 ab1f59fd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 829 additions and 1 deletion
ProcessLib/HeatConduction/Tests.cmake
+ 13
− 0
View file @ d1b043c9
...@@ -42,3 +42,16 @@ AddTest( ...@@ -42,3 +42,16 @@ AddTest(
# EXECUTABLE ogs # EXECUTABLE ogs
# EXECUTABLE_ARGS wedge_1e2_axi_ang_0.02.prj # EXECUTABLE_ARGS wedge_1e2_axi_ang_0.02.prj
# ) # )
# The 25 BHE array benchmark
# test results are compared to 2D simulation result
AddTest(
NAME BHE_Array_2D
PATH Parabolic/T/2D_BHE_array
EXECUTABLE ogs
EXECUTABLE_ARGS bhe2d.prj
TESTER vtkdiff
DIFF_DATA
standard_solution_bhe2d_pcs_0_ts_840_t_72576000.000000.vtu bhe2d_pcs_0_ts_840_t_72576000.000000.vtu temperature temperature 1e-12 0.0
REQUIREMENTS NOT OGS_USE_MPI
)
Tests/Data/Parabolic/T/2D_BHE_array/bhe2d.gml 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
Source diff could not be displayed: it is stored in LFS. Options to address this: view the blob.
Tests/Data/Parabolic/T/2D_BHE_array/bhe2d.prj 0 → 100644
+ 302
− 0
View file @ d1b043c9
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<mesh>bhe2d.vtu</mesh>
<geometry>bhe2d.gml</geometry>
<processes>
<process>
<name>HeatConduction</name>
<type>HEAT_CONDUCTION</type>
<integration_order>2</integration_order>
<thermal_conductivity>lambda</thermal_conductivity>
<heat_capacity>c_p</heat_capacity>
<density>rho</density>
<process_variables>
<process_variable>temperature</process_variable>
</process_variables>
<secondary_variables>
<secondary_variable type="static" internal_name="heat_flux_x" output_name="heat_flux_x"/>
</secondary_variables>
</process>
</processes>
<time_loop>
<processes>
<process ref="HeatConduction">
<nonlinear_solver>basic_picard</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<abstol>1.e-6</abstol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<output>
<variables>
<variable> temperature </variable>
<variable> heat_flux_x </variable>
</variables>
</output>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial> 0.0 </t_initial>
<t_end> 93312000 </t_end>
<timesteps>
<pair>
<repeat>1</repeat>
<delta_t>86400</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<type>VTK</type>
<prefix>bhe2d</prefix>
<timesteps>
<pair>
<repeat> 8 </repeat>
<each_steps> 120 </each_steps>
</pair>
</timesteps>
</output>
</time_loop>
<parameters>
<parameter>
<name>lambda</name>
<type>Constant</type>
<value>2.0</value>
</parameter>
<parameter>
<name>c_p</name>
<type>Constant</type>
<value>1500</value>
</parameter>
<parameter>
<name>rho</name>
<type>Constant</type>
<value>1950.0</value>
</parameter>
<parameter>
<name>T0</name>
<type>Constant</type>
<value>10</value>
</parameter>
<parameter>
<name>p_spatial</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>sourceterm_heat_transfer_coefficient</name>
<type>CurveScaled</type>
<curve>point_Neumann</curve>
<parameter>p_spatial</parameter>
</parameter>
<parameter>
<name>ambient_temperature</name>
<type>Constant</type>
<value>10</value>
</parameter>
</parameters>
<curves>
<curve>
<name>point_Neumann</name>
<coords>0 10368000 10368000.01 31104000 31104000.01 41472000 41472000.01 62208000 62208000.01 72576000 72576000.01 93312000</coords>
<values>
-35 -35 0 0 -35 -35 0 0 -35 -35 0 0
</values>
</curve>
</curves>
<process_variables>
<process_variable>
<name>temperature</name>
<components>1</components>
<order>1</order>
<initial_condition>T0</initial_condition>
<boundary_conditions>
<boundary_condition>
<geometrical_set>geometry</geometrical_set>
<geometry>bottom_left</geometry>
<type>Dirichlet</type>
<parameter>ambient_temperature</parameter>
</boundary_condition>
</boundary_conditions>
<source_terms>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po0</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po1</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po2</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po3</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po4</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po5</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po6</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po7</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po8</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po9</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po10</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po11</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po12</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po13</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po14</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po15</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po16</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po17</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po18</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po19</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po20</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po21</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po22</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po23</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
<source_term>
<geometrical_set>geometry</geometrical_set>
<geometry>po24</geometry>
<type>Nodal</type>
<parameter>sourceterm_heat_transfer_coefficient</parameter>
</source_term>
</source_terms>
</process_variable>
</process_variables>
<nonlinear_solvers>
<nonlinear_solver>
<name>basic_picard</name>
<type>Picard</type>
<max_iter>10</max_iter>
<linear_solver>general_linear_solver</linear_solver>
</nonlinear_solver>
</nonlinear_solvers>
<linear_solvers>
<linear_solver>
<name>general_linear_solver</name>
<lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
<eigen>
<solver_type>CG</solver_type>
<precon_type>DIAGONAL</precon_type>
<max_iteration_step>10000</max_iteration_step>
<error_tolerance>1e-16</error_tolerance>
</eigen>
<petsc>
<prefix>gw</prefix>
<parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
</petsc>
</linear_solver>
</linear_solvers>
</OpenGeoSysProject>
\ No newline at end of file
Tests/Data/Parabolic/T/2D_BHE_array/bhe2d.vtu 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
File added
Tests/Data/Parabolic/T/2D_BHE_array/standard_solution_bhe2d_pcs_0_ts_840_t_72576000.000000.vtu 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
File added
scripts/docker/Dockerfile.gcc.full
+ 2
− 1
View file @ d1b043c9
...@@ -57,13 +57,14 @@ RUN apt-get update \ ...@@ -57,13 +57,14 @@ RUN apt-get update \
graphviz \ graphviz \
libxml2-utils \ libxml2-utils \
nodejs \ nodejs \
pandoc \
pandoc-citeproc \ pandoc-citeproc \
yarn yarn
RUN update-alternatives --install /usr/bin/node node /usr/bin/nodejs 10 RUN update-alternatives --install /usr/bin/node node /usr/bin/nodejs 10
# Hugo # Hugo
RUN curl -L -o hugo.tar.gz https://github.com/gohugoio/hugo/releases/download/v0.32.3/hugo_0.32.3_Linux-64bit.tar.gz \ RUN curl -L -o hugo.tar.gz https://github.com/gohugoio/hugo/releases/download/v0.36.1/hugo_0.36.1_Linux-64bit.tar.gz \
&& tar xf hugo.tar.gz \ && tar xf hugo.tar.gz \
&& mv hugo /usr/local/bin/hugo \ && mv hugo /usr/local/bin/hugo \
&& rm -rf hugo.tar.gz LICENSE.md README.md && rm -rf hugo.tar.gz LICENSE.md README.md
......
web/content/docs/benchmarks/heatconduction/BHE_array_benchmark.pandoc 0 → 100644
+ 115
− 0
View file @ d1b043c9
+++
author = "Shuang Chen and Haibing Shao"
date = "2018-02-21T13:44:00+01:00"
title = "BHE Array 2D"
weight = 123
project = "Parabolic/T/2D_BHE_array/bhe2d.prj"
[menu]
[menu.benchmarks]
parent = "heatconduction"
+++
{{< data-link >}}
## Problem description
When shallow geothermal energy is extracted by using Borehole Heat Exchanger (BHE) for heating of buildings, it causes the decrease of subsurface temperature in the vicinity of BHE. In this benchmark, a 2D numerical model has been constructed to model the above temperature variation. The model is validated against the super-imposed analytical solution. Additionally the impact of mesh density to the accuracy of numerical result is also discussed.
## Analytical Solution
For the temperature change in an infinite homogeneous subsurface caused by one single BHE, it can be calculated by the line-source analytical solution (cf. Stauffer et al. (2013) , section 3.1.3), with the assumption that thermal conduction is the only process and no groundwater flow is present. In this case, the ground temperature $T$ is subject to the following equation.
\begin{equation}
T-T_0 = \frac{q_b}{4\pi \lambda}E_1 \frac{r^2}{4\alpha t}
\end{equation}
where $q_b$ is the heat extraction rate on the BHE and $E_1$ denotes the exponential integral function. $T_0$ refers to the initial ground Temperature, and $r$ is the distance between the observation point and the BHE.
In case multiple BHEs are present, Bayer et al. (2014) proposed to calculate the temperature change by super-impose the line-source model of equation (1). The overall temperature change at an observation point with a local coordinate (i, j) can be then calculated as
\begin{equation}
\Delta \mathop T\nolimits_{i,j} \left( {t,\mathop q\nolimits_{k = 1,...,n} } \right) = \sum\limits_{k = 1}^n {\Delta \mathop T\nolimits_{i,j,k} } \left( {t,\mathop q\nolimits_k } \right).
\end{equation}
where ${\mathop q\nolimits_k }$ is a sequence of heat extraction pulses at t =1, ... ,m time steps. In this benchmark, the time step size of heat extraction pulses is set to 120 days, reflecting a 4-month long heat period and the 8-month long recovery interval every year. Within one time step the heat extraction rate on each BHE remains constant. By combining equation (1) and equation (2), the super-imposed temperature change due to the thermal load ${\mathop q\nolimits_k }$ imposed on each BHE can be calculated as
\begin{equation}\begin{split}
\Delta \mathop T\nolimits_{i,j} \left( {t,\mathop q\nolimits_{k,l = 1,...,m} } \right)= \sum\limits_{k = 1}^m {\frac{{\mathop q\nolimits_l - \mathop q\nolimits_{l - 1} }}{{4\pi L\lambda }}} E_1\left[ {\frac{{{{\left( {i - \mathop x\nolimits_k } \right)}^2} + {{\left( {j - \mathop y\nolimits_k } \right)}^2}}}{{4\alpha \left( {\mathop t\nolimits_m - \mathop t\nolimits_l } \right)}}} \right] \\
= \sum\limits_{l = 1}^m {\sum\limits_{k = l}^n {\frac{{\mathop q\nolimits_{k,l} }}{{4\pi L\lambda }}} } \left( {E_1\left[ {\frac{{{{\left( {i - \mathop x\nolimits_k } \right)}^2} + {{\left( {j - \mathop y\nolimits_k } \right)}^2}}}{{4\alpha \left( {\mathop t\nolimits_m - \mathop t\nolimits_{l - 1} } \right)}}} \right] - E_1\left[ {\frac{{{{\left( {i - \mathop x\nolimits_k } \right)}^2} + {{\left( {j - \mathop y\nolimits_k } \right)}^2}}}{{4\alpha \left( {\mathop t\nolimits_m - \mathop t\nolimits_l } \right)}}} \right]} \right).
\end{split}\end{equation}
where ${\mathop q\nolimits_{k,l} }$ is the heat extraction of the k-*th* BHE at time step *l*. The equation (3) will be used to calculate the analytical solution of the overall temperature change in this model for validating the numerical results. It is written in python code and can be found [here](../bhe_array_analytical_solver.py).
## Numerical model setup
A 2D numerical model was constructed and simulated with the Finite element code OpenGeoSys (OGS). The subsurface was represented by a $100 \times 100~m$ square-shaped domain, inside of which 25 BHEs were installed (cf. Figure 1). The distance between adjacent BHEs is kept at 5 m. The ground temperature variation caused by the operation of the BHE array was simulated over a three years’ period. In the model, a 4-month heating period is assumed every year from January to April, with a constant heat extraction rate of 35 W/m on each BHE. In the rest months, the BHE system is shut down and no heat extraction was imposed. The parameters applied in the numerical model can be found in the table below.
In this model, the quad element was adopted to compose the mesh. The initial temperature of the model domain is set to 10 $^{\circ}$C. For the need of modelling a fixed 10 $^{\circ}$C temperature boundary condition was imposed at coordinate (0 m, 0 m). In the model domain, the locations of the BHEs are identified as red dots as in Figure 1. On each of these points, a sink term was specified in the numerical model, with the specific heat extraction rate as listed in the following table.
| Parameter | Symbol | Value | Unit |
| -------------------------------- |:------------ | -------------------:| ----------------:|
| Soil thermal conductivity | $\lambda$ | 1.720 | $Wm^{-1}K^{-1}$ |
| Soil heat capacity | $\rho c$ | $2.925\times10^{6}$ | $J^{-3} mK^{-1}$ |
| Ground thermal diffusivity | $\alpha$ | $5.7\times10^{-7}$ | $Wm^{-1}K^{-1}$ |
| Initial subsurface temperature | $T_0$ | 10 | $^{\circ}C$ |
| Heat extraction rate of the BHE | $q$ | 35 | $W/m$ |
| Diameter of the BHE | $D$ | 0.15 | $m$ |
{{< img src="../BHE_array_benchmark_figures/figure_1.png" >}}
Figure 1: Model geometry, BHE location, and the observation profile
Different meshes were adopted to analyse the impact of mesh density on the numerical results. According to Diersch et al. (2011) the different element size can affect the accuracy of the numerical result significantly for such type of BHE simulation. The optimal element size $\triangle$ in a 2D model around the BHE node should have the following relationship with respect to the BHE diameter:
\begin{equation}
\begin{split}
\Delta = {\rm{ }}a{r_b}\ \hspace{6mm}
a = \left\{ \begin{array}{l}
4.81 \hspace{2mm} for\hspace{2mm} n=4\\
6.13 \hspace{2mm} for\hspace{2mm} n=6\\
6.16 \hspace{2mm} for\hspace{2mm} n=8
\end{array}\right.
\label{eq_4}
\end{split}
\end{equation}
where $r_b$ is the BHE radius. n denotes the number of surrounding nodes. n = 8 is typical for a squared grid meshes. In this study, the BHE diameter is assumed to be 0.15 m. Based on equation (4) the optimal element size should be set to approximately 0.5 m.
## Numerical modelling results
Figure 2 and 3 show the comparison of the temperature distribution along the observation profile (position see Figure 1) using analytical solution with the numerical results from OGS5 and OGS6 for every 4 months in the whole simulated time. It shows the numerical solution has a very good agreement with the analytical solution.
{{< img src="../BHE_array_benchmark_figures/figure_2.png" width="200">}}
Figure 2: The temperature evolution of the BHEs field along the observation profile
{{< img src="../BHE_array_benchmark_figures/figure_3.png" width="200">}}
Figure 3: The temperature evolution of the BHEs field along the observation profile
In order to investigate the impact of mesh density on the accuracy of numerical result, the simulated temperature profile at the observation point A (53 m, 52.5 m) was plotted and compared against the analytical solution. Figure 3 shows the relative difference of the computed temperature between the analytical and numerical solution by using different mesh size (2.5 m, 1 m, 0.5 m, 0.25 m and 0.2 m). The results show that the difference becomes smaller when the mesh size is approaching 0.5 m, which is expected as the optimal mesh size mentioned in Diersch et al. (2011). From Figure 4, it can be found that the absolute error of temperature values at point A should be less than 2.5e-3 if the mesh size is kept denser than 0.5m.
{{< img src="../BHE_array_benchmark_figures/figure_4.png" width="200">}}
Figure 4: The relative difference of computed temperature at point A between the analytical and numerical solution using different mesh size
{{< img src="../BHE_array_benchmark_figures/figure_5.png" width="200">}}
Figure 5: The absolute difference of computed temperature along the diagonal profile between the analytical and numerical solution using different mesh size
## Summary
In this benchmark, a 2D numerical model has been constructed to simulate the ground temperature variation caused by heat extraction by BHE. The results show a very good agreement against the analytical solution. Additionally the impact of mesh density on the accuracy of numerical result was also investigated, and the optimal mesh size was found to be 0.5 m. In future studies, the pipe-line network system will be implemented into the model and coupled with the BHE.
## References
[1] Bayer, P., de Paly, M., Beck, M., 2014. Strategic optimization of borehole heat exchanger field for seasonal geothermal heating and cooling. Applied Energy 136, 445-453.
URL http://dx.doi.org/10.1016/j.apenergy.2014.09.029
[2] Diersch, H.-J., Bauer, D., Heidemann, W., Ruehaak, W., Sch¨atzl, P., 2011. Finite element modeling of borehole heat exchanger systems: Part 2. Numerical simulation. Computers & Geosciences 37, 1136-1147.
[3] Stauffer, F., Bayer, P., Blum, P., Giraldo, N., Kinzelbach, W., 2013. Thermal Use of Shallow Groundwater. CRC Press, 290.
web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_1.png 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_1.png

131 B

web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_2.png 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_2.png

131 B

web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_3.png 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_3.png

131 B

web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_4.png 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_4.png

131 B

web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_5.png 0 → 100644 LFS
+ 3
− 0
View file @ d1b043c9
web/content/docs/benchmarks/heatconduction/BHE_array_benchmark_figures/figure_5.png

131 B

web/content/docs/benchmarks/heatconduction/bhe_array_analytical_solver.py 0 → 100644
+ 332
− 0
View file @ d1b043c9
# -*- coding: utf-8 -*- {}
import matplotlib.pyplot as plt
import numpy as np
import scipy
import os
import math
from mpmath import *
mp.dps = 25; mp.pretty = True
#part 1:gif_splitfile to txt
#新建文件夹txt，用于存储原文件的分解文件集。
os.mkdir("C://george//PhD//papers//BHE_sc//result//txt")
#os.chdir(path) 方法用于改变当前工作目录到指定的路径。此处创建路径到新建文件夹txt
os.chdir("C://george//PhD//papers//BHE_sc//result//txt")
txtfiles = os.getcwd() #获得当前运行脚本所在目录作为文件储存地址
# define the path of the data file
str_tec_file_path = "C://george//PhD//papers//BHE_sc//result//tec_0_5m.tec"
data_num =0
with open(str_tec_file_path, 'r') as f_in:
i = -1 #此处写-1是因为要使输出文件名第一个必须是以0结尾，不然和part2的文件读入循环对不上。必须要对上是因为\n
#part2数组值从0位开始储存的。
for line in f_in:
if '=' not in line:
txtfiles.write(line)
#读出单位文件中有多少行数据，存入变量data_num
data_num = data_num + 1
elif 'ZONE' in line: #从关键词ZONE开始读取行
data_num = 0
i = i + 1
txtfiles = open('txtfile{}.txt'.format(i), 'w')#格式化命名文件名存储
#记录总共的分出文件数目，用于part 2
numfile = i + 1
#关闭并完整文件释放内存
txtfiles.close()
# plotting T_in T_out curves,建立作图文件夹及及改变工作路径
os.mkdir("C://george//PhD//papers//BHE_sc//result//png")
os.chdir("C://george//PhD//papers//BHE_sc//result//png")
pngfiles = os.getcwd()
#建立source term坐标矩阵，用以计算各热源对referece point的温度影响矩阵
po_x = np.array([40,40,40,40,40],dtype=float).reshape(-1,1) + np.arange(0,25,5)
po_y = np.arange(60,35,-5).reshape(-1,1) + np.zeros(5)
po_dist_to_referencepo = np.zeros([5,5])
Temp_po_to_referencepo = np.zeros([5,5])
#解析解1基础数据项,q是timestep变化量
q1 = -35
q2 = 35
q3 = 0
time_trans = 120*24*60*60 #3个月一输出数据
T0 = 10
T = T0
poro = 10e-20
lamda_f = 0.5984
density_f = 998.2032
cp_f = 4182.0
lamda_sp = 2.0
density_sp = 1950.0
cp_sp = 1500.0
alpha_f = lamda_f/(density_f*cp_f)
alpha_sp = lamda_sp/(density_sp*cp_sp)
lamda_med = (1 - poro)*lamda_sp + poro*lamda_f
alpha_med = (1 - poro)*alpha_sp + poro*alpha_f
#解析解2基础数据项：
qq = np.array([-35,0,0,-35,0,0,-35,0,0,-35,0]).reshape(1,-1)
qq_all = np.repeat(qq,data_num,axis=0)
expandedloads= np.array(qq)
numtimesteps = 11
numbhe = 25
ppo_x = np.array([40,40,40,40,40],dtype=float).reshape(-1,1) + np.arange(0,25,5)
ppo_y = np.arange(60,35,-5).reshape(-1,1) + np.zeros(5)
ppo_x_re = np.reshape(po_x, (-1,1))
ppo_y_re = np.reshape(po_y, (-1,1))
# 建立动态变量名
createVar = locals()
# 建立月份名字表，因为数据记录从1月1日有一组数据，所以多出一个月数据。
month = [" Temperature after 0 months ", " Temperature after 4 months", " Temperature after 8 months", " Temperature after 12 months", " Temperature after 16 months", " Temperature after 20 months",
" Temperature after 24 months", " Temperature after 28 months", " Temperature after 32 months", " Temperature after 36 months"," Temperature after 40 months"," Temperature after 44 months"]
#解析解2项：
txtpath2 = f"C://george//PhD//papers//BHE_sc//result//txt//txtfile0.txt"#f格式化字符串
dist_arrayy = 0.0
with open(txtpath2, 'r') as f_in_txt:
for line in f_in_txt:
# split line
data_line = line.split(" ") #maxsplit
# time_double = 0.0
dist_double2 = float(data_line[0])
cordinate2 = math.sqrt(dist_double2**2/2)
# put into the list
dist_arrayy= np.vstack((dist_arrayy, cordinate2))
dist_arrayy_new = np.delete(dist_arrayy, 0, 0)
numtemppoints = len(dist_arrayy_new)
T2=np.zeros([numtemppoints,numtimesteps])
coeff_all = np.zeros([numtemppoints,numtimesteps])
for currstep in range(0,numtimesteps):
Temp_po_to_referencepo= np.zeros([numtemppoints,numbhe])
po_dist_to_referencepo= np.zeros([numtemppoints,numbhe])
localcoeff_all= np.zeros([numtemppoints,1])
localcoeff= np.zeros([numtemppoints,numbhe])
localcoeff1= np.zeros([numtemppoints,numbhe])
for i in range(0,numbhe):
if(time_trans*(currstep+1)-time_trans*0>0):
for j in range(0,numtemppoints):
po_dist_to_referencepo[j,i] = abs(ppo_x_re[i] - (dist_arrayy_new[j]+1))**2 + abs(ppo_y_re[i] - dist_arrayy_new[j])**2
exp = po_dist_to_referencepo[j,i]/(4*alpha_med*time_trans*(currstep+1))
n = e1(exp)
localcoeff[j,i] = 1/(4*math.pi*lamda_med)*n
if(time_trans*(currstep+1)-time_trans*1>0):
for j in range(0,numtemppoints):
po_dist_to_referencepo[j,i] = abs(ppo_x_re[i] - (dist_arrayy_new[j]+1))**2 + abs(ppo_y_re[i] - dist_arrayy_new[j])**2
exp1 = po_dist_to_referencepo[j,i]/(4*alpha_med*time_trans*currstep)
n1 = e1(exp1)
localcoeff[j,i] = localcoeff[j,i] - 1/(4*math.pi*lamda_med)*n1
localcoeff_all= np.sum(localcoeff,axis=1).reshape(-1,1)
coeff_all[:,1:]=coeff_all[:,:numtimesteps-1]
coeff_all[:,:1]=localcoeff_all
for currstep in range(0,numtimesteps):
T2[:,currstep] = np.sum(coeff_all[:,numtimesteps-1-currstep:]*qq_all[:,:currstep+1],axis=1) +10
# loop over
for m in range(0 , numfile):
#文件读取路径
txtpath = f"C://george//PhD//papers//BHE_sc//result//txt//txtfile{m}.txt"#f格式化字符串
txtpath_ogs6 = f"C://george//PhD//papers//BHE_sc//result//txt_ogs6//txtfile{m}.txt"
#动态变量名每组赋初值0.0
createVar['dist_array'+ str(m)] = 0.0
createVar['temperature_array_ogs5'+ str(m)] = 0.0
createVar['temperature_array_ana'+ str(m)] = 0.0
createVar['dist_array2'+ str(m)] = 0.0
createVar['temperature_array_ana2'+ str(m)] = 0.0
createVar['dist_array_ogs6'+ str(m)] = 0.0
createVar['temperature_array_ogs6'+ str(m)] = 0.0
#解析解1：
# if m == 0 or m ==2 or m ==3 or m ==5 or m ==6 or m ==8 or m ==9 or m ==10:
# q = q2
# elif m == 1 or m ==4 or m ==7:
# q = q1
if m == 0:
with open(txtpath, 'r') as f_in_txt:
for line in f_in_txt:
# split line
data_line = line.split(" ") #maxsplit
# time_double = 0.0
dist_double = float(data_line[0])
cordinate = math.sqrt(dist_double**2/2)
temperature = float(data_line[1])
#解析解数据项
time = time_trans + 10e-6
for i in range(0,5):
for j in range(0,5):
po_dist_to_referencepo[i,j] = abs(po_x[i,j]-(cordinate+1))**2 + abs( po_y[i,j]-cordinate)**2
for i in range(0,5):
for j in range(0,5):
exp = po_dist_to_referencepo[i,j]/(4*alpha_med*time)
n1 = e1(exp)
Temp_po_to_referencepo[i,j] = q3/(4*math.pi*lamda_med)*n1
T = np.sum(Temp_po_to_referencepo) + T0
# put into the list
createVar['dist_array'+ str(m)] = np.vstack((createVar['dist_array'+ str(m)], cordinate))
createVar['temperature_array_ogs5'+ str(m)] = np.vstack((createVar['temperature_array_ogs5'+ str(m)], temperature))
createVar['temperature_array_ana'+ str(m)] = np.vstack((createVar['temperature_array_ana'+ str(m)], T))
f_in_txt.close()
# end of loop
#去除21和22行的0.0值，因为此值被带入了dist_array的list第一行，可查看variable explorer temperature_array观察。
#目的是去除做图时的第一个数据0值。
dist_array_new = np.delete(createVar['dist_array'+ str(m)], 0, 0)
temperature_array_ogs5_new = np.delete(createVar['temperature_array_ogs5'+ str(m)], 0, 0)
temperature_array_ana_new = np.delete(createVar['temperature_array_ana'+ str(m)], 0, 0)
with open(txtpath_ogs6, 'r') as f_in_txt:
for line in f_in_txt:
# split line
data_line = line.split(" ") #maxsplit
# time_double = 0.0
dist_double = float(data_line[0])
cordinate = math.sqrt(dist_double**2/2)
temperature = float(data_line[1])
createVar['dist_array_ogs6'+ str(m)] = np.vstack((createVar['dist_array_ogs6'+ str(m)], cordinate))
createVar['temperature_array_ogs6'+ str(m)] = np.vstack((createVar['temperature_array_ogs6'+ str(m)], temperature))
f_in_txt.close()
dist_array_new_ogs6 = np.delete(createVar['dist_array_ogs6'+ str(m)], 0, 0)
temperature_array_ogs6_new = np.delete(createVar['temperature_array_ogs6'+ str(m)], 0, 0)
#plotting
plt.figure()
plt.plot(scipy.dot(dist_array_new,1.0),scipy.dot(temperature_array_ogs5_new,1.0),color = 'r',ls=':',lw=1, marker='^',markersize=1.5, label= 'OGS5')
plt.plot(scipy.dot(dist_array_new_ogs6,1.0),scipy.dot(temperature_array_ogs6_new,1.0),c='g', ls=':',lw=1, marker='o',markersize=1, label= 'OGS6')
# plt.plot(scipy.dot(dist_array_new,1.0),scipy.dot(temperature_array_ana_new,1.0), color = 'g',label= 'ana1' + months[m])
plt.plot(scipy.dot(dist_array_new,1.0),scipy.dot(temperature_array_ana_new,1.0), "b",label= 'Analytical')
plt.xlim([0,100])
plt.ylim([-10,20])
plt.ylabel('Temperature [$^\circ$C]')
plt.xlabel('x [m]')
plt.legend(loc='best',fontsize=8)
plt.title(month[m],fontsize=12)
# plt.show()
plt.savefig('pngfile{}.png'.format(m),dpi = 300, transparent = False)
else:
#解析解数据项,实际是温度场斜三角矩阵相加
T = np.zeros([12])
with open(txtpath, 'r') as f_in_txt:
for line in f_in_txt:
# split line
data_line = line.split(" ") #maxsplit
# time_double = 0.0
dist_double = float(data_line[0])
cordinate = math.sqrt(dist_double**2/2)
i1 = m +1
for m1 in range(1, m + 1 ):
i1 = i1 - 1
#时间变量，timestep温度随时间叠加矩阵每行重新设为0
time = time_trans * i1 + 10e-6
if m1 == 2 or m1 == 5 or m1 == 8 or m1 == 11:
q = q2
elif m1 == 1 or m1 ==4 or m1 ==7 or m1 ==10:
q = q1
elif m1 == 3 or m1 == 6 or m1 == 9:
q = q3
for i in range(0,5):
for j in range(0,5):
po_dist_to_referencepo[i,j] = abs(po_x[i,j]-(cordinate+1))**2 + abs( po_y[i,j]-cordinate)**2
for i in range(0,5):
for j in range(0,5):
exp = po_dist_to_referencepo[i,j]/(4*alpha_med*time)
n1 = e1(exp)
Temp_po_to_referencepo[i,j] = q/(4*math.pi*lamda_med)*n1
T[m1] = np.sum(Temp_po_to_referencepo)
T_sum = np.sum(T) + T0
createVar['temperature_array_ana'+ str(m)] = np.vstack((createVar['temperature_array_ana'+ str(m)], T_sum))
f_in_txt.close()
#ogs5项：
with open(txtpath, 'r') as f_in_txt:
for line in f_in_txt:
# split line
data_line = line.split(" ") #maxsplit
# time_double = 0.0
dist_double = float(data_line[0])
cordinate = math.sqrt(dist_double**2/2)
temperature = float(data_line[1])
# put into the list
createVar['dist_array'+ str(m)] = np.vstack((createVar['dist_array'+ str(m)], cordinate))
createVar['temperature_array_ogs5'+ str(m)] = np.vstack((createVar['temperature_array_ogs5'+ str(m)], temperature))
f_in_txt.close()
#去除21和22行的0.0值，因为此值被带入了dist_array的list第一行，可查看variable explorer temperature_array观察。
#目的是去除做图时的第一个数据0值。
dist_array_new = np.delete(createVar['dist_array'+ str(m)], 0, 0)
temperature_array_ogs5_new = np.delete(createVar['temperature_array_ogs5'+ str(m)], 0, 0)
temperature_array_ana_new = np.delete(createVar['temperature_array_ana'+ str(m)], 0, 0)
#ogs6项
with open(txtpath_ogs6, 'r') as f_in_txt:
for line in f_in_txt:
# split line
data_line = line.split(" ") #maxsplit
# time_double = 0.0
dist_double = float(data_line[0])
cordinate = math.sqrt(dist_double**2/2)
temperature = float(data_line[1])
# put into the list
createVar['dist_array_ogs6'+ str(m)] = np.vstack((createVar['dist_array_ogs6'+ str(m)], cordinate))
createVar['temperature_array_ogs6'+ str(m)] = np.vstack((createVar['temperature_array_ogs6'+ str(m)], temperature))
f_in_txt.close()
dist_array_new_ogs6 = np.delete(createVar['dist_array_ogs6'+ str(m)], 0, 0)
temperature_array_ogs6_new = np.delete(createVar['temperature_array_ogs6'+ str(m)], 0, 0)
#plotting
plt.figure()
plt.plot(scipy.dot(dist_array_new,1.0),scipy.dot(temperature_array_ogs5_new,1.0),color = 'r',ls=':',lw=1, marker='^',markersize=1.5, label= 'OGS5')
plt.plot(scipy.dot(dist_array_new_ogs6,1.0),scipy.dot(temperature_array_ogs6_new,1.0),c='g', ls=':',lw=1, marker='o',markersize=1, label= 'OGS6')
# plt.plot(scipy.dot(dist_array_new,1.0),scipy.dot(temperature_array_ana_new,1.0), color = 'g',label= 'ana' + month[m])
plt.plot(scipy.dot(dist_array_new,1.0),T2[:,m-1],"b",label= 'Analytical')
plt.xlim([0,100])
plt.ylim([-10,20])
plt.ylabel('Temperature [$^\circ$C]')
plt.xlabel('x [m]')
plt.legend(loc='best',fontsize=8)
plt.title(month[m],fontsize=12)
# plt.show()
plt.savefig('pngfile{}.png'.format(m), dpi = 300, transparent = False)
# end of loop
txtpath.close()
pngfiles.close()
\ No newline at end of file
web/content/docs/benchmarks/heatconduction/bhe_array_benchmark.bib 0 → 100644
+ 41
− 0
View file @ d1b043c9
@article{Bayer,
author = {Bayer, Peter and de Paly, Michael and Beck, Markus},
doi = {10.1016/j.apenergy.2014.09.029},
file = {:C$\backslash$:/george/PhD/Lieteratur/Bayer{\_}AE2014.pdf:pdf},
issn = {03062619},
journal = {Applied Energy},
keywords = {Closed systems,Finite line source,Low-enthalpy geothermics,Seasonal use,Simulation},
mendeley-groups = {PhD/artikel/BHE{\_}benchmark},
pages = {445--453},
publisher = {Elsevier Ltd},
title = {{Strategic optimization of borehole heat exchanger field for seasonal geothermal heating and cooling}},
url = {http://dx.doi.org/10.1016/j.apenergy.2014.09.029},
volume = {136},
year = {2014}
}
@article{Diersch,
author = {Diersch, H.-J.G. and Bauer, D and Heidemann, W and R{\"{u}}haak, W and Sch{\"{a}}tzl, P},
doi = {10.1016/j.cageo.2010.08.002},
file = {:D$\backslash$:/Dropbox/Literature/geothermal{\_}heat{\_}pump/Diersch{\_}FEFLOW{\_}2{\_}2011.pdf:pdf},
issn = {0098-3004},
journal = {Computers {\&} Geosciences},
keywords = {Borehole heat exchanger},
mendeley-groups = {Automatically Imported},
pages = {1136--1147},
title = {{Finite element modeling of borehole heat exchanger systems: Part 2. Numerical simulation}},
volume = {37},
year = {2011}
}
@article{Stauffer,
author = {Stauffer, Fritz and Bayer, Peter and Blum, Philipp and Giraldo, Nelson and Kinzelbach, Wolfgang},
doi = {10.1201/b16239},
file = {:C$\backslash$:/george/PhD/Lieteratur/K15881{\_}CoverDraft{\_}01.pdf.pdf:pdf},
isbn = {978-1-4665-6019-2},
journal = {CRC Press},
pages = {290},
title = {{Thermal Use of Shallow Groundwater}},
year = {2013}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment