J. Пример применения функционала — Метод конечных элементов

Рассмотрим, как функционал метода Бубнова-Галёркина из приложения I применяется на практике: разберём одномерный пример от начала до конца — от пробных функций до системы линейных уравнений в матричной форме.

Решаем уравнение теплопроводности $L[u] = f$ , где параболический оператор $L$ задан в (I.2), на отрезке $[x_0, x_4]$ . На концах отрезка задаются краевые условия первого рода (Дирихле, (1.13)) или второго рода (Нейман, (1.14)). В этом приложении мы доводим до матричной формы сам функционал Бубнова-Галёркина, сохраняя граничный интеграл $\int_S \upsilon \cdot \nabla \upsilon \,dS$ в общем виде; подстановка конкретных граничных значений (условий Дирихле или Неймана) разбирается в отдельных приложениях.

Одномерная сетка из четырёх симплексов-отрезков s₁…s₄ и пяти узлов — **Рис. J.1.** Одномерная сетка примера: четыре симплекса-отрезка $s_1, \ldots, s_4$ и пять узлов $x_0, \ldots, x_4$ ; длины отрезков обозначаются $l_{(i)(i+1)}$ .

Возьмём четыре симплекса-отрезка $s_1, s_2, s_3, s_4$ , три внутренние точки $x_1, x_2, x_3$ и две граничные точки $x_0, x_4$ (рисунок J.1). Пробные функции для симплексов-отрезков примут вид

\begin{split} &\upsilon_{(0)(1)}(x) = q_0 \cdot \phi_0(x) + q_1 \cdot \phi_1(x) \\ &\upsilon_{(1)(2)}(x) = q_1 \cdot \phi_1(x) + q_2 \cdot \phi_2(x) \\ &\upsilon_{(2)(3)}(x) = q_2 \cdot \phi_2(x) + q_3 \cdot \phi_3(x) \\ &\upsilon_{(3)(4)}(x) = q_3 \cdot \phi_3(x) + q_4 \cdot \phi_4(x) \\ \end{split}

Функции «крышки» и их производные для симплексов-отрезков имеют вид

\begin{split} &\phi_0(x) = a_0 + b_0 \cdot x \quad \phi_0'(x) = b_0 \\ &\phi_1(x) = a_1 + b_1 \cdot x \quad \phi_1'(x) = b_1 \\ &\phi_2(x) = a_2 + b_2 \cdot x \quad \phi_2'(x) = b_2 \\ &\phi_3(x) = a_3 + b_3 \cdot x \quad \phi_3'(x) = b_3 \\ &\phi_4(x) = a_4 + b_4 \cdot x \quad \phi_4'(x) = b_4 \\ \end{split}

Коэффициенты принимают следующие значения. Введём обозначение $l_{(i)(i+1)} = x_{i+1} - x_i$ для длины симплекса между точками $x_i$ и $x_{i+1}$ .

Коэффициенты «крышек» для произвольного симплекса между точками $x_i$ и $x_{i+1}$ выводятся так же, как в разделе «Пробные функции», и в обозначении $l_{(i)(i+1)} = x_{i+1} - x_i$ равны

\begin{split} &a_i = \frac{\displaystyle x_{i+1}}{\displaystyle l_{(i)(i+1)}} \quad b_i = \frac{\displaystyle -1}{\displaystyle l_{(i)(i+1)}} \\ &a_{i+1} = \frac{\displaystyle -x_i}{\displaystyle l_{(i)(i+1)}} \quad b_{i+1} = \frac{\displaystyle 1}{\displaystyle l_{(i)(i+1)}} \end{split}

Для краевых «крышек» $\phi_0$ и $\phi_4$ берётся соответствующая половина. Подставив эти коэффициенты, сразу запишем пробные функции на каждом симплексе.

Подставим коэффициенты в пробные функции и выразим их относительно $x$ :

Для симплекса $s_1$ :

\begin{split} \upsilon_{(0)(1)}(x) &= q_0 \cdot \phi_0(x) + q_1 \cdot \phi_1(x) \\ &= q_0 \cdot \left(\frac{\displaystyle x_1}{\displaystyle l_{(0)(1)}} - \frac{\displaystyle x}{\displaystyle l_{(0)(1)}}\right) + q_1 \cdot \left(\frac{\displaystyle -x_0}{\displaystyle l_{(0)(1)}} + \frac{\displaystyle x}{\displaystyle l_{(0)(1)}}\right) \\ &= \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} \left[q_0 \cdot (x_1 - x) + q_1 \cdot (x - x_0)\right] \end{split}

Для симплекса $s_2$ :

\begin{split} \upsilon_{(1)(2)}(x) &= q_1 \cdot \phi_1(x) + q_2 \cdot \phi_2(x) \\ &= q_1 \cdot \left(\frac{\displaystyle x_2}{\displaystyle l_{(1)(2)}} - \frac{\displaystyle x}{\displaystyle l_{(1)(2)}}\right) + q_2 \cdot \left(\frac{\displaystyle -x_1}{\displaystyle l_{(1)(2)}} + \frac{\displaystyle x}{\displaystyle l_{(1)(2)}}\right) \\ &= \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right] \end{split}

Для симплекса $s_3$ :

\begin{split} \upsilon_{(2)(3)}(x) &= q_2 \cdot \phi_2(x) + q_3 \cdot \phi_3(x) \\ &= q_2 \cdot \left(\frac{\displaystyle x_3}{\displaystyle l_{(2)(3)}} - \frac{\displaystyle x}{\displaystyle l_{(2)(3)}}\right) + q_3 \cdot \left(\frac{\displaystyle -x_2}{\displaystyle l_{(2)(3)}} + \frac{\displaystyle x}{\displaystyle l_{(2)(3)}}\right) \\ &= \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right] \end{split}

Для симплекса $s_4$ :

\begin{split} \upsilon_{(3)(4)}(x) &= q_3 \cdot \phi_3(x) + q_4 \cdot \phi_4(x) \\ &= q_3 \cdot \left(\frac{\displaystyle x_4}{\displaystyle l_{(3)(4)}} - \frac{\displaystyle x}{\displaystyle l_{(3)(4)}}\right) + q_4 \cdot \left(\frac{\displaystyle -x_3}{\displaystyle l_{(3)(4)}} + \frac{\displaystyle x}{\displaystyle l_{(3)(4)}}\right) \\ &= \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} \left[q_3 \cdot (x_4 - x) + q_4 \cdot (x - x_3)\right] \end{split}

Найдём производные пробных функций по $x$ :

Для симплекса $s_1$ :

\upsilon'_{(0)(1)}(x) = \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} (q_1 - q_0)

Для симплекса $s_2$ :

\upsilon'_{(1)(2)}(x) = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} (q_2 - q_1)

Для симплекса $s_3$ :

\upsilon'_{(2)(3)}(x) = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} (q_3 - q_2)

Для симплекса $s_4$ :

\upsilon'_{(3)(4)}(x) = \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} (q_4 - q_3)

Вычислим квадраты пробных функций:

Для симплекса $s_1$ :

\begin{split} \upsilon_{(0)(1)}^2(x) &= \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} \left[q_0 \cdot (x_1 - x) + q_1 \cdot (x - x_0)\right]^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} \left[q_0^2(x_1 - x)^2 + 2q_0q_1(x_1 - x)(x - x_0) + q_1^2(x - x_0)^2\right] \end{split}

Для симплекса $s_2$ :

\begin{split} \upsilon_{(1)(2)}^2(x) &= \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right]^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left[q_1^2(x_2 - x)^2 + 2q_1q_2(x_2 - x)(x - x_1) + q_2^2(x - x_1)^2\right] \end{split}

Для симплекса $s_3$ :

\begin{split} \upsilon_{(2)(3)}^2(x) &= \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right]^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \left[q_2^2(x_3 - x)^2 + 2q_2q_3(x_3 - x)(x - x_2) + q_3^2(x - x_2)^2\right] \end{split}

Для симплекса $s_4$ :

\begin{split} \upsilon_{(3)(4)}^2(x) &= \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} \left[q_3 \cdot (x_4 - x) + q_4 \cdot (x - x_3)\right]^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} \left[q_3^2(x_4 - x)^2 + 2q_3q_4(x_4 - x)(x - x_3) + q_4^2(x - x_3)^2\right] \end{split}

Вычислим квадраты производных пробных функций:

Для симплекса $s_1$ :

\begin{aligned} \left(\upsilon'_{(0)(1)}(x)\right)^2 &= \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} (q_1 - q_0)^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}^2} \left(q_0^2 - 2q_0q_1 + q_1^2\right) \end{aligned}

Для симплекса $s_2$ :

\begin{aligned} \left(\upsilon'_{(1)(2)}(x)\right)^2 &= \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} (q_2 - q_1)^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left(q_1^2 - 2q_1q_2 + q_2^2\right) \end{aligned}

Для симплекса $s_3$ :

\begin{aligned} \left(\upsilon'_{(2)(3)}(x)\right)^2 &= \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} (q_3 - q_2)^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \left(q_2^2 - 2q_2q_3 + q_3^2\right) \end{aligned}

Для симплекса $s_4$ :

\begin{aligned} \left(\upsilon'_{(3)(4)}(x)\right)^2 &= \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} (q_4 - q_3)^2 \\ &= \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}^2} \left(q_3^2 - 2q_3q_4 + q_4^2\right) \end{aligned}

Запишем линейную аппроксимацию функции источника $f(x)$ на каждом симплексе согласно (6.36). В нашем примере источник действует только на внутренних симплексах: на крайних он тождественно равен нулю.

Для симплекса $s_1$ (между $x_0$ и $x_1$ ):

Поскольку функция плотности тепловых источников тождественно равна нулю на этом интервале:

\widetilde{f}_{(0)(1)}(x) = 0.

Для симплекса $s_2$ (между $x_1$ и $x_2$ ):

\widetilde{f}_{(1)(2)}(x) = e_1^{(1)(2)} + e_2^{(1)(2)} \cdot x,

где:

\begin{split} &e_1^{(1)(2)} = \frac{\displaystyle f_1 \cdot x_2 - f_2 \cdot x_1}{\displaystyle x_2 - x_1} = \frac{\displaystyle f_1 \cdot x_2 - f_2 \cdot x_1}{\displaystyle l_{(1)(2)}} \\ &e_2^{(1)(2)} = \frac{\displaystyle f_2 - f_1}{\displaystyle x_2 - x_1} = \frac{\displaystyle f_2 - f_1}{\displaystyle l_{(1)(2)}} \end{split}

Для симплекса $s_3$ (между $x_2$ и $x_3$ ):

\widetilde{f}_{(2)(3)}(x) = e_1^{(2)(3)} + e_2^{(2)(3)} \cdot x,

где:

\begin{split} &e_1^{(2)(3)} = \frac{\displaystyle f_2 \cdot x_3 - f_3 \cdot x_2}{\displaystyle x_3 - x_2} = \frac{\displaystyle f_2 \cdot x_3 - f_3 \cdot x_2}{\displaystyle l_{(2)(3)}} \\ &e_2^{(2)(3)} = \frac{\displaystyle f_3 - f_2}{\displaystyle x_3 - x_2} = \frac{\displaystyle f_3 - f_2}{\displaystyle l_{(2)(3)}} \end{split}

Для симплекса $s_4$ (между $x_3$ и $x_4$ ):

Поскольку функция плотности тепловых источников тождественно равна нулю на этом интервале:

\widetilde{f}_{(3)(4)}(x) = 0.

Теперь запишем произведение функции источника на соответствующие пробные функции для каждого симплекса.

Для симплекса $s_1$ :

Поскольку $\widetilde{f}_{(0)(1)}(x) = 0$ , имеем:

\widetilde{f}_{(0)(1)}(x) \cdot \upsilon_{(0)(1)}(x) = 0.

Для симплекса $s_2$ :

\begin{aligned} \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) &= \left(e_1^{(1)(2)} + e_2^{(1)(2)} \cdot x\right) \\ &\quad \cdot \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right]. \end{aligned}

Подставляя значения $e_1^{(1)(2)}$ и $e_2^{(1)(2)}$ :

\begin{aligned} \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) &= \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left[(f_1 \cdot x_2 - f_2 \cdot x_1) + (f_2 - f_1) \cdot x\right] \\ &\quad \cdot \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right]. \end{aligned}

Преобразуем первый множитель:

(f_1 \cdot x_2 - f_2 \cdot x_1) + (f_2 - f_1) \cdot x = f_1 \cdot (x_2 - x) + f_2 \cdot (x - x_1).

Следовательно:

\begin{aligned} \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) &= \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \left[f_1 \cdot (x_2 - x) + f_2 \cdot (x - x_1)\right] \\ &\quad \cdot \left[q_1 \cdot (x_2 - x) + q_2 \cdot (x - x_1)\right]. \end{aligned}

Раскроем скобки:

\begin{aligned} \widetilde{f}_{(1)(2)}(x) \cdot \upsilon_{(1)(2)}(x) = \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}^2} \Big[ &f_1 q_1 (x_2 - x)^2 + f_1 q_2 (x_2 - x)(x - x_1) \\ &+ f_2 q_1 (x - x_1)(x_2 - x) + f_2 q_2 (x - x_1)^2 \Big]. \end{aligned}

Для симплекса $s_3$ :

\begin{aligned} \widetilde{f}_{(2)(3)}(x) \cdot \upsilon_{(2)(3)}(x) &= \left(e_1^{(2)(3)} + e_2^{(2)(3)} \cdot x\right) \\ &\quad \cdot \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right]. \end{aligned}

По аналогии с симплексом $s_2$ :

\begin{aligned} \widetilde{f}_{(2)(3)}(x) \cdot \upsilon_{(2)(3)}(x) &= \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \left[f_2 \cdot (x_3 - x) + f_3 \cdot (x - x_2)\right] \\ &\quad \cdot \left[q_2 \cdot (x_3 - x) + q_3 \cdot (x - x_2)\right]. \end{aligned}

Раскроем скобки:

\begin{aligned} \widetilde{f}_{(2)(3)}(x) \cdot \upsilon_{(2)(3)}(x) = \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}^2} \Big[ &f_2 q_2 (x_3 - x)^2 + f_2 q_3 (x_3 - x)(x - x_2) \\ &+ f_3 q_2 (x - x_2)(x_3 - x) + f_3 q_3 (x - x_2)^2 \Big]. \end{aligned}

Для симплекса $s_4$ :

Поскольку $\widetilde{f}_{(3)(4)}(x) = 0$ , имеем:

\widetilde{f}_{(3)(4)}(x) \cdot \upsilon_{(3)(4)}(x) = 0.

Для левой границы $x_0$ (смежный симплекс $s_1 = [x_0, x_1]$ ) производная по внешней нормали равна

\begin{aligned} \frac{\displaystyle \partial \upsilon(x)}{\displaystyle \partial n} \bigg|_{x=x_0} &= \frac{\displaystyle q_0 - q_1}{\displaystyle x_1 - x_0} \\ &= \frac{\displaystyle q_0 - q_1}{\displaystyle l_{(0)(1)}}. \end{aligned}

Для правой границы $x_4$ (смежный симплекс $s_4 = [x_3, x_4]$ ) производная по внешней нормали равна

\begin{aligned} \frac{\displaystyle \partial \upsilon(x)}{\displaystyle \partial n} \bigg|_{x=x_4} &= \frac{\displaystyle q_4 - q_3}{\displaystyle x_4 - x_3} \\ &= \frac{\displaystyle q_4 - q_3}{\displaystyle l_{(3)(4)}}. \end{aligned}

Интеграл по границе согласно (6.51) принимает вид

\begin{aligned} \int_S \upsilon \cdot \nabla \upsilon \,dS &= q_4 \cdot \frac{\displaystyle \partial \upsilon(x)}{\displaystyle \partial n} \bigg|_{x=x_4} \\ &\quad + q_0 \cdot \frac{\displaystyle \partial \upsilon(x)}{\displaystyle \partial n} \bigg|_{x=x_0}. \end{aligned}

Подставляя найденные производные:

\begin{aligned} \int_S \upsilon \cdot \nabla \upsilon \,dS &= q_4 \cdot \frac{\displaystyle q_4 - q_3}{\displaystyle l_{(3)(4)}} \\ &\quad + q_0 \cdot \frac{\displaystyle q_0 - q_1}{\displaystyle l_{(0)(1)}}. \end{aligned}

Выразим относительно $q$ :

\begin{aligned} \int_S \upsilon \cdot \nabla \upsilon \,dS &= \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} \cdot q_0^2 - \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} \cdot q_0 \cdot q_1 \\ &\quad + \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} \cdot q_4^2 - \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} \cdot q_3 \cdot q_4. \end{aligned}

Интеграл жёсткости для каждого симплекса вычисляется как интеграл квадрата производной пробной функции. Поскольку производные пробных функций являются константами, интегрирование сводится к умножению на длину симплекса.

Согласно локальной формуле из раздела «Матрица жёсткости 1D», для каждого симплекса этот интеграл равен $\frac{1}{l_{(i)(i+1)}}$ .

Общий интеграл жёсткости по всей области $M$ равен сумме интегралов по всем симплексам:

\begin{aligned} \int_M (\nabla \upsilon)^2 \,dM &= \sum_{i=1}^{4} \int_{s_i} \left(\upsilon'(x)\right)^2 \,dx \\ &= \frac{\displaystyle 1}{\displaystyle l_{(0)(1)}} \left(q_0^2 - 2q_0q_1 + q_1^2\right) + \frac{\displaystyle 1}{\displaystyle l_{(1)(2)}} \left(q_1^2 - 2q_1q_2 + q_2^2\right) \\ &\quad + \frac{\displaystyle 1}{\displaystyle l_{(2)(3)}} \left(q_2^2 - 2q_2q_3 + q_3^2\right) + \frac{\displaystyle 1}{\displaystyle l_{(3)(4)}} \left(q_3^2 - 2q_3q_4 + q_4^2\right). \end{aligned}

Интеграл демпфирования для каждого симплекса вычисляется как интеграл квадрата пробной функции.

Согласно локальной формуле из раздела «Матрица демпфирования 1D», интеграл квадрата пробной функции на симплексе равен $\frac{l_{(i)(i+1)}}{3}\left(q_i^2 + q_i q_{i+1} + q_{i+1}^2\right)$ .

Общий интеграл демпфирования по всей области $M$ равен сумме интегралов по всем симплексам:

\begin{aligned} \int_M \upsilon^2 \,dM &= \sum_{i=1}^{4} \int_{s_i} \upsilon^2(x) \,dx \\ &= \frac{\displaystyle l_{(0)(1)}}{\displaystyle 3} \left(q_0^2 + q_0q_1 + q_1^2\right) + \frac{\displaystyle l_{(1)(2)}}{\displaystyle 3} \left(q_1^2 + q_1q_2 + q_2^2\right) \\ &\quad + \frac{\displaystyle l_{(2)(3)}}{\displaystyle 3} \left(q_2^2 + q_2q_3 + q_3^2\right) + \frac{\displaystyle l_{(3)(4)}}{\displaystyle 3} \left(q_3^2 + q_3q_4 + q_4^2\right). \end{aligned}

Интеграл нагрузки для каждого симплекса вычисляется как интеграл произведения функции источника на пробную функцию.

Согласно локальной формуле из раздела «Вектор нагрузки 1D», на крайних симплексах $s_1$ и $s_4$ интеграл равен нулю (источник там тождественно нулевой), а на $s_2$ и $s_3$ — соответствующему вкладу локального вектора нагрузки.

Общий интеграл нагрузки по всей области $M$ равен сумме интегралов по всем симплексам:

\begin{aligned} \int_M \widetilde{f}(x) \cdot \upsilon(x) \,dM &= \sum_{i=1}^{4} \int_{s_i} \widetilde{f}(x) \cdot \upsilon(x) \,dx \\ &= 0 + \frac{\displaystyle l_{(1)(2)}}{\displaystyle 6} \left(2f_1 q_1 + f_1 q_2 + f_2 q_1 + 2f_2 q_2\right) \\ &\quad + \frac{\displaystyle l_{(2)(3)}}{\displaystyle 6} \left(2f_2 q_2 + f_2 q_3 + f_3 q_2 + 2f_3 q_3\right) + 0. \end{aligned}

Упростив:

\begin{aligned} \int_M \widetilde{f}(x) \cdot \upsilon(x) \,dM &= \frac{\displaystyle l_{(1)(2)}}{\displaystyle 6} \left(2f_1 q_1 + f_1 q_2 + f_2 q_1 + 2f_2 q_2\right) \\ &\quad + \frac{\displaystyle l_{(2)(3)}}{\displaystyle 6} \left(2f_2 q_2 + f_2 q_3 + f_3 q_2 + 2f_3 q_3\right). \end{aligned}

Запишем все полученные интегралы в матричной форме. Введём вектор неизвестных $\vec{q} = \begin{pmatrix} q_0 & q_1 & q_2 & q_3 & q_4 \end{pmatrix}^T$ .

Интеграл жёсткости можно записать в матричном виде:

\int_M (\nabla \upsilon)^2 \,dM = \vec{q}^T \cdot K \cdot \vec{q},

где $K$ — матрица жёсткости размером $5 \times 5$ :

\begin{aligned}K = \begin{pmatrix} \frac{1}{l_{(0)(1)}} & -\frac{1}{l_{(0)(1)}} & 0 & 0 & 0 \\ -\frac{1}{l_{(0)(1)}} & \frac{1}{l_{(0)(1)}} + \frac{1}{l_{(1)(2)}} & -\frac{1}{l_{(1)(2)}} & 0 & 0 \\ 0 & -\frac{1}{l_{(1)(2)}} & \frac{1}{l_{(1)(2)}} + \frac{1}{l_{(2)(3)}} & -\frac{1}{l_{(2)(3)}} & 0 \\ 0 & 0 & -\frac{1}{l_{(2)(3)}} & \frac{1}{l_{(2)(3)}} + \frac{1}{l_{(3)(4)}} & -\frac{1}{l_{(3)(4)}} \\ 0 & 0 & 0 & -\frac{1}{l_{(3)(4)}} & \frac{1}{l_{(3)(4)}} \end{pmatrix}.\end{aligned}

Матрица $K$ является симметричной и трёхдиагональной.

Интеграл демпфирования можно записать в матричном виде:

\int_M \upsilon^2 \,dM = \vec{q}^T \cdot D \cdot \vec{q},

где $D$ — матрица демпфирования размером $5 \times 5$ :

\begin{aligned}D = \begin{pmatrix} \frac{l_{(0)(1)}}{3} & \frac{l_{(0)(1)}}{6} & 0 & 0 & 0 \\ \frac{l_{(0)(1)}}{6} & \frac{l_{(0)(1)}}{3} + \frac{l_{(1)(2)}}{3} & \frac{l_{(1)(2)}}{6} & 0 & 0 \\ 0 & \frac{l_{(1)(2)}}{6} & \frac{l_{(1)(2)}}{3} + \frac{l_{(2)(3)}}{3} & \frac{l_{(2)(3)}}{6} & 0 \\ 0 & 0 & \frac{l_{(2)(3)}}{6} & \frac{l_{(2)(3)}}{3} + \frac{l_{(3)(4)}}{3} & \frac{l_{(3)(4)}}{6} \\ 0 & 0 & 0 & \frac{l_{(3)(4)}}{6} & \frac{l_{(3)(4)}}{3} \end{pmatrix}.\end{aligned}

Матрица $D$ также является симметричной и трёхдиагональной.

Интеграл нагрузки можно записать в матричном виде:

\int_M \widetilde{f}(x) \cdot \upsilon(x) \,dM = \vec{F}^T \cdot \vec{q},

где $\vec{F}$ — вектор нагрузки размером $5 \times 1$ :

\begin{gathered}\vec{F} = \begin{pmatrix} 0 \\ \frac{2f_1 \cdot l_{(1)(2)}}{6} + \frac{f_2 \cdot l_{(1)(2)}}{6} \\ \frac{f_1 \cdot l_{(1)(2)}}{6} + \frac{2f_2 \cdot l_{(1)(2)}}{6} + \frac{2f_2 \cdot l_{(2)(3)}}{6} + \frac{f_3 \cdot l_{(2)(3)}}{6} \\ \frac{f_2 \cdot l_{(2)(3)}}{6} + \frac{2f_3 \cdot l_{(2)(3)}}{6} \\ 0 \end{pmatrix}.\end{gathered}

Упростив, получим:

\begin{gathered}\vec{F} = \begin{pmatrix} 0 \\ \frac{l_{(1)(2)}}{6} \left(2f_1 + f_2\right) \\ \frac{f_1 \cdot l_{(1)(2)}}{6} + \frac{f_2 (l_{(1)(2)} + l_{(2)(3)})}{3} + \frac{f_3 \cdot l_{(2)(3)}}{6} \\ \frac{l_{(2)(3)}}{6} \left(f_2 + 2f_3\right) \\ 0 \end{pmatrix}.\end{gathered}

Интеграл по границе можно записать в матричном виде:

\int_S \upsilon \cdot \nabla \upsilon \,dS = \vec{q}^T \cdot B \cdot \vec{q},

где $B$ — матрица граничных условий размером $5 \times 5$ :

\begin{aligned}B = \begin{pmatrix} \frac{1}{l_{(0)(1)}} & -\frac{1}{2l_{(0)(1)}} & 0 & 0 & 0 \\ -\frac{1}{2l_{(0)(1)}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -\frac{1}{2l_{(3)(4)}} \\ 0 & 0 & 0 & -\frac{1}{2l_{(3)(4)}} & \frac{1}{l_{(3)(4)}} \end{pmatrix}.\end{aligned}

Матрица $B$ является симметричной и содержит ненулевые элементы только в углах, соответствующих граничным узлам $q_0$ и $q_4$ .

Стационарное уравнение

Согласно (I.5), в матричной форме имеем:

a^2 \cdot \vec{q}^T \cdot K \cdot \vec{q} - 2 \cdot \vec{F}^T \cdot \vec{q} - a^2 \cdot \vec{q}^T \cdot B \cdot \vec{q} \rightarrow \min.

Вычисляя производную по $\vec{q}$ , получаем систему линейных уравнений:

a^2 (K - B) \cdot \vec{q} = \vec{F}.

Нестационарное уравнение

Согласно (I.6), в матричной форме имеем:

\begin{aligned} &\vec{q}_n^T \cdot D \cdot \vec{q}_n + \Delta t \cdot a^2 \cdot \vec{q}_n^T \cdot K \cdot \vec{q}_n - 2 \cdot \Delta t \cdot \vec{F}^T \cdot \vec{q}_n \\ &- 2 \cdot \vec{q}_{n-1}^T \cdot D \cdot \vec{q}_n - \Delta t \cdot a^2 \cdot \vec{q}_n^T \cdot B \cdot \vec{q}_n \rightarrow \min. \end{aligned}

Вычисляя производную по $\vec{q}_n$ , получаем систему линейных уравнений:

\left[D + \Delta t \cdot a^2 (K - B)\right] \cdot \vec{q}_n = \Delta t \cdot \vec{F} + D \cdot \vec{q}_{n-1}.

Обозначения:

$K$ — матрица жёсткости $5 \times 5$ ,
$D$ — матрица демпфирования $5 \times 5$ ,
$B$ — матрица граничных условий $5 \times 5$ ,
$\vec{F}$ — вектор нагрузки $5 \times 1$ ,
$\vec{q}$ — вектор неизвестных $5 \times 1$ для стационарного случая,
$\vec{q}_n, \vec{q}_{n-1}$ — векторы неизвестных на текущем и предыдущем временных слоях,
$a$ — коэффициент теплопроводности,
$\Delta t$ — временной шаг.

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

Вычислим разность матриц $K - B$ , которая входит в итоговые уравнения для стационарного и нестационарного случаев:

\begin{aligned}K - B = \begin{pmatrix} 0 & -\frac{1}{2l_{(0)(1)}} & 0 & 0 & 0 \\ -\frac{1}{2l_{(0)(1)}} & \frac{1}{l_{(0)(1)}} + \frac{1}{l_{(1)(2)}} & -\frac{1}{l_{(1)(2)}} & 0 & 0 \\ 0 & -\frac{1}{l_{(1)(2)}} & \frac{1}{l_{(1)(2)}} + \frac{1}{l_{(2)(3)}} & -\frac{1}{l_{(2)(3)}} & 0 \\ 0 & 0 & -\frac{1}{l_{(2)(3)}} & \frac{1}{l_{(2)(3)}} + \frac{1}{l_{(3)(4)}} & -\frac{1}{2l_{(3)(4)}} \\ 0 & 0 & 0 & -\frac{1}{2l_{(3)(4)}} & 0 \end{pmatrix}.\end{aligned}

(J.1)

Матрица $K - B$ остаётся симметричной и трёхдиагональной. Характерной особенностью является то, что диагональные элементы на границах (первый и последний) обнуляются, что отражает вклад граничных условий в систему.

Хитрость этой матрицы в том, что диагональные элементы на границах обнулились, и из самой системы их не найти. Вы скажете: ну ок, есть же граничные условия — именно. Чтобы довести пример до конца, осталось разобрать способ учёта граничных условий в системе; этому посвящены следующие приложения.

I. Функционал метода Бубнова-Галёркина K. Учёт граничных условий Дирихле