Research Article
An Approximate Method for Solving Nonlinear Heat Diffusion Problems in Spherical Media
Department of Engineering Mathematics and Physics, Faculty of Engineering, Zagazig University, Zagazig, Egypt
In the field of linear heat transfer, the transient heat diffusion equation is linearized by considering the thermal properties to be independent of temperature, furthermore, the boundary conditions are also taken to be linear. This class of linear transient heat diffusion has been treated in detail[1-5] with exact, approximate and purely numerical methods. However, when the thermophysical properties and/or the volumetric heat source become temperature dependent, the field equation becomes nonlinear. In addition, if the temperature level becomes high, radiation and/or change of phase may occur and as a result, the boundary conditions become nonlinear. The problems of heat diffusion with nonlinear boundary conditions appear in combustion systems[4], where in the pre-ignition heating, the particle entering a furnace and traveling toward a flame front receives heat uniformly by thermal radiation from the furnace walls and losses heat uniformly by convection to the surrounding gases.
The transition from a linear to a nonlinear model introduces additional mathematical and computational difficulties, which should be dealt with. As it appears, solving nonlinear problems analytically is complicated. There is no general theory as yet available to handle all types of nonlinear problems, because each nonlinear problem requires a different technique for its analytic treatment. On the other hand, analytical approaches, when applicable, are advantageous because they provide an understanding of the role of various system parameters affecting heat transfer and because they establish the dominant features of the problem.
The integral transform techniques provide a systematic and straight-forward approach to the solution of a certain class of heat equations. The method is particularly suitable for the solution of both homogeneous and nonhomogeneous boundary value problems of heat conduction. The method will be extended to solve the heat diffusion problem in a spherical finite region subject to nonlinear boundary conditions due to radiation exchange at the interface according to the fourth power law[6-8].
PROBLEM DESCRIPTION
A solid sphere defined in a finite interval 0≤ r ≤ R is considered. A heat transfer from the surface takes place simultaneously by convection to the surrounding gases and by radiation to the enclosure. Figure 1 illustrate a small part of solid sphere of area A and emissivity ε that maintained at temperature T and exchanges energy by convection with a gases at T∞ with a heat transfer coefficient h and by radiation with the enclosure at Tr. The heat loss per unit area of the sphere by convection, is given by:
Also, the heat loss per unit area of the sphere by radiation, is given by:
where, σ is the Stefan-Boltzmann constant {σ=5.6697x10-8 W/(m2 x K4}. The heat loss per unit area of the sphere by combined mechanism of convection and radiation, is given by:
Fig. 1: | Simultaneous convection and radiation from a small area of sold sphere |
The convection coefficient h, the surface emissivity ε and the thermophysical properties of the solid are assumed invariant. Finally, the initial temperature is assumed to be uniform throughout the solid.
GOVERNING EQUATION
The above transient heat conduction problem can be described by the following partial differential equation[3]:
(1) |
subject to the following boundary conditions:
(2) |
(2a) |
and the initial condition:
(2b) |
where, is the thermal diffusivity including thermal conductivity k, specific heat C and density ρ of the solid sphere.
NON-DIMENSIONALIZATION
For convenience, we recast the above governing equation and its auxiliary conditions into a dimensionless form. Redefining the variables as follows:
Introducing these new (dimensionless) variables into the governing equations and the auxiliary conditions to obtain the problem in the more concise form:
(3) |
subject to the following boundary conditions:
(4) |
(4a) |
and the initial condition:
(4b) |
The dimensionless temperature distribution in the solid sphere is obtained by solving the above partial differential equation subject to the displayed conditions. A new dependent variable is defined as:
(5) |
Then Eq. (3) and (4) are transformed into:
(6) |
subject to:
(7) |
(7a) |
(7b) |
THE PROPOSED METHOD
The proposed approach is introduced here by solving the problem discussed in the last section. The method proceeds by treating the nonlinearity term in the boundary condition (7a) as a source in the differential Eq. (6) and keeping other conditions unchanged. The nonlinear term f (ξ, τ, Θ) must appear only at ξ = 1 and hence we make use of the Dirac-delta function to represent the function f (ξ, τ, Θ) as a source at ξ =1. Therefore, we may consider the following equivalent problem:
(8) |
subject to:
(9) |
(9a) |
(9b) |
To show the equivalence of the original problem defined by Eq. (6) and (7) with the proposed problem defined by Eq. (8) and (9), we integrate Eq. (6) w.r.t. ξ from 0 to 1 and by the aid of its auxiliary conditions, we get:
which is identical to the integration of Eq. (8) w.r.t. ξ from 0 to 1 with the aid of its auxiliary conditions. Therefore, the solution of the transient heat-conduction problem defined by Eq. (6) and (7) is equivalence to the solution of the proposed problem defined by Eq. (8) and (9).
The exact solution of the differential equation given by Eq. (8) subject to the conditions described by Eq. (9) is obtained by employing the finite integral transform technique. In the finite integral technique, the integral transform pair needed for the solution of a given problem is developed by considering representation of an arbitrary function in terms of the eigenfunctions corresponding to the given eigenvalue problem. Obtaining the required eigenvalue problem may be accomplished by considering the homogeneous part of the nonhomogeneous field equation and then employing separation of variables to obtain the following eigenvalue problem[3]:
(10) |
subject to:
(11) |
(11a) |
The solution to this problem gives the eigenfunctions:
(12) |
The eigenvalus of this problem is determined from the following transcendental equation:
(13) |
The orthogonality relation associated with the eigenvalue problem is given by Holaman[3]:
(14) |
where the normalization integral N(λn) is given by:
(15) |
We now proceed to find the general solution by following the usual procedure associated with the finite integral transform technique. We begin by defining the transform pair as follows:
Integral Transform:
(16) |
Inversion formula:
(16a) |
Operating on Eq. (8) with we obtain:
(17) |
knowing that:
(18) |
Following the standard transformation procedures[5-7] and with the aid of Eq. (18), (17) reduces to the following first order nonhomogeneous ordinary differential equation:
(19) |
subject to a transformed initial condition given by introducing the integral transform Eq. (16) into the initial condition (9b), therefore;
(20) |
The value of Θ4 (1, τ) can be obtained from the inversion formula (16a), so f (1, τ, Θ) can be evaluated and hence (19) yields:
(21) |
which is a system of nonlinear differential equations whose solution can be obtained using an appropriate numerical integration scheme. Runge-Kutta-4th order method for a finite number (say m) of the differential equations can be performed.
For j = 0, 1, 2, ... , we have:
(22) |
subject to the initial conditions:
(23) |
The procedure of the method is as follows:
Consider, for n = 1, 2,..., m, that
(24) |
(24a) |
(24b) |
(24c) |
where, Δτ is the time step and j = 0, 1, 2,...
Assume
(25) |
Finally, the integral transforms for n = 1, 2,..., m can be obtained as:
(26) |
The calculation is started with j = 0 and the integral transforms are evaluated because the initial integral transforms are available by Eq. (25), knowing the integral transform at the end of the second time step are evaluated by setting j = 1. The procedure is repeated to calculate the integral transform φn at the subsequent time steps. Once φn (τ) is obtained, the temperature distribution can be reconstructed through the use of the inversion formula.
COMPARISON BETWEEN PROPSED AND NUMERICAL METHODS
The transient heat conduction problem subject to a nonlinear boundary condition is investigated for some cases. The first case characterized by no heat exchange by radiation at the boundary, i.e. Ω = 0, so that the problem reduces to linear case. The other case has a stronger nonlinearity when Ω = 2. The study is conducted for a wide range of values of Biot number (β), dimensionless temperature of the surrounding and dimensionless temperature of the enclosure , however, only representative results are displayed in Fig. 2-4. The results of the proposed method (dotted line) are compared to those obtained from a finite difference method (solid line) (Fig. 2-4).
In the first case, Ω = 0, all the heat exchange at the outer-surface is due to convection. Figure 2 represent the dimensionless temperature-thickness variation for different Biot numbers (0.001, 0.01, 0.1), the difference in dimensionless temperature predicted from the two solution methods is within (0.2% to 3%). In the other case, the effect of stronger nonlinearity is considered by taking Ω = 2 as in Fig. 3. A difference within (0.2% to 5%) between the two solutions was noticed.
In Fig. 4a, we assume that the heat loss by radiation is zero, therefore the heat exchange at the outer surface of the sphere is due to convection. We notice that the drop of temperature is slow and increases as Biot number grows.
In Figs. 4b and 4c, the temperature response is dominated by the radiation exchange. In Fig. 4b, we take the temperature of enclosure is small enough that implies rapid heat loss from the surface of the sphere. From this Fig. 4b, one can notice that a higher drop of the temperature at a short time. However, the temperature of surface approaches to be steady state as time progresses.
Fig. 2a: | Dimensionless Temperature-Thickness at (β = 0.001 and Ω = 0) |
Fig. 2b: | Dimensionless Temperature-Thickness at (β = 0.01 and Ω = 0) |
Fig. 2c: | Dimensionless Temperature-Thickness at (β = 0.1 and Ω = 0) |
Fig. 3a: | Dimensionless Temperature-Thickness at (β = 0.01 and Ω = 0.1) |
Fig. 3b: | Dimensionless Temperature-Thickness at (β = 0.1 and Ω = 0.1) |
Fig. 3c: | Dimensionless Temperature-Thickness at (β = 1 and Ω = 0.1) |
Fig. 4a: | Dimensionless temperature-time variaion at outer surface of sphere for different biot numbers (convection only) |
Fig. 4b: | Dimensionless temperature-time variation at outer surface of sphere for different biot numbers (convection and heat loss by radiation) |
Fig. 4c: | Dimensionless temperature-time variation at outer surface of sphere for different biot numbers (convection and heat gain by radiation) |
In Fig. 4c, the temperature of enclosure is large enough, therefore, there is heat gain to the surface of the sphere. One can see a higher increase of the temperature at a short time until the steady state is reached as time progresses. In both Figs. 4a and 4c, the effect of Biot number is not as obvious in Fig. 4a. The above results are similar to the results of Abdel-Hamid[8].
A new approach is applied for solving a transient heat conduction problem in spherical coordinates subject to a nonlinear boundary condition due to coupled convection-radiation heat exchange. Treating the dependent term in the boundary condition as a source is employed to obtain the associated homogeneous problem. Then the problem can be solved by any conventional method. It is shown that the proposed method provides a straight-forward methodology for heat equations subject to nonlinear boundary conditions. This approach yielded similar results as the explicit finite difference solution for all considered cases. It is obvious that the finite integral transform technique is a useful methodology that can be used to solve both linear and nonlinear heat conduction problems.