CC BY
18
УДК 519.633
Д. П. Григорьевых, Н. И. Хохлов, И. Б. Петров
Московский физико-технический институт (государственный университет)
Математическое моделирование трещин в твердых деформируемых телах с использованием гексаэдральных сеток
Рассматриваются возможности моделирования трещин и неоднородностей в твердых деформируемых телах при использовании сеточно-характеристического метода на одной лагранжевой гексаэдральной сетке. Рассматриваются задачи прохождения упругих волн через постоянные флюидо- и газонасыщенные трещины и неоднородности, а также процессы образования новых трещин. За основу модели образования трещин взята широко известная модель Майчена-Сака.
Ключевые слова: моделирование, трещины, кластер трещин, сейсморазведка, гексаэдральные сетки.
D. P. Grigorievykh, N. I. Khokhlov, I. B. Petrov
Moscow Institute of Physics and Technology (State University)
Mathematical modeling of cracks in solid deformable bodies using hexahedral grids
The article discusses a possibility of modeling cracks and irregularities in solid deformable bodies using the grid characteristic method on one of Lagrangian hexahedral grids. The problems of elastic waves passing through fluid or gas-saturated cracks and irregularities are discovered as well as processes of formation of new cracks. The basis of the cracking model is the Maichiev-Saka model.
Key words: modeling, cracks, cluster of cracks, seismic survey, hexahedral grids.
1. Введение
В последнее время в задачах сейсморазведки полезных ископаемых широко используются численные методы и полноволновое моделирование. Например, в [1] решается задача прямого моделирования отражения волнового возмущения от кластера трещин, заполненных углеводородами, на тетраэдральных сетках. В [2] решается аналогичная задача, но уже на криволинейных гексаэральных сетках, что накладывает ограничения на рассматриваемую геометрию, невозможность моделирования динамического роста трещин, а также уменьшает допустимый шаг интегрирования по времени из-за сильно искривленных ячеек. В обеих работах, так же как и в данной, используется сеточно-характеристический метод [3], разработанных для численного решения задач динамического поведения деформируемых твердых тел в [4].
В данной работе предлагается алгоритм моделирования трещин на прямоугольных гек-саэдральных сетках, сочетающий в себе возможности метода на тетраэдральных сетках и очевидное преимущество гексаэральных сеток по скорости расчетов.
2. Определяющие уравнения
Бесконечно малый объем линейно-упругой среды подчиняется следующей системе дифференциальных уравнений:
= Vац, [г,] = х,у,г),
да,
г]
дек1 дЪ ,
= Х5 ^ 5к1 + К^гк + 5ц + ) ,
м
Щк1-
(1)
где р - плотность материала, Л и ц, - параметры Ламе, V^ - ковариантная производная по координате ^ - символ Кронекера, Ьг - компоненты скорости, а^ - компоненты тензора напряжений, е^ - компоненты тензора деформаций.
Продольная и поперечная скорости звука с\ и С2 выражаются через плотность и параметры Ламе следующим образом:
С1
,С2 = Л V р V р
Систему (1) можно записать в следующей канонической форме системы линейных дифференциальных уравнений в частных производных гиперболического типа:
ди
Ж
д и
дх
д и
т + Ах + Ау + я ду
д и
дх
где и = (ь х, V у, Ух, а хх, аХу, а хх, ауу, аух, агг) вектор искомых функций, х,у,х независимые переменные, их,уу, - компоненты скорости V, ахх, &ху, , &уу, , - компоненты тензора напряжений а, Ах, Ау, Аг - матрицы (9 х 9) вида:
Ая
( 0 0 0 _ 1 р 0 0 0 0 0 0
0 0 0 _ 1 р 0 0 0 0 0
0 0 0 0 _ 1 р 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
- X 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
V - X 0 0 0 0 0 0 0 0
0
А,,
0 0 _1 р 0 0 0 0 0 0 0
0 0 _1 р 0 0 0 0 0 0
0 0 0 _1 р 0 0 0 0 0
0 - X 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 - X 0 0 0 0 0 0 0
Аг
0 _ 1 р 0 0 0 0 0 0 0 0
0 _ 1 р 0 0 0 0 0 0 0
0 0 _ 1 р 0 0 0 0 0 0
0 0 —X 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 —X 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 —(X + 2^) 0 0 0 0 0 0
3. Численный метод
Переход от одного временного слоя к другому может быть записан в виде
ип+1 = ^(Ах, Ау, Аг)ип,
где Р - оператор перехода, вид которого определяет конкретную разностную схему. В работе используется следующая схема расщепления:
Р{AX, Ау, Аг) — РХ(АХ^Ру(Ау)Рг(Аг).
Таким образом, решение исходной системы сводится к последовательному решению трех одномерных систем, причем для решения следующей системы используется решение предыдущей как начальное приближение:
9их + А дих х дх
0;
9иу + . диу
М + Ау ду
0; ^ + А ди
т
дх
0,
матрицы Ах, Ау, Аг разрежены и аналитически приводятся к диагональному виду Лх — Лу — Аг — diag(cl, — с1,с2, —с2,с2, —с2, 0, 0, 0). Выпишем решение задачи для и где г = х, у или г:
А^ = .
После замены переменных — &гиг получаем 9 независимых линейных уравнений переноса:
д^ + А дг
0.
Тогда решение на (п + 1)-м временном слое будет иметь вид
и?+1 — О71^?+1.
(2)
Для решения одномерных уравнений в частных производных гиперболического типа применяется сеточно-характеристический метод [3] четвертого порядка точности (индекс г опущен для простоты записи):
ш.
п+1
— — аД — — а(Д? — аА2))),
где
а — Т,
п — 24(—+2
п '2 — 24 (—Шт+2 +
п '3 — 24 (2шт+2 —
п '4 — 24 (шт+2 — '
оп
оп
0п
гп
п
п
I - собственное значение, т - шаг интегрирования по времени, к - размер ячейки. В работе используется постоянный шаг интегрирования по времени, который выбирается достаточно малым, чтобы условие Куранта-Фридрихса-Леви ( и ^ 1) выполнялось на протяжении всего расчета. Так как в данной работе рассматриваются только небольшие деформации и искривления сеток, то было достаточно использовать шаг по времени, в е = 1.5 раза меньший, чем если бы была прямоугольная сетка. При больших деформациях необходимо использовать переменный шаг интегрирования по времени.
Для обеспечения монотонности разностной схемы используется ограничитель:
шп+1 = Гтах(<^_згдп{1)),ш^+1 > тах«,^_sгgn{V)), т \тШ«, ШПт ^а^ < тШ«, иП_8гдп(1)).
На границе области интегрирования система уравнений дополняется тремя граничными условиями:
Вип = Ь, (3)
где В - матрица 3 х 9, Ь - вектор размерности 3.
Решение (2) на верхнем временном слое будет иметь вид
ип+1 = О(гп)^п+1(гп) + ^{оиЬ)ип+1{оиЬ) = ип+1(гп) + О(ои€)^п+1(о«4), (4)
где О(гп') и - прямоугольные матрицы, составленные из тех столбцов матрицы О ,
которые соответствуют входящим или выходящим из области интегрирования характеристикам соответственно. и>п+1(1п^ находится так же, как и шп+1 для внутренних точек, неизвестен лишь ^п+1(оиЬ), который можно выразить из граничных условий (3):
^п+1(ои±) = )_1(Ь - Вип+1(*п)). (5)
Подставляя это в (4), получаем общую формулу для расчета граничных узлов:
ип+1 = ип+1(*п) + О(ои*)(ВО(ои*))_\Ъ - Вип+1(*п)). (6)
Введем вектор п - нормаль к поверхности. Оператором а ® Ъ обозначается тензорное произведение векторов, (а®Ь)у = aibj. В работе использовались 3 типа граничных условий.
1. Для заданной плотности поверхностных сил на границе f (а ■ п = f), скорость и тензор напряжения рассчитываются по формулам
^+1 = ип.+1(гп.) _ С1гп+1 - (С1 - С2)(гп+1 ■ п)п
рс 1С2 '
ап.+1 = ап.+1(гП) - гп+1 ^ п - п ^ гп.+1 -
(гп+1 ■ п)
X + 2/л
где гп+1 = ап+1(,1п) ■ п - f, I = ё1а®(1,1,1).
2. Для заданной скорости V границы:
( XI - 2(X + р)п ® п), (7)
ап+1 = сгп+1(;'-п) _
г}П+1 = у,
р[(п ■ дп+1)((с 1 - 2С2 - Сз)п ® п + С31) +
+С2(дп+1 ® П + п ® дп+1)], (8
где сэ = А/ , 9П+1 = - V.
3. Для смешанных условий (заданы только нормальная составляющая Уп и тангенциальная составляющая fT) находится полная плотность поверхностных сил f, затем используются формулы (7):
f = и + п{рс 1Уп + п ■ (ап+1('п) ■ п - рс^п+1('п))^) . (9)
Расчет контактной границы проводится следующим образом. Рассматриваются 2 типа контакта: полное слипание и свободное скольжение.
При полном слипании скорости на обеих поверхностях контактной границы совпадают при выполнении третьего закона Ньютона:
<+1 — <+1 — «™+1, <+1 • п — о-ь+1 • П, п — па — —пъ, (10)
где индексы а и Ь относятся к разным сторонам.
Окончательная расчетная формула с учетом (8) и (10) имеет вид
— -1-[Ра [(П • Уа;+1(*П) )(Са1 — Са2)п + ] +
Ра Са2 + РЬ СЬ2
+РЬ [(П • уЪа+1(га) )(сы — СЬ2)П + Сь2У'а+1{га)] —
— (аа+1(га) — фп+\(гп)) ^ п —
ъ
Ра (Са1 — Са2) + РЬ (СЬ1 — СЬ2) Ра Са1 + РЬ СЬ1
П,
(П • [РаСа1 »а+1(*а) + РЬСЬ1ЬЪ+1(Ы) — (а^+1('а) — аЪЪ+1(*а)) • п])}. (11)
Условия свободного скольжения имеют вид
п • Уа+1 — п • Vа+1 — Уп+1, аа+1 • п — <та+1 • п, [п, аа+1 • п] — п, аЪЪ+1 • п — 0, п — Па — —пъ, (12)
где [а, Ь] - векторное произведение векторов. Далее вычисляется Уп+1:
УП+1 1
Ра Са2 + РЬ) СЬ2
РаСа1 V^а+1(iа) + рЬСЬ1 — (^"О — ) • п] • П, (13)
после чего используются смешанные граничные условия (9) с fт — 0.
4. Модель бесконечно тонкой двухбереговой трещины
В данной работе используется только дискретная модель разрушений [5].
• У
• У* У У •
• •
« У У А • 1 л
• * / Л \ \ 9 Ч Ь
у Ч Л г ь
Ч г
Рис. 1. Дополнительные расчетные узлы на трещине. Красным отмечены узлы, участвующие при расчете зеленого
В каждом разрушенном узле вводится дополнительный узел таким образом, что при расчете узлов слева от трещины используется один, а справа - другой (см. рис. 1). Сами
же разрушенные узлы обсчитываются как граница сетки, по формулам (7) для газонасыщенной трещины и (12) для флюидонасыщенной.
В качестве первой тестовой задачи было рассмотрено отражение плоской продольной волны от наклонной газонасыщенной трещины (рис. 2). Как видно, расчеты обоими методами качественно совпадают, отличия наблюдаются лишь в хвосте отраженной волны.
0.012 -1-1-1-1-1-1-г
0 50 100 150 200 250 300 350 400
номер временного шага
Рис. 2. Отражение плоской продольной волны от газонасыщенной трещины. На левой половине рисунка - по описываемой модели, а на правой - криволинейной структурной сетке, на которой трещина проходит вдоль границ ячеек. На нижней половине рисунка приведен график модуля скорости от номера расчетного шага в выделенных жирным точках (зеленым - по описываемой модели, а красным - криволинейной сетке)
Как правило, внутри Земной коры трещины располагаются группами, параллельно друг другу, образуя так называемый кластер трещин. В трещинах могут быть полезные ископаемые - природный газ или нефть, которые и представляют интерес для геологов. На основании отклика от кластера трещин можно определить параметры месторождения (обратная задача). На рис. 3 и 4 приведен расчет модельной задачи для сейсморазведки с
использованием модели двухбереговых бесконечно тонких трещин. Подобный расчет был проведен в [1] с использованием треугольных сеток, результаты качественно совпадают.
Рис. 3. Прохождение плоской волны через кластер наклоненных газонасыщенных трещин (разные моменты времени)
Рис. 4. Прохождение плоской волны через кластер наклоненных флюидонасыщенных трещин (разные моменты времени)
Для моделирования динамического роста трещин используется критерий по максимальному напряжению:
тах(<7г) > а0, (14)
где - максимальное напряжение (предел прочности на разрыв). После выполнения критерия в узле образуется трещина, нормаль которой сонаправленна с максимальным напряжением, и вводится новый расчетный узел. На рис. 5 показан пример такого расчета: по массивной плите ударяет небольшой ударник, что вызывает хорошо известный «эффект откола».
5. Модель неоднородности
Стоит заметить, что на основании модели бесконечно тонкой трещины можно рассчитывать и трещины конечной толщины, а также полости (поры). Для этого достаточно сделать замкнутую трещину. Также в дополнительных узлах можно ставить условие слипания (10), а по разные стороны границы задать разные параметры среды, тогда будем иметь неоднородность. На рис. 6 приведен пример такой неоднородности.
Iжг М
Рис. 6. Прохождение плоской волны через неоднородность с меньшей скоростью звука
о о о о о о о о
О О о о < о о о о о ° о
о о о о
О ° О О о о о О Л ° ° ° о о ^ о о о
о о ° о о э о о о о о * о ° о э о о
о °о ° О - о о о о ° ° о ~ о о о
о о о о о о о 00 о О О о о ос о ° о о о О 0 О О О о
Рис. 7. Прохождение плоской волны через пористую среду. На правой половине рисунка представлен момент времени после отражения от нижней границы
В текущей программной реализации расчет идет с постоянным шагом интегрирования по времени, что накладывает ограничения на перепад параметров. Но алгоритмических препятствий для расчетов переменным и кратным шагом здесь нет, это является целью дальнейшей работы.
Отдельным интересом было моделирование пористой среды. На рис. 7 изображена модельная задача пористой среды. Плоская волна шириной порядка 10 диаметров пор проходит через материал и отражается от нижней границы. Как видно, происходит рассеивание и затухание волны с характерным интервалом затухания порядка длины волны.
6. Заключение
В данной работе были предложены довольно простые алгоритмы, позволяющие учитывать трещины и неоднородности в твердых деформируемых телах. Ключевой особенностью является работа с одной структурной сеткой, что дает большую производительность, чем при расчете на тетраэдральных сетках, а также возможность легко организовать распределенные вычисления. Проведена серия тестовых расчетов в двумерном случае, показывающая адекватное поведение данных алгоритмов. Описываемые алгоритмы могут быть обобщены и на трехмерный случай. Результаты работы могут быть применены для моделирования разрушений и прямых задач сейсморазведки.
Работа выполнена при финансовой поддержке РФФИ (проект No 15-07-01931-А).
Литература
1. Голубев В.И., Петров И.Б. Численное моделирование волновых процессов в трещиноватых средах в трехмерной постановке // Вестник Балтийского федерального университета им. И. Канта. 2014. № 4.
2. Голубев В.И., Петров И.Б., Хохлов Н.И., Шульц К И. Численный расчет волновых процессов в трещиноватых средах на гексаэдральных сетках сеточно-характеристическим методом // Ж. вычисл. матем. и матем. физ. 2015. Т. 55, № 3. C. 512-522.
3. Магомедов К.М., Холодов А.С. О построении разностных схем для уравнений гиперболического типа на основе характеристических соотношений // Ж. вычисл. матем. и матем. физ. 1969. Т. 9, № 2. С. 373—386.
4. Петров И.Б., Холодов А.С. Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом // Ж. вычисл. матем. и матем. физ. 1984. Т. 24, № 5. С. 722-739.
5. Майчен Дж., Сак С. Метод расчета «Тензор» // Вычисл. методы в гидродинамике. М.: Мир, 1967. С. 185-211.
6. Петров И.Б., Тормасов А.Г. Численное исследование косого соударения жесткого шарика с двухслойной упругопластической плитой // Матем. моделирование. 1992. Т. 4, № 3. С. 20-27.
7. Иванов В.Д., Кондауров В.И., Петров И.Б., Холодов А.С. Расчет динамического деформирования и разрушения упругопластических тел сеточно-характеристическими методами // Матем. моделирование. 1990. Т. 2, № 11. С. 10-29.
8. Холодов А.С., Холодов Я.А. О критериях монотонности разностных схем для уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. 2006. Т. 46, № 9. С. 1638-1667.
9. Кукуджанов В.Н. Вычислительная механника сплошных сред. М.: Физматлит, 2008. 320 с.
10. Бураго Н.Г., Кукуджанов В.Н. Обзор контактных алгоритмов // Известия РАН. Сер. Механика твердого тела. 2005. № 1. С. 45-87
11. Qianyi Chen, Jingtao Wang, Kaixin Liu Improved CE/SE scheme with particle level set method for numerical simulation of spall fracture due to high-velocity impact // Journal of Computational Physics. 2010. 229. 7503-7519.
12. Gagnon R.E., Wang J. Numerical simulations of a tanker collision with a bergy bit incorporating hydrodynamics, a validated ice model and damage to the vessel // Cold Regions Science and Technology, 2012.
References
1. Golubev V.I., Petrov I.B. Numerical modeling of wave processes in 3D fractured media. Herald Baltic Federal University. Kant. 2014. N 4.
2. Golubev V.I., Petrov I.B., Khokhlov N.I., Shul'ts K.I. Numerical computation of wave propagation in fractured media by applying the grid-characteristic method on hexahedral meshes. Computational Mathematics and Mathematical Physics. 2015. V. 55, I. 3. P. 509518.
3. Magomedov K.M., Kholodov A.S. The construction of difference schemes for hyperbolic equations based on characteristic relations. USSR Computational Mathematics and Mathematical Physics. 1969. V. 9, I. 2. P. 158-176.
4. Petrov I.E., Kholodov A.S. Numerical study of some dynamic problems of the mechanics of a deformable rigid body by the mesh-characteristic method. USSR Computational Mathematics and Mathematical Physics. 1984. V. 24, I. 3. P. 61-73.
5. Machen J., Sachs S. The method of calculating «Tensor». Fundamental methods in hydrodynamics. 1967. P. 185-211.
6. Petrov I.B., Tormasov I.G. Numerical investigation of oblique collision of hard sphere with two-layer elastic-plastic slab. Matem. Mod. 1992. 4:3. 20-27.
7. Ivanov V.D., Kondaurov V.I., Petrov I.B., Kholodov A.S. Calculation of dynamic deformation and distructure of elastic-plastic body by grid-characteristic methods. Matem. Mod. 1990. 2:11. 10-29.
8. Kholodov A.S., Kholodov Ya.A. Monotonicity criteria for difference schemes designed for hyperbolic equations. Computational Mathematics and Mathematical Physics. 2006. V. 46, I. 9. P. 1560-1588.
9. Kukudzhanov V.N. Computational continuum mechanics. M.: Fizmatlit, 2008. 320 c.
10. Burago N.G., Kukudzhanov V.N. Overview of contact algorithms. Izvestiya RAN. Ser. Solid Mechanics. 2005. N 1. P. 45-87.
11. Qianyi Chen, Jingtao Wang, Kaixin Liu Improved CE/SE scheme with particle level set method for numerical simulation of spall fracture due to high-velocity impact. Journal of Computational Physics. 2010. 229. 7503-7519.
12. Gagnon R.E., Wang J. Numerical simulations of a tanker collision with a bergy bit incorporating hydrodynamics, a validated ice model and damage to the vessel. Cold Regions Science and Technology, 2012.
Поступила в редакцию 22.10.2015.
