На правах рукописи
СЕЛИН Илья Александрович
МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
ПРИ ЧИСЛЕННО-АНАЛИТИЧЕСКИХ РЕШЕНИЯХ ЗАДАЧ
ДЛЯ УРАВНЕНИЙ ВОЛНОВОГО ТЕПЛОПЕРЕНОСА
Специальность 05.13.18 – «Математическое моделирование,
численные методы и комплексы программ»
АВТОРЕФЕРАТ
диссертации на соискание учёной степени кандидата физико-математических наук
Москва – 2010
Работа выполнена на кафедре «Вычислительная математика и программирование»
Московского Авиационного Института (Государственного Технического Университета).
Научный руководитель: заслуженный деятель науки, доктор физико-математических наук, профессор Формалёв Владимир Федорович
Официальные оппоненты: доктор физико-математических наук, профессор Волков Игорь Куприянович кандидат физико-математических наук, доцент Жаворонок Сергей Игоревич
Ведущая организация: Московская государственная академия тонкой химической технологии им. М.В.Ломоносова
Защита состоится «24» декабря 2010 года в 12 ч. 00 мин. на заседании диссертационного совета Д 212.125.04 при Московском авиационном институте (государственном техническом университете) по адресу: 125993, г.Москва, А-80, ГСП-3, Волоколамское шоссе, д. 4, Учёный совет МАИ.
С диссертацией можно ознакомиться в библиотеке Московского Авиационного Института (Государственного Технического Университета).
Автореферат разослан «_» 2010 года.
Учёный секретарь диссертационного совета Д 212.125.04, кандидат физико-математических наук М.В. Ротанина
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Объект исследования. Объектом исследования в диссертационной работе является процесс волнового переноса потенциала при высокоинтенсивных воздействиях.
Актуальность работы. При исследовании процессов переноса потенциалов (температуры, давления, концентрации и т.д.) при высокоинтенсивных воздействиях (при взаимодействии лазерного излучения с твёрдой поверхностью, входе гиперзвуковых летательных аппаратов в плотные слои атмосферы, сверхзвуковых контактных взаимодействиях и т.п.) гипотеза о пропорциональности вектора потока вектору градиента потенциала, основанная на феноменологическом подходе, приводит к парадоксу о бесконечной скорости распространения возмущений, что противоречит фундаментальным законам естествознания. Уравнения переноса возмущения потенциала, полученные на основе этих законов, являются уравнениями параболического типа.
Попытка учесть конечную скорость возмущения в рамках феноменологических гипотез не привела к успеху. Позднее, в различных областях естествознания были высказаны гипотезы о распространении возмущений как о процессе, являющемся одновременно волновым и диффузионным. Первым, кто высказал эту гипотезу был Максвелл, развитие идей которого привело к законам, где появилось слагаемое, учитывающие релаксационные явления в средах.
Математическое моделирование явлений переноса потенциала на основе этих законов приводит к уравнениям гиперболического типа, учитывающих внутреннюю вязкость среды и описывающих волновой перенос в условиях высокоинтенсивных воздействий с конечной скоростью.
Математическому моделированию процессов волнового переноса при высокоинтенсивных воздействиях посвящено незначительное число работ. Это связано со сложностью постановки и решения подобных задач, так как фронт волны математически представляет собой сильный разрыв потенциала поля (разрыв первого рода) или слабые разрывы (разрывы первых или вторых производных). Кроме этого с помощью численного моделирования сложно точно выделить фронт волны, а следовательно и решение до и после фронта.
Из работ, в которых наиболее полно рассмотрен волновой перенос можно отметить монографию Шашкова А.Г., Бубнова А.В., Яновского С.Ю. Однако в ней практически полностью отсутствует математический аппарат исследования волновых процессов переноса.
В работах Соболева С.Л. проведён анализ волнового переноса субстанции на основе релаксационных явлений в неравновесных средах и высказано предположение о существовании ударных волн в нелинейных средах. Однако математический аппарат исследования волнового переноса возмущений так же отсутствует.
В работах Самарского А.А. и его школы исследуется волновой перенос в нелинейных одномерных средах на основе градиентных законов. Для одномерного пространства приводится теорема о существовании волнового переноса в случае зависимости коэффициента переноса от потенциала. В работе Баутина С.П. предложен метод функциональных рядов для исследования волнового переноса в одномерных средах.
Вопросами волнового переноса занимались так же Лыков А.В., Корнеев А.С., Рубина Л.И., Формалев В.Ф., Баренблатт Г.И.
В перечисленных работах отсутствует общий подход к математическому моделированию волнового переноса в линейных и нелинейных средах, а так же в средах с анизотропией свойств. Кроме этого, отсутствуют аналитические и численные методы решения задач волнового переноса в многомерных линейных и нелинейных анизотропных средах, а так же математический аппарат исследования ударных волн в твёрдых средах.
Вместе с тем большинство материалов естественного и искусственного происхождения в той или иной степени являются анизотропными, для которых предполагается, что скорость распространения возмущений в различных направлениях различна. Пренебрежение волновым переносом приводит к значительным погрешностям в существенно нестационарных процессах в окрестности начального времени, которые, не затухая, распространяются и на большие времена.
Нерешённость перечисленных актуальных проблем обусловила цель данной диссертации: разработка математического аппарата и численно-аналитических методов моделирования волнового переноса в линейных и нелинейных анизотропных средах и определение новых явлений, сопровождающих волновой перенос.
Методы исследования. В работе использовались методы интегральных преобразований, методы подобия и размерностей, математического анализа, численные методы, методы математического моделирования.
Достоверность результатов диссертационной работы обеспечивается адекватными математическими моделями, точными аналитическими решениями, использованием абсолютно устойчивых численных методов и сравнением результатов с аналитическими решениями.
Научная новизна.
– предложен новый метод построения математических моделей, описывающих процессы переноса при высокоинтенсивных воздействиях, на основе представления законов переноса возмущений в виде ряда по малому параметру.
– доказана теорема об эквивалентности волнового переноса на основе гиперболического уравнения и волнового переноса на основе квазилинейного параболического уравнения.
– предложен метод решения уравнений в частных производных со смешанными производными, с помощью которого впервые получены аналитические решения задачи Коши и краевой задачи волнового переноса в нелинейном анизотропном пространстве.
– впервые показано существование ударных волн в нелинейных пространствах и сформулированы условия их возникновения.
– предложена модификация метода переменных направлений с экстраполяцией для численного решения уравнений волнового теплопереноса по увеличению порядка аппроксимации.
Практическая значимость.
Практическая значимость результатов работы состоит в следующем:
1. Разработан программный комплекс, реализующий расчёты температурных полей при высокоинтенсивных тепловых воздействиях.
2. Предложенный математический аппарат позволил моделировать волновой перенос на основе уравнения параболического типа с запаздывающим аргументом по времени.
3. Полученные аналитические решения волнового переноса в нелинейном анизотропном пространстве позволили установить не только качественные, но и количественные характеристики волн в твёрдых средах.
4. Модифицирован и обоснован метод переменных направлений с экстраполяцией для решения задач волнового теплопереноса в анизотропных средах.
Апробация работы. Полученные результаты докладывались на следующих конференциях: 16-ая международная конференция по вычислительной механике и современным прикладным программным системам (Алушта, 2009); XV Международный симпозиум «Динамические и технологические проблемы механики конструкций и сплошных сред им. А.Г.Горшкова» (Ярополец, 2009); 8-ая международная конференция по неравновесным процессам в соплах и струях (Алушта, 2010); 5-ая Российская национальная конференция по теплообмену (Москва, 2010).
Публикации. Основные результаты диссертации опубликованы в 11 работах, в том числе в 6 журналах, рекомендованных ВАК.
Структура и объём диссертации. Диссертация состоит из введения, четырёх глав, заключения и списка литературы. Объём составляет 141 страницу, 66 рисунков, 93 источника в списке литературы.
СОДЕРЖАНИЕ РАБОТЫ
Во введении обоснована актуальность работы, описаны проблемы волнового переноса, проведён обзор литературы, представлены общие положения процесса, перечислены решаемые задачи в соответствии с целью работы.
В первой главе описан существующий метод построения моделей переноса потенциала при высокоинтенсивных воздействиях. Данный метод основан на обобщённых градиентных законах переноса, которые в общем случае для изотропной среды имеют вид:
где q вектор плотности потока, коэффициент переноса, U потенциал, время релаксации. Уравнение переноса потенциала при высокоинтенсивных воздействиях, полученное на основе этого закона, в одномерном случае имеет вид Для описания процессов переноса в анизотропных средах уравнение (1) записывается в следующем виде где тензор коэффициентов переноса. В соответствии с (3), уравнение, описывающее поведение потенциала в анизотропной среде будет иметь вид Уравнения (2) и (4) относятся к уравнениям гиперболического типа и учитывают как диффузионный так и волновой перенос. При высокоинтенсивных процессах, слагаемое играет существенную роль. На основе уравнения (2) поставлены начально-краевые задачи в одномерных и многомерных областях. Получены аналитические решения для первой, второй и третьей краевых задач.
Так как время релаксации для твёрдых тел имеет порядок 10 13 10 10 c, то можно рассмотреть следующую технологию моделирования: записав (3) в виде можно считать, что слагаемые слева представляют собой первые два члена разложения в ряд Тейлора по степеням вектора плотности. Весь ряд имеет вид Тогда (5) можно записать следующим образом или, вводя новую переменную t = t +, а затем заменяя t на t, что физически означает отставание по времени вектора плотности потока от градиента потенциала на величину. На основе выражения (6) получаются следующие дифференциальные уравнения переноса потенциала:
для одномерного случая для анизотропного трёхмерного случая Уравнения (7) и (8) являются уравнениями параболического типа с запаздывающим аргументом по времени, которые значительно проще реализуются в вычислительных процедурах по сравнению с уравнениями (2) и (4).
Проводя аналогию между законами (1) и (6) и полученными на их основе уравнениями (2) и (7) ясно, что уравнения (7) и (8) можно получить не из закона переноса (6), а непосредственно из уравнений (2) и (4), применив предложенный метод моделирования.
Данный подход в моделирования возможен в случаях, когда время релаксации является малой величиной по сравнению с характерным временем процесса.
При численном решении задач для уравнения (7) оказалось, что при использовании неявной схемы с шагом интегрирования по времени равном времени релаксации, неявная схема трансформируется в явную для классических уравнений переноса.
Для получения значений функции при любом шаге интегрирования по времени предложен новый метод численного решения, в соответствии с которым временная область разбивается на блоки, равные коэффициенту релаксации, а блоки разбиваются на слои равные шагу численного интегрирования по времени. Тогда, при вычислении значений на новом временном слое, необходимые данные берутся не из предыдущего временного слоя, а из предыдущего блока на соответствующем временном слое.
Применение рассмотренных подходов моделирования явлений переноса при высокоинтенсивных воздействиях проанализировано на примере процессов теплопереноса.
При аналитическом решении краевых задач на основе уравнения (2) использованы параметры характерные для процессов теплопереноса в твёрдых телах. Проведено сравнение полученных решений начально-краевых задач для уравнения (2) с решениями краевых задач для уравнения теплопроводности (рис. 1). Во всех случаях решение уравнения (2) имеет затухающий разрыв первого рода – фронт тепловой волны. Получены аналитические выражения закона затухания волны и скорость распространения волны. В данных задачах, по истечении времени ~ 10 решения уравнения (2) и классического уравнения теплопроводности совпадают.
Сравнение аналитического решения первой начально-краевой задачи для уравнения (2) в разные моменты времени {1,3,5} c аналитическим решением первой начально-краевой задачи При решении первой и второй начально-краевых задач в двумерной области для уравнения, полученного из закона (1), показано, что разрыв первого рода существует в направлении обоих осей (рис. 2). Решения получены методом переменных направлений с экстраполяцией. На рисунке 3 показано сравнение решения, полученного в двумерной области, с аналитическим и численным решением в одномерной области.
Рис. 2. Численное решение двумерной задачи: слева – вдоль оси y,справа – вдоль оси x на Рис.3. Сравнение решения первой краевой задачи {3} для уравнения, полученного из закона (1) в двумерной области, с аналитическим {1} и численным {2} решениями первой краевой Из рисунка видно, что полученное численное решение верно.
Сравнение численных решений первой и второй краевых задач для уравнения теплопроводности с запаздывающим аргументом по времени с аналитическими решениями уравнения вида (2) показало, что интегральное расхождение между решениями составляет менее 4% (рис. 4), что подтверждает правомерность предложенного метода моделирования на основе представления законов переноса возмущений в виде ряда по малому параметру и предложенного численного метода решения.
Рис.4. Сравнение численного решения уравнения с запаздывающим аргументом с аналитическим в одинаковых точках в моменты времени 2 10 10 c (слева) и 4 10 10 c (справа).
Во второй главе доказана Теорема1 (об эквивалентности процесса переноса потенциала на основе уравнений гиперболического типа и квазилинейных уравнений параболического типа).
Пусть (U ) = U и > 0, функция U ( x, t ) непрерывна и существуют первая производная по времени и вторая производная по пространственной переменной. Тогда, квазилинейный параболический и гиперболический дифференциальные операторы эквивалентны в смысле волнового переноса.
Эта теорема позволила исследовать волновой перенос потенциала на основе квазилинейных уравнений.
Предложен новый метод аналитического решения квазилинейных дифференциальных уравнений переноса потенциала в анизотропной среде, который заключается во введении последовательной цепочки автомодельных переменных. Таким образом исходное уравнение значительно упрощается с введением каждой новой переменной.
Этим методом получено решение задачи Коши о распространении возмущений в анизотропной нелинейной среде от мгновенного источника:
и – главные компоненты тензора переноса.
Цепочка применяемых преобразований выглядит следующим образом:
Аффинное преобразование x = cos( ) sin( ) Замена переменных q1 = r cos( ), q2 = r sin( ).
Решение имеет вид С помощью этого же метода получено решение следующей первой начально-краевой задачи:
Решение имеет вид где a =, b аналог коэффициента L в задаче Коши.
В обоих случаях установлено, что областями ненулевого решения в различные моменты времени являются области, ограниченные эллипсами, т.е. фронт возмущений имеет форму эллипсов. При увеличении параметра поверхности, образуемые графиками функции U ( x, y, t ), имеют большую область определения (рис. 4).
При ;1 (1; ) производная по времени и первые производные по пространственным переменным существуют только внутри областей, ограниченных эллипсами. Первые производные по времени и по пространству не существуют в точках фронта при > 1, а вторые производные по пространственным переменным не существуют в точках фронта при > пространственным переменным существуют всюду, но в точках фронта нет производных более высокого порядка.
Применяя данные модели к теории теплопроводности, показано, что при степенной зависимости коэффициента теплопроводности от температуры, где показатель степени больше нуля, теплоперенос имеет волновой характер. Из полученных выражений можно определить скорости распространения тепловых волн в разных направлениях.
Правильность решений, полученных с помощью предложенного метода подтверждается сравнением с одномерными решениями А.А. Самарского.
Рис. 5. Волновой теплоперенос в нелинейной анизотропной среде.
Сформулировано и доказано условие существования решения в виде бегущей волны для нелинейного анизотропного двумерного пространства:
Теорема2. Пусть функция U ( x, y, t ) удовлетворяет уравнению Тогда при выполнении неравенства Существует решение, ограниченное по пространственным переменным, т.е. существует фронт волны.
Для одномерного случая это условие имеет вид При совместном использовании двух методов моделирования волнового переноса – на основе релаксационных явлений и на основе использования зависимости коэффициента переноса от потенциала получается следующее уравнение где (U ) = U, > 0, с помощью которого можно моделировать волновой перенос с возрастающей амплитудой. Так, численное решение данного уравнения в рамках теории теплопроводности показало существование тепловой ударной волны в твердом теле при импульсном источнике на границе. Формирование тепловой ударной волны следует из того, что каждый следующий импульс движется по более прогретой области относительно предыдущего импульса. Поскольку (T ) = T и > 0, то с ростом температуры увеличивается скорость каждого следующего импульса. Таким образом, фронт тепловых волн накладывается, образуя тепловую ударную волну (рис.5).
Проведено исследование влияния параметров импульсного источника на формирование тепловых ударных волн. Тепловая ударная волна формируется всегда при выполнении условий: 1 и импульсный источник с характеристикой 2, где d 2 – время существования импульса на границе, а d 1 – время его отсутствия, т.е. отношение амплитуд текущего импульса к предыдущему больше или равно двум.
Показано, что решение, полученное без учёта времени релаксации даёт значительные погрешности не только в начальные моменты времени, а на протяжении всего процесса, и не позволяет получить тепловые ударные волны.
В третьей главе для решения задач волнового теплопереноса на основе закона Максвелла-Каттанео в анизотропной среде используются численные методы. Впервые для уравнений гиперболического типа применён метод переменных направлений с линейной экстраполяцией по времени. Шаблон схемы представлен на рисунке.
Доказана устойчивость данной схемы и определён порядок аппроксимации.
Первая краевая задача решалась в следующей постановке:
где 11, 12, 21, 22 – компоненты тензора теплопроводности.
При увеличении по отношению к увеличивается скорость тепловой волны в направлении оси x при = 0 ; при этом разрыв на фронте тепловой волны уменьшается, что вызвано более интенсивным распространением тепла вдоль оси x и большей диссипацией в этом направлении. При увеличении по отношению к аналогичные результаты проявляются в направлении оси y. При изменении угла происходит смещение тепловой волны: при ;0 смещение изотерм и разрывов первого рода будет происходить в левую сторону, а при 0; – в правую.
Если >, то при увеличении угла ориентации главных осей тензора теплопроводности происходит более интенсивный прогрев. При ; коэффициент при второй производной станет отрицательным и изотермы направляются в другую сторону с отклонением на угол меньший на по отношению к заданному.
Вторая краевая задача решалась в следующей постановке:
При решении этой задачи, для аппроксимации краевого условия со вторым порядком применялась модификация данного численного метода с использованием интегроинтерполяционного метода Самарского А.А., что позволило увеличить порядок аппроксимации краевых условий, содержащих производные.
При увеличении по отношению к увеличивается скорость тепловой волны в направлении оси y. В отличие от первого краевого условия, при увеличении температура на границе снижается. Когда >, то тепло растекается вдоль оси x и тепловая волна имеет более широкий фронт, при этом уровень температуры на границе так же снижается. При изменении угла ориентации главных осей тензора теплопроводности, установлено, что когда >, то максимальные значения температур будут наблюдаться = 30 0.
Из рассмотренных задач в анизотропной области так же установлено, что угол ориентации главных осей тензора теплопроводности O, O влияет одинаково на решение при моделировании теплопереноса с учётом времени релаксации и без учёта.
«