0

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ

УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ ИНДУСТРИАЛЬНЫЙ УНИВЕРСИТЕТ»

(ФГБОУ ВПО «МГИУ»)

Кафедра «Информационные системы и технологии»

ДИПЛОМНАЯ РАБОТА

по направлению «Математическое обеспечение и администрирование информационных систем»

студента Ушакова Ивана Александровича на тему «Разработка математической модели импульсного уплотнения песчано-глинистых литейных форм»

Руководитель работы: проф., д. т. н., Матвеенко Иван Владимирович Студент-дипломник _ И. А. Ушаков Руководитель работы, проф., д. т. н. _ И.В. Матвеенко

ДОПУСКАЕТСЯ К ЗАЩИТЕ

Зав. кафедрой к.ф.-м.н доцент _ Е.А. Роганов МОСКВА Оглавление Аннотация.

Введение

Литературный обзор.

Глава 1. Математическая модель процессов уплотнения.

Глава 2. Реализация математической модели средствами МКЭ.

Глава 3. Компьютерное моделирование ударно-прессово фильтрационного процесса уплотнения.

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

Заключение.

Литература.

Приложения.

Аннотация.

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

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

Работа содержит 55 страниц, 13 иллюстраций, 3 приложения.

Введение.

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

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

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

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

Результатом работы будет программа на языке Cи, которая производит расчёт на ЭВМ поля плотности формовочной смеси с нелинейно изменяющимися свойствами при импульсно-ударном методе уплотнения в сечении опоки с моделью с учётом больших перемещений смеси.

Литературный обзор.

Разработка технологии производства отливок – это одна из важнейших задач в литейном производстве. Самая трудоемкая операция в данном процессе – это операция формовки, получения качественной литейной формы. Одним из способов формовки является ударно-прессово фильтрационное уплотнение песчано-глинистых форм [8]. Конструктивное устройство установки такого типа представлено на рисунке Основные узлы установки: ресивер 1, внутри которого размещён пневмоударник, представляющий собой пневмоцилиндр 2, поршень 3, шток 4, на конце которого закреплена ударная плита 5. Пневмоцилиндр крепится на основании 6, в котором по центру отверстие для штока, ряд отверстий 7 для выпуска воздуха из ресивера в кожух 8, через отверстие 9 в ударной плите, а также каналы 10 соединяющие ресивер с полостью рабочего цилиндра. Ударная плита 5 в нижней части штока крепится гайками 11, а верхняя соединена замковым узлом 18, служащим для фиксации штока в верхнем положении и пуска ударного механизма уплотнения смеси в опоке. Вместе с формовочной установкой на рисунке 1 представлены также модельная плита с вентами 12, модели 13, опоки 14, наполнительная рамка 15 и формовочная смесь 16, прокладка 17.

Работа установки заключается в следующем:

1. Нерабочее положение: Поршень с пресс - плитой находятся в опущенном положении, давления воздуха в ресивере нет. Кран «б», впуска воздуха из магистрали, закрыт.

2. Подготовке к пуску (Подъём плиты): Подъём поршня с ударной плитой осуществляется путём подачи воздуха при открытых кранах «а» и «б». Кран «а» в положении впуска воздуха под поршень. После подъёма поршня с плитой производится их фиксация замковым узлом 18.

';

3. Пуск: Впуск воздуха в ресивер и цилиндр. Кран «б» и «а» открыты. После впуска воздуха в ресивер кран «б» закрывается. Кран «а» остаётся открытым. После чего защёлка фиксатора сбрасывается и происходит удар пресс - плиты по смеси. Воздух из поршневой полости вместе с воздухом из ресивера уходит через смесь в атмосферу.

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

Способ уплотнения литейных форм и стержней включает установку опоки на подмодельную плиту с моделью, засыпку смеси, размещение на опоке камеры воздушно – импульсной головки, уплотнение смеси по средствам фильтрации воздушного потока через формовочную смесь и импульсного удара прессовой колодки. Особенностями данного типа уплотнения являются: отработанный воздух из штоковой полости переходит под давлением в кожух и проходит через смесь в венты, тем самым совершая работу по уплотнению формовочной смеси; наличие вент вокруг модели и по периметру опоки.

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

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

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

Принято известное понятие допредельного и запредельного этапов деформирования. Под понятием «предельное напряжение смеси» понимается такое напряженное состояние, при котором отдельные песчинки сомкнулись друг с другом (без проскальзывания друг относительно друга) до «беспустотного» соприкосновения (рис. 2 б). Предельное состояние наступает при структурной плотности формовочной смеси. За этим пределом осуществляется переход смеси от упруго-вязкого состояния к упругопластическому. [11] В нашей работе учитывалось явление структурной плотности формовочной смеси, состоящее в том, что для уплотнения влажной смеси от насыпной плотности (при просеивании её через сетку) до структурной, с хорошей точностью равной насыпной плотности смеси в сухом состоянии, необходимы весьма малые давления. Это уплотнение происходит при первой фазе ударно-прессово фильтрационного процесса (воздухом). Для дальнейшего уплотнения образца требуется уже значительная ударно-прессовая нагрузка во второй фазе процесса (деформация за счет сжатия).

Рассмотрим второй этап ударно-прессово фильтрационного уплотнения. Он представляет собой удар прессовой колодкой, что схоже с высокоскоростным способом уплотнения. Исследованию и разработке высокоскоростного процесса посвящены работы И. В. Матвеенко [4], М. А. Шалимовой [3] и других авторов, в которых отмечается, что данный способ позволяет получать качественные литейные формы из смесей как высокой, так и низкой прочности с весьма малыми энергозатратами.

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

где k – отношение массы уплотняемой смеси к массе подвижных частей машины, R – коэффициент восстановления при ударе.

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

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

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

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

Исследуя высокоскоростное прессование песчано-глинистых смесей, отмечается появление объёмной волны сжатия, под действием которой частицы смеси начинают колебаться в направлении действия приложенной нагрузки. Увеличение скорости падения ударной колодки влечёт рост среднего напряжения от до МПа, а также остаточной деформации, которая при 3% влажности смеси составляет 0,35 – 0,4, а при влажности 5% 0,4 – 0,52.

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

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

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

с динамической стороны – возникновение, а затем исчезновение больших ударных сил.

Постановка любой задачи об ударе состоит в формулировке инерционных и реологических свойств тел, образующих исследуемую модель (выделение инерционных и безынерционных элементов; использование представлений об абсолютно твёрдых телах и о телах, обладающих упругостью, вязкостью и пластичностью), а также свойств, которые приписываются силам ударной природы (их распределённость или сосредоточенность в пространстве или во времени). [3] Одной из важнейших задач при разработке технологии производства отливок является изготовление качественной формы. Условием получения качественной литейной формы является тесная взаимосвязь трёх параметров технологического процесса: режим машины (т.е. давление прессования, давление сжатого воздуха в ресивере или масса и скорость падения ударной колодки), габариты оснастки (размеры опоки и модели) и реологические параметры смеси (вязкость, упругость и пластичность).

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

Задача определения режима формовки может быть решена с помощью натурных испытаний и при моделировании процессов уплотнения на ЭВМ. Экспериментальное изучение уплотнения смесей связано с большими затратами труда, энергии и материалов. Такого рода вопросы успешно решает проведение вычислительных экспериментов, основанных на создании математической модели процессов, происходящих в смесях при их уплотнении, поведение которых исследуется на ЭВМ гораздо быстрее и дешевле соответствующих лабораторных опытов. Кроме того, компьютерное моделирование позволяет проследить перемещение смеси в ходе её уплотнения и выбрать оптимальный режим формообразования и расположения моделей в оснастке.

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

Это возможно только при решении соответствующей математической модели на ЭВМ с применением современных численных методов.

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

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

В данной работе мы строим модель, описывающую поведение среды, отличающейся от идеализированных тел (формовочная смесь). Мы будем заниматься изучением и прогнозированием поведения формовочной смеси под воздействием внешних или массовых сил. Для решения этих задач, с одной стороны, изучим физические особенности формовочных смесей, а с другой – применим математический аппарат механики сплошной среды. Однако в отличие от этой дисциплины, допускающей, что тела являются сплошными, мы будем рассматривать формовочную смесь как пористую дисперсную среду, способную необратимо изменять свой объём, то есть уплотняться.

Конечной целью механики сплошной среды является математическое описание движения деформируемых тел. Под понятием «сплошная среда» подразумевают модель такого материального тела, которое хотя и состоит из отдельных частиц – атомов, молекул, но заполняет часть пространства непрерывным (сплошным) образом. Гипотеза сплошности (континуума) материального тела позволяет, рассматривая напряжение и деформацию бесконечно малых объёмов и используя аппарат дифференциального исчисления, переходить к изучению напряжённодеформированного состояния всего тела.

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

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

Для корректного решения задач об уплотнении необходим переход к моделям деформируемого тела. В моделях с распределёнными параметрами отражается действительная распределённость инерционных и деформационных свойств реальных тел. Конкретный вид уравнений зависит от формы соударяющихся тел, условий соударения и от принятого закона связи между напряжениями и деформациями. При малых деформациях рассматриваемого тела принимается, что материал следует закону Гука, но в более сложных системах, когда при ударе развиваются напряжения, превосходящие предел упругости материала, необходимо учитывать вязкие и пластические свойства материала. Весьма сложным оказывается явление удара в упругопластических стержнях, если производная увеличивается с возрастанием напряжения.

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

Математическое моделирование – это идеальное научное знаковое формальное моделирование, при котором описание объекта осуществляется на языке математики, а исследование модели проводится с использованием тех или иных математических методов.

При построении математической модели уплотнения формовочных смесей необходимо учитывать следующие специфические свойства, которые существенно усложняют задачу: большие деформации (перемещения слоёв), изменение реологических свойств формовочной смеси в процессе уплотнения, изменение модуля продольной деформации (модуля Юнга) и коэффициента поперечной деформации (коэффициента Пуассона) в зависимости от плотности смеси, переменная ударная нагрузка.

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

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

Существенным является также свойство полноты математической модели, когда модель даёт принципиальную возможность с помощью математических методов получить интересующие нас утверждения. К модели предъявляется требование устойчивости относительно погрешностей в исходных данных. Желательным является свойство наглядности математической модели, то есть ясный содержательный смысл её компонент, который даёт возможность проконтролировать модель, наметить план решения математической задачи, а также ориентировочно предвидеть результат решения, что может существенно ускорить процедуру.

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

Исследование одномерной задачи намного проще двумерной и ею занимались многие учёные. Сущность одномерного расчёта состоит в разбиении уплотняемой смеси на элементарные слои и последовательном определении для каждого из них напряжения сжатия и плотности с учётом граничных условий. Достоинство таких методик – в сравнительной простоте математической модели, основными их недостатками являются невозможность исследование опок с моделями и невозможность учёта силы трения на границах смеси с оснасткой. Данные недостатки существенны, так как все опоки на реальных предприятиях имеют модели (без них нет смысла в изготовлении литейной формы) и очень важным является распределение плотности в околомодельной зоне. Кроме того, в результате силы трения смеси на границе с оснасткой, слои смеси в ходе уплотнения двигаются неравномерно, и это играет большую роль в распределении плотности. Все эти недостатки можно исключить в результате рассмотрения двумерной задачи, трёхмерная задача может послужить лишь уточнением этих явлений, но никаких кардинально новых данных не даст, кроме того её моделирование значительно сложнее, поэтому мы остановимся на исследовании уплотнения в двумерной области.

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

В работе И. В. Матвеенко [4] разработан конечно-разностный метод нахождения двумерного напряжённо-деформированного состояния при воздействии динамической нагрузки, для чего весь уплотняемый объём разбивается на элементарные ячейки, для каждой из которых составляется система линейных уравнений. Задача решается пошагово, с корректировкой реологических свойств смеси и восстановлением координатной сетки на каждом шаге. Результатом решения являются напряжения сжатия, возникающие при уплотнении, скорости перемещения смеси на гранях ячеек и плотность смеси в каждой ячейке. В этом методе сложность заключается в аппроксимации граничных условий и построении нерегулярной конечно-разностной сетки. Это связано с точностью их аппроксимации, а также с видом самой границы: если она криволинейна, то задача представления граничных условий в конечно-разностной форме с сохранением аппроксимации становится проблемой. В этом заключается одно из существенных ограничений использования метода конечных разностей.

Одним из наиболее эффективных методов исследования сложного напряжённого состояния в механике сплошных сред является метод конечных элементов (МКЭ) вследствие простоты принципа, универсальности и удобства реализации на ЭВМ.

Важными характеристиками МКЭ являются хорошая аппроксимация криволинейной области, возможность применения неравномерной сетки разбиения области на элементы. Способ построения численного решения отличается от традиционных разностных схем в первую очередь принципом построения континуального приближённого решения. Так, в разностных схемах обязательно присутствуют этапы дискретизации, а затем уже проводится восполнение полученного дискретного решения до континуального, что порождает неоднозначность (особенно характерную для схем второго и выше порядков аппроксимации). В методе конечных элементов с самого начала построения численного решения ищется наилучшее (в той или иной норме) приближение точного решения в некотором пространстве (обычно это пространство кусочно-гладких функций). Таким образом, в этом подходе отсутствует этап восполнения. Можно считать подобные аппроксимации математически более строгими.

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

Методические принципы, используемые при разработке математических моделей формовочных смесей и теории уплотнения: сочетание феноменологической макроскопической теории, использующей гипотезу сплошности, и статистического подхода, основанного на введении средних для объёма частиц характеристик и изучении законов их распределения, а также доведение математического аппарата и расчётных зависимостей до конкретных инженерных приложений. Целесообразным представляется использование упрощённой математической постановки задачи с последующей разработкой физически обоснованных формул с экспериментально определёнными коэффициентами.

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

Для учёта изменения свойств формовочной смеси в процессе её уплотнения при исследовании сечения опоки с моделью необходимо привлечение результатов испытаний реологических характеристик в двумерных условиях. Такими характеристиками при использовании соотношений нелинейной упругости являются продольный модуль деформирования (модуль Юнга) и коэффициент поперечной деформации (коэффициент Пуассона).

В подавляющем большинстве реологических исследований авторы использовали в своих экспериментах образцы формовочной смеси с различной начальной плотностью. Это является причиной затруднений при попытке сопоставления полученных ими результатов. Между тем, в реальных формовочных процессах смесь начинает уплотняться, как правило, от плотности, полученной ею при свободной засыпке в опоку. Поэтому при численных расчётах процессов уплотнения литейных форм необходимо иметь диаграммы изменения реологических характеристик смесей во всём диапазоне плотностей.

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

Глава 1.

Математическая модель процессов уплотнения.

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

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

Уравнения движения:

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

В уравнениях (1), (2) третий член учитывает инерционность, а четвёртый – вязкость формовочной смеси.

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

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

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

где – величина напряжения, нормального к рассматриваемой границе, – скорость частиц смеси вдоль границы, – эмпирический коэффициент внешнего трения, подобранный из экспериментов работы [4]:

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

Рассмотрим первый этап этого процесса. Уплотнение потоком сжатого воздуха продолжается, пока ударная колодка двигается к верхнему слою смеси.

Величина ускорения для пневмоударного устройства определяется из уравнения движения:

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

где – эмпирический коэффициент трения в подвижных частях установки, для нашей установки (определен из [8]).

Время первой фазы ударно-прессово фильтрационного способа уплотнения рассчитывается следующим образом:

где – ускорение ударной колодки, рассчитанное по формуле (4), – расстояние, которое проходит ударная колодка до верхнего слоя смеси.

Скорость колодки в момент удара о верхний слой смеси:

При фильтрации через формовочную смесь потоком сжатого воздуха имеет место сила межфазного трения:

где – давление воздуха на входе в слой смеси, – давление воздуха на выходе из слоя смеси, – площадь опоки.

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

Сила динамической составляющей воздушного потока:

где – частота вынужденных колебаний смеси.

Таким образом, сила внешней нагрузки выражается следующим образом:

где – площадь опоки, – статический перепад давления, – плотность воздуха, – скорость воздушного потока.

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

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

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

Уравнение состояния для газа, находящегося в ресивере:

Если принять, что истечение воздуха из ресивера в надопочное пространство установки протекает по адиабатическому закону, а температура при этом постоянна, то уравнение теплового баланса при истечении воздуха из полости постоянного объёма (ресивера) имеет вид:

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

– давление в ресивере. Тогда расход воздуха из ресивера в технологическое пространство можно определить по известной зависимости [11]:

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

Формула для градиента давления:

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

Скорость воздушного потока вычисляется при помощи закона Бернулли:

Рассмотрим второй этап ударно-прессово фильтрационного уплотнения.

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

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

где m – масса подвижных частей механизма, P(t) – контактная сила удара.

При неупругом соударении пластическая деформация h пропорциональна силе удара где С – жёсткость кома формовочной смеси.

При решении дифференциального уравнения (12) получено выражение для контактной силы удара:

В нашем случае – скорость движения колодки в момент удара о смесь, которая вычисляется по формуле (6).

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

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

В нашем случае имеют место большие перемещения слоёв смеси при сжатии, обусловленные геометрической нелинейностью процесса уплотнения. В связи с этим в кинематических соотношениях, связывающих деформации с перемещениями U и V по координатным осям, необходимо учитывать не только линейные, но и квадратичные члены градиентов перемещений. Это увеличивает объём и сложность вычислений. Данные формулы были выведены Грином и Сен-Венаном (известно как тензор деформаций Грина) [1]:

где - компоненты тензора деформаций.

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

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

Уравнения состояния, связывающие компоненты тензора напряжений с компонентами тензора деформаций, должны учитывать нелинейное поведение реологических характеристик уплотняемой формовочной смеси и приняты в виде уравнений Генки [2]:

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

Соотношения Генки, устанавливающие связь между компонентами напряжений и деформаций (16), основаны на положениях теории пластичности, основными допущениями которой являются [2]:

деформация формы вызывается девиатором напряжений и не зависит от шарового тензора напряжений, деформация объёма вызывается шаровым тензором напряжений и не зависит от девиатора напряжений;

связь между компонентами напряжений и деформаций остаётся неизменной при любом виде напряжённого состояния;

напряжённое и деформированное состояния тела принимают подобными.

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

В классической теории пластичности принимается, что объёмные деформации упруги и полностью обратимы. Однако объёмные деформации формовочных смесей нелинейны и лишь частично обратимы (упругими деформациями можно пренебречь, так как они составляют небольшой процент от всей деформации).

Учитывая вышесказанное, уравнения (16) были переписаны следующим образом [2]:

Данные соотношения представляют собой закон Гука, используемый в классической теории пластичности механики сплошной среды, с той лишь разницей, что продольный модуль деформирования деформации ( ) являются переменными величинами. Данный факт является результатом физической нелинейности среды. Соотношения, связывающие модуль Юнга и коэффициент Пуассона с объёмной деформацией, зависят от свойств формовочной смеси и являются экспериментальными величинами.

равными нулю. Если при этом считать равным нулю, приходим к уравнениям плоско напряжённого состояния:

Если считать нулевой деформацию, то получим уравнения плоско деформированного состояния:

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

Проведённые эксперименты также показали, что модуль Юнга и коэффициент Пуассона в соотношениях (17) и (18) являются переменными величинами, нелинейно зависящими от достигнутой плотности формовочной смеси.

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

Метод конечных элементов (МКЭ) является сеточным методом, предназначенным для решения задач, в которых модель объекта задаётся системой дифференциальных уравнений в частных производных с заданными краевыми условиями. В основном МКЭ применяется для численного решения задач в механике сплошной среды.

Для таких задач необходимо сформулировать условие, позволяющее минимизировать невязку функционала по всей области. Это можно сделать при помощи метода взвешенным невязок. Однако данный метод обладает следующими недостатками: поскольку решение ищется сразу по всей области, то количество пробных и весовых функций должно быть значительным для обеспечения приемлемой точности, но при этом возникают трудности при вычислении, особенно при решении плоских и объёмных задач, когда требуется вычисление двойных и тройных интегралов по области с криволинейными границами. Поэтому на практике используют метод конечных элементов. Идея метода заключается в том, что непрерывную функцию, которую требуется вычислить (например, вектор перемещений), аппроксимируют дискретной моделью, которая строится на множестве кусочно-непрерывных (обычно линейных) функций. Таким образом, область сплошной среды разбивают на подобласти (конечные элементы), а точность решения задачи обеспечивают использованием большого количества конечных элементов. Обычно конечные элементы делают простой формы и вычисление интегралов по ним не вызывает особых затруднений. Конечные элементы строятся таким образом, что значения функции в любой точке элемента определяются через её значения в узлах (вершинах элемента). Для этого используют специальные базисные функции (функции формы), которые обладают следующими свойствами:

1. в узле аппроксимации функции имеют значение равное единице;

2. отличны от нуля в конечных элементах, содержащих этот узел аппроксимации;

3. во всей остальной области равны нулю.

Уравнения для элементов выводятся путём минимизации либо функционала, связанного с дифференциальным уравнением, либо энергии деформации на основе принципа стационарности (минимума) потенциальной энергии. Именно второй способ использован в данной работе, рассмотрим его подробнее.

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

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

Левая часть данного равенства представляет собой сумму вариации энергии деформации конструкции и вариации потенциальной энергии внешней нагрузки. – полная потенциальная энергия. Это означает, что для обеспечения равновесия полная потенциальная энергия должна принимать стационарное значение, что обеспечивает сходимость МКЭ. Можно показать, что для упругого материала полная потенциальная энергия не только стационарна, но и минимальна.

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

Глава 2.

Реализация математической модели средствами МКЭ.

Представленная математическая модель преобразована в матричное волновое уравнение движения [1]:

где * + – глобальный вектор перемещения узлов сетки конечных элементов, на которые разбито исследуемое сечение опоки с моделью;

[K] – глобальная матрица жёсткости, в построении которой участвуют кинематические соотношения и уравнения состояния;

[C] – глобальная матрица демпфирования, зависящая от вязкости среды;

[M] – глобальная матрица масс, зависящая от плотности среды;

{R} – глобальный вектор нагрузки, определяемый на первом этапе силой воздушного потока, на втором – ударной силой колодки и свойствами модельно-опочной оснастки.

Глобальные матрицы [K], [C] и [M] складываются из локальных матриц, относящихся к отдельным элементам:

элементов, то есть элементы локальных матриц записываются в те ячейки глобальной матрицы, к которым относятся соответствующие узлы конечного элемента (если узел используется несколькими конечными элементами, то в глобальной матрице производится суммирование значений локальных матриц соответствующих данному узлу). Пример составления ансамбля приведён на рисунке 3 [1]:

Используя преобразования, представленные в работе [1], можно показать, что формулы для локальных матриц будут выглядеть следующим образом:

где интегрирование ведётся по объёму V конечного элемента, значок «т» означает транспонирование матрицы, верхний индекс «е» здесь и далее определяет принадлежность соответствующего символа или понятия к элементу.

[N] – матрица формы конечного элемента;

– плотность смеси, заключённой в конечном элементе (постоянная для каждого шага по времени);

– вязкость смеси, отнесённая к единице площади в конечном элементе;

[B] – матрица градиентов конечного элемента, отвечающая за кинематические соотношения связи деформаций с перемещениями в элементе;

[D] – матрица реологических свойств смеси в конечном элементе, отвечающая за связь напряжений с деформациями в элементе.

Вид матриц [N], [B] и [D] приведён ниже.

Из формул (21) – (23) видно, что локальные матрицы жёсткости, демпфирования и масс определяются не только свойствами среды, заключённой в соответствующем конечном элементе, но и геометрией самого элемента.

Таким образом, в уравнении (19) глобальные матрицы [K], [C] и [M] зависят от свойств смеси и геометрии разбитого на конечные элементы сечения литейной формы, а вектор {R} отвечает за действующую в этом сечении нагрузку в зависимости от вида формовочного процесса (на первом этапе импульсно-ударного уплотнения он определяется статической и динамической силами воздушного потока, на втором – ударной силой колодки).

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

На рисунке 5 показан произвольный элемент единичной толщины, узлы которого пронумерованы буквами i, j, k против часовой стрелки начиная с узла i. Перемещение узла i определялось двухкомпонентным вектором-столбцом * + * +, так как задача двумерная (узел i имеет две степени свободы). - перемещения узла i по осям X и Y соответственно.

Рис. 4. Конечно-элементная дискретизация области треугольными элементами.

Аналогично перемещения узлов j и k представлялись вектор столбцами * + и * +. Перемещение элемента в целом задавалось 6-компонентным вектором Перемещение любой точки с координатами x и y внутри элемента определяется аналогично вектор столбцом, также имеющим две компоненты, но в виде функций, зависящих от координат интересуемой точки:

Для аппроксимации перемещения внутри элемента используют базисные функции.

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

где – базисные функции, а – площадь треугольного конечного элемента.

Пример треугольного конечного элемента и значений его базисных функций представлен на рисунке 6 [7]:

В результате для произвольного треугольного конечного элемента можно составить линейный полином, который будет однозначно определять перемещения внутри элемента вдоль оси X и вдоль оси Y:

и определяются аналогично.

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

С учётом определения функций формы перемещение в элементе найдено по зависимости [1]:

При использовании в работе треугольных конечных элементов постоянной толщины формулы вычисления локальных матриц масс (23) и демпфирования (22) преобразованы следующим образом [1]:

демпфирования треугольного элемента единичной толщины и площади [1]:

Большие перемещения формовочной смеси, достигаемые при уплотнении, учтены при помощи матрицы градиентов [B]:

Поскольку в соответствии с кинематическими соотношениями (уравнениями Грина) (14) все компоненты тензора деформаций * + представлены в виде суммы линейной * + и нелинейной * + частей, то и матрица градиентов [B] представлена такой же комбинацией матриц. Эти матрицы служат для выражения линейной и нелинейной частей деформации конечного элемента через его перемещения. Данное соотношение представлено в работе [1]:

формы, составим следующую матрицу для линейных членов в уравнениях для деформации:

, Значения матрицы справедливы, поскольку в формулах (22), (23) от x и y зависят дифференцирования Для удобства практических расчётов нелинейная часть деформаций представлена в виде [1]:

,где матрица [G] составляется аналогично линейной части матрицы градиентов, через функции формы и последующего их дифференцирования.

,Нелинейная часть матрицы градиентов зависит от перемещений узлов конечного элемента (так как в матрице [A] присутствуют перемещения) и соответственно принята в виде [1]:

Напряжения, действующие в элементе, связаны с его деформацией посредством матрицы [D], отражающей свойства формовочной смеси:

Матрица свойств построена при помощи соотношений Генки, которые устанавливают связь между компонентами тензора напряжений и тензора деформаций. Для плосконапряжённого состояния [D] имеет вид:

,Для плоскодеформированного состояния матрица [D] определяется по формуле:

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

Экспериментальное определение реологических характеристик при уплотнении смеси проведено в работе [3] для определённой глинисто-песчаной смеси. Мы будем использовать смесь с такими же составом и характеристиками (структурная плотность, внешний коэффициент трения, вязкость), и воспользуемся полученными соотношениями. Для практических расчётов следует применять диаграммы зависимости реологических коэффициентов смеси от плотности исследуемых объектов, а не от их относительной деформации. Зависимость продольного модуля деформации (модуля Юнга) от плотности является расчётной относительно функций модуля компрессионного сжатия ( ). Модуль компрессионного сжатия установлен в условиях одноосного сжатия смеси в объёмах с боковыми стенками. Вид функций ( ), отображающей затухающий характер развития объёмных зависимости деформаций с ростом всестороннего давления (рис. 7).

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

Данные величины вычисляются отдельно для каждого конечного элемента с учётом достигнутой в нём плотности на i-м временном шаге. В этом случае продольный модуль деформирования формовочной смеси определяется исходя из следующей зависимости:

Рис. 7. Вид зависимости напряжений от деформаций.

В свою очередь формула (21) для вычисления локальной матрицы жёсткости при использовании треугольных конечных элементов преобразована к виду:

где - площадь конечного элемента.

За граничные условия отвечает вектор нагрузок * +, в рассматриваемом случае содержащий количество компонент, равное удвоенному числу узлов треугольной сетки, при этом каждому узлу соответствует вертикальная (y) и горизонтальная (x) составляющие внешней силы на первом и втором этапах ударно-прессово фильтрационного уплотнения. Условия компрессионного сжатия, содержащие запрет на перемещения в определённом направлении узлов, расположенных на поверхностях, соприкасающихся с оснасткой, удовлетворяются обнулением соответствующих компонент вектора перемещений * + при помощи вектора * +.

Вектор * + учитывал также нагрузку от силы трения, возникающей на границе смеси с модельно-опочной оснасткой. Сила трения вычислялась по формуле (3) и относилась к узлам, находящимся на границе с оснасткой. При этом для узлов на вертикальных и горизонтальных границах заполнялись соответственно x и y составляющие вектора * + относящиеся к этим узлам.

Статическая и динамическая силы воздушного потока, уплотняющая ударная нагрузка колодки, а также нагрузка от давления верхнего поршня, действующие по оси y, учитывались компонентами вектора * +, которые соответствуют узлам контрлада формы. Начальная плотность смеси принимается равной структурной (плотность засыпанной в опоку неуплотнённой формовочной смеси).

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

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

Производные вектора перемещений были заменены их разностными аналогами на временном шаге:

где i – номер итерации, а – величина шага по времени.

В соответствии с этим образована система линейных алгебраических уравнений относительно вектора перемещений на каждой итерации:

Разрешим данное матричное уравнение относительно неизвестного вектора перемещений на i+1 шаге:

В квадратных скобках указаны результирующие матрицы, которые получаются из исходных путём линейных преобразований по правилам сложения матриц и умножения матрицы на коэффициент. Таким образом, зная перемещения на i и i- шаге, мы можем получить вектор перемещений узлов конечных элементов на i+ шаге, решив данное матричное уравнение. В работе для этого использован метод Гаусса.

Так как первая и вторая производные от перемещений по времени вычисляются в средней точке разностной схемы, то в этой же точке вычисляются и все значения, ходящие в уравнение (36). В частности, матрицы градиентов [B] и свойств материала [D], которые участвуют в построении глобальной матрицы жёсткости, получаются зависящими от достигнутых значений перемещений и, таким образом, изменяются по времени. Такой же динамический характер имеет вектор нагрузок системы {R}.

Очень важным и сложным вопросом при решении матричного волнового уравнения движения итерационным способом является вопрос выбора размера временного шага, так как от этого в значительной степени зависит точность полученного решения. В рассмотренной литературе нет однозначного мнения на этот счёт. Было принято решение использовать метод, предложенный в работе [3]. Где при вычислениях величина временного шага подбиралась с использованием условия Куранта, согласно которому за промежуток времени распространяющаяся волна должна пройти расстояние, не большее половины минимального размера элемента.

Таким образом где k – коэффициент, учитывающий скорость роста внешней нагрузки, L – наименьшая сторона наименьшего элемента в сетке, a – максимальная скорость звука в среде. В твёрдой среде скорость звука для продольной волны равна где E – продольный модуль деформирования (модуль Юнга), – коэффициент поперечной деформации (коэффициент Пуассона), – плотность среды.

Коэффициент k находится в пределах 0 – 0,5 и выбирается в зависимости от времени возрастания функции нагружения. В рассматриваемом случае быстрого роста внешней нагрузки выбиралось значение k из интервала 0 – 0,3.

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

При этом,а определяется как градиент давления и находится по формуле (10) с использованием конечно-разностной аппроксимации:

где – объём надопочного пространства, – рассчитывается по формуле (9).

Формула скорости уменьшения давления в ресивере:

где – объём ресивера.

Также будем определять скорость воздушного потока на каждом временном шаге при помощи закона Бернулли (11):

Тогда:

где – расстояние от верхнего слоя смеси до основания пневмоцилиндра, – плотность воздуха, – давление в ресивере на i-м шаге по времени:

– давление над смесью на i-м шаге по времени:

Решением уравнения (36) является вектор перемещений узлов конечноэлементного представления рассматриваемого сечения литейной формы. Затем по нему находятся деформации и плотность в элементах среды.

Плотность смеси в элементе на i-м временном шаге вычисляется по формуле:

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

Компоненты тензора деформаций вычислялись при помощи текущего вектора перемещений и матрицы градиентов [B].

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

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

Таким образом, разработанная на основе МКЭ методика расчёта плотности формовочной смеси позволяет учесть нелинейный рост её реологических характеристик и уплотняющей нагрузки при импульсно-ударном процессе с большими пластическими деформациями в опоках с моделью.

Глава 3.

Компьютерное моделирование ударно-прессово фильтрационного процесса уплотнения.

На основе предложенной методики построен алгоритм численного расчёта ударно-прессово фильтрационного уплотнения формовочной смеси.

1. Задание начальных данных: разбиение исследуемого сечения литейной формы на конечные элементы и задание геометрии дискретной области; задание начальной плотности и законов изменения реологических характеристик формовочной смеси в зависимости от достигнутой плотности; задания параметров импульсноударного формовочного процесса (давление в ресивере, диаметр впускного клапана, объём ресивера, массу подвижных частей ударного механизма).

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

2.1. Вычисление локальных матриц жёсткости по формуле (35), локальных матриц демпфирования по формуле (28) и локальных матриц масс по формуле (27) с учётом достигнутой величины плотности в элементе.

2.2. Построение глобальных матриц жесткости, масс и демпфирования из локальных матриц при помощи формул (20) по правилу составления ансамбля в методе конечных элементов.

2.3. Вычисление соответствующего временному шагу вектора нагрузки системы.

Для граничных узлов вычисляется сила внешнего трения по формуле (3), для узлов контрлада формы на первом этапе уплотнения вычисляется сила воздушного потока по формуле (7), на втором – вычисляется сила удара по формуле (13). Время окончания первого этапа вычисляется по формуле (5).

2.4. Учёт ограничений на перемещения смеси на границах с оснасткой при помощи вычёркивания из матричного уравнения (36) строк и столбцов, отвечающих за перемещение в горизонтальном направлении для узлов, прилегающих к вертикальным стенкам опоки и модели, и отвечающих за вертикальные перемещения для узлов, прилегающих к горизонтальным поверхностям оснастки.

2.5. Решение системы линейных алгебраических уравнений (36).

2.6. По найденным значениям перемещений вычисление в каждом элементе деформаций по формуле (29).

2.7. Вычисление напряжений при помощи полученных деформаций по формуле 2.8. Вычисление плотности в соответствии с формулой (42).

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

2.10. Переход к следующему временному шагу.

Блок-схема вычислительной программы представлена на рисунке 8.

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

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

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

Разработанная прикладная программа требует для своей работы дополнительные входные параметры, а именно:

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

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

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

задание нагрузок – это выражается в задании особенностей импульсноударной формовочной установки (давление в ресивере, объём ресивера, диаметр впускного клапана, площадь вент, масса подвижных частей механизма колодки, жёсткость, кроме того необходимо указать коэффициент внешнего трения смеси с оснасткой).

Рис. 8. Блок-схема программы.

В работе использованы координаты X и Y узлов треугольной сетки, расположенные в определённом порядке.

Для обеспечения максимальной простоты и гибкости программы нагрузки вычислялись в виде одного вектора, который также отвечал за силу трения смеси на границах с опокой и моделью, построенный вектор затем непосредственно использован в программе. Поскольку рассматривались нестационарные граничные силы, для формирования вектора нагрузки в зависимости от геометрии системы и времени был использован специальный метод для вычисления значений узловых сил при действии распределённой нагрузки. При составлении вектора внешних нагрузок использовались массивы из входных данных, после чего компоненты вектора, соответствовавшие левым, правым и нижним границам, заполнялись значениями вычисленной силы трения. Компоненты вектора внешних нагрузок, соответствующие узлам контрлада формы, в интервале времени, соответствующем уплотнению фильтрацией потока сжатого воздуха, заполнялись вычисленным в методе значением силы воздушного потока, а затем, в последующем интервале времени уплотнения колодкой, заполнялись значением силы удара колодки. Во время уплотнения потоком сжатого воздуха на каждом шаге по времени вычислялся статический перепад давления по формуле (38), который использовался для нахождения силы статического перепада давления воздушного потока, также вычислялась скорость воздушного потока по формуле (41), которая использовалась для нахождения силы динамической составляющей воздушного потока. После чего высчитывалась сила внешней нагрузки для этапа уплотнения воздухом (7). Ударная сила колодки второго этапа уплотнения определялась по формуле (13). У вектора нагрузки заполнялись только те компоненты, которые отвечали за движение узла в определённом направлении (для узлов на левой, правой и верхней границах по оси Y, для узлов на нижней границе по оси X).

Граничные условия были учтены путём соответствующего изменения системы алгебраических уравнений непосредственно перед её решением. В системе (36) из матрицы в левой части вычёркивались строки и столбцы, отвечающие за перемещение прилегающих к границе узлов в перпендикулярном направлении к поверхностям опоки или модели. Соответственно в искомом векторе и векторе правой части матричного уравнения также вычеркивались данные элементы. Далее полученная система решалась прямым методом исключения Гаусса, который позволяет получить точное (в пределах ошибки округления) решение, а результирующий вектор дополнялся нулевыми значениями для вычеркнутых компонент.

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

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

Последним этапом решения задачи явилось вычисление плотностей для элементов по найденным деформациям при помощи формулы (42) (для этого в программе отводится отдельный одномерный массив).

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

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

Результат данной дискретизации представлен на рисунке 9.

Характеристики расчётной области: 50 узлов, 69 элементов, 10 слоёв.

Шаг по времени определяется при помощи условия Куранта (29), установлено, что условие устойчивости разностной схемы нарушается, если dt превосходит предельное значение порядка 1 мс (для каждой задачи оно определяется экспериментально).

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

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

Рис. 9. Конечно-элементное разбиение сечения опоки с моделью.

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

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

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

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

При ударно-прессово фильтрационном уплотнении принималось, что смесь до впуска сжатого воздуха и начала движения колодки покоится (то есть скорость всех её узлов равна нулю), её плотность и реологические характеристики во всех точках формы одинаковы. Во время продувки воздухом происходит предварительное уплотнение смеси, в момент падения на смесь ударной колодки на верхние слои контрлада формы начинает действовать ударная нагрузка, развивающаяся за время удара от нуля до максимального значения; при этом верхние слои смеси приходят в движение. Следует отметить, что для учёта одинаковых перемещений узлов (вертикальные перемещения на колодке) был применён подход, состоящий в усреднении соответствующих компонент вектора перемещений на текущем временном шаге уплотнения ударной колодкой. Другие подходы, формализующие принцип одновременного перемещения узлов, не были найдены в рассмотренной литературе по этому вопросу.

В ходе моделирования было выявлено, что перемещение узлов в горизонтальном направлении крайне незначительно. Процесс уплотнения формовочной смеси носит больше послойный характер, поскольку движение узлов происходит преимущественно в вертикальном направлении, а перемещение узлов, лежащих на одной высоте, практически одинаково и различается только вблизи стенок, где действует сила трения. Тем не менее, некоторые конечные элементы испытывают деформации сдвига в направлении оси X и даже деформации скручивания (разворачиваясь на определённый угол), но данные виды деформации незначительны по своей величине.

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

Глава 4.

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

Способ ударно-прессово фильтрационного уплотнения песчано-глинистых форм был зарегистрирован в Государственном реестре изобретений Российской Федерации 10 апреля 2010 года Иваном Владимировичем Матвеенко и Сергеем Александровичем Кондратьевым, о чем свидетельствует патент в приложении 3.

Для проверки адекватности разработанного компьютерного метода расчёта поля плотностей в сечении литейной формы, уплотняемой ударно-прессово фильтрационным методом, сравнивались результаты работы программы и экспериментальных данных. К сожалению, мы не располагали установкой по данному процессу, из-за конструктивных особенностей и массы более одной тонны её нет в лаборатории кафедры литейного производства МГИУ. Поэтому для проведения эксперимента нам пришлось сымитировать работу установки ударнопрессово фильтрационного уплотнения на имеющихся машинах.

Данный процесс можно разделить на две части:

1. Предварительное уплотнение фильтрацией потоком сжатого воздуха (кратковременная продувка) 2. Последующее уплотнение плоской прессовой колодкой.

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

Сначала мы вычислили время первой фазы и необходимое давление в ресивере, и произвели уплотнение потоком сжатого воздуха на лабораторной импульсной установке, расположенной в лаборатории кафедры «Литейное производство» в МГИУ.

Принцип работы данной машины (рис. 10) заключается в том, что при открытии клапана 2 поток сжатого воздуха из ресивера 1 через отверстие 3 мгновенно воздействует на формовочную смесь 4, находящуюся в опоке 6, разгоняет её в направлении модели 5. Отработанный сжатый воздух через специальный клапан сброса давления или через венты в модельной оснастке, а частично через неплотности соединений оснастки и машины, уходит в атмосферу. [9] Рис. 10. Схема лабораторной импульсной установки.

После чего перенесли опоку с предварительно уплотненной смесью в лабораторную ударно-прессовую установку, находящуюся в лаборатории кафедры «Литейное производство» в МГИУ, и произвели допрессовку.

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

Рис. 11. Схема лабораторной ударно-прессовой установки.

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

В настоящей работе применялась смесь следующей рецептуры:

Механические свойства:

Для увеличения точности измерения плотности формовочная смесь была просеяна через сито.

Смесь вначале засыпалась в опоку размерами 290X190X150 мм и насыпную рамку высотой 100 мм, на дне опоки располагалась модельная плита с моделью в виде прямоугольного параллелепипеда размерами 195X95X50 мм. Затем на стенде воздушно-импульсного уплотнения было совершено воздействие на формовочную смесь в опоке в течение 0,1 секунды с давлением в ресивере, равным трём атмосферам. Время действия первого этапа уплотнения определялось по формуле (5). Ускорение для пневмоударного устройства высчитывалось по формуле (4). Масса подвижных частей ударного механизма в установке равна 20 кг.

– диаметр рабочего цилиндра;

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

После этого опоку с предварительно уплотненную смесь мы поставили на ударнопрессовый стенд, где производилась допрессовка. Скорость ударной колодки в момент касания верхнего слоя смеси составила 5 м/c. Скорость движения колодки регулировалась давлением в пневмоударнике и была выбрана с учетом вычислений по формуле (6).

Для нахождения распределения плотности смеси по сечению формы использовали отбор образцов смеси шприцеобразными пробоотборниками с диаметром 20 и мм и высотой, соответственно, 25 и 50 мм из различных точек формы. По массе и объёму образца вычислялась их плотность. Масса образцов определялась на аналитических весах OHAUS Scout Pro ISO 9001 с точностью до 0,001 грамма.

В нашей работе насыпная плотность аэрированной влажной смеси составляла а насыпная плотность сухой смеси 1221 За начальную плотность формовочной смеси была принята структурная плотность.

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

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

В целом расхождение вычислительных и натурных экспериментов не превышает 10% и свидетельствует о количественной адекватности разработанной компьютерной модели реальным импульсно-ударным формовочным процессам. Сравнение расчётного и экспериментального распределений плотности формы по высоте опоки с моделью, уплотнённой ударно-прессово фильтрационным способом, представлено на рисунках 12 и 13.

Рис. 12. Расчётное и экспериментальное распределения плотности по высоте в кармане Рис. 13. Расчётное и экспериментальное распределения плотности по высоте в Заключение.

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

Способ ударно-прессово фильтрационного уплотнения песчано-глинистых форм был зарегистрирован в Государственном реестре изобретений Российской Федерации апреля 2010 года Иваном Владимировичем Матвеенко и Сергеем Александровичем Кондратьевым, о чем свидетельствует патент в приложении 3.

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

Конечным результатом работы стала программа, написанная на языке Cи, которая рассчитывает плотность смеси на ЭВМ.

Были проведены эксперименты на установках в лаборатории кафедры «Литейное производство» в МГИУ. Также были проведены расчёты по компьютерной визуализации картин динамического уплотнения литейных форм. Сравнение результатов вычислений с экспериментальными данными показали качественную и количественную адекватность разработанного метода компьютерного моделирования процессов уплотнения реальным условиям изготовления литейных форм.

Данная программа может быть использована на литейных заводах для расчёта плотности формы после уплотнения на соответствующей установке, а также в качестве обучающего материала для студентов кафедры «Литейное производство» в МГИУ.

Литература.

[1] Зенкевич О. Метод конечных элементов в технике. – М., Мир, 1975.

[2] Вялов С. С. Реологические основы механики грунтов. – М., Высшая школа, 1978.

[3] Шалимова М. А. Моделирование высокоскоростных процессов уплотнения литейных форм. Дисс. … канд. техн. наук. – М., 1998.

[4] Матвеенко И. В. Процессы динамического уплотнения литейных форм и выбор параметров формовочных машин. – М., МГИУ, 1979.

[5] Матвеенко И. В. Оборудование литейных цехов. – М., МГИУ, 2003.

[6] Седов Л. И. Механика сплошной среды. – М., Наука, 1970.

[7] Hunter P., Pullan A. Fem/bem notes. – New Zealand, The University of Auckland, 2002.

[8] Матвеенко И. В., Кондратьев С. А., Зубарев А. А. Способ ударно-прессово фильтрационного уплотнения песчано-глинистых форм и установка для его осуществления. – Патент РФ №2385784 от 10.04.2010.

[9] Матвеенко И. В., Исагулов А. З., Дайкер А. А. Динамические и импульсные процессы и машины для уплотнения литейных форм. – Алма-Ата, Наука, 1998.

[10] Иванов Ю. И., Благонравов Б.П., Волкомич А. А., Спицкиц В. П. Сборник научных трудов МАМИ. – М., 2002, с. 43.

[11] Каменский В. В. Теоретические основы разработки технологии получения высококачественных отливок со сложной ребристой поверхностью. Дисс.

доктор техн. наук. – М., 2005.

Приложение 1.

Результаты работы программы.

Средняя плотность конечных элементов, уплотнённых импульсно-ударным методом. Дискретизация сечения опоки с моделью на конечные элементы соответствует разбиению, представленному на рисунке 9. Нумерация элементов ведётся слева направо и снизу вверх.

0: 1495.8 1: 1495.8 2: 1495.8 3: 1495.8 4: 1495.8 5: 1518.8 6: 1518.8 7: 1518.8 8:

1518.8 9: 1518.8 10: 1541.8 11: 1541.8 12: 1541.8 13: 1541.8 14: 1541.8 15: 1564. 16: 1564.84 17: 1564.85 18: 1564.86 19: 1564.83 20: 1595.8 21: 1595.8 22: 1595.8 23:

1595.8 24: 1585.8 25: 1585.84 26: 1585.85 27: 1585.86 28: 1585.83 29: 1609.58 30:

1609.49 31: 1609.41 32: 1609.25 33: 1609.1 34: 1609.38 35: 1609.67 36: 1609.33 37:

1630.57 38: 1630.46 39: 1630.39 40: 1630.24 41: 1630.09 42: 1630.37 43: 1630.64 44:

1630.29 45: 1662.04 46: 1660.69 47: 1659.39 48: 1656.44 49: 1653.48 50: 1658.23 51:

1662.99 52: 1662.42 53: 1682.62 54: 1681.24 55: 1680.05 56: 1677.23 57: 1674.32 58:

1679.06 59: 1683.5 60: 1682.69 61: 1739.52 62: 1739.41 63: 1739.56 64: 1739.72 65:

1739.76 66: 1739.79 67: 1739.42 68: 1739. Приложение 2.

Код программы.

#include #include #include #include #define H 0. #define W 0. #define M 20. #define C #define S_DENSITY #define NUM_NODE #define NUM_ELEM #define B_SIZE #define V #define K 0. #define P_R #define P_N #define F_C 0. #define K_T 0. #define L 0. #define S_0 0. #define RO_V 1. #define MU_V 0. #define S_V 0. #define G 9. #define R 8. #define T #define V_N 0. #define V_R 0. #define M_V 0. #define HI 1. double delta_p;

double p_res;

double p_nadop;

void print(double* x, double* y){ int i;