CC BY
149
УДК 539.37+517.91
АНАЛИЗ ЭФФЕКТИВНОСТИ МЕТОДОВ АДАМСА И ГИРА ПРИ РЕШЕНИИ ЖЕСТКИХ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ В ПАКЕТЕ SPFCC
М.Е. Семенов1,2, С.Н. Колупаева1
Томский государственный архитектурно-строительный университет 2Томский политехнический университет E-mail: sme@tpu.ru; ksn@tsuab.ru
Проведен анализ эффективности использования численных методов Адамса и Гира при решении жесткой системы обыкновенных дифференциальных уравнений в пакете SPFCC. Описан подход, позволяющий обнаруживать участки интервала интегрирования, где численная устойчивость метода накладывает ограничения на длину внутреннего шага интегрирования.
Ключевые слова:
Жесткая система обыкновенных дифференциальных уравнений, численные методы Адамса и Гира, область устойчивости. Key words:
Stiff set of differential equations, numerical Adams and Gear methods, stability region.
В последние годы достижения в области математического моделирования и вычислительного эксперимента как информационной технологии получения новых знаний об окружающем нас мире приобретают все большее значение для различных областей наук. При математическом описании различных явлений природы и многих фундаментальных законов окружающего нас мира используется аппарат дифференциальных уравнений [1-3]. Исследование эволюции деформационной подсистемы, роли различных механизмов генерации и аннигиляции дефектов различного типа в деформационном упрочнении для металлов с гранецентрирован-ной кубической структурой (ГЦКС) и сплавов на их основе в условиях различных воздействий является лишь одним примером задач, при решении которых может эффективно использоваться математическое моделирование с использованием дифференциальных уравнений, в том числе балансового типа.
Однако получение точных решений дифференциальных уравнений возможно лишь для ограниченного класса задач. В связи с этим возникает необходимость в разработке и использовании различных методов, позволяющих получить приближенное решение задачи. Выбор, сборка и разработка численных методов для решения систем обыкновенных дифференциальных уравнений (ОДУ) баланса деформационных дефектов (которые, как правило, являются жесткими) при описании процессов, определяющих пластическое поведение материалов в условиях различных деформирующих воздействий, является одним из необходимых этапов исследований [4, 5]. Среди различных численных методов решения дифференциальных уравнений важную роль играют разностные методы вычислений. Их существенным достоинством является простая алгоритмизация и компьютерная реализация [1].
Математическая модель процессов
пластической деформации
Для исследования закономерностей пластической деформации и эволюции деформационной
дефектной среды в условиях деформации с постоянной скоростью деформирования, при постоянном напряжении либо постоянной нагрузке (в условиях растяжения или сжатия) для металлов и дисперсно-упрочненных материалов с ГЦКС с недеформируемыми частицами второй фазы создан пакет прикладных программ Slip Plasticity of Face-Centered Cubic (SPFCC) [6-9]. Для описания процессов пластической деформации в пакете SPFCC использованы математические модели, включающие уравнения баланса деформационных дефектов различного типа [7, 10, 11]. В настоящей работе анализ эффективности используемых численных алгоритмов приведен на примере математической модели пластической деформации для дисперсно-упрочненных материалов с недеформируемыми частицами [4, 7, 10]:
^ = (1 )D- - -(1 -Ч)pIbmm(r , p-J'i)x
da Db a
x(c2 Q + cQ + cQ) +
i^4p(pp (C1vQ1v + C2vQ2v) + p'pCiQi) + ^
2b ¿1
+ ~(PdCQ + Pd (C1vQ1v + C2vQ2v ))
r
(1)
dpp _<X>$
da
2Л 2pb
-—4Pp'pb(2c2vQ2v + 2CQ) - ^bp cQ , (2)
a a
d PP _<X>S 4a b JppvcQ^ -da 2ЛрЬ a p ' '
2a
- '2ab4ppVp (2c2vQ2v + c1vQ1v )>
a
(3)
d p
1 2b v „ „ 2b v . ^ „ „ . ,
d ^r - -a-pdcQ -—pd(c2vQ + cQX(4)
da Ль ara ar
d p'd da
1 2b ■ 2b ■ ТГ -—pd (c2vQ2v + cQ )- — pdciQ , (5)
Л b air ar.
Сс = дтфп -Са О Г ((1 -Щ )рт + Рр + Рс ) х *ь2д1 + о,ъсъ + +д (си + е2 у )
(6)
Сс,
IV
Са
6О
(
(((1 )Рт +Рр + Рс )Ь2 + С + ^ ) Х
Х Q1vС1v + Си 1
т(а
^ + д )С С2v ,
(7)
Сс~
5дт1
Суп
, (
Са 6О (((1 )Рт + Рр + Ра)Ь2 + С )х
2v + Й4 С2v
Л
хд..с.
+-е^2- (8)
Переменные модели: р - суммарная плотность дислокаций; рт - плотность сдвигообразующих дислокаций; р/, рр' - плотность дислокаций в призматических петлях вакансионного, межузельного типа; р/ и р1! - плотность дислокаций в дипольных конфигурациях вакансионного и межузельного типа; с, с2„ и с„ - концентрация межузельных атомов, бивакансий и вакансий; а - степень деформации сдвига. Параметры модели [12]: 5 - диаметр упрочняющих частиц; Б - диаметр зоны сдвига; Ь - модуль вектора Бюргерса, параметр ^ определяется формой дислокационных петель и их распределением в зоне сдвига; тйуп - напряжение, избыточное над статическим сопротивлением движению дислокаций; ю - доля винтовых дислокаций; Qj - кинетический коэффициент реакции для ¡¡-го типа дефекта (¡=1,V); <х> - средняя величина параметра, характеризующего геометрию дислокаций на частицах; q - параметр, определяющий интенсивность генерации точечных дефектов; а - параметр междислокационных взаимодействий; О - модуль сдвига; х1 - напряжение трения; та - атермическая составляющая напряжения; Ра! - вероятность аннигиляции винтовых дислокаций; га - критический радиус захвата; Лр - расстояние между частицами второй фазы; с=с+с(>) - суммарная концентрация точечных дефектов, с/0' - концентрация термодинамически равновесных точечных дефектов ¡¡-го типа (¡=1,у).
Система уравнений модели (1)-(8) дополнена уравнением, связывающим скорость деформации с приложенным воздействием и характеристиками дефектной среды [13]:
а = ВvDвД,2 х а ^ (1 -Д )Х
^3[((1 -Д )Рт +Рр +РС )(Т-Тд )]1/3
х О4/3Ь1/3(т2 -О2Ь2Дрт)рт1/2 Х
<ехр
0,2 ОЬ3 - (т-Т„)ЛЬ2 " кТ
(9)
где ув - частота Дебая; £ - доля дислокаций леса; г - приложенное напряжение; Д - доля реагирующих дислокаций леса; к - постоянная Больцмана; Т - температура; Л - длина свободного дислокационного сегмента.
В настоящем исследовании анализ проведен для деформации с постоянной скоростью деформирования.
Сложность численного моделирования закономерностей пластической деформации заключается в том, что переменные математической модели (1)-(9) являются разнопорядковыми величинами и на интервале интегрирования изменяются на порядки величины, а интенсивности накопления деформационных дефектов существенно различны. Кроме этого в процессе расчетов при достижении некоторых физических условий происходит трансформация системы уравнений - добавляются в отдельные уравнения дополнительные слагаемые либо к системе уравнений подключаются дополнительные уравнения. Например, уравнения (2), (3) подключаются к системе ОДУ только при условии достижения некоторого физического условия.
Для решения систем ОДУ могут использоваться различные приближенные методы численного интегрирования [2, 14, 15], однако их эффективность для решения определенного класса задач может сильно отличаться. Эффективность конкретного метода применительно к решению конкретной задачи зависит от достижения необходимой точности численного решения и приемлемого времени расчета [2, 14]. Эти два фактора тесно взаимосвязаны, так как увеличение точности влечет за собой увеличение времени расчета. Исходя из этого, определяются методы численного интегрирования, имеющие оптимальное сочетание указанных факторов при решении задачи.
Точность численного решения и расчетное время во многом зависят от выбора шага интегрирования, порядка локальной ошибки и численной устойчивости метода. Под численной устойчивостью метода понимается свойство невозрастания локальной ошибки при переходе к следующему шагу интегрирования [2, 14]. Оптимальная стратегия использования методов подразумевает наличие процедуры автоматического управления порядком метода интегрирования и величиной шага [2, 15].
К настоящему времени не существует универсальных методов решения жестких нелинейных систем дифференциальных уравнений [14]. Надежных простых способов оценки степени жесткости системы ОДУ не существует, и поэтому необходимы численные методы, работающие без проверок на жесткость. Требованиям, предъявляемым к методам решения жестких систем уравнений, удовлетворяют, в частности, неявные алгоритмы Адам-са-Мултона и неявные многошаговые алгоритмы Гира различного порядка [2, 15].
Достаточно высокую устойчивость имеют численные методы прогноза-коррекции. При этом прогнозируемое значение на следующем шаге ин-
с
тегрирования сначала определяется по явному од-ношаговому или многошаговому методу, а уточненное (корректирующее) значение - по неявному методу интегрирования. Например, прогноз осуществляется по методу Адамса-Бешфорта, а коррекция - по методу Адамса-Мултона [14].
В настоящее время существует ряд мощных математических пакетов программ широкого назначения (Maple, MatLab, MathCad и др.), позволяющих решать системы ОДУ различной степени жесткости. Ряд задач связанных с моделированием пластической деформации удается решить и с использованием вышеназванных пакетов программ [16, 17]. Но опыт такого подхода продемонстрировал высокую трудоемкость получения, накопления и обработки данных, особенно при комплексных исследованиях, и, как правило, необходимость некоторого упрощения моделей. Для проведения исследования процессов пластической деформации с использованием математических пакетов пользователь должен иметь представление о методах решения ОДУ, навыки работы с программой и, как правило, программирования на внутреннем языке пакета (например, для обработки физических ограничений). Исходя из этого, использование пакетов программ широкого назначения для неподготовленного пользователя представляется затруднительным или даже невозможным.
Не менее важной проблемой является задача эффективной организации хранения и дальнейшей комплексной обработки полученных результатов проведенных компьютерных экспериментов.
Поскольку система математических моделей пластического поведения материалов развивается на единой концептуальной основе, создаются модификации моделей с различными возможностями и приложениями, разрабатываются модели для материалов различного типа и для описания их поведения при различных воздействиях возникла необходимость создания специализированного программного обеспечения для поддержки моделирования. Одной из основных проблем была проблема разработки вычислительного модуля и, соответственно, выбора или сборки численного метода, позволяющего решать системы различной и изменяющейся на порядки на интервале интегрирования жесткости. Анализ показал, что для численного решения приведенной выше системы уравнений целесообразно использование неявных методов численного интегрирования, позволяющих в широком диапазоне изменять шаг интегрирования, исходя из условий достижения необходимой точности.
Анализ возможностей численного метода
Существенное изменение степени жесткости и, соответственно, поведения системы (1)-(9) на интервале интегрирования потребовало реализации в пакете SPFCC численных методов. Для анализа возможностей численных методов было проведено тестирование вычислительного модуля на класси-
ческих задачах различной степени жесткости [18]. Следующим необходимым шагом был детальный анализ поведения численных методов при решении системы модели (1)-(9): исследование изменения шага и порядка методов в процессе расчетов, изменения числа жесткости, количества вычислений правой части системы уравнений, якобиана системы, количества шагов интегрирования, количества вызовов программы решения системы линейных алгебраических уравнений.
При нахождении численного решения жесткой системы ОДУ необходимо, прежде всего, проводить интегрирование с таким шагом к, который будет обеспечивать устойчивое вычисление для самой быстрозатухающей компоненты решения.
Так, изменение величины шага может вызвать значительное увеличение времени вычислений, поскольку предварительно сохраненные «предшествующие» значения, соответствующие величине шага к, должны быть интерполированы для получения ряда преобразованных предшествующих значений, соответствующих новой величине шага к„0=тах{гг/г-1/г+1}к„, (10)
где ^ - отношение величины шага предлагаемого на следующем шаге кнов и текущей величины шага кп (коэффициент пропорциональности); /=#, д-1, д+1, индекс q означает, что порядок метода на следующем шаге не изменятся, д±1 - порядок метода на следующем шаге увеличивается/уменьшается на единицу; кп - размер шага в точке хп. Заметим, что при выборе максимального значения из гр г?-1,
в (10) необходимо запомнить значение q, которое определяет порядок метода на следующем шаге. Таким образом, размер шага кнов есть некоторая функция от порядка q используемого метода интегрирования.
В условиях, когда быстро быстроизменяющиеся переменные не оказывают значительного влияния на поведение решения, автоматический выбор шага может привести к существенному увеличению шага. Необходимо ограничивать увеличение шага к до величины шага кпол, указанного исследователем, что позволяет сохранить необходимое качество визуального представления полученного решения.
При решении системы (1)-(9) в пакете 8РБСС реализованы жестко устойчивые методы дифференцирования назад (в записи Гира). Для эффективного использования алгоритма управления размером шага интегрирования в программе 8РЁСС используется представление решения в виде вектора Нордсика [2, 14], что позволяет исключить необходимость нахождения решения с новым шагом кнов, достаточно провести интерполяцию полученного решения в точке кп в точку кпол.
Комплекс программ 8РБСС позволяет в диалоговом режиме сформировать математическую модель на основе системы уравнений (1)-(9) и получить ее численное решение при большом диапазоне изменения начальных условий и параметров задачи. При этом в процессе решения от исследователя не требуется проведения дополнительных на-
строек. В качестве иллюстративного примера рассмотрим результаты моделирования пластической деформации для дисперсно-упрочненных материалов с медной, никелевой и алюминиевой матрицей. Выбор материалов для исследования был обусловлен тем, что они достаточно полно представлены теоретическими и экспериментальными результатами других авторов. На рис. 1 приведены в качестве примера результаты вычислительных экспериментов с использованием комплекса программ 8РБСС, при следующих значениях параметров модели [12]: 6=2,5.1010 м-2, /=5, ул=1013е-1, а=0,5, а =0,3, 0=0,14, |=0,5, т/ = 1 МПа, а^-0,33, &>5=0,3, <^>=4 и начальных условиях: рт=1012м-2, Р/=Р/=Р/=Р</=0, с=сгс2=0.
Расчеты проводились до степени деформации сдвига равной единице, шаг пользователя Нпол=0,01, относительная точность решения 10-5. На рис. 1 приведены кривые деформационного упрочнения и плотностей дислокационных призматических петель и диполей, концентрации вакансий и межузельных атомов для температуры 273 К и скорости деформации а=10-4 с-1.
При степени деформации 0,1 было достигнуто значение физического критерия, при котором к системе (1)-(9) подключаются дополнительно два уравнения (2), (3). При этом на кривых, описывающих поведение дефектной подсистемы, наблюдается резкое изменение плотности дислокаций, концентрации вакансий и межузельных атомов. При этом происходит изменения числа жесткости 5 системы ОДУ баланса деформационных дефектов (рис. 2, а). На начальном интервале интегрирования (степень деформации до 0,05) число жесткости 5 системы ОДУ баланса деформационных дефектов экспоненциально возрастает с 106до 109. Это свидетельствует об увеличении степени жесткости системы ОДУ, при этом использование жестко устойчивого неявного метода Гира переменного порядка позволило на этом участке увеличить шаг интегрирования до 0,0001 (рис. 2, б). Далее вплоть до подключения дополнительных двух дифференциальных уравнений (2), (3) наблюдается незначительное изменение поведения кривой течения (рис. 1, а), и расчеты проводятся практически с максимально возможным размером шага (рис. 2, б).
При подключении новых уравнений к системе ОДУ баланса деформационных дефектов происходит переключение на явный метод интегрирования (метод Адамса переменного порядка), и шаг интегрирования уменьшается до минимально возможного (рис. 2, б). После построения с помощью метода Адамса (первого и второго порядков) необходимых точек разгона осуществляется переключение на использование неявного метода Гира. При степени деформации 0,2 число жесткости 5 достигает максимального значения, равного 1011, выходит на стационарное поведение, при этом метод Гира позволяет увеличить шаг интегрирования до максимально возможного шага пользователя Нпол=0,01.
0,0 0,2 а
Рис. 1. Зависимость: а) напряжения течения; б) плотности дислокационных призматических петель рр и плотности дислокационных диполей р; в) концентрации вакансий о и межузельных атомов с, от степени деформации
Анализ, проведенный в [19], показывает, что, всякий раз, когда жесткость системы ОДУ увеличивается, произведение длины шага Н и доминирующего собственного значения матрицы Якоби шах{|Яе(Яг)||, /=1,2,...,я, не выходит заграницы области устойчивости.
а
Рис. 2. Характеристики вычислительного процесса: а) число жесткости б; б) внутренний шаг интегрирования Ь; в) произведение длиныI шага и доминирующего собственного значения матрицыI Якоби 1г<р(Х)
На рис. 2, в приведена зависимость к<(А)= =к-шах{|Ке(А;)|}, г=1,2,...,8 для системы ОДУ баланса деформационных дефектов, где к - шаг интегрирования, шах{|Ке(А;)|| - максимальное значение абсолютной величины вещественной части среди собственных чисел якобиана системы (для всех переменных системы). С помощью этой величины можно оценить область устойчивости используемого численного метода. На начальном ин-
СПИСОК ЛИТЕРАТУРЫ
1. Абрамов В.В. Математическое моделирование движения астероида 99942 АрорЫз на основе методов Адамса с переменным шагом // Вестник Самарского гос. техн. ун-та. Сер. Физ.-мат. науки. - 2008. - Т. 1. - № 16. - С. 144-148.
2. Катаева Л.Ю., Карпухин В.Б. Численное моделирование динамических систем, описываемых обыкновенными дифференциальными уравнениями // Наука и техника транспорта. -2008. - № 2. - С. 57-66.
тервале интегрирования величина к-шах{|Яе(^)|} возрастает (используем жесткоустойчивый метод Гира), но в процессе расчетов может быть достигнуто критическое значение физического условия, определяющего начало образования еще двух типов деформационных дефектов. На рис. 2 критическое значение достигнуто при степени деформации, равной 0,1, при этом происходит перестройка системы ОДУ баланса деформационных дефектов (к системе уравнений подключаются два уравнения). При этом происходит переключение на метод Адамса и величина к<(Х) резко убывает, затем (благодаря переключению на использование метода Гира) она снова возрастает, выходит на стационарное поведение.
После прохождения степени деформации 0,1 порядок метод Гира начинает увеличиваться и достигает максимально возможного порядка (шестого), заложенного в программе 8РБСС. В работе [15] показано, что увеличение порядка метода Гира выше шестого будет нарушать требования жесткой устойчивости метода. Для сохранения устойчивости метода в программе 8РБСС предусмотрена процедура автоматического управления шагом интегрирования. На рис. 2, б, видно, что большая часть интервала интегрирования пройдена с максимально доступным шагом к. Заметим, что несколько раз потребовалось уменьшить шаг интегрирования. Это происходило на фоне увеличения числа жесткости 5. Изменение длины внутреннего шага интегрирования к оказало влияние на поведение кривой, представленной на рис. 2, в.
Выводы
Показано, что использование численных методов Адамса и Гира переменного порядка при решении жесткой системы обыкновенных дифференциальных уравнений математической модели пластической деформации с гранецентрированной кубической структурой материалов является оптимальным (устойчивость, величина шага интегрирования, объем вычислений). Для этих целей может быть использован вычислительный модуль комплекса программ 8РБСС.
Выявлены участки интервала интегрирования, где численная устойчивость метода накладывает ограничения на длину внутреннего шага интегрирования, что позволяет оптимизировать планирование экспериментов по исследованию закономерностей пластической деформации.
3. Шестернева О.В., Мальцева Т.В. Метод построения математической модели линейного динамического объекта // Вестник Сибирского государственного аэрокосмического университета. - 2009. - № 2. - С. 4-8.
4. Данейко О.И., Ковалевская Т.А., Колупаева С.Н., Семенов М.Е., Мелкозерова Н.А. Влияние скорости деформации на деформационное упрочнение и эволюцию дефектной подсистемы в гетерофазных материалах с ГЦК-матрицей // Известия вузов. Физика. - 2009. - Т. 52. - № 9/2. - С. 125-131.
5. Колупаева С.Н., Семенов М.Е. Исследование численного решения жесткой системы обыкновенных дифференциальных уравнений на примере модели пластической деформации // Прикладные задачи математики и механики: Матер. XVIII Междунар. научно-техн. конф., 13-17 сентября 2010 г. - Севастополь: Изд-во СевНТУ, 2010. - С. 15-18.
6. Колупаева С.Н., Семенов М.Е. Пакет прикладных программ для исследования пластической деформации скольжения в г.ц.к. материалах // Вестник Томского государственного архитектурно-строительного университета. - 2005. - № 1. -С. 36-46.
7. Ковалевская Т.А., Данейко О.И., Колупаева С.Н., Семенов М.Е. Математическое моделирование деформационного упрочнения сплавов, содержащих недеформируемые дисперсные частицы // Журнал функциональных материалов. -2007.- Т. 1. - №3. - С. 98-103.
8. Семенов М.Е., Колупаева С.Н. Свидетельство об официальной регистрации программы для ЭВМ № 20055612381 Slip Plasticity of Face-Centered Cubic v1.0 (SPFCC). Зарегистрировано в Реестре программ для ЭВМ 12 сентября 2005 г. Федеральная служба по интеллектуальной собственности, патентам и товарным знакам.
9. Семенов М.Е., Колупаева С.Н. Электронный информационный образовательный ресурс: «Программный комплекс SPFCC для исследования закономерностей пластической деформации в материалах с гранецентрированной кубической структурой» // Хроники объединенного фонда электронных ресурсов «Наука и образование», 2010. - № 10. URL: http://ofernio.ru/portal/newspaper/ofernio/2010/10.doc (дата обращения: 11.01.2011).
10. Колупаева С.Н., Ковалевская Т.А., Данейко О.И., Семенов М.Е., Кулаева Н.А. Моделирование температурной и скоростной зависимости напряжения течения и эволюции деформационной дефектной среды в дисперсно-упрочненных материалах // Известия РАН. Серия физическая. - 2010. - Т. 74. -№11. - С. 1588-1593.
11. Колупаева С.Н., Семенов М.Е., Пуспешева С.И. Математическое моделирование температурной и скоростной зависимо-
сти деформационного упрочнения ГЦК-металлов // Деформация и разрушение материалов. - 2006. - № 4. - С. 40-46.
12. Попов Л.Е., Кобытев В.С., Ковалевская Т.А. Пластическая деформация сплавов. - М.: Металлургия, 1984. - 182 с.
13. Семенов М.Е., Колупаева С.Н., Ковалевская Т.А., Даней-ко О.И. Математическое моделирование деформационного упрочнения и эволюции деформационной дефектной среды в дисперсно-упрочненных материалах // В сб.: Эволюция структуры и свойства металлических материалов / под ред. А.И. Потекаева. - Томск: Изд-во НТЛ, 2007. - С. 5-41.
14. Butcher J.C. Numerical Methods for Ordinary Differential Equations. - Chichester: John Wiley & Sons, Ltd, 2008. - 463 p.
15. Семенов М.Е., Колупаева С.Н. Анализ областей абсолютной устойчивости неявных методов решения систем обыкновенных дифференциальных уравнений // Известия Томского политехнического университета. - 2010. - Т. 317. - № 2. -С. 16-22.
16. Колупаева С.Н., Комарь Е.В., Ковалевская Т.А. Математическое моделирование деформационного упрочнения дисперсно-упрочненных материалов с некогерентной упрочняющей фазой // Физическая мезомеханика. - 2004. - Т. 7. - Спец. выпуск. - Ч. 1. - С. 23-26.
17. Старенченко В.А., Старенченко С.В., Колупаева С.Н., Пантю-хова О.Д. Генерация точечных дефектов в сплавах со сверхструктурой L12 // Известия вузов. Физика. - 2000. - Т. 43. -№1. - С. 66-70.
18. Семенов М.Е., Колупаева С.Н. Численное решение жестких систем обыкновенных дифференциальных уравнений при моделировании пластической деформации / Зимняя школа по механике сплошных сред (пятнадцатая). Сб. статей. В 3-х частях. Ч. 2. - Екатеринбург: УрО РАН, 2007. - С. 196-199.
19. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. - M.: Мир, 1999. - 685 с.
Поступила 11.01.2011 г.
