Известия РАН. Механика твердого тела, 2022, № 2, стр. 14-24

ОБ ОПТИМИЗАЦИИ КИНЕМАТИЧЕСКОЙ СХЕМЫ ФОРМООБРАЗОВАНИЯ ОБТЯЖКОЙ ТОНКИХ УПРУГОПЛАСТИЧЕСКИХ ОБШИВОК

К. С. Бормотин a*, А. А. Буренин b**, А. А. Кривенок c***

a Комсомольский-на-Амуре государственный университет
Комсомольск-на-Амуре, Россия

b Институт машиноведения и металлургии ДВО РАН
Комсомольск-на-Амуре, Россия

c Комсомольский-на-Амуре авиационный завод им. Ю.А. Гагарина
Комсомольск-на-Амуре, Россия

* E-mail: cvmi@knastu.ru
** E-mail: aab@imim.ru
*** E-mail: baikal-190@yandex.ru

Поступила в редакцию 17.10.2021
После доработки 20.10.2021
Принята к публикации 21.10.2021

Полный текст (PDF)

Аннотация

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

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

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

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

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

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

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

2. Формулировка задач оптимального управления при формообразовании тонкостенных конструкций. Пусть $V \subset {{R}^{3}}$ – область деформируемого тела с границей S. Область задания краевых смещений – ${{S}_{b}}$ (рис. 1). Контактная поверхность пуансона с деформируемым телом обозначается через ${{S}_{c}}$ (${{S}_{c}} \subset S$). Обозначим через ${\mathbf{u}} = ({{u}_{1}},{{u}_{2}},{{u}_{3}})$, ${\mathbf{\tilde {u}}} = ({{\tilde {u}}_{1}},{{\tilde {u}}_{2}},{{\tilde {u}}_{3}})$ – векторы текущих, остаточных перемещений деформируемого тела, t – параметр деформирования. Моделирование обтяжки листа на пуансоне выполняется с учетом малых деформаций, но больших перемещений и поворотов (общая Лагранжева формулировка (TL) [6]).

Рис. 1.

Модель обтяжного пресса (1– заготовка из листа, 2 – пуансон, 3 – контрольная точка).

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

(2.1)
${{J}_{1}}({\mathbf{\dot {u}}}) = {{W}_{c}} + a({\mathbf{\dot {u}}},{\mathbf{\dot {u}}})\quad {\text{при}}\quad {\mathbf{\dot {u}}}{{{\text{|}}}_{{{{S}_{b}}}}} = {\mathbf{\dot {u}}}{\kern 1pt} *,\quad J{}_{2}({\mathbf{\dot {\tilde {u}}}}) = {{W}_{c}} + a({\mathbf{\dot {\tilde {u}}}},{\mathbf{\dot {\tilde {u}}}})$
где ${{W}_{c}}$ – контактный потенциал [6, 7]; потенциальные формы
$a({\mathbf{\dot {u}}},\delta {\mathbf{\dot {u}}}) = \int\limits_V {\frac{{\partial E(\dot {u}_{{i,j}}^{{}})}}{{\partial {{{\dot {u}}}_{{i,j}}}}}\delta {{{\dot {u}}}_{{i,j}}}} dV,\quad E(\dot {u}_{{i,j}}^{{}}) = \frac{1}{2}{{c}_{{ijkl}}}{{\dot {\varepsilon }}_{{ij}}}{{\dot {\varepsilon }}_{{kl}}} - {{c}_{{ijkl}}}{{\dot {\varepsilon }}_{{ij}}}\dot {\varepsilon }_{{kl}}^{p} + \frac{1}{2}{{\sigma }_{{ij}}}\dot {u}_{{k,i}}^{{}}\dot {u}_{{k,j}}^{{}}$
$a({\mathbf{\dot {\tilde {u}}}},\delta {\mathbf{\dot {\tilde {u}}}}) = \int\limits_V {\left( {\frac{{\partial{ \tilde {E}}(\dot {\tilde {u}}_{{i,j}}^{{}})}}{{\partial {{{\dot {\tilde {u}}}}_{{i,j}}}}}} \right)\delta {{{\dot {\tilde {u}}}}_{{i,j}}}} dV,\quad \tilde {E}(\dot {\tilde {u}}_{{i,j}}^{{}}) = \frac{1}{2}{{c}_{{ijkl}}}{{\dot {\tilde {\varepsilon }}}_{{ij}}}{{\dot {\tilde {\varepsilon }}}_{{kl}}} - {{c}_{{ijkl}}}{{\dot {\tilde {\varepsilon }}}_{{ij}}}\dot {\varepsilon }_{{kl}}^{p} + \frac{1}{2}{{\tilde {\sigma }}_{{ij}}}\dot {\tilde {u}}_{{k,i}}^{{}}\dot {\tilde {u}}_{{k,j}}^{{}}$
${{c}_{{ijkl}}}$ – компоненты тензора упругих констант; ${{\dot {\varepsilon }}^{p}}_{{ij}}$ – компоненты скоростей пластических деформаций (${{\dot {\varepsilon }}^{p}}_{{ij}} = \lambda \frac{{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{f} }}{{\partial {{\sigma }_{{ij}}}}}$, $\lambda > 0$, уравнение $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{f} = 0$ определяет поверхность текучести в пространстве компонент девиатора тензора напряжений) [6], ${{\dot {\varepsilon }}_{{ij}}}$, ${{\dot {\tilde {\varepsilon }}}_{{ij}}}$ – компоненты скорости текущих и остаточных деформаций Грина–Лагранжа, ${{\sigma }_{{ij}}}$, ${{\tilde {\sigma }}_{{ij}}}$ – компоненты текущего и остаточного второго тензора напряжений Пиола–Кирхгофа, ${{u}_{{i,j}}} = \frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}}$, $i,j,k,l = 1,2,3$.

Второй тензор напряжений Пиола–Кирхгофа может быть определен через скорости деформаций

${{\dot {\sigma }}_{{ij}}} = {{c}_{{ijkl}}}({{\dot {\varepsilon }}_{{kl}}} - \dot {\varepsilon }_{{kl}}^{p})$

Обобщая вариационные принципы с функционалами (2.1) для произвольных деформаций тела с помощью текущей Лагранжевой формулировки или UL-подхода потенциальные формы примут вид

$a({\mathbf{v}},\delta {\mathbf{v}}) = \int\limits_{{{V}^{t}}} {\frac{{\partial W({v}_{{i,j}}^{{}})}}{{\partial {v}_{{i,j}}^{{}}}}\delta {v}_{{i,j}}^{{}}} dV,\quad a({\mathbf{\tilde {v}}},\delta {\mathbf{\tilde {v}}}) = \int\limits_{{{V}^{t}}} {\left( {\frac{{\partial{ \tilde {W}}({\tilde {v}}_{{i,j}}^{{}})}}{{\partial {{{{\tilde {v}}}}_{{i,j}}}}}} \right)\delta {{{{\tilde {v}}}}_{{i,j}}}} dV$
$\begin{gathered} W\left( {{{{v}}_{{i,j}}}} \right) = \frac{1}{2}{{c}^{t}}_{{ijkl}}{{d}_{{ij}}}{{d}_{{kl}}} - c_{{ijkl}}^{t}{{d}_{{ij}}}d_{{kl}}^{p} + \frac{1}{2}{{p}_{{ij}}}{v}_{{k,i}}^{{}}{v}_{{k,j}}^{{}} \\ \tilde {W}({\tilde {v}}_{{i,j}}^{{}}) = \frac{1}{2}c_{{ijkl}}^{t}{{{\tilde {d}}}_{{ij}}}{{{\tilde {d}}}_{{kl}}} - c_{{ijkl}}^{t}{{{\tilde {d}}}_{{ij}}}d_{{kl}}^{p} + \frac{1}{2}{{{\tilde {p}}}_{{ij}}}{\tilde {v}}_{{k,i}}^{{}}{\tilde {v}}_{{k,j}}^{{}} \\ \end{gathered} $
где ${\mathbf{v}}$, ${\mathbf{\tilde {v}}}$ – векторы текущих и остаточных скоростей перемещений, ${v}_{{i,j}}^{{}}$, ${\tilde {v}}_{{i,j}}^{{}}$ – компоненты тензоров градиентов скоростей, ${{d}_{{ij}}}$, ${{\tilde {d}}_{{ij}}}$ – компоненты тензоров скоростей деформаций, ${{p}_{{ij}}}$, ${{\tilde {p}}_{{ij}}}$ – компоненты текущих и остаточных тензоров напряжений Коши. В данном случае используется текущая конфигурация в качестве отсчетной, пространственное дифференцирование, а интеграл берется по текущей конфигурации.

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

(2.2)
${{J}_{3}} = \,\int\limits_S {\sqrt {\sum\limits_i {{{{({{{\tilde {u}}}_{i}}(T) - \tilde {u}_{{_{i}}}^{*})}}^{2}}} } dS} \, \to \inf $
где $\tilde {u}_{i}^{*}$ – заданные компоненты остаточных перемещений листа до пуансона, обеспечивающих прилегание листа на поверхности пуансона.

Задача оптимального управления, включающая уравнения механики деформируемого твердого тела и функционал оптимизации

(2.3)
${{J}_{4}} = \,\int\limits_0^T {\int\limits_V {{{\sigma }_{{ij}}}\dot {\varepsilon }_{{ij}}^{p}dVdt} } \, \to \inf $
определяет кинематическую схему движения зажимов с минимальной поврежденностью материала в пластичности.

При реализации кинематических схем обтяжки с этапами, чередующими нагружение и разгрузку, происходит циклическое упруго-пластическое деформирование. В этом случае возможно появление эффекта Баушингера в материале заготовки [8] и при анализе траекторий деформирования необходимо учитывать анизотропное упрочнение [9, 10]. Идеальный эффект Баушингера определяется принципом Мазинга, который сводится к удвоенному начальному пределу текучести равному разности пределов текучести растяжения и сжатия [8, 9].

Изотропное упрочнение определяется условием пластичности Хубера–Мизеса [9]:

(2.4)
где ${{s}_{{ij}}}$ – компоненты девиатора тензора напряжений, $d\bar {\varepsilon }_{{}}^{p}$ – интенсивность приращений пластических деформаций, функция $\Phi $ определяется по диаграмме растяжения материала.

В случае кинематического (трансляционного) упрочнения уравнение поверхности нагружения определяется в виде:

(2.5)
где ${{\rho }_{{ij}}}$ – компоненты девиатора добавочного напряжения, $({{s}_{{ij}}} - {{\rho }_{{ij}}})$ – компоненты девиатора активного напряжения [9], $\sigma _{T}^{0}$ – предел текучести, который является постоянной величиной. Координаты центра поверхности пластичности зависят от пластической деформации, например, с помощью простейшей зависимости ${{\rho }_{{ij}}} = с\varepsilon _{{ij}}^{p}$ (c – постоянная для определенного материала), предложенной А.Ю. Ишлинским.

Для согласования результатов теории и эксперимента целесообразно комбинировать расширение и жесткое смещение поверхности пластичности (изотропное и трансляционное упрочнение), в этом случае условие пластичности имеет вид [9]:

(2.6)

3. Численный метод оптимизации кинематической схемы формообразования панелей. С учетом дискретизации вариации (2.1) формируются уравнения конечных элементов для решения контактной задачи [6, 7]

(3.1)
${}^{{t + dt}}{{{\mathbf{K}}}^{{(r - 1)}}}\Delta {{{\mathbf{U}}}^{{(r)}}} = {}^{{t + dt}}{\mathbf{R}}_{{}}^{{(r - 1)}},\quad {}^{{t + dt}}{{{\mathbf{\tilde {K}}}}^{{(r - 1)}}}\Delta {{{\mathbf{\hat {U}}}}^{{(r)}}} = {}^{{t + dt}}{\mathbf{\hat {R}}}_{{}}^{{(r - 1)}}$
где ${}^{{t + dt}}{{{\mathbf{K}}}^{{(r - 1)}}}$, ${}^{{t + dt}}{{{\mathbf{\tilde {K}}}}^{{(r - 1)}}}$ – матрицы касательной жесткости (в матрицах включены дополнительные элементы, образующиеся от контактных ограничений), ${}^{{t + dt}}{\mathbf{R}}_{{}}^{{(r - 1)}}$, ${}^{{t + dt}}{\mathbf{\hat {R}}}_{{}}^{{(r - 1)}}$ – векторы внешних и внутренних сил. Верхние индексы величин $t + dt$ – значение параметра нагружения, $(r - 1)$ – номер итерации метода Ньютона–Рафсона. Остаточные узловые перемещения определяются по формуле ${\mathbf{\tilde {U}}} = {\mathbf{U}} + {\mathbf{\hat {U}}}$.

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

${{\sigma }_{T}} = {{\sigma }_{T}}^{0} + H{{\bar {\varepsilon }}^{p}},\quad {{\bar {\varepsilon }}^{p}} = \int {{{{\dot {\bar {\varepsilon }}}}^{p}}dt} = \int {\sqrt {\frac{2}{3}\dot {\varepsilon }_{{ij}}^{p}\dot {\varepsilon }_{{ij}}^{p}} dt} $
где модуль $H = \frac{{\Delta \sigma }}{{\Delta {{\varepsilon }^{p}}}}$ находится из одноосной зависимости напряжения от деформации.

При кинематическом упрочнении добавочные напряжения зависят от текущего напряжения и накопленной интенсивности пластической деформации [11]:

$\Delta {{\rho }_{{ij}}} = \sqrt {\frac{2}{3}} H\Delta {{\varepsilon }^{p}}\frac{{\left( {{{s}_{{ij}}} - {{\rho }_{{ij}}}} \right)}}{{{{{\left( {\left( {{{s}_{{ij}}} - {{\rho }_{{ij}}}} \right)\left( {{{s}_{{ij}}} - {{\rho }_{{ij}}}} \right)} \right)}}^{{\frac{1}{2}}}}}}$

Для комбинированной модели упрочнения вводится параметр $\beta \in [0,1]$ и предел текучести, приращения добавочных напряжений определяются по формулам [11]:

${{\sigma }_{T}} = \sqrt {\frac{2}{3}} (\sigma _{T}^{0} + \left( {1 - \beta } \right)H{{\bar {\varepsilon }}^{p}}),\quad \Delta {{\rho }_{{ij}}} = \sqrt {\frac{2}{3}} \beta H\Delta {{\varepsilon }^{p}}\frac{{({{s}_{{ij}}} - {{\rho }_{{ij}}})}}{{{{{\left( {\left( {{{s}_{{ij}}} - {{\rho }_{{ij}}}} \right)\left( {{{s}_{{ij}}} - {{\rho }_{{ij}}}} \right)} \right)}}^{{\frac{1}{2}}}}}}$

Данная модель при β = 1 соответствует кинематическому упрочнению и при β = 0 изотропному упрочнению.

Для приближенного решения задачи оптимального управления интервал $[0,T]$ разбивается на N частей: $0 = {{t}_{0}} < {{t}_{1}} < {{t}_{2}} < \ldots < {{t}_{N}} = T$. Используя дискретные по времени уравнения пошаговой процедуры интегрирования (3.1), при условии $dt \leqslant \Delta t = {{t}_{{n + 1}}}$tn, вычисляются напряжения, деформации и перемещения.

Расстояние контрольной точки листа до пуансона после деформирования до момента t и разгрузки определяется формулой: $l(t) = \sqrt {\sum\limits_i {{{{({{{\tilde {u}}}_{i}}(t) - \tilde {u}_{i}^{*})}}^{2}}} } $. Приращение длины после деформирования до момента $t + \Delta t$ и разгрузки определяется.

$\Delta l(t) = \frac{{\partial l(t)}}{{\partial t}}\Delta t = \frac{{\sum\limits_i {({{{\tilde {u}}}_{i}}(t) - \tilde {u}_{i}^{*})\Delta {{{\tilde {u}}}_{i}}} }}{{\sqrt {\sum\limits_i {{{{({{{\tilde {u}}}_{i}}(t) - \tilde {u}_{i}^{*})}}^{2}}} } }}$

Таким образом, дискретный аналог функционала (2.2) представляется в виде

(3.2)
${{\bar {J}}_{3}} = \sum\limits_{n = 0}^{N - 1} {\sum\limits_S {\frac{{\left( {\sum\limits_i {\left( {\sum\limits_{m = 0}^{n - 1} {\sum\limits_{t = {{t}_{m}}}^{{{t}_{{m + 1}}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} - {\mathbf{\tilde {U}}}_{i}^{*}} } \right)\sum\limits_{t = {{t}_{n}}}^{{{t}_{{n + 1}}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} } } \right)}}{{{{{\left( {{{{\sum\limits_i {\left( {\sum\limits_{m = 0}^{n - 1} {\sum\limits_{t = {{t}_{m}}}^{{{t}_{{m + 1}}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} - {\mathbf{\tilde {U}}}_{i}^{*}} } \right)} }}^{2}}} \right)}}^{{1/2}}}}}} } \to \inf $
где $\Delta {{{\mathbf{\tilde {U}}}}_{i}}(t)$ – компоненты приращения остаточных перемещений узлов, вызванных от приращения текущих перемещений ($\Delta {\mathbf{\tilde {U}}}(t) = \left\{ {\Delta {{{{\mathbf{\tilde {U}}}}}_{1}}(t),\Delta {{{{\mathbf{\tilde {U}}}}}_{2}}(t),\Delta {{{{\mathbf{\tilde {U}}}}}_{3}}(t)} \right\}$), внутренняя сумма при n = 0 равна нулю. Начальные приращения остаточных перемещений равны остаточным перемещениям, полученным после деформирования до момента ${{t}_{1}}$ и разгрузки $\sum\limits_{t = {{t}_{0}}}^{{{t}_{1}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} = \Delta {{{\mathbf{\tilde {U}}}}_{i}}({{t}_{1}}) = {{{\mathbf{\tilde {U}}}}_{i}}({{t}_{1}})$, остальные $\sum\limits_{t = {{t}_{n}}}^{{{t}_{{n + 1}}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} = \Delta {{{\mathbf{\tilde {U}}}}_{i}}({{t}_{{n + 1}}}) = {{{\mathbf{\tilde {U}}}}_{i}}({{t}_{{n + 1}}}) - {{{\mathbf{\tilde {U}}}}_{i}}({{t}_{n}})$, $n \geqslant 1$.

Дискретный аналог функционала (2.3) принимает вид

(3.3)
$\,{{\bar {J}}_{4}} = \sum\limits_{n = 0}^{N - 1} {\sum\limits_{t = {{t}_{n}}}^{{{t}_{{n + 1}}}} {\,\sum\limits_V {{}^{t}{{\sigma }_{{ij}}}\Delta \varepsilon _{{ij}}^{p}} } } \to \inf $

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

Для применения алгоритма оптимизации удобно кривые траекторий деформирования при обтяжке аппроксимировать дугами эллипса, тогда в декартовой системе координат вектор-функция ${\mathbf{U}}(t) = ({{{\mathbf{U}}}_{x}}(t),{{{\mathbf{U}}}_{y}}(t),0)$ перемещений точек краев листа на границе ${{S}_{b}}$ и зажима задается радиусом и углом в виде

${{{\mathbf{U}}}_{x}}(t) = {{{\mathbf{U}}}_{x}}(r_{x}^{0} + {{f}_{1}}(t) \cdot r_{x}^{*},{{\varphi }^{0}} + {{f}_{2}}(t) \cdot \varphi {\text{*}}),\quad {{{\mathbf{U}}}_{y}}(t) = {{{\mathbf{U}}}_{y}}(r_{y}^{0} + {{f}_{1}}(t) \cdot r_{y}^{*},{{\varphi }^{0}} + {{f}_{2}}(t) \cdot \varphi {\text{*}})$
где ${\mathbf{U}}_{x}^{*} = {{{\mathbf{U}}}_{x}}(r_{x}^{0} + r_{x}^{*},{{\varphi }^{0}} + \varphi _{{}}^{*})$, ${\mathbf{U}}_{y}^{*} = {{{\mathbf{U}}}_{y}}(r_{y}^{0} + r_{y}^{*},{{\varphi }^{0}} + \varphi _{{}}^{*})$ – конечные значения перемещений, $r_{x}^{0}$,$r_{y}^{0}$, ${{\varphi }^{0}}$ – начальные значения радиусов и угла, $r_{x}^{*}$, $r_{y}^{*}$, φ* – заданные значения изменения радиусов и угла, ${{f}_{1}}(t) \in [0,1]$, ${{f}_{2}}(t) \in [0,1]$. Необходимо найти оптимальную функциональную зависимость компонент перемещений ${{{\mathbf{U}}}_{y}} = f({{{\mathbf{U}}}_{x}})$.

Для использования рекуррентных соотношений Беллмана строится сетка в пространствах $({{r}_{x}},\varphi )$ и $({{r}_{y}},\varphi )$. Шаги по радиусам ${{r}_{x}}$,${{r}_{y}}$, заданы и равны $\Delta {{r}_{x}}$, $\Delta {{r}_{y}}$, по углу $\varphi $$\Delta \varphi $.

При решении задачи деформирования листа при перемещении краев с положения $a \in \,{{H}_{n}}$ в $b\, \in \,{{H}_{{n + 1}}}$ по некоторой функции ${{{\mathbf{U}}}_{y}} = f({{{\mathbf{U}}}_{x}})$ с помощью уравнений (3.1) определяются на отрезке времени $[{{t}_{n}},{{t}_{{n + 1}}}]$ перемещения, деформации и напряжения листа. ${{H}_{n}}$ – шкала состояний с M узлами. Уравнения (3.1) и критерий оптимизации (3.2) или (3.3) образуют дискретную задачу оптимального управления. В такой постановке строится функция Беллмана и задача решается методом динамического программирования [12, 13]. Задача

$A(a,b,{\mathbf{U}}(t),{{t}_{{n + 1}}}) = \sum\limits_S {\frac{{\left( {\sum\limits_i {\left( {\sum\limits_{t = {{t}_{0}}}^{{{t}_{n}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} - {\mathbf{\tilde {U}}}_{i}^{*}} \right)\sum\limits_{t = {{t}_{n}}}^{{{t}_{{n + 1}}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} } } \right)}}{{{{{\left( {{{{\sum\limits_i {\left( {\sum\limits_{t = {{t}_{0}}}^{{{t}_{n}}} {\,\Delta {{{{\mathbf{\tilde {U}}}}}_{i}}(t)} - {\mathbf{\tilde {U}}}_{i}^{*}} \right)} }}^{2}}} \right)}}^{{1/2}}}}}} \to \inf $
или
$A(a,b,{\mathbf{U}}(t),{{t}_{{n + 1}}}) = \sum\limits_{t = {{t}_{n}}}^{{{t}_{{n + 1}}}} {\sum\limits_V {{}^{t}{{\sigma }_{{ij}}}\Delta {{\varepsilon }_{{ij}}}^{p}} \,} \to \inf ,\quad a = {\mathbf{U}}({{t}_{n}}),\quad b = {\mathbf{U}}({{t}_{{n + 1}}})$
будет определять элементарную операцию, соединяющую точки a и b. Через ${{\Delta }_{{n + 1}}}(a,b)$ обозначается множество всех управлений ${\mathbf{U}}(t)$ (функций, связывающих точки a и b) для которых выполняется (3.1).

Если все ближайшие точки всех соседних шкал попарно соединены элементарными операциями, то величина

$\sum\limits_{n = 0}^{N - 1} {{{L}_{{n + 1}}}({{{\mathbf{U}}}_{{n{{p}_{n}}}}},{{{\mathbf{U}}}_{{n + 1,{{p}_{{n + 1}}}}}})} $
представляет собой значение критерия на управлении ${\mathbf{U}}(t) = {{{\mathbf{U}}}_{{n + 1}}}(t)$, ${{t}_{n}} \leqslant t \leqslant {{t}_{{n + 1}}}$, $n = 0,1,...,N - 1$, ${\mathbf{U}}(0) = 0$, где ${{L}_{{n + 1}}}(a,b) = \mathop {\inf }\limits_{{\mathbf{U}} \in {{\Delta }_{{n + 1}}}(a,b)} A(a,b,{\mathbf{U}}(t),{{t}_{{n + 1}}})$.

Таким образом, необходимо отыскать минимум по всевозможным наборам точек $({{{\mathbf{U}}}_{{0{{p}_{0}}}}},{{{\mathbf{U}}}_{{1{{p}_{1}}}}},....,{{{\mathbf{U}}}_{{N{{p}_{N}}}}})$, ${{{\mathbf{U}}}_{{np}}} \in \,{{H}_{n}}$, $n = 0,...,N$.

Пусть

${{C}_{s}}(a) = \inf \left\{ {\sum\limits_{n = 1}^s {{{L}_{n}}\left( {{{{\mathbf{U}}}_{{n - 1{{p}_{{n - 1}}}}}},{{{\mathbf{U}}}_{{n{{p}_{n}}}}}} \right)} } \right\}$
где нижняя грань берется по всем наборам точек $({{{\mathbf{U}}}_{{0{{p}_{0}}}}},{{{\mathbf{U}}}_{{1{{p}_{1}}}}},....,{{{\mathbf{U}}}_{{s{{p}_{s}}}}})$, $a = {{{\mathbf{U}}}_{{s{{p}_{s}}}}}$. Таким образом, ${{C}_{s}}(a)$ выражает собой кратчайшее расстояние между точкой $a \in \,{{H}_{s}}$ и шкалой ${{H}_{0}}$. Функция ${{C}_{s}}(a)$ удовлетворяет следующим рекуррентным соотношениям

${{C}_{s}}(b) = \mathop {\inf }\limits_{a \in {{H}_{{s - 1}}}} \left\{ {{{L}_{s}}(a,b) + {{C}_{{s - 1}}}(a)} \right\},\quad s = 1,...,N,\quad {{C}_{0}}(a) = 0$

4. Численные результаты решения задачи оптимизации и экспериментальные данные. Решение задачи моделирования и оптимизации процесса обтяжки листа проводится методом конечных элементов в MSC.Marc, используя созданные пользовательские подпрограммы. В силу симметрии рассматривается только правая часть листа (рис. 1). Заготовка толщиной 3 мм и длинной 1000 мм задается набором 4 угольных элементов для 2D анализа (условия плоской деформации). Контактное взаимодействие заготовки с пуансоном на этапе нагружения определяется с коэффициентом трения равным 0.1. Заготовка имеет свойства материала 1163Т: модуль Юнга Е = 7073 кГ/мм2, коэффициент Пуассона $\nu = 0.34$, предел текучести ${{\sigma }_{Т}} = 30$ кГ/мм2, упрочнение задается степенной функцией ${{\sigma }_{Т}} = A{{({{\varepsilon }_{0}} + {{\bar {\varepsilon }}^{p}})}^{g}}$ (${{\varepsilon }_{0}}$ – деформации, соответствующие начальному пределу текучести, A, g – константы материала).

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

На рис. 2 представлены варианты монотонных траекторий деформирования от точки A до B для правого края листа при обтяжке с сеткой при N = 4, M = 1 (рис. 2, a) и N = M = 2 (рис. 2, b), а также оптимальный путь (выделен жирным) по критерию (2.2). Значения по осям координат на рисунках указаны в мм. На рис. 3 представлено отличие траекторий, заданных дугами эллипса, и траектории, выполняющих изгиб по пуансону без сжатия и растяжения нейтрального слоя листа (кривая с маркерами). Как можно увидеть, оптимальная траектория отклоняется от кривой эллипса (рис. 2, b).

Рис. 2.

Возможные варианты и оптимальная траектория движения зажима пресса при деформировании заготовки: (a) – при N = 4, M = 1; (b) – при N = M = 2.

Рис. 3.

Сравнение траекторий движения зажима пресса при деформировании заготовки для моделирования и эксперимента.

Анализ кинематических схем деформирования среди немонотонных траекторий в пространстве $({{r}_{x}},\varphi )$ и различных вариантов упрочнения представлены на рис. 4.

Рис. 4.

Возможные варианты и оптимальная траектория движения зажима пресса при деформировании заготовки по критерию (2.2) при N = M = 2 – (a), при N = M = 4 – (b), по критерию (2.3) при N = M = 4 – (c).

Решение задачи оптимизации с сеткой N = M = 4 на рис. 4, b и рис. 4, c совпадает результатами при N = M = 3, а по критерию (2.3) и при N = M = 2.

При учете каждой модели пластического упрочнения (2.4), (2.5), (2.6) с коэффициентом $\beta = 0.25$ решения, рассмотренных задач, не изменялись.

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

Натурные эксперименты проведены на обтяжном прессе Т-600 с использованием оснастки для формообразования обшивки детали самолета (рис. 1). Заготовки представляют собой полосы из листа материала 1163Т шириной 200 мм, длиной 2000 мм и толщиной 3 мм. Результаты экспериментов проведены для трех траекторий нагружения (жирная, штриховая и штрихпунктирная кривые на рис. 2, a) и представлены на рис. 5. Данные траектории соответствуют трем схемам нагружения: изгиб–растяжение, растяжение–изгиб–растяжение и растяжение–изгиб. В этих схемах траектория изгиба по пуансону вычислена при условиях отсутствия сжатия и растяжения нейтрального слоя листа. Данная траектория отмечена кривой с маркерами на рис. 3.

Рис. 5.

Образцы заготовок после формообразования на прессе Т-600 по разным схемам обтяжки.

На рис. 6 представлены результаты моделирования задач обтяжки в зависимости от формулировки (общей Лагранжевой и текущей Лагранжевой) и от учета трения при контакте заготовки с пуансоном в сравнении с экспериментальными данными. Для сравнения формулировок TL и UL расчеты проводились с одинаковым количеством постоянных шагов и одинаковыми условиями сходимости. Образующиеся деформации по разным формулировкам незначительно отличаются и не достигают 5%. Согласно рис. 6 результаты моделирования качественно совпадают с экспериментами, а отклонение объясняется отличием траектории изгиба в схемах обтяжки (при моделировании изгиб выполняется по дуге эллипса).

Рис. 6.

Сравнение с экспериментальными данными вычисленной в разных формулировках поверхности разгруженного листа: (a) – после растяжения, изгиба и растяжения; (b) – после растяжения и изгиба (1 – пуансон, 2 – эксперимент, 3 – UL, 4 – TL, 5 – TL с трением, 6 – UL с трением).

В результате сравнения решений, полученные в TL, UL, можно заключить, что для данной задачи предпочтительнее использовать TL. Кроме того, на остаточную форму значительное влияние оказывает трение листа с пуансоном.

5. Заключение. В рассмотренном методе оптимизации исходная задача сводится к последовательности вспомогательных более простых задач минимизации. Данный метод позволяет рассматривать различные траектории нагружения и при этом уменьшает объем вычислений по сравнению с простым перебором всевозможных путей деформирования, так как в процессе расчета не оптимальные траектории исключаются. Согласно результатам моделирования формообразования детали одинарной кривизны из тонкого листа и натурных испытаний оптимальной кинематической схемой формообразования по остаточным перемещениям и поврежденности является изгиб с последующим растяжением. Учет в исследуемом множестве траекторий деформирования с участками разгрузки и в задачах обтяжки моделей материалов, включающих эффект Баушингера, не изменяет оптимальную кинематическую схему. Данные траектории будут иметь место в задачах формообразования панелей двойной кривизны [1, 4, 14] либо толстых панелей.

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

Работа выполнена при финансовой поддержке Российского научного фонда (грант № 21-11-00165).

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

  1. Михеев В.А., Сурудин С.В. Основы расчета процесса формообразования обтяжкой тонких оболочек двойной кривизны // Известия Самарского научного центра РАН. 2017. Т. 19. № 1 (3). С. 555–562.

  2. He J., Xia Z.C., Zhu X., Zeng D., Li Sh. Sheet metal forming limits under stretch-bending with anisotropic hardening // Int. J. Mech. Sci. 2013. V. 75. P. 244–256. https://doi.org/10.1016/j.ijmecsci.2013.07.007

  3. Погарцева М.М. Методика создания управляющих программ для обтяжных прессов, применяемых в авиационной промышленности // Сибирский журнал науки и технологий. 2017. Т. 18. № 2. С. 404–414.

  4. Krivenok A.A., Burenin A.A. On the parallel kinematics of the FET stretching press in the stretch forming operations in the manufacture of parts with complex spatial geometry // J. Sib. Fed. Univ. Math. Phys. 2021. V. 14. № 6. P. 1–11. https://doi.org/10.17516/1997-1397-2021-14-6-1-11

  5. Буренин А.А., Ковтанюк Л.В. Большие необратимые деформации и упругое последействие. Владивосток: Дальнаука, 2013. 312 с.

  6. Коробейников C.H. Нелинейное деформирование твердых тел. Новосибирск: Изд-во СО РАН, 2000. 262 с.

  7. Wriggers P. Computational contact mechanics. Heidelberg: Springer, 2006. https://doi.org/10.1007/978-3-540-32609-0

  8. Москвитин В.В. Пластичность при переменных нагружениях. М.: Изд-во МГУ, 1965. 263 с.

  9. Малинин Н.Н. Прикладная теория пластичности и ползучести. М.: Машиностроение, 1975. 400 с.

  10. Каратушин С.И., Храмова Д.А., Пехов В.А. Эффект Баушингера при различных видах пластической деформации // Известия высших учебных заведений. Машиностроение. 2017. № 12 (693). С. 45–50. https://doi.org/10.18698/0536-1044-2017-12-45-50

  11. Nam-Ho Kim. Introduction to Nonlinear Finite Element Analysis. N. Y.: Springer Science+Business Media, 2015. https://doi.org/10.1007/978-1-4419-1746-1

  12. Васильев Ф.П. Методы оптимизации. М.: Факториал Пресс, 2002. 824 с.

  13. Моисеев Н.Н. Элементы теории оптимальных систем. М: Наука, 1975. 528 с.

  14. Бормотин К.С., Герасимов К.Е., Романютин М.И. Численный метод оптимизации кинематической схемы формообразования панелей двойной кривизны // Ученые записки КнАГТУ. 2020. № VII–1 (47). С. 59–69.

Дополнительные материалы отсутствуют.