Вход

Решение обратной задачи вихретокового контроля

Дипломная работа* по физике
Дата добавления: 15 апреля 2003
Язык диплома: Русский
Word, rtf, 2.8 Мб
Диплом можно скачать бесплатно
Скачать
Данная работа не подходит - план Б:
Создаете заказ
Выбираете исполнителя
Готовый результат
Исполнители предлагают свои условия
Автор работает
Заказать
Не подходит данная работа?
Вы можете заказать написание любой учебной работы на любую тему.
Заказать новую работу
* Данная работа не является научным трудом, не является выпускной квалификационной работой и представляет собой результат обработки, структурирования и форматирования собранной информации, предназначенной для использования в качестве источника материала при самостоятельной подготовки учебных работ.
Очень похожие работы

  1. Техническое задание Разработать алгоритм решения обратной задачи вихретокового контроля (ВТК ). Объе к том контроля (ОК ) являются проводящие немагнитные листы . Объекты контроля подверг а ются термообработке (закалка , отпуск ) или насыщению внешних слоев р азличными вещ е ствами , что приводит к изменению механических , а вследствие этого и электромагнитных свойств материала листа по глубине . Задача заключается в определении , в рамках допустимой погрешности , зависимости электропроводности (ЭП ) от глубины (Н ) в ОК для данного состояния . Метод контроля заключается в измерении определенного количества комплексных значений вносимой ЭДС на различных частотах с помощью накладного вихретокового преобразователя (НВТП ). Необходимо выбрат ь математическую модель задачи , способ аппроксимации искомого решения , рассмотреть алгоритм решения . Используя программную реализацию , исследовать поведение погрешности аппрокс и мации зависимости (Н ) от следующих факторов : 1. От величины приборной погрешности измерения ЭДС 2. От вида зависимости электропроводности от глубины (Н ) 3. От параметров аппроксимации решения 4. От диапазона частот возбуждения ВТП 2. Ана лиз технического задания. Основная задача вихретокового контроля с помощью накладных преобразователей с о стоит из двух подзадач : · Прямой задачи расчета вносимой ЭДС в присутствии немагнитного проводящего листа с произвольной зависимостью ЭП по глубине. · Обратной задачи нахождения зависимости ЭП как функции глубины в немагнитном пр о водящем листе по результатам измерений определенного количества комплексных знач е ний вносимой ЭДС. 2.1 Прямая задача ВТК Полагая зависимость ЭП от глубины известной проведем ее кусочно-постоянную а п проксимацию . Это позволяет свести исходную задачу к расчету ЭДС в многослойном листе , в каждом слое которого ЭП принимает постоянное значение . Как показано в работе [ 50 ] , подобная модель вполне адекватно описывает задачу и дает отлич ное согласование с результатами опытов. Рекуррентные формулы для произвольного количества слоев хорошо известны [ 1-5,36, 42,43,50-52 ]. Таким образом решение прямой задачи в рамках принятой модели затруднений не вызывает . 2.2 Обратная задача ВТК С математи ческой точки зрения обратная задача ВТК относится к классу некорректных задач [ 49 ] и ее решение неустойчиво т.е . при сколь угодно малой погрешности исходных данных ( набора измеренных вносимых ЭДС ) погрешность решения ( рассчитанных локал ь ных значений ЭП ) может быть сколь угодно большой , а одному набору измерений может отвечать много (формально бесконечно много ) распределений ЭП по глубине. При попытке расчета некорректной задачи как корректной , вычислительный процесс за счет неустойчивости сваливается в заведомо худшую сторону . В нашем случае это означ а ет получение распределения ЭП , которое , хотя и обеспечивает требуемое совпадение изм е ренной и вычисленной ЭДС , но является явно нереальным из-за осцилляций . Следует отм е тить , что амплитуда и частота осцилля ций распределения ЭП растут при увеличении числа независимых параметров аппроксимации ЭП ( коэффициентов полинома в случае полином и альной аппроксимации , количества узлов при сплайн-аппроксимации и т.д .). При наличии погрешности измерения вносимой ЭДС , пре вышающей на несколько порядков вычислительную погрешность и на практике составляющей не менее (0.5-1)% от измеряемого сигнала , ситуация значительно осложняется. Учитывая вышеизложенное для выделения из множества допустимых распределений решения , наиболее у довлетворяющего физической реальности , в алгоритмах решения обра т ной задачи необходимо использовать дополнительную априорную информацию . На практ и ке это реализуется введением некоторых критериев , позволяющих отличить решение , отв е чающее практике , от физиче ски нереального . Для решения обратной задачи ВТК предлагались три возможные стратегии [ 46 ] : Решение большого числа прямых задач и табуляция результатов для различных моделей . Измеренные данные с помощью некоторых критериев сравниваются с таблицей . Под ход очень экстенсивный и требующий проведения избыточного числа расчетов , поэтому на практике встречающийся редко. Условная минимизация невязки измеренных и расчитанных данных . Очень мощный и универсальный метод , широко распространен для решения обрат ных задач в различных областях техники [ 41,44,49 ] . Позволяет восстанавливать произвольное распределение ЭП по глубине (вообще говоря произвольное 3 D распределение ), но требуется довольно сложная процедура расчета. Аналитическое инвертирование ядра опе ратора и использование алгоритма , зависящего от ядра уравнения [ 46 ] . Потенциально самый малозатратный метод , однако как и все ан а литические , применим далеко не всегда. В нашем случае остановимся на втором подходе , поскольку он сочетает в себе униве р сальност ь , точность и относительную простоту реализации . В целом процесс решения обратной задачи сводится к итерационному решению прямой задачи для текущей оценки распределения ЭП и внесению изменений в эту оценку в соотве т ствии с величиной невязки. 2.3 Модель за дачи Приведем основные положения , на основе которых будет построена модель нашей з а дачи : · ОК представляет из себя находящуюся в воздухе проводящую пластину толщиной Н состоящую из N плоско-параллельных слоев толщиной b i . · В пределах каждого слоя удел ьная электропроводность имеет постоянное значение т.е . распределение по глубине аппроксимируется кусочно-постоянной зависимостью. · Возбуждающая и измерительная обмотки ВТП заменяются нитевид ными моделями . Следует отметить , что это предположение сказывается лишь на решении прямой зад а чи , а проведя интегрирование можно получить выражения для катушек конечных ра з меров. · Для численного моделирования реальных распределений ЭП применим пять типов а п проксимации : сплайном , кусочно-постоянную , кусочно-линейную , экспоненциальную и гиперболическим тангенсом . В процессе решения прямой задачи с их помощью в ы числяются значения в центральных точках слоев пластины. 2.4 Анализ литературы 2.4.1 Зарубежные методы решения Решению обратной задачи ВТК посвящен ряд работ в зарубежных изданиях . Следует отметить монографию [ 38 ] , в которой рассмотрены случаи импульсного возбуждения , а оп е рируют в частотной и временной областях напряженн остью электрического поля. Подход к решению квазистационарных задач рассмотрен в цикле статей [ 45-51 ] . Он о с нован на интегральной постановке задачи с помощью функций Грина [ 31-34,39 ] . Для илл ю страции рассмотрим решение обратной задачи ВТК согласно [ 49 ] . А . Прямая задача Определим функцию v(r)=( (r) - 0 )/ 0 , где (r) - произвольное распределение пров о димости , а 0 - ее базовая величина . Функция v(r) может представлять собой как описание произвольного распределения проводимости (в этом случае для удобства полагаем (r)= 0 вне некоторого ОК объема V , тогда v(r ) отлична от нуля только в пределах V ) так и некот о рого дефекта (для трещины v(r )=-1 внутри дефекта и равна нулю вне его ). Рассмотрим систему уравнений Максвелла в предположении гармонического возбу ж дения exp(-jwt ) и пренебрегая токами смещения : ( 2.4.1) где P(r)=[ (r)- 0 ] E(r)= 0 v(r) E(r) - может интерпретироваться как плотность дип о лей эффективного тока , причиной которого является вариация (r)- 0 . Решение уравнений Максвелла можно предс тавить в виде ( 2.4.2) где E i (r) - возбуждающее поле , а G(r|r ’ ) - функция Грина , удовлетворяющая уравн е нию G(r|r ’ ) + k 2 G(r|r ’ )= (r-r ’ ) , k 2 =-j 0 0 , (r-r ’ ) - трехмерная дельта-функция. Импеданс ВТП можно выразить как ( 2.4.3) где интеграл берется по измерительной катушке , J(r) - плотность тока в возбуждающей катушке . Применяя теорему взаимности импеданс можно представить через возбуждающее поле : ( 2.4.4) где интеграл берется по объему ОК. В . Обратная задача Пусть v(r ) - оценка истинной функции v true (r ) , Z obs (m) - измеренный импеданс ВТП в точке r 0 на частоте возбуждения , m=(r 0 , ) - вектор в некоторой области определения M , Z[m,v] - оценка величины Z obs (m) на основе решения прямой задачи. Определим функционал невязки измеренных и рассчитанных значений импеданса ВТП как : ( 2.4.5) Предположим , что для решения обратной задачи используется итерационный алгоритм типа метода спуска : v n (r)= v n-1 (r)+ s n (r) . Можно показать , что в случае метода на искорейш е го спуска итерация имеет вид : v n (r)= v n-1 (r)- F[ v n-1 (r) ] , где градиент функционала F[v] можно определить как : ( 2.4.6) где Re обозначает вещественную часть , * обозначает комплексную сопряженность. Требуемый в (2.4.6) градиент импеданса можно определить как : Z(r) = - 0 E(r) E * (r) ( 2.4.7) где E * (r) - решение уравнения ( 2.4.8) С . Аппроксимация при решении обратной задачи Пусть электропроводность моделируется с помощью конечного числа переменных (например узловых значений некоторой аппроксимации ), а вектор р состоит из этих пер е менных . Тогда выражение (2.4.7) принимает вид : ( 2.4.9) где ( Z ) j - j -ая компонента градиента импеданса. Значение j -ой компоненты градиента невяз ки (2.4.6) можно представить как : ( 2.4.10) Следует обратить внимание на то , что в случае дискретного пространства М (конечное число измерений ) интеграл в (2.4.10) заменяется суммой. С уче том приведенных преобразований итерация метода наискорейшего спуска прин и мает вид : p j n = p j n-1 - ( F n-1 ) j ( 2.4.11) где n - номер итерации. D . Пример применения В к ачестве примера рассмотрим функцию v(r ) в виде v(r )= c i i ( r ) , i=1,N , где i ( r ) - множество линейно независимых базовых функций с коэффи циентами c i . Рассматривая к о эффициенты c i в роли параметров аппроксимации ( c i = p i ) получим из (2.4.9) для компонентов градиента импеданса : ( 2.4.12) В случае проводящего ОК , состоящего из N параллельных слоев с проводимостью j распределение электропроводности по глубине можно представить с помощью функций Хевисайда H(z) как (z)= j [ H( z-z j ) - H( z-z j+1 ) ] . Подставляя в (2.4.12) базовые функции вида i ( z )= [H( z-z j )-H( z-z j+1 )] , получим оконч а тельное выражение : ( 2.4.13) Отметим основное преимущество такого решения . Несмотря на определенную сло ж ность вычислений при решении интегральных уравнений (2.4.2-2.4.8) для расчета градиента импеданса НВТП необходимо р ешить только две такие задачи. 2.4.2 Отечественные методы решения Подход , в значительной мере аналогичный работам [ 45-51 ] был предложен в работе [ 41 ] . Из-за небольшого объема в ней уделено недостсточное внимание вопросам практич е ской реализации , объяснены не все обозначения и не приведены результаты численного моделирования . В целом это значительно снижает практическую ценность статьи . Приведем основные положения этой работы. Прямая задача Пусть круговой виток радиусом а с током I находится в точке P=P s (r, ,z), (- , ) вблизи немагнитного ОК , занимающего область V . Пусть ОК обладает электрической пров о димость ю = 0 (Р ) являющейся произвольной функцией координат . Требуется по N изм е рениям величины э.д.с . определить как функцию координат точек P V . Причем i -ое изм е рение э.д.с . будем проводить на i -ом измерительном круговом витке с координатами P i =P i (r, ,z) i=1,N при неизменных частоте и расположении возб уждающего витка. В общем случае напряженность электрического поля Е определяется через векторный магнитный потенциал А , причем А = А 0 + А вн , где А 0 - возбуждающий , а А вн - вносимый п о тенциалы . (2.4.14) Вводя функцию Грина G(p,p 0 ) получим (2.4.15) При этом вносимая напряженность электрического поля E вн = -j A вн (2.4.16) Вносимая э.д.с ., наводи мая в i -ом витке (2.4.17) где функция Грина G(P,P 0 ) имеет вид (2.4.18) В дальнейшем рассмотрим случай , при котором V -пол упространство (r> 0 ,\ \ < ,z<0) с электрической проводимостью = 0 (Р ) , где (Р ) - произвольная функция координаты z . В этом случае выражение (2.4.17) примет вид (2.4.19) где k 2 =jw 0 0 . Для напряженности электрического поля Е (Р ) справедливо представление (2.4.20) где Е 0 (р ) - возбуждающее поле витка. После проведения серии из N измерений величины вн выражение (2.4.19) дает связь между вносимыми ЭДС i и (z)E(r,z) . Чтобы определить непосредственно = ( z ) , находим E(r,z) при известной функции (z)E(r,z) из (2.4 .20 ), после чего исключаем E(r,z) из известн о го. Обозначим x ( p ) =-k 2 (z)E(r,z) , а измеряемую совокупность ЭДС через Fi . Тогда (2.4 .19 ) можно записать в операторной форме как F = Px + (2.4. 21 ) где - погрешность измерения. Обратная задача Построим функционал Ф (х )= ||F-Px|| 2 + ||x-x 0 || 2 , где х 0 - некоторое известное близкое к искомому распредел ение , удовлетворяющее F 0 =Px 0 . Образуем вариацию функционала Ф (х ) , используя определение сопряженного оператора ( Px,y )=( x,P*y ) . Для нахождения мин и мума Ф (х ) приравняем его вариацию Ф нулю. Вводя вспомогательную функцию u=x-x 0 и учитывая F 0 =Px 0 проведем ряд преобраз о ваний . Искомое распределение ( z ) можно найти из равенства (2.4.22) где напряженность электрического поля в точке р для известного распределения (z) имеет вид (2.4. 23 ) (2.4. 24 ) Система алгебраических уравнений для определения коэффициентов С i имеет вид (2.4. 25 ) , j=1,N (2.4. 26 ) 3. Прямая задача ВТК для НВТП 3.1 Уравнение Гельмгольца для векторного потенциала Взаимодейст вие преобразователя с объектом контроля определяется системой уравн е ний Максвелла в дифференциальной форме [ 6 ] : (3.1) где H - вектор напряженности магнитного поля E - вектор напряженнос ти электрического поля B - вектор магнитной индукции - вектор плотности полного тока - вектор плотности токов проводимос ти - удельная электрическая проводимость - вектор плотности токов смещения D - вектор электрического смещения - вектор плотности токов переноса - вектор скорости переноса j стор - вектор плотности сторонних токов Дополним систему (3.1) уравнениями связи : B = 0 H (3 . 2) B = rot A (3.3) где 0 = 4 10 -7 - магни тная постоянная - относительная магнитная проницаемость A - векторный потенциал магнитного поля Преобразуем систему уравнений (3.1) с учетом следующих предположений [4] : * ОК неподвижен относительно электромагнитного по ля т.е . j пер = 0 * среда изотропна и ее параметры не зависят от напряженностей полей * воздействия синусоидальны * последовательность дифференцирования по времени и пространственным координатам можно изменять , а операция дифференцирования линейна (3.4.1) (3.4.2) Поскольку ротор градиента любого скаляра тождественно равен нулю , величину в скобках выражения (3.4.2) можно приравнять градиенту некоторого скаляра , например скалярного потенциала электрического поля : (3.5) Заменяя векторы напряженности магнитного и электрического поля в (3.4.1) через ве к торный потенциал магнитного поля получаем : grad div A - A = - 0 ( + j 0 ) ( grad + j A ) + 0 j стор (3.6) Откуда после очевидных преобразований следует : ( 3.7) где k 2 = 2 0 0 - j 0 (3.8) Поскольку векторный потенциал магнитного поля задан с точностью до градиента н е которого скаляра , а потенциал с точностью до постоянной величины , имеется в озможность положить значение величины в квадратных скобках выражения (3.7) равное нулю (так наз ы ваемая калибровка Лоренца ). В результате получаем уравнение Гельмгольца для векторного потенциала магнитного поля : (3.9) В дальнейших рассуждениях используем следующие предположения : · Поле НВТП квазистационарно в том смысле , что волновыми процессами в воздухе можно пренебречь . Это вполне оправдано т.к . размеры НВТП и ОК обычно много меньше дл ины волны в воздухе , а потери на излучение по сравнению с потерями в ОК малы. · В проводящем теле будем рассматривать только волновые процессы , обусловленные наличием параметров и т.е . токами смещения ( пропорциональными 0 ) как и в воздухе пренебрегаем . Легко показать , что это предположение справедливо не только для металлов , но и для полупроводниковых материалов с удельным сопротивлением до 50 [ Ом см ] . В этом случае выражение (3.8) принимает вид : . 3.2 Поле витка над многослойной средой Введем цилиндрическую систему координат ( r, , z ) . Пусть : · - ток , протек а ющий по нитевидной во з буждающей обмотке с р а ди усом R 1 , находящейся на расстоянии h от N - слойной среды · j стор = I ( z - h ) ( - R 1 ) Отметим , что в силу осевой симметрии системы В цилиндрической системе координат выражение (3.9) имеет следующий вид : (3.10) Применяя к (3.10) преобразование Фурье-Бесселя с ядром в виде функции Бесселя пе р вого порядка имеющее вид : получаем (3.11) Так как на поверхностях раздела слоев ОК должна сохранятся непрерывность танге н циальных составляющих векторов напряженностей магнитного и электрического поля , д о полняем уравнение (3.11) граничными условиями на поверхностях слоев ОК ( граничные условия одинаковы для А и А * ) : (3.12) (3.13) Решив уравнение (3.11) с учетом граничных условий (3.12-3.13) и применяя к решению обратное преобразование Фурье-Бесселя имеющее вид : пол у чаем для полупространства над ОК (3.14) где = ( , , ) - функция граничных условий . 3.3 Воздействие проводящего ОК на НВТП Для большинства инженерных расчетов можно использовать нитевидную модель обм о ток НВТП использованную в (п 3.2). При данном упрощении пол учаем : - напряженность электрического поля (3.15) - ЭДС наводимая в измерительной обмотке с радиусом R 2 и числом витков w 2 (3.16) Анализируя формулу (3.14) можно заметить , что первый интеграл представляет собой векторный потенциал создаваемый возбуждающей обмоткой , а второй интеграл - векторный потенциал вносимый ОК . В практике ВТК обычно анализируются вносимые параметры НВ ТП (напряжение , импеданс ) поэтому получим выражение для вычисления вносимого напряжение кругового трансформаторного накладного ВТП используя (3.15-3.16): (3.17) Подставляя выражение для вн осимого векторного потенциала (3.14) в уравнение (3.17) окончательно получаем : (3.18) где = 2 f - круговая частота тока возбуждения I 0 - магнитная постоянная w и , w в - числа витков измерительной и возбуждающей обмоток НВТП R = (R и R в ) - эквивалентный радиус НВТП K r = (R в /R и ) - параметр НВТП x - переменная интегрирования h * = (h и + h в )/2 - обобщенный зазор J 1 - функция Бесселя 1 рода 1 порядка m - функция граничных условий Функция граничных условий для m -слойного ОК с плоскопараллельными сл оями м о жет быть вычислена по рекуррентной формуле [2] : (3.19) где (3.20) (3.21) (3.22) th(z) - гиперболический тангенс m - относительная магнитная проницаемость m -го слоя b m * = 2 t m / R - относительная толщина m-го слоя t m - толщина m-го слоя q m - обобщенный параметр m-го слоя 1 - функция граничных условий д ля нижнего полубесконечного слоя , для воздуха ( = 1 , = 1 , = 0 ) 1 = 0 При анализе годографов для удобств а используют нормированные зависимости . Для НВТП нормировку производят по модулю максимального вносимого напряжения , которое соответствует идеально проводящему ОК и вычисляется при м = -1 : (3.23) Такая нормировка обобщает полученные результаты , расширяет область их примен е ния и делает их однозначными. Отметим , что для получения часто используемого в ВТК значения импеданса НВТП достаточн о разделить правую часть (3.18) на величину тока возбуждения I . 4. Обратная задача ВТК для НВТП Решение обратной задачи ВТК состоит в нахождении зависимости ( h ) распределения электропроводности по глубине пластины используя набор из N измеренных с помощью НВТП вносимых напряжений . Математически обратную задачу можно представить инт е гральным уравнением (4.1) Поскольку явного метода решения уравнения (4.1) н е существует , применим к нему метод квазирешения (п 5.3.2). В постановке для локального в смысле Чебышева критерия п о лучим корректную задачу минимизации функционала невязки : , i=1,N (4.2) Уч ет априорной информации в обратной задачи ВТК удобно проводить в виде инте р вала [ min , max ] , которому могут принадлежать значения электропроводности . В этом сл у чае можно рассматривать задачу (4.2) как задачу нелинейного программирования вида : (4.3) Заметим , что поскольку ограничения в задаче (4.3) являются линейными , разумным представляется применение метода условного градиент а (п 6.2.1). Рассмотрим процесс реш е ния системы (4.3) в предположении , что электропроводность аппроксимируется по узловым значениям j , j=1,M . (4.4) Линеариз уем функционал Ф в окрестности исследуемой точки 0 разложив его в ряд Тейлора с использованием только первых производных. (4.5) Пусть y = max Ф i ’ = Ф p ’ 0 . В этом случае мы можем свести задачу (4.4) к эквивален т ной задаче линейного программирования , состоящей в условной минимизации функции y . Рассмотрим процесс приведения задачи линейного программирования к каноническому в и ду . Рас крывая модуль в (4.5) получаем систему уравнений (4.6) Рассмотрим выражение под модулем в (4.5) и введем некоторые обозначения (4.7) (4.8) (4.9) (4.10) С учетом системы (4.8 - 4.10) постановка задачи (4.4) принимает вид (4.11) Раскрывая скобки в (4.11) и исключая y из первых 2N неравенств кроме р -го получаем систему неравенств (4.12) Приведем систему неравенств (4.12) к каноническому виду (6.1). Для этого , в с оотве т ствии со стандартным подходом , запишем все неравенства в виде равенств , добавляя в левые части неравенств неотрицательные переменные v . (4.13) В матричном виде полученная система имеет вид Ax = b (4.14) , где искомый вектор-столбец из 2( N+M)+1 элементов имеет вид x = ( y , z 1 , ... , z M , v 1 , ... , v 2N+M ) T . В системе линейных алгебраических уравнений (4.13) параметр м инимизации y определен строкой с номером p , которую в дальнейшем будем называть базовой. Рассмотрим алгоритм симплексного метода для решения задачи (4.14): 1. Выбор начального базиса - допустимого решения (4.14). В нашем случае базис до л жен состоять из 2 N+M переменных . Удобно задать начальный базис , присвоив дополнител ь ным переменным v i значения правых частей b i тех строк , в которых коэффициент матрицы A при них равен 1 . Начальное значение параметра минимизации y равно значению правой ч а сти базовой строк и . Все остальные компоненты искомого вектора х принимаются равными нулю. 2. Определение переменной , которая должна войти в очередной пробный базис . Для этого проводится анализ базовой строки p матрицы A . Из всех положительных элементов строки p , не явля ющихся коэффициентами при базисных переменных , выбирается элемент с наибольшим значением . Переменная , у которой этот элемент является коэффициентом , до л жен войти в очередной пробный базис , т.е . за новую базисную переменную принимается та , которая имеет н аибольший вес в функции y . Если в базовой строке p нет небазисных пер е менных с положительными коэффициентами , то в силу не отрицательности элементов х сл е дует сделать вывод , что оптимальному решению , т.е . минимуму y соответствует выбранный ранее базис . Выч исления завершаются также и при запрете изменения переменных по огр а ничениям. 3. Определение максимальной допустимой величины новой базисной переменной , не выходящей за пределы имеющихся ограничений . Вычисляются отношения значений правых частей (4.14) к соответствующим значениям коэффициентов при новой базисной переменной во всех строках , кроме базовой . При этом не рассматриваются отношения , в которых знам е натель равен нулю или отрицателен , т.е . при положительной правой части подобные случаи соответствую т бесконечным значениям переменных . Определяется номер строки q , где это отношение наименьшее . Новой базисной переменной присваивается значение отношения в строке q . Переменная , входившая в прежний базис и определявшаяся строкой q , исключается из базиса и приравнивается нулю . Если во всех строках , кроме базовой , коэффициенты при новой переменной равны нулю или отрицательны , то в силу не отрицательности элементов х и ограничения базиса ( 2N+M ) переменными , следует признать , что эта переменная не может на данн ом шаге вычислений войти в базис . В этом случае необходимо вернуться к пункту 2 , не рассматривая запрещенную переменную. 4. Преобразование системы (4.14) таким образом , чтобы в строке q коэффициент при вновь введенном параметре был равен 1 , а в остальны х строках - 0 . Это достигается путем линейных преобразований равенств , входящих в (4.14). Т.к . коэффициенты при параметрах , входящих в новое пробное базисное решение , становятся равными 1 и в каждую строку вх о дит только один базисный параметр , то значение нового базиса определяется правой частью уравнений . Далее следует возврат к пункту 2 . Решая систему (4.14) находим вектор min , соответствующий текущему решению зад а чи (4.13). Возвращаясь к методу условного градиента отметим , что направление спуска опр е деляется как -s n = min - 0 , а очередное итерационное решение задачи (4.3) определяется в ы ражением n+1 = n - s n . Для получения окончательного результата требуется определить оптимальную величину шага в направлении s n , что можно осуществить путем одноме р ной минимизации ф ункции ( )= n - s n методом золотого сечения. 5. Некорректные задачи 5.1 Основные понятия . Корректность по Адамару В самом общем виде большинство обратных задач может быть представлено в виде операторного уравнения A x = f , x X , f F ( 5.1 ) где А - оператор , оп ределенный на непустом множестве некоторого метрического простра н ства Х с метрикой X и действующий в метрическое пространство F с метрикой F , а по з а данному элементу f требуется определить решен ие х [ 10-14 ] . Введем в пространстве X норму || x ||= x i 2 и в пространстве F норму || f ||= f j 2 . Зам е тим , что метрики в соответствующих пространствах будут иметь вид (x,y) = \\ x-y \\ . В нашем случае обозначения в (5.1) имеют следующий смысл : А - оператор , согласно которому вычисляется величина относ и тельного напряжения , вносимого пластиной с электрической проводимостью ( h ) х ( h ) - электрическая проводимость пластины как функция глубины f U * вн - величина относительного вносимого напряжения НВТП Согласно классического определения задача (5.1) называется корректной по Адамару если при любой фиксиров анной правой части ее решение : · существует в Х · единственно в Х · непрерывно зависит от f В реальных условиях правая часть (5.1) известна всегда с некоторой погрешностью , т.е . f = f 0 + f , причем обычно f 0 принадлежи т пространству гладких функций , а погрешность f выводит ее из этого класса . Вследствие этого получаем постановку задачи , для решения к о торой невозможно применение обычных методов решения корректных задач , т.к . любой фиксирова нной правой части (5.1) соответствует бесконечное множество наборов исходных данных т.е . возможных распределений ЭП по глубине пластины. 5.2 Корректность по Тихонову Задача (5.1) называется корректной по Тихонову на множестве корректности М X е с ли : · точное решение задачи существует и принадлежит М · принадлежащее М решение единственно для любой правой части · принадлежащее М решение непрерывно относительно правой части В данном подходе к вопросу корректности существовани е решения и его принадле ж ность некоторому множеству не доказывается , а постулируется в самой постановке задачи. Физически гипотеза о принадлежности искомого решения определенному множеству корректности может интерпретироваться для нашей задачи предположен иями : 1. Исследуемая среда устроена не слишком сложно , т.е . ее физические характерист и ки ( , ) являются достаточно гладкими функциями ( т.е . их можно моделировать с п о мощью аппроксимаций типа кусоч но-постоянной , кусочно-линейной и т.п .). Предп о ложение основывается на физическом смысле поверхностной обработки . 2. Значения функций находятся во вполне определенных пределах ( для ( h ) истинность данного предположение не вы зывает сомнения ). 5.3 Вариационные методы решения некорректных задач Вариационные методы решения некорректных задач являются наиболее универсал ь ными из известных способов решения . Практически любая некорректная задача , для которой разработан какой-либо ме тод решения , может быть решена также и вариационным способом [15] . Для выбора подходящего метода решения обратной задачи рассмотрим постановки наиболее распространенных вариационных методов в терминах вычислительной математики и нашей задачи . Пусть фиксир ованный набор данных состоит из измеренных на N частотах N ко м плексных значений вносимой ЭДС U i изм , текущее рассчитанное значение которых U i ( ) . Требуется определить для выбранного типа аппроксимации ЭП значения М параметров а п проксимации ( обычно используются узловые значения ) . 5.3.1 Метод регуляризации Метод основан на стабилизации невязки ( Ax,f ) при помощи вспомогательного неотр и цательного функционала (x) . Идея метода состоит в том , чтобы минимизировать облада ю щий сглаживающими свойставами функционал Ф (x,f) , имеющий следующий вид : , параметр регуляризации 0 (5.2) Используя классический регуляризирующий функционал вида в терминах нашей задачи получаем : (5.3) Основное преимущество метода состоит в регуляризации простейшим способом , в ра м ках использован ия квадратичного функционала . Это позволяет использовать для решения некорректной задачи хорошо известные и легко программируемые методы минимизации квадратичных функционалов [ 17 ]. Оборотной стороной достоинств метода являются его недостатки. Требование м ин и мизации нормы решения и , как следствие , выбор гладкой реализации , в нашем случае будет приходить в противоречие с физикой задачи и в принципе не позволит находить решения с выраженным приповерхностным изменением ЭП. Еще один принципиальный недостаток ме тода состоит в постановке функционала как квадратичного , единого для всех измерений . Его минимум в общем случае не гарантирует минимизацию отклонения для произвольного i -го измерения в следствии нелокальности условия минимизации. Кроме того , следует учитыв ать отсутствие надежных априорных рекомендаций по в ы бору параметра регуляризации . Обычно подходящие значения можно подобрать только после ряда численных экспериментов по решению однотипных зада ч . Изменение характера искомого решения приводит к необходимости поиска нового значения . 5.3.2 Метод квазирешений Метод использует одну из форм критерия невязки и заключается в сведении невязки к минимуму на некотором непусто м множестве P , содержащем подмножество искомых реш е ний. Квазирешением уравнения (5.1) на множестве P X называется всякий элемент y P для которого справедливо равенство F ( Ay , f ) = inf( Ax , f ), x P . Понятие квазирешения обо б щает понятие решения , а для его существования не требуется принадлежность решения множеству P . Исходя из вышеизложенного получаем постановку метода в виде задач и условной м и нимизации функционала Ф ( x,f ) : (5. 4 ) Отметим , что множество Р может иметь простой вид , например интервала [ x min , x max ]. В терминах нашей задачи ВТК постановка задачи (5.4) пр имет вид : (5.5) Для того , чтобы гарантировать минимизацию отклонения для произвольного i -го изм е рения , можно применить к первому выражению в (5.5) локальный в смысле Чебышева кр и терий , в с оответствии с которым получаем окончательное выражение : (5.6) Основное преимущество метода состоит в том , что само понятие квазирешения снимает трудности с требованиями тихоновской коррект ности : первым (вызывающим переопред е ленность задачи ) и третьим (обычно принадлежность приближенной правой части уравнения (5.1) множеству N=AM неизвестна , а критерии этой принадлежности часто сами бывают н е устойчивы ). Кроме этого при рассмотрении задачи в виде (5.6) возможна постановка минимизац и онной задачи как задачи нелинейного программирования с явно заданными ограничениями на искомые переменные . В этом случае нет необходимости искажать исходный функционал регуляризующими членами как в (п 5.3.1), а треб ования к искомому решению можно удовл е творить , управляя ограничениями на параметры минимизации (в нашем случае - узловые значения ЭП ). 5.3.3 Метод невязки Рассмотрим множество Р формальных решений уравнения (5.1) Р = x : F ( Ax , f ) , где f - приближенная правая часть (5.1), известная с погрешностью . В качестве прибли женного решения (5.1) нельзя брать произвольный элемент множ е ства Р , т.к . не гарантируется близость Р к множеству точных решений . Для выбора прибл и женного решения предлагается использовать стабилизирующий функционал (х ) из (п 5.3.1) следующим образом : ( х ) = inf ( х ) , x P . Этот способ приводит к выбору элементов мн о жества Р имеющих минимально допустимую невязку. С учетом этого постановка метода состоит в условной минимизации функционала Ф (х ) : (5.7) Как и для метода регуляризации можно использовать стабилизирующий функционал вида (х )= ||x|| 2 , ч то приводит в обозначениях нашей задачи к системе : (5.8) При использовании локального в смысле Чебышева критерия система (5.8) окончател ь но примет вид : (5.9) 6. Нелинейное программирование Содержание нелинейного программирования составляют теория и методы решения з а дач о нахождении экстремумов нелинейных функций на множествах , определяемых лине й ными и нелинейными ограничени ями (равенствами и неравенствами ) [16-29] . Рассмотрим наиболее распространенные методы решения на примере основной задачи нелинейного программирования вида : (6.1) 6.1 Метод штрафных функций Идея метода состоит в замене экстремальной задачи с ограничениями (6.1) на задачу безусловной минимизации однопараметрической функции , > 0 (6.2) Непрерывную фу нкцию (х ) называют штрафом , если (х ) = 0 для х Х и (х ) > 0 в пр о тивном случае . Функция (х ) должна быть вы брана таким образом , чтобы решение задачи (6.2) сходилось при 0 к решению исходной задачи (6.1) или , по крайней мере , стремилось к нему. Приведем часто используемые выражения для штрафа : , k>0 (6.3) (6.4) (6.5) Наибольшее применение находит штраф (6.3). Выражение (6.5) гарантирует конечность метода при любом k>0 . При численной ре ализации метода штрафных функций возникают проблемы выбора начального значения параметра и способа его изменения . Сложность состоит в том , что выбор достаточно малого увеличивает вероятность сх одимости решения (6.2) к решению (6.1), а скорость сходимости градиентных методов вычисления точек минимума (6.2), как правило , падает с убыванием величины . 6.2 Релаксационные методы Релаксационным методом называют процесс п остроения последовательности точек х k : х k X , ( х k+1 ) ( х k ) ; k=0,1... . Основными представителями этого класса являю т ся мет оды спуска , алгоритм которых состоит из следующих шагов : 1. Выбор начального приближения х 0 2. Выбор в точке х k направления спуска - s k 3. Нахождение очередного приближения х k +1 = х k - k s k , где длина шага k >0 Различия методов состоят в выборе либо направления спуска , либо способа движения вдоль выбранного направления . В последнем случае обычно используют одномерную мин и мизацию функции х k +1 ( ) = х k - s k ( при этом точность вычисления точки минимума фун к ции х k +1 ( ) следует согласовывать с точностью вычисления значений функции (х ) ) или сп о соб удвоения (величина шага удваивается пока выполняется условие ( х k+1 ) ( х k ) ). 6.2.1 Метод условного градиента Идея метода заключается в линеаризации нелинейной функции (х ). В этом методе в ы бор направления спуска осуществляется следующим образом : Линеаризируя функцию (х ) в точке х К получаем Ф (х )= ( х k ) + ( ( х k ) ’ , х - х k ) Минимизируя линейную функцию Ф (х ) на множестве Х находим х min Направление спуска получаем как - s k = х min - х k Таким образом итерация метода имеет вид : x k+1 =x k + k (s k+1 - x k ) , s k+1 =arg min( f(x k ),x). Основное преимущество метода проявляется в случае задания допустимого множества с помощью линейных ограничений . В этом случае получаем задачу линейного программир о вания , решаемую стандартными методами (например симплексным ). 6.2.2 Метод проекции градиента Этот метод является аналогом метода градиентного спуска , используемого в задачах без ограничений . Его идея состоит в про ектировании точек , найденных методом наискоре й шего спуска , на допустимое множество , определяемое ограничениями . Проекцией точки y на множество Х называется точка P(y) Х такая , что || P(y) - y || || x - y || для всех х Х . Задача проектирования формализуется как || x - y || 2 min , x Х. Выбор направления спуска осуществляется следующим образом : Находим т очку r k = х k - ’ ( х k ) Находим проекцию p k точки r k на множество Х Направление спуска получаем как - s k = p k - х k Таким образом итерация метода имеет вид : x k+1 =P X [ x k - k f( x k ) ], где Р X (у ) - ортог о нальная проекция точки у на множество Х . Для отыскания направления спуска s k необходимо решить задачу минимизации квадр а тичной функции || r k - х || 2 на множестве Х . В общем случае эта задача того же порядка сло ж ности , что и исходная , однако для задач , допустимое множество которых имеет простую геометрическую структуру , отыскание проекции значительно упрщается . Например , дл я многомерного параллелепипида Q N = x R N : a x b , отыскание проекции осуществляе т ся путем сравнения n чисел и имеет вид P(x)= a i , x i a i ; x i , x i [ a i ,b i ] ; b i , x i b i . 6.2.3 Метод случайного спуска Метод характеризуется тем , что в качестве направления спуска s K выбирается некот о рая реализация n - мерной случайной величины S с известным законом распределения . Об э ф фективности этого метода судить трудно , однако благодаря использованию быстродейств у ющих ЭВМ он оказывается практически полезным. 6.3 Метод множителей Лагранжа Идея метода состоит в отыскании седловой точки функции Л агранжа задачи (6.1). Для нахождения решения вводится набор переменных i , называемых множители Лагранжа , и составляется функция Лагранжа , имеющая вид : (6.6) Алгоритм метода состоит в следующем : Составление функции Лагранжа Нахождение частных производных функции Лагранжа (6.7) 3. Решение системы из n+m уравнений вида (6.8) Решениями системы (6.8) являются точки , которые могут быть решениями задачи . 4. Выбор точек , в которых достигается экстремум и вычисление функции (х ) в этих точк ах. 7. Линейное программирование Задача линейного программирования в каноническом виде имеет вид [15,16]: (7.1) Приведение к каноническому виду любой задачи линейного программирования ос у щ ествляется путем введения дополнительных неотрицательных переменных , за счет чего ограничения , имеющие вид неравенств , принимают вид эквивалентных им равенств. Любая задача линейного программирования может быть решена за конечное число итераций с помощью симплексного метода [17,18] . Следует отметить , что поскольку этот м е тод разработан для неотрицательных элементов x j , это условие учитывается неявно и в с и стему уравнений (7.1) при численной реализации не входит. 7.1 Алгоритм симплексного метода 1. Приведен ие к каноническому виду 2. Выбор начального базиса 3. Проверка оптимальности базиса Матрицу А можно рассматривать как совокупность столбцов a j т.е . a j x j =b где j=1,N . Не ограничивая общности мож но считать , что базис образуют первые m столбцов , тогда остальные можно представить в виде a k = a j jk , j=1,m где jk . - некоторые числа . Рассмотрим коэффициенты k = c j jk - c k где j=1,m и k=1,N . Заметим , что для баз о вых столбцов k 0 . Проверка на оптимальность осуществляется следующим образом : k 0 , k=1,N - текущий базис оптимален - решение не ограничено сверху - существует другой , более подходящий базис 4. Составление нового базиса 4.1 Выбор элемента для введе ния в базис. В базис вводится любой столбец , для которого k 0 , обозначим его p 4.2 Выбор элемента исключаемого из базиса Из текущего базиса исключается столбец , для которого минимально отношение b i /A ip , i=1,M обозначим его b r /A rp 4.3 Преобразование вектора b и матрицы А по методу Жордана-Гаусса 4.4 Переход к пункту 3 8. Одномерная минимизация Несмотря на кажущуюся простоту , для широкого класса функций решение задачи м и нимиза ция функции одного переменного (х ) сопряжено с некоторыми трудностями . С о д ной стороны , в практических задачах часто неизвестно , является ли функция дифференцир у емой . С другой стороны , задача решения уравнения (х )=0 может на практике оказаться весьма сложной . Ввиду этого существенное значение приобретают методы минимизации , не требующие вычисления производной [15] . Поскольку нас интересует приближенное определение то чки минимума , то для этого исследуют поведение функции в конечном числе точек , способами выбора которых разл и чаются методы одномерной минимизации. К методам , в которых при ограничениях на количество вычислений значений (х ) д о стигается в определенном смысле наилучшая точность , относятся методы Фибоначчи и зол о того сечения [17,18] . Оба метода строятся по единой схеме и различаются способом выбора параметра . Исходными данными для них являются : отрез ок [a 0 ,b 0 ] содержащий точку м и нимума и точность определения точки минимума . 8.1 Алгоритм методов I. h 0 = b 0 - a 0 , k = 1 , (0.5,1) , h 1 = h 0 , h 2 = h 0 - h 1 , c 1 = a 0 + h 2 , d 2 = b 0 - h 2 II. Вычислить текущие значения ( c k ) и ( d k ) и действовать в соответствии с ними : ( c k ) ( d k ) ( c k ) ( d k ) a k = a k-1 c k-1 b k = d k-1 b k-1 d k = c k-1 c k = d k-1 h k+2 = h k -h k-1 h k -h k-1 d k = b k -h k+2 c k = a k +h k+2 III. Если ( h k ) то x min =min ( c k ) , ( d k ) иначе k++ и переход к шаг у II Следует отметить , что на каждом шаге кроме первого , производится только одно в ы числение значения функции ( x) . Легко показать , что для получения оптимальной последовательности отрезков , стяг и вающихся к точке минимума , нео бходимо положить k = F k-1 / F k , где F - число Фибоначчи. 8.2 Метод Фибоначчи Решая вопрос , при каких значениях параметра за конечное число итераций N мы п о лучим отрезок минимальной длины , получ им = N = F N-1 / F N . Иначе говоря , для поиска м и нимума первоначально необходимо найти число Фибоначчи N такое , что F N+1 (b 0 -a 0 )/ F N+2 . 8.3 Метод золотого сечения В реальной ситуации начиная поиск минимума мы не знаем точного числа требуемых итераций . Вместо вычисления будем выдерживать постоянное отношение длин интервало в h k-2 /h k-1 = h k-1 /h k = . При = ( 5+1)/2 = 1.618034 получаем метод золотого сечения. Сравнивая приведенные методы при больших значениях N можно показать , что знач е ние окончательного интервала неопределенности в методе золотого сечения лишь на 17% больше чем в методе Фибоначчи. 9. Результаты численного моделирования 9.1 Аппроксимации при численном моделировании Для построения моделей реальных распределений ЭП возмо жно применение целого ряда аппроксимаций . Все они могут быть разделены на два класса. 1. Аппроксимации , строящиеся по набору из произвольного числа узлов . Наиболее ра с пространенные из них : кусочно-постоянная , кусочно-линейная и сплайном . В условиях нашей з адачи указанные аппроксимации имеют несколько существенных недостатков : · Результаты аппроксимаций слабо согласуются с реальностью . Кусочно-постоянная и кусочно-линейная аппроксимации принципиально являются негладкими , а аппрокс и мация сплайном сглаживает все , в результате чего возникают значительные немон о тонности и всплески. · При увеличении количества узлов аппроксимации быстро нарастает неустойчивость процесса решения обратной задачи , для противодействия которой требуется примен е ние искусственных прием ов , не гарантирующих успеха . · В реальных условиях мы не имеем достоверной априорной информации о величине ЭП в узлах аппроксимации , расположенных в глубине пластины . 2. Аппроксимации , строящиеся по значениям ЭП на верхней и нижней поверхностях пл а стины и нескольким параметрам аппроксимации . Наиболее известные из них : экспоненц и альная , гиперболическим тангенсом и гауссоидой . Аппроксимации имеют вид : - аппроксимация экспоненциальная - аппроксимация гиперболическим тангенсом - аппроксимация гауссоидой где x - координата , равна нулю на нижней поверхност и пластины и единице на верхней 1 - величина электропроводности на верхней поверхности пластины 2 - величина электропроводности на нижней поверхности пластины - коэффициент , характеризующий крутизну экспоненты - коэффициент - коэффициент , характеризующий крутизну ; =0 соответствует случаю слоя с пр о водимостью 1 и толщиной на полупространстве с проводимостью 2 - коэффициент , характеризующий крутизну Для нашей задачи подобные аппроксимаци и являются предпочтительными , поскольку обладают заметными достоинствами : * Аппроксимации являются монотонными и гладкими , что хорошо согласуется с физ и ческой реальностью. * Пользуясь физически обоснованными рассуждениями мы можем получить необход и мую априорную информацию о величинах ЭП в приповерхностных слоях пластины . * Процесс решения обратной задачи существенно более устойчив и осуществляется зн а чительно быстрее Для иллюстрации наших рассуждений приведем пример применения приведенных в ы ше аппро ксимаций к случаю восстановления кусочно-линейной функции . По оси абсцисс отложена относительная глубина , по оси ординат электропроводность (МСм /м ). На графике показаны аппроксимации : кусочно п о стоянная ( SIci ),кусочно л и нейная ( SIli ), сплайн ( SIs ), экспоненциальная ( SIe ), г и перболическим тангенсом ( Sith ), гауссоидой ( SIg ). Легко заметить , что аппроксимация гипербол и ческим тангенсом хорошо описывает приповерхнос т ные изменения (аналогично экспоненциальной при большом показателе эксп о ненты ). Гауссоида может быть легко воспроизведена с помощью экспоненциальной аппроксимации , поэтому в дальнейшем использована не будет. 9.2 Модели реальных распределений электропроводности Модель задачи должна описывать некоторую пластину , подвергну тую поверхностной обработке . Для определенности зададим толщину пластины равной двум сантиметрам . На основе данных из Приложения 2 зададим значения ЭП вблизи нижней и верхней поверхн о стей соответственно 20 (МСм /м ) и 13 (МСм /м ) . Для решения обратной задачи необходимо задать априорную информацию о величине ЭП в узлах аппроксимации . В качестве таковой примем интервал [ 8,25 ] (МСм /м ) , получе н ный внесением 25% отклонения от считаемых истинными значений . Это отклонение мод е лирует неточность априорной информации . Из-за особенностей реализации алгоритма устойчивость решения сильно зависит от точности задания ЭП в узле , соответствующем нижней поверхности пластины , поэтому ограничение в нем зададим интервалом [ 19,21 ] (МСм /м ) . В нашем случае все возможные модели рас пределений ЭП могут быть разделены на два класса . Распределения относящиеся к первому из них условно назовем глубинными . В них ЭП претерпевает существенные изменения на протяжении всей глубины пластины . Второй класс образуют распределения , ЭП в которых за м етно изменяется лишь в приповерхностном слое глубиной порядка четверти пластины ., поэтому назовем эти распределения поверхнос т ными. Критерием отличия восстановленной функции распределения ЭП от модельной будем считать величину относительной погрешности , по скольку сравнение результатов с ее пом о щью вполне адекватно целям нашей работы. Следует отметит , что погрешность восстановления для поверхностных распределений ЭП представляет практический интерес в области , примерно ограниченной глубиной поря д ка четверти пластины , что обусловлено физическим смыслом поверхностной обработки . П о этому для случаев поверхностных распределений основное внимание будем уделят именно указанным глубинам. Для проверки возможности восстановления приповерхностных изменений ЭП рассмо т рим две базовые модели поверхностных распределений. Базовая модель A1 . Аппроксимация экспонентой . Проводимость 2 =20 [ МСм /м ] Проводимость 1 =13 [ МСм /м ] Показатель экспоненты = 25, 38, 120, 200 . Базовая модель A2 Аппроксимация гиперболическим тангенсом . Проводимость 2 =20 [ МСм /м ] Проводимость 1 = 6 [ МСм /м ] Коэффициент = 1 Коэффициент = 0.1, 0.05, 0.02, 0.01 Для проверки возможности восстановления глубинных распределений ЭП рассмотрим две базовые модели глубинных распределений. Базовая модель B1 Аппр оксимация кусочно-линейная. Проводимость задается в узлах с отсчитываемой от дна пл а стины относительной глубиной 0, 0.25, 0.5, 0.75, 1 . Узловые значения проводимости 20, 20, 17.6, 15.3, 13 , 20, 20, 20 , 16.5 , 13 , 20,20, 20 , 20 ,13 [ МСм /м ] . Базовая модель B2 Аппроксимация сплайном. Проводимость задается в узлах с отсчитываемой от дна пл а стины относительной глубиной 0, 0.25, 0.5, 0.75, 1 . Узловые значения проводимости 20, 20, 17.6, 15.3, 13 , 20, 20, 20 , 16.5 , 13 , 20,20, 20 , 20 ,13 [ МСм /м ] . Заметим , что на практике можно осуществить достаточно точное определение велич и ны ЭП приповерхностных слоев путем измерений проводимости традиционными средств а ми , поэтому дополнительно рассмотрим модельные зада чи при условии известной ЭП на верхней , а так же верхней и нижней поверхностях. Поскольку на практике результаты измерений вносимого напряжения имеют опред е ленную погрешность , все модели будем рассчитывать эмулируя погрешность U = 0 , 1 , 2 , 5% . Для исследования зависимости результатов восстановления распределений ЭП от ч а стоты возбуждения разобьем частотный диапазона три части следующим образом ( глубины проникновения приведены для случая постоянной ЭП =13 МСм /м ): Модели FA1L, FB1L Модели FA1M, FB1M Модели FA1H, FB1H f , [ Гц ] h , [m] f , [ КГц ] h , [m] f , [ КГц ] h , [m] 1 0.1396 5 0.001974 55 0.0005952 10 0.04414 10 0.001396 60 0.0005699 20 0.03121 15 0.00114 80 0.0004935 50 0.01974 20 0.000987 90 0.0004653 100 0.01396 25 0.0008828 10 0 0.0004414 200 0.00987 30 0.0008059 200 0.0003121 500 0.006243 35 0.0007461 300 0.0002549 1000 0.004414 40 0.0006979 500 0.0001974 2000 0.003121 45 0.000658 700 0.0001668 5000 0.001974 54 . 1 0.0006001 1000 0.0001396 Для исследования зависимости результатов восстановления распределений ЭП от числа измеряемых вносимых напряжений N рассмотрим случаи N= 5, 10, 15 . Низкие частоты f , [ Гц ] 1, 5, 10, 20, 35, 50, 100, 150, 200, 500, 750, 1000, 2000, 3500, 5000 f , [ Гц ] 1, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000 f , [ Гц ] 1, 20, 100, 500, 2000 Средние частоты f , [ КГц ] 5 , 7.5, 10, 15,17.5, 20, 25, 27.5, 30, 35, 37.5, 40, 45, 50, 54.1 f , [ КГц ] 5 , 10, 15, 20, 25, 30, 35, 40, 45, 54.1 f , [ КГц ] 5 , 15, 25, 35, 45 Высокие частоты f , [ КГц ] 55, 57.5, 60, 80, 85, 90,100, 150, 200, 300, 400, 500, 700, 850, 1000 f , [ КГц ] 55, 60, 80, 90,100, 200, 300, 500, 700, 1000 f , [ КГц ] 55, 80, 100, 300, 700 9.3 Принципиальная возможность восстановления Для ис следования возможности восстановления распределения ЭП рассмотрим резул ь таты , полученные в предположении наличия точных данных (погрешность измерения отсу т ствует ). На графиках в первых четырех пунктах Приложения 3 рассматриваемые зависим о сти показаны красн ым цветом (исходные данные черным ). Исходя из них можно сделать следующие выводы · Восстановление с помощью аппроксимации , использованной при эмуляции измерений (решении прямой задачи ), приводит к погрешности восстановления порядка 0.1% . · Восстановление глубинных распределений с помощью аппроксимаций экспоненциал ь ной и гиперболическим тангенсом возможно с хорошей точностью ( погрешность 2-5% ) для приповерхностных слоев глубиной порядка четверти пластины. · Восстановление глубинных распределений с помощ ью аппроксимаций сплайном и к у сочно-линейной возможно с хорошей точностью ( погрешность 2-3% ). Погрешность во с становления увеличивается с уменьшением глубины. · Восстановление поверхностных распределений с помощью аппроксимаций сплайном и кусочно-линейно й практически невозможно . Имеют место осцилляции , приводящие к погрешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксимаций экспоненц и альной и гиперболическим тангенсом возможно с хорошей точностью (погрешность 2-3%) . Погрешность восстановления увеличивается с уменьшением глубины , занимаемой распределением. 9.4 Восстановление по зашумленным данным Рассмотренные в данном разделе результаты демонстрируют возможность восстано в ления распределений ЭП в реальных условиях . Г рафики представлены в первых четырех пунктах Приложения 3 . На графиках рассматриваемые зависимости показаны цветами : результат восстановл е ния при погрешности данных равной 1% - голубым , результат восстановления при погре ш ности данных равной 2% - коричневы м , результат восстановления при погрешности данных равной 5% - синим. Исходя из них можно сделать следующие выводы : · Восстановление глубинных распределений с помощью аппроксимаций экспоненциал ь ной и гиперболическим тангенсом возможно с хорошей точностью ( погрешность 2-8% ) для приповерхностных слоев глубиной порядка четверти пластины. · Восстановление глубинных распределений с помощью аппроксимаций сплайном и к у сочно-линейной затруднено ( погрешность осциллирует в пределах 0-10% ). Погрешность восстановл ения увеличивается с уменьшением глубины , занимаемой распределением. · Восстановление поверхностных распределений с помощью аппроксимаций сплайном и кусочно-линейной практически невозможно . Имеют место осцилляции , приводящие к погрешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксимаций экспоненц и альной и гиперболическим тангенсом возможно с хорошей точностью (погрешность 3-6% для одноименных аппроксимаций и 7-10% в противном случае ) . Погрешность восст а новления увелич ивается с уменьшением глубины , занимаемой распределением. 9.5 Восстановление с учетом дополнительной информации С целью улучшения результатов восстановления в реальной обстановке , осложненной наличием зашумленных данных , следует использовать более жесткие ограничения на вел и чину ЭП в приповерхностных слоях . Для иллюстрации приведем несколько графиков , представленных в пятом пункте Пр и ложения 3 . На графиках рассматриваемые зависимости показаны цветами : базовые огран и чения - коричневым , жесткие ограничения н а верхней поверхности - голубым , жесткие огр а ничения на обоих поверхностях - малиновым. Исходя из полученных результатов можно сделать следующие выводы · Восстановление глубинных распределений с помощью аппроксимаций экспоненциал ь ной и гиперболическим тан генсом возможно с удовлетворительной точностью для пр и поверхностных слоев глубиной порядка четверти пластины . Дополнительные жесткие ограничения результаты восстановления не улучшают. · Восстановление глубинных распределений с помощью аппроксимаций сплайн ом и к у сочно-линейной затруднено . Дополнительные жесткие ограничения результаты восст а новления не улучшают. · Восстановление поверхностных распределений с помощью аппроксимаций сплайном и кусочно-линейной практически невозможно . Имеют место осцилляции , пр иводящие к погрешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксимаций экспоненц и альной и гиперболическим тангенсом возможно с удовлетворительной точностью (п о грешность 6 -10% ) . Погрешность восстановления уменьшаетс я при использовании д о полнительные ограничений примерно вдвое , особенно в приповерхностном слое толщ и ной порядка 10% . 9.6 Восстановление при различном возбуждении Для выбора необходимого количества измерений U вн * и соответствующих им частот возбуждения ВТП рассмотрим три возможных диапазона частот , в каждом из которых и с следуем случаи пяти , десяти и пятнадцати частот. На графиках в шестом пункте Приложения 3 рассматриваемые зависимости показаны цветами : 5 частот - коричневым , 10 частот - голубым , 15 частот - малиновым. Область низких частот Исходя из полученных результатов можно сделать следующие выводы · Восстановление глубинных распределений с помощью аппроксимаций экспоненциал ь ной и гиперболическим тангенсом возможно с удовлетворительной точностью для п р и поверхностных слоев глубиной порядка четверти пластины . Для улучшения результатов восстановления следует использовать 10 частот в случае погрешности данных не более 2% и 15 частот в противном случае. · Восстановление глубинных распределений с помощью ап проксимаций сплайном и к у сочно-линейной затруднено . Для улучшения результатов восстановления в приповер х ностном слоев глубиной порядка четверти пластины следует использовать 10 частот в случае погрешности данных не более 2% и 15 частот в противном случае. · Восстановление поверхностных распределений с помощью аппроксимаций сплайном и кусочно-линейной практически невозможно . Имеют место осцилляции , приводящие к погрешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксима ций экспоненц и альной и гиперболическим тангенсом возможно с удовлетворительной точностью (п о грешность 6 -8% ) . Для улучшения результатов восстановления следует использовать 10 частот в случае погрешности данных не более 2% и 15 частот в противном случае. Об ласть средних частот Исходя из полученных результатов можно сделать следующие выводы : · Восстановление глубинных распределений с помощью аппроксимаций экспоненциал ь ной и гиперболическим тангенсом возможно с удовлетворительной точностью для пр и поверхностн ых слоев глубиной порядка четверти пластины . Для улучшения результатов восстановления следует использовать 10 частот в случае погрешности данных не более 2% и 15 частот в противном случае. · Восстановление глубинных распределений с помощью аппроксимаций с плайном и к у сочно-линейной практически невозможно . Имеют место осцилляции , приводящие к п о грешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксимаций сплайном и кусочно-линейной практически невозможно . Имеют место осц илляции , приводящие к погрешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксимаций экспоненц и альной и гиперболическим тангенсом возможно с удовлетворительной точностью (п о грешность 6 -8% ) . Для улучшения результатов в осстановления следует использовать 10 частот в случае погрешности данных не более 2% и 15 частот в противном случае. Область высоких частот Исходя из полученных результатов можно сделать следующие выводы : · Восстановление глубинных распределений с помощью аппроксимаций экспоненциал ь ной и гиперболическим тангенсом возможно с удовлетворительной точностью для пр и поверхностных слоев глубиной порядка четверти пластины . Для улучшения результатов восстановления следует использовать 15, что позволяет восстанавлива тьраспределения с погрешностью (2-5)% . · Восстановление глубинных распределений с помощью аппроксимаций сплайном и к у сочно-линейной практически невозможно . Имеют место осцилляции , приводящие к п о грешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксимаций сплайном и кусочно-линейной практически невозможно . Имеют место осцилляции , приводящие к погрешностям , превышающим 10% . · Восстановление поверхностных распределений с помощью аппроксимаций экспоненц и альной и гипербол ическим тангенсом возможно с удовлетворительной точностью . Для улучшения результатов восстановления следует использовать 15 частот. 10 . Заключение По результатам проделанной работы можно сделать следующие выводы : * Существует принципиальная возможность восстановления как поверхностных так и гл у бинных распределений ЭП с погрешностью , не превышающей (2-3)%. Для восстановления поверхностных распределений следует использовать экспоненциальную и гиперболич е скую аппроксимации , а для глубинных сплайн и кусочно- постоянную (возможно использ о вание экспоненциальной и гиперболической аппрксимаций для в приповерхностном слое глубиной порядка четверти пластины ). * Существенное отрицательное влияние на результаты восстановления имеют погрешность измерения U вн * (не сле дует использовать данные с погрешностью измерения более 2% ) и малая глубина распределения ЭП (распределения ЭП сосредоточенные в приповерхнос т ном слое глубиной менее (3-5)% восстанавливаются хуже ). * Использование жестких ограничений на величину ЭП в пр иповерхностных слоях опра в дано при восстановлении поверхностных распределений , причем при наличии данных с погрешностью , превосходящей 2% , или малой глубины распределения предпочтительнее задавать ограничения на обеих поверхностях . При зашумленности данных порядка (1-2)% достаточно задать жесткие ограничения лишь на верхней поверхности. * В наборе частот возбуждения ВТП должны присутствовать низкочастотные составля ю щие , влияние которых особенно заметно при работе с глубинными распределениями и с о ответству ющими аппроксимациями . Рекомендуется использовать порядка десяти частот , равномерно распределенных по частотному диапазону (0.001 70)КГц . В условиях высокой погрешности измерений или отчетливо выраженных приповерхностных измен ениях ЭП заметное положительное влияние оказывает увеличение числа частот возбуждения ВТП (например , до пятнадцати .). В процессе работы над задачей был проведен анализ литературы , выбрана модель зад а чи и способы ее аппроксимации . При помощи программы , разр аботанной согласно предл о женной модели , были проведены расчеты модельных задач и рассмотрены результаты во с становления распределений ЭП в зависимости от основных влияющих факторов . Таким образом , цели , поставленные в техническом задании , решены в полном о бъеме. 11. Литература 1. Неразрушающий контроль качества изделий электромагнитными методами, Герасимов ВГ , 1978,215 2. Вихретоковый контроль накладными преобразователями., Герасимов ВГ ,1985,86 3. Вихретоковые методы и приборы неразрушающего контроля., Рудаков ВН, 1992, 72 4. Накладные и экранные датчики., Соболев ВС, 1967, 144 5. Теория и расчет накладных вихретоковых преобразователей., Дякин ВВ, 1981, 135 6. Основы анализа физических полей.,Покровский АД, 1982, 89 7. Дефектоскопия металлов., Денель АК, 1972, 303 8. Индукционная структуроскопия., Дорофеев АЛ ,1973,177 9. Структура и свойства металлов и сплавов.Справочник., Шматко ОА ,1987,580 10. Некорректные задачи Численные методы и приложения., Гончарский АВ ,1989,198 11. Некорректные задачи матф изики и анализа., Лаврентьев ММ ,1980,286 12. Линейные операторы и некорректные задачи., Лаврентьев ММ ,1991,331 13. Методы решения некорректно поставленных задач Алгоритмич . аспект., Морозов ВА , 1992,320 14. Численные методы решения некорректных задач., Тихонов АН ,1990,230 15. Начала теории вычислительных методов, Крылов ВИ ,1984,260 16. Математическое программирование в примерах и задачах., Акулич ИЛ ,1993,319 17. Математическое программирование., Карманов ВГ ,1986,286 18. Математическое программировани е., Орехова РА ,1992,290 19. Нелинейное программирование Теория и алгоритмы., Базара М ,1982,583 20. Прикладное нелинейное программирование., Химмельблау Д ,1975,534 21. Введение в методы оптимизации., Аоки М ,1977,344 22. Введение в оптимизацию., Поляк БТ ,1983,384 23. Курс методов оптимизации., Сухарев АГ ,1986,326 24. Практическая оптимизация., Гилл Ф ,1985,509 25. Численные методы оптимизации., Полак Э ,1974,367 26. Алгоритмы решения экстремальных задач., Романовский ИВ ,1977,352 27. Методы решения экст ремальных задач., Васильев ФП ,1981,400 28. Методы решения экстремальных задач и их применение в системах оптимизации., Евт у шенко ЮГ , 1982,432 29. Численные методы решения экстремальных задач., Васильев ФП ,1988,549 30. Введение в вычислительную физику., Федоренко РП ,1994,526 31. Методы математической физики., Арсенин ВЯ ,1984,283 32. Уравнения математической физики., Тихонов АН ,1977 33. Уравнения математической физики., Владимиров ВС ,1988,512 34. Метод интегральных уравнений в теории рассеивания., Колт он Д ,1987,311 35. Теория электромагнитного поля., Поливанов КМ ,1975,207 36. Eddy current testing. Manual on eddy current method., Cecco VS,1981,195 37. Optimization methods with applications for PC., Mistree F,1987,168 38. Electromagnetic inverse p rofiling., Tijhuis AG,1987,465 39. Inverse acoustic and electromagnetic scattering theory., Colton D,1992,305 40. Накладной электромагнитный преобразователь над объектом контроля с изменяющим и ся по глубине электрич ескими и магнитными свойствами , Касимов ГА , Кулаев ЮВ , Д е фектоскопия , 1978, № 6, с 81-84 41. Возможн ости применения методов теории синтеза излучающих систем в задачах эле к тромагнитного контроля , Кулаев ЮВ , 1980, тематический сборник Труды МЭИ , в ы пуск 453, с 12-18 42. Analitical solutions to eddy-current probe-coil problems , Deeds WE, Dodd CV, Journal of Applied Phisics , 19 6 8, vol 39 , 3, p 2 8 29 - 2838 43. General analysis of probe coils near stratified conductors , Deeds WE, D odd CV, International Journal of Nondestructive Testing , 1971, vol3, 2, p109-130 44. Tutorial. A review of least-squares inversion and its application to geophysical problems , Lines LR, Treitel S, Geophysical Prospecting , 1984, vol32, 2, p159-186 45. Eddy current calculations using half-space Green ’ s functions , Bowler JR, J ournal of A p plied Phisics , 1987, vol61, 3, p833-839 46. Reconstruction of 3D conductivity variations from eddy current( electromagnetic induction ) data , Nair SM, Rose JH, Inverse Problems , 1990, 6, p1007-1030 47. Electromagnetic induction (eddy-currents) in a conducting half-space in the absence and pre s ence of inhomogeneities: a new formalism , Nair SM, Rose JH, Journal of Applied Phisics , 1990, vol68, 12, p5995-6009 48. Eddy-current p robe impedance due to a volumetric flaw , Bowler JR, Journal of Applied Phisics , 1991, vol70, 3, p1107-1114 49. Theory of eddy current inversion , Bowler JR, Norton SJ, Journal of Applied Phisics , 1993, vol73, 2, p501-512 50. Impedance of coils over layered metals with continuously variable conductivity and permeabi l ity: Theory and experiment , Rose JH, Journal of Applied Phisics , 1993, vol74, 3, p2076 51. Eddy-current interaction with ideal crack , Bowler JR, Journal of Applied Phisics , 1994, vol75, 12, p8128,8138 52. Method of so lution of forward problems in eddy-current testing , Kolyshkin AA, Journal of Applied Phisics , 1995, vol77, 10, p4903-4912 Приложение 1. Программная реализация Программная реализация изложенного метода решения обратной задачи ВТК ос у ществлена при помощи компилятора Borland Pascal 7.0 и состоит из шести модулей : 1. ErIn12.pas - исполняемый файл , осуществляет ос новной цикл программы 2. EData.pas - содержит глобальные данные и осуществляет чтение файла исходных да н ных 3. EFile.pas - содержит вспомогательные функции и иосуществляет сохранение резул ь татов расчетов 4. EMath.pas - осуществляет поддержку операци й с комплексными числами 5. EDirect.pas - осуществляет решение прямой задачи ВТК 6. EMinimum.pas - осуществляет решение обратной задачи ВТК П 1.1 Исходные данные Исходные данные программы хранятся в текстовом файле ( кодировка ASCII , расшир е ние по умолча нию TXT ). HThick - толщина пластины ,[мм ] nPoints - количество узлов аппроксимации электропроводности для PWL,PWC,SPL аппроксимаций . В случае EXP,HTG аппроксимации вычисление значений ЭП в них производится по окончании расчетов nLayers - количество инте рвалов с кусочно-постоянной электропроводностью , на которые разбивается пластина для непосредственного расчета вносимой ЭДС по реккурентным формулам для многослойной пластины nFreqs - количество частот возбуждения гармоник вносимой относительной ЭДС nSta b - число стабилизируемых значащих цифр epsU - погрешность измерения aG - коэффициент сжатия ограничений (aG<=1); НЕ используется при EXP,HTG аппроксимации nApprox - типы аппроксимации прямой и обратной задач si - значения проводимости в узлах аппрокси мации siMin, siMax - ограничения на возможные значения проводимости в узлах аппроксим а ции в процессе решения ОЗ incVal - величина dx для численного дифференцирования maxSteps - максимальное число отрезков интегрирования maxX - верхний предел интегриров ания при расчете Uvn* Eps - погрешность интегрирования при расчете Uvn* dType - тип разностной производной (=1 правая или =2 центральная ) eqlB - толщины слоев пластины одинаковы ( b=hThick/nLayers) если eqlB >0, в противном случае используются координаты слоев из файла П 1.2 Используемые аппроксимации Примечание . Координата Х [0,1] отсчитывается от дна пластины для всех аппроксимаций. Сплайн ( SPL ), кусочно-линейная ( PWL ), кусочно-постоянная ( PWC ) аппроксимации. В процессе расчето в ищутся значения электропроводности в узлах аппроксимации , причем количество узлов увеличивается от едениицы до nPoints в целях сохранения усто й чивости решения. Начальные значения (узловые значения истинной ЭП для эмуляции измерений U * вн ) задаются в столбце si файла исходных данных , начальные значения ограничений на узл о вые значения ЭП в столбцах siMin и siMax (движение по столбцу сверху вниз соотве т ствует изменению координаты от дна пластины до обрабатываемой повехности ) . Экспоненциальная аппроксимация ( EXP) В случ ае задания экспоненциальной аппроксимации зависимость электропрводности от толщины представляется в виде SIGMA = ( siE-siI )*EXP( -alfa*(1-x) ) + siI Варьруемыми параметрами являются эектропроводность на верхней поверхности siЕ , электропроводность "на бес конечности " siI и параметр alfa . В файле исходных данных в таблице из nPoints строк с подзаголовком "si siMin siMax", информация об ограничениях на параметры siE, siI задается в первой и nPoints -строке . Величина и ограничения для пар а метра alfa задаю тся первой строкой в "special approximation parameters". Аппроксимация гиперболическим тангенсом (HTG) В случае задания аппроксимации гиперболическим тангенсом зависимость электр о прводности от толщины представляется в виде SIGMA = si2 + ( si1-si2 )/2* 1 + th( ( beta-x )/gamma ) Величина и ограничения для параметров si2,beta,gamma задаются начиная со второй строки в "special approximation parameters", для si1 аналогично siI. П 1.3 Результаты расчета Результаты расчета помещаются в текстовый файл ( кодир овка ASCII , расширение по умолчанию LST ), при этом результат каждой итерации отбражается строкой вида : 1 <Ф > 0.000353 Rg= 17.003 15.639 9.697 где первая цифра (в данном случае 1) соответствует номеру текущей внутренней итерации , затем посл е текста "<Ф >" идет значение текущей абсолютной среднеквадратичной невязки по всем гармоникам (в данном случае 0.000353), затем после текста "Rg= ", идут искомые текущие значения переменных минимизации . В случае SPL,PWL,PWC аппроксимаций это непосредствен но узловые значения электропроводности для текущего количества узлов , а для EXP,HTG аппроксимаций это параметры siE, siI, Alfa или si1, si2, Beta, Gamma . B качестве последней строки помещаются n Points вычисленных значений э /проводности в ра в номерно расположенных узлах пластины. П 1 . 4 Основная программа ErIn (****************************************************************************) (* ErIn v1.42 *) ( * Eddy current inverse problem solver. *) (* (C) 1999 by Nikita U.Dolgov *) (* Moscow Power Engineering Institute , Introscopy dept. *) ******************* ********************************************************* Program ErIn; 23.02.99 Uses DOS,CRT, EData, EMath, EDirect, EFile, EMinimum; Var m, mLast, i : byte; loop counters procedure about; Let me to introduce myself begin clrscr; GetTime( clk1.H, clk1.M, clk1.S, clk1.S100 ); get start time writeln('***********************************************************'); writeln('* ErIn v1.42 Basic * *'); writeln('***********************************************************'); end; procedure initParameters; var apDT : byte; approximation type for direct task begin apDT := n Approx SHR 4; XXXXYYYY->0000XXXX fHypTg:=(( apDT AND apHypTg ) = apHypTg); if fHypTg then begin si0[ 1 ]:=si[ 1 ]; si1 - conductivity about bottom of slab si0[ 2 ]:=par0[ 2 ]; si2 - conductivity about top of slab si0[ 3 ]:=par0[ 3 ]; Beta - ratio of approx. si0[ 4 ]:=par0[ 4 ]; Gamma- ratio of approx. mCur:=4; end else if(( apDT AND apExp ) = 0 ) then It 's not an EXP approx. begin for i:=1 to nPoints do si0[ i ] :=si [ i ]; SI data from file mCur:=nPoints; end else begin si0[ 1 ]:=si[ 1 ]; siI - conductivity about bottom of slab s i0[ 2 ]:=si[ nPoints ]; siE - conductivity about top of slab si0[ 3 ]:=par0[ 1 ]; Alfa- ratio of approx. mCur:=3; end; setApproximationType( apDT ); approx. type for direct problem setApproximationData ( si0, mCur ); approx. data for direct problem nApprox := ( nApprox AND $0F ); XXXXYYYY->0000YYYY fHypTg := (( nApprox AND apHypTg ) = apHypTg ); fMulti := (( nApprox AND apExp ) = 0 ) AND NOT fHypTg; It's no t an EXP approx. if fMulti then begin for i:=1 to nPoints do begin Gr[ 1,i ]:=SiMax[ i ]; Gr[ 2,i ]:=SiMin[ i ]; Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2; zero estimate of SI Rgs[ i ]:=1E33; biggest integer end; mLast:=nPoints; loop for every node of approx. mCur :=1; to begin from the only node of approx end else if fHypTg then begin Gr[ 1,1 ]:= siMax[ 1 ]; Gr[ 2,1 ]:= siMin[ 1 ]; Rgs[ 1 ]:=1E33; Gr[ 1,2 ]:=parMax[ 2 ]; Gr[ 2,2 ]:=parMin[ 2 ]; Rgs[ 2 ]:=1E33; Gr[ 1,3 ]:=parMax[ 3 ]; Gr[ 2,3 ]:=parMin[ 3 ]; Rgs[ 3 ] :=1E33; Gr[ 1,4 ]:=parMax[ 4 ]; Gr[ 2,4 ]:=parMin[ 4 ]; Rgs[ 4 ]:=1E33; for i:=1 to 4 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2; mLast:=1; mCur:=4; end else begin Gr[ 1,1 ]:= siMax[1]; Gr [2,1]:= siMin[1]; Rgs[ 1 ]:=1E33; Gr[ 1,2 ]:= siMax[nPoints]; Gr[2,2]:= siMin[nPoints]; Rgs[ 2 ]:=1E33; Gr[ 1,3 ]:= parMax[1]; Gr[2,3]:= parMin[1]; Rgs[ 3 ]:=1E33; for i:=1 to 3 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2, i ] )/2; mLast:=1; mCur :=3; end; initConst( nLayers, parMaxH, parMaxX , parEps, parEqlB ); set probe params end; procedure directTask; emulate voltage measurements [with error] begin for i:=1 to nFreqs do begin getVoltage( freqs[i], Umr[ i ], Umi[ i ] ); "measured" Uvn* if ( epsU > 0 ) then add measurement error begin randomize; Umr[ i ]:=Umr[ i ] *( 1 + epsU*( random-0.5 ) ); randomize; Umi[ i ]:=Umi[ i ]*( 1 + epsU*( random-0.5 ) ); end; end; writeln('* Voltage measurements have been emulated'); setApproximationType( nApprox ); approx. type for inver se problem setApproximationData( Rg, mCur ); approx. data for inverse problem end; procedure reduceSILimits; evaluate SI for m+1 points of approx. using aG var x0, x1, xL, dx, Gr1, Gr2 : real; j, k : byte; begin ------------- ---------------- get SI min/max for m+1 points of approximation dx:=1/( nPoints-1 ); for i:=1 to m+1 do begin k:=1; x1:=0; x0:=( i-1 )/m; for j:=1 to nPoints-1 do begin xL:=( j-1 )/( nPoints-1 ); if( ( xL < x0 ) AND ( x0 <= xL+dx ) )then begin k:=j; x1:=xL; end; end; Gr[ 1,i ]:=siMax[ k ] + ( siMax[ k+1 ]-siMax[ k ] )*( x0-x1 )/dx; Gr[ 2,i ]:=siMin[ k ] + ( siMin[ k+1 ]-siMin[ k ] )*( x0-x1 )/dx; end; ------------------------------------- get SI for m+1 points of approximation for i :=1 to m+1 do begin Rg[i]:=getSiFunction( (i-1)/m ); if ( Rg[i] > Gr[1,i] )then Rg[i]:=Gr[1,i]; if ( Rg[i] < Gr[2,i] )then Rg[i]:=Gr[2,i]; if m > 1 then There're more than 1 point of approx. begin Gr1:= Rg[i]+( Gr[1,i]-Rg[i] )*aG; reduce upper bound Gr2:= Rg[i]-( Rg[i]-Gr[2,i] )*aG; reduce lower bound if ( Gr1 < Gr[1,i] )then Gr[1,i]:=Gr1; test o verflow if ( Gr2 > Gr[2,i] )then Gr[2,i]:=Gr2; end; end; setApproximationData( Rg , m+1 ); end; procedure resultMessage; to announce new results begin if fMulti then begin writeln(' current nodal values of conductivity'); write(' si : ');for i:=1 to m do write(Rg[i] :6:3,' ');writeln; write(' max: ');for i:=1 to m do write(Gr[1,i]:6:3,' ');writeln; write(' min: ');for i:=1 to m do write(Gr[2, i]:6:3,' ');writeln; end else begin for i:=1 to nPoints do si[i]:=getSiFunction( ( i-1 )/( nPoints-1 ) ); if fHypTg then saveHypTgResults else saveExpResults; end; end; procedure clockMessage; user-fri e ndly message begin writeln('***********************************************************'); write( '* approximation points number :',m:3,' * Time '); clock; writeln('*************************** ********************************'); end; procedure done; final message begin Sound(222); Delay(111); Sound(444); Delay(111); NoSound; beep write('* Task processing time '); clock; saveTime; writeln('* Status: Inverse problem has been successfully evaluated.'); end; Begin about; loadData; initParameters; directTask; for m:=1 to mLast do begin if fMulti then begin mCur:=m; clockMessage; end; doMinimization; main part of work setApproximationData( Rg, mCur ); set new approx. data resultMessage; if(( fMulti )AND( m < nPoints ))then reduceSILimits; end; done; End. П 1 . 5 Модуль глобальных данных EData Unit EData; Interface Uses DOS; Const maxPAR = 40; nodes of approximation max number maxFUN = 20; excitation frequencies max number maxSPC = 4; support approximation values number iterImax = 50; max number of internal iterations Const apSpline = 1; approximation type identifiers apHypTg = 3; apEx p = 2; apPWCon = 4; apPWLin = 8; Type Parameters = array[ 1..maxPAR ] of real; Si,Mu data Functionals = array[ 1..maxFUN ] of real; Voltage data SpecialPar = array[ 1..maxSPC ] of real; Special data Var hThick : real; thickness of slab nPoints : integer; nodes of approximation number within [ 0,H ] nLayers : integer; number of piecewise constant SI layers within[ 0,H ] nFreqs : integer; number of excitation frequencies nStab : integer; required number of true digits in SI estimation epsU : real; relative error of measurement Uvn* nApprox : byte; approximation type identifier incVal : real; step for numerical differ. parMaxH : intege r; max number of integration steps parMaxX : real; upper bound of integration parEps : real; error of integration derivType: byte; 1 for right; 2 for central Var freqs : Functionals; frequencies of excitment Uvn* Umr, Umi : Functionals; Re(Uvn*),Im(Uvn*) for "measured" data Uer, Uei : Functionals; Re(Uvn*),Im(Uvn*) for estimated data mu : Parameters; relative permeability nodal values si, si0 : Parameters; conductivity approximation nodal values siMin, siMax : Parameters; conductivity nodal values restrictions par0 : SpecialPar; alfa,si2,beta,gamma - for exp&HypTg parMin,parMax: SpecialPar; -||- min/m ax zLayer : Parameters; relative borders of slab layers [0,1] Var aG : real; scale factor for SImin/max Ft : real; current discrepancy functional value fMulti : boolean; TRUE if it isn't an EXP-approximation fHypTg : boolean; TRUE for Hyperbolic tg approximation parEqlB : boolean; TRUE if b[i]=const mCur : integer; current number of approximation nodes inFileName : strin g; data file name outFileName : string; results file name Var Rg : Parameters; current SI estimation RgS : Parameters; previous SI estimation Gr : array [ 1. .2 ,1..maxPAR ] of real; SI max/min Fh : array [ 1..maxPAR , 1..maxFUN ] of real; current discrepancies Type TTime=record H, M, S, S100 : word; hour,min,sec,sec/100 end; Var clk1, clk2 : TTime; start&finish time procedure loadData; load all user-defined data from file Implementation procedure loadData; var i,eqlB : integer; FF : text; begin assign( F F, outFileName ); clear output file rewrite( FF ); close( FF ); assign( FF, inFileName ); read input file reset( FF ); readln( FF ); readln( FF ); readln( FF, h Thick, nPoints, nLayers, nFreqs, nStab, epsU, aG, nApprox ); readln( FF ); for i:=1 to nFreqs do read( FF, freqs[i] ); readln( FF ); readln( FF ); readln( FF ); for i:=1 to nPoints do readln( FF, si[i], siMin[i], siMax[i] ); readln( FF ); readln( FF ); readln( FF , incVal, parMaxH, parMaxX, parEps, derivType, eqlB ); readln( FF ); readln( FF ); for i:=1 to maxSPC do readln( FF, par0[i] , parMin[i] , parMax[i] ); readln( FF ); i f ( eqlB=0 )then begin for i:=1 to nLayers+1 do read( FF, zLayer[i] ); parEqlB:=false; end else parEqlB:=true; close( FF ); for i:=1 to maxPAR do mu[i]:=1; end; Var str : string; Begin if( ParamCount = 1 )then str:=ParamStr(1) else begin write('Enter I/O file name, please: '); readln( str ); end; inFileName :=str+'.txt'; outFileName:=str+'.lst'; End. П 1 . 6 Модуль работы с файлами EFile Unit EFile; Interface Uses DOS, EData; function isStable( ns : integer; var RG1,RG2 ) : boolean; function saveResults( ns,iter : integer ) : boolean; procedure saveExpResults; procedure saveHypTgResults; procedure clock; procedure saveTime; Implementation Var FF : t ext; i : byte; function decimalDegree( n:integer ) : real; 10^n var s:real; i:byte; begin s:=1; for i:=1 to n do s:=s*10; decimalDegree:=s; end; function isStable( ns:integer ; var RG1,RG2 ) : boolean; var m : real; R1 : Parameters absolute RG1; R2 : Parameters absolute RG2; begin isStable:=TRUE; m:=decimalDegree( ns-1 ); for i:=1 to mCur do begin if NOT(( ABS( R2[i]-R1[i] )*m )<=ABS( R2[i]) ) then isStable:=FALSE; RgS [i]:=Rg[i]; end; end; function saveResults( ns , iter : integer ) : boolean; var sum : real; begin sum:=0; for i:=1 to nFreqs do sum:=sum + Fh[1,i]; sum:=SQRT( sum/nFreqs ); assign( FF , outFileName ); append( FF ); write( iter:2, ' <” >', sum:10:7, ' Rg=' ); write( FF , iter:2, ' <” >', sum:10:7, ' Rg='); for i:=1 to mCur do begin write( Rg[i]:6:3, ' '); write( FF , Rg[i]:6:3, ' '); end; writeln; writeln( FF ); close( FF ); saveResults:=isStable( ns , Rgs , Rg ); end; procedure saveExpResults; begin assign( FF , outFileName ); append( FF ); writeln( ' siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3); wr iteln( FF , ' siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3); write( ' SI: '); write( FF , ' SI: '); for i:=1 to nPoints do begin write( si[i]:6:3,' '); write( FF , si[i]:6:3,' '); end; wri teln; writeln( FF ); close( FF ); end; procedure saveHypTgResults; begin assign( FF , outFileName ); append( FF ); writeln( ' si1=',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,' gamma=',Rg[4]:6:3); writeln( FF , ' si1 =',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,' gamma=',Rg[4]:6:3); write( ' SI: '); write( FF , ' SI: '); for i:=1 to nPoints do begin write( si[i]:6:3,' '); write( FF , si[i]:6:3,' '); end; writeln; writeln( FF ); close( FF ); end; procedure clock; t2 = t2-t1 var H1,M1,S1,H2,M2,S2,sec1,sec2 : longint; begin GetTime( clk2.H, clk2.M, clk2.S, clk2.S100 ); current time H2:=clk2.H; M2:=clk2.M; S2:=clk2.S; H1:=clk1.H; M1:=clk1.M; S1:=clk1.S; sec2:= ( H2*60 + M2 )*60 + S2; sec1:= ( H1*60 + M1 )*60 + S1; if( sec2 < sec1 )then sec2:=sec2 + 85020; +23.59.59 sec2:=sec 2 - sec1; clk2.H := sec2 div 3600; sec2:=sec2 - clk2.H*3600; clk2.M := sec2 div 60; sec2:=sec2 - clk2.M*60; clk2.S := sec2; writeln( clk2.H:2, ':', clk2.M:2, ':', clk2.S:2 ); end; procedure saveTime; begin assign( FF , outFileName ); append( FF ); write( FF ,'* Processing time ',clk2.H:2, ':', clk2.M:2, ':', clk2.S:2 ); close( FF ); end; End. П 1 . 7 Модуль решения прямой задачи ВТК для НВТП EDirect ***************************************** *********************************** ERIN submodule : EDirect , 15.02.99, (C) 1999 by Nikita U.Dolgov **************************************************************************** Estimates Uvn* for Eddy current testing of inhomogeneous m ultilayer slab with surface( flat ) probe. It can do it using one of five types of conductivity approximation : Spline, Exponential, Piecewise constant, Piecewise linear,Hyperbolic tangent **************************************************************************** $F+ Unit EDirect; Interface Uses EData, EMath; Type siFunc = function( x:real ) : real; Var getSiFunction : siFunc; for external getting SI estimate procedure in itConst( par1,par2:integer; par3,par4:real; par5:boolean ); procedure getVoltage( freq : real ; var ur,ui : real ); Uvn* = ur + j*ui procedure setApproximationType( approx : byte ); type of approx. procedure setApproximationItem( SIG:real ; N : byte ); set SIGMA[ N ] procedure setApproximationData( var SIG; nVal : byte ); SIGMA[1..nVal] procedure getApproximationData( var SIG ; var N : byte ); get SIGMA[ N ] Implementation Const PI23 = 2000*pi; 2*pi*KHz mu0 = 4*pi*1E-7; magnetic const Var appSigma : Parameters; conductivity approximation data buffer appCount : byte; size of conductivity approximation data buffer appType : byte; conductivity approximation type identifier Type commonInfo=record w : real; cyclical excitation frecuency R : real; equivalent radius of probe H : real; generalized lift-off of probe Kr : real; parameter of probe eps : real; error of integration xMax : real; upper bound of integration steps : integer; current number of integration steps maxsteps: integer; max number of integration steps Nlay : integer; number of layers in slab sigma : Parameters; conductivity of layers m : Parameters; relative permeability of layers b : Parameters; thickness of layers zCentre : Parameters; centre of layer end; procFunc = procedure( x:real; var result:complex); Var siB, siC, siD : Parameters; support for Spline approx. cInfo : commonInfo; one-way access low level info function siSpline( x:real ) : real; Spline approxi mation begin if( appCount = 1 )then siSpline := appSigma[ 1 ] else siSpline:=Seval( appCount, x, appSigma, siB, siC, siD); end; function siExp( x:real ) : real; Exponential approximation begin siExp:=(appSigma[2]-appSigm a[1])*EXP( -appSigma[3]*(1-x) ) + appSigma[1]; end; function siPWConst( x:real ) : real; Piecewise constant approximation var dx, dh : real; i : byte; begin if( appCount = 1 )then siPWConst := appSigma[ 1 ] else begin dh:=1/ ( appCount-1 ); dx:=dh/2; i:=1; while( x > dx ) do begin i:=i + 1; dx:=dx + dh; end; siPWConst:=appSigma[ i ]; end; end; function siPWLinear( x:real ) : real; Piecewise linear approximation var dx, dh : real; i : byte; begin if( appCount = 1 )then siPWLinear := appSigma[ 1 ] else begin dh:=1/( appCount-1 ); dx:=0; i:=1; repeat i:=i + 1; dx:=dx + dh; until( x <= dx ); siPWLinear:=appSigma[i-1]+( appSigma[i]-appSigma[i-1] )*( x/dh+2-i); end; end; function siHyperTg( x:real ) : real; Hyperbolic tangent approximation begin siHyperTg:=appSigma[2]+(appSigma[1]-appSigma[2])*(1+th((appSigma[3]-x)/appSigma[4]))/2; end; procedure setApproximationType( approx : byte ); begin appType := approx; write('* conductivity approximation type : '); case approx of apSpline : begin writeln('SPLINE'); getSiFunction := siSpline; end; apExp : begin wri teln('EXP'); getSiFunction := siExp; end; apPWCon : begin writeln('PIECEWISE CONST'); getSiFunction : = siPWConst; end; apPWLin : begin writeln('PIECEWISE LINEAR'); getSiFunction := siPWLinear; end; apHypTg : begin writeln('HYPERBOLIC TANGENT'); getSiFunction := siHyperTg; end; end; end; procedure setApproximationData( var SIG ; nVal : byte ); var Sigma : Parameters absolute SIG; i:byte; begin appCount := nVal; for i:=1 to nVal do appSigma[ i ]:=Sigma[ i ]; if( appType = apSpline )then Spline( appCount, appSigma, siB, s iC, siD); end; procedure getApproximationData( var SIG ; var N : byte ); var Sigma : Parameters absolute SIG; i:byte; begin N := appCount; for i:=1 to appCount do Sigma[ i ]:=appSigma[ i ]; end; procedure setApproximationItem( SIG:real ; N : byte ); begin appSigma[ N ] := SIG; if( appType = apSpline )then Spline( appCount, appSigma, siB, siC, siD); end; procedure functionFi( x:real ; var result:complex ); get boundary conditions function value var beta : array[ 1..maxPAR ]o f real; q : array[ 1..maxPAR ]of complex; fi : array[ 0..maxPAR ]of complex; th , z1 , z2 , z3 , z4 , z5 , z6 , z7 : complex; i : byte; begin mkComp( 0, 0, fi[0] ); with cInfo do for i:=1 to Nlay do begin beta[i]:=R*sqrt( w*mu0*sigma[i] ); calculation of beta mkComp( sqr(x), sqr(beta[i])*m[i], z7 ); calculation of q, z7=q^2 SqrtC( z7, q[i] ); mulComp( q[i], b[i], z6 ); calculation of th,z6=q*b tanH( z6, th ); th=tanH(q*b) mkComp( sqr(m[i])*sqr(x), 0, z6 ); z6=m2x2 SubC( z6, z7, z5); z5=m2x2-q2 AddC( z6, z7, z4); z4=m2x2+q2 MulC( z5, th, z1); z1=z5*th MulC( z4, th, z2); z2=z4*th mulComp( q[i], 2*x*m[i], z3 ); z3=2xmq SubC( z2, z3, z4 ); MulC( z4, fi[i-1], z5 ); SubC( z1, z5, z6 ); z6=high AddC( z2, z3, z4 ); MulC( z1, fi[i-1], z5 ); SubC( z4, z5, th ); th=low DivC( z6, th, fi[i] ); end; eqlComp( result, fi[ cInfo. Nlay ] ); end; procedure funcSimple( x:real; var result:complex ); intergrand function value var z : complex; begin with cInfo do begin functionFi( x, result ); mulComp( result, exp( -x*H ), z ); mulComp( z, J1( x*Kr ), result ); mulComp( result, J1( x/Kr ), z ); eqlComp( result, z ); end; end; procedure funcMax( x:real; var result:complex ); max value ; When Fi[Nlay]=1 var z1, z2 : complex; begin with cInfo do begin mkComp( 1,0,z1 ); mulComp( z1, exp(-x*H), z2 ); mulComp( z2, J1( x*Kr ), z1 ); mulComp( z1, J1( x/Kr ), result ); end; end; procedure integralBS( func:procFunc ; var result:complex ); integ ral by Simpson var z , y , tmp : complex; hh : real; i, iLast : word; begin with cInfo do begin hh:=xMax/steps; iLast:=steps div 2; nullComp(tmp); func( 0, z ); eqlComp( y, z ); for i:=1 to iLast do begin func( ( 2*i-1 )*hh , z ); deltaComp( tmp, z ); end; mulComp( tmp, 4, z ); deltaComp( y, z ); nullComp( tmp ); iLast:=iLast-1; for i:=1 to iLast do begin func( 2*i*hh ,z ); deltaComp( tmp, z ); end; mulComp( tmp, 2, z ); deltaComp( y, z ); func( xmax, z ); deltaComp(y,z); mulComp( y, hh/3, z ); eqlComp( result, z ); end; end; I = h/3 * [ F0 + 4*sum(F2k-1) + 2*sum(F2k) + Fn ] procedure integral( F:procFunc; var result:complex ); integral with given error var e , e15 : real; flag : boolean; delta , integ1H , integ2H : complex; begin with cInfo do begin e15 :=eps*15; | I2h-I1h |/(2^k -1 )< Eps ; k=4 for Simpson method steps:=20; flag :=false; integralBS( F, integ2H ); repeat eqlComp( integ1H, integ2H ); steps:=steps*2; integralBS( F, integ2H ); SubC( integ2H, integ1H, delta ); e:=Leng( delta ); if( emaxsteps) ); if( flag )then begin eqlComp( result, integ2H ); end else begin writeln('Error: Too big number of integration steps.'); halt(1); end; end; end; procedure initConst( par1, par2 : integer; par3, par4 : real; par5:boolean ); var i : byte; bThick, dl, x : real; const Ri=0.02; hi=0.005; radius and lift-off of excitation coil Rm=0.02; hm=0.005; radius and lift-off of measuring coil begin with cInfo do begin Nlay :=par1; xMax :=par3; max steps:=par2; R :=sqrt( Ri*Rm ); H :=( hi+hm )/R; Kr :=sqrt( Ri/Rm ); eps :=par4; bThick :=hThick*0.002/R; 2*b/R [m] for i:=1 to Nlay do m[i]:= mu[i]; if par5 then begin bThick:=bThick/NLay; for i:=1 to Nlay do b[i]:=bThick; dl:=1/NLay; x:=dl/2; x grows up from bottom of slab to the top for i:=1 to Nlay do begin zCentre[i]:=x; x:=x + dl; end; end else for i:=1 to Nlay do begin b[i]:=( zLayer[i+1]-zLayer[i] )*bThick; zCentre[i]:=( zLayer[i+1]+zLayer[i] )/2; end; end; end; procedure init( f:real ); get current approach of conductivity values var i : byte; begin with cInfo do begin w:=PI23*f; for i:=1 to Nlay do sigma[i]:=getSiFunction( zCentre[i] )*1E6; end; end; procedure getVoltage( freq : real ; var ur,ui : real ); var U, U0, Uvn, tmp : complex; begin init( freq ); integral( funcSimple, U ); U =Uvn integral( funcMax , U0 ); U0=Uvn max divComp( U, Leng(U0), Uvn ); Uvn=U/|U0| mkComp( 0, 1, tmp ); tmp=( 0+j1 ) MulC ( tmp, Uvn, U ); U= j*Uvn = Uvn* ur := U.re; ui := U.im; end; END. П 1 . 8 Модуль решения обратной задачи ВТК для НВТП EMinimum Unit EMinimum; INTERFACE Uses EData, Crt, EFile, EDirect; procedure doMinimization; IMPLEMENTATION procedure getFunctional( Reg : byte ); var ur, ui, dur, dui, Rgt : real; ur2, ui2: Functionals; i, j, k : byte; begin getApproximationData( si , k ); setApproximationData( Rg, mCur ); case Reg of 0 : for i:=1 to nFreqs do get functional F value begin getVoltage( freqs[i], ur, ui ); Uer[ i ]:=ur; we need it for dU Uei[ i ] :=ui; Fh[1,i] := SQR( ur-Umr[i] ) + SQR( ui-Umi[i] ); end; Right:U'(i)= (U(i+1)-U(i))/h 1 : for i:=1 to mCur do get dF/dSI[i] value begin Rgt:= Rg[i]*( 1+incVal ); si[i]=si[i]+dsi[i] setApproximationItem( Rgt, i ); set new si[i] value for j:=1 to nFreqs do begin get dUr/dSI,dUi/dSI getVoltage( freqs[ j ], ur, ui ); dur:=( ur-Uer[j] )/( Rg[i]*incVal ); dui:=( ui-Uei[j] )/( Rg[i]*incVal ); Fh[i,j]:=2*(dur*(Uer[j]-Umr[j])+dui*(Uei[j]-Umi[j])); end; setApproximationItem( Rg[i], i ); restore si[i] value end; Central:U'(i)= (U(i+1)-U(i-1))/2h 2 : for i:=1 to mCur do get dF/dSI[i] value begin Rgt:=Rg[i]*( 1+incVal ); si[i]=si[i]+dsi[i] setApproximationItem( Rgt, i ); set new si[i] value for j:=1 to nFreqs do getVoltage( freqs[j],ur2[j],ui2[j] ); Rgt:=Rg[i]*( 1-incVal ); si[i]=si[i]-dsi[i] setApproximationItem( Rgt, i ); set new si[i] value for j:=1 to nFreqs do begin get dUr/dSI,dUi/dSI getVoltage( freqs[ j ], ur, ui ); dur:=( ur2[j]-ur )/( 2*Rg[i]*incVal ); dui:=( ui2[j]-ui )/( 2*Rg[i]*incVal ); Fh[i,j]:=2*(dur*(Uer[j]-U mr[j])+dui*(Uei[j]-Umi[j])); end; setApproximationItem( Rg[i], i ); restore si[i] value end; end; setApproximationData( si , k ); end; procedure doMinimization; const mp1Max = maxPAR + 1; mp2Max = maxPAR + 2; m2Max = 2*( maxPAR + maxFUN ); m21Max = m2Max + 1; n2Max = 2*maxFUN; m1Max = maxPAR + n2Max; n1Max = n2Max + 1; mm1Max = maxPAR + n1Max; minDh : real = 0.001; criterion of an exit from golden section method var A : array [ 1 .. m1Max , 1 .. m21Max ] of real; B : array [ 1 .. m1Max] of real; Sx: array [ 1 .. m21Max] of real; Zt: array [ 1 .. maxPAR] of real; Nb: array [ 1 .. m1Max] of integer; N0: array [ 1 .. m21Max] of integer; a1, a2, dh, r, tt, tp, tl, cv, cv1, cl, cp : real; n2, n1, mp1, mp2, mm1, m1, m2, m21 : integer; ain : real; apn : real; iq : integer; k1 : integer; n11 : integer; ip : integer; iterI : integer; i,j,k : integer; label 102 ,103 ,104 ,105 ,106 ,107 ,108; begin n2:=2*nFreqs; n1:=n2+1; m1:=mCur+n2; mp1:=mCur+1; mp2:=mCur+2; m m1:=mCur+n1; m2:=2*( mCur+nFreqs ); m21:=m2+1; for k:=1 to m1Max do for i:=1 to m21Max do A[k,i]:=0; iterI:=0; 102: iterI:=iterI+1; getFunctional( 0 ); for i:=1 to nFreqs do b[i]:= -Fh[1,i]; getFunctional( derivType ); for k:=1 to mCur do begin Zt[k]:=Rg[k]; for i:=1 to nFreqs do begin A[i,k+1]:=Fh[k,i]; A [i+nFreqs,k+1]:=-A[i,k+1]; end; for i:=1 to nFreqs do B[i]:=B[i]+Rg[k]*A[i,k+1]; end; for i:=1 to nFreqs do B[i+nFreqs]:=-B[i]; for i:=n1 to m1 do B[i]:=Gr[1,i-n2]-Gr[2,i-n2]; for i:=1 to m1 do begin if i<=n2 then for k:=2 to mp1 do B[i]:=B[i]-A[i,k]*Gr[2,k-1]; A[i,1]:=-1; if i>n2 then begin A[i,1]:=0; for k:=2 to mp1 do if i-n2=k-1 then A[i,k]:=1 else A[i,k]:=0; end; for k:=mp2 to m21 do if k-mp1=i then A[i,k]:=1 else A[i,k]:=0; end; k1:=1; for k:=1 to n2 do if B[k1]>B[k] then k1:=k; for k:=1 to mp1 do A[k1,k]:=-A[k1,k]; A[k1,mCur+1+k1]:=0; B[k1]:=-B[k1]; for i:=1 to n2 do if i<>k1 then begin B[i]:=B[i]+B[k1]; for k:=1 to mm1 do A[i,k]:=A[i,k]+A[k1,k]; end; for i:=mp2 to m21 do begin Sx[i]:=B[i-mp1]; Nb[i-mp1]:=i; end; for i:=1 to mp1 do Sx[i]:=0; Sx[1]:=B[k1]; Sx[mp1+k1]:=0; Nb[k1]:=1; 103: for i:=2 to m21 do N0[i]:=0; 104: for i:=m21 downt o 2 do if N0[i]=0 then n11:=i; for k:=2 to m21 do if ((A[k1,n11]N0[k])) then n11:=k; if A[k1,n11]<=0 then goto 105; iq:=0; for i:=1 to m1 do if i<>k1 then begin if A[i,n11]>0 then begin iq:=iq+1; if iq=1 then begin Sx[n11]:=B[i]/A[i,n11]; ip:=i; end else begin if Sx[n11]>B[i]/A[i,n11] then begin Sx[n11]:=B[i]/A[i,n11]; ip:=i; end; end; end else if iq=0 then begin N0[n11]:=n11; goto 104; end; end; Sx[Nb[ip]]:=0; Nb[ip]:=n11; B[ip]:=B[ip]/A[ip,n11]; apn:=A[ip,n11]; for k:=2 to m21 do A[ip,k]:=A[ip,k]/apn; for i:=1 to m1 do if i<>ip then begin ain:=A[i,n11]; B[i]:=-B[ip]*ain+B[i]; for j:=1 to m21 do A[i,j]:=-ain*A[ip,j]+A[i,j]; end; fo r i:=1 to m1 do Sx[Nb[i]]:=B[i]; goto 103; 105: for k:=1 to mCur do Sx[k+1]:=Sx[k+1]+Gr[2,k]; a1:=0; a2:=1.; dh:=a2-a1; r:=0.618033; tl:=a1+r*r*dh; tp:=a1+r*dh; j:=1; 108: if j=1 then tt:=tl else tt:=tp; 1 06: for i:=1 to mCur do Rg[i]:=Zt[i]+tt*(Sx[i+1]-Zt[i]); getFunctional( 0 ); cv:=abs(Fh[1,1]); if nFreqs>1 then for k:=2 to nFreqs do begin cv1:=abs(Fh[1,k]); if cvcp then begin a1:=tl; dh:=a2-a1; tl:=tp; tp:=a1+r*dh ; tt:=tp; cl:=cp; j:=4; end else begin a2:=tp; dh:=tp-a1; tp:=tl; tl:=a1+r*r*dh; tt:=tl; cp:=cl; j:=3; end; goto 106; 107: if (iterI < iterImax)AND(NOT saveResults( nStab,iterI )) the n goto 102; end; End. Приложение 2 - Удельная электрическая проводимость материалов Приведем сводку справочных данных согласно [7-9] . Материал min , [ МСм / м ] max , [ МСм / м ] Немагнитные стали 0.4 1.8 Бронзы (БрБ , Бр 2, Бр 9) 6.8 17 Латуни (ЛС 59, ЛС 62) 13.5 17.8 Магниевые сплавы (МЛ 5-МЛ 15) 5.8 18.5 Титановые сплавы (ОТ 4, ВТ 3-ВТ 16) 0.48 2.15 Алюминиевые сплавы (В 95, Д 16, Д 19) 15.1 26.9 Приложение 4 - Abst act The i n verse eddy current problem can be described as the task of reconstructing an unknown distribution of electrical conductivity from eddy-current probe voltage measurements recorded as function of excitation frequency. Conductivity variation may be a result of surface processing with substances like hydrogen and carbon or surface heating. Mathematical reasons and supporting softw are for inverse conductivity profiling were developed by us. Inverse problem was solved for layered plane and cylindrical conductors. Because the inverse problem is nonlinear, we propose using an iterative algorithm which can be formalized as the minimizat ion of an error functional related to the difference between the probe voltages theoretically predicted by the direct problem solving and the measured probe voltages. Numerical results were obtained for some models of conductivity distribution. It was show n that inverse problem can be solved exactly in case of correct measurements. Good estimation of the true conductivity distribution takes place also for measurement noise about 2 percents but in case of 5 percent error results are worse.

© Рефератбанк, 2002 - 2024