An Analytic Solution for 2D Heat Conduction Problems with Space–Time-Dependent Dirichlet Boundary Conditions and Heat Sources

Abstract

This study proposes a closed-form solution for the two-dimensional (2D) transient heat conduction in a rectangular cross-section of an infinite bar with space–time-dependent Dirichlet boundary conditions and heat sources. The main purpose of this study is to eliminate the limitations of the previous study and add heat sources to the heat conduction system. The restriction of the previous study is that the values of the boundary conditions and initial conditions at the four corners of the rectangular region should be zero. First, the boundary value problem of 2D heat conduction system is transformed into a dimensionless form. Second, the dimensionless temperature function is transformed so that the temperatures at the four endpoints of the boundary of the rectangular region become zero. Dividing the system into two one-dimensional (1D) subsystems and solving them by combining the proposed shifting function method with the eigenfunction expansion theorem, the complete solution in series form is obtained through the superposition of the subsystem solutions. Three examples are studied to illustrate the efficiency and reliability of the method. For convenience, the space–time-dependent functions used in the examples are considered separable in the space–time domain. The linear, parabolic, and sine functions are chosen as the space-dependent functions, and the sine, cosine, and exponential functions are chosen as the time-dependent functions. The solutions in the literature are used to verify the correctness of the solutions derived using the proposed method, and the results are completely consistent. The parameter influence of the time-dependent function of the boundary conditions and heat sources on the temperature variation is also investigated. The time-dependent function includes exponential type and harmonic type. For the exponential time-dependent function, a smaller decay constant of the time-dependent function leads to a greater temperature drop. For the harmonic time-dependent function, a higher frequency of the time-dependent function leads to a more frequent fluctuation of the temperature change.

1. Introduction

Heat conduction problems commonly exist in various modern mechanical equipment and advanced instruments, and the boundary conditions and the heat sources of these heat conduction problems are often space–time-dependent. Over the years, much research has been devoted to these problems. The main methods to solve these problems include numerical methods, analytic methods, and experimental methods. The literature review below focuses on the study of heat conduction problems with time- and/or space-dependent boundary conditions and/or heat sources.

Some advanced books on heat conduction described some classical techniques for solving heat conduction problems [1,2,3]. These include Laplace transform, Green’s function, and Duhamel’s theorem. Different approximation methods have been used to study one-dimensional (1D) heat conduction problems with space- and/or time-dependent boundary conditions and heat sources, such as the iterative perturbation method used by Holy [4] and the eigenfunction expansion method used by Özisik and Murray [5]. Johansson and Lesnic [6], Chantasiriwan [7], and Young et al. [8] applied the fundamental solutions method to the time-dependent heat conduction problems. On the basis of the boundary element method, Zhu, Liu, and Lu [9] used the Laplace transform, and Bulgakov, Sarler, and Kuhn [10] used the finite difference scheme to find the solution in the time domain. Lee et al. [11] proposed an integration-free-solution method to derive the analytic closed solution to the one-dimensional heat conduction problem with time-dependent boundary conditions. Furthermore, they extended the solution method and successfully applied it to the inverse analysis of heat conduction problems [12,13].

Several numerical techniques have been proposed for two- and three-dimensional (2D and 3D) heat conduction problems with space- and/or time-dependent boundary conditions and/or heat sources. Walker [14] applied the fundamental solution of diffusion combined with time integration to solve the diffusion equation. Chen, Golberg, and Hon [15] applied the modified Helmholtz fundamental solution to the diffusion equation. Zhu [16] and Sutradhar, Paulino, and Gray [17] used the Laplace transform to deal with the time derivative term in the diffusion equation. Burgess and Mahajerin [18] used the fundamental collocation method to the problems of arbitrary shapes subjected to mixed time-dependent boundary conditions and arbitrary initial conditions. Young, Tsai, and Fan [19] solved the nonhomogeneous diffusion problems using the fundamental solution method and the dual reciprocity method. Numerical computation using the Padé scheme was proposed by Siddique [20] to solve two-dimensional diffusion problems. The finite integral transform method was used by Singh et al. [21] to solve the problem of asymmetric heat conduction in a multilayer annulus space with time-dependent boundary conditions and/or heat sources. Hematiyan et al. [22] presented a boundary element method to analyze the 2D and 3D uncoupled thermoelastic problems with space- and time-dependent heat sources. An enhanced state-space method considering laminate approximation theory was introduced by Daneshjou et al. [23] to analyze the non-Fourier heat conduction of an infinite 2D functionally graded hollow cylinder influenced by a time-dependent heat source.

Gu et al. [24] proposed the generalized finite difference method to recover the time-dependent heat sources associated with a 3D heat equation. The eigenfunction-based solutions of various heat conduction models exhibited a mismatch at the boundaries. Biswas et al. [25] homogenized the generalized time-dependent boundary conditions to remove this mismatch from the solution. Akbari et al. [26] employed Duhamel’s theorem to evaluate the non-Fourier heat conduction of the 3D hollow spheres under the space- and time-dependent boundary conditions. Zhou et al. [27] proposed a polygonal boundary element method to solve the nonlinear heat conduction problems with space-dependent heat source and temperature-dependent thermal properties. The authors [28] derived a closed-form solution to the 2D heat conduction problem with the general Dirichlet boundary condition using the shifting function method with the eigenfunction expansion theorem. However, the temperature values of the boundary conditions and initial condition at the four corners of the rectangular region were restricted to zero, which limited the applicability of this study. Moreover, no heat source was considered in this document.

In this paper, a closed solution to the transient heat conduction problems in a rectangular cross-section of an infinite bar with nonhomogeneous space–time-dependent Dirichlet boundary conditions and heat sources was developed. This study addresses the limitations of the previous study [28]; that is, the values of the boundary conditions and initial conditions at the four corners of the rectangular area do not need to be limited to zero and become functions of time. The space–time-dependent heat sources are also considered in this study. The scope of application becomes much larger. This solution has a wide range of applications, such as laser heating, living tissue research, electronic devices, and microstructural applications [12,25,29,30]. The proposed method can accurately solve the heat conduction problem modeled as a rectangular cross-section of an infinite bar with Dirichlet boundary conditions and heat sources. The first step in solving the two-dimensional heat conduction problem is to change the physic system to be a dimensionless form. Afterward, the temperature function is transformed so that the temperature values at the four endpoints of the boundary conditions of the rectangular region become zero. Then, the system is divided into two separate one-dimensional subsystems. The governing differential equations with nonhomogeneous space- and time-dependent boundary conditions are transformed into differential equations with homogeneous boundary conditions by applying the shifting function method which does not require any integral transformation. Under the transformed homogeneous boundary conditions, the closed-form solution of the system can be further deduced by using the eigenfunction expansion theorem. Combining the solutions of the two individual subsystems, the analytic solution of the original two-dimensional heat conduction system can be obtained. Three examples are employed to illustrate the efficiency and reliability of the method. The solution obtained using the proposed method is in full agreement with the literature. The influence of the parameters is also investigated.

The contributions of this paper are as follows:

(1)

The analytic solution to 2D heat conduction problems with the general Dirichlet boundary conditions using the shifting function method with the expansion theorem method was proposed in our previous study [28]. However, there were two restrictions, the temperature values at the four corners of the rectangular area should be zero, and the heat source was also set to zero. The greatest contribution of this work is that an analytical solution is proposed first for the 2D transient heat conduction in a rectangular cross-section of an infinite bar with the space–time-dependent Dirichlet boundary conditions and heat sources. The temperatures at the four corners of the rectangular region can be functions of time.

(2)

The correctness of the solution in this study is verified by comparing the solutions of some cases using the proposed method with those of Young et al. [8], the previous work [28], and Siddique [20]. To the best of the authors’ knowledge, the other cases in this paper have never been presented in past studies. Furthermore, the case studies show that the proposed method has good convergence to the solution using series expansion and can quickly reach the converged value. The parameter influence of the time-dependent function of the boundary conditions and heat sources on the temperature change is also studied.

2. Mathematical Modeling and Dimensionless Form of Physical System

Consider the transient heat conduction for a rectangular cross-section in an infinite bar with the space–time-dependent Dirichlet boundary conditions on its four sides and the heat generation inside the bar. The material of the bar is isotropic and time-independent, implying that the material properties are constants. Figure 1 shows the geometry, heat sources, boundary conditions, and initial conditions of a rectangular cross-section in an infinite bar. Therefore, the governing equation, boundary conditions, and initial conditions of the physical system are given as follows:

$$k\left[\frac{{\partial }^{2}T(x,y,t)}{\partial x^2}+\frac{{\partial }^{2}T(x,y,t)}{\partial y^2}\right]+g(x,y,t)=\rho c\frac{\partial T(x,y,t)}{\partial t} \mathrm{in} 0\le x\le L_x, 0\le y\le L_y, t>0,$$

(1)

$$T(0,y,t)=f_1(y,t), \mathrm{at} x=0, 0\le y\le L_y,$$

(2)

$$T(L_x,y,t)=f_2(y,t), \mathrm{at} x=L_x, 0\le y\le L_y,$$

(3)

$$T(x,0,t)=f_3(x,t), \mathrm{at} y=0, 0\le x\le L_x,$$

(4)

$$T(x,L_y,t)=f_4(x,t), \mathrm{at} y=L_y, 0\le x\le L_x,$$

(5)

$$T(x,y,0)=T_0(x,y), \mathrm{at} t=0, 0\le x\le L_x, 0\le y\le L_y,$$

(6)

where $T(x,y,t)$ is the temperature function, $g(x,y,t)$ denotes the heat source inside the rectangular cross-section, $x$ and $y$ are the two-dimensional space variables, $L_x$ and $L_y$ are the thickness of the rectangular bar in the $x$- and $y$-directions, respectively, and $t$ is the time variable. In addition, $k$, $\rho$, and $c$ represent the thermal conductivity, mass density, and specific heat, respectively. $f_i(y,t)\hspace{1em}i=1,2$ and $f_i(x,t)\hspace{1em}i=3,4$ represent the general case of space–time-dependent temperatures prescribed along the surfaces at the left and right sides and at the bottom and top sides, respectively. It is worth noting that the heat source function $g(x,y,t)$ in Equation (1) is a function of space variables, which means that there are infinitely many heat sources allowed in this function. If multiple heat source functions are required to model the heat sources, the function $g(x,y,t)$ in Equation (1) can be changed to a summation term of these functions. The matching boundary conditions to initial conditions have the following properties:

$$f_1(y, 0)=T_0(0, y), f_2(y,0)=T_0(L_x,y), f_3(x,0)=T_0(x,0), f_4(x,0)=T_0(x,L_y).$$

(7)

For ease of discussion, a dimensionless form of the 2D heat conduction system is first generated. After introducing the dimensionless parameters

$$\begin{gathered} \psi (X,Y,\tau )=\frac{T(x,y,t)}{T_r}, \overline{G}(X,Y,\tau )=\frac{L_y^2 g(x,y,t)}{k{T}_r}, \tau =\frac{\alpha t}{L_y^2}, X=\frac{x}{L_x}, Y=\frac{y}{L_y}, L_r=\frac{L_y}{L_x}, \\ {\overline{F}}_i(Y,\tau )=\frac{f_i(y,t)}{T_r}, i=1,2, {\overline{F}}_i(X,\tau )=\frac{f_i(x,t)}{T_r}, i=3,4, {\psi }_0(X,Y)=\frac{T_0(x,y)}{T_r}, \end{gathered}$$

(8)

the dimensionless form of the boundary value problem becomes

$$\left[L_r^2\frac{{\partial }^2\psi (X,Y,\tau )}{\partial X^2}+\frac{{\partial }^2\psi (X,Y,\tau )}{\partial Y^2}\right]+\overline{G}(X,Y,\tau )=\frac{\partial \psi (X,Y,\tau )}{\partial \tau }, \mathrm{in} 0<X<1,\hspace{1em}0<Y<1, \tau >0,$$

(9)

$$\psi (0,Y,\tau )={\overline{F}}_1(Y,\tau ), \psi (1,Y,\tau )={\overline{F}}_2(Y,\tau ), \psi (X,0,\tau )={\overline{F}}_3(X,\tau ), \psi (X,1,\tau )={\overline{F}}_4(X,\tau ),$$

(10)

$$\psi (X,Y,0)={\psi }_0(X,Y).$$

(11)

The parameter $\alpha =\frac{k}{\rho c}$ in Equation (8) represents the thermal diffusivity, and $T_r$ represents the reference temperature. The matching conditions of the boundary conditions and the initial conditions also become

$${\overline{F}}_1(Y, 0)={\psi }_0(0, Y), {\overline{F}}_2(Y,0)={\psi }_0(1,Y), {\overline{F}}_3(X,0)={\psi }_0(X,0), {\overline{F}}_4(X,0)={\psi }_0(X,1).$$

(12)

3. The Solution Method

The previous study [28] presented a method to derive an analytic solution to a 2D heat conduction problem with general Dirichlet boundary conditions. However, the heat generation was not considered in the model. Moreover, the values of boundary conditions and initial conditions at each endpoint of the rectangular area were limited to zero, which limited the applicability of this method. In this study, the constraint of corner zeros of the boundary conditions was removed by first variable-transforming the temperature, and the space–time-dependent heat source was added to the heat conduction system. The system was then divided into two subsystems, each of which can be solved with a 1D heat conduction problem. By properly choosing the shifting function, the second-order governing differential equation with space–time-dependent boundary conditions was transformed into a differential equation with homogeneous boundary conditions. Due to the homogeneous boundary conditions, the eigenfunction expansion theorem could be applied to solve the closed-form solution of the subsystem. Lastly, the solutions of the two subsystems were added to obtain the analytic solution for 2D heat conduction in a rectangular region with space–time Dirichlet boundary conditions and heat sources

3.1. Temperature Variable Transformation

$a(\tau )$, $b(\tau )$, $c(\tau )$, and $d(\tau )$ represent the temperature at each endpoint of the rectangular section (Figure 2a) and can be expressed as follows.

$$\begin{gathered} a(\tau )={\overline{F}}_1(0, \tau )={\overline{F}}_3(0, \tau ), b(\tau )={\overline{F}}_2(0, \tau )={\overline{F}}_3(1, \tau ), \\ c(\tau )={\overline{F}}_1(1, \tau )={\overline{F}}_4(0, \tau ), d(\tau )={\overline{F}}_2(1, \tau )={\overline{F}}_4(1, \tau ). \end{gathered}$$

(13)

To increase the method applicability of our previous study [28], the nonzero temperatures at the four endpoints can be converted to zero first (Figure 2b). Therefore, variable transformation to achieve this goal can be expressed as follows:

$$\theta (X,Y,\tau )=\psi (X,Y,\tau )-\left\{\begin{array}{l}a(\tau )+X [b(\tau )-a(\tau )]+\\ Y [c(\tau )-a(\tau )+X Y (d(\tau )-c(\tau )-b(\tau )+a(\tau ))]\end{array}\right\}.$$

(14)

Substituting Equation (14) back into Equation (9) can obtain the following nonhomogeneous differential equation:

$$\left[L_r^2\frac{{\partial }^2\theta (X,Y,\tau )}{\partial X^2}+\frac{{\partial }^2\theta (X,Y,\tau )}{\partial Y^2}\right]+G(X,Y,\tau )=\frac{\partial \theta (X,Y,\tau )}{\partial \tau }, \mathrm{in} 0<X<1,\hspace{1em}0<Y<1, \tau >0,$$

(15)

where the transformed dimensionless heat source $G(X,Y,\tau )$ is defined as

$$G(X,Y,\tau )=\overline{G}(X,Y,\tau )-\rho c\left\{\begin{array}{l}\dot{a}(\tau )+X [\dot{b}(\tau )-\dot{a}(\tau )]+\\ Y [\dot{c}(\tau )-\dot{a}(\tau )+X Y [\dot{d}(\tau )-\dot{c}(\tau )-\dot{b}(\tau )+\dot{a}(\tau )]\end{array}\right\},$$

(16)

and the dot is used to indicate the differentiation with respect to the dimensionless time $\tau$.

The dimensionless boundary conditions and initial conditions become

$$\theta (0,Y,\tau )=F_1(Y,\tau )={\overline{F}}_1(Y,\tau )-\left\{a(\tau )+Y [c(\tau )-a(\tau )]\right\}, \mathrm{at} X=0, 0\le Y\le 1,$$

(17)

$$\theta (1,Y,\tau )=F_2(Y,\tau )={\overline{F}}_2(Y,\tau )-\left\{b(\tau )+Y [d(\tau )-b(\tau )]\right\}, \mathrm{at} X=1, 0\le Y\le 1,$$

(18)

$$\theta (X,0,\tau )=F_3(X,\tau )={\overline{F}}_3(X,\tau )-\left\{a(\tau )+X [b(\tau )-a(\tau )]\right\}, \mathrm{at} Y=0, 0\le X\le 1,$$

(19)

$$\theta (X,1,\tau )=F_4(X,\tau )={\overline{F}}_4(X,\tau )-\left\{c(\tau )+X [d(\tau )-c(\tau )]\right\}, \mathrm{at} Y=1, 0\le X\le 1,$$

(20)

$$\begin{gathered} \theta (X,Y,\tau )={\theta }_0(X,Y)={\psi }_0(X,Y)-\left\{\begin{array}{l}a(0)+X [b(0)-a(0)]+\\ Y [c(0)-a(0)]+X Y [d(0)-c(0)-b(0)+a(0)]\end{array}\right\}, \\ \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1. \end{gathered}$$

(21)

In addition, considering the compatibility of the boundary conditions and initial conditions, we have

$$F_1(Y,0)={\theta }_0(0,Y), F_2(Y,0)={\theta }_0(1,Y), F_3(X,0)={\theta }_0(X,0), F_4(X,0)={\theta }_0(X,1).$$

(22)

3.2. Principle of Superposition

Due to the boundary value problem with the linear characteristic, we can divide the physical system into two subsystems $A$ and $B$ in the $X$- and $Y$-directions through the principle of superposition, as shown in Figure 3.

$\theta (X,Y,\tau )$ is spilt into two branches as follows:

$$\theta (X,Y,\tau )={\theta }_a(X,Y,\tau )+{\theta }_b(X,Y,\tau ).$$

(23)

For subsystem $A$, the governing equation, and the boundary and initial conditions of the heat conduction problem are

$$\left[L_r^2\frac{{\partial }^2{\theta }_a(X,Y,\tau )}{\partial X^2}+\frac{{\partial }^2{\theta }_a(X,Y,\tau )}{\partial Y^2}\right]+G_a(X,Y,\tau )=\frac{\partial {\theta }_a(X,Y,\tau )}{\partial \tau }, \mathrm{in} 0\le X\le 1, 0\le Y\le 1,\tau >0,$$

(24)

$${\theta }_a(0,Y,\tau )=F_1(Y,\tau ), {\theta }_a(1,Y,\tau )=F_2(Y,\tau ), {\theta }_a(X,0,\tau )=0, {\theta }_a(X,1,\tau )=0,$$

(25)

$${\theta }_a(X,Y,0)={\theta }_{a0}(X,Y), \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1.$$

(26)

Similarly, for subsystem $B$, the governing equation, and the boundary and initial conditions of the heat conduction problem are

$$\left[L_r^2\frac{{\partial }^2{\theta }_b(X,Y,\tau )}{\partial X^2}+\frac{{\partial }^2{\theta }_b(X,Y,\tau )}{\partial Y^2}\right]+G_b(X,Y,\tau )=\frac{\partial {\theta }_b(X,Y,\tau )}{\partial \tau }, \mathrm{in} 0\le X\le 1, 0\le Y\le 1,\tau >0,$$

(27)

$${\theta }_b(0,Y,\tau )=0, {\theta }_b(1,Y,\tau )=0, {\theta }_b(X,0,\tau )=F_3(X,\tau ), {\theta }_b(X,1,\tau )=F_4(X,\tau ),$$

(28)

$${\theta }_b(X,Y,0)={\theta }_0(X,Y)-{\theta }_{a0}(X,Y)={\theta }_{b0}(X,Y).$$

(29)

Considering the similarity of the two subsystems, for the sake of brevity, subsystem $A$ is solved first, while subsystem $B$ is solved in Appendix A.

3.3. Reduction to One-Dimensional Problem

Considering two homogeneous boundary conditions on opposite sides of a rectangular cross-section, namely, $Y=0$ and $Y=1$, it is reasonable to assume that the temperature ${\theta }_a(X,Y,\tau )$, the heat source $G_a(X,Y,\tau )$, and the dimensionless boundary values $F_i(Y,\tau )\hspace{1em}(i=1,2)$ defined in Equations (24) and (25) become

$${\theta }_a(X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\left[{\theta }_m(X,\tau )\sin (m\pi Y)\right],}$$

(30)

$$G_a(X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\left[G_{am}(X,\tau )\sin (m\pi Y)\right],}$$

(31)

$$F_i(Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\left[{\overline{F}}_{i,m}(\tau )\sin (m\pi Y)\right]}, i=1,2,$$

(32)

where $G_{am}(X,\tau )$ and ${\overline{F}}_{i,m}(\tau )\hspace{1em}(i=1,2),$ are defined as

$$G_{am}(X,\tau )=2{\displaystyle {\int }_0^1\left[G_a(X,Y,\tau )\sin (m\pi Y)\right]} dY,$$

(33)

$${\overline{F}}_{i,m}(\tau )=2{\displaystyle {\int }_0^1{F}_i}(Y,\tau )\sin (m\pi Y)dY, i=1,2.$$

(34)

${\theta }_m(X,\tau )$ in Equation (30) is determined by meeting the governing equation and the boundary conditions on both edges $X=0$ and $X=1$ (the first two terms in Equation (25)). Substituting Equations (30)–(32) into Equations (24)–(26) yields the following result:

$$\frac{\partial {\theta }_m\left(X,\tau \right)}{\partial \tau }-L_r^2\frac{{\partial }^2{\theta }_m\left(X,\tau \right)}{\partial X^2}+m^2{\pi }^2{\theta }_m\left(X,\tau \right)=G_{am}(X,\tau ), \mathrm{in} 0<X<1, \tau >0,$$

(35)

$${\theta }_m(0,\tau )={\overline{F}}_{1,m}(\tau ), {\theta }_m(1,\tau )={\overline{F}}_{2,m}(\tau ),$$

(36)

$${\theta }_m(X,0)=2{\displaystyle {\int }_0^1{\theta }_{a0}(X,Y)}\sin (m\pi Y)dY.$$

(37)

3.4. The Shifting Function Method

3.4.1. Change of Variable

To find the solution of the second-order partial differential Equation (35) with nonhomogeneous boundary conditions (Equation (36)), the following transformation equation can be employed to extend the shifting function method [11,12,13]:

$${\theta }_m(X,\tau )={\overline{\theta }}_m(X,\tau )+{\displaystyle \sum _{i=1}^2[g_{i,m}(X){\overline{F}}_{i,m}(\tau )}],$$

(38)

where ${\overline{\theta }}_m(X,\tau )$ is the transformed function, and $g_{i,m}(X)(i=1,2)$ are the two shifting functions to be determined.

Substituting Equation (38) into Equations (35)–(37) yields

$$\begin{array}{l}{\dot{\overline{\theta }}}_m(X,\tau )+{\displaystyle \sum _{i=1}^2[g_{i,m}(X){\dot{\overline{F}}}_{i,m}(\tau )}]-L_r^2\left\{{{\overline{\theta }}^{″}}_m(X,\tau )+{\displaystyle \sum _{i=1}^2[{g^{″}}_{i,m}(X){\overline{F}}_{i,m}(\tau )}]\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+m^2{\pi }^2\left\{{\overline{\theta }}_m(X,\tau )+{\displaystyle \sum _{i=1}^2[g_{i,m}(X){\overline{F}}_{i,m}(\tau )}]\right\}=G_{am}(X,\tau )\end{array},$$

(39)

where the double primes are used to denote twice differentiation with respect to dimensionless coordinate $X$. The associated boundary conditions become

$${\overline{\theta }}_m(0,\tau )+{\displaystyle \sum _{i=1}^2[g_{i,m}(0){\overline{F}}_{i,m}(\tau )]}={\overline{F}}_{1,m}(\tau ), {\overline{\theta }}_m(1,\tau )+{\displaystyle \sum _{i=1}^2[g_{i,m}(1){\overline{F}}_{i,m}(\tau )}]={\overline{F}}_{2,m}(\tau ).$$

(40)

3.4.2. The Shifting Functions

For the convenience of subsequent analysis, the shifting functions can be chosen to satisfy the following differential equation and boundary conditions:

$${g^{″}}_{i,m}(X)=0, i=1,2, 0<X<1,$$

(41)

$$g_{i,m}(0)={\delta }_{i1}, g_{i,m}(1)={\delta }_{i2},$$

(42)

where ${\delta }_{ij}$ is the Kronecker delta. Two shifting functions are, thus, easily specified as

$$g_{1,m}(X)=1-X, g_{2,m}(X)=X.$$

(43)

Substituting Equations (41)–(43) into Equations (39) and (40), the differential equation for ${\overline{\theta }}_m(X,\tau )$ is simplified as follows:

$${\dot{\overline{\theta }}}_m(X,\tau )-L_r^2{{\overline{\theta }}^{″}}_m(X,\tau )+m^2{\pi }^2{\overline{\theta }}_m(X,\tau )={\overline{G}}_m(X,\tau ),$$

(44)

and the homogeneous boundary conditions are

$${\overline{\theta }}_m(0,\tau )=0, {\overline{\theta }}_m(1,\tau )=0.$$

(45)

${\overline{G}}_m(X,\tau )$ in Equation (44) is defined as

$${\overline{G}}_m(X,\tau )=-{\displaystyle \sum _{i=1}^2\left\{g_{i,m}(X)[{\dot{\overline{F}}}_{i,m}(\tau )+m^2{\pi }^2{\overline{F}}_{i,m}(\tau )]\right\}}+G_{am}(X,\tau ).$$

(46)

In addition, the initial condition is transformed to be

$${\overline{\theta }}_m(X,0)=2{\displaystyle {\int }_0^1{\theta }_{a0}(X,Y)}\sin (m\pi Y)dY-{\displaystyle \sum _{i=1}^2[g_{i,m}(X){\overline{F}}_{i,m}(0)}].$$

(47)

3.4.3. The Series Expansion Theorem

The solution ${\overline{\theta }}_m(X,\tau )$ of Equations (44) and (45) can be solved by applying the method of separation variable. In this method, the solution ${\overline{\theta }}_m(X,\tau )$ is expressed as

$${\overline{\theta }}_m(X,\tau )={\displaystyle \sum _{n=1}^{\infty }[{\overline{\theta }}_{mn}(X) T_{mna}}(\tau )],$$

(48)

where the space variable ${\overline{\theta }}_{mn}(X)$ is solved by the following Sturm–Liouville eigenvalue problem:

$${{\overline{\theta }}^{″}}_{mn}(X)+{\omega }_n^2{\overline{\theta }}_{mn}(X)=0, 0<X<1,$$

(49)

$${\overline{\theta }}_{mn}(0)=0, {\overline{\theta }}_{mn}(1)=0.$$

(50)

The eigenfunctions ${\overline{\theta }}_{mn}(X)\hspace{1em}(n=1, 2, 3,\cdots ).$ and the corresponding eigenvalues are solved as

$${\overline{\theta }}_{mn}(X)=\sin {\omega }_{n}X, {\omega }_n=n\pi , n=1, 2, 3, \cdots$$

(51)

The eigenfunctions constitute an orthogonal set in the interval $[0, 1]$,

$${\displaystyle {\int }_0^1{\overline{\theta }}_{mi}(X) }{\overline{\theta }}_{mj}(X)dX=\left\{\begin{array}{c}0\hspace{1em}\hspace{1em}\mathrm{for}i\ne j\\ \frac{1}{2}\hspace{1em}\hspace{1em}\mathrm{for}i=j\end{array}\right..$$

(52)

Substituting Equation (51) into Equation (48), substituting Equation (48) into Equation (44), multiplying Equation (44) by ${\overline{\theta }}_{mn}(X)$, and integrating from 0 to 1, the following differential equation can be obtained:

$${\dot{T}}_{mna}(\tau )+{\lambda }_{mna}^2{T}_{mna}(\tau )={\gamma }_{mna}(\tau ),$$

(53)

where ${\lambda }_{mna}$ and ${\gamma }_{mna}$ are given as

$${\lambda }_{mna}=\sqrt{m^2+n^2{L}_r^2}\pi ,$$

(54)

$$\begin{array}{l}{\gamma }_{mna}(\tau )=2{\displaystyle {\int }_0^1{\overline{\theta }}_{mn}(X)} {\overline{G}}_m(X,\tau )dX=\frac{-2}{n\pi }\left[{\dot{\overline{F}}}_{1,m}(\tau )-{(-1)}^n{\dot{\overline{F}}}_{2,m}(\tau )\right]\\ \hspace{1em}\hspace{1em}\hspace{1em} -\frac{2m^2\pi }{n}\left[{\overline{F}}_{1,m}(\tau )-{(-1)}^n{\overline{F}}_{2,m}(\tau )\right]+2{\displaystyle {\int }_0^1{\overline{\theta }}_{mn}(X)} G_{am}(X,\tau )dX\end{array}.$$

(55)

$T_{mna}(0)$ is determined from the initial condition of the transformed function defined in Equation (47) as

$$T_{mna}(0)=4{\displaystyle {\int }_0^1\sin (n\pi X)}{\displaystyle {\int }_0^1{\theta }_{a0}(X,Y)}\sin (m\pi Y)dYdX-\frac{2}{n\pi }\left[{\overline{F}}_{1,m}(0)-{(-1)}^n{\overline{F}}_{2,m}(0)\right].$$

(56)

Equations (53) and (56) can be solved, and their general solution is as follows:

$$T_{mna}(\tau )=e^{-{\lambda }_{mna}^2\tau }{T}_{mna}(0)+{\displaystyle {\int }_0^{\tau }{e}^{-{\lambda }_{mna}^2(\tau -\phi )}{\gamma }_{mna}}(\phi )d\phi .$$

(57)

3.4.4. The Analytic Solution

After substituting the solutions for the transformation function (Equation (48)) and shifting function (Equation (43)) back into Equations (38) and (30), the closed-form solution ${\theta }_a(X,Y,\tau )$ for the subsystem $A$ is derived as follows:

$${\theta }_a(X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (n\pi X)T_{mna}(\tau )]+(1-X){\overline{F}}_{1,m}(\tau )+X{\overline{F}}_{2,m}(\tau )}}\}\sin (m\pi Y).$$

(58)

Thanks to the high symmetry with the subsystem $A$, the solution form of the subsystem $B$ can be easily obtained through a similar derivation process (see Appendix A for details):

$${\theta }_b(X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (n\pi Y)T_{mnb}(\tau )]+(1-Y){\overline{F}}_{3,m}(\tau )+Y{\overline{F}}_{4,m}(\tau )}}\}\sin (m\pi X).$$

(59)

The summation of the two subsystem solutions leads to an analytic solution for the heat conduction in the rectangular cross-section with space–time-dependent Dirichlet boundary conditions and heat sources as follows:

$$\begin{array}{l}\theta (X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (n\pi X)T_{mna}(\tau )]}}+(1-X){\overline{F}}_{1,m}(\tau )+X{\overline{F}}_{2,m}(\tau )\}\sin (m\pi Y)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (n\pi Y)T_{mnb}(\tau )]}}+(1-Y){\overline{F}}_{3,m}(\tau )+Y{\overline{F}}_{4,m}(\tau )\}\sin (m\pi X)\end{array}.$$

(60)

Using the relationship shown in Equation (14), the dimensionless temperature $\psi (X,Y,\tau )$ before temperature variable transformation can be calculated. In addition, using the first identity of Equation (8), the exact solution $T(x,y,t)$ with dimension can also be obtained.

3.4.5. The Extreme Case Study

If the heat source value $g(x,y,t)$ is zero, and the temperatures at four corners of the rectangular region $a(\tau )$, $b(\tau )$, $c(\tau )$, and $d(\tau )$, are all zeros, then the parameters $\overline{G}(X,Y,\tau )$, $G(X,Y,\tau )$, $G_a(X,Y,\tau )$, $G_b(X,Y,\tau )$, $G_{am}(X,\tau )$, and $G_{bm}(Y,\tau )$ are all zeros. Therefore, the obtained exact solution $\theta (X,Y,\tau )$ will be the same as that obtained in our previous work [28]. This extreme case is studied in Example 2.

4. Examples and Verification

In order to verify the advantages of the proposed solution method, two cases involving the presence and absence of a heat source in the medium are explored in detail below. For simplicity, we take $L_x=L_y=L_r=1$ in all examples of the 2D heat conduction problem.

4.1. With Zero Heat Source

In our previous study [28], the boundary condition was restricted such that the temperature values at the four corners of the rectangular area were zero. The method of this study is not limited to this restriction. In Example 1, we show how to solve it using the proposed method if there are nonzero temperature values at the four corners of the rectangular region.

Example 1. 

Take the heat source in the bar as  $g(x,y,t)=0$, and assume the space-time-dependent Dirichlet boundary and initial conditions as follows:

$$T(0,y,t)=f_1(y,t)=[\sin (\frac{\pi y}{2})+\cos (\frac{\pi y}{2})+1] \eta (\alpha t), \mathrm{at} x=0, 0\le y\le 1,$$

(61)

$$T(1,y,t)=f_2(y,t)=[\sin (\frac{\pi y}{2})+\cos (\frac{\pi y}{2})+1] \eta (\alpha t), \mathrm{at} x=1, 0\le y\le 1,$$

(62)

$$T(x,0,t)=f_3(y,t)=[\sin (\frac{\pi x}{2})+\cos (\frac{\pi x}{2})+1] \eta (\alpha t), \mathrm{at} y=0, 0\le x\le 1,$$

(63)

$$T(x,1,t)=f_4(y,t)=[\sin (\frac{\pi x}{2})+\cos (\frac{\pi x}{2})+1] \eta (\alpha t), \mathrm{at} y=1, 0\le x\le 1,$$

(64)

$$T(x,y,0)=\sin (\frac{\pi x}{2})+\cos (\frac{\pi x}{2})+\sin (\frac{\pi y}{2})+\cos (\frac{\pi y}{2}), \mathrm{at} t=0, 0\le x\le 1, 0\le y\le 1.$$

(65)

The dimensionless form of boundary and initial conditions, and the dimensionless heat source $\overline{G}(X,Y,\tau )$ are written as

$$\psi (0,Y,\tau )={\overline{F}}_1(Y,\tau )=\frac{\left[\sin (\frac{\pi Y}{2})+\cos (\frac{\pi Y}{2})+1\right] \eta (\tau )}{T_r}, \mathrm{at} X=0, 0\le Y\le 1,$$

(66)

$$\psi (1,Y,\tau )={\overline{F}}_2(Y,\tau )=\frac{[\sin (\frac{\pi Y}{2})+\cos (\frac{\pi Y}{2})+1] \eta (\tau )}{T_r}, \mathrm{at} X=1, 0\le Y\le 1,$$

(67)

$$\psi (X,0,\tau )={\overline{F}}_3(Y,\tau )=\frac{[\sin (\frac{\pi X}{2})+\cos (\frac{\pi X}{2})+1] \eta (\tau )}{T_r}, \mathrm{at} Y=0, 0\le X\le 1,$$

(68)

$$\psi (X,1,\tau )={\overline{F}}_4(Y,\tau )=\frac{[\sin (\frac{\pi X}{2})+\cos (\frac{\pi X}{2})+1] \eta (\tau )}{T_r}, \mathrm{at} Y=1, 0\le X\le 1,$$

(69)

$$\psi (X,Y,0)={\psi }_0(X,Y)=\frac{\sin (\frac{\pi X}{2})+\cos (\frac{\pi X}{2})+\sin (\frac{\pi Y}{2})+\cos (\frac{\pi Y}{2})}{T_r}, \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1,$$

(70)

$$\overline{G}(X,Y,\tau )=0, \mathrm{at} 0\le X\le 1, 0\le Y\le 1.$$

(71)

Considering the matching conditions at the four corners of the rectangular section at the initial moment, we have

$$\eta (0)=1.$$

(72)

In this case, $a(\tau )$, $b(\tau )$, $c(\tau )$, $d(\tau )$, and their differentiations with respect to $\tau$ are derived as

$$a(\tau )=b(\tau )=c(\tau )=d(\tau )=\frac{2\eta (\tau )}{T_r},$$

(73)

$$\dot{a}(\tau )=\dot{b}(\tau )=\dot{c}(\tau )=\dot{d}(\tau )=\frac{2\dot{\eta }(\tau )}{T_r}.$$

(74)

Thus, the transformed temperature $\theta (X,Y,\tau )$ and heat source $G(X,Y,\tau )$ are given as

$$\theta (X,Y,\tau )=\psi (X,Y,\tau )-\frac{2\eta (\tau )}{T_r},$$

(75)

$$G(X,Y,\tau )=-\frac{2\dot{\eta }(\tau )}{T_r}.$$

(76)

It is worth noting that, after transformation, the dimensionless heat source becomes nonzero.

Following the proposed solution procedure, the boundary and initial conditions are derived as

$$F_1(Y,\tau )=\frac{\left[\sin (\frac{\pi Y}{2})+\cos (\frac{\pi Y}{2})-1\right] \eta (\tau )}{T_r}, \mathrm{at} X=0, 0\le Y\le 1,$$

(77)

$$F_2(Y,\tau )=\frac{\left[\sin (\frac{\pi Y}{2})+\cos (\frac{\pi Y}{2})-1\right] \eta (\tau )}{T_r}, \mathrm{at} X=1, 0\le Y\le 1,$$

(78)

$$F_3(X,\tau )=\frac{\left[\sin (\frac{\pi X}{2})+\cos (\frac{\pi X}{2})-1\right] \eta (\tau )}{T_r}, \mathrm{at} Y=0, 0\le X\le 1,$$

(79)

$$F_4(X,\tau )=\frac{\left[\sin (\frac{\pi X}{2})+\cos (\frac{\pi X}{2})-1\right] \eta (\tau )}{T_r}, \mathrm{at} Y=1, 0\le X\le 1,$$

(80)

$${\theta }_0(X,Y)=\frac{\sin (\frac{\pi X}{2})+\cos (\frac{\pi X}{2})+\sin (\frac{\pi Y}{2})+\cos (\frac{\pi Y}{2})-2}{T_r}, \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1.$$

(81)

Next, $G(X,Y,\tau )$ and $\theta (X,Y,\tau )$ are divided into two parts as follows:

$$G_a(X,Y,\tau )=-\frac{\dot{\eta }(\tau )}{T_r}, G_b(X,Y,\tau )=-\frac{\dot{\eta }(\tau )}{T_r},$$

(82)

$${\theta }_{a0}(X,Y)=\frac{\sin (\frac{\pi X}{2})+\cos (\frac{\pi X}{2})-1}{T_r}, {\theta }_{b0}(X,Y)=\frac{\sin (\frac{\pi Y}{2})+\cos (\frac{\pi Y}{2})-1}{T_r}.$$

(83)

The associated dimensionless quantities $G_{am}(X,\tau )$, $G_{bm}(Y,\tau )$, and ${\overline{F}}_{i,m}(\tau )(i=1,2,3,4)$ are derived as

$$G_{am}(X,\tau )=G_{bm}(X,\tau )=\frac{-2\dot{\eta }(\tau )[1-{(-1)}^m]}{T_r m \pi },$$

(84)

$${\overline{F}}_{i,m}(\tau )=\frac{2\eta (\tau )[1-{(-1)}^m]}{T_r m \pi (4m^2-1)}, i=1,2,3,4.$$

(85)

From Equations (54), (55), (A25), and (A26) one has

$${\lambda }_{mna}={\lambda }_{mnb}=\sqrt{m^2+n^2} \pi ,$$

(86)

$${\gamma }_{mna}(\tau )={\gamma }_{mnb}(\tau )=-\frac{4m^2[1-{(-1)}^m][1-{(-1)}^n]}{T_r mn {\pi }^2(4m^2-1)}[4\dot{\eta }(\tau )+{\pi }^2\eta (\tau )].$$

(87)

Likewise, $T_{mna}(0)$ and $T_{mnb}(0)$ are determined from the initial conditions of the transformed functions defined in Equations (56) and (A27) as

$$T_{mna}(0)=T_{mnb}(0)=\frac{16(m^2-n^2)[1-{(-1)}^m][1-{(-1)}^n]}{T_r mn {\pi }^2(4m^2-1)(4n^2-1)}.$$

(88)

Therefore, $T_{mna}(\tau )$ $\mathrm{a}\mathrm{n}\mathrm{d}$ $T_{mnb}(\tau )$ are determined as follows:

$$\begin{array}{l}{T}_{mna}(\tau )=T_{mnb}(\tau )=\frac{2[1-{(-1)}^m][1-{(-1)}^n]}{T_r mn {\pi }^2(4m^2-1)(4n^2-1)}\{8(m^2-n^2) e^{-(m^2+n^2) {\pi }^2}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-4m^2(4n^2-1){\displaystyle {\int }_0^{\tau } e^{-(m^2+n^2) {\pi }^2(\tau -\varphi )}}[4\dot{\eta }(\varphi )+ {\pi }^2 \eta (\varphi )] d\varphi \end{array}.$$

(89)

Considering the time-dependent term of exponential type,

$$\eta (\tau )=e^{-{\pi }^2 \tau /4},$$

(90)

and substituting Equation (90) back into Equations (85) and (89) yields

$${\overline{F}}_{i,m}(\tau )=\frac{2e^{-{\pi }^2 \tau /4}[1-{(-1)}^m]}{T_r m \pi (4m^2-1)}, i=1,2,3,4,$$

(91)

$$T_{mna}(\tau )=T_{mnb}(\tau )=\frac{16(m^2-n^2) [1-{(-1)}^m][1-{(-1)}^n]}{T_r mn {\pi }^2(4m^2-1)(4n^2-1)}{e}^{-(m^2+n^2) {\pi }^2}.$$

(92)

From Equations (60) and (14), the dimensionless solutions $\theta (X,Y,\tau )$ and $\psi (X,Y,\tau )$ become

$$\begin{array}{l}\theta (X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (m\pi X)\sin (n\pi Y)+\sin (m\pi Y)\sin (n\pi X)]T_{mna}(\tau )}}\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{m=1}^{\infty }[\sin (m\pi X)+\sin (m\pi Y)]} {\overline{F}}_{1,m}(\tau )\end{array},$$

(93)

$$\psi (X,Y,\tau )=\theta (X,Y,\tau )+\frac{2e^{-{\pi }^2 \tau /4}}{T_r}.$$

(94)

$\psi (X,Y,\tau )$ is changed back to dimensional form to obtain the exact solution $T(x,y,t)$, and its series form is expressed as follows:

$$\begin{array}{l}T(x,y,t)={\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (m\pi x)\sin (n\pi y)+\sin (m\pi y)\sin (n\pi x)]}}\cdot \\ \hspace{1em}\hspace{1em}\frac{16(m^2-n^2)[1-{(-1)}^m][1-{(-1)}^n]}{mn {\pi }^2(4m^2-1)(4n^2-1)}{e}^{-(m^2+n^2) {\pi }^2}\}\\ \hspace{1em}\hspace{1em}+{\displaystyle \sum _{m=1}^{\infty }[\sin (m\pi x)+\sin (m\pi y)]} \cdot \frac{2e^{-\alpha {\pi }^2 t /4}[1-{(-1)}^m]}{ m \pi (4m^2-1)}+2e^{-\alpha {\pi }^2 t /4}\end{array}.$$

(95)

Another analytic solution in compact form was derived by Young et al. [8] as follows:

$$T(x,y,t)=[\sin (\frac{\pi x}{2})+\cos (\frac{\pi x}{2})+\sin (\frac{\pi y}{2})+\cos (\frac{\pi y}{2})] e^{-\alpha {\pi }^2 t/4}.$$

(96)

Since the exact solution of Equation (95) is in the form of series summation, when the number of expansion terms $m$ and $n$ becomes larger, the numerical result tends to be more accurate.

Table 1. The temperature of the rectangular region at $x=y=0.5$ and at different times [$\eta (t)=e^{-{\pi }^2 t/4}$ ].

$\mathit{t}$	$\mathit{T}(\mathit{x}=0.5, \mathit{y}=0.5,\mathit{t})$
	$\mathbf{Number} \mathbf{of} \mathbf{Expansion} \mathbf{Terms} (\mathit{m}=\mathit{n})$
	1	3	5	10	20	Exact Solution [8]
0	2.849	2.825	2.830	2.829	2.828	2.828
0.1	2.226	2.207	2.211	2.210	2.210	2.210
0.2	1.739	1.724	1.728	1.727	1.727	1.727
0.4	1.062	1.053	1.055	1.054	1.054	1.054
0.6	0.648	0.643	0.644	0.644	0.644	0.644
0.8	0.396	0.392	0.393	0.393	0.393	0.393
1.0	0.242	0.240	0.240	0.240	0.240	0.240
1.2	0.148	0.146	0.146	0.146	0.146	0.146

shows the midpoint temperature of the rectangular area of the bar at different times. It can be seen from Table 1 that the convergence speed is very fast; moreover, when the number of terms $m=n$ is 5, the maximum error is less than 0.1%. At $0\le \tau \le 1.2$, the temperature using $m=n\ge 10$ expansion terms almost converges to the exact solution in the literature [8].

4.2. With Nonzero Heat Sources

Two examples with the space–time-dependent Dirichlet boundary conditions and heat sources are used to demonstrate the proposed method.

Example 2. 

We assume that the space-time-dependent heat source in the medium is given as  $g(x,y,t)=-\rho c(x+y+1)e^{-t}$, and that the space-time-dependent Dirichlet boundary and initial conditions are as follows:

$$T(0,y,t)=f_1(y,t)=[\sin (\pi y)] e^{-\alpha {\pi }^{2}t}+(y+1) e^{-t}, \mathrm{at} x=0, 0\le y\le 1,$$

(97)

$$T(1,y,t)=f_2(y,t)=[\sin (\pi y)] e^{-\alpha {\pi }^{2}t}+(y+2) e^{-t}, \mathrm{at} x=1, 0\le y\le 1,$$

(98)

$$T(x,0,t)=f_3(x,t)=[\sin (\pi x)] e^{-\alpha {\pi }^{2}t}+(x+1) e^{-t}, \mathrm{at} y=0, 0\le x\le 1,$$

(99)

$$T(x,1,t)=f_4(x,t)=[\sin (\pi x)] e^{-\alpha {\pi }^{2}t}+(x+2) e^{-t}, \mathrm{at} y=1, 0\le x\le 1,$$

(100)

$$T(x,y,0)=T_0(x,y)=\sin (\pi x)+\sin (\pi y)+x+y+1, \mathrm{at} t=0, 0\le x\le 1, 0\le y\le 1.$$

(101)

Therefore, the dimensionless form of boundary and initial conditions can be expressed as

$$\psi (0,Y,\tau )={\overline{F}}_1(Y,\tau )=\frac{\sin (\pi Y) e^{-{\pi }^2 \tau }+(Y+1) e^{- \tau /\alpha }}{T_r}, \mathrm{at} X=0, 0\le Y\le 1,$$

(102)

$$\psi (1,Y,\tau )={\overline{F}}_2(Y,\tau )=\frac{\sin (\pi Y) e^{-{\pi }^2 \tau }+(Y+2) e^{- \tau /\alpha }}{T_r}, \mathrm{at} X=1, 0\le Y\le 1,$$

(103)

$$\psi (X,0,\tau )={\overline{F}}_3(Y,\tau )=\frac{\sin (\pi X) e^{-{\pi }^2 \tau }+(X+1) e^{- \tau /\alpha }}{T_r}, \mathrm{at} Y=0, 0\le X\le 1,$$

(104)

$$\psi (X,1,\tau )={\overline{F}}_4(Y,\tau )=\frac{\sin (\pi X) e^{-{\pi }^2 \tau }+(X+2) e^{- \tau /\alpha }}{T_r}, \mathrm{at} Y=1, 0\le X\le 1,$$

(105)

$$\psi (X,Y,0)={\psi }_0(Y,\tau )=\frac{\sin (\pi X)+\sin (\pi Y)+X+Y+1}{T_r}, \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1,$$

(106)

$$\overline{G}(X,Y,\tau )=-\frac{(X+Y+1) e^{- \tau /\alpha }}{\alpha T_r}, \mathrm{at} 0\le X\le 1, 0\le Y\le 1.$$

(107)

In this case, the time-dependent functions at four corners of the rectangular cross-section and their differentiation with respect to $\tau$ are given as follows:

$$a(\tau )=\frac{ e^{- \tau /\alpha }}{T_r}, b(\tau )=\frac{ 2e^{- \tau /\alpha }}{T_r}, c(\tau )=\frac{ 2e^{- \tau /\alpha }}{T_r}, d(\tau )=\frac{ 3e^{- \tau /\alpha }}{T_r},$$

(108)

$$\dot{a}(\tau )=\frac{ -e^{- \tau /\alpha }}{\alpha T_r}, \dot{b}(\tau )=\frac{ -2e^{- \tau /\alpha }}{\alpha T_r}, \dot{c}(\tau )=\frac{ -2e^{- \tau /\alpha }}{\alpha T_r}, \dot{d}(\tau )=\frac{ -3e^{- \tau /\alpha }}{\alpha T_r}.$$

(109)

The transformed variables $\theta (X,Y,\tau )$ and $G(X,Y,\tau )$ using Equations (14) and (16) become

$$\theta (X,Y,\tau )=\psi (X,Y,\tau )-\frac{(X+Y+1) e^{- \tau /\alpha }}{T_r}, G(X,Y,\tau )=0.$$

(110)

In addition, the transformed boundary and initial conditions are given as follows:

$$F_1(Y,\tau )=F_2(Y,\tau )=\frac{\sin (\pi Y) e^{-{\pi }^2 \tau }}{T_r}, \mathrm{at} 0\le Y\le 1,$$

(111)

$$F_3(X,\tau )=F_4(X,\tau )=\frac{\sin (\pi X) e^{-{\pi }^2 \tau }}{T_r}, \mathrm{at} 0\le X\le 1,$$

(112)

$${\theta }_0(X,Y)=\sin (\pi X)+\sin (\pi Y), \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1.$$

(113)

$G(X,Y,\tau )$ and ${\theta }_0(X,Y,\tau )$ can be separated into two parts:

$$G_a(X,Y,\tau )=G_b(X,Y,\tau )=0,$$

(114)

$${\theta }_{a0}(X,Y)=\frac{\sin (\pi Y)}{T_r}, {\theta }_{b0}(X,Y)=\frac{\sin (\pi X)}{T_r}.$$

(115)

In the case, using one-term expansion ($m=n=1$) in Equation (60), the analytic solution can be obtained as

$$\begin{array}{l}\theta (X,Y,\tau )=[\sin (\pi X)T_{11a}(\tau )+(1-X){\overline{F}}_{1,1}(\tau )+X {\overline{F}}_{2,1}(\tau )]\sin (\pi Y)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+[\sin (\pi Y)T_{11b}(\tau )+(1-Y){\overline{F}}_{3,1}(\tau )+Y {\overline{F}}_{4,1}(\tau )]\sin (\pi X)\end{array},$$

(116)

where the quantities ${\overline{F}}_{i,1}(\tau )(i=1,2,3,4)$ are given as

$${\overline{F}}_{i,1}(\tau )=\frac{e^{-{\pi }^2\tau }}{T_r }, i=1,2,3,4.$$

(117)

From Equations (54), (55), (A25), and (A26), one has

$${\lambda }_{11a}={\lambda }_{11b}=\sqrt{2} \pi , {\gamma }_{11a}(\tau )={\gamma }_{11b}(\tau )=0.$$

(118)

Likewise, $T_{11a}(0)$ and $T_{11b}(0)$ are determined from the initial conditions of the transformed functions defined in Equations (56) and (A27) as

$$T_{11a}(0)=T_{11b}(0)=0.$$

(119)

Therefore, one obtains

$$T_{11a}(\tau )=T_{11b}(\tau )=0.$$

(120)

From Equation (116), the solution in dimensionless form is as follows:

$$\theta (X,Y,\tau )=\frac{\sin (\pi X)+\sin (\pi Y)}{T_r} e^{-{\pi }^2 \tau }.$$

(121)

It is worth noting that the solution $\theta (X,Y,\tau )$ of the boundary value problem is the same as the analytical solution for case 1 of Example 1 [28]. After the temperature $\theta (X,Y,\tau )$ is transformed from the dimensionless temperature, the heat source and the associated boundary and initial conditions of the heat conduction system become an extreme case (case 1 of Example 1 in our previous work [28]). This can verify the correctness of the proposed method.

Furthermore, the dimensionless temperature before the transformation of this example can be derived as

$$\psi (X,Y,\tau )=\frac{[\sin (\pi X)+\sin (\pi Y)] e^{-{\pi }^2 \tau }+(X+Y+1) e^{- \tau /\alpha }}{T_r}.$$

(122)

Finally, the analytic solution $T(x,y,t)$ of this example becomes

$$T(x,y,t)=[\sin (\pi x)+\sin (\pi y)] e^{-\alpha {\pi }^2 t}+(x+y+1) e^{- t}.$$

(123)

Example 3. 

The heat source in the medium is given as $g(x,y,t)=-\rho c (x^2+y^2+4) {\eta }_0(\alpha t)$, where ${\eta }_0(\alpha t)$ is a time-dependent function. Assume that the space-time-dependent Dirichlet boundary and initial conditions are as follows:

$$T(0,y,t)=f_1(y,t)=1+y^2 \eta (\alpha t), \mathrm{at} x=0, 0\le y\le 1,$$

(124)

$$T(1,y,t)=f_2(y,t)=1+(1+y^2) \eta (\alpha t), \mathrm{at} x=1, 0\le y\le 1,$$

(125)

$$T(x,0,t)=f_3(x,t)=1+x^2\eta (\alpha t), \mathrm{at} y=0, 0\le x\le 1,$$

(126)

$$T(x,1,t)=f_4(x,t)=1+(1+x^2) \eta (\alpha t), \mathrm{at} y=1, 0\le x\le 1,$$

(127)

$$T(x,y,0)=T_0(x,y)=1+x^2+y^2, \mathrm{at} t=0, 0\le x\le 1, 0\le y\le 1.$$

(128)

The dimensionless form of boundary and initial conditions becomes

$$\psi (0,Y,\tau )={\overline{F}}_1(Y,\tau )=\frac{1+Y^2 \eta (\tau )}{T_r}, \mathrm{at} X=0, 0\le Y\le 1,$$

(129)

$$\psi (1,Y,\tau )={\overline{F}}_2(Y,\tau )=\frac{1+(1+Y^2) \eta (\tau )}{T_r}, \mathrm{at} X=1, 0\le Y\le 1,$$

(130)

$$\psi (X,0,\tau )={\overline{F}}_3(Y,\tau )=\frac{1+X^2 \eta (\tau )}{T_r}, \mathrm{at} Y=0, 0\le X\le 1,$$

(131)

$$\psi (X,1,\tau )={\overline{F}}_4(X,\tau )=\frac{1+(1+X^2) \eta (\tau )}{T_r}, \mathrm{at} Y=1, 0\le X\le 1,$$

(132)

$$\psi (X,Y,0)={\psi }_0(X,Y)=\frac{1+X^2+Y^2}{T_r}, \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1,$$

(133)

$$\overline{G}(X,Y,\tau )=-\frac{(X^2+Y^2+4) {\eta }_0(\tau )}{\alpha T_r}, \eta (0)=1, \mathrm{at} 0\le X\le 1, 0\le Y\le 1.$$

(134)

In this case, $a(\tau )$, $b(\tau )$, $c(\tau )$, $d(\tau )$, and their differentiation with respect to $\tau$ are

$$a(\tau )=\frac{ 1}{T_r}, b(\tau )=\frac{ 1+\eta \left(\tau \right)}{T_r}, c(\tau )=\frac{ 1+\eta \left(\tau \right)}{T_r}, d(\tau )=\frac{ 1+2\eta \left(\tau \right)}{T_r},$$

(135)

$$\dot{a}(\tau )=0, \dot{b}(\tau )=\frac{ \dot{\eta }(\tau )}{T_r}, \dot{c}(\tau )=\frac{ \dot{\eta }(\tau )}{T_r}, \dot{d}(\tau )=\frac{ 2\dot{\eta }(\tau )}{T_r}.$$

(136)

The transformed $\theta (X,Y,\tau )$ and $G(X,Y,\tau )$ are given as

$$\theta (X,Y,\tau )=\psi (X,Y,\tau )-\frac{1+(X+Y) \eta (\tau )}{T_r},$$

(137)

$$G(X,Y,\tau )=-\frac{(X^2+Y^2+4) {\eta }_0(\tau )+\alpha (X+Y) \dot{\eta }(\tau )}{\alpha T_r}.$$

(138)

Following the proposed procedure, one has the space–time-dependent boundary and initial conditions as follows:

$$F_1(Y,\tau )=\frac{(Y^2-Y) \eta (\tau )}{T_r}, \mathrm{at} X=0, 0\le Y\le 1,$$

(139)

$$F_2(Y,\tau )=\frac{(Y^2-Y) \eta (\tau )}{T_r}, \mathrm{at} X=1, 0\le Y\le 1,$$

(140)

$$F_3(X,\tau )=\frac{(X^2-X) \eta (\tau )}{T_r}, \mathrm{at} Y=0, 0\le X\le 1,$$

(141)

$$F_4(X,\tau )=\frac{(X^2-X) \eta (\tau )}{T_r}, \mathrm{at} Y=1, 0\le X\le 1,$$

(142)

$${\theta }_0(X,Y)=\frac{X^2+Y^2-(X+Y)}{T_r}, \mathrm{at} \tau =0, 0\le X\le 1, 0\le Y\le 1.$$

(143)

One separates $G(X,Y,\tau )$ and $\theta (X,Y,\tau )$ into two parts as follows:

$$G_a(X,Y,\tau )=-\frac{(X^2+2) {\eta }_0(\tau )+\alpha X \dot{\eta }(\tau )}{\alpha T_r}, G_b(X,Y,\tau )=-\frac{(Y^2+2) {\eta }_0(\tau )+\alpha Y \dot{\eta }(\tau )}{\alpha T_r},$$

(144)

$${\theta }_{a0}(X,Y)=\frac{X^2-X}{T_r}, {\theta }_{b0}(X,Y)=\frac{Y^2-Y}{T_r}.$$

(145)

The associated dimensionless quantities $G_{am}(X,\tau )$, $G_{bm}(Y,\tau )$, and ${\overline{F}}_{i,m}(\tau )(i=1,2)$ are

$$\begin{gathered} G_{am}(X,\tau )=\frac{-2[(X^2+2){\eta }_0(\tau )+\alpha X \dot{\eta }(\tau )][1-{(-1)}^m]}{\alpha T_r m \pi }, \\ {G}_{bm}(Y,\tau )=\frac{-2[(Y^2+2){\eta }_0(\tau )+\alpha Y \dot{\eta }(\tau )][1-{(-1)}^m]}{\alpha T_r m \pi }, \end{gathered}$$

(146)

$${\overline{F}}_{i,m}(\tau )=\frac{-4\eta (\tau ) }{T_r m^3 \pi { }^3}[1-{(-1)}^m], i=1,2,3,4.$$

(147)

In addition, one can determine $T_{mna}(0)$ and $T_{mnb}(0)$ from the initial conditions of the transformed functions defined in Equations (56) and (A27) as

$$T_{mna}(0)=T_{mnb}(0)=\frac{-8(m^2-n^2)[1-{(-1)}^m][1-{(-1)}^n]}{T_r m^3 n^3 {\pi }^4}.$$

(148)

Likewise, from Equations (54), (55), (A25), and (A26), one has

$${\lambda }_{mna}={\lambda }_{mnb}=\sqrt{m^2+n^2} \pi ,$$

(149)

$$\begin{array}{l}{\gamma }_{mna}(\tau )={\gamma }_{mnb}(\tau )=\frac{4[1-{(-1)}^m][2-(2-m^2 {\pi }^2){(-1)}^n]}{T_r m^3 n {\pi }^4}\dot{\eta }(\tau )\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{8[1-{(-1)}^m][1-{(-1)}^n]}{T_r m n {\pi }^2}\eta (\tau )-\frac{4[1-{(-1)}^m][2(n^2 {\pi }^2-1) -(3n^2 {\pi }^2-2){(-1)}^n]}{\alpha T_r m n^3 {\pi }^4}{\eta }_0(\tau )\end{array}.$$

(150)

Therefore, one can determine $T_{mna}(\tau )$ and $T_{mnb}(\tau )$ as follows:

$$\begin{array}{l}{T}_{mna}(\tau )=T_{mnb}(\tau )=\frac{-8(m^2-n^2)[1-{(-1)}^m][1-{(-1)}^n]}{T_r m^3 n^3 {\pi }^4} e^{-(m^2+n^2) {\pi }^2}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle {\int }_0^{\tau } e^{-(m^2+n^2) {\pi }^2(\tau -\varphi )}}\{\frac{4[1-{(-1)}^m][2-(2-m^2 {\pi }^2){(-1)}^n]}{T_r m^3 n {\pi }^4}\dot{\eta }(\varphi )\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{8[1-{(-1)}^m][1-{(-1)}^n]}{T_r m n {\pi }^2}\eta (\varphi )-\frac{4[1-{(-1)}^m][2(n^2 {\pi }^2-1)-(3n^2 {\pi }^2-2){(-1)}^n]}{\alpha T_r m n^3 {\pi }^4}{\eta }_0(\varphi )\} d\varphi \end{array}.$$

(151)

(Case 1) Consider the time-dependent functions to be of exponential type as follows:

$${\eta }_0(\tau )=e^{-D_0 \tau }, \eta (\tau )=e^{-D \tau },$$

(152)

where $D_0$ and $D$ represent the decay constants for the heat source and boundary conditions, respectively. From Equations (147) and (151), one obtains

$${\overline{F}}_{i,m}(\tau )=\frac{-4e^{-D\tau } }{T_r m^3 \pi { }^3}[1-{(-1)}^m], i=1,2,3,4.$$

(153)

$$\begin{array}{l}{T}_{mna}(\tau )=T_{mnb}(\tau )=\frac{-8(m^2-n^2)[1-{(-1)}^m][1-{(-1)}^n]}{T_r m^3 n^3 {\pi }^4} e^{-(m^2+n^2) {\pi }^2}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{4D [1-{(-1)}^m][2-(2-m^2 {\pi }^2){(-1)}^n]}{T_r m^3 n {\pi }^4[(m^2+n^2) {\pi }^2-D]}[e^{-D \tau }-e^{-(m^2+n^2){\pi }^2 \tau }]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{8 [1-{(-1)}^m][1-{(-1)}^n]}{T_r m n {\pi }^2 [(m^2+n^2) {\pi }^2-D]}[e^{-D \tau }-e^{-(m^2+n^2){\pi }^2 \tau }]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{4[1-{(-1)}^m][2(n^2 {\pi }^2-1)-(3n^2 {\pi }^2-2){(-1)}^n]}{\alpha T_r m n^3 {\pi }^4 [(m^2+n^2) {\pi }^2-D_0]}[e^{-D_0 \tau }-e^{-(m^2+n^2){\pi }^2 \tau }]\hspace{1em}\hspace{1em}\end{array}.$$

(154)

Furthermore, from Equations (60) and (14), the solutions $\theta (X,Y,\tau )$ and $\psi (X,Y,\tau )$ in dimensionless form are obtained as

$$\begin{array}{l}\theta (X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (n\pi X)\sin (m\pi Y)+\sin (n\pi Y)\sin (m\pi X)] T_{mna}(\tau )}}\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{m=1}^{\infty }[\sin (m\pi Y)+\sin (m\pi X)]} {\overline{F}}_{1,m}(\tau )\end{array}.$$

(155)

$$\psi (X,Y,\tau )=\theta (X,Y,\tau )+\frac{1+(X+Y)e^{-D \tau }}{T_r}.$$

(156)

$\psi (X,Y,\tau )$ is returned to the dimensional form as

$$T(x,y,t)=T_r \theta +1+(x+y) e^{-\alpha D t}.$$

(157)

Considering the case $D_0=D=\alpha =1$, it can be verified that the above series solution is the same as that given by Siddique [20] (shown below).

$$T(x,y,t)=1+(x^2+y^2) e^{- t}.$$

(158)

$T_{mna}(\tau )$ and $T_{mnb}(\tau )$ can be re-expressed as

$$\begin{array}{l}{T}_{mna}(\tau )=T_{mnb}(\tau )=\frac{-8(m^2-n^2)[1-{(-1)}^m][1-{(-1)}^n]}{T_r m^3 n^3 {\pi }^4} e^{-(m^2+n^2) {\pi }^2}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{8(m^2-n^2) [1-{(-1)}^m][1-{(-1)}^n]}{T_r m^3 n^3 {\pi }^4[(m^2+n^2) {\pi }^2-1]}[e^{- \tau }-e^{-(m^2+n^2){\pi }^2 \tau }]\end{array}.$$

(159)

The numerical result for the temperature of the rectangular region of the bar at the middle point in different times is shown in

Table 2. The temperature of the rectangular region at $x=y=0.5$ at different times [$\eta (t)={\eta }_0(t)=e^{- t}$ ].

$\mathit{t}$	$\mathit{T}(\mathit{x}=0.5, \mathit{y}=0.5,\mathit{t})$
	$\mathbf{Number} \mathbf{of} \mathbf{Expansion} \mathbf{Terms} (\mathit{m}=\mathit{n})$
	1	3	5	10	20	Exact Solution [20]
0	1.484	1.503	1.499	1.500	1.500	1.500
0.1	1.438	1.455	1.452	1.452	1.452	1.452
0.2	1.396	1.412	1.409	1.409	1.409	1.409
0.4	1.324	1.337	1.335	1.335	1.335	1.335
0.6	1.266	1.276	1.274	1.274	1.274	1.274
0.8	1.218	1.226	1.224	1.225	1.225	1.225
1.0	1.178	1.185	1.184	1.184	1.184	1.184
1.2	1.146	1.152	1.150	1.151	1.151	1.151

. From this table, one can find that the largest error was less than 0.1% when the number of terms $m=n$ was 5. The calculated temperature converged when 10-term expansion is applied. The temperature at $0\le \tau \le 1.2$ was also verified with the results of the literature [20]. It turns out that the two results were identical when applying $m=n\ge 10$ expansion terms to obtain the proposed solution.

To gain insight into the effect of the decay constant for the boundary conditions and the heat source on the temperature in the middle of the rectangular region, two figures are established. First, Figure 4a plots the three temperature curves for the three decay constants of the boundary conditions ($D=1,1.5,2$ and $D_0=1$). All three curves decayed exponentially and converge to 1. It can be found that, for the smallest decay constant $D=1$ of the boundary conditions, the decrease rate of the temperature curve was largest as expected. This could mean that a smaller decay constant for the boundary conditions causes slower heat dissipation at the boundary; thus, the heat dissipation in the middle is faster and its temperature drops more quickly. Second, the temperature curves for the three cases of the different decay constants of the heat source are shown in Figure 4b ($D_0=1,1.5,2$ and $D=1$). It can be seen that, for the smallest decay constant $D_0=1$ of the boundary conditions, the decrease rate of the temperature curve was largest. This can be explained by the fact that a smaller decay constant of the heat source leads to faster heat dissipation and, therefore, a faster temperature drop.

(Case 2) Consider the time-dependent functions to be of periodic type as follows:

$${\eta }_0(\tau )=\cos ({\omega }_0 \tau ), \eta (\tau )=\cos (\omega \tau ),$$

(160)

where ${\omega }_0$ and $\omega$ represent the frequency for the heat source and boundary conditions, respectively. Then, one has

$$\dot{\eta }(\tau )=-\omega \sin (\omega \tau ),$$

(161)

and one can obtain ${\overline{F}}_{i,m}(\tau )$, ${\gamma }_{mna}(\tau )={\gamma }_{mnb}(\tau )$, and $T_{mna}(\tau )=T_{mnb}(\tau )$ as follows:

$${\overline{F}}_{i,m}(\tau )=\frac{-4\cos (\omega \tau ) }{T_r m^3 \pi { }^3}[1-{(-1)}^m], i=1,2,3,4,$$

(162)

$$\begin{array}{l}{\gamma }_{mna}(\tau )={\gamma }_{mnb}(\tau )=\frac{4 [1-{(-1)}^m][2-(2-m^2 {\pi }^2){(-1)}^n]}{T_r m^3 n {\pi }^4}[-\omega \sin (\omega \tau )]\\ \hspace{1em}\hspace{1em}\hspace{1em}+\frac{8[1-{(-1)}^m][1-{(-1)}^n]}{T_r m n {\pi }^2} \cos (\omega \tau )-\frac{4[1-{(-1)}^m][2(n^2 {\pi }^2-1)-(3n^2 {\pi }^2-2) {(-1)}^n]}{\alpha T_r m n^3 {\pi }^4}\cos ({\omega }_0 \tau )\end{array},$$

(163)

$$\begin{array}{l}{T}_{mna}(\tau )=T_{mnb}(\tau )=\frac{-8(m^2-n^2)[1-{(-1)}^m][1-{(-1)}^n]}{T_r m^3 n^3 {\pi }^4} e^{-(m^2+n^2) {\pi }^2}\\ +\frac{4[1-{(-1)}^m][2-(2-m^2 {\pi }^2){(-1)}^n]}{T_r m^3 n {\pi }^4}\cdot \frac{{\omega }^2 [\cos (\omega \tau )-e^{-(m^2+n^2) {\pi }^2 \tau }]-\omega (m^2+n^2) {\pi }^2\sin (\omega \tau )}{{(m^2+n^2)}^2 {\pi }^4+{\omega }^2}\\ +\frac{8[1-{(-1)}^m][1-{(-1)}^n]}{T_r m n {\pi }^2}\cdot \frac{(m^2+n^2) {\pi }^2 [\cos (\omega \tau )-e^{-(m^2+n^2) {\pi }^2\tau }]+\omega \sin (\omega \tau )}{{(m^2+n^2)}^2 {\pi }^4+{\omega }^2}\\ -\frac{4[1-{(-1)}^m][2(n^2 {\pi }^2-1)-(3n^2 {\pi }^2-2){(-1)}^n]}{\alpha T_r m n^3 {\pi }^4}\cdot \frac{(m^2+n^2) {\pi }^2 [\cos ({\omega }_0\tau )-e^{-(m^2+n^2) {\pi }^2\tau }]+{\omega }_0 \sin ({\omega }_0\tau )}{(m^2+n^2){ }^2 {\pi }^4+{\omega }_0{}^2}\end{array}.$$

(164)

Therefore, from Equation (60), the solution in dimensional form is derived as

$$\begin{array}{l}\theta (X,Y,\tau )={\displaystyle \sum _{m=1}^{\infty }\{{\displaystyle \sum _{n=1}^{\infty }[\sin (n\pi X)\sin (m\pi Y)+\sin (n\pi Y)\sin (m\pi X)] T_{mna}(\tau )}}\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{m=1}^{\infty }[\sin (m\pi Y)+\sin (m\pi X)]} {\overline{F}}_{1,m}(\tau )\end{array},$$

(165)

$$\psi (X,Y,\tau )=\theta (X,Y,\tau )+\frac{1+(X+Y)\cos (\omega \tau )}{T_r}.$$

(166)

Recovering it to the dimensional form yields

$$T(x,y,t)=T_r \theta +1+(x+y) \cos (\alpha \omega t).$$

(167)

Temperature variation in the middle of the rectangular region with different decay constants of harmonic type for the boundary conditions or the heat source was studied numerically. Figure 5a shows the temperature variation in the middle of the rectangular region for the three cases, including the frequencies $\omega =\pi ,5,7$ and ${\omega }_0=\pi$. Likewise, the results for the frequencies ${\omega }_0=\pi ,5,7$ and $\omega =\pi$ are shown in Figure 5b. These plots show that all the temperature curves oscillated up and down with the horizontal line $T=1$. When the frequency of the heat source or boundary conditions was higher, the fluctuation of temperature change was more frequent.

5. Conclusions

In this paper, the analytical solution of 2D transient heat conduction in a rectangular cross-section of an infinitely long bar with the space–time-dependent Dirichlet boundary conditions and heat sources was solved using a new analytical solution method that does not require any kind of integral transformation. Three examples were studied to illustrate the efficiency and reliability of the proposed method. Some results were verified to be the same as those in the literature [8,20,28].

The new findings of the present study are as follows:

(1)

The purpose of this study was to complete the future work of our previous study [28], i.e., to remove the limitations of the previous study and add a heat source to the heat conduction system. The restriction of the temperatures of the boundary conditions and initial conditions at the four corners of the rectangular region to zero in the previous study was successfully eliminated. The zero temperature could be replaced by a function of time.

(2)

From the examples, it was found that the convergence speed was very fast, and the maximum error was less than 0.1% when only five terms were used in the series expansion of the solution. Compared with the literature, the temperature could converge to the exact solution.

(3)

The space–time-dependent functions used for the boundary conditions and heat sources in this study were considered separable in the space–time domain. The influence of the time-dependent function of the boundary conditions and heat sources on the temperature variation was investigated. For the exponential time-dependent function, a smaller decay constant ($D_0$ and $D$) of the time-dependent function ($e^{-D_0\tau }$ and $e^{-D\tau }$) led to a greater temperature drop. The temperatures with different decay constants converged to the same value. For the harmonic time-dependent function, a higher frequency (${\omega }_0$, and $\omega$) of the time-dependent function ($\cos \left({\omega }_0\tau \right)$ and $\cos \left(\omega \tau \right)$) led to a more frequent fluctuation of the temperature change. All temperature curves oscillated above and below a horizontal line.