Full Text

Введение Развитие современных высокоэнергетических методов обработки материалов, позволяющих получать изделия с заданными эксплуатационными свойствами, требует построения связанных термоструктурно-деформационных математических моделей. Достаточно широкое распространение численное моделирование получило при исследовании таких высокотемпературных технологических процессов, как сварка [1-5], литье [6-9], закалка [10-12], поверхностное упрочнение [13-15] и др. В работе [13], развитием которой можно считать данную статью, рассматриваются основные особенности построения систем математического моделирования интенсивных температурных и силовых воздействий, характерных для процессов поверхностного упрочнения металлов концентрированными потоками энергии (при электромеханической, лазерной, плазменной обработке и т.д.). При расчётном анализе напряженно-деформированных состояний (НДС) и механического поведения таких систем в каждый момент времени приходится иметь дело, по сути, с новым телом, структура и свойства которого непрерывно трансформируются в ходе технологического процесса под действием эволюционирующих полей температуры, структуры и напряжений. Это требует привлечения моделей смежных с механикой дисциплин, описывающих процессы тепло- и массопереноса, структурно-фазовых превращений, образования дефектов и др. В данной работе в рамках общей постановки контактной связанной задачи механики по определению НДС неоднородных тел при высокоэнергетических воздействиях выделены задачи теплопроводности и изменения фазового состава, связанного с мартенситным превращением в металлическом сплаве, характерные для исследуемого технологического процесса поверхностного упрочнения деталей. Точность решения таких задач во многом определяется обоснованностью и адекватностью учета наиболее значимых нелинейных эффектов при анализе тепловых процессов (зависимость теплофизических свойств от температуры, нелинейность граничных условий), а также связанностью эффектов, возникающих, например, в результате быстротекущей трансформации структуры (скрытая теплота фазовых превращений, движение межфазных границ). Учет таких особенностей может приводить как к количественному, так и к качественному уточнению результатов [16-19]. Кроме того, в ряде случаев накладываются определенные ограничения на методы и параметры решения нестационарного уравнения теплопроводности, например, при выборе схем интегрирования. В рамках данной статьи на основе метода конечных элементов (МКЭ) излагаются особенности решения термоструктурной задачи о поверхностном упрочнении электромеханической обработкой (ЭМО) титанового сплава Ti6Al2V. Проводится анализ влияния зависимости теплофизических свойств от температуры, радиационного теплообмена, скрытой теплоты фазовых превращений при различных способах их учета, а также величины выбранного временного шага на точность решения связанной задачи. При этом задача расчета активных и остаточных напряжений при контактном термосиловом нагружении в ходе поверхностной упрочняющей обработки металлических тел будет рассмотрена авторами дополнительно. 1. Постановка задачи 1.1. Основные уравнения и граничные условия Решение комплексной задачи основывается на системе дифференциальных уравнений, которые содержат в явном виде связанные члены. (1) (2) где - оператор набла в актуальной конфигурации деформируемого тела; σ - тензор напряжений Коши; f - вектор объемных сил; λ - тензор теплопроводности; T - температура; Qдж - джоулево тепло; Qфаз - скрытая теплота фазовых превращений; Qмех - теплота, соответствующая диссипации механической энергии при пластическом деформировании. Основные уравнения дополняются начальными и граничными условиями применительно к импульсной (переменным током) ЭМО детали твердосплавным инструментом в форме ролика (рис. 1). Рис. 1. Расчетная схема неоднородного тела с трансформирующейся структурой и граничные условия при контактном температурно-силовом воздействии (в ходе ЭМО) Fig. 1. Computational scheme for the heterogeneous body with a transformed structure and boundary conditions under thermo-force contact loading (during EMT) Контактный инструмент моделируется абсолютно жёстким эллипсоидным штампом. Силовые граничные условия задаются равенством нулю на свободной поверхности соответствующих нормальных и касательных напряжений, которые дополняются краевыми условиями Герца-Синьорини. На границе структурных областей записываются условия совместности перемещений и поперечных (к поверхности раздела с нормалью в данной точке) нормальных и касательных напряжений. На бесконечности устанавливается равенство перемещений и напряжений нулю. Принимается гипотеза о естественном, ненагруженном состоянии тела в начальный момент времени - напряжения и деформации отсутствуют. (3) где n - вектор нормали к поверхности деформируемого тела в текущей конфигурации; fвн - вектор внешней нагрузки, действующей на ролик; tc - интенсивность нагрузки в контактной зоне; gn - функция расстояния; - контактное давление. Расчетная область вследствие малой кривизны поверхности и локальности воздействия рассматривается в виде полупространства с подобластями (зонами) с различающимися физико-механическими характеристиками. Границы данных зон соответствуют областям металла с разной структурой, формирующейся в ходе термосилового нагружения (например, фрагментам со структурой закалки и исходного состояния; зоны 1, 2 на рис. 1). Их текущее положение устанавливается на каждом расчетном шаге из решения связанной задачи. В зоне контакта эллипсоидного штампа и детали задается эквивалентный равномерно-распределенный поверхностный источник тепла q. Температура в начальный момент времени и на бесконечности равна температуре окружающей среды , а также краевые условия на поверхности для уравнения (2) заданы системой. (4) где h - коэффициент конвективного теплообмена; σ - постоянная Стефана-Больцмана; ε - коэффициент теплового излучения. В случае ЭМО переменным током удельная мощность источника определяется по формуле , (5) где I - действующее значение силы тока; U - действующее значение контактного напряжения; k - коэффициент, учитывающий теплообмен инструмента и детали; ν - частота тока. В данной статье основной акцент делается на методах и особенностях решения уравнения (2) как составной части общей постановки. Решение термоструктурной составляющей рассматривается без учета диссипации механической энергии при пластическом деформировании. Отметим, что при численном решении задачи теплопроводности область действия теплового источника q принималась эллиптической в соответствии с размерами пятна контакта инструмента с деталью, установленными экспериментально. 1.2. Модель структурно-фазовых превращений Математическая модель трансформации микроструктуры при тепловом воздействии в ходе ЭМО титанового псевдо-α-сплава Ti6Al2V представляет собой моделирование фазового превращения при нагреве и (мартенситного) превращения при охлаждении на основе анализа термокинетической диаграммы распада переохлажденной β-фазы [20-23]. В текущий момент времени в рассматриваемой точке температура начала и конца мартенситного превращения, а также фазовый состав сплава определяются по термокинетической диаграмме в зависимости от скорости охлаждения (рис. 2). В связи с высокой скоростью роста мартенсита в процессе трансформации (близкой к скорости звука) [20] принято считать, что превращение происходит мгновенно. В этом случае кинетика образующейся мартенситной фазы в зависимости от температуры может быть описана эмпирическим законом Койстинена - Марбургера [10] , (6) где ψ - объемная доля мартенситной фазы; - температура начала мартенситного превращения; b - эмпирический коэффициент (для Ti6Al2V b = -0,05). При исследовании высокоскоростных процессов термической обработки необходимо учитывать зависимость величины критических температур Tкр (температур фазовых превращений) от скорости нагрева. Для титановых сплавов с ростом скорости нагрева характерно смещение критических точек ΔTкр при превращении за счет подавления роста β-фазы, фазовый переход при этом осуществляется только путем образования новых зародышей в области более высоких значений температур [24]. Максимальное смещение при скорости нагрева порядка 105-107 оС/с составляет 90 оС [20]. t, с t, с Т, °С Т, °С а б Рис. 2. Термокинетическая диаграмма (а) и пример термического цикла при ЭМО Ti6Al2V (б) (α, α’, α” - стабильные и метастабильные фазы, образующиеся при распаде β-фазы в процессе непрерывного охлаждения; Ms, Mf - температуры начала и конца мартенситного превращения; - температура превращения) [20, 21] Fig. 2. Continuous cooling diagram (a) and an example of the thermal cycle of the Ti6Al2V alloy during EMT, (b) (α, α ', α " are stable and metastable phases formed during the decay of β-phase in the process of a continuous cooling; Ms, Mf are the temperatures at the beginning and end of the martensite transformation; is the temperature of transformation) [20, 21] При этом зависимость ΔTкр от скорости нагрева для диапазона 200 оС/с < < < 3,5∙104 оС/с близка к линейной и может быть представлена в виде 2. Методика численного решения термоструктурной задачи 2.1. Решение нестационарного уравнения теплопроводности методом конечных элементов В основе решения лежит вариационная формулировка уравнения теплопроводности, полученная методом взвешенных невязок [25, 26]: (7) После процедуры дискретизации и ансамблирования удобно записать уравнение (6) в матричной форме: , (8) где - матрица теплоемкости; , - матрица теплопроводности; - матрица конвективного теплообмена; , - приведенный вектор тепловой нагрузки, соответствующий поверхностному источнику; - приведенный вектор тепловой нагрузки, соответствующий конвективного теплообмену; - приведенный вектор тепловой нагрузки, соответствующий удельному объемному тепловому источнику; - вектор температуры в узлах КЭ-сетки; - вектор производной температуры по времени. Данные матрицы для соответствующих конечных элементов могут быть вычислены с использованием следующих соотношений: (9) где , - матрица функций и производных функций формы соответственно. Разрешающее уравнение можно получить после интегрирования матричного соотношения (8) по времени. Для этого применяется метод Галеркина на исследуемом временном интервале Δt: (10) Интегрирование уравнения (10) на интервале Δt приводит к системе линейных уравнений на i+1 временном шаге следующего вида: (11) Данная схема является полунеявной и обладает безусловной устойчивостью для линейной постановки, кроме того, она меньше подвержена осцилляциям на первых шагах решения, чем другие схемы [9]. 2.2. Учет нелинейных и связанных эффектов К нелинейным эффектам относятся зависимости членов дифференциальных уравнений и граничных условий от температуры, а к эффектам связанности - зависимости от других переменных и дополнительных условий. Для данной задачи, в частности, принимались в расчет следующие нелинейности и связанности: 1) зависимость физических свойств от температуры ; 2) зависимость коэффициента теплоотдачи от температуры и тепловое излучение в бесконечную среду ; 3) скрытая теплота фазовых превращений . При этом конкретный вид функций, определяющих зависимость соответствующих теплофизических свойств от температуры, находится экспериментальным путем. 2.2.1. Нелинейность физических свойств Зависимость теплофизических коэффициентов для сплава Ti-6Al-2V представлена в виде квадратичной сплайн-интерполяции известных экспериментальных данных (рис. 3) [27] вида , где n - число интерполяционных участков; N - квадратичные базисные функции; ar, br - весовые коэффициенты. При расчете на каждом шаге определяется номер интерполяционного участка и вычисляется значение необходимых коэффициентов в рамках соответствующего интервала температур. Необходимо отметить, что на стадии охлаждения происходит только мартенситное превращение без существенного изменения теплоемкости. На аналогичной зависимости отсутствуют пики (кривая 2, рис. 3), расчет при этом ведется при одинаковом для всей временной области шаге Δt. Рис. 3. Расчетные кривые температурных зависимостей теплофизических констант (1 - коэффициента теплоемкости на стадии нагрева; 2 - коэффициента теплоемкости на стадии охлаждения; 3 - коэффициента теплопроводности) для титанового сплава Ti6Al2V, где Ts и Tf - температура соответственно начала и конца фазового перехода Fig. 3. The calculated temperature dependences of thermal constants (1) is for the coefficient of heat capacity at the heating stage; 2) is for the specific heat coefficient at the cooling stage; 3) is for the coefficient of thermal conductivity) for the Ti6Al2V titanium alloy, where Ts and Tf are the temperatures at the beginning and end of the phase transition respectively 2.2.2. Скрытая теплота фазовых переходов Скрытая теплота фазовых превращений, выделяющаяся в процессе превращения, вводится в решение в явном виде как объемный источник тепла Q и добавляется инкрементально на каждом временном шаге с учетом закона Койстинена-Марбургера (рис. 4). где - энтальпия фазового перехода; ψ - объемная доля образующейся мартенситной фазы. Рис. 4. Алгоритм решения нелинейной задачи на текущем временном шаге Fig. 4. Algorithm of solving the nonlinear problem at the current time step 2.2.3. Нелинейность теплоотдачи в окружающую среду Радиационный теплообмен в бесконечную среду эквивалентен распределенному по поверхности стоку тепла вида Данное выражение приводит к эквивалентному коэффициенту теплообмена, зависящему от температуры . 2.3. Методика решения нелинейной задачи Алгоритм решения нелинейной задачи методом конечных элементов в случае нестационарного процесса представлен на рис. 4. Процедура представляет собой последовательное уменьшение невязки основного разрешающего уравнения в общем виде: В нулевом приближении решение на i+1 временном шаге получается путем квадратичной экстраполяции значений, полученных на трех предыдущих временных шагах (см. рис. 4, блок 1): Далее вычисляется значение вектора температур на интервале временного интегрирования в точке решения, определяемое коэффициентом α и зависящее от схемы интегрирования. В частности, для метода Галеркина (см. рис. 4, блок 2). На следующем шаге вычисляется значение основных матриц в точке решения (см. рис. 4, блок 3), затем рассчитывается норма невязки (см. рис. 4, блок 4) и сравнивается с допустимой величиной . При выполнении условия сходимости осуществляется переход к следующему временному шагу. Если условие сходимости не достигнуто, систему нелинейных уравнений необходимо линеаризовать с использованием одного из итерационных методов. Отметим, что решение поставленной нестационарной задачи требует использования весьма малых шагов Δt даже при неявной схеме интегрирования, что связано с высокоскоростными процессами нагрева/охлаждения, быстротекущей трансформации структуры и др. В рамках одного малого временного шага (порядка 0,2 мс) задача может считаться слабо нелинейной, и достаточным будет использование методов линеаризации первого порядка. К таким методам, обладающим линейной скоростью сходимости, относится итерационный метод Пикара [25], который сводит нелинейную задачу прямым итерационным способом к последовательному решению системы линейных уравнений. Для решения полученной системы линейных уравнений в данной работе использовался метод сопряженных градиентов с предобусловливанием. В результате процесс решения можно представить в виде следующего рекуррентного соотношения (см. рис. 4, блок 5): где - разрешающая матрица, а - вектор свободных членов правой части уравнения (11). Для ускорения сходимости используется формула релаксации (см. рис. 4, блок 6): где γ - параметр релаксации. Вычисление основных матриц на каждой итерации осуществляется согласно уравнениям (9) с учетом уточненного вектора температур в пределах текущего временного шага. Итерационный процесс продолжается пока не выполнится условие сходимости. 2.4. Используемые конечные элементы Для формирования КЭ модели используются тетраэдральные симплекс- и изопараметрические 8-узловые КЭ [28, 29]. Для корректного моделирования полупространства вводятся бесконечные элементы [30] (рис. 5). Для указанных бесконечных элементов используется раздельный набор интерполяционных функций, определяющих геометрию элемента и температуру в узлах. Для аппроксимации температуры используются линейные функции вида где k - текущий номер узла в локальной нумерации; - безразмерные локальные в общем случае криволинейные координаты. Функции геометрии элемента, бесконечного в одном из направлений, приведены ниже: Рис. 5. Отображение КЭ на бесконечность Fig. 5. Infinite element mapping scheme 3. Анализ результатов 3.1. Анализ термических циклов и структурообразования при ЭМО Анализ термических циклов при ЭМО, полученных при решении термоструктурной задачи, производится применительно к титановому псевдо-α-сплаву Ti6Al2V. Плотность тока для всех представленных режимов I = 400 A, U = 0,4 В, k = 0,3, ν = 50 Гц, скорость обработки варьируется: На рис. 6 представлен термический цикл точек, расположенных на пересечении плоскостей и на глубинах 0, 0,25 мм, 0,5 мм и 0,75 мм (кривые 1, 2, 3 и 4 соответственно). Температурное воздействие на материал носит достаточно сложный характер, различные точки материала могут несколько раз проходить стадию закалки и отпуска (см. рис. 6, кривая 1). При этом максимальная температура практически достигает температуры плавления, однако носит локальный характер и быстро затухает. Глубина упрочненной зоны, соответствующей закалке, составляет порядка 200 мкм. Скорость охлаждения при этом составляет Окончательная структура материала определяется последним вторичным пиком температурного цикла от удаляющегося источника. Если этот пик доходит до температуры отжига, то структура считается разупрочненной. Рис. 6. Термический цикл и картина микроструктуры на поверхности сплава типа Ti6Al2V при обработке со скоростью (I - упрочненная зона, II - зона отжига) Fig. 6. The thermal cycle and microstructure on the Ti6Al2V alloy surface when it is processed at the rate of (I is the hardened zone, II is the annealing zone) 3.2. Влияние нелинейных эффектов и способа их учета на решение термоструктурной задачи 3.2.1. Зависимость теплофизических свойств от температуры Важным аспектом решения задачи является способ учета нелинейных эффектов. В данной работе произведен сравнительный анализ двух способов учета: 1) квазилинейный, заключающийся в вычислении значения нелинейных членов на основе распределения температуры, полученного на предыдущем временном шаге (явная схема) и последующем решении системы линейных уравнений; 2) нелинейный - на основе метода Пикара. На рис. 7 представлена зависимость относительного отклонения температуры на поверхности при решении по явной схеме при различных скоростях обработки от решения полностью нелинейной задачи при рабочем временном шаге Δt = 0,2 мс. Заметим, что для выбранного временного шага Δt = 0,2 мс максимальное отклонение составляет не более 1 %, однако с увеличением шага Δt отклонение возрастает до 9 % (рис. 9). Пики относительного отклонения приходятся на временные интервалы, соответствующие локальным температурным максимумам, поэтому целесообразно использовать итерационные процедуры именно в этой временной области. Для этого выгоднее будет использовать несколько другой критерий сходимости итерационного процесса, например так называемую бесконечную норму вектора невязки , которая характеризует именно максимальную величину вектора. Рис. 7. Зависимость относительного отклонения температуры от текущего времени обработки при различных скоростях движения источника тепла для квазилинейного способа учета (1 - 2 - ) Fig. 7. Dependence of the relative temperature deviation on the current processing time at different motion speeds of the heat source for the quasi-linear method of accounting (1) ; (2 ) 3.2.2. Скрытая теплота фазовых переходов Скрытая теплота, выделяющаяся в процессе мартенситного превращения, связана с изменением энтальпии в процессе фазового перехода. Для металлических сплавов изменение энтальпии лежит в диапазоне Как показали численные эксперименты, учет теплоты мартенситного превращения незначительно влияет на температурное поле (в пределах ) и несколько увеличивает количество образующегося мартенсита. Рис. 8. Картина скоростей охлаждения и микроструктуры: (а) ; (б) (в) Fig. 8. Microstructure and cooling rates: (a) ; (b) ; (c) Более значительное влияние скрытая теплота оказывает на кинетику фазового перехода, замедляя процесс превращения. Это хорошо видно по величине скорости охлаждения, которая становится меньше на 35 % для значения энтальпии . Область структурообразования сужается при этом максимально на 10 % (рис. 8). 3.2.3. Радиационный теплообмен Охлаждение металла за счет теплового излучения становится обычно значимым фактором при температурах выше 300 °C [20]. Однако, как показали численные эксперименты в данной работе, относительное отклонение температуры составляет менее , что объясняется малым интервалом времени термического воздействия, особенно в области высоких температур. При таких временных рамках процесс перераспределения тепла в объеме металлического тела происходит за счет отвода тепла вглубь материала, и учет любого поверхностного теплообмена становится малозначимым. 3.3. Выбор временного шага При решении термоструктурной задачи величина временного шага выбиралась из следующих соображений: 1) обеспечение достаточной устойчивости решения; 2) минимальное накопление ошибок в ходе решения; 3) адекватное описание процесса высокоскоростной трансформации структуры. Проведено исследование влияния максимального относительного отклонения температуры на поверхности при текущем шаге от температуры при базовом шаге Δt = 0,2 мс при различных способах учета нелинейных и связанных эффектов во всей области скоростей обработки (рис. 9). Рис. 9. Зависимость относительного отклонения температуры от величины временного шага при различных способах учета: 1 - квазилинейный; 2 - полностью нелинейный Fig. 9. Relative temperature dependence on the value of time step for different methods of accounting (1 is quasi-linear, 2 is fully non-linear) Относительное отклонение температуры при увеличении величины временного шага для полностью нелинейного способа учета невелико (в пределах 3 %), однако в таком случае увеличивается степень нелинейности задачи в рамках шага. Следовательно, увеличивается число необходимых итераций при решении нелинейной задачи, вычислительные затраты при этом остаются приблизительно на уровне базового шага Δt. С другой стороны, слишком большое значение временного шага ведет к ошибкам в определении структурных зон и несвоевременному выделению скрытой теплоты фазового превращения. Как показывают численные эксперименты, допустимая величина временного шага лежит в диапазоне 0,2-0,5 мс. Для учета по явной схеме отклонение температуры быстро возрастает: так, при базовом шаге 0,2 мс величина отклонения составляет менее 1 %, тогда как при величине шага 2 мс - уже порядка 9 %. Выводы 1. Приведена обобщенная постановка термоструктурно-механической задачи при моделировании НДС и механического поведения неоднородного тела со структурой, трансформирующейся в ходе высокоэнергетического воздействия. В рамках данной постановки показана необходимость предварительного анализа эффектов, связанных с решением вспомогательных задач описания действующих температурных полей и формирующейся структуры тела. 2. На основе численных экспериментов установлено, что тепловые поля, формирующиеся под действием на титановый сплав Ti6Al2V синусоидального источника тепла, характерного для ЭМО переменным током, носят нелинейный импульсный характер, в результате чего образуется регулярная неоднородная структура. 3. В результате анализа значимости нелинейных и связанных эффектов на решения термоструктурной задачи установлено, что скрытая теплота мартенситного превращения увеличивает температуру в момент превращения на и снижает скорость охлаждения на 35 %. Учет любого поверхностного теплообмена при этом слабо влияет на решение в пределах . 4. Обоснованы важность процедуры выбора величины временного шага, необходимой для адекватного воспроизведения высокоскоростной трансформации структуры. При этом установлено, что допустимая величина временного шага лежит в диапазоне 0,2-0,5 мс. 5. Установлено, что использование квазилинейного способа учета нелинейных и связанных эффектов допустимо при малой величине шага в пределах 0,5 мс, относительное отклонение температуры при этом составляет менее 2 %.

About the authors

V P Bagmutov

Volgograd State Technical University

D S Denisevich

Volgograd State Technical University

I N Zakharov

Volgograd State Technical University

A Yu Ivannikov

Volgograd State Technical University

References

Statistics

Views

Abstract - 325

PDF (Russian) - 92