CC BY
34НАУКА И ОБРАЗОВАНИЕ, 2016, №2
УДК 550.343; 551.345
Численное моделирование термического состояния криолитозоны в условиях меняющегося климата
П.П. Пермяков*, С.П. Варламов**, П.Н. Скрябин**, Ю.Б. Скачков**
*Институт физико-технических проблем Севера им. В. П. Ларионова СО РАН, г. Якутск **Институт мерзлотоведения им. П. И. Мельникова СО РАН, г. Якутск
Рассмотрена методология повышения достоверности численного моделирования термического состояния мерзлых грунтов. Достоверность моделирования состоит из погрешностей трех видов: погрешности модели за счет неполного учета факторов; погрешности за счет статистической неоднородности влияющих (входных) факторов; инструментальной погрешности за счет неточного измерения влияющих факторов и вычислительного алгоритма. Достоверность численного моделирования термического режима криолитозоны можно повысить, используя математическую модель средней сложности, уточняя входные параметры. Основными входными факторами являются теплообмен с атмосферой, метеоусловия, снежный и растительный покровы, влажность и литологиче-ский состав отложений. Используя натурные замеры, скорректированы входные данные. Проведен численный расчет с учетом ежемесячной динамики снежного покрова, среднесуточной температуры воздуха, неоднородности грунтового массива. Показана удовлетворительная сопоставимость полученных результатов расчетов с экспериментальными данными на трех ключевых участках Центральной Якутии. Предложенную методику численного моделирования следует использовать для расчёта термического режима многолетнемерзлых грунтов в условиях меняющегося климата.
Ключевые слова: мерзлые грунты, математическая модель, численный эксперимент, сравнение, натурные наблюдения.
Numerical Simulation of the Thermal State of Cryolithozone in a Changing Climate
P.P. Permyakov*, S.P. Varlamov**, P.N. Skryabin**, Yu.B. Skachkov**
*Melnikov Permafrost Institute SB RAS, Yakutsk **Larionov Institute of Physical and Technical Problems of the North SB RAS, Yakutsk
A methodology to improve the reliability of numerical modeling offrozen soils thermal state is considered. The accuracy of the simulation consists of three types of errors: a model error due to incomplete accounting factors; error due to statistical heterogeneity affecting input factors; instruments error due to inaccurate measurement of the influencing factors and computational algorithm. The accuracy of the numerical simulation of thermal regime of cryolithozone can be increased using a mathematical model of middle complication, making more accurate input parameters. The main input factors are heat exchange with the atmosphere, weather conditions, snow and vegetation covers, moisture and lithological composition of the sediments. Using field measurements the input data was adjusted. Numerical calculation taking into account the monthly dynamics of snow cover, average daily air temperature and heterogeneity of soils mass was carried out. Satisfactory comparability of obtained numerical results with the experimental data on three key sites of Central Yakutia was shown. The proposed method of numerical modeling should be used to calculate a thermal regime ofperennial frozen grounds in conditions of a changing climate.
Key words: frozen soils, mathematical model, numerical experiment, comparison, field observations.
ПЕРМЯКОВ Петр Петрович - д.ф.-м.н., в.н.с., e-mail: permyakov2005@mail.ru; ВАРЛАМОВ Степан Прокопье-вич - к.г.н., с.н.с., e-mail: vsp@mpi.ysn.ru; СКАЧКОВ Юрий Борисович - к.г.н., с.н.с., e-mail: vsp@mpi.ysn.ru; СКРЯБИН Павел Николаевич - к.г.н., в.н.с., e-mail: vsp@mpi.ysn.ru.
Введение
Оценка эволюции термического режима грунтовых оснований инженерных сооружений в условиях меняющегося климата требует достоверности и надежности математической модели. Расчет температурного режима осуществляется с помощью детерминистских математических моделей с использованием современных численных методов и программных средств. Сущность математического моделирования состоит в замене исходного объекта его «образом» - математической моделью - и в перспективе её завершение с помощью вычислительно-логических алгоритмов, реализуемых на компьютере. Такое математическое представление процесса является универсальным методом научного познания в области региональной и инженерной геокриологии.
Целью представленной работы является усовершенствование методологии количественного геокриологического прогноза на основе современных информационных технологий.
Материалы и методология исследования
Оптимизация процесса математического моделирования и расчета теплового состояния криолитозоны является одной из самых актуальных современных проблем, независимо от метода реализации. Достоверность прогноза AF предусматривает минимальное отклонение результатов математического моделирования от натурного замера и состоит из погрешностей трех видов [1-3]:
AF =
Т - ТЭ
= AMF + А ^ + А ^ ^ тп,
где Т - результаты математического моделирования; Т - данные натурного замера объекта; AмF - погрешность модели за счет неполного учета факторов; AнF- погрешность за счет статистической неоднородности влияющих (входных) факторов; AиF - погрешность за счет инструментальных измерений влияющих факторов и вычислительного алгоритма.
Качественный характер изменения погрешности модели AмF в зависимости от числа учитываемых факторов выглядит следующим образом [2]. Если учтены все влияющие факторы (практически бесконечное число факторов п, т. е. п^<х>), то модель будет идеальной (AмF=0). Если, наоборот, не учтен ни один фактор (п = 0) -гипотетически такой прогноз можно осуществить методом игральной кости - эта ошибка может быть сколь угодно большой, т. е. AмF^<x>. Погрешность модели зависит от структуры модели. Процесс определения вида модели называется структурной идентификацией [3].
Погрешность неоднородности AHF возрастает с увеличением числа факторов п. Его возрастание примерно пропорционально 4п [2]. Таким образом, для любой функции F можно считать, что погрешность неоднородности возрастает с увеличением числа учитываемых факторов п. Каждый учитываемый фактор (параметр) следует определить, исходя из физических свойств исходного объекта. Такой процесс называется параметрической идентификацией модели [3].
Современные измерительные средства и компьютерные вычислительные алгоритмы имеют минимальную погрешность. Характер изменения инструментальной и алгоритмической погрешностей AиF идентичен погрешности неоднородности, но даже значительно меньше и обычно ими можно пренебречь [2-3].
Качественный вид зависимости погрешности достоверности AF от числа учитываемых факторов неоднородности п показан на рис. 1.
Рис. 1. Зависимость точности от погрешности учитываемых факторов неоднородности
Из рис. 1 видно, что погрешность достоверности прогноза AF имеет минимум. Минимум кривой суммарной ошибки AF находится при небольших значениях учитываемых факторов п. Это говорит о том, что требование повышения достоверности расчета делает целесообразным применение менее точных моделей с меньшим числом учитываемых факторов п [2]. Приведенные выше соображения, казалось бы, свидетельствуют о нецелесообразности разработки сложных моделей процессов с привлечением большого числа учитываемых факторов. Однако анализ показывает, что это не так. Выше нами предполагалось, что все - и простые, и сложные модели являются физически правильными, т. е. не искажают реальную природу явлений. Такое предположение не всегда соответствует действительности, так как любая модель в какой-то
А
мере упрощает природу развития процессов. Поэтому роль сложных моделей с большим числом факторов важна с точки зрения более полного выявления физических особенностей и значимости влияния отдельных факторов. Наряду с этим, следует учитывать, что практическое прогнозирование основывается на моделях средней сложности. Более подробно этапы информационно-компьютерного решения различной сложности задач прикладной и вычислительной математики приведены в работе [1]. При этом следует учесть машинное время и объем занимаемой памяти.
Используя предложенную методологию, проведен численный эксперимент по расчету температурного режима грунтов на трех ключевых участках Центральной Якутии. Для достоверности прогноза процесса промерзания - протаи-вания на деятельном слое грунта выбрали известную математическую модель Стефана с фазовым переходом и провели сравнительный анализ с натурными замерами. Стефановский подход относится к классу средней сложности и соответственно занимает небольшой объем памяти и вычислительного времени [4].
Классическая постановка задачи (известная задача Стефана) состоит из уравнения теплопроводности в талой (От) и мерзлой зонах (Ом), и двух граничных условий на подвижной поверхности (£ (г)) раздела фаз [4-6]:
рмсм (Тм) =МК (Тм ^), дг
дТ
рТсТ (ТТ) —- = й\\{Кг (ТТ )УТТ), дг
*еОм, (1)
х еОТ, (2)
дТТ
Аг—--Км
дп
дТм
м
дп
= ,
х
Тт = Тм = Тф,
й- Ь
(3)
(4)
На верхней границе исследуемой области задается граничное условие Ш-го рода, для которого тепловой поток в грунт определяется конвективным теплообменом между поверхностью пород и воздухом: д Т
- А— = а эф (Т)в - Тп) д п
Здесь Т - температура, °С; г - время, с; х = (хх, х2, х3) - пространственная координата; с - удельная теплоемкость, Дж/кг; р - плотность, кг/м3; К - коэффициент теплопроводности, Вт/(м-К); V -оператор Гамильтона; дТ / дп, Vn - производная и скорость перемещения по-
(5)
верхности £(г) по направлению нормали п, м/с; Тф - температура кристаллизации поровой влаги, °С; Твз, Тп - температура воздуха и поверхности грунта соответственно; аэф - коэффициент конвективного теплообмена с воздухом, Вт/(м2-К); индексы: м - мерзлый; т - талый.
Численная реализация поставленной задачи (1-5) осуществлена итерационной неявной разностной схемой [5].
Обычно при прогнозе температурного режима грунтов используются осредненные средне-многолетние метеорологические данные. Для более точного прогноза нами учитываются ежемесячные и ежесуточные метеорологические данные. В качестве исходных данных приняты ключевые материалы круглогодичных геотемпературных наблюдений на экспериментальных площадках стационаров Чабыда (лиственничник) в 20 км к юго-западу от г. Якутска и Якутск (сосняк и разнотравный луг) [8-9]. Площадки при одинаковом литологическом составе на стационаре Якутск отличаются лишь распределением влажности и растительными покровами.
Исходные информации о температурном режиме воздуха за летний Твл(г) и зимний (Твз(г)) периоды аппроксимируем стандартной тригонометрической функцией [7]:
пО„
Т (г)=-
Твл (г) =
2гз
2г
Э1П
V гз У
(г^
-Э1П
г
V л У
(6)
(7)
где Овз, Овл - суммы отрицательных и положительных температур воздуха, К; г, гл - продолжительность, соответственно, зимнего (среднесуточная отрицательная температура воздуха) и теплого периодов года, с; г - время, с.
Аналогично аппроксимируем коэффициент конвективного теплообмена аэф(г) от высоты снежного покрова (табл. 1).
Обсуждение результатов
Основными факторами, определяющими термическое состояние грунтов слоя годовых теп-лооборотов, являются теплообмен с атмосферой, метеоусловия, снежный и растительный покровы, влажность и литологический состав отложений.
Т а б л и ц а 1
Термическое сопротивление снежного покрова, Вт/(м2К)
Высота снега, м 0,0 0,05 0,10 0,15 0,20 0,25 0,30 0,35 0,40 0,45 0,50 0,55
а , Вт/(м2-К) 12,3 5,44 1,84 1,15 0,67 0,52 0,42 0,37 0,31 0,27 0,21 0,19
Лесной покров снижает приток солнечной энергии к поверхности почвы, что способствует уменьшению её нагревания и протаивания. Поэтому в летний сезон лес оказывает почти всегда охлаждающее воздействие на почву. На стационаре Чабыда экспериментальная площадка расположена в лиственничнике с брусничным покровом на супесчаных почвах [810]. Литологический разрез и теплофизические характеристики грунтов представлены в табл. 2.
Численный расчет глубины протаивания, среднегодовых температур на подошве деятельного слоя и годовых теплооборотов охватывает период 1987-2014 гг. (рис. 2 и табл. 3). При
этом учитываются ежемесячная динамика снежного покрова, среднесуточная температура воздуха (данные метеостанции г. Якутска), а также неоднородность грунтового массива.
На стационаре Якутск в сосновом лесу влажность поверхностного почвенного горизонта бывает несколько большей, чем на разнотравном лугу, что объясняется хорошей водоудерживаю-щей способностью лесной постилки. Почвенные же горизонты, расположенные в средней части деятельного слоя (интервал 0,6-1,2 м), оказываются чаще всего более иссушенными [9]. Данный горизонт играет важную роль в распределении температуры по глубине массива.
Т а б л и ц а 2
Теплофизические характеристики грунтов
Грунт Толщина, м Р , кг/м3 ' Влажность, % К, Вт/(м-К) К,, Вт/(м-К) Ст , кДж/(м3-К) См , кДж/(м3-К)
Лиственничник (стационар Чабыда)
ПРС* 0,0-0,08 580 23 0,34 0,43
Супесь 0,08-1,0 1430 18 1,57 1,87 2198 1696
Суглинок 1,0-1,5 1,5-2,3 1400 31 33 1,50 1,64 1,90 2,05 2930 3224 2261 2407
Супесь 2,3-18 1650 25 2,19 2,68 2847 2114
Разнотравный луг, сосняк (стационар Якутск)
ПРС* 0-0,3 780 20 0,42 0,51
Суглинок 0,3-0,6 1420 10 1,2 1,43 1694 1397
Песок 0,6-8,4 1420 25 1,95 2,45 2467 1724
Песок с включениями гальки 8,4-20 1650 20 2,44 3,07 2867 2003
Примечание: ПРС - почвенно-растительный слой
Среднегодовые температуры грунтов (1) на подошве деятельного слоя (2 м) и годовых теплооборотов (10 м), оС
Ландшафт Год Глубина, м Методы
Натурный, t 0С Численный, t 0С
Брусничный лиственничник (стационар Чабыда) 2007 2 -1,3 -0,9
10 -1,9 -1,9
2011 2 -2,9 -1,9
10 -2,3 -2,3
Разнотравный сосняк (стационар Якутск) 2007 2 -2,5 -1,0
10 -3,0 -2,8
2012 2 -3,0 -2,6
10 -3,2 -3,3
Разнотравный луг (стационар Якутск) 2007 2 -1,3 -0,4
10 -1,8 -1,5
2012 2 -1,2 -1,1
10 -1,6 -1,8
Т а б л и ц а 3
1987
1990
1993
1996
Время, год
1999 2002
2005
2008
2011
2014
130
Рис. 2. Динамика глубины протаивания грунтов в лиственничнике на стационаре Чабыда: 1 - натурные наблюдения; 2 -расчет
Аналогичные результаты за 25 лет представлены для площадки в сосновом лесу (рис. 3). Расчеты дают заниженные результаты глубины протаивания по сравнению с натурными замерами.
На рис. 4 представлена сравнительная динамика глубины протаивания деятельного слоя на разнотравном лугу. Видно, что результаты чис-
ленных расчетов удовлетворительно совпадают с натурными наблюдениями.
Приведены сравнительные результаты численных расчетов температурного поля с данными натурных наблюдений в 2007, 2011-2012 гг. (табл. 3). Самым теплым является 2007 год. Результаты численного расчета температуры грунтов на подошве слоя годовых теплооборотов
1987
1990
1993
Время, год
1996 1999 2002 2005
2008
2011
2014
140
1 —■— 2 -в—3
Рис. 3. Динамика глубины протаивания грунтов в сосняке на стационаре Якутск: 1 - натурные наблюдения; расчетные: 2 - нулевая изотерма; 3 - изотерма -0,3 оС
1987 1990 1993 1996 1999 2002 2005 2008 2011 2014
150
2160 и -170
то |180
*®190
£200
210
220
Л/ \ в эемя , го/ * V.......................
к-' 1 »-■ г" кЛ*1—*-з /
Рис. 4. Динамика глубины протаивания грунтов на разнотравном лугу на стационаре Якутск: 1 - натурные наблюдения; 2 - расчет
ШЕСТАКОВА, СПЕКТОР, ТОРГОВКИН, СПЕКТОР
удовлетворительно совпадают с натурными данными, а незначительное отклонение наблюдается на подошве деятельного слоя.
Выводы
Сущность предложенной методологии состоит в замене исходного объекта его «образом» — математической моделью, и в дальнейшем её изучении с помощью реализуемых на компьютере вычислительно-логических алгоритмов. Сама постановка вопроса о математическом моделировании какого-либо объекта состоит из трех этапов: модель — алгоритм — программа.
Построение модели или качественное описание исследуемого процесса должны пройти этапы структурной (конструктивной) и параметрической идентификации факторов неоднородности.
Динамику термического режима криолитозо-ны более реально отражает математическая модель средней сложности.
Предложенную методику численного моделирования следует рекомендовать для расчета термического режима грунтов слоя годовых тепло-оборотов в условиях меняющегося климата.
Литература
1. Об использовании резервов оптимизации вычислений в компьютерных технологиях решения задач прикладной и вычислительной математики с требуемыми значениями характеристик качества / М.Д. Бабич, В.К. Завдирака, В.А. Людвиченко, И.В. Сергиенко // Журн. вычис лит. математики и мат. физики. - 2010. - Т. 50, № 12. - С. 2285-2295.
2. Гречищев С.Е. Оптимальные модели природных процессов погрешности и надежность // Анализ и оценка природных рисков в строительстве: материалы конференции. - М. - 1997.
- С. 14-16.
3. Алифанов О.М. Идентификация процессов теплообмена летательных аппаратов. - М. : Машиностроение, 1979. - 216 с.
4. Пермяков П.П. Идентификация параметров математической модели тепловлагопереноса в мерзлых грунтах. - Новосибирск : Наука, 1989.
- 86 с.
5. Самарский А.А. Устойчивость разностных схем / А.А. Самарский, А.В. Гулин. - М. : Еди-ториал УРСС, 2005. - 384 с.
6. Павлов А.В. Актуальные аспекты моделирования и прогноза термического состояния криолитозоны в условиях меняющегося климата / А.В. Павлов, Г.З. Перльштейн, Г.С. Типенко // Криосфера Земли. - 2010. - Т. 14, № 1. - С. 3-12.
7. Пособие по прогнозу температурного режима грунтов Якутии / Г.М. Фельдман, А.С. Тетельбаум, Н.И. Шендер, Р.И. Гаврильев. -Якутск : Изд-во ИМЗ СО РАН, 1988. - 240 с.
8. Варламов С.П. Температурный режим грунтов мерзлотных ландшафтов Центральной Якутии / С.П. Варламов, Ю.Б. Скачков, П.Н. Скрябин. -Якутск : Изд-во ИМЗ СО РАН, 2002. - 218 с.
9. Павлов А.В. Теплофизика ландшафта. -Новосибирск : Наука, 1979. - 284 с.
10. Скрябин П.Н. Межгодовая изменчивость теплового режима грунтов района Якутска / П.Н. Скрябин, С.П. Варламов, Ю.Б. Скачков. -Новосибирск : Изд-во СО РАН, 1998. - 144 с.
Поступила в редакцию 04.02.2016
УДК 624.131:551.345
Опыт составления инженерно-геологической карты Республики Саха (Якутия) масштаба 1:1 500 000
А.А. Шестакова, В.Б. Спектор, Я.И. Торговкин, В.В. Спектор
Институт мерзлотоведения им. П.И. Мельникова СО РАН, г. Якутск
Инженерно-геологическая карта территории Республики Саха (Якутия) охватывает площадь около 3 млн. км2, составляющую 1/5 часть территории России. На карте показаны грунтовые и гео-
ШЕСТАКОВА Алена Алексеевна - к.г.н., н.с., aashest@mail.ru; СПЕКТОР Владимир Борисович - д.г.-м.н., г.н.с., vbspector@mpi.ysn.ru; ТОРГОВКИН Ярослав Ильич - к.г.н., н.с., torgovkin@mpi.ysn.ru; СПЕКТОР Валентин Владимирович - к.г.н., доцент, зав.лаб., vvspector@mpi.ysn.ru.
