An Efficient Algorithm for the Numerical Solution of a Three-dimensional Thermal Conductivity Issue

Full Text

ВВЕДЕНИЕ

В производстве листового проката одной из важнейших технологических операций является нагрев слябов под прокатку. Слябы в листопрокатном производстве представляют собой металлические заготовки толщиной 0,1–0,3 м и шириной 0,5–1,5 м, разрезаемые на выходе из установки непрерывной разливки стали на куски длиной 6–12 м. В нагревательную печь слябы укладываются поперек печи, размером сляба по вертикали при этом является его толщина. В зависимости от конструкции печи, теплота от сжигаемого топлива подается либо только к верхней поверхности сляба, либо к верхней и нижней, либо еще и к передней и задней граням сляба (подводом теплоты к торцевым поверхностям сляба обычно пренебрегают ввиду малой роли этого фактора в тепловом балансе сляба). Экспериментальное изучение этого процесса в промышленных условиях затруднено, поэтому обычно применяется только на начальной и завершающей стадиях исследования, а наиболее плодотворным методом изучения процесса нагрева сляба является математическое моделирование [1].

1. МЕТОДОЛОГИЧЕСКАЯ БАЗА ПРОВЕДЕННОГО ИССЛЕДОВАНИЯ

В зависимости от цели исследования, успешно используются одномерные, двухмерные и трехмерные модели нагрева сляба с различными граничными условиями [2; 3].

Одномерные модели используют в случаях, когда можно пренебречь неоднородностью температурного поля не только вдоль сляба, но и по его ширине, полагая, что температура меняется только по одной координате (толщине). В допущении о постоянстве теплофизических характеристик, однородности начального условия и об одностороннем характере нагрева, а также чисто конвективном характере нагрева на поверхности постановка задачи теплопроводности состоит из линейного однородного уравнения теплопроводности

$\frac{дТ}{дt}=\alpha \frac{{д}^2Т}{д{х}^2}, 0<х<\delta ;$ (1)

начального условия

T(x, 0) = T_н; (2)

граничного условия адиабатности на нижней поверхности

${\overline{)\frac{дТ}{дх}}}_{х=0}=0;$ (3)

и линейного граничного условия на тепловоспринимающей поверхности

$\lambda {\overline{)\frac{дТ}{дх}}}_{х=\delta }=\alpha ({Т}_0-{Т}_{х=\delta }).$ (4)

В формулах (1)–(4) использованы следующие обозначения: δ – толщина сляба, м; x – координата по толщине, м; t – время от начала процесса, с; T(x, t) – тем- пература, К; T_н – начальная температура, К; T₀ – температура греющей среды, К; a – коэффициент температуропроводности, м²/с; λ – коэффициент теплопроводности, Вт/м/К; α – коэффициент теплоотдачи, Вт/м²/К.

Для численного решения задачи теплопроводности (1)–(4) методом конечных разностей в расчетной области (0 ≤ x ≤ δ) вводится пространственная сетка с шагом Δx = δ/n, узлы которой имеют координаты x_i = (i – 1)Δx, (i = 1, 2, … , n + 1), где n – количество разбиений по координате. Вводится также дискретное время t_k = kΔt, (k = 1, 2, ...) с постоянным шагом Δt. Сеточные значения температуры тела в узловых точках x_i в моменты времени t_{k – 1} (в начале k-го шага) и t_k (в конце k-го шага) обозначают как T_i^{k – 1} и T_i^k, соответственно.

Вокруг каждого из узлов расчетной сетки (кроме крайних) на расстоянии Δx/2 проведем границы контрольных объемов. В итоге вся расчетная область окажется разделенной на непересекающиеся контрольные объемы, внутри каждого из которых располагается по одной узловой точке. При этом внутренние контрольные объемы (i = 2, … n) будут иметь толщину Δx и узловую точку в центре объема, а приграничные объемы (содержащие 1-ю и (n + 1)-ю узловые точки) – толщину Δx/2 и узловую точку, расположенную на краю объема, как это показано на рис. 1. В дальнейшем делается важное допущение о том, что температура внутри каждого из элементарных объемов в каждый момент времени распределена равномерно и равна температуре в представляющей его узловой точке.

 

Рис. 1. Разбиение расчетной области на контрольные объемы

Fig. 1. Splitting the calculation area into control volumes

 

В рамках наиболее физичной разновидности метода конечных разностей – метода контрольного объема (иногда его называют методом баланса) – разностные уравнения представляют собой уравнения элементарного теплового баланса, записанные для каждого элементарного объема на элементарном промежутке (шаге по времени) Δt. Эти уравнения получаются путем приравнивания изменения энтальпии элементарного объема за шаг по времени к разности между количествами теплоты, вошедшим в этот объем и вышедшим из него за то же время [4].

Опуская промежуточные выкладки, итоговое балансовое уравнение для 1-го контрольного объема (прилегающего к адиабатной поверхности) можно представить в виде

$(1+2\mu f)T_1^k-2\mu f{T}_2^k=2(1-\mu )f{T}_2^{k-1}+\left[1-2(1-\mu )f\right]T_1^{k-1},$ (5)

для внутренних объемов (i = 2, 3, ... n) балансовое уравнение может быть представлено как

$-\mu f{T}_{i-1}^k+(1+2\mu f)T_1^k-\mu f{T}_{i+1}^k=(1-\mu )f{T}_{i-1}^{k-1}+\left[1-2(1-\mu )f\right]T_i^{k-1}+(1-\mu )f{T}_{i+1}^{k-1},$ (6)

а для объема, примыкающего к поверхности теплообмена (i = n + 1) как

$-2\mu f{T}_n^k+\left[1+2\mu f(b+1)\right]T_{n+1}^k=2(1-\mu )f{T}_n^{k-1}+\left[1-2(1-\mu )f(b+1)\right]T_{n+1}^{k-1}+2fb{T}_0.$ (7)

В формулах (5)–(7) использованы следующие обозначения: f = aΔt/Δx² – сеточное число Фурье; b = αΔx/λ – сеточное число Био; μ – степень неявности разностной схемы (0 ≤ μ ≤ 1).

Полученная СЛУ (система линейных уравнений) (5)–(7) имеет трехдиагональную матрицу коэффициентов (в которой отличны от 0 только элементы, стоящие на главной диагонали и двух примыкающих к ней диагоналях [5]).

Для решения таких систем наиболее эффективен метод прогонки, представляющий собой результат применения алгоритма Гаусса к СЛУ с трехдиагональной матрицей коэффициентов. Достоинством этого метода является то, что при его использовании не требуется выделения памяти для хранения матрицы коэффициентов СЛУ (пусть даже треугольного вида), а достаточно двух вспомогательных одномерных массивов прогоночных коэффициентов. Поэтому этот алгоритм является очень быстрым и эффективным [4; 5].

Ситуация кардинальным образом меняется при рассмотрении двумерных и трехмерных задач. Например, если рассматривается трехмерная задача теплопроводности с однородными начальными условиями и несимметричными граничными условиями третьего рода, то при необходимости учета зависимости теплофизических характеристик от температуры постановка задачи состоит из нелинейного уравнения теплопроводности:

$\frac{дТ}{дТ}=\frac{1}{pc\left(T\right)}\left[\frac{д}{дх}\left(\lambda \left(Т\right)\frac{дТ}{дх}\right)+\frac{д}{ду}\left(\lambda \left(Т\right)\frac{дТ}{ду}\right)+\frac{д}{дz}\left(\lambda \left(Т\right)\frac{дТ}{дz}\right)\right], \left\{0<x<{\delta }_x, 0<у<{\delta }_{у}, 0<z<{\delta }_z\right\},$ (8)

начального условия:

T(x, y, z, 0) = T_н, (9)

и шести однотипных граничных условий:

$$\begin{gathered} {\overline{)-\left(\lambda \frac{дТ}{дх}\right)}}_{х=0}={\alpha }_{х}^{\uparrow }(T_0-T_{x=0}), \\ {\overline{)\lambda \frac{дТ}{дх}}}_{x={\delta }_x}={\alpha }_{х}^{\downarrow }(T_0-T_{x={\delta }_x}); \end{gathered}$$ (10)

$$\begin{gathered} {\overline{)-\left(\lambda \frac{дТ}{дy}\right)}}_{y=0}={\alpha }_y^{\uparrow }(T_0-T_{y=0}), \\ {\overline{)\lambda \frac{дТ}{дy}}}_{y={\delta }_y}={\alpha }_y^{\downarrow }(T_0-T_{y={\delta }_y}); \end{gathered}$$ (11)

$$\begin{gathered} {\overline{)-\left(\lambda \frac{дТ}{дz}\right)}}_{z=0}={\alpha }_z^{\uparrow }(T_0-T_{z=0}), \\ {\overline{)\lambda \frac{дТ}{дz}}}_{z={\delta }_z}={\alpha }_z^{\downarrow }(T_0-T_{z={\delta }_z}); \end{gathered}$$ (12)

В выражениях (8)–(12), которые «образуют замкнутую постановку дифференциальной задачи теплопроводности, использованы следующие обозначения: T(x, y, z, t) – температура заготовки (сляба), К; δ_x, δ_y, δ_z – размеры заготовки по соответствующим координатным направлениям (ширина, толщина, длина), м; ρ – плотность материала заготовки, кг/м³; с(T) – удельная теплоемкость материала заготовки, Дж/(кг ∙ К); λ(T) – коэффициент теплопроводности материала заготовки, Вт/(м ∙ К); α^↑_x и α^↓_x – коэффициенты теплоотдачи на задней и передней вертикальных поверхностях сляба, Вт/(м² ∙ К); α^↑_y и α^↓_y – коэффициенты теплоотдачи на нижней и верхней поверхностях сляба, Вт/(м² ∙ К); α^↑_z и α^↓_z – коэффициент теплоотдачи на левой и правой торцевых поверхностях сляба, Вт/(м² ∙ К); T₀ – температура греющей среды, К» [8].

Используя для численного решения метод баланса, вводим дискретное время t_k = kΔt, (k = 1, 2, ...) с постоянным шагом Δt и дискретные координаты x_i = iΔx, (i = 0, 1, 2, … , n_x), y_j = jΔy, (j = 0, 1, 2, … , n_y), z_l = lΔz, (l = 0, 1, 2, … , n_z) тоже для простоты с постоянными шагами Δx, Δy и Δz (n_x, n_y и n_z – количества разбиений заготовки вдоль каждого из координатных направлений).

«В результате вся расчетная область разбивается на элементарные объемы, количество которых N = (n_x + 1)(n_y + 1)(n_z + 1). Каждый из этих объемов содержит один узел пространственной сетки, которая задается трехиндексной нумерацией (i, j, l)» [8]. Разностные аналоги уравнения теплопроводности представляют собой балансовые уравнения для внутренних объемов расчетной области, а разностные аналоги граничных условий – балансовые уравнения для объемов, граница (хотя бы одна) которых совпадает с границей расчетной области. Уравнение для внутренних объемов (0 < i < n_x; 0 < j < n_y, 0 < l < n_z) имеет вид

$T_{i,j,l}^k-T_{i,j,l}^{k-1}=f_x^{-}({\overline{T}}_{i-1,j,l}-{\overline{T}}_{i,j,l})+f_x^{+}({\overline{T}}_{i+1,j,l}-{\overline{T}}_{i,j,l})+f_y^{-}({\overline{T}}_{i,j-1,l}-{\overline{T}}_{i,j,l})+f_y^{+}({\overline{T}}_{i,j+1,l}-{\overline{T}}_{i,j,l})+f_z^{-}({\overline{T}}_{i,j,l-1}-{\overline{T}}_{i,j,l})+f_z^{+}({\overline{T}}_{i,j,l+1}-{\overline{T}}_{i,j,l}) ,$ (13)

где величины

$$\begin{gathered} f_x^{-}=\frac{{\lambda }_{i-0,5,j,l}}{{\left(pc\right)}_{i,j,l}}\frac{∆t}{∆x^2} ; f_y^{-}=\frac{{\lambda }_{i,j-0,5,l}}{{\left(pc\right)}_{i,j,l}}\frac{∆t}{∆y^2} ; \\ {f}_z^{-}=\frac{{\lambda }_{i,j,l-0,5}}{{\left(pc\right)}_{i,j,l}}\frac{∆t}{∆z^2} ; \end{gathered}$$

представляют собой левые, а величины

$$\begin{gathered} f_x^{+}=\frac{{\lambda }_{i+0,5,j,l}}{{\left(pc\right)}_{i,j,l}}\frac{∆t}{∆x^2} ; f_y^{+}=\frac{{\lambda }_{i,j+0,5,l}}{{\left(pc\right)}_{i,j,l}}\frac{∆t}{∆y^2} ; \\ {f}_z^{+}=\frac{{\lambda }_{i,j,l+0,5}}{{\left(pc\right)}_{i,j,l}}\frac{∆t}{∆z^2} ; \end{gathered}$$

– правые сеточные числа Фурье для каждого из координатных направлений [4], а в правой части стоят определяющие температуры, представляющие собой в общем случае линейную комбинацию температур в начале и конце шага по времени, в которой весовым коэффициентом является степень неявности разностной схемы μ:

${\overline{T}}_{i,j,k}=\mu T_{i,j,k}^k+(1-\mu )T_{i,j,k}^{k-1} .$ (14)

Важно определить. какое количество неизвестных содержит каждое из уравнений вида (13). Если используется явная разностная схема (μ = 0), то все определяющие температуры, стоящие в правой части уравнения (13), превращаются в узловые значения в начале шага по времени, и неизвестным останется только одно значение T^k_{i, j, l}, стоящее в левой части этого уравнения. При этом СЛУ превращается в набор расчетных формул, не требующих никаких вспомогательных массивов для решения. Однако явная разностная схема дает сходящееся решение только при выполнении определенных ограничений, накладываемых на сеточные числа Фурье. Фактически это ограничения на шаг по времени, зависящие от шагов разбиения по координатам [5]. При наличии зависимости теплофизических характеристик и параметров граничных условий от температуры бывает трудно проконтролировать соблюдения упомянутых ограничений (условий сходимости), поэтому явную разностную схему стараются не использовать [4–6].

В общем случае применения неявных разностных схем (μ ≠ 0) каждое уравнение содержит семь неизвестных (для трехмерных задач), при этом матрица коэффициентов такой СЛУ оказывается очень разреженной и перестает быть трехдиагональной. Это объясняется тем, что каждый объем граничит только с двумя соседними объемами в каждом из координатных направлений. При этом объемы, являющиеся соседними в расчетной области, не всегда имеют близкие номера в сквозной нумерации узлов. В итоге приходится решать СЛУ общего вида, что для таких гигантских разреженных матриц нецелесообразно, поскольку требует очень больших вспомогательных массивов для хранения обратной матрицы СЛУ размером N × N.

Разумной альтернативой представляется метод простой итерации, позволяющий в каждом уравнении трактовать как неизвестную величину только температуру в конце шага по времени в том узле, для которого записано балансовое уравнение, а температуры в конце шага по времени в остальных узлах брать с предыдущей итерации. Этот подход тоже позволяет избежать применения вспомогательных матриц (как и явная схема), однако требует многократных повторов вычисления (итераций) для достижения заданной точности вычислений. Недостатком этого метода является медленная сходимость, а также возможная потеря устойчивости при сильной зависимости теплофизических характеристик от температуры [7].

2. ПРЕДЛАГАЕМЫЙ ПОДХОД К ВЫБОРУ АЛГОРИТМА РЕШЕНИЯ СИСТЕМЫ РАЗНОСТНЫХ УРАВНЕНИЙ МАТЕМАТИЧЕСКОЙ МОДЕЛИ НАГРЕВА СЛЯБА

Целью данной работы является разработка эффективного алгоритма решения системы разностных уравнений математической модели нагрева сляба, позволяющего с требуемой точностью решать систему разностных уравнений трехмерной задачи теплопроводности без чрезмерных затрат времени и компьютерных ресурсов.

Идея предлагаемого алгоритма основана на том, что, поскольку перенос теплоты по толщине сляба является основным направлением (именно на нем сосредоточены одномерные модели), то именно описывающие этот перенос члены балансового уравнения следует трактовать как неизвестные величины, тогда как для остальных направлений значения неизвестных температур могут быть взяты с предыдущей итерации [8]:

$-\mu f_y^{-}{T}_{i,j-1,l}^k+(1+\mu f_{\sum })T_{i,j,l}^k-\mu f_y^{+}{T}_{i,j+1,l}^k=T_{i,j,l}^{k-1}\left[1-(1-\mu )f_{\sum }\right]+f_x^{-}{\overline{T}}_{i-1,j,l}+f_x^{+}{\overline{T}}_{i+1,j,l}+(1-\mu )(f_y^{-}{T}_{i,j-1,l}^{k-1}+f_y^{+}{T}_{i,j+1,l}^{k-1})+f_z^{-}{\overline{T}}_{i,j,l-1}+f_z^{+}{\overline{T}}_{i,j,l+1},$ (15)

где для сокращения записи использовано обозначение:

$f_{\sum }=f_x^{-}+f_x^{+}+f_y^{-}+f_y^{+}+f_z^{-}+f_z^{+}.$

Уравнения вида (15) формально имеют трех-диагональную матрицу коэффициентов, и поэтому они могут быть решены методом прогонки, однако этот метод должен быть применен в рамках итерационной процедуры, поскольку определяющие температуры T̅_{i –1, j, l}, T̅_{i +1, j, l}, T̅_{i, j, l –1}, T̅_{i, j, l + 1}, стоящие в правой части, требуют уточнения. Такая итерационная процедура сходится очень быстро, поскольку распространение теплоты по ширине и длине сляба, учитываемое ею, является существенно менее выраженным, чем явным образом учитываемое в алгоритме прогонки распространение ее по толщине сляба. В дальнейшем изложении будем называть предлагаемый метод послойным.

Предложенный подход программно реализован в среде Builder C++ 6.0, стартовое окно расчетной программы представлено на рис. 2.

 

Рис. 2. Стартовое окно расчетной программы

Fig. 2. The start window of the calculation program

 

Как видно из представленного рисунка, расчетная программа предполагает задание размеров сляба, плотности материала сляба (как константы), коэффициента теплопроводности λ и удельной теплоемкости c (как полиномов второй степени), начальной температуры сляба, продолжительности нагрева, граничных условий для каждой грани сляба, а также параметров разностной схемы (количество разбиений по каждому из направлений, «шаг по времени Δt, степень неявности разностной схемы, допустимая погрешность расчета температуры») [8]. Кроме того, программа позволяет выбрать один из методов решения – послойный метод или метод простой итерации.

С помощью этого инструмента были проведены вариантные расчеты при одних и тех же исходных данных, но при разном количестве разбиений. Результаты таких расчетов оказались идентичными (расхождение между значениями температур, полученных разными методами, не превышало заданного требуемого уровня точности), однако продолжительность расчета, существенно зависящая от количества узлов расчетной сетки, для двух методов решения СЛУ оказалась различной, как это видно из табл. 1.

Поскольку абсолютное значение продолжительности расчета зависит от мощности компьютера, интерес представляют не абсолютные, а относительные значения продолжительности расчета, представленные в последней колонке табл. 1.

В результате сравнения установлено, что с измельчением пространственной сетки эффективность послойного метода по сравнению с методом простой итерации увеличивается.

 

Таблица 1

Сравнение продолжительности расчета послойным методом и методом простой итерации при различных параметрах расчетной сетки [Comparison of the calculation duration by the layer-by-layer method and the method of simple iteration with different parameters of the calculation grid]

Вариант [Option]	Количество разбиений [Number of splits]			Δt, c	Время счета, ч:м:с [Counting time (hours, minutes, second)]		Эффективность (отношение) [Efficiency (ratio)]
	nx	ny	nz		Простая итерация [Simple iteration]	Послойный метод [Layer-by-layer method]
1	10	10	40	10	0:00:17	0:00:11	1,54
2	20	20	80	10	0:03:32	0:01:39	2,14
3	30	30	120	10	0:18:07	0:07:12	2,52
4	40	40	160	10	1:01:55	0:22:35	2,74

 

ЗАКЛЮЧЕНИЕ

В работе разработана и программно реализована в среде Builder C++ 6.0 математическая модель нагрева сляба, представляющая собой трехмерную нестационарную задачу теплопроводности с зависящими от температуры теплофизическими характеристиками, с возможностью выбора метода решения системы разностных уравнений (метод простой итерации или предложенный послойный метод).

Также в работе изучена зависимость времени расчета от шагов разбиения пространственной сетки. В результате было установлено, что с измельчением сетки эффективность предложенного послойного метода возрастает по сравнению с методом простой итерации.

About the authors

Aleksander V. Vargin

National Research Technological University “MISIS”

Author for correspondence.
Email: mr.vargin@yandex.ru
ORCID iD: 0000-0001-8657-8716

PhD student

Russian Federation, Moscow

Igor A. Levitskiy

National Research Technological University “MISIS”

Email: lewwwis@mail.ru
ORCID iD: 0000-0002-9345-3628

Cand. Sci. (Eng.), Associate Professor

Russian Federation, Moscow

References

Supplementary files