CC BY
32
УДК 519.22
К. А. Беклемишева, И. Б. Пет,ров, A.B. Фаворская
Московский физико-технический институт (государственный университет)
Численное моделирование процессов в твердых деформируемых средах при наличии динамических контактов с помощью сеточно-характеристического
метода
Рассматривается ряд задач динамического деформирования при наличии подвижных контактирующих поверхностей. Адаптирован сеточно-характеристический метод, использующий нерегулярные треугольные двухмерные сетки и гибридную схему 1-2 порядка аппроксимации. Предлагается алгоритм расчета сил трения скольжения и покоя. Получено численное решение задач качения колеса и генерации сдвиговых волн.
Ключевые слова: механика деформируемого твердого тела, сеточно-характеристический метод, нерегулярные сетки, триангуляция Делоне, динамический контакт, обнаружение столкновения, сила трения, сейсморазведка, метод падающего груза.
Значительная группа нестационарных процессов, происходящих в твердых деформируемых телах при наличии контактного взаимодействия двух или более тел, может быть промоделирована с помощью полной замкнутой динамической системы уравнений механики сплошных сред (\!('(') с корректным решением задачи контактного разрыва (см., например, [1, 2]). Данная система дифференциальных уравнений, включающая уравнения движения и реологические соотношения, связывающие тензоры напряжений и деформаций (или их временные производные), является гиперболической. Корректное решение задачи динамического контактного разрыва представляет собой самостоятельную проблему при использовании численных методов для решения задач МСС и может быть реализовано только при использовании подходов, учитывающих характеристические свойства определяющей системы уравнений. Как хорошо известно из курсов уравнений математической физики и численных методов (см., например, [3-5]), такими свойствами обладают сеточно-характеристические методы [5-7], используемые в данной работе. Пример корректного численного решения задачи фиксированного контактного разрыва для твердой деформируемой среды предоставлен в работе [8].
Для оценки возможных повреждений и получения информации о внутренней структуре объекта необходимо иметь возможность корректно описывать волновые процессы как внутри компонентов системы, так и на многочисленных свободных границах и контактных поверхностях.
Моделирование динамического контакта при решении различных задач механики деформируемого твердого тела (МДТТ) представляет собой довольно сложную задачу. При ее решении необходимо реализовать поиск контактирующих узлов расчетной сетки, не допускать пересечения сеток («проникания» одного тела в другое), определять, когда необходимо изменять тип условий на контактной границе.
Также в рамках сеточно-характеристического метода возможно исследовать эффекты, связанные с расчетом силой трения, возникающей на динамическом контакте. Разработка алгоритма ее расчета, приведенного в работе, открывает возможность для решения многих прикладных задач, связанных с проблемой динамического контакта. В работе рассмотрены следующие задачи:
1) Генерация сдвиговых волн для сейсморазведки. Сдвиговые волны, которые в последнее время находят широкое практическое применение в исследовании земной коры и
залежей нефти, могут быть сгенерированы различными способами, но особый интерес представляет метод падающего груза. Этот метод эффективен по времени и трудозатратам, но пока что мало кем реализован по причине недостаточности теоретического исследования проблемы.
2) Качение колеса, которое представляет собой актуальную для железных дорог задачу. Моделирование данного процесса дает возможность исследования на прочность колес поездов и рельсов, влияния сколов и неоднородностей на их разрушение.
Для математического моделирования волновых процессов в исследуемых объектах использовалась система динамических уравнений в частных производных теории упругости (см., например, [1, 2]). Данная система является гиперболической и распадается на совокупность уравнений переноса, что подробно описано в [3-5].
В данной работе используется гибридная разностная схема 1-2 порядка аппроксимации, предложенная в [7]. Схема является устойчивой при выполнении условия Куранта и адекватно описывает динамические процессы. Опорными для нее являются схемы Куранта-Изаксона-Риса и Лакса-Вендроффа. Теория построения различных сеточно-характеристических методов изложена в [5].
В расчетах использовались нерегулярные (треугольные) сетки, генерируемые при помощи триангуляции Делоне. Соответственно в данной работе реализован метод, позволивший производить расчеты на подобных сетках с использованием сеточно-характеристического подхода [9, 10].
Данные сетки имеют целый ряд преимуществ. Во-первых, имеется большое количество программных инструментов промышленного качества для построения треугольных сеток в областях произвольной формы и связности. Во-вторых, треугольная сетка позволяет полностью управлять размерами ячеек, добиваясь либо полной однородности, либо при необходимости сгущения или разряжения в определенных зонах. В-третьих, в работе [11] показано, что использование гибридных схем с расщеплением на неструктурированных треугольных сетках вполне оправдано, так как позволяет получать численные решения удовлетворительного качества в областях интегрирования произвольной формы.
Корректное решение задачи в точках, принадлежащих границам области интегрирования, требует постановки граничных условий в количестве, равном числу выходящих за эту область характеристик [3]. Формулы для общего вида линейных граничных условий, использованные для данной задачи, приведены в [6, 12].
Для расчета искомых функций в каждом узле, лежащем на контактной границе двух деформируемых тел, необходимо знать скорости близлежащих узлов, чтобы для их анализа определить максимально допустимый временной шаг, при котором тела еще не пересекутся. Кроме того, необходимо найти все достаточно близко расположенные к исследуемому узлу участки границы и произвести дальнейшие проверки касательно наличия контакта в этом районе. Определение ближайших границ в нашем случае производится на основе специально построенной триангуляции с вершинами в граничных узлах. Вершинами в данной триангуляции выступают все граничные узлы тел; для построения триангуляции можно использовать те же алгоритмы, что и для сетки. Также можно лишь проводить повторную оптимизацию триангуляции после смещения (удаления, добавления) граничных узлов тел, произошедшего в результате перехода к следующему временному слою.
Ограничение шага интегрирования гарантирует, что никакой треугольник не выродится, а, значит, границы тел не пересекутся. Для определения участка границы, с которым взаимодействует некоторый узел, выполняется перебор всех треугольников, смежных с узлом, и рассматриваются их ребра, лежащие напротив этого узла. Время поиска контакта для каждого узла пропорционально количеству инцидентных ему ребер в триангуляции. А поскольку каждое ребро инцидентно лишь двум узлам, то время обработки всех граничных узлов всех тел составляет О(п).
Рассмотрим более подробно алгоритм расчета искомых функций в точках, принадлежащих контактирующим поверхностям. Если оказалось, что расстояние между двумя узлами, принадлежащими разным участкам границы одного или двух тел, меньше, чем заранее определенный порог, то предстоит установить, взаимодействуют ли эти участки границ друг с другом либо происходит потеря контакта, то есть появление свободных границ. При расчете контакта двух узлов, прежде всего, проводится пробный расчет в предположении, что узлы находятся в состоянии полного слипания; при этом вычисляются «пробные» значения для тензора напряжений и вектора скоростей на временном слое п + 1 1 и г>™+1. Из полученных при таком расчете данных мы рассчитываем силы, действующие на контакте.
Прежде всего рассмотрим силу, направленную по нормали к поверхности:
Р _ грП+1
F ПЛ. - г г.
Р.
Если она положительна, на узел тела а действует сила, сонаправленная нормали р. Узлы считаются неконтактирующими, то есть рассчитываются как две свободные границы. Если же сила отрицательна, то она работает как сила реакции опоры N, и возникает сила трения.
Рассмотрим силу, направленную перпендикулярно к нормали:
F™ _ ТП+1 •
Р1.
Условие на эту силу позволяет разграничить силы трения покоя и скольжения. Если рта < кЕпа, мы получаем случай трения покоя: внешняя сила не превышает силы сцепления. В этом случае становится верным предположение о том, что контактирующие узлы находятся в состоянии полного слипания, и можно использовать уже рассчитанные результаты.
В обратном случае, еслиРта > кЕпа, производится новый расчет по формулам для трения скольжения. Вывод формул приводится ниже.
Найдем силу трения скольжения из следующих условий:
г = = /,
/т = к fp,
Уар = уьр , (V£ - V*) Гг-
Граничное условие при заданной внешней силе:
vn+1 _ vin - — Р + -(1 - — I (Р ■ п) п,
( — - — ) (Р ■ п)1 \С2 ClJ
РС2 Р \ С2
Tn+1 _Tin ± (z ® п + п ® Р) ± J-HL (XI - 2 (л + N00),
X + 2ц
где z — тТтп - /, верхний знак относится к левой границе области интегрирования, нижний - к правой.
Такие смешанные ограничения при использовании такого вида задания граничных условий можно привести к граничным условиям с заданной внешней силой. Получается, сила, действующая на тело а:
Iе1 _ - fb _ f _ fPP + (sign (fT) к fp) p1, где
U _ p~y № - ) ■ P) + 2^ (CC - Tin) f ■ P +1 (| -1) mn - Т™) p) ■ f .
1. Качение железнодорожного колеса
Высокая плотность загрузки железнодорожных путей ставит все более высокие требования к скорости и точности диагностики, и старые методы становятся неэффективны из-за низкой скорости перемещения и обработки данных мобильных датчиков. В случае же запуска по российским железным дорогам высокоскоростных экспрессов (до 500 км/ч) потребуется дополнительный контроль состояния рельсов и колес. Поэтому разработка высокоскоростных мобильных датчиков, основанных на улавливании отклика акустического сигнала, генерируемого при движении поезда, является на данный момент важной задачей. Кроме того, анализ волновой картины и характерных зон постоянного высокого напряжения при численном решении может дать оценку причин возникновения трещин и позволит оптимизировать форму колеса и рельсов.
Для численного решения динамической двухмерной задачи о динамическом деформировании колеса, катящегося по рельсу, используется модель линейно-упругой сплошной среды [1] и применяется сеточно-характеристический гибридный метод на нерегулярных (треугольных) расчетных сетках.
Общий вид расчетной области предоставлен на рис. 1, характерный размер сетки h = 0.008 м, шаг по времени, определяемый из условия Куранта, г = 1.33-10-7 с.
На нижней границе рельса поставлено условие жесткого закрепления (скорость равна нулю). На внутренней границе колеса (контакт с осью) поставлено условие, отражающее вес вагона, что выражается в дополнительном нормальном «прижимающем» напряжении на нижней половине границы. Сила тяжести, действующая на отдельные элементы материала, не учитывается. Мы можем ей пренебречь, т.к. вес вагона гораздо больше, чем вес колеса.
Начальные условия - колесо движется с постоянной скоростью to — 10 м/с, расчет проводился при коэффициентах трения к\ = 0.01 = 0.8, а также при полном слипании на контакте (условие на тангенциальную силу не проверяется, всегда считается случай полного слипания).
1.1. Результаты расчета
На рис. 2 слева отображено поле скоростей: скалярное поле отображает модуль скорости, векторы - направление. На рисунке 2 справа отображено поле компоненты тензора напряжений sxx. На каждом из них а) к\ = 0.01, б) к2 = 0.8 и в) полное слипание.
Как следует из приведенных результатов расчетов, при большем коэффициенте трения колесо сильнее тормозится, но и сильнее закручивается. В случае полного слипания в течение нескольких десятков шагов по времени устанавливается качение без проскальзывания. При конечном коэффициенте трения качение без проскальзывания наступает в течение более долгого времени. Для качения колеса с коэффициентом трения к2 = 0.8 проскальзывание исчезает в момент времени t = 0.0288 с, что согласуется с теоретическими расчетами.
Также из этого рисунка следует, что интенсивность компоненты не зависящей, по
большому счету, от прижимного напряжения, увеличивается при увеличении силы трения.
Г,
9
*
Рис. 2. Поле скоростей (слева) в момент времени Ь = 0.0232 с и компоненты тензора напряжений зжж(справа) в момент времени Ь = 0.0033 с; а) к\ = 0.01, б) = 0.8, в) полное слипание
2. Генерация сдвиговых волн
Применение сдвиговых волн (б'-волн) для сейсморазведки нефтяных залежей является одним из наиболее активно развивающихся направлений. Разницу в свойствах распространения продольных и сдвиговых волн можно использовать для разнообразных практических задач: например, определение пористости пласта, изучение структуры литосферы или различение обломочных и известковых пород. Основную сложность при изучении и моделировании представляют анизотропия сдвиговых волн в слоистых средах и отсутствие эффективных источников генерации этих волн. За последние десять лет исследователи перешли от математических исследований, о которых написано в статье [14], к практическому применению, но еще требуется множество исследований, чтобы технология была полностью отработана.
На данный момент в основном применяются три метода генерации: метод взрывных зарядов, метод падающмх) груза и вибрационные сейсмические источники. В данной работе мы рассмотрим генерацию сдвиговых волн методом падающмх) груза. Этот метод достаточно экономичен энергетически, не требует столь до.;п'ой нодх'отовки, как взрывные траншеи и удобен при отсеве приповерхностных волн. Устройство для мх) осуществления и сравнение с аналогами приводится в [13].
Тяжелая (порядка двух тонн) металлическая плита со специальным покрытием, увеличивающим сцепление с поверхностью, закреплена на грузовике. Конструкция доставляется на исследуемую точку, после чего приводится в действие. Плита опускается на грунт, молот, закрепленный в раме, поднимается в верхнее положение. После этшх) молот свободно надает и ударяет но плите горизонтально в нижней точке свой траектории. Через несколько секунд, после прохождения всех сейсмических сигналов, рама с молотом разворачивается на 180°и производится аналогичный удар с противоположной стороны плиты. Таким образом, в породу уходят два сигнала противоположной полярности. Складывая их, мы можем получить акцентированные Р-волны с сократившимися 5-волнами, а вычитая - чистые 5-волны, не загрязненные возмущениями Р-волн. Метод сложения и вычитания сейсмо-
грамм является на данный момент достаточно типичным методом выделения б'-волн и их отражений из общей волновой картины, в частности, в методе взрывных зарядов. Также стоит заметить, что сенсор, установленный на самом устройстве, позволяет автоматически отфильтровывать приповерхностные волны: к тому времени, когда до него доходит отклик от низлежащих слоев, трещин или каверн, приповерхностные волны успевают отойти достаточно далеко, чтобы не вносить возмущения в волновую картину.
Исследование такого метода генерации сдвиговых волн имеет практический интерес, так как отечественные геологи пока не разработали достаточно эффективного аналога такого аппарата. Математическое моделирование может помочь определить параметры, которые наиболее сильно влияют на его технические характеристики.
2.1. Постановка задачи
Геометрия задачи изображена на рисунке 3. Для получения двух сигналов с помощью ударов с противоположных сторон ударник помещается симметрично относительно центра. Небольшому ударнику придается начальная скорость V = 20 м/с, на нижней границе металлической плиты стоит условие трения с коэффициентом к = 0.9, также задается прижимное воздействие на верхний край плиты Р = 107 Па. Характеристики массивной породы соответствуют граниту: р = 2700 кг/м3, с\ = 6000 м/с, с2 = 3500 м/с;
характеристики ударника и плиты соответствуют стали: р = 7800 кг/м3, с\ = 5700 м/с, с2 = 3100 м/с. Рассматривается отражение от нижней грани породы, на которой установлено условие свободной границы. Боковые границы удалены достаточно далеко, чтобы отражение от них не влияло на волновую картину около датчика, расположенного в точке, отмеченной рис. 3.
Рис. 3. Расчетная область для задачи с колесом
2.2. Результаты расчета
На рисунке 4 показано распределение поля скоростей в момент времени £ = 0.000404 с. Скалярное поле показывает модуль скорости. Видно, что волновая картина асимметрична и зависит от направления удара. Также видны сформировавшиеся Р-волна и ¿'-волна одинаковой интенсивности. Стоит заметить, что «окно» сдвиговой волны направлено в сторону от вертикального направления, что позволяет исследовать точку, находящуюся непосредственно под устройством. Также мы видим начинающую оформляться волну Релея.
Рис. 4. Распределение поля скоростей в момент времени £ = 0.000404
На рисунке 5 приведено сравнение с экспериментом. После вычитания сигналов, полученных при ударе справа и при ударе слева, на графике компоненты ьх (сплошная линия на правом рисунке) в явном виде можно наблюдать поперечные волны. Заметим, что при сложении сигналов (прерывистая линия) эти волны сокращаются. Мы видим качественное совпадение профиля еейемоеигнала, уходящего от нашего генератора в горную породу. Дальнейшие возмущения, которые мы видим на правом графике, являются откликом от установленного ниже препятствия. При моделировании в реальных масштабах отклик будет приходить позже, и сейсмограмма будет полностью аналогична экспериментальной.
Рис. 5. Эксперимент и показания на датчике, компонента ух. Слева сигнал, полученный при испытаниях на земле, посередине на асфальте, справа данные, полученные при расчете
3. Заключение
Для численного решения задач динамичеекого деформирования при наличии подвижных контактирующих поверхностей в данной работе адаптированы сеточно-характериетичеекий метод, использующий нерегулярные треугольные двухмерные сетки и гибридную схему 1 2 порядка, и поиск контактирующих узлов методом триангуляции пространства между телами. Разработан корректный алгоритм для расчета сил трения скольжения и покоя.
Получены численное решение задачи качения колеса и картины распределения скоростей и напряжений. Также была рассчитана генерация сдвиговых волн методом падающего груза, проведен анализ волновой картины и данных сенсора. На волновой картине наблюдаются все волны, которые обычно встречаются на практике: продольные волны, поперечные волны и поверхностные волны Релея. Данные сенсора подтверждают эффективность метода сложения сейсмограмм.
Литература
1. Новацкий В.К. Теория упругости. М.: Мир, 1975.
2. Работное Ю.Н. Механика деформируемого твердого тела. М.: Наука, 1988.
3. Рождественский Б.Л., Ясиеико H.H. Системы квазилинейных уравнений. М.: Наука, 1978. 687 с.
4. Годунов С.К., Забродин, A.B.. Ива-нов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976, 400 с.
5. Магомедов A.M., Холодов A.C. Сеточно-характеристические численные методы. М.: Наука, 1988. 288 с.
6. Пет/ров И.В., Холодов A.C. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом /7 Ж. выч. матем. и матем. физ. 1984. Т. 24, № 5. С. 722 739.
7. Петров И.В., Холодов А. С. О регуляризации разрывных численных решений уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. — 1984. — Т. 24, № 8.
- С. 1172-1188.
8. Петров И.В., Тормасов А.Г., Холодов А.С. О численном изучении нестационарных процессов в деформируемых средах многослойной структуры // Изв. АН СССР сер. МТТ. - 1989. - С. 89-95.
9. Иванов В.Д., Пет,ров И.В., Тормасов А.Г., Пашутин Р.А., Холодов А.С. Сеточно-характеристический метод расчета процессов динамического деформирования на нерегулярных расчетных сетках // Мат. моделирование. — 1999. — Т. 11, № 7. — С. 118-127.
10. Левянт В.В., Петров И.В., Челноков Ф.В. Кластерная природа сейсмической энергии, рассеянной от зоны диффузной каверзности и трещиноватости в массивных породах // Геофизика. - 2005. - 5-19.
11. Белоцерковский О.М., Агапов П.П., Пет,ров И.Б. Численное моделирование последствий механического воздействия на мозг человека при черепно-мозговой травме // Ж. вычисл. матем. и матем. физ. — 2006. — Т. 46, № 3. — С. 1711-1720.
12. Петров И.В., Челноков Ф.В. Численное исследование волновых процессов и процессов разрушения в многослойных преградах // Ж. вычисл. матем. и матем. физ. — 2003.
- Т. 43, № 10. - С. 1562-1579.
13. Layotte Р.С. Marthor: an s-wave impulse source // Shear-wave exploration edited by S.N. Domenico and S.H. Danbom, Geophysical Development Series. — V. 1. — P. 79-96. Society of Exploration Geophvsicists.
14. Helbig K. Shear-waves — what they are and how they can be used // Shear-wave exploration edited by S.N. Domenico and S.H. Danbom, Geophysical Development Series. — V. 1. — P. 19-36. Society of Exploration Geophvsicists.
Поступим в редакцию 23.01.2012.
