Побудова результуючої системи лінійних алгебраїчних рівнянь

Підставляючи розклад (22) у слабку форму рівнянь Гальоркіна (24), отримаємо СЛАР відносно невідомих вузлових значень , яка в мат­ричній формі набуває вигляду

, (27)

з коефіцієнтами

, . (28)

Матриця системи (27) називається матрицею жорсткості, а права час­тина - вектором навантаження.

Враховуючи співвідношення (25), бачимо, що у виразі для коефі­цієн­тів матриці жорсткості (28) лише значення індекса рівні , та дають відмінні від нуля значення самих коефіцієнтів . Причому ін­тег­­рал буде містити лише доданки та , що відповідають вкладам СЕ з номерами та , відповідно. Це означає, по-перше, що матриця жорсткості буде трьохдіагональною і симетричною, а, по-друге, що значення коефіцієнта буде складатися з двох доданків: інтегралу , який є внеском СЕ з номе­ром , та інтегралу , який є внеском СЕ з номе­ром . Тому, на практиці, обчислення коефіцієнтів матриць СЛАР (27) здійснюють не безпосередньо за формулами (28), а шляхом обчислення локальних матриць жорсткості та вектора на­ван­та­жень з наступним рознесенням їх у глобальні матрицю жорсткості та вектор навантаження . Такий підхід до формування результуючої СЛАР в МСЕ у науковій літературі отримав назву ассемблювання.

Отже, на кожному СЕ тепер потрібно обчислити локальну мат­ри­цю жорсткості та вектор на­ван­та­ження . Оскільки ми вибрали лінійні СЕ, тобто елементи з двома вузлами, і у кожному вузлі шу­каєть­ся лише одна невідома величина, то локальна мат­ри­ця жорсткості буде мати розмір (і, відповідно, локальний вектор на­ван­та­ження - розмір 2):

, .

Тоді, враховуючи вище сказане і співвідношення (28), будемо мати, що

,

, (29)

, ,

та

, . (30)

Схематично процес ассемблювання глобальної матриці жорсткості та век­то­ра навантаження з локальних матриць та вектора , відпо­від­но, для тестової сітки з чотирьох СЕ, можна зобразити та­ким чином

, (31)

. (32)

Отже, розв’язання крайової задачі (20)-(21) МСЕ з вико­рис­тан­ням одновимірних кусково-лінійних базисних функцій (25) для зве­лося до СЛАР (27) з матрицею виду (31) та правою частиною виду (32).

Програмна реалізація описаного вище процесу побудови локаль­них матриці жорсткості та вектора навантажень, і їх ассемблювання в середовищі MATHCAD наведена на рис. 3-4. Для обчислення локальних матриць жорсткості на елементі з номером за формулами (29) вве­де­но функцію STIFF(ne). Оскільки базисні функції на першому та ос­тан­ньому СЕ визначаються окремо, то для обчислення локальних матриць жорсткості на цих елементах введено окремі функції STIFF_1 та STIFF_n, від­по­відно. Аналогічним чином організовано процес обчислен­ня локальних векторів навантажень за формулами (30) за допомогою функцій LOAD(ne), LOAD_1 і LOAD_n.

Рис.3. Побудова локальних матриць жорсткості та їх ассемблювання

Для рознесення локальних матриць у глобальні введено функцію FindRow(i,ne), яка дозволяє визначити номер рядка глобальної матриці жорсткості (вектора навантаження) для вузла з номером i на СЕ з номером ne,виду

.

Рис.4. Побудова локальних векторів навантаження та їх ассемблювання

Також, тут використовуються функції

, ,

які призначені для отримання конкретного елемента локальної матриці жорсткості та локального вектора навантаження, відповідно.

ВРАХУВАННЯ ГРАНИЧНИХ УМОВ

Розрізняють два типи граничних умов: головні граничні умови та природні граничні умови. Формальну ознаку поділу граничних умов на головні та природні можна сформулювати таким чином: якщо задано диференціальне рівняння порядку , то граничні умови, що містять по­хідні до порядку включно є головними, а граничні умови, що міс­тять похідні порядку та вище називаються природними граничними умовами.

Рис.5. Врахування граничних умов

Врахування головної граничної умови на лівому(правому) кінці по­лягає в тому, що перший(останній) стовпець глобальної матриці жорст­кості, який є стовпцем коефіцієнтів при відомому значенні, дом­но­жається на це значення і переноситься в праву частину системи, при­чо­му відкидається перша(остання) стрічка глобальної матриці жорсткості і перший (останній) елемент глобального вектора навантаження. Для цього (див. рис. 5) введено функцію модифікації глобального вектора навантаження

та функцію модифікації глобальної матриці жорсткості

.

Природня гранична умова на лівому(правому) кінці для рівняння (20) в загальному випадку має вигляд

, . (33)

Це означає, що, як наслідок інтегрування за частинами, у слабкій формі рівнянь Гальоркіна (24) з’явиться доданок , який з врахуванням умови (33) можна перетворити до вигляду . Тоді, якщо , то відповідне рівняння результуючої СЛАР (тобто перше рівняння, якщо природня умова (33) задана на лівому кінці при і останнє рівняння, якщо природня умова (33) задана на правому кінці при ) доповниться вільним членом , який потрібно просто врахувати у векторі навантаження :

,

.

Якщо , то з’являється додатково доданок , який містить невідому функцію , тобто це означає, що потрібно мо­ди­фіку­вати той елемент глобальної матриці жорсткості , який є коефіцієнтом при відповідному вузловому значення невідомої функції . Таким еле­ментом є перший (останній) елемент на головній діагоналі, якщо при­род­ня умова (33) задана на лівому(правому) кінці. Реалізацію цього випадку задання природної граничної умови в середовищі MATHCAD можна здійснити, визначивши аналогічні функції

,

.

Тоді для врахування однорідної природної граничної умови (21), як видно з рис. 5, залишається здійснити такі виклики наведених вище функцій

, .

Рис.6. Розв’язання СЛАР та оцінка точності наближеного розв’язку