Определение НДС цилиндрической оболочки по заданным перемещениям

Еременко Владимир Романович

Определение НДС цилиндрической оболочки по заданным перемещениям

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

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d363b5f1be77c40d58a0e

UUID: f9bddba7-4545-4c6c-baa9-553e07f23e25

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 259

Еременко Владимир Романович

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 1018087 bytes

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ КАФЕДРА ВЫЧИСЛИТЕЛЬНЫХ МЕТОДОВ МЕХАНИКИ ДЕФОРМИРУЕМОГО ТЕЛА Еременко Владимир Романович Выпускная квалификационная работа бакалавра Определение напряженно деформированного состояния цилиндрической оболочки по заданным перемещениям Направление 010400 Прикладная математика и информатика Научный руководитель, кандидат физ.-мат. наук, доцент Кабриц С.А. Санкт-Петербург 2016

Содержание Введение ........................................................................................................ 3 Постановка задачи ....................................................................................... 5 Обзор литературы ........................................................................................ 6 Глава 1. Математическая формулировка задачи .................................... 10 1.1. Уравнения цилиндрической оболочки ....................................... 11 1.1. Решение поставленной задачи ..................................................... 16 Глава 2. Метод численного решения ....................................................... 19 2.1. Метод ортогональной прогонки .................................................. 21 Глава 3. Анализ полученных результатов ............................................... 23 3.1. Алгоритм нахождения значений нагрузок .................................. 23 3.2. Графики напряжений ..................................................................... 26 3.3. Графики напряжений на лицевых поверхностях трубопровода 29 Выводы ........................................................................................................ 32 Заключение ................................................................................................. 33 Список литературы .................................................................................... 34 Приложение ................................................................................................ 36 2

Введение Нефтегазовая отрасль играет ведущую роль в современной экономике, но она сопряжена со многими проблемами, связанными с транспортировкой продукции на дальние расстояния, как по суше, так и морским путем. Известно много случаев разлития нефти в море, или прорывов трубопроводов на суше. Следствием подобных аварий являются большие человеческие жертвы, а также финансовые потери, вызванные разрушением дорогостоящего оборудования, прекращением производства и снабжением продукцией на неопределенные сроки. Роль трубопроводного транспорта в системе нефтегазовой отрасли очень высока, поскольку он транспортировки нефти на является одним из недорогих заводы и в другие регионы. видов Помимо энергетического обеспечения страны, трубопроводы позволяют уменьшить нагрузку на железнодорожный транспорт. Магистральные трубопроводы состоят из нескольких основных элементов, функционирование каждого из которых необходимо для продолжения транспортировки. Однако анализ надежности и систематизация произошедших аварий показывают, что наибольшее влияние на надежность конструкции оказывают её прямолинейные участки. Трубопроводы зачастую пролегают в труднодоступных не населенных районах, что ведет к увеличению возможности отказа. Помимо ущерба, наносимого производству, аварии могут привести к ущербу окружающей среды. Особенно опасно если трубопровод проходит вблизи крупного населенного пункта. Именно поэтому необходимо особое внимание уделять мониторингу технических средств и модернизации способов наблюдения и выявления опасных участков. Во избежание подобных проблем важным становится повышение надежности линейной части на всех этапах создания трубопровода. Также 3

важно адекватно оценить влияние внешних нагрузок на конструкции, то есть необходимо исследовать конструкцию на надежность при помощи анализа напряженно-деформированного состояния (НДС) трубопровода. Расчет НДС предполагает знание всех нагрузок, действующих на трубопровод. Эти нагрузки меняются в процессе многолетней эксплуатации трубопровода. Причинами изменения нагрузок являются подвижки грунта и т.п. Одним из способов контроля деформации трубопровода является периодическая геодезическая съемка, которая позволяет определить смещения контрольных точек по сравнению с предыдущим положением. Но при расчете НДС возникает ряд проблем таких как: 1. Данные получаемые в ходе геодезической сьемки не всегда верны, возможна ошибка в условии поставленной задачи. 2. Данные снимаются только в дискретном наборе точек, и не всегда по ним можно представить полную картину. В частности, если известны перемещения в трех точках, то можно рассмотреть результат воздействия силы с нескольких сторон. Возможно в какой-то точке из-под земли подняло камень, и он в одной точке действует на трубопровод, тогда необходимо считать, что в этой точке действует сосредоточенная сила. Но более вероятно, что в каком-то месте поднялась земля, а в каком-то наоборот вымыло её, и тогда уже следует считать, что сила распределена по некоему закону. Проблема заключается в том, что от выбора способа нагружения зависит распределение напряжений в трубопроводе. В одном случае это будет оценка, сделанная с большим запасом, поскольку будут получены значения напряжений в разы превышающий случай распределенной нагрузки. 3. Как известно [17,18] задача расчета тонких цилиндрических оболочек является плохо обусловленной (краевой эффект), поэтому для решения этих задач необходимо использовать более «тяжеловесные» численные методы. 4

Постановка задачи Необходимо определить изменение НДС фрагмента трубопровода по вертикальным смещениям, полученным посредством геодезической съемки. Тонкостенность рассматриваемого математической модели объекта использовать позволяет линейную в теорию качестве тонких цилиндрических оболочек. Обычно решение задач для цилиндрических оболочек ищут в виде тригонометрических рядов по круговой координате [3,9]. Специфика "вертикального" нагружения позволяет рассматривать только одну гармонику (случай "ветровой нагрузки"). Классическая постановка задач теории упругости и ее упрощенных вариантов - это найти НДС объекта под действием заданных нагрузок. В рассматриваемом случае задача является в некотором смысле обратной - нужно найти нагрузки, которые реализуют, полученные в эксперименте смещения. Решение такой задачи неоднозначно, т.к. одни и те же смещения в дискретном наборе точек могут возникать при различных вариантах нагружения и, соответственно, неоднозначность возникает при определении НДС. Целью работы является получение диапазона возникающих напряжений и проведение оценку прочности элементов трубопровода, как по «худшему», так и по «лучшему» развитию сценария. 5

Обзор литературы В 1822 году все основные уравнения упругости вывел Луи Коши [6]. Для вывода данных уравнений Коши использует Эйлеровы законы динамики и принцип затвердевания. В данных уравнениях полагается, что тело трехмерное и на каждую точку действует один силовой тензор. Также Коши утверждает, что существует линейная связь между тензорами напряжений и деформаций [6]: 𝜏 = 𝐶 ∙∙ 𝜖 Где 𝜖 – тензор деформаций, 𝜏 – тензор напряжений C называется тензором упругости, который имеет 36 независимых компонент в общем случае. Также существенный вклад в данную теорию внес Джордж Грин, опубликовав работу “О законах отражения и преломления света на общей поверхности двух некристаллических сред” в 1839 году. Основная его идея была в том, что необходимо выбрать в качестве основы создания уравнений теории упругости общий физический принцип, и поэтому предложил принцип: «если все действующие внутренние силы помножить на элементы их соответствующих направлений, общая сумма для любой определенной части тела всегда будет полным дифференциалом некоторой функции». То есть этот принцип описывается уравнением: 𝑑𝑈 𝑑2𝑈 𝑑𝑈 = 𝜏 ∙∙ 𝑑𝜖, ⇒ τ= ⇒ C= 𝑑𝜖 𝑑𝜖𝑑𝜖 где U— энергия деформации. Тогда число независимых компонент можно сократить до 21. А позже Новожилов доказал, что можно ограничиться 13 компонентами [16]. Позднее начала интенсивно развиваться теория оболочек, пластин и стержней. Уравнения данных теорий являются следствиями общей теории упругости. Впервые Коши и Пуассон применили разложения по степеням толщинной координаты и ограничились только низшими членами. Но данный 6

метод сразу подвергся критике Сен-Венана из-за неясности его сходимости. Однако за неимением строгих доказательств он активно используется и по сей день. По теории пластин первая работа была опубликована в 1850 году Кирхгофом. Через 9 лет вышла его же работа, но по теории стержней [10]. На данный момент методы Кирхгофа активно используются и применяются. В 1862 году Клебш применил данные методы в своей работе “Теория упругости”. В данной работе впервые рассматривались вместо напряжений усилия и моменты. Но, увы, данные предположения подверглись критике со стороны А Лява. Клебш решил пренебречь малыми перерезывающими усилиями, но Ляв заметил, что данное предположение ошибочно, так как эти усилия входят в уравнения моментов, которые также малы. В 1874 была написана первая работа по теории оболочек Г. Ароном, в которой он допустил определенные неточности. Поэтому основные идеи были предложены Лявом, Бэссетом и Лэмбой [6]. Хотя две последние вытекают из работ Лява [11]. В итоге теория Лява стала очень популярной. Однако она имела несколько недостатков [6]:  «непоследовательное обращение с малыми членами» Новожилов В.В. [17]  Неясная область применимости  Нарушение 6-го уравнения равновесия  Использование необязательного ограничения 𝜖~ℎ𝜘 Параллельно с развитием теории оболочек активно формировались методы решения краевых задач данной теории. Основные достижения в данной области принадлежат российским ученым, так как зарубежные работы включали в себя точные решения краевых задач, что является явным противоречием с эвристическим характером изначальных уравнений. Впервые, данное предположение ввел Штаерман в 1924 году, используя асимптотические методы. Поэтапно это рассматривается в книге Лурье [14] при построении решений и в работах В.В. Новожилова [17] при комплексном 7

преобразовании уравнений теории оболочек. Также асимптотические методы получили развитие в работах Гольдейнвейзера с 1939 года. К.Ф. Черных в своей работе рассмотрел комбинацию асимптотических методов и методов комплексных преобразований. Иной подход был представлен в работах Власова [4], основанный на интуитивных и физических представлениях. Одна из самых значительных работ в теории оболочек вышла в 1940 году [12]. Она принадлежала Лурье и в ней приводилось тензорное представление основных уравнений теории оболочек. В 1943 году, вышла работа Лурье [13], в которой учтены и исправлены все недостатки теории Лява, кроме области применимости. Также в 1943 году вышли многозначительные работы В.В. Новожилова [17,15]. В них он нашел решение задачи о формулировке простейшей теории оболочек, в которой было выполнено шестое уравнение равновесия. Одновременно с ним данную задачу решил Балабух [1], но возымела успех именно работа Новожилова, так как в ней не только сформулирована задача, но и дано ее использование. Позднее установили, что не существует тензорных уравнений, совпадающих с соотношениями Новожилова, без нарушения одного из свойств:  Сохранения изотропности Если уравнения в линиях кривизны записаны для изотропной оболочки, тогда в любой системе координат также должна быть изотропность  Непрерывной зависимости от радиусов кривизны Здесь говорится не только о тензорном представлении уравнений, но и о сложной операции. Эта операция была проделана в книге Черныха ([18], с. 101, уравнение (6.6)). Но на самом деле ранее подобная задача была решена в публикациях Сандерса и Койтера и получила название теории КотейраСандерса. Но эта теория является спорной, так как все вышеупомянутые работы по своей сути повторяют рассуждения Новожилова. Однако в своей работе Черных показал различие соотношений Койтера-Сандерса и 8

Новожилова. Заключительный этап в разработке теории был сделан Гольденвейзером. Он предложил иную формулировку статистических и кинематических гипотез на основе асимптотического анализа уравнений теории упругости, которые отличаются от гипотез Кирхгофа-Лява. Он принял во внимание поперечную сжимаемость оболочки и за счёт этого получил новые уравнения теории упругости, которые содержали слагаемые, потерянные в работе [15]. В настоящей работе рассматривается следующая проблема: предоставлены данные о вертикальных смещениях некоторых точек прямолинейного участка трубопровода. Требуется определить изменение НДС. Данная проблема не является новой, поскольку расчет НДС является одним из важнейших способов оценки прочности конструкции. Так, например, в [7] рассматривается подобная задача, но численная оценка производится, наиболее распространенным в настоящее время способом, методом конечных элементов (МКЭ). Трубопровод можно рассматривать как тонкую цилиндрическую оболочку, и соответственно для решения использовать уравнения цилиндрической оболочки. Итоговая система уравнений, использующаяся в работе, представлена в статьях [8,9]. Причем в [9] получены результаты под действием сосредоточенных нагрузок. Для решения поставленной задачи будет использоваться численный метод, ортогональной прогонки Годунова [5]. 9

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

1.1. Уравнения цилиндрической оболочки Нагрузка, действующая на срединной поверхности цилиндрической оболочки, может быть описана вектором 𝒒 = {−𝑞1 𝒌 + 𝑞2 𝒆𝟐 + 𝑞𝑛 𝒆𝒓 } Для решения поставленной задачи необходимо найти напряжения на внешней и внутренней лицевых поверхностях оболочки. Формулы для их поиска описаны в [8]. Компоненты напряжения усилия 𝑇1 , 𝑇2 , 𝑆, и моменты 𝑀1 , 𝑀2 , 𝐻, и компоненты перемещения 𝑢𝑟 , 𝑢2 , 𝑢𝑧 , 𝜗1 , 𝜗2 , можно представить через тригонометрические соотношения {𝑇1 , 𝑇2 } = {𝑇1,𝑘 (𝑠), 𝑇2,𝑘 (𝑠)} cos 𝑘𝜑 {𝑀1 , 𝑀2 } = {𝑀1,𝑘 (𝑠), 𝑀2,𝑘 (𝑠)} cos 𝑘𝜑 {𝑢𝑟 , 𝑢𝑧 , 𝜗1 } = {𝑢𝑟,𝑘 (𝑠), 𝑢𝑧,𝑘 (𝑠), 𝜗1,𝑘 (𝑠)} cos 𝑘𝜑 {𝑞1 , 𝑞𝑛 } = {𝑞1,𝑘 (𝑠), 𝑞𝑛,𝑘 (𝑠)} cos 𝑘𝜑 (1.1) {𝑆, 𝐻} = {𝑆,𝑘 (𝑠), 𝐻,𝑘 (𝑠)} sin 𝑘𝜑 {𝑢2 , 𝜗2 } = {𝑢2,𝑘 (𝑠), 𝜗2,𝑘 (𝑠)} sin 𝑘𝜑 {𝑞2 } = {𝑞2,𝑘 (𝑠)} sin 𝑘𝜑 Для определения коэффициентов соотношения (1.1) используется система дифференциальных уравнений, полученная в [8] и представленная в [9]. 𝑑𝐹,𝑘 = 𝑇2,𝑘 (1 − 𝑘 2 ) + 𝑅(𝑞2𝑘 − 𝑞𝑛𝑘 ) 𝑑𝑠 𝑑𝑀,𝑘 = 𝐹,𝑘 + 𝑅2 𝑞1𝑘 𝑑𝑠 𝑑𝑄12,𝑘 1 = 𝑘 (𝑇2,𝑘 + 𝑀2,𝑘 ) − R𝑞2,𝑘 𝑑𝑠 𝑅 𝑑𝑀1,𝑘 2 1 = 𝑘 (𝑄12,𝑘 − 𝐻,𝑘 ) + 𝐹,𝑘 𝑑𝑠 𝑅 𝑅 𝑑𝜘𝑡𝜈,𝑘 1 = 𝑘 𝜖1,𝑘 − 𝑘𝜘1,𝑘 𝑑𝑠 𝑅 𝑑𝑢𝑟,𝑘 = −𝜗1,𝑘 𝑑𝑠 (1.2) 11

𝑑𝜗1,𝑘 = 𝜘1,𝑘 𝑑𝑠 𝑑𝜖2,𝑘 1 1 = 𝑘𝜘𝑡𝜈,𝑘 + 𝑘 𝜔,𝑘 + (𝑘 2 − 1) 𝜗1,𝑘 𝑑𝑠 𝑅 𝑅 Неизвестные функции здесь {𝐹,𝑘 , 𝑀,𝑘 , 𝑄12,𝑘 , 𝑀1,𝑘 , 𝜘𝑡𝜈,𝑘 , 𝑢𝑟,𝑘 , 𝜗1,𝑘 , 𝜖2,𝑘 }. (1.3) 𝐹,𝑘 – усилие 𝑀,𝑘 – статический момент 𝑄12,𝑘 – обобщенная прерывающая сила 𝜘𝑡𝜈,𝑘 – компоненты изгибной деформации 𝑢𝑟,𝑘 – вектор перемещения 𝜗1,𝑘 – угол между осью вращения и нормали к срединной поверхности 𝜖2,𝑘 – компоненты тангенциальной деформации Кроме этого используются дополнительные функции {𝑇2,𝑘 , 𝑀2,𝑘 , 𝐻,𝑘 , 𝜘𝑡𝜈,𝑘 , 𝜖1,𝑘 , 𝜔,𝑘 } (1.4) (1.4), могут быть выражены через (1.3), используя закон Гука: 𝑇2,𝑘 = 𝜈𝑇1,𝑘 + 𝐵(1 − 𝜈 2 )𝜖2,𝑘 1 𝑇 − 𝜈𝜖2,𝑘 𝐵 1,𝑘 1−𝜈 𝑆,𝑘 = 𝐵 𝜔,𝑘 2 𝜖1,𝑘 = 𝑀2,𝑘 = 𝜈𝑀1,𝑘 + 𝐷(1 − 𝜈 2 )𝜘2,𝑘 𝜘1,𝑘 = (1.5) 1 𝑀 − 𝜈𝜘2,𝑘 𝐷 1,𝑘 𝐻,𝑘 = 𝐷(1 − 𝜈)𝜏,𝑘 𝐸ℎ 𝐸ℎ3 𝐵= ,𝐷 = 1 − 𝜈2 12(1 − 𝜈 2 ) А также уравнения соотношения напряженности и геометрические соотношения 𝜘2,𝑘 = 𝜖2,𝑘 𝑟 (1.6) 12

𝜘𝑡𝜈,𝑘 = 𝜏,𝑘 − 𝜔,𝑘 𝑟 В (1.5) E – это модуль Юнга, ν – это коэффициент Пуассона, h – толщина оболочки Рассматриваем случай «ветровой нагрузки», в этом случае в соотношении (1.2) можно оставить только члены с номером k=1. Так же заменой переменных (1.4), через основные переменные (1.3), с использованием (1.5), (1.6), получается система (1.7) 𝑑𝐹𝑥 = 𝜋𝑅(𝑞21 − 𝑞𝑛1 ) 𝑑𝑠 𝑑𝑀𝑦 = 𝐹𝑥 + 𝜋𝑅2 𝑞11 𝑑𝑠 𝑑𝑄12,1 = −𝛼1 𝜖2,1 − 𝛼2 𝑀𝑦 − 𝑅𝑞2,1 𝑑𝑠 𝑑𝑀1,1 = 𝛼3 𝐹𝑥 + 𝛼4 𝑄12,1 + 𝛼5 𝜘𝑡𝜈,1 𝑑𝑠 𝑑𝜘𝑡𝜈,1 = 𝛼6 𝑀1,1 + 𝛼7 𝑀𝑦 𝑑𝑠 𝑑𝑢𝑟,1 = −𝜗1,1 𝑑𝑠 𝑑𝜗1,1 𝑀1,1 𝜈𝜖2,1 = − 𝑑𝑠 𝐷 𝑅 𝑑𝜖2,1 = 𝛼8 𝜘𝑡𝜈,1 + 𝛼9 𝑄12,1 𝑑𝑠 (1.7) Где, 𝜈 (𝜋(𝐵𝑅2 + 𝐷)(𝜈 2 − 1)) 𝛼1 = ; 𝛼 = − 2 𝜋𝑅2 𝜋𝑅2 4𝐷 + 𝐵𝑅2 𝜋𝐵𝑅 2𝜋𝐵𝐷(𝜈 − 1) 𝛼3 = ; 𝛼 = ; 𝛼 = 4 5 4𝐷𝑅 + 𝐵𝜋𝑅3 4𝐷 + 𝐵𝜋𝑅2 4𝐷 + 𝐵𝜋𝑅2 𝐷 − 𝐵𝑅2 1 𝛼6 = ; 𝛼 = − 7 𝐷𝐵𝑅2 𝐵𝜋𝑅3 𝐵𝑅 2 𝛼8 = ; 𝛼 = − (𝜈 − 1)(4𝐷 + 𝐵𝑅2 ) 4𝐷 + 𝐵𝑅2 9 𝐹𝑥 = 𝜋𝐹,1 , 𝑀𝑦 = 𝜋𝑀,1 13

Соотношения (1.1) - (1.6) – служат для построения вычислительной модели, позволяющей определить напряжения в цилиндрической, главной, части трубопровода, под действием сил, связанных с поднятием или оседанием грунта. Известны перемещения поперечного сечения. В работе рассматриваются три возможных случая приложения силы. Сосредоточенные нагрузки, и распределенные (случай линейной интерполяции, и полинома Лагранжа) Известно, что перемещения, ортогональны оси трубопровода и имеют направление 𝒊𝑥 = 𝒆𝑟 cos 𝜑 − 𝒆2 sin 𝜑 Кроме того, перемещения можно записать в виде 𝒖 = 𝑢1 𝒆1 + 𝑢2 𝒆2 + 𝑤𝒏 = = 𝑢𝑟 𝒆𝑟 − 𝑢𝑧 𝒌 + 𝑢2 𝒆2 = = 𝑢𝑟,𝑘 cos 𝑘𝜑 𝒆𝑟 − 𝑢𝑧,𝑘 cos 𝑘𝜑 𝒌 + 𝑢2,𝑘 sin 𝑘𝜑 𝒆2 Тогда перемещение в координате x 𝑢𝑥 = (𝒖 ∙ 𝒊𝑥 ) = (𝑢𝑟,𝑘 cos 2 𝑘𝜑 − 𝑢2,𝑘 sin2 𝑘𝜑) Вычисляем усредненное перемещение поперечного сечения 2𝜋 𝑢𝑥средн 1 1 = ∫ 𝑢𝑥 𝑑𝜑 = (𝑢𝑟,𝑘 − 𝑢2,𝑘 ) 2𝜋 2 0 Используя соотношение из [8] 𝑢𝑟,𝑘 + 𝑢2,𝑘 = 𝜖2,𝑘 ; 𝑢2,𝑘 = 𝑟𝜖2,𝑘 − 𝑢𝑟,𝑘 𝑟 Перепишем 𝑢𝑥средн 𝑢𝑥средн = 𝑢𝑟,𝑘 − 𝑟 Будем отождествлять перемещения 𝜖2,𝑘 2 𝑢𝑥средн с полученными из эксперимента (геодезическая съемка) Сосредоточенная нагрузка в i-ом узле может быть записана, например, так: 𝑞𝑖 = 𝑝𝑖 (𝑒𝑟 cos 𝜑 − 𝑒2 sin 𝜑)𝛿(𝑠 − 𝑠𝑖 ), 𝑝𝑖 = 𝑐𝑜𝑛𝑠𝑡 14

Здесь 𝛿 – дельта-функция Дирака. А распределенная может быть записана в виде 𝑞 = 𝑝(𝑠)(𝑒𝑟 cos 𝜑 − 𝑒2 sin 𝜑) Где 𝑝(𝑠) – определяется из способа интерполяции 15

1.2. Решение поставленной задачи Формулирование краевой задачи: Введем обозначения 𝑽 = (𝑉1 , 𝑉2 , 𝑉3 , 𝑉4 , 𝑉5 , 𝑉6 , 𝑉7 , 𝑉8 ), где 𝑉1 = 𝐹𝑥 𝑉2 = 𝑀𝑦 𝑉3 = 𝑄12,1 𝑉4 = 𝑀1,1 𝑉5 = 𝜘𝑡𝜈,1 𝑉6 = 𝑢𝑟,1 𝑉7 = 𝜗1,1 𝑉8 = 𝜖2,1 Систему дифференциальных уравнений (1.7) можно представить в виде: 𝑑𝑽(𝑠) = 𝐴𝑽(𝑠) + 𝑭𝑝𝐿(𝑠), 0 ≤ 𝑠 ≤ 𝑠𝑙 𝑑𝑠 (1.8) Где A – постоянная матрица 𝐴= 0 1 0 𝛼3 0 0 0 0 −𝛼2 0 𝛼7 0 0 0 0 𝛼4 0 0 0 0 0 (0 0 𝛼9 0 0 0 0 𝛼6 0 1 𝐷 0 0 0 0 𝛼5 0 0 0 0 0 0 0 0 0 0 −𝛼1 0 0 0 0 0 0 0 −1 0 𝜈 0 0 − 𝑅 0 0 0 ) 0 𝛼8 С граничными условиями 𝐺1 𝑽(𝑠) = 0, 𝐺2 𝑽(𝑠𝑙 ) = 0 (1.9) 𝐺𝑖 , 𝑖 = 1,2 – прямоугольные матрицы следующего вида 0 0 𝐺1 = 𝐺2 = ( 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 ) 0 1 В (1.8) функция 𝐿𝑖 (𝑠) – меняется в зависимости, от способа приложения силы. 16

Рассматриваются различные варианты нагружения, обеспечивающие одни и те же перемещения (полученные из замеров). В случае сосредоточенной нагрузки 𝐿𝑖 (𝑠) = 𝛿(𝑠 − 𝑠𝑖 ) Для распределенной нагрузки рассмотрены несколько вариантов: a) Кусочно-постоянная аппроксимация, около замеров 1, |𝑠 − 𝑠𝑖 | ≤ 0.05 ∙ 𝑠𝑙 𝐿𝑖 (𝑠) = { 0, иначе b) Линейный сплайн 𝑠 − 𝑠𝑖 , 𝑠 ∈ [𝑠𝑖−1 , 𝑠𝑖 ] 𝑠𝑖 − 𝑠𝑖−1 𝐿𝑖 (𝑠) = { 𝑠 − 𝑠 𝑖+1 , 𝑠 ∈ [𝑠𝑖 , 𝑠𝑖+1 ] 𝑠𝑖+1 − 𝑠𝑖 c) Интерполяционный полином Лагранжа 𝑚 𝐿𝑖 (𝑠) = ∏ 𝑗=0,𝑗≠𝑖 𝑠 − 𝑠𝑗 𝑠𝑖 − 𝑠𝑗 В системе (1.8), 𝑽(𝑠) – вектор основных неизвестных функций, A – постоянная матрица коэффициентов, F постоянный вектор размерности 8, 𝑝𝑖 – значения нагрузки, действующие на поперечное сечение трубопровода в точке 𝑠 = 𝑠𝑖 . 𝐺1 и 𝐺2 матрицы граничных условий, размерности 4×8. Алгоритм решения На первом этапе граничная задача (1.8) – (1.9), решается с коэффициентами 𝑝𝑖 = 1, 𝑝𝑗 = 0, 𝑗 ≠ 𝑖. Из коэффициентов, формируется квадратная матрица U (матрица податливости). Компоненты матрицы, это усредненные перемещения, 𝑢𝑗𝑖 , где j – номер точки в котором перемещение вычисляется при 𝑝𝑖 = 1, 𝑝𝑗 = 0, 𝑗 ≠ 𝑖, 𝑗 = 1,2, … Решая, систему уравнений 𝑈𝑷 = 𝒖, получаем вектор сил 𝑷 = (𝑝𝑖 , 𝑝2 , … , 𝑝𝑚 )′. В данном случае, u – это известные перемещения, полученные посредством геодезической съемки. На втором этапе решаем следующую систем дифференциальных уравнений (1.10)17

𝑚 𝑑𝑽(𝑠) = 𝐴𝑽(𝑠) + 𝑭 ∑ 𝑝𝑖 𝐿𝑖 (𝑠) 𝑑𝑠 𝑖=1 𝐺1 𝑽(𝑠0 ) = 0, 𝐺2 𝑽(𝑠𝑙 ) = 0 (1.11) В результате, решения краевой задачи (1.10), (1.11) получили НДС трубопровода, в котором реализуются перемещения, полученные из геодезической съемки. Линейная краевая задача (1.8) - (1.11) решается численно. Для решения используется метод ортогональной прогонки Годунова, который обеспечивает необходимую, точность вычислений. Так же стоит отметить, что, в случае сосредоточенных нагрузок, дельта функция в правой части системы (1.10), не вызывает дополнительных трудностей при численном решении, методом Рунге-Куты. Просто добавляется 𝐹𝑝 к вектору 𝑉(𝑠𝑖 ) при решении неоднородного уравнения в точках s, совпадающих с 𝑠𝑖 , и продолжается вычисление по формулам метода. 18

Глава 2. Метод численного решения Рассматривается линейная краевая задача вида: 𝑽′ (𝑠) − 𝐴𝑽(𝑠) = 𝒇(𝑠) (2.1) 𝐺1 𝑽(𝑠0 ) = 𝒃 (2.2) 𝐺2 𝑽(𝑠𝑙 ) = 𝒅 (2.3) Здесь 𝑽, 𝒇, 𝒃, 𝒅 – векторы размерностей соответственно 𝑚, 𝑚, 𝑚 − 𝑟, 𝑟. Размерности матриц 𝐴, 𝐺1 , 𝐺2 𝑙 × 𝑙, (𝑙 − 𝑟) × 𝑙, 𝑟 × 𝑙. Для решения краевой задачи известно много методов численного интегрирования, например, метод стрельбы. Но они не работают, если рассматривается плохо обусловленная задача. Краевую задачу называют хорошо обусловленной, если малые возмущения коэффициентов и правых частей уравнения и граничных условий – приводят к столь же малым по порядку изменениям решения задачи. Если V(s) – решение поставленной задачи, 𝑉1 (𝑠) – решение задачи с возмущениями правых частей и коэффициентов, то тогда задача будет хорошо обусловлена, если max ‖𝑉1 (𝑠) − 𝑉(𝑠)‖ ≤ 𝑀𝜀 𝑠0 ≤𝑠≤𝑧 при не очень больших значениях постоянной M. В противном случае задачу называют плохо обусловленной. Необходимость использования метода ортогональной прогонки возникает при решении линейных краевых задач с пограничным слоем. Т.е. когда решение довольно плавно меняется вдали от краев и очень сильно меняются при приближении к границе. Такие задачи возникают в теории оболочек. Например, при рассмотрении длинных тонких цилиндрических оболочек, говорят, что имеет место так называемый, краевой эффект. Фактически это проявление плохой обусловленности уравнений теории оболочек. При решении однородного уравнения возникают проблемы, связанные с плохой обусловленностью. Например, решением однородного уравнения 𝑽’ = 𝐴𝒔 с постоянной матрицей 𝐴 является 19

𝑙 𝑉𝑗 = ∑ 𝑐𝑖𝑗 𝜀𝑖 𝑒 𝜆𝑖(𝑠−𝑠0 ) 𝑖=1 Где 𝑐𝑖𝑗 – произвольные постоянные, 𝜀𝑖 – собственные вектора, 𝜆𝑖 – собственные числа матрицы 𝐴. Пусть 𝑅𝑒𝜆𝑗 ≤ 𝑅𝑒𝜆𝑙−1 < 𝑅𝑒𝜆𝑙 . 𝑙 𝑉𝑗 (𝑠0 ) = 𝑉𝑗0 , 𝑉𝑗0 = ∑ 𝑐𝑖𝑗 𝜀𝑖 𝑖=1 Тогда, решение однородного можно переписать в виде 𝑙−1 𝑉𝑗 = 𝑒 𝜆𝑙 (𝑠−𝑠0 ) (𝑐𝑙𝑗 𝜀𝑙 ) + ∑ 𝑐𝑖𝑗 𝜀𝑖 𝑒 (𝜆𝑖 −𝜆𝑙 )(𝑠−𝑠0 ) 𝑖=1 Но (𝜆𝑖 − 𝜆𝑙 ) < 0 поэтому при больших значениях 𝑠 − 𝑠0 , 𝑒 (𝜆𝑖−𝜆𝑙 )(𝑠−𝑠0 ) будет стремиться к 0, то есть все значения V будут пропорциональны 𝑒 𝜆𝑙(𝑠−𝑠0 ) (𝑐𝑙𝑗 𝜀𝑙 ), иначе говоря, произойдет сплющивание базиса. Что бы избежать данной проблемы очевидным является уменьшение отрезка 𝑠 − 𝑠0 . Одним из методов, базирующихся на этой идее является метод ортогональной прогонки Годунова [2] 20

2.1 Метод Ортогональной прогонки Разобьем отрезок интегрирования на m точек 𝑠0 = 𝑡0 < 𝑡1 < ⋯ < 𝑡𝑚 = 𝑠𝑙 . На левом конце решаем систему линейных алгебраических уравнений (2.2), в результате получим матрицу 𝐺0 решений размерности 𝑙 × 𝑟 + 1, где r первых столбцов – это общее решение однородного уравнения By=0, а (r+1) –ое – частное решение уравнения (2.2). Проводится ортогонализация и нормировка столбцов матрицы 𝐺0 (последний столбец не нормируется) по следующим формулам: 𝑗−1 𝒛𝑗 = 𝒈𝑗 − ∑(𝒈𝑖 , 𝒒𝑖 )𝒒𝑖 𝑖=1 𝒒𝑗 = 𝒛𝑗 √(𝒛𝑗 , 𝒛𝑗 ) , 𝑗 = 1, … , 𝑟, 𝑟 𝒒0 = 𝒈0 − ∑(𝒈0 , 𝒒𝑖 )𝒒𝑖 𝑖=1 Выразим zj в соотношениях (2.2), (2.3) через 𝒒i и перенесем все 𝒒i в левые части; получим j−1 √(𝐳j , 𝐳j )𝒒j + ∑(𝒈i , 𝒒i )𝒒i = 𝒈j , j = 1, … , r, i=1 r ∑(𝒈0 , 𝒒i )𝒒i + 𝒒0 = 𝒈0 . i=1 Эту систему соотношений можно переписать в матричной форме QW = G, где W= w11 0 ⋮ ( w21 w22 0 0 0 0 ⋯ ⋱ ⋯ wr1 wr2 wrr 0 wr+1,1 wr+1,2 , ⋮ wr+1,r 1 ) wjj = √(𝐳j , 𝐳j ), wij = (𝒈j , 𝒒i ) при i < j ≤ r, wr+1,i = (𝒈0 , 𝒒i ); Полученную матрицу назовем 𝑄0 столбцы этой матрицы будут являться 21

начальными условиями для задачи Коши, решая которую будет получена матрица 𝐺1 , размерности 𝑙 × 𝑟 + 1, где первые 𝒈1 , … , 𝒈𝑟 векторов будут решением уравнения 𝑽′ = 𝐴𝒔 а 𝑟 + 1вый вектор 𝒈0 – является решением уравнения (2.1). Далее повторяем проводим ортогонализацию, затем интегрируем и т.д., пока не дойдем до точки 𝑡𝑚 . После решения на отрезке [𝑡𝑚−1 ; 𝑡𝑚 ] получим матрицу 𝐺𝑚 , которую не нужно ортогонализовывать. На каждом из отрезков [t s , t s+1 ] искомое решение задачи записывается в виде r 𝐕(𝑡) = 𝒈s0 (t) + ∑ 𝜷sj 𝒈sj (t), j=1 или иначе 𝐲(t) = Gs (t)𝜷s, Где 𝜷s = (𝛽1𝑠 , … , 𝛽𝑟𝑠 , 1). Так как выполняется 𝒚(𝑡𝑠+1 ) = 𝐺𝑠 (𝑡𝑠+1 )𝜷𝑠 и выполняется условие непрерывности (𝒚(𝑡𝑗−0 ) = 𝒚(𝑡𝑗+0 )), то 𝒚(𝑡𝑠+1 ), так же может быть найдена по следующей формуле: 𝒚(𝑡𝑠+1 ) = 𝐺𝑠+1 (𝑡𝑠+1 )𝜷𝑠+1 , то 𝜷s+1 = Ws+1 𝜷s . Значение 𝜷m может быть найдено из правого граничного условия (2.3) D𝐲(𝑥1 ) = D(Gm (t m )𝜷m ) = 𝐝. А затем последовательно определяются 𝛃m−1 , … , 𝜷0 . 22

Глава 3. Анализ результатов 3.1 Алгоритм нахождения значений нагрузок При единичных нагрузках (𝑝𝑖 = 1, 𝑝𝑗 ≠ 0, 𝑖 = 1. .3, 𝑗 ≠ 𝑖), решается система (1.10) с граничными условиями (1.11). На каждом i-ом шаге получается вектор 𝒒𝒊 (𝑘) = 𝑉6 (𝑘) − 𝑉8 (𝑘) 2 = 𝑢𝑟,1 (𝑘) − 𝜖2,1 (𝑘) 2 номер k – совпадает с номером точки в которой известно перемещение и меняется от 1 до 3. Из векторов q формируется матрица податливости 𝑄 = (𝒒1 , 𝒒2 , 𝒒3 )′. Вектор заданных смещений u и матрица податливости Q позволяют найти значения нагрузок в точках с известными перемещениями. Для этого необходимо решить систему алгебраических уравнений: 𝑄𝒑 = 𝒖. Полученные p используются для решения системы (1.10). Расчет проводился для прямолинейного участка трубопровода длины 𝑙̃ = 𝑙 𝑅 = 45.1022. Величины с тильдой – записаны в безразмерном виде. Вертикальные перемещения поперечного сечения трубопровода были 𝑠 𝑢 измерены в трех точках(𝑠̃ = , 𝑢̃ = ) 𝑅 𝑅 𝑠̃1 = 7.75, 𝑠̃2 = 22.83, 𝑠̃3 = 36.65, 𝑢̃1 = −0.05, 𝑢̃2 = 0.035, 𝑢̃3 = −0.085. (3.1) 𝑢̃𝑖 - перемещение трубопровода в точке 𝑠̃𝑖 . Эти точки отмечены на графике перемещений вдоль трубопровода, рис. 1, цифрами 1, 2 и 3 соответственно. 23

Рис. 1. График полученных перемещений Значения нагрузок, которые обеспечивают перемещения (3.1): А) в случае сосредоточенных нагрузок, 𝑝̃ = 𝑝 𝐸𝑅 : 𝑝̃1 = −0.595 ∙ 10−5 , 𝑝̃2 = 0.452 ∙ 10−5 , 𝑝̃3 = −0.780 ∙ 10−5 Б) в случае кусочно-постоянной аппроксимации нагрузок в окрестности 𝑝 замеров на участках длиной 0.0022𝑙̃, 𝑝̃ = : 𝐸 𝑝̃1 = −0.591 ∙ 10−4 , 𝑝̃2 = 0.452 ∙ 10−4 , 𝑝̃3 = −0.777 ∙ 10−4 𝑝 В) в случае интерполяции линейным сплайном, 𝑝̃ = : 𝐸 𝑝̃1 = −0.286 ∙ 10−5 , 𝑝̃2 = 0.265 ∙ 10−5 , 𝑝̃3 = −0.347 ∙ 10−5 𝑝 Г) в случае интерполяции полиномом Лагранжа, 𝑝̃ = : 𝐸 𝑝̃1 = −0.257 ∙ 10−5 , 𝑝̃2 = 0.165 ∙ 10−5 , 𝑝̃3 = −0.289 ∙ 10−5 Используя полученные нагрузки, находим НДС трубопровода в четырех случаях нагружения. Заметим, что случай Б при уменьшении участков приложения нагрузок должен переходить в случай А. Действительно, если посчитать усилия от распределенных кусочно-постоянных нагрузок (случай Б) около точек 𝑠𝑖 , получим 𝑃𝑖 = 𝑝̃𝑖 ∙ 0.0022𝑙̃ = 0.1𝑝̃𝑖 , что практически совпадает со значениями сосредоточенных усилий в случае 24

А. Этот факт свидетельствует в пользу достоверности получаемых численных результатов. Для решения поставленной задачи создана программа в прикладном пакете MatLab. Текст программы с комментариями приведен в приложении. Программа тестировалась на задаче, рассмотренной в работе [9] (случай сосредоточенных нагрузок). 25

3.2 Графики напряжений Компоненты нагружения усилия T и моменты M определяются через полученный вектор V(s). 𝑀,1 (𝑠) − 𝑀1,1 (𝑠)) ( 𝜋𝑅 𝑇1 = 𝑅 𝑀1 = 𝑀1,1 (𝑠) 𝑇2 = 𝜈𝑀1 + 𝐵(1 − 𝜈 2 )𝜖2,1 (𝑠) (3.2) 𝑀2 = 𝜈𝑀1 + 𝐷(1 − 𝜈 2 )𝜖2,1 (𝑠) Напряжение определяется по формулам: 𝜎 𝜎̃ = 𝐸 𝑇𝑖 𝜎𝑇𝑖 = ℎ 6𝑀𝑖 𝜉 𝜎𝑀𝑖 = ± 2 ∙ ℎ ℎ/2 Где соответствующие 𝑇𝑖 , 𝑀𝑖 – получены по (3.2). Поскольку все данные о перемещениях известны на срединной поверхности, то построим графики напряжений на срединной поверхности. Рис. 2. График изгибных напряжений 𝜎𝑀1 26

Рис. 3. График изгибных напряжений 𝜎𝑀2 Рис. 4. График тангенциальных напряжений 𝜎𝑇1 27

Рис. 5. График тангенциальных напряжений 𝜎𝑇2 Но для анализа НДС недостаточно результатов представленных на графиках, рис. 2-5, необходимо рассмотреть другие данные, а именно напряжения на лицевых поверхностях трубопровода. Формулы для расчета напряжений 𝜎𝑖 = 𝑇𝑖 ℎ + 6𝑀𝑖 ℎ2 ∙ 𝜉 ℎ/2 ℎ ℎ ; 𝑖 = 1,2; 𝜉 ∈ [− ; ] 2 2 𝜉= ℎ 2 для верхней ℎ поверхности и 𝜉 = − для нижней. 2 28

3.3 Графики напряжений на лицевых поверхностях трубопровода Для того что бы оценить напряжение по сравнению с максимально допустимым, необходимо построить грфики напряжений на внешней и внутренней поверхностях трубопровода, для всех случаев нагружения. Так же построим на графике предельные значения напряжений взятые из технического задания. Модуль юнга E = 206000 МПа, максимально допустимое напряжение 𝑅𝑚𝑎𝑥 = 563 МПа, безразмерная величина максимально допустимого напряжения: 𝜎𝑚𝑎𝑥 = Максимально допустимые 𝑅𝑚𝑎𝑥 = 0.00273 𝐸 напряжения на графиках обозначим пунктирными линиями. Рассмотрим полученные графики Мередиальных напряжений. Формулы для рассчета мередиальныхмнапряжений. 𝑆𝑖𝑔𝑚𝑎 1− = 𝜎𝑇1 − 𝜎𝑀1 𝑆𝑖𝑔𝑚𝑎 1+ = 𝜎𝑇1 + 𝜎𝑀1 Рис. 6. График меридиональных(осевых) напряжений на внутренней лицевой поверхности трубы 29

Рис. 7. График меридиональных(осевых) напряжений на внешней лицевой поверхности трубы Из графиков, рис. 6, 7 видно что в зависимости от способа нагружения получаются разные результаты. Так скачкообразный график сосредоточенных нагрузок, выходит за пределы максимально допустимого напряжения во всех трех точках, перемещение которых было известно. График кусочно постоянной, около замеров функции ведет себя лучше, чем график сосредоточенных нагрузок, пики явно меньше, однако он тоже выходит за границы. Графики распределенных нагрузок, ведут себя гораздо лучше, оставаясь в пределах допустимого. Только в точке с координатой 36.65 видны выходы за пределы допустимого при всех способах нагружения. Видимо данной точке следует уделить особое внимание, и в ней возможно повреждение трубопровода. Рассморим полученные графики Окружных напряжений. Формулы для рассчета окружных напряжений 𝑆𝑖𝑔𝑚𝑎 2− = 𝜎𝑇2 − 𝜎𝑀2 𝑆𝑖𝑔𝑚𝑎 2+ = 𝜎𝑇2 + 𝜎𝑀2 30

Рис. 8. График окружных напряжений на внутренней лицевой поверхности трубы Рис. 9. График окружных напряжений на внешней лицевой поверхности трубы На графиках, рис. 8,9 видно, что значения окружных напряжений меньше меридиональных и не выходят за пределы максимально допустимых, кроме случая сосредоточенных нагрузок в точке с координатой 36.65 на внешней лицевой поверхности трубы. 31

Выводы 1. Создана программа решения плохо обусловленных линейных краевых задач методом ортогональной прогонки Годунова. 2. Построен обеспечивающих алгоритм заданные определения (совпадающие параметров с нагружения, экспериментальными) перемещения во фрагменте трубопровода. 3. Определено НДС для различных способов нагружения от "жесткого" (сосредоточенные нагрузки в точках снятия перемещений) до "мягкого"(непрерывно распределенная нагрузка). 4. Показано, что измеренные в некоторых точках трубопровода перемещения не позволяют однозначно судить о величине возникающих в конструкции напряжений. 5. В развитие работы можно рассмотреть не только вертикальные нагрузки, но нагрузки, действующие в других направлениях, которые также могут влиять на НДС фрагмента трубопровода. 32

Заключение Разработана методика расчета изменения НДС фрагмента трубопровода по измеренным вертикальным смещениям в некоторых точках. 33

Список используемой литературы 1. Балабух Л.И. Изгиб и кручение конических оболочек // Труды ЦАГИ. 1946. № 577. 63 с. 2. Бахвалов Н.С. Численные методы (анализ, алгебра, обыкновенные дифференциальные уравнения). Главная редакция физико- математической литературы изд-ва «Наука», 1975. С. 571-576. 3. Бидерман В.Л. Механика тонкостенных конструкций. Изд. Машиностроение, 1977 С. 277-282. 4. Власов В.З. Общая теория оболочек и ее приложения в технике. М.: Гостехиздат, 1949. 784 с. 5. Годунов С.К. Метод ортогональной прогонки для решения систем разностных уравнений // Журнал Вычислительной математики и математической физики, 1962. Т. 2. № 6. С. 972–982. 6. Жилин П.А. Прикладная механика основы теории оболочек // Изд. Политехнического университета, 2006 С.10-14 7. Завьялов А.П., Гусейнов К.П., Егоров С.И., Лопатин А.С. Опыт применения программных комплексов на основе метода конечных элементов для оценки напряженно-деформированного состояния магистральных газопроводов, проложенных в особых климатических условиях // Проблемы сбора, подготовки и транспорта нефти и нефтепродуктов. Институт проблем транспорта энергоресурсов Республики Башкортостан, 2014. С.67 – 70. 8. Кабриц С.А., Шамина В.А. Изгиб оболочки вращения поперечной силой и моментом // Вестник. СПБГУ, 2014. С.261 – 269. 9. Kabrits S.A., Shamina V.A. Cilindrical shell under the action of the ring load // Seventh Polyakhov’s Reading, 2015 International conference on Mechanics, 2015. – P. 1-3. 10.Кирхгоф Г. Механика. М.: Изд-во АН СССР, 1962. 402 с. 11.Love A.E.H. A Treatise on the Mathematocal Theory of Elasticity. V. II. 34

Cambridge, 1893. 327 p. 12.Лурье А.И. Общая теория упругих тонких оболочек // ПММ. 1940. Т. 4. № 2. С. 7–34. 13.Лурье А.И. Равновесие упругой симметрично-нагруженной сферической оболочки // ПММ. 1943. Т.7. № 6. С. 393–404. 14.Лурье А.И. Статика тонкостенных упругих оболочек. М.: Гостехиздат, 1947. 252 с. 15.Новожилов В.В., Финкельштейн Р.М. О погрешности гипотез Кирхгофа–Лява в теории оболочек // ПММ. 1943. Т.7. №. 5. С. 323– 330. 16.Новожилов В.В. Теория упругости. Л.: 1958. 369 с. 17.Новожилов В.В. Теория тонких оболочек. Л.: Судопромгиз, 1962. 432 с. 18.Черных К.Ф. Линейная теория оболочек. Ч.1,2. Л.: Изд-во Ленинградского Университета. Ч.1, 1962. 272с. Ч.2, 1964. 395с. 35

Приложение Первая функция, которую необходимо рассмотреть – это основной модуль ортогональной прогонки. Входными параметрами в данной функции являются число точек интегрирования и длина отрезка интегрирования. function y=MainFunc(m,l) [G0,D,d,s]=Gran(); Задание граничных условий Задание переменных, содержатся в которых координаты точек global xx xnom ; ортогонализации (xx), и номера точек, в [x,xnom]=razbivka(m,l,s); которых известно перемещение (xnom). xnom=[1,xnom,length(x)]; Функция razbivka определяет точки таким xx=x; образом, чтобы точки, в которых известны перемещениями, являлись точками ортогонализации N=size(G0); Прямой ход метода. В переменных V, U GG=G0; содержатся [V,U]=Pro(G0); ортогонализации и матрицы W. матрицы Q в точках Mv=size(V); Mu=size(U); Обратный ход метода. Из V берется for i=1:N(2) последняя G(:,i)=V(:,Mv(2)-N(2)+i); end K=Mu(2)/N(2); матрица, для которой находится 𝛽𝑚 из правого граничного условия. Затем последовательно вычисляются все значения 𝛽𝑖 , 𝑖 = 1. . 𝑚 − 1 и соответственно значения искомых P=D*G; функций в каждой точке ортогонализации t=size(P); 𝑦𝑖 , 𝑖 = 1. . 𝑚. for i=1:t(2)-1 S(:,i)=P(:,i); 36

end d=-P(:,t(2)); be(:,K)=S\d; be(t(2),:)=1; y(:,K)=G*be(:,K); for j=2:(Mv(2)/N(2))-1 for i=1:N(2) G(:,i)=V(:,Mv(2)-j*N(2)+i); W(:,i)=U(:,Mu(2)-(j-1)*N(2)+i); end be(:,K-(j)+1)=W\be(:,K-(j)+2); y(:,(Mv(2)/N(2))-j+1)=G*be(:,K-(j)+1); end for i=1:N(2) W(:,i)=U(:,N(2)+i); end be(:,1)=W\be(:,2); G0=ortogon(G0); y(:,1)=G0*be(:,1); Запись полученных значений y в матрицу выходных y=[x;y]; данных. добавляются Так же координаты в нее точек y=y'; ортогонализации. После этого матрица end транспонируется для более удобного вывода. Рассмотрим теперь следующий модуль, модуль Pro, в котором производится прямой ход прогонки. Входной параметр – это матрица 𝐺0 , полученная при решении левого граничного условия 𝐵𝑦 = 𝑏. function [V,U]=Pro(G0) global xx xnom tip qq sknom; 37

x=xx; Задание начальных данных GG=G0; for i=1:length(xx)-1 sknom=i; [GG,W]=ortogon(GG); if i>1 Ортогонализация и нормировка матрицы входной, и запись её в массив матриц. V=[V,GG]; U=[U,W]; else V=GG; U=W; End SK=[-2*pi;0;1;0;0;0;0;0]; Skachok=0*SK; if (tip==1) Вычисление величины скачка для случая сосредоточенных нагрузок, при условии, что в системе уравнений в неоднородной части q21=-qn1, qn1=1. for j=1:length(xnom)-1 if((i==xnom(j)))%||(i+1==xnom(j))) Skachok=qq(j)*SK; end end end GG=Mfun(GG,[x(i),x(i+1)],Skachok); end [GG,W]=ortogon(GG); V=[V,GG]; Интегрирование на отрезке [x(i);x(i+1)], ортогонализация и запись полученной матрицы в массив матриц. В случае сосредоточенных нагрузок скачок не U=[U,W]; нулевой, и он добавляется к правой части End системы Д.У. при решении неоднородного уравнения. 38

Рассмотрим функция которая получает необходимые напряжения 𝜎𝑇1 , 𝜎𝑇2 , 𝜎𝑀1 , 𝜎𝑀2 . Входными параметрами являются число точек интегрирования, длина отрезка интегрирования, а также номер от 1 до трех, указывающий способ интерполирования (1 – сосредоточенные, 2 – линейный сплайн, 3 – полином Лагранжа) function [up,Y]=Napr_Skach(m,l,t) [h0,R,diam_vneshn,nu,E_yung,D,B]=Geom(); global xnom qq tip; tip=t; q=[0,1,0,0,0; Задание переменных, и единичных нагрузок в точках ортогонализации. 0,0,1,0,0; 0,0,0,1,0]; U=zeros(3); for i=1:3 qq=q(i,:); Формирование матрицы податливости U. y=MainFunc(m,l); for j=1:3 U(i,j)=y(xnom(j+1),9)-y(xnom(j+1),7)/2; end end UU=U'; u=[-0.05*R;0.035*R;-0.085*R]; p=UU\u q=[0,p(1),p(2),p(3),0]; Задание перемещений нагрузок, необходимых и расчет для их реализации qq=q; 39

y=MainFunc(m,l); Вычисление вектора y при полученных нагрузках p si=size(y); Расчет for i=1:si(1) запись их в одну матрицу, в 𝜎𝑇1 , 𝜎𝑇2 , 𝜎𝑀1 , 𝜎𝑀2 и которой первый столбец – T11(i)=(y(i,3)/(pi*R)-y(i,5))/R; координаты M11(i)=y(i,5); T21(i)=nu*M11(i)+B*(1-nu^2)*y(i,7); M21(i)=nu*M11(i)+D*(1-nu^2)*y(i,7)/R; ортогонализации, точек 2-5 напряжения T1..M2, и 6 столбец перемещения end up=[y(:,1),T11'/(h0*E_yung),6*M11'/(h0^2*E_yung),T21'/(h0*E_yung), 6*M21'/(h0^2*E_yung),y(:,9)-y(:,7)/2]; plot(up(:,1),up(:,6)); end Кроме представленных используются и другие функции Функция, которая возвращает данные о геометрических характеристиках оболочки function [h0,R,diam_vneshn,nu,E_yung,D,B,q21,qn1,q11]=Geom() h0=0.0105; diam_vneshn=0.720; R=(diam_vneshn-h0)/2; nu=0.3;E_yung=2060;%206e9% D=E_yung*h0*h0*h0/(12*(1-nu^2)); B=E_yung*h0/(1-nu^2); q21=0; qn1=0; q11=0; end Функция граничных условий, возвращает матрицу решений левого граничного условия и матрицы, которыми описывается правое граничное условие. function [B,D,d,s]=Gran() B=[ 0,0,0,0,1,0,1,0; 40

0,0,0,0,0,1,0,0; 0,0,0,1,0,0,0,0; 0,0,0,0,0,0,0,1]; D=[ 0,0,0,0,1,0,1,0; 0,0,0,0,0,1,0,0; 0,0,0,1,0,0,0,0; 0,0,0,0,0,0,0,1]; d=[0;0;0;0]; b=d; Levkray – Функция, решающая B=Levkray(B,b); систему уравнений 𝐵𝑥 = 𝑏 End Функция, решающая систему уравнений Bx=b для левого граничного условия, и возвращающая фундаментальную систему решений. function G=Levkray(BB,b) B=BB; n=size(B); for i=1:n(2) x(i)=i; y(i)=0; end for i=1:n(1) max=B(i,i); maxi=i; for j=i:n(2) if(abs(B(i,j))>abs(max)) max=B(i,j); maxi=j; 41

end end C=B(:,i); B(:,i)=B(:,maxi); B(:,maxi)=C; c=x(i); x(i)=x(maxi); x(maxi)=c; B(i,:)=B(i,:)/B(i,i); for j=1:n(1) if(j~=i) B(j,:)=B(j,:)-B(i,:)*B(j,i); end end end k=1; for j=(n(1)+1):n(2) for i=1:n(2) y(i)=0; end for i=1:n(1) y(x(i))=-B(i,j); end y(x(j))=1; G(:,k)=y'; k=k+1; end G(:,k)=BB\b; End Функция, ортогонализации и нормировки столбцов матрицы, последний 42

столбец не нормируется. function [Q,W]=ortogon(G) N=size(G); for i=1:N(2) p=0; for j=1:i-1 p=p+dot(G(:,i),Q(:,j))*Q(:,j); end z(:,i)=G(:,i)-p; Q(:,i)=z(:,i)/(sqrt(dot(z(:,i),z(:,i)))); end Q(:,N(2))=z(:,N(2)); W=zeros(N(2)); for i=1:N(2)-1 for j=(i+1):N(2) W(i,j)=dot(G(:,j),Q(:,i)); end W(i,i)=sqrt(dot(z(:,i),z(:,i))); end W(N(2),N(2))=1; End Функция интегрирования системы ДУ, в которой в случае скачка, он добавляется при решении неоднородной части уравнения. function M=Mfun(GG,tspan,Skachok) global Nom Nom=0; N=size(GG); for i=1:N(2)-1 43

[t,y]=RunKut(tspan,GG(:,i)); n=size(y); M(:,i)=y(n(1),:); end Nom=1; [tn, yn] = RunKut(tspan,GG(:,N(2))+Skachok); n=size(yn); M(:,N(2))=yn(n(1),:); End Функция, реализующая метод Рунге-Кута 4 порядка. function [x,y]=RunKut(x_span,y0) X=x_span(2); xn=x_span(1); yn=y0; x=xn; y=yn; h=(X-xn)/100; while (xn<X) if((X-xn)<h) h=X-xn; end [xn,yn]=Step(xn,yn,h); x=[x,xn]; y=[y,yn]; end x=x'; y=y'; end Функция, реализующая шаг, по методу Рунге-Куты 4-ого порядка. 44

function [x,y]=Step(x,y,h) k1=h*RKST(x,y); k2=h*RKST(x+h/2,y+k1/2); k3=h*RKST(x+h/2,y+k2/2); k4=h*RKST(x+h,y+k3); x=x+h; y=y+(k1+2*k2+2*k3+k4)/6; end Функция, в которой содержатся уравнения цилиндрической оболочки, которые решаются методом Рунге-Куты. function z=RKST(x,y) global Nom tip; [h0,R,diam_vneshn,nu,E_yung,D,B,q21,qn1,q11]=Geom(); R2=R^2; R3=R2*R; nu2=nu^2; if ((Nom ==1)&&(tip~=1)) qn=inter(x,tip); qn1=qn; q21=-qn1; q11=0; end g(1)= pi*R*(q21-qn1); g(2)= y(1)-pi*R2*q11; g(3)= -(-B*y(6)*pi*R2+B*nu2*y(6)*pi*R2-nu*y(2)+pi*D*nu2*y(6)- pi*D*y(6)+q21*R3*pi)/pi/R2; g(4)= (- 0.2e1*D*pi*B*y(5)*R+0.2e1*D*pi*B*y(5)*nu*R+y(3)*pi*B*R2+0.4e1*y(1)*D +y(1)*B*R2)/(0.4e1*D+B*R2)/R/pi; g(5)= -(R3 * B * pi * y(4) - D * y(2) + D * y(4) * pi * R) / D / R3 / B / pi; g(6)= (-B * y(5) * R + B * y(5) * nu * R - 2 * y(3)) / (-4 * D + 4 * D * nu 45

B * R2 + B * R2 * nu); g(7)= -(-y(4) * R + D * nu * y(6)) / D / R; g(8)= -y(7); z=g'; end Функция, которая возвращает значение нагрузки в зависимости от способа интерполирования для случая распределенной нагрузки. function qn=inter(x,k) global xx xnom qq sknom; if(k==0) Кусочно-заданная сила for i=1:length(xnom)-1 if(abs(x-xx(xnom(i)))<=0.5) qn=qq(i); break; else qn=0; end end elseif(k==2) Линейный сплайн qn=0; for j=1:length(xnom) if((x>xx(xnom(j)))&&(x<=xx(xnom(j+1)))) qn=qq(j)+(qq(j+1)-qq(j))/(xx(xnom(j+1))-xx(xnom(j)))*(xxx(xnom(j))); end end 46

elseif (k==3) Полином Лагранжа s=0; for i=1:length(xnom) L=1; for j=1:length(xnom) if(j~=i) L=L*(x-xx(xnom(j)))/(xx(xnom(i))-xx(xnom(j))); end end s=s+qq(i)*L; end qn=s; else qn=0; end end Функция разбиения отрезка интегрирования на точки, так, чтобы точки, в которых известны перемещения, являлись точками ортогонализации. function [x,minnom]=razbivka(m,l,s) x=0:l/m:l; min=l; minnom=0; for k=1:length(s) min=l; minnom(k)=2; for i=2:length(x)-1; if(min>abs(x(i)-s(k))) min=abs(x(i)-s(k)); minnom(k)=i; 47

end end x(minnom(k))=s(k); end End 48

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв