Идентификация параметров модели тепломассопереноса при техногенном загрязнении мерзлых грунтов

П.П. Пермяков

ИДЕНТИФИКАЦИЯ ПАРАМЕТРОВ МОДЕЛИ ТЕПЛОМАССОПЕРЕНОСА ПРИ ТЕХНОГЕННОМ ЗАГРЯЗНЕНИИ МЕРЗЛЫХ ГРУНТОВ

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, проект № 03-05-65408.

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

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

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

Для простоты изложения рассмотрим одномерное уравнение теплопроводности с учетом фазового перехода поровой воды [5]:

дТ 1 д Л „ дТ ) <

с —— =-----— \ Xг -г—I , 0 < г < Я, 0 < т < ттах;

дт г" дг ^ дг )

X = 0, г = Я, 0 < т < ттах;

Т

с = с (Т) =

= [сск + С^ + (Св - Сл ) (Т) + Ь ^^Т^) Р,

,W(T) - Ш X = X(T) = Xм + (Xт - Xм) Ш - Ш пс,

п0 ^пс

где с - объемная теплоемкость грунта, Дж/(м3-К); Т -температура, К; т - время, с; ттах - наибольшее значение ременной переменной, с; г - пространст енная координата, м; Я - координата массива, м; Xм, X,. - теплопроводности мерзлых, талых пород, Бт/(м-К);. Ь - объемная теплота фазового перехода, Дж/м3; Т0(г) - начальное распределение температур, К; сск, сл, св -удельные теплоемкости грунта, льда и воды, Дж/(кг-К); р - объемная плотность минерального скелета, кг/м3; Ш0, Шпс - начальная и прочносвязанная влажности, %;

Шнв(Т) - функция количества незамерзшей воды, которая в случае техногенного загрязнения грунта зависит и от засоленности Шнв(Т, С), %; V = 0, 1 - соответствуют декартовой и цилиндрической координатам.

Требуется восстановить одно из следующих граничных условий (при г = 0):

Т

- ^ rVX w=rVq(x)’ 0 < Т < TmaX’

dT

- lim rv X = rva(x) (T - Tc ),0 < т < Tm

r^o dr v ’

где искомым параметром и(т) может быть одна из функций {Т (т), <?(т), а(т)} ; Тп(т) - температура поверхно-сти, К; q(x) - плотность теплового потока, Вт/м2; а(т) - коэффициент теплоотдачи, Вт/(м-К); Тс - температура внешней среды, К.

Для восстановления искомого параметра нужны дополнительные замеры температуры внутри исследуемого образца в точках rt (0 < r1 < r2 < ... < rn < R):

T

Данную задачу сформулируем как задачу оптимального управления: найти функцию и(т)из минимума целевого функционала:

Пт тшах / \2

J (u) = S 1 Р,(т) (Т (г,, т)-Т (г,, т) +

i=l 0

тшах 2

0

|э(т) = sup (Т(м,т))= Тф(э(т)),

те[°.тшах]

где р,(т), р(т) - весовые множители с размерностью К-2-с-1; T(r,, т), T(r, т) - расчетная и замеренная температуры в i-й точке горного массива; Т(|(т)), Тф(|э(т)) - соответственно расчетные и заданные температуры на фронте фазового перехода |э(т); пт - число точек измерения.

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

Минимизация целевого функционала (8) проводится с помощью хорошо известных градиентных методов, характеризующихся эффективным началом итерационного процесса из некоторого «грубого» начального приближения и снижением скорости сходимости при приближении к минимуму [6]. Этим требованиям отвечают градиентные методы спуска, в частности методы скорейшего спуска и сопряженных градиентов. В этом случае приближения вычисляются по формуле

Us+1 (т) = Us (т) - PsPs S = n .

Направление спуска р*(т) и коэффициент в* определяются в зависимости от методов:

1) скорейшего спуска

р* = Ы'(«,). Р*: Ы (и+1 )= ™п Ы (и, -Р '(«,));

2) сопряженных градиентов

Л = 7'(и*)+Т*Р*-l, У0 = 0 У* =

Ц-7' («-1

Р*: 7 («+1 )= т™ 7 («-Р*Р*),

где ... г - норма в гильбертовом пространстве Ь2.

II N¿2

Формулы для градиента функционала соответственно неизвестным параметрам будут иметь вид

д7 д7

д7

- 4(т - Дт)+ + (т- 2Дт)+ ] (т-Ч)+ = [тах (0, т-Ч]3;

Дт = тт / т0 - заданный интервал времени. Решение сопряженной задачи у = у (г, т):

+Азу,0 < г < Я,0 <т<т.

дт дг2 дг2

Xгl

дг

= 0, г = Я, 0<т<тт, у(г,тт) = 0,0 <г<Я

У г- =У

Xr,

дг

-XV

= 2р(т)(т - Тэ), і = 1

0<т<ттах,

П- =п*,

xгv дг

дг

-*г"

дг ^ г 0 <т<т

= 2 р

где

. . . дX 1 д /.

А1 =X, А2 =^—+-----------—

1 ’ 2 дг --V дг V /

дг

зи^" ’амто+ч

д Т Хтах д

-±Г- = Г В„ (т)^-(уу)| йт,

дЬ„ 0 а дг 1Г=0

ух о

д Т Хтах

дг= 1ув(т) ^=ойт,

а 0

д Т Хтах

дЫ = Г Ва (Т)(С - Т )Уг=0 йТ ,0 <Т^Ттах.

а о

а = -1,0,1,..., т0 +1.

Искомым параметром является (т0 + 3)-мерный

Вектор и = Ь = (ь_1, Ьо , ..., Ьт0 +1 ):

т0+1

и

а=-1

где Ва(т) - заданный кубический сплайн по переменной т, который имеет вид

Ва (т) = Во (т-аДт),

Во

. = _!, 5Т )_ 5с _5Т Аз = гу дг\дТГ дг ) дТ Р дт .

Граничные условия при г = 0, соответствующие искомым параметрам, имеют вид

У(т) = 0,0 <т<ттах>

^ 1Т(У^ <т-ттах,

^ ^Г ^У ^) = _аУ’0 < т < ттах.

Величина шага р* методов скорейшего спуска и сопряженных градиентов находится по формуле

Пт ттах / \ ттах / / \ \

£ 1 р1

1=1 0 0

р* =-

Пт ттах / \

£ Р і (т) V2йт + І Р (т) V2 ( (т))^т

где V = У(г, т) - решение краевой задачи для приращения:

д 2V

дV „ с^— = А1 —2

1 ^,.2

дт

дг дXV

дV

+ А2 — + A3V,0 < г < Я,° < т < ттах,

дг = 0, г = Я,(0 <T<Tmax,

V(г, 0) = 0,0 < г < Я,

'І* = ^^1г-

дXV

дг

дкУ_

дг

,/ =1, пт ,0

дXV

дг

дг

0 < т < тт

Граничные условия в зависимости от искомых параметров имеют вид

V = ДТП (т), 0 <т<ттах,

Д д (т),0 < т < ттах,

дг

дг

= а (т) (т 1 - Тс )-а (т)V ,0 <т<ттах.

Численная реализация вышеуказанных краевых задач осуществляется методом конечных разностей. В каждом временном слое для численного определения температуры Т(г, т), сопряженной переменной у(г, т) и приращения V(г, т) будем иметь соответствующие трехдиагональные системы алгебраических уравнений: нелинейные для Т и линейные - у и V. Решение этих систем производится методом прогонки (в сочетании с итерациями по коэффициенту при расчете Т). Итерационный процесс минимизации останавливается по критерию невязки [5]. Аналогичным образом строятся алгоритмы идентификации остальных параметров (теплоемкости, теплопроводности, функции количества незамерзшей воды, коэффициента диффузии и т.д.).

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

2

і=1

+

+

+

+

г

+

Г;

+

ч

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

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

Рассмотрим два иллюстративных примера по расчету плотности тепловых потоков и функции количества незамерзшей воды песка при засолении его хлоридом натрия.

Для восстановления теплового потока рассмотрены три ключевых участка: открытый участок с нарушениями растительного и напочвенного покровов, сосновый и лиственничный леса (рис 1).

Рис. 1. Зависимость удельного теплового потока от времени:

1 - открытый участок с нарушениями растительного и напочвенного покровов; 2 - сосновый лес; 3 - лиственничный лес

Они отличаются друг от друга сложением грунта, растительным покровом, влажностью и температурным режимом. Были проведены многолетние круглогодичные наблюдения в разные времена [7, 8]. Отличительной особенностью лесного массива является его затененность деревьями, что влияет на формирования теплового режима и влагообменного процесса. Чем больше сомкнутость древостоя, тем меньше солнечной радиации проникает под полог леса. В период полного развертывания

листвы затененность деревьями снижает приход солнечной радиации. Особенно в листопадных лесах (лиственничный) проникает лишь 40-50 % радиации (июнь-август). В зимнее время плотности теплового потока в сосновом лесу и открытом грунте отличаются друг от друга из-за высоты и структуры снега.

На рис. 2 представлена пространственная зависимость функции количества незамерзшей воды от температуры и засоленности.

Рис. 2. Зависимость функции незамерзшей воды от температуры и засоленности

У чистого песка фазовый переход порового раство- сунка видно, что с повышением засоленности содержа-

ра происходит при температуре 273 К, а с увеличением ние незамерзшей воды увеличивается и соответственно

засоленности большое количество порового раствора резко меняются все теплофизические и массообменные

замерзает при температуре эвтектики 250,6 К. Из ри- характеристики загрязненного мерзлого грунта.

ЛИТЕРАТУРА

1. Фельдман Г.М. Методы расчета температурного режима мерзлых грунтов. М.: Наука, 1973. 254 с.

2. ЧистотиновЛ.В. Миграция влаги в промерзающих неводонасыщенных грунтах. М.: Наука, 1973. 144 с.

3. Crange B. W., Viskanta R., Stevenson W.H. Diffusion of heat and solute during freezing of salt solution // Int. J. Heat and Mass Transfer. 1976. Vol. 19. № 14. P. 373-384.

4. Ентов В.М., Максимов А.М., Цыпкин Г.Г. Образование двухфазной зоны при промерзании пористой среды. М.: ИПМ АН СССР. Препринт № 269, 1986. 56 с.

5. Пермяков П.П. Идентификация параметров математической модели тепловлагопереноса в мерзлых грунтах. Новосибирск: Наука, 1989. 86 с.

6. Алифанов О.М., Артюхин Е.А., Румянцев С.В. Экстремальные методы решения некорректных задач. М.: Наука, 1988. 288 с.

7. Павлов А.В. Теплообмен почвы с атмосферой в северных и умеренных широтах территории СССР. Якутск: Якутское кн. изд-во, 1975. 302 с.

8. Гаврилова М.К. Тепловой режим земной поверхности и грунтов поля и леса в период протаивания // Материалы VIII Всесоюзного междуведомственного совещания по геокриологии. Якутск, 1966. Вып. 4. С. 194-206.

Статья представлена отделом теплофизики и теплоэнергетики Института физико-технических проблем Севера (г. Якутск), поступила в научную редакци. «Кибернетика и информатика» 26 марта 2004 г.