Численный метод решения динамических задач теории упругости в полярной системе координат типа метода конечных элементов Текст научной статьи по специальности «Физика»

УДК 539.3

В.В. Немчинов, В.К. Мусаев

ФГБОУ ВПО «МГСУ»

ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ДИНАМИЧЕСКИХ ЗАДАЧ ТЕОРИИ УПРУГОСТИ В ПОЛЯРНОЙ СИСТЕМЕ КООРДИНАТ ТИПА МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ

Рассмотрен метод построения схемы решения динамических задач теории упругости в полярной системе координат. Получен численный метод, который определяет с одинаковой точностью значения скоростей и напряжений и точно выполняет заданные граничные условия. В качестве тестового примера рассмотрена задача о дифракции продольной волны на круглом отверстии. Произведена оценка точности и сходимости численного решения в зависимости от величины дискретного шага по координатам и времени.

Ключевые слова: конечно-разностная производная, двухслойная по времени схема, дифракция продольной волны, круговая полость, контурные напряжения, численное моделирование, метод конечных элементов, полярная система координат.

Вывод разностной схемы в полярной системе координат имеет преимущественно методическое значение, иллюстрирующее применение общей методики построения расчетных схем, основанных на динамических уравнениях теории упругости первого порядка в скоростях и напряжениях [1]. В результате для любых систем координат получаются расчетные соотношения, позволяющие находить скорости (перемещения) и напряжения, во всех точках произвольной расчетной области, используя конечные элементы (КЭ) (треугольные с тремя узлами или четырехугольные с четырьмя узлами). Аналогично, для объемной упругой задачи: конечные элементы — 6-точечные тетраэдры или 8-точечная призма.

Численное решение динамических задач теории упругости имеет наибольшую точность, если границы рассматриваемых областей совпадают с координатными линиями выбранной системы координат.

Для динамических уравнений теории упругости первого порядка в полярной систему Рис. 1. Полярная система координат (рис. 1) можно ввести безразмерные к°°рдинат (г,ф) и пл°щадь ш-координаты лилинейного КЭ = гё фёг

Координаты: г = г'Ь; ф = ф

время и скорости: t' = с^Ь; и' = и!ср ; V' = у/Ср ; (1)

напряжения: р = Огг/рС2Р; q = °фф/Рс] ; * = °гф/р

С2:

где Ь — характерный линейный размер; р — плотность; ср — скорость продольной волны.

Безразмерная система уравнений динамической теории упругости первого прядка в скоростях и напряжениях в полярной системе координат имеет следующий вид (штрихи опускаем)

д u dp ds p — q — = — +---—-;

dt dr r дф r dv _ds + dq + 2s dt dr r дф r

dp du

— =— + a dt dr

dq du

— = a— + dt dr

^ dv u л +

^гЭф rj

r dv u л +

ds (dv

~dt

^г5ф ry

du Л

---+ -

dr r rдф

(2.1) (2.2)

(2.3)

(2.4)

(2.5)

Рассмотрим неявную двухслойную по времени разностную схему, аппроксимирующую дифференциальную систему уравнений (2.1)—(2.5), где для первого уравнения

0 П 0

р — а

(3)

u(t + At)-u (t) dp ' ds' p - q

— _ ut _— +---+-

At dr гдф r

Остальные уравнения получаются аналогично.

u (t + At) + u (t) At du

u0 =

A t

' u + ■

At [dp ds p - q

= u +—I —+-+-—-

2 y dr гдф r

2 д t

\

(4)

в = (е,/ер )2, а = 1 — 2р.

Для вывода разностных уравнений применим метод Галёркина [2]. Умножим соотношение (2.1)—(2.5) на базисную функцию КЭ N (г,ц>) и проинтегрируем по площади КЭ, тогда для первого уравнения

\uiNfd фёг = Г— ЫГ фёг + Г-^— Ы^ё фёг з з дг гдф

+

dr

•Р - q

+

(5)

Nrd q>dr.

Меняем порядок дифференцирования (формула Грина):

J ut N rd tydr + [ p0 d—Lrd tydr + [ s0 d—Lrd tydr = s S dr S гдФ

= ф (p0 rd ф - s0 dr)— + J N rd qdr.

I S r

Подставляя в (6) соотношения (4), получаем

(6)

ВЕСТНИК

МГСУ-

7/2013

с г дЫ ■ г дЫ ■

I иг Ы/ёфёг + I р—г-гёфёг + I ^—г-гёфёг J 1 дг ^ г дф

Л/ 2

ди дг

а

дг ду и

\\

Угдф г))

дЫ

дг

гё фёг

ду у ди

дг г г дф ) г дф

дЫ

гё фёг

| (р0 гё ф - ^0 ёг)) РN \гё фёг.

(7)

Остальные уравнения преобразуются аналогично.

Представляем значения скоростей и напряжений через базисные функции на элементе

и = N. и, V = N. р = N. р , q = N. qi, * = N. * ■. Обозначим

• т

дг

(8)

■ г дЫ-

—'-Ы,гё<аёг, Бп = I—N,гё(пёг:

3 дг 1 * 11 {гду 1

сдN. дЫ, сдN. дЫ,

з дг дг 3 г дф дг

■> дг гдт ^ гдф Гдф

5 дг Г дф

(9)

-г-^-гйфёг, Н; = {—-гёфёг;

5 дг г 2 г дф г

I. = фёг, М . = Г N N. гёфёг.

J г *

Получаем соотношения для расчета как внутренних точек, так и граничных точек области (интегралы вычисляются квадратурами Гаусса [3]).

М .и! + Ар1 + В +—\С .и] + аЕ..у] +

1 * Ч о _ ч ч

+

р1 + Г.. у' + Н иии ] = ф (р\ёф - ^°ёт) ) +

(10)

+ I. ( - Ч )•

Преобразуем (10) и остальные соотношения в более простой вид ми = Ги, + ф( р0гёф-50ёг));

I

М ^ = ГУ1 + ф ( ^0 гё ф- д0 ёг));

1

М ур] = ¥р х + ^ ( м0 гё ф- а V0 ёг));

(111) (11.2) (11.3)

= ^ + ф (а и0тйф - V0йт )); (11.4)

1

М^ = + в ф ( V0гаф - и0&)) . (11.5)

1

Соотношения (11.1)—(115) будут иметь явный вид, если матрицу масс

£

М г представить в диагональном виде. Для четырех точечного КЭ М и = т =—.

V

Правые части соотношений (11.1)—(115) (для удобства пишем со знаком минус) имеет следующий вид:

- Еи , = А.р1 + В+ — \Спи] + аЕ„У +

1 У* У 2

+ Р' + Е.у' + Иви'] + !„(р - я>); (12.1)

- »» + ^ +| [С»-" - «V +

+ Ем1 + аБм1 + Е7 V1 + Н » и1 1 + 21.У; (12.2)

V V V V _| V

- Ер = Ли1 + аВ у1 +—\C-.p1 + Е^1 +

г 1 V V 2

+а Ц1$] +а Ен а11 + аI ни]; (12.3)

1 1 1

- Еа . = а А.и + В.. + — Га С.. р + аЕ+

лг у у 2

+ Б,/ + Ег]а] ]+ I .и1; (12.4)

- й, = РV +РЩи' + Pffas' + ' +

+ D'P + F'SJ J-P/'V. (12.5) Используем соотношения типа (4), получаем

= Fui + ф(р • rdy-sdr )N,; (13.1)

Щ = Fvt + §(srdq-qdr)); (13.2)

^Pi = FPi + §{u' rdф-av■ dr)N1; (13.3)

= Fqt + £(au ■ rdф-v■ dr)); (13.4)

^s, = Fs. +efy(yrdy-udr )Nt. (13.5) Тогда

m;u't = ^ut §(ptrd ф-stdr )Nt; (14.1)

i

тХ = Щ +—)); (14.2)

тр = Жр1 +—§(и{гё ф-ау{ёг ; (14.3)

2 ,

ш.

= + ~ у-^ж)); (14.4)

т£ = ^ + р — §(гёц-ы(ёг)). (14.5)

2 I

Рассмотрим четырехугольный КЭ с 4 узлами и базисными функциями [2—4]

^(г;,л)=0,25(1-г;)(1-л); N2 (£,п) =0,25(1 + ^)(1 -Л); (15)

ЛГ3(4,Л) =0,25(1 + 4)(1 + л); N4(£,п) = 0,25(1 Ч)(1 + п) .

Если граница проведена так, чтобы г ёф >0; ёг = 0, то интегрируем контурные интегралы по границе элемента, суммируя в граничной точке.

(16.1)

т = (16.2)

щр = Щ ф; (16.3)

а и(гё Ф; (16.4)

= ^ + ф. (16.5)

Из соотношений (16.1)—(16.5) можно получить 3 также линейно независимые комбинации, чтобы знаменатель выражения не обращался в ноль ни при каком выборе dt, гёф, ёг — шагов по времени и координате

и- Р =—г м ; (171)

т. +--гё ф

г 2

Щ Я*

-тв *=—Л—; (172)

т^ + УР — гёф

а р^- £ =-£-Ъ.. (17.3)

щ

Полученные соотношения подобны соотношениям метода [5], полученным для квадратной дискретной сетки. Для свободной границы имеем

Р1 = 0, si = 0,

и из формул (17.1)—(17.3) можно получить остальные значения переменных на границе ui, V , qi.

В качестве примера рассмотрим известную задачу о дифракции продольной волны на круглом отверстии [4, 6—9] (туннельная выработка), имеющую как аналитическое решение, так и много экспериментальных и расчетных данных. При коэффициенте Пуассона V = 1/ 3 (отношение скоростей поперечной и продольной волн с^ер = 0,5) максимальное значение напряжения на границе круглого отверстия равно -2,97 при ? = 3,5 .

Разбивка расчетной области осуществлялась в полярной системе координат (верхняя часть полуплоскости), шаг по углу dфф равнялся шагу по радиусу dr = М ф, за исключением 6-го варианта, где dr = ^ ф/2.

Результаты расчетов приведены в таблице и на рис. 2, где показаны значение напряжения афф в верхней точки окружности, аналитическое значение афф = — 2,97 . Максимум напряжений достигается примерно при безразмерном времени ? « 3,5.

Максимальные значения напряжения афф в верхней точке окружности в зависимости от числа разбиения дуги полуокружности

Номер графика Шаг по координатам Число расчетных точек на полуокружности Max а фф Относительная точность, %

1. Ar = ЯДф 0,049 17 -2,667 10,2

2. Ar = ЯДф 0,049 17 -2,641 11

3. Ar = ЯДф 0,025 33 -2,735 7,9

4. Ar = КДф 0,016 49 -2,84 4,37

5. Ar = ЯДф 0,012 65 -2,98 0,34

6. Ar = ЛДф/ 2 0,016 49 -2,92 1,7

Рис. 2. Изменение напряжения а во времени в верхней точке окружности

Полученные численные результаты показывают хорошую сходимость данного численного метода, но численные эксперименты показали, что ошибка решения на границе окружности линейно зависит от выбранного шага по координатам.

Библиографический список

1. НемчиновВ.В. Двухслойная разностная схема численного решения плоских динамических задач теории упругости // Вестник МГСУ 20i2. № S. С. i04—iii.

2. Флетчер К. Численные методы на основе метода Галёркина. М. : Мир, i9SS. 352 с.

3. СекуловичМ. Метод конечных элементов. М. : Стройиздат, i993. 664 с.

4. Мусаев В.К. Применение метода конечных элементов к решению плоской нестационарной динамической задачи теории упругости // Механика твердого тела. i9S0. № i. С. i67—i73.

5. Сабодаш П.Ф., Чередниченко Р.А. Применение метода пространственных характеристик к решению задач о распространении волн в упругой полуполосе // Известия АН СССР. Механ. твердого тела. i972. № 6. С. iS0—iS5.

6. Гернет Х., Крузе-Паскаль Д. Неустановившаяся реакция находящегося в упругой среде кругового цилиндра произвольной толщины на действие плоской волны расширения // Прикладная механика. Труды американского общества инженеров-механиков. Сер. Е. i966. Т. 33. № 3. С. 4S—60.

7. Bayandin Yu.V., Naimark O.B., Uvarov S.V. Numerical simulation of spall failure in metals under shock compression // AIP Conf. Proc. of the American Physical Society Topical Group on Shock Compression of Condensed Matter, Nashville, TN, 2S June — 3 July 2009. V. П95. Pр. i093—1096. Режим доступа: http://dx.doi.org/i0.i063/L3294992

S. Burago N.G., Zhuravlev A.B., Nikitin I.S. Models of Multiaxial Fatigue Fracture and Service Life Estimation of Structural Elements // Mechanics of Solids. 20ii. Vol. 46. No. 6. Pp. S2S—S3S.

9. Li Y., Liu G.R., ZhangG.Y. An adaptive NS/ES-FEM approach for 2D contact problems using triangular elements // Finite Elem. Anal. Des. 20ii. V. 47, N. 3. Pр. 256—275. Режим доступа: http://dx.doi.org/i0.i0i6/j.finel.20i0.i0.007

Поступила в редакцию в феврале 2013 г.

Об авторах: Немчинов Владимир Валентинович — кандидат технических наук, профессор кафедры прикладной механики и математики, Мытищинский филиал, ФГБОУ ВПО «Московский государственный строительный университет» (ФГБОУ ВПО «МГСУ»), i4i006, Московская область, г. Мытищи, Олимпийский проспект, д. 50, S(495)602-70-29, vvnemchinov@mail.ru;

Мусаев Вячеслав Кадыр оглы—доктор технических наук, профессор, профессор-консультант, Мытищинский филиал, ФГБОУ ВПО «Московский государственный строительный университет» (ФГБОУ ВПО «МГСУ»), i4i006, Московская область, г. Мытищи, Олимпийский проспект, д. 50, musayev-vk@yandex.ru.

Для цитирования: Немчинов В.В., Мусаев В.К. Численный метод решения динамических задач теории упругости в полярной системе координат типа метода конечных элементов // Вестник МГСУ 20i3. № 7. С. 6S—76

V.V. Nemchinov, V.K. Musayev

NUMERICAL METHOD FOR SOLVING DYNAMIC PROBLEMS OF THE THEORY OF ELASTICITY IN THE POLAR COORDINATE SYSTEM SIMILAR TO THE FINITE

ELEMENT METHOD

The authors consider a dynamic problem solving procedure based on the theory of elasticity in the Cartesian coordinate system. This method consists in the development

of the pattern of numerical solutions to dynamic elastic problems within any coordinate system and, in particular, in the polar coordinate system. Numerical solutions of dynamic problems within the theory of elasticity are the most accurate ones, if the boundaries of the areas under consideration coincide with the coordinate lines of the selected coordinate system.

The first order linear system of differential equations is converted into an implicit difference scheme. The implicit scheme is transformed into the explicit method of numerical solutions. Using the Galerkin method, the authors obtain formulas for the calculation of both the points of the computational domain and the boundary points.

Difference ratios similar to those obtained for a discrete rectangular grid and derived in this paper are suitable to design any geometry, which fact significantly increases the value of the methods considered in this paper.

As a test case, the problem of diffraction of a longitudinal wave in a circular cavity, where maximum stresses are obtained analytically, was considered by the authors. The proposed method demonstrated sufficient accuracy of calculations and convergence of numerical solutions, depending on the size of discrete steps. The problem of diffraction of longitudinal waves in a circular cavity was taken for example; however, the proposed method is applicable to any problems within any computational domain.

The polar coordinate system is the best one for any research into the diffraction of plane longitudinal waves in a circular cavity, since the boundaries of the computational domain coincide with the coordinate lines of the selected system.

Key words: finite-difference derivative, double-layer scheme, diffraction of longitudinal waves, circular cavity contour stress, numerical modeling, finite element method, polar coordinate system.

References

1. 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.

2. Fletcher K. Chislennye metody na osnove metoda Galerkina [Numerical Methods Based on the Galerkin Method]. Moscow, Mir Publ., 1988, 352 p.

3. Sekulovich M. Metod konechnykh elementov [Finite Element Method]. Moscow, Stroy-izdat Publ., 1993, 664 p.

4. Musaev V.K. Primenenie metoda konechnykh elementov k resheniyu ploskoy nestat-sionarnoy dinamicheskoy zadachi teorii uprugosti [Application of the Finite Element Method to the Plane Non-stationary Dynamic Problem of the Theory of Elasticity]. Mekhanika tverdogo tela [Solid Body Mechanics]. 1980, no. 1, pp. 167—173.

5. Sabodash P.F., Cherednichenko R.A. Primenenie metoda prostranstvennykh kharak-teristik k resheniyu zadach o rasprostranenii voln v uprugoy polupolose [Application of Method of 3D Characteristics to Problems of Propagation of Waves in an Elastic Half-strip]. Izvestiya AN SSSP. Mekhan. tverdogo tela [News of the Academy of Sciences of the USSR. Solid Body Mechanics]. 1972, no. 6, pp. 180—185.

6. Gernet Kh., Kruze-Paskal' D. Neustanovivshayasya reaktsiya nakhodyashchegosya v uprugoy srede krugovogo tsilindra proizvol'noy tolshchiny na deystvie ploskoy volny rasshire-niya [Unstable Response of an Arbitrary Thickness Circular Cylinder to the Action of a Plane Expansion Wave]. Prikladnaya mekhanika. Trudy amerikanskogo obshchestva inzhenerov-mekhanikov. Ser. E. [Applied Mechanics. Works of the American Society of Mechanical Engineers. Series E.] 1966, vol. 33, no. 3, pp. 48—60.

7. Bayandin Yu.V., Naimark O.B., Uvarov S.V. Numerical Simulation of Spall Failure in Metals under Shock Compression. AIP Conf. Proc. of the American Physical Society. Topical Group on Shock Compression of Condensed Matter. Nashville, TN, 28 June — 3 July 2009, vol. 1195, pp. 1093—1096.

8. Burago N.G., Zhuravlev A.B., Nikitin I.S. Models of Multiaxial Fatigue Fracture and Service Life Estimation of Structural Elements. Mechanics of Solids. 2011, vol. 46, no. 6, pp. 828—838.

9. Li Y., Liu G.R., Zhang G.Y. An Adaptive NS/ES-FEM Approach for Plane Contact Problems Using Triangular Elements. Finite Elem. Anal. Dec. 2011, vol., 47, no. 3, pp. 256—275.

About the authors: 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, 141006, Russian Federation; vvnemchinov@mail.ru; +7 (495) 602-70-29;

Musayev Vyacheslav Kadyr ogly — Doctor of Technical Sciences, Professor, Consulting Professor, Mytischi Branch, Moscow State University of Civil Engineering (MGSU), 50 Olimpiyskiy prospekt, Mytischi, Moscow Region, 141006, Russian Federation; musayev-vk@ yandex.ru.

For citation: Nemchinov V.V., Musaev V.K. Chislennyy metod resheniya dinamicheskikh zadach teorii uprugosti v polyarnoy sisteme koordinat tipa metoda konechnykh elementov [Numerical Method for Solving Dynamic Problems of the Theory of Elasticity in the Polar Coordinate System Similar to The Finite Element Method]. Vestnik MGSU [Proceedings of Moscow State University of Civil Engineering]. 2013, no. 7, pp. 68—76.