CC BY
38
вестник 812012
УДК 539.3
В.В. Немчинов
ФГБОУВПО «МГСУ»
ДВУХСЛОЙНАЯ РАЗНОСТНАЯ СХЕМА ЧИСЛЕННОГО РЕШЕНИЯ ПЛОСКИХ ДИНАМИЧЕСКИХ ЗАДАЧ
ТЕОРИИ УПРУГОСТИ
Разработан численный метод решения плоских динамических задач теории упругости, использующий конечные элементы для аппроксимации расчетных областей сложной формы, вычисляющий скорости и напряжения в среде и с высокой точностью выполняющий заданные граничные условия.
Ключевые слова: разностная производная, двухслойная разностная схема, дифракция продольной волны, окружность, контурные напряжения, численное моделирование, метод конечных элементов.
Безразмерная система дифференциальных уравнений первого порядка, описывающая распространение плоских упругих волн, имеет вид
д и д р д г 1Ч
— = — +—; (1.1) д , д х д у
5у = 5г + д1; (1.2)
д , д х д у
д р д и д V
-Ё- = — + а— ; (1.3)
д t д х д у
д а д и д V „
— = а— + —; (1.4)
д t д х д у
дт д t И
д и + д VI (15)
д у д х /
где и, V — безразмерные скорости по координате х, у соответственно; —
напряжения стхх, стуу, ст^ с множителем (р е2р); t — безразмерное время, связанное
с размерным временем множителем ср1Ь; р — плотность среды; £ — характерный
линейный размер; а = 1 -2р, Р = (/Ср) , ср, с8 — скорость продольной и поперечной волны.
Аппроксимируем данную систему уравнений по времени неявной разностной схемой:
дР° + дг0 (21) и, =--1--; (2.1)
дг0 дд0
V, =— +—; (2.2)
дх ду
ди0 ду0
Р> + ; (2.3)
дх ду ди0 ду0
д, = а-+-; (2.4)
г -р| ^ + ^ ' ду дх
(2.5)
(3)
где, например,
0 (() + и(( + А/)) и + и ~2 = 2 '
Производные по времени аппроксимируем разностной производной вперед, рассматривая скорость по оси х:
и - и
А'
(4)
Остальные скорости и безразмерные напряжения определяются аналогично. Схему (2) можно преобразовать в явный вид, используя разложение в ряд Тейло-ра по времени до первой производной, например, для скорости и :
О и + и и+ и + At и А'. и —---= иН--и .
(5)
2 2 2 Соотношения (2.1) умножаем на функцию формы N и интегрируем по площади элемента, для первого уравнения получаем
dp 0 dr 0 Jut -Ni -ds= ¡-^—Nids + J-N ids.
s s dx s дУ
Меняем порядок дифференцирования:
Г »г 1 го dNi , Г о dNi , \ut- Nt ■ ds + p —-ds + \r —l-ds = j j dx J dy
(6)
=J
dx
s
d(°Ni) d(r°N i)"
dx
dy
ds =
J(dy
- r°dx)Ni dl.
(7)
Представляем скорости и напряжения внутри КЭ через функции формы:
О 0 /Лг О О /лр 0 0 /Лг 0 0 /Лг 0 0 ¡\т ^
и =и ■'М у, V -V ■> Лу, р = р , q =¿7 ■'М/, г =г -1 N /. (8)
Для разностной производной по времени задаем аналогичные соотношения:
U' — u'Nj , V' = v'Nj , р' = p'Nj , q' = q'Nj , Г' = r/Nj . Имеем для скорости и (остальные аналогичны)
и° — и + 0,5 • d' • и. — и + 0,5 d' | — + — I —
— Ни + 0,5 d'
( dN.
дх
jpj +
дх ду
dN ^ _j
ду
(9)
(10)
Окончательно получаем
J N ■ N ■ ds I uj +I J^ ds I p! +
dx
I сdN. dN, Л . ((dN. dN,
0,5dt I f—^—J- ds IuJ + 0,5dt a I J—L—L I J dx dx ) I s dx dy
ds Iv1 +
(rdN- Л ■ (edN- dN, Л .
+ I f—N. ds I rJ + 0,5 dt pi —1 ds I v^ +
[I dy 1 ) dy dx )
+ 0,5dtplf——ds Iu1 = <J(p°dy - r0dx)N dl .
[ Js dy dy ) ^ 7
(11)
иt =
После применения формулы Грина в правой части получаем контурный интеграл по границе элемента.
Для получения глобальных соотношений выражение (11) и все аналогичные необходимо просуммировать по всем элементам расчетной области.
Для внутренних узлов расчетной области при суммировании контурные интегралы сокращаются, так как они вычисляются по границе элемента с учетом направления обхода (против часовой стрелки).
Обозначим матрицы:
,-дМ- гдМ- гдМ- дМ,
А, =1 -¿N,4. В, ^N,„,0,, ,,
в, ^ , = М^ , =£М ,,
5 дх ду ду дх 5 ду ду
М, = £ N , N } (к.
Получаем:
+ А,р< + В,г< + 0,5 Л [(С, ) +(а В, + Щ ) ] = = £(р0Су - г0ах)м,С1,
о
М, V + А, г + В, Ц +0,5 С[(ро, + Е,У + + (рв, + аЕ,)] =
= £(г0 Су - ц0 Ах )МС1,
о
М,р + А,и1 +аВ,- у' +0,5 Л [С, р1 + (в, + аЕ, ) + аЕ,Ц ] =
= £(и0 Су -ау0 Сх >)М1С1,
о
М, ц' +а А, и + В, у] +0,5 Л [аС, р + (аВ, +Е, ) + Е,ц] ~] =
= £(а и0 Су - у0 Сх)М1С1,
о
М, г1 + Р А, у* + вВ, и + 0,5Сг[рЕ, р1 + (рС, + РЕ ,) + рв, Ц ] = Р$(у°Су - и0Сх)М;с11.
(12)
(13.1)
(13.2)
(13.3)
(13.4)
(13.5)
о
Разностная схема (13) будет явной, если матрицу масс взять в диагональном виде, например для 4-угольного элемента:
»г 5 я
М = 4 я« = т,
где 5 — площадь элемента; Я ^ — символ Кронекера; 1 — дискретная точка. Обозначим для первого и аналогично для остальных уравнений
Ей = -(4 + Би т1 + 0,5А[( ) + (аВу + р^)) ]). (14)
Соотношения (15) примут вид
т и] = Еи + £ (р0Су - г0 Сх) ); (15.1)
о
т{ у] = Еи + £ (г0Су - ц0 Сх) ); (15.2)
т. р\ - Еы + ф (ы0¿у -а V0 ¿х) N; (15.3)
о
д; - Ед + ф (аы0¿у - V0 ¿х)); (15.4)
о
т.г. - Ег + ф( ¿у-вы0 ¿х)) . (15.5)
о
Соотношения (15) характерны тем, что необходимо выделять граничные точки и в (15) вычислять контурные интегралы.
Рассмотрим для примера 4-точечный элемент с базисными функциями (функции формы):
N,(5, „)==№„, N. (5, , (16)
N3 (5, „^р+Ш1^ , N4 (5, „)=№„>.
На границе элементы линейны, поэтому интегралы вычисляются легко и формулу (17) можно преобразовать к более простому виду. Для первого уравнения имеем
т1ы] - Еы + ф(р0 ¿у - г0 ¿х )) - Еы + ф(pdy - rdx)) +
о о (17)
+ 0,5 ¿г ф(р^у - г^х)N¡.
о
Так как граница расчетной области в нашем случае состоит из прямых линий, то все необходимые соотношения можно получить в направленных отрезках, не рассматривая тригонометрические соотношения. Получаем формулы, проектирующие скорости по х и у на оси п и т (внешнюю нормаль и касательную к границе КЭ):
dsU = -udх - уСу; (18.1)
dsV = udy - vdх; (18.2)
ds2 = сЬс2 + dy2. (18.3)
Имеем из выражений (18):
-и1 dх —у, dy = -ЖиСх -Жу dy +Ь \_-(рСу-г, dх)dх-( dy - qtdх)dy~^ =
= -ЖиСх -Жу СУ + Ь ^-ргСхСу+г^х2 - ггСу2 + qtdхdy^ = (19)
= - Жи Сх - Жу Су + ЬЯ1,
где Ь = 0,25 (Ь, + Ь2,; Ь,, Ь2 — стороны элементов, примыкающих к границе; К — напряжение по т;
Жы - - Еы + ф (pdy - Ых) N. - -Еы + 0,5^у1 + ¿у2)р. - 0,5 (¿х1 + ¿х2) г (20.1)
ь
результат суммирования граничных элементов в точке 7 .
Жv - - Еv +ф(у - дdx)N¡ ^^^^ + 0,5(^ + ¿у2)г1-0,5 (( + ¿х2); (20.2)
ь
Ся2К = ( - р*)СхСу + г (сСх2 - Су2);
и Су -у,Сх = ЖиСу -ЖуСх + Ь \_(р,Су - г,Сх )Су -(г,Су - qtdх ~)сСх ] = = Ж1^у-Ж\^х+Ь^р^у2 - rtdxdy - rtdydx + д^х2^ = ЖисЬ — ЖУСЬС + ЬОп где сЫ2С) = рс1х2 - 2гс/хс/у + с/с/у1 — напряжение по нормали п .
Ся2 £ = рСу2 + 2 гСхСу + qdх2; (20.3)
Сч2 Q = рСх2 - 2гСхСу + дСу2; (20.4)
Сх2Я = (д — р)хСу + г[Сх2 —Су2 ). (20.5)
Получаем 2 первых соотношения на границе:
т.П, = Ши + ЬЯ,; (21.1)
т 1 Г(=Ш+Ь()(. (21.2)
2 2 2
Если ввести напряжение по касательной к границе: Сх Б=рСх + 2гСхСу + дСу , то получаем оставшиеся 3 соотношения:
= ШБ + ЬаУ1; (21.3)
т£, =^ + ЬУ,; (21.4)
тД, = ШЯ + $Ьи,, (21.5)
где СчШи = -ШиСх-ШуСу; (22.1)
СчШУ = ШиСу-ШуСх ; (22.2)
ШрСх2 +2 ШгСхСу + ШдСу2- (22.3)
2ШШрСх2 - 2ШгСхСу + ШдСу2; (22.4)
= (Шд - Шр)сСХСу + Жг[х2 - Су1), (22.5)
Q — по нормали (к радиусу полярной системы координат); Б — по касательной (направлена по углу полярной системы координат).
Из соотношений (22) и (23) можно образовать 3 линейно независимые комбинации так, чтобы в знаменателе выражений не могло возникнуть деление на ноль:
1 ШЦ-щД;
V - а = Ы; (23.2)
т 1 + Ь
Ш - аЩ ...
Б,-аQt=--, (23.3)
т1
где т^ уже просуммированная масса в граничной точке; Ь = 0,25 (( + Ь2 — получено при суммировании контурных интегралов, содержащих ускорения.
Формулы (23) используются для точек границы и учитывают граничные условия задачи. Например, для свободной границы имеем:
Q = 0;
Я = 0.
Обратно в декартову систему координат (х,у) можно перейти по формулам:
Счи = - иСу + УСх; (24.1)
Счу = - иСх - УСу; (24.2)
Сч2р = БСх2 - 2ЯСхСу + QСy2; (24.3)
Сч2д = БСу2 + 2ЯСхСу + QСx2; (24.4)
С*2г = [-0)СхСу + я[сх2 -С/). (24.5)
В качестве тестовой задачи для исследования сходимости полученного численного метода рассмотрена известная задача о дифракции продольной волны на круглом отверстии единичного диаметра [1—6].
В табл. приведено изменение максимального граничного напряжения в зависимости от .
Сходимость численного решения (ненулевые граничные безразмерные напряжения) в зависимости от шага разбиения дискретной сетки
h ■ "min Число разбиений полуокружности (p < ) Точное решение Относительная ошибка, %
1. 0,098 16 -2,425 -2,97 18,4
2. 0,049 32 -2,65 -2,97 9,9
3. 0,0327 48 -2,711 -2,97 8,7
4. 0,0314 50 -2,737 -2,97 7,7
5. 0,0224 70 -2,74 -2,97 7,7
6. 0,01382 70 (hr=hr /2) -2,8116 -2,97 5,3
7. — -2,97 — —
Графики и^ (рср) и зависимости от времени показаны на рисунке для различных йт;п. Приведены безразмерные напряжения в граничной точке на диаметре окружности, параллельном фронту набегающей продольной волны.
Сравнение напряжений для различных Ат;п(табл.)
(р ср )
График показывает (рис.), что в отличие от точного решения напряжения на границе круга при уменьшении Ат;п изменяются практически линейно. 1рср)
Полученный численный метод сочетает достоинства МКЭ и метода пространственных характеристик Клифтона [1, 2, 5, 6] — точное выполнение граничных условий. Внутри расчетной области метод имеет точность порядка о (к2) согласно выводу, но на границе, как показывают расчеты, точность падает до о (к). Возможно, использование более точных КЭ, например, 8-точечного изопараметрического конечного элемента, сделает численное решение на границе более точным.
Библиографический список
1. Барон М.Л., Метьюс. Дифракция волны давления относительно цилиндрической полости в упругой среде // Прикладная механика. Серия Е. 1961. № 3. С. 31—38.
2. Клифтон Р.Дж. Разностный метод в плоских задачах динамической упругости // Механика. 1968. № 1 (107). С. 103—122.
3. Мусаев В.К. Применение метода конечных элементов к решению плоской нестационарной динамической задачи теории упругости // Механика твердого тела. 1980. № 1. С. 167.
4. Мусаев В.К Воздействие продольной ступенчатой волны на подкрепленное круглое отверстие в упругой среде // Всесоюз. конф. «Современные проблемы строительной механики и прочности летательных аппаратов» : тез. докл. М. : МАИ, 1983. С. 51.
ВЕСТНИК 8/2Q12
5. СабодашП.Ф., Чередниченко Р.А. Распространение упругих волн в полосе, составленной из двух разнородных материалов // «Избранные проблемы прикладной механики» посвященный 60-тилетию академика В.Н. Челомея. М. : ВИНИТИ, 1974. С. 617—624.
6. Clifnon R.J. A difference method for plane problems in dynamic elasticity. Quart. Appl. Mfth. 1967. Vol 25. № 1, pp. 97—116.
Поступила в редакцию в июне 2012 г.
Об авторе: Немчинов Владимир Валентинович — кандидат технических наук, профессор кафедры прикладной механики и математики, Мытищинский филиал ФГБОУ ВПО «Московский государственный строительный университет», Московская обл., г. Мытищи, Олимпийский проспект, д. 50, 8(495) 583-73-81, vvnemchinov@gmail.com.
Для цитирования: Немчинов В.В. Двухслойная разностная схема численного решения плоских динамических задач теории упругости // Вестник МГСУ. 2012. № 8. С. 104—111.
V.V. Nemchinov
BILAYER DIFFERENCE SCHEME OF A NUMERICAL SOLUTION TO TWO-DIMENSIONAL DYNAMIC PROBLEMS OF ELASTICITY
Numerical modeling of dynamic problems of the theory of elasticity remains a relevant task. A complex network of waves that propagate within solid bodies, including longitudinal, transverse, conical and surface Rayleigh waves, etc., prevents the separation of wave fronts for modeling purposes. Therefore, it is required to apply the so-called "pass-through analysis".
The method applied to resolve dynamic problems of the two-dimensional theory of elasticity employs finite elements to approximate computational domains of complex shapes, whereby the software calculates the speed and voltage in the medium at each step. Pre-set boundary conditions are satisfied precisely.
The resulting method is classified as explicit bilayer difference schemes that form special relationships at the boundary points.
The method is based on an implicit bilayer time-difference scheme based on a system of dynamic equations of the theory of elasticity of the first order, which is converted into an explicit scheme with the help of a Taylor series in time, while basic relations are resolved with the help of the Galerkin method. The author demonstrates that the speed and voltage are calculated with the same accuracy as the one provided by the classical finite element method, whereby determination of stresses has to act as a numerically differentiating displacement.
The author identifies the relations needed to calculate both the internal points of the computational domain and the boundary points. The author has also analyzed the accuracy and convergence of the resulting method having completed a numerical simulation of the well-known problem of diffraction of a longitudinal wave speed in a circular aperture. The problem has an analytical solution.
Key words: finite-difference derivative, two-layer difference scheme, diffraction of a longitudinal wave, circle, contour stress, numerical modeling, finite element method.
References
1. Baron M.L., Matthews. Difraktsiya volny davleniya otnositel'no tsilindricheskoy polosti v upru-goy srede [Diffraction of a Pressure Wave with Respect to a Cylindrical Cavity in an Elastic Medium]. Prikladnaya mekhanika [Applied Mechanics]. A series, no. 3, 1961, pp. 31 —38.
2. Klifton R.Dzh. Raznostnyy metod v ploskikh zadachakh dinamicheskoy uprugosti [Difference Method for Plane Problems of Dynamic Elasticity]. Mekhanika [Mechanics]. 1968, no. 1 (107), pp. 103—122.
3. Musaev V.K. Primenenie metoda konechnykh elementov k resheniyu ploskoy nestatsionarnoy dinamicheskoy zadachi teorii uprugosti [Application of the Finite Element Method to Solve a Transient Dynamic Plane Elasticity Problem]. Mekhanika tverdogo tela [Mechanics of Solids]. 1980, no. 1, p. 167.
4. Musaev V.K. Vozdeystvie prodol'noy stupenchatoy volny na podkreplennoe krugloe otverstie v uprugoy srede [Impact of the Longitudinal Steo-shaped Wave on a Supported Circular Hole in an Elastic Medium]. All-Union Conference "Modern Problems of Structural Mechanics and Strength of Aircrafts." Collected abstracts. Moscow Institute of Aviation, 1983, p. 51.
5. Sabodash P.F, Cherednichenko R.A. Rasprostranenie uprugikh voln v polose, sostavlennoy iz dvukh raznorodnykh materialov [Propagation of Elastic Waves in a Band Composed of Two Dissimilar
Materials]. Collected works on "Selected Problems of Applied Mechanics" dedicated to the 60th Anniversary of Academician V.N. Chelomey. Moscow, VINITI, pp. 617—624.
6. Clifnon R.J. A Difference Method for Plane Problems in Dynamic Elasticity. Quart. Appl. Mfth. 1967, vol. 25, no. 1, pp. 97—116.
About the author: Nemchinov Vladimir Valentinovich — Candidate of Technical Sciences, Professor, Department of Applied Mechanics and Mathematics, Mytischi Branch, Moscow State University of Civil Engineering (MGSU), 50 Olimpiyskiy prospekt, Mytischi, Moscow Region, Russian Federation; vvnemchinov@gmail.com; +7 (495) 583-73-81.
For citation: Nemchinov V.B. Dvukhsloynaya raznostnaya skhema chislennogo resheniya ploskikh dinamicheskikh zadach teorii uprugosti [Bilayer Difference Scheme of a Numerical Solution to Two-Dimensional Dynamic Problems of Elasticity]. Vestnik MGSU [Proceedings of Moscow State University of Civil Engineering]. 2012, no. 8, pp. 104—111.
