Моделирование динамики температурного поля грунтов основания здания в криолитозоне Текст научной статьи по специальности «Физика»

электронное научно-техническое издание

НАУКА и ОБРАЗОВАНИЕ

Эл № ФС 77 - 30569. Государственная регистрация №0421100025. ISSN 1994-040S

Моделирование динамики температурного поля грунтов основания здания в криолитозоне

77-30569/274059

# 12, декабрь 2011

Гласко А.В., Федотов А.А., Сидняев Н.И., Храпов П.В., Мельникова Ю.С.

УДК 551.340: 51-74

МГТУ им. Н.Э. Баумана [email protected] [email protected] si ёпуаеу@уапёех. ги [email protected] [email protected]

1. Введение.

Одна из основных проблем строительства в криолитозоне (на Дальнем Севере) состоит в том, что таяние вечномерзлого грунта, в результате теплового воздействия со стороны установленного на нем отапливаемого здания, ведет к проявлению опасных криогенных процессов и, в результате, к аварийной ситуации. К опасным криогенным процессам относятся, прежде всего, следующие процессы [1].

Пучение грунта - увеличение объёма грунта, в результате его промерзания. Неравномерное поднятие поверхности грунта при пучении ведет к повреждениям построенного на поверхности грунта здания (сооружения).

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

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

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

При резком замерзании грунта, в нем образуются морозобойные трещины - трещины шириной порядка 10-15 см, рассекающие поверхность грунта и уходящие на несколько метров в глубину. Со временем в такие трещины попадает вода, а ее замерзание приводит к их росту.

Просадка - сокращение объема вечномерзлого грунта (т.е. глубинного слоя грунта, находящегося под деятельным слоем и не подверженного процессам сезонного промерзания и оттаивания), потеря им прочностных свойств и рост сжимаемости в результате его оттаивания. Неравномерная деформация многолетнемерзлого грунта при оттаивании ведет к повреждениям установленного на его поверхности здания.

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

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

2. Постановка задачи

Рассмотрим здание, установленное непосредственно на поверхность грунта (без фундамента) и поставим задачу моделирования динамики температурного поля в прямоугольной области грунта под зданием, которую назовем расчетной областью. Будем

считать, что здание отапливается и внутри него поддерживается постоянная температура (например, и = 200 С). На рис.1 представлена часть расчетной области О = {0 < х < хЬ, 0 < у < уЬ, 0 < г < гЬ} (прямоугольный параллелепипед), соответствующая первой координатной четверти ( х > 0, у > 0 ). Плоскость Оху отвечает поверхности земли. Ось Ог направлена вглубь грунта. Закрашенный прямоугольник на верхней границе обасти соответствует зданию, расположенному над этой областью. Предполагается симметрия здания, состава грунта и пр. относительно координатных плоскостей х02 и у02, поэтому моделирование динамики температурного поля достаточно провести в рассматриваемой четверти расчетной области. Область О неоднородна и состоит из нескольких литологических слоев, т.е. областей грунта с различным составом и, как следствие, с различными физическими характеристиками. К физическим характеристикам литологических слоев относятся следующие величины: рл - плотность сухого грунта,

Сй - удельная теплоемкость сухого грунта,

- суммарная влажность грунта в долях к массе сухого грунта, Лм - коэффициент теплопроводности мерзлого грунта, - коэффициент теплопроводности талого грунта,

= (и) - доля незамерзшей воды по отношению к массе сухого грунта (см.

ниже).

Рис. 1. Расчетная область

Для примера, на рис. 1 выделено 3 литологических слоя: В1, В2 и В3 (соответственно, области 1, 2 и 3 на рис. 1), Принадлежность точки М(х, у, г), соответствующей области В{ определяется следующими формулами:

В1 ={УМ(х, у, г) е В: г < 2}, В2 = {УМ(х,у,г) е В :6 < х < 21,2 < г < 6}, В3 = {УМ(х, у, г) е В: г > 6 и 0 < х < 6,2 < г < 6}.

(1)

Верхняя граница области В (2=0) тоже, вообще говоря, не является однородной. Так, в качестве примера, на рис. 2 она разбита на три зоны: 1, 2 и 3. В зоне 2 (заполненный текстурой «кирпичная стена» прямоугольник) располагается здание. Зона 1 - зона прямого контакта грунта с атмосферой (отсутствие снежного покрова или тонкий снежный покров). Зона 3 - зона снежных надувов (сугробы).

Рис. 2. Верхняя граница расчетной области

Динамика температурного поля и( X, у, 2, t) в области Д с учетом фазовых переходов - твердое тело-жидкость, описывается уравнением теплопроводности в виде

* ди д (ср + - и*))— = — дt дх

Я

ди дх

+

д_

ду

Я

ди ду

+

д_ д2

Я

ди д2

+ Б(X, у, 2,t), (2)

где с - удельная теплоемкость; р - плотность;

Я - коэффициент теплопроводности;

и( х, у, 2, t) - температура среды;

*

и - температура фазового перехода (отметим, что температура начала замерзания

грунта, вообще говоря, чуть ниже нуля и зависит от разновидности грунта и концентрации

*

порового раствора в нем [2], но здесь приближенно принимается и = 0 °С); Q - теплота фазового перехода;

S(u — u ) - дельта-функция.

s( x, y, z, t) - мощность внутренних источников тепла.

Внутренние источники тепла в задаче отсутствуют, поэтому s( x, y, z, t) = 0 Начальное условие зададим, фиксировав распределение температуры в грунте в момент времени t=0 (на практике, такое распределение может быть найдено по результатам геотехнического мониторинга):

u(x,y,z,0) = p(x,y,z), M(x,y,z) e D. (3)

Граничные условия зададим следующим образом. На границе z = 0 с температурой u = иъ происходит конвективный теплообмен со средой, имеющей

температуру 6(t). Плотность теплового потока JЪ на этой границе задается через коэффициент теплоотдачи h в виде

Jb = h ■ (0(t) — ub). (4)

На границе z = zL поддерживается постоянная температура, равная температуре окружающей среды, имеющей температуру

00 = const. (5)

Боковые границы области D теплоизолированы

du du

—(0, y, z, t) = 0, —(xL, y, z, t) = 0, dx dx

du . n . n du . . Л

—(x,0, z, t) = 0, — (x, yL, z, t) = 0. (6)

dy dy

Слагаемое cp в левой части уравнения (2) можно написать в следующем виде:

cp =

Pd'

cd + Cc (Wtot — Ww) + cwWw + к

dW

л

w

, u < u ; du ) / (7)

Pd ■ (cd + cwWtot), u > u ,

где Сссе =0.49 ккал/(кг-оС) и =1.006 ккал/(кг-оС) - удельные теплоемкости льда и воды,

соответственно, к = 79.4 ккал/кг - удельная теплота фазового перехода льда. Формула для объемной теплоемкости (7) учитывает фазовые переходы в области отрицательных

температур. Функцией Ww = Ww (u) описывается массовая доля (по отношению к массе сухого грунта) незамершей воды при температуре u .

Вследствие вышесказанного плотность влажного грунта р представляется выражением

Pd +Pd * Wot, U < U;

pd +Pd ■ (Wot - Ww) + Pd ■ Ww, u < u < u*; (8)

P =

Pd + Pd ■ Wt

tot

u > u

где и — температура замерзания грунтовой воды; т.е. при и < и вся влага находится в

*

твердом состоянии (лед), при и > и вся влага находится в жидком состоянии (вода); при и < и < и часть влаги ) находится в жидком состоянии, а часть влаги

Ра — ) — в твердом состоянии.

Зависимость = (и) описывается интерполяционной формулой [3]

Ww (u) =

а

■+rw,

(9)

в — и

где постоянные коэффициенты ОСМ!, и ум! свои для каждого литологического слоя грунта и определяются по трем значениям функции при трех различных значениях температуры и . Рис.3 иллюстрирует зависимости от температуры. Кривые 1 и 2

соответствуют двум различным литологическим слоям.

Рис. 3. Кривые незамершей воды

Коэффициент теплопроводности Л в уравнении (2) определяется так

Л =

Лм, и < и* Л

V' и > и .

Коэффициент теплоотдачи И в формуле (4) записывается в виде

и = 1

1

(11)

+ Я

а

где а - коэффициент конвективного теплообмена, Я- термическое сопротивление [3].

У-Ч *

Теплота фазового перехода О при и = и вычисляется по формуле [2]

О = Р(Ща - (и*)).

3. Решение задачи методом контрольного объема

Краевая задача (2) - (6) решалась численно. Прежде всего, коэффициент в левой части уравнения (2) сглаживается и совершается переход к обычной задаче

теплопроводности [4]. Дельта-функция 3(и — и ) заменяется некоторой функцией 3((и — и ), А), которая отлична от нуля лишь внутри некоторого промежутка сглаживания [—А, А]. В результате вместо решения (2) ищется решение уравнения со сглаженным коэффициентом

* ди д (ср + О 3((и — и*), А)) — = —

дГ дх

КУ1 1/Л

г

Л-и

V дх.

+ ■

_д_ ду

Л

ди ду

+ ■

д_ дх

Л

ди дх

+ 5( х, у, 2,г) (12)

Для аппроксимации дельта-функции используются различные формулы, которые строятся с учетом выполнения условия сохранения теплового баланса на отрезке [—А, А]. В настоящей работе используются следующие формулы [4]: ступенчатая аппроксимация

8(и — и*, А) =

Г 1 и — и * <А,

2А,

0, и — и * > А

(13)

и параболическая аппроксимация

8(и - и*, А) =

3

4А 0,

(и - и*)

и -и

и -и

< А,

> А.

(14)

Можно убедиться в том, что для обеих формул (13) и (14) выполняется условие

/Л

| 8((и - и*), А)Ди = 1.

-А

Параметр А влияет на результаты расчетов, он зависит от используемой расчетной сетки и, как правило, определяется опытным путем в результате методических расчетов.

Задача теплопроводности ((12), (3) - (6)) решалась методом контрольного объема. Суть метода изложена ниже. В начале, для простоты изложения, рассматривается одномерная задача. Потом проводится обобщение на трехмерный случай. В качестве примера рассмотрим следующую задачу: найти численное решение и = и (х, t) уравнения теплопроводности

ди д ( п ди Л Л— дх

ср— = — дt дх

V ил у

+ s(u, х, t) (15)

в ограниченной области О = {0 < х < хЬ,0 < t < 4}, удовлетворяющее начальному усло-

вию

и( х,0) = (( х)

(16)

и граничным условиям

и (0, t) = / ),

и( хЬ, t) = ), (17)

где и(х, t) - температура среды, s(u, х, t) - мощность внутренних источников тепла, ((х), /(1,) и v(t) - заданные функции.

Разобьем отрезок [0, хЬ] на т участков, которые назовем контрольными объемами (КО) (рис.4а). Точки сетки

= {xwl,I = 1,...,т +1}={х^1 = 0 < х^2 < ... < < xwm+1 = хЬ} (18)

определяют границы КО; индекс левой границы КО соответствует его номеру: 1-ый, 2-ой, ..., ш-ый контрольные объемы (КО); Иxi = xwi+1 — xwi, I = 1,..., Ш - I — ый шаг сетки (0М1

или ширина / — ого КО (в общем случае сетка предполагается неравномерной). Точки области В, в которых будем искать решение задачи, назовем расчетными. В качестве внутренних расчетных точек выберем центры КО: х1 =(xwi + xwi+12 , 1=1,.,ш+1;

добавим к этим точкам две граничные точки x0 = 0 и xш+1 = xL, и назовем сеткой ( объединенное множество внутренних и граничных расчетных точек -( = {xl, / = 0,1,..., ш +1} (рис. 4 б).

Рис. 4. Построение расчетной сетки

Фрагмент расчетной сетки С, показан на рис. 5. Точки сетки обозначены заглавными буквами W, Р и Е: Р — рассматриваемая точка (Point), а Е и W -соответственно «восточная» (East) и «западная» (West) соседние с ней точки [5,6]. В

скобках указаны индексы расчетных точек - (i — 1), i и (i + 1), соответственно.

Штриховыми линиями показаны грани КО. Для обозначения граней КО, содержащего точку Р, используются буквы w и e . Расстояние между точками W и Р обозначим как (Sx)w, а между точками Р и Е - как (Sx) . Ширину контрольного объема обозначим

через Ax.

Рис. 5. Фрагмент расчетной сетки

Запишем уравнение (15) в виде

Здесь

ср

ди дЗ дt дх

+ £.

з =

дх

(19)

(20)

- плотность теплового потока.

Проинтегрируем обе части равенства (19) по КО, содержащему точку Р(1), т.е. от w до е:

е ди е дЗ е

Г ср—dx = - Г—dx + Г sdx. J дt дх J

дх

В результате получим

аг

(и Р-и°Р) = З - Зе + Г sdx,

(21)

где иР - известное значение температуры и в точке Р в момент времени t, ир -соответствует неизвестной температуре в точке Р в момент времени t + Аt,

а°р = (ср)р Ах/А? . При интегрировании температуру внутри КО считаем постоянной и

равной температуре ир. Первые два слагаемых в правой части уравнения (21) отвечают

плотностям входящего и выходящего потоков тепла, интеграл выражает суммарную мощность генерации тепла в КО.

Для выражения 3 ^ и через температуры в расчетных точках используем кусочно-линейный профиль, изображенный штрихпунктирной линией на рис.5б. В результате получим

Л Л

=ТТ\~ ■(и* - иР ) ' =тт\~ ■( иР - иЕ ) . (22)

(Зх)„ (Зх)е

Обозначим среднюю мощность генерации тепла в КО через £ . Линеаризуем £ по температуре: £ = £с + £рир. Тогда

е

| £ёх = £Ах = (£с + £рир)Ах . (23)

В случае, когда £ не зависит от температуры, имеем £р = 0, £ = £с . Подставив (22) и (23) в выражение (21), получим дискретный аналог уравнения (19)

в виде

арир = а*и* + аЕиЕ + Ь, (24)

где

а а0 =(Ср^рАХ (ЗхЕ (Зх) ' р А?

ар = а* + аЕ + а°р - £рАх, Ь = £сАх + а(ри°р.

При аппроксимации плотности теплового потока по формулам (22) используются значения коэффициента теплопроводности Л на гранях КО. В задаче же все необходимые величины задаются в расчетных точках - центрах КО. Если теплопроводность задана только в расчетной точке, то разумно предположить, что ее значение остается постоянным по всему КО, окружающему эту точку. Другими словами, каждый КО заполнен материалом с постоянной теплопроводностью (соответствующей ее значению в расчетной точке).

Если зависимость температуры от x аппроксимировать кусочно-линейным профилем, изображенным на рис.5,б сплошной линией, то плотность теплового потока на грани e может быть найдена по следующей формуле

Л =

Лп

(8х) е

( иР - ие ) =

Л

(8х) е

(ие - иЕ)

(25)

Исключив ие из этого выражения, получим

а

Е

(8х ) е- + (8х ) е

Л

Л

-1

( ир - иЕ )

(26)

Сравнивая вторую формулу в (22) с (26), запишем выражение для коэффициента

аЕ =

Л

(8 )е

(8х ) е-+(8х ) е

Л

Л

-1

ЛЕ • Лр

ЛЕ(8х ) е-+ЛР (8х ) е

Аналогично записывается формула для коэффициента а.

ж

Л

а,

Лр • Лж

ж

(8х)w ЛР (8х) (8х),

Из (25) можно также вывести формулу для и :

и„

^рир + рЕиЕ

Р + РЕ

где Рр

Л

р,

Л

(8х 1-' (8х)

(27)

(28)

(29)

В вычислительной схеме мы получаем температуры в расчетных точках и обычно не определяем температуры на гранях контрольных объемов. Когда необходимо получить температуру на грани, нужно использовать выражение (29).

Формула (28) используется для вычисления аЖ , когда индекс / точки р принимает значения от 2 до т , а формула (27) - для I = 1,...,т -1. При I = 1 левая грань первого КО совпадает с точкой Ж(0), поэтому полагая (8х) к-= 0, получаем для аЖ

выражение

Яр

а =

(дх)*>+ (30)

На правой границе области при I = т правая грань т — ого КО совпадает с точкой Е(т + 1), т.е. здесь (бх) е+ = 0 . Для коэффициента аЕ получается формула

а* =7ТУ- (31)

(бх) е—

Введем новые обозначения для коэффициентов в уравнении (24):

а№(г) = Лг, ар(г) = Б1 , аЕ(г) = С1, Ь(г) = Вг.

С учетом построенной сетки (см. формулу (18) и рис.4 и рис.5) можем записать выражения для коэффициентов:

Л 1 = а№ (1) = Iя Л = а№ (г)" ЧЯ—1

кх1 г Я1кх1—1 + Я1_1кх1 г = 2,..., т.

1 ? ?

Сг = аЕ(г) = ,2Я+1Я - . , , Ст = аЕ(т) ^

т •

(32)

Я+1кх1 + Я1кх1+1 г = 1,..., т — 1. т кхт

1 ?

Б г = ар (г) = а№ (г) + аЕ (г) + а°р (г) — sр (г)(Лх)г = Лг + Сг + а°р (г) — ^ (г) • кхг,

а°р(г) = (( ер)1кх1 )/л; В г = Ь(г) = (г)кхг +аРР(г)ы°р(г), иРР(г) = и0.

Здесь и°- значение температуры в точке х = х 1 в предыдущий момент времени, $с (г) и 8р (г) - значения функций 8С (х, Т) и 5*р (х, Т) в точке х = хг в текущий момент времени.

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

Б1и1 — С1и2 = В1 + Л1и0, Л2и + Б2^2 С2^3 = ^2 ,

—Аиг—1 + Бгиг — СгЩ+1 = Вг , (33)

—Лт—1ит—2 + Бт—1ит—1 — Ст—1ит = Вт—1 , —Лтит—1 + Бтит = Вт + С т ит+1 .

Система (33) решается методом прогонки [7]. После нахождения поля температуры делается следующий шаг по времени и вычисления продолжаются до нужного момента времени.

Полученную систему линейных алгеьраических уравнений будем решать методом прогонки. Прямой ход метода прогонки состоит в вычислении прогоночных. коэффициентов ai и Д , (0 < i < m)

a = о, во = о,

«i = Cj( Bt - Д .ам), Д =( D + Д' в i-i)/(Bt - At-ам), i = 1,..., m .

Обратный ход метода прогонки состоит в нахождении значений неизвестных

um =Дт ; U = a ■ U+1 + Д , i = m - 1,m - 2,...,1 .

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

Достаточным условием применимости метода прогонки является преобладание по модулю диагональных элементов матрицы системы суммы модулей соседних элементов

[7]. Коэффициенты aW (i), aE(i), a°(i) положительны, в силу постановки задачи (см. формулы (32)). Нарушение указанных достаточных условий возможно лишь при положительных значениях коэффициента sp (i) (см. выражение для коэффициента ap (i) в формулах (32)). Поэтому функцию s , если это необходимо, нужно линеаризовать так, чтобы величина sp (i) была отрицательной, либо равной нулю.

Метод КО обобщается на случай трехмерной нестационарной краевой задачи для уравнения теплопроводности.

Уравнение теплопроводности

ди д cp— = — dt дх

v дх j

д + —

ду

с

. дУ,

д + —

дг

v д2 j

+ s (34)

представляется в виде

ди ср—

дг

д/ д/ д/

дх ду дг

т ,ди ди ди

где / х = —Я —, / = —Я—, / = —Я— - плотности теплового потока в направлении дх ду дг

осей х, у и г, соответственно. Расчетная область В разбивается на множество КО.

Центр каждого КО есть расчетная точка. Типичный контрольный объем представлен

на рис.6.

Рис. 6. Типичный контрольный объем

Буквы W, Р , Е, S, N, B, T обозначают точки сетки, а именно: Р — рассматриваемая точка (Point), а W, Е, S, N, B и T - соответственно, «западная» (West) , «восточная» (East), «южная» (South), «северная» (North), «нижняя» (Bottom) и «верхняя» (Top) соседние точки. В скобках указаны индексы расчетных точек. Грани КО изображаются штриховыми линиями. Для обозначения граней КО, содержащего точку Р, используются буквы w, e, s, n , b и t. Проинтегрировав уравнение (35) по КО, получим

( cp)p (Up - U°P ) = JwAw - JeAe + JsAs — J nAn + J bAb — JA + *AV (36)

Здесь верхний индекс «0» обозначает известное значение температуры в начале шага по времени At, J — плотность теплового потока через грань КО, на которую указывает нижний индекс; A — площадь соответствующей грани, например, Ae = hy(j) • hz(к); s — усредненный по объему источниковый член, A V = hx(i) • hy (j) • hz (к) — объем КО.

Потоки тепла через грани КО e и w могут быть найдены по формулам

JeAe = De • (UP — UE ) , JwAw = Dw '(UW — UP ) . (37)

Здесь De — проводимость между точками Р и Е, которая находится по значениям Л в этих точках. Проводимость D можно вычислить по формуле (см. формулы (26) и (27))

D = A

e e

2ЛеЛр ЛЕМ{ i) + ЛрМ{ i +1)

(38)

Равенства, подобные (37) и (38) , используются также для расчета соответствующих величин на других гранях. Источниковый член £ линеаризуется по температуре: £ = £с + £рир. Если £ не зависит от температуры, то полагаем £р — 0, £ = £с. Подставив выражения для J и £ в (36), получим дискретный аналог уравнения (34)

арир — аЕиЕ + ашиш + а^и^ + а^и^ + а^и^ + а^и^ + Ь, (39)

где аЕ — ве, а№ — , ам — Бп, а5 — Б,, ат — , ав — Въ, Ь — £сАУ + а°ри°р,

a0P

((cp) p AV У At, ap = aE + aW + aN + aS + aT + aB + a°p — sp A V.

Уравнение вида (39) записывается для каждого КО, содержащего внутреннюю расчетную точку. С помощью граничных условий в приграничных КО выполняются

преобразования, в результате которых, значения и на границе не будут входить явным образом в систему уравнений. В итоге получается замкнутая система линейных алгебраических уравнений, которая может быть решена любым подходящим методом. Для решения системы в настоящей работе используется, метод переменных направлений [5, 6], основанный на методе прогонки.

Суть метода переменных направлений проиллюстрируем на примере двумерной задачи. На рис. 7 изображена расчетная область в двухмерии. Через У ^ обозначен двумерный КО, другие обозначения подобны используемым на рис. 6.

Рис. 7. Расчетная область двумерной задачи

Организуется следующий итерационный процесс. Если в дискретном аналоге нестационарного двумерного уравнения теплопроводности

ари р — а Еи е + ашиш + а^и^ + а^и^ + Ь,

(40)

где аЕ — Ое, а№ — В^, ам — Бп, а5 — , Ь — ^АУ + а°ри°р, а(р —((ср)р АУ)/Аг, ар — аЕ + а№ + ам + а3 + а°р — 5р АУ

в направлении оси у предположить известными иы и и3 (из предыдущей итерации), то в нем останутся только три неизвестные: ир, иЕ и и№ . Построив такие уравнения с тремя неизвестными вдоль линии, параллельной оси х , получим

ар(г,])ир(г,]) — аЕ(1+\,])иЕ(¿+1,;) + аШ(г—!,])%(I—1,у) + Ь(г,]) , * — т , (41)

где Ь0, ]) — (г, ]+1) К (г, ]+1) + а3 (г, ]—1)и^ (г, ]—1) + Ь(г, ]), а и обозначает значение и в

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

Метод переменных направлений состоит в том, что вначале прогонка используется для всех линий, параллельных оси Ох, а потом повторяется для всех линий, параллельных оси Оу. Рассчитанные новые значения и вдоль линии используются в

качестве значений и при решении уравнений для соседней линии. Линии выбираются в произвольной последовательности. В настоящей работе сначала рассматривается линия вдоль оси х сразу над нижней границей (] — 1). Затем все параллельные ей линии перебираются снизу вверх до верхней границы (] — п) и в обратном направлении. После этого выполняется прогонка по всем линиям вдоль оси Оу слева направо и обратно. Такой алгоритм обеспечивает быстрое «проникновение» информации о граничных условиях во внутреннюю часть расчетной области.

В трехмерии добавляется направление вдоль оси г.

4. Пример расчета

В качестве примера, был выполнен расчет динамики температурного поля грунтов основания здания рис. 1 с размерами в горизонтальном сечении 20 м х 8 м. Координаты вершин здания (в метрах), соответственно, равны (0,0), (10,0), (10,4) и (0,4). Размеры области В выбирались равными хЬ =21 м, уЬ =15 м, гЬ =17 м. В расчетной области

выделялось три литологических слоя (как на рис.1). Физические характеристики литологических слоев, необходимые для расчета, приведены в табл. 1 [3].

Температура на нижней границе расчетной области г — 2~Ь (см. рис. 1) полагалась равной и(X, у, 2Ь, ?) — в0 — -1,4. На верхней границе области г — 0 выделялось 3 зоны

(см. рис.2). В зоне 2 (где находится здание) физические характеристики остаются постоянными во времени. Условия в зонах 1 и 3 изменяются по интервалам времени -месяцам года (см. табл. 2). За начало отсчета принимается 01 января текущего года: момент начала эксплуатации здания. Через год эти условия повторяются, т.е. условия на верхней границе являются периодическими функциями времени с периодом Т — 8760 часов (1 год=8760 часов). Начальное распределение температуры по глубине представлено в табл. 3.

Таблица 1

Физические характеристики литологических слоев

Номер Характеристика

слоя Р, Cd, Wtot, Ww (u), доли

кг/м3 ккал/(кг-°С) доли ккал/(м-ч-°С) ккал/(м-ч-°С) u, °С Ww (U)

-0,3 0,14

1 1390 0,22 0,25 1,3 1,15 -1,0 0,12

-10,0 0,08

-0,3 0,12

2 1520 0,24 0,22 1,55 1,55 -1,0 0,08

-10,0 0,05

0 0

3 1500 0,23 0,27 2,35 2,15 0 0

0 0

Решение краевой задачи (2) - (6) методом контрольного объема осуществлялось с помощью программы, написанной на языке Visual Fortran. При этом область D разбивалась сначала на блоки, а затем блоки дополнительно разбивались на более мелкие контрольные объемы (ячейки). Размеры блоков в метрах по направлениям осей следующие:

по оси x : 3, 3, 2, 2, 2, 3, 3, 3;

по оси у : 2, 2, 2, 3, 3, 3; по оси г: 2, 2, 2, 2, 3, 3, 3. Всего блоков 8x6x7=336. Разбиение блоков на ячейки выполнялось следующим образом. В каждом блоке по всем трем направлениям х, у иг вводилась равномерная сетка: длины блоков по х, у и г делились на 5 частей. Таким образом, в результате число контрольных объемов по х, у и г составило, соответственно, 40, 30 и 35, т.е. общее число элементов в расчетной области (это же число равняется числу внутренних расчетных точек в области В) равняется 40x30x35=42000. Минимальный безразмерный шаг по пространству АИт{п = 0.4/Ь=0.019, А^шах = 0.6/Ь=0.029; безразмерный шаг по времени Аt =7.3/Т=0.00083. Безразмерная полуширина сглаживания £-функции А =2.0/и=0.08.

Таблица 2

Физические характеристики на верхней границе расчетной области

Зоны Характеристика Месяц

I II III IV V VI VII VIII IX X XI XII

1 Температура в,°С -24,6 -23,8 -19,2 -9,5 3,7 13,1 18,0 12,4 4,6 -4,6 -16,4 -22,4

Термическое сопротивление я, м2-ч-°С/ккал 2,6 2,7 2,8 3,0 1,7 - - - - 1,3 2,1 2,3

Коэффициент конвективного теплообмена а, ккал/(м2-ч-°С) 13,8 12,3 12,3 12,0 12,1 12,6 13,4 13,0 14,3 13,4 13,5 12,7

2 Температура в, °С в =20°С= сопвг

Термическое сопротивление я, м2-ч-°С/ккал Я=1,8= сот(

Коэффициент конвективного теплообмена а, ккал/(м2-ч-°С) а =20.0= сотг

3 Температура в, °С -24,6 -23,2 -19,2 -9,5 3,7 13,1 18,0 12,4 4,6 -4,6 -16,4 -22,4

Термическое сопротивление я, м2-ч-°С/ккал 4,4 4,4 4,7 5,0 2,82 - - - - 2,14 3,57 3,59

Коэффициент конвективного теплообмена а, ккал/(м2-ч-°С) 13,8 12,3 12,3 12,0 12,1 12,6 13,4 13,0 14,3 13,4 13,5 12,7

Таблица 3

Начальное распределение температуры по глубине

Глубина, 1 3 5 7 9,5 12,5 15,5

м

°С -0,1 -0,4 -1,0 -1,2 -1,4 -1,4 -1,4

5. Результаты расчетов

Ниже представлены результаты моделирования динамики температурного поля грунта в расчетной области. На вертикальных линиях, проходящих через точки наблюдения M1(1.5,1), М2(9,5) и М3(19.5,13.5) (см. рис.2), выведены на печать в определенных точках по глубине зависимости температуры и в этих точках от времени ? (рис.8-10). На рис.8 представлены изменения температуры за 5 лет на вертикальной прямой, проходящей через точку М1: 1 - в точке на глубине 1 м (г = 1); 2 - в точке на глубине 3 м (г = 3); 3 - в точке на глубине 5 м (г = 5) ; 4 - в точке на глубине 7 м ( г = 7 ) и 5 - в точке на глубине 15.5 м ( г = 15.5). На рис.9-10 приведены зависимости от температуры через 5 лет. На рис.9 приведены зависимости в точках вертикальной прямой, проходящей через точку М2 (см. рис.2): 1- г = 1, 2 - г = 3, 3 - г = 7, 4 -г = 15.5, а на рис.10 - в точках вертикальной прямой, проходящей через точку М3 (рис.2): 1 - г = 1, 2 - г = 3, 3 - г = 7, 4 - г = 15.5.

Рис. 8. Зависимость температуры от времени в точках наблюдения (1.5,1,1), (1.5,1,3), (1.5,1,5), (1.5,1,7), (1.5,1,15.5) (кривые 1-5, соответственно)

Рис. 9. Зависимость температуры от времени в точках наблюдения (9,5,1), (9,5,3), (9,5,7),

(9,5,15.5) (кривые 1-4, соответственно)

- 1-2= 3-г= :1м. 2 - 7м, 4-2 =3м =15.5м

- 1 [ н

- V \ \

- \ \

— I -Чг- | | | 1 4

0 1 2 3 4 5 Iгоды

Рис. 10. Зависимость температуры от времени в точках наблюдения (19.5,13.5,1), (19.5,13.5,3), (19.5,13.5,7), (19.5,13.5,15.5) (кривые 1-4, соответственно)

Рис. 11. Распределение температуры по разрезу 5X1 — 5X1' при г = 1, 3, 5, 15.5

(кривые 1-4, соответственно)

12

и,°С

-4

1 / -г=1м

2 - 2=Зм

5 -2=5М

\ 4-2= 7м

2 \ 5 - 2= 15.5м

3

5 ----

.......

у,м

О 4 8 12 16

Рис. 12. Распределение температуры по разрезу 5у1 — 5уГ при г = 1, 3, 5, 7, 15.5

(кривые 1-5, соответственно)

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

на указанных глубинах в разрезе 5х1 — £хГ (уравнение разреза: у — 1) а на рис. 12 -зависимости и — и(у) на разных глубинах в разрезе 5у1 — 5уГ (уравнение разреза: х — 5). Приведенные на рис.11,12 зависимости температуры соответствуют по времени 5 годам после начала эксплуатации здания.

6. Заключение.

Результаты моделирования позволяют наблюдать воздействие, оказываемое зданием на грунт, и анализировать возможность проявления опасных криогенных процессов, в результате такого воздействия.

Так например, по рис. 11 видно, что грунт под зданием через 5 лет после начала эксплуатации будет талым, в результате чего произойдет его осадка и здание окажется в аварийной ситуации. Рис. 8 демонстрирует, что на небольших глубинах, осадка произойдет существенно раньше, чем через 5 лет. Например, на глубине 2=3 м (кривая 2) грунт расстаит уже в первый год эксплуатации здания.

Рис. 9-10 демонстрируют годовую периодичность процессов пучения и осадки грунта вблизи здания (точка М2) и на некотором расстоянии от него (точка М3). Периодический процесс пучения и осадки в точке М2 может привести к разрушительным последствиям.

Резкий спад температуры вдоль оси x на кривой 1 рис. 11 (обусловленный выходом за границу отапливаемого здания) указывает на возможность образования наледей.

По рис.8 видно, что в области грунта, расположенной непосредственно под зданием, по прошествии нескольких лет температура становится положительной даже на больших глубинах (примерно через 2 года - на глубине 5 м, соответствующей кривой 3, через 5 лет - на глубине 7 м, соответствующей кривой 4). Если, например, грунт незасоленный песчаный (температура начала замерзания грунта и* — 0° С), а кровля вечномерзлых грунтов расположена выше глубины 5 м, то возможно таяние многолетнемерзлых грунтов с образованием просадки примерно через 2 года после начала эксплуатации здания. По рис.11 видно, что распределение температуры вдоль оси x на фиксированной глубине под зданием не является постоянным (см., например, кривую 3, соответствующую глубине 5 м), что означает неравномерность просадки и перекос здания.

Благодарности. Работа выполнена при финансовой поддержке ОАО «Газпром»

[1] Алексеев С.И. Основания и фундаменты. Часть 12. СПб. , 2007г. 113с.

[2] СНиП 2.02.04-88. Основания и фундаменты на вечномерзлых грунтах. М., 1990.

[3] РСН 67-87. Инженерные изыскания для строительства. Составление прогноза изменений температурного режима вечномерзлых грунтов численными методами. М., 1987.

[4] Самарский А.А., Вабищевич. Вычислительная теплопередача. М.: Едиториал УРСС, 2003. 784 с.

[5] Патанкар С.В. Численные методы решения задач теплообмена и динамики жидкости. М: Энергоатомиздат, 1984. 152 с. (Перевод с английского. Patankar S.V. Numerical heat transfer and fluid flow. New York: Hemisphere Publishing Corporation, 1980)

[6] Патанкар С.В. Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах. М.: Издательство МЭИ, 2003. 312 с. (Перевод с английского. Patankar S.V. Computation of conduction and duct flow heat transfer. Innovative Research, Inc., 1991)

[7] Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 432 с.

electronic scientific and technical periodical

SCIENCE and EDUCATION

_EL № KS 77 -3()56'J..VaU421100025. ISSN 1994-jMOg_

Modeling the dynamics of temperature field soil base of the building in permafrost

77-30569/274059 # 12, December 2011

Glasko A.V., Fedotov A.A., Sidnyaev N. I., Hrapov P.V., Mel'nikova Yu.S.

Bauman Moscow State Technical University

[email protected] [email protected] si dnyaev@yandex. ru [email protected] [email protected]

We simulated the dynamics of temperature field soil base of the building, installed on a ground floor in a continuous permafrost in areas with seasonal freezing and thawing of the soil. The model is based on boundary-value problem for the heat equation with phase transitions. The equation describes the variation with time of the temperature field in a rectangular area of ground under the building and take into account the thermal effect of heated buildings on the ground and seasonal changes in climatic conditions (air temperature). The problem was solved numerically using method of control volume, the solution is implemented in the environment of Compaq Visual Fortran. The simulation results allow to analyze the possibility of dangerous cryogenic processes, resulting in a building in an emergency situation (heave, settlement, etc.).

Publications with keywords: the heat equation, phase transitions, the base plate, permafrost, cryogenic processes

Publications with words: the heat equation, phase transitions, the base plate, permafrost, cryogenic processes

Reference

1. Alekseev S.I., Bases and foundations, Part 2, SPb., 2007, 113 p.

2. Samarskii A.A., Vabishchevich, Computational heat transfer, Moscow, Editorial URSS, 2003, 784 p.

3. Patankar S.V., Numerical heat transfer and fluid flow, New York, Hemisphere Publishing Corporation, 1980.

4. Patankar S.V., Computation of conduction and duct flow heat transfer, Innovative Research, Inc., 1991.

5. Samarskii A.A., Gulin A.V., Numerical methods, Moscow, Nauka, 1989, 432 p.