Вариационные подходы к построению оптимального управления движением упругого тела Текст научной статьи по специальности «Физика»

____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 157, кн. 3 Физико-математические науки

2015

УДК 517.977.58:539.3

ВАРИАЦИОННЫЕ ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ ДВИЖЕНИЕМ УПРУГОГО ТЕЛА

В.В. Саурин, Г.В. Костин

Аннотация

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

Ключевые слова: оптимальное управление, динамика, линейная теория упругости, вариационные принципы, семи-диекретизация.

Введение

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

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

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

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

122

ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

123

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

Обсуждаемый в настоящей работе подход, который имеет ряд преимуществ и учитывает вышеупомянутые недостатки, присущие вариационным и проекционным методом, а также технике метода наименьших квадратов, получил название метода интегро-дифференциальных соотношений (МИДС) [1]. Суть этого подхода заключается в том, что часть управляющих уравнений выполняется точно, а другие соотношения учитываются в интегральном виде. Соотношения, которые должны быть ослаблены, определяются априори, часто на основе физических соображений. Приближенное решение интегро-дифференциальной задачи находится путем минимизации квадратичного функционала при дифференциальных ограничениях в виде уравнений равновесия, кинематических соотношений и граничных условий. Такая формулировка полностью согласуется с идеями метода наименьших квадратов, но оказывается, что она одновременно является и вариационным принципом.

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

В соответствии с идеями МИДС был разработан проекционный подход как модификация метода Петрова -Галеркина. При использовании полудискретных полиномиальных аппроксимаций и проекционной техники трехмерные статические и динамические задачи теории упругости могут быть решены с высокой точностью.

Подходы, обсуждаемые в настоящей работе на примере оптимизации движений упругого тела, были применены не только к статическим и спектральным задачам [2] теории упругости, но и к прямым и обратным начально-краевым задачам механики деформируемого твердого тела [3], гидро- и термодинамики [4, 5].

1. Вариационная формулировка задачи о движении упругого тела

Рассмотрим упругое тело, занимающее в трехмерном пространстве некоторую ограниченную область V с кусочно-гладкой границей Г. Введем динамические переменные a(t, x), p(t, x) и кинематические переменные w(t, x), e(t, x), характеризующие поведение упругой системы и зависящие от времени t € [0, T], а также вектора x = (xi, Х2, хз) € V пространственных координат. Здесь p и w - векторные функции соответственно плотности импульса и перемещений, а а и е - тензоры второго ранга, определяющие распределение упругих напряжений и деформаций во времени и пространстве. Определим также пространственно-временную область Q = (0,T) х V.

В линейной теории упругости локальные уравнения состояния среды, связывающие скорости точек системы wt с функцией плотности импульса p, а также деформации е с напряжениями а можно записать в виде [6]

v(t, x) = 0, i(t, x)= 0, (t, x) € Q. (1)

Здесь введены соответственно векторная функция невязки по скоростям

v = Wt - p—1 (x)p

(2)

124

В.В. САУРИН, Г.В. КОСТИН

и тензорная функция невязки по деформациям

€ = £ - C-1(x) : а.

(3)

Тензор деформаций £ линейно зависит от вектора перемещений:

£

1

2

(Vw + VwT).

(4)

Объемная плотность тела р и тензор модулей упругости C - заданные функции координат x. Компоненты тензора четвертого ранга C обладают следующими свойствами симметрии: Cjki = Cjik = Cknj. Двоеточие между тензорами означает их скалярное произведение (двойную свертку по индексам). В уравнении (4) используется оператор градиента V = (d/dxi, д/дх2, д/дхз) в пространстве координат x.

Используя соотношения (1) и (4), закон изменения плотности импульса упругого тела можно выразить через тензор напряжений а и вектор плотности импульса p в виде

Считается, что внешние объемные силы отсутствуют. Точка между векторами и тензорами означает свертку по одному индексу.

Рассмотрим случай граничных условий, которые могут быть представлены через компоненты векторов перемещения w и внешней нагрузки q = а ■ n в форме

Здесь n - единичный вектор внешней нормали к границе тела Г, wo и q0 - заданные граничные векторы перемещений и напряжений, Г1 и Г2 - непересекающиеся части границы, то есть Г1 р|Г2 = 0 и TiUГ2 = Г.

Для полного описания движения упругого тела необходимо определить его состояние в начальный момент времени (без ограничения общности при t = 0). Это можно сделать, задавая начальные распределения упругих перемещений w0 и плотности импульса р0 как достаточно гладкие функции координат x:

В предлагаемом подходе локальные уравнения состояния упругого тела (1) заменяются интегральным соотношением, связывающим векторы плотности импульса р и скорости wt, а также тензоры напряжений а и деформаций £.

Была предложена [3] следующая интегро-дифференциальная формулировка начально-краевой задачи о движении упругого тела (1)-(7): найти такие неизвестные поля перемещений w*(t,x), напряжений а* (t,x) и плотности импульса p*(t, х), которые удовлетворяют следующему интегральному соотношению:

при выполнении кинематического уравнения (4), уравнения изменения импульса (5), граничных и начальных условий (6), (7).

Подынтегральная функция р в (8) имеет размерность плотности энергии и неотрицательна. Из последнего свойства следует неотрицательность интеграла Ф для произвольных функций w , а и p , что позволяет привести интегро-дифферен-циальную задачу (4)-(7), (8) к минимизационной формулировке: найти допустимые

Pt(t, x) = V ■ a(t, x), (t, x) e П.

(5)

w(t, x) = wo(t, x), x e Гi;

q(t, x) = qo(t, x), x e Г2.

(6)

w(0, x) = w0(x), p(0, x) = p0(x), x e V.

(7)

V

1

2

(p(x)v ■ v + € : C(x) : €), (8)

ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

125

функции w* , а* и р* , которые доставляют минимальное (нулевое) значение функционалу

Ф[ш*,а*, р*\ = min Ф(м, а,р] = 0 (9)

w,a,p

при выполнении ограничений (4)—(7).

Обозначим действительные и произвольно выбранные допустимые перемещения, напряжения, импульсы через соответственно w*, а*, р* и w, а, р и положим

w = w* + Sw, а = а* + 5а, р = р* + 5р.

Тогда с учетом (9) имеем

Ф(м, а,р] = Ф(м*, а*,р*\ + SwФ + 5аФ + 5рФ + 52Ф = Ф|5м, 5а, 5р],

где SwФ, 5аФ, 5рФ - первые вариации относительно неизвестных w, а, р, а 52Ф > 0 - вторая вариация функционала Ф.

Используя функционал Ф, определенный в (8), для произвольных допустимых полей перемещений w, напряжений а и плотности импульса р, удовлетворяющих ограничениям (5)-(7), можно предложить ряд критериев близости к точному решению. Интегральное качество приближенных функций w, а и р можно оценить по величине безразмерного соотношения

Д = ФФ-1 < 5 < 1. (10)

Здесь 5 - выбранное положительное число, а интеграл по времени от полной механической энергии Ф дается формулой

^ = 2^ (р-1(х)р ■ р + £ : C(x) : е) da (11)

п

Распределение в пространстве и времени точности некоторого допустимого движения w, а и р характеризует функция р^,а, р), определенная в (8).

2. Проекционная формулировка задачи о движении упругого тела

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

Ниже описывается вариант метода Петрова -Галеркина [2], в котором используются интегральные проекции вектора невязки скорости v и тензора невязки деформаций £, введенные в (1). В этом случае постановка задачи линейной упругости заключается в следующем: найти допустимые поля перемещений w, напряжений а и плотности импульсов р, удовлетворяющие закону изменения импульса (5), граничным и начальным условиям (6) и (7) и такие, что выполняется следующие интегральные соотношения: У

У p(x)vt(t, x) ■ r(t, x) d^ = 0 Vr e L2(Q), (12)

п

126

В.В. САУРИН, Г.В. КОСТИН

Рис. 1. Область V, занимаемая упругим телом

J £(t, x) : т(t, x) dO = 0 Vт G L2(О). (13)

Q

Здесь r - вектор виртуальных перемещений, т - тензор виртуальных напряжений, а вектор v и тензор £ определены соответственно в (2) и (3).

3. Задача оптимального управления распределенной системой

Рассмотрим упругое тело (балку), имеющее форму прямоугольного параллелепипеда длиной 2щ и размерами поперечного сечения 2a2 х 2a3 , при этом ai ^ а2 + + а3 (см. рис. 1). Введем декартову систему координат Oxix2x3, начало которой расположено в середине тела, а оси Oxk параллельны ребрам длиной 2ak, к = = 1, 2, 3. Пространственная область задачи определяется как

V = {x : \xi\ < ai , i = 1, 2,3}.

Рассматривается случай, когда длинные стороны балки свободны от нагрузок:

a(t, x) • en = 0, xn = ±an, n = 2, 3. (14)

Предполагаем, что одно из торцевых поперечных сечений свободно от нагрузок

a(t, x) • ei = 0, xi = ai, (15)

а другое сечение остается недеформированным и перемещается в соответствии с заданным законом управления u(t):

w(t, x) = (0,yi (t), 0), yi(t) = u(t), xi = -ai. (16)

Здесь ek = (Sik, &2k, $3k), к = 1, 2,3, - орты системы координат Oxix2x3, нормальные к различным участкам границы тела V, j - символ Кронекера. Перемещение концевого сечения yi (t) вдоль оси Ox2 удовлетворяет следующим начальным условиям:

yi (0) = yi (0) = 0. (17)

Ставится задача нахождения ускорения u* (t) - оптимального управления [8], которое переводит балку из начального состояния покоя

w(0, x) = 0, p(0, x) = 0, x G V, (18)

в терминальное

w(T, x) = (0,yx, 0), p(T, x) = 0, x G V, (19)

и минимизирует функционал качества

J [u* ] = min J [u]. neb2 (q,t )

(20)

ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

127

Рассматривается следующий квадратичный функционал

т

J =1У u2(t)dt + 7^, y > 0, (21)

о

где y - заданный весовой коэффициент, Ф - интеграл по времени от энергии, определяемый (11).

4. Алгоритм дискретизации

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

4.1. Семи-дискретизация по пространственным переменным. Для

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

t

p(t, x) = j V- а( t, x) dt. (22)

о

Подставляя выражение (22) в выражение для вектора v, введенного в (2), получим

t

v = wt — Р-1 J V ■ а( t, x) dt. (23)

о

Тензор невязки по деформации, определенный в (3), после подстановки тензора Коши (4) примет вид

£ =2 (Vw + VwT) — C-1 : а. (24)

Зададим аппроксимации неизвестных полей перемещений и напряжений в виде

N N-1

wi(t, x)= J2 w^1 (t,xi)x\ xl3, ws(t, x)= J2 (t,xi)x\ x3,

k+l = 0 k+l=0

N-1

W2(t, x) = yi(t) + w{2kl) (t, xi)xk xt l3,

k+l=0

N N

(t, x) = = ^2 a1kill(t,xiX ann(t, x) = gn J2 a(kl) nn (t,xi)xk x3, (25)

k+l = 0 k+l = 0

N-1 N-2

(t, x) = 9n ain\t,xi)xk 4 а23(t x) = д2дз E a23l)(t,x1)x2 ^

k+l=0 k+l=0

128

В.В. САУРИН, Г.В. КОСТИН

gn

1 — x,

2

n = 2, 3.

1

xn = a

x

n

Здесь N - заданное целое положительное число, которое определяет степени полиномиального разложения неизвестных функций по безразмерным координатам Ж2 и Хз. Выбранные таким образом приближения автоматически удовлетворяют однородным граничным условиям в напряжениях (14) на боковых поверхностях призмы.

Аппроксимации (25) с использованием проекционных соотношений (12), (13) позволяют построить систему дифференциально-алгебраических уравнений в частных производных относительно времени t и координаты xi.

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

j vt(t, x) ■ r(t, x)d^ = 0 Vr,

(26)

Je, . (X x) ■ Х)<Щ = ° Vs.

(27)

Здесь введены векторы виртуальных импульсов r и напряжений s = т ■ e, в виде

N

N

ri(t, x)= £ r1kl)(t,xi)xk xl3, si(t, x)= £ s(id\t,xi)

™k xl x2 x3,

k+l=0

N1

k+l=0

N1

(28)

rn(t, x)= £ r(kV>(t,xi)xk; xl3, Sn(t, x)= £ Snl) (t,x, )xk; x3, n = 2, 3.

k+l=0

k+l = 0

Уравнения (26) и (27) можно в явном виде разрешить относительно первых производных функций dw^/dx,, 3a(km>/дх,, m = 1, 2, 3. Общее число этих функций согласно (25), (28) равно числу виртуальных функций rmf>, Sm1 и составляет 2Nd, где Nd = (N + 1)(3N + 2)/2. Если подставить выражения для этих производных по координате xi в функционал Ф из (8), то алгебраические соотношения, которые необходимы для нахождения функций напряжений a(2l>, a33l>, a(3l> (общим

числом Na вариации

3 5

2 N2 + 2 N + 2 ), можно получить из условия равенства нулю первой

0^22 Ф + ^23 Ф + 0^3 3 Ф = 0

(29)

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

(kl> (kl> (kl> (kl> (kl>

стовых функций rm , sm и вариаций оа22 , °а33 , ^а2з эквивалентна системе

2Nd + Na линейных уравнений относительно переменных w^kl>, а^пи . Дифференциальный порядок системы как по времени t, так и по координате xi равен 2Nd. Из граничного условия (15) следуют Nd краевых условий в напряжениях:

а(У>(t, ai) = 0, i + j < N; ai^*(t,ai)= ai34(t,ai)=0, k + l < N — 1.

(kl>

(kl>

(30)

Из условия (16) и вида аппроксимаций (25) вытекают Nd однородных условий в перемещениях:

w(i>(t, —ai) = 0, i + j < N ; w2kl>(t, —ai)

w(kl>(t, —ai)=0, k +1 < N — 1. (31)

ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

129

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

xi) = 0, i + j < N; w^1 >(0,xi ) = w^1 'l(0,xi) = 0, k + l < N — 1;

„{kl)

(hi)f

,,lij^(0, xi) = 0, i + j < N; w2,^(0, xi) = w3,3(0,xi) = 0, k + l < N — 1.

= „ !kl)(

(32)

Для получения совместной системы к соотношениям (26), (27), (29)—(32) необходимо добавить дифференциальное уравнение относительно функции перемещения сечения балки yi(t) из (16) вместе с начальными условиями (17).

4.2. Краевая задача на собственные значения. Чтобы решить сформулированную в частных производных конечномерную задачу, представим искомые функции (25) в следующем виде [7]

м -1 м -1

wi(t, x) = ^ Wi,i(x)yi+1(t), w3(t, x) = ^2 w3ii(x)yi+1(t),

i=1 i=1 . .

(33)

м-1 м-1

W2(t, x) = yi(t)+£ W2i(x)yi+i(t), a = ^2 °i(x)yi+i(t),

i=1 i=1

где M - заданное положительное число.

Подстановкой (33) и заменой yj (t) = exp(iwt), j = 1,..., M, задача в частных производных (26), (27), (29)-(32) сводится к системе обыкновенных дифференциальных уравнений относительно x1 с однородными граничными условиями, заданными на торцах балки (задача нахождения собственных частот ш).

Решение задачи на собственные значения позволяет с учетом обозначений введенных в (25) представить решение в следующем виде [2]:

N

W1 = ^2 w[kl\xi)xk xl

k+l=0

N

N-

x2 ^ wn = ^2 Wnll(x 1 )X2 ^

k+l=0

N

x2 x3, xnn ^ аПП (xi)x2 x3,

k+l=0

x i i = Yla iill(xi)xk xl

k+l=0

N- N-2

x1n = 9n d'lknl(x 1)xk ^ x23 = 9293 ^2 x(3l)(x 1)xk X3.

k+l=0

k+l=0

(34)

Здесь njlkl') (xi), alkl\xi) - компоненты соответствующих собственных векторов.

В дальнейшем исследуется частный случай колебаний однородной изотропной балки с квадратным поперечным сечением (о,2 = 0,3). Вводя характерные длину x = 02 и время t = 02 \/р/Е, где Е - модуль Юнга, можно привести все рассмотренные уравнения линейной упругости к безразмерному виду. В системе остается два параметра: относительная длина балки a = 01/02 и коэффициент Пуассона v. Неизвестной величиной является безразмерная частота ш = tw. Для упрощения обозначений знак «тильда» опускается. Далее в расчетах используются следующие значения: v = 0.3, 0 = 20.

Для повышения эффективности алгоритма численного решения можно учесть свойства симметрии задачи. В связи с тем, что поперечное сечение балки имеет две оси симметрии, система уравнений (26), (27) и (29) может быть разделена на четыре независимые подсистемы в соответствии с четностью полиномиальной

130

В.В. САУРИН, Г.В. КОСТИН

Табл. 1

Свойство симметрии базисных функций ~ xN2 xN3

Растяжение вдоль Ox1 Изгиб вокруг Ox2 Изгиб вокруг со 0 Кручение вокруг Ox1

N2 N3 N2 N3 N2 N3 N2 N3

w1, 7jj 2n 2n 2n 2n + 1 2n + 1 2n 2n +1 2n + 1

W2, 712 2n — 1 2n — 2 2n — 1 2n — 1 2n 2n 2n 2n + 1

W3, 713 2n — 2 2n — 1 2n 2n 2n — 1 2n — 1 2n + 1 2n

723 2n — 3 2n — 3 2n — 1 2n — 2 2n — 2 2n — 1 2n 2n

части базисных функций из (34), приближенно описывающие растяжение - сжатие, изгибы относительно осей 0x2 и Охз, а также кручение балки [2].

Для различных собственных движений балки в табл. 1 приведены максимальные степени N2(n) и Ns(n) переменных Х2 и хз, которые присутствуют в аппроксимациях перемещений и напряжений (34). Индекс j принимает значения 1, 2, 3, а целое n > 0 определяет дифференциальный порядок соответствующей краевой задачи. Четность (нечетность) чисел N2 и N3 характеризует свойства симметрии (антисимметрии) функций перемещений и напряжений относительно координатных плоскостей 0x2 и Охз соответственно. Если хотя бы одно из этих чисел меньше нуля, то соответствующие функции отсутствуют.

Порядок системы дифференциальных уравнений равен: (n + 1)(3n + 2) - для задачи растяжения-сжатия; (n + 1)(3n + 4) - для задач изгиба; (n + 1)(3n + 6) -для задачи кручения. Минимально возможные размерности аппроксимаций в представлении (34) будут соответственно равны 2, 4 и 6.

Вследствие того, что в настоящей статье рассматривается случай, когда конец балки при xi = —ai движется согласно (16) в направлении оси 0x2, в системе возникают только изгибные колебания относительно оси 0x3.

В качестве примера рассмотрим задачу нахождения собственных частот изгиб-ных колебаний балки с квадратным сечением при n = 0 . Для этого случая функции из (34) имеют вид

Wi = w(10) (xi)x2, W2 = ^Z’(00)(xi),

7ii = ?(10)(xi)x2, <712 = 5(20)(xi)(1 — x2), (35)

<722 = 7(20)(xi)x2(1 — x2), W3 = <713 = 723 = 733 = 0.

Тогда результирующая система дифференциально-алгебраических уравнений (26), (27) и (29) может быть представлена как

4 da

(10)

11

8,

3 dx1 3

(00)

(00)

г12

2~(10)

4— <а га) 3 1

16 da1T + 16 7(10) =0

45 dx1 + 15722 ,

4 d^(10) 4~ (10) +4~ (10) =

3 dx1 3 11

с граничными условиями

+ — а(2 = о,

25 22 ’

8 da

0, --

(00)

12

3 dx1

2 (00) + 4^2 га2

1*2°°. +4 — 208 „ =0

3 dx1 3 1 75 12

( 0)

A (—ai)

(00)

га2 ч—ai)

а11

(10)(ai)

a(020)(a1) = 0.

4

0

(36)

(37)

ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

131

Табл. 2

Собственные частоты консольной балки ut

i 1 2 3 4

ut 1.269 • 10~3 7.951 • 10~3 2.226 • 10~2 4.363 • 10~2

n = 0, ui 1.266 • 10~3 7.843 • 10~3 2.157 • 10~2 4.123 • 10~2

n =1, ut 1.269 • 10~3 7.861 • 10~3 2.162 • 10~2 4.135 • 10~2

(ut - u°)/u° 0.2% 1.4% 3.2% 5.8%

(ut - u°)/u° 0.21% 0.22% 0.25% 0.29%

At 0.30% 0.31% 0.32% 0.33%

Собственные частоты ш находятся из системы (36) как корни характеристического уравнения

4

3

406

104

А4 4----ш2А2 — 4ш2 4----ш4

75

25

0

(38)

где А - соответствующее волновое число.

Для сравнения приведем характеристическое уравнение для балки Тимошенко с такими же параметрами [9]:

4

3

412

104

А4 + —— ш2А2 — 4ш2 + —— ш4

75

25

0.

(39)

Как видно, уравнения (38) и (39) отличаются только коэффициентами при А2 менее чем на 2%.

Для более точного вычисления собственных частот и форм колебаний балки следует применять более высокие степени полиномиальных аппроксимаций n > 0. В первых трех строках табл. 2 приведены значения четырех низших собственных частот колебаний консольной балки квадратного сечения соответственно для модели Эйлера-Бернулли, а также для предложенной в статье модели при n = 0 и n =1 .В четвертой и пятой строках показано различие частот, полученных в рамках этих приближений. Различие между полученными частотами и частотами, вычисленными согласно модели Эйлера-Бернулли, значительно и достигает 5.8% уже для четвертой моды. Последняя строка отражает относительную точность Д1 вычисления i -й собственной формы колебаний при n =1 согласно интегральному критерию (10).

4.3. Система обыкновенных дифференциальных уравнений по времени. Для того чтобы построить систему обыкновенных дифференциальных уравнений относительно времени t, подставим в аппроксимации (33) вместо функций перемещений wi(x) и напряжений Xi(x) собственные формы колебаний из (34), соответствующие частотам шi, i = 1,...,M — 1. Полученные аппроксимации подставляются в интегральные уравнения (26), (27) и (29). В соотношениях (26), (27) вектор vt и тензор £ проецируются на базисные вектор-функции ri(x) и si (x) с компонентами

N

N +1

r1,i = ^14)(Х1)Ж2 ^ r = V4 Xkl)(x- )rkxi rn,i = / y rn,i (x1)x2x3,

k+l = 0 k+l=0

N N-1

ri,i = riki)(xi)rk ^ :Xnk,i)(x 1)xk x3.

k+l=0 k+l=0

(40)

Неизвестные функции rjkp(xi) и sjkp(x 1), j = 1, 2, 3, находятся как решение сопряженной краевой задачи на собственные значения [10].

132

В.В. САУРИН, Г.В. КОСТИН

Выбор таких проекций приводит интегральные уравнения (26), (27) и (29) к диагональному виду

yj = + Ьм+j u(t), j = 2,...,M, (41)

где wi, i = 1,...,M — 1, — приближенные собственные частоты, полученные для выбранной степени аппроксимации n, а bk, к = M + 1,..., 2M, - коэффициенты управления.

5. Конечномерная задача управления

Учитывая начальные условия (17), (32), перепишем систему уравнений (41), добавив к ней уравнение относительно yi(t) из (16), в виде

yj = Ум+j, j = 1,...,M,

Ум+j = —wjj-iyj + Ьм+j u(t), j = 1,...,M, (42)

yk (0)=0, к =1,..., 2M.

Здесь введены новые переменные ум+j (t), являющимися производными по времени от неизвестных yj (t), а также частота wo = 0 и коэффициент управления

Ьм +i = 1.

Систему (42) можно представить в векторном виде

y(t) = f (У, u) = Ay(t) + bu(t) (43)

с однородными начальными условиями

y(0) = 0 (44)

и вытекающими из (19) терминальными соотношениями

У(Т ) = (ут , 0,..., 0)T, (45)

где y = (yl(t),..., у2м(t))T - вектор фазовых переменных, A G М2мх2м и Ь2м -постоянные матрица и вектор.

Для конечномерной динамической системы (43), которая в итоге приближенно определяет поперечные движения рассматриваемого упругого тела, сформулируем задачу оптимального управления, соответствующую задаче (18)—(20): найти функцию управления u*(t), которая в фиксированный момент времени T переводит линейную систему (43) из нулевого состояния (44) в терминальное состояния покоя (45) и минимизирует функционал качества

J\u*] = min J\u\.

u£L2(0,T)

(46)

Здесь квадратичный интеграл

T

J= j fo(t)dt, fo = 2 u2(t) + - y(t)TWy(t), (47)

0

получен в результате дискретизации функционала J, введенного в (21).

Введением вектора сопряженных переменных z(t) G М2м можно согласно принципу максимума Понтрягина [11] определить гамильтониан системы

H\y, z,u] = —fo + fTz.

(48)

ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

133

Рис. 2. Закон оптимального управления u(t)

Используя гамильтониан (48), запишем сопряженную систему уравнений в виде

Z(t = -dH = yWy(t) - ATz(t), (49)

при этом оптимальное управление как линейная функция сопряженных переменных есть

u = bT z (t). (50)

После подстановки этой функции в уравнение (43) нахождение оптимальных движений сводится к решению краевой задачи (43)-(45), (49), (50).

6. Численные результаты

Для иллюстрации работы предложенного алгоритма нахождения управляемых движений упругого тела призматической формы выбраны следующие размерности аппроксимации n = 1, M = 4, 5 и параметры управления ут = 20, T = 5000, y = = 5 • 10-7. Время управления T несколько превышает период колебаний балки по первой моде Ti ~ 4952, перемещение ут составляет половину безразмерной длины балки a.

На рис. 2 представлен оптимальный закон управления u(t), полученный для параметра приближения M = 5 (учитываются только первые 4 моды колебаний). Эта функция слабо отличается от управления, полученного для M = 4, что демонстрирует рис. 3, на котором представлена разность между этими двумя законами Au(t). Как следует из обоих рисунков, учет дополнительной моды колебаний в приближенной модели добавляет в управление незначительную высокочастотную составляющую.

На рис. 4 представлены перемещения w2 (t, x±) двух точек упругого тела с координатами x± = (±a, 0, 0) в фазовой плоскости. Штриховой кривой показаны положение и скорость точки, лежащей на середине граничного сечения балки (x1 = = -a), которое перемещается как жесткое целое. Сплошной кривой обозначена фазовая траектория точки x1 = a, помещенной в середину граничного сечения, свободного от внешних нагрузок. Как можно заметить по отклонению этих кривых друг от друга, при реализации оптимального закона управления для выбранного набора параметров в системе возбуждаются значительные упругие деформации.

Распределение энергии по модам колебаний показано на рис. 4. Сплошная кривая E1 (t) - это изменение кинетической энергии, которая соответствует поступательному движению балки. Она пропорциональна квадрату от переменной у1 (t).

134

В.В. САУРИН, Г.В. КОСТИН

Рис. 3. Изменение управления Au(t) при повышении размерности задачи

Рис. 4. Перемещение в фазовой плоскости точек граничных сечений х±

Рис. 5. Изменение кинетической энергии низших мод Ei(t)

Штриховая кривая E2 (t) соответствует кинетической энергии первой, а штрих-пунктирная E3 (t) - второй моде колебаний. Величины этих энергий пропорциональны соответственно значениям y\ (t) и y2 (t). Остальные энергии не приведены на графике из-за малости их величин. Основные упругие движения относятся к первой моде и сравнимы по энергии с движением балки как целого.

ПОДХОДЫ К ПОСТРОЕНИЮ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

135

Для более точного построения оптимального закона управления движением упругого тела u*(t) и нахождения соответствующих полей перемещений, напряжений и импульсов необходимо, если позволяют вычислительные ресурсы, согласовано повышать порядок полиномиального разложения n и количество учтенных мод M. Стоит отметить, что повышение размерности модели приводит к появлению высокочастотных колебаний функции управления, что может существенно осложнить реализацию полученных законов в технических приложениях.

Работа выполнена при финансовой поддержке РФФИ (проекты № 13-01-00108, 14-01-00282, 15-01-00827) и грантов поддержки ведущих научных школ НШ-2710.2014.1, НШ-2954.2014.1.

Summary

V.V. Saurin, G.V. Kostin. Variational Approaches to Optimal Control Design for Elastic Body Motions.

Variational and projection statements of an initial-boundary value problem for the direct and inverse dynamics of elastic bodies are proposed. An optimal control problem of 3D linear elastic motions for rectilinear beams with the rectangular cross section is studied. Based on the generalized formulations, a numerical algorithm is developed to design the optimal displacement of such elastic beams. The results of numerical simulation and quality analysis are presented and discussed.

Keywords: optimal control, dynamics, linear theory of elasticity, variational principles, semi-discretization.

Литература

1. Kostin G.V., Saurin V.V. Integrodifferential Relations in Linear Elasticity. De Gruyter Studies in Mathematical Physics 10. - Berlin: De Gruyter, 2012. - 280 p.

2. Костин Г.В., Саурин В.В. Моделирование и анализ собственных колебаний упругой призматической балки на основе проекционного подхода // Прикл. матем. и механика. - 2011. - Т. 75, Вып. 6. - С. 995-1010.

3. Костин Г.В. Построение оптимального управления движением упругих тел методом интегродифференциальных соотношений // Изв. РАН. ТиСУ. - 2007. - № 4. -С. 21-31.

4. Saurin V.V., Kostin G.V., Rauh A., Aschemann H. Variational approach to adaptive control design for distributed heating systems under disturbances // Int. Rev. Mech. Eng. - 2011. - V. 5, No 2. - P. 244-256.

5. Kostin G.V., Saurin V.V., Aschemann H., Rauh A. Integrodifferential approaches to frequency analysis and control design for compressible fluid flow in a pipeline element // Math. Comput. Model. Dyn. Syst. - 2014. - V. 20, No 5, - P. 504-527.

6. Тимошенко С.П., Гудьер Дж. Теория упругости. - М.: Наука, 1979. - 560 с.

7. Тихонов А.Н., Самарский А.А. Уравнения математической физики. - М.: Наука, 1977. - 735 с.

8. Лионс Ж.-Л. Оптимальное управление системами, описываемыми уравнениями с частными производными. - М.: Мир, 1972. - 416 с.

9. Тимошенко С.П. Колебания в инженерном деле. - М.: Наука, 1967. - 444 с.

10. Данфорд Н., Шварц Дж. Линейные операторы. - Т. 2. Спектральная теория. Самосопряженные операторы в гильбертовом пространстве. - М.: Мир, 1966. - 1064 с.

136

В.В. САУРИН, Г.В. КОСТИН

11. Понтрягин Л.С., Болтянский В.Г., Гамкрелидзе Р,В., Мищенко Е.Ф. Математическая теория оптимальных процессов. - М.: Наука, 1976. - 392 с.

Поступила в редакцию 15.06.15

Саурин Василий Васильевич - доктор физико-математических наук, старший научный сотрудник, Институт проблем механики им. А.Ю. Ишлинского РАН, г. Москва, Россия.

E-mail: saurin@ipmnet.ru

Костин Георгий Викторович - доктор физико-математических наук, старший научный сотрудник, Институт проблем механики им. А.Ю. Ишлинского РАН, г. Москва, Россия.

E-mail: kostin@ipmnet.ru