Метод контрольного объема на неструктурированной сетке для моделирования гидродинамических процессов и распространения радиоволн тема автореферата и диссертации по механике, 01.02.05 ВАК РФ

На правах^^шиси

Фирсов Дмитрий Константинович

Метод контрольного объёма на неструктурированной сетке для моделирования гидродинамических процессов и распространения радиоволн

(Специальность 01.02.05 — механика жидкости, газа и плазмы)

АВТОРЕФЕРАТ

диссертации на соискание ученой степени доктора физико-математических наук

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

Научный консультант: доктор физико-математических наук, профессор,

Бубенчиков Алексей Михайлович, ФГБОУ ВПО «Национальный

исследовательский Томский государственный университет».

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

Алтухов Юрий Александрович, профессор кафедры теоретических основ информатики ФГБОУ ВПО «Алтайская государственная педагогическая академия», г.Барнаул;

доктор физико-математических наук, Моисеенко Сергей Григорьевич, ведущий научный сотрудник отдела физики космической плазмы ФГБУН «Институт космических исследований» Российской академии наук, г.Москва;

доктор физико-математических наук, ст.н.с., Сафронов Александр Иванович, профессор кафедры нанотехнологии,

материаловедение и механика ФГБОУ ВПО «Тольятгинский государственный университет».

Ведущая организация: ФГБОУ ВПО «Кемеровский государственный

университет».

Защита диссертации состоится 14 февраля 2014 г. в 12 часов на заседании диссертационного совета Д 212.029.08 при ФГАОУ ВПО «Волгоградский государственный университет» по адресу 400062, г.Волгоград, пр-т Университетский, 100, ауд. 4-01 А.

С диссертацией можно ознакомится в библиотеке ФГАОУ ВПО «Волгоградский государственный университет».

Автореферат разослан « декабря 2013 г. Ученый секретарь

диссертационного совета Д 212.029.08 доктор физико-математических наук

у Михайлова В.А.

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

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

Набор численных методов для уравнений Максвелла имеет свою специфику, и многие методы, применимые для расчета движения сплошных сред, не применяются для расчета уравнений Максвелла. Более того, многие численные методы для моделирования распространения радиоволн основаны на интегральных уравнениях и позволяют рассчитывать задачи только лишь для отдельной частоты, что не является приемлемым для расчета движения токопроводящих сред в магнитном поле. Именно поэтому возникает необходимость обобщения численных методов, предназначенных для расчета движения жидкостей и газов на задачи о распространении радиоволн. Более того, многие численные методы для моделирования движения радиоволн основаны на полуаналитических соотношениях, не обобщаемых на задачи о движении жидкостей и газов. Одним из важных приложений численных методов к расчету уравнений Максвелла является получение форм электромагнитного поля, оптимально распространяющегося через волноводы определенной форму. Если в 60-х годах двадцатого века для расчета волноводов требовалось использовать полуаналитические численные решения [10], специфичные для каждой из задач, то в настоящее время в связи с прогрессом в численных алгоритмах расчета уравнений в частных производных можно получать численные решения для таких задач на основе более универсальных численных методов. Последние достижения в методах расчетов на нерегулярной сетке позволяют моделировать, в том числе и сложные приборы, такие как антенны сотовых телефонов и оборудование сопряжения с ними [12-17]. К таким примерам относятся работы [10,31,32], в которых приведены полуаналитические решения задач о распространения радиоволн в волноводах. В настоящее время для получения подобных результатов доступны сверхмощные ЭВМ, с использованием которых можно численно найти все собственные значения и вектора аппроксимированных операторов уравнений в частных производных, получая тем самым решения задачи о распространении радиоволн в волноводах и сопряженных

устройствах. В результате подбор волноводов с оптимальными параметрами для проектируемых устройств становится проще и нагляднее.

Примерами таких задач служат резонирующие антенные комплексы, такие как [23 24] для расчета которых часто используются алгоритмы, основанные на вычислении характеристик распределения радиоволн одной характерной частоты [25]. Последние алгоритмы зачастую используют интегральное уравнение, решаемое на поверхности устройства без использования объёмных сеток, в которых для ускорения произведения матрицы аппроксимации интегрального уравнения на вектор используется идея быстрого преобразования Фурье с использованием предобуславливателя на псевдо-Теплицевой матрице [33]. Анализ поведения приборов на основе алгоритмов, рассчитывающих распространение радиоволн для отдельных частот, требует проведения целой серии расчетов для каждой из гармоник. Однако, моделирование распространения радиоволн в нелинейных средах должно учитывать переизлучение энергии на некотором спектре частот, отличном от частот источника излучения, что делает расчёт и анализ таких проблем затруднительным. Использование временной зависимости для расчета компонент электромагнитных полей позволяет моделировать поведение нелинейных сред в широком диапазоне частот и для сигналов различной формы, что открывает возможность моделирования движения токопроводящих сред, модель движения которых описывается системой уравнений Максвелла для переменных, зависящих от времени, совместно с уравнениями Навье-Стокса. Последнее позволяет рассчитывать движение плазмы или электропроводящей жидкости под воздействием электромагнитных полей, что делает возможным моделирование сложных химических реакций, магнетроны, магнитогидродинамические генераторы [37] и прочие современные устройства, используемые для создания

наночастиц или пучков плазмы.

Одним из наиболее часто используемых методов решения уравнении в частных производных в инженерных приложениях, в том числе для уравнений Навье-Стокса, является метод контрольных объемов, важным достоинством которого является выполнение как локальных законов сохранения на каждом из расчетных элементов, так и глобсшьных законов сохранения во всей расчетной области. Выполнение таких законов чрезвычайно важно в задачах гидромеханики. Именно по этой причине, для расчета мультифизических задач о движении токопроводящих сред в электромагнитном поле и не нарушения общности подхода к аппроксимированию систем уравнений в частных производных важно изучить подходы к аппроксимированию уравнении Максвелла методом контрольных объёмов. Необходимо отметить, что применение метода контрольных объёмов к решению уравнений Максвелла не было в достаточной мере изучено и первые широко известные работы [21,22] по этому методу датируются 1999-2000 годами. Последнее связано в первую очередь с высокой диссипагивностью метода контрольных объемов. Тем не

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

Уравнения Максвелла, записанные во временных переменных, подобны по своей структуре уравнениям Навье-Стокса для расчёта движения нестационарных сред [11]. При исследовании проблем связанных с движением жидкости и газа и воздействия их на обтекаемые ими твёрдые тела, рассматривается модель сплошной текучей вязкой среды, которая подчиняется уравнениям Навье-Стокса. Одна из простейших моделей гидродинамики - модель несжимаемой ньютоновской жидкости, которая при отсутствии воздействия внешних сил имеет вид:

р(и, + и-Уы)-и-У2и+ Ун=0,

где р - давление, р - плотность, ц. - коэффициент вязкости, / - время, и -вектор скоростей.

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

Решение уравнений Навье-Стокса методом контрольных объёмов основано на применении скалярных законов сохранения, применяемых при аппроксимации дивергенции. В случае применения уравнений Максвелла дивергентные законы сохранения необходимо обобщить на векторные законы сохранения, используемые для аппроксимации уравнений Максвелла на основе оператора ротора, которые можно записать в виде:

1д,Е-Уу.Н+аЕ=-1, \1д,Н+ЧхЕ=О,

где Е - электрическое поле, Н - магнитное поле, а - проводимость среды, е - электрическая проницаемость среды, ¡я - магнитная проницаемость среды. J - источниковый член. Дифференциальный оператор VхФ=го!{Ф) есть ротор вектора Ф. Также обозначив через В = еЕ вектор электрической индукции, а через В=\кН вектор магнитной индукции приведем закон Гаусса для электрического поля У-2>=4лр -электрический заряд является источником электрической индукции и для магнитного поля V-В=0 - отсутствие магнитных зарядов.

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

расширить аппроксимационый аппарат для метода контрольных объёмов.

Одной из мультифизических задач является моделирование движения токопроводящей жидкости под воздействием магнитного поля_ Моделирование таких процессов важно, например, при расчете М1Д электростанций и тока крови по кровеносным сосудам под воздействием электромагнитного поля. Движение таких жидкостей описываются системой уравнений [38]:

ед,£-УхЯ = /, И,еа,я(+Ух£'=о,

1=о{Е-\1еиХН),

р(а,и+иЛ7и)-И/Ли+У /хЯ,

V •«=(),

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

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

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

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

Применение метода контрольных объёмов к уравнениям Максвелла уже обсуждалось в работах [21,22,26,27], но ни в одной из них не применялись схемы высокого порядка точности [20,22а], которые широко используются для решения сложных задач в гидродинамике, требующих применения схем высокого порядка точности [39,40]. Несмотря на высокую диссипацию характерную для метода контрольных объёмов, применение схем высокого порядка точности, в том числе и существенно неосциллирующих [201 для решения уравнений Максвелла даёт существенный прирост точности численного решения по отношению к аналитическому [22а]. Использованный в работе уровень абстракции записи уравнений Максвелла позволяет применять схемы различного порядка точности и анализировать устойчивость явных схем, основанных на векторных законах сохранения

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

Метод контрольных объёмов позволяет использовать сеточные элементы с плоскими гранями для получения приемлемой точности, вплоть до схем третьего порядка точности. Применение же метода конечных элементов и его вариации, таких как метод спектральных элементов или метода на основе разрывного Галеркина [34,35] требует использования больших многоузловых криволинеиных элементов, позволяющих описывать геометрию расчётных областей с высоким порядок точности. Последнее делает метод менее простым в использовании и приводит к большему количеству вычислений По этой причине использование метода контрольных объёмов имеет больший смысл для задач до 3-го порядка точности включительно. Для достижения более высокой точности на минимальном числе узлов имеет смысл переходить к использованию спектральных базисных функций в методах на

основе разрывного Галеркина [34,35].

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

уравнениями в частных производных и интегральными соотношениями потоков на гранях элементов. Последнее переносит весь вес штгерполяции с центров элементов и численной аппроксимации производных в узловых точках в центрах элементов на интерполяционные соотношения п™ на гранях элементов. Из факта некорректности аппроксимации производных [28] мы фактически используем интегральные соотношения как одну из форм интерпретации аппроксимации производных численного решен* на элементах неструктурированной сетки, накрывающей область решения, которая имеет больший физический смысл, чем обычная аппроксимация Гоизводных. В результате применения соотношения Гаусса-Острограцкото мы можем добиться выполнения законов сохранения, основанных на построении взаимоэквивалентных взаимодействий

элементами. Более того, при применении метода контрольных объемов основные вычисления ведутся на основе осреднённых значении на элементах, "о позволяет источить высокочастотные колебания решения, и как результат стабилизировать аппроксимацию вокруг разрывов решения.

Посколь^равнения Максвелла содержат лишь первые производные как по пространству так и по времени, то для их разрешения обычно применшотся^вные схемы, кДые нуждаются в ном исследоваиии

критерия устойчивости. Критерий устойчивости Для явных схем „а структурированных сетках с использованием анализа фон Неймана [т зачалю может бьтть получен аналитически в термин- -мстР" параметров элементов и в ^-^ГГмГта а'^ бщГплощГвсГх —на неструктурированной

^ке д^ векторных законов сохранения, в общем случае приводит к нГобхоГмостГчисленных вычислений С . Тем не менее, для уравнении

Максвелла, уже были получены неплохие ан^итические оцен^и шага п времени для явных схем [21,22] для которых С = 1 и 2 соответственна в

настоящей работе проведен численный анализ критерия устойчивости, ™Гбьш представлен в работе [22а]. В результате был получен пшенный критерий устойчивости для которого

геометрии расчетных элементов. Более того, для стр^рированнои с™ состоящей из кубических элементов численньш результат ( С-3 ) полностью совпал с результатами применения критерия фон Неймана.

С целью проверки эффективности новой оценки шага по времени и схем высокого порядкаРточности был реализован численный алгоритм, = уравнений Максвелла на неструктурированной сетке как с обычными "/о лошающими граничными условиями, так и с интегр^ьными уравн™ в качестве граничных условий. Алгоритм имел точность аппроксимации вплоть до четвертого порядка точности по времени и пятого по пространству Тем не менее, использование порядка точности выше третьего имеет смысл толь"лишь при применении спектральных базисных функций в численных схемах на основе методов конечных элементов.

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

Полученные численные алгоритмы расчета уравнений Навье-Стокса и уравнении Максвелла были применены к стандартному набору задач гидромеханики и электротехники. Численные результаты продемонстрировали высокую эффективность примененных численных схем

Цсльн) данной работы является исследование обобщения метода контрольных объёмов, со скалярных законов сохранения на векторные примененные к уравнениям Навье-Стокса и Максвелла, аппроксимированные на неструктурированной сетке. А так же гибридизация с другими численными методами.

Для достижения поставленной пели быпм р^ень, гп^,,.,,. ,,г„...

1. Обобщен метод контрольного объёма высокого порядка точности на векторные законы сохранения, записанные на основе уравнений Максвелла.

2. С использованием уравнений Максвелла улучшен критерий устойчивости для явных схем на основе метода контрольного объёма записанного для векторных законов сохранения, на неструктурированной сетке.

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

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

5. Предложен эффективный способ оценки обратного интегрального уравнения с быстро убывающим ядром.

Основные положения выносимые на защиту

1. Обобщен метод контрольного объёма высокого порядка точности на векторные законы сохранения, записанные на основе уравнений Максвелла.

2. Получен новый критерий устойчивости явных схем высокого порядка точности на основе векторных законов сохранения записанных для уравнений Максвелла и проверен на различный задачах с использованием локальных шагов по времени.

3. Сделано обобщение принципа максимума для схем высокого_порядка точности со скалярных законов сохранения на векторные Принцип максимума позволяет использовать критерии устойчивости, полученный для схем первого порядка точности, при исследовании устойчивости схем высокого порядка точности. МЮЧ1ИЯ

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

5 Предложен эффективный способ оценки обратного интегрального

метода контрольных

объй^^^я Максвелла, которые требуют применения векторных объемов на уравнен спожн0й структуры взаимодействия между

эл^трически^^ пшмми юучение таких -рок—ий

по« сформировать обобщенные подходы к аппроксими=ию систем позволяет «лу р н частных производных, которые будут

^йчи ости ¿ш равнений Максвелла на неструктурированно« сетке

для метода контрольных объемов, примененного к Р^н у

несжимаемой жидкости. отключается в

Поакгичесмязначимоехь диссертационнои работы заключается в 1фактическая ^ ^ КОНтоольного объёма высокого порядка для

возможности испГзовГ ~^^ использованием векторных

пазмепов Найден способ существенного ускорения и повышения качества оГнки обрГог" оператора к интегральному уравнению с существенным убьшанием ядра, последний способ дал существенное ускорение и

возрастание качества при восстановлении резкости изображений.

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

Использование результатов работы. Была разработана универсальная среда решения уравнений Максвелла на неструктурированной стеке, которая использовалась в Университете Манитобы для решения различных прикладных задач. Для Томского государственного университета было разработано учебное пособие [17а], в котором приводится детальное описание алгоритма решения уравнений Навье-Стокса методом контрольного объёма.

Апробация работы. Основные результаты работы докладывались на 14 российски и международных конференциях:

1. «Механика летательных аппаратов» Томск 1999 г.,

2. "Вычислительные технологии-2000" Новосибирск 2000 г,

3. Конференция молодых ученых, посвященная 10 летию ИВТ СО РАН Новосибирск 2000 г.

4. Конференция посвященная 80-летию академика H.H. Яненко, Новосибирск 2001.

5. Международная конференция по Математике и Механике 16-18 Сентября 2003., Россия, Томск, Международная летняя школа молодых ученых «Итерационные методы и матричные вычисления» 2-9 Июня 2002, Ростов-на-Дону,

6. XI Международный Симпозиум «Оптика Атмосферы и Океана, Оптика, Атмосферная физика» 23-26 Июня, 2004, Томск.

7. FACE, June 19-20, 2006, University of Victoria, Victoria, ВС, Canada.

8. 23rd Annual Review of Progress in Applied Computational Electromagnetics March 19-23, 2007 - Verona, Italy ACES.

9. CAIMS*SCMAI 2007 The annual meeting of the Canadian Applied and Industrial Mathematics Society, May 20-24 2007, Banff Center, Alberta, Canada.

10. Union Radio Scientific Internationale-URSI, Chicago, Illinois, USA on August 07-16, 2008, March 30 - April 4, 2008.

11. 24th Annual Review of Progress in Applied Computational Electromagnetics - Niagara Falls, 2008, ACES.

12. ACES 2006, Miami, Florida March 12-16 20062, March 30 - April 4,

2008. , . , , A .

13. ICIAM 2011 - 7th International Congress on Industrial and Applied

Mathematics, July 18 - 22,2011, Vancouver, ВС, Canada.

14. SLAM Conference on Computational Science and Engineering, February

28 - March 4, 2011, Grand Sierra Resort and Casino, USA.

Вклят! am-opa. Все результаты, связанные с созданием численных алгоритмов и проведении математического анализа численных алгоритмов были выполнены автором. Кроме того, автором работы была проведена существенная работа по математической формализации физических и инженерных постановок задач, использовавшихся для исследования численных схем. Основное направление инженерных и физических исследований для уравнений Максвелла было сформулировано Джузеппе Ло Ветри Шу Хонг Люй сформулировал список проблем по восстановлению изображений, в результате работы над которыми были получены новые результаты по восстановлению поврежденных изображений. Все численные результаты, полученные в работе, в том числе по восстановлению поврежденных изображений, были получены непосредственно автором.

Содержание работы

в пррппй главе работы изложены основные подходы к

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

объемов

Итак' обозначим через Г,еД элемент с объёмом V из разбиения области £2= U£0Г,, обозначим через It=[j\, j\-,jm,J - множество индексов элементов имеющих общую грань с элементом Т,, грани меаду соседними элементами Г, и Ту , j&I, обозначим через дТу-дТ,п8Ту, нормали к которым обозначим через ft¡,, а площади граней через StJ.

Продемонстрируем основную идею аппроксимации методом контрольных объёмов на простейшей нестационарной задаче аппроксимации уравнения в частных производных £«=/-V« на контрольном объеме Т, с

использованием теоремы Стокса:

±ÍT¡2Mdx=±SdTf.i,u(x)dx. (1)

Введя обозначение:

(2)

и,~Т$т и(х)<к,

мы можем переписать формулу (]) в виде:

ТГ=7;Ьт/'"и(х)ах- (3)

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

С использованием ряда Тейлора осреднённое значение может быть переписано в виде:

- V V 1 дяи[х1) I . »

м'=кЯ{х1~х1,)

'Г,

что позволяет получить для связи между осреднёным значением на элементе Т) и истинным значением в центре массы:

«(*,)=»-ОДА*2). (5)

Иными словами осреднённая величина со вторым порядком точности аппроксимирует значение в центре масс контрольного объема.

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

Воспользовавшись методом контрольных объёмов перепишем уравнения Навье-Стокса с использованием (3) в виде:

Цгг. И-нйК=0, С граничными условиями следующего вида:

и=0 , , и=»ь(х), *еГ1П, 8и18п=0 , *еГоц,, (7)

где Г — граница, представляющая собой твёрдую стенку, Гт -границана входе в канал, либо подвижная стенка, ГоЛ - граница на выходе из канала.

Для каждого слагаемого (6) была применена как Ий1 так и аппроксимация высокого порядка точности [20], предложенная 0™'вер™еТ^вторая заключайся в восстановлении значении прГизводаыГпо значениям в элементах, соседних с элементом в котором восстанавливаются производные ряда Тейлора к-то порядка точности.

^ V _д(щ.....т„\

.....

где ]е8>; = и номер элемента, из шаблона, построенного для

восстановления^производных (5|={/1 ). В результате, для восстановления значений производных мы формируем в общем случае переопределенную систему линейных уравнений

ьди=ь, (9)

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

здесь е=Ю~10 некоторая константа, исключающая деление на ноль, /, есть

масштаб элемента Г,- (расстояние от центра масс до одной из вершин), Я,

НГ.ПМЯ невязки для решения (9) без весовых функций. При аппротеимации уравнений Навье-Стокса

соответствующие производной первого порядка ^ , оператору Лапласа Д Ví^,m ди'и

и конвективной части 1 •

Введем обозначения: ^ д

])1 . матрица аппроксимации производной дх, ■

Ь - матрица аппроксимации оператора Лапласа. (7(у) матрица аппроксимации конвективной части.

В каждой из матриц мы поставим нули в строках, отвечающих граничным элементам размерности меньше dim . В качестве дополнительных матриц, которые завершающих построение дискретного аналога уравнений Навье-Стокса рассмотрим:

I - единичная матрица,

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

/г - матрица, содержащая единицы на диагонали в строках, отвечающих элементам размерности dim~\ , относящимся к участку границы Гсд£2 , и нули во всех остальных строках ( /г=/-/п ).

Рг - матрица, содержащая аппроксимации производной по нормали 8/дп в строках, отвечающих элементам размерности dim-1 , относящимся к участку границы Гсдй , и нули во всех остальных строках. Необходимо отметить, что Рда служит для реализации граничных условий др/дп=О для давления на стенках расчетной области.

Обозначим через 0 вектор нулевых значений размерности N , где N -количество всех элементов разбиения области QUdQ . Для производной по времени будем использовать конечно-разностную схему 1-го порядка точности:

Тогда дискретный аналог (6) с граничными условиями (7) можно представить в виде:

j=l,..,dim, (10)

Для решения полученной системы вида

A" D Dv О

была применена стандартная итерационная процедура SIMPLE [29].

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

произведение в матричном виде:

м Фг^з -ф3 Ьг

= ФзЬ, -ф^З

¿3 фА -ф2Ь,

(П)

О -ф3 Ф2

фХй=С(ф)й= Фз о -ф,

[-Ф2 Ф1 0

что позволяет записать го(Ь=ЧхЬ в виде дивергенции матрицы:

Ухь=с(У)ь=(^с(ь))г.

„„,,„„, г (т1гй=-С(<р)и и применив формулу

Воспользовавшись соотношением МФ; " *

Грина к (12) получим:

¡т^гуС{ч,))тах = !вт с(ф)гйл=/а,с(й)фЛ • (13)

„о л т\ хтяинения Максвелла могут быть записаны В терминах нового оператора (12) уравнения

в виде:

Lд,E-{йvC{H))т+oE=-J . (14)

[ цЭ,Я+(сНуС(£))г=0

„.»г*. ения (14) в компактной операторной

Последнее позволяет записать уравнения в

форме:

где а=

О ц/]' . УХЕ \ \

а 0 С=М и — Е

0 0 ' Ы Я

(15)

(16)

/лл\„, элементе Г и введя обозначение Проинтегрировав выражение (14) на элементе 1, Д*

и применив соотношение (13) перепишем систему (17) в виде.

' (б,/РГ (е)"1 С(й) Я(е)-1 сг ¡г/г,(е)"13 • (1?)

О,

из чего выражение (15) в интегральной форме примет вид:

Здесь под матрицей *(») подразумевается симметрическая матрица:

Ä(Ä)=Ä(Ä)r=

О -C(n, C(A) 0

(19)

Введя обозначение *(й)=а_,Я(л) разложим матрицу ОД на сумму матриц с положительными и отрицательными собственными значениями Ä(»)-Ä(Ä)++Ä(ft)- . Если наше пространство однородно и изотропно, то e—diagje, е, е) и (i=diag{n, и матрица А(п) имеет шесть

собственных значений заданных в виде A=diag{c,c,-c,-c,0,0) , где - скорость света [21]. Матрица R(n) являемся Симметричной, следовательно мы можем её расписать в виде R{n)=Q(n)AQ{h)'1, где

матрица Q(h) состоит из единичных собственных векторов столбцов' ек к = 1..6 для матрицы R(fi) : '

-теа -тгb тгя -kb VtИ 0 (20)

ж« —ягА -^ä 0 VJwj

где ||Ä||=||#f||=l , h=hxb. Следовательно, матрица д(й) может

быть разложена в виде:

л(/»)=е(й)лё(«)-'=е(«)л+ё(»)-,+ё(й)л-ё(я)-,=л(«)чд(А)-> (2i)

Г-е/м->-7А-/+^г' A+ = diag{c,c,0,0,0,0) , A-=diag{0,0,-c,-с,0,0] ,

Q(n) -(aß(rt)) -Также как и в работах [21,20а] элементы разложения (21)

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

-сС(п)2 —е-1 С(п) M."*'C(«) —сС(й)2 _

cC(hf -e-'ciÄ) ,Ц"'С(Й) сС(й)2

(22)

для которых выполняются обычные свойства для противопоточных схем:

Л(Й)-1(Д(Й)+|Л(Й)|). Более того, с использованием внутренней структуры матрицы д(п) (20) мы

можем записать

Л(йГ = сЕ72=1ё7(й)[аё'(й)]г .

В результате вся разница между методом контрольных объёмов для скалярных и векторных величин заключается в том, что притивопоточное разложение записано не для скаляров, а для собственных чисел оператора Л (и) . Более того, здесь следует отметить важную связь между операторами *+(й) и (я) , полученную в [21]:

с С2 -си' ц 0 ' 1 с2 -с Ц"' 0

С/ц с С2 — С 0 -1 с 2 с с2 0 с

О ц

о

-1/2

с2-с с с2

-1/2

о

0

При расчете ™Рала по грани в выражении (18) в"1 *(в)«А мы

применим разложение по потокам:

1вГ1аЧКЛ)«.А=^гА(ЛГ«,+*^"»"А' (24)

ДЛЯ хедТ, мы получаем

В обшем случае значения слева и справа от грани могут различаться то есть

можем обозначить поток через плоскую грань ЗГ, между элементами , Т] в виде:

<26)

радиационные граничные условия [21,20а]:

и,.=О(0Г,с02. (27)

№я достижения высокого порядка ^^^^ пространственных производных, значения ич аппроксимируют

ПОР^мГлатГ[чтоана каждом элементе Т, электромагнитные параметры

Будем полагать, что на кажд отмеЧено выше, метод

с и II являются константами. 1акже, кар.

о^о. г —¿г^т=

обобщенных векторов, состоящих составляющих:

и решается в осреднённых значениях на элементах сетки Д

и,=\г)ти(х)<Ьс

(28)

на каждом элементе ^ м аппроксимацию на основе

С использованием обозначения ^ зал ^ использованием

метода контрольных объемов для уравнении хуы ^ явтой схемы Эйлера для интегрирования по времени в виде [21,20а].

»Г'^^Ла-'^а^«;^«^«;), (29)

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

Т -1 и

О

г,*г,л " О

г г;' rf+r;1

, (зо)

В результате при использовании оператора Ь дискретизации пространственных производных, для решения временной задачи (20) мы можем применять явную схему Эйлера:

и"+1=н"-Д/1,н\ (31)

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

Наряду с простыми граничными условиями (27) на внешней границе расчетной области была использована гибридизация с интегральными уравнениями [23а]. Итак, при применении интегральных уравнений, рассматриваемая область, ограниченная поверхностью ,5 , будет содержать все источники и структуры. Далее рассмотрим некоторую точку г , лежащую вне поверхности Б , например на поверхности расчетной области 2 . Тогда электромагнитное поле в точке г в момент времени Г может быть представлено в виде интегрирования значений поля на поверхности 5 в прошедшем времени и по всем точкам « на поверхности 5 . В результате, последнее выражение может быть записано в виде:

8(<-т

e{t.r)=is-bтг

'.V 4лЯ

(32)

где [их} = [и+тд,и], d,=ô/dt , 0(г) есть функция Хевисайда, R={r-s) , где iî=||r-s||2 есть расстояние между двумя точками. Здесь } обозначает прошедшее время t=t-т, т-RIc и обозначение [/]=/(?,s) используется, чтобы обозначить значения функции / в позиции s , полученной в прошлом времени по отношению к моменту времени t.

Тогда для возвращения энергии в расчетную область через грани элементов, лежащих на поверхности 2 были использованы потоки, входящие в расчетную область [23а]:

v IJ' :j \ ij > ij '

где , a значение и"1' вычислено с использованием (32).

Рис. 1: Гибридизация интегральных уравнений с методом контрольных объёмов.

Гибридизация метода

контрольных объёмов с интегральными уравнениями позволила проводить расчеты на существенно меньших

расчетных областях, чем метод контрольных объёмов, с граничными условиями (27), а также проводить расчеты на изолированных подобластях, между которыми организован обмен данными с

использованием (32), что позволило существенно

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

ыо(г) = /5сехР(-||г-8||2/а)и(8)^. (33)

Другими словами, мы можем заменить оператор

я (г- 8)=с0 (II г -5II) = с ехр (-||г - «Ц2/0) на опеРат°Р вида

тогда йа(г) = /5Я(г-8)«(5)й5. При достаточно большом Я мы можем получить любую точность для аппроксимации йа(г)~и„(г) . Поэтому, в работе

[9а] для восстановления изображений был предложен алгоритм на основе оператора (34):

т

2о)

(35)

который требует 0(М) операций, в то время как алгоритмы на основе Е2Е преобразований Фурье требуют 0(К 1оЕ(Ы)) операций и дают худший результат (Рис. 2).

Рис. 2: Пример восстановления резкости

a) Изображение, полученное не в резкости на цифровой камере.

b) Восстановление алгоритмом (35).

c) Восстановление с использованием быстрого преобразования Фурье в MATLAB.

В четвертой главе описана связь метода контрольных объёмов, записанного на структурированной сетке, с методом конечных разностей. В главе так же представлен способ аппроксимации уравнений Навье-Стокса специальной конечно-разностной схемой, не обобщаемой на метод контрольных объёмов, записанный для неструктурированной сетки.

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

Одним из таких примеров являются конечно-разностные схемы, записанные на разнесенной сетке, в которой для аппроксимации каждой из расчетных величин используются индивидуальные узлы (Рис. 3). Такие сетки не позволяют рассчитывать задачи в неортогональной системе координат и не обладают достаточной гибкостью, хотя сами схемы позволяют быстро получать надёжные результаты. Для аппроксимации уравнений Максвелла представлена компактная запись известной конечно-разностной схемы Йе, позволяющей максимально быстро проводить расчеты на структурированной сетке. Последняя схема очень часто применяется совместно с методом контрольных объёмов для быстрого расчета в трехмерных подобластях с одинаковыми электрофизическими свойствами, которые предварительно заполняют кубическими элементами [21].

Тем не менее, конечно-разностные схемы могут быть использованы для расчёта задач в областях сложной формы. Для этого уравнения в частных производных записывают в неортогональной системе координат и решают на сетках, полученных с использованием эллиптических сеточных генераторов, аппроксимируя значения искомых величин на неразнесенной сетке (в одних и тех же узлах) [6а,8а].

/

Для уравнений Максвелла

Для уравнений Навье^ Стокса

Рис. 3: Разнесенная сетка.

В качестве примера нестандартной аппроксимации разрывов методом конечных разностей приведен пример аппроксимации Гауссовой кривизны [10а] для задач восстановления изображений. Идея аппроксимации заключается в максимизации значений производных или использовании информации о величине разрывов в задачах восстановления изображений.

Пятая глава посвящена новому критерию устойчивости для метода контрольных объёмов, записанному для векторных законов сохранения на примере применения к уравнениям Максвелла. Новый критерий позволил получить тот же результат для схем первого порядка на структурированной сетке, что и критерий устойчивости, основанный на анализе фон Неймана, что на 50% лучше предыдущих результатов. В проведенных численных экспериментах на неструктурированной сетке новый критерий устойчивости обычно давал более чем 10% увеличение шага по времени для схем первого порядка точности. Для схем высокого порядка точности было предложено обобщение принципа максимума со скалярных на векторные законы сохранения, что позволило использовать критерий устойчивости схем первого порядка точности для схем высокого порядка точности. В главе так же приведены схемы высокого порядка точности по времени и их приложение к схемам с локальными шагами по времени с различными

способами согласования расчетных данных.

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

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

V

Д/= min ——, (36)

kissw CjA,

где для элементов с индексом ( с, — это скорость света, V, - объём элемента и At - площадь внешней поверхности элемента соответственно. Однако, уже в 2000 году эта оценка была улучшена в два раза [22]:

2V:

А /= min —7. (37)

к«?/ c,Aj

Тем не менее, шаг по времени (37) не является необходимым, поэтому возможно получить оценку для большего шага по времени. Так, например, на структурированной сетке, состоящей из кубических элементов с одинаковыми диэлектрическими свойствами, оценка (37) дает: At=h/(3c) что в 1.5 раза меньше, чем получаемая из анализа фон Неймана [21]:

M=0.5hic. (38)

Именно по этой причине критерий (37) не является необходимым. Целью главы является получение нового критерия сходимости, который как минимум является необходимым на сетках с одинаковыми элементами. Здесь стоит отметить, что доказательство того факта, что критерий является необходимым на неструктурированной сетке является не тривиальным. Также, как и в доказательстве, приведенном в работе [22] вывод базируется на использовании естественных физических ограничений, заключающихся в том, что на сетке, покрывающей расчетную область и не имеющую источников, при переходе с одного шага по времени на другой энергия не нарастает. Энергия электромагнитного поля в области ficR3 , покрытой сеткой А , состоящей из элементов Г, с объёмами К, определяет норму следующим образом: \\u\\=<J5п£Е-Е+\1 Н-Н cbc, дискретный аналог которой можно определить по формуле:

11«11(39)

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

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

Новый критерий получается численно с использованием геометрии элементов сетки, и дает больший шаг по времени, чем результаты опубликованные в предыдущих работах [21,22]. Более того на структурированной сетке критерий дает тот же самый результат, что и критерий Фон Ноймана [21] на структурированной сетке. Другими словами, новая оценка позволяет соединить условия устойчивости на структурированной и не структурированной сетках, что позволяет сократить вычислительные затраты при расчете шага по времени.

Для обобщения численного метода оценки шага по времени на схемы высокого порядка точности по пространству приводится способ оценки максимального наклона градиентов в схеме МШСЬ. Предложенный способ оценки наклона градиентов позволил сформулировать принцип максимума [41] для метода контрольных объёмов, примененного к векторным законам сохранения.

Итак, для схем высокого порядка точности по пространству, применённых к векторным законам сохранения [22а] и схеме Эйлера, имеющей первый порядок точности по времени новый критерий устойчивости будет иметь вид:

нормализированный вектор, полученный из одного из рёбер на грани дТ,л , а коэффициент 1<&,<2 будет зависеть от ограничения на вектора потоки размерности 6 между соседними элементами, записанными для схем высокого порядка точности:

Применение принципа максимума и ограничения на величины потоков между соседними элементами позволяет гарантировать устойчивость схемы в районе расположения разрывов. Использование рассчитанного или введенного ограничения 1<б^<2 на величину потоков между соседними элементами позволяет ввести критерий сходимости для схем высокого

порядка точности Д/*=Д /,/&, при котором величина шага по времени хоть и снижается но схема практически не теряет точности и остаётся стабильной. В случае же применения схем высокого порядка точности к областям решения

Д 1<тт

<?А1|С,Г

(40)

не имеющим разрывов сохраняется критерий сходимости для схем первого порядка точности. Поэтому коэффициент б, рекомендуется брать зависящим от задачи. В задачах в которых решение не имеет разрывов рекомендуется брать д, = 1 . В задачах же со значительным числом разрывов было достаточно ввести ограничение 6, = 1.4 . В случае же отказа от ограничений на разрывы рекомендуется вычислять локальные значения 8,, либо брать шаг по времени в 2 раза меньший чем для схемы первого порядка точности. С целью тестирования эффекта вычислимого шага по времени был проведен ряд численных экспериментов с использованием программы, описанной в работе [20а]. Эти вычисления были проведены для обеих неструктурированной и структурированной сеток при 6=1 . Для неструктурированной сетки было получено 5%-15% увеличение шага по времени по сравнению с критерием, вычисленным по формуле (37) из работы [22].

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

В локальных шагах по времени, нам требуется выполнять несколько промежуточных шагов для достижения уровня согласования с соседними шагами по времени. Последнее накладывает требование расчетов предиктора в промежуточных шагах с использованием значений из соседних областей [30] и как следствие расширение подобластей £2- и £2- на соседние области с другими шагами по времени чем Д г,-.

8 8 8

8

7 —

5 1 6 £ ...........................

4 -!--2Щ--- ► /4 5

4

3 3

1 2

£2®, ЛГ, = А( £2°,Д/2 = 1.5 А1 £2°,ДГ3=2ДГ £2°,Д/4=ЗДГ

Рис. 4: Пример расчетной области для локальных шагов {д?,:Д/тп1=Д/(Л',, 1 </<;¥,) повремени.

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

шагов по времени (д?, : Л/та1=Д/,Лг„Аг,е1Ы,1</<А?,! , которые задают разбиение расчетной области £2= и*: ,£2? и согласуются одновременно при достижении Мтах. При этом V/ на всех элементах подобласти £2° мы имеем один шаг по времени Д/(; подобласть £2,' , соседняя с £2, и набранная из элементов соседних областей служит для расчета входящих потоков в область . Подобласть £2* соседняя с £2' , набранная из элементов соседних областей, служит для расчета входящих потоков в область £2,'. Пример такой расчётной области представлен на рис 4, а алгоритм для получения решения на шаге по времени 1п* -с"=А1тт

выглядит следующим образом: о „ _ п

Шаг 0: Зададим на всей вычислительной области £2=и,=1£2, , и=н", Г,-( ,

а затем рассчитаем Ьи .

Шаг 1: Для /т,„=тш,(/,+ДО сформируем множество индексов

Шаг2: V/бзтт рассчитаем предиктор м-и-Л/^ы на множестве

элементов из области £2°, на элементах из множества Гуе£2,'иЦ

воспользуемся формулой й^и —(ti+hti—tt^J})Luj.

Шаг 3: Д/,.

Шаг 4: На Ц^ £2"и£2,' рассчитаем Ьй .

Шаг 5: V ¡&5т,„ на перерассчитаем и с использованием формулы для корректора и=4(и+й-Д*,/.й), а на элементах из множества £2,°и£2'\ и1&5 £2? воспользуемся формулой

—(Г,— Г/1у))Ьну) •

Шаг 6: Если V/ и и"+1 =н , то алгоритм завершается, иначе на

цЕ5 £2° рассчитаем Ьи , а затем возвратимся на Шаг 1.

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

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

Течение жидкости под воздействием постоянного магнитного поля

В качестве примера расчета сопряженной мультифизической задачи о движении жидкости под воздействие магнитного поля приводится результат моделирования стабилизированного течения токопроводящей жидкости в

трубке под воздействием стационарного магнитного поля. Движение такой жидкости описывается системой уравнений:

Д ф=У-(иХВ), Е=-Чц,1=о{Е-иХВ), д,и+ м-Ум—Яе'-Д и+ V р=На21X В, V и=О,

(42)

в которой под На=В0 И{а1\у)0'' понимается число Гартмана, Ке=р и К/и -характерное для рассматриваемой задачи число Рейнольдса, II - характерная скорость течения, Я - характерный размер, (.1 — коэффициент кинематической вязкости.

Уравнения (42) решались методом установления для потенциала магнитного поля и движения жидкости раздельно, с использованием метода конечных разностей как в ортогональной, цилиндрической системе координат, совмещенной с цилиндрической трубкой так и в неортогональной системе координат. В случае применения неортогональной системы координат использовалась конечно-разностная схема, описанная в Главе 4 диссертационной работы. Для решения уравнений в ортогональной системе координат использовался конечно-разностный, упрощенный метод приближения узкого канала, детально описанный в работах [18,5а,7а].

Рис. 5: Эффект Гартмана для контуров тока. Распределения отвечают участку стабилизированного течения.

Последние результаты для матричной противопоточной схемы для

На=5

О

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

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

В связи с этим была проведена серия расчетов без учета влияния индуцированных электрических полей (с ср=0 ). Как показали вычисления, по крайней мере до значений На = 5 при описании течении электропроводящей жидкости в узких трубках с непроводящими стенками мы можем пользоваться осесимметричной постановкой задачи, опирающейся на следующую систему уравнений:

1дих дих 8и\ др I д I ди\

Яе —-+и

\д<

д(ихг) [ д[игг) дх дг

— и,На,

(43)

в которой распределения их и иг уже не зависят от 0,а рх-рАх->) ■

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

X |£еаиетв!1 Ри*& а! ВаЛ-ЭсаНегикайо! [0.015303,0.028997. -в.8339)0гк1 3«е: 0.75т

■Ага1уИс I Ро|упопУа| Г.

Средняя длинна грани

Количество Размер ячеек требуемой памяти

Количество

шагов по

Аппроксимация

¿2 ошибка Е,

времени

элемента

мша.

мша.

7.40%

мша

Полиномиальная-3

Полиномиальная-3

Полиномиальная-3

Конечные Разности

Заключение

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

Так, в работе метод контрольных объёмов был представлен на основе аппроксимации уравнений Навье-Стокса, описывающих движение несжимаемой жидкости. Данная аппроксимация приводит к возникновению неявной системы уравнений, которая требует специальных итерационных методов для получения надежных и устойчивых решений. Необходимо отметить, что явные схемы не позволяют получать надежные результаты для уравнений Навье-Стокса для несжимаемой жидкости. С использованием опыта аппроксимации уравнений Навье-Стокса методом контрольных объёмов была построена явная схема на основе метода контрольных объёмов для уравнений Максвелла, которая была обобщена на схемы высокого порядка точности. Полученная схема дала возможность существенно повысить качество расчетов уравнений Максвелла методом контрольных объёмов и позволила по качеству полученных результатов и простоте реализации на неструктурированной сетке конкурировать с методом разрывного Галёркина. Более того, решение уравнений Максвелла методом контрольных объёмов на неструктурированной сетке позволяет сформулировать надежный подход к расчету систем уравнений в частных производных на основе векторных законов сохранения и обобщить полученный опыт на решение многофизических задач, таких как движение проводящих жидкостей в магнитном поле.

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

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

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

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

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

Литература

1 Козлов А.Н. Кинетика ионизации и рекомбинации в канале плазменного ускорителя. // Известия РАН, механика жидкости и газа. 2000, № 5, с. 181-188.

2 Козлов А.Н. Модель пристеночной проводимости в окрестности макронеоднородной зеркально отражающей поверхности. // Физика плазмы,

2002, т. 28, №2, с. 180-187.

3 Козлов А.Н. Влияние продольного магнитного поля на эффект Холла в канале плазменного ускорителя. // Известия РАН, механика жидкости и газа,

2003, №4, с. 165-175.

4 Козлов А.Н. Динамика вращающихся потоков в канале плазменного ускорителя с продольным магнитным полем. // Физика плазмы, 2006, т. 32, № 5, с. 413-422.

5 Kozlov A.N., Zaborov A.M. Formation of the current attachments in plasma accelerator channel under influence of the longitudinal magnetic field. // Problems of Atomic Science and Technology. Series: Plasma Physics, 2006, N 6, p. 138-140.

6 Морозов A.M., Козлов А.Н. Эффект самоочищения потока водородной плазмы в ускорителе КСПУ. II В сб. "Физика экстремальных состояний вещества - 2007" под ред. Фортова В.Е. и др. Изд. ИПХФ РАН, Черноголовка, 2007, с. 316-319.

7 Козлов А.Н. Воздействие продольного магнитного поля на компрессионные потоки плазмы. // Препринт ИПМ им. М.В. Келдыша РАН. 2007, № 87.

8 Kozlov A.N. Basis of the quasi-steady plasma accelerator theory in the presence of a longitudinal magnetic field. // J. Plasma Physics, 2008, vol. 74, part 2, p. 261-286.

9 Друкаренко СЛ., Климов H.C., Козлов A.H., Москачева А.А., Подковыров B.JJ. Приэлектродные процессы в квазистационарном плазменном ускорителе (КСПУ) с продольным магнитным полем. // В сб. "Физика экстремальных состояний вещества - 2008" под ред. Фортова В.Е. и др. Изд. ИПХФ РАН, Черноголовка, 2008, с. 262-265.

10 А. С. Ильинский, А. Г. Свешников, "Методы исследования нерегулярных волноводов", Ж. вычисл. матем. и матем. Физ., 8:2 (1968), 363373.

11 К. Флетчер Вычислительные методы в динамике жидкостей // М.: Мир, 1991, т.1,2.

12 Григорьев А.Д., Салимое Р.В. Моделирование волновых

электромагнитных полей методом векторных конечных элементов // Известия высших учебных заведений. Радиоэлектроника. — 2005. Вып. 1.- С. 82-85.

13 Салимое Р. В. Встраивание сосредоточенных элементов и портов в векторный метод конечных элементов // Известия высших учебных заведений. Радиоэлектроника. 2008. Вып. 2. С. 66-69.

14 Григорьев А. Д., Салимое Р. В., Тихонов Р. И. Сравнительный анализ результатов моделирования антенны сотового телефона различнымипрограммными средствами. // Труды конференции, посвященной Дню радио, Изд-во СПбГЭТУ, 2008, с. 28.

15 Григорьев А, Д., Салимое Р. В., Тихонов Р. И. Расчет поля и параметров антенн сотовых телефонов. // Материалы 8-й международной конференции "Актуальные проблемы электронного приборостроения", Саратов: Изд-во СГТУ, 2008, с. 243 249.

16 Григорьев А. Д., Салимое Р. В., Тихонов Р. И. Сравнительный анализ результатов моделирования антенн сотовых телефонов.//. Материалы научн.-техн. семинара "Инновационные разработки в СВЧ технике и электронике". СПб.: Изд-во СПбГЭТУ, с. 16-17, 2008.

17 Григорьев А. Д., Салимое Р. В., Тихонов Р. И. Моделировние антенн сотовых телефонов методом векторных конечных элементов// Радиотехника и электроника, том 57, № 3, Март 2012, С. 261-270.

18 А. М. Бубенчиков, А. В. Клевцова, С. Н. Харламов, "Закрученный поток проводящей жидкости в узких трубах при наличии магнитного поля", Матем. Моделирование, Том 16, № 3 (2004), 109-122.

19 J. G. Charney, R. Fjortoft, J. von Neumann, "Numerical Integration of the Barotropic Vorticity Equation", Tellus 2, (1950), pp 237-254.

20 Carl F. Ollivier-Gooch Quasi-ENO Schemes for Unstructured Meshes Based on Unlimited Data-Dependent Least-Squares Reconstruction, Journal of Computational Physics, Vol. 133 (1), 1997, pp 6-17.

21 P. Bonnet, X. Ferrieres, PL. Michielsen, P. Klotz, "Finite Volume Time Domain Method," in Time Domain Electromagnetics, S.M. Rao (editor), Academic Press, San Diego, 1999.

22 S. Piperno, "L2-Stability of the upwind first order finite volume scheme for the Maxwell equations in two and three dimensions on arbitrary unstructured meshes", Mathematical Modeling and Numerical Analysis,Vol. 34, No. 1, pp. 139158,2000.

23 F. Edelvik, G. Ledfelt, P, Lotstedt, D. J. Riley, "An Unconditionally Stable Subcell Model for Arbitrarily Oriented Thin Wires in the FETD Method," // IEEE Trans. Antennas Propagat., Vol. EMC-51, No 8, pp. 1797-1805, Aug. 2003.

24 F. Edelvik, "A New Technique for Accurate and Stable Modeling of Arbitrarily Oriented Thin Wires in the FDTD Method," // IEEE Trans. Electromagn. Compat., Vol. EMC-45, No. 2, pp. 416-423, May 2003.

25 Gibson, Walton C. "The Method of Moments in Electromagnetics," Chapman & Hall/CRC, 2008.

26 Christophe Fumeaux, Dirk Baumann, Pascal Leuchtmann, Riidiger Vahldieck, "A generalized local time-step scheme for the FVTD method for efficient simulation of microwave antennas" // 33rd European Microwave Conference, 2003.

27 Christophe Fumeaux, Dirk Baumann, Pascal Leuchtmann and Riidiger Vahldieck "A Generalized Local Time-Step Scheme for Efficient FVTD Simulations in Strongly Inhomogeneous Meshes" // IEEE Transactions on Microwave Theory and Techniques, Vol. 52, No. 3, March 2004.

28 Самарский A.A. Теория разностных схем. Москва. «Наука». 1989г. 616 с.

29 С. Патанкар Численные методы решения задач теплообмена и динамики жидкости // М. Энергоагомиздат, 1984.

30 Hua-zhong Tang, Gerald Warnecke "High resolution schemes for conservation laws and convection-diffusion equations with vaiyinc time and space grids" // Journal of Computational Mathematics, Vol.24, No.2, 2006, 121-140.

31 Э. А. Хорошееа, Ю. Г. Смирнов Распространение электромагнитных ТМ-волн в круглых диэлектрических волноводах, заполненных нелинейной средой // Известия высших учебных заведений. Поволжский регион. Естественные науки. — 2006. — № 5. — С. 106—114.

32 Э. А. Хорошееа, М. Ю. Медведик, Ю. Г. Смирнов Численное решение задачи о распространении электромагнитных ТМ-волн в круглых диэлектрических волноводах, заполненных нелинейной средой // Известия высших учебных заведений. Поволжский регион. Физико-математические науки.-20Ю.-№ 1.-С. 2-13.

33 Akhmed Ould Khaoua Toepliz preconditioner for pseudo toeplitz matrices // Mecanica Computationa Volume, Bahia Blanca, Argentina November 2003, pp 1421-1429.

34 Guido Kanschat Discontinuous Galerkin Methods for Viscous Incompressible Flow. Deutscher Universitats-Verlag | GWV Fachverlage GmbH, Wiesbaden 2007,183p.

35 Serge Piperno Symplectic local time-stepping in non-dissipative DGTG methods applied to wave propagation problems // ESAIM: Mathematical Modelling and Numerical Analysis Vol. 40, No 5, 2006, pp. 815-841.

36 М. Е. Hubbard Multidimensional Slope Limiters for MUSCL-Type Finite Volume Schemes on Unstructured Grids // Journal of Computational Physics, 1999, V.155, P.54-74.

37 D.V. Gaitonde Tip-to-tail scramjet simulation with plasma-assisted control // Computer Architecture, 2004. Proceedings. 31st Annual International Symposium on 7-11 June 2004, p. 199 - 205.

38 Пейре Роже, Тейлор Томас Д. Вычислительные методы в задачах механики жидкости. Ленинград. Гидрометеоиздат 1986 г.

39 Z. J. Zhang and С. Р. Т. Groth, Parallel High-Order Anisotropic Block-Based Adaptive Mesh Refinement Finite-Volume Scheme, AIAA Paper 2011-3695, June 2011.

40 S. D. McDonald, M. R J. Charest, and С. P. T. Groth, High-Order CENO Finite-Volume Schemes for Multi-Block Unstructured Mesh, AIAA Paper 20113854, June 2011.

41 T. Barth, M. Ohlberger, Finite volume methods: foundation and analysis, in: E. Stein, R. de Borst, T.J.R. Hughes (Eds.), Encyclopedia of Computational Mechanics, Wiley, 2004, pp. 1—57.

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

la Бубенчиков A.M., Фирсов Д.К. Численное исследование вихревых структур в прямоугольной каверне. Сборник статей «Вычислительная гидродинамика». Томск. 1999г. с.8-14.

2а Фирсов Д.К. Расчет течения несжимаемой жидкости на неортогональной неразнесенной сетке. Труды конференции молодых ученых, посвященной 10-летию ИВТ СО РАН. 2000, Т. II. с. 156-160.

За Фирсов Д.К. Конечно-разностный алгоритм расчета течений несжимаемой вязкой жидкости, основанный на алгебраическом разложении скоростей. Труды международной конференции, посвященной 80-летию академика Н.НЛненко. Новосибирск, Академгородок, 24 - 29 июня 2001 года. с.65 8-665.

4а Бубенчиков A.M., Фирсов Д.К, Алъбрандт Е.В. Гемодинамика крупных кровеносных сосудов с аневризмой // Вестник Томского гос. ун-та. Бюл. опер. науч. инф. 2001. № 4. Численные методы в динамике вязкой жидкости. С. 23-31.

5а Бубенчиков A.M., Кчеецова А.В., Фирсов Д.К. Приближение «узкого канала» для течения токопроводящей жидкости // Вестник Томского гос. унта. Бюл. опер. науч. инф. 2001. № 4. Численные методы в динамике вязкой жидкости. С. 53-63.

6а Бубенчиков A.M., Фирсов Д.К. Разностная схема для интегрирования уравнений Навье - Стокса несжимаемой вязкой жидкости на неразнесенной неоргогональной сетке // Вестник Томского гос. ун-та. Бюл. опер. науч. инф. 2001. № 4. Численные методы в динамике вязкой жидкости. С. 5-22.

7а А. М. Бубенчиков, А. В. Кпевцова, Д. К. Фирсов, "Течение проводящей жидкости в тонких трубках в поперечном магнитном поле", Матем. моделирование, Том 15, № 9 (2003), 75-87.

8а А. М. Бубенчиков, Д. К. Фирсов, М. А. Котовщикова, "Ламинарное вязкое течение в искривленных каналах сложного сечения", Матем. моделирование, Том 16 №11, 2004, с.107-119.

9а D.K. Firsov , S.H. Lui "A fast deblurring algorithm"//Applied Mathematics and Computation 183(1): pp 285-291 (2006).

10a D.K. Firsov , S.H. Lui "Domain decomposition methods in image processing using Gaussian curvature" // J. Comput. Appl. Math. 193 (2006), No. 2, pp. 460-473.

11a Dmitry K. Firsov, Ian Jeffrey, Joe LoVetri and Colin Gilmore "An Object Oriented Finite-Volume Time-Domain Computational Engine", ACES 2006, Miami, Florida March 12-16 2006.

12a Dmitry Firsov, Joe LoVetri, Ian Jeffrey, Vladimir Okhmatovskiy, and Walid Chamma "High-Order FVTD on Unstructured Grids", ACES 2006, Miami, Florida March 12-16 2006.

13a Dmitry К Firsovi Joe LoVetri and Vladimir Okhmatovski "FVTD - Integral Equation Hybrid for Maxwell's Equations", FACE, June 19 - 20, 2006, University of Victoria, Victoria, ВС, Canada.

14a S. Yu. Sadov and P. N. Shivakumar and D. Firsov and S. H. Lui and R. Thulasiram, Mathematical Model of Ice Melting on Transmission Line, J. Math. Model. Algorithms, 2007, c. 273-286.

15a Ian Jeffrey, Dmitry K. Firsov, Colin Gilmore, Vladimir Okhmatovski and Joe LoVetri "Parallel High-Order EM-FVTD on an Unstructured Mesh" 23rd Annual Review of Progress in Applied Computational Electromagnetics March 1923, 2007 - Verona, Italy 2007 ACES, 401-408.

16a Dmitry K. Firsov, Joe LoVetri "Necessary Stability Criterion for Unstructured Mesh Upwinding FVTD Schemes for Maxwell's Equations" 23rd Annual Review of Progress in Applied Computational Electromagnetics March 1923,2007 - Verona, Italy 2007 ACES, pp. 1705-1711.

17a Д. К. Фирсов «Метод контрольного объёма на неструктурированной сетке в вычислительной механике», Учебное пособие, Томский Государственный Университет, 2007 г.

18а Бубенчиков A.M., Фирсов Д.К., Котовщикова М.А. Численное решение плоских задач динамики вязкой жидкости методом контрольных объемов на треугольных сетках // Математическое Моделирование. 2007, Том 19, №6, с.

19а А. М. Бубенчиков, В. С. Попонин, Д. К. Фирсов, "Спектральный метод решения плоских краевых задач на неструктурированной сетке", // Математическое моделирование, 2007, т. 19 № 10, с. 3-14.

20а D. Firsov, J. LoVetri, I. Jeffrey, V. Okhmatovski, C. GИтоге, W. ChammaJ"High-Order FVTD on Unstructured Grids using an Object-Oriented Computational Engine", ACES Journal vol. 22, no.l, pp. 71-82, March 2007.

21a Cameron J. Kaye, Dmitry Firsov, АН M. Mehrabani, Lotfollah Shafai and Joe LoVetri, "Design of Low-Profile Wire Antennas Over an AMC Using FVTD», ACES 2008, Niagra Falls 2008.

22a D. K. Firsov and J. LoVetri, "New Stability Criterion for Unstructured Mesh Upwinding FVTD Schemes for Maxwell's Equations," ACES Journal, Vol. 23, No. 3, pp. 193-199, September 2008.

23a D. K. Firsov and J. LoVetri, "FVTD—Integral Equation Hybrid for Maxwell's Equations," Int. J, of Numerical Modelling: Electronic Networks, Devices and Fields (Special Issue on Frontiers of Applied Computational Electromagnetics), vol. 21, no. 1-2, pp. 29-42, January-April 2008.

24a C. Kaye, C. Gilmore, P. Mojabi, D. Firsov, J. LoVetri Development of a Resonant Chamber Microwave Tomography System //Ultra-Wideband, Short Pulse Electromagnetics 9, ISBN 978-0-387-77844-0. Springer Science+BusinessMedia, LLC, 2010, p. 481-488.

71-85.

Подписано в печать 06.12 2013 г. Формат 60 * 84/16. Бумага офсетная. Гарнитура Times. Усл. печ. л. 2,4. Тираж 100 экз. Заказ 1274.

Отпечатано с готового оригинал-макета в ООО «Волгоградское научное издательство» 400011, г. Волгоград, ул. Электролесовская, 55.

Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «Национальный исследовательский Томский государственный университет»

Метод контрольного объёма на неструктурированной сетке для моделирования гидродинамических процессов и распространения радиоволн

(Специальность 01.02.05 — механика жидкости, газа и плазмы)

Диссертация на соискание ученой степени доктора физико-математических наук

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

иэ^и'1 <гЭиЭ61

Фирсов Дмитрий Константинович

Научный консультант:

Профессор, д.ф.-м.н. А. М. Бубенчиков

Томск - 2013

Оглавление

ВВЕДЕНИЕ.......................................................................................................................................5

ГЛАВА 1. МЕТОД КОНТРОЛЬНЫХ ОБЪЕМОВ НА НЕСТРУКТУРИРОВАННОЙ СЕТКЕ..............................................................................................................................................17

§1.1. Преимущества метода контрольных объемов 18

§ 1.2. Неструктурированная сетка и сеточные генераторы 19

п. 1.1.1. Описание элементов, вычисление площадей граней и объёмов...........................22

п. 1.1.2. Представление сетки в теле программы..............................................................24

§ 1.3. Форма конечного объёма и осреднённые значения 26

§ 1.4. Обобщение метода контрольных объёмов на суперэлементные методы высокого порядка 30

§ 1.5. Параллельные вычисления уравнений в частных производных 34

п. 1.5.1. Алгоритм параллелизации на ЭВМ с распределенной памятью..........................37

п. 1.5.2. Алгоритм параллелизации на многопроцессорной ЭВМ.......................................39

ГЛАВА 2. ПРИМЕНЕНИЕ МЕТОДА КОНТРОЛЬНЫХ ОБЪЁМОВ К УРАВНЕНИЯМ

НАВЬЕ-СТОКСА ДЛЯ НЕСЖИМАЕМОЙ ЖИДКОСТИ.....................................................42

§ 2.1. Аппроксимация операторов градиента и дивергенции 45

§ 2.2. Локальная интерполяция MUSCL второго порядка точности 47

§ 2.3. Аппроксимация диффузии 49

§ 2.4. Аппроксимация конвекции 51

§ 2.5. Аппроксимации высокого порядка точности 53

п. 2.5.1. Построение шаблона схемы высокого порядка точности...................................55

п. 2.5.2. Выбор весов для элементов интерполяционного полинома..................................58

п. 2.5.3. Использование интерполяции высокого порядка точности.................................60

§ 2.6. Построение и решение системы линейных алгебраических уравнений 61

п. 2.6.1. Построение СЛАУ из аппроксимации уравнения Пуассона.................................61

п. 2.6.2. Матрицы аппроксимации дифференциальных операторов в уравнениях Навъе-

Стокса.....................................................................................................................................63

§ 2.7. Решение полученной системы алгебраических уравнений 64

п. 2.7.1. Матрица системы уравнений..................................................................................64

п. 2.7.2. Итерационный метод решения системы уравнений (2.55).................................66

п. 2.7.3. Процедура SIMPLE...................................................................................................69

п. 2.7.4. Понятие разреженной матрицы............................................................................72

п. 2.7.5. Итерационные методы для решения систем линейных уравнений.....................73

ГЛАВА 3. МЕТОД КОНТРОЛЬНЫХ ОБЪЕМОВ ДЛЯ УРАВНЕНИЙ МАКСВЕЛЛ А.... 76

§ 3.1. Метод контрольного объема для векторных законов сохранения 77

§ 3.2. Законы сохранения в теории электромагнетизма 78

п. 3.2.1. Грань между элементами имеющими разные электрофизические свойства.... 8 2

п. 3.2.2. Случай грани на поверхности идеального проводника..........................................83

§ 3.3. Отраженное электромагнитное поле 84

п. 3.3.1. Отраженное электромагнитное поле в диэлектриках........................................84

п. 3.3.2. Отраженное поле на поверхности совершенного проводника............................85

§ 3.4. Альтернативная форма записи уравнений Максвелла 86

п. 3.4.1. Отраженное магнитное поле.................................................................................88

§ 3.5. Моделирование источников электромагнитного излучения 89

п. 3.5.1. Моделирование тонких проводов для метода контрольных объёмов.................90

п. 3.5.2. Противопоточная аппроксимация телеграфных уравнений...............................93

п. 3.5.3. Оконечные колебательные контуры на концах телеграфных уравнений...........97

§ 3.6. Граничные условия для уравнений Максвелла в задачах электромагнетизма 100

п. 3.6.1. Поглощающие граничные условия.........................................................................100

п. 3.6.2. Граничные условия Лиао [70,71]...........................................................................101

Экстраполяционная схема Лиао.......................................................................................................................102

Обобщённая интерполяционная формула для граничных условий...............................................................104

п. 3.6.3. Граничные условия, основанные на интегральных уравнениях..........................105

Быстрый расчет интегральных уравнений......................................................................................................109

ГЛАВА 4. СВЯЗЬ МЕТОДА КОНТРОЛЬНЫХ ОБЪЁМОВ И МЕТОДА КОНЕЧНЫХ РАЗНОСТЕЙ..................................................................................................................................114

§4.1. Некоторые сведения из теории разностных схем 114

§ 4.2. Взаимосвязь метода конечных разностей и метода контрольных объёмов 117

§ 4.3. Конечно-разностная аппроксимация уравнений Навье - Стокса 120

п. 4.3.1. Система определяющих уравнений.......................................................................120

п. 4.3.2. Неявная конечно-разностная аппроксимация уравнений движения.................122

п. 4.3.3. Обратимость матрицы 2.....................................................................................123

п. 4.3.4. Дискретизация конвективной части....................................................................127

§ 4.4. Новый способ дискретизации уравнений Навье - Стокса, основанный на алгебраическом

разложении решения. 128

п. 4.4.1. Способ построения решения на основе разложения давления..........................132

§ 4.5. Применение метода конечных разностей к восстановлению изображений 132

ГЛАВА 5. КРИТЕРИЙ УСТОЙЧИВОСТИ ДЛЯ ЯВНЫХ СХЕМ НА ОСНОВЕ КОНТРОЛЬНЫХ ОБЪЁМОВ ДЛЯ УРАВНЕНИЙ МАКСВЕЛЛА....................................138

§ 5.1. Существующие критерии устойчивости и понятие энергетической нормы 138

§ 5.2. Шаг по времени в терминах энергетической нормы 141

§ 5.3. Шаг по времени для метода контрольных объёмов для уравнений Максвелла 143

п. 5.3.1. Схема первого порядка точности для уравнений Максвелла.............................143

п. 5.3.2. Выражения для максимальной энергии................................................................144

п. 5.3.3. Оптимальный способ вычисления критерия устойчивости..............................145

п. 5.3.4. Критерий устойчивости для схем высокого порядка точности.......................147

п. 5.3.5. Численные эксперименты......................................................................................152

§ 5.4. Явные схемы интегрирования по времени 153

п. 5.4.1. Предиктор-корректор...........................................................................................153

п. 5.4.2. Схемы Рунге-Кутта для интегрирования по времени........................................154

§ 5.5. Явные схемы с локальным шагом по времени 157

п. 5.5.1. Согласованные шаги по времени, набранные по степеням двойки....................159

п. 5.5.2. Произвольные согласования в шагах по времени.................................................160

п. 5.5.3. Некоторые аспекты реализации локальных шагов по времени.........................162

п. 5.5.4. Оптимизация выбора локальных шагов по времени............................................163

ГЛАВА 6. ПРИМЕРЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ ДЛЯ МЕТОДА КОНТРОЛЬНЫХ ОБЪЁМОВ..........................................................................................................................166

§6.1. Течение в каверне 166

п. 6.1.1. Описание численного эксперимента.....................................................................166

п. 6.1.2. Описание полученных результатов......................................................................167

§ 6.2. Течение в каверне с вырезанной фигурой произвольной формы 170

п. 6.2.1. Описание численного эксперимента.....................................................................170

п. 6.2.2. Описание полученных результатов......................................................................170

§ 6.3. Течение в канале с уступом 171

п. 6.3.1. Описание численного эксперимента.....................................................................171

п. 6.3.2. Описание полученных результатов......................................................................172

§ 6.4. Обтекание полуцилиндра в канале 172

п. 6.4.1. Описание численного эксперимента.....................................................................172

п. 6.4.2. Описание полученных результатов......................................................................173

§ 6.5. Течение в области произвольной формы 174

§ 6.6. Расчет движения проводящей жидкости методом конечных разностей под воздействием

магнитного поля 175

§ 6.7. Расчет уравнений Максвелла, модель точечного электрического диполя 178

§ 6.8. Расчет отраженного поля от объектов 182

п. 6.8.1. Отражение от металлической сферы.................................................................182

п. 6.8.2. Отражение от диэлектрического проводящего куба........................................186

§ 6.9. Численные эксперименты с различными граничными условиями для уравнений Максвелла 186

п. 6.9.1. Отражение от металлической сферы, заключенной в кубическую расчетную

область..................................................................................................................................186

п. 6.9.2. Отражение от двух кубических металлических объектов................................188

§ 6.10. Расчеты движения электромагнитного поля на основе моделей тонких проводов 190

п. 6.10.1. Моделирование антенны штырь и антенны диполь.........................................190

п. 6.10.2. Расчет резонансных частот для металлического резонатора.......................192

ЗАКЛЮЧЕНИЕ...........................................................................................................................194

ЛИТЕРАТУРА...............................................................................................................................196

Введение

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

Набор численных методов для уравнений Максвелла имеет свою специфику и многие методы, применимые для расчета движения сплошных сред не применяются для расчета уравнений Максвелла. Более того, многие численные методы для моделирования распространения радиоволн основаны на интегральных уравнениях и позволяют рассчитывать задачи только лишь для отдельной частоты, что не является приемлемым для расчета движения токопроводящих сред в магнитном поле. Именно поэтому возникает необходимость обобщения численных методов, предназначенных для расчета движения жидкостей и газов на задачи о распространении радиоволн. Более того, многие численные методы для моделирования движения радиоволн основаны на полуаналитических соотношениях, не обобщаемых на задачи о движении жидкостей и газов. Одним из важных приложений численных методов к расчету уравнений Максвелла является получение форм электромагнитного поля, оптимально распространяющегося через волноводы определенной формы. Если в 60х годах двадцатого века для расчета волноводов требовалось использовать полуаналитические численные решения [11], специфичные для каждой из задач, то в настоящее время в связи с прогрессом в численных алгоритмах расчета уравнений в частных производных можно получать численные решения для таких задач на основе более универсальных численных методов. Последние достижения в методах расчетов на нерегулярной сетке позволяют моделировать, в том числе и сложные приборы, такие как антенны сотовых телефонов и оборудование сопряжения с ними [12-17]. К примерам таких задач относятся [11,18,19] в которых приведены полуаналитиче-

5

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

Примерами таких задач служат резонирующие антенные комплексы, такие как [20,21], для расчета которых часто используются алгоритмы основанные на вычислении характеристик распределения радиоволн одной характерной частоты [22]. Последние алгоритмы зачастую используют интегральное уравнение, решаемое на поверхности устройства без использования объёмных сеток в которых для ускорения произведения матрицы аппроксимации интегрального уравнения на вектор используется идея быстрого преобразования Фурье с использованием предобуславливателя на псевдо Теплицевой матрице [173]. Анализ поведения приборов на основе алгоритмов, рассчитывающих распространение радиоволн для отдельных частот, требует проведения целой серии расчетов для каждой из гармоник. Однако, моделирование распространения радиоволн в нелинейных средах должно учитывать переизлучение энергии на некотором спектре частот, отличном от частот источника излучения, что делает расчёт и анализ таких проблем затруднительным. Использование временной зависимости для расчета компонент электромагнитных полей позволяет моделировать поведение нелинейных сред в широком диапазоне частот и для сигналов различной формы, что открывает возможность моделирования движения токопроводящих сред, модель движения которых описывается системой уравнений Максвелла для переменных, зависящих от времени, совместно с уравнениями Навье-Стокса. Последнее позволяет рассчитывать движение плазмы или электропроводящей жидкости под воздействием электромагнитных полей, что делает возможным моделирование сложных химических реакций, движения плазмы вокруг летательных аппаратов, магнетроны, магнитогидродинамические генераторы [176] и прочие современные устройства, используемые для создания наночастиц или пучков плазмы.

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

мированию систем уравнений в частных производных важно изучить подходы к аппроксимированию уравнений Максвелла методом контрольных объёмов. Необходимо отметить, что применение метода контрольных объёмов к решению уравнений Максвелла не было в достаточной мере изучено и первые широко известные работы [62,65] по этому методу датируются 1999-2000 годами. Последнее связано в первую очередь с высокой диссипативностью метода контрольных объёмов. Тем не менее обобщение метода контрольных объёмов на уравнение Максвелла привело к необходимости изучения аппроксимации векторных законов сохранения и соответственно расширению набора аппроксимационных способов для решения мультифизических задач.

Уравнения Максвелла, записанные во временных переменных, подобны по своей структуре уравнениям Навье-Стокса для расчёта движения нестационарных сред [26]. При исследовании проблем связанных с движением жидкости и газа и воздействия их на обтекаемые ими твёрдые тела, рассматривается модель сплошной текучей вязкой среды, которая подчиняется уравнениям Навье-Стокса. Одна из простейших моделей гидродинамики - модель несжимаемой ньютоновской жидкости, которая при отсутствии воздействия внешних сил имеет вид:

|р(и,+ й'Уи)-цУ2н+ V р=0,

[У-и=0,

где р - давление, р - плотность, V - коэффициент кинематической вязкости, / - время, и -вектор скоростей.

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

Решение уравнений Навье-Стокса методом контрольных объёмов основано на применении скалярных законов сохранения, применяемых при аппроксимации дивергенции. В случае применения уравнений Максвелла дивергентные законы сохранения необходимо обобщить на векторные законы сохранения, используемые для аппроксимации уравнений Максвелла на основе оператора ротора, которые можно записать в виде:

UдlE-VxH+oE = -J, [ \хд,Н+^хЕ = 0,

где Е - электрическое поле, Н - магнитное поле, а - проводимость среды, е - эл�