CC BY
42
Некоторые способы построения расчетных областей для компьютерного моделирования в прямых задачах
М.М. Немирович-Данченко
Институт геофизики СО РАН, Новосибирск, 630090, Россия
В статье предложено несколько способов построения моделей сред со сложной геометрией для решения двумерных и трехмерных задач сейсморазведки, сейсмологии и геомеханики. Особое внимание уделяется приемам обработки растровых изображений реальных сред и создания на этой основе исходных моделей для дальнейших расчетов. Выполнен расчет напряженного состояния в теле гравитационного оползня, проведено сравнение с результатами физического моделирования.
1. Введение
Появление персональных компьютеров с большой вычислительной мощностью привело в последние годы к резкому увеличению использования численных методов при решении самых разнообразных задач. Однако большинство исследователей, обращаясь к численному моделированию таких сложных природных и техногенных процессов, как оползневые явления, взрывы, горные удары упрощают модель среды, сводя ее к классическим прямоугольным, сферическим или цилиндрическим областям. Тем самым снижается адекватность проводимых расчетов и ухудшается возможность проведения сравнения результатов численного моделирования с натурным экспериментом, физическим моделированием, полевым наблюдением. В предлагаемой статье автор применительно к конечно-разностным методам рассматривает способы построения расчетной сетки.
В первом способе вначале строится элементарная расчетная сетка из прямоугольных ячеек (кубов). Затем численно решается задача об упругом деформировании такой среды до требуемой конфигурации. Это возможно сделать, описывая поведение среды как гипоупругое [1], поскольку при этом допускаются сколь угодно большие деформации. При достижении требуемой конфигурации все компоненты тензоров напряжений и деформаций, а также массовые скорости приравниваются нулю. Полученную область можно считать начальной недеформи-рованной средой и ставить для нее краевые задачи.
Второй способ — геометрическое построение расчетной сетки для исходной модели среды, то есть вы-
числение координат узлов сетки с использованием аналитической или числовой информации о структуре среды.
И в том и в другом случаях необходимо иметь либо точную формулу, описывающую границы расчетной области и границы внутри нее, либо иметь значения координат такого числа точек этих границ, которое достаточно для разбиения области на расчетные ячейки с учетом характерных размеров задачи. Поэтому в статье описан также способ обработки растровых изображений, позволяющий получать координаты точек границ.
2. Формирование исходных расчетных областей
Один из способов построения расчетной области с криволинейными границами, предлагаемый в настоящей статье, — это деформирование исходной классической расчетной области до нужной формы. Применяя инкрементальный закон поведения, можно, не вводя специально условие пластичности, добиться практически любой формы, требуемой для проведения определенного численного исследования. Это обусловлено тем, что среда, описываемая инкрементальным законом поведения, как в случае гипоупругости, ведет себя как идеальная резина, т.е. обратимо упруго деформируется до любых значений компонент тензора деформации.
Рассмотрим пример построения сложной расчетной области.
Исходная расчетная область — прямоугольный параллелепипед размерами 100х 100х 10 м3, всего расчетных ячеек — 100000 (рис. 1). Материал среды упругий,
© Немирович-Данченко М.М., 2004
скорость продольных волн 3 000 м/с, скорость поперечных волн 1 500 м/с, плотность 2 000 кг/м3. Одна из боковых граней АВ закреплена, на другой боковой грани CD задается постоянная по времени скорость перемещения 20 м/с.
В направлении от приложенного усилия к закрепленной грани в параллелепипеде начинают распространяться изгибные колебания. На рис. 2 для грани АС приводится векторное поле скоростей смещений для различных моментов времени. Хорошо видны невозмущенная область с левого торца, перемещающиеся вихревые области и векторы скоростей смещений в направлении изгиба.
Волновые поля на рис. 2 позволяют грубо оценить скорости основных процессов. Вихревые области распространяются со скоростью около 300 м/с, экстремумы изгибных колебаний распространяются со скоростью около 900 м/с.
Через 4 секунды процесс деформирования прерывается. Полученная область приведена на рис. 3. Она может быть использована для дальнейших конкретных расчетов. Например, подобный вид имеют многие склоны, в них может быть рассчитано формирование поля напряжений под действием силы тяжести, распространение сейсмических волн и тому подобное. Нужно только иметь в виду, что процесс деформирования первоначально прямоугольного параллелепипеда прерывается нами, поскольку область приняла некоторый удобный для каких-то дальнейших действий вид, но не потому, что изгибные волны затухли. Поэтому использовать можно только полученные эйлеровы координаты лагранжевых расчетных ячеек. Все остальные величины, если продолжать расчеты, используя одну и ту же программу, должны быть обнулены.
Таким образом, для среды, приведенной на рис. 3, формулируется обычная краевая задача с нулевыми начальными напряжениями и деформациями. Следуя [2], можно сказать, что с полученной средой можно работать в терминах малых приращений напряжений и деформаций, пренебрегая их конечностью в начальный момент времени.
В механике, геологии, геомеханике имеется много объектов с особыми точками (линиями). Это надрезы и надпилы, трещины и расселины, дизъюнктивные нарушения. В этих случаях может быть применена описан-
ная выше процедура, но численный метод, на основе которого выполняется деформирование, должен быть дополнен процедурой расчетверения точек расчетной сетки, обеспечивающей хрупкое разрушение первоначально сплошной среды [1]. Проиллюстрируем применение этого подхода для двумерного случая.
Исходная среда показана на рис. 4. Левый торец модели закреплен, к правому приложено усилие вертикально вниз. В середине модели сверху сделан надрез — задана бесконечно тонкая щель, берега ее свободны от напряжений.
На рис. 5 показан результат создания деформированной области с явно выраженной трещиной. В тех случаях, когда границы исследуемой области могут быть выражены аналитически, можно достаточно просто построить криволинейную в общем случае сетку. Но и в этом случае, как и в результате применения описанных процедур деформирования, все расчетные ячейки будут иметь разные площади, шаги по пространству, используемые при аппроксимации производных, будут также различны. Поэтому разностные схемы для проведения численного моделирования на таких сетках должны корректно описывать пространственные производные. Нами используется некоторый минимальный объем, охватывающий расчетную пространственную точку. Для определения пространственных производных применяется суммирование нормальной составляющей потока по поверхности этого объема. Как и в основной численной схеме Неймана-Рихтмайера, используются два набора сеток. Градиенты тензора напряжений центрируются в узлах сетки, а дивергенция скорости относится к мыс-
Рис. 2. Распространение изгибных колебаний в деформируемой прямоугольной призме: t = 0.233 (а) и 0.266 с (б)
Рис. 1. Исходная расчетная область Рис. 3. Готовая для дальнейшего моделирования расчетная область
/
/
/
Рис. 4. Исходная конфигурация. Стрелкой показано направление деформирования
Рис. 5. Деформированная расчетная область
ленным центрам трехмерной ячейки. Такой подход позволяет корректно аппроксимировать пространственные производные с учетом конечного, и зачастую довольно большого, деформирования расчетной сетки.
3. Об обработке растровых изображений
Очень часто вычислитель сталкивается с ситуацией, когда аналитической зависимости, описывающей внутренние и внешние границы области, не существует, табличное ее задание также неизвестно, но есть изображение этой области. Это могут быть изображения геологических разрезов, сейсмические временные разрезы, результаты геофизического профилирования, фотографические изображения образцов, склонов, срезов и т.п. Все это будем называть растровыми изображениями, поскольку после сканирования таких изображений они становятся растровой (точечной) графикой. Возможность использования информации в растровом виде (в отличие от векторного, т.е. цифрового) находит в последнее время сторонников в среде специалистов по обработке геофизических данных [3-5].
Итак, прежде всего имеющееся изображение разреза необходимо отсканировать с достаточным разрешением (нами обычно используется до 200 dpi). При высоком разрешении файл, полученный в результате сканирования, часто оказывается весьма большим (порядка 200 Мбайт). Стандартные графические пакеты, как правило, либо не могут читать файлы такого размера, либо, прочитав, не допускают возможности редактировать материал. Поэтому важным элементом обработки является создание специальных процедур физического чтения-записи, более здесь не обсуждаемых.
Обработка полученных файлов может проводиться с различными целями, например, для улучшения видимого качества изображения (очистка сгибов и пятен, перенесенных при сканировании с бумаги в растровый
образ), для оцифровки выбранного горизонта (границы) и т.д. При этом, как правило, любая обработка информации основана на фильтрации или фильтрация предшествует обработке. Для того чтобы осуществить фильтрацию растрового изображения, нами вводится основное для работы с растрами понятие амплитуды в окрестности точки растрового изображения.
Расчет значения амплитуды в окрестностях точки ^, j) растрового файла производится по формуле:
г+пх/ 2 у+пу/2
/а у) = % % ф(1 , т), (1)
I=г - пх/ 2 т=у - пу /2
где ф(I, т) — значение, возвращаемое функцией Getpixel(l, т); пх и пу — размеры окна, в котором производится суммирование. Изменение размеров окна приведет, вообще говоря, к изменению значения амплитуды. Легко видеть, что введенная таким образом функция осуществляет фильтрацию изображения (сглаживающее интегрирование с переменным окном). Операция (1) некоммутативна, единичная функция из {/} — это операция (1) при пх = пу = 0.
Поясним суть операции (1) на простейшем примере. Исходное растровое изображение задано на рис. 6, а. Минимальный элемент растра — пиксел. Каждый пиксел изображаем квадратом, имеющим свой цвет. Условимся, что черный цвет имеет значение 0, а белый — 15 (в 16-цветной палитре). Эти значения и возвращаются при обращении к функции Getpixel(i, j).
На рис. 6, а приведены две «границы»: граница 1 — четкая, толщиной в один пиксел и граница 2 — размытая, ее координаты сразу определить нельзя (речь ведь идет не о визуальном определении положения этой границы, а об автоматическом построении разностной сетки для любых границ).
На рис. 6, б показано применение формулы (1). В области границы 1 выберем сканирующее окно размером 1x3 пиксела. По смыслу формулы (1) точка (^ ]') — это точка в центре окна. Nx и Ыу в нашем случае равны
а б
ш
1 ■ ш і ■■■ ■■
—і ■ ■ ■ ■ ■ ■ ■ ■ ■■■■■■■■
— ■■■■
■■■
■■■
2■■■■■■■■■ 2■ ■■
■ ■■■
■■■
Рис. 6. Пример работы оператора фильтрации растрового изображения
1 и 3 соответственно; Nyj 2 = 1, поэтому для каждой текущей (i, 7*)-й точки будем суммировать значения цвета в самой этой точке и в точках выше и ниже текущей. В зависимости от условия, накладываемого на результат суммирования, в центр окна, т.е. в текущую точку, будет записано какое-то число (номер цвета). Применительно к границе 1 потребуем, чтобы результат суммирования по формуле (1) был менее 45. Например, если в окне суммирования оказался хотя бы один черный пиксел, то сумма из трех значений цвета будет равна 30. При таком результате в точку (i, j) записывается пиксел черного цвета, в противном случае — белого. В итоге получаем равномерно размытую границу 1 (рис. 6, б). Размывание границ контактов бывает полезно при описании специального поведения среды на контакте. Граница 2 на рис. 6, б была обработана тем же окном 1x3, но с использованием условия, что результат суммирования должен быть равен нулю (все три пиксела черные), тогда в центр окна суммирования ставится черный пиксел. При этом граница стала четкой, толщиной в один пиксел. Учитывая, что каждый пиксел в растровом изображении имеет свою координату, определенную для формата Bitmap единственным способом, мы получаем границу, заданную целочисленными значениями абсцисс и ординат.
4. Расчет напряженного состояния в теле гравитационного оползня
Для проверки работы численной схемы с непрямоугольными ячейками была решена модельная задача о развитии напряжений в теле неустойчивого оползня.
За основу взяты данные из монографии [6]. В ней подробно рассмотрены различные виды оползней, характерные для Кыргызстана. Для анализа напряжений в теле оползня и оценки его устойчивости в [6] приводятся результаты физического моделирования для моделей, приведенных на рис. 7. Рассмотрено три основных угла склона — 15, 25 и 35 градусов. Методом фотоупругости построены области максимальных касательных напряжений. Показывается, что для углов склона 35 градусов наиболее характерным является смещение области концентрации напряжений от угловой точки вниз по склону (рис. 7, а) и это служит признаком потенциальной неустойчивости склона и его оползнеопаснос-ти.
Для численного моделирования была взята область с углом наклона 35 градусов (рис. 8). Мощность покровных отложений (толщина) на склоне составила 24.5 метра. Упругие свойства среды: модуль Юнга E = 5 МПа, плотность р = 1500 кг/м3, коэффициент Пуассона V = 0.32. Сила тяжести входит в уравнения движения и направлена вертикально вниз. Расчет проводится в рамках модели гипоупругой среды [1], уравнения здесь не приводятся.
Рис. 7. Картина напряженного состояния покровных отложений на горных склонах: а — крутизна склона 35°; б — крутизна склона 25°; в — крутизна склона 15°; 1, 2, 3, 4 — порядок полос “п”; 0 — изотропная область [6]
С
Рис. S. Исходная геометрия для расчета оползневого течения
Рис. 9. Изолинии максимальных касательных напряжений в верхней части расчетной области
Рис. 10. Напряжения в теле неустойчивого оползня
Для получения развитой картины оползневого течения при моделировании сцепление тела оползня с нижележащими породами было принято равным нулю (условие скольжения). На рис. 10 приводится вся расчетная область с теневыми картинами максимальных касательных напряжений.
На рис. 10 хорошо видно, во-первых, описанное в [6] проседание среды вблизи точки D (рис. 8). Кроме того, область концентрации растягивающих напряжений спустилась к нижнему правому концу модели, что соответствует подножию склона. Отметим еще, что максимальное значение напряжений достигается на нижней границе склона и равно 0.45 МПа, на дневной поверхности максимальное значение составляет 0.12 МПа.
5. Заключение
В работе приводится техника построения неклассических расчетных областей для конечно-разностного моделирования в прямых задачах. Рассмотрен один прием обработки растровых изображений, могущий оказаться полезным при построении моделей реальных сред, в частности, в задачах геофизики и геомеханики. Работа численной схемы с непрямоугольными ячейками демонстрируется на задаче о напряженном состоянии в теле гравитационного оползня. Результаты расчета качественно достаточно хорошо соответствуют опубликованным опытным данным и данным натурных наблюдений.
Работа выполнена при финансовой поддержке СО РАН и НАН Кыргызстана (интеграционный проект № 172), СО РАН (интеграционный проект № 129, 2003 г.) и гранта Президента РФ по поддержке ведущих научных школ № НШ-1302.2003.5 (школа академика С.В. Гольдина).
Общая длина модели 250 м, длина АВ = 70 м, толщина BD = 24.5 м, число расчетных ячеек 250x31.
Такое сочетание констант среды — модуль Е и коэффициент Пуассона V — определено использованной литературой. При численном моделировании нами используются скорость продольной волны Ур, скорость поперечной волны У и плотность р. Приведем для удобства формулы, связывающие разные системы констант:
Е Е Х уЕ
К= —:гт> Ц = ^У = -
3(1 - 2у) 2(1 + у) (1 + у)(1 - 2у)
Х + 2ц = 3К Ур = , V, = Ц
(1+V) "у Р • • Ур
Здесь К — модуль объемного сжатия; X и ц — константы Ламэ.
Необходимо также отметить, что при расчетах прямых задач с характерным размером модели порядка десятков и сотен метров удобно пользоваться следующей согласованной системой единиц. Размерность шага по пространству — [м], по времени — [мс], плотность — [кг/м3]. Тогда рассчитываемые напряжения имеют размерность [МПа].
На рис. 9 приводится установившаяся теневая картина максимальных касательных напряжений в верхней части модели. Хорошо видно, что особая точка поля напряжений не совпадает с углом модели, а сдвинута вниз по склону, что соответствует результатам физического моделирования [6].
В [6] приводится качественное описание развития неустойчивости оползня. При этом отмечается, что при развитом процессе наблюдается проседание дневной поверхности на расстоянии от бровки порядка толщины слоя. Затем зона концентрации касательных напряжений (зона растяжения) достигает подножия склона.
Литература
1. Немирович-Данченко М.М. Модель гипоупругой хрупкой среды: применение к расчету деформирования и разрушения горных пород // Физ. мезомех. - 1998. - Т. 1. - № 2. - С. 107-114.
2. Жермен П. Курс механики сплошных сред. Общая теория. - М.: Высшая школа, 1983. - 399 с.
3. Montgomery S.L. Raster logs may be basic for a geological workstation
// Oil&Gas J. - 1997. - No. 14. - P. 84-88.
4. Balch A.H. Color sonagrams — a new dimension in seismic data interpretation // Geophysics. - 1971. - V. 36. - No. 6. - P. 1074-1098.
5. Verm R., Hilterman F. Lithology color-coded seismic sections: the calibration of AVO crossplotting to rock properties // The Leading Edge. -1995. - V. 14. - No. 8. - P. 847-853.
6. Айтматов И.Т., Кожогулов К.Ч., Никольская O.B. Геомеханика оползнеопасных склонов. - Бишкек: Изд-во «Илим», 1999. - 208 с.
On ways of generation of calculation areas for computer simulation of direct problems
M.M. Nemirovich-Danchenko
Institute of Geophysics SB RAS, Novosibirsk, 630090, Russia
The paper proposes a number of ways to construct models of media with complex geometry for solving 2D- and 3D-problems of seismic prospecting, seismology and geomechanics. Particular attention is given to the methods of processing of scans of real media and to the generation of initial models on this basis. Stress state of a gravity landslide is calculated and the results obtained are compared to those of physical simulation.
