«MATLAB В ИНЖЕНЕРНЫХ И НАУЧНЫХ РАСЧЕТАХ Одесса Астропринт 2003 ББК Д УДК 539.3:681.3 Монография посвящена иллюстрации возможностей одной из самых эффективных систем компьютерной математики MATLAB в решении ряда научных и ...»

А. Ф. Дащенко, В. Х. Кириллов, Л. В. Коломиец, В. Ф. Оробей

MATLAB

В ИНЖЕНЕРНЫХ И НАУЧНЫХ РАСЧЕТАХ

Одесса «Астропринт» 2003

ББК

Д

УДК 539.3:681.3

Монография посвящена иллюстрации возможностей одной из самых

эффективных систем компьютерной математики MATLAB в решении ряда

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

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

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

Рецензенты: заведующий кафедрой сопротивления материалов и строительной механики Одесского государственного морского университета, доктор технических наук, профессор В.А. Гришин;

заведующий кафедрой высшей математики Одесского национального политехнического университета, доктор технических наук, профессор А. В. Усов Рекомендовано к изданию ученым советом института машиностроения Одесского национального политехнического университета от 17 декабря 2002 г.

Без объявл.

Д © А.Ф. Дащенко, В. Х. Кириллов, ISBN 966-549- Л. В. Коломиец, В. Ф. Оробей,

ВВЕДЕНИЕ

Данная книга посвящена иллюстрации возможностей одной из самых эффективных систем компьютерного программного обеспечения – пакета универсальных интегрированных программ MATLAB. Любознательному читателю предлагается ознакомиться в первом приближении с основами языка программирования и комплексной визуализации результатов решения ряда научных и инженерных задач. Рассматриваются такие проблемы как табулирование функций, решение нелинейных уравнений, поиск оптимальных решений, решение задач Коши, численное интегрирование и другие задачи, традиционно включаемые в курс численных методов. Алгоритм этих задач хорошо известен и разработчики системы MATLAB (фирма Math Works, Inc., U.S.A.) учли опыт численного решения и программирования задач вычислительной математики за все время существования вычислительной техники.

Поэтому в системе MATLAB по каждой проблеме имеется несколько программ (иногда их более 10), предназначенных для ее решения в зависимости от особенностей данной задачи. Кроме чисто научных задач, средствами MATLAB могут быть успешно решены и довольно сложные инженерные проблемы, такие, как поиск спектра частот собственных колебаний и критических сил потери устойчивости стержневых, пластинчатых и оболочечных систем, решение краевых задач для упругих систем и задач сейсмостойкости сооружений и др. Численные результаты таких задач должны сопровождаться соответствующими эпюрами и формами, т. е. визуализацией расчетов.

MATLAB обеспечивает решение и этой проблемы, что очень привлекательно для учебного процесса, где студенты тратят достаточно много времени на графическую часть курсовых и дипломных проектов. Кроме MATLAB существуют и другие, довольно мощные среды программирования и визуализации, такие как Visual Digital Fortran, Delphi, Visual C++ и т. п. Однако, по мнению авторов, в системе MATLAB получаются наиболее простые и в то же время эффективные программы. Читателю предлагается самому убедиться в этом, сравнив приведенные в данной книге программы с программами других сред. Вероятно, нельзя создать в других средах программирования программы более простые, чем в системе MATLAB.

В качестве инженерных задач в пособии рассматриваются задачи строительной механики – науки о расчетах сооружений на статическую, динамическую нагрузки и устойчивость. Для решения задач строительной механики разработано множество методов – методы сил и перемещений, метод конечных элементов, метод конечных разностей, метод R-функций, метод граничных элементов и др. Здесь для иллюстрации возможностей MATLAB выбран один из наиболее эффективных методов – метод граничных элементов (МГЭ), позволяющий существенно упростить алгоритм решаемых задач [2]. Объединение возможностей MATLAB и МГЭ позволяет на качественно более высоком уровне представлять решения инженерных задач, что существенно повышает научное и практическое содержание учебного процесса. В этой связи, кроме основ языка программирования MATLAB, добавлены краткие сведения о МГЭ.

Слово MATLAB состоит из начальных букв слов МАТrix LABoratory – матричная лаборатория. Название системы полностью отражает ее суть. Это действительно матричная лаборатория, где начальным кирпичиком является не простая переменная или константа, а матрица и ее частные случаи - вектор-строка, вектор-столбец.

Систему MATLAB разработал Молер (С.В. Moler) в 70-х г. г. ХХ века, которая использовалась на больших ЭВМ. В начале 80-х г. г. Джон Литл (John Little) из фирмы Math Works, Inс. Модернизировал эту систему для персональных компьютеров типа IBM PC, VAX и Macintosh. Далее к расширению системы были привлечены крупнейшие ученые и научные школы в математике, программировании и естествознании. Это позволило MATLAB стать признанным лидером в решении различных проблем науки и техники среди других подобных систем. Этому способствовало создание языка программирования, который вобрал в себя преимущества традиционных языков (Fortran, Pascal, Basic, C++) и достаточно мощных средств визуализации и моделирования. Более подробно о преимуществах и возможностях системы MATLAB можно узнать в специализированных изданиях [1].

Данная работа содержит 4 главы. В первой главе представлены основные элементы языка программирования MATLAB и средств визуализации результатов расчетов. Вторая глава посвящена решению в системе MATLAB задач вычислительной математики и варианты заданий. Третья глава содержит подробные алгоритмы и программы решения задач статики, динамики и устойчивости упругих систем методом граничных элементов. Приведены так-же варианты заданий на самостоятельную работу. Четвертая глава посвящена пакетам расширений Symbolic Math., Optimization Toolbox u Simulink, позволяющих выполнять символьные вычисления, решать задачи оптимизации и моделировать динамические системы и устройства. Иллюстрация возможностей этих пакетов выполнена на задачах теоретической механики.

ГЛАВА 1. Основные элементы языка программирования и визуализации 1.1. Алфавит языка программирования В MATLAB, как и в других системах, используются все буквы латинского алфавита от А до Z и арабские цифры от 0 до 9. Как и в С++, большие и малые буквы это разные переменные и константы. Кроме букв латинского алфавита используются все специальные символы клавиатуры компьютера.

1.2. Арифметические и логические операторы Число арифметических операторов в MATLAB значительно расширено и включает в себя матричные и арифметические операции. В таблице 1. приводится список арифметических операторов.

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

Операторы отношения служат для сравнения двух величин, векторов или матриц, все операторы отношения имеют две сравниваемые величины и записываются, как показано в таблице 1.2.

Не равно Меньше Больше Меньше или равно Больше или равно Данные операторы выполняют поэлементное сравнение векторов или матриц одинакового размера и логическое выражение принимает значение (True), если элементы идентичны, и значение 0 (False) в противном случае.

Логические операторы служат для реализации поэлементных логических операций над элементами одинаковых по размеру массивов согласно таблице 1.3.

Логическое И Логическое ИЛИ Логическое НЕ Верно, если все элементы вектора равны нулю any (any (a)) Верно, если все элементы вектора не равны нулю all (all (a)) 1.3. Элементарные функции Набор элементарных функций представим их описанием. В тригонометрических функциях углы измеряются в радианах.

arсch x – арккосинус гиперболический acosh(x) arсcth x – арккотангенс гиперболический acoth(x) arссosech x – арккосеканс гиперболический acsch(x) arсsech x – арксеканс гиперболический asech(x) Следует помнить, что все элементарные функции должны записываться в программах малыми буквами. Существуют также специальные математические функции, на которых мы не будем останавливаться. Их описание см. в [1].

1.4. Понятие о файлах-сценариях и файлах-функциях При загрузке системы MATLAB на мониторе появляется основное окно системы, в котором можно выделить окно команд (Command Window).

Система готова к проведению вычислений и созданию программ в командном режиме. Для этого можно на языке MATLAB записывать программы.

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

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

Щелкнув по пиктограмме стрелкой мышки, получаем окно М-файла, на котором можно записывать, редактировать и отлаживать любые программы решения научных и инженерных задач. Данный М-файл по умолчанию имеет название Untitled (Безымянный). Чтобы дать ему имя, необходимо в меню этого окна File выполнить команду Save as и в другом окне указать папку и имя этого файла. После указания имени и сохранения М-файла он готов для выполнения записанной программы. Для этого необходимо щелкнуть мышкой по пиктограмме Выполнить. Она выполнена в виде страницы со стрелкой, направленной вниз. Результат выполнения программы или сообщения об ошибках появится в окне команд. Описанный процесс называется созданием М-файла сценария сессии. Файл-сценарий, именуемый также Scriptфайлом, имеет весьма простую структуру:

% Основной комментарий, если необходимо.

% Дополнительный комментарий, если необходимо.

Тело программы с любыми выражениями.

Важными являются следующие свойства файлов-сценариев:

1. Они не имеют входных и выходных аргументов.

2. Работают с данными из рабочей области.

3. В процессе выполнения не компилируются.

4. Представляют собой последовательность операций, аналогичную той, что используется в сессии.

Кроме М-файла сценария, в MATLAB существует М-файл функция.

Отличие М-файла функции от сценария состоит в том, что он является аналогом подпрограммы типа function в языке Pascal.

Структура М-файла функции с одним выходным параметром имеет вид:

function var = f _ name (Список параметров) % Основной комментарий, если необходимо.

% Дополнительный комментарий, если необходимо.

Тело программы с любыми выражениями.

var = выражение М-файл функция обладает такими свойствами:

1. Он начинается с ключевого слова function, после которого указывается имя переменной var – выходного параметра, имя самой функции f _ name и список ее входных параметров, отделенных запятой.

Внимание: Имя М-файла функции должно совпадать с самой f _ name (именем самой функции). MATLAB автоматически присваивает данное имя при выполнении команды Save as.

2. Результат выполнения М-файла функции присваивается имени функции, которое может использоваться в математических выражениях подобно функциям sin(x), log(x) и т. п.

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

4. Файл-функция является самостоятельным программным модулем, который связан с другими модулями и головной программой через входные и выходные параметры.

5. При обнаружении файла-функции он компилируется и затем исполняется.

Ниже в главах 2, 3 показаны примеры создания и использования в практических расчетах М-файлов сценариев и функций.

1.5. Основы программирования 1.5.1. Оператор присваивания Программирование, т. е. создание определенного набора команд, в системе MATLAB является средством ее расширения и использовании в решении специфических проблем. Отдельные вопросы программирования изложены выше, здесь рассмотрим правила, дополняющие синтаксис языка MATLAB.

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

Для задания переменным определенных значений используется оператор присваивания, вводимый знаком равенства = Имя _ переменной = Выражение ;

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

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

1.5.3. Ввод и вывод данных В языке MATLAB нет явных операторов ввода вывода данных. Эта проблема решается для ввода данных оператором присваивания и использованием системных констант. Вывод данных осуществляется еще проще.

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

При вычислениях в MATLAB используется режим двойной точности.

Однако, при выводе результатов, по умолчанию выдаются числа с 4 цифрами после десятичной точки в действительной форме. Чтобы изменить данную форму вывода, необходимо в программе перед выводимой величиной использовать команду format name, где name имя формата. Для числовых данных name может быть следующим сообщением:

short короткое представление в фиксированном формате (5 знаков);

short е короткое представление в экспоненциальной форме (5 знаков мантиссы и 3 знака порядка);

long – длинное представление в фиксированном формате (15 знаков);

long е – длинное представление в экспоненциальной форме (15 знаков мантиссы и 3 знака порядка).

В качестве примера рассмотрим вывод вектора, содержащий 2 числа:

format name x = [5/3 1.2783 e 7] В различных форматах вывод вектора х будет иметь следующий вид:

shor t shor t e long long e Задание формата сказывается только на форме вывода чисел. Вычисления же происходят в режиме двойной точности, а ввод чисел осуществляется в любом удобном виде.

1.5.5. Формирование векторов и матриц Описанные правила вычислений распространяются и на более сложные вычисления, которые при использовании обычных языков программирования (типа Pascal, Fortran, C++ и др.) требуют составления специальных программ.

MATLAB специально предназначен для проведения сложных вычислений с векторами и матрицами. При этом по умолчанию предполагается, что каждая переменная – это вектор или матрица. Например, если задано х = 1, то это значит, что х – это вектор с одним элементом, равным 1. Если надо задать вектор из трех элементов, то их значения надо перечислить в квадратных скобках, разделяя пробелами.

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

Задание матрицы требует указания несколько строк. Для разграничения строк используется символ ; (точка с запятой).

Для указания отдельного элемента вектора или матрицы используются выражения вида V(i) или T(i, j). Например:

Если элементу Т(i, j) нужно присвоить новое значение х, то используют оператор присваивания Выражение Т(i) с одним индексом дает доступ к элементам матрицы, развернутым в один столбец. Такая матрица образуется из исходной, если подряд выписать ее столбцы. Например:

Наряду с операциями над отдельными элементами матриц и векторов MATLAB позволяет производить арифметические операции сразу над всеми элементами. Для этого перед знаком операции ставится точка.

Имеются также ряд особых функций для задания векторов и матриц.

Отметим функции ones и zeros. Эти функции служат для создания одномерных и многомерных массивов. Функция ones создает массив с единичными элементами Функция zeros создает массив с нулевыми элементами В краевых задачах поиска частот собственных колебаний и критических сил потери устойчивости упругих систем функция zeros находит широкое применение, что существенно упрощает тексты соответствующих программ (см. главу 3).

1.5.6. Типы данных Структура типов данных MATLAB представлена схемой по рис. 1.1.

int16, uint Типы данных array и numerric являются виртуальными (“кажущимися”), т. к. к ним нельзя отнести какие-либо переменные. Они служат для определения и комплектования некоторых типов данных. Таким образом, в MATLAB определены типы данных, представляющие собой многомерные массивы:

1. single – числовые массивы с числами одинарной точности;

2. double – числовые массивы с числами удвоенной точности;

3. chаr – строчные массивы с элементами – символами;

4. sparsе – разреженные матрицы с элементами – числами двойной точности;

5. cell – массивы ячеек;

6. structure – массивы структур с полями;

7. function – handle – дескрипторы функций;

8. int8, …, uint32 – массивы 8-, 16-, 32-разрядных целых чисел со знаками и без знаков.

1.5.7. Оператор двоеточие :

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

х = Начальное _ значение : Шаг : Конечное _ значение ;

Эта конструкция создает возрастающую последовательность чисел, которая начинается с начального значения, изменяется на заданный шаг и завершается конечным значением. Если шаг не задан, то он принимает значение 1. Если конечное значение указано меньшим, чем начальное значение, то выдается сообщение об ошибке. Примеры:

и т.д.

1.5.8. Оператор разветвления if Условный оператор if в MATLAB записывается в общем виде так:

if Логическое условие Оператор 1 elseif Логическое условие Оператор else Оператор 3 end ;

Эта конструкция имеет несколько частных вариантов:

if Логическое условие Оператор 1 end ;

if Логическое условие Оператор 1 else Оператор 2 end ;

Логическое условие записывается в виде:

Выражение 1 Оператор отношения Выражение В качестве операторов отношения используются операторы: = =,, < =, > =, =. Если логическое условие принимает значение 1(true – истина), то выполняются соответствующие операторы. Если логическое условие принимает значение 0(false – ложь), то операторы, следующие за логическим условием, не выполняются. Оператор end указывает на конец условного оператора if. В понятие Оператор 1 входят один или несколько операторов. В последнем случае они разделяются символами, (запятой) или ; (точкой с запятой).

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

1.5.9. Операторы циклов В MATLAB существует 3 типа операторов цикла. С оператором :

(двоеточие) мы познакомились в п.1.5.7. Следующий оператор for … end используется для организации цикла с фиксированным числом повторений. Он имеет вид:

for var = Выражение Операторы end ;

Здесь var – счетчик цикла – любая переменная, обычно это i, j, k, l, m и т. д. Выражение записывается в виде s : d : e, где s – начальное значение счетчика цикла var, d – шаг изменения и е – конечное значение var. Возможна и запись в виде s : e, тогда d = 1. Список операторов завершается ключевым словом end. Оператор continue передает управление в следующую итерацию цикла, пропуская операции, которые записаны за ним. Оператор break используется для досрочного прерывания цикла. Возможны вложенные циклы >> for i = 1 : 3 for j = 1 : 3 a (i, j) = i j ; end ; end ;

В результате выполнения этого цикла формируется матрица а Циклы типа while … end выполняются до тех пор, пока выполняется заданное условие. Оператор записывается в виде:

while Логическое условие Операторы end ;

1.5.10. Сообщения об ошибках и исправление ошибок Система MATLAB контролирует правильность написания программ и, при наличии ошибок, выдает соответствующее сообщение в окне команд.

При этом указывается номер строки, где допущена ошибка, и характер ошибки. После уяснения сути ошибки ее необходимо исправить в тексте программы, запомнить М-файл командой Save и снова выполнить программу. Перед этим желательно очистить окно команд от сообщения об ошибках (чтобы не загромождать полученную картинку) с помощью команды Clear Command Windows (Очистить окно команд) в меню Edit.

1.5.11. Вычисление определителя квадратной матрицы Для вычисления определителя квадратной матрицы используется функция det(a). Если матрица а содержит только целые числа, то результат – тоже целое число. Определитель вычисляется на основе треугольного разложения методом исключения Гаусса. Пример:

>> det(a) Данная функция широко используется в задачах поиска спектров частот собственных колебаний и критических сил потери устойчивости упругих систем (см. главу 3).

1.5.12. Обращение матриц Обращение матриц – одна из наиболее распространенных операций задач строительной механики и других наук. Обратной называют матрицу, получаемую в результате деления единичной матрицы Е на исходную матрицу х, т. е. х-1 = Е/х. Эту процедуру выполняет функция inv(x), которая вычисляет элементы обратной матрицы для исходной квадратной матрицы х. Выдается предупреждающее сообщение, если матрица х плохо масштабирована или близка к вырожденной. На практике вычисление обратной матрицы не так уж необходимо. Чаще обращение применяют для решения систем линейных алгебраических уравнений вида ах = b. Один из путей решения этой системы – х = inv(a) b, хотя лучше использовать метод исключения Гаусса без формирования обратной матрицы, например х = a/b или х = b/a.

1.6. Основы графической визуализации вычислений Во многих областях науки и техники численное решение задач недостаточно для анализа результатов. Необходима еще графическая интерпретация в виде эпюр параметров напряженно-деформированного состояния элементов упругих систем, формы колебаний и потери устойчивости, поведение решений на заданном интервале и т. п. MATLAB позволяет решать эти задачи достаточно простыми процедурами. Вначале необходимо задать интервал изменения аргумента х от начального значения х0 до конечного хк с шагом х, что осуществляется оператором двоеточие : х0 : х : хк. Далее используется команда построения графика какой-либо функции у = f(x), которая носит имя plot.

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

Если число точек достаточно велико, то полученная кривая воспринимается как вполне истинный график функции у = f(x), при 10 – 20 точках получается ломаная кривая.

1.6.1. Построение графиков отрезками прямых Для построения графика функции у = f(x) необходимо задать совокупность точек х и у. Для аргумента х это выполняется оператором двоеточие, для у – надлежащим программированием выражения для функции, т. е. необходимо применить знаки арифметических операций над массивами:

Для отображения таких функций используется декартовая прямоугольная система координат. Команда построения графика функции у = f(x) plot имеет ряд параметров, которые рассмотрим ниже.

рlot(х, у) – строит график функции у = f(x), координаты точек (х, у) которой берутся из векторов одинакового размера х, у.

рlot(х, у, s) – аналогична команде рlot(х, у), но тип линии графика можно задавать с помощью строковой константы s, значения которой представлены в таблице 1.5.

s = y – желтый Таким образом, с помощью строковой константы s можно менять цвет линии, представлять узловые точки различными отметками и менять тип линии графика. Рассмотрим пример построения графиков трех функций y1 = sin x; y2 = sin2 x; y3 = sin3 x с различным стилем:

Здесь график функции у1 строится сплошной фиолетовой линией, график у строится штрихпунктирной линией с точками в виде знака “+” красного цвета, график у3 строится штриховой линией с кружками черного цвета.

1.6.2. Графики дискретных отсчетов функции В сопротивлении материалов и строительной механике принято графики параметров напряженно-деформированного состояния закрашивать (штриховать). Для построения подобных графиков используется команда stem( … ).

stem(x, y) – строит график с закрашенными частями между нулевой линией и кривой у = f(x).

stem( …, LINESPEC ) – аналогична команде stem(x, y), но имеет спецификацию линий LINESPEC, подобную строковой константе s команды plot(х, у, s).

Stem( …, filled ) – строит график функции у = f(x) с закрашенными маркерами.

1.6.3. Создание массивов данных для трехмерной графики Трехмерные поверхности описываются функцией двух переменных z = f(x, y). Построение трехмерных графиков требует определение для х и у двухмерных массивов – матриц. Для создания таких массивов служит функция meshgrid, которая записывается следующим образом:

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

Строки выходного массива Х являются копиями вектора х, а столбцы Y – копиями вектора у.

Приведем еще пример применения функции meshgrid:

Такой вызов функции создает опорную плоскость для построения трехмерной поверхности при изменении х и у от – 1 до 1 с шагом 0.1.

1.6.4. Построение графиков поверхностей Для построения графиков функции z = f(x, y) используется команды plot3 ( … ), которая является аналогом команды plot ( … ). Она строит аксонометрическое изображение трехмерной поверхности и имеет следующие формы:

рlot3 (х, у, z) – строит массив точек, представленных векторами х, у, z и соединяет их отрезками прямых.

рlot3 (X, Y, Z), где X, Y, Z – три матрицы одинакового размера, строит точки с координатами X(i, : ), Y(i, : ) и Z(i, : ) и соединяет их отрезками прямых. Пример построения графика трехмерной поверхности Z = x2 + y2:

рlot3 (X, Y, Z, S) – обеспечивает построение графика поверхности, но со спецификацией стиля линий и точек, соответствующей спецификации команды plot.

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

grid on – добавляет сетку к текущему графику;

grid off – отключает сетку;

grid – последовательно производит включение и отключение сетки.

1.6.6. Представление нескольких графиков в одном окне Иногда необходимо в одном окне поместить несколько графиков без наложения их друг на друга. Особенно удобно такое представление при построении эпюр напряженно-деформированного состояния элементов стержневых и пластинчатых систем. Для этого служит команда subplot, которую необходимо записать перед командой plot.

subplot (m, n, p) – разбивает графическое окно на m n подокон, при этом m – число подокон по горизонтали, п – число окон по вертикали, р – номер подокна, в которое будет выводиться текущий график.

1.6.7. Ввод текста на график с помощью мыши Для маркировки графиков можно ввести любой текст с помощью мыши командой gtext. Команда помещается после команды plot.

gtext ( string ) – выводит на график перемещаемый мышкой маркер в виде крестика. Установив маркер в нужное место и щелкнув кнопкой мыши, получим текст на графике.

1.6.8. Управление свойствами осей графиков Если не задавать масштаб графика, то он строится командой plot автоматически. Не всегда этот масштаб удовлетворяет пользователя. Команда axis позволяет установить любой масштаб.

axis ([ x min x max y min y max ]) – устанавливает нужный диапазон координат графика по осям х и у.

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

ГЛАВА 2. MATLAB в задачах вычислительной Система MATLAB, как следует из первой главы книги, обладает большими возможностями программирования и комплексной визуализации результатов инженерных расчетов и научных исследований. В этой связи покажем применение богатых возможностей MATLAB в решении задач вычислительной математики. Развитие многих наук привело исследователей к необходимости численного решения различных проблем, т. е. применению численных методов.

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

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

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

Рассмотрим два случая табулирования функции:

1. С постоянным шагом изменения аргументов.

2. С произвольным набором значений аргумента.

Алгоритм реализуется путем организации какого-либо цикла.

Пример 1. Вычислить при R = 4.28 10-2; = 2.87;

хi изменяется с шагом х = 2; хп = 2; хк = 10.

Введем обозначение la = 2.87.

Протокол программы:

% Задается начальное значение х, шаг dx и конечное значение х % Для вывода значения у в конце строки символ ; не ставится!

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

Результаты вычислений:

ans = Пример 2. Вычислить и вывести на экран значения функции при х1 = 12.8; х2 = 23.4; х3 = 27.2; х4 = 17.8; х5 = 16.3; х6 = 14.9; а = 1.35; b = 0.98.

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

Протокол программы:

а = 1.35; b = 0.98; х(1) = 12.8; х(2) = 23.4; х(3) = 27.2; х(4) = 17.8; х(5) = 16.3;

х(6) = 14.9;

% В конце строки вычисления функции у символ ; не ставится.

Данные вычислений можно вывести в виде таблицы, если использовать запись [x; y] без точки с запятой или [x y].

Составить программу вычисления значений функции уi для значений аргумента хi. Данные взять из таблицы 2.1.

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

Пусть задана система п линейных алгебраических уравнений с п неизвестными:

Система уравнений (2.1) в матричной форме представляется следующим образом:

где А – квадратная матрица коэффициентов, размером п п строк и столбцов;

Х – вектор-столбец неизвестных;

В – вектор-столбец правых частей.

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

Недостатком метода является накапливание погрешностей в процессе округления, поэтому метод Гаусса без выбора главных элементов используется обычно для решения сравнительно небольших (п 100) систем уравнений с плотно заполненной матрицей и не близким к нулю определителем. Если матрица А сильно разрежена, а ее определитель не близок к нулю, то метод Гаусса пригоден для решения больших систем уравнений. В MATLAB имеется обширный арсенал методов решения систем уравнений (2.2) методом исключения Гаусса. Для этого применяются следующие операторы inv(A) обращение матрицы А.

Выражения дают решения ряда систем линейных уравнений АХ = В, где А – матрица размером m n, В – матрица размером п к. Более сложные случаи решения систем уравнений (2.2) с плохо обусловленной матрицей А освещены в работе [1].

Решить систему 4-х линейных уравнений:

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

Х1 = b/а Х3 = b* inv(a) Результаты решения Варианты заданий. Решить систему линейных алгебраических уравнений с помощью 4-х операторов. Данные взять из таблицы 2.2.

2.3. Аппроксимация функций Одним из распространенных и практически важных случаев связи между аргументом и функцией является задание этой связи в виде некоторой таблицы {xi ; yi}, например, экспериментальные данные. На практике часто приходится использовать табличные данные для приближенного вычисления у при любом значении аргумента х (из некоторой области). Этой цели служит задача о приближении (аппроксимации) функций: данную функцию f(x) требуется приближенно заменить некоторой функцией g(х) так, чтобы отклонение g(х) от f(x) в заданной области было наименьшим. Функция g(х) при этом называется аппроксимирующей. Если приближение строится на заданном дискретном множестве точек {xi}, то аппроксимация называется точечной. К ней относятся интерполирование, среднеквадратичное приближение и др. При построении приближения на непрерывном множестве точек (например, на отрезке [a, b]) аппроксимация называется непрерывной или интегральной. MATLAB имеет мощные средства точечной и непрерывной аппроксимации с визуализацией результата. Рассмотрим наиболее важную точечную аппроксимацию (обработка экспериментальных данных).

Пример 4. Используя линейную и полиномиальную аппроксимации, получить эмпирические формулы для функции у=f(x), заданной в табличном виде:

Оценить погрешность эмпирических формул.

Протокол программы. В окне команд набираются значения xi и yi. Далее выполняется команда построения графика только узловых точек.

>> x = [0.75, 1.50, 2.25, 3.00, 3.75] ;

>> y = [2.50, 1.20, 1.12, 2.25, 4.28] ;

Появляется окно с символами 0 на месте узловых точек (рис. 2.1).

Внимание. Следует помнить, что при полиномиальной аппроксимации максимальная степень полиноме на 1 меньше числа экспериментальных точек!

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

- показать уравнение аппроксимирующей функции у = g(х);

- выбрать метод подбора: сплайн интерполяции В нашей задаче выбираем линейную и полиномиальную аппроксимации.

В окне графика появляются соответствующие графики разноцветом и формулы аппроксимирующих функций (рис. 2.2).

Чтобы узнать погрешность аппроксимации, надо отметить птичкой параметр График остатка в окне Основной Монтаж, и Показать норму остатков. График погрешностей с нормами можно вынести в отдельное окно, или вместе с графиком и аппроксимирующих функций – суб-график. Норма погрешностей указывает на статистическую оценку среднеквадратической погрешности. Чем она меньше, тем точнее полученная аппроксимирующая функция у = g(х). В нашем примере:

Linear : norm of residuals (норма погрешности) = 2. Quadratic : norm of residuals = 0. Cubic : norm of residuals = 0. 4 th degree : norm of residuals = 9.6305e-015.

График погрешностей можно выводить в виде диаграмм (зоны), линий (линии) или отдельных точек (фрагменты). Сам график погрешностей представляет собой зависимость разности g(х) f(x) в условных точках, соединенных прямыми линиями.

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

Варианты заданий. Получить эмпирические формулы и оценить их погрешность для функции у = f(x), заданной таблично. Данные взять из таблицы 2.3.

10.

11.

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

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

Постановка задачи:

Требуется найти функцию у = у(х), удовлетворяющую уравнению и принимающую при х = х0 заданное значение у0:

При этом решение необходимо получить в интервале х0 х хк. Из теории дифференциальных уравнений известно, что решение у(х) задачи Коши (2.3), (2.4) существует, единственно и является гладкой функцией, если правая часть F(x, y) удовлетворяет некоторым условиям гладкости. Численное решение задачи Коши методом Рунге-Кутта 4-го порядка заключается в следующем. На заданном интервале [х0, хк] выбираются узловые точки. Значение решения в нулевой точке известно у(х0) = у0. В следующей точке у(х1) определяется по формуле здесь т. е. данный вариант метода Рунге-Кутта требует на каждом шаге четырехкратного вычисления правой части уравнения (2.3). Этот алгоритм реализован в программе ode45. Кроме этой программы MATLAB располагает обширным набором аналогичных программ, позволяющих успешно решать обыкновенные дифференциальные уравнения.

Пример 5. Решить задачу Коши Точное решение имеет вид Выполним решение данной задачи с помощью программы ode45. Вначале в М-файл записываем правую часть уравнения (2.7), сам М-файл оформляется как файл – функция, даем ему имя F:

function dydx = F(x, y) dydx = zeros(1,1);

dydx(1) = 2*(x^2+y(1));

Для численного решения задачи Коши в окне команд набираются следующие операторы.

Протокол программы.

% Дескриптор @ обеспечивает связь с файлом – функцией правой части % [0 1] – интервал на котором необходимо получить решение % [1] – начальное значение решения >> % Построение графика численного решения задачи Коши (2.7) % Команда позволяет с помощью мышки нанести на график надпись у(х) >> % Последняя команда выводит таблицу численного решения задачи.

Результаты решения. График решения задачи Коши (2.7) показан на рис.

2.3. Численное решение представлено в таблице 2.4, где приведены только отдельные узловые точки. В программе ode45 по умолчанию интервал разбивается на 40 точек с шагом h = 1/40 = 0.025.

Как следует из таблицы 2.4 численное решение программой ode45 является точным.

Варианты заданий. Построить график и вывести в виде таблицы решение задачи Коши на интервале [0; 1] методом Рунге-Кутта 4-го порядка. Данные взять из таблицы 2.5.

2.5. Приближенное вычисление определенных интегралов К вычислениям определенных интегралов сводятся многие практические задачи физики, химии, экологии, механики и других естественных наук. На практике формулой Ньютона-Лейбница не всегда удается воспользоваться. В этом случае используются методы численного интегрирования. Они основаны на следующих соображениях: с геометрической точки зрения определенный интеграл представляет собой площадь криволинейной трапеции. Идея методов численного интегрирования сводится к разбиению интервала [a; b] на множество меньших интервалов и нахождению искомой площади как совокупности элементарных площадей, полученных на каждом частичном промежутке разбиения. В зависимости от использованной аппроксимации получаются различные формулы численного интегрирования, имеющие различную точность.

Рассмотрим методы трапеций и Симпсона (парабол).

Метод трапеций.

Здесь используется линейная аппроксимация, т. е. график функции y = f(x) представляется в виде ломаной, соединяющей точки yi. Формула трапеций В MATLAB данную формулу реализует программа trapz (x,y).

Метод Симпсона Если подынтегральную функцию заменить параболой, то формула Симпсона с постоянным шагом интегрирования предстанет в виде В MATLAB формула Симпсона реализуется программой quad. Подынтегральная функция может задаваться с помощью дескриптора @, тогда она программируется в файле – функции, или с помощью апострофов, тогда она записывается в самой программе quad. Точность вычисления интегралов по умолчанию принята равной 110-6.

Пример 6. Вычислить и вывести на печать по методам трапеций и Симпсона значения интеграла Протокол программы метода трапеций Результат вычислений Протокол программы метода Симпсона Результат вычислений ans = Точное значение интеграла равно 0.785398163.

Как видно из примера 6 полученные результаты являются практически точными, а сами протоколы весьма просты.

Варианты заданий. Вычислить и вывести на печать значения определенного интеграла методами трапеций и Симпсона. Данные взять из таблицы 2.6.

2.6. Численное решение нелинейных уравнений Задача нахождения корней нелинейных уравнений встречается в различных областях научно-технических исследований. Проблема формулируется следующим образом. Пусть задана непрерывная функция f(x) и требуется найти корень уравнения Будем предполагать, что имеется интервал изменения х [a; b], на котором необходимо исследовать функцию f(x) и найти значение х0, при котором f(x0) равно или весьма мало отличается от нуля.

Данная задача в системе MATLAB может быть решена следующим образом. Вначале необходимо построить график функции f(x) на заданном интервале и убедиться в существовании корня или нескольких корней. Затем применить программы поиска корней. Если существует один корень и график f(x) пересекает ось ох, то можно применить программу fzero. Если f(x) имеет больше одного корня и может касаться и пересекать ось ох, то следует применить более мощную программу fsolve из пакета Optimization Toolbox, которая решает задачу методом наименьших квадратов. Программа fzero использует известные численные методы: деление отрезка пополам, секущей и обратной квадратичной интерполяции.

Пример 7. Найти корень нелинейного уравнения 10х + 2х – 100 = 0 на интервале [1.0; 2.0].

Протокол программы >> % Строим график заданной функции >> x = 1.0 : 0.001 : 2.0; y = 10.0.^x + 2.0*x – 100.0;

Появляется окно с графиком функции 10х + 2х – 100 (см. рис. 2.4), из которого следует, что корень функции на заданном интервале существует. Для точного определения корня применяем fzero и fsolve.

>> X1 = fzero ( (10.0.^x + 2.0*x – 100.0), [1.0 2.0]) Результат решения >> X2 = fsolve ( (10.0.^x + 2.0*x – 100.0), 1.0 : 2.0) Результат решения Варианты заданий. Построить график и найти корень нелинейного уравнения. Данные взять из таблицы 2.7.

2.7. Численное решение оптимизационных задач Под оптимизацией понимают процесс выбора наилучшего варианта из всех возможных. С точки зрения инженерных расчетов методы оптимизации позволяют выбрать наилучший вариант конструкции, наилучшее распределение ресурсов, минимальный урон природной среде и т. п. В процессе решения задачи оптимизации необходимо найти оптимальные значения некоторых параметров, их называют проектными параметрами. Выбор оптимального решения проводится с помощью некоторой функции, называемой целевой функцией. Целевую функцию можно записать в виде где х1, х2, …, хп – проектные параметры.

Можно выделить 2 типа задач оптимизации – безусловные и условные.

Безусловная задача оптимизации состоит в отыскании максимума или минимума функции (2.10) от п действительных переменных и определении соответствующих значений аргументов на некотором множестве G n-мерного пространства. Обычно рассматриваются задачи минимизации; к ним легко сводятся и задачи на поиск максимума путем замены знака целевой функции на противоположный. Условные задачи оптимизации – это такие, при формулировке которых задаются некоторые условия (ограничения) на множестве G. Здесь рассмотрим только безусловные задачи оптимизации.

Поиск минимума функции одной переменной.

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

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

Протокол программы >> рlot(x, y) ; grid on Появляется окно с графиком этой функции (рис. 2.5), где отмечаем наличие минимума.

Далее, для точного определения координаты и значения минимума привлекаем программу fminbnd.

>> [x, y] = fminbnd ( (24.0 – 2* x/3 + x.^2/30), 5.0, 20.0) Результат поиска Варианты заданий. Найти и вывести на печать координату и минимальное значение функции f(x) на [a; b]. Данные взять из таблицы 2.8.

Поиск минимума функций нескольких переменных Данная задача значительно сложнее первой. Рассмотрим ее решение на примере функции двух переменных. Алгоритм может быть распространен на функции большего числа переменных. Для минимизации функций нескольких переменных MATLAB использует симплекс – метод Нелдера-Мида. Данный метод является одним из лучших методов поиска минимума функций многих переменных, где не вычисляются производные или градиент функции. Он сводится к построению симплекса в п-мерном пространстве, заданного п + 1 вершиной. В двухмерном пространстве симплекс является треугольником, а в трехмерном – пирамидой. На каждом шаге итераций выбирается новая точка решения внутри или вблизи симплекса. Она сравнивается с одной из вершин симплекса. Ближайшая к этой точке вершина симплекса заменяется этой точкой. Таким образом, симплекс перестраивается и позволяет найти новое, более точное положение точки решения. Алгоритм поиска повторяется, пока размеры симплекса по всем переменным не станет меньше заданной погрешности решения. Программу, реализующую симплекс-методы Нелдера-Мида, удобно использовать в следующей записи [x, min f] = f min search ( … ), где х – вектор координат локального минимума;

min f – значение целевой функции в точке минимума.

Саму целевую функцию удобно представить с помощью дескриптора @ в М-файле.

Пример 9. Найти и вывести на печать координаты и значение минимума функции двух переменных f(x, y) = (x2 + y2 – 3)2 + (x2 + y2 – 2x – 3)2 + 1, если начальная точка поиска имеет координаты М0 (1; 1). Анализ функции показывает, что min f = 1 x = 0, y = 3 = 1.73205.

Строим трехмерный график этой функции, чтобы убедиться в наличии минимума. Возьмем интервал х є [-1; 1]; y є [1; 3].

Протокол программы >> Z = (X.^2 + Y.^2 – 3).^2 + (X.^2 + Y.^2 – 2*X – 3).^2 + 1 ;

После построения трехмерного графика выполняем поиск минимума. В М-файле программируем целевую функцию Решаем поставленную задачу в окне команд >> [xmin, mi f] = fminsearch ( @ Fxy, [1; 1] ) Результаты поиска Как видно, результаты решения задачи точные.

Варианты заданий. Найти и вывести на печать координаты и минимальное значение функции двух переменных. Поиск начать с точки М0 (х0, у0). Данные взять из таблицы 2.9.

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

ГЛАВА 3. Применение MATLAB для решения задач строительной механики стержневых систем Изложенные выше основные данные о программировании и использовании численных методов в интегрированном пакете математического моделирования MATLAB показывают его большие возможности в решении задач прикладных наук. Одной из наиболее значимых прикладных наук в механико-инженерном образовании является строительная механика, где рассматриваются вопросы расчета машиностроительных, судостроительных, авиационных, строительных и т. п. конструкций. Объединение в одном пакете программирования и визуализации результатов расчетов позволяет существенно обогатить учебный процесс, повысить его методический уровень и наглядность теории, значительно уменьшить время на выполнение индивидуальных заданий и обеспечить глубокое овладение современными компьютерными технологиями. По этим причинам внедрение MATLAB в учебный процесс является актуальной задачей высшей школы. Важным условием высокого научного уровня учебного процесса является также широкое использование современных достижений соответствующих наук. В частности, в строительной механике накоплен большой арсенал методов расчета стержневых конструкций на статическую, динамическую нагрузки и устойчивость. Одним из наиболее эффективных методов строительной механики является метод граничных элементов (МГЭ), разработанный рядом зарубежных и отечественных ученых. В данной книге применен численно-аналитический вариант МГЭ, который с исчерпывающей полнотой изложен в работе [2]. Как ниже будет показано, реализация МГЭ в пакете MATLAB исключает какое-либо программирование численных процессов. Сами программы представляют собой лишь процесс введения исходных матриц и вызов численных процедур, т. е. реализуются в полном объеме идеи структурного программирования, что значительно уменьшает количество ошибок и повышает достоверность результатов. Дальнейшая графическая визуализация результатов расчетов логически завершает постановку задачи и ее решение. Таким образом, объединение современных достижений науки и компьютерных технологий будет способствовать подготовке высококвалифицированных специалистов и магистров. Кратко рассмотрим суть МГЭ в задачах деформирования стержневых систем.

3.1 Основные положения метода граничных элементов Значительное число задач механики упругого стержня сводится к решению линейного неоднородного дифференциального уравнения с постоянными коэффициентами удовлетворяющего начальным условиям Как известно, такая задача определения частного решения уравнения (3.1), удовлетворяющего условием (3.2), называется задачей Коши. Для изгиба, поперечных колебаний, кручения тонкостенных стержней, продольнопоперечного изгиба и т. п. видов деформирования, решение задачи Коши можно записать в матричной форме.

или компактно где Y(x) – вектор параметров напряженно-деформированного состояния стержня в текущей точке;

A(x) – квадратная матрица фундаментальных ортонормированных функций уравнения (3.1);

X(o) – вектор начальных параметров;

B – вектор (матрица-столбец) внешней нагрузки.

Если несколько стержней соединены в единую конструкцию, то для системы стержней можно составить матричное уравнение типа (3.4). Матрица А преобразуется к квазидиагональному виду, а векторы Y, X и B будут содержать параметры состояния всех стержней по следующей структуре Если координате х каждого стержня дать граничное значение l i, то для системы матриц (3.5) можно выполнить достаточное простое преобразование по схеме где конечные граничные параметры матрицы Y переносятся на место нулевых параметров вектора Х. При этом, эти векторы дополняются уравнениями равновесия и совместности перемещений узловых точек и граничными условиями. В конце схемы преобразований (3.6) получается система линейных алгебраических уравнений относительно начальных и конечных параметров всех стержней конструкции. После вычисления начальных параметров стержней их напряженно-деформированное состояние определяется по матричному уравнению (3.3). Таким образом, решение прямых задач строительной механики стержневых систем в МГЭ сводится к решению одной системы линейных алгебраических уравнений и вычислению напряженнодеформированного состояния в промежуточных точках по соотношениям метода начальных параметров. Такая схема решения обеспечивает получение весьма точных и достоверных результатов, которые можно представить средствами MATLAB в виде обычных эпюр, форм свободных колебаний, потерь устойчивости и т. п.

Главной операцией в схеме (3.6) является перенос параметров из Y в X.

Процесс переноса конечных параметров вектора Y в вектор X основан на следующих положениях. Векторы Y, X любой стержневой (и не стержневой) конструкции при граничном значении координаты х = l i будут содержать группы параметров.

Первая группа – это нулевые граничные параметры, что определяется заданными условиями опирания (краевыми условиями).

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

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

Эти параметры условно могут быть названы независимыми. Перенос параметров из вектора Y в вектор X должен компенсироваться ненулевыми элементами матрицы А, иначе нарушается исходное уравнение схемы (3.6).

Очевидно, что независимые параметры вектора Y должны быть перенесены на место нулевых параметров вектора X, а зависимые параметры переносятся в соответствии с уравнениями их связи. Перед операцией переноса параметров необходимо освободить поля матрицы А от элементов, связанных с нулевыми параметрами вектора Х, т. е. обнулить столбцы матрицы А, номера которых равны номерам нулевых строк матрицы Х. Далее в матрицу А вводятся ненулевые компенсирующие элементы и преобразования по схеме (3.6) завершены, поскольку в матрице В изменяются только знаки элементов. Правило для определения величины и положения компенсирующих элементов при переносе параметров включает 3 основных случая.

1-й случай. При переносе независимого параметра вектора Y в вектор X компенсирующий элемент матрицы А равен коэффициенту при переносимом параметре со своим знаком по схеме т. е. компенсирующий элемент матрицы А равен ( а) и должен появиться на месте (k, i), где k – номер строки матрицы Y, где находился параметр, i – номер строки матрицы Х, куда переносится параметр. Другими словами, первый индекс положения компенсирующего элемента указывает на старый адрес переносимого параметра, а второй индекс – новый адрес в матрице Х.

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

Более подробно алгоритм формирования системы линейных алгебраических уравнений типа (3.6) для стержневых систем рассмотрен ниже, при решении конкретных задач.

3.2. Расчет неразрезной балки 3.2.1. Определение граничных параметров при статической Исходные данные расчета балки:

l1 = 4.0; l2 = 6.0; l3 = 3.0; l4 = 1.0; m = 20.0; q = 10.0; f1 = - 40.0; f2 = 60.0;

am = 2.0; af1 = 2.0; af2 = 1.0.

Система линейных уравнений для граничных параметров балки Программа расчета граничных параметров балки % Описание матрицы а(16,16); b(16,1) и Х(16,1);

а = zeros(16,16); ); b = zeros (16,1); Х = zeros (16,1);

% Исходные данные расчета балки l1 = 4.0; l2 = 6.0; l3 = 3.0; l4 = 1.0; m = 20.0; q = 10.0;

f1 = - 40.0; f2 = 60.0; am = 2.0; af1 = 2.0; af2 = 1.0;

% Ввод матриц а(16,16) и b(16,1) с помощью операторов присваивания;

a(1,2) = l1; a(1,4) = -l1^3/6; b(1,1) = -m*(l1-am)^2/2;

a(2,2) = 1.0; a(2,4) = -l1^2/2; a(2,6) = -1.0; b(2,1) = -m*(l1-am);

a(3,4) = l1; a(3,7) = -1.0; b(3,1) = m; a(4,1) = -1.0; a(4,4) = 1.0;

a(5,6) = l2; a(5,7) = -l2^2/2; a(5,8) = -l2^3/6; b(5,1) = -q*l2 ^4/24;

a(6,6) = 1.0; a(6,7) = -l2; a(6,8) = -l2^2/2; a(6,10) = -1.0;

b(6,1) = -q*l2 ^3/6; a(7,7) = 1.0; a(7,8) = l2; a(7,11) = -1.0;

b(7,1) = q*l2^2/2; a(8,3) = -1.0; a(8,8) = 1.0; b(8,1) = q*l2;

a(9,10) = l3; a(9,11) = -l3^2/2; a(9,12) = -l3^3/6;

b(9,1) = -f1*(l3-af1)^3/6; a(10,10) = 1.0; a(10,11) = -l3;

a(10,12) = -l1^2/2; a(10,14) = -1.0; b(10,1) = -f1*(l3-af1)^2/2;

a(11,11) = 1.0; a(11,12) = l3; a(11,15) = -1.0; b(11,1) = f1*(l3-af1);

a(12,5) = -1.0; a(12,12) = 1.0; ; b(12,1) = f1; a(13,9) = -1.0;

a(13,14) = l4; a(13,15) = -l4^2/2; a(13,16) = -l4^3/6;

b(13,1) = -f2*(l4-af2)^3/6; a(14,13) = -1.0; a(14,14) = 1.0;

a(14,15) = -l4; a(14,16) = -l4^2/2; b(14,1) = -f2*(l4-af2)^2/2;

a(15,15) = -1.0; a(15,16) = l4; b(15,1) = f2*(l4-af2); a(16,16) = 1.0;

b(16,1) = f2;

% Решение системы уравнений ах = b и вывод результатов Х = a\ b Результаты вычисления граничных параметров представлены в таблице 3.1.

EI (0o)1 = X (2,1) = 12.7984 кНм 2 EI (2o)3 = X (10,1) = 53.7860 кНм M (1o)2 = X (7,1) = 24.1975 кНм 15 M (3o)4 = X (15,1) = 60.0000 кНм Как видно из этого примера, не требуется переставлять строки матриц А и В для исключения нулевых ведущих элементов, что является дополнительным преимуществом MATLAB по сравнению с непосредственным программированием метода исключения Гаусса в таких средах, как Delphi, Visual Fortran, Visual C++ и т.п.

3.2.2. Построение эпюр прогибов, углов поворота, поперечных сил Эпюры параметров изгиба можно строить для каждого стержня в отдельности, используя соотношения метода начальных параметров. При этом удобно выполнять построение эпюр в одном окне, для чего можно привлекать процедуру subplot. Для соблюдения требуемого масштаба можно использовать процедуру axis. Масштабная сетка эпюр появится, если задействовать процедуру grid on. Выполнение этих требований приводит к следующим операторам построения эпюр subplot (2, 2, 1), plot (x, EI); axis ( [0 2 -150 150] ); grid on subplot (2, 2, 2), plot (x, EIfi); axis ( [0 2 -150 150] ); grid on subplot (2, 2, 3), plot (x, Q); axis ( [0 2 -100 100] ); grid on subplot (2, 2, 4), plot (x, M); axis ( [0 2 -100 100] ); grid on.

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

Положительные направления внешней нагрузки и их координаты показаны на рис. 3. Для балки по рис. 3.1 протоколы построения эпюр примут вид:

Стержень 0 – x = 0 : 0.001 : 2.0;

EI = - (X(2,1).* – X(4,1).*x.^3/6);

EIfi = - (X(2,1) – X(4,1).*x.^2/2);

Q = X(4,1); M = X(4,1).*x;

subplot (2, 2, 1), plot (x, EI); axis ( [0 2 -150 150] ); grid on subplot (2, 2, 2), plot (x, EIfi); axis ( [0 2 -150 150] ); grid on subplot (2, 2, 3), plot (x, Q); axis ( [0 2 -100 100] ); grid on subplot (2, 2, 4), plot (x, M); axis ( [0 2 -100 100] ); grid on.

После построения графиков (нажатия кнопки ВЫПОЛНИТЬ) необходимо продолжить формирование команд для завершения эпюр на интервале (2; 4) в М-файле x = 2.0 : 0.001 : 4.0; m = 20.0;

EI = - (X(2,1).* – X(4,1).*x.^3/6 + m.*(x – 2).^2/2);

EIfi = - (X(2,1) – X(4,1).*x.^2/2 + m.*(x – 2));

Q = X(4,1); M = X(4,1).*x – m;

subplot (2, 2, 1), plot (x, EI); axis ( [2 4 -150 150] ); grid on subplot (2, 2, 2), plot (x, EIfi); axis ( [2 4 -150 150] ); grid on subplot (2, 2, 3), plot (x, Q); axis ( [2 4 -100 100] ); grid on subplot (2, 2, 4), plot (x, M); axis ( [2 4 -100 100] ); grid on.

Стержень 1 – x = 0 : 0.001 : 6.0; q = 10.0;

EI = - (X(6,1).*x – X(7,1).*x.^2/2 – X(8,1).*x.^3/6 + q.*x.^4/24);

EIfi = - (X(6,1) – X(7,1).*x – X(8.1).*x.^2/2 + q.*x.^3/6);

Q = X(8,1) – q.*x; M = X(7,1) + X(8,1).*x – q.*x.^2/2;

subplot (2, 2, 1), plot (x, EI); axis ( [0 6 -150 150] ); grid on subplot (2, 2, 2), plot (x, EIfi); axis ( [0 6 -150 150] ); grid on subplot (2, 2, 3), plot (x, Q); axis ( [0 6 -100 100] ); grid on subplot (2, 2, 4), plot (x, M); axis ( [0 6 -100 100] ); grid on.

Стержень 2 – x = 0 : 0.001 : 2.0;

EI = - (X(10,1).*x – X(11,1).*x.^2/2 – X(12,1).*x.^3/6);

EIfi = - (X(10,1) – X(11,1).*x – X(12,1).*x.^2/2);

Q = X(12,1); M = X(11,1) + X(12,1).*x;

subplot (2, 2, 1), plot (x, EI); axis ( [0 2 -150 150] ); grid on subplot (2, 2, 2), plot (x, EIfi); axis ( [0 2 -150 150] ); grid on subplot (2, 2, 3), plot (x, Q); axis ( [0 2 -100 100] ); grid on subplot (2, 2, 4), plot (x, M); axis ( [0 2 -100 100] ); grid on.

Далее вторая часть эпюр строится по операторам x =2. 0 : 0.001 : 3.0;

EI = - (X(10,1).*х – X(11,1).*x.^2/2 – X(12,1).*x.^3/6 – f1.*(x –2).^3/6);

EIfi = - (X(10,1) – X(11,1).*x – X(12,1).*x.^2/2 – f1.*(x –2).^2/2);

Q = X(12,1) + f1; M = X(11,1) + X(12,1).*x + f1.*(x –2);

subplot (2, 2, 1), plot (x, EI); axis ( [2 3 -150 150] ); grid on subplot (2, 2, 2), plot (x, EIfi); axis ( [2 3 -150 150] ); grid on subplot (2, 2, 3), plot (x, Q); axis ( [2 3 -100 100] ); grid on subplot (2, 2, 4), plot (x, M); axis ( [2 3 -100 100] ); grid on Стержень 3 – x = 0 : 0.001 : 1;

EI = - (X(14,1).*x – X(15,1).*x.^2/2 – X(16,1).*x.^3/6);

EIfi = - (X(14,1) – X(15,1).*x– X(16,1).*x.^2/2);

Q = X(16,1); M = X(15,1) + X(16,1).*x;

subplot (2, 2, 1), plot (x, EI); axis ( [0 1 -150 150] ); grid on subplot (2, 2, 2), plot (x, EIfi); axis ( [0 1 -150 150] ); grid on subplot (2, 2, 3), plot (x, Q); axis ( [0 1 -100 100] ); grid on subplot (2, 2, 4), plot (x, M); axis ( [0 1 -100 100] ); grid on.

Следует помнить, что в окно команд необходимо вводить значения неизвестных начальных параметров балки, т. е. вектор Х. Это легко выполняется с помощью отдельного М-файла, на котором записана процедура решения матричного уравнения неразрезной балки. После вектора Х нужно ввести также значения внешней нагрузки и их координаты. Эпюры EI(х), EI(х), Q(х) и М(х) представлены на рис. 3.1.

3.2.3. Определение частот собственных колебаний Матрица фундаментальных функций поперечных колебаний Частоты определяются как корни трансцендентного уравнения где матрица А() берется из задачи статики с заменой матриц фундаментальных функций изгиба на матрицы функций поперечных колебаний Наиболее просто корни могут быть найдены методом перебора в сочетании с процедурой вычисления определителя матрица А(). Организуется цикл вычисления определителя А(), значения которого вместе с частотой выводятся на печать в формате long e. При просмотре таблицы их значений определяются точки, где определитель изменяет знак. Эти точки и будут являться частотами собственных колебаний. Иногда бывают случаи, когда график функции d = А() касается оси. Для контроля этих точек полезно построить график у = d(), просмотр которого позволяет определить интервалы, где находятся корни и существуют ли точки касания графика оси. Далее корни могут быть уточнены повторными прогонами программы (без чального значения частоты 0 = 0.01l 2, с шагом 0 = 0.01l 2. Число вычислений п = 300 – 500 позволяет надежно определить первую частоту, старшие частоты определяются при расширении интервалов вычисления определителя А(). При решении данной задачи принимается, что l=EI=m=1.

Обозначения переменных, принятых в программе а матрица А(); п – порядок матрицы А(); п1 – число циклов вычисления определителя А(); аm – частота собственных колебаний; dam – шаг изменения частоты ; m – счетчик цикла вычисления определителя; l1, l2, l3, l4 – длины стержней неразрезной балки; lа – коэффициент фундаментальных функций; Х – вектор значений частоты ; Y – вектор значений определителя А().

Текст программы l1 = 4.0; l2 = 6.0; l3 = 3.0; l4 = 1.0; n = 16; n1 = 300; am = 0,01;

dam = 0.01; X = zeros (n1,1); Y = zeros (n1,1);

for m = 1 : n la = sqrt (am); a = zeros (n,n);

a(1,2) = (sinh(la*l1) + sin (la*l1))/(2*la);

a(1,4) = (si h(la*l1) - sin (la*l1))/(2*la^3);

a(2,2)=(cosh(la*l1)+cos(la*l1))/2; a(2,4) = (cos h(la*l1) - cos (la*l1))/(2*la^2);

a(2,6) = 1; a(3,2) = la^ 4*a(1,4); a(3,4) = a(1,2); a(3,7) = 1; a(4,1) = 1;

a(4,2) = la^ 4*a(2,4); a(4,4) = a(2,2);

a(5,6) = (sinh(la*l2) + sin (la*l2))/(2*la);

a(5,7)= (cosh(la*l2) cos(la*l2))/(2*la^2);

a(5,8) = (sinh(la*l2) sin (la*l2))/(2*la^ 3);

a(6,6)=(cosh(la*l2)+cos(la*l2))/2; a(6,7) = a(5,6); a(6,8) = a(5,7);

a(6,10) = 1; a(7,6) = la^ 4*a(5,8); a(7,7) = a(6,6); a(7,8) = a(5,6); a(7,11) = 1;

a(8,3) = 1; a(8,6) = la^ 4*a(5,7); a(8,7) = a(7,6); a(8,8) = a(6,6);

a(9,10) = (sinh(la*l3) + sin (la*l3))/(2*la);

a(9,11)= (cosh(la*l3) cos(la*l3))/(2*la^2);

a(9,12) = (sinh(la*l3) sin (la*l3))/(2*la^ 3);

a(10,10)=(cosh(la*l3)+cos(la*l3))/2; a(10,11) = a(9,10); a(10,12) = a(9,11);

a(10,14) = 1; a(11,10) = la^ 4*a(9,12); a(11,11) = a(10,10); a(11,12) = a(9,10);

a(11,15) = 1; a(12,5) = 1; a(12,10) = la^ 4*a(9,11); a(12,11) = a(11,10);

a(12,12) = a(10,10); a(13,9) = 1; a(13,14) = (sin h(la*l4) + sin (la*l4))/(2*la);

a(13,15)= (cosh(la*l4) cos(la*l4))/(2*la^2);

a(13,16) = (sinh(la*l4) sin (la*l4))/(2*la^ 3); a(14,13) = 1;

a(14,14)=(cosh(la*l4)+cos(la*l4))/2; a(15,14) = la^ 4*a(13,16);

a(15,15) = a(14,14); a(15,16)= a(13,14); a(16,14) = la^ 4*a(13,15);

a(16,15) = a(15,14); a(16,16) = a(14,16);

d = det(a); X(m,1) = am; Y(m,1) = d; am = am + dam; end;

% Построение графика у = d(am) plot (X,Y); grid оn % Вывод на печать таблицы значений am и d format long e [X Y] График определителя d = A() для балки принимает вид (рис. 3.5) Как следует из рис. 3.5, график определителя d = A() не имеет точек разрыва 2-го рода (это большое преимущество метода граничных элементов), не касается оси и имеет точки пересечения с горизонталью. Эти точки отмечены символами 1, 2 и т. д. Уточнение частот приводит к следующим значениям

EI EI EI

3.2.4. Построение форм собственных колебаний Для построения форм собственных колебаний необходимо определить граничные параметры стержней балки. Эта задача может быть решена следующим образом:

1. Используется статический или кинематический способ возбуждения собственных колебаний. Для этого приравнивается единице статический или кинематический параметр матрицы Y( l ), равный нулю. Далее этот параметр переносится в матрицу нагрузки В. Из системы уравнений (3.11) следует, что наиболее удобным параметром является Q3-4( l ) = 1, тогда в (16,1) = 1.

(3.16) Решением линейной системы уравнений (при собственной частоте I) определяются граничные параметры балки, которые нужно нормировать относительно какого-либо перемещения. Для данной балки в качестве нормирующей величины можно взять перемещение консольной части 2. По нормированным начальным параметрам строятся формы собственных колебаний балки средствами MATLAB.

Уравнение для граничных параметров балки отличается от уравнения (3.11) только фундаментальными функциями и вектором правой части (см. (3.16)).

Программа вычисления граничных параметров набирается в отдельном М-файле и принимает вид a = zeros (16,16); b = zeros (16,1); X = zeros (16,1); la = sqrt (0,4055125);

l1 = 4.0; l2 = 6.0; l3 = 3.0; l4 = 1.0;

a(1,2) = (sinh(la*l1) + sin (la*l1))/(2*la);

a(1,4) = (sinh(la*l1) sin (la*l1))/(2*la^3);

a(2,2)=(cosh(la*l1)+cos(la*l1))/2; a(2,4) = (cosh(la*l1) - cos (la*l1))/(2*la^2);

a(2,6) = 1; a(3,2) = la^ 4*a(1,4); a(3,4) = a(1,2); a(3,7) = 1; a(4,1) = 1;

a(4,2) = la^ 4*a(2,4); a(4,4) = a(2,2); a(5,6) = (sinh(la*l2) + sin (la*l2))/(2*la);

a(5,7)= (cosh(la*l2) cos(la*l2))/(2*la^2);

a(5,8) = (sinh(la*l2) sin (la*l2))/(2*la^ 3); a(6,6)=(cosh(la*l2)+cos(la*l2))/2;

a(6,7) = a(5,6); a(6,8) =a(5,7); a(6,10) = 1; a(7,6) =la^ 4*a(5,8); a(7,7) =a(6,6);

a(7,8) = a(5,6); a(7,11) = 1; a(8,3) = 1; a(8,6) = la^ 4*a(5,7); a(8,7) = a(7,6);

a(8,8) = a(6,6); a(9,10) = (sinh(la*l3) + sin (la*l3))/(2*la);

a(9,11)= (cosh(la*l3) cos(la*l3))/(2*la^2);

a(9,12) =(sinh(la*l3)sin(la*l3))/(2*la^3); a(10,10)=(cosh(la*l3)+cos(la*l3))/2;

a(10,11) = a(9,10); a(10,12) = a(9,11); a(10,14) = 1; a(11,10) = la^ 4*a(9,12);

a(11,11) =a(10,10); a(11,12)= a(9,10); a(11,15)= 1; a(12,5)= 1; a(12,11) = a(11,10);

a(12,10)=la^ 4*a(9,11); a(12,12)=a(10,10); a(13,14)=(sinh(la*l4) + sin (la*l4))/(2*la);

a(13,9) = 1; a(13,15)= (cos h(la*l4) cos(la*l4))/(2*la^2);

a(13,16) = (sinh(la*l4) sin (la*l4))/(2*la^ 3); a(14,13) = 1;

a(14,14)=(cosh(la*l4)+cos(la*l4))/2; a(14,15) = - a(13,14); a(14,16) = a(13,15);

a(15,14) = la^ 4*a(13,16); a(15,15) = a(14,14); a(15,16) = a(13,14); a(16,14) = la^ 4*a(13,15);

a(16,15) = a(15,14); a(16,16) = a(14,14); b(16,1) = 1; X = a\ b;

X = X \ X(9,1) Результаты вычислений относительных значений граничных параметров балки при собственных частотах сведены в таблицу 3.2.

№ Граничные Относительные значения граничных параметров при частотах п/п параметры балки Выражение для прогиба (формы колебания) стержня при собственных колебаниях в соответствии с методом начальных параметров имеет вид Протокол построения форм колебаний балки принимает вид х1 = 0 : 0.001 : 4.0; х2 = 0 : 0.001 : 6.0; х3 = 0 : 0.001 : 3.0;

х4 = 0 : 0.001 : 1.0; la = sqrt (0.4055125);

EIv1 = (X(2,1)*(sinh(la*x1)+sin(la*x1))/(2*la) … X(4,1)*(sinh(la*x1)sin(la*x1))/(2*la^3));

EIv2 = (X(6,1)*(sinh(la*x2)+sin(la*x2))/(2*la) … X(7,1)*(cosh(la*x2)cos(la*x2))/(2*la^2) … X(8,1)*(sinh(la*x2)sin(la*x2))/(2*la^3));

EIv3 = (X(10,1)*(sinh(la*x3)+sin(la*x3))/(2*la) … X(11,1)*(cosh(la*x3)cos(la*x3))/(2*la^2) … X(12,1)*(sinh(la*x3)sin(la*x3))/(2*la^3));

EIv4 = (X(14,1)*(sinh(la*x4)+sin(la*x4))/(2*la) … X(15,1)*(cosh(la*x4)cos(la*x4))/(2*la^2) … X(16,1)*(sinh(la*x4)sin(la*x4))/(2*la^3));

subplot (2,2,1), plot (x1, EIv1); axis ([0 4 0.5 0.5]); grid on subplot (2,2,2), plot (x2, EIv2); axis ([0 6 0.6 0.6]); grid on subplot (2,2,3), plot (x3, EIv3); axis ([0 3 0.5 0.5]); grid on Перед выполнением данного протокола из отдельного М-файла в окно команд помещается вектор Х вектор значений граничных параметров балки при i. Выполняется это прогонкой программы, решающей систему уравнений (3.16). Далее отрабатывается протокол построения форм колебаний при собственной частоте i. Первые 5 форм собственных колебаний балки при условии, что EI(3l) 4 = 1, представлены на рис 3.6. Можно рекомендовать также объединить протоколы вычисления вектора Х и построения форм колебаний, но только после четкого усвоения навыка построения графиков и выбора их масштабов.

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

3.2.5. Расчет на вынужденные колебания Расчет балки на вынужденные колебания сводится к решению матричного уравнения АХ = В, где матрица А повторяет матрицу частотного уравнения (см. п.3.2.3), а матрица В представлена в этом разделе. Программа также записывается в отдельном М-файле, куда переносится матрица А задачи определения частот собственных колебаний и дополняется элементами вектора В. Сама программа имеет вид.

a = zeros (16,16); b = zeros (16,1); X = zeros (16,1); la = sqrt (0.1*0,315365);

l1 = 4.0; l2 = 6.0; l3 = 3.0; l4 = 1.0; m = 20.0; q = 10.0; f1 = - 40.0;

f2 = 60.0; am = 2.0; af1 = 2.0; af2 = 1.0;

a(1,2) = (sinh(la*l1) + sin (la*l1))/(2*la);

a(1,4) = (sinh(la*l1) sin (la*l1))/(2*la^3);

a(2,2) = (cosh(la*l1)+cos(la*l1))/2;

a(2,4) = (cosh(la*l1) cos (la*l1))/(2*la^2); a(2,6) = 1;

a(3,2) = la^ 4*a(1,4); a(3,4) = a(1,2); a(3,7) = 1; a(4,1) = 1;

a(4,2) = la^ 4*a(2,4); a(4,4) = a(2,2);

b(1,1) = m (cosh(la*(l1am)) cos (la*(l1am))) / (2*la^2);

b(2,1) = m (sinh(la*(l1am)) + sin (la*(l1am))) / (2*la);

b(3,1) = m (cosh(la*(l1am)) + cos (la*(l1am))) / 2;

b(4,1) = m la*(sinh(la*(l1am)) sin (la*(l1am))) / 2;

a(5,6) = (sinh(la*l2) + sin (la*l2))/(2*la);

a(5,7)= (cosh(la*l2) cos(la*l2))/(2*la^2);

a(5,8) = (sinh(la*l2) sin (la*l2))/(2*la^ 3);

a(6,6)=(cosh(la*l2)+cos(la*l2))/2; a(6,7) = a(5,6); a(6,8) =a(5,7);

a(6,10) = 1; a(7,6) =la^ 4*a(5,8); a(7,7) =a(6,6); a(7,8) = a(5,6);

a(7,11) = 1; a(8,3) = 1; a(8,6) = la^ 4*a(5,7); a(8,7) = a(7,6);

a(8,8) = a(6,6);

b(5,1) = q (cosh(la*l2) + cos (la*l2) 2) / (2*la^4);

b(6,1) = q (sinh(la*l2) sin (la*l2)) / (2*la^3);

b(7,1) = q (cosh(la*l2) cos (la*l2)) / (2*la^2);

b(8,1) = q (sinh(la*l2) + sin (la*l2)) / (2*la);

a(9,10) = (sinh(la*l3) + sin (la*l3))/(2*la);

a(9,11)= (cosh(la*l3) cos(la*l3))/(2*la^2);

a(9,12) =(sinh(la*l3)sin(la*l3))/(2*la^3);

a(10,10)=(cosh(la*l3)+cos(la*l3))/2; a(10,11) = a(9,10); a(10,12) = a(9,11);

a(10,14) = 1; a(11,10) = la^ 4*a(9,12); a(11,11) =a(10,10); a(11,12)= a(9,10);

a(11,15)= 1; a(12,5)= 1; a(12,11) = a(11,10); a(12,10)=la^ 4*a(9,11);

a(12,12)=a(10,10);

b(9,1) = f1 (sinh(la*(l3a f1)) sin (la*(l3a f1))) / (2*la^3);

b(10,1) = f1 (cosh(la*(l3a f1)) cos (la*(l3a f1))) / (2*la^2);

b(11,1) = f1 (sinh(la*(l3a f1)) + sin (la*(l3a f1))) / (2*la);

b(12,1) = f1 (cosh(la*(l3a f1)) + cos (la*(l3a f1))) / 2; a(13,9) = 1;

a(13,14)=(sinh(la*l4) + sin (la*l4))/(2*la);

a(13,15)= (cosh(la*l4) cos(la*l4))/(2*la^2);

a(13,16) = (sinh(la*l4) sin (la*l4))/(2*la^ 3); a(14,13) = 1;

a(14,14)=(cosh(la*l4)+cos(la*l4))/2; a(14,15) = - a(13,14); a(14,16) = a(13,15);

a(15,14) = la^ 4*a(13,16); a(15,15) = a(14,14); a(15,16) = a(13,14);

a(16,14) = la^ 4*a(13,15); a(16,15) = a(15,14); a(16,16) = a(14,14);

b(13,1) = f2 (sinh(la*(l4a f2)) sin (la*(l4a f2))) / (2*la^3);

b(14,1) = f2 (cosh(la* l4a f2)) cos (la*(l4a f2))) / (2*la^2);

b(15,1) = f2 (sinh(la*(l4a f2)) + cos (la*(l4a f2))) / (2*la);

b(16,1) = f2 (cosh(la*(l4a f2)) + cos (la*(l4a f2))) / 2;

X = a\ b Для ухода от резонанса частоту вынужденных колебаний принимаем Тогда Результаты расчета динамических граничных параметров балки 3.2.6. Построение эпюр напряженно-деформированного состояния балки при вынужденных колебаниях Эпюры напряженно-деформированного состояния можно построить для каждого стержня в отдельности, используя соотношения метода начальных параметров.

где выражения от внешней нагрузки имеют вид M, F, q – сосредоточенные момент, сила и распределенная нагрузка.

аМ, аF, аH, аK – координаты внешней нагрузки.

Протоколы построения эпюр EI(x), EI(x), M(x) и Q(x) примут вид х = 0 : 0.0001 : 2.0; la = sqrt (0.1*0.315365);

EIv = (X(2,1)*(sinh(la*x)+sin(la*x))/(2*la) X(4,1)* … (sinh(la*x)sin(la*x))/(2*la^3));

EIfi = (X(2,1)*(cosh(la*x)+cos(la*x)) / 2 X(4,1)* … (cosh(la*x)cos(la*x))/(2*la^2));

Q = X(2,1)*la^4*(cosh(la*x)cos(la*x)) / (2*la^2) + … X(4,1)* (cosh(la*x)+cos(la*x)) / 2;

M = X(2,1)*la^4*(sinh(la*x)sin(la*x))/(2*la^3) + … X(4,1)*(sinh(la*x)+sin(la*x))/(2*la);

subplot (2,2,1), plot (x, EIv); axis ([0 2 30 30]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 2 30 30]); grid on subplot (2,2,3), plot (x, Q); axis ([0 2 10 10]); grid on Перед выполнением этого протокола необходимо поместить в окно команд вектор граничных параметров Х (см. п. 3.2.5), значения нагрузки и их координаты.

х = 2.0 : 0.0001 : 4.0; m = 20.0; am = 2.0; la = sqrt (0.1*0.315365);

EIv = (X(2,1)*(sinh(la*x)+sin(la*x))/(2*la) X(4,1)* … (sinh(la*x)sin(la*x))/(2*la^3)+m*(cosh(la*(x-am)) … cos (la*(x-am))) / (2*la^2));

EIfi = (X(2,1)*(cosh(la*x)+cos(la*x)) / 2 X(4,1)*(cosh(la*x) … cos(la*x))/(2*la^2))+ m*(sinh(la*(x-am))+sin(la*(x-am))) / (2*la));

Q = X(2,1)*la^4*(cosh(la*x)cos(la*x)) / (2*la^2) + … X(4,1)* (cosh(la*x)+cos(la*x)) / 2 m*(sinh(la*(x-am)) … sin(la*(x-am)))* la / 2;

M = X(2,1)*la^4*(sinh(la*x)sin(la*x))/(2*la^3) + X(4,1)* … (sinh(la*x)+sin(la*x))/(2*la)) m*(cosh(la*(x-am)+ … cos(la*(x-am))) / 2;

subplot (2,2,1), plot (x, EIv); axis ([2 4 30 30]); grid on subplot (2,2,2), plot (x, EIfi); axis ([2 4 40 40]); grid on subplot (2,2,3), plot (x, Q); axis ([2 4 2 2]); grid on subplot (2,2,4), plot (x, M); axis ([2 4 30 30]); grid on стержень 1- х = 0 : 0.0001 : 6.0; q = 10.0; la = sqrt (0.1*0.315365);

EIv = (X(6,1)*(sinh(la*x)+sin(la*x))/(2*la) X(7,1)* … (cosh(la*x)cos(la*x))/(2*la^2)) X(8,1)*(sinh(la*x)sin(la*x) … ) / (2*la^3)) + q* (cos(la*x) + cos(la*x) 2) / (2*la^4));

EIfi = (X(6,1)*(cosh(la*x)+cos(la*x)) / 2 X(7,1)* (sinh(la*x) + … sin(la*x)) / (2* la) X(8,1)* (cosh(la*x) cos(la*x)) / (2*la^2) + … q*(sinh(la*x)sin(la*x))/(2*la^3));

Q = X(6,1)*la^4*(cosh(la*x)cos(la*x)) / (2*la^2) + X(7,1)* … la^4 *(sinh(la*x)sin(la*x))/(2*la^3) + X(8,1)* (cosh(la*x) + cos(la*x) … ) / 2 q*(sinh(la*x)+sin(la*x))/(2*la);

M = X(6,1)*la^4*(sinh(la*x)sin(la*x))/(2*la^3) + X(7,1)* … (cosh(la*x)+cos(la*x)) / 2 + X(8,1)*(sinh(la*x)+sin(la*x))/(2*la) … q*(cosh(la*x)cos(la*x)) / (2*la^2);

subplot (2,2,1), plot (x, EIv); axis ([0 6 110 110]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 6 60 60]); grid on subplot (2,2,3), plot (x, Q); axis ([0 6 40 40]); grid on subplot (2,2,4), plot (x, M); axis ([0 6 40 40]); grid on стержень 2- х = 0 : 0.0001 : 2.0; la = sqrt (0.1*0.315365);

EIv = (X(10,1)*(sinh(la*x)+sin(la*x))/(2*la) X(11,1)* ( … cosh(la*x)cos(la*x))/(2*la^2)X(12,1)*(sinh(la*x)sin(la*x)) … /(2*la^3));

EIfi = (X(10,1)*(cosh(la*x)+cos(la*x)) / 2 X(11,1)*(sinh(la*x) + … sin(la*x))/(2*la) X(12,1)*(cosh(la*x)cos(la*x))/(2*la^2));

Q = X(10,1)*(cosh(la*x)cos(la*x))*la^4 / (2*la^2) + X(11,1)* … la^4*(sinh(la*x)sin(la*x))/(2*la^3) + X(12,1)*(cosh(la*x) + … cos(la*x)) / 2;

M = X(10,1)*(sinh(la*x)sin(la*x))*la^4/(2*la^3) + X(11,1)* … (cosh(la*x)+cos(la*x)) / 2 + X(12,1)*(sinh(la*x) +(sinh(la*x))/(2*la);

subplot (2,2,1), plot (x, EIv); axis ([0 2 60 60]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 2 90 90]); grid on subplot (2,2,3), plot (x, Q); axis ([0 2 40 40]); grid on subplot (2,2,4), plot (x, M); axis ([0 2 70 70]); grid on Протоколы эпюр второй части стержня имеют вид:

х = 2.0 : 0.0001 : 3.0; la = sqrt (0.1*0.315365); f1 = 40.0; af1 = 2.0;

EIv = (X(10,1)*(sinh(la*x)+sin(la*x))/(2*la) X(11,1)* ( … cosh(la*x)cos(la*x)) / (2*la^3) X(12,1)*(sinh(la*x)sin(la*x)) … /(2*la^3)+ f1*(sinh(la*(x-a f1))sin(la*(x-a f1))) / (2*la^3));

EIfi = (X(10,1)*(cosh(la*x)+cos(la*x)) / 2 X(11,1)*(sinh(la*x)+ … sin(la*x))/(2*la)) X(12,1)*(cosh(la*x) cosh(la*x)) / (2*la^2) + … f1*(cosh(la*(x-a f1))cos(la*(x-a f1))) / (2*la^2));

Q = X(10,1*(cosh(la*x)cos(la*x))*la^4 / (2*la^2) + X(11,1)* … la^4*(sinh(la*x)sin(la*x))/(2*la^3) + X(12,1)*(cosh(la*x) + … cos(la*x)) / 2 f1*(cosh(la*(x-a f1))+cos(la*(x-a f1))) / 2;

M = X(10,1)*(sinh(la*x)sin(la*x))*la^4/(2*la^3) + X(11,1)* … (cos(la*x)+ cos(la*x)) / 2 + X(12,1)*(sinh(la*x)+sin(la*x))/(2*la) … f1*(sinh(la*(x-a f1))+sin(la*(x-a f1))) / (2*la);

subplot (2,2,1), plot (x, EIv); axis ([2 3 60 60]); grid on subplot (2,2,2), plot (x, EIfi); axis ([2 3 90 90]); grid on subplot (2,2,3), plot (x, Q); axis ([2 3 10 10]); grid on subplot (2,2,4), plot (x, M); axis ([2 3 70 70]); grid on х = 0 : 0.0001 : 1.0; la = sqrt (0.1*0.315365);

EIv = (X(14,1)*(sinh(la*x)+sin(la*x))/(2*la) X(15,1)*(cosh(la*x) … cos(la*x))/(2*la^2) X(16,1)*(sinh(la*x)sin(la*x))/(2*la^3));

EIfi = (X(14,1)*(cosh(la*x)+cos(la*x)) / 2 X(15,1)*(sinh(la*x)+ … sinh(la*x))/(2*la) X(16,1)*(cosh(la*x)cos(la*x))/(2*la^2));

Q = X(14,1)*la^4*(cosh(la*x)cos(la*x)) / (2*la^2) + X(16,1)* la^4* … (sinh(la*x)sin(la*x))/(2*la^3)+X(16,1)*(cosh(la*x)+cos(la*x)) / 2;

M=X(14,1)*la^4*(sinh(la*x)sin(la*x))/(2*la^3)+X(15,1)*(cosh(la*x)+… cos(la*x))/2+X(16,1)*(sinh(la*x)+sin(la*x))/(2*la);

subplot (2,2,1), plot (x, EIv); axis ([0 1 110 110]); grid on subplot (2,2,2), plot (x, EIfi); axis ([0 1 120 120]); grid on subplot (2,2,3), plot (x, Q); axis ([0 1 70 70]); grid on subplot (2,2,4), plot (x, M); axis ([0 1 70 70]); grid on Анализ показывает, что несмотря на большой объем протоколов построения эпюр параметров балки, сам процесс довольно простой. В отдельном М-файле необходимо только дополнять недостающие элементы параметров и менять индексы вектора Х. Эпюры EI(x), EI(x), Q(x) и М(x) представлены на рис. 3.7.

3.2.7. Определение критических сил потери устойчивости Исходные данные: l 01 = 4.0 м. ; l122 = 6.0 м. ; l 23 = 3.0 м.; l 344 = 1.0 м. ; EI = 1.

Матрица фундаментальных функций продольно-поперечного изгиба (задача устойчивости плоских стержневых систем) Критические силы определяются как корни трансцендентного уравнения А(F) = 0, где матрицу А можно взять из задач статики или динамики данной балки с заменой фундаментальных функций. Для неразрезной балки можно не учитывать нормальные силы N, потому матрица А(F) примет размер 16 16 строк и столбцов.

Корни (критические силы) можно определить методом последовательного перебора в сочетании с прямым ходом метода исключения Гаусса. Организуется цикл вычисления определителя А(F), в конце которого величины F и d выводятся в окно команд. При просмотре таблицы значений F и d определяются точки, где изменяется знак определителя d = А(F). Эти точки и есть критические силы потери устойчивости. Рекомендуется начинать вычисления с начального значения сжимающей силы F0 = 0.001 с шагом F = 0.001. Число вычислений определителя d = А(F) п1 = 300 – 500 позволяет надежно и достаточно точно определить первую и старшие критические силы. При решении задачи принимаем, что EI = 1. Матрица устойчивости неразрезной балки принимает вид А(F)= Обозначения переменных, принятых в программе а – матрица А(F); п – порядок матрицы А(F); п1 – число циклов вычисления определителя d = А(F) ; d – величина определителя; f – сжимающая сила F; df – шаг изменения сжимающей силы F; m – счетчик цикла вычисления определителя; l1 l4 длины стержней балки; п2 – коэффициент п фундаментальных функций; X, Y – векторы значений сжимающих сил F и определителя d.

п = 16; п1 = 300; f = 0.01; df = 0.01; X = zeros (n1,1); Y = zeros (n1,1);

l1 = 4.0; l2 = 6.0; l3 = 3.0; l4 = 1.0;

for m = 1 : n1 п2 = sqrt (f); a = zeros (n,n);

a(1,2) = sin (п2*l1) / п2; a(1,4) = (п2*l1sin (п2*l1)) / (п2^3);

a(2,2)=cos (п2*l1); a(2,4) = (1cos (п2*l1)) / (п2^2); a(2,6) = 1;

a(3,2) = п2*sin(п2*l1); a(3,7) = 1; a(4,1) = 1; a(4,4) = 1;

a(5,6) = sin (п2*l2) / п2; a(5,7)= (1cos (п2*l2)) / (п2^2);

a(5,8) = (п2*l2sin (п2*l2)) / (п2^ 3); a(6,6)=cos (п2*l2); a(6,7) = a(5,6);

a(6,8) = a(5,7); a(6,10) = 1; a(7,6) = п2*sin(п2*l2); a(7,7) = a(6,6);

a(7,8) = a(5,6); a(7,11) = 1; a(8,3) = 1; a(8,6) = 1; a(9,10) = (sin (п2*l3) / п2;

a(9,11)= (1cos (п2*l3)) / (п2^2); a(9,12) = (п2*l3sin (п2*l3)) / (п2^ 3);

a(10,10)=cos (п2*l3); a(10,11) = a(9,10); a(10,12) = a(9,11); a(10,14) = 1;

a(11,10) = п2*sin(п2*l3); a(11,11) = a(10,10); a(11,12) = a(9,10); a(11,15) = 1;

a(12,5) = 1; a(12,12) = 1; a(13,9) = 1; a(13,14) = sin (п2*l4) / п2;

a(13,15)= (1cos (la*l4)) / (п2^2); a(13,16) = (п2*l4sin (п2*l4)) / (п2^ 3);

a(14,13) = 1; a(14,14)=cos (п2*l4); a(14,15) = a(13,14); a(14,16) = a(13,15);

a(15,14) = п2*sin(п2*l4); a(15,15) = a(14,14); a(15,16) = a(13,14);

a(16,16) = 1; d = det(a); X(m,1) = f; Y(m,1) = d; f = f + df; end;

% Построение графика определителя матрицы устойчивости plot (X,Y); grid on % Вывод таблицы значений f и d в формате long e format long e [X Y] График зависимости d = А(F) позволяет лишь грубо определить интервалы, где находятся корни (критические силы). Как и в динамике, этот график не имеет точек разрыва 2-го рода (в методах сил и перемещений аналогичный график имеет разрывы 2-го рода) и служит вспомогательным средством при поиске критических сил. Уточнить критические силы можно при повторных прогонах программы (без построения графиков) с новыми интервалами для F.

Результаты поиска критических сил потери устойчивости

EI EI EI

F1 = 0.443685 ; F2 = 0.652565 2 ; F3 = 1.051035 2 ;

F4 = 1.6187795 2 ; F5 = 2.430485 2 кН.

3.2.8. Построение форм потери устойчивости Аналогично задаче динамики (п.3.2.4) примем, что b(16,1) = Q3-4( l ) = 1.

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

Программа вычисления граничных параметров имеет вид а = zeros (16,16); b = zeros (16,1); X = zeros (16,1); п2 =sqrt(0.443685);

l1 = 4.0; l2 = 6.0; l3 = 3.0; l4 = 1.0;

a(1,2) =sin(п2*l1)/п2; a(1,4) = (п2*l1sin (п2*l1)) / (п2^3); a(2,2)=cos (п2*l1);

a(2,4) = (1cos (п2*l1)) / (п2^2); a(2,6) = 1; a(3,2) = п2*sin(п2*l1);

a(3,4) = a(1,2); a(3,7) = 1; a(4,1) = 1; a(4,4) = 1;

a(5,6) = sin (п2*l2) / п2; a(5,7)= (1cos (п2*l2)) / (п2^2);

a(5,8) = (п2*l2sin (п2*l2)) / (п2^ 3); a(6,6)=cos (п2*l2); a(6,7) = a(5,6);

a(6,8) = a(5,7); a(6,10) = 1; a(7,6) = п2*sin(п2*l2); a(7,7) = a(6,6);

a(7,8) = a(5,6); a(7,11) = 1; a(8,3) = 1; a(8,8) = 1;

a(9,10) = sin (п2*l3) / п2; a(9,11)= (1cos (п2*l3)) / (п2^2);

a(9,12) = (п2*l3sin (п2*l3)) / (п2^ 3); a(10,10)=cos (п2*l3);

a(10,11) = a(9,10); a(10,12) = a(9,11); a(10,14) = 1;

a(11,10) = п2*sin(п2*l3); a(11,11) = a(10,10); a(11,12) = a(9,10);

a(11,15) = 1; a(12,5) = 1; a(12,12) = 1; a(13,9) = 1;

a(13,14) = sin (п2*l4) / п2; a(13,15)= (1cos (la*l4)) / (п2^2);

a(13,16) = (п2*l4sin (п2*l4)) / (п2^ 3); a(14,13) = 1;

a(14,14)=cos (п2*l4); a(14,15) = a(13,14); a(14,16) = a(13,15);

a(15,14) = п2*sin(п2*l4); a(15,15) = a(14,14); a(15,16) = a(13,14);

a(16,16) = 1; b (16,1) = 1; X = a\ b; X = X / X(9,1) Результаты вычислений относительных значений граничных параметров балки при потере устойчивости сведены в таблицу 3.4.

(3.25) п/п Выражение для прогиба (формы потери устойчивости) согласно методу начальных параметров имеет вид Протокол построения форм потери устойчивости запишется в виде х4 = 0 : 0.001 : 1.0; п2 = sqrt (0.443685);

EIv1 = (X(2,1)*(sin(п2*x1) / п2 X(4,1)*(п2*x1sin(п2*x1)) … EIv2 = (X(6,1)*sin(п2*x2) / п2 X(7,1)*(1cos(п2*x2)) / п2^2) … X(8,1)*( п2*x2sin(п2*x2)) / п2^3);

EIv3 = (X(10,1)*sin(п2*x3) / п2 X(11,1)*(1cos(п2*x3)) / п2^2) … X(12,1)*( п2*x3sin(п2*x3)) / п2^3);

EIv4 = (X(14,1)*sin(п2*x4) / п2 X(15,1)*(1cos(п2*x4)) / п2^2 … X(16,1)*( п2*x4sin(п2*x4)) / п2^3);

subplot (2,2,1), plot (x1, EIv1); axis ([0 4 2 2]); grid on subplot (2,2,2), plot (x2, EIv2); axis ([0 6 2 2]); grid on subplot (2,2,3), plot (x3, EIv3); axis ([0 3 2 2]); grid on subplot (2,2,4), plot (x4, EIv4); axis ([0 1 2 2]); grid on Перед выполнением данного протокола из отдельного М-файла в окно команд помещается вектор Х вектор относительных значений граничных параметров балки. Выполняется это прогонкой программы, решающей систему уравнений (3.25). Далее нужно задействовать протокол построения форм потери устойчивости балки при разных критических силах Fi. Первые форм потери устойчивости балки при условии, что EI 34 (l ) = 1, представлены на рис 3.12.

3.2.9. Варианты заданий Продолжение рис. 3. варианта заданий 3.3. Расчет плоской рамы 3.3.1. Определение граничных параметров при статической Для построения эпюр напряженно-деформированного состояния рамы формируем систему линейных алгебраических уравнений краевой задачи по МГЭ. Для этого:

1. Разбиваем раму на 4 стержня. Нумеруем узлы и стрелками обозначаем начало и конец каждого стержня, т. е. формируем орграф расчета рамы (рис. 3.14, а).

2. Составляем уравнения равновесия и совместности перемещений узлов рамы. Уравнения равновесия узлов 1 и 2 составляем для недеформированного состояния, а уравнения совместности перемещений в соответствии с деформированным состоянием по рис. 3.14, d.

При расчете рамы полагаем, что стержни нерастяжимы и несжимаемы.

Составленные уравнения равновесия и совместности перемещений узлов рамы помещаем в матрицу конечных параметров Y (l ). Матрицы Х и Y (l ) при учете граничных условий примут вид Из анализа матрицы Х следует, что в матрице А нужно обнулить 1, 3, 5, 6, 7 и 11 столбцы. На место нулевых строк матрицы Х переносим независимые параметры матрицы Y. Зависимые параметры матрицы Y переносим в матрицу Х в соответствии с уравнениями их связи. В матрице А появятся компенсирующие элементы. Разрешающее уравнение задачи статики при l = 1 м примет вид Программа определения граничных параметров рамы запишется следующим образом (в отдельном файле) a = zeros (20,20); b = zeros (20,1); X = zeros (20,1);

a(1,2) = 1; a(1,4) = 1/6; a(2,2) = 1; a(2,4) = 1/6; a(2,17) = 1;

a(3,4) = 1; a(3,8) = 1; a(3,18) = 1; a(4,4) = 1; a(4,9) = 1; a(4,20) = 1;

a(5,10) = 1; a(5,19) = 1; b(1,1) = 5/12; b(2,1) = 5/3; b(3,1) = 5; b(4,1) = 10;

a(6,8) = 1/2; a(6,9) = 1/6; a(6,17) = 1; a(7,8) = 1; a(7,9) = 1/2; a(7,12) = 1;

a(7,17) = 1; a(8,8) = 1; a(8,9) = 1; a(8,13) = 1; a(9,9) = 1; ; a(9,15) = 1;

a(10,10) = 1; a(10,14) = 1; b(6,1) = 5/12; b(7,1) = 5/3; b(8,1) = 5; b(9,1) = 10;