Алгоритмы и программное обеспечение для моделирования импульсных процессов

Размещено на http://www.allbest.ru//

Размещено на http://www.allbest.ru//

Алгоритмы и программное обеспечение для моделирования импульсных процессов

В. Н. Аптуков

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

1. Введение© В. Н. Аптуков, А. В. Фонарев, Л. В. Ландик,

Д. В. Щеголев, 2010

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

Для моделирования этих процессов построена математическая модель и алгоритм решения задач динамического деформирования. Задача может рассматриваться в различных постановках - осесимметричная, плоское напряженное состояние, плоское деформированное состояние.

Алгоритм предусматривает использование различных конечно-элементных сеток (регулярных, нерегулярных, смешанных, сформированных в других пакетах). Элементы могут быть 3-угольные, 4-угольные. Конечно-разностная схема реализована в виде алгоритма "поэлементной сборки", что повышает ее эффективность, а также позволяет применять распараллеливание.

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

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

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

2. Уравнения движения

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

Рис. 1. Схема для осесимметричной задачи

Запишем уравнения сохранения импульса [1-3]:

(1)

Здесь - плотность среды; - компоненты вектора массовой скорости по осям и ; - компоненты тензора напряжений Коши. Система уравнений (1) описывает движение материальной частицы среды в эйлеровой системе координат под действием тензора напряжений .

Компоненты тензора скоростей деформаций связаны со скоростями частиц дифференциальными зависимостями Коши

(2)

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

Рис. 2. Область интегрирования для -го узла

Для получения разностных уравнений интегрируем уравнения (1) по площади (рис. 2).

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

В результате для компонент скорости получены соотношения

(3)

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

Суммирование производится по всем элементам, окружающим узел (рис. 2, 3).

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

Рис. 3. Вклад -го элемента в область интегрирования для -го узла

Для оптимального вычисления правой части в (3) при интегрировании (1) по области используется метод поэлементной сборки по всем элементам . Для этого вводятся понятия узловых усилий и узловых масс:

(4)

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

На рис. 3 для узла 1 показана площадь , а на рис. 4 - соответствующие площади для всех узлов данного элемента : в треугольнике - это , в четырехугольнике -.

Рис. 4. Схема поэлементного вклада при вычислении усилий и масс

Соответствующий вклад в узловые силы и массы, например для локального узла 1 треугольного элемента , имеет вид

В результате, с учетом центрирования по двум временным слоям, обеспечивающим 2-й порядок точности, получены разностные соотношения

- для компонент скорости:

(6)

- для координат узлов:

(7)

Шаг по времени определяется из условия устойчивости [3]

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

3. Физические уравнения

Рассмотрим для краткости только модель упруговязкопластичного материала (металлы).

Пластическое течение и объемная сжимаемость рассчитываются независимо [3]. Объемная сжимаемость, согласно принятому подходу, описывается зависимостями среднего давления от текущей плотности. Взаимовлияние может осуществляться через зависимость условия текучести от давления.

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

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

(8)

Отсюда скорости деформаций

. (9)

Тогда новые значения деформаций

, (10)

где- площадь элемента в момент , - приращения деформаций, - шаг по времени.

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

В элементе находим 2 треугольника по принципу короткой диагонали, определяем в них скорости деформаций, а затем усредняем по площади. Расчеты подтвердили правильность этого подхода.

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

.

Тензор напряжений представим в виде суммы шаровой и девиаторной составляющих

, (11)

где - единичный тензор, - девиатор тензора напряжений, - среднее давление.

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

, (12)

где - интенсивность напряжений.

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

Пластическое течение материала описывается моделью упруговязкопластической среды [1, 2]. Условие текучести записывается в форме Мизеса с пределом текучести , зависящим от деформации, скоростей деформаций и температуры. Согласно [1] запишем эту зависимость в виде

, (13)

где , , - функции, учитывающие деформационное упрочнение, вязкое упрочнение, температурное разупрочнение.

Деформационное упрочнение представляется зависимостью [1]

, (14)

где - параметры модели, подбираемые при описании экспериментальных кривых для статических испытаний при нормальной температуре; - статический предел текучести.

Для вязкого упрочнения применим аппроксимацию [1]

, (15)

где - экспериментальные параметры материала, .

Данные многочисленных испытаний показывают существенное изменение деформационных и прочностных свойств металлов при повышенных температурах. Зависимость предела текучести от температуры аппроксимировали функцией [1]

, (16)

где - абсолютная температура; - температура плавления; - параметры модели.

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

Температура определяется из уравнения баланса энергии за счет пластического деформирования. Влияние объемной сжимаемости для средних скоростей удара (до 1 км/с) не превышает 50-60С. Разностный аналог уравнения баланса энергии имеет вид

, (17)

где - теплоемкость.

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

В результате получим

(18)

где - изменение объемной деформации; - модуль сдвига материала.

В ходе вычислений проверяется условие текучести (12), которое отражает ограниченность девиаторных компонент тензора напряжений . Если в пространстве напряжений мы выходим за пределы упругой области, то вычисленные по закону упругости компоненты девиатора тензора напряжений "переносятся" по нормали на поверхность текучести [3].

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

, (19)

где - модуль объемного сжатия; - модуль и показатель адиабаты для данного материала.

Второе соотношение в (19) справедливо при больших давлениях, достигаемых при распространении ударных волн.

Для сглаживания разрывов в схемах сквозного счета введена искусственная вязкость [3], вычисляемая только при сжатии материала:

(20)

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

В результате получаем полные напряжения

. (21)

4. Граничные условия

Граничные условия в задачах соударения, проникания и действия взрыва на грунтовые массивы и элементы конструкций делятся на кинематические, динамические и смешанные [1-3].

Кинематические условия (на части поверхности заданы перемещения или скорости) реализуются на поверхностях жесткой заделки, осях и/или плоскостях симметрии.

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

Смешанные условия (на части поверхности заданы одновременно и перемещения (скорости) и поверхностные нагрузки) реализуются на контактных границах (деформируемое тело - продукты детонации, деформируемое тело - жесткое тело, деформируемое тело - деформируемое тело).

Температурные эффекты на контактных поверхностях не учитываются.

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

5. Программный комплекс для моделирования импульсных процессов

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

Препроцессор в интерактивном режиме выполняет выбор задачи, математической модели, параметров (физические, экспериментальные, схемные, конечно-элементную сетку), формирует для исполнительной части входной текстовый файл *.tm1. На рис. 5, 6 приведены две начальные формы препроцессора (выбор задачи, параметров).

Рис. 5. Начальная форма препроцессора

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

Конечно-элементные сетки могут быть построены препроцессором (регулярные сетки) или взяты из текстовых сеточных файлов, полученных в ANSYS [4] (2 файла *.txt) или MKE_MТДТ/ Mesh_32 [5] (файл *.grd).

В сеточных файлах предусмотрены все форматы матрицы конечных элементов: старый формат - [i, j, k)]; стандартные форматы ANSYS - [i, j, k, k]; [i, j, k, 0]; [i, j, k, l].

Рис. 6. Форма выбора параметров

Точность решения зависит от качества деформируемой сетки, поэтому в исполнительных модулях предусмотрен оперативный режим перестройки треугольной сетки [1]. Для каждого элемента определяется мера искажения эталонного элемента. Здесь - константы влияния формы и площади. Правильный выбор предельного значения меры искажения , , а также коэффициентов искусственной вязкости для треугольных (CL3,CV3) и четырехугольных (CL4,CV4) элементов позволяет получить решение, достаточно хорошо согласующееся с экспериментом.

Процесс решения сопровождается периодической визуализацией сетки деформируемого образца, записью промежуточных и окончательных результатов (НДС) в текстовые файлы (*.res,*.msh) для последующей графической визуализации постпроцессором. На каждом временном шаге в текстовый файл *.txt записываются номер шага N; текущее время T (мкс); укорочение L/L0 (см); текущая скорость V (м/сек); усилие в зоне контакта для последующего построения графиков и анализа решения задачи.

6. Примеры работы пакета программ

В качестве одной из задач импульсного деформирования рассматривается задача соударения упруговязкопластического цилиндра с жесткой стенкой. Эта задача является классическим тестом двух- и трехмерных схем для решения задач импульсного деформирования тел. Численное моделирование задачи соударения позволяет определять параметры моделей на основе сравнения с экспериментальными данными [6].

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

Сравнивалась величина укорочения цилиндров после удара. Получено совпадение результатов в пределах 7% для всех пяти металлов во всем диапазоне скоростей.

Рис. 7. Начальная и конечная форма цилиндра при ударе с различными скоростями (169м/с, …)

На рис. 8 сравниваются расчетные и экспериментальные [6] величины остаточного укорочения уранового цилиндра. Параметры материала: удельный вес = 18,7 гс/см3; объемный модуль K = 210 ГПа; модуль сдвига G = 80 ГПа; модуль адиабаты A=20 ГПа; показатель адиабаты = 5,8; теплоемкость = 0,96Кдж/кг·К; температура плавления = 14000С; статический предел текучести = 0,9 ГПа.

Рис. 8. Укорочение урановых цилиндров при ударе с различными скоростями

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

Результат расчетов удара на разных сетках

V0 (м/с)	L/Lo (эксперимент)	Расчет L/Lo на различных типах сеток
		^	¦	^+¦
158	0,970	0,966	0,961	0,965
181	0,963	0,958	0,953	0,956
203	0,958	0,953	0,947	0,951
205	0,953	0,947	0,940	0,945
273	0,932	0,925	0,917	0,922

На рис. 9 представлены результаты расчета соударения полой сферы и полого цилиндра с жесткой стенкой. Сетки инкорпорированы из пакета ANSYS [4]. Показаны начальная (слева) и конечная (справа) конфигурации. Рисунок иллюстрирует возможности расчета динамических задач на сетках произвольной структуры, что позволяет рассматривать тела произвольной формы.

Рис. 9. Расчет соударения сферы и цилиндра на сетках, сформированных в ANSYS

Пример расчета динамической осадки неоднородного образца в плоской постановке приведен на рис. 10.

Рис. 10. Расчет динамической осадки образца

Заключение

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

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

Проведены расчеты, которые показали, что каждый тип сеток оптимален для своего класса задач и условий деформирования.

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

динамический деформирование нагружение программный

Аптуков В.Н., Мурзакаев Р.Т., Фонарев А.В. Прикладная теория проникания. М.: Наука, 1992. 104 c.

Кукуджанов В.Н., Кондауров В.Н. Численное решение неодномерных задач динамики твердого тела // Проблемы динамики упругопластических сред. М.: Мир, 1975. C.39-84.

Уилкинс М.Л. Расчет упругопластических течений // Вычислительные методы в гидродинамике. М.: Мир. 1967. C.212-263.

ANSYS Basic Analysis Procedures Guide. ANSYS Release 11.0 / ANSYS Inc.

Аптуков В.Н., Ландик Л.В., Фонарев А.В. Метод конечных элементов и нерегулярные сетки для решения стационарных задач переноса тепла и статики упругих тел: учеб. пособие / Перм. ун-т. Пермь, 2002. 120 с.

Уилкинс М.Л., Гуинан М.У. Удар цилиндра по жесткой преграде // Механика: сб. пер. М.: Мир, 1973. №3. C.112-128.

Размещено на Allbest.ru

...