N NON-STATIONARY HEAT TRANSFER IN ANISOTROPIC BODIES WITH GENERAL ANISOTROPY AND GIVEN DISTRIBUTED HEAT SOURCES

A new analytical solution is obtained for the first initial-boundary problem of heat transfer in an anisotropic plate with general anisotropy (non-zero off-diagonal components of thermal conductivity tensor) and given heat sources distributed over space and time. The tensorial character of heat transfer results in presence of mixed derivatives in the differential equation of thermal conductivity. The solution is obtained by constructing a source function and further use of the basic integral Green’s function. Convergence of the functional series of improper integrals is investigated. We analyze numerical results on the effect of principal components and orientation angles of principal axes of thermal conductivity tensor on nonstationary temperature fields. AMS Subject Classification: 80A20, 74E10, 74F05


Introduction
The mathematical simulation and methods of solving the problems of theory of thermal conduction in anisotropic bodies differ considerably from similar problems in isotropic media, first of all, by tensorial character of heat transfer, which results in presence of mixed derivatives in differential equations.This essentially complicates obtaining analytical as well as numerical solutions.
As a result, it becomes impossible to apply the main method of mathematical physics, namely separation of variables because mixed derivatives with respect to spatial variables cannot be separated.To find analytical solutions of such problems, one could apply the integral methods: Fourier and Laplace transforms or Green's functions method.However, their application requires that the calculating regions are unlimited (at least on one side).
Among the first works on analytical theory of thermal conduction in anisotropic media were those by Padovan [1] and Chang and Tsou [2].Some monographs [3,4,5,6] and journal articles [7,8,9,10,11] in this area have been published by the authors of this work in recent years.The problems of anisotropic thermal conduction in semi-infinite anisotropic regions with moving boundary have been solved analytically in [12,13].
In the present work we found and investigated a new analytical solution of the first initial-boundary value problem of thermal conduction in an anisotropic strip with general anisotropy and given distributed heat sources.The analytical solution was obtained by combined application of the methods of Green's functions, integral Fourier transform and separation of variables.Based on the obtained solution, we analyzed temperature fields in the anisotropic strip with different values of principal components and orientation angles of principal axes of thermal conductivity tensor.

Formulation of the Problem
We consider the following problem of determination of non-stationary temperature field u (x, y, t) in an anisotropic strip of thickness (Fig. 1) at whose boundaries the temperature distributions are µ 1 (x, t) at y = 0 and µ 2 (x, t) at y = s, with distributed heat sources f (x, y, t): where η (z) is the Heaviside unit function (η (z) = 0 at z < 0 and η(z) = 1 otherwise); x, y, t are the spatial coordinates and time, respectively; c, ρ are the heat capacity and density; λ xx , λ xy , λ yy are the components of thermal conductivity tensor in the Cartesian coordinate system; they are determined as [3] λ ξ , λ η are the principal components of thermal conductivity tensor; ϕ is the angle orienting principal axes Oξ, Oη of the thermal conductivity tensor relative to the axes Ox, Oy of the Cartesian coordinate system.
The ratio between the maximal main coefficient of thermal conductivity and the minimal main coefficient λ = λ ξ /λ η (or λ = κ η /λ ξ is called degree of anisotropy.It may take on values from unity (isotropic materials) to 200 [11].
The thermal conductivity tensor with components given in Eq. ( 5) is of arbitrary type at λ xy = λ yx = 0.A body with ϕ = 0 or ϕ = π/2 is orthotropic; it has a zero component λ xy at mixed derivative.In this case, λ xx = λ yy .
To apply the integral methods for solving the formulated problem Eqs.(1) (5), it is necessary to specify the conditions at infinity:

Method of Solution
To solve the above problem Eqs.(1)( 6), we use a source function method.For this it is necessary to construct a Green's function G (x, y, t, ψ, ζ) for a source of unit power placed at a point with coordinates ψ, ζ.According to the general theory of the method of Green's function, the latter satisfies a homogeneous ordinary differential equation and homogeneous boundary conditions that correspond to the appropriate Eqs.(1)(3) and inhomogeneous initial condition with a non-uniformity as Dirac delta function.That is the function G (x, y, ψ, ζ) satisfies the problem where δ (x − ψ) and δ (y − ζ) are delta functions.
Since the anisotropic plate is not limited in direction of variable x, it is possible to apply the integral Fourier transform when solving the problem Eqs.(7) (10).For this absolute integrability of the source function (i.e.existence of the integral ) is necessary.Absolute integrability of all inhomogeneous functions in the above problem is sufficient for existence of this integral.In this case, the improper integral within the infinite limits of delta function equals unity, thus satisfying the requirement of absolute integrability.
By applying the direct Fourier transform to the problem Eqs. ( 7) (10), we arrive at the following problem concerning the transform When calculating the Fourier transforms of x derivatives of the first and second orders in Eq. ( 7), we used the additional conditions from Eq. ( 6).Contrary to Fourier transform of the heat conduction equation for isotropic region, an imaginary constant appears before the first derivative in the case of anisotropic region.
To simplify the problems Eqs.(11) (14), let us substitute and determine the coefficients A and B in such a way that the coefficients at G ω and ∂G ω /∂y in Eq. ( 11) are zero; we get where As a result, we arrive at a one-dimensional in the spatial variable y nonstationary problem on function G 00 (y, t, ψ, ζ) determined by the equality with homogeneous boundary conditions and inhomogeneous initial condition: Now it is possible to apply the well-grounded variable separation method [12] to the problem Eqs. ( 18)(21) as a result, we find the solution as series of eigenfunctions Y k = C k • sin πk s y : By applying inverse Fourier transform to Eq. ( 24), we get solution of the problem Eqs.(7) (10) for source function G (x, y, t, ψ, ζ): Integration with respect to variable ω in Eq. (25) gives According to the Green's basic integral formula [15] with allowance made for the Duhamel's theorem [16], the solution of completely inhomogeneous problem of anisotropic heat conduction Eqs.(1)(5) can be written as Equation ( 27) with allowance made for Eq. ( 26) is the exact solution of the problem Eq. (1)(5) concerning heat conduction in an anisotropic strip with arbitrarily given functions f (x, y, t) , µ 1 (x, t) , µ 2 (x, t) , g (x, y) .Let us consider a special case of the problem Eq. ( 1) (5).Assume that f (x, y, t) = 0, g (x, y) = 0; (28) Then, by substituting the source function Eq. ( 26) in Eq. ( 27) and calculating the corresponding integrals, we get solution of the problem Eq. ( 1), Eqs.(28)(30) and Eq. ( 5)as where If the initial condition is g (x, y) = g = const, then it is necessary to add this value to the solution Eq. (31).The latter can be used for calculating temperature fields in anisotropic bodies as well as for testing of numerical methods.

Analysis of Results
Using Eq. ( 27), we performed mass calculations of non-stationary temperature fields in anisotropic plates, at different boundary and initial conditions and for different principal components λ ξ , λ η and orientation angles of principal axes Oξ, Oη of thermal conductivity tensor.Some results of the above calculations are presented below.The initial data were as follows: Shown in Fig. 2 is the temperature field at the moment t = 50s from the initial piecewise continuous temperature distribution given by the expression g(x, y) = 500, −0.03 ≤ |x| ≤ 0.03, 0.01 ≤ y ≤ 0.05; 300, |x| > 0.03, 0 < y < 0.01, 0.05 < y < s.
The boundary conditions µ 1 and µ 2 were taken constant: µ 1 = µ 2 = 300K.The isothermal lines were ellipses, with major axis oriented along the principal axis with higher coefficient of thermal conductivity (in our case, along the axis Oξ ) and minor axis oriented along the principal axis Oη.So the principal axes Oξ, Oη of thermal conductivity tensor are the symmetry axes of temperature field.Since µ 2 = g (x, y) = const, heating spreads from the boundary y = 0 along the principal axis Oξ with higher coefficient of thermal conductivity.It is not symmetry axis of isothermal lines because temperature gradient on the right of this axis is considerably higher than that on the left of the axis Oξ.Shown in Fig. 4 are similar results in the case of exponential boundary condition at y = 0 as µ 1 (x, t) = µ 2 (x, t) = 200 • exp(−(x/0.02) 2 ) + 300, for g(x, y) = µ 2 = 300K = const, in the moment t = 50s, with orientation angle ϕ = 60 • of the principal axis Oξ.The temperature field does not differ qualitatively from that in Fig. 3, but the speed of isotherm moving is higher than that in Fig. 3. Shown in Fig. 5 is temperature distribution with the boundary conditions µ 1 (x, t) = µ 2 (x, t) = 200 • exp(−(x/0.02) 2 ) + 300, t = 50s, ϕ = 30 • and the initial condition g(x.y) = 300K i.e. in the case of heating from two boundaries.The temperature field is obviously demarcated by a linear separatrix along the principal axis with higher coefficient of thermal conductivity and a curvilinear separatrix toward the principal axis Oη with lower coefficient of thermal conductivity.These separatrices divide temperature fields from different boundaries and they are fixed for stationary boundary conditions.
A body is referred to as orthotropic if the principal axes of thermal conductivity tensor coincide with the axes of Cartesian coordinate system Ox ≡ Oξ, Oy ≡ Oη, with ϕ = 0 or Oy ≡ Oξ, Ox ≡ Oη, with ϕ = 90 • .Shown in Fig. 6 are the results of temperature field calculation for orthotropic bodies in the case µ 1 (x, t) = µ 2 (x, t) = 200 • exp(−(x/0.02) 2 ) + 300, g(x, y) = 300K, t = 50s, at ϕ = 0 • (Fig. 6a) and ϕ = 90 • (Fig. 6b).One can see that at ϕ = 0 • (the principal axis with higher coefficient of thermal conductivity is along the axis ) (Fig. 6a), the temperature field is concentrated in the neighborhood of boundaries; this case may be used to cool the central part of the plate.Contrary to this, at ϕ = 90 • (Fig. 6b), the temperature field is concentrated in the central part of the plate; this may be used to just to localize temperature field.In this case, the peripheral areas remain cold.

Figure 1 :
Figure 1: Calculating region and coordinate system.

Figure 2 :
Figure 2: Temperature field in the case g(x, y) =