4.4.2. Дискретные представления поверхностных и объемных интегралов

За исключением нескольких простых задач, найти аналитическое решение уравнений (4.15) и (4.16) невозможно, и поэтому нужно было разработать численные методы решения. Подчеркнем, что (4.15) и (4.16) являются точным представлением решения задачи. При численном решении любая ошибка в конечном результате

получается исключительно из-за дискретизации интегралов и последующего решения алгебраических уравнений.

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

Рис. 4.3.

Следует отметить, что, хотя внутренние ячейки имеют такой же вид, как в схеме дискретизации, используемой в методе конечных элементов, они представляют лишь удобный способ вычисления влияния распределенных по объему сил Как было указано ранее, порядок системы линейных алгебраических уравнений, получающейся при численной реализации МГЭ, определяется лишь количеством использованных граничных элементов (т. е. подобластей и совершенно не связан с количеством использованных внутренних ячеек (т. е. подобластей V). Более того, дискретизация объемными ячейками в МГЭ может быть совершенно произвольной в том смысле, что она не обязательно соответствует поверхностной дискретизации. Впрочем, для упрощения задачи (особенно с точки зрения подготовки данных) обычно удобно согласовать разбиение границы и внутренней области.

Если предположить, что в пределах каждого элемента имеют постоянные значения, то для граничного элемента можно написать

где — координаты некоторой характерной точки граничного элемента, например его середины, длина граничного элемента, поверхность 1-й ячейки. Используя матричные обозначения, эти уравнения можно записать в более удобной форме:

где причем эти индексы введены лишь для удобства. Если необходимо ввести линейно меняющиеся в пределах каждого граничного элемента 9 и то величины и нельзя выносить из-под знака интеграла, как это было сделано в (4.19) и (4.20).

Рис. 4.4.

Однако если представить и в виде функции узловых значений с использованием базисных функций , то, как показано на рис. 4.4, для граничного элемента с узлами и 1-й ячейки с узлами можно написать

Здесь

матрица размером причем определяют, как показано выше, точку приложения мгновенной нагрузки, длина граничного элемента;

Поэтому для лежащей на границе типичной узловой точки уравнения (4.19) и (4.20) можно переписать в виде

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

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

усилия определяются введением вблизи угла двух узловых точек, взятых на произвольно близком расстоянии одна от другой (обычно на расстоянии, составляющем 0.05 длины элемента) и представляющих в пределе крайние точки двух поверхностей (рис. 4.5).

Рис. 4.5.

Некоторые другие особенности решения задач для областей с ребрами и узлами на границе приводятся в гл. 7.

Уравнения (4.22) и (4.23) можно использовать для получения окончательной системы алгебраических уравнений, но сначала нужно вычислить интегралы, стоящие в квадратных скобках. Если точка поля не лежит где-либо внутри нагруженного граничного элемента (в противном случае интегралы становятся сингулярными), то это делается просто по формуле интегрирования Гаусса (см. приложение В).

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

и

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

где функции выраженные в координатах

В полярной системе координат (рис. 4.6) интегралы в правой части (4.24) и (4.25) легко вычислить точно.

Уравнения (4.24) и (4.25) можно записать следующим образом:

где — результаты интегрирования

Таким образом, получаем

Рис. 4.6.

где уравнения (4.26) и (4.27) верны для любых точек за исключением совпадающих точек поля и точек приложения нагрузки (т. е. при соответственно, как на рис. 4.7, а). В этом случае вычисления можно провести двумя способами:

1. Поместить точку поля вблизи точки приложения нагрузки с внутренней стороны как показано на рис. 4.7 (обычно равно 0.01 длины граничного элемента).

Рис. 4.7.

В этом случае следует отбросить член так как он аппроксимируется непосредственно (из-за близости

2. Поместить точку поля на поверхности и затем вычислить интегралы в соответствии с изложенным в гл. 3.

Разница в значениях коэффициентов, полученных в тестовых расчетах обоими способами, оказалась совершенно незначительной.

В отношении членов в (4.22) и (4.23), содержащих объемные интегралы, следует отметить, что существует три типа интегралов, которые нужно вычислять.

1. Интегралы, в которых точка поля не совпадает с вершиной внутренней ячейки и не лежит на ее стороне (рис. 4.8, в), являются интегралами от непрерывных функций и поэтому вычисляются при помощи формулы Гаусса.

2. Интегралы, в которых точка поля совпадает с вершиной внутренней ячейки, вычисляются с помощью четырехточечной формулы интегрирования Гаусса (рис. 4.8, а).

3. Если точка поля лежит на стороне треугольной ячейки, то сначала ячейка делится на две части прямой, проходящей через точку поля (рис. 4.8, б), и получающиеся интегралы вычисляются по формуле интегрирования Гаусса.

Рис. 4.8.

Затем уравнения (4.22) и (4.23) применяются по очереди к каждому узлу, что при заданных граничных условиях приводит к системе линейных алгебраических уравнений для определения векторов в каждом узле (ниже этот процесс объясняется подробно). Как только все значения найдены, смещения, деформации и напряжения во внутренних точках вычисляются по следующим формулам:

Матрицы имеют размеры а вектор-столбцы и 5 имеют размерность 2, 3 и 3 соответственно, например

Для получения системы уравнений интегралы в уравнениях (4.28) — (4.30) можно вычислить по схеме, аналогичной предложенной выше.

Одна особенность НМГЭ заключается в том, что можно вычислить тензор напряжений во внутренней точке, используя уравнение (4.14) дважды [9]. Сначала ядро вычисляется при выборе внешней нормали в виде единичного вектора в направлении оси (т. е. ), что дает Затем, выбирая (единичная нормаль в направлении оси получаем

Если необходимо вычислить напряжения и деформации для большого числа внутренних точек, то для вычисления точных значений смещений в узловых точках внутренних ячеек более эффективно использовать уравнение (4.28). Это уравнение справедливо во всех точках, включая граничные. Напряжения и деформации вычисляются по смещениям при помощи процедуры, аналогичной используемой в методе конечных элементов или в методе конечных разностей. Поэтому, если шестимерные векторы смещений в узлах некоторой внутренней треугольной ячейки, то вектор смещений в любой точке этой внутренней ячейки записывается в виде [см. (4.21)]

Деформации вычисляются подстановкой (4.31) в (4.3) и могут быть записаны в виде

где матрица размером Соответствующие этим деформациям напряжения получаются по формуле

где матрица упругих постоянных размером

Главный недостаток такого алгоритма заключается в невозможности использования коэффициента Пуассона большего 0.49 для представления несжимаемых материалов Основанная, на уравнениях (4.28) — (4.30) полная процедура МГЭ позволяет конечно, рассматривать любое допустимое значение