CC BY
17
МЕХАНИКА
УДК 517.956
ОБ ОДНОМ ПОДХОДЕ К ВОССТАНОВЛЕНИЮ КОЭФФИЦИЕНТОВ ПЕРЕНОСА
© 2009 г. А. О. Ватульян, С.А. Нестеров
Южный федеральный университет, Southern Federal University,
ул. Мильчакова, 8а, Ростов н/Д, 344090, Milchakov St., 8a, Rostov-on-Don, 344090,
dnjme@math.stedu.ru dnjme@math.stedu.ru
Предложен альтернативный подход идентификации коэффициентов переноса для неоднородного стержня. Обратная задача формулируется в пространстве трансформант. На основе метода линеаризации она сводится к поэтапному решению интегрального уравнения Фредгольма 1-го рода с гладким ядром.
Ключевые слова: обратная задача, реконструкция, уравнение теплопроводности, итерационный процесс, интегральное уравнение, линеаризация, коэффициент теплопроводности, теплоемкость, плотность.
Alternative approach of identification of transport factors for inhomogeneous rod is proposed. Inverse problem is formulated in space of transforms. On the basis of method of linearization inverse problem is reduced to step by step solution ofFredholm integral equation of the first kind with smooth kernel.
Keywords: inverse problem, reconstruction, heat conduction equation, iterative process, integral equation, linearization, heat conductivity factor, heat capacity, density.
Математическая модель теплопроводности с постоянными коэффициентами переноса - основа научно-технических расчетов теплового и тепломеханического состояния тел. При этом коэффициенты теплопроводности и температуропроводности определяются из простых макроэкспериментов. Для многих материалов составлены обширные таблицы. Большое внимание уделено слоистым композитам, а также материалам, свойства которых зависят от температуры [1]. Однако в последние годы в различных областях техники для замены слоистых композитов появляются функционально-градиентные материалы - композиты с непрерывным изменением физических свойств. Для их адекватного описания необходим отказ от гипотезы однородности среды. Эффективность использования неоднородных материалов зависит от точного знания их свойств после изготовления, для чего нужно поставить и решить коэффициентные обратные задачи. Наиболее распространенный метод их решения - сведение на этапе параметрической идентификации к задаче нахождения минимума функционала невязки [1], легко приводимой к конечномерной нелинейной проблеме. При этом применяются как градиентные методы, так и генетические алгоритмы. Решение прямых задач находят, как правило, при помощи методов конечных элементов или конечных разностей.
Коэффициентные обратные задачи теории упругости и теплопроводности по определению коэффициентов дифференциальных операторов являются нелинейными некорректными проблемами, для которых в последние годы интенсивно развиваются как математические аспекты исследования, так и вычислительные средства [2-4]. При решении такого типа задач наи-
большие трудности возникают при построении нелинейных операторных соотношений, связывающих искомые и измеряемые в эксперименте функции. Чтобы обойти эту сложность, в [5, 6] предложен новый подход. Прямая задача сведена к решению интегрального уравнения Фредгольма 2-го рода, обратная - интегрального уравнения Фредгольма 1-го рода. В [6] для восстановления модуля сдвига неоднородного по толщине слоя обратная задача формулировалась в пространстве трансформант Фурье, причем использование метода линеаризации приводит к нахождению поправок неизвестных функций из поэтапного решения интегрального уравнения Фредгольма 1 -го рода с гладким ядром. Оказалось, что этот подход эффективен и для реконструкции теплофизических характеристик неоднородного стержня. Решение прямой задачи теплопроводности сводится к интегральному уравнению Фред-гольма 2-го рода [7]; в [8] приведены метод решения и результаты реконструкции коэффициента теплопроводности, в [9] представлены результаты восстановления коэффициентов переноса для экспоненциальных, квадратичных и кусочно-однородных законов неоднородности. Все расчеты по реконструкции неизвестных коэффициентов переноса проводились при точном знании остальных характеристик без наложения шумов. Однако на практике, как правило, неизвестны все коэффициенты, и при реализации процедуры идентификации нужно учитывать погрешность экспериментальной информации.
В настоящей работе на основе метода линеаризации построен итерационный процесс, представлена система интегральных уравнений Фредгольма 1-го рода в пространстве преобразований Лапласа для нахождения
поправок одновременно ко всем коэффициентам переноса. Проведены вычислительные эксперименты, учитывающие погрешность входных данных.
Постановка прямой задачи теплопроводности неоднородного стержня
Рассматривается задача о распространении тепла в прямолинейном неоднородном изотропном стержне длиной /, на одном торце которого х = 0 поддерживается нулевая температура, на другом х = / действует постоянный тепловой поток д0 при нулевых начальных условиях.
Прямая задача состоит в том, чтобы, зная коэффициенты теплоемкости с(х), плотности р(х) и теплопроводности к(х), найти температуру точек стержня Т как функцию времени.
ох ох dt
71(0,0 = 0, (х,0) = 0 .
ВТ
ох
(1)
(2) (3)
Сведение задачи к интегральному уравнению Фредгольма 2-го рода
Точные аналитические решения начально-краевой задачи (1)-(3) могут быть получены только для некоторых типов неоднородности, в частности для некоторых степенных и экспоненциальных функций. В общем случае законов изменения коэффициентов переноса решение может быть построено лишь численно.
Перейдем к безразмерным параметрам
Г К > к, '
c0
Po
knt kn ~
т =-9—; W{z, t) = Tx-2-\ r(z) = c(z)p(z).
q0l
12
— (k (z)—) = r (z)—. dz dz dr
W( 0, т) = 0 ,
~ dW
-k(l) — (\T) = l,
dz
(4)
Щг,0) = 0 . (5)
Применим к системе (4) преобразование Лапласа. С учетом начального условия (5) имеем краевую задачу с параметром
d ' dW
— (k(z)—) = pr(z)W(z,p); dz dz
Щ0,р) = 0,
(6) (7)
dz p
(8)
или каноническую систему обыкновенных дифференциальных уравнений относительно трансформант: <Ш 1 =
dz
dQ dz
k (z) = pr(z)W(z, p) ,
(9)
с0р0Г
Здесь . с о . р() - некоторые характерные коэффициенты теплопроводности, теплоемкости и плотности соответственно.
Для безразмерных параметров краевая задача (1)-(3) примет вид
д ~ . dW^ , ч dW
где р - параметр преобразования Лапласа; W(2, р), Q(z, р) - безразмерные трансформанты температуры и теплового потока.
Проинтегрируем систему дифференциальных уравнений (9) по 2 на отрезке [0,2]: 2 \ _
О Щ)
^ 2 ^ 0(г, р) = | РЖ + С2 •
о
Граничное условие (8) преобразуется к виду
0(1 ,р) = -~. (Ю)
Р
Найдя постоянные интегрирования из граничных условий (7), (10), исключив из системы (9) Q(2, р), поменяв порядок интегрирования в двойных интегралах, получим интегральное уравнение Фредгольма 2-го рода:
~ 1 ~
Щг, р) = \К(2. £ рЩ + /(г, р) . (11)
о
к{1,е,р)=-рг(£) / —
о к{П) РОЩ)
Для решения интегрального уравнения (11) используем метод коллокаций, заменив интегралы их приближенными значениями по квадратурной формуле трапеций, для чего введем равномерное разбиение отрезка [0,1] на п отрезков точками г, = Дг^—1
(/ = 1.л+ 1), где /\г = 1 /я - шаг разбиения. Конечномерная аппроксимация интегрального уравнения (11) представима в матричном виде: А(р) = Е-В(р), А(р)\¥(р) = Б1(р).
Здесь Е - единичная матрица; В(р) - матрица, полученная дискретизацией ядра интегрального уравнения на основе квадратурной формулы трапеций; Б(р), "(р) - векторы значений правой части и трансформант температуры в узловых точках.
Метод тестировался в случае однородного стержня, для которого известно точное решение; тестирование показало высокую точность метода.
Нахождение температуры на торце стержня
Для нахождения безразмерной температуры на торце г = 1 стержня необходимо решить интегральное уравнение (11) и осуществить обратное преобра-
f(z,p) =--1-
зование Лапласа: 1Г(1. г) -- j W{\, р)ертdr, или
2 777
У Г А
переходя к размерной температуре,
pküt
qoi r+>f
= / W(\,p)ecp»1 dt.
1лк01 r_jœ
(12)
- При £°
d ~ dW 1 ~
f±)=Pr„-1(zW„-1(z,P):
dz dz
W„_l(0,p) = 0,
dW 1 1
dz p
при e1
d ~ dW d ~ dW 1 d (kn_ !( z) ^ К d (K (z) W1 ) = dz dz dz dz
(15)
(16)
= p(rn-\(z)fF„ (*> P) + JV„(0,P) = 0,
dW ~ dW ,
dz
dz
(17)
(18)
Функция W (1, р) имеет полюс в нуле и счетное число отрицательных полюсов первого порядка. Интеграл (12) вычисляется по теории вычетов. Для нахождения вычетов в полюсах используется формула, вытекающая из разложения W (1, р) в ряд Лорана в окрестности простого полюса р0 [6].
Постановка обратной задачи
Будем считать, что температура на торце стержня известна в любой момент времени
г1(/,о = /1(о,/>о. (13)
Рассмотрим постановку обратной задачи: из условия (13) найти г (г) и к (г). Она относится к классу коэффициентных обратных задач [1], нелинейна и в общем случае приводится к решению системы нелинейных операторных уравнений.
При такой постановке обратную задачу можно формулировать в пространстве трансформант.
Будем искать к (г) и г (г) в два этапа.
На 1-м определяется нулевое приближение в классе линейных функций: к0 аг + Ъ , г (г) = сг + ё.
На 2-м уточняются законы к (г) и г(г) на основе метода линеаризации в окрестности некоторого известного состояния. Рассмотрим состояние И/и+1 (1 ,р),
кп+\(г), гп 11 (г). для которого известно
\¥т(\,р) = \¥п_1(\,р) + £\¥п(\,р).
Для этого положим
кп+! (г) = кп_х (г) + (кп (г), ги+1 (г) = гп_х (г) + егп (г), К+1(г,р) = Жп_1(2,р) + £Жп(2,р) . (14)
Здесь е - формальный малый параметр; ]¥т - безразмерная трансформанта торцевой температуры, вычисленная для точного закона.
Подставляя (14) в (6)-(8) и приравнивая коэффициенты при одинаковых степенях параметра е , получим краевые задачи:
Умножая (15) на Wп, (17) - на , вычитая одно из другого, интегрируя результат по частям на отрезке [0, 1], с учетом граничных условий (16), (18) получим интегральное уравнение Фредгольма 1-го рода относительно неизвестных функций:
0
1
1
dz
)zdz + p\rn(z)Wt_l(z,p)dz = о
(19)
Р
Одного уравнения вида (19) для нахождения поправок двух неизвестных функций недостаточно. Нужно получить 2-е уравнение, исходя из другой задачи с новыми граничными условиями следующего вида:
ох ох <3/
дТ0 ох
T2(l,t) = 0,
(20)
Т2(х, 0) = 0,
72(0,0 =/2 (О-
Производя действия, аналогичные описанным выше для задачи (1)-(3), (13), получим 2-е интегральное уравнение:
} кп (Z)(^±)2dz + p)rn (zWl, (z, p)dz = о dz о (21)
= -фт-вп_1\ Р
где 6(z,p) - безразмерная трансформанта температуры, которая находится из решения уравнения Фредгольма 2-го рода:
~ 1 ~ в (z, p) = JA(z, £ р)в (£ p)dÇ + g(z, p) ,
A(z,Ç,p) = -pr{£) }
dr/
(л 1 f
g(z,p) =--—
Pzl
тт {г,£}к(Г}) РгЩ)
При решении обратной задачи необходимо получить значения торцевых температур ТТ (I, /) и Т2т (0,/) в ходе проведения двух тепловых экспериментов при граничных условиях (2) и (20) соответственно. В работе натурный эксперимент заменен вычислительным. Значения функций Тт (I, /) и Тт (0, /) находятся при точных заданных значениях коэффициентов переноса.
Определение поправок из системы интегральных уравнений (19), (21) является некорректной задачей и требует регуляризации. Применена дискретная реализация метода Тихонова и осуществлен следующий
I j yLj : отобра-
принцип выбора набора параметров /] ' жение множества р е |,+ао ^ в отрезок т е [).1
С
отрезка
осу-
2 2
ществлено по правилу т =р равномерного разбиения
2>1 'P >
построение
re |,1-<Г
о
(<5е <¡..0,5^ точками 4 ' вычисление Д/ при помощи обратного отображения набора 4 // •
После нахождения решения кп (2) и гп (2) строятся новые приближения кп ¡\(г) = кп | (г) I ки(г) и -ги-1(2) + ги(2); таким образом осуществляется итерационная схема. Критерием выхода из итерационного процесса является условие на функционал
Л 2 2 невязки: Jn = Д(Г1й -Ч\т) + С/2п -Т2т) \сЛ < ст .
с
Здесь [с, ё] - отрезок времени, на котором вычисляется невязка; Т^п (1,1) и Т2п (0,I) - вычисленные торцевые температуры при приближенных значениях коэффициентов переноса на п-й итерации. Отметим попутно, что поиск параметров а, Ь , с и ё при поиске нулевого приближения осуществляется из условия минимизации функционала невязки; а - параметр, характеризующий качество решения (на практике выход из итерационного процесса осуществляется обычно при стабилизации этого функционала).
Численные результаты
На языке Фортран 90 составлена программа, реализующая представленный выше подход. Решение прямой задачи показало, что с ростом t Т(х,0 выходит на стационарное решение; поэтому при экспериментальной проверке необходимо проводить измерения торцевой температуры до выхода на этот уровень. Проведена серия вычислительных экспериментов по восстановлению коэффициентов переноса неоднородного стержня в классах степенных, экспоненциальных, тригонометрических и кусочно--однородных функций. Исследовано влияние на точность результата и характер итерационного процесса количества разбиений п и т, начального приближения, монотонности функций, возмущения входных данных.
Вначале проведены вычислительные эксперименты по восстановлению одной из характеристик по известной другой. Так, для восстановления коэффициента теплопроводности приходим к необходимости решать уравнение
1 ~ (1\¥ 1 о 1 ~ ~
/*„(*)(—= (22)
0 аг р
Для восстановления г(2) приходим к уравнению
1 1 ~
¡гп (2)Ш1_Х (г,р)сЬ = —(1ГТ- ). (23)
о р
Во всех расчетах принято: количество разбиений в квадратурной формуле п = т = 20, | с. с! \ = [100,120]с,
£/„ = -100000 Вт/м2, к0 = 65 Вт/(м - град), / = ОД м,
с0 - 670 Дж/(кг - град), р0 — 3100 кг/м3.
На начальном этапе вычислительных экспериментов коэффициенты переноса имели один и тот же функциональный закон изменения. Выяснено, что дальнейшее увеличение количества разбиений мало
снижает погрешность восстановления, но сильно увеличивает затраты машинного времени. Немонотонные и кусочно-однородные функции восстанавливаются гораздо хуже монотонных. Восстановление функции
г(2) происходит с большей погрешностью, чем к (2). Начальное приближение находилось в классе линейных функций, минимизирующих функционал невязки. Если искать начальное приближение в виде постоянной, то увеличивается количество итераций и растет погрешность восстановления. Поскольку в реальных экспериментах входные данные определяются с некоторой погрешностью, то с целью проверки устойчивости метода к возмущениям входных данных проведены вычислительные эксперименты при наличии аддитивных шумов, начиная с величины в 2 %. Выяснено, что восстановление г(2) невозможно даже при наличии 2 % погрешности исходных данных. Реконструкция коэффициента теплопроводности к ( 2) дает удовлетворительные результаты даже при 4%-м шуме для монотонных и 2%-м для немонотонных функций. При этом погрешность восстановления не превышает 11 %. В случае восстановления г(2) важно обратить внимание на ядро интегрального уравнения (23). Оно обращается в нуль при г = 0, поэтому следует ожидать плохого восстановления г( 2) на торце г = 0. Кроме того, вычислительные эксперименты показали большую погрешность восстановления к(£) на торце г = 1, что также обусловлено свойством ядра уравнения (22). Поэтому для улучшения качества восстановления желательно знать значения реконструируемых коэффициентов на торцах. Результаты расчетов показали, что при известных торцевых значениях погрешность восстановления значительно снижается. Если коэффициенты переноса имели разные законы неоднородности, реконструкция происходила значительно хуже и за большее число итераций.
Проведены вычислительные эксперименты по восстановлению одновременно двух коэффициентов переноса. Результаты показали, что восстановить две функции удается не для всех законов, но коэффициент теплопроводности в большинстве случаев восстанавливался. Хорошие результаты при реконструкции получены, когда коэффициенты переноса являются экспоненциальными функциями либо одна из функций - постоянна или линейна. В этом случае можно применить итерационный процесс для улучшения восстановления. Процедура реконструкции в этом случае устойчива при наличии 2-4 % аддитивных шумов. Они мало влияют на реконструкцию коэффициента теплопроводности, значительно больше на г(2) .
На рис. 1-3 представлены примеры восстановления коэффициентов переноса при известных торцевых значениях. При этом сплошной линией изображен точный закон, точками - восстановленный. На рис. 1 показан пример реконструкции функции г (г) = соз(г) при невозмущенных входных данных. В расчетах такой же закон неоднородности принимался и для коэффициента теплопроводности. Параметр регуляризации -
а -10 5. Погрешность восстановления не превысила 6 %. На рис. 2 представлен результат восстановления функции г(г) - 2г - 2г +1 при таком же законе для
_п
к (г). Параметр регуляризации - а = 10 . При этом
на 8-й итерации значение невязки - 1.33 -10 4. а погрешность восстановления не превысила 16 %. На рис. 3 приводится пример реконструкции коэффициента теплопроводности к(г) = е 1 (такой же закон имеет и г( 2)) при восстановлении одновременно двух функций при входных данных с 2%-м зашумлением в пространстве трансформант; при этом начальное приближение выбрано в виде г(:) = к(г) = 1 - 0,7г .
г
1
0,9 0,8 -0,7 -0,6 -0,5
0,2
0,4 0,6 Рис. 1
0,8
1
0,9 -0,8 -0,7 -0,6 -0,5
0,2 0,4 0,6 Рис. 2
0,8
1
0,9 -
0,8 -
0,7
0,6
0,5 -\
0,4
0,3
0 0,2 0,4 0,6 0,8 z
Рис. 3
Выводы
Предложен эффективный подход к решению коэффициентных обратных задач теплопроводности для неоднородного стержня. Решена прямая задача путем сведения к интегральному уравнению Фредгольма 2-го рода, разработан метод расчета температуры на торце стержня на основе теории вычетов. При известной торцевой температуре в любой момент времени поставлена в пространстве трансформант и решена на основе метода линеаризации обратная задача о восстановлении коэффициентов переноса.
Построена итерационная схема, основанная на численном решении интегрального уравнения Фредгольма 1-го рода. Сформулирован критерий выхода из итерационного процесса. Проведены вычислительные эксперименты по реконструкции коэффициентов уравнения теплопроводности.
Работа выполнена при частичной поддержке Южного математического института (г. Владикавказ).
Литература
1. Алифанов О.М., Артюхин Е.А., Румянцев С.В. Экстре-
мальные методы решения некорректных задач. М., 1988. 288 с.
2. Денисов А.М. Введение в теорию обратных задач. М.,
1994. 206 с.
3. Самарский А.А., Вабищевич П.Н. Численные методы
решения обратных задач математической физики. М., 2004. 480 с.
4. Ватульян А.О. Обратные задачи в механике деформи-
руемого твердого тела. М., 2007. 224 с.
5. Бочарова О.В., Ватульян А.О. Обратные задачи для уп-
ругого неоднородного стержня // Изв. вузов. Сев. Кавк. регион. Естеств. науки. 2008. № 3. С. 33-37.
6. Ватульян А.О., Сатуновский П.С. О реконструкции
модуля сдвига при анализе колебаний неоднородного по толщине слоя // Теоретическая и прикладная механика. Донецк, 2006. Вып. 42. С. 124-129.
7. Нестеров С.А. Об одном методе решения нестационар-
ной задачи теплопроводности для неоднородного стержня // Мат. моделирование и биомеханика в современном университете : тез. докл. Всерос. школы-семинара. п. Дивноморское, 28 мая - 1 июня 2007. Ростов н/Д, 2007. С. 63-64.
8. Нестеров С.А. О реконструкции коэффициента теплопро-
водности неоднородного стержня // Математическое моделирование и биомеханика в современном университете : тез. докл. Всерос. школы-семинара. п. Дивноморское, 2 -6 июня 2008. Ростов н/Д, 2008. С. 75-76.
9. Нестеров С.А. Об особенностях идентификации тепло-
физических характеристик неоднородных тел // Актуальные вопросы современной науки : сб. науч. тр. Ме-ждунар. интернет-конференции. Таганрог, 3-5 сентября 2008. М., 2008. С. 85-90.
Поступила в редакцию
21 января 2009 г.
0
1
r
z
0
