Оглавление
Введение ................................................................................................................... 4
Глава 1. Обзор научной литературы по вопросам укладки кабеля под водой . 9
1.1 Глобальность задачи об укладке кабеля под водой ................................... 9
1.2 Описание современных способов укладки кабеля под водой ................ 11
1.3 История развития конструкции коммуникационного подводного кабеля
............................................................................................................................. 13
1.4 История развития конструкции силового подводного кабеля ............... 15
1.5 Обзор современных работ по механике деформируемого кабеля ......... 17
Глава 2. Современное состояние математических моделей об укладке кабеля
и смежных задачах о движении кабелей под водой .......................................... 19
2.1 Аналитическое решение задачи Рауса-Аппеля ........................................ 19
2.2 История развития задачи об укладке кабеля ............................................ 22
2.3 Описание смежных задач о швартовании ................................................ 24
2.4 Описание смежных задач о буксировке и совместном движении
системы «судно – кабель – подводный аппарат» .......................................... 25
2.5 Выводы, описание возможных направлений исследования ................... 26
Глава 3. Математическое моделирование динамического равновесия кабеля
при его укладке под водой ................................................................................... 27
3.1 Численное конечно-разностное решение задачи Рауса-Аппеля ............ 27
3.2 Учет растяжимости кабеля согласно закону Гука ................................... 29
3.3 Учет подводного течения с профилем скорости, изменчивым по
глубине ............................................................................................................... 33
3.4 Верификация и апробация математической модели динамического
равновесия кабеля при его укладке под водой............................................... 36
3.5 Оценка формы и натяжения кабеля при его укладке вблизи газового
месторождения Ормен Ланге ........................................................................... 43
3.6 Выводы ......................................................................................................... 45
Глава 4. Математическое моделирование движения кабеля при его укладке
под водой ................................................................................................................ 47
4.1 Описание математической модели движения кабеля при его укладке
под водой ............................................................................................................ 47
2
4.2 Описание интегрированной программы по определению
нестационарного пространственного нагружения кабеля ............................ 50
4.3 Верификация конечно-элементной модели на примере задачи о
натяжении якорной цепи .................................................................................. 52
4.4 Учет микроструктуры кабеля .................................................................... 54
4.4.1 Описание постановки задачи для определения эффективных
коэффициентов жесткости кабеля ............................................................... 54
4.4.2 Описание конструкции силовых подводных кабелей ...................... 56
4.4.3 Описание конечно-элементных моделей силовых подводных
кабелей............................................................................................................ 59
4.4.4 Определение эффективных коэффициентов жесткости кабеля ...... 62
4.4.5 Определение влияния микроструктуры кабеля на характеристики
укладки кабеля ............................................................................................... 65
4.5 Учет рассогласованной укладки кабеля, вызванной сматыванием кабеля
с барабана лебедки ............................................................................................ 67
4.5.1 Определение статического распределения реакций схемы
крепления барабана лебедки ........................................................................ 68
4.5.2 Определение момента силы, необходимого для запуска лебедки .. 73
4.5.3 Определение момента силы и потребляемой мощности
электродвигателя при вращении барабана лебедки .................................. 76
4.5.4 Определение влияния рассогласованного динамического движения
на характеристики укладки кабеля .............................................................. 82
4.6 Выводы ......................................................................................................... 83
Заключение ............................................................................................................ 85
Список литературы ............................................................................................... 86
Приложение А. Комплекс программ для определения формы и натяжения
провисающего участка кабеля при его укладке под водой .............................. 96
Приложение Б. Интегрированная программа по определению
нестационарного пространственного нагружения кабеля .............................. 119
3
Введение
Актуальность. Укладка кабеля под водой – это невероятно важный и
глобальный процесс, в который вовлечены многочисленные страны и
компании. Область применения результатов проложенного под водой кабеля
широка и охватывает такие отрасли как: передача телефонных разговоров,
видеозаписей и интернет данных; электроснабжение островов, морских
нефтяных
и
газовых
электростанций,
платформ,
ветровых,
океанологических
волновых
лабораторий;
и
военные;
приливных
научно-
исследовательские.
На сегодняшний день, прокладка подводного телекоммуникационного
кабеля под водой – это доминирующий и надежный способ установления
стабильных коммуникаций между континентами. Прокладка подводного
силового кабеля – это зачастую, единственный, способ осуществления
электроснабжения городов, станций и лабораторий отдаленных водной
преградой.
Укладка кабеля под водой – это не только важный глобальный процесс,
но и вместе с этим невероятно сложный процесс, ведь средняя глубина
океана достигает 3.7 километров, а протяженность кабельных линий может
достигать тысяч километров.
Работа над укладкой кабеля под водой продолжается и сегодня. Кроме
выполнения ремонтных работ и поддержания эксплуатационного состояния
уже проложенных кабельных систем, также ведутся активные работы по
построению новых. Опишем несколько кабельных систем, планируемых к
реализации в 2020 году. Корпорации Google и RTI осуществляют подготовку
к реализации кабельной системы Гонконг – Гуам (HG-G), длиной 3900 км и
соединяющей г. Ченкуанъоу Нью Таун с дер. Пити. Также, корпорация
Google анонсировала прокладку еще одной кабельной системы: США –
Франция (Dunant), длиной 6400 км и соединяющей г. Верджиния-Бич с
комунной Сент-Илер-де-Рье.
4
Цель исследования состоит в разработке цифровой модели процесса
укладки кабеля под водой, позволяющей учитывать многочисленные
физические явления, имеющие место при проведении реальных морских
работ.
Для достижения поставленной цели рассматриваются следующие
задачи исследования:
Анализ современных математических моделей об укладке кабеля, а также
смежных задач, тесно связанных использованием теории механики кабеля
и движением в сопротивляющейся среде;
Разработка математической модели динамического равновесия кабеля при
его укладке под водой с учетом:
растяжимости кабеля согласно закону Гука;
влияния подводного течения с профилем скорости, изменчивым по
высоте;
Разработка комплекса программ для определения формы и натяжения
провисающего участка кабеля;
Разработка математической модели движения кабеля при его укладке под
водой с учетом:
микроструктуры кабеля;
рассогласованной укладки кабеля, вызванной сматыванием кабеля с
барабана лебедки.
Разработка интегрированной программы нагружения кабеля.
Методы решения. Для решения поставленных задач применялись
различные методы и программное обеспечение. Для решения общих задач
использована конечно-разностная схема и итерационный многомерный
метод Ньютона для решения системы нелинейных дифференциальных
уравнений
на
языке
программирования
Matlab
MathWorks.
Также
применялось прямое конечно-элементное моделирование с использованием
явной схемы решения программного комплекса ABAQUS Explicit, а также
языка программирования Fortran. Для решения некоторых частных задач об
5
определении
статического
равновесия
деформируемых
использовалась неявная схема решения, реализованная в
систем
программных
комплексах ABAQUS Implicit и ANSYS. Для моделирования динамики
абсолютно жестких тел использовались методы GSTIFF и WSTIFF
программного комплекса ADAMS и его модуля Machinery. Для подготовки
геометрических моделей использовались системы автоматизированного
проектирования SolidWorks и Siemens NX.
Достоверность результатов исследования и сделанных выводов
подтверждаются
многочисленными
проведенными
верификационными
проверками, обоснованным применением современных программных систем,
а также адекватностью физического представления о процессе и объекте
исследования. Проведены проверки численных методов, а также проведено
сравнение с известным аналитическим решением Меркина.
Научная новизна состоит в следующем:
Впервые учтен профиль подводного течения для задачи Рауса-Аппеля;
Впервые учтена многокомпонентная гетерогенная микроструктура кабеля
в виде эффективных характеристик кабеля и рассогласованность скорости
схода кабеля, вызванная сматыванием с барабана лебедки;
Проведено исследование влияния скорости укладки кабеля, механических
свойств кабеля, подводного течения, микроструктуры и параметров
рассогласованного движения на форму и натяжения провисающего
участка кабеля.
На защиту выносятся следующие положения:
Цифровые математические модели процесса укладки кабеля под водой;
Программа
по
определению
нестационарного
пространственного
нагружения кабеля, интегрированная в конечно-элементный программный
комплекс;
Результаты определения эффективных коэффициентов жесткости кабелей
различных марок;
Результаты влияния скорости укладки кабеля, механических свойств
6
кабеля, подводного течения, микроструктуры и рассогласованного
движения на форму и натяжение провисающего участка кабеля.
В первой главе диссертации приведены общие сведения о процессе
исследования – укладке кабеля под водой, а также об объекте исследования –
подводном кабеле.
Вторая
современного
глава
диссертации
состояния,
математических
моделей
а
об
посвящена
также
истории
укладке
детальному
открытия
кабеля.
описанию
и
развития
Приведено
описание
современных математических моделей смежных задач о швартовании,
буксировке и движению системы «судно – кабель – подводный аппарат».
Схожесть этих классов задач обусловлена использованием теории механики
нити и учетом гидродинамических сил сопротивления воды.
В третьей главе диссертации рассматривается математическое
моделирование динамического равновесия кабеля при его укладке под водой
в состоянии т.н. «кажущегося покоя». Учтены растяжимость кабеля согласно
закону Гука и влияние подводного течения с профилем скорости,
изменчивым по глубине. Проведены численные проверки и сравнение с
известным аналитическим решением с использованием разработанного
автором комплекса программ. Рассмотрено влияние механических свойств
кабеля и согласованной скорости укладки на форму и натяжение
провисающего участка кабеля. Выполнена оценка формы и натяжения кабеля
вблизи газового месторождения Ормен Ланге.
В заключительной, четвертой главе диссертации рассматривается
математическое моделирование движения кабеля при его укладке под водой
с использованием разработанной автором интегрированной программы
нестационарного
описание
пространственного
интегрированной
действующей
в
специальных
нагружения
программы
условиях.
для
кабеля.
Приведено
формирования
Учтена
силы,
многокомпонентная
гетерогенная микроструктура кабеля в виде эффективных характеристик
7
кабеля и рассогласованность скорости схода кабеля, вызванная сматыванием
с барабана лебедки.
Научные результаты диссертации представлены 5 докладами:
дважды на конференции «Неделя науки СПбПУ» в 2015 и 2017 годах в
Санкт-Петербурге;
на конференции «Инженерные технологии MSC Software для высших
учебных заведений» в 2016 году в МАДИ в Москве;
на научном семинаре «Математическое моделирование процесса укладки
кабеля под водой и анализ работоспособности роторного механизма с
цевочной передачей» на кафедре «Механика и процессы управления» в
ИПММ СПбПУ в Санкт-Петербурге в 2019 году;
на съезде «XII Всероссийский съезд по фундаментальным проблемам
теоретической и прикладной механике» в 2019 году в Уфе.
Работа победила в конкурсе грантов комитета по науке и высшей
школы правительства Санкт-Петербурга для студентов и аспирантов вузов,
отраслевых академических институтов, расположенных на территории
Санкт-Петербурга в 2017 году.
Научные результаты диссертации представлены в 9 публикациях, в том
числе 3 публикации в рецензируемых научных изданиях из перечня ВАК по
научной
специальности
05.13.18
–
Математическое
моделирование,
численные методы и комплексы программ:
статья в журнале «Подводные исследования и робототехника» №1(27)
2019 года (импакт-фактор РИНЦ 2017 0.738);
свидетельство на программу для электронных вычислительных машин
2019 года;
статья в журнале «Морские интеллектуальные технологии» №3(45) том 3
2019 года (импакт-фактор РИНЦ 2017 0.197).
8
Глава 1. Обзор научной литературы по вопросам укладки кабеля под
водой
Первая глава диссертации посвящена описанию: глобальности задачи
об укладке кабеля под водой, современных способов укладки кабеля под
водой, истории развития конструкции подводных кабелей и обзору
современных работ по механике деформируемого кабеля.
1.1 Глобальность задачи об укладке кабеля под водой
Глобальная протяженность всех кабельных систем, проложенных под
водой, составляет более 1 млн. км. Этой длины достаточно, чтобы 25 раз
обогнуть Землю. Количество подводных кабельных систем превышает 350
[1]. Отметим также, что средняя глубина океана составляет около 3.7 км [2].
Для наглядности, представим современную мировую карту подводного
кабеля [3], на которой изображены коммуникационные и силовые подводные
кабели, проложенные по дну морей и океанов на рис. 1.1.
Рис. 1.1 Мировая карта проложенного под водой кабеля [3]
9
Прокладкой кабеля под водой активно занимаются такие страны, как:
Франция, Великобритания, США, Корея, ОАЭ, Германия и многие другие.
На рис. 1.2 приведены примеры суден-кабелеукладчиков [4].
Ile de Batz, Франция, Alcatel-Lucent
Submarine Networks Ltd., грузовместимость
5500 тонн кабеля [4,5]
Pacific Guardian, Великобритания, Global
Marine Systems, грузовместимость 1700
тонн кабеля [4,6]
Dependable, США, Tyco Electronics Subsea
Communications LLC, грузовместимость
5446 тонн кабеля [4,7]
Responder, Корея, KT Submarine,
грузовместимость 5598 тонн кабеля [4,8]
Maersk Connector, Германия, Maersk Supply
Service, грузовместимость 7000 тонн кабеля
[10]
Рис. 1.2 Судно-кабелеукладчик [4]
CS Maram, ОАЭ, E-Marine PJSC,
грузовместимость 2750 тонн кабеля [4,9]
10
1.2 Описание современных способов укладки кабеля под водой
Опишем три, хорошо известных, способа укладки кабеля под водой,
применяемых в различных ситуациях. Выбор способа укладки кабеля зависит
от многочисленных факторов, преимущественное влияние оказывают: режим
судоходства, условия укладки кабеля, профиль морского дна, длина
кабельной линии и многие другие [11].
На сегодняшний день, одним из наиболее распространенных способов
укладки кабеля под водой является укладка кабеля с использованием
подводного робота – телеуправляемого необитаемого подводного аппарата,
представленная на рис. 1.3 [11,12]. Подводный робот служит для контроля
провисающего участка кабеля при его укладке на морском дне. Также в
случае разрыва кабельных линий, подводный робот может использоваться
для проведения ремонтных работ.
Спуск подводного телеуправляемого
Укладка подводного силового кабеля с
необитаемого подводного аппарата в воду
использованием телеуправляемого
[12]
необитаемого подводного аппарата [11]
Рис. 1.3 Укладка кабеля под водой с использованием телеуправляемого необитаемого
подводного аппарата [11,12]
Альтернативным способом укладки кабеля является укладка кабеля под
водой путем его предварительного расположения на поверхности воды,
представленная на рис. 1.4 [11]. После прокладки кабеля, буйки,
поддерживающие его на поверхности, отсоединяются от кабеля.
11
Подводный кабель Basslink, соединяющий
Прокладка подводного силового кабеля
Австралию и Тасманию
судном C/S Skagerrak, Норвегия
Рис. 1.4 Укладка кабеля под водой путем предварительного расположения кабеля на
поверхности воды [11]
Третьим способом укладки кабеля под водой является укладка кабеля в
результате его сматывания с барабана лебедки, представленная на рис. 1.5.
При таком способе укладки, как правило, обеспечивается согласованная
скорость укладки кабеля (скорость движения судна совпадает со скоростью
схода кабеля) путем контролирования натяжения кабеля. Именно этому
способу укладки посвящена настоящая диссертация.
Прокладка высоковольтного подводного
силового кабеля судном C.S. Sovereign,
Великобритания
Рис. 1.5 Укладка кабеля под водой в результате его сматывания с барабана лебедки [11]
Спуск глубоководного оптоволоконного
кабеля для передачи информации в воду
12
1.3 История развития конструкции коммуникационного подводного кабеля
На
сегодняшний
день,
подводный
коммуникационный
кабель
осуществляет передачу миллионов телефонных разговоров в совместности с
бесчисленным количеством видеозаписей и интернет-данных [11,13].
Приведем краткое описание основных исторических событий укладки
подводных
коммуникационных
кабелей
и
событий,
определяющих
направление развития конструкции подводного коммуникационного кабеля к
его современному состоянию [12]:
1840 год: начало прокладки телеграфных кабелей через реки и гавани,
рассчитанных на ограниченный срок службы;
1843-1845 года: использование гуттаперчи, привезенной британцами из
Малайзии, в качестве электрической изоляции кабеля, продлевающей
срок службы кабеля;
1850 год: первая прокладка международного телеграфного кабеля
между Великобританией и Францией, с последующим усилением в
1851 году;
1858 год: первая прокладка трансатлантического кабеля между
Ирландией и о. Ньюфаундленд британским пароходом «Грейт Истерн»,
проработавшего 26 дней до поломки. Впоследствии был перепроложен
в 1866 году;
1884 год: первая прокладка телефонного кабеля между г. СанФранциско и г. Окленд;
1920-1930 года: использование коротковолнового радио в качестве
альтернативного способа передачи информации, однако все же
имеющего свои недостатки, обусловленные ограничением мощности и
чувствительности к атмосферным условиям;
1956 год: применение ретрансляторов в первом трансатлантическом
кабеле
TAT-1
способствовало
началу
быстрой
и
надежной
трансокеанской связи;
1961 год: возникновение высококачественной глобальной сети связи;
13
1986 год: первый международный оптоволоконный кабель между
Бельгией и Великобританией;
1988 год: первая трансокеанская оптоволоконная кабельная система,
соединяющая США, Великобританию и Францию.
Конструкция подводного коммуникационного кабеля различается в
зависимости от используемых материалов и непосредственных условий
укладки кабеля. Приведем типичную конструкцию современного подводного
коммуникационного кабеля [12]:
Рис. 1.6 Конструкция современного подводного коммуникационного кабеля [12]
Центральной конструктивной составляющей современного подводного
коммуникационного кабеля являются оптические волокна из кварцевого
стекла в качестве среды для передачи информации различного рода. Иногда,
для обеспечения жесткости расположения оптических волокон, последние
фиксируются и сепарируются трубками из полиэтилена или стекловолокна.
Следующим конструктивным слоем (через гидрофобный заполнитель из
полиэтилена) является сепаратор-экран, представляющий из себя, как
правило, медную трубку. Заключительный конструктивным слоем (снова
минуя гидрофобный заполнитель из полиэтилена) является защитная броня,
представляющая из себя повив из, как правило, стальных прутьев. В качестве
завершающего наружного слоя, обеспечивающего герметичность, может
использоваться повив из полипропиленовых жгутов с пропиткой из битума
[12].
14
1.4 История развития конструкции силового подводного кабеля
На сегодняшний день подводный силовой кабель используется для
электроснабжения: островов; морских нефтяных и газовых платформ;
ветровых, волновых и приливных электростанций; океанологических
лабораторий [14].
Приведем краткое описание основных исторических событий укладки
подводных силовых кабелей и событий, определяющих направление
развития конструкции подводного силового кабеля к его современному
состоянию [11]:
1811
год:
первая
прокладка
подводного
силового
кабеля,
изолированного натуральным каучуком в Германии;
1924 год: использование свинцового щита в качестве водяного барьера;
1937 год: первый кабель с синтетической изоляцией из бутилкаучука;
1952 год: использование маслонаполненной изоляции в качестве
изоляции кабеля;
1954 год: первая прокладка кабеля постоянного тока высокого
напряжения между о. Готланд и г. Вестервик в Швеции длиной 98 км;
1962 год: первой использование этиленпропиленового каучука в
качестве изоляции кабеля;
1973 год: первой использование сшитого полиэтилена в качестве
изоляции;
1990-2000 года: замена маслонаполненной изоляции кабеля на
пластмассовую изоляцию.
Конструкция подводного силового кабеля различается в зависимости
от используемых материалов и непосредственных условий укладки кабеля.
Принято разделять силовые кабели на два типа: высокого напряжения с
переменным током, используемого для кабельных линий, не превышающих
80 км из-за ограничения по дальности; высокого напряжения с постоянным
током, используемого для кабельных линий большой дальности [11].
15
Приведем типичную конструкцию современного подводного силового
кабеля [15]:
Рис. 1.7 Конструкция современного силового подводного кабеля [15]
Центральной конструктивной составляющей современного подводного
силового кабеля являются токопроводящая жила (иногда несколько, как
правило, медные), предназначенная для осуществления энергоснабжения.
Следующим конструктивным слоем (через экран-сепаратор по жиле)
является изоляция из сшитого полиэтилена или бумаги, заполненной
жидкостью. Иногда в связку из токопроводящих жил устанавливается
оптоволоконная
линия
для
передачи
информации.
Защитным
конструктивным слоем также является защитная броня, представляющая
собой повив из, как правило, стальных прутьев. В качестве завершающего
наружного
слоя,
обеспечивающего
герметичность,
также
может
использоваться повив из полипропиленовых жгутов с пропиткой из битума.
Отметим, что между вышеописанными слоями повсеместно используются
гидрофобный заполнитель или полимерные ленты. Иногда в конструкцию
подводного
силового
кабеля
также
включается
свинцовый
щит,
предназначенный для защиты от воды [11].
В настоящей диссертации рассматривается влияние микроструктуры
подводного силового кабеля различного типа на форму и натяжение
провисающего участка кабеля при его укладке под водой.
16
1.5 Обзор современных работ по механике деформируемого кабеля
Статьи 2008 и 2010 годов Немова А.С., Войнова И.Б., Боровкова А.И. и
др. [16,17] посвящены расчетному определению жесткостных характеристик
кабелей
с
иерархической
структурой.
Обобщенные
эффективные
коэффициенты жесткости определяются исходя из решения серии задач
теории упругости о растяжении, изгибе и кручении кабелей со сложной
структурой с применением метода конечных элементов для канатов
сверхпроводящей катушки тороидального поля токамака ITER.
Статья 2015 года D. Zhang и M. Ostoja-Starzewski [18] посвящена
исследованию изгибной жесткости участка троса. Используется модель
контактного взаимодействия между прутьями с учетом трения Кулона. В
работе
представлены
распределения
модуля
вектора
перемещений,
контактного давления и кривой изгибающего момента силы для различных
коэффициентов трения и натяжения кабеля. Статья 2016 года Y. Yu и др. [19]
посвящена исследованию трехточечного изгиба участка троса. Также
используется модель контактного взаимодействия между прутьями с учетом
трения Кулона. В работе представлены распределения модуля вектора
перемещений и кривая изгибающей силы от прогиба троса для различных
конфигураций троса.
Статья 2016 K.H. Leong и др. [20] посвящена вопросам разрушения
кабеля для дистанционного динамика микрофона. В работе учитывается
сложная гетерогенная микроструктура кабеля. Также используется модель
контактного взаимодействия между прутьями с учетом трения Кулона.
Используется
нелинейная
деформирования.
модель
материалов
с
учетом
кривой
В работе представлены кривые усилий кабеля до
разрушения для задач о растяжении и трехточечном изгибе. Статья 2018 F.F.
Luz и др. [21] посвящена исследованию разрушения троса при растяжении.
Также используется модель контактного взаимодействия между прутьями с
учетом трения Кулона. В работе представлено критическое усилие троса при
разрыве.
17
Работа 2011 года M. Karahan и Ö Kalenderli [22] посвящена
исследованию температурного состояния и токовой нагрузки силового
кабеля с гетерогенной структурой, проложенных в почве. В работе
представлено влияние взаимного расположения нескольких кабелей на
температуру, силу тока и срок службы. Статья 2019 года Z. Xu и др. [23]
посвящена исследованию температурного состояния подводного силового
кабеля с гетерогенной структурой. В работе рассматривается влияние на
температуру кабеля следующих факторов: глубина прокладки кабеля от
поверхности морского дна, сила тока, свойства кабеля и окружающей среды,
температура окружающей среды. Статья 2014 года C. Holyk и др. [24] также
посвящена исследованию температурного состояния силового кабеля с
гетерогенной структурой в зависимости от токовой нагрузки для различных
геометрических конфигураций кабеля.
Статьи 2016 года S. Dubitsky и др. [25] и 2018 года J.C. del-Pino-López
[26] и др. посвящены исследованию вопросов электродинамики силового
кабеля
с
гетерогенной
структурой. В
первой
работе
представлено
распределение плотности Джоулева тепла для кабеля. Также рассматривается
зависимость температуры и магнитного поля от токовой нагрузки и
размещения кабельных линий. Во второй работе представлено распределение
электромагнитного потенциала, магнитного потока, плотности тока для
участка силового кабеля.
В настоящей диссертации определяются обобщенные коэффициенты
жесткости для подводных силовых кабелей различных типов со сложной,
многокомпонентной
межкомпонентного
гетерогенной
контактного
структурой
взаимодействия
[27].
Для
описания
используется
модель
клеевого соединения. Эффективные характеристики кабелей определяются
исходя из решения серии задач теории упругости о растяжении, изгибе и
кручении кабелей. В дальнейшем, эти характеристики используются для
описания
пространственного
поведения
комбинированного нагружения.
18
кабеля
под
действием
Глава 2. Современное состояние математических моделей об укладке
кабеля и смежных задачах о движении кабелей под водой
Вторая
глава
диссертации
посвящена
описанию
современного
состояния задачи об укладке кабеля под водой, известной в литературе как
задача Рауса-Аппеля [28-34]. Приводится история открытия и развития
математической модели задачи Рауса-Аппеля. Задача об укладке кабеля под
водой обладает косвенной взаимосвязью с такими соседними классами задач,
как: буксировка, швартование, поведение системы «судно – кабель –
подводный аппарат». Эта схожесть обусловлена использованием теории
механики нити и наличием гидродинамических сил сопротивления. В связи с
этим, также приводится описание особенностей некоторых математических
моделей движения нити под водой для этих соседних классов задач.
2.1 Аналитическое решение задачи Рауса-Аппеля
Прежде чем говорить об истории открытия и развития математической
модели укладки кабеля под водой и тем более – о ее расширении физически
обоснованными нелинейными характеристиками, необходимо, в первую
очередь познакомиться с задачей Рауса-Аппеля. Данная задача представлена
в том или ином виде в нескольких источниках [28-34] в немного
различающихся трактовках с использованием математических моделей,
отвечающих современным, на тот момент, данным. Для того чтобы ввести
читателя в курс дела, приведем наиболее лаконичную и современную
трактовку этой задачи, следуя за Меркиным [32]. Будем придерживаться
терминологии, используемой Меркиным и представим аналитическое
решение в немного расширенном и алгоритмизированном виде.
Пусть катушка (см. рис. 2.1) движется горизонтально с постоянной
скоростью V . При этом, с катушки сматывается кабель со скоростью схода,
совпадающей со скоростью движения катушки. Кабель укладывается на
неподвижную горизонтальную плоскость глубиной H в сопротивляющейся
среде. При равномерном и прямолинейном движении катушки, с течением
19
времени, провисающий участок кабеля приобретает форму динамического
равновесия, в рамках которого движение кабеля осуществляется вдоль линии
«кажущегося покоя» в подвижной системе координат, прямолинейно
перемещающейся в пространстве вместе с катушкой. Исследуемыми
характеристиками являются форма и натяжение провисающей части кабеля.
Рис. 2.1 Постановка задачи Рауса-Аппеля [34]
В качестве математической модели кабеля будем использовать модель
идеально гибкой нити. Гипотеза о равномерном движении нити вдоль линии
«кажущегося покоя» позволяет записать уравнения контурного движения
нити как уравнения статического равновесия:
*
dT
ds
P 0, f T ,
ds
dl
(2.1)
где T * – кажущееся натяжение нити, P – погонная равнодействующая сил,
действующих на нить после растяжения, s – дуговая координата после
растяжения, l – дуговая координата до растяжения, T – действительное
натяжение нити, f T – закон растяжимости нити.
Спроектируем
внешние
силы,
представленные
погонной
равнодействующей сил P на касательную и главную нормали нити.
Отметим, что в данной постановке задачи погонная равнодействующая сил P
состоит из гидродинамической силы сопротивления воды и веса нити в
воде q . Для случая плоского движения нити проекции уравнения (2.1) на
касательную и главную нормаль имеют вид:
20
dT *
q sin 1 cos ,
ds
T*
q cos n sin 2 ,
(2.2)
где радиус кривизны определяется следующими соотношениями:
ds
1 dx
1 dy
.
d cos d sin d
(2.3)
Подставляя выражения (2.3) в уравнения (2.2), нетрудно получить систему
линейных дифференциальных уравнений для определения кажущегося
натяжения нити T * , абсциссы x , ординаты y и дуговой координаты s :
dT * q sin 1 cos *
T ,
d
q cos n sin 2
dx
cos
T *,
2
d q cos n sin
dy
sin
T *,
d q cos n sin 2
ds
1
T *.
2
d q cos n sin
(2.4)
Граничные условия имеют вид:
x 0 0,
y 0 0,
(2.5)
y макс H ,
s 0 0.
В такой постановке задачи решение можно получить аналитически.
Преобразуя первое уравнение системы (2.4), легко получить уравнение
dT * q sin 1 cos
d ,
T*
q cos n sin 2
(2.6)
решением которого является
T * С1e F ,
q sin 1 cos
F
d ,
q cos n sin 2
0
(2.7)
где угол изменяется в пределах:
q q 2 4
n
0; arccos
2n
21
,
(2.8)
а верхний предел определяется особенной точкой системы (2.4). Таким
образом, зная глубину укладки и зависимость F можно определить
константу интегрирования С1 :
С1 макс
0
Пользуясь
найденной
H
F
e sin
d
q cos n sin 2
константой
(2.9)
.
интегрирования
и
оставшимися
уравнениями системы (2.4) можно записать уравнения, определяющие форму
провисающего участка нити:
e F cos
0 q cos n sin 2 d ,
макс
e F sin
y С1
d ,
2
q
cos
sin
n
0
макс
e F
s С1
d .
2
q
cos
sin
n
0
x С1
макс
(2.10)
Динамическое натяжение нити T определяется выражением
T T * V 2 ,
(2.11)
где – погонная масса нити. В заключение Меркин предлагает учесть
влияние растяжимости нити f T на вес нити в воде q численным методом
на ЭВМ.
2.2 История развития задачи об укладке кабеля
Задача определения формы и натяжения гибкой нити в потоке воды или
воздуха издавна привлекала внимание исследователей. Знание формы нити в
потоке необходимо не только для укладки кабеля, но и в минном деле при
постановке и тралении подводных мин, гидрографам при промере глубин
лотом, в текстильной промышленности, при расчете глубоководных якорей,
буксирных тросов, привязных аэростатов, линий высоковольтных передач и
в технике промышленного рыболовства.
В вышеописанной постановке эта задача, известная в литературе как
задача Рауса-Аппеля, привлекла внимание исследователей еще в конце
22
середины XIX века в связи с укладкой трансатлантического телеграфного
кабеля. Впервые данную задачу сформулировал Э. Раус в 1860 году для
нерастяжимой
однородной
весомой
нити
[28].
В
отсутствие
экспериментальных данных о сопротивлении кабелей, Э. Раус принял
допущение, что сопротивление воды перемещению участка кабеля действует
противоположно направлению перемещения, не зависит от угла атаки и
пропорционально первой степени скорости. В соответствие с этим
допущением на каждую единицу массы нити действуют две равные по
величине силы сопротивления, одна из которых появляется за счет
перемещения участка кабеля вдоль формы «кажущегося покоя», а вторая – за
счет горизонтального движения этой формы.
Полное аналитическое решение этой задачи для нерастяжимой
однородной весомой гибкой нити в таких предположениях П. Аппель
приводит в своем мемуаре в 1893 году после установления целого ряда
теорем [29].
Для случая плоского расположения весомой нерастяжимой гибкой нити
в однородном потоке наиболее законченное решение задачи получил А.Н.
Крылов, занимавшийся этой задачей в 1909 году в приложении к якорному
тросу подводной мины. А.Н. Крылов предположил, что приходящаяся на
единицу массы нити сила гидродинамического сопротивления действует в
плоскости потока перпендикулярно направлению нити и по величине
квадратично пропорциональна скорости с синусом угла атаки нити [35]. За
плоскость потока принимается плоскость, образованная направлением нити в
рассматриваемой точке с направлением потока.
В 1941 году А.П. Минаков расширил аналитическое решение П. Аппеля
учетом растяжимости нити [30].
Экспериментальные исследования сопротивления нитей, канатов и
проволок в потоке проводились с начала XX века, со времени появления
первых аэродинамических труб. Первые продувки тросов проводились на
сравнительно больших углах атаки [36-42]. Продувка канатов на малых углах
23
атаки требовала более тонкой методики эксперимента. Только в 1961-1962
годах
И.Р.
Матросовым
экспериментальные
и
результаты
С.Н.
Чубаровым
определения
были
опубликованы
коэффициентов
лобового
сопротивления и «подъемной силы» для каната при различных углах атаки
[43,44],
подтвердившие
модель
гидродинамического
сопротивления,
предложенную А.Н. Крыловым.
В 1966-1970 годах Н.И. Алексеев, следуя за А.П. Минаковым, получил
аналитическое решение для задачи Рауса-Аппеля для тяжелой однородной
гибкой нити с учетом математической модели силы лобового сопротивления
[31,45], предложенной А.Н. Крыловым и экспериментально подтвержденной
И.Р. Матросовым и С.Н. Чубаровым.
В
своей
монографии
1980
года
Д.Р.
Меркин
[32]
приводит
аналитическое решение задачи Рауса-Аппеля и конечный вид интегральных
уравнений, описывающих форму и натяжение провисающей части тяжелой
нерастяжимой
нити
в
декартовой
системе
координат,
связанной
с
движущимся судном.
2.3 Описание смежных задач о швартовании
В статьях 2015-2016 годов Z. Ren, R. Skjetne и др. [46-49] рассматривают
задачу о численном моделировании швартования судна для фиксации
месторасположения на поверхности воды.
В работе [46] в качестве численного метода исследования применяется
метод конечных элементов. В рамках внешних сил действующих на кабель
учитывается действие силы тяжести, силы Архимеда, силы нормального и
касательного гидродинамического сопротивления, связанные с обтеканием
кабеля подводным течением и влияние присоединенных масс жидкости.
Кроме этого, также стоит отметить, что в этой работе представлены
различные профиля скорости подводного течения, в том числе характерный
профиль подводного течения для Мексиканского залива и Норвежского моря
вблизи газового месторождения Ормен Ланге. В работе [47] в качестве
24
численного метода исследования также применяется метод конечных
элементов и также в рамках внешних сил, действующих на кабель,
учитывается действие силы тяжести, силы Архимеда, силы нормального и
касательного гидродинамического сопротивления, связанные с обтеканием
кабеля подводным течением и влияние присоединенных масс жидкости.
В работах [48,49] форма кабеля принимается совпадающей с цепной
линией и определяется исходя из данных измерений датчиков натяжения
линий швартования.
2.4 Описание смежных задач о буксировке и совместном движении системы
«судно – кабель – подводный аппарат»
В статье 2015 года J. Park и N. Kim [50] рассматривают задачу о
численном моделировании буксировки подводного робота кабелем в
сопротивляющейся среде. В качестве численного метода исследования
применяется формулировка полных узловых координат. В рамках внешних
сил, действующих на кабель, учитывается действие силы тяжести, силы
Архимеда,
силы
нормального
и
касательного
гидродинамического
сопротивления связанного с движением кабеля при буксировке, а также
влияние присоединенных масс жидкости.
В статье 2017 года M. Vu и др. [51] рассматривают совместное движение
системы «судно – подводный робот» в сопротивляющейся среде. Для
решения
системы
стрельбы.
В
математическая
дифференциальных
качестве
модель
математической
нерастяжимой
уравнений
модели
нити.
применяется
кабеля
Подводное
метод
используется
течение
не
учитывается. В рамках внешних сил, действующих на кабель, учитывается
действие силы тяжести, силы Архимеда, силы нормального и касательного
гидродинамического сопротивления, связанные с движением кабеля при
буксировке.
В статьях 2009, 2017-2018 годов В.В. Костенко и др. [52-54]
рассматривают движение подводного робота связанного кабелем с судном. В
25
качестве численного метода применяется программа ZONA для гибкой
нерастяжимой нити. В рамках внешних сил, действующих на кабель,
учитывается действие силы тяжести, силы Архимеда, силы нормального и
касательного гидродинамического сопротивления вызванная обтеканием в
однородном потоке воды.
2.5 Выводы, описание возможных направлений исследования
Проведен обзор научной литературы по истории открытия и развития
математических моделей задачи об укладке кабеля под водой, известной в
литературе как задача Рауса-Аппеля. Приведено современное состояние
особенностей математической модели об укладке кабеля и смежных задач,
связанных с исследуемой использованием теории механики нити и наличием
сил сопротивления среды.
Резюмируя проведенный обзор литературы, приведем описание
возможных направлений расширения математической модели задачи РаусаАппеля. Проведение работ при укладке кабеля в реальных морских условиях
нестационарное и сопровождается следующими физическими явлениями:
пространственным движением кабеля;
неровным профилем морского дна;
отклонением судна-кабелеукладчика от курса и скорости хода;
рассогласованным скоростным режимом укладки (случай, в котором
скорость схода кабеля не совпадает со скоростью движения суднакабелеукладчика);
морским волнением (бортовая и килевая качка);
подводным течением, изменчивым по глубине;
изгибом и кручением кабеля;
невертикальной ориентацией силы Архимеда;
действием присоединенных масс жидкости.
Некоторые из этих физических явлений учитываются в последующих главах
диссертации.
26
Глава 3. Математическое моделирование динамического равновесия
кабеля при его укладке под водой
Третья глава диссертации посвящена математическому моделированию
укладки кабеля как механической системы динамического равновесия в
состоянии, т.н. «кажущегося покоя». В настоящей главе рассматривается
влияние таких нелинейных эффектов, как:
растяжимость кабеля;
подводное течение с профилем скорости, изменчивым по глубине.
Наличие этих нелинейных особенностей при укладке кабеля при
использовании теории механики нити и формализации к виду задача РаусаАппеля вырождается в систему нелинейных дифференциальных уравнений.
Для решения задачи о динамическом равновесии применяется численный
итерационный многомерный метод Ньютона на базе конечно-разностной
схемы [55] с использованием разработанного комплекса программ на языке
программирования Matlab MathWorks [56,57].
3.1 Численное конечно-разностное решение задачи Рауса-Аппеля
Альтернативным
аналитическому
решению
задачи
Рауса-Аппеля,
представленному в п. 2.1, является численное решение системы линейных
дифференциальных уравнений (2.4) с граничными условиями вида (2.5).
Применим конечно-разностную схему к (2.4) и сформируем систему
линейных алгебраических уравнений вида:
Ti*1 Ti* q sin i 1 cos i *
Ti ,
i 1 i
q cos i n sin 2 i
xi 1 xi
cos i
Ti* ,
2
i 1 i q cos i n sin i
yi 1 yi
sin i
Ti* ,
2
i 1 i q cos i n sin i
si 1 si
1
Ti* ,
2
i 1 i q cos i n sin i
(3.1)
где i [1; N 1]. Отметим, что каждое из уравнений системы (3.1) содержит
одинаковое количество конечно-разностных N 1 уравнений, в то время как
27
количество неизвестных различается: N
неизвестных по кажущемуся
натяжению нити T * , N 2 неизвестных по ординате y , N 1 неизвестных по
абсциссе x и дуговой координате s .
Отличительной
особенностью
численного
решения
системы
алгебраических уравнений (3.1) с граничными условиями вида (2.5) является
необходимость первоочередного решения первого и третьего уравнения по
кажущемуся натяжению нити T * и ординаты y совместно. Воспользуемся
численным методом Гаусса:
MФ W,
(3.2)
где Ф – вектор-столбец неизвестных, M – основная матрица системы, W –
вектор-столбец правых частей. Применительно к рассматриваемой задаче, Ф ,
M и W будут иметь следующий вид:
TN* 1
Ф
,
y ( N 2)1
(3.3)
M1 ( N 1)( N ) M2 ( N 1)( N 2)
M 2( N 1)2( N 1)
,
M3
M4
(
N
1
)
(
N
)
(
N
1
)
(
N
2
)
*
T
F1
1
0
0
0
T1
0
0
0
*
T Fi
M1 ( N 1)( N ) 0
0
1
0
T
j
0
0
0
*
T FN 1
0
0
0
0
TN 1
M3 ( N 1)( N )
y F1
T1
0
0
0
0
M2 ( N 1)( N 2) 0,
0
0
0
0
0
0
0
0
0
y Fi
T j
0
0
0
0
0
y FN 1
TN 1
0
28
0
0
0 ,
0
1
(3.5)
(3.6)
0
(3.4)
0
0
0,
0
0
(3.7)
M4 ( N 1)( N 2 )
1 0
1 1
0
0 0
0 0
0 0
0 0
0
0 0 0 0
0 0 0
1 1 0 0 .
0 0
0 0 1 1
0 0 0 1
0
0
0
(3.8)
0 N 1
y1
0
M
,
0
y
N ( N 2)1
где
T*
Fi
– уравнения по кажущемуся натяжению T * ,
(3.9)
y
Fi
– уравнения по
ординате y , которые нетрудно сформировать по системе (3.1).
После определения кажущегося натяжения T * и ординаты y , абсцисса
x
и дуговая координата определяются в явном виде в соответствие с
конечно-разностными уравнениями:
Ti* i 1 i cos i
,
q cos i n sin 2 i
Ti* i 1 i
si 1 si
,
q cos i n sin 2 i
xi 1 xi
(3.10)
где i [1; N 1]. Заключительным этапом является переход от кажущегося
натяжения T * к действительному натяжению T :
Ti Ti* V 2 ,
где i [1; N ].
(3.11)
3.2 Учет растяжимости кабеля согласно закону Гука
Расширим постановку задачи Рауса-Аппеля, представленную в п. 2.1,
учетом влияния растяжимости нити согласно закону Гука на погонный вес
нити в воде q :
q
q0
q0
q0
,
*
T V 2
f T 1 T
1
E
E
29
(3.12)
где – площадь поперечного сечения нити, q0 – погонный вес нерастянутой
нити в воде, E – модуль упругости нити.
Перепишем систему (2.4) с учетом выражения (3.12) и сформируем
систему нелинейных дифференциальных уравнений, учитывающую влияние
растяжимости нити согласно закону Гука:
dT * q0E sin 1 cos E T * V 2 *
T ,
d
q0E cos n sin 2 E T * V 2
dx
cos E T * V 2
T *,
d q0E cos n sin 2 E T * V 2
dy
sin E T * V 2
T *,
2
*
2
d q0E cos n sin E T V
ds
E T * V 2
T *.
2
*
2
d q0E cos n sin E T V
(3.12)
Кроме этого, растяжимость нити также влияет на верхний предел
интегрирования макс , представленный выражением (2.8), т.к. особенная точка
системы (3.12) изменилась по сравнению с системой (2.4):
q E q 2 2 E 2 42 E T V 2 2
0
0
n
yH
0; arccos
.
2
2n E TyH V
(3.13)
Данная нелинейная особенность исключает возможность применения
аналитического подхода для определения кажущегося натяжения нити T * ,
ординаты y , абсциссы x и дуговой координаты s . В связи с чем будет
использован численный подход. Применим конечно-разностную схему к
(3.12) и сформируем систему нелинейных алгебраических уравнений вида:
Ti*1 Ti* q0E sin i 1 cos i E Ti* V 2 *
Ti ,
i 1 i
q0E cos i n sin 2 i E Ti* V 2
xi 1 xi
cos i E Ti* V 2
T *,
i 1 i q0E cos i n sin 2 i E Ti* V 2 i
yi 1 yi
sin i E Ti* V 2
Ti* ,
2
*
2
i 1 i q0E cos i n sin i E Ti V
si 1 si
E Ti* V 2
Ti* ,
2
*
2
i 1 i q0E cos i n sin i E Ti V
(3.14)
где i [1; N 1]. Как и ранее, первое и третье уравнения по кажущемуся
натяжению нити T * и ординаты y необходимо решать совместно.
30
Для
решения
системы
нелинейных
алгебраических
уравнений
необходимо использовать итерационные методы. Воспользуемся численным
многомерным методом Ньютона, который имеет вид:
J ( k )ΔΦ( k ) F ( k ) ,
Ф( k 1) Ф( k ) ΔФ( k ) ,
(3.15)
где J (k ) – матрица Якоби, ΔФ(k ) – поправка к приближению для системы, F (k ) –
невязка системы, Ф (k ) – приближение для системы. Отметим, что верхний
предел выражения (3.13) необходимо обновлять на каждой итерации
численного решения. В последующих записях индекс k опускается в целях
формирования
менее
громоздких
выражений.
Применительно
к
рассматриваемой задаче, Ф , F и J имеют следующий вид:
TN* 1
Ф
,
y ( N 2)1
T*
(3.16)
2
Fi A1Ti * A2Ti *1 A3Ti * A4Ti *Ti *1 ,
2
Fi B1yi B 2 yi 1 B3Ti* B 4Ti* B5 yiTi * B6 yi 1Ti * ,
A1 q0E cos i E V 2 n sin 2 i i 1 i q0E sin i E V 2 1 cos i ,
A2 q0E cos i E V 2 n sin 2 i ,
A3 n sin 2 i i 1 i 1 cos i ,
y
A4 n sin 2 i ,
B1 A2,
B 2 A2,
B3 i 1 i E V 2 sin i ,
B 4 i 1 i sin i ,
B5 A4,
B6 A4,
(3.17)
где i [1; N 1].
J1( N 1)( N ) J2 ( N 1)( N 2)
J 2( N 1)2( N 1)
,
J3 ( N 1)( N ) J4 ( N 1)( N 2)
31
(3.18)
F
T F1
T1
0
*
T Fi
0
T
j
0
0
*
J1 ( N 1)( N )
T*
1
T2
0
0
0
F F
0
T*
T*
i
i
0
T j
0
T j 1
0
0
0
0
0
0
0
0
FN 1
TN 1
0,
T Fi
J2( N 1)( N 2)
y j
0
T*
0
FN 1
TN
T*
, (3.19)
*
J3( N 1)( N )
y F1
T
1
0
y Fi
0
T j
0
0
y F1
y
y2
F2
y
2
0
y Fi
0
y j
0
0
0
0
F
J4( N 1)( N 2 )
где
T*
Fi
(3.20)
0
0
0
0
0
0
0
0
0
0
y Fi
T j
0
0
0
0
0
y FN 1
TN 1
0
0
0
0,
0
0
(3.21)
0
0
0
0
0
0
0
0
y Fi
y j
0
0
y Fi
y j 1
0
0
0
0
y
2
y3
0
0
0
0
y FN 2
y N 2
0
0
0
0
0
– уравнения по кажущемуся натяжению T * ,
0
,
(3.22)
F
y
N 2
y N 1
y FN 1
y N 1
y
Fi
– уравнения по
ординате y , которые нетрудно сформировать по системе (3.17).
Отметим, что важной составляющей численного многомерного метода
Ньютона является необходимость указания начального приближения, в
рамках которого будет использоваться решение для системы линейных
алгебраических уравнений (3.1) без учета растяжимости для идентичной
дискретизации (см. п. 3.1).
32
После определения кажущегося натяжения T * и ординаты y , абсцисса
x
и дуговая координата s определяются в явном виде в соответствие с
конечно-разностными уравнениями:
Ti* i 1 i cos i E Ti* V 2
xi 1 xi
,
q0E cos i n sin 2 i E Ti* V 2
Ti* i 1 i E Ti* V 2
si 1 si
,
q0E cos i n sin 2 i E Ti* V 2
(3.23)
где i [1; N 1]. Заключительным этапом является переход от кажущегося
натяжения T * к действительному натяжению T :
Ti Ti* V 2 ,
где i [1; N ].
(3.24)
3.3 Учет подводного течения с профилем скорости, изменчивым по глубине
Расширим постановку задачи Рауса-Аппеля с учетом растяжимости
нити согласно закону Гука, представленную в п. 3.2, дополнительным учетом
влияния гидродинамического сопротивления воды
~
, обусловленного
подводным течением. На рис. 3.1 представлена расширенная постановка
задачи Рауса-Аппеля.
Рис. 3.1 Постановка задачи Рауса-Аппеля с учетом подводного течения с профилем
скорости, изменчивым по глубине
Перепишем систему плоского движения нити в проекциях на главную и
касательную
нормали
(2.4)
с
учетом
~
сопротивления воды :
33
влияния
гидродинамического
dT *
~
q sin 1 cos cos q sin 1 cos C u y cos ,
ds
T*
~
q cos n sin 2 n sin 2 q cos n sin 2 n u 2 y sin 2 .
(3.25)
Здесь знак «-» соответствует встречному течению, а «+» – попутному.
Перепишем также систему (3.12) с учетом влияния гидродинамического
~
сопротивления воды :
dT * q0E sin 1 cos C u y cos E T * V 2 *
T ,
d
q0E cos n sin 2 nu 2 y sin 2 E T * V 2
dx
cos E T * V 2
T *,
2
2
2
*
2
d q0E cos n sin nu y sin E T V
dy
sin E T * V 2
T *,
2
2
2
*
2
d q0E cos n sin nu y sin E T V
ds
E T * V 2
T *.
d q0E cos n sin 2 nu 2 y sin 2 E T * V 2
Кроме этого, гидродинамическое сопротивление воды
(3.26)
~
также
оказывает влияние на верхний предел интегрирования макс , представленный
выражением (3.13), т.к. особенная точка системы (3.26) изменилась по
сравнению с системой (3.12):
q E q 2 2 E 2 4 u 2 H 2 E T V 2 2
0
0
n
n
yH
0; arccos
.
2
2
2
u
H
E
T
V
n
n
yH
(3.27)
Как и ранее будем использовать численный подход для определения
кажущегося натяжения нити T * , ординаты y , абсциссы x и дуговой
координаты s . Применим конечно-разностную схему к (3.26) и сформируем
систему нелинейных алгебраических уравнений вида:
Ti*1 Ti* q0E sin i 1 cos i C u yi cos i E Ti* V 2 *
Ti ,
i 1 i
q0E cos i n sin 2 i nu 2 yi sin 2 i E Ti* V 2
xi 1 xi
cos i E Ti* V 2
Ti* ,
2
2
2
*
2
i 1 i q0E cos i n sin i nu yi sin i E Ti V
yi 1 yi
sin i E Ti* V 2
T *,
i 1 i q0E cos i n sin 2 i nu 2 yi sin 2 i E Ti* V 2 i
si 1 si
E Ti* V 2
Ti* ,
2
2
2
*
2
i 1 i q0E cos i n sin i nu yi sin i E Ti V
34
(3.28)
где i [1; N 1]. Как и ранее, первое и третье уравнения по кажущемуся
натяжению нити T * и ординаты y необходимо решать совместно.
Для решения данной системы нелинейных алгебраических уравнений
воспользуемся численным многомерным методом Ньютона (3.15) как и
~
ранее. Применительно к рассматриваемой задаче, матрица Якоби J имеет
следующий вид:
~
~
J1( N 1)( N ) J2( N 1)( N 2)
~
J 2( N 1)2( N 1) ~
,
~
J3( N 1)( N ) J4( N 1)( N 2)
*
*
T F~1 T F~1
0
0
0
T2
T1
0*
0
0
* ~
*
~
T F~i
T Fi T Fi
0
0
0
T j
T j 1
T j
0
0
0
* ~
T FN 1
0
0
0
0
TN 1
(3.29)
~
J1( N 1)( N )
0
0*
T F~2
y 2
0
*
T F~i
~
J2 ( N 1)( N 2 )
0
y j
0
0
~
y F1
T1
0
y ~
Fi
~
J3 ( N 1)( N )
0
T j
0
0
0
0
0
0
0
0
0
0
0
0
0
F~
0
T*
i
0
0
0
y j
0
0
0
0
0
~
FN 1
y N 1
0
0
0
0
0
0
0
0
0
0
~
y Fi
T j
0
0
0
0
0~
y FN 1
TN 1
0
35
T*
0
0
0
~
FN 1
TN
T*
,
0
0
0 ,
0
0
, (3.30)
(3.31)
(3.32)
~
y F1
yy~2
F2
y
2
0
y F~i
0
y
j
0
0
0
0
F~
~
J4( N 1)( N 2 )
где
T*
~
Fi
0
0
0
0
0
0
0
0
~
y Fi
y j
0
0~
y Fi
y j 1
0
0
0
0
y
2
y3
0
0
0
0
~
y FN 2
y N 2
0
0
0
0
0
0
,
(3.33)
F~
– уравнения по кажущемуся натяжению T * ,
y
N 2
y N 1
~
y FN 1
y N 1
y
~
Fi
– уравнения по
ординате y , которые нетрудно сформировать по системе (3.28).
После определения кажущегося натяжения T * и ординаты y , абсцисса
x
и дуговая координата s определяются в явном виде в соответствие с
конечно-разностными уравнениями:
Ti* i1 i cos i E Ti* V 2
,
q0E cos i n sin 2 i nu 2 yi sin 2 i E Ti* V 2
Ti* i1 i E Ti* V 2
si1 si
,
q0E cos i n sin 2 i nu 2 yi sin 2 i E Ti* V 2
xi1 xi
(3.34)
где i [1; N 1]. Заключительным этапом является переход от кажущегося
натяжения T * к действительному натяжению T :
Ti Ti* V 2 ,
где i [1; N ].
(3.35)
3.4 Верификация и апробация математической модели динамического
равновесия кабеля при его укладке под водой
Для проведения математического моделирования необходимо ввести в
рассмотрение механические параметры кабеля (плотность кабеля каб ,
модуль
упругости
кабеля
E
и
диаметр
кабеля
d ),
свойства
сопротивляющейся среды (плотность воды и динамическая вязкость ) и
параметры укладки кабеля (согласованная скорость движения судна V и
глубина H ).
36
В табл. 3.1 приведены рассматриваемые механические параметры
кабелей:
№
Обозначение
кабеля
Табл. 3.1 Механические параметры кабелей [50,58-60]
Плотность
Модуль
Диаметр
кабеля каб ,
упругости кабеля
Литература
d
кабеля
,
м
E , Н/м2
кг/м3
1
▬▬▬▬
1300
7∙108
0.041
[50]
2
▬▬▬▬
3112.5
9∙109
0.047
[58]
3
▬▬▬▬
5500
2.55∙1010
0.1003
[59]
4
▬▬▬▬
7850
2.15∙1011
0.00599
[60]
Погонный вес нерастянутого кабеля в воде q0 определяется исходя из
соотношения:
q0 ~кабg каб воды g ,
(3.36)
где ~каб – плотность кабеля в воде, g – ускорение свободного падения, воды –
плотность воды.
Нормальная гидродинамическая сила сопротивления n и касательная
гидродинамическая
сила
сопротивления
определяются
исходя
из
соотношений [50,61]:
V2
d,
2 ,
C V ,
n Cn воды
(3.37)
где Cn – коэффициент нормального гидродинамического сопротивления, C –
коэффициент
касательного
гидродинамического
сопротивления.
Эти
коэффициенты будем определять исходя из соотношений:
Cn 1.1 4 Re 1/ 2 ,
C Nu ,
Vd
Re воды ,
(3.38)
Nu 0.55 Re1/ 2 0.084 Re 2 / 3 ,
где – динамическая вязкость воды, Nu – число Нуссельта,
[46-50,59], 0.0013 кг м с [61].
37
воды 1025 кг м3
Для обоснования разработанной математической модели динамического
равновесия кабеля при его укладке под водой с учетом растяжимости кабеля
согласно закону Гука и влияния подводного течения с профилем скорости
изменчивым по глубине, необходимо провести серию проверок.
Результаты математического моделирования приведены в безразмерном
виде для удобства их экстраполяции на любые глубины укладки кабеля.
Также отметим, что цвета графиков соответствуют обозначениям в табл. 3.1.
Основными
численными
настройками
используемого
численного
многомерного итерационного метода Ньютона являются: дискретизация
континуальной модели кабеля по переменной интегрирования – углу
приращения d ; количество итераций многомерного метода Ньютона k . На
рис. 3.2 приведена зависимость максимального натяжения от количества
разбиений конечно-разностной схемы, а на рис. 3.3 приведена зависимость
максимального натяжения от количества итераций многомерного метода
Ньютона.
Рис. 3.2 Зависимость максимального натяжения от количества разбиений конечноразностной схемы (на примере кабеля №3 при V = 3 узла)
38
Рис. 3.3 Зависимость максимального натяжения от количества итераций многомерного
метода Ньютона (на примере кабеля №1 при V = 3 узла)
Анализируя результаты, изображенные на рис. 3.2 и 3.3, необходимо
отметить, что наблюдается сходимость используемого численного метода,
как по количеству разбиений конечно-разностной схемы, так и по количеству
итераций многомерного метода Ньютона. Этот факт говорит о том, что
существует
определенная
дискретизация
определенное количество итераций
k,
угла
приращения
d
и
демонстрирующие сошедшееся
численное решение.
Очень важно отметить, что кривые сходимости, представленные на рис.
3.2
и
3.3
зависят
от
механических
параметров
кабелей,
свойств
сопротивляющейся среды и параметров укладки кабеля. Поэтому при
проведении
инженерных
расчетов
по
оценке
формы
и
натяжения
провисающей части кабеля с использованием разработанного комплекса
программ, необходимо, в первую очередь, проводить анализ сходимости в
соответствие с предлагаемым шаблоном (см. прил. А, часть 6, №1 и №2).
39
Также, важно понимать, что сходимость по количеству итераций
многомерного метода Ньютона определяется, в том числе и начальным
приближением. В случае если эффекты нелинейной модели оказывают
существенное влияние на решение по сравнению с линейной моделью, то
можно наблюдать отсутствие сходимости. В этом случае рекомендуется
ступенчатое добавление нелинейных эффектов с фиксированием результатов
для возможности их использования в качестве начального приближения.
Другой
проверкой
достоверности
разработанной
математической
модели является ее совпадение с линейной задачей при пренебрежительно
малых, но ненулевых нелинейных характеристиках. На рис. 3.4 и 3.5
представлено сравнение методов решения при механических параметрах
кабеля, обеспечивающих нерастяжимость (на примере модифицированного
кабеля с модулем упругости E = 2.55∙1014 Н/м2): аналитическое решение для
нерастяжимого кабеля (см. п. 2.1); численное решение для нерастяжимого
кабеля (см. п. 3.1); численное решение для растяжимого кабеля (см. п. 3.3).
Рис. 3.4 Сравнение методов решения для малых нелинейных характеристик на примере
формы провисающего участка кабеля (при V = 3 узла))
40
Рис. 3.5 Сравнение методов решения для малых нелинейных характеристик на примере
натяжения провисающего участка кабеля (при V = 3 узла)
Сравнивая, представленные на рис. 3.4-3.5, результаты моделирования
при механических параметрах кабеля обеспечивающего нерастяжимость,
можно наблюдать отличное совпадение как для формы провисающего
участка кабеля, так и для его натяжения. Вместе с этим, любопытна
различная величина совпадения результатов (менее 1% для формы кабеля и
менее 3% для натяжения кабеля). Это различие объясняется отсутствием
граничных условий для натяжения кабеля (см. систему (2.5)).
Удостоверившись
в
сходимости
используемого
численного
итерационного многомерного метода Ньютона, а также в вырождении
нелинейного
решения
к
линейному
при
малых
характеристиках
нелинейности, можно приступать к апробации разработанного комплекса
программ для оценки формы и натяжения провисающей части кабеля в
рамках принятых допущений математической модели. На рис. 3.6-3.7
представлена форма и натяжение провисающего участка кабеля для
различных
механических
характеристик
движения судна.
41
и
согласованных
скоростей
Рис. 3.6 Форма провисающего участка кабеля для различных механических характеристик
и согласованных скоростей движения судна
Рис. 3.7 Натяжение провисающего участка кабеля для различных механических
характеристик и согласованных скоростей движения судна
42
Анализирую результаты, представленные на рис. 3.6-3.7, нетрудно
заметить однозначное влияние согласованной скорости укладки кабеля на
форму
и
натяжение
провисающей
части
кабеля.
Нормальное
гидродинамическое сопротивление воды прямо пропорционально скорости
движения судна во второй степени, что при увеличении последнего приводит
к
существенному
увеличению
нормальной
гидродинамической
силы
сопротивления среды, что в свою очередь приводят к более протяженной
форме кабеля и большему натяжению как следствие.
Анализ влияния механических параметров кабеля на исследуемые
характеристики более затруднителен, в виду их участия в различной
комбинации в системе дифференциальных уравнений (3.28). Приобретаемая
форма кабеля преимущественно определяется весом кабеля в воде и
гидродинамической силой сопротивления. В случае увеличения веса кабеля
протяженность
будет
уменьшаться,
а
в
случае
увеличения
гидродинамической составляющей – наоборот.
3.5 Оценка формы и натяжения кабеля при его укладке вблизи газового
месторождения Ормен Ланге
В
рамках
разработанного
комплекса
программ
для
описания
зависимости скорости профиля подводного течения от глубины u y
используется полиномиальная зависимость вида:
~
~
~
~
u y Ay3 By 2 Cy D,
(3.39)
~ ~ ~
где регулировка профиля течения осуществляется коэффициентами A , B , C и
~
D.
В качестве практической иллюстрации, рассмотрим укладку кабеля в
Норвежском море вблизи газового месторождения Ормен Ланге [62,63]. В
соответствие с опубликованными данными по профилю подводного течения
в данной области [46], рассмотрим граничные условия следующего вида:
43
u y 0 0,
u y 0 0,
(3.40)
u yH u Н ,
u y H 0.
Формула (3.39) с учетом (3.40) приобретает конечную функциональную
зависимость:
u y
где
в качестве скорости
2uН 3 3uН 2
y 2 y,
H3
H
течения
на поверхности
(3.41)
uН
рассмотрим
экстремальную величину скорости (равную 0.24 м/с) в соответствие с
данными о распределении средней геофизической скорости вблизи газового
месторождения Ормен Ланге [62], представленной на рис. 3.8.
Рис. 3.8 Поле средней геофизической скорости течения вблизи газового месторождения
Ормен Ланге [62]
На рис. 3.9 представлена форма провисающего участка кабеля при
укладке кабеля в Норвежском море для случая встречного течения в
сопоставлении с решением, не учитывающим подводное течение.
44
Рис. 3.9 Форма провисающего участка кабеля при укладке кабеля в Норвежском море (на
примере кабеля №4 при V = 3 узла)
Анализируя постановку задачи, представленную на рис. 3.1, в
совместности с системой дифференциальных уравнений (3.25) и результаты
моделирования, представленные на рис. 3.9, нетрудно заметить, что учет
встречного течения формирует более протяженный провисающий участок
кабеля, что сопровождается большим натяжением кабеля.
3.6 Выводы
Численно
решена
задача
Рауса-Аппеля
для
системы
линейных
дифференциальных уравнений с применением конечно-разностной схемы и
численного
метода
Гаусса.
Разработана
математическая
модель
динамического равновесия кабеля при его укладке под водой с учетом
растяжимости кабеля согласно закону Гука и влияния профиля подводного
течения, изменчивого по глубине. Разработан комплекс программ для
инженерной оценки формы и натяжения провисающего участка кабеля на
языке программирования Matlab MathWorks.
45
Достоверность разработанной математической модели и комплекса
программ подтверждается соблюдением сходимости конечно-разностной
схемы по переменной интегрирования – углу приращения d , а также
соблюдением сходимости многомерного итерационного метода Ньютона по
количеству итераций k . Проведена проверка вырождения математической
нелинейной модели к известному линейному решению при малых
характеристиках нелинейности.
Проведена оценка формы и натяжения провисающего участка кабеля
для различных механических параметров кабелей, представленных в
литературе, а также для различных согласованных скоростей движения
судна. Рассмотрена укладка кабеля вблизи газового месторождения Ормен
Ланге в Норвежском море.
Использованием разработанного комплекса программ допускает выбор
пользователем учитываемых нелинейных эффектов, обладающих различной
вычислительной производительностью. В частности, время моделирования
для линейной задачи составляет около 1 минуты, а для нелинейной модели,
расширенной учетом растяжимости кабеля и влияния подводного течения
составляет около 2 часов. Такая сравнительно невысокая вычислительная
сложность допускает применение разработанного комплекса программ для
оперативной оценки укладки кабеля.
Вместе с этим, не все актуальные физические явления, неизбежно
сопровождающие укладку кабеля в реальных морских условиях, которые
представлены в п. 2.5, учтены. Нерассмотренные эффекты, в особенности
решение нестационарной пространственной задачи, формируют существенно
более
сложную
для
реализации
задачу
в
рамках
самостоятельно
разработанного комплекса программ на базе конечно-разностной схемы. Для
учета этих нелинейных особенностей укладки кабеля целесообразно
применить прямое конечно-элементное моделирование, чему посвящена
следующая глава диссертации.
46
Глава 4. Математическое моделирование движения кабеля при его
укладке под водой
Четвертая
глава
диссертации
посвящена
математическому
моделированию укладки кабеля как механической системы динамического
движения. В настоящей главе рассматривается влияние таких нелинейных
эффектов, как:
изгиб и кручение кабеля;
рассогласованное
динамическое
движение,
обусловленное
сматыванием кабеля с барабана лебедки.
Для решения задачи о динамическом движении кабеля под водой
применяется
прямое
конечно-элементное
моделирование
[64-69]
в
программной системе ABAQUS [70] с использованием разработанной
интегрированной программы на языке программирования Fortran.
4.1 Описание математической модели движения кабеля при его укладке под
водой
Для возможности описания нелинейных эффектов, связанных с учетом
изгиба и нестационарного движения кабеля при его укладке под водой,
необходимо расширить уравнение статического равновесия (2.1) учетом
действия сил инерции. Также отметим, что учет изгиба выходит за рамки
математической модели нити.
Для
учета
вышеописанных
эффектов,
рассмотрим
систему
дифференциальных уравнений для изгибных колебаний стержня БернуллиЭйлера при наличии осевой силы и продольных колебаний стержня
относительно прогиба w и удлинения u [71]:
EI
4 ws, t
ws, t
2 ws, t
T
s
,
t
qw s, t ,
s 4
s
s
t 2
2u s, t
2u s, t
E
qu s, t .
s 2
t 2
47
(4.1)
Для
реализации
прямого
конечно-элементного
моделирования
процессов укладки кабеля необходимо разработать конечно-элементную
модель. На рис. 4.1 представлены фрагменты конечно-элементной модели
укладки кабеля под водой.
б) область схода кабеля с
барабана лебедки
а) общий вид
в) радиус барабана лебедки
Рис. 4.1 Фрагменты конечно-элементной модели укладки кабеля
Приведем описание ключевых составляющих разработанной конечноэлементной модели, представленной на рис. 4.1:
Горизонтальная абсолютно жесткая поверхность, иллюстрирующая
поверхность дна;
Область ограниченная поверхностью дна и поверхностью воды,
иллюстрирующая сопротивляющуюся среду жидкости. Только в этой
области к конечным элементам кабеля прикладывается нестационарная
пространственная
сила
гидродинамического
48
сопротивления
в
соответствие с разработанной интегрированной программой для
приложения пользовательской нагрузки согласно текущим данным о
расположении, скорости движения и ориентации элемента кабеля
(более подробное описание см. в п. 4.2 и Прил. Б);
Математическая модель кабеля, а также его исходное положение в
воде. Кроме
этого,
планируемая
для
в
рамках
численного
вытравливания
длина
моделирования
кабеля
вся
расположена
вертикально над поверхностью воды и опускается в воду в
соответствие со скоростным законом схода кабеля;
Вспомогательная
вертикальная
поверхность,
расположенная
над
поверхностью воды, иллюстрирующая судно-кабелеукладчик, которая
перемещается горизонтально в пространстве в соответствие со
скоростным законом движения судна-кабелеукладчика;
Контактные взаимодействия между парами объектов: «кабель –
вспомогательная
вертикальная
поверхность»
и
«кабель
–
горизонтальная поверхность дна».
Процесс решения задачи состоит из двух ключевых этапов, граничные
условия для которых изображены на рис. 4.2:
Первый, предварительный, этап решения заключается в определении
деформированного состояния кабеля до начала движения суднакабелеукладчика и сматывания кабеля;
Второй этап решения моделирует укладку кабеля под водой. Кабель,
расположенный над поверхностью воды опускается в воду со
скоростью схода кабеля, а вспомогательная вертикальная поверхность
перемещается
горизонтально
со
кабелеукладчика.
49
скоростью
движения
судна-
а) первый этап
б) второй этап
Рис. 4.2 Граничные условия конечно-элементной модели укладки кабеля
4.2 Описание интегрированной программы по определению нестационарного
пространственного нагружения кабеля
Разработка
интегрированной
программы
по
определению
нестационарного пространственного нагружения кабеля и включение ее в
математическую модель движения кабеля при его укладке под водой
обоснована отсутствием стандартной силы в программной системе конечноэлементного
анализа
ABAQUS
(или
его
аналогов),
зависящей
от
расположения, скорости и ориентации элементов кабеля.
На рис. 4.3 представлена блок-схема разработанной математической
модели и участие в ней интегрированной программы по определению
пользовательской
нагрузки.
Отметим,
что
сила
гидродинамического
сопротивления прикладывается только на втором этапе решения (см. рис.
4.3).
Причем,
для
обеспечения
сходимости
численного
решения,
пользовательская нагрузка прикладывается постепенно (см. рис. 4.4).
50
Рис. 4.3 Блок-схема математической модели движения кабеля при его укладке под водой
Рис. 4.4 Закон постепенного приложения нагрузки на втором этапе
Отметим,
что
закон
постепенного
приложения
нагрузки,
представленный на рис. 4.4 не является универсальным для всех
механических параметров кабелей, свойств сопротивляющейся среды и
параметров укладки кабеля. В связи с чем, в некоторых частных случаях
может потребоваться его модификация.
51
4.3 Верификация конечно-элементной модели на примере задачи о натяжении
якорной цепи
Для обоснования разработанной конечно-элементной модели кабеля и
определения дискретизации конечно-элементной сетки, обеспечивающей
сошедшееся решение, проведем верификацию на примере известной задачи о
натяжении якорной цепи [32], представленной на рис. 4.5. Для верификации
настроек контактного взаимодействия рассмотрим случай, когда часть
якорной цепи лежит на дне и цепная линия начинается в точке O.
Рис. 4.5 Постановка задачи о натяжении якорной цепи
В виду широкой известности этой задачи в литературе, не будем
описывать аналитическое решение, а приведем только финальное уравнение
цепной линии в общем виде и применительно к вышеописанной задаче
(подробное аналитическое решение и алгоритм расчета приведен в [32]):
y ach
ch
x C1
C2 ,
a
(4.2)
l1
h
1 .
a
a
(4.3)
На рис. 4.6 приведено сопоставление аналитического и численного
подходов на примере формы провисающего участка кабеля, а на рис. 4.7
представлена
зависимость
максимального
натяжения
от
количества
конечных элементов, приходящихся на провисающий участок кабеля.
52
Рис. 4.6 Сопоставление аналитического и численного подходов на примере формы
провисающего участка кабеля
Рис. 4.7 Зависимость максимального натяжения от количества конечных элементов,
приходящихся на провисающий участок кабеля
53
4.4 Учет микроструктуры кабеля
Данный параграф посвящен математическому моделированию укладки
кабеля с учетом нелинейного эффекта связанного с изгибом и кручением
кабеля. Для учета этого эффекта рассматривается серия связанных друг с
другом задач:
Подготовка
полномасштабной
пространственной
геометрической
модели кабеля;
Разработка конечно-элементных моделей нагружения кабеля для задач
теории упругости о растяжении, изгибе и кручении кабеля;
Определение обобщенных коэффициентов жесткости математической
модели кабеля, эквивалентной на макроуровне многокомпонентной
гетерогенной конструкции кабеля;
Определение влияния микроструктуры кабеля на характеристики
укладки кабеля.
В качестве многокомпонентной гетерогенной конструкции кабеля
рассмотрим и сопоставим друг с другом 3 различные типа кабеля,
используемые при различных условиях укладки кабеля под водой:
трехжильный подводный силовой кабель с одиночным бронированием;
одножильный подводный силовой кабель с одиночным бронированием;
одножильный подводный силовой кабель с двойным бронированием.
4.4.1 Описание постановки задачи для определения эффективных коэффициентов
жесткости кабеля
Прямой учет микроструктуры кабеля при его укладке под водой не
представляется возможным, т.к. потребовал бы чрезмерно больших
вычислительных
деформированного
и
временных
состояния
ресурсов
для
описания
многокомпонентной
напряженногетерогенной
конструкции кабеля с целью определения его влияния на укладку кабеля под
водой. В связи с этим, будет проведена аналогия между кабелем и балкой с
обобщенными коэффициентами жесткости.
54
Запишем уравнения статического равновесия упругого каната при его
растяжении и кручении [72,73]:
A С T ,
C B M ,
где T
и
(4.4)
– приложенные к кабелю сила и крутящий момент
M
соответственно; и – продольные и угловые деформации кабеля; A , B и С
– обобщенные коэффициенты жесткости кабеля.
Граничные
условия
задач
теории
упругости
для
определения
эффективных коэффициентов жесткости кабеля представлены на рис. 4.8:
а) растяжение
б) изгиб
в) кручение
Рис. 4.8 Граничные условия задач теории упругости
В результате решения трех вышеописанных задач о растяжении, изгибе
и кручении кабеля можно определить эффективный модуль упругости кабеля
E , эффективный модуль сдвига G эффективный момент инерции кабеля на
изгиб I и эффективный момент инерции кабеля на кручение J .
55
Приведем ниже основные уравнения для вычисления эффективных
коэффициентов жесткости. В частности, эффективный модуль упругости
кабеля E определяется исходя из результатов решения первой задачи теории
упругости о растяжении (см. рис. 4.8а):
E
A 1 FZ
L,
S S UZ
(4.5)
где S – площадь поперечного сечения кабеля, FZ – реакция закрепленного
торца кабеля под действием приложенного перемещения U Z свободного
торца кабеля, а L – рассматриваемый участок кабеля.
Если пренебречь изменением площади поперечного сечения при
растяжении, то эффективный модуль сдвига G определяется исходя из
выражения:
G
E
.
2
(4.6)
Эффективный момент инерции кабеля на изгиб I определяется исходя
из аналогии между решением второй задачи теории упругости об изгибе (см.
рис. 4.8б) и уравнением изогнутой оси:
M Y L2
I
,
2 EU X
(4.7)
где M Y – изгибающий момент силы, а U X – прогиб свободного торца кабеля.
Эффективный момент инерции кабеля на кручение J определяется из
результатов решения третьей задачи теории упругости о кручении (см. рис.
4.8в):
J
B 1 MZ
L,
G G Z
(4.8)
где M Z – реакция закрепленного торца кабеля под действием приложенного
поворота Z свободного торца кабеля.
4.4.2 Описание конструкции силовых подводных кабелей
В научно-технической литературе и каталогах компаний-изготовителей
представлен широкий спектр кабелей различных типов, структуры и
56
габаритов. При выборе конкретной конструкции необходимо учитывать
такие
многочисленные
факторы
как:
целевое
назначение
кабеля;
механические характеристики защищенности; условия укладки кабеля;
стоимость материалов, изготовления и доставки кабеля [74].
Рассмотри распространенные и хорошо представленные в литературе
различные типы кабелей: трехжильный силовой подводный кабель с
одиночным бронированием; одножильный силовой подводный кабель с
одиночным бронированием; одножильный силовой подводный кабель с
двойным бронированием.
На
рис.
4.9а
представлены
фрагменты
полномасштабной
геометрической модели для трехжильного силового подводного кабеля с
одиночным бронированием 2XS2YRAA [75,76]. Этот кабель с диаметром 65
мм часто применяется для электроснабжения морских платформ, островов,
ветровых станций, а также прокладывается через реки и озера.
На
рис.
4.9б
представлены
фрагменты
полномасштабной
геометрической модели для одножильного силового подводного кабеля с
одиночным бронированием ZS-YJQ41 [23]. Этот кабель с диаметром 114 мм
прокладывается в Восточно-Китайском море и округе Пинтань.
На
рис.
4.9в
представлены
фрагменты
полномасштабной
геометрической модели для одножильного силового подводного кабеля с
двойным бронированием GASLMLTV [77,78]. Этот кабель с диаметром 67
мм
часто
применяется
в
условиях
повышенного
риска
вредного
вмешательства вследствие оживленного судоходства.
Основными
конструктивными
составляющими
кабелей,
представленных на рис. 4.9, являются: токопроводящая медная жила;
полиэтиленовая изоляция; экранирующие слои; армирующие повивы из
стальной проволоки; полипропиленовые жгуты.
57
а) трехжильный кабель с одиночным бронированием
б) одножильный кабель с одиночным бронированием
в) одножильный кабель с двойным бронированием
Рис. 4.9 Фрагменты полномасштабной пространственной геометрической модели
подводного силового кабеля
58
Отметим, что шаг свивки для разработанных полномасштабных
геометрических моделей кабелей, представленных на рис. 4.9, составляет 8
диаметров внешней окружности бронирования в соответствие с [79].
4.4.3 Описание конечно-элементных моделей силовых подводных кабелей
Для
определения
эффективных
коэффициентов
жесткости
рассматриваемых кабелей необходимо разработать конечно-элементные
модели рассматриваемых кабелей (см. рис. 4.9) для трех вышеописанных
задач теории упругости (см. 4.8). На рис. 4.10 представлены фрагменты
конечно-элементных моделей рассматриваемых кабелей.
Отметим, что решение задач теории упругости о растяжении, изгибе и
кручении проводилось в допущении квазистационарного нагружения,
вследствие чего применялась неявная схема интегрирования ABAQUS
Implicit.
Фрагменты конечно-элементных моделей рассматриваемых участков
кабелей, представленные на рис. 4.10, содержат до 0.5 млн. конечных
элементов. Средний размер конечно элемента в плоскости сечения кабеля
составляет 1 мм. Количество конечных элементов в осевом направлении,
приходящихся на шаг свивки кабеля, составляет 80 конечных элементов.
Отметим,
что
эффективных
проводились
коэффициентов
многочисленные
жесткости
от
проверки
сходимости
дискретизации
конечно-
элементной модели кабеля. В виду подобных зависимостей, в целях
исключения загромождения представим лишь сходимость метода конечных
элементов в осевом направлении на примере эффективного модуля
упругости E (см. рис. 4.11).
В качестве математической модели контактного взаимодействия для
многочисленных контактных взаимодействий была выбрана модель, не
допускающая разрыва границы разделов контактных поверхностей (модель
клеевого соединения, т.н. tie constraint).
59
а) трехжильный кабель с одиночным бронированием
б) одножильный кабель с одиночным бронированием
в) одножильный кабель с двойным бронированием
Рис. 4.10 Фрагменты конечно-элементной модели подводного силового кабеля
60
В рамках математической модели материала используется упругая
модель согласно закону Гука. Обозначение и свойства материалов для
конечно-элементных моделей кабелей приведены в табл. 4.1.
Табл. 4.1 Описание свойств материалов
№
1
2
3
Цвет
Материал
Медь
Полиэтилен
высокой
плотности
Полиэтилен
низкой
плотности
Модуль
Коэффициент Плотность,
упругости,
Литература
Пуассона
кг/м³
Па
125·109
0.34
8 920
[80]
1.1·109
0.43
950
[81]
0.17·109
0.45
920
[81]
4
Свинец
16·109
0.42
11 340
[82]
5
Битум
0.2·109
0.44
1 050
[83-85]
6
Оцинкованная
сталь
210·109
0.291
7 860
[82,86]
7
Полипропилен
1.45·109
0.41
910
[87]
Рис. 4.11 Зависимость эффективного модуля упругости от количества конечных элементов
в осевом направлении
61
4.4.4 Определение эффективных коэффициентов жесткости кабеля
Удостоверившись в сходимости метода конечных элементов, можно
приступить к определению эффективных коэффициентов жесткости кабеля.
На базе конечно-элементных моделей, представленных на рис. 4.10, были
численно решены три задачи теории упругости о растяжении, изгибе и
кручении (см. рис. 4.8). На рис. 4.12 приведено сопоставление поля осевого
перемещениях U Z для рассматриваемых типов кабелей на примере задачи
теории упругости о растяжении. На рис. 4.13 приведено сопоставление поля
изгибного перемещениях U X для рассматриваемых типов кабелей на примере
задачи теории упругости об изгибе. На рис. 4.14 приведено сопоставление
поля модуля вектора перемещений U для рассматриваемых типов кабелей на
примере задачи теории упругости о кручении.
Рис. 4.12 Поле осевого перемещения U Z для задачи о растяжении: а) трехжильный кабель
с одиночным бронированием; б) одножильный кабель с одиночным бронированием; в)
одножильный кабель с двойным бронированием
62
Рис. 4.13 Поле изгибного перемещения U X для задачи об изгибе: а) трехжильный кабель с
одиночным бронированием; б) одножильный кабель с одиночным бронированием; в)
одножильный кабель с двойным бронированием
Рис. 4.14 Поле модуля вектора перемещений U для задачи о кручении: а) трехжильный
кабель с одиночным бронированием; б) одножильный кабель с одиночным
бронированием; в) одножильный кабель с двойным бронированием
63
Отметим, что приведенные результаты (см. рис. 4.12-4.14) согласуются
с граничными условиями, приведенными на рис. 4.8. В частности, для всех
задач теории упругости конец кабеля, с граничными условиями заделки,
остается неподвижным. Нейтральная линия для задачи теории упругости о
кручении остается неподвижной. Отметим также, что условия контактного
взаимодействия на границах разделов соблюдаются.
В
соответствие
с
выражениями
(4.6)-(4.9)
и
результатами
моделирования задач теории упругости, представленными на рис. 4.12-4.14,
можно сформировать таблицу с обобщенными коэффициентами жесткости:
Табл. 4.2 Обобщенные коэффициенты жесткости кабелей
Обобщенный
коэффициент
жесткости
Модуль упругости
E , Па
Модуль сдвига G ,
Па
Жесткость на изгиб
EI , Н∙ м2
Момент инерции на
изгиб I , м4
Жесткость на
кручение B , Н∙ м2
Момент инерции на
кручение J , м4
Трехжильный
кабель с одиночным
бронированием
Одножильный
кабель с одиночным
бронированием
Одножильный
кабель с двойным
бронированием
8.9∙109
9.9∙109
32∙109
4.5∙109
4.9∙109
16∙109
7.6∙103
99∙103
20∙103
8.5∙10-7
101∙10-7
6.4∙10-7
8∙103
84∙103
12∙103
17.9∙10-7
171∙10-7
7.8∙10-7
Анализируя результаты, приведенные в табл. 4.2, можно увидеть, что
трехжильный силовой подводный кабель с одиночным бронированием
2XS2YRAA
обладает
наименьшими
обобщенными
характеристиками:
модуль упругости E , жесткость на изгиб EI , жесткость на кручение B . Такое
поведение с точки зрения механики, обусловлено сравнительно меньшим
диаметром и наличием только одного повива бронирования. А вот
сопоставление
эффективных
коэффициентов
жесткости
для
силовых
подводных кабелей ZS-YJQ41 и GASLMLTV не так очевидно, т.к. первый
большего диаметра, а второй – с двумя повивами бронирования.
64
Приведем объяснение взаимного расположения величин обобщенных
коэффициентов жесткости для одножильного силового подводного кабеля с
одиночным бронированием ZS-YJQ41 и одножильного силового подводного
кабеля GASLMLTV с двойным бронированием с точки зрения механики.
Кабель GASLMLTV обладает наибольшим эффективным модулем упругости
E . Это объясняется большей концентрацией гомогенных материалов с
высокими показателями жесткости.
Т.е. при растяжении бронирование
играет более важную роль по сравнению с габаритными размерами. Кабель
ZS-YJQ41
обладает
следующими
наибольшими
эффективными
характеристиками: жесткостью на изгиб EI и жесткостью на кручение B .
Это объясняется большей удаленностью материалов от нейтральной линии
кабеля. Т.е. при изгибе и кручении габаритные размеры играют более
важную роль по сравнению с бронированием.
4.4.5 Определение влияния микроструктуры кабеля на характеристики укладки
кабеля
Прежде чем проводить математическое моделирование движение
укладки кабеля как механической системы динамического движения с
учетом нелинейных эффектов связанных с изгибом и кручением кабеля,
выясним влияние микроструктуры кабеля на примере задачи о статическом
провисании кабеля, закрепленного по концам, под водой и сопоставим с
уравнением
цепной
линии
(см.
формулу
(4.3)).
По
результатам,
представленным на рис. 4.15, можно сказать, что только для кабеля ZSYJQ41 значение эффективной жесткости на изгиб EI таково, что для
описания формы этого кабеля в воде при укладке недостаточно использовать
теорию механики нити.
На рис. 4.16 представлены результаты движения кабелей различных
типов для постановки задачи, описанной в п. 4.1. Результаты движения
представлены в виде распределения натяжения кабеля на установившейся
форме кабеля, приобретаемой кабелем с течением времени.
65
Рис. 4.15 Сравнение формы кабеля для кабелей различных типов на примере задачи о
стационарном провисании кабеля в воде
а) трехжильный
б) одножильный
в) одножильный
кабель с
кабель с
кабель с
г) легенда
одиночным
одиночным
двойным
бронированием
бронированием
бронированием
Рис. 4.16 Распределения натяжения кабеля на установившейся форме кабеля
66
4.5 Учет рассогласованной укладки кабеля, вызванной сматыванием кабеля с
барабана лебедки
Данный параграф посвящен математическому моделированию укладки
кабеля с учетом нелинейного эффекта рассогласованного динамического
движения, вызванного сматыванием кабеля с барабана лебедки. Для учета
этого эффекта рассматривается серия связанных друг с другом задач:
Определение статического распределения реакций схемы крепления
барабана лебедки;
Определение момента силы, необходимого для запуска лебедки;
Определение
момента
силы
и
потребляемой
мощности
электродвигателя при вращении барабана лебедки согласно заданному
скоростному режиму работы;
Определение влияния рассогласованного динамического движения на
характеристики укладки кабеля.
В качестве лебедки, обеспечивающей сматывание кабеля, рассмотрим
механизм, представляющий из себя барабан с намотанным на него кабелем,
который опирается на вращающиеся опорные ролики (см. рис. 4.17) [88-90].
Барабан приводится в движение через ведущий вал с использованием
зубчато-цевочной передачи [91-95].
Рис. 4.17 Схема устройства лебедки
67
4.5.1 Определение статического распределения реакций схемы крепления барабана
лебедки
Рассмотрим задачу об определении распределения реакций для
статически
неопределимой
схемы
крепления
барабана
лебедки,
изображенной на рис. 4.18.
Рис. 4.18 Схема крепления барабана лебедки
На
рис.
4.18
введены
следующие
обозначения:
N 0 ,
–
реакция
большого/малого опорного ролика; R0, – радиус большого/малого опорного
ролика; rбар – внутренний радиус барабана.
Для определения распределения реакций воспользуемся методом
штрафных функций с моделью контактного взаимодействия согласно
контактной теории Герца [96-98]. В случае контакта между двумя
цилиндрами с параллельными осями сила нормальной реакции
F
определяется согласно контактной теорией Герца для контактной пары
«цилиндр-цилиндр»:
68
F
Здесь
l
–
длина
E **l
4
контактирующих
(4.9)
d.
цилиндров,
d
–
сближение
соприкасающихся тел. Заметим, что в уравнении (4.9) сила нормальной
реакции F прямо пропорциональна сближению соударяющихся тел d .
Фактически это означает, что соотношение сближений соударяющихся тел
определяет соотношение реакций опорных роликов.
Применим метод штрафных функций к статически неопределимой
схеме крепления барабана лебедки и рассмотрим малое смещение барабана в
направлении опорных роликов (см. рис. 4.19) для определения соотношения
сближений соударяющихся тел по большому и малому опорным роликам.
Рис. 4.19 Применение метода штрафных функций для определения сближения
соударяющихся тел
На рис. 4.19 введены следующие обозначения: d 0, – сближение
большого/малого опорного ролика с барабаном;
точка пересечения внешней поверхности
внутренней поверхностью барабана.
69
Pin,ex
малого
– внутренняя/внешняя
опорного
ролика
с
Будем решать данную задачу двумя независимыми подходами:
аналитически,
путем
явного
определения
искомых
сближений
соударяющихся тел и численно с применением программной системы
MSC.ADAMS [99].
Рассмотрим аналитический подход. В виду того, что фактически,
требуется определить длину отрезка
Pin Pex ,
определяемого взаимным
расположением окружностей и прямых линий, наиболее удобно записать
математические
уравнения,
принадлежащих
этим
описывающие
фигурам
и
найти
расположение
искомый
отрезок
точек,
Pin Pex
алгебраически. Запишем следующие уравнения в системе координат XZ с
началом отчета в смещенном центре барабана, состоящие из: уравнения
прямой линии, проходящей через начало координат и центр малого опорного
ролика; уравнения окружности – внутренней поверхности смещенного
барабана; уравнения окружности – внешней поверхности малого опорного
ролика.
z
r
бар.
r
R cos d0
x,
R sin
2
x z 2 rбар
.,
2
x rбар. R sin z d rбар. R cos
бар.
2
2
(4.10)
R .
2
В результате решения системы уравнений, состоящих из первого и второго
уравнений системы (4.10) можно найти координаты точки
Pin. .
В результате
решения системы уравнений, состоящих из первого и третьего уравнений
системы (4.10) можно найти координаты точки
Pin.
и
Pex.
Pex. .
Зная координаты точек
можно найти расстояние между ними по известной формуле.
Рассмотрим
численный
подход.
На
рис.
4.20
представлена
разработанная MBD (Multibody dynamic) модель. В качестве численного
метода контактного взаимодействия используется Impact-метод с контактной
жесткостью k
E ** l
4
и показателем e 1 согласно контактной теории Герца
для контактной пары «цилиндр-цилиндр».
70
Рис. 4.20 MBD модель для определения распределения реакций схемы крепления барабана
лебедки
В
результате
моделирования
статически
неопределимой
схемы
крепления барабана лебедки определено распределение реакций в опорных
роликах аналитическим и численным подходами.
Проведем анализ чувствительности геометрических параметров схемы
крепления на отношение реакций в опорных роликах. Представим
результаты в безразмерном виде для удобства их экстраполяции к любым
габаритным размерам лебедки.
В первую очередь, проведем верификацию на примере варьирования
угла установки малого опорного ролика (см. рис. 4.21). Отметим, что
наблюдается совпадение результатов, а также наличие высокой зависимости
отношения реакций от угла установки малого опорного ролика .
Рассмотрим
также
чувствительность
решения
от
относительных
размеров опорных роликов и величины сближения по большому опорному
ролику (см. рис. 4.22). Отметим, что наблюдается низкая зависимость
отношения реакций от этих параметров варьирования.
71
Рис. 4.21 Сопоставление аналитического и численного подходов на примере варьирования
угла установки малого опорного ролика
Рис. 4.22 Определение чувствительности решения от относительных размеров опорных
роликов и величины сближения по большому опорному ролику
72
4.5.2 Определение момента силы, необходимого для запуска лебедки
Для эксплуатации лебедки согласно заданному скоростному режиму
работы необходимо, в первую очередь, обеспечить запуск лебедки из
состояния покоя. На рис. 4.23 изображены основные механические
параметры лебедки.
Рис. 4.23 Основные механические параметры лебедки
Как и в п. 4.5.1, будем решать данную задачу аналитически и численно.
Представим
только
некоторые
основные
промежуточные
выражения
аналитического решения. Передаточные соотношения для механизма,
представленного на рис. 4.23, имеют следующий вид:
M вал
Rвал
,
М бар. Rбар. cos
M бар.0 rбар
,
М тр.0
R0
M бар. rбар.
.
М тр.
R
73
(4.11)
Распишем более подробно моменты сил трения на большом и малом
опорных роликах при допущении пропорционального разделения силы
зацепления по опорным роликам:
N0
~
r0 ,
M тр.0 N 0 r0 N 0 Fц . з. sin
N 0 2 N
N
1
~
r .
M тр. N r N Fц. з. sin
N 0 2 N cos
(4.12)
Моментом силы, необходимым для запуска лебедки, является сумма
момента силы трения опорных роликов приведенная к валу и собственного
момента трения:
M пуска
rбар. ~
Rвал rбар. ~
N 0 r0 2
N r mбар. g Fвал rвал ,
Rбар. cos R0
R
(4.13)
где – угол зацепления, – статический коэффициент трения, –
эксцентриситет барабана. Подставляя выражение (4.12) в формулу (4.13)
нетрудно получить финальное выражение для аналитического определения
момента пуска. Отметим несколько особенностей: сила зацепления
Fц .з.
определяется из передаточных соотношений (4.11) исходя из равенства
момента силы со стороны барабана и со стороны приводного вала; сила на
валу
Fвал
определяется как суперпозиция веса вала и силы зацепления.
Представим результаты в безразмерном виде для удобства их
экстраполяции к любым габаритным размерам лебедки:
~
M пуска
M пуска
mбар. grвал
.
(4.14)
На рис. 4.24 представлена зависимость угловой скорости приводного
вала при линейно возрастающем моменте силы на валу, полученная двумя
подходами: аналитически и численно. Отметим, что при численном решении
момент пуска задавался кусочно-линейной функцией с выходом на
постоянную величину момента силы для более точного определения начала
движения и как следствие – определения момента запуска барабана лебедки.
74
Рис. 4.24 Зависимость угловой скорости приводного вала при линейно возрастающем
моменте силы на валу
Анализируя результаты, представленные на рис. 4.24, заметим малую
осцилляцию зависимости угловой скорости приводного вала при времени
около 0.2 с. Этот эффект вызван достижением момента силы на валу
величины собственного диссипативного сопротивления вала, а также
наличием
небольшого
зазора
в
зубчато-цевочной
передаче.
После
преодоления этого зазора приводной вал упирается в барабан лебедки и
останавливается вплоть до преодоления диссипативного сопротивления от
барабана.
Также отметим, что в численной модели использовалась модель
трения, допускающая переход от статического коэффициента трения к
динамическому, соответствующая переходу полужидкостного трения в
жидкостное трение согласно кривой Герси-Штрибека [100,101]. При этом
конкретное
значение
Фогельполя
[102,103].
скорости
Именно
перехода
этим
вычислялось
фактором
по
формуле
определяется
излом
зависимости угловой скорости приводного вала при времени около 3.9 с.
75
4.5.3 Определение момента силы и потребляемой мощности электродвигателя при
вращении барабана лебедки
После
преодоления
статических
сил
сопротивления
возможна
эксплуатация лебедки согласно заданному скоростному режиму работы
барабана лебедки. Рассмотрим в качестве скоростного интервала участок с
линейно возрастающей скоростью. Также как и в предыдущих п. 4.5.1 и 4.5.2
будем решать данную задачу аналитически и численно.
Кинетическая энергия T механизма лебедки представлена выражением
(4.15), а потенциальная П – выражением (4.17).
1
T J вал 2 J бар. A вал ,
i
2
2
где
J вал
(4.15)
– момент инерции вала; – КПД зубчато-цевочной передачи [104]; i
– передаточное отношение; A – коэффициент инерции опорных роликов (см.
формулу (4.17)); вал – угловая скорость вращения приводного вала.
2
2
2 rбар.
1 r
J 2 бар. J 0 ,
A 2
i R
i R0
(4.16)
где J – момент инерции малого опорного ролика; J 0 – момент инерции
большого опорного ролика.
П
mбар. g
1 cos вал П0 ,
i
где вал – угол поворота вала;
П0
(4.17)
– начальный уровень потенциальной
энергии. Подставляя значения T и П в уравнение Лагранжа 2-го рода,
нетрудно выразить внешний момент силы на валу
M внеш . ,
необходимый для
реализации скоростного режима работы:
d T T
П
M внеш . M тр. ,
dt вал вал
вал
где
M тр.
(4.18)
– суммарный момент трения, определяемый формулой
N rбар.
2 N rбар.
M тр. M тр.вал 0
r0
r .
i
R
i
R
0
76
(4.19)
Потребляемая мощность электродвигателя
эксплуатационный
скоростной
режим,
Pвнеш . ,
вычисляется
обеспечивающего
по
определению
мощности в соответствии с выражением для суммарного момента трения
M тр.
:
Pвнеш . M внеш .вал .
(4.20)
Как и ранее, представим приведенные характеристики
~
M внеш .
и
~
Pвнеш .
для
удобства экстраполяции к любым габаритным размерам:
M внеш .
~
M внеш .
,
вал
J бар.
i
~
~
Pвнеш . M внеш .вал .
(4.21)
На рис. 4.25 представлено сравнение различных подходов к определению
приведенного момента силы на валу: в результате численного решения до
сглаживания (расчетная кривая, желтая); в результате численного решения
после
сглаживания
(сглаженная
кривая,
красная);
в
результате
аналитического решения в соответствие с вышеприведенными формулами
(аналитическая оценка, черная).
Рис. 4.25 Сравнение различных подходов к определению приведенного момента силы на
валу
77
Отметим, что способ закрепления барабана лебедки не позволяет
использовать при численном моделировании тривиальные типы передачи
движения (например, связи типа gear или couple). Поэтому передача
движения осуществляется напрямую, через решение непосредственной
контактной взаимосвязи с применением Impact-метода (связь типа contact).
На рис. 4.25 для расчетной кривой видно наличие высокочастотных
осцилляций, связанных с периодическим соударением зубьев о цилиндры в
цевочном зацеплении. Однако амплитуда этих колебаний завышена в виду
следующих факторов: наличие коллокационного, а не распределенного
контакта; использование модели абсолютно жесткого, а не деформируемого
тела; зазор в зубчато-цевочной передачи, формирующий непостоянное
контактное взаимодействие.
В связи с наличием вышеописанных факторов целесообразно проводить
сглаживание результатов, для чего автором была разработана и применена
программа по сглаживанию результатов на языке Matlab MathWorks.
Процедура сглаживания состоит в применении частотного фильтра для
осреднения набора точек вокруг рассматриваемой. Другими словами, чтобы
получить сглаженное значение для кривой, необходимо сложить все
значения расчетной кривой в окрестности рассматриваемой точки и поделить
на их количество. Окрестность выбирается симметрично и характеризуется
т.н. радиусом осреднения – числом точек, входящих в осреднение с каждой
из сторон. Весь интервал невозможно сгладить фиксированным радиусом
осреднения, ведь вблизи границ интервала не существует требуемого
количества точек. В этом случае радиус осреднения выбирается равным
максимально возможному числу точек со стороны границы. Граничные
точки расчетной и сглаженной кривых совпадают. Повторяя эту процедуру
несколько раз, можно получить более плавную кривую. Число повторов
характеризуется глубиной осреднения. Параметры сглаживания (радиус и
глубина осреднения) выбираются, исходя из частоты движения, связанной с
эксцентриситетом и шагом интегрирования.
78
Однако впоследствии были обнаружены аналогичные возможности у
вложенных функций в Matltab MathWorks и в OriginLab. После проведения
процедуры сглаживания численных результатов для расчетной кривой можно
наблюдать
совпадение
получившиеся
с
аналитическим
периодические
решением.
колебания
Отметим,
обусловлены
что
наличием
эксцентриситета у барабана лебедки.
Вариация
существенно
коэффициента
влияет
на
динамического
изменение
трения
скольжения
исследуемых
динамических
характеристик, поэтому целесообразно использовать подшипники качения
вместо подшипников скольжения. Существует широкий диапазон оценок
эффективного коэффициента трения от 0.0018 до 0.02 в литературе,
обусловленный многочисленными факторами: условия нагружения, вид
подшипника, условия смазки, температура, скоростной режим и др.
[105,106].
Будем рассматривать диссипативную модель узла подшипника с
помощью
четырех
различных
подходов:
трение
скольжения
с
коэффициентом трения равным 0 (этот подход соответствует идеальной
модели подшипника); трение скольжения с коэффициентом трения равным
0.0018 (этот подход соответствует оценочной модели современного
подшипника,
работающего
в
комфортных
условиях
эксплуатации);
эмпирическая модель трения подшипников фирмы SKF (этот подход
соответствует модели трения SKF, сформированной на экспериментальных
данных, однако возможно, что эта модель корректно описывает поведение
только подшипников фирмы SKF); модель трения подшипника из модуля
MSC.ADAMS/Machinery,
позволяющая
включить
в
механизм
параметризованные подшипники из вложенной библиотеки реальных
подшипников ведущих мировых производителей (этот подход обладает
максимальной информацией об используемом подшипнике, однако в то же
время,
имеет
варьируемые
коэффициенты
собственно и определяют диссипативную модель).
79
демпфирования,
которые
Эмпирическая модель момента трения подшипника
M подш.
от шведской
фирмы-производителя SKF, представлена формулой [105]:
M подш. нагр. гол.M кач. M ск. M упл. M сопр..
(4.22)
Здесь нагр. – коэффициент, зависящий от нагревания смазки;
коэффициент
голодания;
M упл.
–
уменьшения
M кач. –
смазочной
момент трения качения;
пленки
M ск . –
в
режиме
гол.
–
смазочного
момент трения скольжения;
– момент трения, вызванный наличием уплотнений в подшипнике; M сопр.
момент
трения
сопротивления,
связанного
с
гидродинамическим
движением смазочного материала.
На рис. 4.26 и 4.27 представлено сравнение вышеописанных способов
описания подшипника на примере момента силы и потребляемой мощности
приводного вала лебедки, полученные численно в результате сглаживания
расчетных зависимостей.
Рис. 4.26 Сравнение различных способов описания подшипника качения на примере
приведенного момента силы приводного вала лебедки
80
Рис. 4.27 Сравнение различных способов описания подшипника качения на примере
приведенной потребляемой мощности приводного вала лебедки
Вполне ожидаемо было увидеть как колебания вокруг 0 для модели
трения скольжения с коэффициентом равным 0, так и увеличение
исследуемых динамических характеристик при увеличении коэффициента
трения скольжения . Показания эмпирической модели трения SKF близки к
модели трения скольжения с коэффициентом равным 0.0018. Любопытной
моделью
узла
подшипника
является
модель
из
модуля
MSC.ADAMS/Machinery. Эта модель обладает самыми полными сведениями
о геометрических данных подшипника и, вероятно, целиком воспроизводит
процесс движения подшипника с сопутствующим вращением роликов.
Очевидно, что эта модель также обладает полной информацией об условиях
нагружения и скорости вращения подшипника, но гидродинамические
эффекты, связанные с движением рабочей жидкости не учитываются в этой
модели.
81
4.5.4 Определение влияния рассогласованного динамического движения на
характеристики укладки кабеля
Проведем моделирование укладки кабеля с учетом нелинейного
эффекта
рассогласованного
динамического
движения,
вызванного
сматыванием кабеля с барабана лебедки на базе математическое модели,
описание которой представлено в п. 4.1. Рассмотрим 10% рассогласование
скорости движения судна-кабелеукладчика и скорости схода кабеля. Очень
важно понимать, что рассогласованное движение – нестационарное.
Вследствие этого будет отсутствовать установившаяся форма динамического
равновесия, как это было ранее на рис. 3.7, 3.10 и 4.16. В виду того, что весь
процесс укладки кабеля при рассогласованном движении – переходный
процесс, то форма, приобретаемая кабелем, определяется величиной
рассогласования
и
временем,
в
течение
которого
проводилась
рассогласованная укладка кабеля. На рис. 4.28 представлено влияние
рассогласованного динамического движения на форму кабеля при его
укладке под водой.
Анализируя результаты, представленные на рис. 4.28, можно заметить,
что при замедленном режиме укладки кабеля, форма кабеля постепенно
приобретает вертикальное положение, а в случае ускоренного режима –
вначале увеличивается провисающий участок кабеля, а в дальнейшем кабель
приобретает более вытянутую и линейную форму.
Рис. 4.28 Влияние времени, в течение которого проводилась рассогласованная укладка
кабеля на приобретаемую форму кабеля
82
4.6 Выводы
Разработана математическая модель движения кабеля при его укладке
под
водой
с
учетом
микроструктуры
кабеля
и
рассогласованного
динамического движения, обусловленного сматыванием кабеля с барабана
лебедки.
Разработана
нестационарного
интегрированная
пространственного
программа
нагружения
по
определению
кабеля
на
языке
программирования Fortran.
Достоверность
разработанной
математической
модели
и
интегрированной программы подтверждается соблюдением многочисленных
верификационных проверок на протяжении всей главы диссертации. В
частности,
речь
идет
о:
совпадении
результатов
разработанной
математической модели с известным аналитическим решением о натяжении
якорной цепи; соблюдением сходимости конечно-элементной модели
укладки кабеля по количеству конечных элементов, приходящихся на
провисающий участок кабеля; соблюдением сходимости конечно-элементной
модели кабеля при определении эффективных коэффициентов жесткости
кабеля; совпадением результатов статического провисания кабелей с
известным уравнением цепной линии для случая, когда эффективная
жесткость кабеля на изгиб недостаточна для выхода за рамки теории
механики нити; совпадением результатов аналитического решения с
численным для задач о статическом распределении реакций схемы крепления
барабана лебедки, об определении момента силы, необходимого для запуска
лебедки, об определении момента силы и потребляемой мощности
электродвигателя при вращении барабана согласно заданному скоростному
режиму работы.
Определены обобщенные коэффициенты жесткости для силовых
подводных кабелей следующих марок: трехжильный кабель с одиночным
бронированием
2XS2YRAA;
одножильный
кабель
с
одиночным
бронированием ZS-YJQ41; одножильный кабель с двойным бронированием
83
GASLMLTV. Полученные обобщенные коэффициенты жесткости включены
в математическую модель движения кабеля при его укладке под водой.
Определены: распределение реакций для схемы крепления барабана
лебедки; момент силы необходимый для запуска лебедки; момент силы и
потребляемая мощность электродвигателя при вращении барабана лебедки
согласно скоростному режиму работы. Аналитическое решение данных задач
выявило величину влияния механических параметров на характеристики
электродвигателя.
Получены
результаты
для
рассогласованного
динамического движения кабеля при его укладке, которые могут быть
обусловлены работой лебедки.
Использование разработанной математической модели движения
кабеля при его укладке под водой в совместности с разработанной
интегрированной
программой
нестационарного
пространственного
нагружения кабеля допускает учет таких нелинейных эффектов как: изгиб и
кручение кабеля, а также рассогласованное движение кабеля. Время
моделирования для данной нелинейной модели зависит от планируемой к
вытравливанию длины кабеля. В частности, моделирование укладки 200 км
кабеля
согласно
заданному
скоростному
режиму
потребует
более
длительного времени моделирования – около суток расчетного времени,
причем на более мощных вычислительных станциях по сравнению с
предыдущей главой диссертации.
Вместе с этим, не все физические явления были учтены в настоящей
диссертации, оставляя поле для дальнейших исследований. Для учета этих
явлений рекомендуется расширять математическую модель настоящей главы
диссертации.
84
Заключение
Исследована задача о разработке цифровой модели укладки кабеля под
водой, позволяющей учитывать многочисленные физические явления,
имеющие место при проведении реальных морских работ. Определена форма
и натяжение провисающего участка кабеля при его укладке под водой
согласно одному из существующих способов укладки кабеля. Основные
полученные научные и технические результаты представлены в следующих
положениях:
1. Разработана математическая модель динамического равновесия
кабеля при его укладке под водой с учетом растяжимости кабеля согласно
закону Гука и влияния профиля подводного течения, изменчивого по
глубине.
2. Разработан комплекс программ для инженерной оценки формы и
натяжения провисающего участка кабеля на языке программирования Matlab
MathWorks.
3. Разработана математическая модель движения кабеля при его
укладке под водой с учетом микроструктуры кабеля и рассогласованного
динамического движения, обусловленного сматыванием кабеля с барабана
лебедки.
4.
Разработана
программа
по
определению
нестационарного
пространственного нагружения кабеля на языке программирования Fortran,
интегрированная в конечно-элементный программный комплекс.
В рамках дальнейших исследований наибольший интерес представляет
учет неровного профиля морского дна и пространственного движения кабеля
для математической модели движения кабеля при его укладке под водой.
85
Список литературы
[1] ICPC Ltd., Submarine Cables as a Sustainable Use of the Deep Sea
Environment. 2018. 24 p. URL: www.iscpc.org (дата обращения: 12.09.2019).
[2] ICPC Ltd., Submarine Cables and Biodiversity Beyond National Jurisdiction.
2016. 47 p. URL: www.iscpc.org (дата обращения: 12.09.2019).
[3] www.submarinecablemap.com
[4] www.iscpc.org/information/cableships-of-the-world
[5] web.asn.com
[6] globalmarine.group
[7] www.subcom.com
[8] www.ktsubmarine.co.kr/eng/main.asp
[9] www.emarine.ae
[10] www.maersksupplyservice.com/Pages/Home.aspx
[11] ICPC Ltd., About Submarine Power Cables. 2011. 45 p. URL: www.iscpc.org
(дата обращения: 12.09.2019).
[12] ICPC Ltd., About Submarine Telecommunications Cables. 2011. 52 p. URL:
www.iscpc.org (дата обращения: 12.09.2019).
[13] Пешков И.Б. Подводные кабели: современное состояние и тенденции
развития. Обзор // Кабели и провода. 2013. №5 (342). С. 9-15.
[14] ICPC Ltd., Submarine cables and the oceans: connecting the world. 2009. 68
p. URL: www.iscpc.org (дата обращения: 12.09.2019).
[15] Donaghy R., HV Submarine Cable Systems Design, Testing and Installation.
20010. 32 p. URL:
www.engineersireland.ie/public/cigre/session3_robert_donaghy.pdf (дата
обращения: 12.09.2019).
[16] Немов А.С., Войнов И.Б., Боровков А.И. Расчетное определение
жесткостных характеристик кабелей с иерархической структурой // Научнотехнические ведомости СПбГПУ. 2008. № 6 (70). С. 21-27.
86
[17] Nemov A.S., Voynov I.B., Borovkov A.I., Boso D.P., Schrefler B.A.
Generalized stiffness coefficients for ITER superconducting cables, direct FE
modeling and initial configuration // Cryogenics. 2010. Vol. 50. No. 5. P. 304-314.
[18] Zhang D., Ostoja-Starzewski M. Finite element solutions to the bending
stiffness of a single-layered helically wound cable with internal friction // Journal
of Applied Mechanics. 2016. Vol. 83. No. 3. 8 p.
[19] Yu Y., Wang X., Chen Z. A simplified finite element model for structural
cable bending mechanism // International Journal of Mechanics Sciences. 2016.
Vol. 113. P. 196-210.
[20] Leong K.H., Latiff R.H.A., Yusof F., Ooi. C.C., Rahman M.R.A. Intermittent
audio failure analysis of a remote speaker-microphone for a two-way radio // J.
Fail. Anal. and Preven. 2016. Vol. 16. P. 75-85.
[21] Luz F.F., de Menezes E.A.W., da Silva L.V., Cimini C.A. Jr., Amico S.C.
Strength analysis of composite cables // Lat. Am. j. solids struct. 2018. Vol. 15.
No. 4. 9 p.
[22] Karahan M., Kalenderli Ö. Coupled electrical and thermal analysis of power
cables using finite element method. Heat Transfer – Engineering Applications.
2011. P. 205-230.
[23] Xu Z., Hu Z., Zhao L., Zhang Y., Yang Z., Hu S., Li Y. Application of
temperature field modeling in monitoring of optic-electric composite submarine
cable with insulation degradation // Measurement. 2019. Vol. 133. P. 479-494.
[24] Holyk C., Liess H.-D., Grondel S., Kanbach H., Loos F. Simulation and
measurement of the steady-state temperature in multi-core cables // Electric Power
Systems Research. 2014. Vol. 116. P. 54-66.
[25] Dubitsky S., Greshnyakov G., Korovkin N. Comparison of finite element
analysis to IEC-60287 for predicting underground cable ampacity // EnergyCon
2016 IEEE Int. Conference, Leuven, Belgium 4-8 Apr. 2016. 6 p.
[26] Del-Pino-López J.C., Hatlo M., Cruz-Romero P. On simplified 3D finite
element simulations of three-core armored power cables // Energies. 2018. Vol. 11.
No. 11. 14 p.
87
[27] Керестень И.А., Корнилова Е.В., Михайлов А.А. Конечно-элементное
определение
эффективных
коэффициентов
жесткости
для
силовых
подводных кабелей с гетерогенной структурой // Морские интеллектуальные
технологии. 2019 №3 (45) том 3. С. 208-215.
[28] Routh E. An elementary treatise on the dynamics of a system of rigid bodies.
With numerous examples. London. 1860. 336 p.
[29] Appel P. Traite de mecanique rationnelle, t. I. Paris. 1893. 549 p.
[30] Минаков А.П. Основы механики нити. Научно-исследовательские труды
Московского текстильного института. 1941. Т. 9, Вып. 1. 88 с.
[31] Алексеев Н.И. Статика и установившееся движение гибкой нити. М.:
Легкая индустрия. 1970. 270 с.
[32] Меркин Д.Р. Введение в механику гибкой нити. М.: Наука, 1980. 240 с.
[33] Керестень И.А., Михайлов А.А., Войнов И.Б., Боровков А.И. Численное
моделирование укладки растяжимого кабеля на дно моря с движущегося
судна с учетом гидродинамических сил сопротивления воды // Подводные
исследования и робототехника. 2019 №1 (27). С. 12-20.
[34] Керестень И.А., Михайлов А.А., Войнов И.Б., Боровков А.И. Алгоритм
численного
определения
формы
и
натяжения
провисающей
части
растяжимой гибкой нити при сматывании с движущейся катушки под водой
на
горизонтальную
плоскость
//
XII
Всероссийский
съезд
по
фундаментальным проблемам теоретической и прикладной механики.
Аннотации докладов. – Уфа: РИЦ БашГУ, 2019 С. 19.
[35] Крылов А.Н. О равновесии шаровой мины на течении / А. Крылов. –
Санкт-Петербург : тип. Мор. м-ва. 1909. 27 с.
[36] Foppl O. Mitteilungen aus der Cottinger Modellversuchsanstalt. Zeitschrift für
Flugtechnik und Motorluftschiffahrt, 29 okt., 1910.
[37] Prandtl L. Ergebnisse der Aerodynamischen Versuchsanstalt zu Gottingen,
Lieferung II. München und Berlin, 1923.
[38] Eiffel G. Nouvelles researches sur la resistance de l'air et l'aviation. Paris,
1914.
88
[39] Relf E. and Powell C. Thest on smooth and stranded wires inclined to the
wind direction, and a comparison of results in air and water. Reports and
Memoranda of Britich Advisory Comity for Aeronautics, N 307, 1917.
[40] Кузнецов Б.Я. Лобовое сопротивление тросов, проволок, тандеров и
авиационных лент, Труды ЦАГИ, 1931, вып. 97.
[41] Чесалов А.В. Коэффициенты вредных сопротивлений самолета, Труды
ЦАГИ, 1929, вып. 42.
[42] Кузнецов Б.Я. Аэродинамические исследования цилиндров, Труды
ЦАГИ, 1931, вып. 98.
[43] Матросов И.Р. Теоретические основы для расчета движения судна с
тралом // Рыбное хозяйство. 1961. № 6. C. 41-53.
[44] Чубаров С.Н. О влиянии гидродинамических сил на ваер и выбор тросов
для скоростного и глубоководного траления // Рыбное хозяйство. 1962. № 5.
C. 45-49.
[45] Алексеев Н.И. О натяжении и пространственной форме канатов в потоке
воды // Труды ВНИИ морского рыбного хозяйства и океанографии. 1966. Т.
61. С. 277-285.
[46] Ren, Z. and Skjetne, R An on-site current profile estimation algorithm for a
moored floating structure // International federation of automatic control. 2016.
Vol. 49, No. 23. P. 153-158.
[47] Ren, Z. and Skjetne, R. A tension-based position estimation solution of a
moored structure and its uncertain anchor positions // International federation of
automatic control. 2016. Vol. 49, No. 23. P. 251-257.
[48] Ren, Z., Skjetne, R. and Kjerstad,Ø.K. A tension-based position estimation
approach for moored marine vessels // International federation of automatic
control. 2015. Vol. 48, No. 16. P. 248-253.
[49] Ren, Z., Skjetne, R., and Hassani, V. Supervisory control of line breakage for
thruster-assisted position mooring system // International federation of automatic
control. Vol. 48, No. 16. P. 235-240.
89
[50] Park, J. Kim, N. Dynamics of a semi-submersible autonomous underwater
vehicle with a towfish towed by a cable // International Journal of Naval
Architecture and Ocean Engineering. 2015. Vol. 7, No. 2. P. 409-425.
[51] Vu, M.T., Choi, H.S., Kang, J.I., Ji, D.H., Jeong, S.K. A study of hovering
motion of the underwater vehicle with umbilical cable // Ocean Engineering. 2017.
Vol. 135. P. 137-157.
[52] Костенко В.В., Мокеева И.Г. Исследование влияния кабеля связи на
маневренность
телеуправляемого
подводного
аппарата
//
Подводные
исследования и робототехника. 2009. №1 (7). С. 22-27.
[53] Костенко В.В., Львов О.Ю. Комбинированная система связи и навигации
автономного подводного робота с поплавковым модулем // Подводные
исследования и робототехника. 2017. №1 (23). С. 31-43.
[54] Ваулин Ю.В., Костенко В.В., Мокеева И.Г., Матвиенко Ю.В., Рылов
Н.И. Особенности координирования донных источников навигационных
сигналов с использованием буксируемого антенного модуля // Подводные
исследования и робототехника. 2018. №2 (26). С. 4-11.
[55] Киреев В.И., Пантелеев А.В. Численные методы в примерах и задачах:
Учебное пособие. М.: Изд-во МАИ, 2000. 376 с.
[56] Потемкин В.Г. Введение в MATLAB. М.: Диалог-МИФИ. 2000.
[57] Керестень И.А., Войнов И.Б. Конечно-разностная программа укладки
кабеля под водой // Свидетельство о регистрации программы для ЭВМ №
2019613644, Роспатент, М., 20.03.19.
[58] Huang S. Dynamic analysis of three-dimensional marine cables // Ocean
Engineering. 1994. Vol. 21, No. 6. P. 587-605.
[59] Aamo, O. and Fossen, T. Finite element modeling of moored vessels //
Mathematical and Computer Modeling of Dynamical Systems. 2001. Vol. 7, No. 1.
P. 47-75.
[60] Matulea, I.C., Năstase, A.T., Tălmacia, N., Slămnoiu, G., Gonçalves-Coelho,
A.M. On the equilibrium configuration of mooring and towing cables // Applied
Ocean Research. 2008 Vol. 30, No. 2. P. 81-91.
90
[61] Choc, Y.I. and Casarella, M.J. Hydrodynamic resistance of towed cables //
Journal of Hydronautics. 1971. Vol. 5, No. 4. P. 126-131.
[62] Johannessen J.A., Pripp T., Eldevik T. GOCE studies of mean dynamic
topography and ocean circulation in the Nordic Seas. Progress report for: project
number 212020 – GOCE studies of mean dynamic topography and ocean
circulation in the high latitude and Arctic Ocean (GOCE-MDT). 2013. 8 p.
[63]
Мирзоев
Д.А.,
углеводородных
технологиями
Ибрагимов
ресурсов
//
Разработка
И.Э.,
арктики
и
Архипова
О.Л.
инновационными
эксплуатация
нефтяных
Освоение
подводными
и
газовых
месторождений. 2012. №3. С. 49-53.
[64] Сегерлинд Л. Применение метода конечных элементов: Пер. с англ. М.:
Мир, 1979. 392 с.
[65] Зенкевич О., Морган К., Конечные элементы и аппроксимация: Пер. с
англ. М.: Мир, 1986. 318 с.
[66] Бате К.-Ю. Методы конечных элементов. М.: Физматлит, 2010. 1024 с.
[67] Галлагер Р. Метод конечных элементов. Основы: Пер. с англ. М.: Мир,
1984.428 с.
[68] Зенкевич О. Метод конечных элементов в технике: Пер. с англ. М.: Мир,
1975. 541 с.
[69] Деклу Ф. Метод конечных элементов: Пер. с фр. М.: Мир, 1976. 94 с.
[70] Dassault Systems ABAQUS 6.13 Online documentation. – URL:
http://dsk.ippt.pan.pl/docs/abaqus/v6.13/index.html
(дата
обращения:
12.09.2019).
[71] Тимошенко С.П., Янг Д.Х., Уивер У. Колебания в инженерном деле. М.:
Физматлит, 1985. 474 с.
[72] Глушко М.Ф. Стальные подъемные канаты. Киев: Техника, 1966. 328 с.
[73] Сергеев С.Т. Надежность и долговечность подъемных канатов. Киев:
Техника, 1968. 238 с.
[74] Taormina B., Bald J., Want A.,Thouzeau G., Lejart M., Desroy N., Carlier A.
A review of potential impacts of submarine power cables on the marine
91
environment: knowledge gaps, recommendations and future directions //
Renewable and Sustainable Energy Reviews. 2018. Vol. 96. P. 380-391.
[75]
Nexans
Inc.
Submarine
Power
–
Cables.
URL:
https://www.nexans.ru/Germany/2013/SubmPowCables_FINAL_10jun13_engl.pd
f (дата обращения: 12.09.2019).
[76] Ventikos N.P., Stavrou D.I. Submarine power cables: Laying procedure, the
fleet and reliability analysis // Journal of Marine Engineering and Technology.
2013. Vol. 12. No. 1. P. 13-26.
[77]
Hexatronic
Group.
Submarine
cable
–
systems.
URL:
https://hexatronic.com/media/233536/hexatronic_submarine.pdf (дата обращения:
12.09.2019).
[78] Vise S., Adnitt C., Stanisland R. Review of cabling techniques and
environmental effects applicable to the offshore wind farm industry (BERR
Technical Report). 2008. 164 p.
[79] Белоруссов Н.И., Саакян А.Е., Яковлева А.И. Электрические кабели,
провода
и
шнуры:
справочник
–
5-е
изд.,
перераб.
и
доп.
М.:
Энергоатомиздат, 1987. 536 с.
[80] Золоторевский В.С. Механические свойства металлов. 3-е изд., перераб.
и доп. М.: МИСИС, 1998. 400 с.
[81] Крыжановский В.К., Бурлов В.В., Паниматченко А.Д., Крыжановская
Ю.В. Технические свойства полимерных материалов: учеб-справ. пособие –
2-е изд., испр. и доп. СПб.: Профессия, 2005. 248 с.
[82] Лившиц Б.Г., Крапошин В.С., Липецкий Я.Л. Физические свойства
металлов и сплавов: учебник для ВУЗов. М.: Металлургия, 1980. 320 с.
[83] Гун Р.Б. Нефтяные битумы. М.: Химия, 1973. 432 с.
[84] Ярцев В.П., Ерофеев А.В. Эксплуатационные свойства и долговечность
битумно-полимерных композитов: учебное пособие для студентов. Тамбов:
Изд-во. ФГБОУ ВПО «ТГТУ», 2014. 80 с.
[85] Maher A., Bennet T. Evaluation of Poisson’s Ratio for Use in the Mechanistic
Empirical Pavement Design Guide (MEPDG) (Final Report). 2008. 60 p.
92
[86] Чиркин В.С. Теплофизические свойства материалов ядерной техники:
справочник. М.: Атомиздат, 1968. 485 с.
[87] Инструкция по проектированию технологических трубопроводов из
пластмассовых труб (СН 550-82). Госстрой России. М.: ГУП ЦПП, 2000. 63 с.
[88] Керестень И.А., Войнов И.Б., Михайлов А.А., Боровков А.И.
Рационализация процессов пуска и работы катушечного механизма // Неделя
науки СПбПУ : материалы научного форума с международным участием.
Институт прикладной математики и механики. – СПб. : Изд-во Политехн. унта, 2015. – С. 82-84.
[89] Керестень И.А., Войнов И.Б., Михайлов А.А., Боровков А.И. Численное
моделирование, исследование и анализ процессов пуска и работы роторного
механизма // Неделя науки СПбПУ : материалы научного форума с
международным участием. Лучшие доклады. – СПб. : Изв-во Политехн. унта, 2016. – С. 190-194.
[90] Керестень И.А., Войнов И.Б., Михайлов А.А., Боровков А.И.
Идентификация
параметров
контактной
взаимосвязи
Impact-метода
MSC.ADAMS, исследование распределения реакций схемы крепления
роторного механизма [электронный ресурс] : [сайт] – Москва, 2016. – 6 с. –
Режим доступа: http://docs.mscsoftware.ru/conf/vuz2016/21_Keresten-tesis.pdf
(23.05.2016).
[91] Первицкий Ю.Д. Расчет и конструирование точных механизмов.
Учебное пособие для вузов. 2-e изд. Л.: «Машиностроение», 1976. 456 с.
[92] Гаврилов А.Н., Щедровицкий С.С. Прикладная метрология, методика
расчета, детали и элементы приборов. Т.2. Ч. 1. М.: МАШГИЗ, 1964. 596 с.
[93] Литвин Ф.Л. Проектирование механизмов и деталей приборов. Ф.Л.
Литвин. Л.: «Машиностроение», 1973, 696 с.
[94] Керестень И.А., Плотников Ф.С., Войнов И.Б., Михайлов А.А., Боровков
А.И. Численное моделирование зацепления цевочной передачи роторного
механизма и исследование чувствительности факторов зацепления к
изменению межосевого расстояния // Неделя науки СПбПУ : материалы
93
научной конференции с международным участием. Институт прикладной
математики и механики. – СПб. : Изд-во Политехн. ун-та, 2017. – С. 115-118.
[95] Керестень И.А. Численная идентификация пределов работоспособности
цевочной передачи роторного механизма и исследование чувствительности
факторов зацепления на изменение межосевого расстояния // Двадцать
вторая Санкт-Петербургская Ассамблея молодых ученых и специалистов:
Сборник тезисов. – СПб. : Изд-во СПбГУПТД, 2017. – С. 154-155.
[96]. V.L. Popov, Contact Mechanics and Friction: Physical Principles and
Applications, Springer-Verlag Berlin Heidelberg, 2010. 362 p.
[97] I.N. Sneddon, The Relation between Load and Penetration in the
Axisymmetric Boussinesq Problem for a Punch of Arbitrary Profile. Int. J. Eng.
Sci., 1965, v. 3, pp. 47-57.
[98] Прочность, устойчивость, колебания. Справочник, том 2. Под ред. И.А.
Биргера и Я.Г. Пановко. Издательство «Машиностроение», Москва, 1968. 463
с.
[99] MSC Inc. Adams/Solver C++ Statements. 2014. 408 p.
[100] Махутов Н.А. Конструкционная прочность, ресурс и техногенная
безопасность. В двух частях. Часть 1. Новосибирск.: Изд. “Наука”, 2005.
[101] Хлебалин Н.А., Костиков А.Ю. Библиотека моделей трения в
SIMULINK
(Опыт
создания
и
использования).
Труды
II
научной
конференции “Проектирование инженерных и научных приложений в среде
MATLAB”. Секция 5. Моделирование в SIMULINK, 2004.
[102] Пожбелко В.И. Аналитическая нелинейная скоростная характеристика
трения и оптимизация толщины смазочного слоя и эксцентриситета
гидродинамических подшипников, // Известия высших учебных заведений.
Машиностроение, 2 (611). 2011. 23-30 c.
[103] Ицкович Г.М., Чернавский С.А., Кисилев В.А., Боков К.Н., БончОсмоловский М.А. Сборник задач и примеров расчета по курсу деталей
машин. 3-е изд. М.: «Машиностроение», 1965. 328 с.
94
[104] Иванов А.С., Муркин С.В. Конструирование современных моторредукторов: Электронное учебное издание. М.: Изд-во МГТУ им. Н.Э.
Баумана. 2012. 147 с.
[105] SKF inc. Общий каталог SKF. 2006, 1129 с.
[106] Дукельский А.И. Справочник по кранам. Т. 2. Л.: «Машиностроение»,
1973. 504 с.
95
Приложение А. Комплекс программ для определения формы и
натяжения провисающего участка кабеля при его укладке под водой
1. Часть программы (препроцессор), предназначенная для пересчета
необходимых входных расчетных данных по имеющимся физическим
входным данным. Указана номенклатура с описанием физических данных и
образцы входных данных в соответствие с литературой.
function [sigma,mu,q0,Ct,lmbdn,lmbdt,U2,U3,xin] =
Cable_Properties(d,Roc,E,V,H,Umax,opposite)
%%%%%%%%%%%%%%%%
% NOMENCLATURE %
%%%%%%%%%%%%%%%%
% g - ускорение свободного падения, м/с^2
% Roh - плотность воды, кг/м^3
% nu - динамическая вязкость воды, кг/м*с
% V - рассматриваемая скорость движения кабеля, м/с
% d - диаметр кабеля, м
% Roc - плотность кабеля, кг/м^3
% E - модуль упругости кабеля, Па
% sigma - Площадь сечения кабеля (представлена формула для случая круглого
сечения), м^2
% mu - погонная масса кабеля, кг/м
% Roeff - эффективная плотность кабеля в воде, кг/м^3
% q0 - погонный вес нерастяжимого кабеля в воде (с учетом силы Архимеда), Н/м
% Rn - число Рейнольдса
% Nuss - число Нуссельта
% Cn - коэффициент нормального гидродинамического сопротивления
% Ct - коэффициент касательного гидродинамического сопротивления
% lmbdn - погонное усилие нормального гидродинамического сопротивления, Н/м
% lmbdt - погонное усилие касательного гидродинамического сопротивления, Н/м
% H - глубина прокладки кабеля, м
% Umax - максимальная скорость течения воды на поверхности (в качестве
профиля используется кубический полином), м/с
% opposite - направление течения: попутное=-1; встречное=1
%%%%%%%%%%%%%%%
% COMMON DATA %
%%%%%%%%%%%%%%%
g = 9.80665;
Roh = 1025;
nu = 0.0013;
% V = 0.51444*1 [m/s] == 1 [knot];
%%%%%%%%%%%%%%%%%%%
% JINMO 2015 DATA %
%%%%%%%%%%%%%%%%%%%
% d = 0.041;
% Roc = 1300;
% E = 7e8;
%%%%%%%%%%%%%%%%%%%
% HUANG 1994 DATA %
%%%%%%%%%%%%%%%%%%%
% d = 0.047;
96
% Roc = 3112.5;
% E = 9e9;
%%%%%%%%%%%%%%%%%%
% AAMO 2001 DATA %
%%%%%%%%%%%%%%%%%%
% d = 0.1003;
% Roc = 5500;
% E = 2.55e10;
%%%%%%%%%%%%%%%%%%%%%
% MATULEA 2008 DATA %
%%%%%%%%%%%%%%%%%%%%%
% d = 0.00599;
% Roc = 7850;
% E = 2.15e11;
%%%%%%%%%%%%%%%%%%%%
% COMMON SOLUTIONS %
%%%%%%%%%%%%%%%%%%%%
sigma = (pi*(d^2))/4;
mu = sigma*Roc;
Roeff = Roc - Roh;
q0 = sigma*Roeff*g;
Rn = (Roh*V*d)/(nu);
Nuss = (0.55*(Rn^0.5)) + (0.084*(Rn^(2/3)));
Cn = 1.1 + (4*(Rn^(-0.5)));
Ct = pi*nu*Nuss;
lmbdn = (1/2)*Cn*Roh*d*(V^2);
lmbdt = Ct*V;
U2 = (3*Umax)/(H^2);
U3 = -(2*Umax)/(H^3);
%U2 = Umax/(H^2);
%U3 = 0;
%U2 = 0;
%U3 = Umax/(H^3);
xin = (1/2)*Cn*Roh*d;
end
2.
Часть программы (решатель Merkin). Решатель определяет форму и
натяжение кабеля по входным физическим данным для линейной системы
дифференциальных уравнений, в соответствие с аналитическим решением
Меркина.
function [X,Y,T,N1,C] =
Merkin_Curve(d,Roc,E,V,H,Umax,opposite,step1,draw_figures)
%%%%%%%%%%%%%%%%
% READING DATA %
%%%%%%%%%%%%%%%%
[sigma,mu,q0,Ct,lmbdn,lmbdt,U2,U3,xin] =
Cable_Properties(d,Roc,E,V,H,Umax,opposite);
%DEFINE SINGULAR POINT;
Azero = acos((-q0+(sqrt((q0^2)+(4*(lmbdn^2)))))/(2*lmbdn));
97
Amove = Azero - step1;
alpha1 = 0:step1:Amove;
N1=length(alpha1);
Amax = alpha1(length(alpha1));
%DEFINE INTEGRABLE FUNCTION;
N2=N1;
step2 = (Amax)/(N2-1);
alpha2 = 0:step2:Amax;
F = zeros(1,length(alpha2));
for i=1:length(alpha2)
N3=N1;
step3 = alpha2(i)/(N3-1);
alpha3 = 0:step3:alpha2(i);
for j=1:((length(alpha3))-1)
dF1 = ((q0*sin(alpha3(j))) - (lmbdt*(1(cos(alpha3(j))))))/((q0*cos(alpha3(j))) - lmbdn*((sin(alpha3(j)))^2));
dF2 = ((q0*sin(alpha3(j+1))) - (lmbdt*(1(cos(alpha3(j+1))))))/((q0*cos(alpha3(j+1))) - lmbdn*((sin(alpha3(j+1)))^2));
F(i) = F(i) + (((dF1+dF2)/2)*step3);
end
end
if (draw_figures)
figure;
plot(alpha2,F);
grid on;
hold on;
title('ЗАВИСИМОСТЬ ИНТЕГРАЛЬНОЙ ФУНКЦИИ ОТ ВЕРХНЕГО ПРЕДЕЛА
ИНТЕГРИРОВАНИЯ');
xlabel('ВЕРХНИЙ ПРЕДЕЛ ИНТЕГРИРОВАНИЯ');
ylabel('ЗНАЧЕНИЕ ИНТЕГРАЛЬНОЙ ФУНКЦИИ');
end
%DEFINE INTEGRATION CONSTANT;
N4=N1;
step4=(Amax)/(N4-1);
alpha4=0:step4:Amax;
J = 0;
for i=1:((length(alpha4))-1)
dJ1 = ((exp(F(i)))*sin(alpha4(i)))/((q0*cos(alpha4(i))) lmbdn*((sin(alpha4(i)))^2));
dJ2 = ((exp(F(i+1)))*sin(alpha4(i+1)))/((q0*cos(alpha4(i+1))) lmbdn*((sin(alpha4(i+1)))^2));
J = J + (((dJ1+dJ2)/2)*step4);
end
C=H/J;
%EXPLICIT SOLUTION;
N5=N1;
step5 = (Amax)/(N5-1);
alpha5 = 0:step5:Amax;
X = zeros(1,length(alpha5));
Y = zeros(1,length(alpha5));
L = zeros(1,length(alpha5));
T = zeros(1,length(alpha5));
for i=1:length(alpha5)
N6=N1;
step6 = alpha5(i)/(N6-1);
alpha6 = 0:step6:alpha5(i);
for j=1:((length(alpha6))-1)
dX1 = ((exp(F(j)))*cos(alpha6(j)))/((q0*cos(alpha6(j))) lmbdn*((sin(alpha6(j)))^2));
98
dX2 = ((exp(F(j+1)))*cos(alpha6(j+1)))/((q0*cos(alpha6(j+1))) lmbdn*((sin(alpha6(j+1)))^2));
X(i) = X(i) + (((dX1+dX2)/2)*step6);
dY1 = ((exp(F(j)))*sin(alpha6(j)))/((q0*cos(alpha6(j))) lmbdn*((sin(alpha6(j)))^2));
dY2 = ((exp(F(j+1)))*sin(alpha6(j+1)))/((q0*cos(alpha6(j+1))) lmbdn*((sin(alpha6(j+1)))^2));
Y(i) = Y(i) + (((dY1+dY2)/2)*step6);
dL1 = (exp(F(j)))/((q0*cos(alpha6(j))) - lmbdn*((sin(alpha6(j)))^2));
dL2 = (exp(F(j+1)))/((q0*cos(alpha6(j+1))) lmbdn*((sin(alpha6(j+1)))^2));
L(i) = L(i) + (((dL1+dL2)/2)*step6);
end
T(i) = exp(F(i));
end
X = C*X;
Y = C*Y;
L = C*L;
T = (C*T) + (mu*(V^2));
if (draw_figures)
figure;
plot(X,Y);
grid on;
hold on;
title('ФОРМА НЕРАСТЯЖИМОЙ ГИБКОЙ НИТИ');
xlabel('X, М');
ylabel('Y, М');
figure;
plot(Y,T);
grid on;
hold on;
title('НАТЯЖЕНИЕ НЕРАСТЯЖИМОЙ ГИБКОЙ НИТИ');
xlabel('Y, М');
ylabel('НАТЯЖЕНИЕ, Н');
end
end
3. Часть программы (решатель Implicit). Решатель определяет форму и
натяжение кабеля по входным физическим данным для линейной системы
дифференциальных уравнений прямым решением – с использованием метода
Гаусса.
function [X,Y,T,alpha2] =
Implicit_Curve(d,Roc,E,V,H,Umax,opposite,step1,draw_figures)
%%%%%%%%%%%%%%%%
% READING DATA %
%%%%%%%%%%%%%%%%
[sigma,mu,q0,Ct,lmbdn,lmbdt,U2,U3,xin] =
Cable_Properties(d,Roc,E,V,H,Umax,opposite);
%DEFINE SINGULAR POINT;
Azero = acos((-q0+(sqrt((q0^2)+(4*(lmbdn^2)))))/(2*lmbdn));
Amove = Azero - step1;
alpha1 = 0:step1:Amove;
N1=length(alpha1);
99
Amax = alpha1(length(alpha1));
%DEFINE MATRIX FORM;
N2=N1;
step2 = (Amax)/(N2-1);
alpha2 = 0:step2:Amax;
A = zeros((2*N2)-2,(2*N2)-2);
B = zeros(1,(2*N2)-2);
for i=1:((length(alpha2))-1)
for j=1:(length(alpha2))
if (i == j)
A(i,j) = ((((q0*(sin(alpha2(i))))-(lmbdt*(1(cos(alpha2(i))))))*((alpha2(i+1))-(alpha2(i))))/((q0*(cos(alpha2(i))))(lmbdn*((sin(alpha2(i)))^2))))+1;
end
if (i+1 == j)
A(i,j) = -1;
end
end
end
for i=(length(alpha2)):((2*(length(alpha2)))-2)
for j=1:(length(alpha2))
if (i == (j + N2 - 1))
A(i,j) = ((sin(alpha2(j)))*((alpha2(j+1))(alpha2(j))))/((q0*(cos(alpha2(j))))-(lmbdn*((sin(alpha2(j)))^2)));
end
end
if (i == (length(alpha2)))
B(i) = 0;
end
if (i == ((2*(length(alpha2)))-2))
B(i) = H;
end
end
for i=(length(alpha2)):((2*(length(alpha2)))-2)
for j=((length(alpha2))+1):((2*(length(alpha2)))-2)
if (i + N2 == j + N2 - 1)
A(i,j) = -1;
end
if (i + N2 == j + 1 + N2 - 1)
A(i,j) = 1;
end
end
end
%IMPLICIT SOLUTION;
TY = A\B';
T = zeros(1,N2);
for i=1:N2
T(i) = TY(i) + (mu*(V^2));
end
Y = zeros(1,N2-1);
Y(1) = 0;
Y(N1) = H;
for i=(N1+1):((2*N1)-2)
Y(i-N1+1) = TY(i);
end
X = zeros(1,N2-1);
L = zeros(1,N2-1);
X(1)=0;
L(1)=0;
for i=1:(N1-1)
X(i+1) = X(i) + (((cos(alpha2(i)))*((alpha2(i+1)) (alpha2(i)))*(TY(i)))/((q0*cos(alpha2(i))) - (lmbdn*((sin(alpha2(i)))^2))));
100
L(i+1) = L(i) + ((((alpha2(i+1)) (alpha2(i)))*(TY(i)))/((q0*cos(alpha2(i))) - (lmbdn*((sin(alpha2(i)))^2))));
end
if (draw_figures)
figure;
plot(X,Y);
grid on;
hold on;
title('ФОРМА НЕРАСТЯЖИМОЙ ГИБКОЙ НИТИ');
xlabel('X, М');
ylabel('Y, М');
figure;
plot(Y,T);
grid on;
hold on;
title('НАТЯЖЕНИЕ НЕРАСТЯЖИМОЙ ГИБКОЙ НИТИ');
xlabel('Y, М');
ylabel('НАТЯЖЕНИЕ, Н');
end
end
4. Часть программы (решатель Modified). Решатель определяет форму и
натяжения кабеля по входным физическим данным для нелинейной системы
дифференциальных уравнений с учетом растяжимости кабеля итерационным
решением – с использованием метода Ньютона.
function [X,Y,T,alpha3,TK,YK,TE,YE,XE,LE] =
Modified_Curve(d,Roc,E,V,H,Umax,opposite,step1,NN,draw_figures)
%%%%%%%%%%%%%%%%
% READING DATA %
%%%%%%%%%%%%%%%%
[sigma,mu,q0,Ct,lmbdn,lmbdt,U2,U3,xin] =
Cable_Properties(d,Roc,E,V,H,Umax,opposite);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE INITIAL APPROXIMATION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step1 = 0.00004;
[X,Y,T,alpha2] = Implicit_Curve(d,Roc,E,V,H,Umax,opposite,step1,0);
N2 = length(alpha2);
TYK = [T-(mu*V^2),Y];
%%%%%%%%%%%%%%%%%%%
% NEWTON SOLUTION %
%%%%%%%%%%%%%%%%%%%
% NN = 3;
W = zeros((2*N2)-2,(2*N2)-2);
F = zeros(1,(2*N2)-2);
alpha3 = alpha2;
for k=2:NN
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REDEFINE SINGULAR POINT %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
CFK = lmbdn*((sigma*E) + (TYK(k-1,length(alpha3))) + (mu*(V^2)));
101
AzeroK = acos((((-q0*sigma*E) + sqrt(((q0*sigma*E)^2) +
(4*(CFK^2))))/(2*CFK)));
AmoveK = AzeroK - step1;
alpha3 = 0:((AmoveK-0)/(N2-1)):AmoveK;
%%%%%%%%%%%%%%%%%%%%
% LEFT UPPER BLOCK %
%%%%%%%%%%%%%%%%%%%%
for i=1:((length(alpha3))-1)
for j=1:(length(alpha3))
if (i == j)
W1 = -q0*sigma*E*cos(alpha3(j));
W2 = lmbdn*sigma*E*((sin(alpha3(j)))^2);
W3 = lmbdn*mu*(V^2)*((sin(alpha3(j)))^2);
W4 = -((alpha3(j+1))(alpha3(j)))*((q0*sigma*E*(sin(alpha3(j))))-(sigma*E*(lmbdt*(1(cos(alpha3(j))))))-(mu*(V^2)*(lmbdt*(1-(cos(alpha3(j)))))));
W5 = (lmbdn*((sin(alpha3(j)))^2))+(((alpha3(j+1))(alpha3(j)))*(lmbdt*(1-(cos(alpha3(j))))));
W6 = -lmbdn*((sin(alpha3(j)))^2);
W(i,j) = (W1+W2+W3+W4) + (2*(W5)*(TYK(k-1,j))) + (W6*(TYK(k1,j+1)));
end
if (i+1 == j)
W7 = -W1;
W8 = -W2;
W9 = -W3;
W10 = W6;
W(i,j) = (W7+W8+W9) + (W10*(TYK(k-1,j)));
end
end
F1 = -q0*sigma*E*(cos(alpha3(i)));
F2 = lmbdn*sigma*E*((sin(alpha3(i)))^2);
F3 = lmbdn*mu*(V^2)*((sin(alpha3(i)))^2);
F4 = -((alpha3(i+1))-(alpha3(i)))*((q0*sigma*E*(sin(alpha3(i))))(sigma*E*(lmbdt*(1-(cos(alpha3(i))))))-(mu*(V^2)*(lmbdt*(1(cos(alpha3(i)))))));
F5 = -F1;
F6 = -F2;
F7 = -F3;
F8 = lmbdn*((sin(alpha3(i)))^2);
F9 = ((alpha3(i+1))-(alpha3(i)))*(lmbdt*(1-(cos(alpha3(i)))));
F10 = -F8;
F(i) = ((F1+F2+F3+F4)*(TYK(k-1,i))) + ((F5+F6+F7)*(TYK(k-1,i+1))) +
((F8+F9)*((TYK(k-1,i))^2)) + (F10*(TYK(k-1,i))*(TYK(k-1,i+1)));
end
%%%%%%%%%%%%%%%%%%%
% LEFT DOWN BLOCK %
%%%%%%%%%%%%%%%%%%%
WU1 = lmbdn*((sin(alpha3(1)))^2)*(TYK(k-1,N2+1));
WU2 = -((alpha3(2))(alpha3(1)))*((sigma*E*sin(alpha3(1)))+(mu*(V^2)*(sin(alpha3(1)))));
WU3 = -((alpha3(2))-(alpha3(1)))*(sin(alpha3(1)));
WU4 = -lmbdn*((sin(alpha3(1)))^2);
W(length(alpha3),1) = (WU1+WU2) + (2*WU3*(TYK(k-1,1))) + (WU4*(TYK(k1,N2+2)));
FU1 = -q0*sigma*E*(cos(alpha3(1)));
FU2 = lmbdn*sigma*E*((sin(alpha3(1)))^2);
FU3 = lmbdn*mu*(V^2)*((sin(alpha3(1)))^2);
FU4 = -FU1;
FU5 = -FU2;
FU6 = -FU3;
FU7 = lmbdn*((sin(alpha3(1)))^2)*(TYK(k-1,N2+1));
FU8 = -((alpha3(2))(alpha3(1)))*((sigma*E*(sin(alpha3(1))))+(mu*(V^2)*(sin(alpha3(1)))));
102
FU9 = -((alpha3(2))-(alpha3(1)))*(sin(alpha3(1)));
FU10 = -lmbdn*((sin(alpha3(1)))^2);
F(length(alpha3)) = ((FU1+FU2+FU3)*(TYK(k-1,N2+1))) +
((FU4+FU5+FU6)*(TYK(k-1,N2+2))) + ((FU7+FU8)*(TYK(k-1,1))) + (FU9*((TYK(k1,1))^2)) + (FU10*(TYK(k-1,N2+2))*(TYK(k-1,1)));
for i=((length(alpha3))+1):((2*(length(alpha3)))-3)
for j=2:(length(alpha3)-2)
if (i == (j + N2 - 1))
W11 = -((alpha3(j+1))(alpha3(j)))*((sigma*E*sin(alpha3(j)))+(mu*(V^2)*(sin(alpha3(j)))));
W12 = -((alpha3(j+1))-(alpha3(j)))*(sin(alpha3(j)));
W13 = lmbdn*((sin(alpha3(j)))^2);
W14 = -W13;
W(i,j) = W11 + (2*W12*(TYK(k-1,j))) + (W13*(TYK(k-1,j+N2))) +
(W14*(TYK(k-1,j+1+N2)));
end
end
F11 = -q0*sigma*E*(cos(alpha3(i-(N2-1))));
F12 = lmbdn*sigma*E*((sin(alpha3(i-(N2-1))))^2);
F13 = lmbdn*mu*(V^2)*((sin(alpha3(i-(N2-1))))^2);
F14 = -F11;
F15 = -F12;
F16 = -F13;
F17 = -((alpha3(i-(N2-1)+1))-(alpha3(i-(N21))))*((sigma*E*(sin(alpha3(i-(N2-1)))))+(mu*(V^2)*(sin(alpha3(i-(N2-1))))));
F18 = -((alpha3(i-(N2-1)+1))-(alpha3(i-(N2-1))))*(sin(alpha3(i-(N21))));
F19 = lmbdn*((sin(alpha3(i-(N2-1))))^2);
F20 = -F19;
F(i) = ((F11+F12+F13)*(TYK(k-1,i+1))) + ((F14+F15+F16)*(TYK(k1,i+2))) + (F17*(TYK(k-1,i-(N2-1)))) + (F18*((TYK(k-1,i-(N2-1)))^2)) +
(F19*(TYK(k-1,i+1))*(TYK(k-1,i-(N2-1)))) + (F20*(TYK(k-1,i+2))*(TYK(k-1,i(N2-1))));
end
WU5 = -lmbdn*((sin(alpha3(N2-1)))^2)*(TYK(k-1,N2+N2));
WU6 = -((alpha3(N2))-(alpha3(N2-1)))*((sigma*E*(sin(alpha3(N21))))+(mu*(V^2)*(sin(alpha3(N2-1)))));
WU7 = -((alpha3(N2))-(alpha3(N2-1)))*(sin(alpha3(N2-1)));
WU8 = lmbdn*((sin(alpha3(N2-1)))^2);
W((2*(length(alpha3)))-2,(length(alpha3))-1) = (WU5+WU6) + (2*WU7*(TYK(k1,N2-1))) + (WU8*(TYK(k-1,N2+N2-1)));
FU11 = q0*sigma*E*(cos(alpha3(N2-1)));
FU12 = -lmbdn*sigma*E*((sin(alpha3(N2-1)))^2);
FU13 = -lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2);
FU14 = -FU11;
FU15 = -FU12;
FU16 = -FU13;
FU17 = -lmbdn*((sin(alpha3(N2-1)))^2)*(TYK(k-1,N2+N2));
FU18 = -((alpha3(N2))-(alpha3(N2-1)))*((sigma*E*(sin(alpha3(N21))))+(mu*(V^2)*(sin(alpha3(N2-1)))));
FU19 = -((alpha3(N2))-(alpha3(N2-1)))*(sin(alpha3(N2-1)));
FU20 = lmbdn*((sin(alpha3(N2-1)))^2);
F((2*(length(alpha3)))-2) = ((FU11+FU12+FU13)*(TYK(k-1,N2+N2))) +
((FU14+FU15+FU16)*(TYK(k-1,N2+N2-1))) + ((FU17+FU18)*(TYK(k-1,N2-1))) +
(FU19*((TYK(k-1,N2-1))^2)) + (FU20*(TYK(k-1,N2+N2-1))*(TYK(k-1,N2-1)));
%%%%%%%%%%%%%%%%%%%%
% RIGHT DOWN BLOCK %
%%%%%%%%%%%%%%%%%%%%
WU9 = q0*sigma*E*(cos(alpha3(1)));
WU10 = -lmbdn*sigma*E*((sin(alpha3(1)))^2);
WU11 = -lmbdn*mu*(V^2)*((sin(alpha3(1)))^2);
WU12 = -lmbdn*((sin(alpha3(1)))^2);
W(length(alpha3),(length(alpha3))+1) = (WU9+WU10+WU11) + (WU12*(TYK(k1,1)));
103
for i=((length(alpha3))+1):((2*(length(alpha3)))-3)
for j=((length(alpha3))+1):((2*(length(alpha3)))-2)
if (i == j - 1)
W15 = q0*sigma*E*cos(alpha3(j-N2));
W16 = -lmbdn*sigma*E*((sin(alpha3(j-N2)))^2);
W17 = -lmbdn*mu*(V^2)*((sin(alpha3(j-N2)))^2);
W18 = -lmbdn*((sin(alpha3(j-N2)))^2);
W(i,j) = (W15+W16+W17) + (W18*(TYK(k-1,j-N2)));
end
if (i == j)
W19 = -q0*sigma*E*cos(alpha3(j-N2+1));
W20 = lmbdn*sigma*E*((sin(alpha3(j-N2+1)))^2);
W21 = lmbdn*mu*(V^2)*((sin(alpha3(j-N2+1)))^2);
W22 = lmbdn*((sin(alpha3(j-N2+1)))^2);
W(i,j) = (W19+W20+W21) + (W22*(TYK(k-1,j-N2+1)));
end
end
end
WU13 = -q0*sigma*E*(cos(alpha3(N2-1)));
WU14 = lmbdn*sigma*E*((sin(alpha3(N2-1)))^2);
WU15 = lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2);
WU16 = lmbdn*((sin(alpha3(N2-1)))^2);
W((2*(length(alpha3)))-2,(2*(length(alpha3)))-2) = (WU13+WU14+WU15) +
(WU16*(TYK(k-1,N2-1)));
%%%%%%%%%%%%%%%%%%%%%
% MODIFIED SOLUTION %
%%%%%%%%%%%%%%%%%%%%%
DELTA = W\(-F');
for i=1:N2
TYK(k,i) = TYK(k-1,i) + DELTA(i);
end
TYK(k,N2+1) = TYK(k-1,N2+1);
for i=(N2+2):((2*N2)-1)
TYK(k,i) = TYK(k-1,i) + DELTA(i-1);
end
TYK(k,(2*N2)) = TYK(k-1,(2*N2));
end
TK = zeros(NN,N2);
YK = zeros(NN,N2);
for k=1:NN
for i=1:N2
TK(k,i) = TYK(k,i) + (mu*(V^2));
end
for i=(N2+1):(2*N2)
YK(k,i-N2) = TYK(k,i);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%
% ESTABLISHED SOLUTION %
%%%%%%%%%%%%%%%%%%%%%%%%
TE = zeros(1,N2);
YE = zeros(1,N2);
XE = zeros(1,N2);
LE = zeros(1,N2);
XE(1)=0;
LE(1)=0;
for i=1:N2
TE(i) = TYK(k,i) + (mu*(V^2));
YE(i) = TYK(k,i+N2);
if (i < N2)
numerator = (TYK(k,i))*((alpha3(i+1))(alpha3(i)))*((sigma*E)+(TYK(k,i))+(mu*(V^2)));
denominator = (q0*sigma*E*(cos(alpha3(i)))) (lmbdn*((sin(alpha3(i)))^2)*((sigma*E)+(TYK(k,i))+(mu*(V^2))));
104
XE(i+1) = XE(i) + (((cos(alpha3(i)))*numerator)/denominator);
LE(i+1) = LE(i) + (numerator/denominator);
end
end
%%%%%%%%%
% GRAPH %
%%%%%%%%%
if (draw_figures)
figure;
plot(X,Y);
grid on;
hold on;
plot(XE,YE);
grid on;
hold on;
title('ФОРМА ГИБКОЙ НИТИ');
legend('НЕРАСТЯЖИМОЙ','РАСТЯЖИМОЙ');
xlabel('X, М');
ylabel('Y, М');
figure;
plot(Y,T);
grid on;
hold on;
plot(YE,TE);
grid on;
hold on;
title('НАТЯЖЕНИЕ ГИБКОЙ НИТИ');
legend('НЕРАСТЯЖИМОЙ','РАСТЯЖИМОЙ');
xlabel('Y, М');
ylabel('НАТЯЖЕНИЕ, Н');
end
end
5. Часть программы (решатель Wave). Решатель определяет форму и
натяжения кабеля по входным физическим данным для нелинейной системы
дифференциальных уравнений с учетом растяжимости кабеля при наличии
подводного
течения
с
переменным
по
высоте
профилем
итерационным решением – с использованием метода Ньютона.
function [X,Y,T,alpha3,TK,YK,TE,YE,XE,LE] =
Wave_Curve(d,Roc,E,V,H,Umax,opposite,step1,NN,draw_figures)
%%%%%%%%%%%%%%%%
% READING DATA %
%%%%%%%%%%%%%%%%
[sigma,mu,q0,Ct,lmbdn,lmbdt,U2,U3,xin] =
Cable_Properties(d,Roc,E,V,H,Umax,opposite);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINE INITIAL APPROXIMATION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step1 = 0.00004;
[X,Y,T,alpha2] = Implicit_Curve(d,Roc,E,V,H,Umax,opposite,step1,0);
N2 = length(alpha2);
105
скорости
TYK = [T-(mu*V^2),Y];
%%%%%%%%%%%%%%%%%%%
% NEWTON SOLUTION %
%%%%%%%%%%%%%%%%%%%
% NN = 3;
W = zeros((2*N2)-2,(2*N2)-2);
F = zeros(1,(2*N2)-2);
alpha3 = alpha2;
for k=2:NN
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REDEFINE SINGULAR POINT %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
CFK = -((-lmbdn)(opposite*xin*(((U2^2)*(H^4))+(2*U2*U3*(H^5))+((U3^2)*(H^6)))))*((sigma*E)+(T
YK(k-1,length(alpha3)))+(mu*(V^2)));
AzeroK = acos((((-q0*sigma*E) + sqrt(((q0*sigma*E)^2) +
(4*(CFK^2))))/(2*CFK)));
AmoveK = AzeroK - step1;
alpha3 = 0:((AmoveK-0)/(N2-1)):AmoveK;
%%%%%%%%%%%%%%%%%%%%
% LEFT UPPER BLOCK %
%%%%%%%%%%%%%%%%%%%%
WU1 = (q0*sigma*E*cos(alpha3(1)))+(lmbdn*sigma*E*((sin(alpha3(1)))^2))+(lmbdn*mu*(V^
2)*((sin(alpha3(1)))^2))-(((alpha3(2))(alpha3(1)))*((q0*sigma*E*(sin(alpha3(1))))-(lmbdt*sigma*E*(1(cos(alpha3(1)))))-(lmbdt*mu*(V^2)*(1-(cos(alpha3(1)))))));
WU2 = opposite*((alpha3(2))(alpha3(1)))*((Ct*sigma*E*(cos(alpha3(1))))+(Ct*mu*(V^2)*(cos(alpha3(1)))));
WU3 =
opposite*((xin*((sin(alpha3(1)))^2)*sigma*E)+(xin*((sin(alpha3(1)))^2)*mu*(V^
2)));
WU4 = 2*((lmbdn*((sin(alpha3(1)))^2))-((alpha3(2))-(alpha3(1)))*(lmbdt*(1-(cos(alpha3(1))))));
WU5 = -lmbdn*((sin(alpha3(1)))^2);
WU6 = opposite*2*((alpha3(2))-(alpha3(1)))*Ct*(cos(alpha3(1)));
WU7 = opposite*2*xin*((sin(alpha3(1)))^2);
WU8 = -opposite*xin*((sin(alpha3(1)))^2);
W(1,1) = WU1+(WU2*((U2*((TYK(k-1,N2+1))^2))+(U3*((TYK(k1,N2+1))^3))))+(WU3*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))))+(WU4*(TYK(k-1,1)))+(WU5*(TYK(k1,2)))+(WU6*(TYK(k-1,1))*((U2*((TYK(k-1,N2+1))^2))+(U3*((TYK(k1,N2+1))^3))))+(WU7*(TYK(k-1,1))*(((U2^2)*((TYK(k1,N2+1))^4))+(2*U2*U3*((TYK(k-1,N2+1))^5))+((U3^2)*((TYK(k1,N2+1))^6))))+(WU8*(TYK(k-1,2))*(((U2^2)*((TYK(k1,N2+1))^4))+(2*U2*U3*((TYK(k-1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))));
WU9 = (q0*sigma*E*cos(alpha3(1)))-(lmbdn*sigma*E*((sin(alpha3(1)))^2))(lmbdn*mu*(V^2)*((sin(alpha3(1)))^2));
WU10 = -WU3;
WU11 = WU5;
WU12 = WU8;
W(1,2) = WU9+(WU10*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))))+(WU11*(TYK(k-1,1)))+(WU12*(TYK(k1,1))*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))));
FU1 = (q0*sigma*E*cos(alpha3(1)))+(lmbdn*sigma*E*((sin(alpha3(1)))^2))+(lmbdn*mu*(V^
2)*((sin(alpha3(1)))^2))-(((alpha3(2))(alpha3(1)))*((q0*sigma*E*(sin(alpha3(1))))-(lmbdt*sigma*E*(1(cos(alpha3(1)))))-(lmbdt*mu*(V^2)*(1-(cos(alpha3(1)))))));
FU2 = (q0*sigma*E*cos(alpha3(1)))-(lmbdn*sigma*E*((sin(alpha3(1)))^2))(lmbdn*mu*(V^2)*((sin(alpha3(1)))^2));
FU3 = ((lmbdn*((sin(alpha3(1)))^2))+((alpha3(2))-(alpha3(1)))*(lmbdt*(1(cos(alpha3(1))))));
106
FU4 = -lmbdn*((sin(alpha3(1)))^2);
FU5 = opposite*((alpha3(2))(alpha3(1)))*((Ct*sigma*E*(cos(alpha3(1))))+(Ct*mu*(V^2)*(cos(alpha3(1)))));
FU6 =
opposite*((xin*((sin(alpha3(1)))^2)*sigma*E)+(xin*((sin(alpha3(1)))^2)*mu*(V^
2)));
FU7 = -FU6;
FU8 = opposite*((alpha3(2))-(alpha3(1)))*Ct*(cos(alpha3(1)));
FU9 = opposite*xin*((sin(alpha3(1)))^2);
FU10 = -FU9;
F(1) = (FU1*(TYK(k-1,1)))+(FU2*(TYK(k-1,2)))+(FU3*((TYK(k1,1))^2))+(FU4*(TYK(k-1,1))*(TYK(k-1,2)))+(FU5*(TYK(k-1,1))*((U2*((TYK(k1,N2+1))^2))+(U3*((TYK(k-1,N2+1))^3))))+(FU6*(TYK(k-1,1))*(((U2^2)*((TYK(k1,N2+1))^4))+(2*U2*U3*((TYK(k-1,N2+1))^5))+((U3^2)*((TYK(k1,N2+1))^6))))+(FU7*(TYK(k-1,2))*(((U2^2)*((TYK(k1,N2+1))^4))+(2*U2*U3*((TYK(k-1,N2+1))^5))+((U3^2)*((TYK(k1,N2+1))^6))))+(FU8*((TYK(k-1,1))^2)*((U2*((TYK(k-1,N2+1))^2))+(U3*((TYK(k1,N2+1))^3))))+(FU9*((TYK(k-1,1))^2)*(((U2^2)*((TYK(k1,N2+1))^4))+(2*U2*U3*((TYK(k-1,N2+1))^5))+((U3^2)*((TYK(k1,N2+1))^6))))+(FU10*(TYK(k-1,1))*(TYK(k-1,2))*(((U2^2)*((TYK(k1,N2+1))^4))+(2*U2*U3*((TYK(k-1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))));
for i=2:((length(alpha3))-2)
for j=2:((length(alpha3))-1)
if (i == j)
W1 = q0*sigma*E*cos(alpha3(j))+(lmbdn*sigma*E*((sin(alpha3(j)))^2))+(lmbdn*mu*(V^2
)*((sin(alpha3(j)))^2))-(((alpha3(j+1))(alpha3(j)))*((q0*sigma*E*(sin(alpha3(j))))-(lmbdt*sigma*E*(1(cos(alpha3(j)))))-(lmbdt*mu*(V^2)*(1-(cos(alpha3(j)))))));
W2 = 2*((lmbdn*((sin(alpha3(j)))^2))-((alpha3(j+1))(alpha3(j)))*(-lmbdt*(1-(cos(alpha3(j))))));
W3 = -lmbdn*((sin(alpha3(j)))^2);
W4 = opposite*((alpha3(j+1))(alpha3(j)))*((Ct*sigma*E*(cos(alpha3(j))))+(Ct*mu*(V^2)*(cos(alpha3(j)))));
W5 =
opposite*((xin*((sin(alpha3(j)))^2)*sigma*E)+(xin*((sin(alpha3(j)))^2)*mu*(V^
2)));
W6 = opposite*2*((alpha3(j+1))(alpha3(j)))*Ct*(cos(alpha3(j)));
W7 = opposite*2*xin*((sin(alpha3(j)))^2);
W8 = -opposite*xin*((sin(alpha3(j)))^2);
W(i,j) = W1+(W2*TYK(k-1,j))+(W3*TYK(k1,j+1))+(W4*((U2*((TYK(k-1,j+N2))^2))+(U3*((TYK(k1,j+N2))^3))))+(W5*(((U2^2)*((TYK(k-1,j+N2))^4))+(2*U2*U3*((TYK(k1,j+N2))^5))+((U3^2)*((TYK(k-1,j+N2))^6))))+(W6*(TYK(k-1,j))*((U2*((TYK(k1,j+N2))^2))+(U3*((TYK(k-1,j+N2))^3))))+(W7*(TYK(k-1,j))*(((U2^2)*((TYK(k1,j+N2))^4))+(2*U2*U3*((TYK(k-1,j+N2))^5))+((U3^2)*((TYK(k1,j+N2))^6))))+(W8*(TYK(k-1,j+1))*(((U2^2)*((TYK(k1,j+N2))^4))+(2*U2*U3*((TYK(k-1,j+N2))^5))+((U3^2)*((TYK(k-1,j+N2))^6))));
end
if (i+1 == j)
W9 = q0*sigma*E*cos(alpha3(j))(lmbdn*sigma*E*((sin(alpha3(j)))^2))-(lmbdn*mu*(V^2)*((sin(alpha3(j)))^2));
W10 = W3;
W11 = -W5;
W12 = W8;
W(i,j) = W9+(W10*(TYK(k-1,j)))+(W11*(((U2^2)*((TYK(k1,j+N2))^4))+(2*U2*U3*((TYK(k-1,j+N2))^5))+((U3^2)*((TYK(k1,j+N2))^6))))+(W12*(TYK(k-1,j))*(((U2^2)*((TYK(k1,j+N2))^4))+(2*U2*U3*((TYK(k-1,j+N2))^5))+((U3^2)*((TYK(k-1,j+N2))^6))));
end
end
F1 = (q0*sigma*E*cos(alpha3(i)))+(lmbdn*sigma*E*((sin(alpha3(i)))^2))+(lmbdn*mu*(V^
107
2)*((sin(alpha3(i)))^2))-(((alpha3(i+1))(alpha3(i)))*((q0*sigma*E*(sin(alpha3(i))))-(lmbdt*sigma*E*(1(cos(alpha3(i)))))-(lmbdt*mu*(V^2)*(1-(cos(alpha3(i)))))));
F2 = (q0*sigma*E*cos(alpha3(i)))(lmbdn*sigma*E*((sin(alpha3(i)))^2))-(lmbdn*mu*(V^2)*((sin(alpha3(i)))^2));
F3 = ((lmbdn*((sin(alpha3(i)))^2))-((alpha3(i+1))-(alpha3(i)))*(lmbdt*(1-(cos(alpha3(i))))));
F4 = -lmbdn*((sin(alpha3(i)))^2);
F5 = opposite*((alpha3(i+1))(alpha3(i)))*((Ct*sigma*E*(cos(alpha3(i))))+(Ct*mu*(V^2)*(cos(alpha3(i)))));
F6 =
opposite*((xin*((sin(alpha3(i)))^2)*sigma*E)+(xin*((sin(alpha3(i)))^2)*mu*(V^
2)));
F7 = -F6;
F8 = opposite*((alpha3(i+1))-(alpha3(i)))*Ct*(cos(alpha3(i)));
F9 = opposite*xin*((sin(alpha3(i)))^2);
F10 = -F9;
F(i) = (F1*(TYK(k-1,i)))+(F2*(TYK(k-1,i+1)))+(F3*((TYK(k1,i))^2))+(F4*(TYK(k-1,i))*(TYK(k-1,i+1)))+(F5*(TYK(k-1,i))*((U2*((TYK(k1,N2+i))^2))+(U3*((TYK(k-1,N2+i))^3))))+(F6*(TYK(k-1,i))*(((U2^2)*((TYK(k1,N2+i))^4))+(2*U2*U3*((TYK(k-1,N2+i))^5))+((U3^2)*((TYK(k1,N2+i))^6))))+(F7*(TYK(k-1,i+1))*(((U2^2)*((TYK(k1,N2+i))^4))+(2*U2*U3*((TYK(k-1,N2+i))^5))+((U3^2)*((TYK(k1,N2+i))^6))))+(F8*((TYK(k-1,i))^2)*((U2*((TYK(k-1,N2+i))^2))+(U3*((TYK(k1,N2+i))^3))))+(F9*((TYK(k-1,i))^2)*(((U2^2)*((TYK(k1,N2+i))^4))+(2*U2*U3*((TYK(k-1,N2+i))^5))+((U3^2)*((TYK(k1,N2+i))^6))))+(F10*(TYK(k-1,i))*(TYK(k-1,i+1))*(((U2^2)*((TYK(k1,N2+i))^4))+(2*U2*U3*((TYK(k-1,N2+i))^5))+((U3^2)*((TYK(k-1,N2+i))^6))));
end
WU13 = (-q0*sigma*E*cos(alpha3(N2-1)))+(lmbdn*sigma*E*((sin(alpha3(N21)))^2))+(lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2))-(((alpha3(N2))-(alpha3(N21)))*((q0*sigma*E*(sin(alpha3(N2-1))))-(lmbdt*sigma*E*(1-(cos(alpha3(N21)))))-(lmbdt*mu*(V^2)*(1-(cos(alpha3(N2-1)))))));
WU14 = 2*((lmbdn*((sin(alpha3(N2-1)))^2))-((alpha3(N2))-(alpha3(N21)))*(-lmbdt*(1-(cos(alpha3(N2-1))))));
WU15 = -lmbdn*((sin(alpha3(N2-1)))^2);
WU16 = opposite*((alpha3(N2))-(alpha3(N21)))*((Ct*sigma*E*(cos(alpha3(N2-1))))+(Ct*mu*(V^2)*(cos(alpha3(N2-1)))));
WU17 = opposite*((xin*((sin(alpha3(N21)))^2)*sigma*E)+(xin*((sin(alpha3(N2-1)))^2)*mu*(V^2)));
WU18 = opposite*2*((alpha3(N2))-(alpha3(N2-1)))*Ct*(cos(alpha3(N2-1)));
WU19 = opposite*2*xin*((sin(alpha3(N2-1)))^2);
WU20 = -opposite*xin*((sin(alpha3(N2-1)))^2);
W(N2-1,N2-1) = WU13+(WU14*(TYK(k-1,N2-1)))+(WU15*(TYK(k1,N2)))+(WU16*((U2*((TYK(k-1,(2*N2)-1))^2))+(U3*((TYK(k-1,(2*N2)1))^3))))+(WU17*(((U2^2)*((TYK(k-1,(2*N2)-1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))))+(WU18*(TYK(k-1,N21))*((U2*((TYK(k-1,(2*N2)-1))^2))+(U3*((TYK(k-1,(2*N2)1))^3))))+(WU19*(TYK(k-1,N2-1))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)1))^6))))+(WU20*(TYK(k-1,N2))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))));
WU21 = (q0*sigma*E*cos(alpha3(N2-1)))-(lmbdn*sigma*E*((sin(alpha3(N21)))^2))-(lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2));
WU22 = WU15;
WU23 = -WU17;
WU24 = WU20;
W(N2-1,N2) = WU21+(WU22*(TYK(k-1,N2-1)))+(WU23*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)1))^6))))+(WU24*(TYK(k-1,N2-1))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))));
FU11 = (-q0*sigma*E*cos(alpha3(N2-1)))+(lmbdn*sigma*E*((sin(alpha3(N21)))^2))+(lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2))-(((alpha3(N2))-(alpha3(N2-
108
1)))*((q0*sigma*E*(sin(alpha3(N2-1))))-(lmbdt*sigma*E*(1-(cos(alpha3(N21)))))-(lmbdt*mu*(V^2)*(1-(cos(alpha3(N2-1)))))));
FU12 = (q0*sigma*E*cos(alpha3(N2-1)))-(lmbdn*sigma*E*((sin(alpha3(N21)))^2))-(lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2));
FU13 = ((lmbdn*((sin(alpha3(N2-1)))^2))-((alpha3(N2))-(alpha3(N2-1)))*(lmbdt*(1-(cos(alpha3(N2-1))))));
FU14 = -lmbdn*((sin(alpha3(N2-1)))^2);
FU15 = opposite*((alpha3(N2))-(alpha3(N21)))*((Ct*sigma*E*(cos(alpha3(N2-1))))+(Ct*mu*(V^2)*(cos(alpha3(N2-1)))));
FU16 = opposite*((xin*((sin(alpha3(N21)))^2)*sigma*E)+(xin*((sin(alpha3(N2-1)))^2)*mu*(V^2)));
FU17 = -FU16;
FU18 = opposite*((alpha3(N2))-(alpha3(N2-1)))*Ct*(cos(alpha3(1)));
FU19 = opposite*xin*((sin(alpha3(N2-1)))^2);
FU20 = -FU19;
F(N2-1) = (FU11*(TYK(k-1,N2-1)))+(FU12*(TYK(k-1,N2)))+(FU13*((TYK(k-1,N21))^2))+(FU14*(TYK(k-1,N2-1))*(TYK(k-1,N2)))+(FU15*(TYK(k-1,N21))*((U2*((TYK(k-1,(2*N2)-1))^2))+(U3*((TYK(k-1,(2*N2)1))^3))))+(FU16*(TYK(k-1,N2-1))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)1))^6))))+(FU17*(TYK(k-1,N2))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)1))^6))))+(FU18*((TYK(k-1,N2-1))^2)*((U2*((TYK(k-1,(2*N2)1))^2))+(U3*((TYK(k-1,(2*N2)-1))^3))))+(FU19*((TYK(k-1,N21))^2)*(((U2^2)*((TYK(k-1,(2*N2)-1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))))+(FU20*(TYK(k-1,N2-1))*(TYK(k1,N2))*(((U2^2)*((TYK(k-1,(2*N2)-1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))));
%%%%%%%%%%%%%%%%%%%
% LEFT DOWN BLOCK %
%%%%%%%%%%%%%%%%%%%
WU25 = -((alpha3(2))(alpha3(1)))*((sigma*E*(sin(alpha3(1))))+(mu*(V^2)*(sin(alpha3(1)))));
WU26 = lmbdn*((sin(alpha3(1)))^2);
WU27 = opposite*xin*((sin(alpha3(1)))^2);
WU28 = -2*((alpha3(2))-(alpha3(1)))*(sin(alpha3(1)));
WU29 = -WU14;
WU30 = -WU15;
W(length(alpha3),1) =
WU25+(WU26*(TYK(N2+1)))+(WU27*(TYK(N2+1))*(((U2^2)*((TYK(N2+1))^4))+(2*U2*U3*
((TYK(N2+1))^5))+((U3^2)*((TYK(N2+1))^6))))+(WU28*(TYK(1)))+(WU29*(TYK(N2+2))
)+(WU30*(TYK(N2+2))*(((U2^2)*((TYK(N2+1))^4))+(2*U2*U3*((TYK(N2+1))^5))+((U3^
2)*((TYK(N2+1))^6))));
FU21 = (q0*sigma*E*cos(alpha3(1)))+(lmbdn*sigma*E*((sin(alpha3(1)))^2))+(lmbdn*mu*(V^
2)*((sin(alpha3(1)))^2));
FU22 = -FU21;
FU23 = -((alpha3(2))(alpha3(1)))*((sigma*E*(sin(alpha3(1))))+(mu*(V^2)*(sin(alpha3(1)))));
FU24 = -((alpha3(2))-(alpha3(1)))*(sin(alpha3(1)));
FU25 = lmbdn*((sin(alpha3(1)))^2);
FU26 = -FU25;
FU27 =
opposite*((xin*((sin(alpha3(1)))^2)*sigma*E)+(xin*((sin(alpha3(1)))^2)*mu*(V^
2)));
FU28 = -FU27;
FU29 = opposite*xin*((sin(alpha3(1)))^2);
FU30 = -FU29;
F(N2) = (FU21*(TYK(k-1,N2+1)))+(FU22*(TYK(k-1,N2+2)))+(FU23*(TYK(k1,1)))+(FU24*((TYK(k-1,1))^2))+(FU25*(TYK(k-1,N2+1))*(TYK(k1,1)))+(FU26*(TYK(k-1,N2+2))*(TYK(k-1,1)))+(FU27*(TYK(k1,N2+1))*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))))+(FU28*(TYK(k1,N2+2))*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k-
109
1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))))+(FU29*(TYK(k-1,N2+1))*(TYK(k1,1))*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))))+(FU30*(TYK(k-1,N2+2))*(TYK(k1,1))*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))));
for i=((length(alpha3))+1):((2*(length(alpha3)))-3)
for j=2:(length(alpha3)-2)
if (i == (j + N2 - 1))
W13 = -((alpha3(j+1))(alpha3(j)))*((sigma*E*(sin(alpha3(j))))+(mu*(V^2)*(sin(alpha3(j)))));
W14 = -2*((alpha3(j+1))-(alpha3(j)))*(sin(alpha3(j)));
W15 = lmbdn*((sin(alpha3(j)))^2);
W16 = -W15;
W17 = opposite*xin*((sin(alpha3(j)))^2);
W18 = -W17;
W(i,j) = W13+(W14*(TYK(k-1,j)))+(W15*(TYK(k1,j+N2)))+(W16*(TYK(k-1,j+1+N2)))+(W17*(TYK(k-1,j+N2))*(((U2^2)*((TYK(k1,j+N2))^4))+(2*U2*U3*((TYK(k-1,j+N2))^5))+((U3^2)*((TYK(k1,j+N2))^6))))+(W18*(TYK(k-1,j+1+N2))*(((U2^2)*((TYK(k1,j+N2))^4))+(2*U2*U3*((TYK(k-1,j+N2))^5))+((U3^2)*((TYK(k-1,j+N2))^6))));
end
end
F11 = (-q0*sigma*E*cos(alpha3(i-(N21))))+(lmbdn*sigma*E*((sin(alpha3(i-(N21))))^2))+(lmbdn*mu*(V^2)*((sin(alpha3(i-(N2-1))))^2));
F12 = -F11;
F13 = -((alpha3(i-(N2-1)+1))-(alpha3(i-(N21))))*((sigma*E*(sin(alpha3(i-(N2-1)))))+(mu*(V^2)*(sin(alpha3(i-(N2-1))))));
F14 = -((alpha3(i-(N2-1)+1))-(alpha3(i-(N2-1))))*(sin(alpha3(i-(N21))));
F15 = lmbdn*((sin(alpha3(i-(N2-1))))^2);
F16 = -F15;
F17 = opposite*((xin*((sin(alpha3(i-(N21))))^2)*sigma*E)+(xin*((sin(alpha3(i-(N2-1))))^2)*mu*(V^2)));
F18 = -F17;
F19 = opposite*xin*((sin(alpha3(i-(N2-1))))^2);
F20 = -F19;
F(i) = (F11*(TYK(k-1,i+1)))+(F12*(TYK(k-1,i+2)))+(F13*(TYK(k-1,i-(N21))))+(F14*((TYK(k-1,i-(N2-1)))^2))+(F15*(TYK(k-1,i+1))*(TYK(k-1,i-(N21))))+(F16*(TYK(k-1,i+2))*(TYK(k-1,i-(N2-1))))+(F17*(TYK(k1,i+1))*(((U2^2)*((TYK(k-1,i+1))^4))+(2*U2*U3*((TYK(k1,i+1))^5))+((U3^2)*((TYK(k-1,i+1))^6))))+(F18*(TYK(k1,i+2))*(((U2^2)*((TYK(k-1,i+1))^4))+(2*U2*U3*((TYK(k1,i+1))^5))+((U3^2)*((TYK(k-1,i+1))^6))))+(F19*(TYK(k-1,i+1))*(TYK(k-1,i-(N21)))*(((U2^2)*((TYK(k-1,i+1))^4))+(2*U2*U3*((TYK(k1,i+1))^5))+((U3^2)*((TYK(k-1,i+1))^6))))+(F20*(TYK(k-1,i+2))*(TYK(k-1,i-(N21)))*(((U2^2)*((TYK(k-1,i+1))^4))+(2*U2*U3*((TYK(k1,i+1))^5))+((U3^2)*((TYK(k-1,i+1))^6))));
end
WU31 = -((alpha3(N2))-(alpha3(N2-1)))*((sigma*E*(sin(alpha3(N21))))+(mu*(V^2)*(sin(alpha3(N2-1)))));
WU32 = -lmbdn*((sin(alpha3(N2-1)))^2);
WU33 = -2*((alpha3(N2))-(alpha3(N2-1)))*(sin(alpha3(N2-1)));
WU34 = -WU32;
WU35 = opposite*xin*((sin(alpha3(N2-1)))^2);
WU36 = -WU35;
W((2*(length(alpha3)))-2,(length(alpha3))-1) = WU31+(WU32*(TYK(k1,2*N2)))+(WU33*(TYK(k-1,N2-1)))+(WU34*(TYK(k-1,(2*N2)-1)))+(WU35*(TYK(k1,(2*N2)-1))*(((U2^2)*((TYK(k-1,(2*N2)-1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))))+(WU36*(TYK(k1,2*N2))*(((U2^2)*((TYK(k-1,(2*N2)-1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))));
FU31 = (-q0*sigma*E*cos(alpha3(N2-1)))+(lmbdn*sigma*E*((sin(alpha3(N21)))^2))+(lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2));
110
FU32 = -FU31;
FU33 = -((alpha3(N2))-(alpha3(N2-1)))*((sigma*E*(sin(alpha3(N21))))+(mu*(V^2)*(sin(alpha3(N2-1)))));
FU34 = -((alpha3(N2))-(alpha3(N2-1)))*(sin(alpha3(N2-1)));
FU35 = lmbdn*((sin(alpha3(N2-1)))^2);
FU36 = -FU35;
FU37 = opposite*((xin*((sin(alpha3(N21)))^2)*sigma*E)+(xin*((sin(alpha3(N2-1)))^2)*mu*(V^2)));
FU38 = -FU37;
FU39 = opposite*xin*((sin(alpha3(N2-1)))^2);
FU40 = -FU39;
F((2*(length(alpha3)))-2) = (FU31*(TYK(k-1,(2*N2)-1)))+(FU32*(TYK(k1,2*N2)))+(FU33*(TYK(k-1,N2-1)))+(FU34*((TYK(k-1,N2-1))^2))+(FU35*(TYK(k1,(2*N2)-1))*(TYK(k-1,N2-1)))+(FU36*(TYK(k-1,2*N2))*(TYK(k-1,N21)))+(FU37*(TYK(k-1,(2*N2)-1))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)1))^6))))+(FU38*(TYK(k-1,2*N2))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)1))^6))))+(FU39*(TYK(k-1,(2*N2)-1))*(TYK(k-1,N2-1))*(((U2^2)*((TYK(k1,(2*N2)-1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)1))^6))))+(FU40*(TYK(k-1,2*N2))*(TYK(k-1,N2-1))*(((U2^2)*((TYK(k-1,(2*N2)1))^4))+(2*U2*U3*((TYK(k-1,(2*N2)-1))^5))+((U3^2)*((TYK(k-1,(2*N2)-1))^6))));
%%%%%%%%%%%%%%%%%%%%%
% RIGHT UPPER BLOCK %
%%%%%%%%%%%%%%%%%%%%%
for i=2:((length(alpha3))-2)
for j=((length(alpha3))+1):((2*(length(alpha3)))-2)
if (i == (j - N2 + 1))
W19 = opposite*((alpha3(i+1))(alpha3(i)))*((Ct*sigma*E*(cos(alpha3(i))))+(Ct*mu*(V^2)*(cos(alpha3(i)))));
W20 =
opposite*((xin*((sin(alpha3(i)))^2)*sigma*E)+(xin*((sin(alpha3(i)))^2)*mu*(V^
2)));
W21 = -W20;
W22 = opposite*((alpha3(i+1))(alpha3(i)))*Ct*(cos(alpha3(i)));
W23 = opposite*xin*((sin(alpha3(i)))^2);
W24 = -W23;
W(i,j) = (W19*(TYK(k-1,i))*((2*U2*(TYK(k1,i+N2)))+(3*U3*((TYK(k-1,i+N2))^2))))+(W20*(TYK(k-1,i))*((4*(U2^2)*((TYK(k1,i+N2))^3))+(10*U2*U3*((TYK(k-1,i+N2))^4))+(6*(U3^2)*((TYK(k1,i+N2))^5))))+(W21*(TYK(k-1,i+1))*((4*(U2^2)*((TYK(k1,i+N2))^3))+(10*U2*U3*((TYK(k-1,i+N2))^4))+(6*(U3^2)*((TYK(k1,i+N2))^5))))+(W22*((TYK(k-1,i))^2)*((2*U2*(TYK(k-1,i+N2)))+(3*U3*((TYK(k1,i+N2))^2))))+(W23*((TYK(k-1,i))^2)*((4*(U2^2)*((TYK(k1,i+N2))^3))+(10*U2*U3*((TYK(k-1,i+N2))^4))+(6*(U3^2)*((TYK(k1,i+N2))^5))))+(W24*(TYK(k-1,i))*(TYK(k-1,i+1))*((4*(U2^2)*((TYK(k1,i+N2))^3))+(10*U2*U3*((TYK(k-1,i+N2))^4))+(6*(U3^2)*((TYK(k-1,i+N2))^5))));
end
end
end
WU37 = opposite*((alpha3(N2))-(alpha3(N21)))*((Ct*sigma*E*(cos(alpha3(N2-1))))+(Ct*mu*(V^2)*(cos(alpha3(N2-1)))));
WU38 = opposite*((xin*((sin(alpha3(N21)))^2)*sigma*E)+(xin*((sin(alpha3(N2-1)))^2)*mu*(V^2)));
WU39 = -WU38;
WU40 = opposite*((alpha3(N2))-(alpha3(N2-1)))*Ct*(cos(alpha3(N2-1)));
WU41 = opposite*xin*((sin(alpha3(N2-1)))^2);
WU42 = -WU41;
W(N2-1,(2*N2)-2) = (WU37*(TYK(k-1,N2-1))*((2*U2*(TYK(k-1,(2*N2)1)))+(3*U3*((TYK(k-1,(2*N2)-1))^2))))+(WU38*(TYK(k-1,N21))*((4*(U2^2)*((TYK(k-1,(2*N2)-1))^3))+(10*U2*U3*((TYK(k-1,(2*N2)1))^4))+(6*(U3^2)*((TYK(k-1,(2*N2)-1))^5))))+(WU39*(TYK(k1,N2))*((4*(U2^2)*((TYK(k-1,(2*N2)-1))^3))+(10*U2*U3*((TYK(k-1,(2*N2)-
111
1))^4))+(6*(U3^2)*((TYK(k-1,(2*N2)-1))^5))))+(WU40*((TYK(k-1,N21))^2)*((2*U2*(TYK(k-1,(2*N2)-1)))+(3*U3*((TYK(k-1,(2*N2)1))^2))))+(WU41*((TYK(k-1,N2-1))^2)*((4*(U2^2)*((TYK(k-1,(2*N2)1))^3))+(10*U2*U3*((TYK(k-1,(2*N2)-1))^4))+(6*(U3^2)*((TYK(k-1,(2*N2)1))^5))))+(WU42*(TYK(k-1,N2-1))*(TYK(k-1,N2))*((4*(U2^2)*((TYK(k-1,(2*N2)1))^3))+(10*U2*U3*((TYK(k-1,(2*N2)-1))^4))+(6*(U3^2)*((TYK(k-1,(2*N2)1))^5))));
%%%%%%%%%%%%%%%%%%%%
% RIGHT DOWN BLOCK %
%%%%%%%%%%%%%%%%%%%%
WU43 = (q0*sigma*E*cos(alpha3(1)))-(lmbdn*sigma*E*((sin(alpha3(1)))^2))(lmbdn*mu*(V^2)*((sin(alpha3(1)))^2));
WU44 = opposite*((xin*((sin(alpha3(1)))^2)*sigma*E)+(xin*((sin(alpha3(1)))^2)*mu*(V^
2)));
WU45 = -lmbdn*((sin(alpha3(1)))^2);
WU46 = -opposite*xin*((sin(alpha3(1)))^2);
W(N2,N2+1) = WU43+(WU44*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))))+(WU45*(TYK(k-1,1)))+(WU46*(TYK(k1,1))*(((U2^2)*((TYK(k-1,N2+1))^4))+(2*U2*U3*((TYK(k1,N2+1))^5))+((U3^2)*((TYK(k-1,N2+1))^6))));
for i=((length(alpha3))+1):((2*(length(alpha3)))-3)
for j=((length(alpha3))+1):((2*(length(alpha3)))-2)
if (i == j - 1)
W25 = (q0*sigma*E*cos(alpha3(j-N2)))(lmbdn*sigma*E*((sin(alpha3(j-N2)))^2))-(lmbdn*mu*(V^2)*((sin(alpha3(jN2)))^2));
W26 = -lmbdn*((sin(alpha3(j-N2)))^2);
W27 = -opposite*((xin*((sin(alpha3(jN2)))^2)*sigma*E)+(xin*((sin(alpha3(j-N2)))^2)*mu*(V^2)));
W28 = -opposite*xin*((sin(alpha3(j-N2)))^2);
W(i,j) = W25+(W26*(TYK(k-1,j-N2)))+(W27*(((U2^2)*((TYK(k1,j))^4))+(2*U2*U3*((TYK(k-1,j))^5))+((U3^2)*((TYK(k-1,j))^6))))+(W28*(TYK(k1,j-N2))*(((U2^2)*((TYK(k-1,j))^4))+(2*U2*U3*((TYK(k1,j))^5))+((U3^2)*((TYK(k-1,j))^6))));
end
if (i == j)
W29 = (-q0*sigma*E*cos(alpha3(jN2+1)))+(lmbdn*sigma*E*((sin(alpha3(jN2+1)))^2))+(lmbdn*mu*(V^2)*((sin(alpha3(j-N2+1)))^2));
W30 = lmbdn*((sin(alpha3(j-N2+1)))^2);
W31 = opposite*((xin*((sin(alpha3(jN2+1)))^2)*sigma*E)+(xin*((sin(alpha3(j-N2+1)))^2)*mu*(V^2)));
W32 = -W31;
W33 = opposite*xin*((sin(alpha3(j-N2+1)))^2);
W34 = -W33;
W(i,j) = W29+(W30*(TYK(k-1,j-N2+1)))+(W31*((5*(U2^2)*((TYK(k1,j+1))^4))+(12*U2*U3*((TYK(k-1,j+1))^5))+(7*(U3^2)*((TYK(k1,j+1))^6))))+(W32*(TYK(k-1,j+2))*((4*(U2^2)*((TYK(k1,j+1))^3))+(10*U2*U3*((TYK(k-1,j+1))^4))+(6*(U3^2)*((TYK(k1,j+1))^5))))+(W33*(TYK(k-1,j-N2+1))*((5*(U2^2)*((TYK(k1,j+1))^4))+(12*U2*U3*((TYK(k-1,j+1))^5))+(7*(U3^2)*((TYK(k1,j+1))^6))))+(W34*(TYK(k-1,j+2))*(TYK(k-1,j-N2+1))*((4*(U2^2)*((TYK(k1,j+1))^3))+(10*U2*U3*((TYK(k-1,j+1))^4))+(6*(U3^2)*((TYK(k-1,j+1))^5))));
end
end
end
WU47 = (-q0*sigma*E*cos(alpha3(N2-1)))+(lmbdn*sigma*E*((sin(alpha3(N21)))^2))+(lmbdn*mu*(V^2)*((sin(alpha3(N2-1)))^2));
WU48 = lmbdn*((sin(alpha3(N2-1)))^2);
WU49 = opposite*((xin*((sin(alpha3(N21)))^2)*sigma*E)+(xin*((sin(alpha3(N2-1)))^2)*mu*(V^2)));
WU50 = -WU49;
WU51 = opposite*xin*((sin(alpha3(N2-1)))^2);
112
WU52 = -WU51;
W((2*(length(alpha3)))-2,(2*(length(alpha3)))-2) = WU47+(WU48*(TYK(k1,N2-1)))+(WU49*((5*(U2^2)*((TYK(k-1,(2*N2)-1))^4))+(12*U2*U3*((TYK(k1,(2*N2)-1))^5))+(7*(U3^2)*((TYK(k-1,(2*N2)-1))^6))))+(WU50*(TYK(k1,2*N2))*((4*(U2^2)*((TYK(k-1,(2*N2)-1))^3))+(10*U2*U3*((TYK(k-1,(2*N2)1))^4))+(6*(U3^2)*((TYK(k-1,(2*N2)-1))^5))))+(WU51*(TYK(k-1,N21))*((5*(U2^2)*((TYK(k-1,(2*N2)-1))^4))+(12*U2*U3*((TYK(k-1,(2*N2)1))^5))+(7*(U3^2)*((TYK(k-1,(2*N2)-1))^6))))+(WU52*(TYK(k-1,2*N2))*(TYK(k1,N2-1))*((4*(U2^2)*((TYK(k-1,(2*N2)-1))^3))+(10*U2*U3*((TYK(k-1,(2*N2)1))^4))+(6*(U3^2)*((TYK(k-1,(2*N2)-1))^5))));
%%%%%%%%%%%%%%%%%%%%%
% MODIFIED SOLUTION %
%%%%%%%%%%%%%%%%%%%%%
DELTA = W\(-F');
for i=1:N2
TYK(k,i) = TYK(k-1,i) + DELTA(i);
end
TYK(k,N2+1) = TYK(k-1,N2+1);
for i=(N2+2):((2*N2)-1)
TYK(k,i) = TYK(k-1,i) + DELTA(i-1);
end
TYK(k,(2*N2)) = TYK(k-1,(2*N2));
end
TK = zeros(NN,N2);
YK = zeros(NN,N2);
for k=1:NN
for i=1:N2
TK(k,i) = TYK(k,i) + (mu*(V^2));
end
for i=(N2+1):(2*N2)
YK(k,i-N2) = TYK(k,i);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%
% ESTABLISHED SOLUTION %
%%%%%%%%%%%%%%%%%%%%%%%%
TE = zeros(1,N2);
YE = zeros(1,N2);
XE = zeros(1,N2);
LE = zeros(1,N2);
XE(1)=0;
LE(1)=0;
for i=1:N2
TE(i) = TYK(k,i) + (mu*(V^2));
YE(i) = TYK(k,i+N2);
if (i < N2)
numerator = (TYK(k,i))*((alpha3(i+1))(alpha3(i)))*((sigma*E)+(TYK(k,i))+(mu*(V^2)));
denominator = (q0*sigma*E*(cos(alpha3(i)))) (lmbdn*((sin(alpha3(i)))^2)*((sigma*E)+(TYK(k,i))+(mu*(V^2)))) (opposite*xin*((sin(alpha3(i)))^2)*(((U2^2)*((TYK(k,i+N2))^4))+(2*U2*U3*((TYK
(k,i+N2))^5))+((U3^2)*((TYK(k,i+N2))^6)))*((sigma*E)+(TYK(k,i))+(mu*(V^2))));
XE(i+1) = XE(i) + (((cos(alpha3(i)))*numerator)/denominator);
LE(i+1) = LE(i) + (numerator/denominator);
end
end
%%%%%%%%%
% GRAPH %
%%%%%%%%%
if (draw_figures)
figure;
plot(X,Y);
grid on;
hold on;
113
plot(XE,YE);
grid on;
hold on;
title('ФОРМА ГИБКОЙ НИТИ');
legend('НЕРАСТЯЖИМОЙ','РАСТЯЖИМОЙ');
xlabel('X, М');
ylabel('Y, М');
figure;
plot(Y,T);
grid on;
hold on;
plot(YE,TE);
grid on;
hold on;
title('НАТЯЖЕНИЕ ГИБКОЙ НИТИ');
legend('НЕРАСТЯЖИМОЙ','РАСТЯЖИМОЙ');
xlabel('Y, М');
ylabel('НАТЯЖЕНИЕ, Н');
end
end
6. Часть программы (постпроцессор), предназначенная для построения
графиков исследуемых зависимостей. Представлены различные блоки:
обращения к базе данных; анализа сходимости; верификации; варьирования
различных физических входных данных.
function Modified_Convergence
%%%%%%%%%%%%%%%%%%
% 0 DATA STORAGE %
%%%%%%%%%%%%%%%%%%
save('database.mat');
clear;
load('database.mat');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1 EULERIAN SCHEME FOR CABLE MODEL WITHOUT TENSION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X_I1,Y_I1,T_I1,alpha2_I1] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.06,0);
[X_I2,Y_I2,T_I2,alpha2_I2] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.03,0);
[X_I3,Y_I3,T_I3,alpha2_I3] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.015,0);
[X_I4,Y_I4,T_I4,alpha2_I4] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.0076,0);
[X_I5,Y_I5,T_I5,alpha2_I5] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.0038,0);
[X_I6,Y_I6,T_I6,alpha2_I6] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.0019,0);
[X_I7,Y_I7,T_I7,alpha2_I7] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.00096,0);
[X_I8,Y_I8,T_I8,alpha2_I8] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.00048,0);
114
[X_I9,Y_I9,T_I9,alpha2_I9] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.00012,0);
[X_I10,Y_I10,T_I10,alpha2_I10] =
Implicit_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.00006,0);
N_EULERIAN =
[length(T_I1),length(T_I2),length(T_I3),length(T_I4),length(T_I5),length(T_I6
),length(T_I7),length(T_I8),length(T_I9),length(T_I10)];
T_EULERIAN =
[T_I1(length(T_I1)),T_I2(length(T_I2)),T_I3(length(T_I3)),T_I4(length(T_I4)),
T_I5(length(T_I5)),T_I6(length(T_I6)),T_I7(length(T_I7)),T_I8(length(T_I8)),T
_I9(length(T_I9)),T_I10(length(T_I10))];
figure;
plot(N_EULERIAN,T_EULERIAN/(((pi*(0.1003^2))/4)*(55001025)*9.80665*5000),'black-*','LineWidth',2);
grid on;
hold on;
xlabel('КОЛИЧЕСТВО РАЗБИЕНИЙ МЕТОДА ЭЙЛЕРА N');
ylabel('БЕЗРАЗМЕРНОЕ НАТЯЖЕНИЕ T/q0*H');
title('НАТЯЖЕНИЕ КАБЕЛЯ');
angle_I1 = (180/pi)*(alpha2_I1(2) - alpha2_I1(1));
angle_I2 = (180/pi)*(alpha2_I2(2) - alpha2_I2(1));
angle_I3 = (180/pi)*(alpha2_I3(2) - alpha2_I3(1));
angle_I4 = (180/pi)*(alpha2_I4(2) - alpha2_I4(1));
angle_I5 = (180/pi)*(alpha2_I5(2) - alpha2_I5(1));
angle_I6 = (180/pi)*(alpha2_I6(2) - alpha2_I6(1));
angle_I7 = (180/pi)*(alpha2_I7(2) - alpha2_I7(1));
angle_I8 = (180/pi)*(alpha2_I8(2) - alpha2_I8(1));
angle_I9 = (180/pi)*(alpha2_I9(2) - alpha2_I9(1));
angle_I10 = (180/pi)*(alpha2_I10(2) - alpha2_I10(1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2 NEWTON METHOD FOR CABLE MODEL WITH TENSION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X_I11,Y_I11,T_I11,alpha2_I11,TK_I11,YK_I11,TE_I11,YE_I11,XE_I11,LE_I11] =
Modified_Curve(0.041,1300,7e8,3*0.5144,5000,1,1,0.0002,51,0);
N_NEWTON = 0:1:(size(TK_I11,1)-1);
T_NEWTON = 0:1:(size(TK_I11,1)-1);
for i=1:(size(TK_I11,1))
T_NEWTON(i) = TK_I11(i,size(TK_I11,2));
end
figure;
plot(N_NEWTON,T_NEWTON/(((pi*(0.041^2))/4)*(1300-1025)*9.80665*5000),'blue*','LineWidth',2);
grid on;
hold on;
xlabel('КОЛИЧЕСТВО ИТЕРАЦИЙ МЕТОДА НЬЮТОНА k');
ylabel('БЕЗРАЗМЕРНОЕ НАТЯЖЕНИЕ T/q0*H');
title('НАТЯЖЕНИЕ КАБЕЛЯ');
angle_I11 = (180/pi)*(alpha2_I11(2) - alpha2_I11(1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3 COMPARISON OF DIFFERENT MODELS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X_I12,Y_I12,T_I12,N1_I12,C_I12] =
Merkin_Curve(0.00599,7850,2.15e14,3*0.5144,5000,1,1,0.00003,0);
[X_I13,Y_I13,T_I13,alpha2_I13] =
Implicit_Curve(0.00599,7850,2.15e14,3*0.5144,5000,1,1,0.00003,0);
[X_I14,Y_I14,T_I14,alpha2_I14,TK_I14,YK_I14,TE_I14,YE_I14,XE_I14,LE_I14] =
Modified_Curve(0.00599,7850,2.15e14,3*0.5144,5000,1,1,0.00003,6,0);
time_I12 = [61.335005, 60.009655, 59.152498, 58.452049, 58.893390];
time_AVG_I12 = 59.6;
time_I13 = [373.777274, 247.866681, 268.013729, 243.105090, 210.940058];
time_AVG_I13 = 268.7;
115
time_I14_N1 = [425.768919, 470.984755, 506.827169, 426.985569, 464.118817];
time_AVG_I14_N1 = 458.9;
time_I14_N5 = [1423.001206, 1777.255751, 1149.759327, 1344.212041,
1378.551973];
time_AVG_I14_N5 = 1414.6;
time_I14_N27 = [6326.316935, 6450.634972, 6241.156647, 6428.365258,
6387.143878];
time_AVG_I14_N27 = 6366.7;
figure;
plot(X_I12/5000,Y_I12/5000,'red:','LineWidth',2);
grid on;
hold on;
plot(X_I13/5000,Y_I13/5000,'red--','LineWidth',2);
grid on;
hold on;
plot(XE_I14/5000,YE_I14/5000,'red-','LineWidth',2);
grid on;
hold on;
xlabel('БЕЗРАЗМЕРНАЯ АБСЦИССА X/H');
ylabel('БЕЗРАЗМЕРНАЯ ОРДИНАТА Y/H');
title('ФОРМА КАБЕЛЯ');
legend('АНАЛИТИЧЕСКОЕ РЕШЕНИЕ МЕРКИНА','ЧИСЛЕННОЕ РЕШЕНИЕ ДЛЯ НЕРАСТЯЖИМОГО
КАБЕЛЯ','ЧИСЛЕННОЕ РЕШЕНИЕ ДЛЯ РАСТЯЖИМОГО КАБЕЛЯ');
figure;
plot(Y_I12/5000,T_I12/(((pi*(0.00599^2))/4)*(78501025)*9.80665*5000),'red:','LineWidth',2);
grid on;
hold on;
plot(Y_I13/5000,T_I13/(((pi*(0.00599^2))/4)*(7850-1025)*9.80665*5000),'red-','LineWidth',2);
grid on;
hold on;
plot(YE_I14/5000,TE_I14/(((pi*(0.00599^2))/4)*(78501025)*9.80665*5000),'red','LineWidth',2);
grid on;
hold on;
xlabel('БЕЗРАЗМЕРНАЯ ОРДИНАТА Y/H');
ylabel('БЕЗРАЗМЕРНОЕ НАТЯЖЕНИЕ T/q0*H');
title('НАТЯЖЕНИЕ КАБЕЛЯ');
legend('АНАЛИТИЧЕСКОЕ РЕШЕНИЕ МЕРКИНА','ЧИСЛЕННОЕ РЕШЕНИЕ ДЛЯ НЕРАСТЯЖИМОГО
КАБЕЛЯ','ЧИСЛЕННОЕ РЕШЕНИЕ ДЛЯ РАСТЯЖИМОГО КАБЕЛЯ');
angle_I13 = (180/pi)*(alpha2_I13(2) - alpha2_I13(1));
angle_I14 = (180/pi)*(alpha2_I14(2) - alpha2_I14(1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4 VELOCITY VARIATION FOR CABLE MODEL WITH TENSION %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[X_I15,Y_I15,T_I15,alpha2_I15,TK_I15,YK_I15,TE_I15,YE_I15,XE_I15,LE_I15] =
Modified_Curve(0.047,3112.5,9e9,1*0.5144,5000,0.00018,28,0);
[X_I15,Y_I15,T_I15,alpha2_I15,TK_I15,YK_I15,TE_I15,YE_I15,XE_I15,LE_I15] =
Modified_Curve(0.047,3112.5,9e9,1*0.5144,5000,1,1,0.0002,28,0);
%[X_I16,Y_I16,T_I16,alpha2_I16,TK_I16,YK_I16,TE_I16,YE_I16,XE_I16,LE_I16] =
Modified_Curve(0.047,3112.5,9e9,3*0.5144,5000,0.00018,28,0);
%[X_I17,Y_I17,T_I17,alpha2_I17,TK_I17,YK_I17,TE_I17,YE_I17,XE_I17,LE_I17] =
Modified_Curve(0.047,3112.5,9e9,5*0.5144,5000,0.00018,28,0);
[X_I17,Y_I17,T_I17,alpha2_I17,TK_I17,YK_I17,TE_I17,YE_I17,XE_I17,LE_I17] =
Modified_Curve(0.047,3112.5,9e9,5*0.5144,5000,1,1,0.0002,28,0);
angle_I15 = (180/pi)*(alpha2_I15(2) - alpha2_I15(1));
%angle_I16 = (180/pi)*(alpha2_I16(2) - alpha2_I16(1));
angle_I17 = (180/pi)*(alpha2_I17(2) - alpha2_I17(1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5 MECHANICAL PROPERTIES VARIATION FOR CABLE MODEL WITH TENSION %
116
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X_I18,Y_I18,T_I18,alpha2_I18,TK_I18,YK_I18,TE_I18,YE_I18,XE_I18,LE_I18]
Modified_Curve(0.041,1300,7e8,3*0.5144,5000,1,1,0.0002,28,0);
[X_I19,Y_I19,T_I19,alpha2_I19,TK_I19,YK_I19,TE_I19,YE_I19,XE_I19,LE_I19]
Modified_Curve(0.047,3112.5,9e9,3*0.5144,5000,1,1,0.0002,28,0);
[X_I20,Y_I20,T_I20,alpha2_I20,TK_I20,YK_I20,TE_I20,YE_I20,XE_I20,LE_I20]
Modified_Curve(0.1003,5500,2.55e10,3*0.5144,5000,1,1,0.0002,28,0);
[X_I21,Y_I21,T_I21,alpha2_I21,TK_I21,YK_I21,TE_I21,YE_I21,XE_I21,LE_I21]
Modified_Curve(0.00599,7850,2.15e11,3*0.5144,5000,1,1,0.0002,28,0);
angle_I18 = (180/pi)*(alpha2_I18(2) - alpha2_I18(1));
angle_I19 = (180/pi)*(alpha2_I19(2) - alpha2_I19(1));
angle_I20 = (180/pi)*(alpha2_I20(2) - alpha2_I20(1));
angle_I21 = (180/pi)*(alpha2_I21(2) - alpha2_I21(1));
=
=
=
=
%%%%%%%%%%%%%%%%%%%%%
% 4,5 POSTPROCESSOR %
%%%%%%%%%%%%%%%%%%%%%
figure;
plot(XE_I18/5000,YE_I18/5000,'blue-','LineWidth',2);
grid on;
hold on;
plot(XE_I15/5000,YE_I15/5000,'green:','LineWidth',2);
grid on;
hold on;
plot(XE_I19/5000,YE_I19/5000,'green-','LineWidth',2);
grid on;
hold on;
plot(XE_I17/5000,YE_I17/5000,'green-.','LineWidth',2);
grid on;
hold on;
plot(XE_I20/5000,YE_I20/5000,'black-','LineWidth',2);
grid on;
hold on;
plot(XE_I21/5000,YE_I21/5000,'red-','LineWidth',2);
grid on;
hold on;
xlabel('БЕЗРАЗМЕРНАЯ АБСЦИССА X/H');
ylabel('БЕЗРАЗМЕРНАЯ ОРДИНАТА Y/H');
title('ФОРМА КАБЕЛЯ');
legend('ТИП КАБЕЛЯ №1, 3 УЗЛА','ТИП КАБЕЛЯ №2, 1 УЗЕЛ','ТИП КАБЕЛЯ №2, 3
УЗЛА','ТИП КАБЕЛЯ №2, 5 УЗЛОВ','ТИП КАБЕЛЯ №3, 3 УЗЛА','ТИП КАБЕЛЯ №4, 3
УЗЛА');
figure;
plot(YE_I18/5000,TE_I18/(((pi*(0.041^2))/4)*(1300-1025)*9.80665*5000),'blue','LineWidth',2);
grid on;
hold on;
plot(YE_I15/5000,TE_I15/(((pi*(0.047^2))/4)*(3112.51025)*9.80665*5000),'green:','LineWidth',2);
grid on;
hold on;
plot(YE_I19/5000,TE_I19/(((pi*(0.047^2))/4)*(3112.51025)*9.80665*5000),'green-','LineWidth',2);
grid on;
hold on;
plot(YE_I17/5000,TE_I17/(((pi*(0.047^2))/4)*(3112.51025)*9.80665*5000),'green-.','LineWidth',2);
grid on;
hold on;
plot(YE_I20/5000,TE_I20/(((pi*(0.1003^2))/4)*(55001025)*9.80665*5000),'black-','LineWidth',2);
grid on;
117
hold on;
plot(YE_I21/5000,TE_I21/(((pi*(0.00599^2))/4)*(7850-1025)*9.80665*5000),'red','LineWidth',2);
grid on;
hold on;
xlabel('БЕЗРАЗМЕРНАЯ ОРДИНАТА Y/H');
ylabel('БЕЗРАЗМЕРНОЕ НАТЯЖЕНИЕ T/q0*H');
title('НАТЯЖЕНИЕ КАБЕЛЯ');
legend('ТИП КАБЕЛЯ №1, 3 УЗЛА','ТИП КАБЕЛЯ №2, 1 УЗЕЛ','ТИП КАБЕЛЯ №2, 3
УЗЛА','ТИП КАБЕЛЯ №2, 5 УЗЛОВ','ТИП КАБЕЛЯ №3, 3 УЗЛА','ТИП КАБЕЛЯ №4, 3
УЗЛА');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6 WAVE VARIATION FOR CABLE MODEL WITH TENSION AND FLOW %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X_I22,Y_I22,T_I22,alpha2_I22,TK_I22,YK_I22,TE_I22,YE_I22,XE_I22,LE_I22] =
Wave_Curve(0.00599,7850,2.15e11,3*0.5144,5000,0.24,1,0.00026,201,0);
[X_I23,Y_I23,T_I23,alpha2_I23,TK_I23,YK_I23,TE_I23,YE_I23,XE_I23,LE_I23] =
Modified_Curve(0.00599,7850,2.15e11,3*0.5144,5000,0.24,1,0.00026,201,0);
H_profile = 1;
Umax_profile = 0.24;
U2_profile = (3*Umax_profile)/(H_profile^2);
U3_profile = -(2*Umax_profile)/(H_profile^3);
y_profile = 0:0.01:H_profile;
for i=1:length(y_profile)
V_profile(i) = (U3_profile*((y_profile(i))^3)) +
(U2_profile*((y_profile(i))^2));
end
figure;
plot(XE_I23/5000,YE_I23/5000,'red-','LineWidth',2);
grid on;
hold on;
plot(XE_I22/5000,YE_I22/5000,'cyan-','LineWidth',2);
grid on;
hold on;
plot(3.5-2.5*V_profile,y_profile,'cyan--','LineWidth',2);
grid on;
hold on;
xlabel('БЕЗРАЗМЕРНАЯ АБСЦИССА X/H');
ylabel('БЕЗРАЗМЕРНАЯ ОРДИНАТА Y/H');
title('ФОРМА КАБЕЛЯ');
legend('ЧИСЛЕННОЕ РЕШЕНИЕ ДЛЯ РАСТЯЖИМОГО КАБЕЛЯ','ЧИСЛЕННОЕ РЕШЕНИЕ ДЛЯ
РАСТЯЖИМОГО КАБЕЛЯ ПРИ ВСТРЕЧНОМ ТЕЧЕНИИ','ПРОФИЛЬ ПОДВОДНОГО ТЕЧЕНИЯ');
angle_I22 = (180/pi)*(alpha2_I22(2) - alpha2_I22(1));
angle_I23 = (180/pi)*(alpha2_I23(2) - alpha2_I23(1));
end
118
Приложение Б. Интегрированная программа по определению
нестационарного пространственного нагружения кабеля
Программа, предназначенная для приложения пользовательской нагрузки к
элементу кабеля в каждый момент времени в соответствие с текущими
данными о движении и ориентации элемента в пространстве.
subroutine vdload (
1 nblock, ndim, stepTime, totalTime,
2 amplitude, curCoords, velocity, dirCos, jltyp, sname,
3 value )
include 'vaba_param.inc'
dimension curCoords(nblock,ndim), velocity(nblock,ndim),
1 dirCos(nblock,ndim,ndim), value(nblock)
character*80 sname
REAL*8 x_coor, y_coor, z_coor
C ndim=3 --> for 3D beam elements
C nblock --> number of points to apply load
C jltyp load
C
=21 P1NU - force per unit length (FL^-1) in beam local N1-direction
C
(2 axis C in local CS )
C
=22 P2NU - force per unit length (FL^-1) in beam local N2-direction
C
(3 axis C in local CS )
C
=41 PXNU - force per unit length (FL^-1) in global X-direction
C
=42 PYNU - force per unit length (FL^-1) in global Y-direction
C
=43 PZNU - force per unit length (FL^-1) in global Z-direction
C===============
C INITIAL DATA C
C===============
parameter(
* length = 10.0,
* bottom = -1000.0,
* start = 550.0,
* dt_smooth = 100.d0,
* eps_h = -1.0*length,
* water_coord = -20.d0,
* cp = 1.1,
* ro = 1.e3,
* diam = 20.e-3,
* zero = 0.d0,
* one = 1.d0)
C==================
C USER SUBROUTINE C
C==================
x_coor = 0.0
y_coor = 0.0
z_coor = 0.0
vel_x = 0.0
vel_y = 0.0
vel_z = 0.0
dir_x = 0.0
dir_y = 0.0
dir_z = 0.0
zero_force = bottom - water_coord - eps_h
do km = 1, nblock
value(km) = zero
if (totalTime .ge. start) then
x_coor = curCoords(km,1)
119
y_coor = curCoords(km,2)
z_coor = curCoords(km,3)
vel_x = velocity(km,1)
vel_y = velocity(km,2)
vel_z = velocity(km,3)
dir_x = dirCos(km,2,1)
dir_y = dirCos(km,2,2)
dir_z = dirCos(km,2,3)
vel_n = (vel_x*dir_x + vel_y*dir_y + vel_z*dir_z)
if (y_coor .lt. water_coord) then
if (y_coor .gt. zero_force) then
if (totalTime .le. start + dt_smooth) then
value(km) = (cp*ro*diam/2)*(vel_n*vel_n)*dir_y*dir_y
* *(totalTime - start)/dt_smooth/100.0
end if
C==============================
C APPLYING FORCE STEP-BY-STEP C
C==============================
if (totalTime .gt. start + dt_smooth*2) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/90.0
end if
if (totalTime .gt. start + dt_smooth*3) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/80.0
end if
if (totalTime .gt. start + dt_smooth*4) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/70.0
end if
if (totalTime .gt. start + dt_smooth*5) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/60.0
end if
if (totalTime .gt. start + dt_smooth*6) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/50.0
end if
if (totalTime .gt. start + dt_smooth*7) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/40.0
end if
if (totalTime .gt. start + dt_smooth*8) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/30.0
end if
if (totalTime .gt. start + dt_smooth*9) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y/10.0
end if
C================================================
C END OF STABILIZATION TIME TO APPLY 100% FORCE C
C================================================
if (totalTime .gt. start + dt_smooth*10) then
value(km) = cp*ro*diam/2*(vel_n*vel_n)*dir_y*dir_y
end if
end if
end if
end if
end do
return
end
120
Отзывы:
Авторизуйтесь, чтобы оставить отзыв