CC BY
11
КОНЕЧНОЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЯ ПРОЦЕССА РАСПРОСТРАНЕНИЯ УПРУГИХ ВОЛН В ОДНОМЕРНОЙ СРЕДЕ СО
СФЕРИЧЕСКОЙ СИММЕТРИЕЙ
Гуш Максим Николаевич
аспирант, кафедра прикладной математики Новосибирского государственного технического университета,
РФ, г. Новосибирск E-mail: mgush@mail.ru
FINITE-ELEMENT MODELING OF ELASTIC WAVES IN ONE-DIMENSIONAL MEDIUM WITH SPHERICAL SYMMETRY
Maxim Gush
postgraduate student Department of Applied Mathematics Novosibirsk State Technical University,
Russia, Novosibirsk
АННОТАЦИЯ
В данной статье рассматривается конечноэлементное моделирование процесса распространения упругих волн в одномерной среде, обладающей сферической симметрией. Приводятся вывод эквивалентной вариационной постановки и построение конечноэлементного аналога. Оценивается порядок аппроксимации разработанной схемы.
ABSTRACT
The article discusses the finite element modeling of elastic wave propagation in one-dimensional medium with spherical symmetry. An equivalent variation formulation has been developed and finite element discretization has been carried out. The order of approximation has been estimated.
Ключевые слова: упругие волны; метод конечных элементов; вариационная постановка; дискретный аналог.
Keywords: elastic waves; finite element method; variation formulation; finite element discretization.
Распространение упругих волн в изотропной среде обладающей сферической симметрией описывается одномерным волновым уравнением и
^ created by free version of
S DociFreezer
является наиболее простым случаем распространения упругих волн. При этом в силу симметрии в такой среде будут наблюдаться только продольные упругие волны, в которых смещение направлено вдоль вектора распространения самой волны. В практических задачах среды, обладающие сферической симметрией, встречаются крайне редко, однако, данная задача может быть использована для исследования качественного характера получаемых численных решений в задачах распространения упругих волн и для отладки более сложных процедур (моделирование распространения упругих волн в двух - и трёхмерных средах, решение обратной задачи сейсмической разведки и т. д.).
Математическая модель одномерной задачи в сферических координатах описывается следующим уравнением [1, с. 513]:
д2и и пл
§гай и) + = (1)
где: и(т^) — деформация среды в точке г,
f (г,1) — функция задающая внешние силы, ^ = 2С
1 — 2 V
р — плотность, С — модуль сдвига, V — коэффициент Пуассона. Начальные условия:
ди
Краевые условия:
йи
= 0. (2)
ик = °>л
dn
= 0- (3)
Получим вариационную формулировку в форме уравнения Галёркина [2, с. 84]. Для этого потребуем, чтобы невязка дифференциальных
!!\ сгеа1ес1 Ьу ^ее уегсюп о{
д РооРгеегег
уравнений была ортогональна (в смысле скалярного произведения пространства Н0) некоторому пространству Ф пробных функций V, т. е.
д2и и
Р^ъ — grad и) + = ) = 0,
(4)
Расписывая скалярное произведение, получим:
Г д2и
I У " ^
дг2
— I grad и) - V - + | 2т\ —
-у-йП = 0
п
п
п
По свойствам дивергенции:
(5)
- V = А gradv —
(6)
где: V — скалярное поле,
А — векторное поле. Тогда справедливо:
— § ймА - = § А - — §
По теореме о дивергенции:
№п(Ии(Ау)(1П = §яАу-п <Я9
где п — вектор нормали к поверхности 5. Учитывая (8), соотношение (7) примет вид:
— § &\уа - = § А - gradv — йБ.
Тогда
§ div(;qgтad и) - V - = § ц grad и - grad V - — §V - п йБ. (10)
п
п
^ йп
Подставим полученные соотношения (10) в (5):
(7)
(8)
(9)
- аП + № grad и - grad р-йП — М + (11)
йи йп
аа = о.
Поскольку на границе краевыми условиями не определяется значение Ау - п, то слагаемое ^ Ау - п йБ следует исключить из уравнения, потребовав, чтобы пространство пробных функций Ф содержало только функции, которые принимали бы только нулевые значения на границе . Обозначим их у0 .
Обратим внимание, что в полученное выше уравнение входят производные пробных функций у. Поэтому в качестве пространства пробных функций Ф мы можем выбрать НО — пространство функций, имеющих суммируемые с квадратом производные и равные 0 на границе .
Таким образом, получаем систему вариационных уравнение вида:
д2и
Г о2и Г
I Р^ргр0 I ц grad и-%г&& у0аа
а а
[ &и ^ Г и (12)
— I ц—;-У0-П аБ + I 2ц — - у0аа = 0 ] ап J г2
5 а
Ууо е НО.
В силу (3) приходим к уравнению
Г д2и Г Г и
I р^тУ0 - I ц grad и - grad у0 - йа + I 2ц — - у0 - йа = 0,
) дг2 ) J г2 (13)
а а а
ЧУ О е НО.
Сделаем дискретизацию по времени [2, с. 364]. Будем полагать, что ось времени £ разбита на так называемые временные слои значениями ^, к = 1..К, а значения искомой функции и(г, £) на к-м временном слое обозначим через ик(г,гк) = ик(г)9 которые не зависят от времени, но остаются функциями пространственных координат.
Рассмотрим процедуру построения неявной четырехслойной схемы для решения дифференциального уравнения гиперболического типа.
Представим искомое решение ик(г) на интервале (?к-3^к) в следующем виде:
иг(г, 0 = ик-3(г)г}*(1) + иК-*(г)г}$(0 + ик~г{г)7}^(0 + и* (г>/£((14)
к-2
к-1,
где функции иК 3(г),иК 2(г),иК 1(г),иК(г) являются значениями искомой функции при £ = tk-3,t = гк-2,г = гк-1,г = , а функции
лз(£),л2.(£),л1(£),ло(£) являются кубическими полиномами и имеют вид:
(* — гк-2)(г — гк-1)(1 — гк)
(¿к-3 — ^-2)^-3 — {к-О^к-З — ¿кУ
Лз =
Ю =
(£к-2 — ^к-з)(^к-2 — ^к-1)(^к-2 — ¿кУ а _ & — гк-з)& — гк-2)(г — гк)
т =
(*к-1 — ^к-з)(^к-1 — ^к-2)(^к-1 — ¿кУ (С — Ь-з)^ — гк-2)(г — гк-1)
т =
— ^к-3)(^к — ^к-2)(^к — Ък-1)
(15)
(16)
(17)
(18)
Применим представление (14) для аппроксимации второй производной по времени в уравнениях на временном слое £ = :
д
д2
дг2 иг(г, 1)=-^ (икг-3(гШ1) + икг-2(гШ1) + ик-1(гШ1)
+ ик(г)лШ
д2иг
~дг2
д2т к-з+д_А к-2+дж к-1 + д_А *
дг2 щ + дг2 щ + дг2 щ + дг2 Пг'
Подставим (20) в уравнения (13) и получим:
(19)
(20)
Я
п
I
(дгЩ дг2
■и
Н,к-3
+
д2Ч2 П,к-2 , д_2]1иКк-1 ,
дг2
■и'
дг2
д2Ло к
дг2
+ I иК - gradv0 - +
Г ик п
у0-аа = о,
Ыо е н1
уп - с1& +
(21)
Перейдём к конечноэлементной СЛАУ. При построении численного решения по методу Галёркина минимум невязки ищется не на пространствах Нд и Н1, а на аппроксимирующих их конечномерных подпространствах Удн и
2
У£. При этом конечномерное пространство У*, подпространствами которого являются Уд* и У О, мы определим как линейное пространство, натянутое на базисные функции -ф^Ь = 1..п.
Заменим функции и Е Н} аппроксимирующие их функции и* ЕУ^1, а
функцию у0 Е Н° функцией V* Е У*
А9
И (~\гП.
,д2Чз
^,-2 +д2&иП*-1 + д2£иП*
I
а
Г Г и*,«
+ I ц grad ип,К - %га&уО - йа + I 2ц—— - уО - йа = 0,
УУо Е И}.
• Уп - dа +
(22)
Любая функция у0 Е V* может быть представлена в виде линейной комбинации базисных функций пространства V* :
V* = £ Фи
(23)
1ЕМо
Функции и* Е Уд* так же можно представить в виде линейной комбинации базисных функций пространства Уд*:
п
и
= £
(24)
1=1
причём п — п0 компонент векторов весов, ^ ^.....^ должны быть
фиксированы и могут быть определены из условий
Подставляя (23) и (24) в (22) получим дискретный аналог для (13):
п
ад
ь I I д2п% дф^д-фг 2л \
) = 1 п
дЧ . л-2-Ч. _к-1д2ч1\ (.,. .,. .,„ (26)
7 П
I = 1.. П.
В результате мы получили систему из п уравнений с п неизвестными . Чтобы определить матрицу и вектор правой части полученной конечноэлементной СЛАУ, пронумеруем её уравнения и неизвестные следующим образом.
Систему уравнений (26) можно записать в матричном виде:
Ац = Ъ. (27)
где
Ау = I + ^дтг + Т^К (28)
п
* = + + IМ'"0- (29)
У п
Используемые в МКЭ базисные функции являются финитными, т.е. каждая функция 'фI не равна 0 только на нескольких примыкающих к определяющему её /-му узлу конечных элементах 0к. Следовательно, для фиксированной функции фI только небольшое число функций ф] отличны от 0 в подобласти, где ^¿^О. Тогда из соотношений (28) очевидно, что в каждой строке матрицы Асодержится мало ненулевых элементов, поэтому для хранения и обработки матрицы будет использован разреженный формат хранения.
Поскольку базисные функции кусочно-полиномиальные, то интегралы в (28) практично вычислять, как сумму интегралов по конечным элементам 0к на которые разбита расчётная область и тогда матрица А представляется в виде суммы вкладов Ак от соответствующих конечных элементов 0к:
А =
к
(30)
При этом фактически все ненулевые компоненты матрицы Ак (размера п X и) можно разместить в локальной матрице Акразмера 1к X 1к, где 1к-количество базисных функций, отличных от нуля на конечном элементе П к.
Локальные лагранжевы базисные функции на конечном элементе Пр = [хр,хр+1] представляются в виде [2, с. 68]:
X—+1 X X X—
*1(х)=—£-,Х2(Х) = —-,Их=Хр+1—Хр. (31)
Каждая базисная функция равна единице в узле соответствующем её номеру и нулю в остальных.
Локальные матрицы для линейных лагранжевых базисных функций будем вычислять численно по схеме Гаусс-4.
Протестируем аппроксимацию решения по пространству. Расчётная область: г е [—100,100]; Н = 25; Параметры среды:С = 1е + 9,у = 3.5е — 1, р = 0; Аналитическое решение:
ианалит = г3 + 2г2 + 3 г.
|и™«ит_ий|| =2.13518^ + 2. 1 "Н
|цаналит _ ц^/2 II = 1,04263* + 2. ■ >'¿2
Iцаналит _ цЛ/4|| = 6Л1425е + 1,
■ "¿2
Определим по правилу Рунге порядок аппроксимации:
\\и*/4 — и*\\
^ сгеа!ес1 Ьу йгее уетоп of
ю Роа^геегег
Вывод: Решение с дроблением сетки по пространству сходится к правильному. Порядок пространственной аппроксимации О (к2) совпадает с теоретическим.
Протестируем аппроксимацию решения по времени. Расчётная область: г Е [-100,100]; к = 25; Расчётный интервал времени: £ Е [0,10]; к = 2.5; Параметры среды:С = 1е + 9, и = 3.5е — 1,р = 1000; Аналитическое решение:
Определим по правилу Рунге порядок аппроксимации:
к = 1од2 ,, -гтттт: = 3.39 .
Вывод: Решение с дроблением сетки по времени сходится к правильному. Порядок аппроксимации по времени О (к4) совпадает с теоретическим.
Список литературы:
1. Новацкий В. Теория упругости: учебное пособие. М.: «Мир», 1975. — 872 с.
2. Соловейчик Ю.Г., Рояк М.Э., Персова М.Г. Метод конечных элементов для решения скалярных и векторных задач: учебное пособие. Новосибирск: Изд-во НГТУ, 2007. — 896 с.
