CC BY
10
2017, Т. 159, кн. 2 С. 143-160
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
УДК 519.6:532:533:539.3
РАСЧЕТ ДВУМЕРНЫХ ВОЗМУЩЕНИЙ В УПРУГОМ ТЕЛЕ
А.А. Аганин, Н.А. Хисматуллина
Институт механики и машиностроения КазНЦ РАН, г. Казань, 420111, Россия
Аннотация
В работе исследована возможность повышения эффективности расчета линейных волн в упругом теле за счет применения одной из UNO-модификаций классического метода С.К. Годунова, имеющей второй порядок точности. При построении UNO-схем повышенного порядка точности требование монотонности заменяется условием невозрастания полной вариации. Оценка эффективности предлагаемой UNO-модификации осуществляется путем сравнения результатов ее применения для расчета ряда одномерных и двумерных задач о распространении линейных волн в упругом теле и их взаимодействии с поверхностью тела с точными решениями и результатами расчетов методом С.К. Годунова. Показано, что применение схемы UNO позволяет более чем на порядок снизить затраты компьютерных ресурсов.
Ключевые слова: схема UNO, схема С.К. Годунова, эффективность разностных схем, линейная волна, упругое тело
Введение
Для численного исследования волн в сплошных средах широко применяется классический метод С.К. Годунова [1] (например, в газовой динамике [2, 3], динамике упруго-пластических тел [4-6] и др.). В этом методе внутри каждой ячейки расчетной сетки значения искомых функций задачи (плотности, импульса, энергии и т. д.) предполагаются постоянными (среднеинтегральными по объему). Величины, используемые при разностной аппроксимации пространственных производных на гранях ячеек, определяются в результате решения задачи о распаде разрыва (задачи Римана), обусловленной различием состояния среды в примыкающих к грани ячейках. При выполнении условия Куранта [7] схема С.К. Годунова является устойчивой и монотонной. Вместе с тем она имеет первый порядок аппроксимации, что может привести к сильному размазыванию скачков, контактных разрывов и других особенностей решения в областях с большими градиентами.
В настоящей работе изучается возможность повышения эффективности расчета линейных волн в упругом теле за счет применения применения одной из UNO-модификаций (UNO - Uniformly Non Oscillatory, равномерно безосцилляци-онной) классического метода С.К. Годунова, имеющей второй порядок точности. Идеология построения таких схем изложена в [8]. Реализация одной из них для исследования волновых задач газовой динамики и динамики жидкости представлена в работах [9-11].
При построении UNO-схем повышенного порядка точности требование монотонности заменяется условием убывания (невозрастания) полной вариации (условием TVD - Total Variation Diminishing [12]). При этом, в отличие от TVD-схем, где TVD-условие выполняется строго, в UNO-схемах допускается его нарушение, но лишь
на уровне погрешностей аппроксимации (то есть требуется выполнять асимптотически с измельчением расчетной сетки).
Достигнуть второго порядка точности по времени можно по-разному, например, применяя методы Рунге-Кутты [13]. В настоящей работе с этой целью используются выражения временнЫх производных искомых функций через их пространственные производные [14]. Эти выражения позволяют вычислить значения искомых функций на следующем полуслое по времени в центре ячеек сетки и по обе стороны каждой из их граней. Значения этих функций на самих гранях находятся как решения задачи Римана. Далее рассчитываются численные потоки через границы ячеек. Расчет значений на следующем временном слое осуществляется по явной схеме. При вычислении значений производных по пространственным переменным в ячейке используются ограничители. Применяемый в схеме UNO ограничитель включает в себя аппроксимации производных как первого, так и второго порядков.
Оценка эффективности предлагаемой UNO-модификации осуществляется путем сравнения результатов ее применения для расчета ряда одномерных и двумерных задач о распространении линейных волн в упругом теле и их взаимодействии между собой и с поверхностью тела с результатами расчетов этих задач классическим методом С.К. Годунова.
Изучается возможность повышения эффективности расчета линейных волн в упругом теле по сравнению с классическим методом С.К. Годунова за счет применения одного из вариантов его ИМО-модификации второго порядка точности. Рассматривается случай двумерных волн. В качестве расчетной области тела принимается прямоугольник [х1 < х < хг ] х [уь < у < у г] с границами, параллельными осям декартовой системы координат х, у. Распространение линейных волн в такой области описывается следующими уравнениями:
Здесь £ - время, р - плотность, и, V - компоненты скорости по осям х и у, Бхх, Буу , Бху - компоненты девиатора Б тензора напряжений а, Бхх = ахх + Р, Буу = ауу + Р, Бху = тху, Р - всестороннее (гидростатическое) давление, ахх, ауу, тху - компоненты тензора а, А, ¡л - параметры Ламе, К = А+2/3^ - коэффициент объемного расширения.
Наличие волн в расчетной области обусловлено либо начальным состоянием тела в этой области, либо тем или иным воздействием на границы области. Предполагается, что границами расчетной области могут быть «жесткие стенки», «свободные поверхности» или «искусственные границы».
1. На жесткой стенке задается скорость. Пусть, например, верхняя граница у = уг является жесткой стенкой. Тогда в пределах этой границы
1. Постановка задачи
du daxx dSxy dv dSxy да.
(1)
u(t, x, yt) = U(t, x), v(t, x, yt) = V(t, x),
где U(t,x), V(t,x) - заданные функции.
2. На свободной поверхности задаются нормальные и касательные напряжения. Пусть, например, нижняя граница у = уъ является свободной поверхностью. Тогда в пределах этой границы
ауу (г, х, уъ) = —р(г, х), БХу (г, х, уъ) = Т(г, х),
где р(г,х), т(г,х) - заданные функции.
3. Искусственные границы используются для исключения из рассмотрения удаленных частей тела, влияние которых на расчетную область тела незначительно. Искусственная граница определяется условием отсутствия возмущений, приходящих к этой границе с ее внешней стороны. Такие условия обычно называют неотражающими [15]. Их математическая формулировка будет представлена в следующем разделе.
2.1.
2. Методика расчета Векторная форма уравнений динамики линейных волн и рима-
новы инварианты. Для описания расчетных схем, применяемых в настоящей работе, систему (1) удобно представить в следующем виде:
д а дf де
— +---+ — = 0,
дг дх ду
где
а = (и
> Ях
> Яуу, Яху, Р)
f| О хх Яху 4 2
= — ( —, — , зии, — зии, —ки
р р
I Яху Оуу 24
ё = — —-,—, — о —ку
р р 3 3
(2)
(3)
(4)
Важную роль в предлагаемой методике играет решение одномерной задачи о распаде разрыва, основанное на инвариантах Римана. Кроме того, инварианты Римана используются при аппроксимации краевых условий, а также при получении точных аналитических решений одномерных задач, что полезно для оценки эффективности численного метода.
Если в системе (1) исключить все производные по у, то получится система, описывающая одномерные волны в направлении х. В этом случае инварианты Римана имеют вид
ЗА + 2^
ЗА + 4Ме
11 = -:-Ях
4^
+ Яуу + Р, 12 =
4^
~Яхх + P,
т Охх
/з = и--,
рС1
т , Охх
14 = и +--,
рС1
¡5 = V — ^, 16 = V + Яху.
рС2 рС2
(5)
Для полученной аналогичным образом системы, описывающей одномерные волны в направлении у , инварианты Римана задаются выражениями
ЗА + 4^
ЗА + 2^
11 = Яхх +--;-— Яуу + Р ¡2 = -Яуу + Р 1з = у--~,
4^
4^
Оу
рС1
14 = V + , 15 = и — ^, 1б = и + ^.
рС1 рС2 рС2
(6)
В обоих случаях скорости пространственного перемещения инвариантов Д, ¡2 равны нулю, инварианты 1з, /5 перемещаются со скоростями С1 и С2 соответственно, а инварианты /4, /б - со скоростями — С1 и —С2 соответственно.
т
т
2.2. Аппроксимация уравнений. Расчетная область покрывается равномерной прямоугольной сеткой с шагами Нх в направлении х и ку в направлении у. Значения параметров qгk, II < г < 1Г, Кь < к < К1 с целыми индексами относятся к центрам ячеек, а с полуцелыми (qг±l/2,k, q¿,k±l/2) - соответственно к их вертикальным и горизонтальным граням.
Рассматриваемая ИМО-модификация метода С.К. Годунова для системы (2) состоит из трех шагов.
Шаг 1. Вычисление искомых значений qгk в ячейках расчетной сетки на полуцелом временном слое
^+1/2 = + (двУ ДП
где п - номер временного слоя, ДЬп - шаг по времени.
Для расчета входящей в это соотношение производной по времени используется полученное из (2) выражение этой производной через пространственные производные
% = - (а дЛ + В дЛ
Ы \ дх ду
где А = дf /дq, В = дg/дq, так что
дХ) к+ВГк( ддгУ/гк,
' dq
I
Здесь
ik
A™ = A (qik) , B™ = в (q™) •
Введем обозначения
mm(a, ß) = 2 (sgn(a) + sgn(ß)) min (|a|, \ß|),
A1 = qn qn A1 = qn qn
Ai+1/2,k = qi+1,k — qi,k, Ai-1/2,k = qi,k — qi-1,k,
Ai+1/2,k = mm (qi+1,k - 2qi,k + qi-1,k, qi+2,k - 2qi+1,k + , Ai-1/2,k = mm (qi,k - 2qi-1,k + qi-2,k, qi+1,k - 2qi,k + qi-1,k) ■ Тогда значение производной (dq/dx)^ вычисляется по формуле
^ =rmm Га1+1 /2 k - 1A2+1 /2 k,, Ab „k + 1A2
ОХ ~ V 1+1/2,к - 2 г+1/2,к, ^ г-1/2,к^ 2 ¿-1/2,к
/ гк х V
Аналогично (дq/дx)™k вычисляется и (дq/дy)™k . Шаг 2. Расчет потоков через границы между ячейками
^+1/2 = f (п+1/2 \ п+1/2 = g (п+1/2 \
г+1/2,к 1 ^Чг+1/2,^ , &г,к + 1/2 & ^Чг,к+1/2^ ,
где аргументы функций f и g вычисляются как решение задачи о распаде разрыва,
n+1/2 = n+1/2 _ (f^S^ (7)
qi+1/2ß,k = qi+1,k 2 \dx J k ,
п+1/2 _ п+1/2 + К íдq
%+1/2Ь,к _ %к 1дж
гк
п+1/2 Чг,к+1/2Т
п+1/2
п+1/2 %,к + 1
к^ 2
дq
ду
%,к+1/2в _ %кк ' + л I 1
+ 1/2 , К {"
2 \дУ/гк
(8) (9) (10)
Поток f (3) зависит от четырех аргументов и, V, ахх, Бху. Решение задачи о распаде разрыва для этих величин, используя выражения (5) и опуская все индексы, за исключением Ь и К, запишем в виде
_ 1 ( , (ахх) К (ахх
и _ о \иь + ик +--
2 \ РС1
1 ( . . (Бху )К - (Бху )Ь
V _ о \VL + VК +------
2 V Рс2
&хх _ 2 [(ик - и^)рсс1 + (ахх)я + (axx)L],
Бху _ 2 [^я - VL)pC2 + (Яху)я - (БхуЬ] •
Соответственно, решение задачи о распаде разрыва для четырех аргументов потока й (4), используя выражения (6) и опуская все индексы, за исключением В и Т, получаем в виде
(Бху )Т - (Бху )В
и _ 1 \ит + ив -
V _ 2 I vт + vв +
РС2
(ауу)Т - (ауу)в РС1
ауу _ 2 [(,ит - г°в )РС1 + (ауу)Т + (ауу)в]
Бху _ 2 [(ив - ит)рС2 + (Бху)т + (Бху)в] •
Шаг 3. Расчет искомых значений вектора q™k+1 на новом временном слое по явной конечно-разностной схеме
п+1 п ^+1/2 _ fп+1/2 „п+1/2 _ „п+1/2
qik - qik + г+1 / 2,к 4-1/2,к + &г,к+1/2 £>г,к-1/2
АГ
кх
к
У
2.3. Аппроксимация граничных условий. Для аппроксимации граничных условий вдоль каждой из четырех границ расчетной области вводится примыкающий к ней с внешней стороны слой фиктивных ячеек. Для левой границы х11-1/2,к _ Х1 - это слой ячеек с х-координатой центра хх1-1,к _ щ - кх/2, для правой границы Х1г+1/2,к _ хг - слой ячеек с х-координатой центра Х1г+1кк _ _ хг + кх/2, для нижней границы Угккь-1/2 _ Уь - слой ячеек с у-координатой центра Уг,кь-1 _ Уь - ку/2, для верхней границы Уг,кь+1/2 _ Уг - слой ячеек с у-координатой центра Уг,кг+1 _ Уг + ку/2. Значения параметров в фиктивных ячейках, а также на их гранях, совпадающих с границей расчетной области, вычисляются с использованием соответствующих этой границе краевых условий, а также значений параметров в смежных с фиктивными ячейками реальных ячейках расчетной области.
п
п
1
1. Жесткая стенка. Пусть жесткой стенкой является левая граница xjl_i/2kk = = xi. Тогда параметры, необходимые для расчета потока через грани ячеек, определяются в фиктивных ячейках при t = tn и t = tn+1/2 следующим образом:
uil-i,h = 2U(t,Vk) - uTl,k, vIl_i,fc = 2V(t,Vk) - vIbk,
(axx) Il_1,k = (axx)Il,k , (Sxy )Il_1,k = (Sxy )Il,k-
а параметры на совпадающих с границей гранях ячеек при t = tn+1/2 имеют вид
UIl-l/2,k = U (t,Vk ^ vil-1/2,k = V (t,Vk ^
(axx)Il _1/2,k = [UIl_1/2R,k - U (t,Vk)] PC1 - (&xx)h_1/2R,k, (Sxy )Il_1/2,k = [vIl_1/2R,k - V (t, Vk ^ Pc2 + (Sxy )I_1/2R,k,
где величины с индексом R вычисляются по формуле (7).
2. Свободная поверхность. Пусть свободной поверхностью является верхняя граница Vi,Kt+1/2 = Vt. Тогда параметры, необходимые для расчета потока через грани ячеек, определяются в фиктивных ячейках при t = tn и t = tn+1/2 следующими выражениями:
= = 2 p(t,xi) + (&yy )i,Kt
Ui,Kt + 1 = Ui,Kt, vi,Kt + 1 = vi,Kt - 2-,
PC1
(&yy )i,Kt + 1 = (°yy )i,Kt, (Sxy )i,Kt + 1 = 2t(t,xi) - (Sxy )i,Kt и на совпадающих с границей гранях ячеек при t = tn+1/2 следующими:
(°yy )i,Kt + 1/2 = -P(t,xi), (Sxy )i,Kt + 1/2 = T (t,xi),
T(t,xi) - (Sxy)i,Kt + 1/2B
■.,Kt + 1/2 = Ui,Kt + 1/2B + vi,j = vi,Kt + 1/2B
PC2
P(t,xi) + (°yy )i,Kt + 1/2B
PC1
где величины с индексом B вычисляются по формуле (10).
3. Искусственная граница. Пусть искусственной границей является нижняя сторона расчетной области Vi,Kb_1/2 = Vb. Как отмечалось выше, на искусственных границах ставятся условия отсутствия возмущений, приходящих к ним с их внешней стороны. В случае нижней границы эти условия математически можно представить следующим образом:
(I3)i,Kb_1/2 - (I3)b,0 = ° (h)i,Kb_1/2 - (I5)b,0 = 0,
где индексы b, 0 относятся к значениям параметров, соответствующим начальному (фоновому) состоянию тела со стороны нижней границы.
Отсюда параметры, необходимые для расчета потока через грани ячеек, определяются в фиктивных ячейках при t = tn и t = tn+1/2 по формулам
Ui,Kb_1/2 = 0.5 \ UiKb + Ub,o + vi,Kb_1/2 = 0.5 vi,Kb + vb,o +
(Sxy )i,Kb - (Sxy )b,0
PC2
(°yy)i,Kb - (&yy)b,0 PC1
(ауу)г,Кь-1/2 _ °-5 [РС1(Щ,Кь - + (ауу)г,Кь + (ауу)ь,0] > (Бху )г,Кь-1/2 _ °-5 [РС2(иг,Кь - иЬ,0) + (Бху )г,Кь + (Бху )Ь,0] • На искусственной границе при Ь _ Ьп+1/2 получаем
. (Бху )Кь - (Бху )Ь,0 , 0
иКь-1 _ иь,0 +--, VKь-1 _ -VКь + 2Vi кь-1/2,
РС2
(ауу )Кь-1 _ (ауу )Кь , (Бху )Кь-1 _ (Бху )Кь •
Граничные условия на других жестких стенках, свободных поверхностях и искусственных границах реализуются аналогичным образом.
Следует отметить, что приведенный способ аппроксимации искусственных границ имеет первый порядок точности в отличие от аппроксимаций жестких стенок и свободных поверхностей, которые имеют второй порядок. Такая точность на искусственных границах представляется достаточной, так как по своему смыслу такие границы отделяют несущественное дальнее поле тела от представляющего интерес ближнего поля, рассматриваемого как расчетная область, где динамика тела может определяться, например, воздействием на поверхности тела, которые могут быть как жесткими стенками, так и свободными поверхностями.
3. Расчет одномерных задач
3.1. Распространение волн от поверхности тела. Предполагается, что тело занимает полупространство у > 0. Его состояние при Ь _ 0 определяется заданными значениями инвариантов
11 _ ¡1,0, 12 _ ^2,0, ¡3 _ ¡3,0, ¡4 _ 14,0, ¡5 _ ¡5,0, ¡6 _ 1б,0,
где 13,0, 14,0, ¡5,0, 16,0 - некоторые константы.
Рассматриваются два случая распространения одномерных волн от границы тела У _ 0 .В первом случае граница является жесткой стенкой и
и(Ь, 0) _ иГ(Ь), V(Ь, 0) _ vг(t),
где иг(Ь) и vг(t) - заданные функции. Индекс Г здесь и далее указывает на отнесение к границе. Во втором случае граница тела У _ 0 является свободной поверхностью и
а у (Ь, 0)_ -рг(Ь), Бху (Ь, 0)_ тг(Ь),
где рг(Ь) и тг(Ь) - заданные функции.
Распространение волн в обоих случаях описывается в терминах инвариантов следующим образом:
У > С1Ь : 13 _ ¡3,0, 14 _ ¡4,0, 15 _ ¡5,0, 16 _ 16,0•
С2Ь <У < С1Ь : ¡з _ ¡з,г(Ь - у/с1), ¡4 _ ¡4,0, 15 _ ¡5,0, ¡6 _ ¡6,0,
0 <у < С1Ь : ¡з _ ¡з,г(Ь - у/с\), ¡4 _ ¡4,0, ¡5 _ ¡5,г(Ь - У/С2), ¡6 _ ¡6,0,
где граничные значения инвариантов задаются в случае жесткой стенки выражениями
¡з,г _ -!4,0 + 2vг, ¡5,г _ ¡6,0 + 2иг, а в случае свободной поверхности - выражениями
Т _ Т +2РГ Т _ Т +2ТТ
1 3,г _ Т4 0 +--, т5,г _ Т6,0 +--•
РС1 РС2
Рис. 1. Распределения всестороннего давления P в теле в два момента времени ti и t2 при изменении скорости (a) и давления (b) на поверхности тела y = 0 по законам, указанным на вставках. Черные линии - аналитическое решение, красные - схема UNO, зеленые - классический метод С.К. Годунова
В качестве расчетной области принимается отрезок 0 < у < К, правая граница которого у = К является искусственной, на ней ставятся неотражающие условия. В рассматриваемых промежутках времени они сводятся к приравниванию значений всех необходимых для реализации граничных условий параметров значениям этих параметров в соответствующих приграничных ячейках.
На рис. 1 представлены результаты расчетов двух задач о распространении волн от поверхности тела у = 0 .В первой задаче (рис. 1, а) граница у = 0 является жесткой стенкой,
пг(Ь) = 0, vг(г) = V* > 0. Во второй задаче (рис. 1, Ь) граница у = 0 является свободной поверхностью,
(р* > 0, если 0 < г < г* рг(г) = { , тг(г) = 0.
[0, если г > г*
В обоих случаях в начальный момент времени все инварианты равны нулю, так что при г = 0
и(0, у) = v(0, у) = 0, ^(0, у) = Буу(0, у) = БХу(0, у) = Р(0, у) = 0.
В расчетах использовалась равномерная сетка из 100 ячеек.
На рис. 1 видно, что в обеих задачах схема иМО дает результат, значительно более близкий к точному решению, чем метод С.К. Годунова, который довольно сильно «размазывает» скачки. Со временем из-за накопления вычислительных погрешностей отклонение численных решений от точных увеличивается. В результате во второй задаче метод С.К. Годунова к моменту г = г2 заметно «подрезает» экстремум волны (более чем на 15%).
3.2. Взаимодействие волн с границей. Рассматриваемые в данном разделе задачи и их численная реализация отличаются от задач предыдущего раздела главным образом тем, что в начальный момент времени г = 0 в теле в направлении его границы у = 0 распространяются волны, определяемые начальным распределением инвариантов
/4 = ^4,0 (у), 1б = 1б,о(у).
Здесь /4 0(у), /б,о(у) - заданные функции, определяющие профили волн, распространяющихся со скоростями —С1 и —С2 соответственно.
Рис. 2. Взаимодействие волны с жесткой стенкой и с волной, образованной на этой поверхности: (а) - начальное распределение всестороннего давления Р в теле и закон изменения скорости V на жесткой стенке; (Ь), (с) - распределения Р в моменты времени Ь\ и ¿2 соответственно, ¿2 > > 0. Стрелочки указывают направление движения волн. Обозначения кривых те же, что и на рис. 1
Тогда аналитическое решение данных задач в терминах инвариантов имеет в областях
у>С1 Ь (А), С2Ь < у < С1Ь (В), 0 <у < С2Ь (С)
следующий вид:
I3 — Is,0,
I4 — I4,0 (y + Cit), I
1з — 1з,Г (t - y/в! ), I4 — h,0 (y + c11), I I3 — h,r (t - y/ci), I4 — I4,0 (y + C11), I
— h,0, I6 — I6,0 (y + C2t),
— h,0, I6 — I6,0 (y + C2t),
— I5,r(t - y/c2), I6 — I6,0 (y + C2t).
На рис. 2-4 приведены результаты расчетов рассматриваемых задач при начальных данных, определяемых по формулам
P (0,y)— p(y), Sxy (0,y)— T (y), I10 — 0, I20 — 0, I30 — 0, I40 —0.
С учетом этого при t — 0 имеем
(0>y) — iA+r- p(y), Syy(0,y) — -wh-^'
u(0,y) — TM, v(0,y) — - ШЬ+Зп...
' PC2 ' pci (3A + 2-)
На рис. 2 приведены результаты расчетов для начальных данных
U > 0, если y>R/2, ty) — 0
yyJ \0, если y < R/2,
Рис. 3. Взаимодействие волны с подвижной жесткой стенкой и волной, образованной на этой стенке: (а) - начальное распределение скорости и в теле и закон ее изменения на поверхности тела; (Ь), (с) - распределения и в моменты времени и ¿2 соответственно, ¿2 > tl > 0. Стрелочки указывают направление движения волн. Обозначения кривых те же, что и на рис. 1
и граничных условий
„ I v* > 0, если t < t*,
ur (t)=0, vr 0 '
I 0, если t > t*.
В этой задаче реализуется взаимодействие волны с жесткой стенкой и с образованным на ней импульсом. На рис. 2 видно, что к моменту ti уже практически завершилось взаимодействие возбужденного на поверхности тела импульса с передним фронтом распространяющейся к поверхности тела волны. В промежутке ti < t < t2 указанный фронт достиг поверхности тела и отразился от нее. Кроме того, видно, что схема UNO намного лучше воспроизводит точное решение. Численное решение, полученное с помощью метода С.К. Годунова, значительно больше «размазывает» фронты волн. В результате этого начиная с некоторого времени возникает «подрезка» амплитуды возбужденного на поверхности тела импульса, которая со временем все более увеличивается (в момент t2 она составляет 47%).
На рис. 3 приведены результаты расчетов для начальных данных
|У > 0, если 0.5R <y < 0.8R,
Р(У) = т (y) = <
[0, если y < 0.5R или y > 0.8R
и граничных условий
(-п*(п* > 0), если t<t* ur (t) = ^ , vr (t)=0.
0, если t > t*
Такими условиями описывается взаимодействие волны с жесткой стенкой и образованным на ней импульсом. Из рис. 3 видно, что к моменту ti в результате
Рис. 4. Взаимодействие импульса с неподвижной свободной поверхностью: начальное распределение давления Р в теле (а) и его распределения в один из моментов до начала взаимодействия с поверхностью тела (Ь) и в один из моментов после начала взаимодействия (с). Стрелочки указывают направление движения волн. Обозначения кривых те же, что и на рис. 1
«воздействия» на поверхность тела в ее окрестности сформировалась волна, уходящая от поверхности тела со скоростью о2 навстречу распространяющейся в ее сторону также со скоростью С2 волне, определяемой начальными условиями. Через некоторое время эти волны начинают взаимодействовать. К моменту t2 их взаимодействие уже завершилось. Кроме того, к этому моменту передний фронт распространяющейся к поверхности тела волны уже отразился от нее. Видно, что и в этой задаче численное решение схемы UNO значительно ближе к точному. Численное решение, полученное с помощью метода С.К. Годунова, намного сильнее «размазывает» фронты волн и заметно «подрезает» амплитуду возбужденного на поверхности тела импульса. В момент Ь2 «подрезка» составляет около 33%.
На рис. 4 приведены результаты расчетов для начальных данных
P(y) = 0, т (y) =
и граничных условий
jp* > 0, если 0.5R <y < 0.8R,
\0, если y < 0.5R или y > 0.8R,
Pr(t) = 0, тг (t) = 0.
В этом случае моделируется взаимодействие волны с неподвижной свободной поверхностью. На рис. 4 видно, что взаимодействие начинается в интервале ti < t < t2. Рисунок показывает, что и в этой задаче схема UNO дает по сравнению с методом С.К. Годунова значительно более близкое к точному численное решение. В частности, и здесь в решении UNO-схемы нет накопленного к моменту t2 большого «размазывания» фронтов волны и «подрезки» ее амплитуды, что наблюдается в численном решении, полученном с использованием метода С.К.Годунова, где, в частности, подрезка составляет около 4%.
Рис. 5. Динамика импульсов внутри тела и их взаимодействие с неподвижной жесткой стенкой: начальные распределения скоростей и и V в теле (а) и их распределения в один из моментов до начала их взаимодействия со стенкой (Ь) ив один из моментов после начала взаимодействия (с). Стрелочки указывают направление движения волн. Обозначения кривых те же, что и на рис. 1
На рис. 5 представлены результаты расчетов задачи при следующих начальных
данных:
P (0,y ) = <j
Sxy (0,y) = u(0y) = -
PC2
p* > 0 при 0.67R <y < 0.93R, 0 при y < 0.67R или y > 0.93R,
iO.lp* при 0.33R <y < 0.53R, (0 при y < 0.33R или y > 0.53R,
v(0,y) = P^, Sxx(0,y) = 0, Syy(0,y) = 0
pci
и следующих граничных условиях
ur (t) = 0, vr (t) = 0.
С помощью этой задачи описываются распространение в теле продольной и поперечной волн и их взаимодействие с неподвижной жесткой стенкой. В начальный момент t = 0 продольная волна определяется профилем скорости v, а поперечная волна - профилем скорости u. Обе волны распространяются в направлении поверхности тела со скоростями c\ и c2 соответственно, c\ > c2 .К моменту t\ продольная волна догоняет поперечную. К моменту t2 продольная волна уже отразилась от поверхности тела, тогда как поперечная волна еще находится в процессе взаимодействия с поверхностью. Видно, что и здесь, как и в предыдущих задачах, схема UNO значительно меньше «размазывает» фронты обеих волн и лучше описывает их взаимодействие с поверхностью тела.
Рис. 6. Удар по части поверхности тела (область удара закрашена черным): (a) - изолинии ai, вычисленные на сетке 100 х 100 методом С.К. Годунова (черные кривые) и UNO (красные); (b) - изолинии ai, вычисленные методами С.К. Годунова на сетке 800 х 800 (черные кривые) и UNO - на сетке 100 х 100 (красные); (c) - распределения ai вдоль плоскости симметрии x = 0, вычисленные методом UNO на сетках 100 х 100, 200 х 200, 400 х 400, 800 х 800 и 1600 х 1600 (результаты на трех последних сетках графически совпадают)
4. Расчет двумерной задачи
Для иллюстрации работоспособности и эффективности UNO-схемы в двумерном случае рассматривается следующая задача удара по поверхности тела. При t < 0 тело, представляющее собой, как и ранее, упругое полупространство y > 0, находится в покое при
u(0, y) = «(0, y) = 0, Sxx(0, y) = Syy(0, y) = SXy(0, y) = P(0, y) = 0. Поверхность тела является свободной, при этом
ay(t, 0) = 0, Sxy (t, 0) = 0.
В начальный момент времени по части поверхности тела осуществляется ударное воздействие, реализующееся в повышении давления до p* > 0, что в граничном условии выражается следующим образом
ay (t, 0) = p* при - 1/3R < x < 1/3R.
Данная задача имеет плоскость симметрии x = 0. C учетом этого в качестве расчетной области принимается квадрат [0 < x < R] х [-R < y < 0], верхняя граница которого y = 0 является свободной поверхностью, левая x = 0 - неподвижной жесткой стенкой, а правая x = R и нижняя y = —R - искусственными границами.
Работоспособность и эффективность UNO-схемы иллюстрирует рис. 6, где приведен фрагмент расчетной области с изолиниями интенсивности напряжений ai и распределение ai в этом фрагменте вдоль плоскости симметрии задачи x = 0. Интенсивность напряжений ai рассчитывается по формуле
ai = \J3(aX + a1 + axaz + Tly).
Видно, что рассматриваемая UNO-схема намного экономичнее метода С.К. Годунова. Действительно, метод С.К. Годунова позволяет получить численное решение, аналогичное рассчитанному по UNO-схеме на сетке 100 х 100, лишь на сетке
800 х 800. Это означает, что при решении двумерных задач рассматриваемая UNO-схема более чем в 64 раза эффективнее метода С.К. Годунова.
Рис. 6, c демонстрирует быструю сходимость результатов расчетов по схеме UNO при измельчении сетки.
Заключение
Представлены основные положения методики расчета динамики линейных волн в упругом теле и их взаимодействия с поверхностью тела на основе UNO-модификации классического метода С.К. Годунова, имеющей второй порядок точности по пространству и времени. В отличие от TVD-модификаций, где TVD-ограничение выполняется строго, в данной модификации это ограничение выполняется приближенно с точностью до аппроксимации уравнений динамики волн. В результате второй порядок аппроксимации сохраняется в том числе и в областях локальных экстремумов решения.
Исследованы работоспособность и эффективность предложенной методики на ряде одномерных задач, имеющих точное аналитическое решение. Это задачи о распространении в теле возмущений, образованных на его поверхности, задачи о взаимодействии возмущений с поверхностью тела, а также задачи о взаимодействии возмущений в теле между собой и с поверхностью тела. Поверхность тела моделировалась и как подвижная и неподвижная жесткие стенки, и как свободная поверхность. Наряду с точными решениями, для оценки работоспособности и эффективности предлагаемой методики применялось также и сравнение с результатами расчетов классическим методом С.К. Годунова. Показано, что на аналогичных сетках UNO-схема дает результаты намного ближе к точным решениям. В частности, ширина размазывания резких фронтов волн при использовании UNO-схемы намного меньше. В результате к концу рассмотренных интервалов времени метод С.К. Годунова, как правило, приводил к значительной «подрезке» амплитуды волн, тогда как в численных решениях UNO-схемы этого не наблюдалось.
Работоспособность и эффективность предлагаемой методики в двумерных задачах проиллюстрированы на примере динамики приповерхностного слоя тела при ударном воздействии на часть его поверхности. Показано, что по сравнению с методом С.К. Годунова рассматриваемая UNO-схема позволяет снизить затраты компьютерных ресурсов на расчет этой задачи более чем на порядок (примерно в 60 раз).
Литература
1. Годунов С.К., Рябенький В.С. Разностные схемы. - М.: Наука, 1973. - 400 с.
2. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. - М.: Наука, 1976. - 400 с.
3. Aganin A.A. Dynamics of a small bubble in a compressible fluid // Int. J. Numer. Methods Fluids. - 2000. - V. 33, No 2. - P. 157-174. - doi: 10.1002/(SICI)1097-0363(20000530)33:2<157::AID-FLD6>3.0.C0;2-A.
4. Чебан В.Г., Навал И.К., Сабодаш П.Ф., Чередниченко Р.А. Численные методы решения задач динамической теории упругости. - Кишинев: Штиинца, 1976. - 226 с.
5. Иванов Г.В., Волчков Ю.М., Богульский И.О., Кургузов В.Д., Анисимов С.А. Численное решение динамических задач упругопластического деформирования твердых тел. - Новосибирск: Сиб. унив. изд-во, 2002. - 352 с.
6. Аганин А.А., Хисматуллина Н.А. Ударное воздействие струи жидкости на упруго-пластическое тело // Учен. зап. Казан. ун-та. Сер. физ.-матем. науки - 2014. - Т. 156, кн. 2. - С. 72-86.
7. Курант Р., Фридрихе К., Леви Г. О разностных уравнениях математической физики // Усп. матем. наук. - 1941. - Вып. 8. - С. 125-160.
8. Harten A., Engquist B, Osher S., Chakravarthy S.R. Uniformly high order accurate essentially non-oscillatory schemes, III //J. Comp. Phys. - 1987. - V. 71, No 2. -P. 231-303. - doi: 10.1016/0021-9991(87)90031-3.
9. Аганин А.А., Ильгамов М.А., Халитова Т.Ф. Моделирование сильного сжатия газовой полости в жидкости // Матем. моделирование. - 2008. - Т. 20, № 11. - С. 89-103.
10. Аганин А.А., Халитова Т.Ф., Хисматуллина Н.А. Расчет сильного сжатия сферического парогазового пузырька в жидкости // Вычисл. технологии. - 2008. - Т. 13, № 6. - С. 54-64.
11. Аганин А.А., Халитова Т. Ф., Хисматуллина Н.А. Метод численного решения задач сильного сжатия несферического кавитационного пузырька // Вычисл. технологии. -2010. - Т. 15, № 1. - С. 14-32.
12. Harten A. High Resolution schemes for hyperbolic conservation laws //J. Comput. Phys. - 1983. - V. 49, No 3. - P. 357-393. - doi: 10.1016/0021-9991(83)90136-5.
13. Liu X.-D., Osher S., Chan T. Weighted essentially non-oscillatory schemes // J. Comput. Phys. - 1994. - V. 115, No 1. - P. 200-212
14. Lax P.D., Wendroff B. Systems of conservation laws // Comm. Pure Appl. Math. - 1960. -V. 13, No 2. - P. 217-237. - doi: 10.1002/cpa.3160130205.
15. Ильгамов М.А., Гильманов А.Н. Неотражающие условия на границах расчетной области. - М.: ФИЗМАТЛИТ, 2003. - 240 с.
Поступила в редакцию 15.03.17
Аганин Александр Алексеевич, доктор физико-математических наук, профессор, заведующий лабораторией
Институт механики и машиностроения КазНЦ РАН
ул. Лобачевского, д. 2/31, г. Казань, 420111, Россия E-mail: [email protected]
Хисматуллина Наиля Абдулхаевна, кандидат физико-математических наук, старший научный сотрудник
Институт механики и машиностроения КазНЦ РАН
ул. Лобачевского, д. 2/31, г. Казань, 420111, Россия E-mail: [email protected]
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online) UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)
2017, vol. 159, no. 2, pp. 143-160
Computation of Two-Dimensional Disturbances in an Elastic Body
A.A. Aganin* , Khismatullina N.A. **
Institute of Mechanics and Engineering, Kazan Science Center, Russian Academy of Sciences, Kazan, 420111 Russia E-mail: *[email protected], **[email protected]
Received March 15, 2017 Abstract
The classical Godunov method is widely used for the numerical study of waves in continuous media. If the Courant condition is satisfied, the Godunov scheme is stable and monotonous. However, due to its first order of accuracy it can lead to rather large smearing of jumps, contact discontinuities, and other features of the solution in the domains where the solution gradients are great.
In this paper, the possibility of increasing the efficiency of computation of linear waves in an elastic body by applying one of the second-order accurate UNO modifications of the classical Godunov method has been studied. In that UNO modification, the monotonicity requirement is replaced by the TVD condition, and all the parameters inside each grid cell are assumed linear rather than constant as they are in the Godunov method. The TVD condition is met just on the level of approximation. To derive the second order of accuracy in time, the time derivatives of the unknown functions have been expressed in terms of their spatial derivatives. Those expressions allow to calculate the next half-time-layer values of the unknown functions at the center of the grid cells and on both sides of the cell boundaries. The values of the unknown functions on the boundaries themselves have been found by solving the corresponding Riemann problems. Subsequently, numerical flows across the cell boundaries have been computed. Computation of the values at the next time layer has been carried out by an explicit scheme. Some limiters for the values of the spatial derivatives have been introduced. The present UNO scheme limiters use approximations to both the first and second derivatives.
The efficiency of the proposed UNO modification has been estimated by computing a number of one- and two-dimensional problems on propagation of linear waves in an elastic body and their interaction with each other and with the surface of the body. The results of those computations have been compared with the exact solutions and the results of applying the Godunov method. It has been shown that the UNO scheme considered allows one to reduce computational costs by more than a factor of ten.
Keywords: UNO scheme, Godunov scheme, efficiency of difference schemes, linear wave, elastic body
Figure Captions
Fig. 1. Distributions of the uniform pressure P in the body at two moments of time t1 and t2 upon changes in the velocity (a) and pressure (b) at the body surface y = 0 following the laws on the insertions. Black lines - analytical solution, red lines - UNO scheme, green lines - classical Godunov method.
Fig. 2. The interaction of the wave with the rigid wall and with the wave, which is formed at this surface: (a) - initial distribution of the uniform pressure P in the body and the velocity
defect law v on the rigid wall; (b), (c) - distributions of P at the moments of time ti and t2 , respectively, t2 > t1 > 0. Arrows show the direction in which the waves move. The curves are designated in the same way as in Fig. 1.
Fig. 3. The interaction of the wave with the mobile rigid wall and the wave which is formed in this wall: (a) - initial distribution of the velocity u in the body and the law of its defect at the body surface; (b), (c) - distributions u at the moments of time t1 and t2 , respectively, t2 > t1 > 0. Arrows show the direction in which the waves move. The curves are designated in the same way as in Fig. 1.
Fig. 4. The interaction of the wave with the mobile rigid wall and the wave which is formed on this wall: (a) - initial distribution of the velocity u in the body and the law of its defect at the body surface; (b), (c) - distributions u at the moments of time t1 and t2 , respectively, t2 > t1 > 0. Arrows show the direction in which the waves move. The curves are designated in the same way as in Fig. 1.
Fig. 5. The dynamics of impulses inside the body and their interaction with the non-mobile rigid wall: initial distributions of the velocities u and v in the body (a) and their distributions at one of the moments before their interaction with the wall (b) and at one of the moments after the onset of interaction with (c). Arrows show the direction in which the waves move. The curves are designated in the same way as in Fig. 1.
Fig. 6. Impact on the part of the body surface (the area of the impact is colored in black): (a) - isolines oi, computation on the grid 100 x 100 by the Godunov (black curves) and UNO (red curves) methods; (b) - isolines oi, computed by the Godunov method on the grid 800 x 800 (black curves) and by the UNO method on the grid 100 x 100 (red curves); (c) -distributions oi along the symmetry plane x = 0 calculated by the UNO method on the grids 100 x 100, 200 x 200, 400 x 400, 800 x 800, and 1600 x 1600 (the results on the latter three grids coincide graphically).
References
1. Godunov S.K., Ryaben'kii V.S. Difference Schemes. Moscow, Nauka, 1973. 400 p. (InRus-sian)
2. Godunov S.K., Zabrodin A.V., Ivanov M.Ya., Kraiko A.N., Prokopov G.P. Numerical Solution of Multidimensional Problems of Gas Dynamics. Moscow, Nauka, 1976. 400 p. (In Russian)
3. Aganin A.A. Dynamics of a small bubble in a compressible fluid. Int. J. Nu-mer. Methods Fluids, 2000, vol. 33, no. 2, pp. 157-174. doi: 10.1002/(SICI)1097-0363(20000530)33:2¡157::AID-FLD6¿3.0.CO;2-A.
4. Cheban V.G., Naval I.K., Sabodash P.F., Cherednichenko R.A. Numerical Methods for Solving Problems of Dynamic Elasticity Theory. Chisinau, Shtiintsa, 1976. 226 p. (In Russian)
5. Ivanov G.V., Volchkov Yu.M., Bogul'skii I.O., Kurguzov V.D., Anisimov S.A. Numerical Solution of the Dynamic Problems of Elastoplastic Deformation of Solid Bodies. Novosibirsk, Sib. Univ., 2002. 352 p. (In Russian)
6. Aganin A.A., Khismatullina N.A. Liquid jet impact on an elastic-plastic body. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, 2014, vol. 156, no. 2, pp. 72-86. (In Russian)
7. Courant R., Friedrichs K., Lewy H. On difference equations of mathematical physics. Usp. Mat. Nauk, 1941, no. 8, pp. 125-160. (In Russian)
8. Harten A., Engquist B., Osher S., Chakravarthy S.R. Uniformly high order accurate essentially non-oscillatory schemes, III. J. Comp. Phys., 1987, vol. 71, no. 2, pp. 231-303. doi: 10.1016/0021-9991(87)90031-3.
9. Aganin A.A., Il'gamov M.A., Khalitova T.F. Simulation of strong compression of a gas cavity in liquid. Mat. Model., 2008, vol. 20, no. 11, pp. 89-103. (In Russian)
10. Aganin A.A., Khalitova T.F., Khismatullina N.A. Calculation of strong compression of a nonspherical cavitation bubble in fluid. Vychisl. Tekh., 2008, vol. 13, no. 6, pp. 54-64. (In Russian)
11. Aganin A.A., Khalitova T.F., Khismatullina N.A. The method of numerical solution of strong compression of a nonspherical cavitation bubble. Vychisl. Tekh., 2010, vol. 15, no. 1, pp. 14-32. (In Russian)
12. Harten A. High resolution schemes for hyperbolic conservation laws. J. Comput. Phys., 1983, vol. 49, no. 3, pp. 357-393. doi: 10.1016/0021-9991(83)90136-5.
13. Liu X.-D., Osher S., Chan T. Weighted essentially non-oscillatory schemes. J. Comput. Phys., 1994, vol. 115, no. 1, pp. 200-212
14. Lax P.D., Wendroff B. Systems of conservation laws. Commun. Pure Appl. Math., 1960, vol. 13, no. 2, pp. 217-237. doi: 10.1002/cpa.3160130205.
15. Il'gamov M.A., Gil'manov A.N. Non-Reflecting Conditions on Computational Domain Borders. Moscow, FIZMATLIT, 2003. 240 p. (In Russian)
/ Для цитирования: Аганин А.А., Хисматуллина Н.А. Расчет двумерных возму-( щений в упругом теле // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. - 2017. -\ Т. 159, кн. 2. - С. 143-160.
/ For citation: Aganin A.A., Khismatullina N.A. Computation of two-dimensional ( disturbances in an elastic body. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-\ Matematicheskie Nauki, 2017, vol. 159, no. 2, pp. 143-160. (In Russian)
