МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ
УЧРЕЖДЕНИЕ ОБРАЗОВАНИЯ “ГОМЕЛЬСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ им.Ф. СКОРИНЫ”
Математический факультет
Кафедра ВМ и программирования
Курсовой проект
Суперэлементное моделирование пространственной системы “плита - грунтовое основание ”
Исполнитель
студент группы ПМ-44 Рыжик И.А.
Научный руководитель Цурганова Л.А.
к. т. н.
ГОМЕЛЬ 2001

Содержание
  Введение
1. Системы и методы их исследования. Системный подход
1.1 Основные положения общей теории систем
1.2 Классификация систем
1.3 Структура системы
1.4 Системный подход.
1.5 Методы исследования систем
2. Основные понятия теории упругости
2.1 Напряжения
2.2 Деформации
3. Основная концепция метода конечных элементов
4. Характеристики тетраэдрального элемента
4.1 Функции перемещений
4.2 Матрица деформации
4.3 Матрица упругости
4.4 Матрицы жесткости, напряжений и нагрузок
5. Математическая и дискретная модели
5.1 Математическая модель
5.2 Дискретная модель
6 Алгоритмы построения и решения дискретной модели
7. Описание и инструкция работы с приложением
8. Верификация приложения
Заключение
Список использованных источников
Приложение

Введение

В гражданском и промышленном строительстве все чаще под застройку идут территории, содержащие неоднородный грунт. В связи с этим актуальна задача расчета осадок плиты с учетом неоднородности грунтового основания. Цель таких расчетов выявить ослабленные места в грунтовом основании на этапе проектирования и предложить дополнительные мероприятия по подготовке территории под строительство или изменить форму фундамента.
В настоящей работе рассматривается плита на неоднородном линейно-деформируемом основании. Нагрузка на плиту берется вертикальная, равномерно распределенная. Моделируется расчет осадок фундаментной плиты с учетом сложной структуры неоднородного основания.
Математическая модель системы “плита - грунтовое основание" представляет собой третью краевую задачу математической физики. Она описывает условия равновесия системы. Равновесие систем механики твердого деформируемого тела может быть описано уравнениями равновесия в напряжениях, перемещениях либо каким-то вариационным принципом, например принципом минимума полной энергии системы.
Т.к. МКЭ эффективен лишь сравнительно для небольших систем, то для решения математической модели применяется метод суперэлементов. Метод суперэлементов основывается на той же теоретической базе, что и МКЭ, только предварительно ещё используется метод декомпозиции, т.е. вся расчётная область разбивается на отдельные макроэлементы, называемые суперэлементами.
Суперэлементное моделирование системы “плита - грунтовое основание" включает в себя построение и решение дискретной модели. Построение реализуется алгоритмами построения матрицы жесткости, заданием вектора нагрузок и граничных условий для отдельного суперэлемента.
В разработанных алгоритмах учитываются особенности матриц жесткости отдельного суперэлемента и неоднородность грунтового основания системы.
Для удобства пользователя спроектирован интерфейс вывода исходных данных: размеров нерегулярной решетки, выбор характеристик конечных элементов по слоям XOZ (для каждого суперэлемента), интерфейс вывода результатов в табличной форме.
Приложение моделирования расчета осадок плиты реализуется в интегрированной среде программирования Borland Delphi 5.0.

1. Системы и методы их исследования. Системный подход

1.1 Основные положения общей теории систем

Под системой понимают конечную совокупность элементов, связей между ними и между их свойствами, действующими как целостное образование для достижения единой цели. Элементом называют некоторый объект (биологический, информационный, энергетический материал), обладающий рядом определенных свойств, но внутреннее строение (содержание) которого безотносительно к цели рассмотрения. Элементы будем обозначать через M, их совокупность через {M}. Связью называется важный для рассмотрения обмен между элементами (веществом, энергией, информацией и т.п.), т.е. фактор, связывающий элементы и их свойства в целое. Единичным фактом связи выступает воздействие , где i, j - индексы взаимодействующих элементов , . Связи позволяют по свойствам перехода по ним от элемента к элементу соединить два элемента совокупности. Свойства есть качества параметров объектов. Они могут изменяться в результате действия системы. Свойства дают возможность описывать объекты системы количественно. Любая система характеризуется двумя признаками: связанностью (наличием связи между элементами); функциализацией (свойства системы отличаются от свойств отдельных элементов). Применяя ”картежные” определения системы, символически систему можно записать в следующем виде:
 {{M}, {X}, F} (1), где  - система,
{M} - совокупность элементов, {X} - совокупность связей,
F - функция системы.
Запись типа (1) является наиболее простой и достаточно полной.

1.2 Классификация систем

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

1.3 Структура системы

Система может иметь структурное представление, т.е. может быть расчленена на группы элементов с указанием связи между ними. Такое расчленение называется декомпозицией. Декомпозиция на время изучения сохраняется неизменной. Группы элементов называются модулями системы. Они образуются по принципу общих свойств, а также по характеру связей между ними и по другим признакам. Символически структура может быть записана в следующем виде:
: {{M},{X}}  (2), где  - система,
{M}-совокупность модулей, {X}-совокупность связей.

1.4 Системный подход.

Известно, что свойства системы, как сложного объекта, не обнаруживаются в свойствах её отдельных подсистем. Это значит, что традиционный метод изучения целого путём анализа его частей и последующего объединения (суперпозиции) их свойств непригоден для больших и сложных систем. Решением проблемы становится системный подход, суть которого состоит во взаимосвязанном рассмотрении всех элементов (подсистем) системы. При системном подходе система рассматривается не изолированно, а как подсистема более общей системы (системы более высокого ранга). Основным при системном подходе является определение цели, например, определение способов достижения равновесия деформируемой системы; снижение материалоёмкости конструктивных элементов механизмов и т.п. Для каждой цели должен быть выбран свой надёжный критерий эффективности. Например, для информационных систем это может быть оперативность информации, её полнота, надёжность и прогнозируемость развития процессов, входящих в область интересов системы.
Системный подход при исследовании различных систем, явлений, объектов позволяет с единых позиций строить общую методологию исследования указанных систем и процессов независимо от их природы. Эта методология, как и любая другая, содержит определенные этапы.
Этап 1. Определение системы.
а) Определение исследуемой функции системы.
б) Определение области существования системы вместе с ее границей.
в) Определение краевых условий.
г) Декомпозиция системы вплоть до простых элементов.
д) Определение свойств элементов и модулей системы
е) Нахождение связей между элементами и модулями исходной системы.
Этап 2. Построение математической модели.
а) Формальное описание исследуемой функции.
б) Разработка дискретной модели системы.
в) Разработка алгоритмической модели.
г) Проверка адекватности математической модели системы.
Этап 3. Исследование поведения системы при различных входных воздействиях.

1.5 Методы исследования систем

В настоящее время существует несколько методов исследования систем.
Микроподход.
Суть этого метода сводится к исследованию отдельных элементов системы. Выбор этих элементов не однозначен и определяется задачей исследования или системой. При использовании микро подхода изучается структура каждого из выделенных элементов системы, их функции, совокупность и диапазон возможных изменений параметров, после чего делается попытка понять процесс функционирования системы в целом.
Задачи микроподхода заключаются в следующем:
выявление элементов исследуемой системы;
изучение структуры выделенных элементов;
раскрытие функций каждого из элементов;
выявление связей между элементами.
Макроподход.
При этом методе система рассматривается как “черный ящик”, внутреннее строение которого неизвестно. Такая ситуация может быть, например, при изучении недоступных управляющих систем или исследование систем, структура которых изучена недостаточно. В процессе макро подхода исследователь, воздействуя различным образом на вход системы, анализирует ее реакцию на соответствующие входные воздействия. Имея обширную статистическую информацию вследствие ее анализа делается вывод о структуре системы и принципов ее функционирования.
Физическое моделирование.
Это моделирование осуществляется путем воспроизведения исследуемого процесса на модели, имеющий в общем случае отличную от оригинала природу, но одинаковое математическое описание процесса функционирования. При этом физические процессы, протекающие в модели и оригинале, являются подобными. Физическое моделирование позволяет провести исследование процессов и систем, непосредственный анализ которых затруднен или не возможен. Использование физической модели позволяет определить влияние различных параметров на протекание изучаемых процессов, уточнить структуру системы и понять принцип ее функционирования.
Математическое моделирование.
Математическая модель концентрирует в себе описанную в форме математических соотношений совокупность наших знаний, представлений и гипотез о соответствующих объектах или знаниях.
Т. к. знания никогда не бывают абсолютными, а в гипотезах иногда намеренно не учитываются некоторые эффекты, то модель лишь приближенно описывает поведение реальной системы.
Основное назначение модели - это возможность сделать некоторые выводы о поведении реальной системы. Наблюдения над реальной системой (натурные эксперименты) в лучшем случае могут дать материал лишь для проверки той или иной гипотезы, той или иной модели, т.к они представляют собой источник информации ограниченного объема о прошлом этой системы.
Модель допускает значительно более широкие исследования, результаты которых дают нам информацию для прогнозирования поведения системы. Чтобы обеспечить эти и другие возможности приходится решать проблему соотношения (адекватности) модели и системы, т.е. необходимо проводить дополнительные исследования согласованности результатов моделирования с реальными результатами. Создавая модель, исследователь познает систему, т.е. выделяет ее как объект изучения из окружающей среды и строит ее формальное описание в соответствии с поставленными целями, задачами и имеющимися возможностями.
В дальнейшем, анализируется система через поведение ее модели, в том числе делается, и прогнозирование ее функции во времени. Математические модели строят на основе законов и закономерностей, выявленных фундаментальными науками. Для всякой модели необходимо построить моделирующий алгоритм.
В целом можно сказать, что процесс моделирования сводится к трем объектам:
система (реальная, проектируемая, воображаемая);
математическая модель системы;
алгоритмическая модель системы.
В соответствии с этим возникают следующие задачи:
определение (формирование) системы исследования;
построение математической модели системы;
разработка алгоритм решения модели.

2. Основные понятия теории упругости

2.1 Напряжения

B

A

Будем рассматривать следующий случай: возьмем призматический стержень, который растягивается под действием сил, равномерно распределённых по его концевым сечениям, внутренние силы распределены равномерно по поперечному сечению стержня АВ, напряжения можно найти, разделив полное значение растягивающей силы Р на площадь поперечного сечения F.
В общем случае напряжения по сечению распределены неравномерно, чтобы определить значение напряжения в некоторой точке этой плоскости, возьмём элементарную площадку δF в окрестности данной точки и предположим, что силы, возникающие на этой площадке, сводятся к равнодействующей δP. Если теперь равномерно стянуть элементарную площадку δF, то в пределе получится отношение δP/δF, которое определит величину напряжения, возникающего на плоскости АВ в некоторой точке. Направление этого напряжения будет совпадать с направлением равнодействующей δP. В общем случае напряжение направлено под некоторым углом к площадке δF, на которой оно действует, и обычно раскладывается на две составляющие: нормальное напряжение, перпендикулярное к площадке δF, и касательное напряжение, действующее в плоскости площадки.

2.2 Деформации

При рассмотрении деформации в упругом теле предполагается, что
Существуют ограничения, препятствующие перемещению его как жёсткого тела. Таким образом, какое-либо перемещение частиц тела может происходить лишь за счёт его деформации. Малые перемещения частиц при деформировании тела разложим по составляющим u, v, параллельные соответствующим осям координат x, y. Можно предположить, что эти малые величины непрерывно изменяются по всей площади тела.
Рассмотрим бесконечно малый элемент dxdy вблизи точки О тела.
Можно показать, что относительное удлинение по направлению оси y задается производной.
Рассмотрим теперь изменение угла между отрезками ОА и ОВ, которые до деформирования тела были взаимно перпендикулярны. Если u и v - перемещения точки О в направлениях осей x и y, то перемещения точки А в направлении оси у и точки В в направлении оси х будут соответственно равны. Поэтому первоначально прямой угол АОВ между отрезками ОА и ОВ уменьшается на величину, которая представляет собой деформацию сдвига между осями х и у.

3. Основная концепция метода конечных элементов

Метод конечных элементов (МКЭ) основан на идее аппроксимации непрерывной функции (температуры, давления, перемещений и т.п.) дискретной моделью, которая строится на множестве кусочно-непрерывных функций, определённых на конечном числе подобластей, называемых конечными элементами. В качестве функции элемента, чаще всего, принимается полином. Порядок полинома определяется числом используемых в каждом узле элемента данных о непрерывной функции. В общем случае форма конечного элемента может быть произвольной, но для удобства математических выкладок их принимают правильной геометрической формы. Конечные элементы могут быть линейные и криволинейные, одномерные, двумерные и трехмерные. Количество узлов конечного элемента может быть равно или больше количества вершин. В зависимости от этого качества можно проводить классификацию конечных элементов. Выделяют следующие три группы: симплекс-, комплекс - и мультиплекс-элементы.
Симплекс элементам соответствуют полиномы, содержащие константу и линейные члены:
;
Здесь коэффициентов столько сколько узлов.
Комплекс-элементам соответствуют полиномиальные функции, содержащие константу и члены первого и более высоких порядков. Форма комплекс элемента может быть такой же как и у симплекс-элемента, но комплекс элементы имеют количество узлов больше количества вершин.
Интерполяционный полином для двумерного треугольного комплекс элемента имеет вид:

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

4. Характеристики тетраэдрального элемента

4.1 Функции перемещений

На фигуре 1 изображен тетраэдральный элемент ijpm в системе координат x, y, z.

Y

Z

X

i

j

p

m

Перемещение любой точки определяется тремя компонентами u, v, w в направлениях координат x, y, z. Таким образом, вектор перемещений имеет вид
. (1)
Если для задания линейного закона изменения какой-либо величины в плоском треугольном элементе требовались три узловых значения, то в трехмерном случае необходимо задать четыре узловых значения. По аналогии с представлением (4.3) можно записать, например,
. (2)
Приравнивая эти выражения перемещением узловых точек, получаем четыре уравнения типа
 и т.д. (3)
из которых определяются коэффициенты .
Запишем теперь соотношение (2) в следующей форме, с использованием определителя
 (4), где
 (5а)
Величина V в данном случае представляет собой объем тетраэдра. Коэффициентами  обозначены определители
 (5б)
Остальные коэффициенты получаются циклической перестановкой индексов i, j, p, m.
Как видно из фигуры 1, узлы i, j, p, m пронумерованы в соответствии с правилом правой руки, причем первые три узла пронумерованы по часовой стрелке, если смотреть со стороны последнего узла.
Перемещение элемента определяется двенадцатью компонентами перемещений его узлов:
 (6) где
 и т.д.
Перемещение произвольной точки можно записать в виде
 (7)
где скалярные величины определяются соотношениями
и т.д.
А I - единичная матрица размерности 3*3.
Ясно, что эти функции перемещений будут удовлетворять требованиям непрерывности на границах между элементами. Этот результат является прямым следствием линейного закона изменения перемещений.

4.2 Матрица деформации

В трехмерном случае учитываются все шесть компонент деформации. Используя известные обозначения Тимошенко, запишем матрицу деформаций в виде
 (9)
С помощью соотношений (4) - (7) легко убедиться, что
 (10) где
. (11)
Остальные подматрицы получаются простой перестановкой индексов.
Начальные деформации, такие, как обусловленные тепловым расширением, можно записать обычным образом в виде шестикомпонентного вектора, имеющего, например, для изотропного теплового расширения простой вид:
 (12)
где  - коэффициент линейного расширения, а - средняя по элементу температура.

4.3 Матрица упругости

В случае материала с изотропией свойств матрица [D], связывающая шесть компонент напряжения с компонентами деформации, может содержать не более чем 21 независимую постоянную.
В общем случае
. (13)
Так как такое умножение никогда не выполняется в явном виде, запишем здесь матрицу [D] только для изотропного материала, хотя это нетрудно сделать и для случая произвольной анизотропии. При использовании обычных упругих постоянных: модуля упругости Е и коэффициента Пуассона v - матрица имеет вид
 (14)

4.4 Матрицы жесткости, напряжений и нагрузок

Выражение для матрицы жесткости, определяемой в общем случае соотношением , можно проинтегрировать точно, так как компоненты деформации и напряжения постоянны внутри элемента.
Подматрица с индексами rs матрицы жесткости имеет размерность 3*3 и определяется соотношением
, (15)
где V - объем тетраэдра.
Узловые силы, обусловленные начальной деформацией, записываются в виде
, (16) или для i-ой компоненты
.

5. Математическая и дискретная модели

5.1 Математическая модель

Математическая модель системы включает геометрическую, структурную, механико-математическую модели, краевые условия и условия равновесия системы.
Геометрическая модель представляет собой параллелепипед, размеры которого определяются нулевыми перемещениями на его ребрах.
Механико-математическая модель системы “плита-основание”: для основания s_i=E _ie_i, для плиты s_i=E’e_i, E’>>E_i, где E’, E_i -модули упругости основания и плиты, s_i, e_i -интенсивности напряжений и деформаций.
Краевые условия области определения системы “плита-основание": перемещения на всех ребрах, кроме верхнего равны нулю, на верхнем ребре области определения на поверхности плиты задается внешняя нагрузка.

5.2 Дискретная модель

Процесс дискретизации разделяется на 2 этапа:
Разбиение области на подобласти. Подобласти характеризуются стационарностью определяющих характеристик: свойства материала, прилагаемая нагрузка.
Разбиение подобластей на конечные элементы. Подобласти разбиваются на симплекс-элементы.
Дискретизация производится элементами малых размеров. Деформация и напряжение в любом конечном элементе выражаются через перемещения по известным формулам. В узлах элементов вводятся силы, статистически эквивалентные напряжениям на границе соответствующего элемента и внешним силам, приложенным к нему.
Разбивка на элементы производится так, что в пределах одного элемента участок среды рассматривается как однородный. Любой другой элемент, оставаясь однородным, может характеризоваться свойствами, отличными от соседних элементов. Таким образом, система в целом представляет неоднородную среду.
Применение МКЭ для решения системы “плита-основание” приводит к системе линейных алгебраических уравнений с ленточной симметричной матрицей. Ширина ее полуленты зависит от порядка нумерации узлов и определяется по формуле: 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<j то
;
В том случае, если матрица A самосопряжённая и положительно определённая, то матрицу D можно опустить, так как она будет единичной. Метод осуществляется по следующей схеме:
сначала решаем уравнение S*Y=B
затем уравнение SX=Y, находя решение системы.
Наша работа заключается в решении СЛАУ методом квадратного корня, используя ленточную симметрическую матрицу, компактно упакованную.
Полуленточная матрица системы строиться следующим образом:

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

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

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

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

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

Z

X

Y

1

2

3

4

5

6

7

8

9

10

3

12

Однородная плита располагается вертикально на жёстком основание. Усилие Р равномерно распределено по верхнему основанию плиты. Дискретизация пластины производится путём разбиения ее на конечные элементы (тетраэдры). Узлы и полученные конечные элементы нумеруются. Программное приложение рассчитывает значения перемещений в каждом узле модели.
При введении в качестве параметров модели тестового примера следующих величин: нагрузка на плиту Р = 100кг,
параметры плиты: модуль упругости Е = 360 кг/см², коэффициент Пуассона m=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) or (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 comdivssf (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 comdivssf (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<j then
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;