Определение распределения электрического сопротивления в токопроводящей пленке методом томографии Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 517.54

В.Д. Садков, А.В. Лопаткин

ОПРЕДЕЛЕНИЕ РАСПРЕДЕЛЕНИЯ ЭЛЕКТРИЧЕСКОГО СОПРОТИВЛЕНИЯ В ТОКОПРОВОДЯЩЕЙ ПЛЕНКЕ МЕТОДОМ ТОМОГРАФИИ

Нижегородский государственный технический университет им. Р.Е. Алексеева

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

Ключевые слова: импедансная томография, токопроводящая пленка, распределение электрического сопротивления.

Введение

Предложенный в работе [1] алгоритм определения зависимости электрического сопротивления токопроводящей пленки от координат базируется на методе томографии [2-5], учитывающем конкретные свойства рентгеновских лучей. Линии тока в токопроводящей среде принципиально отличаются по своим свойствам от рентгеновских лучей: их траектории не прямолинейны (в том числе и в однородном диске с контактами по периметру, выбранном в [1]), зависят от неизвестного распределения электрического сопротивления и не могут быть установлены до эксперимента; величина сопротивления R, определяется всей совокупностью линий тока, протекающих между контактами i и^трубкой тока), и несет информацию о распределении электрического сопротивления пленки вдоль плоской фигуры, образуемой трубкой тока с ненулевой шириной 5 в любом поперечном сечении (при 5 = 0 контакты i и j точечные и Ry-бесконечно велико). Перечисленные особенности не учтены в [1]. Задача (даже без учета необходимости ее постановки в произвольной двумерной области) не может быть решена предложенным алгоритмом, ибо в определении нуждаются не только распределение электрического сопротивления вдоль совокупности траекторий , но и сами траектории, зависящие от этого неизвестного распределения.

Математическая постановка задачи

Задача в данном случае заключается в определении гладкой или кусочно-гладкой неотрицательной функции p(x,y), удовлетворяющей в токопроводящей области (рис. 1) следующей системе внутренних обратных (инверсных) краевых задач [6, 7]:

1

V

Р( x, у)

VU (x, у)

=0 (1)

UGk, = Ei , i = 1, 2 •••,n

Jk ,i

\UG = 0 , j = 1, 2, ...,i-1, i + 1, •••,n , (2)

1 k,J

f - 0 ( 3)

on G

1 r 1 OU Ii, j 1

— I--dl - —x- --, (4)

Et Jp(x, y) On Et Ri j ^

k Gk ,i

© Садков В.Д., Лопаткин А.В., 2015.

где р(х, у) - искомое распределение электрического сопротивления; п - число позиций контактов или зондов по периметру исследуемой области; т - число контактов (2 <т<п), устанавливаемых на п позиций; Ок,г и 4,— участок границы, занимаемый г-м контактом, и его длина; О - граница области, свободная от контактов; Е, - разность потенциалов между контактами г и у; Ягу - сопротивление между контактами г и у, вычисленное по измеренному току 1гу.

а., 8г а, ¿1 атч С, ат 4,

Рис. 1

Минимальная размерность системы (1) - (4) Ь=п-1 имеет место при т=п, максимальная - Ь=п(п - 1)/2 при т=2. Рост размерности Ь с уменьшением т сопровождается снижением сложности каждой отдельной задачи системы. Отметим, что при любом количестве контактов 2<т<п общее количество независимых измерений сопротивлений Ягу равно п(п - 1)/2 (при т<п контакты перемещаются по п выбранным позициям).

Система (1)-(4) внутренних обратных краевых задач отличается от известных нетрадиционностью постановки: в качестве исходных данных задаются не значения потенциалов поля [6, 7], а потоки поля через участки границы с контактами. Методы решения таких систем в литературе не рассмотрены.

Система (1)-(4) может быть решена конечно-разностным методом путем варьирова-

2 4

ния параметров рк,; (к= 1, 2, ..., М; 1= 1, 2, ..., Ы; 10 <МЫ<10 ) электрического сопротивления ячеек сетки для удовлетворения условий (4) [5]. Алгоритм такого варьирования, последовательно приближающий параметры среды к искомому распределению, построим на базе алгебраического алгоритма реконструкции [2], устойчивого к малым изменениям входных данных:

М N м N

=Е£ркХ«кЛ, У"7 =ЕЕРк,/«К, (5)

к=\ I=1 к=11=\

где М, N - количество элементов изображения (элизов) по координате х и у соответственно; рк,; - электрическое сопротивление элиза с индексом к, I; ¡=1, 2, ..., п иу=1, 2, ..., ¡-1, г +1, ..., п (с учетом Ягу= Яу,г достаточно ограничиться значениями у= ¡+1, г +2, ..., п); к, - весовая функция для г, у измерения, заменяющая реальную трубку эквивалентной с равномерным распределением тока по сечению путем задания значений кгу =1 и кгу = 0 внутри и вне трубки; п^7 - коэффициент формы участка эквивалентной трубки в элизе к, I при г,у измерении.

Представление величины Я, в форме (5) возможно лишь при условии, что ширина эквивалентной трубки тока 5 меньше размера элиза Н вблизи контакта 5^4 и, значит, 4<Н. В случае узких (5<<Н) трубок тока можно считать п\\ = (длина траектории линии тока в трубке с условной единичной шириной в элизе к, I при г, j измерении) и

мы

Я, = т.т.РкА1 , (6)

''7 к=1=1

Методика решения системы (1) - (4) заключается в следующем:

1. Исследуемую область разделяем на элизы равномерной сеткой с шагом Н.

2. Проводим, в соответствии с рис. 1, п=(п-1)/2 измерений и формируем матрицу Я сопротивлений Я^.

3. В приближении однородной среды с помощью конечно-разностного метода для каждой из Ь краевых задач системы (1) - (4) находим структуру поля потенциалов, формиру-

-(0)

ем исходную проекционную матрицу электрических сопротивлений в элизах р .

4. Определяем геометрию эквивалентных трубок тока. Формируем исходную матрицу

-(0) —(0) коэффициентов формы п или (для узких трубок) длин траекторий и .

5. Проводим цикл реконструкции по указанным в [1, 2] формулам ((30),(31) настоящей работы), используя соотношения (5) или (6). При этом не стремимся к получению результата с высокой точностью ( проводим цикл из 3-5 итераций), учитывая повторение ( после уточнения свойств среды в п.6) циклов реконструкции.

6. Проверяем выполнение условия окончания процесса

Р Р ~1

Ры ~Ррл

<8 ,к=1, 2, ...№; 1=1, 2, ..., N (7)

где рР I, рР;/1 - электрическое сопротивление элиза к, I нар-1 и р циклах реконструкции; 8 -

малая величина, имеющая порядок точности измерений токов /у.

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

—(0)

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

Выбор схемы построения системы сеточных уравнений

Для построения системы сеточных уравнений в среде, свойства которой в каждой ячейке сетки могут быть различны (рис. 2), используем уравнение для поля в интегральной форме [8]. Для всех внутренних узлов области (нумерация слева направо и сверху вниз) получаем однотипные уравнения с погрешностью порядка 0(И ):

кр-\4=(рр,д-\+ар,^12 ; кр^-1=(ар^-1+ар+1^_1)/2 ; кр,д+1=(«р+1+«р,д)/2; </-= 1/Р— <3р,д=1/ррд,' кр+1>д=(Ср+11д-1+Ср+11д)/2.

Погрешность порядка 0(Н3) обеспечивается решением задачи на последовательности сеток с шагами Н1=Н/4и Н2=И/8 (Н—размер элиза). При этом не исключается возможность использования неравномерных сеток. Краевое условие (3) реализуется введением "фиктивных" точек вне исследуемой области на расстоянии Н от ее границы и приравниванием потенциалов "фиктивных" точек потенциалам соседних граничных точек. Система сеточных уравнений решается методом верхней релаксации.

Рис. 2

Получение начального приближения для поля потенциалов

Аналитическое начальное приближение в предположении однородной среды для поля потенциалов двумерной области (рис. 1) с минимальными затратами памяти и машинного времени получаем методом теории функций комплексного переменного путем отображения исследуемой области в плоскости Z на верхнюю полуплоскость Z1 , поле потенциалов в которой для каждой из Ь краевых задач иЛх1 у1) находится так [9]:

и г (уЛ = Яе Г

п « - К-ь)

г =1

-0.5

т-1 г=1

(т-1)-гж + А

(9)

гдеа7, Ьг — координаты начал и концов контактов в плоскости Z1; т — общее число контактов.

Задача отображения многоугольной области решается с помощью интеграла Кристоф-феля—Шварца, односвязной или двухсвязной области произвольной формы — методом тригонометрической интерполяции [10]. Постоянные Ак (к=1, 2, т-1) в (9) определятся из удовлетворения краевых условий (2) на контактах. Получающаяся система разрешима и приводится к виду

1 • А = V , (10)

где

1 =

11,т-2 Г 2 т-2

Г\т-3 Г 2.т-3

Г т-1 т-3

11.0

Г 2,0 1

(11)

А — вектор-столбец коэффициентов [А1, А2, ..., Ат-1]; V — вектор-столбец разностей потенциалов контактов [и2-и1,из-и2,..ит-ит-1].

п

Элементы Fk,i матрицы F - гиперэллиптические интегралы вида

a к+1 m

Fk,i = (-1)к J Ы

-0,5

п (t - ai )(t -bi)

¿=i

t dt , (12)

вычисляемые с использованием квадратурных формул наивысшей степени точности после приведения интервала интегрирования к стандартному (-1,1) и выделения особенностей вида

2 —0 5

(1-х ) ' в весовую функцию [11]. Константа Ат находится по (9) при ^ь(х1, у1) = ик и интегрировании в пределах от 0 до ак для любого к. По соотношениям (9) определяем начальное

—(0)

приближение и для потенциалов сеточных узлов в исходной плоскости Z после вычисления их координат в верхней полуплоскости Z1.

Получение начального приближения для проекционной матрицы

По известной структуре поля потенциалов (9), полученной в предположении однородной среды, рассчитываем коэффициенты формы п^г сопротивлений Я, [12]:

Пк,1г = Хк,1, Пк,рг= (Хк.р-Хк.р-Х), Пк,тг= хк,т-1 ,Хк,д= det[ Р Ек,д ] , (13)

где к=1, 2, ..., т; р=2, 3, ..., т-1; ц=1, 2, ..., т; т - число контактов по периметру структуры; det[ Р ] - определитель матрицы F , описанной в (11), (12); det[ Ек Л ] - определитель матрицы Екполучаемой из матрицы Р путем замены ее ц-й строки (^д т-2 ^Чтт-2 , ..., Fq, о) строкой с элементами (Ек,т-2, Ек,т-2, ..., Ек,о). Элементы Ек I имеют вид

b

п (t - ai )(t-bi)

i=i

-0,5

tldt

(14)

Ек ,i = (-1)k-1 J

ак

и вычисляется аналогично элементам Fkj(12).

-(0)

Начальное приближение ро для проекционной матрицы р получаем следующим образом. В предположении однородности среды в пределах трубок тока определяем удельные сопротивления p^ для каждой из n(n-1)/2 трубок :

plj=RlJlnlJ. (15)

Значение ро получаем усреднением величин р^-по (15):

^ n-1 n

p=nn-i)72 . (16)

Очевидно Pijmin<P0 <p ijmax , где pi',jmin> pi,jmax — минимальное и максимальное из значений pij по (15).

Разрешающая способность и чувствительность

По-видимому, необходимо ввести понятия предельной и реальной разрешающих способностей. Предельная разрешающая способность определяется средней шириной 5 трубок тока и зависит от геометрии области, количества, взаимного расположения и длины контак-тов4. Для заданной области, при выбранном способе расположения контактов известной длины по ее периметру, параметр 5 зависит только от количества контактов и легко вычисляется по соотношениям (12) - (15) через усредненный коэффициент формы ns:

^ n-1 n-1

5 = lij max/ns , ns = ZZ ni,jr , (17)

n(n - 1)/ 2 i=1 j=i+1

где /¿¿тах - наибольшее расстояние между контактами г и¿.

Например, для прямоугольной структуры площадью £=7,5х 11,25мм2 с 30 контактами (п=т=30) длиной /к=0,2 мм, равномерно расположенными по периметру, п=232, /¿¿тах=21,6 мм и 5=0,1 мм. Реальная разрешающая способность всегда ниже и определяется минимальным размером элиза Ят;п=(^/(М¥))0,5=0,4 мм, где М¥<п(п-1)/2=435. Размер Н=0,75 мм>Нт;п , соответствующий целым М и N (при равномерном расположении контактов М^=7,5/11,25 мм, ^1,5М, М^1,5М ) получаем с учетом неравенств Н>5 и Н>(3...4)/к. Для прямоугольной структуры площадью £=5х10 мм с т=2, п=30 и 4=0,5 мм учет простейшего условия Нтщ>(3...4)/к =1,25мм (MхN=4х8=32) показывает, что проведение полного цикла из п(п-1)/2=435 измерений бесполезно. Предварительный анализ сразу выявляет низкую разрешающую способность (п^=2,4 ,/у-тах=11,7 мм) и из условия Нтш>5 4,9 следует, что в структуре может быть не более двух элизов.

В этих условиях необходимо увеличить предельную разрешающую способность, что, как показывают соотношения (12) - (15), возможно только за счет увеличения количества контактов т. Отметим, что снизить трудоемкость решения задачи, обеспечив одновременно высокую разрешающую способность, по крайней мере, в некоторой зоне области (особенно при повторном решении задачи с целью уточнения результата), возможно только при

2 4

уменьшении общего количества элизов до приемлемой величины (10 <MN<10 ) за счет использования неравномерной сетки.

При среднем значении сопротивления между контактами ¿¿Яу=п -р0 изменение сопротивления АЯу, определяющее требования к точности (чувствительности) измерительных приборов, связано с изменением удельного сопротивления Арк,/ в элизе к, / следующим образом: АЯу =Арк,/хН/5 .

Оценка степени неоднородности структуры

Гарантированную оценку р(х, у) сверху (ртах) и снизу (рт;п) находим для наихудшего случая (при оценке ртах считаем, что пленка с р=ртах занимает только один из элизов, через которые проходит ¿, ] - трубка; в остальных элизах - пленка с р=р0):

[1-(/у/аН)(1-ргУт1п/рс)] =ртш<р(х, у)<ртах=

=ро [1+$уН)(рушах/ро-1)], (18)

гдергутт, рг,7тах - минимальное и максимальное значения р, по (15); /, - расстояние между контактами г и ] в структуре; а > 1 - количество элизов с р=ртщ, выбираемое из условия положительности квадратных скобок в (18 ).

Определение геометрии эквивалентных трубок тока

Для решения задачи (рис. 1) достаточно определить по известному на равномерной сетке с шагом к полю потенциалов граничные траектории линий тока. Конечными токами этих траекторий (для трубки между контактами г и ] ) являются точки А и В контакта ¿; направление траекторий в целом - на контакт ¿, в каждой ячейке сетки - в сторону градиента потенциала. В таких границах протекает практически весь ток реальной трубки, граничные линии которой (в отличие от эквивалентной трубки) идут из точек А и В вдоль границы структуры в сторону роста потенциала до точек ф иу с наибольшим потенциалом на участках между контактами ¿, ¿+1 и у, у-1 соответственно.

Для получения компактных и удобных для реализации на ЭВМ соотношений примем, что для траекторий, идущих в ячейке р, ц сетки соответственно вверх или вниз (вправо или влево), вертикальная (горизонтальная) координата точки траектории А определяется относи-

тельно нижележащего или вышележащего (лежащего правее или левее) узла сетки с той же горизонтальной (вертикальной) координатой (рис.2). Положение последующей точки А" в масштабе сеточной ячейки через координаты точки А определяется соотношениями:

A© = (A©, A©) = [[Ak ± (1 -а)(1 -ß)h , A ± (1 -а)(1 -ß)h] , AUh = AUV = 0 [Ah , A ± (1 -a)h] , AUh = 0 , AUV * 0 [Ah ± (1 -ß)h , Av] , AUh * 0 , AUv = 0 [Ah ± (1 - a)ha£, , Av ± (1 -a)h] , a^ < 1- ß , < 1/(1- а)

(19)

[А ± (1 -Р)Л , А ± (1 -р)Л /(о^)] , о$> 1-р , о$> 1/(1-а) ,

где Аин/Аи — отношение горизонтального и вертикального приращений потенциала в ячейке р, q сетки, через которую проходит траектория; о - параметр, определяемый через проводимости соседних с р, q ячеек сетки.

В приведенных формулах знак "+" для А'ь (А'у) ставится при движении траектории в сторону роста индекса р сеточных узлов. Параметр о для траекторий с растущим индексом р и падающим q (рТ, qХ), с (рХ^Х), (рХ^Т) и (рТ^Т) определяется (рис. 2) соответственно так:

a = <

(ap,q +ap +1,q )/(a

p, q + a p, q

(a

+ J/(a +a.

-1) )

p , q p +1, q p , q p , q +1

(ap,q +ap-1,q )/(ap,q +ap,q+1)

(ap, q +a p-1, q )/( a p, q +a p, q-1)'

p,q

p,q-

(20) (21) (22) (23)

Если исходная точка А принадлежит горизонтальному нижнему (верхнему) участку границы, то а=0 в (19) и ор+^= 0 в (20), (21) (ор—1^ = 0 в (22), (23)). Если же точка А принадлежит вертикальному левому (правому) участку границы, то Р= 0 в (19) и ор^—1= 0 в (21), (23) (ор^+1 = 0 в (21), (22)). Аналогичным образом находим направление и ячейку сетки, через которую проходит траектория из точки А " в точку Акоординаты точки А " находим по соотношениям (19) — (23), заменяя в них А на А "и А " на А "). После определения траекторий граничных линий тока из точек А и В контакта - (их конечными точками являются соответственно точки С и Б контакта г), геометрию трубки тока между контактами г и - считаем известной. Отметим, что при т = п для контакта г (г=1, 2, ..., п—1) определяется геометрия (п-г) трубок тока.

Работоспособность методики в приближении однородной среды проверялась на сетках с шагом Н, равным Н/2, Н/4, Н/8; потенциалы в узлах определялись методом теории функций комплексного переменного по формулам (9) — (12). Вычисления показали достаточность шага Н=Н/4-Н/8 для полного выявления геометрии трубок тока.

Определение коэффициента формы

Коэффициент формы и^ определяется следующим образом. Участок трубки тока, описываемый в пределах элиза к, I ломаными линиями (рис. 3), разбивается на N элементарных прямоугольных элементов с неоднородным полем (рис. 4, а) высотой Дхт и поперечными координатами ут,1, ут,г, задаваемыми относительно оси симметрии сечения трубки на нижней, ближайшей к контакту источника, стороне элиза.

Результирующее сопротивление участка трубки тока в элизе к,1 при /--измерении

N

г ,] _ ^

находим суммированием сопротивлении гт всех аппроксимирующих элементов Гк1 — ^ Гт.

т—1

Для получения малой (<0,01 %) погрешности определения г1^ необходимо, чтобы ут /Лхт>

7. Из обобщенной структуры аппроксимирующего элемента (рис. 4, а) с верхним 1 и нижним 2 условными контактами, рассчитываемого методом теории функций комплексного переменного, в зависимости от характера поведения граничных траекторий трубки выявим четыре упрощенных варианта и приведем для них расчетные соотношения.

Н

Рис. 3 Рис. 4

В случае ут,/>ут-1,/ , ут,г>ут-1,г (рис. 3) координаты контакта 2 т-го элемента Ут,иГ ут,и=ут/2, длина прямоугольника ут= ут-11 + ут-1,г, координаты контакта 1 ут,й= ут-11 -уст,

уст= (ут-1,1 ~ут-1,г)/2 , ут,1г= ут-1,г+уст , сопротивление Гт

Гт — Рк ^ (я / 2) /{]п 4 + ар ■- 0,5 ехр [- 2а(1 - (Р + у) / 2)][1 + ехр (- 2ау)]}, (24)

а—пут/(4ахт) , Р=2(ут!1г+ут,и)/ут, у=2(у^,1г-у^,11)/ут.

В случае ут,1<ут-1,1 , ут,г<ут-1,г получаем ут,П=ут,1г=ут/2 , ут=ут,&ут,г , уст= (ут,-ут,Г)/2, ут,и1= ут+1,\ - уст , ут,иг= ут+1>г+ уст, а величина сопротивления гт определяется по (24) с заменой в соотношениях для Р и у индекса I на и.

В случае ут,<ут-1,1 , ут,г>ут-1,г получаем ут,иг=ут, 1Гут/2 , ут= ут, 1+ ут-1,г , уст= (ут,-ут-1,г)/2,

ут,и1=ут+1,1-уст , ут, 1г= ут-1,г+ уст ,

гт — рк г (я / 2) /{1п 4 + ар - 0,5 ехр [- а(8 - 0)][1 - ехр (- 2а0)]}, (25)

8 =2(ут,1г+ ут,ид/ут , 0 =2(ут,1г- ут,и1)/ут .

Наконец, вслучае ут,1>ут-1,1 , ут,г<ут-1,г получаем ут,и1=ут,1г=ут/2 , ут= ут-1,1 + ут,г, уст= (ут-1,1 -ут,г)/2 , ут,иг=ут+1,г +уст , ут,Иг ут-2,1 -уст, а величина гт определяется по (25) с заменой в соотношениях для 8 и 0 индексов I на г.

В случае, когда трубка тока пересекает не противоположные, а соседние стороны элиза, разбиение на элементарные аппроксимирующие прямоугольники ведется от обоих граничных сечений, расположенных под углом я/2, и сопровождается выделением в средней области произвольного дугового участка с углом я/2 и центром в точке О пересечения линий разбиения. Процесс разбиения завершается при достижении точкой О такого положения, при котором для каждого из Р элементарных дуговых аппроксимирующих элементов (рис. 4, б),

Р

получаемых из дугового участка с углами ат<6° ( X ат = % ) и радиусами Я1т, Я2т (т=0,

т=1

1, ..., Р), выполняется неравенство Я2т / Я1т>ехр(5ат). При этом дуговой элемент рис. 4, б

1/2

общего вида сводится конформным отображением £ = (тс/ат)1п[г/( Я1тЯ2т) ] к прямоугольному элементу с параметрами Ахт=1, ут=2/ (Я2т), ут/Ахт > 5, методика расчета которого изложена выше.

По известным величинам Г,/ ,рк,1 и Аик/ (разность потенциалов на границах участка

трубки тока в элизе к, I) находим коэффициент формы «к,/ = г^,/ / рк ^ и величину тока

трубки / = Аик^ / г(погрешность в определении потенциала первого порядка дает погрешность второго порядка при определении величины коэффициента формы). Сравнение найденных и измеренных токов с учетом числа верных знаков в последних позволяет своевременно завершить вычислительный процесс, сохранить корректность постановки задачи [6]. Зная величину 1г- и разности потенциалов на границах участков трубки в других элизах,

находим оставшиеся коэффициенты «к,; (коэффициент формы эквивалентной трубки в предположении однородной среды с погрешностью менее 1% совпадает с коэффициентом реальной трубки по (13)). Отметим, что использование аппроксимирующих элементов с однородным полем ( условные контакты 1 и 2 по всей длине прямоугольника рис. 4, а и дугового элемента рис. 4, б) не позволит решить поставленную задачу, так как суммарное накопление погрешности с ростом N опережает ее уменьшение за счет сближения истинной и аппроксимирующей фигур [10].

Алгоритм реконструкции изображения

В основе предложенного алгоритма реконструкции лежит соотношение (5), вырождающееся в случае узких (Н> 205) трубок тока в соотношение (6). Задача реконструкции после формирования матрицы Я измеренных значений, исходной матрицы коэффициентов формы п ^ и исходной проекционной матрицы р(0) решается по соотношению [1, 2]:

Я = р- п + е . (28)

Работа алгоритма основана на определении величины погрешности е по (29) для некоторого значения Яг- и внесении в проекционную матрицу р поправок, пропорциональных полученному значению погрешности е и коэффициенту формы п в элизах так, чтобы для новых значений р погрешность е стала равна 0. Новое да+1 значение вектора р в элизах рассчитывается из предыдущего т-го значения по формуле

РкГ' = рк? + И" - «У (29)

гдек, I — координаты элизов (к=1, 2, ..., М; I =1, 2, ..., Н']— коэффициент поправки для /--го измерения; Коэффициент поправкиНг,- считается по формуле

М N М N

Ии/ = (Яи-ХХРт - «к/)/XX «й - «к/ . (30)

к=1 I к=1 I

Последовательность решения задачи

Окончательно задача (1)—(4) решается в следующей последовательности. 1. По заданной точности восстановления зависимости р(х,у) определяется минимальный размер элиза Н и общее количество элизов MN.

2. По соотношениям (11) — (14) с учетом конструктивно-технологических ограничений на длину ¡к контактов, их взаимного расположения и количества позиций птах проверяется возможность построения измерительной схемы с требуемыми предельной и реальной разрешающими способностями (Н>5, Н>1к) при п(п—1)/2 >MN и т<п. Результаты анализа позволяют выбрать измерительную схему (конкретные значения т и п, п<птах, т<п), окончательно определить размер элиза Н и количество элизов MN.

3. Проводится дискретизация исследуемой области на элизы с использованием регулярной (реже нерегулярной) сетки, нумерация элизов и позиций контактов.

4. Проводится п(п—1)/2 >MN измерений токов 1г- (г=1, 2, ..., п—1; у=/+1, г+2, ..., п) по схеме рис.1 с погрешностью порядка 0,1...0,5% , формируется матрица Я . Практически Яг-можно получить, обработав на ЭВМ результаты п(п—1)/2 независимых измерений сопротивлений исследуемой структуры.

5. По соотношениям (11)—(14) получаем начальные приближения для поля потенциалов всех краевых задач, начальные приближения для проекционной матрицы р(0). Оцениваем степень неоднородности исследуемой структуры.

6. Определяем геометрию эквивалентных трубок тока, вычисляем коэффициенты

формы Ик\ (для узких трубок — длины траекторий ). Формируем исходную матрицу п(0) (й(0)) указанной в [1, 2] формы.

7. Критерием окончания процесса при использовании матрицы й является условие (7); при использование матрицы п — достижение равенства (в пределах погрешности измерений) матриц рассчитываемых и измеренных токов 11——|<<в+у , где V =0,02...0,03 — погрешность вычислений. При выполнении критерия поставленная задача считается решенной, а матрица р — решением. В противном случае переходим к п.8.

8. Осуществляем цикл реконструкции (3—5 итераций) по соотношениям (29)—(30). Уточняем с помощью конечно-разностного метода поля потенциалов всех краевых задач на последовательности сеток с шагами Н1=Н/4, Н2=Н/8, используя предыдущие значения потенциалов в качестве начального приближения. Уточняем геометрию трубок тока, вычисляем

коэффициенты формы «¿,7 (для узких трубок — длины траекторий ). Корректируем содержимое матриц и (й) и переходим к п. 7. Отметим, что окончательное решение (при известной структуре поля) может быть получено с учетом реального неравномерного распределения тока в трубках.

Экспериментальные результаты

Работоспособность методики проверялась на моделях из резистивной бумаги ЭТБ-2М размерами 300x450 мм с 30 равномерно расположенными по периметру прижимными контактами (да=п=30) длиной ¡к = 8 мм. Размеры элиза Н— 30x30 мм, количество элизов —10x15 =150. Дефекты создавались вклеиванием в исходную бумагу участков бумаг с другими р. Погрешность выполнения условия /,■=/, (порядка 3%) при проведении полного цикла из п(п—1)=870 измерений показала достаточное качество контактов модели.

Обработка результатов измерений осуществлялось на ПЭВМ. Минимальный шаг сетки Н2=Н1=Н/8=3,75 мм, максимальное количество узлов 9600. Проверена работоспособность всех частных методик: определения среднего значения р=р0; оценки степени неоднородности; определения поля потенциалов методами теории функций комплексного переменного и конечно-разностным; определения трубок тока и коэффициентов формы.

Заданное распределение р(х,^) в модели (размеры дефектов и расстояния между ними

равны или кратны размерам элиза) восстанавливалось с погрешностью, не превышающей 7%, причем значительная доля этой погрешности связана, по-видимому, с неоднородностью и анизотропией бумаги («3%), качеством контактов модели («3%), а также погрешностью определения (по величине падения напряжения на дополнительных резисторах 6,2 Ом±1%) токов между контактами. В случае, когда размеры дефектов оказывались меньше размеров элиза, значение р в нем восстанавливалось практически как усредненное по площади элиза.

Расчет с использованием матрицы d приводит к росту погрешности до 12-15%.

Исследовалась и структура с измерительной схемой, отличающейся низкой разрешающей способностью. Подтверждена целесообразность ее использования для определения среднего значения р0 пленки.

Выводы

1. Предложена и экспериментально проверена методика определения распределения электрического сопротивления в двумерной токопроводящей области, базирующаяся на методе томографии и учитывающая нелинейность траекторий линий тока и их зависимость от искомого распределения. Методика обобщается на трехмерные задачи, в которых конечно-разностный метод полностью вытесняет метод теории функций комплексного переменного, эффективный только в двумерных или сводимых к ним задачах. С учетом аналогии методика применима и для определения распределения диэлектрической проницаемости и коэффициента теплопроводности в двух - и трехмерных областях.

2. Предложена методика решения систем внутренних обратных краевых задач с уравнением вида (1), краевыми условиями (2),(3) и исходными данными в постановке (4). Методика допускает обобщение на системы с уравнением вида div[(1/p>)gradU - qU=Fи краевым условием (вместо (2), (3)) aU + раи / dn — v (a> 0, P> 0, v — заданные гладкие или кусочно-гладкие функции, a+P > 0), причем определению (при соответствующем увеличении количества и модификации измерений в (4)) подлежат зависимости р(ху) и q(xy) или р(х, y, U) и q(x, y, Ц).Такое обобщение предполагается темой самостоятельной статьи.

Библиографический список

1. Жилинскас, П.Ю. Определение распределения электрического сопротивления методом реконструктивной томографии // ЖТФ. 1994. Т.64. Вып. 2. С. 83—95.

2. Хермен, Г. Восстановление изображений по проекциям. Основы реконструктивной томографии / Г. Хермен. — М.: Мир, 1983. — 352 с.

3. Reineke, N. Visualisation of transient two-phase flows using multielectrode capacitance sensors for electrical tomography / N. Reineke, D. Mewes // Two-Phase Modeling and Experim. 1995. Р. 739—750.

4. Strizzolo, C.N. Capacitance sensors for measurement of phase volume fraction in two-phase pipelines // IEEE Trahs. oflnstrum. and Meas. 1993. 42. T. 3. Р. 726—729.

5. Сушко, И. А. Визуализация распределений поверхностной проводимости томографического сечения методом зон проводимости // Радиоэлектроника. Изв. вузов. 2013. № 7. С. 60—67.

6. Мацевитый, Ю.М. Идентификация теплофизических свойств твердых тел / Ю.М. Мацеви-тый, С.Ф. Лушпенко. — Киев: Наукова думка, 1990. — 216 с.

7. Коздоба, Л.А. Методы решения обратных задач теплопроводности / Л.А. Коздоба. — Киев: Наукова думка, 1982. — 352 с.

8. Демирчян, К.С. Машинные расчеты электромагнитных полей / К.С. Демирчян, В.Л. Чечурин. — М.: Высш. шк., 1986. — 240 с.

9. Мусхелишвилли, Н.И. Сингулярные интегральные уравнения / Н.И. Мусхелишвилли. — М.: Наука, 1968. — 511 с.

10.Лаврик, В.И. Конформные отображения физико-топологических моделей / В.И. Лаврик, В.П. Фильчакова, А.А. Яшин. - Киев: Наукова думка, 1990. - 376 с.

11.Крылов, В.И. Приближенное вычисление интегралов / В.И. Крылов. - М.: Наука, 1967. - 500 с.

12.Ямпурин, Н.П. Современные проблемы радиоэлектроники с позиций теории конформных отображений / Н.П. Ямпурин, Л.В. Широков, В.Д. Садков. - Арзамас: А1ПИ, 2014. - 209 с.

Дата поступления в редакцию 02.07.2015

V.D. Sadkov, A.V. Lopatkin

DETERMINING THE DISTRIBUTION OF ELECTRICAL RESISTANCE IN THE CONDUCTIVE FILM BY A METHOD OF TOMOGRAPHY

Nizhny Novgorod state technical university n.a. R.E. Alexeev

Technique of determining the distribution of electrical resistance in the conductive film, taking into account the nonlinearity of the trajectories of the streamlines and their dependence on the desired distribution. The method is based on the combined use of the finite-difference method, defining and clarifying the trajectory of streamlines, and reconstruction algorithm that detects the distribution of the electrical resistance along these lines. The method is suitable for solving problems on the distribution of permittivity and conductivity in multidimensional domains.

Key words: impedance tomography, conductive film, distribution of electrical resistance.