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

详细

Для численного решения задач оптимального управления в процессах формообразования элементов конструкций в ползучести рассматривается метод динамического программирования. Предлагается реализация разработанного метода в комплексе программ конечно-элементного анализа. Приводится анализ устойчивости данного метода. Библ. 25. Фиг. 7.

全文:

Введение

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

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

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

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

Пусть $V\subset R^3$ — область деформируемого тела с границей S. Поверхность с заданными кинематическими условиями обозначается через S^* $\left(S^{*}\subset S\right).$ Обозначим через $u=(u_1,u_2,u_3)$ вектор перемещений деформируемого тела.

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

$J\left(\dot{u}\right)=a\left(\dot{u},\dot{u}\right)$ при $\dot{u}{|}_{S^{*}}={\dot{u}}^{*},$ (1.1)

где ${\dot{u}}^{*}$ — заданные скорости перемещений в момент времени t; $t\in \left[\mathrm{0,}T\right]$ — время деформирования тела под нагрузкой; потенциальная форма определяется в виде $a(\dot{u},\delta \dot{u})=\underset{V}{\int }[\partial E({\dot{u}}_{i,j}^{})/\partial {\dot{u}}_{i,j}]\delta {\dot{u}}_{i,j}dV$, где $E\left({\dot{u}}_{i,j}^{}\right)=(1/2)c_{ijpl}{\dot{\epsilon }}_{ij}{\dot{\epsilon }}_{pl}-c_{ijpl}{\dot{\epsilon }}_{ij}{\eta }_{pl}^{}+(1/2){\sigma }_{ij}{\dot{u}}_{p,i}^{}{\dot{u}}_{p,j}^{}$, $c_{ijpl}$ — компоненты тензора упругих констант, ${\dot{\epsilon }}_{ij}^{}=(1/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}$; ${\eta }_{pl}^{}$ — компоненты скоростей деформаций ползучести, определяемые по закону установившейся ползучести (см. [13]):

${\eta }_{ij}=\frac{3}{2}B{\sigma }_e^{n-1}{s}_{ij},$ (1.2)

s_ij — компоненты девиатора тензора напряжений, ${\sigma }_e=\sqrt{\frac{3}{2}{s}_{ij}{s}_{ij}}$ — эффективное напряжение (интенсивность напряжений), B, n — константы материала; точкой сверху обозначены скорости перемещений ${\dot{u}}_i$, $i,j,p,l=\mathrm{1,2,3}$.

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

${\dot{\sigma }}_{ij}=c_{ijpl}({\dot{\epsilon }}_{pl}-{\eta }_{pl}^{}).$

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

$A=\underset{V}{\max }\underset{0}{\overset{T}{\int }}{\sigma }_{ij}{\eta }_{ij}dt\to \inf .$ (1.3)

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

В качестве функций управления принимаются перемещения $u\left(t\right)=f\left(t\right)u^{*},$ заданные на границе S^*, а в качестве функции состояния — перемещения, деформации, напряжения в теле V. Таким образом, определив некоторое решение u^* обратной задачи (см. [15]), решается задача поиска оптимальной функции f(t).

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

${\epsilon }_{\alpha \beta }=\frac{1}{2}\left(u_{\alpha ,\beta }+u_{\beta ,\alpha }\right)+\frac{1}{2}{w}_{,\alpha }{w}_{,\beta },$

${\dot{\epsilon }}_{\alpha \beta }=\frac{1}{2}\left({\dot{u}}_{\alpha ,\beta }+{\dot{u}}_{\beta ,\alpha }\right)+\frac{1}{2}{\dot{w}}_{,\alpha }{w}_{,\beta }+\frac{1}{2}{w}_{,\alpha }{\dot{w}}_{,\beta }\left(\alpha ,\beta =\mathrm{1,2}\right).$ (1.4)

Известно, что для пластин в случае малых прогибов оптимальное деформирование происходит по линейному закону $w\left(t\right)=\frac{t}{T}{w}^{*},$ а в случае больших прогибов оптимальное деформирование происходит по нелинейному закону $w\left(t\right)=\sqrt{\frac{t}{T}}{w}^{*},$ где w^* — прогиб в конечный момент времени (см. [6; 12]).

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

Применяя основные процедуры метода конечных элементов к вариационному уравнению с функционалом (1.1), строятся дискретные уравнения задачи деформирования (см. [13; 17])

${}^{t+dt}{K}^{\left(r-1\right)}\Delta U^{\left(r\right)}{=}^{t+dt}{R}^{\left(r-1\right)},$ (2.1)

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

Наряду с дискретизацией по параметру t, вызванной решением нелинейных задач механики методом конечных элементов, для приближенного решения задачи оптимального управления вводится дополнительная сетка: $0<t_1<t_2<\dots <t_N=T$. Учитывая дискретные по времени уравнения пошаговой процедуры интегрирования (2.1) при условии $dt\le \tau =t_{k+1}-t_k$, минимизируемый функционал (1.3) заменяется формулой

$A=\underset{V}{}\sum _{k=0}^{N-1}{\int }_{t_k}^{t_{k+1}}{\sigma }_{ij}{\eta }_{ij}dt=\underset{V}{}\sum _{k=0}^{N-1}\sum _{t=t_k}^{t_{k+1}}{\sigma }_{ij}\Delta {\epsilon }_{ij}^c\to \inf ,$ (2.2)

где $\Delta {{\epsilon }_{ij}}^c$ — компоненты приращений деформаций ползучести, вычисленные методом конечных элементов. В данном случае дискретная задача оптимального управления будет включать дискретные по времени уравнения пошаговой процедуры интегрирования (2.1) и минимизируемый функционал (2.2). В такой постановке строится функция Беллмана и задача решается методом динамического программирования (см. [18; 19]).

Решение задачи оптимизации траекторий деформирования рассматривается на примере формообразования квадратной пластинки с выступами толщиной 12 мм и длиной стороны 180 мм. Выступы необходимы для расчета остаточной формы пластинки с двойной кривизной (см. [20]). Таким образом, находится прогиб пластинки, моделирующий кручение (см. [20]) в виде узловых перемещений по координате, нормальной к поверхности пластинки (значения максимальных перемещений в углах — 80 мм). Данная величина прогиба выбрана с целью уменьшения сопротивления пластинки изгибу и сведения задачи к нахождению прогиба гибкой мембраны (см. [16]). Для более полного анализа рассматривается объемная постановка задачи (фиг. 1) c 8-узловыми изопараметрическими шестигранными конечными элементами с трилинейной аппроксимацией функций. При вычислениях в системе MSC.Marc подключаются дополнительные параметры, улучшающие характеристики сдвига (или изгиба) и описывающие несжимаемое поведение материала путем использования альтернативной функции интерполяции и альтернативной процедуры интегрирования (см. [21]). Количество элементов в модели определено в соответствии с проведенным ранее анализом (см. [20]) и обеспечивает достаточную точность при относительно небольшом времени расчета.

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

В расчетах деформирования пластины используются характеристики материала АК4–1Т (алюминиевого сплава) (см. [20; 22]). В соответствии с этими данными, материал изотропен с параметрами упругости: модуль Юнга E=7000 кГ/мм², коэффициент Пуассона ν=0.4. При температуре T=200 °C стадия установившейся ползучести при сжатии и при растяжении описывается законом Нортона с разными значениями коэффициентов. В данном случае для упрощения вывода условий сходимости метода динамического программирования в задачах оптимального формообразования принимается описание закона ползучести в течение 260 ч с одинаковыми коэффициентами при сжатии и растяжении: $B=2.5\cdot {10}^{-15}{(êã/ì{ì}^2)}^{-n}{(÷àñ)}^{-1}$, n = 8.

При решении аддитивных задач применяется алгоритм, основное содержание которого состоит в формулировке правил последовательного сжатия множества конкурентоспособных вариантов (см. [18; 19]). Алгоритм представляет собой многошаговый процесс, на каждом шаге которого производится исключение некоторого множества вариантов, о котором в процессе работы алгоритма становится известным, что он не содержит оптимального варианта.

Для разработки алгоритма оптимизации при деформировании заготовки в качестве управляющих параметров вводится вектор-функция перемещений узловых точек тела на границе S^* в виде $U_z\left(t\right)=f\left(t\right){U_z}^{*}$, где ${U_z}^{*}$ — решение обратной задачи с линейной функцией f(t) (см. [15]), обеспечивающее необходимую остаточную форму панели. В этом случае строится сетка в пространстве (t, z). Шаг по аргументу t задан и равен $\tau =t_{k+1}-t_k$, по переменной z — Dz. Узлы сетки обозначим через P_g(k). Индекс k означает номер гиперплоскости S_k при заданном значении t, а индекс g означает номер узла в гиперплоскости S_k. Каждые два узла, лежащие в гиперплоскостях $P_q\left(k\right)$ и $P_g(k+1)$, соединены отрезками, длины этих отрезков обозначаются $l_{qg}\left(k\right)=f_k\left(P_q\right(k),P_g(k+1\left)\right)$ (см. [18]).

В результате таких операций можно получить граф, в котором роль вершин играют узлы P_g(k), и вместо исходной задачи будет рассматриваться задача поиска на этом графе кратчайшего пути, соединяющего гиперплоскости S₀ и S_N. Обозначая через l_q(k) ломаную кратчайшей длины, соединяющую узел P_q(k) с гиперплоскостью S₀, можно прийти к рекуррентному соотношению (см. [18])

$l_g(k+1)=\underset{q}{\min }\left\{l_q\left(k\right)+l_{qg}\left(k\right)\right\}$.

Минимум берется по тем номерам q, для которых узлы лежат в допустимой области G_k и принадлежат гиперплоскости S_k.

Для программной реализации метода динамического программирования и построения функции $U_z\left(t\right)=f\left(t\right){U_z}^{*}$ предлагается следующий способ задания граничных условий для рассматриваемой сетки при решении задач (2.1). Шаги метода динамического программирования вычисляются по формуле

$\tau =t_k-t_{k-1}=\frac{T}{N}$, $k=\mathrm{1,...,}N$, $t_0=0$, $\Delta z=\frac{{U_z}^{*}}{M}$.(2.3)

На каждом интервале $[t_{k-1},t_k]$ при решении задачи уравнениями (2.1) задаются граничные условия на перемещения DU_z по следующему алгоритму:

$\Delta t_{step}={\xi }_k\tau$,

$\Delta U_z=({\zeta }_k-{\zeta }_{k-1})\frac{\Delta z}{\tau }dt$, пока $\Delta t_{step}\le t+dt$, иначе $k=k+1$ и повторное выполнение операций.

Граничные условия на всех интервалах $[t_{k-1},t_k]$ могут быть представлены с помощью системы параметров

$[{\xi }_1,{\zeta }_1;{\xi }_2,{\zeta }_2;\mathrm{...};{\xi }_N,{\zeta }_N]$,

где x_k могут принимать значения 0,1,.., N (при условии ${\xi }_k>{\xi }_{k-1}$), а z_k — значения 0,1,.., M.

Вычисления рекуррентных соотношений выполняются путем построения итераций с различными системами параметров и решения уравнений (2.1) в системе MSC.Marc (см. [21]). Ввод граничных условий и вывод значения критерия оптимизации выполняется с помощью разработанных пользовательских программ.

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

Фиг. 2. Траектория деформирования пластинки при N = 3, M = 9.

Фиг. 3. Траектория деформирования пластинки при N = 4, M = 12.

Фиг. 4. Траектория деформирования пластинки при N = 6, M = 18.

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

3. Исследование сходимости численного метода оптимизации

Пусть сетка Q_m характеризуется шагом t_m по временной переменной и шагом $\Delta z_m$ по пространственной переменной. Последовательность сеток такова, что ${\tau }_m\to 0$, $\Delta z_m\to 0$ при $m\to \infty$. Каждой сетке Q_m можно поставить в соответствие конечное множество траекторий $\left\{f\right(t,m\left)\right\},$ построенных с помощью элементарной операции (см. [18; 19]). Данные траектории представляют собой ломаные, которые проходят через узлы, лежащие на гиперплоскостях ${\sum }_k\left(t=k{\tau }_m\right).$ С помощью метода динамического программирования можно определить ломаную, соединяющую начальную и конечную точки и имеющую минимальную длину. Под длиной понимается величина диссипируемой работы (2.2)

$A\left(f\right(t,m\left)\right)=\underset{V}{\max }\sum _{k=0}^{N_m-1}\underset{t_k}{\overset{t_{k+1}}{\int }}{\sigma }_{ij}{\eta }_{ij}dt,$

где $N_m=T/{\tau }_m$, $t_k=k{\tau }_m$.

Обозначим

${\Delta }_m=\left|A\left(f\left(t,m\right)\right)-A^0\right|,$ (3.1)

где диссипируемая работа A⁰ определяется по оптимальной криволинейной траектории f(t), а A — по оптимальной ломаной траектории f(t, m), найденной с помощью метода динамического программирования.

Согласно лемме из [18], если $\underset{m\to \infty }{\lim }{\Delta }_m=\mathrm{0,}$ то траектория деформирования f(t, m) сходится слабо к оптимальной f(t).

При известном потенциале скоростей ползучести можно записать (см. [23])

${\eta }_{ij}=\frac{\partial \Phi }{\partial {\sigma }_{ij}}.$

Тогда, если удельная мощность рассеянной при ползучести энергии $W=W\left({\sigma }_{ij}\right)={\sigma }_{ij}{\eta }_{ij}$ есть однородная функция напряжений степени n+1, то $\Phi \left({\sigma }_{ij}\right)=\frac{1}{n+1}W\left({\sigma }_{ij}\right)$ (см. [23]).

Отсюда следуют и обратные соотношения

${\sigma }_{ij}=\frac{\partial U}{\partial {\eta }_{ij}^{}}.$

В этом случае можно определить $U\left({\eta }_{ij}\right)=\frac{n}{n+1}W\left({\eta }_{ij}\right)$ (см. [24]). Условие устойчивости в малом эквивалентно условию выпуклости функций $W\left({\sigma }_{ij}\right)$ и $W\left({\eta }_{ij}\right)$ (см. [6]).

В силу выпуклости закона выполняется неравенство для любых двух путей деформирования ${\eta }_{ij}$, ${\eta }_{ij}^0$ (см. [19]):

$W^0-W\ge \frac{n+1}{n}{\sigma }_{ij}({\eta }_{ij}^0-{\eta }_{ij})$ или $W-W^0\le \frac{n+1}{n}{\sigma }_{ij}({\eta }_{ij}-{\eta }_{ij}^0)$.

Далее будем обозначать с помощью индекса «0» все величины, относящиеся к оптимальному пути деформирования, тогда $W-W^0\ge 0$ и разница диссипируемой работы примет вид

$A-A^0=\underset{V}{\max }\underset{0}{\overset{T}{\int }}\left(W-W^0\right)dt=\underset{V}{\max }\underset{0}{\overset{T}{\int }}\left({\sigma }_{ij}{\eta }_{ij}-{\sigma }_{ij}^0{\eta }_{ij}^0\right)dt\le \frac{n+1}{n}\underset{V}{\max }\underset{0}{\overset{T}{\int }}{\sigma }_{ij}\left({\eta }_{ij}-{\eta }_{ij}^0\right)dt.$

Таким образом, (3.1) примет вид

${\Delta }_m\le \frac{n+1}{n}\underset{V}{}\left|\sum _{k=0}^{N_m-1}{\int }_{t_k}^{t_{k+1}}{\sigma }_{ij}\left({\eta }_{ij}-{\eta }_{ij}^0\right)dt\right|\le \frac{n+1}{n}\underset{V}{}\sum _{k=0}^{N_m-1}\left|{\int }_{t_k}^{t_{k+1}}{\sigma }_{ij}\left({\eta }_{ij}-{\eta }_{ij}^0\right)dt\right|.$ (3.2)

Считая, что изменения объема на установившейся стадии ползучести не происходит, т.е. ${\eta }_{11}+{\eta }_{22}+{\eta }_{33}=0$ (см. [23]), тензор скорости деформаций ползучести может рассматриваться как девиатор скорости деформаций ползучести. Так как гидростатическое напряжение работы не совершает, то тензор напряжений может быть заменен девиатором тензора. В этом случае мощность рассеяния энергии при ползучести может быть представлена в виде $W=s_{ij}{\eta }_{ij}$.

Представим мощность удельной рассеянной при ползучести энергии в виде (см. [6])

$W={\sigma }_e{\eta }_e$ , (3.3)

где интенсивность скоростей деформаций ползучести ${\eta }_e=\sqrt{\frac{2}{3}{\eta }_{ij}{\eta }_{ij}}$.

Умножим (1.2) на s_ij и свернем: $W=B{\sigma }_e^{n+1}.$ Сравнивая с (3.3), получим ${\eta }_e=B{\sigma }_e^n,$ откуда ${\eta }_{ij}=B{\sigma }_e^n\frac{\partial {\sigma }_e}{\partial s_{ij}}$ совпадает с (1.2). Тогда $s_{ij}{\eta }_{ij}={\sigma }_e{\eta }_e$ продифференцируем по ${\eta }_{kl}:$

$s_{kl}+{\eta }_{ij}\frac{\partial s_{ij}}{\partial {\eta }_{kl}}={\sigma }_e\frac{\partial {\eta }_e}{\partial {\eta }_{kl}}+{\eta }_e\frac{\partial {\sigma }_e}{\partial {\eta }_{kl}}={\sigma }_e\frac{\partial {\eta }_e}{\partial {\eta }_{kl}}+{\eta }_e\frac{\partial {\sigma }_e}{\partial s_{ij}}\frac{\partial s_{ij}}{\partial {\eta }_{kl}}$, но ${\eta }_e\frac{\partial {\sigma }_e}{\partial s_{ij}}={\eta }_{ij}$.

Откуда получим $s_{kl}={\sigma }_e\frac{\partial {\eta }_e}{\partial {\eta }_{kl}}$. Но ${\sigma }_e^{}=B^{-\frac{1}{n}}{{\eta }_e}^{\frac{1}{n}}$, тогда $s_{kl}=B^{-\frac{1}{n}}{{\eta }_e}^{\frac{1}{n}}\frac{\partial {\eta }_e}{\partial {\eta }_{kl}}=\frac{2}{3}{B}^{-\frac{1}{n}}{{\eta }_e}^{\frac{1}{n}-1}{\eta }_{kl}$. Таким образом, напряжения определяются по скоростям деформаций ползучести с помощью соотношения

$s_{ij}=\frac{2}{3}{B}^{-\frac{1}{n}}{{\eta }_e}^{\frac{1}{n}-1}{\eta }_{ij}$ .(3.4)

Для оценки D_m в (3.2), будет использоваться модель гибкой мембраны, для которой скорости деформаций определяются по (1.4). Мембрана толщиной h занимает в плоскости $x_{1}O{x}_2$ область S. Примем такие же допущения, как и в [6]:

1) если время деформирования велико, то компонентами скоростей упругих деформаций можно пренебречь, таким образом, принимаем ${\dot{\epsilon }}_{\alpha \beta }\approx {\eta }_{\alpha \beta }$;

2) при деформировании тонких пластин остаточная форма в основном определяется прогибом w, так как перемещения u_a в плоскости пластин много меньше w.

Пусть искомый прогиб имеет вид

$w\left(t,x_1,x_2\right)=f\left(t\right)w^{*}\left(x_1,x_2\right).$ (3.5)

Тогда в силу указанных предположений и (1.4) ${\eta }_{\alpha \beta }\approx \frac{1}{2}{\dot{w}}_{,\alpha }{w}_{,\beta }+\frac{1}{2}{w}_{,\alpha }{\dot{w}}_{,\beta }=\dot{f}f\left(w_{,\alpha }^{*}{w}_{,\beta }^{*}\right)$. Используя (3.4) с учетом ${\eta }_e=\dot{f}f\sqrt{\frac{2}{3}{w^{*}}_{,\alpha }{w^{*}}_{,\beta }{w^{*}}_{,\alpha }{w^{*}}_{,\beta }}$, найдем напряжения

$\begin{array}{c}{s}_{\gamma \lambda }=\frac{2}{3}{B}^{-\frac{1}{n}}{\left(\dot{f}f\right)}^{\frac{1-n}{n}}\times \\ \times \dot{f}f{\left(\frac{2}{3}{w^{*}}_{,\alpha }{w^{*}}_{,\beta }{w^{*}}_{,\alpha }{w^{*}}_{,\beta }\right)}^{\frac{1-n}{2n}}\left({w^{*}}_{,\gamma }{w^{*}}_{,\lambda }\right)=\\ ={\left(\frac{2}{3}\right)}^{\frac{1+n}{2n}}{B}^{-\frac{1}{n}}{\left(\dot{f}f\right)}^{\frac{1}{n}}\frac{\left({w^{*}}_{,\gamma }{w^{*}}_{,\lambda }\right)}{{\left({w^{*}}_{,\alpha }{w^{*}}_{,\beta }{w^{*}}_{,\alpha }{w^{*}}_{,\beta }\right)}^{\frac{n-1}{2n}}}\left(\alpha ,\beta ,\gamma ,\lambda =\mathrm{1,2}\right).\end{array}$

В этом случае (3.2) примет вид

${\Delta }_m\le \frac{n+1}{n}\underset{S}{\max }h\times \sum _{k=0}^{N_m-1}\left|\underset{t_k}{\overset{t_{k+1}}{\int }}\left(\dot{f}f-{\dot{f}}^0{f}^0\right)s_{\alpha \beta }{w^{*}}_{,\alpha }{w^{*}}_{,\beta }dt\right|=C\sum _{k=0}^{N_m-1}\left|\underset{t_k}{\overset{t_{k+1}}{\int }}\left(\dot{f}f-{\dot{f}}^0{f}^0\right){\left(\dot{f}f\right)}^{\frac{1}{n}}dt\right|,$ (3.6)

где $C={\left(\frac{2}{3}\right)}^{\frac{1+n}{2n}}\frac{n+1}{n}{B}^{-\frac{1}{n}}h\underset{S}{\max }\left|{\left({w^{*}}_{,\alpha }{w^{*}}_{,\beta }{w^{*}}_{,\alpha }{w^{*}}_{,\beta }\right)}^{\frac{n+1}{2n}}\right|$.

Отрезки, соединяющие точки с соседних гиперповерхностей, описываются уравнениями $f\left(t\right)=f\left(t_k\right)+\frac{f\left(t_{k+1}\right)-f\left(t_k\right)}{{\tau }_m}(t-t_k)$, тогда $\dot{f}\left(t\right)=\frac{f\left(t_{k+1}\right)-f\left(t_k\right)}{{\tau }_m}$.

Обозначив значения функций в моменты времени $f\left(t_k\right)=f_k$, можно записать $\dot{f}\left(t\right)=\frac{f_{k+1}-f_k}{{\tau }_m}$ и $\dot{f}f=\frac{f_{k+1}-f_k}{{\tau }_m}{f}_{k+1}$. Положим, что $f_k={f^0}_k+\Delta f_k$, тогда

$\begin{array}{c}\left|\frac{f_{k+1}-f_k}{{\tau }_m}{f}_{k+1}\right|=\left|\frac{{f^0}_{k+1}+\Delta f_{k+1}-{f^0}_k-\Delta f_k}{{\tau }_m}\left({f^0}_{k+1}+\Delta f_{k+1}\right)\right|\le \\ \le \left|\frac{{f^0}_{k+1}-{f^0}_k}{{\tau }_m}\left({f^0}_{k+1}+\Delta f_{k+1}\right)\right|+\left|\frac{\Delta f_{k+1}-\Delta f_k}{{\tau }_m}\left({f^0}_{k+1}+\Delta f_{k+1}\right)\right|\le \\ \le \frac{1}{{\tau }_m}\left|\left({f^0}_{k+1}-{f^0}_k\right){f^0}_{k+1}\right|+\frac{1}{{\tau }_m}\left|\left({f^0}_{k+1}-{f^0}_k\right)\Delta f_{k+1}\right|+\frac{1}{{\tau }_m}\left|\left(\Delta f_{k+1}-\Delta f_k\right){f^0}_{k+1}\right|+\frac{1}{{\tau }_m}\left|\left(\Delta f_{k+1}-\Delta f_k\right)\Delta f_{k+1}\right|\le \\ \le \frac{a}{{\tau }_m}+\frac{b}{{\tau }_m}\Delta z_m+\frac{c}{{\tau }_m}\Delta z_m+\frac{d}{{\tau }_m}\Delta {z_m}^2=\frac{1}{{\tau }_m}\left(a+e\Delta z_m+d\Delta {z_m}^2\right),\end{array}$

где использованы обозначения $|\Delta f_{k+1}-\Delta f_k|\le \Delta z_m$, $\left|\Delta f_k\right|\le \frac{1}{2}\Delta z_m$, $a=\underset{t}{\max }\left|\right({f^0}_{k+1}-{f^0}_k\left){f^0}_{k+1}\right|$, $b=\underset{t}{\max }\frac{1}{2}|{f^0}_{k+1}-{f^0}_k|$, $c=\underset{t}{\max }\left|{f^0}_{k+1}\right|$, $d=\frac{1}{2}$, $e=b+c$.

Аналогично оценим выражение

$\begin{array}{c}\left|\dot{f}f-{\dot{f}}^0{f}^0\right|=\left|\frac{{f^0}_{k+1}+\Delta f_{k+1}-{f^0}_k-\Delta f_k}{{\tau }_m}\times \right.\left.\left({f^0}_{k+1}+\Delta f_{k+1}\right)-\frac{{f^0}_{k+1}-{f^0}_k}{{\tau }_m}{f^0}_{k+1}\right|=\\ =\left|\frac{{f^0}_{k+1}-{f^0}_k}{{\tau }_m}\Delta f_{k+1}+\frac{\Delta f_{k+1}-\Delta f_k}{{\tau }_m}\times \right.\left.\left({f^0}_{k+1}+\Delta f_{k+1}\right)\right|\le \frac{1}{{\tau }_m}\times \left(b\Delta z_m+c\Delta z_m+d\Delta {z_m}^2\right)=\\ =\frac{1}{{\tau }_m}\left(e\Delta z_m+d\Delta {z_m}^2\right).\end{array}$

Тогда для подынтегрального выражения (3.6) будет выполняться неравенство

$\begin{array}{c}\left|(\dot{f}f-{\dot{f}}^0{f}^0){\left(\dot{f}f\right)}^{\frac{1}{n}}\right|\le \frac{1}{{\tau }_m}\left(e\Delta z_m+d\Delta {z_m}^2\right)\times {\left[\frac{1}{{\tau }_m}\left(a+e\Delta z_m+d\Delta {z_m}^2\right)\right]}^{\frac{1}{n}}=\\ =\frac{1}{{\tau }_m}{\mu }_m{\left[\frac{1}{{\tau }_m}\left(a+{\mu }_m\right)\right]}^{\frac{1}{n}}=\frac{1}{{{\tau }_m}^{\frac{n+1}{n}}}{\left[\left(a{{\mu }_m}^n+{{\mu }_m}^{n+1}\right)\right]}^{\frac{1}{n}},\end{array}$

где ${\mu }_m=e\Delta z_m+d\Delta {z_m}^2$.

В результате (3.2) примет вид

${\Delta }_m\le C{N}_m{\tau }_m\frac{1}{{{\tau }_m}^{\frac{n+1}{n}}}{\left[\left(a{{\mu }_m}^n+{{\mu }_m}^{n+1}\right)\right]}^{\frac{1}{n}}=CT\frac{1}{{{\tau }_m}^{\frac{n+1}{n}}}{\left[\left(a{{\mu }_m}^n+{{\mu }_m}^{n+1}\right)\right]}^{\frac{1}{n}}.$ (3.7)

Сравнивая функцию оптимального деформирования, полученного аналитически в случае больших прогибов, с (3.5), найдем $f^0\left(t\right)=\sqrt{\frac{t}{T}}$. Используя разложение в ряд Тейлора, ограничиваясь членами первого порядка, можно получить $f^0\left(t_{k+1}\right)\approx f^0\left(t_k\right)+\frac{1}{2}\sqrt{\frac{{\tau }_m}{T}}\frac{1}{\sqrt{k}}\approx f^0\left(t_k\right)+\frac{1}{2\sqrt{N}}\frac{1}{\sqrt{k}}.$

В этом случае коэффициенты могут быть заданы значениями

$a=\frac{1}{2\sqrt{N}}$, $b=\frac{1}{4\sqrt{N}}$, $c=1$, $e=\frac{1}{4\sqrt{N}}+1$.

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

$a{\left(e\Delta z_m+d\Delta {z_m}^2\right)}^n+{\left(e\Delta z_m+d\Delta {z_m}^2\right)}^{n+1}=\rho {\tau }_m^{n+1+\epsilon },$ (3.8)

где r, e — произвольные положительные постоянные.

Графики изменения величины D_m, вычисленной по формуле (3.7), в зависимости от значений шагов по пространственной и временной переменным, определенных по N, M в (2.3), представлены на фиг. 5. Константа C в (3.7) выбрана при условии совпадения аналитических и численных значений D_m при N=2.

Фиг. 5. Изменение D_m в зависимости от N, M по формуле (3.7).

По результатам численного решения методом динамического программирования задачи оптимального деформирования пластинки (фиг. 1) с разными значениями N, M вычислены значения D_m по (3.1) и $e={\left(\sum _{k=0}^{N-1}\left(f\right(t_k,m)-f^0(t_k){)}^2\right)}^{\frac{1}{2}}$, которые представлены на фиг. 6, 7.

Фиг. 6. Численные результаты изменения D_m в зависимости от N, M.

Фиг. 7. Численные результаты изменения e в зависимости от N, M.

Заключение

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

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