CC BY
148
ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
УДК 519.63
П. Н. Вабищевич, М. В. Васильева
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЗАДАЧ ТЕРМОУПРУГОСТИ
В работе рассматривается связанная квазистационарная линейная система уравнений для температуры и перемещений, описывающая термоупругое состояние тела. Приводятся априорные оценки для дифференциальной задачи. После конечноэлементной аппроксимации по пространству и аппроксимации по времени с использованием двухслойных разностных схем с весами получаем дискретную задачу. Показана устойчивость получаемой дискретной задачи. Представлены результаты численного расчета модельной задачи термоупругости.
Ключевые слова: термоупругость, двухслойные схемы с весами, устойчивость, априорная оценка, метод конечных элементов, численное моделирование.
P N. Vabishchev, M. V Vasilieva
Numerical simulation of thermoelasticity problem
The paper considers the coupled quasi-stationary linear system of equations for temperature and displacement, which describes the thermoelastic behavior of the body. The stability estimates of the studied differential problem are presented. For the numerical solution, we approximate our system using the finite elements method and for approximation by time we use the standard two-layer weighted schemes. The discrete problem has proven stable. The results of the numerical simulation of the model problem of thermoelasticity are represented.
Key words: thermoelasticity, two-layer weighted schemes, stability, prior estimat, finite elements method, numerical simulation.
Многие прикладные проблемы математического моделирования связаны с необходимостью расчета напряженно-деформированного состояния твердых тел. Во многих случаях деформации обусловлены тепловым расширением. Для их исследования используются модели термоупругости [1-4].
Базовые математические модели термоупругости включают уравнения Ламе для перемещений и
ВАБИЩЕВИЧ Петр Николаевич - д. ф.-м. н., профессор, заведующий лабораторией Института проблем безопасного развития атомной энергетики РАН, ведущий ученый -научный руководитель ЦВТ ИМИ СВФУ им. М.К. Аммосова.
E-mail: [email protected]
ВАСИЛЬЕВА Мария Васильевна - к. ф.-м. н., начальник отдела параллельных технологий Центра вычислительных технологий Института математики и информатики СВФУ им. М.К. Аммосова,
E-mail: [email protected]
уравнение теплопроводности. Принципиальный момент заключается в том, что эта система уравнений завязана, уравнение для перемещений включает объемную силу, пропорциональную градиенту температуры, а уравнение для температуры включает член, описывающий сжимаемость среды.
В работе описан вычислительный алгоритм решения задач термоупругости. Он базируется на конечно-элементной аппроксимации температурного поля и поля перемещений по пространству [5-8]. Квазистационарная задача (стационарная для перемещений и нестационарная для температуры) решается на основе использования двухслойной схемы с весами [9]. Показана безусловная устойчивость приближенного решения по начальным данным и правой части при стандартных ограничениях на вес. Получена априорная оценка приближенного решения, которая согласуется с соответствующей оценкой для решения дифференциальной задачи. Возможности
вычислительного алгоритма проиллюстрированы данными расчетов по модельной двумерной задаче для тела прямоугольного сечения с круговым вырезом.
Постановка задачи. При механических и тепловых воздействиях в упругом теле возникают поля перемещений u, деформаций е и напряжений а. Пусть Т - постоянная абсолютная температура, при которой тело находится в начальном равновесном состоянии, а 0 - приращение температуры. Под механическим воздействием понимают воздействие на тело внешних сил (объемных или поверхностных), а под тепловыми
- процессы теплообмена между поверхностью тела и окружающей средой и выделения или поглощения тепла источниками внутри тела.
Математическая модель термоупругого состояния тела описывается связанной системой уравнений для перемещений u и приращения температуры 0 в области О:
и я = і2(Л) - пространство для перемещений.
Отметим )сновные свойства операторов задачи во введенных пространствах [10-14]. На множестве функций, удовлетворяющих приведенным выше условиям на границе, оіредетим оператор:
4 = — /Дг — ТА + /-О §гас1 сііу V.
Оператор А является самосопряженным и положительно оіфеделенньїм Н
(Ате, и) = 4# Л и), (4и,и) > 8и(и,ге), > 0.
Аналогично в Н опредев м оператор В:
Вв = — сііу 0+ |Е;гас1б).
нАи + (Л + н)graddivu — agradе = 0, (1)
дв дdivu с— + аТ ——---------div(Щgrad0) = /. (2)
Здесь X, ц - постоянные Ламе, к - коэффициент теплопроводности, с - удельная объемная теплоемкость при отсутствии деформгщий,
а — ат( ЗА р 2 ц), гае ит - температурный коэффициент линейного расширения, е -тензор деформаций:
1
£ = !
ист - тензор нержген ий:
а = АЧи1 + ! ц £.
Уравнение (1), (2) дополн=м соответствующими начальным
в (Х, 0) = 0
и граничными дП = Г* + Г+* = Г) + Г+ условиями:
Также определим операто]ры дивергенции Б и градиента О, причем
(СД, иЖ = — (Д,МиЖ.
Рассматриваемую зада.- дх= системы уравнений (1), (2) с использованием введенных операторов
можно здисать в виде задачи Коши для системы дифференциальн омператорных уравнений:
Au + аМП = О)
(3)
— (Св + aDu) + Вв = /(£). (4)
Здесь Св=ев и проведено обезразмесив=ие температуры на T. Для (3), (4) задается начальное условие для температуры:
е(О) в о.
(5)
Пусть Я - самосопряженный положительный
оператор, действующий в Н. Тогда можно ввести
энергетическое пространство НЯ, состоящее из
элементор Н со скалярным произведением
-n =0, x Є Г+*, u = 0, x Є Г)*,
ОН
-к—= 0, ХЄГ+, Н = 0, хЄГ0е.
on
Пуснь Н = Ь2 (І2) -гильбертово пространство для
давления со скалярным произведением и нормой в виде: (_u,v) = ju(rMX)dX, \\и\\ = (и,иї*
(u,v)R = (Ru, v), и нормой ||uHR = (Ru,u)2.
Теорема 1. Для решения задачи (3)-(5) верна априорная о ценка:
\iu(t(\i2A + \iem2c< н«сож +
t
+ lie(O)ll?+iJll/(s)||
о
g-i ds
(6)
б
Доказательство: Умножим скалярно в H уравне-
ние (3) на du/dt, а уравнение (4) на в. Это дает
І йи\ і йи\
(#и'л)+Чсе'л) = 0' (се’§)+а(°%’е)+те)=^е)-
Складываа эти равенства, получим
і аи\ і ав\
і'4 ив$)+{С в^)+()в,в) = (/'в )■
Во спользуемся нераве нстюм
и,ва<\т2в
и получим
du$
dt $
dB$ І
Ou> —) М I "В>—;— ) < —
dt $ 4
в-1.
Используя
edu $ 1 н
" lt'U) = 2 Tt (U’U)’
можно записать
^Аи,и)+ •
Откуда непосредственно и вытекает доказываемая оценка (4).
Апнооксимачия но пространству и времени.
Для численного решения системы уравнений (1),(2) построим конечно-элементную апщюксимацш) по пространству
Ийя всех фушщий vиq5 расчетной области О
шагом т > О : tn = пт, П = 0,1, ...,M
и обозначим
un и u(x, tn), Ип = И(х, tn) .
Воепользуемся ставдартныгми двухслойными схемами с весами и определим сле&ющии билинейные и линейные формы
c(M,%) = I CMcdx,
(о.
d(u,q= = I adiv uqdx,
& а
b(d,q)= I (.grad в, grad q)dx,
Jq.
Kf,q) = (f,q) = & f%dx,
Jq.
a(u,v) = I a(u)£(u)dx,
Jq
g (в,v) = I a(grad в,v)dx = О.
Jq.
Используем вариацион^ю формулировку задачи: шейти OEVg и и Є Vu такие, что выполняются
a(un+1,v) + g(Bn+1,v) = 0, <v eVU, (9)
n+1 „) _
-
-)u)s)v)Hx +
- a(grad в,
v)Hx = О, ("7)
Г da\vu f в)
I % ot Св* — | "Св-ссв^с-Ь
(en+1-en \ (un+1-un \
с----------,q\ + dl----------,q? + b(0n+1,q)
(8)
где
вt
5 5
+ J (8grad U, grad q)dx = J fqdx,
5 5
в E Vg = {вЄ H1: в(х) = О,x Є Гц } и X Vu = [H1(^)]d , d = 2,
Для апщюксимации по време ни определим равномерную для простоты сетку по времени с
= Kfa,q), є Vq.
(i„
Здесь
в
п-вl _ —ап+1
-в + (і — р-в
П
/П+1 _ адп+1 -|- (1 — д)0
и весовой параметр 0 < 7 < 1.
Отметим, что для всех u и V имеет место соотношение
d(u, у) = —д(у, и) .
Теорема 2. При а > 0.5 для решения задачи (9) (10) верна априорная оценка
|и^ + ||0п+1||2<||ип||2 +
-Н^+^Н/Г1!^ ,
где У - | |*;ь является нормой в сопряженном к пространстве.
Доказательство: Подставим
ид+1 = +ип+1 + (1 — а)пп
(11)
в у равнение для перемещений (9)
а ( и
п+1
+ д ( Р,пш1
ип+1 — гр1'
'а ’
" /V "
1П + 1
и положим % = 9™ вурашнении длядавления (10):
1#п+1 — Нп Л
&------------------,0г1 ) + й
/ип+1 — ип
.ег1) +
+ "(ЯГ^Г1) = (/"+1,0£+1 )■
Сложим полученные уржнения
а ( и*+1,
ип+1 - иг
/0пв1-0п '
+ с(---------О п
" /V " ж испг^нгч = е/^+^Г1)-
^п+1 — ^71,
— 1- + (Г — (гп+1 — ГП)
и тем что для симметричной билинейной формы
г{и,у) = г{у,и) имеет4 место
е(и + У,и — У Т = е(и, и Т — е(с, сТ .
С учето м этого пол^^шм соедующче неравенство:
-(||И"+1 И* - \\и"\\2а) + -(\Г+1\|* - ||Щ*) + 1л
+ О7-г)** Н“п+1 _ “п11а +
+ ||вп+1-вп||2)<7||/;
\*,Ь>
в
+
из которого при а > 0.5 и следует доказываемая оценка устойчивости по начальным данным и правой части. Численные результаты. Проведем численное моделирование задачи термоупругости в расчетной области, изображенной на рис. 1. Параметры задачи приведены в табл. (система единиц СИ). Расчетная сетка содержала 5825 треугольных ячеек и была построена с использованием программы Gmsh [15] (рис. 2). Моделирование проводилось на расчетное время, равное одному часу. Распределение температуры на различные моменты времени и перемещений в конечный момент приводятся на рис. 3, 4.
Моделирование проводилось с использованием аппроксимационных полиномов первой степени для температуры и второй степени для перемещений. Для решения возникающей системы линейных уравнений использовался стандартный прямой метод Ши-фак-
Поскольку имеет место нераве+ство
*,ь>
тогда при е=1 по лучи м
а I и
п+1
+
(вп+1-вп \ 1 ^ _
+ & (—п—, вп+- # * * И/.П+1 н,ь-
Воспользуемсе соотношением
П + 1
УП — ОУ
п+1
+ (1 — о)уп —
Рис. 1. Расчетная область
Таблица
Параметры задачи
Название Обозначение Значение
Модуль упругости Е 106
Коэффициент Пуассона ц 0.3
Объемная теплоемкость С 4 106
Теплопроводность к 30
Коэффициент линейного расширения CCj' 10-4
Начальная температура Т. 10
Температура на внутренней границе г, 200
Рис. 2. Расчетная сетка
торизации. Для численного решения использовалась открытая библиотека FEniCS [16,17], а для визуализации полученных результатов использовалась программа визуализации научных данных Paraview [18].
Проведенные расчеты демонстрируют высокую эффективность вычислительного алгоритма, его робастность. Расчеты выполнены на вычислительном кластере «Ариан Кузьмин» СВФУ им. М.К.Аммосова.
Работа выполнена при финансовой поддержке ЗАО «Оптоган» (договор №02.G25.31.0090) и гранта РФФИ (проект №13-01-00719А).
Л и т е р а т у р а
1. Biot M. A. Thermoelasticity and irreversible thermodynamics // Journal of Applied Physics. 1956. Vol. 27, no. 3. Pp. 240-253.
2. Lubliner J. Plasticity theory. Dover Publications, ew York, 2008.
3. Simo J., Hughes T. Computational inelasticity. Springer, New York, 1998.
4. Nowacki W. Dynamic problems of thermoelasticity. Pergamon Press, Warsaw, 1986.
5. Hughes T. J. The finite element method: linear static and
200
160 120 80 j 40 10
Рис. З. Распределение температуры
0.23
-0.2
0.1
0
X Y 1 hour
Рис. 4. Распределение перемещений
dynamic finite element analysis. Dover Publications, New York, 2012.
6. Zienkiewicz O. C., Taylor R. L., Zhu J. Z. The Finite Element Method: Its Basis and Fundamentals. Butterworth-Heinemann, Burlington, 2005.
7. Zienkiewicz O. C., Taylor R. L. The Finite Element Method: Solid Mechan- ics. Butterworth-heinemann, Burlington, 2000.
8. Brezzi F., Fortin M. Mixed and hybrid finite element methods. Springer, NewYork, 1991.
9. Самарский А. А. Теория разностных схем. М., Наука, 1989.
10. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.,Наука, 1978.
11. Gaspar F., Lisbona F., Vabishchevich P. A finite difference analysis of Biot’s consolidation model // Applied numerical mathematics. 2003. Vol. 44, no. 4. Pp. 487-506.
12. Lisbona F. J., Vabishchevich P. N. Operator-splitting Schemes for Solv- ing Unsteady Elasticity Problems // Comput. Methods Appl. Math.2001. Vol. 1, no. 2. Pp. 188-198.
13. Afanas’eva N., Vabishchevich P., Vasil’eva M. Unconditionally stable schemes for convection-diffusion problems // Russian Mathematics. 2013. Vol. 57, no. 3. Pp. 1-11.
14. Vabishchevich P. N., Vasil’eva M. V. Explicit-implicit schemes for convection-diffusion-reaction problems // Numerical Analysis and Appli- cations. 2012. Vol. 5, no. 4. Pp. 297-306.
15. Software package Gmsh. http://geuz.org/gmsh/
16. Anders Logg G. N. W., Kent-Andre Mardal. Automated Solution of Differ- ential Equations by the Finite Element Method.2011.
17. Software package FEniCS. http://fenicsproject.org/
18. Software package Paraview. http://www.paraview.org/
