CC BY
28Учитывая частное решение ^(г) и принимая во внимание, что нам нужна действительная часть для краевой задачи (1), будем иметь
T = Re < e-iLUt
2 ^ Jo (rpneiфп) . ^ ~ 7
I ^ Jo (aPne^) sin T ZJ Í{Z,) SÍn f dZ' + T°
n= 1 o
Выделение действительной части не представляет больших сложностей, так как все функции зата-булированы.
Таким образом, данное решение сведено фактически к вычислению интегралов f^ f(z') sin ^ z' dz'. Во многих практически интересных случаях эти интегралы вычисляются аналитически.
СПИСОК ЛИТЕРАТУРЫ
1. Карслоу Г., Егер Д. Теплопроводность твердых тел. М.: Наука, 1964.
2. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. М.: Наука, 1976.
3. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1972.
4. Ват,сон Г.Н. Теория бесселевых функций. Ч. I. М.: ИЛ, 1949.
Поступила в редакцию 15.04.2013
УДК 536.21
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ПРОМЕРЗАНИЯ
ВЛАЖНОГО ГРУНТА Б. П. Лазарев1
Численно реализован в трехмерной постановке метод решения нестационарных задач переноса тепла во влажном грунте с учетом фазовых переходов в области отрицательных температур. В численной модели учитывается физическая нелинейность работы грунта. Для описания фазовых переходов в промерзающем грунте использовалась модель, состоящая из четырех характерных температурных зон, в которых основные термодинамические параметры грунта изменяются по определенным закономерностям. Численное моделирование задачи производилось при помощи метода конечных элементов для дискретизации краевой задачи в пространстве и разностной схемы для дискретизации по времени. Полученные численные решения сравниваются с аналитическими решениями задачи Стефана.
Ключевые слова: промерзание, теплопроводность, нелинейность, фазовый переход, теплоемкость, метод конечных элементов, задача Стефана.
A numerical method is implemented for the solution of the unsteady three-dimensional heat transfer problem in moist soils with consideration of phase transitions in a region of subzero temperatures. The corresponding numerical model takes into account the physical nonlinearity of the soil. For the description of phase transitions in frozen soils, a model that consists of the four characteristic temperature zones where the basic thermodynamic soil parameters vary-according to certain laws is used. The numerical simulation is carried out using the finite element method to discretize the boundary value problem in space and a finite difference scheme to discretize it in time. The obtained numerical solutions are compared with the analytical solutions to the Stefan problem.
Key words: freezing, thermal conductivity, nonlinearity, phase transition, heat capacity, finite element method, Stefan problem.
1. Перенос тепла. Рассмотрим процесс промерзания влажного массива грунта, в котором фронт промерзания распространяется от его поверхности. Процесс переноса тепла во влажном грунте описывается при помощи нестационарного уравнения теплопроводности
V • (А(Т) • VT) = С(Т) •pd-T-Q, (1)
1 Лазарев Борис Петрович — асп. каф. механики композитов мех.-мат. ф-та МГУ, e-mail: lazarev-bpQyandex.ru.
где A(T) — тензор теплопроводности грунта (Вт/м - К), C(T) — тензор удельной теплоемкости грунта (Дж/кг - К), pd — плотность сухого грунта (кг/м3), Q — мощность внутренних источников тепла (Вт/м3) [1, 2]. В этом уравнении тензоры теплопроводности и теплоемкости являются функциями температуры и изменяются в процессе промерзания грунта.
Влажность грунта явным образом не входит в уравнение теплопроводности. Известно, что суммарная влажность грунта при отрицательных температурах является функцией температуры: Wtot = Wtot(T)• Тогда можем переписать тензор теплоемкости в следующем виде:
C(T) = Cy(T) + L0dW^T), (2)
где Lo — скрытая теплота фазового перехода вода-лед. Первое слагаемое этого выражения Cv(T) является объемной теплоемкостью грунта (талого или мерзлого в зависимости от фазы, в которой находится грунт), а второе слагаемое — количеством тепла, которое необходимо затратить, чтобы перевести всю влагу, содержащуюся в единице объема грунта, из жидкого состояния в твердое.
Суммарную влажность грунта можно представить в виде Wtot (T) = Wfr(T) + Ww (T), где Wfr(T) — влажность грунта, обусловленная наличием свободной воды; Ww(T) — влажность грунта, обусловленная наличием связанной воды (содержание незамерзшей воды). Содержание незамерзшей воды для глинистых грунтов может быть найдено из выражения
Ww(T) = Kw(T) - Wp, (3)
где Wp — влажность на границе раекатывания, Kw(T) — коэффициент содержания незамерзшей воды в мерзлых грунтах [2]. Значения коэффициентов для различных типов грунтов приводятся в [3]. Из-за практически полного отсутствия связанной воды в песчаных грунтах содержание в них незамерзшей воды может быть принято равным нулю.
2. Моделирование процесса промерзания с учетом фазовых переходов. Для моделирования процесса промерзания/оттаивания была выбрана модель И.А. Цытовича (см. [3]), которая содержит четыре характерные температурные зоны и учитывает изменение теплофизических характеристик грунта при изменении температуры. Для каждой из зон выписан соответствующий вид функции теплоемкости (с учетом уравнения (2)) и теплопроводности.
I зона. Грунт находится в талом или переохлажденном состоянии в диапазоне температур T > Tf где Tbf — температура начала замерзания свободной воды, функции теплоемкости и теплопроводности принимаются постоянными и равными теплоемкости и теплопроводности талого грунта соответственно:
Ci (T) = Cth = const; Ai (T) = Ath = const.
II зона. Происходит замерзание (оттаивание) свободной воды, содержащейся в грунтах. Эта зона характеризуется значительными фазовыми переходами в диапазоне температур Tfw <T < Tf гд е Tfw — температура начала замерзания связанной воды, функции теплоемкости и теплопроводности в этом случае имеют вид
^ m ~ (Cth-Cf)(rbf-r) , LoPd(W(T)-Ww(T)) WA-t ) — Ь-th--~-~--1--~-~-,
1 bf — 1 f 1 bf — 1 fw
л и=л (d=Ath -
III зона. Происходит замерзание (оттаивание) связанной воды, содержащейся в грунтах. Количество незамерзшей воды начинает приближаться к количеству прочносвязанной воды в грунте в диапазоне температур Tf < T < Tfw, где Tf — температура замерзания практически всей рыхлосвязанной влаги в грунте. Функции теплоемкости и теплопроводности в этом случае имеют вид
CAT) = CAT) + ср„ = с, + <С"--°'><:-г'> + ь0Р1
AIII = A(T) = Af —
Tbf — Tf urd dT
(Ath — Af )(T — Tf)
Tbf — Tf
IV зона. Грунт находится в практически мерзлом состоянии в диапазоне температур T < Tf. Функции теплоемкости и теплопроводности принимаются постоянными и равными теплоемкости и теплопроводности мерзлого грунта соответственно:
Civ (T) = Cf = const; Aiv (T) = Af = const.
3. Краевые условия. Будем рассматривать граничные условия трех видов:
1) граничное условие первого рода в случае, если известна температура на некоторой поверхности Et-Тогда
T (x, t) = To(t), x e Et ;
2) граничное условие второго рода в случае, если известен тепловой поток на поверхности Eqi. Тогда
-X ■ VT ■ n = q°(t), x е Eqi,
где q0(t) — плотность теплового потока (Вт/м2), n — нормаль к поверхности Eqi;
3) граничное условие третьего рода в случае, если на поверхности Eq2 происходит конвективный теплообмен со средой. Тогда
-X ■ VT ■ n = a(T - Ta), X е Eq2,
где a — коэффициент теплоотдачи (Вт/м2 ■К), Ta — температура окружающей среды (К), n — нормаль
Eq2 t
поля температур во всем объеме грунта:
T (x, t) = th (x), t = 0.
Тогда уравнение (1) вместе с начальными и граничными условиями, а также с учетом функций теплопроводности и теплоемкости, в которые подставлено выражение для содержания незамерзшей воды в мерзлом грунте (3), даст нам уравнение нашей модели.
4. Дискретизация. Для дискретизации задачи в пространстве используется метод конечных элементов [4]. С этой целью запишем вариационную постановку, соответствующую дифференциальной. Умножим уравнение (1) скалярно на пробную функцию v е Gt, где gt является функциональным множеством, которое состоит из элементов пространства Соболева W2,, удовлетворяющих граничным условиям первого рода. Полученное уравнение проинтегрируем по объему V:
J vV ■ (X(T) ■ VT) dV = j vC(T) ■ pd ■ TdV - J vQ dV. (4)
y V V
Преобразуем левую часть уравнения (4) при помощи теоремы Грина с учетом граничных условий. Упростим выкладки, ограничившись граничными условиями второго рода, и положим температуру на границе Et постоянной. Тогда получим
У Vv ■ X(T) ■ vTdV + J vqo dE = j vC (T) ■ pd ■ TdV - J vQ dV. (5)
V S,1 V V
Функция T(x,t) е gt, которая удовлетворяет уравнению (5) для любой функции v(x,t) е gt и для каждого t > 0, будет являться обобщенным решение исходной краевой задачи (1).
Рассмотрим Ж-мерное подпространство G^ е gt и соответствующую функцию TN е GN, которая является приближенным решением уравнения (5) при условии, что для всех vN е GN при любом t > 0 она удовлетворяет следующему уравнению:
У VvN ■ X(T) ■ VTNdV + j vNq0dE = J vNC(T) ■ Pd ■ TNdV - J vNQ dV. (6)
V Sq1 V V
Разобьем область V на конечные элементы Ve и запишем уравнение (6) для совокупности этих элементов:
/
Е yj VvNX(T)VTNdV j
vN qodE I = ^ IJ vN C (T )pd TN dV j - ^ lj vN QdV j . (7)
\Sq1e j e \Ve ) e \Ve
Выразим функции Т(х, Ь) и ь(х, Ь) и их производные по пространственным координатам через их значения в узлах элементов:
NN
TN = Е Ns(x)Ts = NsTs, T N = N%lTs, i = 1, 2, 3; vN = v.s(x)vs = Nsvs, vN = N%lvs, i = 1, 2, 3.
s=1 s=1
Подставляя эти выражения в уравнение (7), получим
//
AijNs,iNhj dV I Tsvh | + ^
Ve
J NhQodZ I vh I = ^ IN CpdNsNh dV I T svh I —
\Sqie / / e We
) ^
Это уравнение можно переписать в следующем виде:
^shTsvh — V Klhfsvh) = V FheVh
где Cesh = f CpdNsNh dV — матрица теплоемкости элемента, Kesh = f AijNs,iNh,j dV — матрица теплоем-
Ve Ve ' '
кости элемента, F^ = — f NhqodY< — J NhQ, dV — вектор-столбец правой части уравнения.
Sqie Ve
Учитывая, что пробная функция выбиралась произвольно, перейдем от локальной нумерации к глобальной и получим систему дифференциальных уравнений для всей области:
[C]T + [K]T = {F}. (8)
Для решения системы (8) использовался неявный метод Эйлера ввиду его абсолютной устойчивости. После ее дискретизации по времени при заданных начальных условиях приходим к следующей системе алгебраических уравнений:
([C]ra+1 + т[K]n+1) {Tn+1} = [C]n+1{Tn} + т{F}n+\
Где т — шаг по времени. Решая эту систему, получаем необходимое нам распределение температуры в момент времени tn+1.
При решении модельной задачи о промерзании грунта столбик грунта разбивался на конечные элементы с восемью узлами, которые представляют собой параллелепипеды.
5. Промерзание столбика грунта. Для проверки корректности работы программы результаты сравнивались с аналитическим решением одномерной задачи Стефана. Рассмотрим задачу о промерзании влажного столбика песчаного грунта с влажностью Wtot = Wfr = 0,2, так как для песчаного грунта можно полагать, что влажность грунта за счет связанной воды Ww = 0. На свободной поверхности грунта задана постоянная температура Tc = —5 °C. Температура грунта в начальный момент времени постоянна по глубине и равняется Th = 0°C. Температура начала промерзания Tbf = —0,2 °C. Общее аналитическое решение задачи приводится в работе Лыкова [1] в виде
„ x f Р
erf —= eric ■
2 y/aft , . , ,
Tf(x, t)=Tc + (Ты - Tc)-y- Tth(x, t)=Tu- (TH - Tbf)-w
erf -A= erfc;
где х — расстояние по глубине; а^ь) — коэффициент температуропроводности мерзлого и талого грунта соответственно; ей(ейс) — функция ошибок (дополнительная функция ошибок); в — коэффициент, который приближенно можно получить из выражения
о 2 • Af • Тс
Р =
То • • р • п • af
Графики распределения температуры в грунте в различные моменты времени приведены на рис. 1.
Температура, °С -3 -2
-1
X
к
ю >.
1
!» - ?
NN
4 ^
Рис. 1. Распределение температуры в промерзающем грунте по годам: 1 1-й год. 2 2-й год. 3 4-й год. 4 6-й год. Штриховой линией показаны соответствующие аналитические решения
Рис. 2. Распределение температуры при двустороннем промерзании грунта по годам: 1 1-й год.
2 2-й год. 8 4-й год. 4 5-й год. 5 6-й год. Штриховой линией показаны соответствующие
аналитические решения
Можно заметить, что численным решением хорошо аппроксимируется аналитическое решение, но первое дает более пологую кривую из-за того, что выделение теплоты фазового перехода происходит не мгновенно, а в интервале температур Ты < Т < Тн- Соответственно при уменьшении этого интервала численное решение будет приближаться к аналитическому.
Рассмотрим задачу о промерзании столбика грунта высотой 10 м, когда на обеих поверхностях задается температура Тс = —5 °С, в остальном условия задачи аналогичны первой. Графики распределения температуры в грунте в различные моменты времени приведены на рис. 2.
Численные результаты и в этом случае показывают хорошую сходимость с аналитическим решением. Таким образом, можно сказать, что численное моделирование процесса промерзания влажного грунта
при помощи метода конечных элементов в нелинейной постановке позволяет достаточно точно описывать процесс промерзания и изменения температуры в грунте.
Работа выполнена под руководством C.B. Шешенина.
СПИСОК ЛИТЕРАТУРЫ
1. Лыков A.B. Теория теплопроводности. М.: Высшая школа, 1967.
2. Карслоу Г., Егер Д. Теплопроводность твердых тел. М.: Наука, 1964.
3. Цытович H.A. Механика мерзлых грунтов (общая и прикладная). М.: Высшая школа, 1973.
4. Зенкевич О., Чанг И. Метод конечных элементов в теории сооружений и в механике сплошных сред. М.: Недра, 1974.
Поступила в редакцию 22.04.2013
