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

УДК 531.132.2

С. А. Борисевич, ассистент (БГТУ)

ИССЛЕДОВАНИЕ ДВИЖЕНИЯ ГИБКОГО СТВОЛА ДЕРЕВА С ЕСТЕСТВЕННОЙ КРИВИЗНОЙ В ТРЕХМЕРНОМ ПРОСТРАНСТВЕ

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

A rod-chain model was used to study the three-dimensional dynamics of a tree stem. A numerical method is developed for simulating the mechanical behavior of a flexible tree stem. The motion of the flexible tree stem is determined under the resistance forces exerting on the tree crown. The modeling strategy of this approach employed the exact nonlinear kinematic relationships in the sense of Cosserat theory, and adopted Bernoulli hypothesis. A deformed configuration of the rod is described by the displacement vector centroid curves and an orthonormal moving frame, rigidly attached to the cross-section of the rod. The position of the moving frame relative to the inertial frame is specified by the rotation matrix, parameterized by rotational quaternions. A falling tree with natural curvature as a simple example has been presented to illustrate the use of developed model.

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

В данной статье рассмотрена численная модель механического поведения гибкого ствола дерева при его движении в трехмерном пространстве, основанная на модели [2].

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

декартову систему координат (х, у, 2) в некотором инерциальном базисе (е1, е 2, е3) (рис. 1). Ориентацию поперечного сечения в произвольном сечении 5 будем задавать ортогональным базисом d 1 (5, t) (/ = 1, 2, 3), называемым подвижным, причем единичный вектор d3 - нормальный к поперечному сечению, а d1 и d2 лежат в плоскости сечения. Тогда движение участка стержня описывается его нейтральной линией г(5, t) и тремя ортогональными единичными векторами di (5, t) (/ = 1, 2,3).

d

ei

Рис. 1. Модель стержня Коссера

В инерциальном декартовом базисе (е1, е 2, е3) можно записать:

г(5, t) = х(5, t)е + у(5, t)е2 + 2(5, t>3.

Движение стержня включает скорость осевой линии дг (5, t) / дt и угловую скорость поперечного сечения w(5, t). Вектор угловой скорости поперечного сечения может быть представлен в виде [3]

1 d dd,

w = — d, x—-2 ' dt

а вектор деформации -

1 d dd

u = — d; x—L

2 г ds

(2)

Компоненты нормированы

IЯ1 = 1.

(5)

Поскольку базис ^ (5, t) (г = 1, 2, 3) естественный для внутреннего описания деформаций, разложим все вектора в нем:

3

Vt) = Х^ ( tН ^ t)

г=1

3

" ^ t) = Хиг ^ tН ( t),

г=1

3

V ( t) = Х ^ ^ tН ( t)•

Параметризация матрицы поворота. Математическое моделирование пространственных поворотов, их параметризация и линеаризация рассматриваются в ряде публикаций, которые появились в последние годы [3, 4]. Эти работы являются развитием предыдущих работ, исследующих движение гибких стержней и нитей в плоскости. Ранее было доказано, что свободные от сингулярности уравнения поворота не могут быть получены при использовании трех независимых переменных, и поэтому было предложено [5] использовать четыре компоненты кватерниона в качестве обобщенных координат. Для параметризации матрицы поворота применим компоненты кватернионов q = q(s) = (q^ q2, q3, q4). Тогда выражения для единичных векторов примут вид

di (q) = (q2 - q22 - q32 + q2, 2qiq2 + 2qзq4,

2qiq3- 2q2 q4),

d 2 (q) = (2qq2- 2qзq4, - q2 + q2- q3 + q42,

2q2 q3- 2qq4),

d3(q) = ( 2qq3 + 2q2 q4,2q2 q3- 2qlq4,

2 2 2 2 \ -qi -q2 + q3 + q4).

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

. Гф + у

ql =sin I ^ I cos I —

(3)

'eV (ф + V

42 = sin| 2 I sin

'eV (Ф + v

43 = cos| - IsinI

'e^ (Ф + v

44 = cosi 21cosI

(4)

Подставляя выражения для единичных векторов из (3) в (1) и учитывая выражение для нормы кватернионов (5), получим проекции угловой скорости в связанной системе координат:

®1 = 2 (ъъ+<?2 Ъ - ъЗъ + q4 Ъ), ®2 = 2 ((к - + ъЗъ + Ъ4Ъ ), ®3 = 2(( - Ъ - ъзЪ4 + Ъ4Ъ )•

Аналогично получим компоненты деформации в связанной системе координат из выражения (2):

и1 = 2 (+ ъ3- ъЗ + Ъ4 ъ'), и2 = 2 (ъ3- ъ4 + ъзъ1 + Ъ4 из = 2 (2- - ъЗ ъ4 + Ъ4 ъ3) •

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

Уравнения движения ствола дерева. Располагая локальную ортогональную систему единичных векторов в конечном

числе точек, получим модель дерева в виде цепочки жестких стержней, соединенных упругими шарнирами. Свяжем с каждым стержнем модели локальную ортогональную систему единичных векторов d2, d3). Положение каждого стержня в некоторой неподвижной системе координат (е1, е2, е3) будем задавать радиус-вектором центра масс стержня:

гс (s, 0 = хс (5 t)el + ус (5 t)е

2 + 2С (5,t)ез

и ориентацию каждой локальной системы координат d2, d3) относительно неподвижной системы координат (е1, е 2, е3) будем определять с помощью выражения (3).

Координаты центров тяжести стержней в неподвижной системе координат найдем следующим образом:

к-1

хск=Е(З(чХ ) Д=

г =1

= ]Г1 ( + 2^4) + 11к (к + 2ъкъ4к),

г=1 2

к-1

Уск =Е(З(ЧХ )у Д5 =

г =1

= 1I ( Ъз - 2ъ1 ъ4 )+21к (к - 2« ),

m

i=1

i=1

k-1

k-1

^=Е(3(чх-) Д=

1=1

ЕI (-(Ч )2 - (ч2 )2 + (Ч3 )2 + (Ч4 )2)+

г =1

+| 4 (()2 - (чк )2 + (чк )2 + (чк )2),

где ¡^ - длина 1-го стержня.

Скорость центров масс стержней равна

2 = • 2 + • 2 + • 2 УСк = ХСк + уСк + 2Ск .

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

Т = ткУСк + 2(( ( ) + ))|,

где юхк, юук, ю2к - проекции угловой скорости на оси подвижной системы координат, связанной с телом.

Потенциальная энергия системы состоит из энергии в поле сил тяжести и энергии упругих связей между стержнями. Потенциальная энергия упругих связей имеет вид

П = 1Е ЕС;к (к -и0к )2,

2 1 =1 к=2

где с к - жесткости упругих связей; и]к - компоненты деформации стержня; и 0к - компоненты деформации, соответствующие естественной кривизне стержня.

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

* д+1- д п„ = д+1- 2% +

1 ' Ч - . 2 ;

Д5 Д5

где Д5 - шаг дискретизации.

Потенциальная энергия жестких цилиндров в поле сил тяжести:

к=1

Сила сопротивления, действующая на крону дерева, принята линейной по скорости Ек = -Руск, так как это сделано в работе [2], введена диссипативная функция Релея:

п . 2

ф = 1РЛ .

1=1 2

Тогда обобщенная сила сопротивления равна

о =-дФ ¿¿к ~ -. •

дЧк

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

Ц*=Е Чт («) -1 (« = 1,2,..., п ),

т

которая фактически является условием выполнения нормы кватерниона для каждого стержня модели. Составим систему уравнений Лагранжа второго рода с множителями (следуя методу неопределенных множителей Лагранжа):

( дТ ^

dt

dqt

k J

dT d4k

дП

дЧк

дФ

дЧк

Я

а=1 дЧк

(k = 1,2,..., n),

которая совместно с s = n уравнениями связей служит для определения 4n + s неизвестных величин:

^(1), Я2(1),..., ?4(1), (2),..., ?4(2),...,

(n),..., q^n), »1,..., »n.

Уравнения движения здесь не выписываются ввиду их чрезвычайной громоздкости. Множители Лагранжа в этих уравнениях не имеют физического смысла, и их равенство нулю говорит о правильности вычислений.

Полученные уравнения решены численно с помощью пакета Maple 11 при заданных начальных углах отклонения стержней и их начальных скоростях, пересчитанных в компонентах кватернионов. С учетом начальных условий получена система линейных дифференциальных уравнений. Из данной системы уравнений находились старшие производные компонент кватернионов при помощи метода Ньютона - Раф-сона [6]. Численное интегрирование уравнений движения производилось с помощью модификации базовой схемы Верле, так называемой полушаговой «leap-frog» схемы [5].

Динамические параметры стержней и характеристики пружин сочленения определялись аналогично тому, как это делалось в работе [2].

Пример расчета. В качестве примера рассмотрим падение ствола дерева переменного сечения (усеченный конус), изогнутого в одной

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

16 14 12 10 8 6 4 2

6

Рис. 2. Формы осевой линии ствола дерева в различные моменты времени

Начальное положение стержней в пространстве задавалось с помощью углов Эйлера, а затем по их значениям находились соответствующие компоненты кватернионов при помощи формул (4). Начальные угловые скорости относительно каждой из осей принималась равными нулю.

На рис. 2 представлены формы осевой линии ствола дерева в различные моменты времени. Из рисунка видно, что в процессе движения проекция осевой линии на горизонтальную плоскость поворачивается, несмотря на то, что начальная угловая скорость равна нулю.

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

Литература

1. Жуков, А. В. Теория лесных машин / А. В. Жуков. - Минск: БГТУ, 2001. - 640 с.

2. Борисевич, С. А. Конечно-разностная схема для исследования падения ствола дерева / С. А. Борисевич // Труды БГТУ. Сер. II, Лесная и деревообраб. пром-сть. - 2008. - Вып. XVI. -

C.104-107.

3. Cao, D. Q. Three-dimensional nonlinear dynamics of slender structures: Cosserat rod element approach / D. Q. Cao, D. Liu, C.H.-T. Wang // International Journal of Solids and Structures. -2006. - Vol. 43, № 3/4. - P. 760-783.

4. Zupan, E. The quaternion-based three-dimensional beam theory / E. Zupan, M. Saje,

D. Zupan // Comput. Methods in Appl. Mech. Engrg. - 2009. - Vol. 198, № 49-52. - P. 39443956.

5. Allen, M. P. Computer Simulation of Liquids / M. P. Allen, D. J. Tildesley. - Oxford: Clarendon press, 1999. - 385 p.

6. Дэннис, Дж. Численные методы безусловной оптимизации и решения нелинейных уравнений / Дж. Дэннис, Р. Шнабель. - М.: Мир, 1988. - 440 с.

Поступила 01.03.2011