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

Г I □ ин

Институт Математического Моделирования Российской Академии Наук

Центральный физико-технический институт Министерства Обороны России

На правах рукописи УДК 517.957+519:62

Шапранов Александр Викторович

МЕТОД ДИНАМИЧЕСКОЙ АДАПТАЦИИ В НЕСТАЦИОНАРНЫХ КРАЕВЫХ ЗАДАЧАХ

05.13.18. — Теоретические основы математического моделирования, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой, степени кандидата физико-математических наук

Москва 1994

Работа выполнялась в Институте математического моделирования Российской Академии Наук и в Центральном физико-техническом институте Министерства Обороны России.

Научный руководитель:

кандидат физ.- мат. наук, старший научный сотрудник В.И. Махукин.

Официальные оппоненты:

доктор физ,- мат. наук, старший научный сотрудник В.А. Гасилов,

кандидат физ.- мат. наук, доцент С.А. Волошин.

Ведущая организация:

Институт математики АН Беларуси.

Защита диссертации состоится ' заседании специализированного

_ 1994 года

003.91.01 при

на заседании специализированного совета К Институте математического моделирования РАН по адресу: 125047, Москва, Миусская площадь, дом 4.

С диссертацией мохно ознакомиться в библиотеке ИММ РАН Автореферат разослан "2^" г.

Ученый секретарь специализированного совета К 003.91.01 при ИММ РАН

к. ф.- м. н.

С.Р. Свирщевский

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

Цель работы. Целью данной диссертации является дальнейшее развитие метода динамической адаптации применительно к нестационарным краевым задачам математической физики.

Научная новизна Основным элементом научной новизны является способ определения оптимальной функции преобразования к нестационарной системе координат в методе динамической адаптации. С этой целью вводится требование квазистационарности процессов в новой системе координат. Разработаны и проанализированы разностные схемы для уравнения Бюргерса и модельной задачи об аннигилирующих импульсах.'

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

Апробация диссертации. Материалы диссертации докладывались на международной конференции Европейского сообщества исследования материалов (Е-МИБ, Франция, Страсбург, 1993); на международном лазерном конгрессе (Франция, Лион, 1993); на научном

семинаре профессора Леванова Б.И.

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

Структура и объем работы. Диссертация состоит из введения, четырех глав, заключения и списка литературы, содержащего 131 наименование. Работа изложена на 115 страницах машинописного текста, содержит 32 рисунка.

СОДЕРЖАНИЕ РАБОТЫ.

Во введении показана актуальность темы диссертации, сформупиро-. вала цель работы и приведено ее краткое содержание по главам.

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

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

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

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

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

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

З^-'Н С)

В данном уравнении Р(и) представляет собой поток величины и вдоль оси X и может быть нелинейной дифференцируемой функцией от и. Правая часть Р [и] в общем случае является нелинейным дифференциальным оператором произвольного порядка. В данной работе рассмотрение ограничивается операторами не выше второго порядка.

Запишем уравнение (1) в произвольной нестационарной системе координат (д, т):

ди Оди . 1£Р(в2_«ы (2)

дх _ _ дх ■ = * = ^

или в консервативном виде:

+ = (з)

Здесь Ф - метрический коэффициент, дх/дт - скорость движения системы координат.

Полученное уравнение (3) отражает динамику исходной задачи (1) в новых переменных (<?, т). Уравнение же (4) является уравнением обратного преобразования переменных (д,г) в исходную систему координат (х, <). Задание функции <2 определяет конкретный вид преобразования координат и является параметром управления движением узлов сетки.

Ранее функция (} задавалась всходя из следующих соображений.

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

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

Ф

» = 0,1,2,...

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

„5Ф дР(и) Э ( |3и|\ д /т|92«

Я = -о^г- - х-—1 - *Ы - а~я1

дд ди дд V. \dq\j дд \ \дq1

(5)

где £>, о - подгоночные коэффициенты. Наличие подгоночных коэффициентов в выражении (5) для функции С}, численные значения которых определяются эмпирическим путем, несколько усложняют ее использование в произвольном случае.

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

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

разностных схем, и их уменьшение позитивно сказывается на этих свойствах.

Исходя из этого, при выборе управляющей функции 0 потребуем, чтобы была выбрана такая нестационарная система координат, в которой временные производные решения были бы равны нулю или достаточно малы. Это требование в дальнейшем будем называть требованием ква-зистационарностЕ. Предельным случаем здесь является равенство нулю временной производной решения во всех точках области определения в системе координат (д, г). Подставляя ди/дт = 0 в уравнение (2), можно найти функцию С}, которая должна обеспечить выполнение этого равенства:

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

Можно показать, что функция в виде (6) обеспечивает выполнение еще одного равенства:

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

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

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

где /( - малый параметр и имеет смысл вязкости.

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

(6)

(7)

_ / д ¡ti\ ¡id2u /ди\

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

1. Оптимальную функцию преобразования Q можно определить из условия стационарности процессов в новой нестационарной системе координат.

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

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

Выводы теоретического анализа подтверждаются вычислительными экспериментами. На рис.1,2 представлено решение задачи Бюргерса с начальным условием,' задаваемым в виде асимметричной синусоиды:

и0 (х) = и (х, 0) = sin (2пх) + 0.5 sin (пх)

Значение параметра ц принималось равным Ю-4.

Вначале задача решалась на сетке с фиксированными узлами с использованием разностной схемы Кранка - Николсона. Расчеты показали, что на сетках с числом узлов менее 100 численное решение фактически неустойчиво. На рис.1 представлено решение полученное на сетке с 103 узлами. Увеличение числа узлов до 104 позволяет уменьшить амплитуду наблюдаемых паразитных осцилляции, но не избавиться от них полностью. Данный эффект свидетельствует о чрезвычайно сильной дисперсии используемой разностной схемы.

На рис.2 представлено решение, полученное на адаптивной сетке с общим числом узлов N — 21, на те же моменты времени, что и на рис.1.

8

М=1(Г4 го=0.0

г, =0.2

1г=0.Ч 1Э = 0.6 и = 0.8

Рис.1. Расчет асимметричной синусоиды на Эйлеровой сетке (1001 узел)

-0.5

-1.0

Рис.2. Расчет асимметричной синусоиды на адаптирующейся сетке (21 узел)

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

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

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

'%Г+«17П = -*/>1Ра, 0<х<Х, t> О а Э <9)

Здесь - скорости распространения импульсов, г.>1 > 0, «2 < 0> & -коэффициент взаимодействия.

Начальные условия:

Р? = Р1

О < X < XI XI < х < X

(10)

х2 < х < X

0 < X < Х2

где: 0 <х1 < х2 < X

Граничные условия : /91(0, £) = 0, рг(Х, ¿) = 0, < > 0. (11)

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

Qn = -(vn + Я>„kplp2(^J п = 1,2

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

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

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

В четвертой главе с помощью метода динамической адаптации моделируется процесс быстрых фазовых переходов при импульсной лазерной обработке сверхпроводящей керамики. Математическое описание процессов лазерного нагрева, плавления и испарения сверхпроводящей керамики производится в рамках совмещенного варианта задачи Стефана, объединяющего классическую и однофазную постановки. Математическая модель совмещенного варианта представляет собой краевую задачу для уравнения теплопроводности с двумя подвижными фазовыми границами Г„| (1) и Г|„ (<). Учет объемного нагрева твердой и жидкой фаз производится с помощью члена дй/дх в уравнении теплопроводности и уравнения переноса излучения:

(12)

+ кв = 0 ^ к = в,1

х0<х< ТВ1 (*)

\

Г», (<) < X < Г,„ (<), (>0

Условия на межфазных границах:

эг, ад

£ = I»/ (I) : ~ = Р'ьту>1'

= Г, = 7} = Тт; х = Г,„(<): = — аТ*г\

Р^ь = р„ (иь - и); Р«г + -Рь + Ръ {1>ь - и)2;

Основная сложность численного решения задач типа Стефана связана с наличием подвижных межфазных границ, положение которых Заранее неизвестно и должно определяться в ходе решения. В связи с тем, что область решения заранее неизвестна и меняется с течением времени, возникают трудности с построением пространственных расчетных сеток. Преодолеть эти трудности можно с помощью метода динамической адаптации. В качестве закона преобразования координат задавалась функция вида С} = —(дФ/дq), которая при достаточной величине коэффициента О обеспечивает равномерное распределение узлов на каждый момент времени.

Моделирование проводилось для двух механизмов диссипации энергии лазерного излучения: для поверхностного и объемного. В случае поверхностного нагрева пространственное распределение температуры Т(х), как в твердой, так и в жидкой фазах, рис.3, мало чем отличается от традиционного распределения в материалах с низкой теплопроводностью при воздействии поверхностного теплового источника. К моменту плавления прогретый слой не превышает 0.1 мкм, что приводит к возникновению больших температурных градиентов в приповерхностных слоях. Наличие больших градиентов температуры обуславливает высокую начальную скорость плавления и,/, которая за короткое время достигает своего максимального значения V™* ~ 6 м/с, а затем под влиянием роста жидкой фазы начинает быстро уменьшаться. Однако, дальнейший

Х(»10"4см)

Рис.3. Пространственный профиль температуры при поверхностном механизме выделения энергии

2500 2000 1500

н ЮОО-500:

1 I

2.0 1.5 1.0 0.5 0 X (*10~4см)

Рис.4. Пространственные профили температуры и интенсивности излучения при объемном выделении энергии

спад скорости плавления замедляется процессом интенсивного поверхностного испарения, максимальная скорость которого достигает величины ь%ах ~ 1.5 м/с.

Увеличение длины свободного пробега кванта до

зоооЛ приводит к

заметной диссипации энергии излучения в объеме материала. Пространственный профиль температуры в этом случае имеет плоскую вершину вблизи поверхности с малым температурным градиентом. При достижении на облучаемой поверхности температуры плавления Тт от поверхности вглубь материала начинает распространяться фронт плавления. Из -за невысоких пространственных градиентов температуры начальная скорость плавления оказывается на несколько порядков меньше начальной скорости плавления при поверхностном источнике. Затем под совмест-'ным воздействием процессов объемного нагрева и плавления в приповерхностных слоях материала начинает формироваться температурный максимум, рнс.4. Наличие приповерхностного максимума температуры означает, что некоторый объем твердой фалы оказывается перегретым относительно равновесной температуры плавления Т,п. Возникающие ме-тастабильные состояния наиболее значительно сказываются на скорости движения фронта плавления. Появление максимума температуры на некотором расстоянии от границы Г.,/ вызывает изменение направления теплового потока \Уц , идущего из твердой фазы. В результате тепловые потоки из твердой и жидкой фаз 1УЙ и И7/ на границе раздела складываются, а не вычитаются, как в традиционных ситуациях. Сложение потоков Ил5 и И'/ приводит к резкому увеличению скорости плавления, максимальное значение которой ~ 14 м/с, превосходит более чем

в два раза соответствующее значение в случае поверхностного ис-

точника. Высокие значения г>„| приводят к быстрому росту жидкой фазы, в которой поглощается все большая доля лазерного излучения. Перегрев твердой фазы при этом уменьшается вплоть до полного исчезновения.

Аналогичное перегретое состояние возникает также и в жидкой фазе.

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

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

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

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

3. С помощью метода динамической адаптации исследовано явление возникновения перегретых метастабильяых состояний, сопровождающих быстрые фазовые переходы при импульсной лазерной обработке сверхпроводящей керамики. Моделирование объемного воздействия с интенсивностью 107 Вт/см2 и длительностью 40 не позволило зафиксировать значения перегрева твердой и жидкой фаз, которые соответственно составляют более ста и нескольких сотен градусов.

Список литературы

1. В.И.Мажукян, A.A.Самарский, А.В.Шалранов Метод динамической адаптации в проблеме Бюргерса. //ДАН, 1993, т.ЗЗЗ, №2, стр.165 -169.

2. В.И.Мажукин, А.А.Самарский, Орландо Кастельянос, А.В.Шапра-нов. Метод динамической адаптации для нестационарных задач с большими градиентами. //Математическое моделирование, 1993, т.5, №4, стр.32 - 56.

3. В.И.Мажукин, И.В.Гусев, А.В.Шалранов. Влияние метастабильных состояний на процесс импульсной лазерной обработки сверхпроводящей керамики. //Математическое моделирование, т.5, N 5, 1993, > стр.30 - 60.