CC BY
2Электронный журнал «Труды МАИ». Выпуск № 40
www.mai.ru/science/trudy/
УДК 539.3
Применение метода ГИУ для решения краевых задач трехмерных динамических теорий вязко- и пороупругости
А. А. Белов, Л.А. Игумнов, И.С. Карелин, С.Ю. Литвинчук
Аннотация
Представлены результаты расчетов динамического состояния вязкоупругих тел на основе метода граничных элементов (МГЭ) в сочетании с методами квадратур сверток и Дурбина. Решение вязкоупругих краевых задач в трехмерной постановке в явном времени, организовано без использования шаговых процедур. Рассматривается модель пористой среды с двухфазной внутренней структурой, предложенная Био. Представлены результаты расчетов динамического состояния конечных пороупругих тел на основе МГЭ. Представлены численные эксперименты. Ключевые слова
граничные интегральные уравнения; граничный элемент; вязкоупругость; пороупругость.
Введение
Нестационарные динамические задачи теории вязкоупругости решаются методом граничных интегральных уравнений (ГИУ). Представлена методика численного решения систем граничных интегральных уравнений прямого подхода в сочетании с преобразованием Лапласа и методом Дурбина для расчета неизвестных волновых полей трехмерных вязкоупругих составных тел, и представлен подход метода граничных элементов с явным учетом переменной времени: использована гранично-элементная техника построения дискретного аналога в сочетании с методом квадратур сверток. Вязкоупругие свойства
материалов описываются соотношениями классических регулярных моделей (Максвелла, Фойгта, стандартного вязкоупругого тела) или слабосингулярной степенной моделью.
Исследование волновых процессов в пороупругих телах представляет научный и практический интерес. Для широкого диапазона насыщенных материалов упругая и вязкоупругая теории являются грубым приближением при исследовании распространения волн. Для учета пористости использована теория М. Био [1, 2]. Вопросами распространения волн в пористых насыщенных средах в последние годы занимались Н.С. Городецкая (1998), Н. Дунин, Д. Михайлов, В. Николаевский (2002), R. Ababou и др. (2002), D.F. Aldridg и др. (2005), J. Jocker и D. Smenlders (2005), G. Chao и др. (2005), H.F. Wang (2000) и многие другие. Вопросы, связанные с построением интегральных представлений решений в полной строгой математической постановке, являются актуальными. Несмотря на то, что в литературе представлены варианты сингулярных граничных интегральных уравнений [3, 4], имеются лишь единичные гранично-элементные решения краевых динамических задач пороупругости. В работе используются модифицированные интегральные представления волновых полей в пороупругих средах, полученные ранее M. Schanz [3]. На основе новых ГИУ пороупругости получены численные гранично-элементные решения прямых трехмерных динамических задач.
1. Постановка задачи
Рассматривается кусочно-однородное тело О в трехмерном евклидовом пространстве R3 с декартовой системой координат Ox1 x2 x3. Границу тела обозначим через Г, границы однородных частей Оk (k = 1,...,K) - через Г. Предполагается, что Оk являются изотропными вязкоупругими телами [5, 6]. Введем следующие обозначения для параметров материала каждой однородной части (подконструкции) Оk : рk - плотность материала, Xk (t) и цk (t ) — функции Ламе материала. Динамическое состояние каждой части тела О k описывается следующей системой дифференциальных уравнений в перемещениях:
перемещений точки х = (х1 х2х3) в момент времени Физические и геометрические соотношения имеют вид:
цk (t ) * Auk ( x, t ) + (Xk (t ) + цk (t )) * grad div uk ( x, t ) = рk uk ( x, t ),
(1)
где символ
» * и
означает свертку Стилтьеса по времени t. В уравнениях (1) uk ( x, t ) - вектор
где g , s , i = 1,3, j = 1,3 - тензоры напряжений и деформаций.
Пусть вектор перемещений и функции Ламе материала удовлетворяют условиям: uk (х,0) = uк (х,0) = 0, ^(t - т) = 0, 1(t - т) = 0, где t < т , limt^0A(t) = 0, limt1 (t) = 0.
Конкретный вид функций p,(t) и 1(t) определяется вязкоупругой моделью материала. Будем рассматривать случай пропорциональных функций памяти, тогда достаточно описать физические соотношения, к примеру, для случая i Ф j :
Gy = 2ц * вij = 2|G(t - T)deij (t), dGt) = 1 - R(t), R(t) = К(t) -(t - т)К(т)^т +...,
0 dt 0
J (t) = 1 + К (t),
где G(t) - функция памяти материала, R(t), К (t) - ядра релаксации и ползучести материала. Кроме того, пусть отношение значения модуля на бесконечности к значению модуля в начальный момент (для регулярных моделей) определяется параметром w = G(^) / G(0) .
Введем вектор напряжений tn (х, t) в точке х на элементарной площадке с единичной нормалью n(х):
tn (х, t) = Tn(х)u(х, t) = n(x)X(t) * div u(х, t) + 2^(t) * д u(X, t) + ^(t) * [n(х) x rot u(х, t)]. (2)
^ ' д n(х)
Если х e Гк , то под n(х) будем понимать единичный вектор внешней (по отношению к О к ) нормали к границе Гк. Вектор напряжений tn (х, t) = (t1,12, t3) на границе Гк, соответствующий этой нормали, обозначим через tk (х, t).
Будем рассматривать следующие типы граничных условий для О к : uk (х, t) = ft (х, t), х еГu пГк ; tf (х, t) = gf (х, t), х еГ^пГк ; uf (х, t) = uS (х, t), tf (х, t) = -tf (х, t), х е Гк5 .
Здесь Гu и Гст - части границы Г тела О, по которым заданы соответственно перемещения и поверхностные силы; Г^ - граница жесткого контакта частей О к и О s.
Функции fi (х, t) и gl (х, t) являются заданными функциями координат и времени. Математическая модель пороупругой среды
Система дифференциальных уравнений в преобразованиях Лапласа (параметр s) для смещения мг- и порового давления pP имеет следующий вид [3]:
К,+(К+3 Дм- (*-ßßb -s hp-ßpfu =-Fi,
в ~ Ф13 ~ , в ~ _ ~ п_ кР/ Ф2 5 2
-Рм--—Р -(а- Р)™1л _ -а , в _
К ' ' Ф2 3 + 32 к (Ра +фр1)'
где л, К - константы упругости, ф - пористость, к - проницаемость, а - эффективный коэффициент напряжений, р, ра, р^ - плотности пористого скелета, присоединенной массы
и жидкой среды, , ~ - плотности источников.
Фундаментальное решение для этой системы построено в [3], однако при построении матриц сингулярных решений допущены ошибки, что, в свою очередь, привело к ошибочному интегральному представлению и ГИУ для порового давления. 2. Интегральная формулировка
Применим к исходным уравнениям интегральное преобразование Лапласа:
<х>
/ (3) _ | / (I)е ,
о
где 3 - параметр преобразования Лапласа.
В качестве метода решения будем использовать метод граничных интегральных уравнений [5], в основе которого лежит сведение краевой задачи для дифференциального уравнения движения к интегральному уравнению относительно граничных функций.
Вектор перемещений во внутренних точках области связан с граничными значениями перемещений и усилий следующим образом.
7 С ~к С 7
й (х, з ) _ ] и у (X, уз )у (у, з )йу£, - ] Ту (х, у, з )й/. (у, з )ёу8 , I _ 1,2,3 , х еО к . (3)
Гк Гк
Здесь иу и - соответственно компоненты тензоров фундаментальных и сингулярных решений уравнения (1). Тензор Т выражается из (2) через тензор и с помощью оператора напряжений Тп
Т(хy, з) _ [Тп(х)и(y,3 )]' , где верхний значок «'» означает транспонирование. Формула (3) дает следующее ГИУ:
Су (х)йк (х, 3 ) +|Ту (х, у, 3 Уйк (у, 3 _ |иу (х - у, 3 (у, 3 , I _ 1,2,3 , х е Гк . (4)
Гк Гк
Интеграл в левой части (4) является сингулярным, т.е. понимается в смысле главного значения по Коши, а Су (х) - известный коэффициент при внеинтегральном члене. Если в
точке х поверхность имеет единственную касательную плоскость, то Су (х) _ 5 у / 2 . ГИУ (4) позволяет разработать эффективную численную методику для определения неизвестных
амплитуд граничных перемещений и поверхностных сил. Решением исходной начально-краевой задачи будет вектор-функция и (х, I), полученная путем применения к решению (3), (4) обратного преобразования Лапласа:
а+гю
/(()=2- I/(*)л*
(5)
Для численного обращения (5) будем использовать алгоритм, предложенный Дурбиным [7]. Граничные интегральные уравнения пороупругости
Интегральное представление прямого подхода имеет следующий вид:
с1Г,
Компоненты матриц - ядер интегрального представления - можно найти в [3, 8].
Ядра интегральных представлений допускают следующее выделение особенностей:
^ 1 0 1 + у
+ 0(г0), и*. = —
4пв г
= I ~и * - Р' с1Г - I ~ * - р' иг
_ р _ Г и/ - Р ? _ _ р _ I Г Т/ - (г _ _р _
Р* = 0(г0), и{ = 0(г0), Р1 = Р/*- + 0(г0), и* = 0_1 +ЛУ г\,г} + 5*(3-4v)}- + 0(г0),
с? * =
8пЕ (1 -V)
1 + V {а(1 - 21/)(г„г; -п*) - 2в(1 - v)(гnг] + п*)} + 0(г0),
8пЕ (1 -V)
2 г
V = р/* 8пв
1 - 2v
1 'г,*
1 -V
Т/ г г г, * + пг ~ ' ^ + 0(г 0),
а + в(1 - 2v) 11 г
1 -V
- ((1 - 2v)Sij + 3г,г г* Уп + (1 - 2v)(r п - г,г п* )
Т * =
4 8п(1 - v)г 2
Итоговая система ГИУ примет вид:
с* (У) 0
+ 0(г 0), (( ' + 0(г0).
4пг
0
с( У)
иг^, х) p(t, х)
+
II
0 Г
Т* ($-Т,У, х) ё* ^-Т,У, х)
II
0Г
и* ^ -т, У, х) - Р* ^ -т, У, х) и { ^ -т, У, х) - Р? ^ -т, У, х)
Т/ ^ -т, У, х) ^ -т, У, х) йГГт.
и (т, х ) P(т, х)
йГГт =
ti(т, х) д(т, х)
На основе этой системы строится дискретный аналог, детальное описание которого можно найти в [5].
3. Численные эксперименты
Задача о торцевой ударной силе по призматическому телу с жестко закрепленным концом
Полагалось р = 1Н/м2 (рис. 1), параметры материала: р = 7850 кг/м3, v = 0, Е = 2,11 -1011 Па. Задача имеет аналитическое решение, и известно ее МГЭ-решение в
а-гю
сочетании с методом квадратур сверток [3]. Использована равномерная сетка с 224 ГЭ. Форма нагрузки: p(t) = H(t), t > 0 .
На рисунках кривыми с маркерами представлены кривые аналитических решений для соответствующих параметров вязкости. В [3] приведены результаты для упругого случая с использованием трех неравномерных гранично-элементных сеток из 324, 112 и 42 треугольных граничных элементов. Расчеты показали, что для достижения не меньшей точности, чем в [3], достаточно взять равномерную сетку с 224 граничными элементами.
Поведение перемещений для модели Кельвина-Фойгта продемонстрировано на рис. 2 (для кривой 1 в = 500 ; для кривой 2 в = 100 ; для кривой 3 в = 3; для кривой 4 в = 0,5 ). С уменьшением характерного времени ползучести тs или ростом в = 1/Т материал все более отчетливо начинает вести себя как упругий материал на длительных модулях. При в ^ <х>(тЕ ^ 0) получаем чисто упругий случай.
Рис. 1.
lime
Рис. 2.
Поведение перемещений для модели стандартного вязкоупругого тела продемонстрировано на рис. 3 (для кривой 1 у=0,01; для кривой 2 у=0,3; для кривой 3 у =10; для кривой 4 у = 100). С уменьшением характерных времен релаксации та или ростом
у = 1/ та материал все более отчетливо начинает вести себя как упругий на длительных
модулях. На рис. 3 показано, как перестраивается картина отклика в перемещениях. Отклик перемещений на длительных модулях имеет и большую амплитуду, и больший период колебаний. Перемещения для степенной модели показаны на рис. 4 (при к = 17 для кривой 1 а = 0,95 ; для кривой 2 а = 0,7; для кривой 3 а = 0,3 ). С приближением наследственного ядра к ядру Больцмана отклик и качественно и количественно описывает упругий отклик.
Соответствующий отклик напряжений в упругом случае, полученный для равномерной ГЭ-сетки с 224 элементами, представлен на рис. 5. Аналитическое решение
(кривая 1) сравнивается с результатами ряда исследований, полученных в [3] на основе метода квадратур сверток и МГЭ, с адаптированной сеткой из 324 треугольных ГЭ (кривая 3) и с результатами, построенными здесь, на равномерной сетке из 224 четырехугольных ГЭ (кривая 2). Влияние вязкости на отклик напряжений, снятых с жестко закрепленного конца, представлено на рис. 6-8 соответственно для моделей Кельвина-Фойгта (для кривой 1 в = 0,05; для кривой 2 в = 3; для кривой 3 в = 100 ; для кривой 4 в = 500 ), стандартного вязкоупругого тела (для кривой 1 у = 0,01; для кривой 2 у = 0,3; для кривой 3 у = 10; для кривой 4 у = 100) и степенной модели (при к = 17 для кривой 1 а = 0,3; для кривой 2 а = 0,7; для кривой 3 а = 0,95 ). Тенденция, выявленная на картинах полей перемещений, прослеживается и в картинах полей напряжений.
Рис. 3.
Рис. 4.
Рис. 5.
Рис. 6.
Рис. 7. Рис. 8.
Решение пороупругой задачи
Рассмотрена задача, изображенная на рис. 13, со следующими параметрами материала: К = 8-109И/ж2; О = 6-109И/ж2; Я = 4,7-108И / ж2; к = 1,9-10-10ж4/Нс; р = 2458кг/ж3; р{ = 1000кг/ж3; ф = 0,19; а = 0,867. Нагрузка, действующая на тело,
^2 = 1И / ж . Граничные условия в области Лапласа имеют вид:
Р2(У2 = 0) = 0 ~2(У2 = 0) = 0 *~2(У2 = 1) = -1/^ Р(У2 = 1) = 0 .
Гранично-элементная сетка, состоит из 504 элементов. Результаты расчетов перемещений и давлений приведены на рис. 9, 10 соответственно. Графики, изображенные сплошной линией, соответствуют точному решению, штриховой - решениям, полученным по новым ядрам, пунктирной - по ядрам из [3].
Рис. 9. Рис. 10.
Задача о действии скачка давления на торец составного призжатического тела Рассмотрим задачу о торцевом ударе силой Р^) = 1Н / ж2 составного призматического тела с жестко закрепленным концом (рис. 11). Эта задача имеет аналитические решения как для упругой так и для вязкоупругих моделей. Для сравнения
численных результатов с аналитическими рассмотрим подобласти с одинаковыми параметрами материала: Е =2,11 • 1011 Н / м2; v = 0,0; р = 7850кг / м3. Задача решается в безразмерных величинах:
х = х / а, и = Еи / (р0 а), р = р / р 0, 1 = с 21 / а, а = а/ р 0, где а=10, р 0=1 Н/м2.
Безразмерные параметры получились следующие: Е = 1;V = 0;р = 0,5 . Гранично-элементная сетка представлена на рис. 12. Каждая из подобластей содержит по 72 элемента и 88 точек на четверти сетки, таким образом, вся геометрическая модель содержит 576 элементов.
Закрепленный
Рис. 11. Рис. 12. Рис. 13.
На рис. 14-15 приведено сравнение аналитического упругого решения (штриховая линия) с численными ГЭ-результатами, полученными на основе метода квадратур сверток (кривая 1) и на основе метода Дурбина (кривая 2). На рис. 14 приведены графики напряжений (в области контакта подобластей), а на рис. 15 графики перемещений (в области контакта подобластей).
1. С »10° С кМ"'
Рис. 14. Рис. 15.
На рис. 16 приведены графики напряжений для модели Кельвина-Фойгта при различных значениях параметра вязкости в в областях заделки, на рис. 17 приведены графики перемещений на свободном конце консоли. На рис. 18 приведены графики напряжений для модели стандартного вязкоупругого тела при различных значениях
параметра вязкости у в областях заделки, на рис. 19 приведены графики перемещений на свободном конце консоли. На рис. 20 приведены графики напряжений в области заделки, а на рис. 21 графики перемещений на свободном конце консоли для степенной модели при фиксированном значении параметра к=17 и при различных значениях параметра а. Аналогичным будет влияние параметра к при фиксированном значении параметра а. По итогам экспериментов, результаты которых приведены на рис. 18-19, можем заключить, что численно продемонстрирован эффект перестройки волновых полей граничных и внутренних перемещений и напряжений, когда свойства вязкоупругого материала изменялись с мгновенных модулей на длительные. В откликах перемещений изменялись (увеличивались) амплитуда и период искомой функции.
Рис. 16. 1 - упругое ГЭ-решение, ГЭ-решение для модели Кельвина-Фойгта 2 - в = 100, 3 - в = 10.
Рис. 17. 1 - упругое ГЭ-решение, ГЭ-решение для модели Кельвина-Фойгта 2 - в = 100, 3 - в = 10.
Рис. 18. 1 - упругое ГЭ-решение, ГЭ-решение для модели стандартного вязкоупугого тела 2 - у = 200, 3 - у = 50. 4 - у = 1.
Рис. 19. 1 - упругое ГЭ-решение, ГЭ-решение для модели стандартного вязкоупугого тела 2 - у = 200, 3 - у = 50. 4 - у = 1.
Рис. 20. 1( ) - упругое ГЭ-решение, ГЭ- Рис. 21. 1( ) - упругое ГЭ-решение, ГЭ-решение для степенной модели при решение для степенной модели при
постоянном к = 17 2(......) - а = 0,3 ; постоянном к = 17 2(......) - а = 0,3 ;
3(-.-.-) - а = 0,7 ; 4(-—) - а = 0,95 . 3(-—-) - а = 0,7 ; 4(-—) - а = 0,95 .
Задача о действии вертикальной силы на поверхность упругого полупространства
Рассмотрим задачу о действии вертикальной силы Р(^) = Р0 /(^), Р0 = 1000И / ж2 на поверхность полупространства (рис. 22).
Параметры материала: Е = 2,5 -108 Н / ж2; V = 0,298; р = 1884кг / ж3 ( с1 = 425ж / с; с2 = 228ж /с, сЯ = 211ж / с ). Исследуется точка А полупространства на расстоянии 15ж от источника силы.
Рассматриваются два варианта нагружения /^) = И(^) и /^) = Н(^) - И(^ - 0,05). Гранично-элементная сетка строится с учетом двух плоскостей симметрии и содержит 864 элемента и 913 точек (на четверти сетки). Геометрическая модель исследуемой задачи показана на рис. 23.
Рис. 22. Рис. 23.
Для случая нагрузки /(^) = Н(^) на рис. 24 дано сравнение численных результатов, полученных методами квадратур сверток и Дурбина с численным решением из [3]. При
решении в [3] использовалась изопараметрическая ГЭ-схема с раскрытием сингулярности на основе подхода из [9].
Для случая нагрузки /(V) = Н(V) - Н(V - 0,05) на рис. 25 дано сравнение численных результатов, полученных методами квадратур сверток и Дурбина.
Проведенное исследование показало, что созданные ГЭ-схемы на основе метода Дурбина и метода квадратур сверток дают близкие результаты. Сравнение с ГЭ-результатом из научной литературы позволяет остановиться на выборе ГЭ-схемы для подобных задач и утверждать, что ГЭ-решение из [3] дает заниженный результат по амплитуде волны Рэлея, а также завышенный нефизичный результат по падению амплитуды сразу за фронтом волны Рэлея.
Рис. 24. Сравнение ГЭ-решения
—•—•— Schanz [3],-ГЭ-решение по
методам Дурбина и квадратур сверток (полное совпадение). Задача о действии вертикальной
Рис. 25.
1 - метод Дурбина, 2 - метод квадратур сверток.
силы на
поверхность вязкоупругого полупространства
Рассмотрим решение задачи при действии нагрузки Р(^) = Р0(Н(V)-Н(V-0,0085)), где Р0 = 1 Н/м2 и параметрами материала: Е = 1,38• 108Н/м2; 1/ = 0,35; р = 1966кг/м3 (с1=335,64м/с ; с2 = 161,24ж/с, ск = 150,5м/с ). В качестве точки наблюдения возьмем точку (2,3333; 2,3333; 0).
На рис. 26 приведены численные результаты для модели Кельвина-Фойгта (Я(го) = Я,/и(ю) = ¡и) при разных значениях параметра вязкости в(кривая 1 соответствует упругому случаю, кривая 2 - в = 100 ; 3 - в = 1; 4 - в = 0,1; 5 - в = 0,01).
На рис. 27 приведены численные результаты для модели стандартного вязкоупругого тела при разных (¡и(<х>) = ¡и,Я(ю) = Я, Я(ю)/ Я(0) = ¡(^)/ ¡(0) = 0,0625) значениях параметра
вязкости у (кривая 1 соответствует упругому случаю, кривая 2 - у = 100; 3 - у = 1; 4 -у = 0,1; 5 -у = 0,01).
Рис. 26. Рис. 27.
На рис. 28, 29 приведены численные результаты для степенной модели при разных значениях параметра вязкости к и а (на рис. 28 кривая 1 соответствует упругому случаю, кривая 2 - к = 5, а = 0,95; 3 - к = 10, а = 0,95; 4 - к = 15, а = 0,95; на рис. 29 кривая 1 соответствует упругому случаю, кривая 2 - к = 5, а = 0,5; 3 - к = 5, а = 0,7 ; 4 - к = 5 , а = 0,95).
Рис. 28. Рис. 29.
Численные ГЭ-эксперименты показали, что если для хевисайдовской нагрузки максимум в отклике перемещений приходится на волну Рэлея, соответствующую скачку подъема силы в начальный момент времени, то для нагрузки вида разности хевисайдовских сил максимум в отклике перемещений приходится на волну Рэлея, соответствующую скачку падения силы в момент I = а . Разница величин амплитуд соответствует величине отклика на постоянную (статическую) часть отклика. Проведены исследования влияния длины импульса нагрузки и вязкости материала на характер поведения отклика перемещений. Использование модели Кельвина-Фойгта может существенно изменить отклик поверхностных волн, однако
имеются диапазоны характерных времен ползучести и релаксации, когда изменения в отклике проявляется лишь в уменьшении амплитуды. Использование модели стандартного вязкоупругого тела и слабосингулярного ядра приводит к изменению скорости отклика перемещений и амплитуды.
Задача о действии вертикальной силы на поверхность полупространства с полостью
Рассматриваются два варианта полости - сферическая и кубическая. В качестве действующей силы рассматривается вертикальная сила Р(/) = Р0 /(/) на площади £ = 1м2 дневной поверхности полупространства. Внутри полупространства на глубине Н = 7,5 м (центр полости) расположена сферическая полость радиуса Я=5м (рис. 30). Исследуются перемещения на поверхности полупространства на расстоянии 15м от границы действия силы.
Гранично-элементная сетка строится с учетом двух плоскостей симметрии. Четверть сетки содержит для полупространства - 864 элемента и 913 точек, для полости - 150 элементов и 171 точку рис. 30.
Рассматривается задача, когда на глубине Н = 7,5 м (центр куба) расположена кубическая полость с длиной ребра 10 м (рис. 31).
Рис. 30. Рис. 31.
Параметры материала выбраны следующими: Е = 2,5 • 108 Н / м2; у = 0,298; р = 1884кг / м2 (^ = 425м / с; с2 = 228м / с). Рассматриваются два варианта нагружения Р(I) = Р0/(I), /(I) = Н(I) и /(I) = Н(I) - Н(I - 0,05)), при Р0 = 1Н / м2.
На рис. 32, 33 приведены графики перемещений во времени, полученные на основе метода Дурбина при нагружении /(/) = Н(/) и /(/) = Н(/) - Н(/ - 0,05) соответственно, причем кривая 1 - решение задачи для полупространства без полости, кривая 2 - решение задачи со сферической полостью и кривая 3 - решение задачи с кубической полостью.
На рис. 34 приведено исследование перемещений для полупространства со сферической полостью, при удалении от источника нагрузки при нагружении /(/) = Н(/). Кривая 1 соответствует удалению 15м, кривая 2 - 19м, кривая 3 - 21,25м, кривая 4 - 23,5м.
Рис. 32. Рис. 33.
На рис. 35 приведено исследование перемещений для полупространства с кубической полостью при удалении от источника нагрузки при нагружении /(^) = И (^). Кривая 1 соответствует удалению 15ж, кривая 2 - 19ж, кривая 3 - 21,25ж, кривая 4 - 23,5ж, кривая 5 -25,5ж, кривая 6 - 27,5ж.
Из графиков видно, что для сферической и кубической полостей форма отклика в сравнении с задачей для полупространства поменялась. Форма полости влияет на форму отклика: для задачи со сферической полостью проявилось два горба в момент прихода поверхностной волны, а для кубической полости проявился временной промежуток с постоянными перемещения.
Рис. 34. Рис. 35.
Задача о штажпе на полупространстве
Рассмотрим задачу о действии вертикальной силы на деформируемый штамп (параметры материала: Е = 3-108 И / ж2, v = 0,2, р = 2000кг / ж3) в форме параллелепипеда (2 х2х1 ж), расположенный на деформируемом полупространстве (параметры материала:
Е = 1,38-108Н/ж2, v = 0,35, р = 1966 кг /ж ). Расчеты проводились для нескольких ГЭ-
сеток, сведения о которых приведены в таблице. Сетка №1 приведена на рис. 36, сетка №4 -на рис. 37.
Таблица
Номер сетки Кол-во элементов
Штамп Полупространство
1 54 90
2 216 360
3 96 864
4 384 3456
Рис. 36. Рис. 37.
Параметры свойств материалов приведены в таблице 1. Рассмотрим случай нагружения РУ) = Р0 (Н(V) - Н(V - 0,0085)); Р0 = 1Н / м2. Задача решается в безразмерных величинах.
В качестве координат исследуемой точки взяты (2,33; 2,33; 0). За начало координат выбран центр контактной грани штампа.
На рис. 38 приведены результаты расчетов для упругого случая, выполненные для сетки 1, причем кривая 1 соответствует расчетам по методу квадратур сверток, кривая 2 - по методу Дурбина, кривая 3 - ГЭ-решение из [10]. Исследовалась сходимость на выбранных сетках. Установлено, что известное из литературы ГЭ-решение дает завышение амплитуды отклика поверхностной волны.
Влияние вязкоупругих свойств на отклик продемонстрировано на рис. 39-41.
Рис. 38.
Рис. 39.
1 - упругое ГЭ-решение, модель Кельвина-Фойгта 2 - в = 100, 3 - в = 1, 4 - в = 0,01
Рис. 40.
1 - упругое ГЭ-решение, модель стандартного вязкоупругого тела 2 - у = 100, 3 - у = 1, 4 - у = 0,01.
Рис. 41.
1 - упругое ГЭ-решение, степенная модель 2 - к = 1;а = 0,5; 3 - к = 1, а = 0,75; 4 - к = 1;а = 0,95 ; 5 - к = 5;а = 0,95 ; 6 - к = 10;а = 0,95.
Использование модели Кельвина-Фойгта может существенно изменить отклик
поверхностных волн, однако имеются диапазоны характерных времен ползучести, когда изменения в отклике проявляются лишь в уменьшении амплитуды. Использование модели стандартного вязкоупругого тела и слабосингулярного ядра приводит к изменению скорости отклика перемещений и амплитуды.
Задача о реакции защитного корпуса атожной станции теплоснабжения на действие ударной силы
Геометрическая модель корпуса атомной станции теплоснабжения и форма нагрузки представлены на рис. 42. Высота корпуса 34,5ж, размеры в плане 15,6х15,6ж, = 30,5ж, Ъ2 = 9,3ж, И3 = 4,215 ж, ЪА = 4,01 ж, Ъъ = 8,445 ж, Ъ6 = 2ж, Ь7 = 2,2ж, 1Х = 9,6ж,
12 = 11,6 м, г = 3,6 м, г2 = 3,8 м, г3 = 4,29 м, г4 = 5,8 м . Параметры материала - р = 2 • 103 кг/м3,
V = 0,2, Е=3 • 1010 Н/м2 (с = 4082,5м / с, с2 = 2041,2м / с).
Задача решается в безразмерных величинах: х = х / а, м = Ей /(р0 а), р = р / р0,
V = с2V/а , а = а/р0, где а = 25, р0 = 4 • 106Н/м2.
Половина ГЭ-сетки АСТ содержит 226 элементов и 257 точек (рис. 43 ), половина ГЭ-сетки полупространства содержит 940 элементов и 982 точки (рис. 44).
Исследуются вертикальные перемещения в точке В. Численные результаты решаемой задачи (кривая 1, рис. 45) сравниваются с результатами, полученными на модели, у которой полупространство моделируется амортизаторами (кривая 2, рис. 43) [5, 11].
0 02 0.04 0.06 о.ое 0.1 0 12 0.14
Г
Рис. 43. Рис. 44. Рис. 45. 1 - с учетом деформаций
полупространства, 2 - без учета
деформаций полупространства.
Расчеты продемонстрировали качественное совпадение и количественные отличия в
искомых граничных полях.
Библиографический список
1. Biot M. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range // J. Acoust. Soc. Am. - 1956. V. 28, № 2. - P. 168-178.
2. Biot M. Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher-frequency range // J. Acoust. Soc. Am. - 1956. V. 28, № 2. - P. 179-191.
3. Schanz M. Wave Propogation in Viscoelastic and Poroelastic Continua / Berlin Springer, 2001. - 170 p.
4. Manolis G.D., Beskos D.E. Integral Formulation and Fundamental Solutions of Dynamic Poroelasticity and Thermoelasticity // Ada Mechanica. 1989. № 76. P. 89-104.
5. Баженов В.Г., Игумнов Л.А. Методы граничных интегральных уравнений и граничных элементов в решении задач трехмерной динамической теории упругости с сопряженными полями. М.: Физматлит, 2008. - 352 с.
6. Локшин А.А., Суворова Ю.В. Математическая теория распространения волн в средах с памятью // Москва. - Изд-во МГУ. - 1982. - С. 152.
7. Durbin F. Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate's method // The Computer Journal. - 1974. - Vol.17, 4. - P. 371-376.
8. Аменицкий А.В., Васильев А.Н., Игумнов Л.А., Карелин И.С. Граничные интегральные уравнения для решения динамических задач трехмерной теории пороупругости // Н.Новгород: Изд-во ННГУ. Проблемы прочности и пластичности. Межвуз. Сборник, 2009. - Вып. 71.,. - С. 164-171.
9. Birgisson B., Siebrits E., Peirce A.P. Elastodynamic direct boundary element methods with Enhanced numerical stability properties // International journal for numerical methods in engineering. - 1999. - Vol. 46. - P. 871-888.
10. Gaul L. [et al.] Boundary Element Methods for the Dynamic Analysis of Elastic, Viscoelastic, and Piezoelectric Solids // Encyclopedia of Computational Mechanics: Edited by E. Stein, R. de Borst and Thomas J. R. Hughes. Solids and Structures. - Jhon Wiley & Sons, Ltd., 2004. -Vol. 2. - P. 751-769.
11. Хуторянский Н.М. Метод гранично-временных элементов в пространственных задачах нестационарной динамики упругих и вязкоупругих тел: автореф. дисс...доктора техн. наук: 01.02.04 / Хуторянский Наум Маркович. - Рига, 1988. - 32с.
Рекомендовано к публикации программным комитетом XVI Международного симпозиума «Динамические и технологические проблемы механики конструкций и сплошных сред» им. А.Г. Горшкова.
Сведения об авторах
Александр Александрович Белов Александр, научный сотрудник НИИ механики Нижегородского государственного университета им. Н.И. Лобачевского, к.ф-м.н., тел. (831)465-76-55, e-mail: belov_a2@mech.unn.ru
Леонид Александрович Игумнов, директор НИИ механики Нижегородского государственного университета им. Н.И. Лобачевского, д.ф-.м.н., тел. (831)465-76-55, e-mail: igumnov @mech.unn.ru
Иван Сергеевич Карелин, младший научный сотрудник НИИ механики Нижегородского государственного университета им. Н.И. Лобачевского, тел. (831)465-76-55, e-mail: igumnov@mech.unn.ru
Светлана Юрьевна Литвинчук, ученый секретарь НИИ механики Нижегородского государственного университета им. Н.И. Лобачевского,к.ф.м.н., тел. (831)465-76-55, e-mail: litvinchuk @mech.unn.ru
