Введение
Задача определения нестационарных волновых полей относится к числу весьма сложных. Отличительной чертой многих практически важных динамических задач механики деформируемого твердого тела является разрывной характер их решений. Таковы случаи, когда рассматриваемые упругие области конечных размеров содержат разрывы в граничных условиях, отверстия с ломаным контуром, вырезы и инородные включения, являющиеся источниками высокой концентрации напряжений. Решение краевых задач для таких областей невозможно без разработки эффективных численных методов. Поэтому центр тяжести проблемы исследования нестационарных волновых движений в сплошных средах все более смещается в сторону разработки и совершенствования разностных схем, позволяющих улучшать результаты расчетов. Продолжают интенсивно разрабатываться приближенные методы в различных модификациях, используемых для численного решения динамических задач однородных и структурно-неоднородных сред, как-то: методы конечных разностей, дробных шагов, пространственных характеристик, конечных элементов, сеточно-характеристический метод Годунова, граничных интегральных уравнений, метод источников и др.
Известными преимуществами обладают конечно-разностные методы, основанные на использовании характеристических поверхностей и соотношений совместности на них. Они максимально приближают область зависимости разностных и исходных дифференциальных уравнений, обеспечивая высокую точность результата как в случае гладких, так и разрывных решений, позволяют корректно рассчитывать границы тела и контактные поверхности. В 1960 г. была предложена явная разностная схема второго порядка точности для гиперболической системы дифференциальных уравнений в частных производных первого порядка относительно трех переменных [6]. Эта схема, основанная на характеристическом подходе, реализована для исследования плоских упругих волн [7]. В последуюшем метод пространственных характеристик развит для решения конкретных задач динамики сплошных сред [1-5].
Постановка задачи. Исследуется плоская деформация упругого тела с прямоугольным поперечным сечением. Сечение имеет размеры 0 ≤ x ≤ l, - L ≤ x 2 ≤ L (рисунок 1). Стороны прямоугольника разбиты соответственно на n1 и n2 частей. При этом шаги разбиения по координатам определены равенствами h 1 = l / n 1 и h 2 = L / n 2. Таким образом, произвольная узловая точка Оi j будет иметь координаты ( x 1i, x 2j). При этом x 1i = i h 1 (i = 0, 1, 2, ... , n 1) и x2j = j h 2 ( j = - n 2, - n 2+1, -n 2 +2, ..., -1, 0, 1, 2, ..., n 2 -1, n 2).
На части L* ≤ x 2 ≤ L** границы x 1 = 0 прямоугольника прикладывается внешняя нормальная нагрузка, распределенная по границе по П-образной форме,
p +q = f (x 2, t) = A sin (wt), 0 ≤ t ≤ t* (1)
τ = 0
Рисунок 1. Исследуемая область
изменяющаяся во времени t и постоянная по поперечной координате x2. Остальная часть границы x1 = 0 прямоугольника считается свободной от воздействия внешних напряжений. В (1) принято, что A - амплитуда внешней нагрузки, а w её частота. Нагрузка действует на ограниченном участке времени t*, определяя период его воздействия и соответствующую длину волны силы возбуждения. В моменты времени, превышающие t *, нагрузка на этом участке границы полностью снимается, т.е. считается
p + q = 0, t = 0 при t ≥ t*. (2)
Граница x 1 = l прямоугольника не нагружена и поэтому считается свободной от каких-либо воздействий, т.е.
p + q = 0, t = 0 при t ≥ 0. (3)
Наконец, границы x 2 = ± L прямоугольника предполагаются закрепленными и на них скорости перемещений равны нулю в любой момент времени, т.е.
v1 (x 1; t) = v2 (x 1; t) = 0 при t ≥ 0. (4)
Перечисленные граничные условия (1) - (4) должны быть дополнены начальными условиями. Предполагается, что в начальный момент времени (t = 0) тело не нагружено и находится в состоянии покоя, т.е.
v1 (x1; x 2; 0) = v 2 (x1; x 2; 0) = p ( x1; x 2; 0) = q (x1; x 2; 0)= τ (x1; x 2; 0) = 0. (5)
Задача заключается в определении внутри прямоугольника полей напряжений и скоростей, вызванных фронтами падающих и многократно дифрагированных упругих волн в момент времени t > 0 .
В условиях плоской деформации волновой процесс во внутренних точках прямоугольника описывается системой динамических уравнений гиперболического типа, содержащей в качестве неизвестных безразмерные напряжения p , q , τ скорости перемещений v1 , v 2
v1, t - p , 1 - q , 1 - t , 2 = 0 ; v2, t - p , 2 + q , 2 - t , 1 = 0 ; (6)
γ 2 ( γ 2 - 1 ) -1 p , t - v1 ,1 - v2 , 2= 0 ; γ 2 q , t - v1 ,1+ v2 , 2 = 0 ;
γ 2 t , t - v1 , 2 - v2 , 1 = 0 .
Следуя работе [9], введены безразмерные независимые переменные и искомые величины:
(7)
где b - характерная длина. В дальнейшем черта над безразмерными параметрами опускается.
При построении численного решения предполагается, что граница прямоугольника совпадает с линией узлов квадратной сетки, которая покрывает исследуемую область.
Таким образом, необходимо найти решение уравнений (6) при сформулированных условиях (1) - (5). Поставленная задача решена методом пространственных характеристик, подробный алгоритм численной реализации которого изложен в [9]. Однако они могут быть использованы только в областях с непрерывным изменением всех входящих параметров. В связи с этим, ниже разработан и приведен алгоритм решения динамических задач в особых точках x 2 = L* и x 2 = L** границы x 1 = 0, в которых входящие параметры терпят разрыв первого рода.
Исследуемая область представлена на рисунке 1. Для определенности ниже рассматривается точка Е1(x 2 = L**) границы x1 = 0 прямоугольника. Окрестности точки Е1 рассматриваются как две угловые точки I и II. Для углов I и II выписываются соответствующие конечно-разностные соотношения, полученные в результате интегрирования по бихарактеристикам и по оси характеристических конусов. Нетрудно видеть, что соотношения для угла I подобны соотношениям правого верхнего угла R рассматриваемой области
δv 1I - δv 2 I + α8 δpI = А1,
δv 1 I + δv 2 I + α2 δqI = А2 (8а)
а для угла II - левого верхнего угла М
δv 1 II + δv 2 II + α8 δpII = А3,
δv 1 II - δv 2 II + α2 δqII = А4. (8б)
Здесь А1, А2, А3, А4 известные на нижнем слое величины.
Слева от точки Е1 всюду и в самой точке Е1 в соответствии с (1) заданы нормальные напряжения. Их приращения в обозначениях (7) можно записать в виде
δpI + δqI = A [ sin (ωt) - sin (ω(t-k) ) ]. (9)
Кроме того, необходимо удовлетворить условиям непрерывности скоростей перемещений и нормальных и касательных напряжений при переходе от точек одного угла к сопряженным точкам другого. При этом
δv1 I = δv1 II, δv2 I = δv2 II,
δpI - δqI = δpII - δqII, δtI = δtII . (10)
Выписанная система уравнений (8), (9) и (10) позволяет единственным образом определить приращения скоростей перемещений δv1 I, δv1 II, δv2 I, δv2 II и приращения напряжений δpI, δqI, δtI, δpII, δqII, δtII в точке Е1 разрыва граничных условий
δv1 = Δ1 /Δ , δv2 = Δ2 /Δ, δpI = Δ3 /Δ, δqI = Δ4 /Δ, dt = 0,
δpII = Δ 5 /Δ , δqII = Δ 6 /Δ. (11)
Введенные здесь определители вычисляются на нижнем слое.
Расчетные формулы (11) используются при вычислениях искомых величин в особой точке Е1. В "левой" особой точке Е2 (x 2 = L* ) границы x1 = 0 прямоугольника точно также устанавливаются расчетные соотношения.
Настоящие расчеты проведены для прямоугольной области 0 ≤ x1 ≤ 100 h1 и | x 2 | ≤ 100 h 2. При этом h1 = h 2 = h = 0.05. Шаг по времени k выбран в соответствии с необходимыми условиями устойчивости [9]
(12)
используемой явной конечно-разностной расчетной схемы. В расчетах он считался равным k = 0.025. Коэффициент А в равенстве (1) принят равным единице, а период приложенной импульсной нагрузки был выбран равным Т = 100k. Таким образом, круговая частота динамической нагрузки w принята равной ω = π / Т = π / 100 k. Область приложения нагрузки принята следующая: x1 = 0 , - 97 h ≤ x 2 ≤ 97 h. Симметричность приложенной нагрузки и характер деформирования исследуемого тела позволяют ограничиться анализом динамических явлений в области только положительных значений переменной x 2 ≥ 0.
На рис. 2 представлена осциллограмма продольной скорости v1 в пяти точках «наблюдения» x2 = 0 (точка наблюдения 1), x2 = 20 h (2), x2 = 40 h (3), x2 = 60 h (4) и x2 = 80 h (5), расположенных на лицевой границе x1 = 0 исследуемого тела. Во всех точках «наблюдения» на первых 100 шагах по времени вид осциллограммы определяется формой заданного на границе x1 = 0 синусоидального импульса (1). Первоначально появляется некоторое искажение в характере поведения параметра (это влияние дифрагированных из точки Е1(x1 = 0,x2 = 97h ) возмущений, распространяющихся со скоростью с1 продольных волн). Это искажение первым наблюдается в точке x2 = 80 h(5), которая является ближайшей к особой точке Е1. Через определенные моменты времени эффект особой граничной точки сказывается последовательно в точках «наблюдения» 4, 3, 2 и 1. Интенсивность влияния особой точки Е1 мала. Более мощное возбуждение обусловлено приходом в точки «наблюдения» дифрагированной из угловой точки R(x1 = 0,x2 = 100 h) волны, распространяющейся со скоростью продольной волны с1 = 1.Обсуждаемая волна также последовательно проходит через точки 5, 4, 3, 2 и 1. К моменту прихода в точку 1 дифрагированных из угловых точек M(x1 = 0 ,x 2 = -100 h ) и R(x1 = 0, x 2 = 100 h ) волн в неё приходит и отраженная от границы x1 =100 h прямоугольника волна. Поэтому на приведенной осциллограмме продольной v1 скорости перемещений в моменты времени, превышающие 400k, можно видеть результат интерференции волн, являющийся наложением волн различного типа.
Таким образом получены расчетные конечно-разностные соотношения динамических задач в особых точках лицевой границы прямоугольника, в которых искомые параметры терпят разрыв первого рода. Путем численной реализации установлена устойчивость расчетных алгоритмов для достаточно большего времени. В результате проведенных исследований можно заключить, что разработанная методика расчета применительно к нестационарным динамическим задачам достаточно правильно передает основные закономерности и особенности протекающих волновых процессов и позволяет проводить исследование напряженно-деформированного состояния в однородных и слоисто-неоднородных средах.
Рецензенты:
- Дасибеков А., доктор технических наук, профессор кафедры прикладной механики Южно- Казахстанского государственного университета имени М. Ауезова, г. Шымкент.
- Нысанов Е. А., доктор физико-математических наук, доцент кафедры «Теория и методика преподавания информатики» Южно-Казахстанского государственного университета имени М. Ауезова, г. Шымкент.
- Антонов Александр Владимирович, д.т.н., профессор, декан факультета кибернетики, Обнинский институт атомной энергетики Национального исследовательского ядерного университета МИФИ Министерства образования и науки РФ, г. Обнинск.