Для чисельного ршення задач нестащонарног теплопередачi розроблений спрощений дискретний аналог, представлений в обезрозмiреному вигля-дi. Цього вдалося досягти внаслидок використан-ня одновимiрноi вихидноп моделi. Використання таког моделi е достаттм для виршення бiльшостi практично важливих задач. Показана стштсть чисельного розрахунку при великих кроках дискре-тизаци за часом i висока точтсть розрахунтв на гранично малих розрахункових Ытках в 3 вузли
Ключовi слова: нестационарна теплопередача, плоска сттка, спрощений чисельний розрахунок,
мала розрахункова ытка
□-□
Для численного решения задач нестационарной теплопередачи разработан упрощенный дискретный аналог, представленный в обезразмеренном виде. Этого удалось достичь вследствие использования одномерной исходной модели. Использование такой модели является достаточным для решения большинства практически важных задач. Показана устойчивость численного расчета при больших шагах дискретизации по времени и высокая точность расчетов на предельно малых расчетных сетках в 3 узла
Ключевые слова: нестационарная теплопередача, плоская стенка, упрощенный численный расчет, малая расчетная сетка -□ □-
UDC 519.63
|DOI: 10.15587/1729-4061.2017.96090]
A SIMPLIFIED METHOD FOR THE NUMERICAL CALCULATION OF NONSTATIONARY HEAT TRANSFER THROUGH A
FLAT WALL
O. Brunetkin
PhD, Associate Professor* E-mail: [email protected] M. Maksymov Doctor of Technical Sciences, Professor, Head of Department* E-mail: [email protected] O. Ly s i u k Postgraduate student* E-mail: [email protected] *Department of thermal power automation processes Odessa National Polytechnic University Shevchenko ave., 1, Odessa, Ukraine, 65044
1. Introduction
The use in energy equipment of fuel with variable composition [1, 2] instead of the certified constant formulation leads to an increase in the number of transient modes. A change in calorific fuel capacity, amount of products of combustion, their thermophysical properties predetermines the non-stationary processes of heating-cooling the elements of power equipment design, as well as nonstationary processes of heat transfer through the heat exchange surfaces. As a consequence, in the processes of energy conversion a more important role is played by its variable accumulation - release in all equipment components in contact with the combustion products. The accumulation of energy (heating - cooling) affects inertness in the processes of heat transfer and, therefore, manageability of the process of energy transformations. In turn, the magnitude of energy accumulation is determined by temperature of the elements of design. Thus, determining the non-stationary temperature and, as a result, the nonstationary magnitude of energy accumulation, is an important element in solving the problem on optimal control of power equipment under conditions of using a non-certified fuel of variable composition.
2. Literature review and problem statement
Analytical solution is a reliable and universal means to calculate parameters of the nonstationary heat exchange
process. But at the same time, its application is cumbersome and covers only particular cases. As a rule, they are confined by the description of the process of heating-cooling the bodies of simple shapes: infinite plate, cylinder, and sphere. Such solutions have been known for a long time, but up to now they have been applied only for a similar kind of processes. Only the object of application has changed in favor of contemporary equipment as, for example, in [3], where the problem on cooling a display is solved. Article [4] examines a process of cooling the plate with an internal heat source. The difference is that the examined plate is a multilayer one. A special feature of paper [5] is the application, instead of the frequently utilized Fourier expansion, the Trefftz functions. In this case, only the approximate solution is obtained. The problem of nonstationary heat transfer is of practical interest and it does not have the analytical solution.
Numerical methods are more universal and can be used to solve any problems on heat transfer. From the above enumeration of simple shapes of bodies, a plate appears the most important one. It has two surfaces, which corresponds to the process of heat transfer. Furthermore, for the case of coaxial cylinders (a pipe) at the ratio between outer and inside diameters of d2/d1<2, heat transfer through the cylindrical wall with an error less than 4 % can be described using the model for a flat wall. Such a relationship of diameters corresponds to the majority of variants of tubular heat exchangers. Heat transfer is predetermined by taking account of the combined processes of heat exchange at the surfaces of a plate. These
©
processes in general are difficult to describe. But in many practically important cases this description is confined by assigning a generalized (integral) magnitude - heat-transfer coefficient. Here arises a contradiction between relative simplicity of problem statement on the nonstationary heat transfer through an infinite plate and the use of sophisticated universal numerical methods. Thus, [6] proposes a transition from the distributed model to the concentrated in the form of a system of ordinary first order differential equations. Its considerable effectiveness is emphasized in this case relative to the finite-difference and finite-element methods. A well-known fact is indirectly highlighted here of the need to maintain certain ratios between the size of computational grid and the calculation step by time: in order to provide for the necessary accuracy, the computational grid is lessened with a simultaneous decrease in step by time to ensure stability of the numerical calculation. This leads to an avalanche-type increase in the number of computations and to the accumulation of errors in calculations. For a practical application in the engineering calculations, in most cases, it is sufficient to consider a problem in the one-dimensional setting. This is justified in view of the smallness of thickness of the heat-transmitting surface in comparison with its other geometric dimensions. A confirmation of this can be found, for example, in [7]. It is demonstrated here that, when solving a nonstationary inverse problem on heat transfer, good results are obtained when using the one-dimensional model in particular.
A solution in the form of functional dependence is one of the advantages of analytical methods over numerical. It is necessary to obtain the solution only once. It can be subsequently used for any arguments. In this case, in the form of this dependence, the interrelation of all phenomena, taken into account in the model, is reflected analytically. Results of numerical calculations lack these properties. In order to somewhat improve this situation, it is possible to use the models and, accordingly, numerical solutions based on them in a generalized (made dimensionless) form. The obtained advantage manifests iteself most vividly in complicated cases. Thus, in [8], when examining the magnetohydrodynamic flows, using the similarity transformations, it was possible to reduce the Navier-Stokes equations to ordinary differential equations. Such a result is possible due to the representation, when rendering dimensionless, all terms of equations in the uniform form and obtaining the possibility of running a fractional analysis. In [9], it was possible in this way to obtain approximated analytical solutions of the systems of ordinary differential equations. However, this approach, as a rule, is not applied to the simpler problems, although it can also yield a positive effect.
Therefore, the nonstationary problem on heat transfer can be solved with the aid of numerical methods. A large number of similar practically important processes can be described using relatively simple one-dimensional models. The problem is in the mismatch between simplicity of the model and complexity of its proposed numerical realization. When employing simple models, additional possibilities to generalize the results of solution by representing the model in a dimensionless form are not used.
3. The aim and tasks of the study
The aim of present work is to develop a simplified discrete analog, designed to calculate the process of nonstation-
ary heat transfer through the infinite plate (one-dimensional model). This will make it possible:
- to represent the discrete analog and, therefore, results of the calculations in the dimensionless form with the possibility of their appropriate generalization;
- to obtain the discrete analog that makes it possible to carry out calculations with the required accuracy on the rough computational grids (with a small number of computational nodes) and at large steps in time retaining the stability of numerical calculations.
To achieve the set aim, the following tasks are to be solved:
- to substantiate the choice of method for constructing the simplified discrete analog;
- to construct a discrete analog in the dimensionless form;
- to prove working ability of the developed analog and evaluate the accuracy of calculations based on it by their comparison with the analytical solutions for the known particular cases;
- to determine the limits in the applicability of the analog at a decrease in the magnitude of computational grid and an increase in the steps of calculation in time.
4. Method for constructing the discrete analog
Underlying the solution of the set problem is the control volume method (CVM) [10, 11]. It combines the benefits of other numerical methods and is deprived of their many shortcomings. Thus, the desired discrete analog in CVM is built as easy as in the finite-difference method. The desired parameters are calculated in the isolated points (nodes). In this case, with the aid of appropriate profiles, similar to the method of finite elements, we consider possible character of change in the calculated parameter between the nodes of a computational grid. Moreover, a discrete analog is built based on the compliance with conservation laws in each particular control volume. This makes it possible to obtain physically noncontradictory results on the grid of any roughness, which, in contrast to the finite-difference method, for example [12], allows using computational grids of small size. The merits of CVM contribute to its application at present by many authors [13, 14].
The use in the development of a discrete analog of the a priori-profiles of change in the desired magnitude, similar to an exponential, agrees with the use of profiles in the form of exponents in the method of integral coefficients [15]. It should be noted that CVM underlies the construction of such a universal programming product as SolidWorks.
[10] emphasizes the development of universal two- and three-dimensional analogs. Only the principle of their construction is explained on the example of one-dimensional analog, but the possibility of solving the one-dimensional variant of the problem based on dimensionless magnitudes is not examined.
Universalism also manifests itself in the fact that even a one-dimensional analog includes the source term that considers sources or sinks of heat fluxes inside the heat exchange surfaces. In the overwhelming majority of cases, in the thermal-power equipment, the sources and the sinks of energy are found outside the body, through which heat exchange is conducted. Therefore, for the simpler algorithm of sources, the term should be excluded from consideration.
Within the framework of CVM, the examined space, similar to the method of finite elements, is split into separate small elements - control volumes. In a three-dimensional orthogonal coordinate system, these are the "orthogonal" elements, in the Cartesian - parallelepipeds, in the one-dimensional model - layers (Fig. 1).
Sxj = const; Axj = const; Sx = Ax.
(1)
ap ■ Tp = aE ■ TE + aw ■ TW + b,
where
_ Ax _ l _ l
a7) — ar- + a^r + p ■c ■ ; at- — ; a^r — ; P E W H At E Sx W Sx
u Ax o
b — p " Tp.
(2)
(3)
Here p, c, X are the density, heat capacity and thermal conductivity of material of the heat exchange surface; At is the step of calculation by time; TP, TE, Tw are the calculated (current) temperatures in the corresponding points; TP is the value of temperature in point P from the previous step
of calculation by time. At the first step of calculation is the value from the initial condition (initial temperature profile).
Let us substitute (3) into (2), divide all terms of the equation by X/Sx. As a result, taking into account (1), we shall obtain
2 + 1(Ax)2
a At
T — T + T +1 (Ax) T P = T + Tw + a ^ At ^
(4)
Fig. 1. Calculation scheme for the one-dimensional discrete analog
Inside these layers are the nodes, which, similar to the finite-difference method, are assigned with the values of calculated magnitudes. Geometric characteristics of the calculated region represented in this way are: Sx are the distances between the grid nodes; Ax are the dimensions of layers (elements), into which the region is split. Positions of the nodes inside different layers may differ. The very dimensions of layers may be different, too. That is why, in a general case, both Sxi and Axi are different from each other. Furthermore, in general, Sx^Ax. In the considered case, at constant thermophysical properties of heat exchange surfaces byn their thickness, in order to simplify the discrete analog, let us assume:
where a=X/( p-c) is the coefficient of thermal diffusivity.
Expression (4) is the one-dimensional discrete analog in dimensional form for solving the problem on nonstationary thermal conductivity. In order to bring it to the dimension-less form, let us assume that thickness of the heat exchange surface is equal to 2l. Let us multiply second term in brackets by (2l)2/(2l)2=1. Designate in this expression Ax/2l= =Sx/2l=A. It is the dimensionless (relative) thickness of layer of the discrete analog. We shall obtain, taking into account further transformations:
1 (Ax)2 _ (2l)2
a At
a At
(A)2 —
1
A(Fo)
■(A)2.
(5)
Here A(Fo) is the dimensionless step of calculation by time - a step of a change in the Fourier number.
After carrying out analogous transformations for the last addend in the right side of expression (4), we shall obtain:
2 + -
1
A(Fo)
(A)2
■T — T + Tw +-
1
A(Fo)
(A)2Tp0. (6)
In CVM, on the boundaries of calculated region, there are the nodes of computational grid (points "1" "N" Fig. 1). They are surrounded by incomplete boundary control volumes. Taking into account (1), their magnitude in the examined case is Ax/2. This must be considered when recording the discrete analog for nodes on the boundaries of calculated region. Furthermore, for these nodes, in contrast to the internal, it is necessary to consider boundary conditions at the corresponding surfaces. Thus, discrete analogs for the internal and boundary nodes do not differ, but they must have common character of their record for the possibility of their regular solution.
Let us examine a discrete analog for the internal points. Based on common constructions from [10] for internal point "Pi", taking into account parameters in points "Wi" and "Ei" (Fig. 1), it is possible to write in order to calculate temperatures:
Expression (6) is the dimensionless discrete analog for internal points of the computed region.
Let us examine a discrete analog for the node of computational grid on the boundary of calculated region, for the certainty in point "1" on the left boundary (Fig. 1). In order to compute temperature in point P1, it is necessary to consider boundary conditions to the left of it, as well as value in point E1. As the boundary conditions, let us examine conditions of the third kind as the most general ones. Based on general constructions from [10], it is possible to write for the boundary point "1"
ap-tp — ae ■ te + b where
(7)
aP — aE + a1 + p-c-—; aE — —;
Ax; At;
u T Ax 1 0
b—«v Tambl +p c—■ -■ TP.
(8)
Here a1, Tambi are the heat transfer coefficient and ambient temperature from the corresponding side of a heat exchange surface.
Similar to the previos case, let us substitute (8) in (7), divide all terms of the equation by X/Sx. As a result, taking into account (1), we shall obtain
Sx 1 (Ax)2 1
1+a--+-■ -—----
1 l a 2 At
T —
P
TSx ™
E +Kr T ^ Tamb''
1 (Ax)2 X To a 2 At P.
Let us multiply in expression (9) the terms, which contain Sx, by (2l)/(2l)=1. The terms that contain (Ax)2 shall be multiplied by (2l)2/(2l)2=1. Consider that Ax/2l=Sx/2l=A. Let us designate as
«r2l _Bi — ~ Bi1
the Biot criterion from the side of node "1" of the computational grid. As a result, we shall obtain:
— 11
1 + Bii(D) + ---—
2 D(Fo)
(D)2
T =
1
= TE + Bil ■ Tambi ■ (D) + -
2 D(Fo)
■(D)2 ■TP.
(10)
Expression (10) is the dimensionless discrete analog for the left boundary point of computed region (Fig. 1).
Upon carrying out transformations, analogous to (7)-(10), but for the right boundary (node "n", Fig. 1), we shall receive the dimensionless discrete analog for the right boundary point of computed region:
— 1 1 1 + Bi2 (D) + -■
2 D(Fo)
(D)2
T =
= Tw + Bi2 ■ TamU ■ (D) + -
1 1
2 D(Fo)
■(D)2 ■TP.
ai ■Tj = bj ■Ti+1 + Cj ■Ti-1 + dj, where
1 1
a = 2 + ■(D)2; bi = 1; ci = 1;di ■(D)2^TP0;
D(Fo) - for the left boundary a1 ■T1 = b1 ■T2 + d1, where
D(Fo)
a1=1+Bi-(D)+2 D(Fo)(D)2; b1=1; d1=^w(D)+ÎD(k)-(D)2-Tp;
- for the right boundary
aN ■ TN = CN ■ TN-1 + dN, (14)
where
aN=1+BÎ2-(D)+2D(Foy(D)2; CN=1
dN = Bi2 Tamb2 (D) + 1 D(Fo)(D)2 Tp.
The solution algorithm is built based on [10] as follows. Direct course of computation - auxiliary coefficients P and Q are computed:
first from (13) P1 and Q1
(11)
P = b =
m _ _ a
1
1 1+Bl1(D)+Îd(fo)-(d)2'
Q d1 Bi1Tamb1(d) + 2d(fo)(D)2 tp0 Q1 = = "
(15)
^ 1 + Bi1 (D) + ;
Here Bi2, Tam^ is the Biot criterion and ambient temperature from the side of node "N" of of the computational grid.
5. Form of the discrete analog and algorithm of solution based on it
Expressions (6), (10) and (11), taken together, represent the dimensionless discrete analog for calculating the non-stationary heat transfer through a flat plate at boundary conditions of the third kind.
The algorithm of solution based on such discrete analog can be realized with the help of TDMA (Tri-diagonal-Matrix Algorithm) - the sweep method. Results of calculation are obtained in one "sweep", without iterations, which simplifies the computation. In order to facilitate the realization of the solution algorithm, let us write discrete analog (6), (10), (11) in the index form, where the indices are counted from the left to the right boundary of calculated region (Fig. 1):
- for internal points
2 D(Fo) then from (12) and (15) Pi and Qi b1
(D)2
Pi= _ D " 1 '
ai Ci i-1 2 + -^ (D)2 - Pi.1 D(Fo)
, ——(D)2TP0 + Qi1
Q = di + CiQi-1 = D(Fo) 7 P ^i-1
Qi „ „ r> 1
(D)2 - Pi-1
(16)
ai -Ci Pi-1 2+ -A-
D(Fo)
Inverse course of computation - for the right boundary bi=0. Consequently, pn=0. From (16):
Qn =
dN + Qi-1 :
aN - Pi-1
(12)
Bi2 ^ Tamb2 ■ (D) + 2 D(Fo) (D)2 TP + Qn-1
= — 1 1 .
1 + Bi2 (D) + ------(D)2 - PN 1
v J 2 D(Fo) N-1
Assume from (17)
Tn = Qn. and then in reverse order
Ti = PTi+1 + Qi.
(17)
(18)
(19)
(13)
Simplicity of the obtained algorithm (15)-(19) makes the realization possible even based on Excel, to say nothing of the more powerful computational tools.
6. Evaluation of adequacy of the obtained results by the available solutions. Possibilities and constraints in the application of the obtained discrete analog
An algorithm must qualitatively and quantitatively correct reflect the processes under investigation. Given this, the estimation of physicalness and adequacy of the developed discrete analog can be carried out in two ways:
1) after appropriate transformation, by the comparison of analog itself with the exact solution in the extreme case of stationary heat transfer through an infinite flat wall (estimation of physicalness - qualitatively correct representation of the examined process);
2) by comparing the results of numerical calculations with the results of analytical solutions in a particular case of the nonstationary process of symmetrical heating (cooling) of an infinite plate (estimation of correctness of quantitative representation of the studied process).
For both these cases, there are precise analytical solutions.
1. Stationary heat transfer through a flat wall
When solving the problem on stationary heat transfer through a flat wall, the heat flux through both surfaces of the wall is assumed identical. For certainty, let us accept (Fig. 2) that the ambient temperature to the left of the wall exceeds temperature to the right, and the heat flux is directed from left to right.
0(2
(Ba)
(Tamb1 Tw! )
Bi,
<X, - Tamb2 ) «1 BV
Fig. 2. Minimal computational grid
This can be written in the form:
- q1 =a1 ■ (Tamb -Tw ) is conditionaly for the left side of the wall;
- q2 = a2 ■ (Tw -Tamb ) is conditionaly for the right side of the wall.
Here:
qi, q2 is the heat flux, respectively, entering from the environment from the left to the wall and exiting the wall from the right to the environment;
a1, a2 are the heat transfer coefficients on the left and right side of the wall, respectively;
Tw1, Tw2 is the temperature of the left and right side of the wall, respectively;
(Tamb1, Tamb2) is the ambient temperature from the left and from the right sides of the plate, respectively.
Taking into account the condition accepted q1=q2, we shall obtain:
«1 ■ (Tamb1 - Tw1 ) = «2 ■ (Tw2 - Tamb2 )
Here BI1, Bi2 are the Biot criteria for the left and right side of the wall.
The discrete analog in the form (17)-(19) is obtained for the case of nonstationary heat transfer. But at constant ambient temperatures from the left and right sides of the wall and sufficiently long course of the heat transfer process, the analog, if correctly constructed (17)-(19), must yield the result, analogous to (20). Furthermore, in the examined case, the profile of a change in temperature from Tw1 to Tw2 (Fig. 2) inside the wall must be of linear character, and the correctly constructed discrete analog also must reflect it. Let us transform the analog (17)-(19) to test the feasibility of these requirements.
Let us split the examined flat wall into three layers (Fig. 2): two near-wall ones and one internal. In this case, temperatures T1, T3 of the analog correspond to temperatures Tw1, Tw2 at the surface of the wall. Temperature T2 is the temperature in the central layer. The near-wall layers in accordance with (Fig. 1) have a half thickness. The choice of only three layers for the splitting is predetermined by the following considerations. The discrete analog for temperatures Ti in central cells (17) is built with the use of temperatures in the adjacent cells Ti-1 and Ti+1. Three layers allow us to construct a similar analog for temperature T2 taking into account temperatures Tw1,Tw2. The near-wall layers make it possible to employ analogs (18), (19) for the left and right boundaries of the plate. Thus, in order to evaluate correctness of representing the process of stationary heat transfer with the help of the developed discrete analog, three calculated layers will suffice and an increase in the number of inner layers at partition adds nothing fundamentally new to the computation.
In the analog (17)-(19), magnitudes Tp°, Tp°, Tp°, utilized in coefficients di, d1, dn, are the temperatures in the corresponding nodes of computational grid from the previous step of calculation by time. In the computation of heat transfer process over a prolonged time interval and its reaching the stationary state, temperatures from the previous step of calculation by time must be equal to the estimated temperatures currently in the corresponding points. We obtain for the case of three layers in question:
Tn° — T; Tn° — T; Tn° — T3.
p, 2 ' p1 1 ' pq 3
(21)
We shall obtain from (17) upon the substitution of expressions for all coefficients:
2 T2+A(k) (A)2 T2—T+T1+A(k) (A)2 T- (22)
Taking into account equality Tp0 — T2 from (21), we shall receive the equality of terms from the left and right sides of expression (22):
1
(A)2 T2 — —■(A)2 ■TP ,
A(Fo) 2 A(Fo)
and after their reduction:
or
2 ■ T2 = T3 + T1 or T2 =
t3+T
(23)
After carrying out analogous transformations for the analogs of left (18) and right (19) boundaries, we shall obtain:
- for the left boundary
T + Bir(D)^1 = T2 + Bi^D)^;
- for the right boundary
T3 + Bï2(D)T3 = T2 + Bi2 ■ (D) ■ Tam..
(24)
(25)
Substitute expression for T2 from (23) to (24) and (25). Carry out transformations and, as a result, we shall obtain: - for the left boundary
^ = Bir(D)T1 - Bi^D)!^; - for the right boundary
= BV(D)Tamb2 - Bi^D)^.
(26)
(27)
In expressions (26) and (27), the left sides are equal. Equate in these expressions the right sides, reduce by A:
Bir(T1 - Tamb1) = Bi2 (Tamb2 - T3) BV(Tamb1- T1) = Bi2 ^(T3 - Tamb2)
we receive as a result:
(Tamb1 - T1) = Bl^ = (T3 - Tamb2 ) Bi1 «1.
(28)
(29)
tg(<M = ^^ = Bi1(Tamb1 - T1),
(30)
- from (Fig. 2) for the section between temperatures (T2-T3) taking into account (25):
tg(^2,3 ) = ^^^ = Bi2^T3 - Tamb2).
(31)
The equality of right sides in expressions (30) and (31) follows from (28). Hence, the left sides are equal in them:
Thus, in the process of stationary heat transfer the sections of profiles of temperatures, represented using the analog (17)-(19), have a common point T2 and identical angles of inclination (32). Consequently, they are one straight line. This agrees with the analytical solution and confirms correct construction of the analog (17)-(19) in this part as well.
2. Symmetrical heating of an infinite plate.
The accuracy of numerical calculations, executed with the help of the discrete analog proposed, can be evaluated by comparing their results with the available analytical solutions. The case of bilateral symmetrical heating of an infinite plate is one of a few that exist. Let us examine a variant when at the initial moment of time (t=0), initial temperature in the plate is distributed evenly. Under these conditions, the proposed analytical solution takes the form of summation of a series. The indicated articles note that at Fo>0.3 a series starts so rapidly converging that the temperature distribution is determined accurately enough by the first term of a series in the form:
0 = 1 --
2 ■ sin
+ sin ■ cos
■ cos^- X) exp(-^42Fo). (33)
Here © is the relative (dimensionless) temperature of a plate; X is the relative (dimensionless) coordinate of the point in question. It is counted from the center of the plate to its surface;
Fo=(a-t)/(S2) - Fourier number;
where a is the coefficient of thermal diffusivity; 5 is the half thickness of the plate.
Relative temperature of plate © is determined by relationship
0 =
T(X) - To T - T
1 amb 1 0
(34)
Comparing the expressions (29) and (20) reveals their concurrence, which confirms correct construction of the analog (17)-(19) in this part.
Let us consider the second part of the test - representation with the help of the proposed discrete analog of the form of temperature profile inside a flat wall. For this purpose, let us determine the tangents of temperature profile angles of inclination in sections (Tj-T2) and (T2-T3):
- from (Fig. 2) for the section between temperatures (Ti-T2) taking into account (24):
where T(X) is the current temperature in the corresponding point of plate; T0 is the initial temperature, evenly distributed in the plate; Tamb is the ambient temperature, to the magnitude of which the plate is heated.
Relative temperature of the plate changes in the range of ©e [0...1].
Relative coordinate is determined from relationship X=x/S and changes in the range of Xe[0...1]. Here x is the absolute coordinate that is counted from the center of the plate to its surface.
Magnitude | is the root of transcendental equation ctg(|)=|/Bi. Equation (33) is the first term of a series. Therefore, it employs | - the first positive root of the transcendental equation.
In order to compare results of the analytical and numerical calculations, one should consider that:
- biot criterion in the analytical solution is computed for the half thickness of the plate (because of its symmetry). However, in the numerical calculation, taking into account the possibility of solving the non-symmetrical problems, it is determined for the full thickness of the plate. Thus, it is necessary to apply relationship
tg(^1,2 ) = tg(^2,3>.
(32)
(BiJnc = (Bi2)nc = 2 Bia
(35)
or
where (Bi1)nc, (Bi2)nc are the Biot criteria for the corresponding sides of the plate in the discrete analog (numerical calculation) for the case of symmetrical heating; (Bi)a is the Biot criterion in analytical calculations;
- the Fourier number due to the above noted reason also differs for the cases of analytical and numerical solutions. Comparing the results, which correspond to the identical moments of dimensionless time, should be conducted based on relationship:
4 ■ Fonc = Foa. (36)
where Foa is the Fourier number in the analytical calculation; Fonc is the corresponding Fourier number in the numerical calculation.
Fig. 3 shows results of computation based on analytical expression (33) and discrete analog (12)-(14). Since the analytical calculations are carried out for the dimension-less temperature (34), then the numerical solutions are obtained based on it. Fig. 3 shows results that cover the broad range of change in the Biot criterion. In the brackets is its value for the numerical calculations based on (35). Each figure presents results for the two moments of time. Magnitudes of the Fo numbers corresponding to them are bound by relationship (36). Maximum relative error e in the numerical results from all estimated points by thickness of the plate is given for each moment of time. All in all, 21estimated points were considered. Relative error was determined relative to the range of change in temperature ©e[0...1].
A comparison of results of the numerical and analytical calculations reveals their good concurrence. However, similar accuracy at a large number of estimated points can be attained by other numerical methods. The employed method of control volumes differs by abiding the conservation laws on computational grids of any accuracy. In order to evaluate an influence of the number of estimated points on the computation error, we performed calculations using the maximally small grids - with three nodes only (Fig. 2). The obtained results and their comparison with the analytical and numerical calculations on a large grid (21 nodes) are given in Table 1. Here, similar to Fig. 3, the magnitudes of Bi and Fo that relate to the numerical calculations are given in brackets. The values of Bi and Fo in the analytical and numerical calculations, similar to the previous case, are bound by relationships (35), (36). The first three lines show the values of relative temperatures for Fo=0.4(0.1). In this case:
- in the first line are the temperatures in the analytical calculations;
- in the second line - in the numerical calculations on a grid with 21 nodes;
- in the third line - in the numerical calculations on a grid with 3 nodes.
The next three lines show relative errors in determining the relative temperatures, represented in the previous lines. In this case:
- in the line, designated by s21, are the errors of results in the numerical calculations on a grid with 21 nodes relative to the analytical calculations;
- in the line, designated by s3, are the errors of results in the numerical calculations on a grid with 3 nodes relative to the analytical calculations;
- in the line, designated by s21-3, are the errors of results in the numerical calculations on a grid with 3 nodes relative to the numerical calculations on a grid with 21 nodes.
Further in the table, the order of results arrangement is analogous.
■-©a. Fo=5 -©a. Fo=l
---©nc, Fo=l.25,e=0.01 ---©nc.Fo=0.25.e=0.01
■—-©a, Fo=30 -©a. Fo=4
---©nc,Fo=7.5,£=0.03 — ©nc.Fo=l.e=0.01
0.8 0.6 0.4 0.2 0.0
0.0 0.5 1.0 0.0 0.5 1.0
a b
-©a. Fo=0.4 -©a. Fo=0.4
---©nc,Fo=0.1,£=0.001 ---©nc.Fo=0.1.e=0.01
-©a. Fo=l -- ©ai Fo=l
Fig. 3. Relative temperature © dependent on the relative X coordinate and relative moment of time Fo for the analytical ©a numerical ©nc calculations: a — Bi=0.004 (0.008); b - Bi=0.5 (1.0); c - Bi=5 (10); d- Bi=50 (100)
Relative errors do not exceed 3.5 % in all given points. The magnitudes of Bi=5 and Bi=50 are selected because of the maximal relative error received in this case. At Bi=0.04 and Bi=0.5, an error in the numerical solution does not exceed 1 %. The magnitudes of the noted errors do not exceed permissible values for the engineering computations (<5 %). This testifies to the applicability of the developed discrete analog for solving the problems on nonstationary heat transfer.
Let us examine an example of solving a problem on the nonstationary heat transfer. The absence of analytical solution is its special feature. Assume that at the initial moment of time body temperature and ambient temperature are in equilibrium. Let us designate it similar to the previous cases, T0. At a certain moment of time, temperature of the environment from one side of the plate rises abruptly to magnitude Tamb From the other side of the plate, ambient temperature remains equal to T0. We shall examine a nonstationary process of change in the temperature inside the plate during heat transfer. We shall perform the computation for relative temperature ©, determined in accordance with (34) using the relative node coordinates of grid X. Some results for the computational grid with 21 nodes at different values of the Fo numbers are shown in Fig. 4. The value of Fo=0 corresponds to the initial moment of time.
Table 1
Results of the analytical and numerical calculations of relative temperature © at symmetrical heating of an infinite plate. Magnitudes of error in the numerical calculations relative to the analytical ones
Bi Fo/e X
0.0 0.5 1.0
0.4 0.842 0.378 0.842
(0.1)21 0.837 0.366 0.837
(0.1)3 0.873 0.393 0.873
e21 0.005 0.012 0.005
e3 0.031 0.015 0.031
5 e21-3 0.036 0.027 0.036
(10) 1.0 0.944 0.779 0.944
(0.25)21 0.941 0.766 0.941
(0.25)3 0.955 0.765 0.955
e21 0.003 0.013 0.003
e3 0.011 0.014 0.011
e21-3 0.014 0.001 0.014
0.4 0.985 0.507 0.985
(0.1)21 0.984 0.487 0.984
(0.1)3 0.991 0.528 0.991
e21 0.001 0.020 0.001
e3 0.006 0.020 0.006
50 e21-3 0.007 0.041 0.007
(100) 1.0 0.996 0.881 0.996
(0.25)21 0.996 0.868 0.996
(0.25)3 0.997 0.851 0.997
e21 <0.001 0.013 <0.001
e3 0.001 0.03 0.001
e21-3 0.001 0.017 0.001
results on a grid with 21 nodes (Table 2). The structure of data here and in Table 1 is analogous.
Table 2
Results of the numerical calculation of relative temperature © at nonstationary heat transfer through an infinite plate. The magnitudes of error in the numerical calculations on a small computational grid
Fo ©/e X
0.0 0.5 1.0
(0.12) ©21 0.710 0.178 0.015
©3 0.718 0.177 0.024
e21-3 0.009 0.001 0.009
(0.26) ©21 0.790 0.336 0.047
©3 0.799 0.332 0.052
e21-3 0.009 0.004 0.005
Fig. 4. Profiles of relative temperatures © along the relative thickness of plate X depending on different Fo numbers
In accordance with (34), for the accepted initial data in this moment of time, the value of relative temperature in all points of the plate is ©=0. The value of Fo=1.7 models reaching the state of stationary heat transfer.
For the estimated calculations of the nonstationary heat transfer processes, of interest is the use of "rough" computational grids. Let us define a possibility of applying the method of control volumes at minimal computational grid, the same as in the case of symmetrical heating (with 3 nodes). For this purpose, let us compare results when using it to the
Errors in the calculations on the "rough" grid in comparison with results on the more detailed one are insignificant. This accounts for its use in many cases. Insignificant errors on such a "rough" grid can be explained based on the following considerations:
- on one hand, the profile of change in the temperature inside the plate, both at its heating and at nonstationary heat transfer, is of exponential character. On the other hand, as noted above, to determine temperature in the nodes, they use the profiles of its change between the nodes. In a general case, the profiles can be any, including as parts of exponents. In the considered case, in order to simplify the form of a discrete analog when using the implicit scheme of its construction, we accepted a spasmodic character of change in temperature between the nodes. Nevertheless, such a simplified character in the representation of change in the temperature profile between the nodes of a computational grid does not in general contradict the character of change in temperature in the real process. In other words, the discrete analog satisfies the physics of examined processes, which predetermines its accuracy on the "rough" computational grids;
- as it was noted above, the discrete analog is built based on abiding the conservation laws on the grids with any "roughness". In combination with the preceding point, this also contributes to improving the accuracy of computations.
It should be noted that in a general case stability and accuracy of numerical methods depend on the relationship between the step of calculation by time and a geometric dimension in the step of a computational grid. Thus, in order to provide for the stability of computation when using the explicit scheme in the finite-difference method, At is determined from relationship [10]:
Dt <
p^ c ■ (Dx)2 21
(37)
Upon performing the transformations, applied to (5), we shall obtain:
= D(Fo) < 1(D)2. 2l V 7 2 W
(38)
In the above examined case of computation in the relative coordinates with 21 nodes of computational
grid and, accordingly, 20 intervals between them, we receive A=0.05. In accordance with (38), when applying the finite-difference method, it would be necessary to obtain A(Fo)<0.0012. In the performed calculations, we employed A(Fo)=0.005, which actually exceeds the noted boundary.
Let us examine an impact of the magnitude of step by time on the stability and accuracy of calculations. Results at different steps are given in Table 3. Calculations were carried out for 21 nodes, but, for visibility, the results are given only for the points at the surfaces and in the center of the plate. In order to compare, we took the same values of Fo=0.12 and Fo=0.26, as those used in Table 2. Taking into account A(Fo)=0.005, the value of temperature at Fo=0.12 was obtained in 24 steps, which is reflected in the appropriate line A(Fo)24. The next lines contain results of the temperature calculation at the same 21 nodes and at the same moment of time Fo=0.12, but at other steps by time and, accordingly, different number of these steps (subscript): A(Fo)8=0.015 and A(Fo)4=0.03.
Table 3
Results of the numerical calculations of relative temperature © in the process of nonstationary heat transfer through an infinite plate at different calculation steps by time
Fo A(Fo) X
0.0 0.5 1.0
(0.12) (0.005)24 0.710 0.178 0.015
(0.015)8 0.704 0.173 0.015
(0.03)4 0.694 0.167 0.015
(0.26) (0.005)52 0.792 0.336 0.047
(0.02)13 0.786 0.328 0.046
(0.065)4 0.774 0.306 0.042
contains relative errors in the rough calculations relative to those accepted as reference. The second part of the table contains analogous data for the moment of time Fo=0.26. The errors calculated demonstrate that when using the discrete analog based on the method of control volumes within the limits of engineering accuracy of computations, it is possible to employ maximally "rough" grids and large steps by time.
Table 4
Results of the numerical calculations of relative temperature
© in the nonstationary heat transfer process through an infinite plate at different combination of the number of nodes in a computational grid and the magnitude of calculation steps by time
Fo A(Fo)/e X
0.0 0.5 1.0
(0.12) (0.005)24 0.7097 0.1775 0.0145
(0.03)4 0.6863 0.1684 0.0230
£ 0.023 0.009 0.009
(0.26) (0.005)22 0.7920 0.3362 0.0470
(0.065)4 0.7797 0.3035 0.0473
£ 0.012 0.033 <0.001
Errors in the calculations at A(Fo)8 and A(Fo)4 relative to the variant at A(Fo)24 are not indicated, though they can be computed easily. The second part of Table 3 contains analogous data, but for the moment of time Fo=0.26 at steps, corresponding by the magnitude, and their number: A(Fo)52=0.005, A(Fo)13=0.02 and A(Fo)4=0.065. A comparison of the given data demonstrates maintaining the high accuracy of calculations also during the partition of the assigned time interval into a maximally small number of calculation steps.
Tables 2, 3 give the results that allow us to argue about the stability of calculations when using the "rough" computational grids and large steps by time. We, however, examined the influence of these factors separately. Of practical interest is the possibility to perform calculations with engineering precision at simultaneous influence of both factors. Based on the example of data, given in Table 4, it is possible to estimate the results of this kind of calculations when comparing to the case of detailed grid and small step by time. Results of the computation on a grid with 21 nodes and at calculation step by time A(Fo)=0.005 are accepted as reference. The first part of the table, as previously, considers time moment Fo=0.12. In order to reach this time at selected step, 24 steps are required. Results for this case are represented in the first line of the table. The second line contains results of the calculation on a grid with 3 nodes and at step by time A(Fo)=0.03, requiring 4 calculation steps to achieve the considered time. The third line
The adequacy of computation on the "rough" grids can be used also in the numerical-analytical calculations of the nonstationary heat transfer processes. The thing is that the discrete analog on the grid with 3 nodes consists of the system of 3 linear algebraic equations. For the current time moment, it is possible to receive a sufficiently simple analytical solution of this system. However, the possibility of using large steps by time allows based on this solution rapid estimation of temperatures at the surfaces of a plate in the nonstationary heat transfer process of and, accordingly, heat fluxes on them. In line with the classical approach, the energy accumulated in the plate is determined by integrating the temperature profile in it. In the examined case, it suffices to find a difference in the heat fluxes at the surfaces of a plate.
7. Conclusions
1. The method of control volumes selected for the discrete analog is based on fulfilling the conservation laws on computational grids of any size, including the smallest ones. This corresponds in the best way to the set problem on developing a simplified discrete analog.
2. The use of the simplified, one-dimensional form of the record of discrete analog made it possible to represent it in the dimensionless form. Such approach contributes not only to a decrease in the number of calculations, but it also facilitates comparing their results with data of the existing analytical calculations. The latter, as a rule, are represented in the dimensionless form.
3. A comparison of results of the calculations based on the developed analog to the available analytical data demonstrates their good recurrence. Relative errors do not exceed, and they are even substantially less than the per-
missible engineering accuracy (<5 %). Furthermore, this is substantially lower than the errors in determining the initial data - heat-transfer coefficients.
4. Based on the numerical calculations, using large grids and small steps by time as the rererence, we evaluated stability and accuracy of calculations on the rougher grids.
We demonstrated maintaining the engineering accuracy on maximally small grids (to 3 nodes) and retention of stability and accuracy at large steps by time, substantially larger than the other numerical methods permit. Such special features of the discrete analog can be used to solve the inverse problems on heat exchange.
References
1. Brunetkin, A. I. Method for determining the composition of combustion gases when burned [Text] / A. I. Brunetkin, M. V. Maksy-mov // Scientific Journal Natsionalnho Mining University. - 2015. - Issue 5. - P. 83-90.
2. Maksymov, M. V. Model and method for determining conditional formula hydrocarbon fuel combustion [Text] / M. V. Maksymov, A. I. Brunetkin, A. V. Bondarenko // Eastern-European Journal of Enterprise Technologies. - 2013. - Vol. 6, Issue 8 (66). -P. 20-27. - Available at: http://journals.uran.ua/eejet/article/view/18702/17074
3. Karvinen, R. Use of Analytical Expressions of Convection in Conjugated Heat Transfer Problems [Text] / R. Karvinen // Journal of Heat Transfer. - 2012. - Vol. 134, Issue 3. - P. 031007. doi: 10.1115/1.4005129
4. Shupikov, A. N. Nonstationary Heat Conduction in Complex-Shape Laminated Plates [Text] / A. N. Shupikov, N. V. Smetankina, Y. V. Svet // Journal of Heat Transfer. - 2007. - Vol. 129, Issue 3. - P. 335. doi: 10.1115/1.2427073
5. Grysa, K. Trefftz Functions Applied to Direct and Inverse Non-Fourier Heat Conduction Problems [Text] / K. Grysa, A. Maciag, J. Adamczyk-Krasa // Journal of Heat Transfer. - 2014. - Vol. 136, Issue 9. - P. 091302. doi: 10.1115/1.4027770
6. Seem, J. E. Transfer Functions for Efficient Calculation of Multidimensional Transient Heat Transfer [Text] / J. E. Seem, S. A. Klein, W. A. Beckman, J. W. Mitchell // Journal of Heat Transfer. - 1989. - Vol. 111, Issue 1. - P. 5. doi: 10.1115/1.3250659
7. Vahabzadeh, A. Analytical investigation of the one dimensional heat transfer in logarithmic various surfaces [Text] / A. Vahabza-deh, M. Fakour, D. D. Ganji, H. Bakhshi // Alexandria Engineering Journal. - 2016. - Vol. 55, Issue 1. - P. 113-117. doi: 10.1016/ j.aej.2015.12.027
8. Mahapatra, T. R. Dual Solutions in Magnetohydrodynamic Stagnation-Point Flow and Heat Transfer Over a Shrinking Surface With Partial Slip [Text] / T. R. Mahapatra, S. K. Nandy, I. Pop // Journal of Heat Transfer. - 2014. - Vol. 136, Issue 10. - P. 104501. doi: 10.1115/1.4024592
9. Zhang, L.-Z. An Analytical Solution to Heat and Mass Transfer in Hollow Fiber Membrane Contactors for Liquid Desiccant Air Dehumidification [Text] / L.-Z. Zhang // Journal of Heat Transfer. - 2011. - Vol. 133, Issue 9. - P. 092001. doi: 10.1115/1.4003900
10. Patankar, S. Chislennye metody resheniya zadach teploobmena i dinamiki zhidkosti [Text] / S. Patankar. - Moscow: Ehnergoat-omizdat, 1984. - 152 p.
11. Brunetkin, A. I. Modified method of controlling the amount used in the solution of nonlinear problems in fluid dynamics with free surface in tanks of complex shapes [Text]: proc. of the VI symp. / A. I. Brunetkin, V. N. Nakosin // The vibrations of elastic structures with a liquid. - Novosibirsk: Siberian Research Institute of Aviation them. S. A. Chaplygin, 1990. - P. 26-30.
12. Kuznetsov, G. V. Numerical Simulation of Convective Heat Transfer Modes in a Rectangular Area With a Heat Source and Conducting Walls [Text] / G. V. Kuznetsov, M. A. Sheremet // Journal of Heat Transfer. - 2010. - Vol. 132, Issue 8. - P. 081401. doi: 10.1115/1.4001303
13. Lapka, P. Immersed Boundary Method for Radiative Heat Transfer Problems in Nongray Media With Complex Internal and External Boundaries [Text] / P. Lapka, P. Furmanski // Journal of Heat Transfer. - 2016. - Vol. 139, Issue 2. - P. 022702. doi: 10.1115/1.4034772
14. He, Y.-L. Numerical Solutions of Nano/Microphenomena Coupled With Macroscopic Process of Heat Transfer and Fluid Flow: A Brief Review [Text] / Y.-L. He, W.-Q. Tao // Journal of Heat Transfer. - 2015. - Vol. 137, Issue 9. - P. 090801. doi: 10.1115/1.4030239
15. Brunetkin, A. I. Integrated approach to solving the fluid dynamics and heat transfer problems [Text] / A. I. Brunetkin // Odes'kyi Politechnichnyi Universytet. Pratsi. - 2014. - Issue 2. - P. 108-115. doi: 10.15276/opu.2.44.2014.21