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

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

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

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

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

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

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

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

дТ 1 д Л „ дТ ) <

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

дт г" дг ^ дг )

X = 0, г = Я, 0 < т < ттах; (2)

Т (г, 0) = Т0(г) , 0 < г < Я, т = 0; (3)

с = с (Т) =

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

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

п0 ^пс

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

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

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

Т(Я, т) = Тп(т), 0 < т < ттах, (4)

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

dT

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

r^o dr v ’

(5)

(6)

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

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

T(r,, т)= Tэ (т), i = imr. (7)

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

Пт тшах / \2

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

i=l 0

тшах 2

+ |Р(т) (т((т))-Тф Т(т))) , 0 <)(т)R, (8)

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

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

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

Во(т) =тАт [(т + 2Дт)++ - 4 (т + Дт)+ + 6(т)+ -(Дт)

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

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

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

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

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

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

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

£ 1 р1 (т) V (Т _Тэ ) йт + Г р (т) V (т ( (т))-Тф ) )

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 г.