рефераты рефераты
Домой
Домой
рефераты
Поиск
рефераты
Войти
рефераты
Контакты
рефераты Добавить в избранное
рефераты Сделать стартовой
рефераты рефераты рефераты рефераты
рефераты
БОЛЬШАЯ ЛЕНИНГРАДСКАЯ БИБЛИОТЕКА
рефераты
 
МЕНЮ
рефераты Решение обратной задачи вихретокового контроля рефераты

БОЛЬШАЯ ЛЕНИНГРАДСКАЯ БИБЛИОТЕКА - РЕФЕРАТЫ - Решение обратной задачи вихретокового контроля

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

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 рефераты