Эффективные методы численного моделирования динамики нелинейных систем абсолютно твёрдых и деформируемых тел (1105385), страница 3
Текст из файла (страница 3)
Производная потенциальнойэнергии по u даёт обобщённую силу ∂П ∂u = − Fuelast (u) .Виртуальная работа внешних сил δW' складывается из работ распределённых сил f(M) и моментов m(M) на возможных перемещениях точек тела.Введём столбец линейных и угловых скоростей& ~&ξ = ⎧⎨ r − ρω + v u ⎫⎬ .ω + ωu ⎭⎩Первая его составляющая выражает собой теорему о сложении скоростей, вторая – угловых скоростей.
Переходя по аналогии к вариациям и учтязависимости (1.3), (1.4), получим выражение для возможных линейных иугловых перемещений через вариации обобщённых координат:~ δω + δv ⎫ ⎧ D δq − ρ~ B δq + D δu ⎫⎧ δr − ρuuδξ = ⎨⎬ = CΦ δq + S δu ,⎬=⎨B δq + B u δu ⎭δω + δω u ⎭ ⎩⎩(1.8)где введены матрицы~⎤⎡I −ρC=⎢⎥,⎣O I ⎦⎡D⎤Φ = ⎢ ⎥,⎣B⎦⎡D ⎤S = ⎢ u⎥.⎣ Bu ⎦Теперь вычислим виртуальную работу сил и моментов F = {fT, mT}T на2T⎡ T⎤ввиду тождества d ⎢ ∂v ⎥ v = ∂v v , которое легко проверяется непосредственно.∂xdt ⎣ ∂x& ⎦15возможном перемещении δξ:δ W ′ = ∫ δ ξ T F dV = ∫ ( δ q T Φ T C T F + δ u T S T F ) d V ,VVа также столбцы обобщённых активных силδW ′= ΦT ⎛⎜ ∫ C T F dV ⎞⎟ = ΦT Fq ,⎝ V⎠δqδW ′= ∫ S T F dV .Q′u =VδuQ′q =(1.9)Аналогично получим обобщённые силы Q ′′ от сосредоточенных сил имоментов реакций Rk (Рис.
1.1), заменив интегрирование по объёму суммированием по числу реакций:Q′q′ = Φ TQ′u′ =(∑ C(∑ SkkTMkTMk)R k = ΦT R q ,)(1.10)Rk = Ru.В этих формулах матрицы C и S вычислены в тех точках тела Mk, в которых действуют реакции Rk.При подстановке выражений (1.6)–(1.10) уравнения (1.5) принимают вид~ ω~ω~ mρ + 2 µω~ ρ& dV ⎫⎞⎛ ⎧ m&r& − mρ& + ∫ µρ&& dV + ωCC∫⎟⎪⎜⎪VVTFRΦ ⎜⎨ ~−−⎬qq⎟ = 0,~ρ~ J ω + 2 µρ~ω~ ρ& dV&&&&&mµdVρrJωρω+++⎟⎜⎪ CCC∫V∫V⎪⎭(1.11)⎠⎝⎩∫V D u ( µ&r& − µρ~ω& + µρ&& + µω~ ω~ ρ + 2 µω~ ρ& ) dV − Q′u − FuTelast− Ru = 0 .Здесь m = ∫ µ dV – масса тела, mρ C = ∫ µρ dV – радиус-вектор центра массVV~ρ~ dV – тентела относительно точки O1, умноженный на массу, J C = − ∫ µρVзор инерции тела относительно точки O1.Чтобы вынести вторые производные из-под знаков интегралов, выразим&& через набор обобщённых координат u и их производные (эту процедуруρиногда называют дискретизацией [55]):& u u& = D u u&& = D u u&& + D&& + a′u′ .ρ(1.12)16Матрица Φ является квадратной (6×6) и невырожденной3, поэтому послеумножения первого уравнения системы (1.11) слева на Φ–T получим окончательные уравнения движения отдельного тела в виде:&& + k q = Fq + R q ,M q w + M qu u&& + k u = Fu + R u .M quT w + M u u(1.13)Использованные обозначения [13]:⎧ &r& ⎫w = ⎨ ⎬,&⎭⎩ω⎡ mEM ( q, u ) = ⎢~⎣m ρCqFu = Q′u + Fuelast ,⎡ ∫ µD u dV ⎤~ T⎤mρVquC⎥ , M u ( u) = ∫ µ DTu D u dV ,⎥ , M ( q, u ) = ⎢~D dV ⎥V⎢ ∫ µρJC ⎦u⎦⎣V~~~&⎧ω⎫′′⎪ ωmρ C + 2 ∫V µωρ dV + ∫V µa u dV ⎪T~~~kq = ⎨~⎬ , k u = ∫V D u µ (ωωρ + 2ωρ& + a′u′ ) dV .~~~⎪⎩ω J C ω + 2 ∫V µ ρ ωρ& dV + ∫V µ ρ a ′u′ dV ⎪⎭В работе [55] уравнения (1.13) получаются на основе принципа возможных перемещений.
Там же отмечено, что в случае «отвердевания» КЭ, когдаρ& (M ) ≡ 0 , и при записи уравнений относительно центра масс, когда ρ C ≡ 0 ,они превращаются в динамические уравнения Ньютона-Эйлера для абсолютно твёрдого тела:m&r& = f акт + f реакт ,~ J ω = mакт + m реакт .& +ωJ ωC1.1.3.(1.14)CПрямой метод формирования уравнений движениясистемы телВ п. 1.1.2 были сформированы уравнения движения отдельного абсолютно твёрдого или деформируемого тела.
Здесь выводятся уравнения движения системы таких тел.Уравнения (1.13) иногда называются полудискретизированными [55],3Кроме отдельных положений, в которых углы ориентации вырождаются. Как отмечено в начале этого параграфа, для избежания вырождения можно использовать четыре параметра Эйлера.17&& были выраженычтобы подчеркнуть, что лишь относительные ускорения ρ&& , формула (1.12).через обобщённые ускорения uВыразим столбец wi ускорений каждого тела через его обобщённые координаты qi и производные от них, подобно тому, как в выражениях (1.3):&& + a′i′⎫⎧ &r& ⎫ ⎧D q&&i + w′i′ .wi = ⎨ i ⎬ = ⎨ i i⎬ = Φi q& i ⎭ ⎩ Bi q&&i + ε′i′ ⎭⎩ω(1.15)Если теперь подставить это значение wi в уравнение (1.13) для i-го телаи выполнить суммирование по всем телам системы, то получим уравнениядвижения системы в обобщённых координатах⎛∑ ⎜⎜i =1 ⎝n⎡ΦiT M iq Φi⎢O⎣O⎤ ⎡ O⎥+⎢O⎦ ⎣M iquT Φi&& i ⎫ n ⎧ΦiT Q iq ⎫ΦiT M iqu ⎤ ⎞⎟ ⎧q⎥ ⎨ ⎬ = ∑ ⎨ u ⎬ , (1.16)&& i ⎭ i =1 ⎩ Qi ⎭M iu ⎦ ⎟⎠ ⎩uгде обозначено:Q iq = Fiq − k iq − M iq w ′i′,Q iu = Fiu − k iu − ∫ D iuT µ i (a′i′ + ~εi′′ρ i ) dVi .ViВ краткой матричной записи уравнения (1.16) можно записать в виде[M q + M qu ] &x& = Q ,(1.17)где x – глобальный столбец обобщённых координат системы, являющийсяобъединением всех столбцов x i = {q iT , u iT }T с учётом глобальной нумерациикоординат (см.
приложение 6.3). Заметим, что силы реакций в эти уравненияне входят, что является следствием идеальности связей в шарнирах.Отметим, что все элементы уравнений движения (1.16) вычисляются сиспользованием только алгебраических матричных операций умножения исложения. Матрицы производных Di, Bi, входящие в соотношения (1.15) идалее через матрицы Φi в (1.16), можно вычислить без использования дифференцирования, используя рекуррентные соотношения для цепочек тел, см.приложение 6.6. У истоков описанного формализма стояли Роберсон и Виттенбург [45], Вукобратович [64], а также Шилен и Кройцер [48].Оценим вычислительную сложность формирования уравнений движенияв этом методе на примере системы в виде цепочки из n тел.
Число алгебраи-18ческих операций для вычисления (трудоёмкость вычисления) составляющихMq и Mqu матрицы масс системы (1.17) равно, соответственно O(n3) и O(n2),как несложно понять по их структуре в уравнении (1.16)4. Трудоёмкость вычисления глобального вектора обобщённых сил также равна O(n2). В последующих пунктах будут рассмотрены более эффективные алгоритмы.1.1.4.Метод составных телПовысить эффективность вычислений можно за счёт более глубокогоанализа внутренней структуры матрицы масс [24].
При этом вновь рассматривается система в виде цепочки n тел5. Матрица масс из (1.16) примет видM i = ΦiT M iq Φi .(1.18)Рассмотрим соотношения между линейными и угловыми скоростямидвух смежных тел цепочки, имеющие вид, аналогичный виду (1.8)Vi +1 = Ci ,i +1Vi + S il +1q& i +1 + Vi′+1 ,(1.19)~⎧ r&i ⎫⎡I −ρi ,i +1 ⎤Vi = ⎨ ⎬ , C i ,i +1 = ⎢– матрица преобразования, ρ i ,i +1 –I ⎥⎦⎩ω i ⎭⎣Oгдерадиус-вектор между центрами масс тел i и (i – 1); S il +1 и q& i +1 – локальнаяматрица Якоби в шарнире (i + 1) и столбец локальных шарнирных скоростей.Анализируя соотношения между скоростями в шарнирах цепочки тел,можно получить замечательные соотношенияΦi =i∑ C k ,i Skl ,k =1[Skl = O ... O S lkO ...
O]для рекуррентного вычисления матриц Φi . Здесь матрица Skl содержит единственную ненулевую подматрицу S lk , расположенную на k-том месте.В конечном итоге переход к рекуррентным соотношениям позволяетснизить трудоёмкость вычисления матрицы масс (1.18) до O(n2) и столбцаобобщённых сил Q до O(n) за счёт изменения порядка суммирования.4Трудоёмкость каждого умножения на матрицу Фi пропорциональна n, так как она имеет размер 6 × nи, вообще говоря, является плотной.
Кроме того, в формуле производится суммирование по n телам.5Здесь для простоты рассматривается цепочка из абсолютно твёрдых тел.191.1.5.Метод отдельных телОтметим, что вычислить матрицу масс за меньшее, чем n2, число операций, невозможно, так как она сама содержит n2 элементов.
Однако Верещагиным [4] был предложен метод формирования уравнений движения сложности O(n), в котором глобальная матрица масс вообще не вычисляется. История его появления6 описывается в работе [49]. Первое удачное применениеэтого формализма принадлежит Айхбергеру в 1993-1994 гг. [22, 23].Изложим суть метода в трактовке работы [77]. Рассмотрим два концевыетела цепочки, изображённые на Рис. 1.2.RnnRn–1nn–1ρnQnn–1Qn–1R* nРис. 1.2. Концевые тела цепочки тел и силы, действующие на нихЗапишем уравнения движения концевого тела n. На него действуютобобщённые силы Qn и силы реакции Rn в шарнире n (приведенные к центрумасс), ср. с уравнениями (1.13):M n w n = Qn + R n .(1.20)Запишем далее то же уравнение для предыдущего тела в цепочке (n – 1)с учётом того, что на него действуют, кроме активных сил Qn–1 и сил реакцийRn–1 в шарнире (n – 1), ещё и реакции R*n со стороны шарнира n:6Развитие этого формализма, как описывают Schwertassek и Rulka в работах 1991 года, не относится квесьма славным эпизодам в истории динамики систем тел.
Метод был опубликован Верещагиным в 1974году, но важность этой работы не была оценена. В середине восьмидесятых формулировка алгоритма былапереоткрыта различными авторами, и только тогда был сделан шаг по его применению к построению эффективных методов моделирования динамики систем тел.20⎡IR ∗n = ⎢~⎣ρnO⎤⋅ (− R n ) = −C Tn R n .⎥I⎦{При приведении столбца R n = f nTm Tn}T сил реакций f(1.21)nи их моментов mn кцентру масс тела (n – 1) необходимо добавить момент от этих сил относительно центра приведения, что и выражается формулой (1.21).Итак, для тела (n – 1) запишем:M n −1w n −1 = Q n −1 + R n −1 − C Tn R n .(1.22)Запишем кинематическое соотношение между ускорениями двух тел,которое является продолжением серии уравнений (1.8) и (1.19):&& n + w ′n .w n = C n w n −1 + S n q(1.23)В случае идеальных связей выполняется соотношениеS Tn R n = 0 ,(1.24)оно вытекает из условия равенства нулю работ на возможном перемещении.Используем условие (1.24) для исключения сил реакций из уравнения(1.20), умножив его слева на S Tn :S Tn M n w n = S Tn Q n .Подставим в это выражение вместо wn его значение из (1.23)&& n + w ′n ) = S Tn Q n .S Tn M n (C n w n −1 + S n qИз последнего уравнения получим важную формулу, связывающую ускорения тела и вторую производную по времени от координат в шарнире:&& n = U n−1S Tn {Q n − M n (C n w n −1 + w ′n )}.q(1.25)Здесь матрица Un имеет видU n = S Tn M n S n .(1.26)Размер этой матрицы равен числу n степеней свободы в шарнире.
Она симметрична и положительно определена. В данном методе требуется обращение лишь подобных ей матриц. И поскольку размер их всегда мал (он не превышает 6), то этим и обусловлена эффективность описываемого метода.Теперь мы можем выразить силы реакции Rn из (1.20) через ускорения21wn, выполнив подстановку выражения (1.25) в (1.23). Так мы исключим реакции Rn из уравнения (1.22) и преобразуем к видуM ∗n −1w n −1 = Q∗n −1 + R n −1 .В этой формуле введены обозначения (n заменено на k)M ∗k −1 = M k −1 + C Tk M k C k − C Tk M k S k U k−1S Tk M k C k ,{(1.27)}Q∗k −1 = Q k −1 − C Tk M k S k U k−1S Tk ( Q k − M k w ′k ) + w ′k − Q k ..(1.28)Таким образом, начав с уравнения (1.20) для тела n, после некоторыхпреобразований мы пришли к сходному уравнению (1.27) для тела (n – 1).Преобразованная матрица масс M ∗n −1 уже не будет блочно-диагональной, нопо-прежнему будет оставаться симметричной и положительно определённой.Теперь мы можем записать уравнение, аналогичное (1.27), для телаk = n – 1, затем для тела k = n – 2 и так далее вплоть до тела k = 1, вычисляяматрицы M ∗k −1 и векторы Q∗k −1 .
При этом в (1.28) вместо величин Mk и Qkбудут входить их преобразованные значения из предыдущего шага рекурсии.Так реализуется обратный ход алгоритма от конца цепочки к её началу.Прямой ход – от начала цепочки (тела 1) до её конца – реализуется с ис-пользованием формул (1.25), (1.23), k = 1,…,n.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
