Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (1185899), страница 29
Текст из файла (страница 29)
Второе уравнение системы имеет вид дУ(з> д1(з> + — =О. ди, дц, >42 и'," =и, 1 >з> и, ==и„ ц(з> =ц 4 Глобальные номера ! 2 3 2 4 3 4 5 3 и'," =и, з ц(м ц з и'," =и . 3' Из (4.27), (4.28) вытекает, что в локальных обозначениях неизвестных это уравнение записывается следующим образом: дД> и',»+д>» и,'»+у<»и'," — ч><»+ +дЯ> и','>+у<~> и',"+д<з> и,'" — >р<~>= О. Обратим внимание, что при записи д>'>'>1ди, использованы коэффициенты первой строки локальной матрицы и первый коэффициент столбца для второго элемента, поскольку температура и, в его локальных обозначениях имеет первый номер: и, = и<'.
Теперь для получения искомого уравнения необходимо в соответствии с (4.32) заменить локальные обозначения неизвестных на глобальные; а>> ">+к з иэ+к2з ив — >р>»+д>>',>из+ +ь>1'*~ и4+к>'-;> и,— >р>1з>= О. (4. ЗЗ) Группируя слагаемые с одинаковыми неизвестными и складывая постоянные коэффициенты, получаем вторую строку глобальной матрицы С и второй коэффициент глобального вектор-столбца Ф в виде С2> (й>1 у>2 +Уй Юйз +д1з Ю1> ~ О) (4.34) Ф >р>» > >р>з> Аналогичным образом можно получить и остальные три уравнения. Алгоритм формирования глобальных матрицы и вектор-столбца.
Полученные выражения (4.34) позволяют изложить принцип формирования л>-го уравнения глобальной системы. Это формирование целесообразно проводить путем постепенного суммирования вкладов от различных элементов. При машинной реализации перед началом формирования массивы, в которых помещаются глобальные С н Ф, обнуляются, а затем к их текущим значениям постепенно добавляются соответствующие коэффициенты локальных матРиц и столбцов.
Ясно, что вклад в т-е уравнение системы дадут только те элементы, у которых в строке индексной матрицы имеется номер и>. Если >и-й узел числится в локальной нумерации какого- либо из этих элементов 1-м (1= 1, 2 или 3), то будет использована 1-я строка локальной матрицы й<"> и 1-й коэффициент локального вектор-столбца «р<" >. Найденный нужный коэффициент локального столбца прибавляется к текущему значению л>-го коэффициента глобального столбца. Коэффициенты выделенной строки локальной матрицы элемента прибавляются к соответствующим коэффициентам и>-й строки глобальной матрицы, имеющим порядковые номера, указанные в строке индексной матрицы, т. е.
первому коэффициенту строки локальной матрицы соответствует первый номер «отсылкнэ в строке индексной матрицы, второму коэффициенту — второй номер, третьему — третий. 143 Описанная процедура лежит в основе алгоритма формирования глобальной матрицы и глобального вектор-столбца. Как было уже отмечено выше, она реализуется путем последовательного перебора элементов следующим образом. Берется первый элемент, анализируется его строка в индексной матрице и устанавливается, в какие три уравнения этот элемент «дает вклад». Далее рассчитываются его локальная матрица и вектор-столбец.
При этом расчете используется информация о наличии у данного элемента граничных сторон, содержащаяся в четвертом столбце индексной матрицы. Пусть локальным номерам 1, 2, 3 соответствуют фактические номера ~', /, й. Тогда первая строка локальной матрицы и первый коэффициент локального вектор-столбца участвуют в формировании ~-й строки глобальной матрицы и 1-го коэффициента глобального вектор- столбца.
Производится сложение найденных локальных коэффициентов йк',), д(,!>, йкф с имеющимися значениями глобальных коэффициентов 6;о 66, 6;ы Затем аналогичная процедура повторяется для второй и третьей строк локальной матрицы и второго и третьего коэффициентов локального столбца. Они участвуют в формировании строк глобальной матрицы и коэффициентов глобального столбца с номерами ~ и й, которые соответствуют локальным номе- рам2иЗ. Изложенный на примере треугольных элементов разбиения метод формирования глобальных матрицы и вектор-столбца, основанный на введении локальной нумерации узлов и неизвестных, легко переносится и на случай более сложных элементов разбиения.
Он является наиболее общим, часто используемым и тем более эффективным, чем сложнее применяемые конечные элементы. Свойства системы разиостных уравнений и методы ее решения. Теперь рассмотрим ряд важных свойств, которыми обладает глобальная матрица. Во-первых, можно доказать, что она является симметричной. Во-вторых, глобальная матрица для задач большой размерности М является сильно разреженной, т. е.
большинство ее элементов — нулевые. Наконец, путем введения разумной нумерации узлов ее можно сделать ленточной. Остановимся на двух последних свойствах матрицы. Очевидно, что коэффициент 6; в гп-й строке глобальной матрицы отличен от нуля, только когда узлы с номерами гп и / являются вершннамн какого-то общего для них элемента. В этом случае в строке индексной матрицы, соответствующей этому общему элементу, будут содержаться номера т и ~. Указанное обстоятельство и объясняет разреженность глобальной матрицы, посколь ку, например, для треугольных элементов при значительном числе треугольников Л' большинство возможных пар узлов т и й не являются вершинами общего треугольника и, следовательно, соответствующий элемент глобальной матрицы 6 „— — О.
Рассмотрим влияние нумерации узлов на структуру глобальной матрицы С. Из сказанного выше вытекает, что расположение нуле- 144 вых элементов в матрице зависит от способа нумерации узлов. Например, в рассмотренном выше конкретном примере при нумерации, указанной на рис. 4.8, глобальная матрица С выглядит так: б„ 6„ бзз 0 0 633 643 0 0 633 бзз бы 0 бзз бзз 634 635 643 643 644 645 0 бзз бм бм (4.35) Если же перенумеровать узлы так, как это сделано на рис. 4.9, тс матрица С примет вид 6„ 0 0 — 651 643 о о б„ 6,3 63, О 63, 63.
634 634 635 О 63 бм б„ б 653 654 655 (4.36) й =шахмат — я(, оо где т, й — номера узлов и-го треугольника. Схематичный вид глобальной ленточной матрицы показан на Рис. 4.10. Символами «х» обозначены ненулевые коэффициенты. ~се коэффициенты, расположенные за пределами полосы, ограни- (4. 37) 145 Общее число нулевых элементов в (4.36) не изменилось по сравнению с (4.35), однако в (4.35) ненулевые элементы расположены лишь на главной диагонали н на двух прилегающих к ней верхних и нижних диагоналях, а в (4.36) эти элементы «разбросаны» по всей матрице.
Таким образом при разумной нумерации узлов глобальная матрица С имеет (з/ ленточный вид, т. е. все ненулевые коэффициенты расположены в пределах полосы, И/ образованной рядом верхних и нижних диагоналей, примыкающих к главной дна- (О гонали. Из симметрии матрицы следует, г что число верхних и нижних диагоналей с отличными от нуля коэффициентами оди- Рис, 4.9 наково. Поскольку для треугольного разбиения коэффициент 6„3 отли' чен от нуля только в случае, когда узлы т и й принадлежат одно- мУ треугольнику, то положение наиболее удаленного от главной диагонали ненулевого элемента матрицы определяется мгксимальной по всем парам общих вершин треугольников разностью номеров узлов, т.
е. величиной ченной линиями, параллельными главной диагонали, равны нулю. В общем случае нулевые коэффициенты встречаются и внутри полосы. Число (Е+1) будем называть шириной полосы (шириной ленты) матрицы. Ленточный характер и симметрия глобальной матрицы позволяют значительно сократить объем машинной памяти, требуемой для ее хранения.
При программировании задачи предусматривают запись матрицы не в виде массива длиной М х М, а в виде массива, содержащего лишь элементы, находящиеся в пределах полосы на главной диагонали и выше. Например, если требуется решить задачу с числом неизвестных узловых темпеЕ+! ратур М=ЗОО, то для записи мат- рицы в общем виде необходимо хра- хох ОЬО нить 300х300=9 10' вещественных ххахаао чисел. Пусть ширина ленты в этой заООхххх О даче равна Е+1 = 40. Тогда запись матрицы в сокращенном виде потребует массив длиной (Е + 1) М == О х О х х х О =- 1,2 1О' элементов, если матрицу О х х х х О х запоминать в виде «прямоугольника» с Е + 1 столбцами и М строками, О О О Оххх, т.
е. для облегчения логики программирования предусматривать место для хранения фиктивных элементов в последних Е строках (заштрихованный треугольник вне матрицы на рис. 4.10). Если же не учитывать фиктивные элементы «хвоста», то потребуется запомнить [(Е + 1) х х М вЂ” Е(Е + 1)!2) =-: 1,12 1О' чисел. Наконец, перейдем к вопросу решения системы уравнений. Для решения систем уравнений МКЭ применяют как прямые, так и итерационные методы. Причем последние обычно используют в тех случаях, когда объем оперативной памяти не позволяет хранить всю глобальную матрицу даже с учетом ленточного симметричного вида. Из прямых методов хорошо зарекомендовал себя на практике и получил широкое распространение метод квадратного корня.
Этот метод пригоден только для систем линейных уравнений с симметричной матрицей и по затратам машинного времени примерно вдвое быстрей метода исключения Гаусса. В математическом обеспечении ЭВМ имеются стандартные программы, реализующие метод квадратного корня. Предусмотрен и случай систем с ленточной матрицей (стандартная подпрограмма МСНВ из математического обеспечения ЕС ЭВМ [15]). В заключение подчеркнем, что использование той или иной стандартной подпрограммы решения системы уравнений требует определенного способа записи глобальной матрицы в одномерный массив. Применяемые способы различны для разных подпрограмм, т.
е. может организовываться запись по 146 строкам или по столбцам, с учетом или без учета «хвоста» ленты. Вто обстоятельство следует учитывать при программировании алго- ритма формирования глобальной матрицы. Е 4АХ ПРОГРАММНАЯ РЕАЛИЗАЦИЯ МКЭ Существенным достоинством МКЭ является возможность составления программ численного расчета полей в областях сложной геометрической конфигурации, которые проще по логической структуре и по заданию исходных данных, чем программы, реализующие метод конечных разностей для таких областей. В данном подразделе рассмотрим в качестве примера структуру программы для решения двумерной задачи (4.1), (4.2) в областях произвольной формы при треугольных элементах разбиения.
В программе решения задачи методом конечных элементов выполняются следующие основные процедуры: 1. Разбиение области на элементы, нумерация элементов, глобальная и локальная нумерации узлов и формирование на их основе индексной матрицы. 2. Формирование глобальных матрицы и вектор-столбца системы алгебраических уравнений, реализуемое на основе расчета локальных матриц и столбцов.
3. Решение системы разностных уравнений. 4. Расчет температур и тепловых потоков в различных точках элементов разбиения, проводимый на основе принятой аппроксимации температурного поля в элементе. Остановимся на особенностях программной реализации первых двух процедур, рассматривая их применительно к решению задачи (4.1), (4.2). Автоматизация разбиения области. Простейший (но наиболее трудоемкий) способ реализации первой процедуры состоит в ручном разбиении области О на треугольные элементы, ручной нумерации Узлов и дальнейшем вводе в качестве исходных данных массивов координат узлов (х ) ы (у„,), и индексной матрицы.