Курс лекций - Математическое моделирование технических объектов, страница 7

Следовательно, с целью повышения точности решений следует уменьшать шаг дискретизации по длине и по времени, что, в свою очередь, ведет к возрастанию количества основных циклов в алгоритме вычислений (пропорционально m x n), то есть – к увеличению времени счета.Другим важным свойством МКР является устойчивость решений. Устойчивость метода означает, что любые ошибки (округления или ограничения) не возрастают в ходе вычисления решения.

Доказано, что указанный процесс сходится и устойчив, если (<0.5) или(<1/[2]). Тем самым накладываются ограничения на выбор шага  по времени при выбранном шаге h разбиения вывода сечениями.Преимущества метода: (а) метод универсален и пригоден для области любой формы,причем источники тепла могут располагаться произвольно; (б) в задаче можно задавать любые граничные и начальные условия – алгоритм решения от этого не меняется; (в) в процессерешения получается полная картина распределения температуры внутри области на каждомшаге по времени, то есть возможен контроль допустимых температур; (г) возможно совершенствование и усложнение модели для получения более точной картины тепловых явлений.4.

Решение системы алгебраических уравнений на ЭВМПри использовании метода конечных разностей и методоа конечных элементов получается система линейных алгебраических уравнений, которая должна быть решена относительно неизвестных узловых параметров. Методы решения этой системы с малым илибольшим числом уравнений мало отличаются друг от друга. Специфика ее решения заключа-27ется в том, что при дискретизации D-области можно контролировать расположение коэффициентов в глобальной матрице жесткости.

В частности, можно получить матрицу ленточноготипа, которая имеет два полезных свойства: а) симметричность относительно главной диагонали и б) положительная определенность, означающая, что все коэффициенты на главнойдиагонали – положительные. Указанные свойства значительно сокращают объем вычислений.

Причем минимизируются ошибки округления. Ленточный характер матрицы позволяетзначительно сократить объем памяти для ее хранения.Одним из эффективных методов решения системы алгебраических уравнений, которыеполучаются при использовании МКЭ, является известный вариант метода исключения Гаусса. На первом этапе исходная матрица преобразуется к треугольному виду, после чего решение получается обратной прогонкой.Рассмотрим систему из 4-х линейных алгебраических уравнений вида:а+111х1 + а+112 х2 + а+113 х3 +а+121 х1 + а+122 х2 + а+123 х3 +а+131 х1 + а+132 х2 + а+133 х3 +а+141 х1 + а+142 х2 + а+143 х3 +а+114 х4 = b+11а+124 х4 = b+12а+134 х4 = b+13а+144 х4 = b+14Выразив из первого уравнения переменную x1, имеем:b+11-a+112-a+113-a+114 [1 x2 x3 x4 ] т[ a+111]+1+1+1a 11a 11a 11Подставляя полученное выражение для х1 во 2-е, 3-е и 4-е уравнения и приводяподобные члены, приходим к системе:x1=а+111х1 +000а+112 х2 +а+113 х3 +а+114 х4 = b+11+ (а22–а21[а12/а11])х2+(а23–а21[а13/а11])х3 +(а24–а21[а14/а11])х4 = b2–b1(а21/а11)+ (а32–а31[а12/а11])х2+(а33–а31[а13/а11])х3 +(а34–а31[а14/а11])х4 = b3–b1(а31/а11)+ (а42–а41[а12/а11])х2+(а43–а31[а14/а11])х3 +(а44–а41[а14/а11])х4 = b4–b1(а41/а11)В трех последних уравнениях все коэффициенты аpq и bp должны иметь верхний индекс (+1), поскольку эти коэффициенты взяты из исходной системы.

Далее указанный индексбудет использован для обозначения номера итерации решения исходной системы. Введемследующие обозначения:apq+(k+1)= apq+k- apk+k( akq+k/ akk+k) ; bp+(k+1)= bp+k- bk+k( apk+k/ akk+k)(14.1)Тогда последнюю систему можно переписать в виде:а+111х1 + а+112 х2 + а+113 х3 + а+114 х4 = b+110+ а+222 х2 + а+223 х3 + а+224 х4 = b+220+ а+232 х2 + а+233 х3 + а+234 х4 = b+230+ а+242 х2 + а+243 х3 + а+244 х4 = b+24Выразив из второго уравнения переменную x2, имеем:X2=[b+22a+222-a+223a+222-a+224a+222][1 x3 x4 ] тПодставляя полученное выражение для х2 в 3-е и 4-е уравнения и приводя подобныечлены, приходим к системе:а+111х1 + а+112 х2 +а+113 х3 +а+114 х4 = b+1128а+222 х2 +а+223 х3 +а+224 х4 = b+22(а33+2–а32+2 [а23+2/а22+2])х3 + (а34+2–а32+2 [а24+2/а22+2])х4 = b+23–b2+2(а32+2/а22+2)(а43+2–а42+2 [а23+2/а22+2])х3 + (а44+2–а42+2 [а24+2/а22+2])х4 = b+24–b2+2(а42+2/а22+2)Коэффициент при неизвестной х3 в третьем уравнении логично было бы обозначитькак а33+3.

Попробуем получить его формально, используя первую формулу в выражении(14.1). С этой целью обозначим p=3 (№ строки) , q=3 (№ столбца), k=2 (номер текущей итерации) и подставим эти индексы в (14.1):a33+(2+1)= a33+2- a32+2( a23+2/ a22+2) ;Получили очевидное совпадение результатов. Вычислим аналогично остальные коэффициенты при неизвестных в третьем и четвертом уравнениях:a34+(2+1)= a34+2- a32+2( a24+2/ a22+2) =a34+3a43+(2+1)= a43+2- a42+2( a23+2/ a22+2) =a43+3a44+(2+1)= a43+2- a42+2( a23+2/ a22+2) =a44+3Непосредственной проверкой можно показать, что правая итерационная формула, спомощью которой вычисляются свободные члены, так же верна. Проводя необходимые вычисления, получаем выражение для исходной системы уравнений после второй итерации выражения неизвестных:где: b3+3= b+23–b2+2(а32+2/а22+2) и b4+3= b4+2 –b2+2(а42+2/а22+2).После третьей итерации система примет вид:а11+1х1+ а12+1х2+ а13+1х3+ а14+1х4 =b1+10+а22+2х2+ а23+2х3+ а24+2х4 =b2+20+0+а33+3х3+ а34+3х4 =b3+3+4+40+0+0+а44 х4 = b4где: a44+4= a44+3–a43+3(а34+3/а33+3) и b4+4= b4+3 –b3+3(а43+3/а33+3).Решение полученной системы выполняем методом обратной прогонки.

Из четвертогоуравнения вычисляем неизвестную Х4:х4 =b4+4 /а44+4Из третьего уравнения вычисляем неизвестную Х3:х3 =[b3+3 -(а34+3х4)]/а33+3Из второго уравнения вычисляем неизвестную Х2:х2 + =[b2+2 – (а24+2х4 + а23+2х3)]/а22+2Наконец, из первого уравнения вычисляем неизвестную Х1:х1 =[b1+1 – (а14+1х4+а13+1х3+а12+1х2)]/а11+1На странице 27 приводится полная программа решения системы из n линейных алгебраических уравнений, в которой блок, реализующий метод обратной прогонки, выделенжирным шрифтом.Однако, непосредственно перед решением система должна быть преобразована, поскольку, как правило, некоторые компоненты неизвестного вектора Ф (в программе – этовектор Xr) узловых значений известны. Так, в большинстве задач теории поля некоторыеграничные значения искомой величины заданы; во всех задачах теории упругости должныбыть фиксированы некоторые перемещения с тем, чтобы исключить перемещение среды жесткого тела.То есть, элементы матриц [K] и {F} системы: [K]{Ф}={F} необходимо преобразовать, чтобы ее решение после преобразование давало правильный результат.

При этом жела-29тельно не изменять программу решения самой системы, поскольку это повлечет за собойтрудности при программировании.Пусть в исходной системе задано фиксированное значение р-й переменной (Хр=Q).Преобразование системы проводим по шагам: коэффициенты р-й строки, кроме диагонального коэффициента, равного аpp+1, приравниваем нулю; свободный член в р-й строке заменяем произведением: (аpp+1Q); уравнения, содержащие переменную Хр, преобразуем, вычитая из обеих частей каждого изних произведение (аqp+1Q), где q – номер строки (qp).Проиллюстрируем это на примере системы уравнений:46,6T1 – 21,7T2 +0+0= 1000- 21,7T1 + 93,2T2 – 21,7T3 +0= 2000(14.2)0- 21,7T2 + 93,2T3 – 21,7T4 = 20000+0– 21,7T3 + 56,6T4 = 1400Здесь, согласно условию задачи, фиксирована одна степень свободы узлового параметра {Т1=150}.

Преобразование системы проводим по шагам: коэффициенты 1-й строки, кроме диагонального коэффициента, равного К11=46.6, приравниваем нулю: 46,6T1+ 0 + 0 + 0 = 1000 свободный член в 1-й строке заменяем произведением: (К11Т1)=6990: 46,6T1+ 0 + 0 + 0 = 6990 переменная Т1 входит еще во второе уравнение, поэтому вычитаем из левой и правой части 2-го уравнения произведение К21Т1=(-21,7150):0 + 93,2T2 – 21,7T3 + 0 = 2000–(-3255)= 5255Таким образом, искомая система для решения примет вид:46,6T1 +0+0+0= 69900+ 93,2T2 – 21,7T3 +0= 52550- 21,7T2 + 93,2T3 – 21,7T4 = 20000+0– 21,7T3 + 56,6T4 = 1400что совпадает с системой из раздела 12.В программе решения системы уравнений, приводимой на стр.27, преобразование выполняется оператором:For i:=1 to n do If defX[i]=1 Then UppCase(i);Собственно преобразование выполняется подпрограммой UppCase, которая в качествепараметра принимает номер фиксированной степени свободы.

Последний выбирается из исходного линейного массива defX, в котором каждый фиксированный параметр должен бытьзаранее помечен единицей. Учебная pascal-программа преобразования и решения системы изn ЛАУ приведена в приложении:5. Метод конечных элементов (МКЭ)5.1. Типы конечных элементов.Используемые в настоящее время численные методы рассматривают ДУ непосредственно в той форме, в которой в которой они были выведены (без каких-либо дальнейших математических преобразований и манипуляций) при помощи: (а) аппроксимации дифференци-30альных операторов конечно-разностными алгебраическими операторами, действующими впоследовательности узлов, находящихся в области; (б) при помощи представления самой области элементами среды, не являющимися бесконечно малыми (то есть конечными элементами), которые в совокупности аппроксимируют реальную систему.

Наиболее громоздкой итрудно программируемой операцией в МКЭ является учет граничных условий задачи, причем точность полученного численного решения полностью зависит от степени измельченности сетки, определяющей узловые точки. Отсюда следует, что в процессе решения задачипрограммисту приходится иметь дело с системами алгебраических уравнений очень высокого порядка.В настоящее время наиболее популярным является второй подход, состоящий в возврате к характерному для физики разбиению ИТО на элементы конечных размеров, причем,чем больше по размерам эти элементы, тем лучше для минимизации числа получаемых уравнений. Поведение каждого элемента приближенно воспроизводит поведение малой областиИТО, которую он представляет.