Моделирование установившихся колебаний пластин, контактирующих с упругим основанием

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

«КУБАНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ»

(ФГБОУ ВПО «КубГУ»)

Кафедра математического моделирования

ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ РАБОТА Бакалавра

моделирование установившихся колебаний пластин, контактирующих с упругим основанием

Работу выполнила Н.А. Рябец

Факультет компьютерных технологий и прикладной математики

Направление 010500.62 - Прикладная математика и информатика

Научный руководитель,

профессор, д-р физ.-мат наук, доцент А.В. Павлова

Нормоконтролер,

доцент, канд. физ.-мат наук М.С. Капустин

Краснодар 2014

РЕФЕРАТ

Выпускная квалификационная работа 51 с., 25 рис., 14 источников, 4 приложения.

ТЕОРИЯ УПРУГОСТИ, ДЕФОРМАЦИИ И ПЕРЕМЕЩЕНИЯ, ПОКРЫТИЕ, ПОДЛОЖКА, УСТАНОВИВШИЕСЯ КОЛЕБАНИЯ, ПРЕОБРАЗОВАНИЕ ФУРЬЕ.

Объектом исследования являются контактные взаимодействия упругой полосы и пластины-покрытия.

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

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

Для реализации расчетов и построения графиков использовались среды Compaq Visual Fortran 6.5 и Maple15.

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

СОДЕРЖАНИЕ

фурье моделирование интегральный

ВВЕДЕНИЕ

1. Основные соотношения теории оболочек и теории упругости

1.1 Общие понятия теории упругости

1.2 Некоторые сведения об оболочках

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

2.1 Задача для покрытия

2.2 Задача для упругой полосы

3. Применение интегрального преобразования Фурье к решению задачи для полосы с покрытием

3.1 Определение и свойства преобразования Фурье

3.2 Применение интегрального преобразования Фурье в задаче для полосы

3.3 Применение интегрального преобразования Фурье в задаче для покрытия

4. Определение интегральных характеристик напряжений и перемещений под покрытием

5. Использованные алгоритмы и результаты вычислительных экспериментов

5.1 Поиск вещественных нулей и полюсов

5.2 Расчет перемещений и напряжений на границе покрытия и подложки

ЗАКЛЮЧЕНИЕ

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

ПРИЛОЖЕНИЕ А

ПРИЛОЖЕНИЕ Б

ПРИЛОЖЕНИЕ В

ПРИЛОЖЕНИЕ Г

ВВЕДЕНИЕ

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

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

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

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

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

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

 

1. Основные соотношения теории оболочек и теории упругости

 

.1 Общие понятия теории упругости

Теория упругости - раздел механики сплошных сред, изучающий деформации упругих твёрдых тел, их поведение при статических и динамических нагрузках.

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

В рамках теории линейной упругости исследуется напряженно-деформированное состояние упругой среды.

Механическое поведение идеально упругой среды во введенной ортогональной декартовой системе координат (х₁, х₂, х₃) описывается с помощью обобщенного закона упругости - закона Гука. В дальнейшем будут использованы обозначения х₁ = x, х₂ = y, х₃ = z. Напряжение, до которого справедлив закон Гука, называют пределом пропорциональности. Для изотропной среды можно считать, что линейные деформации зависят только от нормальных напряжений, а угловые - только от касательных:

,

,

,

, , .

Здесь , ,  - линейные относительные деформации; , ,  - нормальные напряжения; , ,  - угловые деформации (деформации сдвига); , ,  - касательные напряжения; Е - модуль упругости первого рода (модуль Юнга);  - модуль упругости второго рода (модуль сдвига); n - коэффициент поперечной деформации Пуассона. Они связаны между собой линейной зависимостью

Механическое состояние упругого тела, занимающего в начальном состоянии известный объем V с ограничивающей поверхностью S, характеризуется компонентами тензора деформации e_ij и тензора напряжений s_ij. Перемещения в точках тела, под действием заданной системы поверхностных и объемных сил, описываются вектором перемещений u={u₁, u₂, u₃} и представляет собой непрерывные и однозначные функции координат и времени. Механическое состояние упруго тела характеризуется также вектором напряжений t={t₁, t₂, t₃}, возникающих в упругом теле на некоторой элементарной площадке с нормалью . Вектор t выражается через компоненты тензора напряжений t_i=s_ijn_j, i=1,2,3.

В линейной теории упругости перемещения, деформация и напряжения связаны уравнениями движения [1]

, i=1,2,3,…,

соотношениями обобщенного закона Гука

, i,j=1,2,3,…

и геометрическими соотношениями Коши

, i,j=1,2,3.

Здесь  - вектор распределения объемных сил, x={x₁, x₂, x₃} - радиус-вектор заданной точки, r(х) и  заданное распределение плотности и коэффициентов упругости материала. У однородных тел коэффициенты упругости  и плотность r - константы.

В изотропном случае, когда упругие свойства одинаковы во всех направлениях, закон Гука имеет вид

,        (1.1.1)

где ,  - относительное изменение объема,  - символ Кронекера.

Вектор напряжений t, учитывая соотношения (1.1), выражается через перемещения t=Tu, где Т - линейный дифференциальный оператор напряжений. В изотропном случае

.

Колебания однородной упругой среды (в перемещениях) описываются уравнениями Ляме

,        (1.1.2)

где l, m - константы Ляме, r - плотность среды, вектор перемещений .

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

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

,

где

.

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

, i=.

Здесь q_i , i=- плотность заданных поверхностных сил.

При заданных на границе тела перемещениях имеем условия

,,

где j_i(x) - заданные на поверхности S функции.

Если на части поверхности S_s заданы напряжения (плотность поверхностных сил), а на части S_u - перемещения, граничные условия примут вид

, , .

Для свободной поверхности S_l (когда отсутствуют внешние напряжения и не заданы кинематические условия) граничные условия имеют вид

.

Встречаются и другие типы граничных условий.

В зависимости от вида граничных условий различают три типа основных задач теории упругости:

-    Основная задача первого типа состоит в определении компонент напряжений s_ij внутри области V, занятой телом, и компонент u_i(x) вектора перемещения точек внутри области V и точек поверхности S тела по заданным объемным силам F и поверхностным силам q_i.

-       Основная задача второго типа состоит в определении перемещений u_i(x) точек внутри области V и компонент напряжений s_ij по заданным силам F и по заданным перемещениям u_S на поверхности тела.

-       Основная задача третьего типа (смешанная задача) состоит в том, что по заданным поверхностным силам на одной части поверхности S_s и по заданным перемещениям на другой части поверхности тела S_u, требуется определить компоненты напряжений s_ij(x) и перемещений u_i(x). Получив решение данной задачи, можно определить, в частности, усилия связей на , чтобы реализовать заданные перемещения  на этой поверхности, а также можно вычислить перемещения  точек поверхности .

 

.2 Некоторые сведения об оболочках

Тонкой оболочкой называется тело, одно из измерений которого (толщина оболочки) значительно меньше двух других [2, 3]. Толщину оболочки будем предполагать постоянной. Поверхность, делящая пополам толщину оболочки, называется срединной поверхностью. Тонкая оболочка, срединная поверхность которой представляет собой плоскость, называется тонкой пластинкой.

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

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

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

Уравнение поперечных колебаний пластины, учитывающее влияние инерции вращения и деформации поперечного сдвига имеет вид [4]

.   (1.2.1)

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

Из уравнения (1.2.1) как частный случай вытекает уравнения колебаний пластины с учетом сдвига

.

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

.

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

 

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

Рассматривается плоская задача об установившихся колебаниях упругой полосы с покрытием. Колебания предполагаются гармоническими установившимися с частотой . Ось  введенной декартовой системы координат совпадает с верхней границей полосы (рисунок 1).

Рисунок 1 - Упругая полоса с покрытием

Предполагается, что материалы покрытия и подложки являются линейными однородными и изотропными.

 

.1 Задача для покрытия

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

          (2.1.1)

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

При этом заданы следующие граничные условия и условия сопряжения пластины и упругой подложки:

,

.

Далее рассматриваем установившийся (с частотой ) режим колебаний, т.е.

,

,

.

В дальнейшем, в силу линейности задачи, экспоненциальный множитель  можно опустить. При этом уравнение (2.1.1) примет вид

(2.1.2)

Здесь использовано обозначение .

 

.2 Задача для упругой полосы

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

Для установившихся гармонических колебаний инерционный член  в (1.1.2) принимает вид - . Тогда уравнение Ляме для упругой полосы запишется в виде

.

Или в покомпонентной записи

,    (2.2.1)

.  (2.2.2)

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

,

где  - вектор перемещения слоя.

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

.  (2.2.3)

Напряжение на границе покрытия при  в рассматриваемой задаче описываются вектором напряжений

.     (2.2.4)

На бесконечности требуется выполнение условий:

 при ,

и условий излучения. В качестве условия излучения выбран принцип предельного поглощения [5,6], когда в качестве решения задачи для идеально упругой среды берется равномерный предел решения соответствующей задачи для вязкоупругой среды (среды с поглощением) при стремлении вязкости к нулю.

Итак, относительно неизвестных перемещений  имеем краевую задачу (2.2.1) - (2.2.4).

 

3. Применение интегрального преобразования Фурье к решению задачи для полосы с покрытием

 

.1 Определение и свойства преобразования Фурье

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

Если функция , определена на , то ее трансформанта Фурье определяется по формуле

.     (3.1.1)

Ядро может иметь вид , но это не принципиально.

Комплексное преобразование Фурье существует далеко не для всех функций. Очевидно, для существования преобразования Фурье для функции  необходимо, чтобы она была определена на всей числовой оси  и .

Если для функции  существует преобразование Фурье (3.1.1), то обращение его можно определить как

.       (3.1.2)

Преобразования Фурье обладает следующими свойствами [7]:

а)     

б)     

в)     

 

.2 Применение интегрального преобразования Фурье в задаче для полосы

Геометрия задачи позволяет применить преобразование Фурье по переменной  к уравнениям (2.2.1), (2.2.2) и граничным условиям (2.2.3), (2.2.4).

В результате получим систему дифференциальных уравнений

         (3.2.1)

Штрихом здесь и далее обозначается производная по , заглавными буквами обозначены Фурье образы соответствующих компонент амплитуд вектора перемещений

С учетом условий, при  получаем

   (3.2.2)

При

      (3.2.3)

Используя для решения системы (3.2.1) с граничными условиями (3.2.2), (3.2.3) алгоритм, изложенный в [8], выразим искомые перемещения на границе пластины-покрытия и упругой полосы через заданную нагрузку и получим

       (3.2.4)

Матрица  имеет вид , где

, ,

, , , ,

,

Здесь  - плотность, модуль сдвига и коэффициент Пуассона среды полосы,

,

Используя эти соотношения можно построить функциональные соотношения на границе покрытия и подложки.

 

.3 Применение интегрального преобразования Фурье в задаче для покрытия

Аналогично изложенному в 3.2, можно применить преобразование Фурье по переменной x к уравнению (2.1.2).

Тогда уравнение (2.1.2) в трансформантах Фурье можно представить

  (3.3.1)

где  

Иначе можно записать

.  (3.3.2)

где

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

4. Определение интегральных характеристик напряжений и перемещений под покрытием

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

Для упругой полосы было получено соотношение (3.2.4)

Аналогично для покрытия (3.3.2)

Учитывая условия сопряжения на границе покрытия и упругой подложки, получим соотношение для напряжения на границе слоя и покрытия:

.         (4.1.1)

Где  , , , , ,

,

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

(4.4.2)

где  и  описываются соотношениями, приведенными выше.

 

5. Использованные алгоритмы и результаты вычислительных экспериментов

 

.1 Поиск вещественных нулей и полюсов

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

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

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

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

Отдельно выделяются вещественная и мнимая части функции с помощью функций Real и Aimag. Нули и той и другой частей ищутся на интервале, для каждого конкретного значения . Причем  изменяется от 0 до 10 с шагом, равным 0.05, в то время как для  он равен 0.01. Результаты работы программы выводятся в отдельные файлы.

Для поиска нулей использовалась стандартная подпрограмма библиотеки IMSL для поиска вещественных корней вещественной функции ZBREN [9]. Подпрограмма ZBREN применяет метод бисекций, линейную интерполяцию (метод секущих) и обратную квадратичную интерполяцию. На каждом шаге принимается решение, какой из трех методов будет использован для вычисления следующего приближения. Также IMSL использует для определения заданного числа вещественных корней метод Мюллера. Последний, вдобавок, применяется и при вычислении комплексных корней. (Данный метод является трехшаговым и основан примерно на тех же идеях, что и метод обратной квадратичной интерполяции а также метод секущих. Геометрически, на каждом шаге поиска очередного приближения корня используется три точки и строится парабола, проходящая через эти точки. В качестве следующего приближения берётся точка пересечения параболы и оси ).

Суть алгоритма заключается в исследовании отрезка, на котором. Если найден такой промежуток, то запускается стандартная подпрограмма ZBREN, которая имеет вызов:

Call ZBREN(f, errabs,errel,a,b,maxfn).

Параметры подпрограммы ZBREN:

-       Пользвательская функция: ;

-       Входные данные: errabs,errel;

-       Входные/выходные данные: a,b,maxfn.

Здесь errabs - первый критерий завершения вычислений. Корень b принимается, если . Причем, можно задать. Далее, errel - относительная ошибка (второй критерий завершения вычислений). Корень принимается, если относительная разница между приближениями, найденными в двух последовательных аппроксимациях, меньше errel.

Параметры а и b - начало и конец отрезка, на котором выполняется поиск вещественного корня. На выходе они отличны от начального значения, причем в b содержится приближенное значение корня функции f.

Параметр maxfn - на входе задает максимально допустимое число вычислений функции f, на выходе содержит число реальных вызовов функции f.

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

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

Рисунок 2 - Кривые нулей и полюсов функции

На рисунке 2 приведены кривые нулей и полюсов функции  (красным цветом обозначены нули элемента, черным - полюса) при следующих параметрах задачи: ; ; ; ; ; ; ; .

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

Рисунок 3 - Кривые нулей и полюсов функции для значений

; ; ; ; ; ; ;

Рисунок 4 - Кривые нулей и полюсов функции для значений

; ; ; ; ; ; ;

Рисунок 5 - Кривые нулей и полюсов функции для значений

; ; ; ; ; ; ;

Рисунок 6 - Кривые нулей и полюсов функции для значений

; ; ; ; ; ; ;

Рисунок 7 - Кривые нулей и полюсов функции для значений

; ; ; ; ; ; ;

Из рисунков 2 - 7 видно, что наблюдается чередование вещественных нулей и полюсов, с увеличением частоты число вещественных полюсов растет.

При построении графиков варьировались параметры, характеризующие толщину полосы и покрытия, а также их плотность.

 

5.2 Расчет перемещений и напряжений на границе покрытия и подложки

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

Амплитуды напряжений и перемещений могут быть выражены через их Фурье-образы следующим образом:

,

,

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

Для вычисления перемещений была создана программа на языке Fortran, которая производит обратное преобразование Фурье. Использовалась функция DQDAGS из библиотеки IMSL [9], вычисляющая значение интеграла от функции, имеющей конечное число особенностей на заданном промежутке. Она использует адаптивную схему, уменьшающую абсолютную ошибку. Подпрограмма делит отрезок  на подынтервалы и использует 21-точеченое правило Гаусса-Кронрода для оценки интеграла на каждом подынтервале.

DQDAGS имеет вызов:

Call DQDAGS(F,a,b,Errabs,Errel,Res,Errest).

Входными параметрами DQDAGS являются:

-       F - пользовательская функция, интеграл от которой должен быть вычислен;

-       а - нижняя граница интегрирования.

-       b - верхняя граница интегрирования;

-       Errabs - желаемая абсолютная ошибка;

-       Errel - желаемая относительная ошибка.

Выходные параметры:

-       Res - значение вычисленного интеграла;

-       Errest - оценка величины абсолютной ошибки.

Так как подынтегральная функция является комплексной, то вычисления ведутся для вещественной и мнимой части.

На рисунках 8 - 16 представлены графики вещественной части вертикальных смещений пластины для различных значений частоты  и для следующих параметров задачи: ; ; ; ; ; ; ; .

Рисунок 8 - Вещественная часть вертикальных смещений u при =1,9

Рисунок 9 - Вещественная часть вертикальных смещений u при =2,4

Рисунок 10 - Вещественная часть вертикальных смещений u при =3,4

Рисунок 11 - Вещественная часть вертикальных смещений u при =4,4

 

Рисунок 12 - Вещественная часть вертикальных смещений u при =5,4

Рисунок 13 - Вещественная часть вертикальных смещений u при =6,6

 

Рисунок 14 - Вещественная часть вертикальных смещений u при =7,3

Для вычисления напряжений также была написана программа на языке Fortran, которая производит обратное преобразование Фурье. Принцип ее работы аналогичен предыдущему. Ниже приведены графики, построенные по результатам работы программы для различных значений частоты . Значения вычислялись при следующих значениях параметров: ; ; ; ; ; ; ; .

Рисунок 17 - Вещественная часть амплитуды напряжений q при 𝜔=1,4

Рисунок 18 - Вещественная часть амплитуды напряжений q при =2,4

Рисунок 19 - Вещественная часть амплитуды напряжений q при =3,6

Рисунок 20 - Вещественная часть амплитуды напряжений q при =4,2

Рисунок 21 - Вещественная часть амплитуды напряжений q при =5,3

Рисунок 22 - Вещественная часть амплитуды напряжений q при =6,6

Рисунок 23 - Вещественная часть амплитуды напряжений q при =7,6

Рисунок 24 - Вещественная часть амплитуды напряжений q при =8,3

Рисунок 25 - Вещественная часть амплитуды напряжений q при  =9,3

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

 

ЗАКЛЮЧЕНИЕ

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

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

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

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

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

 

1            Демидов, С.П. Теория упругости - М.: Высшая школа, 1979. - 432с.

2       Александров, В.М. Контактные задачи для тел с тонкими покрытиями и прослойками / В.М. Александров, С.М. Мхитарян - М.: Наука, 1983. - 488 с.

3       Вольмир, А.С. Нелинейная динамика пластинок и оболочек - М.: Наука, 1972. - 432с.

4       Григолюк, Э.И. Механика твердых деформируемых тел. Т.5 Неклассическая теория колебаний стержней, пластин и оболочек / Э.И. Григолюк, И.Т. Селезнев - М.: ВИНИТИ, 1973. - 267с.

         Ворович, И.И. Неклассические смешанные задачи теории упругости / И.И. Ворович, В.М. Александров, В.А. Бабешко - М.: Наука, 1974. - 456 с.

         Ворович, И.И. Динамика массивных тел и резонансные явления в деформируемых средах / И.И. Ворович, В.А. Бабешко, О.Д. Пряхина - М.: Научный мир, 1989. - 246с.

         Снеддон, И. Преобразования Фурье - М.: Иностранная литература, 1955. - 668с.

         Глушков, Е.В. Интегральные преобразования в задачах теории упругости / Е.В. Глушков, Н.В. Глушкова - Кубанский государственный университет - М.: Краснодар, 1990. - 72с.

         Бартеньев, О.В. Фортран для профессионалов. Математическая библиотека IMSL. Выпуск 2 - М.: Диалог-МИФИ, 2001. - 309с.

10           Бартеньев, О.В. Современный Фортран - М.: Диалог-МИФИ, 2000. - 449 с.

11     Бабешко, В.А. Динамика неоднородных линейно-упругих сред / В.А. Бабешко, Е.В. Глушков, Ж.Ф. Зинченко - М.: Наука, 1989. - 344с.

12     Варвак, П.М. Справочник по теории упругости / П.М. Варвак, А.Ф. Рябов - Киев: Будiвельник, 1971. - 418с.

         Ворович, И.И. Механика контактных взаимодействий / И.И. Ворович, В.М. Александров - М.: Физматлит, 2001. - 672с.

         Павлова, А.В. Применение интегральных преобразований к решению задач для уравнений в частных производных / А.В. Павлова, С.Е. Рубцов, А.В. Смирнова - Кубанский государственный университет - М.: Краснодар, 2003. - 28с.

ПРИЛОЖЕНИЕ А

Листинг программы на языке Fortran. Поиск нулей и полюсов

program MyDip

dfimsl

real(4)::omega,alpha/global/omega(4), external:: fr(4), external:: fi(4), external:: znamr(4), external:: znami(4)::errabs=0.0,errrel=0.0001 (4)::a,b,nulr,polusr,nuli,polusi,nul,polus(4)::maxfn,k1=0,k2=0,k3=0,flag1,flag11,flag2,flag22,flag3,flag33,flag4,flag44,ff1,ff2

(4, file = 'MOmegaNuli.txt')(5, file = 'MNuli.txt')(6, file = 'MOmegaPolusa.txt')(7, file = 'MPolusa.txt')

!главный циклomega=0.2, 10.0, 0.005

         do alpha=0.0,10.0,0.001

         flag1=0;flag11=0;flag2=0;flag22=0;flag3=0;flag33=0;flag4=0;flag44=0;ff1=0;ff2=0

!поиск нуля вещественной части

                   a=alpha

                   b=alpha+0.001

                   if ((fr(a)*fr(b))<=0.0000) then

                            maxfn=100

                            flag1=1

                            call ZBREN(fr,errabs,errrel,a,b,maxfn)

                                      if (maxfn<100) then

                                      flag11=1

                                      nulr=b

                                      endif

                   end if

!поиск нуля мнимой части

                   a=alpha

                   b=alpha+0.001

                   if ((fi(a)*fi(b))<=0.0000) then

                            maxfn=100

                            flag3=1

                            call ZBREN(fi,errabs,errrel,a,b,maxfn)

                                      if (maxfn<100) then

                                      flag33=1

                                      nuli=b

                   end if

!поиск полюса вещественной части

                   !для функции знаменателя

                   a=alpha

                   b=alpha+0.001

                   if ((znamr(a)*znamr(b))<=0.0000) then

                            maxfn=100

                            flag2=1

                            call ZBREN(znamr,errabs,errrel,a,b,maxfn)

                                      if (maxfn<100) then

                                      flag22=1

                                      polusr=b

                                      endif

                   end if

!поиск полюса мнимой части

                   !для функции знаменателя

                   a=alpha

                   b=alpha+0.001

                   if ((znami(a)*znami(b))<=0.0000) then

                            maxfn=100

                            flag4=1

                            call ZBREN(znami,errabs,errrel,a,b,maxfn)

                                      if (maxfn<100) then

                                      flag44=1

                                      polusi=b

                                      endif

                   end if

                   if ((flag11==1) .and. (flag33==1)) then

                            if (abs(nulr - nuli)<0.0001) then

                            nul = nulr

                            ff1=1

                            endif

                   endif 

                   if ((flag22==1) .and. (flag44==1)) then

                            if (abs(polusr - polusi)<0.0001) then                                   

                            polus = polusr

                            ff2=1

                            endif

                   endif

                   if ((ff1==1) .and. (ff2==1) .and. (nul /= polus) ) then

                            write (4,*) omega !вывод омеги в файл

                            write (5,'(3f14.8)'), nul !форматный вывод корня

                            k1=k1+1

                            write (6,*) omega !вывод омеги в файл

                            write (7,'(3f14.8)'), polus !форматный вывод корня

                            k2=k2+1

                   elseif (ff1==1 .and. (ff2==0)) then

                                      write (4,*) omega !вывод омеги в файл

                                      write (5,'(3f14.8)'), nul !форматный вывод корня

                                      k1=k1+1

                            elseif (ff2==1 .and. (ff1==0)) then

                                      write (6,*) omega !вывод омеги в файл

                                      write (7,'(3f14.8)'), polus !форматный вывод корня

                                      k2=k2+1

                   endif

                  

         end dodo(*,*)k1(*,*)k2

program MyDip

fr(alpha)     !функция

(4), PARAMETER :: i=(0.0,1.0)      !мнимая единица(4) :: fr,alpha

         complex(4)::eps, eps1, eps2, eps3, ka2, del11, del21, k11, k12, k22, psi, gam1, sig1, sig2    !они комплексные

        

!        real::mu1=1, mu2=2.6, nu1=0.16, nu2=0.25, h=0.1, h1=1, ro2=3.9, ro1=1 , omega !старые параметры

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega         !1 набор

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !2 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega        !3 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !4 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=1, ro1=1 , omega  !5 набор

         real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=1, ro1=1 , omega        !6 набор

         common /global/ omega

         eps1=(mu2*(h**3))/(6*(1-nu2))

         eps2=((h**3)/12)*ro2*(omega**2)

         eps3=h*ro2*(omega**2)

         ka2=sqrt((ro1*(omega**2))/mu1)

         eps=(1-2*nu1)/(2-2*nu1)

         psi=eps1*(alpha**4)-eps2*(alpha**2)-eps3

         gam1=(alpha**2)-0.5*(ka2**2)

         sig2=sqrt((alpha**2)-(ka2**2))

         sig1=sqrt((alpha**2)-eps*(ka2**2))

k11=(sig2*(ka2**2)*(sig1*sig2*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)-(alpha**2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)))/(-4)

k12=(sig1*sig2*(gam1+(alpha**2))*(((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)-1)-(gam1*(alpha**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

k22=(sig1*(ka2**2)*(sig2*sig1*((exp(sig2*h1)+exp(-(sig2*h1)))/2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)-(alpha**2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)*((exp(sig1*h1)+exp(-(sig1*h1)))/2)))/(-4)

         del21=(sig1*sig2*((1/4)*(ka2**4)-((gam1+(alpha**2))**2)+((1/4)*(ka2**4)+((gam1+(alpha**2))**2))*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2))-(alpha**2)*((gam1**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

         del11=((exp(sig2*h1)+exp(-(sig2*h1)))/2)

         fr=real(k22/del21)

function fr

fi(alpha)     !функция

(4), PARAMETER :: i=(0.0,1.0)      !мнимая единица(4) :: fi,alpha

         complex(4)::eps, eps1, eps2, eps3, ka2, del11, del21, k11, k12, k22, psi, gam1, sig1, sig2    !они комплексные

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega         !1 набор

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !2 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega        !3 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !4 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=1, ro1=1 , omega  !5 набор

         real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=1, ro1=1 , omega        !6 набор

         common /global/ omega

         eps1=(mu2*(h**3))/(6*(1-nu2))

         eps2=((h**3)/12)*ro2*(omega**2)

         eps3=h*ro2*(omega**2)

         ka2=sqrt((ro1*(omega**2))/mu1)

         eps=(1-2*nu1)/(2-2*nu1)

         psi=eps1*(alpha**4)-eps2*(alpha**2)-eps3

         gam1=(alpha**2)-0.5*(ka2**2)

         sig2=sqrt((alpha**2)-(ka2**2))

         sig1=sqrt((alpha**2)-eps*(ka2**2))

k11=(sig2*(ka2**2)*(sig1*sig2*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)-(alpha**2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)))/(-4)

k12=(sig1*sig2*(gam1+(alpha**2))*(((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)-1)-(gam1*(alpha**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

k22=(sig1*(ka2**2)*(sig2*sig1*((exp(sig2*h1)+exp(-(sig2*h1)))/2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)-(alpha**2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)*((exp(sig1*h1)+exp(-(sig1*h1)))/2)))/(-4)

         del21=(sig1*sig2*((1/4)*(ka2**4)-((gam1+(alpha**2))**2)+((1/4)*(ka2**4)+((gam1+(alpha**2))**2))*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2))-(alpha**2)*((gam1**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

         del11=((exp(sig2*h1)+exp(-(sig2*h1)))/2)

         fi=aimag(k22/del21)

function fi

znamr(alpha)       !функция

(4), PARAMETER :: i=(0.0,1.0)      !мнимая единица(4) :: znamr,alpha

         complex(4)::eps, eps1, eps2, eps3, ka2, del11, del21, k11, k12, k22, psi, gam1, sig1, sig2    !они комплексные

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega         !1 набор

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !2 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega        !3 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !4 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=1, ro1=1 , omega  !5 набор

         real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=1, ro1=1 , omega        !6 набор

         common /global/ omega

         eps1=(mu2*(h**3))/(6*(1-nu2))

         eps2=((h**3)/12)*ro2*(omega**2)

         eps3=h*ro2*(omega**2)

         ka2=sqrt((ro1*(omega**2))/mu1)

         eps=(1-2*nu1)/(2-2*nu1)

         psi=eps1*(alpha**4)-eps2*(alpha**2)-eps3

         gam1=(alpha**2)-0.5*(ka2**2)

         sig2=sqrt((alpha**2)-(ka2**2))

         sig1=sqrt((alpha**2)-eps*(ka2**2))

k11=(sig2*(ka2**2)*(sig1*sig2*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)-(alpha**2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)))/(-4)

k12=(sig1*sig2*(gam1+(alpha**2))*(((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)-1)-(gam1*(alpha**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

k22=(sig1*(ka2**2)*(sig2*sig1*((exp(sig2*h1)+exp(-(sig2*h1)))/2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)-(alpha**2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)*((exp(sig1*h1)+exp(-(sig1*h1)))/2)))/(-4)

         del21=(sig1*sig2*((1/4)*(ka2**4)-((gam1+(alpha**2))**2)+((1/4)*(ka2**4)+((gam1+(alpha**2))**2))*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2))-(alpha**2)*((gam1**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

         del11=((exp(sig2*h1)+exp(-(sig2*h1)))/2)

         znamr=real((k22/del21)*psi+1)

function znamr

znami(alpha)        !функция

(4), PARAMETER :: i=(0.0,1.0)      !мнимая единица(4) :: znami,alpha

         complex(4)::eps, eps1, eps2, eps3, ka2, del11, del21, k11, k12, k22, psi, gam1, sig1, sig2    !они комплексные

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega         !1 набор

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !2 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega        !3 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !4 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=1, ro1=1 , omega  !5 набор

         real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=1, ro1=1 , omega        !6 набор

         common /global/ omega

         eps1=(mu2*(h**3))/(6*(1-nu2))

         eps2=((h**3)/12)*ro2*(omega**2)

         eps3=h*ro2*(omega**2)

         ka2=sqrt((ro1*(omega**2))/mu1)

         eps=(1-2*nu1)/(2-2*nu1)

         psi=eps1*(alpha**4)-eps2*(alpha**2)-eps3

         gam1=(alpha**2)-0.5*(ka2**2)

         sig2=sqrt((alpha**2)-(ka2**2))

         sig1=sqrt((alpha**2)-eps*(ka2**2))

k11=(sig2*(ka2**2)*(sig1*sig2*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)-(alpha**2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)))/(-4)

k12=(sig1*sig2*(gam1+(alpha**2))*(((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2)-1)-(gam1*(alpha**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

k22=(sig1*(ka2**2)*(sig2*sig1*((exp(sig2*h1)+exp(-(sig2*h1)))/2)*((exp(sig1*h1)-exp(-(sig1*h1)))/2)-(alpha**2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2)*((exp(sig1*h1)+exp(-(sig1*h1)))/2)))/(-4)

         del21=(sig1*sig2*((1/4)*(ka2**4)-((gam1+(alpha**2))**2)+((1/4)*(ka2**4)+((gam1+(alpha**2))**2))*((exp(sig1*h1)+exp(-(sig1*h1)))/2)*((exp(sig2*h1)+exp(-(sig2*h1)))/2))-(alpha**2)*((gam1**2)+(sig1**2)*(sig2**2))*((exp(sig1*h1)-exp(-(sig1*h1)))/2)*((exp(sig2*h1)-exp(-(sig2*h1)))/2))/2

         del11=((exp(sig2*h1)+exp(-(sig2*h1)))/2)

         znami=aimag((k22/del21)*psi+1)

end function znami

 

ПРИЛОЖЕНИЕ Б

Листинг программы на языке Fortran. Построение перемещений

Program zer1

use dfimsl/vab/ omega,x1(8),external::f     !,U3,U4(8)::omega,x1,alpha,alpha11,errabs,errel,res,res1,errest,Int1,Int2,a(8)::Int3(4)::i,j,maxfn,maxfn1,ff

!open(1,action='write',file='xf.txt',form='formatted',position='rewind',status='old')        !файл для перемещений(1,action='write',file='yf.txt',form='formatted',position='rewind',status='old')        !файл дл обращений

=0.0d0=0.01d0

                           

                  

                            do x1=0d0,10d0,0.1

                                       

                                       omega=9.3

                                       !omega=2.100009

                                       

                                       do a=-30D0,10D0,1.0

                                       alpha=a

                                                                  alpha11=alpha+0.1

                                                                  call dqdags(f,alpha,alpha11,errabs,errel,res,errest)

                                                                  Int1=Int1+res

                                       end do

                                       

                                       write(1,'(f10.5,2x,15(f8.5,1x))') x1,Int1

                                       Int1=0.0

                            end do

close(1)program zer1

f(alpha)      !функция

/vab/omega,x1(8), PARAMETER :: I=(0.0,1.0) !мнимая единица(8) :: f,alpha,x1

         complex(8)::eps, eps1, eps2, eps3, ka2, del11, del21, k11, k12, k22, psi, gam1, sig1, sig2, elem    !они комплексные

         real(8)::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega         !1 набор

!        real::mu1=1, mu2=2.6, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !2 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=3.8, ro1=1 , omega        !3 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=3.8, ro1=1 , omega        !4 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.1, h1=1, ro2=1, ro1=1 , omega  !5 набор

!        real::mu1=1, mu2=5, nu1=0.3, nu2=0.3, h=0.01, h1=1, ro2=1, ro1=1 , omega        !6 набор

         eps1=(mu2*(h**3))/(6*(1-nu2))

         eps2=((h**3)/12)*ro2*(omega**2)

         eps3=h*ro2*(omega**2)

         ka2=cdsqrt(dcmplx((ro1*(omega**2))/mu1))

         eps=(1-2*nu1)/(2-2*nu1)

         psi=eps1*(alpha**4)-eps2*(alpha**2)-eps3

         gam1=(alpha**2)-0.5*(ka2**2)

         sig2=cdsqrt(dcmplx((alpha**2)-(ka2**2)))

         sig1=cdsqrt(dcmplx((alpha**2)-eps*(ka2**2)))

k11=(sig2*(ka2**2)*(sig1*sig2*((cdexp(sig1*h1)+cdexp(-(sig1*h1)))/2)*((cdexp(sig2*h1)-cdexp(-(sig2*h1)))/2)-(alpha**2)*((cdexp(sig1*h1)-cdexp(-(sig1*h1)))/2)*((cdexp(sig2*h1)+cdexp(-(sig2*h1)))/2)))/(-4)

k12=(sig1*sig2*(gam1+(alpha**2))*(((cdexp(sig1*h1)+cdexp(-(sig1*h1)))/2)*((cdexp(sig2*h1)+cdexp(-(sig2*h1)))/2)-1)-(gam1*(alpha**2)+(sig1**2)*(sig2**2))*((cdexp(sig1*h1)-cdexp(-(sig1*h1)))/2)*((cdexp(sig2*h1)-cdexp(-(sig2*h1)))/2))/2

k22=(sig1*(ka2**2)*(sig2*sig1*((cdexp(sig2*h1)+cdexp(-(sig2*h1)))/2)*((cdexp(sig1*h1)-cdexp(-(sig1*h1)))/2)-(alpha**2)*((cdexp(sig2*h1)-cdexp(-(sig2*h1)))/2)*((cdexp(sig1*h1)+cdexp(-(sig1*h1)))/2)))/(-4)

         del21=(sig1*sig2*((1/4)*(ka2**4)-((gam1+(alpha**2))**2)+((1/4)*(ka2**4)+((gam1+(alpha**2))**2))*((cdexp(sig1*h1)+cdexp(-(sig1*h1)))/2)*((cdexp(sig2*h1)+cdexp(-(sig2*h1)))/2))-(alpha**2)*((gam1**2)+(sig1**2)*(sig2**2))*((cdexp(sig1*h1)-cdexp(-(sig1*h1)))/2)*((cdexp(sig2*h1)-cdexp(-(sig2*h1)))/2))/2

         del11=((cdexp(sig2*h1)+cdexp(-(sig2*h1)))/2)

!        f=real(((k22/del21)/((k22/del21)*psi+1))*cdexp(-I*alpha*x1))       !для перемещений

         f=real(((1)/((k22/del21)*psi+1))*cdexp(-I*alpha*x1))  !для напряжений

end function f

ПРИЛОЖЕНИЕ B

Листинг программы в среде Maple15 (построение графиков нулей и полюсов).

(stats);

:= readdata(`c:\\MOmegaNuli.txt`, 1);:= readdata(`c:\\MNuli.txt`, 1);:= readdata(`c:\\MOmegaPolusa.txt`, 1);:= readdata(`c:\\MPolusa.txt`, 1);[scatterplot](x, y);[scatterplot](z, q);

ПРИЛОЖЕНИЕ Г

Листинг программы в среде Maple15 (построение графиков перемещений и напряжений).

(stats);

[listplot](readdata(`c:\\xf19n.txt`, float, 2), color = red, labels = ["x", "u"], labelfont = [TIMES, BOLDITALIC, 11]);[listplot](readdata(`c:\\yf14.txt`, float, 2), color = red, labels = ["x", "u"], labelfont = [TIMES, BOLDITALIC, 11]);