1. Постановка задачи.
В инженерной практике часто возникает необходимость в решении плоских и пространственных задач теории упругости для тел сложной геометрии и тел, состоящих из материалов с упругими различными характеристиками. Эти задачи позволяют определять напряженно-деформированное состояние тел не только в рамках упругой постановки, но и получать решение упруг опластиче-ских задач, так как многие численные методы их решения сводятся к последовательности упругих задач. К числу таких задач относятся, например, задачи о распределении напряжений вблизи выработок в сплошных горных породах, о давлении на подземные сооружения со стороны окружающего массива, о напряженно-деформированном состоянии гидротехнических сооружений.
Проблема концентрации напряжений в деформируемом твердом теле с различного рода концентраторами напряжений актуальна, поскольку концентрация напряжений является одним из основных факторов, определяющих прочность и долговечность материалов и элементов конструкций независимо от того, вызвана она неоднородностью структуры материала или является следствием наличия геометрических особенностей в рассматриваемых областях.
Пространственные задачи теории упругости для неодносвязных тел относятся к одному из наиболее сложных и недостаточно изученных разделов механики деформируемого тела. Вопросы исследования этих задач связаны с необходимостью построения точных различных моделей упругой среды, учитывающих ее механические и геометрические неоднородности. Последнее обстоятельство в наибольшей степени относится к теории разрушения, а также к механике композитных материалов. С другой стороны, необходимость решения задач этого класса диктуется большой практической важностью исследования концентраций напряжений в окрестности трещин, полостей, включений при расчетах на прочность, при проектировании сооружений.
Характерным отличием для всех упомянутых проблем является сложная геометрия и многосвязность областей, в которых решаются краевые задачи, и, как следствие этого, существенные математические трудности при их решении. В настоящее время наиболее хорошо изучены различные типы краевых задач теории упругости для односвязных тел достаточно простой геометрии.
2. Состояние проблемы.
Методы решений задач теории упругости для многосвязных областей можно разделить на аналитические, численно-аналитические и численные. Широкий круг задач о концентрации напряжений в упругих двумерных телах решается с использованием методов теории функций комплексного переменного. Для плоских областей сложной конфигурации, занятых упругой средой, наиболее эффективным методом является метод Н. И. Мусхелишвили, позволяющий получить аналитическое решение основных задач теории упругости для случаев, когда известны функции, осуществляющие конформное отображение единичного круга или полуплоскости на рассматриваемую область [1−7]. В работе [8] отмечается, что построение этих отображений является наиболее трудоемким моментом при решении задач теории упругости методом Мусхелишвили. Авторы работ [9,10] на основе многократных конформных отображений выводят функции напряжений многих комплексных переменных для многосвязных областей. В итоге задача сводится к бесконечной системе линейных уравнений, решаемой численно.
Среди работ, посвященных изучению напряженного состояния многосвязных областей, важное место занимает метод интегральных уравнений (метод комплексного потенциала). Класс задач об упругом равновесии кольца или плоскости с круговым отверстием решен Галиным Л. А., Лурье И. И., Шерманом с применением метода потенциала. В работах [11−18] линейная задача теории упругости для многосвязных областей и для областей с границей сложной формы сводится к системе сингулярных интегральных уравнений.
Для задач расчета концентраций напряжений в многосвязных областях с успехом применяется метод Фурье разделения переменных, с помощью которого в ряде случаев удается получить точные решения в тригонометрических рядах [19,20].
В.Л. Рвачевым был разработан новый метод (метод R-функций) для решения задач теории упругости в областях сложной конфигурации, позволяющий строить координатные последовательности для областей практически произвольной геометрии и краевых условий сложного вида. Применение аппарата R-функций [21,22] позволяет получить решение в виде формулы, содержащей явную зависимость от геометрических и физических параметров задачи.
Успехи методов теории функций комплексного переменного в двумерных задачах вызвали попытки использовать комплексные функции при решении пространственных задач теории упругости. Теория кватернионных функций развилась в качестве многомерного аналога комплексного анализа. Получены кватернионные представления общего решения уравнений Ламе в виде пространственного аналога формулы Колосова-Мусхелишвили [23]. В работе [24] предложен метод решения ряда краевых задач теории упругости для шара и пространства с шаровой полостью с использованием кватернионных функций.
К числу наиболее популярных и эффективных методов решения краевых пространственных задач для тел сложной формы относятся методы потенциала, в частности, методы граничных интегральных уравнений. Теоретическую основу этих методов составляет разработанная Купрадзе В. Д. и Михлиным теория обобщенных упругих потенциалов и многомерных сингулярных интегральных уравнений. Теория потенциала по своей общности занимает то же место для пространственных задач, что и теория функций комплексного переменного для плоских задач теории упругости. Но в частных задачах для тел сложной формы, применение теории потенциала затруднительно. В работах [25−29] рассматриваются основные краевые задачи теории упругости для полупространства, ослабленного полостью, причем поверхность полости предполагается принадлежащей классу поверхностей Ляпунова. Включения в форме сферы, эллипсоида, цилиндра, параллелепипеда в основную среду рассмотрены в работах [30−34], а в [35] геометрические конфигурации достаточно произвольной формы конструируются из набора типовых фрагментов канонических поверхностей с применением аффинных преобразований и инверсий.
Метод Фурье с успехом применяется для решения задач теории упругости для тел со сферическими, эллипсоидальными, цилиндрическими включениями [36,39], обобщается на задачи для тел, ограниченных одновременно несколькими координатными поверхностями цилиндрической, сферической и сфероидальной системы координат [40,42]. В общей постановке метод Фурье обоснован для областей типа полупространства, цилиндра, тора, сфероида, гиперболоида, параболоида и конуса [43].
Для решения осесимметричных пространственных задач применяется метод декомпозиции путем установления зависимостей между осесимметричным и плоским состояниями. Это делает возможным исследование пространственной задачи с помощью аппарата ТФКП. Такую же связь находят авторы работ [45,46] и тем самым сводят расчет пространственных конструкций к расчету двумерных систем.
При решении задач пространственной теории упругости для тел сложной конфигурации применяются также вариационные принципы Лагранжа и Кас-тильяно. Использование их позволяет построить устойчивую разностную схему после введения равномерной сетки и получения разностного оператора вариационным методом [47,50].
Метод возмущения формы, основанный на разложении по малому параметру потенциала простого слоя, распределенного по двусвязной области [51], распространен и на тела неканонической формы [52,53 ].
В последние годы получает распространение метод геометрического погружения, основная идея которого заключается в сведении исходной краевой задачи для тела сложной геометрии к итерационной последовательности краевых задач для канонической области [54]. Численная реализация метода геометрического погружения возможна на основе вариационных методов, методов конечных и граничных элементов. В качестве канонической области может быть выбрано все евклидово пространство [55], или ненарушенный горными работами породный массив в случае расчета напряженно-деформированного состояния горного массива в окрестности участков выработанного пространства.
Подавляющее большинство задач для плоских и трехмерных областей решается с помощью численных методов, основанных на дискретизации расчетной области, например, методом конечных разностей [58,59]. Непосредственное использование метода конечных разностей для областей сложной геометрии осложняется ухудшением сходимости в случае существенного отличия формы тела от канонической. Наибольшее распространение в практике расчетов для областей сложной геометрии получил метод конечных элементов (МКЭ). Простота основных идей, гибкость, наглядность, доступность математического аппарата для широкого класса пользователей сделали этот метод очень популярным.
В настоящее время наиболее распространенной является схема МКЭ, основанная на аппроксимации смещений с использованием принципа минимума потенциальной энергии. Поскольку деформации в такой постановке определяются после численного дифференцирования поля перемещений, то напряжения вычисляются с меньшей точностью по сравнению с перемещениями. В этой связи есть разработки МКЭ в напряжениях [60].
Применение МКЭ при расчете пространственных задач, даже при малой степени дискретизации расчетной области на трехмерные конечные элементы, приводит к системам линейных алгебраических уравнений, решение которых в классе точных (прямых) методов затруднительно. Причиной этого служит высокий порядок и ширина ленты разрешающей системы уравнений. В самом деле, время решения системы известными прямыми методами растет пропорционально порядку матрицы жесткости и квадрату ширины ее ленты. В этой связи разрабатываются многосеточные методы конечных элементов [61−65], которые являются полуитерационными методами решения пространственной и плоской задачи теории упругости. Они сводятся к следующей последовательности операций: невязка проектируется на крупную (вспомогательную) сетку конечных элементов и система решается в классе точных (прямых) методов, затем полученная поправка интерполируется на мелкую (исходную) сетку конечных элементов и добавляется к начальному приближению решения системы. Полученное решение рассматривается как новое приближение, и итерационный процесс продолжается до достижения требуемой точности.
Обычно при решении пространственных задач теории упругости применяется вариационная постановка, связанная с интегрированием по объему области. Для численной реализации в этом случае применяется МКЭ. Если вариационную постановку связать с интегрированием по поверхности области, то для численной реализации применяется метод граничных элементов (МГЭ). Он является основным конкурентом МКЭ и быстро развивается в последние 10 лет [66−70]. Снижение на единицу размерности задачи и, соответственно, снижение затрат памяти в МГЭ особенно важно в пространственных задачах. Фактором, сдерживающим распространение МГЭ, является сложность математического аппарата (необходимость вычисления сингулярных интегралов), решение несимметричных линейных систем с плотно заполненными матрицами.
В работах [71,72], в качестве альтернативы обычному МКЭ, предложен новый подход к решению задач анализа напряженно-деформируемого состояния массивных пространственных конструкций, названный методом конечных объемов. Он основан на дискретизации деформируемого трехмерного тела в виде ансамбля конечных объемов произвольной топологии с помещением характеристических узлов в центрах масс.
В последние годы получили развитие различные модификации МКЭ. Например:
• в [73] предложен метод внешних конечноэлементных аппроксимаций,.
• для дискретизации случайных полей разрабатывается стохастический метод конечных элементов [74],.
• для описания траектории распространения трещины и разрушения предложены сингулярные конечные элементы [75],.
• для расчета напряженно-деформируемого состояния очень длинных выработок в породном массиве рассматриваются «бесконечные» конечные элементы [76],.
• для конкретных задач теории упругости разрабатываются специальные конечные элементы [77−79].
При анализе сложных пространственных: тел методом конечных элементов приходится рассматривать расчетные модели с несколькими тысячами степеней свободы. Это требует большой скорости счета, а современные процессоры имеют максимальное быстродействие до трех млрд. операций в секунду, поэтому в последние годы появилась методика использования параллельных процессоров [80−83]. Метод применим к конечноэлементным моделям с большой шириной ленты матрицы жесткости, в противном случае эффективность вычислений с использованием параллельных процессоров очень низкая [84].
Одним из наиболее трудоемких этапов решения задач МКЭ является дискретизация расчетной области, т. е. построение сетки конечных элементов в области решения. В связи с этим МКЭ применяется обычно для решения задач или с достаточно простой геометрией области, где строятся конечные элементы, форма которых обусловлена геометрией аппроксимируемой области [8592], или с небольшим количеством узлов и элементов, допускающих подготовку исходных данных вручную. В некоторых работах, например, при расчете напряженно-деформируемого состояния гравитационной плотины [93], с целью избежать нерегулярности сетки конечных элементов, применяется принцип стандартной области. Стандартная область имеет простую форму, покрыта регулярной конечноэлементной сеткой и содержит в себе действительную область, занятую телом. Цель алгоритма состоит в том, чтобы освободить от усилий требуемый контур, что эквивалентно близости к нулю поверхностных напряжений на нем. Трудности дискретизации области существенно осложняют применение МКЭ для решения практически важных задач механики деформируемого твердого тела, в частности, задач механики горных пород, отличающихся сложной геометрией областей решения. Поэтому высок интерес к разработке алгоритмов автоматического построения конечноэлементной сетки.
Алгоритмы компьютерного автоматического построения конечноэле-ментных сеток разделены на полуавтоматические и полностью автоматизированные [94]. В полуавтоматических методах КЗ сетка формируется из треугольников или прямоугольников для двумерной задачи и в виде четырехгранников и шестигранников для трехмерных задач. Если область решаемой задачи имеет сложную конфигурацию, то она предварительно разбивается на простые подобласти. Например, в работе [95] при формировании дискретной модели тела предварительно вручную выделялись составные блоки, т. е. крупные элементы, конфигурация которых позволяет запрограммировать процесс разбивки на конечные элементы, или в [96] производится декомпозиция исходной полигональной области на выпуклые многоугольники, каждый из которых разбивается на четырехугольные КЭ.
Алгоритмы полностью автоматизированного построения КЭ сеток, включая метод Делоне, не требуют предварительного разбиения рассматриваемой области сложной конфигурации, но связаны с вводом достаточно громоздкой дополнительной информации об ее границах [97,99]. Например, в работе [100] предложен эффективный метод для дискретизации большого класса деформируемых твердых тел на топологически простые подобласти. Используются геометрические свойства анализируемых тел, конфигурации их срединных поверхностейподбор подобластей осуществляется из некоторого числа предложенных определенных типов поверхностей. Таким образом формируется КЭ сетка из шестигранных элементов для тел с выпуклыми границами.
В работе [101] описан алгоритм для автоматизированного построения трехмерной КЭ сетки, применимой для конструкций невыпуклой геометрии. Реализация алгоритма начинается с размещения редких узлов на границах геометрической модели и соответствующей трехмерной триангуляцией по методу Делоне. Затем вводятся инкрементальным образом последующие узлы с проверкой геометрии сетки из тетраэдров и топологической совместимости между триангуляцией Делоне и геометрической моделью.
Чтобы улучшить точность конечно-элементного анализа, предложены схемы адаптивного автоматического построения КЭ сетки с применением оценки получающейся погрешности расчета. Поскольку компоненты напряжений должны быть совместны в каждом узле сетки, их несовместность может рассматриваться как критерий ошибки конечно-элементного анализа. Элемент, характеризующийся большой несовместностью напряжений в узловых точках, автоматически делится на более мелкие элементы с введением дополнительных узлов [102,103]. В работе [104] предлагается другой метод адаптивного перепостроения сеток для трехмерного КЭ анализа напряженно-деформируемого состояния массивных пространственных тел. В первом приближении генерируется конечно-элементная модель, следующая начальной плотности распределения узлов по заданию пользователя. Затем реализуется конечно-элементный анализ и производится оценка погрешности расчета. Определяется оптимальная плотность узлов и сетка перестраивается. Процесс повторяется до тех пор, пока точность решения не будет удовлетворять критерию, назначенному пользователем.
Во всех рассмотренных работах конечно-элементная сетка формируется так, чтобы отвечать реальной конфигурации области. В работе [105] при КЭ расчете концентраций напряжений в упругих телах около круглого отверстия, был введен неоднородный конечный элемент, внутри которого проходила граница круглого отверстия, и матрица жесткости такого элемента вычислялась на основе численного интегрирования. Но проверка попадания узла интегрирования в выпуклую область с аналитически заданной границей, в частности, в круг, не представляет трудностей.
В заключение можно сделать вывод: для численного исследования напряженно-деформированного состояния многосвязных тел применяется конечно-элементная сетка, учитывающая геометрические особенности исследуемой области. Построение такой сетки является одним из наиболее трудоемких этапов решения задач методом конечных элементов.
3. Краткая характеристика работы.
Целью данной работы является исследование напряженно-деформированного состояния многосвязных тел сложной геометрии в двумерной и трехмерной постановке. Задачи, которые были решены для достижения поставленной цели, заключаются в следующем.
• Разработаны схемы расчета, основанные на методе конечных элементов в перемещениях, с дискретизацией исследуемой области на прямоугольники (в случае плоской задачи) или на прямоугольные параллелепипеды (в случае пространственной задачи).
• Разработан комплекс программ на языке СИ, реализующих предложенную методику с использованием неоднородных конечных элементов.
• Проведено сравнение численных расчетов с аналитическим решением для двумерных и трехмерных задач теории упругости.
• Проведен двумерный конечно-элементный анализ напряженно-деформированного состояния системы плотина-основание Братской ГЭС в упругой постановке и с учетом нелинейных деформаций бетона. Выполнена оценка длительной прочности сооружения при воздействии периодических нагрузок.
• Построена имитационная модель усталостного разрушения гребня бетонной плотины при воздействии периодических нагрузок.
Структура диссертационной работы в целом соответствует последовательному решению перечисленных выше задач.
В первой главе определяются основные соотношения метода конечных элементов при дискретизации области на прямоугольные КЭ с использованием неоднородных конечных элементов. Описывается процедура проверки попадания в невыпуклую плоскую область узла интегрирования. Для задачи Кирша выполнено сравнение результатов расчета выбранным методом и аналитическим.
Во второй главе проводится расчет напряженно-деформированного состояния глухой секции плотины Братской ГЭС по предложенной методике, и сравниваются полученные результаты с результатами, полученными по общепринятой методике. Сделан расчет напряженно-деформированного состояния секции с полостью расширенного шва в рамках упругой модели. Построены упругопластическая и вероятностная модели разрушения гребня бетонной плотины от периодических силовых воздействий, вызванных движением железнодорожных составов. Исследована устойчивость системы плотина-основание при периодических нагрузках, вызванных колебанием уровня воды в водохранилище.
В третьей главе определяются основные соотношения МКЭ для пространственных задач теории упругости. Рассматриваются тела невыпуклой геометрии с дискретизацией области на прямоугольные параллелепипеды с использованием неоднородных конечных элементов (НКЭ). Описана процедура проверки попадания точки в трехмерную невыпуклую область. На двух примерах показывается сходимость предложенного подхода.
В заключении сформулированы основные результаты и выводы диссертации.
Основные положения и результаты диссертации были доложены: • На научно-технической конференции БрГТУ, г. Братск, 1997 г., 1999 г., 2000 г.
• На 3-й региональной научно-практической конференции «Интеллектуальные и материальные ресурсы Сибири», г. Иркутск, 1999 г.
• На международной научно-практической конференции «Математические методы в образовании, науке и промышленности», г. Тирасполь, 1999 г.
• На всероссийском семинаре «Энергетика: экология, надежность, безопасность», г. Томск, 1999 г.
• На международной конференции «Геодинамика и напряженное состояние земных недр», г. Новосибирск, 1999 г.
Основное содержание диссертации опубликовано в работах [140−155].
Основные результаты и выводы состоят в следующем:
1. Предложен подход к дискретизации области расчета при конечноэлемент-ном моделировании с введением неоднородного конечного элемента. Рассмотрены двумерные и трехмерные области сложной конфигурации Получены соотношения метода конечных элементов для регулярной прямоугольной сетки, не зависящей от геометрических особенностей объекта, в случае плоской задачи теории упругости. Для пространственных объектов сложной конфигурации получены соотношения МКЭ с применением сетки из прямоугольных параллелепипедов.
2. Построен алгоритм численного решения, реализующий предложенную методику с сеткой прямоугольных элементов, и с сеткой из параллелепипедов. Показано удовлетворительное совпадение результатов расчета, выполненного численно, с аналитическим решением. Программа расчета полей напряжений плоской двусвязной области с применением НКЭ зарегистрирована в Реестре программ для ЭВМ (свидетельство № 2 000 610 780 от 22 августа 2000 г.).
3. Исследовано напряженно-деформированное состояние секции плотины с расширенным швом. Расчеты показывают, что при нормальном подпорном уровне в 112 м в теле плотины появляются зоны опасных растягивающих напряжений: Максимальные растягивающие напряжения до 1,2 МПа наблюдаются на внутренней стороне расширенного шва, примыкающей к верховой грани, на высоте около 50 м от подножия плотины.
4. Исследована устойчивость верхней части плотины Братской ГЭС при периодических нагрузках, вызванных прохождением по гребню плотины железнодорожных составов. Построена модель последовательной квазистатической смены состояний объекта: появление зон пластичности, их развитие и разрушение. Исследована кинетика разрушения, показано, что разрушение может прекращать дальнейшее развитие в условиях отсутствия дополнительного увеличения нагрузок.
5. Сделана оценка устойчивости гидротехнического сооружения, подвергающегося переменным нагрузкам, связанным с изменением уровня воды в водохранилище. Оценен запас прочности сооружения, находящегося в районе, подверженном сейсмическим воздействиям со стороны Байкальской рифтовой зоны.
6. Построена вероятностная модель усталостного разрушения гребня бетонной плотины. Установлено вероятное число циклов нагружения, вызывающее потерю несущей способности гребня плотины при различных уровнях нагрузки.
7. Разработан комплекс программ на языке СИ для решения задач пластичности, для моделирования появления и распространения зон разрушения от периодических воздействий, для имитационного моделирования усталостного разрушения тел, имеющих геометрические особенности.
ЗАКЛЮЧЕНИЕ
.
В диссертационной работе предложен подход к дискретизации области расчета при конечно-элементном моделировании с введением неоднородного конечного элемента. Рассмотрены двумерные и трехмерные области сложной конфигурации. По предложенной методике исследована устойчивость плотины Братской ГЭС при воздействии периодических нагрузок.
