автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Вычислительные схемы МКЭ-моделирования трехмерных элeктpомагнитных и тепловых полей в сложных областях

с-' *

Новосибирский государственный технический университет ° #

N %

На правах рукописи

Соловейчик Юрий Григорьевич

Вычислительные схемы МКЭ-моделирования трехмерных электр0!\1агнит11ых н тепловых полей в сложных областях

Специальность 05.13.16. - Применение вычислительной техники, математического моделирования и математических методов в научных исследованиях (в области технических наук)

Автореферат диссертации на соискание ученой степени доктора технических наук

Новосибирск - 1997

Работа выполнена в Новосибирском государственном техническом университете

Официальные оппоненты: доктор технических наук, старший научный

сотрудник Б.М. Глинский

доктор физико-математических наук, профессор В.Г1. Ильин

доктор технических наук, профессор А.Д. Рычков

Ведущая организация: Институт теоретической и прикладной механики СО РАН

Защита состоится , 1997 года в 13 часов на заседании

диссертационного совета Д063.34.03 при Новосибирском государственном техническом университете (630092, Новосибирск-92, пр.К.Маркса, 20).

С диссертацией можно ознакомиться в читальном зале библиотеки НГТУ. Автореферат разослан Си^^е^ЖЯ. 1997 г.

Ученый секретарь

диссьртацнонного /

соьита .У О • к.т.л,доц. Г. П.Чимшьшш

^ / /у* V

Общая характеристика работы

Актуальность темы. Метод конечных элементов (МКЭ) является одним из наиболее мощных инструментов математического моделирования. Теория МКЭ и многие аспекты его применения для решения различных задач науки и техники разрабатывались в основном в работах зарубежных ученых: Зенкевича О., ОбэнаЖ.-П., Стренга Г., СьярлеФ., Темама Р. и др. Более широкому использованию МКЭ при решении различных прикладных задач способствовали работы отечественных ученых: Ильина В.П., Корнеева В.Г., Марчука Г.И., Шайдурова В В. и др. Значительный вклад в развитие МКЭ применительно к решению задач электромагнетизма внесли работы Кулона Ж.-Л., Сабоннадье-ра Ж.-К., Сильвестера Р. и др:

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

В проблеме построения конечноэлементных сеток в качестве наиболее важных можно выделить два аспекта:

• правильное описание сложной-трехмерной геометрии;

• построение сеток, оптимальных с точки зрения минимизаций ошибок аппроксимации при ограниченном числе узлоа, в том числе адаптивных сеток.

Этим вопросам уделяется очень много внимания в большом количестве работ, связанных с использованием МКЭ. Столь значительный интерес к проблеме построения сеток объясняется тем, что при решении любой практической задачи методом конечных элементов исследователь тем или иным образом вынужден разбивать расчетную область на большое число отдельных ячеек определенного вида, и ог того, насколько эффективно ему удастся это сделать, часто зависит сама возможность решения с нужной точностью конкретной практической задачи.

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

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

С другой стороны, имеется возможность использования в качестве базисных функций кусочных полиномов высоких порядков. Конечно, при одинаковой дискретизации расчетной области кусочные полиномы более высоких порядков позволяют существенно снизить ошибку конечноэлементной аппроксимации. Однако, это не означает их безусловного преимущества перед кусочными полиномами более низких порядков с точки зрения суммарных затрат вычислительных ресурсов на получение численного решения с заданной точностью. Действительно, при использовании полиномов высоких порядков даже на коаечноэлемектных сетках с крупными ячейками число узлов может быть не меньше числа узлов в конечноэлементной сетке с гораздо более мелкими ячейками при использовании на ней в качестве базисных функций полиномов более низких порядков (из-за необходимости введения дополнительных узлов на конечном элементе при задании на нем полинома высокого порядка).-Не менее важным является и то, что при использовании в качестве базисных функций кусочных полиномов высоких порядков Существенно увеличивается заполненность матриц конечноэлементных систем линейных алгебраических уравнений (СЛАУ) ненулевыми элементами. А поскольку для решения таких СЛАУ используется чаще всего метод сопряженных градиентов (являющийся одним из наиболее эффективных методов решения конечноэлементных СЛАУ, получаемых в результате аппроксимации трехмерных краевых задач) или какой-либо подобный ему метод, то пропорционально росту заполненности матрицы конечноэлементной СЛАУ растут и вычислительные затраты на решение этой СЛАУ как по памяти, так и по времени счета. Кроме того, конечные элементы с базисными функциями в виде кусочных полиномов высоких порядков могут оказаться очень неэффективными в областях со сложными границами (внешними й внутренними), когда из-за необходимости описания геометрии этих границ с нужной точностью потребуются конечные элементы достаточно малых размеров и используемые на них кусочные полиномы высоких порядков окажутся чрезмерно точными для аппроксимации исследуемого решения. То есть, если для получения численного решения с необходимой для прак-тпчсских целей точностью на конечноэлементной сетке с небольшими размерами ячеек, требуемыми для описания сложной геометрии расчетной области, достаточно было бы кусочных нолиномов низких порядков, то более высокая точность конечных элементов с полиномами высоких порядков окажхчея фактически бес-

полезной, и довольно существенный рост вычислительных затрат при их использовании будет совершенно неоправданным.

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

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

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

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

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

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

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

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

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

Научная ношппа работы состоит в следующем.

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

2. Дано обоснование применимости методики снижения погрешности аппроксимации при выделении полей влияния отдельных объектов путем вычисления разности решений двух краевых задач (без решения краевой задачи в специальной постановке).

3. Предложены подходы, основанные на сочетании полуаналптпчеекпх ме-

годов решения задач электромагнитного зондирования ^Земли с методами численного анализа, использующими аппроксимацию искомых решений на сетках. Эти подходы позволяют соединить преимущества полуаналитических методов при расчетах характеристик поля в средах с простой геометрией с возможностями сеточных методов при учете различных трехмерных неоднородностей.

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

5. Предложена конечноэлементная реализация известной схемы против потока при численном моделировании тепловых полей с учетом конвективного геплонереноса.

6. Предложена вычислительная схема решения систем конечноэлементных уравнений с симметричными и несимметричными матрицами. Эта схема наряду с другими наиболее эффективными методами решения СЛАУ с сильно разреженными матрицами была с успехом использована при решении многих конечноэлементных СЛАУ.

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

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

TELMA и были с успехом использованы при решении большого числа сложных трехмерных задач перечисленных выше типов.

Личный иклад. Все результаты, изложенные в диссертации без ссылок на чужие работы, принадлежат лично автору.

Апробация работы. Основные результаты работы были представлены и докладывались на: VJ1 Всесоюзном семинаре по комплексам программ математической физике (Горький, 1981); семинаре но численным методам механики сплошной среды под руководством академика Н.Н.Яненко (Новосибирск, 1982); Всесоюзных семинарах но комплексам программ математической физики (Новосибирск 1982, Ташкент 1983, Новосибирск 1984, Красноярск, 1988); объединенном заседании секций "Теоретические проблемы генерирования электромагнитной энергии" и "Проблемы теории поля и динамики в электроэнергетических и. электрофизических устройствах" Научного совега по комплексной проблеме "Научные основы электрофизики и электротехники" АН СССР отделения физико-технических проблем энергетики по теме "Численный эксперимент при создании мощных электрических машин" под руководством академика И.А.Глебова (Ленинград, 1982); II Всесоюзном совещании "Автоматизация проектирования и конструирования" (Ленинград, 1983); симпозиуме "Математическое обеспечение для автоматизации исследований, идентификации и планирования эксперимента" (Харьков, 1984); республиканской научно-технической конференции "Современные проблемы энергетики" (Киев, 1985); Всесоюзном научно-техническом семинаре "Опыт применения средств технической диагностики и контроля за состоянием электроэнергетического оборудования" (Суздаль, 1986), VIII Всесоюзной конференции "Планирование и автоматизация эксперимента в научных исследованиях" (Ленинград, 1986); Всесоюзной конференции "Актуальные проблемы вычислительной и прикладной математики" (Новосибирск, 1987); VI Всесоюзном семинаре "Обратные задачи и идентификация процессов теплообмена")Москва, 1987); Всесоюзном научно-техническом совещании "Вопросы проектирования, исследования и производства мощных турбо-, гидрогенераторов и крупных электрических машин" (Ленинград, 1988); Всесоюзных конференциях "Математическое моделирование; нелинейные проблемы и вычислительная математика" (Звенигород, 1988, 1990); Всесоюзном научно-техническом совещании "Автоматизация проектирования и производства в электромашиностроении" (Суздаль, 1989); Всесоюзной научно-технической сессии "Состояние и перспективы газодинамических тепловых исследований в обеспечении повышения температуры газа в стационарных газотурбинных установках" (Москва, 1989); научном семинаре в техническом университете г.Хемниц (Германия, 1990); Всесоюзной конференции по методам численного решения многомерных нестационарных задач магемапце-ской физики (Арзамае-16, 1991); Международной геофизической конференции и выставке по разведочной геофизике SEG-EAOO (Москва, 1993); Международном совещании-семинаре по механике реагирующих сред и экологии (Томск, 1994); First Asian Computational Fluid Dinamics Conference (Hong Kong, 1995); Particle Accelerator Conference and international Conference on High-Energy

Accelerators (Dallas, 1995); Международной конференции PaCT-95 (Санкт-Петербург, 1995); Международной геофизической конференции и выставке Санкт-Петербург'95; Международной конференции "Неклассическая геоэлектрика" (Саратов, 1995); Международной конференции "Математические модели и численные методы механики сплошных сред" (Новосибирск, 1996); Научно-техническом совещании "Геофизические методы при разведке недр и экологических исследованиях" (Томск, 1996); Международной геофизической конференции "Электромагнитные исследования с контролируемыми источниками" (С.-Петербург, 1996); научных семинарах ВЦ СО РАН, ИВТ СО РАН (Новосибирск). Результаты диссертации включены в 6 зарегистрированных от-четов'по НИР НЭТИ и в 2 заключительных отчета СНИИГГиМСа.

Публикации. По результатам выполненных исследований опубликовано 46 печатных работ.

Структура работы. Диссертация состоит из введения, семи глав, заключения и списка использованных источников (208 наименований). Работа изложена на 333 страницах, включая 39 рисунков.

Основное содержание работы

Глава 1. Вычислительные схемы пысокогочного коиечноэлементного моделирования трехмерных стационарных электрических нолей в проводящих средах

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

где X - проводимость среды;/- плотность источников электрического коля, определяемая стекающими в описываемую областью П среду токами; Л'], -внешние границы области С2.

Решение ('краевой задачи (1)-(2) ищется в виде суммы V -Vй т1''*. При этом сначала ищется решение V0 краевой задачи для уравнения

~div(Xgrad V) =/,

0)

>Y=o.

- -I =о,

И,

(2)

в области П с новой проводимостью }?-Хи(х,у,:) и старым (взятым из уравнения (1)) источниковым членом/ Затем решается краевая задача для уравнения

- сЦ>. цгас1Г) = -<1 Ц(Г - Г). (4)

в области П с исходной проводимостью Х--).(:(,у,:), но новым источниковым членом Очевидно, что сумма этих решений К=4'-»К

удовлетворяет уравнению (I). Краевые условия в краевых задачах для уравнений (3) и (4) должны быть абсолютно аналогичны краевым условиям (2) в краевой задаче для уравнения (I)..

Далее в п.1.2 и л. 1.3 диссертационной работы обосновывается, что вычисление решения Г исходной краевой задачи для уравнения (I) в виде суммы основного (аппроксимирующего) поля Г" и добавочного поля V* оказывается очень полезным для существенного снижения вычислительных затрат ири решении задач электроразведки меюдом конечных элементов. '

В п. 1.3 приводится вариационный эквивалент краевой задачи для добавочного поля в виде;

Ji.grас]Г' цгас1ч1 </а = - л)йга<.1Г" цЫу .

(5)

В вариационном уравнении (5) источники поля У* (определяемые по решению краевой задачи для основного поля ('") учитываются в виде объемных интегралов. Далее в этом же пункте рассматривается конечноэлементная аппроксимация задачи для добавочного поля и приводятся формулы вычисления локальных матриц и векторов конечноэлементной СЛАУ.

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

В п.1.5 приводятся эквивалентная вариационная формулировка и конечно-элементная аппроксимация краевой задачи для добавочного поля, в которых все источники добавочного поля учитываются в виде поверхностных интегралов.

В п.1.6 рассматриваются особенности и проводится сравнительный анализ эквивалентных вариационных формулировок задачи для добавочного поля1 с источниками в виде объемных и поверхностных интегралов. Исследуется влияние ошибок представления основного поля на ошибки вычисления добавочно! о поля. Для вариационной формулировки с источниками в виде поверхностных интегралов описываются процедуры коррекции этих источников, уоршшощие

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

В п. 1.7 рассматриваются основные принципы построения близких к оптимальным конечноэлементных сеток при решении задачи для добавочного поля с учетом особенностей локализации его источников.

В п. 1.8 на примере решения достаточно типичной задачи наземно-скважишюй электроразведки продемонстрированы преимущества рассмотренной в первой главе методики конечноэлементного моделирования. Эта задача заключается в следующем. Электрическое поле изучается в среде с проводимостью 0.2 См/м. Исследуемая среда включает в себя две неоднородности: тело формы

. I !оложте.п.миИ источник <?=<>м)

|Дне1шая поверхность

Jäa&auiL-

R=77 5мм .

0.089 м 0.0775 м имеющий I О7 См/м

Ч9()м

JiiLttÜU

! 1ефтяниИ пласт (о=А (ЕСм/ч)

цилиндрической (описывающее нефтяной пласт) высотой 10 м и радиусом 3000 м, расположенное на глубине 1000 м и имеющее проводимость 0.02 См/м, и соосный с ним усеченный полый цилиндр высотой 1050 м с радиусами (внешний) и (внутренний), проводимость (усеченный полый цилиндр описывает обсадную колонну труб). Схема разреза этой среды приведена на рис.1. Возбуждающий электрическое поле ток циркулирует между двумя электродами, расположенными в обсадной колонне труб на глубине 1000 м и на поверхности Земли. Цель исследования состоит в определении влияния нефтяного пласта па электрическое поле, изучаемое на поверхности Земли (то есть на так называемой "дневной" поверхности).

Результаты проведенных численных исследований представлены на рис.2.

. OipUlüüLMI.III.lÜ источник. (г—1000м)

Обсадим ко.юина труб' (ст=10'См/м)

Вмещающая среда . (а=0.1См/м) ;

Рис. 1. Схема разреза среды для модельной задачи с нефтяным пластом

V - V

Кривая I на рис.2 показывает зависимость функции = - -"■-• 100% от

радиальной координаты г (100м < г< 5000м) на дневной поверхности, где У0,/ -решение осесимметрнчиой задачи с нефтяным пластом, а Ут, - решение осе-симметрнчной задачи без нефтяного пласта (то есть в ней нефтяной пласт заменен вмещающей средой). Таким образом, функция 8 У,,,/ фактически показывает влияние нефтяного пласта на исследуемое электрическое поле. Функции ('„,/ и К„„ были рассчитаны с использованием МКЗ на очень подробных треугольных сетках, так что погрешность вычисления 5Г„,/ составляла около 0.2% - 0.3%, и поэтому кривая I может быть принята за эталон.

У*

Кривая 2 на рис.2 показываег зависимость функции = _ -100% от

К*,

радиальной координаты г (на дневной поверхности), где У - добавочное поле, вычисленное как решение трехмерной краевой задачи для уравнения (4), в которой в качестве функции У" была использована функция Уи.„. Расчет поля У+ выполнялся на довольно грубой конечноэлементиой сетке (в качестве конечных элементов были использованы тетраэдры). Как видно из сравнения кривой 2 с эталонной кривой 1 на рис.2, решение трехмерной задачи для добавочного поля даже на довольно грубой сетке позволяет достаточно точно оценить влияние нефтяного пласта на изучаемое электрическое поле,

У' -У

Кривая 3 на рис.2 показывает зависимость функции &Уи), = '"' ^ ■ 100%

от координаты г на дневной поверхности, где У,', - стандартное конечноэле-ментное решение исходной трехмерной краевой задачи для уравнения (1) на той же тетраэдральной сетке, которая была использована при вычислении добавочного поля У*. Кривая 4 на рис.2 показывает зависимость функции у2 -К

5Уа1 = —-100% от г на дневной поверхности, где У\ - стандартное ко-У.е

нечноэлементное решение исходной краевой задачи на тетраэдральной сетке, содержащей в 8 раз большее количество узлов по сравнению с сеткой, использованной при вычислении поля У*, а кривая 5 является графиком функции

ЪУ01 = -г!—.^¿е . Ю0% от г, где У', - стандартное конечно-элементное решение

исходной краевой задачи на тетраэдральной сетке, содержащей более чем в 20 раз большее количество узлов по сравнению с сеткой, на которой вычислялось поле У*.

Очевидно, что использование сечки, на которой было получено решение Г* с точностью, вполне достаточной для оценки влияния нефтяного пласта на изучаемое электрическое поле, дает совершенно неприемлемые результаты при попытке решить эту задачу с применением стандартного подхода к построению конечноэлементного решения (гю кривой 3 на рис.2 совершенно нево ¡можно

сколь-нибудь уверенно судить о влиянии нефтяного пласта на изучаемое электрическое поле). Увеличение числа узлов в сетке в 8 раз, повлекшее за собой увеличение затрат оперативной памяти в 8 раз и времени счета более чем в 20 раз, хотя и дало возможность получить некоторое представление об исследуемой зависимости, но еще очень сильно уступает по точности исследованию, проведенному через вычисление добавочного поля. И только увеличение числа узлов более чем в 20 раз н последовавшее за этим увеличение затрат оперативной памяти более чем в 20 раз и времени счета более чем в 100 раз позволило вычислить влияние нефтяного пласта с точностью, достаточно близкой к той, которая была достигнута при расчете добавочного поля путем решения краевой задачи для уравнения (4) на очень грубой сетке (хотя точность расчета через решение краевой задачи для добавочного поля все-таки заметно выше; выполнить же дальнейшее дробление сетки оказалось невозможным потому, что были полностью исчерпаны ресурсы оперативной памяти компьютера, на котором проводилось данное исследование).

Рис. 2. Результаты вычисления оценки влияния нефтяного пласта на электрическое поле на дневной поверхности при решении трехмерной задачи напрямую в сопоставлении с результатами расчета при использовании методики раздельного вычисления основного и добавочного полей

Кроме того, значительного сокращения времени счета при решении задачи для добавочного поля можно добиться, если выполнить некоторое "размазывание" обсадной колонны труб, заключающееся в увеличении площади поперечного сечения трубы и соответствующей коррекции ее проводимости. Такая процедура существенно улучшает свойства матрицы конечноэлементной СЛАУ, благодаря чему время счета при решении этой СЛАУ методом сопряженных градиентов (с предобусловливанием, использующим неполное разложение Холесского матрицы СЛАУ) сокращается на порядок и более (при одном и том же количестве узлов в конечноэлементной сетке). Но решение многих модельных задач показало, что применение процедуры "размазывания" трубы корректно лишь в задаче для добавочного поля н при условии, что "размазываемая" труба точно (без "размазывания") учтена в осеснмметричной задаче для основного поля. Проведение процедуры "размазывания" трубы в задаче для основного поля или в исходной задаче для искомого поля (решаемой напрямую) очень часто приводит к появлению существенной систематической ошибки в значениях потенциала, возникающей из-за фактического изменения математической модели изучаемого физического процесса, и поэтому приме-пять эту процедуру при решении исходной трехмерной задачи напрямую, вообще говоря, нельзя. Таким образом, эффективность разработанной методики вычисления искомого поля в виде суммы основного и добавочного полей становится в таких случаях еще более высокой.

Далее в п. 1.8 приводятся результаты применения разработанной методики конечноэлементного моделирования электрических полей при решении еще одной модельной задачи и нескольких практических задач наземно-скважинной электроразведки и каротажа. ,

Глава 2. Особенности конечноэлементного решения трехмерных задач электростатики при использовании методики выделения двумерной составляющей потенциала изучаемого электрического ноля

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

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

Вторая особенность заключается в.следующем. В задачах электростатики не всегда удается достаточно просто и естественно-выпилишь сформулирован-

мое п начале п 2.2 диссертационной работы требование о полном включении трехмерной области (I, в которой определена исходная задача, в близкую к ней область Ц) с двумерной (осесиммегричной) геометрией. Связанные с этим трудности могут возникнуть тогда, когда в исходной электростатической задаче проводник, на котором задано некоторое значение потенциала электрического поля, имеет геометрию, близкую к двумерной (осесиммегричной), по никак не может быть учтен п двумерной задаче таким образом, чтобы все подобласти исходной трехмерной области П (в которых нужно получить распределение потенциала) оказались бы внутри области Q« с двумерной геометрией. II тогда для выполнения требования о включении Q в П<> необходимо /(сформировывать область ß,i подобластями, со всех сторон окруженными границами с заданными на mix краевыми условиями первого рода.

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

Потенциал !'" электрического поля в осесиммегричной задаче может быть доопределен внутри отверстий каждого электрода значениями, равными значениям потенциала исходного поля на поверхностях электродов, преле чего этот доопределенный потенцнтт Г" может быть использован в правых частях дифференциального уравнения и краевых условий при решении краевой задачи для добавочного поля Г*. При этом обратим внимание на то, что вид правой части дифференциального уравнения в краевой задаче для добавочного ноля V* дол-ясен быть соответствующим образом изменен, так как после продолжения К0 внутрь каждого электрода в расширенной области определения I'0 (точнее, в окрестности границ нерасширенной области определения К0, попавших внутрь расширенной области определения У0), может не выполняться дифференциальное уравнение для основного поля. Поэтому в краевой задаче для добавочного поля F+ лиффереициальиое уравнение примет вид (сравните с (4)):

- div(A.gmdr)= -div((A.° - Ä.)gradF°) + f, (6)

тле дополнительные источники /р поля V*, появившиеся из-за расширения области определения потенциала V", определяются соотношением

Нетрудно убедиться, что в рассматриваемом случае соотношение (7) фактически определяет поверхностные источники, заданные на границах Л',1' нерасширенной области определения Г°, которые после продолжения потенциала V" попали внутрь расширенной области его определения. Интенсивности этих поверхностных источников определяются значениями производной потенциала I вдоль нормалей к Б'. С учетом этого эквивалентная вариационная постановка для уравнения (6) примет вид

" . „ (8) £1 X'

где б1' - объединение всех границ нерасширенной области определения У", попавших внутрь расширенной области его определения, а Л - внутренняя по отношению к нерасширенной области определения У0 нормаль к .

В п.2.4 на примере решения существенно трехмерной практической задачи (общий вид расчетной области с заданной в пей конечноэлементнои сеткой приведен на рис.3) подтверждается работоспособность и эффективность методики вычисления электрического. поля в виде суммы основного и добавочного полей в задачах электростатики.

На рис.4 представлены результаты расчета изучаемого электрического поля. Кривая 1 на рис.4 показывает распределение компоненты Е* напряженности Ё* - - угас! V * добавочного поля У* вдоль отрезка прямой дг=0, >-0, Кривая 2 на рис.4 показывает распределение вдоль того же отрезка прямой функции ЬЕ = Е, - где Ег вычислена по потенциалу У, полученному при решении исходной трехмерной задачи напрямую с использованием той же сетки, на которой решалась задача дня добавочного поля. Кривые 3 и 4 на рис.4 также показывают распределение функции 6Е= Е. но для этих кривых компонента Ег была получена по значениям потенциала V, вычисленным при решении исходной трехмерной задачи напрямую с использованием сеток, содержащих соответственно в 10 раз (кривая 3) и в 64 раза (кривая 4) большее число узлов, чем сетка, на которой решалась задача для добавочного поля.

Рис.3. Коиечноэлементвая дискретизация расчетной области

Рис. 4. Распределение функций Е* и (Л. — вдоль отрезка прям-ч, определяемо» уравнениям» г=0,_у=0

Таким образом, для получения с необходимой для практических нужд точностью значений напряженности электрического поля при решении исходной задачи напрямую потребовалась тетраэдрачьная сетка, содержащая в 64 раза большее число узлов по сравнению с числом узлов в сетке, использованной для расчета /; через добавочное поле С учетом затрат на решение осесимметрич-ной задачи при вычислен::и Е в виде суммы Е° +■ Е* вычислительные затраты па решение исходной задачи напрямую возросли почти в 30 раз по памяти компьютера и более чем в 100 раз по времени счета по сравнению с суммарными вычислительными затратами, потребовавшимися на вычисление Е в виде суммы Е° + Е* (в данном случае осесимметрнчная задача требует затрат ресурсов компьютера даже больше, чем трехмерная задача для добавочного поля).

Глава 3. Вычислительные c.vcsiu конечноглемсшного моде.тиропанни нестационарных электромагнитных полей в средах с трехмерными неодмо-родиостими

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

В п3.1 приводятся и анализируются различные математические модели, описывающие трехмерные нестационарные электромагнит ные поля.

В ситуации, когда нельзя пренебречь токами смещения или коэффициент магнитной проницаемости среды р отличен or коэффициента магнитной проницаемости вакуума (то есть среда неоднородна по |д) для описания электромагнитного поля используется уравнение

( 1 -Л 5 Л д*Л :

rot ■ rot/1 +о +с -,=■!, (9)

кц J-didi2

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

В = ю14, Н = ~дА ' (10)

д/

( Е - напряженность электрического поля, В - индукция магнитного ноля).

Если же токами смещения можно пренебречь и среда однородна по ц, то для описания электромагнитного поля используются уравнения

-±M + o\d£ + gradA = 3 , (11)

Ро J

-div(agradH)-div^cr^j = -div.7 , (12)

причем напряженность электрического поля /:" определяется соотношением

=- -gradi'. (13)

dt

В 11.3.2 и п.3.3 рассматриваются вычислительные схемы МКЗ-моделирования трехмерных нестационарных электромагнитных полей, оено-

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

О п.3.2 рассматриваются ситуации, когда для описания электромагнитного поля используется уравнение (9). Вектор-потенциал А представляется в виде суммы двух вектор-потенциалов А = ,4° + Я*, где А° - вектор-потенциал, удовлетворяющий уравнению

а А* - вектор-потенциал, удовлетворяющий уравнению

^ > д! Х 1 д!2

(15)

где 3* = ,7 - .7°.

Далее проводится конечно-разностная аппроксимация уравнения (15) по времени с использованием трехслойной полностью неявной схемы и приводится эквивалентная вариационная постановка в виде:

Г1 рЫЧ'гЮ- Г-- — - - цгас! + =

Ju " Ju й .1 ^

i) ч с1

а и

+ |(л - 0^ЛЫ - + ОД" + (укы + , (16)

¡1

где 4 - одна из пространственных переме :ых ( г, у или г); 0' - значение вектор-потенциала А* на}-\\ временном слое, а IV1 - значение вектор-пепетшала А0 на у-м временном слое (го есть 0'(х,у,г) =

№'(х,у,г)= 1с - номер текущего временного слоя; пробная функ-

ция Ч' - произвольный элемент гильбертова пространства Т^1 функций, определенных в трехмерной области С2, имеющих суммируемые с квадратом первые производные и обращающихся в нуль на границе 5 области П; 0; и 0( - пара-

метры, определяемые соотношениями

е,в(о--о)^(0 + (в»-е)^(0. (18)

Входящие в формулы (17Н18) функции 1)1(/), гь(')" ч'О) являются квадратичными полиномами переменной I на интервале (/к 2, 4):

Значения коэффициентов с, 1\ ис| этих полиномов вычисляются так, чтобы выполнились следующие равенства:

= тМ-.М,,, Ч'(0 = 5,,, 1 = 1-----3 , (20)

где 8,; - символ Кронексра.

Далее в п.3.2 диссертационной работы выводятся формулы для вычисления локальных матриц и векторов конечноэлементнон СЛАУ, аппроксимирующей систему вариационных уравнений (16).

В п.3.3 рассматриваются ситуации, когда для описания электромагнитного поля используется система уравнении

(1*1 И12). Решение (Л,!') системы уравнений (11)-(12) ищется в виде суммы двух решений: решения двумерной (осесимметричной) краевой задачи, а также решения трехмерной краевой задачи для системы уравнений

----Л/Г . 4'''' +йга<1!-") =(а"-а/^*-'' + (21)

р0 V д! ) х \ д! )

(а4 - а)д? - Шу[(о° - а)ВЫ1/0] + Г ,

(22)

где Г =-111 у.Г.

Далее проводится конечно-разностная аппроксимация уравнений (21)-(22) гю времени с использованием трехслойной полностью неявной схемы и приводится эквивалентная вариационная постановка н виде:

"»а

01 • цЫ Ч'с/а + |()/71 Ус/п + рте! V • Ч'</Г2 = а и а

= |(.г - 02(7к~' - о//1-2 + е,»7' + 9,и"+ о,(Гл ,

1}

|а£га(Н" £гас1Ч'гЮ + |о/7к угаеПУП = р^Ч'ЛП + '.) 11 И

+ |(-0,{/к,-е,/71:! +0,И"'к +01((:'к-5)ргас1Ч'с/П+ (24)

о

+ |"(а° - а) цЫ Г° цгай ЧШ, а

гдев-а^ЧО, М00"0)-^)'

В заключение пункта 3.3 диссертационной работы приводятся ({юрмулы для вычисления локальных матриц и векторов конечноэлементнон СЛАУ, аппроксимирующей систему вариационных уравнений (23)-(24).

В п.3.4 анализируются недостатки и достоинства рассмотренных подходов, обсуждаются особенности и трудности реализации соответствующих вычислительных схем и проводится сравнение вычислительных затрат.

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

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

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

Глава 4. Вычислительные схемы МКЭ при изучении влияния вихревых токов на нестационарные электромагнитные поля в технических устройствах

Четвертая глава посвящена разработке эквивалентной вариационной формулировки и конечноэлементноб аппроксимации задач изучения влияния вихревых токов на нестационарные электромагнитные поля в технических устройствах.

В п.4.1 рассматривается математическая модель изучаемых электромагнитных процессов в следующем виде. В подобласти О« с нулевой проводимостью индукция магнитного поля представляется в виде градиента скалярного потенциала:

= -gradí/ , (25)

в результате чего автоматически обеспечивается выполнение уравнения гоШ = 0. Уравнение же d¡v/)' = 0 из системы уравнений Максвелла принимает вид

-div(grad(/) = 0. (26)

В остальных подобластях fii,...,Q¡s расчетной области D индукция магнитного поля й представляется через вектор-потенциал А с помощью соотношения

/>' = rol Л. ' (27)

Полагая, что напряженность электрического поля определяется через вектор-потенциал А с помощью второго из соотношений (10), убеждаемся, что при таком представлении В и Н автоматически выполняются уравнения xüXE--dli dt и div/i = 0 из системы уравнений Максвелла. Уравнение же rot 11 - J + ali1 преобразуется к виду :

rotí-rot4| + cf=./. (28)

J di

Таким образом, в подобласти Ц, с нулевой проводимостью и магии гной проницаемостью р« индукция магнитного поля И1' определяется через скалярный потенциал V и магнитное поле описывается уравнением (26), а в подобластях <П],...,£}р с проводимоетями а„#0 и магнитными проиицаемостммн (а, (а=1,...,Р) индуншя В'' определяется через вектор-потенциал А и электрома!-нитное поле описывается уравнением (28). Чтобы индукция Й, определяемая соотношением

[В1 в«,,...,П,, ,

удовлетворяла системе уравнений Максвелла во всей расчетной области £2=иС2„, на границах ^....Л между подобластью По и подобластями

г,-О

(соответственно) должны выполняться следующие соотношение (обеспечивающие непрерывность нормальной составляющей В и тзнгешдиаль-

ных составляющих Н = ): /'

= (зо)

= (31)

и» И

где п - лгоб.п (например, внешняя) нормаль к рассматриваемой границе (а 5 1), а т - произвольный касательный к границе 5а (го есть ортогонзлмшП к се нормали п) вектор.

В и 4.2 приводится эквивалентная взриациош'зя формулпровга для рзс-смотренной пышс математической модели в виде:

("--- Ёгас!^,егас1Ч'\Я2- Г' ЙЛ впкГГа'П +

•Ч1„ •|М„ дх

а, «.

(32)

+ N = Г,.™,,

Лцо'чОу ' Зг V Л 81 )

„■'Ра дУ

с 1 ( би еи

+ I -------п, + --

■ЬтД дх 2 Ъг

\ л »

(33)

Г 1 ь'ггиЫ - [ 1 ^ ^атГГ./О т-

¿Л &

ц

Но

7 дЛ

------ п, + ' п

{ дх ' д: •

( ЗА, 1 ду

ал, ' п ду

вЛк )

' "< - & Ч'

= о,

(35)

1де ,Ч'У,Ч'г 1! 1Г - произвольные элементы соответствующих гильбертовых пространств , Я<1М , и .

Далее в п.4.2 диссертационной работы доказывается, что если вектор-потенциал А и скалярный потенциал (/ удовлетворяют системе вариационных уравнений (32)-(35), ю эго обеспечивает выполнение не только дифференциальных уравнений (26) н (28), но и условий (30) и (31) на границах Л'и стыковки .¿верного потенциала II и вектор-потенциала А .

В ¡1,43 рассматриваются наиболее часто встречающиеся при решении практических задач условия симметрии электромат нмтных нолей в технических конструкциях. Для этих условий симметрии обосновываются процедуры их учета на плоскостях симметрии, параллельных координат ным плоскостям.

В п.4.4 рассматриваются некоторые особенности и связанные с ними модификации описанной в и.4.2 вариационной формулировки. Здесь доказывается, что эта вариационная формулировка обеспечивает выполнение условий не-протекашм токов через траницы между проводящими и непроводящими электрический ток частями конструкции. Далее обсуждаются различные ситуации, связанные с возможными появлениями разрывов в нроводтшостяч внутри тех подобластей, где индукция магнитного поля описывается векторным потенциалом. Кроме того, рассматриваются некоторые возможности существенного сокращения вычислительных затрат для определенною типа задач за счет задания краевых условий для отдельных компонент вектор-потенциала на границах подобластей расчетной области, а также возможности вычисления нолей влиянии отдельных трехмерных объектов. *

И п.4.5 провалится дискретизация по времени и копечпоэлемеитная аппроксимация по пространству системы вариационных уравнений, описывающей нестационарные электромагнитные поля в технических устройствах. Приводятся формулы, позволяющие вычислять вклады в локальные матрицы и векторы от различных, членов аппроксимируемой системы вариационных уравнений.

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

В п 4.7 приводятся результаты решения тестовой залачи, подтверждающей правильность разработанных в четвертой главе подходов к моделированию нестационарных электромагнитных иоле» - 'ечнических конструкциях.

В п.4.8 приводятся результаты исследования влияния вихревых токов иг электромагнитное поле в технической конструкции, обший вид которой пола зап па рис.5. Па рис.6 изображено рас предсление напряженности электрического поля в лпух сечениях расчетноГ области в момент включения тога I обмотках.

Рис. 5 Общей вид технического устройств?

Г

а) дс=0 0078

б) .Г-0.02П

Рис. 6 Распределение напряженности эчек-трического ноля на левой (при г=0.0078) и правой (при г-0.0213) сторонах пластины пакета

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

В пятой главе обобщаются рассмотренные в первых трех главах подходы к вычислению полей влияния трехмерных объектов и обосновывается подход, не требующий разработки специальных формулировок для задач расчета добавочного поля,

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

В п.5.1 рассматривается процедура вычисления полей влияния трехмерных объектов для задач, записанных в операторном виде:

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

Будем считать, что области действия дифференциально-краевых операторов I и ¡}' совпадают. Тогда поле влияния исследуемого объекта удовлетворяет операторному уравнению

Чтобы убедиться в этом, достаточно вычесть уравнение (37) из уравнения (36) и провести элементарные преобразования.

Таким образом, несмотря па то, что аппроксимирующая задача (37), как и исходная задача (36), является трехмерной, поле и* влияния исследуемою объ-еста можег быть найдено с достаточно высокой точностью из решения задачи (3&), если оператор 1° и источннковый член /"в задаче (37) довольно близш к оператору Ь и источниковому члену/из задачи (36). В этом случае источником поля и+ является в основном воздействие разности операторов /,-/," на поле и". Учет этого обстоятельства позволяет строить эффективные сежи как при численном решении аппроксимирующей задачи (37), ьтк и при численном решении

(36)

(37)

(38)

ется слагаемым го!

........ ~ ' '81

.1 (38), определяющей поле и¥ влияния исследуемого объекта. Очевидно, :я решения каждой из трехмерных задач (37) и (38) можно строить соот-_,юшую коиечноэлеменгную (или конечно-разностную) сетку, учитываю-теннфнку решаемой задачи. Эю, конечно, позволяет еще больше пэвы-сигь точность вычисления поля и (либо сократить вычислительные затраты при ...ши задач (37) и (38)), но, с другой стороны, создает и некоторые трудности в реализации алгоритма расчета 1/. Действительно, для решения задачи (38) поле и", полученное при решении задачи (37), должно быть корректно пересчитано на другую трехмерную сетку, используемую при аппроксимации задачи (38). Значительные трудности могут возникнуть н при конечноэлементнои (конечно-разностной) аппроксимации члена - Л°]н°'в задаче (37). Например,

в краевой задаче для уравнения (4) член -(/, - имеет вид

- ).^га(1 V"в краевой задаче для уравнения (15) этот член оиределя-

(( 1 II чо! I о \дАи ( „ чЭ2Л° У " ] 10 у 1° _а7 а7 +18 ' -¡2 > 3 в |Файв01! задаче для системы уравнений (21)-(22) - членами и

(сг°-а) ^ - Ш\|(а'' С"Поэтому имеет смысл рассмотреть

другой, значительно более простой с точки зрения численной реализации способ вычисления поля и* влияния исследуемого объекта при условии существенной трехмерности поля ¡Л Этот способ заключается в следующем.

На одной и той же сетке (конечноэлементной пли конечно-разностной) решаются две задачи: полная задача с учетом влияния исследуемою объекта (задача (36)) и такая же задача, но без учета влияния исследуемого объекта (задача (37)). Затем вычисляется разность решений этих задач г/* = и-и", являющаяся полем влияния исследуемого объекта. Несмотря на то, что ошибка аппроксимации в каждой из этих задач может быть сопоставима по уровню с полем и* (и даже превышать его), точность вычисления поля и' в виде разности полей и и и" оказывается довольно высокой и вполне достаточной для правильной оценки влияния исследуемого обьента. Этот факт можно объяснить взан-моунпчтожением основных частей ошибок аппроксимации задач (36) п (37) при вычислении разности их решений. Однако без теоретического обоснования этого факта и понимания порождающих его причин исследователь не имеет уверенности п правильности полученных им результатов. Рассмотренные в данной работе подходы к математическому моделированию сложных полей позволяют ие только дать теоретическое обоснование факту взаимоуничтожения ошибок аппроксимации в разности решений дифферешшааыю-краепых задач (36) и (37), но и прогнозировать уровень козможнок ошибки в поле влияния исследуемого объекта и даже определять пути ее уменьшения.

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

Будем считать, '¡то задачи (36) и (37) решаются с использованием МКЭ. В результате конечноэлементной аппроксимации этих задач получим две СЛАУ:

где 1. и 1° - копсчиоолсменгныс матрицы, являющиеся дискретными аналогами дифференциально-краевых операторов /. и / и /" - вектора правых частей конечнозлементпмх СЛАУ, а й и и" - вектора, являющиеся численными решениями дифференииалыю-краавых задач (36) и (37) соответственно. Поскольку СЛАУ (39) и (40) получены в результате конечноэлементной аппроксимации задач (36) и (37) на одной и той же сетке, мы можем провести над ними те же операции, которые проводили выше для операторных уравнении (36) и (37). Сначала пышем СЛАУ (40) из СЛАУ (39), а затем из левой и правой частей получившейся системы вычтем вектор /Л?0. После элементарных преобразований получим

где й* = й - й". Сравнивая (41) с (38) убеждаемся, что вектор ¡Г являеют решением СЛАУ, получающейся в результате конечноэлементной аппроксимации дифференциально-краевой задачи (38). Таким образом, при использовании одной и той же сетки для вычисления численных решений и и й° дифференциально-краевых задач (36) и (37) их разность й* является численным решением дифференциально-краевой задачи (38) на той же конечноэлементной сстке.

Итак, мы убедились, что ошибки аппроксимации в поле ¿Г влияния исследуемого объекта, полученного вычитанием численных решений й и и0 краевых задач (36) и (37), фактически являются ошибками аппроксимации, возникающими при численном решении дифферентшально-краевой задачи (38) на той же ссгке, которая была использована'при численном решении задач (36) и (37). Этот факт имеет важное практическое значение. Он позволяет дать обоснование следующей практической рекомендации. Для снижения ошибки аппроксимации

= / ,

¿V = ,

(39)

(40)

(41)

в поле н', получаемого в виде разности численных решении й и г/° задач (36) и (37), следует делать локальные сгущения узлов не только в местах резких изменений градиентов решений и и г/1, но и в местах резких изменений градиентов г/+ даже несмотря на то, что поля и и и° в этих местах могут быть довольно гладкими (поле и* может быть гораздо ниже по уровню полей и п ц", н поэтому даже резкие изменения производных в поле могу г слабо сказаться па гладкости и).

Таким образом, рассмотренные два способа вычисления полей влияния отдельных объектов, возбуждаемых существенно трехмерными источниками, имеют свои недостатки и достоинства Так, при вычислении поля и* путем решения дифференциально-краевой задачи (38) в процессе численного решения задач (37) и (38) могут быть использованы две разные сетки, учитывающие особенности каждой из этих задач. В результате каждая из задач (37) и (38) может быть решена на оптимальной сетке с минимальным количеством узлов, а эго часто позволяет существенно сократить требуемые вычислительные затраты (память и машинное время) без снижения точности вычисления поля и+. К недостаткам этого способа следует отнести то, чго для его применения необходимо разработать процедуры корректного пересчета решения трехмерной задачи с одной трехмерной сетки на другую, а также процедуры "вычисления окладов в правые части конечноэлементных (или конечно-разностных) СЛАУ от члена (¿-/.")|/' операторного уравнения (38), которые во многих случаях являются совсем не тривиальными.

Пторой способ вычисления и* в виде разности решений задач (36) и (37) не требует разработки никаких дополнительных процедур. Но сетка, используемая для численного решения задач (36) и (37), должна одновременно учитывать как особенности решения и" задачи (37), так и особенности решения и+ задачи (38) (хотя задача (38) в этом случае и не решается). Это приводит к тому, что для достижения требуемой точности в решении и* количество узлов в сетке возрастает, что влечет за собой увеличение затрат вычислительных ресурсов (памяти и времени счета). Поэтому выбор того или иного способа вычисления поля и+ (при условии существенной трехмерности поли н°) зависит от конкретной ситуации. если для решения задач (36) и (37) на сетке, учитывающей особенности и достаточно вычислительных ресурсов, го проще вычислять и* в виде разности решений и и и" задач (36> чт (37); если же оптимальные сетки для задачи (37) и задачи (38) существенно различаются, а вычислительных ресурсов для достижения требуемой точности не достаточно, то имеет смысл находить и* как решение задачи (38) со всеми вытекающими отсюда последствиями (разработкой процедур пересчета решения с сетки па сетку и учета члена воздействия разности операторов I, I? на функцию ?<° при решении задачи (38)).

В п.5.2 обсуждаются подходы, основанные на использовании методики выделения полей влияния отдельных обьектов для уточнения решения трехмерных задач. Эти подходы позволяют сочетать преимущества лолуаналигиче-ских методов при расчетах элекгромлпшгных полей в небольшом числе точек геометрически простой среды с возможностями МКЭ учета любых трехмерных

неоднородностей. Они очень удобны тем, что просты в реализации и не требуют никаких дополнительных доработок используемых в них процедур численного анализа.

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

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

Чтобы получить решение и исходной задачи с высокой точностью, применим описанную в предыдущем пунше методику вычисления полей влияния отдельных обьекгов. Для этого решим с использованием МКЭ две злтачи на одной н той же сетке: исходную задачу с учетом всех трехмерных неоднородностей и задачу в среде без трехмерных неоднородное!ей (то сеть вторая задача должна полностью соответствовать той задаче, которая может быть быстро и точно решена с применением какого-либо из иолуаналнтическич методов). Ко-нечноэлсментные решения этих двух задач обозначим соответственно через ¡7 и и" (обратим особое внимание па то, что коиечноэлемеитная сетка, используемая при вычислениях и и должна быть одна и та же). Тогда функция й, определяемая соотношением

й = й" + й - н°, (42)

будет аппроксимировать точное решение г/ с гораздо более высокой точнооью, чем конечно-элементное решение и : погрешность и-и может на порядок и более превосходить погрешность и - и (при условии, что погрешность и" - йя по-луаналитнческого метода решения задачи в среде без трехмерных неоднородностей также тта порядок или более ниже погрешности и-и конечноэлемент-ной аппроксимации исходной задачи при учете в ней всех трехмерных неоднородностей)- Обоснование этого утверждения базируется па фактах, приведенных выше. Действительно, выше было показано, что при использовании одной и той же сетки для вычисления конечноэлементиых решений и и и'1 их разность и - и" фактически аппроксимирует задачу для добавочного поля и* - и-и", и поэтому при не слишком высоком уровне влияния трехмерных неоднородностей основные части погрешности аппроксимации в функциях и и и" будут взаимно уничтожаться при вычитании и" из V .

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

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

В п.5.3 приводятся примеры решения двух модельных задач из различных классов задач электромагнетизма, иллюстрирующие возможности предлагаемой методики выделения полей влияния отдельных объектов при использовании ее для уточнения искомого решения.

В качестве первого примера рассмотрена-достаточно типичная задача на-земно-скяажшшой электроразведки, описанная в главе I (схема разреза среды для этой модельной задачи показана на рис.1). Для вычисления конечиоэле-ментных решений всех трехмерных задач в этом примере использована та же тетраэдральная сетка, что н при решении задачи для добавочного поля. На этой сетке решались как трехмерные две задачи: исходная задача с учетом нефтяного пласта, и аналогичная задача, в которой нефтяной пласт заменен вмещающей средой. Затем из решения трехмерной задачи с нефтяным пластом вычиталось решение трехмерной задачи без нефтяного пласта, после чего эта разность добавлялась к эталонному решению задачи без нефтяного пласта (согласно формуле (42); эталонное решение вычислялось как решение соответствующей осе-симметричнои задачи). Покажем, что полученная в результате этих действий функция действительно является значительно уточненным решением исходной задачи. Для этого сравним точность решения нескольких задач на плоскости г-0 (то есть на так называемой "дневной поверхности") с помощью [рафиков, приведенных на рис.7.

Рис. 7. Роульгагы решения двух трехмерных задач и вычисления оценки влияния нефтяного anacía на электрическое иоле па дневной поверхности в виде их разности в сопоставлении с результатами paciera при использовании методики ра¡дельного вычисления основного и .

добавочного полей

Кривая I на рис. 7 показывает эталонную зависимость функции у -

8К = —--„-—• 100% от координаты г, где У - решение осесиммегричиои задачи

с нефтяным пластом, а У0 - решение осесимметричной задачи без нефтяного пласта (то есть функция 5Г показывает влияние нефтяного пласта на исследуемое электрическое поле; эта функция фактически я пляс гея функцией 5 У,„г, и поэтому кривая 1 на рис.7 полностью совпадает с кривой I на рис.2). Кривая 2 на

У*

рис.7 показывает зависимость функции 51'- 100% ог г, в которой функция

У

Vf вычислена как решение задачи для добавочного поля (кривая 2 на рис.7 также полностью совпадает с кривой 2 на рис.2). Кривая .3 на рис.7 показывает за- Г - У9

виснмоеть функции 51 = - 100%, где У -решение исходной трехмерной

рй __ у О

задачи, а кривая 4 - зависимость функции й!/0 = —^—100% ог г, где Г1 -

стенке трехмерной задачи без нефтяного пласта на той же сетке, на которой ..илучено решение I'. Заметим, что фактически кривая 4 показывает относительную погрешность решения аппроксимирующей задачи на трехмерной сетке. И, наконец, кривая 5 на рнс.7 показывает зависимость функции

8Г = /0 • 100% от г, где V = Г" + V-V" - полученное по изложенной выше

\

методике уточнения решение исходной задачи (см. формулу (42)).

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

Для второго примера, рассмотренного в п.5.3, пол».'«ено столь же значительное увеличение точности расчета нестационарного электромагнитного поли с использованием методики конечноэлементного моделирования, описанной в главе 4.

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

Глава 6. Использование схемы против потока при расчетах методом конечных элементов тепловых полей с учет ом кчнтек-пшиого тендттереко-са

Шестая глава посвящена разработке вычислительной процедуры, реализующей так называемую схему против потока при конечпоэлеменгной аппроксимации краевых задач, описывающих тепловые поля с учетом конвективного тсплопсреноса.

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

д" ди си ,. /, . ч ,

; + г- + — -- 111\(Я.йга£1 и)~ {, (43)

ох су ' д: ' "

где функция и является искомой температурой,/- источником тепловыделения, X - коэффициентом теплопроводности, а вектор-функция w = (»', ,и'г.и'..-) задает характеристики конвективного тетшоперепоса (определяемые плотностью, тсппоемкостыо и скоростью течения жидкости или газа).

¡Сонечноэлементная аппроксимация краевой задачи для уравнения (43) по методу Галеркина часто приводит к негативным (иногда даже физически неправдоподобным) результатам, например, при больших значениях компонент вектора w и недостаточно мелкой сетке. Таких негативных результатов можно избежать, если использовать для аппроксимации уравнения (43) схему против потока (называемую также схемой с разностями против потока, схемой с разностями против течения и т.д.). Эта схема довольно широко и успешно применяется для конечно-разностной аппроксимации уравнения (43).

Рассмотрим возможность применения схемы против потока в ситуации, когда аппроксимация краевой задачи для уравнения (43) выполняется с использованием МКЭ. Будем считать, что трехмерная область, в которой определено уравнение (43), разбита Tía тетраэдральные конечные элементы í\, а в качестве

базисных функций \|/¡ в разложении приближенного решения г/н нс-

j

пользуются линейные на каждом конечном элементе Í2* полиномы:

+ + (44)

Физически неправдоподобных результатов при аппроксимации по методу Галеркина краевой задачи для уравнения (43) не возникает, если выполняются следующие свойства полученной в результате этой аппроксимации СЛАУ: все диагональные элементы матрицы А конечноэлементной СЛАУ являются ноло-жятельнымн, а все внедиагоналыше - неположительными (отрицательными или нулевыми). При использовании некоторой фиксированной сетки выполнение этих свойств матрицы А конечноэлементной СЛАУ, полученной в результате применения стандартной схемы аппроксимации уравнения (43) по методу Галеркина, зависит только ог значения характеристик конвективного тешюие-реноса й' (точнее, от значения отношения й'Д). Действительно, если обозначить через среднее на Пц значение компоненты нг вектора »■, то в результа-

_ ди

тс аппроксимации по методу 1алсркина вклад от члена в1 —- в компоненту а,,

д.х

матрицы А конечноэлементной СЛАУ на конечном элементе может быть . вычислен по формуле

Гй' ' ¥ </П, (45)

3 Лг

тле ц/, - базисная функция, соответствующая /'-му ywy конечноэлементной сет-

34

кн. Учитывая (44), из (45) получаем яклад в компоненту а,у матрицы А конечно-элементной СЛЛУ:

КЛУ(Я2, (46)

.1,1 г- д"

где = (аналогичная запись может оыть получена для членов н>у у или

()и , I , ,1 , 1 ,

К'. -- И в НИХ Ь) =Ьу ИЛИ Ь] =0(Д

Поскольку базисные функции неотрицательны (по определению), то

знак вклада (46) от члена уг в компоненту ач на конечном элементе О* оп-

дх

ределяетея знаками компоненты IV* конвективного теплопереноса и коэффнци-в представлении (44) базисной функции на Дь Тогда вклад (46) при

:аточно большом (по абсолютной величине) м'х мож^г нарушить положительность диагонального элемента а„ матрицы А конечно-элементной СЛАУ, если Л'Л* < 0 на конечном элементе а также неположительность внедиаго-нального элемента я,/, если >0 на Это и является основной причиной

появления физически неправдоподобного решения задачи, а также возможного существенного ухудшения сходимости итерационных методов решения конеч-ноэлеменгных СЛАУ. Аналогичные рассуждения справедливы и для вкладов в

. -^гтллг ди

компоненты а,, матрицы А конечноэлементнои СЛАУ от членов — и и' —

ду " дг

уравнения (43).

Таким образом, чтобы избежать появления физически неправдоподобных результатов и существенного ухудшения свойств конечноэлементнои СЛАУ при выполнении конечноэлементнои аппроксимации краевой задачи для уравнения (43), нужно добиться того, чтобы вклады (46) от члена и>л (и анало-

ах

гичные вклады от членов IV — и и'_ ) уравнения (43) не нарушали положи-

ду ' д:

тельности диагональных элементов и неположительности внедиатональных элементов матрицы Л конечноэлементнои СЛАУ. Для достижения этой цели применим следующую схему.

Обозначим через А",- множество номеров всех конечных элементов, содержащих у-ый узел конечноэлементнои сетки, а через О, - конечный носитель базисной функции то есть = У О, . Для каждого / выберем такой ко-

1 А ьЛ*,

нечнын элемент Ик , А, е К1, для которого

¿М*' >0,

и1'' • /)/ < О для j ф /,

где у - номера всех остальных (помимо /-го) узлов конечного элемента £2, .

Вклады от члена и>т —- уравнения (43) в компоненту я,у матрицы /1 конечноэле-дл

мснтной СЛАУ на каждом конечном элементе 0* для к е. К, (то есть Па принадлежит носителю ПА. базисной функции »)/,) будем вычислять по формуле

где j - номер одного из узлов конечного элемента (а не Пд). Эго означает , что соотношение (48) (в отличие от соотношения (46)) дает вклады в элементы а,У матрицы Л конечноэлементиои СЛАУ с номерами столбцов j, соответствующими номерам узлов конечного элемента Г1t (а не £2*).

Таким образом, рассматриваемая нами схема отличаете,' тг стандартной аппроксимации по методу Г'алеркина тем, что на конечном элементе П* член

ди

уравнения (43) аппроксимируется через средние и узловые значения

функций irt и ди/дх, вычисляемые не на самом конечном элементе Па, а на соседнем с ним конечном элементе П4 (принадлежащем, как и конечному носителю ilKi базисной функции ip,), для которого выполняются соотношения (47).

Нетрудно убедиться, что соотношения (47) гарантированно выполняются для такого конечного элемента, в котором вершина с номером / лежит ниже по потоку вдоль оси дг, чем остальные вершины этого Конечного элемента (то есть если ,vx положительно, то ;-ая вершина должна лежать правее всех остальных, а если wx отрицательно - то левее). Очевидно, это и является основным свойством схемы против потока: для аппроксимации функций и1, и ди'дх при вычислении вклада в ;-ое уравнение конечноэлементной СЛАУ от члена w — выби-

д.х

рается тот конечный элемент, в котором /-ая вершина лежит ниже по течению

вдоль оси л', чем все остальные вершины этого конечного элемента. Аналогично

должны быть подобраны конечные элементы для аппроксимации функций w„

ди ди ди ди

и иг, --- при вычислении вкладов от членов и> — - и ir - —.

ду dz ду " 0:

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

(48)

поток жидкости (газа) резко меняет направление движения. В этих случаях более

точной (для слагаемого и» - — ) является аппроксимация следующего вида: дх

(49)

где и', вычисляется на конечном элементе О*, а не на П(1 (функция и*, может

быть заменена, например, своим ннтерполянтом на Г}*)• Тогда конечный элемент П, должен быть выбран таким, чтобы выполнялись соотношения

Л,4' J^hm/Í} > 0 ,

'Ч

Л'- ¿ 0 для у' *«',

(50)

где, как и раньше, j-- номера остальных (помимо i-го) узлов-конечного элемента ílt , a í\ = U íit - конечный носитель базисной функции у,.

1 ' i el',

В п.6.2 приводятся особенности программной реализации конечноэле-ментной схемы против потока, состоящие в следующем.

Чтобы вычислить вклады в локальные матрицы конечных элементов при проведении поэлементной сборки матрицы конечноэлементной СЛАУ с использованием схемы против потока вида (49), для каждого узла i необходимо знать номера к, конечных элементов О, из конечного носителя Qt базисной

функции v¡i„ которые обеспечивают выполнение соотношения (50) для члена

ди „ ди ди „

\чх - - и аналогичных соотношении для членов w — и w —. Поэтому сборку дх ду ■ dz

магрицы конечноэлементной СЛАУ удобно выполнять за два прохода по конечным элементам. При этом иа первом проходе для всех узлов сетки можно не

_ , ди ди

только определить по три номерам, (для каждого из членов w,—, wy - - и

ди

w, уравнения (43) может быть свой номер к, конечного элемента, для кото" с:

poro выполняются соотношения типа (50)), но и сформировать и запомнить информацию, необходимую для сборки глобальной матрицы конечноэлементной СЛАУ на втором проходе. Такой информацией являются, например, значения

= jw-tv|/(í/n (51)

для каждого узла конечноэлемепшой сетки (аналогично определяются и £„:).

Очевидно, чго для каждого узла конечноэлемепшой сеткн может существовать несколько конечных элементов с номерами к„ для которых выполняются С001 ношения типа (50). Введем критерии

шах

ки К,

шах

кчГ.

Ь*, шт(-л;)

1*'

В части пат

Ь , па г,

J*'

(52)

(53)

критерия (52) вычисляется наихудший с точки

зрения выполнения неравенств (50) вклад ог члена и1, - в компоненту а„ маг-

дх

рицы конечноэлемептной СЛАУ (для каждого конечного элемент и. Ол из конечного носителя базисной функции при положшелыюм значении ** (см. (51)). Поэтому по критерию (52) определяется номер к* конечного элемента, дяз которого наихудший на одном конечном элементе вклад в компонеиш маг-рицы копсчноэлсменгпой СЛАУ является наилучшим среди всех конечных элемепюв нз конечною носителя С1А- базисной функции »¡¡, при положительном

с.' . То есть для конечного элемента с номером к* выполняется соагношеинс

пн^//' ,'1тпп|- Ь*' ^ = шах шщ|л1, тп>|-/>* ^

(54)

По критерию (53) определяется номер конечною элемента, котсры» обладает темп же свойствами, но при отрицательном значении . Это значит, что для конечного элемента с номером к~ выполняется соотношение

тш^-б,1'ттй'" =,пах тш^--//^, тш/?^

(55)

В критерии (52) и (53) можно ввести и средние значения входящих и неравенства (50) сомножителей:

шах

/

V, тш(-/><)] +а Й;

(56)

л;), 1шпл; ]+а -/),'+ ^ л;

(57)

где а - некоторый вес.

Итак, на первом проходе для каждого (/-го) узла конечноэлементной сетки

вычисляются номера тетраэдров к* и к~ (по критериям (52)-(53) или (56)-(57))

и значения интегралов (5)). При этом во всех узлах сетки формируются по три

ди 8 и ди —. >у, — и н>. — дх Уду г д:

уравнения (43). Далее для каждого / формируется три номера конечных элементов. Первый номер к, определяет номер того конечного элемента, который дол-

ди

жен быть использован для аппроксимации члена и», -■- в схеме против потока.

значения к*, к, и Е,, для каждого из членов , и н1. -- исходного

дх

Этот номер определяется соотношением

д. _ [К , если >0, , если с," < 0.

(58)

Аналогично определяются еще два номера конечных элементов для аппрокси-

ди ди мации IV — ни' -.

> Ву - д:

На втором проходе по конечным элементам при аппроксимации члена

ди . ди ди. п

и' - (аналогично и для и' — и и', —) на конечном элементе Ик проверяется, дх ' ду " &

не совпадает ли сформированное после первого прохода значение к, для каждого узла / элемента Па с номером к этого конечного элемента. Если такое совпадение есть, то к элементам а¡, матрицы А конечноэлементной СЛАУ добавляются вклады Ь' ■ £,', где } - номера всех (включая /-ый) узлов конечного элемента О/с, а значение 4* было вычислено на первом проходе по формуле (51).

Таким образом, вычисленные на первом проходе значения с,' используются дважды: сначала для определения номеров по формуле (58) после окончания первого прохода (но не на втором проходе, а перед его началом), а затем для формирования компонент матрицы конечноэлементной СЛАУ на втором проходе по конечным элементам.

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

ди , ди ди,

ции конвективного члена тг — (и членов н>. — и в>_ — ), а достаточно пройти Х6х ' ду ' д:

по узлам сетки и дособрать глобальную матрицу по коэффициентам />'' и значениям (в этом случае вклады от члена -с11\(Х цгас1г^) следует формировать на первом проходе). Однако такой подход, являющийся довольно экономичным относительно затрат машинного времени на генерацию глобальной матрицы, требует намного больше памяти по сравнению с двухпроходовым вариантом. Выбор того или иного варианта реализации схемы против потока определяется тем, какой из вычислительных ресурсов более дефицитен для решаемой задачи - время или память.

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

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

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

Кроме этих достоинств схема против потока имеет еще одно преимущество, играющее гораздо более важную роль при решении трехмерных задач с большим количеством узлов в конечноэлементной сетке. При использовании стандартного метода Галеркина для аппроксимации дифференциального уравнения (43) с достаточно большими значениями характеристик конвективного теплопереноса и-,, и'г или наиболее эффективные итерационные методы, применяемые для решения конечноэлементных СЛАУ с несимметричными матрицами, могут либо довольно плохо сходиться, либо не сходиться вообще (даже и в тех случаях, когда аппроксимация по методу Галеркина дает физически правдоподобные результаты), и поэтому для решения таких СЛАУ приходится использовать либо довольно, неэффективные итерационные, либо прямые методы. Эю может привести либо к офомным затратам вычислительных ресурсов (памяти компьютера или времени счета), либо даже к невозможности

решения исследуемой задачи из-за нехватки вычислительных ресурсов. Применение же рассмотренной выше схемы против потока позволило без каких-либо существенных затруднений решить все конечноэлементные СЛАУ итерацион-годами, такими как ВСй, ОМЯЕЯ и др. '

Глава 7. Итерационный метод решення СЛАУ с симметричными н несимметричными матрицами, получаемыми в результате конечноэлсмснт-нон аппроксимации краевых задач

В седьмой главе приводится описание итерационного метода решения ко-нечноэлементных СЛАУ с симметричными и несимметричными матрицами.

В п.7.1 приводится схема метода, которая состоит В следующем. Для решения СЛАУ

Ах / (59)

используется итерационный процесс:

л" - произвольное начальное приближение,

/•"=/-Л.т", :"=г\ - (60)

(61)

(62)

гк +а1Лгк', (63)

(64)

Р, =

г 1 - г! I - У. / т 2

(65)

■г^РУ+О-Р.У"1. (66)

Смысл шагов (64)-(66), определяющих выбор направления состоит в следующем. Обозначим через //к гиперплоскость, образованную всевозможными линейными комбинациями векторов гы и гк. Тогда, если у* > 0, то лежит в Нк между лыи -гк, а если у ^ < 0, то гк лежит в Нк между гы и гк. При этом направление 2к фиксируется в Як таким образом, что минимизация вдоль него квадрата нормы невязки (гк+|, гкн) приводит к наименьшему среди всех

значений, получающихся в результате минимизации (г^1, Л4) вдоль любого из векторов геНк.

В п.7.2 дается обоснование предложенной схемы, а в п.7.3 обсуждаются ее вычислительные затраты.

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

Конечноэлементная аппроксимация тепловых полей в задачах без конвективного теплопереноса и электромагнитных полей в режиме постоянного тока или в режиме становления поля для задач электромагнитного зондирования Земли приводит к СЛАУ (59) с симметричной положительно определенной матрицей А. На примерах решений этих СЛАУ проводилось сравнение метода (60)-(66) с МСГ, являющимся одним из наиболее эффективных итерационных методов решения такого рода СЛАУ. Для обоих методов (метода (60)-(66) и МСГ) выполнялась одинаковая процедура предобуеловливапия с использованием неполного разложения Холесского. Критерием остановки итерационных процессов для обоих методов являлось выполнение неравенства

Сравнение метода (60)-(66) с МСГ дало следующие результаты. Для всех задач, по которым проводилось сравнение, количество итераций, требуемых для выполнения условия (67), в методе (60)-(б6) очень незначительно отличалось (причем иногда в большую сторону, но чаще в меньшую) от количества итераций, требуемого для выполнения того же условия в МСГ. Таким образом, метод (60)-(66) при решении конечноэлементных СЛАУ с симметричными матрицами показал сходимость не хуже, чем МСГ, являющийся одним из наиболее эффективных методов решения СЛАУ такого типа.

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

Кроме того, во всех решенных задачах расчета нестационарных электро-

(67)

магнитных полей в технических устройствах и моделирования теплообмена с учетом конвективного геплопереноса метод (60)-(66) позволял получать решения конечно-элементных СЛАУ при хранении всех вещественных чисел с одинарной точностью (в 4-байтовых словах), в то время как при использовании метода ВСО для решения многих из этих задач хранение вещественных чисел с одинарной точностью приводило к его расходимости, и только переход к хранению вещественных чисел с двойной точностью (в 8-байтовых словах) позволял для таких задач добиться сходимости метода ВСО.

Таким образом, метод (60)-(66) при решении задач определенного типа (аппроксимация которых приводила к конечноэлементным СЛАУ с несимметричными матрицами) часто позволял добиваться практически двукратного сокращения требуемого для расчета объема оперативной памяти компьютера. Также немаловажным представляется и тот факт, что метод (60)-(66) оказался более устойчивым к негативным влияниям ошибок округления, так как скорость его сходимости практически не менялась при переходе от хранения вещественных чисел с одинарной точностью к их хранению с двойной точностью. Это свойство метода (60)-(66) может дать ему дополнительные преимущества при решении конечноэлемеитиых СЛАУ очень большой размерности.

Основные результаты работы

Па основе проведенных в диссертационной работе исследований получены следующие теоретические и практические результаты.

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

го и добавочного трехмерного).

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

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

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

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

6. На основе предложенной вариационной формулировки разработаны и реалиюваиы алтритмм конечноэлеменгной аппроксимации нестационарных

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

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

8. Разработана вычислительная процедура, реализуют»« .а« называемую схему против потока при конечноэлементной аппроксимации краевых задач, описывающих процессы теплообмена с учетом конвективного теплоиереноса. На примере решения модельной задачи теплообмена в одном фрагменте технической конструкции показана высокая эффективность реализованной вычислительной схемы.

9. Предложен, обоснован и реализован итерационный метод решения конечноэлементных СЛАУ с разреженными матрицами. Этот метод с успехом применялся при решении различных задач электромагнетизма и теплопереноса, конечноэлеменгная аппроксимация которых приводила к СЛАУ как с симметричными матрицами (например, в задачах электростатики, каротажа, наземно-скважинной и индукционной электроразведки, диффузионного теплообмена), так и с несимметричными матрицами (в задачах теплообмена с учетом конвективного теплопереноса, а также в задачах изучения влияния вихревых токов на нестационарные электромагнитные Коля в технических устройствах).

Основные работы, опубликованные по теме диссертации

1. Глебов H.A., Данилевич Я.Б., Беляев СЛ., Брынский Е.А., Журавлев Г.С., Карчов Д.С., Смолин U.M., Соловейчик Ю.Г., Шурина Э.Г(. Температурные поля ротора турбогенератора со сверхпроводящей обмоткой возбуждения И Изв. АН СССР, Энергетика и транспорт. - 19S2, № 5. - с.79-83.

2. Карчов Д.С., Соловейчик Ю Г, Шурина Э.П. Один подход к решению газовой динамики в винтовом канале // Численные методы механики сплошной

t

среды - 1982, т. 13, № 3,- с. 148-149.

3. Карчов Д.С., Шурина Э.П., Соловейчик ЮГ. О программно-математическом обеспечении вычислительною эксперимента по исследованию тепловых процессов в электрических машинах нового типа // Применение ЭВМ в оптимальном планировании и проектировании. - Новосибирск, НЭТИ, 1982. - с.22-33.

4. Соловейчик Ю.Г. Об алгоритмическом обеспечении вычислительного эксперимента по исследованию движения газа в каналах теплообменников // Применение ЭВМ в оптимальном планировании и проектировании. - Новосибирск, НЭТИ, 1982. - с.8 (-90.

5. Шурина Э.П., Карчов Д.С., Соловейчик Ю.Г. Решение сопряженной задачи для ротора К'П//Исследование крупных электрических машин.-Л., Всесоюзный научно-исследовательский институт электромашиностроения, 1982. - с.24-32.

6. Шурина Э.П., Соломахин В.И., Карчов Д.С., Свидченков М.В., Соловейчик Ю.Г. Организация комплекса программ расчета трехмерных тепловых полей в составных областях // Комплексы программ математической физики. - Новосибирск, ВЦ СО АН СССР, 1982. - с.291-296.

7. Соловейчик Ю.Г. Программный комплекс моделирования турбулентного течения в каналах // Алгоритмическое и программное обеспечение оптимального планирования и проектирования. - Новосибирск, НЭП 1. 1983. - с.45-48.

8. Соловейчик Ю.Г., Шурина Э.П. Численное моделирование теплообмена при турбулентном течении газа в каналах // Алгоритмическое и программное обеспечение оптимального планирования и проектирования. - Новосибирск, НЭТИ, 1983. -с.38-45.

9. Глебов H.A., Дашшевич Я.Б., Журавлев Г.С., Иванов С.А., Карчов Д.С., Соловейчик Ю.Г., Шурина Э.П. Теплообмен в роторе сверхпроводникового турбогенератора// Изв. АН СССР, Энергетика и транспорт. - 1984, № 4. - с. 101-106.

10. Шурина ЭЛ., Соловейчик Ю.Г., Соломахин В.И., Рояк М.Э. Триангуляция сложных областей в задачах моделирования теплового состояния электрических машин // Машинные методы оптимизации, моделирования и планирования эксперимента. - Новосибирск, НЭТИ, 1988. - с.86-88.

Н.Глебов H.A., Данилевич Я.Б., Журавлев Г.С., Захаров В.П., Соловейчик Ю.Г., Тугаев В.А., Шурина Э.П. Проблемы совершенствования термодинамических процессов в роторе сверхпроводнпкового турбогенератора и пути их решения // Изв. АН СССР, Энергетика и транспорт - 1988, № 6,- с.74-79.

12.Кади-Оглы И.А., Соловейчик Ю.Г., Шурина Э.П. Математическое моделирование теплового состояния турбогенераторов с водяным охлаждением // Вопросы проектирования, исследования производства мощных турбо-, гидрогенераторов и крупных электрических машин. - Ленинград, Всесоюзный научно-исследовательский инеттут электромашиностроения, 1988,- с.53-54.

13.Соловейчик Ю.Г. Математические модели трехмерных электромагнитных полей в электротехнических устройствах // Оптимальное проектирование, планирование экспериментов и моделирование многофакторных объектов. -Новосибирск, НЭТИ, 1989. - с. 119-125.

14. Карчов /(.С., Соловейчик Ю.Г.. Васьковский Ю Н. Математическое модели-

>

рованне трехмерного электромагнитного поля с помощью пакета программ РЭМПСО // Техническая электродинамика. - 1990, № 6. - с.32-38.

15. Шурина Э.П., Карчов Д.С., Соловейчик Ю.Г. Моделирование тепловоз состояния трехмерного составного объекта. Модели и алгоритмы // Численные методы и оптимизация. - Таллинн, АН Эстонии, 1990. - с.86-93.

16. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Особенности моделирования нелинейных физических процессов в трехмерных областях // Вопросы атомной иауки и техники. Серия: Математическое моделирование физических процессов. - М., 1992, Вып.З. - с.86-87.

17. Шурина Э П., Соловейчик Ю.Г., Рояк М.Э. Моделирование электромагнитных нолей в трехмерных областях И Вычислительные технологии. - Новосибирск, Институт вычислительных технологии СО РАН, 1993, Т.2, №6. - с.48-53.

18 Моисеев B.C., Соловейчик Ю.Г., Рояк М.Э. Математическое моделирование сложиопостроенных сред // Сборник рефератов №2 Международной геофизической конференции и выставки по разведочной геофизике SEG-EAGO. -М., 1993.-с. 15.

19. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Математическое моделирование физических полей, обусловленных локальными возмещениями // Вычислительные технологии. - Новосибирск, Институт вычислительных технологий СО РАН, 1994, Т.З, №8. - с. 143-147.

20. Моисеев B.C. Соловейчик Ю.Г. Рояк М.Э. Тригубович Г.М. Влияние диэлектрической проницаемости среды на распространение электромагнитной волны // НЕКЛАССИЧЕСКАЯ ГЕОЭЛЕКТРИКА: Материалы международной конференции (28 августа - 1 сентября 1995г., г. Саратов). - Саратов: Нижневолжский НИИ геологии и геофизики, 1995. - с.60-61.

21. Gnidiev A., Royak М., Shurina Е.,.Soloveichik Yu, Tiimov M., Vobly P. Mathematical Simulation of 3-D Magnitostatic Fields Using the Complex of Programs MASTAC // Abstracts of ЛМСА-95.-Novosibirsk: NCC Publisher, 1995.-P. 131-132

22. Royak M., Shurina В., Soloveichik Yu. and Malyshkin V. Parallelization of Computer Code MASTAC Three-Dimensional Finite Elements Method Implementing II Proceedings of PaCT-95 Lecture Notes in Computer Science. - Germany: Springer, 1995. -P.304-313.

23. Slnirina E.P., Soloveichik Yu.G., Royak M.E., Vobly P.D., Grudiev A.V., Tiunov M.A Three-Dimensional Non-iinei!i Magnitostatic Fields Mathematical Modelling // The Scientific Conference on Use of Research Conversion Results in the Siberian Institutions of Higher Education for International Cooperation (S1BCONVERS-95).- Tomsk, 1995. -P.30.

24. Shurina E.P., Soloveitcliik J.G., Royak M.E. Three-dimensional Fields Modeling on Irregular Mesh Using Finite Elements Method II Proceedings of the First Asian Computational Fluid Dynamics Conference 16-19 January. - V.3, Hong Kong, 1995.-P.1225-1226.

25. Моисеев B.C., Соловейчик Ю.Г., Рояк М.Э., Тригубович Г.М. Новое в математическом моделировании электромагнитных полей искусственных источников тока// Международная геофизическая конференция "Электромагшп-

ные исследования с контролируемыми источниками": тезисы докладов. С.-Петербург, 1996. -с.22-23.

26. Моисеев B.C., Соловейчик Ю.Г., Рояк М.Э., Васильев А.В. Математическое обеспечение задач электроразведки для сложи о п острое иных сред // Научно-техническое совещание "Геофизические методы при разведке недр и экологических исследованиях", сборник материалов (докладов).-Томск, 1996 - с.65-66.

27. Рояк М.Э., Соловейчик Ю.Г. Алгоритмы построения нерегулярных треугольных и тетраэдральных сеток // Сб. научных трудов НГТУ. - 1Ioboch-бирск, НГТУ, 1996г., №2(4). - с.39-46.

28. Соловейчик Ю.Г., Рояк М.Э. Математическое моделирование трехмерных нестационарных электромагнитных полей в электротехнических устройствах // Труды третьей международной научно-технической конференции "Актуальные проблемы электронного приборостроения" АПЭП-96 в 11 томах./Т.б Моделирование и вычислительная техника. - Новосибирск, НГТУ, 1996т .-с.66-69

29. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Решение трехмерных нелинейных магнитостатических задач с использованием двух потенциалов. - Новосибирск. 1996. - 28 с. (Препринт / РАН. Сиб. отд-ние. ВЦ ; № 1070).

30. Royak М., Shurina Е., Soloveichik Yu. and Gnrdicv A., Tiunov M., Vobiy P, MAS-TAC - new code for solving three-dimensional non-linear magnetostatic problems // Proceedings of Par-95 (Particle Accelerator Conference and International Conference of High-Energy AcceIerators).-Vol.4, USA, Texas, Dallas: IEEE, 1996.-P.2291-2293.

31. Soloveichik Y.G. Iterative method for solving finite element systems of algebraic equations // Computers Math. Applic. - Vol.33, №6, 1996. - P.87-90.

32. Soloveichik Y.G., Royak M.E. The computing schemes of non-stationary electromagnetic fields FEM modeling in mediums with three-dimensional inhomogeneity // Bulletin of the Novosibirsk Computing Center. Series: Computer Science. - Novosibirsk: NCC Publisher, Issue 4, 1996. - P.55-78.