Алгоритм декомпозиции конечно-элементной модели механической системы с односторонними контактами

Развитие методов декомпозиции обусловлено необходимостью решения поставляемых практикой все более сложных задач, требующих компьютерного моделирования инженерных объектов и технологических процессов. Цель разработки метода декомпозиции - корректно заменить решение одной большой задачи решением серии меньших задач и тем самым уменьшить затраты времени как в последовательных, так и в параллельных алгоритмах численного моделирования. Применение методов декомпозиции позволяет расширить класс объектов, моделирование которых может быть выполнено с применением существующей компьютерной техники. Появление многопроцессорных вычислительных систем и мультиядерных процессоров стимулировало эволюцию методов декомпозиции. Развитие методов декомпозиции привело к формированию одного из фундаментальных направлений вычислительной математики. К актуальным проблемам данного направления относится совершенствование алгоритмов численного моделирования с применением метода конечных элементов [6], [7]. В данной статье в развитие работы [3] рассматривается алгоритм декомпозиции конечно-элементной модели механической системы с полукоэрцитивны-ми (односторонними) контактами * . Примером такой модели может быть система «контактная сеть - токоприемники электроподвижного состава» [1], [4]. Объектом исследования является математическая модель механической системы с полукоэрцитивными (односторонними) контактами. Модель получена как физически обоснованное описание контактного взаимодействия подсистем механической системы. Предмет изучения - вычислительная эффективность алгоритма компьютерной реализации модели.

Пусть конечно-элементная модель разбита на п подструктур. Запишем в векторно-матричной форме уравнения движения каждой подструктуры k , к = 1, .„, п.

M ⁽ k ) U ⁽ k ) + K ^{( k} ) U ^{( k} ) + R ^{( k} ) U ^{( k} ) = P ^{( k} ) + C ^{( k} ) N . (1)

Здесь U ^{( k} ) , U ^{( k} ) и U ^{( k} ) - векторы (одномерные массивы) соответственно перемещений, скоростей и ускорений узлов; R ^{( k} ) , K ^{( k} ) и M ^{( k} ) - матрицы соответственно жесткости, демпфирования и масс. P ^{( k )} - вектор внешних воздействий (внешних сил). Соотношение (1), которое можно интерпретировать как систему уравнений равновесия, отличается от общепринятой формулировки слагаемым C ^{( k} ) N , где C ^{( k} ) - матрица коэффициентов в уравнениях равновесия. Элементами вектора N являются силы и моменты (пары сил), появляющиеся при контакт н ом взаимодействии подсистем (см. рисунок). N и U ^{( k} ) заранее не известны, их определение - цель решения задачи.

N yba a

N y ab )

M^ab )

-Nxab )

b

M< ba )   ^{C (} b ) ^N =

N X ba )

_ _ f1

C ⁽ a ) N = 0

V 0

_— C ⁽ a )

1J

f N Xab ⁾) N yab) M' ab ’

V J

fN))

N ab ) M ( ab )

V J

Узлы а и b двух подсистем

Если, например, N (^ab ) = 0 и M^{( ab} ) = 0, то C ⁽ a ) N = = 1 - N y^ab ) , C ^{( b} ) N = - 1 - N yy,^ab ) . Это модель скользящего вдоль оси x шарнирно-подвижного контакта.

Зазоры укажем в векторе D , определив их линейными и угловыми компонентами. Например, для двух контактирующих узлов а и b получим:

D ⁽ ab ) = C a ) T u ⁽ a ) + C b ) T u ⁽ b ) + d 0 ab )

Не указывая верхние индексы, запишем дискретный аналог уравнения (1), используя явную схему с односторонн и ми конечными разностями. С учетом (2) при D ₀ = 0 получим [2]:

^

AU _Z = P + CN ; D _z = C T U _z ,         (3)

где: A = M т - ² + K t - ¹ + R ,

P = P + т - MU - 2 + (2 M t ^{- 2} + K t ^{- 1} ) U _Z -_b

τ – шаг по времени, – номер шага.

Для подструктуры k на шаге получим (индекс не указан):

a ⁽ k ) u ⁽ k ) _— C ⁽ k ) n = p ⁽ k ) ; d ⁽ k ) = c ⁽ k ) T u ⁽ k )

Систему этих соотношений для всех k = 1, …, n запишем в виде блочного векторно-матричного равенства:

' A (1)	0          •••	0	- c (1)	^		" U <1) >		f P <1) ^
0 • • •	a (2)      ••• •••           •••	0 •••	- c (2) •••			u '2' •••	—	P < 2) •••	. (4)
0	0          •••	A ( n )	- C ( n )			U < n )		P <  n )
(1) T V C	- c (2) T   •••	- C ( n ) T	0	>		V N )		-D V D )

Из (4) следует:

U ⁽ k ) = A ^- 1⁽ k ) ( p ⁽ k ) + c ⁽ k ) N ).             (5)

Используя (5) и учитывая, что, согласно (4),

C TU = 2k=1 C( k) T U( k) = D, получим:

2 k = 1 C T ⁽ k ) A ^{- 1(} k ) ( P ⁽ k ) + C ⁽ k ) N ) = D .

Если односторонних связей нет, то D = 0 . Тогда [3]

N = — ( e k = 1 C T ⁽ k ) A ^{- 1(} k ) C ⁽ k ) ) — ¹ 2 k = 1 C T ⁽ k ) A ^{- 1(} k ) P ⁽ k ) . (6)

Подставив (6) в (5), определим U ^{( k )} , k = 1, …, n . Таким обр а зом, задача определения векторов N и U при D = 0 решена. Если связи од н осторонние, то D * о и определение векторов N и U сводится к линейной задаче дополнительности:

D = 2 k = 1 C T ⁽ k ) A ^{- 1(} k ) C ⁽ k ) N + 2 k = 1 C T ⁽ k ) A ^{- 1(} k ) P ⁽ k ) ,

N >  0 , D >  0 , N T D = 0 .

Решив эту задачу, например, с применением энергетического критерия перехода односторонних св язе й в действительное состояние [5], найдем N и D . Затем определим U ^{( k )} ( 5), k = 1, …, n .

Тестирование алгоритма при D = 0 выполнено на компьютере с процессором Intel® Core™2 Duo CPU E8400@3.00 GHz, кэш L1 32 Kb, L2 6144 Kb, память DDR3 667 MHz 2 Gb, Intel P43, ОС Windows XP 32 bit. Модель механической системы [4] разбивалась на n подсистем с примерно одинаковым числом неизвестных. Затраты времени (в секундах) при решении задачи с общим числом уравнений N ₀ и числом подструктур n приведены в таблице.

N 0	n
	1	2	3	4	5	6	7
1000	0,36	0,26	0,15	0,11	0,10	0,09	0,09
2000	2,21	1,98	0,81	0,57	0,41	0,34	0,29
3000	6,52	9,08	3,14	1,51	1,11	0,91	0,72

При N ₀ > 1 выполнялась только декомпозиция, алгоритм вычислений – последовательный. Алгоритм параллельных вычислений обеспечит дальнейшее уменьшение времени счета.

Таким образом, применение представленного метода декомпозиции в задачах определенного класса [2], [4] позволяет уменьшить затраты времени при численном моделировании.

ПРИМЕЧАНИЕ

* В данной статье исправлены опечатки, допущенные в [3] по технической причине.