 |
 |
 |
|
|||||||||
|
 | |||||||||||
|
 |
||||||||
|
МЕНЮ
|
БОЛЬШАЯ ЛЕНИНГРАДСКАЯ БИБЛИОТЕКА - РЕФЕРАТЫ - Решение обратной задачи вихретокового контроляРешение обратной задачи вихретокового контроля9. Результаты численного моделирования 9.1 Аппроксимации при численном моделировании Для построения моделей реальных распределений ЭП возможно применение целого ряда аппроксимаций. Все они могут быть разделены на два класса. 1. Аппроксимации, строящиеся по набору из произвольного числа узлов. Наиболее распространенные из них: кусочно-постоянная, кусочно-линейная и сплайном. В условиях нашей задачи указанные аппроксимации имеют несколько существенных недостатков: · Результаты аппроксимаций слабо согласуются с реальностью. Кусочно-постоянная и кусочно-линейная аппроксимации принципиально являются негладкими, а аппроксимация сплайном сглаживает все, в результате чего возникают значительные немонотонности и всплески. · При увеличении количества узлов аппроксимации быстро нарастает неустойчивость процесса решения обратной задачи, для противодействия которой требуется применение искусственных приемов, не гарантирующих успеха. · В реальных условиях мы не имеем достоверной априорной информации о величине ЭП в узлах аппроксимации, расположенных в глубине пластины. 2. Аппроксимации, строящиеся по значениям ЭП на верхней и нижней поверхностях пластины и нескольким параметрам аппроксимации. Наиболее известные из них: экспоненциальная, гиперболическим тангенсом и гауссоидой. Аппроксимации имеют вид: | |||||
|
- аппроксимация экспоненциальная |
||||||
|
- аппроксимация гиперболическим тангенсом |
||||||
|
- аппроксимация гауссоидой |
где
x
- координата, равна нулю на нижней поверхности пластины и единице на верхней
s1
- величина электропроводности на верхней поверхности пластины
s2
- величина электропроводности на нижней поверхности пластины
a
- коэффициент, характеризующий крутизну экспоненты
b
- коэффициент
g
- коэффициент, характеризующий крутизну; g=0 соответствует случаю слоя с проводимостью s1 и толщиной b на полупространстве с проводимостью s2
d
- коэффициент, характеризующий крутизну
Для нашей задачи подобные аппроксимации являются предпочтительными, поскольку обладают заметными достоинствами:
· Аппроксимации являются монотонными и гладкими, что хорошо согласуется с физической реальностью.
· Пользуясь физически обоснованными рассуждениями мы можем получить необходимую априорную информацию о величинах ЭП в приповерхностных слоях пластины.
· Процесс решения обратной задачи существенно более устойчив и осуществляется значительно быстрее
Для иллюстрации наших рассуждений приведем пример применения приведенных выше аппроксимаций к случаю восстановления кусочно-линейной функции. По оси абсцисс отложена относительная глубина, по оси ординат электропроводность (МСм/м).
На графике показаны аппроксимации: кусочно постоянная(SIci),кусочно линейная(SIli), сплайн(SIs), экспоненциальная(SIe), гиперболическим тангенсом (Sith), гауссоидой(SIg).
Легко заметить, что аппроксимация гиперболическим тангенсом хорошо описывает приповерхностные изменения (аналогично экспоненциальной при большом показателе экспоненты). Гауссоида может быть легко воспроизведена с помощью экспоненциальной аппроксимации, поэтому в дальнейшем использована не будет.
9.2 Модели реальных распределений электропроводности
Модель задачи должна описывать некоторую пластину, подвергнутую поверхностной обработке. Для определенности зададим толщину пластины равной двум сантиметрам. На основе данных из Приложения 2 зададим значения ЭП вблизи нижней и верхней поверхностей соответственно 20 (МСм/м) и 13 (МСм/м).
Для решения обратной задачи необходимо задать априорную информацию о величине ЭП в узлах аппроксимации. В качестве таковой примем интервал [8,25] (МСм/м), полученный внесением 25% отклонения от считаемых истинными значений. Это отклонение моделирует неточность априорной информации.
Из-за особенностей реализации алгоритма устойчивость решения сильно зависит от точности задания ЭП в узле, соответствующем нижней поверхности пластины, поэтому ограничение в нем зададим интервалом [19,21] (МСм/м).
В нашем случае все возможные модели распределений ЭП могут быть разделены на два класса. Распределения относящиеся к первому из них условно назовем глубинными. В них ЭП претерпевает существенные изменения на протяжении всей глубины пластины. Второй класс образуют распределения, ЭП в которых заметно изменяется лишь в приповерхностном слое глубиной порядка четверти пластины., поэтому назовем эти распределения поверхностными.
Критерием отличия восстановленной функции распределения ЭП от модельной будем считать величину относительной погрешности, поскольку сравнение результатов с ее помощью вполне адекватно целям нашей работы.
Следует отметит, что погрешность восстановления для поверхностных распределений ЭП представляет практический интерес в области, примерно ограниченной глубиной порядка четверти пластины, что обусловлено физическим смыслом поверхностной обработки. Поэтому для случаев поверхностных распределений основное внимание будем уделят именно указанным глубинам.
Для проверки возможности восстановления приповерхностных изменений ЭП рассмотрим две базовые модели поверхностных распределений.
Базовая модель A1.
Аппроксимация экспонентой.
Проводимость s2=20 [МСм/м]
Проводимость s1=13 [МСм/м]
Показатель экспоненты a ={ 25, 38, 120, 200 }.
Базовая модель A2
Аппроксимация гиперболическим тангенсом.
Проводимость s2=20 [МСм/м]
Проводимость s1= 6 [МСм/м]
Коэффициент b = 1
Коэффициент g = { 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} [МСм/м].
Заметим, что на практике можно осуществить достаточно точное определение величины ЭП приповерхностных слоев путем измерений проводимости традиционными средствами, поэтому дополнительно рассмотрим модельные задачи при условии известной ЭП на верхней, а так же верхней и нижней поверхностях.
Поскольку на практике результаты измерений вносимого напряжения имеют определенную погрешность, все модели будем рассчитывать эмулируя погрешность dU= 0,1,2,5%.
Для исследования зависимости результатов восстановления распределений ЭП от частоты возбуждения разобьем частотный диапазона три части следующим образом( глубины проникновения приведены для случая постоянной ЭП s=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
100
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 profiling., 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², 1968, vol39, ?3, p2829-2838
43." General analysis of probe coils near stratified conductors " , Deeds WE, Dodd 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, ²Journal of Applied 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 presence of inhomogeneities: a new formalism " , Nair SM, Rose JH, ²Journal of Applied Phisics², 1990, vol68, ?12, p5995-6009
48." Eddy-current probe 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 permeability: 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 solution 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
- количество частот возбуждения гармоник вносимой относительной ЭДС
nStab
- число стабилизируемых значащих цифр
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 качестве последней строки помещаются nPoints вычисленных значений э/проводности в равномерно расположенных узлах пластины.
П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 := nApprox 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}
si0[ 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 not 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 inverse 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 overflow}
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-friendly 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;
apExp = 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 : integer; {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/max
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 : string; {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( FF, outFileName ); {clear output file}
rewrite( FF );
close( FF );
assign( FF, inFileName ); {read input file}
reset( FF );
readln( FF );
readln( FF );
readln( FF, hThick, 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 );
if ( 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 : text;
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);
writeln( 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;
writeln;
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:=sec2 - 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 multilayer 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 initConst( 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 approximation}
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]-appSigma[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
writeln('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, siC, 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 ]of 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 );{integral 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
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( e<e15 )then flag:=true;
until( ( flag ) OR (steps>maxsteps) );
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;
maxsteps:=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 ); 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]-Umr[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; mm1:=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 downto 2 do
if N0[i]=0 then n11:=i;
for k:=2 to m21 do
if ((A[k1,n11]<A[k1,k]) AND (k<>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;
for 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;
106:
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 cv<cv1 then cv:=cv1;
end;
if (j=1) or (j=3) then cl:=cv
else cp:=cv;
if j=1 then
begin
j:=2;
goto 108;
end;
if dh<MinDh then goto 107;
if cl>cp 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 )) then goto 102;
end;
End.
Приложение 2 - Удельная электрическая проводимость материалов
Приведем сводку справочных данных согласно[7-9].
Материал
smin ,[МСм/м]
smax ,[МСм/м]
Немагнитные стали
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 - Abstract
The inverse 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 software 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 minimization 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 shown 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.
 |
НОВОСТИ | |
| Изменения | ||
| Прошла модернизация движка, изменение дизайна и переезд на новый более качественный сервер | ||
|
СЧЕТЧИК | |
| БОЛЬШАЯ ЛЕНИНГРАДСКАЯ БИБЛИОТЕКА | ||
 |
© 2010 | |