Решение задач с односторонними связями и трением Кулона при динамических воздействиях Текст научной статьи по специальности «Физика»

Динамика конструкций и сооружений

РЕШЕНИЕ ЗАДАЧ С ОДНОСТОРОННИМИ СВЯЗЯМИ И ТРЕНИЕМ КУЛОНА ПРИ ДИНАМИЧЕСКИХ ВОЗДЕЙСТВИЯХ

А. А. ЛУКАШЕВИЧ, канд. техн. наук, доцент

Санкт-Петербургский государственный политехнический университет 195251, Санкт-Петербург, ул. Политехническая 29, ПГК, а.105, а. 1^@,Ьк ги

Предлагается численный подход к решению динамических задач с односторонними связями и трением Кулона. Решение строится на основе метода пошагового анализа (по времени) с использованием рамных контактных конечных элементов. Приведенный алгоритм иллюстрируется численным примером.

КЛЮЧЕВЫЕ СЛОВА: контактная задача, односторонняя связь, трение Кулона, пошаговый анализ.

Задачи с односторонними связями и трением Кулона между контактирующими поверхностями довольно часто встречаются при расчете различного рода сооружений. Решение таких задач при действии статических нагрузок рассматривалось, например, в работах [1, 2]. В то же время нередки случаи, когда требуется учитывать динамические воздействия на сооружение. Конструктивная нелинейность здесь будет проявляться в смене рабочих схем сооружения по времени - включение и выключение односторонних связей, как в нормальном, так и касательном направлении. Предполагается, что между двумя последовательными событиями на контакте, т. е. в пределах каждой такой рабочей схемы, характер деформирования системы является линейным.

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

Рассмотрим плоскую динамическую задачу взаимодействия линейно упругих тел V+ и V ~ (которыми могут являться, например, сооружение и основание) с контактирующими поверхностями соответственно и S~ . Для расчета указанной системы будем использовать дискретную расчетную модель МКЭ, для которой справедливо следующее матричное уравнение движения с начальными условиями [3]:

ми *+м + Си *+м + ки = Р ; и| = и 0, и| = и 0, (1)

' 11=0 0' 11=0 0' у '

где М, С и К - соответственно матрицы масс, демпфирования и жесткости системы конечных элементов; и , и , и и Р '+Л' - векторы узловых перемещений, скоростей, ускорений и внешней узловой нагрузки в момент времени t +Дt. Кроме того, на части внешних границ , (соответственно для

V +, V ~) должны быть заданы граничные условия в усилиях, а на , «и - в перемещениях.

Принимая во внимание, что перемещение в момент времени t +At можно представить как Ut+At = Ut + AUt+At, где AU t+At - приращение перемещений на шаге At, преобразуем уравнение (1) к виду, позволяющему решение конструктивно нелинейной динамической задачи свести к решению последовательности линейных динамических задач на основе метода пошагового анализа состояния на контакте. На каждом шаге по времени At будет решаться система линейных дифференциальных уравнений второго порядка вида:

MÜt+At + CÜt+At + KAUt+At = Pt+At - KUt. (2)

Здесь предполагается, что в начальный момент времени (при t = 0) известны векторы перемещений, скоростей и ускорений и необходимо найти решение (2) на интересующем нас интервале времени от 0 до некоторого T.

Для численного интегрирования уравнений движения (2) используется неявная конечно-разностная схема Ньюмарка, в основе которой предположение о линейном изменении ускорений в интервале At. При этом используются следующие зависимости между приращениями перемещений, скоростями и ускорениями для момента времени t +At:

U*+a* =-[au *+a* - AtU 11-f— -11U t;

«(At)2 L J 12« ) (3)

U t+At = U t + [(1 - ß)U t + ß U *+a* ] At.

Здесь а и ß - параметры, определяющие точность и устойчивость процесса интегрирования. Принимаем а = 1/4, ß = 1/2, что соответствует случаю постоянного среднего ускорения на каждом из интервалов At. В этом случае схема интегрирования по Ньюмарку для линейных задач является безусловно устойчивой, т. е. решение при любых начальных условиях не растет неограниченно при больших значениях шага At [3].

Подставив выражения (3) в уравнение (2), тем самым, исключив U t+At и U t+At из числа неизвестных, после несложных преобразований получим следующее матричное уравнение для определения AUt+At:

K AUt+At = P*+At, (4)

где

K = K + a0M + «C;P*+a* = P t+At + (a2Ut + a3Ut) + C(a4Ut + a5Ut)- KUt.

Коэффициенты а0 - а5 зависят от шага At и параметров а и ß:

1 ß 1 1 1 ß л ßAt .

a0 =-al =-; a2 =-; a3 =--1; a4 =--1; a5 =--At.

a(At )2 aAt aAt 2a a 2a

Система алгебраических уравнений (4) решается методом LDLT - факторизации, учитывающим разреженность симметричной матрицы K и ее переменный профиль. После нахождения AUt+A и, соответственно, Ut+A для вычис-

" TT t + At « TT t + At n

ления ускорений U и скоростей U используются уравнения (3). В свою очередь в любой момент времени t' в пределах интервала At (t < t' < t + At) значения ускорений U(t'), скоростей U(t' ) и перемещений U(t') могут быть рассчитаны по следующим формулам:

U(t') = Ut + [U t+At - U11; U(t') = Ut + [Ut + U '];

W At 2 (5)

U(t') = Ut + (t '-1) Ut + (t '-1)2 [Üt + U' 1/4.

Первое из уравнений (5), согласно схеме Ньюмарка, выражает линейный закон изменения ускорения на интервале At, второе и третье получаются из выражений (3) при подстановке в них значений а = 1/4, в = 1/2 и замены величины (г +Аг) на величину t'.

Кроме того, при учете односторонних связей с трением Кулона дополнительно к начальным условиям при t = 0 и граничным условиям на , Би должны выполняться также следующие условия на граничных поверхностях Б+, Б-, записанные для момента времени t [2]:

иП +5^ > 0; аП < 0; (иП ) а 1п = 0, хе Б+, Б;; (6)

<3т

<-/а п;

а,

п

г ■ г

< 0; (

а,

+ /а п) и, = 0 х е «с

.

(7)

Здесь ип, и, - взаимные перемещения контактных поверхностей соответственно по нормали и по касательной в момент времени г; и, = ди, / дг - скорость взаимных касательных перемещений; 5п - заданный начальный зазор между

контактными поверхностями; агп, а, - нормальное и касательное напряжения на площадке контакте; / - коэффициент трения Кулона.

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

данного слоя). Полагая, что контактный слой прикреплен к поверхности Б- и его толщина £ мала по сравнению с длиной контакта, представим расчетную схему контактного взаимодействия тел V + и V-, как показано на рис. 1, а.

б)

и IN

их 7 б> Я 7 А

V

т

▼ и

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

Считая 7-й ККЭ (рис. 1, б) как 7-ю дискретную связь (по нормали и по касательной) между контактными узлами конечно-элементной сетки, расположенными на и Б-, перепишем условия (6), (7) в терминах усилий и перемещений для каждого 7-го контактного элемента (7 е Бс, = Б+ и Б-):

и^ + 5° > 0; Ы; < 0; (и^ + 5° ) Ы* = 0; Ы < |е<| = -/Ы/; д'и^ > 0; (&'-е') и^ = 0.

(8) (9)

Здесь и'п7, и'т7 - взаимные нормальное и касательное смещения противолежащих узлов 7-го ККЭ на и Б-; 1л^7 - скорость взаимного смещения (проскальзыва-

ния) этих узлов, Ы/ , Q■ - продольная и поперечная силы в 7-м ККЭ в момент времени t, причем поперечная сила по направлению противоположна контактному напряжению С; Qtpi - предельная поперечная сила, соответствующая

предельному касательному напряжению, но также противоположная ему по направлению:

Qtpi = ^П ) • ^

Р7

* , С 0 N

Последнее из условий (8) означает, что при контакте (ип + 5п) = 0, N ■ < 0 ; при отрыве Ы- = 0, и1П1 + 5 > 0. Два последних условия (9) имеют следующий физический смысл - они устанавливают соответствие между скоростью взаимного проскальзывания противолежащих узлов на S+ и Sc" и вели-

Qгt < QU

чиной усилия Q| в момент времени При условиях

'■1*7 = 0 имеет

место состояние сцепления (допредельного трения); при Q■ = Qtpi, ^ 0 -

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

случае направления u!ti и Qit должны совпадать).

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

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

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

Считается, что в текущий момент времени * состояние на контакте известно. Для каждой 7-й связи (т. е. 7-го ККЭ) определены значения взаимных смещений и1П7, и , скоростей и и контактных усилий N1, Q¡ , Qtpi . Пусть часть связей (7 е £1с) находится в состоянии сцепления, другая часть (7 е Б2С) - в состоянии контакта со скольжением и, наконец, третья часть связей (7 е 8-$с) - в состоянии отрыва. При этом 8с = £1с и $2с и 3-$с . В начале расчета, при * = 0, перемещения иП 7 и и ° , (также, как скорости и ускорения) принимаются нуле-

выми. Задаются величины , которые, в частности, могут быть равны нулю.

1. Выполняется основной шаг по времени At (в процессе расчета величина

текущего шага может изменяться в соответствии с установленными моментами

наступления событий на контакте). Из решения (4) определяются приращения

- А г+Аг л г+Аг г+Аг г+Аг

перемещений Аип- , Аит- , затем перемещения ип- , ит- , скорости

и{+Аг и усилия N1+Аг, Qti +Аг, Аг для момента времени г +Аг.

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

Таким образом, момент проскальзывания (при срабатывании условия

Qi - Qp у) для каждой г-й связи, находящейся до этого в состоянии сцепления,

вычисляем с помощью следующей итерационной формулы:

Д-1Ч /Д-К ^

у е ^. (10)

т к ,к-1 ¡-.к-1 ,к-1\ - гИ + (г2у - гИ )

QPi щГ) - Qi (С1)

[Qi (г!-1) - (Ть-1)] - [Qpi (г|г-1) - Qpi (гь"1)]

Здесь г1к 1, 1 - начало и конец интервала, в котором с помощью линейной интерполяции находится к-го приближения Тгк . Для первого приближения этот интервал (г-, г°) равен шагу Дг, т. е. г°. - г, г° - г + Аг, Qг (Т-0) - Qгt,

Qг (г°) - Qi +АТ, в дальнейшем, при последующих приближениях, интервал монотонно сужается. Так, например, после нахождения очередного приближения Тгк , с помощью уравнений (5) выполняется пересчет взаимных смещений, а по ним и усилий Qi (Тк), Qpi (Тгк) для момента времени Тк и сопоставление полученных величин с граничными значениями Qi (Т^_1), Qpi (Т1ки Qi (Т^1), Qpi (Г2к_1) для текущего интервала. В результате чего устанавливаются границы

интервала (^, г2,г) для следующего - (к+1)-го приближения. Окончанием итерационного процесса служит достижение относительной разности величин Qi (Т-1с ) и Qpi (Ггк) некоторой заданной точности вычислений -

-рг

Qr (ТГ ) - QVг (Тк)

< е . При этом окончательно принимаем Т- - Т]к .

• г

Qpг (1к )

Аналогично определяем момент зацепления (при срабатывании условия - - 0) для каждой г-й связи, находящейся до этого в состоянии скольжения:

г е ^. (11)

тк - Л-1 + (гк-1 Л-1)

к 1

-Щг(гц )

к 1 к 1 и%г(г2г )-и%г(г1г )

Момент отрыва, т. е. выключения г-й односторонней связи (при срабатыва-

нии условия N! _ 0), находящейся до этого в состоянии контакта:

7k _ tk-1 + (tk-1 tk-l) ti _ t1i + (t2i -t1i )

-N (ti"*)

Ni(t2k-')- Ni(t*-1)

i e Sic, S2c. (12)

Момент контакта, т. е. включения 7-й односторонней связи (при срабатывании условия и1П1 +5 П7 = 0 ), находящейся до этого в состоянии отрыва:

7 k _ + k-1 . k-i v k-i \ ti _ t1i + (t 2 i - t1i )

( i+k-К e 0

-Um (t1i )- 5m

v Um (t2k-1) - Uni (t1k,

i e S3C. (13)

3. Из всех значений ti, найденных с помощью формул (10) - (13) и находящихся в интервале (t, t +At), выбирается наименьшее, соответствующее моменту наступления ближайшего по времени события на контакте:

7 _ min (7), i e Sc. (14)

В случае, если t > t + At - выполняется следующий основной шаг интегрирования по времени продолжительностью At и уже в пределах этого шага производится поиск момента наступления очередного события на контакте (т. е. следует переход к п. 1).

4. Если же t < t < t + At - выполняется перерасчет уточненного таким образом шага с величиной A t _ t - t. При этом изменение предельных сил трения на контакте учитывается приложением к противолежащим узлам i-го ККЭ на S +c, S-c корректирующих сил FT { (с соответствующим знаком), принимаемых:

F i _-AQP1 _-A • (QT-Q'pi), i e S 2 с. (15)

На ту же величину корректируется значение поперечной силы в i-м ККЭ:

Qt _ Qt +AQP1, i e S2с. (16)

В результате перерасчета шага определяются значения u tTi, u tTi, N ti , Q't , Q'pj для текущего момента времени t . Проверяются условия наступления ожидаемого события - в случае проскальзывания это будет условие Q i _ Qtpi, зацепления - U^ _ 0, отрыва - N. _ 0 , контакта - (utni + 5° ) _ 0 . Если соответствующее условие (с учетом заданной погрешности вычислений) не срабатывает, следует еще раз уточнить момент времени t , но уже на интервале (t, 7 ), либо (7, t+At).

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

Продемонстрируем изложенный подход на примере расчета плоской стержневой системы, которая может моделировать, к примеру, участок трубопровода с перепадом рельефа, в условиях динамического нагружения (рис. 2, а). Считается, что система закреплена от боковых смещений (т. е. из плоскости xy) и не разрушается в местах концентрации напряжений. На опорах i = 1, 2 дейст-

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

Изначально сооружение находится в состоянии покоя, затем на левом конце прикладывается горизонтальная импульсная нагрузка Р(г). Закон изменения импульса имеет треугольную форму, продолжительность 0,1 с, амплитуда 100 кН (рис. 2, б).

а) б)

Л

р (т)

м

777

<г

Л1

20 м

->

10 м

->

20 м

100

0 0,05 0,1 Т, с

2

х

Рис. 2. Стержневая система под действием импульсной нагрузки Р (Т)

Продольная жесткость стержней ЕА задавалась равной 9-106 кН, изгибная жесткость Е1 = 2-106 кН м ; погонная масса стержней т = 0,4 т/м. Начальные зазоры на контактах задавались нулевыми, т. е. 5П1 - 5П2 - 0 . Матрица демпфирования системы здесь вычислялась как С = а М (коэффициент а принимался равным 0,1 с-1). Контактное взаимодействие моделировалось с помощью рамных ККЭ, соединяющих опорные узлы стержневой системы с неподвижными опорами.

а)

«т,

А

0,5 -0,4- — 0,3---0,2---0,1---0

-0,1----0,2 ---0,3----0,4 ---0,5---

1 1 1 г1

1 'Л1/1

__J____1___ХУ___I

мм 2,5 --

— ---/ = 0,1

2,0 - -- - / = 0,5

— ----- / = 1,0

1,5 - — -----Г = 1,3

1,0-----Г---т—4—

0,5 --

0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 0,45 0,5 0,55 0,6 0,65 0,7 Т, с

0

Рис. 3. Смещения на опоре 1 в зависимости от времени а - горизонтальные; б - вертикальные (отрыв)

В начальный момент времени t = 0 на всех опорах задавалось состояние сцепления. Нормальные силы взаимодействия № на контакте узлов с соответствующими опорами принимались равными реакциям от собственного веса сооружения. В дальнейшем, в результате действия динамической нагрузки, возможно как проскальзывание, так и последующее зацепление на контакте, а также выключение (отрыв) - включение односторонних связей в течение рассматриваемого промежутка времени. При этом, вследствие геометрии рассматриваемой стержневой системы, нормальные силы взаимодействия Ы- также меняются во времени.

Численные расчеты производились с использованием разработанной автором вычислительной программы [5]. Целью расчетов являлась оценка влияния коэффициента трения (величина f в расчетах варьировалась от 0 до 1,5) на поведение стержневой системы при динамическом воздействии. На рис. 3 показаны горизонтальные и вертикальные смещения на опорах 1, 2 в зависимости от времени при разных значениях коэффициента f (шаг интегрирования здесь принимался А* = 0,0002 с).

На рис. 4 представлены зависимости между амплитудными значениями горизонтальных (скольжения) их и вертикальных (отрыва) ип смещений на опорах и заданной величиной коэффициента трения f на контакте. ит /ип, мм

/V

—1 — 1 1 1 1

1 1 1 1 —1— 1 1 1 1 1 1 1 1 1 1 на опоре _ 1

1 1 1 1 1 ч 1 1 1 1 1--- 1

1 л 1 1 1 1- _1 \ 1 1 1 1 1 1 Ь 1 _ 1 1 Л 1 1 1 1- 1 л 1 1 1

1 1 1 1 —\— 1 1 1 1 1 1 —1----1--- V V 1 т 1 —+- 1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 ит 1 1 п — 4— 1 1 1 г~

- т _ Г 1 п т 1 —

1 --1— 1 —1— 1 —1- 1 1 —1-1— 1 1—1- — 1--1— —1-1— —!—1-—1-1- ~ Т —• -О. —1-1-И»

0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1 1,2 1,3 1,4

Рис. 4. Амплитудные значения смещений ит и ип на опорах для разных значений коэффициента трения /

Как видно из приведенных графиков, с увеличением коэффициента f отрыв сооружения от опор в результате динамического воздействия заметно уменьшается. Причем в рассмотренном здесь примере имеет место некоторое пороговое значение коэффициента f (0,6-0,7), при котором влияние трения на контакте на величину максимального отрыва имеет экстремальный характер.

Для того, чтобы изучить зависимость решения от величины шага интегрирования по времени, было просчитано поведение стержневой системы для разных значений А* в диапазоне от 0,0001 до 0,0064 с (с последовательным удвоением длины шага). Сравнение полученных результатов позволяет заключить, что предлагаемый численный подход показывает удовлетворительную внутреннюю сходимость в достаточно широком диапазоне временных шагов интегрирования. Так, величины горизонтальных и вертикальных смещений на опорах не намного отличаются при назначении основного шага в диапазоне от 0,0001 до 0,0008 с. С дальнейшим увеличением А* наблюдается некоторое ухудшение сходимости, особенно при больших значениях коэффициента трения. Заметим, что в общем случае выбор оптимального шага интегрирования представляет собой достаточно сложную проблему [3]. В этой связи, на основании проведен-

73

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

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

1. Лукашевич А. А. Решение контактных задач при учете податливости и прочности односторонних связей методом пошагового нагружения // Вестник гражданских инженеров. - 2009. - № 3 (20). - С. 18-23.

2. Вовкушевский А. В. Постановка и решение контактной задачи теории упругости с трением при произвольном процессе нагружения // Труды ЛПИ. - 1985. - № 405. -

3. Бате Л., Вилсон. Е. Численные методы анализа и метод конечных элементов. -М.: Стройиздат, 1982. - 448 с.

4. Лукашевич А. А. Построение и реализация схем прямого метода конечных элементов для решения контактных задач // Изв. вузов. Строительство. - 2007. - № 12. -

5. Лукашевич А. А. Решение двумерных задач контактного взаимодействия с учетом односторонних связей и трения Кулона методом конечных элементов / Программа для ЭВМ. Зарегистрирована в Реестре программ для ЭВМ Федеральной службы по интеллектуальной собственности, патентам и товарным знакам РФ № 2008610113, 2008.

THE DECISION OF PROBLEMS WITH UNILATERAL CONSTRAINTS AND CULON'S FRICTION UNDER THE DYNAMIC LOAD

Lukashevich A. A.

The numerical approach to the decision of dynamic problems with unilateral constraints and Culon's friction is offered. The decision is formulated on the basis of step-by-step analysis (in time) method with use of frame contact finite elements. The reduced algorithm is illustrated by numerical example.

KEY WORDS: contact problem, unilateral constraint, Colon's friction, step-by-step

Л и т е р а т у р а

С. 9 - 13.

С. 18 - 23.

analysis.