Assume on a 2-D plane, we have a constant power source \(P_0\) at the origin, and we would like to find a solution of the temperature field with polar cooridnates of \(r>0\).
The Laplacian in a 2-D polar coordinate is \[\nabla ^{2}=\dfrac{\partial}{r \partial r}r\dfrac{\partial}{\partial r}+\dfrac{\partial^2}{r^2 \partial \theta^2}\]
Since the solution is obviously radially symmetric, the PDE can be simplified into \[\dfrac{\partial T}{\partial t}= \dfrac{\lambda}{\rho c}\nabla ^{2}T=\alpha\nabla ^{2}T=\alpha \dfrac{\partial}{r \partial r}r\dfrac{\partial T}{\partial r}\]
We seek a variable \(u=u(r,t)\) that can describe the similarity of the solution. The rationale behind this is that with two variables and one equation, we should be able to find a characteristic variable that can describe the system sufficiently.
Apply the chain rule to both sides we shall have \[\dfrac{\partial u}{\partial t}\dfrac{\partial T}{\partial u}=\dfrac{\alpha}{r} \dfrac{\partial u}{\partial r}\dfrac{\partial}{\partial u}r\dfrac{\partial u}{\partial r}\dfrac{\partial T}{\partial u}\]
Rearranging yields \[\dfrac{r}{\alpha} \dfrac{\dfrac{\partial u}{\partial t}}{\dfrac{\partial u}{\partial r}}\dfrac{\partial T}{\partial u}=-\dfrac{r}{\alpha} \dfrac{\partial r}{\partial t}|_u\dfrac{\partial T}{\partial u}=\dfrac{\partial}{\partial u}r\dfrac{\partial u}{\partial r}\dfrac{\partial T}{\partial u}\]
Because this differential equation should be an ODE of variable u, we would thus have \[\begin{cases} r\dfrac{\partial u}{\partial r}=f(u)\\ r\dfrac{\partial r}{\partial t}|_u=g(u)\end{cases}\]
Note that the selection of u is not unique (any variable other than the chosen u that satisfies \(u'=h(u)\) is also a qualified candidate for the characteristic method). So suppose we choose the form \(u=kr^a t^b\), then we have \[\begin{cases} r\dfrac{\partial u}{\partial r}=a u\\ r\dfrac{\partial r}{\partial t}|_u=-r \dfrac{bu/t}{au/r}=-\dfrac{br^2}{at}\end{cases}\]
Observe that if we have \(a=2\) and \(b=-1\), the equations would become \[\begin{cases} r\dfrac{\partial u}{\partial r}=2 u \\ r\dfrac{\partial r}{\partial t}|_u=\dfrac{u}{2} \end{cases}\]
Then our ODE becomes \[-\dfrac{1}{4\alpha}(u\dfrac{\partial T}{\partial u})=\dfrac{\partial}{\partial u}(u\dfrac{\partial T}{\partial u})\]
Which yields \[(u\dfrac{\partial T}{\partial u})=(\lim _{u\rightarrow 0}(u\dfrac{\partial T}{\partial u}))e^{-\dfrac{u}{4\alpha}}\]
The limit is necessary since \(u=0\) is not in the domain of definition. However, we know that the total power delivered outwards at a given \(r\) should approach \(P_0\) as \(r \rightarrow 0\). This yields \[P_0=-\lim _{r\rightarrow 0}2\pi r q(r)=lim _{r\rightarrow 0}2\pi r \lambda\dfrac{\partial T}{\partial r}=lim _{r\rightarrow 0}4\pi\lambda u \dfrac{\partial T}{\partial u}\\lim _{r\rightarrow 0}u \dfrac{\partial T}{\partial u}=\dfrac{P_0}{4\pi\lambda}\]
Thus we have \[(u\dfrac{\partial T}{\partial u})=\dfrac{P_0}{4\pi\lambda}e^{-\dfrac{u}{4\alpha}}\\T(u)=C+\dfrac{P_0}{4\pi\lambda}\int \dfrac{e^{-\dfrac{u'}{4\alpha}}}{u'}du'\] Rescaling by a factor of \(\dfrac{1}{4\alpha}\) we have \[T(u)=C+\dfrac{P_0}{4\pi\lambda}\int \dfrac{e^{-u'}}{u'}du'\]
For boundary condition \(T(u) \rightarrow 0\) as \(u \rightarrow \infty\) it is easy to see that \[T(r,t)=T(u)=\dfrac{P_0}{4\pi\lambda}\int_u^\infty \dfrac{e^{-u'}}{u'}du'=\dfrac{P_0}{4\pi\lambda}E_1(u)=\dfrac{P_0}{4\pi\lambda}E_1(\dfrac{r^2}{4\alpha t})\]
Since the function \(\dfrac{e^{-\dfrac{u}{4\alpha}}}{u}\) approaches zero at spectrum convergence rate as u approaches infinity.
It is important to realize that the heat sources form a linear space-time invariant (LSTI) system. That is, for a heat source field \(P(\overrightarrow {r},t)\) and initial temperature field \(T(\overrightarrow {r},0)=0\) everywhere, we can derive the temperature field as the following convolution formula: \[T(\overrightarrow {r},t)=\dfrac{1}{4\pi\lambda}\int_{Domain} \dfrac{\partial P(\overrightarrow {r'},\tau)}{\partial \tau}E_1(\dfrac{||\overrightarrow {r}-\overrightarrow {r'}||^2}{4\alpha (t-\tau)}) d\overrightarrow {r'}d\tau\]
As a comparison, it is also possible to derive a LSTI form for given variations of temperature as boundary conditions: \[T(\overrightarrow {r},t)=\int_{Domain} \dfrac{\partial T_{BC}(\overrightarrow {r'},\tau)}{\partial \tau}\dfrac{exp(-\dfrac{||\overrightarrow {r}-\overrightarrow {r'}||^2}{4\alpha (t-\tau)})}{\sqrt{4\pi\alpha(t-\tau)}} d\overrightarrow {r'}d\tau\] Which generalizes the solution and thus can be applied in both assignment three and assignment four at homework 1.
Note that in both cases, if the boundary conditions are constants, the derivative of \(P\) or \(T_{BC}\) would be a delta function with respect to time.
We now formalize the heat equation as the following: \[\dfrac{\partial T}{\partial t}= \dfrac{\lambda}{\rho c}\nabla ^{2}T+\dfrac{P_0}{\rho c}\delta(\overrightarrow r)=\alpha\nabla ^{2}T+\dfrac{P_0}{\rho c}\delta(\overrightarrow r)\]
Taking the Fourier transform of this equation with respect to the spatial variables one shall have \[\dfrac{\partial F[T]}{\partial t}=-4\pi^2\alpha s^2 F[T]+\dfrac{\alpha P_0}{\lambda}\]
Solving it yields \[F[T]=F[T](\overrightarrow s,0)e^{-4\pi^2\alpha s^2 t}+\dfrac{P_0}{4\pi^2\lambda s^2}(1-e^{-4\pi^2\alpha s^2 t})\]
The first term on the right hand side is what we have derived in assignment three of homework one. With our initial condition we can conveniently neglect it.
Let us now rearrange the equations into \(-4\pi^2s^2F[T]=\dfrac{P_0}{\lambda}(e^{-4\pi^2\alpha s^2 t}-1)\) and take the inverse Fourier Transform of it, noticing that the function is radially symmetric: \[\dfrac{\partial}{r \partial r}r\dfrac{\partial T}{\partial r}=\dfrac{P_0}{\lambda}(\dfrac{e^{-\dfrac{r^2}{4\alpha t}}}{4\pi\alpha t}-\delta(r))\\ \dfrac{\partial}{\partial r}r\dfrac{\partial T}{\partial r}=\dfrac{P_0}{\lambda}(r\dfrac{e^{-\dfrac{r^2}{4\alpha t}}}{4\pi\alpha t}-r\delta(r))\\ r\dfrac{\partial T}{\partial r}=\dfrac{P_0}{2\pi\lambda}(1-e^{-\dfrac{r^2}{4\alpha t}})+C_1\\\dfrac{\partial T}{\partial r}=(\dfrac{P_0}{2\pi\lambda}-C_1)\dfrac{1}{r}-\dfrac{P_0}{2\pi\lambda}\dfrac{e^{-\dfrac{r^2}{4\alpha t}}}{r}\\T(r,t)=(\dfrac{P_0}{2\pi\lambda}-C_1)ln(r)-\dfrac{P_0}{4\pi\lambda}\int_{0}^{\dfrac{r^2}{4\alpha t}}\dfrac{e^{-u}}{u}du+C_2\]
If the boundary condition is to be held (that is, the temperature field vanishes at infinity), we will need \(C_1=\dfrac{P_0}{2\pi\lambda}\) and \(C_2=\Gamma(0)\). Then we obtain the same result from the previous section.
In this assignment we plug \(u=\dfrac{r^2_{wall}}{4\alpha t}\) into the equation for a temperature time series when the heat sink is held constant. The initial temperature from assignment four is \(T_m=10.1917\) degree Celsius.
The convolution of the function is computed with the following algorithm: for every given \(t_i\), we find all the variations of heat sink \(P_j\) when \(j<i\). The effect of this \(P_j\) is \(T_{ij}=\dfrac{P_j}{4\pi\lambda}E_1(\dfrac{r^2}{4\alpha (t_i-t_j)})\). We then sum everything up as the following to get the total impact of the heat sink on the borewall: \[T_i=\Sigma_jT_{ij}\]
The resulting code for this scheme is:
lambda=2.5
rho=2300
c=900
alpha=lambda/rho/c
r_wall=0.07
T_m=10
time=seq(0,9)
P=c(0,.7,2.7,10.2,19.7,29.5,35.5,32.3,28.0,23.8)
diff_P=diff(P,1)
T_wall=c()
sum=c()
u=c()
for(i in 1:length(time)){
sum=0
if(i==1){
T_wall=T_m
}
else{
for(j in 1:(i-1)){
u=r_wall^2/4/alpha/((i-j)*30*86400)
sum=sum+diff_P[j]*expint_E1(u)/4/pi/lambda
}
T_wall[i]=T_m-sum
}
}We apply the scheme in the pevious section for a given variation of power demand per unit area. Again the P/A and \(T_m\) are again obtained from assignment 4. It is obvious from the graph below that a constant heat sink is about the moving average of the variable heat sink case.
According to the second law of Thermodynamics, the maximum heat able to transfer from lower temperature to higher temperature is \(Q_{carnot}=Q_{Demand}\dfrac{T_L}{T_H}\). The real available extracted heat is \(Q_{real}=Q_{Demand}-Q_{el}=Q_{Demand}-\dfrac{1}{1+\eta_{rel} \dfrac{T_L}{T_H-T_L}}Q_{Demand}=\dfrac{Q_{Demand}}{1+\dfrac{1}{\eta_{real}}\dfrac{T_H-T_L}{T_L}}\). Note that from the heat transient equations of a given heat source field, the temperature is also a dependent variable of \(Q_{real}\). This means that the equations must be solved dynamically. The easiest way of solving such system is to use a Euler forward scheme, that is \[\begin{cases} T_{t+1}=f(T_t,Q_{real,t})\\Q_{real,t}=g(T_t,Q_{real,t}) \end{cases}\]
Below we show the scheme used. Note that to obtain \(P_{geo}\) at initial time \(t(0)\), we have to define the temperature at a starting time \(t(-1)\).
lambda=2.5
rho=2300
c=900
alpha=lambda/rho/c
r_wall=0.07
eta_rel=.5
T_m=10
T_H=35
time=seq(0,10)-1
P_Demand=c(0,.7,2.7,10.2,19.7,29.5,35.5,32.3,28.0,23.8)
T_wall=c()
P_geo=c()
diff_P_geo=c()
sum=c()
u=c()
for(i in 1:length(time)){
if(i==1){
T_wall[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=0
for(j in 1:(i-1)){
u=r_wall^2/4/alpha/((i-j)*30*86400)
sum=sum+diff_P_geo[j]*expint_E1(u)/4/pi/lambda
}
T_wall[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]
}
}The resistance between temperature at infinity and origin can be written as \[\Delta T=P_0 \Sigma_i R_i=P_0(R_s+R_b)=P_0(\dfrac{E_1(\dfrac{r_b^2}{4\alpha t})}{4\pi\lambda}+R_b)\]
Since the resistance can be seen as connected as a series, we can conveniently find the equivalent radius to compute the temperature evolution when borehole resistance is included. \[E_1(\dfrac{r_a^2}{4\alpha t})=E_1(\dfrac{r_b^2}{4\alpha t})+R_b\\r_a \sim r_be^{-2\pi\lambda R_b}\]
From the graph below, it is obvious that when the heat demand is higher, the temperature difference between two models becomes larger.
The length per living area at the borehole length \(l\) is \(P_{geo}=\dfrac {Q_{geo}}{l}\). As we can see, for a period of thirty years, the minimum temperature occurs at 29.5 years. This means that we find the \(l\) such that \[T_{wall}(t=29.5;l)=0\]
As we can see in the graph below, the relationship is almost linear, and the minimum required length per area is about 0.4 meter.
The additional required power in the heat pump system is \(Q_{add}=Q_{demand}-Q_{geo}\). The annual additional energy required per area is the integral over the entire year of \(Q_{add}\).
In the scenario with borehole resistance, the annual electricity energy required per area grows from \(9.73 \dfrac{kWh}{m^2 a}\) at the beginning to \(10.20 \dfrac{kWh}{m^2 a}\) in the 30th year. This is only \(4.83\%\) of increase, or an annual increase rate of \(0.163\%\). The effect is thus very small.
Below we plot the additional energy required due to the borehole resistance. The effect of the borehole resistance is around \(5.0\sim5.5\%\) of the annual electricity energy demand per area, and remains fairly constant (though slightly decreases) in the time span.
The average COP over the thirty year time span is \(COP_{avg}=\dfrac{\Sigma_i Q_{demand}}{\Sigma_i Q_{el}}=4.9569>\dfrac{P_{el}}{P_{gas}}\).
Since the average COP is greater than the ratio of the average price between electricity and heating gas, it makes economic sense to build the geothermal heat pump.
Similar to what has been done in section h), we plot the annual electricity savings if a longer borehole is applied below. We can see that the savings increase as time goes on.
The total energy saved in the thirty year period is \(11.5731 kWh\), around \(3.82\%\) of the original electricity demand. This yields a financial saving (without considering depreciation) of \(0.28\times 60.0785=3.2405\) Euro per unit area. On the other hand, we have to increse the borehole length for \(\dfrac{1}{15}\) meter per unit area, which results to \(4.6667\) Euro of investment per unit area. So the investment does not pay off very well.
Again we plot the annual electricity savings if better insulation is applied below. We can see in this case the savings decrease as time goes on.
The total energy savings in the thirty year is \(3.1767 kWh\), which is only about \(1.05\%\) of the original electricity demand. This will save around \(0.8895\) Euro per area, but the investment costs \(1.6667\) Euro per area. So such investment also does not pay off.
Taking the results from 5(c), we can plot the graph below. From the graph we can see that \(\Delta T\) has a linear realtionship with \(P_l\).
We define the weighted mean temperature drop as \(\bar {\Delta T}=\dfrac{\Sigma P\Delta T}{\Sigma P}\). The resulting weighted mean temperature drop from 6(a) is \(9.1735\) degree Celsius.
We define average thermal resistance as \(\bar R=\dfrac{\bar {\Delta T}}{\bar P}\), where \(\bar P\) is the average power per length. The resulting average thermal resistance from 6(a) is \(0.5332 \dfrac{mK}{W}\).
We use the following equation to find the fluid temperature: \[T_f(z)=T_m(z)-a\dfrac{\partial T_m}{\partial z}+(T_f(0)-(T_m(0)-a\dfrac{\partial T_m}{\partial z}))e^{-\dfrac{z}{a}}\\a=R\rho_fc_fQ\]
As water, we have \(\rho_f=1000 kg/m^3\) and \(c_f=4184 J/(kg\cdot K)\). These value actuaaly changes as temperature changes, but let’s consider them constant for simplicity.
The power extraction per length at a given depth is \(P_l=\dfrac{T_m-T_f}{R}\).
The total power extraction is the integral of power extraction per length over z, which happens to be \[P_{total}=\rho_fc_fQ(T_f(z)-T_f(0))\]
The constant power demand for the system is \(P_{Demand}=50\times5000/8760=28.5388kW\).From the graph in the previous section, this is the total power extraction when the borehole length is about 1300 meter.
At 1510 meter, the temperature is about 35 degree Celsius. So the minimum borehole length is defined by the minimum required temperature.
Since the required power demand, \(P_{Demand}=50\times5000/8760=28.5388kW\), is a constant, we are looking for a solution to the optimum problem \[min[z^*];\\T_f(z^*)\geq35\\P(z^*)\geq P_{Demand}\]
The equation is over-complicated to solve directly, so we will instead find the minimum depth \(z^*\) that satisfies the minimum for a range of Q and then find a global minimum. As plotted below, we find the optimized flow rate to be 0.680 L/s, and the minimized depth to be 1317 meters.
Temperature reaches 25 degree Celsius at \(z=\dfrac{25-10.5}{0.05}=290\) meters. If we have insulation above this depth, that means \(T_f=25\) (degree Celsius) at this depth. So we can use the equation in 7(a) and the algorithm in 7(e) again, this time setting intital depth \(z_0=290\) (meter). The optimum flow rate is still 0.679 L/s, but the minimized depth is reduced to 1288 meters, with is a reduction of 29 meters (about 2.20%).
To deal with larger living area, we simply increase the power demand. For a living area of 10000 \(m^2\), the optimum flow rate is 1.362 L/s and the minimized depth is 1379 meters. For a living area of 2.724 L/s and 1920 meters.
Following the procedure provided above, we find the minimized required depth for a range of living area. Observe that optimum flow rate occurs at a flow rate of about \[\dfrac{P_{Demand}}{\rho_fc_fQ^*}\sim10\\Q^*\sim\dfrac{P_{Demand}}{10\rho_fc_f}\]
Which makes the computation much faster. We can then compute the LCoE of the geothermal heating system (assuming no depreciation and a 30 year lifetime as in assignment 5). As shown in the graph below, the system is economic viable when the living area is more then 22700 \(m^2\) under these assumptions.