CC BY
36
МЕТОДЫ ГЕОГРАФИЧЕСКИХ ИССЛЕДОВАНИЙ
УДК 551.4.013
В.О. Михайлов1
ТРЕХМЕРНАЯ МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ОБВАЛЬНЫХ ПРОЦЕССОВ
В статье рассмотрены основные принципы построения разработанной автором трехмерной математической модели движения обвальных масс DEBRIS, а также результаты реализации этой модели для случая Аттабадского обвала в долине р. Хунза. Движение обвала в модели описывается траекториями движения центров масс обломков, упруго взаимодействующих между собой и с поверхностью ложа и испытывающих при этом влияние трения.
Ключевые слова: математическое моделирование, обвалы.
Введение. Математическое моделирование потоков вещества, встречающихся в естественных условиях, — один из наиболее важных и необходимых инструментов, используемых для описания динамики и определения характеристик конкретных потоков вещества в различном фазовом состоянии. В частности, моделирование катастрофических обвалов необходимо для проведения защитных мероприятий, а также для прогноза параметров и зон охвата этими катастрофическими процессами при задаваемых условиях.
К достоинствам математического моделирования по сравнению с экспериментальным относится соответствие характерных размеров модельного и реального процессов, а также возможность моделировать процессы различного типа и генезиса, в том числе те, для которых экспериментальное моделирование осуществить затруднительно.
Существуют различные модели движения обвалов. Совокупность этих моделей можно разделить в зависимости от физической основы на две группы — дискретные и непрерывные модели. Непрерывные модели основаны на уравнениях движения и непрерывности [5] и представляют движение обвала как сплошной среды. К непрерывным относятся модели Божинского—Назарова [1], SLOPE/W [14], FLO-2D [8], LAHARZ [13], модель К. Ху [15] и многие другие. Дискретные модели основаны непосредственно на уравнении движения отдельных условных блоков вещества. Примеры дискретных моделей — модели DAN-3D, DAN-W [9] и RAMMS [12], CRSP[6], UDEC [10]. Некоторые модели, например модель Божинского—Назарова [1], дополняются также стохастическими моделями [2], в которых параметры обвала представляются как случайные величины. Автор статьи предложил специфическую модель обвала DEBRIS [3], в основе которой лежит представление обвала в виде совокупности траекторий отдельных взаимодействующих друг с другом обломков в трехмерном пространстве.
Характеристика модели DEBRIS. Модель движения обвальных масс DEBRIS (Digital Elementary
Balls & Relief Interaction Simulation, т.е. цифровая модель взаимодействия элементарных сферических частиц и поверхности) основана на численном решении дифференциального уравнения движения материальной точки, составленного в соответствии со вторым законом Ньютона. Движение обвала в таком случае описывается траекториями движения центров масс обломков. Ложе потока в модели соответствует поверхности цифровой триангуляционной модели рельефа. Частицы потока в модели имеют одинаковый размер, шарообразную форму и рассматриваются в качестве материальных точек. При движении частицы упруго взаимодействуют между собой, а также с поверхностью ложа, испытывая при этом влияние трения. В связи с рассмотрением обломков породы в качестве материальных точек в модели не отражается процесс вращения обломков, при этом под трением подразумевается трение качения.
Если образно представить модель, то движение обвала в ней имитируется движением потока горошин в лотке, форма которого соответствует рельефу выбранной территории. Горошины при этом упруго соударяются между собой, а также со стенками лотка, испытывая при этом влияние трения. Физические характеристики горошин при этом соответствуют физическим характеристикам обломков пород, а физические характеристики лотка — характеристикам ложа потока.
Физическая основа модели. Рассмотрим вывод основных рабочих формул математической модели DEBRIS, представляющих ее физическую основу. Сначала рассмотрим силы, действующие на произвольно выбранную г-ю частицу в потоке (рис. 1). На г-ю частицу действуют следующие силы:
1) сила тяжести, направленная вертикально вниз и определяемая как Fi тяжести = mi g, где g — ускорение свободного падения, тг=4/3)кЯ3рг — масса г-й частицы, рг — плотность г'-й частицы, R — радиус г'-й частицы;
2) сила нормальной реакции опоры, направленная перпендикулярно поверхности, с которой соприкасается г-я частица, и определяемая согласно закону Гука как
1 Московский государственный университет имени М.В. Ломоносова, географический факультет, научно-исследовательская
лаборатория снежных лавин и селей, инженер, аспирант; e-mail: vladg27@rambler.ru
Рис. 1. Силы, действующие на выбранную г-ю частицу: рт — сила тяжести, F¡р — сила нормальной реакции, F¡jупр — сила упругости, ртр — сила трения каче-
реакции Е^г' [(Яг' /Яг' ]пнормали
где Е — модуль Юнга для материала частицы г; — площадь соприкосновения частицы г с поверхностью ложа, зависящая от величины сжатия; Яг — радиус г-й частицы; 4 — расстояние между центром частицы г и поверхностью ложа по нормали к нему, пнормали — единичный вектор нормали к поверхности, на которой расположена г-я частица;
3) силы упругости, действующие на выбранную г-ю частицу со стороны каждой соударяющейся с ней у-й частицы и направленные в противоположную сторону от у-й частицы. Сила определяется согласно закону Гука как
^/упругости Е^у [(Яг + Яу ¿у) / (Яг + ^/)]пг/частицы,
где Е — модуль Юнга для материала частиц г и '; Бу — площадь соприкосновения частиц г и у', зависящая от величины сжатия; Я1 и Яу — радиусы г-й и у-й частиц соответственно; ¿у — расстояние между центрами частиц г и у; Пу частицы — единичный вектор, направленный из центра г-й частицы к центру у-й частицы. Следует заметить, что соударение частиц между собой и с поверхностью ложа принимается в модели абсолютно упругим. В случае относительно малой величины деформации обломков в модели принимается линейная зависимость деформации от напряжения, описываемая законом Гука. Применение закона Гука в случае крайне малых деформаций обломков требует повышенной точности измерения величины относительного сжатия обломков Дйу,, а
У
также использования бесконечно малых значений временного шага в связи с быстрым изменением величины относительного сжатия в течение времени;
4) сила сухого трения качения и вязкого трения частицы о поверхность опоры и о поверхность других частиц, направленная противоположно вектору скорости частицы и определяемая как
Р г трения [Н-г реакции + упругости) +
+ )] скорости,
где ^ — коэффициент (угол) сухого трения качения; Р реакции — модуль силы нормальной реакции опоры; ^Ру упругости — модуль суммарной силы упругости; действующей на г-ю частицу; — коэффициент вязкого трения; V — скорость г-й частицы; н1 скорости — единичный вектор скорости г-й частицы.
Отметим, что силы сухого трения качения, вязкого трения, упругости и нормальной реакции опоры приложены к поверхности частицы, что обусловливает возникновение крутящего момента. Однако ввиду рассмотрения частиц в качестве материальных точек в модели не учитывается их вращение. В модели определяется лишь движение центра массы частицы. Подчеркнем, что силы сухого трения качения, вязкого трения, упругости и нормальной реакции опоры возникают только при соударении частиц между собой или с поверхностью опоры.
Сумма указанных выше сил, действующих на г-ю частицу ^суммарная, примет следующий вид:
р = Р + Р + ХР +
суммарная тяжести р г реакции у упругости
+ Р
^ I трения •
После этого применим второй закон Ньютона для получения уравнения движения выбранной -й
частицы: ^суммарная = тгаг, 1Де Суммарная — векторная
сумма сил; тг — масса г-й частицы; аг — ускорение г-й частицы. С учетом определения ускорения как второй производной от радиус-вектора материальной точки по времени выразим из указанного уравнения ускорение частицы и перепишем полученное уравнение, спроектировав его составляющие на координатные оси X, У и Z в виде:
с1\ _ Р1СуммХ(х, у, г, О
Ж2 т1
<12У1 _ -^1_сумму(*> У-> 0
л2 т{
_^_сумм2(х, у, 0
Л2
(1)
где Р_суммХх У, г, 1), Р_суммУ(х, У, г 1) и Р_умм2<х, У, г, 0 — проекции на оси X, У и Z соответственно суммарной силы, действующей на г-ю частицу.
Полученные уравнения (1) — уравнения движения г-й частицы. В результате численного решения уравнений (1) можно получить траекторию движения г-й частицы, т.е. зависимость координат г-й частицы от времени. При численном решении использован метод конечных разностей [4] с временным шагом 0,01 с. Укажем выражения, определяющие координаты 1-й частицы в последующий момент времени хц+м , Уи+М , гц+м при известных ее координатах в текущий хи, уи , ги и предыдущий х,-д, уи-ы, моменты времени:
м
ТП)
- Ъ сумм7 . Л
, (2)
ГП:
21л+м — 21л-м
-^¿суммг (х> У,2, 0
т.-
ния
Выражения (2) — основные рабочие формулы математической модели DEBRIS.
Исходные данные модели. После рассмотрения физической и численной основы математической модели DEBRIS рассмотрим исходные данные модели. Минимальный набор исходных данных, необходимых для использования модели DEBRIS, составляют следующие данные:
1) цифровая модель поверхности начального положения вещества потока («уровня») в виде массива точек поверхности с известными координатами X,
Y и Z;
2) цифровая модель поверхности ложа потока («ложа») в виде массива точек поверхности с известными координатами X, Y и Z.
Указанные цифровые модели поверхностей строятся на основании имеющихся фактических материалов — топографических карт, результатов дешифрирования аэрофото- и космических снимков, полевых исследований, а также некоторых предположений;
3) физические характеристики вещества потока: плотность (р), коэффициент (угол) сухого трения
и вязкого трения (п), модуль Юнга (E). Значения этих характеристик в большинстве случаев можно получить из справочных таблиц свойств различных материалов.
Алгоритм работы модели. Рассмотрим собственно ход работы модели DEBRIS с заданными исходными данными.
Шаг 1 (создание поверхностей «уровня» и «ложа» потока). Сначала на основании имеющихся полей точек с заданными прямоугольными координатами X,
Y и Z, соответствующих поверхностям уровня и ложа потока, строится нерегулярная триангуляционная сеть. При ее построении используется метод кратчайшего ребра. Поверхности уровня и ложа потока при этом представляют собой совокупности физических треугольников, для каждого из которых определяется вектор нормали «нормали.
Шаг 2 (создание частиц). Далее пространство между двумя поверхностями — уровнем и ложем потока — заполняется частицами с заданными физическими характеристиками, соответствующими свойствам вещества потока. Совокупность частиц при этом представляет собой собственно толщу пришедшего в движение вещества потока. Каждой частице при этом приписывается индивидуальный порядковый номер i. Заполняющие пространство частицы в начальный момент времени упакованы максимально плотно — каждая частица во внутренней области пространства соприкасается с двенадцатью соседними частицами. При заполнении пространства частицам приписываются координаты их положения в пространстве на начальном и предшествующем временных шагах с учетом заданной начальной скорости их движения. Для каждой i-й частицы имеется следующий набор параметров: радиус (R), плотность (рг), коэффици-
ент (угол) сухого трения коэффициент вязкого трения (п), модуль Юнга (Е), координаты на начальном временном шаге (х,, у,, и на предыдущем временном шаге (х,-д* уи_Д( , 2,-д).
Шаг 3 (перемещение частиц). Затем частицы «пускаются» по поверхности слоя, соответствующего поверхности ложа потока (рис. 2). На каждом текущем временном шаге t для каждой выбранной г-й частицы определяются величины составляющих сил F¿тяжести,
^реакцию ^/упругости ^трения по координатн^1м осям X, Г
и 2 согласно приведенным выше выражениям.
Рис. 2. Представление движения частиц по поверхности «ложа» потока в модели DEBRIS
При этом значения параметров, входящих в выражения для определения векторов сил, определяются следующим образом. Определение вектора силы тяжести проводится с учетом заданных параметров. При определении вектора силы нормальной реакции для каждой i-й частицы проверяются все возможные треугольники поверхности ложа потока на предмет соприкосновения с ними i-й частицы.
При определении вектора суммарной силы упругости для каждой выбранной i-й частицы проверяются все остальные j-е частицы на предмет их соударения с выбранной i-й частицей. Соударение произошло, если расстояние между частицами меньше суммы их радиусов. В случае соударения частиц i и j вычисляются и суммируются возникающие при этом векторы сил упругости.
При определении вектора силы трения качения модуль скорости частицы v t и единичный вектор скорости частицы И;скоросхи вычисляются через разность пространственных координат частицы на текущем x,t, yt,t, zt,t и на предыдущем хц-м, y^t, zt,t-&t временном шаге.
Для каждого вычисленного вектора силы, действующей на каждую i-ю выбранную частицу, устанавливаются его проекции на координатные оси X, Y, Z. Затем определяются проекции на координатные оси X, Y, Z результирующей векторной суммы сил, действующих на каждую выбранную i-ю частицу, путем сложения соответствующих проекций сил;
таким образом получают массивы проекций на координатные оси X, Y, Z векторов результирующих сил, действующих на каждую г'-ю частицу. После этого на основании набора данных, полученных на этом шаге работы модели, определяют координаты каждой г'-й частицы на последующем (t+At) временном шаге по формулам (2) с учетом ее координат на текущем (t) и предыдущем (t-At) временных шагах.
Шаг 3 (перемещение частиц) повторяется сначала для каждой частицы. При этом определяются координаты каждой г'-й частицы на последующем шаге (t+2At) по формулам (2) c учетом уже установленных ранее ее координат на новом текущем (t+At) и новом предыдущем (t) временных шагах и т.д. В ходе работы модели осуществляется визуализация смещения частиц с течением времени последовательно на каждом временном шаге.
Математическая модель DEBRIS построена в средах Borland C++ Builder 6 (основная программа) и Microsoft Visual C/C++ 2008 (вспомогательная программа).
Реализация модели для Аттабадского обвала. Для
демонстрации работы модели DEBRIS приведен пример реализации Аттабадского обвала, произошедшего в долине р. Хунза (Каракорум, Пакистан) 4 января 2010 г. В модели использована цифровая модель рельефа участка бассейна р. Хунза, построенная автором с использованием топографической карты масштаба 1:200 000. При моделировании на момент, предшествовавший обвалу, обломки в зоне его зарождения находились на правом борту долины р. Хунза. При реконструкции площади зоны зарождения и аккумуляции использованы космические снимки ALI, опубликованные NASA [11], и наземные фотографии, полученные очевидцем обвала Инайятом Али. Для частиц обвальных масс заданы следующие параметры: радиус R = 10 м, плотность р = 2000 кг/ м3, значения коэффициентов (угла) сухого трения качения ^ = 0,3 и вязкого трения п = 0 кН-с/м, модуль
Юнга Е = 40 ГПа. Начальная скорость всех частиц принималась равной нулю, т.е. частицы вначале принимались неподвижными. Общее число обломков, заданных в модели, — 1799 шт. Движение обвальных масс визуализируется в различные моменты времени Т (рис. 3). Временной шаг At в модели принимался равным 0,01 с. К моменту времени Т = 30 с обвальные массы начали движение вниз по склону и вышли в дно долины р. Хунза, ударившись о противоположный борт долины. После выхода обвальных масс в дно, в условиях резкого уменьшения уклона частицы испытали торможение. Торможению также способствовал удар обвальных масс о противоположный борт долины. К моменту времени Т = 50 с основная часть потока остановилась, сформировав нагромождения обломочного материала на дне долины.
Таким образом, в результате моделирования получены следующие результаты. Общий период движения обвальных масс до полной остановки составил примерно 60 с. Установленная в модели скорость движения фронтальной части обвала колебалась в среднем от 50 до 60 м/с. Установленная максимальная мощность отложений обвала в долине составляла около 120 м, ширина 600 м. Общая длина аккумулятивного тела составила примерно 2—2,5 км. Величина объема отложений в долине р. Хунза в результате моделирования составила 45 млн м3. Распределение обвальных отложений, полученное в результате моделирования, близко к натурным наблюдениям Д. Пэтли [7] и специалистов организации «ФОКУС Гуманитарная помощь», занимающихся ликвидацией последствий обвала. Реальный объем отложений обвала в долине оценивается примерно в 45—50 млн м3, в модели же объем составил 45 млн м3. Параметры аккумулятивного тела в долине в целом схожи с модельными данными. Общее время движения обвала оценивается примерно в 1,5 мин против 1 мин в модели, т.е. скорость движения частиц в модели несколько больше, чем в реальности.
Рис. 3. Результаты реализации Аттабадского обвала в бассейне р. Хунза в модели DEBRIS в различные моменты времени (Т), равные соответственно ноль с (а), 30 с (б) и 50 с (в). Черным показаны частицы, имитирующие обломки породы, серое — след движения частиц по поверхности цифровой модели рельефа. Горизонтали проведены через 100 м
Заключение. В заключение необходимо отметить, что математическая модель потока вещества DEBRIS — способ решения задачи об определении параметров и результатов прохождения конкретного обвала при заданных начальных условиях. В перспективе планируется совершенствование модели, во-первых, путем повышения точности входных данных (цифровая модель рельефа, число структурных частиц) и, во-вторых, путем учета дополнительных факторов (неоднородность формы и размера частиц, фазовые переходы, дробление частиц, поступление в поток дополнительного вещества со стенок ложа).
Для практического осуществления этих задач необходимо как техническое преобразование и доработка программы модели, так и использование вычислительных машин высокой мощности. Планируется проверить работу модели посредством мо-
СПИСОК ЛИТЕРАТУРЫ
1. Божинский А.Н., Назаров А.Н. Динамика двухфазного селевого потока // Вестн. Моск. ун-та. Сер. 5. География. 1999. № 5. С. 15-19.
2. Божинский А.Н., Назаров A.M, Черноус П.А. Вероятностная модель движения снежных лавин // Там же. 2000. № 5. С. 8-12.
3. Михайлов В.О., Черноморец С.С. Трехмерное математическое моделирование ледниковых и селевых катастроф с помощью модели DEBRIS // Мат-лы VII Междунар. науч. конф. «Устойчивое развитие горных территорий в условиях глобальных изменений». Владикавказ, 2010. (CD-ROM).
4. Самарский А.А., Гулин А.В. Численные методы математической физики. М.: Научный мир, 2000. 315 с.
5. Эглит М.Э. Неустановившиеся движения в руслах и на склонах. М.: Изд-во Моск. ун-та, 1986. 96 с.
6. Colorado Department of Transportation, Colorado Rockfall Simulation Program (CRSP), URL: http://www.colo-radodot.info/programs/geotech (дата последнего обращения 20.03.2011)
7. Dave's Landslide Blog. 2010. URL: http://daveslandslide-blog.blogspot.com (дата последнего обращения: 20.03.2011)
8. FLO-2D Software. URL: http://www.flo-2d.com/ (дата последнего обращения 20.03.2011)
делирования крупных катастрофических природных процессов, имеющих различные генезис и фазовое состояние движущихся масс, т.е. не только обвалов, но и оползней и селей. На данный момент проведено моделирование для нескольких случаев реализации катастрофического процесса — Геналдонской ледниковой катастрофы в Северной Осетии, селя по р. Герхожан-Су в Кабардино-Балкарии, Уаскаранской ледниковой катастрофы в Перу, обвала-оползня Игун в Китае и упомянутого выше Аттабадского обвала в Пакистане. Результаты данных реализаций будут представлены в отдельной статье.
Автор выражает искреннюю благодарность сотрудникам научно-исследовательской лаборатории снежных лавин и селей, особенно С.С. Черноморцу,
A.Н. Божинскому, В.А. Светлосанову, Ю.Б. Андрееву,
B.Н. Кудину за ценные дискуссии.
9. Hungr O. A model for the runout analysis of rapid flow slides, debris flows and avalanches. // Canad. Geotechn. J. 1995. Vol. 32, N. 4. P. 610-623.
10. Itasca, UDEC: Universal distinct element code. URL: http://www.itascacg.com/udec/overview.php (дата последнего обращения 20.03.2011)
11. Landslide lake in Northwest Pakistan. NASA Earth observatory. URL: http://earthobservatory.nasa.gov/IOTD/view. php?id=45499 (дата последнего обращения 20.03.2011)
12. RAMMS: Rapid mass movements. URL: http://ramms. slf.ch/ramms/ (дата последнего обращения 20.03.2011)
13. Schilling S.P. LAHARZ; GIS programs for automated mapping of lahar-inundation hazard zones. U.S. Geological survey open-file Report 98-638. 1998, 80 p. URL: http:// vulcan.wr.usgs.gov/Projects/LAHARZ/framework.html (дата последнего обращения 20.03.2011)
14. SLOPE/W 2007: Slope stability analysis. URL: http:// www.w2hk.net/geoslope/slopew.html (дата последнего обращения 20.03.2011)
15. Wei F, Hu K, Lopez J.L., Cui P. Method and its application of the momentum model for debris flow risk zoning // Chin. Scien. Bull. 2003. Vol. 48. N 6. P. 594-598.
Поступила в редакцию 24.06.2010
V.O. Mikhailov
3D MATHEMATICAL MODEL OF LANDFALL PROCESSES
The article deals with the main approaches to the compilation of the DEBRIS 3D mathematical model of landfall material movement, as well as the results of modeling for the Attabad landfall episode in the Khunza River valley. The model describes the landfall movement by the trajectories of the mass centers of particular rick fragments, which are in elastic interaction with one another and with the bottom surface and undergo the influence of adhesion.
Key words: mathematical modeling, landfalls.
