Алгоритмы и программные модули моделирования геометрически нелинейного стержневого конечного элемента

Программные модули используются в процедурах статического и динами­ческого анализа геометрически нелинейных стержневых систем. Подпрограм­мы вычисления матрицы касательной жесткости и вектора упругих узловых сил реакции СКЭ вызываются на каждой итерации метода Ньютона-Рафсона (Нью­тона-Канторовича) [63] в задачах статического расчета, или же вызываются в необходимой последовательности на итерациях полного или модифицирован­ного методов Ньютона-Рафсона при численном интегрировании уравнений движения методом Ньюмарка.

Матрица касательной жесткости СКЭ вычисляется подпрограммой MatrKTan.

CALL MatrKTan( Schema, AL, DispL, EA, Elz, Ely, E, Ip,

Git, Indx, Indy, Indz, bAxForce, Tmp, MLen, N, KTM)

Описание параметров данной подпрограммы дано в табл.2.3

Таблица 2.3. Параметры подпрограммы MatrKTan

В подпрограмме MatrKTan вычисляется матрица геометрической жестко­

сти и матрица жесткости, учитывающая большие перемещения узлов элемента 129

по формулам (2.50) и (2.48). Осевая сила, входящая в (2.50), вычисляется по уточненной нелинейной формуле (2.51), причем имеется возможность вычис­лять ее по общепринятой линеаризованной формуле, которая имеет вид

где- производные функций формы осевых перемещений узлов в на­

чале и конце стержневого элемента по осевой координате х, u_iи и _j- осевые перемещения узлов в начале и конце стержня в его локальной системе коорди­нат.

Интегрирование по дайне элемента в формулах (2.48) и (2.50) выполнялось чис­ленно с использованием квадратурных формул Гаусса-Лежандра в нормализованных оптимальных координатах. Порядок подынтегральных полиномов в указанных фор­мулах для каждого элемента матриц зависит от вида расчетной схемы (плоская или пространственная, число локальных узловых степеней свободы в элементе), способа закрепления концов стержневого элемента (выбора соответствующих функций фор­мы) и способа вычисления осевой силы. Максимальный порядок полиномов - 8, одна­ко, большое количество элементов матриц имеют порядок полиномов 6 или 4, кото­рый может еще далее снижаться при разных условиях.

Поэтому для повышения эффективности вычислений был использован оп­тимальный алгоритм численного интегрирования матриц по длине СКЭ. Вы­брано два множества весовых коэффициентов и нормализованных координат Гаусса-Лежандра с 3-мя и 6-ю оптимальными точками позволяющими точно интегрировать полиномы 6 и 12 порядка. Исходя из условий закрепления кон­цов и способа вычисления осевой силы, в подпрограмме определяется макси­мальный порядок полинома элементов матрицы геометрической жесткости, что позволяет выбрать схему интегрирования с 3-мя или 6-ю точками.

Интегрирование осуществляется только по ненулевым элементам матрицы геометрической жесткости, которые определяются по виду расчетной схемы. Интегрирование матрицы [к] (2.48) в подпрограмме осуществляется поэле­ментно. Исходя из вышеназванных условий, для каждого элемента матрицы [к] определяется максимальный порядок полинома, что позволяет выбрать схему интегрирования с наименьшими вычислительными затратами. В некото­

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

Вектор упругих узловых сил реакции СКЭ вычисляется подпрограммой VectFN. Формат обращения к этой подпрограмме имеет вид:

CALL VectFN( Schema, AL, DispL, EA, Elz, Ely, E, Ip,

Git, Indx, Indy, Indz, N, FN )

где FN (N) - выходной параметр: вектор упругих узловых сил реакции стерж­невого конечного элемента. (REAL* 8). Описание других параметров подпро­граммы дано в табл.2.3.

В подпрограмме VectFN вектор узловых сил реакции вычисляется по фор­муле (2.45). Численное интегрирование по длине элемента выполняется с по­мощью квадратурных формул Гаусса-Лежандра по 6-ти оптимальным точкам. Максимальный порядок подынтегрального полинома в формуле (2.45) равен 8.

Входным параметром подпрограмм MatrKTan и VectFN является вектор текущих перемещений узлов элемента в его текущей локальной CK - DispL. Для вычисления этого вектора используется методика учета больших поворо­тов и перемещений узлов в глобальной СК, разработанная в главе 2.

Необходимо отметить, что положение и ориентация локальной CK стерж­невого элемента изменяются каждый раз, когда изменяются координаты его уз­лов при деформировании стержня, а также при перемещении его как абсолютно твердого тела. Для вычисления вектора перемещений и поворотов узлов стерж­ня в его локальной системе координат служит подпрограмма FNDDISP. Фор­мат обращения к этой подпрограмме следующий:

CALL FNDDISP( Schema, NCoord, NLSIter, BarVZ, Coordl, Coord2, NodeCosl, NodeCos2, BarCos, Tl, T2, ТЗ, T4, DispL, N, NBvz, IErr ) Описание параметров дано в табл.2.4 (см. также табл.2.3)

Таблица 2.4. Параметры подпрограммы FNDDISP

Подпрограмма FNDDISP вычисляет текущую ориентацию локальной CK СКЭ с помощью матриц ориентации узлов в начале и конце элемента и матри­цы ориентации локальной CK стержня в его недеформированном состоянии.

CALL NLNodesRecalc( Schema, NPoint, NCoord, NLSIter, NDFStr,

CoordIt, NodeCos, Fix, Degree, DisVecIt, IErr)

Описание параметров подпрограммы дано в табл.2.5 (и табл.2.3, 2.4)

Таблица 2.5. Параметры подпрограммы NLNodesRecalc

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

При выполнении статических расчетов внешнюю узловую нагрузку на сис­тему необходимо постепенно увеличивать от нуля до заданного значения в те­чение нескольких итераций. Это позволяет удовлетворить условие малости приращений перемещений и поворотов на одной итерации решения системы нелинейных уравнений (2.38), а также позволяет достичь лучшей сходимости итераций и получить более точное решение. Для вычисления вектора невязки внутренних и внешних узловых сил в пакет «COMPASS» был включен про­граммный модуль AsmbNRF. Модуль формирует глобальный вектор внутрен­них сил системы и определяет невязку внутренних и внешних узловых СИЛ C учетом пошагового увеличения значений внешних сил. При этом формула 133

где і - номер текущей итерации; N_li- параметр, равный числу шагов пропор­ционального приращения внешней нагрузки от нуля до заданного уровня; N_cs- параметр, равный числу итераций между приращениями внешней нагрузки. Иными словами, внешняя нагрузка постепенно увеличивается за N_uшагов на величинуна каждом шаге, и после каждого приращения нагрузки вы­

полняется N_csитераций для получения промежуточного решения.

2.8.