Конечноэлементное моделирование процесса распространения упругих волн в одномерной среде со сферической симметрией

Тип работы:
Реферат
Предмет:
Физико-математические науки


Узнать стоимость

Детальная информация о работе

Выдержка из работы

КОНЕЧНОЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЯ ПРОЦЕССА РАСПРОСТРАНЕНИЯ УПРУГИХ ВОЛН В ОДНОМЕРНОЙ СРЕДЕ СО
СФЕРИЧЕСКОЙ СИММЕТРИЕЙ
Гуш Максим Николаевич
аспирант, кафедра прикладной математики Новосибирского государственного технического университета,
РФ, г. Новосибирск 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)
ик = °& gt-л
dn
= 0- (3)
Получим вариационную формулировку в форме уравнения Галёркина [2, с. 84]. Для этого потребуем, чтобы невязка дифференциальных
!! сгеа1ес1 Ьуее уегсюп о{
д РооРгеегег
уравнений была ортогональна (в смысле скалярного произведения пространства Н0) некоторому пространству Ф пробных функций V, т. е.
д2и и
Р^ъ — grad и) + =) = 0,
(4)
Расписывая скалярное произведение, получим:
Г д2и
I У & quot- ^
дг2
— I grad и) — V — + | 2 т —
-у-йП = 0
п
п
п
По свойствам дивергенции:
(5)
— V = А gradv —
(6)
где: V — скалярное поле,
А — векторное поле. Тогда справедливо:
— § ймА — = § А — - §
По теореме о дивергенции:
№п (Ии (Ау)(1П = § яАу-п & lt-Я9
где п — вектор нормали к поверхности 5. Учитывая (8), соотношение (7) примет вид:
— § & amp-уа — = § А — 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 и-%г& amp-&- у0аа
а а
[ & amp-и ^ Г и (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 + и* (г& gt-/?((14)
к-2
к-1,
где функции иК 3(г), иК 2(г), иК 1(г), иК (г) являются значениями искомой функции при? = tk-3,t = гк-2,г = гк-1,г =, а функции
лз (?), л2. (?), л1(?), ло (?) являются кубическими полиномами и имеют вид:
(* - гк-2)(г — гк-1)(1 — гк)
(?к-3 — ^-2)^-3 — {к-О^к-З — ?кУ
Лз =
Ю =
(?к-2 — ^к-з)(^к-2 — ^к-1)(^к-2 — ?кУ, а _ & amp- - гк-з)& amp- - гк-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& amp- +
(21)
Перейдём к конечноэлементной СЛАУ. При построении численного решения по методу Галёркина минимум невязки ищется не на пространствах Нд и Н1, а на аппроксимирующих их конечномерных подпространствах Удн и
2
У?. При этом конечномерное пространство У*, подпространствами которого являются Уд* и У О, мы определим как линейное пространство, натянутое на базисные функции -ф^Ь = 1.п.
Заменим функции и Е Н} аппроксимирующие их функции и* ЕУ1, а
функцию у0 Е Н° функцией V* Е У*
А9
И (~гП.
, д2Чз
,-2 +д2& amp-иП*-1 + д2? иП*
I
а
Г Г и*,"
+ I ц grad ип, К — %га& amp-уО — йа + 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М'-& quot-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. 13 518^ + 2. 1 & quot-Н
|цаналит _ ц^/2 II = 1,4 263* + 2. ¦ & gt-'-?2
Iцаналит _ цЛ/4|| = 6Л1425е + 1,
¦ & quot-?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 с.

ПоказатьСвернуть
Заполнить форму текущей работой