Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (1185899), страница 24
Текст из файла (страница 24)
Однако при записи линейной системы уравнений всем неизвестным надо присвоить сквозную нумерацию и представить их в виде одномерного массива — вектор- столбца. Такая перенумерация позволяет представить систему разностных уравнений в общепринятой матричной форме записи систем линейных алгебраических уравнений и воспользоваться стандартными программами их решения. Выполним перенумерацию по горизонтальным прямым слева направо и снизу вверх. В этом случае неизвестные нижней горизонтальной прямой обозначаются иь и», ..., иэ, неизвестные второй горизонтальной прямой — ичю, / / ! ..., и/»/ и т.
д. Пример такой перенумерации показан на рис. 3.!4. Общая формула пересчета индексов л, т двумерного массива в индекс й одномерного массива имеет вид и„=и„, й = — (т — 1) й/+л. После введения новой нумерации формируется матрица А и столбец свободных членов Р линейной системы Ах) =-Р. 12 Рис. ЗЛ4 Если рассмотреть какую-либо й-ю строку матрицы А, то входящие в нее ненулевые коэффициенты будут соответствовать температурам, участвую«цим в записи уравнения теплового баланса для Ьй элементарной ячейки. Из введенной нумерации и схемы !3.76) вытекает, что в балансе 0-й внутренней ячейки участвуют только температуры узлов с номерами й — 1, й+ 1, й — Ж, й+ А!. В балансах граничных ячеек участвует еще меньше температур. А это означает, что в любой строке матрицы А коэффициенты, отстоящие от диагонали далее чем на Ф, равны нулю: аи ... а„»„0 0...0 а„„' ... а„„...
а„„,»0... 0 аии,ии-и аии,ии Пример заполнения матрицы для области с девятью узлами, представленной на рис. 3.14, показан на рис. 3.15, в котором символом «х» отмечены ненулевые коэффициенты. Таким образом все ненулевые коэффициенты лежат в ленте шириной 2 М + 1, осью которой является главная диагональ матрицы. В пределах ленты также могут встречаться нулевые коэффициенты. Для прямоугольной области структура матрицы внутри ленты является упорядоченной, однако если рассмотреть область более сложной формы, например изображенную на рис.
3.12, то ленточный характер матрицы сохранится, но ненулевые коэффициенты будут более сложным образом расположены внутри ленты. Отметим, что матрица А является симметрич- «З40078 Рис. 336 !!6 ной, и поэтомудостаточноформировать и хранить в памяти машины только верхнюю часть ленты. Симметрия матрицы вытекает из условия согласования тепловых потоков между соседними элементарными объемами в методе баланса. Мы проводили перенумерацию чпо горизонталям». Можно было организовать ее и аппо вертикалям», двигаясь снизу вверх и слева направо. Тогда ширина ленты равнялась бы 2 М+1. Очевидно, что для зкономии машинной памяти, а для многих алгоритмов и машинного времени выгоднее иметь минимальную ширину ленты.
Поэтому следует проводить перенумерацню вдоль той координаты, по которой берется меньшее число узловых точек. Таким образом, при использовании неявной схемы сначала в соответствии с принятой перенумерацией неизвестных проводится формирование ленты матрицы А и столбца свободных членов, а затем с помощью обращения к какой-либо стандартной программе решения линейной системы уравнений с ленточной матрицей находятся искомые значения температуры. Пример использования такой стандартной программы рассматривается в следующей главе применительно к системе уравнений метода конечных элементов, которая также имеет ленточный вид. Мы рассмотрели методику решения линейной задачи.
В случае наличия каких-либо нелинейностей применяют такие же приемы, как и для одномерной задачи (см. й 3.6), и приходится на каждом шаге по времени решать линейную систему уравнений (при наличии итераций — решать многократно) с меняющейся от шага к шагу (и от одной итерации к другой) матрицей. Основной проблемой при реализации описанного подхода является быстрый рост затрат машинного времени с увеличением числа узловых точек в области.
Например, при использовании специальных модификаций метода Гаусса для ленточных матриц число арифметических операций для решения системы уравнений пропорционально КЕ», где К вЂ” общее число узловых точек в области, равное числу неизвестных в системе, Š— ширина ленты матрицы.
Особенно неприятно это для нестационарных нелинейных задач. С целью сокращения затрат машинного времени были разработаны конечно-разностные схемы, у которых эти затраты на каждом шаге по времени пропорциональны числу узловых точек К. Такие схемы называются экономичными. Из экономичных схем, получивших распространение на практике, рассмотрим в следующем параграфе локально-одномерную (24). К ее достоинствам относятся безусловная устойчивость, возможности применения как для двух-, так и для трехмерных задач.
При решении многомерных стационарных задач применяют два подхода. При первом составляется и решается система конечно-разностных уравнений для стационарной задачи. Эта система получается из (3.76) — (3.78) обнулением левых частей уравнений (так 117 дТ как — = О). Структура матрицы системы стационарных разностдт ных уравнений полностью совпадает с рассмотренной выше для неявной схемы, поэтому для решения стационарной задачи требуется число арифметических операций, пропорциональное КЕ«.
В ряде случаев (особенно для трехмерных задач) такой подход приводит к большим затратам машинного времени и памяти, резко возрастающим с увеличением числа узлов К. Другой подход, называемый «счетом на установление», заключается в определении решения стационарной задачи путем моделирования процесса выхода в стационарный режим нестационарного температурного поля, которое рассчитывается по какой-либо экономичной разностной схеме. При этом приходится делать определенное число шагов по времени. Общие затраты машинного времени равны произведению числа шагов по времени l на затраты на одном шаге.
При использовании экономичных схем затраты на расчет поля на одном шаге пропорциональны числу узлов сетки К. Поэтому общие затраты времени с увеличением числа узлов растут медленнее, чем при решении стационарной системы с ленточной матрицей. Кроме того, при счете на установление нет необходимости хранить в памяти матрицу А, содержащую ьК элементов. $3.$. лОкАльнО-ОднОмеРнАя схемА Локально-одномерная схема является «типичным представите- лема широкого класса схем, применяемых для решения многомерных задач и задач расчета совместно протекающих процессов, описываемых несколькими уравнениями (например, уравнениями теплопроводности и диффузии или уравнениями Навье — Стокса и энергии для потока жидкости).
Отличительная особенность этих схем — сочетание сильных сторон явных схем (малые затраты машинного времени на шаге по времени) и неявных схем (безусловная устойчивость). В таких схемах протекание многомерного физического процесса на каждом временном шаге представляется как результат последовательной реализации соответствующих одномерных процессов„ каждый из которых начинается от распределения поля, возникшего после окончания предыдущего одномерного процесса. На основе такого представления, называемого ропцеплением задачи по пространственным переменным, моделирование одномерных процессов проводится с помощью неявных схем, а последовательное действие процессов учитывается по существу явным образом, т.
е. решение многомерной задачи сводится к расчету на каждом шаге по времени набора одномерных задач, решаемых в случае уравнения теплопроводности методом прогонки. Применение неявной аппроксимации одномерных задач обеспечивает устойчивость схемы, а общее число арифметических действий оказывается пропорционально числу Нв дб д«д ч« (3.79) б(х, у, т««) = Т (х, у, тд,), (3. 80) дв д«»» г В д«, т««(т<тм т у (3.81) ш (х, у, тэ,) = б (х, у, тз). (3.82) Пока мы предполагаем, что приближенное решение начинается из точного распределения Т (х, у, т;,). Отметим, что сначала решается уравнение (3.79) с начальным условием (3.80), а затем уравнение (3.81), в котором в качестве начального условия принимается полученное к концу временного интервала распределение б (х, у, т;).
Граничные условия для уравнений (3.79), (3.81) соответствуют граничным условиям исходной задачи по направлениям х и у. ыа узловых точек, поскольку алгоритм прогонки обладает этим свойством (см. ~ 3.4). Аналогичный подход используется и для задач расчета нескольких совместно протекающих процессов, в которых на каждом временном шаге расщепление проводится по физическим процессам, т. е. последовательно решаются отдельные уравнения со своими граничными условиями, а значения величин, определяемых из других уравнений, берутся из уже полученных на данном или предыдущем временном шаге полей.
После расщепления по физическим процессам отдельные многомерные задачи можно далее расщеплять и по пространственным координатам. Описанная методика внешне весьма проста, однако ее конкретное воплощение, связанное с решением вопросов о том, «что можно и как нужно расщеплятык во многих случаях связано с большими трудностями и требует от расчетчика высокой математической квалификации и хорошего понимания физики исследуемого процесса. Рассмотрим обоснование допустимости одного из вариантов расщепления и соответствующую этому варианту локально-одномерную схему применительно к задаче (3.73) — (3.75) при равномерной по пространственным координатам сетке с шагами Й„и Ь».
Для простоты покажем допустимость расщепления уравнения теплопроводности по пространственным переменным при достаточно малом шаге по времени без учета дискретизации пространственной области. С этой целью сопоставим точное решение уравнения (3.73) в конце малого промежутка времени [тэ „т,[ с его приближенным решением ш (х, у, т;), получаемым в результате решения на этом же промежутке времени следующей системы, возникающей при расщеплении: (3.84) 120 Покажем, что для точного решения Т (х, у, тз) и решения ш (х, у, т;) расщепленной задачи справедливо соотношение ш(х, у, т;)= Т(х, у, т;)+0(Лта).
(3. 83) Из (3.73) имеем Т(х,у,тЗ) — Т(х,у,тт,) (д'Т д~Т 1! — 1 до и, следовательно, Т (х, у, тз) = Т (х, у, т;,) + ГдТ дт и — Дтя, +Лта ~ —, + —,~ + — + 0(Лт'), Аналогично из (3.79), учитывая (3.80), получаем Г д'Т 17 — ~ атсо б(х, У,т)) = — Т(х, У,тт,)+ Ата( д, ) -1- +0(Атт). (385) Далее нз (3.81) и (3.85) находим Г д~д ~/ ш(х у тз) =б(х у тз)+Лта( д з 1 д' Т У ' Л'~~~о = Т(х, у, т;,)+Л а( дх, 1 + — '+ +Ат~( д ) +(Ата) ( д д ) +0(Лт'). (3.86) Сравнивая (3.86) и (3.84), убеждаемся в справедливости соотно- шения (3.83). Таким образом, при условии использования в начале интерва- ла (тз „тз) точного распределения (3.80) в конце интервала появ- ляется погрешность порядка 0 (Лт').
Однако при прохождении все- го промежутка [О, т„„х] с шагом Лт только на первом шаге можно будет взять точное начальное распределение, а на последующих шагах в качестве начального распределения для б (х, у, т) на дан- ном шаге придется использовать распределение ш (х, у, тз,) в кон- це предыдущего интервала. Тогда вместо условия (3.80) следует рассмотреть условия вида: при 1 = 1 б (х, у, О) = Т, (х, у), при 1 ) 1 б (х, у, т; 1) = ш (х, у, т) т) (3.8?) При расчете по приближенной методике (3.79), (3.87), (3.81), (3.82) будет происходить накопление погрешности.
Однако накоп- ленная погрешность не может превзойти суммы пошаговых погреш- ностей, поскольку для уравнения теплопроводности искажение решения возмущением начального распределения температуры всегда меньше этого возмущения. Общее число шагов по времени равно т ,„!Лт, и таким образом накопленная к моменту време- ни т „. погрешность составляет 0 (Лт») т,„!Лт — 0 (Лт), т. е.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
