Р РК 218-149-2018 МЕТОДИЧЕСКИЕ РЕКОМЕНДАЦИИ ПО ПРАКТИЧЕСКОМУ ИСПОЛЬЗОВАНИЮ ПРОГРАММЫ ПО ОПРЕДЕЛЕНИЮ НАПРЯЖЕННО-ДЕФОРМИРОВАННОГО СОСТОЯНИЯ МНОГОСЛОЙНОЙ КОНСТРУКЦИИ ДОРОЖНОЙ ОДЕЖДЫ

Приказ Председателя Комитета автомобильных дорог Министерства по инвестициям и развитию Республики Казахстан от 21 декабря 2018 года № 122

Предисловие

1	РАЗРАБОТАНЫ И ВНЕСЕНЫ	Акционерным обществом "Казахстанский дорожный научно-исследовательский институт" (АО "КаздорНИИ")
2	УТВЕРЖДЕНЫ И ВВЕДЕНЫ В ДЕЙСТВИЕ	Приказом Председателя Комитета автомобильных дорог Министерства по инвестициям и развитию Республики Казахстан № 122 от 21 декабря 2018 года
3	СОГЛАСОВАНЫ	Акционерным обществом "НК "ҚазАвтоЖол" № 03/14-2-2591-И от 12 ноября 2018 года
4	СРОК ПЕРВОЙ ПРОВЕРКИ	2023 год
	ПЕРИОДИЧНОСТЬ ПРОВЕРКИ	5 лет
5	ВВЕДЕНЫ ВПЕРВЫЕ

Документ доступен к просмотру в информационно-правовой системе нормативно-правовых актов Республики Казахстан "Әдiлет"

Настоящие рекомендации не могут быть полностью или частично воспроизведены, тиражированы и распространены без разрешения Комитета автомобильных дорог Министерства по инвестициям и развитию Республики Казахстан

Содержание

	Введение
1	Область применения
2	Нормативные ссылки
3	Термины и определения
4	Методика использования программы BASIC_NDS_MKE_8_uzlov
4.1	Общая постановка задачи
4.2	Формирование массивов координат
4.3	Граничные условия
4.4	Определяющие соотношения уравнения равновесия
4.5	Вычисление матрицы градиентов четырехугольного элемента с восемью узлами
4.6	Алгоритм определения деформации и напряжений
4.7	Краткая инструкция по использованию программы BASIC_NDS_MKE_8_uzlov
4.8	Примеры расчетов по программе BASIC_NDS_MKE_8_uzlov
Приложение А (обязательное) Исходный код программы BASIC_NDS_MKE_8_uzlov Приложение Б (обязательное) Исходный код подпрограммы BtDBnds Библиография

Введение

Настоящие методические рекомендации разработаны согласно плану работ, принятому в рамках выполнения темы Б01.02 "Разработка теоретических решений и программ для определения напряженно-деформированного состояния многослойной конструкции дорожной одежды".

Документ содержит постановку типовой задачи по определению напряженно-деформированного состояния многослойной дорожной конструкции, статически нагруженной весом автомобиля, и основные формулы из алгоритма расчета методом конечных элементов. В методе применяется высокоточный квадратичный прямоугольный конечный элемент с восемью узлами. В приложении А приведен текст (листинг) расчетной программы BASIC_NDS_MKE_8_uzlovна языке MATLAB [1-4], снабженный необходимыми комментариями. В приложении Б приводится листинг подпрограммы BDB0, предназначенный для вычисления матрицы жесткости элемента. Приводится также краткая инструкция по использованию программы BASIC_NDS_MKE_8_uzlov. В конце документа приведен список использованной литературы.

1 Область применения

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

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

2 Нормативные ссылки

Для применения настоящих рекомендаций необходимы следующие ссылочные нормативные документы:

СП РК 3.03-103-2014 "Проектирование жестких дорожных одежд"

СП РК 3.03-104-2014 "Проектирование нежестких дорожных одежд"

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

3 Термины и определения

В настоящих рекомендациях применяются следующие термины с соответствующими определениями:

3.1 Дорожная одежда: Многослойная конструкция в пределах проезжей части автомобильной дороги, воспринимающая нагрузку от автотранспортного средства и передающая ее на грунт.

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

3.3 Деформация: Относительная величина, характеризующая изменение линейных размеров тела по отношению к его первоначальным размерам.

3.4 Напряжение: Относительная величина, определяемая нормальной или касательной нагрузкой, приходящейся на единицу площади.

3.5 Дорожная одежда нежесткая: Дорожная одежда со слоями, устроенными из разного вида асфальтобетонов, из материалов и грунтов укрепленных битумом, цементом, известью, комплексными и другими вяжущими, а также из слабосвязных зернистых материалов

3.6 Конструктивный слой: Каждый слой дорожной одежды, состоящий из однородных материалов и отличающийся от соседних слоев видом материалов, его прочностью и составом. Учитывается при расчете прочности дорожной одежды.

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

3.8 Основание дорожное: Нижний несущий слой дорожной одежды, воспринимающий нагрузки от транспортных средств совместно с покрытием и предназначенный для ее распределения на дополнительные слои или непосредственно на грунт земляного полотна.

3.9 Метод конечных элементов: Численный метод решения дифференциальных уравнений математической физики.

4 Методика использования программы BASIC_NDS_MKE_8_uzlov

4.1 Общая постановка задачи

Для построения математической модели автомобильной дороги общего пользования рассмотрим конструкцию, состоящую из двухслойной асфальтобетонной дорожной одежды и основания дорожной одежды, состоящего из щебеночно-песчаной смеси с 8 % цементной добавкой, щебеночной смеси и гравийно-песчаной смеси. Основание дорожной одежды устраивается на грунтовом основании из легкого суглинка (рисунок 1).

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

Поэтому, приведенные в таблице 1 геометрические размеры конструктивных слоев в рассматриваемой математической модели (рисунок 1) и физико-механические свойства материалов конструктивных элементов (в том числе модули упругости) назначены согласно [6].

1 – покрытие, асфальтобетон мелкозернистый плотный; 2 – асфальтобетон крупнозернистый пористый; 3 – щебеночно-песчаная смесь 8 %; 4 – щебеночная смесь; 5 –гравийно-песчаная смесь; 6 – грунт суглинок легкий

Рисунок 1 – Схематический вид дорожной конструкции

Значения физико-механических свойств материалов конструктивных элементов вводятся в 108-110 строках приложения А.

Таблица 1 – Характеристики конструктивных слоев

Материал конструктивного слоя	Модуль упругости Е, МПа	Коэффициент Пуассона,	Высота слоя, м
1 Асфальтобетон мелкозернистый плотный	3200	0,18	0,05
2 Асфальтобетон крупнозернистый пористый	2000	0,27	0,10
3 Щебеночно-песчаная смесь 8%	800	0,30	0,40
4 Щебеночная смесь	275	0,30	0,45
5 Гравийно-песчаная смесь	180	0,30	1,00
6 Грунт суглинок легкий	50	0,35	-

Для компьютерной реализации алгоритма решения задачи о напряженно-деформированном состоянии дорожной конструкции методом конечных элементов разработана расчетная программа BASIC_NDS_MKE_8_uzlovна алгоритмическом языке системы MATLAB.

Исследуемая область разбивается на квадратичные прямоугольные конечные элементы с восемью узлами (рисунок 2). Местная нумерация узлов с 1 по 8 в таком конечном элементе начинается с левого нижнего узла, и осуществляется в направлении против часовой стрелки.

Рисунок 2 – Прямоугольный элемент с восемью узлами

Для создания конечно-элементной сетки применялись переменные шаги и в направлении координатных осей и (рисунок 3).

В каждом горизонтальном ряде (рисунок 3) расположено по 60 элементов (egor=60), а количество горизонтальных рядов элементов – 23 (ever=23). Следовательно, исследуемая область разбита на 1380 конечных элементов, а число узлов равно np=(3*ever+2)*egor+2*ever+1=4307 (строки 12-15 приложения А).

Рисунок 3 – Схематический вид конечно-элементной разбивки исследуемой области

Глобальная нумерация узлов начинается с крайнего левого вертикального ряда узлов, и осуществляется сверху вниз и перемещается слева направо (рисунок 3).

При этом разбивка получается нерегулярной, т.к. в восьмиузловых элементах центральный узел отсутствует. В данном случае, в одном вертикальном ряде имеются 47 узлов, а в следующем – 24 узла.

Требуется определить компоненты вектора перемещений в узловых точках конечно-элементной сетки под действием веса автомобиля, передаваемого через одно колесо и приложенного в виде вертикальных нагрузок в соответствующих узлах (рисунок 4).

L и H – ширина и высота исследуемой области соответственно; P – нагрузка от колеса транспортного средства; u_i, v_i – компоненты перемещения точек на границах в направлении координатных осей x и y; I – VI конструктивные слои дорожной одежды и грунтового основания

Рисунок 4 – Расчетная схема задачи

4.2 Формирование массивов координат

В программе BASIC_NDS_MKE_8_uzlovзаданы массивы dx и dy, содержащие значения шагов по координатам. Наличие таких массивов позволяет задавать переменные шаги по координатам (строки 44-63 приложения А).

В блоке формирования координат программы (строки 68-91 приложения А) глобальные номера узлов текущего конечного элемента обозначены через n1, n2, …, n8, соответствующие локальным номерам 1,2, …, 8, приведенным на рисунке 2. Конкретные значения глобальных номеров узлов, вычисляемые по алгоритмам вычисления n1, n2, …, n8 в программе BASIC_NDS_MKE_8_uzlovопределяются номером элемента в текущем горизонтальном ряде узлов (m=1:egor) и номером элемента в текущем вертикальном ряде узлов (n=1:ever).

4.3 Граничные условия

Из рисунка 4 видно, что граничные условия в перемещениях по бокам и на нижней границе исследуемой области задаются в виде нулевых значений соответствующих компонентов вектора перемещений (u=0 и u=v=0). Чтобы задавать граничные условия в виде известных перемещений некоторых узлов также нужно создать специальный массив Mz, содержащий номера узлов задаваемых известных перемещений (строки 23-43 приложения А).

Граничные условия в нагрузках задаются в виде сосредоточенных вертикальных усилий, рассчитанных от веса автомобиля, приходящегося на одно колесо. Методика расчета граничных условий в нагрузках приведена в приложении А (строки 112-166). Конечный результат расчета представляет собой вектор правой части уравнения равновесия системы (1).

4.4 Определяющие соотношения уравнения равновесия

Уравнение равновесия системы в случае применения квадратичных элементов, как и в случае линейных треугольных конечных элементов, в матричной форме имеет вид [7]:

(1)

где [K] – матрица жесткости системы размерности (2*np,2*np);

– вектор узловых перемещений размерности 2*np;

– вектор узловых нагрузок размерности 2*np;

np – общее число узлов в конечно-элементной сетке.

Порядок вычисления значений компонентов вектора , соответствующие граничным условиям в нагрузках приведены в строках 150-165приложения А. Остальные компоненты этого вектора, в силу условий равновесия узлов, равны нулю.

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

где – общее число элементов в конечно-элементной сетке; – матрица жесткости текущего элемента "".

Матрица жесткости элемента вычисляется с помощью объемного интеграла (2):

(2)

где алгоритм вычисления матрицы упругости приведены в строках 204-207 приложения А.

В выражении (2) так называемая матрица градиентов формируется путем дифференцирования функций формы по координатам .

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

4.5 Вычисление матрицы градиентов четырехугольного элемента с восемью узлами

Для численного интегрирования объемного интеграла в соотношении (2) предлагается [7] переходить к локальной системе координат (рисунок 5). Здесь замену переменных интегрирования можно сделать с помощью соотношения (3):

(3)

где – единичная толщина элемента, а – определитель матрицы преобразования координат Якоби. Тогда объемный интеграл в (2) приводится к виду (4):

(4)

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

Рисунок 5 – Система координат для линейного прямоугольного элемента

4.6 Алгоритм определения деформации и напряжений

Выражение перепишем в следующем виде (5):

(5)

где

,

, ,

,

, (6)

,

, .

Теперь из выражения (5) компоненту тензора деформации определим следующим образом:

(7)

или учитывая равенства (5)

(8)

В выражении (8) координаты и относятся к линейному прямоугольнику (рисунок 5), и им при рассмотрении восьмиугольных прямоугольников будут соответствовать и .

Аналогично определяются и остальные компоненты тензора деформации. Из формулы (5) для компоненты имеем:

или, учитывая (4) и , получим:

(9)

Деформация сдвига определяется аналогичным образом.

Для определения компонентов напряжений воспользуемся формулами вычисления напряжений [7]:

, и , (10)

где, для плоской деформации ;

- модуль упругости;

- коэффициент Пуассона.

Программа численной реализации приведенного алгоритма вычисления компонентов вектора деформации и вектора напряжений приведена в строках 412-527 приложения А настоящего документа. Кроме этого, в приложении А предусмотрены фрагменты программ построения графиков горизонтальных перемещений точек характерных вертикальных сечений (строки 297-312 приложения А), вертикальных перемещений точек горизонтальных сечений, расположенных на различных глубинах (строки 316-408 приложения А) и значений компонент векторов деформаций и напряжений (строки 441-475, 529-602 приложения А).

4.7 Краткая инструкция по использованию программы

BASIC_NDS_MKE_8_uzlov

Целью создания настоящего нормативного документа является разработка расчетной программы на языке MATLAB.

Расчетная программа предназначена для вычисления значений компонентов перемещений, деформаций и напряжений в точках исследуемой дорожной конструкции (рисунок 4).

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

Для запуска программы BASIC_NDS_MKE_8_uzlovдля решения конкретной задачи,в первую очередь необходимо, чтобы в компьютере была установлена одна из версии программного комплекса MATLAB,и в директорию этого комплекса был внесен текст (листинг) программы BASIC_NDS_MKE_8_uzlov. На рабочий стол компьютера выносится ярлык программного комплекса MATLAB (рисунок 6).

Рисунок 6 – Ярлык программного комплекса MATLAB

В приложении А настоящего документа приведен исходный код программы BASIC_NDS_MKE_8_uzlov. Для облегчения обращения к тексту листинг снабжен построчной нумерацией, которую нужно удалить перед установкой программы в директорию программного комплекса MATLAB.

Вызов программного комплекса MATLAB для расчетов осуществляется двойным щелчком на ярлыке (рисунок 6).

На рабочем столе появится главная страница MATLAB (рисунок 7).

Рисунок 7 – Вид главной страницы программы наMATLAB

На рисунке 7 на главной странице MATLAB приведен фрагмент программы BASIC_NDS_MKE_8_uzlov. Если же этот фрагмент будет отсутствовать, а на главной странице присутствует фрагмент другой программы, то поиск необходимой программы BASIC_NDS_MKE_8_uzlov будет произведен простым нажатием на кнопку OPEN, расположенной в верхней части главной страницы. После нажатия кнопки OPEN получим список программ, находящихся в директории MATLAB (рисунок 8).

Рисунок 8 – Вид списка программ в директории

Простым нажатием на название нужной программы BASIC_NDS_MKE_8_uzlovиз приведенного списка вызывается эта программа на главную страницу (рисунок 7).

Рисунок 9 – Строка 109-110 для введения физико-механических характеристик материалов слоев

По умолчанию в программе установлены физико-механические характеристики автомобильной дороги.Для изменения их значений данные необходимо ввести в строки 109-110 (рисунок 9). После введения физико-механических характеристик вводим значение прилагаемой вертикальной силы на поверхность в строку 150,а так же вводим номера поверхностных узлов куда необходимо приложить нагрузку (рисунок10).

Рисунок 10 – Строка для введения значения прилагаемой силы на поверхность и номеров узлов приложения сил

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

Рисунок 11 – Получение графиков вертикальных перемещений узлов на глубине h= 5 см

Рисунок 12 – Получение графиков вертикальных перемещений узлов на глубине h= 15 см

Рисунок 13 – Получение графиков вертикальных перемещений узлов на глубине h= 25 см

Рисунок 14 – Получение графиков вертикальных перемещений узлов на глубине h= 45 см

Рисунок 15 – Получение графиков компонент напряжений SigmaX

Рисунок 16 – Получение графиков компонент напряжений SigmaY

Рисунок 17 – Получение графиков компонент напряжений TauXY

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

Рисунок 18 – Кнопка Run на панели инструментов для запуска программы

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

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

4.8 Примеры расчетов по программе BASIC_NDS_MKE_8_uzlov

Ниже приводятся результаты компьютерной реализации программы BASIC_NDS_MKE_8_uzlovв виде графиков перемещений и напряжений.

1–горизонтальные смещения точек вертикального сечения, отстоящего от вертикальной оси симметрии слева на 30 см; 2–горизонтальные смещения точек вертикальной оси симметрии; 3–горизонтальные смещения точек вертикального сечения, отстоящего от вертикальной оси симметрии справа на 30 см

Рисунок 19 – Горизонтальные смещения точек вертикального сечения

Рисунок 20 – Вертикальные смещения точек горизонтальных сечений, расположенных на глубинах: 1– h=5 см; 2–h=15 см; 3–h=25 см; 4–h=45 см; 5–h=80 см.

Рисунок 21 – Распределение напряжений на глубинах h

Рисунок 22 – Распределение напряжений на глубинах h

Рисунок 23 – Распределение напряжений на глубинах h

Приложение А

(обязательное)

Исходный код программы BASIC_NDS_MKE_8_uzlov

1	% БАЗОВАЯ ВЕРСИЯ ПРОГРАММЫ BASIC_NDS_MKE_8_uzlov
2	% РЕШЕНИЕ ТЕСТОВОЙ ЗАДАЧИ О НДС С ПОМОЩЬЮ
3	%ЧЕТЫРЕХУГОЛЬНОГО
4	% КВАДРАТИЧНОГО ЭЛЕМЕНТА С 8-Ю УЗЛАМИ МЕТОДОМ КОНЕЧНЫХ
5	% ЭЛЕМЕНТОВ
6
7	% Начало разработки 8.07.2017 г
8	% Окончание 1.09.2018 г
9
10	Tic
11
12	ever=23; egor=60; % количество конечных элементов в вертикальном (ever) и
13	% горизонтальном (egor) рядах элементов в конечно-элементной
14	% разбивке
15	np=(3*ever+2)*egor+2*ever+1; % общее число узлов
16	x=zeros(np); y=zeros(np);
17
18	% ПОРЯДОК НУМЕРАЦИИ ЭЛЕМЕНТОВ - СЛЕВА НАПРАВО (РЯДЫ % НАЧИНАЮТСЯ СНИЗУ)
19	% ПОРЯДОК НУМЕРАЦИИ УЗЛОВ - СЛЕВА НАПРАВО. РЯДЫ УЗЛОВ ТАКЖЕ % НАЧИНАЮТСЯ
20	% СНИЗУ. ПРОМЕЖУТОЧНЫЕ РЯДЫ УЗЛОВ, ПРОХОДЯЩИЕ ЧЕРЕЗ СЕРЕДИНУ % ЭЛЕМЕНТОВ
21	% СОДЕРЖАТ egor+1 УЗЛОВ
22
23	% ФОРМИРОВАНИЕ МАССИВА НОМЕРОВ УЗЛОВ С ЗАДАННЫМИ ЗНАЧЕНИЯМИ
24	% ПЕРЕМЕЩЕНИЙ ДЛЯ ЗАДАНИЯ ГРАНИЧНЫХ УСЛОВИЙ В ПЕРЕМЕЩЕНИЯХ
25
26	Mz=zeros(400);
27	np=ever*(3*egor+2)+2*egor+1;
28	for i=1:(2*egor+1)
29	Mz(i)=i;
30	End
31	k=2*egor+2;
32	for i=1:ever
33	Mz(k)=(2*egor+1)+1+(i-1)*(3*egor+2);
34	Mz(k+1)=Mz(k)+egor;
35	Mz(k+2)=Mz(k+1)+1;
36	Mz(k+3)=Mz(k+2)+2*egor;
37	k=k+4;
38	End
39	for i=1:(2*egor+1)
40	Mz(k)=Mz(i)+np;
41	k=k+1;
42	End
43
44	% МАССИВЫ ШАГОВ ПО КООРДИНАТНЫМ ОСЯМ В МЕТРАХ
45	format long
46	dx=zeros(2*egor);
47	dy=zeros(2*ever);
48
49	dx=[
50	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
51	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
52	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
53	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
54	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
55	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
56	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
57	0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 0.250e-01 ...
58	0.250e-01 0.250e-01 0.250e-01 0.250e-01];
59
60	dy=[0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 ...
61	0.1000 0.1000 0.1000 0.1000 0.1000 0.0750 0.0750 ...
62	0.0625 0.0625 0.0500 0.0500 0.0500 0.0500 0.0250 ...
63	0.0250 0.0250];
64
65	yy=[0.000 0.200 0.400 0.600 0.800 1.000 1.200 1.400 1.600 1.800 2.000 2.200 ...
66	2.400 2.550 2.700 2.825 2.950 3.050 3.150 3.250 3.350 3.400 3.450 3.500];
67
68	% БЛОК ФОРМИРОВАНИЯ КООРДИНАТ
69	% Положительное направление оси х - слева направо
70	% Положительное направление оси у - снизу вверх
71	k=1;
72	NN=zeros(ever);
73	for n=1:ever
74	for m=1:egor
75	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
76	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
77	n6=n5-1; n7=n5-2; n8=n4-1;
78
79	x(n1)=(m-1)*2*dx(m); x(n2)=x(n1)+dx(m); x(n3)=x(n1)+2*dx(m);
80	x(n4)=x(n3); x(n5)=x(n3); x(n6)=x(n2); x(n7)=x(n1); x(n8)=x(n1);
81
82	y(n1)=yy(n); y(n2)=y(n1); y(n3)=y(n1); y(n4)=y(n3)+dy(n);
83	y(n5)=y(n3)+2*dy(n); y(n6)=y(n5); y(n7)=y(n5); y(n8)=y(n4);
84	if(m==36)
85	NN(2*k-1)=n3;
86	NN(2*k)=n4;
87	k=k+1;
88	end
89	End
90	End
91	NN(1:2*ever+1);
92
93	% ИСХОДНЫЕ ДАННЫЕ И ОБНУЛЕНИЕ МАТРИЦ
94	np2=2*np;
95	K=zeros(np2,np2);
96	U=zeros(np2+1);
97	U0=zeros(np2+1);
98	F=zeros(np2);
99
100	% НАИМЕНОВАНИЕ МАТЕРИАЛОВ КОНСТРУКТИВНЫХ СЛОЕВ (СВЕРХУ ВНИЗ)
101	% 1. АБ МЗ ПЛОТНЫЙ
102	% 2. АБ КЗ ПОРИСТЫЙ
103	% 3. ЩЕБЕНОЧНО-ПЕСЧАНАЯ СМЕСЬ 8%
104	% 4. ЩЕБЕНОЧНАЯ СМЕСЬ
105	% 5. ГРАВИЙНО-ПЕСЧАНАЯ СМЕСЬ
106	% 6. ГРУНТ СУГЛИНОК ЛЕГКИЙ
107
108	% ФИЗИКО-МЕХАНИЧЕСКИЕ ХАРАКТЕРИСТИКИ МАТЕРИАЛОВ СЛОЕВ(СВЕРХУ % ВНИЗ)
109	E=[6933e+05 4142e+05 1000e+05 250e+05 180e+05 58e+05]; % Модуль упр.кг/м^2
110	nu=[0.18 0.27 0.30 0.30 0.30 0.35]; % КОЭФФИЦИЕНТЫ ПУАССОНА
111
112	% ГРАНИЧНЫЕ УСЛОВИЯ В НАГРУЗКАХ
113
114	% Нагрузка от одиночного колеса
115
116	% ДАННЫЕ ИЗ ИНТЕРНЕТА
117	% 1 Па=0.102 кгс/м^2
118	% 0.6 МПа=6.12 кгс/см^2=60000.0 кгс/м^2
119
120	% РАСПРЕДЕЛЕННАЯ ПО УГЛОВЫМ УЗЛАМ ЭЛЕМЕНТА НАГРУЗКА ОТ % ОДИНОЧНОГО КОЛЕСА
121	% ПРИКЛАДЫВАЕТСЯ ПО ЦЕНТРУ РАСЧЕТНОЙ СХЕМЫ. При этом промежуточные % узлы
122	% игнорируются
123
124	% НАГРУЗКА НА КОЛЕСО РАССЧИТЫВАЕТСЯ ИЗ РАСЧЕТА НОРМАТИВНОЙ % ИНТЕНСИВНОСТИ
125	% НАГРУЗКИ q=0.6 MPa
126	% ЕСЛИ 1 Па=0.102 кгс/м^2 (из Интернета) ТО
127	% 0.6 МПа=6.12 кгс/см^2 = 60000.0 кгс/м^2.
128	% Таким образом, если на 1 погонный метр полотна приходится 60000 кг веса, то
129	% на одно расстояние 0.05 м – 3000 кг веса.
130	% Если ширину колеса принять равной 40 см, а ширину одного конечного
131	% элемента в горизонтальном направлении – 5 см, то транспортная нагрузка
132	% будет распределена по 9 узлам в пределах 8 конечных элементов.
133	% Расстояние между угловыми узлами - 0.05 м (5 cм - ширина конечного элемента).
134	% При этом надо учесть, что
135	% два крайних узла (1-ый и 9-ый узлы) будут нагружены наполовину.
136	% В таком случае, нагрузки в узловых точках будут распределены
137	% следующим образом:
138	% P1=P9=1500 кг, P2=P3= … =P8=3000 кг
139
140	% РАСЧЕТ НАГРУЗКИ НА УЗЕЛ
141	% РАССТОЯНИЕ МЕЖДУ УЗЛАМИ 0.025 М (учитываются и промежуточные узлы)
142	% НА 1 М ПРИХОДИТСЯ 60000 КГ ВЕСА, А НА ОДНО РАССТОЯНИЕ 0.025 М
143	% ПРИХОДИТСЯ 1500 КГ ВЕСА
144
145	% РЕЗУЛЬТАТЫ БУДУТ ПОЛУЧЕНЫ В КГ И МЕТРАХ
146
147	% ВЕЛИЧИНЫ ВЕРТИКАЛЬНЫХ СИЛ, ПРИЛОЖЕННЫХ ПО УГЛОВЫМ УЗЛАМ
148	% В ЦЕНТРАЛЬНОЙ ЧАСТИ ЭПЮРЫ
149
150	P=0.3e+04; % кг
151
152	% НОМЕРА УГЛОВЫХ УЗЛОВ ПРИЛОЖЕНИЯ СИЛ
153
154	nF(1:9)=[4239 4241 4243 4245 4247 4249 4251 4253 4255];
155	for i=1:9
156	j=nF(i);
157	xP(i)=x(j);
158	end
159
160	% ПРИЛОЖЕНИЕ СИЛ К УГЛОВЫМ УЗЛАМ
161	F(nF(1)+np)=-P/2;
162	F(nF(9)+np)=-P/2;
163	for i=1:7
164	F(nF(1+i)+np)=-P;
165	end
166
167	% ФОРМИРОВАНИЕ МАТРИЦЫ ЖЕСТКОСТИ ЭЛЕМЕНТА [ke] И
168	% ГЛОБАЛЬНОЙ МАТРИЦЫ ЖЕСТКОСТИ СИСТЕМЫ [K]
169
170	mpe=[-0.577350 0.577350]; % КООРДИНАТЫ ТОЧЕК ИНТЕГРИРОВАНИЯ
171	mH=[1 1]; % ВЕСОВЫЕ КОЭФФИЦИЕНТЫ
172	K=zeros(np2,np2);
173	for n=1:ever
174	for m=1:egor
175	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
176	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
177	n6=n5-1; n7=n5-2; n8=n4-1;
178
179	X1=x(n1); Y1=y(n1);
180	X2=x(n3); Y2=y(n3);
181	X3=x(n5); Y3=y(n5);
182	X4=x(n7); Y4=y(n7);
183
184	% ФОРМИРОВАНИЕ МАТРИЦЫ УПРУГОСТИ ДЛЯ СЛУЧАЯ ПЛОСКОЙ % ДЕФОРМАЦИИ
185	if((n>=1)&(n<=7))
186	e1=E(6); nu1=nu(6); % ГРУНТ СУГЛИНОК ЛЕГКИЙ
187	end
188	if((n>=8)&(n<=12))
189	e1=E(5); nu1=nu(5); % ГРАВИЙНО-ПЕСЧАНАЯ СМЕСЬ
190	end
191	if((n>=13)&(n<=16))
192	e1=E(4); nu1=nu(4); % ЩЕБЕНОЧНАЯ СМЕСЬ
193	end
194	if((n>=17)&(n<=20))
195	e1=E(3); nu1=nu(3); % ЩЕБЕНОЧНО-ПЕСЧАНАЯ СМЕСЬ 8%
196	end
197	if((n>=21)&(n<=22))
198	e1=E(2); nu1=nu(2); % АБ КЗ ПОРИСТЫЙ
199	end
200	if(n==23)
201	e1=E(1); nu1=nu(1); % АБ МЗ ПЛОТНЫЙ
202	end
203
204	D=zeros(3,3);
205	D(3,3)=e1/(2*(1+nu1)); D(2,2)=2*D(3,3)*(1-nu1)/(1-2*nu1);
206	D(1,1)=D(2,2); D(1,2)=2*D(3,3)*nu1/(1-2*nu1);
207	D(2,1)=D(1,2);
208
209	dJ=abs((X1-X2)*(Y2-Y4))/4; % якобиан
210
211	% ИСПОЛЬЗОВАНИЕ КВАДРАТУРЫ ГАУССА ДЛЯ ЧИСЛЕННОГО % ИНТЕГРИРОВАНИЯ
212	% ДВОЙНОГО ИНТЕГРАЛА В МАТРИЦЕ ЖЕСТКОСТИ ЭЛЕМЕНТА [Ke]
213
214	% BtDBnds1(psi,eta,D,X1,X2,Y1,Y3) - подпрограмма вычисления
215	% матрицы жесткости элемента
216	Ke=zeros(16,16);
217	psi=mpe(1); eta=mpe(1);
218	k11=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
219	psi=mpe(1); eta=mpe(2);
220	k12=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
221	psi=mpe(2); eta=mpe(1);
222	k21=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
223	psi=mpe(2); eta=mpe(2);
224	k22=BtDBnds1(psi,eta,D,X1,X2,Y1,Y3);
225	Ke=dJ*(mH(1)*(mH(1)*k11+mH(2)*k12)+mH(2)*(mH(1)*k21+mH(2)*k22));
226
227	ne8=[n1 n2 n3 n4 n5 n6 n7 n8];
228	for i=1:8
229	for j=1:8
230	K(ne8(i),ne8(j))=K(ne8(i),ne8(j))+Ke(i,j);
231	K(ne8(i),ne8(j)+np)=K(ne8(i),ne8(j)+np)+Ke(i,j+8);
232	K(ne8(i)+np,ne8(j))=K(ne8(i)+np,ne8(j))+Ke(i+8,j);
233	K(ne8(i)+np,ne8(j)+np)=K(ne8(i)+np,ne8(j)+np)+Ke(i+8,j+8);
234	end
235	end
236	end
237	end
238
239	disp('Матрица К сформирована')
240
241	% ПРЕОБРАЗОВАНИЕ СЛАУ
242
243	np2=2*np;
244	i=1;
245	for n=1:np2
246	if(n==Mz(i))
247	F(n)=K(n,n)*U(n);
248	for m=1:np2
249	if(m~=n)
250	K(n,m)=0.0;
251	end
252	end
253	i=i+1;
254	end
255	end
256
257	i=1;
258	for n=1:np2
259	if(n==Mz(i))
260	for m=1:np2
261	if(m~=n)
262	F(m)=F(m)-K(m,n)*U(n);
263	K(m,n)=0.0;
264	end
265	end
266	i=i+1;
267	end
268	end
269
270	% РЕШЕНИЕ СЛАУ
271	U=K\F;
272	disp('Система [K]{U}={F} решена')
273
274	% Номера узлов характерных вертикальных сечений
275	nom12=[
276	49 146 231 328 413 510 595 692 ...
277	777 874 959 1056 1141 1238 1323 1420 ...
278	1505 1602 1687 1784 1869 1966 2051 2148 ...
279	2233 2330 2415 2512 2597 2694 2779 2876 ...
280	2961 3058 3143 3240 3325 3422 3507 3604 ...
281	3689 3786 3871 3968 4053 4150 4235];
282	nom15=[
283	61 152 243 334 425 516 607 698 ...
284	789 880 971 1062 1153 1244 1335 1426 ...
285	1517 1608 1699 1790 1881 1972 2063 2154 ...
286	2245 2336 2427 2518 2609 2700 2791 2882 ...
287	2973 3064 3155 3246 3337 3428 3519 3610 ...
288	3701 3792 3883 3974 4065 4156 4247];
289	nom18=[
290	73 158 255 340 437 522 619 704 ...
291	801 886 983 1068 1165 1250 1347 1432 ...
292	1529 1614 1711 1796 1893 1978 2075 2160 ...
293	2257 2342 2439 2524 2621 2706 2803 2888 ...
294	2985 3070 3167 3252 3349 3434 3531 3616 ...
295	3713 3798 3895 3980 4077 4162 4259];
296
297	% ГРАФИКИ ГОРИЗОНТАЛЬНЫХ ПЕРЕМЕЩЕНИЙ УЗЛОВ В СРЕДНЕМ СЕЧЕНИИ
298	Y=zeros(2*ever+1);
299	X150=zeros(2*ever+1);
300	X120=zeros(2*ever+1);
301	X180=zeros(2*ever+1);
302	for i=1:2*ever+1
303	j120=nom12(i);
304	j150=nom15(i);
305	j180=nom18(i);
306	Y(i)=y(j150);
307	X120(i)=U(j120)-0.0001; % СЛЕВА ОТ ВЕРТИКАЛЬНОЙ ОСИ СИММЕТРИИ НА 30 СМ
308	X150(i)=U(j150)*10; % НА ВЕРТИКАЛЬНОЙ ОСИ СИММЕТРИИ
309	X180(i)=U(j180)+0.0001; % СПРАВА ОТ ВЕРТИКАЛЬНОЙ ОСИ СИММЕТРИИ НА % 30 СМ
310	end
311	hPlot=plot(X150,Y,'b-',X120,Y,'r-',X180,Y,'k-.');grid on;xlabel('Ux, m');ylabel('y, m')
312	set(hPlot,'LineWidth',2);
313
314	figure
315
316	% ФОРМИРОВАНИЕ МАССИВА НОМЕРОВ ГОРИЗОНТАЛЬНЫХ РЯДОВ УЗЛОВ % Nv(I,m)
317
318	% ОБЩЕЕ ЧИСЛО УЗЛОВ В КАЖДОМ ПОЛНОМ ГОРИЗОНТАЛЬНОМ РЯДУ РАВНО % 2*egor+1,
319	% А ЧИСЛО ТАКИХ ГОРИЗОНТАЛЬНЫХ РЯДОВ РАВНО ever+1
320	% СЕЧЕНИЯ, ПРОХОДЯЩИЕ ЧЕРЕЗ СЕРЕДИНУ ЭЛЕМЕНТОВ, НЕ % РАССМАТРИВАЮТСЯ.
321	% ОНИ ИСПОЛЬЗУЕТСЯ ДЛЯ ПОСТРОЕНИЯ ГРАФИКОВ ИСКОМЫХ ВЕЛИЧИН НА % КАЖДОМ
322	% ГОРИЗОНТАЛЬНОМ УРОВНЕ (ГЛУБИНЕ)
323
324	Nv=zeros(ever+1,2*egor+1);
325	for n=1:(ever+1)
326	for m=1:(2*egor+1)
327	I=2*n-1;
328	Nv(I,m)=m+(3*egor+2)*(n-1);
329	end
330	end
331	for n=1:ever
332	for m=1:(egor+1)
333	I=2*n;
334	Nv(I,2*m-1)=m+(2*egor+1)*n+(egor+1)*(n-1);
335	end
336	end
337
338	format short
339	Set=zeros(2*ever+1,4);
340	for n=1:2*ever+1
341	i=Nv(n,1);
342	Set(n,1)=n;
343	Set(n,2)=i;
344	Set(n,3)=y(i);
345	Set(n,4)=3.50-y(i);
346	end
347	for m=1:47
348	% fprintf('%6u %6u %10.4f %10.4f\n',Set(m,1),Set(m,2),Set(m,3),Set(m,4));
349	end
350
351	% ГРАФИКИ ВЕРТИКАЛЬНЫХ ПЕРЕМЕЩЕНИЙ УЗЛОВ
352	% ГООРИЗОНТАЛЬНЫХ СЕЧЕНИЙ НА ГЛУБИНАХ h
353	X=zeros(2*egor+1);
354	Y=zeros(2*egor+1);
355	for i=1:(2*egor+1)
356	j=Nv(45,i); % ГЛУБИНА h=5 СМ
357	X(i)=x(j);
358	Y(i)=U(j+np);
359	end
360	hPlot=plot(X,Y,'-');grid on;xlabel('X, m');ylabel('Uy, m')
361	set(hPlot,'LineWidth',2);
362	hold on
363
364	X=zeros(2*egor+1);
365	Y=zeros(2*egor+1);
366	for i=1:(2*egor+1)
367	j=Nv(41,i); % ГЛУБИНА h=15 СМ
368	X(i)=x(j);
369	Y(i)=U(j+np);
370	end
371	hPlot=plot(X,Y,'-');grid on;xlabel('X, m');ylabel('Uy, m')
372	set(hPlot,'LineWidth',2);
373	hold on
374
375	X=zeros(2*egor+1);
376	Y=zeros(2*egor+1);
377	for i=1:(2*egor+1)
378	j=Nv(39,i); % ГЛУБИНА h=25 СМ
379	X(i)=x(j);
380	Y(i)=U(j+np);
381	Ux325(i)=U(j);
382	Uy325(i)=Y(i);
383	end
384	hPlot=plot(X,Y,'--');grid on;xlabel('X, m');ylabel('Uy, m')
385	set(hPlot,'LineWidth',2);
386	hold on
387
388	for i=1:(2*egor+1)
389	j=Nv(35,i); % ГЛУБИНА h=45 СМ
390	X(i)=x(j);
391	Y(i)=U(j+np);
392	X345(i)=X(i);
393	Ux305(i)=U(j);
394	Uy305(i)=Y(i);
395	end
396	hPlot=plot(X,Y,'--');grid on;xlabel('X, m');ylabel('Uy, m')
397	set(hPlot,'LineWidth',2);
398	hold on
399
400	for i=1:(2*egor+1)
401	j=Nv(29,i); % ГЛУБИНА h=80 СМ
402	X(i)=x(j);
403	Y(i)=U(j+np);
404	Ux270(i)=U(j);
405	Uy270(i)=Y(i);
406	end
407	hPlot=plot(X,Y,'-.');grid on;xlabel('X, m');ylabel('Uy, m')
408	set(hPlot,'LineWidth',2);
409
410	figure
411
412	% ВЫЧИСЛЕНИЕ КОМПОНЕНТОВ ДЕФОРМАЦИИ В СРЕДИННЫХ ТОЧКАХ
413	% ЭЛЕМЕНТОВ (eps=0.0, eta=0.0)
414
415	epsX0=zeros(ever*egor);
416	epsY0=zeros(ever*egor);
417	gamXY0=zeros(ever*egor);
418
419	for n=1:ever
420	for m=1:egor
421	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
422	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
423	n6=n5-1; n7=n5-2; n8=n4-1;
424
425	dJ=abs((x(n1)-x(n3))*(y(n3)-y(n7)))/4; % якобиан
426	i=m+egor*(n-1);
427	epsX0(i)=(U(n4)-U(n8))/(x(n4)-x(n8));
428	epsY0(i)=(-U(n2+np)+U(n6+np))/(y(n6)-y(n2));
429	gamXY0(i)=-(-U(n2)+U(n6))/(x(n4)-x(n8))+(U(n4+np)-U(n8+np))/(y(n6)-y(n2));
430	end
431	end
432
433	% НОМЕРА ЭЛЕМЕНТОВ ПОСТРОЧНО
434	NEv=zeros(ever,egor);
435	for n=1:ever
436	for m=1:egor
437	i=m+egor*(n-1);
438	NEv(n,m)=i;
439	end
440	end
441
442	% НОМЕР СРЕДИННОГО УЗЛА, ДЛЯ КОТОРОГО ВЫЧИСЛЯЮТСЯ ДЕФОРМАЦИИ % И НАПРЯЖЕНИЯ,
443	% СОВПАДАЕТ С НОМЕРОМ ЭЛЕМЕНТА, РАСПОЛОЖЕННОГО В МАССИВЕ NEv(n,m)
444
445	X=zeros(egor); Y25=zeros(egor); Y200=zeros(egor); Y400=zeros(egor);
446	for m=1:egor
447	n=Nv(1,2*m);
448	X(m)=x(n);
449	i25=NEv(23,m);
450	i200=NEv(20,m);
451	i400=NEv(18,m);
452	Y25(m)=epsX0(i25); % 2,5 cм
453	Y200(m)=epsX0(i200); % 20,0 cм
454	Y400(m)=epsX0(i400); % 40,0 см
455	end
456	hPlot=plot(X,Y25,'-',X,Y200,'--',X,Y400,'-.');...
457	grid on;xlabel('X, m');ylabel('epsX')
458	set(hPlot,'LineWidth',2);
459	figure
460
461	X=zeros(egor); Y25=zeros(egor); Y200=zeros(egor); Y400=zeros(egor);
462	for m=1:egor
463	n=Nv(1,2*m);
464	X(m)=x(n);
465	i25=NEv(23,m);
466	i200=NEv(20,m);
467	i400=NEv(18,m);
468	Y25(m)=epsY0(i25); % 2,5 cм
469	Y200(m)=epsY0(i200); % 20,0 cм
470	Y400(m)=epsY0(i400); % 40,0 см
471	end
472	hPlot=plot(X,Y25,'-',X,Y200,'--',X,Y400,'-.');...
473	grid on;xlabel('X, m');ylabel('epsY')
474	set(hPlot,'LineWidth',2);
475	figure
476
477	% ВЫЧИСЛЕНИЕ КОМПОНЕНТОВ НАПРЯЖЕНИЙ Sigma В СРЕДИННЫХ ТОЧКАХ
478	% ЭЛЕМЕНТОВ (eps=0.0, eta=0.0)
479	SigX=zeros(ever); SigY=zeros(ever);
480	k=0;
481	for n=1:ever
482	for m=1:egor
483	n1=2*m-1+(3*egor+2)*(n-1); n2=n1+1; n3=n1+2;
484	n4=n3+2*egor-m+1; n5=n3+(3*egor+2);
485	n6=n5-1; n7=n5-2; n8=n4-1;
486
487	% ФОРМИРОВАНИЕ МАТРИЦЫ УПРУГОСТИ ДЛЯ СЛУЧАЯ ПЛОСКОЙ % ДЕФОРМАЦИИ
488	if((n>=1)&(n<=7))
489	e1=E(6); nu1=nu(6); % ГРУНТ СУГЛИНОК ЛЕГКИЙ
490	end
491	if((n>=8)&(n<=12))
492	e1=E(5); nu1=nu(5); % ГРАВИЙНО-ПЕСЧАНАЯ СМЕСЬ
493	end
494	if((n>=13)&(n<=16))
495	e1=E(4); nu1=nu(4); % ЩЕБЕНОЧНАЯ СМЕСЬ
496	end
497	if((n>=17)&(n<=20))
498	e1=E(3); nu1=nu(3); % ЩЕБЕНОЧНО-ПЕСЧАНАЯ СМЕСЬ 8%
499	end
500	if((n>=21)&(n<=22))
501	e1=E(2); nu1=nu(2); % АБ КЗ ПОРИСТЫЙ
502	end
503	if(n==23)
504	e1=E(1); nu1=nu(1); % АБ МЗ ПЛОТНЫЙ
505	end
506
507	D=zeros(3,3);
508	D(3,3)=e1/(2*(1+nu1)); D(2,2)=2*D(3,3)*(1-nu1)/(1-2*nu1);
509	D(1,1)=D(2,2); D(1,2)=2*D(3,3)*nu1/(1-2*nu1);
510	D(2,1)=D(1,2);
511
512	X1=x(n1); X2=x(n3); Y2=y(n3); Y4=y(n7);
513	dJ=abs((X1-X2)*(Y2-Y4))/4; % Якобиан
514	G=e1/(2*(1+nu1));
515
516	i=m+egor*(n-1);
517	% В МПа
518	SigmaX0(i)=(D(1,1)*epsX0(i)+D(1,2)*epsY0(i))/100000;
519	SigmaY0(i)=(D(2,1)*epsX0(i)+D(2,2)*epsY0(i))/100000;
520	TauXY0(i)=D(3,3)*gamXY0(i)/100000;
521	if(m==30)
522	k=k+1;
523	SigX(n)=SigmaX0(i);
524	SigY(n)=SigmaY0(i);
525	end
526	end
527	end
528
529	% ПОСТРОЕНИЕ ГРАФИКОВ КОМПОНЕНТ НАПРЯЖЕНИЙ SigmaX0, SigmaY0 И
530	% TauXY0 В СРЕДИННЫХ УЗЛАХ ЭЛЕМЕНТОВ
531
532	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
533	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
534	for m=1:egor
535	n=Nv(1,2*m);
536	X(m)=x(n);
537	i23=NEv(23,m); % h=2,5 см
538	i22=NEv(22,m); % h=7,5 см
539	i20=NEv(20,m); % h=20,0 см
540	i18=NEv(18,m); % h=40,0 см
541	i14=NEv(14,m); % h=87,5 см
542	i12=NEv(12,m); % h=120,0 см
543	Y23(m)=SigmaX0(i23); % 2,5 cм
544	Y22(m)=SigmaX0(i22); % 7,5 cм
545	Y20(m)=SigmaX0(i20); % 20,0 cм
546	Y18(m)=SigmaX0(i18); % 40,0 cм
547	Y14(m)=SigmaX0(i14); % 87,5 cм
548	Y12(m)=SigmaX0(i12); % 120,0 cм
549	end
550	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y18,'-.',X,Y14,'-r');...
551	grid on;xlabel('X, m');ylabel('SigmaX, МПа')
552	set(hPlot,'LineWidth',2);
553	figure
554
555	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
556	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
557	for m=1:egor
558	n=Nv(1,2*m);
559	X(m)=x(n);
560	i23=NEv(23,m); % h=2,5 см
561	i22=NEv(22,m); % h=7,5 см
562	i20=NEv(20,m); % h=20,0 см
563	i18=NEv(18,m); % h=40,0 см
564	i14=NEv(14,m); % h=87,5 см
565	i12=NEv(12,m); % h=120,0 см
566	Y23(m)=SigmaY0(i23); % 2,5 cм
567	Y22(m)=SigmaY0(i22); % 7,5 cм
568	Y20(m)=SigmaY0(i20); % 20,0 cм
569	Y18(m)=SigmaY0(i18); % 40,0 cм
570	Y14(m)=SigmaY0(i14); % 87,5 cм
571	Y12(m)=SigmaY0(i12); % 120,0 cм
572	end
573	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y18,'-.',X,Y14,'-r');...
574	grid on;xlabel('X, m');ylabel('SigmaY, МПа')
575	set(hPlot,'LineWidth',2);
576	figure
577
578	X=zeros(egor); Y23=zeros(egor); Y22=zeros(egor); Y20=zeros(egor);
579	Y18=zeros(egor); Y14=zeros(egor); Y12=zeros(egor);
580	for m=1:egor
581	n=Nv(1,2*m);
582	X(m)=x(n);
583	i23=NEv(23,m); % h=2,5 см
584	i22=NEv(22,m); % h=7,5 см
585	i21=NEv(21,m); % h=12,5 см
586	i20=NEv(20,m); % h=20,0 см
587	i18=NEv(18,m); % h=40,0 см
588	i16=NEv(16,m); % h=61,25 см
589	i14=NEv(14,m); % h=87,5 см
590	i12=NEv(12,m); % h=120,0 см
591	Y23(m)=TauXY0(i23); % 2,5 cм
592	Y22(m)=TauXY0(i22); % 7,5 cм
593	Y21(m)=TauXY0(i21); % 12,5 cм
594	Y20(m)=TauXY0(i20); % 20,0 cм
595	Y18(m)=TauXY0(i18); % 40,0 cм
596	Y16(m)=TauXY0(i16); % 61,25 cм
597	Y14(m)=TauXY0(i14); % 87,5 cм
598	Y12(m)=TauXY0(i12); % 120,0 cм
599	end
600	hPlot=plot(X,Y22,'-',X,Y20,'--',X,Y14,'-.');...
601	grid on;xlabel('X, m');ylabel('TauXY, МПа')
602	set(hPlot,'LineWidth',2);
603
604	toc

Приложение Б

(обязательное)

Исходный код подпрограммы BtDBnds

1	function[BDB]=BtDBnds(psi,eta,D,X1,X2,Y1,Y3);
2	A=(-2/(X1-X2));
3	B1=A*(1-eta)*(2*psi+eta)/4; B2= A*(-psi*(1-eta));
4	B3= A*(1-eta)*(2*psi-eta)/4; B4= A*(1-eta^2)/2;
5	B5= A*(1+eta)*(2*psi+eta)/4; B6= A*(-psi*(1+eta));
6	B7= A*(1+eta)*(2*psi-eta)/4; B8= A*(-(1-eta^2)/2);
7	B=(-2/(Y1-Y3));
8	C1=B*(1-psi)*(psi+2*eta)/4; C2=-B*(1-psi^2)/2;
9	C3=B*(1+psi)*(-psi+2*eta)/4; C4=-B*eta*(1+psi);
10	C5=B*(1+psi)*(psi+2*eta)/4; C6=B*(1-psi^2)/2;
11	C7=B*(-(1-psi)*(psi-2*eta)/4); C8=-B*eta*(1-psi);
12
13
14	B=[B1 B2 B3 B4 B5 B6 B7 B8 0 0 0 0 0 0 0 0;
15	0 0 0 0 0 0 0 0 C1 C2 C3 C4 C5 C6 C7 C8;
16	C1 C2 C3 C4 C5 C6 C7 C8 B1 B2 B3 B4 B5 B6 B7 B8];
17	Bt=B';
18	BtD=Bt*D;
19	BDB=BtD*B;

Библиография

[1]Мартынов Н.Н., Иванов А.П. MATLAB 5.Х. Вычисления, визуализация, программирование. – М.:КУДИЦ-ОБРАЗ, 2000.-336 с.

[2] Коткин Г.Л., Черкасский В.С. Компьютерное моделирование физических процессов с использованием MATLAB:Учебное пособие/Новосиб. ун-т. Новосибирск, 2001. 173 с.

[3] Потемкин В.Г. Система инженерных и научных расчетов MATLAB 5.Х. В 2-х томах. –М.:ДИАЛОГ-МИФИ, 1999. 670 с.

[4] Тынкевич М.А. Численные методы.-Кемерово.:КузГТУ.1977.-122 с.

[5] Сегерлинд Л. Применение метода конечных элементов. М: Мир, 1979.– 392 с.

[6] СТ РК 1293-2004 Методы определения модуля упругости дорожных одежд нежесткого типа и их классификация.