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

'"V . . . . •

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

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

Шурина Элла ГТетровна

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

ертация на соискание ученой степени доктора технических наук в форме научного доклада

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

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

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

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

профессор Ю.Е. Воскобойников

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

доктор технических наук,. профессор О.П. Солоненко

Ведущая организация: Институт вычислительных технологий СО РА1

Защита состоится _ %У97 года в .часов на засе

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

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

Диссертация в форме научного доклада разослана '¡/¡0" ¿И^СЦ

Ученый секретарь диссертационного совета . к.т.н, доц. Г.П.Чикильди

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

1.1. Актуальность темы

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

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

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

1.2. Цель работы

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

1.3. Научная новщна

I [аучная иовиша определяется следующими результатами :

• разработаны вычислительные схемы н алгоритмы моделирования теплов состояния мощных турбогенераторов традиционного (газовое, водяное ох ждение) и нетрадиционного (использование явления сверхпроводимости) полнений;

• предложены математическая модель и алгоритм конечно-элементного мо лирования трехмерных тепловых полей с учетом радиационного теплооб на;

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

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

1.4. На защиту выносятся:

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

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

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

• Конечноэлементные модели радиационного теплообмена.

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

• Результаты численных расчетов трехмерных магнитных и тепловых поле реальных технических устройствах.

1.5. Практическая ценность работы

В работе получены следующие практически значимые результаты :

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

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

танных методик конечиоэлементного моделирования для решения сложных практических задач;

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

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

Результаты работы были использованы в различных организациях, в том где, во ВНИИэлектромаш (г. Санкт-Петербург), ИВТАН (г. Москва), ЛПЭО лектросила" (г: Санкт-Петербург), ИЯФ СО РАН (г. Новосибирск). Имеются гы о внедрении результатов моделирования.

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

Личный вклад. Постановка изложенных в работе задач, выбор метода их шпации. принадлежат лично автору. Разработка алгоритмов моделирования нового состояния объектов, состоящих из разномасштабных фрагментов с ¡личными теплофизическнми свойствами, с различными типами тепловых 1тактов; а также алгоритмов решения задач магнитостатики выполнены под човодством автора и при его непосредственном участии.

1.6. Апробация работы

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

Всесоюзных конференциях по автоматизации научных исследований на юве применения ЭВМ (Новосибирск, 1977,1981); Всесоюзных семинарах по индексам программ математической физики (Горький, 1981, Новосибирск, 32,1984, Ташкент 1983, Красноярск, 1988 ); Всесоюзной школе-семинаре по "оматизации научных исследований (Новосибирск, 1985); Второй Всесоюзной тференцин по техническому использованию сверхпроводимости (Ленинград, 33); Всесоюзном научно-техническом семинаре "Опыт применения средств дшческой диагностики и контроля за состоянием электроэнергетического эрудоваиия" (Суздаль, 1986); Всесоюзной конференции "Актуальные про-:мы вычислительной и прикладной математики" (Новосибирск, 1987); VI ¡союзном семинаре "Обратные задачи и идентификация процессов теплооб-на'"(Москва, 1987); Всесоюзном научно-техническом совещании "Вопросы

проектирования, исследования н производства мощных турбо-, гидрогенераторов и крупных электрических машин" (Ленинград, 1988); Всесоюзном научно-техническом совещании "Автоматизация проектирования и производства в электромашиностроении" (Суздаль, 1989); Всесоюзной конференции "Теплообмен в парогенераторах" (Новосибирск, 1988); Всесоюзной конференции "Математическое моделирование: нелинейные проблемы и вычислительная математика" (Звенигород, 1988, 1990); Всесоюзной конференции по методам численного решения многомерных нестационарных задач математической физики (Арзамас-16, 1991); Всесоюзной научно-технической сессии "Состояние и перспективы газодинамических тепловых исследований в обеспечении повышения температуры газа в стационарных газотурбинных установках" (Москва, 1989); Международном совешании-семинаре по механике реагирующих сред и экологии (Томск, 1994); Международной школе-семинаре "Аналитические методы и оптимизация процессов жидкости и газа" САМГОП-94 (Арзамас-16, 1994), First Asian Computational Fluid Dmamics СопГегепсе (Hong Kong, 1995); Particle Accelerator Conference and International Conference on High-Energy Accelerators (Dallas,1995); Международной конференции PaCT-95 (Санкт-Петербург, 1995); Международной конференции "Математические модели и численные методы механики сплошных сред" (Новосибирск, 1996); семинарах ВЦ СО РАН, ИВТ СО РАИ, ИТПМ СО РАН (Новосибирск), ВЦ РАН (Москва), Научного Совета АН СССР но комплексной проблеме "Научные основы электрофизики и Электроэнергетики", секция № I 'Теоретические проблемы генерирования электромагнитной энергии" (Ленинград). Результаты работы были использованы а 8 отчетах по научно-исс^довательской работе.

По тематике диссертации под руководством автора защищены 4 кандидатских диссертации.

1.7. Публикации

Основные материалы научного доклада опубликованы в 68 печашых работах.

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

2.1. Моделирование теплового состоянии энергетических объектов

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

промышленной реализации конструкторско-технологических решений. Рассмотрим тепловое состояние турбогенератора (мощность свыше 500 мВт), характерное для его неподвижной (статор) и подвижной (ротор) частей [3, 6-10, 13, 14, 16,18-20].

Для машин традиционного исполнения вероятны технологические дефекты, приводящие к появлению очагов повышенного тепловыделения в сердечнике статора, которые необходимо идентифицировать на этапе тепловых испытаний. Процедура идентификации источников тепловыделения включает в себя два класса задач: во-первых, создание набора стандартных образов (прямые задачи); во-вторых, обработка результатов измерений и алгоритм для сравнения этих результатов с результатами решения прямых задач [3]. Прямые задачи описываются трехмерными уравнениями теплопроводности, с заданным начальным распределением температуры, краевыми условиями Дирихле и Неймана на различных границах области решения, множеством разнотипных источников тепловыделения [1, 3, 7, 11, 12]. После решения серии прямых задач была определена стратегия проведения стендовых тепловизионных измерений, а именно: измерялись значения трехмерного температурного поля в плоскости (г,<р), но для нескольких моментов времени, что позволило идентифицировать как глубинные, так и пршраничные источники тепловыделения [I, 2,3,4, 7,12, 16, 17). Для прогнозирования теплового состояния стагора я рабочем режиме решалась стационарная задача, в которой функция источника определялась из решения задачи идентификации, что позволяло определить температурное поле в статоре в установившемся режиме и прогнозировать его работоспособность [3].

Термомеханическая надежность вращающегося фрагмента ту рбогенератора (ротора) традиционного исполнения определяется системой охлаждения обмоток. Описание этих задач и некоторые подходы к их решению рассматривались в работах И.Ф.Филиппова, Г.М.Хуторецкого. В [16-19] приведены математические модели процессов массо- и теплопереноса, алгоритмы и результаты численного моделирования для различных вариантов технологического исполнения каналов системы охлаждения ротора турбогенератора. Результаты моделирования позволили сделать вывод об оптимальном соотношении (скорость хладагента, градиент давления, температура хладагента на входе-выходе канала системы охлаждения) и конфигурации каждого элемента системы охлаждения ротора турбогенератора (ТВВ-500-2Е, ТВФ-110).

Цель исследования температурного поля ротора турбогенератора нетрадиционного исполнения со сверхпроводящей обмоткой (СПТГ) заключается в определении оптимального конструкторского решения, минимизирующего тепло-притоки в низкотемпературную зону, расход гелия, обеспечивающий состояние сверхпроводимости обмотки при различных режимах работы объекта. Математическая модель температурного поля ротора СПТГ, соответствующая ей вариационная постановка и алгоритмы моделирования изложены в [4-6, 8-11]. Ис-

следование расходных характеристик хладагента (гелия) в роторе СПТГ пот] бовало решения задач о движении жидкости (газа) в криволинейных враша щнхся каналах. Математические модели, характеризующие эти процессы, вариационные постановки, конечноэлементные вычислительные схемы, резу. таты вычислительных экспериментов приведены в [5, И, 15, 16, 24].

В результате вычислительных экспериментов, реализующих многовариантн расчеты при создании и исследовании мощных электрических машин традишк ного и нетрадиционного исполнении, было предложено технологическое peniei для системы каналов газового охлаждения роторов турбогенераторов ТВВ-5 ТВВ-500-2Е, ТВФ-110, позволяющее осуществить оптимальное соотношеши скорость хладагента - тепловое состояние обмотки ротора. Для СПТГ исследов; влияние размеров поперечного сечения каналов теплообменника на газолинами скне характеристики течения хладагента, влияния различных материалов, при; няемых в конструкции КТГ-20, величины контактных сопротивлений, степени куумирования полостей ротора КТГ-20, что дало возможность выработать ре мендации по выбору проекта конструкции ротора СПТГ. - -

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

2.2.1. Методика расчета трехмерною магнитного поля с использованием двух потенциалов

Одним из наиболее эффективных ^етодов численного решения трехм пых магнитостатических задач является метод, основанный на использован двух потенциалов; полного потенциала у в области Q'1' с ферромагнитными i териалами и неполного потенциала л в области Qp, содержащей токовые обм ки. Впервые, в 1979 году, такой подход был предложен J.Simkin, C.W. Тгс bridge. В нашей стране он использовался в некоторых работах, проводимых i руководством В.П.Ильина. Суть этого подхода состоит в следующем.

В области полного потенциала Q'1' напряженность магнитного поля Н ni ностью определяется градиентом полного потенциала

Я = ЯЧ'=- grady. (

В области неполного потенциала Л определяется суммой двух полей:

H = Hc + Hp = Hc-&adpt (

где Нс - напряженность магнитного поля, создаваемого токовыми обмоткам однородном пространстве (т.е. в пространстве с относительной магнитной п ницаемостью Ц=1), а Нр ~ - grad р -напряженность дополнительного магн

то поля, создаваемого ферромагнитными материалами. Очевидно, что вектор-/нкция И, определяемая уравнениями (1) и (2), удовлетворяет уравнению \Н - 0 в области £У полного потенциала и уравнению т\Н = 3 в области ¡полного потенциала (] - плотность тока в возбуждающих магнитное поле ¡мотках).

Потенциалы у и р могут быть найдены как решение дифференциального ювнения

с1ИрЙ) = 0 (3)

области П=П,|Ш\ Для получения эквивалентной вариационной формулиров-I уравнение (3) домножается на пробную функцию о и интегрируется по об-1сти Г2:

|ШУ(|ЛЯ)ОС/П = 0. (4)

«

Учитывая, что область О является объединением областей £У и П4', пере-жаюшихся по их общей границе Л", и в области р=1 и ¡Iе = 0, из (4) поучаем

^(рЯ'ЧиЛГ + =0. (5)

и41 И*

Применяя к каждому из слагаемых левой части уравнения (5) формулу рина (интегрирования по частям) и учитывая, что на внешней границе области ! либо Н=0, либо и=0 (из принципа большого объема или условий симметрии агнитного поля), получим

- ёгас) ЫГГ + - р/'^и сЮ." + |яр= 0, (6)

о''' Л* а* у*

ае Л"1' и ,4" - одна и та же поверхность 5" (граница между £У и с внеш-ими нормалями п" и пг соответственно, причем п" =-йр, т.е. с1§'г =Яч'(/5", а (V" = -я*.Таким образом, (6) преобразуется к виду

- |ц//'" ёгаёЫП4' - вгас1и<ЯУ + - Нр)пос18" =0, (7)

делзнт.

Учитывая, что на поверхности выполняется равенство р.Н4' = Н" + Не, ■случаем,.что подиитегральное выражение р//'1' - Ир в поверхностном инте-

трале равно Н° ■ Тогда уравнение (7) преобразуется к виду

- ргас! шЛУ - |цгас1 р цгас1 ш/П" + л и<Д''' = 0. (8)

I.» и'' X"

Вариационное уравнение (8) обеспечивает непрерывность нормальной составляющей Й = ц„цЯ на границе У = О '1 ПП?. В конечноэлементной формулировке в качестве и берутся базисные функции <р„ и поверхностный интеграл

^Н'писМ" определяет одну из составляющих вектора правой части конечной-

элементной системы линейных алгебраических уравнений (СЛАУ).

Для того, чтобы однозначно определить решение уравнения (8), на границе области необходимо задать соотношение между потенциалами V)) н р, например, в виде

у=/>+м. V (9)

Функция м-ц/-/? называется скачком (или разрывом) потенциалов. Она строится таким образом, чтобы обеспечить непрерывность тангенциальной составляющей Л, вектора напряженности /7 на границе Л":

Й? = Н? + Н[.

Учитывая (1) и (2), получим

- (grad 41), = -(grad />), + /?;.

Подставим вместо ч* его значение из (9): -(grad (/> + «)), = -(grad/7)t +/?;,

или

- (grad u), =/?;".

(Ю)

Уравнение (10) определяет скачок потенциалов и на границе Л"' таким образом, что тангенциальные составляющие векторов Нчц //'" + 1Г совпадают на

границе 5й между областями £Х" и Пр, если лотенциалы ц/ и р связаны соотношением (9). Тем самым достигается непрерывность тангенциальной составляющей вектора И на S1, необходимая для выполнения уравнения rot Я = J в области Q, непрерывность нормальной составляющей II на Л*', необходимую для выполнения уравнения div Д = 0 в области П, обеспечивает вариационное уравнение (8).

2.2.2. Вычисление скачка потенциалов

Правая часть уравнения (10), определяющего скачо* потенциалов на поверхности Л" = П'' ПО'', содержит касательную составляющую напряженности магнитного поля токовых обмоток П[. Величина Н' почти во всех практических задачах не имеет аналитического представления, но может быть достаточно легко рассчитана в любой точке поверхности X". Поэтому процедуры вычисления скачка и должны учитывать тот факт, что поле /Г известно только в отдельных точках поверхности

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

Ям

= (11).

д1,

где /, - ребра треугольников .V", а Н' ~ проекция Н' на ребро /,.

Одним из возможных методов вычисления скачка потенциалов является метод, основанный на интегрировании уравнений (11) по ребрам /,. Для этого в одном или нескольких узлах, являющихся вершинами треугольников ЯЦ, задаются значения и=0 (как правило, это узлы, в которых потенциалы \у=р=0), и с помощью квадратурных формул по ребрам треугольников , исходящим из этих узлов, интегрируется уравнение (11). В результате в узлах, расположенных на противоположных концах ребер, вычисляются значения скачка^ потенциалов и. Затем определяется следующее множество ребер, исходящих из узлов с уже вычисленными значениями и, и по ним проводится интегрирование уравнения (11). Этот процесс повторяется до тех пор, пока значения и не будут вычислены во всех вершинах треугольников Л'^' на поверхности = О'' ПОр .

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

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

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

Более высокую точность вычисления и могут дать процедуры интегрирования уравнений (И) с осреднением значений к, вычисленных по различным путям. Однако при использовании квадратурных формул для интегрирования (11) введение процедур осреднения делает алгоритм вычисления и очень сложным, поскольку различных путей, ведущих от узлов со значениями м=0 к любым другим вершинам треугольников .S'j', может быть слишком много, и поэтому необходимо разрабатывать способы отбора путей интегрирования для усреднения, а это является совсем нетривиальной задачей.

Устранить сложности, связанные с осреднением результатов интегрирования уравнения (И) по различным путям, позволяет подход, основанный на определении разрыва потенциала и из решения вариационной задачи [38]

/(«) = > II f" + ni\ <", ""in • (12)

"Ii!"

1 I.

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

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

5/ ,

ников Я]?.

Для решения вариационной задачи (12) будем использовать метод конечных элементов (МКЭ). Будем считать, что функция и является линейной комбинацией базисных функций ср,, заданных на ребрах /,:

» = (13)

I

Подставляя (13) в (12), получим

'-siHf'

• I, ■ -

'<«>-2; I WTrjr ^^f^iHHiy К- (14)

<

Дифференцируя (14) по весам </> и приравнивая производные к нулю, получаем СЛАУ для вектора весов д :

.</,'' 1 А

\

Для вычисления вектора правой части /*' = - //¿с//, необходимо за-

I,

дать вид Н,сна ребрах /¡. Даже на не слишком подробных конечноэлемеитных сетках высокую точность вычисления обеспечивает аппроксимация квадратичными полиномами /// на /,. Такая аппроксимация может быть легко построена

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

Для достижения необходимой точности решения задачи (12) при аппроксимации /7е на /, квадратичными полиномами ребра I, могут быть разбиты на несколько конечных элементов. Наиболее удобно контролировать точность вычисления и при разбиении /, на кратное 2 число элементов. Матрица СЛАУ

хранить в разреженном строчном формате. При удвоении числа конечных элементов на всех ребрах /, портрет матрицы модифицируется очень легко, поэтому процесс генерации является нетрудоемким. Гораздо более трудоемким является процесс решения СЛАУ (15) методом сопряженных градиентов с предобу-слойливаннем. Точность вычисления удобно контролировать с помощью относительного отклонения в вершинах треугольников Я? при разбиении ребер /, на 2" и 2л+| конечных элементов.

2.2.3. Реализация метода Ньютона при решении трехмерных

Для вычисления потенциалов *|» (полного) и р (неполного) будем использовать полученное выше вариационное уравнение (8). Будем считать, что расчетная область 0=0'1ипр представлена в виде совокупности согласованных конечных элементов Г2„, с узлами (ц. Каждому узлу О-, соответствует непрерывная кусочно-полиномиальная базисная функция <р,, равная единице в С-, и нулю во всех остальных узлах (<р;зО на всех конечных элементах, не содержащих узел О,). Будем также считать, что поверхность скачка потенциалов Л"=Пч'ПОр представлена в виде совокупности поверхностных элементов являющихся гранями конечных элементов П,,..

имеет очень разреженную структуру, и поэтому ее удобно

нелинейных магннтостатнческнх задач с использованием двух потенциалов

Будем искать полный потенциал у и неполный потенциал р в виде линейных комбинаций базисных функций <р,:

= ХУ41" (16>

е/Ш")

>с/<и')

где /(О4") - множество индексов узлов конечноэлементной сетки, принадлежащих области О4" вместе со всеми ее границами, включая Л"=ПтПО'', а /(О1') -множество индексов узлов сетки, принадлежащих области также вместе со всеми ее границами. Подставляя (16) и (17) в вариационное уравнение (8), получим

- |ц( I/* цШр^цтадиЛП* -

(18)

- Г У^,'' цга(1(р, + [//'и ис/Г = 0.

- ипг/(Ч') у"

Обозначим через /(Л") множество индексов узлов сетки, принадлежащих поверхности •У=0>ТЮр. Тогда из равенства (9) следует

< =</'+",. ¡е1(Х"), (19)

где и, - коэффициенты разложения скачка потенциалов и по базисным функциям ср, на поверхности Л" Обозначим через ¡1 вектор, компоненты которого </, определяются соотношениями

\с]/е/ССТ) . Ч, = Г (2°)

а через ¡1" - вектор, компоненты которого су," определяются соотношениями

?.={ (21) (О, 1« /(Г).

Таким образом, компоненты у, вектора ц определяются значениями ц'" во всех узлах области Г241 вместе с границей У и значениями в тех узлах области Пр, которые не лежат на границе У; вектор ¿¡" имеет ненулевые компоненты, равные скачку и, только в тех узлах сетки, которые принадлежат поверхности б*1. Размерности векторов </ и ¿¡" совпадают и равны полному числу узлов в ко-

нечноэлементной сетке V. Вектор </" является известным из соотношения (21) после вычисления скачка потенциалов г/, а вектор ¿} подлежит определению. Сучетом соотношений (19)и (21)уравиенне (18) преобразуется к виду

- Гш ^¡Г </, рте! <р, I grad и ¿/Г2 +

I! ^ ,„1 1

+ (/; ёгаё grad и<К1+ |/?сп и ¿Я" = О,

(22)

где V - полное число узлов в конечноэлементной сетке. Заменяя в вариационном равенстве (22) пробную функцию и поочередно на все базисные функции Ф= I,..., V, получим для компонент </, вектора д систему из V уравнений

рад ф, grad </, = ^Г ф( grad

</Г +

/ = !,....V

(23)

или, в матричном виде Ац = I- .

Компоненты А

(24)

(25)

матрицы А являются функциями магнитной проницаемости р, которая, в свою очередь, является функцией напряженности магнитного поля Н. А поскольку в области П'г И = - цгж! ц/ = - бгпс1 <р1, то компоненты А/, фактически явля-

117(11»)

ются функциями вектора </ : Ак = Таким образом, конечноэлементная

аппроксимация вариационного уравнения (8) приводит к системе уравнений (24) с матрицей А, компоненты которой (25) зависят от искомого решения ¿¡ . Поэтому система уравнений (24) является нелинейной, и ее удобно записать в виде

у

Рассмотрим решение конечноэлементной системы нелинейных уравненш (26) с использованием метода Ньютона.

На каждом шаге итерационного процесса решения системы нелинейны) уравнений (26) методом Ньютона необходимо выполнять линеаризацию все) уравнений решаемой системы. Обозначим через ¿/" решение, полученное н; предыдущей итерации процесса Ньютона (для первой итерации Г/" - начально! приближение). Выполним линеаризацию ;-го слагаемого левой части /-го урав нения системы (26) с использованием его разложения в ряд Тейлора в окрест ности точки у":

Входящие в (27) частные производные й(/1(|(7)<7,)/37 имеют следующш

вид

Таким образом, для линеаризации системы уравнений (26) необходимо вы числить частные производные компонент (25) матрицы А по переменным ц^ Однако для практической реализации такой подход к линеаризации уравненш конечноэлементной системы (26) очень неудобен, так как сборку глобально! матрицы конечноэлементной системы уравнений легко вести для компонент являющихся числами, и очень трудно для компонент, которые являются неко торы ми функциями вектора решения <7. Поэтому процедуру линеаризации от дельных слагаемых системы (26) лучше перенести на этап генерации локальны: матриц отдельных конечных элементов, и собирать из локальных матриц и век торов уже линеаризованные уравнения конечноэлементной системы.

Будем считать, что расчетная область П представлена совокупностью тет раэдральных конечных элементов Пш, на которых заданы кусочно-линейны базисные функции ч>,. Тогда, учитывая соотношения (1) и (16), получим, 'п внутри каждого конечного элемента ПтеО'|/ напряженность магнитною поля !-, определяется соотношением

(27)

¿4, сЦ,(<7)

ч,. -I = ' >

(28)

(20)

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

гает номера четырех вершин тетраэдра Пт). В этом случае квадрат модуля Н на конечном элементе Йт также является константой и определяется соотношением

(30)

Проинтегрируем обе части соотношения (30) по конечному элементу П,„ и разделим их на объем |ПЯ! этого конечного элемента:

| цгаёфц ёгас¡<р//П

Мг-

(31)

Обозначим через 7т матрицу, компоненты которой определяются соотношением

1 = ^цгас! ф,. цга<1<[>гг/0.

(32)

Тогда на конечном элементе От соотношение (31) с учетом (32) может быть записано в виде следующей квадратичной формы переменных (/,:

Поскольку напряженность магнитного поля Н на конечном ^элементе Пт является константой, то и относительная магнитная, проницаемость р на От также является константой. Поэтому из матрицы 7га легко получить матрицу вклада от конечного элемента в матрицу А системы (26). Действительно, если обозначить через А™ локальную матрицу элемента П„„ расширенную до размерности и глобальной матрицы А конечноэлементной системы (24), то она может быть представлена в виде

Лт=ц7т. (34)

Матрица А конечноэлементной системы уравнений (26) фактически является суммой таких локальных матриц Ат:

(35)

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

понентами локальных матриц конечных элементов П,„:

А(<7Н=]Г Л:(</)</,. (36)

т

Поэтому вместо линеаризации слагаемых Аи ((/)</, системы конечноэле-ментных уравнений (26) можно выполнять линеаризацию всех слагаемых А"('])'!, на этапе генерации локальных матриц, и только после этого собирать глобальную матрицу конечноэлементной системы уравнений. Очевидно, что в результате таких действий мы получим ту же самую линеаризованную систему, а выполнять линеаризацию каждого слагаемого нелинейной конечноэлементной системы уравнений гораздо удобней на этапе генерации локальных матриц.

Линеаризация слагаемых Л,"'(</)'/, выполняется так же, как это было сделано в формуле (27) для слагаемого Л,,(<7)</,:

ЛМ)«,*^),:^^^, ' (37)

г ъ,

где Г}" - решение на предыдущей итерации по нелинейностям, а частные производные д( А"(с})у, )/<?</, имеют вид

д(А;;(<1)<1,)

(38)

Вычислим производные 8А"{Г1)! ёц1. Для этого воспользуемся соотношением (34) и учтем, что компоненты '![" матрицы V" не зависят от (/,:

Магнитная проницаемость ц зависит от опосредованно, ц является функцией только модуля напряженности Н, которая в свою очередь является некоторой функцией переменных Нам удобно будет воспользоваться квадратичной зависимостью квадрата модуля // от переменных (¡¡, и поэтому мы будем рассматривать р как функцию Н1. Тогда соотношение (39) перепишется в виде:

алд?) _ Ф ОН1

Величина (1^с!Нг вычисляется по кривой ц(Н2) или р(//). характеризующей магнитные свойства материалов в области О41, а значения 3//а /5(/; легко получить из соотношения (33):

Подставляя (41) в (40), получаем

дА',"(с/) _ 2 ф

•С ' 1

а, п„.,,//'■•!>•• . <42)

Преобразуем правую часть соотношения (37) с учетом соотношений (38) и

(42):

, и<°

= а;:(Г) ч, + 2$ ¿4 -ч))-. (43)

/ЯП г

где Л'ЩГу) - компонента локальной матрицы, определяемой соотношением (33), в котором магнитная проницаемость ]Ц вычисляется по значению Н2 из формулы (33) при Ц = ч".

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

Обозначим через А" локальную матрицу, а через / - локальный вектор конечного элемента Я,„, полученные после линеаризации вкладов от £)„ в нелинейную систему (26). С помощью соотношения (43) соберем вклады в (т.е. в сомножитель перед от всех нелинейных членов 1-го уравнения

системы, т.е. для всех 7 е/(£1т) (где /(£}„,) - множество всех номеров вершин конечного элемента Г2,„):

Из (43) легко получить и вклады в . Для этого нужно просуммировать по / все слагаемые правой части соотношения (43), не содержащие ¡¡, в виде сомножителей, и затем сменить у них знак, т.к !■]"' находится в другой части линеаризованного уравнения:

2 (/ц

'п-мг

(45)

Соотношения (44) и (45) можно записать в более компактном виде, если обозначить через вектор й'"' произведение матрицы /"" на вектор (/":

= * (46)

С учетом (46) соотношения (44) и (45) можно переписать в следующем, более удобном для вычислений виде:

' п,„: с/и2

(47)

(48)

Таким образом, матрица А и вектор !■' линеаризованной системы

Ад = (49)

получаемой в результате замены в системе нелинейных уравнений (26) отдельных слагаемых ее уравнений на их разложения в ряд Тейлора (27), собираются стандартным образом из локальных матриц А'" и векторов /•'"' конечных элементов От, определяемых формулами (47) и (48). Эго относится к процессу генерации глобальной матрицы и вектора по конечным элементам Пш, расположенным в области £1'*' полного потенциала, содержащей в себе все магнитные материалы. При учете вкладов в А и Г от конечных элементов, находящихся в области неполного потенциала Ор, локальная матрица Л'" полностью совпадает с матрицей 7™ (см. формулу (32)), поскольку в области неполного потенциала ц=1, а локальный вектор К", определяемый объемным интегралом из правой части уравнения (23), вычисляется по формуле

Рт~Ттч4. (50)

Вклады в глобальный вектор правой части /•' линеаризованной системы уравнений (49) от поверхностного интеграла из правой части уравнения (23) удобней учитывать отдельным проходом по всем треугольникам поверхности до или после выполнения рассмотренной выше полной процедуры генерации глобальных матрицы и вектора линеаризованной системы уравнений.

Отметим, что применение рассмотренной процедуры линеаризации может ухудшить свойства СЛАУ по сравнению с методом простой итерации, в котором глобальная матрица конечноэлементной СЛАУ является матрицей /1(г/°) . Действительно, при <]^<1Нг<0, что справедливо для большей части кривой намагничивания, добавки от линеаризуемых членов могут нарушить неотрицательную определенность матрицы, характерную для (см. соотношение (47)). Это может сказаться на свойствах глобальной матрицы линеаризованной СЛАУ и, как следствие, на скорости сходимости итерационного процесса ее решения. Еще более важным является тот факт, что при значительных изменениях на соседних итерациях величины ф/<#/2 итерационный процесс учета нелинейности, обусловленной зависимостью р от И, может войти в колебательный режим, существенно замедляющий его сходимость. Поэтому имеет смысл ввести параметр торможения этого процесса хе[0,1 ] и модифицировать формулы (47) и (48) следующим образом:

к = (51)

р,„\с!Нг 'V 'I

Очевидно, при х=0 мы получим процедуру решения системы нелинейных уравнений (26) методом простой итерации, а при х=1 - процедуру решения системы (26) с линеаризацией в соответствии с методом Ньютона. Конечно, значение % желательно выбрать таким, чтобы суммарное количество итераций метода сопряженных градиентов (МСГ) с учетом количества итераций по нелинейно-стям и количества итераций МСГ при решении линеаризованной системы (49) на одной итерации по нелинейности было минимальным. Ниже в п.2.4.2 будет исследована зависимость вычислительных затрат от коэффициента релаксации 9 и коэффициента торможения % для некоторых практических задач магнитостатики.

Рассмотрим проблему задания кривых намагничивания и вычисления производных ц по Я1.

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

(52)

"у

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

Для реализации описанного выше подхода к решению магнитостатическо! задачи необходимо вычислять производную ц по Я2 (JyüäH1). Кроме того, npi вычислении компонент локальной матрицы АЦ'(<]Ь), входящих в соотношени (51), требуется значение р. Оно может быть получено по значению Я2, доста точно просто вычисляемому по формуле

II* = У [</-'}'7"<7". (53)

которая сразу следует из формулы (33) после замены в ней вектора с/ на векто г/°. Поэтому может показаться, что наиболее удобным для реализации описат ной выше вычислительной схемы является кубический сплайн, интерподирук тин зависимость ¡¡(II2). Однако, это ие совсем так. Дело в том, что кубически сплайн, построенный по набору значений р и может иметь участки иемонс томности, особенно в конце кривой намагничивания, не характерные для peani ной зависимости магнитной проницаемости от величины магнитного поля, в з время как кубический сплайн, построенный по аналогичному набору значены р и //, таких участков немонотонности не имеет, и в этом слу чае зависимое! р(Н~) также не имеет участков немонотонности. Появление участков немон< томности в зависимости р(Я2) приводит ¡^существенному замедлению сходили сти итерационного процесса Ньютона и возникновению в нем плохо затухав щих осцилляции. Поэтому кубический сплайн лучше строить по набору знач* ний р и//, т.к. в этом случае дополнительные вычислительные затраты на этаг генерации локальных матриц являются незначительными, а ускорение сходим* сти итерационного процесса Ньютона становится существенным.

Будем считать, что зависимость магнитной проницаемости от величин магнитного поля задзна кубическим сплайном р(Я), построенным по набо| значений ри, Яа, а=1,..,(}, и для любого значения Н можно вычислить по этод сплайну значение р и его производной с/рУ</Я. Тогда для любого Я<//ц (//¡5 - и следнее значение в исходной таблице (Ри,Я„), определяющей кривую намаги чнвания) не представляет никакого труда получить значение с/ц/</Я\ необход мое для вычисления локальных матриц и векторов (см. формулы (51 )-(52)):

_Ф = 1 (V1

с!Нг 2HdH'

Но очень часто таблица (ц,„//ц) обрывается на том месте, где »сличи üBliwlH (jTfj - магнитная проницаемость вакуума) становится достаточно блт

кой к единице, а при решении практической задачи величина Я может превысить последнее табличное значение Яц. Получим выражения, позволяющие вычислять р и с/рШ для Я>Я|) при условии, что для Н>Н\\. Итак, с одной стороны,

~ = рЯ. (55)

Но

С другой стороны, для 11>В\\ (где йц= роРрЯр)

В = Н„ М! __ »„ ^ ЛВ (56)

Ро Мо I1» Но '

Учитывая, что /уро^РрЯр и Д#/|Л<г '//-//ц, соотношение (56) преобразуется к виду

= Рр'^1 + Я - Яц. (57)

И О

Из (55) и (57) получаем

= -//„+//, (58)

откуда

»1 = 1 + 01,-1)—(59) А продифференцировав (59) по Я, получим и значение е/р с/Я: с/Я р Я

= -(рр - о ,/;. (60)

Таким образом, соотношения (59) и (60) позволяют вычислять значения р и с/р/г/Я при значениях Я, превышающих последнее табличное значение Яр, если при Я>ЯР выполняется соотношение с//?/р0с/Я=1.

2.3, МКЭ-моделированне трехмерных температурных полей с учетом радиационного теплообмена

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

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

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

(61)

заданным в подобластях О, расчетной области П (О = 1_)П(), где р - плотность,

Ср- теплоемкость, к - теплопроводность материалов,/- плотность объемных источников тепловыделения, а также краевыми условиями

(62)

' кВМ =6, (63)

5л |4-1

= 0, (64)

№

заданными на внешних границах 5', З2, расчетной области О, и условиями сопряжения

еи

к 8п

-аоЕ((/^-(У^) = 0, (65)

заданными на поверхностях 5„, являющихся границами подобластей Г), (а0 -константа излучения абсолютно черного тела, е - коэффициент черноты, 0 й е £ I). Условия сопряжения (65) определяют теплообмен излучением между границами и подобластей П, и П, соответственно. В начальный момент времени /о температура (/ в области О определяется функцией На.

В краевой задаче (61)-(65) коэффициенты р, С.,, X могут зависеть от искомой функции и. Однако зависимости эти довольно гладкие, и поэтому при решении краевой задачи (6|)-(65) в процессе движения по временным слоям, как правило, бывает вполне достаточно использовать для вычисления этих коэффициентов значения и с предыдущих временных слоев (возможно с экстраполяцией и некоторой релаксацией). И даже при решении стационарной задачи

' ди '

(— = 0) учет зависимости коэффициента Я от температуры (/ с использовани-81

ем метода простой итерации дает вполне приемлемую сходимость по нелиней-

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

Существенное отличие в технику конечноэлементной аппроксимации краевой задачи (61)-(65) по сравнению с аппроксимацией других краевых задач для дифференциальных уравнений эллиптического или параболического типов вносит наличие в ней нелинейных условий сопряжения (65), описывающих радиационный теплообмен между границами \!>,/ и ¿у подобластей расчетной области П. Рассмотрим некоторые пути учета условий радиационного теплообмена (65) при выполнении конечноэлементной аппроксимации краевой задачи (61)-(65). ■ Будем считать, что поверхности и между которыми задан теплообмен излучением, потги параллельны друг другу, и для всех точек, лежащих на может быть установлено взаимно - однозначное соответствие с точками, лежащими на Б;/, Перепишем уравнения (65) в виде

8и = (66)

X дп

где

.П^е^С/)^^^!,)- ,, (67)

Уравнение (66) может быть учтено в эквивалентной вариационной формулировке для краевой задачи (61)-(64), (66) аналогично учету в ней краевого условия третьего рода (64). В результате эквивалентная вариационная формулировка примет вид „ 81/,

|рС:(, ~ Эг/П + рте! и £гас! 9 сЮ + |р Г/ЭсД" +

£1

|п(г/ - (7)ЗбК = |о!Ш + |/Э</П

° " _ (68) +

и-',

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

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

ментов, имеющих параллельные друг другу основания (например, параллелепипедов, призм, некоторых шестигранников). В этом случае последнее слагаемое в левой части уравнения (68) может быть аппроксимировано на конечном элементе, одно основание которого лежит на поверхности а другое - на соответствующей ей поверхности S,,. Тогда доя элемента поверхности S0 может быть построена локальная матрица, точно такой же размерности, как и для любого конечного элемента iî™ CÎÎ, и поэтому при сборке глобальной матрицы можно использовать одну и ту же процедуру как при обработке конечного элемента П"', так и при обработке одного и того же элемента поверхностей \ и S,,.

Менее удобным вариационное уравнение (68) становится в том случае, когда для его аппроксимации в качестве конечных элементов используются тетраэдры, хотя с точки зрения наиболее точного описания геометрии тетраэдр является одним из лучших конечных элементов. В этой ситуации локальные матрицы для элемента поверхности Su и для тетраэдрального конечного элемента £У" будут иметь разную размерность. Кроме того, разбиение поверхностей .S'y и SJt на треугольники должно быть абсолютно идентичным, что в некоторых ситуациях может существенно ограничивать возможности оптимального размещения узлов в тетраэдральных сетках. Обратим также внимание на то, что при использовании вариационного уравнения (68) тештообмен излучением из каждой точки поверхности SJt идет только в одну точку поверхности т.е. вариационная задача с использованием уравнения (68) описывает излучение только "из точки в точку" для тех точек поверхностей и SJt, между которыми установлено взаимно-однозначное соответствие, и поэтому в данной постановке невозможно провести какое-либо исследование влиянияjia изучаемое температурное поле других типов излучения, например, с рассеиванием. Определенные сложности в этом подходе возникают и в ситуациях, когда поверхности Ли и SJt имеют несколько различающиеся площади. В этом случае для сохранения баланса тепловой энергии необходимо вводить специальные множители в подинтегральные выражения последнего слагаемого левой части уравнения (68) для различных границ S,j подобластей £1,.

Рассмотрим другой подход к конечноэлементному моделированию тепловых нолей с учетом радиационного теплообмена. Этот подход базируется на том, что температура на каждой из поверхностей S:J меняется гораздо слабее, чем между поверхностями Sy и S.,.

Заполним пространство между поверхностями S,t и SJt подобластями В канедой подобласти определим "псевдотеплопроводиость" X по формуле

). = aoE8({/j,t+t/j5. .), (69)

где под Î7L и Î/L понимаются значения температуры U в точках поверхно-

стей и Л}/, ближайших к той точке, в которой определяется значение А., а 8 -расстояние между этими точками, которое может изменяться при движении но поверхностям X,, или Л},, но остается неизменным при движении от 5,/ к Л',, например, вдоль нормали к этим поверхностям. Рассмотрим краевую задачу для уравнения

рСр~-йп[краАи)±/, (70)

определенного в области являющейся объединением исходной области £1 и всех подобластей П,,, со значениями р = р, Ср = Ср, Х-Х, / = / внутри О и р = 0, С =0, X = X, / = 0 (см. (69)) внутри П,у и с краевыми условиями (62)-(64). Нетрудно убедиться, что если на всех поверхностях решение (/ краевой задачи (70), (62)-(64) принимает постоянные значения или все расстояния 8 между поверхностями или очень малы, то функция (У является решением исходной краевой задачи (61)-(65). Если же в любой из подобластей С17 какое-либо из сформулированных выше условий, а именно постоянство и на или малость толщины 8 подобласти не выполняется, то Пеней коэффициент теплопроводности X может быть задан в виде тензора Ь со значением /. • = X (величина X. определяется соотношением (69)) для нормали п к поверхности и значениями /. ■ т = 0 для всех направлений т, перпендикулярных п. Таким образом, введение "псевдотеплопроводности" в виде скаляра X или тензора I позволяет вместо условий сопряжения (65) использовать стандартное уравнение диффузии тепла внутри подобластей Ду со скалярным коэффициентом теплопроводности или с коэффициентом теплопроводности в виде тензора. При этом баланс тепловой энергии сохраняется автоматически, в связи с чем отпадает необходимость его контроля, а исследователь получает возможность изучить влияние рассеянного излучения за счет задания тензора I с ненулевыми значениями /. • х, которые должны лежать в диапазоне от 0 до X.

Подход, основанный на решении краевой задачи для уравнения (70), целесообразен потому, что при его использовании не возникает принципиальных трудностей, связанных с применением тетраэдров для конечноэлементнон аппроксимации этой краевой задачи. Кроме того, не накладывается никаких ограничений и на разбиение поверхностей Яу, сетки на поверхностях и могут быть различными. При этом нестандартной остается только процедура вычисления коэффициента X на каждом конечном элементе из подобластей £2„. Рассмотрим одну из возможных реализаций этой процедуры.

Вначале отметим, что при построении тетраэдральной сетки в расчетной области О внутри подобластей не должно быть ни одной внутренней точки, т.е. во всех подобластях Ц, должен быть только один слой тетраэдров с верши-

нами, лежащими либо на поверхности либо на поверхности Л},.

Для вычисления коэффициента X на каждом тетраэдре из подобласти П/, необходимо иметь возможность получить соответствующее этому тетраэдру значение и на поверхно-1 сти Л*,, и значение и на поверх-' поста Эти значения могут быть легко получены осреднением значений функций С/ в вершинах тетраэдра, принадлежащих одной и той же поверхности. Таким образом, лля-вычисления коэффициента X вполне достаточно уметь определять, какие из вершин тетраэдра, принадлежащего подобласти Ц¡, лежат на поверхности Бу, а какие на поверхности Бу.

Рис. 1. Криволинейный дипольный магнит

2.4. Расчеты трехмерных нелинейных магнитных и тепловых полей в реальных Конструкциях

2.4.1. Оценка точности расчета магнитного поля для 45 - градусного криволинейного днполышго магнита

Точность вычислений магнитного поля проиллюстрируем для реального магнита [33, 34].

Это 45-градусный криволинейный дипольный магнит с радиусом кривизны 1.12 м. и сложной формой фасок. На рис.1 показан его общий вид Разбиение на тетраэдры проводилось с помощью препроцессора пакета МАЯТАС [33, 35] Изображенная на рис. 2 тетраэдральная сетка, содержит около 16000 узлов. Расчеты проводились как на этой сетке, т ак и на более подробной сетке, . Р|,с- 2 Тетраэдральная сетка для моделирования

криволинейного днпольного магнита

содержащей около 64000 узлов. Кроме того, после изготовления магнита были проведены физические измерения поля. На рис.3 приведены результаты сравнения измерений вертикальной составляющей индукции магнитного поля в среднем сечении магнита с результатами расчета на двух сетках. Пунктирной линией на нем обозначены результаты расчета на сетке с 16000 узлов, сплошной линией - на более подробной сетке (64000 узлов), символами "0" - результаты измерений. Абсолютная погрешность измерений индукции магнитного поля не превышает Ю-4 Т. На рис.4 показаны отклонения вычисленного поля от результатов измерения в процентах к значению поля в центре магнита.

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

П(хдатьная ююрлиан (см) Продольная координата (см)

Рис. З Вертикальная составляющая индукции Рис. 4 Отклонения вычисленного поля от магнитного поля в среднем сечении магнита результатов измерения в процентах

к значению поля в центре магнита

2.4.2. Исследование сходимости итерационных процессов учета зависимости р(Я) при решении различных трехмерных задач магнитостатики

с«

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

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

?= 0д + (1-е)?\ (71)

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

|| л л I

\F~Aq,!

Л—-А ^ е, , (72)

, Итерационный прогдесс по нелинейностям прекращался тогда, когда выполнялось неравенство:

ер]

' <е,

(73)

где критерии Б/ выхода из итерационного процесса по нелинейностям во всех проведенных исследованиях принимался равным 104 (под д,да,Г',Л в формулах (71)-(73) понимаются те же векторы и матрицы, что и в п.2.2.3).

Рис. 5. Общий вид вигглера.

Рис. 6. Общий вид андулятора

В таблицах 1-19 приведены результаты исследования, показывающие зависимость от параметров х> 0. ег и от основных параметров решаемой задачи количества итераций по нелинейностям и суммарного числа итераций при решении МСГ-подобным методом линеаризованных систем конечноэлементных уравнений (49) на всех итерациях по нелинейностям.

,г , , Таблица 2. Суммарное число итерации

Табл. I показывает зависимость 1 ' 1

количества итераций по нелинейностям от параметров х и 0 при с,--Ю-1 и суммарном токе в обмотках /=10кА для задачи расчета магнитного поля в 45-градусном криволинейном диполыюм магните. Ток в обмотках величиной ЮкА не приводит к значительному насыщению железа в этом устройстве, и

поэтому количество итераций по пели- Таблица 3. Число итераций по нелнн.

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

Результаты, приведенные в табл.З-4, соответствуют той же задаче (для 45-градусного Криволинейного дипольного Таблица 4. Суммарное число итераций магнита), но при суммарном токе в об- (дипольный магнит. ег=10Л /=20кА)

мотках /=20кА. Такое увеличение тока вызывает существенное насыщение железа во многих местах конструкции, что приводит к заметному росту числа итераций по нелинейностям. Область значений оптимальных параметров X 11 в более стабильна: имеет место незначительное смещение в область увеличения параметра 9.

[3 табл.5-6 представлены результаты расчета магнитного поля в том же 45-градусном диполыюм магните при значеннн тока /=20кА н уровне предельной относительной невязки конечноэлементных СЛАУ £,.=0.05. Столь значительное

Таблица 1. Число итераций по нелин.

(дипольный магнит. е,~|0 4, /=10кА)

0.4 0.5 0 6 0.7 0.8 0.9 • 1.0

0.4 _ — 10 10 9 9 9

(Г 5 10 У 9 8 7 К 9

0 6 9 Я 7 7 7 8 —

0.7 Я 7 7 7 7 9 —

08 7 7 8 9 -- _ —

0,4 9 10 _ — — — —

1.0

(дипольный магнит, £,-10 , /=)0кА)

«V 0.4 0 5 (16 0.7 0 8 0.9 10

0 4 _ _ 1152 1140 1019 975 1031

«5 1155 1037 (024 902 755 916 999

Ой 1014 91 1 785 763 801 932 —

0 7 НИ 753 764 Ш 796 990 —

ОХ 794 748 862 979 _ - —

0 9 Н)Г>| 1104

1 (I

л

(дипольный магнит, сг=10~ , /=20кА)

0\у 0 4 05 06 0.7 0 8 0 9 1.0

0 3 34 29

»4 . — _ — 33 29 26 23

(15 _ 33 30 27 24 21 19

0 6 30 28 25 23 21 18 —

(17 26 24 22 20 (X 16 —

0 8 23 21 19 18 16 _ —

0 9 21 19 1» 16 !5 _ —

1 0 19 Ж 16 15 — — —

0\у 0 4 0.5 0.6 0.7 0.8 0 9 1 0

03 2469 1937

04 — — - 2653 2242 18114 1698

0 5 — 2709 2421 2115 1810 1491 1452

Об 2427 2240 1979 1736 1488 1314 —

0.7 2089 1*82 1689 1499 1320 1249 —

0.8 1813 1631 1475 1320 1520 — _

0 9 1615 1456 1627 1270 1165 — —

1 0 1474 2149 1372 1191 — — —

Таблица 5. Число итераций по нелин

№1С 04 0 5 0 6 0 7 0 8 0.9 10

0* _ 40 36 33 29 26 23

0 5 36 33 30 27 24 21 18

06 30 28 25 23 21 18 —

0.7 26 24 22 20 18 16 —

0.8 23 21 19 18 16 — —

09 21 19 18 16 15 _ —

10 19 18 16 15 — — —

Таблица 6. Суммарное число итераций

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

Табл.7-8 показывают зависимость от параметров х и 0 количества итераций по нешшейностям и суммарного числа итераций в процессе решения линеаризованных систем конечноэле-ментных уравнений для задачи вычисления магнитного поля в вигглере при значении суммарного тока в обмотках /=16кА и параметра с^КГ3. Результаты, приведенные табл.9-10, соответствуют той же задаче, но с измененным значением параметра сг=0.01, а результаты, приведенные в табл.11-12 соответствуют той же задаче с теми же параметрами /=16кА и ег=0.01, но при использовании для конечноэлементной аппроксимации другой, более подробной тетраэдральной сетки (число узлов в ней увеличено на порядок).

«V 0.4 0.5 0.6 0.7 0 8 0.9 1.0

0.4 — 670 666 618 550 452 606

0.5 573 558 554 519 473 439 505

0.6 489 455 481 446 439 373 —

<1.7 463 377 438 403 463 410 —

08 419 395 336 425 425 — —

09 413 374 390 487 390 _ —

1.0 396 441 444 385 — — —

Таблица 7. Число итераций по нелин

Таблица 9. Число итераций по нелин.

0.4 0.5 1 0.6 0.7 0.8 09 1 0

0 4 40 37

0.5 _ — — — .16 32 30

0 6 _ - — 35 30 27 25

0.7 — — 16 31 26 . 24 22

0 8 _ 36 31 26 23 21 —

0 9 37 33 28 24 21 19 —

1 0 34 30 26 22 19 — —

о\х 04 0 5 06 0.7 0.8 0.9 1.0

04 — — _ — - 40 37

0 5 — _ _ — 36 32 30

0.6 — — _ 35 30 27 27

0.7 — _ 36 31 26 23 21

0.8 — 36 31 26 2.3 21 _

0.9 37 33 28 24 21 19 —

1.0 34 30 26 22 19 — —

Таблица 8. Суммарное число итераций

(виггпер, е,= 10 /=16кА)

н 0.4 0.5 06 0.7 0.8 0.9 1 0

0,4 1451 1380

0 5 _ — — — 1359 1141 1110

0 6 — _ _ 1294 ИЗО 958 915

0 7 — — 1302 1136 971 851 813

0 8 — 1263 1115 952 855 9.16 _

0 9 1265 1151 100.3 874 780 664 —

1 (1 1157 1041 919 795 699 — —

Таблица 10. Суммарное число итераций

(вигглер, е,=0 01, /=16кА)

0\/ 0.4 0.5 0.6 0.7 08 0.9 1.0

0.4 1038 873

05 — — — 996 816 703

06 _ -- _ 951 828 691 И 1

0.7 — — 958 838 712 588 504

08 — 940 821 697 626 539 —

0,9 940 856 735 636 572 490 _

10 859 773 679 582 509 — —

По результатам, приведенным в табл.7-12, можно сделать вывод, что и для задачи расчета магнитного поля в вигглере увеличение параметра ег до значе-

ния 0.01 не повлияло на количество итераций по нелинейностям, но заметно уменьшило суммарное число итераций в процессах решения линеаризованных систем конечноэлементных уравнений. Изменение же конечноэлементной сетки очень слабо отразилось на количестве итераций по нелинейностям и несколько повлияло на суммарное число итераций в процессах решения линеаризованных систем конечноэлементных уравнений, что является вполне естественной реакцией на увеличение размерности решаемых СЛАУ, но область оптимальных значений параметров х и 0 изменилась довольно мало, переместившись в среднем на 0.1 в сторону уменьшения параметра 9.

Результаты, приведенные в табл.13-14, соответствуют задаче вычисления магнитного поля в вигглере, но при увеличенном вдвое суммарном токе в обмотках и е,=0.01. Уменьшение величины ег до значения 10~4 не изменило количество итераций по нелинейностям, но практически вдвое увеличило суммарное число итераций при решении систем линеаризованных конечноэлементных уравнений (см. табл.15).

Таблица 11. Число итераций по нелин.

0\у 0 4 0 5 Ой 07 ОХ (19 1 0

0.4 37 34

05 — — — — 37 30 27

0.6 _ — — 18 .11 26 —

0 7 — _ ЗУ 33 27 22 —

ОХ — 34 29 24 — —.

0 9 40 35 31 26 21 — _

1 0 36 32 28 23 — — —

Таблица 13. Число итераций по нелин (вигглер, Ег= 0 01, /~32кА)

о\7 04 0.5 0.6 0.7 0.8 0.9 1.0

0 6 — — — _ 40 36 34

0 7 — — 38 36 34 32 —

0 8 38 35 34 32 29 28 _

0.9 34 32 31 29 27 25 —

1 0 31 29 28 25 25 23 —

Таблица 14. Суммарное число итераций

Таблица 12 Суммарное число итераций (вигглер, подр.сетка £,-10 /=16кА)

(>\у 0 4 0 5 0 6 0,7 0 8 0.9 10

0 4 — — — _ — 1206 1163.

0 5 — _ _ _ 1144 967 935

0 6 — — — (112 956 836 —

0.7 —' — 1101 966 827 924 —

0 8 — 1070 958 853 735 — —

0 9 1082 95Х 869 764 645 — —

10 975 К82 784 679 — — —

0\/ 0.4 0.5 0.6 0.7 0.8 0 9 1.0

0.6 — — — — 957 869 852

07 — — 909 859 807 769 —

08 878 827 804 758 696 674 —

0.9 7Н1_ 751 725 680 640 603 —

1.0 70У 674 652 592 586 547 —

Таблица 1?. Суммарное число итераций

в\г 04 0.5 0.6 0.7 0.8 09 1.0

0 6 — — — — 1799 1657 1641

0.7 _ — 1670 1601 1559 1466 —

0.8 1627 1525 1490 1415 1295 1278 —

0.9 1453 1390 1352 1275 1196 1136 _

10 1321 1255 1217 1102 ПОО Ш42 —

В табл. 16-19 приведены данпые, соответствующие расчету магнитного поля в андуляторе. Табл.16 и табл.18 показывают зависимость от параметров

X и 0 количества итераций по нелинейностям при значении параметра £,, равном 10"4 и 0.05 соответственно. Необходимо отметить, что эта задача является единственной из всех рассмотренных, в которой увеличение значения параметра ег от Ю-1 до 0.05 привело к некоторому росту числа итераций по пе.'ппичшо-

стям. При этом количество итераций по нелинейностям при значении параметра е,.=1(Г3 не отличалось от количества итераций по нелинейностям при значении параметра ег=10"4 для всех значений параметров х и 0 и при ограничении максимального числа итераций по нелинейностям значением 35. С другой стороны, суммарное число итераций в процессе решения систем линеаризованных ко-нечноэлементных уравнений при значениях параметров х " 0, близких к оптимальным, для значения ег=0.05 оказывается почти в четыре раза меньше, чем для е,.=1(см. табл.17 и табл.19 ), и более чем в три раза меньше, чем для ег=1(Г3. Счедовательно, оптимальным среди значений параметра бг=10~*, КГ3, 0.05 по суммарным вычислительным затратам безусловно является значение е,=0.05.

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

Таблица 16. Число итераций по нешш. Таблица 18, Число итераций по недин.

(андулятор, сг"ЮЛ /~19.2кА) _(андулятор, е,=0.05, /=19.2кА)

0У/, 0.4 0.5 0 6 0.7 0.8 0.9 10 о\7, 0.5 0.6 0.7 ОХ 0 9 1.0

0.7 — — — _ — _ — 0.7 __ — — - — —

0 8 — _ — _ 35 33 32 0 8 — _ — _ 34 32

0 4 _ — 35 33 31 30 29 0.9 — • _ — 34 31 29

1 0 35 33 32 30 28 27 27 1.0 — 34 31 29 27

Таблица 17. Суммарное число итераций (андулятор, Е,=10 4, /=19.2кА) Таблица 19. Суммарное число итераций (андулятор, к.=0.05, /=19.2кА)

0 4 0 5 об 0 7 0.8 0.9 1.0 оу 0 5 о.о 0.7 0 8 0.9 1 0

07 _ — — — _ — 0.7 — — — — — —

ох — — — — 5151 4936 4912 0.8 — _ — 1283 1259

0.9 — — 4951 4745 4550 4481 4415 0.9 _ — — 1309 1217 1158

1 0 4784 4591 4513 4298 4097 4043 4143 1.0 — — 1274 1190 11.19 1087

При проведении вычислительных исследований была выявлена еще одна особенность поведения итерационных процессов данного класса нелинейных задач магнитостатики. Зависимости количества итераций по нелинейностям суммарного числа итераций при решении систем линеаризованных конечно-элементных уравнений от параметров х и 0 достаточно гладкая для всех выполненных расчетов. Только при расчете магнитного поля в 45-градусном криволинейном дипольном магните (суммарный ток в обмотках /=20кА) для двух наборов параметров ("// 0.5, 0=1 и /~0.6, 0-0.9) наблюдалось локальное возрастание суммарного количества итераций при решении конечноэлементных СЛАУ, причем это происходило при всех включенных в исследование значениях параметров г.г и /, хотя в количестве -итераций по нелинейностям подобных выбро-

:ов для этих знамений параметров не было. -

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

Для того, чтобы более детально изучить особенности поведения итерационных процессов в окрестности оптимального значения х при фиксированном значении параметра 0, рассмотрим зависимости относ^ельнон погрешности итерационного процесса 5 = - от номера п итерации по нели-

нейностям для некоторых значений параметров % и 0 при решении различных задач.

На рис.7-12 приведены графики зависимости 6(л) при различных значениях параметров 0, х и с, для различных технических устройств.

По результатам вычислительных экспериментов, представленных на рис.712, можно сделать следующие выводы. По поведению кривой 5(п) практически всегда можно определить, с какой стороны от оптимальных значений параметров х и 0 выбраны текущие значения этих параметров. Так, если функция 5(/i) монотонно убывает с ростом п, то с большой долей уверенности можно утверждать, что либо оптимальное значение параметра х больше текущего, либо оптимальное значение параметра 0 больше текущего. Если же функция 5(и) осциллирует, то или оптимальное значение параметра % меньше текущего, или оптимальное значение параметра 0 меньше текущего. Исходя из этого, может быть рекомендована следующая стратегия ускорения сходимости итерационного процесса по нелинейностям. Выбираются значения параметров % и 0 меньше единицы (например, 0.6). Выполняется около 10 итераций по нелинейностям. Если функция 8(п) на последних пяти итерациях монотонно убывает, то значение б увеличивается на 0.1 и делаются очередные 5-6 итераций. Процесс увеличения параметра 0 на 0.1 после каждых 5-6 итераций продолжается до тех пор, пока этот параметр не достигнет значения 1 или функция o(/i) не перестанет монотонно убывать. Если после увеличения параметра 0 на 0.1 функция S(/j)

перестает монотонно убывать, то необходимо вернуться к предыдущему его значению и продолжить итерационный процесс с этим значением 0 до получения решения. При этом одновременно целесообразно контролировать осцилляции функции 5(/i), и в случае их появления уменьшить 6, например на 0.05. Если же и при значении 0=1 функция 5(н) продолжает монотонно убывать в течение не менее чем пяти-шести итераций, то необходимо увеличить на 0.1 параметр Х- Такое увеличение х на 0.1 можно продолжать до тех пор, пока либо не будет достигнуто значение х=1, либо функция 5(и) не начнет осциллировать. В случае появления осцилляций функции 5(п) необходимо вернуться к предыдущему значению

8 1е-1

1е-2

1е-3

1е-4

1е-5

......ÍNJ........к

Л.У:......;.........

i » » i

5 10 15 20 25

а) х = 0.5, кривая 1:6 = 0.8, кривая 2: 0 = 0.9, кривая 3:0 = 10.

5 г ■

1е-1

i »r«U-

1е-2 ' 1i \*

1е-3 1е-4 1е-5

......-

Л \л8- :

- V ;V :

v

5 10 15 20 25

6) х = 0.6, кривая 1:0 = 0.8, кривая 2: 0 = 0.9, кривая 3:6= 1.0.

8 1е-1 1е-2 1е-3 1е-4 1е-5

! '! ! i

\ i : Ф

Ф \

i l \

10

15 20

25

б 1е-1 1е-2 1е-3 1е-4 1е-5

14.

......

г

/•*

...

......

vVW

(2Ц!

®

V

-L

К

ю

15

20 25

в) х = 0.7, кривая I: 0 = 0.7, кривая 2: 0 = 8.0, г)х = 0.8, кривая 1: 0 = 0.7, кривая 2: 0 = 0.8, кривая 3: 0 = 0.9, кривая 4:9 = 1.0. кривая 3: 0 = 0.9, кривая 4:6= 1.0.

5

5

Рис. 7 Графики зависимости 6(/i) для 45-градусного криволинейного ди-польного магнита, ег=Т0Л 1= ЮкА

1е-1 1е-2 1е-3 1 е-4 1е-5

I • _

12\

5 10 15 20 25 "

а) X = 0 6. кривая 1:8 = 0.8, кривая 2:9 = 0.9, кривая 3 8= 1.0.

6 1е-1

1в-2

1е-3

1е-4

1е-5

........!~

........"... ф .......з '' ■ Г ; (2)

10 15 20 25 «

б)х = 0.8, кривая 1:0 = 0.7, кривая 2: 6 = 0 8, кривая 3:0 = 02, кривая 4. 6= 10

6 1е-1 1е-2 1е-3 1е-4 1е-5

Лч -

л;/

5 10 15 20 25 "

в) х = 0 9, в

кривая 1: 6 = 0.6, кривая 2: Э = 0.7, кривая 3:8 = 0.8, кривая 4: 9 = 0.9.

1е-1 1е-2 1е-3 1 е-4 1е-5

ш

\

-Ф

5 10 15 20 25 г)Х = Ю, кривая 1:0 = 0.4, кривая 2:0 = 0.5, кривая 3:9 = 0.6, кривая 4:0 = 0.7.

Рис. 8 Графики зависимости 5(я) для 45-градусного криволинейного дипольио-го магнита, е,^ I О 4, У=20кА

5 1е-1 1е-2 |е-3 1е-4 1е-5

4%

Ш".....Г"""" "Г"

: ! / X,

б

1е-1

1е-2

1е-3

1е-4

1е-5

5 10 15 20 25

а)Х = 0.6, кривая 1:0 = 0.7, кривая 2:8 = 0.8, кривая 3:0 = 0.9, кривая 4:0=1.0.

1

I '

(5

25

5 10 15 20

б)х = 0.8,

кривая 10 = 0.7, кривая 2: в = 0.8, кривая 3: 9 = 0.9, кривая 4: 0 = 1.0.

п

3

1е-1 1е-2 1е-3 1е-4 1е-5

Ж

...... ...<3,1.

А......

ф

5

25

10 15 20

в)х = 0.9, кривая 1:6 = 0.6, кривая 2:9, = 0.7, кривая 3.0 = 0.8, кривая 4:0 = 0.9.

б -''Л

1е-1 1е-2

1е-3 1е-4 1е-5

г ,

• --V

Л, '

X......I.....я

.....

; 12, :

5

25

Ю 15 20 г)* =1.0, кривая 1:6 = 0.4, кривая 2: 9 = 0.5, кривая 3: 8 = 0.6, кривая 4:0 = 0.7.

Рис. 9 Графики зависимости 6(и) для 45-градусного криволинейного диполыш го магнита, ег=0.05, /=20кА

8 е-) е-2 е-3 е-4

6-5

-4)

ка

5 10 15 20 25 " а)х=0.6, кривая 1: 9 = 0.7, кривая 2: 8 = 0.8, кривая 3:0 = 09, кривая 4. 6=1.0.

8 1а-1 1е-2 1е-3 1е-4 1e-S

! | i

'S. ]...! .

..........j...... ........: ; ; ® ; ! i ! <i> @ L.......

5 10 15 20 25

б) X = 0.8, кривая 1:6 = 0.7, кривая 2: S = 0.8, кривая 3: 8 = 0.9, кривая 4:9 = 10.

S le-1 le-2j 1 е-з| 1е-4| ■ 1e-5L—

■V.

Vf

di:

\

(âi • >2i

s

10 15 20 25

в) X = 0.9, кривая 1: 0 = 0.7, кривая 2: 0 = 0.8, кривая 3:8 = 0.9, кривая 4: 8 = 1.0.

Ъ

6 1е-1 1е-2 1е-3 1е-4 1е-5

" 3 . . "//. ].....; ...

Л^Ь i

.........N-Л' """î......... ч» ; ® .........•..........•.

.....

5 10 15 20 25 Г)Х=10, кривая J : в = 0.7, кривая 2: Ô = 0.8, кривая 3: 0 = 0.9.

Рис. 10 Графики зависимости 8(л) для вигглера, ег=0.01, /=16кА

5 1e-1 1e-2 1e-3 1e-4 1e-5

ц i ; fi. i' ;. •

M l1^

i ; [ (S)

5

25

10 15 20

a) X = 0.6, кривая 1:0 = 0.7, кривая 2:0 = 0.8, кривая 3: 9 = 0.9, кривая 4: 0 = 10.

8 1е-1 1е-2 1е-3 le-4 1е-5

»4

4

Ф

%щп.............;

4®.....i.........."Му«;:;*

là è

5

25

10 15 20 6) X = 0.8, кривая 1:0 = 0.7, кривая 2: 8 = 0.8, кривая 3:0 = 0.9, кривая 4. 6 = 1.0.

5.

1е-1 1е-2 1е-3 1е-4 1е~5

•S

i .....

а

; (3) i <g\

5

25

10 15 20

в) X = 0.9, кривая 1:9 = 0.6, кривая 2:0 = 0.7, кривая 3:6 = 0.8, кривая 4:0 = 0.9.

п

5 1е-1

1е-2

1е-3

1е-4

1е-5

, ..... ,

...

: Ф \ ; .......■

: (2) ' .

5 10 15 20 25 r)x=l.O, кривая 1:9 = 0.5, кривая 2 9 = 0 6, кривая 3: 0 = 0.7, кривая 4: 0 = 0.8.

Рис. 11 Графики зависимости 5(п) для витглера, Ег=0.01, /= 32кА

Ф

......

;

15 20 25 "

1...........1......1

5 10

а)х=-0.6, кривая 1:9 = 0.9, кривая 2; 0 = 1.0.

S 1е-1 1е-2 1е-3 1в-4 1е-5

..........

11 г >

s п

б) X = 0.8, кривая 1:0 = О.д, кривая 2: 0 = 1.0.

•s

Ф

■■ .....

:

5 10 15 20 25 в)Х = 0.9, кривая 1:6 = 0.9, кривая 2:9 = 1.0,

5 10 15 20 25

Г)Х = 1Д кривая 1:9 = 0.8, кривая 2: 6 = 0.9, кривая 3: 0 = 1.0,

Рис. 12 Графики зависимости 5(л) для андулятора, ег=0.05,1=19.2кА

При решении одной и той же задачи на разных сетках или при решени класса близких задач, различающихся лишь небольшими изменениями геомс рии (например, при оптимизации фасок) или Значений токов в обмотках, подбс оптимальных параметров 0 и % достаточно выполнить при решении первой з; дачи из этого класса. Все остальные задачи вполне могут быть рассчитаны использованием полученных при решении первой задачи параметров 0 и ; Критерием применимости значений этих параметров служит появление осци. ляций функции б (и). Если при решении очередной задачи из одного класса о цилляшш функции 6(л) будут иметь место, то необходимо уменьшить значен! параметра 9 или %.

Приведенные в таблицах 1-19 и на рис.7-12 данные отражают лишь н большую, но наиболее характерную часть результатов проведенных болыш серий вычислительных экспериментов, целью которых было изучение област< локализации оптимальных параметров х, 0 и ег в зависимости от типов технич ских устройств, геометрии их ферромагнитных сердечников, величин токов обмотках, степени дискретизации расчетной области. Полученные результат показывают, что области локализации оптимальных параметров х, 0 и ег завис в основном лишь от типа технического устройства и почти не меняются при » менениях геометрии ферромагнитных сердечников, связанных с оптймизаци отдельных конструктивных элементов этих технических устройств. Даже суп ственное изменение токов в обмотках, приводящее к значительному изменен! коэффициента магнитной проницаемости (за счет насыщения отдельных э; ментов конструкции), довольно слабо влияет на области локализации от мальных параметров 0 и ег. Практически не меняются оптимальные значен параметров х,0 и ег и при изменении дискретизации расчетной области. По.г ченные закономерности поведения зазнсимостей 8(п) при изменении пара ли ров х, 0 и ег позволили сформулировать конструктивную стратегию ускорен итерационного процесса по нелинейностям.

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

2.4.3. Расчет нестационарного трехмерного поля в тепловом экране

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

Цель исследования заключалась в изучении влияния количества пластин в ране, свойств материалов, из которых изготовлены шпильки и упоры, а также эффициента черноты поверхностей пластин на величину теплового потока, оникаюшего от излучателя тепла в резервуар с жидким азотом. Требовалось <же оценить время выхода процесса теплообмена на стационарный режим, [я этого были проведены расчеты более чем 50 вариантов нестационарных мпературных полей при различных значениях коэффициентов уравнения (61), повий сопряжения (65) и краевых условий (62). При этом температурное поле числялось как решение краевой задачи для уравнения (70), в котором козф-циент теплопроводности X в пространстве между поверхностями пластин ределялся в виде тензора /. с единственным ненулевым элементом X, распо-женным на диагонали его третьей строки. Значение X определялось по фор-ле (69), в которой параметр 5 полагался равным расстоянию между пласти-«и экрана.

На рис.14 приведены графики зависимостей от времени / тепловых потоков, проникающих в резервуар с жидким азотом для варианта расчета с закреплением теплового экрана на резервуаре с жидким азотом при температуре поверхности излучателя тепла, равной 500 К, и температуре поверхности резервуара с жидким азотом, равной 78 К. Начальная температура теплового экрана в этом варианте расчета было принята равной 300 К. Кривая 1 на рис.14 показывает график зависимости от I теплового потока, проникающего в резервуар с жидким азотом за счет излучения тепла с пластины экрана, кривая 2 - за счет теплопроводности упоров, а кривая 3 показывает зависимость от I суммарного теплового потока, проникающего в резервуар с жидким азотом от теплового экрана. Приведенные на рис.14 графики зависимостей тепловых потоков позволяют судить не только о количестве тепла, проникающего за определенный интервал времени в резервуар с жидким азотом, и о времени выхода процесса теплообмена на стационарный режим, но и о том, какая часть суммарного теплового потока с теплового экрана попадает в резервуар с жидким азотом за счет излучения, а какая - за счет теплопроводности упоров. Полученные результаты были использованы при оптимизации конструкции тепловых экранов.

3. Основные результаты и выводы

1. Разработаны алгоритмы математического моделирования теплового состояния турбогенераторов традиционного и нетрадиционного исполнений.

2. Предложена математическая модель радиационного теплообмена.

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

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

ный алгоритм, учитывающий разрыв потенциалов.

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

6. Разработанные модели и алгоритмы использовались при проектировании различных технических устройств таких, как турбогенераторы со сверхпроводящей обмоткой ротора (К'ГГ-20), в результате чего была предложена модификация системы каналов в теплообменниках (изменение числа заходов винтового канала, место посадки теплового экрана). Были проведены расчеты для учета влияния консзрукторско-технологического исполнения (геометрия ферромагнитного сердечника, фаски) криволинейного 45-градусного дипольного магнита на магнитное поле в рабочей зоне; выполнен сравнительный анализ различных реализаций сверхпроводящих микроандуляторов, разрабатываемых в Brookhaven National Laboratory (USA); исследованы магнитные поля в рабочей зоне и поля рассеяния мультиполыюго вигглера и андулятора, разработанных в ИЯФ СО РАН (Новосибирск).

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

1. Кузнецов Ю.А., Соломахин В.И., Шурина Э.П. К зг&аче расчета температурных полей. // Разностные и вариационные методы. 1977. -Новосибирск,-С. 155 - 160.

2. Соломахин В.И., Шурина Э.П. Анализ теплового поля дефектов активной стали статоров турбогенераторов. //: Вопросы надежности, автоматического контродя и защиты мощных синхронных генераторов. - 1978. Ленинград : ВНИИэлектромаш. - С.29 - 41.

3. Глебов И.А., Кыков В.М., Большаков Л.П., Шурина Э.П. Распознавание очагов нагрева в массивных ферромагнитных телах.// Дефектоскопия. -1981. № II.-С. 60-70.

4. Шурина Э.П. Математическое обеспечение вычислительных экспериментов в краевых задачах. //Тезисы докл. Всесоюзной конференции : Автоматизация научных исследований на основе" применения ЭВМ. -1981.Новосибирск. - С. 169.

5. Карчов Д.С., Шурина Э.П. Решение осесимметричных краевых задач методом конечных элементов.// Применение ЭВМ в оптимальном планировании и проектировании. 1981. Новосибирск: НЭТИ. - С.71 - 77.

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

7. Соломахин В.И., Карчов Д.С., Свидченков М.В., Соловейчик Ю.Г., Шурина Э.П. Организация комплекса программ расчета трехмерных тепловых

полей в составных областях.// Комплексы программ математической фи ки. Материалы VII Всесоюзного семинара по комплексам программ ма матической физики. - Новосибирск.1982. - С.291 - 296.

8. Карчов Д.С., Соловейчик Ю.Г., Шурина Э.П. Один подход к решению дачи газовой динамики в винтовом канале.//Численные методы мехаш сплошной Среды. - 1982. T.I3, № 3. - С.148 - 149.

9. Шурина Э.Г7., Карчов Д.С., Соловейчик Ю.Г.Решение сопряженной задг для ротора КТГ. // Исследование крупных электрических машин. -19 Ленинград : ВНИИэлектромаш. - С.24 - 32.

10. Глебов И.А., Данилевич Я.Б., Беляев С.Н., Брынский Е.А., Журавлев Г. Карчов Д.С., Смолин И.М., Соловейчик Ю.Г., Шурина Э.П. Температ ные поля ротора турбогенератора со сверхпроводящей обмоткой возб) денпя. //Известия АН СССР.Энсргетика и транспорт. - 1982,- № 5. - С. 7 83.

И. Данилевич 51.Б., Карчов Д.С., Соловейчик Ю.Г., Шурина Э.П. Програи ное и математическое обеспечение САПР сверхпроводниковых турбоге раторов. // Автоматизация проектирования и конструирования. М. - 19К Ч. 1.-С.101 -102.

12. Соломахин В.И., Шурина Э.П. Математическое обеспечение тепловых р четов в статорах турбогенераторов. // Исследование генераторов с полн водяным охлаждением. - 1983. Ленинград : ВНИИэлектромаш. С. 73 - 82

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

14. Глебов И.А., Данилевич Я.Б., Шурина Э.П. Теплообмен во вращагоще? криволинейном канале.// Доклады АН СССР. - 1983. Т. 269, № 5. - C.lOi 1098.

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

16. Шурина Э.П. Численное моделирование турбулентных течений жидко* (газа) в каналах. // Автоматизация научных исследований. Тезисы до Всесоюзной школы ИАЭ СОАН СССР. Новосибирск. - 1985. С.174 - 175

17. Соломахин В.И., Шурина Э.П. Алгоритм идентификации тепловых ист ников в статоре турбогенератора. // Математическое обеспечение иссле панин стохастических и детерминированных моделей. - 1986. Новосиби: : НЭТИ. - С.8 - 15.

18 Петров Ю.В., Свидченков М.В., Шурина Э.П. Алгоритмы и результ; численною моделирования системы охлаждения обмоток рогора турбс перл юра.// Техническая электродинамика. - 1987'. № 5. - С.75 - 79.

. Петров Ю.В., Свидчеиков M B., Шурина Э.П. Математические модели системы охлаждения роторов турбогенераторов. Деп. в ВИНИТИ № 2306 -B87.1987.IOc.

. Соломахин В.И., Шурина Э.П. Пакет прикладных программ моделирования теплового состояния статоров мощных турбогенераторов. // Комплексы программ математической физики и архитектура ЭВМ. 1988.Красноярск. - С.324 - 328.

. Карчов Д.С., Шурина Э.П. Моделирование и анализ чувствительности теплового состояния сверхпроводникового турбогенератора. //Техническая электродинамика. - 1988. №1. - С.68 - 74.

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

. Свидчеиков М.В., Соловейчик Ю.Г., Шурина Э.П. Планирование эксперимента и оценивание параметров теплопереноса при турбулентном течении хладагента во вращающихся и неподвижных каналах. //Тезисы докл. VI Всесоюзн. семинара :Обратные задачи и идентификация процессов теплообмена. М.,1988. С.47 -48.

. Соловейчик Ю.Г., Шурина Э.П. Численное моделирование электромагнитных полей в трехмерных составных областях.// Angewandte numerical magnetic field berecliining. - 1989. Chemnitz. №1. P. 179 - 189.

. Данилевич Я.Б., Журавлев Г.С., Свидчеиков М.В., Шурина Э.П. Темпера-фрное поле ротора турбогенератора при многострунной газовой системе охлаждения. // Исследование и вопросы проектирования турбо - и гидрогенераторов. - 1989. Ленинград : ВНИИэлектромаш. - С.53 - 58.

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

. Захаров В.П., Соловейчик Ю.Г., Шурина Э.П. Моделирование гидродинамических процессов и теплопереноса во вращающихся криволинейных каналах. // Численные методы и программное обеспечение, -1990. М.:ОВМ АНСССР.-С.41-55. «

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

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

. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Математическое моделнрова-

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

31. Шурина Э.П., Захаров В.П. Об одном численном алгоритме интегрирования уравнений Навье - Стокса. // Вопросы атомной науки и техники. Сер. Математическое моделирование физических процессов. М. - 1994. Вып. 4,-С. 45 - 52.

32. E.P.Shurina, Y.G. Soloveichik, M.E.Rojak. Three-dimensional fields modeling on irregular mesh using finite elements method. Proceedings of First Asian Computational Fluid Dynamics Conference. Hong Kong. 1995. P. 1225 - 1226.

33. A. Grudiev, M. Rojak, E.Shurina, Yu. Soloveichik, M. Tiunov and P. Vobly. Mathematical Simulation of 3-D Magnetostatic Fields Using the Complex of Programs MASTAC. Abstracts of AMCA - 95. Novosibirsk. NCC Publisher. 1995. P. 131 - 132.

34. E.P.Shurina, Y.G. Soloveichik, M.E.Rojak, P.D.Vobly, A.V.Gmdiev, M.A.Tiunov. Three-dimensional non-linear magnetostatic fields mathematical modeling. In Abstracts of SIBCONVERS'95.Tomsk. Russia. 1995.P.30.

35. M. Rojak, E.Shurina, Yu. Soloveichik and V. Malyshkin. Parallelization of Computer Code MASTAC Three-Dimensional Finite Elements Method Implementing. Proceedings of PaCT - 95. Lecture Notes in Computer Science 964. Parallel Computing Technologies. Springer. Germany. 1995. P. 304 - 313.

36. Шурина Э.П., Соловейчик Ю.Г., Рояк М.Э. Вычислительные схемы решения трехмерных нелинейных магнитостатических задач.//Тезисы докл. Международной конференции : Математические модели и численные^ме-годы механики сплошных сред. - 1996. - Новосибирск; СО РАН. - С. 5 ЯР-

37. М. Rojak/ E.Shurina, Yu. Soloveichik and A. Grudiev, M. Tiunov, P. Vobly. MASTAC - new code for solving three-dimensional non-linear magnetostatic problems. Proceedings of Par - 95 (Particle Accelerator Conference and International Conference of Higli-Energy Accelerators). USA. Dallas. Texas. IEEE. 1996. Vol. 4. P.2291 - 2293.

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