Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Сложности решения системы (8) связаны с тем, что W — это малый «векторный» параметр, т.е. система является сингулярно-возмущенной, описывающей взаимодействие процессов с существенно разными характерными временами. Если, например, для переменных Ф время 0.1 -г-1с является большим, то для температуры а вре-
206
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
мя 0.1 -г- 1 с является малым. (Во избежание недоразумений отметим, что все конкретные цифры относятся к тем расчетам, результаты которых будут приведены ниже.)
Рассчитываемый процесс был связан с проработкой проекта научно-исследовательского реактора, в котором при «холодном» состоянии (а0 «0 0C, Ф° « 0) внезапно (поднятием регулирующих стержней) создается надкритическая ситуация (X1« 10 -г 200). Начальные данные Ф° определяются по существу «флуктуационным» фоном. Первый этап процесса (сравнительно длительный) происходит практически при постоянной температуре а0, потоки нейтронов экспоненциально нарастают (Ф(0 « CjeVip1), но пока еще слишком малы, чтобы привести к заметному изменению температуры. На этом этапе из всех компонент ряда (4), входящих в начальный фон Ф°, выделяется первая.
Постепенно Ф(0 достигает такого значения, что начинает изменяться и температура а( t), сначала медленно, потом все быстрее и быстрее. Ho одновременно начинает изменяться и X1. Характер зависимости коэффициентов оператора L(a) таков, что X1 смещается
при росте а в сторону отрицательных значений (отрицательная обратная связь между а и X1). Темп роста Ф (и, следовательно, а) замедляется. Наконец, X1 становится отрицательной, Ф начинает экспоненциально убывать, a a(t) выходит на некоторое стационарное состояние (около 1000 -і- 2000 °С).
Рисунок 20 иллюстрирует сказанное вы- 4 ше. На нем представлены X1 (г), и а(() в точке максимума по пространственным переменным. Период начального разгона (от фона до значений Ф, вызывающих заметное изменение а) не рассчитывается. Вместо этого находится ^1 (главная спектральная задача!) и в
качестве начальных данных Ф° берется JV^1 , где N назначается на
основании простых оценок Л(а°). Такое «волевое решение» приводит к некоторому неопределенному сдвигу математического времени t относительно физического. Этот сдвиг допускает несложную оценку и для приложений не очень важен. Расчету подлежит только переходной процесс, длительность которого порядка 1 с.
Опишем в общих чертах метод, который использовался в расчетах и позволил проинтегрировать уравнение (8) за 20 -ь 30 шагов по времени. Правда, шаг (переход от состояния в момент t к состоянию в момент t + т) предполагает выполнение достаточно сложной операции — решения главной спектральной задачи. В целом этот под-
Рис. 20
§ 16]
ГЛАВНАЯ СПЕКТРАЛЬНАЯ ЗАДАЧА ДЛЯ КРАЕВЫХ ЗАДАЧ
207
ход оказался эффективным именно благодаря тому, что был использован очень быстрый способ нахождения главной собственной функции (многосеточный метод; см. § 14).
Основу метода составляет определение (разумеется, приближенное) искомого решения системы уравнений (11) из «асимптотического уравнения»
da
Tt
= A(Ci)Ar0 ехр
•Фі(0, а(0) = ос°,
где ^p1(Z) — главная собственная функция спектральной задачи
Lfct(Z) ] Ч>= №
A1 — соответствующая точка спектра; N0 — нормирующий множитель, определяемый «начальными данными» (при Z = 0 множитель N0 выбирается таким, чтобы уже началось незначительное, но заметное изменение температуры а). Опуская некоторые технические детали, стандартный шаг численного интегрирования можно представить следующими операциями.
Пусть в некоторый момент времени Z уже получены значения a(Z)
t
и I Ai (t') dt'. Определяются коэффициенты оператора L(a), решается
о
главная спектральная задача, находятся Tp1 (0 и A1(J)- После этого для вычисления a( t + т) делается один шаг по явной схеме Эйлера:
a(Z + т) = а(/) + T^4[a(Z)]-/V0 ехр
$ A1 dt'
Естественно возникает вопрос о правомочности такого приближения. Некоторые физически достаточно убедительные аргументы дает следующее рассуждение. Подставим используемую конструкцию в исходные уравнения задачи. Разумеется, они не будут выполнены: возникает некоторая «невязка». Если она очень мала относительно входящих в уравнение членов, то это в известной мере оправдает вышеизло-
t
женный подход. Обозначая для удобства A(Z) = ^ A1 dt', вычисляем
о
а
at
-L(a)
N0eW) гр, (0 =
= 7У0ел(<)
ац>!
A1^ip1 +% -Jt--------Ца)гр1
208
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
Итак, погрешность метода зависит от соотношения WdylZdt и членов, входящих в L(a)yr (В эти выражения входят коэффициенты, погрешности измерения которых часто не так уж малы; кроме того, сама формулировка исходной задачи, которую мы принимаем за истину, в действительности основана на пренебрежении некоторыми относительно малыми величинами.) Заметим, что собственная функция Tp1 определена с точностью до нормировки и этим фактором тоже можно разумно распорядиться с целью уменьшения погрешности. Можно обыграть два обстоятельства: величины U1 и и2 существенно разные (скорости быстрых и медленных нейтронов связаны соотношением U1 « IO2U2). Следовательно, важнее получить медленное изменение второй компоненты Tp1.