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

Известия Тульского государственного университета Естественные науки. 2010. Вып. 1. С. 123-136

^ Механика

УДК 539.3

Решение краевых задач нелинейной термоупругости *

Д. В. Христич, М. Ю. Соколова

Аннотация. Выполнена постановка термомеханических задач о конечном деформировании изотропных и анизотропных материалов под действием внешних механических и тепловых факторов с учетом их взаимного влияния. Проведен анализ влияния учета геометрической нелинейности на напряженно-деформированное состояние при конечных деформациях упругих полых цилиндров различной толщины.

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

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

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

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

* Работа выполнена при финансовой поддержке РФФИ (проект № 10-01-97501-р_центр_а), ФЦП «Научные и научно-педагогические кадры инновационной России» (контракт П1125) и АВЦП «Развитие научного потенциала высшей школы» (проект №2.1.1/941).

приращений, вызываемых приращениями внешних воздействий в данный момент времени.

Кинетические формы равновесия в настоящее время записываются в дифференциальном и вариационном видах. В работе [1] условия равновесного протекания процесса было предложено формулировать в виде вариационного соотношения, полученного на основе модифицированного принципа Журдена. Это вариационное соотношение записано в текущей конфигурации, что позволяет использовать при решении задач смешанный Эйлеро-Лагранжев подход. Сформулируем вариационное условие сохранения равновесия в отсчетной конфигурации.

Условие равновесия элементарного материального объема йУо сплошной среды заключается в равенстве нулю главного вектора действующих на него поверхностных и массовых сил:

— вектор внешних массовых сил.

Для элементарного объема в начальном базисе условие равновесия имеет

вид

Умножим полученное выражение (3) на отличную от нуля во всех точках рассматриваемого материального объема вариацию скорости 5Р и проинтегрируем по начальному объему:

О

Е = V Р + р,

О

где V = — оператор Гамильтона в исходной системе координат,

Р = (Ф-1)Т  £ = £  Ф-1

— тензор условных напряжений Пиола-Кирхгофа,

О

Ф = V х — аффинор деформаций,

£ = в8 — обобщенный тензор напряжений Коши, 8 — тензор истинных

напряжений Коши,

Є = ТЛТ" — изменение элементарного материального объема,

Дифференцируя уравнение (2) по времени, получим

Преобразуем первое слагаемое в левой части соотношения (4) под знаком интеграла, используя известное из тензорного анализа соотношение

О / О \ о

V

и теорему Остроградского-Гаусса:

/ (V Р  5МУо = / V  (1ь  5у) йУо - / Р  -5(у V)<1У0 =

Ло V / ,/Уо ,/Уо

= / по  Р  5МУо - / Р  -5(у V)йУо = / Ро  5уйХо - / Р  -5(у V)йУа, J'Zo ,/Уо Уе 0 Лс

где Ро = по  I3 — скорость изменения вектора внешней нагрузки, приложенной на внешней поверхности Хо с вектором единичной нормали по.

После подстановки преобразованного слагаемого в выражение (4) получим условие равновесного протекания процесса деформирования в вариационной форме:

/ Р  5

]Уо JE о лУо

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

Вариационное соотношение (5) эквивалентно локальному условию сохраняющегося равновесия в виде

V Р + ^ ёУо

= 0 ^ Ьо.

Если подчинить вектор К начальному условию К = 0, то локальное условие равновесия (2) будет удовлетворено вследствие выполнения условия (6).

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

В вариационное соотношение (5) входит производная тензора напряжений Пиола-Кирхгофа по времени I3. Получим связь между тензором Р и «повернутым» обобщенным тензором напряжений Коши

£к = Я  £  Я-1,

где Я — ортогональный тензор поворота, входящий в полярное разложение аффинора деформаций Ф = и  И,, и = Ит — левая мера искажения, а также между производными по времени указанных тензоров напряжений. Из определений (1) и (7) следует, что

Р = (И-1  £к  И)  Ф-1 = И-1  £к  И-1.

Дифференцируя последнее соотношение по времени, получим

Р = (И-1)  £к  И-1 + И-1  £ к  И-1 + И-1  £к  (И-1) .

Подставим выражение для I3 в вариационное соотношение (5):

/ [(И-1)  £к  И-1 + И-1  £к  И-1 + И-1  £к  (И-1)]   5(ь V)йУа =

= / Ро  Ш£о + / Р  5Шй.

,/Ео JVo

Соотношение (8) представляет собой условие равновесного протекания процесса деформирования, записанное через «повернутый» обобщенный тензор £к в отсчетной конфигурации.

Тепловое воздействие на тело определяется притоком тепловой энергии через поверхность £, ограничивающую объем V, и местными источниками тепла, которые имеют физико-химическую природу и в дальнейшем не учитываются. Общий тепловой поток через поверхность £ за время АЬ представим в следующем виде:

AQ — —

/ д

Уе Уео

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

Используя теорему Остроградского-Гаусса, из выражения (9) получим

AQ — — Ш V -дойУоАі.

V)

Тепловое воздействие на единицу начального объема получим в виде

О

в! Q — — V -доАі.

Вектор теплового потока полагаем связанным с неоднородным температурным полем законом Фурье [2]

О

до

где Ло — тензор теплопроводности, определяемый природой вещества,

О

V Т — градиент температуры.

Для изотропного материала полагают Ло — АЕ, где А — коэффициент

О

теплопроводности. При этом закон Фурье принимает вид до (X, і) — —А УТ. В анизотропном материале тензор Ло отличен от шарового и определяется присущей материалу симметрией термических свойств.

Соотношение (11) позволяет определить скорость притока тепла к еди-

О

нице объема материала в виде Q — —V -до, которая с учетом закона Фурье для теплового потока (12) принимает вид

Q — V

В изотропном материале 0 = АУ2Т, где V2 = V -V — оператор Лапласа. Скорость притока тепла 0 входит и в закон изменения энтропии

0 = ро

где 5, т — соответственно скорость изменения удельной энтропии и скорость диссипации.

Приравнивая правые части (13) и (14), получим уравнение теплопроводности в общем виде

ро

которое должно быть конкретизировано на основании представлений для энтропии 5 и скорости диссипации т в рассматриваемой модели материала.

Умножим левую и правую части уравнения (15) на отличную от нуля во всех точках рассматриваемого материального объема вариацию скорости изменения температуры 5Т и проинтегрируем это уравнение по начальному объему:

/ ро(ТБ — т)5ТвУо — / ІУо Jvo

V ( Ло

5ТвУо.

Преобразуем правую часть выражения (16), используя теорему Остроградского-Гаусса:

/ V •( Ло  V Т ) £ТвУо

.)Уо

V V  ^Ло -РтёТ^ вУо — I (Ло -^7^  б(VТ ) вУо —

= - Уе по  до5Тё£о - ^ ^Ло  V Т^  ^V Т^ dVо.

Подставив преобразованное слагаемое в уравнение (16), получим уравнение теплопроводности в вариационной форме:

/ ро(Т5 - w)5TdVо = -/ по  qо5Td£о -/ (л •VT^  ^VT^ dVо.

^0 ./Ео Яо \ УЧУ

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

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

Основными неизвестными в рассматриваемой термомеханической задаче являются вектор перемещений п, тензор напряжений 8 и температура Т. Используем определяющие уравнения, предложенные в работе [3], которые по форме совпадают с выражениями, введенными для линейной анизотропной термоупругости в книге [4]:

£к = N -М - В

Здесь £к — обобщенный тензор напряжений (7), N В — тензоры упругих

и термических свойств материала, М — неголономная мера деформаций,

которая определяется из уравнения

М = = ¥к = И  W  И-1,

dt

где W = 1 (Рь + — тензор деформации скорости.

В используемой модели выражение для энтропии имеет вид [3]

1 Т

Б =— В -М + с£ 1п

ро —о

а скорость производства диссипации для рассматриваемых обратимых процессов ад — 0. В выражении (20) для энтропии ро — начальная плотность,

е£ — удельная теплоемкость материала. Тогда

1 • 1 Т

Б = — В-М + — В-М + е£ -

Ро Ро -

и уравнение теплопроводности в вариационной форме (17) принимает вид / (В  -МТ + В  МТ + о£роТ) 5ТйУо =

= - ^ по  дЪ^Жо — ^ ^Ло  V ^V т) йУо-

Конкретизируем соотношения (18), (21) для изотропного материала. В изотропном материале тензор упругости имеет вид

N = ^К + 3О^ (ё\ё\ё\ё\ + 62626262 + 63636363) +

+ ^К — зО^ (61616262 + 62626161 +62626363 +63636262 + 63636161 +61616363) +

+20 (61626162 + 61626261 + 62616162 + 62616261 + 62636263 + 62636362+

+63626263 + 63626362 + 63616361 + 63616163 + 61636361 + 61636163),

где К, О — модули объемной упругости и сдвига.

Тензор температурных напряжений имеет вид

В = 3аК (6161 + 6262 + 6363),

где а — линейный коэффициент температурного расширения.

Тогда определяющие соотношения (18) могут быть записаны в виде

Хк = К [в — 3л/3а(Т — То)) 1о + 20М,

где М = М — (М  -1о)1о — девиатор тензора М, 1о = -^= (6161 + 6262 + 6363).

Полагая 13 = 0 ив соответствии с (19) М = Wк, в левой части уравнения теплопроводности (21) получим

В  -МТ + В  -МТ = 3^3аК 1о  •WкT = 3аКвТ,

О ООО

а в правой части Ло  УТ = АУТ, УТ  $(УТ) = У2Т$Т. Тогда для изотропного материала уравнение теплопроводности в вариационной форме (21) принимает вид

/

Ьо Уео -'у0

Эволюционные соотношения для перемещений, напряжений и температуры имеют вид

*

Начальные условия характеризуют состояние тела в начальный момент времени £о:

и

Граничные условия статического типа предусматривают задание в каждой точке поверхности Хр закона изменения внешних сил как функции времени

Ро = Ро

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

и = и*

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

Ро = -Ро*

Поверхности Хр, Хи, Хри не пересекаются: Хо = Хр и Хи и Хри.

Также необходимо задать на части поверхности Ху закон изменения температуры

Т = Т*

на части поверхности Хд — закон изменения теплового потока

д = д*

На части поверхности Хс может происходить свободный теплообмен с

окружающей средой, для описания которого принимается закон Ньютона [2]. Это приводит к граничному условию

О

по  Ло -УТ + а0

где а0 — коэффициент теплообмена, Тн = Тн (х, *) — известная температура

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

дТ

А— + а0

дп

Поверхности Ху, Хд, Хс не пересекаются: Хо = Ху и Хд и Хс.

При задании граничных условий полагаем, что функции (25)—(28) являются дифференцируемыми функциями времени. Тогда выражения (25)—(28) могут быть преобразованы к виду

Ро = Р0*

V = У*

Ро = -Ро*

Т = Т*

Таким образом, постановка связанной краевой задачи нелинейной анизотропной термоупругости включает определения мер деформаций и кинематические соотношения (19), уравнение равновесия сплошной среды в вариационной форме (8), уравнение теплопроводности в вариационной форме (17), соотношения, определяющие связь между напряжениями, конечными деформациями и температурой в нелинейно упругом теле (18), эволюционные соотношения (23), начальные условия (24) и граничные условия (29)—(35).

Интегрирование системы уравнений поставленной краевой задачи с учетом заданных граничных и начальных условий выполняется численными методами конечных элементов и пошагового нагружения [1, 5, 6].

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

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

КмьУь = Рм,

где Уь — обобщенные компоненты скорости, а именно: три компоненты вектора скорости V и скорость изменения температуры Т.

Численное интегрирование системы вариационных уравнений (8), (17) для определения компонент вектора узловых перемещений и величин узло-

вых температур проводится в соответствии с разностной схемой аппроксимации функций пь(Ь) и Ть(Ь) первого порядка по времени Ь. Принятая схема аппроксимации позволяет использовать для описания развития процесса во времени метод пошагового нагружения [1, 5, 6].

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

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

Были проведены расчеты напряженно-деформированного состояния в прямом круговом тонкостенном цилиндре (К/К = 0,1), нагруженном внутренним давлением. Схема нагружения и закрепления торцов цилиндра показана на рисунке 1(а). Процесс нагружения полагался изотермическим и проводился при температуре Т = 273К.

а) б)

Рис. 1. Схемы нагружения тонкостенного

цилиндров

При данной схеме нагружения проведены расчеты для изотропных материалов трех типов с различными жесткостью и сжимаемостью: для матери-

ала 1 задано Е = 105 МПа, V = 0,25 (сталь), для материала 2 — Е = 8 МПа,

V = 0,25 (сжимаемый мягкий материал), для материала 3 — Е = 8 МПа,

V = 0,4 (резина).

Характерный размер конечного элемента при разбиении расчетной области составлял 0,02 толщины стенки цилиндра К. Величина шага нагружения при пошаговом решении задачи выбиралась такой, чтобы приращение деформаций за один шаг не превышало 0,05%. Расчеты проводились до деформаций 20%.

В ходе расчетов предполагалось определить, в какой мере учет геометрической нелинейности влияет на результаты. С этой целью проанализированы зависимости перемещений внутренних точек цилиндра п от прикладываемого давления р, найденные с учетом и без учета геометрической нелинейности. Эти зависимости представлены на рисунке 2. Эти решения совпадают до уровня деформаций 0,6%. При достигнутых деформациях 20% отличия линейного и нелинейного решений составляют 40-50%.

О 0Д1 008 0.12 0,16 0.2

_______________________________________________р, МПа

линейное решение * нелинейное решение

Рис. 2. Зависимости и(р) для материала 3

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

Выполнен расчет напряженно-деформированного состояния в полом прямом круговом толстостенном цилиндре, для которого К/К = 0,25. Схема нагружения и закрепления торцов цилиндра показана на рисунке 1(б).

Расчеты выполнялись для отношений Ь/К = 8 и Ь/К = 12. В качестве материала цилиндра принималась резина (Е = 8 МПа, V = 0,4). Температура полагалась постоянной и равной 273 К.

Характерный размер конечного элемента при разбиении расчетной области составлял 0,1 толщины стенки цилиндра. Величина шага нагружения выбиралась так же, как при решении предыдущей задачи. Расчеты проводились до деформаций 20%.

Наибольшее перемещение в радиальном направлении совершает точка, лежащая на внутренней поверхности цилиндра на одинаковом расстоянии от

оз

0,23 02 0,15 0,1 0,05 0

О 0,02 006 0.1 0.Н 0.18

_____________________________________________р, МПа

— м атериал 1 * м атериал 2 м атериал 3

Рис. 3. Зависимости и(р) для различных материалов

его концов. Графики зависимостей радиального перемещения этой точки от приложенного внутреннего давления приведены на рисунке 4.

0,015 0,01 0,005 0

0 0^ 0,16 021 032 0,1

р, МПа

линейное решение 1  нелинейное решение 1

—*—линейное решение 2 —о— нелинейное решение 2

Рис. 4. Зависимости и(р) для короткого (1) и длинного (2) цилиндров

В данных расчетах интерес представляет изменение формы цилиндра с ростом давления. Геометрически нелинейная постановка задачи позволяет проследить за образованием и развитием «бочкообразности», характер которой зависит от относительной длины первоначального цилиндра. На рисунках 5, 6 изображены формы расчетной области (половины осевого сечения цилиндра) при различных уровнях деформаций, полученные в процессе расчета.

Рис. 5. Формы расчетной области для короткого цилиндра при различных деформациях: а) — 5%; б) — 10%; в) — 15%; г) — 20%

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

Рис. 6. Формы расчетной области для длинного цилиндра при различных деформациях: а) — 5%; б) — 10%; в) — 15%; г) — 20%

деформаций 1,5%. При деформации 20% наибольшая величина радиального перемещения точки на внутренней поверхности цилиндра, рассчитанная по нелинейной теории, превышает соответствующую величину, рассчитанную по линейной теории, на 28,9% для короткого цилиндра и на 26,7% для длинного. Максимальное радиальное перемещение точки короткого цилиндра больше, чем в длинном цилиндре, на 4,8%.

Список литературы

1. Адамов В.И., Маркин А.А. Моделирование процессов обработки давлением осесимметричных изделий // Изв. вузов. Машиностроение. 1989. № 12. С. 104-108.

2. Ильюшин А.А. Механика сплошной среды. М.: Изд-во МГУ, 1990. 310 с.

3. Маркин А.А., Соколова М.Ю. Вариант определяющих соотношений нелинейной термоупругости для анизотропных тел // Прикладная механика и техническая физика. 2003. Т. 44, № 1. С. 170-175.

4. Новацкий В. Теория упругости. М.: Мир, 1975. 872 с.

5. Маркин А.А., Соколова М.Ю. Вариант теории конечного упругопластического деформирования // Прогрессивная технология приборостроения. / Межвуз. сб. науч. тр. М.: Изд-во ВЗМИ, 1987. С. 57-61.

6. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритмы, приложения. М.: Наука, 1986. 231 с.

7. Астапов В.Ф., Маркин А.А., Соколова М.Ю. Определение упругих свойств материалов из опытов на сплошных цилиндрах // Механика твердого тела. 2002. № 1. С. 104-111.

8. Белл Дж.Ф. Экспериментальные основы механики деформируемых твердых тел. Часть II. Конечные деформации. М.: Наука, 1984. 432 с.

9. Васин Р.А., Ильюшин А.А., Моссаковский П.А. Исследование определяющих соотношений и критериев разрушения на сплошных и толстостенных трубчатых цилиндрических образцах // Механика твердого тела. 1994. № 2. С. 177-184.

Христич Дмитрий Викторович (dmitro@tula.net), к. ф.-м. н., доцент, кафедра математического моделирования, Тульский государственный университет.

Соколова Марина Юрьевна (sokolova@tula.net), д. ф.-м. н., профессор, кафедра математического моделирования, Тульский государственный университет.

Nonlinear elasticity boundary value problems solution D.V. Khristich, M.Yu. Sokolova

Abstract. A propounding of thermomechanical problems about finite deformations of isotropic and anisotropic materials under the effect of external mechanical and heat factors with consideration of their mutual influence is fulfilled. An analysis of geometrical nonlinearity consideration influence on stress-strain state at finite deformations of elastic hole cylinders of different thickness is carried out.

Keywords: thermomechanics, thermoelasticity, finite deformations, variational principles, coupled boundary value problems.

Khristich Dmitry (dmitro@tula.net), candidate of physical and mathematical sciences, associate professor, department of mathematical modelling, Tula State University.

Sokolova Marina (sokolova@tula.net), doctor of physical and mathematical sciences, professor, department of mathematical modelling, Tula State University.

Поступила 19.01.2010