Calculation of Some Integrals Arising in Heat Transfer in Geothermics

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.

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. The in 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. Section 2 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 which 𝑇(𝑟,𝑧,𝑡) is the temperature field,𝜕𝑇𝜕𝑡−𝑘∇²𝑄𝑇=^𝑧𝐶𝛿(𝑟)[𝐻(𝑧)−𝐻(𝑧−𝐿)]. (2.1) On the one hand, the thermal parameters of (2.1) are the ground thermal diffusivity 𝑘(𝑚²𝑠⁻¹), the heat flow per length unit 𝑄_𝑧(𝑊𝑚⁻¹), and the volumetric heat capacity of ground 𝐶(𝐽𝑚⁻³𝐾⁻¹). On the other hand, the geometry of the heat source, depicted in Figure 1, has been taken into account introducing the Heaviside function 𝐻 for its finite length 𝐿 and the Dirac delta distribution 𝛿 for its infinitely thin width. Equation (2.1) is subject to to the initial condition𝑇(𝑟,𝑧,0)=𝑇₀+𝑘_geo𝑧,𝑧≥0, (2.2) where 𝑇₀ is the undisturbed ground temperature and 𝑘_geo(𝐾𝑚⁻¹) 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 𝑧=0,𝑇(𝑟,0,𝑡)=𝜓(𝑡),𝑡≥0, (2.3) or to a Neumann boundary condition, which represents a time-dependent heat flux on the surface 𝑧=0,−𝜆𝜕_𝑧𝑇(𝑟,0,𝑡)=𝜓(𝑡),𝑡≥0, (2.4) where 𝜆(𝑊𝑚⁻¹𝐾⁻¹) 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) and where the functions 𝑣₀, 𝑣_𝑠 and 𝑣_𝑑 satisfy the Cauchy problems presented below.

Figure 1 

Borehole heat exchanger model.

2.1. The Problem for 𝑣_𝑑

The partial differential equation that 𝑣_𝑑(𝑟,𝑧,𝑡) satisfies is 𝜕𝑣_𝑑𝜕𝑡−𝑘∇²𝑣_𝑑=𝑄_𝑧𝐶𝛿(𝑟)[𝐻(𝑧)−𝐻(𝑧−𝐿)], (2.6) subject to a homogeneous initial condition,𝑣_𝑑(𝑟,𝑧,0)=0,𝑧≥0, (2.7) and to a homogeneous Dirichlet boundary condition,𝑣_𝑑(𝑟,0,𝑡)=0,𝑡≥0, (2.8) or to a Neumann boundary condition, homogenous as well,𝜕_𝑧𝑣_𝑑(𝑟,0,𝑡)=0,𝑡≥0. (2.9) The solution to (2.6), (2.7), and (2.8) is [4]𝑣_𝑑()=𝑄𝑟,𝑧,𝑡^𝑧8𝜋𝜆^∞_𝑟²_/4𝑘𝑡𝑒𝑑𝑢^−𝑢𝑢×𝑧2erf𝑟√𝑢+erf𝐿−𝑧𝑟√𝑢−erf𝐿+𝑧𝑟√𝑢, (2.10) and the solution to (2.6), (2.7), and (2.9) is 𝑣_𝑑()=𝑄𝑟,𝑧,𝑡^𝑧8𝜋𝜆^∞_𝑟²_/4𝑘𝑡𝑒^−𝑢𝑢erf𝑧+𝐿𝑟√𝑢−erf𝑧−𝐿𝑟√𝑢𝑑𝑢, (2.11) that we will present in a future work.

2.2. The Problem for 𝑣₀

The partial differential equation that 𝑣₀(𝑧,𝑡) satisfies is 𝜕𝑣₀𝜕𝜕𝑡=𝑘²𝑣₀𝜕𝑧², (2.12) subject to the initial condition 𝑣₀(𝑧,0)=𝑇₀+𝑘_geo𝑧,𝑧≥0, (2.13) and to a Dirichlet boundary condition, 𝑣₀(0,𝑡)=0,𝑡≥0, (2.14) or to a Neumann homogeneous boundary condition 𝜕_𝑧𝑣₀(0,𝑡)=0,𝑡≥0. (2.15) The solution to (2.12), (2.13), and (2.14) is [5, Equation 2.4.13],𝑣₀(𝑧,𝑡)=𝑇₀𝑧erf2√𝑘𝑡+𝑘_geo𝑧, (2.16) and the solution to (2.12), (2.13), and (2.15) is [5, Equation 2.9.7] 𝑣₀(𝑧,𝑡)=𝑇₀+𝑘_geo𝑧𝑧erf2√+2√𝑘𝑡𝑘𝑡√𝜋−𝑧exp²,4𝑘𝑡 (2.17) where the error function is defined as (𝑧)2erf∶=√𝜋^𝑧₀𝑒^−𝑢2𝑑𝑢. (2.18)

2.3. The Problem for 𝑣_𝑠

The partial differential equation that 𝑣_𝑠(𝑧,𝑡) satisfies is𝜕𝑣_𝑠𝜕𝜕𝑡=𝑘²𝑣_𝑠𝜕𝑧², (2.19) subject to an homogeneous initial condition𝑣_𝑠(𝑧,0)=0,𝑧≥0, (2.20) and to the time-dependent Dirichlet boundary condition 𝑣_𝑠(0,𝑡)=𝜓(𝑡),𝑡≥0, (2.21) or to the time-dependent Neumann boundary condition −𝜆𝜕_𝑧𝑣_𝑠(0,𝑡)=𝜓(𝑡),𝑡≥0. (2.22) The solution to (2.19), (2.20) and (2.21) is [5, Equation 2.5.1.] 𝑣_𝑠()=2𝑧,𝑡√𝜋^∞√_{𝑧/2𝑘𝑡}𝜓𝑧𝑡−²4𝑘𝑢²𝑒^−𝑢2𝑑𝑢, (2.23) and the solution to (2.19), (2.20) and (2.22) is [5, Equation 2.9.9] 𝑣_𝑠()=√𝑧,𝑡𝑘𝜆√𝜋^𝑡₀𝜓(𝑡−𝑢)√𝑢𝑒^{−𝑧2/4𝑘𝑢}𝑑𝑢. (2.24)

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, 𝜓(𝑡) is a periodic function of a certain period 𝜏, continuous and bounded, that we may expand in Fourier series𝜓(𝑡)=𝑎₀2+^∞_𝑛=1𝑎_𝑛cos(𝑛𝜔𝑡)+𝑏_𝑛sin(𝑛𝜔𝑡), (3.1) where the frequency is 𝜔=2𝜋/𝜏, and the coefficients 𝑎_𝑛 and 𝑏_𝑛 are given by 𝑎₀=2𝜏^𝜏/2_−𝜏/2𝜓(𝑢)𝑎𝑑𝑢,_𝑛=2𝜏^𝜏/2_−𝜏/2𝜓(𝑢)()𝑏cos𝑛𝜔𝑢𝑑𝑢,_𝑛=2𝜏^𝜏/2_−𝜏/2𝜓(𝑢)sin(𝑛𝜔𝑢)𝑑𝑢, (3.2) where 𝑛∈ℕ.

3.1. Dirichlet Boundary Condition

Substituting (3.1) in (2.23), we obtain 𝑣_𝑠()=𝑎𝑧,𝑡⁰2𝑧erfc2√+2𝑘𝑡√𝜋^∞_𝑛=1𝑎_𝑛^∞_{√𝑧/2𝑘𝑡}cos𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2𝑑𝑢+𝑏_𝑛^∞_{√𝑧/2𝑘𝑡}sin𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2𝑑𝑢, (3.3) where erfc(𝑧)=1−erf(𝑧) is the complementary error function. The integrals given in (3.3) for 𝑡≥0 may be rewritten as ^∞_{√𝑧/2𝑘𝑡}cos𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2𝑑𝑢=Re𝑒^{𝑖𝑛𝜔𝑡}𝐼_𝐷,𝑛(𝑧,𝑡)_{∞√𝑧/2𝑘𝑡},sin𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2𝑑𝑢=Im𝑒^{𝑖𝑛𝜔𝑡}𝐼_𝐷,𝑛(𝑧,𝑡), (3.4) where we have defined𝐼_𝐷,𝑛(𝑧,𝑡)∶=^∞_{√𝑧/2𝑘𝑡}exp−𝑖𝑛𝜔𝑧²4𝑘𝑢²−𝑢²𝑑𝑢. (3.5) In order to calculate (3.5), let us apply the integral [6, Equation 7.4.33] exp−𝑎𝑥²−𝑏𝑥²√𝑑𝑥=𝜋4√𝑎𝑒^2√𝑎𝑏erf√√𝑎𝑥+𝑏𝑥+𝑒^{√−2𝑎𝑏}erf√√𝑎𝑥−𝑏𝑥+const, (3.6) taking 𝑎=1, 𝑏=𝑖𝑛𝜔𝑧²/4𝑘, and 𝑧≥0, so that𝐼_𝐷,𝑛()=√𝑧,𝑡𝜋4𝑧√exp𝑖𝑛𝜔√𝑘𝑧erfc2√+√𝑘𝑡√𝑖𝑛𝜔𝑡+exp−𝑧𝑖𝑛𝜔√𝑘𝑧erfc2√−√𝑘𝑡𝑖𝑛𝜔𝑡. (3.7) Finally, substituting (3.7) in (3.4), we obtain for 𝑧,𝑡≥0,^∞_{√𝑧/2𝑘𝑡}cos𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2√𝑑𝑢=𝜋4Re𝑒^{𝑖𝑛𝜔𝑡}𝐴⁺_𝑛(𝑧,𝑡)+𝐴⁻_𝑛(𝑧,𝑡)_{∞√𝑧/2𝑘𝑡},sin𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2√𝑑𝑢=𝜋4Im𝑒^{𝑖𝑛𝜔𝑡}𝐴⁺_𝑛(𝑧,𝑡)+𝐴⁻_𝑛(𝑧,𝑡), (3.8) where we have defined𝐴^±_𝑛()±𝑧√𝑧,𝑡∶=exp𝑖𝑛𝜔√𝑘𝑧erfc2√±√𝑘𝑡𝑖𝑛𝜔𝑡. (3.9) 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 𝑣_𝑠()=𝑎𝑧,𝑡⁰√𝑘2𝜆√𝜋^𝑡₀𝑒^{−𝑧2/4𝑘𝑢}√𝑢+√𝑑𝑢𝑘𝜆√𝜋^∞_𝑛=1𝑎_𝑛^𝑡₀cos[𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}𝑑𝑢+𝑏_𝑛^𝑡₀sin[𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}𝑑𝑢. (3.10)(3.11)(3.12) Integrating by parts and performing the substitution 𝜁=𝑧/2√𝑘𝑢, considering 𝑧,𝑡≥0, we may calculate the integral given in (3.10) ^𝑡₀𝑒^{−𝑧2/4𝑘𝑢}√𝑢𝑑𝑢=2√𝑡exp−𝑧²−√4𝑘𝑡𝜋𝑧√𝑘𝑧erfc2√.𝑘𝑡 (3.13) In order to calculate the integrals given in (3.11) and (3.12), let us perform the substitution 𝑢=𝑣², considering 𝑡≥0, so thatcos[^𝑡₀𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}𝑑𝑢=2^√𝑡₀cos𝑛𝜔𝑡−𝑣²𝑒^{−𝑧2/4𝑘𝑣2}𝑒𝑑𝑣=2Re^{𝑖𝑛𝜔𝑡}𝐼_𝑁,𝑛(𝑧,𝑡), (3.14) and similarlysin[^𝑡₀𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}𝑑𝑢=2Im𝑒^{𝑖𝑛𝜔𝑡}𝐼_𝑁,𝑛(𝑧,𝑡), (3.15) where we have defined𝐼_𝑁,𝑛(𝑧,𝑡)∶=^√𝑡₀exp−𝑖𝑛𝜔𝑣²−𝑧²4𝑘𝑣²𝑑𝑣. (3.16) Taking 𝑎=𝑖𝑛𝜔, 𝑏=𝑧²/4𝑘, and 𝑧≥0 in (3.6), we may rewrite (3.16) as𝐼_𝑁,𝑛()=√𝑧,𝑡𝜋4√√𝑖𝑛𝜔exp−𝑧𝑖𝑛𝜔√𝑘𝑧erfc2√−√𝑘𝑡𝑧√𝑖𝑛𝜔𝑡−exp𝑖𝑛𝜔√𝑘𝑧erfc2√+√𝑘𝑡𝑖𝑛𝜔𝑡. (3.17) Finally, substituting (3.17) in (3.14) and (3.15), for 𝑧,𝑡≥0, we obtain cos[^𝑡₀𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}√𝑑𝑢=𝜋2𝑒Re^{𝑖𝑛𝜔𝑡}√𝐴𝑖𝑛𝜔^−𝑛(𝑧,𝑡)−𝐴⁺_𝑛(𝑧,𝑡)_𝑡0,sin[𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}√𝑑𝑢=𝜋2𝑒Im^{𝑖𝑛𝜔𝑡}√𝐴𝑖𝑛𝜔^−𝑛(𝑧,𝑡)−𝐴⁺_𝑛(𝑧,𝑡). (3.18) 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 ^∞₀cos𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2√𝑑𝑢=𝜋4Re𝑒^{𝑖𝑛𝜔𝑡}𝐴_𝐷,𝑛(𝑧)_∞0,sin𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2√𝑑𝑢=𝜋4Im𝑒^{𝑖𝑛𝜔𝑡}𝐴_𝐷,𝑛(𝑧), (4.1)(4.2) where we have defined𝐴_𝐷,𝑛(𝑧)∶=lim_𝑡→∞𝐴⁺_𝑛(𝑧,𝑡)+𝐴⁻_𝑛(𝑧,𝑡). (4.3) Similarly, (3.18) becomes cos[^∞₀𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}√𝑑𝑢=𝜋2𝑒Re^{𝑖𝑛𝜔𝑡}√𝐴𝑖𝑛𝜔^𝑁,𝑛(𝑧)_∞0,sin[𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}√𝑑𝑢=𝜋2𝑒Im^{𝑖𝑛𝜔𝑡}√𝐴𝑖𝑛𝜔^𝑁,𝑛(𝑧), (4.4) where we have defined𝐴_𝑁,𝑛(𝑧)∶=lim_𝑡→∞𝐴⁻_𝑛(𝑧,𝑡)−𝐴⁺_𝑛(𝑧,𝑡). (4.5)

4.1.1. Dirichlet Case

Taking into account that erfc(∞)=0 and erfc(−∞)=2, we may calculate the limit given in (4.3), so that (4.1) becomes^∞₀cos𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2√𝑑𝑢=𝜋2Re𝑒^{𝑖𝑛𝜔𝑡}√exp−𝑧𝑖𝑛𝜔√𝑘=√𝜋2𝑒^{−𝑧𝑛𝜔/2𝑘}Re𝑒^{𝑖(𝑛𝜔𝑡−𝑧𝑛𝜔/2𝑘)}, (4.6) where we have substitute √𝑖=(1+𝑖)/√2. Thus, for 𝑧,𝑡≥0,^∞₀cos𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2√𝑑𝑢=𝜋2𝑒^{−𝑧𝑛𝜔/2𝑘}𝑧√cos𝑛𝜔𝑡−𝑛𝜔√.2𝑘 (4.7) Similarly, we may obtain^∞₀sin𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2√𝑑𝑢=𝜋2𝑒^{−𝑧𝑛𝜔/2𝑘}𝑧√sin𝑛𝜔𝑡−𝑛𝜔√.2𝑘 (4.8) Notice that (4.7) and (4.8) are sinusoidal in 𝑡, with an amplitude exponentially decreasing in depth 𝑧. In the literature, we may find [7, Equation 3.928.1-2]. ^∞₀𝑒^{−𝑝2𝑥2−𝑞2/𝑥2}sin𝑎²𝑥²+𝑏²𝑥²√𝑑𝑥=𝜋𝑒2𝑟^{−2𝑟𝑠cos(𝐴+𝐵)}×sin(𝐴+2𝑟𝑠sin(𝐴+𝐵))_∞0,𝑒^{−𝑝2𝑥2−𝑞2/𝑥2}cos𝑎²𝑥²+𝑏²𝑥²√𝑑𝑥=𝜋𝑒2𝑟^{−2𝑟𝑠cos(𝐴+𝐵)}×cos(𝐴+2𝑟𝑠sin(𝐴+𝐵)), (4.9) where 𝑎²+𝑝²>0 and 𝑟=⁴𝑎⁴+𝑝⁴,𝑠=⁴𝑏⁴+𝑞⁴,1𝐴=2tan⁻¹𝑎²𝑝²1,𝐵=2tan⁻¹𝑏²𝑞², (4.10) thus, taking 𝑝=1, 𝑞=0, and 𝑎=0,^∞₀𝑒^−𝑥2𝑏sin²𝑥²√𝑑𝑥=𝜋2𝑒^−√2𝑏sin√2𝑏_∞0,𝑒^−𝑥2𝑏cos²𝑥²√𝑑𝑥=𝜋2𝑒^−√2𝑏cos√2𝑏. (4.11) Therefore, rewriting the left side of (4.7) as^∞₀cos𝑛𝜔𝑡−𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2𝑑𝑢=cos(𝑛𝜔𝑡)^∞₀cos𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2()𝑑𝑢+sin𝑛𝜔𝑡^∞₀sin𝑛𝜔𝑧²4𝑘𝑢²𝑒^−𝑢2𝑑𝑢, (4.12) and applying (4.11), taking 𝑏=√𝑛𝜔𝑧/2√𝑘, (𝑧≥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 𝑧,𝑡≥0, arriving at cos[^∞₀𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}√𝑑𝑢=𝜋𝑒^{−𝑧𝑛𝜔/2𝑘}√𝑧√𝑛𝜔cos𝑛𝜔𝑡−𝑛𝜔√−𝜋2𝑘4_∞0,sin[𝑛𝜔(𝑡−𝑢)]√𝑢𝑒^{−𝑧2/4𝑘𝑢}√𝑑𝑢=𝜋𝑒^{−𝑧𝑛𝜔/2𝑘}√𝑧√𝑛𝜔sin𝑛𝜔𝑡−𝑛𝜔√−𝜋2𝑘4. (4.13)(4.14) Once again, (4.13) and (4.7) are sinusoidal in 𝑡, with an amplitude exponentially decreasing in depth 𝑧. In the literature, we may find [7, Equation 3.957.1-2] ^∞₀𝑥^𝜇−1𝑒^{−𝛽2/4𝑥}𝛽cos𝑎𝑥𝑑𝑥=2√𝑎^𝜇𝑒^{−𝑖𝜋𝜇/4}𝐾_𝜇𝛽𝑒^𝑖𝜋/4√𝑎+𝑒^{𝑖𝜋𝜇/4}𝐾_𝜇𝛽𝑒^{−𝑖𝜋/4}√𝑎,_∞0𝑥^𝜇−1𝑒^{−𝛽2/4𝑥}𝛽sin𝑎𝑥𝑑𝑥=𝑖2√𝑎^𝜇𝑒^{−𝑖𝜋𝜇/4}𝐾_𝜇𝛽𝑒^𝑖𝜋/4√𝑎−𝑒^{𝑖𝜋𝜇/4}𝐾_𝜇𝛽𝑒^{−𝑖𝜋/4}√𝑎, (4.15) where Re𝛽>0, Re𝜇<1 and 𝑎>0. Therefore, taking in (4.15),𝜇=1/2, 𝑎=𝑛𝜔, 𝛽=𝑧/√𝑘, (𝑧≥0), and knowing that the Macdonald function of order 1/2 is [8, Equation 5.5.5], 𝐾_1/2(𝑧)=𝜋𝑒2𝑧^−𝑧, (4.16) we eventually get the same results as (4.13) and (4.14).

4.2. Solution on the Ground Surface

Notice that in the special case 𝑧=0, (3.8) becomes trivial, cos(^∞₀𝑛𝜔𝑡)𝑒^−𝑢2√𝑑𝑢=𝜋2cos(𝑛𝜔𝑡_∞0),sin(𝑛𝜔𝑡)𝑒^−𝑢2√𝑑𝑢=𝜋2sin(𝑛𝜔𝑡), (4.17) while (3.18) are reduced to cos[^𝑡₀𝑛𝜔(𝑡−𝑢)]√𝑢𝑑𝑢=√𝑒𝜋Re^{𝑖𝑛𝜔𝑡}√√𝑖𝑛𝜔erf𝑖𝑛𝜔𝑡_𝑡0,sin[𝑛𝜔(𝑡−𝑢)]√𝑢𝑑𝑢=√𝑒𝜋Im^{𝑖𝑛𝜔𝑡}√√𝑖𝑛𝜔erf𝑖𝑛𝜔𝑡. (4.18) Since [8, Equation 2.4.1] erf√𝑖𝑥√𝑖=√2𝐶2𝜋𝑥−𝑖𝑆2𝜋𝑥,𝑥∈ℝ, (4.19) where 𝐶(𝑧)∶=^𝑧₀cos𝜋𝑡²2𝑆(𝑧)𝑑𝑡,∶=^𝑧₀sin𝜋𝑡²2𝑑𝑡 (4.20) are the Fresnel integrals; then (4.18) may be expressed as cos[^𝑡₀𝑛𝜔(𝑡−𝑢)]√𝑢𝑑𝑢=2𝜋()𝐶𝑛𝜔cos𝑛𝜔𝑡2𝑛𝜔𝑡𝜋+sin(𝑛𝜔𝑡)𝑆2𝑛𝜔𝑡𝜋_𝑡0,sin[𝑛𝜔(𝑡−𝑢)]√𝑢𝑑𝑢=2𝜋()𝐶𝑛𝜔sin𝑛𝜔𝑡2𝑛𝜔𝑡𝜋−cos(𝑛𝜔𝑡)𝑆2𝑛𝜔𝑡𝜋. (4.21) The results given in (4.21), agrees with [7, Equation 2.653.1-2] sin𝑥√𝑥𝑑𝑥=√2𝜋𝑆√𝑥,cos𝑥√𝑥𝑑𝑥=√2𝜋𝐶√𝑥. (4.22)

Acknowledgments

The authors wish to thank the financial support received from Generalitat Valenciana under Grant no. GVA 3012/2009 and from Universidad Politécnica de Valencia under Grant no. PAID-06-09. This work has been partially supported by the Structural Funds of the European Regional Development Fund (ERDF).