Перейти до змісту

Віртуальна робота зовнішніх сил і складання глобальних рівнянь

У розділі Дискретизація віртуальної роботи внутрішніх сил ліву частину слабкої форми зведено до вектора внутрішніх сил елемента \(\boldsymbol{q}^e\) (метод UL) або \(\boldsymbol{Q}^e\) (метод TL). У цьому розділі з віртуальної роботи зовнішніх сил вводиться вузловий вектор зовнішніх сил елемента \(\boldsymbol{F}^e\), після чого застосовується операція складання, яка перевпорядковує та сумує вузлові величини елементів за глобальними номерами вузлів, у результаті чого отримується нелінійна система рівнянь відносно вузлових переміщень, яку розв’язує нелінійний структурний аналіз FrontISTR.

Поелементний розклад віртуальної роботи зовнішніх сил

Праву частину рівняння з розділу Принцип віртуальної роботи можна розкласти по елементах на віртуальну роботу зовнішніх сил, що складається з об’ємних сил і заданих поверхневих навантажень на механічних границях. Щоб записати інтерполяцію переміщення, введену в розділі Функції форми та скінченно-елементна апроксимація, у матричній формі, для вузла \(\alpha\) розташуємо функцію форми \(N_\alpha^e\) на діагоналі блока \(d \times d\) \(\boldsymbol{N}_\alpha\) і сформуємо з таких блоків матрицю \(\boldsymbol{N} = [\boldsymbol{N}_1, \ldots, \boldsymbol{N}_{n_e}]\), так що \(\delta\boldsymbol{u} = \boldsymbol{N}\, \delta\boldsymbol{u}^e\). Підстановка цього виразу у віртуальну роботу зовнішніх сил, записану в опорній конфігурації, дає

\[ \delta W^{\mathrm{ext}} = \sum_e \delta\boldsymbol{u}^{eT} \boldsymbol{F}^e, \qquad \boldsymbol{F}^e_\alpha = \int_{\Omega^e_0} \boldsymbol{N}_\alpha^T \rho_0 \boldsymbol{g}\, dV + \int_{\Gamma^e_{0t}} \boldsymbol{N}_\alpha^T \bar{\boldsymbol{t}}_0\, d\Gamma_0 \]

де вузловий вектор зовнішніх сил елемента впорядковано як \(\boldsymbol{F}^e = (\boldsymbol{F}^{eT}_1, \ldots, \boldsymbol{F}^{eT}_{n_e})^T\). Таким чином, віртуальна робота зовнішніх сил зводиться до тієї самої форми «вузловий вектор елемента × тестова функція», що й частина внутрішніх сил (при записі в поточній конфігурації ту саму форму отримуємо замінами \(dV \to dv\), \(\rho_0 \to \rho\) і \(\bar{\boldsymbol{t}}_0 \to \bar{\boldsymbol{t}}\)).

Складання вузлових величин елементів

Вузлові величини \(\boldsymbol{Q}^e_\alpha, \boldsymbol{F}^e_\alpha\), отримані для кожного елемента, сумуються у глобальні вектори, впорядковані за глобальним номером вузла. Нехай глобальний номер вузла, що відповідає локальному номеру вузла \(\alpha\) елемента \(\Omega^e\), дорівнює

\[ \mathrm{gdx}(e, \alpha) = i_g \]

Тоді вузлова величина елемента збігається з відповідною компонентою глобальної вузлової величини (наприклад, \(\boldsymbol{u}^e_\alpha = \boldsymbol{u}_{i_g}\)). Оскільки вузол \(i_g\) зазвичай спільний для кількох елементів, визначимо множину пар \((e, \alpha)\), глобальний номер вузла яких дорівнює \(i_g\), як

\[ \mathcal{E}(i_g) = \{ (e, \alpha) \mid \mathrm{gdx}(e, \alpha) = i_g \} \]

Переписавши за допомогою цієї множини суму як \(\sum_e \sum_\alpha = \sum_{i_g} \sum_{(e,\alpha) \in \mathcal{E}(i_g)}\), отримуємо вузлову внутрішню силу та глобальний вектор внутрішніх сил для всіх \(n_g\) вузлів:

\[ \boldsymbol{Q}_{i_g} = \sum_{(e,\alpha) \in \mathcal{E}(i_g)} \boldsymbol{Q}^e_\alpha, \qquad \boldsymbol{Q} = (\boldsymbol{Q}^T_1, \ldots, \boldsymbol{Q}^T_{n_g})^T \]

Тут \(\boldsymbol{Q}_{i_g}\) є результуючою вузлових внутрішніх сил елементів, що діють у вузлі \(i_g\), і дорівнює \(\boldsymbol{0}\), коли зовнішня сила відсутня та виконується рівновага. Та сама процедура в методі UL дає \(\boldsymbol{q}_{i_g}, \boldsymbol{q}\); оскільки їхні числові значення задовольняють \(\boldsymbol{q} = \boldsymbol{Q}\), нижче використовується позначення \(\boldsymbol{Q}\), крім випадків, коли потрібне розрізнення. Глобальний вектор зовнішніх сил \(\boldsymbol{F}\) отримується таким самим накопиченням.

У реалізації множина \(\mathcal{E}(i_g)\) явно не будується; натомість внески додаються до відповідних компонентів усередині циклу по елементах.

Ініціалізувати глобальний вектор внутрішніх сил Q нулем: Q_{i_g} = 0  (i_g = 1, ..., n_g)
for e = 1 to (кількість елементів)
    for α = 1 to n_e
        i_g = gdx(e, α)
        Q_{i_g} += Q^e_α
    end for
end for

Глобальний вектор зовнішніх сил \(\boldsymbol{F}\) будується за тією самою процедурою. Операція додавання та збереження вузлових величин елементів у вектори й матриці, пронумеровані глобальними номерами вузлів, називається складанням. Для тензорів другого порядку, пов’язаних із двома номерами вузлів (наприклад, матриць жорсткості), аналогічне складання виконується з використанням множини \(\mathcal{E}^2(i_g, i_h) = \{ (e, \alpha, \beta) \mid \mathrm{gdx}(e, \alpha) = i_g\ \mathrm{and}\ \mathrm{gdx}(e, \beta) = i_h \}\) (конкретну побудову див. у розділі Дотична матриця жорсткості).

Нелінійні рівняння, що розв’язуються

Підставляючи складені внутрішні та зовнішні сили у рівняння з розділу Принцип віртуальної роботи і використовуючи той факт, що він виконується для будь-якої тестової функції \(\delta\boldsymbol{u}^n\), яка задовольняє геометричні граничні умови, отримуємо

\[ \boldsymbol{Q}(\boldsymbol{u}^n) - \boldsymbol{F}(\boldsymbol{u}^n) = \boldsymbol{0} \]

У контексті покрокового аналізу (Схема інкрементного аналізу) відновимо часовий нижній індекс \(_{n+1}\) і опустимо верхній індекс \(^n\), що позначає глобальний вузловий вектор. Тоді рівняння, яке потрібно розв’язати, має вигляд

\[ \boldsymbol{Q}(\boldsymbol{u}_{n+1}) - \boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{0} \]

Отже, дискретизована крайова задача знаходження вузлового переміщення \(\boldsymbol{u}_{n+1}\) у момент часу \(t_{n+1}\) зводиться до розв’язання цього нелінійного рівняння відносно переміщення разом із геометричними граничними умовами. Лінеаризацію рівняння та побудову дотичної матриці жорсткості описано в розділі Дотична матриця жорсткості, а ітераційний метод розв’язання — у розділі Метод Ньютона—Рафсона.

Пов’язані теми