Эволюционные изотропные схемы численного решения нелинейных краевых задач квазистатического деформирования часть 1. Особенности сходимости нелинейного численного решения в условиях малых деформаций Текст научной статьи по специальности «Математика»

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФИЗИКО-ХИМИЧЕСКИХ ПРОЦЕССОВ

УДК 539.2:544.2:678.01:519.7:539.3:517.958

ЭВОЛЮЦИОННЫЕ ИЗОТРОПНЫЕ СХЕМЫ ЧИСЛЕННОГО РЕШЕНИЯ НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ КВАЗИСТАТИЧЕСКОГО ДЕФОРМИРОВАНИЯ Часть 1. Особенности сходимости нелинейного численного решения в условиях малых деформаций

АЛЬЕС М.Ю.

Институт механики УрО РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34

АННОТАЦИЯ. Рассматриваются особенности решения систем проекционно-сеточных уравнений неупругого поведения полимерных изделий в условиях малых деформаций, особенности расчета среднего напряжения и относительного изменения объема и решения нелинейных систем конечно-элементных уравнений итерационными методами.

КЛЮЧЕВЫЕ СЛОВА: эластомерные композиты, нелинейное механическое поведение, краевые задачи, метод конечных элементов, особенности сходимости итерационных методов.

ВВЕДЕНИЕ

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

Наиболее распространенными методами решения проекционно-сеточных систем нелинейных уравнений большого порядка являются стационарные итерационные методы [2 - 4]. Практические проблемы сходимости обсудим на примере методов простой итерации и Ньютона. Для простоты изложения и без ущерба существу проблемы ограничимся рассмотрением модели деформирования без учета температурной и полимеризационной (химической усадки) составляющих

S- = -rj 0)

о = KQQ, (2)

где S., Э. - девиаторы тензоров напряжений о" и деформаций s.; o,Q- среднее напряжение и относительное изменение объема; /иэ, Кэ - эффективные (кажущиеся) модули сдвига и объемного растяжения-сжатия; Г. - тензорная функция. Различные представления

нелинейных функций Мэ,Кэ,Гэ. для высоконаполненных эластомеров рассмотрены в [1]. В простейшем случае мэ = М = const Кэ = К = const, Г.. = 0 имеем закон Гука.

t

При ju3 = д = const Кэ = K = const Гэ = Jr(t -£) Э (£) d£, где Г - регулярная часть ядра

0

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

ПРОЕКЦИОННО-СЕТОЧНАЯ ПОСТАНОВКА

Для области V, описываемой ограниченной поверхностью S = Su U Sa, Su П Sa = 0 запишем слабую форму Галеркина для уравнения равновесия с естественными граничными условиями на Sa

Jo4 V mdV - Jpg1 VmdV - J f1 V JS = 0 (3)

V V s,

и уравнения, связывающего относительное изменение объема и среднее напряжение

J(vy-aK-1 )vmdV = 0. (4)

V

Здесь иг - вектор перемещений; р - плотность; g1 - единичный вектор массовых сил; f1 - вектор внешних поверхностных сил; vm - пробные функции; Su, Sa - части поверхности, на которых заданы перемещения и напряжения. Произведя триангуляцию p и представив элементные распределения перемещений и среднего напряжения в виде (vi - базисная функция) иг = VA, a = Vi,, i ёГя из (3), (4) с учетом (1), (2) получим следующую систему алгебраических уравнений относительно неизвестных узловых значений перемещений uim , и среднего напряжения am ( m = 1, M )

m

N

V

m

n=1 v

2 П J Дэ (u/v' Vi + uv1 Vi) + [ 1 - 3 эКэ-1 I Vi,g

2 —1 L, „

V' VydV„ -

N NN

П J(20^3gy + Г/)V'VydV„ - 2 П Jpg1 VydV„ -^ЛL^ J f1 VedS„ = 0, (5)

V n=1 V n=1 S

и f n KJn.

2 от Д и ; - К;> 1а 1) dVn = 0, (6)

п=1 V

у п

где т = 1,М; I,уеГп; РеВп; NМ - число конечных элементов и узлов в триангуляции р; Гп - множество узлов элемента; £ - множество элементов, у которых хотя бы одна грань используется при аппроксимации $ ; Вп - множество узлов элемента, находящихся на $ ;

p > n ^ 1 ^ —' p :

Ч л*-»* ** .wiMiv HJJnm

nm, ЛL - логические процедуры: nm = 1, если узел У^Гп конечного элемента n

соответствует узлу т = 1,М триангуляции р, От = 0 в противоположном случае; ЛП = 1, Уп е П, ЛП = 0 в противоположном случае. Символ е обозначает, что при суммировании по впереди стоящему немому индексу производится последовательный перебор элементов данного множества; индексом п отличается принадлежность п -му конечному элементу. Граничные условия на являются естественными для уравнения равновесия. Граничные условия в перемещениях удовлетворяются непосредственной заменой соответствующего уравнения в (6), например, для узла т, тождеством игт = и

*

где игт - заданное значение перемещения. В уравнениях (5), (6) на параметр Кэ не

накладывается никаких ограничений. Матрица системы, при этом, имеет размерность 4М (в двумерных задачах 3М).

РАСЧЕТ СРЕДНЕГО НАПРЯЖЕНИЯ И ОТНОСИТЕЛЬНОГО ИЗМЕНЕНИЯ ОБЪЕМА

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

котором имеет только стационарное решение

Kf 1д% а + е- Кэ-1а = 0, f = 0, о = 0.

(7)

Здесь £- фиктивное безразмерное время; Kf - параметр алгоритма, размерность которого совпадает с размерностью объемного модуля. Аппроксимируя (7) по неявной схеме с шагом Д£ f > 0 будем иметь

N s+l s+1 s+

J "I - (K- - Kf - Kf:Vl °t

n=l V

¥ydVn = 0, Kf = Kf A£ f.

ff

ff

(8)

Уравнения (8) решаются итерационно совместно с уравнениями (5), в которых под узловыми значениями перемещений и\ и среднего напряжения ае следует понимать их (S +1) -е приближения. Отметим, что система (5), (8), за счет зависимости эффективных модулей |дэ, Кэ от искомого решения, является нелинейной. Особенности решения систем нелинейных сеточных уравнений неупругого деформирования обсуждаются ниже. Здесь для определенности и простоты изложения полагаем, что Цэ, Кэ в рассматриваемый момент времени t + At вычислены по решению, полученному в предыдущий момент времени (t + At) = дэ (t), Кэ (t + At) = Кэ (t) и при итерациях не изменяются. Из анализа последовательности векторов ошибок следует, что достаточным условием сходимости процесса (5), (8) к решению исходной системы (5), (6) является неравенство (H + D) 1D < 1, где H - матрица системы (5), (6); D - диагональная матрица D = diag {... ,0,..., (К-1 )1,..., (К-1 )м}. Итерационные параметры (Kff ) m = 1, M , как это

следует из условия сходимости, должны принадлежать нижней грани области Kf > 2Кэ при Кэ > 0; верхней грани области Kf < 2Кэ при Кэ < 0. Для линейной упругой модели имеем Кэ = К = const > 0 и соответственно Кff > 2К . В нелинейных моделях (см. например, [1]), в которых учитывается относительное изменение объема при сдвиге, величина и знак эффективного модуля Кэ определяются совместным влиянием среднего напряжения и

сдвиговыми деформациями. Например, в случае Кэ = К(1 + o~lKD(lи)neqa) [1] при о> 0 или а< KD(lи )neqa в изменении объема превалирует дилатационный механизм и Кэ > 0, Kff > 2Кэ. При -KD(l и )neqa <а< 0 превалирует дилатансия и Кэ < 0, Kf < 2Кэ. Отметим, что при |о| ^ KD(lи )neqa имеем |Кэ| ^ да, е ^ 0, и решение задачи в перемещениях в этой части области интегрирования невозможно (или некорректно при > 102 ).

Для экономии ресурсов ЭВМ удобно разделить единую сеточную систему на две взаимосвязанные системы уравнений меньшего порядка соответственно относительно перемещений и среднего напряжения. Подставив (S + 1)-е приближение среднего

напряжения в уравнения равновесия и воспользовавшись слабой формулировкой Галеркина будем иметь следующую систему для отыскания (£ + 1) -го приближения перемещений

N

V

т

п=1 V

N

( ,5+1

Д э

и ^

5+1

и\ V1vl +11 - 3 Д-^1 I К#< ^ ]Уг -

V

2

, +1 а '

-£ о: | г-- V,уу¿V - £ Щ | X11 - 3 Д-К-1 IV IV'уу¿V

п=1 V

N

\У

'т

п=1 V

у п

-£ от | р*¿V, - £ л^т I ^^=о,

(9)

ап \-1

где К, = К#К- (К, - Кз Г , К, * Кэ, Х = К (К, - К )" .

Среднее напряжение на (£ + 1)-й итерации отыскивается из системы (8). Размерность системы уравнений можно уменьшить, положив а поэлементно постоянным и сведя, тем самым, систему (8) к группе скалярных уравнений. Расчеты показали, однако, что кусочно-линейные аппроксимации а дают более устойчивое его распределение при решении.

ОСОБЕННОСТИ РЕШЕНИЯ НЕЛИНЕЙНЫХ СИСТЕМ КОНЕЧНО-ЭЛЕМЕНТНЫХ УРАВНЕНИЙ

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

зависящих от реализуемого напряженно-деформированного состояния, предыстории деформирования, характеризуются сильной нелинейностью коэффициентов матрицы и вектора правой части. По величине соседние коэффициенты в зависимости от принадлежности узлов зонам растяжения или сжатия, ветви активного нагружения или разгрузки, отличаются на порядки. Наиболее распространенными на практике методами решения систем нелинейных уравнений большого порядка являются стационарные итерационные методы, такие как, Ньютона, простой итерации, Пикара и т.д. [2 - 4]. Возникающие практические проблемы сходимости обсудим на примере методов простой итерации и Ньютона применительно к модели (1), (2) в которой (Фаррис, [1])

(

Дэ =Д

к.

!

е и

Кэ = К(1 + а-1 КО(1 и) "е^)-1, Г

9 (

= е

ир

I,

и.

-1

Для простоты изложения и без ущерба существу проблемы ограничимся рассмотрением модели деформирования при Гэ = 0. Из условия локальной сходимости метода простой итерации [4] для модели деформирования Фарриса имеем

В9 I-1 + т

1 -

1

I.

\р

и..

||р

к

I'

< 1.

(10)

< 1,

(11)

КО (Iи )" Я ехр(яа)

где £и — скорость изменения интенсивности деформаций. Видно, что этот метод сходится только при определенных сочетаниях материальных параметров В, р, т, К, О, п, я модели Фарриса, интенсивности деформаций 1и в конечных элементах и скорости её изменения £и, среднего напряжения а и объемной деформации 9. В общем случае сходимости нет, и даже при сколько угодно точном начальном приближении удобный в

I

и

о

р

реализации процесс простых итераций может разойтись. Условия (10), (11) накладывают очень жесткие ограничения на параметры модели, реализуемые напряжения и деформации. В частности, при растяжении (а > 0, 9> 0) метод, как это видно из (11), сходится только в

той области, где напряжения а удовлетворяют неравенству еча<(КО(Iи )пч) . Ширина

этой области определяется малой для полимеров величиной К 1. Неприемлемы также ограничения, накладываемые на параметры т, р при разгрузке. Из условия (10) для простейшего режима деформирования (квазистатический режим «нагрузка и разгрузка» с равными скоростями, t0 - момент начала разгрузки, 2t0- окончание цикла «нагрузка-

разгрузка») и при В = 0 следует, что метод сходится, если, например, в момент t t0 выполняется неравенство тр-1(2р +1) < 1, в момент t 2t0 - неравенство т < 1. Реальные же значения параметров т , р модели Фарриса больше единицы ( т > р > 1).

При использовании метода Ньютона, исходя из ниже следующей оценки допустимого промаха е< 2М1М-1 для скалярной функции /(х) = 0 (|/ (х)| > М1, |/"(х)| < М2)

х е (х* -е, х* + е) [6], для модели Фарриса имеем следующие оценки допустимого промаха начального приближения относительно элементных значений интенсивности деформаций I и и среднего напряжения а

е, < 2/1/2-1

(12)

где

/ = В91-1 + т

1-

V

"р

+1.

-1 <н

и

р +1

/ = 1-71 (1 + /1)-В61 -2 - тр-11

( п

О (1 и )" ч2а ехР(^а)

р

"р у

р

"р у

8а <

1 + -

а.

(13)

2 (К 1 + О (1и )пч ехр(ча)

Из (12), (13) следует, что е -окрестность, в зависимости от значений материальных параметров модели, реализуемого напряженно-деформированного состояния и его предыстории, может быть очень узкой. Например, в областях с сопоставимым вкладом в

объемную деформацию механизмов дилатации и дилатансии К-1а ~ О (I и )п ехр(ча) и при

а ^ -4-1, имеем еа ^ 0. При переходе от нагрузки к разгрузке ^ t0) для вышеуказанного

режима деформирования при т = 2 , р = 20 допустимый промах е£ составляет величину не

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

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

I

1

и

и

I

и

и

I

I

и

и

I

I

и

и

Коши. Центральным моментом, при этом, является выбор подходящего, органически вписывающегося в "физику" решаемой задачи, линейного оператора, с помощью которого находятся последовательные приближения. В механике деформируемого твердого тела очевидным и естественным вариантом является выбор изотропного оператора линейной теории упругости Ду = 2^3^. + K9g . Различные реализации этого подхода заключаются далее в том или ином выборе коэффициентов Ламе д, K, являющихся итерационными параметрами задачи. Наиболее теоретически разработаны и практически опробованы идеи этого подхода в задачах теории пластичности. Основа была заложена в работах А.А. Ильюшина (см., например, [7]), предложившего метод упругих решений. Сходимость метода в классе обобщенных решений при активном деформировании впервые доказана в [8]. Развитие метода для решения задач деформирования с достаточно общими нелинейными определяющими соотношениями получило в работах Б.Е. Победри (см., например, [9]). В дальнейшем, частях 2 и 3 данной публикации, воспользовавшись идеями указанного подхода и метода искусственной сжимаемости, а также в продолжение идей, изложенных в [10], будут построены всегда сходящиеся методы решения систем нелинейных проекционно-сеточных уравнений большого порядка рассматриваемых квазистатических задач неупругого деформирования полимерных изделий сложной формы при немонотонных режимах нагружения с учетом эффектов внутренней вязкости материала и микроструктурных повреждений, нелинейного объемного поведения.

СПИСОК ЛИТЕРАТУРЫ

1. Альес М.Ю. Феноменологические описания нелинейного сопротивления деформированию высоконаполненных эластомерных (нано) композитов. Часть 1. Малые деформации // Химическая физика и мезоскопия. 2010. Т. 12, № 1. С. 69-77.

2. Самарский А.А., Гулин А.В. Численные методы. М. : Наука, 1989. 432 с.

3. Оден Дж. Конечные элементы в нелинейной механике сплошных сред. М. : Мир, 1976. 464 с.

4. Ортега Дж., Рейнболдт В. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М. : Мир, 1975. 560 с.

5. Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Л. : Гидрометеоиздат, 1986. 352 с.

6. Самарский А.А. Введение в численные методы. М. : Наука, 1987. 288 с.

7. Ильюшин А.А. Пластичность. Основы общей математической теории. М. : Изд-во АН СССР, 1963. 272 с.

8. Ворович И.И., Красовский Ю.П. О методе упругих решений // Доклады АН СССР. 1959. Т. 126, № 4. С. 740-743.

9. Победря Б.Е. О сходимости метода «упругих решений» в нелинейной вязкоупругости // Доклады АН СССР. 1970. Т. 195, вып. 26. С. 307-310.

10. Альес М.Ю. Эволюционные схемы численного решения нелинейных краевых задач неньютоновских ползущих течений // Химическая физика и мезоскопия. 2013. Т. 15, № 2. С. 193-200.

EVOLUTIONARY ISOTROPIC NUMERICAL SCHEMES FOR ACTION BOUNDARY VALUE PROBLEMS OF NONLINEAR QUASI-STATIC DEFORMATION. PART 1. FEATURES OF THE CONVERGENCE OF THE NUMERICAL SOLUTIONS OF THE NONLINEAR TERMS IN THE SMALL DEFORMATIONS

Alies M.Yu.

Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia

SUMMARY. The features of solutions of systems of projection-difference equations of the inelastic behavior of polymer products in small strains, especially the voltage and calculate the average relative change in the volume and the solutions of nonlinear systems of finite element equations by iterative methods.

KEYWORDS: elastomer composites, the nonlinear mechanical behavior, boundary value problems, finite element method, particularly the convergence of iterative methods.

Альес Михаил Юрьевич, доктор физико-математических наук, профессор, главный научный сотрудник ИМ Уро РАН, тел. 89128563824, е-таг7:айе5ту@,таг7.ги