Алгоритмы эффективной организации взаимосвязей упругой и гидродинамической составляющих конечноэлементного моделирования упругогидродинамики узлов трения Текст научной статьи по специальности «Физика»

Машиностроение и машиноведение

УДК 519.3: 621.822

DOI: 10.12737/article 59bllcba23e7a6.8666367G

М.В.Зернин, Н.Н.Рыбкин

АЛГОРИТМЫ ЭФФЕКТИВНОЙ ОРГАНИЗАЦИИ ВЗАИМОСВЯЗЕЙ УПРУГОЙ И ГИДРОДИНАМИЧЕСКОЙ СОСТАВЛЯЮЩИХ КОНЕЧНОЭЛЕМЕНТНОГО МОДЕЛИРОВАНИЯ УПРУГОГИДРОДИНАМИКИ УЗЛОВ ТРЕНИЯ

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

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

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

M.V. Zernin, N.N. Rybkin

EFFICIENT ORGANIZATION ALGORITHMS OF INTERCONNECTIONS OF ELASTIC AND HYDRODYNAMIC CONSTITUENTS OF FINITE ELEMENT SIMULATION OF FRICTION UNIT ELASTIC HYDRODYNAMICS

Basic algorithms decreasing labor intensity in the solution of a problem connected with elastic hydrodynamics are formulated. Some options for the formation of cost-effective grids of finite elements are presented. It is shown the efficient use of the method of condensation in order that in a non-linear procedure only the units located on operation surfaces of a friction unit should be considered. A cost-effective option for the realization of a condensation method is de-

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

В БГТУ развиваются [1-5] конечно-элементные методы расчета гидродинамики динамически нагруженных подшипников скольжения (ПС) на основе уравнений Рейнольдса. Выполняется подготовка к реализации третьей версии, в которой будут учтены неоднородное поле температур, а также упругие и термоупругие перемещения деталей ПС на основе трехмерного ко-нечноэлементного моделирования [6; 7]. Здесь рассмотрены особенности взаимосвязи гидродинамической задачи (определение параметров течения смазывающей жидкости в зазоре) с задачей упругого деформирования деталей узла трения. Отметим, что в такой упругогидродинамиче-ской задаче (УГД-задаче) появляется дополнительная нелинейность именно из-за

scribed in a short way. The algorithm of the reduction of pressures distributed on a surface to the system of concentrated unit loads is stated thoroughly. The computation algorithm of points displacements of a surface in the direction perpendicular to it is shown.

Key words: elastic hydrodynamics problems, elastic deformations, hydrodynamic pressures, finite element method, grid concentration, reduction method, condensation method, approximating polynomial.

необходимости определения параметров взаимосвязи двух задач в итерационном процессе. Задача упругая - линейна при заданных граничных условиях, задача гидродинамическая - нелинейна и требует организации трех вложенных друг в друга итерационных процедур [5]. Для динамически нагруженных ПС цикл их нагруже-ния разбит на малые временные шаги (обычно связанные с углом относительного поворота деталей). Этот расчет также завершается не за один цикл, так как критерием его корректного завершения является замыкание траектории движения центра вала относительно центра ПС.

Если решать только гидродинамическую задачу (отыскание поля давлений в масляном слое при изотермических усло-

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

упруго деформируются от действия давлений. Требуется многократное решение указанных задач в трехмерной постановке. Фактически гидродинамическая задача должна быть решена многократно, и трудоемкость решения такой УГД-задачи возрастает в десятки раз.

Варианты снижения трудоемкости за счет

Для компьютерной реализации алгоритма должны быть предусмотрены меры по повышению его эффективности, в том числе за счет выбора рационального типа конечных элементов (КЭ) для температурной и упругой (термоупругой) задач. В нашей статье [8] показано, что за счет применения специальных КЭ с квазисовместными аппроксимациями можно снизить трудоемкость решения ТУГД-задачи более чем на порядок. Эффективность применения квазисовместных двадцати-узловых КЭ заключается в точном моделировании полей искомых параметров

рациональных конечноэлементных сеток

сравнительно малым количеством КЭ с большим соотношением сторон (до 100:1 и более). Например, конечноэлементные сетки вала и деталей ПС могут быть такими, как показано на рис. 1а. Но экономичное представление упругих деталей малым количеством КЭ может усложнить учет граничных условий кавитации. Граница области давлений (сжимающих) и области их отсутствия (кавитации) может проходить не по узлам сетки, а по граням достаточно больших КЭ. Учет таких граничных условий возможен, но достаточно сложно реализуется.

Рис. 1. Конечноэлементные сетки деталей подшипника скольжения: а - при использовании квазисовместных конечных элементов; б - при сгущении сетки; в - при применении приема редуцирования для согласования

различающихся сеток

Альтернативные варианты - существенное сгущение сетки КЭ в зоне контактирования деталей (рис. 1 б) или использование метода редуцирования. В последнем случае отдельные конструктивные элементы (баббитовый слой на рис. 1в) или фрагменты деталей (тонкий поверхностный слой вала на рис. 1в) разбиты на мелкие КЭ, а основная часть деталей - на более крупные КЭ. Сетки с различными по размерам КЭ корректно сшиваются по ме-

тоду редуцирования [9]: требуемая гладкость функции на границе обеспечивается за счет соответствующего согласования полиномов, аппроксимирующих эту функцию по разные стороны от границы. Алгоритмы автоматической генерации сеток переменной густоты и алгоритмы редуцирования реализованы в промышленных программных пакетах, например в БЕМАР [10]. Поэтому можно подготавливать требуемые сетки в промышленных программ-

ных пакетах, а процесс итерационного решения УГД-задачи реализовать как от-

Снижение трудоемкости за счет применения

Еще один прием, позволяющий значительно снизить трудоемкость задачи, также часто используется и называется «конденсация». Заключается он в том, что можно выделить группы узлов, которые необходимы для организации итерационных процедур УГД-задачи. Такие узлы расположены на рабочих поверхностях ПС, так как именно податливость поверхности и теплоотвод через поверхности из слоя масла вглубь деталей ПС участвуют во взаимосвязях составляющих задач. Такие узлы в методе конденсации называют внешними или граничными (обозначены индексом все остальные - внутренними (обозначены индексом \п).

Метод конденсации состоит в том, что исключают внутренние узлы. Но исключают их таким образом, что оставшиеся в расчетной схеме граничные узлы отображают свойства всей конструкции (со всеми узлами). И после решения УГД-задачи можно рассмотреть исходные разбивки со всеми узлами и КЭ и определить в них все требуемые характеристики напряженно-деформированного состояния (НДС). Процедура конденсации основана на следующем представлении разрешающей системы линейных алгебраических уравнений (СЛАУ) МКЭ:

№ }=М. (1)

- к к -

vn-vn . vn-gr

к к

_ gr-vn. gr-gr _

дельный внешний модуль.

процедуры конденсации

Если все внутренние узлы с перемещениями {{та} сгруппировать вверху, а все граничные узлы с перемещениями {{} внизу, то систему (1) можно представить в виде следующих блоков:

Такое представление системы является условным. В реальных системах внутренние и граничные узлы расположены вперемешку, и преобразование системы к виду (2) может существенно снизить эффективность хранения матрицы жесткости. Обычно матрица жесткости является ленточной и симметричной. Перестановка элементов этой матрицы может существенно увеличить ширину ленты, возрастут требуемые для хранения матрицы объемы памяти.

Простейшим преобразованием из уравнения (2) можно получить два отдельных уравнения:

{ }+[К^г КЛ={Р }, (3)

К™ К КЫр! (4)

Выразив {{уи} из уравнения (3) и подставив результат в уравнение (4), после преобразований получим:

Ик 1-1"^ 1Гк HK 11/ }={P 1-K 1Гк Хг{р }

U- gr-gr Л V gr-vn X vn-vn J L vn-gr W\ gr) V gr ) L gr-vn± vn-vnj t vn ) .

(5)

Из условия симметрии матрицы жесткости можно записать:

\_Kgr-vn ] = \Kvn-gr ] . Одно из распространенных названий метода конденсации - «метод суперэлементов». Поэтому окончательно можно записать разрешающее уравнение для суперэлемента (индекс эв) в следующем виде:

\К„ Ю=Р»}, (6)

где искомый вектор перемещений {{е }= }; матрица коэффициентов су-

перэлемента

Гк 1 = "к 1_Гк I1 Гк Нк 1

L sel L gr- gr \ L vn - gr J L vn-vn J L vn - gr J

vn - gr J 5

частей

вектор правых

{Pse}= {Pgr}" ГKvn - gr f ГKvn-vn 1"' {Pn }.

Следует отметить, что сокращение количества неизвестных не приводит к ухудшению точности решения задачи. Система уравнений (6) сохраняет все свойства системы уравнений (1). После определения по (6) значений граничных пере-

мещений }={Уя;.} их можно подставить в уравнение (3) и вычислить значения внутренних перемещений {}. Последующее определение параметров НДС здесь не обсуждается.

Остановимся на вариантах реализации метода конденсации. Выше приведено много матричных уравнений, и может сложиться впечатление, что трудоемкость реализации метода конденсации достаточно высока, так как требуется выполнение серии матричных операций. Однако фактически операции метода конденсации являются операциями метода Гаусса решения СЛАУ. Так, получение уравнения (6) из уравнения (1) является незавершенной (неполной) процедурой прямого хода метода Гаусса. Незавершенность состоит в том, что исключены только внутренние узлы (степени свободы). Причем такую процедуру исключения авторы реализовали без изменения формы хранения исходной матрицы СЛАУ (1). После выполнения таких действий можно выбрать коэффициенты матрицы суперэлемента \Ке ] и вектора нагрузок на узлы суперэлемента {Р\е}.

После определения перемещений узлов суперэлемента можно продолжить решение СЛАУ (1) с уже известными значениями }={ия,.}. Фактически это неполная реализация обратного хода метода Гаусса, а именно определение {'иш} при уже известных (из решения уравнения (6)) перемещениях {и:ее}={ияг}.

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

ритмов решения нелинейных задач МКЭ сводятся к методу последовательных приближений и на каждой итерации необходимо решать СЛАУ. Метод конденсации позволяет выделить те уравнения, которые не участвуют в нелинейной схеме, и исключить эти уравнения, выполнив неполный прямой ход метода Гаусса. В нелинейной схеме остается система уравнений значительно меньших размеров. Далее эту нелинейную систему уравнений можно разрешить путем многократного решения СЛАУ. Причем такие СЛАУ необязательно решать методом Гаусса. Многие исследования показали, что эффективнее применять метод квадратного корня Холецко-го [6]. На этом этапе тоже можно экономить. Если требуется многократное решение СЛАУ с одной и той же матрицей коэффициентов, а изменяются только правые части СЛАУ, то в нелинейной процедуре требуется многократное преобразование правых частей и определение неизвестных. Начальные процедуры (обычно наиболее трудоемкие) можно выполнить один раз. Так, при использовании метода Гаусса прямой ход (факторизацию матрицы) можно выполнить один раз. При использовании метода Холецкого один раз выполняется разложение матрицы коэффициентов на верхнюю треугольную и нижнюю треугольную.

Если применять метод конденсации, то в решении УГД-задачи задействованы только узлы рабочих поверхностей узла трения. Трудоемкость может быть снижена на порядки.

Алгоритм моделирования взаимосвязи у]

Между конечноэлементными сетками, используемыми для решения гидродинамической задачи и упругой задачи, непосредственной взаимосвязи нет. Для примера на рис. 2б показана развертка масляного слоя крейцкопфного ПС судового дизеля. На рис. 2а приведены следующие обозначения параметров, используемых при решении ГД-задачи: ё -диаметр вала; Б - диаметр рабочей поверхности

угой и гидродинамической задач

подшипника (расточки); х, у, 2 - декартовы координаты в плоскости, перпендикулярной оси подшипника, и вдоль оси подшипника; в - угловая координата; а - скорость вращения вала в подшипнике; вх и ву -проекции эксцентриситета (удаления центра вала от центра расточки подшипника) на координатные оси; Гх и Гу - проекции на координатные оси нагрузки, действующей на вал. На рис. 2а можно выделить две

системы координат: декартову - х, у, г и цилиндрическую - в, z. Для наглядного представления полученных результатов часто используют также условную систему координат, отличающуюся от реальной цилиндрической системы координат отсчетом угловой координаты вдоль прямой линии. На рис. 2б вместо декартовой коор-

динаты г использована угловая координата в, т.е. показана так называемая развертка масляного слоя. Иначе сетку КЭ, моделирующую масляный слой, изобразить невозможно. Результаты решения ГД-задачи также нагляднее демонстрировать на развертке.

а) б)

Рис. 2. Система координат и геометрические параметры гидродинамической задачи (а) и конечноэлементная сетка развертки масляного слоя крейцкопфного подшипника судового дизеля (б)

Масляный слой дискретизирован на треугольные конечные элементы первого порядка. Вблизи канавок реализовано сгущение сетки. Варианты конечноэлемент-ных сеток для упругих деталей этого же ПС показаны на рис. 1а,б. Сетки упругой и ГД задач не взаимосвязаны. Результатами решения ГД-задачи являются зазоры (толщины масляного слоя h(6,z)) и давления в масляном слое p(0,z). Результатами решения упругой задачи являются в первую очередь перемещения узлов сетки КЭ. При необходимости можно определить другие параметры НДС. Взаимосвязь гидродинамической и упругой задач реализуется посредством учета давлений и изменения геометрии деталей (изменение зазоров между их поверхностями): полученные в результате решения гидродинамической задачи давления должны быть учтены в качестве силовых граничных условий в упругой задаче; полученные в результате решения упругой задачи перемещения рабочих поверхностей должны быть учтены в следующей итерации решения гидродинамической задачи. Основная идея учета податливости заключается в дополнении зависимости для расчета зазоров при не-деформируемых поверхностях ПС слагаемым, учитывающим деформацию поверхностей:

h = Const — ex cos в — ey sin в + L(p) , (7)

где первые три слагаемых определяют зазор (толщину слоя масла) из геометрии цилиндрического ПС с недеформируемы-ми поверхностями (рис. 2а), а Ь(р) - оператор податливостей, связывающий значения упругого радиального перемещения одной или обеих поверхностей с действующими давлениями. В нашем случае значения перемещений определяются из расчетов по МКЭ упругих деталей.

Рассмотрим подробнее, как следует реализовать такие взаимосвязи, на примере цилиндрических ПС. Эпюра давлений (даже на развертке рабочей поверхности) может быть достаточно сложной геометрической фигурой (в зависимости от наличия канавок, отверстий и т.п. элементов). Для примера на рис. 3 а приведена эпюра давлений в крейцкопфном ПС с системой канавок для подачи масла (сетка КЭ приведена на рис. 2б). На рис. 3б приведена эпюра давлений в ПС турбины с эллиптической (лимонной) расточкой при перекосе осей вала и ПС. На этой эпюре давлений имеются две зоны гидродинамического контакта (два масляных клина), расположенные на верхнем и нижнем вкладышах. В обеих эпюрах (рис. 3а,б) можно выделить участки с высокими градиентами давлений. В декартовой и цилиндрической системах координат эпюра давлений будет выглядеть столь сложно, что изобразить ее

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

1. Эпюра давлений - всегда гладкая поверхность, хотя могут быть большие градиенты вблизи канавок и других подобных элементов.

2. Давления всегда направлены по нормали к рабочим поверхностям ПС.

а) б)

Рис. 3. Примеры эпюр давлений: а - в крейцкопфном подшипнике скольжения с системой масляных канавок; б - в турбинном подшипнике скольжения с эллиптической расточкой при перекосе осей

Учитывая эти положения, можно принять, что на грань каждого КЭ, принадлежащую участку рабочей поверхности, действует по нормали эпюра, которую можно описать полиномом (рис. 4а). Любая точка грани удалена от центра декартовой системы координат х-у-2 на расстояние радиуса-вектора я . Учесть эту эпюру в упругом расчете можно посредством приведения такой распределенной по поверхности грани нагрузки к системе нагрузок, приложенных к узлам этой грани, по методике, приведенной в работах [6; 7]. Для этого (рис. 4б) на грани с распределенной нагрузкой вводится локальная

двухмерная система криволинейных координат (5 , а 8 узлов, расположенные на этой грани, повторно нумеруются (первый раз они нумеровались в составе всех 20 узлов КЭ). На декартовых осях откладываются орты е , е2 и еъ, а в точке центра грани - единичная нормаль п. Также изображены орты я,5 и я производных

радиуса-вектора я по локальным координатам % и п. Здесь и далее знак запятой в нижнем индексе означает операцию дифференцирования по указанной следом координате.

а) б) в) г)

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

Координаты точек этой грани в трехмерной декартовой системе координат можно выразить [6; 7] через соответству-

ющие координаты восьми узлов и восемь двухмерных функций формы для этих узлов:

= £ N (5,4)хг; .у = £ N (5,п)У ; * = £ N (5,п)г,

(8)

X

1=1

1=1

2=1

где функции формы взяты из таблицы.

Таблица

Функции формы восьмиузлового элемента и их производные по локальным координатам

Номера узлов / Функции формы N ($,п) Производная по £ N ($,п) Производная по ц N ,($,п) $

1, 2, 3, 4 (1+$$, )(1 + пп, )($$, +ПП-1)/4 (2$$,2 + 2$$Пп + $, П2П2)/4 (2пП +$$Л + 2$$, П + $$ П,2)/4 +1 +1

5, 7 (1 -$2)(1 + пп, )/2 -$(1 + пп) (1 -$2)П /2 0 +1

6, 8 (1 + $$, )(1 -п2)/2 $ (1 -П2)/2 (1 + $$)п +1 0

Далее используются производные от глобальных координат точки грани по локальным координатам £ и щ. Формулы этих

8 8

производных получаем в результате дифференцирования формул (8).

х

$

= е; у,$ = е^$($,п)у; 2,$ = е2;

,=1 8

,=1 8

п)у, ; 2,п=Еып($л)2 ,.

хп = ЕNп($л)х,; у,п = ЕNп($ ,=1 ,=1

Производные от функций формы по локальным координатам приведены в таблице.

В общем случае распределенная нагрузка, приложенная к грани элемента, может быть направлена произвольно и характеризоваться тремя эпюрами: рх ($,п), р^ ($,п) и р2($,п) .Так как в УГД-задаче в

качестве распределенной по грани нагрузки выступает гидродинамическое давле-

ео8(п,е) = А/л/ЕО - ¥2 ; ео8(п,е2) = 5/л/ЕО - ¥2 ; ео8(п,е3) = С/л/ЕО - ¥2 где константы правой части вычисляются по формулам

(9)

(10)

ние, всегда направленное по нормали к поверхности, то его можно характеризовать одной эпюрой р($,п). Но для определения проекций давления на глобальные оси координат х, у, г требуется вычислять направляющие косинусы внешней нормали к точкам грани КЭ. В [6] приведен подробный вывод соответствующих формул. Здесь приведем их окончательный вид в удобной для реализации форме:

(11)

А=у, $2, п- У, п2, $;

Е = х$+ у$ + 2$ ;

5 = х,п2,$ х,$2,п; 2,2,2 О = хп+ У,п + 2п';

С=х,$у,п- хпУ$;

¥=х, $х, п+ у, $у,

п + 2,$2,п .

Подставив в формулы (9) и (10) значения координат узлов, можно вычислить в этих конкретных точках поверхности производные от декартовых координат по локальным координатам. Далее можно вычислить параметры А, В, С, Е, G и ¥ в этих точках.

Компоненты вектора сил, действующие на /-й узел, статически эквивалентные действующей по поверхности грани нагрузке р($,п), можно вычислить по формуле [6]

Р =<

рг

Рх

рг

РУ

рг

Р1

£ £ ЕзN ($,п)р($,п)

(12)

где Е3 - единичная матрица размером 3*3.

Интегралы (12) можно взять аналитически только при простейших видах функции давлений р($,п). Так, в случае приложения к плоской прямоугольной

грани площадью А¥ постоянного распределенного нормального давления р($,п) =p=const суммарная сила Р=рА¥ разносится по восьми узлам грани следующим образом. К узлам, лежащим на цен-

8

,=1

трах сторон (номера 5, 6, 7, 8), прикладываются усилия Р/3, направленные так же, как и давления р. К узлам, расположенным в углах грани (номера 1, 2, 3 , 4), прикладываются усилия Р/12, направленные против действующего давления р.

При всех более сложных вариантах функции давлений интегралы (12) следует определять численно по методу Гаусса -Лежандра:

Р -И

Р1

Рх

Р2

РУ

Р.

Р*

Р.

= £з N (5,1) р(5,1)

]=1 к=1

(13)

где Е3 - единичная матрица размером 3*3; п - количество точек интегрирования по каждому из двух направлений; § и Щк - координаты точек интегрирования; Н) и Нк -соответствующие весовые коэффициенты. Для трех компонент вектора усилий, действующих на /-й узел, можно записать:

РХ. = N (5,1) Р(5,1) А

]=\ к=1

РУ = N 5,1)р(5,1)в

;=1 к=1

А к ;

п п

5!НАк

РР*. = N (5,1)р(5,1)с

]=1 к=1

Порядок интегрирования (и соответствующее число точек интегрирования) выбирается в зависимости от порядка по-динтегральных функций. Используются функции формы N (5,1) второго порядка. Они умножаются на функции давлений р(5,1), которые тоже имеют порядок не ниже второго. Рассмотрим выбор параметров численного интегрирования подробнее на основе [6; 7]. Используют метод одномерного численного интегрирования, известный как квадратура Гаусса - Ле-жандра. В этом методе для известной степени полинома т интегрируемой функции подбираются такое число точек интегри-

рования и такие значения их координат, чтобы интегрирование выполнялось точно. Значения весовых коэффициентов получаются из решения системы уравнений. В [7] показано, что существует такое соответствие степени интегрируемого полинома и количества требуемых точек интегрирования: при т=1 п=1; при т=3 п=2; при т=5 п=3; при т=7 п=4. Значения степени полинома т всегда нечетные. Если нужно интегрировать полином четной степени, количество точек интегрирования следует выбирать для ближайшего большего нечетного.

В случае интегрирования по нескольким локальным координатам применяются мультипликативные формулы. По каждому из направлений применяется одномерный вариант. На рис. 5 приведена информация о координатах и весовых коэффициентах для п=2, 3, 4 в двухмерном случае. В зависимости от показателя степени полинома по обоим независимым направлениям и тп выбирается количество точек интегрирования, обычно соответствующее п2. Так, при ,тщ < 3 п2=4 (рис. 5а), при ,тл < 5 п2=9 (рис. 5б) и при ,тл < 7 п2=16 (рис. 5в).

а) б) в)

Рис. 5. Расположение точек интегрирования и значения весовых коэффициентов:

а -т% ,тп < 3; б - ,тп < 5; в - ,тп < 7

п

п

п

п

В качестве тестовых задач рассчитаны задачи об изгибе консольно закреплённой балки (рис. 6). Если приложена постоянная распределённая нагрузка интенсивностью q (рис. 6а), то аналитическое решение для упругой линии является уравнением четвертой степени относительно продольной координаты. Если к консольно закреплённой балке приложена линейно убывающая распределённая нагрузка (рис. 6б), то уравнение упругой линии является

уравнением пятого порядка относительно продольной координаты. Если к консольно закреплённой балке приложена квадратично убывающая распределённая нагрузка (рис. 6в), то уравнение упругой линии является уравнением шестой степени относительно продольной координаты. Результаты расчетов по МКЭ с нагрузкой, приведенной к узлам (по описанной выше методике), соответствуют аналитическим решениям.

а) б) в)

Рис. 6. Изгиб консольно закрепленной балки под действием распределенной нагрузки: а - равномерной; б - линейной; в - параболической

Варианты получения формул, аппроксимирующих давления

В квадратуре Гаусса - Лежандра для численного интегрирования (12) используется формула для гидродинамических давлений р($,п). Причем эта формула должна быть полиномиальной, так как в зависимости от степени этого полинома выбирается система точек интегрирования. В результате решения ГД-задачи на треугольной сетке КЭ имеем значения давлений в узлах. В пределах каждого треугольного КЭ реализована билинейная интерполяция давлений. На основе этих данных необходимо построить полиномиальную аппроксимацию для эпюры давлений. Далее следует выделить отдельные полиномиальные формулы для давлений на грань каждого из КЭ упругой задачи, поверхность которого является частью рабочей поверхности узла трения. Возможны различные варианты решения таких задач.

В первом варианте строится двухмерный аппроксимирующий полином давлений на развертке поверхности трения, например по методике [11]. Далее при рассмотрении отдельных граней КЭ упругой задачи по полученной аппроксимирующей формуле вычисляются значения давлений

в узлах КЭ упругой задачи. Постулируется, что степень полинома р($,п) не превышает порядка аппроксимации, используемой для геометрии упругого КЭ. Локальная формула для давлений на грань такого КЭ строится с применением функций формы упругого КЭ: 8

р = ЕN ($,п)р,. (14)

,=1

Этот вариант применим в тех случаях, когда всю эпюру давлений, действующую на поверхность, можно представить одной аппроксимирующей формулой. Однако возможно, что на поверхность воздействует несколько разделенных эпюр давлений (рис. 3б). В этих случаях требуется построение нескольких аппроксимирующих формул типа (14). Или же можно использовать двухмерные полиномиальные сплайн-аппроксимации [12]. Причем сплайн строится из участков, соответствующих граням КЭ упругой задачи. При этом обеспечивается необходимая гладкость общей аппроксимирующей формулы и получаются локальные аппроксимации для отдельных граней.

Вычисление перемещений поверхностей по нормали (дополнительных зазоров)

В результате решения упругой задачи имеем перемещения узлов упругих КЭ (по три составляющие в соответствии с декар-

товой системой координат). По этой информации требуется получить для узлов треугольной сетки ГД-задачи значения до-

стей). Перемещения любой точки грани КЭ можно вычислить по формулам

полнительных толщин масляного слоя (дополнительных зазоров от упругого деформирования одной и другой поверхно-

8 8 8

их = 2 N (5,1)их1; иу = 2 Н, (5,1)иу1; { = 2 Н, (5,1)и..

2 = 1 2=1 2=1

Далее, имея три компоненты перемещений, можно вычислить перемещение по нормали в этой точке поверхности. Для

и = и соб(п, е) + и СОБ(П, е) + и СОБ(П, е3) где направляющие косинусы вычисляются по формулам (11). Именно эти значения, вычисленные для точек, соответствующих узлам сетки КЭ гидродинамической зада-

(15)

этого используются направляющие косинусы.

" " (16)

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

Выводы

Рассмотрены алгоритмы снижения трудоемкости решения связанной УГД-задачи. Показано, что можно существенно снизить трудоемкость за счет следующих мероприятий: сгущения сеток и применения метода редуцирования, применения квазисовместных аппроксимаций, процедуры конденсации. Причем процедуру конденсации можно реализовать в экономичном алгоритме как неполную процедуру прямого хода в методе исключения Гаусса. В целом трудоемкость решения УГД-задачи может существенно снизиться.

СПИСОК ЛИТЕРАТУРЫ

1. Мишин, А.В. Конечноэлементная методика расчета динамически нагруженных подшипников скольжения с учетом отклонений формы рабочей поверхности от цилиндрической / А.В.Мишин, М.В.Зернин // Сборка в машиностроении и приборостроении. - 2008. - № 2. -С. 43-54.

2. Зернин, М.В. Методика расчетной оценки предельных размеров дефектов поверхностей подшипников скольжения по критерию влияния их на параметры гидродинамики / М.В.Зернин, А.В.Мишин, Н.Н.Рыбкин // Вестник Брянского государственного технического университета. - 2013. - № 3. - С.14-23.

3. Рыбкин, Н.Н. Реализация методики расчетной оценки параметров гидродинамики подшипников скольжения с учетом радиальной податливости поверхностей / Н.Н.Рыбкин, М.В.Зернин// Вестник Брянского государственного технического университета. - 2013. -№ 4. - С. 59-65.

4. Программа для ЭВМ «Bearing Builder Finite Element Method»/ А.В.Мишин, М.В.Зернин,

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

Н.Н.Рыбкин, С.М.Шалыго, В.В.Сопранцов. -Зарегистр. в Федер. службе по интел. собств. (Роспатент) 09.01.14, Свид. № 2014610341.

5. Зернин, М.В. Гидродинамический анализ подшипников скольжения. Ч. 1. Учет неци-линдричности рабочих поверхностей/ М.В.Зернин, А.В. Мишин, Н.Н.Рыбкин, С.В. Шилько //Трение и износ. - 2014. - Т. 35.- № 5.- С. 584-595.

6. Шабров, Н.Н. Метод конечных элементов в расчётах деталей тепловых двигателей/ Н.Н.Шабров. - Л.: Машиностроение, 1983. -212 с.

7. Зенкевич, О. Конечные элементы и аппроксимация: [пер. с англ.] / О.Зенкевич, К.Морган. -М.: Мир, 1986.-318 с.

8. Зернин, М.В. Трехмерный конечный элемент для моделирования температурной, упругой и термоупругой составляющих в связанной задаче термоупругогидродинамики узлов трения / М.В.Зернин, Н.Н.Рыбкин // Вестник Брянского государственного технического университета . - 2016. - № 5. - С. 23-32.

9. Вороненок,Е.А. Метод редуцированных элементов для расчета конструкций / Е.А.Вороненок, О.М.Палий, С.В.Сочинский. -М.: Судостроение, 1990. - 224 с.

10. Рудаков, К.Н. ив8 Бешар 9.3. Геометрическое и конечноэлементное моделирование конструкций/ К.Н.Рудаков. - М.: ДМК Пресс, 2009. - 296 с.

11. Гордон, Р. Аппроксимация эмпирически полученной поверхности методом наименьших квадратов / Р.Гордон. - Режим доступа: http://elib.ispu.ru (дата обращения: 20.02.2012).

12. Корнишин, М.С. Вычислительная геометрия в задачах механики оболочек/ М.С.Корнишин, В.Н.Паймушин, В.Ф.Снигирев. - М.: Наука, 1989. - 208 с.

1. Mishin, A.V. Finite element procedure for computation of friction bearings loaded dynamically taking into account deviations of operation surface form from cylindrical one / A.V.Mishin, M.V.Zernin // Assemblage in Mechanical Engineering and Instrument Making. - 2008. - № 2. -pp. 43-54.

2. Zernin, M.V. Method for rated computation of surface defect limit dimensions for friction bearings on criterion of their impact upon hydrodynamics parameters / M.V.Zernin, A.V.Mishin, N.N.Rybkin // Bulletin of Bryansk State Technical University. - 2013. - № 3. - pp. 14-23.

3. Rybkin, N.N. Realization of procedure for estimated parameters of friction bearing hydrodynamics taking into account surface radial flexibility / N.N.Rybkin, M.V.Zernin// Bulletin of Bryansk State Technical University. - 2013. - № 4. -pp. 59-65.

4. Computer Program "Bearing Builder Finite Element Method"/ A.V.Mishin, M.V.Zernin, N.N.Rybkin, S.M.Shalygo, V.V.Soprantsev. -Registered in Federal Service of Intellectual Property. (Rospatent) 09.01.14, Certificate № 2014610341.

5. Zernin, M.V. Hydrodynamic analysis of friction bearings. Part 1. Registration of operation surface non-cyclicity/ M.V.Zernin, A.V. Mishin, N.N.Rybkin, S.V. Shilko //Friction and Wear. -2014. - Vol. 35.- № 5.- pp. 584-595.

Сведения об авторах:

Зернин Михаил Викторович, к.т.н., доцент кафедры «Информатика и программное обеспечение» Брянского государственного технического университета, e-mail: zerninmv@mail.ru.

Zernin Mikhail Victorovich, Can. Eng., Assistant Prof. of the Dep. "Informatics and Software", Bryansk State Technical University, e-mail: zerninmv@mail.ru.

6. Shabrov, N.N. Finite Element Method in Computation of Heat-Engine Parts/ N.N.Shabrov. - L.: Mechanical Engineering, 1983. - pp. 212.

7. Zenkevich, O. Finite Elements and Approximation: [transl. from Engl.] / O.Zenkevich, K.Morgan.- M.: Mir, 1986.-pp. 318.

8. Zernin, M.V. 3 D finite element for simulation of temperature resilient and thermo-elastic constituents in coupled problem of friction unit thermodynamics / M.V.Zernin, N.N.Rybkin // Bulletin of Bryansk State Technical University . - 2016. - № 5. - pp. 23-32.

9. Voronenok,E.A. Reduced Element Method for Structure Computation / E.A.Voronenok, O.M.Paly, S..V.Sochinsky. - M.: Shipbuilding, 1990. - pp. 224.

10. Rudakov, K.N. UGS Femap 9.3. Geometrical and Finite Element Simulation of Structures/ K.N.Rudakov. - M.: DMK Press, 2009. - pp. 296.

11. Gordon, R. Approximation of surface obtained empirically by least-squares method / R.Gordon. - Access mode: http://elib.ispu.ru (Address date: 20.02.2012).

12. Kornishin, M.S. Computational Geometry in Problems of Casing Mechanics/ M.S.Kornishin, V.N.Paimushin, V.F.Snigiryov. - M.: Science, 1989. - pp. 208.

Статья поступила в редколлегию 9.03.17. Рецензент: д.т.н., профессор Брянского государственного технического университета

Кеглин Б..Г.

Рыбкин Николай Николаевич, аспирант кафедры «Механика и ДПМ» Брянского государственного технического университета, e-mail:

rnikolai@mks.ru.

Rybkin Nikolay Nikolayevich, Post graduate student of the Dep. "Mechanics and DPM", Bryansk State Technical University, e-mail: rnikolai@mks.ru.