CC BY
44УДК 630.3 Статья
Алгоритм декомпозиции линейной задачи дополнительности и его применение для моделирования соударений балансов в корообдирочном барабане
Г еннадий Н. Колесников
Петрозаводский государственный университет, пр. Ленина 33, 185910 Петрозаводск, Россия;
E-Mail: kolesnikovgn@ya.ru
Tel.: +7(8142)711039; Fax: +7(8142)711000.
Получена: 13 сентября 2013 /Принята: 08 ноября 2013 / Опубликована: 16 ноября 2013
Аннотация. Проблема численного моделирования соударений балансов в корообдирочном барабане рассматривается как источник контактной задачи. Решение данной задачи необходимо для прогнозирования силы соударений балансов друг с другом и с корпусом корообдирочного барабана. С применением конечно-элементной модели контактная задача сведена к поиску решения линейной задачи дополнительности с положительно определенной матрицей коэффициентов. Динамика объекта исследования моделируется как множество его статических состояний, упорядоченных во времени с помощью метода конечных разностей. Обоснован новый алгоритм декомпозиции линейной задачи дополнительности с положительно определенной матрицей коэффициентов. На каждом шаге во времени численная реализация алгоритма сводится к исключениям Гаусса-Жордана с выбором разрешающего элемента по энергетическому критерию. Упомянутый критерий используется для поиска глобального минимума дискретного множества квадратичных функций, каждая из которых имеет одну переменную. С точки зрения механики энергетический критерий используется как критерий перехода односторонних связей из возможного состояния в действительное. Решение указанной выше линейной задачи дополнительности по новому алгоритму определяется за меньшее число шагов, чем при использовании описанных в литературе алгоритмов.
Ключевые слова: негладкая динамика, соударения, корообдирочный барабан, численное моделирование, линейная задача дополнительности, энергетический критерий
Article
Decomposition Algorithm of Linear Complementarity Problem and its Application for Simulation of Pulpwood Collisions in Debarking Drum
Gennady N. Kolesnikov
Petrozavodsk State University, Lenin av. 33, 185910 Petrozavodsk, Russia;
E-Mail: kolesnikovgn@ya.ru
Tel.: +7(8142)711039; Fax: +7(8142)711000
Received: 13 September 2013 /Accepted: 08 November 2013 /Published: 16 November 2013
Abstract: The problem of numerical simulation of pulpwood collisions in a debarking drum is regarded as the source of the contact problem. The solution of this problem will allow to predict the strength of pulpwood collisions with each other and with the drum body. The contact problem reduces to finding a solution of a linear complementarity problem with a positive definite matrix of coefficients if a finite element model is applied. The dynamics of the research object is modeled as a set of its static states and time-ordered by the method of finite differences. A new decomposition algorithm of a linear complementarity problem with a positive definite matrix of coefficients is substantiated. At each time step the algorithm numerical implementation reduces to Gauss-Jordan eliminations with the choice of the matrix main element by the proposed energy criterion. The energy criterion is used to find the global minimum of a discrete set of quadratic functions, where each one has a single variable. In terms of mechanics, the energy criterion at each step of the algorithm is used as the criterion for the transition of unilateral constraints from a possible state to a real state. The solution of the mentioned linear complementarity problem by the new algorithm is found in fewer steps in comparison to test examples known from the literature.
Keywords: non-smooth dynamics, collision, debarking drum, numerical modeling, linear complementarity problem, energy criterion
1. Введение
Необходимость моделирования соударений в механических системах часто появляется в инженерных приложениях [1,2]. Только одно из приложений кратко рассматривается в данной работе, чтобы проиллюстрировать применение предлагаемого алгоритма к исследованию соударений деформируемых тел. Это приложение относится к прикладной механике технологических процессов дробления и разрушения, где использование численного моделирования соударений позволяет получить определенный технический эффект [4]. К данной области прикладных исследований относится многоплановая проблема моделирования соударений при очистке балансов от коры в установках барабанного типа [5]. В таких установках кора разрушается при соударениях балансов и посредством сил трения при их взаимодействии друг с другом и с корпусом барабана, в результате вращения последнего вокруг своей продольной оси. Именно эти два типа взаимодействия и приводят к очистке балансов от коры[6]. При этом технология очистки древесины с помощью корообдирочных барабанах остаётся доминирующей [7,8], развиваясь параллельно с другими технологиями окорки [9-11].
Для решения затронутой проблемы необходимы методы прогнозирования силы соударений балансов [5-7,12]. Если сила соударений избыточно велика, то разрушается не только кора, но и древесина, а также повреждаются торцы балансов [13]. Если уменьшить силу соударений, то увеличится продолжительность процесса или уменьшится степень очистки балансов [14]. Очевидно, величина силы соударений балансов должна находиться в определенном интервале значений [15]. Точное выполнение этого условия практически недостижимо по причинам вариабельности свойств древесины и случайного характера соударений балансов. Однако для определенного вида балансов можно экспериментально подобрать такие значения степени заполнения барабана и скорости его вращения, которые обеспечат наилучшее соответствие данному условию [8,16]. Для расширения появляющихся в этой связи возможностей уточнения технологических параметров окорки необходима методика численного прогнозирования сил контактного взаимодействия при соударениях балансов друг с другом и с корпусом корообдирочного барабана.
Возможности экспериментальных исследований процесса соударений тел в реальных технологических процессах ограничены. Поэтому для определения сил взаимодействия при соударениях целесообразно использовать методы математического моделирования [17,18]. Алгоритмы решения таких задач представлены, например, в обзоре [19]. Новые результаты обсуждаются в статьях [5,20-22 и др.].
В данной работе рассматривается обоснование подхода, в котором задача о соударениях деформируемых тел сводится к поиску решения линейной задачи дополнительности. При этом для решения линейной задачи дополнительности предлагается алгоритм декомпозиции, базирующийся на идее работ [23-25]. В рассматриваемой далее проблеме соударений все или некоторая часть ограничений на перемещения тел являются односторонними [23,26]. Такие
ограничения (связи) в полукоэрцитивных контактах препятствуют взаимопроникновению тел, но не препятствуют появлению зазоров. Появление зазоров изменяет структуру механической системы. Коэрцитивный контакт не приводит к появлению зазоров и не изменяет структуру механической системы (Рис. 1). Более детально эти аспекты рассмотрены, например, в статье [20].
(а) (Ь)
Рисунок 1. (a) Коэрцитивный контакт. (b) Полукоэрцитивный контакт Figure 1. (a) Coertive contact. (b) Semicoercive contact
Для каждой односторонней связи возможны два состояния: «включено» («есть контакт») и «выключено» («нет контакта»). Реализуется только то состояние, в котором потенциальная энергия механической системы минимальна. Если в механической системе число односторонних связей равно , то число возможных состояний модели механической системы равно 2 п. Условие минимума потенциальной энергии позволяет выбрать из всех возможных состояний то действительное состояние, которое является единственным. Поиск действительного состояния односторонних связей представляет собой достаточно сложную проблему, которая исследуется многими авторами [19,20,22 и др.]. Известные в данной области алгоритмы не всегда удобны на практике (см., например, [27]).
Сложность проблемы объясняется тем, что траектория движения каждого тела состоит из участков контактного и бесконтактного движения (рис. 1). Соударения сопровождаются изменениями структуры механической системы и, соответственно, появлением негладкой функции скорости. Решение можно найти, например, методом припасовывания, который сводится к поэтапному интегрированию уравнений движения на отдельно взятых участках, а начальные условия для участка определяются по конечным значениям переменных на предыдущем участке. Однако “практическое применение метода припасовывания осложняется тем, что моменты ударов заранее неизвестны, а само «припасовывание» может быть сложно выполнимым технически” [18, с. 4]. Многообразие методов негладкой динамики представлено в работах [2,19,22,27,28]. В частности, перебор всех возможных ( 2п) вариантов применяются при небольших значениях .
Анализ литературы показывает, что проблема соударений многих тел сохраняет свою актуальность и требует продолжения исследований в инженерных приложениях. Цель данной работы: обоснование эффективного в вычислительном отношении алгоритма декомпозиции линейной задачи дополнительности и его применение в численном моделировании соударений деформируемых тел на примере определения верхних оценок сил контактного взаимодействия при соударениях балансов в корообдирочном барабане.
2. Материалы и методы
2.1. Физическая модель объекта исследования
Определение сил контактного взаимодействия при соударениях балансов в корообдирочном барабане предполагает выбор физической модели объекта исследования. По причине сложности объекта исследования и различий целей работ в литературе [6,9,11,16, 29,30] используются различные подходы к выбору физической модели взаимодействия балансов в корообдирочном барабане. В данной работе используется физическая модель в виде условного штабеля балансов, обоснование которой приведено в работе [15]. Некоторые результаты применения этой модели опубликованы в статье [5]. Предполагается, что: траектория движения каждого тела состоит из участков контактного и бесконтактного движения; соударения рассматриваются с применением модели центрального удара; деформации баланса при ударе существенно меньше его диаметра; масса каждого тела сосредоточена в центре тяжести; точки на границе тел «безмассовые»; перемещение этих точек относительно центра масс тела зависит от жесткости [31]. Изменениями массы балансов в процессе очистки от коры пренебрегаем. При соударениях тела без взаимопроникновения соударяются друг с другом.
Величина зазора Dj (i = 1,2 , 3 для рисунка 2) и сила контактного взаимодействия V связаны соотношением D¿V = 0, где в соответствии с правилом знаков [17] Dj >0, V > 0; равенство = 0 моделирует то обстоятельство, что одна из величин ( или )
обязательно равна нулю. Диссипация энергии при соударениях тел определяется коэффициентами вязкости Лу (_/ = 1,2 для рисунка 2). В фазе бесконтактного движения Лу = 0.
Использование данной модели в работе [5] позволило численно уточнить некоторые детали взаимодействия еловых балансов неодинакового диаметра в корообдирочном барабане. Результаты работы [5] демонстрируют эффективность использованной методики численного моделирования и указывают на целесообразность развития данного направления прикладных исследований. Для обоснования математического описания физической модели необходимы нижеследующие предварительные замечания.
Ur
1
u4
з
U
2
hi
ш
Рисунок 2. (a) Модель соударяющихся тел. (b) Номера узлов и перемещения Figure 1. (a) The model of colliding bodies. (b) The node numbers and movements
2.2. Предварительные замечания
Пусть в некоторой механической системе число возможных контактных пар равно n. В фазе контактного движения в некоторой контактной паре i зазор Dj = 0; при этом сила контактного взаимодействия iVj > 0 . Сила считается неотрицательной в соответствии с правилом знаков, обычно используемым в контактных задачах [17].
Определение 1. Связь i находится в состоянии «включено», если Vj > 0 и D j = 0.
Определение 2. Связь i находится в состоянии «выключено», если Vj = 0 и D j >0.
Следовательно, в любом состоянии односторонней связи i выполняются соотношения
V >0, D j > 0, D jVj = 0.
Утверждение 1. Двустороннюю связь можно интерпретировать как частный случай односторонней связи, которая всегда находится в состоянии «включено».
Утверждение 2. Отсутствие связи можно интерпретировать как частный случай односторонней связи, которая всегда находится в состоянии «выключено».
Из этих утверждений следует: если в механической системе действительное состояние каждой односторонней связи известно, то анализ (расчет) такой системы не отличается от анализа механической системы с двусторонними связями.
Определение 3. Определим элементарное событие, которое с физической точки зрения представляет собой переход односторонней связи из состояния «включено» («выключено») в состояние «выключено» («включено»).
Элементарное событие сопровождается скачкообразным переходом системы из одного состояния в другое состояние, или с одного уровня потенциальной энергии механической системы на другой, что характерно для квантовых систем.
Потенциальная энергия механической системы определяется деформациями тел, состоянием односторонних связей, зависит от внешних и внутренних сил. Условие минимума потенциальной энергии позволяет найти действительное состояние односторонней связи. Если действительное состояние всех односторонних связей найдено, то дальнейший расчет может быть сведён к решению системы линейных алгебраических уравнений [17, с. 155]. Таким образом, основная проблема - определение действительного состояния односторонних связей. На этой проблеме фокусируется внимание в данной статье.
Установление воздействия на механическую систему сопровождается, в общем случае, изменением состояния односторонних связей. В процессе установлении воздействия односторонние связи переходят из возможного состояния в действительное; состояние механической системы трансформируется так, что в любой точке процесса достигается минимум потенциальной энергии. Этот процесс рассмотрим как очередность указанных выше элементарных событий.
Для обоснования предлагаемого далее алгоритма ключевую роль играет гипотеза о том, что элементарные события реализуются в определенной очередности [23]. Это означает, что на некотором достаточно малом отрезке времени осуществляется только одно элементарное событие. Взамен времени может быть использован другой подобный времени параметр, например, параметр нагрузки при однопараметрическом воздействии [23]. Определив эту очередность, мы найдем также действительное состояние каждой односторонней связи. На каждом отрезке времени между элементарными событиями механическая система представляет собой линейно деформируемую конструкцию, что позволяет использовать на данном отрезке хорошо известный метод конечно-элементного анализа механической системы с двусторонними связями, который сводится к решению системы линейных алгебраических уравнений [33,34].
2.3. Математическое описание физической модели
В простейшем случае каждое тело моделируется двумя конечными элементами. Матрица жесткости г тела / (Рис. 2, / = 1,2 ) имеет стандартный вид [33]:
1 - 1 0
г = 5* - 1 2 - 1
0 -1 1
(1)
Физической модели с безмассовыми узлами 1 и 3, 4 и 6 на интерфейсной границе того же тела (Рис. 2) соответствуют нижеследующие диагональные матрицы масс т и демпфирования А:
0 0 0 0 0 0
Л — 0 Лі 0 ; т — 0 Мі 0
0 0 0 0 0 0
В данном случае использована простейшая модель демпфирования. Более детально этот вопрос рассмотрен в [35]; в ряде случаев матрица демпфирования может быть представлена в виде линейной комбинации матриц масс и жесткости [34].
Выполнив в (1) и (2) перестановку строк и столбцов (что потребует изменения нумерации искомых величин), запишем:
(3)
"2 " " Гд. 0 0 Г Мі 0 0
— и " 1 0 ; л — 0 0 0 ; т — 0 0 0
" 0 1 0 0 0 0 0 0
Перепишем (3) в виде блочных матриц:
ГГтт I"тЬ
-ГЪт Г^
_ [ ^тт ІтЬІ л _ [Лтт 0 1 ; т — [ттт 0]
[Гьт Гь ь]; [ 0 0 [ 0 0 1
(4)
Соответственно, перемещения наделённых массой узлов запишем в матрицу-столбец Ут; перемещения безмассовых узлов, находящихся на границе тела, запишем в матрицу-столбец {7^. Для одного тела /, ( / = 1 ,2 по рисунку 2), получим:
№ о-
Гтт _ [ 25і] , ГтЬ _ ГЬ т _ ["5і "5і] ; ГЬ Ь _
Лтт _ [ 2Лі]; М тт _ [Мі] , и —
и,
и
ъ.
Ю -і. ; ит — [ У] ; и —
и 1-і
У+1
(5)
(6)
Для системы двух тел 1 и 2 (Рис. 2):
Гтт _
251 0
0 25^
;
тЪ
—
1Ът
" 1 " 1 00
0
"5 2
0
" 2
Гь ь _
1 0 0 0
0 1 0 0
0 0 2 0
0 0 0 2
(7)
Лтт _
и —
ъ.
и?ТЇ
2
; —
1
Уз
У4
Уб
(8)
Аналогично, для системы п тел запишем:
Гтт —
■251 0 0 " 1 1 1 " 0 0 0 0
0 2-2 ■ 0 7 0 0 " 2 " 2 0 0
; ГтЬ ГЬт
0 0 • 2 5П. _ 0 0 0 0 " "5
(9)
Гь b =
тН 0 0 0 0 0 і
0 і 0 0 0 0 Уэ
0 0 2 0 0 0 i/4
0 0 0 2 0 0 ; U = //б
0 0 0 0 •• S lJn 0 У3п - 2
_0 0 0 0 0 s lJnJ - У3п -
(10)
тН < 0 0 тН § 0 0 2
= 0 2 0 ; M mm = 0 2 0 ; Um = /5
0 0 Лn- 0 0 Mn. э і
U =
U,
u,
b .
(її)
Соответственно, вес каждого тела укажем в матрице-столбце Fm:
Г^і'
(12)
Здесь Fj = Mjg , где д = 9, 8 1 м / с2.
Проекции сил контактного взаимодействия N на оси неподвижной декартовой системы координат запишем в матрицу-столбец Fb = С bN . Здесь N - матрица-столбец сил контактного взаимодействия; С ь - матрица направляющих косинусов, в которой число столбцов к равно числу зазоров, а число строк равно числу указанных выше безмассовых узлов:
=
-1 0 0 0
0 - 1 0 0
0 1 0 0
0 0 - 0
0 0 1 • 0
1 0 0 0 - -
N =
і
2
N,
(13)
Для рассматриваемой конечно-элементной модели уравнение движения можно записать в часто используемой векторно-матричной форме [33,34]:
Здесь
rU + MJ+m U = F.
=
(14)
(15)
Необходимо подчеркнуть, что в стандартном уравнении вида (14) неизвестными являются только элементы вектора и [33,34]. В рассматриваемом случае неизвестными являются элементы векторов и и N. Поэтому, в дополнение к (14), необходимы соотношения между
2
силами контактного взаимодействия N и зазорами D, что усложняет решение проблемы. Кроме того, затруднения при решении задачи появляются по той причине, что, как отмечено
выше, перемещения, зазоры и силы контактного взаимодействия при соударениях тел друг с
другом и с корпусом барабана изменяются скачкообразно; перемещения Um и U й, силы контактного взаимодействия N и зазоры D не являются гладкими функциями времени.
Соотношение между перемещениями узлов U и зазорами D запишем, принимая во внимание геометрические аспекты задачи:
D = С £Ub + D о . (16)
Здесь D о - матрица-столбец начальных зазоров.
Принимая во внимание структуру матриц г, À, m (4), запишем соотношение (14) в виде системы двух равенств:
rmmUm AmmUjn + MmmUm + rmbUb = Fm, (17)
гЬ m Um гЬ bUb С bN . (18)
Поскольку взаимопроникновения тел при соударениях нет, то зазоры (16) не могут быть отрицательными:
D = C[Ub+ Dо > 0. (19)
По физическому смыслу задачи, как пояснено выше,
D > 0; N > 0 ; NTD = 0. (20)
Таким образом, формализованная задача сводится к решению системы равенств и неравенств (17)-(20).
Систему равенств и неравенств (17)-(20) решим численно, применяя метод конечных разностей с шагом по времени т. Для этого скорость Um и ускорение Um на шаге i (т.е., при tj = /т) выразим через перемещения Um на шагах i, i- 1 и i-2 , используя следующие конечноразностные соотношения первого порядка точности:
.Ц (j _^ -i- Ц (j) Ц (j _ 2)-9ÏT (j _^ -1- И (j)
Vj(i) _ ~um “H um _ жт(0 _ m ^ um /оi \
m = т ; m = т2 ' ( )
Конечно-разностная аппроксимация большей точности требует использования перемещений U^ 3) на шаге ( i- 3 ) и т. д., что для негладкой функции может оказаться некорректным. Более детально этот вопрос рассмотрен, например, в [35].
Применение центральных конечно-разностных схем целесообразно при N=0, что имеет
место в фазе бесконтактного движения. Если имеет место соударение тел, то N ^ 0 и поэтому
на шаге в случае применения центральных конечных разностей конечно-разностный аналог системы равенств и неравенств (17)-(20), записанный для момента времени t j = i т, будет содержать неизвестные U(j+1) и N (j) для разных моментов времени. Препятствием к использованию величины N(j+1) взамен N ( j ) служит то обстоятельство, что при tj+1 и tj
могут быть неодинаковыми структура механической системы, фаза движения, а значит и структура матриц коэффициентов в соотношениях (17)-(20).
Резюмируя, заметим, что применение односторонних конечно-разностных схем (21) для численного решения системы равенств и неравенств (17)-(20) позволяет построить алгоритм, который не требует заранее знать, в контактной или в бесконтактной фазе движется любое из рассматриваемого множества тел. При этом контактная и бесконтактная фазы движения моделируются одной и той же системой равенств и неравенств (17)-(20). Предлагаемый подход позволяет моделировать переход односторонних связей из возможного состояния в действительное. Реализуя этот подход, подставим (21) в (17). Получим после преобразований:
Тогда, с учетом (22) и (23), взамен аналитических соотношений (17)-(20) получим их дискретные аналоги для момента времени t ^ = / т:
Обозначим:
(24)
(25)
(26)
(27)
Из (25) найдём:
(28)
(29)
Подставим (25) в (28) и полученное при этом выражение и "подставим в (26). Получим:
(0
(30)
(31)
D С ^ = AN(i ) + Bc (32)
D © > 0 ; N © > 0 ; D ©N © T = 0 . (33)
Соотношения (32) и (33) моделируют все возможные состояния рассматриваемого объекта исследования на шаге i. Матрица-столбец Bc i ) моделирует внешнее воздействие на объект; реакцию объекта на внешнее воздействие моделируют матрицы-столбцы зазоров D С i -* и сил контактного взаимодействия N © ; структуру объекта моделирует ( п Хп )-матрица A ;
структура объекта зависит также от зазоров, которые определены элементами матрицы-
столбца .
Система равенств и неравенств (32) и (33) представляет собой линейную задачу дополнительности (далее по тексту ЛЗД), которая появляется во многих приложениях [1,37] и для которой общепринято обозначение LCP (Linear Complementarity Problem) [36]. В задачах рассматриваемого класса матрица коэффициентов A в (32) положительно определенная. Алгоритмы решения линейной задачи дополнительности исследованы в [37,38 и др.]. Для решения часто используется хорошо известный метод Lemke [39]. Однако в задачах рассматриваемого класса этот метод требует большого объема вычислений, на что обращено внимание, например, в статье [27]. Поэтому продолжаются поиски более экономичных алгоритмов решения линейной задачи дополнительности [38,40-42]. Анализ литературы показал, что остаются недостаточно исследованными перспективы применения декомпозиции и исключений Гаусса-Жордана при решении рассматриваемой задачи.
3. Результаты
3.1. Обоснование алгоритма декомпозиции линейной задачи дополнительности
В ближайшем изложении, используя работы [23, 24, 25], рассмотрим новое
формализованное обоснование алгоритма декомпозиции линейной задачи дополнительности
(32), (33) с положительно определенной матрицей коэффициентов A. Для этого случая доказаны теоремы существования и единственности решения [37,43-45]. Однако продолжаются поиски алгоритмов решения, наиболее эффективных в вычислительном отношении [46,47].
В целях компактности дальнейшего изложения воспользуемся общепринятыми обозначениями [43,45], чтобы записать соотношения (32) и (33) в эквивалентном виде одной из возможных формулировок ЛЗД: для положительно определенной (п Х п) матрицы (М) и ( ) вектора найти векторы и такие, что
w = Mz + q; (34)
w > 0; z > 0; wTz = 0 .
(35)
В силу положительной определенности матрицы М решение задачи (34), (35) для любого q Є Мп существует и единственно [37,43].
Если q > 0, то ш = q; z=0.
Рассмотрим случай, когда вектор q содержит отрицательные элементы.
В силу неотрицательности ш и г решение задачи (34), (35) не изменится, если условие ш7г = 0 заменить требованием г ¿ш £ = 0 [43]:
г^ = г(Маг1 + М£2г2 + ■ ■ ■ + М^гп + q¿) = 0, і = 1, 2,. . . ,п. (36)
В данном случае гі = 0, і = 1 ,2 ,... , п не отвечает условию ш > 0. Следовательно, равны нулю выражения в скобках (36) и для задачи (34), (35) необходимо решить систему следующих уравнений и неравенств:
ш = Мг + q; (37)
ш > 0; г > 0 . (38)
Если Му = Му і (аналогично (32), (33)), то задача (37), (38) эквивалентна следующей задаче поиска координат точки условного минимума
Е = т і п(^г7Мг + г7 q); г > 0 . (39)
Рассмотрим возможные варианты эволюции объекта, для которого сформулирована
задача (39). Моделируемым объектом может быть некоторый процесс или механическая система. Одна из таких систем рассмотрена выше (рис. 2).
Ситуация 1. Предположим, что моделируемый объект принудительно переведен в возможное состояние, такое, что г = 0, т.е. что все гі = 0, і = 1, 2 .... ,п. Тогда по (34) ш = q. При этом, если 0 , то данное состояние объекта является не только возможным, но и действительным. В этом случае объект, предоставленный самому себе, останется в данном состоянии: = 0, = .
Ситуация 2. Пусть объект принудительно переведен в возможное состояние, такое, что = 0. Если при этом некоторые элементов вектора отрицательны, то вектор = не отвечает условию ш >0 (35). Предоставленный самому себе, объект будет
трансформироваться из возможного состояния в действительное. Выполним декомпозицию, разбив процесс трансформации на отдельные стадии. Каждой стадии соответствует один шаг алгоритма.
Пусть на старте шага 1 все г^. = 0, к = 1,2 .... ,п. Пусть на финише шага 1 в векторе г только один элемент гі принимает новое значение, такое, что гі > 0, 1 < і < п. Номер элемента і = і * заранее неизвестен. Значение і* определим, используя (39). Примем во внимание, что, согласно (35), ші = 0; г^. = 0, к = 1, 2....,п, к Ф і.
Тогда Е (39) является функцией одной переменной г ¿:
Е (г) = ті п(^М і ¿г2 + q ¿г) ; г > 0 . (40)
В стационарной точке производная функции (40) равна нулю:
MaZi + ^ = 0 . (41)
Из (41) определим Zi = — q i/Mi i и найдем значение функции (40) в стационарной точке:
q2
f?=— т (42)
В стационарной точке имеет место минимум функции (40), т.к. вторая производная равна Mi i >0 (в задачах рассматриваемого класса матрица M положительно определена). Таким образом, на финише рассматриваемого шага значение функции (40) может только уменьшаться на величину (42), что обеспечивает сходимость метода.
Напомним, как указано выше, на старте шага q i < 0. Если в векторе q больше одного
отрицательного элемента, то для каждого q i < 0 определяется значение локального
2
минимума: f? = — ^~; i = 1,2 ,...,n; q i < 0 . Наименьшее из полученных значений
(глобальный минимум)
f • = m i n ( — ; i = 1, 2,..., n; qi < 0 (43)
соответствует указанному выше индексу = того элемента , который на финише шага принимает новое значение 0, .
Если для различных значений получены одинаковые значения , то используется любое из них (такой случай тестируется ниже).
Нестрогое неравенство 0 означает, что если на старте шага = 0, то на финише того же шага 0 и, следовательно, = 0 . Таким образом, в предлагаемом алгоритме не исключается из рассмотрения случай, когда на финише шага z i = 0 и w i = 0. Этот случай в некоторых алгоритмах затрудняет поиск решения [37].
Очевидно, для выбора указанного выше i = i ? целесообразно использовать взамен (43) формулу
2 ^
р?=— I • (44)
Соответственно, задачу (43) можно записать в эквивалентной форме:
р • = mi n ( — ; i = 1, 2,..., n; q i < 0 (45)
Формула (45) есть не что иное, как критерий определения тех переменных Zi и Wi, которые при данном воздействии переходят из возможного состояния в действительное. При этом действительное состояние моделируемого объекта определяется как состояние, в котором достигается глобальный минимум дискретного множества функций (40).
Дальнейшие вычисления сводятся к преобразованию равенства (34) по схеме одного шага исключений Гаусса-Жордана с разрешающим элементом М£ ¿. Исключения Гаусса-Жордана представляют собой удобную схему реализации тождественных алгебраических
преобразований; в данном случае выполняются преобразования системы уравнений (34).
Применение исключений Гаусса-Жордана в алгоритме моделирования механических систем с односторонними связями рассмотрено в [23]. Другие приложения исследованы в [49].
3.2. Алгоритм решения линейной задачи дополнительности
Один шаг исключений Гаусса-Жордана с разрешающим элементом М£ ^ выполняется по следующей схеме. Из уравнения / системы (34)
и, = Му^ + + ■ ■ ■ + Му„г„ + <7у, у = 1,2 ,. . . , п (46)
выражается г£ и подставляется в остальные уравнения той же системы. Система
уравнений (34) может быть представлена в виде следующей таблицы:
(47)
Из уравнения / выразим г£ и подставим в остальные уравнения. Получим:
г1 ■ • Щ ■ ч
1 МП-М^М1/МН ■ ■ Ми/Ма ' М1„-М1 ¿М1„/М1 1 ql - М1 ¿q ¡/М* 1
-М^/Мн ■ ■ 1 /Ма • 'Мт/Мц ~Ч1/М и
Пп М„1- М„ ¿М11/ М1 1 ■ " Мп1/Мц • " мпп - Мш М1п / М ц Чп~МтЧ1/Ми
(48)
С учетом изложенного выше конкретизируем предлагаемый алгоритм решения ЛЗД.
Шаг 1. Для всех q£ < 0 вычислить критерий /* (45). С использованием данного критерия определяется разрешающий элемент М£ £ для одного шага исключений Гаусса-Жордана (48). Если для различных значений / получены одинаковые значения критерия (45), то используется любое из них.
Шаг 2. Если преобразованный вектор q* (48) содержит отрицательные элементы, то повторяется шаг 1. Если в преобразованном векторе q* все элементы неотрицательны, то преобразованные векторы и и г являются решением задачи. Иначе выполняется переход к шагу 1.
Резюмируя, отметим следующее. Предлагаемый алгоритм является шаговым алгоритмом. На каждом шаге решается выполняются преобразования (41). Таким образом, исходная задача линейной дополнительности (34), (35) с (п X п )-матрицей М редуцирована к последовательности задач определения стационарных значений квадратичной функции с одной переменной (40). Ключевую роль в предлагаемом подходе играет критерий (45), который можно назвать энергетическим критерием выбора разрешающего элемента в матрице коэффициентов линейной задачи дополнительности. Правомерность термина «энергетический критерий» обсуждается ниже.
3.3. Ограничения на область применения алгоритма
Данный алгоритм разработан для решения ЛЗД (34), (35) с симметричной положительно определенной матрицей коэффициентов . Заметим, что к задачам такого класса сводятся многочисленные задачи прикладной механики [22,23,26,27,41], в том числе численное моделирование соударений балансов в корообдирочном барабане [5,31]. Это обстоятельство определяет достаточно большую область применения предлагаемого алгоритма, а также актуальность и практическую значимость выполненного исследования.
3.3. Тестовые примеры
Пример 1. Найти векторы и и г, такие, что и = Мг + q; и > 0, г > 0, итг = 0, где
"4 - 0 0 -
- 4 - 0
= , q =
0 - 4 - -4
_0 0 - 4. -2-
Известное решение = = 0 0 получено в статье [46] по итерационному
алгоритму, который базируется на эквивалентности ЛЗД и нелинейного уравнения ( = 0 При этом потребовалось три итерации.
Рассмотрим применение предложенного нами алгоритма. Исходные данные и результаты вычислений представим в табличной форме по аналогии с (47) и (48). В отдельном столбце приведем значения р * (45).
Принимаем р* = р1 = - 4 . Выполним один шаг исключений (48) с разрешающим элементом 11 =
Находим р* = р3 = - 4 . Выполним один шаг исключений (48) с разрешающим элементом М33 = 4 :
Получаем: и1 = 0; г2 = 0; и3 = 0; г4 = 0; г1 = 1; и2 = 1 ; г3 = 1; и4 = 1.
Таким образом, = 0 0 ; = 0 0 Для решения задачи
потребовалось два шага исключений Гаусса-Жордана с выбором разрешающего элемента по критерию (45).
Пример 2. Найти векторы и , такие, что = ; 0 0 = 0 где
1 00 - - - -1-
- 0 - - -
- - 100 - , q = 3
_ - 4 - - 00. -
42
Известное решение г* = 10 — 0 —I опубликовано в статье [40]. В данной статье
линейная задача дополнительности сформулирована как задача квадратичного программирования и предложен итерационный метод решения сформулированной задачи. Решение г* получено в [40] за 54 итерации.
Рассмотрим применение предложенного нами алгоритма с применением энергетического критерия (45). Исходные данные и результаты вычислений представим в табличной форме по аналогии с (47) и (48). В отдельном столбце приведем значения критерия р* (45). Незначащие нули далее не указываем.
Принимаем р* = р2 = -0,0 8 . Выполним один шаг исключений (48) с разрешающим элементом М22 = 5 0 :
г1 И2 гз г4 Ч * VI
99,92 - 0,04 -3,24 -4,28 0,9 2
2 0,04 0 0 0 0 0,04
Из -3,24 0 99,28 -11,84 2,76
4 -4,28 - 0, 14 -11,84 199,0 2 -4,28 -0,09204
Выполним один шаг исключений (44) с разрешающим элементом М44 = 199,02. Получим:
г! И2 гз 4 Ч * VI
1 99,82796 -0,043 0 1 1 -3,49462 -0,0 2 1 5 1 0,82 79 5 7
2 0,043 0 1 1 0 0 00 0 0 000 0 0,0430 1 1
Из -3,49462 0 99,5762 0 0 2, 50 53 76
г4 0,02 1 50 5 -0,000 70 3 0,059492 0,00502 5 0,02 150 5
Таким образом, г1 = 0; и2 = 0; г3 = 0; и4 = 0;
IV! = 0,82 795 7; г2 = 0,0430 1 1; и3 = 2, 50 53 76; г4 = 0,02 150 5. г = [0 0,0430 1 1 0 0,02 1505]т; и/ = [0,82 795 7 0 2, 50 53 76 0] т .
Для решения задачи потребовалось два шага исключений Гаусса-Жордана с выбором разрешающего элемента по критерию (45).
Пример 3. Найти векторы и и г, такие, что и = Мг + 4; и > 0, г > 0, итг = 0, где
"4 - 2 0 -
- 4 - 0 1
= , Ч = -1.
2 - 4 - -2
_0 0 - 4. 0
Опубликованных решений данной задачи обнаружить не удалось. Особенность решения задачи в том, что, как показано ниже, три пары переменных г^ ( / = 2 , 3,4) равны нулю.
Такой случай затруднителен для некоторых известных алгоритмов [37].
Рассмотрим применение предложенного нами алгоритма. Исходные данные и результаты вычислений представим в табличной форме по аналогии с (47) и (48). В отдельном столбце приведем значения критерия (45).
По (45) находим: p* = pi = - 4 . Выполним один шаг исключений (45) с разрешающим элементом Мп = 4 :
Получаем: и1 = 0; г2 = 0; г3 = 0; г4 = 0; г1 = 1 ; и2 = 0; и3 = 0; г4 = 0.
Таким образом, = 0 0 0 ; = 0 0 0 0
Для решения задачи потребовался один шаг исключений Гаусса-Жордана с выбором разрешающего элемента по критерию (45).
Обсуждение и заключение. В данной работе рассмотрено применение предлагаемого алгоритма в конечно-элементной модели соударений балансов в корообдирочном барабане. Простейший случай показан на рисунке 2. Результаты применения алгоритма для моделирования более сложных технологических ситуаций рассмотрены в статье [5].
В данной работе приведено новое обоснование критерия (45). При этом отличительная особенность заключается в том, что для обоснования данного критерия использована только формулировка линейной задачи дополнительности (34), (35).
Аналогичный критерий, анонсированный в сообщении [24], получен в работе [25]. Однако для его обоснования был использован принцип минимума потенциальной энергии упругой механической системы. Этот принцип находится вне рамок формулировки линейной задачи дополнительности (34), (35). Однако задача (40) эквивалентна задаче поиска условного минимума полной потенциальной энергии.
Поскольку критерий (45) получен как результат поиска глобального минимума, то представляется возможным утверждение о предельной эффективности (простоте) предложенного алгоритма решения линейной задачи дополнительности (34), (35) с положительно определенной матрицей коэффициентов. Это означает, что, предположительно, число шагов (или итераций) в рассмотренных выше тестовых примерах невозможно уменьшить. Правомерность этого предположения позволит проверить публикация полученных результатов в журнале открытого доступа.
Применение полученных результатов для численного моделирования взаимодействия еловых балансов неодинакового диаметра в корообдирочном барабане достаточно подробно рассмотрено в статье [5].
Перспективы дальнейших исследований в затронутом направлении численного моделирования технологических процессов могут быть связаны с прогнозированием степени очистки в зависимости от числа соударений и сил контактного взаимодействия балансов друг с другом и с корпусом барабана. Оценки числа соударений в начальной стадии обрушения массива балансов могут быть определены по результатам численного
моделирования [5]. Получение и анализ подобных данных в целях прогнозирования степени очистки с учетом характеристик древесины балансов (ель, сосна и др.) требует продолжения исследований с учетом известных в данной области результатов [15,52-59].
Представленный выше алгоритм решения линейной задачи дополнительности с положительно определенной матрицей коэффициентов, очевидно, может быть использован при моделировании контактного взаимодействия в двух- и трехмерных системах многих тел. При этом задача идентификации контактных пар может составить предмет отдельного исследования.
Благодарности
Работа выполнена в рамках Программы стратегического развития Петрозаводского государственного университета на 2012-2016 годы по госбюджетной теме 126-12 Министерства образования и науки РФ.
Литература
1. Ferris M. C., Pang J. S. Engineering and economic applications of complementarity problems // Siam Review. - 1997. - Т. 39. - №. 4. - С. 669-713.
2. Glocker C., Studer C. Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics // Multibody System Dynamics. - 2005. - Т. 13. - №. 4.
- С. 447-463.
3. Khulief Y.A. Modeling of impact in multibody systems: An overview // J. Comput. Nonlinear Dyn. - 2013. - Т. 8. - С. 021012.
http: //www .researchgate.net/publication/235708326_Modeling_of_Impact_in_Multibody_Systems_An_Ove rview/file/e0b49518983d093f38.pdf_[Cited 24 September 2013].
4. Линч А. Д. Циклы дробления и измельчения: моделирование, оптимизация,
проектирование и управление: Пер. с англ. - М.: Недра, 1981. - 343 с.
5. Васильев С. Б., Доспехова Н. А., Колесников Г. Н. Численное моделирование взаимодействия еловых балансов неодинакового диаметра в корообдирочном барабане // Resources and Technology. - 2013. - Т. 10. - № 1. - С. 024-038.
6. Бойков С. П. Теория процессов очистки древесины от коры. - Л.: Издательство Ленинградского университета, 1980. - 152 с.
7. Baroth R. Literature review of the latest development of wood debarking. - University of Oulu, Control Engineering Laboratory, 2005. - 29 pp..
8. Oman M. Influence of log characteristics on drum debarking of pulpwood //Scandinavian Journal of Forest Research. - 2000. - Т. 15. - №. 4. - С. 455-463.
9. Газизов А. М., Шапиро В. Я., Григорьев И. В. Моделирование процесса разрушения коры при роторной окорке древесины // Вестник Красноярского государственного аграрного университета. - 2008. - № 5. - С. 271-279.
10. Гаспарян Г. Д. Основы метода и технологии ультразвуковой окорки круглых лесоматериалов // Фундаментальные исследования. - 2013. - № 6 (часть 1). - С. 19-23.
11. Оскерко В. Е. Новый принцип окорки лесоматериалов // Строительные и дорожные машины. - 2007. - № 3. - С. 13-16.
12. Васильев С. Б., Колесников Г. Н., Никонова Ю. В., Раковская М. И. Влияние локальной жесткости корпуса корообдирочного барабана на изменение силы соударений и величину потерь древесины // Ученые записки Петрозаводского государственного университета. Серия: Естественные и технические науки. - 2008. - № 96. - С. 84-91.
13. Васильев С. Б., Колесников Г. Н., Никонова Ю. В., Раковская М. И. Исследование закономерностей изменения силы соударений с целью снижения потерь при окорке древесины в барабане // Известия Санкт-Петербургской лесотехнической академии. -2008. - № 185. - С. 195-202.
14. Девятникова Л. А. Потенциал ресурсосбережения в технологии подготовки круглых лесоматериалов к переработке на щепу // Политематический сетевой электронный научный журнал Кубанского государственного аграрного университета. - 2013. - № 88.
- С. 188-206.
15. Никонова Ю. В. Обоснование конструктивно-технологических параметров корообдирочных барабанов с применением численного моделирования динамического взаимодействия балансов: автореф. дисс. ... канд. техн. наук - Петрозаводск, 2009. - 20 с.
16. Isokangas, A., Hyvonen, A., Pollanen, K., Leiviska, K. Analysis of debarking times using pilot scale drum. In 59th Appita Annual Conference and Exhibition: Incorporating the 13th ISWFPC (International Symposium on Wood, Fibre and Pulping Chemistry), Auckland, New Zealand, 16-19 May 2005: Proceedings (p. 95). Appita Inc.
17. Аргатов И. И., Дмитриев Н. Н. Основы теории упругого дискретного контакта. - Санкт-Петербург: Политехника, 2003. - 233 с.
18. Маркеев А. П. Динамика твердого тела при наличии его соударений с твердой поверхностью // Нелинейная динамика. - 2008. - Т. 4. - № 1. - С. 1-38.
19. Бураго Н. Г., Кукуджанов В. Н. Обзор контактных алгоритмов // Известия Российской академии наук. Механика твердого тела. - 2005. - № 1. - С. 45-87.
20. Dobias, J., S. Ptak, Z. Dostal, D. Horak, R. Kucera and V. Vondrak. Semicoercive contact
problems with large displacements by FETI domain decomposition method // European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2004. P. Neittaanmaki, T. Rossi, K. Majava, and O. Pironneau (eds.), R. Owen and M. Mikkola (assoc. eds.). Jyvaskyla, 24-28 July 2004. - P. 1-12.
http://eng.imamod.ru/~serge/arc/conf/ECC0MAS_2004/ECC0MAS_V1/proceedings/pdf/6.pdf [Cited 24 September 2013].
21. Li Y. et al. A contact analysis approach based on linear complementarity formulation using smoothed finite element methods // Engineering Analysis with Boundary Elements. - 2013. -Т. 37. - №. 10. - С. 1244-1258.
22. Pfeiffer, F. On non-smooth multibody dynamics. Proceedings of the Institution of Mechanical Engineers, Part K. // Journal of Multi-body Dynamics. -2012. - Т. 226(2). - C. 147-177.
23. Колесников Г. Н. Дискретные модели механических и биомеханических систем с
односторонними связями. - Петрозаводск: Петрозаводский государственный
университет, 2004. - 202 с.
24. Колесников Г. Н., Раковская М. И. Энергетический критерий очередности перехода односторонних связей в действительное состояние // Обозрение прикладной и промышленной математики. - 2006. - Т. 13. - С. 652.
25. Колесников Г. Н., Раковская М. И. Энергетический критерий очередности перехода односторонних связей в действительное состояние // Строительная физика в XXI веке: Материалы научно-технической конференции / Под. Ред. И. Л. Шубина. - М.: НИИСФ РААСН. - 2006. - С. 606-608.
26. Ким Т. С., Яцура В. Г. Расчет систем с односторонними связями как задача о дополнительности // Строительная механика и расчет сооружений. - 1989. - № 3. - С. 41-
43.
27. Pfeiffer F., Glocker C. (2000). Contacts in multi-body systems. // Journal of Applied Mathematics and Mechanics. - 2000. - Т. 64(5). - С. 773-782.
28. Zhu S. Q., Hao F., Gao H. T., Han Y. L. Numerical Simulation for Non-Smooth Contact Problem in Mechanical Systems //Applied Mechanics and Materials. - 2013. - Т. 341. - С. 491-495.
29. Galvanetto U., Colombo A. (2013, January). Computation of the Basins of Attraction in Nonsmooth Dynamical Systems. // IUTAM Symposium on Nonlinear Dynamics for Advanced Technologies and Engineering Design (pp. 17-29). - Springer Netherlands.
30. Krabbenhoft K. et al. Granular contact dynamics using mathematical programming methods // Computers and Geotechnics. - 2012. - Т. 43. - С. 165-176.
31. Никонова Ю. В., Раковская М. И. Методика определения жесткости балансов, результаты численных экспериментов и испытаний образцов // Resources and Technology.
- 2010. - № 8. - С. 100-106.
32. Фиалко С. Ю. Прямые методы решения систем линейных уравнений в современных МКЭ комплексах. - М.: Издательство СКАД СОФТ. - 2009. - 161 с.
33. Бидерман В. Л. Теория механических колебаний: Учебник для вузов. - М.: Высшая школа, 1980. - 408 с.
34. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов: Пер. с англ.
- М.: Стройиздат, 1982. - 447 с.
35. Зылев В. Б. Вычислительные методы в нелинейной механике конструкций. - М.: НИЦ «Инженер», 1999. - 144 с.
36. Cottle R. W., Dantzig G. B. Complementary pivot theory of mathematical programming // Linear algebra and its applications. - 1968. - Т. 1. - №. 1. - С. 103-125.
37. Murty, K. G. (1988). Linear complementarity, linear and nonlinear programming. - Berlin: Heldermann Verlag. - 629 pp.
38. Csizmadia Z., Illes T. New criss-cross type algorithms for linear complementarity problems with sufficient matrices // Optimization Methods and Software. - 2006. - Т. 21. - №. 2. - С. 247-266.
39. Lemke C. E. Bimatrix equilibrium points and mathematical programming // Management science. - 1965. - Т. 11. - №. 7. - С. 681-689.
40. Elfoutayeni Y., Khaladi M. Using vector divisions in solving the linear complementarity problem // Journal of Computational and Applied Mathematics. - 2012. - Т. 236. - №. 7. - С. 1919-1925. arXiv:1005.1417v1 [cs.RO] 9 May 2010 [Cited 24 September 2013].
41. Ловцов А. Д. Алгоритмы линейной задачи дополнительности в применении к расчету систем с односторонними связями // Тезисы докладов XX Междунар. конф. «Математическое моделирование в механике сплошных сред. Метод граничных и конечных элементов», 24-26 сент. 2003. СПб.- 2003. - С. 128-129.
42. Kheirfam, B. 2013. A full-Newton step infeasible interior-point algorithm for linear complementarity problems based on a kernel function // Algorithmic Operations Research. - Т. 7(2). - 103-110.
43. Попов Л. Д. Введение в теорию, методы и экономические приложения задач о
дополнительности: Учебное пособие. - Екатеринбург: Издательство Уральского
государственного университета, 2001. - 124 с.
44. Cottle R. W., Dantzig G. B. A generalization of the linear complementarity problem // Journal of Combinatorial Theory. - 1970. - Т. 8. - №. 1. - С. 79-90.
45. Втюрина, М. В., Жадан, В. Г. Барьерно-проективный метод с наискорейшим спуском для линейных задач дополнительности // Журнал вычислительной математики и математической физики. - 2005. - Том 45. - № 5. - С. 792-812.
46. Elfoutayeni Y., Khaladi M. A New Interior Point Method for Linear Complementarity Problem // Applied Mathematical Sciences. - 2010. - Т. 4. - №. 66. - С. 3289-3306. http://www.m-hikari.com/ams/ams-2010/ams-65-68-2010/elfoutayeniAMS65-68-2010.pdf [Cited 24 September 2013].
47. Долгачев М. В. Обзор алгоритмов линейной задачи дополнительности в применении к контактным задачам // Новые идеи нового века: материалы международной научной конференции ФАД ТОГУ = The new Ideas of New Century: The International Scientific Conference Proceedings of FAD PNU. - 2006. - Т. 1. - С. 280-284.
48. Бахвалов, Н. С. Численные методы: Учебное пособие. - М.: Наука, 1993. - 631 с.
49. Yang K., Li Y., Xia Y. A Parallel Method for Matrix Inversion Based on Gauss-Jordan Algorithm // Journal of Computational Information Systems. - 2013. -. 9. №. 14. - С. 55615567.
50. Tasora, A., Anitescu, M. A complementarity-based rolling friction model for rigid contacts // Meccanica. - 2013. - Т. 48(7). - С. 1643-1659.
51. Duchesne I., Nylinder M. Measurement of the bark/wood shear strength: practical methods to evaluate debarking resistance of Norway spruce and Scots pine pulpwood // Forest products journal. - 1996. - Т. 46. - № 11/12.- С. 57-62.
52. Hatton J. V. Debarking of frozen wood // Tappi journal. - 1987. - T. 70(2) - C. 61-66.
53. Isokangas, A., Leiviska, K. Optimization of wood losses in log debarking drum. Paperi ja puu.
2005. - Т. 87(5). - С. 324-328.
54. Karha K., Jylha P., Laitila J. Integrated procurement of pulpwood and energy wood from early thinnings using whole-tree bundling // Biomass and Bioenergy. - 2011. - Т. 35. - № 8. - С. 3389-3396.
55. Lauri P., Kallio M., Schneider U. A. The future development of the use of wood in Russia and its potential impacts on the EU forest sector // Scandinavian Journal of Forest Research. - 2013.
- Т. 28. - №. 3. - С. 291-302.
56. Piggott R. R., Thompson R. A. Drum debarking: key factors for design and performance // Tappi Journal. - 1987. - Т. 70. - №. 8. - С. 37-41.
57. Tamarozzi T., Ziegler P., Eberhard P., Desmet W. Static modes switching in gear contact
simulation // Mechanism and Machine Theory. - 2013. - Т. 63. - С. 89-106.
58. Шегельман И. Р., Колесников Г. Н., Васильев А. С., Никонова Ю. В. Моделирование технологического процесса очистки древесины в корообдирочном барабане с применением метода дискретных элементов // Известия Санкт-Петербургской лесотехнической академии. - 2008. - № 184. - С. 172-179.
59. Мазуркевич Е. О., Петрова Е. Г., Стрекаловский А. С. О численном решении линейной задачи дополнительности // Журнал вычислительной математики и математической физики. - 2009. - Том 49. - № 8. - С. 1385-1398.
References
1. Ferris, M. C., Pang, J. S. 1997.Engineering and economic applications of complementarity problems // Siam Review 39(4): 669-713.
2. Glocker, C., Studer, C. 2005. Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics. Multibody System Dynamics 13(4): 447-463.
3. Khulief, Y. A. 2013. Modeling of impact in multibody systems: An overview. J. Comput. Nonlinear Dyn. 8: 021012.
http://www.researchgate.net/publication/235708326_Modeling_of_Impact_in_Multibody_Syst ems_An_0verview/file/e0b49518983d093f38.pdf_[Cited 24 September 2013].
4. Lynch, A. D. 1981. Crushing and grinding cycles: modeling, optimization, design and management. (In Russian). Moscow, Russia, 343 pp.
5. Vasilyev, S. B., Dospehova, N. A., Kolesnikov G. N. 2013. Simulation of Unequal Diameter Spruce Pulpwood Interaction in Debarking Drum. (In Russian). Resources and Technology 10: 024-038.
6. Boykov, S. P. 1980. The theory of processes of wood debarking. (In Russian). Leningrad State University, Leningrad, USSR, 152 pp.
7. Baroth, R. 2005. Literature review of the latest development of wood debarking. Control Engineering Laboratory University of Oulu, Oulu, Finland, Report A No 27, 29 pp.
8. Oman, M. 2000. Influence of log characteristics on drum debarking of pulpwood. Scandinavian Journal of Forest Research 15(4): 455-463.
9. Gazizov, A. M., Shapiro, V. Ja., Grigoriev, I. V. 2008. Modeling of bark destruction at ring debarking. (In Russian). Vestnik Krasnoyarsk State Agrarian University 5: 271-279.
10. Gasparyan, G. 2013. Basis of the method and technology of ultrasonic bark round wood. (In Russian). Fundamentalnye issledovaniya 6(1): 19-23.
11. Oskerko, V. E. 2007. New debarking principle of wood. (In Russian) Construction and Road Building Machinery 3: 13-16.
12. Vasilyev, S. B., Kolesnikov, G. N., Nikonova, Ju. V., Rakovskaya, M. I. 2008. Influence of Local Stiffness of the Drum on Collisions Forces and Size of Wood Losses. (In Russian). Proceedings of Petrozavodsk State University. Natural and Engineering Sciences 96: 84-91.
13. Vasilyev, S. B., Kolesnikov, G. N., Nikonova, Ju. V., Rakovskaya, M. I. 2008. Research in to change of force in pulpwood collisions in drum for the purpose of decreasing losses at debarking. (In Russian). Izvestiya Saint Petersburg State Forest Technical University 185: 195202.
14. Devyatnikova, L. A. 2013. A potential resource-saving technology of preparation of round forest products to chipping. (In Russian). Scientific Journal of the Kuban State Agrarian University 88(4): 188-206.
15. Nikonova, Ju. V. 2009. Justification of constructive and technological parameters of debarking drums using numerical modeling of dynamic interaction of pulpwood, PhD thesis abstract, Petrozavodsk State University, Petrozavodsk, Russia, 20 pp..
16. Isokangas, A., Hyvonen, A., Pollanen, K., Leiviska, K. 2005. Analysis of debarking times using pilot scale drum. Proceedings of 59th Appita Annual Conference and Exhibition Incorporating the 13th ISWFPC (International Symposium on Wood, Fibre and Pulping Chemistry), Auckland, New Zealand, 16-19 May 2005. Appita Inc., p. 95..
17. Argatov, I. I., Dmitriev, N. N. 2003. Fundamentals of the theory of discrete elastic contact (In Russian). Politehnica, St. Petersburg, Russia, 233 pp.
18. Markeev, A. P. 2008. Rigid body dynamics in the presence of his collision with a solid surface. (In Russian). Nonlinear Dynamics 4(1): 1-38.
19. Burago, N. G., Kukudzanov,V. N. 2005. Revue of contact algorithms. (In Russian). Mechanics of Solids 1: 45-87.
20. Dobiâs, J., Ptâk, S., Dostâl, Z., Horâk, D., Kucera, R., Vondrâk, V. 2004. Semicoercive contact problems with large displacements by FETI domain decomposition method. European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2004. P. Neittaanmaki, T. Rossi, K. Majava, and O. Pironneau (eds.), R. Owen and M. Mikkola (assoc. eds.). Jyvaskyla, 24-28 July 2004, pp. 1-12. http://eng.imamod.ru/~serge/arc/conf/ECCOMAS_2004/ECCOMAS_V1/proceedings/pdf76.pdf [Cited 24 September 2013].
21. Li, Y. et al. 2013. A contact analysis approach based on linear complementarity formulation using smoothed finite element methods. Engineering Analysis with Boundary Elements.-T. 37(10): 1244-1258.
22. Pfeiffer, F. 2012. On non-smooth multibody dynamics. Proceedings of the Institution of Mechanical Engineers, Part K. Journal of Multibody Dynamics 226(2): 147-177.
23. Kolesnikov, G. N. 2004. Discrete models of mechanical and biomechanical systems with unilateral constraints. (In Russian). Petrozavodsk State University, Petrozavodsk, Russia, 202 pp.
24. Kolesnikov, G. N., Rakovskaya, M. I. 2006. Energy criterion of priority transition unilateral constraints in the actual state. (In Russian). Review of Applied and Industrial Mathematics 13: 652.
25. Kolesnikov, G. N., Rakovskaya, M. I. 2006. Energy criterion of priority transition unilateral constraints in the actual state. (In Russian). Building physics in the XXI Century: Materials Science and Technology Conference / Under. Ed. I. L. Shubin. NIISF RAASN, Moscow, Russia, pp. 606-608.
26. Kim, T. S., Yatsura, V. G. 1989.Analysis of systems with unilateral constraints as the problem of complementarity. (In Russian). Stroitelnaya mehanika i raschet sooruzheny 3: 41-43.
27. Pfeiffer, F., Glocker, C. 2000. Contacts in multi-body systems. Journal of Applied Mathematics and Mechanics 64(5): 773-782.
28. Zhu, S. Q., Hao, F., Gao, H. T., Han, Y. L. 2013. Numerical Simulation for Non-Smooth Contact Problem in Mechanical Systems. Applied Mechanics and Materials 341: 491-495.
29. Galvanetto, U., Colombo, A. 2013. Computation of the Basins of Attraction in Non-smooth Dynamical Systems. In IUTAM Symposium on Nonlinear Dynamics for Advanced Technologies and Engineering Design (pp. 17-29). Springer Netherlands.
30. Krabbenhoft, K. et al. 2012. Granular contact dynamics using mathematical programming methods. Computers and Geotechnics 43:165-176.
31. Nikonova, Ju. V., Rakovskaya, M. I. 2010. Cross section stiffness of pulpwood: experiment and numerical modeling. (In Russian). Resources and Technology 8: 100-106..
32. Fialko, S. Y. 2009. Direct methods for solving systems of linear equations in modern FEM complexes. (In Russian). Moscow: Publishing SCAD SOFT, Russia, 161 pp.
33. Biderman, V. L. 1980. The theory of mechanical vibrations: Textbook (In Russian). Moscow, Russia, 408 pp.
34. Bathe, K. J., Wilson, E. L. 1982. Numerical methods of analysis and finite element method (In Russian). Stroyizdat, Moscow, USSR, 447 pp.
35. Zylev, V. B. 1999. Computational methods in nonlinear mechanics of structures. (In Russian). NIC "Inzhener", Moscow, Russia, 144 pp.
36. Cottle, R. W., Dantzig, G. B. 1968. Complementary pivot theory of mathematical programming. Linear algebra and its applications 1(1): 103-125.
37. Murty, K. G. 1988. Linear complementarity, linear and nonlinear programming. Berlin: Heldermann Verlag. 629 pp.
38. Csizmadia, Z., Illes, T. 2006. New criss-cross type algorithms for linear complementarity problems with sufficient matrices. Optimization Methods and Software 21(2): 247-266.
39. Lemke, C. E. 1965. Bimatrix equilibrium points and mathematical programming. Management science 11(7): 681-689.
40. Elfoutayeni, Y., Khaladi, M. 2012. Using vector divisions in solving the linear complementarity problem. Journal of Computational and Applied Mathematics 236(7): 19191925. arXiv:1005.1417v1 [cs.RO] 9 May 2010
41. Lovtsov, A. D. 2003. Algorithms for linear complementarity problems applied to the analysis of systems with unilateral constraints. (In Russian). In: Abstracts of the XX Int. Conf. "Mathematical modeling in continuum mechanics. The method of boundary and finite elements ", Sept. 24 -26. 2003. St. Petersburg, Russia, pp. 128-129.
42. Kheirfam, B. 2013. A full-Newton step infeasible interior-point algorithm for linear complementarity problems based on a kernel function. Algorithmic Operations Research 7(2): 103-110.
43. Popov, L. D. 2001.Introduction to the theory, methods and economic applications of complementarity problems: the manual. Ekaterinburg: Publishing, Ural State University, Russia, 124 pp.
44. Cottle, R. W., Dantzig, G. B. 1970. A generalization of the linear complementarity problem. Journal of Combinatorial Theory 8(1): 79-90.
45. Vtyurina, M. V., Zhadan, V. G. 2005. Barrier-projection method of quickest descent for linear complementarity problems. (In Russian). Computational Mathematics and Mathematical Physics 45(5): 792-812.
46. Elfoutayeni, Y., Khaladi, M. A. 2010. New Interior Point Method for Linear Complementarity
Problem. Applied Mathematical Sciences 4(66): 3289-3306. http://www.m-
hikari.com/ams/ams-2010/ams-65-68-2010/elfoutayeniAMS65-68-2010.pdf [Cited 24 September 2013].
47. Dolgacev, M. V. 2006. Overview of algorithms of linear complementarity problems applied to contact problems. (In Russian). The new Ideas of New Century: The International Scientific Conference Proceedings of FAD TOGU. Vol. 1: 280-284.
48. Bahvalov, N. S. 1993. Numerical methods: textbook (In Russian). Nauka, Moscow, Russia, 631 pp.
49. Yang K., Li Y., Xia Y. 2013. A Parallel Method for Matrix Inversion Based on Gauss-Jordan Algorithm. Journal of Computational Information Systems. 9(14): 5561-5567.
50. Tasora, A., Anitescu, M. 2013. A complementarity-based rolling friction model for rigid contacts. Meccanica 48(7), pp 1643-1659.
51. Duchesne, I., Nylinder, M. 1996. Measurement of the bark/wood shear strength: practical methods to evaluate debarking resistance of Norway spruce and Scots pine pulpwood. Forest products journal 46(11/12): 57-62
52. Hatton, J. V. 1987. Debarking of frozen wood. Tappi journal 70(2): 61-66.
53. Isokangas, A., Leiviska, K. 2005. Optimization of wood losses in log debarking drum. Paperi jaPuu 87(5): 324-328.
54. Karha, K., Jylha, P., Laitila, J. 2011. Integrated procurement of pulpwood and energy wood from early thinnings using whole-tree bundling. Biomass and Bioenergy 35(8): 3389-3396.
55. Lauri, P., Kallio, M., Schneider, U.A. 2013. The future development of the use of wood in Russia and its potential impacts on the EU forest sector. Scandinavian Journal of Forest Research 28(3): 291-302.
56. Piggott, R. R., Thompson, R. A. 1987. Drum debarking: key factors for design and performance. Tappi Journal 70(8): 37-41.
57. Tamarozzi, T., Ziegler, P., Eberhard, P., Desmet, W. 2013. Static modes switching in gear contact simulation. Mechanism and Machine Theory 63: 89-106.
58. Shegelman, I. R., Kolesnikov, G. N., Vasilyev A. S., Nikonova, Yu. V. 2008. Process simulation cleaning wood in debarking drum using discrete element method. Proceedings of the St. Petersburg Forest Technical Academy. 184: 172-179.
59. Mazurkevich, E. O., Petrova, E. G., & Strekalovsky, A. S. (2009). On the numerical solution of the linear complementarity problem. Computational Mathematics and Mathematical Physics, 49(8), 1318-1331.
© 2013 Ko^ecHHKOB r.H.
