# Calculation of Some Integrals Arising in Heat Transfer in Geothermics

Authors: J. L. G. Santander; P. Castañeda Porras; J. M. Isidro; P. Fernández de Córdoba
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

We calculate some integrals involved in the temperature field evaluation of the ground, heated by a borehole heat exchanger. This calculation allows a faster computation of that component of the temperature field which involves the periodic oscillation of the ambient temperature or the ambient heat flux.


---

## Body

## 1. Introduction

Ground-coupled heat pumps (GCHPs) are an attractive choice of system for heating and cooling buildings [1]. By comparison with standard technologies, these pumps offer competitive levels of comfort, reduced noise levels, lower greenhouse gas emissions, and reasonable environmental safety. Furthermore, their electrical consumption and maintenance requirements are lower than those required by conventional systems and, consequently, have lower annual operating cost [2].A GCHP exchanges heat with the ground through a buried U tube loop. Since this exchange strongly depends on the thermal properties of the ground, it is very important to have a knowledge of these properties when designing GCHP air-conditioning systems. Methods to estimate them include literature searches, conducting laboratory experiments on soil/rock samples, and/or performing field tests. Thein situ tests are based on studying the thermal response of the borehole heat exchanger to a constant heat injection or extraction. The outputs of the thermal response test are the inlet and outlet temperature of the heat-carrier fluid as a function of time. From these experimental data, and with an appropriate model describing the heat transfer between the fluid and the ground, the thermal conductivity of the surroundings is inferred.The results presented in this study are based on a model of the underground heat transfer due to a borehole heat exchanger assuming that it works as a line heat source of finite depth inside the ground. We assume as well a constant power of length unit for the heat source considered and a periodic oscillation of the ambient temperature or heat flux on the ground surface.This paper is organized as follows. Section2 describes the equations that model the heat transfer inside the ground due to a borehole heat exchanger. Section 3 solves the integrals given in literature for periodic oscillations of the ambient temperature or heat flux, on the ground surface. Section 4 shows that our results agree with the literature in some particular cases, such as the solution on the ground surface and the solution of the quasistationary regime.


## 2. The Equations

Let us consider the heat equation in cylindrical coordinates, in whichT(r,z,t) is the temperature field,(2.1)∂T∂t-k∇2T=QzCδ(r)[H(z)-H(z-L)].
On the one hand, the thermal parameters of (2.1) are the ground thermal diffusivity k(m2s-1), the heat flow per length unit Qz(Wm-1), and the volumetric heat capacity of ground C(Jm-3K-1). On the other hand, the geometry of the heat source, depicted in Figure 1, has been taken into account introducing the Heaviside function H for its finite length L and the Dirac delta distribution δ for its infinitely thin width. Equation (2.1) is subject to to the initial condition(2.2)T(r,z,0)=T0+kgeoz,z≥0,
where T0 is the undisturbed ground temperature and kgeo(Km-1) is the geothermal gradient. We may consider that (2.1) is subject to a Dirichlet boundary condition, which represents a time-dependent temperature on the surface z=0,(2.3)T(r,0,t)=ψ(t),t≥0,
or to a Neumann boundary condition, which represents a time-dependent heat flux on the surface z=0,(2.4)-λ∂zT(r,0,t)=ψ(t),t≥0,
where λ(Wm-1K-1) is the ground thermal conductivity. We split the Cauchy problem posed in (2.1), (2.2), and (2.3) or (2.4) in three terms [3], so that(2.5)T(r,z,t)=vd(r,z,t)+v0(z,t)+vs(z,t),
and where the functions v0, vs and vd satisfy the Cauchy problems presented below.Figure 1
Borehole heat exchanger model.


### 2.1. The Problem forvd

The partial differential equation thatvd(r,z,t) satisfies is (2.6)∂vd∂t-k∇2vd=QzCδ(r)[H(z)-H(z-L)],
subject to a homogeneous initial condition,(2.7)vd(r,z,0)=0,z≥0,
and to a homogeneous Dirichlet boundary condition,(2.8)vd(r,0,t)=0,t≥0,
or to a Neumann boundary condition, homogenous as well,(2.9)∂zvd(r,0,t)=0,t≥0.
The solution to (2.6), (2.7), and (2.8) is [4](2.10)vd(r,z,t)=Qz8πλ∫r2/4kt∞due-uu×{2erf⁡(zru)+erf⁡(L-zru)-erf⁡(L+zru)},
and the solution to (2.6), (2.7), and (2.9) is (2.11)vd(r,z,t)=Qz8πλ∫r2/4kt∞e-uu{erf⁡(z+Lru)-erf⁡(z-Lru)}du,
that we will present in a future work.


### 2.2. The Problem forv0

The partial differential equation thatv0(z,t) satisfies is (2.12)∂v0∂t=k∂2v0∂z2,
subject to the initial condition (2.13)v0(z,0)=T0+kgeoz,z≥0,
and to a Dirichlet boundary condition, (2.14)v0(0,t)=0,t≥0,
or to a Neumann homogeneous boundary condition (2.15)∂zv0(0,t)=0,t≥0.
The solution to (2.12), (2.13), and (2.14) is [5, Equation 2.4.13],(2.16)v0(z,t)=T0erf⁡(z2kt)+kgeoz,
and the solution to (2.12), (2.13), and (2.15) is [5, Equation 2.9.7] (2.17)v0(z,t)=T0+kgeo{zerf⁡(z2kt)+2ktπexp⁡(-z24kt)},
where the error function is defined as (2.18)erf⁡(z):=2π∫0ze-u2du.


### 2.3. The Problem forvs

The partial differential equation thatvs(z,t) satisfies is(2.19)∂vs∂t=k∂2vs∂z2,
subject to an homogeneous initial condition(2.20)vs(z,0)=0,z≥0,
and to the time-dependent Dirichlet boundary condition (2.21)vs(0,t)=ψ(t),t≥0,
or to the time-dependent Neumann boundary condition (2.22)-λ∂zvs(0,t)=ψ(t),t≥0.
The solution to (2.19), (2.20) and (2.21) is [5, Equation 2.5.1.] (2.23)vs(z,t)=2π∫z/2kt∞ψ(t-z24ku2)e-u2du,
and the solution to (2.19), (2.20) and (2.22) is [5, Equation 2.9.9] (2.24)vs(z,t)=kλπ∫0tψ(t-u)ue-z2/4kudu.


## 2.1. The Problem forvd

The partial differential equation thatvd(r,z,t) satisfies is (2.6)∂vd∂t-k∇2vd=QzCδ(r)[H(z)-H(z-L)],
subject to a homogeneous initial condition,(2.7)vd(r,z,0)=0,z≥0,
and to a homogeneous Dirichlet boundary condition,(2.8)vd(r,0,t)=0,t≥0,
or to a Neumann boundary condition, homogenous as well,(2.9)∂zvd(r,0,t)=0,t≥0.
The solution to (2.6), (2.7), and (2.8) is [4](2.10)vd(r,z,t)=Qz8πλ∫r2/4kt∞due-uu×{2erf⁡(zru)+erf⁡(L-zru)-erf⁡(L+zru)},
and the solution to (2.6), (2.7), and (2.9) is (2.11)vd(r,z,t)=Qz8πλ∫r2/4kt∞e-uu{erf⁡(z+Lru)-erf⁡(z-Lru)}du,
that we will present in a future work.


## 2.2. The Problem forv0

The partial differential equation thatv0(z,t) satisfies is (2.12)∂v0∂t=k∂2v0∂z2,
subject to the initial condition (2.13)v0(z,0)=T0+kgeoz,z≥0,
and to a Dirichlet boundary condition, (2.14)v0(0,t)=0,t≥0,
or to a Neumann homogeneous boundary condition (2.15)∂zv0(0,t)=0,t≥0.
The solution to (2.12), (2.13), and (2.14) is [5, Equation 2.4.13],(2.16)v0(z,t)=T0erf⁡(z2kt)+kgeoz,
and the solution to (2.12), (2.13), and (2.15) is [5, Equation 2.9.7] (2.17)v0(z,t)=T0+kgeo{zerf⁡(z2kt)+2ktπexp⁡(-z24kt)},
where the error function is defined as (2.18)erf⁡(z):=2π∫0ze-u2du.


## 2.3. The Problem forvs

The partial differential equation thatvs(z,t) satisfies is(2.19)∂vs∂t=k∂2vs∂z2,
subject to an homogeneous initial condition(2.20)vs(z,0)=0,z≥0,
and to the time-dependent Dirichlet boundary condition (2.21)vs(0,t)=ψ(t),t≥0,
or to the time-dependent Neumann boundary condition (2.22)-λ∂zvs(0,t)=ψ(t),t≥0.
The solution to (2.19), (2.20) and (2.21) is [5, Equation 2.5.1.] (2.23)vs(z,t)=2π∫z/2kt∞ψ(t-z24ku2)e-u2du,
and the solution to (2.19), (2.20) and (2.22) is [5, Equation 2.9.9] (2.24)vs(z,t)=kλπ∫0tψ(t-u)ue-z2/4kudu.


## 3. Harmonic Analysis of the Boundary Condition

We may consider that on the ground surface we have a periodic temperature or heat flux oscillation, due to daily or annual cycles. Therefore,ψ(t) is a periodic function of a certain period τ, continuous and bounded, that we may expand in Fourier series(3.1)ψ(t)=a02+∑n=1∞ancos⁡(nωt)+bnsin⁡(nωt),
where the frequency is ω=2π/τ, and the coefficients an and bn are given by (3.2)a0=2τ∫-τ/2τ/2ψ(u)du,an=2τ∫-τ/2τ/2ψ(u)cos⁡(nωu)du,bn=2τ∫-τ/2τ/2ψ(u)sin⁡(nωu)du,
where n∈ℕ.


### 3.1. Dirichlet Boundary Condition

Substituting (3.1) in (2.23), we obtain (3.3)vs(z,t)=a02erfc⁡(z2kt)+2π∑n=1∞an∫z/2kt∞cos⁡(nωt-nωz24ku2)e-u2du+bn∫z/2kt∞sin⁡(nωt-nωz24ku2)e-u2du,
where erfc(z)=1-erf(z) is the complementary error function. The integrals given in (3.3) for t≥0 may be rewritten as (3.4)∫z/2kt∞cos⁡(nωt-nωz24ku2)e-u2du=Re[einωtID,n(z,t)],∫z/2kt∞sin⁡(nωt-nωz24ku2)e-u2du=Im⁡[einωtID,n(z,t)],
where we have defined(3.5)ID,n(z,t):=∫z/2kt∞exp⁡(-inωz24ku2-u2)du.
In order to calculate (3.5), let us apply the integral [6, Equation 7.4.33] (3.6)∫exp⁡(-ax2-bx2)dx=π4a[e2aberf⁡(ax+bx)+e-2aberf⁡(ax-bx)]+const,
taking a=1, b=inωz2/4k, and z≥0, so that(3.7)ID,n(z,t)=π4[exp⁡(zinωk)erfc⁡(z2kt+inωt)+exp⁡(-zinωk)erfc⁡(z2kt-inωt)].
Finally, substituting (3.7) in (3.4), we obtain for z,t≥0,(3.8)∫z/2kt∞cos⁡(nωt-nωz24ku2)e-u2du=π4Re{einωt[An+(z,t)+An-(z,t)]},∫z/2kt∞sin⁡(nωt-nωz24ku2)e-u2du=π4Im⁡{einωt[An+(z,t)+An-(z,t)]},
where we have defined(3.9)An±(z,t):=exp⁡(±zinωk)erfc⁡(z2kt±inωt).
As far as we know, the integrals given in (3.8) are not reported in the literature.


#### 3.1.1. Neumann Boundary Condition

Substituting (3.1) in (2.24), we obtain (3.10)vs(z,t)=a0k2λπ∫0te-z2/4kuudu(3.11)+kλπ∑n=1∞an∫0tcos⁡[nω(t-u)]ue-z2/4kudu(3.12)+bn∫0tsin⁡[nω(t-u)]ue-z2/4kudu.
Integrating by parts and performing the substitution ζ=z/2ku, considering z,t≥0, we may calculate the integral given in (3.10) (3.13)∫0te-z2/4kuudu=2texp⁡(-z24kt)-πzkerfc⁡(z2kt).
In order to calculate the integrals given in (3.11) and (3.12), let us perform the substitution u=v2, considering t≥0, so that(3.14)∫0tcos⁡[nω(t-u)]ue-z2/4kudu=2∫0tcos⁡[nω(t-v2)]e-z2/4kv2dv=2Re[einωtIN,n(z,t)],
and similarly(3.15)∫0tsin⁡[nω(t-u)]ue-z2/4kudu=2Im⁡[einωtIN,n(z,t)],
where we have defined(3.16)IN,n(z,t):=∫0texp⁡(-inωv2-z24kv2)dv.
Taking a=inω, b=z2/4k, and z≥0 in (3.6), we may rewrite (3.16) as(3.17)IN,n(z,t)=π4inω[exp⁡(-zinωk)erfc⁡(z2kt-inωt)-exp⁡(zinωk)erfc⁡(z2kt+inωt)].
Finally, substituting (3.17) in (3.14) and (3.15), for z,t≥0, we obtain (3.18)∫0tcos⁡[nω(t-u)]ue-z2/4kudu=π2Re{einωtinω[An-(z,t)-An+(z,t)]},∫0tsin⁡[nω(t-u)]ue-z2/4kudu=π2Im⁡{einωtinω[An-(z,t)-An+(z,t)]}.
As far as we know, the integrals given in (3.18) are not reported in the literature.


## 3.1. Dirichlet Boundary Condition

Substituting (3.1) in (2.23), we obtain (3.3)vs(z,t)=a02erfc⁡(z2kt)+2π∑n=1∞an∫z/2kt∞cos⁡(nωt-nωz24ku2)e-u2du+bn∫z/2kt∞sin⁡(nωt-nωz24ku2)e-u2du,
where erfc(z)=1-erf(z) is the complementary error function. The integrals given in (3.3) for t≥0 may be rewritten as (3.4)∫z/2kt∞cos⁡(nωt-nωz24ku2)e-u2du=Re[einωtID,n(z,t)],∫z/2kt∞sin⁡(nωt-nωz24ku2)e-u2du=Im⁡[einωtID,n(z,t)],
where we have defined(3.5)ID,n(z,t):=∫z/2kt∞exp⁡(-inωz24ku2-u2)du.
In order to calculate (3.5), let us apply the integral [6, Equation 7.4.33] (3.6)∫exp⁡(-ax2-bx2)dx=π4a[e2aberf⁡(ax+bx)+e-2aberf⁡(ax-bx)]+const,
taking a=1, b=inωz2/4k, and z≥0, so that(3.7)ID,n(z,t)=π4[exp⁡(zinωk)erfc⁡(z2kt+inωt)+exp⁡(-zinωk)erfc⁡(z2kt-inωt)].
Finally, substituting (3.7) in (3.4), we obtain for z,t≥0,(3.8)∫z/2kt∞cos⁡(nωt-nωz24ku2)e-u2du=π4Re{einωt[An+(z,t)+An-(z,t)]},∫z/2kt∞sin⁡(nωt-nωz24ku2)e-u2du=π4Im⁡{einωt[An+(z,t)+An-(z,t)]},
where we have defined(3.9)An±(z,t):=exp⁡(±zinωk)erfc⁡(z2kt±inωt).
As far as we know, the integrals given in (3.8) are not reported in the literature.


### 3.1.1. Neumann Boundary Condition

Substituting (3.1) in (2.24), we obtain (3.10)vs(z,t)=a0k2λπ∫0te-z2/4kuudu(3.11)+kλπ∑n=1∞an∫0tcos⁡[nω(t-u)]ue-z2/4kudu(3.12)+bn∫0tsin⁡[nω(t-u)]ue-z2/4kudu.
Integrating by parts and performing the substitution ζ=z/2ku, considering z,t≥0, we may calculate the integral given in (3.10) (3.13)∫0te-z2/4kuudu=2texp⁡(-z24kt)-πzkerfc⁡(z2kt).
In order to calculate the integrals given in (3.11) and (3.12), let us perform the substitution u=v2, considering t≥0, so that(3.14)∫0tcos⁡[nω(t-u)]ue-z2/4kudu=2∫0tcos⁡[nω(t-v2)]e-z2/4kv2dv=2Re[einωtIN,n(z,t)],
and similarly(3.15)∫0tsin⁡[nω(t-u)]ue-z2/4kudu=2Im⁡[einωtIN,n(z,t)],
where we have defined(3.16)IN,n(z,t):=∫0texp⁡(-inωv2-z24kv2)dv.
Taking a=inω, b=z2/4k, and z≥0 in (3.6), we may rewrite (3.16) as(3.17)IN,n(z,t)=π4inω[exp⁡(-zinωk)erfc⁡(z2kt-inωt)-exp⁡(zinωk)erfc⁡(z2kt+inωt)].
Finally, substituting (3.17) in (3.14) and (3.15), for z,t≥0, we obtain (3.18)∫0tcos⁡[nω(t-u)]ue-z2/4kudu=π2Re{einωtinω[An-(z,t)-An+(z,t)]},∫0tsin⁡[nω(t-u)]ue-z2/4kudu=π2Im⁡{einωtinω[An-(z,t)-An+(z,t)]}.
As far as we know, the integrals given in (3.18) are not reported in the literature.


## 3.1.1. Neumann Boundary Condition

Substituting (3.1) in (2.24), we obtain (3.10)vs(z,t)=a0k2λπ∫0te-z2/4kuudu(3.11)+kλπ∑n=1∞an∫0tcos⁡[nω(t-u)]ue-z2/4kudu(3.12)+bn∫0tsin⁡[nω(t-u)]ue-z2/4kudu.
Integrating by parts and performing the substitution ζ=z/2ku, considering z,t≥0, we may calculate the integral given in (3.10) (3.13)∫0te-z2/4kuudu=2texp⁡(-z24kt)-πzkerfc⁡(z2kt).
In order to calculate the integrals given in (3.11) and (3.12), let us perform the substitution u=v2, considering t≥0, so that(3.14)∫0tcos⁡[nω(t-u)]ue-z2/4kudu=2∫0tcos⁡[nω(t-v2)]e-z2/4kv2dv=2Re[einωtIN,n(z,t)],
and similarly(3.15)∫0tsin⁡[nω(t-u)]ue-z2/4kudu=2Im⁡[einωtIN,n(z,t)],
where we have defined(3.16)IN,n(z,t):=∫0texp⁡(-inωv2-z24kv2)dv.
Taking a=inω, b=z2/4k, and z≥0 in (3.6), we may rewrite (3.16) as(3.17)IN,n(z,t)=π4inω[exp⁡(-zinωk)erfc⁡(z2kt-inωt)-exp⁡(zinωk)erfc⁡(z2kt+inωt)].
Finally, substituting (3.17) in (3.14) and (3.15), for z,t≥0, we obtain (3.18)∫0tcos⁡[nω(t-u)]ue-z2/4kudu=π2Re{einωtinω[An-(z,t)-An+(z,t)]},∫0tsin⁡[nω(t-u)]ue-z2/4kudu=π2Im⁡{einωtinω[An-(z,t)-An+(z,t)]}.
As far as we know, the integrals given in (3.18) are not reported in the literature.


## 4. Particular Cases

As a consistency check, let us verify that in certain particular cases the new integrals given in (3.8) and (3.18) are reduced to integrals reported in the literature. These particular cases have to do with the quasistationary regime and the solution on the ground surface.


### 4.1. Quasistationary Regime

Since there is a periodic oscillation on the ground surface, the surroundings of the borehole heat exchanger never reach a stationary regime. However, according to [5, Section 2.6], we may define a quasistationary regime in which the ground temperature field is stabilized periodically. In this quasistationary regime, (3.8) becomes (4.1)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π4Re{einωtAD,n(z)},(4.2)∫0∞sin⁡(nωt-nωz24ku2)e-u2du=π4Im⁡{einωtAD,n(z)},
where we have defined(4.3)AD,n(z):=lim⁡t→∞[An+(z,t)+An-(z,t)].
Similarly, (3.18) becomes (4.4)∫0∞cos⁡[nω(t-u)]ue-z2/4kudu=π2Re{einωtinωAN,n(z)},∫0∞sin⁡[nω(t-u)]ue-z2/4kudu=π2Im⁡{einωtinωAN,n(z)},
where we have defined(4.5)AN,n(z):=lim⁡t→∞[An-(z,t)-An+(z,t)].


#### 4.1.1. Dirichlet Case

Taking into account thaterfc(∞)=0 and erfc(-∞)=2, we may calculate the limit given in (4.3), so that (4.1) becomes(4.6)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π2Re{einωtexp⁡(-zinωk)}=π2e-znω/2kRe{ei(nωt-znω/2k)},
where we have substitute i=(1+i)/2. Thus, for z,t≥0,(4.7)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π2e-znω/2kcos⁡(nωt-znω2k).
Similarly, we may obtain(4.8)∫0∞sin⁡(nωt-nωz24ku2)e-u2du=π2e-znω/2ksin⁡(nωt-znω2k).
Notice that (4.7) and (4.8) are sinusoidal in t, with an amplitude exponentially decreasing in depth z. In the literature, we may find [7, Equation 3.928.1-2]. (4.9)∫0∞e-p2x2-q2/x2sin⁡(a2x2+b2x2)dx=π2re-2rscos⁡(A+B)×sin⁡(A+2rssin⁡(A+B)),∫0∞e-p2x2-q2/x2cos⁡(a2x2+b2x2)dx=π2re-2rscos⁡(A+B)×cos⁡(A+2rssin⁡(A+B)),
where a2+p2>0 and (4.10)r=a4+p44,s=b4+q44,A=12tan-1(a2p2),B=12tan-1(b2q2),
thus, taking p=1, q=0, and a=0,(4.11)∫0∞e-x2sin⁡(b2x2)dx=π2e-2bsin⁡(2b),∫0∞e-x2cos⁡(b2x2)dx=π2e-2bcos⁡(2b).
Therefore, rewriting the left side of (4.7) as(4.12)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=cos⁡(nωt)∫0∞cos⁡(nωz24ku2)e-u2du+sin⁡(nωt)∫0∞sin⁡(nωz24ku2)e-u2du,
and applying (4.11), taking b=nωz/2k, (z≥0), we eventually get the same result as (4.7). Similarly, we may check the result given in (4.8).


#### 4.1.2. Neumann Case

In a similar way as in the previous subsection, we may calculate (4.4) for z,t≥0, arriving at (4.13)∫0∞cos⁡[nω(t-u)]ue-z2/4kudu=πe-znω/2knωcos⁡(nωt-znω2k-π4),(4.14)∫0∞sin⁡[nω(t-u)]ue-z2/4kudu=πe-znω/2knωsin⁡(nωt-znω2k-π4).
Once again, (4.13) and (4.7) are sinusoidal in t, with an amplitude exponentially decreasing in depth z. In the literature, we may find [7, Equation 3.957.1-2]
(4.15)∫0∞xμ-1e-β2/4xcos⁡axdx=(β2a)μ[e-iπμ/4Kμ(βeiπ/4a)+eiπμ/4Kμ(βe-iπ/4a)],∫0∞xμ-1e-β2/4xsin⁡axdx=i(β2a)μ[e-iπμ/4Kμ(βeiπ/4a)-eiπμ/4Kμ(βe-iπ/4a)],
where Reβ>0, Reμ<1 and a>0. Therefore, taking in (4.15),μ=1/2, a=nω, β=z/k, (z≥0), and knowing that the Macdonald function of order 1/2 is [8, Equation 5.5.5], (4.16)K1/2(z)=π2ze-z,
we eventually get the same results as (4.13) and (4.14).


### 4.2. Solution on the Ground Surface

Notice that in the special casez=0, (3.8) becomes trivial, (4.17)∫0∞cos⁡(nωt)e-u2du=π2cos⁡(nωt),∫0∞sin⁡(nωt)e-u2du=π2sin⁡(nωt),
while (3.18) are reduced to (4.18)∫0tcos⁡[nω(t-u)]udu=πRe{einωtinωerf⁡(inωt)},∫0tsin⁡[nω(t-u)]udu=πIm⁡{einωtinωerf⁡(inωt)}.
Since [8, Equation 2.4.1] (4.19)erf⁡(ix)i=2[C(2πx)-iS(2πx)],x∈R,
where (4.20)C(z):=∫0zcos⁡πt22dt,S(z):=∫0zsin⁡πt22dt
are the Fresnel integrals; then (4.18) may be expressed as (4.21)∫0tcos⁡[nω(t-u)]udu=2πnω[cos⁡(nωt)C(2nωtπ)+sin⁡(nωt)S(2nωtπ)],∫0tsin⁡[nω(t-u)]udu=2πnω[sin⁡(nωt)C(2nωtπ)-cos⁡(nωt)S(2nωtπ)].
The results given in (4.21), agrees with [7, Equation 2.653.1-2] (4.22)∫sin⁡xxdx=2πS(x),∫cos⁡xxdx=2πC(x).


## 4.1. Quasistationary Regime

Since there is a periodic oscillation on the ground surface, the surroundings of the borehole heat exchanger never reach a stationary regime. However, according to [5, Section 2.6], we may define a quasistationary regime in which the ground temperature field is stabilized periodically. In this quasistationary regime, (3.8) becomes (4.1)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π4Re{einωtAD,n(z)},(4.2)∫0∞sin⁡(nωt-nωz24ku2)e-u2du=π4Im⁡{einωtAD,n(z)},
where we have defined(4.3)AD,n(z):=lim⁡t→∞[An+(z,t)+An-(z,t)].
Similarly, (3.18) becomes (4.4)∫0∞cos⁡[nω(t-u)]ue-z2/4kudu=π2Re{einωtinωAN,n(z)},∫0∞sin⁡[nω(t-u)]ue-z2/4kudu=π2Im⁡{einωtinωAN,n(z)},
where we have defined(4.5)AN,n(z):=lim⁡t→∞[An-(z,t)-An+(z,t)].


### 4.1.1. Dirichlet Case

Taking into account thaterfc(∞)=0 and erfc(-∞)=2, we may calculate the limit given in (4.3), so that (4.1) becomes(4.6)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π2Re{einωtexp⁡(-zinωk)}=π2e-znω/2kRe{ei(nωt-znω/2k)},
where we have substitute i=(1+i)/2. Thus, for z,t≥0,(4.7)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π2e-znω/2kcos⁡(nωt-znω2k).
Similarly, we may obtain(4.8)∫0∞sin⁡(nωt-nωz24ku2)e-u2du=π2e-znω/2ksin⁡(nωt-znω2k).
Notice that (4.7) and (4.8) are sinusoidal in t, with an amplitude exponentially decreasing in depth z. In the literature, we may find [7, Equation 3.928.1-2]. (4.9)∫0∞e-p2x2-q2/x2sin⁡(a2x2+b2x2)dx=π2re-2rscos⁡(A+B)×sin⁡(A+2rssin⁡(A+B)),∫0∞e-p2x2-q2/x2cos⁡(a2x2+b2x2)dx=π2re-2rscos⁡(A+B)×cos⁡(A+2rssin⁡(A+B)),
where a2+p2>0 and (4.10)r=a4+p44,s=b4+q44,A=12tan-1(a2p2),B=12tan-1(b2q2),
thus, taking p=1, q=0, and a=0,(4.11)∫0∞e-x2sin⁡(b2x2)dx=π2e-2bsin⁡(2b),∫0∞e-x2cos⁡(b2x2)dx=π2e-2bcos⁡(2b).
Therefore, rewriting the left side of (4.7) as(4.12)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=cos⁡(nωt)∫0∞cos⁡(nωz24ku2)e-u2du+sin⁡(nωt)∫0∞sin⁡(nωz24ku2)e-u2du,
and applying (4.11), taking b=nωz/2k, (z≥0), we eventually get the same result as (4.7). Similarly, we may check the result given in (4.8).


### 4.1.2. Neumann Case

In a similar way as in the previous subsection, we may calculate (4.4) for z,t≥0, arriving at (4.13)∫0∞cos⁡[nω(t-u)]ue-z2/4kudu=πe-znω/2knωcos⁡(nωt-znω2k-π4),(4.14)∫0∞sin⁡[nω(t-u)]ue-z2/4kudu=πe-znω/2knωsin⁡(nωt-znω2k-π4).
Once again, (4.13) and (4.7) are sinusoidal in t, with an amplitude exponentially decreasing in depth z. In the literature, we may find [7, Equation 3.957.1-2]
(4.15)∫0∞xμ-1e-β2/4xcos⁡axdx=(β2a)μ[e-iπμ/4Kμ(βeiπ/4a)+eiπμ/4Kμ(βe-iπ/4a)],∫0∞xμ-1e-β2/4xsin⁡axdx=i(β2a)μ[e-iπμ/4Kμ(βeiπ/4a)-eiπμ/4Kμ(βe-iπ/4a)],
where Reβ>0, Reμ<1 and a>0. Therefore, taking in (4.15),μ=1/2, a=nω, β=z/k, (z≥0), and knowing that the Macdonald function of order 1/2 is [8, Equation 5.5.5], (4.16)K1/2(z)=π2ze-z,
we eventually get the same results as (4.13) and (4.14).


## 4.1.1. Dirichlet Case

Taking into account thaterfc(∞)=0 and erfc(-∞)=2, we may calculate the limit given in (4.3), so that (4.1) becomes(4.6)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π2Re{einωtexp⁡(-zinωk)}=π2e-znω/2kRe{ei(nωt-znω/2k)},
where we have substitute i=(1+i)/2. Thus, for z,t≥0,(4.7)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=π2e-znω/2kcos⁡(nωt-znω2k).
Similarly, we may obtain(4.8)∫0∞sin⁡(nωt-nωz24ku2)e-u2du=π2e-znω/2ksin⁡(nωt-znω2k).
Notice that (4.7) and (4.8) are sinusoidal in t, with an amplitude exponentially decreasing in depth z. In the literature, we may find [7, Equation 3.928.1-2]. (4.9)∫0∞e-p2x2-q2/x2sin⁡(a2x2+b2x2)dx=π2re-2rscos⁡(A+B)×sin⁡(A+2rssin⁡(A+B)),∫0∞e-p2x2-q2/x2cos⁡(a2x2+b2x2)dx=π2re-2rscos⁡(A+B)×cos⁡(A+2rssin⁡(A+B)),
where a2+p2>0 and (4.10)r=a4+p44,s=b4+q44,A=12tan-1(a2p2),B=12tan-1(b2q2),
thus, taking p=1, q=0, and a=0,(4.11)∫0∞e-x2sin⁡(b2x2)dx=π2e-2bsin⁡(2b),∫0∞e-x2cos⁡(b2x2)dx=π2e-2bcos⁡(2b).
Therefore, rewriting the left side of (4.7) as(4.12)∫0∞cos⁡(nωt-nωz24ku2)e-u2du=cos⁡(nωt)∫0∞cos⁡(nωz24ku2)e-u2du+sin⁡(nωt)∫0∞sin⁡(nωz24ku2)e-u2du,
and applying (4.11), taking b=nωz/2k, (z≥0), we eventually get the same result as (4.7). Similarly, we may check the result given in (4.8).


## 4.1.2. Neumann Case

In a similar way as in the previous subsection, we may calculate (4.4) for z,t≥0, arriving at (4.13)∫0∞cos⁡[nω(t-u)]ue-z2/4kudu=πe-znω/2knωcos⁡(nωt-znω2k-π4),(4.14)∫0∞sin⁡[nω(t-u)]ue-z2/4kudu=πe-znω/2knωsin⁡(nωt-znω2k-π4).
Once again, (4.13) and (4.7) are sinusoidal in t, with an amplitude exponentially decreasing in depth z. In the literature, we may find [7, Equation 3.957.1-2]
(4.15)∫0∞xμ-1e-β2/4xcos⁡axdx=(β2a)μ[e-iπμ/4Kμ(βeiπ/4a)+eiπμ/4Kμ(βe-iπ/4a)],∫0∞xμ-1e-β2/4xsin⁡axdx=i(β2a)μ[e-iπμ/4Kμ(βeiπ/4a)-eiπμ/4Kμ(βe-iπ/4a)],
where Reβ>0, Reμ<1 and a>0. Therefore, taking in (4.15),μ=1/2, a=nω, β=z/k, (z≥0), and knowing that the Macdonald function of order 1/2 is [8, Equation 5.5.5], (4.16)K1/2(z)=π2ze-z,
we eventually get the same results as (4.13) and (4.14).


## 4.2. Solution on the Ground Surface

Notice that in the special casez=0, (3.8) becomes trivial, (4.17)∫0∞cos⁡(nωt)e-u2du=π2cos⁡(nωt),∫0∞sin⁡(nωt)e-u2du=π2sin⁡(nωt),
while (3.18) are reduced to (4.18)∫0tcos⁡[nω(t-u)]udu=πRe{einωtinωerf⁡(inωt)},∫0tsin⁡[nω(t-u)]udu=πIm⁡{einωtinωerf⁡(inωt)}.
Since [8, Equation 2.4.1] (4.19)erf⁡(ix)i=2[C(2πx)-iS(2πx)],x∈R,
where (4.20)C(z):=∫0zcos⁡πt22dt,S(z):=∫0zsin⁡πt22dt
are the Fresnel integrals; then (4.18) may be expressed as (4.21)∫0tcos⁡[nω(t-u)]udu=2πnω[cos⁡(nωt)C(2nωtπ)+sin⁡(nωt)S(2nωtπ)],∫0tsin⁡[nω(t-u)]udu=2πnω[sin⁡(nωt)C(2nωtπ)-cos⁡(nωt)S(2nωtπ)].
The results given in (4.21), agrees with [7, Equation 2.653.1-2] (4.22)∫sin⁡xxdx=2πS(x),∫cos⁡xxdx=2πC(x).



---

*Source: 784794-2010-12-20.xml*