CC BY
188ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
CЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИЙ МЕТОД ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ВОЛНОВЫХ ПРОЦЕССОВ В ТРЕХМЕРНЫХ ЗАДАЧАХ ДИНАМИЧЕСКОГО НАГРУЖЕНИЯ СЛОЖНЫХ КОНСТРУКЦИЙ
Петров И.Б., Фаворская А.В., Хохлов Н.И., Миряха В.А., Санников А.В., Беклемышева К.А., Голубев В.И.
Московский физико-технический институт, http://mipt.ru 141700 Московская область, г. Долгопрудный, Российская Федерация Поступила в редакцию 06.11.2014
Рассматривается численный метод для компьютерного моделирования сложных пространственных динамических процессов в гетерогенных средах, в частности, распространения упругих волн в анизотропной среде композитных материалов. В качестве математической модели рассматривается гиперболическая система уравнений сплошной линейно-упругой среды с анизотропией, используется сеточно-характеристический метод на криволинейных структурированных сетках. Реализована точная постановка граничных условий заданной внешней силы, заданной скорости границы, смешанные граничные условия и неотражающие граничные условия, контактное условие полного слипания, контактное условие свободного скольжения и контактное условие с динамическим трением и трением покоя. Метод использован для численного моделирования задач дефектоскопии элементов железнодорожного пути. Эффекттивность расчетов подтверждена сравнением с результатами, полученными разрывным методом Галеркина на неструктурированных тетраэдральных сетках. В серии численных экспериментов показана возможность диагностировать повреждения в рельсе на раннем этапе их формирования и следить за динамикой их развития.
Ключевые слова: компьютерное моделирование, сеточно-характеристический метод, структурированные тетраэдральные сетки, дефектоскопия железнодорожных путей
PACS: 02.60.CB, 02.70.-C, 02.70.DH, 89.20.FF_
Содержание
1. Введение (34)
2. особенности численного моделирования (36)
3. Сеточно-характеристический метод (36)
4. Разрывный метод галеркина на неструктурированных треугольных сетках
(37)
5. Динамическая диагностика элементов пути (38)
5.1. Распространение упругих волн в рельсе (39)
5.2. Распространение упругих волн в подвижной системе «колесо-рельс» (40)
5.3. Диагностика дефектов в рельсе на раннем этапе формирования (41)
5.4. Численное моделирование волновой картины в поврежденном рельсе (41)
5.5. Распространение ультразвуковых волн в профиле рельса с горизонтальным расслоением головки (43)
5.6. Определение наличия карстового включения в грунте под насыпью (45)
6. Заключение (45) Литература (45)
1. ВВЕДЕНИЕ
Решение современных задач моделирования пространственных динамических процессов в сложных гетерогенных средах, в том числе, в частности, задач о распространении упругих волн в анизотропной среде композитных материалов в состоянии динамического нагружения, требует введения все более усложняющихся механико-математических моделей. Так как система уравнений математической модели состояния сплошной линейно-упругой среды с анизотропией [1, 2] является гиперболической и требуется высокоточный расчет волновых процессов, оптимальным является применение сеточно-характеристического метода [3, 4, 5]. В отличие от метода конечных элементов и разрывного метода Галеркина, которые используются в большинстве прикладных пакетов для моделирования задач деформируемого твердого тела, данный метод учитывает характеристические свойства
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
исходной системы уравнений. Это позволяет моделировать распространение волн в объеме, взаимодействие волновых фронтов с границами детали, а также получать полное решение нестационарных контактных задач, что обеспечивает учет влияния всех внутренних контактных границ между средами с различной реологией. Метод позволяет получить высокое временное и пространственное разрешение и рассчитать компоненты тензора напряжений и вектора скоростей деформации в любой точке рассматриваемой конструкции. Также он позволяет исследовать влияние разрушенных зон конструкции на волновую картину распространения в ней напряжений. В рамках данного метода можно применять различные критерии разрушения и учитывать одновременно различные механизмы разрушения материала. Также метод позволяет использовать различные модели материала (приближение упругого тела, вязкоупругого, упруго-пластического и вязко-упруго-пластического), в том числе анизотропные. Используется также алгоритм расчета силы трения при решении пространственных динамических задач с использованием сеточно-характеристического метода.
Метод позволяет использовать
монотонные разностные схемы повышенного порядка точности, строить корректные численные алгоритмы на границах областей интегрирования и на контактных границах. Использование иерархических сеток
позволяет учитывать большое количество негомогенных включений в моделируемой среде. При этом используется интерполяция высоких порядков на неструктурированных треугольных и тетраэдральных сетках [6, 7], в том числе интерполяция кусочно-параболическая, интерполяция с ограничителем и гибридной монотонной. Применение сеточно-характеристического метода [8-10], а также разрывного метода Галеркина [11] для упругих сред позволяет реализовать точную постановку граничных условий заданной скорости границы, заданной внешней силы и контактных условий полного слипания и свободного скольжения. Примеры использования
сеточно-характеристического метода с
квадратичной интерполяцией с ограничителем на неструктурированных треугольных сетках для решения задач в неоднородных пространственных конструкциях можно найти в [12], с интерполяцией высоких порядков на тетраэдральных сетках — в [13].
При переходе от двумерных к пространственным задачам увеличиваются и объемы данных. Поэтому требуется применение высокопроизводительных вычислительных систем. Разработанный алгоритм был распараллелен на вычислительный кластер с оптимальным использованием системных ресурсов.
Детальное описание всех волновых процессов вблизи всех имеющихся в поставленной задаче неоднородностей требует использования достаточно подробной сетки. Чем меньше исследуемые неоднородные включения, тем больше требуется шагов по времени, а также операций на каждом временном слое. Однако в большинстве случаев неоднородности локализованы в небольшом объеме внутри области интегрирования. При данной постановке задач оптимальным является применение иерархических сеток, сгущающихся в местах расположения неоднородностей.
Но для сеточно-характеристических методов размер шага интегрирования по времени напрямую зависит от размера минимального шага по пространству. Поэтому использование обычных иерархических сеток не сократит количество шагов по времени, уменьшится только число операций на каждом временном слое.
Как показали проведенные теоретические и численные исследования, применение именно сеточно-характеристических методов позволяет использовать специальные иерархические сетки с кратным шагом. Можно делать кратный шаг не только по пространству, но и по времени и, таким образом, сокращать не только число операций, необходимых для интегрирования задач на каждом временном слое, но и число выполняемых шагов по времени в той части области интегрирования, где нет неоднородностей и используется более крупная сетка.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
2. ОСОБЕННОСТИ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ
Состояние бесконечно малого объема сплошной линейно-упругой анизотропной среды подчиняется нестационарным уравнениям теории упругости для случая трех переменных в некоторой ортонормированной системе координат (х1, х2, х3):
Р^г = ^Рц ,
р = Цук1е к! + ^ у . (1)
Здесь р — плотность среды, V — компонента: вектора скорости смещения, а. и е.. — компоненты тензоров напряжений Коши и деформации, V,. - ковариантная производная по ] координате, К — добавочная правая часть. Вид компонент тензора 4-го порядка . определяется реологией среды. Для линейно-упругого случая они имеют вид
. = Я81)8к1 + К848,1 +
В этом соотношении, которое обобщает закон Гука, Я и ¡х — параметры Ляме, а 8;. — символ Кронекера.
Первая строка в системе уравнений (1) представляет три уравнения движения, вторая — шесть реологических соотношений. Вектор искомых функций, состоящий из девяти компонент, имеет вид
и = ^ ^ аи, а12, а1з, а22, а2з, азз} Т.
Тогда перечисленные модели твердого тела допускают запись системы уравнений (1) динамики деформируемого твердого тела в матричном виде [14]:
ди ^ . ди — = / А-,
j=i
dx.
(2)
A =
L — матрицы размера 9х9:
0 0 0 1 0 0 0 0 0
p 1
0 0 0 0 0 0 0 0
p 1
0 0 0 0 0 0 0 0
Р
-(Я+2ц) 0 0 0 0 0 0 0 0
0 -ц 0 0 0 0 0 0 0
0 0 -ц 0 0 0 0 0 0
-Я 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
-Я 0 0 0 0 0 0 0 0
A =
0 0 0 0 -— 0 0 0 0
p
0 0 0 0 0 -— 0 0 0
p
0 0 0 0 0 0 -— 0 0
p
0 -Я 0 0 0 0 0 0 0
-ц 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 -(Я + 2ц) 0 0 0 0 0 0 0
0 0 -ц 0 0 0 0 0 0
0 0 0 0 0 0 0
A =
0 -Я 0 0
0
0 0 0
0 0 -— 0 0 0
p
0 0 0 -— 0 0
p
1
0 0 0 0 0 0 0 Р 0
0 0 -Я 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
-ц 0 0 0 0 0 0 0 0
0 0 -Я 0 0 0 0 0 0
0 -ц 0 0 0 0 0 0 0
0 0 -(Я+2ц) 0 0 0 0 0 0
Система уравнений (2) является гиперболической, поэтому каждая из матриц А. имеет 9 вещественных собственных значений и базис из собственных векторов.
Данная запись является канонической формой записи системы уравнений, принятой в вычислительной математике для построения сеточно-характеристических разностных схем.
3. СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИЙ МЕТОД
При численном моделировании задач
динамики деформируемого твердого тела
применение сеточно-характеристического
метода начинается с применения метода
расщепления по пространственным
координатам, в результате чего имеем три
одномерные системы
ди „ ди , „ „
— = А-, у = 1,2,3.
дГ дх
а дху (3)
Каждая из этих систем является гиперболической и обладает полным набором
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
собственных секторов с действительными собственными значениями, поэтому каждую из систем можно переписать в виде ди ^ А ди
— = П/Л,. ПЛ,-,
д/ 1 11 1 дх .
где матрица ^ — матрица, составленная из собственных векторов, Л — диагональная матрица, элементами которой являются собственные значения. Для всех координат матрица Л имеет вид
Л = diag{ с, -с, с, -с, с, -с, 0, 0, 0},
ср =у1 (Л + 2/и) / р где р * — продольная скорость
звука в среде, с ="^и/ р — поперечная скорость звука.
После замены переменных V = ^и каждая из систем (3) распадается на девять независимых скалярных уравнений переноса (индекс ] далее опускается, где это возможно)
5у . ду
— + Л— = 0. д/ дх
Одномерные уравнения переноса решаются с помощью метода характеристик, либо обычными конечно-разностными схемами.
После того как все компоненты V перенесены, восстанавливается само решение
и
п+1 =
«Г = К -ст(А1 -ст(Д2 -ст(Аэ -оА4) ) ),
Д1 = ^4(-2«:+2 + 1б«т+1 - 16м:-1 + 2С2),
Д2 = 24(-м:+2+!б«:+1 - 30«: +16«:-1 --2), д з = ¿(2«:+2 - 4«:+1 + 4«: - 2«:-2), Д4 = ¿(м:+2 - 4«:+1 + 6«: - 4«^ +
Кроме того, используется сеточно-характеристический критерий монотонности [18]. В случае положительного значения X критерий монотонности имеет следующий вид:
Ш1И
ш {«:,«:-1}< «г < шах {«:,«:-1}.
Для отрицательных значений X он будет симметричен. В простейшей реализации применяется ограничитель вида
ми+1 = ■ :
ип {«
«"„+1 > шах {«:,
ип {«
В данной работе использовались Т^О-разностные схемы [15] 2-го порядка точности. В программе реализовано 15 различных лимитеров [16] — ограничителей наклона, обеспечивающих монотонность решения с сильными разрывами. В расчетах в основном использовался ограничитель зиретЬев,
предложенный в [17]: фл (г) = шах[0, ш1п(2г, 1), ш1п(г, 2)].
Также использовались сеточно-
характеристические монотонные разностные схемы, принцип построения которых описан в [18]. В программе реализованы схемы 2-4 порядка точности, большинство расчетов проводилось с использованием схемы 4-го порядка точности. Приведем ее для численного решения одномерного линейного уравнения упругости
и + Хи = 0, о = Хт/ Ь\
В [18] показано, что данный ограничитель сохраняет 4-й порядок в областях, где решение ведет себя достаточно гладко (выполняется характеристический критерий), в случае больших градиентов решения порядок схемы понижается до 3-го.
4. РАЗРЫВНЫЙ МЕТОД ГАЛЕРКИНА НА НЕСТРУКТУРИРОВАННЫХ ТРЕУГОЛЬНЫХ СЕТКАХ
Для достаточно большого класса задач дефектоскопии, травматологии, сейсморазведки и др. при численном моделировании волновых процессов в контактирующих акустической и упругой средах, а также постановке контактных условий между ними возможно использование на неструктурированных тетраэдральных сетках [19] разрывного метода Галеркина (развитие классического метода конечных элементов для решений с резкими скачками или разрывами), который известен как численный метод с высоким порядком сходимости по пространственным координатам и времени. Метод Галеркина был программно реализован под условия задач настоящей работы, показал хорошие показатели точности и производительности и явился критерием, по которому оценивалась эффективность сеточно-характеристического метода.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
При этом использовался подход [20], суть которого заключается в аналитическом решении задачи Римана [21] на контакте акустической и упругой сред. Разрывный метод Галеркина позволяет использовать полиномы высокого порядка, тем самым добиваясь высокого порядка сходимости по пространственным координатам, что является критичным при моделировании волновых процессов в гетерогенных средах В качестве системы базисных полиномов использовались полиномы Лагранжа [22]. Интегрирование по времени проводилось с помощью 7-стадийного интегратора Дорманда-Принса [23, 24] 5-го порядка точности (разновидность классического метода Рунге-Кутты 4-го порядка) — адаптивного алгоритма с автоматическим изменением порядка и шага по времени. Он рассчитывает два решения в конце шага и по разнице выдает рекомендацию по уменьшению/увеличению шага интегрирования.
5. ДИНАМИЧЕСКАЯ ДИАГНОСТИКА ЭЛЕМЕНТОВ ПУТИ
В настоящей работе сеточно-характеристический метод положен в основу численного решения такой трехмерной задачи динамического нагружения сложной конструкции как задача дефектоскопии элементов железнодорожного пути [25].
В работах [26, 27] приведены обзоры дефектов в рельсах и существующих методов дефектоскопии. К таковым методам относятся: ультразвуковой метод (ultrasonic), магнитодинамический метод (magneticfluxleakage), метод вихревых
токов (pulsed eddy current technology), автоматизированный визуальный метод (automated visual inspection) и рентгеновская дефектоскопия (radiography). Возможно также комбинирование перечисленных подходов для повышения качества и достоверности обнаружения дефектов, например, совмещение автоматического визуального метода с методом вихревых токов.
Системы дефектоскопии, основанные на автоматизированном визуальном методе, способны оценивать качество профиля рельса, измерять степень износа, стыковой зазор,
изменение местоположения шпал, отсутствие балласта, отсутствие болтов, поверхностные повреждения, включая усталостные повреждения из-за контакта «колесо-рельс», волнообразные деформации головки рельса. Скорость работы систем такого класса варьируется в широком диапазоне (от 1 до 320 км/ч) в зависимости от требуемого разрешения изображений и типа детектируемых повреждений. Однако при использовании визуального метода не представляется возможным обнаружение внутренних дефектов.
Наиболее распространен ультразвуковой метод, зарекомендовавший себя как качественное средство высокоскоростного (до 70 км/ч) обнаружения глубоких поверхностных (глубже, чем 4 мм) и внутренних дефектов, особенно в головке и шейке рельса. Данный метод продолжает активно развиваться, и на настоящий момент известны следующие его модификации: EMAT (electro-magnetic acoustic transducer), long-range ultrasonic, ultrasonic phased array и laser ultrasonic. В его развитие большой вклад вносит компьютерное моделирование физических процессов (распространения упругих волн), происходящих в материале рельса. Необходимо отметить, что существуют разнообразные подходы к проведению численного эксперимента. Так, например, в работе [28] для моделирования использовался метод конечных элементов. Сопоставление полученных результатов с последующим натурным экспериментом позволило
подтвердить возможность оценки остаточных напряжений в сварных соединениях. В работе [29] с использованием масс-пружинной модели упругого тела производилось моделирование распространения упругих волн для одномерного, двумерного и трехмерного случая. Авторами [30] предложены модификации метода конечных разностей во временной области. Они используют разделение продольных и поперечных волн для наглядного отображения вычисленных динамических полей скоростей. В работе [31] численно получен отклик от поперечного дефекта в головке рельсы. Рассмотрена его зависимость от размера и ориентации дефекта. Авторами [32] использовался полуаналитический
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
конечно-элементный метод (semi-analytical finite element method, SAFE) для моделирования распространения волн в волноводах произвольного поперечного сечения. В работе [33] проведено моделирование распространения упругих волн в рельсе, получены динамические картины скоростей и напряжений при наличии различных дефектов. В работе [34] с применением метода SAFE проведено моделирование распространения упругих волн в головке рельсы, генерируемых лазером. Анализ результатов позволил выделить отдельные моды, по всей видимости, чувствительные к некоторым типам дефектов головки рельса. Авторами [35] предложен метод анализа результатов моделирования распространения упругих волн в подошве рельса для детектирования точечных коррозий. В работе [36] проведено численное исследование специального типа волн - critically refracted longitudinal waves. Оно позволило подобрать оптимальные параметры генератора для полевых измерений.
Необходимо, однако, отметить, что в перечисленных работах компьютерное моделирование проводилось с использованием коммерческого закрытого программного обеспечения (ABAQUS, ANSYS и т.д.). Отсутствие детального понимания применяемых вычислительных алгоритмов, а также возможности оценки корректности расчетов существенно снижает ценность полученных практических результатов. В настоящей работе изложен сеточно-характеристический подход к моделированию процесса ультразвуковой дефектоскопии, на основе которого разработан программный комплекс, позволяющий провести компьютерный эксперимент. Приводятся результаты расчетов [25] с его использованием.
5.1. Распространение упругих волн в рельсе
На рис. 1 изображена криволинейная структурная расчетная сетка для расчета распространения упругих волн в рельсе. На рис. 2 показана волновая картина, возникающая при прохождении упругой волны через рельс. Расчет был проведен с помощью сеточно-характеристического метода. На рис. 2а отображено начальное возмущение, на рис. 2b,c,d продемонстрировано дальнейшее распространение упругой волны в профиле рельса.
Рис. 1. Рельс. Криволинейная структурная расчетная сетка.
Профиль.
На рис. 3 изображена характерная волновая картина, возникающая при распространении в рельсе упругой волны. На рис. 2-4 цветом показан модуль скорости.
Также было проведено сравнение расчетов сеточно-характеристическим методом на криволинейных структурных сетках и методом Галеркина на неструктурированных треугольных сетках с теорией о распространении упругих
a b
¥ I ?
c d
Рис. 2. Распространение упругих волн в рельсе. Профиль.
Рис. 3. Распространение упругих волн в рельсе. Характреная волновая картина.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Рис. 4. Распространение упругих волн в рельсе. Сравнение методов: сверху - сеточно-характеристический метод на криволинейных структурированных сетках; внизу - метод
Галеркина на неструктурированных сетках. волн в стержнях. На рис. 4 приведено сравнение
данных методов между собой.
На рис. 4 (сверху) изображена волновая картина, полученная сеточно-характеристическим методом на криволинейной структурной сетке, ниже — разрывным методом Галеркина на неструктурированной треугольной сетке. Для обоих расчетов использовались криволинейная структурная и неструктурированная треугольная сетки соответственно, состоящие из 80000 узлов, границы области интегрирования свободные. Для обоих методов расхождение расчетных скоростей распространения звука в рельсах с теоретическими значениями составляет менее 4%.
Приведем здесь с некоторыми дополнениями обзор результатов [25] численного моделирования качения колеса, прохождения волн через рельсы с горизонтальным и вертикальным расслоением головки, поперечной трещиной в головке и трещиной в шейке рельса, а также сравнения распространения упругих волн в железнодорожном полотне и насыпи в случае наличия в насыпи карстового включения и без него.
5.2. Распространение упругих волн в подвижной системе «колесо-рельс»
На рис. 5 изображена расчетная треугольная сетка для численного моделирования
Рис. 6. Распространение упругих волн в подвижной системе
"рельс — колесо". распространения упругих волн в подвижной системе «колесо-рельс». Более детально показано место контакта.
На рис. 6 изображена волновая картина, полученная при моделировании распространения упругих волн в подвижной системе «колесо-рельс». Результат был получен с помощью сеточно-характеристического метода на неструктурированных
тетраэдральных сетках с использованием гибридной монотонной интерполяции. Оттенками серого показана компонента тензора напряжений о .
На рис. 7 приведено распределение охх в зоне динамического контакта между колесом и рельсом. Данный контакт рассчитывается с помощью контактного условия трения.
На рис. 8 представлены зависимости модуля скорости от времени (сейсмотрассы) для
шт
Рис. 5. Подвижная система "рельс-колесо". Расчетная сетка.
Рис. 7. Расчет динамического контакта между колесом и рельсом.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Рис. 8. Зависимости модуля скорости от времени (сейсмотрассы), зафиксированные в неповрежденном и
поврежденных рельсах. неповрежденного рельса, рельса с вертикальным расслоением головки, горизонтальным расслоением головки и поперечной трещиной в головке. Использовался сеточно-характеристический метод. Сплошной линией изображена сейсмотрасса для неповрежденного рельса, пунктирной — при вертикальном расслоении головки, штрихпунктирной — при горизонтальном, мелкой пунктирной — при поперечной трещине в головке.
5.3. Диагностика дефектов в рельсе на раннем этапе формирования
Интересна возможность диагностировать повреждения в рельсе [37] на раннем этапе их формирования и следить за динамикой их развития. Это позволяет точнее строить прогнозы по поводу плановых замен объектов инфраструктуры и детектировать небольшие дефекты, которые могут оказаться опасными при увеличении скоростей движения железнодорожного состава.
Были рассмотрены два типа дефектов: кластер трещин длиной 0.5 мм в головке рельса, состоящий из 300 равномерно распределенных трещин и кластер пор диаметром 0.5 мм в
Рис 10. Волновая картина, возникающая при прохождении упругих волн в рельсе с кластером пор.
головке рельса, состоящий из 300 равномерно
распределенных пор.
На рис. 9 изображена тетраэдральная сетка,
сгущающаяся вокруг пор, а на рис. 10 — волновая
картина, возникающая при прохождении
упругих волн через кластер пор. На рис. 10
градацией серого показан модуль скорости.
Данная серия численных экспериментов была
выполнена с помощью разрывного метода
Галеркина на неструктурированных треугольных
сетках.
5.4. Численное моделирование волновой картины в поврежденном рельсе
На рис. 11, 12 приведены фотографии четырех типов повреждений, возникающих в рельсах. Данные взяты из книги [37], где дан более подробный перечень дефектов и повреждений рельсов.
На рис. 11 (слева) изображено вертикальное расслоение головки, код дефекта 30в 1-2. На рис. 11 (справа) представлено горизонтальное расслоение головки, код дефекта 30г 1-2. На рисунке 12 (слева) изображена поперечная трещина в головке, код дефекта 20 1-2. На рис. 12
Поперечное сечение
Поперечное сечение
Ж Ж
Общий вид
Рис. 9. Треугольная сетка, сгущающаяся вокруг пор <
Рис. 11. Расслоение головки: слева - вертикальное, код дефекта 30в 1-2, справа - горизонтальное, код дефекта 30г 1-2.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Рис. 12. Поперечная трещина в головке (слева), код дефекта 20 1-2; трещина в шейке рельса (справа), код дефекта 53 1-2. (справа) показана трещина в шейке рельса, код дефекта 53 1-2.
На рис. 13 приведены результаты численного моделирования отражения упругих волн от вертикального расслоения
головки в начальной стадии формирования — рис. 13a,b,c,d; от горизонтального — рис 13e,f,g,h. На рис. 14a,b,c,d — от поперечной трещины в головке, также в начальной стадии формирования. На рис. 14e,f,g,h для сравнения показана волновая картина в рельсе без повреждений в те же моменты времени. На рис. 13а,е и 14а изображено отражение упругой волны от дефекта; на рисунках 13 b,f и 14b — дальнейшее формирование первого отклика; на рисунках 13 c,g и 14с — формирование второго отклика от отраженных от края рельса упругих волн; на рисунках 13d,h и 14d — дальнейшее формирование второго отклика; на рис. 13-14 цветом показан модуль скорости.
На рис. 15 приведены зависимости v(t) (сейсмотрассы), зафиксированные в
"о
Я 1 ' Я 1
"ТЯк
о 1 1 у
ЯЯ г
£ ь Рис. 13. Распространение упругих волн в рельсе с расслоением
головки: вертикальным - а,Ь,е,^, горизонтальным - е,/,£,Ь; с
поперечной трещиной в головке - ^1,1,т.
V ' -, 5 * У • •л
и и
Г. .,
£ ь Рис. 14. Распространение упругих волн в рельсе с поперечной
трещиной в головке - ^1,1,т; в неповрежденном рельсе - п,о,р,г.
b
a
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
повреждение выражено уже ярче. Цветом показан модуль скорости.
5.5. Распространение ультразвуковых волн в профиле рельса с горизонтальным расслоением головки
Проведено численное моделирование распространения ультразвуковых волн в профиле железнодорожного рельса Р-65. Предполагалось наличие дефекта с кодом 30Г [37] (горизонтальное расслоение металла головки рельса) ровно посередине головки рельса.
Пьезоэлемент, возбуждающий сигнал, моделировался с помощью приложения к площадке размером 15 мм на горизонтальной поверхности головки рельса внешней силы, изменяющейся по синусоидальному закону (длина волны излучателя составляет 1 мм, частота излучателя равняется 6.25 МГц). Амплитуда внешней силы — 1 Н. Начальная фаза - 0. Материал рельса — сталь, с параметрами X = 146.1 х 109 Па, ^ = 79.3Х109 Па, р = 7.8Х103 кг/м3, ср = 6250 м/с. Гармонический сигнал модулировался прямоугольным меандром длительностью, равной 10 периодам.
Приемник сигнала предполагался точечным и расположенным на поверхности головки рельса, на оси симметрии. Приемник сигнала измерял компоненты скорости. Расчет выполнялся в двумерной постановке. Использовалась неструктурированная треугольная расчетная сетка с 104 ячеек, изображенная на рис. 18. Сетка построена с использованием библиотеки triangle. Применялся разрывный метод Галеркина с полиномами Лагранжа 6-го порядка.
Рис. 15. Зависимости v(t) (сейсмограммы), зафиксированные в поврежденных рельсах.
поврежденных рельсах. Красным цветом изображена сейсмотрасса для неповрежденного рельса, зеленым — при вертикальном расслоении головки, фиолетовым — при горизонтальном, синим — при поперечной трещине в головке. Видно, что на начальных этапах сейсмотрассы существенно
различаются для неповрежденного рельса и рельса с повреждениями.
На рис. 16 изображено формирование отклика от трещины в шейке рельса. На верхнем рисунке оно соответствует начальной стадии повреждения, а на нижнем
Рис. 16. Формирование отклика от трещины в шейке рельса: сверху - начальная стадия повреждения, ниже - развитая стадия повреждения
На рис. 17 изображены зависимости горизонтальной компоненты скорости vJJ) (сейсмотрассы): синим цветом — для неповрежденного рельса, зеленым — при наличии повреждения в начальной стадии, красным - при наличии ярче выраженного повреждения.
Рис. 17. Зависимости у(1) (сейсмограммы), зафиксированные в неповрежденном рельсе и рельсах с трещиной в шейке разной длины.
Рис. 18. Распространение ультразвуковых волн в рельсе с горизонтальным расслоением головки. Расчетная сетка.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Интегрирование по времени проводилось с помощью 7-стадийного метода Дорманда-Принса 5-го порядка с адаптивным шагом по времени. Расслоение располагалось на расстоянии 22.5 мм от поверхности головки и моделировалось с помощью свободной границы.
На рис. 19а,Ь,с изображено формирование диагностирующей волны, момент времени — 1 мс. На рис. 19й,в,/ приведены формирующиеся сферическиепоперечные волнывмоментвремени 3 мс и показано дальнейшее распространение диагностирующей волны. На рис. 19g,h,i показано формирование отклика от горизонтального расслоения головки в момент времени 4 мс. На рис. 19]^,1 изображено распространение отклика
Рис. 20. Распространение ультразвуковых волн в рельсе с горизонтальным расслоением головки. Зависимость (сейсмотрасса).
к датчику, расположенному на поверхности рельса, в момент времени 7 мс. На рис. 19т,пр приведено отражение отклика от поверхности рельса в момент времени 8 мс. Рис. 19 q,r,s изображает сформировавшуюся отраженную от поверхности рельса волну в момент времени 10 мс. На рис. цветом показан модуль
скорости, на рис. 19Ь,е^^,п,г — горизонтальная компонента скорости, на рис. 19с/,г,1р^ — вертикальная компонента скорости.
На рис. 20 приведена зависимость V ($) (сейсмотрасса), а на рис. 21 — v(j). Отклик следовало ожидать спустя время 2х22.5/ср = 7.2 х10-6 с, что согласуется с расчетными данными.
Обратим внимание на два факта. Во-первых, амплитуда отраженного от дефекта сигнала на приемнике вдвое больше, чем у испущенного сигнала. Это связано с тем, что скорость не меняет фазы при отражении от свободной поверхности и интерферирует с собой же. Во-вторых, в интервале от 2.5х10-6 с до 4х10-6 с наблюдается возмущение, сравнимое по амплитуде с сигналом источника.
|па*т
q г s
Рис. 19. Распространение ультразвуковых волн в рельсе с горизонтальным расслоением головки: а,Ь,с - формирование волны, сС,е,/ - распространение волны, - формирование отклика, - распространение отклика, т,пр - отражение от границы, q,г,s - распространение отраженной волны.
Г
Рис. 21. Распространение ультразвуковых волн в рельсе с горизонтальным расслоением головки. Зависимость V (сейсмотрасса).
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Это связано с возникающей дифракционной картиной на краях пьезоэлемента. 5.6. Определение наличия карстового включения в грунте под насыпью
На рис. 22 приведено начальное возмущение для расчета волновой картины в насыпи с карстовым включением (а) и без него (Ь). На рис. 22а-Ь цветом показан модуль скорости. На рис. 22е,с1 отображено дальнейшее распространение упругих волн в рельсе, насыпи и грунте. На рис. 22в,/ отображено формирование первого отклика от карстового включения, а на рис. 22^,Ь — второго.
с иид ■
WW
го 43
g h
Рис. 22. Распространениеупругих волн в насыпи: a, b - начальное возмущение, с каверной (a), без каверныг (b); c, d - распространение возмущения, с каверной (с), без каверныг (d); e, f - формирование первого отклика, с каверной (e), тот же момент времени без каверныг f);g, h - формирование второго отклика, с каверной (g), тот же момент времени без каверныг (h).
На рис. 23, 24 приведены зависимости v(t) (сейсмотрассы), полученные с левого и правого рельса соответственно. На обоих графиках
200 2» 2«
Рис. 23. Зависимости v(t) (сейсмограммы), зафиксированные на левом рельсе.
Рис. 24. Зависимости v(t) (сейсмограммы), зафиксированныге
на правом рельсе красный цвет соответствует расчету без карстового включения, а синий — с карстовым включением.
4. ЗАКЛЮЧЕНИЕ
Сеточно-характеристический метод высокой точности на тетраэдральных иерархических сетках с использованием интерполяции высоких порядков и кратным шагом по времени развит и использован для численного моделирования волновых процессов в трехмерных задачах динамического нагружения сложных
контрукций. Изложены основные элементы метода. Разработанный алгоритм реализован на высокопроизводительных вычислительных системах. Проведены серии численных экспериментов по использованию метода для решения задач дефектоскопии железнодорожных путей. Показано, что метод позволяет получить высокое временное и пространственное разрешение и рассчитать компоненты тензора напряжений и вектора скоростей деформации в любой точке рассматриваемой конструкции. Также он позволяет исследовать влияние разрушенных зон на волновую картину повреждений. Применение накопленного в данной области опыта позволит существенно ускорить разработку эффективных методов неразрушающей дефектоскопии.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 13-07-13169 офи_м_РЖД.
ЛИТЕРАТУРА
1. Новацкий ВК. Теорияупругости. М., Мир, 1975, 872 с.
2. Кондауров ВИ, Фортов В.Е. Основы термомеханики конденсированной среды. М., МФТИ, 2002, 123 с.
3. Магомедов КМ. Холодов АС. О построении разностных схем для уравнений гиперболического типа на основе характеристических соотношений. Журнал вычислительной математики
46 петров KR, фаворская А.к, хохлов ки., миряхл ВА миФПРМЛММПММиР ТРУМПППГИИ
САННИКОВ А.В., БЕКЛЕМЫШЕВА К.А., ГОЛУБЕВ В.И. ИНФОРМАЦИОННЫЕ 1ЬЛНОПО1ИИ
и математической физики, 1969, 9(2):373-386.
4. Магомедов КМ, Холодов АС. Сеточно-характеристические численные методы. М., Наука, 1988, 123 с.
5. Иванов ВД, Кондауров ВИ, Петров ИБ, Холодов АС. Расчет динамического деформирования и разрушения упругопластических тел сеточно-характеристическими методами. Матем. моделирование, 1990, 2(11):10-29.
6. Петров ИБ, Лобанов АИ. Лекции по вычислительной математике. М., Интернет-ун-т информац. технологий, 2006, 123 с.
7. Петров ИБ, Фаворская АБ. Библиотека по интерполяции высоких порядков на неструктурированных треугольных и тетраэдральных сетках. Информационные технологии, 2011, 9:30-32.
8. Петров ИБ, Фаворская АБ, Санников АВ, Квасов ИЕ. Сеточно-характеристический метод с использованием интерполяции высоких порядков на тетраэдральных иерархических сетках с кратным шагом по времени. Математическое моделирование, 2013, 25(2):42-52.
9. Муратов МВ, Петров ИБ, Санников АВ, Фаворская АВ. Сеточно-характеристический метод на неструктурированных тетраэдральных сетках. Журнал вычислительной математики и математической физики, 2014, 54(5):821-832.
10. Петров ИБ, Хохлов НИ. Моделирование задач 3D сейсмики на высокопроизводительных вычислительных системах. Матем. моделирование,
2014, 26(1):83-95.
11. Käser M, Dumbser M. An arbitrary high-oder discontinuous Galerkin method for elastic waves on unstructured meshes - I.The two-dimensional isotropic case with external source terms. Geophysical Journal International, 2006, 66(2):855-877.
12. Квасов ИЕ, Петров ИБ, Челноков ФБ. Расчет волновых процессов в неоднородных пространственных конструкциях. Математическое моделирование, 2009, 21(5):3-9.
13. Квасов ИЕ, Петров ИБ, Санников АВ, Фаворская АВ. Сеточно-характеристический метод высокой точности на тетраэдральных иерархических сетках с кратным шагом по времени. Компьютерные исследования и моделирование, 2012, 3(1):161-171.
14. Петров ИБ, Холодов АС. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом. Журнал вычислительной математики и математической физики, 1984, 24(5):722-739.
15. Ami Harten. High resolution schemes for hyperbolic conaervation laws. J. of Comp. Physics, 1997, 135(2):337-365.
16. Петров ИБ, Хохлов НИ. Сравнение TVD лимитеров для численного решения уравнений динамики дефомируемого твердого тела сеточно-
характеристическим методом. Математические модели и задачиуправления. М., МФТИ, 2011, с. 104-111.
17. Roe PL. Characteristic-Based Schemes for the Euler Equation. Annual Review of Fluid Mechanics, 1986, 18:337-365.
18. Холодов АС, Холодов ЯА. О критериях монотонности разностных схем для уравнения гиперболического типа. Журнал вычислительной математики и математической физики, 2006, 46(9):1560-1588.
19. Миряха ВА, Санников АВ, Петров ИБ. Численное моделирование волновых процессов в гидроупругих задачах разрывным методом Галеркина на неструктурированных треугольных сетках. ВестникБалтийского федеральногоуниверситета им. И. Канта. 2014, 10:16-20.
20. Wilcox LC, Stadler G, Burstedde C. et al. A highorder discontinuous Galerkin method for wave propagation through coupled elastic—acoustic media. J. of Comp.Phys, 2010, 229:9373-9396.
21. LeVeque RL. Finite volume methods for hyperbolic problems. Cambridge, 2002.
22. Hesthaven JS, Warburton T. Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Texts in Applied Mathematics. Springer, 2008, vol. 54.
23. DormandJR, Prince PJ. A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 1980, 6(1):19-26, doi:10.1016/0771-050X(80)90013-3
24. Hairer E, N0rsett SP, Wanner G. Solving ordinary differential equations I: Nonstiff problems. Berlin, New York: Springer-Verlag, 2008.
25. Петров ИБ, Голубев ВИ, Миряха ВА, Хохлов НИ, Фаворская АВ, Санников АВ, Беклемышева КА. Динамическая диагностика элементов пути. Техника железных дорог, 2013, 4:64-77.
26. Cannon DF, Edel K-O, Grassie SL, Sawley K. Rail defects: an overview. Fatigue & Fracture of Engineering Materials & Structures, 2003, 26(10):865-886.
27. Papaelias MPh, Roberts C, Davis CL. A review on non-destructive evaluation of rails: State-of-the-art and future development. Journal of Rail and Rapid Transit, 2008, 222(4):367-384.
28. Джавади Я, Наджафабади МА, Ахлаги М. Оценка остаточных напряжений в сварных соединениях из разнородных компонент с использованием моделирования методом конечных элементов и измерения головной ультразвуковой волны. Дефектоскопия, 2012, 9:48-61.
29. Сыч ТВ, Герасимов СИ, Кулешов ВК. Моделирование распространения акустических волн методом конечных элементов. Дефектоскопия, 2012, 3:3-9.
30. Бархатов ВА. Моделирование ультразвуковых волн методом конечных разностей во временной области. Двумерная задача. Оптимальные алгоритмы. Анализ погрешностей. Поглощающие области вблизи границ сеток. Дефектоскопия, 2009, 6:76-82.
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
31. Bartoli I, di Scalea FL, Fateh M, Viola E. Modeling guided wave propagation with application to the longrange defect detection in railroad tracks.NDT&E International, 2005, 38(5):325-334.
32. Bartoli I, Marzani A, di Scale'a FL, Viola E. Modeling wave propagation in damped waveguides of arbitrary cross-section. NDT8cE International, 2006, 295(3-5):685-707.
33. Zumpano G, Meo M. A new damage detection technique based on wave propagation for rails. Intern. Journal of Solids and Structures, 2006, 43:1023-1046.
34. Coccia S, Bartoli I, Marzani A, di Scalea FL, Salamone S, Fateh M. Numerical and experimental study of guided waves for detection of defects in the rail head. NDT&E International, 2011, 44:93-100.
35. Cerniglia D, Pantano A, Vento MA. Guided wave propagation in a plate edge and application to NDI of rail base. Journal Nondestuct Eval., 2012, 31:245-252.
36. Chaki S, Ke W, Demouveau H. Numerical and experimental analysis of the critically refracted longitudinal beam. Ultrasonics, 2013, 53(1):65-69.
37. Грицык В.И. Дефекты рельсов железнодорожного пути. М., Маршрут, 2005, 80 с.
Петров Игорь Борисович
д.ф.-м.н, проф., чл.-корр. РАН и РАЕН, зав. тф. информатики
МФТИ, ф-т управления и прикладной математики
141701 Моск. обл., г.Долгопрудный, Россия
8 495 408-6695, petrov@mipt.ru
Фаворская Алена Владимировна
ассистент кафедры информатики
GRID-CHARACTERISTIC METHOD
МФТИ, ф-т управления и прикладной математики 141701 Моск. обл., г. Долгопрудный, Россия 8 967 0083986, aleanera@yandex.ru Хохлов Николай Игоревич
к.ф.-м.н, заведующий лабораторией
МФТИ, ф-т управления и прикладной математики
141701 Моск. обл., г. Долгопрудный, Россия
8 916 0105212, k_h@inbox.ru
Миряха Владислав Андреевич
ассистент кафедры информатики
МФТИ, ф-т управления и прикладной математики 141701 Моск. обл., г. Долгопрудный, Россия 8 915 1810630, vlad.miryaha@gmail.com Санников Александр Владимирович
ассистент кафедры информатики
МФТИ, ф-т управления и прикладной математики 141701 Моск. обл., г. Долгопрудный, Россия 8 906 1096720, donxenapo@gmail.com Беклемышева Катерина Алексеевна
к.ф.-м.н., ассистент кафедры информатики МФТИ, ф-т управления и прикладной математики 141701 Моск. обл., г. Долгопрудный, Россия 8 963 6429147, amisto@yandex.ru Голубев Василий Иванович ассистент кафедры информатики
МФТИ, ф-т управления и прикладной математики 141701 Моск. обл., г. Долгопрудный, Россия 8 926 4572707, w.golubev@mail.ru
FOR NUMERICAL MODELING
OF WAVE PROCESSES IN THREE-DIMENSIONAL PROBLEMS OF DYNAMIC LOADING OF COMPLEX STRUCTURES
Petrov Igor B., Favorskaya Alena V., Khokhlov Nikolay I., Miryakha Vladislav A., Sannikov Alexander V., Beklemysheva Katerina A., Golubev Vasiliy I.
Moscow Institute of Physics and Technology, Department of Computer Science, http://www.mipt.ru
141701 Dolgoprudny, Moscow region, Russian Federation
petrov@mipt.ru
This paper is to inform about the numerical method of computer modeling of wave propagation in three-dimensional problems of dynamic loading of complex structures. As a method of modeling uses grid-characteristic method. This method uses unstructured tetrahedral hierarchical meshes, a multiple time step and the high-order interpolation, it has the precise formulation of contact conditions. The use of this grid-characteristic method makes it possible to use the multiple time step and thereby increase productivity and significantly reduce the computation time. This method is used here for modeling the railway defectoscopy for security and timely detection of defects. The comparison of wave patterns that occur during the transmission of elastic waves in the rail obtained by grid-characteristic method on curvilinear structured grids and discontinuous Galerkin method on unstructured tetrahedral grids was done.
Keywords: grid-characteristic method, computer simulation, discontinuous Galerkin method, defectoscopy of railway.
PACS: 02.60.Cb, 02.70.-c, 02.70.Dh, 89.20.Ff
Bibliography - 37 references Received 06.11.2014 RENSIT, 2015, 7(1):34-47_DOI: 10.17725/RENSITe.2015.07.034
