CC BY
94
УДК 539.3
АЛГОРИТМ РЕАЛИЗАЦИИ ЗАДАЧ ТЕОРИИ УПРУГОСТИ И ПЛАСТИЧНОСТИ ВАРИАЦИОННО-РАЗНОСТНЫМ МЕТОДОМ. ЧАСТЬ II
В.Н. Барашков
НИИ ПММ при Томском государственном университете E-mail: ger@mail.tomsknet.ru
Представлена численная методика определения дву- и трехмерного упругопластического напряженно-деформированного состояния твердого деформируемого тела с помощью вариационно-разностного метода, реализующего вариационный принцип Лагранжа методом конечных разностей. Физические соотношения принимаются согласно теории малых упругопластических деформаций, а геометрические соотношения берутся в виде уравнений Коши. Физически нелинейная задача решается методом переменных параметров упругости. На примере задач о деформировании тела вращения оживальной формы и цилиндрического сектора проводится сравнение прямого и итерационного методов решения системы линейных алгебраических уравнений большого порядка, к которой приводит использование необходимого условия экстремума сеточного аналога функционала полной потенциальной энергии системы.
4. Оценка времени решения СЛАУ итерационными и прямым методами
Проведем оценку времени счета системы уравнений (4) итерационными и прямым методами, а также отметим недостатки и преимущества этих методов. В литературе существует достаточно много работ, в которых рассматривается этот вопрос. Причем анализ времени решения проводится зачастую в контексте количества арифметических операций, необходимых для реализации одного итерационного шага или решения СЛАУ прямым методом.
Методом Гаусса"... плотная система ^уравнений с N неизвестными может быть решена с помощью примерно №/3 умножений плюс другие арифметические операции, которыми обычно пренебрегают при грубых оценках количества операций" [1].
В работе [2] отмечается, что методом Зейделя "... реализация одного итерационного шага осуществляется за 2АР-М арифметических действий (М-1 сложенней, М-1 умножений и одно деление). Если в каждой строке матрицы А отлично от нуля лишь т элементов, а именно эта ситуация имеет место для сеточных эллиптических уравнений, то на реализацию итерационного шага потребуется 2тМ-Марифметических действий, т.е. число действий пропорционально числу неизвестных М'. Следует заметить, что такая оценка верна лишь для случая т<М, что бывает достаточно редко.
В работе [3] считается, что"... Число ^арифметических операций, необходимых для реализации метода Гаусса, определяется следующей формулой [4, 5]: ЛГ=2и(л+1)(и+2)/3+й(и-1), где п - число неизвестных. Таким образом, время, необходимое для выполнения арифметических операций при решении СЛАУ методом Гаусса, примерно пропорционально кубу числа неизвестных". И далее у Н.В. Коп-ченовой и И.А. Марона читаем:"... Если метод итераций сходится, он дает следующие преимущества по сравнению с методами, рассмотренными выше (имеются в виду прямые методы):
1) Если итерации сходятся достаточно быстро, т.е. если для решения СЛАУ требуется менее п итераций, то получаем выигрыш во времени, так как число арифметических действий, необходимых для
одной итерации, пропорционально л2, а общее число арифметических действий в методе Гаусса, например, пропорционально и3.
2) Погрешность округления в методе итераций сказываются значительно меньше, чем в методе Гаусса. Кроме того метод итераций является самоисправляющимся, т.е. отдельная ошибка, допущенная в вычислениях, не отражается на окончательном результате, так как ошибочное приближение можно рассматривать как новый начальный вектор.
3) Метод итераций становиться особенно выгодным при решении систем, у которых значительное число коэффициентов равно нулю.
4) Процесс итераций приводит к выполнению однообразных операций и сравнительно легко программируется на ЭВМ".
В работе [1] также отмечается тот факт, что итерационные методы"... в отличие от прямых методов имеют тенденцию быть самокорректирующимися и, следовательно, минимизируют ошибки округления".
H.H. Калиткин [6] пишет:"... Для систем небольшого порядка п < 200 применяются практически только прямые методы'. Итерационные методы выгодны для систем специального вида, со слабо заполненной матрицей очень большого порядка и »103... 105. Метод Гаусса с выбором главного элемента выгоден для систем общего вида с плотно заполненной матрицей. Он требует примерно л2 ячеек в оперативной памяти. При вычислениях производится ~2л3/3 арифметических действий; из них половина сложений, половина умножений и п делений".
Сделать подобные оценки для МВР с ОКР не представляется возможным, ибо сходимость итерационного процесса зависит от выбранной величины Можно лишь с уверенностью сказать, что время решения СЛАУ этим методом меньше, нежели итерациями Зейделя.
В итоге можно сделать вывод о том, что время вычисления методом Зейделя пропорционально л2, а методом Гаусса - и3. В случае, если итерационный процесс сходится за fc итераций (k<n), время вычислений в методе Зейделя становится меньше. Но эти
данные о реализации итерационного и прямого методов относятся кплотно заполненной квадратной матрице (их«). Использование ВРМ в задачах механики деформируемого твердого тела приводит, как говорилось выше, к симметричной положительно определенной матрице ленточной структуры с шириной ленты И. В этом случае необходимо ориентироваться на другие оценки относительно времени реализации решения СЛАУ. Для системы п уравнений с /г членами в каждом время одной итерации будет пропорционально л/г. Если решение с заданной точностью достигается за к итераций, тогда время итерационного метода решения пропорционально кпк Практически не бывает, чтобы итерационный метод дал решение с удовлетворительной точностью за количество итераций к меньшее /г, ибо шаг ленты И для двумерной осесим-метричной задачи зависит от количества узлов сетки по толщине конструкции, которое в большинстве случаев меньше аналогичных параметров вдоль длины тела. Метод Гаусса для симметричной матрицы ленточной структуры дает решение за время, пропорциональное пИ2/2 [7,8].
Для сравнения методов решения системы уравнений (4) была решена квазистатическая задача определения осесимметричного упругопластического НДС образованного двумя усеченными эллипсоидами толстостенного тела вращения оживальной формы при действии массовых сил большой интенсивности [7]. Дискретизация расчетной области проводится с помощью четырехугольных ячеек. Методом Зейделя для со = 1 и метод верхней релаксации с оптимальным коэффициентом релаксации (йор, = 1,9445 (для /=250, Д=50) была просчитана система алгебраических уравнений со 120 неизвестными и шириной ленты А=10 матрицы {А}. Судя по количеству проделанных итераций, необходимых для выполнения заданной точности вычисления перемещений 8 = 0,01 (см. табл. 1), об ускорении сходимости итерационного процесса говорить не приходится. Но результаты значительно различаются между собой.
Верным из них является тот, при котором величина энергии 3 имеет меньшее значение, ибо ищется ее минимум. О близости к истинному решению также можно судить по степени выполнения теоремы Клапейрона, согласно которой
Р = 2, (9)
или соотношения
э^и-^+^и-ги^и, (Ю)
получаемому подстановкой (9) в выражение для функционала энергии (2). МВР с в>ор1 дает решение ближе к истинному, нежели решение, полученное методом Зейделя.
В табл. 2 приведены результаты исследования сходимости решения СЛАУ МВР с выбором ОКР в зависимости от точности е для 500 неизвестных и ширины лены й= 12 матрицы {Л}. Здесь введены следующие обозначения: t - время счета задачи в условных временных единицах (УВЕ), необходимое
для получения решения с точностью s; ^ - среднеквадратичная невязка решения системы уравнений. Во всех трех вариантах (аор1= 1,9487 (/= 600, Д = 50). МВР с ОКР значительно ускоряет процесс поиска решения СЛАУ. Для сравнения был проведен расчет НДС методом Зейделя. Для достижения точности s = Ю-4 потребовалось 360 УВЕ (« 22000 итераций) вместо 45 УВЕ (2687 итераций) при использовании МВР с аор1.
Время счета сокращается примерно вдвое, если коэффициенты ау матрицы {А} не вычислять на каждой итерации, а хранить в оперативной памяти ЭВМ во время решения каждой упругой задачи. Пересчет же коэффициентов проводить только для вновь полученных величин функции пластичности A.A. Ильюшина (т.е. для решения следующей упругой задачи). Но в этом случае теряется одно из преимуществ итерационных методов перед прямыми методами - возможность реализации СЛАУ большого порядка.
Свободным от недостатков, связанных с большими затратами машинного времени, является метод Гаусса решения СЛАУ, с помощью которого для получения аналогичных результатов понадобилась 1 УВЕ (см. нижнюю строку табл. 2). Из таблицы также следует, что теорема Клапейрона (9) и соотношение (10) выполняются точно, а для определения упругопластического НДС тела вращения потребовалось сделать четыре приближения.
Таким образом, результаты численного эксперимента дают право утверждать, что для реализуемых ВРМ двумерных упругопластических задач, время решения СЛАУ МВР с аор1 меньше, нежели итерациями Зейделя и больше, чем методом Гаусса. Сравнивая результаты решения задачи определения НДС тела вращения двумя методами (МВР и Гаусса), обнаруживаем совпадение величин перемещений до четырех значащих цифр. Уменьшая величину погрешности £ в МВР при решении упругопластичес-
Таблица 1. Результаты решения задачи методом Зейделя и
MBPc(i>opl
U к Слагаемые функционала энергии, Дж
э и (4+4)/2
1,0000 682 -6001 3891 4955
1,9445 692 -6290 6060 6200
Таблица2. Результаты решения задачи МВР с (Oopt дляраз-личных значений точности е
Е к t, УВЕ El Слагаемые функционала энергии, Дж
э и {А + Аг)П
10_i 1508 27 IO"6 -6300,7 6148,2 6224,5
Ю-4 2687 45 10"' -6302,4 6305,1 6303,5
Ю-4 3855 66 10^ -6302,4 6302,3 6302,5
10J 4 1 10^ -6302,4 6302,4 6302,4
кой задачи, можно добиться большего совпадения результатов. Итак, для задач данного класса метод Гаусса с точки зрения времени реализации является предпочтительнее по сравнению с итерационными методами.
Анализ методов и времени решения СЛАУ для трехмерной задачи теории упругости проводился на решении задачи о деформировании цилиндрического сектора, представляющего собой треть толстостенного цилиндра с внутренним отверстием (труба) и нагруженного на торце равномерным давлением.
Следует отметить резкое увеличение объемов вычислений (как вспомогательных - тех же коэффициентов ав матрицы {А}, так и при решении СЛАУ) для трехмерной задачи по сравнению с двумерной задачей. Получающиеся при варьировании перемещений в одном узле три уравнения содержат подлежащие вычислению 243 коэффициента и три свободных члена, содержащие внешнюю нагрузку и граничные геометрические условия. В двумерном же случае для четырехугольных ячеек система двух уравнений для варьируемого узла имеет 36 коэффициентов и два свободных члена. Ширина к ленточной матрицы {А} СЛАУ зависит от порядка нумерации (направления обхода) узлов конечно-разностной сетки при составлении сеточных уравнений и для трехмерной задачи гораздо больше величины к для двумерной задачи, а размер матрицы даже на грубой сетке значительно превышает аналогичную величину для двумерной задачи. В этой ситуации в некоторой степени выручает свойство матрицы решаемой СЛАУ - ее симметрия относительно главной диагонали. Поэтому количество коэффициентов а., подлежащих определению в одном варьируемом узле при использовании метода Гаусса, уменьшается до 123.
Вследствие симметрии матрицы {А} при реализации СЛАУ (4) методом Гаусса использовались лишь ее главная диагональ и верхняя треугольная матрица. Это позволило почти вдвое сэкономить необходимую для матрицы память и решить методом Гаусса систему уравнений на сетке (г'х_/><&) = = (5x9x10) относительно л=1350 неизвестных с шириной ленты А=189 за 4 УВЕ. Назовем это решение точным (см. табл. 3, N=2).
При решении СЛАУ итерационными методами были использованы метод Зейделя и МВР с выбором соор1. В качестве критерия окончания итерационного процесса брались условия выполнения с за-
данной точностью е5 теоремы Клапейрона (9) и интегрального соотношения (10), или равенства искомых величин перемещений в двух соседних итерациях. Применение МВР с ОКР по сравнению с методом Зейделя не дало сколько-нибудь заметного выигрыша во времени решения, который имел место при реализации двумерных задач.
Время решения задачи при использовании итерационных методов существенно зависит от выбора нулевого приближения по перемещениям и задаваемой точности. Так, для получения величин деформаций и напряжений с погрешностью 6 % на той же сетке потребовалось выполнить 29 итераций за 1 УВЕ (табл. 3, N=1). Для вычисления этих же величин с погрешностью 0,5 % при удачно выбранном нулевом приближении по перемещениям было проделано уже 102 итерации за 3,5 УВЕ (N=2).
Таким образом, время решения для данной конкретной задачи теории упругости в трехмерной постановке прямым и итерационными методами примерно одинаково. В приведенных выше результатах упоминалось об относительно небольшом количестве итераций, необходимых для определения параметров НДС с заданной точностью. Практически одинаковые затраты машинного времени объясняются многократным вычислением одних и тех же коэффициентов ау (их количество для рассмотренной СЛАУ составляет порядка 110000) матрицы {А} системы уравнений при решении итерационными методами, в то время как при использовании метода Гаусса вычисление этих коэффициентов делается только один раз при наборе матрицы. При запоминании коэффициентов а9 время решения системы уравнений итерационным методом значительно уменьшается.
Описанную выше задачу о деформировании цилиндрического сектора удалось реализовать методом итераций на больших сетках: (9x9x10), я=2430, к=297 {N=4) и (9x9x19), и=4617, А=297, (#=5), что с помощью метода Гаусса сделать не удалось в связи с недостатком оперативной памяти. При использовании метода Гаусса только для хранения главной диагонали и верхней треугольной матрицы системы уравнений для варианта N=5 потребовалось бы 1327293 ячейки памяти.
Следует отметить, что хотя варианты ЛГ=1, 4 и считались МВР (как и все остальные, кроме N=3) с определением ОКР, но вследствие относительно большой величины заданной погрешности вычисления итерационный процесс для них завершился
Таблица 3. Резуль тэты решения трехмерной задачи на разных сетках с различной точностью
N со к УВЕ Слагаемые функционала энергии, Дж
Э и {А+Аг)!2
1 1,0000 0,10 29 1,0 -2210,769 2471,163 2340,967
2 1,6931 0,01 102 3,5 -2218,455 2240,017 2229,236
3 - - - 4,0 -2218,514 2218,514 2218,514
4 1,0000 0,05 71 4,6 -2221,114 2339,614 2275,864
5 1,8049 0,05 110 15,0 -2221,024 2329,915 2275,470
до момента определения со^ Поэтому в табл. 3 для этих вариантов ш = 1 и СЛАУ решалась, по сути дела, методом Зейделя.
Для реализации СЛАУ методом Гаусса в дву- и трехмерных задачах использовалась программа, описанная в [9], которая подвергалась корректировке в каждом конкретном случае.
Анализируя полученные численные результаты, можно сделать некоторые выводы. К перечисленным выше из работы [3] достоинствам итерационных методов следует добавить также следующие моменты:
1) они не требуют обязательного хранения матрицы {А} коэффициентов в памяти ПЭВМ;
2) возможность с их помощью решать СЛАУ большего порядка по сравнению с прямыми методами;
3) они более алгоритмичны и требуют меньшей квалификации программиста;
4) при их использовании проще реализовывать геометрические граничные условия.
К недостаткам итерационных методов решения СЛАУ можно отнести необходимость:
1) наличия нулевого приближения, определение которого в отдельных случаях может представлять самостоятельную задачу (хотя метод итераций сойдется к решению даже тогда, если в качестве нулевого приближения взять просто нули);
2) определения точности итераций е.
В то же время, в методе Гаусса даже на относительно небольшой сетке матрицу {А} коэффициентов системы уравнений, имеющую от нескольких десятков и сотен тысяч элементов и больше, необходимо хранить в оперативной памяти ПЭВМ. Обращение к внешним запоминающим устройствам может в некоторой степени устранить этот недостаток, делая при этом процесс реализации задачи более длительным во времени и все более зависящим от технического состояния ПЭВМ и уровня математического обеспечения. Метод Гаусса требует более высокой квалификации программиста, так как вычислительный алгоритм в этом случае сложнее, нежели в итерационных методах. Следует отметить большую чувствительность метода Гаусса к способу аппроксимации дифференциальных соотношений разностными. К негативной стороне применения метода Гаусса можно отнести также и некоторую не-алгоритмичность составления матрицы в случае изменения граничных условий, формы и пр. Итерационные методы лишены этого недостатка, ибо в этом случае узел конечно-разностной сетки, имеющий отношение, например, к граничным условиям, просто не варьируется, т.е. он обходится во время итерационного процесса реализации СЛАУ.
Несколько слов о сходимости метода итераций. П.М. Сосис в работе [8], ссылаясь на [4,10,11], пишет: "... Было доказано, что существуют универсаль-
ные итерационные алгоритмы, обладающие неизбежной сходимостью. Например, процесс Зейделя обладает неизбежной сходимостью для симметричных положительно определенных матриц, а этим свойством обладают матрицы методов сил и деформаций (при правильной постановке задачи). Однако необходимо еще, чтобы в системе было диагональное преобладание (т.е. главные коэффициенты больше побочных). В противном случае процесс будет сходиться настолько медленно, что решение трудно будет получить".
Тем не менее, несмотря на перечисленные недостатки, метод Гаусса, применительно к рассматриваемому классу двумерных упругопластических задач, является более эффективным прежде всего с точки зрения времени решения СЛАУ. Очевидно, что в задачах оптимального проектирования [12], где рассматриваются десятки и сотни проектов, метод Гаусса, как способ решения системы алгебраических уравнений в задаче определения упругопластичес-кого напряженно-деформированного состояния, являясь более быстродействующим, предпочтителен по сравнению с итерационными методами.
Что касается решения СЛАУ в трехмерных задачах теории упругости (и пластичности), то здесь, анализируя результаты численного сравнения методов решения СЛАУ, трудно отдать явное предпочтение какому-либо методу. Итерационные методы, технически проще реализуемые, с гораздо меньшими временными затратами позволяют проводить оценочные расчеты и практически за то же время, что и прямой метод Гаусса, решать задачу теории упругости с достаточно высокой точностью. Кроме того, итерационные методы позволяют получать решения на значительно больших конечно-разно-стных сетках.
5. Выводы
Для решения задачи определения статического (и квазистатического) упругопластического напряженно-деформированного состояния тела при действии внешних объемных сил и поверхностных нагрузок в качестве основы взят вариационный принцип Лагранжа. Применение метода переменных параметров упругости сводит решение физически нелинейной задачи к решению последовательности упругих задач с уточняемыми на каждом шаге параметрами нелинейности. Реализация каждой отдельной упругой задачи проводится вариационно-разностным методом. Аппроксимация производных при дискретизации задачи в двумерных задачах проводится с помощью известных соотношений Ноха. Для трехмерной задачи используются выражения для пространственных производных, которые в трехмерном пространстве определяются через интеграл по замкнутой поверхности. Использование необходимого условия экстремума функции многих переменных (полной потенциальной энергии системы) сводит задачу определения НДС к решению системы линейных алгебраических уравнений отно-
системы линейных алгебраических уравнений относительно искомых перемещений в узлах конечно-разностной сетки. Получающаяся матрица коэффициентов СЛАУ симметрична, положительно определена и имеет ленточную структуру, что значительно облегчает поиск решения.
В заключение можно сделать некоторые выводы по методам решения системы линейных алгебраических уравнений большого порядка в задачах теории упругости (и пластичности) - наиболее важном и сложном этапе реализации ВРМ. Проведенный численный анализ методов реализации СЛАУ показал преимущества прямого метода (метод Гаусса) перед итерационными при решении двумерных уп-ругопластических задач. При реализации трехмерных задач однозначно нельзя выделить тот, или иной метод. Для проведения оценочных расчетов пред-
почтение следует отдать итерационным методам, а для определения НДС при примерно одинаковых временных затратах и одинаковых по точности результатах итерационные методы позволяют реали-зовывать решение задач на значительно больших ко-нечно-разностных сетках. В то же время использование обоих методов может служить одним из способов проверки достоверности полученных результатов.
Современные компьютеры позволяют решать упругопластическую задачу на сетках порядка (300x300). С помощью созданных методик удалось получить физически достоверные результаты для ряда задач на относительно грубых сетках. Этим определяется возможность применения разработанных методик для решения задач оптимального проектирования машиностроительных конструкций.
СПИСОКЛИТЕРАТУРЫ
1. Вазов В., Форсайт Дж. Разностные методы решения дифференциальных уравнений в частных производных. - М.: Иностранная литература, 1963. - 488 с.
2. Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978. - 592 с.
3. Копченова Н.В., Марон И.А. Вычислительная математика в примерах и задачах. - М.: Наука, 1972. -367 с.
4. Березин И.С., Жидков Н.П. Методы вычислений. -Т. 2. - М.: Физматгиз, 1959. - 620 с.
5. Демидович Б.П., Марон И.А. Основы вычислительной математики. - М.: Наука, 1970. - 664 с.
6. Калиткин H.H. Численные методы. - М.: Наука, 1978.- 512 с.
7. Барашков В.Н., Люкшин Б.А. К реализации вариа-ционно-разностного метода для осесимметричных
задач теории упругости и пластичности / Томск, ун-т. - Томск, 1983. - 14 с. - Деп. в ВИНИТИ 13.03.83, № 1335-83.
8. Сосис П.М. Статически неопределимые системы. -Киев: Буд1вельник, 1968. - 311 с.
9. Сосис П.М. Алгоритмический язык АЛГОЛ-бО и его применение в строительной механике. - Киев: Бущвельник, 1965. - 172 с.
10. Березин И.С., Жидков Н.П. Методы вычислений. -Т. 1. - М.: Физматгиз, 1959. - 464 с.
11. Фаддеев Д.К., Фадцеева В.Н. Вычислительные методы линейной алгебры. - М.: Физматгиз, 1963. - 734 с.
12. Барашков В.Н., Люкшин Б.А. Алгоритм прочностного проектирования осесимметричных упругопла-стических конструкций с использованием вариационно-разностного метода // Прикладные проблемы прочности и пластичности. Анализ и оптимизация деформируемых систем: Всесоюз. межвуз. сб. - Горький, 1988. - С. 91-97.
