Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (1185899), страница 30
Текст из файла (страница 30)
Однако в реальных двумерных (и тем более трехмерных) задачах число узлов и элементов может составлять несколько сотен, а иногда и тысяч, и поэтому построение расчетной сетки вручную и ввод больших массивов чисел в качестве исходных данных нецелесообразны из-за значительных затрат времени на их подготовку и большой веРоятности появления ошибок. Следовательно, возникает задача автоматизации процедуры разбиения области на элементы, нумерации элементов и узлов и формирования индексной матрицы. При эгом требуется в качестве входной информации для соответствую'цей подпрограммы задавать сравнительно небольшое число данных, описывающих геометрию области сложной формы и густоту сетки, а на ее выходе получать массивы координат узлов и индекс- 147 ной матрицы, позволяющие перейти к формированию глобальной матрицы системы уравнений МКЭ.
Для областей типа прямоугольника, треугольника или некоторых более сложных, но конкретных конфигураций можно составить частные подпрограммы, реализующие заполнение массивов координат узлов н индексной матрицы на основе данных о размерах области и шагах по осям х, у. Например, можно покрыть область прямоугольной сеткой, а затем построить треугольники путем разбиения элементарных прямоугольников диагоналями, как это сделано на рис.
4.1!. Программирование таких процедур сводится к организации нескольких циклов с переменными пределами. 2 Ф х Рис. 4.11 Рис. 4.12 Однако для практических приложений больший интерес представляют универсальные программы автоматического разбиения областей различной сложной формы. В литературе предложены различные способы описания геометрии и алгоритмы дискретизации областей для МКЭ. Наибольшее распространение получил следующий подход. Область сложной формы разбивается вручную на подобласти, которые называются «макроэлементамн».
Эти подобласти должны достаточно хорошо описывать геометрию расчетной области. Обычно макроэлементы выбирают в форме треугольников и выпуклых четырехугольников, но иногда используют и подобласти, ограниченные кривыми второго порядка. Число таких макроэлементов обычно невелико (несколько единиц или десятков для сложных областей), и поэтому это разбиение можно описать путем задания координат узлов макроэлементов и некоторой условной нумерации макроэлементов и их узлов. Затем реализуется автоматическое разбиение каждого из макро- элементов на элементарные треугольные элементы. Для этого в исходных данных лишь указывают параметры, характеризующие густоту сетки в каждом из макроэлементов.
Приведем два примера автоматизации разбиения. В первом используются макроэлементы в виде выпуклых четырехугольников (рис. 4. (2). Для каждого макрочетырехугольника задается следующая информация: номера узлов в вершинах и признаки принадлежности сторон границе, числа дроблений на отрезки по двум 148 смежным сторонам, массивы координат вершин четырехугольников.
Каждый из макрочетырехугольников автоматически разбивается на маленькие четырехугольники прямыми, проходящими через точки разбиения противоположных сторон на равные отрезки. Полученные четырехугольники в свою очередь разбиваются короткими диагоналями на элементарные треугольники.
В результате число треугольников в каждом макроэлементе равно 2 Й,я«, где й„й,— числа отрезков на смежных сторонах. Заметим, что при таком способе для общих сторон макроэлементов должна задаваться одинаковая кратность дробления. Другой способ автоматизации разбиения иллюстрирует рис. 4.13. Здесь в ка- «=1 честве макроэлементов взяты треугольники. Информация о них задается почти в таком же виде, как и для «элементар- х 0 ных» треугольников: массивы координат «=г вершин н индексная матрица, но с одним отличием.
В строке индексной матрицы для каждого макроэлемента содержится еще одно число — кратность дробления й. Если я =- О, то макроэлемент не дробится и принимается в качестве конечного элемента. При й = ! путем соединения центров сторон проводится разбиение макроэлемента на четыре подобных треугольника (рис. 4ДЗ). При й = 2 каждый из полученных четырех треугольников еще раз разбивается на четыре подобных и т. д., т. е. число полученных из макроэлемента треугольников равно 4".
Кратность дробления соседних макроэлементов может различаться не более чем на единицу. При этом, чтобы избежать появления «лишних» узлов на стороне треугольника с меньшей кратностью дробления автоматически проводится построение еще нескольких треугольников. Для этого узел, лежащий на стороне треугольника, соединяется с противоположной вершиной, как это показано на рис. 4.13 пунктирными линиями. Достоинством данного способа разбиения является возможность резко сгущать сетку в областях с большими градиентами температур, используя при этом сравнительно небольшое число макроэлементов. Приведенные на «словесном» уровне описания методик автоматического разбиения двумерной области могут создать впечатление о простоте их программной реализации. Однако на самом деле программное «воплощение» этих алгоритмов требует довольно значительных усилий.
В процессе автоматического разбиения области проводится и нумерация получившихся элементов разбиения. После этого возникает задача оптимальной перенумерации узлов с целью уменьшения ширины ленты матрицы. В областях, имеющих форму, близкую к 149 прямоугольной, нумерацию узлов целесообразно проводить последовательно вдоль направлений, параллельных короткой стороне и содержащих меньшее число узлов. В данном случае ситуация аналогична той, что была рассмотрена ранее (см. 5 3.7) применительно к решению двумерной задачи методом конечных разностей. В некоторых программах нумерацию проводят последовательно в пределах каждого макроэлемента в порядке их обхода.
При этом в треугольниках, лежащих у границ макроэлементов, могут возникать большие разности номеров узлов. Можно предложить и другие подходы. Например, хорошо зарекомендовал себя для решения мно- гих практических задач следующий 12 эвристический способ перенумерации узлов. Узел с номером 1 выбирают где-либо на границе области.
Номера 7 2, 3 и т.д. присваивают узлам, кото- рые соединены с узлом номер 1 сто- 10 роной элементарного треугольника б (рис. 4.14). После того, как пронумерованы все узлы, соединяющиеся с узлом номер 1, в качестве «опорного» г берется узел номер 2 и продолжается нумерация узлов, соединенных с ним Рис.
4л4 и не пронумерованных ранее. Подоб- ная процедура последовательной нумерации у ~лов, соединенных с наименьшим из ранее присвоенных номеров, продолжается до полной нумерации всех узлов сетки. После такой перенумерации подсчитывается получившаяся ширина ленты матрицы, которая также является выходным параметром программы разбиения и перенумерации. Поскольку программная реализация описанных автоматических процедур довольно сложна, примеры соответствующих подпрограмм здесь не приводятся. Формирование глобальных матрицы и вектор-столбца. Перейдем к рассмотрению реализации алгоритмов формирования глобальных матрицы и вектор-столбца.
Обычно они строятся на основе методики, изложенной в предыдущем параграфе. В этом случае оформляется подпрограмма для расчета локальных матриц и столбцов. Входными данными для нее являются координаты узловых точек, передаваемые в соответствии с локальной нумерацией узлов, информация о сторонах элемента, лежащих на границе области, и значения параметров Х„ Хч, д„ д,,а. Последние могут вычисляться с помощью специальных подпрограмм, в которых входными параметрами являются координаты узлов элемента, а выходными — значения Х„Ха, д„, д„и для этого элемента.
В приводимом ниже примере (рис. 4.15) предусмотрено использование двух таких подпрограмм: для расчета Х„, Х„, д, (подпрограмма ААА) 150 и а, д, (подпрограмма ВВВ). Обращения к этим подпрограммам имеют вид: СА1.1. ААА(ХС, УС, АХ, АУ, ЯЧ) и СА1.1. ВВВ(ХС, Ъ'С, А).г, ЯЯ) Эти подпрограммы должны составляться пользователем для конкретной задачи. При формировании матрицы происходит обращение к ним с текущими координатами ХС, УС центра элемента (для ААА) или центра стороны элемента (для ВВВ). Прием, аналогичный описанному, использовался в главе 3 при решении одномерной задачи теплопроводности методом конечных разностей.
Напомним, что при выводе уравнений МКЭ мы считали свойства и мощности постоянными в пределах элемента, поэтому в случае разрывных функций желательно, чтобы линии разрыва совпадали с границами элементов. Входными данными для подпрограммы формирования глобальных матрицы и столбца, приведенной на рис. 4.15, являются: И— число треугольных элементов, М вЂ” число узлов, МЬ вЂ” ширина ленты матрицы, Х, У вЂ” массивы координат узлов х„, у длиной М, 1й10 — индексная матрица — массив длиной 4эМ. Все эти данные получаются в результате выполнения процедуры разбиения и перенумерации узлов. В индексной матрице для каждого треугольного элемента (и = — 1, ..., У) записываются четыре целых числа. Первые три числа соответствуют глобальным номерам 1, 1', й узлов данного треугольника, записанным в порядке обхода против часовой стрелки, начиная с любой вершины.
Четвертое число — это признак, с помощью которого задается внешняя граница области. Напомним, что согласно принятому нами способу этот признак принимает значения: 0— если ни одна из сторон треугольника не принадлежит границе, !в если граничной является сторона между узлами 1 н 1, 2 — если граничной являются две стороны между узлами 1 н 1 и между узлами 1 и Ф. Другие ситуации с расположением граничных сторон можно исключить соответствующим выбором порядка записи номеРов 1,1, й, В приводимой подпрограмме из-за малой размерности локальных матрицы н столбца для расчета их коэффициентов не используется специальная подпрограмма, а они вычисляются непосредственно в ней по формулам (4.29), (4.30) (операторы 83 — 122).