Алгоритм компактного хранения и решения СЛАУ высокого порядка (стр. 3 из 5)
2 МЕТОДЫ КОМПАКТНОГО ХРАНЕНИЯ МАТРИЦЫ ЖЕСТКОСТИ
Матрица жесткости, получающаяся при применении МКЭ, обладает симметричной структурой, что позволяет в общем случае хранить только верхнюю треугольную часть матрицы. Однако для задач с большим количеством неизвестных это так же приводит к проблеме нехватки памяти. Предлагаемый в данной работе метод, позволяет хранить только ненулевые члены матрицы жесткости. Суть его заключается в следующем.
Первоначально, с целью выявления связей каждого узла с другими, производится анализ структуры дискретизации области на КЭ. Например, для КЭ - сетки, изображенной на рис. 1, соответствующая структура связей будет иметь вид:
| № узла | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Связи | 1, 2, 5, 6, 7 | 1, 2, 3, 6 | 2, 3, 4, 6 | 3, 4, 5, 6, 7 | 1, 4, 5, 7 | 1, 2, 3, 4, 6, 7 | 1, 4, 5, 6, 7 |
 |
Тогда, для хранения матрицы жесткости необходимо построчно запоминать информацию о коэффициентах, соответствующих узлам, с которыми связан данный узел. На рис. 2 приведены матрица жесткости и ее компактное представление для сетки изображенной на рис 1 [9].
 |
Текст подпрограммы, реализующий предложенный алгоритм анализа структуры КЭ-разбиения тела, приведен в Приложении 1.
Данный способ компактного хранения матрицы жесткости позволяет легко его использовать совместно с каким-нибудь численным методом. Наиболее удобным для этой цели представляется использование вышеизложенного итерационного метода Ланцоша, так как на каждой итерации требуется только перемножать матрицу коэффициентов СЛАУ и заданный вектор. Следовательно, для использования предложенного метода компактного хранения СЛАУ необходимо построить прямое и обратное преобразование в первоначальную квадратную матрицу.
Пусть
– элементпервоначальнойквадратной матрицы размерностью 
, а 
- ее компактное представление. Тогда для обратного преобразования будут справедливы следующие соотношения: 
, (*)где m – количество степеней свободы (m=1,2,3).
Для прямого преобразования будут справедливы соотношения, обратные к соотношениям (*).
3 ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ
Для проверки предлагаемого метода компактного хранения матрицы жесткости была решена задача о контактном взаимодействии оболочечной конструкции и ложемента [12] (рис. 4).
 |
Данная задача часто возникает на практике при транспортировке или хранении с горизонтальным расположением оси оболочечные конструкции устанавливаются на круговые опоры - ложементы. Взаимодействие подкрепленных оболочечных конструкций и ложементов осуществляется через опорные шпангоуты, протяженность которых вдоль оси оболочки соизмерима с шириной ложементов и много меньше радиуса оболочки и величины зоны контакта.
Данная задача решалась методом конечных элементов при помощи системы FORL [5]. Дискретная модель ложемента (в трехмерной постановке) представлена на Рис. 5.
 |
При построении данной КЭ-модели было использовано 880 узлов и 2016 КЭ в форме тетраэдра. Полный размер матрицы жесткости для такой задачи составляет
байт, что приблизительно равно 2,7 Мбайт оперативной памяти. Размер упакованного представления составил около 315 Кбайт.Данная задача решалась на ЭВМ с процессором Pentium 166 и 32 МБ ОЗУ двумя способами – методом Гаусса и методом Ланцоша. Сопоставление результатов решения приведено в Таблице 1.
Таблица 1.
| Время решения (сек) |  |  |  |  |  |  | |
| Метод Гаусса | 280 | 2.2101 | -2.4608 | 1.3756 | -5.2501 | 1.7406 | -2.3489 |
| Метод Ланцоша | 150 | 2.2137 | -2.4669 | 1.3904 | -5.2572 | 1.7433 | -2.3883 |
Из Таблицы 1 легко видеть, что результаты решения СЛАУ методом Гаусса и методом Ланцоша хорошо согласуются между собой, при этом время решения вторым способом почти в два раза меньше, чем в случае использования метода Гаусса.
ВЫВОДЫ.
В данной работе были рассмотрены способы компактного хранения матрицы коэффициентов системы линейных алгебраических уравнений (СЛАУ) и методы ее решения. Разработан алгоритм компактного хранения матрицы жесткости, позволяющий в несколько раз (иногда более чем в десятки раз) сократить объем требуемой памяти для хранения такой матрицы жесткости.
Классические методы хранения, учитывающие симметричную и ленточную структуру матриц жесткости, возникающих при применении метода конечных элементов (МКЭ), как правило, не применимы при решении контактных задач, так как при их решении матрицы жесткости нескольких тел объединяются в одну общую матрицу, ширина ленты которой может стремиться к порядку системы.
Предложенная в работе методика компактного хранения матриц коэффициентов СЛАУ и использования метода Ланцоша позволили на примере решения контактных задач добиться существенной экономии процессорного времени и затрат оперативной памяти.
СПИСОК ССЫЛОК.
1. Зенкевич О., Морган К. Конечные методы и аппроксимация // М.: Мир, 1980
2. Зенкевич О., Метод конечных элементов // М.: Мир., 1975
3. Стрэнг Г., Фикс Дж. Теория метода конечных элементов // М.: Мир, 1977
4. Бахвалов Н.С.,Жидков Н.П., Кобельков Г.М. Численные методы // М.: наука, 1987
5. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления // М.:Наука, 1984
6. Бахвалов Н.С. Численные методы // М.: Наука, 1975
7. Годунов С.К. Решение систем линейных уравнений // Новосибирск: Наука, 1980
8. Гоменюк С.И., Толок В.А. Инструментальная система анализа задач механики деформируемого твердого тела // Приднiпровський науковий вiсник – 1997. – №4.
9. F.G. Gustavson, “Some basic techniques for solving sparse matrix algorithms”, // editer by D.J. Rose and R.A.Willoughby, Plenum Press, New York, 1972
10. А.Джордж, Дж. Лиу, Численное решение больших разреженных систем уравнений // Москва, Мир, 1984
11. D.J. Rose, “A graph theoretic study of the numerical solution of sparse positive definite system of linear equations” // New York, Academic Press, 1972
12. Мосаковский В.И., Гудрамович В.С., Макеев Е.М., Контактные задачи теории оболочек и стержней // М.:”Машиностроение”, 1978
ПРИЛОЖЕНИЕ 1
Исходный текст программы, реализующий анализ структуры КЭ-разбиения объекта.
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <fstream.h>
#include "matrix.h"
#define BASE3D_4 4
#define BASE3D_8 8
#define BASE3D_10 10
const double Eps = 1.0E-10;
DWORD CurrentType = BASE3D_10;
void PrintHeader(void)
{
printf("Command format: AConvert -<switch> <file1.in> <file2.in> <file3.out> [/Options]\n");
printf("Switch: -t10 - Tetraedr(10)\n");
printf(" -c8 - Cube(8)\n");
printf(" -s4 - Shell(4)\n");
printf(" -s8 - Shell(8)\n\n");
printf("Optins: /8 - convert Tetraedr(10)->8*Tetraedr(4)\n");
printf(" /6 - convert Cube(8)->6*Tetraedr(4)\n");
}
bool Output(char* fname,Vector<double>& x,Vector<double>& y,Vector<double>& z, Matrix<DWORD>& tr, DWORD n,
DWORD NumNewPoints,DWORD ntr,Matrix<DWORD>& Bounds,DWORD CountBn)
{
char* Label = "NTRout";
int type = CurrentType,
type1 = (type==BASE3D_4 || type==BASE3D_10) ? 3 : 4;
DWORD NewSize,
i,
j;
ofstream Out;
if (type==BASE3D_4) type1 = 3;
else if (type==BASE3D_8) type1 = 4;
else if (type==BASE3D_10) type1 = 6;
Out.open(fname,ios::out | ios:: binary);
if (Out.bad()) return true;
Out.write((const char*)Label,6 * sizeof(char));
if (Out.fail()) return true;
Out.write((const char*)&type,sizeof(int));
if (Out.fail()) return true;
Out.write((const char*)&CountBn,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
Out.write((const char*)&(NewSize = n + NumNewPoints),sizeof(DWORD));
if (Out.fail()) return true;
Out.write((const char*)&(NumNewPoints),sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
for (DWORD i = 0; i < n; i++)
{
Out.write((const char*)&x[i],sizeof(double));
Out.write((const char*)&y[i],sizeof(double));
Out.write((const char*)&z[i],sizeof(double));
if (Out.fail())
{
Out.close();
return true;
}
}
for (i = 0; i < NumNewPoints; i++)
{
Out.write((const char*)&x[n + i],sizeof(double));
Out.write((const char*)&y[n + i],sizeof(double));
if (Out.fail())
{
Out.close();
return true;
}
}
Out.write((const char*)&(ntr),sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
for (i = 0; i < ntr; i++)
for (j = 0; j < (DWORD)type; j++)
{
DWORD out = tr[i][j];
Out.write((const char*)&out,sizeof(DWORD));
if (Out.fail())
{
Out.close();
return true;
}
}
for (i = 0; i < CountBn; i++)
for (j = 0; j < (DWORD)type1; j++)
{
DWORD out = Bounds[i][j];
Out.write((const char*)&out,sizeof(DWORD));
if (Out.fail())