Федоренко Р.П. Введение в вычислительную физику (1185915), страница 42
Текст из файла (страница 42)
Начальные данные Ф» определяются по существу «флуктуационным» фоном, Первый этап процесса (сравнительно длительный) происходит практически при постоянной температуре а», потоки нейтронов экспоненцнально нарастают (Ф(1) ж с,е" ' р,), но пока еще слишком малы, чтобы привести к заметному изменению температуры. На этом этапе из всех компонент ряда (4), входящих в начальный фон Фо, выделяется первая. Постепенно Ф(1) достигает такого значения, что начинает изменяться и температура а(1), сначала медленно, потом все быстрее и быстрее.
Но одновременно начинает изменяться и )ь,. Характер зависимости коэффициентов оператора Е(а) таков, что Х, смещается при росте а в сторону отрицательных значе«сп ний (отрицательная обратная связь между а и !ь,). Темп роста Ф (и, следовательно, а) зв замедляегся. Наконец, Х~ становится отриЮ6» цательной, Ф начинает экспоненциально убывать, а а(1) выходит на некоторое стаци«з «з с «онарное состояние (около 1000 -ь 2000 'С). Рисунок 20 иллюстрирует сказанное вы~до ше. на нем представлены х,(е), и а(1) в точке максимума по пространственным переменным, Период начального разгона (от фона до значений Ф, вызывающих заметное изменение а) не рассчитывается.
Вместо этого находится р, (главная спектральная задача!) и в качестве начальных данных Ф» берется Д1~р,, где й' назначается на основании простых оценок А(а»). Такое «волевое решение» приводит к некоторому неопределенному сдвигу математического времени 1 относительно физического. Этот сдвиг допускает несложную оценку и для приложений не очень важен. Расчету подлежит только переходной процесс, длительность которого порядка 1 с. Опишем в общих чертах метод, который использовался в расчетах и позволил проинтегрировать уравнение (8) за 20 —: 30 шагов по времени. Правда, шаг (переход от состояния в момент ! к состоянию в момент 1+ т) предполагает выполнение достаточно сложной операции — решения главной спектральной задачи. В целом этот под- 5 161 гллвн»я спскт»Ассов»я злдАчл для кглввых з»длч эоз ход оказался эффективным именно благодаря тому, что был использован очень быстрый способ нахождения главной собственной функции (многосеточный метод; см.
3 !4). Основу метода составляет определение (разумеется, приближенное) искомого решения системы уравнений (11) из «асимптотиче- ского уравнения» ! —, = А(а)МоехР ~ Лс(С') асС' !Рс(с), а(0) = ао, о где зрс(с) — главная собственная функция спектральной задачи 14а(сН р= Лй'р; Л! — соответствующая точка спектра; Дсо — нормируюший множитель, определяемый «начальными данными» (при С= О множитель дс выбирается таким, чтобы уже началось незначительное, но заметное изменение температуры а). Опуская некоторые технические детали, стандартный шаг численного интегрирования можно представить следующими операциями. Пусть в некоторый момент времени с уже получены значения а(с) и ~ Л,(с') асс'.
Определяются коэффициенты оператора 1.(а), решается о главная спектральная задача, находятся !рс(!) и Л,(с). После этого для вычисления а(С + т) делается один шаг по явной схеме Эйлера: а(С + т) = а(С) + тАСа(С)1М» ЕхР ~ Л, о!С' Р,(С), о с+! ! пересчитывается интеграл ~ = ~ + тЛ,(с), и т.д. о о Естественно возникает вопрос о правомочности такого приближения.
Некоторые физически достаточно убедительные аргументы дает следующее рассуждение. Подставим используемую конструкцию в исходные уравнения задачи. Разумеется, они не будут выполнены: возникает некоторая «невязка». Если она очень мала относительно входящих в уравнение членов, то это в известной мере оправдает вышеизлос женный подход. Обозначая для удобства Л(с) = ~ Л, сй', вычисляем о с »с (а) оо ""с(с) дс елс!'! Л изр .1 ~ ' ца)„р у ехЯ й" пгивлижвнныв ми«оды вычислительной виэики 1ч. п Итак, погрешность метода зависит от соотношения й'д~р~/д~ и членов, входящих в Ь(а) 1рг (В этн выражения входят коэффициенты, погрешности измерения которых часто не так уж малы; кроме того, сама формулировка исходной задачи, которую мы принимаем за истину, в действительности основана на пренебрежении некоторыми относительно малыми величинами.) Заметим, что собственная функция ~р, определена с точностью до нормировки и этим фактором тоже можно разумно распорядиться с целью уменьшения погрешности.
Можно обыграть дза обстоятельства: величины о, и о существенно разные (скорости быстрых и медленных нейтронов связаны соотношением о~ ж 1Озоз). Следовательно, важнее получать медленное изменение второй компоненты 1рг Другое обстоятельство, которое можно было предвидеть на основе опыта расчетов в этой области: компоненты главной собственной функции ~р, (в зависимости от пространственных координат, которые в принятых обозначениях опущены) — просто устроенные функции «колоколообразной» формы.
С течением времени (при изменении а(г)) их форма меняется не очень существенно. Заметно меняется лишь отношение их амплитуд (т.е. соотношение между потоками быстрых и медленных нейтронов в холодном и горячем реакторе). Из сказанного выше следует рецепт нормировки: постоянной полагается гильбертова норма второй компоненты ~рг При численном решении задачи оценивалась и величина идгыр,!дп оиа фактически оказалась достаточно малой по сравнению с другими членами (около 1 —: 2 %). 5 17. Жесткие системы обыкновенных днфференцнапьных уравнений В пятидесятых годах при решении задач Коши для систем, описывающих кинетику реагирующих друг с другом химических веществ, вычислители столкнулись с крайне неприятным явлением.
Расчеты проводились с помощью хорошо отработанных программ, в которых использовались методы Рунге — Кутты и надежные алгоритмы автоматического выбора шага. Эти алгоритмы очень быстро вырабатывали шаг численного интегрирования, столь малый, что часто ие было никакой возможности рассчитать процесс на требуемом для приложений отрезке времени, даже используя наиболее мощные ЭВМ той эпохи. Визуальный анализ правых частей, казалось, не давал оснований для каких-то опасений.
Типичная система уравнений химической кинетики выглядела (технические подробности опускаем) примерно так: х'="', А!х~+~Х" А! х'х", й у, й 1, 2, ..., М, (1) » !7! жксткив системы одт где А,', А'.х — константы, характеризующие скорости протекания тех нли иных реакций. Бросалась в глаза, правда, существенная разница в их значениях: они отличались друг от друга часто на много порядков. В то же время исключить какие-то «малые» члены из (1), оставив только самые большие, было нельзя, ориентируясь лишь на значения А. Существенными были и концентрации разных веществ х'! они могли очень сильно изменяться с течением времени. Затрачивая значительное машинное время, удавалось получать начальные отрезки траекторий и провести анализ ситуации.
Он выявил следующую характерную картину. В начале процесса происходит сильное изменение х(() и выбираемый программой шаг численного интегрирования вполне разумен: он очень мал, но так и должно быть для аккуратного инте! рирования столь быстро меняющихся функций. Через небольшое времв г харак- х((! тер траектории резко меняется, она стано- и х х вится гладкой, медленно меняющейся, но программа этого «не замечает» и выбирает такой же малый шаг, Попытки «подска- ~Ч(»с~х((! зать» программе выбор существенно боль- Я$~ ч(~~~(~х шего шага, согласованного с гладкостью Хс(((7 ~ф~х(П решения, немедленно приводили к вычислительной катастрофе. В з 5, 7 специально указывалось„что прн оценке вычислитель- РиС.
2! ной сложности задачи Коши для системы обыкновенных' дифференциальных уравнений х = У(х) существенны два фактора: строение поля траекторий в окрестности интегрируемой траектории и свойства матрицы /„(х). Анализ поля направлений таких систем, получивших название «жестких», дал характерную картину, качественно представленную рис, 21. Траектория х(() состоит из короткого участка быстрого ее изменения (так называемого «пограничного слоя») и длительного участка очень медленной ее эволюции (иногда его называют «квазистационарным режимом»), Основные трудности связаны именно с расчетом последнего. Пограничный слой интегрируется очень малым шагом, но он настолько краток, что число шагов интегрирования вполне приемлемо. На рис.
21 при помощи «микроскопа» с последовательно увеличивающимся разрешением показана структура поля направлений в окрестности к((). Сначала видны траектории, отвесно падающие на х(г). При следующем увеличении видно, что, приближаясь к х((), они поворачивают, стремясь двигаться «параллельно» х(Г). И лишь при еще большем увеличении видна стандартная картина практически параллельных линий.
числвнныв мвтолы млтвмлтичвскоа оизлки 1ч и зш Если из точки х(г) траектории сдвинуться по касательной в точку х' = х(Г) + тх(1), то, хотя расстояние х(~ + т) — х' = О(тз) ничтожно, У(х ) не имеет ничего общего с х(» + с), направление г'(х') скорее напоминаег перпендикуляр к л, 1я1 траектории х(г). То же самое поо Лг лУчаетсл и пРи опРеДелении х' отрезком ряда Тейлора из трех- четырех н более членов при том о l о ЯсЛ значения т, которое хотелось бы использовать для численного интегрирования квазистационарною режима.