Математическое моделирование разрушений в твердых деформируемых телах сеточно-характеристическим методом (1187403)
Текст из файла
Московский физико-технический институт(государственный университет)Факультет управления и прикладной математикиКафедра информатикиМатематическое моделирование разрушений втвердых деформируемых телахсеточно-характеристическим методомДипломная работа на степень магистрастудента 973 группыГригорьевых Д.П.Научный руководитель:к.ф.-м.н. Хохлов Н.И.Работа выполнена в МФТИ,г.Долгопрудный,2015г.Содержание1 Введение32 Определяющие уравнения43 Обзор методов моделирования разрушений64 Програмный комплекс rect75 Сеточно-характеристический метод на прямоугольных сетках85.1 Покоординатное расщепление .
. . . . . . . . . . . . . . . .85.2 Решение одномерных систем . . . . . . . . . . . . . . . . .95.3 Расчетная сетка. . . . . . . . . . . . . . . . . . . . . . . .105.4 Граничные корректоры . . . . . . . . . . . . . . . . . . . .115.5 Контактные корректоры . . . . . . . . .
. . . . . . . . . . .126 Модели разрушения146.1 Модель сдвигового разрушения . . . . . . . . . . . . . . . .156.1.1Сейсмостойкость купола АЭС . . . . . . . . . . . . .166.1.2Сейсмостойкость двухэтажного жилого дома . . . .176.1.3Взрыв круга из центра . . . . . . . . . . . . . . . . .176.2 Модель однобереговых трещин . . . . . . . . . . .
. . . . .186.2.1Соударение стальной и стеклянных пластин . . . .6.2.2Соударение стального шарика и стеклянной пластины 216.2.3Соударение двух шариков о стеклянную пластину .236.2.4Рост трещин при растяжение пластины . . . . . . .236.2.5Отражение плоской волны от однобереговой трещины 241206.3 Модель двухбереговых трещин . .
. . . . . . . . . . . . . .6.3.1Кластер двухбереговых трещин . . . . . . . . . . . .6.3.2Прохождение плоской волны через флюидонасыщенные поры . . . . . . . . . . . . . . . . . . . . . . . .7 Заключение2526282921ВведениеКомпьютерное моделирование широко применяется в настоящее вре-мя практически во всех областях: строительстве, георазведке, медицинеи прочее. Подходов к моделированию тоже множество, от лучевых методов, до численного решения системы определяющих уравнений. Моделирование разрушений одна из основных задач в этой области.В данной работе рассматриваются задачи дискретного динамического разрушения твердых деформируемых тел.
Приведен краткий обзорсуществующих алгоритмов разрушения, описание програмного комплеса, численного метода и результаты расчетов. Для учета разрушений использовались модели сдвигового разрушения и образования трещин, заоснову которой взята модель Майчена-Сака [1]. Для численного решения рассматриваемых задач используется сеточно-характеристическийметод [2], разработанный для численного решения задач динамическогоповедения деформируемых твердых тел в [3].В данной работе рассматривались задачи о соударении стальных шариков со стеклянной пластиной с откольными разрушениями, сейсмостойкости наземных сооружений, отражение плоской волны от кластератрещин и другие.32Определяющие уравненияБесконечно малый объем линейно-упругой среды подчиняется следу-ющей системе дифференциальных уравнений:∂vi= ∇j σij , (i, j = x, y, z),∂t∂σij∂kl= qijkl+ Fij ,∂t∂t= λδij δkl + µ(δik δjl + δil + δjk ),ρqijkl(1)где ρ - плотность материала, λ и µ - параметры Ламе, ∇j - ковариантнаяпроизводная по координате j, δij - символ Кронекера, vi - компонентыскорости, σij - компоненты тензора напряжений, ij - компоненты тензорадеформаций, Fij - не обязательная правая часть.Продольная и поперечная скорости звука c1 и c2 и модуль всестороннего сжатия K выражаются через плотность и параметры Ламе следующим образом:sλ + 2µ,ρrµc2 =,ρc1 =2K = λ + µ.3(2)(3)(4)Производную тензора напряжений можно приближенно выразить через скорость:ij1≈ (∇i vj + ∇j vi ).∂t2(5)Тогда систему (1) можно переписать в следующей форме:∂u∂u∂u∂u+ Ax+ Ay+ Az= 0,∂t∂x∂y∂z4(6)где u = (vx , vy , vz , σxx , σxy , σxz , σyy , σyz , σzz )T - вектор искомых функций,x, y, z - независимые переменные, vx , vy , vz - компоненты скорости v, σxx ,σxy , σxz , σyy , σyz , σzz - компоненты тензора напряжений σ, Ax , Ay , Az матрицы (9 × 9) вида:000000000−(λ + 2µ) 00Ax = 0−µ 000 −µ−λ00000−λ0000− ρ1 000 000 0−λ0Ay = 00−µ 000 0 −(λ + 2µ) 0 00−µ0−λ05− ρ1000 00− ρ100 0000000 00000 00000 00000 00000 0000− ρ1 0 0000− ρ1000− ρ1 000000,0000 0 00 0 00 0 00 0 00 0 00 0 0,0 0 00 0 00 0 0000000000000000000 0 0 0(7)(8)− ρ1000 00− ρ10 000− ρ1 00−λ0Az = 000 0−µ 000 00−λ0 0 −µ0000 −(λ + 2µ) 030 0 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00000.0000(9)Обзор методов моделирования разрушенийМетод частиц.
Подробно метод описан в [4]. Среда моделирует-ся множеством частиц, связанных друг с другом потенциальным поле.Разрушение происходит когда расстояние между частицами превышаетпотенциал взаимодействия, таким образом это естественный для моделипроцесс. Интересные результаты по моделированию откольных разрушений и пробоя тонкой пластины приведены в [5]. Описание континуальноймодели дано в [6]. Метод позволяет легко учитывать разлет осколков инаиболее подходит для расчетов больших деформаций и пробоев.
Из минусов - низкий порядок. Большинство работ по разрушению выполняютименно этим методом.Тетраэдральные сетки. Данный метод хорошо подходит для моделиривания стационарных трещин, как в [7] и [8], за счет использованияиерархических сеток. Динамический рост разрушения имеет схожие проблемы и идеи, как и при расчетах на прямоугольных сетках, описанныев данной работе.64Програмный комплекс rectНа кафедре информатики разрабатывается програмный комплекс длярешения уравнений упругости на прямоугольных сетках. Комплекс состоит из множества модулей и имеет легко расширяемую архитектуру.Далее список основных элементов:1.
grid – Описание сетки, может наследовать параметры от другой сетки.2. grid.schema – Используемая одномерная схема на 5-ти точечномшаблоне. Реализовано несколько десятков различных схем.3. grid.f iller – Функция, применяемая к сетке перед переходом на новый временной слой. Например, используется для копирования граничных узлов на 2 мнимых слоя при создании условия поглощения.4. grid.corrector – Функция, применяемая к сетке после перехода нановый временной слой.
Это наиболее широкая часть. Здесь реализован учет правой части уравнений (например, ускорение свободногопадения), сила на границе и, собственно, реализованный в даннойработе учет разрушения.5. contact – Функция, применяемая одновременно к двум сеткам. Реализованы условия слипания и скольжения.6. initial – Начальные условия. Плоская волна, взрыв, землетрясениеи другие.7. saver – Как сохранять результаты расчетов (vtk, cvs).7Комплекс может работать с произвольным количеством сеток, но наданный момент расчет в несколько потоков с использованием MPI возможен только на одной сетке.
Поэтому отдельный интерес представлетисследование алгоритмов, использующих всего одну прямоугольную сетку.Скорость расчета получается в десятки раз выше, чем в программах, работающих на тетраэдральных сетках, но проблематично считатьсложные структуры.5Сеточно-характеристический метод на прямоугольных сетках5.1Покоординатное расщеплениеПереход от одного временного слоя к другому может быть записан ввиде:un+1 = F (Ax , Ay , Az )un ,(10)где F - оператор перехода, определяющий используемую разностную схему.В работе используется слудующая схема расщепления:F (Ax , Ay , Az ) = Fx (Ax )Fy (Ay )Fz (Az ).(11)В результате решение исходной системы сводится к последовательному решению трех одномерных систем:∂u∂u∂u∂u∂u∂u+ Ax= 0;+ Ay= 0;+ Az= 0.∂t∂x∂t∂y∂t∂z8(12)5.2Решение одномерных системA разрежены и аналитически приводятся к диагональному виду Λ =diag(c1 , −c1 , c2 , −c2 , c2 , −c2 , 0, 0, 0).
Далее индекс у A опускается, подразумевая, что мы решаем последовательно задачу для x, y и z.A = Ω−1 ΛΩ.(13)После замены переменных ω = Ωu получаем 9 независимых линейныхуравнений переноса.∂ωω+Λ= 0,∂t∂α(14)где α - соответствующая ось координат (x, y или z). Далее восстанавливается решение на (n + 1)-ом временном слое:un+1 = Ω−1 ω n+1 .(15)Для решения одномерных уравнений в частных производных гиперболического типа применяется сеточно-характеристический метод 4-огопорядка точности:n+1nωm= ωm− σ(∆n1 − σ(∆n2 − σ(∆n3 − σ∆n4 ))),гдеσ = λτh,1nnnn∆n1 = 24(−2ωm+2+ 16ωm+1− 16ωm−1+ 2ωm−2),1nnnnn(−ωm+2+ 16ωm+1− 30ωm+ 16ωm−1− ωm−2),∆n2 = 241nnnn∆n3 = 24(2ωm+2− 4ωm+1+ 4ωm− 2ωm−2),nnnnn∆n4 = 1 (ωm+2− 4ωm+1+ 6ωm− 4ωm−1+ ωm−2).249(16)Для обеспечения монотонности разностной схемы используется ограничитель:n+1ωm=n+1nnmax(ω n , ω n> max(ωm, ωm−sign(λ))mm−sign(λ) ), ωmn+1nnmin(ω n , ω n< min(ωm, ωm−sign(λ))mm−sign(λ) ), ωm5.3Расчетная сеткаВ расчетах используется гексаэдральная структурная лагранжеваярасчетная сетка.
В общем случае область интегрирования разбиваетсяна подобласти, каждая из которых имеет свою регулярную ijk-сетку (см.рис. 1). На каждом шаге по времени криволинейная сетка отображаетсяна прямоугольную.Рис. 1: 5 расчетных сеток для круга.Корректировка координат узлов сетки на каждом шаге по времени,для каждого узла осуществляется следующим образом:r n+1 = r n + v n τ,(17)где r n - координаты узла, v n - скорость среды, τ - шаг интегрированияпо времени.105.4Граничные корректорыНа границе области интегрирования система уравнений дополняетсятремя граничными условиями:Bun = b,(18)где B - матрица 3 × 9, b - вектор размерности 3.Решение (15) на верхнем временном слое будет иметь вид:un+1 = Ω(in) ω n+1(in) + Ω(out) ω n+1(out) = un+1(in) + Ω(out) ω n+1(out) ,(19)где Ω(in) и Ω(out) - прямоугольные матрицы, составленные из тех столбцов матрицы Ω−1 , которые соответствуют входящим или выходящимиз области интегрирования характеристикам, соответственно.
ω n+1(in)находится так же, как и ω n+1 для внутренних точек, неизвестен лишьω n+1(out) , который можно выразить из граничных условий (18):ω n+1(out) = (BΩ(out) )−1 (b − Bun+1(in) ).(20)Подставляя это в (19) получаем общую формулу для расчета граничныхузлов:un+1 = un+1(in) + Ω(out) (BΩ(out) )−1 (b − Bun+1(in) ).(21)Введем вектор n - нормаль к поверхности. Оператором a ⊗ b обозначается тензорное произведение векторов, (a ⊗ b)ij = ai bj . В работеиспользовались 3 типа граничных условий.1. Для заданной плотности поверхностных сил на границе f (σn =11f ), скорость и тензор напряжения расчитывается по формулам:c1 z n+1 − (c1 − c2 )(z n+1 · n)n,ρc1 c2= σ n+1(in) − z n+1 ⊗ n − n ⊗ z n+1 −(z n+1 · n)−(λI − 2(λ + µ)n ⊗ n),λ + 2µun+1 = un+1(in) −σ n+1(22)где z n+1 = σ n+1(in) n − f , I = diag(1, 1, 1).2. Для заданной скорости V границы:v n+1 = V ,σ n+1 = σ n+1(in) − ρ[(n · g n+1 )((c1 − 2c2 − c3 )n ⊗ n + c3 I) ++c2 (g n+1 ⊗ n + n ⊗ g n+1 )],где c3 =qλ2ρ(λ+2µ) ,(23)g n+1 = v n+1(in) − V .3.
Для смешанных условий (заданы только нормальная состовляющая Vn и тангенциальная составляющая fτ ) находится полная плотностьповерхностных сил f , затем используются формулы (22):n+1(in)n+1(in)f = fτ + n ρc1 Vn + n · (σn − ρc1 v) .5.5(24)Контактные корректорыРасчет контактной границы проводится следующим образом. Рассматриваются 2 типа контакта: полное слипание и свободное скольжение.При полном слипании скорости на обеих поверхностях контакта совпадают при выполняении третьего закона Ньютона:van+1 = vbn+1 = v n+1 , σan+1 na = −σbn+1 nb , na = −nbгде индексы a и b относятся к разным сторонам контакта.12(25)Необходимо так подобрать скорость контакта v n+1 , чтобы на контактной границе выполнялись условия (23) и (25). Окончательная расчетная формула имеет вид:v n+1 =1{ρa [(n · van+1(in) )(ca1 − ca2 )n + ca2 van+1(in) ] +ρa ca2 + ρb cb2n+1(in)+ρb [(n · vbn+1(in))(cb1 − cb2 )n + cb2 vb]−n+1(in)−(σan+1(in) − σb)n −ρa (ca1 − ca2 ) + ρb (cb1 − cb2 )−nρa ca1 + ρb cb1n+1(in)(n · [ρa ca1 van+1(in) + ρb cb1 vbn+1(in)− (σan+1(in) − σb)n])}.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.