Декомпозиция конечно-элементной модели

Размерность задач, для решения которых применяется метод конечных элементов, имеет устойчивую тенденцию к росту. Увеличение быстродействия вычислительных систем решает только часть появляющихся в этой связи проблем. Приоритетным направлением в обеспечении эффективного использования возможностей вычислительной техники является развитие методов декомпозиции конечно-элементных алгоритмов. Обзоры работ в данной области исследований, актуальные проблемы и новые методы их решения представлены в [4], [7]. Для решения задач определенного класса может быть использован предлагаемый далее метод декомпозиции конечно-элементной модели.

Запишем в векторно-матричной форме уравнение движения конечно-элементной модели некоторой механической системы:

MU + KU + RU = P + CN . (1)

Здесь U , U и U - векторы (одномерные массивы) перемещений, скоростей и ускорений узлов соответственно; R , K и M - матрицы жесткости, демпфирования и масс соответственно; P - вектор внешних воздействий (внешних сил). В терминах метода перемещений строительной механики соотношение (1) можно интерпретировать как систему уравнений равновесия конечно-элементной модели [5; 10-15]. Уравнение (1) отличается от общ епринятой формулировки [6] компонентом CN , где C - матрица коэффициентов в уравнениях равновесия. Компонентами вектора N являются силы и моменты (пары сил), появляющиеся в физической модели как результат контактного взаимодействия конечных элементов. Векторы N и U заранее не известны и подлежат определению как результат решения задачи.

Результатом внешнего воздействия и указанного выше контактного взаимодействия являются, в частности, деформации конечных элемен- тов, линейные и угловые перемещения в узлах. Эти перемещения ведут к появлению соответствующих зазоров, если в физической модели нет связей, которые не допускали бы изменений линейных и угловых расстояний между каждой парой узлов контактирующих конечных элементов. Текущие и начальные значения зазоров (D и D0 соответственно) связаны с перемещениями узлов U соотношением

D = C T и + D ₀. (2)

Из физического смысла контактного взаимодействия в узлах конечных элементов следует, что каждый столбец матрицы C содержит только два ненулевых элемента, один из которых равен 1, другой - (-1). Заметим, что существует класс механических систем, в которых при некоторых воздействиях возможно появление зазоров, что рассмотрено, например, в [2], [3]. Далее рассматривается задача, по физическому смыслу которой зазоры равны нулю, то есть D = 0 , D ₀ = 0 в формуле (2).

Дискретизация уравнения движения (1) по времени с использованием явной схемы с односторонними конечными разностями приводит на шаге i с учетом равенства (2) к следующим соотношениям [2]:

AU , = P , + CN ; C T U , = 0 , (3) где A = M r ^{- 2} + K t ^{- 1} + R , P i = P ■ т ^-2 MU _{i - 2} + (2 M r ^-2 + + K t "¹) U i _ ₁ , т - шаг по времени, i - номер шага.

Пусть конечно-элементная модель разбита на n подструктур. Каждая подструктура содержит m конечных элементов, m >  1. Рассматривая независимо от других каждую подструктуру к , к = 1, .., n , сформируем по аналогии с (3) матрицы A , ^{( к} ) , C ,^к ) , вектор P '^к ) . Не указывая индекс i , запишем соотношения вида (3) для подструктуры к . Получим:

A < к ) и ( к ) = P ( к ) ₊ C T ( к ) N ; C к ) и ⁽ к ) = 0

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

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

U ⁽ k ) = A ^{- 1(} k ) ( P ⁽ k ) + C ⁽ k > N ).               (5)

Касаясь вопроса о существовании указанной в (5) обратной матрицы, отметим, что система с распределенными параметрами под действием движущейся сосредоточенной массы рассмотрена в статье [6], автором которой, в частности, доказано, что матрица коэффициентов при старших производных в уравнении вида (1) всегда обратима.

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

C TU = 2k=1 CT(k 'U1k) = 0, получим:

• 2 k = 1 C T ⁽ k ) a ^{- 1(} k ) ( P ⁽ k ) + C ⁽ k ) N ) = 0 .

Отсюда находим:

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

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

Из физического смысла контактного взаимодействия конечных элементов (или подс т руктур) следует, что каждый столбец матрицы C содержит только два ненулевых элемента, один из которых равен 1, другой – (-1). При этом каждый столбец блока C ^{( k )} содержит только один ненулевой элемент, равный +1 или -1. Поэтому матрица H = 2 k = i ( C T ⁽ k ) A ^-1( k ) C ⁽ k ) ) является диагональной. Учет особенностей структуры блока C ^{( k )} ведет к существенному упрощению вычислений по формулам (6) и (5).

Применение представленного метода декомпозиции конечно-элементной модели в алгоритмах расчета позволяет уменьшить затраты времени и объем требуемой оперативной памяти, что особенно важно при итерационном решении нелинейных задач [1], [3]. Кроме того, создаются новые возможности для параллельных вычислений.