Method of additional boundary conditions in the problem of heat transfer for non-newtonian fluid moving in laminar mode in circular pipe

Full Text

Введение. В ряде современных технических устройств для эффективного отвода тепла используются элементы конструкций с каналами, по которым движется жидкий теплоноситель: системы охлаждения ядерных реакторов [1], зеркал технологических лазеров [2], лопаток современных газотурбинных установок [3] и др. С точки зрения технологической реализации таких элементов, по-видимому, наиболее целесообразно каналы в них создавать за счёт системы трубок, упакованных в связующий материал, так как трубки при этом выполняют как роль каналов, так и роль армирующих элементов, которые можно укладывать по достаточно сложным траекториям [4, 5]. Как показано в работах [6, 7], для математического описания процессов теплопереноса в таких трубчато-армированных конструкциях предварительно нужно решить задачу о теплопереносе в отдельно взятой тепловой трубке, причём это решение должно быть достаточно простым по структуре, чтобы оно было удобным для дальнейшего анализа проблемы в целом и при этом достаточно точным с точки зрения последующего инженерного приложения. Этим условиям удовлетворяет метод дополнительных граничных условий, предложенный в [8] для приближённого решения задачи Гретца-Нуссельта и использованный в [6, 7] для определения температурного поля в отдельно взятой тепловой трубке при стационарном и нестационарном теплопереносе. Однако решения, полученные в [6-8], основывались на предположении о том, что по трубе в стабилизированном ламинарном режиме протекает жидкость, распределение скорости течения которой в поперечном сечении подчиняется закону Пуазеля-Гагена и диссипация механической энергии которой не учитывается. В действительности же в качестве жидкого теплоносителя могут быть использованы достаточно вязкие вещества (например, масла [9]), профиль скорости которых в круглой трубе при стабилизированном ламинарном режиме течения не соответствует закону Пуазеля-Гагена, т. е. такой теплоноситель представляет собой неньютоновскую жидкость. Кроме того, при исследовании тепловых процессов в таком теплоносителе целесообразно учитывать и диссипацию механической энергии, вызванную неравномерностью профиля скорости по сечению трубы [10]. В связи с этим настоящее исследование посвящено построению приближённого решения задачи теплопереноса для неньютоновской жидкости, прокачиваемой в стабилизированном ламинарном режиме по круглой трубе, методом дополнительных граничных условий с учётом диссипации механической энергии в теплоносителе. Из сравнения решений, приведённых в [6, 7] для ньютоновской жидкости, прокачиваемой по тепловой трубке, можно заключить, что в случае стационарного теплообмена методом дополнительных граничных условий можно получить более точное решение, чем в нестационарной задаче, без увеличения порядка разрешающего дифференциального уравнения теплопереноса внутри трубки. В силу этого в настоящем исследовании также раздельно изучаются стационарный и нестационарный случаи. 1. Нестационарный случай. Рассматривается стабилизированный ламинарный режим течения несжимаемой нетермочувствительной неньютоновской жидкости в круглой трубе, профиль скорости которой подчиняется степенному закону Освальда-де Виля. При этом уравнение теплового баланса жидкости с учётом диссипации механической энергии в цилиндрической системе координат в предположении об осевой симметрии задачи имеет вид [10] ∂T ∂T ∂2T 1 ∂T ∂2T k dv + v(r) =a + + + ∂t ∂z ∂z 2 r ∂r ∂r2 ρcp dr 0 r r0 , 0 z L, t 0, n+1 , (1) 579 Я н к о в с к и й А. П. где a= λ = const, ρcp r r0 v(r) = v0 1 - (n+1)/n ; (2) λ, ρ, cp - коэффициент теплопроводности, плотность и удельная теплоемкость (при постоянстве давления) жидкости; k, n - реологические постоянные жидкости: k - постоянная консистентности (k 0), n - индекс течения (0 < n 1); v - скорость жидкости в продольном (осевом) направлении z; v0 - скорость жидкости на оси трубы (при r = 0); r - полярный радиус; r0 - радиус внутренней стенки трубы; L - длина трубы; T - температура жидкости; t - время. При n = 1 согласно (2) получаем профиль течения жидкости, соответствующий закону Пуазеля-Гагена. Уравнение (1) с учётом (2) целесообразно переписать так: ∂2T 1 ∂T v(r) ∂T r + = + D(T ) - W 2 ∂r r ∂r a ∂z r0 0 r r0 , 0 z L, t 0, (n+1)/n , (3) где k n + 1 v0 n+1 1 ∂( · ) ∂ 2 ( · ) - ; (4) = const, D( · ) = λ n r0 a ∂t ∂z 2 D - параболический оператор одномерной задачи теплопроводности. Последние слагаемые в правых частях уравнений (1), (3) порождены диссипацией механической энергии в жидкости [10]. Для приближённого решения начально-краевой задачи, соответствующей уравнению (3), используем метод начальных функций, совмещённый с методом дополнительных граничных условий. Согласно идее метода начальных функций [11], разложим функцию T в ряд Маклорена по переменной r: W = ∞ T (t, z, r) = T0 (t, z) + m=1 rm ∂ m T m! ∂rm r=0 , 0 r r0 , (5) где T0 (t, z) ≡ T (t, z, 0) - температура жидкости в точках, лежащих на оси трубы («начальная» функция). Так как распределение температуры T в трубе считается осесимметричным, в (5) ∂mT lim = 0, m = 1, 3, 5, 7, . . . (6) r→0 ∂r m Ограничиваясь частичной суммой ряда (5), с учётом (6) получим M T (t, z, r) = T0 (t, z) + T2m (t, z) m=1 0 r r0 , 0 z L, t r r0 0, 2m , (7) где T0 , T2m - функции, подлежащие определению и зависящие только от времени t и от осевой координаты z. 580 Метод дополнительных граничных условий в задаче теплопереноса . . . Выразим в (7) все функции T2m через T0 . Для этого будем использовать уравнение (3) и его производные по r при r → 0. Так как левая часть уравнения (3) полностью совпадает с левой частью равенства (9) в работе [7], то, повторяя почти дословно рассуждения (14)-(26) из [7], с учётом (4), (6) получим T2 (t, z) = 2 r0 ¯ D (T0 ) , 4 T4 (t, z) = 4 r0 ¯ 2 D (T0 ) , 64 0 < n < 1, (8) где v0 ∂( · ) 1 ∂( · ) ∂ 2 ( · ) v0 ∂( · ) ¯ D( · ) = D( · ) + = - + . a ∂z a ∂t ∂z 2 a ∂z (9) Второе равенство (8) справедливо только при выполнении строгого неравенства n < 1. В случае же n = 1 выражение для T4 содержит еще одно слагаемое (см. формулу (26) в [7]). Так как при n = 1 профиль скорости соответствует закону Пуазеля-Гагена (стабилизированному ламинарному режиму течения ньютоновской жидкости), этот случай далее не рассматривается. Согласно (8) с учётом (9) функции T2 , T4 определяются через производные от начальной функции T0 (t, z). Аналогично можно получить выражения и для функций T2m (3 m M ) в разложении (7) через производные от T0 . На основании (8), (9) эти зависимости можно представить в таком виде: T2m (t, z) = (m) (T0 (t, z); W ), 1 m M, (10) где (m) - дифференциальный оператор порядка m по времени t и порядка 2m по пространственной переменной z, который при m 3 может зависеть и от величины W , характеризующей диссипацию механической энергии в жидкости. Подставив (10) в (7), получим выражение для температуры жидкости T (t, z, r) = (M ) (r; 0 0 r r0 , z T0 (t, z); W ), L, t (11) 0, где (M ) - линейный дифференциальный оператор порядка M по времени t и порядка 2M по пространственной переменной z от начальной функции T0 , коэффициенты которого зависят от переменной r и величины W . Для определения функции T0 необходимо использовать основное граничное условие на стенке трубы (при r = r0 ), например, граничное условие I рода: T (t, z, r0 ) = Tw (t, z), 0 z L, t 0, (12) где Tw - заданная температура стенки трубы. Подставив (11) в (12), получим разрешающее дифференциальное уравнение параболического типа порядка 2M : (M ) (r0 ; T0 (t, z); W ) = Tw (t, z), 0 z L, t 0, (13) где правая часть - известная функция. 581 Я н к о в с к и й А. П. Для однозначного интегрирования уравнения (13) для начальной функции T0 необходимо задать соответствующее количество начальных (при t = 0) и граничных условий на входе (z = 0) и выходе (z = L) жидкости из трубы. Таким образом строится приближённое решение уравнения (3) (или (1)) в рамках метода начальных функций. Повысить точность приближённого решения (7) без увеличения порядка разрешающего уравнения можно за счёт совмещения метода начальных функций с методом дополнительных граничных условий [8]. С этой целью откажемся от выражений (11) при m 3 и воспользуемся дополнительным граничным условием, которое получается из (3) с учётом (2), (4) при r = r0 : ∂2T 1 ∂T + 2 ∂r r ∂r r=r0 = D(Tw ) - W, 0 z L, t 0, (14) где в правой части использовано основное граничное условие (12). Подставив в (12), (14) разложение (7) и ограничившись рассмотрением случая M = 4, получим систему уравнений 4 T2m (t, z) = Tw (t, z), m=0 4 (15) r2 m T2m (t, z) = 0 D(Tw ) - W . 4 2 m=0 Из системы (15), являющейся системой линейных алгебраических уравнений (СЛАУ), с учётом (8) определим функции 1 r2 - 0 [D(Tw ) - W ] + 16(Tw - T0 )- 7 4 15 2 ¯ 3 4¯ - r0 D(T0 ) - r0 D2 (T0 ) , 4 16 2 1 r0 [D(Tw ) - W ] - 9(Tw - T0 )+ T8 (t, z) = 7 4 5 4¯ 2¯ +2r0 D(T0 ) + r0 D2 (T0 ) . 64 T6 (t, z) = (16) Таким образом, функции T6 , T8 , как и T2 , T4 (см. (8)), также удалось выразить через начальную функцию T0 . Для определения функции T0 используем интегральное уравнение теплового баланса жидкости в трубе [8], т. е. уравнение (3) проинтегрируем по r от 0 до r0 с весом r: v0 a r0 1- 0 r r0 (n+1)/n ∂T r dr+ ∂z r0 0 r0 = 0 582 r0 r r0 D(T )r dr - W + 0 ∂2T 1 ∂T + rdr = ∂r2 r ∂r r0 0 (n+1)/n r dr = ∂ ∂T ∂T r dr = r0 ∂r ∂r ∂r r=r0 , (17) Метод дополнительных граничных условий в задаче теплопереноса . . . где учтён предельный переход (6) при m = 1. Подставим в левую и правую части цепочки равенств (17) выражение (7) при M = 4, тогда с учётом (8), (9), (16) получим 6 6 v0 r0 r0 1 ∂ 2 D3 (T0 ) + D (T0 )+ a3 + 2688 2a 672 ∂z v2 r6 1 37 4 2 ∂2 + 0 2 0 a3 + D(T0 ) + r D (T0 )+ 2 a 2688 ∂z 1120 0 4 v0 r0 47 ∂ v 3 r6 a3 ∂ 3 T0 27 2 + a2 + D(T0 ) + 0 03 + r0 D(T0 )+ 2a 560 ∂z 2a ∂z 3 35 2 4 2 v0 r0 1 ∂ 2 T0 v0 r0 a1 3 ∂T0 + + + - a2 + 2a2 56 ∂z 2 a 2 7 ∂z 4 r4 v0 r0 a5 ∂ - 0 D2 (Tw ) - D(Tw )+ 1120 2a ∂z 2 3 2 v0 r0 a4 ∂Tw 24 + r0 D(Tw ) + - (Tw - T0 ) = 35 2a ∂z 7 1 n 2 - W, (18) = r0 3n + 1 14 где 24(n + 1)(127n2 + 22n + 1) a1 = , 35(3n + 1)(9n + 1)(11n + 1) 3(n + 1)(1415n2 + 216n + 9) a2 = , 560(5n + 1)(9n + 1)(11n + 1) (n + 1)(189n2 + 26n + 1) (19) a3 = , 1344(7n + 1)(9n + 1)(11n + 1) (n + 1)(139n + 11) , a4 = 35(9n + 1)(11n + 1) (n + 1)(19n + 1) a5 = , 0 < n < 1. 560(9n + 1)(11n + 1) Последнее слагаемое в левой части уравнения (18) содержит разность (Tw - T0 ), характеризующую величину температурного напора от стенки к жидкости. Согласно (4), разрешающее уравнение (18) с учётом (19) содержит высшие производные от T0 по времени t третьего порядка (∂ 3 T0 /∂t3 ) и шестого порядка по z (∂ 6 T0 /∂t6 ), поэтому для однозначного интегрирования уравнения (18) необходимо задать три начальных условия (для T0 , ∂T0 /∂t, ∂ 2 T0 /∂t2 при t = 0) и по три граничных условия на входе и выходе жидкости из трубы (для T0 , ∂T0 /∂z, ∂ 2 T0 /∂z 2 или их линейных комбинаций при z = 0 и z = L). Получим сначала все необходимые начальные условия для функции T0 . Пусть в момент времени t = 0 в трубе задано начальное распределение температуры T (0, z, r) = T 0 (z, r), 0 r r0 , 0 z L, (20) где T 0 - заданная функция. Так как в настоящем исследовании рассматривается лишь конечное разложение температуры T (см. формулу (7)), в общем случае равенство (20) 583 Я н к о в с к и й А. П. выполнить точно не удается, поэтому начальное условие (20) можно удовлетворить лишь приближённо, используя либо метод коллокаций [12], либо метод наименьших квадратов [13], либо метод взвешенных невязок [14]. При использовании последнего метода вместо равенства (20) применим, например, следующие интегральные соотношения: r0 T (0, z, r) 0 r r0 2m r dr = r0 = 0 T 0 (z, r) r r0 2m 0 r dr ≡ θm (z), m = 0, 1, 2, (21) 0 0 0 где θ0 , θ1 , θ3 - известные функции переменной z ∈ [0, L]. Подставим в (21) разложение (7) при M = 4 и учтём (8), (9), (16), тогда после приведения подобных слагаемых получим: 24 27 2 ¯ r4 ¯ T0 + r0 D(T0 ) + 0 D2 (T0 ) = 35 560 1344 2 0 = 2 θ0 (z) + r0 2 r0 ¯ 11 4 ¯ 2 9 T0 + D(T0 ) + r D (T0 ) = 35 42 26880 0 2 0 = 2 θ1 (z) + r0 4 20 11 2 ¯ r0 ¯ 2 T0 + r D(T0 ) + D (T0 ) = 147 784 0 3920 2 0 = 2 θ2 (z) + r0 2 r0 11 [D(Tw ) - W ] - Tw , 560 35 2 r0 17 [D(Tw ) - W ] - Tw , 840 70 (22) 2 r0 29 [D(Tw ) - W ] - Tw , 1176 147 где в правые части перенесены значения всех известных в момент времени t = 0 функций и величин. Из системы (22), как из СЛАУ, в каждой точке z ∈ [0, L] в начальный мо¯ ¯ мент времени можем определить значения функций T0 , D(T0 ), D2 (T0 ). Тогда с учётом (9) получаем следующие начальные условия: ¯ T0 (0, z) = T0 (z), ¯ ¯ ∂T0 ∂ 2 T0 ∂ T0 ¯ ¯ = aD(T0 ) + a 2 - v0 , ∂t t=0 ∂z ∂z t=0 2¯ ¯ ∂ 2 T0 ∂ 4 T0 2 ∂ T0 ¯ ¯ = a2 D2 (T0 ) - a2 - v0 + ∂t2 t=0 ∂z 4 ∂z 2 ¯ ∂ 2 ∂T0 ∂ ∂T0 ∂ 3 T0 +2a 2 - 2v0 + 2v0 a 3 ∂z ∂t ∂z ∂t ∂z (23) t=0 , ¯ ¯ ¯ ¯ ¯ где T0 , D(T0 ), D2 (T0 ) - значения функций, известных из решения системы (22) в каждой точке z ∈ [0, L]. Последовательно определяя из (23) T0 , ∂T0 /∂t, ∂ 2 T0 /∂t2 , получим все необходимые начальные условия для интегрирования уравнения (18) с учётом (19). (При таком последовательном вычислении правые части в (23) - известные функции переменной z.) 584 Метод дополнительных граничных условий в задаче теплопереноса . . . Определим теперь граничные условия для уравнения (18). Пусть на входе (z = z∗ = 0) и на выходе (z = z∗∗ = L) жидкости из трубы заданы, например, распределения температуры по поперечному сечению T (t, z∗ , r) = T∗ (t, r), T (t, z∗∗ , r) = T∗∗ (t, r), 0 r r0 , t 0, (24) где T∗ , T∗∗ - известные функции. По-прежнему в силу использования конечного разложения температуры (7) в общем случае граничные условия (24) выполнить не удается, поэтому, как и в случае начальных условий, равенства (24) можно удовлетворить лишь приближённо, используя, например, метод коллокаций, метод наименьших квадратов или метод взвешенных невязок [12-14]. В случае применения последнего метода вместо (24) используем интегральные равенства, аналогичные (21): r0 r r0 T (t, z∗ , r) 0 2m r dr = r0 = T∗ (t, r) 0 r r0 2m ∗ rdr ≡ θm (t), m = 0, 1, 2 (∗ → ∗∗), (25) ∗ ∗∗ где θm и θm (m = 0, 1, 2) - известные функции переменной t. Подставим в (25) разложение (7) при M = 4, тогда с учётом (8), (9), (16) получим систему, аналогичную (22), в которой следует сделать замены 0 ∗ θm (z) → θm (t) (z = z∗ ), 0 ∗∗ θm (z) → θm (t) (z = z∗∗ ), m = 0, 1, 2. (26) ∗ ∗∗ Так как функции θm (t) и θm (t) (m = 0, 1, 2) известны (в силу правых равенств в (25)), а зависимость температуры стенки Tw от t и z предполагается заданной (см. (12)), правые части в (22) с учётом замены (26) также известны. Следовательно, из (22) с учётом (26), как из СЛАУ, можем определить ¯ ¯ значения функций T0 , D(T0 ), D2 (T0 ) (в сечениях z = z∗ , z = z∗∗ ), зависящих только от времени t. При этом с учётом (9) в точке z = z∗ (или z = z∗∗ ) получаем требуемые граничные условия: ¯ T0 (t, z∗ ) = T0 (t), ¯ ∂T0 ∂ T0 ∂ 2 T0 ¯ ¯ - v0 = - aD(T0 ) , 2 ∂z ∂z ∂t z=z∗ 2 ∂ 4 T0 ∂ 3 T0 2 ∂ T0 a2 - 2v0 a 3 + v0 = ∂z 4 ∂z ∂z 2 2 ¯ ∂ T0 ∂ ¯ ¯ ¯ ¯ = a2 D2 (T0 ) + - 2a D(T0 ) 2 ∂t ∂t a (27) z=z∗ (z∗ → z∗∗ ), ¯ ¯ ¯ ¯ ¯ где функции T0 (t), D(T0 ), D2 (T0 ) известны из решения системы (22) с учётом замены (26) в точке z = z∗ (или z = z∗∗ ) в каждый момент времени t 0. Последовательно определяя в (27) правые части, получим в точках z = z∗ и z = z∗∗ по три граничных условия, необходимых для однозначного интегрирования разрешающего уравнения (18), если в сечениях z = z∗ и z = z∗∗ заданы граничные условия I рода для температуры жидкости (24). 585 Я н к о в с к и й А. П. Таким образом, исходя из начальных (20) и граничных (24) условий, заданных для температуры жидкости T , удалось определить начальные (23) и граничные (27) условия для начальной функции T0 , необходимые для однозначного интегрирования разрешающего уравнения (18). После решения начально-краевой задачи (18), (23), (27) определим начальную функцию T0 (t, z) (t 0, 0 z L), затем по формулам (8), (16) с учётом (9), (4) можем вычислить и функции T2m (t, z) (m = 1, 4), подставив которые в разложение (7), получим приближённое распределение температуры жидкости в трубе T (t, z, r), найденное совместным применением метода начальных функций и метода дополнительных граничных условий в нестационарном случае с учётом диффузионной составляющей теплопереноса в осевом направлении z. Уточнить полученное приближённое решение можно лишь за счёт использования равенства (10) при m 3, но в конечном итоге это приведёт к увеличению порядка разрешающего уравнения по сравнению с (18). 2. Стационарный случай. Предполагается, что в трубе установился стационарный режим теплопереноса. Кроме того, как это обычно делается [8, 10], будем пренебрегать диффузионной составляющей теплопереноса в жидкости вдоль трубы по сравнению с её конвективной составляющей. В этом случае по сравнению с ранее изложенным (см. п. 1) за счёт еще одного дополнительного граничного условия на стенке трубы (r = r0 ) удаётся уточнить решение рассматриваемой задачи без повышения порядка разрешающего дифференциального уравнения. Согласно сделанным предположениям, в уравнении (1) не учитываем слагаемые, содержащие производные ∂T /∂t и ∂ 2 T /∂z 2 , т. е. вместо (3) используем уравнение ∂2T 1 ∂T v(r) ∂T r + = -W ∂r2 r ∂r a ∂z r0 0 r r0 , z 0, (n+1)/n , (28) где T = T (z, r); обозначения W , a, v(r) даны в (2), (4). Приближённое решение уравнения (28) по-прежнему разыскивается в виде (5)-(7). При этом на основании (8), (9) с учётом сделанных выше предположений получаем1 T2 (z) = 2 r0 v0 dT0 , 4a dz T4 (z) = 4 2 r0 v0 d2 T0 , 64a2 dz 2 T0 = T0 (z). (29) Справедливыми остаются рассуждения (10)-(13), где все функции не зависят от переменной t. Для дальнейшего уточнения решения рассматриваемой задачи теплопереноса вновь помимо основного граничного условия (12) используем дополнительные граничные условия на стенке трубы. С этой целью рассмотрим уравнение (28) при r = r0 , тогда вместо (14) с учётом (2) получим ∂2T 1 ∂T + 2 ∂r r ∂r r=r0 = -W = const, z 0. (30) 1 Соотношения (29), как и (8), имеют место при выполнении неравенств 0 < n < 1; случай n = 1 по-прежнему не рассматривается. 586 Метод дополнительных граничных условий в задаче теплопереноса . . . Для получения второго дополнительного граничного условия продифференцируем (28) по r: ∂3T 1 ∂2T 1 ∂T + - 2 = ∂r3 r ∂r2 r ∂r v(r) ∂ ∂T = a ∂z ∂r - n + 1 v0 r an r0 r0 1/n ∂T ∂z - 1/n W n+1 r n r0 r0 . (31) Рассмотрим уравнение (31) при r = r0 . С учётом (2) получим следующее дополнительное граничное условие: ∂3T 1 ∂2T 1 ∂T + - 2 3 2 ∂r r ∂r r ∂r r=r0 =- n + 1 v0 dTw n+1W - , an r0 dz n r0 z 0, (32) где в правой части использовано основное граничное условие (12). Подставим в равенства (12), (30), (32) разложение (7), в котором в отличие от нестационарного случая (см. п. 1) примем M = 5. Тогда после элементарных преобразований с учётом (29) получим 2 r4 v 2 d2 T0 r0 v0 dT0 2 + c2m 0 20 + d2m r0 W + a dz a dz 2 n + 1 2 v0 dTw + e2m r + W + f2m Tw , m = 3, 4, 5, (33) n 0 a dz T2m (z) = a2m T0 + b2m где 225 72 200 , a8 = -f8 = , a10 = -f10 = - , 47 47 47 87 46 57 b6 = - , b8 = , b10 = - , 94 47 188 93 13 57 , c8 = , c10 = - , c6 = - 1504 3008 1504 13 41 15 d6 = , d8 = - , d10 = , 94 188 188 9 1 7 e6 = - , e8 = , e10 = - . 752 47 752 a6 = -f6 = - (34) Для определения функции T0 (z) по-прежнему используем интегральное уравнение теплового баланса жидкости в трубе, т. е. уравнение (17), в котором, согласно сделанным предположениям, следует принять D(T ) ≡ 0, тогда с учётом (33), (34), (29) вместо (18) получим A1 6 3 r0 v0 d3 T0 r4 v 2 d2 T0 r2 v0 dT0 + A2 0 2 0 + A3 0 + 240T0 + a3 dz 3 a dz 2 a dz r4 v 2 d2 Tw r2 v0 dTw 2 + A4 0 2 0 + A5 0 - 240Tw + A6 r0 W = 0, a dz 2 a dz z 0, (35) 587 Я н к о в с к и й А. П. где n+1 47 57 93 13 - + - , 64 3(7n + 1) 2(9n + 1) 5(11n + 1) 3(13n + 1) 1 47 87 368 57 = 3 + (n + 1) - + - 8 5n + 1 9n + 1 5(11n + 1) 3(13n + 1) 47 50 45 12 = 24 + (n + 1) - + - , 3n + 1 9n + 1 11n + 1 13n + 1 9 32 7 (n + 1)2 - + - , = 32n 2(9n + 1) 5(11n + 1) 3(13n + 1) 50 45 12 1 + - + , = (n + 1) - 2n 9n + 1 11n + 1 13n + 1 94n n+1 + . = 11 - 2n 3n + 1 A1 = A2 A3 A4 A5 A6 , (36) Разрешающее уравнение (35) с учётом (36) определяет начальную функцию T0 (z) при заданной температуре стенки Tw (z). Так как в рассматриваемом случае диффузионная составляющая теплопроводности жидкости в осевом направлении z не учитывалась, для однозначного интегрирования уравнения (35) необходимо задать соответствующие краевые условия только на входе трубы (z = z∗ = 0). В силу того что уравнение (35) содержит третью производную от T0 (z), в точке z = z∗ нужно получить значения для T0 , dT0 /dz, d2 T0 /dz 2 , которые определим из традиционного краевого условия, задающего распределение температуры T во входном сечении трубы [8, 10] (см. первое равенство (24) при условии ∂T∗ /∂t ≡ 0). Вновь используя метод взвешенных невязок для приближённого выполнения первого равенства (24), ∗ получим, что справедливыми остаются равенства (25), где θm (m = 0, 1, 2) - постоянные величины. Подставим в (25) разложение (7) при M = 5 и учтём равенства (33), (34), (29), тогда после приведения подобных слагаемых будем иметь 2 2 r4 v 2 d2 T0 n + 1 r0 v0 dTw r0 v0 dT0 + Bm3 0 20 - Bm4 + Bm5 Tw + a dz a dz 2 n a dz n+1 2 ∗ 2 + Bm6 + Bm7 r0 W = 2 θm , m = 0, 1, 2, z = 0, (37) n r0 Bm1 T0 + Bm2 где B01 = 30 , 47 B02 = B04 = B06 = - B11 = 75 , 329 B14 = B16 = B21 = B12 = 1 , 5640 115 , 987 B24 = B26 = - 588 13 , B05 = 45120 B15 = B22 = 5 , B25 = 42112 219 , 5640 17 , 47 89 , 4935 179 , 658 107 , 10528 214 , 987 B03 = B07 = B13 = B17 = B23 = B27 = 43 , 90240 1 , 235 103 , 421120 107 , 39480 61 , 421120 59 . 31584 (38) Метод дополнительных граничных условий в задаче теплопереноса . . . Если, как и предполагалось, температура стенки Tw задана, то в (37) все слагаемые, содержащие сомножители W , Tw , dTw /dz (при z = z∗ = 0), извест∗ ны и их следует перенести в правую часть, где θm , согласно правым равенствам в (25), также известные величины. Таким образом, преобразованная система (37) с учётом (38) представляет собой замкнутую СЛАУ, состоящую из трёх уравнений и содержащую в качестве неизвестных величин T0 , dT0 /dz, d2 T0 /dz 2 при z = z∗ = 0. Решив эту систему, определим все необходимые краевые условия для однозначного интегрирования уравнения (35). Если температура жидкости на входе в трубу постоянна по сечению z = = z∗ = 0, т. е. (см. первые равенства (24)) ¯ T (z∗ , r) = T∗ (r) = T∗ = const, (39) то выражения для величин в правых частях (37) согласно (25) имеют вид ∗ θm = 2 r0 ¯ T∗ = const, 2(m + 1) m = 0, 1, 2. (40) В случае, когда для аппроксимации краевого условия на входе трубы (см. (39) или первое равенство (24)) используется метод наименьших квадратов, следует рассмотреть минимизацию функции r0 ¯ ¯ ¯ J(T0 , T0 , T0 ) = 2 ¯ ¯ ¯ T (z, r; T0 , T0 , T0 ) - T∗ (r) dr, z = z∗ = 0 (41) 0 трёх переменных dT0 d2 T0 ¯ ¯ T0 ≡ , T0 ≡ , z = z∗ = 0. (42) dz dz 2 Необходимое условие минимума функции (41) приводит к трём равенствам ¯ T0 ≡ T0 (z), r0 0 0 r0 r0 0 ∂T ¯ ¯ ¯ T (z, r; T0 , T0 , T0 ) - T∗ (r) ¯ dr = 0, ∂ T0 ∂T ¯ ¯ ¯ T (z, r; T0 , T0 , T0 ) - T∗ (r) ¯ dr = 0, ∂T ∂T ¯ ¯ ¯ T (z, r; T0 , T0 , T0 ) - T∗ (r) ¯ dr = 0, ∂ T0 (43) z = z∗ = 0. Подставляя в (43) разложение (7) и учитывая зависимости функций T2m (1 m M ) от величин (42) (см. (29), (33)), после элементарных преобразований будем иметь систему M Dlm T2m (z) = Tl∗ , l = 0, 1, 2, z = z∗ = 0, M = 5, (44) m=0 где Tl∗ r0 M 1 = r0 Clm 0 M Dlm = j=0 m=0 r r0 Clj , 2(m + j) + 1 2m T∗ (r)dr, (45) l = 0, 1, 2, M = 5; 589 Я н к о в с к и й А. П. ненулевые значения коэффициентов Clj следующие: C00 C05 C15 C25 = C11 = C22 = 47, C03 = -200, C04 = 225, = -72, C13 = -174, C14 = 184, = -57, C23 = -114, C24 = 93, = -26. (46) Подставим в (44) выражения (29), (33) и соберём подобные слагаемые, после чего с учётом (42), (45), (46) окончательно получим r 2 v0 ¯ r4 v2 ¯ ¯ Fl1 T0 + Fl2 0 T0 + Fl3 0 20 T0 = a a M = Tl∗ 2 Dlm d2m r0 W + e2m - m=3 n + 1 2 v0 dTw r + W + f2m Tw , n 0 a dz l = 0, 1, 2, z = z∗ = 0, M = 5, (47) где Fl1 = Fl2 = Fl3 = 1 47Dl0 + 47 5 C0m Dlm , m=3 5 1 47Dl1 + 188 (48) C1m Dlm , m=3 5 1 47Dl2 + 3008 C2m Dlm , l = 0, 1, 2; m=3 коэффициенты d2m , e2m , f2m определены в (34). В правые части равенств (47) перенесены все известные в сечении z = z∗ = 0 величины. Система трёх уравнений (47) с учётом (34), (45), (46), (48) образует замк¯ ¯ ¯ нутую СЛАУ относительно неизвестных T0 , T0 , T0 . После решения этой системы с учётом обозначений (42) получим в рамках метода наименьших квадратов (41) все необходимые краевые условия для однозначного интегрирования разрешающего уравнения (35). 3. Модельная задача. В качестве модельной задачи рассмотрим случай стационарного теплопереноса при постоянной температуре стенки Tw (z) = Tw = const, z 0. (49) Для удобства дальнейшего анализа результатов расчетов введём в рассмотрение безразмерные переменные и величины r2 W 3n + 1 β T0 (z) - Tw ¯ , W = ¯ 0 = , ¯∗ - Tw 2n 0.097 T T∗ - Tw a z v r0 ¯ n+1 Z = 2z = , Pe = , v= ¯ v0 , r0 Pe a 3n + 1 v r0 ¯ θ0 (Z) = (50) ¯ где v - средняя скорость жидкости по сечению [10]; Pe - число Пекле; T∗ - ¯ характерное значение температуры жидкости на входе в трубу; β - безразмерный параметр, характеризующий влияние диссипации механической энергии на теплообмен и введённый в [10]; числовой множитель 1/0.097 введён 590 Метод дополнительных граничных условий в задаче теплопереноса . . . для того, чтобы значения параметра β в настоящем исследовании и в [10] полностью совпадали. Разрешающее уравнение (35) с учётом (49), (50) примет вид E3 d2 θ0 dθ0 d3 θ0 ¯ + E2 2 + E1 + 240θ0 = -A6 W = const, 3 dZ dZ dZ Z 0, (51) где E1 = N A3 , E2 = N 2 A2 , E3 = N 3 A1 , N= 3n + 1 ; n+1 (52) постоянные Ai определены в (36). Предположим, что на входе в трубу (z = z∗ = 0) температура жидкости постоянна по сечению (см. (39)), т. е. рассматривается аналог задачи Гретца- Нуссельта. Тогда на основании метода взвешенных невязок для однозначного интегрирования уравнения (51) из системы (37) с учётом (40), (49), (50), (52) получим следующие краевые условия при Z = 0: Bm3 N 2 dθ0 d2 θ0 + Bm2 N + Bm1 θ0 = 2 dZ dZ 1 n+1 ¯ = - Bm6 + Bm7 W , m+1 n m = 0, 1, 2, (53) где значения коэффициентов Bmi определены в (38). Если же необходимые краевые условия для уравнения (51) определяются методом наименьших квадратов, то из системы (47) с учётом (39), (45), (49), (50), (52) получим Fl3 N 2 d2 θ0 dθ0 + Fl3 N + Fl1 θ0 = dZ 2 dZ n+1 n+1 Dl3 - 164 - 16 Dl4 + = Dl0 - 104 - 9 n n ¯ n+1 W + 60 - 7 Dl5 , l = 0, 1, 2, Z = 0, (54) n 752 где значения коэффициентов Fli , Dlj определены в (45), (48). На основании (7), (29), (33), (49), (50), (52) безразмерная температура жидкости в трубе имеет выражение M θ2m (Z)R2m , θ(Z, R) = θ0 (Z) + 0 R 1, Z 0, M = 5, (55) r N dθ0 N 2 d2 θ0 , θ2 (Z) = , θ4 (Z) = , r0 4 dZ 64 dZ 2 dθ0 d2 θ0 n+1 ¯ θ2m (Z) = a2m θ0 + b2m N + c2m N 2 2 + d2m + e2m W , dZ dZ n m = 3, 4, 5; (56) m=1 где T - Tw θ= ¯ , T∗ - Tw R= 591 Я н к о в с к и й А. П. коэффициенты a2m , b2m , c2m , d2m , e2m определены в (34). Согласно (50), (52) уравнение (51) является обыкновенным неоднородным дифференциальным уравнением третьего порядка с постоянными коэффициентами и постоянной правой частью, поэтому его решение имеет вид ¯ A6 W θ0 (Z) = - + 240 3 Gm exp(αm Z), Z 0, (57) m=1 где αm - характеристические числа уравнения (51); Gm - постоянные интегрирования, определяемые из краевых условий, например, (53) или (54). В качестве конкретного примера рассмотрим течение расплавленного полиэтилена высокого давления П153-01К (n = 0.463 [10]). При этом значении n корни характеристического уравнения, соответствующего (51), с учётом (52) имеют следующие значения (см. (57)): α1 = -3.9547, α2 = -23.2203, α3 = -42.7902. (58) На рис. 1 пунктирные линии характеризуют заданное распределение температуры жидкости на входе в трубу (Z = 0), а сплошные кривые - аппроксимации этих температур, полученные разными методами. Горизонтальная пунктирная прямая 1 на рис. 1, а соответствует безразмерной температуре жидкости, заданной постоянной по входному сечению трубы (см. (39) с учётом (56)). Кривая 2 получена при Z = 0, β = 0 (отсутствие диссипации) по формулам (53), (55)-(57) с учётом (58), т. е. при использовании метода взвешенных невязок для аппроксимации краевого условия (39). Кривая 3 определена при Z = 0, β = 0 на основании соотношений (54)-(58), т. е. при использовании метода наименьших квадратов для аппроксимации краевого условия для температуры жидкости на входе в трубу. Сопоставление кривых 2 и 3 с линией 1 показывает, что применение метода взвешенных невязок (в его варианте (53) или (37), (38)) оказалось в данной задаче неудачным (кривая 2). Метод же наименьших квадратов (кривая 3) позволяет удовлетворительно аппроксимировать температуру жидкости, заданную на входе в трубу. В работе [10] для приближённого решения рассматриваемой задачи использовался другой аналитический (а именно, вариационный) метод с применением преобразования Лапласа по продольной координате z. Решение, приведённое в [10], имеет следующую структуру: 3 3 (γmi + βδmi ) exp(χi Z)(1 - R2m ), θ(Z, R) = m=1 i=1 0 R 1, Z (59) 0, где [10] χ1 = -4.0008, χ2 = -14.8547, χ3 = -145.3775; (60) числовые значения коэффициентов γmi , δmi приведены в формуле (2.34) в [10]. Сравнение значений (58), (60) показывает, что первые характеристические числа, полученные разными методами, в этой задаче практически совпадают, а вторые и третьи характеристические числа (отвечающие за распределение температурного поля на входном термическом участке [8]) различаются 592 Метод дополнительных граничных условий в задаче теплопереноса . . . a b Рис. 1. Аппроксимация температуры жидкости на входе в трубу: а) случай без учета температуры стенки; b) случай с учетом температуры стенки [Figure 1. An approximation of the fluid temperature at the inlet of the pipe in case without the wall temperature (a), and in case with the wall temperature (b); dashed lines describe the predetermined fluid temperature at the inlet of the pipe; solid lines describe the fluid temperature approximations (they were obtained by various methods); line 1 in Fig. 1, a describe the dimensionless fluid temperature at the inlet of the pipe (look at Eq. (39) with considering (56)); line 2 in Fig. 1, a is obtained by Eqs. (53), (55)-(57) with considering Eqs. (58) when Z = 0, β = 0 (the parameter β characterizes the mechanical energy dissipation, look at Eq. (50)); line 3 in Fig. 1, a is obtained by Eqs. (54)-(58) when Z = 0, β = 0; line 4 in Fig. 1, a is obtained by Eqs. (59), (60) when Z = 0, β = 0; line 1 in Fig. 1, b corresponds to the boundary condition (61) when η = 10; line 2 in Fig. 1, b approximates the line 1 (Fig. b) when β = 0 (look at Eq. (50)); line 3 in Fig. 1, b approximates the line 1 (Fig. b) when β = 1; dashed line 4 in Fig. 1, b describes the relationship (59) obtained when β = 1 with considering (65); solid line 4 in Fig. 1, b approximate the line 4 (Fig. b) by the least squares method (look at Eq. (54) with considering Eqs. (62), (63), (65))] 593 Я н к о в с к и й А. П. существенно. Кроме того, сопоставление приближённых решений (55), (59) с учётом (56), (57) указывает на то, что оба решения представляют собой линейную комбинацию трёх экспоненциальных функций, зависящих от продольной координаты Z, но решение (55) по сравнению с (59) имеет на два порядка более высокую аппроксимацию по переменной R, т. е. более точно описывает распределение температуры в поперечных сечениях трубы. Кривая 4 на рис. 1, а определена по формулам (59), (60) при Z = 0, β = 0. Сравнение кривых 3 и 4 с линией 1 показывает, что кривая 3 лучше аппроксимирует пунктирную прямую 1, чем линия 4. Дополнительные расчёты показали, что при наличии диссипации (например, при β = 1) расчётные кривые визуально практически не отличаются от кривых 3 и 4. Строго говоря, краевое условие (39) не согласовано с граничным усло¯ вием на стенке (T∗ (r0 ) = T∗ = Tw ), это приводит к необходимости построения обобщённого (разрывного) решения задачи теплопереноса в жидкости, протекающей по трубе [15]. Приближённые же методы, разработанные в настоящем исследовании или в [8, 10], позволяют построить решения, близкие лишь к классическим решениям соответствующих краевых задач, в которых все краевые условия согласованы. Исходя из этого зададим на входе в трубу краевое условие для температуры жидкости, согласованное с граничным условием на стенке трубы. С этой целью вместо (39) используем, например, следующее краевое условие: ¯ ¯ T (z∗ , r) = T∗ (r) = T∗ - (T∗ - Tw ) exp(-η(r0 - r)), (61) 0 r r0 , η = const > 0. Из равенства (61) следует, что T (z∗ , r0 ) = T∗ (r0 ) = Tw , т. е. условие согласования выполняется. Однако при подстановке (61) в (6) равенства в последнем строго выполняться не будут, но при достаточно больших значениях параметра η можно считать, что условия (6) выполняются приближённо с приемлемой точностью, по крайней мере при m = 1. Так, на рис. 1, b пунктирная кривая 1 соответствует краевому условию (61) при η = 10, а кривые 2 и 3 аппроксимируют линию 1 при β = 0 и β = 1 соответственно. Кривые 2 и 3 определены методом наименьших квадратов, т. е. на основе системы (54), из правых частей которой согласно (61) нужно вычесть величины Pl (l = 0, 1, 2), определяемые так: P0 = 47R0 - 200R6 + 225R8 - 72R10 , P1 = 47R2 - 174R6 + 184R8 - 57R10 , P2 = 47R4 - 114R6 + 93R8 - 26R10 , (62) где [16] Ri = 1 eη 1 Ri eηR dR = 0 = 1 Ri eηR + eη η i (-1)k k=1 i(i - 1) . . . (i - k + 1) i-k R η k+1 1 . (63) 0 Из сравнения кривых 2 и 3 с линией 1 на рис. 1, b видно, что в случае задания согласованного краевого условия на входе в трубу метод дополнительных граничных условий позволяет хорошо аппроксимировать входную 594 Метод дополнительных граничных условий в задаче теплопереноса . . . температуру жидкости, причём при наличии диссипации механической энергии (см. кривую 3) аппроксимация улучшается. Если потребовать, чтобы краевое условие на входе в трубу было согласовано не только с основным граничным условием на стенке (T∗ (r0 ) = Tw ), но и с дополнительным граничным условием (30), то после подстановки (61) в (30) с учётом (50) получим ¯ (r0 η)2 + r0 η = W . (64) Так как согласно (61) интерес представляет только значение η > 0, из (64) имеем √ -1 + D ¯ , D = 1 + 4W , (65) η= 2r0 ¯ откуда следует, что при W > 0 (β > 0) выполняется строгое неравенство η > 0. На рис. 1, b пунктирная кривая 4 характеризует зависимость (61), полученную при β = 1 с учётом (65), а сплошная линия 5 соответствует аппроксимации кривой 4 по методу наименьших квадратов (см. (54) с учётом (62), (63), (65)). Как видно из рис. 1, b, кривые 4 и 5 визуально практически не различаются, что свидетельствует об очень хорошей степени аппроксимации краевого условия для температуры на входе в трубу. На рис. 2, 3 изображены профили температур в разных сечениях трубы, определённые без учёта (β = 0, рис. 2) и с учётом (β = 1, рис. 3) механической диссипации энергии. На рис. 2, а и 3, а приведены результаты расчётов, выполненных на основе метода дополнительных граничных условий, а на рис. 2, б и 3, б для сравнения приведены результаты расчёта по формулам из [10] (см. (59), (60)). Пунктирная горизонтальная прямая на рис. 2, 3 характеризует заданную постоянную безразмерную температуру на входе в трубу (см. (39), (56)), остальные же линии определены на основе указанных методов. Кривые 1 определены при Z = 0, линии 2 - при Z = 0.01, кривые 3 - при Z = 0.05, линии 4 - при Z = 0.1, кривые 5 - при Z = 0.5, линии 6 - при Z = 1, т. е. кривые 1 аппроксимируют горизонтальные пунктирные прямые. Согласно этому, линия 1 на рис. 2, а совпадает с кривой 3 на рис. 1, а, а линия 1 на рис. 2, b совпадает с кривой 4 на рис. 1, a. Из сравнения кривых с одинаковыми номерами на рис. 2, а и 2, b следует, что при отсутствии диссипации (β = 0) расчёт по формулам из [10] (рис. 2, b) приводит к более быстрому выравниванию температуры жидкости и стенки вдоль трубы, чем в расчётах по методу дополнительных граничных условий (рис. 2, а). Ещё большее различие в расчётах наблюдается при учёте диссипации (β = 1). Действительно, в этом случае расчёт по методу дополнительных граничных условий приводит к тому, что температурное поле вдоль трубы достаточно быстро стабилизируется и уже при Z 0.5 почти не изменяется (см. кривые 5, 6 на рис. 3, a, которые визуально неразличимы), причём на оси трубы (R = 0) температура практически не отличается от температуры на входе (θ(Z, 0) ≈ 1). Согласно же расчёту по формулам, приведённым в [10] 595 Я н к о в с к и й А. П. a b Рис. 2. Профили температуры жидкости в разных сечениях трубы, рассчитанные без учёта диссипации механической энергии (β = 0) по методу дополнительных граничных условий (a); по методу, предложенному в [10] (b) [Figure 2. Temperature profiles for fluid at the different sections of the pipe calculated without taking into account mechanical energy dissipation (when β = 0) by the additional boundary conditions method (look at Eqs. (54)-(58)) (a), and by the method proposed in [10] (look at Eqs. (59), (60)) (b); lines with the label 1 determined at Z = 0; lines with the label 2 determined at Z = 0.01; lines with the label 3 determined at Z = 0.05; lines with the label 4 determined at Z = 0.1; lines with the label 5 determined at Z = 0.5; lines with the label 5 determined at Z = 1; dashed lines characterize the given dimensionless constant temperature at the inlet of the pipe (look at Eqs. (39), (56)); line 1 in Fig. 2, a coincides with the line 3 in Fig. 1, a, and line 1 in Fig. 2, b coincides with the line 4 in Fig. 1, a] 596 Метод дополнительных граничных условий в задаче теплопереноса . . . a b Рис. 3. Профили температуры жидкости в разных сечениях трубы, рассчитанные c учётом диссипации механической энергии (β = 1) по методу дополнительных граничных условий (a); по методу, предложенному в [10] (b) [Figure 3. Temperature profiles for fluid at the different sections of the pipe calculated with taking into account mechanical energy dissipation (when β = 1) by the additional boundary conditions method (look at Eqs. (54)-(58)) (a), and by the method proposed in [10] (look at Eqs. (59), (60)) (b); lines with the label 1 determined at Z = 0; lines with the label 2 determined at Z = 0.01; lines with the label 3 determined at Z = 0.05; lines with the label 4 determined at Z = 0.1; lines with the label 5 determined at Z = 0.5; lines with the label 5 determined at Z = 1] (рис. 3, b), даже при наличии существенной диссипации механической энергии температура жидкости вдоль трубы существенно изменяется, постепенно выравниваясь с температурой стенки. Это существенное различие в профилях температуры, приведённых на рис. 3, a и 3, b, объясняется, очевидно, тем, что кривые, изображённые на 597 Я н к о в с к и й А. П. рис. 3, a, получены по формулам, которые имеют на два порядка более высокую аппроксимацию температуры по R, чем формулы, по которым рассчитаны кривые на рис. 3, b (см. соответственно (55)-(58) и (59), (60)).

About the authors

Andrey P Yankovskii

Khristianovich Institute of Theoretical and Applied Mechanics, Siberian Branch of the Russian Academy of Sciences

Email: lab4nemir@rambler.ru
(Dr. Phys. & Math. Sci.; lab4nemir@rambler.ru), Leading Research Scientist, Lab. of Fast Processes Physic 4/1, Institutskaya st., Novosibirsk, 630090, Russian Federation

References

Supplementary files