«ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ УПРУГОГО ДЕФОРМИРОВАНИЯ КОМПОЗИТНЫХ ОБОЛОЧЕК ВРАЩЕНИЯ ...»

РОССИЙСКАЯ АКАДЕМИЯ НАУК

СИБИРСКОЕ ОТДЕЛЕНИЕ

Институт вычислительных технологий

На правах рукописи

УДК 539.3

Юрченко Андрей Васильевич

ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ

УПРУГОГО ДЕФОРМИРОВАНИЯ

КОМПОЗИТНЫХ ОБОЛОЧЕК ВРАЩЕНИЯ

05.13.18 — математическое моделирование, численные методы и комплексы программ Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель к.ф.-м.н., доцент С. К. Голушко Новосибирск - Оглавление Введение 1 Методы решения краевых задач механики композитных пластин и оболочек вращения 1.1 О сведении двумерных краевых задач к одномерным... 1.2 Особенности систем дифференциальных уравнений при решении краевых задач..................... 1.3 Методы решения краевых задач для систем обыкновенных дифференциальных уравнений................ 1.3.1 Метод начальных параметров............. 1.3.2 Метод дискретной ортогонализации......... 1.3.3 Метод сплайн-коллокации............... 2 Алгоритм решения краевых задач для систем обыкновенных дифференциальных уравнений 2.1 Проблемы вычисления векторов начальных данных и решения многоточечных задач................. 2.2 Обеспечение устойчивости расчетов............. 2.3 Схема алгоритма........................ 2.4 Анализ эффективности алгоритма при решении задачи изгиба слоистых кольцевых пластин.............. 2.5 Обеспечение точности расчетов с использованием неравномерных сеток.......................... 3 Определяющие соотношения статики упругих композитных оболочек вращения 3.1 Моделирование свойств и критерии прочности полиармированных композитов..................... 3.1.1 Модель В.В. Болотина................. 3.1.2 Модели Ю.В. Немировского.............. 3.1.3 Критерии прочности и начального разрушения композиционных материалов............... 3.2 Сравнительный анализ расчетных характеристик композиционных материалов с экспериментальными данными.. 3.3 Исходные и разрешающие системы уравнений теории оболочек вращения......................... 3.3.1 Исходные уравнения и соотношения......... 3.3.2 Разрешающие системы уравнений.......... 4 Напряженно-деформированное состояние рефлектора параболической антенны 4.1 Постановка задачи....................... 4.2 Рефлектор под действием осесимметричного нагружения собственным весом и температурой.............. 4.3 Рефлектор под действием собственного веса и ветровой нагрузки.............................. 4.4 Рефлектор под действием собственного веса, температурной и ветровой нагрузок.................... 4.5 Анализ достоверности численных решений......... 5 Особенности поведения и начальное разрушение армированных куполов 5.1 Постановка задачи....................... 5.2 Купол под действием осесимметричного нагружения собственным весом......................... 5.3 Купол под действием собственного веса и давления ветра. 5.4 Купол под действием собственного веса, ветровой и температурной нагрузок....................... 5.5 Анализ достоверности численных решений......... 6 Влияние анизотропии материала на деформирование резинокордной тороидальной оболочки 6.1 Влияние выбора модели композиционного материала и теории оболочек на результаты расчетов напряженно–деформированного состояния оболочки............................. 6.2 Влияние анизотропии и неоднородности материала на поведение оболочки........................ 6.3 Об использовании несимметричных схем армирования... Заключение Литература Приложение Введение Фундаментальная задача научных исследований — выявление причинно– следственных связей, общих тенденций и закономерностей. Так как проведение натурных экспериментов затрудняется их дороговизной и сложностью, проблемами при обеспечении исследователя желаемым количеством измеряемых параметров, а в ряде случаев невозможностью реализации, то моделирование процессов становится одним из наиболее распространенных методов исследования объектов и явлений различной природы. Особая роль при этом отводится вычислительному эксперименту.

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

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

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

В качестве связующего используются полимерные, углеродные, металлические или керамические матрицы. Армирующими элементами являются стеклянные, углеродные, борные или органические, а также ряд металлических волокон. Перспективно, также, использование в качестве арматуры нитевидных кристаллов (усов). В отличие от обычных волокон, обладающих прочностью не превышающей 1-3% от модуля упругости E, усы, вследствие малых поперечных размеров и высокой степени совершенства, достигают прочности в 5-15% E. В таблице 1 приведены механические свойства характерных представителей разных классов матриц и волокон [10, 66, 72], используемые в работе. Здесь – удельная плотность, E – модуль упругости, – коэффициент линейного температурного расширения материала, а – предел прочности (текучести).

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

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

Методами математического моделирования задачи определения НДС композитных элементов конструкций приводятся к решению пространственной задачи механики деформируемого твердого тела, сформулированной в виде краевых задач для систем дифференциальных уравнений с частными производными. Решение таких задач в аналитическом виде возможно лишь в простейших случаях, а наиболее распространенным подходом является применение численных алгоритмов. Среди универсальных алгоритмов, позволяющих решать в пространственной постановке задачи механики сплошных сред, в том числе механики композитных конструкций, выделяются методы конечных и граничных элементов. Они завоевывают все большую популярность, вследствие существенного прогресса, достигнутого в разработке высокопроизводительной вычислительной техники. Однако их использование приводит к очень большим объемам вычислений, особенно при решении задач механики композитных конструкций в пространственной постановке. Это обусловлено эффектами анизотропии материала, наличием в системах уравнений, описывающих НДС, малых и больших параметров, что приводит к появлению ярко выраженных краевых эффектов и заставляет существенно измельчать сетки элементов в областях погранслоев. Получаемые объемы вычислений столь высоки, что даже на современных высокопроизводительных системах для решения обозначенных задач требуется существенное количество времени, что, в частности, накладывает существенные ограничения на проведение вариантных расчетов. При этом, достоверность получаемых в результате расчета чисел бывает очень сомнительна [64].

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

Цель диссертационной работы заключается в:

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

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

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

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

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

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

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

Исследования выполнялись в соответствии с планами научно-исследовательских работ Института вычислительных технологий СО РАН по теме "Теоретические исследования моделей и разработка эффективных численных методов решения нелинейных задач математической физики" (номер государственной регистрации 01. 2. 00 313336), поддерживались грантами: Федеральной целевой программы "Интеграция" (грант № 274); Президента РФ для поддержки молодых российских ученых и ведущих научных школ РФ (№ НШ–2314. 2003. 1).

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

Апробация работы. Основные положения диссертационной работы докладывались и обсуждались на следующих конференциях: V Всероссийской научно–технической конференции "Механика летательных аппаратов и современные материалы" (Томск, 1998); II и V Сибирских школах–семинарах "Математические проблемы механики сплошных сред" (Новосибирск, 1998; 2001); XXXVII, XXXVIII Международных научных конференциях "Студент и научно–технический прогресс" (Новосибирск, 1999; 2000); V, VI и VII научных конференциях "Современные методы математического моделирования природных и антропогенных катастроф" (Красноярск, 1999; 2001; 2003); XVI школе–семинаре "Информационные технологии в задачах математического моделирования" в рамках научных мероприятий "Вычислительные технологии — 2000" (Новосибирск, 2000); Конференции молодых ученых, посвященной 10–летию ИВТ СО РАН (Новосибирск, 2000); Международной конференции "Современные проблемы прикладной математики и механики: теория, эксперимент и практика", посвященной 80-летию академика Н.Н. Яненко (Новосибирск, 2001); XVII, XVIII, XIX Межреспубликанских конференциях по численным методам решения задач теории упругости и пластичности (Новосибирск, 2001; Кемерово, 2003; Бийск, 2005); Международных конференциях молодых ученых по математике, математическому моделированию и информатике (Новосибирск, 2001;

2002); Международной конференции "Вычислительные и информационные технологии в науке, технике и образовании" (Казахстан, Алматы, 2004).

В полном объеме материалы диссертации докладывались и обсуждались на Объединенном семинаре "Информационно–вычислительные технологии" Института вычислительных технологий СО РАН, Новосибирского государственного университета и Новосибирского государственного технического университета (руководители — академик Ю.И. Шокин и д.ф.-м.н., профессор В.М. Ковеня; Новосибирск, 2005); семинаре "Проблемы математического и численного моделирования" Института вычислительного моделирования СО РАН (руководитель — чл.-корр.

РАН В.В. Шайдуров; Красноярск, 2005).

Публикации. По теме диссертации опубликовано 14 статей в научных журналах и сборниках трудов конференций, а также тезисы докладов на научных конференциях [17]–[31], [74], [75], [80].

Структура и объем диссертации. Диссертация состоит из введения, 6 глав, заключения и приложения. Общий объем диссертации составляет 164 стр., включая 43 рисунка и 13 таблиц. Список литературы содержит 83 наименования.

Содержание работы.

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

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

Вторая глава посвящена разработке алгоритма решения многоточечных краевых задач для систем обыкновенных дифференциальных уравнений (ОДУ) на основе МДО. Предложено решение проблемы построения векторов начальных данных для многоточечных краевых задач. Описана соответствующая модификация МДО. Численно исследована проблема обеспечения устойчивости счета при использовании МДО.

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

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

Четвертая глава посвящена моделированию поведения главного зеркала параболической антенны. Исследовано влияние структурных и механических характеристик КМ на НДС параболического рефлектора, выполненного в виде тонкой композитной оболочки и подверженного действию собственного веса, ветровой и температурной нагрузок и их комбинаций. Показано, что спектральные характеристики матриц систем ОДУ, полученные из исходной системы с частными производными методом разделения переменных, существенно зависят от механических и структурных параметры материала, номера рассчитываемых гармоник в разложении решения в ряды Фурье. Проведены параметрические исследования характера поведения рефлектора в зависимости от структуры КМ и механических характеристик фазовых материалов. Проведен анализ влияния вида нагружения на НДС рефлектора. Достоверность численных результатов обеспечена совпадением результатов, полученных двумя принципиально разными численными методами.

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

В шестой главе, изучено влияние неоднородности и анизотропии КМ на деформирование резинокордной тороидальной оболочки, влияние выбора модели КМ и теории оболочек на результаты расчетов напряженно– деформированного состояния оболочки. Показана необходимость учета нелинейных слагаемых при расчете резинокордных тороидальных оболочек. Показано существенное влияние на результаты расчетов выбора модели КМ. Рассмотрены схемы армирования оболочки с несимметричным, относительно меридиана, строением. Показано, что их применение позволяет реализовать значительно более широкий спектр НДС в оболочке, в сравнении симметричными схемами.

В заключении сформулированы основные результаты работы.

Глава Методы решения краевых задач механики композитных пластин и оболочек вращения Вопросы численного решения краевых задач определения НДС пластин и оболочек рассматриваются, в той или иной мере, во многих научных трудах, посвященных теориям оболочек и пластин [5, 8, 11, 33, 39, 43, 47, 67, 68, 71] и др. Это связано в первую очередь с тем, что получение решений соответствующих краевых задач аналитическими методами сопряжено с существенными, порою непреодолимыми трудностями.

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

Широкое применение к задачам определения НДС оболочек и пластин получили вариационно–проекционные методы (метод конечных элементов, метод граничных элементов) [11, 47, 67, 68, 69, 71] и др. Их применение сопряжено с теми же проблемами, что и в случае решения пространственных задач теории упругости: большим объемом вычислений, трудностями при решении задач механики композитных конструкций.

Применение проекционно–разностного подхода позволяет в ряде случаев решать задачи классической теории оболочек и уточненной теории типа Тимошенко [8, 41, 44, 77] и др. Среди существенных ограничений, накладываемых указанными подходами, отметим трудности при решении задач по уточненным теориям оболочек, учитывающих поперечные сдвиги на основе гипотез о их нелинейном распределении по толщине оболочки [5, 33] и др. Системы дифференциальных уравнений, возникающие при использовании таких теорий, имеют более высокий порядок, чем в случае классической теории и уточненной теории типа Тимошенко, а в структурах их решений присутствуют быстроизменяющиеся компоненты, например, экспоненциального характера или типа модифицированных функций Бесселя. При решении соответствующих краевых задач вариационно–проекционными, проекционно- и конечно–разностными методами указанные факторы могут приводить либо к неоправданно высоким вычислительным затратам, либо к невозможности получения приемлемой точности расчетов, либо к непреодолимой неустойчивости счета в принципе.

Важное место среди подходов к решению двумерных краевых задач теории пластин и оболочек занимают различные методы понижения размерности задачи, например, метод сплайн–аппроксимации функций в одном из координатных направлений [37, 38], для задач расчета замкнутых в окружном направлении оболочек вращения и круглых пластин — метод разделения переменных с применением тригонометрического представления функций в окружном направлении [38, 39, 43, 54] и др. С помощью такого подхода двумерные краевые задачи определения НДС пластин и оболочек сводятся к решению краевых задач для систем ОДУ. Порядок систем дифференциальных уравнений, при этом, возрастает кратно, пропорционально таким параметрам как количество учитываемых гармоник в разложениях функций в ряды Фурье, порядок используемых сплайнов и д.п. В ряде практически важных случаев краевая задача для систем ОДУ высокого порядка может распадаться на подзадачи меньшего порядка, что существенно упрощает их численное решение.

При численном интегрированим краевых задач для систем ОДУ, возникающих либо после понижения размерности исходной задачи, либо при определении одномерного НДС пластин и оболочек, наибольшее распространение получил метод дискретной ортогонализации [8, 33, 38, 39, 43, 49] и др. Среди других методов решения одномерных краевых задач выделим метод сплайн–коллокации [45, 76]. Особо отметим метод инвариантного погружения [3, 5], разработанный специально для решения задач определения НДС многослойных пластин и оболочек в рамках уточненной теории [5].

Нелинейные задачи обычно линеаризуются на разных этапах решения исходной задачи. Для этого используются различные методы линеаризации [42], модификации метода Ньютона [48, 73], метод последовательных нагружений [63] и др.

1.1 О сведении двумерных краевых задач к одномерным Рассмотрим двумерную краевую задачу для системы дифференциальных уравнений с частными производными, разрешенной относительно производных по одному из направлений вида:

где к немому индексу m применяется правило суммирования. К системам такого вида можно привести, например, исходные соотношения определяющие НДС композитных оболочек. При этом порядок системы уравнений и то, какие из производных по второму направлению будут входить в систему, зависит от используемой теории оболочек. Так для классической теории Кирхгофа — Лява порядок системы уравнений равняется 8–ми и в правой части присутствуют производные первого и второго порядков; для теории Тимошенко порядок системы равен 10–ти, а в правой части присутствуют производные только первого порядка; для теории Андреева — Немировского порядок разрешенной системы дифференциальных уравнений равен 12; для теории "ломаной линии" Григолюка — Куликова — 4K + 6, где K — число слоев оболочки.

Рассмотрим несколько практически важных классов краевых задач расчета замкнутых в окружном направлении оболочек вращения и круговых пластин. Первый — класс осесимметричных задач, когда граничные условия и вектор f не зависят от окружной координаты, а все компоненты матриц Am, m > 0 равны нулю. В таком случае решение краевой задачи также не зависит от окружной координаты: y = y(s), где s — координата, отсчитываемая по меридиану.

Другой важный класс задач — линейные задачи неосесимметричной деформации оболочек осесимметричного строения, что соответствует системам уравнений вида:

где — координата в окружном направлении. Граничные условия по окружной координате — условия "склейки", условия по меридиональной координате запишем в виде:

где sl s1 s2 · · · sS sr ; sl, sr — границы области решения задачи, S — порядок системы уравнений. Для этого класса краевых задач опишем метод разделения переменных для сведения двумерной краевой задачи к одномерным.

Итак, если участвующие в постановке краевой задачи (1.2), (1.3) вектор–функция свободных элементов системы f (s, ) и функция gp () допускают разложение в ряды Фурье вида:

то обобщенное решение соответствующей краевой задачи также можно представить в виде ряда Фурье:

Такое представление решения автоматически удовлетворяет условиям "склейки" по окружной координате.

Подставляя разложения (1.4), (1.5) в систему уравнений (1.2) и граничные условия (1.3) получаем:

Домножая скалярно в пространстве L2 выражения (1.6), (1.7) на элементы тригонометрического базиса sin(n) и cos(n) получаем ряд одномерных краевых задач для определения амплитудных факторов:

dyn (s) = (1)m n2m A2m (s)yn (s) + (1)m+1 n2m+1 A2m+1 (s)yn (s)+ dyn (s) = (1)m n2m+1 yn (s)A2m+1 (s) + (1)m n2m yn (s)A2m (s)+ Отметим, что если существует Am (s) = 0, где m — нечетное, то краевые задачи (1.9) и (1.10) для определения yn (s) и yn (s) являются связанными, и, тем самым, порядок разрешающей системы ОДУ удваивается.

Вследствие линейности исходной двумерной и разрешающих одномерных краевых задач, при использовании изложенного подхода необходимо рассчитывать лишь те гармоники, для которых отличны от нуля коэффициенты f n (s), f n (s) или gp,n. Это существенно сокращает объем необходимых для решения задачи вычислений.

1.2 Особенности систем дифференциальных уравнений при решении краевых задач К решению одномерных краевых задач сводится широкий класс задач математической физики и других областей науки, в частности — задачи расчета и анализа НДС композитных оболочечных конструкций. При этом, получаемые системы уравнений могут иметь высокий порядок, содержать большие и малые параметры, вследствие чего решения задач имеют ярко выраженный характер погранслоев. Проблемы численного интегрирования задач Коши для жестких систем ОДУ разработаны достаточно подробно, в отличие от вопросов решения краевых задач. Рассмотрим основные отличия систем уравнений в краевых задачах от систем уравнений в задачах Коши.

В работах, посвященных проблемам численного интегрирования задач Коши для жестких систем ОДУ [46, 60, 65] и др., рассматриваются задачи, решение которых устойчиво относительно начальных данных.

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

Так на рис. 1.1а представлены области устойчивости методов типа Рунге — Кутты второго (пунктирная линия), третьего (штриховая линия) и Рунге — Кутты — Мерсона четвертого (сплошная линия) порядка [60].

В отличие от задач Коши, класс краевых задач с устойчивыми решениями включает в себя задачи, спектры матриц систем которых содержат и отрицательные, и положительные действительные значения, в том числе очень большие. К таким задачам относятся, в частности, задачи статики тонкостенных оболочечных конструкций. Матрицы систем дифференциальных уравнений для краевых задач, возникающих при расчете НДС композитных оболочек вращения на основе классической теории Кирхгофа — Лява или теорий типа Тимошенко, содержат положительные и отрицательные собственные числа с величинами порядка (обозначены цифрами 2 на рис. 1.1б) и, кроме того, очень малые по модулю числа, в том числе ноль. В случае использования уточненных теорий [5, 33] величины действительных собственных чисел достигают значений на два порядка более высоких — от 100 до 1000 и выше. На рис. 1.1 б цифрами 1 обозначены ненулевые действительные собственные значения матрицы системы уравнений, описывающей изгиб длинной трехслойной прямоугольной пластины [5]. Для сравнения приведена область устойчивости метода Рунге — Кутты — Мерсона (обозначена цифрой 3).

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

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

1.3 Методы решения краевых задач для систем обыкновенных дифференциальных уравнений Методы решения одномерных краевых задач можно условно разделить на два больших класса: методы, сводящие краевую задачу к одной или ряду задач Коши и методы, сводящие дифференциальную краевую задачу непосредственно к системе алгебраических уравнений. К первому типу относятся, например, методы начальных параметров [7, 8], ортогональной прогонки [1, 15], дискретной ортогонализации [16], инвариантного погружения [5]. Ко второму классу относятся конечно–разностные [8, 41] и др., конечно–элементные методы [67, 68, 69, 71] и др., методы коллокации, включая метод сплайн–коллокации [45, 76].

1.3.1 Метод начальных параметров Изложим простой способ сведения решения одномерных краевых задач к решению задач Коши. Рассмотрим линейную краевую задачу для системы ОДУ вида где y = y(x) — искомая вектор-функция порядка S; A(x) — матрица системы порядка S S; f (x) — вектор свободных элементов системы порядка S; xl и xr — левая и правая границы области решения задачи;

Gl, Gr — матрицы граничных условий порядка Sl S и Sr S ранга Sl и Sr соответствено, а gl, gr — свободные векторы граничных условий порядка Sl и Sr, Sl + Sr = S.

Общее решение системы (1.11) представимо в виде где Y(x) — матричная функция, вектор-столбцы yi (x), i = 1,..., S которой образуют фундаментальную систему решений (ФСР) однородной части системы (1.11) а y0 (x) — частное решение системы (1.11); c — вектор произвольных постоянных.

В силу единственности решения задачи Коши yi (x) и y0 (x) являются решениями соответствующих задач Коши с начальными данными yi (xl ), y0 (xl ). Кроме того, любое невырожденное линейное преобразование L столбцов набора Y(x) не меняет его статуса фундаментальной системы решений (ФСР) системы (1.15). Из вышесказанного, в частности, следует, что набор {yi (xl )} должен быть линейно независим, т.е. образовывать базис S-мерного векторного пространства, при том, что это может быть произвольный базис.

Таким образом, решение краевой задачи (1.11)–(1.13) представимо в виде (1.14), где Y(x) — матричная функция, составленная из векторстолбцов yi (x), которые являются решениями задач Коши для системы (1.15) с начальными данными yi (xl ) = yi,l, где {yi,l } — произвольный базис S-мерного векторного пространства; y0 (x) — решение задачи Коши для системы (1.11) с произвольными начальными данными y0 (xl ) = y0,l.

Вектор произвольных постоянных c определяется из граничных условий (1.12), (1.13) как решение системы линейных алгебраических уравнений (СЛАУ) Набор вектор-столбцов Y разбиваем на два поднабора Y1, Y2, вектор c на две части c1 и c2, где Y1 — первые Sl столбцов Y, Y2 — оставшиеся Sr столбцов, c1 — первые Sl элементов c, c2 — остальные Sr элементов. Теперь подсистему (1.16) можно записать в виде где Y1,l = Y1 (xl ), Y2,l = Y2 (xl ). Т.к. Yl = Y(xl ) выбирается как произвольный линейно-независимый набор вектор-столбцов, то его можно построить так, чтобы Y2,l Gl, Y1,l = Gl. Если при этом выбрать y0,l : Gl · y0,l = gl, что допустимо, то система (1.18) приводится к виду Gl · Gl · c1 = 0. Вследствие невырожденности Gl · Gl, единственным решением такой системы становится c1 = 0.

Вектор–функции из построенного таким образом набора Y1 не имеют вклада в решение задачи (1.11)–(1.13) и их можно не искать, а неизвестная часть c2 вектора c определяется из граничного условия на правом крае как решение системы Такой подход позволяет при решении краевой задачи (1.11)–(1.13) не искать полную ФСР системы (1.15), ограничиваясь лишь Sr независимыми функциями, построенными особым образом. Т.к. не имеет принципиального значения какую из границ принять за "левый", а какую за "правый" край, общее количество решаемых задач Коши всегда будет не больше [S/2] + 1.

Алгоритм, основанный на изложенном подходе — построение специальным образом начальных данных Y2,l, y0,l, решение (Sr + 1) задач Коши, отыскание неизвестной части c2 вектора свободных коэффициентов, получение решения краевой задачи (1.11)–(1.13) в виде называют методом начальных параметров или методом стрельбы. Он является одним из самых простых методов решения линейных краевых задач для систем ОДУ.

При решении методом начальных параметров краевых задач для систем ОДУ с малыми и большими параметрами возникает две проблемы. Одна из них — описанная ранее проблема неограниченного роста погрешности численного интегрирования задач Коши, в случае, когда в спектре матрицы системы есть значения с большой положительной действительной частью. Другая проблема заключается в том, что в некоторый момент векторы-столбцы матрицы Y могут стать линейно зависимыми. Начиная с такого момента, во всех последующих точках, наборы векторов Y также будут вырожденными, в итоге система (1.19) становится либо неразрешимой, либо имеет бесконечное количество решений, что недопустимо.

По изложенным причинам метод начальных параметров практически неприменим к решению задач расчета НДС тонкостенных композитных конструкций даже при использовании классической теории Кирхгофа — Лява и теорий типа Тимошенко.

В дальнейшем индекс 2 в обозначениях Y2, c2 будет опускаться 1.3.2 Метод дискретной ортогонализации Метод дискретной ортогонализации [16] предназначен для решения двухточечных краевых задач для систем линейных обыкновенных дифференциальных уравнений и является модификацией метода начальных параметров. Он разработан в начале 60-х годов, практически в одно время с двумя другими, схожими методами решения краевых задач для систем ОДУ: методом факторизации, разработанном И.М. Гельфандом и О.В. Локуциевским первоначально для уравнения второго порядка [15], и методом ортогональной прогонки А.А. Абрамова [1]. Все три метода, как и метод начальных параметров, сводят решение краевой задачи к решению задач Коши. Однако, в последних двух методах, строится новая система ОДУ, где искомой является матричная функция и размерность задачи Коши при этом увеличивается кратно. В методе дискретной ортогонализации, так же как и в методе начальных параметров, увеличивается количество решаемых задач Коши при сохранении самой системы уравнений.

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

Изложим основы алгоритма метода дискретной ортогонализации. В любой точке и на любом интервале представление (1.20) можно преобразовать следующим образом где Y(x) = Y(x) · L, = L1 · (c + l), y0 (x) = y0 (x) Y(x) · l; L — произвольное невырожденное линейное преобразование порядка Sr, а l — произвольный вектор порядка Sr. Так как вид представлений (1.20) и (1.21) совпадает, эти преобразования можно продолжать рекурсивно.

Интервал решения разбивается на J подинтервалов узлами xj (xl, xr ), j = 1... J 1, xj xj+1. Начальные данные Yl, y0,l строятся аналогично методу начальных параметров. При достижении заданной точки xj во время численного интегрирования задач Коши набор Y(xj ), y0 (xj ) преобразуется по формуле (1.21), что будет соответствовать применению преобразования на интервале [xj, xr ]. В качестве L выбирается оператор Wj, ортогонализующий и нормирующий Y(xj ), а вектор l = Wj · wj : y0 (xj ) Y(xj ). Для построения Wj и wj используется модифицированная процедура ортогонализации Грамма-Шмидта.

Таким образом, на каждом подынтервале [xj1, xj ], будем иметь где Y(x) = Y(x) · j1 Wi, = cj1, y0 (x) = y0 (x) j1 Y(x) ( k=1 Wk ) · wi, а cj1 — результат применения рекуррентных соотношений В результате, на последнем подынтервале получаем для решения задачи представление вида (1.22) при j = J, которое аналогично (1.20), что позволяет, как и в методе начальных параметров, вычислить неизвестный вектор cJ1 разрешив систему вида (1.19). Получая оставшиеся неизвестные векторы cj по рекуррентной формуле, обратной (1.23), избавляемся от всех неизвестных в правых частях представлений (1.22) для всех подынтервалов [xj1, xj ], построив, тем самым, решение краевой задачи (1.11)–(1.13) на всей области [xl, xr ].

Опишем конкретный вид процедур, применяемых в алгоритме МДО:

процедуры ортогонализации Грамма — Шмидта (см., например, [8]) с необходимыми модификациями и процедуры численного интегрирования задачи Коши для системы ОДУ. Процедура ортогонализации набора векторов Y, y0 и построения матрицы W и вектора w для соответствующего линейного преобразования является одной из основных процедур, применяемых в методе. Пусть на входе имеется набор векторов Y, y0, необходимо построить ортонормированный набор Y = Y · W и ортогональный ему вектор y0 = y0 Y · W · w.

В качестве базисного вектора выбирается первый столбец y1 матрицы Y, он нормируется y1 = y1 / y1, где норма порождена обычным скалярным произведением: y = (y, y), это и будет первый столбец матрицы Y. Первый элемент первого столбца W1 матрицы W при этом равен W1,1 = 1/ y1, остальные равны нулю Wj,1 = 0, j = 2,..., Sr.

Далее i-ый столбец yi, i = 2,..., Sr матрицы Y строится следующим образом. Вектор yi раскладывается на проекцию yi = j=1 wj yj и ортогональное дополнение yi = yi yi к подпространству, натянутому на уже построенные вектора yj, j = 1,..., i 1, здесь wj = (yi, yj ). Теi1) Построенная таким образом матрица W будет иметь верхнетреугольный вид с ненулевыми диагональными элементами, если же в некоторый момент yi = 0, то это означает, что набор Y линейно зависим и процедура останавливается, извещая об ошибке. В случае успешного выполнения процедуры ортогонализации можно определить вектор w = Y · y0. При численной реализации условие равенства нулю нормы вектора yi необходимо заменить приближенным из-за существования и накопления ошибок округления.

Другая важная процедура — численное интегрирование задач Коши для систем ОДУ вида (1.11) и (1.15) на интервале [xj1, xj ], j = 1,..., J.

Литература, посвященная проблемам численного интегрирования задач Коши, в том числе для жестких систем, достаточно обширна [7, 46, 60, 65] и др. В ней описано множество подходов и методов — неявных, явных, типа Рунге — Кутты, типа Розенброка и т.д. Среди требований, предъявляемых к процедуре численного интегрирования, выделим достаточную устойчивость и простоту реализации при высокой точности. Явные методы уступают неявным методам и методам типа Розенброка в вопросе обеспечения устойчивости счета, однако существенно превосходят их по простоте реализации. Среди явных схем выделим одношаговые методы типа Рунге — Кутты, так как применение многошаговых схем к задачам с ярко выраженными локальными и краевыми эффектами может приводить к сглаживанию и обрезанию "пиков" решения. Среди методов Рунге — Кутты наилучшим образом зарекомендовала себя пятистадийная схема четвертого порядка точности, построенная Мерсоном, обладающая достаточно широкой областью устойчивости и высоким порядком сходимости при относительной простоте реализации [60].

Итак, если необходимо решить задачу Коши для системы ОДУ вида то решение y(x1 ) строим с помощью процедуры Рунге — Кутты — Мерсона 4-го порядка по формуле где 1.3.3 Метод сплайн-коллокации Запишем нелинейную краевую задачу для системы обыкновенных дифференциальных уравнений в виде:

Sp ; xp [xl, xr ], xp+1 >= xp.

Для решения задачи (1.27), область определения разбивается на J подинтервалов точками j (j = 0,..., J) так, чтобы {xp } {j }. На каждом подинтервале [j1, j ] выбираются Q + 1 точек коллокации jq (j = 1,..., J, j0 = j1, jQ = j ). Решение y(x) представляется непрерывной кусочно-многочленной функцией u(x) порядка Q+1, точно удовлетворяющей системе (1.27) в точках коллокации и граничным условиям, что соответствует системе Система (1.28) в общем случае нелинейна. Для ее линеаризации используется метод Ньютона — Рафсона. Решение ищется в виде u(x) = lim u(k) (x), где u(k) (x) — кусочно-многочленная функция, являющаяся решением линейной задачи Здесь vi (x) = f /ui (x, u(k1) (x)), ui (x) — i-ая компонента вектор-функции u(x).

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

Здесь V = (Q + 1)J + 1 — размерность пространства; i — вектор коэффициентов размерности S; Bsi — i-й B-сплайн. После подстановки (1.30) в (1.29) получается линейная система алгебраических уравнений для определения коэффициентов i :

где Матрица системы уравнений (1.31) почти блочно-диагональная. Для ее решения применяется метод исключения Гаусса с частичным выбором главного элемента.

Метод сплайн-коллокации реализован в пакете прикладных программ COLSYS [70, 76, 81].

Глава Алгоритм решения краевых задач для систем обыкновенных дифференциальных уравнений Как уже отмечалось, метод дискретной ортогонализации получил широкое распространение при решении задач теории пластин и оболочек. Он хорошо зарекомендовал себя и показал высокую эффективность при расчете осесимметричного НДС резервуаров и сосудов высокого давления [17]–[19]. Использование при этом классической теории оболочек Кирхгофа — Лява позволило применять метод на равномерных сетках и без особых модификаций. Однако, переход к использованию уточненных теорий, учитывающих деформации поперечного сдвига на основе гипотез о нелинейном распределении вектора перемещений по толщине оболочки, заставил искать пути по адаптации алгоритма МДО к решению краевых задач для систем ОДУ, в спектрах матрицы Якоби которых содержатся сильно удаленные от мнимой оси величины из обеих комплексных полуплоскостей. Отметим, что в работе [16] представлен общий алгоритм МДО, но отсутствует описание ряда необходимых при его реализации процедур. Тем самым, у исследователей имеется определенная свобода в реализации алгоритма метода. Среди известных модификаций МДО отметим модификации с сохранением промежуточной информации о структуре решения (векторы y0 и матрицы Y, W) лишь в суммарном виде в точках выдачи результатов [35] и без сохранения этой информации с решением на обратном ходу алгоритма задачи Коши [6]. Отметим, что последняя модификация неприменима к решению задач, у которых в спектрах матриц Якоби систем уравнений содержатся большие отрицательные собственные числа.

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

2.1 Проблемы вычисления векторов начальных данных и решения многоточечных задач Начальные данные Yl, y0,l для задач Коши строятся на первом этапе алгоритма метода и должны удовлетворять подсистеме (1.16), что, по сути, эквивалентно отысканию общего решения системы При этом набор векторов Yl является общим решением однородной части системы (2.1), а вектор y0,l — частным решением последней. Соответственно, таких решений можно построить бесконечно много. Однако, с учетом специфики решаемой далее задачи, в соответствии с излагаемым методом, желательно выполнить следующие условия: во-первых, набор Yl должен быть ортонормированным, а во-вторых, вектор y0,l должен быть ортогонален набору Yl.

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

Суть алгоритма состоит в том, что матрица Yl строится как ортогональное дополнение к матрице Gl путем разложения элементов произвольного базиса пространства на проекции и ортогональные дополнения, после чего находится и частное решение y0,l системы (2.1). Удобно априори считать, что набор векторов-строк Gl невырожден и ортонормирован.

Добиться ортонормированности невырожденного набора можно применением описанной выше процедуры ортогонализации:

тогда выражения (2.1) можно переписать в виде Выражения (2.1) и (2.2) эквивалентны вследствие невырожденности W и их общий вид совпадает, что позволяет положить заранее, что требование ортонормированности выполняется.

Для того, чтобы построить ортогональное дополнение к набору векторов Gl выбирается произвольный базис E = {ei }, i = 1,..., S векторного пространства (достаточно использовать обычный базис из единичных ортов). Последовательно все ei раскладываются в сумму проекции i e и ортогонального дополнения i к подпространству, натянутому на векe тора набора Gl и j, j < i. Ненулевые i, дополнительно нормированные = i / i, образуют искомое ортогональное дополнение Yl. Отметим, ei e e что условие неравенства нулю векторов i, эквивалентное отличию от нуля их нормы i, при численной реализации нужно заменить приe ближенным.

Для получения последнего искомого начального вектора y0,l необходимо найти решение системы Так как Gl — ортонормирован, то выбирая y0,l = Gl · gl и подставляя его в систему (2.3) получаем тождественное равенство при этом вектор y0,l будет ортогонален Yl по построению.

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

Рассмотрим многоточечную краевую задачу для системы (1.11) с граничными условиями вида:

где xp xp+1, сумма рангов матриц Gp равна размерности задачи S, а ранг каждой такой матрицы равен Sp. Для удобства будем полагать что наборы векторов-строк Gp ортонормированы, чего можно добиться описанным выше способом.

В качестве начальных векторов–решений yi, i > 0 в первом узле интервала решения задачи xl выбираем произвольный базис S–мерного векторного пространства, например, состоящий из единичных ортов, а начальный вектор для неоднородной задачи Коши выбираем нулевым y0 = 0. Далее решаем построенные задачи Коши, аналогично алгоритмам методов дискретной ортогонализации и начальных параметров. При достижении точки xp имеем: вектор y0 и набор векторов–столбцов Y, которые являются решением соответствующих задач Коши. Количество столбцов матрицы Y, при этом, должно равняться Sp = p >=p Sp. Требуется построить набор из Sp Sp вектор-столбцов Y и вектор y0 такие, что для любого вектора c размерности Sp Sp выполняется тождество:

выполнив, как и в случае двухточечной краевой задачи, условия ортонормированности набора Y и ортогональности к нему вектора y0. Кроме того, необходимо построить зависимость c = c(c ), где размерность вектора c равна Sp и для любого вектора c выполняется тождество:

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

Итак, пусть существует такое, что выполняется (2.7), тогда для любого вектора c имеем:

где · c = c |c. Подставляя (2.8) в (2.6) получаем:

Так как по построению Y Gp и так как Gp — ортонормированный набор векторов, то выражение принимает вид:

При этом из (2.5) получаем y0 = Gp · gp. Таким образом определяются y0 и связь между c и c и остается описать алгоритм построения Y и.

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

Опишем соответствующий алгоритм по шагам:

1) задаются начальные значения i = 1, J = 2) определяются проекция yi и ортогональное дополнение yi вектора yi на пространство натянутое на вектор-столбцы Gj, j = 1,..., Sp матрицы Gp и построенные векторы yj, j = 1,..., J:

при этом элементы матрицы находящиеся на пересечении i-го столбца и j-ой (j = 1,..., Sp + J) строки равны соответствующим коэффициентам в разложении:

3) если норма yi не равна нулю, то переходим к пункту 4) полагаем элементы матрицы равными увеличиваем i на единицу и переходим к пункту 5) полагаем вектор yJ+1 и элементы матрицы равными увеличиваем i и J на единицу 6) если i = Sp + 1 и J = Sp Sp + 1, то процедура завершена успешно, иначе если i = Sp + 1 либо J = Sp Sp + 1, то процедура завершена неуспешно, иначе переходим к пункту 2.

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

В узле xp|p=P данный обобщенный алгоритм восстанавливает решение задачи, которое будет равно вектору y0, и если xP < xr, то дальнейшее решение задачи эквивалентно решению задачи Коши. В процессе обратного хода при прохождении через точки xp восстанавливается вектор c(c ) и по соответствующим формулам вычисляется решение, аналогично алгоритмам методов дискретной ортогонализации и начальных параметров.

2.2 Обеспечение устойчивости расчетов Для решения методом дискретной ортогонализации краевых задач, возникающих при определении напряженно–деформированного состояния оболочек вращения в рамках классической теории оболочек Кирхгофа — Лява либо теории типа Тимошенко, достаточно использовать равномерные сетки узлов ортогонализации и интервалов интегрирования. При этом типичное количество узлов ортогонализации, достаточное для устойчивого счета, относительно невелико — до 20, а количество интервалов интегрирования, необходимое для достижение приемлемой относительной точности (порядка 105 ) составляет от 200 до 400. Однако при переходе к использованию уточненных теорий пластин и оболочек, учитывающих поперечные сдвиги, таких значений оказывается явно недостаточно, и появляется необходимость в дополнительном контроле точности и устойчивости расчетов и в их обеспечении.

В методе дискретной ортогонализации, в дополнение к шагу интегрирования, появляется ещё один управляющий параметр — расстояние между узлами ортогонализации. Совместно, эти параметры позволяют контролировать и управлять устойчивостью численного расчета, его точностью и объемом вычислений. Дополняя алгоритм метода дискретной ортогонализации процедурами выбора шага интегрирования и расстояния между узлами ортогонализации, можно построить алгоритм решения краевых задач для систем ОДУ с малыми и большими параметрами, наподобии тех, что возникают при определении НДС пластин и оболочек в рамках уточненных теорий [5, 33].

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

Для этого рассмотрим задачу изгиба слоистой длинной прямоугольной пластины, нагруженной равномерно распределенным давлением. Обезразмеренная разрешающая система ОДУ, описывающая поведение такой пластины в уточненной постановке, имеет вид [5]:

Здесь где h, l — полная толщина и длина пластины; hk, Ek, k — толщина, модуль Юнга и коэффициент Пуассона k-го слоя пластины; K — количество слоев; W,, U — безразмерные обобщенные перемещения.

Выбирая в качестве искомой вектор-функции y(x) = W, W, W, W,,, U, U, представляем систему (2.12) в виде (1.11). Ненулевые коэффициенты матрицы A и вектора f будут при этом равны где c1 = a2 a3 a1 a5, c2 = a1 a6 a2, c3 = a2 a1 a4. Разрешающая система уравнений является системой с постоянными коэффициентами. Шесть собственных значений матрицы системы 1 = 2 = 3 = 4 = 5 = 6 = 0, оставшиеся два 7 = 8 =, где = A65 + A46 A64 R [5]. При этом, величина пропорциональна отношению l/h.

С учетом введенных обозначений, общее решение системы (2.12) имеет вид где b = (24C f4 )/A46, C = f4 A65 /(242 ), C1,..., C8 — произвольные постоянные. Полученные аналитические решения, наряду с полиномиальными, содержат экспоненциальные компоненты вида exp((x 1)) (кривая 1 на рис. 2.1а), exp(x) (кривая 2 на рис. 2.1а) (в случае, приведенном на рис. 2.1а, = 155.66).

Жесткое защемление пластины соответствует граничным условиям вида В этом случае произвольные постоянные будут равны 3 A65 /(6A46 A64 ).

Рассмотрим трехслойную пластину, внешние слои которой имеют толщину h1 и выполнены из материала с модулем Юнга E1 и коэффициентом Пуассона 1. Внутренний слой, толщины h2, выполнен из материала с модулем Юнга E2 и коэффициентом Пуассона 2. Из рис. 2.1б видно, что спектральный радиус матрицы системы существенно зависит от механически характеристик материалов и структуры пластины, и эта зависимость нелинейная.

На рис. 2.2 приведен вид решения задачи изгиба жестко защемленной на обоих краях трехслойной пластины в разрешающих функциях, нормированных в равномерной метрике (геометрические и механические параметры пластины: l/h = 15, h2 /h1 = 8, E1 /E2 = 10, 1 = 2 = 0.3).

Сплошные линии на рис. 2.2 соответствуют W (x), U (x); штриховые — W (x), U (x); пунктирные — W (x), (x); штрих-пунктирные — W,. Видно, что получаемые решения имеют ярко выраженные краевые эффекты.

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

Исследовалось влияние расстояния между узлами ортогонализации и количества интервалов интегрирования между ними на точность расчетов, которая использовалась как критерий их устойчивости. Проведенные численные эксперименты показали, что с увеличением жесткости задачи начинает проявляться проблема "размазанности" "области перехода" от неустойчивых расчетов к устойчивым при уменьшении расстояния между узлами ортогонализации. При этом, данная проблема проявляется сильнее при большем числе интервалов интегрирования между соседними узлами ортогонализации. Так на рис. 2.3 представлены зависимости погрешности вычислений (логарифмированной по основанию 10) от числа узлов ортогонализации J для различного числа интервалов интегрирующей процедуры: штриховые кривые соответствуют трем интервалам, пунктирные — четырем, сплошные — пяти. Зависимости получены для описанной выше задачи при значении l/h = (рис. 2.3а) и l/h = 500 (рис. 2.3б), спектральные радиусы матриц системы при этом равняются = 1037.73 и 5188.63 соответственно. Из рисунка видно, что при меньшем спектральном радиусе прослеживается четкая граница перехода от неустойчивых решений к устанавливающимся, т.е. таким, погрешность которых стабильно уменьшается при уменьшении шага сетки. При больших спектральных радиусах очень существенной становится проблема сдерживания экспоненциального роста погрешности численного интегрирования, что проявляется в возникновении осциляций погрешности на графиках, соответствующих большему числу интервалов интегрирования между ортогонализациями. Отметим, что при значении l/h = 1000 осциляции возникают уже при использовании трех интервалов интегрирования. Таким образом, чем больше уравнений, тем меньшее количество интегрирований допускается между узлами ортогонализации.

Рассмотрим теперь какое расстояние между узлами ортогонализации необходимо для устойчивого расчета. Для этого построим зависимость числа узлов ортогонализации J при котором расчет является устойчивым (по признаку стабильно небольшой погрешности) от спектрального радиуса системы. На рис. 2.4а представлены два варианта таа б ких зависимостей: при нефиксированном (больше трех, подбираемом в зависимости от результатов) количестве интегрирований между узлами ортогонализации (пунктирная кривая) и при фиксированном (три интегрирования, сплошная кривая). Во втором случае фиксировался и критерий устойчивости — достижение погрешности вычислений порядка 103. Видно, что использование числа интервалов интегрирования больше трех целесообразно при относительно небольших, иначе это приводит к нестабильности в поведении метода, о чем сказано выше, а выигрыш в количестве ортогонализаций при этом составляет 0, где — коэффициент Пуассона, общий для всех слоев, E (k) — модуль упругости k-го слоя, k = d,..., d, hd1 = h/2, hd,..., hd = h/2 — координаты внешних поверхностей и границ раздела слоев, z — координата, отсчитываемая по нормали к срединной поверхности пластины.

Получить разрешающую систему уравнений можно двумя путями:

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

и переписать исходную систему уравнений в виде:

Из (2.24) следует:

Подставляя это выражение в (2.25), с учетом (2.23), выводим:

Подобным образом получается и разрешающая система для функций (2.22):

Разрешающие системы ОДУ (2.28) и (2.23), (2.26), (2.27) имеют 6–ой порядок и могут быть записаны в виде:

Для разрешающего вектора y = 0, 1, 2, w0, w1, w2 ненулевые компоненты матрицы системы A и вектора свободных элементов f равны:

где приняты обозначения:

В этом случае два ненулевых собственных значения матрицы A равны Если в качестве разрешающего вектора выбрать y = 0, 1, 2, w0, w1, w2, то ненулевые компоненты матрицы A и вектора f равны:

В этом случае собственные значения матрицы A равны 1 = 0, 2 = Таким образом, получены две разрешающие системы ОДУ с различной структурой матрицы A — спектр матрицы с компонентами (2.30) состоит из четырех нулей и двух взаимнообратных действительных значений, не зависящих от координаты r. В спектре матрицы A с компонентами (2.31) содержатся положительные и отрицательные действительные величины, стремящиеся к бесконечности пропорционально r1 при r 0. При этом общие решения данных систем дифференциальных уравнений имеют схожие структуры, отличающиеся лишь множителями вида r и r1.

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

либо При q3 = P = const общее решение построенной системы ОДУ, выписанное в функциях w0, w1, 0, имеет вид:

где Im, Km, — модифицированные функции Бесселя первого и второго рода, а Ci, i = 1,..., 6 — свободные коэффициенты. Соответственно, решения исходной системы (2.20) и системы (2.28) выражаются через w0, w1, 0 следующим образом:

Системы алгебраических уравнений для определения свободных коэффициентов Ci, i = 1,..., 6 получаются подстановкой общего решения в соотношения граничных условий. Разрешая их и подставляя Ci в общие решения, получаем аналитическое решение соответствующих краевых задач.

Рассмотрим вид решения задачи изгиба трехслойных жестко защемленных кольцевых пластин. На рис. 2.5 представлен вид функций w,, 0 для нагруженных равномерным давлением P = 106 Па пластин с различными внутренними радиусами R0 = 0.5, 0.1, 0.01 м, внешним радиусом R1 = 10 м, толщинами слоев h0 = 0.6, h1 = 0.2 м, модулями упругости материалов слоев E0 = 40, E1 = 200 ГПа. Наибольший интерес представляет вид компонент решения, 0 вследствие наличия модифицированных функций Бесселя в структуре их аналитического представления, из–за чего само решение имеет ярко выраженный характер погранслоев. Из рис. 2.5 б видно, что в решении краевой задачи для системы (2.23), (2.26), (2.27) (задача № 1, пунктирные кривые) краевые эффекты больше проявляются на внешней кромке пластины, а в решении краевой задачи для системы (2.28) (задача № 2, сплошные кривые) — на внутренней и, при этом, величина краевого эффекта пропорциональна На рис. 2.6 можно видеть влияние на вид функции других параметров задачи — соотношений h1 /h0 (рис. 2.6 а) и E1 /E0 (рис. 2.6 б).

Цифрам на рис. 2.6 а соответствуют соотношения толщин слоев 0.5h0 = 9h1 (кривая 1), h0 = 8h1 (кривая 2), h0 = 3h1 (кривая 3), 3h0 = 4h1 (кривая 4), при полной толщине пластины h = 1.0 м. Остальные параметры пластины равны: R0 = 0.1, R1 = 10 м, E0 = 20, E1 = 200 ГПа. Видно, что соотношение толщин слоев существенно влияет на амплитуды деформаций поперечного сдвига, при этом, чем меньше толщина внешних, более жестких слоев, по отношению к толщине внутреннего слоя, тем более ярко выраженными становятся краевые эффекты на кромках пластины.

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

"1" для E1 = 100E0, "2" для E1 = 10E0, "3" для E1 = 2E0, "4" для E1 = E0, "5" для 2E1 = E0, во всех случаях принято E1 = 200. Другие параметры пластины при этом равны: R0 = 0.1, R1 = 10 м, h0 = 0.6, h1 = 0.2 м. При уменьшении соотношения E1 /E0 величины экстремумов поперечных сдвигов вблизи кромок пластины ведут себя немонотонно:

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

Аналитические решения, имеющие ярко выраженный характер погранслоев, переменные коэффициенты систем дифференциальных уравнений и наличие в них параметров, управляющих величиной и характером краевых эффектов, позволяют исследовать эффективность разработанного алгоритма численного решения одномерных краевых задач при решении описанной задачи изгиба кольцевой пластины. Два варианта разрешенных систем уравнений, имеющих отличные друг от друга структуры матриц Якоби и, как следствие, аналитических решений, удобно использовать для оценки работы алгоритма при решении задач, в которых краевые эффекты более ярко выражены либо на левой границе интервала решения (задача № 2), либо на правой (задача № 1).

В табл. 2.3 представлены результаты работы алгоритма при решении задачи изгиба трехслойных кольцевых пластин разного радиуса, жестко защемленных на краях. Радиус внутреннего отверстия при этом равен R0 = 0.1; толщины и модули упругости материалов слоев: h1 = 0.2, h0 = 0.6 м, E1 = 200, E0 = 40 ГПа. В таблице приведены максимальные приведенные погрешности, получаемые при определении, и w, получаемые при определении w; общее число вызовов процедуры численного интегрирования задачи Коши J и число интервалов в итоговой сетке интегрирующей процедуры J ; число узлов ортогонализации J. Видно, что алгоритм успешно справляется с решением обеих краевых задач — получаемые погрешности не превосходят 0.1% кроме одного случая. Узлы ортогонализации расположены достаточно равномерно по интервалу решения задачи везде, кроме области прилегающей к левой границе, поэтому их число растет пропорционально длине интервала, и, при этом, совпадает для обоих вариантов задачи.

При изменении радиуса внутреннего отверстия можно видеть, как проявляется наличие в спектре матрицы системы дифференциальных уравнений задачи № 2 величин, обратно пропорциональных r. В табл. 2. представлены результаты работы алгоритма при решении задач с разными значениями параметра R0. При этом R1 = 10 м, а остальные параметры задачи идентичны предыдущему случаю. Из таблицы видно, что при решении задачи № 2 может возникнуть необходимость в дополнительных узлах ортогонализации вблизи левой границы. Однако, вследствие того, что область, в которой собственные числа матрицы Якоби много больше единицы, относительно узка, то необходимости сильно мельчить сетку узлов ортогонализации не возникает.

Для обеих задач в спектре матриц Якоби их систем ОДУ содержатся величины, пропорциональные h1 при h 0. При этом для задачи № слагаемые, вносящие такой вклад в собственные значения, не зависят от r, поэтому, когда h мало, спектральные радиусы матрицы Якоби много больше единицы на всем интервале решения задачи. В табл. 2.5 приведеТаблица 2. ны результаты применения разработанного алгоритма решения краевых задач для систем ОДУ к задаче изгиба трехслойных кольцевых пластин различной толщины h с параметрами: R0 = 0.1, R1 = 10 м, E0 = 40, E1 = 200 ГПа, h0 = 3h1.

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

после интегрирования задач Коши на всех участках [xj, xj+1 ]. Конечно, можно фиксировать единую относительную погрешность для всех задач Коши, тогда итоговая относительная погрешность должна будет иметь тот же порядок. Однако такой подход может привести к неоправданному измельчению сетки интегрирующей процедуры и, соответственно, росту количества вызовов последней. Альтернативный вариант заключается в выполнении предварительного расчета с выполнением лишь условий устойчивости и дальнейшее уточнение решения при необходимости.

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

На рис. 2.7 представлен вид погрешности для функции, получаемой при решении описанной выше задачи изгиба длинной слоистой пластины методом дискретной ортогонализации на равномерной сетке интегрирующей процедуры, содержащей 300 элементов (а), а также, для сравнения, погрешности, получаемой при использовании пакета COLSYS (б). Наблюдается "взрывное" накопление погрешностей вблизи границ, которое можно попробовать нивелировать или хотя бы сгладить. Для этого построим неравномерную сетку на основе экспоненциальных функций вида (1+exp(l (x0l x)))1, (1+exp(r (x0r x))1, где l = r = 40, x0l = 0.12, x0r = 0.88 (параметры сетки подобраны экспериментальным путем). Вид шага (x) такой сетки, содержащей 300 элементов, приведен нии этой сетки погрешности для функции. В табл. 2.6, 2.7 приведены относительные погрешности, полученные на равномерных сетках и сетках, построенных по вышеприведенному принципу соответственно.

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

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

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

При использовании процедуры Рунге — Кутты — Мерсона для численного интегрирования задач Коши, можно воспользоваться следующей формулой приближенного определения погрешности вычислений:

При построении неравномерной сетки воспользуемся принципом эквираспределения [82], выбирая в качестве управляющей функцию, зависящую от погрешности:

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

Рассмотрим применение изложенного подхода. Для начала оценим точность формулы (2.35). На рис. 2.9 представлены точная (сплошные кривые) и приближенно вычисленная (штриховые линии) погрешности, полученные на задаче изгиба длинной слоистой пластины при l/h = 15, для последовательностей сеток из 50, 100 элементов (слева) и 200, элементов (справа). Видно, что во-втором случае погрешность определяется более чем удовлетворительно, как по характеру, так и по значениям.

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

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

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

На рис. 2.10, 2.11 представлены величины шагов сеток интегрирующей процедуры и вид получаемых при их использовании погрешностей для тестовых задач с параметрами l/h = 15 и l/h = 150 соответственно.

Сплошные кривые соответствуют формуле (2.36), пунктирные — формуле (2.37), штриховые — формуле (2.38). Из рис. 2.10б, 2.11б видно, что при использовании формулы (2.38) погрешность в области краевых эффектов "размазывается" лучше, чем при использовании формулы (2.36).

Однако формула (2.36) строит более разреженную сетку в области интервала решения где краевые эффекты отсутствуют, что совершенно не увеличивает максимальные полученные погрешности, однако позволяет существенно сократить вычислительные затраты.

Для того, чтобы сравнить эффективность использования строящихся сеток между собой и с сетками равномерными, приведем данные о количестве элементов в каждой из сеток, необходимого для достижения заданной погрешности (табл. 2.8 и 2.9 для тех же задач, соответственно). Использование формулы (2.36) позволяет сократить количество выТаблица 2. Тип сетки количество интервалов интегрирования Тип сетки количество интервалов интегрирования числений по сравнению с равномерной сеткой от трех до двадцати раз в зависимости от задачи и требуемой точности, формулы (2.38) — от трех до девяти раз. При этом, чем выше требуемая точность, чем больше спектральный радиус матрицы системы уравнений в задаче, тем эффективнее использование этих формул. Формула (2.37) показывает себя неэффективно во всех случаях.

Таким образом можно сказать, что применение неравномерных сеток, строящихся по формулам (2.36), (2.38) позволяет в несколько раз сокращать время расчета, необходимое для достижения заданного по точности результата. Однако наиболее эффективным будет применение их сочетания: формулы (2.38) в областях краевых (локальных) эффектов и формулы (2.36) в областях, где краевые эффекты выражены слабо, либо отсутствуют. Конечно, при построении соответствующих алгоритмов необходимо учитывать "цену" построения таких сеток, т.е. время, которое на это затрачивается.

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

Классическая теория оболочек основана на гипотезах Кирхгофа — Лява сохранения и недеформируемости нормального к отсчетной поверхности элемента (см. [39, 61] и др.). Использование достаточно простых соотношений классической теории Кирхгофа — Лява позволяет в ряде практических случаев получить удовлетворительные результаты. Однако повышение требований к прочности и надежности современных оболочечных конструкций приводит к необходимости рассмотрения теорий, основанных на менее жестких предположениях, чем гипотеза сохранения нормального элемента.

Большинство уточненных, по сравнению с теорией Кирхгофа — Лява, теорий, описывающих НДС оболочек, получены с использованием гипотез о характере распределения напряжений, деформаций и перемещений по толщине оболочки. Преимуществом такого подхода является относительная простота разрешающих соотношений. Однако такой подход не обладает возможностью уточнения полученных на его основе результатов.

При построении уточненных теорий основным фактором является учет деформаций сдвига во всех или отдельных слоях оболочки. Гипотеза о прямолинейном элементе для всего пакета в целом стала одной из самых распространенных, используемая при построении теории оболочек типа Тимошенко (см. [5, 39] и др.). Хотя для многослойных конструкций с существенно различными механическими параметрами слоев принятие гипотезы прямой линии для всего пакета в целом может вносить существенную погрешность в получаемые результаты.

Использование гипотезы прямой линии для каждого слоя в отдельности (см., например, [33]) позволяет находить более точные решения, но приводит к уравнениям, порядок которых зависит от количества слоев, что затрудняет получение конкретных результатов. Другой способ построения уточненных моделей заключается в задании нелинейного закона распределения поперечных напряжений по толщине. Значительный интерес представляет теория, разработанная в [4, 5], позволяющая рассчитывать НДС многослойных анизотропных оболочек с учетом поперечного сдвига в каждом слое. Порядок полученной системы уравнений при этом не зависит от количества слоев.

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

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

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

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

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

3.1 Моделирование свойств и критерии прочности полиармированных композитов К настоящему времени разработано большое число структурных моделей композитов. Первые модели, представляющие собой простое правило смеси, были предложены Фойгтом для осреднения матрицы жесткости и Рейссом для осреднения матрицы податливости. Далее этот подход развивали Д.С. Аболиньш, В.В. Болотин, Ф.Я. Булавс, А.К. Малмейстер, Ю.В. Немировский, А.Л. Рабинович, Дж. Сендецки, А.М. Скудра, В.П. Тамуж, Ю.М. Тарнопольский, Г.А. Тетерс и др., применявшие гипотезы Фойгта и Рейсса в различных комбинациях относительно направления армирования. Гипотеза об однородном НДС в структурных компонентах — основное предположение для этих моделей.

Более сложные модели упругого поведения композиционных материалов, учитывающие неоднородность полей структурных напряжений и деформаций, построены в работах Н.С. Бахвалова, Г.А. Ванина, Т. Ишикавы, Р. Кристенсена, С. Кобаяши, Б.Е. Победри, Ю.В. Соколкина, З. Хашина, Р. Хилла, Л.П. Хорошуна, Т.Д. Шермергора, С. Штрикмана и др.

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

Связь между осредненными напряжениями, 3 и деформациями e, 3 в полиармированном слое (обобщенный закон Гука) имеет вид:

где — приращение температуры. Соотношения (3.1) называют также соотношениями термоупругости или, при отсутствии температурных слагаемых, просто соотношениями упругости.

3.1.1 Модель В.В. Болотина Структурная модель армированного волокнами композита, описанная в [9] (МБ), стала основанием большого количества исследований в данном направлении и широко используется при моделировании композитных конструкций. Модель основывается на предположениях об однородности НДС в изотропных упругих волокнах и во всем объеме изотропной идеально упругой матрицы, совместности деформаций волокон и связующего в направлении армирования, равенства напряжений в волокнах и матрице по остальным компонентам.

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

где отброшены члены, имеющие порядок квадрата от коэффициентов Пуассона. Здесь E1, E2 — эффективные модули упругости в направлении армирования и поперек его, G — эффективный модуль сдвига, — эффективный коэффициент Пуассона в плоскости слоя; E,, с индексами «c» и «a» — модули упругости, коэффициенты Пуассона и объемные содержания связующего и армирующих волокон, соответственно, при этом c + a = 1.

Из соображений симметрии тензора податливости имеем:

Формулы для эффективных коэффициентов теплового расширения имеют вид:

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

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

Невыписанные величины получаются по правилу симметрии, либо полагаются равными нулю. Здесь и далее полагаем, = 1, 2 и =.

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

a33 = az Gc + z где z и az = 1 z — удельные толщины слоев содержащих арматуру и прослоек связующего, l1 () = cos(), l2 () = sin(), а p() — плотность распределения волокон по углам, нормированная на отрезке [0, ].

Плотность является обобщенной функцией : она может содержать особенности типа дельта-функции.

На основе исходных предположений модели в работе [2] для плоского напряженного состояния получены формулы, связывающие осредненные напряжения и деформации с напряжениями и деформациями в компонентах армированного слоя: в волокнах арматуры (добавочный индекс «a») и в связующем (добавочный индекс «c1») Величины, соответствуют системе координат, связанной с направлением армирования в рассматриваемом слое, и определяются по формулам:

Там же выписана более полная формула для эффективного модуля E2, в которой учитываются слагаемые порядка квадрата от коэффициентов Пуассона:

Наряду с описанной моделью рассмотрим базирующуюся на ней модель армированного слоя [33] (МБГ), которая позиционируется авторами для описания свойств резинокордных материалов. В качестве величин характеризующих удельное содержания компонент принимаются диаметр D круглого волокна и плотность укладки — характерное количество нитей n на ширину l сечения представительного элемента плоскостью, перпендикулярной к направлению укладки волокон. Описанные характеристики позволяют вычислить объемное содержание арматуры в слое:

где h — толщина слоя. Для определения эффективного модуля E2 используется формула (3.12). Для учета сдвигов в слое и поперечных сдвигов используются формулы: (3.4) для A1212 G12 и A1313 G13, (3.2) для A2323 G23. В данной модели нет описания температурных компонент. При восстановлении напряжений в волокнах их напряженнодеформированное состояния рассматривается как одномерное:

а напряжения в связующем не восстанавливаются. Здесь e11 = e11 l1 () + 2e12 l1 ()l2 () + e22 l2 () — продольные деформации в волокне.

3.1.2 Модели Ю.В. Немировского Модели КМ с одномерными и двумерными волокнами являются развитием идей заложенных в МБ. Пренебрежение “многомерностью” НДС волокна и рассмотрение его как нити позволило построить существенно более простые, с точки зрения математики, модели КМ с одномерными волокнами (МОВ, МОВУ) [55, 56]. Соотношения этих моделей, в отличие от МБ, линейны относительно таких характеристик КМ как объемные содержания связующего и армирующих волокон и имеют простой вид, поэтому их удобно применять при решении задач рационального проектирования композитных конструкций. Модель с двумерными волокнами (МДВ) [57] является, в свою очередь, уточнением МБ с более четкими и точными математическими формулировками и выводами формул для характеристик армированного материала.

Ниже будут выписаны соотношения для структурной модели с двумерными волокнами. Уравнения для перечисленных выше моделей с одномерными волокнами, а также нитяной модели (НМ) [62] вытекают из них как частные случаи. Материал армированного слоя оболочки полагается идеально упругим, с плоскостью симметрии, параллельной отсчетной поверхности оболочки. Приведем основные предположения для рассматриваемых моделей [57].

1) Полиармированный слой (рис. 3.1) представляет собой изотропное упругое однородное связующее с внедренной в него регулярной сетью однонаправленных упругих изотропных волокон.

2) Число армирующих волокон достаточно велико, так что армированный слой можно считать квазиоднородным.

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

4) В процессе деформирования удлинения и сдвиги остаются малыми по сравнению с единицей.

5) Элементы композиции соединены идеально: отсутствует проскальзывание между армирующими элементами и связующим.

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

7) (a) Армированный слой представляет собой сеть одномерных волокон; нагрузка воспринимается только армирующими волокнами, а напряжения в связующем не учитываются (рис. 3.2 a).

(b) Армированный слой представляет собой сеть одномерных волокон; поперечные сдвиговые напряжения воспринимаются только слоями связующего, а армированные слои являются абсолютно жесткими на сдвиг (рис. 3.2 b,c).

(c) Армирующие волокна имеют прямоугольное поперечное сечение (рис. 3.2 d), а их НДС рассматривается как двумерное.

Предположения 1 — 6 общие для всех рассматриваемых моделей.

Предположение 7a соответствует нитяной модели, предположения 7b, 7c — моделям с одномерными и с двумерными волокнами соответственно.

Компоненты тензоров эффективных жесткостей и эффективных тепловых жесткостей выпишем в виде [5, 57]:

ЗдесьEa, — модуль Юнга и коэффициент Пуассона материала n-го семейства арматуры; n, z — интенсивности армирования в поверхn ности и в направлении толщины оболочки для n-го семейства арматуры;

c = 1 n для МДВ и МОВУ, а для МОВ и НМ c = 0; n — угол укладки n-го семейства арматуры.

Соотношения МДВ [57] получаются при = 1; моделям с одномерными волокнами [56] соответствует = 0; формулы для НМ [62] получаем при = 0, Ec = 0.

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

Для нитяной модели напряжения и деформации арматуры n-го слоя будут иметь вид [62]:

Для модели с одномерными волокнами к соотношениям (3.16) добавляются выражения для напряжений в прослойках связующего [58]:

Для уточненной модели с одномерными волокнами к напряжениям и деформациям, определяемым соотношениями (3.16), (3.17) нужно добавить выражения для напряжений и деформаций связующего в армирующем слое:

Невыписанные компоненты напряжений для волокон и связующего для МОВ и НМ полагаются равными нулю.

Для МДВ напряжения в арматуре, в связующем армирующего слоя и в связующем между армирующими слоями будут иметь вид [5]:

Здесь, = 1, 2 и используется правило суммирования по немым индексам.

3.1.3 Критерии прочности и начального разрушения композиционных материалов Проблемы прочности КМ разрабатывались многими авторами и получили в литературе широкое освещение. Так же как и при построении моделей КМ, при установлении критериев прочности можно выделить два основных похода — феноменологический и структурный. В рамках первого из них КМ рассматривается как квазиоднородная упругая среда, для которой постулируется критерий прочности. Параметры, входящие в его математическую формулировку, определяются из экспериментальных данных. Среди феноменологических критериев прочности важное место занимает тензорно–полиномиальный критерий [53], который обобщает практически все известные феноменологические критерии. Однако следует отметить, что даже для относительно простых видов напряженного состояния требуется реализовать весьма трудоемкие программы экспериментов и математической обработки полученных данных. Другой недостаток таких критериев — их формулировка в терминах средних напряжений, что не позволяет выявить механизм возникновения начального разрушения и исследовать направление его дальнейшего развития.