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

Методы и программы для численного решения параболических уравнений развиваются на протяжении многих десятилетий. Большое внимание к уравнениям и системам уравнений данного типа вызвано тем, что они описывают множество процессов в природе, технике, экономике и социологии [1]. Алгоритмы решения параболических уравнений входят в состав многих современных интегрированных систем автоматизированного проектирования, инженерного анализа и подготовки производства [25], например, SDRC, Unigraphics, ADAMS, T-FLEX, ИСПА. Укажем также на современные пакеты прикладных программ, применяемые для исследований в технике и для разработки новых технологий — ANS YS, Р ATRAN, NASTRAN, COSMOS, PHOENICS, N3S, FLUENT, STAR-CD, CFX, NUMECA, и другие. В ИПМ им. M.B. Келдыша РАН и МММ РАН для научных и прикладных исследований были созданы такие комплексы, как ТЕКОН [3] и ACTA [4].

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

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

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

Для решения нелинейных параболических задач широко используются конечно-разностные [1−8] и конечно-элементные вычислительные технологии [9 — 13, 16, 18- 24].

В реализации вышеупомянутых методик можно выделить три основных этапа:

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

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

— обработка и анализ результатов численного моделирования.

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

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

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

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

Для дискретизации сложных областей в работе также используется метод тетраэдризации [26]. Одним из наиболее эффективных методов тетраэдризации является фронтальный метод. Данный метод позволяет управлять формой и размерами элементов в процессе построения сетки.

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

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

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

В зависимости от типа прикладных задач для их численного решения применяют различные конечно-элементные подходы [2, 3, 7, 9, 19, 22, 24, 50]. Задачи, рассматриваемые в данной работе, в общем виде могут быть сформулированы следующим образом. В некоторой трехмерной расчетной области О с кусочно-гладкой границей Г рассматривается система дифференциальных уравнений:

V = -к&аАср где ф — искомая функция,.

Ж — потока,.

Су, к, <2 — соответственно, коэффициенты уравнений и функция плотности источников, зависящие от координат и искомого решения.

На различных участках границы области Г могут задаваться граничные условия следующих типов:

— первого рода ср = Тр ,.

— второго рода д (р!дп = -Ц,.

— третьего рода адср /дп + /3(р = у, а1 + р1 ф 0.

В исходной постановке подразумевается, что неизвестными являются и поток Ж, и функция ф. Для восполнения ф и ^ по значениям в узлах в пределах конечного элемента можно использовать либо элементы одного типа, либо элементы различных типов. Также можно использовать элементы одного типа, но разных порядков аппроксимации, как, например, принято для аппроксимации разных переменных при моделировании уравнения Навье-Стокса [5, 7, 51]. Это фактически приводит к понятию смешанной конечно-элементной формулировки неизвестных функций. Гибридные и смешанные конечные элементы широко применяются для решения задач механики сплошных сред [50]. В разнородной среде с сильно переменными (разрывными) свойствами использование подобного метода позволяет сохранить точность в областях разрыва свойств. Отметим что, в 1980;х в Институте Прикладной Математики им. М. В. Келдыша и в Институте Математического Моделирования РАН под руководством академика A.A. Самарского разрабатывались, так называемые, вариационно-потоковые разностные схемы для решения эллиптических и параболических уравнений [3, 61]. Математическую основу этих схем составляет вариационный принцип Био [40]. Название «вариационно-потоковые» данные схемы получили вследствие того, что, согласно методу Био, решение соответствующей задачи представляется как минимизация некоторого функционала, зависящего как от искомой функции (скаляра) так и от компонент потока, зависящего от градиента искомой функции. Сеточные уравнения, соответствующие вариационно-потоковым схемам для уравнения теплопроводности [3], записываются в терминах компонент тепловых потоков, нормальных к граням ячеек разностной сетки. Вариационно-разностные схемы можно трактовать как разновидность смешанных конечно-элементных аппроксимаций Равьярта-Томаса минимального порядка точности.

Зачастую на практике систему исходных уравнений (1) с двумя неизвестными функциями сводят к уравнению более высокого порядка с одной неизвестной функцией. Из системы двух дифференциальных уравнений (1) путем исключения переменной W получают одно дифференциальное уравнение вида:

С^ = с1Мк. ёгас1(р) + 0, (2) дt и решают его относительно ф.

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

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

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

— автоматическое построение геометрической модели расчетной области и ее дискретизации (генерации сетки) для трехмерных областей сложной формы,.

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

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

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

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

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

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

Согласно нормативно-технической документации (ГОСТ 17−102) дефектом называется каждое отдельное несоответствие продукции требованиям. При применении средств неразрушающего контроля под дефектом понимается любое отклонение измеряемого параметра от проектно-конструкторской документации. По происхождению дефекты подразделяются на:

— производственно-технические, возникающие в процессе проектирования и изготовления изделия, его монтажа, установки;

— эксплуатационные, возникающие в процессе эксплуатации изделия, в результате процессов деградации, а также в результате неправильной эксплуатации изделия.

Поскольку номенклатура изделий весьма обширна, на практике выделяют несколько групп, для которых характерны однотипные дефекты [66]:

— силовые металлоконструкции;

— сосуды, теплообменные аппараты, трубопроводы;

— механизмы и машинное оборудование;

— трубопроводный транспорт, корпуса систем под давлением, парогенераторы, гидросистемы и т. п.;

— контрольно-измерительные приборы и автоматика, оборудование систем управления;

— кабельное оборудование;

— электронное оборудование;

— объекты, содержащие радиоактивные вещества;

— конструкции строительных сооружений.

Выбор метода неразрушающего контроля должен быть основан помимо знания о характере дефекта на таких факторах, как:

— условия работы изделия;

— форма и размеры изделия;

— физические свойства материала изделия;

— условия контроля и наличие подходов к проверяемому объекту;

— технические условия на изделия, содержащие количественные критерии недопустимости дефектов и зачастую нормирующие применение методов контроля на конкретном изделии;

— чувствительность методов диагностики.

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

Из инженерной практики известны различные методы НК [59−71, 76]:

— радиационные,.

— ультразвуковые,.

— интроскопия,.

— магнитные,.

— капиллярные,.

— тепловые,.

— рентгеновские методы,.

— вихретоковые,.

— электрические,.

— оптические,.

— вибрационные,.

— акустические, и др.

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

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

Дефекты конструкций существенно отличаются по степени опасности. Определенные формы дефектов вообще не являются концентраторами механических напряжений и не представляют опасности для работы конструкции. И наоборот, некоторые виды концентрации механических напряжений не обнаруживаются дефектоскопами, хотя очень опасны. Дефектоскопы не показывают фактических значений напряжений, а, следовательно, и степени опасности дефекта. К тому же, области остаточных напряжений и зоны с остаточной пластической деформацией многие известные методы не обнаруживают. Для значительного количества дефектоскопов требуется вспомогательная обработка поверхности в зоне проведения замера, что снижает оперативность работы, повышает ее стоимость, а так же негативно сказывается на состоянии испытуемого изделия. Таким образом, все традиционные технологии диагностики ориентированны на поиск дефектов и измерение их размеров, а об опасности найденных дефектов судят по результатам расчетов, исходя из размеров и топологии дефекта. Следовательно, достоверность оценок полностью зависит от модели и метода расчета [61]. Данные расчеты [61] основаны на инженерных методиках, использующих статистическую информацию замеров, проводимых на большом количестве образцов. К сожалению, в силу того, что измерения проводятся в конечном числе точек, результаты мониторинга не дают информации о распределениях напряжений в труднодоступных областях изделия, особенно если оно испытывает сложно распределенные нагрузки. В связи с этим, существенно возрастает роль математического моделирования как средства оценки непосредственного состояния испытуемого объекта. А также, по полученным результатам моделирования можно сделать выводы о сроках и экономической целесообразности эксплуатации в определенных эксплуатационных режимах данного изделия.

При проведении дефектоскопии объектов, выполненных из ферромагнитных материалов, для оценки механических напряжений часто используют магнитные методы. Среди известных магнитных методов на практике применяется метод «магнитной памяти металла», заключающийся в регистрации изменения нормальной к поверхности изделия, находящегося в магнитном поле, составляющей напряженности магнитного поля. [72]. Физические основы данного метода заключаются в том, что в областях концентрации напряжений и деформаций возникают локальные изменения магнитных свойств материала. Основными достоинствами данного метода являются: простота проведения замеров, отсутствие необходимости проведения предварительной подготовки испытуемого изделия, возможность проведения замеров на сравнительно большом количестве изделий, а так же контроль протяженных конструкций в сжатые календарные сроки. Данный метод позволяет производить оценку состояния контролируемых объектов с учетом материала, условий эксплуатации, особенностей элементов конструкций. Зоны концентрации деформаций при использовании этого метода предложено фиксировать на поверхности труб по изменению нормальной составляющей напряженности магнитного поля.

Зависимость между магнитным состоянием тел и механическими нагрузками описывается в теории ферромагнетизма [73, 74, 76, 37] как один из частных случаев общей параметрической зависимости магнитной проницаемости ферромагнетиков от температуры, напряженности и направления внешнего магнитного поля, внутренних механических напряжений, химического состава металлов. Результаты экспериментальных исследований [23, 72] показывают, что у ферромагнетика, находящегося в магнитном поле, под действием механических нагрузок возникает анизотропия магнитных свойств. Сложная многопараметрическая зависимость магнитных свойств во многих случаях затрудняет анализ результатов измерений. Поэтому одних только экспериментальных данных недостаточно для построения количественной теории магнитомеханических явлений. В то же время наличие такой теории является необходимым условием создания точных диагностических средств. В связи с этим актуальна разработка математической модели намагничивания и проведение численных исследований магнитомеханических явлений.

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

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

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

Диссертация состоит из введения, трех глав и заключения.

Основные результаты диссертации опубликованы в следующих работах:

1. Гасилов В. А., Карагичев А. Б., Карташева E. JL, Тарасов Д. С. Анализ топологии и декомпозиция полиэдров. Фундаментальные физико-математические проблемы и моделирование технико-технологичеких систем. Выпуск 2. М.: МГТУ «Станкин» 1999. 165−177с.

2. A.S. Boldarev, V.A. Gasilov, E.L. Kartasheva, A.B. Karagichev, D.S. Tarassov. Modeling of transport and deposition of particles in two-phase turbulent flows. Фундаментальные физико-математические проблемы и моделирование технико-технологичеких систем. Выпуск 2. М.: МГТУ «Станкин» 1999. 74−83с.

3. Гасилов В. А., Захаров C.B., Круковский А. Ю., Новикова Т. П. Тарасов Д.С. Алгоритм расчета радиационного энергообмена в импульсной магнитоускоряемой плазме. Фундаментальные физико-математические проблемы и моделирование технико-технологичеких систем. Выпуск 3. М.: МГТУ «Станкин» 2000. 73−86с.

4. Гасилов В. А., Карагичев А. Б., Карташева Е. Л. Дутанова A.B., Тарасов Д. С., Фрязинов О. В. Интерактивная система генерации тетраэдральных сеток ITERA. В сборнике трудов «Четвертая мсеждународная конференция по математическому моделированию» М.: МГТУ «Станкин» 2001. Т.2. 84−91с.

5. Гасилов В. А., Гуськов С. Ю., Карагичев А. Б., Карташева Е. Л., Кузнецов Е. В., Тарасов Д. С. Физико-математическия модель для анализа магнитомеханического эффекта. Препринт № 14 М.- Физич. Инст. им. П. Н. Лебедева РАН, 2002 16с.

Заключение

.

Сформулируем основные результаты данной работы.

Создан комплекс программ подготовки геометрических данных для решения трехмерных задач математической физики. Комплекс включает генераторы сеток, созданные на основе методов:

• объединения сечений (кинематический),.

• фронтальной тетраэдризации,.

• объединения гранично-согласованных сеток.

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

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