CC BY
38
UDC 519.633.6, 519.688
Computational Scheme for Solving Heat Conduction Problem in Multilayer Cylindrical Domain
A. S. Ayriyan*, J. BuSa Jr.*, E. E. Donets*, H. Grigorian*S, J. PribiS*
* Laboratory of Information Technologies Joint Institute for Nuclear Research Joliot-Curie str. 6, Dubna, Moscow region, Russia, 141980 ^ Department of Mathematics and Theoretical Informatics Technical University of Kosice Letna str. 9, Kosice, Slovak Republic, 04001 * Laboratory of High Energies Joint Institute for Nuclear Research Joliot-Curie str. 6, Dubna, Moscow region, Russia, 141980 S Department of Theoretical Physics Yerevan State University Alek Manoogian str. 1, Yerevan, Armenia, 0025
The computational scheme for solving heat conduction problem with periodic source function in multilayer cylindrical domain is suggested. The domain has a non-trivial geometry and the thermal coefficients are non-linear functions of temperature and have discontinuity of the first kind at the borders of the layers. The computational scheme is based on an algorithm for solving difference problem using the explicit-implicit method. The OpenCL realization of the suggested algorithm for calculations performed on a GPU is also compared to calculations performed using a CPU. It is shown that the scheme can be successfully applied to simulations of thermal processes in pulsed cryogenic cell, which is intended for pulse feeding the working gases into the working space of the ion source within the millisecond range. The results are given for a simulation of one of the particular cell structures, which is assumed to correspond to the practical realization. The computational scheme can be used for the optimization problem of the cell model parameters.
Key words and phrases: heat evolution, periodical heating source, multilayer cylindrical structure, finite-difference scheme, OpenCL realization.
Introduction
In modern science and technology of thermal conductivity, there is a very common phenomena for the study of objects with complex geometric and physical structure. The main goal of this work is to suggest a computational scheme for solving heat conduction problem with a periodic source in multilayer cylindrical domain. This problem describes the thermal processes inside a pulsed cryogenic cell in frame of previously suggested model [1,2]. The function of this cell is pulse feeding (in the millisecond range) the working gases into the working space of the ion source [3].
4-
z 'o
max 0
Figure 1. Schematic view of the object slice. The slice of the object: 0 — cooler, 1 — electrical insulator, 2 — heat source (conductive layer), 3 — external insulator, 4 — liquid helium temperature terminal with T = 4.2 K
2
1
r
z
*
*
*
0
r
r
r
0
2
z
Received 13th January, 2015.
The work was supported by RFBR grants 14-01-00628 and 14-01-31227.
The thermal processes can be described by the following system of parabolic partial differential equations with temperature depended coefficients [4]:
^»mf =1 I K<T> §) + I (X~<T> i)+X-(T, '>, «
where r G [0, rmax( 2)], z G [0,2max(f)] (or (r, z) G Q) and i > 0. The domain consists of different layers with different densities and thermal coefficients; thus, the index m is corresponds to a layer. The source function in Eq. (1) is Xm(T) = 0 for the layers m = 0,1, and 3 (there is no source) and has a periodical time dependence.
In a common case, the thermal coefficients are nonlinear functions of the temperature and the spatial coordinates with discontinuities of the first kind (for m +1 layers, there are m points of discontinuities). The initial and the boundary conditions are taken as
dT
— = 0 V (r, z) G 5tt \{(r, z):z = Zmaxj,
d n
T = To
V (r, z) G {(r, Z) : Z = Zmax},
(2)
T(r,z,t = 0) = To, V (r, z) G n,
where SQ is the boundary of the Q, and n is the normal vector of 5Q. The temperature at the right side is always equal to T0 because of contact with thermostat.
The parameter p and the functions cv, A and Xi = X(Ti) have discontinuities of the first kind at the following surfaces with radii: r*, r*, and r* in the interval [0...rmax\. Conjugation conditions between materials are considered to be ideal,i.e. at the discontinuity points, the temperature and the heat currents from left side and right side are equal.
1. Numerical Sheme
The initial-boundary value problem Eqs. (1)-(2) has been approximated by the following mixed explicit-implicit finite difference scheme [5,6]:
pij ^ T%'3 T%'3 = A [ Titj ] + Aj [ Th] ] + xit.
(3)
where Ti,j is temperature on the next time step, Ti,j - temperature on the current time step, is time-step.
Numerical solutions of Eq. (3) can be obtained using a special non-uniform grid [5]:
u = {(t,x, z) | 0 < TO, ti = k ■ ht, k G N0;
0 < r < fmax, Ti+1 = Ti + hi+i, i = 0, . . . , Nj - 1; (4) 0 < Z < Zmax, Zj + i = Zj + T]j+i, j = 0, . . . , Mi - 1}.
The spatial finite difference operator is:
A [T"] = r- 5
ri+2X+2-^--ri-ïï i-1 —
i-1,3
A ] = ^: 'h
Xi 4
Ti
i,j+2 "
i+l
T
— 2
hi
Ti,j Ti,j-1
(5)
(6)
2 V3+1 2 m
where i = l,...,Nj — 1, j = l,...,Mi — 1, hi = ri — ri-1, r]j = Zj — Zj-1 hi = (hi+i + hi) /2, ti j = ( 77^+1 + r]j ) /2, Ti,j = T (ri, Zj, tk ), = T (n, Zj, tk+i),
Ai,j = Xm(Ti,j), cvi,j = cvm(Ti,j), Xi,j = Xm(Ti,j), ri± 1 = (n + n±i)/2, K± 2 ,j = ^m(Ti,j + Ti±i,j)/2, Xi,j± 1 = \m(Ti,j + Ti,i±i)/2.
The indices « and j are global for the whole computational domain. The index to, which marks the material is chosen correspondingly to the domain region, from which the pair (i, j) is taken.
The difference equation (Eq. (3)) could be soved by the Thomas algorithm [7]. For that, one needs to initialize values a0, fi0, and TN.,j for the forward sweep and back sweep of the method respectively. These values have to be initialized to satisfy the boundary conditions (2):
ao = 1, fa = 0.
~ -i (7)
1Ni ,j = Z-.
3 1 - aN,-i
Note that in the case of a non-uniform grid or discontinuities of the first kind of the thermal coefficients, Scheme (3) has the first order difference approximation via spatial coordinatesi [4,5]. Therefore, for approximations of the boundary and conjugation conditions, it is enough to use the difference approximation of the first order.
The coefficients for the forward sweep of the Thomas algorithm ai* and pi* at the discontinuity point (on the border between layers) are given as follows:
( \ * h ■
Am+iai*
ai* =
Pi* =
Xmhi* + i@i*-i
Kn+l^ hi* + ^m^i* +i (1 — ai*-i)
(8)
The difference scheme (3)-(8) has unconditional stability related to spatial step hi and conditional stability related to ^ [6]:
min YnA t < -j121 min
p(T, r, z)cv(T, r, z)
X(T, r, z)
(9)
Generally, the conditional stability could be a strong limitation for the practical usage of a difference scheme. However, our scheme is practical in the cases where spatial step is sufficiently large in one direction, especially, where the step in other direction is too small; and that is the case of our model.
Of course, one could discuss the Alternating Direction Implicit (ADI) method [4,7] which is absolutely stable with relation to the choice of spatial step. The motivation for our choice of method derives from the relative simplicity of technical realization when compared to ADI and the natural affinity for parallel computing (parallelization in direction when it has conditional stability).
The numerical algorithm, described above, we implemented by means of the OpenCL language, and this realization is based on the following idea. In each time step, the cycle for j-index from 1 to M — 1 is parallelized. Each called thread simultaneously calculates the sought-after function by the Thomas algorithm, see Fig. 2. In the figure, we show the discretization of the function domain. Particularly, we group a set of points corresponding to one jth thread. We also show the points involved in the calculation of the given (i,j) point (bold point and crossed points on the Fig. 2).
1 Actually, the order of difference approximation depends on the choice of a norm [4,5].
Figure 2. Schematic representation of discretization of the function domain
2. Results of Numerical Simulations and Discussion
In this work we discuss the results of the numerical simulations only for one particular choice of the configuration of the object (Fig. 1). Its geometrical characteristics have been taken as follows: r0 = 0.12 cm, n = 0.125 cm, r2 = 0.13 cm, rmax = 0.1301cm, z0 = 4 cm, zmax = 5 cm. In the frame of this work, physical and engineering needs of the object geometry are not discussed.
The temperature dependencies of thermal coefficients, cym specific heat capacity and Xm thermal conductivity, for each materials are given as it is in [1,2]. For the chosen materials the corresponding data points have been taken from [8]. The densities of the chosen materials are p0 = 8.92, pi = p2 = 2 and p3 = 2.5 in units g/cm3. For this concrete configuration, the value 10-5 s is suitable for time-step t (9)). The period of source switching is tpvd = 25 ms, where the heating time is tsvc = 1 ms. The critical value of temperature is taken as Tcrit = 42 K (temperature of evaporation of working gases).
The initial temperature has been taken to be equal T0 = 4.2 K (the temperature of liquid helium). Note that because of the structural features of the object, especially because of existence of tiny layers covering the core cylinder, the choice of spatial step in the radial direction has to be smaller in comparison to the size of the layers (at least in order of magnitude) to guarantee the stability of the solution. Therefore, our choice of the difference scheme (see Section 1), which is suitable for the technical realization is justified.
40'
30.
£20. I—
10. 0
0.1
0.05
r [cm] 0 0
z [cm]
a At t = 1 ms after switching on the heat source at the very beginning
40
30
£20 I—
10 0
0.1 0.05 r [cm]
z [cm]
25 ms before switching on the
b At t
heat source of the next period
5
5
Figure 3. Temperature profile T(r, z) at different times
In Fig. 3 we show the temperature profiles at the very beginning. The temperature inside the object at t = 1 ms is shown in the left (Fig. 3,a). At the same moment, the source is switching off in the first period of source function. Because the radius of the cylinder is much smaller than its length, the heat first flows to the central axis and then to the right border, where the cryostat (the liquid helium temperature terminal) is located at zmax = 5 cm.
I- 20
10 -
ol-1-1-1-
0 5 10 15 20
t [s]
Figure 4. Time evolution of the surface temperature at z = 0 (solid line), critical temperature for evaporator and condensation of working gas (dashed line)
In the Fig. 3,b we show the temperature distribution at t = 25 ms (after first period). One period of the source switching Tpvd is enough time to equilibrate the temperature in radial direction (see Fig. 3,b). Such behavior repeats for approximately 5 ~ 10 sec (setting time) until a stable periodic regime is achieved (see Fig. 4). The setting time varies depending on the design of the object and switching time of the source.
The time calculations for different N x M are given in the Table 1 (here and in
the table: N = max Nj and M = max MA To demonstrate the results of OpenCL Vj J Vi
algorithm, the calculations for temperature evolution in time interval t E [0, 0.0765] with time-step t = 10-6 s have been carried out. In the table we use the following notations: CPU - Core i7 3517U (Ubuntu 11.0) and GPU - GF GTX 470 (Core
2 Duo, Debian 6.0). During the compilation of programs -O2 optimization flag has been used. It is shown that there is an interval of increasing number of points of discretization in axial direction (M) where the calculation time using GPU remains the same. In the Table 1 we compare the calculation times for CPU and GPU using
3 different grids with the same N = 431 and two with the same M = 401. Therefore, the choice of our algorithm allows us increase the density of our computational grid in the axial direction, practically without loss of calculation time.
Table 1
Calculation time (in seconds) of OpenCL implementation for different grid size
N x M T cpu T GPU TCPU/TGPU
431 x 101 309 345 0.896
431 x 201 607 349 1.739
431 x 401 1315 357 3.683
631 x 401 1968 509 3.866
3. Conclusion
We have suggested a computational scheme for solving heat conduction problem with periodic right side (source) in multilayer cylindrical domain. The scheme is based on explicit-implicit method. Because of the complex srtucure of the domain, the choice of the method is suitable for the technical realization. Its OpenCL realization has been developed. It is shown that there is an interval of increasing number of points of discretization in axial direction, where the calculation time using GPU remains the same.
The scheme has been applied for simulation of heat conduction process with a periodical source in criogenic cell. The results show that the temperature regime on the surface of the cryogenic cell has non-negligible setting mode (around 5 ~ 10 seconds) which is larger in comparison to the working period of the device. The key
characteristics of the device is working thermal process that has been achieved by a
particular choice of the model parameters.
References
1. A. Ayriyan, J. Pribis, Mathematical Simulation of Heat Conductivity in Composite Object with Cylindrical Symmetry, Matem. model. 24 (12) (2012) 113-118, in Russian.
2. A. Ayriyan, et al., Numerical Simulation of Heat Conductivity in Composite Object with Cylindrical Symmetry, Lecture Notes in Computer Science 7125 (2012) 264269.
3. D. E. Donets, et al., Physics Research and Technology Developments of Electron String Ion Sources, Review of Scientific Instruments 83 (2), art 02A512.
4. A. A. Samarskii, P. N. Vabishchevich, Computational Heat Transfer, Vol. 1, John Wiley & Sons Ltd., Chichester, England, 1995.
5. A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker Inc., New York, USA, 2001.
6. N. N. Yanenko, Fractional Step Methods for Solution of Multidimensional Problems of Mathematical Physics, Nauka, Novosibirsk, 1967, in Russian.
7. W. H. Press, et al., Numerical Recipes, 3rd Edition, Vol. 1, Cambridge University Press, New York, USA, 2007.
8. National Institute of Standards and Technology. URL http://www.nist.com
УДК 519.633.6, 519.688 Вычислительная схема решения задачи теплопроводности в многослойной цилиндрической области
А. С. Айриян*, Я. Буша мл.*, Е. Е. Донец*, О. Григорян*§,
Я. Прибиш*
* Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри 6, Дубна, Московская область, Россия, 141980 ^ Кафедра математики и теоретической информатики Технический университет г. Кошице ул. Летна 6, Кошице, Словакия, 04001 * Лаборатория физики высоких энергий Объединённый институт ядерных исследований ул. Жолио-Кюри 6, Дубна, Московская область, Россия, 141980 § Кафедра теоретической физики Ереванский государственный университет ул. Алека Манукяна 1, Ереван, Армения, 0025
В работе предложена вычислительная схема решения задачи теплопроводности с периодической функцией источника в многослойной цилиндрической области. Область имеет нетривиальную геометрию и теплофизические коэффициенты являются нелинейными функциями, зависящими от температуры, и имеют разрывы первого рода на границах слоёв. Вычислительная схема основана на алгоритме решения разностной задачи с явно-неявной пространственно-временной аппроксимацией. Обсуждается реализация алгоритма на языке OpenCL для проведения расчётов на ГПУ, дано сравнение с вычислениями, выполненными на ЦПУ. Показано, что схема может быть успешно применена для моделирования тепловых процессов в импульсной криогенной камере, предназначенной для импульсной подачи рабочих газов в рабочую область ионного источника в миллисекундном диапазоне. Результаты приведены для моделирования конкретной структуры криогенной камеры, которая удовлетворяет требованиям для его практической реализации.
Ключевые слова: уравнение теплопроводности, периодический источник, многослойная цилиндрическая структура, конечно-разностная схема, OpenCL.
Литература
1. Айриян А. С., Прибиш Я. Моделирование процесса теплопроводности в составном образце с цилиндрической симметрией // Математическое моделирование. — 2012. — Т. 24, № 12. — С. 113-118.
2. Ayriyan A. et al. Numerical Simulation of Heat Conductivity in Composite Object with Cylindrical Symmetry // Lecture Notes in Computer Science. — 2012. — Vol. 7125. — Pp. 264-269.
3. Donets D. E. et al. Physics Research and Technology Developments of Electron String Ion Sources // Review of Scientific Instruments. — 2012. — Vol. 83, No 2. — Art 02A512.
4. Samarskii A. A., Vabishchevich P. N. Computational Heat Transfer. — Chichester, England: John Wiley & Sons Ltd., 1995. — Vol. 1.
5. Samarskii A. A. The Theory of Difference Schemes. — New York, USA: Marcel Dekker Inc., 2001.
6. Яненко Н. Н. Метод дробных шагов решения многомерных задач математической физики. — Новосибирск: Наука, 1967.
7. Press W. H. et al. Numerical Recipes. — 3rd edition. — New York, USA: Cambridge University Press, 2007. — Vol. 1.
8. National Institute of Standards and Technology. — http://www.nist.com.
