CC BY
10Вестник Сыктывкарского университета. Сер. 1.Вып. 4.2001
УДК 539.3
Разреженные матрицы при решении задач теории оболочек
В. Л. Никит, енков
Рассмотрена структура разреженных матриц большой раз-мерйости, получаемых при решении линейного приближения для краевых задач теории оболочек двумерным методом сплайн-коллокации. Обсуждаются вопросы упаковки-распаковки и определения программной среды при реализации метода Гаусса с выбором главного элемента..
При решении нелинейных краевых задач теории оболочек методом сплайн-коллокации[1, 2] появляется необходимость реализации итерационной процедуры для операторного уравнения второго рода, составной частью которой является СЛАУ с сильно разреженной матрицей. Размерность матрицы такова, что не позволяет для ее хранения использовать оперативную и даже дисковую память. В данной работе проводится анализ структуры матрицы СЛАУ, предложен простой способ упаковки-распаковки и программная реализация метода исключения Гаусса, построенная на его основе х.
1. Нелинейный вариант разрешающих уравнений почти цилиндрических оболочек
Параметризация с помощью круговой цилиндической оболочки а0 радиуса Е0, позволяет выписать основные геометические соотношения для произвольной почти цилиндической оболочки а [2, 3, 4, 5].
Пусть поверхность отсчета сг0 отображается на искомую поверхность а с помощью векторного равенства (см. рис. 1)
г(а!,а2) = г(аьо2) + Н(аиа2)п (1)
*В работе частично использованы материалы дипломных работ выпускниц математического факультета 1999 года В.И. Казаниной и Ю.В. Басацкой.
функция Н(а-[,а-2) является параметром
Рис.1
Введем обозначения
#1
дн
да1'
Н2 =
дН
дао
(2)
Тогда геометрические параметры исходной поверхности (компоненты метрического тензора аг) и тензора кривизны Ь^, символы Кристоф-феля Гкг] и др.) будут функциями от Я, #ь Я2[4, 5].
Будем предполагать, что криволинейные координаты на исходной поверхности а ортогональны, т.е. (т.е. либо #1 = 0 , либо Я2 = 0 ). При #1 = 0 имеем цилиндрическую поверхность с произвольной формой поперечного сечения. При Я2 = О имеем круговую цилиндрическую оболочку переменного радиуса (оболочка вращения).
Далее, следуя [6, 7], можно получить нелинейную (в соответствии с принятой квадратичной аппроксимацией) систему разрешающих уравнений равновесия оболочки относительно трех перемещений гц, и2, в матричной форме записи для указанных оболочек
¿(о1,а2,Я,Я,-)и + С(«1, а2, Я, Я,, ги) = я,
(3)
Разреженные матрицы при решении задач теории оболочек... 187 где
Ь = |||| . ¿,7 = 1:3- матрица из операторов дифференцирования соответствующих линейной теории оболочек; в (го) = [6*1 («>), О2{и'), 6'„('ш)] - вектор нелинейных слагаемых; Ч — [91- Я2' Чп] ~ вектор приведенных нагрузок; и = [г<], к»]- вектор смещений.
2. Численный метод
В рамках итерационного подхода [8] к системе (3) одним из этапов является решение линейной системы
Цаиа2,Н,Нг) и<*> = Ч<*>. (4)
Условие периодичности решения (4) по окружной переменной а 2 будет выполнено, если перемещения в узлах некоторой сетки ((сц\с*2к)а € (0. /V — [),к Е (0,М)) искать в виде двумерных В-сп л айнов [9, 10]
А'-1 М+2 N-1 М+2
щ ~ Е Е хыВк(а12)Ьг{а{) = ^ ^ хкгВ0(а2 а2)Ь0(а\ - о^),
к=03=-2 к=0з=-2
/V—1 М+2 Ы-1 М+2
= Е Е ^гВк(а12)Ь;(а\) = Е Е йиВоМ - ак2)Ъ0(а\ - а'),
А,—0 ] = -2 к=0 з = -2
Л' —1 М+2 ЛГ-1 М+2
(С = Е Е = Е Е хьгВ0(а12 - а£)Ьо(а^ - "!),
/г=0 2 ¿=0 ¿=-2
(5)
где
Дз(а) - 27г-периодический В-сплайн,
Ь0(а) - В-сплайн, удовлетворяющий граничным условиям по
После подстановки (5) в (4) система, примет вид (индекс к опускаем)
' ¿и (Е^о1 хкгВк(а'2Ща\)) + Ь12 ¿Ма^ьЩ +
+¿13 = д1{а{ла12),
¿л (е^ГО Й-2.г,,Л,(П'2)/,;(П;)) + ¿и (е^1 +
+¿23 (Е^о1 Ей-" )) = ^(п(,о'2).
¿3. (Е^Го1 Е; = -2 -а, «,(<>',)/,,(О0) + ¿32 (Е^О1 Е^-2 +
+¿13 (Е^о1 = </"(а{,а<).
(6)
Добавляя в систему необходимое количество граничных условий (для доопределения сплайна /)0(«1) за пределами сетки воспользуемся однородными условиями на перемещения), приходим к СЛАУ
А[п.п\ ■ А"[п] (7)
относительно п = ЗАГ(Л/ + 5) переменных Х[п] = (||-1'и|!рь'||), расположенных в указанном порядке.
Структура матрицы А приведена на рис.2, где выделены области расположения ненулевых элементов.
первый регулярный блок_
О
I сдвиг на ЗК элементов
О
3. Структура матрицы и степень разреженности
Данная матрица в области ненулевых элементов состоит из прямоугольных блоков размером ЗА ■ 1ЬN элементов. Два первых и два последних блока составлены из бЛг строк, порожденных граничными условиями, и Г2Л; строк, относящихся к доопределению сплайна Ь0(аг). Остальные блоки включают в себя по ЗДГ строк, соответствующих точкам .сетки по окружной координате ог2 для трех уравнений разрешающей системы (3). В состав прямоугольного блока входят пять прямоугольных блоков ЗА' • ЗА\ каждый из которых образован представителями групп неизвестных (Ц^ь'Ц) в каждом из трех уравнений системы (см. рис. 3). И, наконец, прямоугольный блок N ■ N (рис. 4) соответствует одной из групп неизвестных (например, (||.гь||) в одном из уравнений системы. Этот блок (назовем его атомарным) имеет структуру почти пятидиагоналыюй матрицы, нарушаемую (в основаниях побочной диагонали) за счет условий периодичности используемых /?-сплайнов.
^71 1''
' г
---[— |
4 Н-.....
- ь„ ц, г
- - и,
-
I - и. и. Ци
-н- -
Рис.3
Подсчитаем теперь число ненулевых элементов матрицы. Атомарный блок имеет в каждой строке пять ненулевых элементов. В прямоугольный блок входит пятнадцать атомарных. Следовательно, в каждой строке матрицы А[п,п] ровно 75 ненулевых из 3N( M + 5) элементов. Степень заполненности матрицы (в процентах) равна,
Sn = 75/7* • 100% (8)
и убывает пропорционально п. 'Гак 5150 ~ 16.6% (N = M = 10), •Siобоо = 0.25% {N = M = 100), a при (Л' = M - 1000) ,S'10o5ooo = 0.0025% (один из 40000 элементов - ненулевой). Если учесть число элементов матрицы хотя бы при N — M = 100 (99225 • 104), то объем памяти для ее хранения составит примерно 15 Гб. Это обстоятельство не позволяет использовать обычные алгоритмы решения СЛАУ, а требует реализации методов, работающих с упакованными матрицами.
4. Упаковка и реализация метода Гаусса
Регулярная структура атомарного блока позволяет произвести плотную упаковку исходной матрицы без использования дополнительных разметочных массивов и получить упакованную матрицу размером п-75 (в рассматриваемом примере она содержит 2362500 элементов и требует для хранения примерно 37 Мб памяти). Ниже приведены фрагменты Pascal-процедур, реализующих упаковку (Pack) и распаковку (Unpack) строки.
{процедура упаковки строки b ---> у,
1 - номер строки в атомарном блоке) Procedure Pk(l : integer; var b : Vector; var y : Pvector ); var i,j : int eger ;
for i:»l to 15 do
for j:-l to 5 do y[5*(i-l)+j]:=b[n*(i-l)+l-3+j];
У[76]:=b[mr];{правая часть} у[0]:=1;{номер строки} end{of Pk};
{процедура распаковки строки у ---> Ь}
Procedure Upk(l : integer; var b : Vector; var у : Pvector );
var i,j : integer;
begin
for i: = l to nr do b[i]:=0; 1:=Round(y [0]);
else for i:=l to 15 do
for j:=1 to 5 do b[n*(i-l)+l-3+j]:=y[5*(i-l)+j];
b [nr] : =y [76] ; end{of Upk};
Метод Гаусса с выбором главного элемента использует данные процедуры упаковки и распаковки внутри своих итераций следующим образом:
1. Выбор главного элемента
(a) Вызов очередной строки (Unpack)
(b) Выбор максимального по абсолютной величине ведущего элемента
2. Итерация исключения
(a) Вызов ведущей строки (Unpack)
(b) Преобразование ведущей строки (Pack)
(c) Вызов текущей строки (Unpack)
(d) Преобразование текущей строки (Pack)
3. Обратный ход
(a) Вызов текущей строки (Unpack)
(b) Вычисление очередной компоненты решения (Pack)
Процедуры исключения метода Гаусса будут различаться для следующих частей матрицы А (см. рис. 2 и фрагмент программы ниже).
{процедура исключения для всей матрицы (прямой ход)}
Procedure GAUSS.D;
begin
BoardGl; {первый блок граничных условий) BoardG2; {второй блок граничных условий} for nb:=0 to m-3 do RegG(nb); {регулярные блоки} NonRegGl;
NonRegG2; {предпоследнийй блок граничных условий} BoardGm; {последний блок граничных условий} end{of GAUSS.D};
Процедура обратного хода пишется для прямоугольного блока и имеет особенности применения для регулярных и нерегулярных блоков.
{процедура обратного хода для прямоугольного блока}
Procédure GR(iO : integer ; ni : integer);
begin
i:=n0; repeat
As(ni,y,A);{строка матрицы A в y} Upk(с,y);{упаковка y в с} s : =0 ;
for j:=i0+l to nr-1 do s:=s+c[j]*yl[nl+j-iO]; ylCnl] :=c[nr]-s; ni :=nl-l; iO:=i0-l; i:=i-l; until i<0; end{of GR};
Pascal-nporpaMMa метода исключения Гаусса для упакованной матрицы была реализована в Turbo-pascal 7.0 и протестирована на примерах с матрицами подобной структуры и известными решениями.
5. Выбор среды программной реализации
При расчете конструкций в рамках полученной модели (например, овальных корпусов аппаратов давления или корпусов с локальными вмятинами и выпучиваниями[11, 12, 13]) требуется довольно значительное число узлов сетки по каждой из координат (TV, M ~ 5 • 102 — 103). Так, при прочностном анализе подопорной области автоклава АП16 — 2 х 40.4 (длина корпуса 40.4 м., диаметр 2 м.) необходимо, чтобы по осевой координате на рассматриваемую область приходилось хотя бы десять точек сетки. Осевой размер подопорной области - 0.5 м. Тогда сетка по координате «i должна содержать не менее 20 х 40.4 = 808 точек. По окружной координате наиболее интересным в прочностном отношении является консольный участок опоры протяженностью примерно 3° (или ~ 0.05 м.). Это приводит к еще большему числу точек сетки по координате а2. Сетка должна быть еще мельче, если анализируется поведение параметров НДС (напряженно-деформированного состояния) вблизи локальных вмятин или мест выпучивания корпуса. Это приводит к СЛАУ с матрицами, содержащими миллионы переменных и требующими для хранения (в упакованном виде) памяти порядка от сотен мегабайт до десятков гигабайт. Ясно, что при выборе программной среды следует остановиться на тех средствах создания приложений для Windows, в которых достаточно элементарно осуществляется доступ к динамическому выделению необходимой дисковой
памяти (лучше, если о таком доступе заботиться вообще не надо). Такими возможностями обладает Delphi, тем более что Pascal-процедуры встраиваются туда практически без изменений.
Литература
1. Никитенков B.JT. Деформационный вариант нелинейных уравнений статики ребристых оболочек. -"Автоклавы. Расчет, проектио-вание, опыт эксплуатации//!1. Всесоюзного семинара "Автоматизация инженерных расчетов при проектиовании автоклавов" (под ред. проф. Е.И.Михайловского), Сыктывкар, 1992.С. 168-197.
2. Развитие численных методов нелинейной теории оболочек// Отчет о НИР (рук. B.JI. Никитенков). Сыктывкар, 1998. 114 с.
3. Корнишин М.С., Паймушин В.Н., Снигирев В.Ф. Вычислительная геометрия в задачах механики оболочек. М: Наука, 1989. 207 с.
4. Казанина. В.И. Нелинейные уравнения для овальной цилиндрической оболочки// Дипломная работа (рук В Л. Никитенков). Сык-тывкар:СыктГУ, 1999. 36 с.
5. Басацкая Ю.В. Нелинейные уравнения почти цилиндрической оболочки вращения// Дипломная работа (рук B.JI. Никитенков). Сыктывкар:СыктГУ, 1999. 41 с.
6. Черных К.Ф. Нелинейная теория упругости в машиностроительных расчетах. - Л.: Машиностроение, Денингр. отд-ние, 1986.336с.
7. Новожилов В.В., Черных К.Ф., Михайловский Е.И. Линейная теория тонких оболочек. - Л.: Политехника, 1991.656с.
8. Nikitenkov V. L. Many-layer iterative procedures for solving the problems of mechanics of deformed rigid body//Transactions of SPASP. V. 1. 1997. P.101-115.
9. Малоземов B.H., Певный A.B. Полиномиальные сплайны. Л.:Изд-во Ленингр. ун-та, 1986. 120 с.
10. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.:Наука, 1980. 352 с.
11. Никитенков B.JI. Вопросы прочности и проектирования тяжелых горизонтальных аппаратов давления//Автореф. дисс. докт. техн. наук. С-Пб.:1996. 39 с.
12. Доценко В.Д. Исследование влияния овальности на. напряженно-деформированное состояние автоклавов./ Автореф. дисс. канд. техн. наук, Гатчина, 1986, 14 с.
13. Никитенков B.JI., Белых Д.Г. Система АВГОР// Автоматизация проектиования, 1997. 4. С.46-50.
Summary
Nikitenkov V.L. Rarefied matrixes in problems of shell theory
In this article the sructure of rarefied matrixes is considered. These matrixes received as a result of solution of linear boundary value problems of shell theory by two-dimensional spline-collocation. Also questions of packing and program realization of Gauss method are discussed.
Сыктывкарский университет
Поступила 21.10.2000
