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

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

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

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

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

Поступила в редакцию 21.05.2022
После доработки 22.05.2022
Принята к публикации 23.05.2022

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

Аннотация

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

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

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

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

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

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

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

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

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

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

(2.1)
${{J}_{1}}({\mathbf{\dot {\bar {u}}}},{\mathbf{\dot {u}}}) = {{\dot {W}}_{c}} + a({\mathbf{\dot {u}}},{\mathbf{\dot {u}}})\quad {\text{при}}\quad {\mathbf{\dot {\bar {u}}}}{{|}_{{{{S}_{c}}}}} = {\mathbf{\dot {\bar {u}}}}{\kern 1pt} *$
где ${\mathbf{\dot {\bar {u}}}}{\kern 1pt} *$ – заданные скорости перемещений контактных тел в момент времени t; $t \in [0,T]$ – время деформирования тела под нагрузкой; Wc – соотношения, полученные наложением контактных условий на движения тел методом множителей Лагранжа или методом штрафных функций; потенциальная форма определяется в виде $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$, где $E(\dot {u}_{{i,j}}^{{}}) = \frac{1}{2}{{c}_{{ijpl}}}{{{{\dot {\varepsilon }}}}_{{ij}}}{{{{\dot {\varepsilon }}}}_{{pl}}} - {{c}_{{ijpl}}}{{{{\dot {\varepsilon }}}}_{{ij}}}{{\dot {\varepsilon }}}_{{pl}}^{c} + \frac{1}{2}{{{{\sigma }}}_{{ij}}}\dot {u}_{{p,i}}^{{}}\dot {u}_{{p,j}}^{{}}$, ${{c}_{{ijpl}}}$ – компоненты тензора упругих констант, ${{\dot {\varepsilon }}}_{{ij}}^{{}} = {{(1} \mathord{\left/ {\vphantom {{(1} {2)}}} \right. \kern-0em} {2)}}(\dot {u}_{{i,j}}^{{}} + \dot {u}_{{j,i}}^{{}} + \dot {u}_{{p,i}}^{{}}u_{{p,j}}^{{}} + u_{{p,i}}^{{}}\dot {u}_{{p,j}}^{{}})$ – компоненты скоростей деформаций, ${{u}_{{i,j}}} = \frac{{\partial {{u}_{i}}}}{{\partial {{x}_{j}}}}$; ${{\dot {\varepsilon }}}_{{pl}}^{c}$ – компоненты скоростей деформаций ползучести, определяемые по закону установившейся ползучести; точкой сверху обозначены скорости перемещений ${{\dot {u}}_{i}}$, ${{\dot {\bar {u}}}_{i}}$, $i,j,p,l = 1,2,3$.

Компоненты скорости второго тензора напряжений Пиола–Кирхгофа определяются соотношениями

${{{{\dot {\sigma }}}}_{{ij}}} = {{c}_{{ijpl}}}({{{{\dot {\varepsilon }}}}_{{pl}}} - {{\dot {\varepsilon }}}_{{pl}}^{c})$

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

(2.2)
$A = \int\limits_0^T {\int\limits_V {{{{{\sigma }}}_{{ij}}}{{\dot {\varepsilon }}}_{{ij}}^{c}dVdt} } \to \inf $

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

В качестве функций управления принимаются перемещения ${\mathbf{\bar {u}}}(t) = f(t){\mathbf{\bar {u}}}{\kern 1pt} *$ точек контактных тел на границе Sc, а функции состояния – перемещения, деформации, напряжения в теле V. Таким образом, определив некоторое решение ${\mathbf{\bar {u}}}{\kern 1pt} *$ обратной задачи [19], решается задача поиска оптимальной функции $f(t)$.

В данной работе предлагается анализ следующих моделей определяющих соотношений в ползучести:

1. Неассоциированный закон установившейся ползучести для изотропных сред, имеющих разные характеристики при сжатии и растяжении [20]:

(2.3)
${{{{\dot {\varepsilon }}}}^{c}}_{{ij}} = {{\gamma }}{{s}_{{ij}}},\quad {{\gamma }}({{{{\sigma }}}_{e}},{{\theta }}) \equiv \frac{1}{2}\left[ {{{{{\gamma }}}_{1}} + {{{{\gamma }}}_{2}} + ({{{{\gamma }}}_{2}} - {{{{\gamma }}}_{1}})\sin 3{{\theta }}} \right]$
где
(2.4)
${{{{\gamma }}}_{1}}({{{{\sigma }}}_{e}}) \equiv \frac{3}{2}{{B}^{ + }}{{\sigma }}_{e}^{{{{n}_{ + }} - 1}},\quad {{{{\gamma }}}_{2}}({{{{\sigma }}}_{e}}) \equiv \frac{3}{2}{{B}^{ - }}{{\sigma }}_{e}^{{{{n}_{ - }} - 1}}$
σe – эффективное напряжение, которое определяется через компоненты девиатора тензора напряжений: ${{{{\sigma }}}_{e}} = \sqrt {\frac{3}{2}{{s}_{{ij}}}{{s}_{{ij}}}} $, θ – угол напряженного состояния.

Для материала АК4-1 при температуре T = 200°С стадия установившейся ползучести при сжатии и при растяжении в течение 250 ч, описывается законом Нортона с разными значениями коэффициентов [12, 13]:

– растяжение: ${{B}^{ + }} = 0.5 \times {{10}^{{ - 14}}}$ (кГ/мм2${{)}^{{ - {{n}_{ + }}}}}$ (час)–1, ${{n}_{ + }} = 8$

– сжатие: ${{B}^{ - }} = 0.25 \times {{10}^{{ - 14}}}$ (кГ/мм2${{)}^{{ - {{n}_{ - }}}}}$ (час)–1, ${{n}_{ - }} = 8$

Для данного материала характеристики упругости одинаковы при растяжении и сжатии: модуль Юнга E = 7000 кГ/мм2, коэффициент Пуассона ${{\nu }} = 0.4$.

2. Неассоциированный закон установившейся ползучести для изотропных сред разносопротивляющихся растяжению и сжатию с учетом параметра поврежденности имеет вид (2.3) с функциями:

(2.5)
${{{{\gamma }}}_{1}}({{{{\sigma }}}_{e}}) \equiv \frac{3}{2}\frac{{{{B}^{ + }}{{\sigma }}_{e}^{{{{n}_{ + }} - 1}}}}{{A_{*}^{{{{m}_{ + }}}}{{{(1 - {{\omega }})}}^{{{{m}_{ + }}}}}}},\quad {{{{\gamma }}}_{2}}({{{{\sigma }}}_{e}}) \equiv \frac{3}{2}\frac{{{{B}^{ - }}{{\sigma }}_{e}^{{{{n}_{ - }} - 1}}}}{{A_{*}^{{{{m}_{ - }}}}{{{(1 - {{\omega }})}}^{{{{m}_{ - }}}}}}}$

Параметр поврежденности ${{\omega }}$ определяется из уравнения:

(2.6)
${{\dot {\omega }}}({{{{\sigma }}}_{e}},{{\theta }}) \equiv \frac{1}{2}\left[ {{{{{{\dot {\omega }}}}}_{1}} + {{{{{\dot {\omega }}}}}_{2}} + ({{{{{\dot {\omega }}}}}_{2}} - {{{{{\dot {\omega }}}}}_{1}})\sin 3{{\theta }}} \right]$
где ${{{{\dot {\omega }}}}_{1}} = \frac{{{{B}^{ + }}{{\sigma }}_{e}^{{{{n}_{ + }} + 1}}}}{{A_{*}^{{{{m}_{ + }} + 1}}{{{(1 - {{\omega }})}}^{{{{m}_{ + }}}}}}}$, ${{{{\dot {\omega }}}}_{2}} = \frac{{{{B}^{ - }}{{\sigma }}_{e}^{{{{n}_{ - }} + 1}}}}{{A_{*}^{{{{m}_{ - }} + 1}}{{{(1 - {{\omega }})}}^{{{{m}_{ - }}}}}}}$.

Параметры данного закона ползучести для материала АК4-1 при температуре T = = 200°С [12, 13]:

растяжение: ${{B}^{ + }} = 0.45 \times {{10}^{{ - 14}}}\;{{({{{\text{кГ}}} \mathord{\left/ {\vphantom {{{\text{кГ}}} {{\text{м}}{{{\text{м}}}^{{\text{2}}}}}}} \right. \kern-0em} {{\text{м}}{{{\text{м}}}^{{\text{2}}}}}})}^{{{{m}_{ + }} - {{n}_{ + }}}}}$ (час)–1, ${{n}_{ + }} = 8$, ${{m}_{ + }} = 7$;

сжатие: ${{B}^{ - }} = 0.29 \times {{10}^{{ - 14}}}\;{{({{{\text{кГ}}} \mathord{\left/ {\vphantom {{{\text{кГ}}} {{\text{м}}{{{\text{м}}}^{{\text{2}}}}}}} \right. \kern-0em} {{\text{м}}{{{\text{м}}}^{{\text{2}}}}}})}^{{{{m}_{ - }} - {{n}_{ - }}}}}$ (час)–1, ${{n}_{ - }} = 8$, ${{m}_{ - }} = 2$.

Параметр поврежденности в случае одноосного деформирования есть отношение текущей работы рассеяния A к ее величине на момент разрушения $A_{*}^{{}}$, ($A_{*}^{{}} = 1$ кГ/мм2) т.е. приведенная работа ${{\omega }} = \frac{A}{{A_{*}^{{}}}}$ [12, 13].

3. Неассоциированный закон установившейся ползучести для анизотропных сред разносопротивляющихся растяжению и сжатию имеет вид (2.3) с учетом следующих функций

(2.7)
${{{{\gamma }}}_{1}}({{{{\sigma }}}_{{ij}}}) \equiv T_{1}^{{{{n}_{ + }}}},\quad {{{{\gamma }}}_{2}}({{{{\sigma }}}_{{ij}}}) \equiv T_{2}^{{{{n}_{ - }}}}$
где квадратичные формы компонент тензора напряжений [15]:

(2.8)
$\begin{gathered} {{T}_{1}} = (A_{{11}}^{ + }{{({{{{\sigma }}}_{{22}}} - {{{{\sigma }}}_{{33}}})}^{2}} + A_{{22}}^{ + }{{({{{{\sigma }}}_{{33}}} - {{{{\sigma }}}_{{11}}})}^{2}} + A_{{33}}^{ + }{{({{{{\sigma }}}_{{11}}} - {{{{\sigma }}}_{{22}}})}^{2}} + \\ \, + 2A_{{12}}^{ + }{{{{\sigma }}}_{{12}}}^{2} + 2A_{{23}}^{ + }{{{{\sigma }}}_{{23}}}^{2} + 2A_{{31}}^{ + }{{{{\sigma }}}_{{31}}}^{2}{{)}^{{\frac{1}{2}}}} \\ \end{gathered} $
(2.9)
$\begin{gathered} {{T}_{2}} = (A_{{11}}^{ - }{{({{{{\sigma }}}_{{22}}} - {{{{\sigma }}}_{{33}}})}^{2}} + A_{{22}}^{ - }{{({{{{\sigma }}}_{{33}}} - {{{{\sigma }}}_{{11}}})}^{2}} + A_{{33}}^{ - }{{({{{{\sigma }}}_{{11}}} - {{{{\sigma }}}_{{22}}})}^{2}} + \\ \, + 2A_{{12}}^{ - }{{{{\sigma }}}_{{12}}}^{2} + 2A_{{23}}^{ - }{{{{\sigma }}}_{{23}}}^{2} + 2A_{{31}}^{ - }{{{{\sigma }}}_{{31}}}^{2}{{)}^{{\frac{1}{2}}}} \\ \end{gathered} $
$A_{{11}}^{ + } = \frac{1}{2}\left( {{{{(B_{{22}}^{ + })}}^{{\frac{2}{{{{n}_{ + }}}}}}} + {{{(B_{{33}}^{ + })}}^{{\frac{2}{{{{n}_{ + }}}}}}} - {{{(B_{{11}}^{ + })}}^{{\frac{2}{{{{n}_{ + }}}}}}}} \right),\quad 2A_{{12}}^{ + } = 4{{(B_{{12}}^{ + })}^{{\frac{2}{{{{n}_{ + }}}}}}} - A_{{11}}^{ + } - A_{{22}}^{ + }$
$A_{{11}}^{ - } = \frac{1}{2}\left( {{{{(B_{{22}}^{ - })}}^{{\frac{2}{{{{n}_{ - }}}}}}} + {{{(B_{{33}}^{ + })}}^{{\frac{2}{{{{n}_{ - }}}}}}} - {{{(B_{{11}}^{ + })}}^{{\frac{2}{{{{n}_{ - }}}}}}}} \right),\quad 2A_{{12}}^{ - } = 4{{(B_{{12}}^{ - })}^{{\frac{2}{{{{n}_{ - }}}}}}} - A_{{11}}^{ - } - A_{{22}}^{ - }$

Остальные компоненты $A_{{ij}}^{ + }$, $A_{{ij}}^{ - }$ (i, j = 1, 2, 3) получаются циклической перестановкой индексов.

На основе экспериментов на растяжение и сжатие сплошных цилиндрических образцов, а также на кручение тонкостенных трубчатых образцов, вырезанных в направлении нормали к плите толщиной 42 мм из сплава АК4-1 и в продольном направлении, при температуре T = 200°С в течение 400 ч найдены параметры закона ползучести для трансверсально-изотропного материала [11, 13]:

${{n}_{ + }} = {{n}_{ - }} = n = 12,\quad B_{{23}}^{ + } = B_{{31}}^{ + } = 2.976 \times {{10}^{{ - 35}}}\;{\text{МП}}{{{\text{а}}}^{{ - n}}} \cdot {{{\text{c}}}^{{ - 1}}}$
$\begin{gathered} B_{{11}}^{ + } = B_{{22}}^{ + } = B_{{33}}^{ + } = B_{{12}}^{ + } = 8.935 \times {{10}^{{ - 35}}}\;{\text{МП}}{{{\text{а}}}^{{ - n}}} \cdot {{{\text{c}}}^{{ - 1}}} \\ B_{{23}}^{ - } = B_{{31}}^{ - } = 0.811 \times {{10}^{{ - 35}}}\;{\text{МП}}{{{\text{а}}}^{{ - n}}} \cdot {{{\text{c}}}^{{ - 1}}} \\ \end{gathered} $
$B_{{11}}^{ - } = B_{{22}}^{ - } = B_{{33}}^{ - } = B_{{12}}^{ - } = 1.805 \times {{10}^{{ - 35}}}\;{\text{МП}}{{{\text{а}}}^{{ - n}}} \cdot {{{\text{c}}}^{{ - 1}}}$
или

$\begin{gathered} B_{{23}}^{ + } = B_{{31}}^{ + } = {\text{8}}{\text{.48}} \times {{10}^{{ - 20}}}\;{{({\text{кГ/м}}{{{\text{м}}}^{{\text{2}}}})}^{{ - n}}}\;{{{\text{ч}}}^{{ - 1}}} \\ B_{{11}}^{ + } = B_{{22}}^{ + } = B_{{33}}^{ + } = B_{{12}}^{ + } = {\text{2}}{\text{.545}} \times {{10}^{{ - 19}}}\;{{({\text{кГ/м}}{{{\text{м}}}^{{\text{2}}}})}^{{ - n}}}\;{{{\text{ч}}}^{{ - 1}}} \\ \end{gathered} $
$\begin{gathered} B_{{23}}^{ - } = B_{{31}}^{ - } = {\text{2}}{\text{.31}} \times {{10}^{{ - 20}}}\;{{({\text{кГ/м}}{{{\text{м}}}^{{\text{2}}}})}^{{ - n}}}\;{{{\text{ч}}}^{{ - 1}}} \\ B_{{11}}^{ - } = B_{{22}}^{ - } = B_{{33}}^{ - } = B_{{12}}^{ - } = {\text{5}}{\text{.14}} \times {{10}^{{ - 20}}}\;{{({\text{кГ/м}}{{{\text{м}}}^{{\text{2}}}})}^{{ - n}}}\;{{{\text{ч}}}^{{ - 1}}} \\ \end{gathered} $

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

(2.10)
${{{{\dot {\varepsilon }}}}^{c}}_{{ij}} = \frac{{\partial \phi }}{{\partial {{{{\sigma }}}_{{ij}}}}}$
где $\phi $ – скалярная потенциальная тензорная функция напряжений [13, 14]:

(2.11)
$\begin{gathered} \phi ({{\sigma }_{e}},\theta ) \equiv \frac{1}{2}[{{\phi }_{1}} + {{\phi }_{2}} + ({{\phi }_{2}} - {{\phi }_{1}})\sin 3\theta ] \\ {{\phi }_{1}}(T) = {{T}_{1}}^{{{{n}_{ + }} + 1}}{\text{/}}{{n}_{ + }} + 1,\quad {{\phi }_{2}}(T) = {{T}_{2}}^{{{{n}_{ - }} + 1}}{\text{/}}{{n}_{ - }} + 1 \\ \end{gathered} $

В этом случае в квадратичных формах (2.8), (2.9) коэффициенты находятся по формулам:

$A_{{11}}^{ + } = \frac{1}{2}\left( {{{{(B_{{22}}^{ + })}}^{{\frac{2}{{{{n}_{ + }} + 1}}}}} + {{{(B_{{33}}^{ + })}}^{{\frac{2}{{{{n}_{ + }} + 1}}}}} - {{{(B_{{11}}^{ + })}}^{{\frac{2}{{{{n}_{ + }} + 1}}}}}} \right),\quad 2A_{{12}}^{ + } = 4{{(B_{{12}}^{ + })}^{{\frac{2}{{{{n}_{ + }} + 1}}}}} - A_{{11}}^{ + } - A_{{22}}^{ + }$
$A_{{11}}^{ - } = \frac{1}{2}\left( {{{{(B_{{22}}^{ - })}}^{{\frac{2}{{{{n}_{ - }} + 1}}}}} + {{{(B_{{33}}^{ + })}}^{{\frac{2}{{{{n}_{ - }} + 1}}}}} - {{{(B_{{11}}^{ + })}}^{{\frac{2}{{{{n}_{ - }} + 1}}}}}} \right),\quad 2A_{{12}}^{ - } = 4{{(B_{{12}}^{ - })}^{{\frac{2}{{{{n}_{ - }} + 1}}}}} - A_{{11}}^{ - } - A_{{22}}^{ - }$

Остальные компоненты $A_{{ij}}^{ + }$, $A_{{ij}}^{ - }$ (i, j = 1, 2, 3) получаются циклической перестановкой индексов.

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

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

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

Наряду с дискретизацией по времени t, вызванной решением нелинейных задач механики методом конечных элементов, для приближенного решения задачи оптимального управления вводится дополнительная сетка: $0 < {{t}_{1}} < {{t}_{2}} < \ldots < {{t}_{N}} = T$. Учитывая дискретные по времени уравнения пошаговой процедуры интегрирования (3.1) при условии $dt \leqslant {{t}_{{k + 1}}} - {{t}_{k}}$ минимизируемый функционал (2.2) заменяется формулой

(3.2)
$\bar {A} = \sum\limits_{k = 0}^{N - 1} {\sum\limits_{t = {{t}_{k}}}^{{{t}_{{k + 1}}}} {\int\limits_V {{{\sigma }_{{ij}}}\Delta {{\varepsilon }_{{ij}}}^{c}dV} } } \to \inf $

В данном случае, дискретная задача оптимального управления будет включать дискретные по времени уравнения пошаговой процедуры интегрирования (3.1) и минимизируемый функционал (3.2). В такой постановке строится функция Беллмана и задача решается методом динамического программирования [2225].

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

Для разработки алгоритма оптимизации при деформировании заготовки в реконфигурируемом пуансоне в качестве управляющих параметров вводится вектор-функция перемещений узловых точек контактных тел на границе ${{S}_{c}}$ (стержней) в виде ${{{\mathbf{\bar {U}}}}_{z}}(t) = f(t){\mathbf{\bar {U}}}_{z}^{*}$, где ${\mathbf{\bar {U}}}_{z}^{*}$ – заданные конечные перемещения. В этом случае строится сетка в пространстве (t, z). Шаг по аргументу t задан и равен $\Delta t$, по переменной z${{\Delta }_{z}}$. Узлы сетки обозначим через ${{P}_{m}}(n)$. Индекс n означает номер гиперплоскости Σn при заданном значении t, а индекс m означает номер узла в гиперплоскости Σn. Каждые два узла, лежащие в гиперплоскостях ${{P}_{q}}(n)$ и ${{P}_{m}}(n + 1)$, соединены отрезками, длины этих отрезков обозначаются ${{l}_{{qm}}}(n) = {{f}_{n}}({{P}_{q}}(n),{{P}_{m}}(n + 1))$ [23].

В результате таких операций строится граф с вершинами ${{P}_{m}}(n)$, и вместо исходной задачи будет рассматриваться задача поиска на этом графе кратчайшего пути, соединяющего гиперплоскости Σ0 и ΣN. Обозначая через ${{l}_{m}}(n)$ ломанную кратчайшей длины, соединяющую узел ${{P}_{m}}(n)$ с гиперплоскостью Σ0, можно прийти к рекуррентному соотношению [23]:

${{l}_{s}}(n + 1) = \mathop {\min }\limits_m \{ {{l}_{m}}(n) + {{l}_{{ms}}}(n)\} $

Минимум берется по тем номерам m, для которых узлы лежат в допустимой области ${{G}_{n}}$ и принадлежат гиперплоскости ${{\Sigma }_{n}}$.

Расчет функций ${{{\mathbf{\bar {U}}}}_{z}}(t) = f(t){\mathbf{\bar {U}}}_{z}^{*}$ выполняется методом динамического программирования, реализованным в системе Marc, по заданной сетке. Размеры сетки или шаги метода динамического программирования вычисляются по формуле:

$\Delta t = {{t}_{k}} - {{t}_{{k - 1}}} = \frac{T}{N},\quad k = 1,...,N,\quad {{t}_{0}} = 0,\quad {{\Delta }_{z}} = \frac{{{\mathbf{\bar {U}}}_{z}^{*}}}{M}$

На каждом интервале $[{{t}_{{k - 1}}},{{t}_{k}}]$ при решении задачи уравнениями (3.1) определяются граничные условия на перемещения контактных тел $\Delta {{{\mathbf{U}}}_{z}}$ по заданному алгоритму [25].

В случае деформирования анизотропных материалов законы движения стержней могут быть различны. Задача оптимизации приобретает d-ю размерность, где d – количество независимо движущихся стержней. Увеличение размерности задачи оптимизации приводит к поиску набора функций, обеспечивающих экстремальное значение критерия, т.е. рассматриваются несколько графов с вершинами $P_{{{{m}_{1}}}}^{1}({{n}_{1}}),...,P_{{{{m}_{d}}}}^{d}({{n}_{d}})$, соответствующих нескольким функциям ${\mathbf{\bar {U}}}_{z}^{1}(t) = {{f}^{1}}(t){\mathbf{\bar {U}}}_{z}^{{1*}}$, …, ${\mathbf{\bar {U}}}_{z}^{d}(t) = {{f}^{d}}(t){\mathbf{\bar {U}}}_{z}^{{d*}}$ для разных контактных тел. Решение задачи оптимизации находится определением на этих графах кратчайших путей с помощью рекуррентного соотношения, причем для каждого узла $P_{{{{m}_{{i - 1}}}}}^{{i - 1}}(n)$ на гиперплоскости ${{\Sigma }_{n}}^{{i - 1}}$ перебираются все узлы $P_{{{{m}_{i}}}}^{i}(n)$ на гиперплоскости ${{\Sigma }_{n}}^{i}$.

В частности, при расчете 2-х оптимальных траекторий (размерность задачи оптимизации d = 2) методом динамического программирования для каждого варианта 1-й траектории анализируются все варианты 2-й траектории.

4. Численные результаты моделирования и оптимизации. Тестирование реализованных в программе MSC.Marc моделей определяющих соотношений (2.3)–(2.11) проводилось при решении задач растяжения и сжатия кубического образца в течение T = = 260 ч постоянными нагрузками (табл. 1). Размеры образца: 1 × 1 × 1 мм.

Таблица 1.

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

Модель Деформации и параметр поврежденности Результаты при растяжении: ${{\sigma }_{x}} = 15$ кГ/мм2 Результаты при сжатии: ${{\sigma }_{x}} = - 15$ кГ/мм2
Модель разносопротивляемости в ползучести (2.3), (2.4) $\varepsilon _{x}^{c}$ 0.0032 –0.0017
$\bar {A}$, кГ/мм2 0.0483 0.0257
Модель разносопротивляемости в ползучести с учетом поврежденности (2.3), (2.5), (2.6) $\varepsilon _{x}^{c}$ 0.0035 –0.002
$\omega $ 0.0519 0.0307
$\bar {A}$, кГ/мм2 0.0519 0.0307
Ассоциированный и неассоциированный закон ползучести для анизотропных разносопротивляющихся сред (2.3), (2.7)–(2.11) $\varepsilon _{x}^{c}$ 0.0079 –0.0018
$\bar {A}$, кГ/мм2 0.1189 0.0271

Согласно численным данным решения задач параметры поврежденности, учитываемые в модели определяющих соотношений (2.6) и вычисленные по результатам (3.2), совпадают. Кроме того, для приближенной оценки параметра поврежденности достаточно использовать модель определяющих соотношений (2.3), (2.4). Завышенные результаты, полученные при ползучести для анизотропных разносопротивляющихся сред, можно объяснить аппроксимацией экспериментальных данных, полученных при большем времени деформирования.

Для оценки особенностей деформирования материала с разносопротивляемостью при растяжении и сжатии рассматривается задача кручения квадратной пластинки [20]. Квадратная пластина имеет размеры 200 мм × 200 мм и толщину 20 мм. Начало прямоугольной системы координат находится в центре пластины, а ось z направлена перпендикулярно пластине (пространственные координаты пластины: 100 мм ≤ х ≤ ≤ 100 мм, 100 мм ≤ у ≤ 100 мм, и 10 мм ≤ z ≤ 10 мм). В каждом углу пластинки для задания распределенной нагрузки установлены квадратные площадки размером 20 мм × × 20 мм × 10 мм. К четырем площадкам в течение 260 ч приложены два набора распределённых сил одинаковой величины и противоположного направления, эквивалентные сосредоточенной силе по величине 1850 кГ. В результате деформирования образуется седловидная поверхность пластины (рис. 1), и данное испытание имеет сложное напряженное состояние, в котором одновременно существуют растягивающие и сжимающие напряжения. Полученные величины прогибов (перемещения по оси z) для разных моделей определяющих соотношений представлены в табл. 2.

Рис. 1.

Конечно-элементная модель исходной и деформированной панели с распределением перемещений по оси z.

Таблица 2.

Значения прогибов в момент времени T = 260 ч, полученные при использовании двух типов обобщенного закона Нортона

Тип обобщенного закона Нортона Прогиб (мм) в узле с координатами (35.35, 35.35, 0) Прогиб (мм) в узле с координатами (90,90,0)
Неассоциированный (2.3), (2.7)–(2.9) 2.774 18.303
Ассоциированный (2.10), (2.11) 2.752 18.086

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

Решение задачи оптимизации траекторий деформирования рассматривается на примере многоточечного формообразования квадратной пластинки в установке с верхней и нижней матрицами, включающими по четыре стержня. Размеры пластинки 500 × 500 × 45 мм. Перемещение жестких тел задается и выполняется на 30 мм. В случае, если имеется только конечная форма панели, то расчет конечного положения стержней, обеспечивающих необходимую упреждающую форму панели, может быть выполнен итерационным методом решения обратной задачи [19].

Анализ решений данной задачи с разными определяющими соотношениями позволяет сделать следующие выводы:

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

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

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

1. закон установившейся ползучести для изотропных сред с одинаковыми характеристиками при сжатии и растяжении (2.3), (2.4): ${{B}^{ + }} = {{B}^{ - }}$;

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

3. неассоциированный закон установившейся ползучести для анизотропных сред разносопротивляющихся растяжению и сжатию (2.3), (2.7), (2.8), (2.9).

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

Рис. 2.

Деформированная конфигурация пластинки и максимальное значение энергии рассеяния.

Известно, что для пластин в случае малых прогибов оптимальное деформирование происходит по линейному закону ${{{\mathbf{\bar {U}}}}_{z}}(t) = \frac{t}{T}{\mathbf{\bar {U}}}_{z}^{*}$, а в случае больших прогибов оптимальное деформирование происходит по нелинейному закону ${{{\mathbf{\bar {U}}}}_{z}}(t) = \sqrt {\frac{t}{T}} {\mathbf{\bar {U}}}_{z}^{*}$ [8–10].

В случае одинакового движения всех контактных тел (d = 1) оптимальная функция представлена на рис. 3, 4 (жирная сплошная кривая – численные результаты, штрихпунктирная кривая – аналитические данные для изгиба пластин [810]). Численное решение задачи оптимизации траектории деформирования, одинаковой для всех стержней, сводится к перебору вариантов при каждом параметре tk. Набор функций f(t) задается ломанными линиями, проходящими от точки O к точке B (рис. 3, 4). В результате, оптимальное решение, полученное методом динамического программирования при N = 3, M = 9 и N = 4, M = 12 для трех случаев определяющих соотношений, приближается к аналитической кривой и не совпадает с линейной функцией.

Рис. 3.

Варианты законов движения контактных тел при N = 3, M = 9 и оптимальная функция.

Рис. 4.

Варианты законов движения контактных тел при N = 4, M = 12 и оптимальная функция.

В случае различного движения контактных тел, в частности, при d = 2 группы контактных жестких тел 1 и 2 (рис. 2) будут перемещаться по разным функциям. Для модели материала, соответствующей закону Нортона (т.е. без учета разносопротивляемости, анизотропии), вычисленная оптимальная система траекторий состоит из двух одинаковых функций, которые совпадают с траекторией на рис. 3, 4 (жирная сплошная линия). Вычисленная оптимальная система траекторий с законом установившейся ползучести для изотропных сред, имеющих разные характеристики при сжатии и растяжении (2.3), (2.4), и с законом установившейся ползучести для анизотропных сред, разносопротивляющихся растяжению и сжатию (2.3), (2.7), (2.8), (2.9), состоит из двух разных функций для групп контактных тел 1 и 2 (рис. 5, 6 соответственно).

Рис. 5.

Варианты законов движения контактных тел при N = 3, M = 9 и оптимальные функции для модели материала (2.3), (2.4).

Рис. 6.

Варианты законов движения контактных тел при N = 3, M = 9 и оптимальные функции для модели материала (2.7)–(2.9).

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

Однако использование данного численного метода требует значительных вычислительных ресурсов. Так, если обозначить через r – время расчета пластинки для произвольного пути деформирования, соединяющего соседние гиперплоскости, то можно оценить время расчета методом динамического программирования и простым перебором всевозможных путей деформирования (табл. 3). При увеличении N, M время расчета сокращается значительно, например, при N = 4, M = 9, d = 1 время расчета методом динамического программирования 10 ⋅ 10 ⋅ 2r + 10 ⋅ 10 ⋅ 3r + 10 ⋅ 4r = 540r, а простым перебором траекторий 10 ⋅ 10 ⋅ 10 ⋅ 4r = 4000r.

Таблица 3.

Сравнение времени расчета методом динамического программирования при N = 3, M = 9 с перебором траекторий

Размерность, d Время расчета методом динамического программирования Время расчета простым перебором траекторий
1 10*10*2r + 10*3r = 230r 10*10*3r = 300r
2 10*10*10*10*2r +10*10*3r =20300 r 10*10*10*10*3r =30000r
4 (10*10)4 *2r + (10)4*3r =20003*104r (10*10)4*3r = 30000*104 r

Вычисления рекуррентных соотношений выполняются путем построения итераций и решения уравнений (3.1) в системе MSC.Marc. Ввод граничных условий и вывод значения критерия оптимизации выполняется с помощью пользовательских программ. Программная реализация метода оптимизации выполнена с учетом алгоритма распределенных вычислений в Visual C, в которой удаленная связь между машинами и потоками устанавливается с помощью протокола распределенной модели COM (Distributed COM). Программа позволяет распределять вычисления на потоки и вычислительные машины.

Результаты параллельного многопоточного конечно-элементного расчета на многопроцессорной вычислительной машине (Core i5-10400F), в зависимости от заданного количества потоков, представлены в табл. 4. Фактические значения времени расчета свидетельствуют о перспективности увеличения вычислительных потоков.

Таблица 4.

Сравнение времени расчета методом динамического программирования при N = 2, M = 9, d = 2 в зависимости от количества потоков на сервере

Время расчета методом динамического программирования Время расчета при заданном количестве потоков, мин
1 2 10 18 20
10*10*2r = 200 r (или 100 задач) 1644 877 291 217 197

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

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

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

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

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

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

  1. Аннин Б.Д., Олейников А.И., Бормотин К.С. Моделирование процессов формообразования панелей крыла самолета SSJ-100 // ПМТФ. 2010. Т. 51. № 4. С. 155–165. https://doi.org/10.1007/s10808-010-0074-2

  2. Ribeiro F.C., Marinho E.P., Inforzato D.J., Costa P.R., Batalha G.F. Creep age forming: a short review of fundaments and applications // J. Achiev. Mater. Manuf. Eng. 2010. V. 43. № 1. P. 353–361.

  3. Adachi T., Kimura S., Nagayama T., Takehisa H., Shimanuki M. Age forming technology for aircraft wing skin // Proc. of the 9th Int. Conf. on Aluminium Alloys. Ed. by J.F. Nie, A.J. Morton and B.C. Muddle. Melbourne: Inst. of Mater. Eng. Australasia Ltd., 2004. P. 202–207.

  4. Кривенок А.А. Моделирование в системе MSC.Marc процесса формообразования деталей в режиме термодеформационного старения с учетом усадки материала // Уч. зап. КнАГТУ. 2013. № III-1(15). С. 4–10. https://doi.org/10.17084/2013.III-1(15).1

  5. Simon D., Kern L., Wagner J., Reinhart G. A reconfigurable tooling system for producing plastic shields // Proc. CIRP. 2014. V. 17. P. 853–858. https://doi.org/10.1016/j.procir.2014.01.095

  6. Walczyk D.F., Lakshmikanthan J., Kirk D.R. development of a reconfigurable tool for forming aircraft body panels // J. Manuf. Sys. 1998. V. 17. № 4. P. 287–296. https://doi.org/10.1016/S0278-6125(98)80076-9

  7. Su S.Z., Li M.Z., Liu C.G., Ji C.Q., Setchi R., Larkiola J., Panteleev I., Stead I., Lopez R. Flexible tooling system using reconfigurable multi-point thermoforming technology for manufacturing freeform panels // Key Eng. Mater. 2012. V. 504–506. P. 839–844. https://doi.org/10.4028/www.scientific.net/KEM.504-506.839

  8. Цвелодуб И.Ю. Об оптимальных путях деформирования в условиях ползучести. Некоторые приложения к задачам обработки материалов давлением // Изв. АН СССР. МТТ. 1987. № 6. С. 128–136.

  9. Цвелодуб И.Ю. Постулат устойчивости и его приложения в теории ползучести металлических материалов. Новосибирск: Ин-т гидродинамики СО АН СССР, 1991. 133 с.

  10. Бормотин К.С., Олейников А.И. Вариационные принципы и оптимальные решения обратных задач изгиба пластин при ползучести // ПМТФ. 2012. Т. 53. № 5. С. 136–146. https://doi.org/10.1134/S0021894412050148

  11. Горев Б.В., Масанов И.Ж. Особенности деформирования листовых конструкционных алюминиевых сплавов и плит в режимах ползучести // Технол. машиностр. 2009. № 7. С. 13–20.

  12. Соснин О.В., Горев Б.В., Рубанов В.В. Кручение квадратной пластинки из материала, разносопротивляющегося растяжению и сжатию при ползучести // Расчеты прочности судовых конструкций и механизмов. Сборник трудов НИИВТа № 117. Новосибирск: Новосиб. ин-т инженеров вод. трансп., 1976. С. 78–88.

  13. Банщикова И.А. Ползучесть изотропных и ортотропных сплавов и длительная прочность элементов конструкций. Дисс. … доктора физ.-мат. наук. Новосибирск, 2020. 338 с.

  14. Банщикова И.А. Построение определяющих уравнений для ортотропных при ползучести материалов с различными свойствами при растяжении и сжатии // ПМТФ. 2020. № 1. С. 102–117. https://doi.org/10.15372/PMTF20200110

  15. Соснин О.В. Об анизотропной ползучести материалов // ПМТФ. 1965. № 6. С. 99–104.

  16. Yuhao Guo, Gang Liu, Yi Huang. A complemented multiaxial creep constitutive model for materials with different properties in tension and compression // Eur. J. Mech. A/Solids. 2022. V. 93. P. 104510. https://doi.org/10.1016/j.euromechsol.2022.104510

  17. Yan Yu, Wan Min, Wang Haibo, Huang Lin. Design and optimization of press bend forming path for producing aircraft integral panels with compound curvatures // Chinese J. Aeronaut. 2010. V. 23(2). P. 274−282. https://doi.org/10.1016/S1000-9361(09)60216-8

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

  19. Бормотин К.С. Итеративный метод решения геометрически нелинейных обратных задач формообразования элементов конструкций в режиме ползучести // Ж. выч. мат. мат. физ. 2013. Т. 53. № 12. С. 145–153. https://doi.org/10.1134/S0965542513120026

  20. Коробейников С.Н., Олейников А.И., Горев Б.В., Бормотин К.С. Математическое моделирование процессов ползучести металлических изделий из материалов, имеющих разные свойства при растяжении и сжатии // Выч. мет. программ. 2008. Т. 9. № 1. С. 346–365.

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

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

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

  24. Бормотин К.С., Вин Аунг. Метод динамического программирования в задачах оптимального деформирования панели в режиме ползучести // Выч. мет. программ. 2018. Т. 19. С. 470–478. https://doi.org/10.26089/NumMet.v19r442

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

  26. Marc 2021, Vol A: Theory and User Information, MSC.Software Corporation URL: http://www.mscsoftware.com/product/marc.

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