Диссертация (1145283), страница 22
Текст из файла (страница 22)
Как отмечалось в главе 1, характерное значение числа Рейнольдса для рассматриваемых задач о транспортировки газа составляет величину порядка 108 , поэтомузначения αo в рассматриваемых задачах велики. При больших значения αoграничное условие (2.55) переходит в граничное условие первого рода:t > t0 ,r=R:T1 = T0 (t).˜(2.55)Модель ЛЛ представляет собой обобщение модели Л1, рассмотренной впараграфе 2.4, на цилиндр, имеющий два слоя обшивок внешней поверхности. Как и в модели Л1, модифицированное условие Стефана (2.63) содержитдополнительный приток тепла к фронту замерзания, пропорциональный скорости нарастания морского льда с эффективным параметром α.
Модель ЛЛочевидным образом обобщается на большее число слоев.Величины λ4 , c4 , ρ4 , T∗ , γ, α, q3 в модели ЛЛ, как и в модели Л1, счи-таются постоянными. На сегодняшний день нет общепризнанной теорииформирования и изменения солености морского льда. Поэтому используютсякачественные оценки.Все характеристики нарастающего льда для модели ЛЛ рассчитывалисьпо методу, изложенному в п.2.4.2 параграфа 2.4.
Этот метод предполагаетпроведение численных экспериментов по модели ЛЛ, в результате которыхдля каждого набора параметров определяется зависимость от времени скорости оледенения W (t). По этой скорости с помощью формулы Цурикова(2.3) рассчитывается средняя соленость растущего морского льда sc , кото-132рая затем сравнивается с выбранной величиной солености, при которой были рассчитаны все теплофизические характеристика льда. Если отклонениевелико, процесс выбора sc продолжается.
Изменение скорости оледененияв модели ЛЛ, как и в модели Л1, может быть достигнуто за счет изменения эффективного параметра α. Расчеты по модели ЛЛ показали, что приα > 466235 (кДж/м3 ) рассчитанная средняя соленость нарастающего льдаоказывается ниже известных нам ее качественных оценок. Поэтому было выбрано меньшее значение α, равное: α = 330048.81 кДж/м3 .В результате численного эксперимента по модели ЛЛ был найден наборпараметров 5, при котором достигнуто удовлетворительное согласие с формулой Цурикова и с ожидаемым значением средней солености нарастающегольда, равной sc = 80 /00 . В этом наборе при α = 330048.81 кДж/м3 средняяпо слою соленость льда оказалась равной sc = 90 /00 . Остальные характеристики морского льда, определенные при солености sc = 90 /00 в интервалетемператур от T∗ = 271.236 K до T = 265.236 K были вычислены по форму-лам и таблицам параграфа 2.2.
В указанных ранее единицах измерения ихзначения равны: γ = 303000, ρ4 = 931, c4 = 2100, λ4 = 2.15. Следующийсогласованный набор параметров модели ЛЛ назовемНабор параметров 5q3 = 30.4447 Дж/(c · м2 ),T∗ = 271.236 K, T0 (t) = (T∗ − 6) К, R = 0.51 м,δ2 = 0.12 м, ρ2 = 2300 кг/м3 , c1 = 450 Дж/(кг · K), c2 = 924 Дж/(кг · K),λ1 = 24 Дж/(c · м · K), λ2 = 1.7 Дж/(c · м · K), αo = 500 Вт/(м2 · К),ρ1 = 104 кг/м3 , м, δ1 = 0.04 м, γ = 303000 Дж/кг, ρ4 = 931 кг/м3 ,c4 = 2100 Дж/(кг · K), λ4 = 2.15 Дж/(c · м · K), α = 330048.81 кДж/м3 .Этот набор соответствует цилиндру, имеющему два слоя обшивки: первыйслой, состоящий из стали, второй — из бетона.Набор параметров 5 представляет собой лишь один из возможных вариантов согласованного набора параметров процесса оледенения цилиндра сослоями в морской воде заданной солености.133Для адекватного моделирования реального морского газопровода выборэффективного параметра α и средней солености нарастающего льда должныосновываться на экспериментальных данных, полученных в реальных условиях функционирования газопровода.В модели ЛЛ уравнения (2.53), (2.57) и (2.60) – одномерные уравнениятеплопроводности в слоях обшивок и в толще льда соответственно; (2.55)– граничное условие третьего рода на внутренней поверхности газопровода;(2.54), (2.58), (2.61) – начальные условия; (2.56), (2.59) – условия равенстватемператур и потоков тепла на стыке слоев; (2.62) – условие неизменноститемпературы на фронте замерзания морской воды; (2.63) – модифицированное условие Стефана, введенное в модели Л1 параграфа 2.4 (2.21).Алгоритм численного решениясистемы уравнений модели ЛЛАлгоритм численного решения уравнений модели ЛЛ является обобщением алгоритма параграфа 2.4.
Расчет максимальной толщины y∗ слоя льдана внешней многослойной стенке газопровода приведен далее, в пункте 5.2,уравнение (2.5.2.∗ ). Численное решение задачи (2.53)–(2.63) с явным выде-лением фронта оледенения строится по неявной разностной схеме решенияуравнений теплопроводности (2.53), (2.57), (2.60) с использованием итерационной процедуры решения возникающей нелинейной системы уравнений. Каки в алгоритме п. 2.4.1 параграфа 2.4, переменный шаг по времени τ [n + 1]на (n + 1)-м временном слое определяется таким образом, чтобы слой льдаза время τ [n + 1] увеличился на заданную постоянную величину h, равнуюшагу сетки по радиусу.Введем равномерную сетку по радиусу в области первого и второго слоевобшивки цилиндра и в области слоя льда.
В области слоя льда r4 [j] = R2 +jh, j = 0, 1, . . . , N4 . Обозначим n — номер шага по времени; число N4 вмомент времени n равно: N4 = n. Максимальное значение n — величина N∗— определяется как целая часть от величины (y∗ /h), где y∗ — максимальнаятолщина слоя льда при выбранных параметрах модели ЛЛ, найденная далеев пункте 5.2 (уравнение (2.5.2.∗ )).Толщина слоя льда за n шагов по времени станет равной (y0 +nh); вводим134массив шагов по времени τ [n], n = 1, . . . , N∗ ; τ s — величина шага по временина (n + 1)-м временном слое в s-й итерации; после окончания итерационногопроцесса искомому шагу τ [n + 1] присваивается значение τ s+1 в последнейитерации; tn+1 = t0 + τ [1] + τ [2] + . .
. + τ [n + 1] — суммарное время нарастанияльда толщиной y = y0 + (n + 1)h.T4s [j] — температура льда в j-м узле на (n + 1)-м временном слое в s-йитерации, T4n+1 [j] — температура льда в j-м узле на (n+1)-м временном слое,T4n+1 [j] присваивается значение T4s+1 [j], вычисленное в последней итерациипри выполнении условия окончания итерационного процесса, величины T4n [j],τ [n] считаются известными, величины T4n+1 [j], τ [n + 1] рассчитываются поитерационному алгоритму (j = 0, 1, .
. . , N4 ).К алгоритму расчета температуры льда и величины шага по времени,приведенному в п. 2.4.1 параграфа 2.4, добавляется расчет по неявной схемеметодом прогонки массивов температур в первом и втором слоях обшивкигазопровода. Это не вносит принципиальных изменений в общую схему численного решения задачи с явным выделением фронта оледенения, однакоусложняет расчет температуры в слое льда на каждой итерации.Введем дополнительно следующие обозначения. ri [j] = Ri−1 + jh, j =0, 1, . . . , Ni , (i = 1, 2, 4) — текущее значение пространственной переменнойв слоях обшивки и в слое льда (здесь полагается, что R0 = R); величиныNi определяются как целая часть от δi /h, i = 1, 2; N4 для (n + 1)-го слояпо времени равно: N4 = n + 1, y∗ = N∗ h — максимальная толщина льда висследуемой задаче, ее расчет приведен в п.
2.5.2 (уравнение (5.2.∗ )).Tin+1 [j], (i = 1, 2, 4) — температура в j-м узле на (n + 1)-м временномслое; Tis [j], (i = 1, 2, 4) — температура в j-м узле на (n + 1)-м временномслое в s-й итерации. Величины Tin [j], j = 0, 1, . . . , Ni , (i = 1, 2, 4) и τn считаются известными, искомые величины Tin+1 [j], j = 0, 1, . . . , Ni , (i = 1, 2, 4)и τn+1 находятся по алгоритму, являющемуся обобщением алгоритма п. 2.4.1параграфа 2.4.В нулевом приближении шаг τ (0) на (n + 1)-м слое по времени задаетсяравным шагу τn на n-м временном слое: τ (0) = τn .В качестве нулевого приближения первого шага по времени τ10 можно135взять величину, определенную из разностного аналога модифицированногоусловия Стефана (2.63), при условии, что массив температуры льда равен:T41 [0] = T20 [N 2], T41 [1] = T∗ .
Из этого условия находим:τ10 =(γρ + α)h2.(λ4 (T∗ − T20 [N 2]) − q3 h)Затем значение τ10 уточняется по аналогу пунктов 9) – 12). Расчет по моделиЛЛ начинается в момент времени t0 , в который выполняются необходимыеи достаточные условия начала оледенения внешней многослойной стенки цилиндра в морской воде. Как отмечалось в главе 1, эти условия имеют вид:t < t0 → T2 (R2 , t) > T∗ , t > t0 → T2 (R2 , t) < T∗ ,∂T2 λ2> q3 .∂r R2 ,t02) Уравнения теплопроводности (2.53), (2.57), (2.60) аппроксимируются по неявной монотонной схеме Самарского и решаются методом прогонкианалогично схеме п.
2.4.1 параграфа 2.4.Разностные аналоги этих уравнений можно представить в видеλihTis [j] − Tin [j]=((1+)(Tis [j + 1]−s2τρi c i h2ri [j]−Tis [j]) − (1 −h)(Tis [j] − Tis [j − 1])),2ri [j]j = 0, 1, . . . , Ni ,i = 1, 2, 4.(2.64)Уравнения (2.64) приводятся к виду, используемому в методе прогонки,Ai [j] Tis [j − 1] − Bi [j] Tis [j] + Ci [j] Tis [j + 1] = Di [j], j = 0, 1, . .
. , Ni , i = 1, 2, 4.(2.65)Коэффициенты Ai [j], Bi [j], Ci [j], Di [j] зависят от величины шага по времени τ s , они рассчитываются по формулам:Ai [j] = η(1 −h),2ri [j]Bi [j] = 1 + 2η,λi τ sη=,ρi c i h 2h),2ri [j]Di [j] = −Tin [j],i = 1, 2, 4.(2.66)Ci [j] = η(1 +j = 0, 1, . . . , Ni ,136Рассчитываются первые прогоночные коэффициенты для первого слояобшивки газопровода.
Для граничного условия (2.55) они имеют следующийвид:C1 [1],P1 [2] =Z(ν/(1 + ν))A1 [1] T0 (tsn+1 ) − D1 [1]Q1 [2] =,ZZ = B1 [1] − A1 [1]/(1 + ν),ν = αo h/λ1 .˜ первые прогоночные коэффициенты дляДля граничного условия (2.55)первого слоя обшивки газопровода записываются в виде:C1 [1],P1 [2] =B1 [1]A1 [1] T1s [0] − D1 [1]Q1 [2] =,B1 [1]T1s [0] = T0 (tsn+1 ).Далее в цикле по j = 2, . . . , (N1 − 1) рассчитываются остальные прого-ночные коэффициенты для первого слоя по формулам:P1 [j + 1] =C1 [j],B1 [j] − A1 [j] P1 [j]Q1 [j + 1] =A1 [j] Q1 [j] − D1 [j].B1 [j] − A1 [j] P1 [j]Первые прогоночные коэффициенты для второго слоя определяются из разностного аналога условий (2.56):T1s [N1 ] = T2s [0],(2.67)λ2 s(T2 [1] − T2s [0]),λ1sT1 [N1 − 1] = P1 [N1 ]T1s [N1 ] + Q1 [N1 ],T1s [N1 ] − T1s [N1 − 1] =T2s [0] =λ2λ1 (1 − P1 [N1 ] +λ2λ1 )T2s [1] +(2.68)(2.69)Q1 [N1 ](1 − P1 [N1 ] +λ2λ1 ),T2s [0] = P2 [1]T2s [1] + Q2 [1].(2.70)(2.71)Из системы уравнений (2.67) – (2.71) находятся выражения для первыхпрогоночных коэффициентов в расчете температуры второго слоя:P2 [1] =λ2λ1 (1 − P1 [N1 ] +λ2λ1 ),Q2 [1] =Q1 [N1 ](1 − P1 [N1 ] +λ2λ1 ).(2.72)137Далее в цикле по j = 1, .
. . , (N2 − 1) рассчитываются остальные прого-ночные коэффициенты для второго слояP2 [j + 1] =C2 [j],B2 [j] − A2 [j] P2 [j]Q2 [j + 1] =A2 [j] Q2 [j] − D2 [j].B2 [j] − A2 [j] P2 [j](2.73)Первые прогоночные коэффициенты для слоя льда определяются из разностного аналога условий (2.59):T2s [N2 ] = T4s [0],(2.74)λ4 s(T4 [1] − T4s [0]),λ2sT2 [N2 − 1] = P2 [N2 ]T2s [N2 ] + Q2 [N2 ],T2s [N2 ] − T2s [N2 − 1] =T4s [0] =λ4λ2 (1 − P2 [N2 ] +λ4λ2 )T4s [1] +(2.75)(2.76)Q2 [N2 ](1 − P2 [N2 ] +λ4λ2 ),T4s [0] = P4 [1]T4s [1] + Q4 [1].(2.77)(2.78)Из системы уравнений (2.74) – (2.78) находятся выражения для первыхпрогоночных коэффициентов в расчете температуры слоя льда:P4 [1] =λ4λ2 (1 − P2 [N2 ] +λ4λ2 ),Q4 [1] =Q2 [N2 ](1 − P2 [N2 ] +λ4λ2 ).(2.79)Далее в цикле по j = 1, . . . , N4 − 1 рассчитываются остальные прогоноч-ные коэффициенты для слоя льда (N4 = n + 1)P4 [j + 1] =C4 [j],B4 [j] − A4 [j] P4 [j]Q4 [j + 1] =A4 [j] Q4 [j] − D4 [j].B4 [j] − A4 [j] P4 [j](2.80)Температуре льда на фронте оледенения присваивается значение температуры фазового перехода:T4s [N4 ] = T∗ .(2.81)По найденным прогоночным коэффициентам для слоя льда, второго ипервого слоев обшивок и из условий равенства температур на стыках лед–второй слой, второй слой–первый слой, в обратных циклах рассчитываютсямассивы температур во всех слоях, включая слой льда, в s-й итерацииTis [j − 1] = Pi [j] Tis [j] + Qi [j],i = 4, 2, 1.(2.82)138После того, как массив температуры в слое льда в s-й итерации рассчитан,осуществляется расчет величины шага по времени в (s + 1)-й итерации.