Математическое и программное обеспечение расчета сильфонных компенсаторов на прочность

Аннотация

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

Полный текст

Сильфонные компенсаторы (рисунок 1) находят все более широкое применение в различных отраслях промышленности. Их применение для компенсации температурных деформаций элементов конструкций дает возможность повысить рабочие параметры оборудования, упростить структуру производственных объектов, сократить габариты, получить в результате существенный технический и экономический эффект. Рисунок 1 – Сильфонные компенсаторы В большинстве случаев сильфонные компенсаторы подвергаются повторным воздействиям нагрузки с ограниченным числом циклов (102—104). Число циклов нагружения определяется количеством пусков, остановок и резких изменений режима работы оборудования за время эксплуатации. В этих условиях ограничение работы компенсатора областью упругих деформаций не позволяет полностью использовать его компенсирующую способность. Для увеличения хода сильфонного компенсатора целесообразно и необходимо допустить пластические деформации в наиболее нагруженных зонах оболочки. Поскольку циклическое упругопластическое деформирование сильфона сопровождается обычно его малоцикловым разрушением, одна из основных задач заключается в обеспечении необходимой долговечности конструкции. Затруднения, связанные с разработкой строгих методов расчета сопротивления конструкций малоцикловому разрушению, и острая необходимость решения важных практических задач обусловили использование инженерных методов, основанных на сочетании результатов упругого расчета с экспериментальными данными. Так, например, расчетная методика, принятая в стандарте ASME по котлам и сосудам давления, основана на принципе деформационного эквивалента, согласно которому условные напряжения, определяемые путем расчета в предположении упругого поведения материала, являются мерой усталостного разрушения [1]. Аналогичный подход к инженерным расчетам малоцикловой прочности ряда конструкций, в том числе и компенсаторов, принят и в нормах расчета [2]. Для практической реализации инженерных методов расчетной оценки малоцикловой прочности сильфонных компенсаторов необходимо располагать надежными и удобными в применении методами и алгоритмами автоматизированного расчета упругого напряженно-деформированного состояния сильфонов. Основным элементом сильфонного компенсатора (рисунок 2) является тонкостенная осесимметричная оболочка (сильфон), состоящая из участков тороидальных оболочек положительной и отрицательной гауссовой кривизны, плавно сопряженных кольцевыми пластинками или участками конических оболочек. Рисунок 2 – Расчетная схема сильфонного компенсатора Присоединительные патрубки представляют собой цилиндрические элементы толщиной h1 и длиной L. Толщина стенки сильфона h и физико-механические характеристики материала в общем случае переменны вдоль меридиана оболочки и являются заданными функциями дуги меридиана s. Сильфон нагружен внутренним давлением q и нагрет до температуры T. Задано также относительное осевое перемещение торцов компенсатора λ. Для конечно-элементного моделирования осесимметричной составной оболочечной конструкции, работающей в условиях осесимметричного термомеханического нагружения, используем оболочечный элемент (рисунок 3). Оболочечный элемент ограничен двумя узлами i и j (i < j), каждый из которых имеет три степени свободы: осевое перемещение u, радиальное перемещение v, угол поворота нормали . Начало координат (s = 0) помещается в узел i. Материал оболочечного элемента изотропный, упругий, подчиняется закону Гука. На рассматриваемый элемент действуют распределенная по срединной поверхности элемента нагрузка qn, нормальная к этой поверхности, а также распределенная по срединной поверхности элемента нагрузка qτ, направленная по касательной к меридиану. Кроме того, элемент может быть нагрет до температуры T = T0 + ky (здесь y – расстояние от срединной поверхности элемента). Температура срединной поверхности изменяется вдоль меридиана по заданному закону . По толщине стенки температура изменяется по линейному закону с заданным градиентом , где ΔT(s) – перепад температуры по толщине стенки. Рисунок 3 – Оболочечный элемент Введем понятие вектора состояния, , где: U – осевое усилие в оболочке, V – радиальное усилие, M – меридиональный изгибающий момент, ζ – осевое перемещение точек координатной поверхности оболочки, ξ – радиальное перемещение точек координатной поверхности оболочки, – угол поворота нормали к срединной поверхности оболочки, r – радиус параллельного круга, |т – символ транспонирования. Напряженно-деформированное состояние оболочечного элемента описывается системой линейных дифференциальных уравнений шестого порядка: , (1) где: – матрица коэффициентов системы, – вектор, учитывающий поверхностную нагрузку и температурные деформации. Основная задача при конечно-элементном анализе составной оболочечной конструкции заключается в формировании разрешающей системы уравнений равновесия узлов системы: , где: – глобальная матрица жесткости всей системы, – вектор узловых перемещений, – вектор внешних узловых усилий, – глобальный вектор узловых усилий всей системы, обусловленных силовым и температурным воздействием. Для построения глобальной матрицы жесткости системы необходимо сформировать матрицы жесткости оболочечных элементов и шпангоутов. Для построения вектора необходимо сформировать для каждого оболочечного элемента векторы краевых обобщенных усилий , обусловленные поверхностной нагрузкой и температурным воздействием. В силу линейности системы дифференциальных уравнений (1) можно установить однозначную зависимость между краевыми усилиями на торцах оболочечного элемента и краевыми перемещениями: . Здесь – матрица жесткости оболочечного элемента размерностью 6*6, зависящая от геометрических и механических характеристик. Столбцы матрицы жесткости представляют собой обобщенные усилия на торцах элемента i и j , вызываемые единичными перемещениями этих торцов при отсутствии поверхностной нагрузки. Вектор является вектором краевых обобщенных усилий, обусловленных поверхностной нагрузкой и температурным воздействием на оболочечный элемент при нулевых смещениях торцов. Алгоритм расчета сильфонного компенсатора содержит следующие операции. 1. Формируем матрицы жесткости для оболочечных элементов системы. С этой целью для каждого элемента решаем шесть однородных краевых задач вида: (2) с граничными условиями, соответствующими единичным перемещениям торцов элемента. 2. Определяем векторы краевых обобщенных усилий для оболочечных элементов, выполняя для каждого элемента решение неоднородной краевой задачи (1) с нулевыми граничными условиями. Для решения систем дифференциальных уравнений (1), (2) применяем метод ортогональной прогонки. 3. Формируем глобальную матрицу жесткости для составной оболочечной конструкции. Для формирования матрицы используем матрицы жесткости оболочечных элементов системы. Размерность матрицы жесткости равна 3 NU*3 NU, где NU – общее число узлов в системе. 4. Формируем глобальный вектор узловых усилий из векторов краевых обобщенных усилий, построенных для оболочечных элементов системы. Размерность вектора равна 3*NU. 5. По заданным внешним усилиям формируем вектор внешних узловых усилий для всей системы. Размерность вектора равна 3*NU. Отметим, что вектор содержит не только заданные значения внешних силовых факторов, но и реакции наложенных связей (реакции опор). 6. Выполняя операцию составления ансамбля [3], формируем разрешающую систему уравнений метода конечных элементов: , где: – глобальный вектор нагрузки; – матрица жесткости всей системы; – вектор узловых перемещений; – вектор внешних узловых усилий; – глобальный вектор узловых усилий, обусловленных силовым и температурным воздействием на рассматриваемую конструкцию. 7. Модифицируем матрицу жесткости системы и глобальный вектор нагрузки с учетом граничных условий. Если на узел n наложена связь с номером i, (i =1,2,3), то все элементы строки j = матрицы заменяем нулями, кроме диагонального элемента, который заменяем единицей. Элемент j = вектора узловой нагрузки заменяем нулем. В результате получаем модифицированную матрицу жесткости и модифицированный вектор нагрузки . 8. Выполняем решение системы уравнений: . В результате решения системы уравнений находим все компоненты вектора узловых перемещений . 9. Определяем реакции связей, наложенных на систему. С этой целью возвращаемся к немодифицированной матрице жесткости системы . Если на узел n наложена связь с номером i, (i =1,2,3), то формируем вектор из элементов строки j = 3(n-1)+i этой матрицы. Размерность вектора равна 3*NU. Вычисляем произведение векторов P = . Реакцию связи вычисляем по формуле: , где: – j-я компонента глобального вектора внешних воздействий на оболочечную конструкцию. 10. Выполняем расчет напряженно-деформированного состояния оболочечных элементов конструкции. С этой целью для каждого оболочечного элемента методом ортогональной прогонки выполняем решение краевой задачи для системы дифференциальных уравнений (1) с граничными условиями: при s = 0: ζ = ζi , ξ = ξi , , при s = sk: ζ = ζj , ξ = ξj , , где: ζi , ξi , – перемещения i-го узла оболочечного элемента; ζi , ξi , – перемещения j-го узла оболочечного элемента. Все узловые перемещения определены на предыдущем этапе. Численная реализация предложенного метода и алгоритма компьютерного анализа напряженно-деформированного состояния компенсаторов сильфонного типа осуществлена в виде программного обеспечения. Программный комплекс «Compensator» имеет модульную структуру, функционирует в операционных системах Windows XP/7, предоставляет пользователю удобный графический интерфейс. Предъявляет умеренные требования к аппаратному обеспечению. Позволяет выполнять графический синтез и численный анализ напряженно-деформированного состояния исследуемых конструкций, осуществлять поиск оптимальных конструктивных параметров по критерию минимума интенсивности напряжений. Вывод результатов численного анализа осуществляется в табличной и графической формах. В качестве примера выполним с помощью программы «Compensator» расчет типовой конструкции сильфонного компенсатора = 400 мм по следующим исходным данным: = 400 мм; = 460 мм; = = 10 мм; L = 50 мм; = 5 мм. Материал сильфона – аустенитная нержавеющая сталь Х18Н10Т. Величина внутреннего давления q = 0,5 МПа. Осевой ход компенсатора Δ = 50 мм. Компенсатор установлен в тепловой контур с предварительным натягом λ0 = Δ/2 = 25 мм. Результаты численного анализа представлены на рисунке 4. Рисунок 4 – Интенсивность напряжений в исходной конструкции компенсатора Показаны графики интенсивности условных упругих напряжений в точках внутренней и наружной поверхности сильфонного компенсатора при относительном смещении торцов λ=-25 мм и внутреннем давлении q = 0,5 МПа. Результаты расчета показывают, что напряжения в сильфоне распределены крайне неравномерно, при этом напряжения на внутренней поверхности оболочки существенно превышают напряжения на наружной поверхности. Максимальные значения интенсивности условных упругих напряжений составляют 1232 МПа во впадинах гофров сильфона. Для оценки долговечности рассматриваемой конструкции компенсатора используем зависимость [4] (3) где: – амплитуда интенсивности условных упругих напряжений, МПа, в наиболее нагруженном сечении сильфона; – разрушающее число циклов нагружения. Для аустенитной нержавеющей стали m = 0,35; C = 9400. В данном случае при = 1232 МПа долговечность компенсатора N = 332 цикла. С целью увеличения ресурса изделия оптимизируем его конструкцию. Увеличиваем на 10 мм ширину кольцевых пластинок, соединяющих тороидальные оболочечные элементы положительной и отрицательной гауссовой кривизны. Остальные конструктивные параметры компенсатора оставляем без изменений. Высота гофров при этом составит 40 мм. Наружный диаметр сильфона станет равным 480 мм. Результаты численного анализа оптимизированной конструкции сильфонного компенсатора представлены на рисунке 5. Рисунок 5 – Интенсивность напряжений в оптимизированной конструкции компенсатора Максимальные значения интенсивности условных упругих напряжений существенно снизились и составляют теперь 893 МПа на внутренней поверхности тороидальных элементов отрицательной гауссовой кривизны. Выполненные конструктивные изменения позволили в 2,5 раза увеличить ресурс изделия. В соответствии с формулой (3) долговечность оптимизированной конструкции компенсатора при амплитуде интенсивности условных упругих напряжений = 893 составляет = 832 цикла.