Математическая модель и алгоритмы решения программного комплекса НЕ АТ-3 для расчета температурных полей

В настоящее время существует развитая теория разностных схем [1, 2], которая предлагает достаточно широкий выбор аппроксимаций для уравнения теплопроводности. При создании программного комплекса НЕАТ:3 для расчета температурных полей основными требованиями являлись универсальность программ и возможность реализации комплекса на ЭВМ малой и средней мощности. Под универсальностью понимается возможность решения задач в плоской, осесимметричной и объемной постановках при задании граничных условий I, II и III рода, многообразии граничных условий в одном расчетном узле, наличии в конструкции зазоров и прослоек, нелинейности свойств материалов, внутренних тепловыделениях и т. д. Под ЭВМ малой и средней мощности подразумеваются машины типа ЕС-1022, ЕС-1033.
В настоящее время основными методами численного решения уравнения теплопроводности являются метод конечных разностей и метод конечных элементов. При выборе класса методов авторы статьи руководствовались следующими соображениями относительно методики конечных элементов. Основой метода конечных элементов является минимизация некоторой интегральной величины, которая строится по исходному дифференциальному уравнению.
Таким образом, метод конечных элементов безусловно эффективен на простых областях. Если же область неоднородна (имеются зазоры, прослойки), то для корректного описания требуется дифференциальное уравнение другой структуры. Кроме того, небольшие геометрические размеры вызывают необходимость мелкого разбиения, а это, в свою очередь, увеличивает обусловленность матрицы системы уравнений и приводит к существенному увеличению времени расчета. Следует также указать на трудность подготовки исходных данных для объемной задачи. Однако самым существенным недостатком в рамках поставленной задачи является неэффективность вычислительных процедур МКЭ при решении нестационарных нелинейных задач, вызванная необходимостью расчета матрицы жесткости на каждом шаге.
Далее приводится достаточно универсальный алгоритм, базирующийся на методе конечных разностей, который позволил решить поставленную задачу.
Существенным является то, что при выводе уравнения (1) не использовались какие-либо упрощающие допущения о форме элемента, что позволяет использовать прямоугольные и треугольные сетки, единообразно описывать осесимметричные и объемные задачи. Таким образом, снимаются ограничения на геометрию расчетной области, свойственные ряду разностных схем.
Следовательно, можно задать граничные условия I и III рода. Существенно то обстоятельство, что интерполяционный метод позволяет единообразно описать достаточно широкий класс задач.
Учет лучистого теплообмена с окружающей средой можно провести, воспользовавшись методикой, описанной в работе [3] (приведение закона Стефана — Больцмана к граничным условиям III рода).
Для решения задач, связанных с термической обработкой, в рамках предлагаемой методики можно воспользоваться подходом, описанным в работе [1], т. е. использовать графики изменения свойств материалов в зависимости от температуры.
Таким образом, система (1) описывает практически весь спектр задач при определении температурного состояния, возникающих в энергомашиностроении.
Однако, несмотря на универсальность и другие положительные качества метода, решающими факторами являются возможность реализации «средней» задачи на конкретной ЭВМ и скорость ее решения.
Проанализируем последние результаты в области численных методов решения систем обыкновенных дифференциальных уравнений. Еще в 1952 г. в работе [4] было указано на существование противоречия между визуально медленным изменением компонент решений некоторых классов дифференциальных уравнений и очень малым шагом дискретности при численном интегрировании этих уравнений явными разностными схемами. Такое свойство систем дифференциальных уравнений называется «жесткостью».
Несмотря на обилие работ в этой области [5—7], устоявшейся терминологии в настоящее время не существует. Суть явления согласно современным представлениям заключается в наличии сильно различающихся собственных чисел в спектре матрицы Якоби
системы обыкновенных дифференциальных уравнений. В этом случае в( решении системы появляются компоненты с существенно разными временами затухания. Ограничение на шаг интегрирования для явных линейных методов, как известно нз работы [8], имеет характер. В случае, когда моделируемый интервал очень велик по сравнению с Н , это ограничение становится определяющим для возможности решения задачи за разумное время. Увеличить шаг интегрирования при затухании «быстрых» компонент невозможно из-за ошибок округления, неизбежно появляющихся при реализации метода интегрирования на ЭВМ.
Был предложен целый ряд методов для решения жестких систем (введено понятие «жесткая устойчивость», изобретены «жестко устойчивые» методы [5]). Самым универсальным и надежным средством для решения жестких систем уравнений следует, по-видимому, считать неявные разностные схемы.
Численные эксперименты показали, что система (1) для типовых задач, как правило, проявляет свойства жесткости, поэтому для ее решения используется неявный метод Эйлера. Для удобства дальнейшего изложения перепишем систему (1) в векторно-матричных обозначениях:
Однако, несмотря на то, что порядок аппроксимаций у этого метода выше, чем у неявного метода Эйлера, при ступенчатых изменениях граничных условий решение получается осциллирующим (на этот же факт указано в работе).
В настоящее время имеются стандартные программы решения жестких систем дифференциальных уравнений. Они, как правило, включают гребёнку методов с автоматическим переходом по какому-либо эвристическому алгоритму, однако рассчитаны на «малые системы».
Понятие «малая система», конечно, условно и зависит от типа применяемой ЭВМ. Будем считать систему малой» если #2<Мм0зу» где N — размерность системы уравнений, М — объем оперативного запоминающего устройства в словах. Другими словами, малая система — это система, обратную матрицу которой можно разместить в МОЗУ.
Так, в одной из распространенных программ для решения ОДУ [7] высокая эффективность достигнута, в основном, благодаря двум факторам: обратная матрица типа (М—НА) вычисляется не на каждом шаге, а лишь в случае плохой сходимости итераций и хранится в МОЗУ; применяется схема «предиктор — корректор», так что предсказанное значение вектора неизвестных приводит к малым невязкам системы типа (3), и число итераций для уточнения решения невелико.
Применительно к системе (2) для больших N эффективность метода прямо связана со скоростью решения системы алгебраических уравнений (3).
Для решения системы (3) возможно использование как прямых, так и итерационных методов. Матрица (М—НА) обладает теми же свойствами, что и матрица А, а именно, диагональным преобладанием и симметричностью, что позволяет использовать прямые методы решения линейных систем, основанных на идее исключения Гаусса.
В программном комплексе НЕАТ-3 используется итерационный метод. Выбор этого класса методов линейной алгебры обусловлен стремлением снизить требования к минимальному объему памяти, который в методе Гаусса определяется шириной ленты матрицы.
Для организации итерационного метода используется разложение матрицы (М—НА) вида
(М—Н А)=Э3+{М—Н Л—Я3),
где 1>3 — трехдиагональная матрица. Таким образом итерации производятся по формуле:
1-й шаг ОаТ1к%\ = (д3-М + М)г|+1 +
МТъ + НЬ (ть+з.) = Ь;
где верхний индекс означает номер итерации.
Для интерпретации этого метода с физической точки зрения следует рассмотреть алгоритм нумерации элементов. Если задача одномерная, то при последовательной нумерации матрица (М-*-Н А)=03 и итерации (4) сходятся за один шаг. При принятых алгоритме разбиения и нумерации исходной области (разбиение на слои; последовательная нумерация слева направо в каждом слое и сверху вниз) итерационный процесс (4) можно интерпретировать как расчленение области на ряд одномерных задач с учетом их взаимодействия по предыдущей итерации.
Подобный алгоритм используется в работе [9], но только при подготовке и хранении исходных данных. Такая нумерация позволяет также экономным образом организовать программы обработки результатов.
Практическая проверка алгоритма показала, что для большинства задач при использовании «предиктора» второго порядка точности процесс сходится за 3—4 итерации.
На ЭВМ ЕС-1022 с емкостью МОЗУ 512 кбайт программный комплекс НЕАТ-3 позволяет довести число элементов расчетной схемы до 3000. Применение комплекса в течение 2 лет для решения широкого круга задач показало его надежность.