Метод численного решения краевых задач нелинейной теории упругости Текст научной статьи по специальности «Физика»

УПРАВЛЯЮЩИЕ И ИЗМЕРИТЕЛЬНЫЕ СИСТЕМЫ

УДК 519.3:539.3:624.131

МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ НЕЛИНЕЙНОЙ ТЕОРИИ УПРУГОСТИ

К.С. КУРОЧКА

Большинство методов решения краевых задач классической теории упругости базируется на предположениях физической и геометрической линейности рассматриваемого твёрдого тела. Причём, под геометрической линейностью подразумевается линейный характер зависимости между деформациями и перемещениями, а под физической - линейный характер зависимости между напряжениями и деформациями. Однако большинство реальных материалов имеют нелинейную природу, обладают рядом свойств: анизотропия, дилатация и другие. Использование линейных моделей для таких материалов приводит к значительным погрешностям [1,2].

Учёт одновременно геометрической и физической нелинейности является сложной математической задачей [1, 2]. Поэтому в дальнейшем будем анализировать только физическую нелинейность.

Для решения краевых задач линейной теории упругости наиболее целесообразно применять методы математического моделирования на основе теории систем и системного подхода, метод конечных элементов и (или) суперэлементов [2]. Поэтому рассмотрим методы моделирования нелинейного деформирования твёрдых тел с использованием метода конечных элементов.

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

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

Существующие методы решения краевых задач нелинейной теории упругости

Решение линейных задач теории упругости методом конечных элементов сводится к решению системы линейных алгебраических уравнений вида [2]:

Учреждение образования «Гомельский государственный технический университет имени П.О. Сухого», Республика Беларусь

Введение

(1)

• где {к} - вектор узловых перемещений.

Пусть, например, при условии

матрица жёсткости, (г?! - вектор узловых

(М) <0 или (р2 (И})< 0,

(2)

• где {сг} - вектор напряжений; {е} - вектор деформаций, деформирование области до определенного предела определяется линейно-упругим законом:

(3)

• где [Е0 ] - линейный модуль упругости, а при невыполнении условия (2) деформирование определяется нелинейным законом

F (Ы Н)=0. (4)

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

М=ММ-к })+К } (5)

• где [о] - матрица упругости, содержащая механические характеристики

материала; {к0} и (<т0} - начальные деформации и начальные напряжения.

В случае использования уравнения состояния (3) или (5), вектор [Я] в (1) определяется внешними действующими силами и предварительным напряжённо-деформируемым состоянием.

Можно утверждать, что, если удастся найти такое решение уравнения (1), что при соответствующем подборе одного или нескольких входящих в (5) параметров [о], к} или {<г0} это уравнение и соотношение (4) удовлетворяются при одинаковых значениях напряжения и деформации, то полученное решение будет искомым.

Очевидно, что при решении следует руководствоваться некоторыми алгоритмами для определения [о], {к } К }. при этом можно ограничиться вычислением только одной из указанных величин, зафиксировав другие.

Вопрос этого выбора зависит от метода решения линейной задачи и уравнения состояния. В любом случае указанные параметры определяются итерационно, т. е. постепенно уточняются. Если определяется [о]. то получается так называемый метод переменной жесткости. Это название предопределено тем, что [о] определяет жесткость материала. Если же определяются к0} или {<г0}, то имеем так называемые методы

начальных деформаций или начальных напряжений [1,2].

Метод переменной жёсткости [1, 2] используется тогда, когда уравнение состояния (4) известно в явном виде. По аналогии с (5), уравнение состояния для нелинейных деформаций можно записать таким образом:

М=0№1-к Мы}. (6)

Используя (6) для построения уравнения (1), получим:

[к (ЫШ= {Я}. (7)

Начальные напряжения и деформации из (6) вошли в вектор {я} из (7). Система (7) решается итерационно.

1. На основании физических данных определяются упругие характеристики тела и строится система (7). Она будет эквивалентна системе (1).

2. Решается полученная система и определяется вектор {§■}.

3. Для всех конечных элементов, для которых выполняется условие (2) на основании (6), исходя из уравнения состояния (4), вычисляется [£(Ы0)] (вектор {§■}(, - первое приближение решения) и матрица [к({^}0)].

4. Решается система (7) с полученной матрицей жесткости.

Пункты 3 и 4 повторяются до тех пор, пока {§'}г+1 - {§} < {^}, где {^} - вектор допустимой точности решения.

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

Метод начальных напряжений [1,2] применим, если определяющее уравнение (4) разрешимо относительно напряжений:

Ын }= /, ({«})• (8)

• тогда соотношение (5) для упругого материала можно привести к виду:

{ы}= ф({80}), (9)

• задавая соответствующим образом {<г0}. Так как {<г0} влияет на вектор внешних

сил {Я} , то приходим к решению уравнения

[к Ы={Я({?})}. (10)

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

н=[к т, (11)

• где {я}0 соответствует приложенным нагрузкам.

Определяются напряжения {<г0 }1, необходимые для приведения упругого решения в соответствие с реальными напряжениями при достигнутых деформациях. Далее, с учетом начального напряжения {<г0 }1, находится {Я}1 и определяется

Ы1 =[к ]1{я}1 и т. д.

Процесс останавливается аналогично методу переменной жесткости.

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

а{£}1 =[к ]-1 А{Я}1 и т. д.

Итерации продолжаются до тех пор, пока величина А^} не станет достаточно близкой к нулю.

В случае, когда закон нелинейного деформирования задаётся в виде

Ы = Аё”~" к, (12)

• вектор {<г0} остаётся неизвестным. Поэтому возникает задача по известной интенсивности напряжений <у1 найти компоненты вектора {<г0}.

Предположим, что существует некоторое линейное тело, для которого

(<т„ } = Е '[Е ](£'|,

• причём вектор {к'} известен из линейного решения задачи, тогда

е' е'

где {ал } - определяется из линейного закона деформирования,

Г

_ н

{а11 } - определяется из нелинейного закона деформирования. Заметим, что

' 1 1 ,

{ал } = Ео [Е ]{е'}, ±- {а о} = [ Е ]{ е'} , — {ал } = [Е]{е'},

• отсюда

Е '

{ао}= — {О1}. (13)

Ео

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

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

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

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

М = Л ({а}). (14)

Совпадение соотношений (14) и (4) может быть достигнуто при соответствующем выборе {е}. Уравнение

[К о Ь}=№»

• опять решается итерационным методом, но теперь упругие деформации, полученные на каждом шаге, сравниваются с деформациями, соответствующими определяющему соотношению (14), и их разность используется для оценки невязки силы Д{Я},. В дальнейшем процесс идентичен описанному выше, и, в частности, матрица жесткости остается постоянной на любом шаге.

Недостатки существующих методов решения краевых задач нелинейной теории

упругости

При комплексном моделировании сложных пространственных объектов (например, зданий и сооружений), когда система линейных алгебраических уравнений (1) содержит сотню тысяч неизвестных и более и требуются огромные затраты машинного времени на её решение, становится весьма важным, сколько необходимо провести итераций для

нахождения решения. В рассмотренных методах количество итераций составляет порядка 6-8 [3]. Поэтому для практического применения при решении таких задач данные методы будут малопригодными.

Предлагаемый метод простой линеаризации предлагается применять следующим образом:

1. Решаем упругую задачу методом конечных элементов (1).

2. Для каждого конечного элемента находим ({а0},{е0}).

3. Если условие (2) не выполняется, то предположим, что найденное линейное решение отличается от точного на Ае,, т. е.

Сг = 80 + Аег , (15)

• где е, - интенсивность нелинейных деформаций; е0 - интенсивность линейных деформаций.

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

/ е,)=f (8;)+/е >ч +... (16)

Из разложения в ряд Тейлора возьмём первых два члена:

/ (е,) * /е)+ /'(<)Ае,.

Тогда а */Г? Х' +,//°)_ V? , т. е. нелинейный закон деформирования (4)

А С

приближённо заменим выражением

аН * Ае, + С. (17)

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

Предположим, что существует некоторое линейно-деформируемое тело, для которого закон Гука имеет вид:

аг = Е сек ег,

ан 0

• где есек = г , еi - интенсивность деформации, полученная из решения задачи

е г

(1); аН - соответствующая этим деформациям интенсивность напряжений, определённая

из нелинейного закона деформирования по формуле (8). Закон Гука для этого тела можно переписать в виде:

СЛ }= Есек Е ]{е}.

Компоненты вектора {е}, соответствующие е0, известны из решения задачи (1), следовательно из закона Гука можно определить и компоненты вектора а*. В точке е0 значения а л и а^ совпадают. Тогда положим {тл }= {г11} в точке е0.

Рассмотрим (17). Пусть а,л = а,н - С . Предположим, что существует некоторое линейное тело, для которого закон Гука будет ал = Ае, или С* }=[Е ]{е^}. Теперь, вместо (17), исходный нелинейный закон деформирования будем аппроксимировать следующим выражением:

С }=Е‘]{е}+[с*]. (18)

Потребуем, чтобы в точке е0 (18) совпадало с (8). Тогда из (18) определим матрицу [С *]:

[С * ]=Сн }-[£ * ]{е}.

Деформации конечного элемента можно выразить через перемещения

И=№), (19)

• где [5] - координатная матрица конечного элемента.

Подставим (19) в (18):

{а) = [Г ][в]{я}+[с * ]. (20)

Полная энергия деформируемой системы минимальна, если она находится в состоянии статического равновесия:

дП = 0, (21)

3{я}

где

п = 2 М № -{я}' Я. (22)

2 V

Подставив (19) и (20) в (22), интегрируя, получим:

П = 2 {я}' {В}' [е• ][в]{я}г + 2{я}' {В}' [с• ]г - {я}' {Я}, (23)

• где V - объём конечного элемента.

Подставим (23) в (21) и продифференцируем:

/ V}' [е*][V] {я}= {Я}- V{В}' [с*]. (24)

Уравнение (24) - это и есть основное уравнение метода конечных элементов. Решаем это уравнение и для каждого конечного элемента находим ({а}, {е}).

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

Верификация предложенного метода и разработанного программного обеспечения

К

Для верификации рассмотрим моделирование деформаций фрагмента перекрытия, состоящего из двух многопустотных плит ПК 63.15.8АТ800А-8 с отверстием в середине пролёта и ПК 63.15.6АТ800АТ-2 с замоноличенными стыками с верхней связью плит [4].

Многопустотные плиты рассматривались как неоднородные нелинейно-деформируемые объёмные тела, определённые в и-связанной области.

Расчётный пролёт плиты ПК 63.15.8АТ800А-8 L0 = 6,2 м, ширина плиты В = 1,49м, рабочая арматура 6014АТ8ОО. Плита имеет 7 отверстий диаметром 159 мм, защитный слой бетона 20 мм, в середине пролёта выполнено сквозное отверстие с размерами 400^1700 мм для сантехнических устройств, модуль упругости бетона 22,4 МПа.

Расчётный пролёт плиты ПК 63.15.6АТ800АТ-2 L0 = 6,2 м, ширина плиты В = 1,49м,

рабочая арматура 5012АТ800. Плита имеет 7 отверстий диаметром 159 мм, защитный слой бетона 20 мм, модуль упругости бетона 22,5 МПа.

Модуль упругости бетона замоноличивания 29,5 МПа для первого стыка, 32,7 МПа для второго и третьего стыков. Опорные столики плит перекрытия были приведены к закладным деталям в стеновых панелях, а плиты перекрытия соединены между собой в

верхней зоне металлическими пластинами. Загружение проводилось от 0^/2 до

/ м

5740^2 .

/ м

При аппроксимации конечными элементами сквозные горизонтальные отверстия заменялись равновеликими параллелепипедами и исключались из рассмотрения. Цилиндрическая арматура заменялась также равновеликими параллелепипедами. Нелинейный закон деформирования определялся на основании экспериментальных данных [4] в виде (12). Точность решений 5 = 0,001.

Граничные условия задавались следующим образом: по боковым плоскостям запрещены перемещения вдоль оси Y, торцы плиты жёстко закреплены.

При моделировании использовались конечные элементы в форме параллелепипедов [8]. Вдоль оси X каждая многопустотная плита разбивалась на 30 конечных элементов, вдоль оси Y - 20, вдоль оси Z - 15. Найденные значения прогибов сравнивались с экспериментальными данными [4]. Результаты расчётов различными методами приведены в табл. 1 и табл. 2. Вычислительный эксперимент проводился на компьютере с процессором АМО АШ1оп 2500 МГц и ОЗУ 1024 Мб.

Таблица 1

Прогибы плиты ПК 63.15.8АТ800А-8, мм

Нагр узка, / м Линейное решение Переменной жёсткости Начальных напряжений Начальных деформаций Простой линеаризации Экспериментальные данные [4]

Прогиб, мм Время расчёта, с Прогиб, мм й и я а р е т 8 о - л о К Время расчёта, с Прогиб, мм Кол-во итераций Время расчёта, с Прогиб, мм й и ц а р е т 8 о в - л о К Время расчёта, с Прогиб, мм й и ц а р е т 8 о в - л о К Время расчёта, с

1380 0,27 447 0,54 3 1341 0,62 4 1780 0,49 3 1362 0,69 2 902 0,71

2761 0,68 447 0,70 3 1367 0,69 4 1793 0,72 3 1371 0,76 2 897 0,75

4141 0,72 447 0,74 4 1802 0,73 6 2694 0,77 3 1365 0,82 2 888 0,81

5522 1,08 447 2,28 7 3147 2,16 6 2691 2,22 5 2235 2,41 2 904 2,54

6902 1,72 447 3,44 8 3612 3,38 7 3175 3,54 5 2237 3,57 3 1355 3,61

8283 2,33 447 7,19 12 5365 7,15 8 3612 7,31 7 3190 7,44 4 1785 7,69

Таблица 2

Прогибы плиты ПК 63.15.6АТ800АТ-2, мм

Нагр узка, Н/ /м2 Линейное решение Переменной жёсткости Начальных напряжений Начальных деформаций Простой линеаризации Экспериментальные данные [4]

Прогиб, мм Время расчёта, с Прогиб, мм й и я а р е т 8 о - л о К Время расчёта, с Прогиб, мм Кол-во итераций Время расчёта, с Прогиб, мм й и ц а р е т 8 о в - л о К Время расчёта, с Прогиб, мм й и ц а р е т 8 о в - л о К Время расчёта, с

1380 0,12 447 0,14 3 1341 0,09 4 1780 0,03 3 1362 0,05 2 902 0,05

2761 0,20 447 0,24 3 1367 0,27 4 1793 0,18 3 1371 0,22 2 897 0,24

4141 0,37 447 0,41 4 1802 0,39 6 2694 0,47 3 1365 0,43 2 888 0,45

5522 0,98 447 1,28 7 3147 1,76 6 2691 1,72 5 2235 1,75 2 904 1,83

6902 1,68 447 4,12 8 3612 4,08 7 3175 4,28 5 2237 4,25 3 1355 4,29

8283 2,22 447 10,70 12 5365 11,15 8 3612 11,01 7 3190 11,05 4 1785 11,12

Проведенный анализ результатов вычислительного эксперимента

позволяет сделать следующие выводы:

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

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

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

Литература

1. Ильюшин, А.А. Механика сплошной среды: учебник. - 3-е изд. /А.А. Ильюшин. -Москва: Издательство МГУ, 1990. - 310 с.

2. Постнов, В.А. Метод конечных элементов в расчётах судовых конструкций /В.А. Постнов, И.Я. Хархурим. - Л.: Судостроение, 1974. - 342 с.

3. Курочка, К.С. Оценка эффективности некоторых методов решения краевых задач нелинейной теории упругости: тезисы докладов международной научно-технической конференции «Актуальные проблемы развития транспортных систем» /К.С. Курочка -Гомель: БелГУТ, 1998. - С. 204-205.

4. Золотухин, Ю.Д. Результаты натурных испытаний многопустотных плит перекрытия

экспериментального жилого дома с широким шагом несущих железобетонных поперечных стен в г. Речица /Ю.Д. Золотухин, В.С. Кульбицкий //Пространственные конструктивные системы зданий и сооружений, методы расчёта, конструирования и технология возведения: тр. междунар. науч.-техн. конф. В 2 т. Минск,

10-12 окт. 2001 г.- Минск, 2001.-Т. 1.-С. 59-70.

5. Тимошенко, С.П. Пластинки и оболочки /С.П. Тимошенко, С. Войновский-Кригер. -Москва: Физматгиз, 1963. - 636 с.

6. Arnold D.N., Falk R.S. Asymptotic analysis of the boundary layer for the Reissner-Mindlin plate model, SIAM J. Math. Anal., 27:486-514, 1999.

7. Reissner E. Reflections on the theory of elastic plates, Applied Mech., Reviews, 38:14531464, 1985.

8. Курочка, К.С. Расчёт деформаций неоднородных плит геометрически сложной структуры /К.С. Курочка. - Вестник ГГТУ им. П.О. Сухого. - 2003. - № 1. -С. 39-47.

- Получено 02.04 2004 г.