Полак_и_др__Вычислительные_методы_в_химической_кинетике (972296), страница 59
Текст из файла (страница 59)
ВТ(РРС (И, Т, У, Н, НЕ, ЕР5, МР), где И; МР— целые, Т, Н, НЕ, ЕРВ— действительные, У вЂ” одномерный массив из И элементов. Подпрограмма осуществляет численное интегрирование системы из И обыкновенных дифференциальных уравнений 1-го порядка, Т вЂ” текущее время, перад обращением к подпрограмме этой переменной должно быть присвоено ' Программа разработана А.А. Левицким. 237 начальное значение текущего времени, при выходе иэ подпрограммы Т содержит значение времени, соответствующее концу интервала интегрирования. У вЂ” значения интегрируемых функций.
На входе в подпрограмму в этот массив засылаются начальные условия для всех переменных, на выходе из подпрограммы в У' содержатся решения системы ОДУ для соответствующего момента времени. Н на входе в подпрограмму должно содержать значение начального шага интегрирования. Так как в начале интегрирование системы ОДУ осуществляется методом 1-го порядка, то рекомендуется задевать Н достаточно малым, в дальнейшем подпрограмма 8Т!РР увеличивает порядок метода, "разгоняется", и значение шага интегрирования быстро растет.
На выходе в Н содержит я значение текущего шага интегрирования. НŠ— конечный момент времени, до которого необходима интегрировать систему ОДУ. ЕРБ — относительная точность интегрирования. МР— параметр, определяющий особенности решения уравнения для корректора. Если МР= 20, то при решении используется "метод с неподвижной точкой", МР = 21-вычисление якобизма осуществляется по аналитическим формулам, запрограммированным в подпрограмме РЕОЕВЧ, МР = 22 — вычисление якобиана производится путем чиа ленного дифференцирования, МР= 23 — вместо якобиана при решении уравнения для корректора методом Ньютона-Рафсана используется только главная диагональ якобиане, вычисленная численным дифференцированием [67].
..~ Подпрограмма 8Т1РРС вызывает подпрограммы 8Т1РР и 1РЙЗН Пакет 87!ГЕ, являющийся переработкой пакета 6ЕАВ [282], подробно описан в работе [67]. 8 приводимой здесь версии пакет 8Т1РР несколько сокращен за счет исключения не употрабляемых в данном случае подпрограмм для !ибаты с ленточными матрицами и исключения Ваэмахечостм интегри. рования методом Адамса. Перечислим лишь функции подпрограмм пакета 871РР: СО8ЕТ вычисляет коэффициенты метода интегрирования соответствующего порядка; ВЕЗСА(. осуществляет перемасштабирование решений системы и их производных при изменении шага интегрирования; МАТОЕС осуществляет (. (7-разложение матрицы Š— (36 А, где А — якобиан системы, Š— единичная матрица, И вЂ” текущий шаг интегрирования; ((— константа (см. главу Я]; МАТ801 решает соответствующую систему линейных уравнений РВЕО1С вычисляет значение прадиктора; подпрограмма ЗТ1РР управляет процессом интегрирования.
Кроме подпрограммы 8Т(РГ, подпрограмма 8Т(РРС вызывает целую функцию!ГЙ8Н (1у, ТР, Н, У. ЛТАВТ],где (у, .БТАВТ-целые, ТР, Н— действительные, У вЂ” двумерный массив, (У- жало уравнений,ЮТАВТ— порядок метода интегрирования, ТР— текущее значение времени, Н вЂ” текущее значение шага интегриравения, первый столбец массива У содержит значения решений системы ОДУ в данный моэаент времени, второй — нормированные на шаг значения первых производных решений и тд. (ЕЙ8Н служит для управления печатью решений в процессе интегрирования. Для вывода значений решений в заданные моменты времени производится интерполяция с использованием сохраняемых нормированных значений производных вектора решений.
Прх желании печать можно организовать не для определенных моментов времени, а для опраделанмых значений концентрации той или иной компоненты. В случае, воли требуется продолжать интегрирование 1РЙ8Н, надо присвоить значение О; если по каким. либо причинам необходимо выйти иэ 8Т! РРС, то! РЙЕН следует присвоить значение 1. ,О1РР0Й ((у, Т, У. Р], где (у — целое, Т-действительное, У, Р— одномерные массивы длиной (У. Эта подпрограмма вычисляет значения правых 238 Реакция СН, -СН, +Н СН, +Н -СН, -Н, 5 СН,+СН, С Н,, 4 С Н, .С,Н„+Н, 5 С,Н. ° СН, +СНЗ а с',н, -О.Й, +н, 7 С,Н,.-С~.С+Н, 9 н+Й+н,-н,+н, Аг я, Ег, икая/мель 13,7 0 93,0 13,7 0 12,9 129 0 0 13,9 0 69,0 14,7 0 79.3 9,5 0 40,0 0,2 0 г 30,0 19В -1 0 Исходные данные к программе задаютса в следующем виде. 1.
И, М, М7., ИТМ, МЕ,. Н, ЕРЗ, ТК в формате (613, ЗЕ8.11. Дла рассматриваемого примера перфокарта исходных денных имеет вид 0080080080040211.00Š— 101.00Š— 023000.000. Число компонент И = 8, число реакций М = В.М7. целесообразно выбирать так, чтобы ЭВМ оперировала со значениями концентраций порядка единицы. В процессе решения истиннье'значения концентраций умножаются на 10ьг ', поэтому следует взать М7. б. Число точек, в которых необходимо печатать решения, ИТМ = 4. Так как при аналитическом вычислении акобиана режим работы программы обеспечивает максимальное быстродействие, то МЕ= 21. Начальный шаг И следует выбирать малым по сравнению с характерными временами происходящих процессов, Н = 10 удовлетворяет этому частей системы ОДУ по значениам времени Т и решений У, соответствующих этому моменту времени, и засылает их в массив Е.
РЕОЕВ17 (И, Т, У, УУУОIМ!, где И, ИУОIМ вЂ” целые, Т вЂ” действительное, У вЂ” одномерный массив длиной И. Эта подпрограмма вычисляет значения акобиана интегрируемой системы ОДУ по значениям времени Т и Решений У, соответствующих этому ыоменту времени, и засылает элементы я кобиана в массив общего блока ВТСОМ4.
Обмен информацией в программе осуществляетса через общий блок КIИЕТ, переменные которого имеют следующие значения: С вЂ” массив, содержащий значения констант скорости химических реакций; РŠ— массив. содержащий нормированные значения натуральных логерифьюв пред- экспонентов констант скорости; ТИ вЂ” массив, содержащий показатели степени в температурных множителях констант скорости; ЕА — массив, содержащий значения энергий активации констант скорости; гу — массив, содержащий значения скоростей реакций. вычисляемых в процессе решения; 7. 77 — целый массив, содержащий коды химических реакций. Каждан реакция кодируется девяткой целых чисел: первое — число веществ в тевой части уравнения хиьмческой реакции, второе — число веществ в п1хгвой части этого уравнения, далее следуют новара веществ, участвующих в реакции, записанные слева направоУР— рабочий массив, используемый для печати; ч77И вЂ” массив, содержащий значения моментов времени, в которые необходимо печатать решение; ТК вЂ” температура, К; ТЕМ— температура, нкап/мольгчАМŠ— масштабный множитель;УИ вЂ” число ком.
поиенг в кинэтической схеме; М вЂ” число реакций; УМ7. — десятичный логарифм АМ11 7ТМ вЂ” текущее значение индекса массива ТМ; !СИ вЂ” целый г/ массив, содержащий наименования компонент кинетической схемы. Для иллюстрации задания исходных данных к программе рассмотрим в качестве примера следующую модельную задачу: рассчитать концентрации компонент СН4, СгН4, СгН4, СгНг, Нг, С, СНз. Н для значений време.
ни г = 10 ", 1О, 10 г, 1 с. Начальная концентрация метана 10 1 моль/сыт, начальные концентрации остальных компонент равны нулю, температура пповедения реакции 3000 К. Кинетическая схема: условию. Как правило, для большинства кинетических схем достаточной является относительная погрешность на шаге интегрирования ЕРЕ= 10 ', однако можно проводить два расчета, например с ЕРЯ= 10 з и 10 з, и сравнивать результаты. Температура, для которой проводится расчет, ТК = 3000 К.
2. Наименования компонент, присутствующих в кинетической схеме, перфорируются по 10 штук на карте в формате (1ОАВ) . Дпя рассматриваемого примера все наименования компонент умещаются на одной карте:- СН4 С2Н6 С2Н4,, С2Н2 Н2 С нСНЗ Н 3. Записи уравнений химических реакций пербюрируются по одной карте на реакцию в формате (ЗГБ.2,7 (2А4,А1). Последовательно задаются: десятичный логарифм предзкспонента, показатель степени в температурном множителе константы скорости реакции, знергия активации в ккал/моль, символическая запись реакции. Для рассматриваемого примера необходимо отперфори ровать следующие перфокарты: 13.
700000093. ООСН4 =СНЗ „, +Н 13. 700000012. 9ОСН4» з +Ни и н н =СНЗ»и+Н2 12.900000000000СНЗ н н» +СНЗ»»».»»» =С2НБ 13. 900000069. ООС2Н6» ° =С2Н4 н +Н2 14. 700000079. ЗОС2НБ»»».» =СНЗ ч»» ~ н» +СНЗ 08. 600000040. ООС2Н4»» н =С2Н2 н н»» +Н2 06. 200000030. ООС2Н2 и» н» =Сини и ни +Снниинн н+Н2 18.80 — 1.
0000000Н» ванн» н н+Н» ин» ч»н +Н2» ни л»и = Н2»и» нн +Н2 4. Начальные значения концентраций задвотся в формате (10Е8.1! по 10 чисел на карте. Для рассматриваемого примера данные должны быть отперфорированы в виде 1. ООЕ-ОБ Б. Значения моментов времени, в которые следует печатать решения также задаются в формате (10ЕВ. 1), по 10 чисел на карте. Дпя рассматриваемого примера данные должны быть отперфорированы в виде 1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
