Диссертация (1145283), страница 23
Текст из файла (страница 23)
Дляэтого выполняются пункты 9) – 11) алгоритма п. 2.4.1 параграфа 2.4. В принятых здесь обозначениях они записываются следующим образом.9) По найденному массиву температуры T4s [j], j = 0, . . . , N4 рассчитываются в s-й итерации потоки тепла от фронта оледенения в лед (qI ) и отольда во второй слой (qII ):qI = λ4 (T4s [N4 ] − T4s [N4 − 1])/h,qII = λ4 (T4s [1] − T4s [0])/h.10) Рассчитывается толщина прироста слоя льда hs за время шага τ s :hs = τ s (qI − q3 )/(γρ4 + α).11) Определяется величина нового шага τ s+1 по величине шага τ s и повеличинам потоков qI , qII :τ s+1 = τ s + (γρ4 + α)(h − hs )/(qII − q3 ).Поток q3 , входящий в условие Стефана (2.63), считается, как отмечалось взамечании 1 параграфа 2.2 этой главы, заданным и неизменным.12) Проверяется выполнение условия окончания итерационного процессарасчета шага по времени τ [n + 1] на (n + 1)-м временном слое.
Итерационныйпроцесс считается завершенным при выполнени неравенства: s+1τ− τ s ≤ ε,(2.83)ε — заданная величина.Если неравенство τ s+1 − τ s ≤ ε не выполняется, то цикл по итераци-ям не завершен, т.е. повторяется расчет массивов температур во всех слоях,включая слой льда, с новым шагом по времени τ s+1 и далее расчет по алгоритму пунктов 9) – 12). В противном случае шагу по времени τ [n+1] на новомвременном слое присваивается значение τ s+1 в последней итерации, а массивам температур Tin+1 [j], j = 0, 1, . .
. , Ni ,i = 1, 2, 4 на новом временномслое присваиваются массивы температур, найденные в последней итерации.139По представленному алгоритму была отлажена программа «ЛЛ» расчетанестационарной задачи оледенения многослойной боковой стенки цилиндра вморской воде. Программа «ЛЛ» позволяет рассчитать нестационарные распределения температуры во всех слоях, включая слой льда, толщину слояльда, время оледенения, величину шага по времени, скорость оледенения ивеличину суммарного потока тепла от воды к фронту замерзания. О выбореэффективныех параметров модели ЛЛ говорилось выше. В тестовых расчетах использовался набор параметров 5. Как следует из решения приведенногодалее уравнения (5.2.∗ ), максимальная (равновесная) толщина льда для на-бора параметров 5 равна: y∗ = 18.2378 см. В таблице 2.10 приведена частьрезультатов расчета динамики оледенения многослойного цилиндра по вычислительному алгоритму этого параграфа.Таблица 2.10y,см0.30.5135t,час13.06021.12142.144140.317265.898τ ,часммW, час3.9354.0334.2875.4717.0180.2540.2480.2330.1830.142Продолжение таблицы 2.10y,смt,часτ ,часммW, час791112427.989 641.103 930.573 1117.6059.12012.13216.81122.4720.1100.0820.0600.040По величинам скорости оледенения W (t) из таблицы 2.10 рассчитаны солености морского льда s4 (W ) по формуле Цурикова (2.3), результаты расчетаприведены в таблице 2.11.Таблица 2.11.y, м0.003 0.005W , мм/час0.254 0.248 0.233 0.183 0.142 0.110 0.082 0.060s4 (W ), 0/008.938.850.018.6450.037.880.057.1350.076.440.095.700.115.0140Как следует из таблицы 2.11, средняя соленость нарастающего морскогольда по Цурикову для W (t), вычисленной по набору 5, составила величинуsc = 7.30 /00 .
Таким образом, выбор средней солености нарастающего льдаs4 = 90 /00 можно признать удовлетворительным. В параграфе 2.7 приведеносравнение расчета динамики оледенения по модели ЛЛ для разных наборовпараметров и исследована чувствительность модели ЛЛ к их вариациям.2.5.2. Приближенные нестационарные и квазистационарныемодели ЛЛ.11, ЛЛ.1, алгоритмы численного решенияИз анализа расчетов для набора 5 по алгоритму решения нестационарной задачи, приведенному выше, следует, что в начале процесса оледененияраспределение температуры в слое бетона отличается от квазистационарногораспределения.
Это отличие сказывается и на динамике нарастания льда. Ниже приведены уравнения квазистационарной модели. Расчет по нестационарной модели, представленный в таблице 2.10, сравнивался с расчетом этого жеварианта по квазистационарной модели (решение по ней, как показано далее,сводится к решению обыкновенного дифференциального уравнения (2.96)),результат сравнения приведен в таблице 2.12.Таблица 2.12y, см0.30.51.03.05.0t, час13.060 21.121 42.144 140.317 265.898y1 , см0.3090.5091.0103.0005.004t1 , час12.0020.0041.00138.00263.00Здесь в первых двух строчках приведены данные из таблицы 2.10, а втретьей и четвертой строках — данные, полученные в результате расчета того же процесса оледенения по квазистационарной модели (2.96).
Нетрудновидеть, что квазистационарная модель дает большую скорость нарастанияльда, чем нестационарная модель. Как показал анализ, это связано с рядомфакторов, в частности с тем, что в квазистационарной модели температурана внешней границе второго слоя (состоящего из бетона) ниже, чем в нестационарной модели. Коэффициент температуропроводности бетона для набора141параметров 5 равен: a2 = 0.79997 · 10−6 ( мс ), он меньше, чем коэффициент2температуропроводности льда, равный: a4 = 0.10997 · 10−5 ( мс ). Это приводит к тому, что в начале динамики оледенения в слое бетона даже при2медленном процессе нарастания льда не успевает сформироваться квазистационарное распределение температуры.Рассмотрим, как изменится ситуация при медленном понижении температуры газа T0 (t).
С этой целью вместо условия T0 = 265.236 (K) былрассмотрен более реалистичный для морских газопроводов закон поведениятемпературы в сечениях, удаленных от входа:T0 (t) =n1+ n3 .n2 + t(2.84)В законе (2.84) время t размерное (в секундах). Рассматривались разные значения параметров n1 , n2 , n3 . Если температура газа, начиная с 273.15Кельвина, понижается по закону (2.84), опускаясь до значения 265.236 Кельвина за одни сутки, то, как показывают расчеты, процесс оледенения наступает раньше, чем успевает установиться квазистационарное распределениетемпературы в слое бетона.
Первый слой (слой стали) обладает сравнитель2но высоким коэффициентом температуропроводности: a1 = 0.5333 ·10−5 ( мс ),поэтому в нем квазистационарное распределение температуры наступаетраньше, чем в слоях бетона и льда. При n1 = 194634.723 (c · K),n2 =21834.7233 (c), n3 = 264.236 (K) температура газа понижается по закону(2.84), опускаясь до значения 265.236 Кельвина за двое суток. Как показанодалее, при таком поведении температуры газа практически во всех слоях кмоменту начала оледенения успевает установиться квазистационарное распределение температуры, и при медленном нарастании слоя льда расчет понестационарной модели совпадает с расчетом по квазистационарной моделиЛЛ.1, приведенной далее.Для ряда режимов возможны промежуточные ситуации, в которых нарастание слоя льда и распределение в нем температуры достаточно точномоделируется квазистационарной моделью, а распределения температуры вслоях обшивки далеки от квазистационарных.
В этих ситуациях справедливаследующая модель ЛЛ.11142Приближенная нестационарная модель ЛЛ.11Тепловые процессы в слоях обшивки газопровода моделируются как вмодели ЛЛ нестационарными уравнениями теплопроводности (2.53), (2.57)с начальными и граничными условиями, аналогичными соответствующимусловиям в модели ЛЛ, а модель динамики оледенения (2.60) – (2.63) заменяется следующим квазистационарным вариантом, который для t > t0 имеетвид:T4 (r, t) = A4 (t) + B4 (t) ln r,T4 (R2 , t) = T2 (R2 , t), T4 (R2 + y, t) = T∗ ,dy∂T4 ,−q=γ̃λ43∂r R2 +y(t)dtγ̃ = γρ4 + α.Расчет толщины слоя льда в этой модели при известном значении T2 (R2 , t) —температуры внешней поверхности газопровода — осуществляется по обыкновенному дифференциальному уравнению, следующему из условия Стефанадля квазистационарного распределения температуры в слое льда:dyλ4 (T∗ − T2 (R2 , t))q3=− .ydtγ̃(R2 + y) ln(1 + R2 ))γ̃(2.85)Тепловые уравнения в слоях обшивок решаются по неявным разностнымсхемам численно методом прогонки с фиксированным шагом по времени τ .Алгоритм численного решения системы уравнений модели ЛЛ.11 существенно проще, чем алгоритм численного решения системы уравнений модели ЛЛ.(n+1)Приведем уравнение для расчета температуры T2[N 2] на границевторого слоя и слоя льда (здесь используются те же обозначения, что и вчисленном алгоритме модели ЛЛ).(n+1)T2[N 2] =ξ=Q2 [N 2] + ξ T∗,(1 − P2 [N 2] + ξ)h λ4,λ2 R2 ln (1 + y (n) /R2 )(2.86)(2.87)143где h — величина шага по радиусу, P2 [N 2] и Q2 [N 2] — прогоночные коэффициенты во втором слое в граничном узле N 2.В равенствах (2.86–2.87) толщина льда y (n) берется на предыдущем временном слое.
Новая толщина льда y (n+1) находится из численного решения(n+1)уравнения (2.85), в котором T2 (R2 , t) полагается равным T2[N 2].Проведенные расчеты показали, что при достаточно малом шаге τ повремени для численного решения уравнения (2.85) можно воспользоваться,например, методом Эйлера. В этом случае y (n+1) определяется равенством:y (n+1) = y (n) + τ(n+1)λ4 (T∗ − T2[N 2])(n)γ̃(R2 + y (n) ) ln(1 + yR2 )−q3γ̃!.(2.88)Так же, как при решении уравнения (2.39), для начала счета по этому алгоритму необходимо задать начальную толщину слоя льда y0 , отличную отнуля, поскольку при y0 = 0 знаменатель дроби в (2.88) обращается в ноль.Значение толщины слоя льда y0 , образовавшегося за промежуток времениот t0 до (t0 + τ ), можно найти из приближенного равенства, следующего изусловия Стефана (2.63):γ̃λ4 (T∗ − T2 (t0 ))y0− q3 → y02 + e1 y0 − e2 = 0,=τy0q3 τλ4 (T∗ − T2 (t0 ))τ, e2 =.γ̃γ̃Физический смысл имеет положительный корень этого квадратного уравнеe1 =ния, равный:re1e21y0 = − ++ e2 .(2.89)24Алгоритм численного решения уравнений модели ЛЛ.11 реализован ввиде программы, вошедшей к программный комплекс «Лед» [70].
Проведенные расчеты показали, что с увеличением времени t0 достижения условийвозникновения оледенения уменьшается значение потока тепла от воды q3 ,соответственно уменьшается максимальная толщина слоя льда y∗ . Для набо-ра параметров 5 при законе поведения температуры газа (2.84) с параметрами n1 , n2 , n3 равными: n1 = 194634.723 (c · K), n2 = 21834.7233 (c), n3 =144264.236 (K), при условии равенства в начальный момент времени температуры во всех слоях обшивки газопровода температуре окружающей соленойводы T ∗ = 272.15 K возникновение оледенения наступают через t0 = 5.41167часа, в этот момент времени поток от воды равен: q3 = 15.87534 Дж/(c · м2 ).Максимальная толщина слоя льда y∗ при этих параметрах, рассчитаннаяпо приведенному далее уравнению (5.2.∗ ), равна: y∗ = 0.10395 м.