Перемещения и деформации в пластине

Министерство образования и науки России

Федеральное государственное бюджетное 
образовательное учреждение высшего образования
 «Сыктывкарский государственный университет имени Питирима Сорокина»
(ФГБОУ ВО «СГУ им. Питирима Сорокина»)
Институт точных наук и информационных технологий
Кафедра математического моделирования и кибернетики



Допустить к защите
Зав. кафедрой математического
моделирования и кибернетики
___________________________ к.ф.-м.н. Беляев Ю.Н. 
«____»______________ 2017 г.









Выпускная квалификационная работа

Метод обобщённой реакции











Исполнитель: обучающийся 149 группы Грабовский А.А.







Научный руководитель: к.ф.-м.н., доцент, Ермоленко А.В.

















Сыктывкар 2017
Содержание

Введение………………………………………………………………………..3

Обзор теории пластин…………………………………………………………4
2.1.  Основные уравнения теории упругости…...…...……………………….4
2.2.  Основные понятия и гипотезы…………………………………………..6
2.3.  Перемещения и деформации в пластине………………………………..7
2.4.  Напряжения в пластине………………………………………………….9
2.5.  Усилия в пластине………………………………………………………..9
2.6.  Дифференциальное уравнение изогнутой срединной поверхности
        пластины…………………………………………………………………11

Обзор необходимых сведений из численных методов…………………….13
3.1.  Численное дифференцирование………………………………………..13
3.1.  Построение обратной матрицы……………………………………...…14

Метод обобщённой реакции…………………………………………………16

Задача о расчёте пластины по классической теории с постоянной нагрузкой……………………………………………………………………...17

Задача о нахождении прогиба прямоугольной пластины с переменной нагрузкой, используя обратную матрицу…………………………………..20
6.1.  Постановка задачи………………………………………………………20
6.2.  Решение по классической теории……………………………………...20
6.3.  Решение по теории Кармана……………………………………….......20

Заключение…………………………………………………………………...23

Источники…………………………………………………………………….24

Приложение…………………………………………………………………..25




Введение

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

Объект исследования – это прямоугольная пластина, жёстко закреплённая по контуру. Предметом исследования является прогиб пластины при малых деформациях.

Цель – создать программное обеспечение для построения  поверхностей функций прогиба и реакции прямоугольной пластины по классической  и обобщённой теориям.

Задачи:

Изучение литературы

Нахождение способа вычисления численного решения

Составление итерационных схем для классической и обобщённой теорий

Выбор языка программирования

Написание программного кода

Тестирование программы

Сравнение результатов, построенных по классической и  обобщённой теории для одинаковых входных данных

Введение раскрывает актуальность, определяет объект, предмет,  цель и задачи исследования. Во втором пункте выводится уравнение Софи–Жермен–Лагранжа, используя основные уравнения теории упругости. В третьем пункте затрагиваются некоторые сведения из численных методов. В четвёртом пункте описывается метод обобщённой реакции. В пятом пункте рассматривается численное решение уравнения  Софи–Жермен–Лагранжа для постоянной нагрузки. В шестом пункте поясняется нахождение численного решения по классической и обобщённой теориям. В заключении сравниваются результаты и делаются выводы о результатах расчётов.




2. Обзор теории пластин

2.1 Основные уравнения теории упругости

Рассмотрим основные группы формул из теории упругости.

Группа статических уравнений.

В эту группу входят дифференциальные уравнения равновесия:

                                                                                (1.1)

и условия на поверхности (граничные условия):

                                                                                 (1.2)

Группа геометрических уравнений.

В эту группу входят формулы Коши:

                                                                                 (1.3)

Группа физических уравнений.

В эту группу входят формулы закона Гука:

                                                                    

                                                                  (1.4)

                                     

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

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

шесть компонент тензора напряжений:




шесть компонент тензора деформаций:





три вектора перемещения:



Для нахождения этих неизвестных функций в нашем распоряжении есть пятнадцать уравнений: три дифференциальных уравнения равновесия (1.1), шесть формул Коши (1.3) и шесть формул Гука (1.4). Таким образом, с математической точки зрения задача может быть решена и сводится к интегрированию этих пятнадцати уравнений при удовлетворении условий на поверхности (1.2).

Решения указанных уравнений находятся в зависимости от того, какие величины приняты за основные неизвестные:

Решение в перемещениях, если основными неизвестными являются три вектора перемещения:



Решение в напряжениях, если основными неизвестными являются шесть компонент тензора напряжений:




Решение в смешанной форме, если основными неизвестными являются некоторые векторы перемещений и некоторые компоненты тензора напряжений.







2.2 Основные понятия и гипотезы

Пластиной называется тело, ограниченное двумя параллельными плоскостями, расстояние между которыми, называемое толщиной пластины , мало по сравнению с его другими размерами.

Плоскость, делящая пластину пополам по толщине, называется срединной плоскостью. При изгибе пластины срединная плоскость превращается в изогнутую срединную поверхность пластины.

Линия пересечения боковой поверхности пластины со срединной плоскостью называется контуром пластины.

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

Тонкими называются пластины, имеющие отношения толщины к наименьшему характерному размеру в следующих пределах:



и величину ожидаемых прогибов не более ? h.

Пластины, у которых , рассчитываются по теории толстых плит, а пластины, имеющие прогибы более ? h рассчитываются по геометрически нелинейной теории гибких пластин или мембран.

Теория тонких пластин основана на следующих гипотезах, предложенных Кирхгофом:

Гипотеза прямых нормалей: любой линейный элемент, нормальный к срединной плоскости пластины, остаётся прямолинейным и нормальным к срединной поверхности после деформации и длина его не изменяется.

Любой линейный элемент, нормальный к срединной плоскости, направлен вдоль оси z, и, следовательно, первая часть гипотезы предполагает, что прямые углы между этим элементом и осями x и y остаются прямыми, т.е. сдвиги в указанных плоскостях отсутствуют:

                                                                                         (1.5)

Допущение о сохранении длины прямолинейного элемента предполагает, что линейная деформация в направлениях оси z отсутствует:

                                                                                                          (1.6)

Гипотеза о нерастяжимости срединной плоскости: в срединной плоскости отсутствуют линейные относительные деформации и деформации сдвига, т.е. срединная плоскость является нейтральной. Следовательно, в срединной плоскости перемещения отсутствуют

                                                                                             (1.7)

Гипотеза об отсутствии давления между слоями пластины. Ввиду малости давления между слоями пластины, параллельными срединной плоскости, напряжением  по сравнению с напряжениями  и  можно пренебрегать.





2.3 Перемещения и деформации в пластине

Рассмотрим пластину с нормальной нагрузкой к срединной плоскости. Под действием этой нагрузки пластина получит перемещения. Для их определения воспользуемся гипотезами Кирхгофа.

Согласно первой гипотезе линейная деформация в направлении оси z равна нулю (1.6). Подставляя это условие в третью формулу Коши (1.3) получаем:



откуда следует, что прогибы пластины w не зависят от координат z, т.е.



Это означает, что все точки пластины, лежащие на одной вертикали, получают одинаковые прогибы. Следовательно, достаточно определить прогибы срединной плоскости пластины, чтобы узнать прогибы всех её точек.

Рассматривая условия для сдвигов (1.5), из формул Коши (1.3) получаем:

                                                         ,

                                                         ,

откуда находим производные от составляющих перемещения u и v по координате z:

                                                                 ,

                                                                 .

Интегрируя эти уравнения по z, получаем:

                                                      .                                 (1.8)

Для вычисления функций  и  воспользуемся гипотезой о нерастяжимости срединной плоскости. Согласно этой гипотезе составляющие перемещения  и  на срединной плоскости при z=0 равны нулю. Подставляя эти условия в формулы (1.8), получаем:

                                                            ,

                                                            .

Тогда формулы (1.8) примут следующий вид:

                                                               .                                          (1.9)

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

Составляющие деформации в пластинке, отличные от нуля, найдем с помощью формул Коши (1.3), подставляя в них значения составляющих перемещения (1.9):

                                                                      (1.10)





2.4 Напряжения в пластине

Для вычисления нормальных напряжений  и  возьмём две первые формулы закона Гука (1.4) и на основании третьей гипотезы отбросим напряжение  по сравнению с напряжениями  и . Тогда получим:

                                              ,

                                             

откуда с учётом формул (1.10) находим:

                                              ,                                 (1.11)

                                              ,                                 (1.12)

Четвёртая формула закона Гука (1.4) после подстановки угловой деформации  из формул (1.10) примет такой вид:

                                       .                                   (1.13)





2.5 Усилия в пластине

Исследуем, какие усилия создаются напряжениями в сечениях пластины, нормальных к её срединной плоскости. Подсчитаем изгибающий момент. Обозначим через  погонный (приходящийся на единицу ширины сечения) изгибающий момент в сечении с нормалью x. Изгибающий момент в рассматриваемом сечении создаётся нормальными напряжениями . Равнодействующая этих напряжений на площадке толщиной dz и шириной, равной единице, равна , а изгибающий момент . Суммируя моменты от напряжения  на всех таких площадках по толщине пластины, получаем выражение для погонного изгибающего момента в сечении с нормалью х:

                                            .

Подставляя сюда значение нормального напряжения  из формулы (1.11) и вынося за знак интеграла величины, не зависящие от координаты z, находим:

                                         .

После интегрирования получаем:

                                                   .                                (1.14)

Входящая сюда величина

                                                                                                       (1.15)

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

Аналогично найдём погонный изгибающий момент в сечении с нормалью y, подставляя формулу (1.12):

                                           .                               (1.16)

Погонный крутящий момент в сечении с нормалями х и y равен:

                                              .

После подстановки касательного напряжения  из формулы (1.13)  и интегрирования находим:

                                                   .                                     (1.17)

2.6 Дифференциальное уравнение изогнутой срединной поверхности пластины

Вырежем из срединной плоскости пластины бесконечно малый элемент dxdy и покажем действующие на него нагрузки (Рис. 1). 



Рис. 1

Рассматриваемый элемент срединной плоскости находится в равновесии, следовательно, должно выполняться условие равновесия, т.е. сумма всех внешних сил равна нулю и сумма всех моментов внешних сил равна нулю. А именно шесть условий равновесия: три уравнения проекций на координатные оси и три уравнения моментов относительно этих осей.

Спроектируем все силы на ось z:

  .

После упрощения получаем:

                                                       .                                          (1.18)

Уравнение моментов всех сил относительно оси y дает:

                        

                        

                         .

После упрощения получаем:

                                                      .                                            (1.19)

Аналогично из уравнения моментов относительно оси х получаем:

                                                      .                                            (1.20)

Исключим из уравнений (1.18) – (1.20) поперечные силы . В результате получим:

                                           .                                  (1.21)

Подставив в полученное уравнение (1.21) моменты из формул (1.14), (1.16) и (1.17), найдем:

              .

Откуда после упрощения получим:

                                       ,                                   (1.22)

или

                                                          .                                               (1.23)

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

Уравнение (1.23) должно быть дополнено граничными условиями. Условия на контуре пластины зависят от характера закрепления ее краёв.



3. Обзор необходимых сведений из численных методов

Численное дифференцирование

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

    .     (3.1)

Пусть определены значения  функции  в точках сетки. Существуют три основные разностные отношения: левая конечная разность

                             ,                          (3.2)

правая конечная разность

                             ,                          (3.3)

и центральная конечная разность

                        .                       (3.4)

Все остальные разностные отношения выводятся из комбинаций уравнений (3.2) – (3.4). Например, получим разностное отношение для производной второго порядка

                       .

Выпишем разностные отношения, которые будем использовать в дальнейшем

                                                                                       (3.5)

                                                                                       (3.6)

                                                              (3.7)

                                                                 (3.8)

                                                                 (3.9)

                   

                                                (3.10)



Построение обратной матрицы

Квадратная матрица  порядка  называется обратной, если существует невырожденная матрица  того же порядка такая, что выполняется условие

 ,

где  – единичная матрица порядка. Невырожденность матрицы означает, что её определитель не равен нулю. Обозначается следующим образом . Обратные матрицы не существуют для неквадратных матриц. 

Далее нам понадобиться понятие элементарного преобразования строк (столбцов), к которым относятся: перестановка местами любых двух строк матрицы, умножение на ненулевую константу любой строки матрицы и прибавление к любой строке матрицы другой строки, умноженной на ненулевое число. Есть различные способы нахождения обратной матрицы. Рассмотрим три основных метода.

В методе Жордана – Гаусса нужно взять исходную матрицу  и единичную матрицу  одной размерности. Затем с помощью элементарных преобразований строк нужно привести матрицу  к единичному виду. Те же самые операции необходимо применить и к единичной матрице . В результате матрица  станет единичной, а матрица  – обратной. Сложность алгоритма оценивается как .

Во втором методе используется следующая формула

                                                         ,                                             (3.11)

где  – присоединённая матрица, т.е. матрица, которая получается с помощью алгебраических дополнений для соответствующих элементов транспонированной матрицы. Сложность алгоритма -  .

	Третьим методом является LU – разложение. Методом исключения Гаусса расширенная матрица  приводится к диагональному виду. В результате получаем расширенную матрицу , где  – это верхняя треугольная матрица с единицами на главной диагонали, а  – нижняя треугольная матрица. При умножении этих матриц получается исходная матрица . Обратная матрица находится по следующей формуле

                                                .                                          (3.12)

После умножения уравнения (3.17) на , а затем на  получим два уравнения

                                                       ,                                               (3.13)

                                                       .                                               (3.14)

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




Метод обобщённой реакции

Метод обобщённой реакции решения контактных задач со свободной границей был предложен и обоснован в Сыктывкарском государственном университете Михайловским Евгением Ильичом. Этот метод представляет собой итерационную схему, общий вид которой для уравнения Софи–Жермен–Лагранжа записывается следующим образом

                                                 (4.1)

                                                                    (4.2)

где  – это функция Грина,  – некоторая константа,  – зазор между пластиной.

                                                        (4.3)

Формула (4.2) называется положительной срезкой функции, которой эквиваленты, следующие условия

                                       – условие непроникновения,                     (4.4)

                                – условие односторонности связи,                    (4.5)

               – условие дополняющей нежёсткости.           (4.6)

Аналогично для системы Кармана можно записать итерационную схему в следующем виде

                    (4.7)

                                 ,                                (4.8)

где

                 (4.9)

                                                                  (4.10)

                                                    (4.11)



Задача о расчёте пластины по классической теории c постоянной нагрузкой

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

                                                ,                                                 (5.1)

где w – искомый прогиб в любой точке пластины;  – нормальная нагрузка;  – цилиндрическая жесткость пластины; h–толщина пластины; 
? – коэффициент Пуассона; E–модуль Юнга; 
 – бигармонический оператор.

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

                    .              (5.2)

Решение исходной краевой задачи (5.2) будем отыскивать в виде таблицы значений в точках сетки

                                        .                                  (5.3)

На пластину нанесена сетка. Для аппроксимации производных будем использовать формулы (3.8)–(3.10).

Аппроксимируем условия закрепления пластины, воспользовавшись формулами (3.2) – (3.3)



(5.4)

(5.4)

Подставим в краевую задачу (5.2) формулы (3.2)–(3.10). Получим разностную схему:

 .                                                                          (5.5)

Уравнение (5.5) можно представить в виде сеточного шаблона, изображённого на рисунке (5.1).



Рис. 5.1

Этот шаблон нужно применить ко всем внутренним узлам (узлы нумеруются по порядку слева направо и сверху вниз), при этом центр шаблона должен совпадать с номером узла, к которому он применяется. Коэффициенты для узлов берутся из разностного уравнения (5.5). Например, для первого узла при  и 9 внутренних узлах получим первую строку матрицы следующего вида

                                           ,

для второго узла 2-ю строку

                                        .

Таким образом, применив данный шаблон к каждому внутреннему узлу сетки, получим матрицу коэффициентов A размера . Например для 9 внутренних узлов данная матрица имеет следующий вид:



A   = 

A   = 

где

                                              

                                              

                                              

 Решив систему  можно найти значение прогиба в любом узле сетки.


Задача о нахождении прогиба прямоугольной пластины с переменной нагрузкой, используя обратную матрицу

Постановка задачи

Аналогично пункту 5 требуется найти прогиб прямоугольной пластины жёстко закреплённой по контуру. Искать решение будем с помощью обратной матрицы, которая определяется по методу Жордана–Гаусса.





Решение по классической теории

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

Итерационная схема, записанная в виде формул (4.1) – (4.3) может быть представлена в следующем виде

                                                                                         (6.1)

                                                                    (6.2)                         

                                                                                       (6.3)





Решение по теории Кармана

Система Кармана записывается в следующем виде

                       ,                         (6.4)

где  – это функция напряжений Эри, получающаяся из условий

                                           .                        (6.5)

Граничные условия функции напряжений записываются аналогично функции прогиба

                                 (6.6)

Из формулы (6.6) после аппроксимации получается



(6.7)

(6.7)                                         

Запишем итерационную схему, используя формулы (4.7) – (4.11)
                                                                        (6.8)

                                                                (6.9)  (6.10)

                                                             (6.11)

В формулах (6.10) и (6.11) вторые частные производные от некоторой функции заменяются матричными записями

                                (6.12)

                                                                              (6.13)

                                                                              (6.14)

где  или  и



(6.15)

(6.15),



(6.16)

(6.16),



(6.17)

(6.17).

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




Заключение

Была написана программа для построения поверхностей функций прогиба и реакции прямоугольной пластины с помощью языка программирования Python 3 и библиотеки matplotlib. Эта программа вызывает скомпилированный код, написанный на языке программирования C++, чтобы создать JSON-файлы для соответствующих функций, по которым в дальнейшем и строятся данные поверхности. Для создания JSON-файла использовалась библиотека rapidjson. Графический интерфейс программы был написан с помощью библиотеки PyQt 5 и изображён на Рис. 7.1.



Рис. 7.1

Если указанное число итераций равно нулю, то итерационная схема выполняется до тех пор, пока выполнено следующее условие

                                             .                                                 (7.1)

При вычислении нормы была использована евклидова норма. В качестве  было выбрано значение . 




Источники

Михайловский	 Е.И., Бадокин К.В., Ермоленко А.В. Теория изгиба пластин типа Кармана без гипотез Кирхгофа // Вестн. Сыктывкарского ун-та Сер. 1. Вып. 3. 1999

Львов Г.И. Основы теории пластин и оболочек // Учебные пособия: Учебник - , 2014.

Теребушко О.И. Основы теории упругости и пластичности.–М: Наука. Главная редакция физико-математической литературы, 1984.–320 с.

Кузнецова Е.В. Изгиб пластин: Учебное-методическое пособие к решению задач и лабораторному практикуму по исследованию прогибов при нагружении прямоугольных и круглых пластин / Пермь. Гос. Техн. ун-т. Пермь, 2006. 32с.

Ефимов В.В. Динамика и прочность авиационных конструкций. Часть I: учебное пособие. - М.: МГТУ ГА, 2013. – 72 с., 29 ил., лит.: 9 наим.

А. А. Лахтин. Расчёт прямоугольной пластины методом конечных разностей: конспект лекций. Екатеринбург:   Изд–во УрГУПС, 2013. 15с.

А.А. Самарский, А.В. Гулин. Численные методы: Учеб. пособие для вузов. – М.: Наука. Гл. ред. физ-мат. лит., 1989. 432 с.

Прохорёнок Н.А., Дронов В.А. Python 3 и PyQt 5. Разработка приложений.–Спб.: БХВ-Петербург, 2016.–832 с.: ил.

https://matplotlib.org/gallery.html

https://github.com/miloyip/rapidjson

https://docs.python.org/3/




Приложение

Заголовочный файл C++ Helper.h:

#pragma once
#include "iostream"
using namespace std;


class Helper {

public:
    Helper(void); ~Helper(void);
    void set_args(int argc, char* argv[], int &n, double &a, double &b, 
        double &h, double &q0, double &E, double &nu);

    void init_matr_idx(int n, int** &matr_idx);

    void init_A1(int n, double hx, double hy, int** &matr_idx, double** &A1);

    void init_A2(int max_idx, double** &A2);
    void init_Axx_Ayy_Axy(int n, double hx, double hy, int** &matr_idx,
        double** &Axx, double** &Ayy, double** &Axy);

    void front_process(int max_idx, double** &A1, double** &A2);

    void back_process(int max_idx, double** &A1, double** &A2);
    void init_w_r(int n, double* &w, double* &r, double** &W, double** &R);
    void save_matrix_to_json(int n, double** matrix, string matrix_name);

};


Вспомогательный класс C++ Helper.cpp для нахождения обратной матрицы методом Жордана – Гаусса:

#include "Helper.h"

#include "math.h"
#include "rapidjson/writer.h"
#include "rapidjson/stringbuffer.h"
#include "fstream"
using namespace rapidjson;

Helper::Helper(void) {}

Helper::~Helper(void) {}

void Helper::set_args(int argc, char* argv[], int &n, double &a, double &b,
        double &h, double &q0, double &E, double &nu) {
    if (argc > 1) { n = atoi(argv[1]); } if (argc > 2) { a = atof(argv[2]); }
    if (argc > 3) { b = atof(argv[3]); } if (argc > 4) { h = atof(argv[4]); }
    if (argc > 5) { q0 = atof(argv[5]); } if (argc > 6) { E = atof(argv[6]); }
    if (argc > 7) { nu = atof(argv[7]); }
}



void Helper::init_matr_idx(int n, int** &matr_idx) {
    int k = 1;
    for (int i = 2; i <= n - 2; i++)
        for (int j = 2; j <= n - 2; j++) { matr_idx[i][j] = k++; }
}


void Helper::init_A1(int n,double hx,double hy,int** &matr_idx,double** &A1) {
    double hx4 = pow(hx, 4), hy4 = pow(hy, 4), hx2hy2 = pow(hx, 2) * pow(hy, 2);
    int k = 0;
    for (int i = 2; i <= n - 2; i++) {
        for (int j = 2; j <= n - 2; j++) {
            A1[k][matr_idx[i][j] - 1] = 6 / hx4 + 6 / hy4 + 8 / hx2hy2;
            if (matr_idx[i - 2][j] != 0) { A1[k][matr_idx[i - 2][j] - 1] = 1 / hy4; }
            if (matr_idx[i + 2][j] != 0) { A1[k][matr_idx[i + 2][j] - 1] = 1 / hy4; }
            if (matr_idx[i][j - 2] != 0) { A1[k][matr_idx[i][j - 2] - 1] = 1 / hx4; }
            if (matr_idx[i][j + 2] != 0) { A1[k][matr_idx[i][j + 2] - 1] = 1 / hx4; }
            if (matr_idx[i-1][j] != 0) A1[k][matr_idx[i-1][j]-1] = -4 / hy4 -4 / hx2hy2;
            if (matr_idx[i+1][j] != 0) A1[k][matr_idx[i+1][j]-1] = -4 / hy4 -4 / hx2hy2;
            if (matr_idx[i][j-1] != 0) A1[k][matr_idx[i][j-1]-1] = -4 / hx4 -4 / hx2hy2;
            if (matr_idx[i][j+1] != 0) A1[k][matr_idx[i][j+1]-1] = -4 / hx4 -4 / hx2hy2;
            if (matr_idx[i-1][j-1] != 0) A1[k][matr_idx[i - 1][j - 1] - 1] = 2 / hx2hy2;
            if (matr_idx[i-1][j+1] != 0) A1[k][matr_idx[i - 1][j + 1] - 1] = 2 / hx2hy2;
            if (matr_idx[i+1][j-1] != 0) A1[k][matr_idx[i + 1][j - 1] - 1] = 2 / hx2hy2;
            if (matr_idx[i+1][j+1] != 0) A1[k][matr_idx[i + 1][j + 1] - 1] = 2 / hx2hy2;
            k++;
}}}

void Helper::init_A2(int max_idx, double** &A2) {
    for (int i = 0; i < max_idx; i++) { A2[i][i] = 1.0; }
}



void Helper::init_Axx_Ayy_Axy(int n, double hx, double hy, int** &matr_idx,
        double** &Axx, double** &Ayy, double** &Axy) {
    double hx2 = pow(hx, 2), hy2 = pow(hy, 2); int k = 0;
    for (int i = 2; i <= n - 2; i++) {
        for (int j = 2; j <= n - 2; j++) {
            Axx[k][matr_idx[i][j] - 1] = -2 / hx2; Ayy[k][matr_idx[i][j] - 1] = -2 / hy2;
            Axy[k][matr_idx[i][j] - 1] = 1 / (hx * hy);
            if (matr_idx[i - 1][j] != 0) {
                Ayy[k][matr_idx[i - 1][j] - 1] = 1 / hy2;
                Axy[k][matr_idx[i - 1][j] - 1] = -1 / (2 * hx * hy); }
            if (matr_idx[i + 1][j] != 0) {
                Ayy[k][matr_idx[i + 1][j] - 1] = 1 / hy2;
                Axy[k][matr_idx[i + 1][j] - 1] = -1 / (2 * hx * hy); }
            if (matr_idx[i][j - 1] != 0) {
                Axx[k][matr_idx[i][j - 1] - 1] = 1 / hx2;
                Axy[k][matr_idx[i][j - 1] - 1] = -1 / (2 * hx * hy); }
            if (matr_idx[i][j + 1] != 0) {
                Axx[k][matr_idx[i][j + 1] - 1] = 1 / hx2;
                Axy[k][matr_idx[i][j + 1] - 1] = -1 / (2 * hx * hy); }
            if (matr_idx[i - 1][j + 1] != 0) { 
                Axy[k][matr_idx[i - 1][j + 1] - 1] = 1 / (2 * hx * hy); }
            if (matr_idx[i + 1][j - 1] != 0) { 
                Axy[k][matr_idx[i + 1][j - 1] - 1] = 1 / (2 * hx * hy); }
            k++;
}}}

void Helper::front_process(int max_idx, double** &A1, double** &A2) {
    double tmp;
    for (int k = 0; k < max_idx; k++) {
        if (A1[k][k] == 0.0) {
            for (int j = 0; j < max_idx; j++) {
                tmp = A1[k][j]; A1[k][j] = A1[k + 1][j]; A1[k + 1][j] = tmp;
                tmp = A2[k][j]; A2[k][j] = A2[k + 1][j]; A2[k + 1][j] = tmp;
        }}
        if (A1[k][k] != 0.0 && A1[k][k] != 1.0) {
            tmp = A1[k][k];
            for (int j = k; j < max_idx; j++) { A1[k][j] /= tmp; }
            for (int j = 0; j < max_idx; j++) { A2[k][j] /= tmp; }
            if (k == max_idx - 1) { break; }
        }
        for (int i = k + 1; i < max_idx; i++) {
            tmp = A1[i][k];
            for (int j = k; j < max_idx; j++) { A1[i][j] -= A1[k][j] * tmp; }
            for (int j = 0; j < max_idx; j++) { A2[i][j] -= A2[k][j] * tmp; }
}}}

void Helper::back_process(int max_idx, double** &A1, double** &A2) {
    double tmp;
    for (int k = max_idx - 1; k > 0; k--) {
        for (int i = k - 1; i > -1; i--) {
            tmp = A1[i][k];
            for (int j = k; j > i; j--) { A1[i][j] -= A1[k][j] * tmp; }
            for (int j = 0; j < max_idx; j++) { A2[i][j] -= A2[k][j] * tmp; }
}}}

void Helper::init_w_r(int n, double* &w, double* &r, 
        double** &W, double** &R) {
    int k = 0;
    for (int i = 2; i <= n - 2; i++) {
        for (int j = 2; j <= n - 2; j++) {
            W[i][j] = w[k]; R[i][j] = r[k]; k++;
}}}

void Helper::save_matrix_to_json(int n, double** matrix, string matrix_name) {
    ofstream fout((matrix_name + ".json").c_str());
    StringBuffer sb; Writer writer(sb);
    writer.StartObject(); writer.Key(matrix_name.c_str()); writer.StartArray();
    for (int i = 0; i < n + 1; i++) {
        writer.StartArray();
        for (int j = 0; j < n + 1; j++) { writer.Double(matrix[i][j]); }
        writer.EndArray(); }
    writer.EndArray(); writer.EndObject();
    fout << sb.GetString() << endl; fout.close();
}

Файл C++ ct_const.cpp для построения JSON-файла матрицы прогиба по классической теории с постоянной нагрузкой:
#include "math.h"
#include "Helper.h"

int **matr_idx = NULL; double **A1 = NULL, **A2 = NULL;
double *_w = NULL; double **w = NULL;

void alloc_arrays(int n, int max_idx, int** &matr_idx,
        double** &A1, double** &A2, double* &_w, double** &w) {
    matr_idx = new int*[n + 1]; w = new double*[n + 1];
 .......................