CC BY
57
УДК 539.3
О решении квазистатических задач микро- и мезомеханики в динамической постановке
В.А. Романова1, P.P. Балохонов1, Е.Е. Батухтина1, Е.С. Емельянова2, М.В. Сергеев2
1 Институт физики прочности и материаловедения СО РАН, Томск, 634055, Россия 2 Национальный исследовательский Томский государственный университет, Томск, 634050, Россия
При моделировании характерных процессов на мезоуровне расчетная область содержит большое количество структурных элементов (зерен, включений, пор и т.д.), которые необходимо аппроксимировать расчетной сеткой с достаточной степенью подробности. Известно, что требования к вычислительным мощностям нелинейно возрастают с увеличением количества расчетных элементов. В связи с этим актуальной проблемой при решении задач микромеханики является минимизация вычислительных затрат без потери информативности и точности решения. Существенно минимизировать вычислительные затраты позволяет решение задач квазистатики в динамической постановке. В настоящей работе исследованы условия применимости динамических подходов для моделирования процессов квазистатического деформирования в задачах микро- и мезомеханики, где микроструктура рассматривается в явном виде. Конечно-элементные расчеты квазистатического растяжения в динамической и статической постановках задачи проводились для модельных материалов, нечувствительных к скорости деформации. Показано, что основным параметром, влияющим на совпадение решения динамической задачи со статическим решением, является время, за которое скорость деформирования достигает амплитудного значения. При условии линейного возрастания скорости в течение времени, за которое упругая волна проходит через образец более двух раз, погрешность отклонения динамического решения от статического не превышает 0.1 %. При этом требуемая оперативная память и время вычислений статической задачи превышают на порядок аналогичные затраты на расчеты в динамической постановке. Таким образом, на примере модельных структур показана возможность и целесообразность использования явных динамических методов для решения квазистатических задач микро- и мезомеханики.
Ключевые слова: мезомеханика, микроструктурные модели, численное моделирование, квазистатические процессы, динамические задачи, вычислительные затраты
DOI 10.24411/1683-805X-2018-12007
On the solution of quasi-static micro- and mesomechanical problems
in a dynamic formulation
V.A. Romanova1, R.R. Balokhonov1, E.E. Batukhtina1, E.S. Emelyanova2, and M.V. Sergeev2
1 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634055, Russia 2 National Research Tomsk State University, Tomsk, 634050, Russia
The modeling of characteristic processes at the mesoscale involves the use of a computational domain containing a large number of structural elements (grains, inclusions, voids, etc.), which must be approximated by a detailed computational mesh. The requirements for computational power are known to increase nonlinearly with the increase in the number of calculated elements. For micromechanical problems, this fact makes it necessary to minimize computational costs without loss of information and solution accuracy. The computational costs can be substantially minimized by solving quasi-static problems in a dynamic formulation. In this paper, we study the conditions for the applicability of dynamic approaches for simulating quasi-static deformation in micro- and mesomechanical problems, in which the microstructure is considered explicitly. Finite element calculations of quasi-static tension in the dynamic and static formulation of the problem are carried out for strain-rate insensitive model materials. It has been shown that the coincidence between the dynamic and static solutions is mainly determined by the time during which the strain rate increases up to the amplitude value. If the strain rate increases linearly during the time when the elastic wave passes through the specimen more than two times, the error in the deviation of the dynamic solution from the static solution does not exceed 0.1%. In so doing, the required operational memory and computational time of the static problem exceed by an order of magnitude the similar costs for dynamic calculations. Thus, we have shown using model structures as an example that explicit dynamic methods can be effectively applied for solving quasi-static micro- and mesomechanical problems.
Keywords: mesomechanics, microstructural models, numerical modeling, quasi-static processes, dynamic problems, computational costs
© Романова В.А., Балохонов P.P., Батухтина Е.Е., Емельянова Е.С., Сергеев М.В., 2018
1. Введение
В рамках современных представлений механики деформируемого твердого тела, физики пластичности и материаловедения важнейшая роль в формировании деформационного поведения материалов под нагрузкой отводится микроструктуре [1-6]. Экспериментальные и теоретические данные свидетельствуют о том, что в микрообьемах материала, сопоставимых по размеру с масштабом микроструктурных элементов (например зерен, включений, микропор и т.д.), напряженно-де формированное состояние характеризуется существенной неоднородностью. Локальные характеристики напряженно-деформированного состояния могут на порядок отличаться от среднего уровня. Именно в местах концентрации напряжений происходит накопление необратимой деформации и зарождение микротрещин.
Одним из перспективных подходов, позволяющих исследовать механизмы деформации и разрушения на микро- и мезоуровнях, является численное моделирование с явным учетом микроструктуры [4-10]. В общем случае в пределах структурного элемента свойства материала считаются постоянными и изменяются при переходе через границу раздела. Решение таких задач предьявляет высокие требования к вычислительным мощностям. С одной стороны, для реалистичного воспроизведения деформационных процессов на микро- и мезоуровнях рассматриваемый микрообьем должен содержать достаточное количество структурных элементов. С другой стороны, структурные элементы и приграничные области должны быть аппроксимированы с достаточной степенью подробности для обеспечения приемлемой точности решения задачи. Еще более высокие требования к разрешению расчетных сеток предьявляют задачи о зарождении и распространении трещин, где для строгости и единообразия описания структурных элементов разной формы и обьема рекомендуется применять регулярные сетки. Необходимость использования подробных сеток с большим количеством элементов делает актуальной проблему минимизации вычислительных затрат без потери информативности и точности решения.
Одним из подходов к решению задач квазистатики, позволяющих существенно минимизировать требования к оперативной памяти, дисковому пространству и быстродействию, является замена уравнений равновесия на уравнения движения [11-17]. Таким образом может быть осуществлен переход от неявных схем интегрирования к явным, имеющим существенные преимущества с точки зрения вычислительных затрат. Явные схемы первоначально были разработаны для моделирования процессов высокоскоростной деформации и ударного взаимодействия, где инерционные члены играют определяющую роль [13, 18, 19]. Расчет изменения напряженно-деформированного состояния производится на каждом шаге по времени, при этом механическое
возмущение (волна) передается от элемента к элементу. Условие устойчивости накладывает ограничения на величину шага по времени таким образом, чтобы расстояние, пройденное волной за один шаг по времени, не превысило минимальное расстояние между соседними узлами сетки. Скорость распространения волн малой амплитуды близка к скорости звука. Шаг по времени, обеспечивающий выполнение условия устойчивости, достаточно мал и для моделирования квазистатических процессов требуется большое количество шагов по времени.
Неявные схемы интегрирования в отличие от явных являются безусловно устойчивыми. При решении статических линейных задач точность решения не зависит от количества расчетных шагов, за которое достигается заданная деформация, т.е. конечное деформированное состояние может быть достигнуто за один шаг. При решении нелинейных задач число шагов существенно увеличивается, что связано с необходимостью описать области нелинейного изменения напряженно-деформированного состояния во времени и пространстве с достаточной степенью подробности. Нелинейность задачи может быть обусловлена различными факторами, включая нелинейный отклик материала (например пластичность), геометрическую нелинейность расчетной области, нелинейные граничные условия и др. Кроме того, изменение напряженно-деформированного состояния от начального к конечному состоянию определяется путем итераций, число которых также зависит от степени нелинейности задачи. Таким образом, при решении нелинейных задач преимущества неявных методов, связанные с безусловной устойчивостью, во многих случаях практически исчезают.
Известно, что решения задач в динамической и статической постановках будут сходиться при определенных условиях [16, 17]. В общем случае это условие минимизации волновых эффектов в динамической постановке путем плавного приложения нагрузки и использование моделей материалов, не чувствительных к скорости нагружения. Вместе с тем целенаправленного исследования этих условий для моделирования квазистатического нагружения в материалах с внутренними границами раздела не проводилось. Этим обусловлено ограниченное применение динамических подходов для моделирования квазистатического нагружения в задачах микро- и мезомеханики.
Настоящая работа посвящена исследованию условий применимости динамических подходов для моделирования процессов квазистатического деформирования в микрообьемах структурно-неоднородных материалов. Расчеты одноосного растяжения в динамической и квазистатической постановках проводились для идеализированной двухкомпонентной модели с прямолинейными границами раздела и модели с поликристаллической структурой.
2. Статические и динамические краевые задачи упругопластического деформирования
Общая система уравнений в динамической постановке в декартовой системе координат, описывающая движение точек деформируемого твердого тела в полях внешних воздействий без учета массовых сил, включает уравнения движения
РХ = °Ц, j, (1)
уравнение неразрывности V
V
--x,-,- = 0
и соотношения для скоростей полных деформаций
е - = Ь, - + )' (3)
где — координаты; V = р/ р0 — относительный обь-ем; р0 и р—текущая и начальная плотность материала; о- и е- — компоненты тензоров напряжений и полных деформаций; точка над символом обозначает материальную производную по времени, запятая перед нижним индексом — производную по координатам, i = 1...3.
В общем случае относительное движение точек сплошной среды может быть разложено на деформацию и жесткий поворот [16, 20]. С точки зрения описания поведения материала основной интерес представляет поле деформаций, поскольку оно связано определяющими соотношениями с полем напряжений, в то время как жесткий поворот не вызывает изменения напряженного состояния. В случае малых деформаций для корректировки напряженного состояния с учетом поправки на поворот используется производная Яуманна [13, 16].
Статическая постановка задачи предполагает, что система находится в состоянии равновесия. При отсутствии массовых сил уравнения равновесия имеют вид
о-- = 0. (4)
В случае малых упругих деформаций и малых поворотов, что является характерным для металлов, справедливы соотношения Коши
Ev = + Uji ).
(5)
Как в динамической, так и в статической постановке задачи для замыкания систем уравнений (1)-(3) и (4), (5) требуется сформулировать определяющие соотношения, задающие связь между компонентами тензоров напряжений и деформаций. В общем случае упруго-пластической среды в качестве определяющих соотношений используется обобщенный закон Гука. Для замыкания системы уравнений статики (4), (5) закон Гука записывается в виде
<ij = Cijkl ЕИ,
(6)
где е- — компоненты тензора упругих деформаций; С-м — матрица упругих модулей. В динамической постановке задачи обобщенный закон Гука формули-
руется в скоростной форме:
°у = СуЫ & к1 ■ (7)
Для упругопластической среды тензор полных деформаций удобно представить в виде суммы упругой е- и пластической ер составляющих:
е =ee. +eP j j j •
С учетом (8) закон Гука (6) примет вид
aij = Cijkî (еи Eld) (2) или в скоростной форме
а ii
= Cijklkl Е pl).
(8) (9) (10)
Для определения компонент тензора пластических деформаций построено большое количество различных моделей, выбор которых зависит от материала, условий нагружения, характерных процессов, протекающих в материале, и т.д.
Граничные условия, задающие внешнее воздействие, могут быть либо кинематическими, когда на поверхности тела заданы перемещения или скорости, либо статическими, когда на поверхности тела заданы поверхностные силы. Граничные условия в напряжениях имеют вид
aijnj I о
V j Ра
= т, (11)
где Sо — поверхность или ее часть; Т и п- — компоненты векторов поверхностных сил и нормали к поверхности Sо соответственно. Кинематические граничные условия в перемещениях записываются в виде
и-1^ = и*, (12)
*
где и — компоненты вектора перемещений, заданные на поверхности Би ■
Кинематические граничные условия в случае динамической постановки задачи задаются в скоростях:
й^ = V*, (13)
где V* — заданные значения компонент вектора скорости на поверхности Sv ■ Как правило, кинематические граничные условия отражают условия закрепления или движения точек поверхности относительно заданного направления.
Очевидно, что в условиях развития неравновесных динамических процессов, например резонансных явлений, поверхностных волн и т.д. [18], возникает вопрос о корректности решения квазистатической задачи в динамической постановке. Физическая постановка задач микромеханики предполагает, что исследуемый микрообьем «вырезается» из макроскопического образца в форме прямоугольного бруска или куба, линейные размеры которого по трем направлениям отличаются несущественно. Выбор такой геометрии исключает возникновение резонансных явлений и волновых эффектов, которые могут возникать при определенной комби-
нации геометрических параметров, свойств материалов и условий нагружения.
Вопрос о формулировке граничных условий, наиболее реалистично отражающих условия нагружения микроструктур, на сегодняшний день остается открытым. Даже для микрообъемов, выходящих на поверхность, задание реальных условий нагружения представляется сложно реализуемой задачей. Как правило, в расчетах задаются идеализированные граничные условия. Для микрообъемов, выходящих на свободную поверхность образца, на соответствующей грани задаются условия отсутствия внешних сил. На поверхностях нагружения формулируются граничные условия, задающие растяжение, сжатие, сдвиг и т.п. На поверхностях микрообъемов, находящихся внутри образца, часто задают условия симметрии или периодические граничные условия, имитирующие влияние окружающего материала, или условия свободных поверхностей.
3. Особенности численного решения динамических и статических задач
Методы численного интегрирования дифференциальных уравнений могут быть условно разделены на явные и неявные в зависимости от выбора алгебраических соотношений, аппроксимирующих точное решение системы или производные, входящие в исходную систему [16]. Решение статической задачи (4)-(6), (11), (12) возможно только в рамках неявных подходов. Приближенные алгебраические соотношения, построенные для системы дифференциальных уравнений в рамках неявных подходов, предполагают нахождение значений искомых функций через функции, также неизвестные на данном шаге интегрирования. Такие схемы являются безусловно устойчивыми, но в численной реализации требуют привлечения итерационных методов. Кроме того, моделирование эволюции напряженно-деформированного состояния в области нелинейности и быстрого изменения параметров может потребовать существенного увеличения расчетных шагов.
Динамическая задача (1)-(3) может быть решена явными методами [13, 16]. В рамках явных подходов значения искомых функций на данном шаге интегрирования определяются через значения функций, известные с предыдущего шага. Такие методы являются условно устойчивыми. Условие устойчивости накладывает существенные ограничения на величину шага интегрирования по времени. При решении динамической задачи шаг интегрирования по времени связан со скоростью распространения упругих волн в материале. Это наиболее быстропротекающий физический процесс, описываемый уравнениями (1)-(3). Условие устойчивости задает связь между шагом интегрирования по времени и наименьшим шагом по пространству в виде
Дt = ^ (14)
С1
где ^ < 1 — число Куранта; — минимальный шаг расчетной сетки; С1 — скорость продольной звуковой волны. Согласно (14), упругое возмущение за один шаг по времени не должно проходить расстояние, превышающее минимальный шаг по пространству. Очевидно, что расчет квазистатического нагружения с реальными скоростями деформации 10-5—10—4 с-1 в рамках явных подходов является нецелесообразным, поскольку требует огромного количества шагов по времени.
Проведем сравнительный анализ вычислительных затрат при решении динамической и статической задач методом конечных элементов [16]. Основными характеристиками вычислительных затрат при решении задачи численными методами являются требуемая оперативная память, дисковое пространство и время вычислений. В общем случае эти характеристики зависят от целого ряда факторов, включая количество искомых функций, метод решения и размерность задачи, разрешение расчетной сетки и тип элементов. Однозначно определить такую зависимость довольно сложно, однако приблизительные оценки могут быть приведены на основании имеющихся литературных данных [16]. Известно, что оперативная память и время расчета при решении задач явными методами пропорционально возрастают с увеличением числа элементов. Кроме того, время расчета обратно пропорционально наименьшему размеру элементов. Например, в случае трехмерной конечно-элементной модели с кубическими элементами увеличение числа элементов по каждому из трех направлений в два раза при сохранении геометрических размеров расчетной области приведет к увеличению времени расчета в 16 раз в результате увеличения числа элементов в 8 раз и уменьшения размера элементов в два раза. Дисковое пространство и оперативная память не зависят от размера элементов, поэтому эти параметры возрастут только в 8 раз.
Гораздо сложнее предсказать, насколько изменится производительность вычислений при использовании неявных методов. По грубым оценкам для многих задач вычислительные требования возрастают пропорционально квадрату количества степеней свободы. Количество степеней свободы, в свою очередь, определяется типом и количеством элементов, а также постановкой и размерностью задачи. Так, например, каждый узел кубического элемента имеет 3 степени свободы, связанные со смещениями по трем координатам. Увеличение количества элементов в трехмерной конечно-элементной модели с кубическими элементами по трем направлениям в 2 раза приводит к увеличению количества степеней свободы примерно в 23, при этом вычислительные требования (время расчета, дисковое пространство и оперативная память) возрастают примерно в (23)2, т.е. в 64 раза.
Точные оценки необходимых вычислительных ресурсов можно дать только в каждом конкретном случае.
Для примера на рис. 1 приведено сравнение вычислительных затрат на решение упругой задачи в динамической и статической постановках на персональном компьютере с 4-ядерным процессором Intel i7 и 16 ГБ оперативной памяти с использованием конечно-элементного пакета ABAQUS/Standard и ABAQUS/Explicit (6.12) [16]. В первой серии расчетов оценивалось максимально возможное количество элементов, допустимых при данном объеме оперативной памяти компьютера, при решении задачи с использованием от одного до четырех процессоров. Геометрическая модель представляла собой кубический образец, аппроксимированный регулярной сеткой с линейными кубическими элементами (C3D8R [16]). Проведенные оценки показали, что максимальное количество элементов Nmax, которое может быть задано в динамической постановке, при прочих равных условиях (размерность задачи, модель изотропного упругого материала, условия нагружения, параметры вывода) в 8-10 раз превышает значение Nmax в статических расчетах.
В обоих случаях Nmax нелинейно увеличивается с увеличением количества ядер (рис. 1, а). Так, использование двух ядер позволяет увеличить Nmax в статической задаче на 63 %, а в динамической — на 72 % по сравнению с расчетами без параллельных вычислений. С увеличением количества ядер эта разница нелинейно уменьшается. Например, для трех и четырех ядер разница в Nmax составляет 16 % в динамических и 11 % в статических расчетах.
Для сравнения временных затрат на решение задачи в динамической и статической постановках необходимо определить условия, при которых эти решения совпадут. Рассмотрим задачу одноосного растяжения изотропного линейно-упругого материала. Такая задача имеет аналитическое решение, с которым будем сравнивать результаты моделирования в динамике и статике. При одноосном растяжении стержня только одна компонента тензо-
ра напряжений, направленная вдоль оси растяжения, отлична от нуля и определяется законом Гука:
стп = Ееп, (15)
где Е — модуль Юнга.
Конечно-элементная модель представляла собой прямоугольный образец 5x5x5 см3, аппроксимированный регулярной сеткой с кубическими элементами. Модули упругости, заданные в расчетах, соответствовали характеристикам поликристаллического титана: Е = 112 ГПа, V = 0.32. Оценки показали, что характеристики напряженно-деформированного состояния, полученные численно из решения задачи в статической постановке, совпадают с аналитическими значениями с относительной погрешностью ~10-5.
Характеристики напряженно-деформированного состояния, полученные в статическом и динамическом расчетах, сравнивались при деформации ~0.001. В статической постановке такая деформация достигается за один расчетный шаг. Для достижения заданной деформации в динамической постановке необходимо сделать достаточно большое количество шагов по времени. При этом каждое смещение узлов на поверхности расчетной области генерирует упругую волну, амплитуда которой напрямую зависит от амплитуды заданных смещений. Граничные условия в динамической постановке задачи формулируются в скоростях (13). Минимизация волновых эффектов, неизбежно присутствующих в динамике, возможна путем плавного наращивания скорости растяжения до амплитудного значения V*. Инерционные члены в уравнении движения (1) связаны с ускорением. Чем больше время наращивания скорости t*, тем меньше ускорение и, соответственно, меньше проявляются динамические эффекты. После достижения заданной амплитуды процесс нагружения продолжается при постоянной скорости. В расчетах скорость наращивалась до амплитудного значения V* по линейному закону в течение времени t * (рис. 2, а).
0.04
12 3 4
Количество параллельных процессов
Статика Динамика
-1-1-1-1-г
12 3 4
Количество параллельных процессов
Рис. 1. Зависимость максимально возможного количества элементов (а) и времени вычислений (б) динамической и статической задач от количества параллельных процессов
0.10-
с Ё Ь
I
| 0.05 Н
в ь
0.00-
Время
0
Тб
Скорость деформации, с 1 □ 104 о 102 ■ 10"1 • 10"3
Аппроксимационная кривая
10
20 30
МКС
40
Рис. 2. Схема задания скорости нагружения в динамических расчетах (а) и отклонение экстремальных значений полей напряжений, полученных в динамических расчетах, от среднего уровня, соответствующего статическому решению, в зависимости от амплитудных значений скорости V* и времени наращивания скорости (б)
В результате волновых эффектов напряженно-деформированное состояние при моделировании одноосного растяжения в динамической постановке является неоднородным (рис. 3, б) — все компоненты тензоров напряжений и деформаций отклоняются от равновесного состояния в некотором диапазоне. При этом средний уровень соответствует статическому значению (рис. 3, а). Очевидно, что чем меньше диапазон отклонения, тем меньше влияние волновых эффектов и тем более точно решение динамической задачи приближается к решению статической задачи. Для определения параметров, влияющих на точность совпадения динамического и статического решений, была проведена серия динамических расчетов, в которых варьировались амплитуда V * и время наращивания скорости растяжения £ ■
Зависимость отклонения экстремальных значений полей напряжений, полученных в динамических расчетах с разными V* и £*, от статического значения приведена на рис. 2, б. Анализ результатов показал, что
основным параметром, влияющим на минимизацию волновых эффектов, является время наращивания скорости до амплитудных значений. Точки на рис. 2, б, полученные при одних и тех же значениях £ * и разных амплитудных значениях скорости растяжения V*, практически совпадают. В то же время диапазон отклонения экстремальных значений напряжений от среднего уровня уменьшается экспоненциально с увеличением £*■ Сопоставление зависимости, приведенной на рис. 2, б, с геометрическими параметрами образца и упругими характеристиками материала, определяющими скорость распространения упругих возмущений, показало, что для значений £ *, при которых упругая волна дважды проходит через образец, волновые эффекты существенно снижаются и относительная погрешность отклонения экстремальных значений от среднего не превышает 0.01 %. Таким образом, приемлемое совпадение статического и динамического решений с относительной погрешностью порядка 0.01 % достигается при выполнении условия
о
Рис. 3. Поля эквивалентных напряжений, полученные в статическом (а) и динамическом (б) расчетах при одноосном растяжении
где L — длина образца вдоль оси растяжения; С — скорость продольной звуковой волны в материале.
Для оценки временных затрат на решение задачи в динамической и статической постановках одни и те же расчеты были проведены на одном и том же компьютере с использованием от 1 до 4 ядер. Время наращивания скорости до амплитудного значения в динамических расчетах было выбрано из условия (16). Полученные зависимости приведены на рис. 1, б. Сравнение показало, что время, затраченное на расчет задачи в статической постановке, в 7-10 раз превышает время расчета динамическим методом. При этом выполнение условия (16) обеспечивает совпадение динамического и статического решений с относительной погрешностью ~0.01 %.
Проведенные оценки убедительно показывают целесообразность решения задач квазистатики в рамках динамического подхода. В обоих случаях время расчетов уменьшается с увеличением количества ядер. Однако если использование двух ядер обеспечивает почти двукратный выигрыш по времени по сравнению с расчетами на одном ядре, то разница для трех и четырех ядер составляет 22 % в динамических и 13 % в статических расчетах. Такая нелинейная зависимость связана с увеличением времени обмена информацией между параллельными процессами при увеличении количества ядер.
4. Моделирование квазистатического нагружения структурно-неоднородных материалов в рамках динамического подхода
4.1. Квазистатическое растяжение двухкомпонент-ной структуры с прямолинейными границами раздела
Границы раздела между разнородными средами являются источниками отражения и преломления упругих
волн, неизбежно возникающих при решении динамической задачи. Если в однородном материале границами раздела являются только поверхности образца, то в неоднородном материале все внутренние границы раздела влияют на изменение волновой динамики в той или иной степени. Вопрос о влиянии внутренних границ раздела на решение задачи квазистатики в динамической постановке был исследован на примере двухком-понентной модели 5x5x5 см3, состоящей из двух изотропных упругих материалов, разделенных прямолинейными границами (рис. 4, а). По плотности и упругим характеристикам материалы соответствовали титану и алюминию.
В динамической постановке скорость растяжения наращивалась до амплитудного значения V* в течение времени t*, а затем поддерживалась постоянной (рис. 2, а). Амплитудное значение скорости движения узлов составляло 2.5 м/с, что соответствовало скорости деформации 102 с-1. Время нарастания скорости С варьировалось в расчетах от 5 • 10-6 до 5 • 10-5 с. Решения задач в динамической и квазистатической постановках сравнивались при деформации 0.125 %.
В отличие от случая одноосного растяжения однородного материала все компоненты тензоров напряжений и деформаций в двухкомпонентной структуре отличны от нуля. В качестве иллюстрации на рис. 4, б, в приведены распределения компонент тензора напряжений а22 и ст33, полученные в квазистатическом расчете. В однородном материале в условиях одноосного растяжения вдоль оси Х1 эти напряжения равны нулю. В двухкомпонентной структуре напряжения а22 и ст33 отличны от нуля и демонстрируют как положительные, так и отрицательные значения в диапазоне ±13 МПа. Максимальные ненулевые напряжения наблюдаются вблизи границ раздела. Следует отметить, что средние по объему значения а22 и ст33 равны нулю, что соответствует равновесному состоянию.
0
СТ33
Рис. 4. Модельная двухкомпонентная структура (а) и распределение напряжений о22 (б) и о33 (в), полученные в квазистатических расчетах при растяжении, е = 0.125 %
оеф МПа Я 150 - 125 100 75
Рис. 5. Распределение эквивалентных напряжений в модельном образце с прямолинейными границами раздела при растяжении до е = 0.125 % для динамической (а, б) и статической (в) постановки задачи. Время нагружения в динамической задаче £ * = = 5 • 10-6 (а) и 5-10"5 с (б)
Основной вклад в эквивалентные напряжения вносит о11? действующее в направлении растяжения и принимающее значения от 75 до 150 МПа при деформации 0.125 %. Распределение напряжений для динамического и статического расчета приведены на рис. 5. Сравнительный анализ показал, что приближение решения динамической задачи к статическому решению существенно зависит от значения £ *■ Так, при £ * = 5 • 10 5 с относительная погрешность отклонения статического и динамического решений не превышает 0.05 % (ср. рис. 5, б и в). За это время упругие волны проходят через образец более 8 раз. При уменьшении времени нагружения до £ * = 5 • 106 с волновые эффекты существенно влияют на картину распределения напряжений (ср. рис. 5, а и в). При таких условиях моделируемый процесс не может рассматриваться как квазистатический. Динамические расчеты, проведенные для разных значений показали, что для модели с прямолинейными границами раздела совпадение статического и динамического решений с относительной погрешностью, не превышающей 0.05 %, достигается при выполнении условия (16), как и в случае однородного материала.
4.2. Моделирование квазистатического растяжения поликристаллической структуры
Поскольку кривизна границ раздела является одним из определяющих факторов, влияющих на значения локальных характеристик напряженно-деформированного состояния, важную роль в задачах микро- и мезо-механики играет реалистичное воспроизведение геометрии структурных элементов. В связи с этим анализ условий применимости динамического подхода для моделирования квазистатического растяжения был проведен для поликристаллической модели, представленной на рис. 6, а. Генерация трехмерной структуры осуществлялась методом пошагового заполнения [21]. В качестве начальных условий в обьеме, дискретизиро-ванном регулярной сеткой размером 70 х 70 х 70 с шагом 10 мкм, распределялись центры зарождения, координаты которых определялись случайным образом с помощью генератора случайных чисел. Для всех зерен задавались сферический закон и одинаковая скорость роста.
По определению статическая формулировка задачи означает, что материал находится в равновесии, т.е. внешние приложенные силы и внутренние силы, выз-
100 80
<3 60 С
= 40
20 0
та
- Статика, Ас = 0.00375
■ Статика, Ас = 0.00075
■ Динамика
0.00 0.02 0.04 0.06
Рис. 6. Модель поликристалла (а), распределение эквивалентных напряжений при растяжении до 7.5 % (б) и о-е-диаграммы, полученные в квазистатических и динамическом расчетах (в)
ванные напряжениями в материале, уравновешивают друг друга. В случае вязкоупругого или вязкопластичес-кого материала напряжения зависят от скорости деформации. В этом случае баланс между приложенной нагрузкой и откликом материала достигается за определенное время, зависящее от скорости нагружения и вязких свойств материала. Очевидно, что одним из условий применимости динамического подхода для моделирования квазистатических процессов является нечувствительность материала к скорости нагружения. Только в этом случае динамическая постановка задачи, где скорости приложенной деформации на несколько порядков выше скоростей квазистатических процессов, обеспечит решение, совпадающее со статическим.
В настоящей работе для описания упругопластичес-кого отклика зерен использовалась модель упругоплас-тического изотропного материала с линейным изотропным упрочнением [13]. В рамках этой модели упруго-пластический переход описывается критерием текучести Мизеса, согласно которому компоненты тензора пластических деформаций в (9) и (10) равны нулю и напряжения определяются по упругому закону при условии
—Sj.Sj,- < а.
2 J 'J
(17)
где Sj — компоненты девиатора напряжений; а = = ст0 + _Depq — предел текучести с учетом деформационного упрочнения. В области пластического деформирования компоненты тензора напряжений сносятся по нормали на поверхность текучести. Упругопласти-ческие характеристики зерен варьировались в пределах 20 % относительно средних значений: E = 112 ГПа, v = = 0.32, а0 = 100 МПа [22]. Коэффициент деформационного упрочнения D для всех зерен равнялся 50 МПа.
Граничные условия задавали одноосное растяжение вдоль оси Xj (рис. 6, а). В динамической постановке задачи скорость движения поверхностей нагружения линейно наращивалась до амплитудного значения v*, а затем поддерживалась постоянной (рис. 2, а). С учетом условия (16) время нарастания нагрузки С задавалось равным 2 • 10-6 с, что более чем в десять раз превышает время, за которое упругая волна проходит через образец. Амплитудное значение скорости деформации в динамическом расчете составляло 1.4 •Ю2 с-1. В статической постановке задачи граничные условия на поверхностях нагружения задавались в смещениях вдоль оси Xj. Напряженно-деформированное состояние, полученное в динамической и статической постановках, сравнивалось при деформации 7.5 %. В квазистатических расчетах варьировалось количество расчетных шагов, за которое достигалась заданная степень деформации 7.5 %. В первом случае количество шагов равнялось 20, а величина приращения деформации за один шаг состав-
ляла Де = 0.00375. Во втором случае та же деформация достигалась за 100 шагов величиной Де = 0.00075.
Как и в случае двухкомпонентной модели с прямолинейными границами раздела, все компоненты тензоров напряжений и деформаций в поликристаллической структуре отличны от нуля уже на упругой стадии нагру-жения и вносят вклад в формирование деформационного отклика материала. Однако наибольший интерес представляет анализ инвариантных характеристик напряженно-деформированного состояния. Проведем сравнительный анализ эквивалентных напряжений и деформаций, полученных в динамической и статической постановках задачи. Поля напряжений во всех трех расчетах качественно практически совпадают. В качестве иллюстрации распределение эквивалентных напряжений, характерное для всех случаев, приведено на рис. 6, б. Также практически полностью совпадают осредненные кривые нагружения (рис. 6, в). По оси абсцисс отложена деформация образца
е = -
L - Ln
(18)
где L и L0 — начальная и конечная длина образца вдоль оси растяжения. По оси ординат отложена величина эквивалентных напряжений, осредненных по всем расчетным элементам
(19)
где N — количество элементов в конечно-элементной модели.
Расчеты показали, что существенное влияние на распределение пластических деформаций в статической постановке задачи оказывает величина Де. Этот вывод проиллюстрирован на рис. 7, а, б, где поля эквивалентных пластических деформаций, полученные в статических расчетах при значениях Де = 0.00375 и 0.00075, приведены при одной и той же степени деформации образца. Приведенные распределения пластических деформаций (рис. 7, а, б) демонстрируют существенные качественные и количественные отличия. Очевидно, что чем меньше задается Де в области пластического течения, тем более точное решение получается в статической постановке. При этом важно отметить, что распределения пластических деформаций, полученные в динамической постановке задачи (рис. 7, в), согласуются со статическими расчетами для Де = 0.00075 и качественно, и количественно (ср. рис. 7, б и в).
Для более точного количественного сравнения на рис. 8 и 9 приведены кривые отклонения друг от друга соответствующих характеристик напряженно-деформированного состояния, полученных в разных расчетах. По оси абсцисс отложена величина относительной погрешности
0.03
Рис. 7. Распределение пластической деформации в модельном поликристалле при растяжении до 7.5 % для квазистатической постановки при Де = 0.00375 (а) и 0.00075 (б) и динамической постановки (в)
бо еч(и) =
п)
п)
-1
бе%(п)=
• 100%,
(п )
(п )
(20)
•100%, п = 1, ..., N
где N — количество конечных элементов в модели, символами * и ** обозначены величины, полученные в разных расчетах. По оси ординат отложена относительная обьемная доля материала, демонстрирующая соответствующие значения бо^ и бе
Сравнение значений пластической деформации (рис. 8, а), полученных в динамическом и статическом расчете с шагом Де = 0.00075, показало, что для подавляющего большинства элементов величина 8ерп не
еч
превышает 0.1 %. Лишь в единичных элементах, главным образом расположенных в областях локализации пластической деформации, эта величина принимает значения от 0.5 до 4.5 %. Разница между эквивалент-
ными напряжениями, полученными в динамике и статике с шагом Де = 0.00075, не превышает сотые доли процента во всей расчетной области (рис. 8, б).
Существенно более сильное расхождение наблюдается между характеристиками напряженно-деформированного состояния, полученными в статических расчетах с разным шагом по деформации Де (рис. 9). Для подавляющего количества элементов погрешность между значениями ерч, полученными в расчетах с разными Де, составляет более 10 % и в единичных локальных областях достигает 80 %. Погрешность отклонения полей напряжений составляет порядка 1 % и в локальных областях достигает ~1.3 %, что для оценки напряжений, ограниченных поверхностью текучести, является существенным.
Проведенные оценки показали, что решение задачи в статической постановке с шагом Де = 0.00375 не обеспечивает достаточной точности описания напряженно-деформированного состояния на микроуровне.
Рис. 8. Сравнение характеристик напряженно-деформированного состояния, полученных при решении задачи о растяжении поликристалла в динамической постановке (а) и в статическом расчете с шагом Де = 0.00075 (б)
20 40 60
8eL %
Рис. 9. Сравнение характеристик напряженно-деформированного состояния, полученных при решении задачи о растяжении поликристалла в статической постановке с шагом Де = 0.00375 (а) и 0.00075 (б)
При этом время, затраченное на вычисления, составляло 3 ч 20 мин. Адекватное описание деформационных процессов на микроуровне было получено при решении квазистатической задачи, где заданная степень деформации достигалась за 100 шагов величиной Де = = 0.00075, а также динамической задачи при выполнении условия (16). При этом при равных вычислительных мощностях расчет задачи в динамической постановке занял пятьдесят минут, а в квазистатической — более 33 ч. Таким образом, для рассмотренного случая поликристаллического материала, не чувствительного к скорости деформации, можно сделать вывод о целесообразности использования динамических подходов и явных методов интегрирования для моделирования квазистатических процессов.
5. Заключение
Проблема минимизации вычислительных затрат является крайне актуальной для решения задач микро- и мезомеханики, где структура материала вводится в рассмотрение в явном виде. Для моделирования характерных процессов на мезоуровне расчетная область должна содержать достаточно большое количество структурных элементов (зерен, включений, пор и т.д.), которые, в свою очередь, должны быть аппроксимированы расчетной сеткой с достаточной степенью подробности. Необходимость использования подробных расчетных сеток обуславливает высокие требования к вычислительным мощностям, в ряде случаев делающих численный анализ нецелесообразным.
Существенно минимизировать вычислительные затраты (требуемую оперативную память, дисковое пространство и время расчетов) позволяет решение квазистатических задач в динамической постановке. В настоящей работе исследованы условия применимости динамических подходов для моделирования процессов квазистатического деформирования в задачах микро-и мезомеханики, где микроструктура рассматривается
в явном виде. Динамические и квазистатические конечно-элементные расчеты одноосного растяжения проводились для однородного материала двухкомпонентной структуры с прямолинейными границами раздела и поликристаллической модели материала, не чувствительного к скорости нагружения.
Анализ результатов расчетов показал, что решение динамической задачи совпадает с высокой точностью со статическим решением при выполнении нескольких условий, одним из которых является минимизация волновых эффектов в динамических расчетах путем плавного наращивания скорости деформирования. При этом основным параметром, влияющим на согласие характеристик напряженно-деформированного состояния в динамических и статических расчетах, является не амплитудное значение скорости деформации, а время его достижения. Показано, что при наращивании скорости деформации в течение времени, за которое упругая волна проходит через образец более двух раз, погрешность отклонения динамического решения от статического составляет менее 0.1 %.
Другим условием применимости динамического подхода для моделирования квазистатического нагружения является выбор модели материала, не чувствительной к скорости деформации. Это означает, что баланс между приложенной нагрузкой и откликом материала достигается за бесконечно малое время при любой скорости нагружения. Это условие накладывает определенные ограничения на класс статических задач, которые можно решать в динамической постановке. За пределами этого класса оказываются задачи вязкоупругости и вяз-копластичности, где реальная реология играет определяющую роль. В этом случае вопрос о применимости динамической постановки для моделирования квазистатических процессов требует отдельного изучения.
Проведенные расчеты показали, что при прочих равных условиях время решения квазистатической задачи и объем требуемой оперативной памяти более чем в
10-40 раз превышали соответствующие вычислительные затраты на решение динамической задачи. При этом отклонении оба решения совпадали с высокой точностью. Таким образом, показана целесообразность использования динамического подхода для моделирования квазистатических процессов в задачах микро- и мезомеханики.
Работа выполнена в рамках государственного задания на 2017-2020 гг. Результаты раздела 4.2 получены в рамках проекта РФФИ № 17-08-00643А.
Литература
1. Панин В.Е., Егорушкин В.Е. Деформируемое твердое тело как нелинейная иерархически организованная система // Физ. мезо-мех. - 2011. - Т. 14. - № 3. - С. 7-26.
2. Поверхностные слои и внутренние границы раздела в гетерогенных материалах / Под ред. В.Е. Панина. - Новосибирск: Изд-во СО РАН, 2006. - 520 с.
3. Панин В.Е., Егорушкин В.Е. Основы физической мезомеханики пластической деформации и разрушения твердых тел как нелинейных иерархически организованных систем // Физ. мезомех. -2015. - Т. 18. - № 5. - С. 100-113.
4. McDowell D.L. A perspective on trends in multiscale plasticity // Int. J. Plasticity. - 2010. - V. 26. - P. 1280-1309.
5. Bieler T.R., Eisenlohr P., Roters F., Kumar D., Mason D.E., Crimp M.A., Raabe D. The role of heterogeneous deformation on damage nucleation at grain boundaries in single phase metals // Int. J. Plasticity. - 2009. - V. 25. - P. 1655-1683.
6. Kovalevskaya Zh.G., Ivanov Yu.F., Perevalova O.B., Klimenov V.A., Uvarkin P. V. Study of microstructure of surface layers of low-carbon steel after turning and ultrasonic finishing // Phys. Met. Metall. -2013. - V. 114. - No. 1. - P. 41-53.
7. Psakhie S., Ovcharenko V., Baohai Yu, Shilko E., Astafurov S., Ivanov Yu., Byeli A., Mokhovikov A. Influence of features of interphase boundaries on mechanical properties and fracture pattern in metal-ceramic composites // J. Mater. Sci. Technol. - 2013. - V 29. -No. 11.- P. 1025-1034.
8. Chawla N., Sidhu R.S., Ganesh V.V. Three-dimensional visualization and microstructure-based modeling of deformation in particle-reinforced composites // Acta Mater. - 2006. - V. 54. - No. 6. - P. 15411548.
9. Diard O., Leclercq S., Rousselier G., Cailletaud G. Evaluation of finite element based analysis of 3D multicrystalline aggregates plasticity. Application to crystal plasticity model identification and the study
of stress and strain fields near grain boundaries // Int. J. Plasticity. -2005.- V. 21. - P. 691-722.
10. Makarov P.V.. Schmauder S., Cherepanov O.I., Smolin Yu.I., Romanova V.A., Balokhonov R.R., Saraev D.Yu., Soppa E., Kizler P., Fischer G., Hu S., Ludwig M. Simulation of elastic plastic deformation and fracture of materials at micro-, meso- and macrolevels // Theor. Appl. Frac. Mech. - 2001. - V. 37. - No. 1-3. - P. 183-244.
11. Naghibi H., Janssen D., Khoshgoftar M., Sprengers A., Semih E., Van Den Boogaard T., Verdonschot N. A comparison between dynamic implicit and explicit finite element simulations of the native knee joint // Med. Eng. Phys. - 2016. - V 38. - P. 1123-1130.
12. Rojek J., Zienkiewicz O.C., Onate E., Postek E. Advances in FE explicit formulation for simulation of metalforming processes // J. Mater. Proc. Tech. - 2001. - V. 119. - P. 41-47.
13. Wilkins M.L. Computer Simulation of Dynamic Phenomena. - Berlin: Springer-Verlag, 1999. - 265 p.
14. Lindgren L.-E., Edberg J. Explicit versus implicit finite element formulation in simulation of rolling // J. Mater. Proc. Tech. - 1990. -V. 24. - P. 85-94.
15. Zhang H., Dong X., Wang Q., Zeng Z. An effective semi-implicit integration scheme for rate dependent crystal plasticity using explicit finite element codes // Comput. Mater. Sci. - 2012. - V 54. - P. 208218.
16. Getting Started with Abaqus: Keywords Edition. - ABAQUS 6.12 PDF Documentation, 2012.
17. Romanova V.A., Soppa E., Schmauder S., Balokhonov V.A. Meso-mechanical analysis of the elasto-plastic behavior of a 3D composite-structure under tension // Comput. Mech. - 2005. - V. 36. - P. 475483.
18. Smolin A.Yu., Konovalenko Ig.S., Psakhie S.G. Identification of elastic waves generated in the contact zone of a friction couple // Tech. Phys. Lett. - 2007. - V. 33. - No. 7. - P. 600-603.
19. Balokhonov R.R., Makarov P.V., Romanova V.A., Smolin I.Yu. Simulation of crystal plasticity under dynamic loading // Comput. Mater. Sci. - 1999. - V. 16. - P. 355-361.
20. Трусов П.В., Швейкин А.И., Янц А.Ю. О разложении движения, независимых от выбора системы отсчета производных и определяющих соотношениях при больших градиентах перемещений: взгляд с позиций многоуровневого моделирования // Физ. мезо-мех. - 2016. - Т. 19. - № 2. - С. 49-65.
21. Romanova V., Balokhonov R., Makarov P., Schmauder S., Soppa E. Simulation of elasto-plastic behavior of an artificial 3D-structure under dynamic loading // Comput. Mater. Sci. - 2003. - V. 28. - P. 518528.
22. Kovalevskaya Z.G., Khimich M.A., Belyakov A.V., Shulepov I.A. Evaluation of physical and mechanical properties of structural components of Ti-Nb alloy // Adv. Mater. Res. - 2014. - V. 1040. - C. 3942.
Поступила в редакцию 04.07.2017 г.
Сведения об авторах
Романова Варвара Александровна, д.ф.-м.н., внс ИФПМ СО РАН, [email protected] Балохонов Руслан Ревович, д.ф.-м.н., внс ИФПМ СО РАН, [email protected] Батухтина Екатерина Евгеньевна, аспирант ИФПМ СО РАН, [email protected] Емельянова Евгения Сергеевна, магистрант ТГУ, [email protected] Сергеев Максим Владимирович, магистрант ТГУ, [email protected]
