МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФИЗИКО-ХИМИЧЕСКИХ ПРОЦЕССОВ
УДК 539.2:544.2:678.01:519.7:539.3:517.958
ЭВОЛЮЦИОННЫЕ ИЗОТРОПНЫЕ СХЕМЫ ЧИСЛЕННОГО РЕШЕНИЯ НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ КВАЗИСТАТИЧЕСКОГО ДЕФОРМИРОВАНИЯ. Часть 4. Проекционно-сеточная постановка при конечных деформациях
АЛЬЕС М.Ю.
Институт механики УрО РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. Исследуются эволюционные схемы численного решения систем нелинейных проекционно-сеточных уравнений неупругого поведения полимерных изделий в условиях конечных деформаций. Дается проекционно-сеточная постановка задачи.
КЛЮЧЕВЫЕ СЛОВА: эластомерные композиты, нелинейное механическое поведение, конечные деформации, краевые задачи, метод конечных элементов, эволюционные численные схемы.
ДИФФЕРЕНЦИАЛЬНАЯ ПОСТАНОВКА
Сформулируем квазистатическую задачу нелинейной вязкоупругости при конечных деформациях. Пусть тело, имеющее в момент t0, соответствующий началу рассматриваемого процесса нагружения, однородную конфигурацию р, к моменту времени t в результате деформирования переходит в новую конфигурацию р. Будем полагать, что в момент t0 известны конфигурация р и значения всех параметров, характеризующих предысторию деформирования. Приняв в качестве меры деформации некоторой окрестности физической точки тела меру Коши-Грина [1] (5* — символ Кронекера; gtk — метрический тензор; V — набла оператор)
О = F^kF¡ gtк, (1)
где
^ = 5* + V и, (2)
запишем соотношения для тензора конечных деформаций Грина
а, = 0,5(О - gг] К (3)
или
а, = 0,5^ +V и +V1Uk V и* ). (4)
В (1) - (4) компоненты векторов и тензоров, ковариантные производные определены в базисе отсчетной конфигурации р. При решении задачи в р -базисе удобно перейти к тензору деформации Альманзи [1]
а, = 0,50^, Д. V£к). (5)
Нижними индексами у векторов и тензоров обозначены ковариантные компоненты, верхними - контравариантные. Крышечкой сверху обозначено, что компоненты берутся в
р -базисе.
В рамках теории малых деформаций последними слагаемыми в правых частях в (4), (5), пренебрегается, тензоры Грина и Альманзи равны друг другу. При этом отпадает необходимость в рассмотрении более чем одной конфигурации р . Перемещения же в геометрически линейной теории рассматриваются не как разность между векторами места частицы тела в р и р конфигурациях, а лишь как векторное поле, определенное в евклидовом пространстве [2].
Перейдем к уравнениям равновесия в объеме К и на части границы , подверженной воздействию поверхностных сил. Наиболее просто и физически наглядно записывается баланс сил в базисе актуальной конфигурации р ( Р - плотность; п; - внешняя нормаль к поверхности тела; Г — вектор внешних поверхностных сил; gi — вектор массовых сил; — части поверхности, на которых задаётся плотность внешних поверхностных сил)
+ р^] = 0, (6) = Гi. (7)
Трудности в (6), (7) не фигурируют в явном виде и заключаются в том, что, как правило, сама эта конфигурация и ее базис наперед неизвестны, тогда как отсчетная конфигурация входит в состав исходных данных. Сравнение эйлерова и лагранжева подхода рассматриваются во многих работах (см., например, [3]). Отметим только, что в задачах с определяющими уравнениями, зависящими от истории деформирования, преимущества имеет постановка в терминах отсчетной конфигурации. В этом случае интегрирование определяющих соотношений производится наиболее простым способом.
Введение симметричного тензора Пиола-Кирхгоффа р, определенного в
отсчетной конфигурации и связанного с тензором ст., простым соотношением [ 1 ] (13(А) = det [ А] — третий инвариант симметричного тензора второго ранга)
а -Р, (8)
>/?з(ё)
вместо (6), (7) дает [4]:
V; (Ра(5]а + Уаи>)) + pg] = 0, (9)
Ра(51 +Уаи)п\ = Г . (10)
Выражение для компонент вектора внешних поверхностных сил Г3, обусловленных действием поверхностного давления р , выглядит следующим образом [4]
Г3 = р413(0)0'ап1 (31 +Уаи). (11)
В (9) — (11) компоненты векторов и тензоров, ковариантные производные определены в базисе отсчетной конфигурации. При рассмотрении малых деформаций различие между лагранжевым и эйлеровым описанием исчезает, напряженное состояние определяется
единственным тензором О/ и уравнения равновесия в объеме тела V и на части границы принимают вид О +рg} = 0, ауп} = Г .
На части поверхности тела 5и могут быть известны перемещения. Кинематические граничные условия имеют вид
и.\ = и .
1 15 1
(12)
Замыкается постановка краевой задачи заданием связи между тензором напряжений Пиола-Кирхгоффа р и тензором меры деформаций Коши-Грина О,. Для
моделей деформирования (см. [5]) Свансона (— девиатор тензора р,
„ „ Г г Л1/р
\гЛр =1 | IР (Е)^Е — норма Лебега, Iи — интенсивность деформаций)
и Фарриса
запишем
t дЭ (Е)
= I (I и ,|| I и|| ® )| ^-£)—
I(1 и , II1Л ®) 1Л ®)
= е
Б—
1и
2/
Г I Лт
дЕ
г Г
1 - С
V V Г
лл
1-
Э +
VII и ||р J
1 -
V
I
I щ-Е)йЭц (Е)
р, = 2/ э + к -
Г~ 3//12 -1Г Л
у О I О 7" / 6 У
V
3 + 2/1 (а)
J
(13)
(14)
(15)
где — = (О) -1, ^(а) — первый инвариант симметричного тензора второго ранга;
~ 3К в +1 в-у
Кэ =---
э 3 + 2/1 (а) в 1 + в-г
(16)
Г = Э ^) + гтп •
1 Э Этп V ) ^ 1 Э ;
(17)
/Э, Г Э - определяются следующими соотношениями
для модели Свансона (13):
М (II1 Л®)
Г Г £ лл
1 _ 1 и
ГЭ=10 (11Л®)
1 - С
V
Г Г
1 - С
V V
I
V II и |I® J J
Л t
1 - 1 и
{г^ -Е)Э (ЕМЕ;
для модели Фарриса (14):
/э = /
Б—
4
1-
• +
т
Vм и||р J
(18) (19)
(20)
— (
гЭ = I
1 _ 1и
|гсг -ЕРУ- (ЕМЕ.
(21)
*
0
0
I
и
I
I
и
0
и I ®
0
и
®
I
и
и I ®
и 1I® I 0
При записи соотношений (15) в уравнениях (13), (14) было использовано интегрирование по частям для свертки Больцмана
t t
/Я(/(£) = 2рЭц (I)-Щ. (£)*£, т = -т, 2р = Я(0) (22)
0 0
Предполагается, что функция Э, (£) £ е [0, I] имеет интегрируемую производную Э, (£). Точка обозначает производную по аргументу.
Таким образом, решение задачи о напряженно-деформированном состоянии тела, имеющего в отсчетный момент времени t0 конфигурацию р, в постановке с учетом конечности деформаций сводится к определению в заданный момент времени t
конфигурации р(, перемещений и , деформаций , напряжений Р, удовлетворяющих уравнениям (4), (9), (15), при соответствующих граничных условиях (10), (12) и известных параметрах предыстории деформирования. Модели для приращения объема за счет порообразования приведены в [5]:
I! (а) = 3К~~~~~, /, (Р = 3/Г1 (О)(д/7з(О)1 (а) - ),
1 + в-у
где - девиатор тензора Р,; через у обозначено приращение объема, связанное с порообразованием в = а(К - О)- + у.
Аппроксимации для функции у можно принять исходя из моделей
Фарриса у(£ и ,0) = D(l и )"е°
или В.В.Мошева у(ат ,0) = С (е„ )аео,
отражающих закономерности объемных изменений в общем виде и не связанных по своей структуре (в части вкладов у) с понятиями "малые" или "конечные" деформации.
Отметим, что тензор Пиола-Кирхгоффа Р, - удобная вспомогательная величина,
введенная в данном случае для упрощения построения сеточных систем уравнений и вычислительных алгоритмов, и непосредственно не определяющая напряженного состояния, определение которого (при известном Р, ) требует возвращения к тензору Коши о, (8).
ПРОЕКЦИОННО-СЕТОЧНАЯ АППРОКСИМАЦИЯ
Для области V , ограниченной поверхностью 8 = 8и и запишем слабую форму Галеркина для уравнения равновесия (9) с естественным граничными условиями (10) на 8>а
| Р'а 5 + уу -\pgvJV -| = 0, (23)
V V
где Р, ^ - определяются соотношениями (15), (11).
Представив элементные распределения перемещений в виде
к, = УЛ 1 ег„,
получим следующую систему нелинейных алгебраических уравнений относительно неизвестных узловых значений перемещений к,т ( т = 1,М )
!(к)к + ^ (к) - G - рк) = 0, (24)
N
Цк)к = ^от//Зэ (к^'у + к;V>у - 2к^у, (25)
N с
F (к) = 2 от / ^, (26)
п=1 V'
^ =ФУ- + Р^ЧУ
Фу = З (к;VУlKfcsVу - 3Vawsg,1) + + Кэв^ -ГЭ - 2 (3 + 2/1 (а))-1 (3//э1 к - Гэ ) gy,
= & (к^у + к; vaу - 2 к\vуt ga + + ккVУiкksV>s -ткыVЧккVу^'" ) + + -Га -2(3 + 2/1(а))-1 [З]лэ1 к -Гэ)g,^
N
с = ^От , (27)
N
р(к)=^л ^от / 1У^п, (28)
п=1 е
/ = р (1 + 0)( gа + к^'у + к; Vау + кй VУtкks Vау )(# + к^„у) (т = 1,М; I,к,Г,г-Гп; в<-в)
Здесь у - базисная функция; N, М - число конечных элементов и узлов в триангуляции р; Гп - множество узлов элемента; £ - множество элементов, у которых хотя бы одна грань используется при аппроксимации £ет; Вп - множество узлов элемента, находящихся на £ет; От , Л £ - логические процедуры: От = 1, если узел У^Гп конечного элемента п соответствует узлу т = 1,М триангуляции р; О7т = 0 в противоположном случае; Л £ = 1, V« е £, Л £ = 0 в противоположном случае. Символ е обозначает, что при суммировании по впереди стоящему немому индексу производится последовательный перебор элементов данного множества; индексом п отличается принадлежность п -у конечному элементу.
Заметим, что структура уравнений (23) такова, что позволяет построить проекционно-сеточную систему по типу "шахматная доска", когда клетки одного цвета образуют группу уравнений относительно одной компоненты перемещений в узлах сеточной области, другого цвета - относительно другой компоненты. Использование такой структуры уравнений представляется интересным с позиций распараллеливания вычислений.
Таким образом, решение исходной дифференциальной задачи сводится к решению системы нелинейных алгебраических уравнений большого порядка относительно неизвестных узловых значений перемещений uim (m = 1,M ). Очевидно, что по сравнению с
подобными системами для малых деформаций, к трудностям, обусловленным нелинейным механическим поведением, добавились существенные трудности, связанные с изменением метрики континуума. Особенности решения таких систем и устойчивые всегда сходящиеся алгоритмы будут приведены последующих публикациях.
СПИСОК ЛИТЕРАТУРЫ
1. Лурье А.И. Нелинейная теория упругости. М. : Наука, 1980. 515 с.
2. Лурье А.И. Теория упругости. М. : Наука, 1970. 939 с.
3. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритмы, приложения. М. : Наука, 1986. 231 с.
4. Оден Дж. Конечные элементы в нелинейной механике сплошных сред. М. : Мир, 1976. 464 с.
5. Альес М.Ю. Феноменологические описания нелинейного сопротивления деформированию высоконаполненных эластомерных (нано) композитов. Часть 2. Конечные деформации // Химическая физика и мезоскопия. 2010. Т. 12, № 2. С. 219-223.
EVOLUTIONARY ISOTROPIC NUMERICAL SCHEMES FOR ACTION BOUNDARY VALUE PROBLEMS OF NONLINEAR QUASI-STATIC DEFORMATION. PART 4. PROJECTION-GRID PRODUCTION AT FINITE DEFORMATIONS
Alies M.Yu.
Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. We investigate the evolutionary scheme of numerical solutions of systems of nonlinear equations projection-grid inelastic behavior of polymeric products under finite deformations. Given projection-difference formulation of the problem.
KEYWORDS: elastomer composites, nonlinear mechanical behavior, finite deformation, boundary value problems, finite element method, numerical evolutionary scheme.
Альес Михаил Юрьевич, доктор физико-математических наук, профессор, главный научный сотрудник ИМ УрО РАН, тел. 89128563824, e-mail: [email protected]