CC BY
112
2014 Математика и механика № 3(29)
УДК 519.6
Е.С. Шерина, А.В. Старченко
РАЗНОСТНЫЕ СХЕМЫ НА ОСНОВЕ МЕТОДА КОНЕЧНЫХ ОБЪЁМОВ ДЛЯ ЗАДАЧИ ЭЛЕКТРОИМПЕДАНСНОЙ ТОМОГРАФИИ1
Получены разностные схемы для решения задачи электроимпедансной томографии. Вывод схем осуществляется с помощью метода конечных объёмов на неструктурированных сетках. Численное сравнение для конечных объёмов разной формы выполнено на двух тестовых задачах с аналитическим решением. Полученные результаты также сравниваются с решением по методу конечных элементов.
Ключевые слова: электроимпедансная томография, метод конечных объёмов, метод конечных элементов, разностные схемы, неструктурированные сетки.
1. Введение
Электроимпедансная томография (ЭИТ) - это быстро развивающийся неинва-зивный метод визуализации со спектром применения, включающим медицинскую томографию, неразрушающий контроль материалов и контроль промышленных процессов [1-5]. Метод ЭИТ восстанавливает неизвестное распределение или оценивает несколько параметров электрической проводимости (или полного импеданса) внутри объекта, если измерено электрическое напряжение или сила тока на его границе. По реконструкции даётся оценка внутренней структуры объекта. ЭИТ позволяет наблюдать динамические процессы и работать со статическими изображениями.
Отмечают несколько преимуществ ЭИТ перед другими видами томографии: компактный размер оборудования, простота устройства, невысокая стоимость и быстрый сбор данных. Кроме того, в медицинских приложениях ценится безопасность для биологических тканей. Направление ЭИТ представляется перспективным в разных приложениях, но его широкое внедрение ограничено невысоким качеством получаемых изображений. Наравне с необходимостью усовершенствования аппаратных средств томографии следует уделять внимание разработке математической теории ЭИТ и новых численных методов решения такой задачи.
Восстановление проводимости с помощью ЭИТ является нелинейной и некорректной обратной коэффициентной задачей. Сложность задачи вызвана нарушением условия устойчивости в классическом определении корректности по Адама-ру [2]. Как следствие, реконструкция статических изображений ЭИТ очень чувствительна к ошибкам измерений прибора и погрешности аппроксимации дифференциальной задачи разностной, а также качеству физического моделирования рассматриваемого объекта. Данное исследование направлено на изучение особенностей численного моделирования прямой задачи ЭИТ, являющейся важной компонентой обратной задачи. Отдельное внимание уделено вопросу уменьшения погрешности аппроксимации, которая появляется в ходе дискретизации задачи. Все
1 Работа выполнена при финансовой поддержке Минобрнауки РФ в рамках государственного задания № 2014/223 (код проекта 2382).
результаты получены для двумерной статической задачи и при необходимости могут быть перенесены на трёхмерный случай или использоваться в динамической реконструкции.
В ЭИТ на поверхность объекта прикрепляется набор электродов (рис. 1, а), затем выполняется серия измерений, в которых через объект пропускается электрический ток небольшой силы и высокой частоты. Обычно, сила тока выбирается в интервале от 1 до 5 мА с частотой от 10 до 150 кГц. Схема проведения эксперимента - схема подведения тока и сбора данных - определяется исследователем. Между электродами граница объекта контактирует с воздухом. Результирующее напряжение измеряется на электродах. Вся серия данных «электрический ток -напряжение» используется для реконструкции неизвестного распределения электрической проводимости внутри объекта.
В данной статье использована модель ЭИТ с идеально проводящими электродами, которая называется в научной литературе gap model [6].
Предполагается, что учитываемые моделью физические процессы описываются электрическим и магнитным полями, которые удовлетворяют стационарным уравнениям Максвелла. В данной работе рассматривается случай применения ЭИТ к объектам биологической природы. Ткани проводят электричество, потому что они содержат ионы, которые переносят электрический заряд. Проводимость биологических тканей обладает определенной спецификой и зависит от частоты тока. Когда электромагнитные волны проходят через биологический объект, интенсивность магнитного поля становится пренебрежимо малой [1, 2], поэтому магнитная компонента в уравнениях Максвелла не рассматривается. В этом случае система уравнений Максвелла может быть с использованием закона Ома сведена к уравнению эллиптического типа в частных производных. В силу отсутствия внутри объекта источников электрического поля, а также с учётом 1-го правила Кирхгофа для электрической цепи, сумма втекающих и вытекающих токов равняется нулю - выполняется закон сохранения заряда. Также предполагается, что на электроде, где инжектируется ток, известно значение потенциала электрического тока.
С учётом принятых условий математическая модель ЭИТ для двумерного случая (рис. 1, а) может быть представлена краевой задачей [2, 7]:
2. Постановка задачи ЭИТ
2.1. Физическая постановка задачи
2.2. Математическая постановка задачи
V-(cVcp) = 0, (x,y)eQ;
(1)
ст— = Jk, (xy) e ek, k = l,..^L ; dn
dp
(2)
ст—= 0, (x,y)e dQ/ U ek ;
dn k=1
(3)
где о(ху) - электрическая проводимость, ф(х,у) - потенциал электрического поля, п - единичный вектор внешней нормали, Jk - плотность электрического тока на
k-м электроде, ek - k-й электрод, L - число электродов. Функция электрического потенциала фе H1 (Q), плотность тока j =ст • Уф е L2 (Q).
По принятой в литературе терминологии задача (1) - (3) с функцией электрического потенциала ф(х,у) в качестве неизвестного будет в дальнейшем называться прямой задачей ЭИТ. Обратной задачей ЭИТ именуется задача (1) - (3) по поиску неизвестной функции проводимости o(x,y) в области ß, дополненная данными измерений электрического напряжения на электродах.
3. Построение численных схем для задачи ЭИТ
Чувствительность реконструкции ЭИТ к ошибкам измерений и погрешности аппроксимации дифференциальных операторов ограничивает круг подходящих численных методов, которыми можно с хорошим качеством решать обратную задачу. Как правило, дифференциальные задачи приводят к дискретному аналогу каким-либо сеточным методом. Среди них можно выделить следующие группы: конечно-разностные (МКР), конечно-элементные (МКЭ), гранично-элементные (МГЭ), конечно-объёмные (МКО) и другие. Задачу ЭИТ изначально ассоциировали с электродинамикой и моделированием цепей, поэтому был опробован, а впоследствии и наиболее распространился подход по МКЭ. МКЭ обладает гибкостью, когда необходимо быстро улучшить качество аппроксимации без сложных изменений алгоритма. Но в [8] отмечается, что МКЭ не удовлетворяет условию непрерывности электрического тока, тем самым не соблюдается закон сохранения на элементе сетки и области в целом. С МКЭ ищется решение не конкретной дифференциальной задачи, а обобщённое решение вариационной задачи, выведенной из исходной и записанной в интегральной форме. В данной работе приводятся некоторые результаты построения конечно-объёмных схем для прямой задачи ЭИТ, которые сравниваются с решением МКЭ по нескольким показателям.
3.1. Дискретизация расчётной области
На начальном этапе решения задачи выполняется дискретизация расчётной области rah, и определяется набор сеточных функций электрического потенциала uh и проводимости ch Создание расчётных сеток - отдельная область исследования. В настоящее время доступны свободно распространяемые и коммерческие программные продукты для генерации сеток (GMSH, GAMBIT, QMG, Netgen, Tgrid, FEMLAB и др.).
В данном исследовании численные эксперименты проводились на треугольных неструктурированных сетках (M ячеек, N узлов), построенных в программах GAMBIT и GMSH. Пример расчётной сетки приведён на рис. 1, б. Заметим, что сетка детализована вдоль границы области, потому что градиент потенциала быстро изменяется вблизи контактов электродов. Для оценки качества построенных сеток с треугольными ячейками обычно вводится какой-нибудь числовой показатель. Чтобы осуществлять контроль качества построенных сеток была выбрана группа параметров: коэффициент деформации элемента сетки, угловой коэффициент, соотношение сторон треугольного элемента, соотношение радиуса вписанной к радиусу описанной окружности для элемента сетки. Как правило, точность расчётов снижается, если элементы треугольной сетки сильно отличаются от правильных треугольников. Следует заметить, что решающим фактором при выборе сетки будет не только набор показателей качества и размер сетки, но и её способ-
ность обеспечивать работоспособность конкретного применяемого численного метода в целом. Должна быть обеспечена вычислительная экономичность расчётов и достаточная точность приближенного решения полученных сеточных уравнений, поэтому, по-видимому, следует говорить о том, насколько хорошо сетка подходит для применяемого численного метода.
«5
----
( V
\ Ю У
У&15
1э\/ V V у\
Рис. 1. Двумерная модель области реконструкции О (а), на границе сЮ обозначены электроды ек, к = 1,...,16; расчётная сетка со сгущением к электродам (б), размер сетки -М = 4168 элементов, N = 2181 узел; расчётная сетка для тестовой задачи с точечными электродами 1 и 13 (в), размер сетки - М = 136 элементов, N = 81 узел [8]
3.2. Варианты конечных объёмов
Когда дифференциальная задача с помощью сеточного метода преобразуется в разностную, то выбирается способ соответствия, по которому неизвестные значения сеточной функции связываются с расчётной сеткой. Функции электрического потенциала ф(х,х) и проводимости а(х,,у) заменяются дискретными аналогами, т.е. сеточными функциями ик и <зк соответственно, иъ е иъ,аъ е!ъ, Ц и Т - гильбертовы пространства. На триангуляции области обычно рассматриваются два варианта - значение сеточных функций определяется в некоторой точке внутри ячейки сетки (рис. 2, а) или в узле сетки (рис. 2, б, в). Конечный объём в обоих случаях строится вокруг точки, с которой связываются значения сеточной функции. В настоящей работе обсуждается реализация МКО на трёх конечных объёмах: треугольный конечный объём (рис. 2, а), барицентрический конечный объём (рис. 2, б) и конечный объём - ячейка Дирихле - Вороного (рис. 2, в).
Р2
Рис. 2. Фрагменты расчётной сетки: а - треугольный конечный объём с центром в точке Р0; б - барицентрический конечный объём, построенный вокруг вершины Р0; в - конечный объём - ячейка Дирихле - Вороного, - построенный вокруг вершины Р0
7
10
0
5
Р
Р
Барицентрическая ячейка строится вокруг вершины Р0 по точкам пересечения медиан и точкам середин сторон треугольников. В ячейке Дирихле - Вороного вместо барицентров выбираются точки пересечения срединных перпендикуляров. Треугольные конечные объёмы соответствуют ячейкам сетки, в них положение точки Р0 задается в центре масс, точке пересечения срединных перпендикуляров или биссектрис.
3.3.Вывод численных схем
Согласно принципу МКО для каждого конечного объёма записывается ин-тегро-интерполяционный баланс. Уравнение (1) интегрируется по конечному объёму 0.т,
|У-(стУф) " = 0. (4)
"т
По формуле Грина уравнение (4) преобразуется к равенству нулю криволинейного интеграла по границе от плотности тока
стдф <Ц = 0. (5)
д0 дп
д"т
Дальнейшая оценка интеграла в левой части (5) зависит от формы конечного объёма и способа аппроксимации градиента потенциала в подынтегральном выражении. Функция проводимости задается кусочно-постоянной по треугольным элементам сетки, 0 < с < о(х,у) < С, с и С - константы.
На треугольных конечных объёмах (рис. 2, а) интеграл в правой части (5) вычисляется по сторонам текущего треугольника 0.т по интерполяционной формуле, в данном случае - средних прямоугольников. Опробовано несколько вариантов расположения точки Р0 внутри конечного объёма, с которой ассоциировались сеточные функции проводимости и электрического потенциала ик. Дискретные значения определялись в точке пересечения медиан, срединных перпендикуляров или биссектрис треугольника. Производные аппроксимировались в серединах сторон конечного объёма, причём оценка выполнялась с учётом условий сопряжения на границе двух смежных ячеек, например для грани ' (рис. 2, а):
стдфТ° г дф"
дп )к V дп
К =(ф)^
чтобы схема стала чувствительной к изменению свойств токопроводящей среды. Потоки через границу конечного объёма могут быть аппроксимированы разными соотношениями. Один из предложенных подходов рассматривает аппроксимацию
(дфУ0 ик - иР (дфУ* иР - ия вида I — I « —1-- и I — I и —1-L, где ЛпР(1 - длина отрезка перпенди-
V дп ДпР(| V дп ДпР1
куляра из точки Р0 до стороны АпР1 - от точки Р1 до стороны а иР0 = ф(хР0,уР0). Преобразования с использованием условий сопряжения приводят к системе сеточных уравнений:
АР1 (иР1 - иР0 ) Д/у + АР2 (иР2 - иР0 ) Д' + АР3 (иР3 - иР0 ) Д1кг = 0 (6)
где
а Р а Р
Ро Рч
Ам =
Ч аР ДпР +аР ДпР
Ч = 1,2,3,
О Ч Ч о
для текущего конечного объёма; Л/у- обозначена длина стороны г/. Заметим, что при таком выборе конечного объёма граничные условия Неймана (2) - (3) записываются без погрешности аппроксимации.
На барицентрических конечных объёмах [9] (рис. 2, б) функция потенциала
приближается на множестве базисных функций {Тп (х,у)}}= , где N - количество
узлов сетки, полиномом
'(х,У) = ЁипТп (х,у),
(7)
в котором ип - значение сеточной функции и в п-й вершине сетки.
Каждая базисная функция, как и в МКЭ, задается локально на треугольной ячейке сетки по узловым значениям, причем её значение не равно нулю только в одной вершине треугольника. Рассмотрим треугольник АРоРтРт+1, т = 1, для каждой вершины запишется своя базисная функция:
1
'(х, У ) = - 1
Т (т) Ро
(х,У) =
1
1
У Ур
1 т
Урт+
Т(т) I
2S„
ЛР о
х
УРо У
Урт+
Т
(т) Р
т+1
(х У) =
2Sm
1 хРо УРо _ 1 = 2 1 хРо УРо
1 Рт УРт , Sm 1 Рт УРт
1 х У 1 хРт+1 УРт+1
(8)
где Sm - площадь треугольника ДРоРтРт+1, т = 1, 2,...МР, МР - число треугольников вокруг точки Ро. Для одного треугольника базисные функции в сумме дают 1. Тогда распределение потенциала на каждом треугольнике интерполируется линейной функцией, принимающей значения потенциала в его вершинах. Для ДРоРтРт+1 в выражении (7) останутся только 3 ненулевые базисные функции,
ик (х,у) = ир Тт (х,у) + и
Р Тр) (х,у) + ир
т
Т(рт) (х,у). Кроме того, градиент по-
тенциала принимает постоянное значение на АРоРтРт+1. Если подставить градиент в интеграл по контуру (5) и проинтегрировать по участкам границы многоугольного объёма, которые проходят в ДРоРтРт+1, то будет получен вклад в разностную формулу от этого треугольного элемента с постоянной проводимостью ст. Схема для конечного объёма вокруг узла Ро (рис. 2, б) запишется в виде
МР а
X ТтК ((УРт - УРт+1)2 + (хР
т=1 т
- х-,
)2)+
)(хРт+1 - хРт ))
+ (9)
+иРт ((УРт+1 - УРо)(УРт - УРт+1)+(хРо - х +иРт+1 ((УРо - УРт )(УРт ~ УРт +1 ) + (^т - хРо )(^т+1 - ^т )И = о,
где суммирование выполняется по всем треугольным элементам сетки с общей вершиной Ро, причем, когда значение индекса т в (9) больше Мр, то т = 1.
п=1
На объёмах Дирихле - Вороного (рис. 2, в) используется другой способ аппроксимации производной в (5). Функция потенциала задается линейной функцией ик(х, у) = ат(х-хР0) + Ьт(у-уро)+иро по т-му треугольному элементу сетки:
ат = ((иРт—иР0)(УРт+\—УР0) — (uPm+1—uP0)(УPm—УP0))/(2Sm), Ьт = ((иРт+1-иР0)(хРт-хР0) — (uPm—uP0)(xPm+1—xP0))/(2Sm), &т = ((хРт хР0)(уРт+1 -уР0) — (хРт+1-хР0)(уРт—уР0))/2.
Градиент линейного потенциала принимает постоянное значение на АР0РтРт+1. Тогда подстановка градиента и значения вектора внешней нормали к текущему отрезку границы конечного объёма Дирихле в (5) даёт вклад в коэффициенты для вершины Р0 от АР0РтРт+1. Схема для конечного объёма (рис. 2, в) вокруг узла Р0 имеет вид
мР
Хст»
т=1
- иР
Р0
АКтСт
«/(хР - хР )2 + (уР - уР )2
. V Рт Р0 ' у Рт •ЛР0/ (10)
(
- иР<>
\Ст^-т+\ I
:Рт+1 - хР0)2 +(урт+1 - ур0)2
= 0,
где через Ст обозначена точка пересечения срединных перпендикуляров в АР0РтРт+1; Ят, Ят+1 - середины сторон Р0Рт и Р0Рт+1 соответственно. Под |-| подразумевается длина отрезка.
Интегральные соотношения (5) для граничных узлов на барицентрических объёмах и объёмах Дирихле - Вороного аппроксимируются с использованием значений производных от функции ик(х,у), взятых из краевых условий (2), (3). Каждое соотношение вносит вклад только в соответствующую компоненту правой части СЛАУ. В схеме для треугольных конечных объёмов краевые условия (2), (3) подставляются в криволинейный интеграл для граничных ячеек.
Поскольку в сеточных построителях нельзя гарантировать «правильность» треугольных ячеек, МКО на ячейках Дирихле - Вороного может давать неоднозначные результаты. Если точка пересечения срединных перпендикуляров одного из треугольников не лежит внутри, то в коэффициентах для вершин накапливается погрешность.
Решение задачи Неймана для эллиптического уравнения (1) - (3) определяется с точностью до некоторой аддитивной постоянной [10], при практической реализации у решаемых СЛАУ отсутствует диагональное преобладание. Однако в измерительной системе ЭИТ всегда есть электрод заземления, который также считается точкой отсчёта распределения электрического потенциала в объекте. Этот факт можно использовать для вычислительных целей и задавать на одном электроде условие Дирихле.
3.4. Исследование свойств численных схем
Запись соотношения (6), (9) или (10) на всей расчётной сетке приводит к разреженной симметричной системе линейных алгебраических уравнений (СЛАУ) с диагональным преобладанием относительно неизвестных сеточных значений потенциала при использовании условия Дирихле на одном из электродов. Числа обусловленности систем зависят от размера сетки, порядка нумерации узлов, свойств треугольных элементов сетки, распределения проводимости и т.п. На тес-
товой сетке (рис. 1, в) число обусловленности равно 133 для конечно-объёмных схем (9) и (10) и конечно-элементной системы, 1596 - для конечно-объёмной схемы (6) на однородной задаче (o(x,y) = 1 См/м в Q); 1428 и 7647 - на неоднородной задаче (o(x,y) = 5 См/м в круге радиуса r0 = 5 см и o(x,y) = 1 См/м вне круга) соответственно. На тестовой сетке (рис. 1, б) получено число обусловленности 2,4104 и 6,1104 на однородной задаче, 3,9104 и 2,6105 - на неоднородной задаче соответственно.
Построенные конечно-объёмные схемы локально консервативны, так как на любой грани, разделяющей соседние конечные объёмы, выходящий поток равен входящему потоку. Локальная консервативность построенных разностных схем позволяет говорить и об их интегральной консервативности: выполнении интегрального закона сохранения тока для всей области Именно выполнение этого условия в разностном виде обеспечивает разрешимость полученных систем сеточных уравнений [11]. Устойчивость разностных схем (6), (9) и (10) доказывается с использованием мажоранты Гершгорина и принципа максимума [12].
3.5. Выбор численного метода решения СЛАУ
Для решения полученной СЛАУ были опробованы различные итерационные методы: Якоби, Гаусса - Зейделя, верхней релаксации (SOR), методы подпространств Крылова BiCG и BiCGSTAB. В качестве критерия завершения итерационного процесса использовались условия малости невязки и ошибки (<е = 10-10). Вычислительную устойчивость и высокую скорость сходимости на МКО- и МКЭ-системах сеточных уравнений продемонстрировал SOR со значением параметра релаксации ю = 1,9. Этот метод выбран для дальнейших вычислений. Заметим, что для МКО- и МКЭ-схем нет необходимости хранить матрицы коэффициентов целиком. Методы решения СЛАУ были модифицированы для работы со сжатыми данными, в которых записаны только ненулевые компоненты матриц.
4. Тестирование на задаче с аналитическим решением
Проведена серия экспериментов на задаче с аналитическим решением, в которых оценивалась точность полученных разностных методов. Рассматривалась двумерная область О радиуса Я с круговым включением радиуса г0 в центре (рис. 1, а). Аналитическое решение для модели ЭИТ (1) - (3) с неоднородностью в центре круга записывается в полярных координатах через ряд Фурье [13]:
u (r ,9) =
W 1
R 2 -Cn
n=1 n
í
W 1
R 2 -
n=1 n
An I —
R i
R
(an cos(n9) + bn sin(n9)), 0 < r < r0
-Bj -
,R
(11)
(( cos(n9) + bn sin(n9)), r0 < r < R;
Cn =-
1
An =-(1 + *)Cn, Bn = An -1;
s2n ' n 2
1+*-(1 - *) R
1 2n 1 2n
an =- f j(9)cos(n9)d9, bn =- f j(9)sin(n9)d9,
П 0 П 0
где j(0) - распределение плотности тока по поверхности круга, * = Oi/a2, о1 - элек-
трическая проводимость включения, с2 - фоновая электрическая проводимость круга.
Относительная погрешность вычислений определена через аналитическое решение по формуле
\\uh - uh
11 num an|1 , x 100 %,
h h где u num - элекгрическии потенциал, найденный численным методом, u an - аналитическое решение (11), uan_max, uan_min - максимальное и минимальное значение аналитического потенциала в области Q соответственно.
Тест 1. Первая группа тестов выполнена на сетке (рис. 1, в), представленной в статье G. Dong и др. [8], где авторы провели сравнение МКЭ- с МКО-подходом, разработанным на барицентрических конечных объёмах. Авторы не упоминают использование базисных функций (8), но говорят о линейной интерполяции потенциалов по треугольным элементам сетки. Вывод формул конечно-объёмной схемы и детали аналитического решения для тестов в [8], в основном, упущены. Отдельным пунктом выделено различие обработки граничных условий в МКО- и МКЭ-подходах, когда сила тока задается на точечных электродах. Численные эксперименты выполнены на однородном круге с постоянной проводимостью. Относительная погрешность находится в пределах 0,3 % для МКО и 0,5 % для МКЭ в большинстве граничных и внутренних узлов. Точки 1 и 13, где задавались граничные условия, были исключены из сравнения со ссылкой на сингулярность аналитического решения в этих точках и бесконечное значение потенциала в них. Максимальная относительная ошибка вблизи этих точек возрастает до 2,56 % для МКО и 7,4 % для МКЭ [8].
В настоящей работе приведены результаты сравнения МКО с тремя вариантами выбора конечного объёма (рис. 2) и МКЭ на однородном объекте. Тесты выполнены для круга радиуса R = 13 см. Для однородного случая электрическая проводимость задавалась с = 1 См/м. Электрический ток инжектируется на паре узлов 1-13 (рис. 1, в), плотность тока J = 1 А/м2, J3 = -1 А/м2. Вычисленная функция электрического потенциала сравнивается с аналитическим решением (11) в граничных узлах 1-24 и промежуточных узлах вдоль оси Ох между точками 1 и 13.
Конечно-объёмные схемы (9) и (10) дают схожее приближенное решение, которое проиллюстрировано графиками (рис. 3). На графиках приводится аналитическое и численное решение в граничных узлах (рис. 3, а), а также вдоль сечения области осью Ох (рис. 3, б). Относительная погрешность методов находится в пределах 0,25 % в граничных точках и 1 % - во внутренних точках вдали от электродов. Аппроксимация граничных условий приводит к погрешности до 2,75 % в узлах, относящихся к электродам. Подобные результаты объясняются тем, что распределение потенциала на каждом элементе сетки использует линейную интерполяцию потенциала в узлах. Матрицы коэффициентов схожи на одинаковой сетке, СЛАУ отличаются правой частью, поскольку граничные условия обрабатываются по-разному.
Также были проведены расчёты МКЭ для однородного круга. Относительная погрешность метода сохраняется на том же уровне. Вариационная формулировка МКЭ для задачи Неймана (1) - (3) получена по [14, 15] на конечномерном подпространстве Uh е H1 (Q), где приближенное решение uh может быть представ-
лено в виде линейной комбинации кусочно-линейных базисных функций (7). Решая составленные сеточные уравнение для метода Галеркина относительно иь, находим приближение электрического потенциала в узлах расчётной сетки.
Приближённое решение по МКО на треугольных конечных объёмах получено в центрах масс треугольных ячеек. Для сравнения с аналитическим решением на границе области д^ численное решение интерполировалось из центров треугольников в граничные точки по направлению перпендикуляра с использованием разложения в ряд Тейлора. На графиках (рис. 3) метод по схеме (6) демонстрирует относительную погрешность до 0,75 % в граничных точках, где нет контакта с электродами, и 5,06 %о - на электродах. Электрический потенциал вдоль сечения 1-13 интерполируется по значениям, которые вычислены в паре смежных объёмов, имеющих общую сторону на оси Ох. Относительная погрешность вдоль сечения составляет 0,17 %о во внутренних узлах, а максимальное значение ошибки 1,49 % находится вблизи электродов 1 и 13.
6т
Угол, рад
6
4-
-Аналитическое решение
• • • МКО, барицентрические ячейки ооо МКО, треугольные ячейки
-15 -10 -5 0 5 10 15
х9 см
ей И
ю 1
о
л ¡2 2-о
о
И н О
°1 О
—9 1 • 1 • 1 '
3 4 5
Угол, рад
ей И
ю
§ 4 4-
л ¡2 2-
8 о
о
И н О
-15 -10 -5 0 5 10 15
х, см
Рис. 3. Сравнение приближенного МКО-решения с аналитическим на круге для случая однородного объекта с точечными электродами: а, б - распределение приближенного и аналитического решения по внешней границе круга и по сечению у = 0; в, г - относительная погрешность вычислений
6
6
в
г
0
6
7
Тест 2. В другой серии экспериментов конечно-объёмные разностные схемы (6), (9) и (10) протестированы на модели с электродами размером 2 см (рис. 1, а) для случая однородного и неоднородного объекта. Тесты выполнены для круга радиуса Я = 13 см с внутренним включением размером г0 = 5 см. Когда рассматривался однородный случай, электрическая проводимость задавалась как в тесте 1. В экспериментах с неоднородным объектом проводимость включения и фона задавалась а1 = 5 См/м и а2 = 1 См/м соответственно. Расчёты выполнены на сетке 4186 элемент (рис. 1, б). Плотность тока задана на паре диаметрально размещённых электродов в\ и е9, = 1 А/м2, З13 = -1 А/м2.
На графиках (рис. 4) аналитическое решение задачи сравнивается с численным МКО-решением для схемы на барицентрических объёмах. Приближенное решение по схемам (9) и (10) незначительно отличается друг от друга и от МКЭ решения в граничных узлах, расположенных вдали от активной пары электродов, что может быть проиллюстрировано графиком (рис. 4, а). Относительная погрешность в граничных узлах находится на уровне 0,01 % с отклонением до 0,34 % в узлах, относящихся к электродам. Вдоль сечения объекта осью Ох аналитическое решение сравнивается с вычисленной функцией потенциала, которая интерполируется в 100 промежуточных точках (рис. 4, б). Относительная ошибка достигает 0,8 % внутри объекта и увеличивается до 4,5 % на границе, что можно объяснить неточностью интерполяции.
Результаты для конечно-объёмной разностной схемы (6) показали, что в граничных точках относительная погрешность возрастает до 0,9 % в местах контакта активных электродов и держится на уровне 0,2-0,4 % на участках, граничащих с воздухом. На поперечном сечении круга относительная ошибка - до 1,5 % во внутренних точках и 5 % вблизи электродов.
Угол, рад х, см
Рис. 4. Сравнение приближенного МКО-решения с аналитическим решением на круге для случая однородного и неоднородного объекта при пропускании тока через электроды е1 - е9: а - распределение потенциала по внешней границе круга; б - распределение потенциала вдоль сечения у = 0
При проведении расчётов для прямой задачи ЭИТ с неоднородным кругом видно, что профиль функции потенциала изменяется, когда электрический ток
проходит через неоднородность. Математическая модель (1) - (3) пригодна для изучения изотропных объектов и ограничивается случаем идеального контакта электродов.
Тест 3. Выполнены расчёты на группе измельчающихся сеток одинакового качества на задаче с неоднородной проводимостью. Более детальные сетки получены разделением треугольных элементов на четыре части. Чтобы сохранить качество сетки, в исходном треугольнике соединялись середины сторон, и образовывалось 4 новых треугольных элемента. Результаты для конечно-объёмных схем и метода конечных элементов приведены в таблице и говорят о численной сходимости методов на последовательности сеток с уменьшающимся шагом.
Относительная ошибка конечно-объемных разностных схем и метода конечных элементов
Номер Размер сетки, Относительная ошибка, %
схемы* ячеек По всей области По граничным узлам
Максимум Среднее Максимум Среднее
790 3,202 0,256 3,202 0,300
1 3160 2,382 0,211 2,382 0,292
12640 1,795 0,159 1,795 0,230
790 1,156 0,060 1,156 0,101
2 3160 0,347 0,059 0,347 0,067
12640 0,292 0,064 0,292 0,064
790 1,157 0,059 1,157 0,101
3 3160 0,345 0,058 0,345 0,066
12640 0,290 0,063 0,290 0,063
790 1,156 0,060 1,156 0,101
4 3160 0,347 0,059 0,347 0,067
12640 0,292 0,064 0,292 0,064
* Обозначения: конечно-объёмная схема для треугольных конечных объёмов - 1, для барицентрических конечных объёмов - 2, для конечных объёмов Дирихле - Вороного - 3, метод конечных элементов - 4.
6. Заключение
Для численного решения прямой задачи электроимпедансной томографии на неструктурированных сетках построены конечно-объёмные аппроксимации. В качестве конечных объёмов выбирались барицентрические ячейки, ячейки Дирихле - Вороного, треугольные конечные объёмы. Исследование аппроксимационных свойств полученных разностных схем затруднено вследствие нерегулярной структуры вычислительной сетки. Устойчивость разностных схем доказывалась с использованием принципа максимума и мажоранты Гершгорина. Сходимость проверялась численно на последовательности сеток с кратно уменьшающимися размерами ячеек. На основе сравнительного анализа получаемых численных решений по построенным разностным схемам с точным решением модельной задачи для однородного и неоднородного дисков и с численным решением на основе метода конечных элементов установлено, что высокую точность расчётов (сравнимую с МКЭ) обеспечивают разностные схемы на барицентрических ячейках, ячейках Дирихле - Вороного. Несколько худшие результаты, тем не менее с высокой точностью, показывает разностная схема на треугольных конечных объё-
мах, которая позволяет без погрешности аппроксимации учесть граничные условия Неймана.
Для решения прямой задачи ЭИТ разработан комплекс программ на языке C/C++.
ЛИТЕРАТУРА
1. Пеккер Я.С., Бразовский К.С., Усов В.Ю. и др. Электроимпедансная томография. Томск: Изд-во НТЛ, 2004. 192 с.
2. Lionheart W., Polydorides N., Borsic A. Electrical Impedance Tomography: Methods, History and Applications. Manchester, 2004. 62 p.
3. Brown B.H. Electrical impedance tomography (EIT): a review // J. Med. Eng. Technol. 2003. V. 27. No. 3. P. 97-108. doi: 10.1080/0309190021000059687.
4. Holder D.S. Medical impedance tomography and process impedance tomography: a brief review // Meas. Sci. Technol. 2001. V. 12. P. 991-996.
5. Barber D.C., Brown B.H. Applied Potential Tomography // J. Phys. E: Sci. Instrum. 1984. V. 17. No. 9. P. 723-733. doi:10.1088/0022-3735/17/9/002.
6. Somersalo E., Cheney M., Isaacson D. Existence and uniqueness for electrode models for electric current computed tomography // SIAM J. Appl. Math. 1992. V. 52. No. 4. P. 10231040.
7. Шерина Е.С., Старченко А.В. Численный метод реконструкции распределения электрического импеданса внутри биологических объектов по измерениям тока на границе // Вестник Томского государственного университета. Математика и механика. 2012. № 4(20). С. 36-49.
8. Dong G., Zou J., Bayford R.H., et al. The comparison between FVM and FEM for EIT forward problem // IEEE Trans. Magnetics. 2005. V. 41. No. 5. P. 1468-1471. doi: 10.1109/ TMAG.2005.844558.
9. Li J., Yuan Y. Numerical simulation and analysis of generalized difference method on triangular networks for electrical impedance tomography // App. Math. Modelling. 2009. V. 33. No. 5. P. 2175-2186.
10. ТихоновА.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 735 с.
11. МарчукГ.И. Методы вычислительной математики. М.: Наука, 1989. 608 с.
12. Крылов В.И., Бобков В.В., П Монастырный.И. Вычислительные методы. М.: Наука, 1977. 399 с.
13. Isaacson D. Distinguishability of conductivities by electric current computed tomography // IEEE Trans. Medical Imaging. 1986. V. 5. No. 2. P. 91-95.
14. Vauhkonen M. Electrical impedance tomography and prior information. PhD. Kuopio, 1997. 110 p.
15. Рояк М.Э., Соловейчик Ю.Г., Шурина Э.П. Сеточные методы решения краевых задач математической физики. Новосибирск: Изд-во НГТУ, 1998. 120 с.
Статья поступила 14.04.2014 г.
Sherina E.S., Starchenko A.V. FINITE VOLUME SCHEMES FOR THE ELECTRICAL IMPEDANCE TOMOGRAPHY PROBLEM. In electrical impedance tomography, the currents are applied on electrodes placed on the surface of an object. The electrical conductivity is reconstructed in the interior of the object using voltage measurements on its surface. Knowing the conductivity distribution provides information about the internal object's structure which can be useful for medical and industry applications. Instability of the EIT problem causes difficulties challenging a successful reconstruction. Since a static EIT imaging is sensitive to the measurement noise and approximation errors, this paper addresses the problem of reducing the latter. The finite volume method is presented for solving the EIT forward problem, which is a significant part of the inverse problem. Finite volume schemes were obtained for unstructured grids. The schemes were derived for three kinds of finite volumes, which can be considered relative to triangulation
38
Е.С. WepHHa, A.B. frapneHKO
of the domain. The approaches are based on vertex-centered and cell-centered discretization, where the numerical solution is associated with mesh vertices or mesh elements, respectively. In the first case, a finite volume approximation was introduced on barycentric volumes and Dirichlet - Voronoi volumes on the assumption of a linear distribution of electric potential. Triangular finite volumes were utilized for approximation based on cell-centered discretization. Both cases suggested a piecewise constant conductivity distribution over grid cells. Numerical comparison for the obtained finite volume schemes was carried out on a test problem that can be solved analytically. The results were compared to a solution by the finite element method.
Keywords: electrical impedance tomography, finite volume method, finite element method, difference schemes, unstructured meshes
SherinaEkaterina Sergeevna (M. Sc., Tomsk State University, Tomsk, Russian Federation) E-mail: sherina@math.tsu.ru
Starchenko Aleksandr Vasilyevish (Doctor of Physics and Mathematics, Prof., Tomsk State University, Tomsk, Russian Federation) E-mail: starch@math.tsu.ru
REFERENCES
1. Pekker Ya.S., Brazovskiy K.S.,. Usov V.Yu, Plotnikov M.P., Umanskiy O.S. Elektroim-pedansnaya tomografiya. Tomsk, NTL Publ., 2004. 192 p. (in Russian)
2. Lionheart W., Polydorides N., Borsic A. Electrical Impedance Tomography: Methods, History and Applications. Manchester, 2004. 62 p.
3. Brown B.H. Electrical impedance tomography (EIT): a review (2003) J. Med. Eng. Technol. V. 27. No 3, pp. 97-108. doi: 10.1080/0309190021000059687.
4. Holder D.S. Medical impedance tomography and process impedance tomography: a brief review (2001) Meas. Sci. Technol. V. 12, pp. 991-996.
5. Barber D.C., Brown B.H. Applied Potential Tomography (1984) J. Phys. E: Sci. Instrum. V. 17. No. 9, pp. 723-733. doi:10.1088/0022-3735/17/9/002.
6. Somersalo E., Cheney M., Isaacson D. Existence and uniqueness for electrode models for electric current computed tomography (1992) SIAM J. Appl. Math. V. 52. No. 4, pp. 10231040.
7. Sherina E.S., Starchenko A.V. Chislennyy metod rekonstruktsii raspredeleniya elek-tricheskogo impedansa vnutri biologicheskikh ob"ektov po izmereniyam toka na granitse (2012) Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika. No. 4(20), pp. 36-49. (in Russian)
8. Dong G., Zou J., Bayford R.H., Xinshan M., Shangkai G., Weili Y., Manling G. The comparison between FVM and FEM for EIT forward problem (2005) IEEE Trans. Magnetics. V. 41. No. 5, pp. 1468-1471. doi: 10.1109/TMAG.2005.844558.
9. Li J., Yuan Y. Numerical simulation and analysis of generalized difference method on triangular networks for electrical impedance tomography (2009) App. Math. Modelling. V. 33. No. 5, pp. 2175-2186.
10. Tikhonov A.N., Samarskiy A.A. Uravneniya matematicheskoy fiziki. Moscow, Nauka Publ., 1977. 735 p.
11. Marchuk G.I. Metody vychislitel'noy matematiki. Moscow, Nauka Publ., 1989. 608 p. (in Russian)
12. Krylov V.I., Bobkov V.V., Monastyrnyy P.I. Vychislitel'nye metody. Moscow, Nauka Publ., 1977. 399 p.
13. Isaacson D. Distinguishability of Conductivities by Electric Current Computed Tomography (1986) IEEE Trans. Medical Imaging. V. 5. No. 2, pp. 91-95.
14. Vauhkonen M. Electrical impedance tomography and prior information. PhD. Kuopio, 1997. 110 p.
15. Royak M.E., Soloveychik Yu.G., Shurina E.P. Setochnye metody resheniya kraevykh zadach matematicheskoy fiziki. Novosibirsk: NGTU Publ., 1998. 120 p. (in Russian)
