The total length of the tubes per square meter living area is \(L=\dfrac{1.5\times16}{15}=1.6\ (m^{-1})\). This is a much higher length per area than deep borehole geothermal systems, which also means smaller power per length.
| Month | Power.per.Length.per.Area |
|---|---|
| Units | W/m^3 |
| 1 | 0.145833333333333 |
| 2 | 0.5625 |
| 3 | 2.125 |
| 4 | 4.10416666666667 |
| 5 | 6.14583333333333 |
| 6 | 7.39583333333333 |
| 7 | 6.72916666666667 |
| 8 | 5.83333333333333 |
| 9 | 4.95833333333333 |
\(\lambda\) for PE is \(0.48 \dfrac{W}{mK}\). The thermal resistance of a radially symmetric is \[R_{th}=\int^{r_{out}}_{r_{in}}\dfrac{\nabla T}{q(r)}dr=\int^{r_{out}}_{r_{in}}\dfrac{ \nabla T}{\lambda A(r) \nabla T} dr=\int^{r_{out}}_{r_{in}}\dfrac{dr}{2\pi \lambda r}\\R_{th}=\dfrac{1}{2\pi\lambda}ln(\dfrac{r_{out}}{r_{in}})\]
So the thermal resistance is \(R_{th}=\dfrac{1}{2\pi\lambda}ln(\dfrac{16}{16-3})=0.06885 (\dfrac{m}{WK})\). The apparent radius is \(r_a=r_be^{-2\pi\lambda R_{th}}=5.4256 mm\).
Note that we obtain a lower temperature than the answer. If we assume \(\lambda_{PE}\) to equal to \(\lambda_{rock}\), we will obtain the same temperature as the answer, as shown below. This is also true in 8(d).
In the assignment it was not clarified temperature from which tube should we calculate. So we solve the general space-time convolution problem here.
L=16/15
lambda=2.5
rho=2300
c=900
alpha=lambda/rho/c
lambda_pe=0.48
r_b=0.016
R_th=1/lambda_pe/2/pi*log(r_b/(r_b-0.003))
r_a=r_b*exp(-2*pi*lambda*R_th)
eta_rel=.5
T_H=35
time=seq(0,9)
P_Demand=c(0,.7,2.7,10.2,19.7,29.5,35.5,32.3,28.0,23.8)*(1/3)/L
T_m=10
x=seq(0,15,.1)
y=seq(0,3,.1)
T_field=array(dim=c(length(x),length(y),length(time)))
T_wall=matrix(nrow=16,ncol=length(time))
T_wall[,1]=T_m
P_geo=matrix(nrow=16,ncol=length(time))
diff_P_geo=matrix(nrow=16,ncol=length(time))
dis=array(dim=c(length(x),length(y),16))
for(i in 1:length(x)){
for(j in 1:length(y)){
for(k in 1:16){
dis[i,j,k]=(x[i]-k+1)^2+(y[j]-1.5)^2
if(dis[i,j,k]==0){
dis[i,j,k]=r_a^2
}
}
}
}
dis_vir=array(dim=c(length(x),length(y),16))
for(i in 1:length(x)){
for(j in 1:length(y)){
for(k in 1:16){
dis_vir[i,j,k]=(x[i]-k+1)^2+(y[j]+1.5)^2
}
}
}
sum_0=matrix(0,nrow=length(x),ncol=length(y))
sum=sum_0
u=array(dim=c(length(x),length(y),16))
u_vir=array(dim=c(length(x),length(y),16))
for(i in 1:length(time)){
if(i==1){
T_wall[,i]=T_m
T_field[,,i]=T_m
P_geo[,i]=P_Demand[i]/(1+(T_H-T_wall[,i])/(273.15+T_wall[,i])/eta_rel)
diff_P_geo[,i]=0
}
else{
sum=sum_0
for(j in 1:(i-1)){
u=dis/4/alpha/((i-j)*30*86400)
u_vir=dis_vir/4/alpha/((i-j)*30*86400)
for(k in 1:16){
sum=sum+diff_P_geo[k,j]*(expint_E1(u[,,k])-expint_E1(u_vir[,,k]))/4/pi/lambda
}
}
T_wall[,i]=T_m-sum[(c(1:16)-1)/.1+1,1.5/0.1+1]
T_field[,,i]=T_m-sum
P_geo[,i]=P_Demand[i]/(1+(T_H-T_wall[,i])/(273.15+T_wall[,i])/eta_rel)
diff_P_geo[,i]=P_geo[,i]-P_geo[,i-1]
}
}Below we plot the temperature field after nine months. The mean wall temperature is \(5.22\) degree Celsius.
As shown in the previous section, the minimum Temperature occurs at the center (the 8th and the 9th tubes). This implies that we calculate only the temperature evolution there.
Since the heat equation is linear, we can decompose it to a harmonic part (obtained in assignment four) and a convolution part which is calculated in the previous solution. A small catch with this method is that the extractable power must be adjusted accordingly to the changing temperature.
It is easy to see that in winter, there exist risks that the water would freeze.
With an algorithm similar to that in 5(h), we obtain an annual additional electricity need of \(8.93 kWh\). This is \(8.8\%\) less than the longer borehole length case in 5(h).
With a denser configuration, the power per length decreases and thus the temperature drop is lower at the wall. This reduces the additional required electricity than the original system although some of this effect is offset by the fact that the borehole in the center is now affected by more other boreholes.
From the first law of thermodynamics we have \(Q=Q_{in}-Q_{out}=\rho_wQc_w(T_H-T_L)\). From the second law we have \(\eta_{Carnot}=1-\dfrac{T_L}{T_H}\) and \(W=\eta Q\). So there is a linear relationship between flow rate and output. At \(300 \dfrac{L}{s}\), we get \(16.66 MW\) for maximum work and \(94.02 MW\) for thermal power.
The pressure at one of the doublet well is \[p=\dfrac{\mu Q}{2\pi kl}ln(\dfrac{r_e}{r_i})\] and the difference is two times of this. Meanwhile, the required work from a pump to make the water circulate is \[P=\Delta p Q=\dfrac{\mu Q^2}{\pi kl}ln(\dfrac{r_e}{r_i})\] To compute these values we assumed constant viscosity, which is of course not the case. Like density \(\rho\), viscosity \(\mu\) is heavily sensitive to temperature change.
Below we see that pressure has a linear relationship with flow rate and pump work has quadratic relationship with flow rate.
As shown below, the net power is an quadratic function to flow rate. Optimizing yields \(Q^*=175 \dfrac{L}{s}\).
The diffusion velocity is \(v=-\dfrac{k}{\mu}\nabla p\), on the straight line we will have \[\nabla p=\dfrac{\partial p}{\partial r}=-\dfrac{\mu Q}{2\pi kl}(\dfrac{1}{r_e}+\dfrac{1}{r_i})\\v=\dfrac{Q}{2\pi l}(\dfrac{1}{r_e}+\dfrac{1}{r_i})\\v_a=\dfrac{\rho_f c_f}{\rho_m c_m+\phi \rho_f c_f} v\]
Where we will convenient take the approximation \(v_a \sim 1.5 v\) here.
Below we see that the advection velocity is maximized at the wells (about \(1.1606\times10^{-3}\dfrac{m}{s}\) or \(36625 \dfrac{m}{yr}\)) and minimized at the center (about \(3.7136\times10^{-7}\dfrac{m}{s}\) or \(11.7190 \dfrac{m}{yr}\)).
\[T_{residue}=\int^{d}_{0}\dfrac{dr}{v_a(r)}=\dfrac{2\pi l}{Q}\int^{d}_{0}\dfrac{r(d-r)dr}{d}\\=\dfrac{2\pi l}{Qd}\int^{d}_{0}r(d-r)dr=\dfrac{2\pi l}{Qd}(\dfrac{r^2d}{2}-\dfrac{r^3}{3})|^d_0=\dfrac{\pi l d^2}{3Q}\]
Which is about \(4.0392\times10^9\) seconds, or \(127.997\) years, a pretty safe value for a project with lifespan of 30 to 40 years.
Note that this is the minimum residual time of the water. Choosing other paths yields a longer residual time. The expected residual time should be computed with conformal mapping techniques (setting \(w=\dfrac{1}{z-\dfrac{d}{2}}-\dfrac{1}{z+\dfrac{d}{2}}\) will transform the function onto a complex plane).
If we have a shorter distance for the doublet, we have the same maximum power output but less required pumpwork due to less pressure drop. This will result in higher net power output at a new optimum flow rate \(Q^*=183 \dfrac{L}{s}\).
The advection velocity at the wells will increase to around \(v_a=1.2137\times10^{-3}\dfrac{m}{s}\) or \(38301 \dfrac{m}{yr}\), and the advection velocity at the center is increased to \(v_a=5.8251\times10^{-7}\dfrac{m}{s}\) or \(18.3821 \dfrac{m}{yr}\).
The residual time along the straight line will be reduced to \(1.7167\times10^9\) seconds, or \(54.401\) years. This is the most significant change of the results, but still a safe value when considering the project.