49829 (Суперэлементное моделирование пространственной системы "плита – грунтовое основание"), страница 3

Применение МКЭ для решения системы “плита-основание” приводит к системе линейных алгебраических уравнений с ленточной симметричной матрицей. Ширина ее полуленты зависит от порядка нумерации узлов и определяется по формуле: B= (R+1) Q, где R - максимальная разность разностей номеров узлов конечных элементов, Q - число неизвестных (степеней свободы) в каждом узле.

6. Алгоритмы построения и решения дискретной модели

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

Второй этап- дискретизация расчетной области, учитывающая особенности структуры грунтового основания. В результате строится нерегулярная решетка с массивами шагов по координатным осям. Каждый параллелепипед дискретной решетки делится на шесть тетраэдральных элементов.

Для каждого конечного элемента (тетраэдра) необходимо задать характеристики: модуль упругости, коэффициент Пуассона.

Третий этап - задание краевых условий. Граничные условия расчетной области определяются системой внешних сил и выбором размеров расчетной области (этап 1). Система внешних сил задается в виде вектора нагрузок, определенного для всех узлов расчетной области. С каждым узлом связано три значения нагрузки: одно по направлению оси OX, второе по направлению оси OY, третье по направлению оси OZ. Вектор нагрузок задается на верхнем ребре. На всех остальных обычно задаются нулевые перемещения. Четвертый этап - формирование матрицы жесткости. Построение матрицы жесткости производится с учетом ее особенностей: симметричности, ленточности. Матрица жесткости (МЖ) размещается в ОП упакованной в прямоугольник, т.е. хранится верхняя полулента. Для построения МЖ используется аналитический алгоритм построения [1].

Согласно которому матрица жесткости имеет вид:

где

где i - номер узла, связанного с узлами j; j=1,2,3,4;

Пятый этап - учет граничных условий в МЖ. Используется вектор усилий и вектор корректировки, с помощью которого описываются задаваемые граничные значения перемещений. Учёт граничных условий приводит к изменению матрицы жёсткости [K] и векторов узловых сил и перемещений. Матрица [K] уже не будет сингулярной.

Шестой этап - решение системы линейных алгебраических уравнений. На этом этапе используется метод квадратного корня, учитывающий упаковку МЖ в прямоугольник.

Этот метод состоит в следующем:

Если матрица симметрическая, то её можно представить следующим образом:

A=S*DS,

Где S - верхняя треугольная матрица с положительными элементами на главной диагонали; D - диагональная матрица, с элементами +1 или -1 на главной диагонали; S* - нижняя треугольная матрица. Коэффициенты и вычисляются по формулам:

i=j то, ;

;

i

;

В том случае, если матрица A самосопряжённая и положительно определённая, то матрицу D можно опустить, так как она будет единичной. Метод осуществляется по следующей схеме:

сначала решаем уравнение S*Y=B

затем уравнение SX=Y, находя решение системы.

Наша работа заключается в решении СЛАУ методом квадратного корня, используя ленточную симметрическую матрицу, компактно упакованную.

Полуленточная матрица системы строиться следующим образом:

В методе квадратного корня используется функция, с помощью которой меняются оба индекса.

7. Описание и инструкция работы с приложением

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

Пользователю имеет возможность задавать характеристики каждого отдельного симплекс-элемента (тетраэдра). Благодаря этому система “плита-грунтовое основание" может быть рассмотрена как неоднородная.

Ввод вектора узловых сил осущестляется путем ввода величины силы узлу к которому она прилагается. Вектор перемещений сразу предполагает перемещения по трем осям. Однако пользователь имеет возможность запретить перемещения по какой-либо оси для каждого узла.

Выходными характеристиками приложения являются перемещениями в необходимых узлах. По анализу которых делается вывод об осадке плиты.

8. Верификация приложения

Рассмотрим следующий пример:

Однородная плита располагается вертикально на жёстком основание. Усилие Р равномерно распределено по верхнему основанию плиты. Дискретизация пластины производится путём разбиения ее на конечные элементы (тетраэдры). Узлы и полученные конечные элементы нумеруются. Программное приложение рассчитывает значения перемещений в каждом узле модели.

При введении в качестве параметров модели тестового примера следующих величин: нагрузка на плиту Р = 100кг,

параметры плиты: модуль упругости Е = 360 кг/см², коэффициент Пуассона =0.2, h=100см, l =100см, приращения по Ox,Oy, Oz =50см.

Разрешены только вертикальные перемещения.

Полученный результат имеет вид:

вертикальные перемещения в узлах:

1.09651.03951.06241.0996

0.53560.53080.53420.5345

0.00000.00000.00000.0000

Заключение

В курсовой работе реализовано моделирование расчета осадок большеразмерной плиты на основании системного подхода и метода конечных элементов. Разработаны алгоритмы построения матрицы жесткости упакованной в прямоугольник, решения системы линейных алгебраических уравнений для упакованной МЖ. Спроектирован удобный интерфейс ввода исходных данных и вывода результатов. Создан программный продукт моделирования расчета осадок в среде Delphi 5.0.

Проведена верификация программного продукта на основе задачи, имеющей аналитическое решение. Она показала хорошее совпадение результатов с точностью примерно 90-95%.

Разработанное приложение может быть использовано для предварительных расчетов оснований фундаментов плит с учетом сложной структуры основания в инженерной практике на этапе проектирования.

Список использованных источников

1. Быховцев В.Е. Компактный алгоритм построения матрицы жесткости в МКЭ. - Известия АН БССР, серия физ. - матем. наук, №1, 1983, с.34-37.

2. Быховцев В.Е., Ермашов В.П., Богданова Т.Г. Влияние формы фундамента на его осадки. - Фундаменты на искусственных основаниях в условиях Белорусской ССР, сб. научных трудов, Минск: БелНИИС, 1986, с.47-55.

3. Винокуров Е.Ф. Расчёт оснований и фундаментов. - Минск: АН БССР, 1960. - 294 с.

4. Галлагер Р. Метод конечных элементов: основы. - М.: "Мир", 1984. - 428с.

5. Ильюшин А.А. Механика сплошной среды: Учебник. - 3-е изд. - Москва: Издательство МГУ, 1990. - 310 с.

6. Зенкевич О. Метод конечных элементов в технике: Перев. с англ. - Москва: "Мир", 1975. - 544с.

7. Лурье А.И. Нелинейная теория упругости. - Москва: "Наука". Главная редакция физико-математической литературы, 1980. - 512 с.

8. Партон В.З. Перлин П.И. Методы математической теории упругости: Учебное пособие. - Москва: Наука. Главная редакция физико-математической литературы, 1981. - 688 с.

9. Сесков В.Е., Быховцев В.Е., Лях В.Н., Цурганова Л.А. Определение несущей способности и осадки микросвайных фундаментов в выштампованных скважинах методом вычислительного и физического экспериментов. - Основания и фундаменты, сб. н. трудов, Минск: НПТО "Белстройнаука", 1986, с.26-35.

10. Сивцова Е.П. Расчёт осадки одиночной сваи с учётом работы острия. - Сб. трудов НИИ оснований №53, М., 1963.

11. Цытович H.А. Механика грунтов. - Москва: Госстройиздат, 1963. - 636 с.

Приложение

(основные функции, процедуры и алгоритмы приложения)

функции подсчета коэффициентов b, c, d используемых для формирования матрицы жесткости function det (a11,a12,a13,a21,a22,a23,a31,a32,a33: integer): real;

begin

det: =a11*a22*a33-a11*a23*a32+a12*a23*a31-

a12*a21*a33+a13*a21*a32-a13*a22*a31;

end;

function formb (a,ntetr: integer): real;

begin

if a=cells [ntetr,1] then formb: =-det (1,kordy [ntetr,2],kordz [ntetr,2],1,kordy [ntetr,3],kordz [ntetr,3],1,kordy [ntetr,4],kordz [ntetr,4]);

if a=cells [ntetr,2] then formb: =det (1,kordy [ntetr,3],kordz [ntetr,3],1,kordy [ntetr,4],kordz [ntetr,4],1,kordy [ntetr,1],kordz [ntetr,1]);

if a=cells [ntetr,3] then formb: =-det (1,kordy [ntetr,4],kordz [ntetr,4],1,kordy [ntetr,1],kordz [ntetr,1],1,kordy [ntetr,2],kordz [ntetr,2]);

if a=cells [ntetr,4] then formb: =det (1,kordy [ntetr,1],kordz [ntetr,1],1,kordy [ntetr,2],kordz [ntetr,2],1,kordy [ntetr,3],kordz [ntetr,3]);

end;

function formc (a,ntetr: integer): real;

begin

if a=cells [ntetr,1] then formc: =det (1,kordx [ntetr,2],kordz [ntetr,2],1,kordx [ntetr,3],kordz [ntetr,3],1,kordx [ntetr,4],kordz [ntetr,4]);

if a=cells [ntetr,2] then formc: =-det (1,kordx [ntetr,3],kordz [ntetr,3],1,kordx [ntetr,4],kordz [ntetr,4],1,kordx [ntetr,1],kordz [ntetr,1]);

if a=cells [ntetr,3] then formc: =det (1,kordx [ntetr,4],kordz [ntetr,4],1,kordx [ntetr,1],kordz [ntetr,1],1,kordx [ntetr,2],kordz [ntetr,2]);

if a=cells [ntetr,4] then formc: =-det (1,kordx [ntetr,1],kordz [ntetr,1],1,kordx [ntetr,2],kordz [ntetr,2],1,kordx [ntetr,3],kordz [ntetr,3]);

end;

function formd (a,ntetr: integer): real;

begin

if a=cells [ntetr,1] then formd: =-det (1,kordx [ntetr,2],kordy [ntetr,2],1,kordx [ntetr,3],kordy [ntetr,3],1,kordx [ntetr,4],kordy [ntetr,4]);

if a=cells [ntetr,2] then formd: =det (1,kordx [ntetr,3],kordy [ntetr,3],1,kordx [ntetr,4],kordy [ntetr,4],1,kordx [ntetr,1],kordy [ntetr,1]);

if a=cells [ntetr,3] then formd: =-det (1,kordx [ntetr,4],kordy [ntetr,4],1,kordx [ntetr,1],kordy [ntetr,1],1,kordx [ntetr,2],kordy [ntetr,2]);

if a=cells [ntetr,4] then formd: =det (1,kordx [ntetr,1],kordy [ntetr,1],1,kordx [ntetr,2],kordy [ntetr,2],1,kordx [ntetr,3],kordy [ntetr,3]);

end;

процедура формирования матрицы жесткости

procedure formprmatr (a,b,k: integer);

var ro,G,lya,Mu,E,vv: extended;

begin

Mu: =0.2; E: =360;

G: =E/ (2* (1+Mu));

lya: = (2*Mu*G) / (1-2*Mu);

ro: =2*G+lya;

vv: =1/ (360*V [ (k div 7) +1]);

prmatr [3*a-2,3*b-3*a+1]: =prmatr [3*a-2,3*b-3*a+1] +vv* (formb (a,k) *formb (b,k) *ro+G* (formc (a,k) *formc (b,k) +formd (a,k) *formd (b,k)));

prmatr [3*a-2,3*b-3*a+2]: =prmatr [3*a-2,3*b-3*a+2] +vv* (formb (a,k) *formc (b,k) *lya+formc (a,k) *formb (b,k) *G);

prmatr [3*a-2,3*b-3*a+3]: =prmatr [3*a-2,3*b-3*a+3] +vv* (formb (a,k) *formd (b,k) *lya+formd (a,k) *formb (b,k) *G);

if (3*a-1<=3*b-2) then prmatr [3*a-1,3*b-3*a]: =prmatr [3*a-1,3*b-3*a] +vv* (formc (a,k) *formb (b,k) *lya+formb (a,k) *formc (b,k) *G);

prmatr [3*a-1,3*b-3*a+1]: =prmatr [3*a-1,3*b-3*a+1] +vv* (formc (a,k) *formc (b,k) *ro+G* (formb (a,k) *formb (b,k) +formd (a,k) *formd (b,k)));

prmatr [3*a-1,3*b-3*a+2]: =prmatr [3*a-1,3*b-3*a+2] +vv* (formc (a,k) *formd (b,k) *lya+formd (a,k) *formc (b,k) *G);

if (3*a<=3*b-2) then prmatr [3*a,3*b-3*a-1]: =prmatr [3*a,3*b-3*a-1] +vv* (formd (a,k) *formb (b,k) *lya+formb (a,k) *formd (b,k) *G);

if (3*a<=3*b-1) then prmatr [3*a,3*b-3*a]: =prmatr [3*a,3*b-3*a] +vv* (formd (a,k) *formc (b,k) *lya+formc (a,k) *formd (b,k) *G);

prmatr [3*a,3*b-3*a+1]: =prmatr [3*a,3*b-3*a+1] +vv* (formd (a,k) *formd (b,k) *ro+G* (formc (a,k) *formc (b,k) +formb (a,k) *formb (b,k)));

end;

функция для получения необходимого элемента в прямоугольной матрице

function value (i: integer; j: integer): real;

begin

if (j=i+m) then value: =0;

if (i>j) then value: =prmatr [j, i-j+1] ;

if (i<=j) then value: =prmatr [i,j-i+1] ;

end;

функция для получения элементов прямоугольной матрицы с учетом ее сжатия

function value2 (i,j: integer): real;

var k,n1: integer;

begin

n1: =kx*ky*kz*3;

for k: =1 to n1 do

begin

if (P [k] =1) and (i>=k) then inc (i);

if (P [k] =1) and (j>=k) then inc (j);

end;

value2: =value (i,j);

end;

процедура сжатия силового вектора

procedure compressf (k: integer);

var i: integer;

begin

for i: =k to 3*kx*ky*kz do F [i]: =F [i+1] ;

inc (Count);

end;

функция, возвращающая значение сигнума от числа

function sign (f: real): shortint;

begin

if f=0 then sign: =0

else sign: =round (abs (f) /f);

end;

алгоритм занесения координат и номеров узлов тетраздров

point: =1;

kp: =ky*kz;

ntetr: =1;

x: =0;

for i: =1 to kx-1 do

begin

y: =0;

for j: =1 to ky-1 do

begin

z: =0;

for k: =1 to kz-1 do

begin

k1: =point+kp;

k2: =point+kz+1;

z: =z+hz [k] ;

cells [ntetr,1]: =k2; kordx [ntetr,1]: =x; kordy [ntetr,1]: =y+hy [j] ; kordz [ntetr,1]: =z;

cells [ntetr,2]: =point; kordx [ntetr,2]: =x; kordy [ntetr,2]: =y; kordz [ntetr,2]: =z-hz [k] ;

cells [ntetr,3]: =point+1; kordx [ntetr,3]: =x; kordy [ntetr,3]: =y; kordz [ntetr,3]: =z;

cells [ntetr,4]: =k1; kordx [ntetr,4]: =x+hx [i] ; kordy [ntetr,4]: =y; kordz [ntetr,4]: =z-hz [k] ;

cells [ntetr+1,1]: =k2; kordx [ntetr+1,1]: =x; kordy [ntetr+1,1]: =y+hy [j] ; kordz [ntetr+1,1]: =z;

cells [ntetr+1,2]: =point+1; kordx [ntetr+1,2]: =x; kordy [ntetr+1,2]: =y; kordz [ntetr+1,2]: =z;

cells [ntetr+1,3]: =k1+1; kordx [ntetr+1,3]: =x+hx [i] ; kordy [ntetr+1,3]: =y; kordz [ntetr+1,3]: =z;

cells [ntetr+1,4]: =k1; kordx [ntetr+1,4]: =x+hx [i] ; kordy [ntetr+1,4]: =y; kordz [ntetr+1,4]: =z-hz [k] ;

cells [ntetr+2,1]: =k2; kordx [ntetr+2,1]: =x; kordy [ntetr+2,1]: =y+hy [j] ; kordz [ntetr+2,1]: =z;

cells [ntetr+2,2]: =k1+1; kordx [ntetr+2,2]: =x+hx [i] ; kordy [ntetr+2,2]: =y; kordz [ntetr+2,2]: =z;

cells [ntetr+2,3]: =k1+kz+1; kordx [ntetr+2,3]: =x+hx [i] ; kordy [ntetr+2,3]: =y+hy [j] ; kordz [ntetr+2,3]: =z;

cells [ntetr+2,4]: =k1; kordx [ntetr+2,4]: =x+hx [i] ; kordy [ntetr+2,4]: =y; kordz [ntetr+2,4]: =z-hz [k] ;

cells [ntetr+3,1]: =k2; kordx [ntetr+3,1]: =x; kordy [ntetr+3,1]: =y+hy [j] ; kordz [ntetr+3,1]: =z;

cells [ntetr+3,2]: =k1+kz+1; kordx [ntetr+3,2]: =x+hx [i] ; kordy [ntetr+3,2]: =y+hy [j] ; kordz [ntetr+3,2]: =z;

cells [ntetr+3,3]: =k1+kz; kordx [ntetr+3,3]: =x+hx [i] ; kordy [ntetr+3,3]: =y+hy [j] ; kordz [ntetr+3,3]: =z-hz [k] ;

cells [ntetr+3,4]: =k1; kordx [ntetr+3,4]: =x+hx [i] ; kordy [ntetr+3,4]: =y; kordz [ntetr+3,4]: =z-hz [k] ;

cells [ntetr+4,1]: =k2; kordx [ntetr+4,1]: =x; kordy [ntetr+4,1]: =y+hy [j] ; kordz [ntetr+4,1]: =z;

cells [ntetr+4,2]: =k1+kz; kordx [ntetr+4,2]: =x+hx [i] ; kordy [ntetr+4,2]: =y+hy [j] ; kordz [ntetr+4,2]: =z-hz [k] ;

cells [ntetr+4,3]: =point+kz; kordx [ntetr+4,3]: =x; kordy [ntetr+4,3]: =y+hy [j] ; kordz [ntetr+4,3]: =z-hz [k] ;

cells [ntetr+4,4]: =k1; kordx [ntetr+4,4]: =x+hx [i] ; kordy [ntetr+4,4]: =y; kordz [ntetr+4,4]: =z-hz [k] ;

cells [ntetr+5,1]: =k2; kordx [ntetr+5,1]: =x; kordy [ntetr+5,1]: =y+hy [j] ; kordz [ntetr+5,1]: =z;

cells [ntetr+5,2]: =point+kz; kordx [ntetr+5,2]: =x; kordy [ntetr+5,2]: =y+hy [j] ; kordz [ntetr+5,2]: =z-hz [k] ;

cells [ntetr+5,3]: =point; kordx [ntetr+5,3]: =x; kordy [ntetr+5,3]: =y; kordz [ntetr+5,3]: =z-hz [k] ;

cells [ntetr+5,4]: =k1; kordx [ntetr+5,4]: =x+hx [i] ; kordy [ntetr+5,4]: =y; kordz [ntetr+5,4]: =z-hz [k] ;

V [k+ (j-1) * (kz-1) + (i-1) * (kz-1) * (ky-1)]: = (hz [k] *hy [j] *hx [i]) /6;

inc (point);

inc (ntetr,6);

end;

inc (point);

y: =y+hy [j] ;

end;

inc (point,kz);

x: =x+hx [i] ;

end;

алгоритм построения матрицы жесткости и ее компрессии

begin

n: =kx*ky*kz*3;

m: = (ky*kz+1) *3;

for k: =1 to koltetr do

for i: =1 to 4 do

for j: =1 to 4 do

if (cells [k, i] <=cells [k,j]) then formprmatr (cells [k, i],cells [k,j],k);

Count: =0;

for i: =1 to 3*kx*ky*kz do if P [i] =1 then compressf (i-Count);

dec (n,count);

end;

алгоритм решения системы уравнений

begin

s [1,1]: =sqrt (abs (value2 (1,1))); d [1]: =sign (value2 (1,1));

for j: =2 to n do s [1,j]: =value2 (1,j) / (d [1] *s [1,1]);

for i: =2 to n do

begin

for j: =2 to n do

begin

if i>j then s [i,j]: =0;

sum: =0;

for k: =1 to i-1 do sum: =sum+s [k, i] *s [k, i] *d [k] ;

d [i]: =sign (value2 (i, i) - sum);

if i

begin

sum: =0;

for k: =1 to i-1 do sum: =sum+s [k, i] *s [k,j] *d [k] ;

s [i,j]: = (value2 (i,j) - sum) / (s [i, i] *d [i]);

end;

sum: =0;

for k: =1 to j-1 do sum: =sum+d [k] *s [k,j] *s [k,j] ;

s [j,j]: =sqrt (abs (value2 (j,j) - sum));

end;

end;

y [1]: =F [1] / (s [1,1] *d [1]);

for k: =2 to n do

begin

sum: =0;

for i: =1 to k-1 do sum: =sum+s [i,k] *y [i] *d [i] ;

y [k]: = (F [k] -sum) / (s [k,k] *d [k]);

end;

x [n]: =y [n] /s [n,n] ;

for k: =n-1 downto 1 do

begin

sum: =0;

for i: =k+1 to n do sum: =sum+s [k, i] *x [i] ;

x [k]: = (y [k] -sum) /s [k,k] ;

end;

end;