CC BY
62
ВЕСТНИК ПНИПУ
2014 Механика № 2
УДК 539.3
Н.В. Пустовой, В.Е. Левин, Д.А. Красноруцкий
Новосибирский государственный технический университет, Новосибирск, Россия
АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ НЕЛИНЕЙНОЙ КРАЕВОЙ ЗАДАЧИ ДИНАМИЧЕСКОГО ДЕФОРМИРОВАНИЯ ТОНКОГО СТЕРЖНЯ
В статье представлен алгоритм разработанной подпрограммы решения двухточечной краевой задачи для системы нелинейных дифференциальных уравнений первого порядка. Новый алгоритм ранее опубликованной подпрограммы KLPALG объединил в себе основные идеи подпрограммы BVPFD (DD14AD, PASVA3) и PASSIN, реализующей методику продолжения решения по параметру. Кроме того, приведены обобщенные результаты трудов авторов в задаче о нелинейном динамическом деформировании тонкого пространственного криволинейного стержня при расчете по его дифференциальной модели.
Неизвестные функции, входящие в уравнения движения, разыскиваются в дискретных точках. По методам прямого интегрирования производные по времени выражаются через текущие координаты и найденные на предыдущих шагах по времени. Первая производная по координате заменяется конечной разностью, добавляются краевые условия. Полученная система нелинейных алгебраических уравнений решается с помощью метода Ньютона с контролем длины шага из условия сходимости. Матрица Якоби этой системы имеет блочную трехдиагональную структуру, которая поддается эффективному LU-разложению. Такая декомпозиция матрицы Якоби позволяет быстро решать соответствующие системы линейных алгебраических уравнений больших размеров. Если условие сходимости метода Ньютона дает слишком маленький шаг, тогда применяется техника продолжения решения по параметру (псевдодлина дуги). После того как решена основная система нелинейных уравнений, для уточнения узловых значений вычисляемых функций применяется так называемый метод отсроченной коррекции (deferred correction method). Этот метод позволяет вычесть из получаемого решения ошибку, внесенную аппроксимацией производной по методу конечных разностей на начальном этапе численного решения. Получаемое таким образом численное решение имеет назначенную точность. Такая методика реализована в виде подпрограммы KLPALG, алгоритм которой представлен в данной статье.
Ключевые слова: нелинейная краевая задача, тонкий криволинейный стержень, дифференциальная модель, численное решение, KLPALG, BVPFD, DD14AD, PASVA3, PASSIN.
N.V. Pustovoy, V.E. Levin, D.A. Krasnorutskiy
Novosibirsk State Technical University, Novosibirsk, Russian Federation
THE NUMERICAL ALGORITHM FOR SOLVING
NONLINEAR BOUNDARY PROBLEM OF THIN ROD'S DYNAMIC DEFORMATIONS
At this paper the algorithm of the subprogram for solving a two-point boundary value problem for system of the nonlinear differential equations of the first order is presented. The new algorithm of the subprogram named KLPALG united in itself the main ideas of the subprograms BVPFD (DD14AD, PASVA3) and PASSIN realizing the technique of continuation of solution by parameter. Besides, the generalized results of works of authors in a problem of nonlinear dynamic deformation of a thin spatial curvilinear rod calculated by its differential model are presented.
The unknown functions in the equations of motion are calculated at discrete mesh points. The methods of direct integration allow us to express time derivatives by the current coordinates and coordinates and velocities calculated in the previous time steps. The first derivative of coordinate is replaced by finite difference; boundary conditions are added. The obtained system of nonlinear algebraic equations is solved by Newton method with the step length control of the convergence conditions. The Jacobi matrix of this system is of the block-tridiagonal structure which lends itself to efficient LUdecomposition. This decoupling of the Jacobi matrix allows you to quickly solve the corresponding system of linear algebraic equations of the big sizes. If the condition of convergence of Newton's method gives too small step, then used the technique of continuation of the solution on the parameter a (pseudo arc-length). As soon as the system of nonlinear equations is solved, to refine the nodal values of the calculated functions we use the deferred correction method. This method subtracts from the received solution the mistakes made by the approximation derived by the method of finite differences in the initial phase of the numerical solution. Thus obtained numerical solution is of accuracy appointed by user. This method is implemented in KLPALG subroutine which algorithm is presented in this paper.
Keywords: nonlinear boundary value problem, thin curvilinear rod, differential model, numerical solution, KLPALG, BVPFD, DD14AD, PASVA3, PASSIN.
Моделирование динамических процессов в ряде прикладных задач нелинейной механики деформируемого твердого тела удается построить на базе математической модели тонкого стержня [1-3] - упругого тела, один линейный размер которого существенно больше двух других. Введением ряда гипотез [4-6] деформирование стержня, по существу, заменяется деформированием его осевой линии, проходящей через центры жесткости поперечных сечений, ортогональных к этой осевой линии. Положение точки на стержне определяется естественной координатой - длиной дуги осевой линии. В геометрически нелинейной постановке задачи такая модель стержня позволяет описать множество разнообразных пространственных форм статического и динамического равновесия. В задаче о деформировании стержня наиболее распространены два подхода: на основе метода конечных элементов
и на основе решения краевых задач для систем дифференциальных уравнений. В данной работе используется дифференциальная модель стержня [7-10]. Преимущество используемой модели стержня перед другими известными дифференциальными моделями [2, 3, 11, 12] состоит в том, что в уравнения не входит начальная кривизна осевой линии стержня, что позволяет напрямую описывать стержни с изломами осевой линии и со скачками кривизны. При расчете динамического деформирования учитывается инерция, в отличие, например, от работы [12], где динамика представляется в виде набора статических конфигураций, а форма стержня представляет собой спираль. Задача о динамическом нелинейном деформировании тонкого криволинейного стержня представляет собой краевую задачу для системы нелинейных дифференциальных уравнений в частных производных. Решение этой задачи с помощью методов прямого интегрирования [13-15] сводится к решению последовательности двухточечных нелинейных краевых задач для системы нелинейных дифференциальных уравнений первого порядка. Для решения этих задач существует ряд подпрограмм [16-18], базирующихся на разработках прошлого столетия. Решение разыскивается в дискретных точках, производные в уравнениях заменяются конечными разностями [19, 18], строится итерационный процесс по методу Ньютона-Рафсона [20], также возможны и другие подходы, например [21]. При решении задачи о динамическом деформировании тонкого стержня на определенном этапе развития дифференциальной модели [8, 10] у авторов данной статьи появилась необходимость внести определенные изменения в текст подпрограммы БУРБО [17, 18], которая хорошо подходила для решения поставленной задачи о деформировании стержня до этого момента. Кроме внесения изменений в формирование нелинейных уравнений при применении численной методики развитие дифференциальной модели стержня [8, 10] потребовало добавить в подпрограмму метод продолжения решения по параметру в виде [22, 23]. В работах [18, 22-26], а также связанных с ними описаны основные идеи, реализованные в подпрограммах БУРБО [17], РАБУАЗ [18], РАББШ [22], но не приводится подробных алгоритмов, по которым можно было бы составить программу. После анализа указанных работ авторами данной статьи была разработана своя подпрограмма КЬРЛЬО [10], которая объединила в себе основные идеи указанных выше работ и позволила решать более сложные краевые задачи [8], чем
БУРЕБ [17], РАБУАЗ [18] (саму подпрограмму РАББШ [22] найти не удалось). После этого алгоритм КЬРЛЬО претерпел существенные улучшения, которые приблизили его по быстроте исполнения кода к БУРББ. Этот новый алгоритм КЬРЛЬО и представлен в данной статье.
1. Математическая модель тонкого стержня
Точки недеформированной осевой линии стержня задаются радиусом-вектором
е®=** ОН, (1)
где ^ - произвольный параметр, связанный с естественной координатой - длиной стержня э; ^ 2 3 - орты глобальной (неподвижной) системы координат; х12 3 - координаты точек осевой линии, здесь и далее
по повторяющимся индексам ведется суммирование от 1 до 3.
Для описания ориентации поперечных сечений стержня до деформирования используется «матрица поворота» (матрица начальной геометрии) р(^). Направляющие векторы главных осей инерции поперечных сечений недеформированного стержня записываются следующим образом:
е; ® = (!;Н, (2)
е3 - направлен по касательной к осевой линии стержня, здесь и далее
неуказанные индексы пробегают значения от 1 до 3.
Направляющие векторы главных осей поперечного сечения после деформации имеют следующее выражение через компоненты матрицы начальной геометрии в и матрицы поворота при деформировании X :
е;* =Р;* I. (3)
Жесткостные характеристики стержня задаются следующими величинами: EJl2 - изгибные жесткости; 033 - жесткость на кручение; Е¥ - жесткость на растяжение-сжатие. Жесткости также могут зависеть от текущей деформированной конфигурации - текущей кривизны осевой линии или удлинения.
Инерционные характеристики стержня: р^- погонная плотность; рJ1 2 3 - произведение плотности и моментов инерции поперечных сечений. Если ось жесткости стержня не проходит через центры инерции поперечных сечений [27], тогда задаются координаты центра тяжести поперечного сечения относительно центра жесткости Д12 . Вектор положения центров масс относительно центров жесткости
Й = Д1ё; + Д 2ё2 =Д1р1к "кпёп +Д2р2к^кИ!И , К =(( + Д 2^2 к )) • (4)
На стержень могут действовать внешние распределенные нагрузки и моменты. В общем случае они могут зависеть не только от текущей точки на осевой линии, но и от времени I, пространственного положения, текущей кривизны и других параметров:
Я = я(,ид,...), ш = ш(,ид,...). (5)
В качестве внешней распределенной нагрузки может выступать, например, собственный вес, аэродинамическая нагрузка от набегающего потока воздуха [7], силы от контактного взаимодействия витков осевой линии [28], силы инерции от пробегающего внутреннего потока жидкости и др. Кроме распределенных силовых факторов на стержень могут действовать сосредоточенные силы и моменты. Сосредоточенные нагрузки на концах стержня задаются через краевые условия. Задание сосредоточенных сил и моментов в других точках приводит к многоточечной нелинейной краевой задаче, которая может быть решена в рамках численного алгоритма решения двухточечной нелинейной краевой задачи с некоторыми модификациями в алгоритме.
Мгновенные угловые ускорения и скорости поворота поперечного сечения стержня относительно главных осей е12 3 имеют следующий вид [8]:
^ 1 = —кт„в3ткк„Р2к — кт„в3тккпв2к , = ^3^2 =—кт„в3тккпв2к ,
^2 = кт„в3ткк„в1к + ктпв3ткк„в1к , ^2 = ^зХ = кт„в3ткк„в1к , (6)
^3 = кт„в1ткк„в2к + кт„в1тккпв2к , ^3 = = ктпв1тккпв2к ,
где точка обозначает производную по времени I.
Система уравнений движения тонкого упругого стержня, описывающая любые повороты поперечных сечений [8] при динамическом деформировании, имеет вид
и1Л =(1+е)х; {К ]г- ; К^ = (КкгА<»р,5 -Ккр^,5) ¿4= КшвыС^кКк№п
С =
кз-1 0 0
0 КЗ 2-1 0
0 0 ОЗ.
-1
3 у
(7а)
(7б)
а л=
Мг ,5 =
р^ (5).((г + Я,)-дг (5,...)] ^(1+е);
л
' рЗп (5)12 пвпк К ^ +
(5). г((р + я р)-(( г + я г )я
V
(1+4
(7в)
г = 1,2,3, Р = 2,3,1, (7г) г = 3,1,2;
(5,->,5 + Х],5К¡гО-р -Х;,5К;р°
е=х;,5Ка\_ЕР(5)-5,5] 1 +«АГ , ^,5 = ,5хк,5)1/2,
где а - коэффициент температурного удлинения; АТ - изменение температуры; нижний индекс 5 после запятой означает производную по этой переменной.
Система (7) состоит из 18 нелинейных дифференциальных уравнений первого порядка по 5 и второго порядка по времени Разрешающими функциями системы (7) являются 18 функций: и123 (5,*) -проекции вектора перемещений точек осевой линии стержня на оси глобальной системы координат; К у (5,*) - компоненты матрицы, описывающей поворот главных осей инерции поперечных сечений при деформировании; О123 (5,*) - глобальные проекции вектора внутренних усилий; М12 3 (5,*) - глобальные проекции внутренних моментов.
На девять компонентов матрицы поворота X (5,*) накладываются
дополнительные зависимости, которые вытекают из общих свойств матриц поворота. Поэтому на практике вместо 9 компонент можно использовать только 3 функции в качестве параметров поворота, что
уменьшает число неизвестных функций и число уравнений с 18 до 12. В качестве параметров поворота могут выступать, например, углы Эйлера, самолетные углы, компоненты вектора конечного поворота [29, 30] и др. При использовании вектора конечного поворота угол поворота тройки ортов, связанных с поперечным сечением стержня, ограничивается величиной 2п [11] в общем случае. Для использования вектора конечного поворота в качестве параметра для любых поворотов при деформировании разработана специальная численная методика [8], суть которой будет представлена ниже. Если деформирование стержня укладывается в рамки ограничения поворота при деформировании 2п, тогда группа уравнений (7б) заменяется тремя уравнениями следующего вида:
Э1ти d^J в » в = МРК1 кр ds
Э®] dc¡ Е3Х d\
дЬтеА.* 1 в - к1кр ds Э®] dc¡ Ш2 dc¡
> в = МрРзк1 кр ds эю; в1т 1 кпв2к 03ъ ^ ,
где й - №. ^ - вектор конечного поворота. Таким образом, имеем 12 дифференциальных уравнений и 12 неизвестных функций и123, ю12 3, 2 3, М12 3. Такое же количество уравнений и неизвестных получается при применении численной методики [8], но при этом ограничения вектора конечного поворота (при углах поворота кратных 2п определитель СЛАУ (8) равен нулю) не мешают получить численное решение краевой задачи.
Уравнения движения (7) являются системой нелинейных дифференциальных уравнений в частных производных. Для решения можно применять одношаговые и многошаговые методы прямого неявного интегрирования по времени [13, 14, 9] (методы Ньюмарка, Хаболта, Парка и др.). Основная идея численного интегрирования уравнений динамики стержня заключается в том, что производные по времени (скорости и ускорения) аппроксимируются специальными выражениями через координаты (перемещения и повороты). В эти выражения
также могут входить известные координаты, скорости и ускорения с предыдущих шагов по времени. Применение таких методов позволяет снизить размерность исходной задачи нелинейной динамики стержня, сводя ее, по существу, к задаче нелинейной статики.
Конфигурация стержня в каждый момент времени описывается следующими векторами-столбцами:
Х(^)={и12,з ),Ш1Я3 ),а,2,з ),М12,3 )}Т, XX), X).
После проведения ряда тестов, авторы статьи пришли к выводу, что наиболее подходящим (по точности и устойчивости) для решения обозначенной задачи из рассмотренных методов прямого интегрирования является метод Парка [15, 13, 9]. Для расчетов используется модификация [9] этого метода:
X (^ ) = аХ^_3) + а2 Х(Г_2) + а3 Х^^) + а4 Х^0) + а5 Х(^), (9)
где
1 дЦ 1 dL4 1 dL5 1 1 dL. 1 dL5
a,---, a2 =--4+---, a3 =--1+--4+--1
14 д/1 2 4 д/ 4 д/2' 3 2 д/ 4 д/2 4 д/з
1 дЦ, 1 dL. 1 дЦ5 1 dL3 1 дЦ4 1 дЦ5
a. =--1+--4+--a5 =--3+--4+---.
4 2 д/2 4 д/3 4 д/4' 5 2 д/з 4 д/4 4 д/
л
Здесь Ln =Цft), Ln(х) = £ / П
х - X,
j=1 ^ i=0, j*i Xj Xi
- полином Лагранжа для
интерполяции по n точкам.
После замены скоростей и ускорений выражениями (9) в системе (7) остаются только первые производные по параметру ^ . Для замыкания задачи необходимо добавить краевые условия - 6 на левом конце и 6 на правом конце стержня.
2. Численная методика решения нелинейной краевой задачи
Представленная здесь численная методика решения нелинейной краевой задачи базируется на методике [18], которая в свое время породила несколько версий подпрограмм. Из числа известных и доступных подпрограмм это BVPFD [17] из библиотеки IMSL (International Mathematical and Statistical Library) и DD14AD из HSL Mathematical
Software Library. Помимо пользовательского интерфейса эти подпрограммы различаются алгоритмом работы в части критерия расходимости итераций метода Ньютона. Существенная разница в работе этих подпрограмм может быть устранена изменением текста DD14AD. Для поставленной задачи был разработан другой алгоритм, который соединил в себе основные идеи [18] и алгоритм продолжения решения по параметру [22, 23] (подпрограмма PASSIN), что оказалось незаменимым при реализации численной методики, позволяющей использовать вектор конечного поворота при любых деформациях осевой линии [8, 10]. Первая версия такого алгоритма представлена в работе [10], однако после сравнительного анализа процесса решения тестовых краевых задач подпрограммами KLPALG [10] и BVPFD [18, 17] наш алгоритм претерпел существенные улучшения. Новая версия подпрограммы KLPALG, алгоритм которой представлен в этой статье, выполняет одинаковое с BVPFD количество вызовов функций расчета невязок и факторизации матрицы Якоби на тестовых примерах, где возможно использование BVPFD. На более сложных примерах, где BVPFD выдает расхождение метода Ньютона, в подпрограмме KLPALG включается алгоритм продолжения решения по параметру [22, 23], который приводит к искомому решению.
Согласно [18] неизвестные разыскиваются в дискретных узлах, а дифференциальные уравнения аппроксимируются по правилу трапеций:
Wj +1 - Wj = 1 h} 2
f ((j, W;)+f (;+1,w;+1 )1, j=1,2..N, (10)
где в нашем случае системы уравнении стержня:
0 = Ь}1 <£,2 <...<£,ы+1 - узлы, в которых разыскиваются значения
неизвестных функций; Ъ^ =Ъ,+1 — ^- шаг сетки;
W ={и1,и2,из,Ы1,Ы2,Щ,0>1&2£з,М1,М2,Мз}Т - вектор-столбец 12 узловых неизвестных;
f (W) - правая часть системы дифференциальных уравнений движения тонкого упругого стержня (7)-(8) с учетом (9).
Краевые условия представляются в виде g1 (\¥1) = 0, g 2 (\¥1, WN+1 ) = 0 , где g12- векторы-столбцы, g1 имеет р строк, g2 имеет (12 - р) строк, q из которых связывают W1 и WN+1. Если начало стержня не связано с концом, тогда р = 6, q = 0, следовательно, g1 и g 2 имеют по 6 строк (по 6 краевых условий на каждом конце стержня), таким образом, в дискретной форме краевая задача может быть представлена в следующем виде:
^) = 0,
(11)
где
W :
W2
W
N+1.
\Ж 1 1,1 ^(1)
Ж "2,1 и 21)
Ж ''12,1 =■ = ■= м 30)
Ж ''1,2 и1(2)
Ж /'12^+1 _ м 3N+1)
FЯ(W)
gl (^)
- Wl - \ ((1 + * 2 )/2
Wз - - ¿2 (2 + *3 )/2
WN+1 WN ^ ( + fN+1 )/2|
82 (Wl,WN+1 )
где *; =*(;,W;), = $; +1 Ч;, ] = ^ .
Решение разыскивается итерационным способом:
W^ = Wу + ^ • AWу = Wу + ^ • ^(Wу)-1 Fп(W*), (12)
где FЖ (W- матрица Якоби системы уравнений (11) в точке Wу для
V -го приближения; ^ - параметр контроля шага, определяется особым образом [18]. Матрица Якоби имеет следующую структуру:
Fw ( >
ёщ (^ )
-I)/2 0
ё 2Щ +1 )
I-((2Щ^1 )/2 -I -(2щ^2 )/2
0 0 0
0
1 ((+\)Щ^Ы )/2
ё2ЩЖ+! С^1,^^+1 )
Для определения неизвестных узловых значений функций используется система нелинейных уравнений (11). Численная методика [8] предполагает заменить в (11) уравнения для узловых значений проекций вектора конечного поворота
юк, 1 +1 -юк, 1
-ь] (( + У3+к,1+1 )/2 = 0, к = 1,2,3, 1 = 1,2...#, (13) на уравнения
"к1п,1 +1 -к1п,1 - ^ (( + /П,! +1 )) (к 2 п,1 + к 2 п, 1 +1 )/2 = 0
x, 1+1 - кт,1 - ^ ((+/; 1+1)/ 2] (к 3п,1 +к 3п,1+1))=0, к 2 п,1+1 -к 2 п,1 - ^ (( + ))(3п,1 +к 3п,1+1 )/2 = 0
(14)
где по индексу п ведется суммирование от 1 до 3, а по индексу 1 нет суммирования. Этот индекс указывает на номер блока уравнений 1,2...N.
В итерационном процессе (12) для поиска приращения ЛWУ на каждом шаге решается система линейных алгебраических уравнений:
¥щ ( ^ *). (15)
Матрица Якоби Fщ (Wv) имеет следующую блочную трехдиаго-нальную структуру:
0
0
0
?№ (>А ^
г А1 С1 0 . .. 0 0
В2 А2 С 2 . .. 0 0
0 Вз Аз . .. 0 0
0 0 0 . . Аы С N
v Б 0 0 . . ВЫ+1 А ^ы+1у
(16)
где Ак, Вк, Ск - матрицы 12*12, нулями обозначены матрицы с нулевыми элементами, матрица Б = g ш (W1,WJV+1) имеет q ненулевых элементов, если начало стержня связано с концом. Матрица Вк имеет р ненулевых строк, Ск имеет (12 -р) ненулевых строк:
Вк
...
v 0 у
Ск
Г 0 ^
V--/
Если А V Ф 0, тогда матрица (16) может быть представлена в виде разложения
г I 0 0 ... 0 0 > г а 1 С1 0 . .. 0 0
в 2 I 0 ... 0 0 0 а2 С2 . .. 0 0
и = 0 вз I ... 0 0 0 0 аз . .. 0 0
0 0 0 ... I 0 0 0 0 . . аы С N
v 51 5 2 5 з .. (в N+1 + 5N ) Iу v 0 0 0 . .. 0 а N+1у
где матрицы ак и вк удовлетворяют следующим уравнениям:
«1 = А!, а. = А. - в. С61а1 = Б, . = 2,3. ..Ы+1;
в ..а .-1 = В., 5. а. =-б;-1С;-1, . = 2,3...N+1. (17)
Решение системы (15) сводится к решению двух систем
Ь• У = -Б;(*), и• AW* = У . (18)
Алгоритм подпрограммы ЬЦ-факторизации матрицы Якоби следующий.
На входе в подпрограмму имеются массивы аг = Аг, в= Вг, Сг, 51 = Б, они являются рабочими массивами и в них остается результат.
1. Начало цикла, г = 1.
2. Ш-факторизация матрицы аг:
2.1. переставляются столбцы на первых Р строках так, чтобы на главной диагонали стояли максимальные по модулю элементы;
2.2. переставляются строки оставшейся части матрицы аг и матрицы в г так, чтобы на главной диагонали стояли максимальные по модулю элементы (при этом строки могут меняться в матрицах Сг и аг+1 );
2.3. рассчитывается Ш-разложение матрицы аг, на главной диагонали которой стоят максимальные по модулю элементы.
3. Вычисляется 6г = 6а6г+1 =-6гСг , вг+1 = вi+1а-1, аж =
= аг +1 - вг +1Сг .
4. г = г + 1, если г < N+1 идем п. 2.
5. Вычисляется а N+1 = а N+1 + 6 N+1.
6. Производится Ш-разложение аN+1, как в п. 2.3, при этом
предварительно переставляются столбцы так, чтобы на главной диагонали стояли максимальные по модулю элементы.
На выходе из подпрограммы LU-факторизации матрицы Якоби имеем измененные массивы аг, вг, Сг, 6г и массивы перестановок
столбцов и строк.
Для внесения в расчетную схему дополнительных сосредоточенных сил и моментов, приложенных не на концах стержня, необходимо добавить псевдоузел в точке приложения силы. То есть необходимо в вектор (11) внести дополнительные матричное уравнение, ко-
торое будет стыковать старый узел, в котором приложен сосредоточенный силовой фактор и новый псевдоузел, по расположению совпадающий со старым. Дополнительное матричное уравнение должно содержать равенство всех функций слева и справа, за исключением внутренних усилий или моментов, которые должны иметь скачок, равный приложенному внешнему сосредоточенному силовому фактору.
Увеличение числа дополнительно приложенных сил и моментов не изменит вид (16), просто возрастет количество неизвестных, а следовательно, в алгоритме решения полученной таким образом нелинейной системы ничего не изменится. Другими словами, процесс добавления сосредоточенных силовых факторов отразится на формировании нелинейных уравнений (11). По сути, такой способ эквивалентен разбиению стержня в местах приложения сил и моментов на несколько составных частей.
Стратегия выбора шага модифицированного метода Ньютона (12) следующая. Длина шага ^ выбирается максимальной из ряда
{1,1/ 2,1/'4,...}, для которой выполняется условие сходимости
подпрограмме DD14AD. В новой версии подпрограммы KLPALG выбрана величина c = 0,11 при сравнении процессов расчета с BVPFD.
После того как решена система (11), полученное решение нужно скорректировать (уточнить) с учетом той ошибки, которую внесла аппроксимация производной по методу конечных разностей (10). Для этого используется метод отсроченной коррекции (deferred correction method) [18, 25, 26]. Рассмотрим кратко применение этого метода к нашей задаче (более подробные выкладки представлены в [10]).
Пусть W* - точное решение, тогда так называемая локальная ошибка аппроксимации (local truncation error) [18], обусловленная использованием приближенной формулы (10) c учетом W*J = f(j,Wj), имеет следующее выражение:
Необходимо отметить, что локальная ошибка аппроксимации (19) записана для точек, находящихся посередине между узловыми точками ^у, у = 1,2.+1, шаг сетки в общем случае Ъ1 = £,/+1 /, у = 1,2..
r(Wv)-r(Wv + ^vAWv)>c• r(Wv),
где r(x) = ||X|2, ||X||n = П)), c = 1 в [18] и c = 0,5 в [24, с. III.3-2]
и в
) , j = 1,2..N, где ak = ----Г1 + (-1)-] . (19)
Получим такую аппроксимацию (19), какая предлагается в [25]. Пусть задана гладкая функция у(^) и произвольная сетка с узлами ^., - -1,2...#+1. Определим следующую функцию:
т у (к)(£)
¿0^)-, (20)
к-0 к!
где М^) такая функция, что ;+1/2 ) = л(((;+1 - )/2);+1 Ч-, --1,2...#.
Справедливо следующее равенство:
т(^)- Ь( *£)+О(т+1).
Аппроксимируем (20) линейной комбинацией значений функции у(^) в узловых точках:
т+1
М(у)-), (21)
.-1
где г - заданное целое значение (от которого зависит, какие узлы участвуют в аппроксимации); wi - пока неизвестные весовые коэффициенты.
Пусть необходимо достигнуть точности порядка О(кт+1). Разложим (21) в ряд Тейлора:
т+1 - у (-)(£)
м (у)-1 ^ 2 ^5%+г чг -
.-1 - (22) - 2 Ж-^т1+о(/.^)т+1),
1-0 V.-1 у
где а. -((.+гЧ)/М().
Сравнивая (22) с (20), получаем СЛАУ относительно неизвестных весовых коэффициентов wi:
т+1
2^. а- - а-, - - 0,1...т . (23)
.-1
Находя wi из (23), получим М(у) = £(у,^) + о(и(^)т+1^
Система уравнений (23) является Вандермондовой системой уравнений [25]. Быстрый и эффективный алгоритм решения этой системы описан в [31], там же приведен текст подпрограммы на языке «Алгол». Таким образом, находя весовые коэффициенты из (23), можно по узловым значениям W строить аппроксимацию локальной ошибки аппроксимации (19).
Метод отсроченной коррекции [25] предполагает приближенный расчет локальной ошибки аппроксимации (19) приближенно, по сути, с помощью метода конечных разностей (производные функции выражаются через узловые значения этой функции). Рассмотрим алгоритм расчета 8 к (W) - конечно-разностной аппроксимации (19).
Разложение локальной ошибки аппроксимации (19) содержит члены, начиная со второй производной функции {(W*+1/2), и имеет
только четные производные, следовательно, аппроксимация 81 (W) будет составляться по 4 узловым точкам, 8 2 (W) - по 6, 8 3 (W) - по 8, 8к - по (2к + 2) узловым значениям. Поэтому корректное построение 8к (W) требует определенного числа разбиений по длине
Для нахождения 8 к (W) в аппроксимации (21) следует положить т = 2к +1 и
стержня.
г (] ) = ш1и [# _ 2к _ 1, тах(0, ] _ к _ 1)], ] = 1,2...#,
тогда коэффициенты аг. в (23) имеют вид
аг (У) = + ^ )/2]). (24)
(24)
Конечно-разностное представление 8к (W) локальной ошибки аппроксимации (19) для (14) имеет следующий вид:
8 к («3 ) = [к 3т, ] +к 3т,1 +1 > к 2т )/2.
Пусть W0, рассчитанное с точностью 0(к2) решение уравнения, а (W0) - конечно-разностное представление локальной ошибки аппроксимации с точностью 0(к2). Тогда после решения СЛАУ
Fщ )А = -Ь1 (W0) (26)
получаем аппроксимацию О (к2) глобальной ошибки Wз - W0 (где
Wз - точное решение), т.е.
А = W3 - W0 + 0(к2).
Следует отметить, что для нахождения решения (26) не требуется много вычислений, так как после решения FJt(W) = 0 уже имеется ЬЦ-
декомпозиция (17) матрицы Якоби Fщ (W0).
Далее, решая нелинейную систему уравнений
Fп(W) = 81 (W0), (27)
получаем аппроксимацию решения с ошибкой 0(к4), т.е. решение нелинейной проблемы W1 удовлетворяет следующему равенству:
W1 - W0 = 0(к4).
Это первый шаг метода отсроченной коррекции. Далее аналогично (26) для оценки текущей глобальной ошибки решаем СЛАУ:
Fщ (W1 )А = (W0)-Ь2 (W1) (28)
и получаем аппроксимацию 0(к4) глобальной ошибки:
А = W3 - W2 + 0(к4).
Другой не менее важный аспект численной методики решения нелинейной краевой задачи - это алгоритм продолжения решения по параметру. В некоторых задачах, в частности при использовании в ка-
честве параметров поворота при деформировании на большие углы компоненты вектора конечного поворота [8, 10], сходимость метода Ньютона может быть очень медленной. Зачастую время расчета такого итерационного процесса при выборе длины ньютоновского шага (12) из условия сходимости метода стремится к бесконечности. Встает вопрос о выборе хорошего начального приближения. Подробно методика изложена в [23, 22, 10]. Здесь приведем основной алгоритм.
Процедура продолжения решения по параметру [10] имеет следующий вид.
1) Пусть дано решение ,Yk]. Вычисляем--, решая СЛАУ
dу
. dW,
Fw (Wt ,Y,)—k = -Fy(Wk, Yk ) = -F;(W0).
(29)
dY TV
2) Определяем знак Якобиана sfk = sign [FW (Wk, Yk ) (он меняется
при переходе через точки бифуркации), вычисляем Wk и уk по следующим формулам:
Y k
sf0 • Sfk
f
i+e
dW* dW,
Wk =Y k
-1
dWk d y
(30)
ё у ёу
где sf0 - выбирается сначала при Y = Y0 = 0 так, чтобы у0 > 0, так как нам необходимо изменять параметр Y от 0 до 1; I * (о0) такое, что
I * (о 0 )1 (о 0 ) = ||1 (о0 )||\ 0е (0,1) - коэффициент.
3) Рассчитываем начальное приближение (делаем прогноз):
( W 0 >
yyk+1 vYk+i j
(Wk ^ ( W, ^ (Wk ^ (AW;
vYk j
+
vYk j
Ao =
vYk j
+
0
ay;
v" Ik j
4) Для получения уточнения
( w V+1 ^
TTk+1 YV+1
. Ik+1
(коррекции) имеем следую-
щие выражения:
Fw(Wk,Yk)• *V = -F(W;+1 ,YV+1), FW(Wk,Yk)• yV = -FT(W^,yV+1 ),
Ayv+1
AII/V+l V . A„,V+1 V
, AW = z +Ay y ,
WV+1 n jV . hijV+1 V+1 V . a v+1
= W* + AW , Y V+1 = Y k + aY .
(31)
Здесь
NN3 (W(C1 ),у(^1 ),01 ) = 0- W * (o 0 )[w(o)-W(o 0)] +
+(1- 0) • у(°o )[y(o)- Y(oo)] - (o- o0), О = o0 + Ao
или
NV =0^ w; [w;+1 -Wt ] + (1 -0) Yk [y k+1 - Yk ]-Ao
Описанная процедура встраивается в основную процедуру решения нелинейной проблемы особым образом, алгоритм подпрограммы представлен ниже.
3. Алгоритм подпрограммы
На входе в подпрограмму имеется: начальное приближение решения W0; N - число разбиений по длине стержня; е - требуемая точность решения задачи; min д - минимальный коэффициент ньютоновского шага (0 < min^< 1), ниже которого включается алгоритм продолжения решения по параметру; 0е (0,1) - коэффициент для нормализации псевдодлины дуги (pseudoarclength normalization) [23, 22].
k = 0 - счетчик алгоритма метода отсроченной коррекции (deferred corrections method);
PM = 0 - алгоритм продолжения решения по параметру выключен;
Y 0 = 0, y1 = 1 - переменные для алгоритма продолжения решения по параметру;
PE = ложь - метка стадии завершения работы алгоритма продолжения решения по параметру;
Инициализация рабочих переменных
S1 = S1 (W0 ) = 0 - аппроксимация локальной ошибки усечения (local truncation error), вектор-столбец размерностью 12(N+1);
S0 - вспомогательный рабочий вектор-столбец 12(N+1).
brl = га (некое большое число, например 10 ) - для контроля расхождения метода Ньютона;
W = W0 - рабочий вектор-столбец размерностью 12(N+1);
W 2 - вспомогательный рабочий вектор-столбец 12(N+1). ITNEW = 0 - счетчик итераций метода Ньютона; NEPS = k • h2, где k - маленькая константа, h - шаг сетки по переменной ^;
bTOL = - для контроля расхождения метода отсроченной коррекции;
CASI = ложь - метка о выключении перерасчета Якобиана (квазиньютоновский метод).
Цикл итераций метода Ньютона
1. Если PM > 0, тогда D = -Fn(W)-(у1 -1)Fn(W0), иначе
D = -Fn(W) + S1; r = ||D||2, r1 = 4T;
2. Если r1 <е, тогда {
2.1. Если PM = 0, тогда идем в 25;
2.2. Если PM > 1, тогда { PM = 3 ; CASI = ложь } ()
3. Если CASI = правда, тогда идем в 5;
4. Расчет Якобиана FW (W), градиента G(W), факторизация
f- (W);
5. Расчет AW=F- (W )• D;
6. Если ERREV, тогда идем в 18;
7. Если PM > 0, тогда идем в пункт 33;
8. Если ITNEW > 0, тогда br1 = r1;
9. br 2 = r1;
10. rG = -(GT •AW);
11. Если CASI = ложь и (rG < 0 или |r -rG| > r ), тогда AW = -G;
12. W 2 = W+AW, ц = 1;
13. Расчет Fn(W2), D = -Fn(W2) + S1, r = ||D||2, r1 = jr ;
14. Если (br1 - r1)< 0,11ц • br 2, тогда { (контроль шага)
14.1. | = |0.5;
14.2. W2 = W+ |AW ;
14.3. Если |< min |, тогда {
Если k > 0, тогда {
Метод отсроченной коррекции разошелся, нужно уточнить сетку, возврат с ошибкой IERR=4} Иначе {
Запускаем алгоритм продолжения решения по параметру PM = 1, Y = 0; Идем в 1 . } }
14.4. Идем в 13. }
15. W = W 2;
16. Если ITNEW < ITMAX, тогда { ITNEW = ITNEW +1, идем в пункт 2. }
17. Достигнут предел по количеству итераций, метод Ньютона разошелся или имеет слишком маленькую скорость сходимости, выход с ошибкой IERR=3.
Оценка точности решения краевой задачи
18. ERREV = ложь ;
19. Расчет cTOL = ||AW||2, cTOL = VcTOL ;
20. Если cTOL < TOL, тогда { проблема решена успешно IERR=0, выход. }
21. Если cTOL > bTOL , тогда {
Метод отсроченной коррекции разошелся, нужно уточнить сетку, возврат с ошибкой IERR=4.}
22. bTOL = cTOL ;
23. ITNEW=ITNEW+1;
24. Идем в пункт 1.
Метод отсроченной коррекции (Deferred Correction Method)
25. Сохраняем S0 = S1;
26. S1 = Sk+1 (W) по формулам (25);
27. D = S0 - S1;
28. ERREV = истина; (отправка на оценку достигнутой точности);
29. к = к +1;
30. ITNEW = 0, br1 = 10300 , bTOL = 10300 ;
31. CASI = истина ;
32. Идем в пункт 3.
Метод продолжения решения по параметру
33. Решаем СЛАУ (29) FW (W)^ = -Fn(W0);
34. Знак Якобиана sf = sign(FW);
Корректор
35. Если PM = 2, тогда {
35.1. Вычисляем нормализацию Nv по формуле (31)
N V=e-W; [W-W 2] + (1 -e)-Y 0 [y,-Y 0 ]-До; Вычисляем ДY1 по формулам (31):
dW'
аъ =- |> v+e- w0*a] /
(i-e) о-e- W0
0 d у
a = a+äyj dW; W = W+A, Yi =Yi +A^; d y
=IAL• где Kl=njSlXr;
35.2. Если r2 > r1 или Ay1 >Ay0 или (y1 <y0 и РЕ = ложь ), тогда {
Отключаем стадию завершения алгоритма РЕ = ложь ; Уменьшаем шаг Ao = 0,5Ao;
AYi =Y0 • Ao, Yi =Y0 +AYi; Если Ayi < minAy , тогда
{слишком маленький шаг по параметру продолжения , IERR=3, выход.} W = W 2+W0 Ao;
Идем в 1.} =r2;
35.4. Идем в пункт 1.}
35.3. ri = r2;
Предиктор
36. Если РМ - 3, тогда {
36.1. Вычисляем Ж0 и у0 по формулам (30)
У 0 - ^0 •
(
1+0
dW dW
1
-1/2
, ^ -у,
d у
АУ 0-У 0 Ао;
d у d у
36.2. Если Ау0 < 0, тогда $/0 - —л/0, идем в пункт 36.1;
36.3. у 2 -У1+АУ0;
36.3. Если ( РЕ - ложь и у2 > 1) или РЕ - истина, тогда { Ставим алгоритм на завершение РЕ - истина Целим в единицу: Ао - (1 — у1) / у0;
АУ0-У0Ао, у 2-У1+АУ0;
Если |Ау0 < 8, тогда у1 -1, РМ - 0, идем в пункт 25. }
36.4. у0 71 -72;
36.5. Запоминаем текущий вектор W2 - W;
36.6. Делаем прогноз: W - W+^"0Ао;
36.7. Включаем алгоритм на коррекцию РМ - 2 ;
36.8. г1 -— (некое большое число, например 10300);
36.9. Идем в пункт 1; }
Предиктор, первый шаг
37. Если РМ -1 , тогда {
37.1. Вычисляем Цу0 и у0 по формулам (30)
У 0 -
(
1+0
dW_ d у dу
V
1
1/2
, ^ -у,
d у
37.2. Выбираем самый первый шаг Ао-1.0/у0 (условно из расчета, что потребуется один шаг для достижения у-1);
37.3. АУ0 -У0 -А°;
37.4. Если Ау0 < 0, тогда sf0 -—1, Ау0 -— Ау0; иначе -1;
37.5. Запоминаем текущий вектор W2 - W ;
37.5. Делаем прогноз: W - W+^Т0Ао, у1 -у0 +Ау0;
37.6. Ставим на коррекцию алгоритм продолжения решения
РМ = 2;
37.7. г1 (некое большое число, например 10300);
37.8. Идем в пункт 1. }
Конец алгоритма.
Основное отличие ББ14ЛБ от БУРЕБ - это наличие в последней условия как в пункте 8 приведенного выше алгоритма. Это условие, по сути, отменяет проверку сходимости после первого шага метода Ньютона.
Выводы
Тестирование данного алгоритма проводилось на достаточно большом количестве задач. Например, были повторены расчеты, опубликованные в работах [28, 32]. Точных тестов на скорость не проводилось, но предварительно можно заключить, что разработанная новая подпрограмма КЬРЛЬО, алгоритм которой представлен в данной статье, не уступает ББ14ЛБ и БУРЕБ по скорости исполнения кода и способна решать более широкий круг нелинейных краевых задач. Новизна представленных в данной статье результатов заключается: во-первых, в обобщении результатов многолетней работы авторов данной статьи в области механики тонких криволинейных стержней (представлена краткая сводка основных формул с описанием); во-вторых, представлен новый алгоритм подпрограммы КЬРЛЬО и описана используемая численная методика, что может помочь читателю в сжатые сроки ознакомиться, понять и реализовать в виде подпрограммы большой массив информации, представленный в прошлом столетии в ряде публикаций на английском языке, в том числе с ошибками, опечатками и недомолвками.
В представленной методике не рассмотрен вопрос о выборе оптимальной сетки, этот алгоритм [18] планируется внедрить в КЬРЛЬО. Другой нерассмотренный аспект подпрограмм БУРЕБ и ББ14ЛБ -это процедура продолжения решения, когда решение предыдущей краевой задачи экстраполируется вперед по методу Эйлера и при этом происходит программный выбор длины такого шага [33, 34]. Такой подход оказался неработоспособен в задачах [8, 10], поэтому в подпрограмму КЬРЛЬО был внедрен другой алгоритм продолжения решения по параметру [23, 22, 10].
Работа выполнена при финансовой поддержке Министерства образования и науки РФ по государственному заданию №2014/138, проект № 435.
Библиографический список
1. Светлицкий В.А. Механика стержней: учебник для втузов: в 2 ч. Ч. 1. Статика. - М.: Высш. шк., 1987. - 320 с.
2. Жилин П.А. Прикладная механика. Теория тонких упругих стержней: учеб. пособие. - СПб.: Изд-во Политехн. ун-та, 2007. - 101 с.
3. Светлицкий В.А. Строительная механика машин. Механика стержней: в 2 т. Т. 2. Динамика. - М.: Физматлит, 2009. - 383 с.
4. Kirchhoff G.R. Ueber das Gleichgewicht und die Bewegung einer elastischen Stäben // Crelle Journal fuer die reine und angewandte Mathematik. - 1858. - Bd. 56. - S. 285-313.
5. Кирхгоф Г. Механика. - М.: Изд-во АН СССР, 1962. - 402 с.
6. Ляв А. Математическая теория упругости.- М.; Л.: Объединение научно-технических издательств, 1935. - 674 с.
7. Пустовой Н.В., Левин В.Е., Красноруцкий Д.А. Применение геометрически нелинейных уравнений стержня к расчету статики и динамики тросов. Ч. 1 // Научный вестник Новосиб. гос. техн. ун-та, 2012. - № 1 (46). - С. 83-92.
8. Пустовой Н.В., Левин В.Е., Красноруцкий Д.А. Методика вычисления параметров больших поворотов поперечных сечений гибкого стержня при расчетах в рамках его дифференциальной модели. Ч. 1 // Научный вестник Новосиб. гос. техн. ун-та. - 2013. - № 2 (51). -С.155-164.
9. Красноруцкий Д.А., Левин В.Е., Пустовой Н.В. Нелинейная динамика тонких упругих стержней // Нелинейные колебания механических систем: труды IX Всерос. науч. конф. (Нижний Новгород, 2429 сентября 2012 г.) / под ред. Д.В. Баландина, В.И. Ерофеева, И.С. Павлова. - Н. Новгород: Наш дом, 2012. - С. 557-565.
10. Пустовой Н.В., Левин В.Е., Красноруцкий Д.А. Методика вычисления параметров больших поворотов поперечных сечений гибкого стержня при расчетах в рамках его дифференциальной модели. Ч. 2 // Научный вестник Новосиб. гос. техн. ун-та. - 2013. - № 3(52). -С.146-159.
11. Сорокин Ф.Д. Прямое тензорное представление уравнений больших перемещений гибкого стержня с использованием вектора конечного поворота // Известия РАН. МТТ. - 1994. - № 1. - С. 164-168.
12. Nonlinear dynamic deformation simulation for helical rod like objects / H. Du, W. Xiong, H. Wang, Z. Wang, B. Yuan // Engineering Review. - 2013. - Vol. 33. - Iss. 3. - P. 233-238.
13. Bathe K.J. Finite Element Procedures. Englewood Cliffs. - NY: Prentice Hall, 1996. - 1037 p.
14. Newmark N.M. A Method of Computation for Structural Dynamics // Journal of Engineering Mechanics Division, ASCE. - 1959. - Vol. 85. -No. EM3. - P. 67-94.
15. Park К.С. An improved stiffly stable method for direct integration of nonlinear structural dynamic equations // Journal of Applied Mechanics, ASME. - 1975. - Vol.42. - Iss. 2. - P. 464-470.
16. Shampine L.F., Muir P.H., Xu H. A User-Friendly Fortran BVP Solver // Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM). - 2006. - Vol. 1. - No. 2. - P. 201-217.
17. IMSL: Fortran Numerical Library. User's Guide. Math Library. Version. 7.0, available at: http://www.roguewave.com/documents.aspx7en-tryid= 519&comma nd=core_download. - Date 27.05.2014. - Title from screen.
18. Pereyra V. Pasva3: An adaptive finite difference fortran program for first order nonlinear, ordinary boundary problems // Lecture Notes in Computer Science. - 1979. - Vol. 76. - P. 67-88.
19. Rashidinia J. Finite difference methods for a class of two-point boundary value problems // IUST International Journal of Engineering Science. - 2008. - Vol. 19. - No. 5-2. - P. 67-72.
20. Вайнберг A.M. Математическое моделирование процессов переноса. Решение нелинейных краевых задач. - М.; Иерусалим, 2009. - 209 с.
21. Dinkar Sharma, Ram Jiwari, Sheo Kumar. Numerical Solution of Two Point Boundary Value Problems Using Galerkin-Finite Element Method // International Journal of Nonlinear Science. - 2012. - Vol. 13. -No. 2. - P. 204-210.
22. Lentini M. Boundary Value Problems over Semi-Infinite Intervals: Ph.D. Thesis / Cal. Inst, of Technology. - 1978. - 123 p.
23. Keller H.B. Constructive Methods for Bifurcation and Nonlinear Eigenvalue Problems // Lecture Notes in Mathematics, 704. - SpringerVerlag Berlin Heidelberg, New York, 1979. - P. 241-251.
24. Pereyra V., Keller H.B. Finite Difference Solution of Two-Point Boundary Value Problems: Preprint 69 / Dept. Math., Univ. - Southern California. - Los Angeles, 1976. - 130 p.
25. Pereyra V. High Order Finite Difference Solution of Differential Equations. - Stanford Univ. Comp. Sci. Report STAN-CS-73-348, 1973. - 86 p.
26. Lentini M., Pereyra V. An adaptive finite difference solver for nonlinear two point boundary problems with mild boundary layers // SIAM J. Numer. Anal. - 1977. - Vol. 14. - No. 1. - P. 91-111.
27. Красноруцкий Д.А. Развитие модели тонкого упругого стержня для расчета изгибно-крутильных колебаний авиационных лопастей // Наука. Промышленность. Оборона: тр. 13 Всерос. науч.-техн. конф.; Новосиб. гос. техн. ун-т. - Новосибирск, 2012. - С. 328-332.
28. Пустовой Н.В., Левин В.Е., Красноруцкий Д.А. Математическое моделирование контактного взаимодействия витков гибкого стержня при петлеобразовании // Прикладные задачи математики: материалы 21-й Междунар. науч.-техн. конф. (Севастополь, 16-20 сент. 2013 г.). - Севастополь: Изд-во Севастоп. нац. техн. ун-та, 2013. -С. 47-51.
29. Аргирис Дж. Современные достижения в методах расчета конструкций с применением матриц. - М.: Изд-во лит. по строительству, 1968. - 242 с.
30. Argyris J.H. An excursion into large rotations // Comp. Meth. Appl. Mech. Eng. - 1982. - Vol. 32. - No. 1. - P. 85-155.
31. Björck A., Pereyra V. Solution of Vandermonde Systems of Equations // Mathematics of computation. - 1970. - Vol. 24. - No. 112. -P. 893-903.
32. Пустовой Н.В., Левин В.Е., Красноруцкий Д.А. Применение геометрически нелинейных уравнений стержня к расчету статики и динамики тросов. Ч. 2 // Научный вестник Новосиб. гос. техн. ун-та. -2012. - № 2. - С. 106-116.
33. Ortega J.M. and Rheinboldt W.C. Iterative Solution of Nonlinear Equations in Several Variables. - New York: Academic Press, 1970. -572 p.
34. Deuflhard P. A Stepsize Control for Continuation Methods and its Special Application to Multiple Shooting Techniques // Mathematik. -1979. - P. 115-146.
References
1. Svetlitskii V.A. Mekhanika sterzhnei [Mechanics of rods]. Statika. Moscow: Vysshaia shkola, 1987. 320 p.
2. Zhilin P.A. Prikladnaia mekhanika. Teoriia tonkikh uprugikh sterzhnei: uchebnoe posobie [Applied mechanics. The theory of thin elastic rods]. Sankt-Peterburgskiy politekhnicheskiy universitet, 2007. 101 p.
3. Svetlitskii V.A. Stroitelnaya mekhanika mashin. Mekhanika sterzhnei [Construction mechanics of machines. Mechanics of rods]. Dinamika. Moscow: Fizmatlit, 2009. 383 p.
4. Kirchhoff G.R. Ueber das Gleichgewicht und die Bewegung einer elastischen Staben [About the equilibrium and motion of an elastic rods]. Crelle Journal fuer die reine und angewandte Mathematik, 1858, Bd. 56, pp. 285-313.
5. Kirkhgof G. Mekhanika [Mechanics]. Akademiia nauk SSSR. Moscow, 1962. 402 p.
6. Liav A. Matematicheskaia teoriia uprugosti [A treatise on the mathematical theory of elasticity]. Moscow; Leningrad: Otdelenie nauchno-tekhnicheskikh izdatelstv, 1935. 674 p.
7. Pustovoy N.V., Levin V.E., Krasnorutskiy D.A. Primenenie ge-ometricheski nelineinykh uravnenii sterzhnia k raschetu statiki i dinamiki trosov. Chast 1 [Applying of geometrically nonlinear equations of a rod to calculate statics and dynamics of cables. Part 1]. NSTU Scientific Bulletin, 2012, no. 1 (46), pp. 83-92.
8. Pustovoy N.V., Levin V.E., Krasnorutskiy D.A. Metodika vy-chisleniia parametrov bol'shikh povorotov poperechnykh sechenii gibkogo sterzhnia pri raschetakh v ramkakh ego differentsial'noi modeli. Chast 1 [The method of calculation of parameters of big rotations of cross-sections of a flexible rod using its differential model. Part 1]. NSTU Scientific Bulletin, 2013, no. 2 (51), pp. 155-164.
9. Krasnorutskiy D.A., Levin V.E., Pustovoy N.V. Nelineinaia di-namika tonkikh uprugikh sterzhnei [Nonlinear dynamics of thin elastic rods]. Nelineinye kolebaniya mekhanicheskikh sistem: trudy IX Vserossiy-skoy nauchnoi konferentsii (Nizhniy Novgorod, 24-29 sentyabrya 2012 g.). N. Novgorod: Nash dom, 2012, pp. 557-565.
10. Pustovoy N.V., Levin V.E., Krasnorutskiy D.A. Metodika vy-chisleniia parametrov bolshikh povorotov poperechnykh sechenii gibkogo sterzhnia pri raschetakh v ramkakh ego differentsial'noi modeli. Chast 2. [The method of calculation of parameters of big rotations of cross-sections of a flexible rod using its differential model. Part 2]. NSTU Scientific Bulletin, 2013, no. 3(52), pp. 146-159.
11. Sorokin F.D. Priamoe tenzornoe predstavlenie uravnenii bol'shikh peremeshchenii gibkogo sterzhnia s ispol'zovaniem vektora konechnogo povorota [Direct tensor representation of the equations large displacements of an elastic rod with using the vector of finite rotation]. MTT, 1994, no. 1, pp. 164-168.
12. Du H., Xiong W., Wang H., Wang Z., Yuan B. Nonlinear dynamic deformation simulation for helical rod like objects. Engineering Review, 2013, vol. 33, iss. 3, pp. 233-238.
13. Bathe K.J. Finite Element Procedures. Englewood Cliffs. NY: Prentice Hall, 1996. 1037 p.
14. Newmark N.M. A Method of Computation for Structural Dynamics. Journal of Engineering Mechanics Division, ASCE, July 1959, vol. 85, no. EM3, pp. 67-94.
15. Park K.S. An improved stiffly stable method for direct integration of nonlinear structural dynamic equations. Journal of Applied Mechanics, ASME, June 1975, vol. 42, iss. 2, pp. 464-470.
16. Shampine L.F., Muir P.H., Xu H. A User-Friendly Fortran BVP Solver. Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM), 2006, vol. 1, no. 2, pp. 201-217.
17. IMSL [electronic resource]: Fortran Numerical Library. User's Guide. Math Library. Version. 7.0, available at: http://www.rogue-wave.com/documents.aspx? entryid=519&command=core_download.
18. Pereyra V. Pasva3: An adaptive finite difference fortran program for first order nonlinear, ordinary boundary problems. Lecture Notes in Computer Science, 1979, vol. 76, pp 67-88.
19. Rashidinia J. Finite difference methods for a class of two-point boundary value problems. IUST International Journal of Engineering Science, 2008, vol. 19, no. 5-2, pp. 67-72.
20. Vainberg A.M. Matematicheskoe modelirovanie protsessov pere-nosa. Reshenie nelineinykh kraevykh zadach [Computer-aided simulation of transfer processes. Solving of a nonlinear boundary-value problems]. Moskow-Jerusalaem, 2009. 209 p.
21. Dinkar Sharma, Ram Jiwari, Sheo Kumar. Numerical Solution of Two Point Boundary Value Problems Using Galerkin-Finite Element Method. International Journal of Nonlinear Science, 2012, vol. 13, no. 2, pp. 204-210.
22. Lentini M. Boundary Value Problems over Semi-Infinite Intervals: Ph.D. Thesis, Cal. Inst, of Technology, 1978, 123 p.
23. Keller H.B. Constructive Methods for Bifurcation and Nonlinear Eigenvalue Problems. Lecture Notes in Mathematics, 704. Springer-Verlag Berlin Heidelberg, New York, 1979, pp. 241-251.
24. Pereyra V., Keller H.B. Finite Difference Solution of Two-Point Boundary Value Problems: Preprint 69. Dept. Math., Univ. Southern California, Los Angeles, 1976, 130 p.
25. Pereyra V. High Order Finite Difference Solution of Differential Equations, Stanford Univ. Comp. Sci. Report STAN-CS-73-348, 1973, 86 p.
26. Lentini M., Pereyra V. An adaptive finite difference solver for nonlinear two point boundary problems with mild boundary layers. SIAM J. Numer. Anal, 1977, vol. 14, no. 1, pp. 91-111.
27. Krasnorutskiy D.A. Razvitie modeli tonkogo uprugogo sterzhnia dlia rascheta izgibno-krutil'nykh kolebanii aviatsionnykh lopastei [Development of a thin elastic rod model for the calculation of flexural-torsional vibration of aircraft blades]. Trudy 13 Vserosiiskoi nauchno-tekhnicheskoi konferentsii "Nauka. Promyshlennost. Oborona". Novosibirskiy gosudarst-vennyi tekhnicheskiy universitet, 2012, pp. 328-332.
28. Pustovoy N.V., Levin V.E., Krasnorutskiy D.A. Matematicheskoe modelirovanie kontaktnogo vzaimodeistviia vitkov gibkogo sterzhnia pri petleobrazovanii. [Mathematical modeling of self-contact at a looping process of flexible rod]. Materialy 21 mezhdunarodnoi nauchno-tekhnicheskoi konferentsii "Prikladnye zadachi matematiki" (Sevastopol', 16-20 Sent. 2013 g.). Sevastopolskiy natsionalnyi tekhnicheskiy universitet, 2013, pp. 47-51.
29. Argiris Dzh. Sovremennye dostizheniia v metodakh rascheta kon-struktsii s primeneniem matrits [Recent developments in the methods of calculation of structures using matrixes]. Moscow: Izdatelstvo literatury po stroitelstvu, 1968. 242 p.
30. Argyris J.H. An excursion into large rotations. Comp. Meth. Appl. Mech. Eng, 1982, vol. 32, no. 1, pp. 85-155.
31. Björck A., Pereyra V. Solution of Vandermonde Systems of Equations. Mathematics of computation, 1970, vol. 24, no. 112, pp. 893-903
32. Pustovoy N.V., Levin V.E., Krasnorutskiy D.A. Primenenie ge-ometricheski nelineinykh uravnenii sterzhnia k raschetu statiki i dinamiki trosov. Chast 2 [Applying of geometrically nonlinear equations of a rod to calculate statics and dynamics of cables. Part 2]. NSTU Scientific Bulletin, 2012, no. 2, pp. 106-116.
33. Ortega J.M. and Rheinboldt W.C. Iterative Solution of Nonlinear Equations in Several Variables. New York: Academic Press, 1970. 572 p.
34. Deuflhard P. A Stepsize Control for Continuation Methods and its Special Application to Multiple Shooting Techniques. Mathematik, 1979, pp. 115-146
Об авторах
Пустовой Николай Васильевич (Новосибирск, Россия) - доктор технических наук, профессор, заведующий кафедрой «Прочность летательных аппаратов», ректор Новосибирского государственного технического университета (630073, г. Новосибирск, пр. К. Маркса, 20, e-mail: [email protected]).
Левин Владимир Евгеньевич (Новосибирск, Россия) - доктор технических наук, профессор, заместитель заведующего кафедрой «Прочность летательных аппаратов» Новосибирского государственного технического университета (630073, г. Новосибирск, пр. К. Маркса, 20, e-mail: [email protected]).
Красноруцкий Дмитрий Александрович (Новосибирск, Россия) - кандидат технических наук, доцент кафедры «Прочность летательных аппаратов» Новосибирского государственного технического университета (630073, г. Новосибирск, пр. К. Маркса, 20, e-mail: [email protected]).
About the authors
Nikolay V. Pustovoy (Novosibirsk, Russian Federation) - Doctor of Technical Sciences, Professor, Head of the Department of Aircraft Strength, rector of Novosibirsk State Technical University (20, K. Marks av., 630073, Novosibirsk, Russian Federation, e-mail: [email protected]).
Vladimir E. Levin (Novosibirsk, Russian Federation) - Doctor of Technical Sciences, Professor, Deputy head of the Department of Aircraft Strength of Novosibirsk State Technical University (20, K. Marks av., 630073, Novosibirsk, Russian Federation, e-mail: [email protected]).
Dmitry A. Krasnorutskiy (Novosibirsk, Russian Federation) -Ph. D. in Technical Sciences, Associate Professor of the Department of Aircraft Strength of Novosibirsk State Technical University (20, K. Marks av., 630073, Novosibirsk, Russian Federation, e-mail: [email protected]).
Получено 04.04.2014
Просьба ссылаться на эту статью в русскоязычных источниках следующим образом:
Пустовой Н.В., Левин В.Е., Красноруцкий Д. А. Алгоритм численного решения нелинейной краевой задачи динамического деформирования тонкого стержня // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2014. - № 2. - С. 168-199.
Please cite this article in English as:
Pustovoy N.V., Levin V.E., Krasnorutskiy D.A. The numerical algorithm for solving nonlinear boundary problem of thin rod's dynamic deformations. PNRPUMechanics Bulletin. 2014. No. 2. P. 168-199.
