Единое окно доступа к образовательным ресурсам

Метод конечных элементов

Голосов: 4

Метод конечных элементов - один из наиболее эффективных численных методов решения математических задач, описывающих состояние физических систем сложной структуры. На простых примерах пояснены основы метода конечных элементов и проиллюстрированы его основные достоинства.

Приведенный ниже текст получен путем автоматического извлечения из оригинального PDF-документа и предназначен для предварительного просмотра.
Изображения (картинки, формулы, графики) отсутствуют.
                                                                          МАТЕМАТИКА

                                               МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ

                                                                               Л. А. РОЗИН
                                            Санкт-Петербургский государственный технический университет



                                                                                                   В науке и технике постоянно приходится сталки-
                                                                                               ваться с проблемой расчета систем, имеющих сложную
                           THE FINITE ELEMENTS METHOD                                          геометрическую конфигурацию и нерегулярную физи-
                           L. A. ROZIN                                                         ческую структуру. Компьютеры позволяют выполнять
                                                                                               такие расчеты при помощи приближенных численных
                                                                                               методов. Метод конечных элементов (МКЭ) является
                           The finite elements method is one of the most                       одним из них. В последние десятилетия он занял веду-
                           effective numerical methods for solving of                          щее положение и получил широкое применение. В ста-
                           mathematical problems, characterizing the                           тье на простых примерах мы рассмотрим сущность ме-
                           state of physical systems with complicated                          тода конечных элементов и отметим его основные
                                                                                               достоинства.
                           structure. The foundation of the finite ele-
                                                                                                   Предположим, что состояние системы описывает-
                           ments method are illustrated by simple exam-                        ся некоторой функцией. Пусть эта функция является
                           ples, which demonstrate its important advan-                        единственным решением математической задачи, сфор-
                           tages.                                                              мулированной на основе физических законов. Реше-
                                                                                               ние состоит в отыскании из бесконечного множества
                           Метод конечных элементов – один из наи-                             функций такой, которая удовлетворяет уравнениям
                                                                                               задачи. Если задача достаточно сложная, то ее точное
                           более эффективных численных методов ре-
                                                                                               решение невозможно. Вместо того чтобы искать тре-
                           шения математических задач, описываю-                               буемую функцию среди бесконечного множества раз-
                           щих состояние физических систем сложной                             нообразных функций, задача упрощается. Рассматри-
                           структуры. На простых примерах поясне-                              вается некоторое семейство функций, определяемых
                           ны основы метода конечных элементов и                               конечным числом параметров. Как правило, среди та-
                                                                                               ких функций нет точного решения задачи. Однако со-
                           проиллюстрированы его основные досто-                               ответствующим подбором параметров можно попы-
                           инства.                                                             таться приближенно удовлетворить уравнениям задачи
                                                                                               и тем самым построить ее приближенное решение.
                                                                                               Такой общий подход характерен для многих прибли-
                                                                                               женных методов. Специфическим в методе конечных
                                                                                               элементов является построение семейства функций,
                                                                                               определяемых конечным числом параметров.
                                                                                                   Допустим, требуется построить такое семейство
                                                                                               функций u(x) при a x b. Интервал ab разбивается на
                                                                                               конечное число частей (элементов), соединяющихся
                                                                                               между собой и с концами интервала в узловых точках
                                                                                               (узлах) xi (рис. 1). В пределах каждого элемента задается
      © Розин Л.А., 2000




                                                                                               функция, например в виде линейного полинома. Она
                                                                                               определяется своими значениями u(xi) в узлах на кон-
                                                                                               цах элемента. Если отыскиваемая функция является
                                                                                               непрерывной, то значения ее в каждом узле для сосед-
                                                    www.issep.rssi.ru                          них элементов совпадают. В результате имеем семейст-
                                                                                               во кусочно-линейных непрерывных функций, которые



120                                              С О Р О С О В С К И Й О Б РА З О В АТ Е Л Ь Н Ы Й Ж У Р Н А Л , Т О М 6 , № 4 , 2 0 0 0


                                              МАТЕМАТИКА
       u     u1                                                        мента к их полному набору особенно удобен для геоме-
                                                                       трически и физически сложных систем.
                    u2
                          u3                                               2. Каждое отдельное алгебраическое уравнение,
                                  u4
                                                                       полученное на основе метода конечных элементов, со-
                                       u5                              держит незначительную часть узловых неизвестных от
                                                u6                     общего их числа. Другими словами, многие коэффици-
                                                                       енты в уравнениях алгебраической системы равны ну-
             a   1    2    3    4    5 b                               лю, что значительно облегчает ее решение.
              x1   x2   x3   x4   x5   x6 x
                                                                           3. Задачи, решение которых описывается функци-
                         Рис. 1                                        ями, удовлетворяющими функциональным уравнени-
                                                                       ям, носят название континуальных. В отличие от них
                                                                       решение так называемых дискретных задач точно опре-
изображаются в виде ломаных и определяются конеч-
                                                                       деляется конечным числом параметров, удовлетворяю-
ным числом параметров – своими узловыми значения-
                                                                       щих соответствующей системе алгебраических уравне-
ми. На рис. 1 показана одна из функций такого семей-
                                                                       ний. Метод конечных элементов, так же как и другие
ства. Здесь 5 элементов, 6 узлов и 6 узловых параметров
                                                                       численные методы, по существу приближенно заменя-
u(xi) = ui . В случае нескольких переменных схема мето-
                                                                       ет континуальную задачу на дискретную. В методе ко-
да конечных элементов в принципе не меняется. Таким
                                                                       нечных элементов вся процедура такой замены имеет
образом, метод конечных элементов заменяет задачу
                                                                       простой физический смысл. Это позволяет более пол-
отыскания функции на задачу отыскания конечного
                                                                       но представить себе весь процесс решения задачи, из-
числа ее приближенных значений в отдельных точках-
                                                                       бежать многих возможных ошибок и правильно оце-
узлах. При этом если исходная задача относительно
                                                                       нить получаемые результаты.
функции состоит из функционального уравнения, на-
пример дифференциального уравнения с соответству-                          4. Помимо континуальных задач схема метода ко-
ющими граничными условиями, то задача метода ко-                       нечных элементов применяется для соединения эле-
нечных элементов относительно ее значений в узлах                      ментов и формирования алгебраических уравнений
представляет собой систему алгебраических уравнений.                   при решении непосредственно дискретных задач. Это
                                                                       расширяет сферу применения метода.
    С уменьшением максимального размера элементов                          Первая работа, где рассматривалась схема типа
увеличивается число узлов и неизвестных узловых па-                    метода конечных элементов, принадлежит известному
раметров. Вместе с этим повышается возможность бо-                     математику Р. Куранту [2]. Построение метода с ис-
лее точно удовлетворить уравнениям задачи и тем са-                    пользованием физических соображений и его название
мым приблизиться к искомому решению. В настоящее                       “метод конечных элементов” содержатся в статье, на-
время уже изучены многие вопросы, касающиеся схо-                      писанной инженерами [3]. Такое сочетание специаль-
димости приближенного решения методом конечных                         ностей авторов характерно для работ по методу конеч-
элементов к точному. Для линейных задач, когда неиз-                   ных элементов. В последующем было опубликовано
вестные функции и операции над ними входят во все                      много статей и книг, посвященных этому методу и его
соотношения задачи только в первой степени, метод                      различным модификациям. Некоторое представление
конечных элементов получил достаточно полное мате-                     об этом можно получить из списка литературы, напри-
матическое обоснование [1]. В дальнейшем будем рас-                    мер в [4]. Метод конечных элементов реализован в
сматривать только линейные задачи, решение которых                     больших универсальных компьютерных пакетах про-
метод конечных элементов сводит к решению систем ли-                   грамм, которые имеют широкое применение.
нейных алгебраических уравнений. Отметим несколько
важных достоинств метода конечных элементов.                                ПРИМЕР ДИСКРЕТНОЙ ЗАДАЧИ
    1. Метод конечных элементов позволяет построить                    Обратимся к дискретной задаче, состояние которой
удобную схему формирования системы алгебраических                      точно определяется конечным числом параметров.
уравнений относительно узловых значений искомой                        Рассмотрим упругий стержень в виде прямого кругово-
функции. Приближенная аппроксимация решения при                        го цилиндра, длина которого значительно больше его
помощи простых полиномиальных функций и все не-                        диаметра. Это позволяет отождествить стержень с его
обходимые операции выполняются на отдельном типо-                      осью. Пусть три таких стержня расположены на оси x и
вом элементе. Затем производится объединение элемен-                   соединены между собой в точках 1 и 2 (рис. 2, а). Точки
тов, что приводит к требуемой системе алгебраических                   a и b закреплены, что условно изображено на рисунке.
уравнений. Такой алгоритм перехода от отдельного эле-                  К осям стержней вдоль x приложим внешнюю нагрузку.



                                       Р О З И Н Л . А . М Е ТО Д К О Н Е Ч Н Ы Х Э Л Е М Е Н ТО В                               121


                                                                                 МАТЕМАТИКА
                     a                     б                  в                  г                                                                               (r)
                                                                                                                                           (r)       ( r ) du
                                       i                                                                                               N          = c --------- ,
                                                                                                                                                               -                                                   (49)
                 a
                                                                                                                                                            dx
                                                                  N (r)      a
                                                                                     (1)          где c(r) > 0 носит название продольной жесткости
                           (1)
                                                                                                  стержня и определяется из опыта. Пусть c(r) = const для
                 1                             l(r) q(r)(x)               dx 1       1            элемента r. Подставляя (2) в (1) получим задачу относи-
                                 q(r)(x) (r)
                                                                                     (2)
                                                                                                  тельно u(r)(x) в виде дифференциального уравнения и
                                                                                                  граничных условий
                                                                             2       3
                           (2)                                                                                                                         2 (r)
                                       j                                                                                                  (r) d u          (r)
                                                        (r)
                                                      N + dN          (r)
                                                                                     (3)
                                                                                                                                       – c ------------ = q ,
                                                                                                                                                    2
                                                                                                                                                      -
                                                                                                                                               dx
                                           x
                                                                                                                          (r)              (r)              (r)        (r)               (r)
                                                                             3       2                                   u ( 0 ) = ui ,                    u (l ) = u j .                                          (50)
                 2
                                                                                     (4)              В нашем примере для случая дискретной задачи по-
                           (3)                                                                    ложим q(r) = 0 и будем считать, что на стержневую систе-
                                                                             4       4
                                                                                                  му действуют только сосредоточенные силы P1 и P2 со-
                                                                                                  ответственно в узлах 1 и 2. Тогда решение (3) примет вид
                 b                                                                   (5)
                                                                                                                          (r)    (r)
                                                                                                       (r)      u j – ui                   (r)                       (r)             (r)       (r)         (r)
                           x                                                 5       5
                                                                                                      u ( x ) = -------------------- x + u i ,
                                                                                                                         (r)
                                                                                                                                   -                             N           = ж ( u j – u i ),
                                                                                                                       l
                                                                                                                                                                                                                   (51)
                                                                                     (6)                                                                                           (r)
                                                                                                                                                 (r)               (r)     c
                                                                             b
                                                                                                                            0     x         l ,                ж         = ------ .
                                                                                                                                                                             (r)
                                                                                                                                                                                -
                                                                                                                                                                           l

                                                                                     x                На основании (4) можно заключить, что состояние
                                                                                                  типового элемента r, то есть u(r)(x), N(r), точно определя-
                                               Рис. 2                                             ется двумя параметрами – перемещениями его узлов
                                                                                                      (r)    (r)
                                                                                                  u i , u j , что делает задачу дискретной.
      Очевидно, точки на осях стержней перемещаются                                                                                                                                 (r)        (r)
      вдоль x. Силы и перемещения считаются положитель-                                              Рассмотрим внутренние силы f i , f j , действую-
      ными, если они направлены в положительном направ-                                           щие в узлах i, j на элемент r. Поскольку имеет место ли-
      лении x. Задача состоит в определении перемещений                                                                                                                                              (r)     (r)
                                                                                                  нейная задача, то они линейно зависят от u i , u j :
      точек, принадлежащих осям стержней и продольных
                                                                                                                                 (r)             (r) (r)           (r) (r)
      внутренних сил в поперечных сечениях стержней.                                                                            f i = f ii u i + f ij u j ,
                                                                                                                                 (r)             (r) (r)           (r) (r)
                                                                                                                                                                                                                   (52)
          Согласно методу конечных элементов, представим                                                                        f j = f ji u i + f jj u j .
      стержневую систему в виде элементов, соединенных в
                                                                                                                   (r)                                                                                      (r)
      узлах. В качестве элементов примем отдельные стерж-                                         Здесь f lt (l = i, j; t = i, j) есть внутренняя сила f l , дей-
      ни, а узлов – точки 1 и 2. На рис. 2, а в скобках указаны                                   ствующая на элемент r в узле l и возникающая от еди-
      номера элементов. Обратимся к типовому для данной                                           ничного перемещения узла t. При этом перемещение
      системы элементу r. На элемент r с узлами i, j (рис. 2, б )                                                                                                                                                   (r)
                                                                                                  другого узла равно нулю. На рис. 3, а показана сила f ji
      может действовать распределенная нагрузка интенсив-
                                                                                                                                                           (r)                 (r)
                                                 (r) (r)
      ности q(r)(x) и перемещения его узлов u i , u j . Примем                                    для сжатого элемента при u i = 1, u j = 0.
      x = 0 в узле i и обозначим длину элемента l (r). Получим                                       Соотношения (5) можно представить в матричной
      задачу для функции u(r)(x) перемещений точек оси r.                                         форме. Введем столбцы f (r), u(r) и матрицу K(r)
      Бесконечно малая часть элемента dx находится в рав-
                                                                                                              (r)                     (r)                                         (r)                   (r) 
      новесии под действием нагрузки q(r)(x)dx и продольных
                                                                                                            = i          ,    u(r) =  i                 ,                      =  ii
                                                                                                      (r)      f                         u                                   (r)       f                   f ij 
                                                                                                  f                                                                      K                                         .
                                                                                                                                                                                                                   (53)
      внутренних сил N (r)(x), действующих так, как показано                                                  (r)                     (r)                                         (t)                   (r) 
                                                                                                              fj                      uj                                          f ji                f jj 
      на рис. 2, в. Из условия равновесия dx имеем
                                                                                                  Тогда (5) можно записать в виде
                                                              (r)
                     (r)         (r)                  dN
               dN          + q d x = 0,               ----------- = – q ( r ) .
                                                                -                          (48)                                        f (r) = K(r)u(r).                                                           (54)
                                                         dx
                                                                                                     Для упругой пружины коэффициент пропорцио-
          Согласно закону Гука, для упругого стержня                                              нальности между силой и перемещением называется



122                                               С О Р О С О В С К И Й О Б РА З О В АТ Е Л Ь Н Ы Й Ж У Р Н А Л , Т О М 6 , № 4 , 2 0 0 0


                                                                                      МАТЕМАТИКА
                       а                           б                                                       стемы в целом определяется двумя узловыми переме-
                               (r)                                                                         щениями u1 , u2 и рассматриваемая задача является дис-
                   i       u   i     = 1
                                                                                                           кретной.
                                                       (1)
                                                                                                               Для всей системы можно записать соотношения
                                                        (1)
                                                       f 11                                                типа (5) относительно суммарных для смежных эле-
                           (r)
                                               1       u1 = 1                                              ментов внутренних сил в узлах 1, 2. Обозначим их f1 , f2 .
                                                        (2)                                                Очевидно, как и в случае типового элемента, они долж-
                                                       f 11
                   j       u
                               (r)
                               j     = 0                                                                   ны линейно зависеть от u1 , u2 :
                                                       (2)
                                                                                                                                                       f1 = f11u1 + f12u2 ,
                 (r)           (r)
                                               2       u2 = 0
             f ji = f j                                                                                                                                f2 = f21u1 + f22u2 .

                                                       (3)                                                    Введем столбцы f, u и матрицу жесткости всей сис-
                                                                                                           темы K по формулам

                                                                                                                  f                                 u                      f                  f 12   
                                                                                                             f =  1 ,                          u =  1 ,               K =  11                         . (58)
                                                       f 11 = f
                                                                    (1)
                                                                    11    +f
                                                                               (2)
                                                                               11    = f1                         f2                                u2                     f 21               f 22   

                                                                                                           Тогда матричное соотношение типа (7) для всей систе-
                                               Рис. 3
                                                                                                           мы будет
коэффициентом жесткости пружины. Аналогично K(r)                                                                                                             f = Ku.                                         (59)
носит название матрицы жесткости элемента r.                                                               Здесь flt (l = 1, 2; t = 1, 2) есть суммарная для смежных
    От типового элемента перейдем к отдельным эле-                                                         элементов в узле l внутренняя сила fl , возникающая от
ментам данной системы. Для элемента 2 с двумя узлами                                                       единичного перемещения узла t. При этом перемеще-
1, 2 справедливы все зависимости (4)–(7), где следует                                                      ние другого узла равно нулю. Эти суммарные силы оп-
положить r = 2, i = 1, j = 2. Поскольку точки a, b непо-                                                   ределяются через узловые силы в смежных элементах
движны, то состояние элемента 1 определяется пере-
                                                                                                                                     (1)         (2)                     (1)        (2)          (2)
мещением узла 1, а элемента 3 – перемещением узла 2.                                                               f 11 = f 11 + f 11 ,                         f 12 = f 12 + f 12 = f 12 ,
                                                                                                                                                                                                             (60)
На основании (4) будем иметь для элементов 1 и 3 соот-                                                                               (2)
                                                                                                                   f 21 = f 21 + f 21 = f 21 ,
                                                                                                                                                 (3)          (2)                    (2)
                                                                                                                                                                         f 22 = f 22 + f 22 .
                                                                                                                                                                                                 (3)

ношения
                                                                                                                                                                                   (1)     (2)
                             u1
                                     (1)                                                                       На рис. 3, б показаны силы f 11 , f 11 при u1 = 1,
                     (1)                               (1)          (1) (1)
                 u         = ------- x,
                                (1)
                                                   N          =ж u ,    1                                  u2 = 0. При этом элемент 1 растянут, а элемент 2 сжат. В
                              l
                                                                                               (55)                (1)                     (3)
                   u2
                           (3)                                                                             (13) f 12 = 0, f 21 = 0, поскольку узел 2 не принадлежит
           (3)                   (3)                         (3)            (3) (3)
       u         = ------- x + u 2 ,
                      (3)
                                                       N           = –ж u .     2                          элементу 1, а узел 1 не принадлежит элементу 3. Из (13)
                    l
                                                                                                           следует, что матрица жесткости системы строится на
В (8) координата x для элемента 1 равна нулю в точке a,                                                    основе коэффициентов жесткости для отдельных эле-
а для элемента 3 – в узле 2. Зависимости (5) для элемен-                                                   ментов. Алгоритмически выполнить это можно по-раз-
тов 2 и 3 примут соответственно вид                                                                        ному [5]. Например, можно для всех элементов строить
                  (1)                (1) (1)            (3)          (3) (3)                               матрицы жесткости одинаковой размерности равной
             f 1 = f 11 u 1 ,                          f 2 = f 22 u 2 .                        (56)        размерности матрицы K, основываясь на столбце u пе-
Сравнивая (5) или (7) для элемента 2 с (9) для элемен-                                                     ремещений всех узлов системы. Это возможно, по-
                                                                                                                               (r)
тов 1 и 3, можно заключить, что вместо матрицы жест-                                                       скольку f lt = 0, если по крайней мере один из узлов l
кости для двуузлового элемента 2 фигурируют коэффи-                                                        или t не принадлежит элементу r. В данном примере бу-
циенты жесткости для одноузловых элементов 1 и 3.                                                          дем иметь
Теперь все известно о каждом отдельном элементе сис-
темы. Следующим шагом является соединение элемен-                                                                                 (1)                                          (2)              (2) 
                                                                                                                                                                              =  11
                                                                                                                         (1)                           0                (2)       f              f 12 
тов в узлах на основе условий                                                                                      K           =  f 11                    ,       K                                    ,
                                                                                                                                  0                                            (2)              (2) 
            (1)                (2)                      (2)          (3)
                                                                                                                                                       0                         f 21           f 22 
         u1 = u1 = u1 ,                            u2 = u2 = u2 ,                              (57)
                                                                                                                                                       (3)     0         0 
где u1 , u2 – перемещения узлов системы 1, 2. Отсюда                                                                                              K          =                
                                                                                                                                                                           (3)
следует, что состояние соединенных элементов или си-                                                                                                           0        f 22 



                                                                           Р О З И Н Л . А . М Е ТО Д К О Н Е Ч Н Ы Х Э Л Е М Е Н ТО В                                                                               123


                                                                                            МАТЕМАТИКА
      и, согласно (13), K = K(1) + K(2) + K(3). Из условия равно-                                            фициенты fij образуют квадратную матрицу, состоящую
                                                            (r)                 (r)             (r)          из n строк и n столбцов:
      весия элемента r следует f                            lt    = ±N                  при u   l     = 0,
       (r)
      u t = 1, где N (r) > 0 при растяжении и N (r) < 0 при сжа-                                                                f 11 f 12 … f 1n 
                                                                                                                           K =  .............................. .
      тии. В результате на основании (4) и (8) будем иметь                                                                                                    
                                                                                                                                f n1 f n2 … f nn 
                     (1)       (1)                    (2)           (2)               (2)
                   f 11 = ж ,                    f 11 = f 22 = ж ,
                    (2)       (2)               (2)                  (3)               (3)
                                                                                                                   Если обозначить столбец неизвестных u, а столбец
                   f 12 = f 21 = – ж ,                            f 22 = ж ,                                 свободных членов P, то (16) принимает матричную
      а матрица K примет вид                                                                                 форму (15). Система алгебраических уравнений долж-
                                                                                                             на быть невырожденной, то есть иметь единственное
                     (1)     (2)                                   (2)                                     решение. Казалось бы, дальнейшее ясно. Можно вос-
                  K= ж +ж                                   –ж
                                                                                 .                   (61)   пользоваться для решения (16), например, методом ис-
                      –ж
                          (2)
                                                       ж
                                                            (2)
                                                                  +ж
                                                                          (3)                               ключения Гаусса. Однако при применении приближен-
                                                                                                             ных методов обычно приходится иметь дело с системами
          Из условия равновесия узлов 1, 2 следует f1 = P1 ,                                                 большого порядка n, и матрица, вообще говоря, может
      f2 = P2 или для столбцов f = P, где P – столбец из P1 , P2 .                                           иметь такую структуру, которая затрудняет получение
      Подставляя сюда вместо f его выражение согласно (12),                                                  решения. При этом на точности результата в той или
      окончательно получим систему алгебраических урав-                                                      иной степени сказываются неизбежные в процессе вы-
      нений относительно u1 , u2                                                                             числений ошибки округления. Одним из важных до-
                                          (1)           (2)                     (2)                          стоинств метода конечных элементов является то, что
                                     (ж         + ж )u 1 – ж u 2 = P1 ,                                      он обычно приводит к таким системам алгебраических
         Ku = P     или               (2)                     (2)           (3)
                                                                                                      (62)
                                    –ж u1 + ( ж                     + ж )u 2 = P2 .                          уравнений, матрицы K которых позволяют эффектив-
                                                                                                             но строить решение.
      После определения u1 , u2 в результате решения (15) на-                                                      Выясним, какой желательно иметь матрицу K в
      ходятся u(r)(x), N (r) во всех элементах системы при по-                                               (16). Пределом мечты была бы система (16) с диаго-
      мощи (4) и (8).                                                                                        нальной матрицей K, когда все fij = 0 при i j. В этом
          Таким образом, схема метода конечных элементов                                                     случае (16) распадается на отдельные уравнения
      для дискретных задач состоит из представления систе-                                                   fiiui = Pi . Такое может быть, только если в физической
      мы в виде совокупности отдельных элементов, исполь-                                                    системе, рассчитываемой методом конечных элемен-
      зования точного решения для типового элемента и со-                                                    тов, узлы между собой не связаны, то есть по существу
      единения элементов в систему. Матрица жесткости                                                        системы не существует. Однако теперь уже ясно, к чему
      всей системы определяется посредством матриц жест-                                                     надо стремиться: следует так выполнять процесс пост-
      кости отдельных элементов и является матрицей систе-                                                   роения алгебраической системы уравнений, чтобы мат-
      мы алгебраических уравнений относительно неизвест-                                                     рица по возможности содержала больше нулевых ко-
      ных узловых перемещений. Наличие точного решения                                                       эффициентов и была близка к диагональной, другими
      для типового элемента, зависящего от конечного числа                                                   словами, желательно, чтобы в каждое уравнение входи-
      параметров – узловых перемещений, делает задачу                                                        ло относительно небольшое число неизвестных в со-
      дискретной.                                                                                            седних узлах.
                                                                                                                   Матрицы, близкие к диагональным, обычно имеют
             СИСТЕМА ЛИНЕЙНЫХ                                                                                ленточную структуру, когда все ненулевые и некоторые
             АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ                                                                        нулевые коэффициенты находятся между двумя лини-
      Метод конечных элементов сводит решение линейной                                                       ями, параллельными главной диагонали. Например,
      задачи к решению системы линейных алгебраических                                                                                        t1
      уравнений
                                                                                                                                   *     0        *   0    0   
                       f11u1 + f12u2 + … + f1nun = P1 ,                                                                                                        
                                                                                                                          t        0     *        *   *    0   
                                                                                                                       K= 2                                    ,   (64)
                   .............................................                                      (63)                         *     *        *   0    0   
                                                                                                                                   0     *        0   *    *   
                      fn1u1 + fn2u2 + … + fnnun = Pn .                                                                                                         
                                                                                                                                   0     0        0   *    *   
      Здесь ui (i = 1, 2, …, n) – неизвестные, Pi (i = 1, 2, …, n) –
      заданные свободные члены, fij (i, j = 1, 2, …, n) – коэф-                                              где знак * заменяет коэффициенты, отличные от нуля,
      фициенты при неизвестных. По аналогии с (11) коэф-                                                     а главная диагональ и параллельные ей линии указаны



124                                               С О Р О С О В С К И Й О Б РА З О В АТ Е Л Ь Н Ы Й Ж У Р Н А Л , Т О М 6 , № 4 , 2 0 0 0


                                               МАТЕМАТИКА
пунктиром. Ленточную матрицу характеризует ширина                                                      (1)             (2)
                                                                                                   ж +ж                            P1
ленты t = t1 + t2 + 1, равная наибольшему числу коэф-                                        u 2 = ----------------------- u 1 – --------,
                                                                                                               (2)
                                                                                                                         -           (2)
                                                                                                           ж                     ж
фициентов в строке в пределах ленты. В данном случае                                                     (2)
                                                                                                                                                                               (65)
t1 = t2 = 2 и t = 5. Для диагональной матрицы t = 1. При                                             ж                               P2
                                                                                       u 2 = -----------------------) u 1 + -----------------------) .
                                                                                                 (2)             (3
                                                                                                                   -            (2)             (3
                                                                                                                                                  -
решении системы уравнений с ленточной матрицей                                               ж +ж                           ж +ж
участвуют только те коэффициенты, которые располо-
жены в пределах ленты. Число арифметических опера-                        В прямоугольной системе координат u1 , u2 на рис. 4
ций, необходимых для решения системы алгебраичес-                     уравнение прямой будет u2 = u1 tgα + g, где α – угол
ких уравнений с полностью заполненной матрицей                        между прямой и положительным направлением оси u1 ,
методом Гаусса, при больших n имеет порядок n3. В то                  g – отрезок отсекаемый прямой на оси u2 . Уравнения
же время для ленточной матрицы при t1 = t2 и t1 n он                  (18) описывают две прямые на рис. 4, а решение (18)
составляет nt 1 .
                 2                                                    представляет собой координаты точки пересечения
                                                                      этих прямых. Здесь
    Для примера ленточной матрицы обратимся к за-
                                                                                                   (1)               (2)                                        (2)
дачам предыдущего раздела, но с пятью узлами и шес-                                       ж +ж                                                            ж
                                                                                 tg α 1 = ----------------------- ,
                                                                                                      (2)
                                                                                                                -                        tg α 2 = -----------------------) .
                                                                                                                                                      (2)             (3
                                                                                                                                                                        -
тью элементами на рис. 2, г. Аналогично (11) матрица K                                            ж                                               ж +ж
будет иметь коэффициенты flt . По смыслу flt они отлич-
ны от нуля только для тех узлов l, где перемещение узла                   Если α1 = α2 и прямые параллельны, то решение
t вызывает отличную от нуля силу при условии, что ос-                 системы (18) не существует и она является вырожден-
тальные узлы, кроме t, неподвижны. Отсюда при нуме-                   ной. Если α1 и α2 различаются мало, то система близка
рации узлов, показанной на рис. 2, г, слева от оси x имеем            к вырожденной. При этом незначительные изменения
                                                                      углов α1 , α2 сильно скажутся на координатах точки пе-
                                                                    ресечения прямых, то есть на решении. Таким образом,
                    *    *   0   0   0                              плохая обусловленность объясняется тем, что система
                    *    *   *   0   0   
                                                                      является почти вырожденной.
               K=                        .
                    0    *   *   *   0   
                    0    0   *   *   *                                  В качестве примера обратимся к (18). Пусть
                                                                    ж(1) = ж(3) = ж и ж(2)   ж, то есть элемент 2 на рис. 2, а
                    0    0   0   *   *   
                                                                      значительно более жесткий, чем элементы 1 и 3. При
                                                                      этом tgα1 ≈ tgα2 и система (18) почти вырожденная. В
Здесь t = 3 и матрица K является трехдиагональной.
                                                                      данном случае разумно изменить постановку задачи и
    При применении метода конечных элементов ши-                      считать элемент 2 абсолютно жестким по сравнению с
рина полосы ленточной матрицы зависит от нумера-                      элементами 1 и 3. Это позволяет объединить узлы 1 и 2
ции узлов. Например, если пронумеровать узлы так,                     в один узел, который обозначим 12, и приложить к не-
как показано на рис. 2, г справа от оси x, то K примет                му суммарную силу P12 = P1 + P2 . Если в (18) положить
вид (17). Вообще если элементы имеют несколько уз-                    u1 = u2 = u12 , вычесть из первого уравнения второе и по-
лов, то при t1 = t2 величина t1 равна максимальной по                 сле преобразований пренебречь ж по сравнению с ж(2),
элементам величине наибольшей разности между но-                      то задача сведется к одному уравнению 2жu12 = P12 .
мерами узлов в отдельном элементе. В первом случае
нумерации узлов слева на рис. 2, г    t1 = 1, а при нуме-                                                       u2
рации справа t1 = 2.
    В некоторых случаях исходная постановка задачи
может оказаться настолько плохой, что даже метод ко-                                                                            P2
                                                                                                                       ------------------------)
                                                                                                                           (1)             (3
нечных элементов не может помочь. И надо ее менять.                                                                    ж +ж
При этом имеет место система алгебраических уравне-                                              α2                                                α1
ний, в которой малые изменения коэффициентов или                                                                       0                                               u1
                                                                                                           P1
свободных членов приводят к значительному измене-                                                        --------)
                                                                                                             (2
                                                                                                                -
                                                                                                         ж
нию решения. Такие системы уравнений носят назва-
ние плохо обусловленных. Выясним, в чем причина
плохой обусловленности на примере системы (15), ко-
торую перепишем в виде                                                                                                     Рис. 4




                                      Р О З И Н Л . А . М Е ТО Д К О Н Е Ч Н Ы Х Э Л Е М Е Н ТО В                                                                                     125


                                                                     МАТЕМАТИКА
          ПРИМЕР КОНТИНУАЛЬНОЙ ЗАДАЧИ                                                          узлах и определяемых своими узловыми значениями ui ,
                                                                                               i = 1, 2, …, n. На концах интервала 0, 1 они обращаются
      Обратимся к задаче (3) для одного элемента. В общем
                                                                                               в нуль. Каждую из таких функций можно изобразить в
      случае задания q(r)(x) она является континуальной зада-                                  виде ломаной линии.
                                                                       (r)        (r)
      чей. Для простоты положим с(r) = 1, l (r) = 1, u i = u j = 0
                                                                                                   Остается определить ui в (20). Это можно сделать
      и опустим индекс r, тогда задача (3) будет                                               по-разному путем приближенного удовлетворения
                  −u" = q(x),              u(0) = u(1) = 0,                             (66)   уравнению в (19). Однако, поскольку уравнение в (19)
                                                                                               содержит u", а уже u' в (20) терпит разрывы непрерыв-
      где штрихи означают дифференцирование по x. Со-                                          ности в узлах, воспользуемся следующим приемом.
      гласно схеме метода конечных элементов разобьем ин-                                      Обозначим R(x) = u"(x) + q(x) невязку уравнения в (19).
      тервал 0, 1 на элементы, соединенные в узлах xi (i = 0, 1,                               Точное решение дает R(x) = 0, и, следовательно,
      …, n + 1) (рис. 5). Будем разыскивать приближенное ре-
                                                                                                                     1
      шение (19) среди функций семейства с конечным чис-
      лом параметров в виде                                                                                          ∫ [ u" ( x ) + q ( x ) ]ϕ ( x ) dx = 0                                  (68)
      u(x) = u0ϕ0(x) + u1ϕ1(x) + … + unϕn(x) + un + 1ϕn + 1(x). (67)
                                                                                                                     0


                                                                                               для любых функций ϕ(x), которые носят название тес-
      Здесь u(x) приближенно представлена линейной ком-
                                                                                               товых. Поскольку разыскивается приближенное реше-
      бинацией некоторых функций ϕi(x) с коэффициентами
                                                                                               ние в форме (20) и для него, как правило, R(x) 0, то
      (параметрами) ui = u(xi) – неизвестными значениями
                                                                                               выполнение тестового условия (21) на базе (20) невоз-
      искомой функции в узлах xi . Для того чтобы в (20)
                                                                                               можно. Смягчим выполнение условия (21), потребо-
      u(xi) = ui во всех узлах xi , функции ϕi(x) должны удов-
                                                                                               вав, чтобы оно выполнялось только для n функций
      летворять условиям ϕi(xi) = 1 и ϕi(xj) = 0 для всех узлов xj
                                                                                               ϕj(x), которые совпадают с пробными. Такой прием но-
      при j i. Кроме того, чтобы выполнялись граничные
                                                                                               сит название метода Галёркина. Выполним в (21) инте-
      условия (19), следует в (20) положить u0 = un + 1 = 0. В ос-
                                                                                               грирование по частям при условии ϕ(x) = ϕj(x) и ϕj(0) =
      тальном функции ϕi(x), которые носят название проб-
                                                                                               = ϕj(1) = 0, тогда вместо (21) получим
      ных, можно выбирать в довольно широких пределах.
      Общие требования к ним состоят в возможности вы-                                                 1
      полнить процесс построения приближенного решения
      и на основе (20) при n        ∞ осуществить сколь угодно                                         ∫ ( – u'ϕ' + qϕ ) dx = 0,
                                                                                                                     j             j                       j = 1, 2, …, n.                   (69)
                                                                                                       0
      точно соответствующую аппроксимацию любой функ-
      ции, среди которых разыскивается решение задачи.                                             Теперь уже в задачу (22) входит u' и можно подста-
      Очевидно, выбор ϕi(x) играет важнейшую роль как в от-                                    вить u из (20) в (22), что дает систему линейных алгеб-
      ношении трудоемкости расчета, так и точности резуль-                                     раических уравнений относительно ui вида (16) с коэф-
      тата. Метод конечных элементов оперирует в качестве                                      фициентами fij и свободными членами Pi
      ϕi(x) кусочно-полиномиальными функциями, отлич-
                                                                                                                           1                                1
      ными от нуля в пределах небольшого числа элементов
      вблизи узла xi . Именно это делает метод максимально                                                               i ∫
                                                                                                                f ij = ϕ 'ϕ 'j d x,                         ∫
                                                                                                                                                      P i = qϕ j d x.                        (70)
      эффективным. Поскольку u(x) по своему физическому                                                                    0                                0

      смыслу должна быть непрерывной функцией, выберем
                                                                                               Здесь fij = fji и матрица K симметричная, что характерно
      ϕi(x) в виде кусочно-линейных функций-“домиков”,
                                                                                               для метода Галёркина. Для простоты примем длину
      отличных от нуля на двух элементах (см. рис. 5). Каж-
                                                                                               элементов одинаковой и равной h. Согласно рис. 5, на-
      дая такая функция ϕi(x), i = 1, 2, …, n, равна единице в
                                                                                               клон ϕ i' функции ϕi равен 1/ h на интервале xi − 1 , xi и
      xi и нулю во всех остальных узлах. При этом набор
                                                                                               −1/h на интервале xi , xi + 1 .
      функций u(x) в (20) будет состоять из непрерывных
      функций линейных в пределах элементов с изломами в                                           Кроме того, произведение ϕ i'ϕ 'j отлично от нуля
                                                                                               только при j = i, j = i 1, когда соответствующие два
                             ϕi − 1        ϕi         ϕi + 1                                   элемента, которые несут на себе функции ϕi и ϕj , пере-
                                                                                               крываются (см. рис. 5). В противном случае ϕ 'ϕ 'j = 0.
                                                                                                                                             i
                                            1                                                  Если i = j, то
                         h            h           h            h                       x
                                                                                                                xi + 1                     xi                    xi + 1
      x0 = 0 x1              xi − 1        xi         xi + 1                 xn xn + 1 = 1                                                             2
                                                                                                                                                    -- d x +
                                                                                                                                                     1                        1 2
                                                                                                                                                                           – -- d x = 2
                                                                                                                 ∫       ( ϕ 'i ) d x =    ∫                      ∫
                                                                                                                               2
                                                                                                       f ij =                                         -                        -        --
                                                                                                                                                                                         -
                                                                                                                                                    h                    h          h
                                          Рис. 5                                                                xi – 1                    xi – 1                  xi




126                                             С О Р О С О В С К И Й О Б РА З О В АТ Е Л Ь Н Ы Й Ж У Р Н А Л , Т О М 6 , № 4 , 2 0 0 0


                                            МАТЕМАТИКА
Аналогично fi j = −1/ h при j = i 1. Следовательно,                  ментов для континуальной и дискретной задач в основ-
матрица K в данном случае оказывается трехдиаго-                     ном совпадают. Здесь, так же как и в случае дискретной
нальной:                                                             задачи, можно выполнить построение матрицы K(r) для
                                                                     типового элемента и из них в процессе соединения
                2    –1                                            элементов в систему сформировать матрицу K для всей
                                                                   системы. Аналогично формируются и свободные чле-
                –1   2    –1         0 
        K = --                          .
            1                                                        ны уравнений. Алгоритм метода конечных элементов
             -                                           (71)
            h                                                      особенно эффективен для решения двух- и трехмер-
                0         –1   2     –1                            ных задач, где проявляются основные преимущества
                                        
                               –1    2                             этого метода.

Согласно (23), интегрирование в Pi совершается только                     ЛИТЕРАТУРА
на двух соседних элементах. Решение полученной сис-                  1. Стренг Г., Фикс Дж. Теория метода конечных элементов.
темы алгебраических уравнений дает ui и позволяет                    М.: Мир, 1977. 349 с.
представить приближенное решение в форме (20).                       2. Courant R. // Bull. Amer. Math. Soc. 1943. Vol. 49. P. 1–43.
    В данном примере непрерывность пробных функ-                     3. Turner M., Clough R., Martin H., Topp L. // J. Aeronaut Sci.
ций ϕi позволила воспользоваться (22). Кроме того, все               1956. Vol. 23, № 9. P. 805–823.
функции ϕi отличны от нуля на разных интервалах 2h,                  4. Зинкевич О., Морган К. Конечные элементы и аппроксима-
что делает их существенно различными и построенную                   ция. М.: Мир, 1986. 318 с.
на их основе при помощи (23) систему линейных ал-                    5. Розин Л.А. Стержневые системы как системы конечных
                                                                     элементов. Л.: Изд-во ЛГУ, 1976. 232 с.
гебраических уравнений невырожденной. Более того,
матрица (24) оказалась ленточной и каждое уравнение
                                                                                  Рецензент статьи Ю.Г. Мартыненко
связывает не более трех неизвестных в соседних узлах.
Полученное приближенное решение (20) в виде лома-
                                                                                                   ***
ной линии хорошо аппроксимирует решение задачи
при достаточно больших n.                                            Леонид Александрович Розин, доктор физико-матема-
                                                                     тических наук, профессор, зав. кафедрой строитель-
    Таким образом, для континуальной задачи метод
                                                                     ной механики и теории упругости Санкт-Петербург-
конечных элементов осуществляет приближенный пе-                     ского государственного технического университета,
реход к дискретной задаче на основе (20) и соответству-              заслуженный деятель науки и техники РФ. Область
ющих кусочно-полиномиальных функций ϕi , отличных                    научных интересов – численные методы решения за-
от нуля на нескольких соседних элементах, содержащих                 дач механики деформируемых систем. Автор более
узел xi . Дальнейшие процедуры метода конечных эле-                  160 статей и семи монографий.




                                     Р О З И Н Л . А . М Е ТО Д К О Н Е Ч Н Ы Х Э Л Е М Е Н ТО В                                       127



    
Яндекс цитирования Яндекс.Метрика